|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
|
#include "GeneratorInterface/GenFilters/plugins/PythiaProbeFilter.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include <iostream>
#include <memory>
using namespace edm;
using namespace std;
using namespace Pythia8;
PythiaProbeFilter::PythiaProbeFilter(const edm::ParameterSet& iConfig)
: token_(consumes<edm::HepMCProduct>(
edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
MomID(iConfig.getUntrackedParameter("MomID", 0)),
GrandMomID(iConfig.getUntrackedParameter("GrandMomID", 0)),
chargeconju(iConfig.getUntrackedParameter("ChargeConjugation", true)),
nsisters(iConfig.getUntrackedParameter("NumberOfSisters", 0)),
naunts(iConfig.getUntrackedParameter("NumberOfAunts", 0)),
minptcut(iConfig.getUntrackedParameter("MinPt", 0.)),
maxptcut(iConfig.getUntrackedParameter("MaxPt", 14000.)),
minetacut(iConfig.getUntrackedParameter("MinEta", -10.)),
maxetacut(iConfig.getUntrackedParameter("MaxEta", 10.)),
countQEDCorPhotons(iConfig.getUntrackedParameter("countQEDCorPhotons", false)) {
//now do what ever initialization is needed
vector<int> defID;
defID.push_back(0);
exclsisIDs = iConfig.getUntrackedParameter<vector<int> >("SisterIDs", defID);
exclauntIDs = iConfig.getUntrackedParameter<vector<int> >("AuntIDs", defID);
identicalParticle = false;
for (unsigned int ilist = 0; ilist < exclsisIDs.size(); ++ilist) {
if (fabs(exclsisIDs[ilist]) == fabs(particleID))
identicalParticle = true;
}
// create pythia8 instance to access particle data
edm::LogInfo("PythiaProbeFilter::PythiaProbeFilter") << "Creating pythia8 instance for particle properties" << endl;
if (!fLookupGen.get())
fLookupGen = std::make_unique<Pythia>();
}
PythiaProbeFilter::~PythiaProbeFilter() {
// do anything here that needs to be done at desctruction time
// (e.g. close files, deallocate resources etc.)
}
bool PythiaProbeFilter::AlreadyExcludedCheck(std::vector<unsigned int> excludedList, unsigned int current_part) const {
bool result = false;
for (unsigned int checkNow : excludedList) {
if (current_part != checkNow)
continue;
result = true;
break;
}
return result;
}
//
// member functions
//
// ------------ method called to produce the data ------------
bool PythiaProbeFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
using namespace edm;
bool accepted = false;
Handle<HepMCProduct> evt;
iEvent.getByToken(token_, evt);
const HepMC::GenEvent* myGenEvent = evt->GetEvent();
//access particles
for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
++p) {
//select tag particle
if (fabs((*p)->pdg_id()) == fabs(particleID))
if ((*p)->pdg_id() != particleID && !chargeconju)
continue;
if (fabs((*p)->pdg_id()) != fabs(particleID) && chargeconju)
continue;
//kinematic properties of tag
if ((*p)->momentum().perp() < minptcut || (*p)->momentum().perp() > maxptcut)
continue;
if ((*p)->momentum().eta() < minetacut || (*p)->momentum().eta() > maxetacut)
continue;
//remove probe side particles
bool excludeTagParticle = false;
if (naunts == 0) {
if ((*p)->production_vertex()) {
for (HepMC::GenVertex::particle_iterator anc = (*p)->production_vertex()->particles_begin(HepMC::parents);
anc != (*p)->production_vertex()->particles_end(HepMC::parents);
++anc) {
if (fabs((*anc)->pdg_id()) != fabs(MomID) && chargeconju)
continue;
else if ((*anc)->pdg_id() != MomID && !chargeconju)
continue;
int nsis = 0;
int exclsis = 0;
std::vector<unsigned int> checklistSis;
if ((*anc)->end_vertex()) {
for (HepMC::GenVertex::particle_iterator sis = (*anc)->end_vertex()->particles_begin(HepMC::children);
sis != (*anc)->end_vertex()->particles_end(HepMC::children);
++sis) {
//identify the tag particle in the decay
if ((*p)->pdg_id() == (*sis)->pdg_id() && (identicalParticle || !chargeconju))
continue;
if (fabs((*p)->pdg_id()) == fabs((*sis)->pdg_id()) && !identicalParticle && chargeconju)
continue;
//remove QED photons
if ((*sis)->pdg_id() == 22 && !countQEDCorPhotons)
continue;
nsis++;
//check if this sis is excluded already
for (unsigned int ilist = 0; ilist < exclsisIDs.size(); ++ilist) {
if (AlreadyExcludedCheck(checklistSis, ilist)) {
continue;
}
if (fabs(exclsisIDs[ilist]) == fabs((*sis)->pdg_id()) && chargeconju) {
exclsis++;
checklistSis.push_back(ilist);
}
if (exclsisIDs[ilist] == (*sis)->pdg_id() && !chargeconju) {
exclsis++;
checklistSis.push_back(ilist);
}
}
}
}
if (nsis == exclsis && nsis == nsisters) {
excludeTagParticle = true;
break;
}
}
}
} else if (naunts > 0) {
//now take into account that we have up 2 generations in the decay
if ((*p)->production_vertex()) {
for (HepMC::GenVertex::particle_iterator anc = (*p)->production_vertex()->particles_begin(HepMC::parents);
anc != (*p)->production_vertex()->particles_end(HepMC::parents);
++anc) {
if (fabs((*anc)->pdg_id()) != fabs(MomID) && chargeconju)
continue;
else if ((*anc)->pdg_id() != MomID && !chargeconju)
continue;
int nsis = 0;
int exclsis = 0;
std::vector<unsigned int> checklistSis;
int naunt = 0;
int exclaunt = 0;
std::vector<unsigned int> checklistAunt;
if ((*anc)->end_vertex()) {
for (HepMC::GenVertex::particle_iterator sis = (*anc)->end_vertex()->particles_begin(HepMC::children);
sis != (*anc)->end_vertex()->particles_end(HepMC::children);
++sis) {
//identify the particle under study in the decay
if ((*p)->pdg_id() == (*sis)->pdg_id() && (identicalParticle || !chargeconju))
continue;
if (fabs((*p)->pdg_id()) == fabs((*sis)->pdg_id()) && !identicalParticle && chargeconju)
continue;
//remove QED photons
if ((*sis)->pdg_id() == 22 && !countQEDCorPhotons)
continue;
nsis++;
for (unsigned int ilist = 0; ilist < exclsisIDs.size(); ++ilist) {
if (AlreadyExcludedCheck(checklistSis, ilist))
continue;
if (fabs(exclsisIDs[ilist]) == fabs((*sis)->pdg_id()) && chargeconju) {
exclsis++;
checklistSis.push_back(ilist);
}
if (exclsisIDs[ilist] == (*sis)->pdg_id() && !chargeconju) {
exclsis++;
checklistSis.push_back(ilist);
}
}
}
}
//check sisters
if (nsis != exclsis || nsis != nsisters)
break;
if ((*anc)->production_vertex()) {
for (HepMC::GenVertex::particle_iterator granc =
(*anc)->production_vertex()->particles_begin(HepMC::parents);
granc != (*anc)->production_vertex()->particles_end(HepMC::parents);
++granc) {
if (fabs((*granc)->pdg_id()) != fabs(GrandMomID) && chargeconju)
continue;
else if ((*granc)->pdg_id() != GrandMomID && !chargeconju)
continue;
for (HepMC::GenVertex::particle_iterator aunt = (*granc)->end_vertex()->particles_begin(HepMC::children);
aunt != (*granc)->end_vertex()->particles_end(HepMC::children);
++aunt) {
if ((*aunt)->pdg_id() == (*anc)->pdg_id())
continue;
if ((*aunt)->pdg_id() == 22 && !countQEDCorPhotons)
continue;
naunt++;
for (unsigned int ilist = 0; ilist < exclauntIDs.size(); ++ilist) {
if (AlreadyExcludedCheck(checklistAunt, ilist))
continue;
if (fabs(exclauntIDs[ilist]) == fabs((*aunt)->pdg_id()) && chargeconju) {
exclaunt++;
checklistAunt.push_back(ilist);
}
if (exclauntIDs[ilist] == (*aunt)->pdg_id() && !chargeconju) {
exclaunt++;
checklistAunt.push_back(ilist);
}
}
}
}
}
//check aunts
if (naunt == exclaunt && naunt == naunts) {
excludeTagParticle = true;
break;
}
}
}
}
if (excludeTagParticle)
continue;
accepted = true;
break;
}
if (accepted) {
return true;
} else {
return false;
}
}
|