gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp
Go to the documentation of this file.
1 
2 #include "gambit/cmake/cmake_variables.hpp"
3 #ifndef EXCLUDE_ROOT
4 #ifndef EXCLUDE_RESTFRAMES
5 
6 #include <vector>
7 #include <cmath>
8 #include <memory>
9 #include <algorithm>
10 #include <iomanip>
11 
12 #include "gambit/ColliderBit/analyses/Analysis.hpp"
13 #include "gambit/ColliderBit/ATLASEfficiencies.hpp"
14 #include "gambit/ColliderBit/mt2_bisect.h"
15 
16 #include "RestFrames/RestFrames.hh"
17 #include "TLorentzVector.h"
18 
19 using namespace std;
20 
21 /* The ATLAS 13 TeV 3 lepton low mass recursive jigsaw search
22 
23  Based on code kindly supplied by Abhishek Sharma
24 
25  Note that use of ROOT is compulsory for the RestFrames package
26 
27  Based on: https://arxiv.org/pdf/1806.02293.pdf
28 
29  Code adapted by Martin White
30 
31  KNOWN ISSUES
32 
33  1) Need to check overlap removal step when the paper comes out. For now, have assumed it is the same as the stop analysis.
34 
35 */
36 
37 namespace Gambit {
38  namespace ColliderBit {
39 
40  bool sortByPT_RJ3L(HEPUtils::Jet* jet1, HEPUtils::Jet* jet2) { return (jet1->pT() > jet2->pT()); }
41  bool sortLepByPT_RJ3L(HEPUtils::Particle* lep1, HEPUtils::Particle* lep2) { return (lep1->pT() > lep2->pT());}
42  //bool sortByMass(HEPUtils::Jet* jet1, HEPUtils::Jet* jet2) { return (jet1->mass() > jet2->mass()); }
43 
44  bool SortLeptons(const pair<TLorentzVector,int> lv1, const pair<TLorentzVector,int> lv2)
45  //bool VariableConstruction::SortLeptons(const lep lv1, const lep lv2)
46  {
47  return lv1.first.Pt() > lv2.first.Pt();
48  }
49 
50  bool SortJets(const TLorentzVector jv1, const TLorentzVector jv2)
51  {
52  return jv1.Pt() > jv2.Pt();
53  }
54 
55 
56  class Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb : public Analysis {
57 
58  protected:
59  // Numbers passing cuts
60  int _num2L2JHIGH;
61  int _num2L2JINT;
62  int _num2L2JLOW;
63  int _num2L2JCOMP;
64  int _num3LHIGH;
65  int _num3LINT;
66  int _num3LLOW;
67  int _num3LCOMP;
68 
69  private:
70 
71  vector<int> cutFlowVector;
72  vector<string> cutFlowVector_str;
73  int NCUTS; //=16;
74 
75  // Recursive jigsaw objects (using RestFrames)
76 
80 
81  unique_ptr<RestFrames::LabRecoFrame> LAB_2L2J;
82  unique_ptr<RestFrames::DecayRecoFrame> C1N2_2L2J;
83  unique_ptr<RestFrames::DecayRecoFrame> C1a_2L2J;
84  unique_ptr<RestFrames::DecayRecoFrame> N2b_2L2J;
85 
86  unique_ptr<RestFrames::DecayRecoFrame> Wa_2L2J;
87  unique_ptr<RestFrames::DecayRecoFrame> Zb_2L2J;
88 
89  unique_ptr<RestFrames::VisibleRecoFrame> J1_2L2J;
90  unique_ptr<RestFrames::VisibleRecoFrame> J2_2L2J;
91  unique_ptr<RestFrames::VisibleRecoFrame> L1_2L2J;
92  unique_ptr<RestFrames::VisibleRecoFrame> L2_2L2J;
93 
94  unique_ptr<RestFrames::InvisibleRecoFrame> X1a_2L2J;
95  unique_ptr<RestFrames::InvisibleRecoFrame> X1b_2L2J;
96 
97  unique_ptr<RestFrames::InvisibleGroup> INV_2L2J;
98 
99  unique_ptr<RestFrames::SetMassInvJigsaw> X1_mass_2L2J;
100  unique_ptr<RestFrames::SetRapidityInvJigsaw> X1_eta_2L2J;
101 
102  unique_ptr<RestFrames::ContraBoostInvJigsaw> X1X1_contra_2L2J;
103 
107 
108  unique_ptr<RestFrames::LabRecoFrame> LAB_3L;
109  unique_ptr<RestFrames::DecayRecoFrame> C1N2_3L;
110  unique_ptr<RestFrames::DecayRecoFrame> C1a_3L;
111  unique_ptr<RestFrames::DecayRecoFrame> N2b_3L;
112 
113  // unique_ptr<RestFrames::DecayRecoFrame> Wa_3L;
114  // unique_ptr<RestFrames::DecayRecoFrame> Zb_3L;
115 
116  unique_ptr<RestFrames::VisibleRecoFrame> L1a_3L;
117  unique_ptr<RestFrames::VisibleRecoFrame> L1b_3L;
118  unique_ptr<RestFrames::VisibleRecoFrame> L2b_3L;
119 
120  unique_ptr<RestFrames::InvisibleRecoFrame> X1a_3L;
121  unique_ptr<RestFrames::InvisibleRecoFrame> X1b_3L;
122 
123  unique_ptr<RestFrames::InvisibleGroup> INV_3L;
124 
125  unique_ptr<RestFrames::SetMassInvJigsaw> X1_mass_3L;
126  unique_ptr<RestFrames::SetRapidityInvJigsaw> X1_eta_3L;
127 
128  unique_ptr<RestFrames::ContraBoostInvJigsaw> X1X1_contra_3L;
129 
130  // combinatoric (transverse) tree
131  // for jet assignment
132  unique_ptr<RestFrames::LabRecoFrame> LAB_comb;
133  unique_ptr<RestFrames::DecayRecoFrame> CM_comb;
134  unique_ptr<RestFrames::DecayRecoFrame> S_comb;
135  unique_ptr<RestFrames::VisibleRecoFrame> ISR_comb;
136  unique_ptr<RestFrames::VisibleRecoFrame> J_comb;
137  unique_ptr<RestFrames::VisibleRecoFrame> L_comb;
138  unique_ptr<RestFrames::InvisibleRecoFrame> I_comb;
139  unique_ptr<RestFrames::InvisibleGroup> INV_comb;
140  unique_ptr<RestFrames::SetMassInvJigsaw> InvMass_comb;
141  unique_ptr<RestFrames::CombinatoricGroup> JETS_comb;
142  unique_ptr<RestFrames::MinMassesCombJigsaw> SplitJETS_comb;
143 
144  // 2L+NJ tree (Z->ll + W/Z->qq)
145  unique_ptr<RestFrames::LabRecoFrame> LAB_2LNJ;
146  unique_ptr<RestFrames::DecayRecoFrame> CM_2LNJ;
147  unique_ptr<RestFrames::DecayRecoFrame> S_2LNJ;
148  unique_ptr<RestFrames::VisibleRecoFrame> ISR_2LNJ;
149 
150  unique_ptr<RestFrames::DecayRecoFrame> Ca_2LNJ;
151  unique_ptr<RestFrames::DecayRecoFrame> Z_2LNJ;
152  unique_ptr<RestFrames::VisibleRecoFrame> L1_2LNJ;
153  unique_ptr<RestFrames::VisibleRecoFrame> L2_2LNJ;
154 
155  unique_ptr<RestFrames::DecayRecoFrame> Cb_2LNJ;
156  unique_ptr<RestFrames::SelfAssemblingRecoFrame> JSA_2LNJ;
157  unique_ptr<RestFrames::VisibleRecoFrame> J_2LNJ;
158 
159  unique_ptr<RestFrames::InvisibleRecoFrame> Ia_2LNJ;
160  unique_ptr<RestFrames::InvisibleRecoFrame> Ib_2LNJ;
161 
162  unique_ptr<RestFrames::InvisibleGroup> INV_2LNJ;
163  unique_ptr<RestFrames::SetMassInvJigsaw> InvMass_2LNJ;
164  unique_ptr<RestFrames::SetRapidityInvJigsaw> InvRapidity_2LNJ;
165  unique_ptr<RestFrames::ContraBoostInvJigsaw> SplitINV_2LNJ;
166  unique_ptr<RestFrames::CombinatoricGroup> JETS_2LNJ;
167 
168  // 2L+1L tree (Z->ll + Z/W->l)
169  unique_ptr<RestFrames::LabRecoFrame> LAB_2L1L;
170  unique_ptr<RestFrames::DecayRecoFrame> CM_2L1L;
171  unique_ptr<RestFrames::DecayRecoFrame> S_2L1L;
172  unique_ptr<RestFrames::VisibleRecoFrame> ISR_2L1L;
173 
174  unique_ptr<RestFrames::DecayRecoFrame> Ca_2L1L;
175  unique_ptr<RestFrames::DecayRecoFrame> Z_2L1L;
176  unique_ptr<RestFrames::VisibleRecoFrame> L1_2L1L;
177  unique_ptr<RestFrames::VisibleRecoFrame> L2_2L1L;
178 
179  unique_ptr<RestFrames::DecayRecoFrame> Cb_2L1L;
180  unique_ptr<RestFrames::VisibleRecoFrame> Lb_2L1L;
181 
182  unique_ptr<RestFrames::InvisibleRecoFrame> Ia_2L1L;
183  unique_ptr<RestFrames::InvisibleRecoFrame> Ib_2L1L;
184 
185  unique_ptr<RestFrames::InvisibleGroup> INV_2L1L;
186  unique_ptr<RestFrames::SetMassInvJigsaw> InvMass_2L1L;
187  unique_ptr<RestFrames::SetRapidityInvJigsaw> InvRapidity_2L1L;
188  unique_ptr<RestFrames::ContraBoostInvJigsaw> SplitINV_2L1L;
189 
190  // Debug histos
191 
192  void JetLeptonOverlapRemoval(vector<HEPUtils::Jet*> &jetvec, vector<HEPUtils::Particle*> &lepvec, double DeltaRMax) {
193  //Routine to do jet-lepton check
194  //Discards jets if they are within DeltaRMax of a lepton
195 
196  vector<HEPUtils::Jet*> Survivors;
197 
198  for(unsigned int itjet = 0; itjet < jetvec.size(); itjet++) {
199  bool overlap = false;
200  HEPUtils::P4 jetmom=jetvec.at(itjet)->mom();
201  for(unsigned int itlep = 0; itlep < lepvec.size(); itlep++) {
202  HEPUtils::P4 lepmom=lepvec.at(itlep)->mom();
203  double dR;
204 
205  dR=jetmom.deltaR_eta(lepmom);
206 
207  if(fabs(dR) <= DeltaRMax) overlap=true;
208  }
209  if(overlap) continue;
210  Survivors.push_back(jetvec.at(itjet));
211  }
212  jetvec=Survivors;
213 
214  return;
215  }
216 
217  void LeptonJetOverlapRemoval(vector<HEPUtils::Particle*> &lepvec, vector<HEPUtils::Jet*> &jetvec, double DeltaRMax) {
218  //Routine to do lepton-jet check
219  //Discards leptons if they are within DeltaRMax of a jet
220 
221  vector<HEPUtils::Particle*> Survivors;
222 
223  for(unsigned int itlep = 0; itlep < lepvec.size(); itlep++) {
224  bool overlap = false;
225  HEPUtils::P4 lepmom=lepvec.at(itlep)->mom();
226  for(unsigned int itjet= 0; itjet < jetvec.size(); itjet++) {
227  HEPUtils::P4 jetmom=jetvec.at(itjet)->mom();
228  double dR;
229 
230  dR=jetmom.deltaR_eta(lepmom);
231 
232  if(fabs(dR) <= DeltaRMax) overlap=true;
233  }
234  if(overlap) continue;
235  Survivors.push_back(lepvec.at(itlep));
236  }
237  lepvec=Survivors;
238 
239  return;
240  }
241 
242 
243  public:
244 
245  // Required detector sim
246  static constexpr const char* detector = "ATLAS";
247 
248  Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb() {
249 
250  set_analysis_name("ATLAS_13TeV_RJ3L_lowmass_36invfb");
251  set_luminosity(36.);
252 
253  _num2L2JHIGH=0;
254  _num2L2JINT=0;
255  _num2L2JLOW=0;
256  _num2L2JCOMP=0;
257  _num3LHIGH=0;
258  _num3LINT=0;
259  _num3LLOW=0;
260  _num3LCOMP=0;
261 
262  NCUTS=70;
263 
264  for(int i=0;i<NCUTS;i++){
265  cutFlowVector.push_back(0);
266  cutFlowVector_str.push_back("");
267  }
268 
269 
270  // Recursive jigsaw stuff
271  #pragma omp critical (init_ATLAS_13TeV_RJ3L_lowmass_36invfb)
272  {
273 
274  // // DEBUG:
275  // RestFrames::SetLogPrint(RestFrames::LogDebug, true);
276  // RestFrames::SetLogPrint(RestFrames::LogVerbose, true);
277 
278  LAB_2L2J.reset(new RestFrames::LabRecoFrame("LAB_2L2J","lab2L2J"));
279  C1N2_2L2J.reset(new RestFrames::DecayRecoFrame("C1N2_2L2J","#tilde{#chi}^{ #pm}_{1} #tilde{#chi}^{ 0}_{2}"));
280  C1a_2L2J.reset(new RestFrames::DecayRecoFrame("C1a_2L2J","#tilde{#chi}^{ #pm}_{1}"));
281  N2b_2L2J.reset(new RestFrames::DecayRecoFrame("N2b_2L2J","#tilde{#chi}^{ 0}_{2}"));
282 
283  Wa_2L2J.reset(new RestFrames::DecayRecoFrame("Wa_2L2J","W_{a}"));
284  Zb_2L2J.reset(new RestFrames::DecayRecoFrame("Zb_2L2J","Z_{b}"));
285 
286  J1_2L2J.reset(new RestFrames::VisibleRecoFrame("J1_2L2J","#it{j}_{1}"));
287  J2_2L2J.reset(new RestFrames::VisibleRecoFrame("J2_2L2J","#it{j}_{2}"));
288  L1_2L2J.reset(new RestFrames::VisibleRecoFrame("L1_2L2J","#it{l}_{1}"));
289  L2_2L2J.reset(new RestFrames::VisibleRecoFrame("L2_2L2J","#it{l}_{2}"));
290 
291  X1a_2L2J.reset(new RestFrames::InvisibleRecoFrame("X1a_2L2J","#tilde{#chi}^{ 0}_{1 a}"));
292  X1b_2L2J.reset(new RestFrames::InvisibleRecoFrame("X1b_2L2J","#tilde{#chi}^{ 0}_{1 b}"));
293 
294  LAB_2L2J->SetChildFrame(*C1N2_2L2J);
295 
296  C1N2_2L2J->AddChildFrame(*C1a_2L2J);
297  C1N2_2L2J->AddChildFrame(*N2b_2L2J);
298 
299  C1a_2L2J->AddChildFrame(*Wa_2L2J);
300  C1a_2L2J->AddChildFrame(*X1a_2L2J);
301 
302  N2b_2L2J->AddChildFrame(*Zb_2L2J);
303  N2b_2L2J->AddChildFrame(*X1b_2L2J);
304 
305  Wa_2L2J->AddChildFrame(*J1_2L2J);
306  Wa_2L2J->AddChildFrame(*J2_2L2J);
307 
308  Zb_2L2J->AddChildFrame(*L1_2L2J);
309  Zb_2L2J->AddChildFrame(*L2_2L2J);
310 
311 
312  if(!LAB_2L2J->InitializeTree())
313  {
314  str errmsg;
315  errmsg = "Some problem occurred when calling LAB_2L2J->InitializeTree() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
316  piped_errors.request(LOCAL_INFO, errmsg);
317  }
318 
319 
321  //Setting the invisible
323  INV_2L2J.reset(new RestFrames::InvisibleGroup("INV_2L2J","#tilde{#chi}_{1}^{ 0} Jigsaws"));
324  INV_2L2J->AddFrame(*X1a_2L2J);
325  INV_2L2J->AddFrame(*X1b_2L2J);
326 
327  // Set di-LSP mass to minimum Lorentz-invariant expression
328  X1_mass_2L2J.reset(new RestFrames::SetMassInvJigsaw("X1_mass_2L2J", "Set M_{#tilde{#chi}_{1}^{ 0} #tilde{#chi}_{1}^{ 0}} to minimum"));
329  INV_2L2J->AddJigsaw(*X1_mass_2L2J);
330 
331  // Set di-LSP rapidity to that of visible particles
332  X1_eta_2L2J.reset(new RestFrames::SetRapidityInvJigsaw("X1_eta_2L2J", "#eta_{#tilde{#chi}_{1}^{ 0} #tilde{#chi}_{1}^{ 0}} = #eta_{2jet+2#it{l}}"));
333  INV_2L2J->AddJigsaw(*X1_eta_2L2J);
334  X1_eta_2L2J->AddVisibleFrames(C1N2_2L2J->GetListVisibleFrames());
335 
336 
337  X1X1_contra_2L2J.reset(new RestFrames::ContraBoostInvJigsaw("X1X1_contra_2L2J","Contraboost invariant Jigsaw"));
338  INV_2L2J->AddJigsaw(*X1X1_contra_2L2J);
339  X1X1_contra_2L2J->AddVisibleFrames(C1a_2L2J->GetListVisibleFrames(), 0);
340  X1X1_contra_2L2J->AddVisibleFrames(N2b_2L2J->GetListVisibleFrames(), 1);
341  X1X1_contra_2L2J->AddInvisibleFrame(*X1a_2L2J, 0);
342  X1X1_contra_2L2J->AddInvisibleFrame(*X1b_2L2J, 1);
343 
344  if(!LAB_2L2J->InitializeAnalysis())
345  {
346  str errmsg;
347  errmsg = "Some problem occurred when calling LAB_2L2J->InitializeAnalysis() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
348  piped_errors.request(LOCAL_INFO, errmsg);
349  }
350 
351 
352  LAB_3L.reset(new RestFrames::LabRecoFrame("LAB_3L","lab"));
353  C1N2_3L.reset(new RestFrames::DecayRecoFrame("C1N2_3L","#tilde{#chi}^{ #pm}_{1} #tilde{#chi}^{ 0}_{2}"));
354  C1a_3L.reset(new RestFrames::DecayRecoFrame("C1a_3L","#tilde{#chi}^{ #pm}_{1}"));
355  N2b_3L.reset(new RestFrames::DecayRecoFrame("N2b_3L","#tilde{#chi}^{ 0}_{2}"));
356 
357  L1a_3L.reset(new RestFrames::VisibleRecoFrame("L1a_3L","#it{l}_{1a}"));
358  L1b_3L.reset(new RestFrames::VisibleRecoFrame("L1b_3L","#it{l}_{1b}"));
359  L2b_3L.reset(new RestFrames::VisibleRecoFrame("L2b_3L","#it{l}_{2b}"));
360 
361  X1a_3L.reset(new RestFrames::InvisibleRecoFrame("X1a_3L","#tilde{#chi}^{ 0}_{1 a} + #nu_{a}"));
362  X1b_3L.reset(new RestFrames::InvisibleRecoFrame("X1b_3L","#tilde{#chi}^{ 0}_{1 b}"));
363 
364 
365  LAB_3L->SetChildFrame(*C1N2_3L);
366 
367  C1N2_3L->AddChildFrame(*C1a_3L);
368  C1N2_3L->AddChildFrame(*N2b_3L);
369 
370  C1a_3L->AddChildFrame(*L1a_3L);
371  C1a_3L->AddChildFrame(*X1a_3L);
372 
373  N2b_3L->AddChildFrame(*L1b_3L);
374  N2b_3L->AddChildFrame(*L2b_3L);
375  N2b_3L->AddChildFrame(*X1b_3L);
376 
377 
378  if(!LAB_3L->InitializeTree())
379  {
380  str errmsg;
381  errmsg = "Some problem occurred when calling LAB_3L->InitializeTree() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
382  piped_errors.request(LOCAL_INFO, errmsg);
383  }
384 
385 
386  //setting the invisible components
387  INV_3L.reset(new RestFrames::InvisibleGroup("INV_3L","Invisible system LSP mass Jigsaw"));
388  INV_3L->AddFrame(*X1a_3L);
389  INV_3L->AddFrame(*X1b_3L);
390 
391 
392  // Set di-LSP mass to minimum Lorentz-invariant expression
393  X1_mass_3L.reset(new RestFrames::SetMassInvJigsaw("X1_mass_3L", "Set M_{#tilde{#chi}_{1}^{ 0} #tilde{#chi}_{1}^{ 0}} to minimum"));
394  INV_3L->AddJigsaw(*X1_mass_3L);
395 
396  // Set di-LSP rapidity to that of visible particles and neutrino
397  X1_eta_3L.reset(new RestFrames::SetRapidityInvJigsaw("X1_eta_3L", "#eta_{#tilde{#chi}_{1}^{ 0} #tilde{#chi}_{1}^{ 0}} = #eta_{3#it{l}}"));
398  INV_3L->AddJigsaw(*X1_eta_3L);
399  X1_eta_3L->AddVisibleFrames(C1N2_3L->GetListVisibleFrames());
400 
401 
402  X1X1_contra_3L.reset(new RestFrames::ContraBoostInvJigsaw("X1X1_contra_3L","Contraboost invariant Jigsaw"));
403  INV_3L->AddJigsaw(*X1X1_contra_3L);
404  X1X1_contra_3L->AddVisibleFrames(C1a_3L->GetListVisibleFrames(),0);
405  X1X1_contra_3L->AddVisibleFrames(N2b_3L->GetListVisibleFrames(),1);
406  X1X1_contra_3L->AddInvisibleFrames(C1a_3L->GetListInvisibleFrames(),0);
407  X1X1_contra_3L->AddInvisibleFrames(N2b_3L->GetListInvisibleFrames(),1);
408 
409  if(!LAB_3L->InitializeAnalysis())
410  {
411  str errmsg;
412  errmsg = "Some problem occurred when calling LAB_3L->InitializeAnalysis() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
413  piped_errors.request(LOCAL_INFO, errmsg);
414  }
415 
416 
418  // RestFrames stuff
419 
420  // combinatoric (transverse) tree
421  // for jet assignment
422  LAB_comb.reset(new RestFrames::LabRecoFrame("LAB_comb","LAB"));
423  CM_comb.reset(new RestFrames::DecayRecoFrame("CM_comb","CM"));
424  S_comb.reset(new RestFrames::DecayRecoFrame("S_comb","S"));
425  ISR_comb.reset(new RestFrames::VisibleRecoFrame("ISR_comb","ISR"));
426  J_comb.reset(new RestFrames::VisibleRecoFrame("J_comb","Jets"));
427  L_comb.reset(new RestFrames::VisibleRecoFrame("L_comb","#it{l}'s"));
428  I_comb.reset(new RestFrames::InvisibleRecoFrame("I_comb","Inv"));
429 
430  LAB_comb->SetChildFrame(*CM_comb);
431  CM_comb->AddChildFrame(*ISR_comb);
432  CM_comb->AddChildFrame(*S_comb);
433  S_comb->AddChildFrame(*L_comb);
434  S_comb->AddChildFrame(*J_comb);
435  S_comb->AddChildFrame(*I_comb);
436 
437  if(!LAB_comb->InitializeTree())
438  {
439  str errmsg;
440  errmsg = "Some problem occurred when calling LAB_comb->InitializeTree() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
441  piped_errors.request(LOCAL_INFO, errmsg);
442  }
443 
444  // 2L+NJ tree (Z->ll + W/Z->qq)
445  LAB_2LNJ.reset(new RestFrames::LabRecoFrame("LAB_2LNJ","LAB"));
446  CM_2LNJ.reset(new RestFrames::DecayRecoFrame("CM_2LNJ","CM"));
447  S_2LNJ.reset(new RestFrames::DecayRecoFrame("S_2LNJ","S"));
448  ISR_2LNJ.reset(new RestFrames::VisibleRecoFrame("ISR_2LNJ","ISR"));
449  Ca_2LNJ.reset(new RestFrames::DecayRecoFrame("Ca_2LNJ","C_{a}"));
450  Z_2LNJ.reset(new RestFrames::DecayRecoFrame("Z_2LNJ","Z"));
451  L1_2LNJ.reset(new RestFrames::VisibleRecoFrame("L1_2LNJ","#it{l}_{1}"));
452  L2_2LNJ.reset(new RestFrames::VisibleRecoFrame("L2_2LNJ","#it{l}_{2}"));
453  Cb_2LNJ.reset(new RestFrames::DecayRecoFrame("Cb_2LNJ","C_{b}"));
454  JSA_2LNJ.reset(new RestFrames::SelfAssemblingRecoFrame("JSA_2LNJ", "J"));
455  J_2LNJ.reset(new RestFrames::VisibleRecoFrame("J_2LNJ","Jets"));
456  Ia_2LNJ.reset(new RestFrames::InvisibleRecoFrame("Ia_2LNJ","I_{a}"));
457  Ib_2LNJ.reset(new RestFrames::InvisibleRecoFrame("Ib_2LNJ","I_{b}"));
458 
459  LAB_2LNJ->SetChildFrame(*CM_2LNJ);
460  CM_2LNJ->AddChildFrame(*ISR_2LNJ);
461  CM_2LNJ->AddChildFrame(*S_2LNJ);
462  S_2LNJ->AddChildFrame(*Ca_2LNJ);
463  S_2LNJ->AddChildFrame(*Cb_2LNJ);
464  Ca_2LNJ->AddChildFrame(*Z_2LNJ);
465  Ca_2LNJ->AddChildFrame(*Ia_2LNJ);
466  Cb_2LNJ->AddChildFrame(*JSA_2LNJ);
467  Cb_2LNJ->AddChildFrame(*Ib_2LNJ);
468  Z_2LNJ->AddChildFrame(*L1_2LNJ);
469  Z_2LNJ->AddChildFrame(*L2_2LNJ);
470  JSA_2LNJ->AddChildFrame(*J_2LNJ);
471 
472  if(!LAB_2LNJ->InitializeTree())
473  {
474  str errmsg;
475  errmsg = "Some problem occurred when calling LAB_2LNJ->InitializeTree() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
476  piped_errors.request(LOCAL_INFO, errmsg);
477  }
478 
479 
480  // 2L+1L tree (Z->ll + Z/W->l)
481  LAB_2L1L.reset(new RestFrames::LabRecoFrame("LAB_2L1L","LAB"));
482  CM_2L1L.reset(new RestFrames::DecayRecoFrame("CM_2L1L","CM"));
483  S_2L1L.reset(new RestFrames::DecayRecoFrame("S_2L1L","S"));
484  ISR_2L1L.reset(new RestFrames::VisibleRecoFrame("ISR_2L1L","ISR"));
485  Ca_2L1L.reset(new RestFrames::DecayRecoFrame("Ca_2L1L","C_{a}"));
486  Z_2L1L.reset(new RestFrames::DecayRecoFrame("Z_2L1L","Z"));
487  L1_2L1L.reset(new RestFrames::VisibleRecoFrame("L1_2L1L","#it{l}_{1}"));
488  L2_2L1L.reset(new RestFrames::VisibleRecoFrame("L2_2L1L","#it{l}_{2}"));
489  Cb_2L1L.reset(new RestFrames::DecayRecoFrame("Cb_2L1L","C_{b}"));
490  Lb_2L1L.reset(new RestFrames::VisibleRecoFrame("Lb_2L1L","#it{l}_{b}"));
491  Ia_2L1L.reset(new RestFrames::InvisibleRecoFrame("Ia_2L1L","I_{a}"));
492  Ib_2L1L.reset(new RestFrames::InvisibleRecoFrame("Ia_2L1L","I_{b}"));
493 
494  LAB_2L1L->SetChildFrame(*CM_2L1L);
495  CM_2L1L->AddChildFrame(*ISR_2L1L);
496  CM_2L1L->AddChildFrame(*S_2L1L);
497  S_2L1L->AddChildFrame(*Ca_2L1L);
498  S_2L1L->AddChildFrame(*Cb_2L1L);
499  Ca_2L1L->AddChildFrame(*Z_2L1L);
500  Ca_2L1L->AddChildFrame(*Ia_2L1L);
501  Z_2L1L->AddChildFrame(*L1_2L1L);
502  Z_2L1L->AddChildFrame(*L2_2L1L);
503  Cb_2L1L->AddChildFrame(*Lb_2L1L);
504  Cb_2L1L->AddChildFrame(*Ib_2L1L);
505 
506  if(!LAB_2L1L->InitializeTree())
507  {
508  str errmsg;
509  errmsg = "Some problem occurred when calling LAB_2L1L->InitializeTree() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
510  piped_errors.request(LOCAL_INFO, errmsg);
511  }
512 
513 
515 
516  // combinatoric (transverse) tree
517  // for jet assignment
518  INV_comb.reset(new RestFrames::InvisibleGroup("INV_comb","Invisible System"));
519  INV_comb->AddFrame(*I_comb);
520 
521  InvMass_comb.reset(new RestFrames::SetMassInvJigsaw("InvMass_comb", "Invisible system mass Jigsaw"));
522  INV_comb->AddJigsaw(*InvMass_comb);
523 
524  JETS_comb.reset(new RestFrames::CombinatoricGroup("JETS_comb","Jets System"));
525  JETS_comb->AddFrame(*ISR_comb);
526  JETS_comb->SetNElementsForFrame(*ISR_comb, 1);
527  JETS_comb->AddFrame(*J_comb);
528  JETS_comb->SetNElementsForFrame(*J_comb, 0);
529 
530  SplitJETS_comb.reset(new RestFrames::MinMassesCombJigsaw("SplitJETS_comb", "Minimize M_{ISR} and M_{S} Jigsaw"));
531  JETS_comb->AddJigsaw(*SplitJETS_comb);
532  SplitJETS_comb->AddCombFrame(*ISR_comb, 0);
533  SplitJETS_comb->AddCombFrame(*J_comb, 1);
534  SplitJETS_comb->AddObjectFrame(*ISR_comb,0);
535  SplitJETS_comb->AddObjectFrame(*S_comb,1);
536 
537  if(!LAB_comb->InitializeAnalysis())
538  {
539  str errmsg;
540  errmsg = "Some problem occurred when calling LAB_comb->InitializeAnalysis() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
541  piped_errors.request(LOCAL_INFO, errmsg);
542  }
543 
544  // 2L+NJ tree (Z->ll + W/Z->qq)
545  INV_2LNJ.reset(new RestFrames::InvisibleGroup("INV_2LNJ","Invisible System"));
546  INV_2LNJ->AddFrame(*Ia_2LNJ);
547  INV_2LNJ->AddFrame(*Ib_2LNJ);
548 
549  InvMass_2LNJ.reset(new RestFrames::SetMassInvJigsaw("InvMass_2LNJ", "Invisible system mass Jigsaw"));
550  INV_2LNJ->AddJigsaw(*InvMass_2LNJ);
551  InvRapidity_2LNJ.reset(new RestFrames::SetRapidityInvJigsaw("InvRapidity_2LNJ", "Set inv. system rapidity"));
552  INV_2LNJ->AddJigsaw(*InvRapidity_2LNJ);
553  InvRapidity_2LNJ->AddVisibleFrames(S_2LNJ->GetListVisibleFrames());
554  SplitINV_2LNJ.reset(new RestFrames::ContraBoostInvJigsaw("SplitINV_2LNJ", "INV -> I_{a}+ I_{b} jigsaw"));
555  INV_2LNJ->AddJigsaw(*SplitINV_2LNJ);
556  SplitINV_2LNJ->AddVisibleFrames(Ca_2LNJ->GetListVisibleFrames(), 0);
557  SplitINV_2LNJ->AddVisibleFrames(Cb_2LNJ->GetListVisibleFrames(), 1);
558  SplitINV_2LNJ->AddInvisibleFrame(*Ia_2LNJ, 0);
559  SplitINV_2LNJ->AddInvisibleFrame(*Ib_2LNJ, 1);
560 
561  JETS_2LNJ.reset(new RestFrames::CombinatoricGroup("JETS_comb","Jets System"));
562  JETS_2LNJ->AddFrame(*J_2LNJ);
563  JETS_2LNJ->SetNElementsForFrame(*J_2LNJ, 0);
564 
565  if(!LAB_2LNJ->InitializeAnalysis())
566  {
567  str errmsg;
568  errmsg = "Some problem occurred when calling LAB_2LNJ->InitializeAnalysis() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
569  piped_errors.request(LOCAL_INFO, errmsg);
570  }
571 
572  // 2L+1L tree (Z->ll + Z/W->l)
573  INV_2L1L.reset(new RestFrames::InvisibleGroup("INV_2L1L","Invisible System"));
574  INV_2L1L->AddFrame(*Ia_2L1L);
575  INV_2L1L->AddFrame(*Ib_2L1L);
576 
577  InvMass_2L1L.reset(new RestFrames::SetMassInvJigsaw("InvMass_2L1L", "Invisible system mass Jigsaw"));
578  INV_2L1L->AddJigsaw(*InvMass_2L1L);
579  InvRapidity_2L1L.reset(new RestFrames::SetRapidityInvJigsaw("InvRapidity_2L1L", "Set inv. system rapidity"));
580  INV_2L1L->AddJigsaw(*InvRapidity_2L1L);
581  InvRapidity_2L1L->AddVisibleFrames(S_2L1L->GetListVisibleFrames());
582  SplitINV_2L1L.reset(new RestFrames::ContraBoostInvJigsaw("SplitINV_2L1L", "INV -> I_{a}+ I_{b} jigsaw"));
583  INV_2L1L->AddJigsaw(*SplitINV_2L1L);
584  SplitINV_2L1L->AddVisibleFrames(Ca_2L1L->GetListVisibleFrames(), 0);
585  SplitINV_2L1L->AddVisibleFrames(Cb_2L1L->GetListVisibleFrames(), 1);
586  SplitINV_2L1L->AddInvisibleFrame(*Ia_2L1L, 0);
587  SplitINV_2L1L->AddInvisibleFrame(*Ib_2L1L, 1);
588 
589  if(!LAB_2L1L->InitializeAnalysis())
590  {
591  str errmsg;
592  errmsg = "Some problem occurred when calling LAB_2L1L->InitializeAnalysis() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
593  piped_errors.request(LOCAL_INFO, errmsg);
594  }
595 
596  }
597  }
598 
599 
600  void run(const HEPUtils::Event* event) {
601 
602  // Clear
603  LAB_2L2J->ClearEvent();
604  LAB_comb->ClearEvent();
605  LAB_2LNJ->ClearEvent();
606  LAB_2L1L->ClearEvent();
607  LAB_3L->ClearEvent();
608 
609  // Missing energy
610  HEPUtils::P4 ptot = event->missingmom();
611  TVector3 ETMiss;
612  ETMiss.SetXYZ(ptot.px(),ptot.py(),0.0);
613 
614  // Baseline lepton objects
615  vector<HEPUtils::Particle*> baselineElectrons, baselineMuons, baselineTaus;
616 
617  for (HEPUtils::Particle* electron : event->electrons()) {
618  if (electron->pT() > 10. && electron->abseta() < 2.47) baselineElectrons.push_back(electron);
619  }
620 
621  // Apply electron efficiency
622  ATLAS::applyElectronEff(baselineElectrons);
623 
624  for (HEPUtils::Particle* muon : event->muons()) {
625  if (muon->pT() > 10. && muon->abseta() < 2.4) baselineMuons.push_back(muon);
626  }
627 
628  // Apply muon efficiency
629  ATLAS::applyMuonEff(baselineMuons);
630 
631  // Photons
632  vector<HEPUtils::Particle*> signalPhotons;
633  for (HEPUtils::Particle* photon : event->photons()) {
634  signalPhotons.push_back(photon);
635  }
636 
637 
638  // No taus used in 13 TeV analysis?
639  //for (HEPUtils::Particle* tau : event->taus()) {
640  //if (tau->pT() > 10. && tau->abseta() < 2.47) baselineTaus.push_back(tau);
641  //}
642  //ATLAS::applyTauEfficiencyR1(baselineTaus);
643 
644 
645  // Jets
646  vector<HEPUtils::Jet*> bJets;
647  vector<HEPUtils::Jet*> nonBJets;
648  vector<HEPUtils::Jet*> trueBJets; //for debugging
649 
650  // Get b jets
652  const std::vector<double> a = {0,10.};
653  const std::vector<double> b = {0,10000.};
654  const std::vector<double> c = {0.77}; // set b-tag efficiency to 77%
655  HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
656  for (HEPUtils::Jet* jet : event->jets()) {
657  bool hasTag=has_tag(_eff2d, jet->eta(), jet->pT());
658  if (jet->pT() > 20. && fabs(jet->eta()) < 2.4) {
659  if(jet->btag() && hasTag && fabs(jet->eta()) < 2.4 && jet->pT() > 20.){
660  bJets.push_back(jet);
661  } else {
662  nonBJets.push_back(jet);
663  }
664  }
665  }
666 
667 
668  // Overlap removal
669  vector<HEPUtils::Particle*> signalElectrons;
670  vector<HEPUtils::Particle*> signalMuons;
671  vector<HEPUtils::Particle*> signalLeptons;
672  vector<HEPUtils::Particle*> electronsForVeto;
673  vector<HEPUtils::Particle*> muonsForVeto;
674 
675  vector<HEPUtils::Jet*> signalJets;
676  vector<HEPUtils::Jet*> signalBJets;
677  vector<HEPUtils::Jet*> signalNonBJets;
678 
679  // Overlap removal is the same as the 8 TeV analysis
680  JetLeptonOverlapRemoval(nonBJets,baselineElectrons,0.2);
681  LeptonJetOverlapRemoval(baselineElectrons,nonBJets,0.4);
682  LeptonJetOverlapRemoval(baselineElectrons,bJets,0.4);
683  LeptonJetOverlapRemoval(baselineMuons,nonBJets,0.4);
684  LeptonJetOverlapRemoval(baselineMuons,bJets,0.4);
685 
686 
687  // Also we have already sorted jets by their b tag properties, so reset the b tag variable for each jet to the right category
688  // i.e. this was previously 100% true for true b jets then the efficiency map was applied above
689  for (HEPUtils::Jet* jet : bJets) {
690  jet->set_btag(true);
691  signalJets.push_back(jet);
692  signalBJets.push_back(jet);
693  }
694 
695  for (HEPUtils::Jet* jet : nonBJets) {
696  if(jet->pT() > 20. && fabs(jet->eta()) < 2.4) {
697  jet->set_btag(false);
698  signalJets.push_back(jet);
699  signalNonBJets.push_back(jet);
700  }
701  }
702 
703  //Put signal jets in pT order
704  std::sort(signalJets.begin(), signalJets.end(), sortByPT_RJ3L);
705  std::sort(signalBJets.begin(), signalBJets.end(), sortByPT_RJ3L);
706  std::sort(signalNonBJets.begin(), signalNonBJets.end(), sortByPT_RJ3L);
707 
708  for (HEPUtils::Particle* electron : baselineElectrons) {
709  signalElectrons.push_back(electron);
710  signalLeptons.push_back(electron);
711  }
712 
713  for (HEPUtils::Particle* muon : baselineMuons) {
714  signalMuons.push_back(muon);
715  signalLeptons.push_back(muon);
716  }
717 
718  std::sort(signalLeptons.begin(), signalLeptons.end(), sortLepByPT_RJ3L);
719 
720  // We now have the signal electrons, muons, jets and b jets- move on to the analysis
721  bool m_is2Lep=false;
722  bool m_is2Lep2Jet=false;
723  bool m_is2L2JInt=false;
724 
725  bool m_is3Lep=false;
726  bool m_is3LInt=false;
727  // bool m_is3Lep2Jet=false;
728  // bool m_is3Lep3Jet=false;
729 
730  // bool m_is4Lep=false;
731  // bool m_is4Lep2Jet=false;
732  // bool m_is4Lep3Jet=false;
733 
734  bool m_foundSFOS=false;
735 
736  // double m_H2PP_visible = -999.;
737  // double m_H2PP_invisible = -999.;
738  // double m_IaPP = -999.;
739  // double m_IbPP = -999.;
740  // double m_IaPa = -999.;
741  // double m_IbPb = -999.;
742  // double m_IaLAB = -999;
743  // double m_IbLAB = -999;
744  // double m_H4PP_Lept1A = -999.;
745  // double m_H4PP_Lept1B = -999.;
746  // double m_H4PP_Lept2B = -999.;
747  // double m_mu = -999;
748  // double m_pileUp_weight = -999;
749 
750 
752  // int m_nBaselineLeptons = -999;
753  // int m_nSignalLeptons = -999;
754 
755  double m_lept1Pt = -999;
756  // double m_lept1Eta = -999;
757  // double m_lept1Phi =-999;
758  double m_lept1sign=-999;
759  // double m_lept1origin = -999;
760  // double m_lept1type = -999;
761 
762  double m_lept2Pt =-999;
763  // double m_lept2Eta=-999;
764  // double m_lept2Phi =-999;
765  double m_lept2sign =-999;
766  // double m_lept2origin = -999;
767  // double m_lept2type = -999;
768 
769  double m_lept3Pt =-999;
770  // double m_lept3Eta =-999;
771  // double m_lept3Phi =-999;
772  double m_lept3sign =-999;
773  // double m_lept3origin = -999;
774  // double m_lept3type = -999;
775 
776  // double m_lept4Pt =-999;
777  // double m_lept4Eta =-999;
778  // double m_lept4Phi =-999;
779  // double m_lept4sign =-999;
780  // double m_lept4origin = -999;
781  // double m_lept4type = -999;
782  // double m_Zlep1Pt = -999;
783  // double m_Zlep1Phi = -999;
784  // double m_Zlep1Eta = -999;
785  // double m_Zlep1No = -999;
786  // double m_Zlep1sign = -999;
787 
788  // double m_Zlep2Pt = -999;
789  // double m_Zlep2sign = -999;
790  // double m_Zlep2Phi = -999;
791  // double m_Zlep2Eta = -999;
792  // double m_Zlep2No = -999;
793 
794  // double m_WlepPt = -999;
795  // double m_WlepPhi = -999;
796  // double m_WlepEta = -999;
797  // double m_WlepNo = -999;
798  // double m_Wlepsign = -999;
799 
800  // VR setup
801  // double m_lept1Pt_VR = -999;
802  // double m_lept1Eta_VR = -999;
803  // double m_lept1Phi_VR = -999;
804  // double m_lept1sign_VR = -999;
805 
806  // double m_lept2Pt_VR = -999;
807  // double m_lept2Eta_VR = -999;
808  // double m_lept2Phi_VR = -999;
809  // double m_lept2sign_VR = -999;
810 
811  //Jet Variables
812  int m_nJets=0;
813  // int m_nBtagJets=0;
814 
815  double m_jet1Pt =-999;
816  // double m_jet1Eta =-999;
817  // double m_jet1Phi =-999;
818  // double m_jet1M=-999;
819  // double m_jet1origin=-999;
820  // double m_jet1type=-999;
821 
822  double m_jet2Pt=-999;
823  // double m_jet2Eta=-999;
824  // double m_jet2Phi=-999;
825  // double m_jet2M=-999;
826  // double m_jet2origin=-999;
827  // double m_jet2type=-999;
828 
829  // double m_jet3Pt=-999;
830  // double m_jet3Eta=-999;
831  // double m_jet3Phi=-999;
832  // double m_jet3M=-999;
833  // double m_jet3origin=-999;
834  // double m_jet3type=-999;
835 
836  // double m_jet4Pt=-999;
837  // double m_jet4Eta=-999;
838  // double m_jet4Phi=-999;
839  // double m_jet4M=-999;
840  // double m_jet4origin=-999;
841  // double m_jet4type=-999;
842 
843  //Di-Lepton System: Calculated for OS Pairs
844  double m_mll=-999;
845  // double m_mt2=-999;
846  // double m_dRll=-999;
847  // double m_ptll=-999;
848  // double m_Zeta=-999;
849 
850  //Tri-Lepton System:
851  // double m_mlll=-999;
852  // double m_ptlll=-999;
853  double m_mTW=-999;
854  // double m_mTW_alt = -999;
855  // double m_mll_alt = -999;
856  //Di-Jet system: Calculated for the Two Leading Jets
857  double m_mjj=-999;
858  // double m_dRjj=-999;
859  // double m_ptjj=-999;
860  // double m_mj2j3 = -999;
861  //calculation of overall jet mass
862  // double m_mJ=-999;
863  // double m_mjjW=-999;//closest to the W-boson mass
864 
865  //Cleaning Variable: If MET is in the same direction as the Jet
866  double m_minDphi=-999;
867  // Some lab frame angles and stuff
868  // double m_dphill = -999;
869  // double m_dphilep1MET = -999;
870  // double m_dphilep2MET = -999;
871  // double m_dphilep3MET = -999;
872  // double m_dphiJMET = -999;
873  // double m_dphilll = -999;
874  // double m_dphilllMET = -999;
875  // double m_dphillMET = -999;
876  // double m_dphijj = -999;
877  // double m_dphijet1MET = -999;
878  // double m_dphijet2MET = -999;
879  // double m_dphijjMET = -999;
880  // double m_dphil3MET = -999;
881  // double m_MET=-999;
882  // double m_MET_phi = -999;
883  // double m_METTST = -999;
884  // double m_METTST_phi = -999;
885  // double m_Meff=-999;
886  // double m_LT=-999;
887 
888  // double m_MDR=-999;
889  // double m_PP_VisShape=-999;
890  // double m_gaminvPP=-999;
891  // double m_MP=-999;
892 
893  // double m_mC1=-999;
894  // double m_mN2=-999;
895 
896  // double m_mTW_Pa=-999;
897  // double m_mTW_PP=-999;
898 
899  // double m_mTZ_Pb=-999;
900  // double m_mTZ_PP=-999;
901 
902  // 3L CA
903  // double m_min_mt = -999;
904  // double m_pt_lll = -999;
905  // double m_mTl3 = -999;
906  //##############################//
907  //# Recursive Jigsaw Variables #//
908  //##############################//
909 
910  //Scale Variables
911  double m_H2PP=-999;
912  // double m_HT2PP=-999;
913  // double m_H3PP=-999;
914  // double m_HT3PP=-999;
915  double m_H4PP=-999;
916  double m_HT4PP=-999;
917  double m_H5PP=-999;
918  double m_HT5PP=-999;
919  // double m_H6PP=-999;
920  // double m_HT6PP=-999;
921 
922  double m_H2Pa=-999;
923  double m_H2Pb=-999;
924  double m_minH2P=-999;
925  // double m_R_H2Pa_H2Pb=-999;
926  double m_H3Pa=-999;
927  double m_H3Pb=-999;
928  double m_minH3P=-999;
929  // double m_R_H3Pa_H3Pb=-999;
930  double m_R_minH2P_minH3P=-999;
931  // double m_minR_pT2i_HT3Pi=-999;
932  // double m_maxR_H1PPi_H2PPi=-999;
933 
934  //Anglular Variables
935  // double m_cosPP=-999;
936  // double m_cosPa=-999;
937  // double m_cosPb=-999;
938  double m_dphiVP=-999;
939  // double m_dphiPPV=-999;
940  // double m_dphiPC1=-999;
941  // double m_dphiPN2=-999;
942 
943  // double m_sangle=-999;
944  // double m_dangle=-999;
945 
946  //Ratio Variables
947  // double m_RPZ_HT4PP=-999;
948  double m_RPT_HT4PP=-999;
949  // double m_R_HT4PP_H4PP=-999;
950 
951  // double m_RPZ_HT5PP=-999;
952  double m_RPT_HT5PP=-999;
953  // double m_R_HT5PP_H5PP=-999;
954  // double m_W_PP = -999;
955  // double m_WZ_PP = -999;
956 
958  double m_PTCM=-999;
959  double m_PTISR=-999;
960  double m_PTI=-999;
961  double m_RISR=-999;
962  // double m_cosCM=-999;
963  // double m_cosS=-999;
964  // double m_MISR=-999;
965  // double m_dphiCMI=-999;
966  // double m_dphiSI=-999;
967  double m_dphiISRI=-999;
968  // double m_HN2S=-999;
969  // double m_R_Ib_Ia=-999;
970  // double m_H11S = -999.;
971  // double m_HN1Ca = -999.;
972  // double m_HN1Cb = -999.;
973  // double m_H11Ca = -999.;
974  // double m_H11Cb = -999.;
975  // double m_cosC = -999.;
976  // double m_Is_Z = -999.;
977  // double m_Is_OS = -999;
978  double m_MZ = -999.;
979  double m_MJ = -999.;
980  // double m_mTWComp =-999.;
981  // double m_cosZ = -999.;
982  // double m_cosJ = -999.;
983  int m_NjS = 0;
984  int m_NjISR = 0;
985  int m_NbS = 0;
986  int m_NbISR = 0;
987 
988  // double m_MZ_VR = -999;
989  // double m_MJ_VR = -999;
990  // double m_PTCM_VR = -999;
991  // double m_PTISR_VR = -999;
992  // double m_PTI_VR = -999;
993  // double m_RISR_VR = -999;
994  // double m_dphiISRI_VR = -999;
995  // int m_NjS_VR = 0;
996  // int m_NjISR_VR = 0;
997 
998 
999  // double m_H2PP_VR = -999;
1000  // double m_H5PP_VR = -999;
1001  // double m_HT5PP_VR = -999;
1002  // double m_RPT_HT5PP_VR = -999;
1003  // double m_dphiVP_VR = -999;
1004  // double m_R_minH2P_minH3P_VR=-999;
1005 
1006  // double m_DPhi_METW = -999;
1007  //compressed
1008  // double m_WmassOnZ = -999;
1009  // double m_WptOnZ = -999;
1010  // double m_DPhi_METZ = -999;
1011  // double m_NonWJet_pT = -999;
1012  // double m_DPhi_METJetLeading = -999;
1013  // double m_DR_WOnZ2Jet = -999;
1014  // double m_DPhi_METNonWJet = -999;
1015  // double m_DPhi_METWonZ = -999;
1016 
1017  // Testing for low mass 3L
1018  // double m_M_I = -999;
1019  // double m_p_z_I = -999;
1020  // double m_p_z_Ia = -999;
1021  // double m_p_z_Ib = -999;
1022  // double m_boostx = -999;
1023  // double m_boosty = -999;
1024  // double m_boostz = -999;
1025 
1026  // Classify events
1027 
1028  m_nJets = signalJets.size();
1029 
1030  //if(signalLeptons.size()==2)std::cout << "m_nJets " << m_nJets << " signalLeptons.size() " << signalLeptons.size() << " pt1 " << signalLeptons[0]->pT() << " pt2 " << signalLeptons[1]->pT() << std::endl;
1031 
1032  if (signalLeptons.size()==2) m_is2Lep = true;
1033  else if (signalLeptons.size()==3) {m_is3Lep = true; //cout << "3L here" << endl;
1034  }
1035  // else if (signalLeptons.size()==4) m_is4Lep = true;
1036  //else return;
1037 
1038  if(m_is2Lep && m_nJets>1 ) m_is2Lep2Jet = true;
1039  if(m_is2Lep && m_nJets>2 && m_nJets<5) m_is2L2JInt = true;
1040  if(m_is3Lep && m_nJets>0 ) m_is3LInt = true;
1041  // if(m_is3Lep && m_nJets>1) m_is3Lep2Jet = true; //
1042  // if(m_is3Lep && m_nJets>2) m_is3Lep3Jet = true; //
1043  // if(m_is4Lep && m_nJets>1) m_is4Lep2Jet = true; //
1044  // if(m_is4Lep && m_nJets>2) m_is4Lep3Jet = true; //
1045 
1046  if(signalLeptons.size()==3)m_is3Lep=true;
1047 
1048  TLorentzVector metLV;
1049  //TLorentzVector bigFatJet;
1050  metLV.SetPxPyPzE(ptot.px(),ptot.py(),0.,sqrt(ptot.px()*ptot.px()+ptot.py()*ptot.py()));
1051 
1052  //Put the Jets in a more useful form
1053  vector<TLorentzVector> myJets;
1054  for(unsigned int ijet=0; ijet<signalJets.size();ijet++)
1055  {
1056  TLorentzVector tmp;
1057  tmp.SetPtEtaPhiM(signalJets[ijet]->pT(),signalJets[ijet]->eta(),signalJets[ijet]->phi(),signalJets[ijet]->mass());
1058  myJets.push_back(tmp);
1059  }
1060 
1061 
1062  //Put the Leptons in a more useful form
1063  vector<pair<TLorentzVector,int> > myLeptons;
1064  //vector<lep> myLeptons;
1065  for(unsigned int ilep=0; ilep<signalLeptons.size(); ilep++)
1066  {
1067  pair<TLorentzVector,int> temp;
1068  TLorentzVector tlv_temp;
1069 
1070  tlv_temp.SetPtEtaPhiM(signalLeptons[ilep]->pT(),signalLeptons[ilep]->eta(),signalLeptons[ilep]->phi(),0.0);
1071  temp.first = tlv_temp;
1072  int lepton_charge=0;
1073  if(signalLeptons[ilep]->pid()<0)lepton_charge=-1;
1074  if(signalLeptons[ilep]->pid()>0)lepton_charge=1;
1075  temp.second = lepton_charge;
1076  //temp.third = lepton_origin->at(lep_signal_index[ilep]);
1077  //temp.fourth = lepton_type->at(lep_signal_index[ilep]);
1078  //temp = make_tuple(tlv_temp,lepton_charge->at(lep_signal_index[ilep]),lepton_origin->at(lep_signal_index_[ilep]),lepton_type->at(lepton_signal_index[ilep]));
1079  myLeptons.push_back(temp);
1080  }
1081 
1082  sort(myJets.begin(), myJets.end(), SortJets);
1083  sort(myLeptons.begin(), myLeptons.end(), SortLeptons);
1084 
1085  if(m_is2Lep2Jet)
1086  {
1087 
1088  // if(myLeptons[0].first.Pt()<25.0 || myLeptons[1].first.Pt()<25.0) return;
1089 
1090  //Setting the Standard Variables
1091  //Di-Lepton System:
1092  m_lept1Pt = myLeptons[0].first.Pt();
1093  // m_lept1Eta = myLeptons[0].first.Eta();
1094  // m_lept1Phi = myLeptons[0].first.Phi();
1095  m_lept1sign = myLeptons[0].second;
1096 
1097  m_lept2Pt = myLeptons[1].first.Pt();
1098  // m_lept2Eta = myLeptons[1].first.Eta();
1099  // m_lept2Phi = myLeptons[1].first.Phi();
1100  m_lept2sign = myLeptons[1].second;
1101 
1102  m_mll = (myLeptons[0].first+myLeptons[1].first).M();
1103  // m_ptll = (myLeptons[0].first+myLeptons[1].first).Pt();
1104  // m_dRll = myLeptons[0].first.DeltaR(myLeptons[1].first);
1105  // m_Zeta = fabs(myLeptons[0].first.Eta() - myLeptons[1].first.Eta());
1106 
1107  vector<TLorentzVector> vleptons;
1108  vleptons.push_back(myLeptons[0].first);
1109  vleptons.push_back(myLeptons[1].first);
1110  //m_mt2 = myTool.GetMt2(vleptons,metLV);
1111 
1112  //min{d#phi}
1113  double mindphi=100000;
1114  double dphi=0;
1115  TLorentzVector tempjet;
1116  for(unsigned int ijet=0; ijet<signalJets.size();ijet++)
1117  {
1118  tempjet.SetPtEtaPhiM(signalJets[ijet]->pT(),signalJets[ijet]->eta(),signalJets[ijet]->phi(),signalJets[ijet]->mass());
1119 
1120  dphi = fabs(metLV.DeltaPhi(tempjet));
1121 
1122  if(dphi<mindphi) mindphi=dphi;
1123  }
1124 
1125  m_minDphi = mindphi;//cleaning variable for missmeasured jets;
1126 
1127  //just use the two leading jets
1128  int indexJ1=0;
1129  int indexJ2=1;
1130 
1131  //Di-Jet System: Here we decide which jets to use as output. The leading and sub-leading jet pair, or the jet pair with invariant mass closest to the W-Mass
1132  //jet closest to the W-boson mass
1133  m_jet1Pt = myJets[indexJ1].Pt();
1134  // m_jet1Eta = myJets[indexJ1].Eta();
1135  // m_jet1Phi = myJets[indexJ1].Phi();
1136  // m_jet1M = myJets[indexJ1].M();
1137 
1138  m_jet2Pt = myJets[indexJ2].Pt();
1139  // m_jet2Eta = myJets[indexJ2].Eta();
1140  // m_jet2Phi = myJets[indexJ2].Phi();
1141  // m_jet2M = myJets[indexJ2].M();
1142 
1143  // if(m_nJets>2) {
1144  // // m_jet3Pt = myJets[2].Pt();
1145  // // m_jet3Eta = myJets[2].Eta();
1146  // // m_jet3Phi = myJets[2].Phi();
1147  // // m_jet3M = myJets[2].M();
1148  // // m_mj2j3 = (myJets[1] + myJets[2]).M();
1149  // if(m_nJets>3) {
1150  // // m_jet4Pt = myJets[3].Pt();
1151  // // m_jet4Eta = myJets[3].Eta();
1152  // // m_jet4Phi = myJets[3].Phi();
1153  // // m_jet4M = myJets[3].M();
1154  // }
1155  // }
1156 
1157  m_mjj = (myJets[indexJ1]+myJets[indexJ2]).M();
1158  // m_ptjj = (myJets[indexJ1]+myJets[indexJ2]).Pt();
1159  // m_dRjj = myJets[indexJ1].DeltaR(myJets[indexJ2]);
1160 
1161 
1163  //Variables for the conventional approach
1164  // m_DPhi_METW = fabs((myJets[indexJ1]+myJets[indexJ2]).DeltaPhi(metLV));
1165  //for the comrpessed tree
1166  // double min_dPhi = 1000;
1167  // int WindexJ1 = -999;
1168  // int WindexJ2 = 999;
1169  // for (int j0=0;j0<myJets.size();j0++) {
1170 
1171  // double my_min_dphi= fabs(myJets[j0].DeltaPhi(metLV+myLeptons[0].first+myLeptons[1].first));
1172  // if (min_dPhi>my_min_dphi) {
1173  // min_dPhi = my_min_dphi;
1174  // WindexJ1 = j0;
1175  // }
1176  // }
1177 
1178  // min_dPhi = 1000;
1179  // for (int j1=0;j1<myJets.size();j1++) {
1180  // double my_min_dphi= fabs(myJets[j1].DeltaPhi(metLV+myLeptons[0].first+myLeptons[1].first));
1181 
1182  // if (min_dPhi>my_min_dphi) {
1183  // if (j1!=WindexJ1) {
1184  // min_dPhi = my_min_dphi;
1185  // WindexJ2 = j1;
1186  // }
1187  // }
1188  // }
1189 
1190  // m_WmassOnZ=(myJets[WindexJ1]+myJets[WindexJ2]).M();
1191  // m_WptOnZ=(myJets[WindexJ1]+myJets[WindexJ2]).Pt();
1192  // m_DPhi_METZ=fabs((myLeptons[0].first+myLeptons[1].first).DeltaPhi(metLV));
1193  // TLorentzVector nonWjetsLV;
1194  // for(int kjet=0;kjet<myJets.size();kjet++) {
1195  // if(kjet!=WindexJ1 && kjet!=WindexJ2) {
1196  // nonWjetsLV+=myJets[kjet];
1197  // }
1198  // }
1199  // if(m_nJets>2) {
1200  // // m_NonWJet_pT=nonWjetsLV.Pt();
1201  // // m_DPhi_METJetLeading = fabs(myJets[0].DeltaPhi(metLV));
1202  // // m_DR_WOnZ2Jet = myJets[WindexJ1].DeltaR(myJets[WindexJ2]);
1203  // // m_DPhi_METNonWJet = fabs(nonWjetsLV.DeltaPhi(metLV));
1204  // // m_DPhi_METWonZ = fabs((myJets[WindexJ1]+myJets[WindexJ2]).DeltaPhi(metLV));
1205  // }
1207 
1208 
1209  L1_2L2J->SetLabFrameFourVector(myLeptons[0].first); // Set lepton 4-vectors
1210  L2_2L2J->SetLabFrameFourVector(myLeptons[1].first);
1211  J1_2L2J->SetLabFrameFourVector(myJets[indexJ1]); // Set jets 4-vectors
1212  J2_2L2J->SetLabFrameFourVector(myJets[indexJ2]);
1213  TVector3 MET = ETMiss; // Get the MET
1214  MET.SetZ(0.);
1215  INV_2L2J->SetLabFrameThreeVector(MET); // Set the MET in reco tree
1216  TLorentzVector lep1;
1217  TLorentzVector lep2;
1219  //Lotentz vectors have been set, now do the boosts
1221 
1222  // Analyze the event
1223  if(!LAB_2L2J->AnalyzeEvent())
1224  {
1225  str errmsg;
1226  errmsg = "Some problem occurred when calling LAB_2L2J->AnalyzeEvent() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
1227  piped_warnings.request(LOCAL_INFO, errmsg);
1228  return;
1229  }
1230 
1231 
1232  //cout << L1_2L2J->GetFourVector(*LAB_2L2J).Pt() << endl;
1233  if (L1_2L2J->GetFourVector(*LAB_2L2J).Pt() > L2_2L2J->GetFourVector(*LAB_2L2J).Pt()){
1234  // m_Zlep1Pt = L1_2L2J->GetFourVector(*LAB_2L2J).Pt();
1235  // m_Zlep1sign = myLeptons[0].second;
1236  // m_Zlep1No = 0;
1237  // m_Zlep2Pt = L2_2L2J->GetFourVector(*LAB_2L2J).Pt();
1238  // m_Zlep2sign = myLeptons[1].second;
1239  // m_Zlep2No = 1;
1240  lep1 = L1_2L2J->GetFourVector(*LAB_2L2J);
1241  lep2 = L2_2L2J->GetFourVector(*LAB_2L2J);
1242  }
1243  else {
1244  // m_Zlep1Pt = L2_2L2J->GetFourVector(*LAB_2L2J).Pt();
1245  // m_Zlep1sign = myLeptons[1].second;
1246  // m_Zlep1No = 1;
1247  // m_Zlep2Pt = L1_2L2J->GetFourVector(*LAB_2L2J).Pt();
1248  // m_Zlep2sign = myLeptons[0].second;
1249  // m_Zlep2No = 0;
1250 
1251  lep1 = L2_2L2J->GetFourVector(*LAB_2L2J);
1252  lep2 = L1_2L2J->GetFourVector(*LAB_2L2J);
1253  }
1254  // set the jet lab frame 4-vector
1255  TLorentzVector jet1 = J1_2L2J->GetFourVector(*LAB_2L2J);
1256  TLorentzVector jet2 = J2_2L2J->GetFourVector(*LAB_2L2J);
1257 
1258  // Some lab frame stuff
1259  // m_dphill = lep1.DeltaPhi(lep2);
1260  // m_dphilep1MET = fabs(lep1.DeltaPhi(metLV));
1261  // m_dphilep2MET = fabs(lep2.DeltaPhi(metLV));
1262  // m_dphillMET = fabs((lep1 + lep2).DeltaPhi(metLV));
1263  // m_dphijet1MET = fabs(jet1.DeltaPhi(metLV));
1264  // m_dphijet2MET = fabs(jet2.DeltaPhi(metLV));
1265  // m_dphijj = fabs(jet1.DeltaPhi(jet2));
1266  // m_dphijjMET = fabs((jet1 + jet2).DeltaPhi(metLV));
1267  //... then by setting the Variables
1268  TLorentzVector vP_V1aPP = J1_2L2J->GetFourVector(*C1N2_2L2J);
1269  TLorentzVector vP_V2aPP = J2_2L2J->GetFourVector(*C1N2_2L2J);
1270  TLorentzVector vP_V1bPP = L1_2L2J->GetFourVector(*C1N2_2L2J);
1271  TLorentzVector vP_V2bPP = L2_2L2J->GetFourVector(*C1N2_2L2J);
1272  TLorentzVector vP_IaPP = X1a_2L2J->GetFourVector(*C1N2_2L2J);
1273  TLorentzVector vP_IbPP = X1b_2L2J->GetFourVector(*C1N2_2L2J);
1274 
1275  TLorentzVector vP_V1aPa = J1_2L2J->GetFourVector(*C1a_2L2J);
1276  TLorentzVector vP_V2aPa = J2_2L2J->GetFourVector(*C1a_2L2J);
1277  TLorentzVector vP_IaPa = X1a_2L2J->GetFourVector(*C1a_2L2J);
1278  TLorentzVector vP_V1bPb = L1_2L2J->GetFourVector(*N2b_2L2J);
1279  TLorentzVector vP_V2bPb = L2_2L2J->GetFourVector(*N2b_2L2J);
1280  TLorentzVector vP_IbPb = X1b_2L2J->GetFourVector(*N2b_2L2J);
1281 
1282 
1283  //Variables w/ 4 objects
1284  //Four vector sum of all visible objets + four vector sum of inv objects
1285  m_H2PP = (vP_V1aPP + vP_V2aPP + vP_V1bPP + vP_V2bPP).P() + (vP_IaPP+vP_IbPP).P();//H(1,1)PP
1286  // m_HT2PP = (vP_V1aPP + vP_V2aPP + vP_V1bPP + vP_V2bPP).Pt() + (vP_IaPP+vP_IbPP).Pt();//HT(1,1)PP
1287  //Scalar sum of all visible objects + vector sum of invisible momenta
1288  m_H5PP = vP_V1aPP.P() + vP_V2aPP.P() + vP_V1bPP.P() + vP_V2bPP.P() + (vP_IaPP + vP_IbPP).P();//H(4,1)PP
1289  m_HT5PP = vP_V1aPP.Pt() + vP_V2aPP.Pt() + vP_V1bPP.Pt() + vP_V2bPP.Pt() + (vP_IaPP + vP_IbPP).Pt();//HT(4,1)PP
1290  //scalar sum of all objects
1291  // m_H6PP = vP_V1aPP.P() + vP_V2aPP.P() + vP_V1bPP.P() + vP_V2bPP.P() + vP_IaPP.P() + vP_IbPP.P();//H(4,2)PP
1292  // m_HT6PP = vP_V1aPP.Pt() + vP_V2aPP.Pt() + vP_V1bPP.Pt() + vP_V2bPP.Pt() + vP_IaPP.Pt() + vP_IbPP.Pt();
1293 
1294  m_H2Pa = (vP_V1aPa + vP_V2aPa).P() + vP_IaPa.P();
1295  m_H2Pb = (vP_V1bPb + vP_V2bPb).P() + vP_IbPb.P();
1296  m_H3Pa = vP_V1aPa.P() + vP_V2aPa.P() + vP_IaPa.P();
1297  m_H3Pb = vP_V1bPb.P() + vP_V2bPb.P() + vP_IbPb.P();
1298  m_minH2P = std::min(m_H2Pa,m_H2Pb);
1299  m_minH3P = std::min(m_H3Pa,m_H3Pb);
1300  // m_R_H2Pa_H2Pb = m_H2Pa/m_H2Pb;
1301  // m_R_H3Pa_H3Pb = m_H3Pa/m_H3Pb;
1302  m_R_minH2P_minH3P = m_minH2P/m_minH3P;
1303  // std::cout << " m_R_minH2P_minH3P " << m_R_minH2P_minH3P << " " << m_minH2P << " " << m_minH3P << std::endl;
1304  // double H3PTa = vP_V1aPa.Pt() + vP_V2aPa.Pt() + vP_IaPa.Pt();
1305 
1306  // m_minR_pT2i_HT3Pi = std::min(vP_V1aPa.Pt()/H3PTa,vP_V2aPa.Pt()/H3PTa);
1307 
1308 
1309  // m_R_HT5PP_H5PP = m_HT5PP/m_H5PP;
1310 
1311  // Invisible in the PP frame, Px frames and lab frame
1312  TLorentzVector vP_IaLAB = X1a_2L2J->GetFourVector(*LAB_2L2J);
1313  TLorentzVector vP_IbLAB = X1b_2L2J->GetFourVector(*LAB_2L2J);
1314  // m_IaLAB = vP_IaLAB.P();
1315  // m_IbLAB = vP_IbLAB.P();
1316  // m_IaPP = vP_IaPP.P();
1317  // m_IbPP = vP_IbPP.P();
1318  // m_IaPa = vP_IaPa.P();
1319  // m_IbPb = vP_IbPb.P();
1320 
1321  // double jetMetphiP = (vP_V1aPa+vP_V2aPa).DeltaPhi(vP_IaPa);
1322  // m_mTW_Pa = sqrt(2*(vP_V1aPa+vP_V2aPa).Pt()*vP_IaPa.Pt()*(1-cos(jetMetphiP)));
1323 
1324  // double jetMetphiPP = (vP_V1aPP+vP_V2aPP).DeltaPhi(vP_IaPP+vP_IbPP);
1325  // m_mTW_PP = sqrt(2*(vP_V1aPP+vP_V2aPP).Pt()*(vP_IaPP+vP_IbPP).Pt()*(1-cos(jetMetphiPP)));
1326 
1327  // double dilepMetphiP = (vP_V1bPb+vP_V2bPb).DeltaPhi(vP_IbPb);
1328  // m_mTZ_Pb = sqrt(2*(vP_V1bPb+vP_V2bPb).Pt()*vP_IbPb.Pt()*(1-cos(dilepMetphiP)));
1329 
1330  // double dilepMetphiPP = (vP_V1bPP+vP_V2bPP).DeltaPhi(vP_IaPP+vP_IbPP);
1331  // m_mTZ_PP = sqrt(2*(vP_V1bPP+vP_V2bPP).Pt()*(vP_IaPP+vP_IbPP).Pt()*(1-cos(dilepMetphiPP)));
1332 
1333 
1334  // double H1PPa = (vP_V1aPP + vP_V2aPP).P();
1335  // double H1PPb = (vP_V1bPP + vP_V2bPP).P();
1336  // double H2PPa = vP_V1aPP.P() + vP_V2aPP.P();
1337  // double H2PPb = vP_V1bPP.P() + vP_V2bPP.P();
1338  // m_maxR_H1PPi_H2PPi = std::max(H1PPa/H2PPa,H1PPb/H2PPb);
1339 
1340  // signal variables
1341  TLorentzVector vP_Va = C1a_2L2J->GetVisibleFourVector(*C1a_2L2J);
1342  TLorentzVector vP_Vb = N2b_2L2J->GetVisibleFourVector(*N2b_2L2J);
1343  // m_MP = (vP_Va.M2()-vP_Vb.M2())/(2.*(vP_Va.E()-vP_Vb.E()));
1344 
1345  // double P_P = C1a_2L2J->GetMomentum(*C1N2_2L2J);
1346 
1347  // double MPP = 2.*sqrt(P_P*P_P + m_MP*m_MP);
1348  TVector3 vP_PP = C1N2_2L2J->GetFourVector(*LAB_2L2J).Vect();
1349  double Pt_PP = vP_PP.Pt();
1350  // double Pz_PP = fabs(vP_PP.Pz());
1351  m_RPT_HT5PP = Pt_PP / (Pt_PP + m_HT5PP);
1352  // m_RPZ_HT5PP = Pz_PP / (Pz_PP + m_HT5PP);
1353 
1354  // m_PP_VisShape = C1N2_2L2J->GetVisibleShape();
1355 
1356  // m_gaminvPP = 2.*m_MP/MPP;
1357  // m_MDR = m_PP_VisShape*C1N2_2L2J->GetMass();
1358 
1359  // m_mC1 = C1a_2L2J->GetMass();
1360  // m_mN2 = N2b_2L2J->GetMass();
1361 
1362 
1363  //Angular properties of the sparticles system
1364  // m_cosPP = C1N2_2L2J->GetCosDecayAngle(); //decay angle of the PP system
1365  // m_cosPa = C1a_2L2J->GetCosDecayAngle(*X1a_2L2J);//decay angle of the C1a system
1366  // m_cosPb = N2b_2L2J->GetCosDecayAngle(*X1b_2L2J);//decay angle of the N2b system
1367 
1368  //difference in azimuthal angle between the total sum of visible ojects in the C1N2 frame
1369  // m_dphiPPV = C1N2_2L2J->GetDeltaPhiBoostVisible();
1370  m_dphiVP = C1N2_2L2J->GetDeltaPhiDecayVisible();
1371 
1372  //hemisphere variables
1373  // m_dphiPC1 = C1a_2L2J->GetDeltaPhiDecayPlanes(*Wa_2L2J);
1374  // m_dphiPN2 = N2b_2L2J->GetDeltaPhiDecayPlanes(*Zb_2L2J);
1375 
1376  // m_sangle =(m_cosPa+(m_dphiVP-acos(-1.)/2.)/(acos(-1.)/2.))/2.;
1377  // m_dangle =(m_cosPa-(m_dphiVP-acos(-1.)/2.)/(acos(-1.)/2.))/2.;
1378  }//end is 2L2J event
1379 
1380  if(m_is3Lep){
1381 
1382  // bool m_pass3L_presel;
1383 
1384  // if(myLeptons[0].first.Pt()<25.0 || myLeptons[1].first.Pt()<25.0 || myLeptons[2].first.Pt()<20.0) {
1385  // m_pass3L_presel=false;
1386  // }
1387  // else {
1388  // m_pass3L_presel=true;
1389  // }
1390 
1391  //if(!m_pass3L_presel)return;
1392 
1393  //Tri-Lepton System
1394  //Here we choose leptons based on where they "come from"
1395  //lept1 and lept2 are the lepton pair with invariant mass closest to the Z-Mass
1396  //lept3 is the remaining lepton
1397  //This is meant to emulate lept1 and lept2 being produced by the Z, while lept3 is produced by the W
1398 
1399  double diff = 10000000000.0;
1400  int Zlep1 = -99;
1401  int Zlep2 = -99;
1402  double Zmass = -999.0;
1403  bool foundSFOS = false;
1404 
1405  for(unsigned int i=0; i<myLeptons.size(); i++)
1406  {
1407  for(unsigned int j=i+1; j<myLeptons.size(); j++)
1408  {
1409  //Opposite-Sign
1410  if(myLeptons[i].second*myLeptons[j].second<0)
1411  {
1412  //Same-Flavor
1413  if(abs(myLeptons[i].second)==abs(myLeptons[j].second))
1414  {
1415  double mass = (myLeptons[i].first+myLeptons[j].first).M();
1416  double massdiff = fabs(mass-91.1876);
1417  if(massdiff<diff)
1418  {
1419  diff=massdiff;
1420  Zmass=mass;
1421  Zlep1 = i;
1422  Zlep2 = j;
1423  foundSFOS = true;
1424  }
1425  }
1426  }
1427  }
1428  }
1429 
1430  if(!foundSFOS) {
1431  m_foundSFOS=false;
1432  }
1433  else {
1434  m_foundSFOS=true;
1435  }
1436 
1437  if(m_foundSFOS){
1438 
1439  int Wlep1 = -999;
1440  if( (Zlep1==0 && Zlep2==1) || (Zlep1==1 && Zlep2==0) ) Wlep1=2;
1441  else if( (Zlep1==0 && Zlep2==2) || (Zlep1==2 && Zlep2==0) ) Wlep1=1;
1442  else if((Zlep1==1 && Zlep2==2) || (Zlep1==2 && Zlep2==1) ) Wlep1=0;
1443 
1444  //Knowing the indices, we perform assignments
1445  m_lept1Pt = myLeptons[0].first.Pt();
1446  m_lept1sign = myLeptons[0].second;
1447 
1448  m_lept2Pt = myLeptons[1].first.Pt();
1449  m_lept2sign = myLeptons[1].second;
1450 
1451  m_lept3Pt = myLeptons[2].first.Pt();
1452  m_lept3sign = myLeptons[2].second;
1453 
1454 
1455  m_mll = Zmass; //based on mass minimization
1456 
1457  vector<TLorentzVector> Leptons;
1458  Leptons.push_back(myLeptons[Wlep1].first);
1459  Leptons.push_back(myLeptons[Zlep1].first);
1460  Leptons.push_back(myLeptons[Zlep2].first);
1461 
1462 
1463  double wlepMetphi = myLeptons[Wlep1].first.DeltaPhi(metLV);
1464 
1465  m_mTW = sqrt(2*myLeptons[Wlep1].first.Pt()*metLV.Pt()*(1-cos(wlepMetphi)));
1466 
1467  INV_3L->SetLabFrameThreeVector(ETMiss); //set the MET in the event
1468 
1469 
1470  L1a_3L->SetLabFrameFourVector(Leptons[0]); // Set lepton from W
1471  L1b_3L->SetLabFrameFourVector(Leptons[1]); // Set lepton1 from Z
1472  L2b_3L->SetLabFrameFourVector(Leptons[2]); // Set lepton2 from Z
1473 
1474  if(!LAB_3L->AnalyzeEvent())
1475  {
1476  str errmsg;
1477  errmsg = "Some problem occurred when calling LAB_3L->AnalyzeEvent() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
1478  piped_warnings.request(LOCAL_INFO, errmsg);
1479  return;
1480  }
1481 
1482 
1483  TLorentzVector l1;
1484  TLorentzVector l2;
1485  TLorentzVector l3 = L1a_3L->GetFourVector(*LAB_3L);
1486 
1487  //if(DEBUG) cout << "WlepPt: " << m_WlepPt << " Wlepsign: " << m_Wlepsign << endl;
1488  if (L1b_3L->GetFourVector(*LAB_3L).Pt() > L2b_3L->GetFourVector(*LAB_3L).Pt()){
1489  l1 = L1b_3L->GetFourVector(*LAB_3L);
1490  l2 = L2b_3L->GetFourVector(*LAB_3L);
1491  }
1492  else {
1493  l2 = L1b_3L->GetFourVector(*LAB_3L);
1494  l1 = L2b_3L->GetFourVector(*LAB_3L);
1495  }
1496 
1497  // More lab frame stuff
1498 
1499  //if(DEBUG) cout << "Zlep1: " << m_Zlep1Pt << " " << m_Zlep1sign << " Zlep2Pt: " << m_Zlep2Pt << " " << m_Zlep2sign << endl;
1500  TLorentzVector vP_V1aPP = L1a_3L->GetFourVector(*C1N2_3L);
1501  TLorentzVector vP_V1bPP = L1b_3L->GetFourVector(*C1N2_3L);
1502  TLorentzVector vP_V2bPP = L2b_3L->GetFourVector(*C1N2_3L);
1503  TLorentzVector vP_I1aPP = X1a_3L->GetFourVector(*C1N2_3L);
1504  TLorentzVector vP_I1bPP = X1b_3L->GetFourVector(*C1N2_3L);
1505 
1506  TLorentzVector vP_V1aPa = L1a_3L->GetFourVector(*C1a_3L);
1507  TLorentzVector vP_I1aPa = X1a_3L->GetFourVector(*C1a_3L);
1508 
1509  TLorentzVector vP_V1bPb = L1b_3L->GetFourVector(*N2b_3L);
1510  TLorentzVector vP_V2bPb = L2b_3L->GetFourVector(*N2b_3L);
1511  TLorentzVector vP_I1bPb = X1b_3L->GetFourVector(*N2b_3L);
1512 
1513 
1514 
1515  //Variables w/ 4 objects
1516 
1518  //Four vector sum of all visible objets + four vector sum of inv objects
1519 
1520  //Scalar sum of all visible objects + vector sum of invisible momenta
1521  m_H4PP = vP_V1aPP.P() + vP_V1bPP.P() + vP_V2bPP.P() + (vP_I1aPP + vP_I1bPP).P();//H(3,1)PP
1522  m_HT4PP = vP_V1aPP.Pt() + vP_V1bPP.Pt() + vP_V2bPP.Pt() + (vP_I1aPP + vP_I1bPP).Pt();//HT(3,1)PP
1523 
1524  // Invisible components again
1525  TLorentzVector vP_IaLAB = X1a_3L->GetFourVector(*LAB_3L);
1526  TLorentzVector vP_IbLAB = X1b_3L->GetFourVector(*LAB_3L);
1527 
1528 
1529  // Testing for low mass 3L
1530  TLorentzVector p_Ia_Lab = X1a_3L->GetFourVector(*LAB_3L);
1531  TLorentzVector p_Ib_Lab = X1b_3L->GetFourVector(*LAB_3L);
1532  TVector3 lab_to_pp = C1N2_3L->GetBoostInParentFrame();
1533 
1535  m_H2Pa = (vP_V1aPa).P() + (vP_I1aPa).P(); //H(1,1)P
1536  m_H2Pb = (vP_V1bPb + vP_V2bPb).P() + vP_I1bPb.P();//H(1,1)P
1537 
1538  m_H3Pa = vP_V1aPa.P() + vP_I1aPa.P();//H(1,1)P
1539  m_H3Pb = vP_V1bPb.P() + vP_V2bPb.P() + vP_I1bPb.P();//H(2,1)P
1540 
1541  m_minH2P = std::min(m_H2Pa,m_H2Pb);
1542  m_minH3P = std::min(m_H3Pa,m_H3Pb);
1543  // m_R_H2Pa_H2Pb = m_H2Pa/m_H2Pb;
1544  // m_R_H3Pa_H3Pb = m_H3Pa/m_H3Pb;
1545  m_R_minH2P_minH3P = m_H2Pb/m_H3Pb;
1546 
1547  // double H1PPa = (vP_V1aPP).P();
1548  // double H1PPb = (vP_V1bPP + vP_V2bPP).P();
1549  // double H2PPa = vP_V1aPP.P() + vP_I1aPP.P();
1550  // double H2PPb = (vP_V1bPP+vP_V2bPP).P() + vP_I1bPP.P();
1551  // m_maxR_H1PPi_H2PPi = std::max(H1PPa/H2PPa,H1PPb/H2PPb);
1552 
1554  //m_dRll_I_PP = (vP_V1bPP+vP_V1bPP).DeltaR(vP_I1bPP);
1555  //m_R_Ib_Ia = (vP_V1bPP + vP_V2bPP + vP_I1bPP).P()/(vP_V1aPP+vP_I1aPP).P();
1556 
1557  // signal variables
1558  TLorentzVector vP_Va = C1a_3L->GetVisibleFourVector(*C1a_3L);
1559  TLorentzVector vP_Vb = N2b_3L->GetVisibleFourVector(*N2b_3L);
1560 
1561  TVector3 vP_PP = C1N2_3L->GetFourVector(*LAB_3L).Vect();
1562  double Pt_PP = vP_PP.Pt();
1563  m_RPT_HT4PP = Pt_PP / (Pt_PP + m_HT4PP);
1564 
1565 
1566  // mt_min here
1567 
1568  /*double min0 = -999;
1569  double min1 = -999;
1570  double min2 = -999;
1571  double lepmetphi0 = myLeptons[0].first.DeltaPhi(metLV);
1572  double lepmetphi1 = myLeptons[1].first.DeltaPhi(metLV);
1573  double lepmetphi2 = myLeptons[2].first.DeltaPhi(metLV);
1574 
1575  if (myLeptons[0].second == -myLeptons[1].second) min0 = sqrt(2*myLeptons[2].first.Pt()*metLV.Pt()*(1-cos(lepmetphi2)));
1576  if (myLeptons[0].second == -myLeptons[2].second) min1 = sqrt(2*myLeptons[1].first.Pt()*metLV.Pt()*(1-cos(lepmetphi1)));
1577  if (myLeptons[1].second == -myLeptons[2].second) min2 = sqrt(2*myLeptons[0].first.Pt()*metLV.Pt()*(1-cos(lepmetphi0)));
1578 
1579  if (min0 > 0 && min1 > 0) m_min_mt = min(min0,min1);
1580  else if (min0 > 0 && min2 > 0) m_min_mt = min(min0,min2);
1581  else if (min1 > 0 && min2 > 0) m_min_mt = min(min1,min2);
1582  else if (min0 > 0 && min1 < 0 && min2 < 0) m_min_mt = min0;
1583  else if (min1 > 0 && min0 < 0 && min2 < 0) m_min_mt = min1;
1584  else if (min2 > 0 && min0 < 0 && min1 <0) m_min_mt = min2;*/
1585 
1586  } // end of if(m_foundSFOS)
1587  } // end of m_is3Lep
1588 
1589  if(m_is3LInt || m_is2L2JInt) {
1590 
1591  //min{d#phi}
1592  double mindphi=100000;
1593  double dphi=0;
1594  TLorentzVector tempjet;
1595  for(unsigned int ijet=0; ijet<signalJets.size();ijet++)
1596  {
1597  tempjet.SetPtEtaPhiM(signalJets[ijet]->pT(),signalJets[ijet]->eta(),signalJets[ijet]->phi(),signalJets[ijet]->mass());
1598 
1599  dphi = fabs(metLV.DeltaPhi(tempjet));
1600 
1601  if(dphi<mindphi) mindphi=dphi;
1602  }
1603 
1604  m_minDphi = mindphi;//cleaning variable for missmeasured jets;
1605  //if( fabs(mindphi)<0.4) return;
1606 
1607 
1608 
1609  vector<RestFrames::RFKey> jetID;
1610  for(int i = 0; i < int(myJets.size()); i++){
1611 
1612  TLorentzVector jet = myJets[i];
1613 
1614  // transverse view of jet 4-vectors
1615  jet.SetPtEtaPhiM(jet.Pt(),0.0,jet.Phi(),jet.M());
1616  jetID.push_back(JETS_comb->AddLabFrameFourVector(jet));
1617  }
1618 
1619  TLorentzVector lepSys(0.,0.,0.,0.);
1620  for(int i = 0; i < int(myLeptons.size()); i++){
1621  TLorentzVector lep1;
1622  // transverse view of mu 4-vectors
1623  lep1.SetPtEtaPhiM(myLeptons[i].first.Pt(),0.0,myLeptons[i].first.Phi(),myLeptons[i].first.M());
1624  lepSys = lepSys + lep1;
1625  }
1626  L_comb->SetLabFrameFourVector(lepSys);
1627 
1628  INV_comb->SetLabFrameThreeVector(ETMiss);
1629 
1630  if(!LAB_comb->AnalyzeEvent())
1631  {
1632  str errmsg;
1633  errmsg = "Some problem occurred when calling LAB_comb->AnalyzeEvent() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
1634  piped_warnings.request(LOCAL_INFO, errmsg);
1635  return;
1636  }
1637 
1638  for(int i = 0; i < int(signalJets.size()); i++){
1639  if(JETS_comb->GetFrame(jetID[i]) == *J_comb){
1640  m_NjS++;
1641  if(signalJets[i]->btag()) m_NbS++;
1642  } else {
1643  m_NjISR++;
1644  if(signalJets[i]->btag()) m_NbISR++;
1645  }
1646  }
1647 
1648  // 2LNJ analysis
1649  if(m_is2L2JInt){
1650 
1651  // put jets in their place
1652  int NJ = jetID.size();
1653  TLorentzVector vISR(0.,0.,0.,0.);
1654  for(int i = 0; i < NJ; i++){
1655  if(JETS_comb->GetFrame(jetID[i]) == *J_comb){
1656  JETS_2LNJ->AddLabFrameFourVector(myJets[i]);
1657  } else {
1658  vISR += myJets[i];
1659  }
1660  }
1661 
1662  ISR_2LNJ->SetLabFrameFourVector(vISR);
1663 
1664  // put leptons in their place
1665  L1_2LNJ->SetLabFrameFourVector(myLeptons[0].first);
1666  L2_2LNJ->SetLabFrameFourVector(myLeptons[1].first);
1667 
1668  INV_2LNJ->SetLabFrameThreeVector(ETMiss);
1669 
1670  if(!LAB_2LNJ->AnalyzeEvent())
1671  {
1672  str errmsg;
1673  errmsg = "Some problem occurred when calling LAB_2LNJ->AnalyzeEvent() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
1674  piped_warnings.request(LOCAL_INFO, errmsg);
1675  return;
1676  }
1677 
1678  }
1679 
1680  // 2L1L analysis
1681  if(m_is3LInt){
1682 
1683  // put jets in their place
1684  int NJ = jetID.size();
1685  TLorentzVector vISR(0.,0.,0.,0.);
1686  for(int i = 0; i < NJ; i++){
1687  if(JETS_comb->GetFrame(jetID[i]) != *J_comb) vISR += myJets[i];
1688  }
1689 
1690  ISR_2L1L->SetLabFrameFourVector(vISR);
1691 
1692  // put leptons in their place
1693  // find min mass OS pair
1694  pair<int,int> iSFOS;
1695  double mSFOS = -1.;
1696  for(int i = 0; i < 2; i++){
1697  for(int j = i+1; j < 3; j++){
1698  if((signbit(myLeptons[i].second) && !signbit(myLeptons[j].second)) || (!signbit(myLeptons[i].second) && signbit(myLeptons[j].second))){
1699  if(mSFOS < 0. ||
1700  (myLeptons[i].first+myLeptons[j].first).M() < mSFOS){
1701  mSFOS = (myLeptons[i].first+myLeptons[j].first).M();
1702  iSFOS.first = i;
1703  iSFOS.second = j;
1704  }
1705  }
1706  }
1707  }
1708 
1709  for(int i = 0; i < 3; i++){
1710  if(i == iSFOS.first)
1711  L1_2L1L->SetLabFrameFourVector(myLeptons[i].first);
1712  if(i == iSFOS.second)
1713  L2_2L1L->SetLabFrameFourVector(myLeptons[i].first);
1714  if(i != iSFOS.first && i != iSFOS.second) {
1715  Lb_2L1L->SetLabFrameFourVector(myLeptons[i].first);
1716  //calculate the mTWComp with the remaining lepton
1717  TLorentzVector themetLV;
1718  themetLV.SetPxPyPzE(ETMiss.X(),ETMiss.Y(),0.,sqrt(ETMiss.X()*ETMiss.X()+ETMiss.Y()*ETMiss.Y()));
1719  // double wlepMetphi = myLeptons[i].first.DeltaPhi(themetLV);
1720  // m_mTWComp = sqrt(2*myLeptons[i].first.Pt()*themetLV.Pt()*(1-cos(wlepMetphi)));
1721  }
1722  }
1723 
1724  INV_2L1L->SetLabFrameThreeVector(ETMiss);
1725 
1726  if(!LAB_2L1L->AnalyzeEvent())
1727  {
1728  str errmsg;
1729  errmsg = "Some problem occurred when calling LAB_2L1L->AnalyzeEvent() from the Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb analysis class.\n";
1730  piped_warnings.request(LOCAL_INFO, errmsg);
1731  return;
1732  }
1733 
1734  }
1735 
1736  TLorentzVector vP_CM;
1737  TLorentzVector vP_ISR;
1738  TLorentzVector vP_I;
1739 
1740  if(m_is2L2JInt){
1741 
1742  vP_CM = CM_2LNJ->GetFourVector();
1743  vP_ISR = ISR_2LNJ->GetFourVector();
1744  vP_I = (*Ia_2LNJ+*Ib_2LNJ).GetFourVector();
1745 
1746  // m_cosCM = CM_2LNJ->GetCosDecayAngle();
1747  // m_cosS = S_2LNJ->GetCosDecayAngle();
1748  // m_MISR = ISR_2LNJ->GetMass();
1749  // m_dphiCMI = acos(-1.)-fabs(CM_2LNJ->GetDeltaPhiBoostVisible());
1750  // m_dphiSI = acos(-1.)-fabs(S_2LNJ->GetDeltaPhiBoostVisible());
1751 
1752  // m_HN2S = //Z_2LNJ->GetFourVector(*S_2LNJ).E() +
1753  // L1_2LNJ->GetFourVector(*S_2LNJ).E()+
1754  // L2_2LNJ->GetFourVector(*S_2LNJ).E()+
1755  // J_2LNJ->GetFourVector(*S_2LNJ).E() +
1756  // Ia_2LNJ->GetFourVector(*S_2LNJ).P() +
1757  // Ib_2LNJ->GetFourVector(*S_2LNJ).P();
1758  // m_H11S = 2.*(*Ia_2LNJ+*Ib_2LNJ).GetFourVector(*S_2LNJ).P();
1759  // m_HN1Ca = Z_2LNJ->GetFourVector(*Ca_2LNJ).E()+
1760  Ia_2LNJ->GetFourVector(*Ca_2LNJ).P();
1761  // m_HN1Cb = J_2LNJ->GetFourVector(*Cb_2LNJ).E()+
1762  // Ib_2LNJ->GetFourVector(*Cb_2LNJ).P();
1763  // m_H11Ca = 2.*Ia_2LNJ->GetFourVector(*Ca_2LNJ).P();
1764  // m_H11Cb = 2.*Ib_2LNJ->GetFourVector(*Cb_2LNJ).P();
1765  // m_cosC = Ca_2LNJ->GetCosDecayAngle();
1766 
1767  // if((signbit(myLeptons[0].second) && !signbit(myLeptons[1].second)) || (!signbit(myLeptons[0].second) && signbit(myLeptons[1].second))) m_Is_OS = 1;
1768  // if(myLeptons[0].second+myLeptons[1].second == 0) m_Is_Z = 1;
1769  m_MZ = Z_2LNJ->GetMass();
1770  m_MJ = J_2LNJ->GetMass();
1771 
1772  // m_cosZ = Z_2LNJ->GetCosDecayAngle();
1773  //if(m_NjS > 1)
1774  // m_cosJ = JSA_2LNJ->GetCosDecayAngle();
1775  // m_dphiJMET = fabs(J_2LNJ->GetFourVector(*LAB_2LNJ).DeltaPhi(metLV));
1776  }
1777 
1778  if(m_is3LInt){
1779 
1780  vP_CM = CM_2L1L->GetFourVector();
1781  vP_ISR = ISR_2L1L->GetFourVector();
1782  vP_I = (*Ia_2L1L+*Ib_2L1L).GetFourVector();
1783 
1784  // m_cosCM = CM_2L1L->GetCosDecayAngle();
1785  // m_cosS = S_2L1L->GetCosDecayAngle();
1786  // m_MISR = ISR_2L1L->GetMass();
1787  // m_dphiCMI = acos(-1.)-fabs(CM_2L1L->GetDeltaPhiBoostVisible());
1788  // m_dphiSI = acos(-1.)-fabs(S_2L1L->GetDeltaPhiBoostVisible());
1789 
1790  // m_HN2S = //Z_2L1L->GetFourVector(*S_2L1L).E() +
1791  // L1_2L1L->GetFourVector(*S_2L1L).E() +
1792  // L2_2L1L->GetFourVector(*S_2L1L).E() +
1793  // Lb_2L1L->GetFourVector(*S_2L1L).E() +
1794  // Ia_2L1L->GetFourVector(*S_2L1L).P() +
1795  // Ib_2L1L->GetFourVector(*S_2L1L).P();
1796  // m_H11S = 2.*(*Ia_2L1L+*Ib_2L1L).GetFourVector(*S_2L1L).P();
1797  // m_HN1Ca = Z_2L1L->GetFourVector(*Ca_2L1L).E()+
1798  Ia_2L1L->GetFourVector(*Ca_2L1L).P();
1799  // m_HN1Cb = Lb_2L1L->GetFourVector(*Cb_2L1L).E()+
1800  // Ib_2L1L->GetFourVector(*Cb_2L1L).P();
1801  // m_H11Ca = 2.*Ia_2L1L->GetFourVector(*Ca_2L1L).P();
1802  // m_H11Cb = 2.*Ib_2L1L->GetFourVector(*Cb_2L1L).P();
1803  // m_cosC = Ca_2L1L->GetCosDecayAngle();
1804  // m_Is_OS = 1;
1805  // if(myLeptons[0].second+myLeptons[1].second == 0 ||
1806  // myLeptons[0].second+myLeptons[2].second == 0 ||
1807  // myLeptons[1].second+myLeptons[2].second == 0) m_Is_Z=1;
1808  m_MZ = Z_2L1L->GetMass();
1809  // m_cosZ = Z_2L1L->GetCosDecayAngle();
1810  }
1811 
1812  m_PTCM = vP_CM.Pt();
1813 
1814  TVector3 boostZ = vP_CM.BoostVector();
1815  boostZ.SetX(0.);
1816  boostZ.SetY(0.);
1817 
1818  vP_CM.Boost(-boostZ);
1819  vP_ISR.Boost(-boostZ);
1820  vP_I.Boost(-boostZ);
1821 
1822  TVector3 boostT = vP_CM.BoostVector();
1823  vP_ISR.Boost(-boostT);
1824  vP_I.Boost(-boostT);
1825 
1826  TVector3 vPt_ISR = vP_ISR.Vect();
1827  TVector3 vPt_I = vP_I.Vect();
1828  vPt_ISR -= vPt_ISR.Dot(boostZ.Unit())*boostZ.Unit();
1829  vPt_I -= vPt_I.Dot(boostZ.Unit())*boostZ.Unit();
1830 
1831  m_PTISR = vPt_ISR.Mag();
1832  m_RISR = -vPt_I.Dot(vPt_ISR.Unit()) / m_PTISR;
1833  m_PTI = vPt_I.Mag();
1834  m_dphiISRI = fabs(vPt_ISR.Angle(vPt_I));
1835 
1836  }//end INTERMEDIATE
1837 
1838  // Cutflow check
1839 
1840  cutFlowVector_str[0] = "No cuts ";
1841  cutFlowVector_str[1] = "3LLOW: Preselection ";
1842  cutFlowVector_str[2] = "3LLOW: 75 GeV < mll < 105 GeV ";
1843  cutFlowVector_str[3] = "3LLOW: mTW > 100 GeV ";
1844  cutFlowVector_str[4] = "3LLOW: m_HT4PP/m_H4PP > 0.9 ";
1845  cutFlowVector_str[5] = "3LLOW: m_H4PP > 250 GeV ";
1846  cutFlowVector_str[6] = "3LLOW: pT_PP/(pT_PP + HT_PP(3,1)) ";
1847  cutFlowVector_str[7] = "2L2JLOW: Preselection ";
1848  cutFlowVector_str[8] = "2L2JLOW: mll ";
1849  cutFlowVector_str[9] = "2L2JLOW: mjj ";
1850  cutFlowVector_str[10] = "2L2JLOW: HT_PP(1,1)/HT_PP(4,1) ";
1851  cutFlowVector_str[11] = "2L2JLOW: pT_PP/(pT_PP + HT_PP(4,1)) ";
1852  cutFlowVector_str[12] = "2L2JLOW: minDPhi ";
1853  cutFlowVector_str[13] = "2L2JLOW: HPP(4,1) ";
1854  cutFlowVector_str[14] = "2L2JINT: Preselection ";
1855  cutFlowVector_str[15] = "2L2JINT: mll ";
1856  cutFlowVector_str[16] = "2L2JINT: mjj ";
1857  cutFlowVector_str[17] = "2L2JINT: HT_PP(1,1)/HT_PP(4,1) ";
1858  cutFlowVector_str[18] = "2L2JINT: pT_PP/(pT_PP + HT_PP(4,1)) ";
1859  cutFlowVector_str[19] = "2L2JINT: minDPhi ";
1860  cutFlowVector_str[20] = "2L2JINT: HPP(4,1) ";
1861  cutFlowVector_str[21] = "2L2JHIGH: Preselection ";
1862  cutFlowVector_str[22] = "2L2JHIGH: mll ";
1863  cutFlowVector_str[23] = "2L2JHIGH: mjj ";
1864  cutFlowVector_str[24] = "2L2JHIGH: m_R_minH2P_minH3P>0.8";
1865  cutFlowVector_str[25] = "2L2JHIGH: m_RPT_HT5PP < 0.05 ";
1866  cutFlowVector_str[26] = "2L2JHIGH: 0.3 < minDPhiVP > 2.8 ";
1867  cutFlowVector_str[27] = "2L2JHIGH: m_H5PP>800. ";
1868  cutFlowVector_str[28] = "2L2JCOMP: Preselection ";
1869  cutFlowVector_str[29] = "2L2JCOMP: mZ ";
1870  cutFlowVector_str[30] = "2L2JCOMP: mJ ";
1871  cutFlowVector_str[31] = "2L2JCOMP: dPhi_ISR_I ";
1872  cutFlowVector_str[32] = "2L2JCOMP: R_ISR ";
1873  cutFlowVector_str[33] = "2L2JCOMP: p_ISRT ";
1874  cutFlowVector_str[34] = "2L2JCOMP: p_IT ";
1875  cutFlowVector_str[35] = "2L2JCOMP: pT_CM ";
1876  cutFlowVector_str[36] = "3LHIGH: Preselection ";
1877  cutFlowVector_str[37] = "3LHIGH: mll ";
1878  cutFlowVector_str[38] = "3LHIGH: mTW ";
1879  cutFlowVector_str[39] = "3LHIGH: m_HT4PP/m_H4PP ";
1880  cutFlowVector_str[40] = "3LHIGH: HPb(1,1)/HPb(2,1) ";
1881  cutFlowVector_str[41] = "3LHIGH: m_H4PP ";
1882  cutFlowVector_str[42] = "3LHIGH: pT_PP/(pT_PP + HT_PP(3,1)) ";
1883  cutFlowVector_str[43] = "3LCOMP: Preselection ";
1884  cutFlowVector_str[44] = "3LCOMP: mll ";
1885  cutFlowVector_str[45] = "3LCOMP: mTW ";
1886  cutFlowVector_str[46] = "3LCOMP: dPhi_ISRI ";
1887  cutFlowVector_str[47] = "3LCOMP: R_ISR ";
1888  cutFlowVector_str[48] = "3LCOMP: p_ISRT ";
1889  cutFlowVector_str[49] = "3LCOMP: p_IT ";
1890  cutFlowVector_str[50] = "3LCOMP: pT_CM ";
1891  cutFlowVector_str[51] = "3LINT: Preselection ";
1892  cutFlowVector_str[52] = "3LINT: mll ";
1893  cutFlowVector_str[53] = "3LINT: mTW ";
1894  cutFlowVector_str[54] = "3LINT: m_HT4PP/m_H4PP ";
1895  cutFlowVector_str[55] = "3LINT: HPb(1,1)/HPb(2,1) ";
1896  cutFlowVector_str[56] = "3LINT: m_H4PP ";
1897  cutFlowVector_str[57] = "3LINT: pT_PP/(pT_PP + HT_PP(3,1)) ";
1898 
1899  //std::cout << " m_is3Lep " << m_is3Lep << " m_is2Lep2Jet " << m_is2Lep2Jet << " m_is2L2JInt " << m_is2L2JInt << " m_is3LInt " << m_is3LInt << " m_is3Lep2Jet " << m_is3Lep2Jet << " m_is3Lep3Jet " << m_is3Lep3Jet << " m_is4Lep2Jet " << m_is4Lep2Jet << " m_is4Lep3Jet " << m_is4Lep3Jet << std::endl;
1900 
1901  //if(m_is3Lep)std::cout << " m_is3Lep " << m_is3Lep << " m_lept1sign " << m_lept1sign << " m_lept2sign " << m_lept2sign << " m_lept1Pt " << m_lept1Pt << " m_lept2Pt " << m_lept2Pt << " m_lept3Pt " << m_lept3Pt << " signalBJets.size() " << signalBJets.size() << " signalJets.size() " << signalJets.size() << std::endl;
1902 
1903  for(int j=0;j<NCUTS;j++){
1904 
1905  if( (j==0) ||
1906 
1907  (j==1 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>40. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()==0) ||
1908 
1909  (j==2 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>40. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()==0 && m_mll>75. && m_mll<105.) ||
1910 
1911  (j==3 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>40. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()==0 && m_mll>75. && m_mll<105. && m_mTW>100.) ||
1912 
1913  (j==4 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>40. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()==0 && m_mll>75. && m_mll<105. && m_mTW>100. && m_HT4PP/m_H4PP > 0.9) ||
1914 
1915  (j==5 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>40. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()==0 && m_mll>75. && m_mll<105. && m_mTW>100. && m_HT4PP/m_H4PP > 0.9 && m_H4PP > 250.) ||
1916 
1917  (j==6 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>40. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()==0 && m_mll>75. && m_mll<105. && m_mTW>100. && m_HT4PP/m_H4PP > 0.9 && m_H4PP > 250. && m_RPT_HT4PP < 0.05) ||
1918 
1919  /*cutFlowVector_str[7] = "2L2JLOW: Preselection ";
1920  cutFlowVector_str[8] = "2L2JLOW: 80 GeV < mll < 100 GeV";
1921  cutFlowVector_str[9] = "2L2JLOW: 70 GeV < mjj < 90 GeV ";
1922  cutFlowVector_str[10] = "2L2JLOW: HT_PP(1,1)/HT_PP(4,1) ";
1923  cutFlowVector_str[11] = "2L2JLOW: pT_PP/(pT_PP + HT_PP(4,1)) ";
1924  cutFlowVector_str[12] = "2L2JLOW: minDPhi ";
1925  cutFlowVector_str[13] = "2L2JLOW: HPP(4,1) ";*/
1926 
1927  (j==7 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()==2) ||
1928 
1929  (j==8 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()==2 && m_mll>80. && m_mll<100.) ||
1930 
1931  (j==9 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()==2 && m_mll>80. && m_mll<100. && m_mjj>70. && m_mjj<90.) ||
1932 
1933  (j==10 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()==2 && m_mll>80. && m_mll<100. && m_mjj>70. && m_mjj<90. && m_H2PP/m_H5PP>0.35 && m_H2PP/m_H5PP<0.6) ||
1934 
1935  (j==11 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()==2 && m_mll>80. && m_mll<100. && m_mjj>70. && m_mjj<90. && m_H2PP/m_H5PP>0.35 && m_H2PP/m_H5PP<0.6 && m_RPT_HT5PP<0.05) ||
1936 
1937  (j==12 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()==2 && m_mll>80. && m_mll<100. && m_mjj>70. && m_mjj<90. && m_H2PP/m_H5PP>0.35 && m_H2PP/m_H5PP<0.6 && m_RPT_HT5PP<0.05 && m_minDphi>2.4) ||
1938 
1939  (j==13 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()==2 && m_mll>80. && m_mll<100. && m_mjj>70. && m_mjj<90. && m_H2PP/m_H5PP>0.35 && m_H2PP/m_H5PP<0.6 && m_RPT_HT5PP<0.05 && m_minDphi>2.4 && m_H5PP>400.) ||
1940 
1941  /*cutFlowVector_str[14] = "2L2JINT: Preselection ";
1942  cutFlowVector_str[15] = "2L2JINT: mll ";
1943  cutFlowVector_str[16] = "2L2JINT: mjj ";
1944  cutFlowVector_str[17] = "2L2JINT: HT_PP(1,1)/HT_PP(4,1) ";
1945  cutFlowVector_str[18] = "2L2JINT: pT_PP/(pT_PP + HT_PP(4,1)) ";
1946  cutFlowVector_str[19] = "2L2JINT: minDPhi ";
1947  cutFlowVector_str[20] = "2L2JINT: HPP(4,1) ";*/
1948 
1949  (j==14 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2) ||
1950 
1951  (j==15 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100.) ||
1952 
1953  (j==16 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100.) ||
1954 
1955  (j==17 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100. && m_R_minH2P_minH3P>0.8 ) ||
1956 
1957  (j==18 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100. && m_R_minH2P_minH3P>0.8 && m_RPT_HT5PP<0.05) ||
1958 
1959  (j==19 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100. && m_R_minH2P_minH3P>0.8 && m_RPT_HT5PP<0.05 && m_dphiVP>0.6) ||
1960 
1961  (j==20 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100. && m_R_minH2P_minH3P>0.8 && m_RPT_HT5PP<0.05 && m_dphiVP>0.6 && m_H5PP>600.) ||
1962 
1963  /* cutFlowVector_str[21] = "2L2JHIGH: Preselection ";
1964  cutFlowVector_str[22] = "2L2JHIGH: mll ";
1965  cutFlowVector_str[23] = "2L2JHIGH: mjj ";
1966  cutFlowVector_str[24] = "2L2JHIGH: m_R_minH2P_minH3P>0.8";
1967  cutFlowVector_str[25] = "2L2JHIGH: m_RPT_HT5PP < 0.05 ";
1968  cutFlowVector_str[26] = "2L2JHIGH: 0.3 < minDPhiVP > 2.8 ";
1969  cutFlowVector_str[27] = "2L2JHIGH: m_H5PP>800. ";*/
1970 
1971  (j==21 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2) ||
1972 
1973  (j==22 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100.) ||
1974 
1975  (j==23 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100.) ||
1976 
1977  (j==24 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100. && m_R_minH2P_minH3P>0.8) ||
1978 
1979  (j==25 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100. && m_R_minH2P_minH3P>0.8 && m_RPT_HT5PP< 0.05) ||
1980 
1981  (j==26 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100. && m_R_minH2P_minH3P>0.8 && m_RPT_HT5PP< 0.05 && m_dphiVP>0.3 && m_dphiVP<2.8 ) ||
1982 
1983  (j==27 && m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100. && m_R_minH2P_minH3P>0.8 && m_RPT_HT5PP< 0.05 && m_dphiVP>0.3 && m_dphiVP<2.8 && m_H5PP>800.) ||
1984 
1985  /*cutFlowVector_str[28] = "2L2JCOMP: Preselection ";
1986  cutFlowVector_str[29] = "2L2JCOMP: mZ ";
1987  cutFlowVector_str[30] = "2L2JCOMP: mJ ";
1988  cutFlowVector_str[31] = "2L2JCOMP: dPhi_ISR_I ";
1989  cutFlowVector_str[32] = "2L2JCOMP: R_ISR ";
1990  cutFlowVector_str[33] = "2L2JCOMP: p_ISRT ";
1991  cutFlowVector_str[34] = "2L2JCOMP: p_IT ";
1992  cutFlowVector_str[35] = "2L2JCOMP: pT_CM ";*/
1993 
1994  (j==28 && m_is2L2JInt && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && m_NjS==2 && m_NjISR>0) ||
1995 
1996  (j==29 && m_is2L2JInt && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && m_NjS==2 && m_NjISR>0 && m_MZ>80. && m_MZ<100.) ||
1997 
1998  (j==30 && m_is2L2JInt && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && m_NjS==2 && m_NjISR>0 && m_MZ>80. && m_MZ<100. && m_MJ>50. && m_MJ<110.) ||
1999 
2000  (j==31 && m_is2L2JInt && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && m_NjS==2 && m_NjISR>0 && m_MZ>80. && m_MZ<100. && m_MJ>50. && m_MJ<110. && m_dphiISRI>2.8) ||
2001 
2002  (j==32 && m_is2L2JInt && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && m_NjS==2 && m_NjISR>0 && m_MZ>80. && m_MZ<100. && m_MJ>50. && m_MJ<110. && m_dphiISRI>2.8 && m_RISR > 0.40 && m_RISR < 0.75 ) ||
2003 
2004  (j==33 && m_is2L2JInt && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && m_NjS==2 && m_NjISR>0 && m_MZ>80. && m_MZ<100. && m_MJ>50. && m_MJ<110. && m_dphiISRI>2.8 && m_RISR > 0.40 && m_RISR < 0.75 && m_PTISR > 180.) ||
2005 
2006  (j==34 && m_is2L2JInt && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && m_NjS==2 && m_NjISR>0 && m_MZ>80. && m_MZ<100. && m_MJ>50. && m_MJ<110. && m_dphiISRI>2.8 && m_RISR > 0.40 && m_RISR < 0.75 && m_PTISR > 180. && m_PTI > 100.) ||
2007 
2008  (j==35 && m_is2L2JInt && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && m_NjS==2 && m_NjISR>0 && m_MZ>80. && m_MZ<100. && m_MJ>50. && m_MJ<110. && m_dphiISRI>2.8 && m_RISR > 0.40 && m_RISR < 0.75 && m_PTISR > 180. && m_PTI > 100. && m_PTCM < 20.) ||
2009 
2010 
2011  /* cutFlowVector_str[36] = "3LHIGH: Preselection ";
2012  cutFlowVector_str[37] = "3LHIGH: 75 GeV < mll < 105 GeV ";
2013  cutFlowVector_str[38] = "3LHIGH: mTW ";
2014  cutFlowVector_str[39] = "3LHIGH: m_HT4PP/m_H4PP ";
2015  cutFlowVector_str[40] = "3LHIGH: HPb(1,1)/HPb(2,1) ";
2016  cutFlowVector_str[41] = "3LHIGH: m_H4PP ";
2017  cutFlowVector_str[42] = "3LHIGH: pT_PP/(pT_PP + HT_PP(3,1)) ";*/
2018 
2019 
2020  (j==36 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>60. && m_lept3Pt>40. && signalBJets.size()==0 && signalJets.size()<3) ||
2021 
2022  (j==37 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>60. && m_lept3Pt>40. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105.) ||
2023 
2024  (j==38 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>60. && m_lept3Pt>40. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>150.) ||
2025 
2026  (j==39 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>60. && m_lept3Pt>40. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>150. && m_HT4PP/m_H4PP > 0.75) ||
2027 
2028  (j==40 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>60. && m_lept3Pt>40. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>150. && m_HT4PP/m_H4PP > 0.75 && m_R_minH2P_minH3P>0.8) ||
2029 
2030  (j==41 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>60. && m_lept3Pt>40. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>150. && m_R_minH2P_minH3P>0.8 && m_HT4PP/m_H4PP > 0.75 && m_H4PP > 550. ) ||
2031 
2032  (j==42 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>60. && m_lept3Pt>40. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>150. && m_R_minH2P_minH3P>0.8 && m_HT4PP/m_H4PP > 0.75 && m_H4PP > 550. && m_RPT_HT4PP < 0.2) ||
2033 
2034  /*cutFlowVector_str[43] = "3LCOMP: Preselection ";
2035  cutFlowVector_str[44] = "3LCOMP: mll ";
2036  cutFlowVector_str[45] = "3LCOMP: mTW ";
2037  cutFlowVector_str[46] = "3LCOMP: dPhi_ISRI ";
2038  cutFlowVector_str[47] = "3LCOMP: R_ISR ";
2039  cutFlowVector_str[48] = "3LCOMP: p_ISRT ";
2040  cutFlowVector_str[49] = "3LCOMP: p_IT ";
2041  cutFlowVector_str[50] = "3LCOMP: pT_CM ";*/
2042 
2043  (j==43 && m_is3LInt && ((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign))) && m_lept1Pt>25. && m_lept2Pt>25. && m_lept3Pt>20. && signalBJets.size()==0 && signalJets.size()<4) ||
2044 
2045  (j==44 && m_is3LInt && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>25. && m_lept2Pt>25. && m_lept3Pt>20. && signalBJets.size()==0 && signalJets.size()<4 && m_mll>75. && m_mll<105.) ||
2046 
2047  (j==45 && m_is3LInt && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>25. && m_lept2Pt>25. && m_lept3Pt>20. && signalBJets.size()==0 && signalJets.size()<4 && m_mll>75. && m_mll<105. && m_mTW>100.) ||
2048 
2049  (j==46 && m_is3LInt && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>25. && m_lept2Pt>25. && m_lept3Pt>20. && signalBJets.size()==0 && signalJets.size()<4 && m_mll>75. && m_mll<105. && m_mTW>100. && m_dphiISRI>2.0) ||
2050 
2051  (j==47 && m_is3LInt && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>25. && m_lept2Pt>25. && m_lept3Pt>20. && signalBJets.size()==0 && signalJets.size()<4 && m_mll>75. && m_mll<105. && m_mTW>100. && m_dphiISRI>2.0 && m_RISR>0.55 && m_RISR<1.0) ||
2052 
2053  (j==48 && m_is3LInt && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>25. && m_lept2Pt>25. && m_lept3Pt>20. && signalBJets.size()==0 && signalJets.size()<4 && m_mll>75. && m_mll<105. && m_mTW>100. && m_dphiISRI>2.0 && m_RISR>0.55 && m_RISR<1.0 && m_PTISR>100.) ||
2054 
2055  (j==49 && m_is3LInt && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>25. && m_lept2Pt>25. && m_lept3Pt>20. && signalBJets.size()==0 && signalJets.size()<4 && m_mll>75. && m_mll<105. && m_mTW>100. && m_dphiISRI>2.0 && m_RISR>0.55 && m_RISR<1.0 && m_PTISR>100. && m_PTI>80.) ||
2056 
2057  (j==50 && m_is3LInt && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>25. && m_lept2Pt>25. && m_lept3Pt>20. && signalBJets.size()==0 && signalJets.size()<4 && m_mll>75. && m_mll<105. && m_mTW>100. && m_dphiISRI>2.0 && m_RISR>0.55 && m_RISR<1.0 && m_PTISR>100. && m_PTI>80. && m_PTCM<25.) ||
2058 
2059  /* cutFlowVector_str[51] = "3LINT: Preselection ";
2060  cutFlowVector_str[52] = "3LINT: mll ";
2061  cutFlowVector_str[53] = "3LINT: mTW ";
2062  cutFlowVector_str[54] = "3LINT: m_HT4PP/m_H4PP ";
2063  cutFlowVector_str[55] = "3LINT: HPb(1,1)/HPb(2,1) ";
2064  cutFlowVector_str[56] = "3LINT: m_H4PP ";
2065  cutFlowVector_str[57] = "3LINT: pT_PP/(pT_PP + HT_PP(3,1)) ";*/
2066 
2067  (j==51 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>50. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()<3) ||
2068 
2069  (j==52 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>50. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105.) ||
2070 
2071  (j==53 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>50. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>130.) ||
2072 
2073  (j==54 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>50. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>130. && m_HT4PP/m_H4PP > 0.8 ) ||
2074 
2075  (j==55 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>50. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>130. && m_HT4PP/m_H4PP > 0.8 && m_R_minH2P_minH3P>0.75) ||
2076 
2077  (j==56 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>50. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>130. && m_HT4PP/m_H4PP > 0.8 && m_R_minH2P_minH3P>0.75 && m_H4PP > 450. ) ||
2078 
2079  (j==57 && m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>50. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>130. && m_HT4PP/m_H4PP > 0.8 && m_R_minH2P_minH3P>0.75 && m_H4PP > 450. && m_RPT_HT4PP < 0.15)
2080 
2081  )cutFlowVector[j]++;
2082  }
2083 
2084  // Now apply the signal region cuts
2085 
2086  if(m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100. && m_R_minH2P_minH3P>0.8 && m_RPT_HT5PP< 0.05 && m_dphiVP>0.3 && m_dphiVP<2.8 && m_H5PP>800.)_num2L2JHIGH++;
2087 
2088  if(m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()>=2 && m_mll>80. && m_mll<100. && m_mjj>60. && m_mjj<100. && m_R_minH2P_minH3P>0.8 && m_RPT_HT5PP<0.05 && m_dphiVP>0.6 && m_dphiVP<2.6 && m_H5PP>600.)_num2L2JINT++;
2089 
2090 
2091  if(m_is2Lep2Jet && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && signalJets.size()==2 && m_mll>80. && m_mll<100. && m_mjj>70. && m_mjj<90. && m_H2PP/m_H5PP>0.35 && m_H2PP/m_H5PP<0.6 && m_RPT_HT5PP<0.05 && m_minDphi>2.4 && m_H5PP>400.)_num2L2JLOW++;
2092 
2093  if(m_is2L2JInt && m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign) && m_lept1Pt>25. && m_lept2Pt>25. && m_jet1Pt>30. && m_jet2Pt>30. && signalBJets.size()==0 && m_NjS==2 && m_NjISR>0 && m_MZ>80. && m_MZ<100. && m_MJ>50. && m_MJ<110. && m_dphiISRI>2.8 && m_RISR > 0.40 && m_RISR < 0.75 && m_PTISR > 180. && m_PTI > 100. && m_PTCM < 20.)_num2L2JCOMP++;
2094 
2095  if(m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>60. && m_lept3Pt>40. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>150. && m_HT4PP/m_H4PP > 0.75 && m_R_minH2P_minH3P>0.8 && m_H4PP > 550. && m_RPT_HT4PP < 0.2)_num3LHIGH++;
2096 
2097  if(m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>50. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()<3 && m_mll>75. && m_mll<105. && m_mTW>130. && m_HT4PP/m_H4PP > 0.8 && m_R_minH2P_minH3P>0.75 && m_H4PP > 450. && m_RPT_HT4PP < 0.15)_num3LINT++;
2098 
2099  if(m_is3LInt && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>25. && m_lept2Pt>25. && m_lept3Pt>20. && signalBJets.size()==0 && signalJets.size()<4 && m_mll>75. && m_mll<105. && m_mTW>100. && m_dphiISRI>2.0 && m_RISR>0.55 && m_RISR<1.0 && m_PTISR>100. && m_PTI>80. && m_PTCM<25.)_num3LCOMP++;
2100 
2101  if(m_is3Lep && (((m_lept1sign*m_lept2sign<0 && abs(m_lept1sign)==abs(m_lept2sign)) || (m_lept1sign*m_lept3sign<0 && abs(m_lept1sign)==abs(m_lept3sign)) || (m_lept2sign*m_lept3sign<0 && abs(m_lept2sign)==abs(m_lept3sign)))) && m_lept1Pt>60. && m_lept2Pt>40. && m_lept3Pt>30. && signalBJets.size()==0 && signalJets.size()==0 && m_mll>75. && m_mll<105. && m_mTW>100. && m_HT4PP/m_H4PP > 0.9 && m_H4PP > 250. && m_RPT_HT4PP < 0.05)_num3LLOW++;
2102 
2103  return;
2104 
2105  } // end analyze method
2106 
2108  void combine(const Analysis* other)
2109  {
2110  const Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb* specificOther
2111  = dynamic_cast<const Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb*>(other);
2112 
2113  if (NCUTS != specificOther->NCUTS) NCUTS = specificOther->NCUTS;
2114 
2115  for (int j=0; j<NCUTS; j++)
2116  {
2117  cutFlowVector[j] += specificOther->cutFlowVector[j];
2118  cutFlowVector_str[j] = specificOther->cutFlowVector_str[j];
2119  }
2120 
2121  _num3LLOW+= specificOther->_num3LLOW;
2122  }
2123 
2124 
2125  virtual void collect_results() {
2126 /*
2127  double scale_by=1.;
2128  cout << "------------------------------------------------------------------------------------------------------------------------------ "<<endl;
2129  cout << "CUT FLOW: ATLAS 13 TeV 3 lep low mass RJ signal region "<<endl;
2130  cout << "------------------------------------------------------------------------------------------------------------------------------"<<endl;
2131  cout << right << setw(40) << "CUT" << setw(20) << "RAW" << setw(20) << "SCALED"
2132  << setw(20) << "%" << setw(20) << "clean adj RAW"<< setw(20) << "clean adj %" << endl;
2133  for (int j=0; j<NCUTS; j++) {
2134  cout << right << setw(40) << cutFlowVector_str[j].c_str() << setw(20)
2135  << cutFlowVector[j] << setw(20) << cutFlowVector[j]*scale_by << setw(20)
2136  << 100.*cutFlowVector[j]/cutFlowVector[0] << "%" << setw(20)
2137  << cutFlowVector[j]*scale_by << setw(20) << 100.*cutFlowVector[j]/cutFlowVector[0]<< "%" << endl;
2138  }
2139  cout << "------------------------------------------------------------------------------------------------------------------------------ "<<endl;
2140 */
2141  // add_result(SignalRegionData("SR label", n_obs, {s, s_sys}, {b, b_sys}));
2142  add_result(SignalRegionData("2L2JHIGH", 0, {_num2L2JHIGH, 0.}, {1.9, 0.8}));
2143  add_result(SignalRegionData("2L2JINT", 1, {_num2L2JINT, 0.}, {2.4, 0.9}));
2144  add_result(SignalRegionData("2L2JLOW", 19, {_num2L2JLOW, 0.}, {8.4, 5.8}));
2145  add_result(SignalRegionData("2L2JCOMP", 11, {_num2L2JCOMP, 0.}, {2.7, 2.7}));
2146  add_result(SignalRegionData("3LHIGH", 2, {_num3LHIGH, 0.}, {1.1, 0.5}));
2147  add_result(SignalRegionData("3LINT", 1, {_num3LINT, 0.}, {2.3, 0.5}));
2148  add_result(SignalRegionData("3LLOW", 20, {_num3LLOW, 0.}, {10., 2.0}));
2149  add_result(SignalRegionData("3LCOMP", 12, {_num3LCOMP, 0.}, {3.9, 1.0}));
2150 
2151  return;
2152  }
2153 
2154 
2155  protected:
2156  void analysis_specific_reset() {
2157  _num2L2JHIGH=0;
2158  _num2L2JINT=0;
2159  _num2L2JLOW=0;
2160  _num2L2JCOMP=0;
2161  _num3LHIGH=0;
2162  _num3LINT=0;
2163  _num3LLOW=0;
2164  _num3LCOMP=0;
2165 
2166  std::fill(cutFlowVector.begin(), cutFlowVector.end(), 0);
2167  }
2168 
2169  }; // end class Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb
2170 
2171  DEFINE_ANALYSIS_FACTORY(ATLAS_13TeV_RJ3L_lowmass_36invfb)
2172 
2173 
2174  //
2175  // Derived analysis class for the RJ3L_lowmass SRs
2176  //
2177  class Analysis_ATLAS_13TeV_RJ3L_2Lep2Jets_36invfb : public Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb {
2178  public:
2179  Analysis_ATLAS_13TeV_RJ3L_2Lep2Jets_36invfb() {
2180  set_analysis_name("ATLAS_13TeV_RJ3L_2Lep2Jets_36invfb");
2181  }
2182  virtual void collect_results() {
2183  add_result(SignalRegionData("2L2JHIGH", 0, {_num2L2JHIGH, 0.}, {1.9, 0.8}));
2184  add_result(SignalRegionData("2L2JINT", 1, {_num2L2JINT, 0.}, {2.4, 0.9}));
2185  add_result(SignalRegionData("2L2JLOW", 19, {_num2L2JLOW, 0.}, {8.4, 5.8}));
2186  add_result(SignalRegionData("2L2JCOMP", 11, {_num2L2JCOMP, 0.}, {2.7, 2.7}));
2187  }
2188 
2189  };
2190  // Factory fn
2191  DEFINE_ANALYSIS_FACTORY(ATLAS_13TeV_RJ3L_2Lep2Jets_36invfb)
2192 
2193  class Analysis_ATLAS_13TeV_RJ3L_3Lep_36invfb : public Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb {
2194  public:
2195  Analysis_ATLAS_13TeV_RJ3L_3Lep_36invfb() {
2196  set_analysis_name("ATLAS_13TeV_RJ3L_3Lep_36invfb");
2197  }
2198  virtual void collect_results() {
2199  add_result(SignalRegionData("3LHIGH", 2, {_num3LHIGH, 0.}, {1.1, 0.5}));
2200  add_result(SignalRegionData("3LINT", 1, {_num3LINT, 0.}, {2.3, 0.5}));
2201  add_result(SignalRegionData("3LLOW", 20, {_num3LLOW, 0.}, {10., 2.0}));
2202  add_result(SignalRegionData("3LCOMP", 12, {_num3LCOMP, 0.}, {3.9, 1.0}));
2203  }
2204 
2205  };
2206  // Factory fn
2207  DEFINE_ANALYSIS_FACTORY(ATLAS_13TeV_RJ3L_3Lep_36invfb)
2208 
2209 
2210  } // end namespace ColliderBit
2211 } // end namespace Gambit
2212 
2213 #endif
2214 #endif
Gambit::piped_errors
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
Gambit::ColliderBit::has_tag
bool has_tag(const HEPUtils::BinnedFn2D< double > &effmap, double eta, double pt)
Randomly get a tag result (can be anything) from a 2D |eta|-pT efficiency map.
Definition: Utils.hpp:173
Gambit::ColliderBit::sortLepByPT_RJ3L
bool sortLepByPT_RJ3L(HEPUtils::Particle *lep1, HEPUtils::Particle *lep2)
Definition: Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp:41
DEFINE_ANALYSIS_FACTORY
#define DEFINE_ANALYSIS_FACTORY(ANAME)
For analysis factory function definition.
Definition: Analysis.hpp:124
Gambit::ColliderBit::ATLAS::applyMuonEff
void applyMuonEff(std::vector< HEPUtils::Particle *> &muons)
Randomly filter the supplied particle list by parameterised muon efficiency.
Definition: ATLASEfficiencies.hpp:83
combine_hdf5.pid
pid
Definition: combine_hdf5.py:273
Gambit::Piped_exceptions::request
void request(std::string origin, std::string message)
Request an exception.
Definition: exceptions.cpp:527
LOCAL_INFO
#define LOCAL_INFO
Definition: local_info.hpp:34
Gambit::piped_warnings
Piped_exceptions piped_warnings
Global instance of Piped_exceptions class for warnings.
run
void run(bool, bool, bool, int, double, double, int, int, int, int, int, double, char[], int, int[], bool, bool, bool, bool, double, int, double(*)(double *, int, int, void *), void(*)(int, int, int, double *, double *, double *, double, double, double, void *), void *)
std
STL namespace.
Gambit::ColliderBit::sortByPT_RJ3L
bool sortByPT_RJ3L(HEPUtils::Jet *jet1, HEPUtils::Jet *jet2)
Definition: Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp:40
b
START_MODEL b
Definition: demo.hpp:235
Gambit::ColliderBit::SortLeptons
bool SortLeptons(const pair< TLorentzVector, int > lv1, const pair< TLorentzVector, int > lv2)
Definition: Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp:44
mt2_bisect.h
Gambit::ColliderBit::Analysis
A class for collider analyses within ColliderBit.
Definition: Analysis.hpp:41
Gambit::str
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
Gambit::ColliderBit::SignalRegionData
A simple container for the result of one signal region from one analysis.
Definition: AnalysisData.hpp:54
int
DS5_MSPCTM DS_INTDOF int
Definition: DarkBit_rollcall.hpp:142
M
START_MODEL M
Definition: CosmoNuisanceModels.hpp:45
colouring.combine
def combine(region_file, pip_file)
Definition: colouring.py:169
Gambit::ColliderBit::SortJets
bool SortJets(const TLorentzVector jv1, const TLorentzVector jv2)
Definition: Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp:50
Gambit::ColliderBit::ATLAS::applyElectronEff
void applyElectronEff(std::vector< HEPUtils::Particle *> &electrons)
Randomly filter the supplied particle list by parameterised electron efficiency.
Definition: ATLASEfficiencies.hpp:60
Gambit
TODO: see if we can use this one:
Definition: Analysis.hpp:33
Analysis.hpp
Class for ColliderBit analyses.
ATLASEfficiencies.hpp
Functions that do super fast ATLAS detector simulation based on four-vector smearing.