3 #include "Math/GenVector/VectorUtil.h" 7 #include "TLorentzVector.h" 13 using namespace HepMC;
17 ptSeedThr(iConfig.getParameter<double>(
"PtSeedThr")),
18 etaSeedThr(iConfig.getParameter<double>(
"EtaSeedThr")),
19 ptGammaThr(iConfig.getParameter<double>(
"PtGammaThr")),
20 etaGammaThr(iConfig.getParameter<double>(
"EtaGammaThr")),
21 ptTkThr(iConfig.getParameter<double>(
"PtTkThr")),
22 etaTkThr(iConfig.getParameter<double>(
"EtaTkThr")),
23 ptElThr(iConfig.getParameter<double>(
"PtElThr")),
24 etaElThr(iConfig.getParameter<double>(
"EtaElThr")),
25 dRTkMax(iConfig.getParameter<double>(
"dRTkMax")),
26 dRSeedMax(iConfig.getParameter<double>(
"dRSeedMax")),
27 dPhiSeedMax(iConfig.getParameter<double>(
"dPhiSeedMax")),
28 dEtaSeedMax(iConfig.getParameter<double>(
"dEtaSeedMax")),
29 dRNarrowCone(iConfig.getParameter<double>(
"dRNarrowCone")),
30 pTMinCandidate1(iConfig.getParameter<double>(
"PtMinCandidate1")),
31 pTMinCandidate2(iConfig.getParameter<double>(
"PtMinCandidate2")),
32 etaMaxCandidate(iConfig.getParameter<double>(
"EtaMaxCandidate")),
33 invMassMin(iConfig.getParameter<double>(
"InvMassMin")),
34 invMassMax(iConfig.getParameter<double>(
"InvMassMax")),
35 energyCut(iConfig.getParameter<double>(
"EnergyCut")),
36 nTkConeMax(iConfig.getParameter<
int>(
"NTkConeMax")),
37 nTkConeSum(iConfig.getParameter<
int>(
"NTkConeSum")),
38 acceptPrompts(iConfig.getParameter<
bool>(
"AcceptPrompts")),
39 promptPtThreshold(iConfig.getParameter<double>(
"PromptPtThreshold")) {
52 std::vector<const GenParticle*> seeds;
56 std::vector<const GenParticle*>
egamma;
60 std::vector<const GenParticle*>
stable;
65 for(HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
p != myGenEvent->particles_end(); ++
p) {
68 ((*p)->status()==1&&(*p)->pdg_id() == 22) ||
69 ((*p)->status()==1&&
abs((*p)->pdg_id()) == 11))
73 if ((*p)->momentum().perp() >
ptSeedThr &&
81 if ((*p)->status() == 1) {
84 if (
abs((*p)->pdg_id()) == 211 ||
85 abs((*p)->pdg_id()) == 321 ||
86 abs((*p)->pdg_id()) == 11 ||
87 abs((*p)->pdg_id()) == 13 ||
88 abs((*p)->pdg_id()) == 15) {
90 if ((*p)->momentum().perp() >
ptTkThr &&
91 fabs((*p)->momentum().eta()) <
etaTkThr) {
97 if ((*p)->pdg_id() == 22 &&
100 egamma.push_back(*
p);
102 if (
abs((*p)->pdg_id()) == 11 &&
103 (*p)->momentum().perp() >
ptElThr &&
104 fabs((*p)->momentum().eta()) <
etaElThr) {
105 egamma.push_back(*
p);
110 if (seeds.size() < 2)
return accepted;
125 std::vector<TLorentzVector> candidate;
128 std::vector<TLorentzVector> candidateNarrow, candidateSeed;
130 std::vector<const GenParticle*>::iterator itSeed;
134 int first_different_id;
136 for(itSeed = seeds.begin(); itSeed != seeds.end(); ++itSeed) {
138 TLorentzVector energy, narrowCone, temp1, temp2, tempseed;
140 tempseed.SetXYZM((*itSeed)->momentum().px(), (*itSeed)->momentum().py(), (*itSeed)->momentum().pz(), 0);
141 for(
auto itEn = egamma.begin(); itEn != egamma.end(); ++itEn) {
142 temp1.SetXYZM((*itEn)->momentum().px(), (*itEn)->momentum().py(), (*itEn)->momentum().pz(), 0);
144 double DR = temp1.DeltaR(tempseed);
145 double DPhi = temp1.DeltaPhi(tempseed);
146 double DEta = (*itEn)->momentum().eta()-(*itSeed)->momentum().eta();
147 if(DPhi<0) DPhi=-
DPhi;
148 if(DEta<0) DEta=-DEta;
163 if ( energy.Et() != 0. ) {
166 temp2.SetXYZM(energy.Px(), energy.Py(), energy.Pz(), 0);
169 for(
auto itStable = stable.begin(); itStable != stable.end(); ++itStable) {
170 temp1.SetXYZM((*itStable)->momentum().px(), (*itStable)->momentum().py(), (*itStable)->momentum().pz(), 0);
171 double DR = temp1.DeltaR(temp2);
182 this_id = (*itSeed)->pdg_id();
184 while (mom->pdg_id() == this_id) {
186 const GenParticle* mother = mom->production_vertex() ?
187 *(mom->production_vertex()->particles_in_const_begin()) :
nullptr;
190 if (mom ==
nullptr) {
195 first_different_id = mom->pdg_id();
197 if (mom->status() == 2 &&
std::abs(first_different_id)>100) isPrompt=
false;
200 if(isPrompt) counter=0;
206 candidate.push_back(energy);
207 candidateSeed.push_back(tempseed);
208 candidateNarrow.push_back(narrowCone);
209 nTracks.push_back(counter);
212 if (candidate.size() <2)
return accepted;
214 TLorentzVector minvMin, minvMax;
223 for(
unsigned int i=0;
i<candidate.size()-1; ++
i) {
225 if (candidate[
i].Energy() <
energyCut)
continue;
229 for(
unsigned int j=
i+1; j<candidate.size(); ++j) {
232 if (candidate[j].Energy() <
energyCut)
continue;
237 if (nTracks[
i] + nTracks[j] >
nTkConeSum)
continue;
240 if (candidate[
i].
Pt() > candidate[j].Pt()) {
253 minvMin = candidate[
i] + candidate[j];
256 minvMax = candidate[
i] + candidate[j];
const unsigned int nTracks(const reco::Vertex &sv)
bool filter(const HepMC::GenEvent *myGenEvent) override
~PythiaHepMCFilterGammaGamma() override
const double etaMaxCandidate
const double pTMinCandidate2
Abs< T >::type abs(const T &t)
PythiaHepMCFilterGammaGamma(const edm::ParameterSet &)
const double promptPtThreshold
bool accepted(std::vector< std::string_view > const &, std::string_view)
const double pTMinCandidate1
const double dRNarrowCone