ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AnaTutorialECCE.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AnaTutorialECCE.cc
1 #include "AnaTutorialECCE.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 <g4vertex/GlobalVertex.h>
25 //PID includes
29 #include <g4eval/JetEvalStack.h>
30 #include <g4eval/SvtxEvalStack.h>
31 
33 #include <HepMC/GenEvent.h>
34 #include <HepMC/GenVertex.h>
37 
41 #include <g4main/PHG4Hit.h>
42 #include <g4main/PHG4Particle.h>
44 #include <phool/PHCompositeNode.h>
45 #include <phool/getClass.h>
46 
48 #include <TFile.h>
49 #include <TH1.h>
50 #include <TH2.h>
51 #include <TMath.h>
52 #include <TNtuple.h>
53 #include <TTree.h>
54 
56 #include <cassert>
57 #include <sstream>
58 #include <string>
59 
60 using namespace std;
61 
72 AnaTutorialECCE::AnaTutorialECCE(const std::string &name, const std::string &filename)
73  : SubsysReco(name)
74  , m_outfilename(filename)
75  , m_hm(nullptr)
76  , m_minjetpt(5.0)
77  , m_mincluspt(0.25)
78  , m_analyzeTracks(true)
79  , m_analyzeClusters(true)
80  , m_analyzeJets(true)
81  , m_analyzeTruth(false)
82 {
87 }
88 
93 {
94  delete m_hm;
95  delete m_hepmctree;
96  delete m_truthjettree;
97  delete m_recojettree;
98  delete m_tracktree;
99 }
100 
105 {
106  if (Verbosity() > 5)
107  {
108  cout << "Beginning Init in AnaTutorialECCE" << endl;
109  }
110 
111  m_outfile = new TFile(m_outfilename.c_str(), "RECREATE");
112 
113  m_phi_h = new TH1D("phi_h", ";Counts;#phi [rad]", 50, -6, 6);
114  m_eta_phi_h = new TH2F("phi_eta_h", ";#eta;#phi [rad]", 10, -1, 1, 50, -6, 6);
115 
116  return 0;
117 }
118 
124 {
125  if (Verbosity() > 5)
126  {
127  cout << "Beginning process_event in AnaTutorialECCE" << endl;
128  }
130  if (m_analyzeTruth)
131  {
132  getHEPMCTruth(topNode);
133  getPHG4Truth(topNode);
134  }
135 
137  if (m_analyzeTracks)
138  {
139  getTracks(topNode);
140  }
142  if (m_analyzeJets)
143  {
144  getTruthJets(topNode);
145  getReconstructedJets(topNode);
146  }
147 
149  if (m_analyzeClusters)
150  {
151  getEMCalClusters(topNode);
152  }
153 
155 }
156 
162 {
163  if (Verbosity() > 1)
164  {
165  cout << "Ending AnaTutorialECCE analysis package" << endl;
166  }
167 
169  m_outfile->cd();
170 
172  if (m_analyzeTracks)
173  m_tracktree->Write();
174 
176  if (m_analyzeJets)
177  {
178  m_truthjettree->Write();
179  m_recojettree->Write();
180  }
181 
183  if (m_analyzeTruth)
184  {
185  m_hepmctree->Write();
186  m_truthtree->Write();
187  }
188 
190  if (m_analyzeClusters)
191  {
192  m_clustertree->Write();
193  }
194 
196  m_phi_h->Write();
197  m_eta_phi_h->Write();
198 
200  m_outfile->Write();
201  m_outfile->Close();
202 
204  delete m_outfile;
205 
206  if (Verbosity() > 1)
207  {
208  cout << "Finished AnaTutorialECCE analysis package" << endl;
209  }
210 
211  return 0;
212 }
213 
221 {
223  PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
224 
226  if (!hepmceventmap)
227  {
228  cout << PHWHERE
229  << "HEPMC event map node is missing, can't collected HEPMC truth particles"
230  << endl;
231  return;
232  }
233 
235  if (Verbosity() > 1)
236  {
237  cout << "Getting HEPMC truth particles " << endl;
238  }
239 
242  for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
243  eventIter != hepmceventmap->end();
244  ++eventIter)
245  {
247  PHHepMCGenEvent *hepmcevent = eventIter->second;
248 
249  if (hepmcevent)
250  {
252  HepMC::GenEvent *truthevent = hepmcevent->getEvent();
253  if (!truthevent)
254  {
255  cout << PHWHERE
256  << "no evt pointer under phhepmvgeneventmap found "
257  << endl;
258  return;
259  }
260 
262  HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
263 
265  m_partid1 = pdfinfo->id1();
266  m_partid2 = pdfinfo->id2();
267  m_x1 = pdfinfo->x1();
268  m_x2 = pdfinfo->x2();
269 
271  m_mpi = truthevent->mpi();
272 
274  m_process_id = truthevent->signal_process_id();
275 
276  if (Verbosity() > 2)
277  {
278  cout << " Iterating over an event" << endl;
279  }
281  for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
282  iter != truthevent->particles_end();
283  ++iter)
284  {
286  m_truthenergy = (*iter)->momentum().e();
287  m_truthpid = (*iter)->pdg_id();
288 
289  m_trutheta = (*iter)->momentum().pseudoRapidity();
290  m_truthphi = (*iter)->momentum().phi();
291  m_truthpx = (*iter)->momentum().px();
292  m_truthpy = (*iter)->momentum().py();
293  m_truthpz = (*iter)->momentum().pz();
295 
297  m_hepmctree->Fill();
299  }
300  }
301  }
302 }
303 
311 {
313  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
314 
315  if (!truthinfo)
316  {
317  cout << PHWHERE
318  << "PHG4TruthInfoContainer node is missing, can't collect G4 truth particles"
319  << endl;
320  return;
321  }
322 
325 
327  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
328  iter != range.second;
329  ++iter)
330  {
332  const PHG4Particle *truth = iter->second;
333 
335  m_truthpx = truth->get_px();
336  m_truthpy = truth->get_py();
337  m_truthpz = truth->get_pz();
339  m_truthenergy = truth->get_e();
340 
342 
343  m_truthphi = atan2(m_truthpy, m_truthpx);
344 
345  m_trutheta = atanh(m_truthpz / m_truthenergy);
347  if (m_trutheta != m_trutheta)
348  m_trutheta = -99;
349  m_truthpid = truth->get_pid();
350 
352  m_truthtree->Fill();
353  }
354 }
355 
361 {
363  SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "TrackMap");
364  EICPIDParticleContainer *pidcontainer = findNode::getClass<EICPIDParticleContainer>(topNode, "EICPIDParticleMap");
365 
366  if (Verbosity() > 1 and pidcontainer == nullptr)
367  {
368  cout << "EICPIDParticleContainer named EICPIDParticleMap does not exist. Skip saving PID info" << endl;
369  }
370 
371  if (!trackmap)
372  {
373  cout << PHWHERE
374  << "TrackMap node is missing, can't collect tracks"
375  << endl;
376  return;
377  }
378 
380  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
381 
382  if (Verbosity() > 1)
383  {
384  cout << "Get the tracks" << endl;
385  }
386  for (SvtxTrackMap::Iter iter = trackmap->begin();
387  iter != trackmap->end();
388  ++iter)
389  {
390  SvtxTrack *track = iter->second;
391 
393  m_tr_px = track->get_px();
394  m_tr_py = track->get_py();
395  m_tr_pz = track->get_pz();
396  m_tr_p = sqrt(m_tr_px * m_tr_px + m_tr_py * m_tr_py + m_tr_pz * m_tr_pz);
397 
398  m_tr_pt = sqrt(m_tr_px * m_tr_px + m_tr_py * m_tr_py);
399 
400  // Make some cuts on the track to clean up sample
401  if (m_tr_pt < 0.5)
402  continue;
403 
404  m_tr_phi = track->get_phi();
405  m_tr_eta = track->get_eta();
406 
407  m_charge = track->get_charge();
408  m_chisq = track->get_chisq();
409  m_ndf = track->get_ndf();
410  m_dca = track->get_dca();
411  m_tr_x = track->get_x();
412  m_tr_y = track->get_y();
413  m_tr_z = track->get_z();
414 
416  SvtxTrack_FastSim *temp = dynamic_cast<SvtxTrack_FastSim *>(iter->second);
417  if (!temp)
418  {
419  if (Verbosity() > 0)
420  std::cout << "Skipping non fast track sim object..." << std::endl;
421  continue;
422  }
423 
425  PHG4Particle *truthtrack = truthinfo->GetParticle(temp->get_truth_track_id());
426  if (truthtrack)
427  {
428  m_truth_is_primary = truthinfo->is_primary(truthtrack);
429 
430  m_truthtrackpx = truthtrack->get_px();
431  m_truthtrackpy = truthtrack->get_py();
432  m_truthtrackpz = truthtrack->get_pz();
434 
435  m_truthtracke = truthtrack->get_e();
436 
438  m_truthtrackphi = atan2(m_truthtrackpy, m_truthtrackpx);
439  m_truthtracketa = atanh(m_truthtrackpz / m_truthtrackp);
440  m_truthtrackpid = truthtrack->get_pid();
441  }
442 
443  // match to PIDparticles
444  if (pidcontainer)
445  {
446  // EICPIDParticle are index the same as the tracks
447  const EICPIDParticle *pid_particle =
448  pidcontainer->findEICPIDParticle(track->get_id());
449 
450  if (pid_particle)
451  {
452  // top level log likelihood sums.
453  // More detailed per-detector information also available at EICPIDParticle::get_LogLikelyhood(EICPIDDefs::PIDCandidate, EICPIDDefs::PIDDetector)
457  }
458  }
459 
460  m_tracktree->Fill();
461  }
462 }
463 
468 {
469  if (Verbosity() > 1)
470  {
471  cout << "get the truth jets" << endl;
472  }
473 
475  JetMap *truth_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Truth_r05");
476 
478  JetMap *reco_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_r05");
479  if (!m_jetEvalStack)
480  {
481  m_jetEvalStack = new JetEvalStack(topNode, "AntiKt_Tower_r05",
482  "AntiKt_Truth_r05");
483  }
484  m_jetEvalStack->next_event(topNode);
486 
487  if (!truth_jets)
488  {
489  cout << PHWHERE
490  << "Truth jet node is missing, can't collect truth jets"
491  << endl;
492  return;
493  }
494 
496  for (JetMap::Iter iter = truth_jets->begin();
497  iter != truth_jets->end();
498  ++iter)
499  {
500  Jet *truthJet = iter->second;
501 
502  m_truthjetpt = truthJet->get_pt();
503 
504  std::set<PHG4Particle *> truthjetcomp =
505  trutheval->all_truth_particles(truthJet);
506  int ntruthconstituents = 0;
507  //loop over the constituents of the truth jet
508  for (std::set<PHG4Particle *>::iterator iter2 = truthjetcomp.begin();
509  iter2 != truthjetcomp.end();
510  ++iter2)
511  {
512  //get the particle of the truthjet
513  PHG4Particle *truthpart = *iter2;
514  if (!truthpart)
515  {
516  cout << "no truth particles in the jet??" << endl;
517  break;
518  }
519 
520  ntruthconstituents++;
521  }
522 
523  if (ntruthconstituents < 3)
524  continue;
526  if (m_truthjetpt < m_minjetpt)
527  continue;
528 
529  m_truthjeteta = truthJet->get_eta();
530  m_truthjetpx = truthJet->get_px();
531  m_truthjetpy = truthJet->get_py();
532  m_truthjetpz = truthJet->get_pz();
533  m_truthjetphi = truthJet->get_phi();
534  m_truthjetp = truthJet->get_p();
535  m_truthjetenergy = truthJet->get_e();
536 
537  m_recojetpt = 0;
538  m_recojetid = 0;
539  m_recojetpx = 0;
540  m_recojetpy = 0;
541  m_recojetpz = 0;
542  m_recojetphi = 0;
543  m_recojetp = 0;
544  m_recojetenergy = 0;
545  m_dR = -99;
546  float closestjet = 9999;
548  for (JetMap::Iter recoIter = reco_jets->begin();
549  recoIter != reco_jets->end();
550  ++recoIter)
551  {
552  const Jet *recoJet = recoIter->second;
553  m_recojetpt = recoJet->get_pt();
554  if (m_recojetpt < m_minjetpt - m_minjetpt * 0.4)
555  continue;
556 
557  m_recojeteta = recoJet->get_eta();
558  m_recojetphi = recoJet->get_phi();
559 
560  if (Verbosity() > 1)
561  {
562  cout << "matching by distance jet" << endl;
563  }
564 
565  float dphi = m_recojetphi - m_truthjetphi;
566  if (dphi > TMath::Pi())
567  dphi -= TMath::Pi() * 2.;
568  if (dphi < -1 * TMath::Pi())
569  dphi += TMath::Pi() * 2.;
570 
571  float deta = m_recojeteta - m_truthjeteta;
574  m_dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
575 
578  if (m_dR < truth_jets->get_par() && m_dR < closestjet)
579  {
580  // Get reco jet characteristics
581  m_recojetid = recoJet->get_id();
582  m_recojetpx = recoJet->get_px();
583  m_recojetpy = recoJet->get_py();
584  m_recojetpz = recoJet->get_pz();
585  m_recojetphi = recoJet->get_phi();
586  m_recojetp = recoJet->get_p();
587  m_recojetenergy = recoJet->get_e();
588  }
589  }
590 
592  m_truthjettree->Fill();
593  }
594 }
595 
600 {
602  JetMap *reco_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_r05");
604  JetMap *truth_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Truth_r05");
605 
606  if (!m_jetEvalStack)
607  {
608  m_jetEvalStack = new JetEvalStack(topNode, "AntiKt_Tower_r05",
609  "AntiKt_Truth_r05");
610  }
611  m_jetEvalStack->next_event(topNode);
613  if (!reco_jets)
614  {
615  cout << PHWHERE
616  << "Reconstructed jet node is missing, can't collect reconstructed jets"
617  << endl;
618  return;
619  }
620 
621  if (Verbosity() > 1)
622  {
623  cout << "Get all Reco Jets" << endl;
624  }
625 
627  for (JetMap::Iter recoIter = reco_jets->begin();
628  recoIter != reco_jets->end();
629  ++recoIter)
630  {
631  Jet *recoJet = recoIter->second;
632  m_recojetpt = recoJet->get_pt();
633  if (m_recojetpt < m_minjetpt)
634  continue;
635 
636  m_recojeteta = recoJet->get_eta();
637 
638  // Get reco jet characteristics
639  m_recojetid = recoJet->get_id();
640  m_recojetpx = recoJet->get_px();
641  m_recojetpy = recoJet->get_py();
642  m_recojetpz = recoJet->get_pz();
643  m_recojetphi = recoJet->get_phi();
644  m_recojetp = recoJet->get_p();
645  m_recojetenergy = recoJet->get_e();
646 
647  if (Verbosity() > 1)
648  {
649  cout << "matching by distance jet" << endl;
650  }
651 
653  m_truthjetid = 0;
654  m_truthjetp = 0;
655  m_truthjetphi = 0;
656  m_truthjeteta = 0;
657  m_truthjetpt = 0;
658  m_truthjetenergy = 0;
659  m_truthjetpx = 0;
660  m_truthjetpy = 0;
661  m_truthjetpz = 0;
662 
663  Jet *truthjet = recoeval->max_truth_jet_by_energy(recoJet);
664  if (truthjet)
665  {
666  m_truthjetid = truthjet->get_id();
667  m_truthjetp = truthjet->get_p();
668  m_truthjetpx = truthjet->get_px();
669  m_truthjetpy = truthjet->get_py();
670  m_truthjetpz = truthjet->get_pz();
671  m_truthjeteta = truthjet->get_eta();
672  m_truthjetphi = truthjet->get_phi();
673  m_truthjetenergy = truthjet->get_e();
675  }
676 
678  else if (truth_jets)
679  {
682  float closestjet = 9999;
683  for (JetMap::Iter truthIter = truth_jets->begin();
684  truthIter != truth_jets->end();
685  ++truthIter)
686  {
687  const Jet *truthJet = truthIter->second;
688 
689  float thisjetpt = truthJet->get_pt();
690  if (thisjetpt < m_minjetpt)
691  continue;
692 
693  float thisjeteta = truthJet->get_eta();
694  float thisjetphi = truthJet->get_phi();
695 
696  float dphi = m_recojetphi - thisjetphi;
697  if (dphi > TMath::Pi())
698  dphi -= TMath::Pi() * 2.;
699  if (dphi < -1. * TMath::Pi())
700  dphi += TMath::Pi() * 2.;
701 
702  float deta = m_recojeteta - thisjeteta;
705  m_dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
706 
709  if (m_dR < reco_jets->get_par() && m_dR < closestjet)
710  {
711  m_truthjetid = -9999;
712  m_truthjetp = truthJet->get_p();
713  m_truthjetphi = truthJet->get_phi();
714  m_truthjeteta = truthJet->get_eta();
715  m_truthjetpt = truthJet->get_pt();
716  m_truthjetenergy = truthJet->get_e();
717  m_truthjetpx = truthJet->get_px();
718  m_truthjetpy = truthJet->get_py();
719  m_truthjetpz = truthJet->get_pz();
720  closestjet = m_dR;
721  }
722  }
723  }
724  m_recojettree->Fill();
725  }
726 }
727 
735 {
739  RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_BECAL");
740 
741  if (!clusters)
742  {
743  cout << PHWHERE
744  << "EMCal cluster node is missing, can't collect EMCal clusters"
745  << endl;
746  return;
747  }
748 
750  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
751  if (!vertexmap)
752  {
753  cout << "AnaTutorialECCE::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;
754  assert(vertexmap); // force quit
755 
756  return;
757  }
758 
759  if (vertexmap->empty())
760  {
761  cout << "AnaTutorialECCE::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;
762  return;
763  }
764 
765  GlobalVertex *vtx = vertexmap->begin()->second;
766  if (vtx == nullptr)
767  return;
768 
770  CaloTriggerInfo *trigger = findNode::getClass<CaloTriggerInfo>(topNode, "CaloTriggerInfo");
771 
773  if (trigger)
774  {
775  m_E_4x4 = trigger->get_best_EMCal_4x4_E();
776  }
777  RawClusterContainer::ConstRange begin_end = clusters->getClusters();
779 
781  for (clusIter = begin_end.first;
782  clusIter != begin_end.second;
783  ++clusIter)
784  {
786  const RawCluster *cluster = clusIter->second;
787 
792  CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
793  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster, vertex);
794  m_clusenergy = E_vec_cluster.mag();
795  m_cluseta = E_vec_cluster.pseudoRapidity();
796  m_clustheta = E_vec_cluster.getTheta();
797  m_cluspt = E_vec_cluster.perp();
798  m_clusphi = E_vec_cluster.getPhi();
799 
800  m_phi_h->Fill(m_clusphi);
802 
803  if (m_cluspt < m_mincluspt)
804  continue;
805 
806  m_cluspx = m_cluspt * cos(m_clusphi);
807  m_cluspy = m_cluspt * sin(m_clusphi);
809 
810  //fill the cluster tree with all emcal clusters
811  m_clustertree->Fill();
812  }
813 }
814 
820 {
821  m_recojettree = new TTree("jettree", "A tree with reconstructed jets");
822  m_recojettree->Branch("m_recojetpt", &m_recojetpt, "m_recojetpt/D");
823  m_recojettree->Branch("m_recojetid", &m_recojetid, "m_recojetid/I");
824  m_recojettree->Branch("m_recojetpx", &m_recojetpx, "m_recojetpx/D");
825  m_recojettree->Branch("m_recojetpy", &m_recojetpy, "m_recojetpy/D");
826  m_recojettree->Branch("m_recojetpz", &m_recojetpz, "m_recojetpz/D");
827  m_recojettree->Branch("m_recojetphi", &m_recojetphi, "m_recojetphi/D");
828  m_recojettree->Branch("m_recojeteta", &m_recojeteta, "m_recojeteta/D");
829  m_recojettree->Branch("m_recojetenergy", &m_recojetenergy, "m_recojetenergy/D");
830  m_recojettree->Branch("m_truthjetid", &m_truthjetid, "m_truthjetid/I");
831  m_recojettree->Branch("m_truthjetp", &m_truthjetp, "m_truthjetp/D");
832  m_recojettree->Branch("m_truthjetphi", &m_truthjetphi, "m_truthjetphi/D");
833  m_recojettree->Branch("m_truthjeteta", &m_truthjeteta, "m_truthjeteta/D");
834  m_recojettree->Branch("m_truthjetpt", &m_truthjetpt, "m_truthjetpt/D");
835  m_recojettree->Branch("m_truthjetenergy", &m_truthjetenergy, "m_truthjetenergy/D");
836  m_recojettree->Branch("m_truthjetpx", &m_truthjetpx, "m_truthjetpx/D");
837  m_recojettree->Branch("m_truthjetpy", &m_truthjetpy, "m_truthjyetpy/D");
838  m_recojettree->Branch("m_truthjetpz", &m_truthjetpz, "m_truthjetpz/D");
839  m_recojettree->Branch("m_dR", &m_dR, "m_dR/D");
840 
841  m_truthjettree = new TTree("truthjettree", "A tree with truth jets");
842  m_truthjettree->Branch("m_truthjetid", &m_truthjetid, "m_truthjetid/I");
843  m_truthjettree->Branch("m_truthjetp", &m_truthjetp, "m_truthjetp/D");
844  m_truthjettree->Branch("m_truthjetphi", &m_truthjetphi, "m_truthjetphi/D");
845  m_truthjettree->Branch("m_truthjeteta", &m_truthjeteta, "m_truthjeteta/D");
846  m_truthjettree->Branch("m_truthjetpt", &m_truthjetpt, "m_truthjetpt/D");
847  m_truthjettree->Branch("m_truthjetenergy", &m_truthjetenergy, "m_truthjetenergy/D");
848  m_truthjettree->Branch("m_truthjetpx", &m_truthjetpx, "m_truthjetpx/D");
849  m_truthjettree->Branch("m_truthjetpy", &m_truthjetpy, "m_truthjetpy/D");
850  m_truthjettree->Branch("m_truthjetpz", &m_truthjetpz, "m_truthjetpz/D");
851  m_truthjettree->Branch("m_dR", &m_dR, "m_dR/D");
852  m_truthjettree->Branch("m_recojetpt", &m_recojetpt, "m_recojetpt/D");
853  m_truthjettree->Branch("m_recojetid", &m_recojetid, "m_recojetid/I");
854  m_truthjettree->Branch("m_recojetpx", &m_recojetpx, "m_recojetpx/D");
855  m_truthjettree->Branch("m_recojetpy", &m_recojetpy, "m_recojetpy/D");
856  m_truthjettree->Branch("m_recojetpz", &m_recojetpz, "m_recojetpz/D");
857  m_truthjettree->Branch("m_recojetphi", &m_recojetphi, "m_recojetphi/D");
858  m_truthjettree->Branch("m_recojeteta", &m_recojeteta, "m_recojeteta/D");
859  m_truthjettree->Branch("m_recojetenergy", &m_recojetenergy, "m_recojetenergy/D");
860  m_tracktree = new TTree("tracktree", "A tree with svtx tracks");
861  m_tracktree->Branch("m_tr_px", &m_tr_px, "m_tr_px/D");
862  m_tracktree->Branch("m_tr_py", &m_tr_py, "m_tr_py/D");
863  m_tracktree->Branch("m_tr_pz", &m_tr_pz, "m_tr_pz/D");
864  m_tracktree->Branch("m_tr_p", &m_tr_p, "m_tr_p/D");
865  m_tracktree->Branch("m_tr_pt", &m_tr_pt, "m_tr_pt/D");
866  m_tracktree->Branch("m_tr_phi", &m_tr_phi, "m_tr_phi/D");
867  m_tracktree->Branch("m_tr_eta", &m_tr_eta, "m_tr_eta/D");
868  m_tracktree->Branch("m_charge", &m_charge, "m_charge/I");
869  m_tracktree->Branch("m_chisq", &m_chisq, "m_chisq/D");
870  m_tracktree->Branch("m_ndf", &m_ndf, "m_ndf/I");
871  m_tracktree->Branch("m_dca", &m_dca, "m_dca/D");
872  m_tracktree->Branch("m_tr_x", &m_tr_x, "m_tr_x/D");
873  m_tracktree->Branch("m_tr_y", &m_tr_y, "m_tr_y/D");
874  m_tracktree->Branch("m_tr_z", &m_tr_z, "m_tr_z/D");
875  m_tracktree->Branch("m_tr_pion_loglikelihood", &m_tr_pion_loglikelihood, "m_tr_pion_loglikelihood/F");
876  m_tracktree->Branch("m_tr_kaon_loglikelihood", &m_tr_kaon_loglikelihood, "m_tr_kaon_loglikelihood/F");
877  m_tracktree->Branch("m_tr_proton_loglikelihood", &m_tr_proton_loglikelihood, "m_tr_proton_loglikelihood/F");
878  m_tracktree->Branch("m_truth_is_primary", &m_truth_is_primary, "m_truth_is_primary/I");
879  m_tracktree->Branch("m_truthtrackpx", &m_truthtrackpx, "m_truthtrackpx/D");
880  m_tracktree->Branch("m_truthtrackpy", &m_truthtrackpy, "m_truthtrackpy/D");
881  m_tracktree->Branch("m_truthtrackpz", &m_truthtrackpz, "m_truthtrackpz/D");
882  m_tracktree->Branch("m_truthtrackp", &m_truthtrackp, "m_truthtrackp/D");
883  m_tracktree->Branch("m_truthtracke", &m_truthtracke, "m_truthtracke/D");
884  m_tracktree->Branch("m_truthtrackpt", &m_truthtrackpt, "m_truthtrackpt/D");
885  m_tracktree->Branch("m_truthtrackphi", &m_truthtrackphi, "m_truthtrackphi/D");
886  m_tracktree->Branch("m_truthtracketa", &m_truthtracketa, "m_truthtracketa/D");
887  m_tracktree->Branch("m_truthtrackpid", &m_truthtrackpid, "m_truthtrackpid/I");
888 
889  m_hepmctree = new TTree("hepmctree", "A tree with hepmc truth particles");
890  m_hepmctree->Branch("m_partid1", &m_partid1, "m_partid1/I");
891  m_hepmctree->Branch("m_partid2", &m_partid2, "m_partid2/I");
892  m_hepmctree->Branch("m_x1", &m_x1, "m_x1/D");
893  m_hepmctree->Branch("m_x2", &m_x2, "m_x2/D");
894  m_hepmctree->Branch("m_mpi", &m_mpi, "m_mpi/I");
895  m_hepmctree->Branch("m_process_id", &m_process_id, "m_process_id/I");
896  m_hepmctree->Branch("m_truthenergy", &m_truthenergy, "m_truthenergy/D");
897  m_hepmctree->Branch("m_trutheta", &m_trutheta, "m_trutheta/D");
898  m_hepmctree->Branch("m_truthphi", &m_truthphi, "m_truthphi/D");
899  m_hepmctree->Branch("m_truthpx", &m_truthpx, "m_truthpx/D");
900  m_hepmctree->Branch("m_truthpy", &m_truthpy, "m_truthpy/D");
901  m_hepmctree->Branch("m_truthpz", &m_truthpz, "m_truthpz/D");
902  m_hepmctree->Branch("m_truthpt", &m_truthpt, "m_truthpt/D");
903  m_hepmctree->Branch("m_numparticlesinevent", &m_numparticlesinevent, "m_numparticlesinevent/I");
904  m_hepmctree->Branch("m_truthpid", &m_truthpid, "m_truthpid/I");
905 
906  m_truthtree = new TTree("truthg4tree", "A tree with truth g4 particles");
907  m_truthtree->Branch("m_truthenergy", &m_truthenergy, "m_truthenergy/D");
908  m_truthtree->Branch("m_truthp", &m_truthp, "m_truthp/D");
909  m_truthtree->Branch("m_truthpx", &m_truthpx, "m_truthpx/D");
910  m_truthtree->Branch("m_truthpy", &m_truthpy, "m_truthpy/D");
911  m_truthtree->Branch("m_truthpz", &m_truthpz, "m_truthpz/D");
912  m_truthtree->Branch("m_truthpt", &m_truthpt, "m_truthpt/D");
913  m_truthtree->Branch("m_truthphi", &m_truthphi, "m_truthphi/D");
914  m_truthtree->Branch("m_trutheta", &m_trutheta, "m_trutheta/D");
915  m_truthtree->Branch("m_truthpid", &m_truthpid, "m_truthpid/I");
916 
917  m_clustertree = new TTree("clustertree", "A tree with emcal clusters");
918  m_clustertree->Branch("m_clusenergy", &m_clusenergy, "m_clusenergy/D");
919  m_clustertree->Branch("m_cluseta", &m_cluseta, "m_cluseta/D");
920  m_clustertree->Branch("m_clustheta", &m_clustheta, "m_clustheta/D");
921  m_clustertree->Branch("m_cluspt", &m_cluspt, "m_cluspt/D");
922  m_clustertree->Branch("m_clusphi", &m_clusphi, "m_clusphi/D");
923  m_clustertree->Branch("m_cluspx", &m_cluspx, "m_cluspx/D");
924  m_clustertree->Branch("m_cluspy", &m_cluspy, "m_cluspy/D");
925  m_clustertree->Branch("m_cluspz", &m_cluspz, "m_cluspz/D");
926  m_clustertree->Branch("m_E_4x4", &m_E_4x4, "m_E_4x4/D");
927 }
928 
935 {
936  m_outfile = new TFile();
937  m_phi_h = new TH1F();
938  m_eta_phi_h = new TH2F();
939 
940  m_partid1 = -99;
941  m_partid2 = -99;
942  m_x1 = -99;
943  m_x2 = -99;
944  m_mpi = -99;
945  m_process_id = -99;
946  m_truthenergy = -99;
947  m_trutheta = -99;
948  m_truthphi = -99;
949  m_truthp = -99;
950  m_truthpx = -99;
951  m_truthpy = -99;
952  m_truthpz = -99;
953  m_truthpt = -99;
954  m_numparticlesinevent = -99;
955  m_truthpid = -99;
956 
957  m_tr_px = -99;
958  m_tr_py = -99;
959  m_tr_pz = -99;
960  m_tr_p = -99;
961  m_tr_pt = -99;
962  m_tr_phi = -99;
963  m_tr_eta = -99;
964  m_charge = -99;
965  m_chisq = -99;
966  m_ndf = -99;
967  m_dca = -99;
968  m_tr_x = -99;
969  m_tr_y = -99;
970  m_tr_z = -99;
974  m_truth_is_primary = -99;
975  m_truthtrackpx = -99;
976  m_truthtrackpy = -99;
977  m_truthtrackpz = -99;
978  m_truthtrackp = -99;
979  m_truthtracke = -99;
980  m_truthtrackpt = -99;
981  m_truthtrackphi = -99;
982  m_truthtracketa = -99;
983  m_truthtrackpid = -99;
984 
985  m_recojetpt = -99;
986  m_recojetid = -99;
987  m_recojetpx = -99;
988  m_recojetpy = -99;
989  m_recojetpz = -99;
990  m_recojetphi = -99;
991  m_recojetp = -99;
992  m_recojetenergy = -99;
993  m_recojeteta = -99;
994  m_truthjetid = -99;
995  m_truthjetp = -99;
996  m_truthjetphi = -99;
997  m_truthjeteta = -99;
998  m_truthjetpt = -99;
999  m_truthjetenergy = -99;
1000  m_truthjetpx = -99;
1001  m_truthjetpy = -99;
1002  m_truthjetpz = -99;
1003  m_dR = -99;
1004 }