File indexing completed on 2024-04-06 12:23:33
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
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
0078
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
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
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
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
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235
0236 for (size_t m = 0; m < wv.particles->size() && !foundPriority; ++m) {
0237
0238 const GenParticle& aParton = *(wv.particles->at(m).get());
0239 int flavour = abs(aParton.pdgId());
0240
0241
0242
0243
0244 if (doPriority) {
0245 vector<int>::const_iterator ipriority = find(priorityList.begin(), priorityList.end(), flavour);
0246
0247 if (ipriority != priorityList.end()) {
0248
0249
0250 double dist = DeltaR(theJet.p4(), aParton.p4());
0251 if (dist <= coneSizeToAssociate) {
0252 tempParticle = m;
0253 foundPriority = true;
0254 }
0255 }
0256 }
0257
0258
0259 else {
0260 foundPriority = false;
0261 }
0262
0263
0264
0265
0266
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;
0289 wv.theHardest = -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
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323 int JetPartonMatcher::fillPhysicsDefinition(const Jet& theJet, WorkingVariables& wv) const {
0324 float TheBiggerConeSize = 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
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352
0353 for (size_t m = 0; m < wv.particles->size() && !foundPriority; ++m) {
0354
0355 const GenParticle& aParticle = *(wv.particles->at(m).get());
0356 int flavour = abs(aParticle.pdgId());
0357
0358
0359
0360
0361 if (doPriority) {
0362 vector<int>::const_iterator ipriority = find(priorityList.begin(), priorityList.end(), flavour);
0363
0364 if (ipriority != priorityList.end()) {
0365
0366
0367 double dist = DeltaR(theJet.p4(), aParticle.p4());
0368 if (dist <= coneSizeToAssociate) {
0369 tempParticle = m;
0370 foundPriority = true;
0371 }
0372 }
0373 }
0374
0375
0376 else {
0377 foundPriority = false;
0378 }
0379
0380
0381 if (!foundPriority) {
0382
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
0395 tempParticle = m;
0396 nInTheCone++;
0397 }
0398
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
0409 if (!foundPriority) {
0410 wv.theNearest3 = tempNearest;
0411
0412 if (nInTheCone != 1)
0413 return -1;
0414 if (theContaminations.empty())
0415 return tempParticle;
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;
0423 if (initialPartonFlavour == 4) {
0424 if (contaminatingFlavour == 4)
0425 continue;
0426 tempParticle = -1;
0427 return tempParticle;
0428 }
0429 }
0430 }
0431
0432 else {
0433 wv.theHeaviest = tempParticle;
0434 wv.theNearest2 = -1;
0435 wv.theNearest3 = -1;
0436 wv.theHardest = -1;
0437 }
0438
0439 return tempParticle;
0440 }
0441
0442
0443 DEFINE_FWK_MODULE(JetPartonMatcher);