Core/src/EmbeddingHepMCFilter.cc

0001 #include "GeneratorInterface/Core/interface/EmbeddingHepMCFilter.h"
0002
0003 #include "boost/algorithm/string.hpp"
0004 #include "boost/algorithm/string/trim_all.hpp"
0005
0006 EmbeddingHepMCFilter::EmbeddingHepMCFilter(const edm::ParameterSet &iConfig)
0007     : ZPDGID_(iConfig.getParameter<int>("BosonPDGID")),
0008       includeDY_(iConfig.existsAs<bool>("IncludeDY") ? iConfig.getParameter<bool>("IncludeDY") : false) {
0009   // Defining standard decay channels
0010   ee.fill(TauDecayMode::Electron);
0011   ee.fill(TauDecayMode::Electron);
0012   mm.fill(TauDecayMode::Muon);
0013   mm.fill(TauDecayMode::Muon);
0014   hh.fill(TauDecayMode::Hadronic);
0015   hh.fill(TauDecayMode::Hadronic);
0016   em.fill(TauDecayMode::Electron);
0017   em.fill(TauDecayMode::Muon);
0018   eh.fill(TauDecayMode::Electron);
0019   eh.fill(TauDecayMode::Hadronic);
0020   mh.fill(TauDecayMode::Muon);
0021   mh.fill(TauDecayMode::Hadronic);
0022
0023   // Filling CutContainers
0024
0025   std::string cut_string_elel = iConfig.getParameter<std::string>("ElElCut");
0026   std::string cut_string_mumu = iConfig.getParameter<std::string>("MuMuCut");
0027   std::string cut_string_hadhad = iConfig.getParameter<std::string>("HadHadCut");
0028   std::string cut_string_elmu = iConfig.getParameter<std::string>("ElMuCut");
0029   std::string cut_string_elhad = iConfig.getParameter<std::string>("ElHadCut");
0030   std::string cut_string_muhad = iConfig.getParameter<std::string>("MuHadCut");
0031
0032   std::vector<std::string> use_final_states = iConfig.getParameter<std::vector<std::string> >("Final_States");
0033
0034   for (std::vector<std::string>::const_iterator final_state = use_final_states.begin();
0035        final_state != use_final_states.end();
0036        ++final_state) {
0037     if ((*final_state) == "ElEl")
0038       fill_cuts(cut_string_elel, ee);
0039     else if ((*final_state) == "MuMu")
0040       fill_cuts(cut_string_mumu, mm);
0041     else if ((*final_state) == "HadHad")
0042       fill_cuts(cut_string_hadhad, hh);
0043     else if ((*final_state) == "ElMu")
0044       fill_cuts(cut_string_elmu, em);
0045     else if ((*final_state) == "ElHad")
0046       fill_cuts(cut_string_elhad, eh);
0047     else if ((*final_state) == "MuHad")
0048       fill_cuts(cut_string_muhad, mh);
0049     else
0050       edm::LogWarning("EmbeddingHepMCFilter")
0051           << (*final_state)
0052           << " this decay channel is not supported. Please choose on of (ElEl,MuMu,HadHad,ElMu,ElHad,MuHad)";
0053   }
0054 }
0055
0056 EmbeddingHepMCFilter::~EmbeddingHepMCFilter() {}
0057
0058 bool EmbeddingHepMCFilter::filter(const HepMC::GenEvent *evt) {
0059   //Reset DecayChannel_ and p4VisPair_ at the beginning of each event.
0060   DecayChannel_.reset();
0061   std::vector<reco::Candidate::LorentzVector> p4VisPair_;
0062
0063   // Going through the particle list. Mother particles are allways before their children.
0064   // One can stop the loop after the second tau is reached and processed.
0065   for (HepMC::GenEvent::particle_const_iterator particle = evt->particles_begin(); particle != evt->particles_end();
0066        ++particle) {
0067     int mom_id = 0;           // no particle available with PDG ID 0
0068     bool isHardProc = false;  // mother is ZPDGID_, or is lepton from hard process (DY process qq -> ll)
0069     int pdg_id = std::abs((*particle)->pdg_id());
0070     HepMC::GenVertex *vertex = (*particle)->production_vertex();
0071     if (vertex != nullptr) {  // search for the mom via the production_vertex
0072       HepMC::GenVertex::particles_in_const_iterator mom = vertex->particles_in_const_begin();
0073       if (mom != vertex->particles_in_const_end()) {
0074         mom_id = std::abs((*mom)->pdg_id());  // mom was found
0075       }
0076       if (mom_id == ZPDGID_) {
0077         isHardProc = true;  // intermediate boson
0078       } else if (includeDY_ && 11 <= pdg_id && pdg_id <= 16 && mcTruthHelper_.isFirstCopy(**particle) &&
0079                  mcTruthHelper_.fromHardProcess(**particle)) {
0080         edm::LogInfo("EmbeddingHepMCFilter") << (*particle)->pdg_id() << " with mother " << (*mom)->pdg_id();
0081         isHardProc = true;  // assume Drell-Yan qq -> ll without intermediate boson
0082       }
0083     }
0084
0085     if (!isHardProc) {
0086       continue;
0087     } else if (pdg_id == tauonPDGID_) {
0088       reco::Candidate::LorentzVector p4Vis;
0089       decay_and_sump4Vis((*particle), p4Vis);  // recursive access to final states.
0090       p4VisPair_.push_back(p4Vis);
0091     } else if (pdg_id == muonPDGID_) {  // Also handle the option when Z-> mumu
0092       reco::Candidate::LorentzVector p4Vis = (reco::Candidate::LorentzVector)(*particle)->momentum();
0093       DecayChannel_.fill(TauDecayMode::Muon);  // take the muon cuts
0094       p4VisPair_.push_back(p4Vis);
0095     } else if (pdg_id == electronPDGID_) {  // Also handle the option when Z-> ee
0096       reco::Candidate::LorentzVector p4Vis = (reco::Candidate::LorentzVector)(*particle)->momentum();
0097       DecayChannel_.fill(TauDecayMode::Electron);  // take the electron cuts
0098       p4VisPair_.push_back(p4Vis);
0099     }
0100   }
0101   // Putting DecayChannel_ in default convention:
0102   // For mixed decay channels use the Electron_Muon, Electron_Hadronic, Muon_Hadronic convention.
0103   // For symmetric decay channels (e.g. Muon_Muon) use Leading_Trailing convention with respect to Pt.
0104   if (p4VisPair_.size() == 2) {
0105     sort_by_convention(p4VisPair_);
0106     edm::LogInfo("EmbeddingHepMCFilter") << "Quantities of the visible decay products:"
0107                                          << "\tPt's: "
0108                                          << " 1st " << p4VisPair_[0].Pt() << ", 2nd " << p4VisPair_[1].Pt()
0109                                          << "\tEta's: "
0110                                          << " 1st " << p4VisPair_[0].Eta() << ", 2nd " << p4VisPair_[1].Eta()
0111                                          << " decay channel: " << return_mode(DecayChannel_.first)
0112                                          << return_mode(DecayChannel_.second);
0113   } else {
0114     edm::LogError("EmbeddingHepMCFilter") << "Size less non equal two";
0115   }
0116
0117   return apply_cuts(p4VisPair_);
0118 }
0119
0120 void EmbeddingHepMCFilter::decay_and_sump4Vis(HepMC::GenParticle *particle, reco::Candidate::LorentzVector &p4Vis) {
0121   bool decaymode_known = false;
0122   for (HepMC::GenVertex::particle_iterator daughter = particle->end_vertex()->particles_begin(HepMC::children);
0123        daughter != particle->end_vertex()->particles_end(HepMC::children);
0124        ++daughter) {
0125     bool neutrino = (std::abs((*daughter)->pdg_id()) == tauon_neutrino_PDGID_) ||
0126                     (std::abs((*daughter)->pdg_id()) == muon_neutrino_PDGID_) ||
0127                     (std::abs((*daughter)->pdg_id()) == electron_neutrino_PDGID_);
0128
0129     // Determining DecayMode, if particle is tau lepton.
0130     // Asuming, that there are only the two tauons from the Z-boson.
0131     // This is the case for the simulated Z->tautau event constructed by EmbeddingLHEProducer.
0132     if (std::abs(particle->pdg_id()) == tauonPDGID_ && !decaymode_known) {
0133       // use these bools to protect againt taus that aren't the last copy (which "decay" to a tau and a gamma)
0134       bool isGamma = std::abs((*daughter)->pdg_id()) == 22;
0135       bool isTau = std::abs((*daughter)->pdg_id()) == 15;
0136       if (std::abs((*daughter)->pdg_id()) == muonPDGID_) {
0137         DecayChannel_.fill(TauDecayMode::Muon);
0138         decaymode_known = true;
0139       } else if (std::abs((*daughter)->pdg_id()) == electronPDGID_) {
0140         DecayChannel_.fill(TauDecayMode::Electron);
0141         decaymode_known = true;
0142       } else if (!neutrino && !isGamma && !isTau) {
0143         DecayChannel_.fill(TauDecayMode::Hadronic);
0144         decaymode_known = true;
0145       }
0146     }
0147     // Adding up all visible momentum in recursive way.
0148     if ((*daughter)->status() == 1 && !neutrino)
0149       p4Vis += (reco::Candidate::LorentzVector)(*daughter)->momentum();
0150     else if (!neutrino)
0151       decay_and_sump4Vis((*daughter), p4Vis);
0152   }
0153 }
0154
0155 void EmbeddingHepMCFilter::sort_by_convention(std::vector<reco::Candidate::LorentzVector> &p4VisPair) {
0156   bool mixed_false_order =
0157       (DecayChannel_.first == TauDecayMode::Hadronic && DecayChannel_.second == TauDecayMode::Muon) ||
0158       (DecayChannel_.first == TauDecayMode::Hadronic && DecayChannel_.second == TauDecayMode::Electron) ||
0159       (DecayChannel_.first == TauDecayMode::Muon && DecayChannel_.second == TauDecayMode::Electron);
0160
0161   if (DecayChannel_.first == DecayChannel_.second && p4VisPair[0].Pt() < p4VisPair[1].Pt()) {
0162     edm::LogVerbatim("EmbeddingHepMCFilter") << "Changing symmetric channels to Leading_Trailing convention in Pt";
0163     edm::LogVerbatim("EmbeddingHepMCFilter") << "Pt's before: " << p4VisPair[0].Pt() << " " << p4VisPair[1].Pt();
0164     std::reverse(p4VisPair.begin(), p4VisPair.end());
0165     edm::LogVerbatim("EmbeddingHepMCFilter") << "Pt's after: " << p4VisPair[0].Pt() << " " << p4VisPair[1].Pt();
0166   } else if (mixed_false_order) {
0167     edm::LogVerbatim("EmbeddingHepMCFilter") << "Swapping order of mixed channels";
0168     edm::LogVerbatim("EmbeddingHepMCFilter") << "Pt's before: " << p4VisPair[0].Pt() << " " << p4VisPair[1].Pt();
0169     DecayChannel_.reverse();
0170     edm::LogVerbatim("EmbeddingHepMCFilter")
0171         << "DecayChannel: " << return_mode(DecayChannel_.first) << return_mode(DecayChannel_.second);
0172     std::reverse(p4VisPair.begin(), p4VisPair.end());
0173     edm::LogVerbatim("EmbeddingHepMCFilter") << "Pt's after: " << p4VisPair[0].Pt() << " " << p4VisPair[1].Pt();
0174   }
0175 }
0176
0177 bool EmbeddingHepMCFilter::apply_cuts(std::vector<reco::Candidate::LorentzVector> &p4VisPair) {
0178   for (std::vector<CutsContainer>::const_iterator cut = cuts_.begin(); cut != cuts_.end(); ++cut) {
0179     if (DecayChannel_.first == cut->decaychannel.first &&
0180         DecayChannel_.second == cut->decaychannel.second) {  // First the match to the decay channel
0181       edm::LogInfo("EmbeddingHepMCFilter")
0182           << "Cut pt1 = " << cut->pt1 << " pt2 = " << cut->pt2 << " abs(eta1) = " << cut->eta1
0183           << " abs(eta2) = " << cut->eta2 << " decay channel: " << return_mode(cut->decaychannel.first)
0184           << return_mode(cut->decaychannel.second);
0185
0186       if ((cut->pt1 == -1. || (p4VisPair[0].Pt() > cut->pt1)) && (cut->pt2 == -1. || (p4VisPair[1].Pt() > cut->pt2)) &&
0187           (cut->eta1 == -1. || (std::abs(p4VisPair[0].Eta()) < cut->eta1)) &&
0188           (cut->eta2 == -1. || (std::abs(p4VisPair[1].Eta()) < cut->eta2))) {
0189         edm::LogInfo("EmbeddingHepMCFilter") << "This cut was passed (Stop here and take the event)";
0190         return true;
0191       }
0192     }
0193   }
0194   return false;
0195 }
0196
0197 void EmbeddingHepMCFilter::fill_cuts(std::string cut_string, EmbeddingHepMCFilter::DecayChannel &dc) {
0198   edm::LogInfo("EmbeddingHepMCFilter") << return_mode(dc.first) << return_mode(dc.second) << "Cut : " << cut_string;
0199   boost::trim_fill(cut_string, "");
0200   std::vector<std::string> cut_paths;
0201   boost::split(cut_paths, cut_string, boost::is_any_of("||"), boost::token_compress_on);
0202   for (unsigned int i = 0; i < cut_paths.size(); ++i) {
0203     // Translating the cuts of a path into a struct which is later accessed to apply them on a event.
0204     CutsContainer cut;
0205     fill_cut(cut_paths[i], dc, cut);
0206     cuts_.push_back(cut);
0207   }
0208 }
0209
0210 void EmbeddingHepMCFilter::fill_cut(std::string cut_string,
0211                                     EmbeddingHepMCFilter::DecayChannel &dc,
0212                                     CutsContainer &cut) {
0213   cut.decaychannel = dc;
0214
0215   boost::replace_all(cut_string, "(", "");
0216   boost::replace_all(cut_string, ")", "");
0217   std::vector<std::string> single_cuts;
0218   boost::split(single_cuts, cut_string, boost::is_any_of("&&"), boost::token_compress_on);
0219   for (unsigned int i = 0; i < single_cuts.size(); ++i) {
0220     std::string pt1_str, pt2_str, eta1_str, eta2_str;
0221     if (dc.first == dc.second) {
0222       pt1_str = return_mode(dc.first) + "1" + ".Pt" + ">";
0223       pt2_str = return_mode(dc.second) + "2" + ".Pt" + ">";
0224       eta1_str = return_mode(dc.first) + "1" + ".Eta" + "<";
0225       eta2_str = return_mode(dc.second) + "2" + ".Eta" + "<";
0226     } else {
0227       pt1_str = return_mode(dc.first) + ".Pt" + ">";
0228       pt2_str = return_mode(dc.second) + ".Pt" + ">";
0229       eta1_str = return_mode(dc.first) + ".Eta" + "<";
0230       eta2_str = return_mode(dc.second) + ".Eta" + "<";
0231     }
0232
0233     if (boost::find_first(single_cuts[i], pt1_str)) {
0234       boost::erase_first(single_cuts[i], pt1_str);
0235       cut.pt1 = std::stod(single_cuts[i]);
0236     } else if (boost::find_first(single_cuts[i], pt2_str)) {
0237       boost::erase_first(single_cuts[i], pt2_str);
0238       cut.pt2 = std::stod(single_cuts[i]);
0239     } else if (boost::find_first(single_cuts[i], eta1_str)) {
0240       boost::erase_first(single_cuts[i], eta1_str);
0241       cut.eta1 = std::stod(single_cuts[i]);
0242     } else if (boost::find_first(single_cuts[i], eta2_str)) {
0243       boost::erase_first(single_cuts[i], eta2_str);
0244       cut.eta2 = std::stod(single_cuts[i]);
0245     }
0246   }
0247 }