Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrackingPerformanceCheck.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrackingPerformanceCheck.C
2 
3 #include <phool/getClass.h>
5 
7 #include <g4main/PHG4Particle.h>
8 #include <g4main/PHG4VtxPoint.h>
9 
10 #include <g4hough/SvtxVertexMap.h>
11 #include <g4hough/SvtxVertex.h>
12 #include <g4hough/SvtxTrackMap.h>
13 #include <g4hough/SvtxTrack.h>
14 #include <g4hough/SvtxClusterMap.h>
15 #include <g4hough/SvtxCluster.h>
16 
17 #include <g4eval/SvtxEvalStack.h>
18 #include <g4eval/SvtxTrackEval.h>
19 #include <g4eval/SvtxVertexEval.h>
20 #include <g4eval/SvtxTruthEval.h>
21 
22 #include "TMath.h"
23 #include <TH1F.h>
24 #include <TH2F.h>
25 #include <TH1D.h>
26 #include <TH2D.h>
27 #include "TString.h"
28 
29 #include <iostream>
30 
31 using namespace std;
32 
34  SubsysReco(name),
35  fLayerTPCBegins(7),
36  fReconstructable_TPCHits(30), // 40/60 or 30/40
37  fFair_ClustersContribution(20), // 25/67 or 20/47
38  fGTrack_Chi2NDF(2.0),
39  fGTrack_TPCClusters(20) // 25/60 or 20/40
40 {
41 }
42 
44  // register histograms
46 
47  const int nptbinsPONE = 52;
48  Int_t nptbins = nptbinsPONE - 1;
49  Double_t ptbins[nptbinsPONE];
50  ptbins[0] = 0.2;
51  for(int i=1; i!=9; ++i)
52  ptbins[i] = 0.2 + i*0.2;
53  for(int i=0; i!=11; ++i)
54  ptbins[9+i] = 2.0 + i*0.5;
55  for(int i=0; i!=10; ++i)
56  ptbins[20+i] = 8.0 + i*1.0;
57  for(int i=0; i!=12; ++i)
58  ptbins[30+i] = 18.0 + i*2.0;
59  for(int i=0; i!=10; ++i)
60  ptbins[42+i] = 45.0 + i*5.0;
61 
62  for(int i=0; i!=nptbinsPONE; ++i) std::cout << i << " " << ptbins[i] << std::endl;
63 
64  fHNEvents = new TH1F("Events","Events",1,-0.5,0.5);
66 
67  TString sTname[5] = {"AllG4Particles","Primaries","Embedded","Reconstructables","RecoMatched"};
68  for(int i=0; i!=5; ++i) {
69  fHTN[i] = new TH1F(Form("Truths%d_N",i),Form("Number of %s",sTname[i].Data()),50,-0.5,49.5);
70  fHTPt[i] = new TH1F(Form("Truths%d_Pt",i),Form("%s Pt",sTname[i].Data()),nptbins,ptbins);
71  fHTPhi[i] = new TH1F(Form("Truths%d_Phi",i),Form("%s Phi",sTname[i].Data()),100,0,TMath::TwoPi());
72  fHTEta[i] = new TH1F(Form("Truths%d_Eta",i),Form("%s Eta",sTname[i].Data()),100,-10,+10);
73  se->registerHisto(fHTN[i]);
74  se->registerHisto(fHTPt[i]);
75  se->registerHisto(fHTPhi[i]);
76  se->registerHisto(fHTEta[i]);
77  if(i<2) continue;
78  fHTNHits[i] = new TH1F(Form("Truths%d_Hits",i),Form("Number of TPCHits in %s",sTname[i].Data()),61,-0.5,60.5);
79  se->registerHisto(fHTNHits[i]);
80  if(i<4) continue;
81  fHTChi2[i] = new TH2F(Form("Truths%d_Chi2NDF",i),Form("%s Chi2NDF;PT;CHI2NDF",sTname[i].Data()),nptbins,ptbins,100,0,5);
82  fHTDca2D[i] = new TH2F( Form("Truths%d_Dca2D",i),Form("%s Dca2D;PT;DCA2D",sTname[i].Data()),nptbins,ptbins,200,-0.003,+0.003);
83  fHTNClustersContribution[i] = new TH2F(Form("Truths%d_NClusters",i), Form("%s Number of Clusters Matched",sTname[i].Data()), nptbins,ptbins,70,-0.5,69.5);
84  fHTPtResolution[i] = new TH2F(Form("Truths%d_ResPt",i), Form("%s PtResolution;PT;REL DIFF",sTname[i].Data()), nptbins,ptbins,400,-1.1,+1.1);
85  fHTPtResolution2[i] = new TH2F(Form("Truths%d_Res2Pt",i), Form("%s PtResolution2;PT;DIFF",sTname[i].Data()), nptbins,ptbins,400,-0.2,+0.2);
86  fHTPhiResolution[i] = new TH2F(Form("Truths%d_ResPhi",i), Form("%s PhiResolution;PT;REL DIFF",sTname[i].Data()), nptbins,ptbins,200,-0.1,+0.1);
87  fHTEtaResolution[i] = new TH2F(Form("Truths%d_ResEta",i), Form("%s EtaResolution;PT;REL DIFF",sTname[i].Data()), nptbins,ptbins,200,-0.1,+0.1);
88  se->registerHisto(fHTChi2[i]);
89  se->registerHisto(fHTDca2D[i]);
95  }
96 
97  TString sRname[4] = {"AllTracks","GoodTracks","GTracksM2R","GTracksExcess"};
98  for(int i=0; i!=4; ++i) {
99  fHRN[i] = new TH1F(Form("Tracks%d_N",i),Form("Number of %s",sRname[i].Data()),100,-0.5,1999.5);
100  fHRPt[i] = new TH1F(Form("Tracks%d_Pt",i),Form("%s Pt",sRname[i].Data()),nptbins,ptbins);
101  fHRPhi[i] = new TH1F(Form("Tracks%d_Phi",i),Form("%s Phi",sRname[i].Data()),100,0,TMath::TwoPi());
102  fHREta[i] = new TH1F(Form("Tracks%d_Eta",i),Form("%s Eta",sRname[i].Data()),100,-1.5,+1.5);
103  fHRTPCClus[i] = new TH2F(Form("Tracks%d_TPCClus",i),Form("%s TPCCLus;PT;TPCClusters",sRname[i].Data()),nptbins,ptbins,70,-0.5,69.5);
104  fHRChi2[i] = new TH2F(Form("Tracks%d_Chi2NDF",i),Form("%s Chi2NDF;PT;CHI2NDF",sRname[i].Data()),nptbins,ptbins,100,0,5);
105  fHRDca2D[i] = new TH2F( Form("Tracks%d_Dca2D",i),Form("%s Dca2D;PT;DCA2D",sRname[i].Data()),nptbins,ptbins,200,-0.003,+0.003);
106  se->registerHisto(fHRN[i]);
107  se->registerHisto(fHRPt[i]);
108  se->registerHisto(fHRPhi[i]);
109  se->registerHisto(fHREta[i]);
110  se->registerHisto(fHRTPCClus[i]);
111  se->registerHisto(fHRChi2[i]);
112  se->registerHisto(fHRDca2D[i]);
113  if(i<2) continue;
114  fHRPtResolution[i] = new TH2F(Form("Tracks%d_ResPt",i), Form("%s PtResolution;PT;REL DIFF",sRname[i].Data()), nptbins,ptbins,400,-1.5,+1.5);
115  fHRPtResolution2[i] = new TH2F(Form("Tracks%d_Res2Pt",i), Form("%s PtResolution2;PT;DIFF",sRname[i].Data()), nptbins,ptbins,400,-0.2,+0.2);
116  fHRPhiResolution[i] = new TH2F(Form("Tracks%d_ResPhi",i), Form("%s PhiResolution;PT;REL DIFF",sRname[i].Data()), nptbins,ptbins,200,-0.1,+0.1);
117  fHREtaResolution[i] = new TH2F(Form("Tracks%d_ResEta",i), Form("%s EtaResolution;PT;REL DIFF",sRname[i].Data()), nptbins,ptbins,200,-0.1,+0.1);
118  fHRNClustersContribution[i] = new TH2F(Form("Tracks%d_NClustersContri",i), Form("%s Number of Clusters Contribution",sRname[i].Data()), nptbins,ptbins,70,-0.5,69.5);
124  }
125 
126  fHNVertexes = new TH1F("Vertexes_N","Number of Vertexes",10,-0.5,9.5);
128  return 0;
129 }
130 
132  //std::cout << "TrackingPerformanceCheck [ENTER]" << std::endl;
133  fHNEvents->Fill(0);
134  fEmbedded.clear();
135  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
136  SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
137  SvtxClusterMap *clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
138  SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
139  if (!truthinfo) {
140  cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
141  exit(-1);
142  }
143  if (!trackmap) {
144  cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap" << endl;
145  exit(-1);
146  }
147  if (!vertexmap) {
148  cerr << PHWHERE << " ERROR: Can't find SvtxVertexMap" << endl;
149  exit(-1);
150  }
151 
152  // We need a method to distinguish between embedded and reconstructables. The current flag for embeded is an int with zero and one, so naturally this can be extended to hold values higher than two. However the current list cannot be modified, because the only accessor returns const_iterator, so we either change that in the base class (PH4TruthInfoContainer) or we clone the list (what we currently do here). Since this list only concerns embedded particles, it is going to be small
153  std::pair< std::map<int,int>::const_iterator, std::map<int,int>::const_iterator > rangeE = truthinfo->GetEmbeddedTrkIds();
154  for(std::map<int,int>::const_iterator iterE=rangeE.first; iterE!=rangeE.second; ++iterE) {
155  int tid = (*iterE).first;
156  int tem = (*iterE).second;
157  if (tem > 0) //embeded signal
158  fEmbedded.insert( std::make_pair(tid,tem) );
159  }
160  //std::cout << "TrackingPerformanceCheck ==> Embedded: " << fEmbedded.size() << std::endl;
161 
162  SvtxEvalStack svtxevalstack(topNode);
163  SvtxTruthEval *trutheval = svtxevalstack.get_truth_eval();
164  SvtxTrackEval *trackeval = svtxevalstack.get_track_eval();
165  //PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
166  PHG4TruthInfoContainer::Range range = truthinfo->GetParticleRange();
167  int nA=0;
168  int nB=0;
169  int nC=0;
170  int nD=0;
171  for(PHG4TruthInfoContainer::ConstIterator iter = range.first;
172  iter != range.second;
173  ++iter) {
174  PHG4Particle* g4particle = iter->second;
175  if(!g4particle) continue;
176  nA++; // all
177  float px = g4particle->get_px();
178  float py = g4particle->get_py();
179  float pz = g4particle->get_pz();
180  float pt = TMath::Sqrt( px*px + py*py );
181  float p = TMath::Sqrt( px*px + py*py + pz*pz );
182  float phi = TMath::Pi()+TMath::ATan2(-py,-px);
183  float eta = 1.e30;
184  if(p != TMath::Abs(pz)) eta = 0.5*TMath::Log((p+pz)/(p-pz));
185  fHTPt[0]->Fill(pt);
186  fHTPhi[0]->Fill(phi);
187  fHTEta[0]->Fill(eta);
188  if( iter->first < 0 ) continue;
189  nB++; // primaries
190  fHTPt[1]->Fill(pt);
191  fHTPhi[1]->Fill(phi);
192  fHTEta[1]->Fill(eta);
193  int id = g4particle->get_track_id();
194  std::map<int,int>::iterator embeddedflag = fEmbedded.find(id);
195  if(embeddedflag==fEmbedded.end()) continue; // not found
196  if((*embeddedflag).second==0) continue; // not embedded
197  //====> EMBEDDED
198  nC++; // embedded
199  //need a method to count how many layers where crossed by truth. Solution is to loop over g4hits and flag layers as hits belong to them.
200  std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(g4particle);
201  bool hit[70];
202  for(int i=0; i!=70; ++i) hit[i] = false;
203  for(std::set<PHG4Hit*>::iterator iter = g4hits.begin();
204  iter != g4hits.end();
205  ++iter) {
206  PHG4Hit *candidate = *iter;
207  unsigned int lyr = candidate->get_layer();
208  hit[lyr] = true;
209  }
210  int nhits = 0;
211  for(int i=fLayerTPCBegins; i!=70; ++i) if(hit[i]) nhits++; // counts only tpc hits
212  fHTPt[2]->Fill(pt);
213  fHTPhi[2]->Fill(phi);
214  fHTEta[2]->Fill(eta);
215  fHTNHits[2]->Fill( nhits );
216  //====> RECONSTRUCTABLES
217  if(nhits<fReconstructable_TPCHits) continue;
218  (*embeddedflag).second = 2;
219  nD++; // reconstructable
220  fHTPt[3]->Fill(pt);
221  fHTPhi[3]->Fill(phi);
222  fHTEta[3]->Fill(eta);
223  fHTNHits[3]->Fill( nhits );
224  SvtxTrack *best = trackeval->best_track_from(g4particle); // get the best track (max nclus) for this truth
225  if(!best) continue;
226  //====> MATCHED TO TRACK
227  float rpx = best->get_px();
228  float rpy = best->get_py();
229  float rpz = best->get_pz();
230  float rpt = TMath::Sqrt( rpx*rpx + rpy*rpy );
231  float rp = TMath::Sqrt( rpx*rpx + rpy*rpy + rpz*rpz );
232  float rphi = TMath::Pi()+TMath::ATan2(-rpy,-rpx);
233  float reta = 1.e30;
234  if(rp != TMath::Abs(rpz)) reta = 0.5*TMath::Log((rp+rpz)/(rp-rpz));
235  int nclusterscontrib = int( trackeval->get_nclusters_contribution(best,g4particle) ); // computes nclus
236  if(nclusterscontrib<fFair_ClustersContribution) continue;
237  float rdca2d = best->get_dca2d();
238  int ndf = int(best->get_ndf());
239  float chi2ndf = -1;
240  if(ndf!=0) chi2ndf = best->get_chisq() / ndf;
241  float ptreldiff = (rpt-pt)/pt;
242  float ptdiff = rpt-pt;
243  float phireldiff = (rphi-phi)/phi;
244  float etareldiff = (reta-eta)/eta;
245  fHTPt[4]->Fill(pt);
246  fHTPhi[4]->Fill(phi);
247  fHTEta[4]->Fill(eta);
248  fHTNHits[4]->Fill( nhits );
249  fHTChi2[4]->Fill(chi2ndf,ptreldiff);
250  fHTDca2D[4]->Fill(pt,rdca2d);
251  fHTPtResolution[4]->Fill(pt,ptreldiff);
252  fHTPtResolution2[4]->Fill(pt,ptdiff);
253  fHTPhiResolution[4]->Fill(pt,phireldiff);
254  fHTEtaResolution[4]->Fill(pt,etareldiff);
255  fHTNClustersContribution[4]->Fill(pt, nclusterscontrib );
256  }
257  fHTN[0]->Fill(nA);
258  fHTN[1]->Fill(nB);
259  fHTN[2]->Fill(nC);
260  fHTN[3]->Fill(nD);
261 
262  // loop over all reco particles
263  nA = 0;
264  nB = 0;
265  nC = 0;
266  nD = 0;
267  //std::cout << "TRACKS!" << std::endl;
268  for(SvtxTrackMap::Iter iter = trackmap->begin();
269  iter != trackmap->end(); ++iter) {
270  //===> ALL TRACKS
271  nA++;
272  SvtxTrack *track = iter->second;
273  if (!track) continue;
274  // counting tpcclusters
275  int tpcclus=0;
276  for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
277  iter != track->end_clusters(); ++iter) {
278  unsigned int cluster_id = *iter;
279  SvtxCluster* cluster = clustermap->get(cluster_id);
280  if(!cluster) continue;
281  int lyr = cluster->get_layer();
282  if(lyr>=fLayerTPCBegins) tpcclus++;
283  }
284  float rpx = track->get_px();
285  float rpy = track->get_py();
286  float rpz = track->get_pz();
287  float rpt = TMath::Sqrt( rpx*rpx + rpy*rpy );
288  float rp = TMath::Sqrt( rpx*rpx + rpy*rpy + rpz*rpz );
289  float rphi = TMath::Pi()+TMath::ATan2(-rpy,-rpx);
290  float reta = 1.e30;
291  if(rp != TMath::Abs(rpz)) reta = 0.5*TMath::Log((rp+rpz)/(rp-rpz));
292  float rdca2d = track->get_dca2d();
293  int ndf = int(track->get_ndf());
294  float chi2ndf = -1;
295  if(ndf!=0) chi2ndf = track->get_chisq() / ndf;
296  fHRPt[0]->Fill(rpt);
297  fHRPhi[0]->Fill(rphi);
298  fHREta[0]->Fill(reta);
299  fHRTPCClus[0]->Fill(rpt,tpcclus);
300  fHRChi2[0]->Fill(rpt,chi2ndf);
301  fHRDca2D[0]->Fill(rpt,rdca2d);
302  if(chi2ndf>fGTrack_Chi2NDF) continue;
303  if(tpcclus<fGTrack_TPCClusters) continue;
304  //===> SELECTED
305  nB++;
306  fHRPt[1]->Fill(rpt);
307  fHRPhi[1]->Fill(rphi);
308  fHREta[1]->Fill(reta);
309  fHRTPCClus[1]->Fill(rpt,tpcclus);
310  fHRChi2[1]->Fill(rpt,chi2ndf);
311  fHRDca2D[1]->Fill(rpt,rdca2d);
312  // check if matches to a MCTruth
313  PHG4Particle* g4particle = trackeval->max_truth_particle_by_nclusters(track); // get the most probable (max nclus) truth
314  if(!g4particle) continue;
315  //SvtxTrack *best = trackeval->best_track_from(g4particle); // get the best track (max nclus) for this truth
316  //if(!best) continue;
317  //std::cout << track->get_id() << " ==> " << g4particle->get_track_id() << " ==> " << best->get_id() << std::endl;
318  //if(best!=track) continue; // skippping low quality duplicates
319  // check if it was embedded
320  int id = g4particle->get_track_id();
321  std::map<int,int>::iterator embeddedflag = fEmbedded.find(id);
322  if(embeddedflag==fEmbedded.end()) continue; // not found
323  //std::cout << " " << track->get_id() << " ==> " << g4particle->get_track_id() << "(" << (*embeddedflag).second << ")" << " ==> " << best->get_id() << std::endl;
324  if((*embeddedflag).second==0) continue; // not embedded
325  // check if it was reconstructable
326  if((*embeddedflag).second==1) continue; // not reconstructable
327  if((*embeddedflag).second>1) (*embeddedflag).second++; // increase counter
328  //std::cout << " " << track->get_id() << " ==> " << g4particle->get_track_id() << "(" << (*embeddedflag).second << ")" << " ==> " << best->get_id() << std::endl;
329  //===> MATCHED TO EMBEDDED RECONSTRUCTABLE
330  nC++;
331  float px = g4particle->get_px();
332  float py = g4particle->get_py();
333  float pz = g4particle->get_pz();
334  float pt = TMath::Sqrt( px*px + py*py );
335  float p = TMath::Sqrt( px*px + py*py + pz*pz );
336  float phi = TMath::Pi()+TMath::ATan2(-py,-px);
337  float eta = 1.e30;
338  if(p != TMath::Abs(pz)) eta = 0.5*TMath::Log((p+pz)/(p-pz));
339  int nclusterscontrib = int( trackeval->get_nclusters_contribution(track,g4particle) ); // computes nclus
340  float ptreldiff = (rpt-pt)/pt;
341  float ptdiff = rpt-pt;
342  float phireldiff = (rphi-phi)/phi;
343  float etareldiff = (reta-eta)/eta;
344  fHRPt[2]->Fill(rpt);
345  fHRPhi[2]->Fill(rphi);
346  fHREta[2]->Fill(reta);
347  fHRTPCClus[2]->Fill(rpt,tpcclus);
348  fHRChi2[2]->Fill(rpt,chi2ndf);
349  fHRDca2D[2]->Fill(rpt,rdca2d);
350  fHRPtResolution[2]->Fill(pt,ptreldiff);
351  fHRPtResolution2[2]->Fill(pt,ptdiff);
352  fHRPhiResolution[2]->Fill(pt,phireldiff);
353  fHREtaResolution[2]->Fill(pt,etareldiff);
354  fHRNClustersContribution[2]->Fill(pt, nclusterscontrib);
355  if((*embeddedflag).second<4) continue; // 2->3 is okay, higher is excess
356  //===> FAKES (counters are the only ones meaningful)
357  nD++;
358  fHRPt[3]->Fill(rpt);
359  fHRPhi[3]->Fill(rphi);
360  fHREta[3]->Fill(reta);
361  fHRTPCClus[2]->Fill(rpt,tpcclus);
362  fHRChi2[3]->Fill(rpt,chi2ndf);
363  fHRDca2D[3]->Fill(rpt,rdca2d);
364  fHRPt[3]->Fill(rpt);
365  fHRPhi[3]->Fill(rphi);
366  fHREta[3]->Fill(reta);
367  fHRChi2[3]->Fill(rpt,chi2ndf);
368  fHRDca2D[3]->Fill(rpt,rdca2d);
369  fHRPtResolution[3]->Fill(pt,ptreldiff);
370  fHRPtResolution2[3]->Fill(pt,ptdiff);
371  fHRPhiResolution[3]->Fill(pt,phireldiff);
372  fHREtaResolution[3]->Fill(pt,etareldiff);
373  fHRNClustersContribution[3]->Fill(pt, nclusterscontrib );
374  }
375  fHRN[0]->Fill( nA );
376  fHRN[1]->Fill( nB );
377  fHRN[2]->Fill( nC );
378  fHRN[3]->Fill( nD );
379 
380  // get leading vertex
381  //SvtxVertexEval *vertexeval = svtxevalstack.get_vertex_eval();
382  //SvtxVertex* maxvertex = NULL;
383  //unsigned int maxtracks = 0;
384  nA=0;
385  for (SvtxVertexMap::Iter iter = vertexmap->begin();
386  iter != vertexmap->end();
387  ++iter) {
388  nA++;
389  //SvtxVertex* vertex = iter->second;
390  //if (vertex->size_tracks() > maxtracks) {
391  // maxvertex = vertex;
392  // maxtracks = vertex->size_tracks();
393  //}
394  }
395  fHNVertexes->Fill(nA);
396 
397  /*
398  if(maxvertex) {
399  PHG4VtxPoint* point = vertexeval->max_truth_point_by_ntracks(maxvertex);
400  if (point) {
401  _dx_vertex->Fill(maxvertex->get_x() - point->get_x());
402  _dy_vertex->Fill(maxvertex->get_y() - point->get_y());
403  _dz_vertex->Fill(maxvertex->get_z() - point->get_z());
404  }
405  }
406  */
407  //std::cout << "TrackingPerformanceCheck [EXIT]" << std::endl;
408  return 0;
409 }
410 
412  return 0;
413 }