Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

made EM enriched filter more indipendent from the generator used and add... #6976

Merged
merged 4 commits into from Jan 25, 2015
Merged

made EM enriched filter more indipendent from the generator used and add... #6976

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7e5a227
made EM enriched filter more indipendent from the generator used and …
mmachet Dec 18, 2014
71ca6b0
loop on neutral mesons removed and prompt photons definition changed
mmachet Jan 19, 2015
2d38fa7
changed the code structure for prompt photons definition
mmachet Jan 22, 2015
e4ebcf7
fixed compilation error
mmachet Jan 24, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
13 changes: 3 additions & 10 deletions GeneratorInterface/GenFilters/interface/PythiaFilterGammaGamma.h
View file
Open in desktop
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
Open in desktop
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