Tau/src/HLTTauMCProducer.cc

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 numOtherParticles = 0;
0122
0123       for (GenParticleRefVector::const_iterator pit = decayProducts.begin(); pit != decayProducts.end(); ++pit) {
0124         int pdg_id = abs((*pit)->pdgId());
0125         if (pdg_id == 11)
0126           numElectrons++;
0127         else if (pdg_id == 13)
0128           numMuons++;
0129         else if (pdg_id == 211 || pdg_id == 321)
0130           numChargedPions++;  // Count both pi+ and K+
0131         else if (pdg_id == 111 || pdg_id == 130 || pdg_id == 310)
0132           numNeutralPions++;  // Count both pi0 and K0_L/S
0133         else if (pdg_id == 12 || pdg_id == 14 || pdg_id == 16) {
0134           if (pdg_id == 16) {
0135             Neutrino.SetPxPyPzE((*pit)->px(), (*pit)->py(), (*pit)->pz(), (*pit)->energy());
0136           }
0137         } else if (pdg_id == 22)
0138           numPhotons++;
0139         else {
0140           numOtherParticles++;
0141         }
0142
0143         if (pdg_id != 12 && pdg_id != 14 && pdg_id != 16) {
0144           TauDecayProduct.SetPxPyPzE((*pit)->px(), (*pit)->py(), (*pit)->pz(), (*pit)->energy());
0145           Visible_Taus += TauDecayProduct;
0146         }
0147         //        cout<< "This has to be the same: " << (*pit)->pdgId()
0148         //<< " "<< (*pit)->status()<< " mother: "<< (*pit)->mother()->pdgId() <<
0149         // endl;
0150       }
0151
0152       int tauDecayMode = kOther;
0153
0154       if (numOtherParticles == 0) {
0155         if (numElectrons == 1) {
0156           //--- tau decays into electrons
0157           tauDecayMode = kElectron;
0158         } else if (numMuons == 1) {
0159           //--- tau decays into muons
0160           tauDecayMode = kMuon;
0161         } else {
0162           //--- hadronic tau decays
0163           switch (numChargedPions) {
0164             case 1:
0165               if (numNeutralPions != 0) {
0166                 tauDecayMode = kOther;
0167                 break;
0168               }
0169               switch (numPhotons) {
0170                 case 0:
0171                   tauDecayMode = kOneProng0pi0;
0172                   break;
0173                 case 2:
0174                   tauDecayMode = kOneProng1pi0;
0175                   break;
0176                 case 4:
0177                   tauDecayMode = kOneProng2pi0;
0178                   break;
0179                 default:
0180                   tauDecayMode = kOther;
0181                   break;
0182               }
0183               break;
0184             case 3:
0185               if (numNeutralPions != 0) {
0186                 tauDecayMode = kOther;
0187                 break;
0188               }
0189               switch (numPhotons) {
0190                 case 0:
0191                   tauDecayMode = kThreeProng0pi0;
0192                   break;
0193                 case 2:
0194                   tauDecayMode = kThreeProng1pi0;
0195                   break;
0196                 default:
0197                   tauDecayMode = kOther;
0198                   break;
0199               }
0200               break;
0201           }
0202         }
0203       }
0204
0205       //          cout<< "So we have a: " << tauDecayMode <<endl;
0206       if (tauDecayMode == kElectron) {
0207         if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
0208              Visible_Taus.phi() < phiMax_) &&
0209             (Visible_Taus.pt() > ptMinMCElectron_)) {
0210           product_Electrons->push_back(Visible_Taus);
0211           product_Leptons->push_back(Visible_Taus);
0212         }
0213       } else if (tauDecayMode == kMuon) {
0214         if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
0215              Visible_Taus.phi() < phiMax_) &&
0216             (Visible_Taus.pt() > ptMinMCMuon_)) {
0217           product_Muons->push_back(Visible_Taus);
0218           product_Leptons->push_back(Visible_Taus);
0219         }
0220       } else if (tauDecayMode == kOneProng0pi0 || tauDecayMode == kOneProng1pi0 || tauDecayMode == kOneProng2pi0) {
0221         if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
0222              Visible_Taus.phi() < phiMax_) &&
0223             (Visible_Taus.pt() > ptMinMCTau_)) {
0224           product_OneProng->push_back(Visible_Taus);
0225           product_OneAndThreeProng->push_back(Visible_Taus);
0226           product_Neutrina->push_back(Neutrino);
0227         }
0228       } else if (tauDecayMode == kThreeProng0pi0 || tauDecayMode == kThreeProng1pi0) {
0229         if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
0230              Visible_Taus.phi() < phiMax_) &&
0231             (Visible_Taus.pt() > ptMinMCTau_)) {
0232           product_ThreeProng->push_back(Visible_Taus);
0233           product_OneAndThreeProng->push_back(Visible_Taus);
0234           product_Neutrina->push_back(Neutrino);
0235         }
0236       } else if (tauDecayMode == kOther) {
0237         if ((Visible_Taus.eta() > etaMin_ && Visible_Taus.eta() < etaMax_ && Visible_Taus.phi() > phiMin_ &&
0238              Visible_Taus.phi() < phiMax_) &&
0239             (Visible_Taus.pt() > ptMinMCTau_)) {
0240           product_Other->push_back(Visible_Taus);
0241         }
0242       }
0243     }
0244   }
0245
0246   iEvent.put(std::move(product_Leptons), "LeptonicTauLeptons");
0247   iEvent.put(std::move(product_Electrons), "LeptonicTauElectrons");
0248   iEvent.put(std::move(product_Muons), "LeptonicTauMuons");
0249   iEvent.put(std::move(product_OneProng), "HadronicTauOneProng");
0250   iEvent.put(std::move(product_ThreeProng), "HadronicTauThreeProng");
0251   iEvent.put(std::move(product_OneAndThreeProng), "HadronicTauOneAndThreeProng");
0252   iEvent.put(std::move(product_Other), "TauOther");
0253   iEvent.put(std::move(product_Neutrina), "Neutrina");
0254   iEvent.put(std::move(product_MET), "MET");
0255   iEvent.put(std::move(product_Mothers), "Mothers");
0256 }
0257
0258 // Helper Function
0259
0260 void HLTTauMCProducer::getGenDecayProducts(const GenParticleRef &mother,
0261                                            GenParticleRefVector &products,
0262                                            int status,
0263                                            int pdgId) const {
0264   const GenParticleRefVector &daughterRefs = mother->daughterRefVector();
0265
0266   for (GenParticleRefVector::const_iterator d = daughterRefs.begin(); d != daughterRefs.end(); ++d) {
0267     if ((status == 0 || (*d)->status() == status) && (pdgId == 0 || std::abs((*d)->pdgId()) == pdgId)) {
0268       products.push_back(*d);
0269     } else
0270       getGenDecayProducts(*d, products, status, pdgId);
0271   }
0272 }