JetMCAlgos/plugins/TauGenJetProducer.cc

0001 #include "PhysicsTools/JetMCAlgos/plugins/TauGenJetProducer.h"
0002
0003 #include "FWCore/Framework/interface/ESHandle.h"
0004
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/Exception.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008
0009 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0010 #include "DataFormats/JetReco/interface/GenJet.h"
0011 #include "DataFormats/Common/interface/RefToPtr.h"
0012
0013 #include "PhysicsTools/HepMCCandAlgos/interface/GenParticlesHelper.h"
0014 #include "PhysicsTools/JetMCUtils/interface/JetMCTag.h"
0015 #include "DataFormats/TauReco/interface/PFTau.h"
0016
0017 using namespace std;
0018 using namespace edm;
0019 using namespace reco;
0020
0021 namespace {
0022   //Map to convert names of decay modes to integer codes
0023   const std::map<std::string, int> decayModeStringToCodeMap = {{"null", PFTau::kNull},
0024                                                                {"oneProng0Pi0", PFTau::kOneProng0PiZero},
0025                                                                {"oneProng1Pi0", PFTau::kOneProng1PiZero},
0026                                                                {"oneProng2Pi0", PFTau::kOneProng2PiZero},
0027                                                                {"threeProng0Pi0", PFTau::kThreeProng0PiZero},
0028                                                                {"threeProng1Pi0", PFTau::kThreeProng1PiZero},
0029                                                                {"electron", PFTau::kRareDecayMode + 1},
0030                                                                {"muon", PFTau::kRareDecayMode + 2},
0031                                                                {"rare", PFTau::kRareDecayMode},
0032                                                                {"tau", PFTau::kNull - 1}};
0033 }  // namespace
0034
0035 TauGenJetProducer::TauGenJetProducer(const edm::ParameterSet& iConfig)
0036     : tokenGenParticles_(consumes<GenParticleCollection>(iConfig.getParameter<InputTag>("GenParticles"))),
0037       includeNeutrinos_(iConfig.getParameter<bool>("includeNeutrinos")),
0038       verbose_(iConfig.getUntrackedParameter<bool>("verbose", false)) {
0039   produces<GenJetCollection>();
0040 }
0041
0042 TauGenJetProducer::~TauGenJetProducer() {}
0043
0044 void TauGenJetProducer::produce(edm::StreamID, Event& iEvent, const EventSetup& iSetup) const {
0045   Handle<GenParticleCollection> genParticles = iEvent.getHandle(tokenGenParticles_);
0046
0047   auto pOutVisTaus = std::make_unique<GenJetCollection>();
0048
0049   using namespace GenParticlesHelper;
0050
0051   GenParticleRefVector allStatus2Taus;
0052   findParticles(*genParticles, allStatus2Taus, 15, 2);
0053
0054   for (IGR iTau = allStatus2Taus.begin(); iTau != allStatus2Taus.end(); ++iTau) {
0055     // look for all status 1 (stable) descendents
0056     GenParticleRefVector descendents;
0057     findDescendents(*iTau, descendents, 1);
0058     if (descendents.empty()) {
0059       edm::LogWarning("NoTauDaughters") << "Tau p4: " << (*iTau)->p4() << " vtx: " << (*iTau)->vertex()
0060                                         << " has no daughters";
0061
0062       math::XYZPoint vertex;
0063       GenJet::Specific specific;
0064       Jet::Constituents constituents;
0065
0066       constituents.push_back(refToPtr(*iTau));
0067       GenJet jet((*iTau)->p4(), vertex, specific, constituents);
0068       jet.setCharge((*iTau)->charge());
0069       jet.setStatus(decayModeStringToCodeMap.at("tau"));
0070       pOutVisTaus->emplace_back(std::move(jet));
0071       continue;
0072     }
0073     // CV: skip status 2 taus that radiate-off a photon
0074     //    --> have a status 2 tau lepton in the list of descendents
0075     GenParticleRefVector status2TauDaughters;
0076     findDescendents(*iTau, status2TauDaughters, 2, 15);
0077     if (!status2TauDaughters.empty())
0078       continue;
0079
0080     // loop on descendents, and take all except neutrinos
0081     math::XYZTLorentzVector sumVisMom;
0082     Particle::Charge charge = 0;
0083     Jet::Constituents constituents;
0084
0085     if (verbose_)
0086       cout << "tau " << (*iTau) << endl;
0087
0088     for (IGR igr = descendents.begin(); igr != descendents.end(); ++igr) {
0089       int absPdgId = abs((*igr)->pdgId());
0090
0091       // neutrinos
0092       if (!includeNeutrinos_) {
0093         if (absPdgId == 12 || absPdgId == 14 || absPdgId == 16)
0094           continue;
0095       }
0096
0097       if (verbose_)
0098         cout << "\t" << (*igr) << endl;
0099
0100       charge += (*igr)->charge();
0101       sumVisMom += (*igr)->p4();
0102
0103       // need to convert the vector of reference to the constituents
0104       // to a vector of pointers to build the genjet
0105       constituents.push_back(refToPtr(*igr));
0106     }
0107
0108     math::XYZPoint vertex;
0109     GenJet::Specific specific;
0110
0111     GenJet jet(sumVisMom, vertex, specific, constituents);
0112
0113     if (charge != (*iTau)->charge())
0114       edm::LogError("TauGenJetProducer") << " charge of Tau: " << (*iTau)
0115                                          << " not equal to charge of sum of charge of all descendents.\n"
0116                                          << " Tau's charge: " << (*iTau)->charge() << " sum: " << charge
0117                                          << " # descendents: " << constituents.size() << "\n";
0118
0119     jet.setCharge(charge);
0120     // determine tau decay mode and set it as jet status
0121     if (auto search = decayModeStringToCodeMap.find(JetMCTagUtils::genTauDecayMode(jet));
0122         search != decayModeStringToCodeMap.end())
0123       jet.setStatus(search->second);
0124     else
0125       jet.setStatus(decayModeStringToCodeMap.at("null"));
0126     pOutVisTaus->push_back(jet);
0127   }
0128   iEvent.put(std::move(pOutVisTaus));
0129 }
0130
0131 DEFINE_FWK_MODULE(TauGenJetProducer);