File indexing completed on 2023-05-10 03:54:28
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
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
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
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
0062
0063
0064
0065 for (GenParticleCollection::const_iterator p = genParticles->begin(); p != genParticles->end(); ++p) {
0066
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
0072
0073 break;
0074 }
0075 }
0076
0077
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 }
0084
0085
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
0091 if (!genP.isPromptDecayed())
0092 continue;
0093
0094 if (std::abs(genP.pdgId()) == 15) {
0095 GenParticleRef genRef(genParticles, index);
0096
0097 GenParticleRefVector daugTaus;
0098 getGenDecayProducts(genRef, daugTaus, 0, 15);
0099 if (daugTaus.empty())
0100 allTaus.push_back(genRef);
0101 }
0102 }
0103
0104
0105 for (GenParticleRefVector::const_iterator t = allTaus.begin(); t != allTaus.end(); ++t) {
0106
0107 GenParticleRefVector decayProducts;
0108 getGenDecayProducts(*t, decayProducts, 1);
0109
0110
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++;
0131 else if (pdg_id == 111 || pdg_id == 130 || pdg_id == 310)
0132 numNeutralPions++;
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
0148
0149
0150 }
0151
0152 int tauDecayMode = kOther;
0153
0154 if (numOtherParticles == 0) {
0155 if (numElectrons == 1) {
0156
0157 tauDecayMode = kElectron;
0158 } else if (numMuons == 1) {
0159
0160 tauDecayMode = kMuon;
0161 } else {
0162
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
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
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 }