TruthJetTagging/src/TruthJetTagging.cc

0001 //____________________________________________________________________________..
0002
0003 #include "TruthJetTagging.h"
0004
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006 #include <fun4all/Fun4AllServer.h>
0007 #include <fun4all/PHTFileServer.h>
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/PHIODataNode.h>    // for PHIODataNode
0010 #include <phool/PHNode.h>          // for PHNode
0011 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0012 #include <phool/PHObject.h>        // for PHObject
0013 #include <phool/getClass.h>
0014 #include <phool/phool.h>                    // for PHWHERE
0015 #include <g4main/PHG4Utils.h>
0016
0017 #pragma GCC diagnostic push
0018 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0019 #include <HepMC/GenEvent.h>
0020 #include <HepMC/GenVertex.h>
0021 #pragma GCC diagnostic pop
0022
0023 #include <phhepmc/PHHepMCGenEvent.h>
0024 #include <phhepmc/PHHepMCGenEventMap.h>
0025 #include <TDatabasePDG.h>
0026
0027 #include <g4main/PHG4Particle.h>
0028 #include <g4main/PHG4TruthInfoContainer.h>
0029 #include <TMath.h>
0030
0031 #include <g4jets/JetMap.h>
0032 #include <g4jets/Jet.h>
0033 #include <g4jets/Jetv1.h>
0034 #include <sstream>
0035
0036 //____________________________________________________________________________..
0037 TruthJetTagging::TruthJetTagging(const std::string &name):
0038  SubsysReco(name)
0039  , _algorithms()
0040  , _radii()
0041  , _do_photon_tagging()
0042  , _do_hf_tagging()
0043  , _embedding_id(1)
0044 {
0045
0046 }
0047
0048 //____________________________________________________________________________..
0049 TruthJetTagging::~TruthJetTagging()
0050 {
0051
0052 }
0053
0054 //____________________________________________________________________________..
0055 int TruthJetTagging::Init(PHCompositeNode *topNode)
0056 {
0057   return Fun4AllReturnCodes::EVENT_OK;
0058 }
0059
0060 //____________________________________________________________________________..
0061 int TruthJetTagging::InitRun(PHCompositeNode *topNode)
0062 {
0063   return Fun4AllReturnCodes::EVENT_OK;
0064 }
0065
0066 //____________________________________________________________________________..
0067 int TruthJetTagging::process_event(PHCompositeNode *topNode)
0068 {
0069 PHG4TruthInfoContainer* truthcontainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0070 if (!truthcontainer)
0071   {
0072     std::cout
0073       << "MyJetAnalysis::process_event - Error can not find DST truthcontainer node "
0074         << std::endl;
0075 exit(-1);
0076 }
0077   PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0078   if (!geneventmap)
0079   {
0080     std::cout << PHWHERE << " - Fatal error - missing node PHHepMCGenEventMap" << std::endl;
0081     return Fun4AllReturnCodes::ABORTRUN;
0082   }
0083
0084   PHHepMCGenEvent* genevt = geneventmap->get(_embedding_id);
0085   if (!genevt)
0086   {
0087     std::cout << PHWHERE << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID " << _embedding_id;
0088     std::cout << ". Print PHHepMCGenEventMap:";
0089     geneventmap->identify();
0090     return Fun4AllReturnCodes::ABORTRUN;
0091   }
0092
0093
0094
0095 HepMC::GenEvent* theEvent = genevt->getEvent();
0096
0097   int n_algo = _algorithms.size();
0098   int n_radii = _radii.size();
0099   if (_do_hf_tagging)
0100     {
0101       if (n_radii != n_algo)
0102     {
0103       std::cout << "TruthJetTagging::process_event - errorr unequal number of jet radii and algoirthms specified" << std::endl;
0104       exit(-1);
0105     }
0106     }
0107   for (int algoiter = 0;algoiter < n_algo;algoiter++)
0108     {
0109       JetMap* tjets = findNode::getClass<JetMap>(topNode, _algorithms.at(algoiter).c_str());
0110       if (!tjets)
0111     {
0112       std::cout
0113         << "MyJetAnalysis::process_event - Error can not find DST JetMap node "
0114         << std::endl;
0115       exit(-1);
0116     }
0117       for (JetMap::Iter iter = tjets->begin(); iter != tjets->end(); ++iter)
0118     {
0119           Jet* tjet = iter->second;
0120       assert(tjet);
0121
0122       if (_do_photon_tagging)
0123         {
0124           float photon_fraction =  TruthJetTagging::TruthPhotonTagging ( truthcontainer , tjet );
0125           tjet->set_property(Jet::prop_gamma,photon_fraction);
0126         }
0127       if (_do_hf_tagging)
0128         {
0129
0130           float jet_radius = 0.4;
0131           int jet_flavor = TruthJetTagging::hadron_tagging(tjet, theEvent, jet_radius);
0132           tjet->set_property(Jet::prop_JetHadronFlavor,jet_flavor);
0133         }
0134
0135     }
0136     }
0137   return Fun4AllReturnCodes::EVENT_OK;
0138 }
0139
0140   float TruthJetTagging::TruthPhotonTagging ( PHG4TruthInfoContainer* truthnode, Jet* tjet)
0141   {
0142     assert(truthnode);
0143     assert(tjet);
0144     float jetpt = tjet->get_pt();
0145     float max_gamma_pt = 0;
0146     for (Jet::ConstIter comp = tjet->begin_comp(); comp != tjet->end_comp(); ++comp)
0147       {
0148     PHG4Particle* particle = truthnode->GetParticle((*comp).second);
0149     if (particle->get_pid() == 22)
0150       {
0151         float particle_pt = TMath::Sqrt(TMath::Power(particle->get_px(),2) + TMath::Power(particle->get_py(),2));
0152         if (particle_pt > max_gamma_pt)
0153           {
0154         max_gamma_pt = particle_pt;
0155           }
0156       }
0157       }
0158     float ratio = max_gamma_pt/jetpt;
0159     return ratio;
0160   }
0161
0162
0163
0164
0165
0166
0167 int TruthJetTagging::hadron_tagging(Jet* this_jet, HepMC::GenEvent* theEvent,
0168                                       const double match_radius)
0169 {
0170   float this_pt = this_jet->get_pt();
0171   float this_phi = this_jet->get_phi();
0172   float this_eta = this_jet->get_eta();
0173
0174   int jet_flavor = 0;
0175   double jet_parton_zt = 0;
0176
0177   for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin();
0178        p != theEvent->particles_end(); ++p)
0179   {
0180     float dR = deltaR((*p)->momentum().pseudoRapidity(), this_eta,
0181                       (*p)->momentum().phi(), this_phi);
0182     if (dR > match_radius)
0183       continue;
0184
0185     int pidabs = 0;
0186     TParticlePDG* pdg_p = TDatabasePDG::Instance()->GetParticle((*p)->pdg_id());
0187     if (pdg_p)
0188     {
0189       if (TString(pdg_p->ParticleClass()).BeginsWith("B-"))
0190       {
0191         pidabs = TDatabasePDG::Instance()->GetParticle("b")->PdgCode();
0192       }
0193       else if (TString(pdg_p->ParticleClass()).BeginsWith("Charmed"))
0194       {
0195         pidabs = TDatabasePDG::Instance()->GetParticle("c")->PdgCode();
0196       }
0197     }
0198
0199     const double zt = (*p)->momentum().perp() / this_pt;
0200
0201     if (pidabs == TDatabasePDG::Instance()->GetParticle("c")->PdgCode()      //
0202         or pidabs == TDatabasePDG::Instance()->GetParticle("b")->PdgCode())  // handle heavy quarks only. All other favor tagged as default 0
0203     {
0204       if (pidabs > abs(jet_flavor))  // heavy quark found
0205       {
0206         jet_parton_zt = zt;
0207         jet_flavor = pidabs;
0208       }
0209       else if (pidabs == abs(jet_flavor))  // same quark mass. next compare zt
0210       {
0211         if (zt > jet_parton_zt)
0212         {
0213           jet_parton_zt = zt;
0214           jet_flavor = pidabs;
0215         }
0216       }
0217
0218       // if (pidabs == TDatabasePDG::Instance()->GetParticle("b")->PdgCode())
0219       // {
0220       //   // if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0221       //   // {
0222       //   //   std::cout << __PRETTY_FUNCTION__
0223       //   //             << " --BOTTOM--> pt / eta / phi = "
0224       //   //             << (*p)->momentum().perp() << " / "
0225       //   //             << (*p)->momentum().pseudoRapidity() << " / "
0226       //   //             << (*p)->momentum().phi() << " with hadron ";
0227       //   //   pdg_p->Print();
0228       //   // }
0229       // }
0230       // else if (pidabs == TDatabasePDG::Instance()->GetParticle("c")->PdgCode())
0231       // {
0232       //   // if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0233       //   // {
0234       //   //   std::cout << __PRETTY_FUNCTION__
0235       //   //             << " --CHARM --> pt / eta / phi = "
0236       //   //             << (*p)->momentum().perp() << " / "
0237       //   //             << (*p)->momentum().pseudoRapidity() << " / "
0238       //   //             << (*p)->momentum().phi() << " with hadron ";
0239       //   //   pdg_p->Print();
0240       //   // }
0241       // }
0242     }
0243   }  //       for (HepMC::GenEvent::particle_const_iterator p =
0244
0245
0246   this_jet->set_property(Jet::prop_JetHadronZT,
0247              jet_parton_zt);
0248   // this_jet->identify();
0249
0250   // if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0251   //   std::cout << __PRETTY_FUNCTION__ << " jet_flavor = " << jet_flavor
0252   //             << std::endl;
0253
0254   return jet_flavor;
0255 }
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277 //____________________________________________________________________________..
0278 int TruthJetTagging::ResetEvent(PHCompositeNode *topNode)
0279 {
0280   return Fun4AllReturnCodes::EVENT_OK;
0281 }
0282
0283 //____________________________________________________________________________..
0284 int TruthJetTagging::EndRun(const int runnumber)
0285 {
0286   return Fun4AllReturnCodes::EVENT_OK;
0287 }
0288
0289 //____________________________________________________________________________..
0290 int TruthJetTagging::End(PHCompositeNode *topNode)
0291 {
0292   return Fun4AllReturnCodes::EVENT_OK;
0293 }
0294
0295 //____________________________________________________________________________..
0296 int TruthJetTagging::Reset(PHCompositeNode *topNode)
0297 {
0298   return Fun4AllReturnCodes::EVENT_OK;
0299 }
0300
0301 //____________________________________________________________________________..
0302 void TruthJetTagging::Print(const std::string &what) const
0303 {
0304   std::cout << "TruthJetTagging::Print(const std::string &what) const Printing info for " << what << std::endl;
0305 }