Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:17 

0001 /** class
0002  *
0003  * \author Stephen Mrenna, FNAL
0004  *
0005  * \version $Id: FlavorHistoryProducer.cc,v 1.0
0006  *
0007  */
0008 
0009 // -------------------------------------------------------------
0010 // Identify the ancestry of the Quark
0011 //
0012 //
0013 // Matrix Element:
0014 //    Status 3 parent with precisely 2 "grandparents" that
0015 //    is outside of the "initial" section (0-5) that has the
0016 //    same ID as the status 2 parton in question.
0017 //
0018 // Flavor excitation:
0019 //    Almost the same as the matrix element classification,
0020 //    but has only one outgoing parton product instead of two.
0021 //
0022 // Gluon splitting:
0023 //    Parent is a quark of a different flavor than the parton
0024 //    in question, or a gluon. Can come from either ISR or FSR.
0025 //
0026 // True decay:
0027 //    Decays from a resonance like top, Higgs, etc.
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   /// constructor
0057   FlavorHistoryProducer(const edm::ParameterSet &);
0058 
0059 private:
0060   /// process one event
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_;         // GenParticles source collection name
0069   const edm::EDGetTokenT<reco::CandidateView> matchedSrcToken_;  // matched particles source collection name
0070   const double matchDR_;                                         // delta r to match matched particles
0071   const int pdgIdToSelect_;                                      // pdg of hf partons to select
0072   const double ptMinParticle_;                                   // minimum pt of the partons
0073   const double ptMinShower_;                                     // minimum pt of the shower
0074   const double etaMaxParticle_;                                  // max eta of the parton
0075   const double etaMaxShower_;                                    // max eta of the shower
0076   const std::string flavorHistoryName_;                          // name to give flavor history
0077   const bool verbose_;                                           // verbose flag
0078 };
0079 
0080 // #include "DataFormats/Common/interface/ValueMap.h"
0081 // #include <iterators>
0082 
0083 // -------------------------------------------------------------
0084 // Identify the ancestry of the Quark
0085 //
0086 //
0087 // Matrix Element:
0088 //    Status 3 parent with precisely 2 "grandparents" that
0089 //    is outside of the "initial" section (0-5) that has the
0090 //    same ID as the status 2 parton in question.
0091 //
0092 // Flavor excitation:
0093 //    If we find only one outgoing parton.
0094 //
0095 // Gluon splitting:
0096 //    Parent is a quark of a different flavor than the parton
0097 //    in question, or a gluon.
0098 //    Can come from either ISR or FSR.
0099 //
0100 // True decay:
0101 //    Decays from a resonance like top, Higgs, etc.
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   // Get a handle to the particle collection (OwnVector)
0127   Handle<CandidateView> particlesViewH;
0128   evt.getByToken(srcToken_, particlesViewH);
0129 
0130   Handle<CandidateView> matchedView;
0131   evt.getByToken(matchedSrcToken_, matchedView);
0132 
0133   // Copy the View to an vector for easier iterator manipulation convenience
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   // List of indices for the partons to add
0140   vector<int> partonIndices;
0141   // List of progenitors for those partons
0142   vector<int> progenitorIndices;
0143   // List of sisters for those partons
0144   vector<int> sisterIndices;
0145   // Flavor sources
0146   vector<FlavorHistory::FLAVOR_T> flavorSources;
0147 
0148   // Make a new flavor history vector
0149   auto flavorHistoryEvent = std::make_unique<FlavorHistoryEvent>();
0150 
0151   // ------------------------------------------------------------
0152   // Loop over partons
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     // Get the candidate
0158     const Candidate *p = particles[j];
0159     // Set up indices that we'll need for the flavor history
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     // Check if we have a status 2 or 3 particle, which is a parton before the string.
0169     // Only consider quarks.
0170     if (idabs <= FlavorHistory::tQuarkId && p->status() == 2 && idabs == pdgIdToSelect_) {
0171       // Ensure the parton in question has daughters
0172       if (nDa > 0) {
0173         // Ensure the parton passes some minimum kinematic cuts
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           //Main (but simple) workhorse to get all ancestors
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           // Identify the ancestry of the Quark
0194           //
0195           //
0196           // Matrix Element:
0197           //    Status 3 parent with precisely 2 "grandparents" that
0198           //    is outside of the "initial" section (0-5) that has the
0199           //    same ID as the status 2 parton in question.
0200           //
0201           // Flavor excitation:
0202           //    If we find only one outgoing parton.
0203           //
0204           // Gluon splitting:
0205           //    Parent is a quark of a different flavor than the parton
0206           //    in question, or a gluon.
0207           //    Can come from either ISR or FSR.
0208           //
0209           // True decay:
0210           //    Decays from a resonance like top, Higgs, etc.
0211           // -------------------------------------------------------------
0212           vector<Candidate const *>::size_type a_size = allParents.size();
0213 
0214           //
0215           // Loop over all the ancestors of this parton and find the progenitor.
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             // Get the index of the progenitor candidate
0222             progenitorIndex = found - particles.begin();
0223 
0224             int aParentId = std::abs(aParent->pdgId());
0225 
0226             // -----------------------------------------------------------------------
0227             // Here we examine particles that were produced after the collision
0228             // -----------------------------------------------------------------------
0229             if (progenitorIndex > 5) {
0230               // Here is where we have a matrix element process.
0231               // The signature is to have a status 3 parent with precisely 2 parents
0232               // that is outside of the "initial" section (0-5) that has the
0233               // same ID as the status 2 parton in question.
0234               // NOTE: If this parton has no sister, then this will be classified as
0235               // a flavor excitation. Often the only difference is that the matrix element
0236               // cases have a sister, while the flavor excitation cases do not.
0237               // If we do not find a sister below, this will be classified as a flavor
0238               // excitation.
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               // Here we have a true decay. Parent is not a quark or a gluon.
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             }  // end if progenitorIndex > 5
0256           }  // end loop over parents
0257 
0258           // Otherwise, classify it as gluon splitting. Take the first status 3 parton with 1 parent
0259           // that you see as the progenitor
0260           if (!foundProgenitor) {
0261             if (flavorSource == FlavorHistory::FLAVOR_NULL)
0262               flavorSource = FlavorHistory::FLAVOR_GS;
0263             // Now get the parton with only one parent (the proton) and that is the progenitor
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               // Get the index of the progenitor candidate
0268               progenitorIndex = found - particles.begin();
0269 
0270               if (aParent->numberOfMothers() == 1 && aParent->status() == 3) {
0271                 foundProgenitor = true;
0272               }
0273             }  // end loop over parents
0274           }  // end if we haven't found a progenitor, and if we haven't found a status 3 parton of the same flavor
0275           // (i.e. end if examining gluon splitting)
0276         }  // End if this parton passes some minimal kinematic cuts
0277       }  // End if this particle is status 2 (has strings as daughters)
0278 
0279       // Make sure we've actually found a progenitor
0280       if (!foundProgenitor)
0281         progenitorIndex = j_max;
0282 
0283       // We've found the particle, add to the list
0284 
0285       partonIndices.push_back(partonIndex);
0286       progenitorIndices.push_back(progenitorIndex);
0287       flavorSources.push_back(flavorSource);
0288       sisterIndices.push_back(-1);  // set below
0289 
0290     }  // End if this particle was a status==2 parton
0291   }  // end loop over particles
0292 
0293   // Now add sisters.
0294   // Also if the event is preliminarily classified as "matrix element", check to
0295   // make sure that they have a sister. If not, it is flavor excitation.
0296 
0297   if (verbose_)
0298     cout << "Making sisters" << endl;
0299   // First make sure nothing went terribly wrong:
0300   if (partonIndices.size() == progenitorIndices.size() && !partonIndices.empty()) {
0301     // Now loop over the candidates
0302     for (unsigned int ii = 0; ii < partonIndices.size(); ++ii) {
0303       // Get the iith particle
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       // Get the iith progenitor
0309       // Now loop over the other flavor history candidates and
0310       // attempt to find a sister to this one
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           // They should be opposite in pdgid and have the same status, and
0318           // the same progenitory.
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       // Here, ensure that there is a sister. Otherwise this is "flavor excitation"
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       // Get the closest match in the matched object collection
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       // Get the closest sister in the matched object collection, if sister is found
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       // Now create the flavor history object
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   // If we got any partons, cache them and then write them out
0390   if (flavorHistoryEvent->size() > 0) {
0391     // Calculate some nice variables for FlavorHistoryEvent
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   // Now add the flavor history to the event record
0406   evt.put(std::move(flavorHistoryEvent), flavorHistoryName_);
0407 }
0408 
0409 // Helper function to get all ancestors of this candidate
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!