Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
psTOFTimezeroEval.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file psTOFTimezeroEval.C
1 #include <g4eval/SvtxEvalStack.h>
2 
3 #include <phool/PHCompositeNode.h>
4 #include <fun4all/Fun4AllReturnCodes.h>
5 #include <phool/getClass.h>
6 
7 #include <g4hough/SvtxVertexMap.h>
8 #include <g4hough/SvtxVertex.h>
9 #include <g4hough/SvtxTrackMap.h>
10 #include <g4hough/SvtxTrack.h>
11 #include <g4hough/SvtxClusterMap.h>
12 #include <g4hough/SvtxCluster.h>
13 #include <g4hough/SvtxHitMap.h>
14 #include <g4hough/SvtxHit.h>
15 
16 #include <g4main/PHG4Hit.h>
17 #include <g4main/PHG4Particle.h>
18 #include <g4main/PHG4VtxPoint.h>
19 #include <g4main/PHG4TruthInfoContainer.h>
20 
21 #include <g4detectors/PHG4Cell.h>
22 
23 #include <TFile.h>
24 #include <TNtuple.h>
25 
26 #include <iostream>
27 #include <set>
28 #include <cmath>
29 #include <cassert>
30 
31 #include "../../psTOF/psTOFTimezeroEval/psTOFTimezeroEval.h"
32 
33 using namespace std;
34 
35 psTOFTimezeroEval::psTOFTimezeroEval(const string &name, const string &filename) :
36  SubsysReco("psTOFTimezeroEval"),
37  _ievent(0),
38  _svtxevalstack(nullptr),
39  _strict(false),
40  _errors(0),
41  _do_vertex_eval(true),
42  _do_gpoint_eval(true),
43  _do_g4hit_eval(true),
44  _do_hit_eval(true),
45  _do_cluster_eval(true),
46  _do_gtrack_eval(true),
47  _do_track_eval(true),
48  _scan_for_embedded(false),
49  _ntp_vertex(nullptr),
50  _ntp_gpoint(nullptr),
51  _ntp_g4hit(nullptr),
52  _ntp_hit(nullptr),
53  _ntp_cluster(nullptr),
54  _ntp_gtrack(nullptr),
55  _ntp_track(nullptr),
56  _filename(filename),
57  _tfile(nullptr) {
58  verbosity = 0;
59 }
60 
61 int psTOFTimezeroEval::Init(PHCompositeNode *topNode) {
62 
63  _ievent = 0;
64 
65  _tfile = new TFile(_filename.c_str(), "RECREATE");
66 
67  if (_do_vertex_eval) _ntp_vertex = new TNtuple("ntp_vertex","vertex => max truth",
68  "event:vx:vy:vz:ntracks:"
69  "gvx:gvy:gvz:gvt:gntracks:"
70  "nfromtruth");
71 
72  if (_do_gpoint_eval) _ntp_gpoint = new TNtuple("ntp_gpoint","g4point => best vertex",
73  "event:gvx:gvy:gvz:gvt:gntracks:"
74  "vx:vy:vz:ntracks:"
75  "nfromtruth");
76 
77  if (_do_g4hit_eval) _ntp_g4hit = new TNtuple("ntp_g4hit","g4hit => best svtxcluster",
78  "event:g4hitID:gx:gy:gz:gt:gedep:"
79  "glayer:gtrackID:gflavor:"
80  "gpx:gpy:gpz:gvx:gvy:gvz:"
81  "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
82  "gembed:gprimary:nclusters:"
83  "clusID:x:y:z:e:adc:layer:size:"
84  "phisize:zsize:efromtruth");
85 
86  if (_do_hit_eval) _ntp_hit = new TNtuple("ntp_hit","svtxhit => max truth",
87  "event:hitID:e:adc:layer:"
88  "cellID:ecell:"
89  "g4hitID:gedep:gx:gy:gz:gt:"
90  "gtrackID:gflavor:"
91  "gpx:gpy:gpz:gvx:gvy:gvz:"
92  "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
93  "gembed:gprimary:efromtruth");
94 
95  if (_do_cluster_eval) _ntp_cluster = new TNtuple("ntp_cluster","svtxcluster => max truth",
96  "event:hitID:x:y:z:ex:ey:ez:ephi:"
97  "e:adc:layer:size:phisize:"
98  "zsize:trackID:g4hitID:gx:"
99  "gy:gz:gt:gtrackID:gflavor:"
100  "gpx:gpy:gpz:gvx:gvy:gvz:"
101  "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
102  "gembed:gprimary:efromtruth");
103 
104  if (_do_gtrack_eval) _ntp_gtrack = new TNtuple("ntp_gtrack","g4particle => best svtxtrack",
105  "event:gtrackID:gflavor:gnhits:"
106  "gpx:gpy:gpz:"
107  "gvx:gvy:gvz:gvt:"
108  "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
109  "gembed:gprimary:"
110  "trackID:px:py:pz:charge:quality:chisq:ndf:nhits:layers:"
111  "dca2d:dca2dsigma:pcax:pcay:pcaz:nfromtruth:layersfromtruth:"
112  "nmaps:nintt:ntpc");
113 
114  if (_do_track_eval) _ntp_track = new TNtuple("ntp_track","svtxtrack => max truth",
115  "event:trackID:px:py:pz:charge:"
116  "quality:chisq:ndf:nhits:layers:"
117  "dca2d:dca2dsigma:pcax:pcay:pcaz:"
118  "presdphi:presdeta:prese3x3:prese:"
119  "cemcdphi:cemcdeta:cemce3x3:cemce:"
120  "hcalindphi:hcalindeta:hcaline3x3:hcaline:"
121  "hcaloutdphi:hcaloutdeta:hcaloute3x3:hcaloute:"
122  "gtrackID:gflavor:gnhits:"
123  "gpx:gpy:gpz:"
124  "gvx:gvy:gvz:gvt:"
125  "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
126  "gembed:gprimary:nfromtruth:layersfromtruth:"
127  "nmaps:nintt:ntpc");
128 
129  return Fun4AllReturnCodes::EVENT_OK;
130 }
131 
132 int psTOFTimezeroEval::InitRun(PHCompositeNode *topNode) {
133  return Fun4AllReturnCodes::EVENT_OK;
134 }
135 
136 int psTOFTimezeroEval::process_event(PHCompositeNode *topNode) {
137 
138  if ((verbosity > 0)&&(_ievent%100==0)) {
139  cout << "psTOFTimezeroEval::process_event - Event = " << _ievent << endl;
140  }
141 
142  if (!_svtxevalstack) {
143  _svtxevalstack = new SvtxEvalStack(topNode);
144  _svtxevalstack->set_strict(_strict);
145  _svtxevalstack->set_verbosity(verbosity+1);
146  } else {
147  _svtxevalstack->next_event(topNode);
148  }
149 
150  //-----------------------------------
151  // print what is coming into the code
152  //-----------------------------------
153 
154  printInputInfo(topNode);
155 
156  //---------------------------
157  // fill the Evaluator NTuples
158  //---------------------------
159 
160  fillOutputNtuples(topNode);
161 
162  //--------------------------------------------------
163  // Print out the ancestry information for this event
164  //--------------------------------------------------
165 
166  printOutputInfo(topNode);
167 
168  ++_ievent;
169  return Fun4AllReturnCodes::EVENT_OK;
170 }
171 
172 int psTOFTimezeroEval::End(PHCompositeNode *topNode) {
173 
174  _tfile->cd();
175 
176  if (_ntp_vertex) _ntp_vertex->Write();
177  if (_ntp_gpoint) _ntp_gpoint->Write();
178  if (_ntp_g4hit) _ntp_g4hit->Write();
179  if (_ntp_hit) _ntp_hit->Write();
180  if (_ntp_cluster) _ntp_cluster->Write();
181  if (_ntp_gtrack) _ntp_gtrack->Write();
182  if (_ntp_track) _ntp_track->Write();
183 
184  _tfile->Close();
185 
186  delete _tfile;
187 
188  if (verbosity > 0) {
189  cout << "========================= psTOFTimezeroEval::End() ============================" << endl;
190  cout << " " << _ievent << " events of output written to: " << _filename << endl;
191  cout << "===========================================================================" << endl;
192  }
193 
194  _errors += _svtxevalstack->get_errors();
195 
196  if (verbosity > -1) {
197  if ((_errors > 0)||(verbosity > 0)) {
198  cout << "psTOFTimezeroEval::End() - Error Count: " << _errors << endl;
199  }
200  }
201 
202  delete _svtxevalstack;
203 
204  return Fun4AllReturnCodes::EVENT_OK;
205 }
206 
207 void psTOFTimezeroEval::printInputInfo(PHCompositeNode *topNode) {
208 
209  if (verbosity > 1) cout << "psTOFTimezeroEval::printInputInfo() entered" << endl;
210 
211  if (verbosity > 3) {
212 
213  // event information
214  cout << endl;
215  cout << PHWHERE << " INPUT FOR EVENT " << _ievent << endl;
216 
217  cout << endl;
218  cout << "---PHG4HITS-------------" << endl;
219  _svtxevalstack->get_truth_eval()->set_strict(_strict);
220  std::set<PHG4Hit*> g4hits = _svtxevalstack->get_truth_eval()->all_truth_hits();
221  unsigned int ig4hit = 0;
222  for(std::set<PHG4Hit*>::iterator iter = g4hits.begin();
223  iter != g4hits.end();
224  ++iter) {
225  PHG4Hit *g4hit = *iter;
226  cout << ig4hit << " of " << g4hits.size();
227  cout << ": PHG4Hit: " << endl;
228  g4hit->identify();
229  ++ig4hit;
230  }
231 
232  cout << "---SVTXCLUSTERS-------------" << endl;
233  SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
234  if (clustermap) {
235  unsigned int icluster = 0;
236  for (SvtxClusterMap::Iter iter = clustermap->begin();
237  iter != clustermap->end();
238  ++iter) {
239  SvtxCluster* cluster = iter->second;
240  cout << icluster << " of " << clustermap->size();
241  cout << ": SvtxCluster: " << endl;
242  cluster->identify();
243  ++icluster;
244  }
245  }
246 
247  cout << "---SVXTRACKS-------------" << endl;
248  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
249  if (trackmap) {
250  unsigned int itrack = 0;
251  for (SvtxTrackMap::Iter iter = trackmap->begin();
252  iter != trackmap->end();
253  ++iter) {
254  cout << itrack << " of " << trackmap->size();
255  SvtxTrack *track = iter->second;
256  cout << " : SvtxTrack:" << endl;
257  track->identify();
258  cout << endl;
259  }
260  }
261 
262  cout << "---SVXVERTEXES-------------" << endl;
263  SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
264  if (vertexmap) {
265  unsigned int ivertex = 0;
266  for (SvtxVertexMap::Iter iter = vertexmap->begin();
267  iter != vertexmap->end();
268  ++iter) {
269  cout << ivertex << " of " << vertexmap->size();
270  SvtxVertex *vertex = iter->second;
271  cout << " : SvtxVertex:" << endl;
272  vertex->identify();
273  cout << endl;
274  }
275  }
276  }
277 
278  return;
279 }
280 
281 void psTOFTimezeroEval::printOutputInfo(PHCompositeNode *topNode) {
282 
283  if (verbosity > 1) cout << "psTOFTimezeroEval::printOutputInfo() entered" << endl;
284 
285  //==========================================
286  // print out some useful stuff for debugging
287  //==========================================
288 
289  if (verbosity > 0) {
290 
291  SvtxTrackEval* trackeval = _svtxevalstack->get_track_eval();
292  SvtxClusterEval* clustereval = _svtxevalstack->get_cluster_eval();
293  SvtxTruthEval* trutheval = _svtxevalstack->get_truth_eval();
294 
295  // event information
296  cout << endl;
297  cout << PHWHERE << " NEW OUTPUT FOR EVENT " << _ievent << endl;
298  cout << endl;
299 
300  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
301 
302  PHG4VtxPoint *gvertex = truthinfo->GetPrimaryVtx( truthinfo->GetPrimaryVertexIndex() );
303  float gvx = gvertex->get_x();
304  float gvy = gvertex->get_y();
305  float gvz = gvertex->get_z();
306 
307  float vx = NAN;
308  float vy = NAN;
309  float vz = NAN;
310 
311  SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
312  if (vertexmap) {
313  if (!vertexmap->empty()) {
314  SvtxVertex* vertex = (vertexmap->begin()->second);
315 
316  vx = vertex->get_x();
317  vy = vertex->get_y();
318  vz = vertex->get_z();
319  }
320  }
321 
322  cout << "===Vertex Reconstruction=======================" << endl;
323  cout << "vtrue = (" << gvx << "," << gvy << "," << gvz << ") => vreco = (" << vx << "," << vy << "," << vz << ")" << endl;
324  cout << endl;
325 
326  cout << "===Tracking Summary============================" << endl;
327  unsigned int ng4hits[100] = {0};
328  std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits();
329  for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
330  iter != g4hits.end();
331  ++iter) {
332  PHG4Hit *g4hit = *iter;
333  ++ng4hits[g4hit->get_layer()];
334  }
335 
336  SvtxHitMap* hitmap = findNode::getClass<SvtxHitMap>(topNode,"SvtxHitMap");
337  unsigned int nhits[100] = {0};
338  if (hitmap) {
339  for (SvtxHitMap::Iter iter = hitmap->begin();
340  iter != hitmap->end();
341  ++iter) {
342  SvtxHit* hit = iter->second;
343  ++nhits[hit->get_layer()];
344  }
345  }
346 
347  SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
348  unsigned int nclusters[100] = {0};
349  if (clustermap) {
350  for (SvtxClusterMap::Iter iter = clustermap->begin();
351  iter != clustermap->end();
352  ++iter) {
353  SvtxCluster* cluster = iter->second;
354  ++nclusters[cluster->get_layer()];
355  }
356  }
357 
358  for (unsigned int ilayer = 0; ilayer < 100; ++ilayer) {
359  cout << "layer " << ilayer << ": nG4hits = " << ng4hits[ilayer]
360  << " => nHits = " << nhits[ilayer]
361  << " => nClusters = " << nclusters[ilayer] << endl;
362  }
363 
364  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
365 
366  cout << "nGtracks = " << std::distance(truthinfo->GetPrimaryParticleRange().first,
367  truthinfo->GetPrimaryParticleRange().second);
368  cout << " => nTracks = ";
369  if (trackmap) cout << trackmap->size() << endl;
370  else cout << 0 << endl;
371 
372  // cluster wise information
373  if (verbosity > 1) {
374 
375  for(std::set<PHG4Hit*>::iterator iter = g4hits.begin();
376  iter != g4hits.end();
377  ++iter) {
378  PHG4Hit *g4hit = *iter;
379 
380  cout << endl;
381  cout << "===PHG4Hit===================================" << endl;
382  cout << " PHG4Hit: "; g4hit->identify();
383 
384  std::set<SvtxCluster*> clusters = clustereval->all_clusters_from(g4hit);
385 
386  for (std::set<SvtxCluster*>::iterator jter = clusters.begin();
387  jter != clusters.end();
388  ++jter) {
389  SvtxCluster *cluster = *jter;
390  cout << "===Created-SvtxCluster================" << endl;
391  cout << "SvtxCluster: "; cluster->identify();
392  }
393  }
394 
395  PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
396  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
397  iter != range.second;
398  ++iter) {
399 
400  PHG4Particle *particle = iter->second;
401 
402  // track-wise information
403  cout << endl;
404 
405  cout << "=== Gtrack ===================================================" << endl;
406  cout << " PHG4Particle id = " << particle->get_track_id() << endl;
407  particle->identify();
408  cout << " ptrue = (";
409  cout.width(5); cout << particle->get_px();
410  cout << ",";
411  cout.width(5); cout << particle->get_py();
412  cout << ",";
413  cout.width(5); cout << particle->get_pz();
414  cout << ")" << endl;
415 
416  cout << " vtrue = (";
417  cout.width(5); cout << truthinfo->GetVtx(particle->get_vtx_id())->get_x();
418  cout << ",";
419  cout.width(5); cout << truthinfo->GetVtx(particle->get_vtx_id())->get_y();
420  cout << ",";
421  cout.width(5); cout << truthinfo->GetVtx(particle->get_vtx_id())->get_z();
422  cout << ")" << endl;
423 
424  cout << " pt = " << sqrt(pow(particle->get_px(),2)+pow(particle->get_py(),2)) << endl;
425  cout << " phi = " << atan2(particle->get_py(),particle->get_px()) << endl;
426  cout << " eta = " << asinh(particle->get_pz()/sqrt(pow(particle->get_px(),2)+pow(particle->get_py(),2))) << endl;
427 
428  cout << " embed flag = " << truthinfo->isEmbeded(particle->get_track_id()) << endl;
429 
430  cout << " ---Associated-PHG4Hits-----------------------------------------" << endl;
431  std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(particle);
432  for(std::set<PHG4Hit*>::iterator jter = g4hits.begin();
433  jter != g4hits.end();
434  ++jter) {
435  PHG4Hit *g4hit = *jter;
436 
437  float x = 0.5*(g4hit->get_x(0)+g4hit->get_x(1));
438  float y = 0.5*(g4hit->get_y(0)+g4hit->get_y(1));
439  float z = 0.5*(g4hit->get_z(0)+g4hit->get_z(1));
440 
441  cout << " #" << g4hit->get_hit_id() << " xtrue = (";
442  cout.width(5); cout << x;
443  cout << ",";
444  cout.width(5); cout << y;
445  cout << ",";
446  cout.width(5); cout << z;
447  cout << ")";
448 
449  std::set<SvtxCluster*> clusters = clustereval->all_clusters_from(g4hit);
450  for (std::set<SvtxCluster*>::iterator kter = clusters.begin();
451  kter != clusters.end();
452  ++kter) {
453 
454  SvtxCluster *cluster = *kter;
455 
456  float x = cluster->get_x();
457  float y = cluster->get_y();
458  float z = cluster->get_z();
459 
460  cout << " => #" << cluster->get_id();
461  cout << " xreco = (";
462  cout.width(5); cout << x;
463  cout << ",";
464  cout.width(5); cout << y;
465  cout << ",";
466  cout.width(5); cout << z;
467  cout << ")";
468  }
469 
470  cout << endl;
471  }
472 
473  if (trackmap&&clustermap) {
474 
475  std::set<SvtxTrack*> tracks = trackeval->all_tracks_from(particle);
476  for (std::set<SvtxTrack*>::iterator jter = tracks.begin();
477  jter != tracks.end();
478  ++jter) {
479 
480  SvtxTrack *track = *jter;
481 
482  float px = track->get_px();
483  float py = track->get_py();
484  float pz = track->get_pz();
485 
486  cout << "===Created-SvtxTrack==========================================" << endl;
487  cout << " SvtxTrack id = " << track->get_id() << endl;
488  cout << " preco = (";
489  cout.width(5); cout << px;
490  cout << ",";
491  cout.width(5); cout << py;
492  cout << ",";
493  cout.width(5); cout << pz;
494  cout << ")" << endl;
495  cout << " quality = " << track->get_quality() << endl;
496  cout << " nfromtruth = " << trackeval->get_nclusters_contribution(track,particle) << endl;
497 
498  cout << " ---Associated-SvtxClusters-to-PHG4Hits-------------------------" << endl;
499 
500  for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
501  iter != track->end_clusters();
502  ++iter) {
503  unsigned int cluster_id = *iter;
504  SvtxCluster* cluster = clustermap->get(cluster_id);
505 
506  float x = cluster->get_x();
507  float y = cluster->get_y();
508  float z = cluster->get_z();
509 
510  cout << " #" << cluster->get_id() << " xreco = (";
511  cout.width(5); cout << x;
512  cout << ",";
513  cout.width(5); cout << y;
514  cout << ",";
515  cout.width(5); cout << z;
516  cout << ") =>";
517 
518  PHG4Hit* g4hit = clustereval->max_truth_hit_by_energy(cluster);
519  if ((g4hit) && (g4hit->get_trkid() == particle->get_track_id())) {
520 
521  x = 0.5*(g4hit->get_x(0)+g4hit->get_x(1));
522  y = 0.5*(g4hit->get_y(0)+g4hit->get_y(1));
523  z = 0.5*(g4hit->get_z(0)+g4hit->get_z(1));
524 
525  cout << " #" << g4hit->get_hit_id()
526  << " xtrue = (";
527  cout.width(5); cout << x;
528  cout << ",";
529  cout.width(5); cout << y;
530  cout << ",";
531  cout.width(5); cout << z;
532  cout << ") => Gtrack id = " << g4hit->get_trkid();
533  } else {
534  cout << " noise hit";
535  }
536  }
537 
538  cout << endl;
539  }
540  }
541  }
542  }
543 
544  cout << endl;
545 
546  } // if verbosity
547 
548  return;
549 }
550 
551 void psTOFTimezeroEval::fillOutputNtuples(PHCompositeNode *topNode) {
552 
553  if (verbosity > 1) cout << "psTOFTimezeroEval::fillOutputNtuples() entered" << endl;
554 
555  SvtxVertexEval* vertexeval = _svtxevalstack->get_vertex_eval();
556  SvtxTrackEval* trackeval = _svtxevalstack->get_track_eval();
557  SvtxClusterEval* clustereval = _svtxevalstack->get_cluster_eval();
558  SvtxHitEval* hiteval = _svtxevalstack->get_hit_eval();
559  SvtxTruthEval* trutheval = _svtxevalstack->get_truth_eval();
560 
561  //-----------------------
562  // fill the Vertex NTuple
563  //-----------------------
564 
565  if (_ntp_vertex) {
566  //cout << "Filling ntp_vertex " << endl;
567  SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
568  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
569  if (vertexmap && truthinfo) {
570  for (SvtxVertexMap::Iter iter = vertexmap->begin();
571  iter != vertexmap->end();
572  ++iter) {
573  SvtxVertex* vertex = iter->second;
574  PHG4VtxPoint* point = vertexeval->max_truth_point_by_ntracks(vertex);
575 
576  float vx = vertex->get_x();
577  float vy = vertex->get_y();
578  float vz = vertex->get_z();
579  float ntracks = vertex->size_tracks();
580 
581  float gvx = NAN;
582  float gvy = NAN;
583  float gvz = NAN;
584  float gvt = NAN;
585  float gntracks = truthinfo->GetNumPrimaryVertexParticles();
586  float nfromtruth = NAN;
587 
588  if (point) {
589  gvx = point->get_x();
590  gvy = point->get_y();
591  gvz = point->get_z();
592  gvt = point->get_t();
593  gntracks = truthinfo->GetNumPrimaryVertexParticles();
594  nfromtruth = vertexeval->get_ntracks_contribution(vertex,point);
595  }
596 
597  float vertex_data[11] = {(float) _ievent,
598  vx,
599  vy,
600  vz,
601  ntracks,
602  gvx,
603  gvy,
604  gvz,
605  gvt,
606  gntracks,
607  nfromtruth
608  };
609 
610  /*
611  cout << "vertex: "
612  << " ievent " << vertex_data[0]
613  << " vx " << vertex_data[1]
614  << " vy " << vertex_data[2]
615  << " vz " << vertex_data[3]
616  << endl;
617  */
618 
619  _ntp_vertex->Fill(vertex_data);
620  }
621  }
622  }
623 
624  //-----------------------
625  // fill the gpoint NTuple
626  //-----------------------
627 
628  if (_ntp_gpoint) {
629  //cout << "Filling ntp_gpoint " << endl;
630  SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
631  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
632 
633  if (vertexmap && truthinfo) {
634 
635  PHG4VtxPoint* point = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
636 
637  if (point) {
638 
639  SvtxVertex* vertex = vertexeval->best_vertex_from(point);
640 
641  float gvx = point->get_x();
642  float gvy = point->get_y();
643  float gvz = point->get_z();
644  float gvt = point->get_t();
645  float gntracks = truthinfo->GetNumPrimaryVertexParticles();
646  float vx = NAN;
647  float vy = NAN;
648  float vz = NAN;
649  float ntracks = NAN;
650  float nfromtruth = NAN;
651 
652  if (vertex) {
653  vx = vertex->get_x();
654  vy = vertex->get_y();
655  vz = vertex->get_z();
656  ntracks = vertex->size_tracks();
657  nfromtruth = vertexeval->get_ntracks_contribution(vertex,point);
658  }
659 
660  float gpoint_data[11] = {(float) _ievent,
661  gvx,
662  gvy,
663  gvz,
664  gvt,
665  gntracks,
666  vx,
667  vy,
668  vz,
669  ntracks,
670  nfromtruth
671  };
672 
673  _ntp_gpoint->Fill(gpoint_data);
674  }
675  }
676  }
677 
678  //---------------------
679  // fill the G4hit NTuple
680  //---------------------
681 
682  if (_ntp_g4hit) {
683  //cout << "Filling ntp_g4hit " << endl;
684  std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits();
685  for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
686  iter != g4hits.end();
687  ++iter) {
688 
689  PHG4Hit *g4hit = *iter;
690  PHG4Particle *g4particle = trutheval->get_particle(g4hit);
691 
692  float g4hitID = g4hit->get_hit_id();
693  float gx = g4hit->get_avg_x();
694  float gy = g4hit->get_avg_y();
695  float gz = g4hit->get_avg_z();
696  float gt = g4hit->get_avg_t();
697  float gedep = g4hit->get_edep();
698  float glayer = g4hit->get_layer();
699 
700  float gtrackID = g4hit->get_trkid();
701 
702  float gflavor = NAN;
703  float gpx = NAN;
704  float gpy = NAN;
705  float gpz = NAN;
706 
707  float gvx = NAN;
708  float gvy = NAN;
709  float gvz = NAN;
710 
711  float gembed = NAN;
712  float gprimary = NAN;
713 
714  float gfpx = 0.;
715  float gfpy = 0.;
716  float gfpz = 0.;
717  float gfx = 0.;
718  float gfy = 0.;
719  float gfz = 0.;
720 
721  if (g4particle) {
722 
723  if (_scan_for_embedded) {
724  if (trutheval->get_embed(g4particle) <= 0) continue;
725  }
726 
727  gflavor = g4particle->get_pid();
728  gpx = g4particle->get_px();
729  gpy = g4particle->get_py();
730  gpz = g4particle->get_pz();
731 
732  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
733 
734  if (vtx) {
735  gvx = vtx->get_x();
736  gvy = vtx->get_y();
737  gvz = vtx->get_z();
738  }
739 
740  PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
741  if (outerhit) {
742  gfpx = outerhit->get_px(1);
743  gfpy = outerhit->get_py(1);
744  gfpz = outerhit->get_pz(1);
745  gfx = outerhit->get_x(1);
746  gfy = outerhit->get_y(1);
747  gfz = outerhit->get_z(1);
748  }
749 
750  gembed = trutheval->get_embed(g4particle);
751  gprimary = trutheval->is_primary(g4particle);
752  } // if (g4particle)
753 
754  std::set<SvtxCluster*> clusters = clustereval->all_clusters_from(g4hit);
755  float nclusters = clusters.size();
756 
757  // best cluster reco'd
758  SvtxCluster* cluster = clustereval->best_cluster_from(g4hit);
759 
760  float clusID = NAN;
761  float x = NAN;
762  float y = NAN;
763  float z = NAN;
764  float e = NAN;
765  float adc = NAN;
766  float layer = NAN;
767  float size = NAN;
768  float phisize = NAN;
769  float zsize = NAN;
770  float efromtruth = NAN;
771 
772  if (cluster) {
773  clusID = cluster->get_id();
774  x = cluster->get_x();
775  y = cluster->get_y();
776  z = cluster->get_z();
777  e = cluster->get_e();
778  adc = cluster->get_adc();
779  layer = cluster->get_layer();
780  size = cluster->size_hits();
781  phisize = cluster->get_phi_size();
782  zsize = cluster->get_z_size();
783  if (g4particle) {
784  efromtruth = clustereval->get_energy_contribution(cluster,g4particle);
785  }
786  }
787 
788  float g4hit_data[36] = {(float) _ievent,
789  g4hitID,
790  gx,
791  gy,
792  gz,
793  gt,
794  gedep,
795  glayer,
796  gtrackID,
797  gflavor,
798  gpx,
799  gpy,
800  gpz,
801  gvx,
802  gvy,
803  gvz,
804  gfpx,
805  gfpy,
806  gfpz,
807  gfx,
808  gfy,
809  gfz,
810  gembed,
811  gprimary,
812  nclusters,
813  clusID,
814  x,
815  y,
816  z,
817  e,
818  adc,
819  layer,
820  size,
821  phisize,
822  zsize,
823  efromtruth
824  };
825 
826  _ntp_g4hit->Fill(g4hit_data);
827  }
828  }
829 
830  //--------------------
831  // fill the Hit NTuple
832  //--------------------
833 
834  if (_ntp_hit) {
835  //cout << "Filling ntp_hit " << endl;
836  // need things off of the DST...
837  SvtxHitMap* hitmap = findNode::getClass<SvtxHitMap>(topNode,"SvtxHitMap");
838  if (hitmap) {
839 
840  for (SvtxHitMap::Iter iter = hitmap->begin();
841  iter != hitmap->end();
842  ++iter) {
843 
844  SvtxHit* hit = iter->second;
845  PHG4Hit* g4hit = hiteval->max_truth_hit_by_energy(hit);
846  PHG4Cell* g4cell = hiteval->get_cell(hit);
847  PHG4Particle* g4particle = trutheval->get_particle(g4hit);
848 
849  float event = _ievent;
850  float hitID = hit->get_id();
851  float e = hit->get_e();
852  float adc = hit->get_adc();
853  float layer = hit->get_layer();
854  float cellID = hit->get_cellid();
855  float ecell = g4cell->get_edep();
856 
857  float g4hitID = NAN;
858  float gedep = NAN;
859  float gx = NAN;
860  float gy = NAN;
861  float gz = NAN;
862  float gt = NAN;
863  float gtrackID = NAN;
864  float gflavor = NAN;
865  float gpx = NAN;
866  float gpy = NAN;
867  float gpz = NAN;
868  float gvx = NAN;
869  float gvy = NAN;
870  float gvz = NAN;
871  float gfpx = NAN;
872  float gfpy = NAN;
873  float gfpz = NAN;
874  float gfx = NAN;
875  float gfy = NAN;
876  float gfz = NAN;
877  float gembed = NAN;
878  float gprimary = NAN;
879 
880  float efromtruth = NAN;
881 
882  if (g4hit) {
883  g4hitID = g4hit->get_hit_id();
884  gedep = g4hit->get_edep();
885  gx = g4hit->get_avg_x();
886  gy = g4hit->get_avg_y();
887  gz = g4hit->get_avg_z();
888  gt = g4hit->get_avg_t();
889 
890  if (g4particle) {
891 
892  if (_scan_for_embedded) {
893  if (trutheval->get_embed(g4particle) <= 0) continue;
894  }
895 
896  gtrackID = g4particle->get_track_id();
897  gflavor = g4particle->get_pid();
898  gpx = g4particle->get_px();
899  gpy = g4particle->get_py();
900  gpz = g4particle->get_pz();
901 
902  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
903 
904  if (vtx) {
905  gvx = vtx->get_x();
906  gvy = vtx->get_y();
907  gvz = vtx->get_z();
908  }
909 
910  PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
911  if (outerhit) {
912  gfpx = outerhit->get_px(1);
913  gfpy = outerhit->get_py(1);
914  gfpz = outerhit->get_pz(1);
915  gfx = outerhit->get_x(1);
916  gfy = outerhit->get_y(1);
917  gfz = outerhit->get_z(1);
918  }
919  gembed = trutheval->get_embed(g4particle);
920  gprimary = trutheval->is_primary(g4particle);
921  } // if (g4particle){
922  }
923 
924  if (g4particle) {
925  efromtruth = hiteval->get_energy_contribution(hit,g4particle);
926  }
927 
928  float hit_data[33] = {
929  event,
930  hitID,
931  e,
932  adc,
933  layer,
934  cellID,
935  ecell,
936  g4hitID,
937  gedep,
938  gx,
939  gy,
940  gz,
941  gt,
942  gtrackID,
943  gflavor,
944  gpx,
945  gpy,
946  gpz,
947  gvx,
948  gvy,
949  gvz,
950  gfpx,
951  gfpy,
952  gfpz,
953  gfx,
954  gfy,
955  gfz,
956  gembed,
957  gprimary,
958  efromtruth
959  };
960 
961  _ntp_hit->Fill(hit_data);
962  }
963  }
964  }
965 
966  //------------------------
967  // fill the Cluster NTuple
968  //------------------------
969 
970  //cout << "check for ntp_cluster" << endl;
971 
972  if (_ntp_cluster && !_scan_for_embedded) {
973  //cout << "Filling ntp_cluster 1 " << endl;
974  // need things off of the DST...
975  SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
976  if (clustermap) {
977 
978  for (SvtxClusterMap::Iter iter = clustermap->begin();
979  iter != clustermap->end();
980  ++iter) {
981 
982  SvtxCluster* cluster = iter->second;
983  SvtxTrack* track = trackeval->best_track_from(cluster);
984 
985  PHG4Hit *g4hit = clustereval->max_truth_hit_by_energy(cluster);
986  PHG4Particle *g4particle = trutheval->get_particle(g4hit);
987 
988  float hitID = cluster->get_id();
989  float x = cluster->get_x();
990  float y = cluster->get_y();
991  float z = cluster->get_z();
992 
993  float ex = sqrt(cluster->get_error(0,0));
994  float ey = sqrt(cluster->get_error(1,1));
995  float ez = cluster->get_z_error();
996 
997  float ephi = cluster->get_rphi_error();
998 
999  float e = cluster->get_e();
1000  float adc = cluster->get_adc();
1001  float layer = cluster->get_layer();
1002  float size = cluster->size_hits();
1003  float phisize = cluster->get_phi_size();
1004  float zsize = cluster->get_z_size();
1005 
1006  float trackID = NAN;
1007  if (track) trackID = track->get_id();
1008 
1009  float g4hitID = NAN;
1010  float gx = NAN;
1011  float gy = NAN;
1012  float gz = NAN;
1013  float gt = NAN;
1014  float gtrackID = NAN;
1015  float gflavor = NAN;
1016  float gpx = NAN;
1017  float gpy = NAN;
1018  float gpz = NAN;
1019  float gvx = NAN;
1020  float gvy = NAN;
1021  float gvz = NAN;
1022  float gfpx = NAN;
1023  float gfpy = NAN;
1024  float gfpz = NAN;
1025  float gfx = NAN;
1026  float gfy = NAN;
1027  float gfz = NAN;
1028  float gembed = NAN;
1029  float gprimary = NAN;
1030 
1031  float efromtruth = NAN;
1032 
1033  if (g4hit) {
1034  g4hitID = g4hit->get_hit_id();
1035  gx = g4hit->get_avg_x();
1036  gy = g4hit->get_avg_y();
1037  gz = g4hit->get_avg_z();
1038  gt = g4hit->get_avg_t();
1039 
1040  if (g4particle) {
1041 
1042  gtrackID = g4particle->get_track_id();
1043  gflavor = g4particle->get_pid();
1044  gpx = g4particle->get_px();
1045  gpy = g4particle->get_py();
1046  gpz = g4particle->get_pz();
1047 
1048  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1049  if (vtx) {
1050  gvx = vtx->get_x();
1051  gvy = vtx->get_y();
1052  gvz = vtx->get_z();
1053  }
1054 
1055  PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
1056  if (outerhit) {
1057  gfpx = outerhit->get_px(1);
1058  gfpy = outerhit->get_py(1);
1059  gfpz = outerhit->get_pz(1);
1060  gfx = outerhit->get_x(1);
1061  gfy = outerhit->get_y(1);
1062  gfz = outerhit->get_z(1);
1063  }
1064 
1065  gembed = trutheval->get_embed(g4particle);
1066  gprimary = trutheval->is_primary(g4particle);
1067  } // if (g4particle){
1068  } // if (g4hit) {
1069 
1070  if (g4particle){
1071  efromtruth = clustereval->get_energy_contribution(cluster,g4particle);
1072  }
1073 
1074  float cluster_data[38] = {(float) _ievent,
1075  hitID,
1076  x,
1077  y,
1078  z,
1079  ex,
1080  ey,
1081  ez,
1082  ephi,
1083  e,
1084  adc,
1085  layer,
1086  size,
1087  phisize,
1088  zsize,
1089  trackID,
1090  g4hitID,
1091  gx,
1092  gy,
1093  gz,
1094  gt,
1095  gtrackID,
1096  gflavor,
1097  gpx,
1098  gpy,
1099  gpz,
1100  gvx,
1101  gvy,
1102  gvz,
1103  gfpx,
1104  gfpy,
1105  gfpz,
1106  gfx,
1107  gfy,
1108  gfz,
1109  gembed,
1110  gprimary,
1111  efromtruth};
1112 
1113  _ntp_cluster->Fill(cluster_data);
1114  }
1115  }
1116 
1117  } else if (_ntp_cluster && _scan_for_embedded) {
1118 
1119  //cout << "Filling ntp_cluster 2 " << endl;
1120 
1121  // if only scanning embedded signals, loop over all the tracks from
1122  // embedded particles and report all of their clusters, including those
1123  // from other sources (noise hits on the embedded track)
1124 
1125  // need things off of the DST...
1126  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
1127  SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
1128  if (trackmap) {
1129 
1130  for (SvtxTrackMap::Iter iter = trackmap->begin();
1131  iter != trackmap->end();
1132  ++iter) {
1133 
1134  SvtxTrack* track = iter->second;
1135  PHG4Particle* truth = trackeval->max_truth_particle_by_nclusters(track);
1136  if (truth) {
1137  if (trutheval->get_embed(truth) <= 0) continue;
1138  }
1139 
1140  for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
1141  iter != track->end_clusters();
1142  ++iter) {
1143 
1144  unsigned int cluster_id = *iter;
1145  SvtxCluster* cluster = clustermap->get(cluster_id);
1146 
1147  PHG4Hit *g4hit = clustereval->max_truth_hit_by_energy(cluster);
1148  PHG4Particle *g4particle = trutheval->get_particle(g4hit);
1149 
1150  float hitID = cluster->get_id();
1151  float x = cluster->get_x();
1152  float y = cluster->get_y();
1153  float z = cluster->get_z();
1154 
1155  float ex = sqrt(cluster->get_error(0,0));
1156  float ey = sqrt(cluster->get_error(1,1));
1157  float ez = cluster->get_z_error();
1158 
1159  float ephi = cluster->get_rphi_error();
1160 
1161  float e = cluster->get_e();
1162  float adc = cluster->get_adc();
1163  float layer = cluster->get_layer();
1164  float size = cluster->size_hits();
1165  float phisize = cluster->get_phi_size();
1166  float zsize = cluster->get_z_size();
1167 
1168  float trackID = NAN;
1169  if (track) trackID = track->get_id();
1170 
1171  float g4hitID = NAN;
1172  float gx = NAN;
1173  float gy = NAN;
1174  float gz = NAN;
1175  float gt = NAN;
1176  float gtrackID = NAN;
1177  float gflavor = NAN;
1178  float gpx = NAN;
1179  float gpy = NAN;
1180  float gpz = NAN;
1181  float gvx = NAN;
1182  float gvy = NAN;
1183  float gvz = NAN;
1184  float gfpx = NAN;
1185  float gfpy = NAN;
1186  float gfpz = NAN;
1187  float gfx = NAN;
1188  float gfy = NAN;
1189  float gfz = NAN;
1190  float gembed = NAN;
1191  float gprimary = NAN;
1192 
1193  float efromtruth = NAN;
1194 
1195  if (g4hit) {
1196  g4hitID = g4hit->get_hit_id();
1197  gx = g4hit->get_avg_x();
1198  gy = g4hit->get_avg_y();
1199  gz = g4hit->get_avg_z();
1200  gt = g4hit->get_avg_t();
1201 
1202  if (g4particle) {
1203 
1204  gtrackID = g4particle->get_track_id();
1205  gflavor = g4particle->get_pid();
1206  gpx = g4particle->get_px();
1207  gpy = g4particle->get_py();
1208  gpz = g4particle->get_pz();
1209 
1210  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1211  if (vtx) {
1212  gvx = vtx->get_x();
1213  gvy = vtx->get_y();
1214  gvz = vtx->get_z();
1215  }
1216 
1217  PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
1218  if (outerhit) {
1219  gfpx = outerhit->get_px(1);
1220  gfpy = outerhit->get_py(1);
1221  gfpz = outerhit->get_pz(1);
1222  gfx = outerhit->get_x(1);
1223  gfy = outerhit->get_y(1);
1224  gfz = outerhit->get_z(1);
1225  }
1226 
1227  gembed = trutheval->get_embed(g4particle);
1228  gprimary = trutheval->is_primary(g4particle);
1229  } // if (g4particle){
1230  } // if (g4hit) {
1231 
1232  if (g4particle){
1233  efromtruth = clustereval->get_energy_contribution(cluster,g4particle);
1234  }
1235 
1236  float cluster_data[38] = {(float) _ievent,
1237  hitID,
1238  x,
1239  y,
1240  z,
1241  ex,
1242  ey,
1243  ez,
1244  ephi,
1245  e,
1246  adc,
1247  layer,
1248  size,
1249  phisize,
1250  zsize,
1251  trackID,
1252  g4hitID,
1253  gx,
1254  gy,
1255  gz,
1256  gt,
1257  gtrackID,
1258  gflavor,
1259  gpx,
1260  gpy,
1261  gpz,
1262  gvx,
1263  gvy,
1264  gvz,
1265  gfpx,
1266  gfpy,
1267  gfpz,
1268  gfx,
1269  gfy,
1270  gfz,
1271  gembed,
1272  gprimary,
1273  efromtruth};
1274 
1275  _ntp_cluster->Fill(cluster_data);
1276  }
1277  }
1278  }
1279  }
1280 
1281  //------------------------
1282  // fill the Gtrack NTuple
1283  //------------------------
1284 
1285  // need things off of the DST...
1286 
1287  //cout << "check for ntp_gtrack" << endl;
1288 
1289  if (_ntp_gtrack) {
1290  //cout << "Filling ntp_gtrack " << endl;
1291 
1292  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
1293  SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
1294  if (truthinfo) {
1295 
1296  PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
1297  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
1298  iter != range.second;
1299  ++iter) {
1300 
1301  PHG4Particle* g4particle = iter->second;
1302 
1303  if (_scan_for_embedded) {
1304  if (trutheval->get_embed(g4particle) <= 0) continue;
1305  }
1306 
1307  float gtrackID = g4particle->get_track_id();
1308  float gflavor = g4particle->get_pid();
1309 
1310  std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(g4particle);
1311  float ng4hits = g4hits.size();
1312  float gpx = g4particle->get_px();
1313  float gpy = g4particle->get_py();
1314  float gpz = g4particle->get_pz();
1315 
1316  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1317  float gvx = vtx->get_x();
1318  float gvy = vtx->get_y();
1319  float gvz = vtx->get_z();
1320  float gvt = vtx->get_t();
1321 
1322  float gfpx = 0.;
1323  float gfpy = 0.;
1324  float gfpz = 0.;
1325  float gfx = 0.;
1326  float gfy = 0.;
1327  float gfz = 0.;
1328 
1329  PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
1330  if (outerhit) {
1331  gfpx = outerhit->get_px(1);
1332  gfpy = outerhit->get_py(1);
1333  gfpz = outerhit->get_pz(1);
1334  gfx = outerhit->get_x(1);
1335  gfy = outerhit->get_y(1);
1336  gfz = outerhit->get_z(1);
1337  }
1338 
1339  float gembed = trutheval->get_embed(g4particle);
1340  float gprimary = trutheval->is_primary(g4particle);
1341 
1342  SvtxTrack* track = trackeval->best_track_from(g4particle);
1343 
1344  float trackID = NAN;
1345  float charge = NAN;
1346  float quality = NAN;
1347  float chisq = NAN;
1348  float ndf = NAN;
1349  float nhits = NAN;
1350  unsigned int layers = 0x0;
1351  float dca2d = NAN;
1352  float dca2dsigma = NAN;
1353  float px = NAN;
1354  float py = NAN;
1355  float pz = NAN;
1356  float pcax = NAN;
1357  float pcay = NAN;
1358  float pcaz = NAN;
1359 
1360  float nfromtruth = NAN;
1361  float layersfromtruth = NAN;
1362 
1363  float nmaps = 0;
1364  float nintt = 0;
1365  float ntpc = 0;
1366 
1367  if (track) {
1368  trackID = track->get_id();
1369  charge = track->get_charge();
1370  quality = track->get_quality();
1371  chisq = track->get_chisq();
1372  ndf = track->get_ndf();
1373  nhits = track->size_clusters();
1374 
1375  for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
1376  iter != track->end_clusters();
1377  ++iter) {
1378  unsigned int cluster_id = *iter;
1379  SvtxCluster* cluster = clustermap->get(cluster_id);
1380  unsigned int layer = cluster->get_layer();
1381  if (layer < 32) layers |= (0x1 << layer);
1382 
1383 
1384  if(layer >= 0 and layer <= 2) ++nmaps;
1385  if(layer >= 3 and layer <= 6) ++nintt;
1386  if(layer >= 7 and layer <= 66) ++ntpc;
1387  }
1388 
1389  dca2d = track->get_dca2d();
1390  dca2dsigma = track->get_dca2d_error();
1391  px = track->get_px();
1392  py = track->get_py();
1393  pz = track->get_pz();
1394  pcax = track->get_x();
1395  pcay = track->get_y();
1396  pcaz = track->get_z();
1397 
1398  nfromtruth = trackeval->get_nclusters_contribution(track,g4particle);
1399  layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track,g4particle);
1400  }
1401 
1402  float gtrack_data[39] = {(float) _ievent,
1403  gtrackID,
1404  gflavor,
1405  ng4hits,
1406  gpx,
1407  gpy,
1408  gpz,
1409  gvx,
1410  gvy,
1411  gvz,
1412  gvt,
1413  gfpx,
1414  gfpy,
1415  gfpz,
1416  gfx,
1417  gfy,
1418  gfz,
1419  gembed,
1420  gprimary,
1421  trackID,
1422  px,
1423  py,
1424  pz,
1425  charge,
1426  quality,
1427  chisq,
1428  ndf,
1429  nhits,
1430  (float) layers,
1431  dca2d,
1432  dca2dsigma,
1433  pcax,
1434  pcay,
1435  pcaz,
1436  nfromtruth,
1437  layersfromtruth,
1438  nmaps,
1439  nintt,
1440  ntpc
1441  };
1442 
1443  /*
1444  cout << " ievent " << _ievent
1445  << " gtrackID " << gtrackID
1446  << " gflavor " << gflavor
1447  << " ng4hits " << ng4hits
1448  << endl;
1449  */
1450 
1451  _ntp_gtrack->Fill(gtrack_data);
1452 
1453  }
1454  }
1455  }
1456 
1457  //------------------------
1458  // fill the Track NTuple
1459  //------------------------
1460 
1461 
1462 
1463  if (_ntp_track) {
1464  //cout << "Filling ntp_track " << endl;
1465 
1466  // need things off of the DST...
1467  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
1468  SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
1469  if (trackmap) {
1470 
1471  for (SvtxTrackMap::Iter iter = trackmap->begin();
1472  iter != trackmap->end();
1473  ++iter) {
1474 
1475  SvtxTrack* track = iter->second;
1476 
1477  float trackID = track->get_id();
1478  float charge = track->get_charge();
1479  float quality = track->get_quality();
1480  float chisq = track->get_chisq();
1481  float ndf = track->get_ndf();
1482  float nhits = track->size_clusters();
1483 
1484  float nmaps = 0;
1485  float nintt = 0;
1486  float ntpc = 0;
1487 
1488  unsigned int layers = 0x0;
1489  for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
1490  iter != track->end_clusters();
1491  ++iter) {
1492  unsigned int cluster_id = *iter;
1493  SvtxCluster* cluster = clustermap->get(cluster_id);
1494  unsigned int layer = cluster->get_layer();
1495  if (layer < 31) layers |= (0x1 << layer);
1496 
1497  if(layer >= 0 and layer <= 2) ++nmaps;
1498  if(layer >= 3 and layer <= 6) ++nintt;
1499  if(layer >= 7 and layer <= 66) ++ntpc;
1500  }
1501 
1502  float dca2d = track->get_dca2d();
1503  float dca2dsigma = track->get_dca2d_error();
1504  float px = track->get_px();
1505  float py = track->get_py();
1506  float pz = track->get_pz();
1507  float pcax = track->get_x();
1508  float pcay = track->get_y();
1509  float pcaz = track->get_z();
1510 
1511  float presdphi = track->get_cal_dphi(SvtxTrack::PRES);
1512  float presdeta = track->get_cal_deta(SvtxTrack::PRES);
1513  float prese3x3 = track->get_cal_energy_3x3(SvtxTrack::PRES);
1514  float prese = track->get_cal_cluster_e(SvtxTrack::PRES);
1515 
1516  float cemcdphi = track->get_cal_dphi(SvtxTrack::CEMC);
1517  float cemcdeta = track->get_cal_deta(SvtxTrack::CEMC);
1518  float cemce3x3 = track->get_cal_energy_3x3(SvtxTrack::CEMC);
1519  float cemce = track->get_cal_cluster_e(SvtxTrack::CEMC);
1520 
1521  float hcalindphi = track->get_cal_dphi(SvtxTrack::HCALIN);
1522  float hcalindeta = track->get_cal_deta(SvtxTrack::HCALIN);
1523  float hcaline3x3 = track->get_cal_energy_3x3(SvtxTrack::HCALIN);
1524  float hcaline = track->get_cal_cluster_e(SvtxTrack::HCALIN);
1525 
1526  float hcaloutdphi = track->get_cal_dphi(SvtxTrack::HCALOUT);
1527  float hcaloutdeta = track->get_cal_deta(SvtxTrack::HCALOUT);
1528  float hcaloute3x3 = track->get_cal_energy_3x3(SvtxTrack::HCALOUT);
1529  float hcaloute = track->get_cal_cluster_e(SvtxTrack::HCALOUT);
1530 
1531  float gtrackID = NAN;
1532  float gflavor = NAN;
1533  float ng4hits = NAN;
1534  float gpx = NAN;
1535  float gpy = NAN;
1536  float gpz = NAN;
1537  float gvx = NAN;
1538  float gvy = NAN;
1539  float gvz = NAN;
1540  float gvt = NAN;
1541  float gfpx = NAN;
1542  float gfpy = NAN;
1543  float gfpz = NAN;
1544  float gfx = NAN;
1545  float gfy = NAN;
1546  float gfz = NAN;
1547  float gembed = NAN;
1548  float gprimary = NAN;
1549 
1550  float nfromtruth = NAN;
1551  float layersfromtruth = NAN;
1552 
1553  PHG4Particle* g4particle = trackeval->max_truth_particle_by_nclusters(track);
1554 
1555  if (g4particle) {
1556 
1557  if (_scan_for_embedded) {
1558  if (trutheval->get_embed(g4particle) <= 0) continue;
1559  }
1560 
1561  gtrackID = g4particle->get_track_id();
1562  gflavor = g4particle->get_pid();
1563 
1564  std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(g4particle);
1565  ng4hits = g4hits.size();
1566  gpx = g4particle->get_px();
1567  gpy = g4particle->get_py();
1568  gpz = g4particle->get_pz();
1569 
1570  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1571  gvx = vtx->get_x();
1572  gvy = vtx->get_y();
1573  gvz = vtx->get_z();
1574  gvt = vtx->get_t();
1575 
1576  PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
1577  if (outerhit) {
1578  gfpx = outerhit->get_px(1);
1579  gfpy = outerhit->get_py(1);
1580  gfpz = outerhit->get_pz(1);
1581  gfx = outerhit->get_x(1);
1582  gfy = outerhit->get_y(1);
1583  gfz = outerhit->get_z(1);
1584  }
1585  gembed = trutheval->get_embed(g4particle);
1586  gprimary = trutheval->is_primary(g4particle);
1587 
1588  nfromtruth = trackeval->get_nclusters_contribution(track,g4particle);
1589  layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track,g4particle);
1590  }
1591 
1592  float track_data[55] = {(float) _ievent,
1593  trackID,
1594  px,
1595  py,
1596  pz,
1597  charge,
1598  quality,
1599  chisq,
1600  ndf,
1601  nhits,
1602  (float) layers,
1603  dca2d,
1604  dca2dsigma,
1605  pcax,
1606  pcay,
1607  pcaz,
1608  presdphi,
1609  presdeta,
1610  prese3x3,
1611  prese,
1612  cemcdphi,
1613  cemcdeta,
1614  cemce3x3,
1615  cemce,
1616  hcalindphi,
1617  hcalindeta,
1618  hcaline3x3,
1619  hcaline,
1620  hcaloutdphi,
1621  hcaloutdeta,
1622  hcaloute3x3,
1623  hcaloute,
1624  gtrackID,
1625  gflavor,
1626  ng4hits,
1627  gpx,
1628  gpy,
1629  gpz,
1630  gvx,
1631  gvy,
1632  gvz,
1633  gvt,
1634  gfpx,
1635  gfpy,
1636  gfpz,
1637  gfx,
1638  gfy,
1639  gfz,
1640  gembed,
1641  gprimary,
1642  nfromtruth,
1643  layersfromtruth,
1644  nmaps,
1645  nintt,
1646  ntpc
1647  };
1648 
1649  /*
1650  cout << "ievent " << _ievent
1651  << " trackID " << trackID
1652  << " nhits " << nhits
1653  << " px " << px
1654  << " py " << py
1655  << " pz " << pz
1656  << " gembed " << gembed
1657  << " gprimary " << gprimary
1658  << endl;
1659  */
1660  _ntp_track->Fill(track_data);
1661  }
1662  }
1663  }
1664 
1665  return;
1666 }