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>
33 #include <HepMC/GenEvent.h>
34 #include <HepMC/GenVertex.h>
74 , m_outfilename(filename)
78 , m_analyzeTracks(
true)
79 , m_analyzeClusters(
true)
81 , m_analyzeTruth(
false)
108 cout <<
"Beginning Init in AnaTutorialECCE" << endl;
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);
127 cout <<
"Beginning process_event in AnaTutorialECCE" << endl;
165 cout <<
"Ending AnaTutorialECCE analysis package" << endl;
208 cout <<
"Finished AnaTutorialECCE analysis package" << endl;
223 PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode,
"PHHepMCGenEventMap");
229 <<
"HEPMC event map node is missing, can't collected HEPMC truth particles"
237 cout <<
"Getting HEPMC truth particles " << endl;
243 eventIter != hepmceventmap->
end();
252 HepMC::GenEvent *truthevent = hepmcevent->
getEvent();
256 <<
"no evt pointer under phhepmvgeneventmap found "
262 HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
267 m_x1 = pdfinfo->x1();
268 m_x2 = pdfinfo->x2();
271 m_mpi = truthevent->mpi();
278 cout <<
" Iterating over an event" << endl;
281 for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
282 iter != truthevent->particles_end();
289 m_trutheta = (*iter)->momentum().pseudoRapidity();
318 <<
"PHG4TruthInfoContainer node is missing, can't collect G4 truth particles"
328 iter != range.second;
363 SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode,
"TrackMap");
364 EICPIDParticleContainer *pidcontainer = findNode::getClass<EICPIDParticleContainer>(topNode,
"EICPIDParticleMap");
366 if (
Verbosity() > 1 and pidcontainer ==
nullptr)
368 cout <<
"EICPIDParticleContainer named EICPIDParticleMap does not exist. Skip saving PID info" << endl;
374 <<
"TrackMap node is missing, can't collect tracks"
384 cout <<
"Get the tracks" << endl;
387 iter != trackmap->
end();
420 std::cout <<
"Skipping non fast track sim object..." << std::endl;
448 pidcontainer->findEICPIDParticle(track->
get_id());
471 cout <<
"get the truth jets" << endl;
475 JetMap *truth_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Truth_r05");
478 JetMap *reco_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Tower_r05");
490 <<
"Truth jet node is missing, can't collect truth jets"
497 iter != truth_jets->
end();
500 Jet *truthJet = iter->second;
504 std::set<PHG4Particle *> truthjetcomp =
506 int ntruthconstituents = 0;
508 for (std::set<PHG4Particle *>::iterator iter2 = truthjetcomp.begin();
509 iter2 != truthjetcomp.end();
516 cout <<
"no truth particles in the jet??" << endl;
520 ntruthconstituents++;
523 if (ntruthconstituents < 3)
546 float closestjet = 9999;
549 recoIter != reco_jets->end();
552 const Jet *recoJet = recoIter->second;
562 cout <<
"matching by distance jet" << endl;
566 if (dphi > TMath::Pi())
567 dphi -= TMath::Pi() * 2.;
568 if (dphi < -1 * TMath::Pi())
569 dphi += TMath::Pi() * 2.;
574 m_dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
578 if (m_dR < truth_jets->get_par() &&
m_dR < closestjet)
602 JetMap *reco_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Tower_r05");
604 JetMap *truth_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Truth_r05");
616 <<
"Reconstructed jet node is missing, can't collect reconstructed jets"
623 cout <<
"Get all Reco Jets" << endl;
628 recoIter != reco_jets->
end();
631 Jet *recoJet = recoIter->second;
649 cout <<
"matching by distance jet" << endl;
682 float closestjet = 9999;
684 truthIter != truth_jets->end();
687 const Jet *truthJet = truthIter->second;
689 float thisjetpt = truthJet->
get_pt();
693 float thisjeteta = truthJet->
get_eta();
694 float thisjetphi = truthJet->
get_phi();
697 if (dphi > TMath::Pi())
698 dphi -= TMath::Pi() * 2.;
699 if (dphi < -1. * TMath::Pi())
700 dphi += TMath::Pi() * 2.;
705 m_dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
709 if (m_dR < reco_jets->get_par() &&
m_dR < closestjet)
739 RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode,
"CLUSTER_BECAL");
744 <<
"EMCal cluster node is missing, can't collect EMCal clusters"
750 GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode,
"GlobalVertexMap");
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;
759 if (vertexmap->
empty())
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;
770 CaloTriggerInfo *trigger = findNode::getClass<CaloTriggerInfo>(topNode,
"CaloTriggerInfo");
781 for (clusIter = begin_end.first;
782 clusIter != begin_end.second;
795 m_cluseta = E_vec_cluster.pseudoRapidity();
821 m_recojettree =
new TTree(
"jettree",
"A tree with reconstructed jets");
841 m_truthjettree =
new TTree(
"truthjettree",
"A tree with truth jets");
860 m_tracktree =
new TTree(
"tracktree",
"A tree with svtx tracks");
889 m_hepmctree =
new TTree(
"hepmctree",
"A tree with hepmc truth particles");
906 m_truthtree =
new TTree(
"truthg4tree",
"A tree with truth g4 particles");
917 m_clustertree =
new TTree(
"clustertree",
"A tree with emcal clusters");