File indexing completed on 2024-09-07 04:37:19
0001
0002 #include <memory>
0003
0004
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
0069
0070
0071 if (hflav == 5) {
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) {
0105 if (ncs > 1) {
0106 if (ncFromGSP > 0)
0107 return 411;
0108 else
0109 return 410;
0110 } else
0111 return 400;
0112 } else {
0113 if (!alltaus.empty()) {
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
0188 void BTVMCFlavourTableProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0189 auto jets = iEvent.getHandle(src_);
0190
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) {
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
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
0273 DEFINE_FWK_MODULE(BTVMCFlavourTableProducer);