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