HepMCCandidate/src/FlavorHistoryEvent.cc

0001 #include "DataFormats/HepMCCandidate/interface/FlavorHistoryEvent.h"
0002 #include "TLorentzVector.h"
0003
0004 #include <iostream>
0005
0006 using namespace reco;
0007 using namespace std;
0008
0009 // Loop over flavor histories, count number of genjets with
0010 // flavor b, c, or l
0011 void FlavorHistoryEvent::cache() {
0012   bool verbose = false;
0013
0014   if (verbose)
0015     cout << "----- Caching Flavor History Event  -----" << endl;
0016   // Set cached to false
0017   cached_ = false;
0018   nb_ = nc_ = 0;
0019   highestFlavor_ = 0;
0020   dR_ = 0.0;
0021   flavorSource_ = FlavorHistory::FLAVOR_NULL;
0022
0023   // get list of flavor --> type --> dR.
0024   // Will sort this later to determine event classification
0025   vector<helpers::FlavorHistoryEventHelper> classification;
0026
0027   // get iterators to the history vector
0028   const_iterator i = histories_.begin(), ibegin = histories_.begin(), iend = histories_.end();
0029   // loop over the history vector and count the number of
0030   // partons of flavors "b" and "c" that have a matched genjet.
0031   for (; i != iend; ++i) {
0032     FlavorHistory const& flavHist = *i;
0033     if (verbose)
0034       cout << "  Processing flavor history: " << i - ibegin << " = " << endl << flavHist << endl;
0035     CandidatePtr const& parton = flavHist.parton();
0036     flavor_type flavorSource = flavHist.flavorSource();
0037
0038     // Now examine the matched jets to see what the classification should be.
0039     int pdgId = -1;
0040     if (parton.isNonnull())
0041       pdgId = abs(parton->pdgId());
0042     ShallowClonePtrCandidate const& matchedJet = flavHist.matchedJet();
0043     // Only count events with a matched genjet
0044     if (matchedJet.masterClonePtr().isNonnull()) {
0045       TLorentzVector p41(matchedJet.px(), matchedJet.py(), matchedJet.pz(), matchedJet.energy());
0046       if (pdgId == 5)
0047         nb_++;
0048       if (pdgId == 4)
0049         nc_++;
0050
0051       // Get the sister genjet
0052       ShallowClonePtrCandidate const& sisterJet = i->sisterJet();
0053       TLorentzVector p42(sisterJet.px(), sisterJet.py(), sisterJet.pz(), sisterJet.energy());
0054
0055       // Now check the source.
0056       double dR = -1;
0057       if (sisterJet.masterClonePtr().isNonnull()) {
0058         dR = p41.DeltaR(p42);
0059       }
0060       // Add to the vector to be sorted later
0061       if (verbose)
0062         cout << "Adding classification: pdgId = " << pdgId << ", flavorSource = " << flavorSource << ", dR = " << dR
0063              << endl;
0064       classification.push_back(helpers::FlavorHistoryEventHelper(pdgId, flavorSource, dR));
0065     } else {
0066       if (verbose)
0067         cout << "No matched jet found, not adding to classification list" << endl;
0068     }
0069   }
0070
0071   // Sort by priority
0072
0073   // Priority is:
0074   //
0075   //  1. flavor (5 > 4)
0076   //  2. type:
0077   //      2a. Flavor decay
0078   //      2b. Matrix element
0079   //      2c. Flavor excitation
0080   //      2d. Gluon splitting
0081   //  3. delta R (if applicable)
0082   if (!classification.empty()) {
0083     std::sort(classification.begin(), classification.end());
0084
0085     if (verbose) {
0086       cout << "Writing out list of classifications" << endl;
0087       copy(classification.begin(), classification.end(), ostream_iterator<helpers::FlavorHistoryEventHelper>(cout, ""));
0088     }
0089
0090     helpers::FlavorHistoryEventHelper const& best = *(classification.rbegin());
0091     dR_ = best.dR;
0092     highestFlavor_ = best.flavor;
0093     flavorSource_ = best.flavorSource;
0094   } else {
0095     dR_ = -1.0;
0096     highestFlavor_ = 0;
0097     flavorSource_ = FlavorHistory::FLAVOR_NULL;
0098   }
0099
0100   // now we're cached, can return values quickly
0101   cached_ = true;
0102 }