gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb.cpp
Go to the documentation of this file.
1 //
2 // Created by dsteiner on 31/07/18.
3 // Amended by Martin White on 08/03/19.
4 //
5 // Based on https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/SUSY-2012-10/ (arXiv:1209.2102)
6 
7 //#include <gambit/ColliderBit/colliders/SpecializablePythia.hpp>
12 
13 namespace Gambit {
14  namespace ColliderBit {
15 
16  using namespace std;
17  using namespace HEPUtils;
18 
19  class Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb : public Analysis {
20 public:
21 
22 // Required detector sim
23 static constexpr const char* detector = "ATLAS";
24 
25 
26 #define CUTFLOWMAP(X) \
27  X(Total_events) \
28  X(electron_eq_4_jets) \
29  X(electron_met_gt_40) \
30  X(electron_eq_2_bjets) \
31  X(electron_sr) \
32  X(muon_eq_4_jets) \
33  X(muon_met_gt_40) \
34  X(muon_eq_2_bjets) \
35  X(muon_sr) \
36  X(twoLep_met_gt_40) \
37  X(twoLep_gt_1_bjet) \
38  X(mll_lt_81) \
39  X(mll_gt_30_lt_81) \
40  X(num_2lsr1) \
41  X(num_2lsr2)
42 #define f(x) x,
43 #define g(x) #x,
44  const std::vector<std::string> cutflowNames = {CUTFLOWMAP(g)};
45  enum cutflowEnum {CUTFLOWMAP(f)};
46 #undef g
47 #undef f
48 #undef CUTFLOWMAP
49 
50 #define VARMAP(X) \
51  X(mTopHad) \
52  X(oneLepSqrtS) \
53  X(mllKey) \
54  X(twoLepSqrtS)
55 #define f(x) x,
56 #define g(x) #x,
57  const std::vector<std::string> varNames = {VARMAP(g)};
58  enum varEnum {VARMAP(f)};
59 #undef g
60 #undef f
61 #undef VARMAP
62 
63  std::map<std::string, std::vector<double>> varResults;
64 
65  std::map<std::string, int> cutflows;
66 
67  double num1LSR=0;
68  double num2LSR1=0;
69  double num2LSR2=0;
70 
71  std::vector<double> calcNuPz(double Mw, P4 metMom, P4 lepMom)
72  {
73  double mu = sqr(Mw) / 2.0 + metMom.px() * lepMom.px() + metMom.py() * lepMom.py();
74  double denom = lepMom.E2() - lepMom.pz2();
75  double a = mu * lepMom.pz() / denom;
76  double a2 = sqr(a);
77  double b = (lepMom.E2() * metMom.E2() - sqr(mu)) / denom;
78  double delta = a2 - b;
79  if (delta < 0)
80  {
81  return {a};
82  }
83  else
84  {
85  return {a + std::sqrt(delta), a - std::sqrt(delta)};
86  }
87  }
88 
89  P4 getBestHadronicTop(
90  std::vector<const Jet *> bJets,
91  std::vector<const Jet *> lightJets,
92  const P4& leptonMom,
93  const P4& metMom,
94  double width,
95  double mean
96  )
97  {
98  // gaussian probability density function
99  auto prob = [&width, &mean](P4 particle) {
100  return 1 - std::erf(1.0 * std::abs(particle.m() - mean) / (std::sqrt(2.0) * width));
101  };
102 
103  double pTotal = 0.0;
104  P4 bestHadronicTop;
105  std::vector<double> nuPzChoices = calcNuPz(80.0, metMom, leptonMom);
106  P4 nu;
107  for (double nuPz : nuPzChoices)
108  {
109  double nuE = std::sqrt(sqr(metMom.px()) + sqr(metMom.py()) + sqr(nuPz));
110  nu.setPE(metMom.px(), metMom.py(), nuPz, nuE);
111  P4 WLep = leptonMom + nu;
112  // go through every bJet
113  for (const Jet* firstBJet : bJets)
114  {
115  for (const Jet* secondBJet : bJets)
116  {
117  if (firstBJet == secondBJet)
118  {
119  continue;
120  }
121  P4 topLep = *firstBJet + WLep;
122  // go through every combination of two light jets
123  for (const Jet* firstLightJet : lightJets)
124  {
125  for (const Jet* secondLightJet : lightJets)
126  {
127  // don't want to use a light jet with itself
128  if (firstLightJet == secondLightJet)
129  {
130  continue;
131  }
132  P4 WHad = *firstLightJet + *secondLightJet;
133  P4 topHad = *secondBJet + WHad;
134  // calculate a new probability
135  double newPTotal = prob(topHad) * prob(WHad) * prob(topLep) * prob(WLep);
136  if (newPTotal > pTotal)
137  {
138  // update the best values
139  pTotal = newPTotal;
140  bestHadronicTop = topHad;
141  }
142  }
143  }
144  }
145  }
146  }
147  return bestHadronicTop;
148  }
149 
150  double calcMt(P4 metVec, P4 lepVec)
151  {
152  double Met = metVec.pT();
153  double pTLep = lepVec.pT();
154  return std::sqrt(2 * pTLep * Met - 2 * AnalysisUtil::dot2D(lepVec, metVec));
155  }
156 
157 
158  double calcSqrtSSubMin(P4 visibleSubsystem, P4 invisbleSubsystem)
159  {
160  double visiblePart = std::sqrt(sqr(visibleSubsystem.m()) + sqr(visibleSubsystem.pT()));
161  double invisiblePart = invisbleSubsystem.pT();
162  double twoDimensionalVecSum =
163  sqr(visibleSubsystem.px() + invisbleSubsystem.px())
164  + sqr(visibleSubsystem.py() + invisbleSubsystem.py());
165  return std::sqrt(sqr(visiblePart + invisiblePart) - twoDimensionalVecSum);
166  }
167 
168  void getBJets(
169  std::vector<const Jet*>& jets,
170  std::vector<const Jet*>* bJets,
171  std::vector<const Jet*>* lightJets)
172  {
174  const std::vector<double> a = {0, 10.};
175  const std::vector<double> b = {0, 10000.};
176  const std::vector<double> c = {0.60};
177  BinnedFn2D<double> _eff2d(a,b,c);
178  for (const Jet* jet : jets)
179  {
180  bool hasTag = has_tag(_eff2d, jet->eta(), jet->pT());
181  if(jet->btag() && hasTag && jet->abseta() < 2.5)
182  {
183  bJets->push_back(jet);
184  }
185  else
186  {
187  lightJets->push_back(jet);
188  }
189  }
190  }
191 
192 
196 Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb()
197  {
198  set_analysis_name("ATLAS_7TeV_1OR2LEPStop_4_7invfb");
199  set_luminosity(4.7);
200  }
201 
206  void run(const HEPUtils::Event* event)
207  {
208  // TODO: take log of plots and constrain the plot range
209  //HEPUtilsAnalysis::analyze(event);
210  //cout << "Event number: " << num_events() << endl;
211  incrementCut(Total_events);
212  std::vector<const Particle*> electrons = event->electrons();
213  std::vector<const Particle*> muons = event->muons();
214  std::vector<const Jet*> jets = event->jets();
215 
216  electrons = AnalysisUtil::filterPtEta(electrons, 20, 2.47);
217  muons = AnalysisUtil::filterPtEta(muons, 10, 2.4);
218  jets = AnalysisUtil::filterPtEta(jets, 20, 4.5);
219 
220  std::vector<const Jet*> bJets, lightJets;
221  getBJets(jets, &bJets, &lightJets);
222 
223  jets = AnalysisUtil::jetLeptonOverlapRemoval(jets, electrons, 0.2);
224  electrons = AnalysisUtil::leptonJetOverlapRemoval(electrons, jets, 0.4);
225  muons = AnalysisUtil::leptonJetOverlapRemoval(muons, jets, 0.4);
226 
227  jets = AnalysisUtil::filterMaxEta(jets, 2.5);
228 
230 
231  std::vector<const Particle*> leptons = AnalysisUtil::getSortedLeptons({electrons, muons});
232  std::sort(electrons.begin(), electrons.end(), AnalysisUtil::sortParticlesByPt);
233  std::sort(muons.begin(), muons.end(), AnalysisUtil::sortParticlesByPt);
234  std::sort(jets.begin(), jets.end(), AnalysisUtil::sortJetsByPt);
235  std::sort(bJets.begin(), bJets.end(), AnalysisUtil::sortJetsByPt);
236  std::sort(lightJets.begin(), lightJets.end(), AnalysisUtil::sortJetsByPt);
237 
238  size_t
239  nLeptons = leptons.size(),
240  nJets = jets.size(),
241  nBJets = bJets.size(),
242  nLightJets = lightJets.size();
243 
244  double Met = event->met();
245  const P4& metVec = event->missingmom();
246 
247  if (nLeptons == 1)
248  {
249  if (!AnalysisUtil::muonFilter7TeV(muons) && muons.size() == 1)
250  {
251  return;
252  }
253  cutflowEnum a, b, c;
254  if (electrons.size() == 1) a = electron_eq_4_jets, b = electron_met_gt_40, c = electron_eq_2_bjets;
255  if (muons.size() == 1) a = muon_eq_4_jets, b = muon_met_gt_40, c = muon_eq_2_bjets;
256  if (nJets == 4)
257  {
258  incrementCut(a);
259  {
260  if (Met > 40)
261  {
262  incrementCut(b);
263  if (nBJets == 2)
264  {
265  incrementCut(c);
266  }
267  }
268  }
269  }
270  }
271 
272 
273 
274 
275  // minimal selection requirements for single lepton
276  if (nLeptons == 1 && nBJets >= 2 && nLightJets >= 2 && Met > 40)
277  {
278  double mT = calcMt(metVec, leptons[0]->mom());
279 
280  auto isValidTop = [](double mean, double width, double mass) {return mass < mean - 0.5 * width;};
281 
282  if (mT > 30)
283  {
284  double mean = 0.0, width = 0.0;
285  P4 hadronicTop;
286 
287  P4 visibleSubsystem = *leptons[0] + *lightJets[0] + *lightJets[1] + *bJets[0] + *bJets[1];
288  double sqrtSsubMin = calcSqrtSSubMin(visibleSubsystem, metVec);
289  bool isOneLep = false;
290  // e-channel
291  if (electrons.size() == 1 && electrons[0]->pT() > 25)
292  {
293  mean = 168.4, width = 18.0;
294  hadronicTop = getBestHadronicTop(bJets, lightJets, *electrons[0], metVec, width, mean);
295  isOneLep = true;
296  }
297 
298  // mu-channel
299  if (muons.size() == 1 && muons[0]->pT() > 20)
300  {
301  mean = 168.2, width = 18.6;
302  hadronicTop = getBestHadronicTop(bJets, lightJets, *muons[0], metVec, width, mean);
303  isOneLep = true;
304  }
305 
306 
307  bool validTop = isValidTop(mean, width, hadronicTop.m());
308 
309  if (isOneLep)
310  {
311  varResults[varNames[mTopHad]].push_back(hadronicTop.m());
312  varResults[varNames[oneLepSqrtS]].push_back(sqrtSsubMin);
313  }
314 
315  // check if we are in the 1LSR signal region
316  if (isOneLep && validTop && sqrtSsubMin < 250)
317  {
318  num1LSR += event->weight();
319  if (electrons.size() == 1) incrementCut(electron_sr);
320  if (muons.size() == 1) incrementCut(muon_sr);
321  }
322  }
323  }
324 
325  if (nLeptons == 2 && Met > 40 && AnalysisUtil::oppositeSign(leptons[0], leptons[1]) && nJets >= 2)
326  {
327  P4 ll = *leptons[0] + *leptons[1];
328  double mll = ll.m();
329  incrementCut(twoLep_met_gt_40);
330  {
331  if (nBJets >= 1)
332  {
333  incrementCut(twoLep_gt_1_bjet);
334  if (mll < 81)
335  {
336  incrementCut(mll_lt_81);
337  }
338  if (mll < 81 && mll > 30)
339  {
340  incrementCut(mll_gt_30_lt_81);
341  }
342  }
343  }
344  }
345 
346  if (nLeptons == 2 && AnalysisUtil::oppositeSign(leptons[0], leptons[1]) && Met > 40 && nJets >= 2 && nBJets >= 1)
347  {
348  P4 ll = *leptons[0] + *leptons[1];
349  double mll = ll.m();
350 
351  P4 visibleSubsystem = *leptons[0] + *leptons[1] + *jets[0] + *jets[1];
352 
353  double sqrtSsubMin = calcSqrtSSubMin(visibleSubsystem, metVec);
354 
355  double mlljj = visibleSubsystem.m();
356  varResults[varNames[mllKey]].push_back(mll);
357  if (mll > 30 && mll < 81)
358  {
359  bool isTwoLeptonEvent = false;
360 
361  // ee channel
362  if (electrons.size() == 2 && electrons[0]->pT() > 25)
363  {
364  isTwoLeptonEvent = true;
365  }
366 
367  // mu-mu channel
368  if (muons.size() == 2 && muons[0]->pT() > 20)
369  {
370  isTwoLeptonEvent = true;
371  }
372 
373  // e-mu channel
374  if (electrons.size() == 1 && muons.size() == 1 && (electrons[0]->pT() > 25 || muons[0]->pT() > 20))
375  {
376  isTwoLeptonEvent = true;
377  }
378 
379  if (isTwoLeptonEvent)
380  {
381  varResults[varNames[twoLepSqrtS]].push_back(sqrtSsubMin);
382 
383  if (sqrtSsubMin < 225)
384  {
385  num2LSR1 += event->weight();
386  incrementCut(num_2lsr1);
387  }
388  if (sqrtSsubMin < 235 && mlljj < 140)
389  {
390  num2LSR2 += event->weight();
391  incrementCut(num_2lsr2);
392  }
393  }
394  }
395  }
396  }
397 
403 void combine(const Analysis* other)
405 {
406  const Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb* specificOther = dynamic_cast<const Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb*>(other);
407  num1LSR += specificOther->num1LSR;
408  num2LSR1 += specificOther->num2LSR1;
409  num2LSR2 += specificOther->num2LSR2;
410 }
411 
412 
413  void collect_results()
414  {
415  //saveCutFlow();
416 
417  add_result(SignalRegionData("1LSR", 50, {num1LSR, 0.}, {38., 7.}));
418  add_result(SignalRegionData("2LSR1", 123, {num2LSR1, 0.}, {115., 15.}));
419  add_result(SignalRegionData("2LSR2", 47, {num2LSR2, 0.}, {46., 7.}));
420 
421  //cout << "1LSR: " << num1LSR << ", 2LSR1: " << num2LSR1 << ", 2LSR2: " << num2LSR2 << endl;
422 
423  /*for (std::pair<std::string, std::vector<double>> entry : varResults)
424  {
425  cout << "SAVE_START:" << entry.first << endl;
426  for (double value : entry.second)
427  {
428  cout << value << endl;
429  }
430  cout << "SAVE_END" << endl;
431  }*/
432  }
433 
434 protected:
435 
436 void analysis_specific_reset()
437 {
438  num1LSR = 0;
439  num2LSR1 = 0;
440  num2LSR2 = 0;
441  for (std::string varName : varNames)
442  {
443  varResults[varName] = {};
444  }
445 }
446 
447 /*void scale(double factor)
448  {
449  HEPUtilsAnalysis::scale(factor);
450  cout << "SAVE_XSEC:" << xsec() << endl;
451  auto save = [](double value, std::string name)
452  {
453  cout << "SAVE_START:" << name << endl;
454  cout << value << endl;
455  cout << "SAVE_END" << endl;
456  };
457  save(num1LSR, "num1LSR");
458  save(num2LSR1, "num2LSR1");
459  save(num2LSR2, "num2LSR2");
460  }*/
461 
462 
463 void incrementCut(int cutIndex)
464 {
465  cutflows[cutflowNames[cutIndex]]++;
466 }
467 
468 void saveCutFlow()
469 {
470  double scale_by = 1.0;
471  cout << "SAVE_START:cuts" << endl;
472  cout << "CUT;RAW;SCALED;%" << endl;
473  double initialCut = cutflows[cutflowNames[Total_events]];
474  double thisCut;
475  for (std::string name : cutflowNames) {
476  thisCut = cutflows[name];
477  cout << name.c_str() << ";"
478  << thisCut << ";"
479  << thisCut * scale_by << ";"
480  << 100. * thisCut / initialCut
481  << endl;
482  }
483  cout << "SAVE_END" << endl;
484 }
485 };
486 
487 DEFINE_ANALYSIS_FACTORY(ATLAS_7TeV_1OR2LEPStop_4_7invfb)
488 }
489 }
490 
static std::vector< const HEPUtils::Jet * > filterPtEta(std::vector< const HEPUtils::Jet *> jets, double pT, double absEta)
Utils for ColliderBit analyses.
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
static std::vector< const HEPUtils::Jet * > jetLeptonOverlapRemoval(std::vector< const HEPUtils::Jet *> jets, std::vector< const HEPUtils::Particle *> leptons, double dR)
static bool sortParticlesByPt(const HEPUtils::Particle *particle1, const HEPUtils::Particle *particle2)
Funk delta(std::string arg, double pos, double width)
Definition: daFunk.hpp:902
static bool muonFilter7TeV(const std::vector< const HEPUtils::Particle *> &muons)
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 *)
STL namespace.
START_MODEL b
Definition: demo.hpp:270
static bool sortJetsByPt(const HEPUtils::Jet *jet1, const HEPUtils::Jet *jet2)
EXPORT_SYMBOLS double sqr(double a)
returns square of double - saves tedious repetition
static std::vector< const HEPUtils::Jet * > filterMaxEta(const std::vector< const HEPUtils::Jet *> &jets, double maxAbsEta)
The Cutflow and Cutflows classes.
const double mu
Definition: SM_Z.hpp:42
static std::vector< const HEPUtils::Particle * > getSortedLeptons(const std::vector< std::vector< const HEPUtils::Particle *>> allLeptons)
void applyTightIDElectronSelection(std::vector< const HEPUtils::Particle *> &electrons)
Efficiency function for Tight ID electrons.
static double dot2D(HEPUtils::P4 mom1, HEPUtils::P4 mom2)
static bool oppositeSign(const HEPUtils::Particle *a, const HEPUtils::Particle *b)
def combine(region_file, pip_file)
Definition: colouring.py:169
static std::vector< const HEPUtils::Particle * > leptonJetOverlapRemoval(std::vector< const HEPUtils::Particle *> leptons, std::vector< const HEPUtils::Jet *> jets, double dR)
DEFINE_ANALYSIS_FACTORY(ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb)
TODO: see if we can use this one:
Definition: Analysis.hpp:33
Class for ColliderBit analyses.
Functions that do super fast ATLAS detector simulation based on four-vector smearing.