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