ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AnaTutorial.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AnaTutorial.cc
1 #include "AnaTutorial.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>
24 
26 #include <g4eval/JetEvalStack.h>
27 #include <g4eval/SvtxEvalStack.h>
28 
30 #include <HepMC/GenEvent.h>
31 #include <HepMC/GenVertex.h>
34 
38 #include <g4main/PHG4Hit.h>
39 #include <g4main/PHG4Particle.h>
41 #include <phool/PHCompositeNode.h>
42 #include <phool/getClass.h>
43 
45 #include <TFile.h>
46 #include <TH1.h>
47 #include <TH2.h>
48 #include <TMath.h>
49 #include <TNtuple.h>
50 #include <TTree.h>
51 
53 #include <cassert>
54 #include <sstream>
55 #include <string>
56 
57 using namespace std;
58 
69 AnaTutorial::AnaTutorial(const std::string &name, const std::string &filename)
70  : SubsysReco(name)
71  , m_outfilename(filename)
72  , m_hm(nullptr)
73  , m_minjetpt(5.0)
74  , m_mincluspt(0.25)
75  , m_analyzeTracks(true)
76  , m_analyzeClusters(true)
77  , m_analyzeJets(true)
78  , m_analyzeTruth(false)
79 {
84 }
85 
90 {
91  delete m_hm;
92  delete m_hepmctree;
93  delete m_truthjettree;
94  delete m_recojettree;
95  delete m_tracktree;
96 }
97 
102 {
103  if (Verbosity() > 5)
104  {
105  cout << "Beginning Init in AnaTutorial" << endl;
106  }
107 
108  m_outfile = new TFile(m_outfilename.c_str(), "RECREATE");
109 
110  m_phi_h = new TH1D("phi_h", ";Counts;#phi [rad]", 50, -6, 6);
111  m_eta_phi_h = new TH2F("phi_eta_h", ";#eta;#phi [rad]", 10, -1, 1, 50, -6, 6);
112 
113  return 0;
114 }
115 
121 {
122  if (Verbosity() > 5)
123  {
124  cout << "Beginning process_event in AnaTutorial" << endl;
125  }
127  if (m_analyzeTruth)
128  {
129  getHEPMCTruth(topNode);
130  getPHG4Truth(topNode);
131  }
132 
134  if (m_analyzeTracks)
135  {
136  getTracks(topNode);
137  }
139  if (m_analyzeJets)
140  {
141  getTruthJets(topNode);
142  getReconstructedJets(topNode);
143  }
144 
146  if (m_analyzeClusters)
147  {
148  getEMCalClusters(topNode);
149  }
150 
152 }
153 
159 {
160  if (Verbosity() > 1)
161  {
162  cout << "Ending AnaTutorial analysis package" << endl;
163  }
164 
166  m_outfile->cd();
167 
169  if (m_analyzeTracks)
170  m_tracktree->Write();
171 
173  if (m_analyzeJets)
174  {
175  m_truthjettree->Write();
176  m_recojettree->Write();
177  }
178 
180  if (m_analyzeTruth)
181  {
182  m_hepmctree->Write();
183  m_truthtree->Write();
184  }
185 
187  if (m_analyzeClusters)
188  {
189  m_clustertree->Write();
190  }
191 
193  m_phi_h->Write();
194  m_eta_phi_h->Write();
195 
197  m_outfile->Write();
198  m_outfile->Close();
199 
201  delete m_outfile;
202 
203  if (Verbosity() > 1)
204  {
205  cout << "Finished AnaTutorial analysis package" << endl;
206  }
207 
208  return 0;
209 }
210 
218 {
220  PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
221 
223  if (!hepmceventmap)
224  {
225  cout << PHWHERE
226  << "HEPMC event map node is missing, can't collected HEPMC truth particles"
227  << endl;
228  return;
229  }
230 
232  if (Verbosity() > 1)
233  {
234  cout << "Getting HEPMC truth particles " << endl;
235  }
236 
239  for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
240  eventIter != hepmceventmap->end();
241  ++eventIter)
242  {
244  PHHepMCGenEvent *hepmcevent = eventIter->second;
245 
246  if (hepmcevent)
247  {
249  HepMC::GenEvent *truthevent = hepmcevent->getEvent();
250  if (!truthevent)
251  {
252  cout << PHWHERE
253  << "no evt pointer under phhepmvgeneventmap found "
254  << endl;
255  return;
256  }
257 
259  HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
260 
262  m_partid1 = pdfinfo->id1();
263  m_partid2 = pdfinfo->id2();
264  m_x1 = pdfinfo->x1();
265  m_x2 = pdfinfo->x2();
266 
268  m_mpi = truthevent->mpi();
269 
271  m_process_id = truthevent->signal_process_id();
272 
273  if (Verbosity() > 2)
274  {
275  cout << " Iterating over an event" << endl;
276  }
278  for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
279  iter != truthevent->particles_end();
280  ++iter)
281  {
283  m_truthenergy = (*iter)->momentum().e();
284  m_truthpid = (*iter)->pdg_id();
285 
286  m_trutheta = (*iter)->momentum().pseudoRapidity();
287  m_truthphi = (*iter)->momentum().phi();
288  m_truthpx = (*iter)->momentum().px();
289  m_truthpy = (*iter)->momentum().py();
290  m_truthpz = (*iter)->momentum().pz();
292 
294  m_hepmctree->Fill();
296  }
297  }
298  }
299 }
300 
308 {
310  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
311 
312  if (!truthinfo)
313  {
314  cout << PHWHERE
315  << "PHG4TruthInfoContainer node is missing, can't collect G4 truth particles"
316  << endl;
317  return;
318  }
319 
322 
324  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
325  iter != range.second;
326  ++iter)
327  {
329  const PHG4Particle *truth = iter->second;
330 
332  m_truthpx = truth->get_px();
333  m_truthpy = truth->get_py();
334  m_truthpz = truth->get_pz();
336  m_truthenergy = truth->get_e();
337 
339 
340  m_truthphi = atan(m_truthpy / m_truthpx);
341 
342  m_trutheta = atanh(m_truthpz / m_truthenergy);
344  if (m_trutheta != m_trutheta)
345  m_trutheta = -99;
346  m_truthpid = truth->get_pid();
347 
349  m_truthtree->Fill();
350  }
351 }
352 
359 {
361  SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
362 
363  if (!trackmap)
364  {
365  cout << PHWHERE
366  << "SvtxTrackMap node is missing, can't collect tracks"
367  << endl;
368  return;
369  }
370 
372  if(!m_svtxEvalStack)
373  {
374  m_svtxEvalStack = new SvtxEvalStack(topNode);
376  }
377 
378  m_svtxEvalStack->next_event(topNode);
379 
382 
384  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
385 
386  if (Verbosity() > 1)
387  {
388  cout << "Get the SVTX tracks" << endl;
389  }
390  for (SvtxTrackMap::Iter iter = trackmap->begin();
391  iter != trackmap->end();
392  ++iter)
393  {
394  SvtxTrack *track = iter->second;
395 
397  m_tr_px = track->get_px();
398  m_tr_py = track->get_py();
399  m_tr_pz = track->get_pz();
400  m_tr_p = sqrt(m_tr_px * m_tr_px + m_tr_py * m_tr_py + m_tr_pz * m_tr_pz);
401 
402  m_tr_pt = sqrt(m_tr_px * m_tr_px + m_tr_py * m_tr_py);
403 
404  // Make some cuts on the track to clean up sample
405  if (m_tr_pt < 0.5)
406  continue;
407 
408  m_tr_phi = track->get_phi();
409  m_tr_eta = track->get_eta();
410 
411  m_charge = track->get_charge();
412  m_chisq = track->get_chisq();
413  m_ndf = track->get_ndf();
414  m_dca = track->get_dca();
415  m_tr_x = track->get_x();
416  m_tr_y = track->get_y();
417  m_tr_z = track->get_z();
418 
420  PHG4Particle *truthtrack = trackeval->max_truth_particle_by_nclusters(track);
421  m_truth_is_primary = truthinfo->is_primary(truthtrack);
422 
423  m_truthtrackpx = truthtrack->get_px();
424  m_truthtrackpy = truthtrack->get_py();
425  m_truthtrackpz = truthtrack->get_pz();
427 
428  m_truthtracke = truthtrack->get_e();
429 
431  m_truthtrackphi = atan(m_truthtrackpy / m_truthtrackpx);
432  m_truthtracketa = atanh(m_truthtrackpz / m_truthtrackp);
433  m_truthtrackpid = truthtrack->get_pid();
434 
435  m_tracktree->Fill();
436  }
437 }
438 
443 {
444  if (Verbosity() > 1)
445  {
446  cout << "get the truth jets" << endl;
447  }
448 
450  JetMap *truth_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Truth_r04");
451 
453  JetMap *reco_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_r04");
454  if(!m_jetEvalStack)
455  {
456  m_jetEvalStack = new JetEvalStack(topNode, "AntiKt_Tower_r04",
457  "AntiKt_Truth_r04");
458  }
459  m_jetEvalStack->next_event(topNode);
461 
462  if (!truth_jets)
463  {
464  cout << PHWHERE
465  << "Truth jet node is missing, can't collect truth jets"
466  << endl;
467  return;
468  }
469 
471  for (JetMap::Iter iter = truth_jets->begin();
472  iter != truth_jets->end();
473  ++iter)
474  {
475  Jet *truthJet = iter->second;
476 
477  m_truthjetpt = truthJet->get_pt();
478 
479  std::set<PHG4Particle *> truthjetcomp =
480  trutheval->all_truth_particles(truthJet);
481  int ntruthconstituents = 0;
482  //loop over the constituents of the truth jet
483  for (std::set<PHG4Particle *>::iterator iter2 = truthjetcomp.begin();
484  iter2 != truthjetcomp.end();
485  ++iter2)
486  {
487  //get the particle of the truthjet
488  PHG4Particle *truthpart = *iter2;
489  if (!truthpart)
490  {
491  cout << "no truth particles in the jet??" << endl;
492  break;
493  }
494 
495  ntruthconstituents++;
496  }
497 
498  if(ntruthconstituents < 3)
499  continue;
501  if (m_truthjetpt < m_minjetpt)
502  continue;
503 
504  m_truthjeteta = truthJet->get_eta();
505  m_truthjetpx = truthJet->get_px();
506  m_truthjetpy = truthJet->get_py();
507  m_truthjetpz = truthJet->get_pz();
508  m_truthjetphi = truthJet->get_phi();
509  m_truthjetp = truthJet->get_p();
510  m_truthjetenergy = truthJet->get_e();
511 
512  m_recojetpt = 0;
513  m_recojetid = 0;
514  m_recojetpx = 0;
515  m_recojetpy = 0;
516  m_recojetpz = 0;
517  m_recojetphi = 0;
518  m_recojetp = 0;
519  m_recojetenergy = 0;
520  m_dR = -99;
521  float closestjet = 9999;
523  for (JetMap::Iter recoIter = reco_jets->begin();
524  recoIter != reco_jets->end();
525  ++recoIter)
526  {
527  const Jet *recoJet = recoIter->second;
528  m_recojetpt = recoJet->get_pt();
530  continue;
531 
532  m_recojeteta = recoJet->get_eta();
533  m_recojetphi = recoJet->get_phi();
534 
535  if (Verbosity() > 1)
536  {
537  cout << "matching by distance jet" << endl;
538  }
539 
540  float dphi = m_recojetphi - m_truthjetphi;
541  if (dphi > TMath::Pi())
542  dphi -= TMath::Pi() * 2.;
543  if (dphi < -1 * TMath::Pi())
544  dphi += TMath::Pi() * 2.;
545 
546  float deta = m_recojeteta - m_truthjeteta;
549  m_dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
550 
553  if (m_dR < truth_jets->get_par() && m_dR < closestjet)
554  {
555  // Get reco jet characteristics
556  m_recojetid = recoJet->get_id();
557  m_recojetpx = recoJet->get_px();
558  m_recojetpy = recoJet->get_py();
559  m_recojetpz = recoJet->get_pz();
560  m_recojetphi = recoJet->get_phi();
561  m_recojetp = recoJet->get_p();
562  m_recojetenergy = recoJet->get_e();
563 
564  }
565  }
566 
567 
569  m_truthjettree->Fill();
570  }
571 }
572 
577 {
579  JetMap *reco_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_r04");
581  JetMap *truth_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Truth_r04");
582 
583  if(!m_jetEvalStack)
584  {
585  m_jetEvalStack = new JetEvalStack(topNode, "AntiKt_Tower_r04",
586  "AntiKt_Truth_r04");
587  }
588  m_jetEvalStack->next_event(topNode);
590  if (!reco_jets)
591  {
592  cout << PHWHERE
593  << "Reconstructed jet node is missing, can't collect reconstructed jets"
594  << endl;
595  return;
596  }
597 
598  if (Verbosity() > 1)
599  {
600  cout << "Get all Reco Jets" << endl;
601  }
602 
604  for (JetMap::Iter recoIter = reco_jets->begin();
605  recoIter != reco_jets->end();
606  ++recoIter)
607  {
608  Jet *recoJet = recoIter->second;
609  m_recojetpt = recoJet->get_pt();
610  if (m_recojetpt < m_minjetpt)
611  continue;
612 
613  m_recojeteta = recoJet->get_eta();
614 
615  // Get reco jet characteristics
616  m_recojetid = recoJet->get_id();
617  m_recojetpx = recoJet->get_px();
618  m_recojetpy = recoJet->get_py();
619  m_recojetpz = recoJet->get_pz();
620  m_recojetphi = recoJet->get_phi();
621  m_recojetp = recoJet->get_p();
622  m_recojetenergy = recoJet->get_e();
623 
624  if (Verbosity() > 1)
625  {
626  cout << "matching by distance jet" << endl;
627  }
628 
630  m_truthjetid = 0;
631  m_truthjetp = 0;
632  m_truthjetphi = 0;
633  m_truthjeteta = 0;
634  m_truthjetpt = 0;
635  m_truthjetenergy = 0;
636  m_truthjetpx = 0;
637  m_truthjetpy = 0;
638  m_truthjetpz = 0;
639 
640  Jet *truthjet = recoeval->max_truth_jet_by_energy(recoJet);
641  if(truthjet){
642  m_truthjetid = truthjet->get_id();
643  m_truthjetp = truthjet->get_p();
644  m_truthjetpx = truthjet->get_px();
645  m_truthjetpy = truthjet->get_py();
646  m_truthjetpz = truthjet->get_pz();
647  m_truthjeteta = truthjet->get_eta();
648  m_truthjetphi = truthjet->get_phi();
649  m_truthjetenergy = truthjet->get_e();
652  }
653 
655  else if(truth_jets)
656  {
659  float closestjet = 9999;
660  for (JetMap::Iter truthIter = truth_jets->begin();
661  truthIter != truth_jets->end();
662  ++truthIter)
663  {
664  const Jet *truthJet = truthIter->second;
665 
666  float thisjetpt = truthJet->get_pt();
667  if (thisjetpt < m_minjetpt)
668  continue;
669 
670  float thisjeteta = truthJet->get_eta();
671  float thisjetphi = truthJet->get_phi();
672 
673  float dphi = m_recojetphi - thisjetphi;
674  if (dphi >TMath::Pi())
675  dphi -= TMath::Pi() * 2.;
676  if (dphi < -1. * TMath::Pi())
677  dphi += TMath::Pi() * 2.;
678 
679  float deta = m_recojeteta - thisjeteta;
682  m_dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
683 
686  if (m_dR < reco_jets->get_par() && m_dR < closestjet)
687  {
688  m_truthjetid = -9999;
689  m_truthjetp = truthJet->get_p();
690  m_truthjetphi = truthJet->get_phi();
691  m_truthjeteta = truthJet->get_eta();
692  m_truthjetpt = truthJet->get_pt();
693  m_truthjetenergy = truthJet->get_e();
694  m_truthjetpx = truthJet->get_px();
695  m_truthjetpy = truthJet->get_py();
696  m_truthjetpz = truthJet->get_pz();
697  closestjet = m_dR;
698  }
699  }
700  }
701  m_recojettree->Fill();
702  }
703 }
704 
712 {
716  RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
717 
718  if (!clusters)
719  {
720  cout << PHWHERE
721  << "EMCal cluster node is missing, can't collect EMCal clusters"
722  << endl;
723  return;
724  }
725 
727  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
728  if (!vertexmap)
729  {
730  cout << "AnaTutorial::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;
731  assert(vertexmap); // force quit
732 
733  return;
734  }
735 
736  if (vertexmap->empty())
737  {
738  cout << "AnaTutorial::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;
739  return;
740  }
741 
742  GlobalVertex *vtx = vertexmap->begin()->second;
743  if (vtx == nullptr)
744  return;
745 
747  CaloTriggerInfo *trigger = findNode::getClass<CaloTriggerInfo>(topNode, "CaloTriggerInfo");
748 
750  if(trigger)
751  {
752  m_E_4x4 = trigger->get_best_EMCal_4x4_E();
753  }
754  RawClusterContainer::ConstRange begin_end = clusters->getClusters();
756 
758  for (clusIter = begin_end.first;
759  clusIter != begin_end.second;
760  ++clusIter)
761  {
763  const RawCluster *cluster = clusIter->second;
764 
769  CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
770  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster, vertex);
771  m_clusenergy = E_vec_cluster.mag();
772  m_cluseta = E_vec_cluster.pseudoRapidity();
773  m_clustheta = E_vec_cluster.getTheta();
774  m_cluspt = E_vec_cluster.perp();
775  m_clusphi = E_vec_cluster.getPhi();
776 
777  if (m_cluspt < m_mincluspt)
778  continue;
779 
780  m_cluspx = m_cluspt * cos(m_clusphi);
781  m_cluspy = m_cluspt * sin(m_clusphi);
783 
784  //fill the cluster tree with all emcal clusters
785  m_clustertree->Fill();
786  }
787 }
788 
794 {
795  m_recojettree = new TTree("jettree", "A tree with reconstructed jets");
796  m_recojettree->Branch("m_recojetpt", &m_recojetpt, "m_recojetpt/D");
797  m_recojettree->Branch("m_recojetid", &m_recojetid, "m_recojetid/I");
798  m_recojettree->Branch("m_recojetpx", &m_recojetpx, "m_recojetpx/D");
799  m_recojettree->Branch("m_recojetpy", &m_recojetpy, "m_recojetpy/D");
800  m_recojettree->Branch("m_recojetpz", &m_recojetpz, "m_recojetpz/D");
801  m_recojettree->Branch("m_recojetphi", &m_recojetphi, "m_recojetphi/D");
802  m_recojettree->Branch("m_recojeteta", &m_recojeteta, "m_recojeteta/D");
803  m_recojettree->Branch("m_recojetenergy", &m_recojetenergy, "m_recojetenergy/D");
804  m_recojettree->Branch("m_truthjetid", &m_truthjetid, "m_truthjetid/I");
805  m_recojettree->Branch("m_truthjetp", &m_truthjetp, "m_truthjetp/D");
806  m_recojettree->Branch("m_truthjetphi", &m_truthjetphi, "m_truthjetphi/D");
807  m_recojettree->Branch("m_truthjeteta", &m_truthjeteta, "m_truthjeteta/D");
808  m_recojettree->Branch("m_truthjetpt", &m_truthjetpt, "m_truthjetpt/D");
809  m_recojettree->Branch("m_truthjetenergy", &m_truthjetenergy, "m_truthjetenergy/D");
810  m_recojettree->Branch("m_truthjetpx", &m_truthjetpx, "m_truthjetpx/D");
811  m_recojettree->Branch("m_truthjetpy", &m_truthjetpy, "m_truthjyetpy/D");
812  m_recojettree->Branch("m_truthjetpz", &m_truthjetpz, "m_truthjetpz/D");
813  m_recojettree->Branch("m_dR", &m_dR, "m_dR/D");
814 
815  m_truthjettree = new TTree("truthjettree", "A tree with truth jets");
816  m_truthjettree->Branch("m_truthjetid", &m_truthjetid, "m_truthjetid/I");
817  m_truthjettree->Branch("m_truthjetp", &m_truthjetp, "m_truthjetp/D");
818  m_truthjettree->Branch("m_truthjetphi", &m_truthjetphi, "m_truthjetphi/D");
819  m_truthjettree->Branch("m_truthjeteta", &m_truthjeteta, "m_truthjeteta/D");
820  m_truthjettree->Branch("m_truthjetpt", &m_truthjetpt, "m_truthjetpt/D");
821  m_truthjettree->Branch("m_truthjetenergy", &m_truthjetenergy, "m_truthjetenergy/D");
822  m_truthjettree->Branch("m_truthjetpx", &m_truthjetpx, "m_truthjetpx/D");
823  m_truthjettree->Branch("m_truthjetpy", &m_truthjetpy, "m_truthjetpy/D");
824  m_truthjettree->Branch("m_truthjetpz", &m_truthjetpz, "m_truthjetpz/D");
825  m_truthjettree->Branch("m_dR", &m_dR, "m_dR/D");
826  m_truthjettree->Branch("m_recojetpt", &m_recojetpt, "m_recojetpt/D");
827  m_truthjettree->Branch("m_recojetid", &m_recojetid, "m_recojetid/I");
828  m_truthjettree->Branch("m_recojetpx", &m_recojetpx, "m_recojetpx/D");
829  m_truthjettree->Branch("m_recojetpy", &m_recojetpy, "m_recojetpy/D");
830  m_truthjettree->Branch("m_recojetpz", &m_recojetpz, "m_recojetpz/D");
831  m_truthjettree->Branch("m_recojetphi", &m_recojetphi, "m_recojetphi/D");
832  m_truthjettree->Branch("m_recojeteta", &m_recojeteta, "m_recojeteta/D");
833  m_truthjettree->Branch("m_recojetenergy", &m_recojetenergy, "m_recojetenergy/D");
834  m_tracktree = new TTree("tracktree", "A tree with svtx tracks");
835  m_tracktree->Branch("m_tr_px", &m_tr_px, "m_tr_px/D");
836  m_tracktree->Branch("m_tr_py", &m_tr_py, "m_tr_py/D");
837  m_tracktree->Branch("m_tr_pz", &m_tr_pz, "m_tr_pz/D");
838  m_tracktree->Branch("m_tr_p", &m_tr_p, "m_tr_p/D");
839  m_tracktree->Branch("m_tr_pt", &m_tr_pt, "m_tr_pt/D");
840  m_tracktree->Branch("m_tr_phi", &m_tr_phi, "m_tr_phi/D");
841  m_tracktree->Branch("m_tr_eta", &m_tr_eta, "m_tr_eta/D");
842  m_tracktree->Branch("m_charge", &m_charge, "m_charge/I");
843  m_tracktree->Branch("m_chisq", &m_chisq, "m_chisq/D");
844  m_tracktree->Branch("m_ndf", &m_ndf, "m_ndf/I");
845  m_tracktree->Branch("m_dca", &m_dca, "m_dca/D");
846  m_tracktree->Branch("m_tr_x", &m_tr_x, "m_tr_x/D");
847  m_tracktree->Branch("m_tr_y", &m_tr_y, "m_tr_y/D");
848  m_tracktree->Branch("m_tr_z", &m_tr_z, "m_tr_z/D");
849  m_tracktree->Branch("m_truth_is_primary", &m_truth_is_primary, "m_truth_is_primary/I");
850  m_tracktree->Branch("m_truthtrackpx", &m_truthtrackpx, "m_truthtrackpx/D");
851  m_tracktree->Branch("m_truthtrackpy", &m_truthtrackpy, "m_truthtrackpy/D");
852  m_tracktree->Branch("m_truthtrackpz", &m_truthtrackpz, "m_truthtrackpz/D");
853  m_tracktree->Branch("m_truthtrackp", &m_truthtrackp, "m_truthtrackp/D");
854  m_tracktree->Branch("m_truthtracke", &m_truthtracke, "m_truthtracke/D");
855  m_tracktree->Branch("m_truthtrackpt", &m_truthtrackpt, "m_truthtrackpt/D");
856  m_tracktree->Branch("m_truthtrackphi", &m_truthtrackphi, "m_truthtrackphi/D");
857  m_tracktree->Branch("m_truthtracketa", &m_truthtracketa, "m_truthtracketa/D");
858  m_tracktree->Branch("m_truthtrackpid", &m_truthtrackpid, "m_truthtrackpid/I");
859 
860  m_hepmctree = new TTree("hepmctree", "A tree with hepmc truth particles");
861  m_hepmctree->Branch("m_partid1", &m_partid1, "m_partid1/I");
862  m_hepmctree->Branch("m_partid2", &m_partid2, "m_partid2/I");
863  m_hepmctree->Branch("m_x1", &m_x1, "m_x1/D");
864  m_hepmctree->Branch("m_x2", &m_x2, "m_x2/D");
865  m_hepmctree->Branch("m_mpi", &m_mpi, "m_mpi/I");
866  m_hepmctree->Branch("m_process_id", &m_process_id, "m_process_id/I");
867  m_hepmctree->Branch("m_truthenergy", &m_truthenergy, "m_truthenergy/D");
868  m_hepmctree->Branch("m_trutheta", &m_trutheta, "m_trutheta/D");
869  m_hepmctree->Branch("m_truthphi", &m_truthphi, "m_truthphi/D");
870  m_hepmctree->Branch("m_truthpx", &m_truthpx, "m_truthpx/D");
871  m_hepmctree->Branch("m_truthpy", &m_truthpy, "m_truthpy/D");
872  m_hepmctree->Branch("m_truthpz", &m_truthpz, "m_truthpz/D");
873  m_hepmctree->Branch("m_truthpt", &m_truthpt, "m_truthpt/D");
874  m_hepmctree->Branch("m_numparticlesinevent", &m_numparticlesinevent, "m_numparticlesinevent/I");
875  m_hepmctree->Branch("m_truthpid", &m_truthpid, "m_truthpid/I");
876 
877  m_truthtree = new TTree("truthg4tree", "A tree with truth g4 particles");
878  m_truthtree->Branch("m_truthenergy", &m_truthenergy, "m_truthenergy/D");
879  m_truthtree->Branch("m_truthp", &m_truthp, "m_truthp/D");
880  m_truthtree->Branch("m_truthpx", &m_truthpx, "m_truthpx/D");
881  m_truthtree->Branch("m_truthpy", &m_truthpy, "m_truthpy/D");
882  m_truthtree->Branch("m_truthpz", &m_truthpz, "m_truthpz/D");
883  m_truthtree->Branch("m_truthpt", &m_truthpt, "m_truthpt/D");
884  m_truthtree->Branch("m_truthphi", &m_truthphi, "m_truthphi/D");
885  m_truthtree->Branch("m_trutheta", &m_trutheta, "m_trutheta/D");
886  m_truthtree->Branch("m_truthpid", &m_truthpid, "m_truthpid/I");
887 
888  m_clustertree = new TTree("clustertree", "A tree with emcal clusters");
889  m_clustertree->Branch("m_clusenergy", &m_clusenergy, "m_clusenergy/D");
890  m_clustertree->Branch("m_cluseta", &m_cluseta, "m_cluseta/D");
891  m_clustertree->Branch("m_clustheta", &m_clustheta, "m_clustheta/D");
892  m_clustertree->Branch("m_cluspt", &m_cluspt, "m_cluspt/D");
893  m_clustertree->Branch("m_clusphi", &m_clusphi, "m_clusphi/D");
894  m_clustertree->Branch("m_cluspx", &m_cluspx, "m_cluspx/D");
895  m_clustertree->Branch("m_cluspy", &m_cluspy, "m_cluspy/D");
896  m_clustertree->Branch("m_cluspz", &m_cluspz, "m_cluspz/D");
897  m_clustertree->Branch("m_E_4x4", &m_E_4x4, "m_E_4x4/D");
898 }
899 
906 {
907  m_outfile = new TFile();
908  m_phi_h = new TH1F();
909  m_eta_phi_h = new TH2F();
910 
911  m_partid1 = -99;
912  m_partid2 = -99;
913  m_x1 = -99;
914  m_x2 = -99;
915  m_mpi = -99;
916  m_process_id = -99;
917  m_truthenergy = -99;
918  m_trutheta = -99;
919  m_truthphi = -99;
920  m_truthp = -99;
921  m_truthpx = -99;
922  m_truthpy = -99;
923  m_truthpz = -99;
924  m_truthpt = -99;
925  m_numparticlesinevent = -99;
926  m_truthpid = -99;
927 
928  m_tr_px = -99;
929  m_tr_py = -99;
930  m_tr_pz = -99;
931  m_tr_p = -99;
932  m_tr_pt = -99;
933  m_tr_phi = -99;
934  m_tr_eta = -99;
935  m_charge = -99;
936  m_chisq = -99;
937  m_ndf = -99;
938  m_dca = -99;
939  m_tr_x = -99;
940  m_tr_y = -99;
941  m_tr_z = -99;
942  m_truth_is_primary = -99;
943  m_truthtrackpx = -99;
944  m_truthtrackpy = -99;
945  m_truthtrackpz = -99;
946  m_truthtrackp = -99;
947  m_truthtracke = -99;
948  m_truthtrackpt = -99;
949  m_truthtrackphi = -99;
950  m_truthtracketa = -99;
951  m_truthtrackpid = -99;
952 
953  m_recojetpt = -99;
954  m_recojetid = -99;
955  m_recojetpx = -99;
956  m_recojetpy = -99;
957  m_recojetpz = -99;
958  m_recojetphi = -99;
959  m_recojetp = -99;
960  m_recojetenergy = -99;
961  m_recojeteta = -99;
962  m_truthjetid = -99;
963  m_truthjetp = -99;
964  m_truthjetphi = -99;
965  m_truthjeteta = -99;
966  m_truthjetpt = -99;
967  m_truthjetenergy = -99;
968  m_truthjetpx = -99;
969  m_truthjetpy = -99;
970  m_truthjetpz = -99;
971  m_dR = -99;
972 }
973 
974