gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp
Go to the documentation of this file.
1 #include <vector>
2 #include <cmath>
3 #include <memory>
4 #include <iomanip>
5 
6 #include "gambit/ColliderBit/analyses/Analysis.hpp"
7 #include "gambit/ColliderBit/ATLASEfficiencies.hpp"
8 #include "gambit/ColliderBit/mt2_bisect.h"
9 
10 // #define CHECK_CUTFLOW
11 
12 
13 using namespace std;
14 
30 
31 namespace Gambit {
32  namespace ColliderBit {
33 
34  bool sortByPT_jet(const HEPUtils::Jet* jet1, const HEPUtils::Jet* jet2) { return (jet1->pT() > jet2->pT()); }
35  bool sortByPT_lep(const HEPUtils::Particle* lep1, const HEPUtils::Particle* lep2) { return (lep1->pT() > lep2->pT()); }
36 
37 
38  class Analysis_ATLAS_13TeV_PhotonGGM_36invfb : public Analysis {
39  private:
40 
41  // Numbers passing cuts
42  std::map<string, EventCounter> _counters = {
43  {"SRaa_SL", EventCounter("SRaa_SL")},
44  {"SRaa_SH", EventCounter("SRaa_SH")},
45  {"SRaa_WL", EventCounter("SRaa_WL")},
46  {"SRaa_WH", EventCounter("SRaa_WH")},
47  {"SRaj_L", EventCounter("SRaj_L")},
48  {"SRaj_L200", EventCounter("SRaj_L200")},
49  {"SRaj_H", EventCounter("SRaj_H")},
50  };
51 
52 
53  // Cut Flow
54  #ifdef CHECK_CUTFLOW
55  vector<int> cutFlowVector;
56  vector<double> cutFlowVector_ATLAS;
57  vector<string> cutFlowVector_str;
58  int NCUTS;
59  #endif
60 
61 
62  // Overlap removal -- discard from first list if deltaR < deltaR1 and
63  // discard from the second list deltaR1 < deltaR < deltaR2
64  /*
65  void JetParticleOverlapRemoval2(vector<const HEPUtils::Jet*> &jetvec, vector<const HEPUtils::Particle*> &particlevec, double deltaR1, double deltaR2, bool use_rapidity=false) {
66 
67  vector<const HEPUtils::Jet*> keep_jets;
68  vector<const HEPUtils::Jet*> kill_jets;
69 
70  vector<const HEPUtils::Particle*> keep_particles;
71  vector<const HEPUtils::Particle*> kill_particles;
72 
73  for(const HEPUtils::Jet* j : jetvec) {
74 
75  for(const HEPUtils::Particle* p : particlevec) {
76 
77  const double dR = (use_rapidity) ? j->mom().deltaR_rap(p->mom()) : j->mom().deltaR_eta(p->mom());
78 
79  if (fabs(dR) <= deltaR1) {
80  // If jet j is not in kill_jets, add it.
81  if ( find(kill_jets.begin(), kill_jets.end(), j) == kill_jets.end() ) kill_jets.push_back(j);
82  }
83  else if ((fabs(dR) > deltaR1) && (fabs(dR) <= deltaR2)) {
84  // If particle p is not in kill_particles, add it.
85  if ( find(kill_particles.begin(), kill_particles.end(), p) == kill_particles.end() ) kill_particles.push_back(p);
86  }
87  }
88  }
89 
90  // All jets that are not in kill_jets should go into keep_jets
91  for(const HEPUtils::Jet* j : jetvec) {
92  if ( find(kill_jets.begin(), kill_jets.end(), j) == kill_jets.end() ) keep_jets.push_back(j);
93  }
94 
95  // All particles that are not in kill_particles should go into keep_particles
96  for(const HEPUtils::Particle* p : particlevec) {
97  if ( find(kill_particles.begin(), kill_particles.end(), p) == kill_particles.end() ) keep_particles.push_back(p);
98  }
99 
100  jetvec = keep_jets;
101  particlevec = keep_particles;
102 
103  return;
104  }
105  */
106 
107  public:
108 
109  // Required detector sim
110  static constexpr const char* detector = "ATLAS";
111 
112  Analysis_ATLAS_13TeV_PhotonGGM_36invfb() {
113 
114  set_analysis_name("ATLAS_13TeV_PhotonGGM_36invfb");
115  set_luminosity(36.1);
116 
117  #ifdef CHECK_CUTFLOW
118  NCUTS= 51;
119  for(int i=0;i<NCUTS;i++){
120  cutFlowVector.push_back(0);
121  cutFlowVector_ATLAS.push_back(-1.0);
122  cutFlowVector_str.push_back("");
123  }
124  #endif
125 
126  }
127 
128  void run(const HEPUtils::Event* event) {
129 
130  // Missing energy w/ smearing
131  double ht = 0;
132  for (const HEPUtils::Particle* p : event->visible_particles()) ht += p->pT();
133  HEPUtils::P4 pmiss = event->missingmom();
134  ATLAS::smearMET(pmiss, ht);
135  const double met = pmiss.pT();
136 
137  // Baseline lepton objects
138  vector<const HEPUtils::Particle*> baselineElectrons, baselineMuons;
139 
140  // Baseline electrons
141  for (const HEPUtils::Particle* electron : event->electrons()) {
142  bool crack = (electron->abseta() > 1.37) && (electron->abseta() < 1.52);
143  if (electron->pT() > 25. && electron->abseta() < 2.47 && !crack) baselineElectrons.push_back(electron);
144  }
145 
146  // Apply electron efficiency
147  ATLAS::applyElectronEff(baselineElectrons);
148 
149  // Apply tight electron selection
150  ATLAS::applyTightIDElectronSelection(baselineElectrons);
151 
152  for (const HEPUtils::Particle* muon : event->muons()) {
153  if (muon->pT() > 25. && muon->abseta() < 2.7) baselineMuons.push_back(muon);
154  }
155 
156  // Apply muon efficiency
157  ATLAS::applyMuonEff(baselineMuons);
158 
159  // Photons
160  vector<const HEPUtils::Particle*> baselinePhotons;
161  for (const HEPUtils::Particle* photon : event->photons()) {
162  bool crack = (photon->abseta() > 1.37) && (photon->abseta() < 1.52);
163  if (photon->pT() > 25. && photon->abseta() < 2.37 && !crack) baselinePhotons.push_back(photon);
164  }
165  ATLAS::applyPhotonEfficiencyR2(baselinePhotons);
166 
167  // Jets
168  vector<const HEPUtils::Jet*> jets28;
169  vector<const HEPUtils::Jet*> jets28_nophooverlap;
170  for (const HEPUtils::Jet* jet : event->jets()) {
171  if (jet->pT() > 30. && fabs(jet->eta()) < 2.8) {
172  jets28.push_back(jet);
173  jets28_nophooverlap.push_back(jet);
174  }
175  }
176 
177 
178  // Overlap removal
179  const bool use_rapidity=true;
180  removeOverlap(baselinePhotons,baselineElectrons, 0.01, use_rapidity); // <-- taken from ATLAS code snippets on HEPData
181 
182  removeOverlap(jets28, baselineElectrons, 0.2, use_rapidity);
183  removeOverlap(jets28_nophooverlap, baselineElectrons, 0.2, use_rapidity);
184  removeOverlap(jets28, baselinePhotons, 0.2, use_rapidity);
185  removeOverlap(baselineElectrons, jets28_nophooverlap, 0.4, use_rapidity);
186  removeOverlap(baselineElectrons, jets28, 0.4, use_rapidity);
187  removeOverlap(baselinePhotons, jets28, 0.4, use_rapidity);
188  removeOverlap(baselineMuons, jets28, 0.4, use_rapidity);
189  removeOverlap(baselineMuons, jets28_nophooverlap, 0.4, use_rapidity);
190 
191 
192  // Make |eta| < 2.5 jets
193  vector<const HEPUtils::Jet*> jets25;
194  for (const HEPUtils::Jet* jet : jets28){
195  if (fabs(jet->eta()) < 2.5) jets25.push_back(jet);
196  }
197 
198  // Put objects in pT order
199  sortByPt(jets25);
200  sortByPt(jets28);
201  sortByPt(jets28_nophooverlap);
202  sortByPt(baselineElectrons);
203  sortByPt(baselineMuons);
204  sortByPt(baselinePhotons);
205 
206 
207  // Function used to get b jets
208  vector<const HEPUtils::Jet*> bJets25;
209  vector<const HEPUtils::Jet*> bJets28;
210 
211  const std::vector<double> a = {0,10.};
212  const std::vector<double> b = {0,10000.};
213  const std::vector<double> c = {0.77};
214  HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
215  for (const HEPUtils::Jet* jet : jets25) {
216  bool hasTag=has_tag(_eff2d, jet->abseta(), jet->pT());
217  if(jet->btag() && hasTag) bJets25.push_back(jet);
218  }
219 
220  for (const HEPUtils::Jet* jet : jets28) {
221  bool hasTag=has_tag(_eff2d, jet->abseta(), jet->pT());
222  if(jet->btag() && hasTag) bJets28.push_back(jet);
223  }
224 
225  // Multiplicities
226  int nLep = baselineElectrons.size() + baselineMuons.size();
227  int nJets25 = jets25.size();
228  int nPhotons = baselinePhotons.size();
229 
230 
231  // Pre-selection
232  bool preSelection2a=false;
233  if(nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75.) preSelection2a=true;
234 
235  bool preSelectionSRLaj = false;
236  if(nPhotons==1 && baselinePhotons[0]->pT() > 145.) preSelectionSRLaj=true;
237 
238  bool preSelectionSRHaj = false;
239  if(nPhotons==1 && baselinePhotons[0]->pT() > 400.) preSelectionSRHaj=true;
240 
241  // Useful variables
242  // "Photon-enhanced" HT, with no overlap removal of photons-jets
243  double HT = 0.;
244  for(const HEPUtils::Particle* photon : baselinePhotons) {
245  HT += photon->pT();
246  }
247  for(const HEPUtils::Particle* electron : baselineElectrons) {
248  HT += electron->pT();
249  }
250  for(const HEPUtils::Particle* muon : baselineMuons) {
251  HT += muon->pT();
252  }
253  for(const HEPUtils::Jet* jet : jets28_nophooverlap) {
254  HT += jet->pT();
255  }
256 
257  // meff
258  double meff = met;
259  for(const HEPUtils::Particle* photon : baselinePhotons) {
260  meff += photon->pT();
261  }
262  for(const HEPUtils::Particle* electron : baselineElectrons) {
263  meff += electron->pT();
264  }
265  for(const HEPUtils::Particle* muon : baselineMuons) {
266  meff += muon->pT();
267  }
268 
269  // Note that meff is only used for aj signal regions -> |jet eta| < 2.5
270  for(const HEPUtils::Jet* jet : jets25) {
271  meff += jet->pT();
272  }
273 
274  // dphimin(a,met)
275  double dphimin_amet = 999.;
276  if (nPhotons == 1) {
277  dphimin_amet = pmiss.deltaPhi(baselinePhotons.at(0)->mom());
278  }
279  else if (nPhotons >= 2) {
280  dphimin_amet = std::min( pmiss.deltaPhi(baselinePhotons.at(0)->mom()), pmiss.deltaPhi(baselinePhotons.at(1)->mom()) );
281  }
282 
283 
284  double dphimin_j25met = 999.;
285  if (jets25.size() == 1) {
286  dphimin_j25met = pmiss.deltaPhi(jets25.at(0)->mom());
287  }
288  else if (jets25.size() >= 2) {
289  dphimin_j25met = std::min( pmiss.deltaPhi(jets25.at(0)->mom()), pmiss.deltaPhi(jets25.at(1)->mom()) );
290  }
291 
292 
293  double dphimin_j28met = 999.;
294  if (jets28.size() == 1) {
295  dphimin_j28met = pmiss.deltaPhi(jets28.at(0)->mom());
296  }
297  else if (jets28.size() >= 2) {
298  dphimin_j28met = std::min( pmiss.deltaPhi(jets28.at(0)->mom()), pmiss.deltaPhi(jets28.at(1)->mom()) );
299  }
300 
301 
302  // RT4
303  // Only used in aj regions -> use |jet eta| < 2.5
304  double RT4 = 0.;
305  if(jets25.size() > 3){
306  RT4 = jets25[0]->pT() + jets25[1]->pT() + jets25[2]->pT() + jets25[3]->pT();
307  }
308  double denom=0.;
309  for(const HEPUtils::Jet* jet : jets25){
310  denom += jet->pT();
311  }
312  RT4=RT4/denom;
313 
314 
315 
316  // All variables are now done
317  // Increment signal region variables
318  // 2a regions
319  if(preSelection2a && met > 150. && HT > 2750 && dphimin_j28met > 0.5) _counters.at("SRaa_SL").add_event(event);
320  if(preSelection2a && met > 250. && HT > 2000 && dphimin_j28met > 0.5 && dphimin_amet > 0.5) _counters.at("SRaa_SH").add_event(event);
321  if(preSelection2a && met > 150. && HT > 1500 && dphimin_j28met > 0.5) _counters.at("SRaa_WL").add_event(event);
322  if(preSelection2a && met > 250. && HT > 1000 && dphimin_j28met > 0.5 && dphimin_amet > 0.5) _counters.at("SRaa_WH").add_event(event);
323 
324  // aj regions
325  if(preSelectionSRLaj && nJets25 >=5 && nLep == 0 && met > 300. && meff > 2000. && RT4 < 0.90 && dphimin_j25met > 0.5 && dphimin_amet > 0.5) _counters.at("SRaj_L").add_event(event);
326  if(preSelectionSRLaj && nJets25 >=5 && nLep == 0 && met > 200. && meff > 2000. && RT4 < 0.90 && dphimin_j25met > 0.5 && dphimin_amet > 0.5) _counters.at("SRaj_L200").add_event(event);
327  if(preSelectionSRHaj && nJets25 >=3 && nLep == 0 && met > 400. && meff > 2400. && dphimin_j25met > 0.5 && dphimin_amet > 0.5) _counters.at("SRaj_H").add_event(event);
328 
329 
330  #ifdef CHECK_CUTFLOW
331 
332  /* */
333  /*********************************************************/
334  cutFlowVector_str[0] = "Total ";
335  /*---------------------------------------*/
336  cutFlowVector_str[1] = "SBL: trigger && 2 photons";
337  cutFlowVector_str[2] = "SBL: PhotonsPt";
338  cutFlowVector_str[3] = "SBL: MET";
339  cutFlowVector_str[4] = "SBL: HT";
340  cutFlowVector_str[5] = "SBL: dPhiMin(jet,met)";
341  cutFlowVector_str[6] = "SBL: dPhiMin(gamma,met)";
342 
343  cutFlowVector_str[7] = "SBH: trigger && 2 photons";
344  cutFlowVector_str[8] = "SBH: PhotonsPt";
345  cutFlowVector_str[9] = "SBH: MET";
346  cutFlowVector_str[10] = "SBH: HT";
347  cutFlowVector_str[11] = "SBH: dPhiMin(jet,met)";
348  cutFlowVector_str[12] = "SBH: dPhiMin(gamma,met)";
349 
350  cutFlowVector_str[13] = "WBL: trigger && 2 photons";
351  cutFlowVector_ATLAS[13] = 26.6;
352  cutFlowVector_str[14] = "WBL: PhotonsPt";
353  cutFlowVector_ATLAS[14] = 21.3;
354  cutFlowVector_str[15] = "WBL: MET";
355  cutFlowVector_ATLAS[15] = 16.9;
356  cutFlowVector_str[16] = "WBL: HT";
357  cutFlowVector_ATLAS[16] = 14.7;
358  cutFlowVector_str[17] = "WBL: dPhiMin(jet,met)";
359  cutFlowVector_ATLAS[17] = 11.0;
360  cutFlowVector_str[18] = "WBL: dPhiMin(gamma,met)";
361  cutFlowVector_ATLAS[18] = 11.0;
362 
363  cutFlowVector_str[19] = "WBH: trigger && 2 photons";
364  cutFlowVector_ATLAS[19] = 19.6;
365  cutFlowVector_str[20] = "WBH: PhotonsPt";
366  cutFlowVector_ATLAS[20] = 19.2;
367  cutFlowVector_str[21] = "WBH: MET";
368  cutFlowVector_ATLAS[21] = 15.6;
369  cutFlowVector_str[22] = "WBH: HT";
370  cutFlowVector_ATLAS[22] = 15.6;
371  cutFlowVector_str[23] = "WBH: dPhiMin(jet,met)";
372  cutFlowVector_ATLAS[23] = 14.8;
373  cutFlowVector_str[24] = "WBH: dPhiMin(gamma,met)";
374  cutFlowVector_ATLAS[24] = 14.6;
375 
376  cutFlowVector_str[25] = "SRL: trigger && 1 photon";
377  cutFlowVector_str[26] = "SRL: lepton veto";
378  cutFlowVector_str[27] = "SRL: pT_gamma";
379  cutFlowVector_str[28] = "SRL: met";
380  cutFlowVector_str[29] = "SRL: Njets";
381  cutFlowVector_str[30] = "SRL: dphimin(jet,met)";
382  cutFlowVector_str[31] = "SRL: dphimin(gamma,met)";
383  cutFlowVector_str[32] = "SRL: meff";
384  cutFlowVector_str[33] = "SRL: RT4";
385 
386  cutFlowVector_str[34] = "SRL200: trigger && 1 photon";
387  cutFlowVector_str[35] = "SRL200: lepton veto";
388  cutFlowVector_str[36] = "SRL200: pT_gamma";
389  cutFlowVector_str[37] = "SRL200: met";
390  cutFlowVector_str[38] = "SRL200: Njets";
391  cutFlowVector_str[39] = "SRL200: dphimin(jet,met)";
392  cutFlowVector_str[40] = "SRL200: dphimin(gamma,met)";
393  cutFlowVector_str[41] = "SRL200: meff";
394  cutFlowVector_str[42] = "SRL200: RT4";
395 
396  cutFlowVector_str[43] = "SRH: trigger && 1 photon";
397  cutFlowVector_str[44] = "SRH: lepton veto";
398  cutFlowVector_str[45] = "SRH: pT_gamma";
399  cutFlowVector_str[46] = "SRH: met";
400  cutFlowVector_str[47] = "SRH: Njets";
401  cutFlowVector_str[48] = "SRH: dphimin(jet,met)";
402  cutFlowVector_str[49] = "SRH: dphimin(gamma,met)";
403  cutFlowVector_str[50] = "SRH: meff";
404 
405  for(int j=0;j<NCUTS;j++){
406  if(
407  (j==0) ||
408 
409  /*
410  cutFlowVector_str[1] = "SBL: trigger && 2 photons";
411  cutFlowVector_str[2] = "SBL: PhotonsPt";
412  cutFlowVector_str[3] = "SBL: MET";
413  cutFlowVector_str[4] = "SBL: HT";
414  cutFlowVector_str[5] = "SBL: dPhiMin(jet,met)";
415  cutFlowVector_str[6] = "SBL: dPhiMin(gamma,met)";
416  */
417 
418  (j==1 && nPhotons==2 && baselinePhotons[0]->pT() > 35. && baselinePhotons[1]->pT() > 25.) ||
419  (j==2 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75.) ||
420  (j==3 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 150.) ||
421  (j==4 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 150. && HT > 2750.) ||
422  (j==5 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 150. && HT > 2750. && dphimin_j28met > 0.5) ||
423  (j==6 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 150. && HT > 2750. && dphimin_j28met > 0.5) || // No extra cut in this case
424 
425 
426  /*
427  cutFlowVector_str[7] = "SBH: trigger && 2 photons";
428  cutFlowVector_str[8] = "SBH: PhotonsPt";
429  cutFlowVector_str[9] = "SBH: MET";
430  cutFlowVector_str[10] = "SBH: HT";
431  cutFlowVector_str[11] = "SBH: dPhiMin(jet,met)";
432  cutFlowVector_str[12] = "SBH: dPhiMin(gamma,met)";
433  */
434 
435  (j==7 && nPhotons==2 && baselinePhotons[0]->pT() > 35. && baselinePhotons[1]->pT() > 25.) ||
436  (j==8 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75.) ||
437  (j==9 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 250.) ||
438  (j==10 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 250. && HT > 2000.) ||
439  (j==11 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 250. && HT > 2000. && dphimin_j28met > 0.5) ||
440  (j==12 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 250. && HT > 2000. && dphimin_j28met > 0.5 && dphimin_amet > 0.5) ||
441 
442 
443  /*
444  cutFlowVector_str[13] = "WBL: trigger && 2 photons";
445  cutFlowVector_str[14] = "WBL: PhotonsPt";
446  cutFlowVector_str[15] = "WBL: MET";
447  cutFlowVector_str[16] = "WBL: HT";
448  cutFlowVector_str[17] = "WBL: dPhiMin(jet,met)";
449  cutFlowVector_str[18] = "WBL: dPhiMin(gamma,met)";
450  */
451 
452  (j==13 && nPhotons==2 && baselinePhotons[0]->pT() > 35. && baselinePhotons[1]->pT() > 25.) ||
453  (j==14 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75.) ||
454  (j==15 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 150.) ||
455  (j==16 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 150. && HT > 1500.) ||
456  (j==17 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 150. && HT > 1500. && dphimin_j28met > 0.5) ||
457  (j==18 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 150. && HT > 1500. && dphimin_j28met > 0.5) || // no additional cut in this case
458 
459 
460  /*
461  cutFlowVector_str[19] = "WBH: trigger && 2 photons";
462  cutFlowVector_str[20] = "WBH: PhotonsPt";
463  cutFlowVector_str[21] = "WBH: MET";
464  cutFlowVector_str[22] = "WBH: HT";
465  cutFlowVector_str[23] = "WBH: dPhiMin(jet,met)";
466  cutFlowVector_str[24] = "WBH: dPhiMin(gamma,met)";
467  */
468 
469  (j==19 && nPhotons==2 && baselinePhotons[0]->pT() > 35. && baselinePhotons[1]->pT() > 25.) ||
470  (j==20 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75.) ||
471  (j==21 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 250.) ||
472  (j==22 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 250. && HT > 1000.) ||
473  (j==23 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 250. && HT > 1000. && dphimin_j28met > 0.5) ||
474  (j==24 && nPhotons==2 && baselinePhotons[0]->pT() > 75. && baselinePhotons[1]->pT() > 75. && met > 250. && HT > 1000. && dphimin_j28met > 0.5 && dphimin_amet > 0.5) || // no additional cut in this case
475 
476 
477  // --------
478 
479 
480  /*
481  cutFlowVector_str[25] = "SRL: trigger && 1 photon";
482  cutFlowVector_str[26] = "SRL: lepton veto";
483  cutFlowVector_str[27] = "SRL: pT_gamma";
484  cutFlowVector_str[28] = "SRL: met";
485  cutFlowVector_str[29] = "SRL: Njets";
486  cutFlowVector_str[30] = "SRL: dphimin(jet,met)";
487  cutFlowVector_str[31] = "SRL: dphimin(gamma,met)";
488  cutFlowVector_str[32] = "SRL: meff";
489  cutFlowVector_str[33] = "SRL: RT4";
490  */
491 
492  (j==25 && nPhotons==1 && baselinePhotons[0]->pT() > 140.) ||
493  (j==26 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 140.) ||
494  (j==27 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145.) ||
495  (j==28 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 300.) ||
496  (j==29 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 300. && nJets25 >= 5) ||
497  (j==30 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 300. && nJets25 >= 5 && dphimin_j25met > 0.4) ||
498  (j==31 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 300. && nJets25 >= 5 && dphimin_j25met > 0.4 && dphimin_amet > 0.4) ||
499  (j==32 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 300. && nJets25 >= 5 && dphimin_j25met > 0.4 && dphimin_amet > 0.4 && meff > 2000.) ||
500  (j==33 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 300. && nJets25 >= 5 && dphimin_j25met > 0.4 && dphimin_amet > 0.4 && meff > 2000. && RT4 < 0.90) ||
501 
502 
503  /*
504  cutFlowVector_str[25] = "SRL200: trigger && 1 photon";
505  cutFlowVector_str[26] = "SRL200: lepton veto";
506  cutFlowVector_str[27] = "SRL200: pT_gamma";
507  cutFlowVector_str[28] = "SRL200: met";
508  cutFlowVector_str[29] = "SRL200: Njets";
509  cutFlowVector_str[30] = "SRL200: dphimin(jet,met)";
510  cutFlowVector_str[31] = "SRL200: dphimin(gamma,met)";
511  cutFlowVector_str[32] = "SRL200: meff";
512  cutFlowVector_str[33] = "SRL200: RT4";
513  */
514 
515  (j==34 && nPhotons==1 && baselinePhotons[0]->pT() > 140.) ||
516  (j==35 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 140.) ||
517  (j==36 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145.) ||
518  (j==37 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 200.) ||
519  (j==38 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 200. && nJets25 >= 5) ||
520  (j==39 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 200. && nJets25 >= 5 && dphimin_j25met > 0.4) ||
521  (j==40 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 200. && nJets25 >= 5 && dphimin_j25met > 0.4 && dphimin_amet > 0.4) ||
522  (j==41 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 200. && nJets25 >= 5 && dphimin_j25met > 0.4 && dphimin_amet > 0.4 && meff > 2000.) ||
523  (j==42 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 145. && met > 200. && nJets25 >= 5 && dphimin_j25met > 0.4 && dphimin_amet > 0.4 && meff > 2000. && RT4 < 0.90) ||
524 
525 
526  /*
527  cutFlowVector_str[34] = "SRH: trigger && 1 photon";
528  cutFlowVector_str[35] = "SRH: lepton veto";
529  cutFlowVector_str[36] = "SRH: pT_gamma";
530  cutFlowVector_str[37] = "SRH: met";
531  cutFlowVector_str[38] = "SRH: Njets";
532  cutFlowVector_str[39] = "SRH: dphimin(jet,met)";
533  cutFlowVector_str[40] = "SRH: dphimin(gamma,met)";
534  cutFlowVector_str[41] = "SRH: meff";
535  */
536 
537  (j==43 && nPhotons==1 && baselinePhotons[0]->pT() > 140.) ||
538  (j==44 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 140.) ||
539  (j==45 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 400.) ||
540  (j==46 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 400. && met > 400.) ||
541  (j==47 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 400. && met > 400. && nJets25 >= 3) ||
542  (j==48 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 400. && met > 400. && nJets25 >= 3 && dphimin_j25met > 0.4) ||
543  (j==49 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 400. && met > 400. && nJets25 >= 3 && dphimin_j25met > 0.4 && dphimin_amet > 0.4) ||
544  (j==50 && nPhotons==1 && nLep==0 && baselinePhotons[0]->pT() > 400. && met > 400. && nJets25 >= 3 && dphimin_j25met > 0.4 && dphimin_amet > 0.4 && meff > 2400.)
545 
546  )cutFlowVector[j]++;
547 
548  }
549 
550  #endif // end #ifdef CHECK_CUTFLOW
551 
552  return;
553 
554  }
555 
557  void combine(const Analysis* other)
558  {
559  const Analysis_ATLAS_13TeV_PhotonGGM_36invfb* specificOther
560  = dynamic_cast<const Analysis_ATLAS_13TeV_PhotonGGM_36invfb*>(other);
561 
562  #ifdef CHECK_CUTFLOW
563  if (NCUTS != specificOther->NCUTS) NCUTS = specificOther->NCUTS;
564  for (int j=0; j<NCUTS; j++) {
565  cutFlowVector[j] += specificOther->cutFlowVector[j];
566  cutFlowVector_str[j] = specificOther->cutFlowVector_str[j];
567  }
568  #endif
569 
570  for (auto& pair : _counters) { pair.second += specificOther->_counters.at(pair.first); }
571 
572  }
573 
574 
575  void collect_results() {
576 
577  #ifdef CHECK_CUTFLOW
578  double scale_by= 70.8 / 1.0e4;
579  cout << "------------------------------------------------------------------------------------------------------------------------------ "<<endl;
580  cout << "CUT FLOW: ATLAS_13TeV_PhotonGGM_36invfb "<<endl;
581  cout << "------------------------------------------------------------------------------------------------------------------------------"<<endl;
582  cout << right << setw(40) << "CUT," << setw(20) << "RAW," << setw(20) << "SCALED,"
583  << setw(20) << "%," << setw(20) << "ATLAS," << setw(20) << "GAMBIT (scaled) /ATLAS" << endl;
584  for (int j=0; j<NCUTS; j++) {
585  cout << right << setw(40) << cutFlowVector_str[j].c_str() << "," << setw(20)
586  << cutFlowVector[j] << "," << setw(20) << cutFlowVector[j]*scale_by << "," << setw(20)
587  << 100.*cutFlowVector[j]/cutFlowVector[0] << "%," << setw(20) << cutFlowVector_ATLAS[j] << "," << setw(20) << (cutFlowVector[j]*scale_by / cutFlowVector_ATLAS[j]) << endl;
588  }
589  cout << "------------------------------------------------------------------------------------------------------------------------------ "<<endl;
590  #endif
591 
592  // add_result(SignalRegionData("SR label", n_obs, {n_sig_MC, n_sig_MC_sys}, {n_bkg, n_bkg_err}));
593 
594  add_result(SignalRegionData(_counters.at("SRaa_SL"), 0., { 0.50, 0.30}));
595  add_result(SignalRegionData(_counters.at("SRaa_SH"), 0., { 0.48, 0.30}));
596  add_result(SignalRegionData(_counters.at("SRaa_WL"), 6., { 3.7, 1.1}));
597  add_result(SignalRegionData(_counters.at("SRaa_WH"), 1., { 2.05, 0.65}));
598  add_result(SignalRegionData(_counters.at("SRaj_L"), 4., { 1.33, 0.54}));
599  add_result(SignalRegionData(_counters.at("SRaj_L200"), 8., { 2.68, 0.64}));
600  add_result(SignalRegionData(_counters.at("SRaj_H"), 3., { 1.14, 0.61}));
601 
602  return;
603  }
604 
605 
606  protected:
607  void analysis_specific_reset() {
608 
609  for (auto& pair : _counters) { pair.second.reset(); }
610 
611  #ifdef CHECK_CUTFLOW
612  std::fill(cutFlowVector.begin(), cutFlowVector.end(), 0);
613  #endif
614  }
615 
616  };
617 
618  // Factory function
619  DEFINE_ANALYSIS_FACTORY(ATLAS_13TeV_PhotonGGM_36invfb)
620 
621  }
622 }
Gambit::ColliderBit::has_tag
bool has_tag(const HEPUtils::BinnedFn2D< double > &effmap, double eta, double pt)
Randomly get a tag result (can be anything) from a 2D |eta|-pT efficiency map.
Definition: Utils.hpp:207
DEFINE_ANALYSIS_FACTORY
#define DEFINE_ANALYSIS_FACTORY(ANAME)
For analysis factory function definition.
Definition: Analysis.hpp:124
Gambit::ColliderBit::ATLAS::smearMET
void smearMET(HEPUtils::P4 &pmiss, double set)
Randomly smear the MET vector by parameterised resolutions.
Definition: ATLASEfficiencies.hpp:266
Gambit::ColliderBit::sortByPT_lep
bool sortByPT_lep(const HEPUtils::Particle *lep1, const HEPUtils::Particle *lep2)
Definition: Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp:35
Gambit::ColliderBit::ATLAS::applyPhotonEfficiencyR2
void applyPhotonEfficiencyR2(std::vector< const HEPUtils::Particle *> &photons)
Definition: ATLASEfficiencies.hpp:152
Gambit::ColliderBit::removeOverlap
void removeOverlap(MOMPTRS1 &momstofilter, const MOMPTRS2 &momstocmp, double deltaRMax, bool use_rapidity=false, double pTmax=DBL_MAX, double btageff=0)
Overlap removal – discard from first list if within deltaRMax of any from the second list Optional a...
Definition: Utils.hpp:299
run
void run(bool, bool, bool, int, double, double, int, int, int, int, int, double, char[], int, int[], bool, bool, bool, bool, double, int, double(*)(double *, int, int, void *), void(*)(int, int, int, double *, double *, double *, double, double, double, void *), void *)
std
STL namespace.
b
START_MODEL b
Definition: demo.hpp:270
Gambit::ColliderBit::sortByPt
void sortByPt(ParticlePtrs &particles)
Definition: Utils.hpp:389
Gambit::ColliderBit::EventCounter
A simple class for counting events of type HEPUtils::Event.
Definition: EventCounter.hpp:26
mt2_bisect.h
Gambit::ColliderBit::Analysis
A class for collider analyses within ColliderBit.
Definition: Analysis.hpp:41
Gambit::ColliderBit::sortByPT_jet
bool sortByPT_jet(const HEPUtils::Jet *jet1, const HEPUtils::Jet *jet2)
Definition: Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp:34
Gambit::ColliderBit::SignalRegionData
A simple container for the result of one signal region from one analysis.
Definition: AnalysisData.hpp:57
Gambit::ColliderBit::ATLAS::applyTightIDElectronSelection
void applyTightIDElectronSelection(std::vector< const HEPUtils::Particle *> &electrons)
Efficiency function for Tight ID electrons.
Definition: ATLASEfficiencies.hpp:446
Gambit::ColliderBit::ATLAS::applyMuonEff
void applyMuonEff(std::vector< const HEPUtils::Particle *> &muons)
Randomly filter the supplied particle list by parameterised muon efficiency.
Definition: ATLASEfficiencies.hpp:85
colouring.combine
def combine(region_file, pip_file)
Definition: colouring.py:169
Gambit::ColliderBit::ATLAS::applyElectronEff
void applyElectronEff(std::vector< const HEPUtils::Particle *> &electrons)
Randomly filter the supplied particle list by parameterised electron efficiency.
Definition: ATLASEfficiencies.hpp:62
Gambit
TODO: see if we can use this one:
Definition: Analysis.hpp:33
Analysis.hpp
Class for ColliderBit analyses.
ATLASEfficiencies.hpp
Functions that do super fast ATLAS detector simulation based on four-vector smearing.