gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
Analysis_ATLAS_8TeV_0LEPStop_20invfb.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 0 lepton direct stop analysis (20fb^-1) - `heavy stop'.
13 
14  Based originally on: ATLAS-CONF-2013-024
15  Now updated to arXiv:1406.1122
16  Code by Martin White, Sky French.
17  Known errors:
18  a) The isolation is already applied in the simulation rather than after overlap removal -> the electron and muon veto technically require a veto on base-line electrons/muons not overlapping with jets
19  b) ETmiss-track cuts ignored...
20 
21  Known features:
22  a) Have removed the tau veto since it disagrees with ATLAS cutflow (we lose ~2/3 of the events whilst ATLAS see a minimal effect. So ditching it is fine. In fact the implementation below is incorrect since it does not take into account the number of tracks associated to the jets.
23  b) Using mt2 bisect method from H. Cheng, Z. Han, arXiv:0810.5178 for mT2 calculation
24 
25 */
26 
27 namespace Gambit {
28  namespace ColliderBit {
29 
30  bool sortByPT(const HEPUtils::Jet* jet1, const HEPUtils::Jet* jet2) { return (jet1->pT() > jet2->pT()); }
31 
32  class Analysis_ATLAS_8TeV_0LEPStop_20invfb : public Analysis {
33  private:
34 
35  // Numbers passing cuts
36  double _numSRA1, _numSRA2, _numSRA3, _numSRA4;
37  double _numSRC1, _numSRC2, _numSRC3;
38 
39  vector<int> cutFlowVector;
40  vector<string> cutFlowVector_str;
41  int NCUTS; //=16;
42 
43  // Debug histos
44 
45  void JetLeptonOverlapRemoval(vector<const HEPUtils::Jet*> &jetvec, vector<const HEPUtils::Particle*> &lepvec, double DeltaRMax) {
46  //Routine to do jet-lepton check
47  //Discards jets if they are within DeltaRMax of a lepton
48 
49  vector<const HEPUtils::Jet*> Survivors;
50 
51  for(unsigned int itjet = 0; itjet < jetvec.size(); itjet++) {
52  bool overlap = false;
53  HEPUtils::P4 jetmom=jetvec.at(itjet)->mom();
54  for(unsigned int itlep = 0; itlep < lepvec.size(); itlep++) {
55  HEPUtils::P4 lepmom=lepvec.at(itlep)->mom();
56  double dR;
57 
58  dR=jetmom.deltaR_eta(lepmom);
59 
60  if(fabs(dR) <= DeltaRMax) overlap=true;
61  }
62  if(overlap) continue;
63  Survivors.push_back(jetvec.at(itjet));
64  }
65  jetvec=Survivors;
66 
67  return;
68  }
69 
70  void LeptonJetOverlapRemoval(vector<const HEPUtils::Particle*> &lepvec, vector<const HEPUtils::Jet*> &jetvec, double DeltaRMax) {
71  //Routine to do lepton-jet check
72  //Discards leptons if they are within DeltaRMax of a jet
73 
74  vector<const HEPUtils::Particle*> Survivors;
75 
76  for(unsigned int itlep = 0; itlep < lepvec.size(); itlep++) {
77  bool overlap = false;
78  HEPUtils::P4 lepmom=lepvec.at(itlep)->mom();
79  for(unsigned int itjet= 0; itjet < jetvec.size(); itjet++) {
80  HEPUtils::P4 jetmom=jetvec.at(itjet)->mom();
81  double dR;
82 
83  dR=jetmom.deltaR_eta(lepmom);
84 
85  if(fabs(dR) <= DeltaRMax) overlap=true;
86  }
87  if(overlap) continue;
88  Survivors.push_back(lepvec.at(itlep));
89  }
90  lepvec=Survivors;
91 
92  return;
93  }
94 
95 
96  public:
97 
98  // Required detector sim
99  static constexpr const char* detector = "ATLAS";
100 
101  Analysis_ATLAS_8TeV_0LEPStop_20invfb() {
102 
103  set_analysis_name("ATLAS_8TeV_0LEPStop_20invfb");
104  set_luminosity(20.1);
105 
106  _numSRA1 = 0 ; _numSRA2 = 0; _numSRA3 = 0; _numSRA4 = 0;
107  _numSRC1 = 0 ; _numSRC2 = 0; _numSRC3 = 0;
108  NCUTS=23;
109 
110  for(int i=0;i<NCUTS;i++){
111  cutFlowVector.push_back(0);
112  cutFlowVector_str.push_back("");
113  }
114 
115  }
116 
117 
118 
119  void run(const HEPUtils::Event* event) {
120 
121  // Missing energy
122  HEPUtils::P4 ptot = event->missingmom();
123  double met = event->met();
124 
125 
126  // Baseline lepton objects
127  vector<const HEPUtils::Particle*> baselineElectrons, baselineMuons, baselineTaus;
128  for (const HEPUtils::Particle* electron : event->electrons()) {
129  if (electron->pT() > 10. && electron->abseta() < 2.47) baselineElectrons.push_back(electron);
130  }
131 
132  // Apply electron efficiency
133  ATLAS::applyElectronEff(baselineElectrons);
134 
135  for (const HEPUtils::Particle* muon : event->muons()) {
136  if (muon->pT() > 10. && muon->abseta() < 2.4) baselineMuons.push_back(muon);
137  }
138 
139  // Apply muon efficiency
140  ATLAS::applyMuonEff(baselineMuons);
141 
142  for (const HEPUtils::Particle* tau : event->taus()) {
143  if (tau->pT() > 10. && tau->abseta() < 2.47) baselineTaus.push_back(tau);
144  }
145  ATLAS::applyTauEfficiencyR1(baselineTaus);
146 
147 
148  // Jets
149  vector<const HEPUtils::Jet*> baselineJets;
150  vector<const HEPUtils::Jet*> bJets;
151  vector<const HEPUtils::Jet*> nonBJets;
152  vector<const HEPUtils::Jet*> trueBJets; //for debugging
153 
154  // Get b jets
156  const std::vector<double> a = {0,10.};
157  const std::vector<double> b = {0,10000.};
158  const std::vector<double> c = {0.7};
159  HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
160  for (const HEPUtils::Jet* jet : event->jets()) {
161  if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) baselineJets.push_back(jet);
162  bool hasTag=has_tag(_eff2d, fabs(jet->eta()), jet->pT());
163  if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) {
164  if(jet->btag() && hasTag && fabs(jet->eta()) < 2.5 && jet->pT() > 20.){
165  bJets.push_back(jet);
166  } else {
167  nonBJets.push_back(jet);
168  }
169  }
170  }
171 
172 
173  // Overlap removal
174  vector<const HEPUtils::Particle*> signalElectrons;
175  vector<const HEPUtils::Particle*> signalMuons;
176  vector<const HEPUtils::Particle*> electronsForVeto;
177  vector<const HEPUtils::Particle*> muonsForVeto;
178 
179  vector<const HEPUtils::Jet*> signalJets;
180  vector<const HEPUtils::Jet*> signalBJets;
181  vector<const HEPUtils::Jet*> signalNonBJets;
182 
183  JetLeptonOverlapRemoval(nonBJets,baselineElectrons,0.2);
184  LeptonJetOverlapRemoval(baselineElectrons,nonBJets,0.4);
185  LeptonJetOverlapRemoval(baselineElectrons,bJets,0.4);
186  LeptonJetOverlapRemoval(baselineMuons,nonBJets,0.4);
187  LeptonJetOverlapRemoval(baselineMuons,bJets,0.4);
188 
189  for (const HEPUtils::Jet* jet : bJets) {
190  if (jet->pT() > 35. && fabs(jet->eta()) < 2.8) {
191  signalJets.push_back(jet);
192  signalBJets.push_back(jet);
193  }
194  }
195 
196  for (const HEPUtils::Jet* jet : nonBJets) {
197  if (jet->pT() > 35. && fabs(jet->eta()) < 2.8){
198  signalJets.push_back(jet);
199  signalNonBJets.push_back(jet);
200  }
201  }
202 
203  //Put signal jets in pT order
204  std::sort(signalJets.begin(), signalJets.end(), sortByPT);
205  std::sort(signalBJets.begin(), signalBJets.end(), sortByPT);
206  std::sort(signalNonBJets.begin(), signalNonBJets.end(), sortByPT);
207 
208  for (const HEPUtils::Particle* electron : baselineElectrons) {
209  signalElectrons.push_back(electron);
210  }
211 
212  for (const HEPUtils::Particle* muon : baselineMuons) {
213  signalMuons.push_back(muon);
214  }
215 
216  // We now have the signal electrons, muons, jets and b jets- move on to the analysis
217 
218  // Calculate common variables and cuts first
219  int nElectrons = signalElectrons.size();
220  int nMuons = signalMuons.size();
221  int nJets = signalJets.size();
222 
223  //Lepton veto
224  bool cut_LeptonVeto=true;
225  if((nElectrons + nMuons)>0.)cut_LeptonVeto=false;
226 
227  //Calculate dphi(jet,met) for the three leading jets
228  bool cut_dPhiJets=false;
229  bool cut_dPhiJet3=false;
230  bool cut_dPhiJet2=false;
231  bool cut_dPhiJet1=false;
232  double dphi_jetmet1=9999;
233  if(nJets>0)dphi_jetmet1=std::acos(std::cos(signalJets.at(0)->phi()-ptot.phi()));
234  double dphi_jetmet2=9999;
235  if(nJets>1)dphi_jetmet2=std::acos(std::cos(signalJets.at(1)->phi()-ptot.phi()));
236  double dphi_jetmet3=9999;
237  if(nJets>2)dphi_jetmet3=std::acos(std::cos(signalJets.at(2)->phi()-ptot.phi()));
238 
239  if(dphi_jetmet3>3.14/fabs(5.0))cut_dPhiJet3=true;
240  if(dphi_jetmet2>3.14/fabs(5.0))cut_dPhiJet2=true;
241  if(dphi_jetmet1>3.14/fabs(5.0))cut_dPhiJet1=true;
242  if(cut_dPhiJet1 && cut_dPhiJet2 && cut_dPhiJet3)cut_dPhiJets=true;
243 
244  //Number of b jets
245  bool passBJetCut=false;
246  if(signalBJets.size()>=2)passBJetCut=true;
247  //MET > 150 GeV
248  bool cut_METGt150=false;
249  if(met>150.)cut_METGt150=true;
250 
251  //Calculate dphi(b,met) for the closest b-jet in phi to MET
252  double dphi_bjetmet_min=9999.;
253  double minphi =9999.;
254  int whichb=0;
255  for(size_t j=0; j<signalBJets.size(); j++) {
256  if(fabs(std::acos(std::cos(signalBJets.at(j)->phi()-ptot.phi())))<minphi) {
257  minphi = fabs(std::acos(std::cos(signalBJets.at(j)->phi()-ptot.phi())));
258  dphi_bjetmet_min = minphi;
259  whichb=j;
260 
261  }
262  }
263 
264 
265  double mT_bjetmet_min = 0;
266  if(passBJetCut) mT_bjetmet_min = sqrt(2*signalBJets.at(whichb)->pT()*met*(1-std::cos(dphi_bjetmet_min)));
267 
268  bool cut_mTbjetmetGt175=false;
269  if(mT_bjetmet_min>175.)cut_mTbjetmetGt175=true;
270 
271  //Calculate dphi(b,met) for the furthest b-jet in phi to MET
272  double dphi_bjetmet_max=0.;
273  double maxphi =0.;
274  int whichb_max=0;
275  for(size_t j=0; j<signalBJets.size(); j++) {
276 
277  if(fabs(std::acos(std::cos(signalBJets.at(j)->phi()-ptot.phi())))>maxphi) {
278  maxphi = fabs(std::acos(std::cos(signalBJets.at(j)->phi()-ptot.phi())));
279  dphi_bjetmet_max = maxphi;
280  whichb_max=j;
281  }
282  }
283 
284  double mT_bjetmet_max = 0;
285 
286  if(passBJetCut) mT_bjetmet_max = sqrt(2*signalBJets.at(whichb_max)->pT()*met*(1-std::cos(dphi_bjetmet_max)));
287 
288  //Common preselection for all signal regions in the fully resolved case
289  bool passJetCutSRA=false;
290 
291  if(nJets>=6){
292  if(signalJets[0]->pT() > 80.
293  && signalJets[1]->pT() > 80.
294  && signalJets[2]->pT() > 35.
295  && signalJets[3]->pT() > 35.
296  && signalJets[4]->pT() > 35.
297  && signalJets[5]->pT() > 35.)passJetCutSRA=true;
298  }
299 
300  //mjjj combinations
301 
302  HEPUtils::P4 mbjj0, mbjj1;
303 
304  double mindphi_12 = 9999.;
305 
306  HEPUtils::P4 W1;
307  HEPUtils::P4 W2;
308  HEPUtils::P4 T1;
309  HEPUtils::P4 T2;
310  HEPUtils::P4 jet1;
311  HEPUtils::P4 jet2;
312  HEPUtils::P4 jet3;
313  HEPUtils::P4 jet4;
314  HEPUtils::P4 jet5;
315  HEPUtils::P4 jet6;
316 
317  //Need to form top quark four vectors from jets
318  //Use the two leading b jets as the b jets (a slight departure from ATLAS which uses the two jets with the highest b weight)
319 
320  vector<const HEPUtils::Jet*> selectBJets;
321  vector<const HEPUtils::Jet*> selectNonBJets;
322  int bjetcount=0;
323 
324  for (const HEPUtils::Jet* jet : signalBJets) {
325  if(bjetcount<2){
326  bjetcount++;
327  selectBJets.push_back(jet);
328  }
329  }
330 
331  //Now take any remaining jets in b jet collection plus the non b jets
332  int i=0;
333  for (const HEPUtils::Jet* jet : signalBJets) {
334  i++;
335  if(i>2)selectNonBJets.push_back(jet);
336  }
337 
338  for (const HEPUtils::Jet* jet : signalNonBJets){
339  selectNonBJets.push_back(jet);
340  }
341 
342  if(nJets>=6 and bjetcount==2) {
343  unsigned int j1 = 0 ; unsigned int j2 = 0; //unsigned int j4 = 0; unsigned int j5 = 0; //int j6 = 0;
344  unsigned int b1 = 0;
345  for(unsigned int k=0; k<selectNonBJets.size(); k++) {
346  for(unsigned int l=k+1; l<selectNonBJets.size(); l++) {
347  jet1.setXYZE(selectNonBJets[k]->mom().px(),selectNonBJets[k]->mom().py(),selectNonBJets[k]->mom().pz(),selectNonBJets[k]->E());
348  jet2.setXYZE(selectNonBJets[l]->mom().px(),selectNonBJets[l]->mom().py(),selectNonBJets[l]->mom().pz(),selectNonBJets[l]->E());
349 
350  if(jet1.deltaR_eta(jet2)<mindphi_12) {
351  j1 = k;
352  j2 = l;
353  mindphi_12 = jet1.deltaR_eta(jet2);
354  W1 = jet1+jet2;
355 
356  }
357 
358  }
359  }
360 
361 
362  double mindphi_w1j3 = 9999.;
363  for(unsigned int p=0; p<selectBJets.size(); p++) {
364 
365  jet3.setXYZE(selectBJets[p]->mom().px(),selectBJets[p]->mom().py(),selectBJets[p]->mom().pz(),selectBJets[p]->E());
366  if(jet3.deltaR_eta(W1)<mindphi_w1j3) {
367  b1 = p;
368  mindphi_w1j3 = jet3.deltaR_eta(W1);
369  T1 = W1+jet3;
370 
371  }
372  }
373 
374  double mindphi_45 = 9999.;
375  for(unsigned int k=0; k<selectNonBJets.size(); k++) {
376  for(unsigned int l=k; l<selectNonBJets.size(); l++) {
377  if(k!=j1 && k!=j2 && l!=j1 && l!=j2) {
378 
379  jet4.setXYZE(selectNonBJets[k]->mom().px(),selectNonBJets[k]->mom().py(),selectNonBJets[k]->mom().pz(),selectNonBJets[k]->E());
380  jet5.setXYZE(selectNonBJets[l]->mom().px(),selectNonBJets[l]->mom().py(),selectNonBJets[l]->mom().pz(),selectNonBJets[l]->E());
381 
382  if(jet4.deltaR_eta(jet5)<mindphi_45) {
383  //j4 = k;
384  //j5 = l;
385  mindphi_45 = jet4.deltaR_eta(jet5);
386  W2 = jet4+jet5;
387 
388  }
389  }
390  }
391  }
392  double mindphi_w2j6 = 9999.;
393  for(unsigned int p=0; p<selectBJets.size(); p++) {
394  if(p!=b1) {
395 
396  jet6.setXYZE(selectBJets[p]->mom().px(),selectBJets[p]->mom().py(),selectBJets[p]->mom().pz(),selectBJets[p]->E());
397 
398  if(jet6.deltaR_eta(W2)<mindphi_w2j6) {
399  //j6 = p;
400  mindphi_w2j6 = jet6.deltaR_eta(W2);
401  T2 = W2+jet6;
402 
403  }
404  }
405  }
406  mbjj0 = T1;
407  mbjj1 = T2;
408  }
409 
410  bool cut_tau=true;
411 
412  /*
413  //Tau Veto
414  for (int j=0; j<signalNonBJets.size(); j++) {
415  if(std::acos(std::cos(signalNonBJets.at(j)->phi()-ptot.phi()))<0.2*3.14)
416  cut_tau=false;
417  }*/
418 
419  //Cutflow flags
420  //bool cut_mjjj0=false;
421  //bool cut_mjjj1=false;
422 
423  bool cut_6jets=false;
424  bool cut_Btag=false;
425 
426  bool cut_METGt130=false;
427  //bool cut_METGt200=false;
428  bool cut_METGt250=false;
429  bool cut_METGt300=false;
430  bool cut_METGt350=false;
431 
432 
433  if(passJetCutSRA)cut_6jets=true;
434  if(passBJetCut)cut_Btag=true;
435  if(met>130.)cut_METGt130=true;
436  if(met>250.)cut_METGt250=true;
437  if(met>300.)cut_METGt300=true;
438  if(met>350.)cut_METGt350=true;
439  //if(nJets>=6) {
440  //if(mbjj0.m()<270 && mjjj0.m()>80) cut_mjjj0=true;
441  //if(mjjj1.m()<270 && mjjj1.m()>80) cut_mjjj1=true;
442  //}
443 
444  //Calculate min transverse mass between signal jets and ptmiss
445 
446  double mtMin=9999;
447  for(const HEPUtils::Jet* jet : signalJets){
448  double dphi_jetmet=std::acos(std::cos(jet->phi()-ptot.phi()));
449  double mT=sqrt(2.*jet->pT()*met*(1. - cos(dphi_jetmet)));
450  if(mT<mtMin)mtMin=mT;
451  }
452 
453  bool isSRA1=false;
454  bool isSRA2=false;
455  bool isSRA3=false;
456  bool isSRA4=false;
457 
458  if(cut_LeptonVeto && cut_Btag && cut_METGt150 && cut_dPhiJets && cut_mTbjetmetGt175 && cut_6jets && mbjj0.m() < 225. && mbjj1.m() < 250. && cut_tau && cut_METGt150)isSRA1=true;
459 
460  if(cut_LeptonVeto && cut_Btag && cut_METGt150 && cut_dPhiJets && cut_mTbjetmetGt175 && cut_6jets && mbjj0.m() < 225. && mbjj1.m() < 250. && cut_tau && cut_METGt250)isSRA2=true;
461 
462  if(cut_LeptonVeto && cut_Btag && cut_METGt150 && cut_dPhiJets && cut_mTbjetmetGt175 && cut_6jets && mbjj0.m() > 50. && mbjj0.m() < 250. && mbjj1.m() > 50. && mbjj1.m() < 400. && mtMin > 50. && cut_tau && cut_METGt300)isSRA3=true;
463 
464  if(cut_LeptonVeto && cut_Btag && cut_METGt150 && cut_dPhiJets && cut_mTbjetmetGt175 && cut_6jets && mbjj0.m() > 50. && mbjj0.m() < 250. && mbjj1.m() > 50. && mbjj1.m() < 400. && mtMin > 50. && cut_tau && cut_METGt350)isSRA4=true;
465 
466  //Now do the mixed regions
467 
468  //Find highest pT b jet
469  //Should no longe be necessary due to sorting of b jet collection
470  /*double leadBJetPt=0;
471  double leadBJetID=0;
472 
473  if(passBJetCut){
474  for(int j=0; j<nJets; j++) {
475  if(signalJets[j]->btag()&&signalJets[j]->pT()>leadBJetPt){
476  leadBJetPt=signalJets[j]->pT();
477  leadBJetID=j;
478  }
479  }
480  }
481 
482  //Find sub-leading pT b jet
483  double subBJetPt=0;
484  double subBJetID=0;
485 
486  if(passBJetCut){
487  for(int j=0; j<nJets; j++) {
488  if(signalJets[j]->btag()&&signalJets[j]->pT()>subBJetPt && j!=leadBJetPt){
489  subBJetPt=signalJets[j]->pT();
490  subBJetID=j;
491  }
492  }
493  }*/
494 
495  //Work out dPhi between B jets
496  double dPhiBB=0;
497  double leadBJetID=0;
498  double subBJetID=1;
499  if(passBJetCut){
500  dPhiBB=std::acos(std::cos(signalBJets.at(leadBJetID)->phi()-signalBJets.at(subBJetID)->phi()));
501  }
502 
503  bool passJetCutSRC=false;
504 
505  /*std::cout << "JET PT CHECK ";
506  for(const HEPUtils::Jet* jet : signalJets){
507  std::cout << jet->pT() << " ";
508  }
509  std::cout << endl;*/
510 
511  if(nJets==5){
512  if(signalJets[0]->pT() > 80.
513  && signalJets[1]->pT() > 80.
514  && signalJets[2]->pT() > 35.
515  && signalJets[3]->pT() > 35.
516  && signalJets[4]->pT() > 35.)passJetCutSRC=true;
517  }
518 
519  bool isSRC1=false;
520  bool isSRC2=false;
521  bool isSRC3=false;
522 
523  if(cut_LeptonVeto && cut_Btag && cut_METGt150 && cut_dPhiJets && cut_mTbjetmetGt175 && cut_tau){
524 
525  if(passJetCutSRC && mT_bjetmet_min > 185. && mT_bjetmet_max>205. && met>160. && dPhiBB > (0.2*3.14))isSRC1=true;
526 
527  if(passJetCutSRC && mT_bjetmet_min > 200. && mT_bjetmet_max>290. && met>160. && dPhiBB > (0.2*3.14))isSRC2=true;
528 
529  if(passJetCutSRC && mT_bjetmet_min > 200. && mT_bjetmet_max>325. && met>215. && dPhiBB > (0.2*3.14))isSRC3=true;
530 
531  }
532 
533  cutFlowVector_str[0] = "No cuts ";
534  cutFlowVector_str[1] = "MET > 130 GeV ";
535  cutFlowVector_str[2] = "Lepton veto ";
536  cutFlowVector_str[3] = "MET > 150 GeV ";
537  cutFlowVector_str[4] = "Jet multiplicity and pT ";
538  cutFlowVector_str[5] = "dPhi(jet,MET) > pi/5 ";
539  cutFlowVector_str[6] = ">=2 b jets ";
540  cutFlowVector_str[7] = "tau veto ";
541  cutFlowVector_str[8] = "mT(b,MET) > 175 ";
542  cutFlowVector_str[9] = "SRA1 ";
543  cutFlowVector_str[10] = "SRA2 ";
544  cutFlowVector_str[11] = "SRA3 ";
545  cutFlowVector_str[12] = "SRA4 ";
546  cutFlowVector_str[13] = "SRC: exactly 5 jets ";
547  cutFlowVector_str[14] = "SRC: dPhi(jet,MET) ";
548  cutFlowVector_str[15] = "SRC: >=2 b jets ";
549  cutFlowVector_str[16] = "SRC: tau veto ";
550  cutFlowVector_str[17] = "SRC: dPhi(b,b) ";
551  cutFlowVector_str[18] = "SRC1";
552  cutFlowVector_str[19] = "SRC2";
553  cutFlowVector_str[20] = "SRC3";
554 
555  for(int j=0;j<NCUTS;j++){
556  if(
557  (j==0) ||
558 
559  (j==1 && cut_METGt130) ||
560 
561  (j==2 && cut_METGt130 && cut_LeptonVeto) ||
562 
563  (j==3 && cut_METGt150 && cut_LeptonVeto) ||
564 
565  (j==4 && cut_METGt150 && cut_LeptonVeto && cut_6jets) ||
566 
567  (j==5 && cut_METGt150 && cut_LeptonVeto && cut_6jets && cut_dPhiJets) ||
568 
569  (j==6 && cut_METGt150 && cut_LeptonVeto && cut_6jets && cut_dPhiJets && cut_Btag) ||
570 
571  (j==7 && cut_METGt150 && cut_LeptonVeto && cut_6jets && cut_dPhiJets && cut_Btag && cut_tau) ||
572 
573  (j==8 && cut_METGt150 && cut_LeptonVeto && cut_6jets && cut_dPhiJets && cut_Btag && cut_tau && cut_mTbjetmetGt175) ||
574 
575  (j==9 && cut_METGt150 && cut_LeptonVeto && cut_6jets && cut_dPhiJets && cut_Btag && cut_tau && cut_mTbjetmetGt175 && isSRA1) ||
576 
577  (j==10 && cut_METGt150 && cut_LeptonVeto && cut_6jets && cut_dPhiJets && cut_Btag && cut_tau && cut_mTbjetmetGt175 && isSRA2) ||
578 
579  (j==11 && cut_METGt150 && cut_LeptonVeto && cut_6jets && cut_dPhiJets && cut_Btag && cut_tau && cut_mTbjetmetGt175 && isSRA3) ||
580 
581  (j==12 && cut_METGt150 && cut_LeptonVeto && cut_6jets && cut_dPhiJets && cut_Btag && cut_tau && cut_mTbjetmetGt175 && isSRA4) ||
582 
583  (j==13 && cut_METGt150 && cut_LeptonVeto && passJetCutSRC) ||
584 
585  (j==14 && cut_METGt150 && cut_LeptonVeto && passJetCutSRC && cut_dPhiJets) ||
586 
587  (j==15 && cut_METGt150 && cut_LeptonVeto && passJetCutSRC && cut_dPhiJets && cut_Btag) ||
588 
589  (j==16 && cut_METGt150 && cut_LeptonVeto && passJetCutSRC && cut_dPhiJets && cut_Btag && cut_tau) ||
590 
591  (j==17 && cut_METGt150 && cut_LeptonVeto && passJetCutSRC && cut_dPhiJets && cut_Btag && cut_tau && dPhiBB > (0.2*3.14)) ||
592 
593  (j==18 && isSRC1) ||
594 
595  (j==19 && isSRC2) ||
596 
597  (j==20 && isSRC3)
598 
599  ){
600 
601  cutFlowVector[j]++;
602 
603  }
604  }
605 
606 
607  /*for(int j=0;j<NCUTS;j++){
608  if(
609  (j==0) ||
610 
611  (j==1 && cut_MuonVeto) ||
612 
613  (j==2 && cut_ElectronVeto && cut_MuonVeto) ||
614 
615  (j==3 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130) ||
616 
617  (j==4 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets) ||
618 
619  (j==5 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets && cut_dPhiJets) ||
620 
621  (j==6 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets && cut_dPhiJets && cut_tau) ||
622 
623  (j==7 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets && cut_dPhiJets && cut_tau && cut_Btag) ||
624 
625  (j==8 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets && cut_dPhiJets && cut_tau && cut_Btag && cut_mTbjetmetGt175) ||
626 
627  (j==9 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets && cut_dPhiJets && cut_tau && cut_Btag && cut_mTbjetmetGt175 && cut_mjjj0) ||
628 
629  (j==10 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets && cut_dPhiJets && cut_tau && cut_Btag && cut_mTbjetmetGt175 && cut_mjjj0 && cut_mjjj1) ||
630 
631  (j==11 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets && cut_dPhiJets && cut_tau && cut_Btag && cut_mTbjetmetGt175 && cut_mjjj0 && cut_mjjj1 && cut_METGt150) ||
632 
633  (j==12 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets && cut_dPhiJets && cut_tau && cut_Btag && cut_mTbjetmetGt175 && cut_mjjj0 && cut_mjjj1 && cut_METGt200) ||
634 
635  (j==13 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets && cut_dPhiJets && cut_tau && cut_Btag && cut_mTbjetmetGt175 && cut_mjjj0 && cut_mjjj1 && cut_METGt250) ||
636 
637  (j==14 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets && cut_dPhiJets && cut_tau && cut_Btag && cut_mTbjetmetGt175 && cut_mjjj0 && cut_mjjj1 && cut_METGt300) ||
638 
639  (j==15 && cut_ElectronVeto && cut_MuonVeto && cut_METGt130 && cut_6jets && cut_dPhiJets && cut_tau && cut_Btag && cut_mTbjetmetGt175 && cut_mjjj0 && cut_mjjj1 && cut_METGt350) )
640 
641 
642  cutFlowVector[j]++;
643  }*/
644 
645  //We're now ready to apply the cuts for each signal region
646  //_numSR1, _numSR2, _numSR3;
647 
648  if(isSRA1) _numSRA1 += event->weight();
649  if(isSRA2) _numSRA2 += event->weight();
650  if(isSRA3) _numSRA3 += event->weight();
651  if(isSRA4) _numSRA4 += event->weight();
652 
653  if(isSRC1) _numSRC1 += event->weight();
654  if(isSRC2) _numSRC2 += event->weight();
655  if(isSRC3) _numSRC3 += event->weight();
656 
657  return;
658 
659  }
660 
662  void combine(const Analysis* other)
663  {
664  const Analysis_ATLAS_8TeV_0LEPStop_20invfb* specificOther
665  = dynamic_cast<const Analysis_ATLAS_8TeV_0LEPStop_20invfb*>(other);
666 
667  if (NCUTS != specificOther->NCUTS) NCUTS = specificOther->NCUTS;
668  for (int j=0; j<NCUTS; j++) {
669  cutFlowVector[j] += specificOther->cutFlowVector[j];
670  cutFlowVector_str[j] = specificOther->cutFlowVector_str[j];
671  }
672  _numSRA1 += specificOther->_numSRA1;
673  _numSRA2 += specificOther->_numSRA2;
674  _numSRA3 += specificOther->_numSRA3;
675  _numSRA4 += specificOther->_numSRA4;
676  _numSRC1 += specificOther->_numSRC1;
677  _numSRC2 += specificOther->_numSRC2;
678  _numSRC3 += specificOther->_numSRC3;
679  }
680 
681 
682  void collect_results() {
683 
684  // add_result(SignalRegionData("SR label", n_obs, {n_sig_MC, n_sig_MC_sys}, {n_bkg, n_bkg_err}));
685 
686  add_result(SignalRegionData("SRA1", 11., {_numSRA1, 0.}, {15.8, 1.9}));
687  add_result(SignalRegionData("SRA2", 4., {_numSRA2, 0.}, {4.1, 0.8}));
688  add_result(SignalRegionData("SRA3", 5., {_numSRA3, 0.}, {4.1, 0.9}));
689  add_result(SignalRegionData("SRA4", 4., {_numSRA4, 0.}, {2.4, 0.7}));
690  add_result(SignalRegionData("SRC1", 59., {_numSRC1, 0.}, {68., 7.}));
691  add_result(SignalRegionData("SRC2", 30., {_numSRC2, 0.}, {34., 5.}));
692  add_result(SignalRegionData("SRC3", 15., {_numSRC3, 0.}, {20.3, 3.}));
693 
694  return;
695  }
696 
697 
698  protected:
699  void analysis_specific_reset() {
700  _numSRA1 = 0 ; _numSRA2 = 0; _numSRA3 = 0; _numSRA4 = 0;
701  _numSRC1 = 0 ; _numSRC2 = 0; _numSRC3 = 0;
702 
703  std::fill(cutFlowVector.begin(), cutFlowVector.end(), 0);
704  }
705 
706  };
707 
708  DEFINE_ANALYSIS_FACTORY(ATLAS_8TeV_0LEPStop_20invfb)
709 
710 
711  }
712 }
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 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
bool sortByPT(const HEPUtils::Jet *jet1, const HEPUtils::Jet *jet2)
A class for collider analyses within ColliderBit.
Definition: Analysis.hpp:41
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.
def combine(region_file, pip_file)
Definition: colouring.py:169
void applyElectronEff(std::vector< const HEPUtils::Particle *> &electrons)
Randomly filter the supplied particle list by parameterised electron efficiency.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
void applyTauEfficiencyR1(std::vector< const HEPUtils::Particle *> &taus)
Randomly filter the supplied particle list by parameterised Run 1 tau efficiency. ...
Class for ColliderBit analyses.
double E(const double x, const double y)
Functions that do super fast ATLAS detector simulation based on four-vector smearing.