JetMCAlgos/plugins/JetFlavourIdentifier.cc

0001 //
0002 // Translation of BTag MCJetFlavour tool to identify real flavour of a jet
0003 // work with CaloJet objects
0004 // Store Infos by Values in JetFlavour.h
0005 // Author: Attilio
0006 // Date: 05.10.2007
0007 //
0008 //
0009 // \class JetFlavourIdentifier
0010 //
0011 // \brief Interface to pull out the proper flavour identifier from a
0012 //        jet->parton matching collection.
0013 //
0014 // In detail, the definitions are as follows:
0015 //
0016 // Definitions:
0017 // The default behavior is that the "definition" is NULL_DEF,
0018 // so the software will fall back to either "physics" or "algorithmic" definition
0019 // as per the "physDefinition" switch.
0020 //
0021 // If the user specifies "definition", then that definition is taken. However,
0022 // if the requested definition is not defined, the software reverts back to
0023 // either "physics" or "algorithmic" based on the "physDefinition" switch.
0024 // For example, if the user specifies "heaviest" as a flavor ID, and there
0025 // are no bottom, charm, or top quarks in the event that match to the jet,
0026 // then the software will fall back to the "physics" or "algorithmic" definition.
0027 //
0028 // Modifications:
0029 //
0030 //     09.03.2008: Sal Rappoccio.
0031 //                 Added capability for all methods of identification, not just
0032 //                 "physics" or "algorithmic". If the requested method does not exist
0033 //                 (i.e. is unphysical), the "physics" or "algorithmic" definition
0034 //                 is defaulted to.
0035
0036 //=======================================================================
0037
0038 // user include files
0039 #include "FWCore/Framework/interface/global/EDProducer.h"
0040 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0041 #include "FWCore/ParameterSet/interface/ParameterSetfwd.h"
0042 #include "FWCore/Utilities/interface/InputTag.h"
0043
0044 #include "FWCore/Framework/interface/Event.h"
0045 #include "FWCore/Framework/interface/EventSetup.h"
0046 #include "FWCore/Framework/interface/MakerMacros.h"
0047 #include "FWCore/Framework/interface/ESHandle.h"
0048 #include "FWCore/Framework/interface/makeRefToBaseProdFrom.h"
0049 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0050
0051 //#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0052
0053 #include "DataFormats/JetReco/interface/Jet.h"
0054 #include "DataFormats/JetReco/interface/JetCollection.h"
0055 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0056 #include "DataFormats/JetMatching/interface/JetFlavour.h"
0057 #include "DataFormats/JetMatching/interface/JetFlavourMatching.h"
0058 #include "DataFormats/JetMatching/interface/MatchedPartons.h"
0059 #include "DataFormats/JetMatching/interface/JetMatchedPartons.h"
0060
0061 #include "DataFormats/Common/interface/Ref.h"
0062 #include "DataFormats/Candidate/interface/Candidate.h"
0063 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0064
0065 #include "DataFormats/Math/interface/Point3D.h"
0066 #include "DataFormats/Math/interface/LorentzVector.h"
0067
0068 #include <memory>
0069 #include <string>
0070 #include <iostream>
0071 #include <vector>
0072 #include <Math/VectorUtil.h>
0073 #include <TMath.h>
0074
0075 using namespace std;
0076 using namespace reco;
0077 using namespace edm;
0078 using namespace ROOT::Math::VectorUtil;
0079
0080 namespace reco {
0081   namespace modules {
0082
0083     //--------------------------------------------------------------------------
0084     //
0085     //--------------------------------------------------------------------------
0086     class JetFlavourIdentifier : public edm::global::EDProducer<> {
0087     public:
0088       enum DEFINITION_T { PHYSICS = 0, ALGO, NEAREST_STATUS2, NEAREST_STATUS3, HEAVIEST, N_DEFINITIONS, NULL_DEF };
0089
0090       JetFlavourIdentifier(const edm::ParameterSet&);
0091       ~JetFlavourIdentifier() override;
0092
0093     private:
0094       void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0095
0096       JetFlavour::Leptons findLeptons(const GenParticleRef&) const;
0097       std::vector<const reco::Candidate*> findCandidates(const reco::Candidate*,
0098                                                          int,
0099                                                          math::XYZTLorentzVector const& thePartonLV) const;
0100       void fillLeptons(const std::vector<const reco::Candidate*>&,
0101                        JetFlavour::Leptons&,
0102                        int,
0103                        int,
0104                        math::XYZTLorentzVector const& thePartonLV) const;
0105       static int heaviestFlavour(int);
0106
0107       EDGetTokenT<JetMatchedPartonsCollection> sourceByReferToken_;
0108       bool physDefinition;
0109       bool leptonInfo_;
0110       DEFINITION_T definition;
0111     };
0112   }  // namespace modules
0113 }  // namespace reco
0114 using reco::modules::JetFlavourIdentifier;
0115
0116 //=========================================================================
0117
0118 JetFlavourIdentifier::JetFlavourIdentifier(const edm::ParameterSet& iConfig) {
0119   produces<JetFlavourMatchingCollection>();
0120   sourceByReferToken_ = consumes<JetMatchedPartonsCollection>(iConfig.getParameter<InputTag>("srcByReference"));
0121   physDefinition = iConfig.getParameter<bool>("physicsDefinition");
0122   leptonInfo_ = iConfig.exists("leptonInfo") ? iConfig.getParameter<bool>("leptonInfo") : false;
0123   // If we have a definition of which parton to identify, use it,
0124   // otherwise we default to the "old" behavior of either "physics" or "algorithmic".
0125   // Furthermore, if the specified definition is not sensible for the given jet,
0126   // then the "physDefinition" switch is used to identify the flavour of the jet.
0127   if (iConfig.exists("definition")) {
0128     definition = static_cast<DEFINITION_T>(iConfig.getParameter<int>("definition"));
0129   } else {
0130     definition = NULL_DEF;
0131   }
0132 }
0133
0134 //=========================================================================
0135
0136 JetFlavourIdentifier::~JetFlavourIdentifier() {}
0137
0138 // ------------ method called to produce the data  ------------
0139
0140 void JetFlavourIdentifier::produce(StreamID, Event& iEvent, const EventSetup& iEs) const {
0141   // Get the JetMatchedPartons
0142   Handle<JetMatchedPartonsCollection> theTagByRef;
0143   iEvent.getByToken(sourceByReferToken_, theTagByRef);
0144
0145   // Create a JetFlavourMatchingCollection
0146   JetFlavourMatchingCollection* jfmc;
0147   if (!theTagByRef->empty()) {
0148     RefToBase<Jet> jj = theTagByRef->begin()->first;
0149     jfmc = new JetFlavourMatchingCollection(edm::makeRefToBaseProdFrom(jj, iEvent));
0150   } else {
0151     jfmc = new JetFlavourMatchingCollection();
0152   }
0153   std::unique_ptr<reco::JetFlavourMatchingCollection> jetFlavMatching(jfmc);
0154
0155   // Loop over the matched partons and see which match.
0156   for (JetMatchedPartonsCollection::const_iterator j = theTagByRef->begin(); j != theTagByRef->end(); j++) {
0157     // Consider this match.
0158     const MatchedPartons aMatch = (*j).second;
0159
0160     // This will hold the 4-vector, vertex, flavour and the leptonian decays (0: no lepton, xyz: x leptons in layer 2, y in layer 1 and z in layer 0) of the requested object.
0161     math::XYZTLorentzVector thePartonLorentzVector(0, 0, 0, 0);
0162     math::XYZPoint thePartonVertex(0, 0, 0);
0163     int thePartonFlavour = 0;
0164     JetFlavour::Leptons theLeptons;
0165
0166     // get the partons based on which definition to use.
0167     switch (definition) {
0168       case PHYSICS: {
0169         const GenParticleRef& aPartPhy = aMatch.physicsDefinitionParton();
0170         if (aPartPhy.isNonnull()) {
0171           thePartonLorentzVector = aPartPhy.get()->p4();
0172           thePartonVertex = aPartPhy.get()->vertex();
0173           thePartonFlavour = aPartPhy.get()->pdgId();
0174           if (leptonInfo_)
0175             theLeptons = findLeptons(aPartPhy);
0176         }
0177         break;
0178       }
0179       case ALGO: {
0180         const GenParticleRef& aPartAlg = aMatch.algoDefinitionParton();
0181         if (aPartAlg.isNonnull()) {
0182           thePartonLorentzVector = aPartAlg.get()->p4();
0183           thePartonVertex = aPartAlg.get()->vertex();
0184           thePartonFlavour = aPartAlg.get()->pdgId();
0185           if (leptonInfo_)
0186             theLeptons = findLeptons(aPartAlg);
0187         }
0188         break;
0189       }
0190       case NEAREST_STATUS2: {
0191         const GenParticleRef& aPartN2 = aMatch.nearest_status2();
0192         if (aPartN2.isNonnull()) {
0193           thePartonLorentzVector = aPartN2.get()->p4();
0194           thePartonVertex = aPartN2.get()->vertex();
0195           thePartonFlavour = aPartN2.get()->pdgId();
0196           if (leptonInfo_)
0197             theLeptons = findLeptons(aPartN2);
0198         }
0199         break;
0200       }
0201       case NEAREST_STATUS3: {
0202         const GenParticleRef& aPartN3 = aMatch.nearest_status3();
0203         if (aPartN3.isNonnull()) {
0204           thePartonLorentzVector = aPartN3.get()->p4();
0205           thePartonVertex = aPartN3.get()->vertex();
0206           thePartonFlavour = aPartN3.get()->pdgId();
0207           if (leptonInfo_)
0208             theLeptons = findLeptons(aPartN3);
0209         }
0210         break;
0211       }
0212       case HEAVIEST: {
0213         const GenParticleRef aPartHeaviest = aMatch.heaviest();
0214         if (aPartHeaviest.isNonnull()) {
0215           thePartonLorentzVector = aPartHeaviest.get()->p4();
0216           thePartonVertex = aPartHeaviest.get()->vertex();
0217           thePartonFlavour = aPartHeaviest.get()->pdgId();
0218           if (leptonInfo_)
0219             theLeptons = findLeptons(aPartHeaviest);
0220         }
0221         break;
0222       }
0223       // Default case is backwards-compatible
0224       default: {
0225         if (physDefinition) {
0226           const GenParticleRef& aPartPhy = aMatch.physicsDefinitionParton();
0227           if (aPartPhy.isNonnull()) {
0228             thePartonLorentzVector = aPartPhy.get()->p4();
0229             thePartonVertex = aPartPhy.get()->vertex();
0230             thePartonFlavour = aPartPhy.get()->pdgId();
0231             if (leptonInfo_)
0232               theLeptons = findLeptons(aPartPhy);
0233           }
0234         } else {
0235           const GenParticleRef& aPartAlg = aMatch.algoDefinitionParton();
0236           if (aPartAlg.isNonnull()) {
0237             thePartonLorentzVector = aPartAlg.get()->p4();
0238             thePartonVertex = aPartAlg.get()->vertex();
0239             thePartonFlavour = aPartAlg.get()->pdgId();
0240             if (leptonInfo_)
0241               theLeptons = findLeptons(aPartAlg);
0242           }
0243         }
0244       } break;
0245     }  // end switch on definition
0246
0247     // Now make sure we have a match. If the user specified "heaviest", for instance,
0248     // and there is no b- or c-quarks, then fall back to the "physDefinition" switch.
0249
0250     if (thePartonFlavour == 0) {
0251       if (physDefinition) {
0252         const GenParticleRef& aPartPhy = aMatch.physicsDefinitionParton();
0253         if (aPartPhy.isNonnull()) {
0254           thePartonLorentzVector = aPartPhy.get()->p4();
0255           thePartonVertex = aPartPhy.get()->vertex();
0256           thePartonFlavour = aPartPhy.get()->pdgId();
0257           if (leptonInfo_)
0258             theLeptons = findLeptons(aPartPhy);
0259         }
0260       } else {
0261         const GenParticleRef& aPartAlg = aMatch.algoDefinitionParton();
0262         if (aPartAlg.isNonnull()) {
0263           thePartonLorentzVector = aPartAlg.get()->p4();
0264           thePartonVertex = aPartAlg.get()->vertex();
0265           thePartonFlavour = aPartAlg.get()->pdgId();
0266           if (leptonInfo_)
0267             theLeptons = findLeptons(aPartAlg);
0268         }
0269       }
0270     }
0271
0272     /*
0273      std::cout << "Leptons of " <<thePartonFlavour << " Jet: " << std::endl;
0274      std::cout << "  electrons: " <<theLeptons.electron << std::endl;
0275      std::cout << "  muons    : " <<theLeptons.muon << std::endl;
0276      std::cout << "  tau      : " <<theLeptons.tau << std::endl;
0277 */
0278     // Add the jet->flavour match to the map.
0279     (*jetFlavMatching)[(*j).first] = JetFlavour(thePartonLorentzVector, thePartonVertex, thePartonFlavour, theLeptons);
0280   }  // end loop over jets
0281
0282   // Put the object into the event.
0283   iEvent.put(std::move(jetFlavMatching));
0284 }
0285
0286 JetFlavour::Leptons JetFlavourIdentifier::findLeptons(const GenParticleRef& parton) const {
0287   JetFlavour::Leptons theLeptons;
0288
0289   auto const& thePartonLV = parton->p4();
0290
0291   ///first daughter of the parton should be an MC particle (pdgId==92,93)
0292   const reco::Candidate* mcstring = parton->daughter(0);
0293   int partonFlavour = std::abs(parton->pdgId());
0294   //  std::cout << "parton DeltaR: " << DeltaR(thePartonLV, parton->p4()) << std::endl;
0295
0296   ///lookup particles with parton flavour and weak decay
0297   std::vector<const reco::Candidate*> candidates = findCandidates(mcstring, partonFlavour, parton->p4());
0298   //  std::cout << "Candidates are:" << std::endl;
0299   //  for(unsigned int j = 0; j < candidates.size(); j++) std::cout << "   --> " << candidates[j]->pdgId() << std::endl;
0300
0301   ///count leptons of candidates
0302   fillLeptons(candidates, theLeptons, 1, partonFlavour, thePartonLV);
0303
0304   return theLeptons;
0305 }
0306
0307 std::vector<const reco::Candidate*> JetFlavourIdentifier::findCandidates(
0308     const reco::Candidate* cand, int partonFlavour, math::XYZTLorentzVector const& thePartonLV) const {
0309   std::vector<const reco::Candidate*> cands;
0310   if (!cand)
0311     return cands;
0312
0313   for (unsigned int i = 0; i < cand->numberOfDaughters(); i++) {
0314     /*
0315     std::cout << "DeltaR - " << partonFlavour << " ";
0316     if (DeltaR(thePartonLV, cand->daughter(i)->p4()) > 0.7) std::cout << "(";
0317     std::cout << cand->daughter(i)->pdgId() << ": " << DeltaR(thePartonLV, cand->daughter(i)->p4());
0318     if (DeltaR(thePartonLV, cand->daughter(i)->p4()) > 0.7) std::cout << ")";
0319     std::cout << std::endl;
0320 */
0321     if (DeltaR(thePartonLV, cand->daughter(i)->p4()) < 0.7) {
0322       int pdgId = std::abs(cand->daughter(i)->pdgId());
0323       int flavour = heaviestFlavour(pdgId);
0324       if (flavour == partonFlavour || (flavour >= 10 && partonFlavour >= 10)) {
0325         //        std::cout << "<------- " << std::endl;
0326         std::vector<const reco::Candidate*> newcands = findCandidates(cand->daughter(i), partonFlavour, thePartonLV);
0327         //        std::cout << " ------->" << std::endl;
0328         std::copy(newcands.begin(), newcands.end(), std::back_inserter(cands));
0329       }
0330       if (partonFlavour >= 10)
0331         cands.push_back(cand->daughter(i));
0332     }
0333   }
0334
0335   if (cands.empty() && std::abs(cand->pdgId()) > 110 &&
0336       !(partonFlavour >= 4 && partonFlavour < 10 && heaviestFlavour(cand->pdgId()) < 4))
0337     cands.push_back(cand);
0338
0339   return cands;
0340 }
0341
0342 void JetFlavourIdentifier::fillLeptons(const std::vector<const reco::Candidate*>& cands,
0343                                        JetFlavour::Leptons& leptons,
0344                                        int rank,
0345                                        int flavour,
0346                                        math::XYZTLorentzVector const& thePartonLV) const {
0347   for (unsigned int j = 0; j < cands.size(); j++) {
0348     for (unsigned int i = 0; i < cands[j]->numberOfDaughters(); i++) {
0349       int pdgId = std::abs(cands[j]->daughter(i)->pdgId());
0350
0351       //      for(int z = 1; z <= rank; z *= 10) std::cout << " ------ ";
0352       //      std::cout << pdgId << std::endl;
0353
0354       ///test for neutrinos because of conversions and dalitz pions
0355       if (pdgId == 12)
0356         leptons.electron += rank;
0357       else if (pdgId == 14)
0358         leptons.muon += rank;
0359       else if (pdgId == 16)
0360         leptons.tau += rank;
0361       else {
0362         int heaviest = heaviestFlavour(pdgId);
0363         int heaviest_ = heaviest < 10 ? heaviest : 0;
0364         if (!heaviest || (flavour < 4 ? (heaviest_ < 4) : (heaviest >= 4))) {
0365           std::vector<const reco::Candidate*> newcands = findCandidates(cands[j]->daughter(i), heaviest, thePartonLV);
0366           if (pdgId <= 110)
0367             newcands.push_back(cands[j]->daughter(i));
0368           fillLeptons(newcands, leptons, rank * 10, std::max(heaviest_, flavour), thePartonLV);
0369         }
0370       }
0371     }
0372   }
0373 }
0374
0375 int JetFlavourIdentifier::heaviestFlavour(int pdgId) {
0376   int flavour = 0;
0377
0378   pdgId = std::abs(pdgId) % 100000;
0379   if (pdgId > 110) {
0380     while (pdgId % 10 > 0 && pdgId % 10 < 6) {
0381       pdgId /= 10;
0382       if (pdgId % 10 > flavour)
0383         flavour = pdgId % 10;
0384     }
0385   } else
0386     flavour = pdgId;
0387
0388   return flavour;
0389 }
0390
0391 //define this as a plug-in
0392 DEFINE_FWK_MODULE(JetFlavourIdentifier);