Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SvtxSimPerformanceCheckReco.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SvtxSimPerformanceCheckReco.C
1 
2 #include "SvtxSimPerformanceCheckReco.h"
3 
4 #include <phool/getClass.h>
5 #include <fun4all/Fun4AllServer.h>
6 
7 #include <g4main/PHG4TruthInfoContainer.h>
8 #include <g4main/PHG4Particle.h>
9 #include <g4main/PHG4VtxPoint.h>
10 
11 #include <g4hough/SvtxVertexMap.h>
12 #include <g4hough/SvtxVertex.h>
13 #include <g4hough/SvtxTrackMap.h>
14 #include <g4hough/SvtxTrack.h>
15 
16 #include <g4eval/SvtxEvalStack.h>
17 #include <g4eval/SvtxTrackEval.h>
18 #include <g4eval/SvtxVertexEval.h>
19 #include <g4eval/SvtxTruthEval.h>
20 
21 #include <TH1D.h>
22 #include <TH2D.h>
23 
24 #include <iostream>
25 
26 using namespace std;
27 
28 SvtxSimPerformanceCheckReco::SvtxSimPerformanceCheckReco(const string &name)
29  : SubsysReco(name),
30  _event(0),
31  _nlayers(7),
32  _inner_layer_mask(0x7),
33  _truept_dptoverpt(NULL),
34  _truept_dca(NULL),
35  _truept_particles_leaving7Hits(NULL),
36  _truept_particles_recoWithExactHits(NULL),
37  _truept_particles_recoWithin1Hit(NULL),
38  _truept_particles_recoWithin2Hits(NULL),
39  _truept_particles_recoWithExactInnerHits(NULL),
40  _truept_particles_recoWithin1InnerHit(NULL),
41  _truept_particles_recoWithin2InnerHits(NULL),
42  _truept_particles_recoWithin3Percent(NULL),
43  _truept_particles_recoWithin4Percent(NULL),
44  _truept_particles_recoWithin5Percent(NULL),
45  _recopt_tracks_all(NULL),
46  _recopt_tracks_recoWithExactHits(NULL),
47  _recopt_tracks_recoWithin1Hit(NULL),
48  _recopt_tracks_recoWithin2Hits(NULL),
49  _recopt_tracks_recoWithExactInnerHits(NULL),
50  _recopt_tracks_recoWithin1InnerHit(NULL),
51  _recopt_tracks_recoWithin2InnerHits(NULL),
52  _recopt_tracks_recoWithin3Percent(NULL),
53  _recopt_tracks_recoWithin4Percent(NULL),
54  _recopt_tracks_recoWithin5Percent(NULL),
55  _recopt_quality(NULL),
56  _dx_vertex(NULL),
57  _dy_vertex(NULL),
58  _dz_vertex(NULL),
59  _truept_quality_particles_recoWithin4Percent(NULL),
60  _recopt_quality_tracks_all(NULL),
61  _recopt_quality_tracks_recoWithin4Percent(NULL) {
62 }
63 
64 int SvtxSimPerformanceCheckReco::Init(PHCompositeNode *topNode) {
65 
66  // register histograms
67  Fun4AllServer *se = Fun4AllServer::instance();
68 
69  _truept_dptoverpt = new TH2D("truept_dptoverpt",
70  "truept_dptoverpt",
71  40,0.0,40.0,
72  200,-0.5,0.5);
73 
74  _truept_dca = new TH2D("truept_dca",
75  "truept_dca",
76  20,0.0,10.0,
77  200,-0.1,0.1);
78 
79  _truept_particles_leaving7Hits = new TH1D("truept_particles_leaving7Hits",
80  "truept_particles_leaving7Hits",
81  20,0.0,10.0);
82 
83  _truept_particles_recoWithExactHits = new TH1D("truept_particles_recoWithExactHits",
84  "truept_particles_recoWithExactHits",
85  20,0.0,10.0);
86 
87  _truept_particles_recoWithin1Hit = new TH1D("truept_particles_recoWithin1Hit",
88  "truept_particles_recoWithin1Hit",
89  20,0.0,10.0);
90 
91  _truept_particles_recoWithin2Hits = new TH1D("truept_particles_recoWithin2Hits",
92  "truept_particles_recoWithin2Hits",
93  20,0.0,10.0);
94 
95  _truept_particles_recoWithExactInnerHits = new TH1D("truept_particles_recoWithExactInnerHits",
96  "truept_particles_recoWithExactInnerHits",
97  20,0.0,10.0);
98 
99  _truept_particles_recoWithin1InnerHit = new TH1D("truept_particles_recoWithin1InnerHit",
100  "truept_particles_recoWithin1InnerHit",
101  20,0.0,10.0);
102 
103  _truept_particles_recoWithin2InnerHits = new TH1D("truept_particles_recoWithin2InnerHits",
104  "truept_particles_recoWithin2InnerHits",
105  20,0.0,10.0);
106 
107 
108  _truept_particles_recoWithin3Percent = new TH1D("truept_particles_recoWithin3Percent",
109  "truept_particles_recoWithin3Percent",
110  20,0.0,10.0);
111 
112  _truept_particles_recoWithin4Percent = new TH1D("truept_particles_recoWithin4Percent",
113  "truept_particles_recoWithin4Percent",
114  20,0.0,10.0);
115 
116  _truept_particles_recoWithin5Percent = new TH1D("truept_particles_recoWithin5Percent",
117  "truept_particles_recoWithin5Percent",
118  20,0.0,10.0);
119 
120  _recopt_tracks_all = new TH1D("recopt_tracks_all",
121  "recopt_tracks_all",
122  20,0.0,10.0);
123 
124  _recopt_tracks_recoWithExactHits = new TH1D("recopt_tracks_recoWithExactHits",
125  "recopt_tracks_recoWithExactHits",
126  20,0.0,10.0);
127 
128  _recopt_tracks_recoWithin1Hit = new TH1D("recopt_tracks_recoWithin1Hit",
129  "recopt_tracks_recoWithin1Hit",
130  20,0.0,10.0);
131 
132  _recopt_tracks_recoWithin2Hits = new TH1D("recopt_tracks_recoWithin2Hits",
133  "recopt_tracks_recoWithin2Hits",
134  20,0.0,10.0);
135 
136  _recopt_tracks_recoWithExactInnerHits = new TH1D("recopt_tracks_recoWithExactInnerHits",
137  "recopt_tracks_recoWithExactInnerHits",
138  20,0.0,10.0);
139 
140  _recopt_tracks_recoWithin1InnerHit = new TH1D("recopt_tracks_recoWithin1InnerHit",
141  "recopt_tracks_recoWithin1InnerHit",
142  20,0.0,10.0);
143 
144  _recopt_tracks_recoWithin2InnerHits = new TH1D("recopt_tracks_recoWithin2InnerHits",
145  "recopt_tracks_recoWithin2InnerHits",
146  20,0.0,10.0);
147 
148  _recopt_tracks_recoWithin3Percent = new TH1D("recopt_tracks_recoWithin3Percent",
149  "recopt_tracks_recoWithin3Percent",
150  20,0.0,10.0);
151 
152  _recopt_tracks_recoWithin4Percent = new TH1D("recopt_tracks_recoWithin4Percent",
153  "recopt_tracks_recoWithin4Percent",
154  20,0.0,10.0);
155 
156  _recopt_tracks_recoWithin5Percent = new TH1D("recopt_tracks_recoWithin5Percent",
157  "recopt_tracks_recoWithin5Percent",
158  20,0.0,10.0);
159 
160  _recopt_quality = new TH2D("recopt_quality",
161  "recopt_quality",
162  20,0.0,10.0,
163  100,0.0,5.0);
164 
165  _dx_vertex = new TH1D("dx_vertex",
166  "dx_vertex",
167  200,-0.01,0.01);
168 
169  _dy_vertex = new TH1D("dy_vertex",
170  "dy_vertex",
171  200,-0.01,0.01);
172 
173  _dz_vertex = new TH1D("dz_vertex",
174  "dz_vertex",
175  200,-0.01,0.01);
176 
177  _truept_quality_particles_recoWithin4Percent = new TH2D("truept_quality_particles_recoWithin4Percent",
178  "truept_quality_particles_recoWithin4Percent",
179  20,0.0,10.0,
180  100,0.0,5.0);
181 
182  _recopt_quality_tracks_all = new TH2D("recopt_quality_tracks_all",
183  "recopt_quality_tracks_all",
184  20,0.0,10.0,
185  100,0.0,5.0);
186 
187  _recopt_quality_tracks_recoWithin4Percent = new TH2D("recopt_quality_tracks_recoWithin4Percent",
188  "recopt_quality_tracks_recoWithin4Percent",
189  20,0.0,10.0,
190  100,0.0,5.0);
191 
192 
193  se->registerHisto(_truept_dptoverpt);
194 
195  se->registerHisto(_truept_dca);
196 
197  se->registerHisto(_truept_particles_leaving7Hits);
198 
199  se->registerHisto(_truept_particles_recoWithExactHits);
200  se->registerHisto(_truept_particles_recoWithin1Hit);
201  se->registerHisto(_truept_particles_recoWithin2Hits);
202 
203  se->registerHisto(_truept_particles_recoWithExactInnerHits);
204  se->registerHisto(_truept_particles_recoWithin1InnerHit);
205  se->registerHisto(_truept_particles_recoWithin2InnerHits);
206 
207  se->registerHisto(_truept_particles_recoWithin3Percent);
208  se->registerHisto(_truept_particles_recoWithin4Percent);
209  se->registerHisto(_truept_particles_recoWithin5Percent);
210 
211  se->registerHisto(_recopt_tracks_all);
212 
213  se->registerHisto(_recopt_tracks_recoWithExactHits);
214  se->registerHisto(_recopt_tracks_recoWithin1Hit);
215  se->registerHisto(_recopt_tracks_recoWithin2Hits);
216 
217  se->registerHisto(_recopt_tracks_recoWithExactInnerHits);
218  se->registerHisto(_recopt_tracks_recoWithin1InnerHit);
219  se->registerHisto(_recopt_tracks_recoWithin2InnerHits);
220 
221  se->registerHisto(_recopt_tracks_recoWithin3Percent);
222  se->registerHisto(_recopt_tracks_recoWithin4Percent);
223  se->registerHisto(_recopt_tracks_recoWithin5Percent);
224 
225  se->registerHisto(_recopt_quality);
226 
227  se->registerHisto(_dx_vertex);
228  se->registerHisto(_dy_vertex);
229  se->registerHisto(_dz_vertex);
230 
231  se->registerHisto(_truept_quality_particles_recoWithin4Percent);
232  se->registerHisto(_recopt_quality_tracks_all);
233  se->registerHisto(_recopt_quality_tracks_recoWithin4Percent);
234 
235  return 0;
236 }
237 
238 int SvtxSimPerformanceCheckReco::process_event(PHCompositeNode *topNode) {
239 
240  ++_event;
241 
242  // need things off of the DST...
243  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
244  if (!truthinfo) {
245  cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
246  exit(-1);
247  }
248 
249  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
250  if (!trackmap) {
251  cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap" << endl;
252  exit(-1);
253  }
254 
255  SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
256  if (!vertexmap) {
257  cerr << PHWHERE << " ERROR: Can't find SvtxVertexMap" << endl;
258  exit(-1);
259  }
260 
261  // create SVTX eval stack
262  SvtxEvalStack svtxevalstack(topNode);
263 
264  SvtxVertexEval* vertexeval = svtxevalstack.get_vertex_eval();
265  SvtxTrackEval* trackeval = svtxevalstack.get_track_eval();
266  SvtxTruthEval* trutheval = svtxevalstack.get_truth_eval();
267 
268  // loop over all truth particles
269  PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
270  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
271  iter != range.second;
272  ++iter) {
273  PHG4Particle* g4particle = iter->second;
274  if (trutheval->get_embed(g4particle) <= 0) continue;
275 
276  std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(g4particle);
277  float ng4hits = g4hits.size();
278 
279  float truept = sqrt(pow(g4particle->get_px(),2)+pow(g4particle->get_py(),2));
280 
281  // examine truth particles that leave 7 detector hits
282  if (ng4hits == _nlayers) {
283  _truept_particles_leaving7Hits->Fill(truept);
284 
285  SvtxTrack* track = trackeval->best_track_from(g4particle);
286 
287  if (!track) {continue;}
288 
289  unsigned int nfromtruth = trackeval->get_nclusters_contribution(track,g4particle);
290  float recopt = track->get_pt();
291 
292  unsigned int ndiff = abs((int)nfromtruth-(int)_nlayers);
293  if (ndiff <= 2) {
294  _truept_particles_recoWithin2Hits->Fill(truept);
295  }
296  if (ndiff <= 1) {
297  _truept_particles_recoWithin1Hit->Fill(truept);
298  }
299  if (ndiff == 0) {
300  _truept_particles_recoWithExactHits->Fill(truept);
301  }
302 
303  unsigned int layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track,g4particle);
304  unsigned int innerhits = (layersfromtruth & _inner_layer_mask);
305  unsigned int ninnerhitsfromtruth = 0;
306  unsigned int ninnerlayers = 0;
307  for (unsigned int i = 0; i < 32; ++i) {
308  ninnerhitsfromtruth += (0x1 & (innerhits >> i));
309  ninnerlayers += (0x1 & (_inner_layer_mask >> i));
310  }
311 
312  ndiff = abs((int)ninnerhitsfromtruth-(int)ninnerlayers);
313  if (ndiff <= 2) {
314  _truept_particles_recoWithin2InnerHits->Fill(truept);
315  }
316  if (ndiff <= 1) {
317  _truept_particles_recoWithin1InnerHit->Fill(truept);
318  }
319  if (ndiff == 0) {
320  _truept_particles_recoWithExactInnerHits->Fill(truept);
321  }
322 
323  float diff = fabs(recopt-truept)/truept;
324  if (diff < 0.05) {
325  _truept_particles_recoWithin5Percent->Fill(truept);
326  }
327  if (diff < 0.04) {
328  _truept_particles_recoWithin4Percent->Fill(truept);
329  _truept_quality_particles_recoWithin4Percent->Fill(truept,track->get_quality());
330  }
331  if (diff < 0.03) {
332  _truept_particles_recoWithin3Percent->Fill(truept);
333  }
334  }
335  }
336 
337  // loop over all reco particles
338  for (SvtxTrackMap::Iter iter = trackmap->begin();
339  iter != trackmap->end();
340  ++iter) {
341 
342  SvtxTrack* track = iter->second;
343  float recopt = track->get_pt();
344 
345  PHG4Particle* g4particle = trackeval->max_truth_particle_by_nclusters(track);
346  float truept = sqrt(pow(g4particle->get_px(),2)+pow(g4particle->get_py(),2));
347 
348  if (trutheval->get_embed(g4particle) > 0) {
349  // embedded results (quality or preformance measures)
350  _truept_dptoverpt->Fill(truept,(recopt-truept)/truept);
351  _truept_dca->Fill(truept,track->get_dca2d());
352 
353  _recopt_quality->Fill(recopt,track->get_quality());
354  } else {
355  // non-embedded results (purity measures)
356  _recopt_tracks_all->Fill(recopt);
357  _recopt_quality_tracks_all->Fill(recopt,track->get_quality());
358 
359  unsigned int nfromtruth = trackeval->get_nclusters_contribution(track,g4particle);
360 
361  unsigned int ndiff = abs((int)nfromtruth-(int)_nlayers);
362  if (ndiff <= 2) {
363  _recopt_tracks_recoWithin2Hits->Fill(recopt);
364  }
365  if (ndiff <= 1) {
366  _recopt_tracks_recoWithin1Hit->Fill(recopt);
367  }
368  if (ndiff == 0) {
369  _recopt_tracks_recoWithExactHits->Fill(recopt);
370  }
371 
372  unsigned int layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track,g4particle);
373  unsigned int innerhits = (layersfromtruth & _inner_layer_mask);
374  unsigned int ninnerhitsfromtruth = 0;
375  unsigned int ninnerlayers = 0;
376  for (unsigned int i = 0; i < 32; ++i) {
377  ninnerhitsfromtruth += (0x1 & (innerhits >> i));
378  ninnerlayers += (0x1 & (_inner_layer_mask >> i));
379  }
380 
381  ndiff = abs((int)ninnerhitsfromtruth-(int)ninnerlayers);
382  if (ndiff <= 2) {
383  _recopt_tracks_recoWithin2InnerHits->Fill(recopt);
384  }
385  if (ndiff <= 1) {
386  _recopt_tracks_recoWithin1InnerHit->Fill(recopt);
387  }
388  if (ndiff == 0) {
389  _recopt_tracks_recoWithExactInnerHits->Fill(recopt);
390  }
391 
392  float diff = fabs(recopt-truept)/truept;
393  if (diff < 0.05) {
394  _recopt_tracks_recoWithin5Percent->Fill(recopt);
395  }
396  if (diff < 0.04) {
397  _recopt_tracks_recoWithin4Percent->Fill(recopt);
398  _recopt_quality_tracks_recoWithin4Percent->Fill(recopt,track->get_quality());
399  }
400  if (diff < 0.03) {
401  _recopt_tracks_recoWithin3Percent->Fill(recopt);
402  }
403  }
404  }
405 
406  // get leading vertex
407  SvtxVertex* maxvertex = NULL;
408  unsigned int maxtracks = 0;
409  for (SvtxVertexMap::Iter iter = vertexmap->begin();
410  iter != vertexmap->end();
411  ++iter) {
412  SvtxVertex* vertex = iter->second;
413 
414  if (vertex->size_tracks() > maxtracks) {
415  maxvertex = vertex;
416  maxtracks = vertex->size_tracks();
417  }
418  }
419 
420  if (maxvertex) {
421  PHG4VtxPoint* point = vertexeval->max_truth_point_by_ntracks(maxvertex);
422 
423  if (point) {
424  _dx_vertex->Fill(maxvertex->get_x() - point->get_x());
425  _dy_vertex->Fill(maxvertex->get_y() - point->get_y());
426  _dz_vertex->Fill(maxvertex->get_z() - point->get_z());
427  }
428  }
429 
430  return 0;
431 }
432 
433 int SvtxSimPerformanceCheckReco::End(PHCompositeNode *topNode) {
434  return 0;
435 }