HepMCCandAlgos/plugins/GenVisTauProducer.cc

0001 #include "FWCore/Framework/interface/global/EDProducer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0006 #include "DataFormats/Common/interface/View.h"
0007 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0008 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0009 #include "DataFormats/JetReco/interface/GenJet.h"
0010 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0011 #include "PhysicsTools/JetMCUtils/interface/JetMCTag.h"
0012 #include "DataFormats/TauReco/interface/PFTau.h"
0013 #include "DataFormats/Math/interface/deltaR.h"
0014 #include "CommonTools/Utils/interface/PtComparator.h"
0015
0016 #include <vector>
0017 #include <iostream>
0018
0019 class GenVisTauProducer : public edm::global::EDProducer<> {
0020 public:
0021   GenVisTauProducer(const edm::ParameterSet& params)
0022       : src_(consumes<reco::GenJetCollection>(params.getParameter<edm::InputTag>("src"))),
0023         srcGenParticles_(consumes<reco::GenParticleCollection>(params.getParameter<edm::InputTag>("srcGenParticles"))),
0024         pTComparator_() {
0025     produces<reco::GenParticleCollection>();
0026   }
0027
0028   ~GenVisTauProducer() override {}
0029
0030   void produce(edm::StreamID id, edm::Event& evt, const edm::EventSetup& es) const override {
0031     edm::Handle<reco::GenJetCollection> genTauJets;
0032     evt.getByToken(src_, genTauJets);
0033
0034     edm::Handle<reco::GenParticleCollection> genParticles;
0035     evt.getByToken(srcGenParticles_, genParticles);
0036     size_t numGenParticles = genParticles->size();
0037
0038     auto genVisTaus = std::make_unique<reco::GenParticleCollection>();
0039
0040     for (const auto& genTauJet : *genTauJets) {
0041       std::string decayMode_string = JetMCTagUtils::genTauDecayMode(genTauJet);
0042       // CV: store hadronic tau decays only
0043       if (decayMode_string == "electron" || decayMode_string == "muon")
0044         continue;
0045       int decayMode = reco::PFTau::kNull;
0046       if (decayMode_string == "oneProng0Pi0")
0047         decayMode = reco::PFTau::kOneProng0PiZero;
0048       else if (decayMode_string == "oneProng1Pi0")
0049         decayMode = reco::PFTau::kOneProng1PiZero;
0050       else if (decayMode_string == "oneProng2Pi0")
0051         decayMode = reco::PFTau::kOneProng2PiZero;
0052       else if (decayMode_string == "threeProng0Pi0")
0053         decayMode = reco::PFTau::kThreeProng0PiZero;
0054       else if (decayMode_string == "threeProng1Pi0")
0055         decayMode = reco::PFTau::kThreeProng1PiZero;
0056       else
0057         decayMode = reco::PFTau::kRareDecayMode;
0058
0059       int pdgId = (genTauJet.charge() > 0) ? -15 : +15;
0060
0061       // CV: store decayMode in status flag of GenParticle object
0062       reco::GenParticle genVisTau(genTauJet.charge(), genTauJet.p4(), genTauJet.vertex(), pdgId, decayMode, true);
0063
0064       // CV: find tau lepton "mother" particle
0065       for (size_t idxGenParticle = 0; idxGenParticle < numGenParticles; ++idxGenParticle) {
0066         const reco::GenParticle& genTau = (*genParticles)[idxGenParticle];
0067         if (abs(genTau.pdgId()) == 15 && genTau.status() == 2) {
0068           reco::Candidate::LorentzVector daughterVisP4;
0069           for (const reco::GenParticleRef& daughter : genTau.daughterRefVector()) {
0070             int abs_pdgId = abs(daughter->pdgId());
0071             // CV: skip neutrinos
0072             if (abs_pdgId == 12 || abs_pdgId == 14 || abs_pdgId == 16)
0073               continue;
0074             daughterVisP4 += daughter->p4();
0075           }
0076           double dR2 = deltaR2(daughterVisP4, genVisTau);
0077           if (dR2 < 1.e-4) {
0078             genVisTau.addMother(reco::GenParticleRef(genParticles, idxGenParticle));
0079             break;
0080           }
0081         }
0082       }
0083
0084       genVisTaus->push_back(genVisTau);
0085     }
0086
0087     std::sort(genVisTaus->begin(), genVisTaus->end(), pTComparator_);
0088     evt.put(std::move(genVisTaus));
0089   }
0090
0091   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0092     edm::ParameterSetDescription desc;
0093     desc.add<edm::InputTag>("src")->setComment("collection of visible gen taus (as reco::GenJetCollection)");
0094     desc.add<edm::InputTag>("srcGenParticles")->setComment("collections of gen particles");
0095     descriptions.add("genVisTaus", desc);
0096   }
0097
0098 private:
0099   const edm::EDGetTokenT<reco::GenJetCollection> src_;
0100   const edm::EDGetTokenT<reco::GenParticleCollection> srcGenParticles_;
0101   const GreaterByPt<reco::GenParticle> pTComparator_;
0102 };
0103
0104 #include "FWCore/Framework/interface/MakerMacros.h"
0105 DEFINE_FWK_MODULE(GenVisTauProducer);