Skip to content

Commit

Permalink
Merge pull request cms-sw#6976 from mmachet/CMSSW_7_4_X
Browse files Browse the repository at this point in the history 
made EM enriched filter more indipendent from the generator used and add...
  • Loading branch information
@cmsbuild @bendavid
cmsbuild authored and bendavid committed Jan 25, 2015
1 parent 119684e commit e39617c
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 168 deletions.
13 changes: 3 additions & 10 deletions GeneratorInterface/GenFilters/interface/PythiaFilterGammaGamma.h
View file
Expand Up @@ -26,8 +26,6 @@ class PythiaFilterGammaGamma : public edm::EDFilter {
explicit PythiaFilterGammaGamma(const edm::ParameterSet&);
~PythiaFilterGammaGamma();

//void writeFile();

virtual bool filter(edm::Event&, const edm::EventSetup&);
private:

Expand All @@ -39,20 +37,15 @@ class PythiaFilterGammaGamma : public edm::EDFilter {
double minetacut;
double maxetacut;
int maxEvents;
int nSelectedEvents, nGeneratedEvents;
int nSelectedEvents, nGeneratedEvents, counterPrompt;

double ptSeedThr, etaSeedThr, ptGammaThr, etaGammaThr, ptTkThr, etaTkThr;
double ptElThr, etaElThr, dRTkMax, dRSeedMax, dPhiSeedMax, dEtaSeedMax, dRNarrowCone, pTMinCandidate1, pTMinCandidate2, etaMaxCandidate;
double invMassWide, invMassNarrow;
double invMassMin, invMassMax;
double energyCut;
int nTkConeMax, nTkConeSum;
bool acceptPrompts;
double promptPtThreshold;

//std::string fileName;
//TH1D *hPtSeed[2], *hEtaSeed[2], *hMassNarrow, *hMassWide;
//TH1I *hPidSeed[2], *hNTk[2], *hSel, *hNTkSum;
//TH1D *hPtCandidate[2], *hEtaCandidate[2];
//TH1I *hPidCandidate[2];

};
#endif
207 changes: 49 additions & 158 deletions GeneratorInterface/GenFilters/src/PythiaFilterGammaGamma.cc
View file
Expand Up @@ -16,7 +16,6 @@ using namespace HepMC;

PythiaFilterGammaGamma::PythiaFilterGammaGamma(const edm::ParameterSet& iConfig) :
label(iConfig.getUntrackedParameter<std::string>("moduleLabel",std::string("generator"))),
//fileName(iConfig.getUntrackedParameter<std::string>("fileName", std::string("plots.root"))),
maxEvents(iConfig.getUntrackedParameter<int>("maxEvents", 100000)),
ptSeedThr(iConfig.getUntrackedParameter<double>("PtSeedThr")),
etaSeedThr(iConfig.getUntrackedParameter<double>("EtaSeedThr")),
Expand All @@ -34,94 +33,31 @@ PythiaFilterGammaGamma::PythiaFilterGammaGamma(const edm::ParameterSet& iConfig)
pTMinCandidate1(iConfig.getUntrackedParameter<double>("PtMinCandidate1")),
pTMinCandidate2(iConfig.getUntrackedParameter<double>("PtMinCandidate2")),
etaMaxCandidate(iConfig.getUntrackedParameter<double>("EtaMaxCandidate")),
invMassWide(iConfig.getUntrackedParameter<double>("InvMassWide")),
invMassNarrow(iConfig.getUntrackedParameter<double>("InvMassNarrow")),
invMassMin(iConfig.getUntrackedParameter<double>("InvMassMin")),
invMassMax(iConfig.getUntrackedParameter<double>("InvMassMax")),
energyCut(iConfig.getUntrackedParameter<double>("EnergyCut")),
nTkConeMax(iConfig.getUntrackedParameter<int>("NTkConeMax")),
nTkConeSum(iConfig.getUntrackedParameter<int>("NTkConeSum")),
acceptPrompts(iConfig.getUntrackedParameter<bool>("AcceptPrompts")),
promptPtThreshold(iConfig.getUntrackedParameter<double>("PromptPtThreshold")) {

//cout<<"ptSeedThr "<<ptSeedThr<<endl;
//cout<<"etaSeedThr "<<etaSeedThr<<endl;
//cout<<"ptGammaThr "<<ptGammaThr<<endl;
//cout<<"etaGammaThr "<<etaGammaThr<<endl;
//cout<<"ptTkThr "<<ptTkThr<<endl;
//cout<<"etaTkThr "<<etaTkThr<<endl;
//cout<<"ptElThr "<<ptElThr<<endl;
//cout<<"etaElThr "<<etaElThr<<endl;
//cout<<"dRTkMax "<<dRTkMax<<endl;
//cout<<"dRSeedMax "<<dRSeedMax<<endl;
//cout<<"dPhiSeedMax "<<dPhiSeedMax<<endl;
//cout<<"dEtaSeedMax "<<dEtaSeedMax<<endl;
//cout<<"dRNarrowCone "<<dRNarrowCone<<endl;
//cout<<"pTMinCandidate1 "<<pTMinCandidate1<<endl;
//cout<<"pTMinCandidate2 "<<pTMinCandidate2<<endl;
//cout<<"etaMaxCandidate "<<etaMaxCandidate<<endl;
//cout<<"invMassWide "<<invMassWide<<endl;
//cout<<"invMassNarrow "<<invMassNarrow<<endl;
//cout<<"nTkConeMax "<<nTkConeMax<<endl;
//cout<<"nTkConeSum "<<nTkConeSum<<endl;
//cout<<"acceptPrompts "<<acceptPrompts<<endl;
//cout<<"promptPtThreshold "<<promptPtThreshold<<endl;

nSelectedEvents = 0;
nGeneratedEvents = 0;

//char a[20];
//for(int i=0; i<2; i++) {
// sprintf(a, "PT Seed %d", i+1);
// hPtSeed[i] = new TH1D(a, a, 100, 0, 200);
// sprintf(a, "Eta Seed %d", i+1);
// hEtaSeed[i]= new TH1D(a, a, 100, -3., 3.);
// sprintf(a, "Pid Seed %d", i+1);
// hPidSeed[i]= new TH1I(a, a, 50, 0, 500);
// sprintf(a, "PT Candidate %d", i+1);
// hPtCandidate[i] = new TH1D(a, a, 100, 0, 200);
// sprintf(a, "Eta Candidate %d", i+1);
// hEtaCandidate[i]= new TH1D(a, a, 100, -3., 3.);
// sprintf(a, "Pid Candidate %d", i+1);
// hPidCandidate[i]= new TH1I(a, a, 50, 0, 500);
// sprintf(a, "NTk Iso %d", i+1);
// hNTk[i]= new TH1I(a, a, 50, 0, 50);
//}
//hMassNarrow = new TH1D("Mass narr.", "Mass narr.", 100, 0, 1000);
//hMassWide= new TH1D("Mass wide", "Mass wide", 100, 0, 1000);
//hNTkSum= new TH1I("NTk Sum Iso", "NTk Sum Iso", 50, 0, 50);
counterPrompt = 0;

}

PythiaFilterGammaGamma::~PythiaFilterGammaGamma()
{
//writeFile();
cout << "Number of Selected Events: " << nSelectedEvents << endl;
cout << "Number of Generated Events: " << nGeneratedEvents << endl;
cout << "Number of Prompt photons: " << counterPrompt << endl;
}

//void PythiaFilterGammaGamma::writeFile() {

//TFile* file = new TFile (fileName.c_str(), "recreate");

//for(int i=0; i<2; i++) {
// hPtSeed[i]->Write();
// hEtaSeed[i]->Write();
// hPidSeed[i]->Write();
// hPtCandidate[i]->Write();
// hEtaCandidate[i]->Write();
// hPidCandidate[i]->Write();
// hNTk[i]->Write();
//}
//hMassNarrow->Write();
//hMassWide->Write();
//hNTkSum->Write();

//file->Close();

//}

bool PythiaFilterGammaGamma::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {


if(nSelectedEvents >= maxEvents) {
//writeFile();
throw cms::Exception("endJob")<<"We have reached the maximum number of events...";
}

Expand All @@ -136,49 +72,22 @@ bool PythiaFilterGammaGamma::filter(edm::Event& iEvent, const edm::EventSetup& i
std::vector<const GenParticle*> seeds, egamma, stable;
std::vector<const GenParticle*>::const_iterator itPart, itStable, itEn;

// Loop on egamma
for(HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {

if (
((*p)->status()==1&&(*p)->pdg_id() == 22) || // gamma
((*p)->status()==1&&abs((*p)->pdg_id()) == 11) || // electron
(*p)->pdg_id() == 111 || // pi0
abs((*p)->pdg_id()) == 221 || // eta
abs((*p)->pdg_id()) == 331 || // eta prime
abs((*p)->pdg_id()) == 113 || // rho0
abs((*p)->pdg_id()) == 223) // omega
((*p)->status()==1&&abs((*p)->pdg_id()) == 11)) // electron

{
// check for eta and pT threshold for seed in gamma, el
if ((*p)->momentum().perp() > ptSeedThr &&
fabs((*p)->momentum().eta()) < etaSeedThr) {

// check if found is daughter of one already taken
bool isUsed = false;

const GenParticle* mother = (*p)->production_vertex() ?
*((*p)->production_vertex()->particles_in_const_begin()) : 0;
const GenParticle* motherMother = (mother != 0 && mother->production_vertex()) ?
*(mother->production_vertex()->particles_in_const_begin()) : 0;
const GenParticle* motherMotherMother = (motherMother != 0 && motherMother->production_vertex()) ?
*(motherMother->production_vertex()->particles_in_const_begin()) : 0;

for(itPart = seeds.begin(); itPart != seeds.end(); itPart++) {

if ((*itPart) == mother ||
(*itPart) == motherMother ||
(*itPart) == motherMotherMother) {
isUsed = true;
break;
}
}

if (!isUsed) seeds.push_back(*p);

seeds.push_back(*p);
}
}
}

if (seeds.size() < 2) return accepted;
}

for(HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {

if ((*p)->status() == 1) {

Expand Down Expand Up @@ -209,10 +118,16 @@ bool PythiaFilterGammaGamma::filter(edm::Event& iEvent, const edm::EventSetup& i
}
}

if (seeds.size() < 2) return accepted;

std::vector<int> nTracks;
std::vector<TLorentzVector> candidate, candidateNarrow, candidateSeed;
std::vector<const GenParticle*>::iterator itSeed;

const GenParticle* mom;
int this_id;
int first_different_id;

for(itSeed = seeds.begin(); itSeed != seeds.end(); ++itSeed) {

TLorentzVector energy, narrowCone, temp1, temp2, tempseed;
Expand All @@ -236,7 +151,6 @@ bool PythiaFilterGammaGamma::filter(edm::Event& iEvent, const edm::EventSetup& i
}

int counter = 0;
//bool isIsolated = true;

if ( energy.Et() != 0. ) {
if (fabs(energy.Eta()) < etaMaxCandidate) {
Expand All @@ -250,57 +164,54 @@ bool PythiaFilterGammaGamma::filter(edm::Event& iEvent, const edm::EventSetup& i
}

if(acceptPrompts) {
bool isPrompt=false;
if((*itSeed)->status() == 1&&(*itSeed)->pdg_id() == 22) {
const GenParticle* mother = (*itSeed)->production_vertex() ?
*((*itSeed)->production_vertex()->particles_in_const_begin()) : 0;
if(mother) {
if(mother->pdg_id()>=-22&&mother->pdg_id()<=22) {
const GenParticle* motherMother = (mother != 0 && mother->production_vertex()) ?
*(mother->production_vertex()->particles_in_const_begin()) : 0;
if(motherMother) {
if(motherMother->pdg_id()>=-22&&motherMother->pdg_id()<=22) {
if((*itSeed)->momentum().perp()>promptPtThreshold) {
isPrompt=true;
}
}
}

if ((*itSeed)->momentum().perp()>promptPtThreshold)
{
bool isPrompt=true;
this_id = (*itSeed)->pdg_id();
mom = (*itSeed);
while (mom->pdg_id() == this_id) {
const GenParticle* mother = mom->production_vertex() ?
*(mom->production_vertex()->particles_in_const_begin()) : 0;

mom = mother;
if (mom == 0) {
break;
}
}

first_different_id = mom->pdg_id();

if (mom->status() == 2 && std::abs(first_different_id)>100) isPrompt=false;


if(isPrompt) counter=0;
if(isPrompt) counterPrompt++;
}
}
if(isPrompt) counter=0;
}
// check number of tracks
//if (counter <= nTkConeMax) isIsolated = true;
}
}

// check pt candidate, check nTrack, check eta
//if (isIsolated) {
candidate.push_back(energy);
candidateSeed.push_back(tempseed);
candidateNarrow.push_back(narrowCone);
nTracks.push_back(counter);
//++itSeed;
//}
}
}

if (candidate.size() <2) return accepted;

TLorentzVector minv, minvNarrow;


//bool filled=false;
TLorentzVector minvMin, minvMax;

int i1, i2;
for(unsigned int i=0; i<candidate.size()-1; ++i) {

if (candidate[i].Energy() <1.) continue;
if (candidate[i].Energy() < energyCut) continue;
if(nTracks[i]>nTkConeMax) continue;
if (fabs(candidate[i].Eta()) > etaMaxCandidate) continue;

for(unsigned int j=i+1; j<candidate.size(); ++j) {
if (candidate[j].Energy() <1.) continue;
if (candidate[j].Energy() < energyCut) continue;
if(nTracks[j]>nTkConeMax) continue;
if (fabs(candidate[j].Eta()) > etaMaxCandidate) continue;

Expand All @@ -317,34 +228,14 @@ bool PythiaFilterGammaGamma::filter(edm::Event& iEvent, const edm::EventSetup& i

if (candidate[i1].Pt() < pTMinCandidate1 || candidate[i2].Pt() < pTMinCandidate2) continue;

minv = candidate[i] + candidate[j];
if (minv.M() < invMassWide) continue;
minvMin = candidate[i] + candidate[j];
if (minvMin.M() < invMassMin) continue;

minvNarrow = candidateNarrow[i] + candidateNarrow[j];
if (minvNarrow.M() > invMassNarrow) continue;
minvMax = candidate[i] + candidate[j];
if (minvMax.M() > invMassMax) continue;

accepted = true;

//if(!filled) {
//hMassWide->Fill(minv.M());
//hMassNarrow->Fill(minvNarrow.M());
//hNTkSum->Fill(nTracks[i] + nTracks[j]);
//hPtCandidate[0]->Fill(candidate[i1].Pt());
//hPtCandidate[1]->Fill(candidate[i2].Pt());
//hEtaCandidate[0]->Fill(candidate[i1].Eta());
//hEtaCandidate[1]->Fill(candidate[i2].Eta());
//hPidCandidate[0]->Fill(seeds[i1]->pdg_id());
//hPidCandidate[1]->Fill(seeds[i2]->pdg_id());
//hPtSeed[0]->Fill(candidateSeed[i1].Pt());
//hPtSeed[1]->Fill(candidateSeed[i2].Pt());
//hEtaSeed[0]->Fill(candidateSeed[i1].Eta());
//hEtaSeed[1]->Fill(candidateSeed[i2].Eta());
//hPidSeed[0]->Fill(seeds[i1]->pdg_id());
//hPidSeed[1]->Fill(seeds[i2]->pdg_id());
//hNTk[0]->Fill(nTracks[i1]);
//hNTk[1]->Fill(nTracks[i2]);
//filled=true;
//}
}
}

Expand Down

0 comments on commit e39617c

Please sign in to comment.