Reference Manual

Go to the documentation of this file.
00001 #include "GeneratorInterface/GenFilters/interface/PythiaFilterGammaGamma.h"
00002 #include "DataFormats/Math/interface/LorentzVector.h"
00003 #include "Math/GenVector/VectorUtil.h"
00004 //#include "CLHEP/HepMC/GenParticle.h"
00005 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00006 
00007 #include "TFile.h"
00008 #include "TLorentzVector.h"
00009 
00010 #include <iostream>
00011 
00012 using namespace edm;
00013 using namespace std;
00014 using namespace HepMC;
00015 
00016 
00017 PythiaFilterGammaGamma::PythiaFilterGammaGamma(const edm::ParameterSet& iConfig) :
00018   label(iConfig.getUntrackedParameter<std::string>("moduleLabel",std::string("generator"))),
00019   //fileName(iConfig.getUntrackedParameter<std::string>("fileName", std::string("plots.root"))),
00020   maxEvents(iConfig.getUntrackedParameter<int>("maxEvents", 100000)),
00021   ptSeedThr(iConfig.getUntrackedParameter<double>("PtSeedThr")),
00022   etaSeedThr(iConfig.getUntrackedParameter<double>("EtaSeedThr")),
00023   ptGammaThr(iConfig.getUntrackedParameter<double>("PtGammaThr")),
00024   etaGammaThr(iConfig.getUntrackedParameter<double>("EtaGammaThr")),
00025   ptTkThr(iConfig.getUntrackedParameter<double>("PtTkThr")),
00026   etaTkThr(iConfig.getUntrackedParameter<double>("EtaTkThr")),
00027   ptElThr(iConfig.getUntrackedParameter<double>("PtElThr")),
00028   etaElThr(iConfig.getUntrackedParameter<double>("EtaElThr")),
00029   dRTkMax(iConfig.getUntrackedParameter<double>("dRTkMax")),
00030   dRSeedMax(iConfig.getUntrackedParameter<double>("dRSeedMax")),
00031   dPhiSeedMax(iConfig.getUntrackedParameter<double>("dPhiSeedMax")),
00032   dEtaSeedMax(iConfig.getUntrackedParameter<double>("dEtaSeedMax")),
00033   dRNarrowCone(iConfig.getUntrackedParameter<double>("dRNarrowCone")),
00034   pTMinCandidate1(iConfig.getUntrackedParameter<double>("PtMinCandidate1")),
00035   pTMinCandidate2(iConfig.getUntrackedParameter<double>("PtMinCandidate2")),
00036   etaMaxCandidate(iConfig.getUntrackedParameter<double>("EtaMaxCandidate")),
00037   invMassWide(iConfig.getUntrackedParameter<double>("InvMassWide")),
00038   invMassNarrow(iConfig.getUntrackedParameter<double>("InvMassNarrow")),
00039   nTkConeMax(iConfig.getUntrackedParameter<int>("NTkConeMax")),
00040   nTkConeSum(iConfig.getUntrackedParameter<int>("NTkConeSum")),
00041   acceptPrompts(iConfig.getUntrackedParameter<bool>("AcceptPrompts")), 
00042   promptPtThreshold(iConfig.getUntrackedParameter<double>("PromptPtThreshold")) {
00043 
00044   //cout<<"ptSeedThr "<<ptSeedThr<<endl;
00045   //cout<<"etaSeedThr "<<etaSeedThr<<endl;
00046   //cout<<"ptGammaThr "<<ptGammaThr<<endl;
00047   //cout<<"etaGammaThr "<<etaGammaThr<<endl;
00048   //cout<<"ptTkThr "<<ptTkThr<<endl;
00049   //cout<<"etaTkThr "<<etaTkThr<<endl;
00050   //cout<<"ptElThr "<<ptElThr<<endl;
00051   //cout<<"etaElThr "<<etaElThr<<endl;
00052   //cout<<"dRTkMax "<<dRTkMax<<endl;
00053   //cout<<"dRSeedMax "<<dRSeedMax<<endl;
00054   //cout<<"dPhiSeedMax "<<dPhiSeedMax<<endl;
00055   //cout<<"dEtaSeedMax "<<dEtaSeedMax<<endl;
00056   //cout<<"dRNarrowCone "<<dRNarrowCone<<endl;
00057   //cout<<"pTMinCandidate1 "<<pTMinCandidate1<<endl;
00058   //cout<<"pTMinCandidate2 "<<pTMinCandidate2<<endl;
00059   //cout<<"etaMaxCandidate "<<etaMaxCandidate<<endl;
00060   //cout<<"invMassWide "<<invMassWide<<endl;
00061   //cout<<"invMassNarrow "<<invMassNarrow<<endl;
00062   //cout<<"nTkConeMax "<<nTkConeMax<<endl;
00063   //cout<<"nTkConeSum "<<nTkConeSum<<endl;
00064   //cout<<"acceptPrompts "<<acceptPrompts<<endl;
00065   //cout<<"promptPtThreshold "<<promptPtThreshold<<endl;
00066   
00067   nSelectedEvents = 0;
00068   nGeneratedEvents = 0;
00069 
00070   //char a[20];
00071   //for(int i=0; i<2; i++) {
00072   //  sprintf(a, "PT Seed %d", i+1);
00073   //  hPtSeed[i] = new TH1D(a, a, 100, 0, 200);
00074   //  sprintf(a, "Eta Seed %d", i+1);
00075   //  hEtaSeed[i]= new TH1D(a, a, 100, -3., 3.);
00076   //  sprintf(a, "Pid Seed %d", i+1);
00077   //  hPidSeed[i]= new TH1I(a, a, 50, 0, 500); 
00078   //  sprintf(a, "PT Candidate %d", i+1);
00079   //  hPtCandidate[i] = new TH1D(a, a, 100, 0, 200);
00080   //  sprintf(a, "Eta Candidate %d", i+1);
00081   //  hEtaCandidate[i]= new TH1D(a, a, 100, -3., 3.);
00082   //  sprintf(a, "Pid Candidate %d", i+1);
00083   //  hPidCandidate[i]= new TH1I(a, a, 50, 0, 500);
00084   //  sprintf(a, "NTk Iso %d", i+1);
00085   //  hNTk[i]= new TH1I(a, a, 50, 0, 50);
00086   //}
00087   //hMassNarrow = new TH1D("Mass narr.", "Mass narr.", 100, 0, 1000);
00088   //hMassWide= new TH1D("Mass wide", "Mass wide", 100, 0, 1000);
00089   //hNTkSum= new TH1I("NTk Sum Iso", "NTk Sum Iso", 50, 0, 50);
00090   
00091 }
00092 
00093 PythiaFilterGammaGamma::~PythiaFilterGammaGamma() 
00094 {  
00095   //writeFile();
00096   cout << "Number of Selected Events: " << nSelectedEvents << endl;
00097   cout << "Number of Generated Events: " << nGeneratedEvents << endl;
00098 }
00099 
00100 //void PythiaFilterGammaGamma::writeFile() {
00101 
00102   //TFile* file = new TFile (fileName.c_str(), "recreate");
00103   
00104   //for(int i=0; i<2; i++) {
00105   //  hPtSeed[i]->Write();
00106   //  hEtaSeed[i]->Write();
00107   //  hPidSeed[i]->Write(); 
00108   //  hPtCandidate[i]->Write();
00109   //  hEtaCandidate[i]->Write();
00110   //  hPidCandidate[i]->Write();
00111   //  hNTk[i]->Write();
00112   //}
00113   //hMassNarrow->Write();
00114   //hMassWide->Write();
00115   //hNTkSum->Write();
00116   
00117   //file->Close();
00118   
00119 //}
00120 
00121 bool PythiaFilterGammaGamma::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
00122 
00123   if(nSelectedEvents >= maxEvents) {
00124     //writeFile();
00125     throw cms::Exception("endJob")<<"We have reached the maximum number of events...";
00126   }
00127 
00128   nGeneratedEvents++;
00129 
00130   bool accepted = false;
00131 
00132   Handle<HepMCProduct> evt;
00133   iEvent.getByLabel(label, evt);
00134   myGenEvent = evt->GetEvent();
00135 
00136   std::vector<const GenParticle*> seeds, egamma, stable; 
00137   std::vector<const GenParticle*>::const_iterator itPart, itStable, itEn;
00138 
00139   for(HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
00140 
00141     if (
00142         ((*p)->status()==1&&(*p)->pdg_id() == 22) ||                   // gamma
00143         ((*p)->status()==1&&abs((*p)->pdg_id()) == 11) ||              // electron
00144         (*p)->pdg_id() == 111 ||                  // pi0
00145         abs((*p)->pdg_id()) == 221 ||             // eta
00146         abs((*p)->pdg_id()) == 331 ||             // eta prime
00147         abs((*p)->pdg_id()) == 113 ||             // rho0 
00148         abs((*p)->pdg_id()) == 223)               // omega
00149       {       
00150         // check for eta and pT threshold for seed in gamma, el
00151         if ((*p)->momentum().perp() > ptSeedThr &&
00152             fabs((*p)->momentum().eta()) < etaSeedThr) {
00153         
00154           // check if found is daughter of one already taken
00155           bool isUsed = false;
00156           
00157           const GenParticle* mother = (*p)->production_vertex() ?
00158             *((*p)->production_vertex()->particles_in_const_begin()) : 0;
00159           const GenParticle* motherMother = (mother != 0  && mother->production_vertex()) ?
00160             *(mother->production_vertex()->particles_in_const_begin()) : 0;
00161           const GenParticle* motherMotherMother = (motherMother != 0 && motherMother->production_vertex()) ?
00162             *(motherMother->production_vertex()->particles_in_const_begin()) : 0;
00163 
00164           for(itPart = seeds.begin(); itPart != seeds.end(); itPart++) {
00165             
00166             if ((*itPart) == mother ||
00167                 (*itPart) == motherMother ||
00168                 (*itPart) == motherMotherMother) {
00169               isUsed = true;
00170               break;
00171             }
00172           }
00173           
00174           if (!isUsed) seeds.push_back(*p);
00175         }
00176       }  
00177   } 
00178    
00179   if (seeds.size() < 2) return accepted;
00180 
00181   for(HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
00182 
00183     if ((*p)->status() == 1) {
00184 
00185       // save charged stable tracks
00186       if (abs((*p)->pdg_id()) == 211 ||
00187           abs((*p)->pdg_id()) == 321 ||
00188           abs((*p)->pdg_id()) == 11 ||
00189           abs((*p)->pdg_id()) == 13 ||
00190           abs((*p)->pdg_id()) == 15) {
00191         // check if it passes the cut
00192         if ((*p)->momentum().perp() > ptTkThr &&
00193             fabs((*p)->momentum().eta()) < etaTkThr) {
00194           stable.push_back(*p);
00195         }
00196       }
00197 
00198       // save egamma for candidate calculation
00199       if ((*p)->pdg_id() == 22 &&
00200           (*p)->momentum().perp() > ptGammaThr &&
00201           fabs((*p)->momentum().eta()) < etaGammaThr) {
00202         egamma.push_back(*p);
00203       }
00204       if (abs((*p)->pdg_id()) == 11 &&
00205           (*p)->momentum().perp() > ptElThr &&
00206           fabs((*p)->momentum().eta()) < etaElThr) {
00207         egamma.push_back(*p); 
00208       }
00209     }
00210   }
00211 
00212   std::vector<int> nTracks;
00213   std::vector<TLorentzVector> candidate, candidateNarrow, candidateSeed;
00214   std::vector<const GenParticle*>::iterator itSeed;
00215 
00216   for(itSeed = seeds.begin(); itSeed != seeds.end(); ++itSeed) {
00217 
00218     TLorentzVector energy, narrowCone, temp1, temp2, tempseed;
00219 
00220     tempseed.SetXYZM((*itSeed)->momentum().px(), (*itSeed)->momentum().py(), (*itSeed)->momentum().pz(), 0);   
00221     for(itEn = egamma.begin(); itEn != egamma.end(); ++itEn) {
00222       temp1.SetXYZM((*itEn)->momentum().px(), (*itEn)->momentum().py(), (*itEn)->momentum().pz(), 0);        
00223         
00224       double DR = temp1.DeltaR(tempseed);
00225       double DPhi = temp1.DeltaPhi(tempseed);
00226       double DEta = (*itEn)->momentum().eta()-(*itSeed)->momentum().eta();
00227       if(DPhi<0) DPhi=-DPhi;
00228       if(DEta<0) DEta=-DEta;
00229         
00230       if (DR < dRSeedMax || (DPhi<dPhiSeedMax&&DEta<dEtaSeedMax)) {
00231         energy += temp1;
00232       }
00233       if (DR < dRNarrowCone) {
00234         narrowCone += temp1;
00235       }
00236     }
00237 
00238     int counter = 0;
00239     //bool isIsolated = true;
00240 
00241     if ( energy.Et() != 0. ) {
00242       if (fabs(energy.Eta()) < etaMaxCandidate) {
00243 
00244         temp2.SetXYZM(energy.Px(), energy.Py(), energy.Pz(), 0);        
00245         
00246         for(itStable = stable.begin(); itStable != stable.end(); ++itStable) {  
00247           temp1.SetXYZM((*itStable)->momentum().px(), (*itStable)->momentum().py(), (*itStable)->momentum().pz(), 0);        
00248           double DR = temp1.DeltaR(temp2);
00249           if (DR < dRTkMax) counter++;        
00250         }
00251 
00252         if(acceptPrompts) {
00253           bool isPrompt=false;
00254           if((*itSeed)->status() == 1&&(*itSeed)->pdg_id() == 22) {
00255             const GenParticle* mother = (*itSeed)->production_vertex() ?
00256               *((*itSeed)->production_vertex()->particles_in_const_begin()) : 0;
00257             if(mother) {
00258               if(mother->pdg_id()>=-22&&mother->pdg_id()<=22) {
00259                 const GenParticle* motherMother = (mother != 0  && mother->production_vertex()) ?
00260                   *(mother->production_vertex()->particles_in_const_begin()) : 0;
00261                 if(motherMother) {
00262                   if(motherMother->pdg_id()>=-22&&motherMother->pdg_id()<=22) {
00263                     if((*itSeed)->momentum().perp()>promptPtThreshold) {
00264                       isPrompt=true;
00265                     }
00266                   }
00267                 }
00268               }
00269             }
00270           }
00271           if(isPrompt) counter=0;
00272         }
00273         // check number of tracks
00274         //if (counter <= nTkConeMax) isIsolated = true;
00275       }
00276     }
00277 
00278     // check pt candidate, check nTrack, check eta
00279     //if (isIsolated) {
00280     candidate.push_back(energy);
00281     candidateSeed.push_back(tempseed);
00282     candidateNarrow.push_back(narrowCone);
00283     nTracks.push_back(counter);
00284     //++itSeed;
00285     //} 
00286   }
00287 
00288   if (candidate.size() <2) return accepted;
00289 
00290   TLorentzVector minv, minvNarrow;
00291   
00292   
00293   //bool filled=false;
00294 
00295   int i1, i2;
00296   for(unsigned int i=0; i<candidate.size()-1; ++i) {
00297     
00298     if (candidate[i].Energy() <1.) continue;
00299     if(nTracks[i]>nTkConeMax) continue;
00300     if (fabs(candidate[i].Eta()) > etaMaxCandidate) continue;
00301     
00302     for(unsigned int j=i+1; j<candidate.size(); ++j) { 
00303       if (candidate[j].Energy() <1.) continue;
00304       if(nTracks[j]>nTkConeMax) continue;
00305       if (fabs(candidate[j].Eta()) > etaMaxCandidate) continue;
00306           
00307       if (nTracks[i] + nTracks[j] > nTkConeSum) continue;
00308 
00309       if (candidate[i].Pt() > candidate[j].Pt()) {
00310         i1 = i;
00311         i2 = j;
00312       } 
00313       else {
00314         i1 = j;
00315         i2 = i;
00316       }
00317 
00318       if (candidate[i1].Pt() < pTMinCandidate1 || candidate[i2].Pt() < pTMinCandidate2) continue;
00319 
00320       minv = candidate[i] + candidate[j];
00321       if (minv.M() < invMassWide) continue;
00322         
00323       minvNarrow = candidateNarrow[i] + candidateNarrow[j];
00324       if (minvNarrow.M() > invMassNarrow) continue;
00325 
00326       accepted = true;
00327 
00328       //if(!filled) {
00329       //hMassWide->Fill(minv.M());
00330       //hMassNarrow->Fill(minvNarrow.M());
00331       //hNTkSum->Fill(nTracks[i] + nTracks[j]);
00332       //hPtCandidate[0]->Fill(candidate[i1].Pt());
00333       //hPtCandidate[1]->Fill(candidate[i2].Pt());
00334       //hEtaCandidate[0]->Fill(candidate[i1].Eta());
00335       //hEtaCandidate[1]->Fill(candidate[i2].Eta());
00336       //hPidCandidate[0]->Fill(seeds[i1]->pdg_id());
00337       //hPidCandidate[1]->Fill(seeds[i2]->pdg_id());
00338       //hPtSeed[0]->Fill(candidateSeed[i1].Pt());
00339       //hPtSeed[1]->Fill(candidateSeed[i2].Pt());
00340       //hEtaSeed[0]->Fill(candidateSeed[i1].Eta());
00341       //hEtaSeed[1]->Fill(candidateSeed[i2].Eta());
00342       //hPidSeed[0]->Fill(seeds[i1]->pdg_id());
00343       //hPidSeed[1]->Fill(seeds[i2]->pdg_id());
00344       //hNTk[0]->Fill(nTracks[i1]);
00345       //hNTk[1]->Fill(nTracks[i2]);
00346       //filled=true;
00347       //}
00348     }
00349   }
00350   
00351   if (accepted) nSelectedEvents++;
00352   return accepted;
00353 }
00354
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/GeneratorInterface/GenFilters/src/PythiaFilterGammaGamma.cc
