Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrkrNtuplizer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrkrNtuplizer.cc
1 #include "TrkrNtuplizer.h"
2 
3 #include <trackbase/ActsGeometry.h>
4 #include <trackbase/ClusterErrorPara.h>
5 #include <trackbase/TrkrCluster.h>
6 #include <trackbase/TrkrDefs.h>
7 #include <trackbase/TrkrHit.h>
8 #include <trackbase/TrkrHitSetContainer.h>
9 #include <trackbase/TrackFitUtils.h>
10 
11 #include <trackbase_historic/TrackSeedContainer_v1.h>
12 
13 #include <trackbase/TpcDefs.h>
14 
15 #include <trackbase/TrkrClusterContainer.h>
16 #include <trackbase/TrkrClusterHitAssoc.h>
17 #include <trackbase/TrkrClusterIterationMapv1.h>
18 #include <trackbase/TrkrHitSet.h>
19 
20 #include <globalvertex/GlobalVertex.h>
21 #include <globalvertex/GlobalVertexMap.h>
22 #include <globalvertex/SvtxVertex.h>
23 #include <globalvertex/SvtxVertexMap.h>
24 
25 #include <trackbase_historic/TrackAnalysisUtils.h>
26 #include <trackbase_historic/ActsTransformations.h>
27 #include <trackbase_historic/SvtxTrack.h>
28 #include <trackbase_historic/SvtxTrackMap.h>
29 #include <trackbase_historic/TrackSeed.h>
30 
31 #include <g4detectors/PHG4TpcCylinderGeom.h>
32 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
33 #include <trackermillepedealignment/HelicalFitter.h>
34 
35 #include <fun4all/Fun4AllReturnCodes.h>
36 #include <fun4all/SubsysReco.h>
37 
38 #include <phool/PHTimer.h>
39 #include <phool/getClass.h>
40 #include <phool/phool.h>
41 #include <phool/recoConsts.h>
42 
43 #include <TFile.h>
44 #include <TNtuple.h>
45 #include <TVector3.h>
46 
47 #include <cmath>
48 #include <iomanip>
49 #include <iostream>
50 #include <iterator>
51 #include <map>
52 #include <memory> // for shared_ptr
53 #include <set> // for _Rb_tree_cons...
54 #include <utility>
55 #include <vector>
56 
57 using namespace std;
58 
59 TrkrNtuplizer::TrkrNtuplizer(const string& /*name*/, const string& filename, const string& trackmapname,
60  unsigned int nlayers_maps,
61  unsigned int nlayers_intt,
62  unsigned int nlayers_tpc,
63  unsigned int nlayers_mms)
64  : SubsysReco("TrkrNtuplizer")
65  , _ievent(0)
66  , _iseed(0)
67  , m_fSeed(NAN)
68  , _do_info_eval(true)
69  , _do_vertex_eval(true)
70  , _do_hit_eval(true)
71  , _do_cluster_eval(true)
72  , _do_clus_trk_eval(true)
73  , _do_track_eval(true)
74  , _do_tpcseed_eval(false)
75  , _do_siseed_eval(false)
76  , _nlayers_maps(nlayers_maps)
77  , _nlayers_intt(nlayers_intt)
78  , _nlayers_tpc(nlayers_tpc)
79  , _nlayers_mms(nlayers_mms)
80  , _ntp_info(nullptr)
81  , _ntp_vertex(nullptr)
82  , _ntp_hit(nullptr)
83  , _ntp_cluster(nullptr)
84  , _ntp_clus_trk(nullptr)
85  , _ntp_track(nullptr)
86  , _ntp_tpcseed(nullptr)
87  , _ntp_siseed(nullptr)
88  , _filename(filename)
89  , _trackmapname(trackmapname)
90  , _tfile(nullptr)
91  , _timer(nullptr)
92 {
93 }
94 
95 TrkrNtuplizer::~TrkrNtuplizer()
96 {
97  delete _timer;
98 }
99 
100 int TrkrNtuplizer::Init(PHCompositeNode* /*topNode*/)
101 {
102  _ievent = 0;
103 
104  _tfile = new TFile(_filename.c_str(), "RECREATE");
105  _tfile->SetCompressionLevel(7);
106  if (_do_info_eval)
107  {
108  _ntp_info = new TNtuple("ntp_info", "event info",
109  "event:seed:"
110  "occ11:occ116:occ21:occ216:occ31:occ316:"
111  "ntrk:"
112  "nhitmvtx:nhitintt:nhittpot:"
113  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
114  }
115 
116  if (_do_vertex_eval)
117  {
118  _ntp_vertex = new TNtuple("ntp_vertex", "vertex => max truth",
119  "event:seed:vertexID:vx:vy:vz:ntracks:chi2:ndof:"
120  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
121  }
122 
123  if (_do_hit_eval)
124  {
125  _ntp_hit = new TNtuple("ntp_hit", "svtxhit => max truth",
126  "event:seed:hitID:e:adc:layer:phielem:zelem:"
127  "cellID:ecell:phibin:tbin:phi:x:y:z:"
128  "nhitmvtx:nhitintt:nhittpot:"
129  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
130  }
131 
132  if (_do_cluster_eval)
133  {
134  _ntp_cluster = new TNtuple("ntp_cluster", "svtxcluster => max truth",
135  "event:seed:hitID:locx:locy:x:y:z:r:phi:eta:theta:phibin:tbin:ex:ey:ez:ephi:pez:pephi:"
136  "e:adc:maxadc:layer:phielem:zelem:size:phisize:zsize:"
137  "trackID:niter:"
138  "nhitmvtx:nhitintt:nhittpot:"
139  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
140  }
141  if (_do_clus_trk_eval)
142  {
143  _ntp_clus_trk = new TNtuple("ntp_clus_trk", "cluster on track",
144  "event:seed:locx:locy:x:y:z:r:phi:eta:theta:phibin:tbin:ex:ey:ez:ephi:pez:pephi:"
145  "e:adc:maxadc:layer:phielem:zelem:size:phisize:zsize:"
146  "resphi:resz:"
147  "trackID:niter:pt:eta:phi:X0:Y0:Z0:charge:"
148  "nhitmvtx:nhitintt:nhittpot:"
149  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
150  }
151 
152  if (_do_track_eval)
153  {
154  _ntp_track = new TNtuple("ntp_track", "svtxtrack => max truth",
155  "event:seed:trackID:crossing:px:py:pz:pt:eta:phi:deltapt:deltaeta:deltaphi:charge:"
156  "quality:chisq:ndf:nhits:nmaps:nintt:ntpc:nmms:ntpc1:ntpc11:ntpc2:ntpc3:nlmaps:nlintt:nltpc:nlmms:layers:"
157  "vertexID:vx:vy:vz:dca2d:dca2dsigma:dca3dxy:dca3dxysigma:dca3dz:dca3dzsigma:pcax:pcay:pcaz:"
158  "nhitmvtx:nhitintt:nhittpot:"
159  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
160  }
161 
162  if (_do_tpcseed_eval)
163  {
164  _ntp_tpcseed = new TNtuple("ntp_tpcseed", "seeds from truth",
165  "event:seed:ntrk:gx:gy:gz:gr:geta:gphi:"
166  "glayer:"
167  "gpx:gpy:gpz:gtpt:gtphi:gteta:"
168  "gvx:gvy:gvz:"
169  "gembed:gprimary:gflav:"
170  "dphiprev:detaprev:"
171  "nhitmvtx:nhitintt:nhittpot:"
172  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
173  }
174  if (_do_siseed_eval)
175  {
176  _ntp_tpcseed = new TNtuple("ntp_siseed", "seeds from truth",
177  "event:seed:ntrk:gx:gy:gz:gr:geta:gphi:"
178  "glayer:"
179  "gpx:gpy:gpz:gtpt:gtphi:gteta:"
180  "gvx:gvy:gvz:"
181  "gembed:gprimary:gflav:"
182  "dphiprev:detaprev:"
183  "nhitmvtx:nhitintt:nhittpot:"
184  "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
185  }
186  _timer = new PHTimer("_eval_timer");
187  _timer->stop();
188 
189  return Fun4AllReturnCodes::EVENT_OK;
190 }
191 
192 int TrkrNtuplizer::InitRun(PHCompositeNode*)
193 {
194  return Fun4AllReturnCodes::EVENT_OK;
195 }
196 
197 int TrkrNtuplizer::process_event(PHCompositeNode* topNode)
198 {
199  if ((Verbosity() > 1) && (_ievent % 100 == 0))
200  {
201  cout << "TrkrNtuplizer::process_event - Event = " << _ievent << endl;
202  }
203 
204  recoConsts* rc = recoConsts::instance();
205  if (rc->FlagExist("RANDOMSEED"))
206  {
207  _iseed = rc->get_IntFlag("RANDOMSEED");
208  m_fSeed = _iseed;
209  }
210  else
211  {
212  _iseed = 0;
213  m_fSeed = NAN;
214  }
215  if(_trackmap == nullptr)
216  _trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname.c_str());
217 
218  if (Verbosity() > 1)
219  {
220  cout << "TrkrNtuplizer::process_event - Seed = " << _iseed << endl;
221  }
222  //-----------------------------------
223  // print what is coming into the code
224  //-----------------------------------
225 
226  printInputInfo(topNode);
227 
228  //---------------------------
229  // fill the Evaluator NTuples
230  //---------------------------
231 
232  fillOutputNtuples(topNode);
233 
234  //--------------------------------------------------
235  // Print out the ancestry information for this event
236  //--------------------------------------------------
237 
238  // printOutputInfo(topNode);
239 
240  ++_ievent;
241  return Fun4AllReturnCodes::EVENT_OK;
242 }
243 
244 int TrkrNtuplizer::End(PHCompositeNode* /*topNode*/)
245 {
246  _tfile->cd();
247 
248  if (_ntp_info)
249  {
250  _ntp_info->Write();
251  }
252  if (_ntp_vertex)
253  {
254  _ntp_vertex->Write();
255  }
256  if (_ntp_hit)
257  {
258  _ntp_hit->Write();
259  }
260  if (_ntp_cluster)
261  {
262  _ntp_cluster->Write();
263  }
264  if (_ntp_clus_trk)
265  {
266  _ntp_clus_trk->Write();
267  }
268  if (_ntp_track)
269  {
270  _ntp_track->Write();
271  }
272  if (_ntp_tpcseed)
273  {
274  _ntp_tpcseed->Write();
275  }
276  if (_ntp_siseed)
277  {
278  _ntp_siseed->Write();
279  }
280 
281  _tfile->Close();
282 
283  delete _tfile;
284 
285  if (Verbosity() > 1)
286  {
287  cout << "========================= TrkrNtuplizer::End() ============================" << endl;
288  cout << " " << _ievent << " events of output written to: " << _filename << endl;
289  cout << "===========================================================================" << endl;
290  }
291 
292  return Fun4AllReturnCodes::EVENT_OK;
293 }
294 
295 void TrkrNtuplizer::printInputInfo(PHCompositeNode* topNode)
296 {
297  if (Verbosity() > 1)
298  {
299  cout << "TrkrNtuplizer::printInputInfo() entered" << endl;
300  }
301 
302  if (Verbosity() > 3)
303  {
304  // event information
305  cout << endl;
306  cout << PHWHERE << " INPUT FOR EVENT " << _ievent << endl;
307 
308  cout << endl;
309 
310  cout << "---SVTXCLUSTERS-------------" << endl;
311  TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
312  if (!clustermap)
313  {
314  clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
315  }
316 
317  if (clustermap != nullptr)
318  {
319  unsigned int icluster = 0;
320  for (const auto& hitsetkey : clustermap->getHitSetKeys())
321  {
322  auto range = clustermap->getClusters(hitsetkey);
323  for (auto iter = range.first; iter != range.second; ++iter)
324  {
325  TrkrDefs::cluskey cluster_key = iter->first;
326  cout << icluster << " with key " << cluster_key << " of " << clustermap->size();
327  cout << ": SvtxCluster: " << endl;
328  iter->second->identify();
329  ++icluster;
330  }
331  }
332  }
333 
334  cout << "---SVXTRACKS-------------" << endl;
335 
336  if (_trackmap)
337  {
338  unsigned int itrack = 0;
339  for (SvtxTrackMap::Iter iter = _trackmap->begin();
340  iter != _trackmap->end();
341  ++iter)
342  {
343  cout << itrack << " of " << _trackmap->size();
344  SvtxTrack* track = iter->second;
345  cout << " : SvtxTrack:" << endl;
346  track->identify();
347  cout << endl;
348  ++itrack;
349  }
350  }
351 
352  cout << "---SVXVERTEXES-------------" << endl;
353  SvtxVertexMap* vertexmap = nullptr;
354  vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs"); // Acts vertices
355 
356  if (vertexmap)
357  {
358  unsigned int ivertex = 0;
359  for (SvtxVertexMap::Iter iter = vertexmap->begin();
360  iter != vertexmap->end();
361  ++iter)
362  {
363  cout << ivertex << " of " << vertexmap->size();
364  SvtxVertex* vertex = iter->second;
365  cout << " : SvtxVertex:" << endl;
366  vertex->identify();
367  cout << endl;
368  }
369  }
370  }
371 
372  return;
373 }
374 
375 void TrkrNtuplizer::printOutputInfo(PHCompositeNode* topNode)
376 {
377  if (Verbosity() > 1)
378  {
379  cout << "TrkrNtuplizer::printOutputInfo() entered" << endl;
380  }
381 
382  //==========================================
383  // print out some useful stuff for debugging
384  //==========================================
385 
386  if (Verbosity() > 100)
387  {
388  // event information
389  cout << endl;
390  cout << PHWHERE << " NEW OUTPUT FOR EVENT " << _ievent << endl;
391  cout << endl;
392 
393  float vx = NAN;
394  float vy = NAN;
395  float vz = NAN;
396 
397  SvtxVertexMap* vertexmap = nullptr;
398 
399  vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs"); // Acts vertices
400 
401  if (vertexmap)
402  {
403  if (!vertexmap->empty())
404  {
405  SvtxVertex* vertex = (vertexmap->begin()->second);
406 
407  vx = vertex->get_x();
408  vy = vertex->get_y();
409  vz = vertex->get_z();
410  }
411  }
412 
413  cout << "===Vertex Reconstruction=======================" << endl;
414  cout << "vreco = (" << vx << "," << vy << "," << vz << ")" << endl;
415  cout << endl;
416 
417  cout << "===Tracking Summary============================" << endl;
418 
419  TrkrHitSetContainer* hitsetmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
420 
421  TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
422  if (!clustermap)
423  {
424  clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
425  }
426 
427  unsigned int nclusters[100] = {0};
428  unsigned int nhits[100] = {0};
429 
430  ActsGeometry* tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
431 
432  if (!tgeometry)
433  {
434  std::cout << PHWHERE << "No Acts geometry on node tree. Can't continue."
435  << std::endl;
436  }
437 
438  if (hitsetmap)
439  {
440  TrkrHitSetContainer::ConstRange all_hitsets = hitsetmap->getHitSets();
441  for (TrkrHitSetContainer::ConstIterator hitsetiter = all_hitsets.first;
442  hitsetiter != all_hitsets.second;
443  ++hitsetiter)
444  {
445  // we have a single hitset, get the layer
446  unsigned int layer = TrkrDefs::getLayer(hitsetiter->first);
447 
448  // count all hits in this hitset
449  TrkrHitSet::ConstRange hitrangei = hitsetiter->second->getHits();
450  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
451  hitr != hitrangei.second;
452  ++hitr)
453  {
454  ++nhits[layer];
455  }
456  auto range = clustermap->getClusters(hitsetiter->first);
457  for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
458  {
459  const auto cluskey = clusIter->first;
460  // const auto cluster = clusIter->second;
461  // unsigned int layer = TrkrDefs::getLayer(cluskey);
462  nclusters[TrkrDefs::getLayer(cluskey)]++;
463  }
464  }
465  }
466 
467  PHG4TpcCylinderGeomContainer* geom_container =
468  findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
469  {
470  if (!geom_container)
471  {
472  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
473  }
474  return;
475  }
476 
477  for (unsigned int ilayer = 0; ilayer < _nlayers_maps + _nlayers_intt + _nlayers_tpc; ++ilayer)
478  {
479  PHG4TpcCylinderGeom* GeoLayer = geom_container->GetLayerCellGeom(ilayer);
480 
481  cout << "layer " << ilayer
482  << " => nHits = " << nhits[ilayer]
483  << " => nClusters = " << nclusters[ilayer] << endl;
484  if (ilayer >= _nlayers_maps + _nlayers_intt && ilayer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
485  {
486  cout << "layer " << ilayer
487  << " => nphi = " << GeoLayer->get_phibins()
488  << " => nz = " << GeoLayer->get_zbins()
489  << " => ntot = " << GeoLayer->get_phibins() * GeoLayer->get_zbins()
490  << endl;
491  }
492  }
493 
494  cout << " => nTracks = ";
495  if (_trackmap)
496  {
497  cout << _trackmap->size() << endl;
498  }
499  else
500  {
501  cout << 0 << endl;
502  }
503 
504  cout << endl;
505 
506  } // if Verbosity()
507 
508  return;
509 }
510 
511 void TrkrNtuplizer::fillOutputNtuples(PHCompositeNode* topNode)
512 {
513  if (Verbosity() > 1)
514  {
515  cout << "TrkrNtuplizer::fillOutputNtuples() entered" << endl;
516  }
517 
518  ActsGeometry* tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
519  if (!tgeometry)
520  {
521  std::cout << "No Acts geometry on node tree. Can't continue."
522  << std::endl;
523  return;
524  }
525 
526  float nhit_tpc_all = 0;
527  float nhit_tpc_in = 0;
528  float nhit_tpc_mid = 0;
529  float nhit_tpc_out = 0;
530 
531  float nhit_intt = 0;
532  float nhit_maps = 0;
533  float nhit_mms = 0;
534 
535  float nclus_all = 0;
536  float nclus_tpc = 0;
537  float nclus_intt = 0;
538  float nclus_maps = 0;
539  float nclus_mms = 0;
540  float nhit[100];
541  for (float& i : nhit) i = 0;
542  float occ11 = 0;
543  float occ116 = 0;
544  float occ21 = 0;
545  float occ216 = 0;
546  float occ31 = 0;
547  float occ316 = 0;
548 
549  TrkrHitSetContainer* hitmap_in = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
550  if (hitmap_in)
551  {
552  TrkrHitSetContainer::ConstRange all_hitsets = hitmap_in->getHitSets();
553  for (TrkrHitSetContainer::ConstIterator hitsetiter = all_hitsets.first;
554  hitsetiter != all_hitsets.second;
555  ++hitsetiter)
556  {
557  // we have a single hitset, get the layer
558  unsigned int layer = TrkrDefs::getLayer(hitsetiter->first);
559 
560  TrkrHitSet::ConstRange hitrangei = hitsetiter->second->getHits();
561  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
562  hitr != hitrangei.second;
563  ++hitr)
564  {
565  nhit[layer]++;
566 
567  if(layer<_nlayers_maps){
568  nhit_maps++;
569  }
570  if(layer>=_nlayers_maps&&layer<_nlayers_maps + _nlayers_intt){
571  nhit_intt++;
572  }
573  if ((float) layer >= _nlayers_maps + _nlayers_intt &&
574  (float) layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
575  nhit_tpc_all++;
576  if((float) layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
577  nhit_mms++;
578  if ((float) layer == _nlayers_maps + _nlayers_intt)
579  {
580  nhit_tpc_in++;
581  }
582  if ((float) layer == _nlayers_maps + _nlayers_intt + _nlayers_tpc - 1)
583  {
584  nhit_tpc_out++;
585  }
586  if ((float) layer == _nlayers_maps + _nlayers_intt + _nlayers_tpc / 2 - 1)
587  {
588  nhit_tpc_mid++;
589  }
590  }
591  }
592  }
593  /**********/
594 
595  PHG4TpcCylinderGeomContainer* geom_container =
596  findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
597  if (!geom_container)
598  {
599  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
600  return;
601  }
602  float nhits = 0;
603  /*
604  PHG4TpcCylinderGeom* GeoLayer;
605  int layer = _nlayers_maps + _nlayers_intt;
606  GeoLayer = geom_container->GetLayerCellGeom(layer);
607  int nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
608  float nhits = nhit[layer];
609  occ11 = nhits / nbins;
610  if (Verbosity() > 1)
611  {
612  cout << " occ11 = " << occ11
613  << " nbins = " << nbins
614  << " nhits = " << nhits
615  << " layer = " << layer
616  << endl;
617  }
618  layer = _nlayers_maps + _nlayers_intt + 15;
619  GeoLayer = geom_container->GetLayerCellGeom(layer);
620  nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
621  nhits = nhit[layer];
622  occ116 = nhits / nbins;
623 
624  layer = _nlayers_maps + _nlayers_intt + 16;
625  GeoLayer = geom_container->GetLayerCellGeom(layer);
626  nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
627  nhits = nhit[layer];
628  occ21 = nhits / nbins;
629  layer = _nlayers_maps + _nlayers_intt + 31;
630  GeoLayer = geom_container->GetLayerCellGeom(layer);
631  nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
632  nhits = nhit[layer];
633  occ216 = nhits / nbins;
634  layer = _nlayers_maps + _nlayers_intt + 32;
635  GeoLayer = geom_container->GetLayerCellGeom(layer);
636  nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
637  nhits = nhit[layer];
638  occ31 = nhits / nbins;
639  layer = _nlayers_maps + _nlayers_intt + 47;
640  GeoLayer = geom_container->GetLayerCellGeom(layer);
641  nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
642  nhits = nhit[layer];
643  occ316 = nhits / nbins;
644 
645  if (Verbosity() > 1)
646  {
647  cout << " occ11 = " << occ11
648  << " occ116 = " << occ116
649  << " occ21 = " << occ21
650  << " occ216 = " << occ216
651  << " occ31 = " << occ31
652  << " occ316 = " << occ316
653  << endl;
654  }
655  */
656  TrkrClusterContainer* clustermap_in = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
657  if (!clustermap_in)
658  {
659  clustermap_in = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
660  }
661 
662  if (clustermap_in)
663  {
664  nclus_all = clustermap_in->size();
665 
666  for (const auto& hitsetkey : clustermap_in->getHitSetKeys())
667  {
668  auto range = clustermap_in->getClusters(hitsetkey);
669  for (auto iter_cin = range.first; iter_cin != range.second; ++iter_cin)
670  {
671  TrkrDefs::cluskey cluster_key = iter_cin->first;
672  unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
673  if (_nlayers_maps > 0)
674  {
675  if (layer_local < _nlayers_maps)
676  {
677  nclus_maps++;
678  }
679  }
680  if (_nlayers_intt > 0)
681  {
682  if ((layer_local >= _nlayers_maps) && (layer_local < (_nlayers_maps + _nlayers_intt)))
683  {
684  nclus_intt++;
685  }
686  }
687  if (_nlayers_tpc > 0)
688  {
689  if (layer_local >= (_nlayers_maps + _nlayers_intt) && layer_local < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
690  {
691  nclus_tpc++;
692  }
693  }
694  if (_nlayers_mms > 0)
695  {
696  if (layer_local >= (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
697  {
698  nclus_mms++;
699  }
700  }
701  }
702  }
703  }
704  //-----------------------
705  // fill the info NTuple
706  //-----------------------
707  if (_ntp_info)
708  {
709  if (Verbosity() > 1)
710  {
711  cout << "Filling ntp_info " << endl;
712  }
713  float ntrk = 0;
714 
715  if (_trackmap)
716  {
717  ntrk = (float) _trackmap->size();
718  }
719  if (Verbosity() > 0)
720  {
721  cout << "EVENTINFO SEED: " << m_fSeed << endl;
722  cout << "EVENTINFO NHIT: " << setprecision(9) << nhit_tpc_all << endl;
723  cout << "EVENTINFO CLUSTPC: " << nclus_tpc << endl;
724  cout << "EVENTINFO NTRKREC: " << ntrk << endl;
725  }
726  float info_data[] = {(float) _ievent, m_fSeed,
727  occ11, occ116, occ21, occ216, occ31, occ316,
728  ntrk,
729  nhit_maps,
730  nhit_intt,
731  nhit_mms,
732  nhit_tpc_all,
733  nhit_tpc_in,
734  nhit_tpc_mid,
735  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
736 
737  _ntp_info->Fill(info_data);
738  }
739 
740  //-----------------------
741  // fill the Vertex NTuple
742  //-----------------------
743  bool doit = true;
744  if (_ntp_vertex && doit)
745  {
746  if (Verbosity() > 1)
747  {
748  cout << "Filling ntp_vertex " << endl;
749  cout << "start vertex time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
750  _timer->restart();
751  }
752 
753  // SvtxVertexMap* vertexmap = nullptr;
754 
755  // vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs"); // Acts vertices
756 
757  float vx = NAN;
758  float vy = NAN;
759  float vz = NAN;
760  float ntracks = NAN;
761 
762  if (Verbosity() > 1)
763  {
764  std::cout << " adding vertex data " << std::endl;
765  }
766 
767  float vertex_data[] = {(float) _ievent, m_fSeed,
768  vx,
769  vy,
770  vz,
771  ntracks,
772  nhit_maps,
773  nhit_intt,
774  nhit_mms,
775  nhit_tpc_all,
776  nhit_tpc_in,
777  nhit_tpc_mid,
778  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
779 
780  _ntp_vertex->Fill(vertex_data);
781  }
782  _timer->stop();
783  cout << "vertex time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
784 
785  //--------------------
786  // fill the Hit NTuple
787  //--------------------
788 
789  if (_ntp_hit)
790  {
791  auto m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
792 
793  if (Verbosity() >= 1)
794  {
795  cout << "Filling ntp_hit " << endl;
796  _timer->restart();
797  }
798  // need things off of the DST...
799  TrkrHitSetContainer* hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
800 
801  if (hitmap)
802  {
803  TrkrHitSetContainer::ConstRange all_hitsets = hitmap->getHitSets();
804  for (TrkrHitSetContainer::ConstIterator iter = all_hitsets.first;
805  iter != all_hitsets.second;
806  ++iter)
807  {
808  const TrkrDefs::hitsetkey hitset_key = iter->first;
809  TrkrHitSet* hitset = iter->second;
810  // get all hits for this hitset
811  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
812  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
813  hitr != hitrangei.second;
814  ++hitr)
815  {
816  TrkrDefs::hitkey hit_key = hitr->first;
817  TrkrHit* hit = hitr->second;
818  float event = _ievent;
819  float hitID = hit_key;
820  float e = hit->getEnergy();
821  float adc = hit->getAdc();
822  float layer_local = TrkrDefs::getLayer(hitset_key);
823  float sector = TpcDefs::getSectorId(hitset_key);
824  float side = TpcDefs::getSide(hitset_key);
825  float cellID = 0;
826  float ecell = hit->getAdc();
827 
828  float phibin = NAN;
829  float tbin = NAN;
830  float phi = NAN;
831  float phi_center = NAN;
832  float x = NAN;
833  float y = NAN;
834  float z = NAN;
835 
836  if (TrkrDefs::getTrkrId(hitset_key) == TrkrDefs::TrkrId::tpcId)
837  {
838  PHG4TpcCylinderGeom* GeoLayer_local = geom_container->GetLayerCellGeom(layer_local);
839  double radius = GeoLayer_local->get_radius();
840  phibin = (float) TpcDefs::getPad(hit_key);
841  tbin = (float) TpcDefs::getTBin(hit_key);
842  phi = GeoLayer_local->get_phicenter(phibin);
843 
844  double zdriftlength = tbin * m_tGeometry->get_drift_velocity() * AdcClockPeriod;
845  // convert z drift length to z position in the TPC
846  // cout << " tbin: " << tbin << " vdrift " <<m_tGeometry->get_drift_velocity() << " l drift: " << zdriftlength <<endl;
847  unsigned short NTBins = (unsigned short) GeoLayer_local->get_zbins();
848  double m_tdriftmax = AdcClockPeriod * NTBins / 2.0;
849  double clusz = (m_tdriftmax * m_tGeometry->get_drift_velocity()) - zdriftlength;
850  if (side == 0)
851  {
852  clusz = -clusz;
853  }
854  z = clusz;
855  phi_center = GeoLayer_local->get_phicenter(phibin);
856  x = radius * cos(phi_center);
857  y = radius * sin(phi_center);
858  }
859 
860  float hit_data[] = {
861  event,
862  (float) _iseed,
863  hitID,
864  e,
865  adc,
866  layer_local,
867  sector,
868  side,
869  cellID,
870  ecell,
871  (float) phibin,
872  (float) tbin,
873  phi,
874  x,
875  y,
876  z,
877  nhit_maps,
878  nhit_intt,
879  nhit_mms,
880  nhit_tpc_all,
881  nhit_tpc_in,
882  nhit_tpc_mid,
883  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
884 
885  _ntp_hit->Fill(hit_data);
886  }
887  }
888  }
889  if (Verbosity() >= 1)
890  {
891  _timer->stop();
892  cout << "hit time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
893  }
894  }
895 
896  //------------------------
897  // fill the Cluster NTuple
898  //------------------------
899 
900  if (Verbosity() >= 1)
901  {
902  cout << "check for ntp_cluster" << endl;
903  _timer->restart();
904  }
905 
906  if (_ntp_cluster)
907  {
908  if (Verbosity() > 1)
909  {
910  cout << "Filling ntp_cluster (all of them) " << endl;
911  }
912  // need things off of the DST...
913  TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
914  if (!clustermap)
915  {
916  clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
917  }
918 
919  TrkrHitSetContainer* hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
920  // TrkrClusterIterationMapv1* _iteration_map = findNode::getClass<TrkrClusterIterationMapv1>(topNode, "CLUSTER_ITERATION_MAP");
921  ClusterErrorPara ClusErrPara;
922 
923  if (Verbosity() > 1)
924  {
925  if (clustermap != nullptr)
926  {
927  cout << "got clustermap" << endl;
928  }
929  else
930  {
931  cout << "no clustermap" << endl;
932  }
933 
934  if (hitsets != nullptr)
935  {
936  cout << "got hitsets" << endl;
937  }
938  else
939  {
940  cout << "no hitsets" << endl;
941  }
942  }
943 
944  if (clustermap && hitsets)
945  {
946  for (const auto& hitsetkey : clustermap->getHitSetKeys())
947  {
948  int hitsetlayer = TrkrDefs::getLayer(hitsetkey);
949  auto range = clustermap->getClusters(hitsetkey);
950  for (auto iter = range.first; iter != range.second; ++iter)
951  {
952  TrkrDefs::cluskey cluster_key = iter->first;
953  TrkrCluster* cluster = clustermap->findCluster(cluster_key);
954  SvtxTrack* track = nullptr;//best_track_from(cluster_key);
955  float niter = 0;
956 
957  float hitID = (float) cluster_key;
958  // auto trkrid = TrkrDefs::getTrkrId(cluster_key);
959  Acts::Vector3 cglob;
960  cglob = tgeometry->getGlobalPosition(cluster_key, cluster);
961  float x = cglob(0);
962  float y = cglob(1);
963  float z = cglob(2);
964  float locx = cluster->getLocalX();
965  float locy = cluster->getLocalY();
966  TVector3 pos(x, y, z);
967  float r = pos.Perp();
968  float phi = pos.Phi();
969  float eta = pos.Eta();
970  float theta = pos.Theta();
971 
972  float ex = 0;
973  float ey = 0;
974  float ez = 0;
975  float ephi = 0;
976  float pez = 0;
977  float pephi = 0;
978  float size = 0;
979  float phisize = 0;
980  float zsize = 0;
981  float maxadc = -999;
982 
983  auto para_errors = ClusErrPara.get_clusterv5_modified_error(cluster,r ,cluster_key);
984 
985  phisize = cluster->getPhiSize();
986  zsize = cluster->getZSize();
987  // double clusRadius = r;
988  ez = sqrt(para_errors.second);
989  ephi = sqrt(para_errors.first);
990  maxadc = cluster->getMaxAdc();
991 
992  float e = cluster->getAdc();
993  float adc = cluster->getAdc();
994  float layer_local = (float) TrkrDefs::getLayer(cluster_key);
995  float sector = TpcDefs::getSectorId(cluster_key);
996  float side = TpcDefs::getSide(cluster_key);
997  phisize = cluster->getPhiSize();
998  zsize = cluster->getZSize();
999  size = cluster->getSize();
1000  // double clusRadius = r;
1001  ez = sqrt(para_errors.second);
1002  ephi = sqrt(para_errors.first);
1003  maxadc = cluster->getMaxAdc();
1004  // count all hits for this cluster
1005  const TrkrDefs::hitsetkey hitsetkey_local = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
1006  int hitsetlayer2 = TrkrDefs::getLayer(hitsetkey_local);
1007  if (hitsetlayer != layer_local)
1008  {
1009  cout << "WARNING hitset layer " << hitsetlayer << "| " << hitsetlayer2 << " layer " << layer_local << endl;
1010  }
1011  /*else{
1012  cout << "Good hitset layer " << hitsetlayer << "| " << hitsetlayer2 << " layer " << layer << endl;
1013  }
1014  */
1015 
1016  float trackID = NAN;
1017  if (track != nullptr)
1018  {
1019  trackID = track->get_id();
1020  }
1021  float phibin = NAN;
1022  float tbin = NAN;
1023  if (TrkrDefs::getTrkrId(hitsetkey_local) == TrkrDefs::TrkrId::tpcId)
1024  {
1025  PHG4TpcCylinderGeom* GeoLayer_local = geom_container->GetLayerCellGeom(layer_local);
1026  phibin = GeoLayer_local->get_pad_float(phi,side);
1027  tbin = GeoLayer_local->get_tbin_float(locy-39.6);
1028  }
1029  else{
1030  phibin = locx;
1031  tbin = locy;
1032  }
1033  float cluster_data[] = {(float) _ievent,
1034  (float) _iseed,
1035  hitID,
1036  locx,
1037  locy,
1038  x,
1039  y,
1040  z,
1041  r,
1042  phi,
1043  eta,
1044  theta,
1045  phibin,
1046  tbin,
1047  ex,
1048  ey,
1049  ez,
1050  ephi,
1051  pez,
1052  pephi,
1053  e,
1054  adc,
1055  maxadc,
1056  layer_local,
1057  sector,
1058  side,
1059  size,
1060  phisize,
1061  zsize,
1062  trackID,
1063  niter,
1064  nhit_maps,
1065  nhit_intt,
1066  nhit_mms,
1067  nhit_tpc_all,
1068  nhit_tpc_in,
1069  nhit_tpc_mid,
1070  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1071 
1072  _ntp_cluster->Fill(cluster_data);
1073  }
1074  }
1075  }
1076  }
1077 
1078  if (Verbosity() >= 1)
1079  {
1080  _timer->stop();
1081  cout << "cluster time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
1082  }
1083  if (_ntp_clus_trk)
1084  {
1085  // if (Verbosity() > 1)
1086  {
1087  cout << "Filling ntp_clus_trk " << endl;
1088  }
1089  // need things off of the DST...
1090  TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
1091  if (!clustermap)
1092  {
1093  clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1094  }
1095 
1096 
1097  ClusterErrorPara ClusErrPara;
1098 
1099  if (Verbosity() > 1)
1100  {
1101  if (clustermap != nullptr) cout << "got clustermap" << endl;
1102  else cout << "no clustermap" << endl;
1103  }
1104 
1105  TrackSeedContainer *_tpc_seeds = findNode::getClass<TrackSeedContainer>(topNode, _clustrackseedcontainer);
1106  TrackSeedContainer *_tpcseeds = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
1107  TrackSeedContainer *_silseeds = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
1108  if (!_tpc_seeds)
1109  {
1110  cerr << PHWHERE << " ERROR: Can't find " << _clustrackseedcontainer << endl;
1111  return ;
1112  }
1113 
1114 
1115  for (unsigned int phtrk_iter = 0; phtrk_iter < _tpc_seeds->size(); ++phtrk_iter){
1116 
1117  TrackSeed *tpcseed = _tpc_seeds->get(phtrk_iter);
1118  if(!tpcseed) {
1119  continue;
1120  }
1121  TrackSeed *silseed = nullptr;
1122  if(_clustrackseedcontainer.find("SvtxTrackSeed") != std::string::npos)
1123  {
1124  unsigned int tpcindex = tpcseed->get_tpc_seed_index();
1125  unsigned int silindex = tpcseed->get_silicon_seed_index();
1126  tpcseed = _tpcseeds->get(tpcindex);
1127  silseed = _silseeds->get(silindex);
1128  }
1129 
1130 
1131  std::vector<Acts::Vector3> clusterPositions;
1132  std::vector<TrkrDefs::cluskey> clusterKeys;
1133  clusterKeys.insert(clusterKeys.end(), tpcseed->begin_cluster_keys(),
1134  tpcseed->end_cluster_keys());
1135  if(silseed)
1136  {
1137  clusterKeys.insert(clusterKeys.end(), silseed->begin_cluster_keys(),
1138  silseed->end_cluster_keys());
1139  }
1140  TrackFitUtils::getTrackletClusters(tgeometry, clustermap,
1141  clusterPositions, clusterKeys);
1142  std::vector<float> fitparams = TrackFitUtils::fitClusters(clusterPositions, clusterKeys);
1143  if(fitparams.size() == 0)
1144  {
1145  continue;
1146  }
1147  nclus_all = clusterKeys.size();
1148  float charge = NAN;
1149  if(tpcseed->get_qOverR()>0)
1150  { charge = 1; }
1151  else
1152  { charge = -1; }
1153  // "pt:eta:phi:X0:Y0:charge:nhits:"
1154  float tpt = tpcseed->get_pt();
1155  float teta = tpcseed->get_eta();
1156  float tphi = tpcseed->get_phi(clustermap,tgeometry);
1157  float tX0 = tpcseed->get_X0();
1158  float tY0 = tpcseed->get_Y0();
1159  float tZ0 = tpcseed->get_Z0();
1160 
1161  float nhits_local = 0;
1162  nhits_local += tpcseed->size_cluster_keys();
1163  if(silseed)
1164  {
1165  nhits_local += silseed->size_cluster_keys();
1166  }
1167  nhit_maps = 0;
1168  nhit_intt = 0;
1169  nhit_tpc_all = 0;
1170  nhit_mms = 0;
1171 
1172  for (size_t i = 0; i < clusterPositions.size(); i++)
1173  { /*
1174  for (SvtxTrack::ConstClusterKeyIter iter_local = tpcseed->begin_cluster_keys();
1175  iter_local != tpcseed->end_cluster_keys();
1176  ++iter_local){
1177  */
1178 
1179  const TrkrDefs::cluskey cluster_key = clusterKeys[i]; //*iter_local;
1180  // TrkrCluster* cluster = clustermap->findCluster(cluster_key);
1181  // unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
1182  switch(TrkrDefs::getTrkrId(cluster_key)) {
1183  case TrkrDefs::TrkrId::mvtxId:
1184  nhit_maps++;
1185  break;
1186  case TrkrDefs::TrkrId::inttId:
1187  nhit_intt++;
1188  break;
1189  case TrkrDefs::TrkrId::tpcId:
1190  nhit_tpc_all++;
1191  break;
1192  case TrkrDefs::TrkrId::micromegasId:
1193  nhit_mms++;
1194  break;
1195  default:
1196  break;
1197  }
1198 
1199  Acts::Vector3 position = clusterPositions[i];
1200  Acts::Vector3 pca = TrackFitUtils::get_helix_pca(fitparams, position);
1201  float cluster_phi = atan2(position(1), position(0));
1202  float pca_phi = atan2(pca(1), pca(0));
1203  float dphi = cluster_phi - pca_phi;
1204  if (dphi > M_PI)
1205  {
1206  dphi = 2 * M_PI - dphi;
1207  }
1208  if (dphi < -M_PI)
1209  {
1210  dphi = 2 * M_PI + dphi;
1211  }
1212  float dz = position(2) - pca(2);
1213 
1214  TrkrCluster* cluster = clustermap->findCluster(cluster_key);
1215  Acts::Vector3 cglob;
1216  cglob = tgeometry->getGlobalPosition(cluster_key, cluster);
1217  float x = cglob(0);
1218  float y = cglob(1);
1219  float z = cglob(2);
1220  float locx = cluster->getLocalX();
1221  float locy = cluster->getLocalY();
1222  TVector3 pos(x, y, z);
1223  float r = pos.Perp();
1224  float phi = pos.Phi();
1225  float eta = pos.Eta();
1226  float theta = pos.Theta();
1227 
1228  float ex = 0;
1229  float ey = 0;
1230  float ez = 0;
1231  float ephi = 0;
1232  float pez = 0;
1233  float pephi = 0;
1234  float size = 0;
1235  float phisize = 0;
1236  float zsize = 0;
1237  float maxadc = -999;
1238 
1239  auto para_errors = ClusErrPara.get_clusterv5_modified_error(cluster, r, cluster_key);
1240 
1241  phisize = cluster->getPhiSize();
1242  zsize = cluster->getZSize();
1243  // double clusRadius = r;
1244  ez = sqrt(para_errors.second);
1245  ephi = sqrt(para_errors.first);
1246  maxadc = cluster->getMaxAdc();
1247 
1248  float e = cluster->getAdc();
1249  float adc = cluster->getAdc();
1250  float layer_local = (float) TrkrDefs::getLayer(cluster_key);
1251  float sector = TpcDefs::getSectorId(cluster_key);
1252  float side = TpcDefs::getSide(cluster_key);
1253 
1254  float trackID = phtrk_iter;
1255  float phibin = NAN;
1256  float tbin = NAN;
1257  const TrkrDefs::hitsetkey hsk = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
1258  if (TrkrDefs::getTrkrId(hsk) == TrkrDefs::TrkrId::tpcId)
1259  {
1260  PHG4TpcCylinderGeom* GeoLayer_local = geom_container->GetLayerCellGeom(layer_local);
1261  phibin = GeoLayer_local->get_pad_float(phi, side);
1262  tbin = GeoLayer_local->get_tbin_float(locy - 39.6);
1263  }
1264  else
1265  {
1266  phibin = locx;
1267  tbin = locy;
1268  }
1269  /*
1270  "pt:eta:phi:X0:Y0:charge:nhits:"
1271  */
1272  float clus_trk_data[] = {(float) _ievent,
1273  (float) phtrk_iter,
1274  locx,
1275  locy,
1276  x,
1277  y,
1278  z,
1279  r,
1280  phi,
1281  eta,
1282  theta,
1283  phibin,
1284  tbin,
1285  ex,
1286  ey,
1287  ez,
1288  ephi,
1289  pez,
1290  pephi,
1291  e,
1292  adc,
1293  maxadc,
1294  layer_local,
1295  sector,
1296  side,
1297  size,
1298  phisize,
1299  zsize,
1300  dphi,
1301  dz,
1302  trackID,
1303  0,
1304  tpt,
1305  teta,
1306  tphi,
1307  tX0,
1308  tY0,
1309  tZ0,
1310  charge,
1311  nhit_maps,
1312  nhit_intt,
1313  nhit_mms,
1314  nhits_local,
1315  nhit_tpc_all,
1316  nhit_tpc_in,
1317  nhit_tpc_mid,
1318  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1319 
1320  _ntp_clus_trk->Fill(clus_trk_data);
1321  }
1322  }
1323  }
1324 
1325 
1326  //------------------------
1327  // fill the Track NTuple
1328  //------------------------
1329 
1330  if (_ntp_track)
1331  {
1332  if (Verbosity() > 1)
1333  {
1334  cout << "Filling ntp_track " << endl;
1335  _timer->restart();
1336  }
1337 
1338  // need things off of the DST...
1339 
1340  GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
1341 
1342  if (_trackmap)
1343  {
1344  for (auto& iter : *_trackmap)
1345  {
1346  SvtxTrack* track = iter.second;
1347  if(!track) continue;
1348  float trackID = track->get_id();
1349  TrackSeed* tpcseed = track->get_tpc_seed();
1350  TrackSeed* silseed = track->get_silicon_seed();
1351  short int crossing_int = track->get_crossing();
1352  float crossing;
1353  if (crossing_int == SHRT_MAX)
1354  {
1355  crossing = NAN;
1356  }
1357  else
1358  {
1359  crossing = (float) crossing_int;
1360  }
1361  float charge = track->get_charge();
1362  float quality = track->get_quality();
1363  float chisq = track->get_chisq();
1364  float ndf = track->get_ndf();
1365  float nhits_local = 0;
1366  if (tpcseed)
1367  {
1368  nhits_local += tpcseed->size_cluster_keys();
1369  }
1370  if (silseed)
1371  {
1372  nhits_local += silseed->size_cluster_keys();
1373  }
1374  unsigned int layers = 0x0;
1375  int maps[_nlayers_maps];
1376  int intt[_nlayers_intt];
1377  int tpc[_nlayers_tpc];
1378  int mms[_nlayers_mms];
1379  if (_nlayers_maps > 0)
1380  {
1381  for (unsigned int i = 0; i < _nlayers_maps; i++)
1382  {
1383  maps[i] = 0;
1384  }
1385  }
1386  if (_nlayers_intt > 0)
1387  {
1388  for (unsigned int i = 0; i < _nlayers_intt; i++)
1389  {
1390  intt[i] = 0;
1391  }
1392  }
1393  if (_nlayers_tpc > 0)
1394  {
1395  for (unsigned int i = 0; i < _nlayers_tpc; i++)
1396  {
1397  tpc[i] = 0;
1398  }
1399  }
1400  if (_nlayers_mms > 0)
1401  {
1402  for (unsigned int i = 0; i < _nlayers_mms; i++)
1403  {
1404  mms[i] = 0;
1405  }
1406  }
1407 
1408  float nmaps = 0;
1409  float nintt = 0;
1410  float nmms = 0;
1411  float ntpc = 0;
1412  float ntpc1 = 0;
1413  float ntpc11 = 0;
1414  float ntpc2 = 0;
1415  float ntpc3 = 0;
1416  float nlmaps = 0;
1417  float nlintt = 0;
1418  float nltpc = 0;
1419  float nlmms = 0;
1420 
1421  if (tpcseed)
1422  {
1423  for (SvtxTrack::ConstClusterKeyIter iter_local = tpcseed->begin_cluster_keys();
1424  iter_local != tpcseed->end_cluster_keys();
1425  ++iter_local)
1426  {
1427  TrkrDefs::cluskey cluster_key = *iter_local;
1428  // TrkrCluster* cluster = clustermap->findCluster(cluster_key);
1429  unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
1430 
1431  if (_nlayers_maps > 0 && layer_local < _nlayers_maps)
1432  {
1433  maps[layer_local] = 1;
1434  nmaps++;
1435  }
1436  if (_nlayers_intt > 0 && layer_local >= _nlayers_maps && layer_local < _nlayers_maps + _nlayers_intt)
1437  {
1438  intt[layer_local - _nlayers_maps] = 1;
1439  nintt++;
1440  }
1441  if (_nlayers_mms > 0 && layer_local >= _nlayers_maps + _nlayers_intt + _nlayers_tpc && layer_local < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
1442  {
1443  mms[layer_local - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
1444  nmms++;
1445  }
1446  if (_nlayers_tpc > 0 && layer_local >= (_nlayers_maps + _nlayers_intt) && layer_local < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
1447  {
1448  tpc[layer_local - (_nlayers_maps + _nlayers_intt)] = 1;
1449  ntpc++;
1450  if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 8)
1451  {
1452  ntpc11++;
1453  }
1454 
1455  if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 16)
1456  {
1457  ntpc1++;
1458  }
1459  else if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 32)
1460  {
1461  ntpc2++;
1462  }
1463  else if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 48)
1464  {
1465  ntpc3++;
1466  }
1467  }
1468  }
1469  }
1470 
1471  if (silseed)
1472  {
1473  for (SvtxTrack::ConstClusterKeyIter iter_local = silseed->begin_cluster_keys();
1474  iter_local != silseed->end_cluster_keys();
1475  ++iter_local)
1476  {
1477  TrkrDefs::cluskey cluster_key = *iter_local;
1478  // TrkrCluster* cluster = clustermap->findCluster(cluster_key);
1479  unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
1480 
1481  if (_nlayers_maps > 0 && layer_local < _nlayers_maps)
1482  {
1483  maps[layer_local] = 1;
1484  nmaps++;
1485  }
1486  if (_nlayers_intt > 0 && layer_local >= _nlayers_maps && layer_local < _nlayers_maps + _nlayers_intt)
1487  {
1488  intt[layer_local - _nlayers_maps] = 1;
1489  nintt++;
1490  }
1491  if (_nlayers_mms > 0 && layer_local >= _nlayers_maps + _nlayers_intt + _nlayers_tpc && layer_local < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
1492  {
1493  mms[layer_local - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
1494  nmms++;
1495  }
1496  if (_nlayers_tpc > 0 && layer_local >= (_nlayers_maps + _nlayers_intt) && layer_local < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
1497  {
1498  tpc[layer_local - (_nlayers_maps + _nlayers_intt)] = 1;
1499  ntpc++;
1500  if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 8)
1501  {
1502  ntpc11++;
1503  }
1504 
1505  if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 16)
1506  {
1507  ntpc1++;
1508  }
1509  else if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 32)
1510  {
1511  ntpc2++;
1512  }
1513  else if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 48)
1514  {
1515  ntpc3++;
1516  }
1517  }
1518  }
1519  }
1520 
1521  if (_nlayers_maps > 0)
1522  {
1523  for (unsigned int i = 0; i < _nlayers_maps; i++)
1524  {
1525  nlmaps += maps[i];
1526  }
1527  }
1528  if (_nlayers_intt > 0)
1529  {
1530  for (unsigned int i = 0; i < _nlayers_intt; i++)
1531  {
1532  nlintt += intt[i];
1533  }
1534  }
1535  if (_nlayers_tpc > 0)
1536  {
1537  for (unsigned int i = 0; i < _nlayers_tpc; i++)
1538  {
1539  nltpc += tpc[i];
1540  }
1541  }
1542  if (_nlayers_mms > 0)
1543  {
1544  for (unsigned int i = 0; i < _nlayers_mms; i++)
1545  {
1546  nlmms += mms[i];
1547  }
1548  }
1549  layers = nlmaps + nlintt + nltpc + nlmms;
1550 
1551  float dca3dxy = NAN, dca3dz = NAN,
1552  dca3dxysigma = NAN, dca3dzsigma = NAN;
1553  float dca2d = NAN, dca2dsigma = NAN;
1554 
1555  int vertexID = track->get_vertex_id();
1556  float vx = NAN;
1557  float vy = NAN;
1558  float vz = NAN;
1559  if(vertexmap)
1560  {
1561  GlobalVertex* vertex = vertexmap->get(vertexID);
1562 
1563  if (vertex)
1564  {
1565  vx = vertex->get_x();
1566  vy = vertex->get_y();
1567  vz = vertex->get_z();
1568  Acts::Vector3 vert(vx,vy,vz);
1569  auto dcapair = TrackAnalysisUtils::get_dca(track, vert);
1570  dca3dxy = dcapair.first.first;
1571  dca3dxysigma = dcapair.first.second;
1572  dca3dz = dcapair.second.first;
1573  dca3dzsigma = dcapair.second.second;
1574  }
1575  }
1576  float px = track->get_px();
1577  float py = track->get_py();
1578  float pz = track->get_pz();
1579  TVector3 v(px, py, pz);
1580  float pt = v.Pt();
1581  float eta = v.Eta();
1582  float phi = v.Phi();
1583  float CVxx = track->get_error(3, 3);
1584  float CVxy = track->get_error(3, 4);
1585  float CVxz = track->get_error(3, 5);
1586  float CVyy = track->get_error(4, 4);
1587  float CVyz = track->get_error(4, 5);
1588  float CVzz = track->get_error(5, 5);
1589  float deltapt = sqrt((CVxx * px * px + 2 * CVxy * px * py + CVyy * py * py) / (px * px + py * py));
1590  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)));
1591  float deltaphi = sqrt((CVyy * px * px - 2 * CVxy * px * py + CVxx * py * py) / ((px * px + py * py) * (px * px + py * py)));
1592  float pcax = track->get_x();
1593  float pcay = track->get_y();
1594  float pcaz = track->get_z();
1595 
1596  float track_data[] = {(float) _ievent, m_fSeed,
1597  trackID,
1598  crossing,
1599  px,
1600  py,
1601  pz,
1602  pt,
1603  eta,
1604  phi,
1605  deltapt,
1606  deltaeta,
1607  deltaphi,
1608  charge,
1609  quality,
1610  chisq,
1611  ndf,
1612  nhits_local, nmaps, nintt, ntpc, nmms,
1613  ntpc1, ntpc11, ntpc2, ntpc3,
1614  nlmaps, nlintt, nltpc, nlmms,
1615  (float) layers,
1616  (float) vertexID,
1617  vx,
1618  vy,
1619  vz,
1620  dca2d,
1621  dca2dsigma,
1622  dca3dxy,
1623  dca3dxysigma,
1624  dca3dz,
1625  dca3dzsigma,
1626  pcax,
1627  pcay,
1628  pcaz,
1629  nhit_maps,
1630  nhit_intt,
1631  nhit_mms,
1632  nhit_tpc_all,
1633  nhit_tpc_in,
1634  nhit_tpc_mid,
1635  nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1636 
1637  if (Verbosity() >= 1)
1638  {
1639  cout << "ievent " << _ievent
1640  << " trackID " << trackID
1641  << " nhits " << nhits
1642  << " px " << px
1643  << " py " << py
1644  << " pz " << pz
1645  << endl;
1646  }
1647 
1648  _ntp_track->Fill(track_data);
1649  }
1650  }
1651  if (Verbosity() > 1)
1652  {
1653  _timer->stop();
1654  cout << "track time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
1655  }
1656  }
1657 
1658  //---------------------
1659  // fill the Gseed NTuple
1660  //---------------------
1661 
1662  if (_ntp_tpcseed)
1663  {
1664  if (Verbosity() > 1)
1665  {
1666  cout << "Filling ntp_tpcseed " << endl;
1667  _timer->restart();
1668  }
1669 
1670  if (Verbosity() > 1)
1671  {
1672  _timer->stop();
1673  cout << "tpcseed time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
1674  }
1675  }
1676 
1677  if (_ntp_siseed)
1678  {
1679  if (Verbosity() > 1)
1680  {
1681  cout << "Filling ntp_tpcseed " << endl;
1682  _timer->restart();
1683  }
1684 
1685  if (Verbosity() > 1)
1686  {
1687  _timer->stop();
1688  cout << "tpcseed time: " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
1689  }
1690  }
1691  return;
1692 }
1693 
1694 std::vector<TrkrDefs::cluskey> TrkrNtuplizer::get_track_ckeys(SvtxTrack* track)
1695 {
1696  std::vector<TrkrDefs::cluskey> cluster_keys;
1697  TrackSeed *tpcseed = track->get_tpc_seed();
1698  TrackSeed *silseed = track->get_silicon_seed();
1699  if(silseed)
1700  {
1701  for(auto iter = silseed->begin_cluster_keys();
1702  iter!= silseed->end_cluster_keys();
1703  ++iter)
1704  { cluster_keys.push_back(*iter); }
1705  }
1706  if(tpcseed)
1707  {
1708  for(auto iter = tpcseed->begin_cluster_keys();
1709  iter!= tpcseed->end_cluster_keys();
1710  ++iter)
1711  { cluster_keys.push_back(*iter); }
1712  }
1713 
1714  return cluster_keys;
1715 }
1716 
1717 SvtxTrack* TrkrNtuplizer::best_track_from(TrkrDefs::cluskey cluster_key)
1718 {
1719  std::map<TrkrDefs::cluskey, SvtxTrack*>::iterator find_iter =
1720  _cache_best_track_from_cluster.find(cluster_key);
1721  if (find_iter != _cache_best_track_from_cluster.end())
1722  {
1723  return find_iter->second;
1724  }
1725 
1726  SvtxTrack* best_track = nullptr;
1727  float best_quality = FLT_MAX;
1728 
1729  std::set<SvtxTrack*> tracks = all_tracks_from(cluster_key);
1730  // loop over all SvtxTracks
1731  for (std::set<SvtxTrack*>::iterator iter = tracks.begin();
1732  iter != tracks.end();
1733  ++iter)
1734  {
1735  SvtxTrack* candidate = *iter;
1736  if (candidate->get_quality() < best_quality)
1737  {
1738  best_quality = candidate->get_quality();
1739  best_track = candidate;
1740  }
1741  }
1742 
1743  _cache_best_track_from_cluster.insert(make_pair(cluster_key, best_track));
1744  return best_track;
1745 }
1746 
1747 std::set<SvtxTrack*> TrkrNtuplizer::all_tracks_from(TrkrDefs::cluskey cluster_key)
1748 {
1749  std::set<SvtxTrack*> tracks;
1750 
1751  if (_cache_track_from_cluster_exists == false) create_cache_track_from_cluster();
1752  std::map<TrkrDefs::cluskey, std::set<SvtxTrack*> >::iterator find_iter =
1753  _cache_all_tracks_from_cluster.find(cluster_key);
1754  if (find_iter != _cache_all_tracks_from_cluster.end())
1755  {
1756  return find_iter->second;
1757  }
1758  else
1759  {
1760  return tracks;
1761  }
1762 
1763  // loop over all SvtxTracks
1764  for (SvtxTrackMap::Iter iter = _trackmap->begin();
1765  iter != _trackmap->end();
1766  ++iter)
1767  {
1768  SvtxTrack* track = iter->second;
1769  std::vector<TrkrDefs::cluskey> cluster_keys = get_track_ckeys(track);
1770 
1771  // loop over all clusters
1772  for (const auto& candidate : cluster_keys)
1773  {
1774 
1775  // if (_strict)
1776  // {
1777  // assert(candidate);
1778  // }
1779  // else if (!candidate)
1780  // {
1781  // ++_errors;
1782  // continue;
1783  // }
1784 
1785  if (cluster_key == candidate)
1786  {
1787  tracks.insert(track);
1788  }
1789  }
1790  }
1791 
1792  _cache_all_tracks_from_cluster.insert(make_pair(cluster_key, tracks));
1793 
1794  return tracks;
1795 }
1796 
1797 void TrkrNtuplizer::create_cache_track_from_cluster()
1798 {
1799 
1800  if(!_trackmap) return;
1801 
1802  // loop over all SvtxTracks
1803  for (SvtxTrackMap::Iter iter = _trackmap->begin();
1804  iter != _trackmap->end();
1805  ++iter)
1806  {
1807  SvtxTrack* track = iter->second;
1808  std::vector<TrkrDefs::cluskey> cluster_keys = get_track_ckeys(track);
1809 
1810  // loop over all clusters
1811  for (const auto& candidate_key : cluster_keys)
1812  {
1813  //check if cluster has an entry in cache
1814  std::map<TrkrDefs::cluskey, std::set<SvtxTrack*> >::iterator cliter =
1815  _cache_all_tracks_from_cluster.find(candidate_key);
1816  if (cliter != _cache_all_tracks_from_cluster.end())
1817  { //got entry
1818  cliter->second.insert(track); //add track to list;
1819  }
1820  else
1821  {
1822  std::set<SvtxTrack*> tracks;
1823  tracks.insert(track);
1824  _cache_all_tracks_from_cluster.insert(make_pair(candidate_key, tracks));
1825  }
1826  }
1827  }
1828  _cache_track_from_cluster_exists = true;
1829 
1830  return;
1831 }
1832 
1833 TMatrixF TrkrNtuplizer::calculateClusterError(TrkrCluster* c, float& clusphi)
1834 {
1835  TMatrixF localErr(3, 3);
1836  localErr[0][0] = 0.;
1837  localErr[0][1] = 0.;
1838  localErr[0][2] = 0.;
1839  localErr[1][0] = 0.;
1840  localErr[1][1] = c->getActsLocalError(0, 0);
1841  localErr[1][2] = c->getActsLocalError(0, 1);
1842  localErr[2][0] = 0.;
1843  localErr[2][1] = c->getActsLocalError(1, 0);
1844  localErr[2][2] = c->getActsLocalError(1, 1);
1845 
1846  TMatrixF ROT(3, 3);
1847  ROT[0][0] = cos(clusphi);
1848  ROT[0][1] = -sin(clusphi);
1849  ROT[0][2] = 0.0;
1850  ROT[1][0] = sin(clusphi);
1851  ROT[1][1] = cos(clusphi);
1852  ROT[1][2] = 0.0;
1853  ROT[2][0] = 0.0;
1854  ROT[2][1] = 0.0;
1855  ROT[2][2] = 1.0;
1856  TMatrixF ROT_T(3, 3);
1857  ROT_T.Transpose(ROT);
1858 
1859  TMatrixF err(3, 3);
1860  err = ROT * localErr * ROT_T;
1861  return err;
1862 }