Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:21:05

0001 #include "HLTriggerOffline/Tau/interface/HLTTauMCProducer.h"
0002 
0003 using namespace edm;
0004 using namespace std;
0005 using namespace reco;
0006 
0007 HLTTauMCProducer::HLTTauMCProducer(const edm::ParameterSet &mc)
0008     : MC_{consumes<GenParticleCollection>(mc.getUntrackedParameter<edm::InputTag>("GenParticles"))},
0009       MCMET_{consumes<GenMETCollection>(mc.getUntrackedParameter<edm::InputTag>("GenMET"))},
0010       ptMinMCTau_{mc.getUntrackedParameter<double>("ptMinTau", 5.)},
0011       ptMinMCElectron_{mc.getUntrackedParameter<double>("ptMinElectron", 5.)},
0012       ptMinMCMuon_{mc.getUntrackedParameter<double>("ptMinMuon", 2.)},
0013       m_PDG_{mc.getUntrackedParameter<std::vector<int>>("BosonID")},
0014       etaMin_{mc.getUntrackedParameter<double>("EtaMin", -2.5)},
0015       etaMax_{mc.getUntrackedParameter<double>("EtaMax", 2.5)},
0016       phiMin_{mc.getUntrackedParameter<double>("PhiMin", -3.15)},
0017       phiMax_{mc.getUntrackedParameter<double>("PhiMax", 3.15)} {
0018   // One Parameter Set per Collection
0019 
0020   produces<LorentzVectorCollection>("LeptonicTauLeptons");
0021   produces<LorentzVectorCollection>("LeptonicTauElectrons");
0022   produces<LorentzVectorCollection>("LeptonicTauMuons");
0023   produces<LorentzVectorCollection>("HadronicTauOneProng");
0024   produces<LorentzVectorCollection>("HadronicTauThreeProng");
0025   produces<LorentzVectorCollection>("HadronicTauOneAndThreeProng");
0026   produces<LorentzVectorCollection>("TauOther");
0027   produces<LorentzVectorCollection>("Neutrina");
0028   produces<LorentzVectorCollection>("MET");
0029   produces<std::vector<int>>("Mothers");
0030 }
0031 
0032 void HLTTauMCProducer::produce(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iES) const {
0033   // All the code from HLTTauMCInfo is here :-)
0034 
0035   unique_ptr<LorentzVectorCollection> product_Electrons(new LorentzVectorCollection);
0036   unique_ptr<LorentzVectorCollection> product_Muons(new LorentzVectorCollection);
0037   unique_ptr<LorentzVectorCollection> product_Leptons(new LorentzVectorCollection);
0038   unique_ptr<LorentzVectorCollection> product_OneProng(new LorentzVectorCollection);
0039   unique_ptr<LorentzVectorCollection> product_ThreeProng(new LorentzVectorCollection);
0040   unique_ptr<LorentzVectorCollection> product_OneAndThreeProng(new LorentzVectorCollection);
0041   unique_ptr<LorentzVectorCollection> product_Other(new LorentzVectorCollection);
0042   unique_ptr<LorentzVectorCollection> product_Neutrina(new LorentzVectorCollection);
0043   unique_ptr<LorentzVectorCollection> product_MET(new LorentzVectorCollection);
0044   unique_ptr<std::vector<int>> product_Mothers(new std::vector<int>);
0045 
0046   edm::Handle<GenParticleCollection> genParticles;
0047   iEvent.getByToken(MC_, genParticles);
0048 
0049   if (!genParticles.isValid())
0050     return;
0051 
0052   // Look for MET
0053   edm::Handle<reco::GenMETCollection> genMet;
0054   iEvent.getByToken(MCMET_, genMet);
0055   LorentzVector MET(0., 0., 0., 0.);
0056   if (genMet.isValid()) {
0057     MET = LorentzVector(genMet->front().px(), genMet->front().py(), 0, genMet->front().pt());
0058   }
0059   product_MET->push_back(MET);
0060 
0061   // Look for primary bosons
0062   // It is not guaranteed that primary bosons are stored in event history.
0063   // Is it really needed when check if taus from the boson is removed?
0064   // Kept for backward compatibility
0065   for (GenParticleCollection::const_iterator p = genParticles->begin(); p != genParticles->end(); ++p) {
0066     // Check the PDG ID
0067     bool pdg_ok = false;
0068     for (size_t pi = 0; pi < m_PDG_.size(); ++pi) {
0069       if (abs((*p).pdgId()) == m_PDG_[pi] && ((*p).isHardProcess() || (*p).status() == 3)) {
0070         pdg_ok = true;
0071         // cout<<" Bsoson particles: "<< (*p).pdgId()<< " " <<(*p).status() << "
0072         // "<< pdg_ok<<endl;
0073         break;
0074       }
0075     }
0076 
0077     // Check if the boson is one of interest and if there is a valid vertex
0078     if (pdg_ok) {
0079       product_Mothers->push_back((*p).pdgId());
0080 
0081       TLorentzVector Boson((*p).px(), (*p).py(), (*p).pz(), (*p).energy());
0082     }
0083   }  // End of search for the bosons
0084 
0085   // Look for taus
0086   GenParticleRefVector allTaus;
0087   unsigned index = 0;
0088   for (GenParticleCollection::const_iterator p = genParticles->begin(); p != genParticles->end(); ++p, ++index) {
0089     const GenParticle &genP = *p;
0090     // accept only isPromptDecayed() particles
0091     if (!genP.isPromptDecayed())
0092       continue;
0093     // check if it is tau, i.e. if |pdgId|=15
0094     if (std::abs(genP.pdgId()) == 15) {
0095       GenParticleRef genRef(genParticles, index);
0096       // check if it is the last tau in decay/radiation chain
0097       GenParticleRefVector daugTaus;
0098       getGenDecayProducts(genRef, daugTaus, 0, 15);
0099       if (daugTaus.empty())
0100         allTaus.push_back(genRef);
0101     }
0102   }
0103 
0104   // Find stable tau decay products and build visible taus
0105   for (GenParticleRefVector::const_iterator t = allTaus.begin(); t != allTaus.end(); ++t) {
0106     // look for all stable (status=1) decay products
0107     GenParticleRefVector decayProducts;
0108     getGenDecayProducts(*t, decayProducts, 1);
0109 
0110     // build visible taus and recognize decay mode
0111     if (!decayProducts.empty()) {
0112       LorentzVector Visible_Taus(0., 0., 0., 0.);
0113       LorentzVector TauDecayProduct(0., 0., 0., 0.);
0114       LorentzVector Neutrino(0., 0., 0., 0.);
0115 
0116       int numElectrons = 0;
0117       int numMuons = 0;
0118       int numChargedPions = 0;
0119       int numNeutralPions = 0;
0120       int numPhotons = 0;
0121       int numNeutrinos = 0;
0122       int numOtherParticles = 0;
0123 
0124       for (GenParticleRefVector::const_iterator pit = decayProducts.begin(); pit != decayProducts.end(); ++pit) {
0125         int pdg_id = abs((*pit)->pdgId());
0126         if (pdg_id == 11)
0127           numElectrons++;
0128         else if (pdg_id == 13)
0129           numMuons++;
0130         else if (pdg_id == 211 || pdg_id == 321)
0131           numChargedPions++;  // Count both pi+ and K+
0132         else if (pdg_id == 111 || pdg_id == 130 || pdg_id == 310)
0133           numNeutralPions++;  // Count both pi0 and K0_L/S
0134         else if (pdg_id == 12 || pdg_id == 14 || pdg_id == 16) {
0135           numNeutrinos++;
0136           if (pdg_id == 16) {
0137             Neutrino.SetPxPyPzE((*pit)->px(), (*pit)->py(), (*pit)->pz(), (*pit)->energy());
0138           }
0139         } else if (pdg_id == 22)
0140           numPhotons++;
0141         else {
0142           numOtherParticles++;
0143         }
0144 
0145         if (pdg_id != 12 && pdg_id != 14 && pdg_id != 16) {
0146           TauDecayProduct.SetPxPyPzE((*pit)->px(), (*pit)->py(), (*pit)->pz(), (*pit)->energy());
0147           Visible_Taus += TauDecayProduct;
0148         }
0149         //        cout<< "This has to be the same: " << (*pit)->pdgId()
0150         //<< " "<< (*pit)->status()<< " mother: "<< (*pit)->mother()->pdgId() <<
0151         // endl;
0152       }
0153 
0154       int tauDecayMode = kOther;
0155 
0156       if (numOtherParticles == 0) {
0157         if (numElectrons == 1) {
0158           //--- tau decays into electrons
0159           tauDecayMode = kElectron;
0160         } else if (numMuons == 1) {
0161           //--- tau decays into muons
0162           tauDecayMode = kMuon;
0163         } else {
0164           //--- hadronic tau decays
0165           switch (numChargedPions) {
0166             case 1:
0167               if (numNeutralPions != 0) {
0168                 tauDecayMode = kOther;
0169                 break;
0170               }
0171               switch (numPhotons) {
0172                 case 0:
0173                   tauDecayMode = kOneProng0pi0;
0174                   break;
0175                 case 2:
0176                   tauDecayMode = kOneProng1pi0;
0177                   break;
0178                 case 4:
0179                   tauDecayMode = kOneProng2pi0;
0180                   break;
0181                 default:
0182                   tauDecayMode = kOther;
0183                   break;
0184               }
0185               break;
0186             case 3:
0187               if (numNeutralPions != 0) {
0188                 tauDecayMode = kOther;
0189                 break;
0190               }
0191               switch (numPhotons) {
0192                 case 0:
0193                   tauDecayMode = kThreeProng0pi0;
0194                   break;
0195                 case 2:
0196                   tauDecayMode = kThreeProng1pi0;
0197                   break;
0198                 default:
0199                   tauDecayMode = kOther;
0200                   break;
0201               }
0202               break;
0203           }
0204         }
0205       }
0206 
0207       //          cout<< "So we have a: " << tauDecayMode <<endl;
0208       if (tauDecayMode == kElectron) {
0209         if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
0210              Visible_Taus.phi() < phiMax_) &&
0211             (Visible_Taus.pt() > ptMinMCElectron_)) {
0212           product_Electrons->push_back(Visible_Taus);
0213           product_Leptons->push_back(Visible_Taus);
0214         }
0215       } else if (tauDecayMode == kMuon) {
0216         if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
0217              Visible_Taus.phi() < phiMax_) &&
0218             (Visible_Taus.pt() > ptMinMCMuon_)) {
0219           product_Muons->push_back(Visible_Taus);
0220           product_Leptons->push_back(Visible_Taus);
0221         }
0222       } else if (tauDecayMode == kOneProng0pi0 || tauDecayMode == kOneProng1pi0 || tauDecayMode == kOneProng2pi0) {
0223         if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
0224              Visible_Taus.phi() < phiMax_) &&
0225             (Visible_Taus.pt() > ptMinMCTau_)) {
0226           product_OneProng->push_back(Visible_Taus);
0227           product_OneAndThreeProng->push_back(Visible_Taus);
0228           product_Neutrina->push_back(Neutrino);
0229         }
0230       } else if (tauDecayMode == kThreeProng0pi0 || tauDecayMode == kThreeProng1pi0) {
0231         if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
0232              Visible_Taus.phi() < phiMax_) &&
0233             (Visible_Taus.pt() > ptMinMCTau_)) {
0234           product_ThreeProng->push_back(Visible_Taus);
0235           product_OneAndThreeProng->push_back(Visible_Taus);
0236           product_Neutrina->push_back(Neutrino);
0237         }
0238       } else if (tauDecayMode == kOther) {
0239         if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
0240              Visible_Taus.phi() < phiMax_) &&
0241             (Visible_Taus.pt() > ptMinMCTau_)) {
0242           product_Other->push_back(Visible_Taus);
0243         }
0244       }
0245     }
0246   }
0247 
0248   iEvent.put(std::move(product_Leptons), "LeptonicTauLeptons");
0249   iEvent.put(std::move(product_Electrons), "LeptonicTauElectrons");
0250   iEvent.put(std::move(product_Muons), "LeptonicTauMuons");
0251   iEvent.put(std::move(product_OneProng), "HadronicTauOneProng");
0252   iEvent.put(std::move(product_ThreeProng), "HadronicTauThreeProng");
0253   iEvent.put(std::move(product_OneAndThreeProng), "HadronicTauOneAndThreeProng");
0254   iEvent.put(std::move(product_Other), "TauOther");
0255   iEvent.put(std::move(product_Neutrina), "Neutrina");
0256   iEvent.put(std::move(product_MET), "MET");
0257   iEvent.put(std::move(product_Mothers), "Mothers");
0258 }
0259 
0260 // Helper Function
0261 
0262 void HLTTauMCProducer::getGenDecayProducts(const GenParticleRef &mother,
0263                                            GenParticleRefVector &products,
0264                                            int status,
0265                                            int pdgId) const {
0266   const GenParticleRefVector &daughterRefs = mother->daughterRefVector();
0267 
0268   for (GenParticleRefVector::const_iterator d = daughterRefs.begin(); d != daughterRefs.end(); ++d) {
0269     if ((status == 0 || (*d)->status() == status) && (pdgId == 0 || std::abs((*d)->pdgId()) == pdgId)) {
0270       products.push_back(*d);
0271     } else
0272       getGenDecayProducts(*d, products, status, pdgId);
0273   }
0274 }