Reference Manual

Go to the documentation of this file.
00001 #include "GeneratorInterface/GenFilters/interface/PythiaFilterEMJetHeep.h"
00002 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00003 //#include "CLHEP/HepMC/GenParticle.h"
00004 
00005 #include <iostream>
00006 #include<list>
00007 #include<vector>
00008 #include<cmath>
00009 #include "TFile.h"
00010 
00011 // using namespace edm;
00012 // using namespace std;
00013 // using namespace HepMC;
00014 
00015 
00016 //namespace{
00017 
00018 double PythiaFilterEMJetHeep::deltaR(double eta0, double phi0, double eta, double phi){
00019     double dphi=phi-phi0;
00020     if(dphi>M_PI) dphi-=2*M_PI;
00021     else if(dphi<=-M_PI) dphi+=2*M_PI;
00022     return sqrt(dphi*dphi+(eta-eta0)*(eta-eta0));
00023 }
00024 
00025 struct ParticlePtGreater{
00026       bool operator()(const HepMC::GenParticle * a,
00027                       const HepMC::GenParticle * b) const{
00028       return a->momentum().perp() > b->momentum().perp();
00029     }
00030 };
00031 
00032 //}
00033 
00034 PythiaFilterEMJetHeep::PythiaFilterEMJetHeep(const edm::ParameterSet& iConfig) :
00035 label_(iConfig.getUntrackedParameter("moduleLabel",std::string("source"))),
00036 //
00037 minEventPt(iConfig.getUntrackedParameter<double>("MinEventPt",40.)),
00038 etaMax(iConfig.getUntrackedParameter<double>("MaxEta", 2.8)),
00039 cone_clust(iConfig.getUntrackedParameter<double>("ConeClust",0.10)), 
00040 cone_iso(iConfig.getUntrackedParameter<double>("ConeIso",0.50)), 
00041 nPartMin(iConfig.getUntrackedParameter<unsigned int>("NumPartMin", 2)),
00042 drMin(iConfig.getUntrackedParameter<double>("dRMin", 0.4)),
00043 //
00044 eventsProcessed(0),
00045 maxnumberofeventsinrun(iConfig.getUntrackedParameter<int>("MaxEvents",1000)),
00046 outputFile_(iConfig.getUntrackedParameter("outputFile",std::string("rootOutputFile"))),
00047 minbias(iConfig.getUntrackedParameter<bool>("Minbias",false)),
00048 debug(iConfig.getUntrackedParameter<bool>("Debug",true)) { 
00049 theNumberOfSelected = 0;   
00050 //  
00051 }
00052 
00053 PythiaFilterEMJetHeep::~PythiaFilterEMJetHeep(){}
00054 
00055 void PythiaFilterEMJetHeep::beginJob() {
00056 
00057   // parametarizations of presel. criteria:
00058 
00059   ptSeedMin_EB =  6.0 + (minEventPt - 80.) * 0.035;
00060   ptSeedMin_EE =  5.5 + (minEventPt - 80.) * 0.033;
00061   fracConePtMin_EB = 0.60 - (minEventPt - 80.) * 0.0009;
00062   fracConePtMin_EE = fracConePtMin_EB;
00063   fracEmPtMin_EB = 0.30 + (minEventPt - 80.) * 0.0017;
00064   if(minEventPt >=225) fracEmPtMin_EB = 0.40 + (minEventPt - 230.) * 0.00063;
00065   fracEmPtMin_EE = 0.30 + (minEventPt - 80.) * 0.0033;
00066   if(minEventPt >=165) fracEmPtMin_EE = 0.55 + (minEventPt - 170.) * 0.0005;
00067   fracTrkPtMax_EB = 0.80 - (minEventPt - 80.) * 0.001;
00068   fracTrkPtMax_EE = 0.70 - (minEventPt - 80.) * 0.001;
00069   isoConeMax_EB = 0.35;
00070   isoConeMax_EE = 0.40;
00071   ptHdMax_EB = 40;
00072   ptHdMax_EB = 45;
00073   ntrkMax_EB = 4;
00074   ntrkMax_EE = 4;
00075 
00076   if( minbias ) {
00077 
00078     std::cout <<" ... Minbias preselection ... " << std::endl;
00079 
00080     minEventPt = 1.0;
00081     ptSeedMin_EB =  1.5; 
00082     ptSeedMin_EE =  1.5; 
00083     fracConePtMin_EB = 0.20; 
00084     fracConePtMin_EE = fracConePtMin_EB;
00085     fracEmPtMin_EB = 0.20; 
00086     fracEmPtMin_EE = 0.20; 
00087     fracTrkPtMax_EB = 0.80; 
00088     fracTrkPtMax_EE = 0.80; 
00089     isoConeMax_EB = 1.0;
00090     isoConeMax_EE = 1.0;
00091     ptHdMax_EB = 7;
00092     ptHdMax_EB = 7;
00093     ntrkMax_EB = 2;
00094     ntrkMax_EE = 2;
00095 
00096   }
00097 
00098 
00099   std::cout <<" minEventPt = " << minEventPt << std::endl;
00100   std::cout <<" etaMax = " << etaMax << std::endl;
00101   std::cout <<" nPartMin = " << nPartMin << std::endl;
00102   std::cout <<" cone_clust = " << cone_clust << std::endl;
00103   std::cout <<" cone_iso = " << cone_iso << std::endl;
00104   std::cout <<" drMin = " << drMin << std::endl << std::endl;
00105   std::cout <<" ptSeedMin_EB = " << ptSeedMin_EB << std::endl;
00106   std::cout <<" ptSeedMin_EE = " << ptSeedMin_EE << std::endl;
00107   std::cout <<" fracConePtMin_EB = " << fracConePtMin_EB << std::endl;
00108   std::cout <<" fracConePtMin_EE = " << fracConePtMin_EE << std::endl;
00109   std::cout <<" fracEmPtMin_EB = " << fracEmPtMin_EB << std::endl;
00110   std::cout <<" fracEmPtMin_EE = " << fracEmPtMin_EE << std::endl;
00111   std::cout <<" fracTrkPtMax_EB = " << fracTrkPtMax_EB << std::endl;
00112   std::cout <<" fracTrkPtMax_EE = " << fracTrkPtMax_EE << std::endl;
00113   std::cout <<" ntrkMax_EB = " << ntrkMax_EB << std::endl;
00114   std::cout <<" ntrkMax_EE = " << ntrkMax_EE << std::endl;
00115   std::cout <<" ptHdMax_EB = " << ptHdMax_EB << std::endl;
00116   std::cout <<" ptHdMax_EE = " << ptHdMax_EE << std::endl;
00117   std::cout <<" isoConeMax_EB = " << isoConeMax_EB << std::endl;
00118   std::cout <<" isoConeMax_EE = " << isoConeMax_EE << std::endl;
00119 
00120 }
00121 
00122 // ------------ method called to produce the data  ------------
00123 bool PythiaFilterEMJetHeep::filter(edm::Event& iEvent, const edm::EventSetup& iSetup){
00124 
00125   if(theNumberOfSelected>=maxnumberofeventsinrun)   {
00126     throw cms::Exception("endJob")<<"we have reached the maximum number of events ";
00127   }
00128 
00129   accepted = false;
00130   edm::Handle<edm::HepMCProduct> evt;
00131   iEvent.getByLabel(label_, evt);
00132 
00133   const HepMC::GenEvent * myGenEvent = evt->GetEvent();
00134 
00135   eventsProcessed++;
00136 
00137   std::list<const HepMC::GenParticle *> seeds;
00138   std::list<const HepMC::GenParticle *> candidates;
00139  
00140 // collect the seeds
00141   for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();   p != myGenEvent->particles_end(); ++p ) {
00142    
00143     double pt=(*p)->momentum().perp();
00144     double eta=(*p)->momentum().eta();
00145 
00146     int pdgid = (*p)->pdg_id(); 
00147     if ( !(pdgid==22 || abs(pdgid)==11 || abs(pdgid)==211 || abs(pdgid)==321) ) continue;
00148     //if ( !(pdgid==22 || abs(pdgid)==11) ) continue;
00149  
00150 //  Selection #1: seed particle ETA
00151     if ( std::abs(eta) > etaMax ) continue;
00152     bool EB = false;
00153     bool EE = false;
00154     if (fabs(eta)<1.5)                         EB = true;
00155     else if (fabs(eta)>=1.5&&fabs(eta)<etaMax) EE = true;   
00156     if (debug) std::cout << " Selection 1 passed " << std::endl;
00157 
00158 //  Selection #2: seed particle pT
00159     if ( EB && pt < ptSeedMin_EB ) continue;
00160     if ( EE && pt < ptSeedMin_EE ) continue;
00161     if (debug) std::cout << " Selection 2 passed " << std::endl;
00162 
00163     seeds.push_back(*p);   
00164   }
00165   
00166   if (seeds.size() < nPartMin ) return false;
00167 
00168   seeds.sort(ParticlePtGreater());
00169 
00170 // select the proto-clusters
00171   std::list<const HepMC::GenParticle*>::iterator itSeed;
00172 
00173   for(itSeed = seeds.begin(); itSeed != seeds.end(); ++itSeed) {
00174 
00175     double pt=(*itSeed)->momentum().perp();
00176     double eta=(*itSeed)->momentum().eta();
00177     double phi=(*itSeed)->momentum().phi();  
00178 
00179     bool EB = false;
00180     bool EE = false;
00181     if (fabs(eta)<1.5)                         EB = true;
00182     else if (fabs(eta)>=1.5&&fabs(eta)<etaMax) EE = true;     
00183 
00184     float setCone_iso=0;
00185     float setCone_clust=0;
00186     float setEM=0;
00187     float setHAD=0;
00188     float ptMaxHadron=0;
00189     float setCharged=0;
00190     unsigned int Ncharged = 0;
00191     
00192     for ( HepMC::GenEvent::particle_const_iterator pp = myGenEvent->particles_begin();  pp != myGenEvent->particles_end(); ++pp ) {
00193     
00194       if ( (*pp)->status()!=1 ) continue;  // select just stable particles
00195       int pid= (*pp)->pdg_id();
00196       int apid= std::abs(pid);
00197       if(apid==310) apid = 22; 
00198       if (apid > 11 && apid < 20) continue; //get rid of muons and neutrinos
00199 
00200       double pt_ =(*pp)->momentum().perp();
00201       double eta_=(*pp)->momentum().eta();
00202       double phi_=(*pp)->momentum().phi();
00203 
00204       float dr = deltaR(eta_, phi_, eta, phi);
00205       if ( dr <= cone_iso )  setCone_iso += pt_;
00206 
00207       bool charged = false;
00208       if ( apid ==211 || apid == 321 || apid == 2212 || apid == 11 ) charged = true;
00209       if ( dr <= cone_clust ) {
00210          setCone_clust += pt_;
00211          if ( apid == 22 || apid == 11) setEM += pt_;
00212          if ( apid>100 )   setHAD += pt_;
00213          if ( apid>100 && pt_ > ptMaxHadron ) ptMaxHadron = pt_;
00214          if ( charged && pt_ > 1 ) { Ncharged++; setCharged += pt_;}
00215       }     
00216     }
00217     setCone_iso -= setCone_clust;
00218 
00219     if ( pt / setCone_clust < 0.15 ) continue;
00220 
00221     // Selection #3: min. pT of all particles in the proto-cluster   
00222     if (EB && setCone_clust < fracConePtMin_EB*minEventPt ) continue;
00223     if (EE && setCone_clust < fracConePtMin_EE*minEventPt ) continue;
00224     if (debug) std::cout << " Selection 3 passed " << std::endl;
00225 
00226      // Selections #4: min/max pT fractions of e/gamma in the proto-cluster from the total pT 
00227     if ( EB && setEM / setCone_clust  < fracEmPtMin_EB ) continue;
00228     if ( EE && setEM / setCone_clust  < fracEmPtMin_EE ) continue;
00229     if (debug) std::cout << " Selection 4 passed " << std::endl;
00230 
00231     // Selection 5: max. track pT fractions and number of tracks in the proto-cluster
00232     if ( (EB && setCharged / setCone_clust > fracTrkPtMax_EB) ||  Ncharged > ntrkMax_EB) continue;
00233     if ( (EE && setCharged / setCone_clust > fracTrkPtMax_EE) ||  Ncharged > ntrkMax_EE) continue;
00234     if (debug) std::cout << " Selection 5 passed " << std::endl;
00235 
00236     // Selection #6: max. pT of the hadron in the proto-cluster
00237     if ( EB && ptMaxHadron > ptHdMax_EB ) continue;
00238     if ( EE && ptMaxHadron > ptHdMax_EE ) continue;
00239     if (debug) std::cout << " Selection 6 passed " << std::endl;
00240 
00241     // Selection #7: max. pT fraction in the dR=cone_iso around the proto-cluster
00242     if ( EB && setCone_iso / setCone_clust > isoConeMax_EB) continue;
00243     if ( EE && setCone_iso / setCone_clust > isoConeMax_EE) continue;
00244     if (debug) std::cout << " Selection 7 passed " << std::endl;       
00245 
00246     if (debug) {
00247       std::cout <<  "(*itSeed)->pdg_id() = " << (*itSeed)->pdg_id() << std::endl;
00248       std::cout <<  "pt = " << pt << std::endl;
00249       std::cout <<  "setCone_clust = " << setCone_clust << std::endl;
00250       std::cout <<  "setEM = " << setEM << std::endl;
00251       std::cout <<  "ptMaxHadron = " << ptMaxHadron << std::endl;
00252       std::cout <<  "setCone_iso = " << setCone_iso << std::endl;
00253       std::cout <<  "setCone_iso / setCone_clust = " << setCone_iso / setCone_clust << std::endl;    
00254     }
00255 
00256     // Selection #8: any two objects should be separated by dr > drMin
00257     bool same = false;
00258     std::list<const HepMC::GenParticle*>::iterator itPart;
00259     for(itPart = candidates.begin(); itPart != candidates.end(); ++itPart) {
00260       float eta_ = (*itPart)->momentum().eta();
00261       float phi_ = (*itPart)->momentum().phi();
00262       float dr = deltaR(eta_, phi_, eta, phi);    
00263       if (dr < drMin ) same = true;
00264     }
00265     if (same) continue;
00266     if (debug) std::cout << " Selection 8 passed " << std::endl;       
00267 
00268     candidates.push_back(*itSeed);
00269 
00270   }  
00271 
00272   if (candidates.size() >= nPartMin ) {
00273     accepted=true;
00274     theNumberOfSelected++;
00275   }
00276 
00277   if (debug) std::cout<<" Event preselected "<<theNumberOfSelected<<" Proccess ID "<<myGenEvent->signal_process_id()<<std::endl;
00278 
00279   return accepted; 
00280   
00281 }
00282 
00283 void PythiaFilterEMJetHeep::endJob() {
00284 
00285   std::cout<<"\n *************************\n "<<std::endl;
00286   std::cout<<" Events processed "<< eventsProcessed <<std::endl;
00287   std::cout<<" Events preselected "<< theNumberOfSelected <<std::endl;
00288   std::cout<<"\n ************************* \n"<<std::endl;
00289 
00290 }
00291
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/GeneratorInterface/GenFilters/src/PythiaFilterEMJetHeep.cc
