ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
eIDMLInterface.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file eIDMLInterface.cc
1 #include "eIDMLInterface.h"
2 
4 #include <calobase/RawCluster.h>
5 #include <calobase/RawClusterContainer.h>
6 #include <calobase/RawClusterUtility.h>
7 #include <calobase/RawTower.h>
8 #include <calobase/RawTowerContainer.h>
9 #include <calobase/RawTowerGeom.h>
10 #include <calobase/RawTowerGeomContainer.h>
11 #include <calotrigger/CaloTriggerInfo.h>
12 
14 #include <g4jets/Jet.h>
15 #include <g4jets/JetMap.h>
16 
18 #include <calobase/RawCluster.h>
19 #include <calobase/RawClusterContainer.h>
20 #include <calobase/RawTower.h>
21 #include <calobase/RawTowerContainer.h>
22 #include <calobase/RawTowerGeomContainer.h>
23 #include <g4vertex/GlobalVertex.h>
24 #include <g4vertex/GlobalVertexMap.h>
25 #include <trackbase_historic/SvtxTrack.h>
26 #include <trackbase_historic/SvtxTrackMap.h>
27 #include <trackbase_historic/SvtxTrack_FastSim.h>
28 #include <trackbase_historic/SvtxVertex.h>
29 #include <trackbase_historic/SvtxVertexMap.h>
30 //PID includes
31 #include <eicpidbase/EICPIDParticle.h>
32 #include <eicpidbase/EICPIDParticleContainer.h>
34 #include <g4eval/JetEvalStack.h>
35 #include <g4eval/SvtxEvalStack.h>
36 
38 #include <HepMC/GenEvent.h>
39 #include <HepMC/GenVertex.h>
40 #include <phhepmc/PHHepMCGenEvent.h>
41 #include <phhepmc/PHHepMCGenEventMap.h>
42 
44 #include <fun4all/Fun4AllHistoManager.h>
45 #include <fun4all/Fun4AllReturnCodes.h>
46 #include <g4main/PHG4Hit.h>
47 #include <g4main/PHG4Particle.h>
48 #include <g4main/PHG4TruthInfoContainer.h>
49 #include <phool/PHCompositeNode.h>
50 #include <phool/getClass.h>
51 
53 #include <TFile.h>
54 #include <TH1.h>
55 #include <TH2.h>
56 #include <TMath.h>
57 #include <TNtuple.h>
58 #include <TTree.h>
59 #include <TVector3.h>
60 
62 #include <cassert>
63 #include <cmath>
64 #include <sstream>
65 #include <string>
66 
67 using namespace std;
68 
79 eIDMLInterface::eIDMLInterface(const std::vector<std::string> &names, const std::string &filename)
80  : SubsysReco("eIDMLInterface_" + names[0])
81  , _calo_names(names)
82  // , m_TTree_Tower_dEta(m_sizeTowerPatch * m_sizeTowerPatch, 0)
83  // , m_TTree_Tower_dPhi(m_sizeTowerPatch * m_sizeTowerPatch, 0)
84  // , m_TTree_Tower_iEta_patch(m_sizeTowerPatch * m_sizeTowerPatch, 0)
85  // , m_TTree_Tower_iPhi_patch(m_sizeTowerPatch * m_sizeTowerPatch, 0)
86  // , m_TTree_Tower_E(m_sizeTowerPatch * m_sizeTowerPatch, 0)
87  , m_outfilename(filename)
88  , m_hm(nullptr)
89  , m_minjetpt(5.0)
90  , m_mincluspt(0.25)
91  , m_analyzeTracks(true)
92  , m_analyzeClusters(true)
93  , m_analyzeJets(true)
94  , m_analyzeTruth(false)
95 {
98  initializeVariables();
99  initializeTrees();
100 }
101 
105 eIDMLInterface::~eIDMLInterface()
106 {
107 }
108 
112 int eIDMLInterface::Init(PHCompositeNode *topNode)
113 {
114  if (Verbosity() > 5)
115  {
116  cout << "Beginning Init in eIDMLInterface" << endl;
117  }
118 
119  m_outfile = new TFile(m_outfilename.c_str(), "RECREATE");
120 
121  m_phi_h = new TH1D("phi_h", ";Counts;#phi [rad]", 50, -6, 6);
122  m_eta_phi_h = new TH2F("phi_eta_h", ";#eta;#phi [rad]", 10, -1, 1, 50, -6, 6);
123 
124  return 0;
125 }
126 
131 int eIDMLInterface::process_event(PHCompositeNode *topNode)
132 {
133  if (Verbosity() > 5)
134  {
135  cout << "Beginning process_event in eIDMLInterface" << endl;
136  }
137 
139  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
140  SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "TrackMap");
141 
142  if (!truthinfo)
143  {
144  cout << PHWHERE
145  << "PHG4TruthInfoContainer node is missing, can't collect G4 truth particles"
146  << endl;
147  return Fun4AllReturnCodes::ABORTRUN;
148  }
149 
151  PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
153  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
154  iter != range.second;
155  ++iter)
156  {
157  initializeVariables();
158 
160  const PHG4Particle *truth = iter->second;
161 
163  m_truthpx = truth->get_px();
164  m_truthpy = truth->get_py();
165  m_truthpz = truth->get_pz();
166  m_truthp = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy + m_truthpz * m_truthpz);
167  m_truthenergy = truth->get_e();
168 
169  m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
170 
171  m_truthphi = atan2(m_truthpy, m_truthpx);
172 
173  m_trutheta = atanh(m_truthpz / m_truthenergy);
174 
175  // eta cut
176  if (m_trutheta < m_etaRange.first or m_trutheta > m_etaRange.second) continue;
177 
179  if (m_trutheta != m_trutheta)
180  m_trutheta = -99;
181  m_truthpid = truth->get_pid();
182 
183  SvtxTrack_FastSim *track = nullptr;
184 
185  if (Verbosity() > 1) cout << __PRETTY_FUNCTION__ << "TRACKmap size " << trackmap->size() << std::endl;
186  for (SvtxTrackMap::ConstIter track_itr = trackmap->begin();
187  track_itr != trackmap->end();
188  track_itr++)
189  {
190  //std::cout << "TRACK * " << track_itr->first << std::endl;
191  SvtxTrack_FastSim *temp = dynamic_cast<SvtxTrack_FastSim *>(track_itr->second);
192  if (!temp)
193  {
194  if (Verbosity() > 1)
195  {
196  std::cout << "PHG4TrackFastSimEval::fill_track_tree - ignore track that is not a SvtxTrack_FastSim:";
197  track_itr->second->identify();
198  }
199  continue;
200  }
201  if (Verbosity() > 1) cout << __PRETTY_FUNCTION__ << " PARTICLE!" << std::endl;
202 
203  if ((temp->get_truth_track_id() - truth->get_track_id()) == 0)
204  {
205  track = temp;
206 
207  break;
208  }
209  }
210 
211  if (track)
212  {
213  // track matched
214 
216  m_tr_px = track->get_px();
217  m_tr_py = track->get_py();
218  m_tr_pz = track->get_pz();
219  m_tr_p = sqrt(m_tr_px * m_tr_px + m_tr_py * m_tr_py + m_tr_pz * m_tr_pz);
220 
221  m_tr_pt = sqrt(m_tr_px * m_tr_px + m_tr_py * m_tr_py);
222 
223  m_tr_phi = track->get_phi();
224  m_tr_eta = track->get_eta();
225 
226  m_charge = track->get_charge();
227  m_chisq = track->get_chisq();
228  m_ndf = track->get_ndf();
229  m_dca = track->get_dca();
230  m_tr_x = track->get_x();
231  m_tr_y = track->get_y();
232  m_tr_z = track->get_z();
233 
234  for (const std::string calo_name : _calo_names)
235  {
236  std::string detector(calo_name);
237  int detector_layer = -1;
238 
239  // special treatment for layered HCAL
240  if (calo_name.find("LFHCAL") == 0)
241  {
242  detector = "LFHCAL";
243  detector_layer = std::stoi(
244  calo_name.substr(
245  calo_name.find("_") + 1));
246 
247  static set<string> known_calo_name;
248  if (known_calo_name.find(calo_name) == known_calo_name.end())
249  {
250  known_calo_name.insert(calo_name);
251 
252  cout << __PRETTY_FUNCTION__ << " : Found new LFHCAL " << calo_name << " which is parsed as "
253  << "detector = " << detector << ", detector_layer = " << detector_layer << endl;
254  }
255  }
256 
257  std::string towernodename = "TOWER_CALIB_" + detector;
258  // Grab the towers
259  RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode,
260  towernodename);
261  if (!towers)
262  {
263  std::cout << PHWHERE << ": Could not find node " << towernodename
264  << std::endl;
265  return Fun4AllReturnCodes::ABORTRUN;
266  }
267  std::string towergeomnodename = "TOWERGEOM_" + detector;
268  RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(
269  topNode, towergeomnodename);
270  if (!towergeom)
271  {
272  std::cout << PHWHERE << ": Could not find node " << towergeomnodename
273  << std::endl;
274  return Fun4AllReturnCodes::ABORTRUN;
275  }
276 
277  bool has_projection(false);
278  // find projections
279  for (SvtxTrack::ConstStateIter trkstates = track->begin_states();
280  trkstates != track->end_states();
281  ++trkstates)
282  {
283  if (Verbosity() > 1) cout << __PRETTY_FUNCTION__ << " checking " << trkstates->second->get_name() << endl;
284  if (trkstates->second->get_name().compare(detector) == 0)
285  {
286  if (Verbosity() > 1) cout << __PRETTY_FUNCTION__ << " found " << trkstates->second->get_name() << endl;
287  has_projection = true;
288 
289  // setting the projection (xyz and pxpypz)
290  for (int i = 0; i < 3; i++)
291  {
292  m_CaloDataMap[calo_name].m_TTree_proj_vec[i] = trkstates->second->get_pos(i);
293  m_CaloDataMap[calo_name].m_TTree_proj_p_vec[i] = trkstates->second->get_mom(i);
294  }
295  // fourth element is the path length
296  m_CaloDataMap[calo_name].m_TTree_proj_vec[3] = trkstates->first;
297  }
298  } // for (SvtxTrack::ConstStateIter trkstates = track->begin_states();
299 
300  // projection match to calorimeter
301  if (has_projection)
302  {
303  TVector3 vec_proj(m_CaloDataMap[calo_name].m_TTree_proj_vec[0],
304  m_CaloDataMap[calo_name].m_TTree_proj_vec[1],
305  m_CaloDataMap[calo_name].m_TTree_proj_vec[2]);
306 
307  const double eta_proj = vec_proj.Eta();
308  const double phi_proj = vec_proj.Phi();
309 
310  double min_tower_r2 = 1000;
311  RawTowerDefs::keytype central_tower_key = -1;
312  const RawTowerGeom *central_tower(nullptr);
313  int maxBinPhi = 0;
314  int minBinPhi = 10000;
315 
316  RawTowerGeomContainer::ConstRange range = towergeom->get_tower_geometries();
318  for (RawTowerGeomContainer::ConstIterator titer = range.first;
319  titer != range.second;
320  ++titer)
321  {
322  const RawTowerGeom *tower_geom = titer->second;
323  assert(tower_geom);
324 
325  // special layer selection for LFHCAL
326  if (detector_layer >= 0)
327  {
328  const int layer = tower_geom->get_binl();
329  if (layer != detector_layer) continue;
330  }
331 
332  // const int bineta = RawTowerDefs::decode_index1(titer->first);
333  const int binphi = RawTowerDefs::decode_index2(titer->first);
334  if (maxBinPhi < binphi) maxBinPhi = binphi;
335  if (minBinPhi > binphi) minBinPhi = binphi;
336 
337  TVector3 vec_tower(
338  tower_geom->get_center_x(),
339  tower_geom->get_center_y(),
340  tower_geom->get_center_z());
341 
342  const double deta = eta_proj - vec_tower.Eta();
343  double dphi = phi_proj - vec_tower.Phi();
344 
345  if (dphi > M_PI) dphi -= 2 * M_PI;
346  if (dphi < -M_PI) dphi += 2 * M_PI;
347 
348  const double r2 = dphi * dphi + deta * deta;
349 
350  if (r2 < min_tower_r2)
351  {
352  min_tower_r2 = r2;
353  central_tower_key = titer->first;
354  central_tower = titer->second;
355  }
356 
357  } // for (RawTowerGeomContainer::ConstIterator titer = range.first;
358 
359  if (central_tower == nullptr) continue;
360  if (central_tower_key < 0) continue;
361 
362  if (Verbosity() > 1)
363  {
364  cout << __PRETTY_FUNCTION__ << " / " << calo_name << " found central tower " << central_tower_key << ": ";
365  cout << " min_tower_r2 = " << min_tower_r2;
366  cout << " decode_index1 = " << RawTowerDefs::decode_index1(central_tower_key);
367  cout << " decode_index2 = " << RawTowerDefs::decode_index2(central_tower_key);
368  cout << " decode_index1v2 = " << RawTowerDefs::decode_index1v2(central_tower_key);
369  cout << " decode_index2v2 = " << RawTowerDefs::decode_index2v2(central_tower_key);
370  cout << " decode_index3v2 = " << RawTowerDefs::decode_index3v2(central_tower_key);
371  cout << " minBinPhi = " << minBinPhi;
372  cout << " maxBinPhi = " << maxBinPhi;
373  central_tower->identify();
374  }
375 
376  // print tower patch
377  // const int central_tower_eta = RawTowerDefs::decode_index1(central_tower_key);
378  // const int central_tower_phi = RawTowerDefs::decode_index2(central_tower_key);
379  const int central_tower_eta = central_tower->get_bineta();
380  const int central_tower_phi = central_tower->get_binphi();
381  const int size_half_tower_patch = (m_CaloDataMap[calo_name].m_sizeTowerPatch - 1) / 2; // 7x7
382  size_t tower_index_patch = 0;
383 
384  m_CaloDataMap[calo_name].m_centralTowerBinEta = central_tower_eta;
385  m_CaloDataMap[calo_name].m_centralTowerBinPhi = central_tower_phi;
386 
387  for (int ieta_patch = -size_half_tower_patch;
388  ieta_patch <= +size_half_tower_patch;
389  ++ieta_patch)
390  {
391  const int bin_eta = central_tower_eta + ieta_patch;
392 
393  for (int iphi_patch = -size_half_tower_patch;
394  iphi_patch <= +size_half_tower_patch;
395  ++iphi_patch)
396  {
397  assert(tower_index_patch < m_CaloDataMap[calo_name].m_TTree_Tower_E.size());
398 
399  int bin_phi = central_tower_phi + iphi_patch;
400 
401  if (bin_phi < minBinPhi) bin_phi = bin_phi - minBinPhi + maxBinPhi + 1;
402  if (bin_phi > maxBinPhi) bin_phi = bin_phi - maxBinPhi + minBinPhi - 1;
403 
404  m_CaloDataMap[calo_name].m_TTree_Tower_iEta_patch[tower_index_patch] = ieta_patch;
405  m_CaloDataMap[calo_name].m_TTree_Tower_iPhi_patch[tower_index_patch] = iphi_patch;
406 
407  // compose tower key:
408  RawTowerDefs::keytype tower_key = 0;
409 
410  if (detector_layer == -1)
411  {
412  if (bin_eta >= 0xFFF or bin_phi >= 0xFFF or bin_eta < 0 or bin_phi < 0)
413  {
414  if (Verbosity())
415  {
416  cout << __PRETTY_FUNCTION__ << " / " << calo_name << " invalid tower geom " << central_tower_key << ": ";
417  cout << " bin_eta = " << bin_eta;
418  cout << " bin_phi = " << bin_phi;
419  cout << " central_tower_eta = " << central_tower_eta;
420  cout << " central_tower_phi = " << central_tower_phi;
421  cout << " central_tower_key = " << central_tower_key;
422  cout << " min_tower_r2 = " << min_tower_r2;
423  cout << " decode_index1 = " << RawTowerDefs::decode_index1(central_tower_key);
424  cout << " decode_index2 = " << RawTowerDefs::decode_index2(central_tower_key);
425  cout << " minBinPhi = " << minBinPhi;
426  cout << " maxBinPhi = " << maxBinPhi;
427  central_tower->identify();
428  }
429 
430  return Fun4AllReturnCodes::ABORTEVENT;
431  // continue;
432  }
433 
434  tower_key = RawTowerDefs::encode_towerid(
435  towergeom->get_calorimeter_id(), bin_eta, bin_phi);
436  }
437  else // (detector_layer > -1)
438  {
439  // LFHCAL special
440 
441  if (bin_eta >= 0x3FF or bin_phi >= 0x3FF or bin_eta < 0 or bin_phi < 0)
442  {
443  if (Verbosity())
444  {
445  cout << __PRETTY_FUNCTION__ << " / " << calo_name << " invalid tower geom " << central_tower_key << ": ";
446  cout << " bin_eta = " << bin_eta;
447  cout << " bin_phi = " << bin_phi;
448  cout << " central_tower_eta = " << central_tower_eta;
449  cout << " central_tower_phi = " << central_tower_phi;
450  cout << " central_tower_key = " << central_tower_key;
451  cout << " min_tower_r2 = " << min_tower_r2;
452  cout << " decode_index1 = " << RawTowerDefs::decode_index1(central_tower_key);
453  cout << " decode_index2 = " << RawTowerDefs::decode_index2(central_tower_key);
454  cout << " minBinPhi = " << minBinPhi;
455  cout << " maxBinPhi = " << maxBinPhi;
456  central_tower->identify();
457  }
458 
459  return Fun4AllReturnCodes::ABORTEVENT;
460  // continue;
461  }
462 
463  tower_key = RawTowerDefs::encode_towerid(
464  towergeom->get_calorimeter_id(), bin_eta, bin_phi, detector_layer);
465  } //else // (detector_layer > -1)
466 
467  const RawTowerGeom *tower_geom = towergeom->get_tower_geometry(tower_key);
468  if (tower_geom)
469  {
470  TVector3 vec_tower(
471  tower_geom->get_center_x(),
472  tower_geom->get_center_y(),
473  tower_geom->get_center_z());
474 
475  const double deta = eta_proj - vec_tower.Eta();
476  double dphi = phi_proj - vec_tower.Phi();
477 
478  m_CaloDataMap[calo_name].m_TTree_Tower_dEta[tower_index_patch] = deta;
479  m_CaloDataMap[calo_name].m_TTree_Tower_dPhi[tower_index_patch] = dphi;
480 
481  if (Verbosity() > 2)
482  {
483  cout << __PRETTY_FUNCTION__ << " / " << calo_name << " process tower geom " << tower_key << ": ";
484  cout << " ieta_patch = " << ieta_patch;
485  cout << " iphi_patch = " << iphi_patch;
486  cout << " bin_eta = " << bin_eta;
487  cout << " bin_phi = " << bin_phi;
488  cout << " deta = " << deta;
489  cout << " dphi = " << dphi;
490  tower_geom->identify();
491  }
492 
493  const RawTower *tower = towers->getTower(tower_key);
494 
495  if (tower)
496  {
497  const double energy = tower->get_energy();
498 
499  if (abs(ieta_patch) <= 1 and abs(iphi_patch) <= 1) m_CaloDataMap[calo_name].m_E3x3 += energy;
500  if (abs(ieta_patch) <= 2 and abs(iphi_patch) <= 2) m_CaloDataMap[calo_name].m_E5x5 += energy;
501  if (abs(ieta_patch) <= 3 and abs(iphi_patch) <= 3) m_CaloDataMap[calo_name].m_E7x7 += energy;
502 
503  m_CaloDataMap[calo_name].m_TTree_Tower_E[tower_index_patch] = energy;
504 
505  if (Verbosity() > 2)
506  {
507  cout << __PRETTY_FUNCTION__ << " / " << calo_name << " process tower " << tower_key << ": ";
508  cout << " ieta_patch = " << ieta_patch;
509  cout << " iphi_patch = " << iphi_patch;
510  cout << " bin_eta = " << bin_eta;
511  cout << " bin_phi = " << bin_phi;
512  cout << " energy = " << energy;
513  tower->identify();
514  }
515 
516  } // if (tower)
517 
518  ++tower_index_patch;
519  } // if (tower_geom)
520 
521  } // for (int iphi_patch = central_tower_phi - size_half_tower_patch;
522 
523  } // for (int ieta_patch = central_tower_eta - size_half_tower_patch;
524 
525  // const int bineta = RawTowerDefs::decode_index1(central_tower_key);
526  // const int binphi = RawTowerDefs::decode_index2(central_tower_key);
527 
528  } // if (has_projection)
529 
530  } // iterate detectors
531 
532  } // if (track)
533 
535  assert(m_truthtree);
536  m_truthtree->Fill();
537  }
538 
539  return Fun4AllReturnCodes::EVENT_OK;
540 }
541 
546 int eIDMLInterface::End(PHCompositeNode *topNode)
547 {
548  if (Verbosity() > 1)
549  {
550  cout << "Ending eIDMLInterface analysis package" << endl;
551  }
552 
554  m_outfile->cd();
555  m_truthtree->Write();
557  m_phi_h->Write();
558  m_eta_phi_h->Write();
559 
561  m_outfile->Write();
562  m_outfile->Close();
563 
565  delete m_outfile;
566 
567  if (Verbosity() > 1)
568  {
569  cout << "Finished eIDMLInterface analysis package" << endl;
570  }
571 
572  return 0;
573 }
574 
581 void eIDMLInterface::getHEPMCTruth(PHCompositeNode *topNode)
582 {
584  PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
585 
587  if (!hepmceventmap)
588  {
589  cout << PHWHERE
590  << "HEPMC event map node is missing, can't collected HEPMC truth particles"
591  << endl;
592  return;
593  }
594 
596  if (Verbosity() > 1)
597  {
598  cout << "Getting HEPMC truth particles " << endl;
599  }
600 
603  for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
604  eventIter != hepmceventmap->end();
605  ++eventIter)
606  {
608  PHHepMCGenEvent *hepmcevent = eventIter->second;
609 
610  if (hepmcevent)
611  {
613  HepMC::GenEvent *truthevent = hepmcevent->getEvent();
614  if (!truthevent)
615  {
616  cout << PHWHERE
617  << "no evt pointer under phhepmvgeneventmap found "
618  << endl;
619  return;
620  }
621 
623  HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
624 
626  m_partid1 = pdfinfo->id1();
627  m_partid2 = pdfinfo->id2();
628  m_x1 = pdfinfo->x1();
629  m_x2 = pdfinfo->x2();
630 
632  m_mpi = truthevent->mpi();
633 
635  m_process_id = truthevent->signal_process_id();
636 
637  if (Verbosity() > 2)
638  {
639  cout << " Iterating over an event" << endl;
640  }
642  for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
643  iter != truthevent->particles_end();
644  ++iter)
645  {
647  m_truthenergy = (*iter)->momentum().e();
648  m_truthpid = (*iter)->pdg_id();
649 
650  m_trutheta = (*iter)->momentum().pseudoRapidity();
651  m_truthphi = (*iter)->momentum().phi();
652  m_truthpx = (*iter)->momentum().px();
653  m_truthpy = (*iter)->momentum().py();
654  m_truthpz = (*iter)->momentum().pz();
655  m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
656 
658  m_hepmctree->Fill();
659  m_numparticlesinevent++;
660  }
661  }
662  }
663 }
664 
671 void eIDMLInterface::getPHG4Truth(PHCompositeNode *topNode)
672 {
674  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
675 
676  if (!truthinfo)
677  {
678  cout << PHWHERE
679  << "PHG4TruthInfoContainer node is missing, can't collect G4 truth particles"
680  << endl;
681  return;
682  }
683 
685  PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
686 
688  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
689  iter != range.second;
690  ++iter)
691  {
693  const PHG4Particle *truth = iter->second;
694 
696  m_truthpx = truth->get_px();
697  m_truthpy = truth->get_py();
698  m_truthpz = truth->get_pz();
699  m_truthp = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy + m_truthpz * m_truthpz);
700  m_truthenergy = truth->get_e();
701 
702  m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
703 
704  m_truthphi = atan2(m_truthpy, m_truthpx);
705 
706  m_trutheta = atanh(m_truthpz / m_truthenergy);
708  if (m_trutheta != m_trutheta)
709  m_trutheta = -99;
710  m_truthpid = truth->get_pid();
711 
713  m_truthtree->Fill();
714  }
715 }
716 
721 void eIDMLInterface::getTracks(PHCompositeNode *topNode)
722 {
724  SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "TrackMap");
725  EICPIDParticleContainer *pidcontainer = findNode::getClass<EICPIDParticleContainer>(topNode, "EICPIDParticleMap");
726 
727  if (Verbosity() > 1 and pidcontainer == nullptr)
728  {
729  cout << "EICPIDParticleContainer named EICPIDParticleMap does not exist. Skip saving PID info" << endl;
730  }
731 
732  if (!trackmap)
733  {
734  cout << PHWHERE
735  << "TrackMap node is missing, can't collect tracks"
736  << endl;
737  return;
738  }
739 
741  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
742 
743  if (Verbosity() > 1)
744  {
745  cout << "Get the tracks" << endl;
746  }
747  for (SvtxTrackMap::Iter iter = trackmap->begin();
748  iter != trackmap->end();
749  ++iter)
750  {
751  SvtxTrack *track = iter->second;
752 
754  m_tr_px = track->get_px();
755  m_tr_py = track->get_py();
756  m_tr_pz = track->get_pz();
757  m_tr_p = sqrt(m_tr_px * m_tr_px + m_tr_py * m_tr_py + m_tr_pz * m_tr_pz);
758 
759  m_tr_pt = sqrt(m_tr_px * m_tr_px + m_tr_py * m_tr_py);
760 
761  // Make some cuts on the track to clean up sample
762  if (m_tr_pt < 0.5)
763  continue;
764 
765  m_tr_phi = track->get_phi();
766  m_tr_eta = track->get_eta();
767 
768  m_charge = track->get_charge();
769  m_chisq = track->get_chisq();
770  m_ndf = track->get_ndf();
771  m_dca = track->get_dca();
772  m_tr_x = track->get_x();
773  m_tr_y = track->get_y();
774  m_tr_z = track->get_z();
775 
777  SvtxTrack_FastSim *temp = dynamic_cast<SvtxTrack_FastSim *>(iter->second);
778  if (!temp)
779  {
780  if (Verbosity() > 0)
781  std::cout << "Skipping non fast track sim object..." << std::endl;
782  continue;
783  }
784 
786  PHG4Particle *truthtrack = truthinfo->GetParticle(temp->get_truth_track_id());
787  if (truthtrack)
788  {
789  m_truth_is_primary = truthinfo->is_primary(truthtrack);
790 
791  m_truthtrackpx = truthtrack->get_px();
792  m_truthtrackpy = truthtrack->get_py();
793  m_truthtrackpz = truthtrack->get_pz();
794  m_truthtrackp = sqrt(m_truthtrackpx * m_truthtrackpx + m_truthtrackpy * m_truthtrackpy + m_truthtrackpz * m_truthtrackpz);
795 
796  m_truthtracke = truthtrack->get_e();
797 
798  m_truthtrackpt = sqrt(m_truthtrackpx * m_truthtrackpx + m_truthtrackpy * m_truthtrackpy);
799  m_truthtrackphi = atan2(m_truthtrackpy, m_truthtrackpx);
800  m_truthtracketa = atanh(m_truthtrackpz / m_truthtrackp);
801  m_truthtrackpid = truthtrack->get_pid();
802  }
803 
804  // match to PIDparticles
805  if (pidcontainer)
806  {
807  // EICPIDParticle are index the same as the tracks
808  const EICPIDParticle *pid_particle =
809  pidcontainer->findEICPIDParticle(track->get_id());
810 
811  if (pid_particle)
812  {
813  // top level log likelihood sums.
814  // More detailed per-detector information also available at EICPIDParticle::get_LogLikelyhood(EICPIDDefs::PIDCandidate, EICPIDDefs::PIDDetector)
815  m_tr_pion_loglikelihood = pid_particle->get_SumLogLikelyhood(EICPIDDefs::PionCandiate);
816  m_tr_kaon_loglikelihood = pid_particle->get_SumLogLikelyhood(EICPIDDefs::KaonCandiate);
817  m_tr_proton_loglikelihood = pid_particle->get_SumLogLikelyhood(EICPIDDefs::ProtonCandiate);
818  }
819  }
820 
821  m_tracktree->Fill();
822  }
823 }
824 
828 void eIDMLInterface::getTruthJets(PHCompositeNode *topNode)
829 {
830  if (Verbosity() > 1)
831  {
832  cout << "get the truth jets" << endl;
833  }
834 
836  JetMap *truth_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Truth_r05");
837 
839  JetMap *reco_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_r05");
840  if (!m_jetEvalStack)
841  {
842  m_jetEvalStack = new JetEvalStack(topNode, "AntiKt_Tower_r05",
843  "AntiKt_Truth_r05");
844  }
845  m_jetEvalStack->next_event(topNode);
846  JetTruthEval *trutheval = m_jetEvalStack->get_truth_eval();
847 
848  if (!truth_jets)
849  {
850  cout << PHWHERE
851  << "Truth jet node is missing, can't collect truth jets"
852  << endl;
853  return;
854  }
855 
857  for (JetMap::Iter iter = truth_jets->begin();
858  iter != truth_jets->end();
859  ++iter)
860  {
861  Jet *truthJet = iter->second;
862 
863  m_truthjetpt = truthJet->get_pt();
864 
865  std::set<PHG4Particle *> truthjetcomp =
866  trutheval->all_truth_particles(truthJet);
867  int ntruthconstituents = 0;
868  //loop over the constituents of the truth jet
869  for (std::set<PHG4Particle *>::iterator iter2 = truthjetcomp.begin();
870  iter2 != truthjetcomp.end();
871  ++iter2)
872  {
873  //get the particle of the truthjet
874  PHG4Particle *truthpart = *iter2;
875  if (!truthpart)
876  {
877  cout << "no truth particles in the jet??" << endl;
878  break;
879  }
880 
881  ntruthconstituents++;
882  }
883 
884  if (ntruthconstituents < 3)
885  continue;
887  if (m_truthjetpt < m_minjetpt)
888  continue;
889 
890  m_truthjeteta = truthJet->get_eta();
891  m_truthjetpx = truthJet->get_px();
892  m_truthjetpy = truthJet->get_py();
893  m_truthjetpz = truthJet->get_pz();
894  m_truthjetphi = truthJet->get_phi();
895  m_truthjetp = truthJet->get_p();
896  m_truthjetenergy = truthJet->get_e();
897 
898  m_recojetpt = 0;
899  m_recojetid = 0;
900  m_recojetpx = 0;
901  m_recojetpy = 0;
902  m_recojetpz = 0;
903  m_recojetphi = 0;
904  m_recojetp = 0;
905  m_recojetenergy = 0;
906  m_dR = -99;
907  float closestjet = 9999;
909  for (JetMap::Iter recoIter = reco_jets->begin();
910  recoIter != reco_jets->end();
911  ++recoIter)
912  {
913  const Jet *recoJet = recoIter->second;
914  m_recojetpt = recoJet->get_pt();
915  if (m_recojetpt < m_minjetpt - m_minjetpt * 0.4)
916  continue;
917 
918  m_recojeteta = recoJet->get_eta();
919  m_recojetphi = recoJet->get_phi();
920 
921  if (Verbosity() > 1)
922  {
923  cout << "matching by distance jet" << endl;
924  }
925 
926  float dphi = m_recojetphi - m_truthjetphi;
927  if (dphi > TMath::Pi())
928  dphi -= TMath::Pi() * 2.;
929  if (dphi < -1 * TMath::Pi())
930  dphi += TMath::Pi() * 2.;
931 
932  float deta = m_recojeteta - m_truthjeteta;
935  m_dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
936 
939  if (m_dR < truth_jets->get_par() && m_dR < closestjet)
940  {
941  // Get reco jet characteristics
942  m_recojetid = recoJet->get_id();
943  m_recojetpx = recoJet->get_px();
944  m_recojetpy = recoJet->get_py();
945  m_recojetpz = recoJet->get_pz();
946  m_recojetphi = recoJet->get_phi();
947  m_recojetp = recoJet->get_p();
948  m_recojetenergy = recoJet->get_e();
949  }
950  }
951 
953  m_truthjettree->Fill();
954  }
955 }
956 
960 void eIDMLInterface::getReconstructedJets(PHCompositeNode *topNode)
961 {
963  JetMap *reco_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_r05");
965  JetMap *truth_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Truth_r05");
966 
967  if (!m_jetEvalStack)
968  {
969  m_jetEvalStack = new JetEvalStack(topNode, "AntiKt_Tower_r05",
970  "AntiKt_Truth_r05");
971  }
972  m_jetEvalStack->next_event(topNode);
973  JetRecoEval *recoeval = m_jetEvalStack->get_reco_eval();
974  if (!reco_jets)
975  {
976  cout << PHWHERE
977  << "Reconstructed jet node is missing, can't collect reconstructed jets"
978  << endl;
979  return;
980  }
981 
982  if (Verbosity() > 1)
983  {
984  cout << "Get all Reco Jets" << endl;
985  }
986 
988  for (JetMap::Iter recoIter = reco_jets->begin();
989  recoIter != reco_jets->end();
990  ++recoIter)
991  {
992  Jet *recoJet = recoIter->second;
993  m_recojetpt = recoJet->get_pt();
994  if (m_recojetpt < m_minjetpt)
995  continue;
996 
997  m_recojeteta = recoJet->get_eta();
998 
999  // Get reco jet characteristics
1000  m_recojetid = recoJet->get_id();
1001  m_recojetpx = recoJet->get_px();
1002  m_recojetpy = recoJet->get_py();
1003  m_recojetpz = recoJet->get_pz();
1004  m_recojetphi = recoJet->get_phi();
1005  m_recojetp = recoJet->get_p();
1006  m_recojetenergy = recoJet->get_e();
1007 
1008  if (Verbosity() > 1)
1009  {
1010  cout << "matching by distance jet" << endl;
1011  }
1012 
1014  m_truthjetid = 0;
1015  m_truthjetp = 0;
1016  m_truthjetphi = 0;
1017  m_truthjeteta = 0;
1018  m_truthjetpt = 0;
1019  m_truthjetenergy = 0;
1020  m_truthjetpx = 0;
1021  m_truthjetpy = 0;
1022  m_truthjetpz = 0;
1023 
1024  Jet *truthjet = recoeval->max_truth_jet_by_energy(recoJet);
1025  if (truthjet)
1026  {
1027  m_truthjetid = truthjet->get_id();
1028  m_truthjetp = truthjet->get_p();
1029  m_truthjetpx = truthjet->get_px();
1030  m_truthjetpy = truthjet->get_py();
1031  m_truthjetpz = truthjet->get_pz();
1032  m_truthjeteta = truthjet->get_eta();
1033  m_truthjetphi = truthjet->get_phi();
1034  m_truthjetenergy = truthjet->get_e();
1035  m_truthjetpt = sqrt(m_truthjetpx * m_truthjetpx + m_truthjetpy * m_truthjetpy);
1036  }
1037 
1039  else if (truth_jets)
1040  {
1043  float closestjet = 9999;
1044  for (JetMap::Iter truthIter = truth_jets->begin();
1045  truthIter != truth_jets->end();
1046  ++truthIter)
1047  {
1048  const Jet *truthJet = truthIter->second;
1049 
1050  float thisjetpt = truthJet->get_pt();
1051  if (thisjetpt < m_minjetpt)
1052  continue;
1053 
1054  float thisjeteta = truthJet->get_eta();
1055  float thisjetphi = truthJet->get_phi();
1056 
1057  float dphi = m_recojetphi - thisjetphi;
1058  if (dphi > TMath::Pi())
1059  dphi -= TMath::Pi() * 2.;
1060  if (dphi < -1. * TMath::Pi())
1061  dphi += TMath::Pi() * 2.;
1062 
1063  float deta = m_recojeteta - thisjeteta;
1066  m_dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
1067 
1070  if (m_dR < reco_jets->get_par() && m_dR < closestjet)
1071  {
1072  m_truthjetid = -9999;
1073  m_truthjetp = truthJet->get_p();
1074  m_truthjetphi = truthJet->get_phi();
1075  m_truthjeteta = truthJet->get_eta();
1076  m_truthjetpt = truthJet->get_pt();
1077  m_truthjetenergy = truthJet->get_e();
1078  m_truthjetpx = truthJet->get_px();
1079  m_truthjetpy = truthJet->get_py();
1080  m_truthjetpz = truthJet->get_pz();
1081  closestjet = m_dR;
1082  }
1083  }
1084  }
1085  m_recojettree->Fill();
1086  }
1087 }
1088 
1095 void eIDMLInterface::getEMCalClusters(PHCompositeNode *topNode)
1096 {
1100  RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_BECAL");
1101 
1102  if (!clusters)
1103  {
1104  cout << PHWHERE
1105  << "EMCal cluster node is missing, can't collect EMCal clusters"
1106  << endl;
1107  return;
1108  }
1109 
1111  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
1112  if (!vertexmap)
1113  {
1114  cout << "eIDMLInterface::getEmcalClusters - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << endl;
1115  assert(vertexmap); // force quit
1116 
1117  return;
1118  }
1119 
1120  if (vertexmap->empty())
1121  {
1122  cout << "eIDMLInterface::getEmcalClusters - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << endl;
1123  return;
1124  }
1125 
1126  GlobalVertex *vtx = vertexmap->begin()->second;
1127  if (vtx == nullptr)
1128  return;
1129 
1131  CaloTriggerInfo *trigger = findNode::getClass<CaloTriggerInfo>(topNode, "CaloTriggerInfo");
1132 
1134  if (trigger)
1135  {
1136  m_E_4x4 = trigger->get_best_EMCal_4x4_E();
1137  }
1138  RawClusterContainer::ConstRange begin_end = clusters->getClusters();
1139  RawClusterContainer::ConstIterator clusIter;
1140 
1142  for (clusIter = begin_end.first;
1143  clusIter != begin_end.second;
1144  ++clusIter)
1145  {
1147  const RawCluster *cluster = clusIter->second;
1148 
1153  CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
1154  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster, vertex);
1155  m_clusenergy = E_vec_cluster.mag();
1156  m_cluseta = E_vec_cluster.pseudoRapidity();
1157  m_clustheta = E_vec_cluster.getTheta();
1158  m_cluspt = E_vec_cluster.perp();
1159  m_clusphi = E_vec_cluster.getPhi();
1160 
1161  m_phi_h->Fill(m_clusphi);
1162  m_eta_phi_h->Fill(m_clusphi, m_cluseta);
1163 
1164  if (m_cluspt < m_mincluspt)
1165  continue;
1166 
1167  m_cluspx = m_cluspt * cos(m_clusphi);
1168  m_cluspy = m_cluspt * sin(m_clusphi);
1169  m_cluspz = sqrt(m_clusenergy * m_clusenergy - m_cluspx * m_cluspx - m_cluspy * m_cluspy);
1170 
1171  //fill the cluster tree with all emcal clusters
1172  m_clustertree->Fill();
1173  }
1174 }
1175 
1180 void eIDMLInterface::initializeTrees()
1181 {
1182  m_outfile = new TFile();
1183  m_phi_h = new TH1F();
1184  m_eta_phi_h = new TH2F();
1185  // m_recojettree = new TTree("jettree", "A tree with reconstructed jets");
1186  // m_recojettree->Branch("m_recojetpt", &m_recojetpt, "m_recojetpt/D");
1187  // m_recojettree->Branch("m_recojetid", &m_recojetid, "m_recojetid/I");
1188  // m_recojettree->Branch("m_recojetpx", &m_recojetpx, "m_recojetpx/D");
1189  // m_recojettree->Branch("m_recojetpy", &m_recojetpy, "m_recojetpy/D");
1190  // m_recojettree->Branch("m_recojetpz", &m_recojetpz, "m_recojetpz/D");
1191  // m_recojettree->Branch("m_recojetphi", &m_recojetphi, "m_recojetphi/D");
1192  // m_recojettree->Branch("m_recojeteta", &m_recojeteta, "m_recojeteta/D");
1193  // m_recojettree->Branch("m_recojetenergy", &m_recojetenergy, "m_recojetenergy/D");
1194  // m_recojettree->Branch("m_truthjetid", &m_truthjetid, "m_truthjetid/I");
1195  // m_recojettree->Branch("m_truthjetp", &m_truthjetp, "m_truthjetp/D");
1196  // m_recojettree->Branch("m_truthjetphi", &m_truthjetphi, "m_truthjetphi/D");
1197  // m_recojettree->Branch("m_truthjeteta", &m_truthjeteta, "m_truthjeteta/D");
1198  // m_recojettree->Branch("m_truthjetpt", &m_truthjetpt, "m_truthjetpt/D");
1199  // m_recojettree->Branch("m_truthjetenergy", &m_truthjetenergy, "m_truthjetenergy/D");
1200  // m_recojettree->Branch("m_truthjetpx", &m_truthjetpx, "m_truthjetpx/D");
1201  // m_recojettree->Branch("m_truthjetpy", &m_truthjetpy, "m_truthjyetpy/D");
1202  // m_recojettree->Branch("m_truthjetpz", &m_truthjetpz, "m_truthjetpz/D");
1203  // m_recojettree->Branch("m_dR", &m_dR, "m_dR/D");
1204  //
1205  // m_truthjettree = new TTree("truthjettree", "A tree with truth jets");
1206  // m_truthjettree->Branch("m_truthjetid", &m_truthjetid, "m_truthjetid/I");
1207  // m_truthjettree->Branch("m_truthjetp", &m_truthjetp, "m_truthjetp/D");
1208  // m_truthjettree->Branch("m_truthjetphi", &m_truthjetphi, "m_truthjetphi/D");
1209  // m_truthjettree->Branch("m_truthjeteta", &m_truthjeteta, "m_truthjeteta/D");
1210  // m_truthjettree->Branch("m_truthjetpt", &m_truthjetpt, "m_truthjetpt/D");
1211  // m_truthjettree->Branch("m_truthjetenergy", &m_truthjetenergy, "m_truthjetenergy/D");
1212  // m_truthjettree->Branch("m_truthjetpx", &m_truthjetpx, "m_truthjetpx/D");
1213  // m_truthjettree->Branch("m_truthjetpy", &m_truthjetpy, "m_truthjetpy/D");
1214  // m_truthjettree->Branch("m_truthjetpz", &m_truthjetpz, "m_truthjetpz/D");
1215  // m_truthjettree->Branch("m_dR", &m_dR, "m_dR/D");
1216  // m_truthjettree->Branch("m_recojetpt", &m_recojetpt, "m_recojetpt/D");
1217  // m_truthjettree->Branch("m_recojetid", &m_recojetid, "m_recojetid/I");
1218  // m_truthjettree->Branch("m_recojetpx", &m_recojetpx, "m_recojetpx/D");
1219  // m_truthjettree->Branch("m_recojetpy", &m_recojetpy, "m_recojetpy/D");
1220  // m_truthjettree->Branch("m_recojetpz", &m_recojetpz, "m_recojetpz/D");
1221  // m_truthjettree->Branch("m_recojetphi", &m_recojetphi, "m_recojetphi/D");
1222  // m_truthjettree->Branch("m_recojeteta", &m_recojeteta, "m_recojeteta/D");
1223  // m_truthjettree->Branch("m_recojetenergy", &m_recojetenergy, "m_recojetenergy/D");
1224  //
1225  // m_tracktree = new TTree("tracktree", "A tree with svtx tracks");
1226  //
1227  // m_hepmctree = new TTree("hepmctree", "A tree with hepmc truth particles");
1228  // m_hepmctree->Branch("m_partid1", &m_partid1, "m_partid1/I");
1229  // m_hepmctree->Branch("m_partid2", &m_partid2, "m_partid2/I");
1230  // m_hepmctree->Branch("m_x1", &m_x1, "m_x1/D");
1231  // m_hepmctree->Branch("m_x2", &m_x2, "m_x2/D");
1232  // m_hepmctree->Branch("m_mpi", &m_mpi, "m_mpi/I");
1233  // m_hepmctree->Branch("m_process_id", &m_process_id, "m_process_id/I");
1234  // m_hepmctree->Branch("m_truthenergy", &m_truthenergy, "m_truthenergy/D");
1235  // m_hepmctree->Branch("m_trutheta", &m_trutheta, "m_trutheta/D");
1236  // m_hepmctree->Branch("m_truthphi", &m_truthphi, "m_truthphi/D");
1237  // m_hepmctree->Branch("m_truthpx", &m_truthpx, "m_truthpx/D");
1238  // m_hepmctree->Branch("m_truthpy", &m_truthpy, "m_truthpy/D");
1239  // m_hepmctree->Branch("m_truthpz", &m_truthpz, "m_truthpz/D");
1240  // m_hepmctree->Branch("m_truthpt", &m_truthpt, "m_truthpt/D");
1241  // m_hepmctree->Branch("m_numparticlesinevent", &m_numparticlesinevent, "m_numparticlesinevent/I");
1242  // m_hepmctree->Branch("m_truthpid", &m_truthpid, "m_truthpid/I");
1243 
1244  m_truthtree = new TTree("T", "A tree one enetry for a truth g4 particles matched with reco objects");
1245  m_truthtree->Branch("m_truthenergy", &m_truthenergy, "m_truthenergy/D");
1246  m_truthtree->Branch("m_truthp", &m_truthp, "m_truthp/D");
1247  m_truthtree->Branch("m_truthpx", &m_truthpx, "m_truthpx/D");
1248  m_truthtree->Branch("m_truthpy", &m_truthpy, "m_truthpy/D");
1249  m_truthtree->Branch("m_truthpz", &m_truthpz, "m_truthpz/D");
1250  m_truthtree->Branch("m_truthpt", &m_truthpt, "m_truthpt/D");
1251  m_truthtree->Branch("m_truthphi", &m_truthphi, "m_truthphi/D");
1252  m_truthtree->Branch("m_trutheta", &m_trutheta, "m_trutheta/D");
1253  m_truthtree->Branch("m_truthpid", &m_truthpid, "m_truthpid/I");
1254 
1255  m_truthtree->Branch("m_tr_px", &m_tr_px, "m_tr_px/D");
1256  m_truthtree->Branch("m_tr_py", &m_tr_py, "m_tr_py/D");
1257  m_truthtree->Branch("m_tr_pz", &m_tr_pz, "m_tr_pz/D");
1258  m_truthtree->Branch("m_tr_p", &m_tr_p, "m_tr_p/D");
1259  m_truthtree->Branch("m_tr_pt", &m_tr_pt, "m_tr_pt/D");
1260  m_truthtree->Branch("m_tr_phi", &m_tr_phi, "m_tr_phi/D");
1261  m_truthtree->Branch("m_tr_eta", &m_tr_eta, "m_tr_eta/D");
1262  m_truthtree->Branch("m_charge", &m_charge, "m_charge/I");
1263  m_truthtree->Branch("m_chisq", &m_chisq, "m_chisq/D");
1264  m_truthtree->Branch("m_ndf", &m_ndf, "m_ndf/I");
1265  m_truthtree->Branch("m_dca", &m_dca, "m_dca/D");
1266  m_truthtree->Branch("m_tr_x", &m_tr_x, "m_tr_x/D");
1267  m_truthtree->Branch("m_tr_y", &m_tr_y, "m_tr_y/D");
1268  m_truthtree->Branch("m_tr_z", &m_tr_z, "m_tr_z/D");
1269 
1270  const string xyzt[] = {"x", "y", "z", "t"};
1271 
1272  for (std::string calo_name : _calo_names)
1273  {
1274  for (int i = 0; i < 4; i++)
1275  {
1276  string bname = calo_name + "_proj_" + xyzt[i];
1277  string bdef = bname + "/F";
1278 
1279  // fourth element is the path length
1280  if (i == 3)
1281  {
1282  bdef = calo_name + "_proj_path_length" + "/F";
1283  }
1284 
1285  m_truthtree->Branch(bname.c_str(), &m_CaloDataMap[calo_name].m_TTree_proj_vec[i], bdef.c_str());
1286  }
1287 
1288  for (int i = 0; i < 3; i++)
1289  {
1290  string bname = calo_name + "_proj_p" + xyzt[i];
1291  string bdef = bname + "/F";
1292  m_truthtree->Branch(bname.c_str(), &m_CaloDataMap[calo_name].m_TTree_proj_p_vec[i], bdef.c_str());
1293  }
1294 
1295  // static const int nTowerInPatch = m_sizeTowerPatch * m_sizeTowerPatch;
1296  m_truthtree->Branch((calo_name + "_Tower_E3x3").c_str(), &m_CaloDataMap[calo_name].m_E3x3, (calo_name + "_Tower_E3x3/F").c_str());
1297  m_truthtree->Branch((calo_name + "_Tower_E5x5").c_str(), &m_CaloDataMap[calo_name].m_E5x5, (calo_name + "_Tower_E5x5/F").c_str());
1298  m_truthtree->Branch((calo_name + "_Tower_E7x7").c_str(), &m_CaloDataMap[calo_name].m_E7x7, (calo_name + "_Tower_E7x7/F").c_str());
1299  m_truthtree->Branch((calo_name + "_centralTowerBinEta").c_str(), &m_CaloDataMap[calo_name].m_centralTowerBinEta, (calo_name + "_centralTowerBinEta/I").c_str());
1300  m_truthtree->Branch((calo_name + "_centralTowerBinPhi").c_str(), &m_CaloDataMap[calo_name].m_centralTowerBinPhi, (calo_name + "_centralTowerBinPhi/I").c_str());
1301  m_truthtree->Branch((calo_name + "_nTowerInPatch").c_str(), &m_CaloDataMap[calo_name].nTowerInPatch, (calo_name + "_nTowerInPatch/I").c_str());
1302  m_truthtree->Branch((calo_name + "_Tower_dEta").c_str(), m_CaloDataMap[calo_name].m_TTree_Tower_dEta.data(), (calo_name + "_Tower_dEta[" + calo_name + "_nTowerInPatch]/F").c_str());
1303  m_truthtree->Branch((calo_name + "_Tower_dPhi").c_str(), m_CaloDataMap[calo_name].m_TTree_Tower_dPhi.data(), (calo_name + "_Tower_dPhi[" + calo_name + "_nTowerInPatch]/F").c_str());
1304  m_truthtree->Branch((calo_name + "_Tower_iEta_patch").c_str(), m_CaloDataMap[calo_name].m_TTree_Tower_iEta_patch.data(), (calo_name + "_Tower_iEta_patch[" + calo_name + "_nTowerInPatch]/I").c_str());
1305  m_truthtree->Branch((calo_name + "_Tower_iPhi_patch").c_str(), m_CaloDataMap[calo_name].m_TTree_Tower_iPhi_patch.data(), (calo_name + "_Tower_iPhi_patch[" + calo_name + "_nTowerInPatch]/I").c_str());
1306  m_truthtree->Branch((calo_name + "_Tower_E").c_str(), m_CaloDataMap[calo_name].m_TTree_Tower_E.data(), (calo_name + "_Tower_E[" + calo_name + "_nTowerInPatch]/F").c_str());
1307 
1308  // m_clustertree = new TTree("clustertree", "A tree with emcal clusters");
1309  // m_clustertree->Branch("m_clusenergy", &m_clusenergy, "m_clusenergy/D");
1310  // m_clustertree->Branch("m_cluseta", &m_cluseta, "m_cluseta/D");
1311  // m_clustertree->Branch("m_clustheta", &m_clustheta, "m_clustheta/D");
1312  // m_clustertree->Branch("m_cluspt", &m_cluspt, "m_cluspt/D");
1313  // m_clustertree->Branch("m_clusphi", &m_clusphi, "m_clusphi/D");
1314  // m_clustertree->Branch("m_cluspx", &m_cluspx, "m_cluspx/D");
1315  // m_clustertree->Branch("m_cluspy", &m_cluspy, "m_cluspy/D");
1316  // m_clustertree->Branch("m_cluspz", &m_cluspz, "m_cluspz/D");
1317  // m_clustertree->Branch("m_E_4x4", &m_E_4x4, "m_E_4x4/D");
1318  } // for (std::string calo_name : _calo_names)
1319 }
1320 
1321 void eIDMLInterface::CaloData::initializeVariables()
1322 {
1323  std::fill(m_TTree_proj_vec.begin(), m_TTree_proj_vec.end(), -9999);
1324  std::fill(m_TTree_proj_p_vec.begin(), m_TTree_proj_p_vec.end(), -9999);
1325 
1326  std::fill(m_TTree_Tower_dEta.begin(), m_TTree_Tower_dEta.end(), 0);
1327  std::fill(m_TTree_Tower_dPhi.begin(), m_TTree_Tower_dPhi.end(), 0);
1328  std::fill(m_TTree_Tower_iEta_patch.begin(), m_TTree_Tower_iEta_patch.end(), 0);
1329  std::fill(m_TTree_Tower_iPhi_patch.begin(), m_TTree_Tower_iPhi_patch.end(), 0);
1330  std::fill(m_TTree_Tower_E.begin(), m_TTree_Tower_E.end(), 0);
1331 
1332  m_centralTowerBinEta = -9999;
1333  m_centralTowerBinPhi = -9999;
1334  m_E3x3 = 0;
1335  m_E5x5 = 0;
1336  m_E7x7 = 0;
1337 }
1343 void eIDMLInterface::initializeVariables()
1344 {
1345  m_partid1 = -99;
1346  m_partid2 = -99;
1347  m_x1 = -99;
1348  m_x2 = -99;
1349  m_mpi = -99;
1350  m_process_id = -99;
1351  m_truthenergy = -99;
1352  m_trutheta = -99;
1353  m_truthphi = -99;
1354  m_truthp = -99;
1355  m_truthpx = -99;
1356  m_truthpy = -99;
1357  m_truthpz = -99;
1358  m_truthpt = -99;
1359  m_numparticlesinevent = -99;
1360  m_truthpid = -99;
1361 
1362  m_tr_px = -99;
1363  m_tr_py = -99;
1364  m_tr_pz = -99;
1365  m_tr_p = -99;
1366  m_tr_pt = -99;
1367  m_tr_phi = -99;
1368  m_tr_eta = -99;
1369  m_charge = -99;
1370  m_chisq = -99;
1371  m_ndf = -99;
1372  m_dca = -99;
1373  m_tr_x = -99;
1374  m_tr_y = -99;
1375  m_tr_z = -99;
1376  m_tr_pion_loglikelihood = -99;
1377  m_tr_kaon_loglikelihood = -99;
1378  m_tr_proton_loglikelihood = -99;
1379 
1380  m_truthtrackpy = -99;
1381  m_truthtrackpz = -99;
1382  m_truthtrackp = -99;
1383  m_truthtracke = -99;
1384  m_truthtrackpt = -99;
1385  m_truthtrackphi = -99;
1386  m_truthtracketa = -99;
1387  m_truthtrackpid = -99;
1388 
1389  m_recojetpt = -99;
1390  m_recojetid = -99;
1391  m_recojetpx = -99;
1392  m_recojetpy = -99;
1393  m_recojetpz = -99;
1394  m_recojetphi = -99;
1395  m_recojetp = -99;
1396  m_recojetenergy = -99;
1397  m_recojeteta = -99;
1398  m_truthjetid = -99;
1399  m_truthjetp = -99;
1400  m_truthjetphi = -99;
1401  m_truthjeteta = -99;
1402  m_truthjetpt = -99;
1403  m_truthjetenergy = -99;
1404  m_truthjetpx = -99;
1405  m_truthjetpy = -99;
1406  m_truthjetpz = -99;
1407  m_dR = -99;
1408 
1409  for (std::string calo_name : _calo_names)
1410  {
1411  m_CaloDataMap[calo_name].initializeVariables();
1412  }
1413 }