Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:14 

0001 #include <memory>
0002 #include <string>
0003 #include <vector>
0004 #include <map>
0005 
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/stream/EDProducer.h"
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 
0010 #include "FWCore/Common/interface/Provenance.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/Utilities/interface/EDMException.h"
0013 #include "FWCore/Utilities/interface/InputTag.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0016 
0017 #include "AnalysisDataFormats/TopObjects/interface/TopGenEvent.h"
0018 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0019 
0020 /**
0021    \class   TopDecaySubset TopDecaySubset.h "TopQuarkAnalysis/TopEventProducers/interface/TopDecaySubset.h"
0022 
0023    \brief   Module to produce the subset of generator particles directly contained in top decay chains
0024 
0025    The module produces the subset of generator particles directly contained in top decay chains. The
0026    particles are saved as a collection of reco::GenParticles. Depending on the configuration of the module,
0027    the 4-vector kinematics can be taken from the status-3 particles (ME before parton showering) or from
0028    the status-2 particles (after parton showering), additionally radiated gluons may be considered during
0029    the creation of the subset or not.
0030 */
0031 
0032 class TopDecaySubset : public edm::stream::EDProducer<> {
0033 public:
0034   /// supported modes to fill the new vectors
0035   /// of gen particles
0036   enum FillMode { kStable, kME };
0037   /// classification of potential shower types
0038   enum ShowerModel { kStart = -1, kNone, kPythia, kHerwig, kPythia8, kSherpa };
0039   /// supported modes to run the code
0040   enum RunMode { kRun1, kRun2 };
0041 
0042   /// default constructor
0043   explicit TopDecaySubset(const edm::ParameterSet& cfg);
0044   /// write output into the event
0045   void produce(edm::Event& event, const edm::EventSetup& setup) override;
0046 
0047 private:
0048   /// find top quarks in list of input particles
0049   std::vector<const reco::GenParticle*> findTops(const reco::GenParticleCollection& parts);
0050   /// find primal top quarks (top quarks from the hard interaction)
0051   /// for Pythia6 this is identical to findDecayingTops
0052   std::vector<const reco::GenParticle*> findPrimalTops(const reco::GenParticleCollection& parts);
0053   /// find decaying top quarks (quarks that decay to qW)
0054   /// for Pythia6 this is identical to findPrimalTops
0055   std::vector<const reco::GenParticle*> findDecayingTops(const reco::GenParticleCollection& parts);
0056   /// find W bosons that come from top quark decays
0057   /// for Pythia6 this is identical to findDecayingW
0058   const reco::GenParticle* findPrimalW(const reco::GenParticle* top);
0059   /// find W bosons that come from top quark decays and decay themselves (end of the MC chain)
0060   /// for Pythia6 this is identical to findPrimalW
0061   //  const reco::GenParticle* findDecayingW(const reco::GenParticle* top);
0062   /// find the last particle in a (potentially) long chain of state transitions
0063   /// e.g. top[status==22]-> top[status==44 -> top[status==44] ->
0064   /// top[status==44] -> top[status==62]
0065   /// this function would pick the top with status 62
0066   const reco::GenParticle* findLastParticleInChain(const reco::GenParticle* p);
0067   /// check the decay chain for the used shower model
0068   ShowerModel checkShowerModel(const std::vector<const reco::GenParticle*>& tops) const;
0069   ShowerModel checkShowerModel(edm::Event& event);
0070   /// check whether W bosons are contained in the original gen particle listing
0071   void checkWBosons(std::vector<const reco::GenParticle*>& tops) const;
0072   /// fill output vector for full decay chain
0073   void fillListing(const std::vector<const reco::GenParticle*>& tops, reco::GenParticleCollection& target);
0074   /// fill output vector for full decay chain
0075   void fillListing(const std::vector<const reco::GenParticle*>& primalTops,
0076                    const std::vector<const reco::GenParticle*>& decayingTops,
0077                    reco::GenParticleCollection& target);
0078 
0079   /// clear references
0080   void clearReferences();
0081   /// fill references for output vector
0082   void fillReferences(const reco::GenParticleRefProd& refProd, reco::GenParticleCollection& target);
0083   /// calculate lorentz vector from input (dedicated to top reconstruction)
0084   reco::Particle::LorentzVector p4(const std::vector<const reco::GenParticle*>::const_iterator top, int statusFlag);
0085   /// calculate lorentz vector from input
0086   reco::Particle::LorentzVector p4(const reco::GenParticle::const_iterator part, int statusFlag);
0087   /// recursively fill vector for all further decay particles of a given particle
0088   void addDaughters(int& idx,
0089                     const reco::GenParticle::const_iterator part,
0090                     reco::GenParticleCollection& target,
0091                     bool recursive = true);
0092   /// fill vector including all radiations from quarks originating from W/top
0093   void addRadiation(int& idx, const reco::GenParticle::const_iterator part, reco::GenParticleCollection& target);
0094   void addRadiation(int& idx, const reco::GenParticle* part, reco::GenParticleCollection& target);
0095 
0096 private:
0097   /// input tag for the genParticle source
0098   edm::EDGetTokenT<reco::GenParticleCollection> srcToken_;
0099   /// input tag for the genEventInfo source
0100   edm::EDGetTokenT<GenEventInfoProduct> genEventInfo_srcToken_;
0101   /// add radiation or not?
0102   bool addRadiation_;
0103   /// print the whole list of input particles or not?
0104   /// mode of decaySubset creation
0105   FillMode fillMode_;
0106   /// parton shower mode (filled in checkShowerModel)
0107   ShowerModel showerModel_;
0108   /// run mode (Run1 || Run2)
0109   RunMode runMode_;
0110 
0111   /// index in new evt listing of parts with daughters;
0112   /// has to be set to -1 in produce to deliver consistent
0113   /// results!
0114   int motherPartIdx_;
0115   /// management of daughter indices for fillRefs
0116   std::map<int, std::vector<int> > refs_;
0117 };
0118 
0119 /// default constructor
0120 TopDecaySubset::TopDecaySubset(const edm::ParameterSet& cfg)
0121     : srcToken_(consumes<reco::GenParticleCollection>(cfg.getParameter<edm::InputTag>("src"))),
0122       genEventInfo_srcToken_(mayConsume<GenEventInfoProduct>(edm::InputTag("generator"))),
0123       addRadiation_(cfg.getParameter<bool>("addRadiation")),
0124       showerModel_(kStart),
0125       runMode_(kRun1) {
0126   // mapping of the corresponding fillMode; see FillMode
0127   // enumerator of TopDecaySubset for available modes
0128   std::string mode = cfg.getParameter<std::string>("fillMode");
0129   if (mode == "kME")
0130     fillMode_ = kME;
0131   else if (mode == "kStable")
0132     fillMode_ = kStable;
0133   else
0134     throw cms::Exception("Configuration") << mode << " is not a supported FillMode!\n";
0135 
0136   mode = cfg.getParameter<std::string>("runMode");
0137   if (mode == "Run1")
0138     runMode_ = kRun1;
0139   else if (mode == "Run2")
0140     runMode_ = kRun2;
0141   else
0142     throw cms::Exception("Configuration") << mode << " is not a supported RunMode!\n";
0143 
0144   // produces a set of GenParticles following
0145   // the decay branch of top quarks to the first level of
0146   // stable decay products
0147   produces<reco::GenParticleCollection>();
0148 }
0149 
0150 /// write output into the event
0151 void TopDecaySubset::produce(edm::Event& event, const edm::EventSetup& setup) {
0152   // create target vector
0153   std::unique_ptr<reco::GenParticleCollection> target(new reco::GenParticleCollection);
0154 
0155   // get source collection
0156   edm::Handle<reco::GenParticleCollection> src;
0157   event.getByToken(srcToken_, src);
0158 
0159   // find out what generator we are dealing with
0160   if (showerModel_ == kStart && runMode_ == kRun2) {
0161     showerModel_ = checkShowerModel(event);
0162   }
0163 
0164   // find top quarks in list of input particles
0165   std::vector<const reco::GenParticle*> tops;
0166   if (runMode_ == kRun1)
0167     tops = findTops(*src);
0168   else
0169     tops = findPrimalTops(*src);
0170 
0171   // determine shower model (only in first event)
0172   if (showerModel_ == kStart && runMode_ == kRun1)
0173     showerModel_ = checkShowerModel(tops);
0174 
0175   if (showerModel_ != kNone) {
0176     // check sanity of W bosons
0177     if (runMode_ == kRun1)
0178       checkWBosons(tops);
0179     else {
0180       // nothing for the moment
0181     }
0182 
0183     // get ref product from the event
0184     const reco::GenParticleRefProd ref = event.getRefBeforePut<reco::GenParticleCollection>();
0185     // clear existing refs
0186     clearReferences();
0187     if (runMode_ == kRun1) {
0188       // fill output
0189       fillListing(tops, *target);
0190       // fill references
0191       fillReferences(ref, *target);
0192     } else {
0193       std::vector<const reco::GenParticle*> decaying_tops = findDecayingTops(*src);
0194       // fill output
0195       fillListing(tops, decaying_tops, *target);
0196       // fill references
0197       fillReferences(ref, *target);
0198     }
0199   }
0200 
0201   // write vectors to the event
0202   event.put(std::move(target));
0203 }
0204 
0205 /// find top quarks in list of input particles
0206 std::vector<const reco::GenParticle*> TopDecaySubset::findTops(const reco::GenParticleCollection& parts) {
0207   std::vector<const reco::GenParticle*> tops;
0208   for (reco::GenParticleCollection::const_iterator t = parts.begin(); t != parts.end(); ++t) {
0209     if (std::abs(t->pdgId()) == TopDecayID::tID && t->status() == TopDecayID::unfrag)
0210       tops.push_back(&(*t));
0211   }
0212   return tops;
0213 }
0214 
0215 /// find primal top quarks (top quarks from the hard interaction)
0216 /// for Pythia6 this is identical to findDecayingTops
0217 std::vector<const reco::GenParticle*> TopDecaySubset::findPrimalTops(const reco::GenParticleCollection& parts) {
0218   std::vector<const reco::GenParticle*> tops;
0219   for (reco::GenParticleCollection::const_iterator t = parts.begin(); t != parts.end(); ++t) {
0220     if (std::abs(t->pdgId()) != TopDecayID::tID)
0221       continue;
0222 
0223     bool hasTopMother = false;
0224     for (unsigned idx = 0; idx < t->numberOfMothers(); ++idx) {
0225       if (std::abs(t->mother(idx)->pdgId()) == TopDecayID::tID)
0226         hasTopMother = true;
0227     }
0228 
0229     if (hasTopMother)  // not a primal top
0230       continue;
0231     tops.push_back(&(*t));
0232   }
0233 
0234   return tops;
0235 }
0236 
0237 /// find decaying top quarks (quarks that decay to qW)
0238 /// for Pythia6 this is identical to findPrimalTops
0239 std::vector<const reco::GenParticle*> TopDecaySubset::findDecayingTops(const reco::GenParticleCollection& parts) {
0240   std::vector<const reco::GenParticle*> tops;
0241   for (reco::GenParticleCollection::const_iterator t = parts.begin(); t != parts.end(); ++t) {
0242     if (std::abs(t->pdgId()) != TopDecayID::tID)
0243       continue;
0244 
0245     bool hasTopDaughter = false;
0246     for (unsigned idx = 0; idx < t->numberOfDaughters(); ++idx) {
0247       if (std::abs(t->daughter(idx)->pdgId()) == TopDecayID::tID)
0248         hasTopDaughter = true;
0249     }
0250 
0251     if (hasTopDaughter)  // not a decaying top
0252       continue;
0253     tops.push_back(&(*t));
0254   }
0255 
0256   return tops;
0257 }
0258 
0259 /// find W bosons that come from top quark decays
0260 /// for Pythia6 this is identical to findDecayingW
0261 const reco::GenParticle* TopDecaySubset::findPrimalW(const reco::GenParticle* top) {
0262   unsigned int w_index = 0;
0263   for (unsigned idx = 0; idx < top->numberOfDaughters(); ++idx) {
0264     if (std::abs(top->daughter(idx)->pdgId()) == TopDecayID::WID) {
0265       w_index = idx;
0266       break;
0267     }
0268   }
0269   return static_cast<const reco::GenParticle*>(top->daughter(w_index));
0270 }
0271 
0272 /// find W bosons that come from top quark decays and decay themselves (end of the MC chain)
0273 /// for Pythia6 this is identical to findPrimalW
0274 //const reco::GenParticle* TopDecaySubset::findDecayingW(
0275 //      const reco::GenParticle* top) {
0276 //  const reco::GenParticle* decaying_W = findLastParticleInChain(findPrimalW(top));
0277 //  return findLastParticleInChain(findPrimalW(top));
0278 //}
0279 
0280 /// find the last particle in a (potentially) long chain of state transitions
0281 /// e.g. top[status==22]-> top[status==44 -> top[status==44] ->
0282 /// top[status==44] -> top[status==62]
0283 /// this function would pick the top with status 62
0284 const reco::GenParticle* TopDecaySubset::findLastParticleInChain(const reco::GenParticle* p) {
0285   int particleID = std::abs(p->pdgId());
0286   bool containsItself = false;
0287   unsigned int d_idx = 0;
0288   for (unsigned idx = 0; idx < p->numberOfDaughters(); ++idx) {
0289     if (std::abs(p->daughter(idx)->pdgId()) == particleID) {
0290       containsItself = true;
0291       d_idx = idx;
0292     }
0293   }
0294 
0295   if (!containsItself)
0296     return p;
0297   else {
0298     if (showerModel_ == kPythia) {
0299       // Pythia6 has a weird idea of W bosons (and maybe other particles)
0300       // W (status == 3) -> q qbar' W. The new W is status 2 and has no daughters
0301       if (p->status() == 3)
0302         return p;
0303     }
0304     return findLastParticleInChain(static_cast<const reco::GenParticle*>(p->daughter(d_idx)));
0305   }
0306 }
0307 
0308 /// check the decay chain for the exploited shower model
0309 TopDecaySubset::ShowerModel TopDecaySubset::checkShowerModel(const std::vector<const reco::GenParticle*>& tops) const {
0310   for (std::vector<const reco::GenParticle*>::const_iterator it = tops.begin(); it != tops.end(); ++it) {
0311     const reco::GenParticle* top = *it;
0312     // check for kHerwig type showers: here the status 3 top quark will
0313     // have a single status 2 top quark as daughter, which has again 3
0314     // or more status 2 daughters
0315     if (top->numberOfDaughters() == 1) {
0316       if (top->begin()->pdgId() == top->pdgId() && top->begin()->status() == TopDecayID::stable &&
0317           top->begin()->numberOfDaughters() > 1)
0318         return kHerwig;
0319     }
0320     // check for kPythia type showers: here the status 3 top quark will
0321     // have all decay products and a status 2 top quark as daughters
0322     // the status 2 top quark will be w/o further daughters
0323     if (top->numberOfDaughters() > 1) {
0324       bool containsWBoson = false, containsQuarkDaughter = false;
0325       for (reco::GenParticle::const_iterator td = top->begin(); td != top->end(); ++td) {
0326         if (std::abs(td->pdgId()) < TopDecayID::tID)
0327           containsQuarkDaughter = true;
0328         if (std::abs(td->pdgId()) == TopDecayID::WID)
0329           containsWBoson = true;
0330       }
0331       if (containsQuarkDaughter && containsWBoson)
0332         return kPythia;
0333     }
0334   }
0335   // if neither Herwig nor Pythia like
0336   if (tops.empty())
0337     edm::LogInfo("decayChain") << " Failed to find top quarks in decay chain. Will assume that this a \n"
0338                                << " non-top sample and produce an empty decaySubset.\n";
0339   else
0340     throw edm::Exception(edm::errors::LogicError,
0341                          " Can not find back any of the supported hadronization models. Models \n"
0342                          " which are supported are:                                            \n"
0343                          " Pythia  LO(+PS): Top(status 3) --> WBoson(status 3), Quark(status 3)\n"
0344                          " Herwig NLO(+PS): Top(status 2) --> Top(status 3) --> Top(status 2)  \n");
0345   return kNone;
0346 }
0347 
0348 /// check the embedded MC information for the shower model
0349 TopDecaySubset::ShowerModel TopDecaySubset::checkShowerModel(edm::Event& event) {
0350   edm::Handle<GenEventInfoProduct> genEvtInfoProduct;
0351   event.getByToken(genEventInfo_srcToken_, genEvtInfoProduct);
0352 
0353   std::string moduleName = "";
0354   if (genEvtInfoProduct.isValid()) {
0355     const edm::StableProvenance& prov = event.getStableProvenance(genEvtInfoProduct.id());
0356     moduleName = edm::moduleName(prov, event.processHistory());
0357     if (moduleName == "ExternalGeneratorFilter") {
0358       moduleName = edm::parameterSet(prov, event.processHistory()).getParameter<std::string>("@external_type");
0359       edm::LogInfo("SpecialModule") << "GEN events are produced by ExternalGeneratorFilter, "
0360                                     << "which is a wrapper of the original module: " << moduleName;
0361     }
0362   }
0363 
0364   ShowerModel shower(kStart);
0365   if (moduleName.find("Pythia6") != std::string::npos)
0366     shower = kPythia;
0367   else if (moduleName.find("Pythia8") != std::string::npos)
0368     shower = kPythia8;
0369   else if (moduleName.find("Herwig6") != std::string::npos)
0370     shower = kHerwig;
0371   else if (moduleName.find("ThePEG") != std::string::npos)
0372     // Herwig++
0373     shower = kHerwig;
0374   else if (moduleName.find("Sherpa") != std::string::npos)
0375     shower = kSherpa;
0376   else
0377     shower = kNone;
0378   return shower;
0379 }
0380 /// check whether the W boson is contained in the original gen particle listing
0381 void TopDecaySubset::checkWBosons(std::vector<const reco::GenParticle*>& tops) const {
0382   unsigned nTops = tops.size();
0383   for (std::vector<const reco::GenParticle*>::iterator it = tops.begin(); it != tops.end();) {
0384     const reco::GenParticle* top = *it;
0385     bool isContained = false;
0386     bool expectedStatus = false;
0387     if (showerModel_ != kPythia && top->begin() == top->end())
0388       throw edm::Exception(edm::errors::LogicError, "showerModel_!=kPythia && top->begin()==top->end()\n");
0389     for (reco::GenParticle::const_iterator td = ((showerModel_ == kPythia) ? top->begin() : top->begin()->begin());
0390          td != ((showerModel_ == kPythia) ? top->end() : top->begin()->end());
0391          ++td) {
0392       if (std::abs(td->pdgId()) == TopDecayID::WID) {
0393         isContained = true;
0394         if (((showerModel_ == kPythia) ? td->status() == TopDecayID::unfrag : td->status() == TopDecayID::stable)) {
0395           expectedStatus = true;
0396           break;
0397         }
0398       }
0399     }
0400     if (!expectedStatus) {
0401       it = tops.erase(it);
0402       if (isContained)
0403         edm::LogInfo("decayChain") << " W boson does not have the expected status. This happens, e.g.,      \n"
0404                                    << " with MC@NLO in the case of additional ttbar pairs from radiated     \n"
0405                                    << " gluons. We hope everything is fine, remove the correspdonding top   \n"
0406                                    << " quark from our list since it is not part of the primary ttbar pair  \n"
0407                                    << " and try to continue.                                                \n";
0408     } else
0409       it++;
0410   }
0411   if (tops.empty() && nTops != 0)
0412     throw edm::Exception(edm::errors::LogicError,
0413                          " Did not find a W boson with appropriate status for any of the top   \n"
0414                          " quarks in this event. This means that the hadronization of the W    \n"
0415                          " boson might be screwed up or there is another problem with the      \n"
0416                          " particle listing in this MC sample. Please contact an expert.       \n");
0417 }
0418 
0419 /// fill output vector for full decay chain
0420 void TopDecaySubset::fillListing(const std::vector<const reco::GenParticle*>& tops,
0421                                  reco::GenParticleCollection& target) {
0422   unsigned int statusFlag;
0423   // determine status flag of the new
0424   // particle depending on the FillMode
0425   fillMode_ == kME ? statusFlag = 3 : statusFlag = 2;
0426 
0427   for (std::vector<const reco::GenParticle*>::const_iterator it = tops.begin(); it != tops.end(); ++it) {
0428     const reco::GenParticle* t = *it;
0429     // if particle is top or anti-top
0430     std::unique_ptr<reco::GenParticle> topPtr(
0431         new reco::GenParticle(t->threeCharge(), p4(it, statusFlag), t->vertex(), t->pdgId(), statusFlag, false));
0432     target.push_back(*topPtr);
0433     ++motherPartIdx_;
0434     // keep the top index for the map to manage the daughter refs
0435     int iTop = motherPartIdx_;
0436     std::vector<int> topDaughters;
0437     // define the W boson index (to be set later) for the map to
0438     // manage the daughter refs
0439     int iW = 0;
0440     std::vector<int> wDaughters;
0441     // sanity check
0442     if (showerModel_ != kPythia && t->begin() == t->end())
0443       throw edm::Exception(edm::errors::LogicError, "showerModel_!=kPythia && t->begin()==t->end()\n");
0444     //iterate over top daughters
0445     for (reco::GenParticle::const_iterator td = ((showerModel_ == kPythia) ? t->begin() : t->begin()->begin());
0446          td != ((showerModel_ == kPythia) ? t->end() : t->begin()->end());
0447          ++td) {
0448       if (td->status() == TopDecayID::unfrag && std::abs(td->pdgId()) <= TopDecayID::bID) {
0449         // if particle is beauty or other quark
0450         std::unique_ptr<reco::GenParticle> bPtr(
0451             new reco::GenParticle(td->threeCharge(), p4(td, statusFlag), td->vertex(), td->pdgId(), statusFlag, false));
0452         target.push_back(*bPtr);
0453         // increment & push index of the top daughter
0454         topDaughters.push_back(++motherPartIdx_);
0455         if (addRadiation_) {
0456           addRadiation(motherPartIdx_, td, target);
0457         }
0458       }
0459       // sanity check
0460       if (showerModel_ != kPythia && td->begin() == td->end())
0461         throw edm::Exception(edm::errors::LogicError, "showerModel_!=kPythia && td->begin()==td->end()\n");
0462       reco::GenParticle::const_iterator buffer = (showerModel_ == kPythia) ? td : td->begin();
0463       if (buffer->status() == TopDecayID::unfrag && std::abs(buffer->pdgId()) == TopDecayID::WID) {
0464         // if particle is a W boson
0465         std::unique_ptr<reco::GenParticle> wPtr(new reco::GenParticle(
0466             buffer->threeCharge(), p4(buffer, statusFlag), buffer->vertex(), buffer->pdgId(), statusFlag, true));
0467         target.push_back(*wPtr);
0468         // increment & push index of the top daughter
0469         topDaughters.push_back(++motherPartIdx_);
0470         // keep the W idx for the map
0471         iW = motherPartIdx_;
0472         if (addRadiation_) {
0473           addRadiation(motherPartIdx_, buffer, target);
0474         }
0475         if (showerModel_ != kPythia && buffer->begin() == buffer->end())
0476           throw edm::Exception(edm::errors::LogicError, "showerModel_!=kPythia && buffer->begin()==buffer->end()\n");
0477         // iterate over W daughters
0478         for (reco::GenParticle::const_iterator wd =
0479                  ((showerModel_ == kPythia) ? buffer->begin() : buffer->begin()->begin());
0480              wd != ((showerModel_ == kPythia) ? buffer->end() : buffer->begin()->end());
0481              ++wd) {
0482           // make sure the W daughter is of status unfrag and not the W itself
0483           if (wd->status() == TopDecayID::unfrag && !(std::abs(wd->pdgId()) == TopDecayID::WID)) {
0484             std::unique_ptr<reco::GenParticle> qPtr(new reco::GenParticle(
0485                 wd->threeCharge(), p4(wd, statusFlag), wd->vertex(), wd->pdgId(), statusFlag, false));
0486             target.push_back(*qPtr);
0487             // increment & push index of the top daughter
0488             wDaughters.push_back(++motherPartIdx_);
0489             if (wd->status() == TopDecayID::unfrag && std::abs(wd->pdgId()) == TopDecayID::tauID) {
0490               // add tau daughters if the particle is a tau pass
0491               // the daughter of the tau which is of status 2
0492               //addDaughters(motherPartIdx_, wd->begin(), target);
0493               // add tau daughters if the particle is a tau pass
0494               // the tau itself, which may add a tau daughter of
0495               // of status 2 to the listing
0496               addDaughters(motherPartIdx_, wd, target);
0497             }
0498           }
0499         }
0500       }
0501       if (addRadiation_ && buffer->status() == TopDecayID::stable &&
0502           (buffer->pdgId() == TopDecayID::glueID || std::abs(buffer->pdgId()) < TopDecayID::bID)) {
0503         // collect additional radiation from the top
0504         std::unique_ptr<reco::GenParticle> radPtr(new reco::GenParticle(
0505             buffer->threeCharge(), buffer->p4(), buffer->vertex(), buffer->pdgId(), statusFlag, false));
0506         target.push_back(*radPtr);
0507         // increment & push index of the top daughter
0508         topDaughters.push_back(++motherPartIdx_);
0509       }
0510     }
0511     // add potential sisters of the top quark;
0512     // only for top to prevent double counting
0513     if (t->numberOfMothers() > 0 && t->pdgId() == TopDecayID::tID) {
0514       for (reco::GenParticle::const_iterator ts = t->mother()->begin(); ts != t->mother()->end(); ++ts) {
0515         // loop over all daughters of the top mother i.e.
0516         // the two top quarks and their potential sisters
0517         if (std::abs(ts->pdgId()) != t->pdgId() && ts->pdgId() != t->mother()->pdgId()) {
0518           // add all further particles but the two top quarks and potential
0519           // cases where the mother of the top has itself as daughter
0520           reco::GenParticle* cand =
0521               new reco::GenParticle(ts->threeCharge(), ts->p4(), ts->vertex(), ts->pdgId(), ts->status(), false);
0522           std::unique_ptr<reco::GenParticle> sPtr(cand);
0523           target.push_back(*sPtr);
0524           if (ts->begin() != ts->end()) {
0525             // in case the sister has daughters increment
0526             // and add the first generation of daughters
0527             addDaughters(++motherPartIdx_, ts->begin(), target, false);
0528           }
0529         }
0530       }
0531     }
0532     // fill the map for the administration
0533     // of daughter indices
0534     refs_[iTop] = topDaughters;
0535     refs_[iW] = wDaughters;
0536   }
0537 }
0538 
0539 void TopDecaySubset::fillListing(const std::vector<const reco::GenParticle*>& primalTops,
0540                                  const std::vector<const reco::GenParticle*>& decayingTops,
0541                                  reco::GenParticleCollection& target) {
0542   std::vector<const reco::GenParticle*>::const_iterator top_start;
0543   std::vector<const reco::GenParticle*>::const_iterator top_end;
0544   if (fillMode_ == kME) {
0545     top_start = primalTops.begin();
0546     top_end = primalTops.end();
0547   } else {
0548     top_start = decayingTops.begin();
0549     top_end = decayingTops.end();
0550   }
0551   for (std::vector<const reco::GenParticle*>::const_iterator it = top_start; it != top_end; ++it) {
0552     const reco::GenParticle* t = *it;
0553     // summation might happen here
0554     std::unique_ptr<reco::GenParticle> topPtr(
0555         new reco::GenParticle(t->threeCharge(), t->p4(), t->vertex(), t->pdgId(), t->status(), false));
0556     target.push_back(*topPtr);
0557     ++motherPartIdx_;
0558 
0559     // keep the top index for the map to manage the daughter refs
0560     int iTop = motherPartIdx_;
0561     std::vector<int> topDaughters;
0562     // define the W boson index (to be set later) for the map to
0563     // manage the daughter refs
0564     int iW = 0;
0565     std::vector<int> wDaughters;
0566     const reco::GenParticle* final_top = findLastParticleInChain(t);
0567 
0568     //iterate over top daughters
0569     for (reco::GenParticle::const_iterator td = final_top->begin(); td != final_top->end(); ++td) {
0570       if (std::abs(td->pdgId()) <= TopDecayID::bID) {
0571         // if particle is beauty or other quark
0572         std::unique_ptr<reco::GenParticle> qPtr(
0573             new reco::GenParticle(td->threeCharge(), td->p4(), td->vertex(), td->pdgId(), td->status(), false));
0574         target.push_back(*qPtr);
0575         // increment & push index of the top daughter
0576         topDaughters.push_back(++motherPartIdx_);
0577         if (addRadiation_) {
0578           // for radation to be added we first need to
0579           // pick the last quark in the MC chain
0580           const reco::GenParticle* last_q = findLastParticleInChain(static_cast<const reco::GenParticle*>(&*td));
0581           addRadiation(motherPartIdx_, last_q, target);
0582         }
0583       } else if (std::abs(td->pdgId()) == TopDecayID::WID) {
0584         // ladies and gentleman, we have a W boson
0585         std::unique_ptr<reco::GenParticle> WPtr(
0586             new reco::GenParticle(td->threeCharge(), td->p4(), td->vertex(), td->pdgId(), td->status(), false));
0587         target.push_back(*WPtr);
0588         // increment & push index of the top daughter
0589         topDaughters.push_back(++motherPartIdx_);
0590         iW = motherPartIdx_;
0591 
0592         // next are the daughers of our W boson
0593         // for Pythia 6 this is wrong as the last W has no daughters at all!
0594         // instead the status 3 W has 3 daughters: q qbar' and W (WTF??!)
0595         const reco::GenParticle* decaying_W = findLastParticleInChain(static_cast<const reco::GenParticle*>(&*td));
0596         for (reco::GenParticle::const_iterator wd = decaying_W->begin(); wd != decaying_W->end(); ++wd) {
0597           if (!(std::abs(wd->pdgId()) == TopDecayID::WID)) {
0598             std::unique_ptr<reco::GenParticle> qPtr(
0599                 new reco::GenParticle(wd->threeCharge(), wd->p4(), wd->vertex(), wd->pdgId(), wd->status(), false));
0600             target.push_back(*qPtr);
0601             // increment & push index of the top daughter
0602             wDaughters.push_back(++motherPartIdx_);
0603             const reco::GenParticle* last_q = findLastParticleInChain(static_cast<const reco::GenParticle*>(&*wd));
0604             addRadiation(motherPartIdx_, last_q, target);
0605             if (std::abs(wd->pdgId()) == TopDecayID::tauID) {
0606               // add tau daughters
0607               // currently it adds k-mesons etc as well, which
0608               // is not what we want.
0609               addDaughters(motherPartIdx_, wd, target);
0610             }
0611           }
0612         }
0613 
0614       } else {
0615         if (addRadiation_ && (td->pdgId() == TopDecayID::glueID || std::abs(td->pdgId()) < TopDecayID::bID)) {
0616           // collect additional radiation from the top
0617           std::unique_ptr<reco::GenParticle> radPtr(
0618               new reco::GenParticle(td->threeCharge(), td->p4(), td->vertex(), td->pdgId(), td->status(), false));
0619           target.push_back(*radPtr);
0620         }
0621         //other top daughters like Zq for FCNC
0622         // for pythia 6 many gluons end up here
0623         //std::cout << "other top daughters: to be implemented"
0624         //      << std::endl;
0625       }
0626     }
0627 
0628     // fill the map for the administration
0629     // of daughter indices
0630     refs_[iTop] = topDaughters;
0631     refs_[iW] = wDaughters;
0632   }
0633 }
0634 
0635 /// calculate lorentz vector from input
0636 reco::Particle::LorentzVector TopDecaySubset::p4(const std::vector<const reco::GenParticle*>::const_iterator top,
0637                                                  int statusFlag) {
0638   // calculate the four vector for top/anti-top quarks from
0639   // the W boson and the b quark plain or including all
0640   // additional radiation depending on switch 'plain'
0641   if (statusFlag == TopDecayID::unfrag) {
0642     // return 4 momentum as it is
0643     return (*top)->p4();
0644   }
0645   reco::Particle::LorentzVector vec;
0646   for (reco::GenParticle::const_iterator p = (*top)->begin(); p != (*top)->end(); ++p) {
0647     if (p->status() == TopDecayID::unfrag) {
0648       // descend by one level for each
0649       // status 3 particle on the way
0650       vec += p4(p, statusFlag);
0651     } else {
0652       if (std::abs((*top)->pdgId()) == TopDecayID::WID) {
0653         // in case of a W boson skip the status
0654         // 2 particle to prevent double counting
0655         if (std::abs(p->pdgId()) != TopDecayID::WID)
0656           vec += p->p4();
0657       } else {
0658         // add all four vectors for each stable
0659         // particle (status 1 or 2) on the way
0660         vec += p->p4();
0661         if (vec.mass() - (*top)->mass() > 0) {
0662           // continue adding up gluons and qqbar pairs on the top
0663           // line untill the nominal top mass is reached; then
0664           // break in order to prevent picking up virtualities
0665           break;
0666         }
0667       }
0668     }
0669   }
0670   return vec;
0671 }
0672 
0673 /// calculate lorentz vector from input (dedicated to top reconstruction)
0674 reco::Particle::LorentzVector TopDecaySubset::p4(const reco::GenParticle::const_iterator part, int statusFlag) {
0675   // calculate the four vector for all top daughters from
0676   // their daughters including additional radiation
0677   if (statusFlag == TopDecayID::unfrag) {
0678     // return 4 momentum as it is
0679     return part->p4();
0680   }
0681   reco::Particle::LorentzVector vec;
0682   for (reco::GenParticle::const_iterator p = part->begin(); p != part->end(); ++p) {
0683     if (p->status() <= TopDecayID::stable && std::abs(p->pdgId()) == TopDecayID::WID) {
0684       vec = p->p4();
0685     } else {
0686       if (p->status() <= TopDecayID::stable) {
0687         // sum up the p4 of all stable particles
0688         // (of status 1 or 2)
0689         vec += p->p4();
0690       } else {
0691         if (p->status() == TopDecayID::unfrag) {
0692           // if the particle is unfragmented (i.e.
0693           // status 3) descend by one level
0694           vec += p4(p, statusFlag);
0695         }
0696       }
0697     }
0698   }
0699   return vec;
0700 }
0701 
0702 /// fill vector including all radiations from quarks originating from W/top
0703 void TopDecaySubset::addRadiation(int& idx,
0704                                   const reco::GenParticle::const_iterator part,
0705                                   reco::GenParticleCollection& target) {
0706   std::vector<int> daughters;
0707   int idxBuffer = idx;
0708   for (reco::GenParticle::const_iterator daughter = part->begin(); daughter != part->end(); ++daughter) {
0709     if (daughter->status() <= TopDecayID::stable && daughter->pdgId() != part->pdgId()) {
0710       // skip comment lines and make sure that
0711       // the particle is not double counted as
0712       // daughter of itself
0713       std::unique_ptr<reco::GenParticle> ptr(new reco::GenParticle(
0714           daughter->threeCharge(), daughter->p4(), daughter->vertex(), daughter->pdgId(), daughter->status(), false));
0715       target.push_back(*ptr);
0716       daughters.push_back(++idx);  //push index of daughter
0717     }
0718   }
0719   if (!daughters.empty()) {
0720     refs_[idxBuffer] = daughters;
0721   }
0722 }
0723 
0724 void TopDecaySubset::addRadiation(int& idx, const reco::GenParticle* part, reco::GenParticleCollection& target) {
0725   std::vector<int> daughters;
0726   int idxBuffer = idx;
0727   for (reco::GenParticle::const_iterator daughter = part->begin(); daughter != part->end(); ++daughter) {
0728     // either pick daughters as long as they are different
0729     // to the initial particle
0730     if (daughter->pdgId() != part->pdgId()) {
0731       std::unique_ptr<reco::GenParticle> ptr(new reco::GenParticle(
0732           daughter->threeCharge(), daughter->p4(), daughter->vertex(), daughter->pdgId(), daughter->status(), false));
0733       target.push_back(*ptr);
0734       daughters.push_back(++idx);  //push index of daughter
0735     }
0736   }
0737   if (!daughters.empty()) {
0738     refs_[idxBuffer] = daughters;
0739   }
0740 }
0741 
0742 /// recursively fill vector for all further decay particles of a given particle
0743 void TopDecaySubset::addDaughters(int& idx,
0744                                   const reco::GenParticle::const_iterator part,
0745                                   reco::GenParticleCollection& target,
0746                                   bool recursive) {
0747   std::vector<int> daughters;
0748   int idxBuffer = idx;
0749   for (reco::GenParticle::const_iterator daughter = part->begin(); daughter != part->end(); ++daughter) {
0750     std::unique_ptr<reco::GenParticle> ptr(new reco::GenParticle(
0751         daughter->threeCharge(), daughter->p4(), daughter->vertex(), daughter->pdgId(), daughter->status(), false));
0752     target.push_back(*ptr);
0753     // increment & push index of daughter
0754     daughters.push_back(++idx);
0755     // continue recursively if desired
0756     if (recursive) {
0757       addDaughters(idx, daughter, target);
0758     }
0759   }
0760   if (!daughters.empty()) {
0761     refs_[idxBuffer] = daughters;
0762   }
0763 }
0764 
0765 /// clear references
0766 void TopDecaySubset::clearReferences() {
0767   // clear vector of references
0768   refs_.clear();
0769   // set idx for mother particles to a start value
0770   // of -1 (the first entry will raise it to 0)
0771   motherPartIdx_ = -1;
0772 }
0773 
0774 /// fill references for output vector
0775 void TopDecaySubset::fillReferences(const reco::GenParticleRefProd& ref, reco::GenParticleCollection& sel) {
0776   int idx = 0;
0777   for (reco::GenParticleCollection::iterator p = sel.begin(); p != sel.end(); ++p, ++idx) {
0778     //find daughter reference vectors in refs_ and add daughters
0779     std::map<int, std::vector<int> >::const_iterator daughters = refs_.find(idx);
0780     if (daughters != refs_.end()) {
0781       for (std::vector<int>::const_iterator daughter = daughters->second.begin(); daughter != daughters->second.end();
0782            ++daughter) {
0783         reco::GenParticle* part = dynamic_cast<reco::GenParticle*>(&(*p));
0784         if (part == nullptr) {
0785           throw edm::Exception(edm::errors::InvalidReference, "Not a GenParticle");
0786         }
0787         part->addDaughter(reco::GenParticleRef(ref, *daughter));
0788         sel[*daughter].addMother(reco::GenParticleRef(ref, idx));
0789       }
0790     }
0791   }
0792 }
0793 
0794 #include "FWCore/Framework/interface/MakerMacros.h"
0795 DEFINE_FWK_MODULE(TopDecaySubset);
[ 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!