ATLAS Offline Software
TtresChi2.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
3  */
4 
5 /* **********************************************************************
6 * Class : TTBarLeptonJetsBuilder_chi2
7 * Author : LPC-ClermontFerrand team
8 ************************************************************************/
9 
10 #include "TopEventReconstructionTools/TtresChi2.h"
11 #include "TopEventReconstructionTools/TtresNeutrinoBuilder.h"
12 
13 //_________________________________________________________________________________________________
14 TtresChi2::TtresChi2(std::string units) {
15  m_debug = 0;
16  m_WReco = ROTATION;
17  if (units == "MeV") m_Units = 1000.;
18  else if (units == "GeV") m_Units = 1.;
19  else std::cerr << "WARNING in TtresChi2: units different from GeV or MeV" << std::endl;
20  m_UsePtDiff = PTDIFF;
21  m_Btag = AFFECTBTAG;
22  m_RunMode = RUNDEFAULT;
23  m_NeutrinoBuilder = new TtresNeutrinoBuilder(units);
24  Init(DATA2012SUMMER2013PT100);
25  res_chi2All = -1.;
26  res_chi2WH = -1.;
27  res_chi2TopH = -1.;
28  res_chi2TopL = -1.;
29  res_Mtl = -1.;
30  res_Mwh = -1.;
31  res_Mth = -1.;
32  res_Mtt = -1.;
33  m_WhChi2Value = -1;
34  m_ThWhChi2Value = -1;
35  m_TlChi2Value = -1;
36  m_PtDiffChi2Value = -1;
37  m_category = -1;
38 }
39 
40 //_________________________________________________________________________________________________
41 TtresChi2::~TtresChi2() {
42  if (m_debug > 0) std::cout << "in destructor " << std::endl;
43  if (m_NeutrinoBuilder) delete m_NeutrinoBuilder;
44 }
45 
46 //_________________________________________________________________________________________________
47 void TtresChi2::Init(Chi2Version version, double highJetMass) {
48  m_highJetMass = highJetMass;
49 
50  switch (version) {
51  case DATA2012SUMMER2013PT100:
52  //Rel17.2: values obtained for Z' 0.5TeV-2 TeV including njets>=5 events
53  // LC jets + Jet Area correction
54  // Applying a cut at 100 GeV on the pT of the leading jet
55  MjjP = 82.74 * m_Units;
56  SMjjP = 9.524 * m_Units;
57 
58  m_TopMinusW_had_mean = 89.99 * m_Units;
59  m_TopMinusW_had_sigma = 13.91 * m_Units;
60 
61  m_Top_lep_mean = 166.6 * m_Units;
62  m_Top_lep_sigma = 17.34 * m_Units;
63 
64  m_PtDiff_mean = -1.903 * m_Units;
65  m_PtDiff_sigma = 47.24 * m_Units;
66 
67  m_PtDiffRel_mean = 0.002678;
68  m_PtDiffRel_sigma = 0.07662;
69 
70  m_PtDiffRelMass_mean = 0.001315;
71  m_PtDiffRelMass_sigma = 0.05092;
72 
73  MTHJJ = 171.8 * m_Units;
74  STHJJ = 16.04 * m_Units;
75  break;
76 
77  case DATA2012SUMMER2013:
78  //Rel17.2: values obtained for Z' 0.5TeV-2 TeV including njets>=5 events
79  // LC jets + Jet Area correction
80  MjjP = 82.41 * m_Units;
81  SMjjP = 9.553 * m_Units;
82 
83  m_TopMinusW_had_mean = 89.01 * m_Units;
84  m_TopMinusW_had_sigma = 15.67 * m_Units;
85 
86  m_Top_lep_mean = 166 * m_Units;
87  m_Top_lep_sigma = 17.48 * m_Units;
88 
89  m_PtDiff_mean = 0.4305 * m_Units;
90  m_PtDiff_sigma = 46.07 * m_Units;
91 
92  m_PtDiffRel_mean = 0.00144;
93  m_PtDiffRel_sigma = 0.07837;
94 
95  m_PtDiffRelMass_mean = 0.0004532;
96  m_PtDiffRelMass_sigma = 0.05112;
97 
98  MTHJJ = 171.8 * m_Units;
99  STHJJ = 16.04 * m_Units;
100  break;
101 
102  case DATA2012AUTOMN2012:
103  //Rel17.2: values obtained for Z' 0.5TeV-2 TeV including njets>=5 events
104  // LC jets
105  MjjP = 83.93 * m_Units;
106  SMjjP = 10.76 * m_Units;
107 
108  m_TopMinusW_had_mean = 91.05 * m_Units;
109  m_TopMinusW_had_sigma = 14.23 * m_Units;
110 
111  m_Top_lep_mean = 168.2 * m_Units;
112  m_Top_lep_sigma = 20.57 * m_Units;
113 
114  m_PtDiff_mean = -8.709 * m_Units;
115  m_PtDiff_sigma = 55. * m_Units;
116 
117  MTHJJ = 173.5 * m_Units;
118  STHJJ = 16.25 * m_Units;
119  break;
120 
121  case DATA2011SUMMER2012:
122  //Rel17.0: values obtained for Z' 1TeV-2 TeV including njets>=5 events
123  // EM+JES jets
124  MjjP = 83.18 * m_Units;
125  SMjjP = 10.67 * m_Units;
126 
127  m_TopMinusW_had_mean = 90.87 * m_Units;
128  m_TopMinusW_had_sigma = 12.78 * m_Units;
129 
130  m_Top_lep_mean = 167.6 * m_Units;
131  m_Top_lep_sigma = 20.53 * m_Units;
132 
133  m_PtDiff_mean = -7.381 * m_Units;
134  m_PtDiff_sigma = 63.96 * m_Units;
135 
136  MTHJJ = 172.2 * m_Units;
137  STHJJ = 16.8 * m_Units;
138  break;
139 
140  default:
141  std::cerr << "TtresChi2::Init chi2 is not defined!!" << std::endl;
142  break;
143  }
144  return;
145 }
146 
147 //_________________________________________________________________________________________________
148 void TtresChi2::Reset() {
149  res_chi2All = -1.;
150  res_chi2WH = -1.;
151  res_chi2TopH = -1.;
152  res_chi2TopL = -1.;
153  res_Mwh = -1.;
154  res_Mth = -1.;
155  res_Mtl = -1.;
156  res_Mtt = -1.;
157  m_WhChi2Value = -1;
158  m_ThWhChi2Value = -1;
159  m_TlChi2Value = -1;
160  m_PtDiffChi2Value = -1;
161  m_category = -1;
162  m_nChi2Values = 0;
163  m_chi2Values.clear();
164  m_chi2Wh_Values.clear();
165  m_chi2Th_Values.clear();
166  m_chi2ThWh_Values.clear();
167  m_chi2Tl_Values.clear();
168  m_i_Wq1.clear();
169  m_i_Wq2.clear();
170  m_i_Thb.clear();
171  m_i_Tlb.clear();
172  m_i_n.clear();
173  m_chi2PtDiff_Values.clear();
174  m_PtDiff_Values.clear();
175  m_Wh_Values.clear();
176  m_ThWh_Values.clear();
177  m_Th_Values.clear();
178  m_Tl_Values.clear();
179  m_chi2Values.reserve(10000);
180  m_chi2Wh_Values.reserve(10000);
181  m_chi2Th_Values.reserve(10000);
182  m_chi2ThWh_Values.reserve(10000);
183  m_chi2Tl_Values.reserve(10000);
184  m_chi2PtDiff_Values.reserve(10000);
185  m_i_Wq1.reserve(10000);
186  m_i_Wq2.reserve(10000);
187  m_i_Thb.reserve(10000);
188  m_i_Tlb.reserve(10000);
189  m_i_n.reserve(10000);
190  m_PtDiff_Values.reserve(10000);
191  m_Wh_Values.reserve(10000);
192  m_ThWh_Values.reserve(10000);
193  m_Th_Values.reserve(10000);
194  m_Tl_Values.reserve(10000);
195  return;
196 }
197 
198 //_________________________________________________________________________________________________
199 std::vector<TLorentzVector*> TtresChi2::GetNeutrinos(TLorentzVector* L, TLorentzVector* MET) {
200  std::vector<TLorentzVector*> neutrinoVector;
201  if (m_WReco ==
202  ROTATION) neutrinoVector = m_NeutrinoBuilder->candidatesFromWMass_Rotation(L, MET, false /*m_useSmallestPz*/);
203  else if (m_WReco ==
204  REALPART) neutrinoVector =
205  m_NeutrinoBuilder->candidatesFromWMass_RealPart(L, MET, false /*m_useSmallestPz*/);
206  else if (m_WReco == DECREASING) neutrinoVector = m_NeutrinoBuilder->candidatesFromWMass_Scaling(L, MET);
207  return neutrinoVector;
208 }
209 
210 //_________________________________________________________________________________________________
211 bool TtresChi2::findMinChiSquare_HMelseLM(TLorentzVector* L, const std::vector<TLorentzVector*>* jetVector,
212  const std::vector<bool>* isJetBtagged, TLorentzVector* MET, int& i_q1_W,
213  int& i_q2_W, int& i_b_had, int& i_b_lep, int& ign1, double& chi2ming1,
214  double& chi2ming1H, double& chi2ming1L) {
215  if (m_debug > 0) std::cout << "entering TtresChi2::findMinChiSquare_HMelseLM()" << std::endl;
216  bool status = false;
217 
218  status = this->findMinChiSquare_HighMass(L, jetVector, isJetBtagged, MET, i_q1_W, i_q2_W, i_b_lep, ign1, chi2ming1,
219  chi2ming1H, chi2ming1L);
220 
221  if (status && TMath::Log10(chi2ming1) < 0.9) return status; // have found a good combination for HM
222 
223  // Otherwise try the LM chi2
224  status = this->findMinChiSquare(L, jetVector, isJetBtagged, MET, i_q1_W, i_q2_W, i_b_had, i_b_lep, ign1, chi2ming1,
225  chi2ming1H, chi2ming1L);
226 
227  return status;
228 }
229 
230 //_________________________________________________________________________________________________
231 bool TtresChi2::findMinChiSquare_LMelseHM(TLorentzVector* L, const std::vector<TLorentzVector*>* jetVector,
232  const std::vector<bool>* isJetBtagged, TLorentzVector* MET, int& i_q1_W,
233  int& i_q2_W, int& i_b_had, int& i_b_lep, int& ign1, double& chi2ming1,
234  double& chi2ming1H, double& chi2ming1L) {
235  if (m_debug > 0) std::cout << "entering TtresChi2::findMinChiSquare_LMelseHM()" << std::endl;
236  bool status = false;
237 
238  status = this->findMinChiSquare(L, jetVector, isJetBtagged, MET, i_q1_W, i_q2_W, i_b_had, i_b_lep, ign1, chi2ming1,
239  chi2ming1H, chi2ming1L);
240 
241  if (status && TMath::Log10(chi2ming1) < 0.9) return status; // have found a good combination for LM
242 
243  // Otherwise try the HM chi2
244  status = this->findMinChiSquare_HighMass(L, jetVector, isJetBtagged, MET, i_q1_W, i_q2_W, i_b_lep, ign1, chi2ming1,
245  chi2ming1H, chi2ming1L);
246 
247  return status;
248 }
249 
250 //_________________________________________________________________________________________________
251 bool TtresChi2::findMinChiSquare_AllRanges(TLorentzVector* L, const std::vector<TLorentzVector*>* jetVector,
252  const std::vector<bool>* isJetBtagged, TLorentzVector* MET, int& i_q1_W,
253  int& i_q2_W, int& i_b_had, int& i_b_lep, int& ign1, double& chi2ming1,
254  double& chi2ming1H, double& chi2ming1L) {
255  if (m_debug > 0) std::cout << "entering TtresChi2::findMinChiSquare_AllRanges()" << std::endl;
256  bool status = false;
257  float massMax = 0.;
258 
259  for (unsigned int i = 0; i < jetVector->size();
260  ++i) if (jetVector->at(i)->M() > massMax) massMax = jetVector->at(i)->M();
261 
262  if (massMax >= m_highJetMass * m_Units) status = this->findMinChiSquare_HighMass(L, jetVector, isJetBtagged, MET, i_q1_W, i_q2_W, i_b_lep, ign1, chi2ming1, chi2ming1H, chi2ming1L);
263  else status = this->findMinChiSquare(L, jetVector, isJetBtagged, MET, i_q1_W, i_q2_W, i_b_had, i_b_lep, ign1, chi2ming1, chi2ming1H, chi2ming1L);
264 
265  return status;
266 }
267 
268 //_________________________________________________________________________________________________
269 bool TtresChi2::findMinChiSquare(TLorentzVector* L, const std::vector<TLorentzVector*>* jetVector, const std::vector<bool>* isJetBtagged, TLorentzVector* MET, int& i_q1_W, int& i_q2_W, int& i_b_had, int& i_b_lep, int& ign1, double& chi2ming1, double& chi2ming1H, double& chi2ming1L) {
270  if (m_debug > 0) std::cout << "entering TtresChi2::findMinChiSquare()" << std::endl;
271  //_________________________________________________________________________________________________
272  i_q1_W = -1;
273  i_q2_W = -1;
274  i_b_had = -1;
275  i_b_lep = -1;
276  ign1 = -1;
277  chi2ming1 = 1e7;
278  chi2ming1H = 1e7;
279  chi2ming1L = 1e7;
280  double chi2ming1WH = 1e7;
281  int n_jets = jetVector->size();
282  int n_bjets = 0;
283 
284  this->Reset();
285 
286  if (L == NULL) {
287  std::cerr << "ERROR : TtresChi2::findMinChiSquare: Lepton is NULL" << std::endl;
288  return false;
289  }
290 
291  if (m_Btag == AFFECTBTAG) {
292  for (unsigned int ib = 0; ib < (unsigned int) isJetBtagged->size(); ++ib) {
293  if (isJetBtagged->at(ib)) n_bjets++;
294  }
295  }
296 
297  //_________________________________________________________________________________________________
298  std::vector<TLorentzVector*> neutrinoVector = GetNeutrinos(L, MET);
299 
300  //_________________________________________________________________________________________________
301  for (unsigned int i = 0; i < (unsigned int) n_jets; i++) {
302  const TLorentzVector* I = jetVector->at(i);
303  for (unsigned int j = i + 1; j < (unsigned int) n_jets; j++) {
304  const TLorentzVector* J = jetVector->at(j);
305  for (unsigned int k = 0; k < (unsigned int) n_jets; k++) {
306  if (k != i && k != j) {
307  const TLorentzVector* K = jetVector->at(k);
308  TLorentzVector TopH = (*I) + (*J) + (*K);
309  TLorentzVector Whad = (*I) + (*J);
310  for (unsigned int n = 0; n < neutrinoVector.size(); n++) {
311  TLorentzVector* N = neutrinoVector[n];
312 
313  TLorentzVector Wlv = (*N) + (*L);
314  for (unsigned int m = 0; m < (unsigned int) n_jets; m++) {
315  if (m != i && m != j && m != k) {
316  const TLorentzVector* M = jetVector->at(m);
317  TLorentzVector Tlv = Wlv + (*M);
318  TLorentzVector Tt = Tlv + TopH;
319  double chi2WH = pow((Whad.M() - MjjP) / SMjjP, 2);
320  double chi2H = chi2WH + pow((TopH.M() - Whad.M() - m_TopMinusW_had_mean) / m_TopMinusW_had_sigma, 2);
321  double chi2L = pow((Tlv.M() - m_Top_lep_mean) / m_Top_lep_sigma, 2);
322 
323  double chi2tempg1 = chi2H + chi2L;
324  if (m_UsePtDiff == PTDIFF) {
325  double chi2Diff = pow((TopH.Pt() - Tlv.Pt() - m_PtDiff_mean) / m_PtDiff_sigma, 2);
326  chi2tempg1 += chi2Diff;
327  } else if (m_UsePtDiff == PTDIFFREL) {
328  double chi2DiffRel = pow(((TopH.Pt() - Tlv.Pt()) / (TopH.Pt() + Tlv.Pt()) - m_PtDiffRel_mean) / m_PtDiffRel_sigma, 2);
329  chi2tempg1 += chi2DiffRel;
330  } else if (m_UsePtDiff == PTDIFFMASS) {
331  double chi2DiffRelMass = pow(((TopH.Pt() - Tlv.Pt()) / (Tt.M()) - m_PtDiffRelMass_mean) / m_PtDiffRelMass_sigma, 2);
332  chi2tempg1 += chi2DiffRelMass;
333  }
334 
335  // does this combination contain a b
336  // k -> hadronic b
337  // m -> leptonic b
338  int n_bJetsComb = 0;
339  if (isJetBtagged->at(k)) {
340  n_bJetsComb++;
341  }
342  // second b
343  if (isJetBtagged->at(m)) {
344  n_bJetsComb++;
345  }
346 
347  bool passBtag = false;
348  if (m_Btag == STDBTAG || m_Btag == NO_BTAGHM) {
349  if (n_bJetsComb > 0) {
350  passBtag = true;
351  }
352  } else if (m_Btag == AFFECTBTAG) {
353  if (n_bjets == 0) {
354  passBtag = true;
355  } else if (n_bjets == 1) {
356  if (n_bJetsComb == 1) {
357  passBtag = true;
358  }
359  } else if (n_bjets >= 2) {
360  if (n_bJetsComb == 2) {
361  passBtag = true;
362  }
363  }
364  } else if (m_Btag == NO_BTAG) {
365  passBtag = true;
366  }
367 
368  if (passBtag && m_RunMode == RUNSTUDY) {
369  m_chi2Values.push_back(chi2tempg1);
370  m_chi2Wh_Values.push_back(chi2WH);
371  m_chi2ThWh_Values.push_back(chi2H - chi2WH);
372  m_chi2Tl_Values.push_back(chi2L);
373  m_chi2PtDiff_Values.push_back(chi2tempg1 - chi2H - chi2L);
374  m_i_Wq1.push_back(i);
375  m_i_Wq2.push_back(j);
376  m_i_Thb.push_back(k);
377  m_i_Tlb.push_back(m);
378  m_i_n.push_back(n);
379  if (m_UsePtDiff == PTDIFF) m_PtDiff_Values.push_back((TopH.Pt() - Tlv.Pt()) / m_Units);
380  else if (m_UsePtDiff == PTDIFFMASS) m_PtDiff_Values.push_back((TopH.Pt() - Tlv.Pt()) / (Tt.M()));
381  m_Wh_Values.push_back(Whad.M() / m_Units);
382  m_ThWh_Values.push_back((TopH.M() - Whad.M()) / m_Units);
383  m_Tl_Values.push_back((Tlv.M()) / m_Units);
384  m_nChi2Values++;
385  }
386 
387  //#################
388  //Original chi2 method
389  //#################
390  if (chi2tempg1 < chi2ming1) {
391  if (passBtag) {
392  chi2ming1 = chi2tempg1;
393  chi2ming1H = chi2H;
394  chi2ming1L = chi2L;
395  m_WhChi2Value = chi2WH;
396  m_ThWhChi2Value = chi2H - chi2WH;
397  m_TlChi2Value = chi2L;
398  m_PtDiffChi2Value = chi2tempg1 - chi2H - chi2L;
399 
400  i_q1_W = i;
401  i_q2_W = j;
402  i_b_had = k;
403  i_b_lep = m;
404  ign1 = n;
405  res_Mtl = Tlv.M();
406  res_Mwh = Whad.M();
407  res_Mth = TopH.M();
408  res_Mtt = Tt.M();
409  res_Tt = Tt;
410 
411  //============================
412  // bjet category splitting
413  //============================
414  if (n_bjets == 0) {
415  m_category = 0;
416  } else if (n_bjets == 1) {
417  if (isJetBtagged->at(k)) {//b-jet on the hadronic side
418  m_category = 2;
419  } else if (isJetBtagged->at(m)) {//b-jet on the leptonic side
420  m_category = 3;
421  } else {
422  std::cerr << "<!> In TtresChi2::findMinChiSquare : cannot find the corresponding category." << std::endl;
423  }
424  } else if (n_bjets >= 2) {
425  if (isJetBtagged->at(k)) {//b-jet on the hadronic side
426  if (isJetBtagged->at(m)) {//b-jet on the leptonic side too
427  m_category = 1;
428  } else {//b-jet only on the hadronic side
429  m_category = 2;
430  }
431  } else if (isJetBtagged->at(m)) {
432  m_category = 3;
433  } else {
434  std::cerr << "<!> In TtresChi2::findMinChiSquare : cannot find the corresponding category." << std::endl;
435  }
436  }
437  }
438  } // end of case a minimal chisquare was found
439  } // end of case unique combination
440  } //end of loop over jets
441  } // end of loop over all neutrino solutions
442  } // end of case this is the electron channel
443  } // end of loop over jets k
444  } // end of loop over jets
445  } // end of loop over jets
446 
447  bool status = false;
448  if (i_q1_W != -1 && i_q2_W != -1 && i_b_had != -1 && i_b_lep != -1 && ign1 != -1) {
449  status = true;
450  }
451 
452  res_chi2All = chi2ming1;
453  res_chi2WH = chi2ming1WH;
454  res_chi2TopH = chi2ming1H;
455  res_chi2TopL = chi2ming1L;
456 
457  for (unsigned int i = 0; i < neutrinoVector.size(); i++) {
458  delete neutrinoVector[i];
459  }
460  neutrinoVector.clear();
461  return status;
462 }
463 
464 //_________________________________________________________________________________________________
465 bool TtresChi2::findMinChiSquare_HighMass(TLorentzVector* L, const std::vector<TLorentzVector*>* jetVector, const std::vector<bool>* isJetBtagged, TLorentzVector* MET, int& i_q1_W, int& i_q2_W, int& i_b_lep, int& ign1, double& chi2ming1, double&
466  chi2ming1H, double& chi2ming1L) {
467  if (m_debug > 0) std::cout << "entering TtresChi2::findMinChiSquare_HighMass()" << std::endl;
468  //_________________________________________________________________________________________________
469  i_q1_W = -1;
470  i_q2_W = -1;
471  i_b_lep = -1;
472  ign1 = -1;
473  chi2ming1 = 1e7;
474  chi2ming1H = 1e7;
475  chi2ming1L = 1e7;
476  double HMtop = -1e6;
477  int n_jets = jetVector->size();
478  int n_bjets = 0;
479 
480  this->Reset();
481 
482  if (L == NULL) {
483  std::cerr << "ERROR : TtresChi2::findMinChiSquare: Lepton is NULL" << std::endl;
484  return false;
485  }
486 
487  if (m_Btag == AFFECTBTAG) {
488  for (unsigned int ib = 0; ib < (unsigned int) isJetBtagged->size(); ++ib) {
489  if (isJetBtagged->at(ib)) n_bjets++;
490  }
491  }
492 
493  //_________________________________________________________________________________________________
494  std::vector<TLorentzVector*> neutrinoVector = GetNeutrinos(L, MET);
495 
496  //_________________________________________________________________________________________________
497  for (unsigned int i = 0; i < (unsigned int) n_jets; i++) {
498  const TLorentzVector* I = jetVector->at(i);
499  for (unsigned int j = i + 1; j < (unsigned int) n_jets; j++) {
500  const TLorentzVector* J = jetVector->at(j);
501  if (I->M() > m_highJetMass * m_Units || J->M() > m_highJetMass * m_Units) {
502  TLorentzVector TopH = (*I) + (*J);
503 
504  for (unsigned int n = 0; n < neutrinoVector.size(); n++) {
505  TLorentzVector* N = neutrinoVector[n];
506 
507  TLorentzVector Wlv = (*N) + (*L);
508  for (unsigned int m = 0; m < (unsigned int) n_jets; m++) {
509  if (m != i && m != j) {
510  const TLorentzVector* M = jetVector->at(m);
511  TLorentzVector Tlv = Wlv + (*M);
512  TLorentzVector Tt = Tlv + TopH;
513 
514  double HMtoptemp = I->M();
515  double chi2H = pow((TopH.M() - MTHJJ) / STHJJ, 2);
516  double chi2L = pow((Tlv.M() - m_Top_lep_mean) / m_Top_lep_sigma, 2);
517 
518  double chi2tempg1 = chi2H + chi2L;
519  if (m_UsePtDiff == PTDIFF) {
520  double chi2Diff = pow((TopH.Pt() - Tlv.Pt() - m_PtDiff_mean) / m_PtDiff_sigma, 2);
521  chi2tempg1 += chi2Diff;
522  } else if (m_UsePtDiff == PTDIFFREL) {
523  double chi2DiffRel = pow(((TopH.Pt() - Tlv.Pt()) / (TopH.Pt() + Tlv.Pt()) - m_PtDiffRel_mean) / m_PtDiffRel_sigma, 2);
524  chi2tempg1 += chi2DiffRel;
525  } else if (m_UsePtDiff == PTDIFFMASS) {
526  double chi2DiffRelMass = pow(((TopH.Pt() - Tlv.Pt()) / (Tt.M()) - m_PtDiffRelMass_mean) / m_PtDiffRelMass_sigma, 2);
527  chi2tempg1 += chi2DiffRelMass;
528  }
529 
530  // does this combination contain a b
531  // k -> hadronic b
532  // m -> leptonic b
533  int n_bJetsComb = 0;
534  // check whether one of the taggers is passed
535  // first b
536  if (isJetBtagged->at(i) || isJetBtagged->at(j)) {
537  n_bJetsComb++;
538  }
539  // second b
540  if (isJetBtagged->at(m)) {
541  n_bJetsComb++;
542  }
543 
544  bool passBtag = false;
545  if (m_Btag == STDBTAG) {
546  if (n_bJetsComb > 0) {
547  passBtag = true;
548  }
549  } else if (m_Btag == AFFECTBTAG) {
550  if (n_bjets == 0) {
551  passBtag = true;
552  } else if (n_bjets == 1) {
553  if (n_bJetsComb == 1) {
554  passBtag = true;
555  }
556  } else if (n_bjets >= 2) {
557  if (n_bJetsComb == 2) {
558  passBtag = true;
559  }
560  }
561  } else if (m_Btag == NO_BTAGHM || m_Btag == NO_BTAG) {
562  passBtag = true;
563  }
564 
565  if (passBtag && m_RunMode == RUNSTUDY) {
566  m_chi2Values.push_back(chi2tempg1);
567  m_chi2Th_Values.push_back(chi2H);
568  m_chi2Tl_Values.push_back(chi2L);
569  m_chi2PtDiff_Values.push_back(chi2tempg1 - chi2H - chi2L);
570  m_i_Wq1.push_back(i);
571  m_i_Wq2.push_back(j);
572  m_i_Tlb.push_back(m);
573  m_i_n.push_back(n);
574  if (m_UsePtDiff == PTDIFF) m_PtDiff_Values.push_back((TopH.Pt() - Tlv.Pt()) / m_Units);
575  else if (m_UsePtDiff == PTDIFFMASS) m_PtDiff_Values.push_back((TopH.Pt() - Tlv.Pt()) / (Tt.M()));
576  m_Th_Values.push_back((TopH.M()) / m_Units);
577  m_Tl_Values.push_back((Tlv.M()) / m_Units);
578  m_nChi2Values++;
579  }
580 
581  //#################
582  //Original chi2 method
583  //#################
584  if (chi2tempg1 < chi2ming1 && HMtoptemp >= HMtop) {
585  if (passBtag) {
586  HMtop = HMtoptemp;
587  chi2ming1 = chi2tempg1;
588  chi2ming1H = chi2H;
589  chi2ming1L = chi2L;
590  m_ThWhChi2Value = chi2H;
591  m_TlChi2Value = chi2L;
592  m_PtDiffChi2Value = chi2tempg1 - chi2H - chi2L;
593 
594  i_q1_W = i;
595  i_q2_W = j;
596  i_b_lep = m;
597  ign1 = n;
598  res_Mtl = Tlv.M();
599  //res_Mwh is not defined here
600  res_Mth = TopH.M();
601  res_Mtt = Tt.M();
602  res_Tt = Tt;
603 
604  //============================
605  // bjet category splitting
606  //============================
607  if (n_bjets == 0) {
608  m_category = 0;
609  } else if (n_bjets == 1) {
610  if (isJetBtagged->at(i) || isJetBtagged->at(j)) {//b-jet on the hadronic side
611  m_category = 2;
612  } else if (isJetBtagged->at(m)) {//b-jet on the leptonic side
613  m_category = 3;
614  } else {
615  std::cerr << "<!> In TtresChi2::findMinChiSquare : cannot find the corresponding category." << std::endl;
616  }
617  } else if (n_bjets >= 2) {
618  if (isJetBtagged->at(i) || isJetBtagged->at(m)) {//b-jet on the hadronic side
619  if (isJetBtagged->at(m)) {//b-jet on the leptonic side too
620  m_category = 1;
621  } else {//b-jet only on the hadronic side
622  m_category = 2;
623  }
624  } else if (isJetBtagged->at(m)) {
625  m_category = 3;
626  } else {
627  std::cerr << "<!> In TtresChi2::findMinChiSquare : cannot find the corresponding category." << std::endl;
628  }
629  }
630  } //end of btag requirement
631  } // end of case a minimal chisquare was found
632  } // end of case unique combination
633  } //end of loop over jets
634  } // end of loop over all neutrino solutions
635  } // end of case this is the electron channel
636  } // end of loop over jets
637  } // end of loop over jets
638 
639  bool status = false;
640  if (i_q1_W != -1 && i_q2_W != -1 && i_b_lep != -1 && ign1 != -1) {
641  status = true;
642  }
643 
644  res_chi2All = chi2ming1;
645  res_chi2TopH = chi2ming1H;
646  res_chi2TopL = chi2ming1L;
647 
648 
649  for (unsigned int i = 0; i < neutrinoVector.size(); i++) {
650  delete neutrinoVector[i];
651  }
652  neutrinoVector.clear();
653  return status;
654 }
655 
656 //_________________________________________________________________________________________________
657 bool TtresChi2::findMinChiSquare_VeryHighMass(TLorentzVector* L, const std::vector<TLorentzVector*>* jetVector, const std::vector<bool>* isJetBtagged, TLorentzVector* MET, int& i_q2_W, int& i_b_lep, int& ign1, double& chi2ming1, double& chi2ming1H, double& chi2ming1L) {
658  if (m_debug > 0) std::cout << "entering TtresChi2::findMinChiSquare_VeryHighMass() " << std::endl;
659  //_________________________________________________________________________________________________
660  i_q2_W = -1;
661  i_b_lep = -1;
662  ign1 = -1;
663  chi2ming1 = 1e7;
664  chi2ming1H = 1e7;
665  chi2ming1L = 1e7;
666  double HMtop = -1e6;
667  int n_jets = jetVector->size();
668  this->Reset();
669 
670  if (L == NULL) {
671  std::cerr << "ERROR : TtresChi2::findMinChiSquare: Lepton is NULL" << std::endl;
672  return false;
673  }
674 
675  std::vector<TLorentzVector*> neutrinoVector = GetNeutrinos(L, MET);
676 
677  //_________________________________________________________________________________________________
678  for (unsigned int i = 0; i < (unsigned int) n_jets; i++) {
679  const TLorentzVector* I = jetVector->at(i);
680  TLorentzVector TopH = (*I);
681 
682  if (I->M() > 150. * m_Units) {
683  for (unsigned int n = 0; n < neutrinoVector.size(); n++) {
684  TLorentzVector* N = neutrinoVector[n];
685 
686  TLorentzVector Wlv = (*N) + (*L);
687  for (unsigned int m = 0; m < (unsigned int) n_jets; m++) {
688  if (m != i) {
689  const TLorentzVector* M = jetVector->at(m);
690  TLorentzVector Tlv = Wlv + (*M);
691  TLorentzVector Tt = Tlv + TopH;
692 
693  double HMtoptemp = I->M();
694  double chi2H = 1.;
695  double chi2L = pow((Tlv.M() - m_Top_lep_mean) / m_Top_lep_sigma, 2);
696  double chi2Diff = pow((TopH.Pt() - Tlv.Pt() - m_PtDiff_mean) / m_PtDiff_sigma, 2);
697  double chi2tempg1 = chi2L;
698 
699  if (m_UsePtDiff) chi2tempg1 += chi2Diff;
700 
701  //#################
702  //Original chi2 method
703  //#################
704  if (chi2tempg1 < chi2ming1 && HMtoptemp >= HMtop) {
705  // does this combination contain a b
706  // k -> hadronic b
707  // m -> leptonic b
708  int n_bJets = 0;
709  // check whether one of the taggers is passed
710  // first b
711  if (isJetBtagged->at(i)) {
712  n_bJets++;
713  }
714  // second b
715  if (isJetBtagged->at(m)) {
716  n_bJets++;
717  }
718 
719  if (n_bJets >= 1) {
720  HMtop = HMtoptemp;
721  chi2ming1 = chi2tempg1;
722  chi2ming1H = chi2H;
723  chi2ming1L = chi2L;
724  i_q2_W = i;
725  i_b_lep = m;
726  ign1 = n;
727  res_Mtl = Tlv.M();
728  //res_Mwh is not defined here
729  res_Mth = TopH.M();
730  res_Mtt = Tt.M();
731  res_Tt = Tt;
732  }
733  } // end of case a minimal chisquare was found
734  } // end of case unique combination
735  } //end of loop over jets
736  } // end of loop over all neutrino solutions
737  } // end of high mass
738  } // end of loop over jets
739 
740  bool status = false;
741  if (i_q2_W != -1 && i_b_lep != -1 && ign1 != -1) {
742  status = true;
743  }
744 
745  res_chi2All = chi2ming1;
746  res_chi2TopH = chi2ming1H;
747  res_chi2TopL = chi2ming1L;
748 
749  for (unsigned int i = 0; i < neutrinoVector.size(); i++) {
750  delete neutrinoVector[i];
751  }
752  neutrinoVector.clear();
753  return status;
754 }
TtresChi2::m_i_Thb
std::vector< int > m_i_Thb
Definition: TtresChi2.h:71
TtresChi2::TtresChi2
TtresChi2(std::string units="GeV")
Definition: TtresChi2.cxx:14
TtresChi2::m_PtDiffRelMass_sigma
double m_PtDiffRelMass_sigma
Definition: TtresChi2.h:41
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
TtresChi2::m_TlChi2Value
double m_TlChi2Value
Definition: TtresChi2.h:42
TtresChi2::REALPART
@ REALPART
Definition: TtresChi2.h:21
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
TtresChi2::m_i_Wq1
std::vector< int > m_i_Wq1
Definition: TtresChi2.h:69
TtresChi2::m_PtDiffRelMass_mean
double m_PtDiffRelMass_mean
Definition: TtresChi2.h:41
TtresChi2::MTHJJ
double MTHJJ
Definition: TtresChi2.h:41
TtresChi2::Reset
void Reset()
Definition: TtresChi2.cxx:148
TtresChi2::PTDIFFREL
@ PTDIFFREL
Definition: TtresChi2.h:27
TtresChi2::m_chi2Tl_Values
std::vector< double > m_chi2Tl_Values
Definition: TtresChi2.h:60
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
TtresChi2::m_i_n
std::vector< int > m_i_n
Definition: TtresChi2.h:73
TtresChi2::m_highJetMass
double m_highJetMass
Definition: TtresChi2.h:50
TtresChi2::m_PtDiff_mean
double m_PtDiff_mean
Definition: TtresChi2.h:40
TtresChi2::res_chi2WH
double res_chi2WH
Definition: TtresChi2.h:53
TtresChi2::findMinChiSquare
bool findMinChiSquare(TLorentzVector *, const std::vector< TLorentzVector * > *, const std::vector< bool > *, TLorentzVector *, int &, int &, int &, int &, int &, double &, double &, double &)
Definition: TtresChi2.cxx:269
PlotCalibFromCool.ib
ib
Definition: PlotCalibFromCool.py:419
TtresChi2::res_Mth
double res_Mth
Definition: TtresChi2.h:53
TtresChi2::findMinChiSquare_AllRanges
bool findMinChiSquare_AllRanges(TLorentzVector *, const std::vector< TLorentzVector * > *, const std::vector< bool > *, TLorentzVector *, int &, int &, int &, int &, int &, double &, double &, double &)
Definition: TtresChi2.cxx:251
TtresChi2::DATA2012AUTOMN2012
@ DATA2012AUTOMN2012
Definition: TtresChi2.h:24
TtresChi2::res_Mtl
double res_Mtl
Definition: TtresChi2.h:53
TtresChi2::NO_BTAG
@ NO_BTAG
Definition: TtresChi2.h:30
TtresChi2::NO_BTAGHM
@ NO_BTAGHM
Definition: TtresChi2.h:30
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
TtresChi2::m_PtDiffRel_sigma
double m_PtDiffRel_sigma
Definition: TtresChi2.h:41
TtresChi2::m_ThWh_Values
std::vector< double > m_ThWh_Values
Definition: TtresChi2.h:65
TtresChi2::Chi2Version
Chi2Version
Definition: TtresChi2.h:23
TtresChi2::m_TopMinusW_had_mean
double m_TopMinusW_had_mean
Definition: TtresChi2.h:40
TtresChi2::m_chi2Wh_Values
std::vector< double > m_chi2Wh_Values
Definition: TtresChi2.h:57
TtresChi2::findMinChiSquare_VeryHighMass
bool findMinChiSquare_VeryHighMass(TLorentzVector *, const std::vector< TLorentzVector * > *, const std::vector< bool > *, TLorentzVector *, int &, int &, int &, double &, double &, double &)
Definition: TtresChi2.cxx:657
TtresChi2::ROTATION
@ ROTATION
Definition: TtresChi2.h:21
TtresChi2::RUNSTUDY
@ RUNSTUDY
Definition: TtresChi2.h:33
TtresNeutrinoBuilder
Definition: TtresNeutrinoBuilder.h:16
TtresChi2::m_PtDiffRel_mean
double m_PtDiffRel_mean
Definition: TtresChi2.h:41
TtresChi2::~TtresChi2
virtual ~TtresChi2()
Definition: TtresChi2.cxx:41
TtresChi2::findMinChiSquare_LMelseHM
bool findMinChiSquare_LMelseHM(TLorentzVector *, const std::vector< TLorentzVector * > *, const std::vector< bool > *, TLorentzVector *, int &, int &, int &, int &, int &, double &, double &, double &)
Definition: TtresChi2.cxx:231
TtresChi2::m_debug
int m_debug
Definition: TtresChi2.h:43
TtresChi2::m_Units
double m_Units
Definition: TtresChi2.h:48
TtresChi2::STHJJ
double STHJJ
Definition: TtresChi2.h:41
TtresChi2::m_Wh_Values
std::vector< double > m_Wh_Values
Definition: TtresChi2.h:64
TtresChi2::m_WReco
WReco m_WReco
Definition: TtresChi2.h:46
TtresChi2::m_chi2Values
std::vector< double > m_chi2Values
Definition: TtresChi2.h:56
lumiFormat.i
int i
Definition: lumiFormat.py:92
TtresChi2::m_nChi2Values
int m_nChi2Values
Definition: TtresChi2.h:75
TtresChi2::DATA2011SUMMER2012
@ DATA2011SUMMER2012
Definition: TtresChi2.h:24
beamspotman.n
n
Definition: beamspotman.py:731
TtresChi2::m_category
int m_category
Definition: TtresChi2.h:51
TtresChi2::m_chi2PtDiff_Values
std::vector< double > m_chi2PtDiff_Values
Definition: TtresChi2.h:61
TCS::MET
@ MET
Definition: Trigger/TrigT1/L1Topo/L1TopoCommon/L1TopoCommon/Types.h:16
TtresChi2::res_Tt
TLorentzVector res_Tt
Definition: TtresChi2.h:54
TtresChi2::MjjP
double MjjP
Definition: TtresChi2.h:40
TtresChi2::m_Tl_Values
std::vector< double > m_Tl_Values
Definition: TtresChi2.h:67
TtresChi2::m_PtDiff_sigma
double m_PtDiff_sigma
Definition: TtresChi2.h:41
TtresChi2::m_Btag
Btag m_Btag
Definition: TtresChi2.h:45
TtresChi2::m_NeutrinoBuilder
TtresNeutrinoBuilder * m_NeutrinoBuilder
Definition: TtresChi2.h:49
TtresChi2::m_TopMinusW_had_sigma
double m_TopMinusW_had_sigma
Definition: TtresChi2.h:40
TtresChi2::PTDIFFMASS
@ PTDIFFMASS
Definition: TtresChi2.h:27
TtresChi2::m_Th_Values
std::vector< double > m_Th_Values
Definition: TtresChi2.h:66
TtresChi2::m_i_Wq2
std::vector< int > m_i_Wq2
Definition: TtresChi2.h:70
TtresChi2::m_UsePtDiff
PtDiff m_UsePtDiff
Definition: TtresChi2.h:44
TtresChi2::m_RunMode
RunMode m_RunMode
Definition: TtresChi2.h:47
TtresChi2::RUNDEFAULT
@ RUNDEFAULT
Definition: TtresChi2.h:33
TtresChi2::SMjjP
double SMjjP
Definition: TtresChi2.h:40
TtresChi2::m_ThWhChi2Value
double m_ThWhChi2Value
Definition: TtresChi2.h:42
TtresChi2::PTDIFF
@ PTDIFF
Definition: TtresChi2.h:27
TtresChi2::findMinChiSquare_HMelseLM
bool findMinChiSquare_HMelseLM(TLorentzVector *, const std::vector< TLorentzVector * > *, const std::vector< bool > *, TLorentzVector *, int &, int &, int &, int &, int &, double &, double &, double &)
Definition: TtresChi2.cxx:211
TtresChi2::GetNeutrinos
std::vector< TLorentzVector * > GetNeutrinos(TLorentzVector *, TLorentzVector *)
Definition: TtresChi2.cxx:199
TtresChi2::m_chi2Th_Values
std::vector< double > m_chi2Th_Values
Definition: TtresChi2.h:58
TtresChi2::Init
void Init(Chi2Version version=DATA2012SUMMER2013, double highJetMass=60.0)
Definition: TtresChi2.cxx:47
TtresChi2::m_chi2ThWh_Values
std::vector< double > m_chi2ThWh_Values
Definition: TtresChi2.h:59
TtresChi2::DATA2012SUMMER2013
@ DATA2012SUMMER2013
Definition: TtresChi2.h:24
perfmonmt-refit.units
string units
Definition: perfmonmt-refit.py:77
TtresChi2.h
get_generator_info.version
version
Definition: get_generator_info.py:33
TtresChi2::res_chi2All
double res_chi2All
Definition: TtresChi2.h:53
TtresChi2::findMinChiSquare_HighMass
bool findMinChiSquare_HighMass(TLorentzVector *, const std::vector< TLorentzVector * > *, const std::vector< bool > *, TLorentzVector *, int &, int &, int &, int &, double &, double &, double &)
Definition: TtresChi2.cxx:465
TtresChi2::m_Top_lep_sigma
double m_Top_lep_sigma
Definition: TtresChi2.h:40
TtresChi2::AFFECTBTAG
@ AFFECTBTAG
Definition: TtresChi2.h:30
TtresChi2::m_WhChi2Value
double m_WhChi2Value
Definition: TtresChi2.h:42
TtresChi2::m_Top_lep_mean
double m_Top_lep_mean
Definition: TtresChi2.h:40
TtresChi2::res_chi2TopH
double res_chi2TopH
Definition: TtresChi2.h:53
TtresChi2::DATA2012SUMMER2013PT100
@ DATA2012SUMMER2013PT100
Definition: TtresChi2.h:24
TtresNeutrinoBuilder::candidatesFromWMass_Scaling
bool candidatesFromWMass_Scaling(const TLorentzVector *, double, Double_t, std::vector< TLorentzVector * > &)
TtresNeutrinoBuilder::candidatesFromWMass_Rotation
std::vector< TLorentzVector * > candidatesFromWMass_Rotation(const TLorentzVector *, const Double_t, const Double_t, const bool)
Definition: TtresNeutrinoBuilder.cxx:206
TtresChi2::m_PtDiffChi2Value
double m_PtDiffChi2Value
Definition: TtresChi2.h:42
merge.status
status
Definition: merge.py:17
TtresChi2::res_Mtt
double res_Mtt
Definition: TtresChi2.h:53
TtresNeutrinoBuilder::candidatesFromWMass_RealPart
std::vector< TLorentzVector * > candidatesFromWMass_RealPart(const TLorentzVector *, const Double_t, const Double_t, const bool)
Definition: TtresNeutrinoBuilder.cxx:288
TtresChi2::res_chi2TopL
double res_chi2TopL
Definition: TtresChi2.h:53
I
#define I(x, y, z)
Definition: MD5.cxx:116
TtresChi2::m_PtDiff_Values
std::vector< double > m_PtDiff_Values
Definition: TtresChi2.h:63
TtresChi2::STDBTAG
@ STDBTAG
Definition: TtresChi2.h:30
TtresChi2::m_i_Tlb
std::vector< int > m_i_Tlb
Definition: TtresChi2.h:72
TtresNeutrinoBuilder.h
TtresChi2::res_Mwh
double res_Mwh
Definition: TtresChi2.h:53
fitman.k
k
Definition: fitman.py:528
TtresChi2::DECREASING
@ DECREASING
Definition: TtresChi2.h:21
Generated on Tue Jul 2 2024 21:32:06 for ATLAS Offline Software by   doxygen 1.8.18