Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-09 02:22:30 

0001 // system include files
0002 #include <memory>
0003 
0004 // user include files
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/stream/EDProducer.h"
0007 
0008 #include "FWCore/Framework/interface/Event.h"
0009 
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0013 
0014 #include "DataFormats/PatCandidates/interface/Jet.h"
0015 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0016 #include "DataFormats/Math/interface/deltaR.h"
0017 
0018 #include "DataFormats/NanoAOD/interface/FlatTable.h"
0019 
0020 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0021 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0022 
0023 class BTVMCFlavourTableProducer : public edm::stream::EDProducer<> {
0024 public:
0025   explicit BTVMCFlavourTableProducer(const edm::ParameterSet& iConfig)
0026       : name_(iConfig.getParameter<std::string>("name")),
0027         src_(consumes<std::vector<pat::Jet> >(iConfig.getParameter<edm::InputTag>("src"))),
0028         genParticlesToken_(consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genparticles"))) {
0029     produces<nanoaod::FlatTable>();
0030   }
0031 
0032   ~BTVMCFlavourTableProducer() override{};
0033   int jet_flavour(const pat::Jet& jet,
0034                   const std::vector<reco::GenParticle>& gToBB,
0035                   const std::vector<reco::GenParticle>& gToCC,
0036                   const std::vector<reco::GenParticle>& neutrinosLepB,
0037                   const std::vector<reco::GenParticle>& neutrinosLepB_C,
0038                   const std::vector<reco::GenParticle>& alltaus,
0039                   bool usePhysForLightAndUndefined) {
0040     int hflav = abs(jet.hadronFlavour());
0041     int pflav = abs(jet.partonFlavour());
0042     int physflav = 0;
0043     if (!(jet.genJet())) {
0044       if (pflav == 0)
0045         return 999;
0046       else
0047         return 1000;
0048     }
0049     if (jet.genParton())
0050       physflav = abs(jet.genParton()->pdgId());
0051     std::size_t nbs = jet.jetFlavourInfo().getbHadrons().size();
0052     std::size_t ncs = jet.jetFlavourInfo().getcHadrons().size();
0053 
0054     unsigned int nbFromGSP(0);
0055     for (const reco::GenParticle& p : gToBB) {
0056       double dr2(reco::deltaR2(jet, p));
0057       if (dr2 < jetR_ * jetR_)
0058         ++nbFromGSP;
0059     }
0060 
0061     unsigned int ncFromGSP(0);
0062     for (const reco::GenParticle& p : gToCC) {
0063       double dr2(reco::deltaR2(jet, p));
0064       if (dr2 < jetR_ * jetR_)
0065         ++ncFromGSP;
0066     }
0067 
0068     //std::cout << " jet pt = " << jet.pt() << " hfl = " << hflav << " pfl = " << pflav << " genpart = " << physflav
0069     //  << " nbFromGSP = " << nbFromGSP << " ncFromGSP = " << ncFromGSP
0070     //  << " nBhadrons " << nbs << " nCHadrons " << ncs << std::endl;
0071     if (hflav == 5) {  //B jet
0072       if (nbs > 1) {
0073         if (nbFromGSP > 0)
0074           return 511;
0075         else
0076           return 510;
0077       } else if (nbs == 1) {
0078         for (std::vector<reco::GenParticle>::const_iterator it = neutrinosLepB.begin(); it != neutrinosLepB.end();
0079              ++it) {
0080           if (reco::deltaR2(it->eta(), it->phi(), jet.eta(), jet.phi()) < 0.4 * 0.4) {
0081             return 520;
0082           }
0083         }
0084         for (std::vector<reco::GenParticle>::const_iterator it = neutrinosLepB_C.begin(); it != neutrinosLepB_C.end();
0085              ++it) {
0086           if (reco::deltaR2(it->eta(), it->phi(), jet.eta(), jet.phi()) < 0.4 * 0.4) {
0087             return 521;
0088           }
0089         }
0090         return 500;
0091       } else {
0092         if (usePhysForLightAndUndefined) {
0093           if (physflav == 21)
0094             return 0;
0095           else if (physflav == 3)
0096             return 2;
0097           else if (physflav == 2 || physflav == 1)
0098             return 1;
0099           else
0100             return 1000;
0101         } else
0102           return 1000;
0103       }
0104     } else if (hflav == 4) {  //C jet
0105       if (ncs > 1) {
0106         if (ncFromGSP > 0)
0107           return 411;
0108         else
0109           return 410;
0110       } else
0111         return 400;
0112     } else {                   //not a heavy jet
0113       if (!alltaus.empty()) {  //check for tau in a simplistic way
0114         bool ishadrtaucontained = true;
0115         for (const auto& p : alltaus) {
0116           size_t ndau = p.numberOfDaughters();
0117           for (size_t i = 0; i < ndau; i++) {
0118             const reco::Candidate* dau = p.daughter(i);
0119             int daupid = std::abs(dau->pdgId());
0120             if (daupid == 13 || daupid == 11) {
0121               ishadrtaucontained = false;
0122               break;
0123             }
0124             if (daupid != 12 && daupid != 14 && daupid != 16 && reco::deltaR2(*dau, jet) > jetR_ * jetR_) {
0125               ishadrtaucontained = false;
0126               break;
0127             }
0128           }
0129         }
0130         if (ishadrtaucontained)
0131           return 600;
0132       }
0133       if (std::abs(pflav) == 4 || std::abs(pflav) == 5 || nbs || ncs) {
0134         if (usePhysForLightAndUndefined) {
0135           if (physflav == 21)
0136             return 0;
0137           else if (physflav == 3)
0138             return 2;
0139           else if (physflav == 2 || physflav == 1)
0140             return 1;
0141           else
0142             return 1000;
0143         } else
0144           return 1000;
0145       } else if (usePhysForLightAndUndefined) {
0146         if (physflav == 21)
0147           return 0;
0148         else if (physflav == 3)
0149           return 2;
0150         else if (physflav == 2 || physflav == 1)
0151           return 1;
0152         else
0153           return 1000;
0154       } else {
0155         if (pflav == 21)
0156           return 0;
0157         else if (pflav == 3)
0158           return 2;
0159         else if (pflav == 2 || pflav == 1)
0160           return 1;
0161         else
0162           return 1000;
0163       }
0164     }
0165   }
0166 
0167   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0168     edm::ParameterSetDescription desc;
0169     desc.add<edm::InputTag>("src")->setComment("input Jet collection");
0170     desc.add<edm::InputTag>("genparticles")->setComment("input genparticles info collection");
0171     desc.add<std::string>("name")->setComment("name of the genJet FlatTable we are extending with flavour information");
0172     descriptions.add("btvMCTable", desc);
0173   }
0174 
0175 private:
0176   void produce(edm::Event&, edm::EventSetup const&) override;
0177 
0178   std::string name_;
0179 
0180   edm::EDGetTokenT<std::vector<pat::Jet> > src_;
0181   constexpr static double jetR_ = 0.4;
0182 
0183   constexpr static bool usePhysForLightAndUndefined = false;
0184   edm::EDGetTokenT<reco::GenParticleCollection> genParticlesToken_;
0185 };
0186 
0187 // ------------ method called to produce the data  ------------
0188 void BTVMCFlavourTableProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0189   auto jets = iEvent.getHandle(src_);
0190   // const auto& jetFlavourInfosProd = iEvent.get(genParticlesToken_);
0191   edm::Handle<reco::GenParticleCollection> genParticlesHandle;
0192   iEvent.getByToken(genParticlesToken_, genParticlesHandle);
0193   std::vector<reco::GenParticle> neutrinosLepB;
0194   std::vector<reco::GenParticle> neutrinosLepB_C;
0195 
0196   std::vector<reco::GenParticle> gToBB;
0197   std::vector<reco::GenParticle> gToCC;
0198   std::vector<reco::GenParticle> alltaus;
0199 
0200   unsigned int nJets = jets->size();
0201 
0202   std::vector<unsigned> jet_FlavSplit(nJets);
0203   for (const reco::Candidate& genC : *genParticlesHandle) {
0204     const reco::GenParticle& gen = static_cast<const reco::GenParticle&>(genC);
0205     if (abs(gen.pdgId()) == 12 || abs(gen.pdgId()) == 14 || abs(gen.pdgId()) == 16) {
0206       const reco::GenParticle* mother = static_cast<const reco::GenParticle*>(gen.mother());
0207       if (mother != nullptr) {
0208         if ((abs(mother->pdgId()) > 500 && abs(mother->pdgId()) < 600) ||
0209             (abs(mother->pdgId()) > 5000 && abs(mother->pdgId()) < 6000)) {
0210           neutrinosLepB.emplace_back(gen);
0211         }
0212         if ((abs(mother->pdgId()) > 400 && abs(mother->pdgId()) < 500) ||
0213             (abs(mother->pdgId()) > 4000 && abs(mother->pdgId()) < 5000)) {
0214           neutrinosLepB_C.emplace_back(gen);
0215         }
0216       } else {
0217         std::cout << "No mother" << std::endl;
0218       }
0219     }
0220 
0221     int id(std::abs(gen.pdgId()));
0222     int status(gen.status());
0223 
0224     if (id == 21 && status >= 21 && status <= 59) {  //// Pythia8 hard scatter, ISR, or FSR
0225       if (gen.numberOfDaughters() == 2) {
0226         const reco::Candidate* d0 = gen.daughter(0);
0227         const reco::Candidate* d1 = gen.daughter(1);
0228         if (std::abs(d0->pdgId()) == 5 && std::abs(d1->pdgId()) == 5 && d0->pdgId() * d1->pdgId() < 0 &&
0229             reco::deltaR2(*d0, *d1) < jetR_ * jetR_)
0230           gToBB.push_back(gen);
0231         if (std::abs(d0->pdgId()) == 4 && std::abs(d1->pdgId()) == 4 && d0->pdgId() * d1->pdgId() < 0 &&
0232             reco::deltaR2(*d0, *d1) < jetR_ * jetR_)
0233           gToCC.push_back(gen);
0234       }
0235     }
0236 
0237     if (id == 15 && false) {
0238       alltaus.push_back(gen);
0239     }
0240   }
0241   for (unsigned i_jet = 0; i_jet < nJets; ++i_jet) {
0242     // from DeepNTuples
0243     const auto& jet = jets->at(i_jet);
0244 
0245     jet_FlavSplit[i_jet] =
0246         jet_flavour(jet, gToBB, gToCC, neutrinosLepB, neutrinosLepB_C, alltaus, usePhysForLightAndUndefined);
0247   }
0248   auto newtab = std::make_unique<nanoaod::FlatTable>(nJets, name_, false, true);
0249   newtab->addColumn<int>("FlavSplit",
0250                          jet_FlavSplit,
0251                          "Flavour of the jet, numerical codes: "
0252                          "isG: 0, "
0253                          "isUD: 1, "
0254                          "isS: 2, "
0255                          "isC: 400, "
0256                          "isCC: 410, "
0257                          "isGCC: 411, "
0258                          "isB: 500, "
0259                          "isBB: 510, "
0260                          "isGBB: 511, "
0261                          "isLeptonicB: 520, "
0262                          "isLeptonicB_C: 521, "
0263                          "isTAU: 600, "
0264                          "isPU: 999,"
0265                          "isUndefined: 1000. "
0266                          "May be combined to form coarse labels for tagger training and flavour dependent attacks "
0267                          "using the loss surface.");
0268   iEvent.put(std::move(newtab));
0269 }
0270 
0271 #include "FWCore/Framework/interface/MakerMacros.h"
0272 //define this as a plug-in
0273 DEFINE_FWK_MODULE(BTVMCFlavourTableProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!