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