Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-03-24 03:15:00

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 JetPartonMatcher
0010 //
0011 // \brief Interface to create a specified "matched" parton to jet match based
0012 //        on the user interpretation.
0013 //
0014 // Algorithm for definitions:
0015 //
0016 // If the particle is matched to any item in the priority list (given by user),
0017 // then count that prioritized particle as the match.
0018 //
0019 // If not resort to default behavior:
0020 //
0021 //
0022 // 1) Algorithmic Definition:
0023 //
0024 // A particle is a parton if its daughter is a string(code=92) or a cluster(code=93)
0025 // If (parton is within the cone defined by coneSizeToAssociate) then:
0026 //           if (parton is a b)                                   then associatedParton is the b
0027 //      else if (associatedParton =! b and parton is a c)         then associatedParton is the c
0028 //      else if (associatedParton =! b and associatedParton =! c) then associatedParton is the one with the highest pT
0029 // associatedParton can be -1 --> no partons in the cone
0030 // True Flavour of the jet --> flavour of the associatedParton
0031 //
0032 // ToDo: if more than one b(c) in the cone --> the selected parton is not always the one with highest pT
0033 //
0034 //
0035 // 2) Physics Definition:
0036 //
0037 // A particle is a parton if its daughter is a string(code=92) or a cluster(code=93)
0038 // if( only one initialParticle within the cone defined by coneSizeToAssociate) associatedInitialParticle is the initialParticle
0039 // TheBiggerConeSize = 0.7 --> it's hard coded!
0040 // if( a parton not coming from associatedInitialParticle in TheBiggerConeSize is a b or a c) reject the association
0041 //
0042 // associatedInitialParticle can be -1 --> no initialParticle in the cone or rejected association
0043 // True Flavour of the jet --> flavour of the associatedInitialParticle
0044 //
0045 //
0046 // Modifications:
0047 //
0048 //     09.03.2008: Sal Rappoccio.
0049 //                 Added capability for "priority" list.
0050 //
0051 
0052 //=======================================================================
0053 
0054 // user include files
0055 #include "FWCore/Framework/interface/global/EDProducer.h"
0056 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0057 #include "FWCore/ParameterSet/interface/ParameterSetfwd.h"
0058 #include "FWCore/Utilities/interface/InputTag.h"
0059 
0060 #include "FWCore/Framework/interface/Event.h"
0061 #include "FWCore/Framework/interface/EventSetup.h"
0062 #include "FWCore/Framework/interface/MakerMacros.h"
0063 #include "FWCore/Framework/interface/ESHandle.h"
0064 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0065 
0066 #include "DataFormats/JetReco/interface/Jet.h"
0067 #include "DataFormats/JetReco/interface/JetCollection.h"
0068 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0069 #include "DataFormats/JetMatching/interface/JetFlavour.h"
0070 #include "DataFormats/JetMatching/interface/JetFlavourMatching.h"
0071 #include "DataFormats/JetMatching/interface/MatchedPartons.h"
0072 #include "DataFormats/JetMatching/interface/JetMatchedPartons.h"
0073 
0074 #include "DataFormats/Common/interface/View.h"
0075 #include "DataFormats/Common/interface/Ref.h"
0076 #include "DataFormats/Common/interface/getRef.h"
0077 //#include "DataFormats/Candidate/interface/Candidate.h"
0078 //#include "DataFormats/Candidate/interface/CandidateFwd.h"
0079 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0080 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0081 
0082 #include <memory>
0083 #include <string>
0084 #include <iostream>
0085 #include <vector>
0086 #include <Math/VectorUtil.h>
0087 #include <TMath.h>
0088 
0089 using namespace std;
0090 using namespace reco;
0091 using namespace edm;
0092 using namespace ROOT::Math::VectorUtil;
0093 
0094 class JetPartonMatcher : public edm::global::EDProducer<> {
0095 public:
0096   JetPartonMatcher(const edm::ParameterSet&);
0097   ~JetPartonMatcher() override;
0098 
0099 private:
0100   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0101 
0102   struct WorkingVariables {
0103     Handle<GenParticleRefVector> particles;
0104     int theHeaviest = 0;
0105     int theNearest2 = 0;
0106     int theNearest3 = 0;
0107     int theHardest = 0;
0108   };
0109 
0110   int fillAlgoritDefinition(const Jet&, WorkingVariables&) const;
0111   int fillPhysicsDefinition(const Jet&, WorkingVariables&) const;
0112 
0113   edm::EDGetTokenT<edm::View<reco::Jet> > m_jetsSrcToken;
0114   edm::EDGetTokenT<GenParticleRefVector> m_ParticleSrcToken;
0115   double coneSizeToAssociate;
0116   bool physDefinition;
0117   bool doPriority;
0118   vector<int> priorityList;
0119 };
0120 
0121 //=========================================================================
0122 
0123 JetPartonMatcher::JetPartonMatcher(const edm::ParameterSet& iConfig) {
0124   produces<JetMatchedPartonsCollection>();
0125   m_jetsSrcToken = consumes<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>("jets"));
0126   m_ParticleSrcToken = consumes<GenParticleRefVector>(iConfig.getParameter<edm::InputTag>("partons"));
0127   coneSizeToAssociate = iConfig.getParameter<double>("coneSizeToAssociate");
0128   if (iConfig.exists("doPriority")) {
0129     doPriority = iConfig.getParameter<bool>("doPriority");
0130     priorityList = iConfig.getParameter<vector<int> >("priorityList");
0131   } else {
0132     doPriority = false;
0133     priorityList.clear();
0134   }
0135 }
0136 
0137 //=========================================================================
0138 
0139 JetPartonMatcher::~JetPartonMatcher() {}
0140 
0141 // ------------ method called to produce the data  ------------
0142 
0143 void JetPartonMatcher::produce(edm::StreamID, Event& iEvent, const EventSetup& iEs) const {
0144   WorkingVariables wv;
0145   edm::Handle<edm::View<reco::Jet> > jets_h;
0146   iEvent.getByToken(m_jetsSrcToken, jets_h);
0147   iEvent.getByToken(m_ParticleSrcToken, wv.particles);
0148 
0149   edm::LogVerbatim("JetPartonMatcher") << "=== Partons size:" << wv.particles->size();
0150 
0151   for (size_t m = 0; m < wv.particles->size(); ++m) {
0152     const GenParticle& aParton = *(wv.particles->at(m).get());
0153     edm::LogVerbatim("JetPartonMatcher") << aParton.status() << " " << aParton.pdgId() << " " << aParton.pt() << " "
0154                                          << aParton.eta() << " " << aParton.phi() << endl;
0155   }
0156 
0157   auto jetMatchedPartons = std::make_unique<JetMatchedPartonsCollection>(reco::JetRefBaseProd(jets_h));
0158 
0159   for (size_t j = 0; j < jets_h->size(); j++) {
0160     const int theMappedPartonAlg = fillAlgoritDefinition((*jets_h)[j], wv);
0161     const int theMappedPartonPhy = fillPhysicsDefinition((*jets_h)[j], wv);
0162 
0163     GenParticleRef pHV;
0164     GenParticleRef pN2;
0165     GenParticleRef pN3;
0166     GenParticleRef pPH;
0167     GenParticleRef pAL;
0168 
0169     if (wv.theHeaviest >= 0)
0170       pHV = wv.particles->at(wv.theHeaviest);
0171     if (wv.theNearest2 >= 0)
0172       pN2 = wv.particles->at(wv.theNearest2);
0173     if (wv.theNearest3 >= 0)
0174       pN3 = wv.particles->at(wv.theNearest3);
0175     if (theMappedPartonPhy >= 0)
0176       pPH = wv.particles->at(theMappedPartonPhy);
0177     if (theMappedPartonAlg >= 0)
0178       pAL = wv.particles->at(theMappedPartonAlg);
0179 
0180     (*jetMatchedPartons)[jets_h->refAt(j)] = MatchedPartons(pHV, pN2, pN3, pPH, pAL);
0181   }
0182 
0183   iEvent.put(std::move(jetMatchedPartons));
0184 }
0185 
0186 //
0187 // Algorithmic Definition:
0188 // Output: define one associatedParton
0189 // Loop on all particle.
0190 //
0191 // (Note: This part added by Salvatore Rappoccio)
0192 // [begin]
0193 // If the particle is matched to any item in the priority list (given by user),
0194 // then count that prioritized particle as the match.
0195 //
0196 // If not resort to default behavior:
0197 // [end]
0198 //
0199 // A particle is a parton if its daughter is a string(code=92) or a cluster(code=93)
0200 // If (parton is within the cone defined by coneSizeToAssociate) then:
0201 //           if (parton is a b)                                   then associatedParton is the b
0202 //      else if (associatedParton =! b and parton is a c)         then associatedParton is the c
0203 //      else if (associatedParton =! b and associatedParton =! c) then associatedParton is the one with the highest pT
0204 // associatedParton can be -1 --> no partons in the cone
0205 // True Flavour of the jet --> flavour of the associatedParton
0206 //
0207 // ToDo: if more than one b(c) in the cone --> the selected parton is not always the one with highest pT
0208 //
0209 int JetPartonMatcher::fillAlgoritDefinition(const Jet& theJet, WorkingVariables& wv) const {
0210   int tempParticle = -1;
0211   int tempPartonHighestPt = -1;
0212   int tempNearest = -1;
0213   float maxPt = 0;
0214   float minDr = 1000;
0215   bool foundPriority = false;
0216 
0217   // Loop over the particles in question until we find a priority
0218   // "hit", or if we find none in the priority list (or we don't want
0219   // to consider priority), then loop through all particles and fill
0220   // standard definition.
0221 
0222   //
0223   // Matching:
0224   //
0225   // 1) First try to match by hand. The "priority list" is given
0226   //    by the user. The algorithm finds any particles in that
0227   //    "priority list" that are within the specified cone size.
0228   //    If it finds one, it counts the object as associating to that
0229   //    particle.
0230   //    NOTE! The objects in the priority list are given in order
0231   //    of priority. So if the user specifies:
0232   //    6, 24, 21,
0233   //    then it will first look for top quarks, then W bosons, then gluons.
0234   // 2) If no priority items are found, do the default "standard"
0235   //    matching.
0236   for (size_t m = 0; m < wv.particles->size() && !foundPriority; ++m) {
0237     //    const Candidate & aParton = *(wv.particles->at(m).get());
0238     const GenParticle& aParton = *(wv.particles->at(m).get());
0239     int flavour = abs(aParton.pdgId());
0240 
0241     // "Priority" behavoir:
0242     // Associate to the first particle found in the priority list, regardless
0243     // of delta R.
0244     if (doPriority) {
0245       vector<int>::const_iterator ipriority = find(priorityList.begin(), priorityList.end(), flavour);
0246       // Check to see if this particle is in our priority list
0247       if (ipriority != priorityList.end()) {
0248         // This particle is on our priority list. If it matches,
0249         // we break, since we take in order of priority, not deltaR
0250         double dist = DeltaR(theJet.p4(), aParton.p4());
0251         if (dist <= coneSizeToAssociate) {
0252           tempParticle = m;
0253           foundPriority = true;
0254         }
0255       }
0256     }
0257     // Here we do not want to do priority matching. Ensure the "foundPriority" swtich
0258     // is turned off.
0259     else {
0260       foundPriority = false;
0261     }
0262 
0263     // Default behavior:
0264     // Look for partons before the color string to associate.
0265     // Do this if we don't want to do priority matching, or if
0266     // we didn't find a priority object.
0267 
0268     if (!foundPriority && aParton.status() != 3) {
0269       double dist = DeltaR(theJet.p4(), aParton.p4());
0270       if (dist <= coneSizeToAssociate) {
0271         if (dist < minDr) {
0272           minDr = dist;
0273           tempNearest = m;
0274         }
0275         if (tempParticle == -1 && (flavour == 4))
0276           tempParticle = m;
0277         if (flavour == 5)
0278           tempParticle = m;
0279         if (aParton.pt() > maxPt) {
0280           maxPt = aParton.pt();
0281           tempPartonHighestPt = m;
0282         }
0283       }
0284     }
0285   }
0286 
0287   if (foundPriority) {
0288     wv.theHeaviest = tempParticle;  // The priority-matched particle
0289     wv.theHardest = -1;             //  set the rest to -1
0290     wv.theNearest2 = -1;            // "                  "
0291   } else {
0292     wv.theHeaviest = tempParticle;
0293     wv.theHardest = tempPartonHighestPt;
0294     wv.theNearest2 = tempNearest;
0295     if (tempParticle == -1)
0296       tempParticle = tempPartonHighestPt;
0297   }
0298   return tempParticle;
0299 }
0300 
0301 //
0302 // Physics Definition:
0303 // A initialParticle is a particle with status=3
0304 // Output: define one associatedInitialParticle
0305 // Loop on all particles
0306 //
0307 // (Note: This part added by Salvatore Rappoccio)
0308 // [begin]
0309 // If the particle is matched to any item in the priority list (given by user),
0310 // then count that prioritized particle as the match.
0311 //
0312 // If not resort to default behavior:
0313 // [end]
0314 //
0315 // A particle is a parton if its daughter is a string(code=92) or a cluster(code=93)
0316 // if( only one initialParticle within the cone defined by coneSizeToAssociate) associatedInitialParticle is the initialParticle
0317 // TheBiggerConeSize = 0.7 --> it's hard coded!
0318 // if( a parton not coming from associatedInitialParticle in TheBiggerConeSize is a b or a c) reject the association
0319 //
0320 // associatedInitialParticle can be -1 --> no initialParticle in the cone or rejected association
0321 // True Flavour of the jet --> flavour of the associatedInitialParticle
0322 //
0323 int JetPartonMatcher::fillPhysicsDefinition(const Jet& theJet, WorkingVariables& wv) const {
0324   float TheBiggerConeSize = 0.7;  // In HepMC it's 0.3 --> it's a mistake: value has to be 0.7
0325   int tempParticle = -1;
0326   int nInTheCone = 0;
0327   int tempNearest = -1;
0328   float minDr = 1000;
0329   bool foundPriority = false;
0330 
0331   vector<const reco::Candidate*> theContaminations;
0332   theContaminations.clear();
0333 
0334   // Loop over the particles in question until we find a priority
0335   // "hit", or if we find none in the priority list (or we don't want
0336   // to consider priority), then loop through all particles and fill
0337   // standard definition.
0338 
0339   //
0340   // Matching:
0341   //
0342   // 1) First try to match by hand. The "priority list" is given
0343   //    by the user. The algorithm finds any particles in that
0344   //    "priority list" that are within the specified cone size.
0345   //    If it finds one, it counts the object as associating to that
0346   //    particle.
0347   //    NOTE! The objects in the priority list are given in order
0348   //    of priority. So if the user specifies:
0349   //    6, 24, 21,
0350   //    then it will first look for top quarks, then W bosons, then gluons.
0351   // 2) If no priority items are found, do the default "standard"
0352   //    matching.
0353   for (size_t m = 0; m < wv.particles->size() && !foundPriority; ++m) {
0354     //    const Candidate & aParticle = *(wv.particles->at(m).get());
0355     const GenParticle& aParticle = *(wv.particles->at(m).get());
0356     int flavour = abs(aParticle.pdgId());
0357 
0358     // "Priority" behavoir:
0359     // Associate to the first particle found in the priority list, regardless
0360     // of delta R.
0361     if (doPriority) {
0362       vector<int>::const_iterator ipriority = find(priorityList.begin(), priorityList.end(), flavour);
0363       // Check to see if this particle is in our priority list
0364       if (ipriority != priorityList.end()) {
0365         // This particle is on our priority list. If it matches,
0366         // we break, since we take in order of priority, not deltaR
0367         double dist = DeltaR(theJet.p4(), aParticle.p4());
0368         if (dist <= coneSizeToAssociate) {
0369           tempParticle = m;
0370           foundPriority = true;
0371         }
0372       }
0373     }
0374     // Here we do not want to do priority matching. Ensure the "foundPriority" swtich
0375     // is turned off.
0376     else {
0377       foundPriority = false;
0378     }
0379 
0380     // Default behavior:
0381     if (!foundPriority) {
0382       // skipping all particle but udscbg (is this correct/enough?!?!)
0383       bool isAParton = false;
0384       if (flavour == 1 || flavour == 2 || flavour == 3 || flavour == 4 || flavour == 5 || flavour == 21)
0385         isAParton = true;
0386       if (!isAParton)
0387         continue;
0388       double dist = DeltaR(theJet.p4(), aParticle.p4());
0389       if (aParticle.status() == 3 && dist < minDr) {
0390         minDr = dist;
0391         tempNearest = m;
0392       }
0393       if (aParticle.status() == 3 && dist <= coneSizeToAssociate) {
0394         //cout << "particle in small cone=" << aParticle.pdgId() << endl;
0395         tempParticle = m;
0396         nInTheCone++;
0397       }
0398       // Look for heavy partons in TheBiggerConeSize now
0399       if (aParticle.numberOfDaughters() > 0 && aParticle.status() != 3) {
0400         if (flavour == 1 || flavour == 2 || flavour == 3 || flavour == 21)
0401           continue;
0402         if (dist < TheBiggerConeSize)
0403           theContaminations.push_back(&aParticle);
0404       }
0405     }
0406   }
0407 
0408   // Here's the default behavior for assignment if there is no priority.
0409   if (!foundPriority) {
0410     wv.theNearest3 = tempNearest;
0411 
0412     if (nInTheCone != 1)
0413       return -1;  // rejected --> only one initialParton requested
0414     if (theContaminations.empty())
0415       return tempParticle;  //no contamination
0416     int initialPartonFlavour = abs((wv.particles->at(tempParticle).get())->pdgId());
0417 
0418     vector<const Candidate*>::const_iterator itCont = theContaminations.begin();
0419     for (; itCont != theContaminations.end(); itCont++) {
0420       int contaminatingFlavour = abs((*itCont)->pdgId());
0421       if ((*itCont)->numberOfMothers() > 0 && (*itCont)->mother(0) == wv.particles->at(tempParticle).get())
0422         continue;  // mother is the initialParton --> OK
0423       if (initialPartonFlavour == 4) {
0424         if (contaminatingFlavour == 4)
0425           continue;         // keep association --> the initialParton is a c --> the contaminated parton is a c
0426         tempParticle = -1;  // all the other cases reject!
0427         return tempParticle;
0428       }
0429     }
0430   }
0431   // If there is priority, then just set the heaviest to priority, the rest are -1.
0432   else {
0433     wv.theHeaviest = tempParticle;  // Set the heaviest to tempParticle
0434     wv.theNearest2 = -1;            //  Set the rest to -1
0435     wv.theNearest3 = -1;            // "                  "
0436     wv.theHardest = -1;             // "                  "
0437   }
0438 
0439   return tempParticle;
0440 }
0441 
0442 //define this as a plug-in
0443 DEFINE_FWK_MODULE(JetPartonMatcher);