Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tracklet.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Tracklet.h
1 #ifndef TRACKLET_H
2 #define TRACKLET_H
3 
4 #include <TTree.h>
5 
6 #include "GenHadron.h"
7 #include "Hit.h"
8 #include "Utilities.h"
9 
10 using namespace std;
11 
12 class Tracklet : public TObject
13 {
14  public:
15  Tracklet();
16  Tracklet(Hit *hit1, Hit *hit2);
17  ~Tracklet();
18 
19  void SetHits(Hit *hit1, Hit *hit2);
20  Hit *Hit1();
21  Hit *Hit2();
22  float dEta();
23  float dPhi();
24  float dR();
25  float Eta();
26  float Phi();
27  float tklVtxZ();
28  int LayerComb(); // Layer combination: can only be 12, 23, 13
29  void SetMatchedGenHardon();
30  bool IsMatchedGenHadron();
31  void SetGenHadron(GenHadron *genhad);
32  GenHadron *MatchedGenHadron();
33 
34  private:
35  Hit *_hit1;
36  Hit *_hit2;
37  float _deta;
38  float _dphi;
39  bool _matched_genhadron;
40  GenHadron *_genhad;
41 };
42 
43 Tracklet::Tracklet()
44 {
45  _hit1 = nullptr;
46  _hit2 = nullptr;
47  _matched_genhadron = false;
48  _genhad = nullptr;
49 }
50 
51 Tracklet::Tracklet(Hit *hit1, Hit *hit2)
52 {
53  _hit1 = hit1;
54  _hit2 = hit2;
55  _deta = _hit2->Eta() - _hit1->Eta();
56  _dphi = deltaPhi(_hit2->Phi(), _hit1->Phi());
57  _matched_genhadron = false;
58  _genhad = nullptr;
59 }
60 
61 Tracklet::~Tracklet() {}
62 
63 void Tracklet::SetHits(Hit *hit1, Hit *hit2)
64 {
65  _hit1 = hit1;
66  _hit2 = hit2;
67  _deta = _hit2->Eta() - _hit1->Eta();
68  _dphi = deltaPhi(_hit2->Phi(), _hit1->Phi());
69  _matched_genhadron = false;
70  _genhad = nullptr;
71 }
72 
73 Hit *Tracklet::Hit1() { return (_hit1); }
74 
75 Hit *Tracklet::Hit2() { return (_hit2); }
76 
77 float Tracklet::dEta() { return (_deta); }
78 
79 float Tracklet::dPhi() { return (_dphi); }
80 
81 float Tracklet::dR() { return (sqrt(_deta * _deta + _dphi * _dphi)); }
82 
83 float Tracklet::Eta() { return (0.5 * (_hit1->Eta() + _hit2->Eta())); }
84 
85 float Tracklet::Phi() { return (0.5 * (_hit1->Phi() + _hit2->Phi())); }
86 
87 float Tracklet::tklVtxZ() { return (_hit1->posZ() - _hit1->rho() * (_hit2->posZ() - _hit1->posZ()) / (_hit2->rho() - _hit1->rho())); }
88 
89 int Tracklet::LayerComb() { return (_hit1->Layer()) * 10 + (_hit2->Layer()); }
90 
91 void Tracklet::SetMatchedGenHardon() { _matched_genhadron = true; }
92 
93 bool Tracklet::IsMatchedGenHadron() { return _matched_genhadron; }
94 
95 void Tracklet::SetGenHadron(GenHadron *genhad)
96 {
97  _matched_genhadron = true;
98  _genhad = genhad;
99 }
100 
101 GenHadron *Tracklet::MatchedGenHadron() { return _genhad; }
102 
103 /*------------------------------------------------------------------------------------------*/
104 
105 class TrackletData
106 {
107  public:
108  bool isdata;
109  vector<vector<Hit *>> layers = {{}, {}};
110  vector<Tracklet *> ProtoTkls, RecoTkls, RecoTkls_GenMatched;
111  vector<GenHadron *> GenHadrons;
112 
113  int event, NClusLayer1, NPrototkl, NRecotkl_Raw, NRecotkl_GenMatched, NGenHadron;
114  float PV_x, PV_y, PV_z, TruthPV_x, TruthPV_y, TruthPV_z;
115  float Centrality_mbd, Centrality_mbdquantity;
116  bool pu0_sel, trig;
117  int process; // single diffractive process
118  vector<float> clusphi, cluseta, clusphisize, cluszsize;
119  vector<float> prototkl_eta, prototkl_phi, prototkl_deta, prototkl_dphi, prototkl_dR;
120  vector<float> recotklraw_eta, recotklraw_phi, recotklraw_deta, recotklraw_dphi, recotklraw_dR;
121  vector<float> recotklgm_eta, recotklgm_phi, recotklgm_deta, recotklgm_dphi, recotklgm_dR;
122  vector<float> GenHadron_Pt, GenHadron_eta, GenHadron_phi, GenHadron_E, GenHadron_matched_Pt, GenHadron_matched_eta, GenHadron_matched_phi, GenHadron_matched_E;
123 };
124 
125 void SetMinitree(TTree *outTree, TrackletData &tkldata)
126 {
127  outTree->Branch("event", &tkldata.event);
128  outTree->Branch("NClusLayer1", &tkldata.NClusLayer1);
129  outTree->Branch("NPrototkl", &tkldata.NPrototkl);
130  outTree->Branch("NRecotkl_Raw", &tkldata.NRecotkl_Raw);
131  outTree->Branch("Centrality_mbd", &tkldata.Centrality_mbd);
132  outTree->Branch("Centrality_mbdquantity", &tkldata.Centrality_mbdquantity);
133  outTree->Branch("PV_x", &tkldata.PV_x);
134  outTree->Branch("PV_y", &tkldata.PV_y);
135  outTree->Branch("PV_z", &tkldata.PV_z);
136  outTree->Branch("trig", &tkldata.trig);
137  outTree->Branch("clusPhi", &tkldata.clusphi);
138  outTree->Branch("clusEta", &tkldata.cluseta);
139  outTree->Branch("clusPhiSize", &tkldata.clusphisize);
140  outTree->Branch("clusZSize", &tkldata.cluszsize);
141  outTree->Branch("prototkl_eta", &tkldata.prototkl_eta);
142  outTree->Branch("prototkl_phi", &tkldata.prototkl_phi);
143  outTree->Branch("prototkl_deta", &tkldata.prototkl_deta);
144  outTree->Branch("prototkl_dphi", &tkldata.prototkl_dphi);
145  outTree->Branch("prototkl_dR", &tkldata.prototkl_dR);
146  outTree->Branch("recotklraw_eta", &tkldata.recotklraw_eta);
147  outTree->Branch("recotklraw_phi", &tkldata.recotklraw_phi);
148  outTree->Branch("recotklraw_deta", &tkldata.recotklraw_deta);
149  outTree->Branch("recotklraw_dphi", &tkldata.recotklraw_dphi);
150  outTree->Branch("recotklraw_dR", &tkldata.recotklraw_dR);
151 
152  if (!tkldata.isdata)
153  {
154  outTree->Branch("NRecotkl_GenMatched", &tkldata.NRecotkl_GenMatched);
155  outTree->Branch("TruthPV_x", &tkldata.TruthPV_x);
156  outTree->Branch("TruthPV_y", &tkldata.TruthPV_y);
157  outTree->Branch("TruthPV_z", &tkldata.TruthPV_z);
158  outTree->Branch("pu0_sel", &tkldata.pu0_sel);
159  outTree->Branch("process", &tkldata.process);
160  outTree->Branch("recotklgm_eta", &tkldata.recotklgm_eta);
161  outTree->Branch("recotklgm_phi", &tkldata.recotklgm_phi);
162  outTree->Branch("recotklgm_deta", &tkldata.recotklgm_deta);
163  outTree->Branch("recotklgm_dphi", &tkldata.recotklgm_dphi);
164  outTree->Branch("recotklgm_dR", &tkldata.recotklgm_dR);
165  outTree->Branch("GenHadron_Pt", &tkldata.GenHadron_Pt);
166  outTree->Branch("GenHadron_eta", &tkldata.GenHadron_eta);
167  outTree->Branch("GenHadron_phi", &tkldata.GenHadron_phi);
168  outTree->Branch("GenHadron_E", &tkldata.GenHadron_E);
169  outTree->Branch("GenHadron_matched_Pt", &tkldata.GenHadron_matched_Pt);
170  outTree->Branch("GenHadron_matched_eta", &tkldata.GenHadron_matched_eta);
171  outTree->Branch("GenHadron_matched_phi", &tkldata.GenHadron_matched_phi);
172  outTree->Branch("GenHadron_matched_E", &tkldata.GenHadron_matched_E);
173  }
174 }
175 
176 void ResetVec(TrackletData &tkldata)
177 {
178  for (size_t i = 0; i < tkldata.layers.size(); i++)
179  {
180  CleanVec(tkldata.layers[i]);
181  }
182  CleanVec(tkldata.ProtoTkls);
183  CleanVec(tkldata.RecoTkls);
184  CleanVec(tkldata.prototkl_eta);
185  CleanVec(tkldata.prototkl_phi);
186  CleanVec(tkldata.prototkl_deta);
187  CleanVec(tkldata.prototkl_dphi);
188  CleanVec(tkldata.prototkl_dR);
189  CleanVec(tkldata.recotklraw_eta);
190  CleanVec(tkldata.recotklraw_phi);
191  CleanVec(tkldata.recotklraw_deta);
192  CleanVec(tkldata.recotklraw_dphi);
193  CleanVec(tkldata.recotklraw_dR);
194  CleanVec(tkldata.recotklgm_eta);
195  CleanVec(tkldata.recotklgm_phi);
196  CleanVec(tkldata.recotklgm_deta);
197  CleanVec(tkldata.recotklgm_dphi);
198  CleanVec(tkldata.recotklgm_dR);
199  CleanVec(tkldata.GenHadrons);
200  CleanVec(tkldata.GenHadron_Pt);
201  CleanVec(tkldata.GenHadron_eta);
202  CleanVec(tkldata.GenHadron_phi);
203  CleanVec(tkldata.GenHadron_E);
204  CleanVec(tkldata.GenHadron_matched_Pt);
205  CleanVec(tkldata.GenHadron_matched_eta);
206  CleanVec(tkldata.GenHadron_matched_phi);
207  CleanVec(tkldata.GenHadron_matched_E);
208 }
209 
210 bool compare_dR(Tracklet *a, Tracklet *b) { return a->dR() < b->dR(); }
211 
212 void ProtoTracklets(TrackletData &tkldata, float cutdr)
213 {
214  float Cut_dR = cutdr; // Nominal: 0.5
215 
216  int iComb = 0;
217  // Get all combinations
218  for (auto ihitl1 : tkldata.layers[0])
219  {
220  for (auto ihitl2 : tkldata.layers[1])
221  {
222  // iComb++;
223  // if (iComb % 500000 == 0)
224  // fprintf(stderr, "processing %d of %d combinations (%.3f %%)\n", iComb, nhits1 * nhits2, (float)iComb / (nhits1 * nhits2) * 100.);
225 
226  float dEta = ihitl2->Eta() - ihitl1->Eta();
227  float dPhi = deltaPhi(ihitl2->Phi(), ihitl1->Phi());
228  float dR = sqrt(dEta * dEta + dPhi * dPhi);
229 
230  if (dR < Cut_dR)
231  {
232  Tracklet *tmptkl = new Tracklet(ihitl1, ihitl2);
233  tkldata.ProtoTkls.push_back(tmptkl);
234 
235  tkldata.prototkl_eta.push_back(tmptkl->Hit1()->Eta());
236  tkldata.prototkl_phi.push_back(tmptkl->Hit1()->Phi());
237  tkldata.prototkl_deta.push_back(tmptkl->dEta());
238  tkldata.prototkl_dphi.push_back(tmptkl->dPhi());
239  tkldata.prototkl_dR.push_back(tmptkl->dR());
240  }
241  else
242  continue;
243  }
244  }
245 
246  tkldata.NPrototkl = tkldata.ProtoTkls.size();
247 }
248 
249 void RecoTracklets(TrackletData &tkldata)
250 {
251  sort(tkldata.ProtoTkls.begin(), tkldata.ProtoTkls.end(), compare_dR);
252 
253  int itkl = 0;
254  for (auto &tkl : tkldata.ProtoTkls)
255  {
256  if (tkl->Hit1()->IsMatchedTkl() || tkl->Hit2()->IsMatchedTkl())
257  {
258  continue;
259  }
260  else
261  {
262  tkldata.RecoTkls.push_back(tkl);
263  tkl->Hit1()->SetMatchedTkl();
264  tkl->Hit2()->SetMatchedTkl();
265 
266  tkldata.recotklraw_eta.push_back(tkl->Hit1()->Eta());
267  tkldata.recotklraw_phi.push_back(tkl->Hit1()->Phi());
268  tkldata.recotklraw_deta.push_back(tkl->dEta());
269  tkldata.recotklraw_dphi.push_back(tkl->dPhi());
270  tkldata.recotklraw_dR.push_back(tkl->dR());
271  }
272  }
273 
274  tkldata.NRecotkl_Raw = tkldata.RecoTkls.size();
275 }
276 
277 void GenMatch_Recotkl(TrackletData &tkldata)
278 {
279  for (auto &tkl : tkldata.RecoTkls)
280  {
281  if (tkl->IsMatchedGenHadron())
282  continue;
283  for (auto &ghadron : tkldata.GenHadrons)
284  {
285  if (ghadron->IsMatchedToRecotkl() || tkl->IsMatchedGenHadron())
286  continue;
287  // Matching criteria
288  if (deltaR(tkl->Eta(), tkl->Phi(), ghadron->Eta(), ghadron->Phi()) > 0.1)
289  continue;
290  else
291  {
292  tkl->SetMatchedGenHardon();
293  ghadron->SetMatchedToRecotkl();
294  tkl->SetGenHadron(ghadron);
295  tkldata.RecoTkls_GenMatched.push_back(tkl);
296  tkldata.recotklgm_eta.push_back(tkl->Hit1()->Eta());
297  tkldata.recotklgm_phi.push_back(tkl->Hit1()->Phi());
298  tkldata.recotklgm_deta.push_back(tkl->dEta());
299  tkldata.recotklgm_dphi.push_back(tkl->dPhi());
300  tkldata.recotklgm_dR.push_back(tkl->dR());
301  tkldata.GenHadron_matched_Pt.push_back(ghadron->Pt());
302  tkldata.GenHadron_matched_eta.push_back(ghadron->Eta());
303  tkldata.GenHadron_matched_phi.push_back(ghadron->Phi());
304  tkldata.GenHadron_matched_E.push_back(ghadron->E());
305  }
306  }
307  }
308 
309  tkldata.NRecotkl_GenMatched = tkldata.RecoTkls_GenMatched.size();
310 }
311 
312 #endif