Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:29:38

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