Skip to content

Commit 3d15312

Browse files
author
Michele Mormile
committed
Moved to CMSSW_14
1 parent 27bec5f commit 3d15312

File tree

3 files changed

+221
-0
lines changed

3 files changed

+221
-0
lines changed

‎GeneratorInterface/GenFilters/plugins/PhotonGenFilter.cc‎

Lines changed: 169 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,169 @@
1+
#include "GeneratorInterface/GenFilters/plugins/PhotonGenFilter.h"
2+
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
3+
#include <iostream>
4+
5+
using namespace edm;
6+
using namespace std;
7+
8+
PhotonGenFilter::PhotonGenFilter(const edm::ParameterSet &iConfig)
9+
: token_(consumes<edm::HepMCProduct>(
10+
edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))) {
11+
// Constructor implementation
12+
ptMin = iConfig.getUntrackedParameter<double>("MinPt", 20.);
13+
etaMin = iConfig.getUntrackedParameter<double>("MinEta", -2.4);
14+
etaMax = iConfig.getUntrackedParameter<double>("MaxEta", 2.4);
15+
drMin = iConfig.getUntrackedParameter<double>("drMin", 0.1);
16+
ptThreshold = iConfig.getUntrackedParameter<double>("ptThreshold", 2.);
17+
}
18+
19+
PhotonGenFilter::~PhotonGenFilter() {
20+
// Destructor implementation
21+
}
22+
23+
bool PhotonGenFilter::filter(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
24+
using namespace edm;
25+
Handle<HepMCProduct> evt;
26+
iEvent.getByToken(token_, evt);
27+
const HepMC::GenEvent *myGenEvent = evt->GetEvent();
28+
bool accepted_event = false;
29+
30+
for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
31+
++p) { // loop through all particles
32+
if ((*p)->pdg_id() == 22 && !hasAncestor(*p, [](int x) {
33+
return x > 30 && x != 2212;
34+
})) { // loop through all photons that don't come from hadrons
35+
if ((*p)->momentum().perp() > ptMin && (*p)->status() == 1 && (*p)->momentum().eta() > etaMin &&
36+
(*p)->momentum().eta() < etaMax) { // check if photon passes pt and eta cuts
37+
bool good_photon = true;
38+
double phi = (*p)->momentum().phi();
39+
double eta = (*p)->momentum().eta();
40+
double pt = (*p)->momentum().perp();
41+
double frixione_isolation_coefficient = pt / (1 - cos(drMin));
42+
vector<double> particles_pt, particles_deltar;
43+
for (HepMC::GenEvent::particle_const_iterator q = myGenEvent->particles_begin();
44+
q != myGenEvent->particles_end();
45+
++q) { // loop through all particles to compute frixione isolation
46+
if (&p != &q) { // don't compare the photon to itself
47+
if ((*q)->momentum().perp() > ptThreshold && (*q)->pdg_id() != 22 && (*q)->pdg_id() != 12 &&
48+
(*q)->pdg_id() != 14 && (*q)->pdg_id() != 16 &&
49+
(*q)->status() == 1) { // check if particle passes pt and status cuts and is not a neutrino or photon
50+
double phi2 = (*q)->momentum().phi();
51+
double deltaphi = fabs(phi - phi2);
52+
if (deltaphi > M_PI)
53+
deltaphi = 2. * M_PI - deltaphi;
54+
double eta2 = (*q)->momentum().eta();
55+
double deltaeta = fabs(eta - eta2);
56+
double deltaR = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
57+
if (deltaR < drMin) // check if particle is within drMin of photon, for isolation
58+
{
59+
particles_pt.push_back((*q)->momentum().perp());
60+
particles_deltar.push_back(deltaR);
61+
}
62+
}
63+
}
64+
}
65+
//order particles_pt and particles_deltar according to increasing particles_deltar
66+
for (unsigned int i = 0; i < particles_deltar.size(); i++) {
67+
for (unsigned int j = i + 1; j < particles_deltar.size(); j++) {
68+
if (particles_deltar[i] > particles_deltar[j]) {
69+
double temp = particles_deltar[i];
70+
particles_deltar[i] = particles_deltar[j];
71+
particles_deltar[j] = temp;
72+
temp = particles_pt[i];
73+
particles_pt[i] = particles_pt[j];
74+
particles_pt[j] = temp;
75+
}
76+
}
77+
}
78+
//calculate frixione isolation
79+
double total_pt = 0;
80+
for (unsigned int i = 0; i < particles_deltar.size(); i++) {
81+
total_pt += particles_pt[i];
82+
if (total_pt > frixione_isolation_coefficient * (1 - cos(particles_deltar[i]))) {
83+
good_photon =
84+
false; // if for some delta R the isolation condition is not satisfied, the photon is not good
85+
break;
86+
}
87+
}
88+
if (good_photon) {
89+
if (hasAncestor(*p, [](int x) { return x == 11 || x == 13 || x == 15; }) ||
90+
hasAncestor(
91+
*p,
92+
[](int x) { return x == 24 || x == 5; },
93+
true,
94+
false)) { // check if photon is from decay and defines the event as acceptable
95+
accepted_event = true;
96+
}
97+
if (!hasAncestor(*p, [](int x) { return x == 11 || x == 13 || x == 15; }) &&
98+
hasAncestor(
99+
*p,
100+
[](int x) { return x == 24 || x == 5; },
101+
false,
102+
true)) { // check if photon comes from the hard process and discards it, in case it is
103+
return false;
104+
}
105+
}
106+
}
107+
}
108+
}
109+
110+
// Implementation for event filtering
111+
return accepted_event; // Accept if it has found at least one good photon from decay and none from hard process
112+
}
113+
114+
void PhotonGenFilter::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
115+
edm::ParameterSetDescription desc;
116+
desc.addUntracked<double>("MaxEta", 2.4);
117+
desc.addUntracked<double>("MinEta", -2.4);
118+
desc.addUntracked<double>("MinPt", 20.);
119+
desc.addUntracked<double>("drMin", 0.1);
120+
desc.addUntracked<double>("ptThreshold", 2.);
121+
122+
descriptions.add("PhotonGenFilter", desc);
123+
}
124+
125+
bool PhotonGenFilter::hasAncestor(
126+
HepMC::GenParticle *particle,
127+
function<bool(int)> check,
128+
bool isWorBFromDecayCheck,
129+
bool isWorBPromptCheck) const { // function to check if a particle has a certain ancestor
130+
if (!particle)
131+
return false;
132+
133+
HepMC::GenVertex *prodVertex = particle->production_vertex();
134+
135+
// If the particle doesn't have a production vertex, it has no parents.
136+
if (!prodVertex)
137+
return false;
138+
139+
// Loop over all parents (incoming particles) of the vertex
140+
for (auto parent = prodVertex->particles_begin(HepMC::parents); parent != prodVertex->particles_end(HepMC::parents);
141+
++parent) {
142+
int pdgId = abs((*parent)->pdg_id());
143+
144+
// Check if the PDG ID respects a check
145+
if (check(pdgId)) {
146+
if (isWorBFromDecayCheck) {
147+
return hasAncestor(
148+
*parent,
149+
[](int x) { return x == 6; },
150+
false,
151+
false); // if the photon has a W or b quark ancestor, check if it comes from a top quark
152+
} else if (isWorBPromptCheck) {
153+
return !hasAncestor(
154+
*parent,
155+
[](int x) { return x == 6; },
156+
false,
157+
false); // if the photon has a W or b quark ancestor, check that it doesn't come from a top quark decay (therefore it comes form the hard process)
158+
} else
159+
return true;
160+
}
161+
162+
return hasAncestor(
163+
*parent, check, isWorBFromDecayCheck, isWorBPromptCheck); // Recursively check the ancestry of the parent
164+
}
165+
166+
return false;
167+
}
168+
169+
DEFINE_FWK_MODULE(PhotonGenFilter);

‎GeneratorInterface/GenFilters/plugins/PhotonGenFilter.h‎

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
#ifndef PhotonGenFilter_h
2+
#define PhotonGenFilter_h
3+
4+
#include <memory>
5+
#include <iostream>
6+
7+
#include "FWCore/Framework/interface/Frameworkfwd.h"
8+
#include "FWCore/Framework/interface/global/EDFilter.h"
9+
10+
#include "FWCore/Framework/interface/Event.h"
11+
#include "FWCore/Framework/interface/MakerMacros.h"
12+
13+
#include "FWCore/ParameterSet/interface/ParameterSet.h"
14+
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
15+
16+
namespace edm {
17+
class HepMCProduct;
18+
}
19+
20+
class PhotonGenFilter : public edm::global::EDFilter<> {
21+
public:
22+
explicit PhotonGenFilter(const edm::ParameterSet&);
23+
~PhotonGenFilter() override;
24+
25+
bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
26+
bool hasAncestor(HepMC::GenParticle* particle,
27+
std::function<bool(int)> check,
28+
bool isWorBFromDecayCheck = false,
29+
bool isWorBPromptCheck = false) const;
30+
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
31+
32+
private:
33+
// Private member variables and functions
34+
35+
const edm::EDGetTokenT<edm::HepMCProduct> token_;
36+
double ptMin;
37+
double etaMin;
38+
double etaMax;
39+
double drMin;
40+
double ptThreshold;
41+
};
42+
43+
#endif

‎GeneratorInterface/GenFilters/python/PhotonGenFilter_cfi.py‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
import FWCore.ParameterSet.Config as cms
2+
3+
photonGenFilter = cms.EDFilter('PhotonGenFilter',
4+
MaxEta = cms.untracked.double(2.4),
5+
MinEta = cms.untracked.double(-2.4),
6+
MinPt = cms.untracked.double(20.),
7+
drMin = cms.untracked.double(0.1),
8+
ptThreshold = cms.untracked.double(2.)
9+
)

0 commit comments

Comments
 (0)