CMS 3D CMS Logo

PythiaHepMCFilterGammaGamma.cc
Go to the documentation of this file.
3 #include "Math/GenVector/VectorUtil.h"
4 //#include "CLHEP/HepMC/GenParticle.h"
6 
7 #include "TLorentzVector.h"
8 
9 #include <iostream>
10 
11 using namespace edm;
12 using namespace std;
13 using namespace HepMC;
14 
15 
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")) {
40 
41 }
42 
44 {
45 }
46 
48 
49  bool accepted = false;
50 
51  // electron/photon seeds
52  std::vector<const GenParticle*> seeds;
53 
54  // other electrons/photons to be added to seeds
55  // to form candidates
56  std::vector<const GenParticle*> egamma;
57 
58  // charged tracks to be taken into account in the isolation cones
59  // around candidates
60  std::vector<const GenParticle*> stable;
61 
62  //----------
63  // 1. find electron/photon seeds
64  //----------
65  for(HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
66 
67  if (
68  ((*p)->status()==1&&(*p)->pdg_id() == 22) || // gamma
69  ((*p)->status()==1&&abs((*p)->pdg_id()) == 11)) // electron
70 
71  {
72  // check for eta and pT threshold for seed in gamma, el
73  if ((*p)->momentum().perp() > ptSeedThr &&
74  fabs((*p)->momentum().eta()) < etaSeedThr) {
75 
76  seeds.push_back(*p);
77  }
78  }
79 
80 
81  if ((*p)->status() == 1) {
82 
83  // save charged stable tracks
84  if (abs((*p)->pdg_id()) == 211 || // charged pion
85  abs((*p)->pdg_id()) == 321 || // charged kaon
86  abs((*p)->pdg_id()) == 11 || // electron
87  abs((*p)->pdg_id()) == 13 || // muon
88  abs((*p)->pdg_id()) == 15) { // tau
89  // check if it passes the cut
90  if ((*p)->momentum().perp() > ptTkThr &&
91  fabs((*p)->momentum().eta()) < etaTkThr) {
92  stable.push_back(*p);
93  }
94  }
95 
96  // save egamma for candidate calculation
97  if ((*p)->pdg_id() == 22 &&
98  (*p)->momentum().perp() > ptGammaThr &&
99  fabs((*p)->momentum().eta()) < etaGammaThr) {
100  egamma.push_back(*p);
101  }
102  if (abs((*p)->pdg_id()) == 11 &&
103  (*p)->momentum().perp() > ptElThr &&
104  fabs((*p)->momentum().eta()) < etaElThr) {
105  egamma.push_back(*p);
106  }
107  }
108  }
109 
110  if (seeds.size() < 2) return accepted;
111 
112  //----------
113  // 2. loop over seeds to build candidates
114  //
115  // (adding nearby electrons/photons
116  // to the seed electrons/photons to obtain the total
117  // electromagnetic energy)
118  //----------
119 
120  // number of tracks around each of the candidates
121  std::vector<int> nTracks;
122 
123  // the candidates (four momenta) formed from the
124  // seed electrons/photons and nearby electrons/photons
125  std::vector<TLorentzVector> candidate;
126 
127  // these are filled but then not used afterwards (could be removed)
128  std::vector<TLorentzVector> candidateNarrow, candidateSeed;
129 
130  std::vector<const GenParticle*>::iterator itSeed;
131 
132  const GenParticle* mom;
133  int this_id;
134  int first_different_id;
135 
136  for(itSeed = seeds.begin(); itSeed != seeds.end(); ++itSeed) {
137 
138  TLorentzVector energy, narrowCone, temp1, temp2, tempseed;
139 
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);
143 
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;
149 
150  // accept if within cone or within rectangular region around seed
151  if (DR < dRSeedMax || (DPhi<dPhiSeedMax&&DEta<dEtaSeedMax)) {
152  energy += temp1;
153  }
154  if (DR < dRNarrowCone) {
155  narrowCone += temp1;
156  }
157  }
158 
159  // number of stable charged particles found within dRTkMax
160  // around candidate
161  int counter = 0;
162 
163  if ( energy.Et() != 0. ) {
164  if (fabs(energy.Eta()) < etaMaxCandidate) {
165 
166  temp2.SetXYZM(energy.Px(), energy.Py(), energy.Pz(), 0);
167 
168  // count number of stable particles within cone around candidate
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);
172  if (DR < dRTkMax) counter++;
173  }
174 
175  if(acceptPrompts) {
176 
177  if ((*itSeed)->momentum().perp()>promptPtThreshold)
178  {
179  // check if *itSeed is a prompt particle
180 
181  bool isPrompt=true;
182  this_id = (*itSeed)->pdg_id();
183  mom = (*itSeed);
184  while (mom->pdg_id() == this_id) {
185 
186  const GenParticle* mother = mom->production_vertex() ?
187  *(mom->production_vertex()->particles_in_const_begin()) : nullptr;
188 
189  mom = mother;
190  if (mom == nullptr) {
191  break;
192  }
193  }
194 
195  first_different_id = mom->pdg_id();
196 
197  if (mom->status() == 2 && std::abs(first_different_id)>100) isPrompt=false;
198 
199  // ignore charged particles around prompt particles
200  if(isPrompt) counter=0;
201  }
202  }
203  }
204  }
205 
206  candidate.push_back(energy);
207  candidateSeed.push_back(tempseed);
208  candidateNarrow.push_back(narrowCone);
209  nTracks.push_back(counter);
210  }
211 
212  if (candidate.size() <2) return accepted;
213 
214  TLorentzVector minvMin, minvMax;
215 
216  //----------
217  // 3. perform further checks on candidates
218  //
219  // (energy, charged isolation requirements etc.)
220  //----------
221 
222  int i1, i2;
223  for(unsigned int i=0; i<candidate.size()-1; ++i) {
224 
225  if (candidate[i].Energy() < energyCut) continue;
226  if(nTracks[i]>nTkConeMax) continue;
227  if (fabs(candidate[i].Eta()) > etaMaxCandidate) continue;
228 
229  for(unsigned int j=i+1; j<candidate.size(); ++j) {
230 
231  // check features of second candidate alone
232  if (candidate[j].Energy() < energyCut) continue;
233  if(nTracks[j]>nTkConeMax) continue;
234  if (fabs(candidate[j].Eta()) > etaMaxCandidate) continue;
235 
236  // check requirement on sum of tracks in both isolation cones
237  if (nTracks[i] + nTracks[j] > nTkConeSum) continue;
238 
239  // swap candidates to have pt[i1] >= pt[i2]
240  if (candidate[i].Pt() > candidate[j].Pt()) {
241  i1 = i;
242  i2 = j;
243  }
244  else {
245  i1 = j;
246  i2 = i;
247  }
248 
249  // require minimum pt on leading and subleading candidate
250  if (candidate[i1].Pt() < pTMinCandidate1 || candidate[i2].Pt() < pTMinCandidate2) continue;
251 
252  // apply requirements on candidate pair mass
253  minvMin = candidate[i] + candidate[j];
254  if (minvMin.M() < invMassMin) continue;
255 
256  minvMax = candidate[i] + candidate[j];
257  if (minvMax.M() > invMassMax) continue;
258 
259  accepted = true;
260 
261  }
262  }
263 
264  return accepted;
265 }
266 
const unsigned int nTracks(const reco::Vertex &sv)
bool filter(const HepMC::GenEvent *myGenEvent) override
static const int stable
Definition: TopGenEvent.h:11
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
PythiaHepMCFilterGammaGamma(const edm::ParameterSet &)
bool accepted(std::vector< std::string_view > const &, std::string_view)
HLT enums.