gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
Analysis_ATLAS_13TeV_2LEPStop_36invfb.cpp
Go to the documentation of this file.
1 #include <vector>
2 #include <cmath>
3 #include <memory>
4 #include <iomanip>
5 
9 
10 using namespace std;
11 
12 /* The ATLAS 2 lepton direct stop analysis (36.1fb^-1) - `heavy stop'.
13 
14  Based on: arXiv:1708.03247
15  Author: Yang Zhang
16 
17  Known errors:
18  1. Jet overlap removal with muon is not done becasue miss trac imformation
19  2. Applied specific efficiencies on muon for validation
20 
21 */
22 
23 namespace Gambit {
24  namespace ColliderBit {
25 
26 
27  bool sortByPT_j(const HEPUtils::Jet* jet1, const HEPUtils::Jet* jet2) { return (jet1->pT() > jet2->pT()); }
28  bool sortByPT_l(const HEPUtils::Particle* lep1, const HEPUtils::Particle* lep2) { return (lep1->pT() > lep2->pT()); }
29 
30  class Analysis_ATLAS_13TeV_2LEPStop_36invfb : public Analysis {
31  private:
32 
33  // Numbers passing cuts
34  std::map<string, EventCounter> _counters = {
35  {"SRASF120", EventCounter("SRASF120")},
36  {"SRADF120", EventCounter("SRADF120")},
37  {"SRASF140", EventCounter("SRASF140")},
38  {"SRADF140", EventCounter("SRADF140")},
39  {"SRASF160", EventCounter("SRASF160")},
40  {"SRADF160", EventCounter("SRADF160")},
41  {"SRASF180", EventCounter("SRASF180")},
42  {"SRADF180", EventCounter("SRADF180")},
43  {"SRBSF120", EventCounter("SRBSF120")},
44  {"SRBDF120", EventCounter("SRBDF120")},
45  {"SRBSF140", EventCounter("SRBSF140")},
46  {"SRBDF140", EventCounter("SRBDF140")},
47  {"SRCSF110", EventCounter("SRCSF110")},
48  {"SRCDF110", EventCounter("SRCDF110")},
49  {"SR4b", EventCounter("SR4b")},
50  };
51 
52  // Cut Flow
53  vector<int> cutFlowVector;
54  vector<string> cutFlowVector_str;
55  int NCUTS;
56 
57  // // debug
58  // ofstream Savelep1;
59  // ofstream Savelep2;
60 
61  // Jet overlap removal
62  void JetLeptonOverlapRemoval(vector<const HEPUtils::Jet*> &jetvec, vector<const HEPUtils::Particle*> &lepvec, double DeltaRMax) {
63  //Routine to do jet-lepton check
64  //Discards jets if they are within DeltaRMax of a lepton
65 
66  vector<const HEPUtils::Jet*> Survivors;
67 
68  for(unsigned int itjet = 0; itjet < jetvec.size(); itjet++) {
69  bool overlap = false;
70  HEPUtils::P4 jetmom=jetvec.at(itjet)->mom();
71  for(unsigned int itlep = 0; itlep < lepvec.size(); itlep++) {
72  HEPUtils::P4 lepmom=lepvec.at(itlep)->mom();
73  double dR;
74 
75  dR=jetmom.deltaR_eta(lepmom);
76 
77  if(fabs(dR) <= DeltaRMax) overlap=true;
78  }
79  if(overlap) continue;
80  Survivors.push_back(jetvec.at(itjet));
81  }
82  jetvec=Survivors;
83 
84  return;
85  }
86 
87  // Lepton overlap removal
88  void LeptonJetOverlapRemoval(vector<const HEPUtils::Particle*> &lepvec, vector<const HEPUtils::Jet*> &jetvec, double DeltaRMax) {
89  //Routine to do lepton-jet check
90  //Discards leptons if they are within DeltaRMax of a jet
91 
92  vector<const HEPUtils::Particle*> Survivors;
93 
94  for(unsigned int itlep = 0; itlep < lepvec.size(); itlep++) {
95  bool overlap = false;
96  HEPUtils::P4 lepmom=lepvec.at(itlep)->mom();
97  for(unsigned int itjet= 0; itjet < jetvec.size(); itjet++) {
98  HEPUtils::P4 jetmom=jetvec.at(itjet)->mom();
99  double dR;
100 
101  dR=jetmom.deltaR_eta(lepmom);
102 
103  if(fabs(dR) <= DeltaRMax) overlap=true;
104  }
105  if(overlap) continue;
106  Survivors.push_back(lepvec.at(itlep));
107  }
108  lepvec=Survivors;
109 
110  return;
111  }
112 
113 
114  public:
115 
116  // Required detector sim
117  static constexpr const char* detector = "ATLAS";
118 
119  Analysis_ATLAS_13TeV_2LEPStop_36invfb() {
120 
121  set_analysis_name("ATLAS_13TeV_2LEPStop_36invfb");
122  set_luminosity(36.1);
123 
124  NCUTS= 66;
125 
126  // //debug
127  // Savelep1.open("lep1.txt");
128  // Savelep2.open("lep2.txt");
129 
130  for(int i=0;i<NCUTS;i++){
131  cutFlowVector.push_back(0);
132  cutFlowVector_str.push_back("");
133  }
134 
135  }
136 
137  void run(const HEPUtils::Event* event) {
138 
139  // Missing energy
140  double met = event->met();
141  HEPUtils::P4 ptot = event->missingmom();
142 
143  // Baseline lepton objects
144  vector<const HEPUtils::Particle*> blElectrons, blMuons; // Used for SR-2body and SR-3body
145  vector<const HEPUtils::Particle*> baselineElectrons, baselineMuons; // Used for SR-4body
146  for (const HEPUtils::Particle* electron : event->electrons()) {
147  // Same with the code snippet, not the experimental report
148  if (electron->pT() > 10. && electron->abseta() < 2.47) blElectrons.push_back(electron);
149  if (electron->pT() > 7. && electron->abseta() < 2.47) baselineElectrons.push_back(electron);
150  }
151 
152  // Apply electron efficiency
153  ATLAS::applyElectronEff(blElectrons);
154  ATLAS::applyElectronEff(baselineElectrons);
155 
156  // Apply loose electron selection
158  ATLAS::applyLooseIDElectronSelectionR2(baselineElectrons);
159 
160  const std::vector<double> a = {0,10.};
161  const std::vector<double> b = {0,10000.};
162  const vector<double> cMu={0.89};
163  HEPUtils::BinnedFn2D<double> _eff2dMu(a,b,cMu);
164  for (const HEPUtils::Particle* muon : event->muons()) {
165  bool hasTrig=has_tag(_eff2dMu, muon->abseta(), muon->pT());
166  // Same with the code snippet, not the experimental report
167  if (muon->pT() > 10. && muon->abseta() < 2.5 && hasTrig) blMuons.push_back(muon);
168  if (muon->pT() > 7. && muon->abseta() < 2.5 && hasTrig) baselineMuons.push_back(muon);
169  }
170 
171  // Apply muon efficiency
172  ATLAS::applyMuonEff(blMuons);
173  ATLAS::applyMuonEff(baselineMuons);
174 
175  // Jets
176  vector<const HEPUtils::Jet*> blJets; // Used for SR-2body and SR-3body
177  vector<const HEPUtils::Jet*> baselineJets; // Used for SR-4body
178  for (const HEPUtils::Jet* jet : event->jets()) {
179  if (jet->pT() > 20. && fabs(jet->eta()) < 2.8) blJets.push_back(jet);
180  if (jet->pT() > 20. && fabs(jet->eta()) < 2.8) baselineJets.push_back(jet);
181  }
182 
183  // Overlap removal
184  JetLeptonOverlapRemoval(blJets,blElectrons,0.2);
185  LeptonJetOverlapRemoval(blElectrons,blJets,0.4);
186  LeptonJetOverlapRemoval(blMuons,blJets,0.4);
187  JetLeptonOverlapRemoval(baselineJets,baselineElectrons,0.2);
188  LeptonJetOverlapRemoval(baselineElectrons,baselineJets,0.4);
189  LeptonJetOverlapRemoval(baselineMuons,baselineJets,0.4);
190 
191  //Signal Jet
192  vector<const HEPUtils::Jet*> sgJets; // Used for SR-2body and SR-3body
193  vector<const HEPUtils::Jet*> signalJets; // Used for SR-4body
194  for (const HEPUtils::Jet* jet : blJets) {
195  if (jet->pT() > 20. && fabs(jet->eta()) < 2.5) {
196  sgJets.push_back(jet);
197  }
198  }
199  for (const HEPUtils::Jet* jet : baselineJets) {
200  if (jet->pT() > 25. && fabs(jet->eta()) < 2.5) {
201  signalJets.push_back(jet);
202  }
203  }
204 
205 
206  //Signal Leptons
207  vector<const HEPUtils::Particle*> sgElectrons, sgLeptons; // Used for SR-2body and SR-3body
208  vector<const HEPUtils::Particle*> signalLeptons; // Used for SR-4body
209  for (const HEPUtils::Particle* electron : blElectrons) {
210  if (electron->pT() > 10. && fabs(electron->eta()) < 2.47){
211  sgElectrons.push_back(electron);
212  }
213  }
215  for (const HEPUtils::Particle* electron : sgElectrons) {
216  sgLeptons.push_back(electron);
217  }
218  for (const HEPUtils::Particle* muon : blMuons) {
219  if (muon->pT() > 10. && fabs(muon->eta()) < 2.4){
220  sgLeptons.push_back(muon);
221  }
222  }
223  ATLAS::applyMediumIDElectronSelectionR2(baselineElectrons);
224  for (const HEPUtils::Particle* electron : baselineElectrons) {
225  signalLeptons.push_back(electron);
226  }
227  for (const HEPUtils::Particle* muon : baselineMuons) {
228  signalLeptons.push_back(muon);
229  }
230 
231 
232  //Put signal jets/leptons in pT order
233  std::sort(signalJets.begin(), signalJets.end(), sortByPT_j);
234  std::sort(signalLeptons.begin(), signalLeptons.end(), sortByPT_l);
235  std::sort(sgJets.begin(), sgJets.end(), sortByPT_j);
236  std::sort(sgLeptons.begin(), sgLeptons.end(), sortByPT_l);
237 
238  // Function used to get b jets
239  vector<const HEPUtils::Jet*> sgbJets; // Used for SR-2body and SR-3body
240  vector<const HEPUtils::Jet*> sgJetsGt25; // Used for SR-2body
241  //const std::vector<double> a = {0,10.};
242  //const std::vector<double> b = {0,10000.};
243  const std::vector<double> c = {0.77};
244  HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
245  for (const HEPUtils::Jet* jet :sgJets) {
246  if (jet->pT() > 25.) {
247  sgJetsGt25.push_back(jet);
248  bool hasTag=has_tag(_eff2d, jet->abseta(), jet->pT());
249  if(jet->btag() && hasTag && jet->pT() > 25.) sgbJets.push_back(jet);
250  }
251  }
252  int nbjet = sgbJets.size(); // Used for SR-2body and SR-3body
253  int njet25= sgJetsGt25.size(); // Used for SR-2body
254  int njet = sgJets.size(); // Used for SR-2body
255 
256  // We now have the signal electrons, muons, jets and b jets- move on to the analysis
257  /*********************************************************/
258  /* */
259  /* SIGNAL REGIONS SR-2body */
260  /* */
261  /*********************************************************/
262  // Flags for SR-2body
263  bool cABC_TriggerOS =false;
264  bool cABC_SF =false;
265  // For SRA
266  bool cA_mllGt111 =false;
267  bool cA_nobjet =false;
268  bool CA_R2l2j =false;
269  bool cA_deltaX =false;
270  bool cA_MT2120 =false;
271  bool cA_MT2140 =false;
272  bool cA_MT2160 =false;
273  bool cA_MT2180 =false;
274  // For SRB
275  bool cBC_mllExMz =false;
276  bool cBC_nbnj =false;
277  bool cB_DelBoost =false;
278  bool cB_MT2120 =false;
279  bool cB_MT2140 =false;
280  // For SRC
281  bool cC_njGt2 =false;
282  bool cC_R2lGt1o2 =false;
283  bool cC_METGt200 =false;
284  bool cC_MT2110 =false;
285 
286  //Lepton Num
287  if(sgLeptons.size() == 2){
288 
289  // Opposite sign leptons, pT(l1,l2)>25,20GeV, mll>20GeV
290  HEPUtils::P4 lepton0=sgLeptons.at(0)->mom();
291  HEPUtils::P4 lepton1=sgLeptons.at(1)->mom();
292  double Mll= (lepton0+lepton1).m();
293 
294  // Savelep1 << sgLeptons[0]->pT() << endl;
295  // Savelep2 << sgLeptons[1]->pT() << endl;
296 
297  if (sgLeptons[0]->pid()*sgLeptons[1]->pid()<0. && sgLeptons[0]->pT() > 25. && sgLeptons[1]->pT() > 20. && Mll>20.){
298  cABC_TriggerOS = true;
299  if (sgLeptons[0]->pid()+sgLeptons[1]->pid()==0) cABC_SF = true;
300  /********* SRA-2body *********/
301  if (Mll>111.2) cA_mllGt111 = true;
302  if (nbjet==0) cA_nobjet = true;
303 
304  double ptJet1=0.; if ( njet >= 1) ptJet1 = sgJets.at(0)->pT();
305  double ptJet2=0.; if ( njet >= 2) ptJet2 = sgJets.at(1)->pT();
306  double R2l2j=met / (met + lepton0.pT() + lepton1.pT() + ptJet1 + ptJet2);
307  if (R2l2j>0.3) CA_R2l2j = true;
308 
309  double DX = fabs((2 * (lepton0.pz() + lepton1.pz())) / 13000.);
310  if (DX<0.07) cA_deltaX = true;
311 
312  //Calculate MT2
313  double mt2ll=0;
314  double pa_a[3] = { 0, sgLeptons[0]->mom().px(), sgLeptons[0]->mom().py() };
315  double pb_a[3] = { 0, sgLeptons[1]->mom().px(), sgLeptons[1]->mom().py() };
316  double pmiss_a[3] = { 0, ptot.px(), ptot.py() };
317  double mn_a = 0.;
318 
319  mt2_bisect::mt2 mt2_event_a;
320  mt2_event_a.set_momenta(pa_a,pb_a,pmiss_a);
321  mt2_event_a.set_mn(mn_a);
322  mt2ll = mt2_event_a.get_mt2();
323 
324  if(mt2ll>120 && mt2ll<140) cA_MT2120 = true;
325  if(mt2ll>140 && mt2ll<160) cA_MT2140 = true;
326  if(mt2ll>160 && mt2ll<180) cA_MT2160 = true;
327  if(mt2ll>180) cA_MT2180 = true;
328 
329  /********* SRB-2body *********/
330  if (Mll>111.2 || Mll<71.2) cBC_mllExMz = true;
331  if (nbjet>0 && njet25>1) cBC_nbnj = true;
332  //Calculate Delta phi_boost
333  HEPUtils::P4 pbll_TLV = lepton0 + lepton1 + ptot;
334  double dPhiEtmisspbll = fabs(pbll_TLV.deltaPhi(ptot));
335  if (dPhiEtmisspbll<1.5) cB_DelBoost = true;
336  if(mt2ll>120 && mt2ll<140) cB_MT2120 = true;
337  if(mt2ll>140 ) cB_MT2140 = true;
338 
339  /********* SRC-2body *********/
340  if (njet25>2) cC_njGt2 = true;
341  double R2L=met / (lepton0.pT() + lepton1.pT());
342  if (R2L>1.2) cC_R2lGt1o2 = true;
343  if (met>200) cC_METGt200 = true;
344  if (mt2ll>110) cC_MT2110 = true;
345  }
346  }
347 
348  /*********************************************************/
349  /* */
350  /* SIGNAL REGIONS SR-4body */
351  /* */
352  /*********************************************************/
353  // Flags
354  bool c4_METOSlepton =false;
355  bool c4_mllGt10 =false;
356  bool c4_SoftLepton =false;
357  bool c4_njetGt2 =false;
358  bool c4_Jet1PtGt150 =false;
359  bool c4_Jet3PtMET =false;
360  bool c4_R2l4j =false;
361  bool c4_R2l =false;
362  bool c4_2bjetveto =false;
363 
364  //MET trigger and Lepton Num
365  if(signalLeptons.size() == 2 && met>200 ){
366 
367  // Opposite sign leptons
368  if (signalLeptons[0]->pid()*signalLeptons[1]->pid()<0) c4_METOSlepton=true;
369 
370  // mll
371  HEPUtils::P4 lep0=signalLeptons.at(0)->mom();
372  HEPUtils::P4 lep1=signalLeptons.at(1)->mom();
373  double mll= (lep0+lep1).m();
374  if(mll>10) c4_mllGt10=true;
375  //Soft lepton
376  if(lep0.pT()<80 && lep1.pT()<35) c4_SoftLepton=true;
377 
378  //Number of jet
379  int nJets = signalJets.size();
380  if(nJets>=2){
381 
382  c4_njetGt2=true;
383  double ptJet1=signalJets.at(0)->pT();
384  double ptJet2=signalJets.at(1)->pT();
385  double ptJet3=0.;
386  if ( nJets >= 3) ptJet3 = signalJets.at(2)->pT();
387  double ptJet4=0.;
388  if ( nJets >= 4) ptJet4 = signalJets.at(3)->pT();
389 
390  //PT(j1)>150GeV
391  if (ptJet1>150) c4_Jet1PtGt150 =true;
392  //PT(j3)/MET<0.14
393  if (nJets>=3 && ptJet3/met < 0.14) c4_Jet3PtMET=true;
394  if (nJets<3) c4_Jet3PtMET=true;
395 
396  //R_{2l4j}>0.35
397  double R2l4j=met / (met + lep0.pT() + lep1.pT() + ptJet1 + ptJet2 + ptJet3 + ptJet4);
398  if (R2l4j>0.35) c4_R2l4j=true;
399  //R_{2l}>12
400  double R2l=met / (lep0.pT() + lep1.pT());
401  if (R2l>12.) c4_R2l=true;
402 
403  //veto b-jet on j1 and j2
404  // Function used to get b jets
405  /*const std::vector<double> a = {0,10.};
406  const std::vector<double> b = {0,10000.};
407  const std::vector<double> c = {0.7};
408  HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);*/
409  bool j1Tag=has_tag(_eff2d, signalJets.at(0)->abseta(), signalJets.at(0)->pT());
410  bool j2Tag=has_tag(_eff2d, signalJets.at(1)->abseta(), signalJets.at(1)->pT());
411  if (!(signalJets.at(0)->btag()&&j1Tag&&signalJets.at(1)->btag()&&j2Tag)) c4_2bjetveto=true;
412  }
413  }
414 
415  /*********************************************************/
416  /* */
417  /* Cut Flow */
418  /* */
419  /*********************************************************/
420  cutFlowVector_str[0] = "Total ";
421  /*---------------------------------------*/
422  cutFlowVector_str[1] = "SR2A--trigger && 2 OS lepton";
423  cutFlowVector_str[2] = "SR2ASF--Same flavour";
424  cutFlowVector_str[3] = "SR2ASF--mll>111GeV";
425  cutFlowVector_str[4] = "SR2ASF--n_{b-jets}=0";
426  cutFlowVector_str[5] = "SR2ASF--R_{2l2j}>0.3";
427  cutFlowVector_str[6] = "SR2ASF--Delta x<0.07";
428  cutFlowVector_str[7] = "SR2ASF--120<MT2<140";
429  cutFlowVector_str[8] = "SR2ASF--140<MT2<160";
430  cutFlowVector_str[9] = "SR2ASF--160<MT2<180";
431  cutFlowVector_str[10] = "SR2ASF--180<MT2";
432 
433  cutFlowVector_str[11] = "SR2ADF--Different falvour";
434  cutFlowVector_str[12] = "SR2ADF--mll>111GeV(only SF)";
435  cutFlowVector_str[13] = "SR2ADF--n_{b-jets}=0";
436  cutFlowVector_str[14] = "SR2ADF--R_{2l2j}>0.3(only SF)";
437  cutFlowVector_str[15] = "SR2ADF--Delta x<0.07";
438  cutFlowVector_str[16] = "SR2ADF--120<MT2<140";
439  cutFlowVector_str[17] = "SR2ADF--140<MT2<160";
440  cutFlowVector_str[18] = "SR2ADF--160<MT2<180";
441  cutFlowVector_str[19] = "SR2ADF--180<MT2";
442  /*---------------------------------------*/
443  cutFlowVector_str[20] = "SR2BC--trigger && 2 OS lepton";
444  cutFlowVector_str[21] = "SR2BSF--Same flavour";
445  cutFlowVector_str[22] = "SR2BSF--mll>111GeV or mll<71GeV";
446  cutFlowVector_str[23] = "SR2BSF--n_{b-jets}>0 && n_{jets}>1";
447  cutFlowVector_str[24] = "SR2BSF--Delta phi_{boost}<1.5";
448  cutFlowVector_str[25] = "SR2BSF--120<MT2<140";
449  cutFlowVector_str[26] = "SR2BSF--140<MT2";
450 
451  cutFlowVector_str[27] = "SR2BDF--Different flavour";
452  cutFlowVector_str[28] = "SR2BDF--mll>111GeV or mll<71GeV(only SF)";
453  cutFlowVector_str[29] = "SR2BDF--n_{b-jets}>0 && n_{jets}>1";
454  cutFlowVector_str[30] = "SR2BDF--Delta phi_{boost}<1.5";
455  cutFlowVector_str[31] = "SR2BDF--120<MT2<140";
456  cutFlowVector_str[32] = "SR2BDF--140<MT2";
457 
458  cutFlowVector_str[33] = "SR2CSF--n_{b-jets}>0 && n_{jets}>1";
459  cutFlowVector_str[34] = "SR2CSF--n_{jets}>2";
460  cutFlowVector_str[35] = "SR2CSF--R_{2l}>1.2";
461  cutFlowVector_str[36] = "SR2CSF--E_T^{miss}>200GeV";
462  cutFlowVector_str[37] = "SR2CSF--110<MT2";
463 
464  cutFlowVector_str[38] = "SR2CDF--n_{b-jets}>0 && n_{jets}>1";
465  cutFlowVector_str[39] = "SR2CDF--n_{jets}>2";
466  cutFlowVector_str[40] = "SR2CDF--R_{2l}>1.2";
467  cutFlowVector_str[41] = "SR2CDF--E_T^{miss}>200GeV";
468  cutFlowVector_str[42] = "SR2CDF--110<MT2";
469 
470  /*---------------------------------------*/
471  cutFlowVector_str[57] = "SR4b--MET trigger && 2 OS Leptons ";
472  cutFlowVector_str[58] = "SR4b--m_{ll}>10GeV ";
473  cutFlowVector_str[59] = "SR4b--PT(l1)<80GeV && PT(l1)<35GeV ";
474  cutFlowVector_str[60] = "SR4b--n_{jets}>2 ";
475  cutFlowVector_str[61] = "SR4b--PT(j1)>150GeV ";
476  cutFlowVector_str[62] = "SR4b--PT(j3)/MET<0.14 ";
477  cutFlowVector_str[63] = "SR4b--R_{2l4j}>0.35 ";
478  cutFlowVector_str[64] = "SR4b--R_{2l}>12 ";
479  cutFlowVector_str[65] = "SR4b--veto on j1 and j2 ";
480 
481  for(int j=0;j<NCUTS;j++){
482  if(
483  (j==0) ||
484  /********* SRA-2body *********/
485  (j==1 && cABC_TriggerOS) ||
486  //Same Flavour
487  (j==2 && cABC_SF) ||
488  (j==3 && cABC_SF && cA_mllGt111) ||
489  (j==4 && cABC_SF && cA_mllGt111 && cA_nobjet) ||
490  (j==5 && cABC_SF && cA_mllGt111 && cA_nobjet && CA_R2l2j) ||
491  (j==6 && cABC_SF && cA_mllGt111 && cA_nobjet && CA_R2l2j && cA_deltaX) ||
492  (j==7 && cABC_SF && cA_mllGt111 && cA_nobjet && CA_R2l2j && cA_deltaX && cA_MT2120) ||
493  (j==8 && cABC_SF && cA_mllGt111 && cA_nobjet && CA_R2l2j && cA_deltaX && cA_MT2140) ||
494  (j==9 && cABC_SF && cA_mllGt111 && cA_nobjet && CA_R2l2j && cA_deltaX && cA_MT2160) ||
495  (j==10 && cABC_SF && cA_mllGt111 && cA_nobjet && CA_R2l2j && cA_deltaX && cA_MT2180) ||
496  //different Flavour
497  (j==11 && cABC_TriggerOS && (!cABC_SF)) ||
498  (j==12 && cABC_TriggerOS && (!cABC_SF)) ||
499  (j==13 && (!cABC_SF) && cA_nobjet) ||
500  (j==14 && (!cABC_SF) && cA_nobjet) ||
501  (j==15 && (!cABC_SF) && cA_nobjet && cA_deltaX) ||
502  (j==16 && (!cABC_SF) && cA_nobjet && cA_deltaX && cA_MT2120) ||
503  (j==17 && (!cABC_SF) && cA_nobjet && cA_deltaX && cA_MT2140) ||
504  (j==18 && (!cABC_SF) && cA_nobjet && cA_deltaX && cA_MT2160) ||
505  (j==19 && (!cABC_SF) && cA_nobjet && cA_deltaX && cA_MT2180) ||
506  /********* SRB-2body *********/
507  (j==20 && cABC_TriggerOS) ||
508  //Same Flavour
509  (j==21 && cABC_SF) ||
510  (j==22 && cABC_SF && cBC_mllExMz) ||
511  (j==23 && cABC_SF && cBC_mllExMz && cBC_nbnj) ||
512  (j==24 && cABC_SF && cBC_mllExMz && cBC_nbnj && cB_DelBoost) ||
513  (j==25 && cABC_SF && cBC_mllExMz && cBC_nbnj && cB_DelBoost && cB_MT2120) ||
514  (j==26 && cABC_SF && cBC_mllExMz && cBC_nbnj && cB_DelBoost && cB_MT2140) ||
515  //different Flavour
516  (j==27 && cABC_TriggerOS && (!cABC_SF)) ||
517  (j==28 && cABC_TriggerOS && (!cABC_SF)) ||
518  (j==29 && (!cABC_SF) && cBC_nbnj) ||
519  (j==30 && (!cABC_SF) && cBC_nbnj && cB_DelBoost) ||
520  (j==31 && (!cABC_SF) && cBC_nbnj && cB_DelBoost && cB_MT2120) ||
521  (j==32 && (!cABC_SF) && cBC_nbnj && cB_DelBoost && cB_MT2140) ||
522  /********* SRC-2body *********/
523  //Same Flavour
524  (j==33 && cABC_SF && cBC_mllExMz && cBC_nbnj) ||
525  (j==34 && cABC_SF && cBC_mllExMz && cBC_nbnj && cC_njGt2) ||
526  (j==35 && cABC_SF && cBC_mllExMz && cBC_nbnj && cC_njGt2 && cC_R2lGt1o2) ||
527  (j==36 && cABC_SF && cBC_mllExMz && cBC_nbnj && cC_njGt2 && cC_R2lGt1o2 && cC_METGt200) ||
528  (j==37 && cABC_SF && cBC_mllExMz && cBC_nbnj && cC_njGt2 && cC_R2lGt1o2 && cC_METGt200 && cC_MT2110) ||
529  //different Flavour
530  (j==38 && (!cABC_SF) && cBC_nbnj) ||
531  (j==39 && (!cABC_SF) && cBC_nbnj && cC_njGt2) ||
532  (j==40 && (!cABC_SF) && cBC_nbnj && cC_njGt2 && cC_R2lGt1o2) ||
533  (j==41 && (!cABC_SF) && cBC_nbnj && cC_njGt2 && cC_R2lGt1o2 && cC_METGt200) ||
534  (j==42 && (!cABC_SF) && cBC_nbnj && cC_njGt2 && cC_R2lGt1o2 && cC_METGt200 && cC_MT2110) ||
535  /********* SR-4body *********/
536  (j==57 && c4_METOSlepton ) ||
537  (j==58 && c4_METOSlepton && c4_mllGt10) ||
538  (j==59 && c4_METOSlepton && c4_mllGt10 && c4_SoftLepton) ||
539  (j==60 && c4_METOSlepton && c4_mllGt10 && c4_SoftLepton && c4_njetGt2) ||
540  (j==61 && c4_METOSlepton && c4_mllGt10 && c4_SoftLepton && c4_Jet1PtGt150) ||
541  (j==62 && c4_METOSlepton && c4_mllGt10 && c4_SoftLepton && c4_Jet1PtGt150 && c4_Jet3PtMET) ||
542  (j==63 && c4_METOSlepton && c4_mllGt10 && c4_SoftLepton && c4_Jet1PtGt150 && c4_Jet3PtMET && c4_R2l4j) ||
543  (j==64 && c4_METOSlepton && c4_mllGt10 && c4_SoftLepton && c4_Jet1PtGt150 && c4_Jet3PtMET && c4_R2l4j && c4_R2l) ||
544  (j==65 && c4_METOSlepton && c4_mllGt10 && c4_SoftLepton && c4_Jet1PtGt150 && c4_Jet3PtMET && c4_R2l4j && c4_R2l && c4_2bjetveto)
545  )cutFlowVector[j]++;
546 
547  }
548  // signal region
549 
550  if ( cABC_SF && cA_mllGt111 && cA_nobjet && CA_R2l2j && cA_deltaX && cA_MT2120 ) _counters.at("SRASF120").add_event(event);
551  if ( (!cABC_SF) && cA_nobjet && cA_deltaX && cA_MT2120 ) _counters.at("SRADF120").add_event(event);
552  if ( cABC_SF && cA_mllGt111 && cA_nobjet && CA_R2l2j && cA_deltaX && cA_MT2140 ) _counters.at("SRASF140").add_event(event);
553  if ( (!cABC_SF) && cA_nobjet && cA_deltaX && cA_MT2140 ) _counters.at("SRADF140").add_event(event);
554  if ( cABC_SF && cA_mllGt111 && cA_nobjet && CA_R2l2j && cA_deltaX && cA_MT2160 ) _counters.at("SRASF160").add_event(event);
555  if ( (!cABC_SF) && cA_nobjet && cA_deltaX && cA_MT2160 ) _counters.at("SRADF160").add_event(event);
556  if ( cABC_SF && cA_mllGt111 && cA_nobjet && CA_R2l2j && cA_deltaX && cA_MT2180 ) _counters.at("SRASF180").add_event(event);
557  if ( (!cABC_SF) && cA_nobjet && cA_deltaX && cA_MT2180 ) _counters.at("SRADF180").add_event(event);
558 
559  if ( cABC_SF && cBC_mllExMz && cBC_nbnj && cB_DelBoost && cB_MT2120 ) _counters.at("SRBSF120").add_event(event);
560  if ( (!cABC_SF) && cBC_nbnj && cB_DelBoost && cB_MT2120 ) _counters.at("SRBDF120").add_event(event);
561  if ( cABC_SF && cBC_mllExMz && cBC_nbnj && cB_DelBoost && cB_MT2140 ) _counters.at("SRBSF140").add_event(event);
562  if ( (!cABC_SF) && cBC_nbnj && cB_DelBoost && cB_MT2140 ) _counters.at("SRBDF140").add_event(event);
563 
564  if ( cABC_SF && cBC_mllExMz && cBC_nbnj && cC_njGt2 && cC_R2lGt1o2 && cC_METGt200 && cC_MT2110 ) _counters.at("SRCSF110").add_event(event);
565  if ( (!cABC_SF) && cBC_nbnj && cC_njGt2 && cC_R2lGt1o2 && cC_METGt200 && cC_MT2110 ) _counters.at("SRCDF110").add_event(event);
566 
567 
568  if (c4_METOSlepton && c4_mllGt10 && c4_SoftLepton && c4_Jet1PtGt150 && c4_Jet3PtMET && c4_R2l4j && c4_R2l && c4_2bjetveto) _counters.at("SR4b").add_event(event);
569  return;
570 
571  }
572 
574  void combine(const Analysis* other)
575  {
576  const Analysis_ATLAS_13TeV_2LEPStop_36invfb* specificOther
577  = dynamic_cast<const Analysis_ATLAS_13TeV_2LEPStop_36invfb*>(other);
578 
579  for (auto& pair : _counters) { pair.second += specificOther->_counters.at(pair.first); }
580 
581  if (NCUTS != specificOther->NCUTS) NCUTS = specificOther->NCUTS;
582  for (int j=0; j<NCUTS; j++)
583  {
584  cutFlowVector[j] += specificOther->cutFlowVector[j];
585  cutFlowVector_str[j] = specificOther->cutFlowVector_str[j];
586  }
587  }
588 
589 
590  void collect_results() {
591 
592  // double scale_by=1.;
593  // cout << "------------------------------------------------------------------------------------------------------------------------------ "<<endl;
594  // cout << "CUT FLOW: ATLAS 13 TeV 2 lep stop paper "<<endl;
595  // cout << "------------------------------------------------------------------------------------------------------------------------------"<<endl;
596  // cout<< right << setw(40) << "CUT" << "," << setw(20) << "RAW" << "," << setw(20) << "SCALED"
597  // << "," << setw(20) << "%" << "," << setw(20) << "clean adj RAW"<< "," << setw(20) << "clean adj %" << endl;
598  // for (int j=0; j<NCUTS; j++) {
599  // cout << right << setw(40) << cutFlowVector_str[j].c_str() << "," << setw(20)
600  // << cutFlowVector[j] << "," << setw(20) << cutFlowVector[j]*scale_by << "," << setw(20)
601  // << 100.*cutFlowVector[j]/cutFlowVector[0] << "%" << "," << setw(20)
602  // << cutFlowVector[j]*scale_by << "," << setw(20) << 100.*cutFlowVector[j]/cutFlowVector[0]<< "%" << endl;
603  // }
604  // cout << "------------------------------------------------------------------------------------------------------------------------------ "<<endl;
605 
606 
607  // add_result(SignalRegionData("SR label", n_obs, {n_sig_MC, n_sig_MC_sys}, {n_bkg, n_bkg_err}));
608 
609  // signal regin 2-body A
610  add_result(SignalRegionData(_counters.at("SRASF120"), 22., { 20.0, 4.6}));
611  add_result(SignalRegionData(_counters.at("SRADF120"), 27., { 23.8, 4.2}));
612  add_result(SignalRegionData(_counters.at("SRASF140"), 6., { 11.0, 2.5}));
613  add_result(SignalRegionData(_counters.at("SRADF140"), 6., { 10.8, 2.1}));
614  add_result(SignalRegionData(_counters.at("SRASF160"), 10., { 5.6, 1.8}));
615  add_result(SignalRegionData(_counters.at("SRADF160"), 7., { 6.4, 1.3}));
616  add_result(SignalRegionData(_counters.at("SRASF180"), 16., { 12.3, 2.6}));
617  add_result(SignalRegionData(_counters.at("SRADF180"), 8., { 5.4, 1.7}));
618  // signal regin 2-body B
619  add_result(SignalRegionData(_counters.at("SRBSF120"), 17., { 16.3, 6.2}));
620  add_result(SignalRegionData(_counters.at("SRBDF120"), 13., { 16.1, 5.3}));
621  add_result(SignalRegionData(_counters.at("SRBSF140"), 9., { 7.4, 1.1}));
622  add_result(SignalRegionData(_counters.at("SRBDF140"), 7., { 4.8, 1.0}));
623  // signal regin 2-body C
624  add_result(SignalRegionData(_counters.at("SRCSF110"), 11., { 5.3, 1.8}));
625  add_result(SignalRegionData(_counters.at("SRCDF110"), 7., { 3.8, 1.5}));
626  // signal regin 4-body
627  add_result(SignalRegionData(_counters.at("SR4b"), 30., { 28., 6.}));
628 
629  return;
630  }
631 
632  protected:
633  void analysis_specific_reset() {
634  for (auto& pair : _counters) { pair.second.reset(); }
635  std::fill(cutFlowVector.begin(), cutFlowVector.end(), 0);
636  }
637 
638  };
639 
640 
641 
642  DEFINE_ANALYSIS_FACTORY(ATLAS_13TeV_2LEPStop_36invfb)
643 
644 
645  }
646 }
bool sortByPT_l(const HEPUtils::Particle *lep1, const HEPUtils::Particle *lep2)
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 DEFINE_ANALYSIS_FACTORY(ANAME)
For analysis factory function definition.
Definition: Analysis.hpp:124
void applyLooseIDElectronSelectionR2(std::vector< const HEPUtils::Particle *> &electrons)
Efficiency function for Loose ID electrons.
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
void set_momenta(double *pa0, double *pb0, double *pmiss0)
Definition: mt2_bisect.cpp:62
A simple class for counting events of type HEPUtils::Event.
A class for collider analyses within ColliderBit.
Definition: Analysis.hpp:41
double get_mt2()
Definition: mt2_bisect.cpp:50
A simple container for the result of one signal region from one analysis.
void applyMuonEff(std::vector< const HEPUtils::Particle *> &muons)
Randomly filter the supplied particle list by parameterised muon efficiency.
void set_mn(double mn)
Definition: mt2_bisect.cpp:130
def combine(region_file, pip_file)
Definition: colouring.py:169
bool sortByPT_j(const HEPUtils::Jet *jet1, const HEPUtils::Jet *jet2)
void applyElectronEff(std::vector< const HEPUtils::Particle *> &electrons)
Randomly filter the supplied particle list by parameterised electron efficiency.
void applyMediumIDElectronSelectionR2(std::vector< const HEPUtils::Particle *> &electrons)
Efficiency function for Loose ID electrons.
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.