ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SvtxEvaluator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SvtxEvaluator.cc
1 #include "SvtxEvaluator.h"
2 
3 #include "SvtxEvalStack.h"
4 
5 #include "SvtxClusterEval.h"
6 #include "SvtxHitEval.h"
7 #include "SvtxTrackEval.h"
8 #include "SvtxTruthEval.h"
9 #include "SvtxVertexEval.h"
10 
11 #include <trackbase/TrkrCluster.h>
12 #include <trackbase/TrkrHit.h>
13 #include <trackbase/TrkrDefs.h>
14 #include <trackbase/ActsTrackingGeometry.h>
15 #include <trackbase/ActsSurfaceMaps.h>
16 
17 #include <tpc/TpcDefs.h>
18 
19 #include <trackbase/TrkrClusterContainer.h>
20 #include <trackbase/TrkrHitSet.h>
21 #include <trackbase/TrkrHitSetContainer.h>
22 #include <trackbase/TrkrClusterHitAssoc.h>
23 #include <trackbase/TrkrClusterIterationMapv1.h>
24 
25 #include <trackbase_historic/SvtxTrack.h>
26 #include <trackbase_historic/SvtxTrackMap.h>
27 #include <trackbase_historic/SvtxVertex.h>
28 #include <trackbase_historic/SvtxVertexMap.h>
29 #include <trackbase_historic/ActsTransformations.h>
30 
31 #include <g4main/PHG4Hit.h>
32 #include <g4main/PHG4Particle.h>
33 #include <g4main/PHG4TruthInfoContainer.h>
34 #include <g4main/PHG4VtxPoint.h>
35 
36 #include <g4detectors/PHG4CylinderCellGeom.h>
37 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
38 
39 #include <fun4all/Fun4AllReturnCodes.h>
40 #include <fun4all/SubsysReco.h>
41 
42 #include <phool/PHTimer.h>
43 #include <phool/getClass.h>
44 #include <phool/phool.h>
45 #include <phool/recoConsts.h>
46 
47 #include <TFile.h>
48 #include <TNtuple.h>
49 #include <TVector3.h>
50 
51 #include <cmath>
52 #include <iostream>
53 #include <iomanip>
54 #include <iterator>
55 #include <map>
56 #include <memory> // for shared_ptr
57 #include <set> // for _Rb_tree_cons...
58 #include <utility>
59 #include <vector>
60 
61 using namespace std;
62 
63 SvtxEvaluator::SvtxEvaluator(const string& /*name*/, const string& filename, const string& trackmapname,
64  unsigned int nlayers_maps,
65  unsigned int nlayers_intt,
66  unsigned int nlayers_tpc,
67  unsigned int nlayers_mms)
68  : SubsysReco("SvtxEvaluator")
69  , _ievent(0)
70  , _iseed(0)
71  , m_fSeed(NAN)
72  , _svtxevalstack(nullptr)
73  , _strict(false)
74  , _errors(0)
75  , _do_info_eval(true)
76  , _do_vertex_eval(true)
77  , _do_gpoint_eval(true)
78  , _do_g4hit_eval(true)
79  , _do_hit_eval(true)
80  , _do_cluster_eval(true)
81  , _do_g4cluster_eval(true)
82  , _do_gtrack_eval(true)
83  , _do_track_eval(true)
84  , _do_gseed_eval(false)
85  , _do_track_match(true)
86  , _do_eval_light(true)
87  , _do_vtx_eval_light(true)
88  , _scan_for_embedded(false)
89  , _scan_for_primaries(false)
90  , _nlayers_maps(nlayers_maps)
91  , _nlayers_intt(nlayers_intt)
92  , _nlayers_tpc(nlayers_tpc)
93  , _nlayers_mms(nlayers_mms)
94  , _ntp_info(nullptr)
95  , _ntp_vertex(nullptr)
96  , _ntp_gpoint(nullptr)
97  , _ntp_g4hit(nullptr)
98  , _ntp_hit(nullptr)
99  , _ntp_cluster(nullptr)
100  , _ntp_g4cluster(nullptr)
101  , _ntp_gtrack(nullptr)
102  , _ntp_track(nullptr)
103  , _ntp_gseed(nullptr)
104  , _filename(filename)
105  , _trackmapname(trackmapname)
106  , _tfile(nullptr)
107  , _timer(nullptr)
108 {
109 }
110 
111 SvtxEvaluator::~SvtxEvaluator()
112 {
113  delete _timer;
114 }
115 
116 int SvtxEvaluator::Init(PHCompositeNode* /*topNode*/)
117 {
118  _ievent = 0;
119 
120  _tfile = new TFile(_filename.c_str(), "RECREATE");
121  _tfile->SetCompressionLevel(0);
122  if (_do_info_eval) _ntp_info = new TNtuple("ntp_info", "event info",
123  "event:seed:"
124  "occ11:occ116:occ21:occ216:occ31:occ316:"
125  "gntrkall:gntrkprim:ntrk:"
126  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
127 
128  if (_do_vertex_eval) _ntp_vertex = new TNtuple("ntp_vertex", "vertex => max truth",
129  "event:seed:vx:vy:vz:ntracks:chi2:ndof:"
130  "gvx:gvy:gvz:gvt:gembed:gntracks:gntracksmaps:"
131  "gnembed:nfromtruth:"
132  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
133 
134 
135 
136  if (_do_gpoint_eval) _ntp_gpoint = new TNtuple("ntp_gpoint", "g4point => best vertex",
137  "event:seed:gvx:gvy:gvz:gvt:gntracks:gembed:"
138  "vx:vy:vz:ntracks:"
139  "nfromtruth:"
140  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
141 
142  if (_do_g4hit_eval) _ntp_g4hit = new TNtuple("ntp_g4hit", "g4hit => best svtxcluster",
143  "event:seed:g4hitID:gx:gy:gz:gt:gpl:gedep:geta:gphi:"
144  "gdphi:gdz:"
145  "glayer:gtrackID:gflavor:"
146  "gpx:gpy:gpz:"
147  "gvx:gvy:gvz:"
148  "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
149  "gembed:gprimary:nclusters:"
150  "clusID:x:y:z:eta:phi:e:adc:layer:size:"
151  "efromtruth:dphitru:detatru:dztru:drtru:"
152  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
153 
154  if (_do_hit_eval) _ntp_hit = new TNtuple("ntp_hit", "svtxhit => max truth",
155  "event:seed:hitID:e:adc:layer:phielem:zelem:"
156  "cellID:ecell:phibin:zbin:phi:z:"
157  "g4hitID:gedep:gx:gy:gz:gt:"
158  "gtrackID:gflavor:"
159  "gpx:gpy:gpz:gvx:gvy:gvz:gvt:"
160  "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
161  "gembed:gprimary:efromtruth:"
162  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
163 
164  if (_do_cluster_eval) _ntp_cluster = new TNtuple("ntp_cluster", "svtxcluster => max truth",
165  "event:seed:hitID:x:y:z:r:phi:eta:theta:ex:ey:ez:ephi:"
166  "e:adc:maxadc:layer:phielem:zelem:size:"
167  "trackID:niter:g4hitID:gx:"
168  "gy:gz:gr:gphi:geta:gt:gtrackID:gflavor:"
169  "gpx:gpy:gpz:gvx:gvy:gvz:gvt:"
170  "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
171  "gembed:gprimary:efromtruth:nparticles:"
172  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
173 
174  if (_do_g4cluster_eval) _ntp_g4cluster = new TNtuple("ntp_g4cluster", "g4cluster => max truth",
175  "event:layer:gx:gy:gz:gt:gedep:gr:gphi:geta:gtrackID:gflavor:gembed:gprimary:gphisize:gzsize:gadc:nreco:x:y:z:r:phi:eta:ex:ey:ez:ephi:adc");
176 
177  if (_do_gtrack_eval) _ntp_gtrack = new TNtuple("ntp_gtrack", "g4particle => best svtxtrack",
178  "event:seed:gntracks:gtrackID:gflavor:gnhits:gnmaps:gnintt:gnmms:"
179  "gnintt1:gnintt2:gnintt3:gnintt4:"
180  "gnintt5:gnintt6:gnintt7:gnintt8:"
181  "gntpc:gnlmaps:gnlintt:gnltpc:gnlmms:"
182  "gpx:gpy:gpz:gpt:geta:gphi:"
183  "gvx:gvy:gvz:gvt:"
184  "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
185  "gembed:gprimary:"
186  "trackID:px:py:pz:pt:eta:phi:deltapt:deltaeta:deltaphi:"
187  "charge:quality:chisq:ndf:nhits:layers:nmaps:nintt:ntpc:nmms:ntpc1:ntpc11:ntpc2:ntpc3:nlmaps:nlintt:nltpc:nlmms:"
188  "dca2d:dca2dsigma:dca3dxy:dca3dxysigma:dca3dz:dca3dzsigma:pcax:pcay:pcaz:nfromtruth:nwrong:ntrumaps:ntruintt:ntrutpc:ntrumms:ntrutpc1:ntrutpc11:ntrutpc2:ntrutpc3:layersfromtruth:"
189  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
190 
191  if (_do_track_eval) _ntp_track = new TNtuple("ntp_track", "svtxtrack => max truth",
192  "event:seed:trackID:px:py:pz:pt:eta:phi:deltapt:deltaeta:deltaphi:charge:"
193  "quality:chisq:ndf:nhits:nmaps:nintt:ntpc:nmms:ntpc1:ntpc11:ntpc2:ntpc3:nlmaps:nlintt:nltpc:nlmms:layers:"
194  "dca2d:dca2dsigma:dca3dxy:dca3dxysigma:dca3dz:dca3dzsigma:pcax:pcay:pcaz:"
195  "presdphi:presdeta:prese3x3:prese:"
196  "cemcdphi:cemcdeta:cemce3x3:cemce:"
197  "hcalindphi:hcalindeta:hcaline3x3:hcaline:"
198  "hcaloutdphi:hcaloutdeta:hcaloute3x3:hcaloute:"
199  "gtrackID:gflavor:gnhits:gnmaps:gnintt:gntpc:gnmms:gnlmaps:gnlintt:gnltpc:gnlmms:"
200  "gpx:gpy:gpz:gpt:geta:gphi:"
201  "gvx:gvy:gvz:gvt:"
202  "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
203  "gembed:gprimary:nfromtruth:nwrong:ntrumaps:ntruintt:"
204  "ntrutpc:ntrumms:ntrutpc1:ntrutpc11:ntrutpc2:ntrutpc3:layersfromtruth:"
205  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
206 
207  if (_do_gseed_eval) _ntp_gseed = new TNtuple("ntp_gseed", "seeds from truth",
208  "event:seed:ntrk:gx:gy:gz:gr:geta:gphi:"
209  "glayer:"
210  "gpx:gpy:gpz:gtpt:gtphi:gteta:"
211  "gvx:gvy:gvz:"
212  "gembed:gprimary:gflav:"
213  "dphiprev:detaprev:"
214  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
215 
216  _timer = new PHTimer("_eval_timer");
217  _timer->stop();
218 
219  return Fun4AllReturnCodes::EVENT_OK;
220 }
221 
222 int SvtxEvaluator::InitRun(PHCompositeNode* /*topNode*/)
223 {
224  return Fun4AllReturnCodes::EVENT_OK;
225 }
226 
227 int SvtxEvaluator::process_event(PHCompositeNode* topNode)
228 {
229  if ((Verbosity() > 1) && (_ievent % 100 == 0))
230  {
231  cout << "SvtxEvaluator::process_event - Event = " << _ievent << endl;
232  }
233 
234  recoConsts *rc = recoConsts::instance();
235  if (rc->FlagExist("RANDOMSEED"))
236  {
237  _iseed = rc->get_IntFlag("RANDOMSEED");
238  m_fSeed = _iseed;
239  }
240  else
241  {
242  _iseed = 0;
243  m_fSeed = NAN;
244  }
245 
246  if(Verbosity() > 1)
247  {
248  cout << "SvtxEvaluator::process_event - Seed = " << _iseed << endl;
249  }
250 
251  if (!_svtxevalstack)
252  {
253  _svtxevalstack = new SvtxEvalStack(topNode);
254  _svtxevalstack->set_strict(_strict);
255  _svtxevalstack->set_verbosity(Verbosity());
256  _svtxevalstack->set_use_initial_vertex(_use_initial_vertex);
257  _svtxevalstack->set_use_genfit_vertex(_use_genfit_vertex);
258  _svtxevalstack->next_event(topNode);
259  }
260  else
261  {
262  _svtxevalstack->next_event(topNode);
263  }
264 
265  //-----------------------------------
266  // print what is coming into the code
267  //-----------------------------------
268 
269  printInputInfo(topNode);
270 
271  //---------------------------
272  // fill the Evaluator NTuples
273  //---------------------------
274 
275  fillOutputNtuples(topNode);
276 
277  //--------------------------------------------------
278  // Print out the ancestry information for this event
279  //--------------------------------------------------
280 
281  //printOutputInfo(topNode);
282 
283  ++_ievent;
284  return Fun4AllReturnCodes::EVENT_OK;
285 }
286 
287 int SvtxEvaluator::End(PHCompositeNode* /*topNode*/)
288 {
289  _tfile->cd();
290 
291  if (_ntp_info) _ntp_info->Write();
292  if (_ntp_vertex) _ntp_vertex->Write();
293  if (_ntp_gpoint) _ntp_gpoint->Write();
294  if (_ntp_g4hit) _ntp_g4hit->Write();
295  if (_ntp_hit) _ntp_hit->Write();
296  if (_ntp_cluster) _ntp_cluster->Write();
297  if (_ntp_g4cluster) _ntp_g4cluster->Write();
298  if (_ntp_gtrack) _ntp_gtrack->Write();
299  if (_ntp_track) _ntp_track->Write();
300  if (_ntp_gseed) _ntp_gseed->Write();
301 
302  _tfile->Close();
303 
304  delete _tfile;
305 
306  if (Verbosity() > 1)
307  {
308  cout << "========================= SvtxEvaluator::End() ============================" << endl;
309  cout << " " << _ievent << " events of output written to: " << _filename << endl;
310  cout << "===========================================================================" << endl;
311  }
312 
313  if(_svtxevalstack)
314  {
315  _errors += _svtxevalstack->get_errors();
316 
317  if (Verbosity() > 1)
318  {
319  if ((_errors > 0) || (Verbosity() > 1))
320  {
321  cout << "SvtxEvaluator::End() - Error Count: " << _errors << endl;
322  }
323  }
324 
325  delete _svtxevalstack;
326  }
327 
328  return Fun4AllReturnCodes::EVENT_OK;
329 }
330 
331 void SvtxEvaluator::printInputInfo(PHCompositeNode* topNode)
332 {
333  if (Verbosity() > 1) cout << "SvtxEvaluator::printInputInfo() entered" << endl;
334 
335  if (Verbosity() > 3)
336  {
337  // event information
338  cout << endl;
339  cout << PHWHERE << " INPUT FOR EVENT " << _ievent << endl;
340 
341  cout << endl;
342  cout << "---PHG4HITS-------------" << endl;
343  _svtxevalstack->get_truth_eval()->set_strict(_strict);
344  std::set<PHG4Hit*> g4hits = _svtxevalstack->get_truth_eval()->all_truth_hits();
345  unsigned int ig4hit = 0;
346  for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
347  iter != g4hits.end();
348  ++iter)
349  {
350  PHG4Hit* g4hit = *iter;
351  cout << ig4hit << " of " << g4hits.size();
352  cout << ": PHG4Hit: " << endl;
353  g4hit->identify();
354  ++ig4hit;
355  }
356 
357  cout << "---SVTXCLUSTERS-------------" << endl;
358  TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
359  if(!clustermap)
360  clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
361 
362  TrkrHitSetContainer* hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
363 
364  if (clustermap!=nullptr&&hitsets!=nullptr){
365  unsigned int icluster = 0;
366  auto hitsetrange = hitsets->getHitSets(TrkrDefs::TrkrId::inttId);
367  for (auto hitsetitr = hitsetrange.first;
368  hitsetitr != hitsetrange.second;
369  ++hitsetitr){
370  auto range = clustermap->getClusters(hitsetitr->first);
371  for( auto iter = range.first; iter != range.second; ++iter ){
372  TrkrDefs::cluskey cluster_key = iter->first;
373  cout << icluster << " with key " << cluster_key << " of " << clustermap->size();
374  cout << ": SvtxCluster: " << endl;
375  iter->second->identify();
376  ++icluster;
377  }
378  }
379  }
380 
381  cout << "---SVXTRACKS-------------" << endl;
382  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname.c_str());
383  if (trackmap)
384  {
385  unsigned int itrack = 0;
386  for (SvtxTrackMap::Iter iter = trackmap->begin();
387  iter != trackmap->end();
388  ++iter)
389  {
390  cout << itrack << " of " << trackmap->size();
391  SvtxTrack* track = iter->second;
392  cout << " : SvtxTrack:" << endl;
393  track->identify();
394  cout << endl;
395  ++itrack;
396  }
397  }
398 
399  cout << "---SVXVERTEXES-------------" << endl;
400  SvtxVertexMap* vertexmap = nullptr;
401  if(_use_initial_vertex)
402  vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
403  else if (_use_genfit_vertex)
404  vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapRefit");
405  else
406  vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs"); // Acts vertices
407 
408  if (vertexmap)
409  {
410  unsigned int ivertex = 0;
411  for (SvtxVertexMap::Iter iter = vertexmap->begin();
412  iter != vertexmap->end();
413  ++iter)
414  {
415  cout << ivertex << " of " << vertexmap->size();
416  SvtxVertex* vertex = iter->second;
417  cout << " : SvtxVertex:" << endl;
418  vertex->identify();
419  cout << endl;
420  }
421  }
422  }
423 
424  return;
425 }
426 
427 void SvtxEvaluator::printOutputInfo(PHCompositeNode* topNode)
428 {
429  if (Verbosity() > 1) cout << "SvtxEvaluator::printOutputInfo() entered" << endl;
430 
431  //==========================================
432  // print out some useful stuff for debugging
433  //==========================================
434 
435  if (Verbosity() > 100)
436  {
437  SvtxTrackEval* trackeval = _svtxevalstack->get_track_eval();
438  SvtxClusterEval* clustereval = _svtxevalstack->get_cluster_eval();
439  SvtxTruthEval* trutheval = _svtxevalstack->get_truth_eval();
440 
441  // event information
442  cout << endl;
443  cout << PHWHERE << " NEW OUTPUT FOR EVENT " << _ievent << endl;
444  cout << endl;
445 
446  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
447 
448  PHG4VtxPoint* gvertex = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
449  float gvx = gvertex->get_x();
450  float gvy = gvertex->get_y();
451  float gvz = gvertex->get_z();
452 
453  float vx = NAN;
454  float vy = NAN;
455  float vz = NAN;
456 
457  SvtxVertexMap* vertexmap = nullptr;
458  if(_use_initial_vertex)
459  vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
460  else if (_use_genfit_vertex)
461  vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapRefit");
462  else
463  vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs"); // Acts vertices
464 
465  if (vertexmap)
466  {
467  if (!vertexmap->empty())
468  {
469  SvtxVertex* vertex = (vertexmap->begin()->second);
470 
471  vx = vertex->get_x();
472  vy = vertex->get_y();
473  vz = vertex->get_z();
474  }
475  }
476 
477  cout << "===Vertex Reconstruction=======================" << endl;
478  cout << "vtrue = (" << gvx << "," << gvy << "," << gvz << ") => vreco = (" << vx << "," << vy << "," << vz << ")" << endl;
479  cout << endl;
480 
481  cout << "===Tracking Summary============================" << endl;
482  unsigned int ng4hits[100] = {0};
483  std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits();
484  for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
485  iter != g4hits.end();
486  ++iter)
487  {
488  PHG4Hit* g4hit = *iter;
489  ++ng4hits[g4hit->get_layer()];
490  }
491 
492  TrkrHitSetContainer *hitsetmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
493 
494  TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
495  if(!clustermap)
496  clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
497 
498  unsigned int nclusters[100] = {0};
499  unsigned int nhits[100] = {0};
500 
501  ActsTrackingGeometry *tgeometry = findNode::getClass<ActsTrackingGeometry>(topNode,"ActsTrackingGeometry");
502  ActsSurfaceMaps *surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,"ActsSurfaceMaps");
503  if(!tgeometry or !surfmaps)
504  {
505  std::cout << PHWHERE << "No Acts geometry on node tree. Can't continue."
506  << std::endl;
507  }
508 
509  if (hitsetmap)
510  {
511  TrkrHitSetContainer::ConstRange all_hitsets = hitsetmap->getHitSets();
512  for (TrkrHitSetContainer::ConstIterator hitsetiter = all_hitsets.first;
513  hitsetiter != all_hitsets.second;
514  ++hitsetiter){
515  // we have a single hitset, get the layer
516  unsigned int layer = TrkrDefs::getLayer(hitsetiter->first);
517 
518  // count all hits in this hitset
519  TrkrHitSet::ConstRange hitrangei = hitsetiter->second->getHits();
520  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
521  hitr != hitrangei.second;
522  ++hitr)
523  {
524  ++nhits[layer];
525  }
526  auto range = clustermap->getClusters(hitsetiter->first);
527  for( auto clusIter = range.first; clusIter != range.second; ++clusIter ){
528  const auto cluskey = clusIter->first;
529  // const auto cluster = clusIter->second;
530  unsigned int layer = TrkrDefs::getLayer(cluskey);
531  nclusters[layer]++;
532  }
533  }
534  }
535 
536  PHG4CylinderCellGeomContainer* geom_container =
537  findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
538  {
539  if (!geom_container)
540  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
541  return;
542  }
543 
544 
545  for (unsigned int ilayer = 0; ilayer < _nlayers_maps + _nlayers_intt + _nlayers_tpc; ++ilayer)
546  {
547  PHG4CylinderCellGeom* GeoLayer = geom_container->GetLayerCellGeom(ilayer);
548 
549  cout << "layer " << ilayer << ": nG4hits = " << ng4hits[ilayer]
550  << " => nHits = " << nhits[ilayer]
551  << " => nClusters = " << nclusters[ilayer] << endl;
552  if(ilayer>=_nlayers_maps + _nlayers_intt && ilayer < _nlayers_maps + _nlayers_intt + _nlayers_tpc){
553  cout << "layer " << ilayer
554  << " => nphi = " << GeoLayer->get_phibins()
555  << " => nz = " << GeoLayer->get_zbins()
556  << " => ntot = " << GeoLayer->get_phibins() * GeoLayer->get_zbins()
557  << endl;
558  }
559  }
560 
561  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname.c_str());
562 
563  cout << "nGtracks = " << std::distance(truthinfo->GetPrimaryParticleRange().first, truthinfo->GetPrimaryParticleRange().second);
564  cout << " => nTracks = ";
565  if (trackmap)
566  cout << trackmap->size() << endl;
567  else
568  cout << 0 << endl;
569 
570  // cluster wise information
571  if (Verbosity() > 1)
572  {
573  for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
574  iter != g4hits.end();
575  ++iter)
576  {
577  PHG4Hit* g4hit = *iter;
578 
579  cout << endl;
580  cout << "===PHG4Hit===================================" << endl;
581  cout << " PHG4Hit: ";
582  g4hit->identify();
583 
584  std::set<TrkrDefs::cluskey> clusters = clustereval->all_clusters_from(g4hit);
585 
586  for (std::set<TrkrDefs::cluskey>::iterator jter = clusters.begin();
587  jter != clusters.end();
588  ++jter)
589  {
590  TrkrDefs::cluskey cluster_key = *jter;
591  cout << "===Created-SvtxCluster================" << endl;
592  cout << "SvtxCluster: ";
593  TrkrCluster *cluster = clustermap->findCluster(cluster_key);
594  cluster->identify();
595  }
596  }
597 
598  PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
599  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
600  iter != range.second;
601  ++iter)
602  {
603  PHG4Particle* particle = iter->second;
604 
605  // track-wise information
606  cout << endl;
607 
608  cout << "=== Gtrack ===================================================" << endl;
609  cout << " PHG4Particle id = " << particle->get_track_id() << endl;
610  particle->identify();
611  cout << " ptrue = (";
612  cout.width(5);
613  cout << particle->get_px();
614  cout << ",";
615  cout.width(5);
616  cout << particle->get_py();
617  cout << ",";
618  cout.width(5);
619  cout << particle->get_pz();
620  cout << ")" << endl;
621 
622  cout << " vtrue = (";
623  cout.width(5);
624  cout << truthinfo->GetVtx(particle->get_vtx_id())->get_x();
625  cout << ",";
626  cout.width(5);
627  cout << truthinfo->GetVtx(particle->get_vtx_id())->get_y();
628  cout << ",";
629  cout.width(5);
630  cout << truthinfo->GetVtx(particle->get_vtx_id())->get_z();
631  cout << ")" << endl;
632 
633  cout << " pt = " << sqrt(pow(particle->get_px(), 2) + pow(particle->get_py(), 2)) << endl;
634  cout << " phi = " << atan2(particle->get_py(), particle->get_px()) << endl;
635  cout << " eta = " << asinh(particle->get_pz() / sqrt(pow(particle->get_px(), 2) + pow(particle->get_py(), 2))) << endl;
636 
637  cout << " embed flag = " << truthinfo->isEmbeded(particle->get_track_id()) << endl;
638 
639  cout << " ---Associated-PHG4Hits-----------------------------------------" << endl;
640  std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(particle);
641  for (std::set<PHG4Hit*>::iterator jter = g4hits.begin();
642  jter != g4hits.end();
643  ++jter)
644  {
645  PHG4Hit* g4hit = *jter;
646 
647  float x = 0.5 * (g4hit->get_x(0) + g4hit->get_x(1));
648  float y = 0.5 * (g4hit->get_y(0) + g4hit->get_y(1));
649  float z = 0.5 * (g4hit->get_z(0) + g4hit->get_z(1));
650 
651  cout << " #" << g4hit->get_hit_id() << " xtrue = (";
652  cout.width(5);
653  cout << x;
654  cout << ",";
655  cout.width(5);
656  cout << y;
657  cout << ",";
658  cout.width(5);
659  cout << z;
660  cout << ")";
661 
662  std::set<TrkrDefs::cluskey> clusters = clustereval->all_clusters_from(g4hit);
663  for (std::set<TrkrDefs::cluskey>::iterator kter = clusters.begin();
664  kter != clusters.end();
665  ++kter)
666  {
667  TrkrDefs::cluskey cluster_key = *kter;
668 
669  TrkrCluster *cluster = clustermap->findCluster(cluster_key);
670  ActsTransformations transformer;
671  auto glob = transformer.getGlobalPosition(cluster,surfmaps,tgeometry);
672  float x = glob(0);
673  float y = glob(1);
674  float z = glob(2);
675 
676  cout << " => #" << cluster_key;
677  cout << " xreco = (";
678  cout.width(5);
679  cout << x;
680  cout << ",";
681  cout.width(5);
682  cout << y;
683  cout << ",";
684  cout.width(5);
685  cout << z;
686  cout << ")";
687  }
688 
689  cout << endl;
690  }
691 
692  if (trackmap && clustermap)
693  {
694  std::set<SvtxTrack*> tracks = trackeval->all_tracks_from(particle);
695  for (std::set<SvtxTrack*>::iterator jter = tracks.begin();
696  jter != tracks.end();
697  ++jter)
698  {
699  SvtxTrack* track = *jter;
700 
701  float px = track->get_px();
702  float py = track->get_py();
703  float pz = track->get_pz();
704 
705  cout << "===Created-SvtxTrack==========================================" << endl;
706  cout << " SvtxTrack id = " << track->get_id() << endl;
707  cout << " preco = (";
708  cout.width(5);
709  cout << px;
710  cout << ",";
711  cout.width(5);
712  cout << py;
713  cout << ",";
714  cout.width(5);
715  cout << pz;
716  cout << ")" << endl;
717  cout << " quality = " << track->get_quality() << endl;
718  cout << " nfromtruth = " << trackeval->get_nclusters_contribution(track, particle) << endl;
719 
720  cout << " ---Associated-SvtxClusters-to-PHG4Hits-------------------------" << endl;
721 
722  for (SvtxTrack::ConstClusterKeyIter iter = track->begin_cluster_keys();
723  iter != track->end_cluster_keys();
724  ++iter)
725  {
726  TrkrDefs::cluskey cluster_key = *iter;
727  TrkrCluster *cluster = clustermap->findCluster(cluster_key);
728  ActsTransformations transformer;
729  auto glob = transformer.getGlobalPosition(cluster,surfmaps,tgeometry);
730  float x = glob(0);
731  float y = glob(1);
732  float z = glob(2);
733 
734  cout << " #" << cluster_key << " xreco = (";
735  cout.width(5);
736  cout << x;
737  cout << ",";
738  cout.width(5);
739  cout << y;
740  cout << ",";
741  cout.width(5);
742  cout << z;
743  cout << ") =>";
744 
745  PHG4Hit* g4hit = clustereval->max_truth_hit_by_energy(cluster_key);
746  if ((g4hit) && (g4hit->get_trkid() == particle->get_track_id()))
747  {
748  x = 0.5 * (g4hit->get_x(0) + g4hit->get_x(1));
749  y = 0.5 * (g4hit->get_y(0) + g4hit->get_y(1));
750  z = 0.5 * (g4hit->get_z(0) + g4hit->get_z(1));
751 
752  cout << " #" << g4hit->get_hit_id()
753  << " xtrue = (";
754  cout.width(5);
755  cout << x;
756  cout << ",";
757  cout.width(5);
758  cout << y;
759  cout << ",";
760  cout.width(5);
761  cout << z;
762  cout << ") => Gtrack id = " << g4hit->get_trkid();
763  }
764  else
765  {
766  cout << " noise hit";
767  }
768  }
769 
770  cout << endl;
771  }
772  }
773  }
774  }
775 
776  cout << endl;
777 
778  } // if Verbosity()
779 
780  return;
781 }
782 
783 void SvtxEvaluator::fillOutputNtuples(PHCompositeNode* topNode)
784 {
785  if (Verbosity() > 1) cout << "SvtxEvaluator::fillOutputNtuples() entered" << endl;
786 
787  ActsTrackingGeometry *tgeometry = findNode::getClass<ActsTrackingGeometry>(topNode,"ActsTrackingGeometry");
788  ActsSurfaceMaps *surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,"ActsSurfaceMaps");
789  if(!tgeometry or !surfmaps)
790  {
791  std::cout << "No Acts geometry on node tree. Can't continue."
792  << std::endl;
793  return;
794  }
795 
796  ActsTransformations actsTransformer;
797 
798  SvtxVertexEval* vertexeval = _svtxevalstack->get_vertex_eval();
799  SvtxClusterEval* clustereval = _svtxevalstack->get_cluster_eval();
800  SvtxTrackEval* trackeval = _svtxevalstack->get_track_eval();
801  SvtxHitEval* hiteval = _svtxevalstack->get_hit_eval();
802  SvtxTruthEval* trutheval = _svtxevalstack->get_truth_eval();
803 
804  float nhit_tpc_all = 0;
805  float nhit_tpc_in = 0;
806  float nhit_tpc_mid = 0;
807  float nhit_tpc_out = 0;
808  float nclus_all = 0;
809  float nclus_tpc = 0;
810  float nclus_intt = 0;
811  float nclus_maps = 0;
812  float nclus_mms = 0;
813  float nhit[100];
814  for(int i = 0; i<100;i++)nhit[i] = 0;
815  float occ11 = 0;
816  float occ116 = 0;
817  float occ21 = 0;
818  float occ216 = 0;
819  float occ31 = 0;
820  float occ316 = 0;
821 
822  TrkrHitSetContainer* hitmap_in = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
823  if (hitmap_in)
824  {
825  TrkrHitSetContainer::ConstRange all_hitsets = hitmap_in->getHitSets();
826  for (TrkrHitSetContainer::ConstIterator hitsetiter = all_hitsets.first;
827  hitsetiter != all_hitsets.second;
828  ++hitsetiter)
829  {
830  // we have a single hitset, get the layer
831  unsigned int layer = TrkrDefs::getLayer(hitsetiter->first);
832  if(layer >= _nlayers_maps + _nlayers_intt && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
833  {
834  // count all hits in this hitset
835  TrkrHitSet::ConstRange hitrangei = hitsetiter->second->getHits();
836  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
837  hitr != hitrangei.second;
838  ++hitr)
839  {
840  nhit[layer]++;
841  nhit_tpc_all++;
842  if ((float) layer == _nlayers_maps + _nlayers_intt) nhit_tpc_in++;
843  if ((float) layer == _nlayers_maps + _nlayers_intt + _nlayers_tpc - 1) nhit_tpc_out++;
844  if ((float) layer == _nlayers_maps + _nlayers_intt + _nlayers_tpc / 2 - 1) nhit_tpc_mid++;
845  }
846  }
847  }
848  }
849  /**********/
850  PHG4CylinderCellGeomContainer* geom_container =
851  findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
852  if (!geom_container)
853  {
854  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
855  return;
856  }
857  PHG4CylinderCellGeom* GeoLayer;
858  int layer = _nlayers_maps + _nlayers_intt;
859  GeoLayer = geom_container->GetLayerCellGeom(layer);
860  int nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
861  float nhits = nhit[layer];
862  occ11 = nhits/nbins;
863  if(Verbosity() > 1)
864  {
865  cout << " occ11 = " << occ11
866  << " nbins = " << nbins
867  << " nhits = " << nhits
868  << " layer = " << layer
869  << endl;
870  }
871  layer = _nlayers_maps + _nlayers_intt + 15;
872  GeoLayer = geom_container->GetLayerCellGeom(layer);
873  nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
874  nhits = nhit[layer];
875  occ116 = nhits/nbins;
876 
877  layer = _nlayers_maps + _nlayers_intt + 16;
878  GeoLayer = geom_container->GetLayerCellGeom(layer);
879  nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
880  nhits = nhit[layer];
881  occ21 = nhits/nbins;
882  layer = _nlayers_maps + _nlayers_intt + 31;
883  GeoLayer = geom_container->GetLayerCellGeom(layer);
884  nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
885  nhits = nhit[layer];
886  occ216 = nhits/nbins;
887  layer = _nlayers_maps + _nlayers_intt + 32;
888  GeoLayer = geom_container->GetLayerCellGeom(layer);
889  nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
890  nhits = nhit[layer];
891  occ31 = nhits/nbins;
892  layer = _nlayers_maps + _nlayers_intt + 47;
893  GeoLayer = geom_container->GetLayerCellGeom(layer);
894  nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
895  nhits = nhit[layer];
896  occ316 = nhits/nbins;
897 
898  if(Verbosity() > 1)
899  {
900  cout << " occ11 = " << occ11
901  << " occ116 = " << occ116
902  << " occ21 = " << occ21
903  << " occ216 = " << occ216
904  << " occ31 = " << occ31
905  << " occ316 = " << occ316
906  << endl;
907  }
908  TrkrHitSetContainer* hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
909 
910  TrkrClusterContainer* clustermap_in = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
911  if(!clustermap_in)
912  clustermap_in = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
913 
914  nclus_all = clustermap_in->size();
915 
916  auto hitsetrange = hitsets->getHitSets(TrkrDefs::TrkrId::inttId);
917  for (auto hitsetitr = hitsetrange.first;
918  hitsetitr != hitsetrange.second;
919  ++hitsetitr){
920  auto range = clustermap_in->getClusters(hitsetitr->first);
921  for( auto iter_cin = range.first; iter_cin != range.second; ++iter_cin ){
922  TrkrDefs::cluskey cluster_key = iter_cin->first;
923  unsigned int layer = TrkrDefs::getLayer(cluster_key);
924  if (_nlayers_maps > 0)
925  if (layer < _nlayers_maps) nclus_maps++;
926  if (_nlayers_intt > 0)
927  if (layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt) nclus_intt++;
928  if (_nlayers_tpc > 0)
929  if (layer >= _nlayers_maps + _nlayers_intt && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc) nclus_tpc++;
930  if (_nlayers_mms > 0)
931  if (layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc) nclus_mms++;
932  }
933  }
934  //-----------------------
935  // fill the info NTuple
936  //-----------------------
937  if (_ntp_info)
938  {
939  if (Verbosity() > 1)
940  {
941  cout << "Filling ntp_info " << endl;
942  }
943  float ntrk = 0;
944  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname.c_str());
945  if (trackmap)
946  ntrk = (float) trackmap->size();
947  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
948  int nprim = truthinfo->GetNumPrimaryVertexParticles();
949  if (Verbosity() > 0){
950  cout << "EVENTINFO SEED: " << m_fSeed << endl;
951  cout << "EVENTINFO NHIT: " << setprecision(9) << nhit_tpc_all << endl;
952  cout << "EVENTINFO NTRKGEN: " << nprim << endl;
953  cout << "EVENTINFO NTRKREC: " << ntrk << endl;
954 
955  }
956  float info_data[] = {(float) _ievent,m_fSeed,
957  occ11,occ116,occ21,occ216,occ31,occ316,
958  (float) nprim,
959  0,
960  ntrk,
961  nhit_tpc_all,
962  nhit_tpc_in,
963  nhit_tpc_mid,
964  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
965 
966  _ntp_info->Fill(info_data);
967  }
968 
969 
970  //-----------------------
971  // fill the Vertex NTuple
972  //-----------------------
973  bool doit = true;
974  if (_ntp_vertex && doit)
975  {
976  if (Verbosity() > 1)
977  {
978  cout << "Filling ntp_vertex " << endl;
979  cout << "start vertex time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
980  _timer->restart();
981  }
982 
983  SvtxVertexMap* vertexmap = nullptr;
984  if(_use_initial_vertex)
985  vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
986  else if (_use_genfit_vertex)
987  vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapRefit");
988  else
989  vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs"); // Acts vertices
990 
991  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
992 
993  if (vertexmap && truthinfo)
994  {
995  const auto prange = truthinfo->GetPrimaryParticleRange();
996  map<int, unsigned int> embedvtxid_particle_count;
997  map<int, unsigned int> embedvtxid_maps_particle_count;
998  map<int, unsigned int> vertex_particle_count;
999 
1000  if (_do_vtx_eval_light == false)
1001  {
1002  for (auto iter = prange.first; iter != prange.second; ++iter) // process all primary paricle
1003  {
1004  const int point_id = iter->second->get_vtx_id();
1005  int gembed = truthinfo->isEmbededVtx(iter->second->get_vtx_id());
1006  ++vertex_particle_count[point_id];
1007  ++embedvtxid_particle_count[gembed];
1008  PHG4Particle* g4particle = iter->second;
1009 
1010  if (_scan_for_embedded && gembed <= 0) continue;
1011 
1012  std::set<TrkrDefs::cluskey> g4clusters = clustereval->all_clusters_from(g4particle);
1013  unsigned int nglmaps = 0;
1014 
1015  int lmaps[_nlayers_maps + 1];
1016  if (_nlayers_maps > 0)
1017  {
1018  for (unsigned int i = 0; i < _nlayers_maps; i++)
1019  {
1020  lmaps[i] = 0;
1021  }
1022  }
1023  for (const TrkrDefs::cluskey g4cluster : g4clusters)
1024  {
1025  unsigned int layer = TrkrDefs::getLayer(g4cluster);
1026  //cout<<__LINE__<<": " << _ievent <<": " <<gtrackID << ": " << layer <<": " <<g4cluster->get_id() <<endl;
1027  if (_nlayers_maps > 0 && layer < _nlayers_maps)
1028  {
1029  lmaps[layer] = 1;
1030  }
1031  }
1032  if (_nlayers_maps > 0)
1033  {
1034  for (unsigned int i = 0; i < _nlayers_maps; i++)
1035  {
1036  nglmaps += lmaps[i];
1037  }
1038  }
1039  float gpx = g4particle->get_px();
1040  float gpy = g4particle->get_py();
1041  float gpz = g4particle->get_pz();
1042  float gpt = NAN;
1043  float geta = NAN;
1044 
1045  if (gpx != 0 && gpy != 0)
1046  {
1047  TVector3 gv(gpx, gpy, gpz);
1048  gpt = gv.Pt();
1049  geta = gv.Eta();
1050  // gphi = gv.Phi();
1051  }
1052 
1053  if (nglmaps == 3 && fabs(geta) < 1.0 && gpt > 0.5)
1054  ++embedvtxid_maps_particle_count[gembed];
1055  }
1056  }
1057 
1058  auto vrange = truthinfo->GetPrimaryVtxRange();
1059  map<int, bool> embedvtxid_found;
1060  map<int, int> embedvtxid_vertex_id;
1061  map<int, PHG4VtxPoint*> embedvtxid_vertex;
1062  for (auto iter = vrange.first; iter != vrange.second; ++iter) // process all primary vertexes
1063  {
1064  const int point_id = iter->first;
1065  int gembed = truthinfo->isEmbededVtx(point_id);
1066 
1067  if (_scan_for_embedded && gembed <= 0) continue;
1068 
1069  auto search = embedvtxid_found.find(gembed);
1070  if (search != embedvtxid_found.end())
1071  {
1072  embedvtxid_vertex_id[gembed] = point_id;
1073  embedvtxid_vertex[gembed] = iter->second;
1074  }
1075  else
1076  {
1077  if (vertex_particle_count[embedvtxid_vertex_id[gembed]] < vertex_particle_count[point_id])
1078  {
1079  embedvtxid_vertex_id[gembed] = point_id;
1080  embedvtxid_vertex[gembed] = iter->second;
1081  }
1082  }
1083  embedvtxid_found[gembed] = false;
1084  }
1085 
1086  unsigned int ngembed = 0;
1087  for (std::map<int, bool>::iterator iter = embedvtxid_found.begin();
1088  iter != embedvtxid_found.end();
1089  ++iter)
1090  {
1091  if (iter->first >= 0 || iter->first != iter->first) continue;
1092  ++ngembed;
1093  }
1094 
1095  for (SvtxVertexMap::Iter iter = vertexmap->begin();
1096  iter != vertexmap->end();
1097  ++iter)
1098  {
1099  SvtxVertex* vertex = iter->second;
1100  PHG4VtxPoint* point = vertexeval->max_truth_point_by_ntracks(vertex);
1101  float vx = vertex->get_x();
1102  float vy = vertex->get_y();
1103  float vz = vertex->get_z();
1104  float ntracks = vertex->size_tracks();
1105  float chi2 = vertex->get_chisq();
1106  float ndof = vertex->get_ndof();
1107  float gvx = NAN;
1108  float gvy = NAN;
1109  float gvz = NAN;
1110  float gvt = NAN;
1111  float gembed = NAN;
1112  float gntracks = truthinfo->GetNumPrimaryVertexParticles();
1113  float gntracksmaps = NAN;
1114  float gnembed = NAN;
1115  float nfromtruth = NAN;
1116  if (point)
1117  {
1118  const int point_id = point->get_id();
1119  gvx = point->get_x();
1120  gvy = point->get_y();
1121  gvz = point->get_z();
1122  gvt = point->get_t();
1123  gembed = truthinfo->isEmbededVtx(point_id);
1124  gntracks = embedvtxid_particle_count[(int) gembed];
1125  if (embedvtxid_maps_particle_count[(int) gembed] > 0 && fabs(gvt) < 2000. && fabs(gvz) < 13.0)
1126  gntracksmaps = embedvtxid_maps_particle_count[(int) gembed];
1127  gnembed = (float) ngembed;
1128  nfromtruth = vertexeval->get_ntracks_contribution(vertex, point);
1129  embedvtxid_found[(int) gembed] = true;
1130  }
1131 
1132  float vertex_data[] = {(float) _ievent,m_fSeed,
1133  vx,
1134  vy,
1135  vz,
1136  ntracks,
1137  chi2,
1138  ndof,
1139  gvx,
1140  gvy,
1141  gvz,
1142  gvt,
1143  gembed,
1144  gntracks,
1145  gntracksmaps,
1146  gnembed,
1147  nfromtruth,
1148  nhit_tpc_all,
1149  nhit_tpc_in,
1150  nhit_tpc_mid,
1151  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps,nclus_mms};
1152 
1153  _ntp_vertex->Fill(vertex_data);
1154  }
1155 
1156  if (!_scan_for_embedded)
1157  {
1158  for (std::map<int, bool>::iterator iter = embedvtxid_found.begin();
1159  iter != embedvtxid_found.end();
1160  ++iter)
1161  {
1162  if (embedvtxid_found[iter->first]) continue;
1163 
1164  float vx = NAN;
1165  float vy = NAN;
1166  float vz = NAN;
1167  float ntracks = NAN;
1168 
1169  float gvx = NAN;
1170  float gvy = NAN;
1171  float gvz = NAN;
1172  float gvt = NAN;
1173  float gembed = iter->first;
1174  float gntracks = NAN;
1175  float gntracksmaps = NAN;
1176  float gnembed = NAN;
1177  float nfromtruth = NAN;
1178 
1179  PHG4VtxPoint* point = embedvtxid_vertex[gembed];
1180 
1181  if (point)
1182  {
1183  const int point_id = point->get_id();
1184  gvx = point->get_x();
1185  gvy = point->get_y();
1186  gvz = point->get_z();
1187  gvt = point->get_t();
1188  gembed = truthinfo->isEmbededVtx(point_id);
1189  gntracks = embedvtxid_particle_count[(int) gembed];
1190  if (embedvtxid_maps_particle_count[(int) gembed] > 0 && fabs(gvt) < 2000 && fabs(gvz) < 13.0)
1191  gntracksmaps = embedvtxid_maps_particle_count[(int) gembed];
1192  gnembed = (float) ngembed;
1193  // nfromtruth = vertexeval->get_ntracks_contribution(vertex,point);
1194  }
1195 
1196  if(Verbosity() > 1)
1197  std::cout << " adding vertex data " << std::endl;
1198 
1199  float vertex_data[] = {(float) _ievent,m_fSeed,
1200  vx,
1201  vy,
1202  vz,
1203  ntracks,
1204  gvx,
1205  gvy,
1206  gvz,
1207  gvt,
1208  gembed,
1209  gntracks,
1210  gntracksmaps,
1211  gnembed,
1212  nfromtruth,
1213  nhit_tpc_all,
1214  nhit_tpc_in,
1215  nhit_tpc_mid,
1216  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1217 
1218  _ntp_vertex->Fill(vertex_data);
1219  }
1220  }
1221 
1222  }
1223  if (Verbosity() > 1)
1224  {
1225  _timer->stop();
1226  cout << "vertex time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
1227  }
1228  }
1229 
1230  //-----------------------
1231  // fill the gpoint NTuple
1232  //-----------------------
1233 
1234  if (_ntp_gpoint)
1235  {
1236  if (Verbosity() > 1)
1237  {
1238  cout << "Filling ntp_gpoint " << endl;
1239  _timer->restart();
1240  }
1241  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1242 
1243  if (truthinfo)
1244  {
1245  auto vrange = truthinfo->GetPrimaryVtxRange();
1246  const auto prange = truthinfo->GetPrimaryParticleRange();
1247 
1248  map<int, unsigned int> vertex_particle_count;
1249  for (auto iter = prange.first; iter != prange.second; ++iter) // process all primary paricle
1250  {
1251  ++vertex_particle_count[iter->second->get_vtx_id()];
1252  }
1253 
1254  for (auto iter = vrange.first; iter != vrange.second; ++iter) // process all primary vertexes
1255  {
1256  const int point_id = iter->first;
1257  PHG4VtxPoint* point = iter->second;
1258 
1259  // PHG4VtxPoint* point = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
1260 
1261  if (point)
1262  {
1263  SvtxVertex* vertex = vertexeval->best_vertex_from(point);
1264 
1265  float gvx = point->get_x();
1266  float gvy = point->get_y();
1267  float gvz = point->get_z();
1268  float gvt = point->get_t();
1269  float gntracks = vertex_particle_count[point_id];
1270 
1271  float gembed = truthinfo->isEmbededVtx(point_id);
1272  float vx = NAN;
1273  float vy = NAN;
1274  float vz = NAN;
1275  float ntracks = NAN;
1276  float nfromtruth = NAN;
1277 
1278  if (vertex)
1279  {
1280  vx = vertex->get_x();
1281  vy = vertex->get_y();
1282  vz = vertex->get_z();
1283  ntracks = vertex->size_tracks();
1284  nfromtruth = vertexeval->get_ntracks_contribution(vertex, point);
1285  }
1286 
1287  float gpoint_data[] = {(float) _ievent,m_fSeed,
1288  gvx,
1289  gvy,
1290  gvz,
1291  gvt,
1292  gntracks,
1293  gembed,
1294  vx,
1295  vy,
1296  vz,
1297  ntracks,
1298  nfromtruth,
1299  nhit_tpc_all,
1300  nhit_tpc_in,
1301  nhit_tpc_mid,
1302  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1303 
1304  _ntp_gpoint->Fill(gpoint_data);
1305  }
1306  }
1307  }
1308  if (Verbosity() > 1)
1309  {
1310  _timer->stop();
1311  cout << "gpoint time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
1312  }
1313  }
1314 
1315  //---------------------
1316  // fill the G4hit NTuple
1317  //---------------------
1318 
1319  if (_ntp_g4hit)
1320  {
1321  if (Verbosity() > 1)
1322  {
1323  cout << "Filling ntp_g4hit " << endl;
1324  _timer->restart();
1325  }
1326  std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits();
1327  for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
1328  iter != g4hits.end();
1329  ++iter)
1330  {
1331  PHG4Hit* g4hit = *iter;
1332  PHG4Particle* g4particle = trutheval->get_particle(g4hit);
1333 
1334  float g4hitID = g4hit->get_hit_id();
1335  float gx = g4hit->get_avg_x();
1336  float gy = g4hit->get_avg_y();
1337  float gz = g4hit->get_avg_z();
1338  TVector3 vg4(gx, gy, gz);
1339  float gt = g4hit->get_avg_t();
1340  float gpl = g4hit->get_path_length();
1341  TVector3 vin(g4hit->get_x(0), g4hit->get_y(0), g4hit->get_z(0));
1342  TVector3 vout(g4hit->get_x(1), g4hit->get_y(1), g4hit->get_z(1));
1343  float gdphi = vin.DeltaPhi(vout);
1344  float gdz = fabs(g4hit->get_z(1) - g4hit->get_z(0));
1345  float gedep = g4hit->get_edep();
1346  float glayer = g4hit->get_layer();
1347 
1348  float gtrackID = g4hit->get_trkid();
1349 
1350  float gflavor = NAN;
1351  float gpx = NAN;
1352  float gpy = NAN;
1353  float gpz = NAN;
1354  TVector3 vec(g4hit->get_avg_x(), g4hit->get_avg_y(), g4hit->get_avg_z());
1355  float geta = vec.Eta();
1356  float gphi = vec.Phi();
1357  float gvx = NAN;
1358  float gvy = NAN;
1359  float gvz = NAN;
1360 
1361  float gembed = NAN;
1362  float gprimary = NAN;
1363 
1364  float gfpx = 0.;
1365  float gfpy = 0.;
1366  float gfpz = 0.;
1367  float gfx = 0.;
1368  float gfy = 0.;
1369  float gfz = 0.;
1370 
1371  if (g4particle)
1372  {
1373  if (_scan_for_embedded)
1374  {
1375  if (trutheval->get_embed(g4particle) <= 0) continue;
1376  }
1377 
1378  gflavor = g4particle->get_pid();
1379  gpx = g4particle->get_px();
1380  gpy = g4particle->get_py();
1381  gpz = g4particle->get_pz();
1382 
1383  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1384 
1385  if (vtx)
1386  {
1387  gvx = vtx->get_x();
1388  gvy = vtx->get_y();
1389  gvz = vtx->get_z();
1390  }
1391  PHG4Hit* outerhit = nullptr;
1392  if (_do_eval_light == false)
1393  outerhit = trutheval->get_outermost_truth_hit(g4particle);
1394 
1395  if (outerhit)
1396  {
1397  gfpx = outerhit->get_px(1);
1398  gfpy = outerhit->get_py(1);
1399  gfpz = outerhit->get_pz(1);
1400  gfx = outerhit->get_x(1);
1401  gfy = outerhit->get_y(1);
1402  gfz = outerhit->get_z(1);
1403  }
1404 
1405  gembed = trutheval->get_embed(g4particle);
1406  gprimary = trutheval->is_primary(g4particle);
1407  } // if (g4particle)
1408 
1409  std::set<TrkrDefs::cluskey> clusters = clustereval->all_clusters_from(g4hit);
1410  float nclusters = clusters.size();
1411 
1412  // best cluster reco'd
1413  TrkrDefs::cluskey cluster_key = clustereval->best_cluster_from(g4hit);
1414 
1415  float clusID = NAN;
1416  float x = NAN;
1417  float y = NAN;
1418  float z = NAN;
1419  float eta = NAN;
1420  float phi = NAN;
1421  float e = NAN;
1422  float adc = NAN;
1423  float layer = NAN;
1424  float size = NAN;
1425  float efromtruth = NAN;
1426  float dphitru = NAN;
1427  float detatru = NAN;
1428  float dztru = NAN;
1429  float drtru = NAN;
1430 
1431  TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
1432  if(!clustermap)
1433  clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1434 
1435  TrkrCluster *cluster = clustermap->findCluster(cluster_key);
1436 
1437  if (cluster)
1438  {
1439  clusID = cluster_key;
1440  auto global = actsTransformer.getGlobalPosition(cluster,surfmaps,tgeometry);
1441  x = global(0);
1442  y = global(1);
1443  z = global(2);
1444  TVector3 vec2(x, y, z);
1445  eta = vec2.Eta();
1446  phi = vec2.Phi();
1447  e = cluster->getAdc();
1448  adc = cluster->getAdc();
1449  layer = (float) TrkrDefs::getLayer(cluster_key);
1450  size = 0.0;
1451  // count all hits for this cluster
1452 
1453  TrkrClusterHitAssoc *cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
1454  std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
1455  hitrange = cluster_hit_map->getHits(cluster_key);
1456  for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
1457  clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
1458  {
1459  ++size;
1460  }
1461 
1462  dphitru = vec2.DeltaPhi(vg4);
1463  detatru = eta - geta;
1464  dztru = z - gz;
1465  drtru = vec2.DeltaR(vg4);
1466  if (g4particle)
1467  {
1468  efromtruth = clustereval->get_energy_contribution(cluster_key, g4particle);
1469  }
1470  }
1471 
1472  float g4hit_data[] = {(float) _ievent,m_fSeed,
1473  g4hitID,
1474  gx,
1475  gy,
1476  gz,
1477  gt,
1478  gpl,
1479  gedep,
1480  geta,
1481  gphi,
1482  gdphi,
1483  gdz,
1484  glayer,
1485  gtrackID,
1486  gflavor,
1487  gpx,
1488  gpy,
1489  gpz,
1490  gvx,
1491  gvy,
1492  gvz,
1493  gfpx,
1494  gfpy,
1495  gfpz,
1496  gfx,
1497  gfy,
1498  gfz,
1499  gembed,
1500  gprimary,
1501  nclusters,
1502  clusID,
1503  x,
1504  y,
1505  z,
1506  eta,
1507  phi,
1508  e,
1509  adc,
1510  layer,
1511  size,
1512  efromtruth,
1513  dphitru,
1514  detatru,
1515  dztru,
1516  drtru,
1517  nhit_tpc_all,
1518  nhit_tpc_in,
1519  nhit_tpc_mid,
1520  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1521 
1522  _ntp_g4hit->Fill(g4hit_data);
1523  }
1524  if (Verbosity() > 1)
1525  {
1526  _timer->stop();
1527  cout << "g4hit time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
1528  }
1529  }
1530 
1531  //--------------------
1532  // fill the Hit NTuple
1533  //--------------------
1534 
1535  if (_ntp_hit)
1536  {
1537  if (Verbosity() >= 1)
1538  {
1539  cout << "Filling ntp_hit " << endl;
1540  _timer->restart();
1541  }
1542  // need things off of the DST...
1543  TrkrHitSetContainer *hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1544  PHG4CylinderCellGeomContainer* geom_container =
1545  findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
1546  if (!geom_container)
1547  {
1548  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
1549  return;
1550  }
1551 
1552  if (hitmap)
1553  {
1554  TrkrHitSetContainer::ConstRange all_hitsets = hitmap->getHitSets();
1555  for (TrkrHitSetContainer::ConstIterator iter = all_hitsets.first;
1556  iter != all_hitsets.second;
1557  ++iter)
1558  {
1559  TrkrDefs::hitsetkey hitset_key = iter->first;
1560  TrkrHitSet *hitset = iter->second;
1561 
1562  // get all hits for this hitset
1563  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
1564  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
1565  hitr != hitrangei.second;
1566  ++hitr)
1567  {
1568  TrkrDefs::hitkey hit_key = hitr->first;
1569  TrkrHit *hit = hitr->second;
1570  PHG4Hit* g4hit = hiteval->max_truth_hit_by_energy(hit_key);
1571  PHG4Particle* g4particle = trutheval->get_particle(g4hit);
1572  float event = _ievent;
1573  float hitID = hit_key;
1574  float e = hit->getEnergy();
1575  float adc = hit->getAdc();
1576  float layer = TrkrDefs::getLayer(hitset_key);
1577  float sector = TpcDefs::getSectorId(hitset_key);
1578  float side = TpcDefs::getSide(hitset_key);
1579  float cellID = 0;
1580  float ecell = hit->getAdc();
1581 
1582  float phibin = NAN;
1583  float zbin = NAN;
1584  float phi = NAN;
1585  float z = NAN;
1586 
1587  if (layer >= _nlayers_maps + _nlayers_intt && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc )
1588  {
1589  PHG4CylinderCellGeom* GeoLayer = geom_container->GetLayerCellGeom(layer);
1590  phibin = (float) TpcDefs::getPad(hit_key);
1591  zbin = (float) TpcDefs::getTBin(hit_key);
1592  phi = GeoLayer->get_phicenter(phibin);
1593  z = GeoLayer->get_zcenter(zbin);
1594  }
1595 
1596  float g4hitID = NAN;
1597  float gedep = NAN;
1598  float gx = NAN;
1599  float gy = NAN;
1600  float gz = NAN;
1601  float gt = NAN;
1602  float gtrackID = NAN;
1603  float gflavor = NAN;
1604  float gpx = NAN;
1605  float gpy = NAN;
1606  float gpz = NAN;
1607  float gvx = NAN;
1608  float gvy = NAN;
1609  float gvz = NAN;
1610  float gvt = NAN;
1611  float gfpx = NAN;
1612  float gfpy = NAN;
1613  float gfpz = NAN;
1614  float gfx = NAN;
1615  float gfy = NAN;
1616  float gfz = NAN;
1617  float gembed = NAN;
1618  float gprimary = NAN;
1619 
1620  float efromtruth = NAN;
1621 
1622  if (g4hit)
1623  {
1624  g4hitID = g4hit->get_hit_id();
1625  gedep = g4hit->get_edep();
1626  gx = g4hit->get_avg_x();
1627  gy = g4hit->get_avg_y();
1628  gz = g4hit->get_avg_z();
1629  gt = g4hit->get_avg_t();
1630 
1631  if (g4particle)
1632  {
1633  if (_scan_for_embedded)
1634  {
1635  if (trutheval->get_embed(g4particle) <= 0) continue;
1636  }
1637 
1638  gtrackID = g4particle->get_track_id();
1639  gflavor = g4particle->get_pid();
1640  gpx = g4particle->get_px();
1641  gpy = g4particle->get_py();
1642  gpz = g4particle->get_pz();
1643 
1644  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1645 
1646  if (vtx)
1647  {
1648  gvx = vtx->get_x();
1649  gvy = vtx->get_y();
1650  gvz = vtx->get_z();
1651  gvt = vtx->get_t();
1652  }
1653 
1654  PHG4Hit* outerhit = nullptr;
1655  if (_do_eval_light == false)
1656  outerhit = trutheval->get_outermost_truth_hit(g4particle);
1657  if (outerhit)
1658  {
1659  gfpx = outerhit->get_px(1);
1660  gfpy = outerhit->get_py(1);
1661  gfpz = outerhit->get_pz(1);
1662  gfx = outerhit->get_x(1);
1663  gfy = outerhit->get_y(1);
1664  gfz = outerhit->get_z(1);
1665  }
1666  gembed = trutheval->get_embed(g4particle);
1667  gprimary = trutheval->is_primary(g4particle);
1668  } // if (g4particle){
1669  }
1670 
1671  if (g4particle)
1672  {
1673  efromtruth = hiteval->get_energy_contribution(hit_key, g4particle);
1674  }
1675 
1676  float hit_data[] = {
1677  event,
1678  (float) _iseed,
1679  hitID,
1680  e,
1681  adc,
1682  layer,
1683  sector,
1684  side,
1685  cellID,
1686  ecell,
1687  (float) phibin,
1688  (float) zbin,
1689  phi,
1690  z,
1691  g4hitID,
1692  gedep,
1693  gx,
1694  gy,
1695  gz,
1696  gt,
1697  gtrackID,
1698  gflavor,
1699  gpx,
1700  gpy,
1701  gpz,
1702  gvx,
1703  gvy,
1704  gvz,
1705  gvt,
1706  gfpx,
1707  gfpy,
1708  gfpz,
1709  gfx,
1710  gfy,
1711  gfz,
1712  gembed,
1713  gprimary,
1714  efromtruth,
1715  nhit_tpc_all,
1716  nhit_tpc_in,
1717  nhit_tpc_mid,
1718  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1719 
1720  _ntp_hit->Fill(hit_data);
1721  }
1722  }
1723  }
1724  if (Verbosity() >= 1)
1725  {
1726  _timer->stop();
1727  cout << "hit time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
1728  }
1729  }
1730 
1731  //------------------------
1732  // fill the Cluster NTuple
1733  //------------------------
1734 
1735  if (Verbosity() >= 1)
1736  {
1737  cout << "check for ntp_cluster" << endl;
1738  _timer->restart();
1739  }
1740 
1741  if (_ntp_cluster && !_scan_for_embedded)
1742  {
1743  if (Verbosity() > 1) cout << "Filling ntp_cluster (all of them) " << endl;
1744  // need things off of the DST...
1745  TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
1746  if(!clustermap)
1747  clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1748 
1749  TrkrClusterHitAssoc* clusterhitmap = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
1750  TrkrHitSetContainer* hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1751  TrkrClusterIterationMapv1* _iteration_map = findNode::getClass<TrkrClusterIterationMapv1>(topNode, "CLUSTER_ITERATION_MAP");
1752 
1753  if (Verbosity() > 1){
1754  if (clustermap != nullptr)
1755  cout << "got clustermap" << endl;
1756  else
1757  cout << "no clustermap" << endl;
1758  if (clusterhitmap != nullptr)
1759  cout << "got clusterhitmap" << endl;
1760  else
1761  cout << "no clusterhitmap" << endl;
1762 
1763  if (hitsets != nullptr)
1764  cout << "got hitsets" << endl;
1765  else
1766  cout << "no hitsets" << endl;
1767  }
1768 
1769  if (clustermap != nullptr && clusterhitmap != nullptr && hitsets != nullptr){
1770  auto hitsetrange = hitsets->getHitSets();
1771  for (auto hitsetitr = hitsetrange.first;
1772  hitsetitr != hitsetrange.second;
1773  ++hitsetitr){
1774  int hitsetlayer = TrkrDefs::getLayer(hitsetitr->first);
1775  auto range = clustermap->getClusters(hitsetitr->first);
1776  for( auto iter = range.first; iter != range.second; ++iter ){
1777  TrkrDefs::cluskey cluster_key = iter->first;
1778  TrkrCluster *cluster = clustermap->findCluster(cluster_key);
1779  SvtxTrack* track = trackeval->best_track_from(cluster_key);
1780  PHG4Particle* g4particle = clustereval->max_truth_particle_by_cluster_energy(cluster_key);
1781  float niter = 0;
1782  if(_iteration_map!=NULL)
1783  niter = _iteration_map->getIteration(cluster_key);
1784  float hitID = (float) cluster_key;
1785  auto cglob = actsTransformer.getGlobalPosition(cluster,surfmaps,tgeometry);
1786  float x = cglob(0);
1787  float y = cglob(1);
1788  float z = cglob(2);
1789  TVector3 pos(x, y, z);
1790  float r = pos.Perp();
1791  float phi = pos.Phi();
1792  float eta = pos.Eta();
1793  float theta = pos.Theta();
1794  auto globerr = calculateClusterError(cluster,phi);
1795  float ex = sqrt(globerr[0][0]);
1796  float ey = sqrt(globerr[1][1]);
1797  float ez = cluster->getZError();
1798  float ephi = cluster->getRPhiError();
1799 
1800  float e = cluster->getAdc();
1801  float adc = cluster->getAdc();
1802  float layer = (float) TrkrDefs::getLayer(cluster_key);
1803  float sector = TpcDefs::getSectorId(cluster_key);
1804  float side = TpcDefs::getSide(cluster_key);
1805  float size = 0;
1806  float maxadc = -999;
1807  // count all hits for this cluster
1808  TrkrDefs::hitsetkey hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
1809  int hitsetlayer2 = TrkrDefs::getLayer(hitsetkey);
1810  if(hitsetlayer!=layer){
1811  cout << "WARNING hitset layer " << hitsetlayer << "| " << hitsetlayer2 << " layer " << layer << endl;
1812  }
1813  /*else{
1814  cout << "Good hitset layer " << hitsetlayer << "| " << hitsetlayer2 << " layer " << layer << endl;
1815  }
1816  */
1817  float sumadc = 0;
1818  TrkrHitSetContainer::Iterator hitset = hitsets->findOrAddHitSet(hitsetkey);
1819  std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
1820  hitrange = clusterhitmap->getHits(cluster_key);
1821  for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
1822  clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
1823  {
1824  TrkrHit* hit = hitset->second->getHit(clushititer->second);
1825  ++size;
1826  sumadc += (hit->getAdc() - 70);
1827  if((hit->getAdc()-70)>maxadc)
1828  maxadc = (hit->getAdc()-70);
1829  }
1830  e = sumadc;
1831 
1832  float trackID = NAN;
1833  if (track!=NULL) trackID = track->get_id();
1834 
1835  float g4hitID = NAN;
1836  float gx = NAN;
1837  float gy = NAN;
1838  float gz = NAN;
1839  float gr = NAN;
1840  float gphi = NAN;
1841  //float gedep = NAN;
1842  float geta = NAN;
1843  float gt = NAN;
1844  float gtrackID = NAN;
1845  float gflavor = NAN;
1846  float gpx = NAN;
1847  float gpy = NAN;
1848  float gpz = NAN;
1849  float gvx = NAN;
1850  float gvy = NAN;
1851  float gvz = NAN;
1852  float gvt = NAN;
1853  float gfpx = NAN;
1854  float gfpy = NAN;
1855  float gfpz = NAN;
1856  float gfx = NAN;
1857  float gfy = NAN;
1858  float gfz = NAN;
1859  float gembed = NAN;
1860  float gprimary = NAN;
1861 
1862  float efromtruth = NAN;
1863 
1864  if(Verbosity() > 1)
1865  {
1866  TrkrDefs::cluskey reco_cluskey = cluster->getClusKey();
1867  std::cout << PHWHERE << " **** reco: layer " << layer << std::endl;
1868  cout << " reco cluster key " << reco_cluskey << " r " << r << " x " << x << " y " << y << " z " << z << " phi " << phi << " adc " << adc << endl;
1869  }
1870  float nparticles = NAN;
1871 
1872  // get best matching truth cluster from clustereval
1873  std::shared_ptr<TrkrCluster> truth_cluster = clustereval->max_truth_cluster_by_energy(cluster_key);
1874  if(truth_cluster)
1875  {
1876  if(Verbosity() > 1)
1877  {
1878  TrkrDefs::cluskey truth_cluskey = truth_cluster->getClusKey();
1879  cout << "Found matching truth cluster with key " << truth_cluskey << " for reco cluster key " << cluster_key << " in layer " << layer << endl;
1880  }
1881 
1882  g4hitID = 0;
1883  gx=truth_cluster->getX();
1884  gy=truth_cluster->getY();
1885  gz=truth_cluster->getZ();
1886  efromtruth = truth_cluster->getError(0,0);
1887 
1888  TVector3 gpos(gx, gy, gz);
1889  gr = gpos.Perp(); // could also be just the center of the layer
1890  gphi = gpos.Phi();
1891  geta = gpos.Eta();
1892 
1893  if (g4particle)
1894  {
1895  gtrackID = g4particle->get_track_id();
1896  gflavor = g4particle->get_pid();
1897  gpx = g4particle->get_px();
1898  gpy = g4particle->get_py();
1899  gpz = g4particle->get_pz();
1900 
1901  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1902  if (vtx)
1903  {
1904  gvx = vtx->get_x();
1905  gvy = vtx->get_y();
1906  gvz = vtx->get_z();
1907  gvt = vtx->get_t();
1908  }
1909 
1910  PHG4Hit* outerhit = nullptr;
1911  if (_do_eval_light == false)
1912  outerhit = trutheval->get_outermost_truth_hit(g4particle);
1913  if (outerhit)
1914  {
1915  gfpx = outerhit->get_px(1);
1916  gfpy = outerhit->get_py(1);
1917  gfpz = outerhit->get_pz(1);
1918  gfx = outerhit->get_x(1);
1919  gfy = outerhit->get_y(1);
1920  gfz = outerhit->get_z(1);
1921  }
1922 
1923  gembed = trutheval->get_embed(g4particle);
1924  gprimary = trutheval->is_primary(g4particle);
1925  } // if (g4particle){
1926 
1927  if(Verbosity() > 1)
1928  {
1929  TrkrDefs::cluskey ckey = truth_cluster->getClusKey();
1930  cout << " truth cluster key " << ckey << " gr " << gr << " gx " << gx << " gy " << gy << " gz " << gz << " gphi " << gphi << " efromtruth " << efromtruth << endl;
1931  }
1932  } // if (truth_cluster) {
1933 
1934  if (g4particle)
1935  {
1936 
1937  }
1938  nparticles = clustereval->all_truth_particles(cluster_key).size();
1939 
1940  float cluster_data[] = {(float) _ievent,
1941  (float) _iseed,
1942  hitID,
1943  x,
1944  y,
1945  z,
1946  r,
1947  phi,
1948  eta,
1949  theta,
1950  ex,
1951  ey,
1952  ez,
1953  ephi,
1954  e,
1955  adc,
1956  maxadc,
1957  layer,
1958  sector,
1959  side,
1960  size,
1961  trackID,
1962  niter,
1963  g4hitID,
1964  gx,
1965  gy,
1966  gz,
1967  gr,
1968  gphi,
1969  geta,
1970  gt,
1971  gtrackID,
1972  gflavor,
1973  gpx,
1974  gpy,
1975  gpz,
1976  gvx,
1977  gvy,
1978  gvz,
1979  gvt,
1980  gfpx,
1981  gfpy,
1982  gfpz,
1983  gfx,
1984  gfy,
1985  gfz,
1986  gembed,
1987  gprimary,
1988  efromtruth,
1989  nparticles,
1990  nhit_tpc_all,
1991  nhit_tpc_in,
1992  nhit_tpc_mid,
1993  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1994 
1995  _ntp_cluster->Fill(cluster_data);
1996  }
1997  }
1998  }
1999  }
2000  else if (_ntp_cluster && _scan_for_embedded)
2001  {
2002  if (Verbosity() > 1)
2003  cout << "Filling ntp_cluster (embedded only) " << endl;
2004 
2005  // if only scanning embedded signals, loop over all the tracks from
2006  // embedded particles and report all of their clusters, including those
2007  // from other sources (noise hits on the embedded track)
2008 
2009  // need things off of the DST...
2010  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname.c_str());
2011 
2012  TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
2013  if(!clustermap)
2014  clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
2015 
2016  TrkrClusterHitAssoc* clusterhitmap = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
2017  TrkrHitSetContainer* hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
2018  TrkrClusterIterationMapv1* _iteration_map = findNode::getClass<TrkrClusterIterationMapv1>(topNode, "CLUSTER_ITERATION_MAP");
2019 
2020  if (trackmap != nullptr && clustermap != nullptr && clusterhitmap != nullptr && hitsets != nullptr){
2021  for (SvtxTrackMap::Iter iter = trackmap->begin();
2022  iter != trackmap->end();
2023  ++iter)
2024  {
2025  SvtxTrack* track = iter->second;
2026 
2027  PHG4Particle* truth = trackeval->max_truth_particle_by_nclusters(track);
2028  if (truth)
2029  {
2030  if (trutheval->get_embed(truth) <= 0) continue;
2031  }
2032 
2033  for (SvtxTrack::ConstClusterKeyIter iter = track->begin_cluster_keys();
2034  iter != track->end_cluster_keys();
2035  ++iter)
2036  {
2037  TrkrDefs::cluskey cluster_key = *iter;
2038  TrkrCluster* cluster = clustermap->findCluster(cluster_key);
2039  if(!cluster) continue; // possible to be missing from corrected clusters if cluster mover fails
2040 
2041  PHG4Hit* g4hit = clustereval->max_truth_hit_by_energy(cluster_key);
2042  PHG4Particle* g4particle = trutheval->get_particle(g4hit);
2043 
2044  //float hitID = cluster_key;
2045  float niter = 0;
2046  if(_iteration_map!=NULL)
2047  niter = _iteration_map->getIteration(cluster_key);
2048  float hitID = (float) TrkrDefs::getClusIndex(cluster_key);
2049  auto glob = actsTransformer.getGlobalPosition(cluster,surfmaps,tgeometry);
2050  float x = glob(0);
2051  float y = glob(1);
2052  float z = glob(2);
2053  TVector3 pos(x, y, z);
2054  float r = pos.Perp();
2055  float phi = pos.Phi();
2056  float eta = pos.Eta();
2057  float theta = pos.Theta();
2058  auto globerr = calculateClusterError(cluster,phi);
2059  float ex = sqrt(globerr[0][0]);
2060  float ey = sqrt(globerr[1][1]);
2061  float ez = cluster->getZError();
2062 
2063  float ephi = cluster->getRPhiError();
2064 
2065  float e = cluster->getAdc();
2066  float adc = cluster->getAdc();
2067  float layer = (float) TrkrDefs::getLayer(cluster_key);
2068  float sector = TpcDefs::getSectorId(cluster_key);
2069  float side = TpcDefs::getSide(cluster_key);
2070  // count all hits for this cluster
2071 
2072  float size = 0;
2073  float maxadc = -999;
2074  // count all hits for this cluster
2075  TrkrDefs::hitsetkey hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
2076  TrkrHitSetContainer::Iterator hitset = hitsets->findOrAddHitSet(hitsetkey);
2077  std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
2078  hitrange = clusterhitmap->getHits(cluster_key);
2079  for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
2080  clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
2081  {
2082  TrkrHit* hit = hitset->second->getHit(clushititer->second);
2083  ++size;
2084  if(hit->getAdc()>maxadc)
2085  maxadc = hit->getAdc();
2086  }
2087 
2088  float trackID = NAN;
2089  trackID = track->get_id();
2090 
2091  float g4hitID = NAN;
2092  float gx = NAN;
2093  float gy = NAN;
2094  float gz = NAN;
2095  float gr = NAN;
2096  float gphi = NAN;
2097  //float gedep = NAN;
2098  float geta = NAN;
2099  float gt = NAN;
2100  float gtrackID = NAN;
2101  float gflavor = NAN;
2102  float gpx = NAN;
2103  float gpy = NAN;
2104  float gpz = NAN;
2105  float gvx = NAN;
2106  float gvy = NAN;
2107  float gvz = NAN;
2108  float gvt = NAN;
2109  float gfpx = NAN;
2110  float gfpy = NAN;
2111  float gfpz = NAN;
2112  float gfx = NAN;
2113  float gfy = NAN;
2114  float gfz = NAN;
2115  float gembed = NAN;
2116  float gprimary = NAN;
2117 
2118  float efromtruth = NAN;
2119 
2120  //cout << "Look for truth cluster to match reco cluster " << cluster_key << endl;
2121 
2122  // get best matching truth cluster from clustereval
2123  std::shared_ptr<TrkrCluster> truth_cluster = clustereval->max_truth_cluster_by_energy(cluster_key);
2124  if(truth_cluster)
2125  {
2126  if(Verbosity() > 1)
2127  {
2128  TrkrDefs::cluskey truth_cluskey = truth_cluster->getClusKey();
2129  cout << " Found matching truth cluster with key " << truth_cluskey << " for reco cluster key " << cluster_key << " in layer " << layer << endl;
2130  }
2131 
2132  g4hitID = 0;
2133  gx=truth_cluster->getX();
2134  gy=truth_cluster->getY();
2135  gz=truth_cluster->getZ();
2136  efromtruth = truth_cluster->getError(0,0);
2137 
2138  TVector3 gpos(gx, gy, gz);
2139  gr = gpos.Perp();
2140  gphi = gpos.Phi();
2141  geta = gpos.Eta();
2142 
2143  if (g4particle)
2144  {
2145  gtrackID = g4particle->get_track_id();
2146  gflavor = g4particle->get_pid();
2147  gpx = g4particle->get_px();
2148  gpy = g4particle->get_py();
2149  gpz = g4particle->get_pz();
2150 
2151  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
2152  if (vtx)
2153  {
2154  gvx = vtx->get_x();
2155  gvy = vtx->get_y();
2156  gvz = vtx->get_z();
2157  gvt = vtx->get_t();
2158  }
2159  PHG4Hit* outerhit = nullptr;
2160  if (_do_eval_light == false)
2161  outerhit = trutheval->get_outermost_truth_hit(g4particle);
2162  if (outerhit)
2163  {
2164  gfpx = outerhit->get_px(1);
2165  gfpy = outerhit->get_py(1);
2166  gfpz = outerhit->get_pz(1);
2167  gfx = outerhit->get_x(1);
2168  gfy = outerhit->get_y(1);
2169  gfz = outerhit->get_z(1);
2170  }
2171 
2172  gembed = trutheval->get_embed(g4particle);
2173  gprimary = trutheval->is_primary(g4particle);
2174  } // if (g4particle){
2175  } // if (g4hit) {
2176 
2177  float nparticles = clustereval->all_truth_particles(cluster_key).size();
2178 
2179  float cluster_data[] = {(float) _ievent,
2180  (float) _iseed,
2181  hitID,
2182  x,
2183  y,
2184  z,
2185  r,
2186  phi,
2187  eta,
2188  theta,
2189  ex,
2190  ey,
2191  ez,
2192  ephi,
2193  e,
2194  adc,
2195  maxadc,
2196  layer,
2197  sector,
2198  side,
2199  size,
2200  trackID,
2201  niter,
2202  g4hitID,
2203  gx,
2204  gy,
2205  gz,
2206  gr,
2207  gphi,
2208  geta,
2209  gt,
2210  gtrackID,
2211  gflavor,
2212  gpx,
2213  gpy,
2214  gpz,
2215  gvx,
2216  gvy,
2217  gvz,
2218  gvt,
2219  gfpx,
2220  gfpy,
2221  gfpz,
2222  gfx,
2223  gfy,
2224  gfz,
2225  gembed,
2226  gprimary,
2227  efromtruth,
2228  nparticles,
2229  nhit_tpc_all,
2230  nhit_tpc_in,
2231  nhit_tpc_mid,
2232  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
2233 
2234  _ntp_cluster->Fill(cluster_data);
2235  }
2236  }
2237  }
2238  }
2239  if (Verbosity() >= 1){
2240  _timer->stop();
2241  cout << "cluster time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
2242  }
2243 
2244  // fill the truth cluster NTuple
2245  //-----------------------------------
2246 
2247  if (Verbosity() > 1)
2248  {
2249  cout << "check for ntp_g4cluster" << endl;
2250  _timer->restart();
2251  }
2252 
2253  if (_ntp_g4cluster)
2254  {
2255  if (Verbosity() >= 1)
2256  cout << "Filling ntp_g4cluster " << endl;
2257 
2258  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
2259  PHG4TruthInfoContainer::ConstRange range = truthinfo->GetParticleRange();
2260  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
2261  iter != range.second;
2262  ++iter)
2263  {
2264 
2265  PHG4Particle* g4particle = iter->second;
2266 
2267  if (_scan_for_embedded)
2268  {
2269  if (trutheval->get_embed(g4particle) <= 0) continue;
2270  }
2271 
2272  float gtrackID = g4particle->get_track_id();
2273  float gflavor = g4particle->get_pid();
2274  float gembed = trutheval->get_embed(g4particle);
2275  float gprimary = trutheval->is_primary(g4particle);
2276 
2277  if(Verbosity() > 1)
2278  cout << PHWHERE << " PHG4Particle ID " << gtrackID << " gflavor " << gflavor << " gprimary " << gprimary << endl;
2279 
2280  // Get the truth clusters from this particle
2281  std::map<unsigned int, std::shared_ptr<TrkrCluster> > truth_clusters = trutheval->all_truth_clusters(g4particle);
2282 
2283  // loop over layers and add to ntuple
2284  for ( auto it = truth_clusters.begin(); it != truth_clusters.end(); ++it )
2285  {
2286  unsigned int layer = it->first;
2287  std::shared_ptr<TrkrCluster> gclus = it->second;
2288 
2289  float gx = gclus->getX();
2290  float gy = gclus->getY();
2291  float gz = gclus->getZ();
2292  float gt = NAN;
2293  float gedep = gclus->getError(0,0);
2294  float gadc = (float) gclus->getAdc();
2295 
2296  TVector3 gpos(gx, gy, gz);
2297  float gr = sqrt(gx*gx+gy*gy);
2298  float gphi = gpos.Phi();
2299  float geta = gpos.Eta();
2300 
2301  if(Verbosity() > 1)
2302  {
2303  TrkrDefs::cluskey ckey = gclus->getClusKey();
2304  std::cout << PHWHERE << " **** truth: layer " << layer << std::endl;
2305  cout << " truth cluster key " << ckey << " gr " << gr << " gx " << gx << " gy " << gy << " gz " << gz
2306  << " gphi " << gphi << " gedep " << gedep << " gadc " << gadc << endl;
2307  }
2308 
2309  float gphisize = gclus->getSize(1,1);
2310  float gzsize = gclus->getSize(2,2);
2311 
2312  // Find the matching TrkrCluster, if it exists
2313 
2314  float x = NAN;
2315  float y = NAN;
2316  float z = NAN;
2317  float r = NAN;
2318  float phi = NAN;
2319  float eta = NAN;
2320  float ex = NAN;
2321  float ey = NAN;
2322  float ez = NAN;
2323  float ephi = NAN;
2324  float adc = NAN;
2325 
2326  float nreco = 0;
2327 
2328  TrkrCluster *reco_cluster = clustereval->reco_cluster_from_truth_cluster(gclus);
2329  if(reco_cluster)
2330  {
2331  nreco = 1;
2332  auto glob = actsTransformer.getGlobalPosition(reco_cluster,surfmaps,tgeometry);
2333  x = glob(0);
2334  y = glob(1);
2335  z = glob(2);
2336 
2337  TVector3 pos(x, y, z);
2338  r = sqrt(x*x+y*y);
2339  phi = pos.Phi();
2340  eta = pos.Eta();
2341  auto globerr = calculateClusterError(reco_cluster,phi);
2342  ex = sqrt(globerr[0][0]);
2343  ey = sqrt(globerr[1][1]);
2344  ez = reco_cluster->getZError();
2345  ephi = reco_cluster->getRPhiError();
2346 
2347  adc = reco_cluster->getAdc();
2348 
2349  if(Verbosity() > 1)
2350  {
2351  TrkrDefs::cluskey reco_cluskey = reco_cluster->getClusKey();
2352  cout << " reco cluster key " << reco_cluskey << " r " << r << " x " << x << " y " << y << " z " << z << " phi " << phi << " adc " << adc << endl;
2353  }
2354  }
2355  if(nreco == 0 && Verbosity() > 1)
2356  {
2357  if(Verbosity() > 1)
2358  cout << " ----------- Failed to find matching reco cluster " << endl;
2359  }
2360 
2361  // add this cluster to the ntuple
2362 
2363  float g4cluster_data[] = {(float) _ievent,
2364  (float) layer,
2365  gx,
2366  gy,
2367  gz,
2368  gt,
2369  gedep,
2370  gr,
2371  gphi,
2372  geta,
2373  gtrackID,
2374  gflavor,
2375  gembed,
2376  gprimary,
2377  gphisize,
2378  gzsize,
2379  gadc,
2380  nreco,
2381  x,
2382  y,
2383  z,
2384  r,
2385  phi,
2386  eta,
2387  ex,
2388  ey,
2389  ez,
2390  ephi,
2391  adc };
2392  _ntp_g4cluster->Fill(g4cluster_data);
2393  }
2394  }
2395  }
2396 
2397  //------------------------
2398  // fill the Gtrack NTuple
2399  //------------------------
2400 
2401  // need things off of the DST...
2402 
2403  //cout << "check for ntp_gtrack" << endl;
2404 
2405  if (_ntp_gtrack)
2406  {
2407  if (Verbosity() > 1)
2408  {
2409  cout << "Filling ntp_gtrack " << endl;
2410  _timer->restart();
2411  }
2412 
2413  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
2414  if (truthinfo)
2415  {
2416 
2417  PHG4TruthInfoContainer::ConstRange range = truthinfo->GetParticleRange();
2418  if(_scan_for_primaries){
2419  range = truthinfo->GetPrimaryParticleRange();
2420  }
2421 
2422  Float_t gntracks = (Float_t) truthinfo->GetNumPrimaryVertexParticles();
2423  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
2424  iter != range.second;
2425  ++iter)
2426  {
2427 
2428  PHG4Particle* g4particle = iter->second;
2429 
2430  if (_scan_for_embedded)
2431  {
2432  if (trutheval->get_embed(g4particle) <= 0) continue;
2433  }
2434 
2435  float gtrackID = g4particle->get_track_id();
2436  float gflavor = g4particle->get_pid();
2437 
2438  std::set<TrkrDefs::cluskey> g4clusters = clustereval->all_clusters_from(g4particle);
2439 
2440  float ng4hits = g4clusters.size();
2441  unsigned int ngmaps = 0;
2442  unsigned int ngmms = 0;
2443  unsigned int ngintt = 0;
2444  unsigned int ngintt1 = 0;
2445  unsigned int ngintt2 = 0;
2446  unsigned int ngintt3 = 0;
2447  unsigned int ngintt4 = 0;
2448  unsigned int ngintt5 = 0;
2449  unsigned int ngintt6 = 0;
2450  unsigned int ngintt7 = 0;
2451  unsigned int ngintt8 = 0;
2452  unsigned int ngtpc = 0;
2453  unsigned int nglmaps = 0;
2454  unsigned int nglintt = 0;
2455  unsigned int ngltpc = 0;
2456  unsigned int nglmms = 0;
2457 
2458  int lmaps[_nlayers_maps + 1];
2459  if (_nlayers_maps > 0)
2460  for (unsigned int i = 0; i < _nlayers_maps; i++) lmaps[i] = 0;
2461 
2462  int lintt[_nlayers_intt + 1];
2463  if (_nlayers_intt > 0)
2464  for (unsigned int i = 0; i < _nlayers_intt; i++) lintt[i] = 0;
2465 
2466  int ltpc[_nlayers_tpc + 1];
2467  if (_nlayers_tpc > 0)
2468  for (unsigned int i = 0; i < _nlayers_tpc; i++) ltpc[i] = 0;
2469 
2470  int lmms[_nlayers_mms + 1];
2471  if (_nlayers_mms > 0)
2472  for (unsigned int i = 0; i < _nlayers_mms; i++) lmms[i] = 0;
2473 
2474  for (const TrkrDefs::cluskey g4cluster : g4clusters)
2475  {
2476  unsigned int layer = TrkrDefs::getLayer(g4cluster);
2477  //cout<<__LINE__<<": " << _ievent <<": " <<gtrackID << ": " << layer <<": " <<g4cluster->get_id() <<endl;
2478  if (_nlayers_maps > 0 && layer < _nlayers_maps)
2479  {
2480  lmaps[layer] = 1;
2481  ngmaps++;
2482  }
2483  if (_nlayers_mms > 0 && layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
2484  {
2485  lmms[layer - _nlayers_tpc - _nlayers_intt - _nlayers_maps] = 1;
2486  ngmms++;
2487  }
2488  if (_nlayers_intt > 0 && layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt)
2489  {
2490  lintt[layer - _nlayers_maps] = 1;
2491  ngintt++;
2492  }
2493 
2494  if (_nlayers_intt > 0 && layer == _nlayers_maps && layer < _nlayers_maps + _nlayers_intt)
2495  {
2496  ngintt1++;
2497  }
2498 
2499  if (_nlayers_intt > 1 && layer == _nlayers_maps + 1 && layer < _nlayers_maps + _nlayers_intt)
2500  {
2501  ngintt2++;
2502  }
2503 
2504  if (_nlayers_intt > 2 && layer == _nlayers_maps + 2 && layer < _nlayers_maps + _nlayers_intt)
2505  {
2506  ngintt3++;
2507  }
2508 
2509  if (_nlayers_intt > 3 && layer == _nlayers_maps + 3 && layer < _nlayers_maps + _nlayers_intt)
2510  {
2511  ngintt4++;
2512  }
2513 
2514  if (_nlayers_intt > 4 && layer == _nlayers_maps + 4 && layer < _nlayers_maps + _nlayers_intt)
2515  {
2516  ngintt5++;
2517  }
2518 
2519  if (_nlayers_intt > 5 && layer == _nlayers_maps + 5 && layer < _nlayers_maps + _nlayers_intt)
2520  {
2521  ngintt6++;
2522  }
2523 
2524  if (_nlayers_intt > 6 && layer == _nlayers_maps + 6 && layer < _nlayers_maps + _nlayers_intt)
2525  {
2526  ngintt7++;
2527  }
2528 
2529  if (_nlayers_intt > 7 && layer == _nlayers_maps + 7 && layer < _nlayers_maps + _nlayers_intt)
2530  {
2531  ngintt8++;
2532  }
2533  if (_nlayers_tpc > 0 && layer >= _nlayers_maps + _nlayers_intt && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
2534  {
2535  ltpc[layer - (_nlayers_maps + _nlayers_intt)] = 1;
2536  ngtpc++;
2537  }
2538  }
2539  if (_nlayers_maps > 0)
2540  for (unsigned int i = 0; i < _nlayers_maps; i++) nglmaps += lmaps[i];
2541  if (_nlayers_intt > 0)
2542  for (unsigned int i = 0; i < _nlayers_intt; i++) nglintt += lintt[i];
2543  if (_nlayers_tpc > 0)
2544  for (unsigned int i = 0; i < _nlayers_tpc; i++) ngltpc += ltpc[i];
2545  if (_nlayers_mms > 0)
2546  for (unsigned int i = 0; i < _nlayers_mms; i++) nglmms += lmms[i];
2547 
2548  float gpx = g4particle->get_px();
2549  float gpy = g4particle->get_py();
2550  float gpz = g4particle->get_pz();
2551  float gpt = NAN;
2552  float geta = NAN;
2553  float gphi = NAN;
2554  if (gpx != 0 && gpy != 0)
2555  {
2556  TVector3 gv(gpx, gpy, gpz);
2557  gpt = gv.Pt();
2558  geta = gv.Eta();
2559  gphi = gv.Phi();
2560  }
2561  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
2562  float gvx = vtx->get_x();
2563  float gvy = vtx->get_y();
2564  float gvz = vtx->get_z();
2565  float gvt = vtx->get_t();
2566 
2567  float gfpx = 0.;
2568  float gfpy = 0.;
2569  float gfpz = 0.;
2570  float gfx = 0.;
2571  float gfy = 0.;
2572  float gfz = 0.;
2573 
2574  PHG4Hit* outerhit = nullptr;
2575  if (_do_eval_light == false)
2576  outerhit = trutheval->get_outermost_truth_hit(g4particle);
2577 
2578  if (outerhit)
2579  {
2580  gfpx = outerhit->get_px(1);
2581  gfpy = outerhit->get_py(1);
2582  gfpz = outerhit->get_pz(1);
2583  gfx = outerhit->get_x(1);
2584  gfy = outerhit->get_y(1);
2585  gfz = outerhit->get_z(1);
2586  }
2587 
2588  float gembed = trutheval->get_embed(g4particle);
2589  float gprimary = trutheval->is_primary(g4particle);
2590 
2591  float trackID = NAN;
2592  float charge = NAN;
2593  float quality = NAN;
2594  float chisq = NAN;
2595  float ndf = NAN;
2596  float nhits = NAN;
2597  float nmaps = 0;
2598  float nintt = 0;
2599  float ntpc = 0;
2600  float nmms = 0;
2601  float ntpc1 = 0;
2602  float ntpc11 = 0;
2603  float ntpc2 = 0;
2604  float ntpc3 = 0;
2605  float nlintt = 0;
2606  float nlmaps = 0;
2607  float nltpc = 0;
2608  float nlmms = 0;
2609  unsigned int layers = 0x0;
2610  float dca2d = NAN;
2611  float dca2dsigma = NAN;
2612  float dca3dxy = NAN;
2613  float dca3dxysigma = NAN;
2614  float dca3dz = NAN;
2615  float dca3dzsigma = NAN;
2616  float px = NAN;
2617  float py = NAN;
2618  float pz = NAN;
2619  float pt = NAN;
2620  float eta = NAN;
2621  float phi = NAN;
2622  float deltapt = NAN;
2623  float deltaeta = NAN;
2624  float deltaphi = NAN;
2625  float pcax = NAN;
2626  float pcay = NAN;
2627  float pcaz = NAN;
2628 
2629  float nfromtruth = NAN;
2630  float nwrong = NAN;
2631  float ntrumaps = NAN;
2632  float ntruintt = NAN;
2633  float ntrumms = NAN;
2634  float ntrutpc = NAN;
2635  float ntrutpc1 = NAN;
2636  float ntrutpc11 = NAN;
2637  float ntrutpc2 = NAN;
2638  float ntrutpc3 = NAN;
2639  float layersfromtruth = NAN;
2640 
2641  if (_do_track_match)
2642  {
2643  SvtxTrack* track = trackeval->best_track_from(g4particle);
2644 
2645  if (track)
2646  {
2647  trackID = track->get_id();
2648  charge = track->get_charge();
2649  quality = track->get_quality();
2650  chisq = track->get_chisq();
2651  ndf = track->get_ndf();
2652  nhits = track->size_cluster_keys();
2653 
2654  vector <int> maps(_nlayers_maps, 0);
2655  vector <int> intt(_nlayers_intt, 0);
2656  vector <int> tpc(_nlayers_tpc, 0);
2657  vector <int> mms(_nlayers_mms, 0);
2658 
2659  if (_nlayers_maps > 0)
2660  {
2661  for (unsigned int i = 0; i < _nlayers_maps; i++) maps[i] = 0;
2662  }
2663  if (_nlayers_intt > 0)
2664  {
2665  for (unsigned int i = 0; i < _nlayers_intt; i++) intt[i] = 0;
2666  }
2667  if (_nlayers_tpc > 0)
2668  {
2669  for (unsigned int i = 0; i < _nlayers_tpc; i++) tpc[i] = 0;
2670  }
2671  if (_nlayers_mms > 0)
2672  {
2673  for (unsigned int i = 0; i < _nlayers_mms; i++) mms[i] = 0;
2674  }
2675 
2676  for (SvtxTrack::ConstClusterKeyIter iter = track->begin_cluster_keys();
2677  iter != track->end_cluster_keys();
2678  ++iter)
2679  {
2680  TrkrDefs::cluskey cluster_key = *iter;
2681  //TrkrCluster* cluster = clustermap->findCluster(cluster_key);
2682  unsigned int layer = TrkrDefs::getLayer(cluster_key);
2683  if (_nlayers_maps > 0 && layer < _nlayers_maps)
2684  {
2685  maps[layer] = 1;
2686  nmaps++;
2687  }
2688  if (_nlayers_intt > 0 && layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt)
2689  {
2690  intt[layer - _nlayers_maps] = 1;
2691  nintt++;
2692  }
2693  if (_nlayers_tpc > 0 &&
2694  layer >= (_nlayers_maps + _nlayers_intt) &&
2695  layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
2696  {
2697  tpc[layer - (_nlayers_maps + _nlayers_intt)] = 1;
2698  ntpc++;
2699 
2700  if((layer - (_nlayers_maps + _nlayers_intt))<8){
2701  ntpc11++;
2702  }
2703 
2704  if((layer - (_nlayers_maps + _nlayers_intt))<16){
2705  //std::cout << " tpc1: layer " << layer << std::endl;
2706  ntpc1++;
2707  }
2708  else if((layer - (_nlayers_maps + _nlayers_intt))<32){
2709  //std::cout << " tpc2: layer " << layer << std::endl;
2710  ntpc2++;
2711  }
2712  else if((layer - (_nlayers_maps + _nlayers_intt))<48){
2713  //std::cout << " tpc3: layer " << layer << std::endl;
2714  ntpc3++;
2715  }
2716  }
2717  if (_nlayers_mms > 0 && layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
2718  {
2719  mms[layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
2720  nmms++;
2721  }
2722  }
2723  if (_nlayers_maps > 0)
2724  for (unsigned int i = 0; i < _nlayers_maps; i++) nlmaps += maps[i];
2725  if (_nlayers_intt > 0)
2726  for (unsigned int i = 0; i < _nlayers_intt; i++) nlintt += intt[i];
2727  if (_nlayers_tpc > 0)
2728  for (unsigned int i = 0; i < _nlayers_tpc; i++) nltpc += tpc[i];
2729  if (_nlayers_mms > 0)
2730  for (unsigned int i = 0; i < _nlayers_mms; i++) nlmms += mms[i];
2731 
2732  layers = nlmaps + nlintt + nltpc + nlmms;
2733  /* cout << " layers " << layers
2734  << " nmaps " << nmaps
2735  << " nintt " << nintt
2736  << " ntpc " << ntpc
2737  << " nlmaps "<< nlmaps
2738  << " nlintt " << nlintt
2739  << " nltpc " << nltpc
2740  << endl;
2741  */
2742  dca2d = track->get_dca2d();
2743  dca2dsigma = track->get_dca2d_error();
2744  dca3dxy = track->get_dca3d_xy();
2745  dca3dxysigma = track->get_dca3d_xy_error();
2746  dca3dz = track->get_dca3d_z();
2747  dca3dzsigma = track->get_dca3d_z_error();
2748  px = track->get_px();
2749  py = track->get_py();
2750  pz = track->get_pz();
2751  TVector3 v(px, py, pz);
2752  pt = v.Pt();
2753  eta = v.Eta();
2754  phi = v.Phi();
2755  float CVxx = track->get_error(3,3);
2756  float CVxy = track->get_error(3,4);
2757  float CVxz = track->get_error(3,5);
2758  float CVyy = track->get_error(4,4);
2759  float CVyz = track->get_error(4,5);
2760  float CVzz = track->get_error(5,5);
2761  deltapt = sqrt((CVxx*px*px+2*CVxy*px*py+CVyy*py*py)/(px*px+py*py));
2762  deltaeta = sqrt((CVzz*(px*px+py*py)*(px*px+py*py)+pz*(-2*(CVxz*px+CVyz*py)*(px*px+py*py)+CVxx*px*px*pz+CVyy*py*py*pz+2*CVxy*px*py*pz))/((px*px+py*py)*(px*px+py*py)*(px*px+py*py+pz*pz)));
2763  deltaphi = sqrt((CVyy*px*px-2*CVxy*px*py+CVxx*py*py)/((px*px+py*py)*(px*px+py*py)));
2764  pcax = track->get_x();
2765  pcay = track->get_y();
2766  pcaz = track->get_z();
2767 
2768  nfromtruth = trackeval->get_nclusters_contribution(track, g4particle);
2769  nwrong = trackeval->get_nwrongclusters_contribution(track, g4particle);
2770 
2771  if (_nlayers_maps == 0)
2772  {
2773  ntrumaps = 0;
2774  }
2775  else
2776  {
2777  ntrumaps = trackeval->get_layer_range_contribution(track, g4particle, 0, _nlayers_maps);
2778  }
2779  if (_nlayers_intt == 0)
2780  {
2781  ntruintt = 0;
2782  }
2783  else
2784  {
2785  ntruintt = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps, _nlayers_maps + _nlayers_intt);
2786  }
2787  if (_nlayers_mms == 0)
2788  {
2789  ntrumms = 0;
2790  }
2791  else
2792  {
2793  ntrumms = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + _nlayers_tpc, _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms);
2794  }
2795  ntrutpc = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
2796  ntrutpc1 = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 16);
2797  ntrutpc11 = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 8);
2798  ntrutpc2 = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt+16, _nlayers_maps + _nlayers_intt + 32);
2799  ntrutpc3 = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt+32, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
2800 
2801  layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track, g4particle);
2802  }
2803  }
2804 
2805  float gtrack_data[] = {(float) _ievent,m_fSeed,
2806  gntracks,
2807  gtrackID,
2808  gflavor,
2809  ng4hits,
2810  (float) ngmaps,
2811  (float) ngintt,
2812  (float) ngmms,
2813  (float) ngintt1,
2814  (float) ngintt2,
2815  (float) ngintt3,
2816  (float) ngintt4,
2817  (float) ngintt5,
2818  (float) ngintt6,
2819  (float) ngintt7,
2820  (float) ngintt8,
2821  (float) ngtpc,
2822  (float) nglmaps,
2823  (float) nglintt,
2824  (float) ngltpc,
2825  (float) nglmms,
2826  gpx,
2827  gpy,
2828  gpz,
2829  gpt,
2830  geta,
2831  gphi,
2832  gvx,
2833  gvy,
2834  gvz,
2835  gvt,
2836  gfpx,
2837  gfpy,
2838  gfpz,
2839  gfx,
2840  gfy,
2841  gfz,
2842  gembed,
2843  gprimary,
2844  trackID,
2845  px,
2846  py,
2847  pz,
2848  pt,
2849  eta,
2850  phi,
2851  deltapt,
2852  deltaeta,
2853  deltaphi,
2854  charge,
2855  quality,
2856  chisq,
2857  ndf,
2858  nhits,
2859  (float) layers,
2860  nmaps,
2861  nintt,
2862  ntpc,
2863  nmms,
2864  ntpc1,
2865  ntpc11,
2866  ntpc2,
2867  ntpc3,
2868  nlmaps,
2869  nlintt,
2870  nltpc,
2871  nlmms,
2872  dca2d,
2873  dca2dsigma,
2874  dca3dxy,
2875  dca3dxysigma,
2876  dca3dz,
2877  dca3dzsigma,
2878  pcax,
2879  pcay,
2880  pcaz,
2881  nfromtruth,
2882  nwrong,
2883  ntrumaps,
2884  ntruintt,
2885  ntrutpc,
2886  ntrumms,
2887  ntrutpc1,
2888  ntrutpc11,
2889  ntrutpc2,
2890  ntrutpc3,
2891  layersfromtruth,
2892  nhit_tpc_all,
2893  nhit_tpc_in,
2894  nhit_tpc_mid,
2895  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
2896 
2897  /*
2898  cout << " ievent " << _ievent
2899  << " gtrackID " << gtrackID
2900  << " gflavor " << gflavor
2901  << " ng4hits " << ng4hits
2902  << endl;
2903  */
2904 
2905  _ntp_gtrack->Fill(gtrack_data);
2906 
2907  }
2908 
2909  }
2910  if (Verbosity() > 1)
2911  {
2912  _timer->stop();
2913  cout << "gtrack time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
2914  }
2915  }
2916 
2917  //------------------------
2918  // fill the Track NTuple
2919  //------------------------
2920 
2921  if (_ntp_track)
2922  {
2923  if (Verbosity() > 1)
2924  {
2925  cout << "Filling ntp_track " << endl;
2926  _timer->restart();
2927  }
2928 
2929  // need things off of the DST...
2930  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname.c_str());
2931  if (trackmap)
2932  {
2933  for (SvtxTrackMap::Iter iter = trackmap->begin();
2934  iter != trackmap->end();
2935  ++iter)
2936  {
2937  SvtxTrack* track = iter->second;
2938  float trackID = track->get_id();
2939  float charge = track->get_charge();
2940  float quality = track->get_quality();
2941  float chisq = track->get_chisq();
2942  float ndf = track->get_ndf();
2943  float nhits = track->size_cluster_keys();
2944  unsigned int layers = 0x0;
2945  int maps[_nlayers_maps];
2946  int intt[_nlayers_intt];
2947  int tpc[_nlayers_tpc];
2948  int mms[_nlayers_mms];
2949  if (_nlayers_maps > 0)
2950  {
2951  for (unsigned int i = 0; i < _nlayers_maps; i++) maps[i] = 0;
2952  }
2953  if (_nlayers_intt > 0)
2954  {
2955  for (unsigned int i = 0; i < _nlayers_intt; i++) intt[i] = 0;
2956  }
2957  if (_nlayers_tpc > 0)
2958  {
2959  for (unsigned int i = 0; i < _nlayers_tpc; i++) tpc[i] = 0;
2960  }
2961  if (_nlayers_mms > 0)
2962  {
2963  for (unsigned int i = 0; i < _nlayers_mms; i++) mms[i] = 0;
2964  }
2965 
2966  float nmaps = 0;
2967  float nintt = 0;
2968  float nmms = 0;
2969  float ntpc = 0;
2970  float ntpc1 = 0;
2971  float ntpc11 = 0;
2972  float ntpc2 = 0;
2973  float ntpc3 = 0;
2974  float nlmaps = 0;
2975  float nlintt = 0;
2976  float nltpc = 0;
2977  float nlmms = 0;
2978 
2979  for (SvtxTrack::ConstClusterKeyIter iter = track->begin_cluster_keys();
2980  iter != track->end_cluster_keys();
2981  ++iter)
2982  {
2983  TrkrDefs::cluskey cluster_key = *iter;
2984  //TrkrCluster* cluster = clustermap->findCluster(cluster_key);
2985  unsigned int layer = TrkrDefs::getLayer(cluster_key);
2986 
2987  if (_nlayers_maps > 0 && layer < _nlayers_maps)
2988  {
2989  maps[layer] = 1;
2990  nmaps++;
2991  }
2992  if (_nlayers_intt > 0 && layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt)
2993  {
2994  intt[layer - _nlayers_maps] = 1;
2995  nintt++;
2996  }
2997  if (_nlayers_mms > 0 && layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
2998  {
2999  mms[layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3000  nmms++;
3001  }
3002  if (_nlayers_tpc > 0 && layer >= (_nlayers_maps + _nlayers_intt) && layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
3003  {
3004  tpc[layer - (_nlayers_maps + _nlayers_intt)] = 1;
3005  ntpc++;
3006  if((layer - (_nlayers_maps + _nlayers_intt))<8){
3007  ntpc11++;
3008  }
3009 
3010  if((layer - (_nlayers_maps + _nlayers_intt))<16){
3011  ntpc1++;
3012  }
3013  else if((layer - (_nlayers_maps + _nlayers_intt))<32){
3014  ntpc2++;
3015  }
3016  else if((layer - (_nlayers_maps + _nlayers_intt))<48){
3017  ntpc3++;
3018  }
3019  }
3020  }
3021  if (_nlayers_maps > 0)
3022  for (unsigned int i = 0; i < _nlayers_maps; i++) nlmaps += maps[i];
3023  if (_nlayers_intt > 0)
3024  for (unsigned int i = 0; i < _nlayers_intt; i++) nlintt += intt[i];
3025  if (_nlayers_tpc > 0)
3026  for (unsigned int i = 0; i < _nlayers_tpc; i++) nltpc += tpc[i];
3027  if (_nlayers_mms > 0)
3028  for (unsigned int i = 0; i < _nlayers_mms; i++) nlmms += mms[i];
3029  layers = nlmaps + nlintt + nltpc + nlmms;
3030  float dca2d = track->get_dca2d();
3031  float dca2dsigma = track->get_dca2d_error();
3032  float dca3dxy = track->get_dca3d_xy();
3033  float dca3dxysigma = track->get_dca3d_xy_error();
3034  float dca3dz = track->get_dca3d_z();
3035  float dca3dzsigma = track->get_dca3d_z_error();
3036  float px = track->get_px();
3037  float py = track->get_py();
3038  float pz = track->get_pz();
3039  TVector3 v(px, py, pz);
3040  float pt = v.Pt();
3041  float eta = v.Eta();
3042  float phi = v.Phi();
3043  float CVxx = track->get_error(3,3);
3044  float CVxy = track->get_error(3,4);
3045  float CVxz = track->get_error(3,5);
3046  float CVyy = track->get_error(4,4);
3047  float CVyz = track->get_error(4,5);
3048  float CVzz = track->get_error(5,5);
3049  float deltapt = sqrt((CVxx*px*px+2*CVxy*px*py+CVyy*py*py)/(px*px+py*py));
3050  float deltaeta = sqrt((CVzz*(px*px+py*py)*(px*px+py*py)+pz*(-2*(CVxz*px+CVyz*py)*(px*px+py*py)+CVxx*px*px*pz+CVyy*py*py*pz+2*CVxy*px*py*pz))/((px*px+py*py)*(px*px+py*py)*(px*px+py*py+pz*pz)));
3051  float deltaphi = sqrt((CVyy*px*px-2*CVxy*px*py+CVxx*py*py)/((px*px+py*py)*(px*px+py*py)));
3052  float pcax = track->get_x();
3053  float pcay = track->get_y();
3054  float pcaz = track->get_z();
3055 
3056  float presdphi = track->get_cal_dphi(SvtxTrack::PRES);
3057  float presdeta = track->get_cal_deta(SvtxTrack::PRES);
3058  float prese3x3 = track->get_cal_energy_3x3(SvtxTrack::PRES);
3059  float prese = track->get_cal_cluster_e(SvtxTrack::PRES);
3060 
3061  float cemcdphi = track->get_cal_dphi(SvtxTrack::CEMC);
3062  float cemcdeta = track->get_cal_deta(SvtxTrack::CEMC);
3063  float cemce3x3 = track->get_cal_energy_3x3(SvtxTrack::CEMC);
3064  float cemce = track->get_cal_cluster_e(SvtxTrack::CEMC);
3065 
3066  float hcalindphi = track->get_cal_dphi(SvtxTrack::HCALIN);
3067  float hcalindeta = track->get_cal_deta(SvtxTrack::HCALIN);
3068  float hcaline3x3 = track->get_cal_energy_3x3(SvtxTrack::HCALIN);
3069  float hcaline = track->get_cal_cluster_e(SvtxTrack::HCALIN);
3070 
3071  float hcaloutdphi = track->get_cal_dphi(SvtxTrack::HCALOUT);
3072  float hcaloutdeta = track->get_cal_deta(SvtxTrack::HCALOUT);
3073  float hcaloute3x3 = track->get_cal_energy_3x3(SvtxTrack::HCALOUT);
3074  float hcaloute = track->get_cal_cluster_e(SvtxTrack::HCALOUT);
3075 
3076  float gtrackID = NAN;
3077  float gflavor = NAN;
3078  float ng4hits = NAN;
3079  unsigned int ngmaps = 0;
3080  unsigned int ngintt = 0;
3081  unsigned int ngmms = 0;
3082  unsigned int ngtpc = 0;
3083  unsigned int nglmaps = 0;
3084  unsigned int nglintt = 0;
3085  unsigned int nglmms = 0;
3086  unsigned int ngltpc = 0;
3087  float gpx = NAN;
3088  float gpy = NAN;
3089  float gpt = NAN;
3090  float geta = NAN;
3091  float gphi = NAN;
3092  float gpz = NAN;
3093  float gvx = NAN;
3094  float gvy = NAN;
3095  float gvz = NAN;
3096  float gvt = NAN;
3097  float gfpx = NAN;
3098  float gfpy = NAN;
3099  float gfpz = NAN;
3100  float gfx = NAN;
3101  float gfy = NAN;
3102  float gfz = NAN;
3103  float gembed = NAN;
3104  float gprimary = NAN;
3105 
3106  float nfromtruth = NAN;
3107  float nwrong = NAN;
3108  float ntrumaps = NAN;
3109  float ntruintt = NAN;
3110  float ntrumms = NAN;
3111  float ntrutpc = NAN;
3112  float ntrutpc1 = NAN;
3113  float ntrutpc11 = NAN;
3114  float ntrutpc2 = NAN;
3115  float ntrutpc3 = NAN;
3116  float layersfromtruth = NAN;
3117 
3118  if (_do_track_match)
3119  {
3120  PHG4Particle* g4particle = trackeval->max_truth_particle_by_nclusters(track);
3121  if (g4particle)
3122  {
3123  if (_scan_for_embedded)
3124  {
3125  if (trutheval->get_embed(g4particle) <= 0) continue;
3126  }
3127 
3128  gtrackID = g4particle->get_track_id();
3129  gflavor = g4particle->get_pid();
3130 
3131  std::set<TrkrDefs::cluskey> g4clusters = clustereval->all_clusters_from(g4particle);
3132  ng4hits = g4clusters.size();
3133  gpx = g4particle->get_px();
3134  gpy = g4particle->get_py();
3135  gpz = g4particle->get_pz();
3136 
3137  int lmaps[_nlayers_maps + 1];
3138  if (_nlayers_maps > 0)
3139  for (unsigned int i = 0; i < _nlayers_maps; i++) lmaps[i] = 0;
3140 
3141  int lintt[_nlayers_intt + 1];
3142  if (_nlayers_intt > 0)
3143  for (unsigned int i = 0; i < _nlayers_intt; i++) lintt[i] = 0;
3144 
3145  int ltpc[_nlayers_tpc + 1];
3146  if (_nlayers_tpc > 0)
3147  for (unsigned int i = 0; i < _nlayers_tpc; i++) ltpc[i] = 0;
3148 
3149  int lmms[_nlayers_mms + 1];
3150  if (_nlayers_mms > 0)
3151  for (unsigned int i = 0; i < _nlayers_mms; i++) lmms[i] = 0;
3152 
3153  for (const TrkrDefs::cluskey g4cluster : g4clusters)
3154  {
3155  unsigned int layer = TrkrDefs::getLayer(g4cluster);
3156  if (_nlayers_maps > 0 && layer < _nlayers_maps)
3157  {
3158  lmaps[layer] = 1;
3159  ngmaps++;
3160  }
3161 
3162  if (_nlayers_intt > 0 && layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt)
3163  {
3164  lintt[layer - _nlayers_maps] = 1;
3165  ngintt++;
3166  }
3167 
3168  if (_nlayers_tpc > 0 && layer >= _nlayers_maps + _nlayers_intt && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
3169  {
3170  ltpc[layer - (_nlayers_maps + _nlayers_intt)] = 1;
3171  ngtpc++;
3172  }
3173 
3174  if (_nlayers_mms > 0 && layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
3175  {
3176  lmms[layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3177  ngmms++;
3178  }
3179 
3180  }
3181  if (_nlayers_maps > 0)
3182  for (unsigned int i = 0; i < _nlayers_maps; i++) nglmaps += lmaps[i];
3183  if (_nlayers_intt > 0)
3184  for (unsigned int i = 0; i < _nlayers_intt; i++) nglintt += lintt[i];
3185  if (_nlayers_tpc > 0)
3186  for (unsigned int i = 0; i < _nlayers_tpc; i++) ngltpc += ltpc[i];
3187  if (_nlayers_mms > 0)
3188  for (unsigned int i = 0; i < _nlayers_mms; i++) nglmms += lmms[i];
3189 
3190  TVector3 gv(gpx, gpy, gpz);
3191  gpt = gv.Pt();
3192  geta = gv.Eta();
3193  gphi = gv.Phi();
3194  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
3195  gvx = vtx->get_x();
3196  gvy = vtx->get_y();
3197  gvz = vtx->get_z();
3198  gvt = vtx->get_t();
3199 
3200  PHG4Hit* outerhit = nullptr;
3201  if (_do_eval_light == false)
3202  outerhit = trutheval->get_outermost_truth_hit(g4particle);
3203  if (outerhit)
3204  {
3205  gfpx = outerhit->get_px(1);
3206  gfpy = outerhit->get_py(1);
3207  gfpz = outerhit->get_pz(1);
3208  gfx = outerhit->get_x(1);
3209  gfy = outerhit->get_y(1);
3210  gfz = outerhit->get_z(1);
3211  }
3212  gembed = trutheval->get_embed(g4particle);
3213  gprimary = trutheval->is_primary(g4particle);
3214 
3215  nfromtruth = trackeval->get_nclusters_contribution(track, g4particle);
3216  nwrong = trackeval->get_nwrongclusters_contribution(track, g4particle);
3217  if (_nlayers_maps == 0)
3218  {
3219  ntrumaps = 0;
3220  }
3221  else
3222  {
3223  ntrumaps = trackeval->get_layer_range_contribution(track, g4particle, 0, _nlayers_maps);
3224  }
3225  if (_nlayers_intt == 0)
3226  {
3227  ntruintt = 0;
3228  }
3229  else
3230  {
3231  ntruintt = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps, _nlayers_maps + _nlayers_intt);
3232  }
3233  if (_nlayers_mms == 0)
3234  {
3235  ntrumms = 0;
3236  }
3237  else
3238  {
3239  ntrumms = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + _nlayers_tpc, _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms);
3240  }
3241  ntrutpc = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
3242  ntrutpc1 = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 16);
3243  ntrutpc11 = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 8);
3244  ntrutpc2 = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt+16, _nlayers_maps + _nlayers_intt + 32);
3245  ntrutpc3 = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt+32, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
3246  layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track, g4particle);
3247  }
3248  }
3249 
3250  float track_data[] = {(float) _ievent,m_fSeed,
3251  trackID,
3252  px,
3253  py,
3254  pz,
3255  pt,
3256  eta,
3257  phi,
3258  deltapt,
3259  deltaeta,
3260  deltaphi,
3261  charge,
3262  quality,
3263  chisq,
3264  ndf,
3265  nhits, nmaps, nintt, ntpc,nmms,
3266  ntpc1,ntpc11,ntpc2,ntpc3,
3267  nlmaps, nlintt, nltpc,nlmms,
3268  (float) layers,
3269  dca2d,
3270  dca2dsigma,
3271  dca3dxy,
3272  dca3dxysigma,
3273  dca3dz,
3274  dca3dzsigma,
3275  pcax,
3276  pcay,
3277  pcaz,
3278  presdphi,
3279  presdeta,
3280  prese3x3,
3281  prese,
3282  cemcdphi,
3283  cemcdeta,
3284  cemce3x3,
3285  cemce,
3286  hcalindphi,
3287  hcalindeta,
3288  hcaline3x3,
3289  hcaline,
3290  hcaloutdphi,
3291  hcaloutdeta,
3292  hcaloute3x3,
3293  hcaloute,
3294  gtrackID,
3295  gflavor,
3296  ng4hits,
3297  (float) ngmaps,
3298  (float) ngintt,
3299  (float) ngtpc,
3300  (float) ngmms,
3301  (float) nglmaps,
3302  (float) nglintt,
3303  (float) ngltpc,
3304  (float) nglmms,
3305  gpx,
3306  gpy,
3307  gpz,
3308  gpt,
3309  geta,
3310  gphi,
3311  gvx,
3312  gvy,
3313  gvz,
3314  gvt,
3315  gfpx,
3316  gfpy,
3317  gfpz,
3318  gfx,
3319  gfy,
3320  gfz,
3321  gembed,
3322  gprimary,
3323  nfromtruth,
3324  nwrong,
3325  ntrumaps,
3326  ntruintt,
3327  ntrutpc,
3328  ntrumms,
3329  ntrutpc1,
3330  ntrutpc11,
3331  ntrutpc2,
3332  ntrutpc3,
3333  layersfromtruth,
3334  nhit_tpc_all,
3335  nhit_tpc_in,
3336  nhit_tpc_mid,
3337  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
3338 
3339  if(Verbosity() >= 1)
3340  cout << "ievent " << _ievent
3341  << " trackID " << trackID
3342  << " nhits " << nhits
3343  << " px " << px
3344  << " py " << py
3345  << " pz " << pz
3346  << " gembed " << gembed
3347  << " gprimary " << gprimary
3348  << endl;
3349 
3350  _ntp_track->Fill(track_data);
3351  }
3352  }
3353  if (Verbosity() > 1)
3354  {
3355  _timer->stop();
3356  cout << "track time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
3357  }
3358  }
3359 
3360  //---------------------
3361  // fill the Gseed NTuple
3362  //---------------------
3363 
3364  if (_ntp_gseed)
3365  {
3366  if (Verbosity() > 1)
3367  {
3368  cout << "Filling ntp_gseed " << endl;
3369  _timer->restart();
3370  }
3371 
3372  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
3373 
3374  float gx = NAN;
3375  float gy = NAN;
3376  float gz = NAN;
3377  float gr = NAN;
3378  float geta = NAN;
3379  float gphi = NAN;
3380  float glayer = NAN;
3381  float gpx = NAN;
3382  float gpy = NAN;
3383  float gpz = NAN;
3384  float gtpt = NAN;
3385  float gtphi = NAN;
3386  float gteta = NAN;
3387  float gvx = NAN;
3388  float gvy = NAN;
3389  float gvz = NAN;
3390  float gembed = NAN;
3391  float gprimary = NAN;
3392  float gflav = NAN;
3393  float dphiprev = NAN;
3394  float detaprev = NAN;
3395 
3396  float xval[_nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms];
3397  float yval[_nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms];
3398  float zval[_nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms];
3399  if (truthinfo)
3400  {
3401  int ntrk = 0;
3402  PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
3403  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
3404  iter != range.second;
3405  ++iter)
3406  {
3407  ntrk++;
3408  PHG4Particle* g4particle = iter->second;
3409  for (unsigned int i = 0; i < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms; i++)
3410  {
3411  xval[i] = 0;
3412  yval[i] = 0;
3413  zval[i] = 0;
3414  }
3415  std::set<PHG4Hit*> truth_hits = trutheval->all_truth_hits(g4particle);
3416  for (std::set<PHG4Hit*>::iterator iter = truth_hits.begin();
3417  iter != truth_hits.end();
3418  ++iter)
3419  {
3420  PHG4Hit* g4hit = *iter;
3421  unsigned int layer = g4hit->get_layer();
3422  //cout << " g4hit " << g4hit->get_hit_id() << " layer = " << layer << endl;
3423  if (layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
3424  {
3425  //cout << PHWHERE << " skipping out of bounds detector id " << layer << endl;
3426  continue;
3427  }
3428  xval[layer] = g4hit->get_avg_x();
3429  yval[layer] = g4hit->get_avg_y();
3430  zval[layer] = g4hit->get_avg_z();
3431  }
3432 
3433  for (unsigned int i = 0; i < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms; i++)
3434  {
3435  gx = xval[i];
3436  gy = yval[i];
3437  gz = zval[i];
3438  if (gx == 0 && gy == 0) continue;
3439 
3440  TVector3 vg4(gx, gy, gz);
3441  glayer = i;
3442  gr = vg4.Perp();
3443  geta = vg4.Eta();
3444  gphi = vg4.Phi();
3445  gpx = g4particle->get_px();
3446  gpy = g4particle->get_py();
3447  gpz = g4particle->get_pz();
3448  TVector3 vg4p(gpx, gpy, gpz);
3449 
3450  gtpt = vg4p.Perp();
3451  gtphi = vg4p.Phi();
3452  gteta = vg4p.Eta();
3453 
3454  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
3455 
3456  if (vtx)
3457  {
3458  gvx = vtx->get_x();
3459  gvy = vtx->get_y();
3460  gvz = vtx->get_z();
3461  }
3462 
3463  gembed = trutheval->get_embed(g4particle);
3464  gprimary = trutheval->is_primary(g4particle);
3465  gflav = g4particle->get_pid();
3466  if (i >= 1)
3467  {
3468  if (xval[i - 1] != 0 && yval[i - 1] != 0)
3469  {
3470  TVector3 vg4prev(xval[i - 1], yval[i - 1], zval[i - 1]);
3471  dphiprev = vg4.DeltaPhi(vg4prev);
3472  detaprev = geta - vg4prev.Eta();
3473  }
3474  }
3475 
3476  float ntrk_f = ntrk;
3477  float _ievent_f = _ievent;
3478  float gseed_data[] = {_ievent_f,
3479  ntrk_f,
3480  gx,
3481  gy,
3482  gz,
3483  gr,
3484  geta,
3485  gphi,
3486  glayer,
3487  gpx,
3488  gpy,
3489  gpz,
3490  gtpt,
3491  gtphi,
3492  gteta,
3493  gvx,
3494  gvy,
3495  gvz,
3496  gembed,
3497  gprimary,
3498  gflav,
3499  dphiprev,
3500  detaprev,
3501  nhit_tpc_all,
3502  nhit_tpc_in,
3503  nhit_tpc_mid,
3504  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
3505 
3506  _ntp_gseed->Fill(gseed_data);
3507  }
3508  }
3509  }
3510 
3511  if (Verbosity() > 1)
3512  {
3513  _timer->stop();
3514  cout << "g4hit time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
3515  }
3516  }
3517  return;
3518 
3519 }
3520 
3521 TMatrixF SvtxEvaluator::calculateClusterError(TrkrCluster* c, float& clusphi)
3522 {
3523  TMatrixF localErr(3,3);
3524  localErr[0][0] = 0.;
3525  localErr[0][1] = 0.;
3526  localErr[0][2] = 0.;
3527  localErr[1][0] = 0.;
3528  localErr[1][1] = c->getActsLocalError(0,0);
3529  localErr[1][2] = c->getActsLocalError(0,1);
3530  localErr[2][0] = 0.;
3531  localErr[2][1] = c->getActsLocalError(1,0);
3532  localErr[2][2] = c->getActsLocalError(1,1);
3533 
3534  TMatrixF ROT(3,3);
3535  ROT[0][0] = cos(clusphi);
3536  ROT[0][1] = -sin(clusphi);
3537  ROT[0][2] = 0.0;
3538  ROT[1][0] = sin(clusphi);
3539  ROT[1][1] = cos(clusphi);
3540  ROT[1][2] = 0.0;
3541  ROT[2][0] = 0.0;
3542  ROT[2][1] = 0.0;
3543  ROT[2][2] = 1.0;
3544  TMatrixF ROT_T(3,3);
3545  ROT_T.Transpose(ROT);
3546 
3547  TMatrixF err(3,3);
3548  err = ROT * localErr * ROT_T;
3549  return err;
3550 }