gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
Analysis_ATLAS_13TeV_0LEPStop_36invfb.cpp
Go to the documentation of this file.
1 #include <vector>
2 #include <cmath>
3 #include <memory>
4 #include <iomanip>
5 
9 
10 // #include "RestFrames/RestFrames.hh"
11 
12 using namespace std;
13 
14 /* The ATLAS 0 lepton direct stop analysis
15 
16  Based on: https://arxiv.org/abs/1709.04183
17 
18  Code by Martin White
19 
20  KNOWN ISSUES
21 
22  a) Should apply Very Loose selection to electron candidates
23  b) Cannot apply requirement that the ETmiss calculated from tracking information is aligned in phi with that calculated from the calo system.
24  c) We do not apply the tau veto. Could approximate by removing events with tagged taus in?
25 
26 Note: have removed RJ code for now to save time (the cutflows are too divergent to make this code usable, probably due to the Pythia ISR modelling).
27 
28 */
29 
30 namespace Gambit {
31  namespace ColliderBit {
32 
33  // Need two different functions here for use with std::sort
34  bool sortByPT13(const HEPUtils::Jet* jet1, const HEPUtils::Jet* jet2) { return (jet1->pT() > jet2->pT()); }
35  bool sortByPT13_sharedptr(std::shared_ptr<HEPUtils::Jet> jet1, std::shared_ptr<HEPUtils::Jet> jet2) { return sortByPT13(jet1.get(), jet2.get()); }
36 
37  // Need two different functions here for use with std::sort
38  bool sortByMass(const HEPUtils::Jet* jet1, const HEPUtils::Jet* jet2) { return (jet1->mass() > jet2->mass()); }
39  bool sortByMass_sharedptr(std::shared_ptr<HEPUtils::Jet> jet1, std::shared_ptr<HEPUtils::Jet> jet2) { return sortByMass(jet1.get(), jet2.get()); }
40 
41  double calcMT(HEPUtils::P4 jetMom,HEPUtils::P4 metMom){
42 
43  //std::cout << "metMom.px() " << metMom.px() << " jetMom PT " << jetMom.pT() << std::endl;
44 
45  double met=sqrt(metMom.px()*metMom.px()+metMom.py()*metMom.py());
46  double dphi = metMom.deltaPhi(jetMom);
47  double mt=sqrt(2*jetMom.pT()*met*(1-std::cos(dphi)));
48 
49  return mt;
50 
51  }
52 
53 
54  class Analysis_ATLAS_13TeV_0LEPStop_36invfb : public Analysis {
55  private:
56 
57  // Numbers passing cuts
58 
59  std::map<string, EventCounter> _counters = {
60  {"SRA_TT", EventCounter("SRA_TT")},
61  {"SRA_TW", EventCounter("SRA_TW")},
62  {"SRA_T0", EventCounter("SRA_T0")},
63  {"SRB_TT", EventCounter("SRB_TT")},
64  {"SRB_TW", EventCounter("SRB_TW")},
65  {"SRB_T0", EventCounter("SRB_T0")},
66  {"SRC1", EventCounter("SRC1")},
67  {"SRC2", EventCounter("SRC2")},
68  {"SRC3", EventCounter("SRC3")},
69  {"SRC4", EventCounter("SRC4")},
70  {"SRC5", EventCounter("SRC5")},
71  {"SRD_low", EventCounter("SRD_low")},
72  {"SRD_high", EventCounter("SRD_high")},
73  {"SRE", EventCounter("SRE")},
74  };
75 
76  vector<int> cutFlowVector;
77  vector<string> cutFlowVector_str;
78  int NCUTS; //=16;
79 
80  // Define RestFrames objects
81 
82  /*unique_ptr<RestFrames::LabRecoFrame> LAB;
83  unique_ptr<RestFrames::DecayRecoFrame> CM;
84  unique_ptr<RestFrames::DecayRecoFrame> S;
85  unique_ptr<RestFrames::VisibleRecoFrame> ISR;
86  unique_ptr<RestFrames::VisibleRecoFrame> V;
87  unique_ptr<RestFrames::InvisibleRecoFrame> I;
88  unique_ptr<RestFrames::InvisibleGroup> INV;
89  unique_ptr<RestFrames::CombinatoricGroup> VIS;
90  unique_ptr<RestFrames::SetMassInvJigsaw> InvMass;
91  unique_ptr<RestFrames::MinMassesCombJigsaw> SplitVis;*/
92 
93 
94 
95  void LeptonLeptonOverlapRemoval(vector<const HEPUtils::Particle*> &lep1vec, vector<const HEPUtils::Particle*> &lep2vec, double DeltaRMax) {
96 
97  //Routine to do jet-lepton check
98  //Discards jets if they are within DeltaRMax of a lepton
99 
100  vector<const HEPUtils::Particle*> Survivors;
101 
102  for(unsigned int itlep1 = 0; itlep1 < lep1vec.size(); itlep1++) {
103  bool overlap = false;
104  HEPUtils::P4 lep1mom=lep1vec.at(itlep1)->mom();
105  for(unsigned int itlep2 = 0; itlep2 < lep2vec.size(); itlep2++) {
106  HEPUtils::P4 lep2mom=lep2vec.at(itlep2)->mom();
107  double dR;
108 
109  dR=lep1mom.deltaR_eta(lep2mom);
110 
111  if(fabs(dR) <= DeltaRMax) overlap=true;
112  }
113  if(overlap) continue;
114  Survivors.push_back(lep1vec.at(itlep1));
115  }
116  lep1vec=Survivors;
117 
118  return;
119  }
120 
121  void JetLeptonOverlapRemoval(vector<const HEPUtils::Jet*> &jetvec, vector<const HEPUtils::Particle*> &lepvec, double DeltaRMax) {
122  //Routine to do jet-lepton check
123  //Discards jets if they are within DeltaRMax of a lepton
124 
125  vector<const HEPUtils::Jet*> Survivors;
126 
127  for(unsigned int itjet = 0; itjet < jetvec.size(); itjet++) {
128  bool overlap = false;
129  HEPUtils::P4 jetmom=jetvec.at(itjet)->mom();
130  for(unsigned int itlep = 0; itlep < lepvec.size(); itlep++) {
131  HEPUtils::P4 lepmom=lepvec.at(itlep)->mom();
132  double dR;
133 
134  dR=jetmom.deltaR_eta(lepmom);
135 
136  if(fabs(dR) <= DeltaRMax) overlap=true;
137  }
138  if(overlap) continue;
139  Survivors.push_back(jetvec.at(itjet));
140  }
141  jetvec=Survivors;
142 
143  return;
144  }
145 
146  void LeptonJetOverlapRemoval(vector<const HEPUtils::Particle*> &lepvec, vector<const HEPUtils::Jet*> &jetvec, double DeltaRMax) {
147  //Routine to do lepton-jet check
148  //Discards leptons if they are within DeltaRMax of a jet
149 
150  vector<const HEPUtils::Particle*> Survivors;
151 
152  for(unsigned int itlep = 0; itlep < lepvec.size(); itlep++) {
153  bool overlap = false;
154  HEPUtils::P4 lepmom=lepvec.at(itlep)->mom();
155  for(unsigned int itjet= 0; itjet < jetvec.size(); itjet++) {
156  HEPUtils::P4 jetmom=jetvec.at(itjet)->mom();
157  double dR;
158 
159  dR=jetmom.deltaR_eta(lepmom);
160 
161  if(fabs(dR) <= DeltaRMax) overlap=true;
162  }
163  if(overlap) continue;
164  Survivors.push_back(lepvec.at(itlep));
165  }
166  lepvec=Survivors;
167 
168  return;
169  }
170 
171 
172  public:
173 
174  // Required detector sim
175  static constexpr const char* detector = "ATLAS";
176 
177  Analysis_ATLAS_13TeV_0LEPStop_36invfb() {
178 
179  set_analysis_name("ATLAS_13TeV_0LEPStop_36invfb");
180  set_luminosity(36.);
181 
182  NCUTS=120;
183 
184  for(int i=0;i<NCUTS;i++){
185  cutFlowVector.push_back(0);
186  cutFlowVector_str.push_back("");
187  }
188 
189  // RestFrames initialisation
190 
191  /*LAB.reset(new RestFrames::LabRecoFrame("LAB","lab"));
192  CM.reset(new RestFrames::DecayRecoFrame("CM","cm"));
193  S.reset(new RestFrames::DecayRecoFrame("S","s"));
194  ISR.reset(new RestFrames::VisibleRecoFrame("ISR","isr"));
195  V.reset(new RestFrames::VisibleRecoFrame("V","v"));
196  I.reset(new RestFrames::InvisibleRecoFrame("I","i"));
197 
198  // Connect the frames
199  LAB->SetChildFrame(*CM);
200  CM->AddChildFrame(*ISR);
201  CM->AddChildFrame(*S);
202  S->AddChildFrame(*V);
203  S->AddChildFrame(*I);
204 
205  // Initialize the tree
206  LAB->InitializeTree();
207 
208  // Define groups
209  INV.reset(new RestFrames::InvisibleGroup("INV","inv"));
210  INV->AddFrame(*I);
211  VIS.reset(new RestFrames::CombinatoricGroup("VIS","vis"));
212  VIS->AddFrame(*ISR);
213  VIS->SetNElementsForFrame(*ISR,1,false);
214  VIS->AddFrame(*V);
215  VIS->SetNElementsForFrame(*V,0,false);
216 
217  // set the invisible system mass to zero
218  InvMass.reset(new RestFrames::SetMassInvJigsaw("InvMass","kSetMass"));
219  INV->AddJigsaw(*InvMass);
220 
221  // define the rule for partitioning objects between "ISR" and "V"
222  SplitVis.reset(new RestFrames::MinMassesCombJigsaw("CombPPJigsaw", "kMinMasses"));
223  VIS->AddJigsaw(*SplitVis);
224  // "0" group (ISR)
225  SplitVis->AddFrame(*ISR, 0);
226  // "1" group (V + I)
227  SplitVis->AddFrame(*V,1);
228  SplitVis->AddFrame(*I,1);
229 
230  LAB->InitializeAnalysis();*/
231 
232  }
233 
234 
235 
236  void run(const HEPUtils::Event* event) {
237 
238  // Missing energy
239  HEPUtils::P4 metVec = event->missingmom();
240  double Met = event->met();
241 
242 
243  // Baseline lepton objects
244  vector<const HEPUtils::Particle*> baselineElectrons, baselineMuons, baselineTaus;
245 
246  for (const HEPUtils::Particle* electron : event->electrons()) {
247  if (electron->pT() > 7. && electron->abseta() < 2.47) baselineElectrons.push_back(electron);
248  }
249  for (const HEPUtils::Particle* muon : event->muons()) {
250  if (muon->pT() > 6. && muon->abseta() < 2.7) baselineMuons.push_back(muon);
251  }
252 
253  // Apply lepton efficiencies
254  ATLAS::applyElectronEff(baselineElectrons);
255  ATLAS::applyMuonEff(baselineMuons);
256 
257  // Photons
258  vector<const HEPUtils::Particle*> signalPhotons;
259  for (const HEPUtils::Particle* photon : event->photons()) {
260  signalPhotons.push_back(photon);
261  }
262 
263 
264  // No taus used in 13 TeV analysis?
265  //for (const HEPUtils::Particle* tau : event->taus()) {
266  //if (tau->pT() > 10. && tau->abseta() < 2.47) baselineTaus.push_back(tau);
267  //}
268  //ATLAS::applyTauEfficiencyR1(baselineTaus);
269 
270 
271  // Jets
272  vector<const HEPUtils::Jet*> bJets;
273  vector<const HEPUtils::Jet*> nonBJets;
274  vector<const HEPUtils::Jet*> trueBJets; //for debugging
275 
276  // Get b jets
278  const std::vector<double> a = {0,10.};
279  const std::vector<double> b = {0,10000.};
280  const std::vector<double> c = {0.77}; // set b-tag efficiency to 77%
281  HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
282  for (const HEPUtils::Jet* jet : event->jets())
283  {
284  bool hasTag=has_tag(_eff2d, fabs(jet->eta()), jet->pT());
285  if (jet->pT() > 20. && fabs(jet->eta()) < 2.8)
286  {
287  if(jet->btag() && hasTag && fabs(jet->eta()) < 2.5 && jet->pT() > 20.)
288  {
289  bJets.push_back(jet);
290  }
291  else
292  {
293  nonBJets.push_back(jet);
294  }
295  }
296  }
297 
298  // Overlap removal
299  // Note: use paper description instead of code snippet
300  JetLeptonOverlapRemoval(nonBJets,baselineElectrons,0.2);
301  LeptonJetOverlapRemoval(baselineElectrons,nonBJets,0.4);
302  LeptonJetOverlapRemoval(baselineElectrons,bJets,0.4);
303  LeptonJetOverlapRemoval(baselineMuons,nonBJets,0.4);
304  LeptonJetOverlapRemoval(baselineMuons,bJets,0.4);
305 
306  // Fill a jet-pointer-to-bool map to make it easy to check
307  // if a given jet is treated as a b-jet in this analysis
308  map<const HEPUtils::Jet*,bool> analysisBtags;
309  for (const HEPUtils::Jet* jet : bJets) analysisBtags[jet] = true;
310  for (const HEPUtils::Jet* jet : nonBJets) analysisBtags[jet] = false;
311 
312  // Signal object containers
313  vector<const HEPUtils::Particle*> signalElectrons;
314  vector<const HEPUtils::Particle*> signalMuons;
315  vector<const HEPUtils::Particle*> electronsForVeto;
316  vector<const HEPUtils::Particle*> muonsForVeto;
317 
318  vector<const HEPUtils::Jet*> signalJets;
319  vector<const HEPUtils::Jet*> signalBJets;
320  vector<const HEPUtils::Jet*> signalNonBJets;
321 
322  // It seems that there are no extra signal jet requirements (unlike 8 TeV analysis)
323  for (const HEPUtils::Jet* jet : bJets) {
324  signalJets.push_back(jet);
325  signalBJets.push_back(jet);
326  }
327 
328  for (const HEPUtils::Jet* jet : nonBJets) {
329  signalJets.push_back(jet);
330  signalNonBJets.push_back(jet);
331  }
332 
333  //Put signal jets in pT order
334  std::sort(signalJets.begin(), signalJets.end(), sortByPT13);
335  std::sort(signalBJets.begin(), signalBJets.end(), sortByPT13);
336  std::sort(signalNonBJets.begin(), signalNonBJets.end(), sortByPT13);
337 
338  for (const HEPUtils::Particle* electron : baselineElectrons) {
339  signalElectrons.push_back(electron);
340  }
341 
342  for (const HEPUtils::Particle* muon : baselineMuons) {
343  signalMuons.push_back(muon);
344  }
345 
346  // Need to recluster jets at this point (R=0.8 and R=1.2)
347  vector<std::shared_ptr<HEPUtils::Jet>> fatJetsR8=get_jets(signalJets,0.8);
348  vector<std::shared_ptr<HEPUtils::Jet>> fatJetsR12=get_jets(signalJets,1.2);
349 
350  //Put 1_2 signal jets in decreasing pT order
351  std::sort(fatJetsR12.begin(), fatJetsR12.end(), sortByPT13_sharedptr);
352 
353  //Put 0_8 signal jets in pT order
354  std::sort(fatJetsR8.begin(), fatJetsR8.end(), sortByPT13_sharedptr);
355 
356  // We now have the signal electrons, muons, jets and b jets- move on to the analysis
357 
358  // The following code follow the ATLAS public snippet closely
359  float DRBB = 0;
360  int NBJets = signalBJets.size();
361 
362  float AntiKt8M_0 = 0;
363  float AntiKt8M_1 = 0;
364  float AntiKt12M_0 = 0;
365  float AntiKt12M_1 = 0;
366  // float MtTauCand = -1 ;
367  float MtBMin = 0 ;
368  float MtBMax = 0 ;
369 
370  if (fatJetsR8.size()>0) AntiKt8M_0 = fatJetsR8[0]->mass() ;
371  if (fatJetsR8.size()>1) AntiKt8M_1 = fatJetsR8[1]->mass() ;
372  if (fatJetsR12.size()>0) AntiKt12M_0 = fatJetsR12[0]->mass() ;
373  if (fatJetsR12.size()>1) AntiKt12M_1 = fatJetsR12[1]->mass() ;
374  if (NBJets>1) DRBB = signalBJets[1]->mom().deltaR_eta(signalBJets[0]->mom());
375 
376  double dPhi_min = 1000.;
377  double dPhi_max = 0.;
378  if (signalBJets.size()>=2) {
379  for (const HEPUtils::Jet* jet : signalBJets) {
380  double dphi = fabs(metVec.deltaPhi(jet->mom()));
381  if (dphi<dPhi_min) {
382  dPhi_min=dphi;
383  MtBMin=calcMT(jet->mom(),metVec);
384  }
385  if (dphi>dPhi_max) {
386  dPhi_max=dphi;
387  MtBMax=calcMT(jet->mom(),metVec);
388  }
389  }
390  }
391 
392  float realWMass = 80.385;
393  float realTopMass = 173.210;
394 
395  //Chi2 method
396  double Chi2min = 99999999999999999.;
397  int W1j1_low = -1,W1j2_low = -1,W2j1_low = -1,W2j2_low = -1,b1_low = -1,b2_low = -1;
398 
399  double m_mt2Chi2 = 0;
400 
401  if (signalJets.size()>=4 && signalBJets.size()>=2 && signalNonBJets.size()>=2)
402  {
403  for(int W1j1=0; W1j1<(int)signalNonBJets.size(); W1j1++) {// <------------------This lines has to be replaced
404  for(int W2j1=0;W2j1<(int)signalNonBJets.size(); W2j1++) {// <------------------This lines has to be replaced
405  if (W2j1==W1j1) continue;// <------------------This lines has to be added
406  // for(int W1j1=0; W1j1<(int)ljets.size()-1; W1j1++) {
407  // for(int W2j1=W1j1+1;W2j1<(int)ljets.size(); W2j1++) {
408  for(int b1=0;b1<(int)signalBJets.size();b1++){
409  for(int b2=0;b2<(int)signalBJets.size();b2++){
410  if(b2==b1) continue;
411  double chi21, chi22, mW1, mW2, mt1, mt2;
412 
413  if(W2j1>W1j1){
414 
415  mW1 = signalNonBJets[W1j1]->mass();
416  mW2 = signalNonBJets[W2j1]->mass();
417  mt1 = (signalNonBJets[W1j1]->mom()+signalBJets[b1]->mom()).m();
418  mt2 = (signalNonBJets[W2j1]->mom()+signalBJets[b2]->mom()).m();
419 
420  chi21 = (mW1-realWMass)*(mW1-realWMass)/realWMass + (mt1-realTopMass)*(mt1-realTopMass)/realTopMass;
421  chi22 = (mW2-realWMass)*(mW2-realWMass)/realWMass + (mt2-realTopMass)*(mt2-realTopMass)/realTopMass;
422 
423  if(Chi2min > (chi21 + chi22)){
424  Chi2min = chi21 + chi22;
425  if(chi21 < chi22){
426  W1j1_low = W1j1;
427  W1j2_low = -1;
428  W2j1_low = W2j1;
429  W2j2_low = -1;
430  b1_low = b1;
431  b2_low = b2;
432  }
433  else{
434  W2j1_low = W1j1;
435  W2j2_low = -1;
436  W1j1_low = W2j1;
437  W1j2_low = -1;
438  b2_low = b1;
439  b1_low = b2;
440  }
441  }
442  }
443 
444  if (signalNonBJets.size()<3)
445  continue;
446 
447  for(int W1j2=W1j1+1;W1j2 < (int)signalNonBJets.size(); W1j2++) {
448  if(W1j2==W2j1) continue;
449 
450  //try bll,bl top candidates
451  mW1 = (signalNonBJets[W1j1]->mom() + signalNonBJets[W1j2]->mom()).m();
452  mW2 = (signalNonBJets[W2j1])->mass();
453  mt1 = (signalNonBJets[W1j1]->mom() + signalNonBJets[W1j2]->mom() + signalBJets[b1]->mom()).m();
454  mt2 = (signalNonBJets[W2j1]->mom() + signalBJets[b2]->mom()).m();
455  chi21 = (mW1-realWMass)*(mW1-realWMass)/realWMass + (mt1-realTopMass)*(mt1-realTopMass)/realTopMass;
456  chi22 = (mW2-realWMass)*(mW2-realWMass)/realWMass + (mt2-realTopMass)*(mt2-realTopMass)/realTopMass;
457  if(Chi2min > (chi21 + chi22)){
458  Chi2min = chi21 + chi22;
459  if(chi21 < chi22){
460  W1j1_low = W1j1;
461  W1j2_low = W1j2;
462  W2j1_low = W2j1;
463  W2j2_low = -1;
464  b1_low = b1;
465  b2_low = b2;
466  }
467  else{
468  W2j1_low = W1j1;
469  W2j2_low = W1j2;
470  W1j1_low = W2j1;
471  W1j2_low = -1;
472  b2_low = b1;
473  b1_low = b2;
474  }
475  }
476  if(signalNonBJets.size() < 4)continue;
477  //try bll, bll top candidates
478  for(int W2j2=W2j1+1;W2j2<(int)signalNonBJets.size(); W2j2++){
479  if((W2j2==W1j1) || (W2j2==W1j2)) continue;
480  if(W2j1<W1j1) continue; //runtime reasons, we don't want combinations checked twice <--------------------This line should be added
481  mW1 = (signalNonBJets[W1j1]->mom() + signalNonBJets[W1j2]->mom()).m();
482  mW2 = (signalNonBJets[W2j1]->mom() + signalNonBJets[W2j2]->mom()).m();
483  mt1 = (signalNonBJets[W1j1]->mom() + signalNonBJets[W1j2]->mom() + signalBJets[b1]->mom()).m();
484  mt2 = (signalNonBJets[W2j1]->mom() + signalNonBJets[W2j2]->mom() + signalBJets[b2]->mom()).m();
485  chi21 = (mW1-realWMass)*(mW1-realWMass)/realWMass + (mt1-realTopMass)*(mt1-realTopMass)/realTopMass;
486  chi22 = (mW2-realWMass)*(mW2-realWMass)/realWMass + (mt2-realTopMass)*(mt2-realTopMass)/realTopMass;
487  if(Chi2min > (chi21 + chi22)){
488  Chi2min = chi21 + chi22;
489  if(chi21 < chi22){
490  W1j1_low = W1j1;
491  W1j2_low = W1j2;
492  W2j1_low = W2j1;
493  W2j2_low = W2j2;
494  b1_low = b1;
495  b2_low = b2;
496  }
497  else{
498  W2j1_low = W1j1;
499  W2j2_low = W1j2;
500  W1j1_low = W2j1;
501  W1j2_low = W2j2;
502  b2_low = b1;
503  b1_low = b2;
504  }
505  }
506  }
507  }
508  }
509  }
510  }
511  }
512 
513  HEPUtils::P4 WCand0=signalNonBJets[W1j1_low]->mom();
514  if (W1j2_low != -1) WCand0 +=signalNonBJets[W1j2_low]->mom();
515  HEPUtils::P4 topCand0 = WCand0 + signalBJets[b1_low]->mom();
516 
517  HEPUtils::P4 WCand1 = signalNonBJets[W2j1_low]->mom();
518  if(W2j2_low != -1) WCand1 += signalNonBJets[W2j2_low]->mom();
519  HEPUtils::P4 topCand1 = WCand1 + signalBJets[b2_low]->mom();
520 
521  HEPUtils::P4 tempTop0=HEPUtils::P4::mkEtaPhiMPt(0.,topCand0.phi(),173.210,topCand0.pT());
522  HEPUtils::P4 tempTop1=HEPUtils::P4::mkEtaPhiMPt(0.,topCand1.phi(),173.210,topCand1.pT());
523 
524  // Note that the first component here is the mass
525  // This must be the top mass (i.e. mass of our vectors) and not zero!
526 
527  double pa_a[3] = { tempTop0.m() , tempTop0.px(), tempTop0.py() };
528  double pb_a[3] = { tempTop1.m() , tempTop1.px(), tempTop1.py() };
529  double pmiss_a[3] = { 0, metVec.px(), metVec.py() };
530  double mn_a = 0.;
531 
532  mt2_bisect::mt2 mt2_event_a;
533 
534  mt2_event_a.set_momenta(pa_a,pb_a,pmiss_a);
535  mt2_event_a.set_mn(mn_a);
536 
537  m_mt2Chi2 = mt2_event_a.get_mt2();
538 
539  }
540 
541  float MT2Chi2 = m_mt2Chi2;
542 
543  // RestFrames stuff
544 
545  /*double CA_PTISR=0;
546  double CA_MS=0;
547  double CA_NbV=0;
548  double CA_NjV=0;
549  double CA_RISR=0;
550  double CA_dphiISRI=0;
551  double CA_pTjV4=0;
552  double CA_pTbV1=0;
553 
554  int m_NjV(0);
555  int m_NbV(0);
556  int m_NbISR(0);
557  double m_pTjV4(0.);
558  double m_pTbV1(0);
559  double m_PTISR(0.);
560  double m_MS(0.);
561  double m_RISR(0.);
562  double m_dphiISRI(0.);
563 
564  LAB->ClearEvent();
565 
566  if (!(Met<250 || (baselineElectrons.size()+baselineMuons.size())>0 || signalJets.size()<4 || signalBJets.size()<1 || signalJets[3]->pT()<40)){
567 
568  std::vector<RestFrames::RFKey> jetID;
569 
570  for(size_t i=0;i<signalJets.size();i++){
571 
572  TLorentzVector tmpJet;
573  tmpJet.SetPtEtaPhiM(signalJets[i]->pT(),0.,signalJets[i]->phi(),signalJets[i]->mass());
574 
575  jetID.push_back(VIS->AddLabFrameFourVector(tmpJet));
576 
577  }
578 
579  TVector3 ETMiss;
580  ETMiss.SetXYZ(metVec.px(),metVec.py(),0.0);
581  INV->SetLabFrameThreeVector(ETMiss);
582 
583  if(!LAB->AnalyzeEvent()) std::cout << "Something went wrong..." << std::endl;
584 
585  for(size_t i = 0; i < signalJets.size(); i++){
586  if (VIS->GetFrame(jetID[i]) == *V){ // sparticle group
587  m_NjV++;
588  if (m_NjV == 4) m_pTjV4 = signalJets[i]->pT();
589  if ( analysisBtags.at(signalJets[i]) && fabs(signalJets[i]->eta())<2.5) {
590  m_NbV++;
591  if (m_NbV == 1) m_pTbV1 = signalJets[i]->pT();
592  }
593  } else {
594  if ( analysisBtags.at(signalJets[i]) && fabs(signalJets[i]->eta())<2.5)
595  m_NbISR++;
596  }
597  }
598 
599  // need at least one jet associated with MET-side of event
600  if(m_NjV >= 1)
601  {
602  TVector3 vP_ISR = ISR->GetFourVector(*CM).Vect();
603  TVector3 vP_I = I->GetFourVector(*CM).Vect();
604 
605  m_PTISR = vP_ISR.Mag();
606  m_RISR = fabs(vP_I.Dot(vP_ISR.Unit())) / m_PTISR;
607 
608  m_MS = S->GetMass();
609 
610  m_dphiISRI = fabs(vP_ISR.DeltaPhi(vP_I));
611 
612  CA_PTISR=m_PTISR;
613  CA_MS=m_MS;
614  CA_NbV=m_NbV;
615  CA_NjV=m_NjV;
616  CA_RISR=m_RISR;
617  CA_dphiISRI=m_dphiISRI;
618  CA_pTjV4=m_pTjV4;
619  CA_pTbV1=m_pTbV1;
620  }
621  }*/
622 
623 
624  bool isSRA_TT=false;
625  bool isSRA_TW=false;
626  bool isSRA_T0=false;
627  bool isSRB_TT=false;
628  bool isSRB_TW=false;
629  bool isSRB_T0=false;
630  bool isSRC1=false;
631  bool isSRC2=false;
632  bool isSRC3=false;
633  bool isSRC4=false;
634  bool isSRC5=false;
635  bool isSRD_low=false;
636  bool isSRD_high=false;
637  bool isSRE=false;
638 
639  cutFlowVector_str[0] = "No cuts ";
640  cutFlowVector_str[1] = "Derivation skim";
641  cutFlowVector_str[2] = "Lepton veto ";
642  cutFlowVector_str[3] = "Njets >= 4 ";
643  cutFlowVector_str[4] = "Nbjets >= 1 ";
644  cutFlowVector_str[5] = "met > 250 GeV ";
645  cutFlowVector_str[6] = "dPhi(jet,MET) > 0.4 ";
646  cutFlowVector_str[7] = "pT jet 1 > 80 GeV ";
647  cutFlowVector_str[8] = "pT jet 3 > 40 GeV ";
648  cutFlowVector_str[9] = "m jet0, R=1.2 > 120 GeV ";
649  cutFlowVector_str[10] = "SRA-TT: m jet1, R=1.2 > 120 GeV";
650  cutFlowVector_str[11] = "SRA-TT: met > 400 GeV";
651  cutFlowVector_str[12] = "SRA-TT: m jet0, R=0.8 > 60 GeV ";
652  cutFlowVector_str[13] = "SRA-TT: mT(b,MET) min > 200 ";
653  cutFlowVector_str[14] = "SRA-TT: deltaR(b,b) > 1 ";
654  cutFlowVector_str[15] = "SRA-TT: mT2 > 400 GeV";
655  cutFlowVector_str[16] = "SRA-TT: Nbjets >=2 ";
656  cutFlowVector_str[17] = "SRA-TW: m jet1, R=1.2 < 120 GeV";
657  cutFlowVector_str[18] = "SRA-TW: m jet1, R=1.2 > 60 GeV";
658  cutFlowVector_str[19] = "SRA-TW: met > 500 GeV ";
659  cutFlowVector_str[20] = "SRA-TW: m jet0, R=0.8 > 60 GeV";
660  cutFlowVector_str[21] = "SRA-TW: mT(b,MET) min > 200 GeV";
661  cutFlowVector_str[22] = "SRA-TW: mT2 > 400 GeV ";
662  cutFlowVector_str[23] = "SRA-TW: Nbjets >=2 ";
663  cutFlowVector_str[24] = "SRA-T0: m jet1, R=1.2 < 60 GeV";
664  cutFlowVector_str[25] = "SRA-T0: m jet0, R=0.8 > 60 GeV";
665  cutFlowVector_str[26] = "SRA-T0: met > 550 GeV ";
666  cutFlowVector_str[27] = "SRA-T0: mT(b,MET) min > 200 GeV";
667  cutFlowVector_str[28] = "SRA-T0: mT2 > 500 GeV ";
668  cutFlowVector_str[29] = "SRA-T0: Nbjets >=2 ";
669  cutFlowVector_str[30] = "SRB-TT: m jet1, R=1.2 > 120 GeV";
670  cutFlowVector_str[31] = "SRB-TT: deltaR(b,b) > 1.2";
671  cutFlowVector_str[32] = "SRB-TT: mT(b,MET) max > 200 GeV";
672  cutFlowVector_str[33] = "SRB-TT: mT(b,MET) min > 200 GeV";
673  cutFlowVector_str[34] = "SRB-TT: Nbjets >=2 ";
674  cutFlowVector_str[35] = "SRB-TW: m jet1, R=1.2 < 120 GeV";
675  cutFlowVector_str[36] = "SRB-TW: m jet1, R=1.2 > 60 GeV";
676  cutFlowVector_str[37] = "SRB-TW: deltaR(b,b) > 1.2";
677  cutFlowVector_str[38] = "SRB-TW: mT(b,MET) max > 200 GeV";
678  cutFlowVector_str[39] = "SRB-TW: mT(b,MET) min > 200 GeV";
679  cutFlowVector_str[40] = "SRB-TW: Nbjets >=2 ";
680  cutFlowVector_str[41] = "SRB-T0: m jet1, R=1.2 < 60 GeV";
681  cutFlowVector_str[42] = "SRB-T0: mT(b,MET) min > 200 GeV";
682  cutFlowVector_str[43] = "SRB-T0: deltaR(b,b) > 1.2";
683  cutFlowVector_str[44] = "SRB-T0: mT(b,MET) max > 200 GeV";
684  cutFlowVector_str[45] = "SRB-T0: met > 250 GeV ";
685  cutFlowVector_str[46] = "SRB-T0: Nbjets >=2 ";
686 
687  // Cutflow for SRD
688  cutFlowVector_str[47] = "SRD-high: No cuts ";
689  cutFlowVector_str[48] = "SRD-high: Derivation skim";
690  cutFlowVector_str[49] = "SRD-high: Lepton veto ";
691  cutFlowVector_str[50] = "SRD-high: Njets >= 4 ";
692  cutFlowVector_str[51] = "SRD-high: Nbjets >= 1 ";
693  cutFlowVector_str[52] = "SRD-high: met > 250 GeV ";
694  cutFlowVector_str[53] = "SRD-high: dPhi(jet,MET) > 0.4 ";
695  cutFlowVector_str[54] = "SRD-high: pT jet 1 > 80 GeV ";
696  cutFlowVector_str[55] = "SRD-high: pT jet 3 > 40 GeV ";
697  cutFlowVector_str[56] = "SRD-high: Njets >= 5 ";
698  cutFlowVector_str[57] = "SRD-high: pT jet 1 > 150 ";
699  cutFlowVector_str[58] = "SRD-high: pT jet 3 > 80 ";
700  cutFlowVector_str[59] = "SRD-high: pT jet 4 > 60 ";
701  cutFlowVector_str[60] = "SRD-high: mT(b,MET) min > 350 GeV ";
702  cutFlowVector_str[61] = "SRD-high: mT(b,MET) max > 450 GeV ";
703  cutFlowVector_str[62] = "SRD-high: Nbjets >=2 ";
704  cutFlowVector_str[63] = "SRD-high: met > 250 GeV ";
705  cutFlowVector_str[64] = "SRD-high: deltaR(b,b) > 0.8";
706  cutFlowVector_str[65] = "SRD-high: pT0b + pT1b > 400 GeV";
707  cutFlowVector_str[66] = "SRD-low: Njets >=5";
708  cutFlowVector_str[67] = "SRD-low: NBjets >=2";
709  cutFlowVector_str[68] = "SRD-low: met > 250 GeV";
710  cutFlowVector_str[69] = "SRD-low: mT(b,MET) min > 250 GeV ";
711  cutFlowVector_str[70] = "SRD-low: mT(b,MET) max > 300 GeV ";
712  cutFlowVector_str[71] = "SRD-low: deltaR(b,b) > 0.8";
713  cutFlowVector_str[72] = "SRD-low: pT jet 1 > 150 GeV ";
714  cutFlowVector_str[73] = "SRD-low: pT jet 3 > 100 GeV ";
715  cutFlowVector_str[74] = "SRD-low: pT jet 4 > 60 GeV ";
716  cutFlowVector_str[75] = "SRD-low: pT0b + pT1b > 300 GeV";
717 
718  // Cutflow for SRE
719  cutFlowVector_str[76] = "SRE: met > 550 GeV";
720  cutFlowVector_str[77] = "SRE: m jet0, R = 0.8 > 120 GeV";
721  cutFlowVector_str[78] = "SRE: m jet1, R = 0.8 > 80 GeV";
722  cutFlowVector_str[79] = "SRE: HT > 800 GeV";
723  cutFlowVector_str[80] = "SRE: met/sqrt(HT) > 18 GeV^1/2";
724  cutFlowVector_str[81] = "SRE: mT(b,MET) min > 200 GeV";
725  cutFlowVector_str[82] = "SRE: NBjets >=2";
726 
727  // Cutflow for SRC1
728 
729  cutFlowVector_str[83] = "SRC: Derivation skim";
730  cutFlowVector_str[84] = "SRC: Lepton veto ";
731  cutFlowVector_str[85] = "SRC: Njets >= 4 ";
732  cutFlowVector_str[86] = "SRC: Nbjets >= 1 ";
733  cutFlowVector_str[87] = "SRC: met > 250 GeV ";
734  cutFlowVector_str[88] = "SRC: dPhi(jet,MET) > 0.4 ";
735  cutFlowVector_str[89] = "SRC: pT jet 1 > 80 GeV ";
736  cutFlowVector_str[90] = "SRC: pT jet 3 > 40 GeV ";
737  cutFlowVector_str[91] = "SRC: NSbjet >=1";
738  cutFlowVector_str[92] = "SRC: NSjet >=5";
739  cutFlowVector_str[93] = "SRC: pT0sb > 40";
740  cutFlowVector_str[94] = "SRC: mS > 300";
741  cutFlowVector_str[95] = "SRC: dPhi(ISR,met) > 3";
742  cutFlowVector_str[96] = "SRC: pTISR > 400";
743  cutFlowVector_str[97] = "SRC: pT4S > 50";
744  cutFlowVector_str[98] = "SRC1: 0.30 <= R_ISR <= 0.40";
745  cutFlowVector_str[99] = "SRC2: 0.40 <= R_ISR <= 0.50";
746  cutFlowVector_str[100] = "SRC3: 0.50 <= R_ISR <= 0.60";
747  cutFlowVector_str[101] = "SRC4: 0.60 <= R_ISR <= 0.70";
748  cutFlowVector_str[102] = "SRC5: 0.70 <= R_ISR <= 0.80";
749 
750  int nElectrons=signalElectrons.size();
751  int nMuons=signalMuons.size();
752  int nJets=signalJets.size();
753 
754  bool cut_LeptonVeto=true;
755  if((nElectrons + nMuons)>0.)cut_LeptonVeto=false;
756 
757  double Ht=0;
758 
759  for(size_t jet=0;jet<signalJets.size();jet++)Ht+=signalJets[jet]->pT();
760 
761  double HtSig = Met/sqrt(Ht);
762 
763  bool devSkim = false;
764 
765  if( (Ht > 150.) || (signalElectrons.size() > 0 && signalElectrons[0]->pT() > 100.) || (signalElectrons.size() > 1 && signalElectrons[0]->pT() > 20. && signalElectrons[1]->pT() > 20.) || (signalMuons.size() > 0 && signalMuons[0]->pT() > 100.) || (signalMuons.size() > 1 && signalMuons[0]->pT() > 20. && signalMuons[1]->pT() > 20.) || (signalPhotons.size() > 0 && signalPhotons[0]->pT() > 100.) || (signalPhotons.size() > 1 && signalPhotons[0]->pT() > 50. && signalPhotons[1]->pT() > 50.))devSkim=true;
766 
767  bool cut_dPhiJetsPresel=false;
768  bool cut_dPhiJet2=false;
769  bool cut_dPhiJet1=false;
770  double dphi_jetmet1=9999;
771  if(nJets>0)dphi_jetmet1=std::acos(std::cos(signalJets.at(0)->phi()-metVec.phi()));
772  double dphi_jetmet2=9999;
773  if(nJets>1)dphi_jetmet2=std::acos(std::cos(signalJets.at(1)->phi()-metVec.phi()));
774  if(dphi_jetmet2>0.4)cut_dPhiJet2=true;
775  if(dphi_jetmet1>0.4)cut_dPhiJet1=true;
776  if(cut_dPhiJet1 && cut_dPhiJet2)cut_dPhiJetsPresel=true;
777 
778  bool cut_dPhiJet3=false;
779  bool cut_dPhiJets_AB=false;
780  double dphi_jetmet3=9999;
781  if(nJets>2)dphi_jetmet3=std::acos(std::cos(signalJets.at(2)->phi()-metVec.phi()));
782  if(dphi_jetmet3>0.4)cut_dPhiJet3=true;
783  if(cut_dPhiJetsPresel && cut_dPhiJet3)cut_dPhiJets_AB=true;
784 
785  //if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && DRBB > 1.2 && MtBMax > 200. && MtBMin > 200.)std::cout << "Met " << Met << " AntiKt12M_0 " << AntiKt12M_0 << " AntiKt12M_1 " << AntiKt12M_1 << " DRBB " << DRBB << " MtBMax " << MtBMax << " MtBMin " << MtBMin << std::endl;
786 
787 
788  for(int j=0;j<NCUTS;j++){
789  if(
790  (j==0) ||
791 
792  (j==1 && devSkim) ||
793 
794  (j==2 && devSkim && cut_LeptonVeto) ||
795 
796  (j==3 && devSkim && cut_LeptonVeto && signalJets.size()>3) ||
797 
798  (j==4 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0) ||
799 
800  (j==5 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250.) ||
801 
802  (j==6 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB) ||
803 
804  (j==7 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80.) ||
805 
806  (j==8 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40.) ||
807 
808  (j==9 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120.) ||
809 
810  // SRA-TT
811 
812  (j==10 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120.) ||
813 
814  (j==11 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 400. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120.) ||
815 
816  (j==12 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 400. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && AntiKt8M_0>60.) ||
817 
818  (j==13 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 400. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && AntiKt8M_0>60. && MtBMin > 200.) ||
819 
820  (j==14 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 400. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && AntiKt8M_0>60. && MtBMin > 200. && DRBB > 1.) ||
821 
822  (j==15 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 400. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && AntiKt8M_0>60. && MtBMin > 200. && DRBB > 1. && MT2Chi2>400.) ||
823 
824  (j==16 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 400. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && AntiKt8M_0>60. && MtBMin > 200. && DRBB > 1. && MT2Chi2>400. && NBJets>=2) ||
825 
826  // SRA-TW
827 
828  /* cutFlowVector_str[17] = "SRA-TW: m jet1, R=1.2 < 120 GeV";
829  cutFlowVector_str[18] = "SRA-TW: m jet1, R=1.2 > 60 GeV";
830  cutFlowVector_str[19] = "SRA-TW: met > 500 GeV ";
831  cutFlowVector_str[20] = "SRA-TW: m jet0, R=0.8 > 60 GeV";
832  cutFlowVector_str[21] = "SRA-TW: mT(b,MET) min > 200 GeV";
833  cutFlowVector_str[22] = "SRA-TW: mT2 > 400 GeV ";*/
834 
835  (j==17 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120.) ||
836 
837  (j==18 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. && AntiKt12M_1>60.) ||
838 
839  (j==19 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. && AntiKt12M_1>60. && Met > 500.) ||
840 
841  (j==20 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. && AntiKt12M_1>60. && Met > 500. && AntiKt8M_0>60.) ||
842 
843  (j==21 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. && AntiKt12M_1>60. && Met > 500. && AntiKt8M_0>60. && MtBMin > 200.) ||
844 
845  (j==22 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. && AntiKt12M_1>60. && Met > 500. && AntiKt8M_0>60. && MtBMin > 200. && MT2Chi2>400.) ||
846 
847  (j==23 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. && AntiKt12M_1>60. && Met > 500. && AntiKt8M_0>60. && MtBMin > 200. && MT2Chi2>400.) ||
848 
849  /* cutFlowVector_str[24] = "SRA-T0: m jet1, R=1.2 < 60 GeV";
850  cutFlowVector_str[25] = "SRA-T0: m jet0, R=0.8 > 60 GeV";
851  cutFlowVector_str[26] = "SRA-T0: met > 550 GeV ";
852  cutFlowVector_str[27] = "SRA-T0: mT(b,MET) min > 200 GeV";
853  cutFlowVector_str[28] = "SRA-T0: mT2 > 500 GeV ";
854  cutFlowVector_str[29] = "SRA-T0: Nbjets >=2 "; */
855 
856 
857  (j==24 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0 > 120. && AntiKt12M_1<60.) ||
858 
859  (j==25 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0 > 120. && AntiKt12M_1<60. && AntiKt8M_0>60.) ||
860 
861  (j==26 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0 > 120. && AntiKt12M_1<60. && AntiKt8M_0>60. && Met > 550.) ||
862 
863  (j==27 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0 > 120. && AntiKt12M_1<60. && AntiKt8M_0>60. && Met > 550. && MtBMin > 200.) ||
864 
865  (j==28 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0 > 120. && AntiKt12M_1<60. && AntiKt8M_0>60. && Met > 550. && MtBMin > 200. && MT2Chi2 > 500.) ||
866 
867  (j==29 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0 > 120. && AntiKt12M_1<60. && AntiKt8M_0>60. && Met > 550. && MtBMin > 200. && MT2Chi2 > 500.) ||
868 
869  /* cutFlowVector_str[30] = "SRB-TT: m jet1, R=1.2 > 120 GeV";
870  cutFlowVector_str[31] = "SRB-TT: deltaR(b,b) > 1.2";
871  cutFlowVector_str[32] = "SRB-TT: mT(b,MET) max > 200 GeV";
872  cutFlowVector_str[33] = "SRB-TT: mT(b,MET) min > 200 GeV";
873  cutFlowVector_str[34] = "SRB-TT: Nbjets >=2 ";*/
874 
875  (j==30 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120.) ||
876 
877  (j==31 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && DRBB > 1.2) ||
878 
879  (j==32 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && DRBB > 1.2 && MtBMax > 200.) ||
880 
881  (j==33 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && DRBB > 1.2 && MtBMax > 200. && MtBMin > 200.) ||
882 
883  (j==34 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && DRBB > 1.2 && MtBMax > 200. && MtBMin > 200.) ||
884 
885  /* cutFlowVector_str[35] = "SRB-TW: m jet1, R=1.2 < 120 GeV";
886  cutFlowVector_str[36] = "SRB-TW: m jet1, R=1.2 > 60 GeV";
887  cutFlowVector_str[37] = "SRB-TW: deltaR(b,b) > 1.2";
888  cutFlowVector_str[38] = "SRB-TW: mT(b,MET) max > 200 GeV";
889  cutFlowVector_str[39] = "SRB-TW: mT(b,MET) min > 200 GeV";
890  cutFlowVector_str[40] = "SRB-TW: Nbjets >=2 ";*/
891 
892  (j==35 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120.) ||
893 
894  (j==36 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. &&AntiKt12M_1>60.) ||
895 
896  (j==37 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. &&AntiKt12M_1>60. && DRBB > 1.2) ||
897 
898  (j==38 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. &&AntiKt12M_1>60. && DRBB > 1.2 && MtBMax > 200.) ||
899 
900  (j==39 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. &&AntiKt12M_1>60. && DRBB > 1.2 && MtBMax > 200. && MtBMin > 200.) ||
901 
902  (j==40 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. &&AntiKt12M_1>60. && DRBB > 1.2 && MtBMax > 200. && MtBMin > 200.) ||
903 
904  /* cutFlowVector_str[41] = "SRB-T0: m jet1, R=1.2 < 60 GeV";
905  cutFlowVector_str[42] = "SRB-T0: mT(b,MET) min > 200 GeV";
906  cutFlowVector_str[43] = "SRB-T0: deltaR(b,b) > 1.2";
907  cutFlowVector_str[44] = "SRB-T0: mT(b,MET) max > 200 GeV";
908  cutFlowVector_str[45] = "SRB-T0: met > 250 GeV ";
909  cutFlowVector_str[46] = "SRB-T0: Nbjets >=2 ";*/
910 
911  (j==41 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<60.) ||
912 
913  (j==42 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<60. && MtBMin > 200.) ||
914 
915  (j==43 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<60. && MtBMin > 200. && DRBB > 1.2) ||
916 
917  (j==44 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<60. && MtBMin > 200. && DRBB > 1.2 && MtBMax > 200.) ||
918 
919  (j==45 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<60. && MtBMin > 200. && DRBB > 1.2 && MtBMax > 200.) ||
920 
921  (j==46 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<60. && MtBMin > 200. && DRBB > 1.2 && MtBMax > 200.) ||
922 
923  // Cutflow for SRD
924  /*cutFlowVector_str[47] = "SRD-high: No cuts ";
925  cutFlowVector_str[48] = "SRD-high: Derivation skim";
926  cutFlowVector_str[49] = "SRD-high: Lepton veto ";
927  cutFlowVector_str[50] = "SRD-high: Njets >= 4 ";
928  cutFlowVector_str[51] = "SRD-high: Nbjets >= 1 ";
929  cutFlowVector_str[52] = "SRD-high: met > 250 GeV ";
930  cutFlowVector_str[53] = "SRD-high: dPhi(jet,MET) > 0.4 ";
931  cutFlowVector_str[54] = "SRD-high: pT jet 1 > 80 GeV ";
932  cutFlowVector_str[55] = "SRD-high: pT jet 3 > 40 GeV ";
933  cutFlowVector_str[56] = "SRD-high: Njets >= 5 ";
934  cutFlowVector_str[57] = "SRD-high: pT jet 1 > 150 ";
935  cutFlowVector_str[58] = "SRD-high: pT jet 3 > 80 ";
936  cutFlowVector_str[59] = "SRD-high: pT jet 4 > 60 ";
937  cutFlowVector_str[60] = "SRD-high: mT(b,MET) min > 350 GeV ";
938  cutFlowVector_str[61] = "SRD-high: mT(b,MET) max > 450 GeV ";
939  cutFlowVector_str[62] = "SRD-high: Nbjets >=2 ";
940  cutFlowVector_str[63] = "SRD-high: met > 250 GeV ";
941  cutFlowVector_str[64] = "SRD-high: deltaR(b,b) > 0.8";
942  cutFlowVector_str[65] = "SRD-high: pT0b + pT1b > 400 GeV";*/
943 
944  (j==47) ||
945 
946  (j==48 && devSkim) ||
947 
948  (j==49 && devSkim && cut_LeptonVeto) ||
949 
950  (j==50 && devSkim && cut_LeptonVeto && signalJets.size()>3) ||
951 
952  (j==51 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0) ||
953 
954  (j==52 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250.) ||
955 
956  (j==53 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB) ||
957 
958  (j==54 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80.) ||
959 
960  (j==55 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. ) ||
961 
962  (j==56 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. ) ||
963 
964  (j==57 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>40. ) ||
965 
966  (j==58 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>80. ) ||
967 
968  (j==59 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>80. && signalJets[4]->pT()>60.) ||
969 
970  (j==60 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>80. && signalJets[4]->pT()>60. && MtBMin > 350.) ||
971 
972  (j==61 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>80. && signalJets[4]->pT()>60. && MtBMin > 350. && MtBMax > 450.) ||
973 
974  (j==62 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>80. && signalJets[4]->pT()>60. && MtBMin > 350. && MtBMax > 450.) ||
975 
976  (j==63 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>80. && signalJets[4]->pT()>60. && MtBMin > 350. && MtBMax > 450.) ||
977 
978  (j==64 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>80. && signalJets[4]->pT()>60. && MtBMin > 350. && MtBMax > 450. && DRBB > 0.8) ||
979 
980  (j==65 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>80. && signalJets[4]->pT()>60. && MtBMin > 350. && MtBMax > 450. && DRBB > 0.8 && ( (signalBJets[0]->pT() + signalBJets[1]->pT())>400.)) ||
981 
982  /*cutFlowVector_str[66] = "SRD-low: Njets >=5";
983  cutFlowVector_str[67] = "SRD-low: NBjets >=2";
984  cutFlowVector_str[68] = "SRD-low: met > 250 GeV";
985  cutFlowVector_str[69] = "SRD-low: mT(b,MET) min > 250 GeV ";
986  cutFlowVector_str[70] = "SRD-low: mT(b,MET) max > 300 GeV ";
987  cutFlowVector_str[71] = "SRD-low: deltaR(b,b) > 0.8";
988  cutFlowVector_str[72] = "SRD-low: pT jet 1 > 150 GeV ";
989  cutFlowVector_str[73] = "SRD-low: pT jet 3 > 100 GeV ";
990  cutFlowVector_str[74] = "SRD-low: pT jet 4 > 60 GeV ";
991  cutFlowVector_str[75] = "SRD-low: pT0b + pT1b > 300 GeV";*/
992 
993 
994  (j==66 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40.) ||
995 
996  (j==67 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40.) ||
997 
998  (j==68 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40.) ||
999 
1000  (j==69 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && MtBMin > 250.) ||
1001 
1002  (j==70 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && MtBMin > 250. && MtBMax > 300.) ||
1003 
1004  (j==71 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && MtBMin > 250. && MtBMax > 300. && DRBB > 0.8) ||
1005 
1006  (j==72 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>40. && MtBMin > 250. && MtBMax > 300. && DRBB > 0.8) ||
1007 
1008  (j==73 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>100. && MtBMin > 250. && MtBMax > 300. && DRBB > 0.8) ||
1009 
1010  (j==74 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>100. && signalJets[4]->pT()>60. && MtBMin > 250. && MtBMax > 300. && DRBB > 0.8) ||
1011 
1012  (j==75 && devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>100. && signalJets[4]->pT()>60. && MtBMin > 250. && MtBMax > 300. && DRBB > 0.8 && ( (signalBJets[0]->pT() + signalBJets[1]->pT())>300.)) ||
1013 
1014  /* cutFlowVector_str[76] = "SRE: met > 550 GeV";
1015  cutFlowVector_str[77] = "SRE: m jet0, R = 0.8 > 120 GeV";
1016  cutFlowVector_str[78] = "SRE: m jet1, R = 0.8 > 80 GeV";
1017  cutFlowVector_str[79] = "SRE: HT > 800 GeV";
1018  cutFlowVector_str[80] = "SRE: met/sqrt(HT) > 18 GeV^1/2";
1019  cutFlowVector_str[81] = "SRE: mT(b,MET) min > 200 GeV";
1020  cutFlowVector_str[82] = "SRE: NBjets >=2";*/
1021 
1022  (j==76 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 550. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40.) ||
1023 
1024  (j==77 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 550. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt8M_0 > 120.) ||
1025 
1026  (j==78 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 550. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt8M_0 > 120. && AntiKt8M_1 > 80.) ||
1027 
1028  (j==79 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 550. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt8M_0 > 120. && AntiKt8M_1 > 80. && Ht > 800.) ||
1029 
1030  (j==80 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 550. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt8M_0 > 120. && AntiKt8M_1 > 80. && Ht > 800. && HtSig > 18.) ||
1031 
1032  (j==81 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 550. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt8M_0 > 120. && AntiKt8M_1 > 80. && Ht > 800. && HtSig > 18. && MtBMin > 200.) ||
1033 
1034  (j==82 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && Met > 550. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt8M_0 > 120. && AntiKt8M_1 > 80. && Ht > 800. && HtSig > 18. && MtBMin > 200.)
1035 
1036  /*cutFlowVector_str[83] = "SRC: Derivation skim";
1037  cutFlowVector_str[84] = "SRC: Lepton veto ";
1038  cutFlowVector_str[85] = "SRC: Njets >= 4 ";
1039  cutFlowVector_str[86] = "SRC: Nbjets >= 1 ";
1040  cutFlowVector_str[87] = "SRC: met > 250 GeV ";
1041  cutFlowVector_str[88] = "SRC: dPhi(jet,MET) > 0.4 ";
1042  cutFlowVector_str[89] = "SRC: pT jet 1 > 80 GeV ";
1043  cutFlowVector_str[90] = "SRC: pT jet 3 > 40 GeV ";
1044  cutFlowVector_str[91] = "SRC: NSbjet >=1";
1045  cutFlowVector_str[92] = "SRC: NSjet >=5";
1046  cutFlowVector_str[93] = "SRC: pT0sb > 40";
1047  cutFlowVector_str[94] = "SRC: mS > 300";
1048  cutFlowVector_str[95] = "SRC: dPhi(ISR,met) > 3";
1049  cutFlowVector_str[96] = "SRC: pTISR > 400";
1050  cutFlowVector_str[97] = "SRC: pT4S > 50";
1051  cutFlowVector_str[98] = "SRC1: 0.30 <= R_ISR <= 0.40";
1052  cutFlowVector_str[99] = "SRC2: 0.40 <= R_ISR <= 0.50";
1053  cutFlowVector_str[100] = "SRC3: 0.50 <= R_ISR <= 0.60";
1054  cutFlowVector_str[101] = "SRC4: 0.60 <= R_ISR <= 0.70";
1055  cutFlowVector_str[102] = "SRC5: 0.70 <= R_ISR <= 0.80";*/
1056 
1057  /*(j==83 && devSkim) ||
1058 
1059  (j==84 && devSkim && cut_LeptonVeto) ||
1060 
1061  (j==85 && devSkim && cut_LeptonVeto && signalJets.size()>3) ||
1062 
1063  (j==86 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0) ||
1064 
1065  (j==87 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250.) ||
1066 
1067  (j==88 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB) ||
1068 
1069  (j==89 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80.) ||
1070 
1071  (j==90 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. ) ||
1072 
1073  (j==91 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1) ||
1074 
1075  (j==92 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5) ||
1076 
1077  (j==93 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40) ||
1078 
1079  (j==94 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300) ||
1080 
1081  (j==95 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00) ||
1082 
1083  (j==96 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400) ||
1084 
1085  (j==97 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400 && CA_pTjV4 > 50) ||
1086 
1087  (j==98 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400 && CA_pTjV4 > 50 && CA_RISR >= 0.30 && CA_RISR <= 0.4) ||
1088 
1089  (j==99 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400 && CA_pTjV4 > 50 && CA_RISR >= 0.40 && CA_RISR <= 0.5) ||
1090 
1091  (j==100 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400 && CA_pTjV4 > 50 && CA_RISR >= 0.50 && CA_RISR <= 0.6) ||
1092 
1093  (j==101 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400 && CA_pTjV4 > 50 && CA_RISR >= 0.60 && CA_RISR <= 0.7) ||
1094 
1095  (j==102 && devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400 && CA_pTjV4 > 50 && CA_RISR >= 0.70 && CA_RISR <= 0.8) */
1096 
1097 
1098  ){
1099 
1100  cutFlowVector[j]++;
1101  }
1102 
1103  }
1104 
1105 
1106  if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 400. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && AntiKt8M_0>60. && MtBMin > 200. && DRBB > 1. && MT2Chi2>400. && NBJets>=2)isSRA_TT=true;
1107 
1108  if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. && AntiKt12M_1>60. && Met > 500. && AntiKt8M_0>60. && MtBMin > 200. && MT2Chi2>400.)isSRA_TW=true;
1109 
1110  if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0 > 120. && AntiKt12M_1<60. && AntiKt8M_0>60. && Met > 550. && MtBMin > 200. && MT2Chi2 > 500.)isSRA_T0=true;
1111 
1112  if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1>120. && DRBB > 1.2 && MtBMax > 200. && MtBMin > 200.)isSRB_TT=true;
1113 
1114  if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<120. &&AntiKt12M_1>60. && DRBB > 1.2 && MtBMax > 200. && MtBMin > 200.)isSRB_TW=true;
1115 
1116  if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt12M_0>120. && AntiKt12M_1<60. && MtBMin > 200. && DRBB > 1.2 && MtBMax > 200.)isSRB_T0=true;
1117 
1118  /*if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400 && CA_pTjV4 > 50 && CA_RISR >= 0.30 && CA_RISR <= 0.4)isSRC1=true;
1119 
1120  if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400 && CA_pTjV4 > 50 && CA_RISR >= 0.40 && CA_RISR <= 0.5)isSRC2=true;
1121 
1122  if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400 && CA_pTjV4 > 50 && CA_RISR >= 0.50 && CA_RISR <= 0.6)isSRC3=true;
1123 
1124  if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400 && CA_pTjV4 > 50 && CA_RISR >= 0.60 && CA_RISR <= 0.7)isSRC4=true;
1125 
1126  if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>0 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && CA_NbV >= 1 && CA_NjV >= 5 && CA_pTbV1 > 40 && CA_MS > 300 && CA_dphiISRI > 3.00 && CA_PTISR > 400 && CA_pTjV4 > 50 && CA_RISR >= 0.70 && CA_RISR <= 0.8)isSRC5=true;*/
1127 
1128  if( devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>80. && signalJets[4]->pT()>60. && MtBMin > 350. && MtBMax > 450. && DRBB > 0.8 && ( (signalBJets[0]->pT() + signalBJets[1]->pT())>400.))isSRD_high=true;
1129 
1130  if(devSkim && cut_LeptonVeto && signalJets.size()>4 && signalBJets.size()>1 && Met > 250. && cut_dPhiJets_AB && signalJets[1]->pT()>150. && signalJets[3]->pT()>100. && signalJets[4]->pT()>60. && MtBMin > 250. && MtBMax > 300. && DRBB > 0.8 && ( (signalBJets[0]->pT() + signalBJets[1]->pT())>300.))isSRD_low=true;
1131 
1132  if(devSkim && cut_LeptonVeto && signalJets.size()>3 && signalBJets.size()>1 && Met > 550. && cut_dPhiJets_AB && signalJets[1]->pT()>80. && signalJets[3]->pT()>40. && AntiKt8M_0 > 120. && AntiKt8M_1 > 80. && Ht > 800. && HtSig > 18. && MtBMin > 200.)isSRE=true;
1133 
1134 
1135  if(isSRA_TT) _counters.at("SRA_TT").add_event(event);
1136  if(isSRA_TW) _counters.at("SRA_TW").add_event(event);
1137  if(isSRA_T0) _counters.at("SRA_T0").add_event(event);
1138  if(isSRB_TT) _counters.at("SRB_TT").add_event(event);
1139  if(isSRB_TW) _counters.at("SRB_TW").add_event(event);
1140  if(isSRB_T0) _counters.at("SRB_T0").add_event(event);
1141  if(isSRC1) _counters.at("SRC1").add_event(event);
1142  if(isSRC2) _counters.at("SRC2").add_event(event);
1143  if(isSRC3) _counters.at("SRC3").add_event(event);
1144  if(isSRC4) _counters.at("SRC4").add_event(event);
1145  if(isSRC5) _counters.at("SRC5").add_event(event);
1146  if(isSRD_low) _counters.at("SRD_low").add_event(event);
1147  if(isSRD_high) _counters.at("SRD_high").add_event(event);
1148  if(isSRE) _counters.at("SRE").add_event(event);
1149 
1150  return;
1151 
1152  }
1153 
1155  void combine(const Analysis* other)
1156  {
1157  const Analysis_ATLAS_13TeV_0LEPStop_36invfb* specificOther
1158  = dynamic_cast<const Analysis_ATLAS_13TeV_0LEPStop_36invfb*>(other);
1159 
1160  for (auto& pair : _counters) { pair.second += specificOther->_counters.at(pair.first); }
1161 
1162  if (NCUTS != specificOther->NCUTS) NCUTS = specificOther->NCUTS;
1163  for (int j=0; j<NCUTS; j++) {
1164  cutFlowVector[j] += specificOther->cutFlowVector[j];
1165  cutFlowVector_str[j] = specificOther->cutFlowVector_str[j];
1166  }
1167  }
1168 
1169 
1170  void collect_results() {
1171 
1172  // double scale_by=1.;
1173  // cout << "------------------------------------------------------------------------------------------------------------------------------ "<<endl;
1174  // cout << "CUT FLOW: ATLAS 13 TeV 0 lep stop paper "<<endl;
1175  // cout << "------------------------------------------------------------------------------------------------------------------------------"<<endl;
1176  // cout<< right << setw(40) << "CUT" << setw(20) << "RAW" << setw(20) << "SCALED"
1177  // << setw(20) << "%" << setw(20) << "clean adj RAW"<< setw(20) << "clean adj %" << endl;
1178  // for (int j=0; j<NCUTS; j++) {
1179  // cout << right << setw(40) << cutFlowVector_str[j].c_str() << setw(20)
1180  // << cutFlowVector[j] << setw(20) << cutFlowVector[j]*scale_by << setw(20)
1181  // << 100.*cutFlowVector[j]/cutFlowVector[0] << "%" << setw(20)
1182  // << cutFlowVector[j]*scale_by << setw(20) << 100.*cutFlowVector[j]/cutFlowVector[0]<< "%" << endl;
1183  // }
1184  // cout << "------------------------------------------------------------------------------------------------------------------------------ "<<endl;
1185 
1187 
1188  add_result(SignalRegionData(_counters.at("SRA_TT"), 11, {8.6, 2.1}));
1189  add_result(SignalRegionData(_counters.at("SRA_TW"), 9, {9.3, 2.2}));
1190  add_result(SignalRegionData(_counters.at("SRA_T0"), 18, {18.7, 2.7}));
1191  add_result(SignalRegionData(_counters.at("SRB_TT"), 38, { 39.3, 7.6}));
1192  add_result(SignalRegionData(_counters.at("SRB_TW"), 53, {52.4, 7.4}));
1193  add_result(SignalRegionData(_counters.at("SRB_T0"), 206, { 179., 26.}));
1194 
1195  // MJW removes the recursive jigsaw signal regions for the Feb 2018 SUSY scans
1196  // The ISR modelling in Pythia does not give reliable answers
1197 
1198  /*
1199  add_result(SignalRegionData(_counters.at("SRC1"), 20, { 20.6, 6.5}));
1200  add_result(SignalRegionData(_counters.at("SRC2"), 22, { 27.6, 4.9}));
1201  add_result(SignalRegionData(_counters.at("SRC3"), 22, { 18.9, 3.4}));
1202  add_result(SignalRegionData(_counters.at("SRC4"), 1, { 7.7, 1.2}));
1203  add_result(SignalRegionData(_counters.at("SRC5"), 0, { 0.91, 0.73}));
1204  */
1205 
1206  add_result(SignalRegionData(_counters.at("SRD_low"), 27, { 25.1, 6.2}));
1207  add_result(SignalRegionData(_counters.at("SRD_high"), 11, { 8.5,1.5}));
1208  add_result(SignalRegionData(_counters.at("SRE"), 3, { 3.64,0.79}));
1209 
1210  return;
1211  }
1212 
1213 
1214  protected:
1215  void analysis_specific_reset() {
1216  for (auto& pair : _counters) { pair.second.reset(); }
1217 
1218  std::fill(cutFlowVector.begin(), cutFlowVector.end(), 0);
1219  }
1220 
1221  };
1222 
1223 
1224  DEFINE_ANALYSIS_FACTORY(ATLAS_13TeV_0LEPStop_36invfb)
1225 
1226 
1227  }
1228 }
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
bool sortByMass_sharedptr(std::shared_ptr< HEPUtils::Jet > jet1, std::shared_ptr< HEPUtils::Jet > jet2)
#define DEFINE_ANALYSIS_FACTORY(ANAME)
For analysis factory function definition.
Definition: Analysis.hpp:124
bool sortByPT13_sharedptr(std::shared_ptr< HEPUtils::Jet > jet1, std::shared_ptr< HEPUtils::Jet > jet2)
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 *)
const double mt
Definition: topness.h:39
bool sortByPT13(const HEPUtils::Jet *jet1, const HEPUtils::Jet *jet2)
STL namespace.
START_MODEL b
Definition: demo.hpp:270
bool sortByMass(const HEPUtils::Jet *jet1, const HEPUtils::Jet *jet2)
double calcMT(HEPUtils::P4 jetMom, HEPUtils::P4 metMom)
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.
DS5_MSPCTM DS_INTDOF int
void set_mn(double mn)
Definition: mt2_bisect.cpp:130
std::vector< std::shared_ptr< HEPUtils::Jet > > get_jets(const std::vector< MOM *> &moms, double R, double ptmin=0 *GeV, FJNS::JetAlgorithm alg=FJNS::antikt_algorithm)
Run jet clustering from any P4-compatible momentum type.
Definition: Utils.hpp:269
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
Class for ColliderBit analyses.
Functions that do super fast ATLAS detector simulation based on four-vector smearing.