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.
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
double EM enrichment filter support for HepMCfilter #22980
Changes from 16 commits
4e442ba
920bd47
87cc06c
ece1cb0
8e59d8d
61c5490
5ffcaa4
90b2fc3
501e132
c7cef8b
df3514f
0b3182c
ebcdf80
08210c8
a934bbb
c579a4d
05c3cd6
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,40 @@ | ||
import FWCore.ParameterSet.Config as cms | ||
|
||
doubleEMenrichingHepMCfilterParams = cms.PSet( | ||
# seed thresholds | ||
PtSeedThr = cms.untracked.double(5.0), | ||
EtaSeedThr = cms.untracked.double(2.8), | ||
|
||
# photon thresholds | ||
PtGammaThr = cms.untracked.double(0.0), | ||
EtaGammaThr = cms.untracked.double(2.8), | ||
|
||
# electron threshold | ||
PtElThr = cms.untracked.double(2.0), | ||
EtaElThr = cms.untracked.double(2.8), | ||
|
||
dRSeedMax = cms.untracked.double(0.0), | ||
dPhiSeedMax = cms.untracked.double(0.2), | ||
dEtaSeedMax = cms.untracked.double(0.12), | ||
dRNarrowCone = cms.untracked.double(0.02), | ||
|
||
PtTkThr = cms.untracked.double(1.6), | ||
EtaTkThr = cms.untracked.double(2.2), | ||
dRTkMax = cms.untracked.double(0.2), | ||
|
||
PtMinCandidate1 = cms.untracked.double(15.0), | ||
PtMinCandidate2 = cms.untracked.double(15.0), | ||
EtaMaxCandidate = cms.untracked.double(3.0), | ||
|
||
NTkConeMax = cms.untracked.int32(2), | ||
NTkConeSum = cms.untracked.int32(4), | ||
|
||
# mass 80..Inf GeV | ||
InvMassMin = cms.untracked.double(80.0), | ||
InvMassMax = cms.untracked.double(14000.0), | ||
|
||
EnergyCut = cms.untracked.double(1.0), | ||
|
||
AcceptPrompts = cms.untracked.bool(True), | ||
PromptPtThreshold = cms.untracked.double(15.0), | ||
) |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,49 @@ | ||
#ifndef PYTHIAHEPMCFILTERGAMMAGAMMA_h | ||
#define PYTHIAHEPMCFILTERGAMMAGAMMA_h | ||
|
||
// | ||
// Package: GeneratorInterface/GenFilters | ||
// Class: PythiaHepMCFilterGammaGamma | ||
// | ||
// Original Author: Matteo Sani | ||
// | ||
// | ||
|
||
#include "FWCore/Framework/interface/Frameworkfwd.h" | ||
#include "GeneratorInterface/Core/interface/BaseHepMCFilter.h" | ||
#include "FWCore/Framework/interface/EDFilter.h" | ||
|
||
#include "FWCore/Framework/interface/Event.h" | ||
#include "FWCore/Framework/interface/MakerMacros.h" | ||
#include "FWCore/ParameterSet/interface/ParameterSet.h" | ||
|
||
namespace edm { | ||
class HepMCProduct; | ||
} | ||
|
||
#include "TH1D.h" | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. where is this used? In any case a filter does not seem to me the right place to fill histograms, even though this is present in a couple of them in GenFilters... There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. thanks for spotting this, this was a leftover from the original PythiaFilterGammaGamma.h in CMSSW_9_3_1. We removed this in 05c3cd6 . |
||
#include "TH1I.h" | ||
|
||
class PythiaHepMCFilterGammaGamma : public BaseHepMCFilter { | ||
public: | ||
explicit PythiaHepMCFilterGammaGamma(const edm::ParameterSet&); | ||
~PythiaHepMCFilterGammaGamma() override; | ||
|
||
/** @return true if this GenEvent passes the double EM enrichment | ||
criterion */ | ||
bool filter(const HepMC::GenEvent* myGenEvent) override; | ||
private: | ||
|
||
const edm::EDGetTokenT<edm::HepMCProduct> token_; | ||
const int maxEvents; | ||
|
||
const double ptSeedThr, etaSeedThr, ptGammaThr, etaGammaThr, ptTkThr, etaTkThr; | ||
const double ptElThr, etaElThr, dRTkMax, dRSeedMax, dPhiSeedMax, dEtaSeedMax, dRNarrowCone, pTMinCandidate1, pTMinCandidate2, etaMaxCandidate; | ||
const double invMassMin, invMassMax; | ||
const double energyCut; | ||
const int nTkConeMax, nTkConeSum; | ||
const bool acceptPrompts; | ||
const double promptPtThreshold; | ||
|
||
}; | ||
#endif |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,225 @@ | ||
#include "GeneratorInterface/Core/interface/PythiaHepMCFilterGammaGamma.h" | ||
#include "DataFormats/Math/interface/LorentzVector.h" | ||
#include "Math/GenVector/VectorUtil.h" | ||
//#include "CLHEP/HepMC/GenParticle.h" | ||
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" | ||
|
||
#include "TLorentzVector.h" | ||
|
||
#include <iostream> | ||
|
||
using namespace edm; | ||
using namespace std; | ||
using namespace HepMC; | ||
|
||
|
||
PythiaHepMCFilterGammaGamma::PythiaHepMCFilterGammaGamma(const edm::ParameterSet& iConfig) : | ||
maxEvents(iConfig.getUntrackedParameter<int>("maxEvents", 0)), | ||
ptSeedThr(iConfig.getUntrackedParameter<double>("PtSeedThr")), | ||
etaSeedThr(iConfig.getUntrackedParameter<double>("EtaSeedThr")), | ||
ptGammaThr(iConfig.getUntrackedParameter<double>("PtGammaThr")), | ||
etaGammaThr(iConfig.getUntrackedParameter<double>("EtaGammaThr")), | ||
ptTkThr(iConfig.getUntrackedParameter<double>("PtTkThr")), | ||
etaTkThr(iConfig.getUntrackedParameter<double>("EtaTkThr")), | ||
ptElThr(iConfig.getUntrackedParameter<double>("PtElThr")), | ||
etaElThr(iConfig.getUntrackedParameter<double>("EtaElThr")), | ||
dRTkMax(iConfig.getUntrackedParameter<double>("dRTkMax")), | ||
dRSeedMax(iConfig.getUntrackedParameter<double>("dRSeedMax")), | ||
dPhiSeedMax(iConfig.getUntrackedParameter<double>("dPhiSeedMax")), | ||
dEtaSeedMax(iConfig.getUntrackedParameter<double>("dEtaSeedMax")), | ||
dRNarrowCone(iConfig.getUntrackedParameter<double>("dRNarrowCone")), | ||
pTMinCandidate1(iConfig.getUntrackedParameter<double>("PtMinCandidate1")), | ||
pTMinCandidate2(iConfig.getUntrackedParameter<double>("PtMinCandidate2")), | ||
etaMaxCandidate(iConfig.getUntrackedParameter<double>("EtaMaxCandidate")), | ||
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")) { | ||
|
||
if (maxEvents != 0) edm::LogInfo("PythiaFilterGammaGamma::PythiaFilterGammaGamma") << "WARNING, ignoring unsuported option, maxEvents = " << maxEvents << endl; | ||
|
||
} | ||
|
||
PythiaHepMCFilterGammaGamma::~PythiaHepMCFilterGammaGamma() | ||
{ | ||
} | ||
|
||
bool PythiaHepMCFilterGammaGamma::filter(const HepMC::GenEvent* myGenEvent) { | ||
|
||
bool accepted = false; | ||
|
||
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 | ||
|
||
{ | ||
// check for eta and pT threshold for seed in gamma, el | ||
if ((*p)->momentum().perp() > ptSeedThr && | ||
fabs((*p)->momentum().eta()) < etaSeedThr) { | ||
|
||
seeds.push_back(*p); | ||
} | ||
} | ||
|
||
|
||
if ((*p)->status() == 1) { | ||
|
||
// save charged stable tracks | ||
if (abs((*p)->pdg_id()) == 211 || | ||
abs((*p)->pdg_id()) == 321 || | ||
abs((*p)->pdg_id()) == 11 || | ||
abs((*p)->pdg_id()) == 13 || | ||
abs((*p)->pdg_id()) == 15) { | ||
// check if it passes the cut | ||
if ((*p)->momentum().perp() > ptTkThr && | ||
fabs((*p)->momentum().eta()) < etaTkThr) { | ||
stable.push_back(*p); | ||
} | ||
} | ||
|
||
// save egamma for candidate calculation | ||
if ((*p)->pdg_id() == 22 && | ||
(*p)->momentum().perp() > ptGammaThr && | ||
fabs((*p)->momentum().eta()) < etaGammaThr) { | ||
egamma.push_back(*p); | ||
} | ||
if (abs((*p)->pdg_id()) == 11 && | ||
(*p)->momentum().perp() > ptElThr && | ||
fabs((*p)->momentum().eta()) < etaElThr) { | ||
egamma.push_back(*p); | ||
} | ||
} | ||
} | ||
|
||
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; | ||
|
||
tempseed.SetXYZM((*itSeed)->momentum().px(), (*itSeed)->momentum().py(), (*itSeed)->momentum().pz(), 0); | ||
for(itEn = egamma.begin(); itEn != egamma.end(); ++itEn) { | ||
temp1.SetXYZM((*itEn)->momentum().px(), (*itEn)->momentum().py(), (*itEn)->momentum().pz(), 0); | ||
|
||
double DR = temp1.DeltaR(tempseed); | ||
double DPhi = temp1.DeltaPhi(tempseed); | ||
double DEta = (*itEn)->momentum().eta()-(*itSeed)->momentum().eta(); | ||
if(DPhi<0) DPhi=-DPhi; | ||
if(DEta<0) DEta=-DEta; | ||
|
||
if (DR < dRSeedMax || (DPhi<dPhiSeedMax&&DEta<dEtaSeedMax)) { | ||
energy += temp1; | ||
} | ||
if (DR < dRNarrowCone) { | ||
narrowCone += temp1; | ||
} | ||
} | ||
|
||
int counter = 0; | ||
|
||
if ( energy.Et() != 0. ) { | ||
if (fabs(energy.Eta()) < etaMaxCandidate) { | ||
|
||
temp2.SetXYZM(energy.Px(), energy.Py(), energy.Pz(), 0); | ||
|
||
for(itStable = stable.begin(); itStable != stable.end(); ++itStable) { | ||
temp1.SetXYZM((*itStable)->momentum().px(), (*itStable)->momentum().py(), (*itStable)->momentum().pz(), 0); | ||
double DR = temp1.DeltaR(temp2); | ||
if (DR < dRTkMax) counter++; | ||
} | ||
|
||
if(acceptPrompts) { | ||
|
||
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()) : nullptr; | ||
|
||
mom = mother; | ||
if (mom == nullptr) { | ||
break; | ||
} | ||
} | ||
|
||
first_different_id = mom->pdg_id(); | ||
|
||
if (mom->status() == 2 && std::abs(first_different_id)>100) isPrompt=false; | ||
|
||
|
||
if(isPrompt) counter=0; | ||
} | ||
} | ||
} | ||
} | ||
|
||
candidate.push_back(energy); | ||
candidateSeed.push_back(tempseed); | ||
candidateNarrow.push_back(narrowCone); | ||
nTracks.push_back(counter); | ||
} | ||
|
||
if (candidate.size() <2) return accepted; | ||
|
||
TLorentzVector minvMin, minvMax; | ||
|
||
int i1, i2; | ||
for(unsigned int i=0; i<candidate.size()-1; ++i) { | ||
|
||
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() < energyCut) continue; | ||
if(nTracks[j]>nTkConeMax) continue; | ||
if (fabs(candidate[j].Eta()) > etaMaxCandidate) continue; | ||
|
||
if (nTracks[i] + nTracks[j] > nTkConeSum) continue; | ||
|
||
if (candidate[i].Pt() > candidate[j].Pt()) { | ||
i1 = i; | ||
i2 = j; | ||
} | ||
else { | ||
i1 = j; | ||
i2 = i; | ||
} | ||
|
||
if (candidate[i1].Pt() < pTMinCandidate1 || candidate[i2].Pt() < pTMinCandidate2) continue; | ||
|
||
minvMin = candidate[i] + candidate[j]; | ||
if (minvMin.M() < invMassMin) continue; | ||
|
||
minvMax = candidate[i] + candidate[j]; | ||
if (minvMax.M() > invMassMax) continue; | ||
|
||
accepted = true; | ||
|
||
} | ||
} | ||
|
||
return accepted; | ||
} | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
why adding a dependence of the GeneratorInterface to root??
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This was carried over from the original
GeneratorInterface/GenFilters/src/PythiaFilterGammaGamma.cc
which usesTLorentzVector
.Should we replace this with a different class instead ?
We moved the code in
PythiaFilterGammaGamma.cc
to a new class inGeneratorInterface/Core
because the code needs to be used both from the originalGeneratorInterface/GenFilters/src/PythiaFilterGammaGamma.cc
(for backwards compatibility) and also fromGeneratorInterface/Core/src/HepMCFilterDriver.cc
.We could therefore not leave it in
GeneratorInterface/GenFilters
because we would have to add a dependency ofGeneratorInterface/Core
onGeneratorInterface/GenFilters
which would create a circular dependency.