analysis/ResonanceJetTaggingOutputs/BuildResonanceJetTaggingTree.cc

0001 #include "BuildResonanceJetTaggingTree.h"
0002
0003 #include <resonancejettagging/ResonanceJetTagging.h>
0004
0005 #include <phool/phool.h>
0006
0007 /// Jet includes
0008 #include <jetbase/Jet.h>
0009 #include <jetbase/JetContainerv1.h>
0010 #include <jetbase/Jetv2.h>
0011
0012 /// Tracking includes
0013 #include <trackbase_historic/SvtxPHG4ParticleMap_v1.h>
0014 #include <trackbase_historic/SvtxTrack.h>
0015 #include <trackbase_historic/SvtxTrackMap.h>
0016
0017 #include <g4eval/SvtxEvalStack.h>   // for SvtxEvalStack
0018 #include <g4eval/SvtxTrackEval.h>   // for SvtxTrackEval
0019
0020 /// HEPMC truth includes
0021 #pragma GCC diagnostic push
0022 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0023 #include <HepMC/GenEvent.h>
0024 #include <HepMC/GenVertex.h>
0025 #pragma GCC diagnostic pop
0026
0027 #include <phhepmc/PHHepMCGenEvent.h>
0028 #include <phhepmc/PHHepMCGenEventMap.h>
0029 #include <phhepmc/PHGenIntegral.h>
0030
0031 /// Fun4All includes
0032 #include <fun4all/Fun4AllReturnCodes.h>
0033
0034 #include <phool/PHCompositeNode.h>
0035 #include <phool/getClass.h>
0036
0037 #include <g4main/PHG4Particle.h>            // for PHG4Particle
0038 #include <g4main/PHG4TruthInfoContainer.h>  // for PHG4TruthInfoContainer
0039 #include <g4main/PHG4VtxPoint.h>            // for PHG4VtxPoint
0040
0041 #include <kfparticle_sphenix/KFParticle_truthAndDetTools.h>
0042
0043 #include <KFParticle.h>
0044 #include <kfparticle_sphenix/KFParticle_Container.h>
0045
0046 /// ROOT includes
0047 #include <TFile.h>
0048 #include <TTree.h>
0049 #include <TH1I.h>
0050
0051 /// C++ includes
0052 #include <cassert>
0053 #include <sstream>
0054 #include <string>
0055
0056 /**
0057  * BuildResonanceJetTaggingTree is a class developed to reconstruct jets containing a D-meson
0058  * The class can be adapted to tag jets using any kind of particle
0059  * Author: Antonio Silva (antonio.sphenix@gmail.com)
0060  * Contributor: Jakub Kvapil (jakub.kvapil@cern.ch)
0061  */
0062
0063 /**
0064  * Constructor of module
0065  */
0066 BuildResonanceJetTaggingTree::BuildResonanceJetTaggingTree(const std::string &name, const std::string &filename, const ResonanceJetTagging::TAG tag)
0067   : SubsysReco(name)
0068   , m_outfilename(filename)
0069   , m_tagcontainer_name("")
0070   , m_jetcontainer_name("")
0071   , m_truth_jetcontainer_name("")
0072   , m_dorec(true)
0073   , m_dotruth(false)
0074   , m_stable_mother(true)
0075   , m_nDaughters(0)
0076   , m_svtx_evalstack(nullptr)
0077   , m_trackeval(nullptr)
0078   , m_tag_particle(tag)
0079   , m_tag_pdg(0)
0080   , m_outfile(nullptr)
0081   , m_eventcount_h(nullptr)
0082   , m_taggedjettree(nullptr)
0083 {
0084   /// Initialize variables and trees so we don't accidentally access
0085   /// memory that was never allocated
0086   initializeVariables();
0087   initializeTrees();
0088
0089   switch (m_tag_particle) {
0090     case ResonanceJetTagging::TAG::D0:
0091       m_tag_pdg = 421;
0092       m_nDaughters = 2;
0093       break;
0094     case ResonanceJetTagging::TAG::D0TOK3PI:
0095       m_tag_pdg = 421;
0096       m_nDaughters = 4;
0097       break;
0098     case ResonanceJetTagging::TAG::DPLUS:
0099       m_tag_pdg = 411;
0100       m_nDaughters = 3;
0101       break;
0102     case ResonanceJetTagging::TAG::DSTAR:
0103       m_tag_pdg = 413;
0104       m_nDaughters = 0;
0105       break;
0106     case ResonanceJetTagging::TAG::JPSI:
0107       m_tag_pdg = 443;
0108       m_nDaughters = 2;
0109       break;
0110     case ResonanceJetTagging::TAG::K0:
0111       m_tag_pdg = 310;
0112       m_nDaughters = 2;
0113       break;
0114     case ResonanceJetTagging::TAG::GAMMA:
0115       m_tag_pdg = 22;
0116       m_nDaughters = 0;
0117       break;
0118     case ResonanceJetTagging::TAG::ELECTRON:
0119       m_tag_pdg = 11;
0120       m_nDaughters = 0;
0121       break;
0122     case ResonanceJetTagging::TAG::LAMBDAC:
0123       m_tag_pdg = 4122;
0124       m_nDaughters = 3;
0125       break;
0126     case ResonanceJetTagging::TAG::LAMBDAS:
0127       m_tag_pdg = 3122;
0128       m_nDaughters = 2;
0129       break;
0130   }
0131 }
0132
0133 /**
0134  * Destructor of module
0135  */
0136 BuildResonanceJetTaggingTree::~BuildResonanceJetTaggingTree()
0137 {
0138
0139 }
0140
0141 /**
0142  * Initialize the module and prepare looping over events
0143  */
0144 int BuildResonanceJetTaggingTree::Init(PHCompositeNode */*topNode*/)
0145 {
0146   if (Verbosity() > 5)
0147   {
0148     std::cout << "Beginning Init in BuildResonanceJetTaggingTree" << std::endl;
0149   }
0150
0151   return 0;
0152 }
0153
0154 /**
0155  * Main workhorse function where each event is looped over and
0156  * data from each event is collected from the node tree for analysis
0157  */
0158 int BuildResonanceJetTaggingTree::process_event(PHCompositeNode *topNode)
0159 {
0160   if (Verbosity() > 5)
0161   {
0162     std::cout << "Beginning process_event in BuildResonanceJetTaggingTree" << std::endl;
0163   }
0164
0165   if(m_nDaughters == 0)
0166   {
0167     std::cout<<"ERROR: Number of Decay Daughters Not Set, ABORTING!";
0168     return Fun4AllReturnCodes::ABORTRUN;
0169   }
0170
0171   m_eventcount_h->Fill(0);
0172
0173   switch (m_tag_particle) {
0174     case ResonanceJetTagging::TAG::D0:
0175       [[fallthrough]];
0176     case ResonanceJetTagging::TAG::D0TOK3PI:
0177       [[fallthrough]];
0178     case ResonanceJetTagging::TAG::DPLUS:
0179       [[fallthrough]];
0180     case ResonanceJetTagging::TAG::LAMBDAC:
0181       [[fallthrough]];
0182     case ResonanceJetTagging::TAG::LAMBDAS:
0183       [[fallthrough]];
0184     case ResonanceJetTagging::TAG::K0:
0185       return loopHFHadronic(topNode);
0186       break;
0187     default:
0188       std::cout<<"ERROR: Fill Tree Function Not Set, ABORTING!";
0189       return Fun4AllReturnCodes::ABORTRUN;
0190       break;
0191   }
0192
0193   return Fun4AllReturnCodes::EVENT_OK;
0194 }
0195
0196 /**
0197  * End the module and finish any data collection. Clean up any remaining
0198  * loose ends
0199  */
0200 int BuildResonanceJetTaggingTree::End(PHCompositeNode *topNode)
0201 {
0202   if (Verbosity() > 1)
0203   {
0204     std::cout << "Ending BuildResonanceJetTaggingTree analysis package" << std::endl;
0205   }
0206
0207   /// Change to the outfile
0208   m_outfile->cd();
0209
0210   if(m_dotruth) {
0211     PHGenIntegral *phgen = findNode::getClass<PHGenIntegral>(topNode, "PHGenIntegral");
0212     if(phgen)
0213       {
0214     m_intlumi = phgen->get_Integrated_Lumi();
0215     m_nprocessedevents = phgen->get_N_Processed_Event();
0216     m_nacceptedevents = phgen->get_N_Generator_Accepted_Event();
0217     m_xsecprocessedevents = phgen->get_CrossSection_Processed_Event();
0218     m_xsecacceptedevents = phgen->get_CrossSection_Generator_Accepted_Event();
0219     m_runinfo->Fill();
0220       }
0221     m_runinfo->Write();
0222   }
0223
0224   m_taggedjettree->Write();
0225   m_eventcount_h->Write();
0226
0227   m_outfile->Close();
0228
0229   delete m_outfile;
0230
0231   if (Verbosity() > 1)
0232   {
0233     std::cout << "Finished BuildResonanceJetTaggingTree analysis package" << std::endl;
0234   }
0235
0236   return 0;
0237 }
0238
0239 int BuildResonanceJetTaggingTree::loopHFHadronic(PHCompositeNode *topNode)
0240 {
0241   KFParticle_Container* kfContainer = nullptr;
0242
0243   if(m_dorec)
0244   {
0245     m_taggedJetContainer = getJetContainerFromNode(topNode, m_jetcontainer_name);
0246     if(!m_taggedJetContainer) return Fun4AllReturnCodes::ABORTEVENT;
0247
0248     kfContainer = getKFParticleContainerFromNode(topNode, m_tagcontainer_name);
0249     if(!kfContainer) return Fun4AllReturnCodes::ABORTEVENT;
0250   }
0251
0252   HepMC::GenEvent *hepMCGenEvent = nullptr;
0253
0254   if(m_dotruth)
0255   {
0256     m_truth_taggedJetContainer = getJetContainerFromNode(topNode, m_truth_jetcontainer_name);
0257     if(!m_truth_taggedJetContainer) return Fun4AllReturnCodes::ABORTEVENT;
0258
0259     hepMCGenEvent = getGenEventFromNode(topNode, "PHHepMCGenEventMap");
0260     if(!hepMCGenEvent) return Fun4AllReturnCodes::ABORTEVENT;
0261
0262     if (!m_svtx_evalstack)
0263     {
0264       m_svtx_evalstack = new SvtxEvalStack(topNode);
0265
0266       m_trackeval = m_svtx_evalstack->get_track_eval();
0267     }
0268     else
0269     {
0270       m_svtx_evalstack->next_event(topNode);
0271     }
0272
0273   }
0274
0275   m_eventcount_h->Fill(1);
0276
0277   Jet *recTagJet = nullptr;
0278   Jet *genTagJet = nullptr;
0279
0280   KFParticle *recTag = nullptr;
0281   HepMC::GenParticle *genTag = nullptr;
0282
0283   std::vector<int> recDaughtersID;
0284   std::vector<int> recJetIndex;
0285
0286   if(m_dorec)
0287   {
0288     for(auto recJet : *m_taggedJetContainer)
0289     {
0290       recTagJet = recJet;
0291
0292       if(!recTagJet) continue;
0293
0294       for(auto comp : recTagJet->get_comp_vec())
0295       {
0296         if(comp.first == Jet::SRC::VOID)
0297         {
0298           recTag = kfContainer->get(comp.second);
0299         }
0300       }
0301
0302       if(!recTag)
0303       {
0304         continue;
0305       }
0306
0307       recDaughtersID = recTag->DaughterIds();
0308
0309       resetTreeVariables();
0310
0311       m_reco_tag_px = recTag->Px();
0312       m_reco_tag_py = recTag->Py();
0313       m_reco_tag_pz = recTag->Pz();
0314       m_reco_tag_pt = recTag->GetPt();
0315       m_reco_tag_eta = recTag->GetEta();
0316       m_reco_tag_phi = recTag->GetPhi();
0317       m_reco_tag_m = recTag->GetMass();
0318       m_reco_tag_e = recTag->E();
0319
0320       m_reco_jet_px = recTagJet->get_px();
0321       m_reco_jet_py = recTagJet->get_py();
0322       m_reco_jet_pz = recTagJet->get_pz();
0323       m_reco_jet_pt = recTagJet->get_pt();
0324       m_reco_jet_eta = recTagJet->get_eta();
0325       m_reco_jet_phi = recTagJet->get_phi();
0326       m_reco_jet_m = recTagJet->get_mass();
0327       m_reco_jet_e = recTagJet->get_e();
0328
0329       genTagJet = nullptr;
0330       genTag = nullptr;
0331
0332       if(m_dotruth)
0333       {
0334         findMatchedTruthD0(topNode, genTagJet, genTag, recDaughtersID);
0335
0336         if((genTagJet) && (genTag))
0337         {
0338           recJetIndex.push_back(genTagJet->get_id());
0339
0340           m_truth_tag_px = genTag->momentum().px();
0341           m_truth_tag_py = genTag->momentum().py();
0342           m_truth_tag_pz = genTag->momentum().pz();
0343           m_truth_tag_pt = std::sqrt(m_truth_tag_px * m_truth_tag_px + m_truth_tag_py * m_truth_tag_py);
0344           m_truth_tag_eta = genTag->momentum().pseudoRapidity();
0345           m_truth_tag_phi = atan2(m_truth_tag_py, m_truth_tag_px);
0346           m_truth_tag_m = genTag->momentum().m();
0347           m_truth_tag_e = genTag->momentum().e();
0348
0349           m_truth_jet_px = genTagJet->get_px();
0350           m_truth_jet_py = genTagJet->get_py();
0351           m_truth_jet_pz = genTagJet->get_pz();
0352           m_truth_jet_pt = genTagJet->get_pt();
0353           m_truth_jet_eta = genTagJet->get_eta();
0354           m_truth_jet_phi = genTagJet->get_phi();
0355           m_truth_jet_m = genTagJet->get_mass();
0356           m_truth_jet_e = genTagJet->get_e();
0357
0358         }
0359       }
0360
0361       m_taggedjettree->Fill();
0362     }
0363   }
0364
0365   if(m_dotruth)
0366   {
0367     resetTreeVariables();
0368
0369     for(auto mcJet : *m_truth_taggedJetContainer)
0370     {
0371       genTagJet = mcJet;
0372
0373       if(!genTagJet) continue;
0374
0375       //Check if truth was matched to reconstructed
0376       if(m_dorec)
0377       {
0378         if(isReconstructed(genTagJet->get_id(), recJetIndex))
0379         {
0380           continue;
0381         }
0382       }
0383
0384       for(auto comp : genTagJet->get_comp_vec())
0385       {
0386         if(comp.first == Jet::SRC::VOID)
0387         {
0388           genTag = hepMCGenEvent->barcode_to_particle(comp.second);
0389         }
0390       }
0391
0392       if(!genTag)
0393       {
0394         continue;
0395       }
0396
0397       m_truth_tag_px = genTag->momentum().px();
0398       m_truth_tag_py = genTag->momentum().py();
0399       m_truth_tag_pz = genTag->momentum().pz();
0400       m_truth_tag_pt = std::sqrt(m_truth_tag_px * m_truth_tag_px + m_truth_tag_py * m_truth_tag_py);
0401       m_truth_tag_eta = genTag->momentum().pseudoRapidity();
0402       m_truth_tag_phi = atan2(m_truth_tag_py, m_truth_tag_px);
0403       m_truth_tag_m = genTag->momentum().m();
0404       m_truth_tag_e = genTag->momentum().e();
0405
0406       m_truth_jet_px = genTagJet->get_px();
0407       m_truth_jet_py = genTagJet->get_py();
0408       m_truth_jet_pz = genTagJet->get_pz();
0409       m_truth_jet_pt = genTagJet->get_pt();
0410       m_truth_jet_eta = genTagJet->get_eta();
0411       m_truth_jet_phi = genTagJet->get_phi();
0412       m_truth_jet_m = genTagJet->get_mass();
0413       m_truth_jet_e = genTagJet->get_e();
0414
0415       /// iterate over all constituents and add them to the tree
0416       for(auto comp : genTagJet->get_comp_vec())
0417       {
0418         HepMC::GenParticle* constituent = hepMCGenEvent->barcode_to_particle(comp.second);
0419         /// Don't include the tagged particle
0420         if(constituent == genTag)
0421         {
0422           continue;
0423         }
0424
0425         m_truthjet_const_px.push_back(constituent->momentum().px());
0426         m_truthjet_const_py.push_back(constituent->momentum().py());
0427         m_truthjet_const_pz.push_back(constituent->momentum().pz());
0428         m_truthjet_const_e.push_back(constituent->momentum().e());
0429       }
0430
0431       m_taggedjettree->Fill();
0432     }
0433   }
0434
0435   return Fun4AllReturnCodes::EVENT_OK;
0436 }
0437
0438 void BuildResonanceJetTaggingTree::findMatchedTruthD0(PHCompositeNode *topNode, Jet *&mcTagJet, HepMC::GenParticle *&mcTag, std::vector<int> decays)
0439 {
0440   m_truth_taggedJetContainer = getJetContainerFromNode(topNode, m_truth_jetcontainer_name);
0441   if(!m_truth_taggedJetContainer) return;
0442
0443   HepMC::GenEvent *hepMCGenEvent = getGenEventFromNode(topNode, "PHHepMCGenEventMap");
0444   if(!hepMCGenEvent) return;
0445
0446   PHG4Particle *g4particle = nullptr;
0447   PHG4Particle *g4parent = nullptr;
0448   std::vector<HepMC::GenParticle*> mcTags(m_nDaughters);
0449   std::vector<int> daughtersID(m_nDaughters);
0450
0451   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0452
0453   // Truth map
0454   SvtxPHG4ParticleMap_v1 *dst_reco_truth_map = findNode::getClass<SvtxPHG4ParticleMap_v1>(topNode, "SvtxPHG4ParticleMap");
0455
0456   if(m_stable_mother)
0457   {
0458     if (dst_reco_truth_map)
0459     {
0460       for (int idecay = 0; idecay < m_nDaughters; idecay++)
0461       {
0462         std::map<float, std::set<int>> truth_set = dst_reco_truth_map->get(decays[idecay]);
0463         const auto &best_weight = truth_set.rbegin();
0464         int best_truth_id = *best_weight->second.rbegin();
0465         g4particle = truthinfo->GetParticle(best_truth_id);
0466         mcTags[idecay] = getMother(topNode, g4particle);
0467         if (mcTags[idecay] == nullptr)
0468         {
0469           return;
0470         }
0471       }
0472     }
0473     else
0474     {
0475       SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0476       if(!trackmap) return;
0477
0478       for (int idecay = 0; idecay < m_nDaughters; idecay++)
0479       {
0480         SvtxTrack *track = trackmap->get(decays[idecay]);
0481         if(!track) return;
0482         g4particle = m_trackeval->max_truth_particle_by_nclusters(track);
0483
0484         if(!g4particle)
0485         {
0486           return;
0487         }
0488
0489         g4parent = truthinfo->GetParticle(g4particle->get_primary_id());
0490
0491         if(g4parent == nullptr)
0492         {
0493           return;
0494         }
0495
0496         mcTags[idecay] = hepMCGenEvent->barcode_to_particle(g4parent->get_barcode());
0497         if(mcTags[idecay] == nullptr)
0498         {
0499           return;
0500         }
0501       }
0502     }
0503
0504     // check is decays are from the same mother, otherwise it is background
0505     for (int idecay = 1; idecay < m_nDaughters; idecay++)
0506     {
0507       if (mcTags[idecay]->barcode() != mcTags[idecay - 1]->barcode())
0508       {
0509         return;
0510       }
0511     }
0512
0513     mcTag = mcTags[0];
0514   }
0515   else
0516   {
0517     SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0518     if(!trackmap) return;
0519
0520     for (int idecay = 0; idecay < m_nDaughters; idecay++)
0521     {
0522       SvtxTrack *track = trackmap->get(decays[idecay]);
0523       if(!track) return;
0524       g4particle = m_trackeval->max_truth_particle_by_nclusters(track);
0525
0526       if(!g4particle)
0527       {
0528         return;
0529       }
0530
0531       daughtersID[idecay] = g4particle->get_barcode();
0532     }
0533   }
0534
0535
0536
0537
0538
0539   for(auto mcJet : *m_truth_taggedJetContainer)
0540   {
0541     for(auto comp : mcJet->get_comp_vec())
0542     {
0543       if(m_stable_mother)
0544       {
0545         if(comp.first == Jet::SRC::VOID && int(comp.second) == mcTag->barcode())
0546         {
0547           mcTagJet = mcJet;
0548         }
0549       }
0550       else
0551       {
0552         if(comp.first == Jet::SRC::VOID)
0553         {
0554           HepMC::GenParticle* TheMother = hepMCGenEvent->barcode_to_particle(int(comp.second));
0555           HepMC::GenVertex *TagVertex = TheMother->end_vertex();
0556           int nMatches = 0;
0557           //if HEPMC vertex exists
0558           if(TagVertex)
0559           {
0560             for (HepMC::GenVertex::particle_iterator it = TagVertex->particles_begin(HepMC::descendants); it != TagVertex->particles_end(HepMC::descendants); ++it)
0561             {
0562               for(unsigned int idau = 0; idau < daughtersID.size(); idau++)
0563               {
0564                 if(daughtersID[idau] == (*it)->barcode())
0565                 {
0566                   nMatches++;
0567                   daughtersID[idau] = -99999;
0568                 }
0569               }
0570             }
0571           }
0572           if(nMatches == m_nDaughters)
0573           {
0574             mcTag = TheMother;
0575             mcTagJet = mcJet;
0576           }
0577         }
0578       }
0579     }
0580
0581     if(mcTagJet)
0582     {
0583       break;
0584     }
0585   }
0586   return;
0587 }
0588
0589 HepMC::GenParticle *BuildResonanceJetTaggingTree::getMother(PHCompositeNode *topNode, PHG4Particle *g4daughter)
0590 {
0591   PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0592
0593   /// If the node was not properly put on the tree, return
0594   if (!hepmceventmap)
0595   {
0596     std::cout << PHWHERE
0597               << "HEPMC event map node is missing, can't collected HEPMC truth particles"
0598               << std::endl;
0599     return nullptr;
0600   }
0601
0602   PHHepMCGenEvent *hepmcevent = hepmceventmap->get(1);
0603   if (!hepmcevent)
0604   {
0605     std::cout << "no hepmcevent!!!" << std::endl;
0606     return nullptr;
0607   }
0608
0609   HepMC::GenEvent *hepMCGenEvent = hepmcevent->getEvent();
0610   if (!hepMCGenEvent)
0611   {
0612     return nullptr;
0613   }
0614
0615   HepMC::GenParticle *mcDaughter = nullptr;
0616
0617   if (g4daughter->get_barcode() > 0)
0618   {
0619     mcDaughter = hepMCGenEvent->barcode_to_particle(g4daughter->get_barcode());
0620   }
0621   if (!mcDaughter)
0622   {
0623     return nullptr;
0624   }
0625
0626   HepMC::GenVertex *TagVertex = mcDaughter->production_vertex();
0627   for (HepMC::GenVertex::particle_iterator it = TagVertex->particles_begin(HepMC::ancestors); it != TagVertex->particles_end(HepMC::ancestors); ++it)
0628   {
0629     if (std::abs((*it)->pdg_id()) == m_tag_pdg)
0630     {
0631       return (*it);
0632     }
0633   }
0634
0635   return nullptr;
0636 }
0637
0638 bool BuildResonanceJetTaggingTree::isReconstructed(int index, std::vector<int> indexRecVector)
0639 {
0640   for(auto indexRec : indexRecVector)
0641   {
0642     if(index == indexRec) return true;
0643   }
0644   return false;
0645 }
0646
0647 void BuildResonanceJetTaggingTree::initializeTrees()
0648 {
0649   delete m_runinfo;
0650   m_runinfo = new TTree("m_runinfo","A tree with the run information");
0651   m_runinfo->Branch("m_intlumi",&m_intlumi, "m_intlumi/F");
0652   m_runinfo->Branch("m_nprocessedevents", &m_nprocessedevents, "m_nprocessedevents/F");
0653   m_runinfo->Branch("m_nacceptedevents", &m_nacceptedevents, "m_nacceptedevents/F");
0654   m_runinfo->Branch("m_xsecprocessedevents", &m_xsecprocessedevents, "m_xsecprocessedevents/F");
0655   m_runinfo->Branch("m_xsecacceptedevents", &m_xsecacceptedevents, "m_xsecacceptedevents/F");
0656
0657   delete m_taggedjettree;
0658   m_taggedjettree = new TTree("m_taggedjettree", "A tree with Tagged-Jet info");
0659   m_taggedjettree->Branch("m_reco_tag_px", &m_reco_tag_px, "m_reco_tag_px/F");
0660   m_taggedjettree->Branch("m_reco_tag_py", &m_reco_tag_py, "m_reco_tag_py/F");
0661   m_taggedjettree->Branch("m_reco_tag_pz", &m_reco_tag_pz, "m_reco_tag_pz/F");
0662   m_taggedjettree->Branch("m_reco_tag_pt", &m_reco_tag_pt, "m_reco_tag_pt/F");
0663   m_taggedjettree->Branch("m_reco_tag_eta", &m_reco_tag_eta, "m_reco_tag_eta/F");
0664   m_taggedjettree->Branch("m_reco_tag_phi", &m_reco_tag_phi, "m_reco_tag_phi/F");
0665   m_taggedjettree->Branch("m_reco_tag_m", &m_reco_tag_m, "m_reco_tag_m/F");
0666   m_taggedjettree->Branch("m_reco_tag_e", &m_reco_tag_e, "m_reco_tag_e/F");
0667   m_taggedjettree->Branch("m_reco_jet_px", &m_reco_jet_px, "m_reco_jet_px/F");
0668   m_taggedjettree->Branch("m_reco_jet_py", &m_reco_jet_py, "m_reco_jet_py/F");
0669   m_taggedjettree->Branch("m_reco_jet_pz", &m_reco_jet_pz, "m_reco_jet_pz/F");
0670   m_taggedjettree->Branch("m_reco_jet_pt", &m_reco_jet_pt, "m_reco_jet_pt/F");
0671   m_taggedjettree->Branch("m_reco_jet_eta", &m_reco_jet_eta, "m_reco_jet_eta/F");
0672   m_taggedjettree->Branch("m_reco_jet_phi", &m_reco_jet_phi, "m_reco_jet_phi/F");
0673   m_taggedjettree->Branch("m_reco_jet_m", &m_reco_jet_m, "m_reco_jet_m/F");
0674   m_taggedjettree->Branch("m_reco_jet_e", &m_reco_jet_e, "m_reco_jet_e/F");
0675
0676   m_taggedjettree->Branch("m_truth_tag_px", &m_truth_tag_px, "m_truth_tag_px/F");
0677   m_taggedjettree->Branch("m_truth_tag_py", &m_truth_tag_py, "m_truth_tag_py/F");
0678   m_taggedjettree->Branch("m_truth_tag_pz", &m_truth_tag_pz, "m_truth_tag_pz/F");
0679   m_taggedjettree->Branch("m_truth_tag_pt", &m_truth_tag_pt, "m_truth_tag_pt/F");
0680   m_taggedjettree->Branch("m_truth_tag_eta", &m_truth_tag_eta, "m_truth_tag_eta/F");
0681   m_taggedjettree->Branch("m_truth_tag_phi", &m_truth_tag_phi, "m_truth_tag_phi/F");
0682   m_taggedjettree->Branch("m_truth_tag_m", &m_truth_tag_m, "m_truth_tag_m/F");
0683   m_taggedjettree->Branch("m_truth_tag_e", &m_truth_tag_e, "m_truth_tag_e/F");
0684   m_taggedjettree->Branch("m_truth_jet_px", &m_truth_jet_px, "m_truth_jet_px/F");
0685   m_taggedjettree->Branch("m_truth_jet_py", &m_truth_jet_py, "m_truth_jet_py/F");
0686   m_taggedjettree->Branch("m_truth_jet_pz", &m_truth_jet_pz, "m_truth_jet_pz/F");
0687   m_taggedjettree->Branch("m_truth_jet_pt", &m_truth_jet_pt, "m_truth_jet_pt/F");
0688   m_taggedjettree->Branch("m_truth_jet_eta", &m_truth_jet_eta, "m_truth_jet_eta/F");
0689   m_taggedjettree->Branch("m_truth_jet_phi", &m_truth_jet_phi, "m_truth_jet_phi/F");
0690   m_taggedjettree->Branch("m_truth_jet_m", &m_truth_jet_m, "m_truth_jet_m/F");
0691   m_taggedjettree->Branch("m_truth_jet_e", &m_truth_jet_e, "m_truth_jet_e/F");
0692   m_taggedjettree->Branch("m_truthjet_const_px", &m_truthjet_const_px);
0693   m_taggedjettree->Branch("m_truthjet_const_py", &m_truthjet_const_py);
0694   m_taggedjettree->Branch("m_truthjet_const_pz", &m_truthjet_const_pz);
0695   m_taggedjettree->Branch("m_truthjet_const_e", &m_truthjet_const_e);
0696
0697 }
0698 void BuildResonanceJetTaggingTree::initializeVariables()
0699 {
0700   delete m_outfile;
0701   m_outfile = new TFile(m_outfilename.c_str(), "RECREATE");
0702
0703   delete m_eventcount_h;
0704   m_eventcount_h = new TH1I("eventcount_h", "Event Count", 2, -0.5, 1.5);
0705   m_eventcount_h->GetXaxis()->SetBinLabel(1,"N raw ev anal");
0706   m_eventcount_h->GetXaxis()->SetBinLabel(2,"N ev anal");
0707 }
0708
0709 void BuildResonanceJetTaggingTree::resetTreeVariables()
0710 {
0711   // Tagged-Jet reconstructed variables
0712   m_reco_tag_px = NAN;
0713   m_reco_tag_py = NAN;
0714   m_reco_tag_pz = NAN;
0715   m_reco_tag_pt = NAN;
0716   m_reco_tag_eta = NAN;
0717   m_reco_tag_phi = NAN;
0718   m_reco_tag_m = NAN;
0719   m_reco_tag_e = NAN;
0720   m_reco_jet_px = NAN;
0721   m_reco_jet_py = NAN;
0722   m_reco_jet_pz = NAN;
0723   m_reco_jet_pt = NAN;
0724   m_reco_jet_eta = NAN;
0725   m_reco_jet_phi = NAN;
0726   m_reco_jet_m = NAN;
0727   m_reco_jet_e = NAN;
0728   //Truth info
0729   m_truth_tag_px = NAN;
0730   m_truth_tag_py = NAN;
0731   m_truth_tag_pz = NAN;
0732   m_truth_tag_pt = NAN;
0733   m_truth_tag_eta = NAN;
0734   m_truth_tag_phi = NAN;
0735   m_truth_tag_m = NAN;
0736   m_truth_tag_e = NAN;
0737   m_truth_jet_px = NAN;
0738   m_truth_jet_py = NAN;
0739   m_truth_jet_pz = NAN;
0740   m_truth_jet_pt = NAN;
0741   m_truth_jet_eta = NAN;
0742   m_truth_jet_phi = NAN;
0743   m_truth_jet_m = NAN;
0744   m_truth_jet_e = NAN;
0745
0746   m_truthjet_const_px.clear();
0747   m_truthjet_const_py.clear();
0748   m_truthjet_const_pz.clear();
0749   m_truthjet_const_e.clear();
0750 }
0751
0752 JetContainerv1* BuildResonanceJetTaggingTree::getJetContainerFromNode(PHCompositeNode *topNode, const std::string &name)
0753 {
0754   JetContainerv1 *jetcont = findNode::getClass<JetContainerv1>(topNode, name);
0755
0756   if (!jetcont)
0757   {
0758     std::cout << PHWHERE
0759          << "JetContainerv1 node is missing, can't collect jets"
0760          << std::endl;
0761     return 0;
0762   }
0763
0764   return jetcont;
0765 }
0766 KFParticle_Container* BuildResonanceJetTaggingTree::getKFParticleContainerFromNode(PHCompositeNode *topNode, const std::string &name)
0767 {
0768   KFParticle_Container *cont = findNode::getClass<KFParticle_Container>(topNode, name);
0769
0770   if (!cont)
0771   {
0772     std::cout << PHWHERE
0773          << "KFParticle_Container node is missing, can't collect HF particles"
0774          << std::endl;
0775     return 0;
0776   }
0777
0778   return cont;
0779 }
0780 HepMC::GenEvent* BuildResonanceJetTaggingTree::getGenEventFromNode(PHCompositeNode *topNode, const std::string &name)
0781 {
0782   PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, name);
0783
0784   /// If the node was not properly put on the tree, return
0785   if (!hepmceventmap)
0786   {
0787     std::cout << PHWHERE
0788          << "HEPMC event map node is missing, can't collected HEPMC truth particles"
0789          << std::endl;
0790     return 0;
0791   }
0792
0793   PHHepMCGenEvent *hepmcevent = hepmceventmap->get(1);
0794
0795   if (!hepmcevent)
0796   {
0797     std::cout << PHWHERE
0798          << "PHHepMCGenEvent node is missing, can't collected HEPMC truth particles"
0799          << std::endl;
0800     return 0;
0801   }
0802
0803   HepMC::GenEvent *hepMCGenEvent = hepmcevent->getEvent();
0804
0805   if (!hepmceventmap)
0806   {
0807     std::cout << PHWHERE
0808          << "HepMC::GenEvent node is missing, can't collected HEPMC truth particles"
0809          << std::endl;
0810     return 0;
0811   }
0812
0813   return hepMCGenEvent;
0814
0815 }