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>
30 #include <HepMC/GenEvent.h>
31 #include <HepMC/GenVertex.h>
71 , m_outfilename(filename)
75 , m_analyzeTracks(
true)
76 , m_analyzeClusters(
true)
78 , m_analyzeTruth(
false)
105 cout <<
"Beginning Init in AnaTutorial" << endl;
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);
124 cout <<
"Beginning process_event in AnaTutorial" << endl;
162 cout <<
"Ending AnaTutorial analysis package" << endl;
205 cout <<
"Finished AnaTutorial analysis package" << endl;
220 PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode,
"PHHepMCGenEventMap");
226 <<
"HEPMC event map node is missing, can't collected HEPMC truth particles"
234 cout <<
"Getting HEPMC truth particles " << endl;
240 eventIter != hepmceventmap->
end();
249 HepMC::GenEvent *truthevent = hepmcevent->
getEvent();
253 <<
"no evt pointer under phhepmvgeneventmap found "
259 HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
264 m_x1 = pdfinfo->x1();
265 m_x2 = pdfinfo->x2();
268 m_mpi = truthevent->mpi();
275 cout <<
" Iterating over an event" << endl;
278 for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
279 iter != truthevent->particles_end();
286 m_trutheta = (*iter)->momentum().pseudoRapidity();
315 <<
"PHG4TruthInfoContainer node is missing, can't collect G4 truth particles"
325 iter != range.second;
361 SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxTrackMap");
366 <<
"SvtxTrackMap node is missing, can't collect tracks"
388 cout <<
"Get the SVTX tracks" << endl;
391 iter != trackmap->
end();
446 cout <<
"get the truth jets" << endl;
450 JetMap *truth_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Truth_r04");
453 JetMap *reco_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Tower_r04");
465 <<
"Truth jet node is missing, can't collect truth jets"
472 iter != truth_jets->
end();
475 Jet *truthJet = iter->second;
479 std::set<PHG4Particle *> truthjetcomp =
481 int ntruthconstituents = 0;
483 for (std::set<PHG4Particle *>::iterator iter2 = truthjetcomp.begin();
484 iter2 != truthjetcomp.end();
491 cout <<
"no truth particles in the jet??" << endl;
495 ntruthconstituents++;
498 if(ntruthconstituents < 3)
521 float closestjet = 9999;
524 recoIter != reco_jets->end();
527 const Jet *recoJet = recoIter->second;
537 cout <<
"matching by distance jet" << endl;
541 if (dphi > TMath::Pi())
542 dphi -= TMath::Pi() * 2.;
543 if (dphi < -1 * TMath::Pi())
544 dphi += TMath::Pi() * 2.;
549 m_dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
553 if (m_dR < truth_jets->get_par() &&
m_dR < closestjet)
579 JetMap *reco_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Tower_r04");
581 JetMap *truth_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Truth_r04");
593 <<
"Reconstructed jet node is missing, can't collect reconstructed jets"
600 cout <<
"Get all Reco Jets" << endl;
605 recoIter != reco_jets->
end();
608 Jet *recoJet = recoIter->second;
626 cout <<
"matching by distance jet" << endl;
659 float closestjet = 9999;
661 truthIter != truth_jets->end();
664 const Jet *truthJet = truthIter->second;
666 float thisjetpt = truthJet->
get_pt();
670 float thisjeteta = truthJet->
get_eta();
671 float thisjetphi = truthJet->
get_phi();
674 if (dphi >TMath::Pi())
675 dphi -= TMath::Pi() * 2.;
676 if (dphi < -1. * TMath::Pi())
677 dphi += TMath::Pi() * 2.;
682 m_dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
686 if (m_dR < reco_jets->get_par() &&
m_dR < closestjet)
716 RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode,
"CLUSTER_CEMC");
721 <<
"EMCal cluster node is missing, can't collect EMCal clusters"
727 GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode,
"GlobalVertexMap");
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;
736 if (vertexmap->
empty())
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;
747 CaloTriggerInfo *trigger = findNode::getClass<CaloTriggerInfo>(topNode,
"CaloTriggerInfo");
758 for (clusIter = begin_end.first;
759 clusIter != begin_end.second;
772 m_cluseta = E_vec_cluster.pseudoRapidity();
795 m_recojettree =
new TTree(
"jettree",
"A tree with reconstructed jets");
815 m_truthjettree =
new TTree(
"truthjettree",
"A tree with truth jets");
834 m_tracktree =
new TTree(
"tracktree",
"A tree with svtx tracks");
860 m_hepmctree =
new TTree(
"hepmctree",
"A tree with hepmc truth particles");
877 m_truthtree =
new TTree(
"truthg4tree",
"A tree with truth g4 particles");
888 m_clustertree =
new TTree(
"clustertree",
"A tree with emcal clusters");