File indexing completed on 2023-03-17 11:15:42
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 #include "FWCore/Framework/interface/global/EDProducer.h"
0031 #include <string>
0032 #include <vector>
0033 #include <set>
0034 #include <utility>
0035 #include <algorithm>
0036
0037 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0038 #include "DataFormats/Common/interface/Ptr.h"
0039 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0040 #include "DataFormats/Candidate/interface/ShallowClonePtrCandidate.h"
0041 #include "DataFormats/HepMCCandidate/interface/FlavorHistory.h"
0042 #include "DataFormats/HepMCCandidate/interface/FlavorHistoryEvent.h"
0043 #include "DataFormats/Common/interface/Handle.h"
0044 #include "FWCore/Framework/interface/Event.h"
0045 #include "FWCore/Framework/interface/EventSetup.h"
0046 #include "FWCore/Utilities/interface/EDMException.h"
0047 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0048 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0049 #include <fstream>
0050
0051 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0052 #include "DataFormats/Math/interface/deltaR.h"
0053
0054 class FlavorHistoryProducer : public edm::global::EDProducer<> {
0055 public:
0056
0057 FlavorHistoryProducer(const edm::ParameterSet &);
0058
0059 private:
0060
0061 void produce(edm::StreamID, edm::Event &e, const edm::EventSetup &) const override;
0062
0063 void getAncestors(const reco::Candidate &c, std::vector<reco::Candidate const *> &moms) const;
0064
0065 reco::CandidateView::const_iterator getClosestJet(edm::Handle<reco::CandidateView> const &pJets,
0066 reco::Candidate const &parton) const;
0067
0068 const edm::EDGetTokenT<reco::CandidateView> srcToken_;
0069 const edm::EDGetTokenT<reco::CandidateView> matchedSrcToken_;
0070 const double matchDR_;
0071 const int pdgIdToSelect_;
0072 const double ptMinParticle_;
0073 const double ptMinShower_;
0074 const double etaMaxParticle_;
0075 const double etaMaxShower_;
0076 const std::string flavorHistoryName_;
0077 const bool verbose_;
0078 };
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104 using namespace std;
0105 using namespace reco;
0106 using namespace edm;
0107
0108 FlavorHistoryProducer::FlavorHistoryProducer(const ParameterSet &p)
0109 : srcToken_(consumes<CandidateView>(p.getParameter<InputTag>("src"))),
0110 matchedSrcToken_(consumes<CandidateView>(p.getParameter<InputTag>("matchedSrc"))),
0111 matchDR_(p.getParameter<double>("matchDR")),
0112 pdgIdToSelect_(p.getParameter<int>("pdgIdToSelect")),
0113 ptMinParticle_(p.getParameter<double>("ptMinParticle")),
0114 ptMinShower_(p.getParameter<double>("ptMinShower")),
0115 etaMaxParticle_(p.getParameter<double>("etaMaxParticle")),
0116 etaMaxShower_(p.getParameter<double>("etaMaxShower")),
0117 flavorHistoryName_(p.getParameter<string>("flavorHistoryName")),
0118 verbose_(p.getUntrackedParameter<bool>("verbose")) {
0119 produces<FlavorHistoryEvent>(flavorHistoryName_);
0120 }
0121
0122 void FlavorHistoryProducer::produce(StreamID, Event &evt, const EventSetup &) const {
0123 if (verbose_)
0124 cout << "---------- Hello from FlavorHistoryProducer! -----" << endl;
0125
0126
0127 Handle<CandidateView> particlesViewH;
0128 evt.getByToken(srcToken_, particlesViewH);
0129
0130 Handle<CandidateView> matchedView;
0131 evt.getByToken(matchedSrcToken_, matchedView);
0132
0133
0134 vector<const Candidate *> particles;
0135 for (CandidateView::const_iterator p = particlesViewH->begin(); p != particlesViewH->end(); ++p) {
0136 particles.push_back(&*p);
0137 }
0138
0139
0140 vector<int> partonIndices;
0141
0142 vector<int> progenitorIndices;
0143
0144 vector<int> sisterIndices;
0145
0146 vector<FlavorHistory::FLAVOR_T> flavorSources;
0147
0148
0149 auto flavorHistoryEvent = std::make_unique<FlavorHistoryEvent>();
0150
0151
0152
0153
0154 vector<const Candidate *>::size_type j;
0155 vector<const Candidate *>::size_type j_max = particles.size();
0156 for (j = 0; j < j_max; ++j) {
0157
0158 const Candidate *p = particles[j];
0159
0160 vector<Candidate const *>::size_type partonIndex = j;
0161 vector<Candidate const *>::size_type progenitorIndex = j_max;
0162 bool foundProgenitor = false;
0163 FlavorHistory::FLAVOR_T flavorSource = FlavorHistory::FLAVOR_NULL;
0164
0165 int idabs = std::abs((p)->pdgId());
0166 int nDa = (p)->numberOfDaughters();
0167
0168
0169
0170 if (idabs <= FlavorHistory::tQuarkId && p->status() == 2 && idabs == pdgIdToSelect_) {
0171
0172 if (nDa > 0) {
0173
0174 if ((p)->pt() > ptMinShower_ && fabs((p)->eta()) < etaMaxShower_) {
0175 if (verbose_)
0176 cout << "--------------------------" << endl;
0177 if (verbose_)
0178 cout << "Processing particle " << j << ", status = " << p->status() << ", pdg id = " << p->pdgId() << endl;
0179
0180
0181 vector<Candidate const *> allParents;
0182 getAncestors(*p, allParents);
0183
0184 if (verbose_) {
0185 cout << "Parents : " << endl;
0186 vector<Candidate const *>::const_iterator iprint = allParents.begin(), iprintend = allParents.end();
0187 for (; iprint != iprintend; ++iprint)
0188 cout << " status = " << (*iprint)->status() << ", pdg id = " << (*iprint)->pdgId()
0189 << ", pt = " << (*iprint)->pt() << endl;
0190 }
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212 vector<Candidate const *>::size_type a_size = allParents.size();
0213
0214
0215
0216
0217 for (vector<Candidate const *>::size_type i = 0; i < a_size && !foundProgenitor; ++i) {
0218 const Candidate *aParent = allParents[i];
0219 vector<Candidate const *>::const_iterator found = find(particles.begin(), particles.end(), aParent);
0220
0221
0222 progenitorIndex = found - particles.begin();
0223
0224 int aParentId = std::abs(aParent->pdgId());
0225
0226
0227
0228
0229 if (progenitorIndex > 5) {
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239 if (aParent->numberOfMothers() == 2 && aParent->pdgId() == p->pdgId() && aParent->status() == 3) {
0240 if (verbose_)
0241 cout << "Matrix Element progenitor" << endl;
0242 if (flavorSource == FlavorHistory::FLAVOR_NULL)
0243 flavorSource = FlavorHistory::FLAVOR_ME;
0244 foundProgenitor = true;
0245 }
0246
0247 else if ((aParentId > pdgIdToSelect_ && aParentId < FlavorHistory::gluonId) ||
0248 aParentId > FlavorHistory::gluonId) {
0249 if (verbose_)
0250 cout << "Flavor decay progenitor" << endl;
0251 if (flavorSource == FlavorHistory::FLAVOR_NULL)
0252 flavorSource = FlavorHistory::FLAVOR_DECAY;
0253 foundProgenitor = true;
0254 }
0255 }
0256 }
0257
0258
0259
0260 if (!foundProgenitor) {
0261 if (flavorSource == FlavorHistory::FLAVOR_NULL)
0262 flavorSource = FlavorHistory::FLAVOR_GS;
0263
0264 for (vector<Candidate const *>::size_type i = 0; i < a_size && !foundProgenitor; ++i) {
0265 const Candidate *aParent = allParents[i];
0266 vector<Candidate const *>::const_iterator found = find(particles.begin(), particles.end(), aParent);
0267
0268 progenitorIndex = found - particles.begin();
0269
0270 if (aParent->numberOfMothers() == 1 && aParent->status() == 3) {
0271 foundProgenitor = true;
0272 }
0273 }
0274 }
0275
0276 }
0277 }
0278
0279
0280 if (!foundProgenitor)
0281 progenitorIndex = j_max;
0282
0283
0284
0285 partonIndices.push_back(partonIndex);
0286 progenitorIndices.push_back(progenitorIndex);
0287 flavorSources.push_back(flavorSource);
0288 sisterIndices.push_back(-1);
0289
0290 }
0291 }
0292
0293
0294
0295
0296
0297 if (verbose_)
0298 cout << "Making sisters" << endl;
0299
0300 if (partonIndices.size() == progenitorIndices.size() && !partonIndices.empty()) {
0301
0302 for (unsigned int ii = 0; ii < partonIndices.size(); ++ii) {
0303
0304 const Candidate *candi = particles[partonIndices[ii]];
0305 if (verbose_)
0306 cout << " Sister candidate particle 1: " << ii << ", pdgid = " << candi->pdgId()
0307 << ", status = " << candi->status() << endl;
0308
0309
0310
0311 for (unsigned int jj = 0; jj < partonIndices.size(); ++jj) {
0312 if (ii != jj) {
0313 const Candidate *candj = particles[partonIndices[jj]];
0314 if (verbose_)
0315 cout << " Sister candidate particle 2: " << jj << ", pdgid = " << candj->pdgId()
0316 << ", status = " << candj->status() << endl;
0317
0318
0319 if (candi->pdgId() == -1 * candj->pdgId() && candi->status() == candj->status()) {
0320 sisterIndices[ii] = partonIndices[jj];
0321 if (verbose_)
0322 cout << "Parton " << partonIndices[ii] << " has sister " << sisterIndices[ii] << endl;
0323 }
0324 }
0325 }
0326
0327
0328 if (sisterIndices[ii] < 0) {
0329 if (verbose_)
0330 cout << "No sister. Classified as flavor excitation" << endl;
0331 flavorSources[ii] = FlavorHistory::FLAVOR_EXC;
0332 }
0333
0334 if (verbose_)
0335 cout << "Getting matched jet" << endl;
0336
0337 CandidateView::const_iterator matched = getClosestJet(matchedView, *candi);
0338 CandidateView::const_iterator matchedBegin = matchedView->begin();
0339 CandidateView::const_iterator matchedEnd = matchedView->end();
0340
0341 if (verbose_)
0342 cout << "Getting sister jet" << endl;
0343
0344 CandidateView::const_iterator sister =
0345 (sisterIndices[ii] >= 0 && static_cast<unsigned int>(sisterIndices[ii]) < particles.size())
0346 ? getClosestJet(matchedView, *particles[sisterIndices[ii]])
0347 : matchedEnd;
0348
0349 if (verbose_)
0350 cout << "Making matched shallow clones : ";
0351 ShallowClonePtrCandidate matchedCand;
0352 if (matched != matchedEnd) {
0353 if (verbose_)
0354 cout << " found" << endl;
0355 matchedCand = ShallowClonePtrCandidate(CandidatePtr(matchedView, matched - matchedBegin));
0356 } else {
0357 if (verbose_)
0358 cout << " NOT found" << endl;
0359 }
0360
0361 if (verbose_)
0362 cout << "Making sister shallow clones : ";
0363 ShallowClonePtrCandidate sisterCand;
0364 if (sister != matchedEnd) {
0365 if (verbose_)
0366 cout << " found" << endl;
0367 sisterCand = ShallowClonePtrCandidate(CandidatePtr(matchedView, sister - matchedBegin));
0368 } else {
0369 if (verbose_)
0370 cout << " NOT found" << endl;
0371 }
0372
0373 if (verbose_)
0374 cout << "Making history object" << endl;
0375
0376 FlavorHistory history(flavorSources[ii],
0377 particlesViewH,
0378 partonIndices[ii],
0379 progenitorIndices[ii],
0380 sisterIndices[ii],
0381 matchedCand,
0382 sisterCand);
0383 if (verbose_)
0384 cout << "Adding flavor history : " << history << endl;
0385 flavorHistoryEvent->push_back(history);
0386 }
0387 }
0388
0389
0390 if (flavorHistoryEvent->size() > 0) {
0391
0392 if (verbose_)
0393 cout << "Caching flavor history event" << endl;
0394 flavorHistoryEvent->cache();
0395
0396 if (verbose_) {
0397 cout << "Outputting pdg id = " << pdgIdToSelect_ << " with nelements = " << flavorHistoryEvent->size() << endl;
0398 vector<FlavorHistory>::const_iterator i = flavorHistoryEvent->begin(), iend = flavorHistoryEvent->end();
0399 for (; i != iend; ++i) {
0400 cout << *i << endl;
0401 }
0402 }
0403 }
0404
0405
0406 evt.put(std::move(flavorHistoryEvent), flavorHistoryName_);
0407 }
0408
0409
0410 void FlavorHistoryProducer::getAncestors(const Candidate &c, vector<Candidate const *> &moms) const {
0411 if (c.numberOfMothers() == 1) {
0412 const Candidate *dau = &c;
0413 const Candidate *mom = c.mother();
0414 while (dau->numberOfMothers() != 0) {
0415 moms.push_back(dau);
0416 dau = mom;
0417 mom = dau->mother();
0418 }
0419 }
0420 }
0421
0422 CandidateView::const_iterator FlavorHistoryProducer::getClosestJet(Handle<CandidateView> const &pJets,
0423 reco::Candidate const &parton) const {
0424 double dr = matchDR_;
0425 CandidateView::const_iterator j = pJets->begin(), jend = pJets->end();
0426 CandidateView::const_iterator bestJet = pJets->end();
0427 for (; j != jend; ++j) {
0428 double dri = deltaR(parton.p4(), j->p4());
0429 if (dri < dr) {
0430 dr = dri;
0431 bestJet = j;
0432 }
0433 }
0434 return bestJet;
0435 }
0436
0437 #include "FWCore/Framework/interface/MakerMacros.h"
0438
0439 DEFINE_FWK_MODULE(FlavorHistoryProducer);