HF-Jet/TruthGeneration/HFJetTruthTrigger.cc

0001 #include "HFJetTruthTrigger.h"
0002
0003 #include "HFJetDefs.h"
0004
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006 #include <fun4all/Fun4AllServer.h>
0007 #include <phool/getClass.h>
0008
0009 #include <phool/PHCompositeNode.h>
0010
0011 #include <g4main/PHG4Hit.h>
0012 #include <g4main/PHG4Particle.h>
0013 #include <g4main/PHG4TruthInfoContainer.h>
0014
0015 #include <TDatabasePDG.h>
0016 #include <TFile.h>
0017 #include <TH2D.h>
0018 #include <TLorentzVector.h>
0019 #include <TString.h>
0020 #include <TTree.h>
0021 #include <g4jets/Jet.h>
0022 #include <g4jets/JetMap.h>
0023 #pragma GCC diagnostic push
0024 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0025 #include <HepMC/GenEvent.h>
0026 #include <HepMC/GenVertex.h>
0027 #include <phhepmc/PHHepMCGenEvent.h>
0028 #include <phhepmc/PHHepMCGenEventMap.h>
0029 #pragma GCC diagnostic pop
0030 #include <cmath>
0031 #include <cstddef>
0032 #include <iostream>
0033
0034 HFJetTruthTrigger::HFJetTruthTrigger(std::string filename, int flavor,
0035                                      std::string jet_node, int maxevent)
0036   : SubsysReco("HFJetTagger_" + jet_node)
0037   , _verbose(0)
0038   , _ievent(0)
0039   , _total_pass(0)
0040   , _f(nullptr)
0041   , _h2(nullptr)
0042   , _h2all(nullptr)
0043   , _h2_b(nullptr)
0044   , _h2_c(nullptr)
0045   , _embedding_id(1)
0046 {
0047   _foutname = filename;
0048
0049   _flavor = flavor;
0050
0051   _maxevent = maxevent;
0052
0053   _pt_min = 20;
0054   _pt_max = 100;
0055
0056   _eta_min = -.7;
0057   _eta_max = +.7;
0058
0059   _jet_name = jet_node;
0060
0061   _rejection_action = Fun4AllReturnCodes::DISCARDEVENT;
0062 }
0063
0064 int HFJetTruthTrigger::Init(PHCompositeNode* topNode)
0065 {
0066   _verbose = true;
0067
0068   _ievent = 0;
0069
0070   _total_pass = 0;
0071
0072   _f = new TFile(_foutname.c_str(), "RECREATE");
0073
0074   _h2 = new TH2D("h2", "", 100, 0, 100.0, 40, -2, +2);
0075   _h2_b = new TH2D("h2_b", "", 100, 0, 100.0, 40, -2, +2);
0076   _h2_c = new TH2D("h2_c", "", 100, 0, 100.0, 40, -2, +2);
0077   _h2all = new TH2D("h2all", "", 100, 0, 100.0, 40, -2, +2);
0078
0079   return 0;
0080 }
0081
0082 int HFJetTruthTrigger::process_event(PHCompositeNode* topNode)
0083 {
0084   PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0085   if (!geneventmap)
0086   {
0087     std::cout << PHWHERE << " - Fatal error - missing node PHHepMCGenEventMap" << std::endl;
0088     return Fun4AllReturnCodes::ABORTRUN;
0089   }
0090
0091   PHHepMCGenEvent* genevt = geneventmap->get(_embedding_id);
0092   if (!genevt)
0093   {
0094     std::cout << PHWHERE << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID " << _embedding_id;
0095     std::cout << ". Print PHHepMCGenEventMap:";
0096     geneventmap->identify();
0097     return Fun4AllReturnCodes::ABORTRUN;
0098   }
0099
0100   HepMC::GenEvent* theEvent = genevt->getEvent();
0101   //theEvent->print();
0102
0103   JetMap* truth_jets = findNode::getClass<JetMap>(topNode, _jet_name);
0104   if (!truth_jets)
0105   {
0106     std::cout << PHWHERE << " - Fatal error - node " << _jet_name << " JetMap missing." << std::endl;
0107     return Fun4AllReturnCodes::ABORTRUN;
0108   }
0109   const double jet_radius = truth_jets->get_par();
0110
0111   if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0112     std::cout << __PRETTY_FUNCTION__ << ": truth jets has size "
0113               << truth_jets->size() << " and R = " << jet_radius << std::endl;
0114
0115   int ijet_t = 0;
0116   bool pass_event = false;
0117   for (JetMap::Iter iter = truth_jets->begin(); iter != truth_jets->end();
0118        ++iter)
0119   {
0120     Jet* this_jet = iter->second;
0121
0122     float this_pt = this_jet->get_pt();
0123     float this_eta = this_jet->get_eta();
0124
0125     if (this_pt < 10 || fabs(this_eta) > 5)
0126       continue;
0127
0128     _h2all->Fill(this_jet->get_pt(), this_eta);
0129
0130     if (this_pt > _pt_min && this_pt < _pt_max && (this_eta) > _eta_min && (this_eta) < _eta_max)
0131     {
0132       //pass_event = true;
0133       _h2->Fill(this_jet->get_pt(), this_eta);
0134       if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0135         this_jet->identify();
0136     }
0137     else
0138     {
0139       continue;
0140     }
0141
0142     const int jet_flavor = parton_tagging(this_jet, theEvent, jet_radius);
0143
0144     if (abs(jet_flavor) == _flavor)
0145     {
0146       pass_event = true;
0147       if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0148       {
0149         this_jet->identify();
0150         std::cout << " --> this is flavor " << jet_flavor
0151                   << " like I want " << std::endl;
0152       }
0153     }
0154     else
0155     {
0156       if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0157       {
0158         this_jet->identify();
0159         std::cout << " --> this is flavor " << jet_flavor
0160                   << " which I do NOT want " << std::endl;
0161       }
0162     }
0163
0164     hadron_tagging(this_jet, theEvent, jet_radius);
0165
0166     ijet_t++;
0167   }
0168
0169   if (pass_event && _total_pass >= _maxevent)
0170   {
0171     if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_SOME)
0172       std::cout << __PRETTY_FUNCTION__ << " --> FAIL due to max events = "
0173                 << _total_pass << std::endl;
0174     _ievent++;
0175     return _rejection_action;
0176   }
0177   else if (pass_event)
0178   {
0179     _total_pass++;
0180     if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_SOME)
0181       std::cout << __PRETTY_FUNCTION__ << " --> PASS, total = " << _total_pass
0182                 << std::endl;
0183     _ievent++;
0184     return Fun4AllReturnCodes::EVENT_OK;
0185   }
0186   else
0187   {
0188     if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_SOME)
0189       std::cout << __PRETTY_FUNCTION__ << " --> FAIL " << std::endl;
0190     _ievent++;
0191     return _rejection_action;
0192   }
0193 }
0194
0195 int HFJetTruthTrigger::End(PHCompositeNode* topNode)
0196 {
0197   if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_SOME)
0198     std::cout << __PRETTY_FUNCTION__ << " DVP PASSED " << _total_pass
0199               << " events" << std::endl;
0200
0201   _f->cd();
0202   _f->Write();
0203   //_f->Close();
0204
0205   return 0;
0206 }
0207
0208 //! tag jet flavor by parton matching, like PRL 113, 132301 (2014)
0209 int HFJetTruthTrigger::parton_tagging(Jet* this_jet, HepMC::GenEvent* theEvent,
0210                                       const double match_radius)
0211 {
0212   float this_pt = this_jet->get_pt();
0213   float this_phi = this_jet->get_phi();
0214   float this_eta = this_jet->get_eta();
0215
0216   int jet_flavor = 0;
0217   double jet_parton_zt = 0;
0218
0219   //std::cout << " truth jet #" << ijet_t << ", pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << ", checking flavor" << std::endl;
0220
0221   //TODO: lack taggign scheme of gluon splitting -> QQ_bar
0222   for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin();
0223        p != theEvent->particles_end(); ++p)
0224   {
0225     float dR = deltaR((*p)->momentum().pseudoRapidity(), this_eta,
0226                       (*p)->momentum().phi(), this_phi);
0227     if (dR > match_radius)
0228       continue;
0229
0230     int pidabs = abs((*p)->pdg_id());
0231     const double zt = (*p)->momentum().perp() / this_pt;
0232
0233     if (pidabs == TDatabasePDG::Instance()->GetParticle("c")->PdgCode()      //
0234         or pidabs == TDatabasePDG::Instance()->GetParticle("b")->PdgCode()   //
0235         or pidabs == TDatabasePDG::Instance()->GetParticle("t")->PdgCode())  // handle heavy quarks only. All other favor tagged as default 0
0236     {
0237       if (pidabs > abs(jet_flavor))  // heavy quark found
0238       {
0239         jet_parton_zt = zt;
0240         jet_flavor = (*p)->pdg_id();
0241       }
0242       else if (pidabs == abs(jet_flavor))  // same quark mass. next compare zt
0243       {
0244         if (zt > jet_parton_zt)
0245         {
0246           jet_parton_zt = zt;
0247           jet_flavor = (*p)->pdg_id();
0248         }
0249       }
0250
0251       if (pidabs == TDatabasePDG::Instance()->GetParticle("b")->PdgCode())
0252       {
0253         if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0254           std::cout << __PRETTY_FUNCTION__
0255                     << " --BOTTOM--> pt / eta / phi = "
0256                     << (*p)->momentum().perp() << " / "
0257                     << (*p)->momentum().pseudoRapidity() << " / "
0258                     << (*p)->momentum().phi() << std::endl;
0259       }
0260       else if (pidabs == TDatabasePDG::Instance()->GetParticle("c")->PdgCode())
0261       {
0262         if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0263           std::cout << __PRETTY_FUNCTION__
0264                     << " --CHARM --> pt / eta / phi = "
0265                     << (*p)->momentum().perp() << " / "
0266                     << (*p)->momentum().pseudoRapidity() << " / "
0267                     << (*p)->momentum().phi() << std::endl;
0268       }
0269     }
0270   }  //       for (HepMC::GenEvent::particle_const_iterator p =
0271
0272   if (abs(jet_flavor) == TDatabasePDG::Instance()->GetParticle("b")->PdgCode())
0273   {
0274     _h2_b->Fill(this_jet->get_pt(), this_eta);
0275   }
0276   else if (abs(jet_flavor) == TDatabasePDG::Instance()->GetParticle("c")->PdgCode())
0277   {
0278     _h2_c->Fill(this_jet->get_pt(), this_eta);
0279   }
0280
0281   this_jet->set_property(static_cast<Jet::PROPERTY>(prop_JetPartonFlavor),
0282                          jet_flavor);
0283   this_jet->set_property(static_cast<Jet::PROPERTY>(prop_JetPartonZT),
0284                          jet_parton_zt);
0285   //          this_jet->identify();
0286
0287   return jet_flavor;
0288 }
0289
0290 int HFJetTruthTrigger::hadron_tagging(Jet* this_jet, HepMC::GenEvent* theEvent,
0291                                       const double match_radius)
0292 {
0293   float this_pt = this_jet->get_pt();
0294   float this_phi = this_jet->get_phi();
0295   float this_eta = this_jet->get_eta();
0296
0297   int jet_flavor = 0;
0298   double jet_parton_zt = 0;
0299
0300   //std::cout << " truth jet #" << ijet_t << ", pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << ", checking flavor" << std::endl;
0301
0302   for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin();
0303        p != theEvent->particles_end(); ++p)
0304   {
0305     float dR = deltaR((*p)->momentum().pseudoRapidity(), this_eta,
0306                       (*p)->momentum().phi(), this_phi);
0307     if (dR > match_radius)
0308       continue;
0309
0310     int pidabs = 0;
0311     TParticlePDG* pdg_p = TDatabasePDG::Instance()->GetParticle((*p)->pdg_id());
0312     if (pdg_p)
0313     {
0314       if (TString(pdg_p->ParticleClass()).BeginsWith("B-"))
0315       {
0316         pidabs = TDatabasePDG::Instance()->GetParticle("b")->PdgCode();
0317       }
0318       else if (TString(pdg_p->ParticleClass()).BeginsWith("Charmed"))
0319       {
0320         pidabs = TDatabasePDG::Instance()->GetParticle("c")->PdgCode();
0321       }
0322     }
0323
0324     const double zt = (*p)->momentum().perp() / this_pt;
0325
0326     if (pidabs == TDatabasePDG::Instance()->GetParticle("c")->PdgCode()      //
0327         or pidabs == TDatabasePDG::Instance()->GetParticle("b")->PdgCode())  // handle heavy quarks only. All other favor tagged as default 0
0328     {
0329       if (pidabs > abs(jet_flavor))  // heavy quark found
0330       {
0331         jet_parton_zt = zt;
0332         jet_flavor = pidabs;
0333       }
0334       else if (pidabs == abs(jet_flavor))  // same quark mass. next compare zt
0335       {
0336         if (zt > jet_parton_zt)
0337         {
0338           jet_parton_zt = zt;
0339           jet_flavor = pidabs;
0340         }
0341       }
0342
0343       if (pidabs == TDatabasePDG::Instance()->GetParticle("b")->PdgCode())
0344       {
0345         if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0346         {
0347           std::cout << __PRETTY_FUNCTION__
0348                     << " --BOTTOM--> pt / eta / phi = "
0349                     << (*p)->momentum().perp() << " / "
0350                     << (*p)->momentum().pseudoRapidity() << " / "
0351                     << (*p)->momentum().phi() << " with hadron ";
0352           pdg_p->Print();
0353         }
0354       }
0355       else if (pidabs == TDatabasePDG::Instance()->GetParticle("c")->PdgCode())
0356       {
0357         if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0358         {
0359           std::cout << __PRETTY_FUNCTION__
0360                     << " --CHARM --> pt / eta / phi = "
0361                     << (*p)->momentum().perp() << " / "
0362                     << (*p)->momentum().pseudoRapidity() << " / "
0363                     << (*p)->momentum().phi() << " with hadron ";
0364           pdg_p->Print();
0365         }
0366       }
0367     }
0368   }  //       for (HepMC::GenEvent::particle_const_iterator p =
0369
0370   this_jet->set_property(static_cast<Jet::PROPERTY>(prop_JetHadronFlavor),
0371                          jet_flavor);
0372   this_jet->set_property(static_cast<Jet::PROPERTY>(prop_JetHadronZT),
0373                          jet_parton_zt);
0374   //          this_jet->identify();
0375
0376   if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
0377     std::cout << __PRETTY_FUNCTION__ << " jet_flavor = " << jet_flavor
0378               << std::endl;
0379
0380   return jet_flavor;
0381 }