CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_9_patch3/src/RecoJets/JetProducers/plugins/InputGenJetsParticleSelector.cc

Go to the documentation of this file.
00001 /* \class GenJetInputParticleSelector
00002 *
00003 *  Selects particles that are used as input for the GenJet collection.
00004 *  Logic: select all stable particles, except for particles specified in
00005 *  the config file that come from
00006 *  W,Z and H decays, and except for a special list, which can be used for
00007 *  unvisible BSM-particles.
00008 *  It is also possible to only selected the partonic final state, 
00009 *  which means all particles before the hadronization step.
00010 *
00011 *  The algorithm is based on code of Christophe Saout.
00012 *
00013 *  Usage: [example for no resonance from nu an mu, and deselect invisible BSM 
00014 *         particles ]
00015 *
00016 *  module genJetParticles = InputGenJetsParticleSelector {
00017 *                InputTag src = "genParticles"
00018 *                bool partonicFinalState = false  
00019 *                bool excludeResonances = true   
00020 *                vuint32 excludeFromResonancePids = {13,12,14,16}
00021 *                bool tausAsJets = false
00022 *                vuint32 ignoreParticleIDs = {   1000022, 2000012, 2000014,
00023 *                                                2000016, 1000039, 5000039,
00024 *                                                4000012, 9900012, 9900014,
00025 *                                                9900016, 39}
00026 *        }
00027 *
00028 *
00029 * \author: Christophe Saout, Andreas Oehler
00030 *
00031 */
00032 
00033 #include "InputGenJetsParticleSelector.h"
00034 #include "FWCore/Framework/interface/MakerMacros.h"
00035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00036 //#include <iostream>
00037 #include <memory>
00038 #include "CommonTools/CandUtils/interface/pdgIdUtils.h"
00039 
00040 using namespace std;
00041 
00042 InputGenJetsParticleSelector::InputGenJetsParticleSelector(const edm::ParameterSet &params ):
00043   inTag(params.getParameter<edm::InputTag>("src")),
00044   partonicFinalState(params.getParameter<bool>("partonicFinalState")),
00045   excludeResonances(params.getParameter<bool>("excludeResonances")),
00046   tausAsJets(params.getParameter<bool>("tausAsJets")),
00047   ptMin(0.0){
00048   if (params.exists("ignoreParticleIDs"))
00049     setIgnoredParticles(params.getParameter<std::vector<unsigned int> >
00050                         ("ignoreParticleIDs"));
00051   setExcludeFromResonancePids(params.getParameter<std::vector<unsigned int> >
00052                         ("excludeFromResonancePids"));
00053 
00054   produces <reco::GenParticleRefVector> ();
00055       
00056 }
00057 
00058 InputGenJetsParticleSelector::~InputGenJetsParticleSelector(){}
00059 
00060 void InputGenJetsParticleSelector::setIgnoredParticles(const std::vector<unsigned int> &particleIDs)
00061 {
00062   ignoreParticleIDs = particleIDs;
00063   std::sort(ignoreParticleIDs.begin(), ignoreParticleIDs.end());
00064 }
00065 
00066 void InputGenJetsParticleSelector::setExcludeFromResonancePids(const std::vector<unsigned int> &particleIDs)
00067 {
00068   excludeFromResonancePids = particleIDs;
00069   std::sort( excludeFromResonancePids.begin(), excludeFromResonancePids.end());
00070 }
00071 
00072 bool InputGenJetsParticleSelector::isParton(int pdgId) const{
00073   pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
00074   return (pdgId > 0 && pdgId < 6) || pdgId == 7 ||
00075     pdgId == 9 || (tausAsJets && pdgId == 15) || pdgId == 21;
00076   // tops are not considered "regular" partons
00077   // but taus eventually are (since they may hadronize later)
00078 }
00079 
00080 bool InputGenJetsParticleSelector::isHadron(int pdgId)
00081 {
00082   pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
00083   return (pdgId > 100 && pdgId < 900) ||
00084     (pdgId > 1000 && pdgId < 9000);
00085 }
00086 
00087 bool InputGenJetsParticleSelector::isResonance(int pdgId)
00088 {
00089   // gauge bosons and tops
00090   pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
00091   return (pdgId > 21 && pdgId <= 42) || pdgId == 6 || pdgId == 8 ;  //BUG! war 21. 22=gamma..
00092 }
00093 
00094 bool InputGenJetsParticleSelector::isIgnored(int pdgId) const
00095 {
00096   pdgId = pdgId > 0 ? pdgId : -pdgId;
00097   std::vector<unsigned int>::const_iterator pos =
00098     std::lower_bound(ignoreParticleIDs.begin(),
00099                      ignoreParticleIDs.end(),
00100                      (unsigned int)pdgId);
00101   return pos != ignoreParticleIDs.end() && *pos == (unsigned int)pdgId;
00102 }
00103 
00104 bool InputGenJetsParticleSelector::isExcludedFromResonance(int pdgId) const
00105 {
00106   pdgId = pdgId > 0 ? pdgId : -pdgId;
00107   std::vector<unsigned int>::const_iterator pos =
00108     std::lower_bound(excludeFromResonancePids.begin(),
00109                      excludeFromResonancePids.end(),
00110                      (unsigned int)pdgId);
00111   return pos != excludeFromResonancePids.end() && *pos == (unsigned int)pdgId;
00112  
00113 }
00114 
00115 static unsigned int partIdx(const InputGenJetsParticleSelector::ParticleVector &p,
00116                             //const reco::GenParticle *particle)
00117                             const reco::GenParticle *particle)
00118 {
00119   InputGenJetsParticleSelector::ParticleVector::const_iterator pos =
00120     std::lower_bound(p.begin(), p.end(), particle);
00121   if (pos == p.end() || *pos != particle)
00122     throw cms::Exception("CorruptedData")
00123       << "reco::GenEvent corrupted: Unlisted particles"
00124       " in decay tree." << std::endl;
00125 
00126   return pos - p.begin();
00127 }
00128     
00129 static void invalidateTree(InputGenJetsParticleSelector::ParticleBitmap &invalid,
00130                            const InputGenJetsParticleSelector::ParticleVector &p,
00131                            const reco::GenParticle *particle)
00132 {
00133   unsigned int npart=particle->numberOfDaughters();
00134   if (!npart) return;
00135 
00136   for (unsigned int i=0;i<npart;++i){
00137     unsigned int idx=partIdx(p,dynamic_cast<const reco::GenParticle*>(particle->daughter(i)));
00138     if (invalid[idx])
00139       continue;
00140     invalid[idx] = true;
00141     //cout<<"Invalidated: ["<<setw(4)<<idx<<"] With pt:"<<particle->daughter(i)->pt()<<endl;
00142     invalidateTree(invalid, p, dynamic_cast<const reco::GenParticle*>(particle->daughter(i)));
00143   }
00144 }
00145   
00146   
00147 int InputGenJetsParticleSelector::testPartonChildren
00148 (InputGenJetsParticleSelector::ParticleBitmap &invalid,
00149  const InputGenJetsParticleSelector::ParticleVector &p,
00150  const reco::GenParticle *particle) const
00151 {
00152   unsigned int npart=particle->numberOfDaughters();
00153   if (!npart) {return 0;}
00154 
00155   for (unsigned int i=0;i<npart;++i){
00156     unsigned int idx = partIdx(p,dynamic_cast<const reco::GenParticle*>(particle->daughter(i)));
00157     if (invalid[idx])
00158       continue;
00159     if (isParton((particle->daughter(i)->pdgId()))){
00160       return 1;
00161     }
00162     if (isHadron((particle->daughter(i)->pdgId()))){
00163       return -1;
00164     }
00165     int result = testPartonChildren(invalid,p,dynamic_cast<const reco::GenParticle*>(particle->daughter(i)));
00166     if (result) return result;
00167   }
00168   return 0;
00169 }
00170 
00171 InputGenJetsParticleSelector::ResonanceState
00172 InputGenJetsParticleSelector::fromResonance(ParticleBitmap &invalid,
00173                                                            const ParticleVector &p,
00174                                                            const reco::GenParticle *particle) const
00175 {
00176   unsigned int idx = partIdx(p, particle);
00177   int id = particle->pdgId();
00178 
00179   if (invalid[idx]) return kIndirect;
00180       
00181   if (isResonance(id) && particle->status() == 3){
00182     return kDirect;
00183   }
00184 
00185   
00186   if (!isIgnored(id) && (isParton(id)))
00187     return kNo;
00188   
00189   
00190   
00191   unsigned int nMo=particle->numberOfMothers();
00192   if (!nMo)
00193     return kNo;
00194   
00195 
00196   for(unsigned int i=0;i<nMo;++i){
00197     ResonanceState result = fromResonance(invalid,p,dynamic_cast<const reco::GenParticle*>(particle->mother(i)));
00198     switch(result) {
00199     case kNo:
00200       break;
00201     case kDirect:
00202       if (dynamic_cast<const reco::GenParticle*>(particle->mother(i))->pdgId()==id || isResonance(id))
00203         return kDirect;
00204       if(!isExcludedFromResonance(id))
00205         break;
00206     case kIndirect:
00207       return kIndirect;
00208     }
00209   }
00210 return kNo;
00211 }
00212 
00213     
00214 bool InputGenJetsParticleSelector::hasPartonChildren(ParticleBitmap &invalid,
00215                                                      const ParticleVector &p,
00216                                                      const reco::GenParticle *particle) const {
00217   return testPartonChildren(invalid, p, particle) > 0;
00218 }
00219     
00220 //######################################################
00221 //function NEEDED and called per EVENT by FRAMEWORK:
00222 void InputGenJetsParticleSelector::produce (edm::Event &evt, const edm::EventSetup &evtSetup){
00223 
00224  
00225   std::auto_ptr<reco::GenParticleRefVector> selected_ (new reco::GenParticleRefVector);
00226     
00227   edm::Handle<reco::GenParticleCollection> genParticles;
00228   //  evt.getByLabel("genParticles", genParticles );
00229   evt.getByLabel(inTag, genParticles );
00230     
00231     
00232   ParticleVector particles;
00233   for (reco::GenParticleCollection::const_iterator iter=genParticles->begin();iter!=genParticles->end();++iter){
00234     particles.push_back(&*iter); 
00235   }
00236   
00237   std::sort(particles.begin(), particles.end());
00238   unsigned int size = particles.size();
00239   
00240   ParticleBitmap selected(size, false);
00241   ParticleBitmap invalid(size, false);
00242 
00243   for(unsigned int i = 0; i < size; i++) {
00244     const reco::GenParticle *particle = particles[i];
00245     if (invalid[i])
00246       continue;
00247     if (particle->status() == 1)
00248       selected[i] = true;
00249     if (partonicFinalState && isParton(particle->pdgId())) {
00250           
00251       if (particle->numberOfDaughters()==0 &&
00252           particle->status() != 1) {
00253         // some brokenness in event...
00254         invalid[i] = true;
00255       }
00256       else if (!hasPartonChildren(invalid, particles,
00257                                   particle)) {
00258         selected[i] = true;
00259         invalidateTree(invalid, particles,particle); //this?!?
00260       }
00261     }
00262         
00263   }
00264  unsigned int count=0;
00265   for(size_t idx=0;idx<genParticles->size();++idx){ 
00266     const reco::GenParticle *particle = particles[idx];
00267     if (!selected[idx] || invalid[idx]){
00268       continue;
00269     }
00270         
00271     if (excludeResonances &&
00272         fromResonance(invalid, particles, particle)) {
00273       invalid[idx] = true;
00274       //cout<<"[INPUTSELECTOR] Invalidates FROM RESONANCE!: ["<<setw(4)<<idx<<"] "<<particle->pdgId()<<" "<<particle->pt()<<endl;
00275       continue;
00276     }
00277         
00278     if (isIgnored(particle->pdgId())){
00279       continue;
00280     }
00281 
00282    
00283     if (particle->pt() >= ptMin){
00284       edm::Ref<reco::GenParticleCollection> particleRef(genParticles,idx);
00285       selected_->push_back(particleRef);
00286       //cout<<"Finally we have: ["<<setw(4)<<idx<<"] "<<setw(4)<<particle->pdgId()<<" "<<particle->pt()<<endl;
00287       count++;
00288     }
00289   }
00290   evt.put(selected_);
00291 }
00292       
00293       
00294   
00295 //define this as a plug-in
00296 DEFINE_FWK_MODULE(InputGenJetsParticleSelector);