ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CaloEvaluator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CaloEvaluator.cc
1 #include "CaloEvaluator.h"
2 
3 #include "CaloEvalStack.h"
4 #include "CaloRawClusterEval.h"
5 #include "CaloRawTowerEval.h"
6 #include "CaloTruthEval.h"
7 
8 #include <g4main/PHG4Particle.h>
9 #include <g4main/PHG4Shower.h>
10 #include <g4main/PHG4TruthInfoContainer.h>
11 #include <g4main/PHG4VtxPoint.h>
12 
13 #include <g4vertex/GlobalVertex.h>
14 #include <g4vertex/GlobalVertexMap.h>
15 
16 #include <calobase/RawCluster.h>
17 #include <calobase/RawClusterContainer.h>
18 #include <calobase/RawClusterUtility.h>
19 #include <calobase/RawTower.h>
20 #include <calobase/RawTowerContainer.h>
21 #include <calobase/RawTowerGeom.h>
22 #include <calobase/RawTowerGeomContainer.h>
23 
24 #include <fun4all/Fun4AllReturnCodes.h>
25 #include <fun4all/SubsysReco.h>
26 
27 #include <phool/getClass.h>
28 #include <phool/phool.h>
29 
30 #include <TFile.h>
31 #include <TNtuple.h>
32 #include <TTree.h>
33 
34 #include <CLHEP/Vector/ThreeVector.h>
35 
36 #include <cmath>
37 #include <cstdlib>
38 #include <iostream>
39 #include <set>
40 #include <utility>
41 
42 using namespace std;
43 
44 CaloEvaluator::CaloEvaluator(const string& name, const string& caloname, const string& filename)
45  : SubsysReco(name)
46  , _caloname(caloname)
47  , _ievent(0)
48  , _towerID_debug(0)
49  , _ieta_debug(0)
50  , _iphi_debug(0)
51  , _eta_debug(0)
52  , _phi_debug(0)
53  , _e_debug(0)
54  , _x_debug(0)
55  , _y_debug(0)
56  , _z_debug(0)
57  , _truth_trace_embed_flags()
58  , _truth_e_threshold(0.0)
59  , // 0 GeV before reco is traced
60  _reco_e_threshold(0.0)
61  , // 0 GeV before reco is traced
62  _caloevalstack(nullptr)
63  , _strict(false)
64  , _do_gpoint_eval(true)
65  , _do_gshower_eval(true)
66  , _do_tower_eval(true)
67  , _do_cluster_eval(true)
68  , _ntp_gpoint(nullptr)
69  , _ntp_gshower(nullptr)
70  , _ntp_tower(nullptr)
71  , _tower_debug(nullptr)
72  , _ntp_cluster(nullptr)
73  , _filename(filename)
74  , _tfile(nullptr)
75 {
76 }
77 
78 int CaloEvaluator::Init(PHCompositeNode* /*topNode*/)
79 {
80  _ievent = 0;
81 
82  _tfile = new TFile(_filename.c_str(), "RECREATE");
83 
84  if (_do_gpoint_eval) _ntp_gpoint = new TNtuple("ntp_gpoint", "primary vertex => best (first) vertex",
85  "event:gvx:gvy:gvz:"
86  "vx:vy:vz");
87 
88  if (_do_gshower_eval) _ntp_gshower = new TNtuple("ntp_gshower", "truth shower => best cluster",
89  "event:gparticleID:gflavor:gnhits:"
90  "geta:gphi:ge:gpt:gvx:gvy:gvz:gembed:gedep:"
91  "clusterID:ntowers:eta:x:y:z:phi:e:efromtruth");
92 
93  //Barak: Added TTree to will allow the TowerID to be set correctly as integer
94  if (_do_tower_eval)
95  {
96  _ntp_tower = new TNtuple("ntp_tower", "tower => max truth primary",
97  "event:towerID:ieta:iphi:eta:phi:e:x:y:z:"
98  "gparticleID:gflavor:gnhits:"
99  "geta:gphi:ge:gpt:gvx:gvy:gvz:"
100  "gembed:gedep:"
101  "efromtruth");
102 
103  //Make Tree
104  _tower_debug = new TTree("tower_debug", "tower => max truth primary");
105 
106  _tower_debug->Branch("event", &_ievent, "event/I");
107  _tower_debug->Branch("towerID", &_towerID_debug, "towerID/I");
108  _tower_debug->Branch("ieta", &_ieta_debug, "ieta/I");
109  _tower_debug->Branch("iphi", &_iphi_debug, "iphi/I");
110  _tower_debug->Branch("eta", &_eta_debug, "eta/F");
111  _tower_debug->Branch("phi", &_phi_debug, "phi/F");
112  _tower_debug->Branch("e", &_e_debug, "e/F");
113  _tower_debug->Branch("x", &_x_debug, "x/F");
114  _tower_debug->Branch("y", &_y_debug, "y/F");
115  _tower_debug->Branch("z", &_z_debug, "z/F");
116  }
117 
118  if (_do_cluster_eval) _ntp_cluster = new TNtuple("ntp_cluster", "cluster => max truth primary",
119  "event:clusterID:ntowers:eta:x:y:z:phi:e:"
120  "gparticleID:gflavor:gnhits:"
121  "geta:gphi:ge:gpt:gvx:gvy:gvz:"
122  "gembed:gedep:"
123  "efromtruth");
124 
125  return Fun4AllReturnCodes::EVENT_OK;
126 }
127 
128 int CaloEvaluator::process_event(PHCompositeNode* topNode)
129 {
130  if (!_caloevalstack)
131  {
132  _caloevalstack = new CaloEvalStack(topNode, _caloname);
133  _caloevalstack->set_strict(_strict);
134  _caloevalstack->set_verbosity(Verbosity() + 1);
135  }
136  else
137  {
138  _caloevalstack->next_event(topNode);
139  }
140 
141  //-----------------------------------
142  // print what is coming into the code
143  //-----------------------------------
144 
145  printInputInfo(topNode);
146 
147  //---------------------------
148  // fill the Evaluator NTuples
149  //---------------------------
150 
151  fillOutputNtuples(topNode);
152 
153  //--------------------------------------------------
154  // Print out the ancestry information for this event
155  //--------------------------------------------------
156 
157  printOutputInfo(topNode);
158 
159  ++_ievent;
160 
161  return Fun4AllReturnCodes::EVENT_OK;
162 }
163 
164 int CaloEvaluator::End(PHCompositeNode* /*topNode*/)
165 {
166  _tfile->cd();
167 
168  if (_do_gpoint_eval) _ntp_gpoint->Write();
169  if (_do_gshower_eval) _ntp_gshower->Write();
170  if (_do_tower_eval)
171  {
172  _ntp_tower->Write();
173  _tower_debug->Write();
174  }
175  if (_do_cluster_eval) _ntp_cluster->Write();
176 
177  _tfile->Close();
178 
179  delete _tfile;
180 
181  if (Verbosity() > 0)
182  {
183  cout << "========================= " << Name() << "::End() ============================" << endl;
184  cout << " " << _ievent << " events of output written to: " << _filename << endl;
185  cout << "===========================================================================" << endl;
186  }
187 
188  if (_caloevalstack) delete _caloevalstack;
189 
190  return Fun4AllReturnCodes::EVENT_OK;
191 }
192 
193 void CaloEvaluator::printInputInfo(PHCompositeNode* topNode)
194 {
195  if (Verbosity() > 2) cout << "CaloEvaluator::printInputInfo() entered" << endl;
196 
197  // print out the truth container
198 
199  if (Verbosity() > 1)
200  {
201  cout << endl;
202  cout << PHWHERE << " NEW INPUT FOR EVENT " << _ievent << endl;
203  cout << endl;
204 
205  // need things off of the DST...
206  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
207  if (!truthinfo)
208  {
209  cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
210  exit(-1);
211  }
212 
213  cout << Name() << ": PHG4TruthInfoContainer contents: " << endl;
214 
215  PHG4TruthInfoContainer::Range truthrange = truthinfo->GetParticleRange();
216  for (PHG4TruthInfoContainer::Iterator truthiter = truthrange.first;
217  truthiter != truthrange.second;
218  ++truthiter)
219  {
220  PHG4Particle* particle = truthiter->second;
221 
222  cout << truthiter->first << " => pid: " << particle->get_pid()
223  << " pt: " << sqrt(pow(particle->get_px(), 2) + pow(particle->get_py(), 2)) << endl;
224  }
225  }
226 
227  return;
228 }
229 
230 void CaloEvaluator::printOutputInfo(PHCompositeNode* topNode)
231 {
232  if (Verbosity() > 2) cout << "CaloEvaluator::printOutputInfo() entered" << endl;
233 
234  CaloRawClusterEval* clustereval = _caloevalstack->get_rawcluster_eval();
235  CaloTruthEval* trutheval = _caloevalstack->get_truth_eval();
236 
237  //==========================================
238  // print out some useful stuff for debugging
239  //==========================================
240 
241  if (Verbosity() > 1)
242  {
243  // event information
244  cout << endl;
245  cout << PHWHERE << " NEW OUTPUT FOR EVENT " << _ievent << endl;
246  cout << endl;
247 
248  // need things off of the DST...
249  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
250  if (!truthinfo)
251  {
252  cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
253  exit(-1);
254  }
255 
256  // need things off of the DST...
257  GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
258 
259  PHG4VtxPoint* gvertex = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
260  float gvx = gvertex->get_x();
261  float gvy = gvertex->get_y();
262  float gvz = gvertex->get_z();
263 
264  float vx = NAN;
265  float vy = NAN;
266  float vz = NAN;
267  if (vertexmap)
268  {
269  if (!vertexmap->empty())
270  {
271  GlobalVertex* vertex = (vertexmap->begin()->second);
272 
273  vx = vertex->get_x();
274  vy = vertex->get_y();
275  vz = vertex->get_z();
276  }
277  }
278 
279  cout << "vtrue = (" << gvx << "," << gvy << "," << gvz << ") => vreco = (" << vx << "," << vy << "," << vz << ")" << endl;
280 
281  PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
282  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
283  iter != range.second;
284  ++iter)
285  {
286  PHG4Particle* primary = iter->second;
287 
288  cout << endl;
289 
290  cout << "===Primary PHG4Particle=========================================" << endl;
291  cout << " particle id = " << primary->get_track_id() << endl;
292  cout << " flavor = " << primary->get_pid() << endl;
293  cout << " (px,py,pz,e) = (";
294 
295  float gpx = primary->get_px();
296  float gpy = primary->get_py();
297  float gpz = primary->get_pz();
298  float ge = primary->get_e();
299 
300  cout.width(5);
301  cout << gpx;
302  cout << ",";
303  cout.width(5);
304  cout << gpy;
305  cout << ",";
306  cout.width(5);
307  cout << gpz;
308  cout << ",";
309  cout.width(5);
310  cout << ge;
311  cout << ")" << endl;
312 
313  float gpt = sqrt(gpx * gpx + gpy * gpy);
314  float geta = NAN;
315  if (gpt != 0.0) geta = asinh(gpz / gpt);
316  float gphi = atan2(gpy, gpx);
317 
318  cout << "(eta,phi,e,pt) = (";
319  cout.width(5);
320  cout << geta;
321  cout << ",";
322  cout.width(5);
323  cout << gphi;
324  cout << ",";
325  cout.width(5);
326  cout << ge;
327  cout << ",";
328  cout.width(5);
329  cout << gpt;
330  cout << ")" << endl;
331 
332  PHG4VtxPoint* vtx = trutheval->get_vertex(primary);
333  float gvx = vtx->get_x();
334  float gvy = vtx->get_y();
335  float gvz = vtx->get_z();
336 
337  cout << " vtrue = (";
338  cout.width(5);
339  cout << gvx;
340  cout << ",";
341  cout.width(5);
342  cout << gvy;
343  cout << ",";
344  cout.width(5);
345  cout << gvz;
346  cout << ")" << endl;
347 
348  cout << " embed = " << trutheval->get_embed(primary) << endl;
349  cout << " edep = " << trutheval->get_shower_energy_deposit(primary) << endl;
350 
351  std::set<RawCluster*> clusters = clustereval->all_clusters_from(primary);
352  for (std::set<RawCluster*>::iterator clusiter = clusters.begin();
353  clusiter != clusters.end();
354  ++clusiter)
355  {
356  RawCluster* cluster = (*clusiter);
357 
358  float ntowers = cluster->getNTowers();
359  float x = cluster->get_x();
360  float y = cluster->get_y();
361  float z = cluster->get_z();
362  float phi = cluster->get_phi();
363  float e = cluster->get_energy();
364 
365  float efromtruth = clustereval->get_energy_contribution(cluster, primary);
366 
367  cout << " => #" << cluster->get_id() << " (x,y,z,phi,e) = (";
368  cout.width(5);
369  cout << x;
370  cout << ",";
371  cout.width(5);
372  cout << y;
373  cout << ",";
374  cout.width(5);
375  cout << z;
376  cout << ",";
377  cout.width(5);
378  cout << phi;
379  cout << ",";
380  cout.width(5);
381  cout << e;
382  cout << "), ntowers = " << ntowers << ", efromtruth = " << efromtruth << endl;
383  }
384  }
385  cout << endl;
386  }
387 
388  return;
389 }
390 
391 void CaloEvaluator::fillOutputNtuples(PHCompositeNode* topNode)
392 {
393  if (Verbosity() > 2) cout << "CaloEvaluator::fillOutputNtuples() entered" << endl;
394 
395  CaloRawClusterEval* clustereval = _caloevalstack->get_rawcluster_eval();
396  CaloRawTowerEval* towereval = _caloevalstack->get_rawtower_eval();
397  CaloTruthEval* trutheval = _caloevalstack->get_truth_eval();
398 
399  //----------------------
400  // fill the Event NTuple
401  //----------------------
402 
403  if (_do_gpoint_eval)
404  {
405  // need things off of the DST...
406  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
407  if (!truthinfo)
408  {
409  cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
410  exit(-1);
411  }
412 
413  // need things off of the DST...
414  GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
415 
416  PHG4VtxPoint* gvertex = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
417  float gvx = gvertex->get_x();
418  float gvy = gvertex->get_y();
419  float gvz = gvertex->get_z();
420 
421  float vx = NAN;
422  float vy = NAN;
423  float vz = NAN;
424  if (vertexmap)
425  {
426  if (!vertexmap->empty())
427  {
428  GlobalVertex* vertex = (vertexmap->begin()->second);
429 
430  vx = vertex->get_x();
431  vy = vertex->get_y();
432  vz = vertex->get_z();
433  }
434  }
435 
436  float gpoint_data[7] = {(float) _ievent,
437  gvx,
438  gvy,
439  gvz,
440  vx,
441  vy,
442  vz};
443 
444  _ntp_gpoint->Fill(gpoint_data);
445  }
446 
447  //------------------------
448  // fill the Gshower NTuple
449  //------------------------
450 
451  if (_ntp_gshower)
452  {
453  if (Verbosity() > 1) cout << Name() << " CaloEvaluator::filling gshower ntuple..." << endl;
454 
455  GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
456 
457  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
458  if (!truthinfo)
459  {
460  cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
461  exit(-1);
462  }
463 
464  PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
465  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
466  iter != range.second;
467  ++iter)
468  {
469  PHG4Particle* primary = iter->second;
470 
471  if (primary->get_e() < _truth_e_threshold) continue;
472 
473  if (!_truth_trace_embed_flags.empty())
474  {
475  if (_truth_trace_embed_flags.find(trutheval->get_embed(primary)) ==
476  _truth_trace_embed_flags.end()) continue;
477  }
478 
479  float gparticleID = primary->get_track_id();
480  float gflavor = primary->get_pid();
481 
482  PHG4Shower* shower = trutheval->get_primary_shower(primary);
483  float gnhits = NAN;
484  if (shower)
485  gnhits = shower->get_nhits(trutheval->get_caloid());
486  else
487  gnhits = 0.0;
488  float gpx = primary->get_px();
489  float gpy = primary->get_py();
490  float gpz = primary->get_pz();
491  float ge = primary->get_e();
492 
493  float gpt = sqrt(gpx * gpx + gpy * gpy);
494  float geta = NAN;
495  if (gpt != 0.0) geta = asinh(gpz / gpt);
496  float gphi = atan2(gpy, gpx);
497 
498  PHG4VtxPoint* vtx = trutheval->get_vertex(primary);
499  float gvx = vtx->get_x();
500  float gvy = vtx->get_y();
501  float gvz = vtx->get_z();
502 
503  float gembed = trutheval->get_embed(primary);
504  float gedep = trutheval->get_shower_energy_deposit(primary);
505 
506  RawCluster* cluster = clustereval->best_cluster_from(primary);
507 
508  float clusterID = NAN;
509  float ntowers = NAN;
510  float eta = NAN;
511  float x = NAN;
512  float y = NAN;
513  float z = NAN;
514  float phi = NAN;
515  float e = NAN;
516 
517  float efromtruth = NAN;
518 
519  if (cluster)
520  {
521  clusterID = cluster->get_id();
522  ntowers = cluster->getNTowers();
523  x = cluster->get_x();
524  y = cluster->get_y();
525  z = cluster->get_z();
526  phi = cluster->get_phi();
527  e = cluster->get_energy();
528 
529  efromtruth = clustereval->get_energy_contribution(cluster, primary);
530 
531  // require vertex for cluster eta calculation
532  if (vertexmap)
533  {
534  if (!vertexmap->empty())
535  {
536  GlobalVertex* vertex = (vertexmap->begin()->second);
537 
538  eta =
539  RawClusterUtility::GetPseudorapidity(
540  *cluster,
541  CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
542  }
543  }
544  }
545 
546  float shower_data[] = {(float) _ievent,
547  gparticleID,
548  gflavor,
549  gnhits,
550  geta,
551  gphi,
552  ge,
553  gpt,
554  gvx,
555  gvy,
556  gvz,
557  gembed,
558  gedep,
559  clusterID,
560  ntowers,
561  eta,
562  x,
563  y,
564  z,
565  phi,
566  e,
567  efromtruth};
568 
569  _ntp_gshower->Fill(shower_data);
570  }
571  }
572 
573  //----------------------
574  // fill the Tower NTuple
575  //----------------------
576 
577  if (_do_tower_eval)
578  {
579  if (Verbosity() > 1) cout << "CaloEvaluator::filling tower ntuple..." << endl;
580 
581  string towernode = "TOWER_CALIB_" + _caloname;
582  RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(topNode, towernode.c_str());
583  if (!towers)
584  {
585  cerr << PHWHERE << " ERROR: Can't find " << towernode << endl;
586  exit(-1);
587  }
588 
589  string towergeomnode = "TOWERGEOM_" + _caloname;
590  RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnode.c_str());
591  if (!towergeom)
592  {
593  cerr << PHWHERE << " ERROR: Can't find " << towergeomnode << endl;
594  exit(-1);
595  }
596 
597  RawTowerContainer::ConstRange begin_end = towers->getTowers();
598  RawTowerContainer::ConstIterator rtiter;
599  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
600  {
601  RawTower* tower = rtiter->second;
602 
603  if (tower->get_energy() < _reco_e_threshold) continue;
604 
605  RawTowerGeom* tower_geom = towergeom->get_tower_geometry(tower->get_id());
606  if (!tower_geom)
607  {
608  cerr << PHWHERE << " ERROR: Can't find tower geometry for this tower hit: ";
609  tower->identify();
610  exit(-1);
611  }
612 
613  //cout<<"Tower ID = "<<tower->get_id()<<" for bin(j,k)= "<<tower->get_bineta()<<","<<tower->get_binphi()<<endl; //Added by Barak
614  const float towerid = tower->get_id();
615  const float ieta = tower->get_bineta();
616  const float iphi = tower->get_binphi();
617  const float eta = tower_geom->get_eta();
618  const float phi = tower_geom->get_phi();
619  const float e = tower->get_energy();
620  const float x = tower_geom->get_center_x();
621  const float y = tower_geom->get_center_y();
622  const float z = tower_geom->get_center_z();
623 
624  //Added by Barak
625  _towerID_debug = tower->get_id();
626  _ieta_debug = tower->get_bineta();
627  _iphi_debug = tower->get_binphi();
628  _eta_debug = tower_geom->get_eta();
629  _phi_debug = tower_geom->get_phi();
630  _e_debug = tower->get_energy();
631  _x_debug = tower_geom->get_center_x();
632  _y_debug = tower_geom->get_center_y();
633  _z_debug = tower_geom->get_center_z();
634 
635  PHG4Particle* primary = towereval->max_truth_primary_particle_by_energy(tower);
636 
637  float gparticleID = NAN;
638  float gflavor = NAN;
639  float gnhits = NAN;
640  float gpx = NAN;
641  float gpy = NAN;
642  float gpz = NAN;
643  float ge = NAN;
644 
645  float gpt = NAN;
646  float geta = NAN;
647  float gphi = NAN;
648 
649  float gvx = NAN;
650  float gvy = NAN;
651  float gvz = NAN;
652 
653  float gembed = NAN;
654  float gedep = NAN;
655 
656  float efromtruth = NAN;
657 
658  if (primary)
659  {
660  gparticleID = primary->get_track_id();
661  gflavor = primary->get_pid();
662 
663  PHG4Shower* shower = trutheval->get_primary_shower(primary);
664  if (shower)
665  gnhits = shower->get_nhits(trutheval->get_caloid());
666  else
667  gnhits = 0.0;
668  gpx = primary->get_px();
669  gpy = primary->get_py();
670  gpz = primary->get_pz();
671  ge = primary->get_e();
672 
673  gpt = sqrt(gpx * gpx + gpy * gpy);
674  if (gpt != 0.0) geta = asinh(gpz / gpt);
675  gphi = atan2(gpy, gpx);
676 
677  PHG4VtxPoint* vtx = trutheval->get_vertex(primary);
678 
679  if (vtx)
680  {
681  gvx = vtx->get_x();
682  gvy = vtx->get_y();
683  gvz = vtx->get_z();
684  }
685 
686  gembed = trutheval->get_embed(primary);
687  gedep = trutheval->get_shower_energy_deposit(primary);
688 
689  efromtruth = towereval->get_energy_contribution(tower, primary);
690  }
691 
692  float tower_data[] = {(float) _ievent,
693  towerid,
694  ieta,
695  iphi,
696  eta,
697  phi,
698  e,
699  x,
700  y,
701  z,
702  gparticleID,
703  gflavor,
704  gnhits,
705  geta,
706  gphi,
707  ge,
708  gpt,
709  gvx,
710  gvy,
711  gvz,
712  gembed,
713  gedep,
714  efromtruth};
715 
716  _ntp_tower->Fill(tower_data);
717  _tower_debug->Fill(); //Added by Barak (see above for explanation)
718  }
719  }
720 
721  //------------------------
722  // fill the Cluster NTuple
723  //------------------------
724 
725  if (_do_cluster_eval)
726  {
727  if (Verbosity() > 1) cout << "CaloEvaluator::filling gcluster ntuple..." << endl;
728 
729  GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
730 
731  string clusternode = "CLUSTER_" + _caloname;
732  RawClusterContainer* clusters = findNode::getClass<RawClusterContainer>(topNode, clusternode.c_str());
733  if (!clusters)
734  {
735  cerr << PHWHERE << " ERROR: Can't find " << clusternode << endl;
736  exit(-1);
737  }
738 
739  // for every cluster
740 
741  for (const auto& iterator : clusters->getClustersMap())
742  {
743  RawCluster* cluster = iterator.second;
744 
745  // for (unsigned int icluster = 0; icluster < clusters->size(); icluster++)
746  // {
747  // RawCluster* cluster = clusters->getCluster(icluster);
748 
749  if (cluster->get_energy() < _reco_e_threshold) continue;
750 
751  float clusterID = cluster->get_id();
752  float ntowers = cluster->getNTowers();
753  float x = cluster->get_x();
754  float y = cluster->get_y();
755  float z = cluster->get_z();
756  float eta = NAN;
757  float phi = cluster->get_phi();
758  float e = cluster->get_energy();
759 
760  // require vertex for cluster eta calculation
761  if (vertexmap)
762  {
763  if (!vertexmap->empty())
764  {
765  GlobalVertex* vertex = (vertexmap->begin()->second);
766 
767  eta =
768  RawClusterUtility::GetPseudorapidity(
769  *cluster,
770  CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
771  }
772  }
773 
774  PHG4Particle* primary = clustereval->max_truth_primary_particle_by_energy(cluster);
775 
776  float gparticleID = NAN;
777  float gflavor = NAN;
778 
779  float gnhits = NAN;
780  float gpx = NAN;
781  float gpy = NAN;
782  float gpz = NAN;
783  float ge = NAN;
784 
785  float gpt = NAN;
786  float geta = NAN;
787  float gphi = NAN;
788 
789  float gvx = NAN;
790  float gvy = NAN;
791  float gvz = NAN;
792 
793  float gembed = NAN;
794  float gedep = NAN;
795 
796  float efromtruth = NAN;
797 
798  if (primary)
799  {
800  gparticleID = primary->get_track_id();
801  gflavor = primary->get_pid();
802 
803  PHG4Shower* shower = trutheval->get_primary_shower(primary);
804  if (shower)
805  gnhits = shower->get_nhits(trutheval->get_caloid());
806  else
807  gnhits = 0.0;
808  gpx = primary->get_px();
809  gpy = primary->get_py();
810  gpz = primary->get_pz();
811  ge = primary->get_e();
812 
813  gpt = sqrt(gpx * gpx + gpy * gpy);
814  if (gpt != 0.0) geta = asinh(gpz / gpt);
815  gphi = atan2(gpy, gpx);
816 
817  PHG4VtxPoint* vtx = trutheval->get_vertex(primary);
818 
819  if (vtx)
820  {
821  gvx = vtx->get_x();
822  gvy = vtx->get_y();
823  gvz = vtx->get_z();
824  }
825 
826  gembed = trutheval->get_embed(primary);
827  gedep = trutheval->get_shower_energy_deposit(primary);
828 
829  efromtruth = clustereval->get_energy_contribution(cluster,
830  primary);
831  }
832 
833  float cluster_data[] = {(float) _ievent,
834  clusterID,
835  ntowers,
836  eta,
837  x,
838  y,
839  z,
840  phi,
841  e,
842  gparticleID,
843  gflavor,
844  gnhits,
845  geta,
846  gphi,
847  ge,
848  gpt,
849  gvx,
850  gvy,
851  gvz,
852  gembed,
853  gedep,
854  efromtruth};
855 
856  _ntp_cluster->Fill(cluster_data);
857  }
858  }
859 
860  return;
861 }