ATLAS Offline Software
MissingMassProb.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // Class handling the probability calculation of the MissingMassCalculator
6 // author Michael Huebner <michael.huebner@no.spam.cern.ch>
7 
8 // Local include(s):
11 #include <cmath>
12 
14 
15 namespace {
16  constexpr double GEV = 1000.0;
17 }
18 
19 using namespace DiTauMassTools;
20 using ROOT::Math::PtEtaPhiMVector;
22 
23 // The wrapper functions make use of the ignore template defined in HelperFunctions.h.
24 // This is to avoid warnings during compilation.
25 // All wrapper functions need to have the same structure such that they can be
26 // handled within the probList vectors and the "user" from the Calculator only
27 // has to call Prob->apply() to have a common handling of all probability calculations.
28 
29 double MissingMassProb::MetProbabilityWrapper( MissingMassProb* prob, MissingMassInput& preparedInput, const int & tau_type1, const int & tau_type2, const PtEtaPhiMVector & tauvec1, const PtEtaPhiMVector & tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector & nuvec2 ){
30  ignore(tau_type1);
31  ignore(tau_type2);
32  ignore(tauvec1);
33  ignore(tauvec2);
34  ignore(nuvec1);
35  ignore(nuvec2);
36  return prob->MetProbability(preparedInput, 1.,1.,1.,1.);
37 }
38 
39 double MissingMassProb::mEtAndTauProbabilityWrapper( MissingMassProb* prob, MissingMassInput& preparedInput, const int & tau_type1, const int & tau_type2, const PtEtaPhiMVector & tauvec1, const PtEtaPhiMVector & tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector & nuvec2 ){
40  ignore(tau_type1);
41  ignore(tau_type2);
42  ignore(tauvec1);
43  ignore(tauvec2);
44  ignore(nuvec1);
45  ignore(nuvec2);
46  return prob->mEtAndTauProbability(preparedInput);
47 }
48 
49 double MissingMassProb::dTheta3d_probabilityFastWrapper( MissingMassProb* prob, MissingMassInput& preparedInput, const int & tau_type1, const int & tau_type2, const PtEtaPhiMVector & tauvec1, const PtEtaPhiMVector & tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector & nuvec2 ){
50  double Prob = 1.;
51  if (tau_type1>=0) {
52  PtEtaPhiMVector totalTau1;
53  totalTau1+=tauvec1;
54  totalTau1+=nuvec1;
55  const double tau1_tmpp = totalTau1.P();
56  const double angle1 = Angle(nuvec1,tauvec1);
57  Prob*=prob->dTheta3d_probabilityFast(preparedInput, tau_type1, angle1, tau1_tmpp);
58  }
59  if (tau_type2>=0) {
60  PtEtaPhiMVector totalTau2;
61  totalTau2+=tauvec2;
62  totalTau2+=nuvec2;
63  const double tau2_tmpp = totalTau2.P();
64  const double angle2 = Angle(nuvec2,tauvec2);
65  Prob*=prob->dTheta3d_probabilityFast(preparedInput, tau_type2, angle2, tau2_tmpp);
66  }
67  return Prob;
68 }
69 
70 double MissingMassProb::TauProbabilityWrapper( MissingMassProb* prob, MissingMassInput& preparedInput, const int & tau_type1, const int & tau_type2, const PtEtaPhiMVector & tauvec1, const PtEtaPhiMVector & tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector & nuvec2 ){
71  if( prob->GetUseHT() || (preparedInput.m_tauTypes==TauTypes::hh) ) {
72  return prob->TauProbability(preparedInput, tau_type1, tauvec1, nuvec1, tau_type2, tauvec2, nuvec2, preparedInput.m_MetVec.R()); // customized prob for Njet25=0
73  } else {
74  return prob->TauProbability(preparedInput, tau_type1, tauvec1, nuvec1, tau_type2, tauvec2, nuvec2);
75  }
76 }
77 
78 double MissingMassProb::MnuProbabilityWrapper( MissingMassProb* prob, MissingMassInput& preparedInput, const int & tau_type1, const int & tau_type2, const PtEtaPhiMVector & tauvec1, const PtEtaPhiMVector & tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector & nuvec2 ){
79  ignore(tauvec1);
80  ignore(tauvec2);
81  if(prob->GetUseMnuProbability()==1){
82  if(preparedInput.m_tauTypes==TauTypes::ll) return prob->MnuProbability(preparedInput, nuvec1.M())*prob->MnuProbability(preparedInput, nuvec2.M()); // lep-lep
83  else if(tau_type1==8 && (tau_type2>=0 && tau_type2<=5)) return prob->MnuProbability(preparedInput, nuvec1.M()); // lep-had: tau1==lepton
84  else if((tau_type1>=0 && tau_type1<=5) && tau_type2==8) return prob->MnuProbability(preparedInput, nuvec2.M()); // lep-had: tau2==lepton
85  else {
86  Warning("DiTauMassTools", "something went really wrong in MNuProb...");
87  return 1.;
88  }
89  } else {
90  return 1.;
91  }
92 }
93 
94 double MissingMassProb::dTheta3d_probabilityNewWrapper( MissingMassProb* prob, MissingMassInput& preparedInput, const int & tau_type1, const int & tau_type2, const PtEtaPhiMVector & tauvec1, const PtEtaPhiMVector & tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector & nuvec2 ){
95  ignore(preparedInput);
96  double Prob = 1.;
97  if (tau_type1>=0) {
98  PtEtaPhiMVector totalTau1;
99  totalTau1+=tauvec1;
100  totalTau1+=nuvec1;
101  const double angle1 = Angle(nuvec1,tauvec1);
102  double prob_tmp = 1e-10;
103  if (angle1!=0.) prob_tmp=prob->GetFormulaAngle1()->Eval(angle1);
104  if (prob_tmp<=0.) prob_tmp=1e-10;
105  Prob*=prob_tmp;
106  }
107  if (tau_type2>=0) {
108  PtEtaPhiMVector totalTau2;
109  totalTau2+=tauvec2;
110  totalTau2+=nuvec2;
111  const double angle2 = Angle(nuvec2,tauvec2);
112  double prob_tmp = 1e-10;
113  if (angle2!=0.) prob_tmp=prob->GetFormulaAngle2()->Eval(angle2);
114  if (prob_tmp<=0.) prob_tmp=1e-10;
115  Prob*=prob_tmp;
116  }
117  // very rare cases where parametrisation flips
118  if (std::isnan(Prob)) Prob = 0.;
119  return Prob;
120 }
121 
122 double MissingMassProb::TauProbabilityNewWrapper( MissingMassProb* prob, MissingMassInput& preparedInput, const int & tau_type1, const int & tau_type2, const PtEtaPhiMVector & tauvec1, const PtEtaPhiMVector & tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector & nuvec2 ){
123  ignore(preparedInput);
124  ignore(tau_type1);
125  ignore(tau_type2);
126  double Prob = 1.;
127  double R1 = nuvec1.P()/tauvec1.P();
128  double R2 = nuvec2.P()/tauvec2.P();
129  Prob*=prob->GetFormulaRatio1()->Eval(R1);
130  Prob*=prob->GetFormulaRatio2()->Eval(R2);
131  // not observed, just a protection
132  if (std::isnan(Prob)) Prob = 0.;
133  return Prob;
134 }
135 
136 double MissingMassProb::MnuProbabilityNewWrapper( MissingMassProb* prob, MissingMassInput& preparedInput, const int & tau_type1, const int & tau_type2, const PtEtaPhiMVector & tauvec1, const PtEtaPhiMVector & tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector & nuvec2 ){
137  ignore(tauvec1);
138  ignore(tauvec2);
139  double Prob = 1.;
140  if(prob->GetUseMnuProbability()==1){
141  if(preparedInput.m_tauTypes==TauTypes::ll) Prob*=(prob->GetFormulaNuMass()->Eval(nuvec1.M())*prob->GetFormulaNuMass()->Eval(nuvec2.M()));
142  else if(tau_type1==8 && (tau_type2>=0 && tau_type2<=5)) Prob*=prob->GetFormulaNuMass()->Eval(nuvec1.M());
143  else if((tau_type1>=0 && tau_type1<=5) && tau_type2==8) Prob*=prob->GetFormulaNuMass()->Eval(nuvec2.M());
144  else {
145  Warning("DiTauMassTools", "something went really wrong in MNuProb...");
146  }
147  }
148  // not observed, just a protection
149  if (std::isnan(Prob)) Prob = 0.;
150  return Prob;
151 }
152 
154  if (m_paramVectorNuMass.size() > 0){
155  for(int i=0; i<m_formulaNuMass->GetNpar(); i++){
156  m_formulaNuMass->SetParameter(i, m_paramVectorNuMass[0]->GetParameter(i));
157  }
158  }
159 }
160 
161 void MissingMassProb::setParamAngle(const PtEtaPhiMVector& tauvec, int tau, int tautype) {
162  double Pt_tau = tauvec.Pt();
163  int type = tautype;
164  if (tautype > 4 && tautype < 8) type = 4;
165  if (tautype <= 5){
166  if (m_paramVectorAngle.size() > 0){
167  for(int i=0; i<=3; i++){
168  double par = m_paramVectorAngle[i+(type)*4]->Eval(Pt_tau);
169  if (tau==1){
170  m_formulaAngle1->SetParameter(i, par);
171  } else {
172  m_formulaAngle2->SetParameter(i, par);
173  }
174  }
175  }
176  } else {
177  if (m_paramVectorAngleLep.size() > 0){
178  for(int i=0; i<=3; i++){
179  double par = m_paramVectorAngleLep[i]->Eval(Pt_tau);
180  if (tau==1){
181  m_formulaAngle1->SetParameter(i, par);
182  } else {
183  m_formulaAngle2->SetParameter(i, par);
184  }
185  }
186  }
187  }
188 }
189 
190 void MissingMassProb::setParamRatio(int tau, int tautype) {
191  int type = tautype;
192  if(tautype > 4 && tautype < 8) type = 4;
193  if (tautype <= 5){
194  if(tau==1){
196  for(int i=0; i<m_formulaRatio1->GetNpar(); i++){
197  m_formulaRatio1->SetParameter(i, m_paramVectorRatio[5*(tau-1)+type]->GetParameter(i));
198  }
199  } else {
201  for(int i=0; i<m_formulaRatio2->GetNpar(); i++){
202  m_formulaRatio2->SetParameter(i, m_paramVectorRatio[5*(tau-1)+type]->GetParameter(i));
203  }
204  }
205  } else {
206  if(tau==1){
208  for(int i=0; i<m_formulaRatio1->GetNpar(); i++){
209  m_formulaRatio1->SetParameter(i, m_paramVectorRatioLep[0]->GetParameter(i));
210  }
211  } else {
213  for(int i=0; i<m_formulaRatio2->GetNpar(); i++){
214  m_formulaRatio2->SetParameter(i, m_paramVectorRatioLep[0]->GetParameter(i));
215  }
216  }
217  }
218 }
219 
220 // Default Constructor
221 MissingMassProb::MissingMassProb(MMCCalibrationSet::e aset, const std::string& paramFilePath) {
222  m_mmcCalibrationSet = aset;
223  m_allowUseHT = false;
224  m_UseHT = false;
225 
226  m_fParams = NULL;
227  if (!paramFilePath.empty()){
228  std::string total_path = "DiTauMassTools/"+paramFilePath;
229  m_fParams = TFile::Open( (const char*) PathResolverFindCalibFile(total_path).c_str() ,"READ");
230  }
232  m_probListConstant.push_back( std::bind(&mEtAndTauProbabilityWrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5, std::placeholders::_6, std::placeholders::_7) );
233  m_probListOneTau.push_back( std::bind(&dTheta3d_probabilityNewWrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5, std::placeholders::_6, std::placeholders::_7) );
234  m_probListTwoTau.push_back( std::bind(&TauProbabilityNewWrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5, std::placeholders::_6, std::placeholders::_7) );
235  m_probListTwoTau.push_back( std::bind(&MnuProbabilityNewWrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5, std::placeholders::_6, std::placeholders::_7) );
236  if (m_fParams){
238  }
239  }
240  else {
241  m_probListConstant.push_back( std::bind(&mEtAndTauProbabilityWrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5, std::placeholders::_6, std::placeholders::_7) );
242  m_probListOneTau.push_back( std::bind(&dTheta3d_probabilityFastWrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5, std::placeholders::_6, std::placeholders::_7) );
243  m_probListTwoTau.push_back( std::bind(&TauProbabilityWrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5, std::placeholders::_6, std::placeholders::_7) );
244  m_probListTwoTau.push_back( std::bind(&MnuProbabilityWrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5, std::placeholders::_6, std::placeholders::_7) );
245  }
246 
247  if (m_fParams) if (m_fParams->IsOpen()) m_fParams->Close();
248 
250  // [tau_type][parLG][par]
251  // leptonic tau
252  //-par[0]
253  s_fit_param[0][0][0][0]=-9.82013E-1;
254  s_fit_param[0][0][0][1]=9.09874E-1;
255  s_fit_param[0][0][0][2]=0.0;
256  s_fit_param[0][0][0][3]=0.0;
257  //-par[1]
258  s_fit_param[0][0][1][0]=9.96303E1;
259  s_fit_param[0][0][1][1]=1.68873E1;
260  s_fit_param[0][0][1][2]=3.23798E-2;
261  s_fit_param[0][0][1][3]=0.0;
262  //-par[2]
263  s_fit_param[0][0][2][0]=0.0;
264  s_fit_param[0][0][2][1]=0.0;
265  s_fit_param[0][0][2][2]=0.0;
266  s_fit_param[0][0][2][3]=0.3; // fit value is 2.8898E-1, I use 0.3
267  //-par[3]
268  s_fit_param[0][0][3][0]=9.70055E1;
269  s_fit_param[0][0][3][1]=6.46175E1;
270  s_fit_param[0][0][3][2]=3.20679E-2;
271  s_fit_param[0][0][3][3]=0.0;
272  //-par[4]
273  s_fit_param[0][0][4][0]=1.56865;
274  s_fit_param[0][0][4][1]=2.28336E-1;
275  s_fit_param[0][0][4][2]=1.09396E-1;
276  s_fit_param[0][0][4][3]=1.99975E-3;
277  //-par[5]
278  s_fit_param[0][0][5][0]=0.0;
279  s_fit_param[0][0][5][1]=0.0;
280  s_fit_param[0][0][5][2]=0.0;
281  s_fit_param[0][0][5][3]=0.66;
282  //-----------------------------------------------------------------
283  // 1-prong hadronic tau
284  //-par[0]
285  s_fit_param[0][1][0][0]=-2.42674;
286  s_fit_param[0][1][0][1]=7.69124E-1;
287  s_fit_param[0][1][0][2]=0.0;
288  s_fit_param[0][1][0][3]=0.0;
289  //-par[1]
290  s_fit_param[0][1][1][0]=9.52747E1;
291  s_fit_param[0][1][1][1]=1.26319E1;
292  s_fit_param[0][1][1][2]=3.09643E-2;
293  s_fit_param[0][1][1][3]=0.0;
294  //-par[2]
295  s_fit_param[0][1][2][0]=1.71302E1;
296  s_fit_param[0][1][2][1]=3.00455E1;
297  s_fit_param[0][1][2][2]=7.49445E-2;
298  s_fit_param[0][1][2][3]=0.0;
299  //-par[3]
300  s_fit_param[0][1][3][0]=1.06137E2;
301  s_fit_param[0][1][3][1]=6.01548E1;
302  s_fit_param[0][1][3][2]=3.50867E-2;
303  s_fit_param[0][1][3][3]=0.0;
304  //-par[4]
305  s_fit_param[0][1][4][0]=4.26079E-1;
306  s_fit_param[0][1][4][1]=1.76978E-1;
307  s_fit_param[0][1][4][2]=1.43419;
308  s_fit_param[0][1][4][3]=0.0;
309  //-par[5]
310  s_fit_param[0][1][5][0]=0.0;
311  s_fit_param[0][1][5][1]=0.0;
312  s_fit_param[0][1][5][2]=0.0;
313  s_fit_param[0][1][5][3]=0.4;
314  //-----------------------------------------------------------------
315  // 3-prong hadronic tau
316  //-par[0]
317  s_fit_param[0][2][0][0]=-2.43533;
318  s_fit_param[0][2][0][1]=6.12947E-1;
319  s_fit_param[0][2][0][2]=0.0;
320  s_fit_param[0][2][0][3]=0.0;
321  //-par[1]
322  s_fit_param[0][2][1][0]=9.54202;
323  s_fit_param[0][2][1][1]=2.80011E-1;
324  s_fit_param[0][2][1][2]=2.49782E-1;
325  s_fit_param[0][2][1][3]=0.0;
326  //-par[2]
327  s_fit_param[0][2][2][0]=1.61325E1;
328  s_fit_param[0][2][2][1]=1.74892E1;
329  s_fit_param[0][2][2][2]=7.05797E-2;
330  s_fit_param[0][2][2][3]=0.0;
331  //-par[3]
332  s_fit_param[0][2][3][0]=1.17093E2;
333  s_fit_param[0][2][3][1]=4.70000E1;
334  s_fit_param[0][2][3][2]=3.87085E-2;
335  s_fit_param[0][2][3][3]=0.0;
336  //-par[4]
337  s_fit_param[0][2][4][0]=4.16557E-1;
338  s_fit_param[0][2][4][1]=1.58902E-1;
339  s_fit_param[0][2][4][2]=1.53516;
340  s_fit_param[0][2][4][3]=0.0;
341  //-par[5]
342  s_fit_param[0][2][5][0]=0.0;
343  s_fit_param[0][2][5][1]=0.0;
344  s_fit_param[0][2][5][2]=0.0;
345  s_fit_param[0][2][5][3]=0.95;
346 
347 
349  // [tau_type][parLG][par]
350  // leptonic tau
351  //-par[0]
352  s_fit_param[1][0][0][0]=-9.82013E-1;
353  s_fit_param[1][0][0][1]=9.09874E-1;
354  s_fit_param[1][0][0][2]=0.0;
355  s_fit_param[1][0][0][3]=0.0;
356  s_fit_param[1][0][0][4]=0.0;
357  //-par[1]
358  s_fit_param[1][0][1][0]=9.96303E1;
359  s_fit_param[1][0][1][1]=1.68873E1;
360  s_fit_param[1][0][1][2]=3.23798E-2;
361  s_fit_param[1][0][1][3]=0.0;
362  s_fit_param[1][0][1][4]=0.0;
363  //-par[2]
364  s_fit_param[1][0][2][0]=0.0;
365  s_fit_param[1][0][2][1]=0.0;
366  s_fit_param[1][0][2][2]=0.0;
367  s_fit_param[1][0][2][3]=0.3; // fit value is 2.8898E-1, I use 0.3
368  s_fit_param[1][0][2][4]=0.0;
369  //-par[3]
370  s_fit_param[1][0][3][0]=9.70055E1;
371  s_fit_param[1][0][3][1]=6.46175E1;
372  s_fit_param[1][0][3][2]=3.20679E-2;
373  s_fit_param[1][0][3][3]=0.0;
374  s_fit_param[1][0][3][4]=0.0;
375  //-par[4]
376  s_fit_param[1][0][4][0]=1.56865;
377  s_fit_param[1][0][4][1]=2.28336E-1;
378  s_fit_param[1][0][4][2]=1.09396E-1;
379  s_fit_param[1][0][4][3]=1.99975E-3;
380  s_fit_param[1][0][4][4]=0.0;
381  //-par[5]
382  s_fit_param[1][0][5][0]=0.0;
383  s_fit_param[1][0][5][1]=0.0;
384  s_fit_param[1][0][5][2]=0.0;
385  s_fit_param[1][0][5][3]=0.66;
386  s_fit_param[1][0][5][4]=0.0;
387  //-----------------------------------------------------------------
388  //-----------------------
389  // Only hadronic tau's were re-parametrized in MC12. The parameterization now
390  // goes to P(tau)=1 TeV.
391  //-----------------------------------------------------------------
392  // 1-prong hadronic tau
393  //---par[0]
394  s_fit_param[1][1][0][0]=0.7568;
395  s_fit_param[1][1][0][1]=-0.0001469;
396  s_fit_param[1][1][0][2]=5.413E-7;
397  s_fit_param[1][1][0][3]=-6.754E-10;
398  s_fit_param[1][1][0][4]=2.269E-13;
399  //---par[1]
400  s_fit_param[1][1][1][0]=-0.0288208;
401  s_fit_param[1][1][1][1]=0.134174;
402  s_fit_param[1][1][1][2]=-142.588;
403  s_fit_param[1][1][1][3]=-0.00035606;
404  s_fit_param[1][1][1][4]=-6.94567E-20;
405  //---par[2]
406  s_fit_param[1][1][2][0]=-0.00468927;
407  s_fit_param[1][1][2][1]=0.0378737;
408  s_fit_param[1][1][2][2]=-260.284;
409  s_fit_param[1][1][2][3]=0.00241158;
410  s_fit_param[1][1][2][4]=-6.01766E-7;
411  //---par[3]
412  s_fit_param[1][1][3][0]=-0.170424;
413  s_fit_param[1][1][3][1]=0.135764;
414  s_fit_param[1][1][3][2]=-50.2361;
415  s_fit_param[1][1][3][3]=0.00735544;
416  s_fit_param[1][1][3][4]=-7.34073E-6;
417  //---par[4]
418  s_fit_param[1][1][4][0]=-0.0081364;
419  s_fit_param[1][1][4][1]=0.0391428;
420  s_fit_param[1][1][4][2]=-141.936;
421  s_fit_param[1][1][4][3]=0.0035034;
422  s_fit_param[1][1][4][4]=-1.21956E-6;
423  //-par[5]
424  s_fit_param[1][1][5][0]=0.0;
425  s_fit_param[1][1][5][1]=0.0;
426  s_fit_param[1][1][5][2]=0.0;
427  s_fit_param[1][1][5][3]=0.624*0.00125; // multiplied by a bin size
428  s_fit_param[1][1][5][4]=0.0;
429  //-----------------------------------------------------------------
430  // 3-prong hadronic tau
431  //---par[0]
432  s_fit_param[1][2][0][0]=0.7562;
433  s_fit_param[1][2][0][1]=-1.168E-5;
434  s_fit_param[1][2][0][2]=0.0;
435  s_fit_param[1][2][0][3]=0.0;
436  s_fit_param[1][2][0][4]=0.0;
437  //---par[1]
438  s_fit_param[1][2][1][0]=-0.0420458;
439  s_fit_param[1][2][1][1]=0.15917;
440  s_fit_param[1][2][1][2]=-80.3259;
441  s_fit_param[1][2][1][3]=0.000125729;
442  s_fit_param[1][2][1][4]=-2.43945E-18;
443  //---par[2]
444  s_fit_param[1][2][2][0]=-0.0216898;
445  s_fit_param[1][2][2][1]=0.0243497;
446  s_fit_param[1][2][2][2]=-63.8273;
447  s_fit_param[1][2][2][3]=0.0148339;
448  s_fit_param[1][2][2][4]=-4.45351E-6;
449  //---par[3]
450  s_fit_param[1][2][3][0]=-0.0879411;
451  s_fit_param[1][2][3][1]=0.110092;
452  s_fit_param[1][2][3][2]=-75.4901;
453  s_fit_param[1][2][3][3]=0.0116915;
454  s_fit_param[1][2][3][4]=-1E-5;
455  //---par[4]
456  s_fit_param[1][2][4][0]=-0.0118324;
457  s_fit_param[1][2][4][1]=0.0280538;
458  s_fit_param[1][2][4][2]=-85.127;
459  s_fit_param[1][2][4][3]=0.00724948;
460  s_fit_param[1][2][4][4]=-2.38792E-6;
461  //-par[5]
462  s_fit_param[1][2][5][0]=0.0;
463  s_fit_param[1][2][5][1]=0.0;
464  s_fit_param[1][2][5][2]=0.0;
465  s_fit_param[1][2][5][3]=0.6167*0.00125; // multiplied by a bin size
466  s_fit_param[1][2][5][4]=0.0;
467 
468  // TER parameterization, for now based on p1130, to be checked and updated with new tag
469  // [tau_type][eta_bin][parameter]
470  // for 1-prongs
471  s_ter_sigma_par[0][0][0]=0.311717;
472  s_ter_sigma_par[0][0][1]=0.0221615;
473  s_ter_sigma_par[0][0][2]=0.859698;
474  s_ter_sigma_par[0][1][0]=0.290019;
475  s_ter_sigma_par[0][1][1]=0.0225794;
476  s_ter_sigma_par[0][1][2]=0.883407;
477  s_ter_sigma_par[0][2][0]=0.352312;
478  s_ter_sigma_par[0][2][1]=0.0196381;
479  s_ter_sigma_par[0][2][2]=0.629708;
480  s_ter_sigma_par[0][3][0]=0.342059;
481  s_ter_sigma_par[0][3][1]=0.0275107;
482  s_ter_sigma_par[0][3][2]=0.48065;
483  s_ter_sigma_par[0][4][0]=0.481564;
484  s_ter_sigma_par[0][4][1]=0.0197219;
485  s_ter_sigma_par[0][4][2]=0.0571714;
486  s_ter_sigma_par[0][5][0]=0.41264;
487  s_ter_sigma_par[0][5][1]=0.0233964;
488  s_ter_sigma_par[0][5][2]=0.515674;
489  s_ter_sigma_par[0][6][0]=0.20112;
490  s_ter_sigma_par[0][6][1]=0.0339914;
491  s_ter_sigma_par[0][6][2]=0.944524;
492  s_ter_sigma_par[0][7][0]=0.0892094;
493  s_ter_sigma_par[0][7][1]=0.0210225;
494  s_ter_sigma_par[0][7][2]=1.34014;
495  s_ter_sigma_par[0][8][0]=0.175554;
496  s_ter_sigma_par[0][8][1]=0.0210968;
497  s_ter_sigma_par[0][8][2]=0.813925;
498  s_ter_sigma_par[0][9][0]=0.0;
499  s_ter_sigma_par[0][9][1]=0.0340279;
500  s_ter_sigma_par[0][9][2]=1.30856;
501  // for 3-prongs
502  s_ter_sigma_par[1][0][0]=0.303356;
503  s_ter_sigma_par[1][0][1]=0.0299807;
504  s_ter_sigma_par[1][0][2]=1.25388;
505  s_ter_sigma_par[1][1][0]=0.358106;
506  s_ter_sigma_par[1][1][1]=0.0229604;
507  s_ter_sigma_par[1][1][2]=1.02222;
508  s_ter_sigma_par[1][2][0]=0.328643;
509  s_ter_sigma_par[1][2][1]=0.025684;
510  s_ter_sigma_par[1][2][2]=1.02594;
511  s_ter_sigma_par[1][3][0]=0.497332;
512  s_ter_sigma_par[1][3][1]=0.0215113;
513  s_ter_sigma_par[1][3][2]=0.30055;
514  s_ter_sigma_par[1][4][0]=0.4493;
515  s_ter_sigma_par[1][4][1]=0.0280311;
516  s_ter_sigma_par[1][4][2]=0.285793;
517  s_ter_sigma_par[1][5][0]=0.427811;
518  s_ter_sigma_par[1][5][1]=0.0316536;
519  s_ter_sigma_par[1][5][2]=0.457286;
520  s_ter_sigma_par[1][6][0]=0.165288;
521  s_ter_sigma_par[1][6][1]=0.0376361;
522  s_ter_sigma_par[1][6][2]=1.3913;
523  s_ter_sigma_par[1][7][0]=0.289798;
524  s_ter_sigma_par[1][7][1]=0.0140801;
525  s_ter_sigma_par[1][7][2]=0.83603;
526  s_ter_sigma_par[1][8][0]=0.186823;
527  s_ter_sigma_par[1][8][1]=0.0213053;
528  s_ter_sigma_par[1][8][2]=0.968934;
529  s_ter_sigma_par[1][9][0]=0.301673;
530  s_ter_sigma_par[1][9][1]=0.0145606;
531  s_ter_sigma_par[1][9][2]=0.514022;
532 
533 
534 }
535 
536 // Default Destructor
538 }
539 
540 double MissingMassProb::apply(MissingMassInput& preparedInput, const int & tau_type1, const int & tau_type2, const PtEtaPhiMVector & tauvec1, const PtEtaPhiMVector & tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector & nuvec2, bool constant, bool oneTau, bool twoTau) {
541  double prob = 1.;
542  if (constant == true) {
543  for (auto& f: m_probListConstant) {
544  prob*=f(preparedInput, tau_type1, tau_type2, tauvec1, tauvec2, nuvec1, nuvec2);
545  }
546  } else if (oneTau == true) {
547  for (auto& f: m_probListOneTau) {
548  prob*=f(preparedInput, tau_type1, tau_type2, tauvec1, tauvec2, nuvec1, nuvec2);
549  }
550  } else if (twoTau == true) {
551  for (auto& f: m_probListTwoTau) {
552  prob*=f(preparedInput, tau_type1, tau_type2, tauvec1, tauvec2, nuvec1, nuvec2);
553  }
554  }
555  return prob;
556 }
557 
558 
559 double MissingMassProb::MetProbability(MissingMassInput& preparedInput, const double & met1,const double & met2,const double & MetSigma1,const double & MetSigma2) {
560 
561 
562  double metprob;
563  if(MetSigma1>1.0 && MetSigma2>1.0) // it doesn't make sense if MET resolution sigma is <1 GeV
564  {
565  //SpeedUp
566  metprob=exp(-0.5*(met1*met1/(MetSigma1*MetSigma1)+met2*met2/(MetSigma2*MetSigma2)))/(MetSigma1*MetSigma2*2*TMath::Pi());
567  }
568  else
569  {
570  if(preparedInput.m_fUseVerbose==1) Warning("DiTauMassTools", "MissingMassCalculator::MetProbability: either MetSigma1 or MetSigma2 are <1 GeV--- too low, returning prob=1");
571  metprob=1.;
572  }
573 
574 
575  return metprob;
576 
577 
578 }
579 
581 {
582  double proba=1.;
583  double metprob;
584 
585  // deltaMEt is the difference between the neutrino sum and the MEt (or -HT if useHT),
586  //corrected possibly from tau E scanning
587 
588  double deltaMetX=preparedInput.m_metVec.X()-preparedInput.m_inputMEtX;
589  double deltaMetY=preparedInput.m_metVec.Y()-preparedInput.m_inputMEtY;
590 
591  // possibly correct for subtract tau scanning here
592  // const double met_smearL=(deltaMetVec.X()*m_metCovPhiCos+deltaMetVec.Y()*m_metCovPhiSin;
593  // const double met_smearP=-deltaMetVec.X()*m_metCovPhiSin+deltaMetVec.Y()*m_metCovPhiCos;
594 
595  // rotate into error ellipse axis
596  const double met_smearL=deltaMetX*cos(preparedInput.m_METcovphi)+deltaMetY*sin(preparedInput.m_METcovphi);
597  const double met_smearP=-deltaMetX*sin(preparedInput.m_METcovphi)+deltaMetY*cos(preparedInput.m_METcovphi);
598 
599 
600  if (m_UseHT)
601  {
602  if ( preparedInput.m_tauTypes==TauTypes::hh ) {
603 
604  metprob=MHtProbabilityHH(preparedInput, met_smearL,met_smearP,preparedInput.m_inputMEtT,preparedInput.m_MEtT,preparedInput.m_htOffset); // for had-had
605  }
606  else {
607  metprob=MHtProbability(preparedInput, met_smearL,met_smearP,preparedInput.m_inputMEtT,preparedInput.m_MEtT,preparedInput.m_htOffset); // for lep-had Winter 2012 analysis
608  }
609  }
610  else {
611  metprob=MetProbability(preparedInput, met_smearL,met_smearP,preparedInput.m_METsigmaL,preparedInput.m_METsigmaP);
612  }
613 
614  proba=metprob;
615 
616  return proba;
617 }
618 
619 //------------------- simple TauProbability for LFV
620 double MissingMassProb::TauProbabilityLFV(MissingMassInput& preparedInput, const int & type1, const PtEtaPhiMVector & vis1, const PtEtaPhiMVector & nu1)
621 {
622  double prob=1.0;
623  if(m_fUseTauProbability==0) return prob; // don't apply TauProbability
624  double prob1=1.0;
625  const double mtau=ParticleConstants::tauMassInMeV / GEV;
626  const double R1=nu1.E()/vis1.E();
627  //--- dealing with 1st tau
628  double m1=nu1.M();
629  double m2=vis1.M();
630  double E1=0.5*(mtau*mtau+m1*m1-m2*m2)/mtau;
631  double E2=mtau-E1;
632  if(E1<=m1 || E1>=mtau)
633  {
634  if(preparedInput.m_fUseVerbose==1) Warning("DiTauMassTools", "Warning in MissingMassCalculator::TauProbability: bad E1, returning 0 ");
635  return 0.0;
636  }
637  if(E2<=m2 || E2>=mtau)
638  {
639  if(preparedInput.m_fUseVerbose==1) Warning("DiTauMassTools", "Warning in MissingMassCalculator::TauProbability: bad E2, returning 0 ");
640  return 0.0;
641  }
642  preparedInput.m_tlv_tmp.SetPxPyPzE(0.,0.,0.,0.);
643  preparedInput.m_tlv_tmp+=nu1;
644  preparedInput.m_tlv_tmp+=vis1;
645  // double p=(nu1+vis1).P();
646  double p=preparedInput.m_tlv_tmp.P();
647  double V=p/sqrt(p*p+mtau*mtau);
648  double p0;
649  if(type1==8) p0=sqrt(E2*E2-m2*m2); // leptonic tau
650  else p0=E1; // hadronic tau
651  prob1=0.5*mtau/(p0*V*pow(R1+1.0,2));
652  // avoid too large values
653  prob=std::min(prob1,1.);
654  return prob;
655 }
656 
657 double MissingMassProb::TauProbability(MissingMassInput& preparedInput, const int & type1, const PtEtaPhiMVector & vis1, const PtEtaPhiMVector & nu1,
658  const int & type2, const PtEtaPhiMVector & vis2, const PtEtaPhiMVector & nu2)
659 {
660  double prob=1.0;
661  if(m_fUseTauProbability==0) return prob; // don't apply TauProbability
662  double prob1=1.0;
663  double prob2=1.0;
664  const double mtau=ParticleConstants::tauMassInMeV / GEV;
665  const double R1=nu1.E()/vis1.E();
666  const double R2=nu2.E()/vis2.E();
667  //--- dealing with 1st tau
668  double m1=nu1.M();
669  double m2=vis1.M();
670  double E1=0.5*(mtau*mtau+m1*m1-m2*m2)/mtau;
671  double E2=mtau-E1;
672  if(E1<=m1 || E1>=mtau)
673  {
674  if(preparedInput.m_fUseVerbose==1) Warning("DiTauMassTools", "Warning in MissingMassCalculator::TauProbability: bad E1, returning 0 ");
675  return 0.0;
676  }
677  if(E2<=m2 || E2>=mtau)
678  {
679  if(preparedInput.m_fUseVerbose==1) Warning("DiTauMassTools", "Warning in MissingMassCalculator::TauProbability: bad E2, returning 0 ");
680  return 0.0;
681  }
682  preparedInput.m_tlv_tmp.SetPxPyPzE(0.,0.,0.,0.);
683  preparedInput.m_tlv_tmp+=nu1;
684  preparedInput.m_tlv_tmp+=vis1;
685  // double p=(nu1+vis1).P();
686  double p=preparedInput.m_tlv_tmp.P();
687  double V=p/sqrt(p*p+mtau*mtau);
688  double p0;
689  if(type1==8) p0=sqrt(E2*E2-m2*m2); // leptonic tau
690  else p0=E1; // hadronic tau
691  prob1=0.5*mtau/(p0*V*pow(R1+1.0,2));
692  // avoid too large values
693  prob1=std::min(prob1,1.);
694 
695 
696  //--- dealing with 2nd tau
697  m1=nu2.M();
698  m2=vis2.M();
699  E1=0.5*(mtau*mtau+m1*m1-m2*m2)/mtau;
700  E2=mtau-E1;
701  if(E1<=m1 || E1>=mtau)
702  {
703  if(preparedInput.m_fUseVerbose==1) Warning("DiTauMassTools", "Warning in MissingMassCalculator::TauProbability: bad E1, returning 0 ");
704  return 0.0;
705  }
706  if(E2<=m2 || E2>=mtau)
707  {
708  if(preparedInput.m_fUseVerbose==1) Warning("DiTauMassTools", "Warning in MissingMassCalculator::TauProbability: bad E2, returning 0 ");
709  return 0.0;
710  }
711  preparedInput.m_tlv_tmp.SetPxPyPzE(0.,0.,0.,0.);
712  preparedInput.m_tlv_tmp+=nu2;
713  preparedInput.m_tlv_tmp+=vis2;
714  // p=(nu2+vis2).P();
715  p=preparedInput.m_tlv_tmp.P();
716 
717 
718  V=p/sqrt(p*p+mtau*mtau);
719  if(type2==8) p0=sqrt(E2*E2-m2*m2); // leptonic tau
720  else p0=E1; // hadronic tau
721  prob2=0.5*mtau/(p0*V*pow(R2+1.0,2));
722  // avoid too large values
723  prob2=std::min(prob2,1.);
724  prob=prob1*prob2;
725  return prob;
726 }
727 
728 
729 // --------- Updated version of TauProbability for lep-had events with Njet25=0, takes into account Winter-2012 analysis cuts
730 double MissingMassProb::TauProbability(MissingMassInput& preparedInput, const int & type1, const PtEtaPhiMVector & vis1, const PtEtaPhiMVector & nu1,
731  const int & type2, const PtEtaPhiMVector & vis2, const PtEtaPhiMVector & nu2, const double & detmet) {
732  double prob=1.0;
733 
734  if(m_fUseTauProbability==0) return prob; // don't apply TauProbability
735 
736  if(m_UseHT)
737  {
738  if(detmet<20.0) // low MET Njet=0 category
739  {
740  const double R1=nu1.P()/vis1.P();
741  const double R2=nu2.P()/vis2.P();
742  const double lep_p1[4]={0.417,0.64,0.52,0.678};
743  const double lep_p2[4]={0.23,0.17,0.315,0.319};
744  const double lep_p3[4]={0.18,0.33,0.41,0.299};
745  const double lep_p4[4]={0.033,0.109,0.129,0.096};
746  const double lep_p5[4]={0.145,0.107,0.259,0.304};
747  int ind=3;
748  int indT=3;
749  const double n_1pr[4]={-0.15,-0.13,-0.25,-0.114};
750  const double s_1pr[4]={0.40,0.54,0.62,0.57};
751  const double n_3pr[4]={-1.08,-1.57,-0.46,-0.39};
752  const double s_3pr[4]={0.53,0.85,0.61,0.53};
753  double Ptau=0.0;
754  double Plep=0.0;
755  if(type1>=0 && type1<=5)
756  {
757  Ptau=(nu1+vis1).P();
758  Plep=(nu2+vis2).P();
759  }
760  if(type2>=0 && type2<=5)
761  {
762  Ptau=(nu2+vis2).P();
763  Plep=(nu1+vis1).P();
764  }
765  if(Plep<50.0 && Plep>=45.0) ind=2;
766  if(Plep<45.0 && Plep>=40.0) ind=1;
767  if(Plep<40.0) ind=0;
768  if(Ptau<50.0 && Ptau>=45.0) indT=2;
769  if(Ptau<45.0 && Ptau>=40.0) indT=1;
770  if(Ptau<40.0) indT=0;
771  if(type1==8) prob=prob*(lep_p5[ind]*TMath::Gaus(R1,lep_p1[ind],lep_p2[ind])+TMath::Landau(R1,lep_p3[ind],lep_p4[ind]))/(1+lep_p5[ind]);
772  if(type2==8) prob=prob*(lep_p5[ind]*TMath::Gaus(R2,lep_p1[ind],lep_p2[ind])+TMath::Landau(R2,lep_p3[ind],lep_p4[ind]))/(1+lep_p5[ind]);
773 
774  if(type1>=0 && type1<=2) prob=prob*TMath::Gaus(R1,n_1pr[indT],s_1pr[indT]);
775  if(type2>=0 && type2<=2) prob=prob*TMath::Gaus(R2,n_1pr[indT],s_1pr[indT]);
776  if(type1>=3 && type1<=5) prob=prob*TMath::Gaus(R1,n_3pr[indT],s_3pr[indT]);
777  if(type2>=3 && type2<=5) prob=prob*TMath::Gaus(R2,n_3pr[indT],s_3pr[indT]);
778 
779  }
780  else // high MET Njet=0 category
781  {
782  const double R1=nu1.P()/vis1.P();
783  const double R2=nu2.P()/vis2.P();
784  const double lep_p1[4]={0.441,0.64,0.79,0.8692};
785  const double lep_p2[4]={0.218,0.28,0.29,0.3304};
786  const double lep_p3[4]={0.256,0.33,0.395,0.4105};
787  const double lep_p4[4]={0.048,0.072,0.148,0.1335};
788  const double lep_p5[4]={0.25,0.68,0.10,0.2872};
789  int ind=3;
790  const double p_1prong=-3.706;
791  const double p_3prong=-5.845;
792  double Ptau=0.0;
793  double Plep=0.0;
794  if(type1>=0 && type1<=5)
795  {
796  Ptau=(nu1+vis1).P();
797  Plep=(nu2+vis2).P();
798  }
799  if(type2>=0 && type2<=5)
800  {
801  Ptau=(nu2+vis2).P();
802  Plep=(nu1+vis1).P();
803  }
804  if(Plep<50.0 && Plep>=45.0) ind=2;
805  if(Plep<45.0 && Plep>=40.0) ind=1;
806  if(Plep<40.0) ind=0;
807  const double scale1prong=Ptau>45.0 ? 1.0 : -1.019/((Ptau*0.0074-0.059)*p_1prong);
808  const double scale3prong=Ptau>40.0 ? 1.0 : -1.24/((Ptau*0.0062-0.033)*p_3prong);
809  if(type1==8) prob=prob*(lep_p5[ind]*TMath::Gaus(R1,lep_p1[ind],lep_p2[ind])+TMath::Landau(R1,lep_p3[ind],lep_p4[ind]))/(1+lep_p5[ind]);
810  if(type2==8) prob=prob*(lep_p5[ind]*TMath::Gaus(R2,lep_p1[ind],lep_p2[ind])+TMath::Landau(R2,lep_p3[ind],lep_p4[ind]))/(1+lep_p5[ind]);
811 
812  if(type1>=0 && type1<=2) prob=prob*exp(p_1prong*R1*scale1prong)*std::abs(p_1prong*scale1prong)*0.02; // introduced normalization to account for sharpening of probability at low E(tau)
813  if(type2>=0 && type2<=2) prob=prob*exp(p_1prong*R2*scale1prong)*std::abs(p_1prong*scale1prong)*0.02;
814  if(type1>=3 && type1<=5) prob=prob*exp(p_3prong*R1*scale3prong)*std::abs(p_3prong*scale3prong)*0.02;
815  if(type2>=3 && type2<=5) prob=prob*exp(p_3prong*R2*scale3prong)*std::abs(p_3prong*scale3prong)*0.02;
816  }
817  }
818  //----------------- had-had channel ---------------------------------------
819  if( preparedInput.m_tauTypes==TauTypes::hh )
820  {
821 
822  if(m_UseHT) // Events with no jets
823  {
824 
825  const double R[2]={nu1.P()/vis1.P(),nu2.P()/vis2.P()};
826  const double E[2]={(nu1+vis1).E(),(nu2+vis2).E()};
827  const int tau_type[2]={type1,type2};
828  int order1= vis1.Pt()>vis2.Pt() ? 0 : 1;
829  int order2= vis1.Pt()>vis2.Pt() ? 1 : 0;
830 
831  double par_1p[2][6]; // P(tau)-scaling; lead, sub-lead
832  double par_3p[2][6]; // P(tau)-scaling; lead, sub-lead
833 
834  par_1p[0][0]=0.1273; par_1p[0][1]=-0.2479; par_1p[0][2]=1.0; par_1p[0][3]=-43.16; par_1p[0][4]=0.0; par_1p[0][5]=0.0;
835  par_1p[1][0]=0.3736; par_1p[1][1]=-1.441; par_1p[1][2]=1.0; par_1p[1][3]=-29.82; par_1p[1][4]=0.0; par_1p[1][5]=0.0;
836  par_3p[0][0]=0.1167; par_3p[0][1]=-0.1388; par_3p[0][2]=1.0; par_3p[0][3]=-44.77; par_3p[0][4]=0.0; par_3p[0][5]=0.0;
837  par_3p[1][0]=0.3056; par_3p[1][1]=-2.18; par_3p[1][2]=1.0; par_3p[1][3]=-19.09; par_3p[1][4]=0.0; par_3p[1][5]=0.0;
838  // parameters for sub-leading tau
839  const double C1p=0.062;
840  const double C3p=0.052;
841  const double G1p=1.055;
842  const double G3p=1.093;
843  // Probability for leading tau
844 
845  if( tau_type[order1]>=0 && tau_type[order1]<=2 ) // 1-prong
846  {
847  //double x=std::min(300.,std::max(E[order1],45.0));
848  // 50 to be safe. TO be finalised.
849  // double x=std::min(300.,std::max(E[order1],50.0));
850  double x=std::min(E[order1],300.0);
851  const double slope=par_1p[0][0]+par_1p[0][1]/(par_1p[0][2]*x+par_1p[0][3])+par_1p[0][4]*x > 0.01 ?
852  par_1p[0][0]+par_1p[0][1]/(par_1p[0][2]*x+par_1p[0][3])+par_1p[0][4]*x : 0.01;
853  prob=prob*exp(-R[order1]/slope)*0.04/std::abs(slope);
854  }
855  if( tau_type[order1]>=3 && tau_type[order1]<=5 ) // 3-prong
856  {
857  //double x=std::min(300.,std::max(E[order1],45.0));
858  // double x=std::min(300.,std::max(E[order1],50.0));
859  double x=std::min(E[order1],300.0);
860  const double slope=par_3p[0][0]+par_3p[0][1]/(par_3p[0][2]*x+par_3p[0][3])+par_3p[0][4]*x > 0.01 ?
861  par_3p[0][0]+par_3p[0][1]/(par_3p[0][2]*x+par_3p[0][3])+par_3p[0][4]*x : 0.01;
862  prob=prob*exp(-R[order1]/slope)*0.04/std::abs(slope);
863  }
864  // Probability for sub-leading tau
865  if( tau_type[order2]>=0 && tau_type[order2]<=2 ) // 1-prong
866  {
867  const double par[4]={0.1147,-0.09675,-35.0,3.711E-11};
868  double x=std::min(300.,std::max(E[order2],30.0));
869  // double x=std::min(300.,std::max(E[order2],50.0));
870  const double sigma=G1p*(par_1p[1][0]+par_1p[1][1]/(par_1p[1][2]*x+par_1p[1][3])+par_1p[1][4]*x+par_1p[1][5]*x*x) > 0.01 ?
871  G1p*(par_1p[1][0]+par_1p[1][1]/(par_1p[1][2]*x+par_1p[1][3])+par_1p[1][4]*x+par_1p[1][5]*x*x) : 0.01;
872  if(x<36.0) x=36.0;
873  const double mean=par[0]+par[1]/(x+par[2])+par[3]*pow(x,4);
874  prob=prob*C1p*TMath::Gaus(R[order2],mean,sigma);
875  }
876  if( tau_type[order2]>=3 && tau_type[order2]<=5 ) // 3-prong
877  {
878  double x=std::min(300.,std::max(E[order2],20.0));
879  // double x=std::min(300.,std::max(E[order2],50.0));
880  const double sigma=G3p*(par_3p[1][0]+par_3p[1][1]/(par_3p[1][2]*x+par_3p[1][3])+par_3p[1][4]*x+par_3p[1][5]*x*x) > 0.01 ?
881  G3p*(par_3p[1][0]+par_3p[1][1]/(par_3p[1][2]*x+par_3p[1][3])+par_3p[1][4]*x+par_3p[1][5]*x*x) : 0.01;
882  const double par[4]={0.2302,-2.012,-36.08,-0.000373};
883  if(x<37.0) x=37.0;
884  const double mean=par[0]+par[1]/(x+par[2])+par[3]*x;
885  prob=prob*C3p*TMath::Gaus(R[order2],mean,sigma);
886  }
887  }
888  if(!m_UseHT) // Events with jets
889  {
890  const double R1=nu1.P()/vis1.P();
891  const double R2=nu2.P()/vis2.P();
892  const double E1=(nu1+vis1).E();
893  const double E2=(nu2+vis2).E();
894  int order1= vis1.Pt()>vis2.Pt() ? 0 : 1;
895  int order2= vis1.Pt()>vis2.Pt() ? 1 : 0;
896  const double slope_1p[2]={-3.185,-2.052}; // lead, sub-lead
897  const double slope_3p[2]={-3.876,-2.853}; // lead, sub-lead
898  double par_1p[2][3]; // scaling below 100 GeV; lead, sub-lead
899  double par_3p[2][3]; // scaling below 100 GeV; lead, sub-lead
900  par_1p[0][0]=-0.3745; par_1p[0][1]=0.01417; par_1p[0][2]=-7.285E-5; // [0][i] is always adjusted to match slope at 100 GeV
901  par_1p[1][0]=-0.3683; par_1p[1][1]=0.01807; par_1p[1][2]=-9.514E-5;
902  par_3p[0][0]=-0.3055; par_3p[0][1]=0.01149; par_3p[0][2]=-5.855E-5;
903  par_3p[1][0]=-0.3410; par_3p[1][1]=0.01638; par_3p[1][2]=-9.465E-5;
904  double scale1;
905  double scale2;
906  if(type1>=0 && type1<=2) // 1-prong
907  {
908  scale1=E1>100.0 ? 1.0 : 1.0/std::abs((par_1p[order1][0]+par_1p[order1][1]*E1+par_1p[order1][2]*E1*E1)*slope_1p[order1]);
909  if(scale1<1.0) scale1=1.0;
910  if(scale1>100.0) scale1=100.0;
911  prob=prob*std::abs(slope_1p[order1])*scale1*exp(slope_1p[order1]*scale1*R1)*0.04;
912  }
913  if(type1>=3 && type1<=5) // 3-prong
914  {
915  scale1=E1>100.0 ? 1.0 : 1.0/std::abs((par_3p[order1][0]+par_3p[order1][1]*E1+par_3p[order1][2]*E1*E1)*slope_3p[order1]);
916  if(scale1<1.0) scale1=1.0;
917  if(scale1>100.0) scale1=100.0;
918  prob=prob*std::abs(slope_3p[order1])*scale1*exp(slope_3p[order1]*scale1*R1)*0.04;
919  }
920  if(type2>=0 && type2<=2) // 1-prong
921  {
922  scale2=E2>100.0 ? 1.0 : 1.0/std::abs((par_1p[order2][0]+par_1p[order2][1]*E2+par_1p[order2][2]*E2*E2)*slope_1p[order2]);
923  if(scale2<1.0) scale2=1.0;
924  if(scale2>100.0) scale2=100.0;
925  prob=prob*std::abs(slope_1p[order2])*scale2*exp(slope_1p[order2]*scale2*R2)*0.04;
926  }
927  if(type2>=3 && type2<=5) // 3-prong
928  {
929  scale2=E2>100.0 ? 1.0 : 1.0/std::abs((par_3p[order2][0]+par_3p[order2][1]*E2+par_3p[order2][2]*E2*E2)*slope_3p[order2]);
930  if(scale2<1.0) scale2=1.0;
931  if(scale2>100.0) scale2=100.0;
932  prob=prob*std::abs(slope_3p[order2])*scale2*exp(slope_3p[order2]*scale2*R2)*0.04;
933 
934  }
935  }
936 
937  }
938  // prob=std::min(prob,1.); // Sasha commented out this line, it was introduced by David. Have to ask about its purpose.
939 
940  return prob;
941 }
942 
943 
944 // returns Mnu probability according pol6 parameterization
945 double MissingMassProb::MnuProbability(MissingMassInput& preparedInput, double mnu, double binsize)
946 {
947  double prob=1.0;
948  double norm=4851900.0;
949  double p[7];
950  p[0]=-288.6/norm; p[1]=6.217E4/(2.0*norm); p[2]=2.122E4/(3.0*norm); p[3]=-9.067E4/(4.0*norm);
951  p[4]=1.433E5/(5.0*norm); p[5]=-1.229E5/(6.0*norm); p[6]=3.434E4/(7.0*norm);
952  double int1=0.0;
953  double int2=0.0;
954  double x1= mnu+0.5*binsize < 1.777-0.113 ? mnu+0.5*binsize : 1.777-0.113;
955  double x2= mnu-0.5*binsize > 0.0 ? mnu-0.5*binsize : 0.0;
956  for(int i=0; i<7; i++)
957  {
958  int1=p[i]*pow(x1,i+1)+int1;
959  int2=p[i]*pow(x2,i+1)+int2;
960  }
961  prob=int1-int2;
962  if(prob<0.0)
963  {
964  if(preparedInput.m_fUseVerbose==1) Warning("DiTauMassTools", "Warning in MissingMassCalculator::MnuProbability: negative probability!!! ");
965  return 0.0;
966  }
967  if(prob>1.0)
968  {
969  if(preparedInput.m_fUseVerbose==1) Warning("DiTauMassTools", "Warning in MissingMassCalculator::MnuProbability: probability > 1!!! ");
970  return 1.0;
971  }
972  return prob;
973 }
974 
975 // returns Mnu probability according pol6 parameterization
976 double MissingMassProb::MnuProbability(MissingMassInput& preparedInput, double mnu)
977 {
978  if(m_fUseMnuProbability==0) return 1.0;
979  double prob=1.0;
980  const double norm=4851900.0;
981  double p[7];
982  p[0]=-288.6; p[1]=6.217E4; p[2]=2.122E4; p[3]=-9.067E4;
983  p[4]=1.433E5; p[5]=-1.229E5; p[6]=3.434E4;
984  double int1=0.0;
985  for(int i=0; i<7; i++)
986  {
987  int1+=p[i]*pow(mnu,i);
988  }
989  prob=int1/norm;
990  if(prob<0.0)
991  {
992  if(preparedInput.m_fUseVerbose==1) Warning("DiTauMassTools", "Warning in MissingMassCalculator::MnuProbability: negative probability!!! ");
993  return 0.0;
994  }
995  if(prob>1.0)
996  {
997  if(preparedInput.m_fUseVerbose==1) Warning("DiTauMassTools", "Warning in MissingMassCalculator::MnuProbability: probability > 1!!! ");
998  return 1.0;
999  }
1000  return prob;
1001 }
1002 
1003 double MissingMassProb::MHtProbability(MissingMassInput& preparedInput, const double & d_mhtX, const double & d_mhtY, const double & mht,
1004  const double & trueMetGuess, const double & mht_offset) {
1005  // all param except trueMetguess unchanged in one event. So can cache agaisnt this one.
1006  //disable cache if (trueMetGuess==trueMetGuesscache) return mhtprobcache;
1007  double mhtprob;
1008  // if(MHtSigma1>0.0 && MHtSigma2>0.0 && MHtGaussFr>0.0)
1009 
1010  // if RANDOMNONUNIF MET already follow the double gaussian parameterisation. So weight should not include it to avoid double counting
1011  // formula to be checked IMHO the two gaussian should be correlated
1012  mhtprob=exp(-0.5*pow(d_mhtX/preparedInput.m_MHtSigma1,2))+preparedInput.m_MHtGaussFr*exp(-0.5*pow(d_mhtX/preparedInput.m_MHtSigma2,2));
1013  mhtprob*=(exp(-0.5*pow(d_mhtY/preparedInput.m_MHtSigma1,2))+preparedInput.m_MHtGaussFr*exp(-0.5*pow(d_mhtY/preparedInput.m_MHtSigma2,2)));
1014 
1015  const double n_arg=(mht-trueMetGuess-mht_offset)/preparedInput.m_MHtSigma1;
1016  mhtprob*=exp(-0.25*pow(n_arg,2)); // assuming sqrt(2)*sigma
1017  return mhtprob;
1018 }
1019 
1020 double MissingMassProb::MHtProbabilityHH(MissingMassInput& preparedInput, const double & d_mhtX, const double & d_mhtY, const double & mht,
1021  const double & trueMetGuess, const double & mht_offset) {
1022  double prob=1.0;
1023 
1024  // updated for final cuts, May 21 2012
1025  // should be merged
1026  prob=prob*(0.0256*TMath::Gaus(d_mhtX,0.0,preparedInput.m_MHtSigma1)+0.01754*TMath::Gaus(d_mhtX,0.0,preparedInput.m_MHtSigma2));
1027  prob=prob*(0.0256*TMath::Gaus(d_mhtY,0.0,preparedInput.m_MHtSigma1)+0.01754*TMath::Gaus(d_mhtY,0.0,preparedInput.m_MHtSigma2));
1028  const double n_arg=(mht-trueMetGuess-mht_offset)/5.7; // this sigma is different from MHtSigma1; actually it depends on dPhi
1029  prob=prob*exp(-0.5*pow(n_arg,2))/(5.7*sqrt(2.0*TMath::Pi())); // assuming sigma from above line
1030 
1031  return prob;
1032 }
1033 
1034 //SpeedUp static instantation
1035 // first index is the calibration set : 0: MMC2011, 1:MMC2012
1036 // second index is the decay 0 : lepton, 1 : 1 prong, 2 3 prong
1037 thread_local double MissingMassProb::s_fit_param[2][3][6][5];
1038 // first parameter: 0- for 1-prong; 1- for 3-prong
1039 thread_local double MissingMassProb::s_ter_sigma_par[2][10][3];
1040 
1041 // returns dTheta3D probability based on ATLAS parameterization
1042 double MissingMassProb::dTheta3d_probabilityFast(MissingMassInput& preparedInput, const int & tau_type,const double & dTheta3d,const double & P_tau) {
1043  double prob=1.0E-10;
1044  int tau_code; // 0: l, 1:1-prong, 2:3-prong
1045  if(tau_type==8) tau_code = 0;
1046  else if(tau_type>=0 && tau_type<=2) tau_code = 1;
1047  else if(tau_type>=3 && tau_type<=5) tau_code = 2;
1048  else if(tau_type==6) return prob;
1049  else
1050  {
1051  Warning("DiTauMassTools", "---- WARNING in MissingMassCalculator::dTheta3d_probabilityFast() ----");
1052  Warning("DiTauMassTools", "%s", ("..... wrong tau_type="+std::to_string(tau_type)).c_str());
1053  Warning("DiTauMassTools", "%s", ("..... returning prob="+std::to_string(prob)).c_str());
1054  Warning("DiTauMassTools", "____________________________________________________________");
1055  return prob;
1056  }
1057 
1058 
1059  double myDelThetaParam[6];
1060 
1061  for (int i=0;i<6;++i)
1062  {
1067  myDelThetaParam[i]=dTheta3Dparam(i,tau_code,P_tau,s_fit_param[1][tau_code][i]);
1068  }
1069  double dTheta3dVal=dTheta3d;
1070 
1071  if (tau_type==8) prob=myDelThetaLepFunc(&dTheta3dVal, myDelThetaParam);
1072  else prob=myDelThetaHadFunc(&dTheta3dVal, myDelThetaParam);
1073 
1074  if (false)
1075  {
1076 
1077  if( preparedInput.m_fUseVerbose==1 && (prob>1.0 || prob<0.0))
1078  {
1079  Warning("DiTauMassTools", "---- WARNING in MissingMassCalculator::dTheta3d_probabilityFast() ----");
1080  Warning("DiTauMassTools", "%s", ("..... wrong probability="+std::to_string(prob)).c_str());
1081  Warning("DiTauMassTools", "%s", ("..... debugging: tau_type="+std::to_string(tau_type)+"dTheta3d="+std::to_string(dTheta3d)+" P_tau="+std::to_string(P_tau)).c_str());
1082  Warning("DiTauMassTools", "____________________________________________________________");
1083  prob=1.0E-10;
1084  }
1085  }
1086 
1087  return prob;
1088 }
1089 
1090 // dTheta probability density function for hadronic taus
1091 double MissingMassProb::myDelThetaHadFunc(double *x, double *par)
1092 {
1093  double fitval=1.0E-10;
1094  if(x[0]>TMath::Pi() || x[0]<0.0) return fitval;
1095  const double arg=x[0];
1096  const double arg_L=arg;
1097  const double mean=par[1];
1098  const double sigmaG=par[2];
1099  const double mpv=par[3];
1100  const double sigmaL=par[4];
1101 
1106  const double norm=sqrt(2.0*TMath::Pi());
1107  const double g1=TMath::Gaus(arg,mean,sigmaG)/norm;
1108  const double g2=TMath::Landau(arg_L,mpv,sigmaL)/norm;
1109  fitval=par[0]*g1/sigmaG+par[5]*g2/sigmaL;
1110  }
1111 
1112  if(fitval<0.0) return 0.0;
1113  return fitval;
1114 }
1115 
1116 // dTheta probability density function for leptonic taus
1117 double MissingMassProb::myDelThetaLepFunc(double *x, double *par)
1118 {
1119  double fitval=1.0E-10;
1120  if(x[0]>TMath::Pi() || x[0]<0.0) return fitval;
1121  double arg=x[0]/par[1];
1122 
1123  double normL=par[5];
1124  if(normL<0.0) normL=0.0;
1125 
1126  if(arg<1) arg=sqrt(std::abs(arg));
1127  else arg=arg*arg;
1128  const double arg_L=x[0];
1129  const double mean=1.0;
1130  const double sigmaG=par[2];
1131  const double mpv=par[3];
1132  const double sigmaL=par[4];
1133  const double g1=normL*TMath::Gaus(arg,mean,sigmaG);
1134  const double g2=TMath::Landau(arg_L,mpv,sigmaL);
1135  fitval=par[0]*(g1+g2)/(1.0+normL);
1136  if(fitval<0.0) return 0.0;
1137  return fitval;
1138 }
1139 
1140 // returns parameters for dTheta3D pdf
1141 double MissingMassProb::dTheta3Dparam(const int & parInd, const int & tau_type, const double & P_tau, const double *par) {
1142  //tau_type 0: l, 1:1-prong, 3:3-prong
1143  if(P_tau<0.0) return 0.0;
1144 
1145 
1146  if(parInd==0) {
1151  return (par[0]+par[1]*P_tau+par[2]*pow(P_tau,2)+par[3]*pow(P_tau,3)+par[4]*pow(P_tau,4))*0.00125;
1152  }
1153  }
1154  else { // parInd==0
1159  if(tau_type==0) return par[0]*(exp(-par[1]*P_tau)+par[2]/P_tau)+par[3]+par[4]*P_tau;
1160  else return par[0]*(exp(-par[1]*sqrt(P_tau))+par[2]/P_tau)+par[3]+par[4]*P_tau;
1161  }
1162  }
1163  return 0.;
1164 }
1165 
1167  // compute MET resolution (eventually use HT)
1168  if(preparedInput.m_METsigmaP<0.1 || preparedInput.m_METsigmaL<0.1)
1169  {
1171  {
1172  if(preparedInput.m_fUseVerbose==1) { Info("DiTauMassTools", "Attempting to set LFV MMC settings"); }
1173  double mT1 = mT(preparedInput.m_vistau1,preparedInput.m_MetVec);
1174  double mT2 = mT(preparedInput.m_vistau2,preparedInput.m_MetVec);
1175  int sr_switch = 0;
1176  if(preparedInput.m_vistau1.M()<0.12 && preparedInput.m_vistau2.M()<0.12) // lep-lep case
1177  {
1178  if(preparedInput.m_LFVmode==0) // e+tau(->mu) case
1179  {
1180  if(preparedInput.m_vistau1.M()<0.05 && preparedInput.m_vistau2.M()>0.05)
1181  {
1182  if(mT1>mT2) sr_switch = 0; // SR1
1183  else sr_switch = 1; // SR2
1184  }
1185  else
1186  {
1187  if(mT1>mT2) sr_switch = 1; // SR2
1188  else sr_switch = 0; // SR1
1189  }
1190  }
1191  if(preparedInput.m_LFVmode==1) // mu+tau(->e) case
1192  {
1193  if(preparedInput.m_vistau1.M()>0.05 && preparedInput.m_vistau2.M()<0.05)
1194  {
1195  if(mT1>mT2) sr_switch = 0; // SR1
1196  else sr_switch = 1; // SR2
1197  }
1198  else
1199  {
1200  if(mT1>mT2) sr_switch = 1; // SR2
1201  else sr_switch = 0; // SR1
1202  }
1203  }
1204  }
1205  if((preparedInput.m_vistau1.M()<0.12 && preparedInput.m_vistau2.M()>0.12) || (preparedInput.m_vistau2.M()<0.12 && preparedInput.m_vistau1.M()>0.12)) // lep-had case
1206  {
1207  if(preparedInput.m_vistau1.M()<0.12 && preparedInput.m_vistau2.M()>0.12)
1208  {
1209  if(mT1>mT2) sr_switch = 0; // SR1
1210  else sr_switch = 1; // SR2
1211  }
1212  else
1213  {
1214  if(mT1>mT2) sr_switch = 1; // SR2
1215  else sr_switch = 0; // SR1
1216  }
1217  }
1218 
1219  m_UseHT = false;
1220  if(preparedInput.m_Njet25==0) // 0-jet
1221  {
1222  double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1223  double METresScale;
1224  double METoffset;
1225  if(sr_switch==0)
1226  {
1227  METresScale = 0.41*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1228  METoffset = 7.36*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1229  }
1230  else
1231  {
1232  METresScale = 0.34*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1233  METoffset = 6.61*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1234  }
1235  if(preparedInput.m_fUseVerbose==1) {
1236  Info("DiTauMassTools", "%s", ("SumEt = "+std::to_string(preparedInput.m_SumEt)).c_str());
1237  Info("DiTauMassTools", "%s", ("METoffset = "+std::to_string(METoffset)).c_str());
1238  Info("DiTauMassTools", "%s", ("METresScale = "+std::to_string(METresScale)).c_str());
1239  }
1240 
1241  double sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1242  preparedInput.m_METsigmaP = sigma;
1243  preparedInput.m_METsigmaL = sigma;
1244  if(preparedInput.m_fUseVerbose==1) {
1245  Info("DiTauMassTools", "%s", ("=> METsigmaP = "+std::to_string(preparedInput.m_METsigmaP)).c_str());
1246  Info("DiTauMassTools", "%s", ("=> METsigmaL = "+std::to_string(preparedInput.m_METsigmaL)).c_str());
1247  }
1248  }
1249  if(preparedInput.m_Njet25>0) // Inclusive 1-jet
1250  {
1251  double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1252  double sigma = 0.;
1253  double METresScale;
1254  double METoffset;
1255  if(sr_switch==0)
1256  {
1257  METresScale = 0.38*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1258  METoffset = 7.96*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1259  }
1260  else
1261  {
1262  METresScale = 0.39*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1263  METoffset = 6.61*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1264  }
1265 
1266  // MET resolution can't be perfect in presence of other objects (i.e., electrons, jets, taus), so assume minSumEt = 5.0 for now
1267  sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1268  preparedInput.m_METsigmaP = sigma;
1269  preparedInput.m_METsigmaL = sigma;
1270  } // Njet25>0
1271  }
1272  else //LFV
1273  {
1274  //DRDRMERGE end addition
1275 
1276  if(preparedInput.m_METScanScheme==1) // default for Winter 2012 and further
1277  {
1278  //LEP-HAD
1279  if ( preparedInput.m_tauTypes==TauTypes::lh ) // lephad case
1280  {
1281  //0-jet
1282  if(preparedInput.m_Njet25==0)//0-jet
1283  {
1284  // placeholder for 2019 tune
1288  if(preparedInput.m_MetVec.R()<20.0) // 0-jet low MET case
1289  {
1290  if(std::abs(preparedInput.m_DelPhiTT)>2.95 && m_allowUseHT) // use mHt only if dPhi(lep-tau)>2.95
1291  {
1292  m_UseHT = true;
1293  // giving priority to external settings
1294  if(preparedInput.m_MHtSigma1<0.0) preparedInput.m_MHtSigma1 = 4.822;
1295  if(preparedInput.m_MHtSigma2<0.0) preparedInput.m_MHtSigma2 = 10.31;
1296  if(preparedInput.m_MHtGaussFr<0.0) preparedInput.m_MHtGaussFr = 6.34E-5;
1297  }
1298  else
1299  {
1300  m_UseHT = false;
1301  double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1302  double METresScale = 0.32*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1303  double METoffset = 5.38*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1304  double sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1305  preparedInput.m_METsigmaP = sigma;
1306  preparedInput.m_METsigmaL = sigma;
1307  }
1308  }
1309  else // 0-jet high MET case
1310  {
1311  if(std::abs(preparedInput.m_DelPhiTT)>2.8 && m_allowUseHT) // use mHt only if dPhi(lep-tau)>2.8
1312  {
1313  m_UseHT = true;
1314  // giving priority to external settings
1315  if(preparedInput.m_MHtSigma1<0.0) preparedInput.m_MHtSigma1 = 7.5;
1316  if(preparedInput.m_MHtSigma2<0.0) preparedInput.m_MHtSigma2 = 13.51;
1317  if(preparedInput.m_MHtGaussFr<0.0) preparedInput.m_MHtGaussFr = 6.81E-4;
1318  preparedInput.m_METsigmaP = preparedInput.m_MHtSigma2; // sigma of 2nd Gaussian for missing Ht resolution
1319  preparedInput.m_METsigmaL = preparedInput.m_MHtSigma2;
1320  }
1321  else
1322  {
1323  m_UseHT = false;
1324  double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1325  double METresScale = 0.87*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1326  double METoffset = 4.16*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1327  double sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1328  preparedInput.m_METsigmaP = sigma;
1329  preparedInput.m_METsigmaL = sigma;
1330  }
1331  } // high MET
1332  } // MMC2016MC15C or MMC2019
1333  // 2015 high-mass tune; avergare MET resolution for Mh=600,1000 mass points
1335  {
1336  m_UseHT = false;
1337  double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1338  double METresScale = 0.65*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1339  double METoffset = 5.0*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1340  double sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1341  preparedInput.m_METsigmaP = sigma;
1342  preparedInput.m_METsigmaL = sigma;
1343  } // MMC2015HIGHMASS
1344  } // 0 jet
1345  //1-jet
1346  else if(preparedInput.m_Njet25>0) // Inclusive 1-jet and VBF lep-had categories for Winter 2012
1347  {
1348  double sigmaSyst=0.10; // 10% systematics for now (be conservative)
1349  double sigma=0.;
1350  // 2015 high-mass tune; average MET resolution for Mh=400,600 mass points (they look consistent);
1352  {
1353  double METresScale=0.86*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1354  double METoffset=3.0*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1355  // MET resolution can't be perfect in presence of other objects (i.e., electrons, jets, taus), so assume minSumEt=5.0 for now
1356  sigma= preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1357  }
1358  //2016 mc15c or 2019
1362  {
1363  double x = preparedInput.m_DelPhiTT;
1364  double dphi_scale = x > 0.3 ? 0.9429 - 0.059*x + 0.054*x*x : 0.728;
1365  double METoffset = 1.875*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1366  double METresScale1 = 8.914*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1367  double METresScale2 = -8.53*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1368 
1369 
1370  sigma = preparedInput.m_SumEt > 80.0 ? METoffset + METresScale1*TMath::Log(sqrt(preparedInput.m_SumEt)+METresScale2) : 5.0;
1371  sigma = sigma * dphi_scale;
1372  }
1373  //
1374 
1375  preparedInput.m_METsigmaP=sigma;
1376  preparedInput.m_METsigmaL=sigma;
1377  } // Njet25>0
1378 
1379  } // lep-had
1380 
1381  //HAD-HAD
1382  else if(preparedInput.m_tauTypes==TauTypes::hh) // had-had events
1383  {
1384  if(preparedInput.m_Njet25==0 && m_mmcCalibrationSet==MMCCalibrationSet::MMC2015HIGHMASS) //0-jet high mass hadhad
1385  {
1386  // 2015 high-mass tune; average of all mass points
1387  // double METresScale=-1.;
1388  // double METoffset=-1.;
1389  double sigmaSyst=0.10; // 10% systematics for now (be conservative)
1390 
1391  double METresScale=0.9*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1392  double METoffset=-1.8*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1393  double sigma= preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : std::abs(METoffset);
1394  preparedInput.m_METsigmaP=sigma;
1395  preparedInput.m_METsigmaL=sigma;
1396 
1397  }
1398  else if(preparedInput.m_Njet25==0 &&
1402  {
1403  double sigmaSyst=0.10; // 10% systematics for now (be conservative)
1404  double x = preparedInput.m_DelPhiTT;
1405  double dphi_scale = x > 2.5 ? 11.0796 - 4.61236*x + 0.423617*x*x : 2.;
1406  double METoffset = -8.51013*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1407  double METresScale1 = 8.54378*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1408  double METresScale2 = -3.97146*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1409  double sigma= preparedInput.m_SumEt>80.0 ? METoffset+METresScale1*TMath::Log(sqrt(preparedInput.m_SumEt)+METresScale2) : 5.;
1410  sigma = sigma*dphi_scale;
1411 
1412  preparedInput.m_METsigmaP=sigma;
1413  preparedInput.m_METsigmaL=sigma;
1414 
1415  }
1416  else if(preparedInput.m_Njet25==0 && m_allowUseHT) // 0-jet high MET had-had category for Winter 2012
1417  {
1418 
1419  m_UseHT=true; // uncomment this line to enable HT also for HH (crucial)
1420  // updated for final cuts, may 21 2012
1421  if(preparedInput.m_MHtSigma1<0.0) preparedInput.m_MHtSigma1=5.972;
1422  if(preparedInput.m_MHtSigma2<0.0) preparedInput.m_MHtSigma2=13.85;
1423  // if(MHtGaussFr<0.0) MHtGaussFr=0.4767; // don't directly use 2nd fraction
1424  }
1425  //1-jet
1426  else // Inclusive 1-jet and VBF categories
1427  {
1428  double METresScale=-1.;
1429  double METoffset=-1.;
1430  double sigmaSyst=0.10; // 10% systematics for now (be conservative)
1431 
1432  // previous value in trunk
1434  // 2015 high-mass tune; average of all mass points
1435  METresScale = 1.1*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1436  METoffset = -5.0*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1437  }
1438  // MET resolution can't be perfect in presence of other objects (i.e., electrons, jets, taus), so assume minSumEt=5.0 for now
1439  double sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : std::abs(METoffset);
1440 
1444  double x = preparedInput.m_DelPhiTT;
1445  double dphi_scale = x > 0.6 ? 1.42047 - 0.666644*x + 0.199986*x*x : 1.02;
1446  METoffset = 1.19769*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1447  double METresScale1 = 5.61687*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1448  double METresScale2 = -4.2076*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1449  sigma= preparedInput.m_SumEt>115.0 ? METoffset+METresScale1*TMath::Log(sqrt(preparedInput.m_SumEt)+METresScale2) : 12.1;
1450  sigma = sigma*dphi_scale;
1451  } //for hh 2016 mc15c
1452 
1453  preparedInput.m_METsigmaP = sigma;
1454  preparedInput.m_METsigmaL = sigma;
1455 
1456  }// 1 jet
1457  } // had-had
1458  //LEP-LEP
1459  else if(preparedInput.m_tauTypes==TauTypes::ll) // setup for LEP-LEP channel
1460  {
1461  if(m_mmcCalibrationSet==MMCCalibrationSet::MMC2015HIGHMASS) // placeholder for 2015 high-mass tune; for now it is the same as 2012
1462  {
1463  m_UseHT = false;
1464  double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1465  double METresScale = -1.0;
1466  double METoffset = -1.0;
1467  double sigma = 5.0;
1468  // tune is based on cuts for Run-1 paper analysis
1469  if(preparedInput.m_Njet25==0)
1470  {
1471  // use tune for emebedding
1472  METresScale=-0.4307*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1473  METoffset=7.06*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1474  double METresScale2=0.07693*(1.0+preparedInput.m_METresSyst*sigmaSyst); // quadratic term
1475  // this is a tune for Higgs125
1476  // METresScale=-0.5355*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1477  // METoffset=11.5*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1478  // double METresScale2=0.07196*(1.0+preparedInput.m_METresSyst*sigmaSyst); // quadratic term
1479  sigma= preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt)+METresScale2*preparedInput.m_SumEt : METoffset;
1480  }
1481  if(preparedInput.m_Njet25>0)
1482  {
1483  // use tune for embedding
1484  METresScale=0.8149*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1485  METoffset=5.343*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1486  // this is a tune for Higgs125
1487  // METresScale=0.599*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1488  // METoffset=8.223*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1489  sigma= preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1490  }
1491  preparedInput.m_METsigmaP = sigma;
1492  preparedInput.m_METsigmaL = sigma;
1493  } // end of MMC2015HIGHMASS
1494 
1498  // 2016 MC15c + 2019 leplep
1499  {
1500  m_UseHT=false;
1501  double sigmaSyst=0.10; // 10% systematics for now (be conservative)
1502  double METresScale=-1.0;
1503  double METoffset=-1.0;
1504  double sigma=5.0;
1505  double min_sigma = 2.0;
1506  // tune is based on cuts for Run-1 paper analysis
1507  if(preparedInput.m_Njet25==0)
1508  {
1509  // Madgraph Ztautau MET param
1510  METoffset = 4.22581*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1511  METresScale = 0.03818*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1512  double METresScale2= 0.12623;
1513  sigma= preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt)+METresScale2*preparedInput.m_SumEt : min_sigma;
1514  if (m_fUseDphiLL) {
1515  double p0 = 2.60131;
1516  double p1const = 1.22427;
1517  double p2quad = -1.71261;
1518  double DphiLL = std::abs(Phi_mpi_pi(preparedInput.m_vistau1.Phi()-preparedInput.m_vistau2.Phi()));
1519  sigma *= (DphiLL < p0) ? p1const : p1const+
1520  p2quad*p0*p0 - 2*p2quad*p0*DphiLL+p2quad*DphiLL*DphiLL;
1521  }
1522  if (sigma < min_sigma) sigma = min_sigma;
1523  }
1524  if(preparedInput.m_Njet25>0)
1525  {
1526  // Madgraph Ztautau MET param
1527  METoffset = 5.42506*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1528  METresScale = 5.36760*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1529  double METoffset2 = -4.86808*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1530  if (preparedInput.m_SumEt > 0.0) {
1531  double x = sqrt(preparedInput.m_SumEt);
1532  sigma = (x+METoffset2 > 1) ? METoffset+METresScale*log(x+METoffset2) : METoffset;
1533  } else {
1534  sigma = METoffset;
1535  }
1536  if (m_fUseDphiLL) {
1537  double p0 = 2.24786;
1538  double p1const = 0.908597;
1539  double p2quad = 0.544577;
1540  double DphiLL = std::abs(Phi_mpi_pi(preparedInput.m_vistau1.Phi()-preparedInput.m_vistau2.Phi()));
1541  sigma *= (DphiLL < p0) ? p1const : p1const+
1542  p2quad*p0*p0 - 2*p2quad*p0*DphiLL+p2quad*DphiLL*DphiLL;
1543  }
1544  if (sigma < min_sigma) sigma = min_sigma;
1545  }
1546  preparedInput.m_METsigmaP=sigma;
1547  preparedInput.m_METsigmaL=sigma;
1548  }//2016 mc15c
1549 
1550  } // lep-lep
1551 
1552  } //preparedInput.METScanScheme
1553 
1554  if(preparedInput.m_METScanScheme==0) // old scheme with JER
1555  {
1556  if(preparedInput.m_dataType==0 || preparedInput.m_dataType==1) preparedInput.SetMetScanParamsUE(preparedInput.m_SumEt,preparedInput.m_METcovphi,preparedInput.m_dataType);
1557  else preparedInput.SetMetScanParamsUE(preparedInput.m_SumEt,preparedInput.m_METcovphi,0);
1558  }
1559  }
1560  } // end else LFV
1561 
1562  return;
1563 }
DiTauMassTools::MMCCalibrationSet::MMC2024
@ MMC2024
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::MissingMassProb::m_paramVectorAngle
std::vector< TF1 * > m_paramVectorAngle
Definition: MissingMassProb.h:110
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
MissingMassProb.h
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:215
ParticleConstants::PDG2011::tauMassInMeV
constexpr double tauMassInMeV
the mass of the tau (in MeV)
Definition: ParticleConstants.h:32
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
DiTauMassTools::MissingMassProb::setParamAngle
void setParamAngle(const PtEtaPhiMVector &tauvec, int tau, int tautype)
Definition: MissingMassProb.cxx:161
DiTauMassTools::TauTypes::lh
@ lh
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:53
DiTauMassTools::MissingMassProb::m_allowUseHT
bool m_allowUseHT
Definition: MissingMassProb.h:120
DiTauMassTools::MissingMassProb::MetProbabilityWrapper
static double MetProbabilityWrapper(MissingMassProb *prob, MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2)
Definition: MissingMassProb.cxx:29
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
mean
void mean(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition: dependence.cxx:254
scale2
#define scale2
Definition: JetAttributeHisto.cxx:42
DiTauMassTools::MMCCalibrationSet::MMC2015HIGHMASS
@ MMC2015HIGHMASS
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::TauTypes::hh
@ hh
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:53
DiTauMassTools::MissingMassProb::m_fUseTauProbability
bool m_fUseTauProbability
Definition: MissingMassProb.h:122
DiTauMassTools::MissingMassProb::m_mmcCalibrationSet
MMCCalibrationSet::e m_mmcCalibrationSet
Definition: MissingMassProb.h:118
DiTauMassTools::MissingMassInput::m_METsigmaP
double m_METsigmaP
Definition: MissingMassInput.h:62
DiTauMassTools::MissingMassProb::m_formulaAngle1
TF1 * m_formulaAngle1
Definition: MissingMassProb.h:101
DiTauMassTools::MissingMassInput::m_METScanScheme
int m_METScanScheme
Definition: MissingMassInput.h:74
DiTauMassTools::MissingMassProb::m_formulaAngle2
TF1 * m_formulaAngle2
Definition: MissingMassProb.h:102
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
DiTauMassTools::MissingMassInput::SetMetScanParamsUE
void SetMetScanParamsUE(double sumEt, double phi_scan=0.0, int data_code=0)
Definition: MissingMassInput.cxx:89
DMTest::P
P_v1 P
Definition: P.h:23
DiTauMassTools::MissingMassInput::m_MHtGaussFr
double m_MHtGaussFr
Definition: MissingMassInput.h:70
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
DiTauMassTools::MissingMassInput::m_metVec
XYVector m_metVec
Definition: MissingMassInput.h:79
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:217
DiTauMassTools::MissingMassProb::MET
void MET(MissingMassInput &preparedInput)
Definition: MissingMassProb.cxx:1166
DiTauMassTools::MissingMassInput::m_inputMEtT
double m_inputMEtT
Definition: MissingMassInput.h:81
DiTauMassTools::MissingMassInput::m_MHtSigma1
double m_MHtSigma1
Definition: MissingMassInput.h:68
DiTauMassTools::MissingMassInput::m_METsigmaL
double m_METsigmaL
Definition: MissingMassInput.h:63
DiTauMassTools::MissingMassProb::m_fUseDphiLL
bool m_fUseDphiLL
Definition: MissingMassProb.h:124
DiTauMassTools::MissingMassInput::m_vistau1
PtEtaPhiMVector m_vistau1
Definition: MissingMassInput.h:54
DiTauMassTools::MissingMassProb::MissingMassProb
MissingMassProb(MMCCalibrationSet::e aset, const std::string &paramFilePath)
Definition: MissingMassProb.cxx:221
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
DiTauMassTools::MissingMassProb::m_UseHT
bool m_UseHT
Definition: MissingMassProb.h:121
DiTauMassTools::MissingMassProb::m_fParams
TFile * m_fParams
Definition: MissingMassProb.h:117
DiTauMassTools::MissingMassProb::m_formulaNuMass
TF1 * m_formulaNuMass
Definition: MissingMassProb.h:109
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
covarianceTool.prob
prob
Definition: covarianceTool.py:678
Phi_mpi_pi
__HOSTDEV__ double Phi_mpi_pi(double)
Definition: GeoRegion.cxx:7
DiTauMassTools::MissingMassProb::s_fit_param
static thread_local double s_fit_param[2][3][6][5]
Definition: MissingMassProb.h:97
x
#define x
DiTauMassTools::MMCCalibrationSet::LFVMMC2012
@ LFVMMC2012
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
DiTauMassTools::MissingMassProb::TauProbabilityWrapper
static double TauProbabilityWrapper(MissingMassProb *prob, MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2)
Definition: MissingMassProb.cxx:70
DiTauMassTools::MissingMassProb::dTheta3Dparam
double dTheta3Dparam(const int &parInd, const int &tau_type, const double &P_tau, const double *par)
Definition: MissingMassProb.cxx:1141
DiTauMassTools::MissingMassInput::m_METcovphi
double m_METcovphi
Definition: MissingMassInput.h:61
DiTauMassTools::MissingMassProb::s_ter_sigma_par
static thread_local double s_ter_sigma_par[2][10][3]
Definition: MissingMassProb.h:98
DiTauMassTools::MissingMassInput::m_tauTypes
TauTypes::e m_tauTypes
Definition: MissingMassInput.h:85
DiTauMassTools::MissingMassProb::apply
double apply(MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2, bool constant=false, bool oneTau=false, bool twoTau=false)
Definition: MissingMassProb.cxx:540
DiTauMassTools::MissingMassInput::m_Njet25
int m_Njet25
Definition: MissingMassInput.h:66
DiTauMassTools::ignore
void ignore(T &&)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:58
DiTauMassTools::MissingMassProb::mEtAndTauProbabilityWrapper
static double mEtAndTauProbabilityWrapper(MissingMassProb *prob, MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2)
Definition: MissingMassProb.cxx:39
DiTauMassTools
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:24
DiTauMassTools::MissingMassProb::myDelThetaHadFunc
double myDelThetaHadFunc(double *x, double *par)
Definition: MissingMassProb.cxx:1091
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
python.changerun.m1
m1
Definition: changerun.py:30
DiTauMassTools::MissingMassInput::m_htOffset
double m_htOffset
Definition: MissingMassInput.h:83
DiTauMassTools::MissingMassInput::m_SumEt
double m_SumEt
Definition: MissingMassInput.h:64
fitman.g1
g1
Definition: fitman.py:619
DiTauMassTools::MissingMassProb::MnuProbability
double MnuProbability(MissingMassInput &preparedInput, double mnu, double binsize)
Definition: MissingMassProb.cxx:945
ParticleConstants.h
lumiFormat.i
int i
Definition: lumiFormat.py:85
DiTauMassTools::readInParams
void readInParams(TDirectory *dir, MMCCalibrationSet::e aset, std::vector< TF1 * > &lep_numass, std::vector< TF1 * > &lep_angle, std::vector< TF1 * > &lep_ratio, std::vector< TF1 * > &had_angle, std::vector< TF1 * > &had_ratio)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:173
DiTauMassTools::MissingMassInput::m_MHtSigma2
double m_MHtSigma2
Definition: MissingMassInput.h:69
DiTauMassTools::MissingMassProb::dTheta3d_probabilityFastWrapper
static double dTheta3d_probabilityFastWrapper(MissingMassProb *prob, MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2)
Definition: MissingMassProb.cxx:49
DiTauMassTools::MissingMassProb::m_formulaRatioHad1
TF1 * m_formulaRatioHad1
Definition: MissingMassProb.h:107
DiTauMassTools::MissingMassInput::m_MEtT
double m_MEtT
Definition: MissingMassInput.h:82
DiTauMassTools::MissingMassProb::dTheta3d_probabilityNewWrapper
static double dTheta3d_probabilityNewWrapper(MissingMassProb *prob, MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2)
Definition: MissingMassProb.cxx:94
DiTauMassTools::MMCCalibrationSet::MMC2019
@ MMC2019
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::mT
double mT(const VectorType &vec, const XYVector &met_vec)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:72
DiTauMassTools::MissingMassProb::dTheta3d_probabilityFast
double dTheta3d_probabilityFast(MissingMassInput &preparedInput, const int &tau_type, const double &dTheta3d, const double &P_tau)
Definition: MissingMassProb.cxx:1042
DiTauMassTools::MissingMassProb::m_formulaRatioLep2
TF1 * m_formulaRatioLep2
Definition: MissingMassProb.h:106
hist_file_dump.f
f
Definition: hist_file_dump.py:140
GEV
#define GEV
Definition: PrintPhotonSF.cxx:25
DiTauMassTools::MMCCalibrationSet::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
dumpNswErrorDb.constant
def constant
Definition: dumpNswErrorDb.py:28
DiTauMassTools::MissingMassInput::m_dataType
int m_dataType
Definition: MissingMassInput.h:60
DiTauMassTools::MissingMassInput::m_fUseVerbose
bool m_fUseVerbose
Definition: MissingMassInput.h:78
DiTauMassTools::MissingMassProb::setParamRatio
void setParamRatio(int tau, int tautype)
Definition: MissingMassProb.cxx:190
create_dcsc_inputs_sqlite.arg
list arg
Definition: create_dcsc_inputs_sqlite.py:48
DiTauMassTools::MissingMassProb::m_paramVectorRatio
std::vector< TF1 * > m_paramVectorRatio
Definition: MissingMassProb.h:112
DiTauMassTools::MissingMassInput::m_tlv_tmp
PtEtaPhiMVector m_tlv_tmp
Definition: MissingMassInput.h:80
DiTauMassTools::MissingMassProb::m_paramVectorRatioLep
std::vector< TF1 * > m_paramVectorRatioLep
Definition: MissingMassProb.h:113
DiTauMassTools::MissingMassInput::m_MetVec
XYVector m_MetVec
Definition: MissingMassInput.h:53
fitman.g2
g2
Definition: fitman.py:624
PathResolver.h
DiTauMassTools::MissingMassProb::m_probListConstant
std::list< std::function< double(MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2)> > m_probListConstant
Definition: MissingMassProb.h:92
DiTauMassTools::MissingMassProb::m_formulaRatioLep1
TF1 * m_formulaRatioLep1
Definition: MissingMassProb.h:105
DiTauMassTools::MissingMassProb::MetProbability
double MetProbability(MissingMassInput &preparedInput, const double &met1, const double &met2, const double &MetSigma1, const double &MetSigma2)
Definition: MissingMassProb.cxx:559
DiTauMassTools::MissingMassInput::m_DelPhiTT
double m_DelPhiTT
Definition: MissingMassInput.h:67
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:28
VP1PartSpect::E
@ E
Definition: VP1PartSpectFlags.h:21
DiTauMassTools::MissingMassProb::MHtProbabilityHH
double MHtProbabilityHH(MissingMassInput &preparedInput, const double &d_mhtX, const double &d_mhtY, const double &mht, const double &trueMetGuess, const double &mht_offset)
Definition: MissingMassProb.cxx:1020
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
DiTauMassTools::MissingMassProb::m_formulaRatio1
TF1 * m_formulaRatio1
Definition: MissingMassProb.h:103
PathResolverFindCalibFile
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
Definition: PathResolver.cxx:283
DiTauMassTools::MissingMassProb::TauProbabilityNewWrapper
static double TauProbabilityNewWrapper(MissingMassProb *prob, MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2)
Definition: MissingMassProb.cxx:122
DiTauMassTools::MissingMassInput::m_vistau2
PtEtaPhiMVector m_vistau2
Definition: MissingMassInput.h:55
DiTauMassTools::MissingMassProb::m_formulaRatioHad2
TF1 * m_formulaRatioHad2
Definition: MissingMassProb.h:108
DiTauMassTools::MissingMassProb
Definition: MissingMassProb.h:28
DiTauMassTools::MissingMassProb::m_fUseMnuProbability
bool m_fUseMnuProbability
Definition: MissingMassProb.h:123
DiTauMassTools::MissingMassProb::setParamNuMass
void setParamNuMass()
Definition: MissingMassProb.cxx:153
DiTauMassTools::MissingMassProb::m_probListTwoTau
std::list< std::function< double(MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2)> > m_probListTwoTau
Definition: MissingMassProb.h:94
DiTauMassTools::MissingMassProb::myDelThetaLepFunc
double myDelThetaLepFunc(double *x, double *par)
Definition: MissingMassProb.cxx:1117
DiTauMassTools::MissingMassProb::TauProbabilityLFV
double TauProbabilityLFV(MissingMassInput &preparedInput, const int &type1, const PtEtaPhiMVector &vis1, const PtEtaPhiMVector &nu1)
Definition: MissingMassProb.cxx:620
DiTauMassTools::MissingMassInput::m_inputMEtX
double m_inputMEtX
Definition: MissingMassInput.h:81
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
DiTauMassTools::MissingMassInput
Definition: MissingMassInput.h:21
DiTauMassTools::MMCCalibrationSet::MMC2016MC15C
@ MMC2016MC15C
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::MissingMassProb::m_paramVectorAngleLep
std::vector< TF1 * > m_paramVectorAngleLep
Definition: MissingMassProb.h:111
DiTauMassTools::MissingMassInput::m_METresSyst
int m_METresSyst
Definition: MissingMassInput.h:75
DiTauMassTools::MissingMassProb::TauProbability
double TauProbability(MissingMassInput &preparedInput, const int &type1, const PtEtaPhiMVector &vis1, const PtEtaPhiMVector &nu1, const int &type2, const PtEtaPhiMVector &vis2, const PtEtaPhiMVector &nu2)
Definition: MissingMassProb.cxx:657
DiTauMassTools::MissingMassProb::m_probListOneTau
std::list< std::function< double(MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2)> > m_probListOneTau
Definition: MissingMassProb.h:93
DiTauMassTools::MissingMassProb::MnuProbabilityNewWrapper
static double MnuProbabilityNewWrapper(MissingMassProb *prob, MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2)
Definition: MissingMassProb.cxx:136
scale1
#define scale1
Definition: JetAttributeHisto.cxx:41
DiTauMassTools::MissingMassProb::MHtProbability
double MHtProbability(MissingMassInput &preparedInput, const double &d_mhtX, const double &d_mhtY, const double &mht, const double &trueMetGuess, const double &mht_offset)
Definition: MissingMassProb.cxx:1003
DiTauMassTools::MMCCalibrationSet::UPGRADE
@ UPGRADE
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
DiTauMassTools::MissingMassProb::MnuProbabilityWrapper
static double MnuProbabilityWrapper(MissingMassProb *prob, MissingMassInput &preparedInput, const int &tau_type1, const int &tau_type2, const PtEtaPhiMVector &tauvec1, const PtEtaPhiMVector &tauvec2, const PtEtaPhiMVector nuvec1, const PtEtaPhiMVector &nuvec2)
Definition: MissingMassProb.cxx:78
DiTauMassTools::Angle
double Angle(const VectorType1 &vec1, const VectorType2 &vec2)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:61
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
python.SystemOfUnits.m2
float m2
Definition: SystemOfUnits.py:107
DiTauMassTools::MissingMassInput::m_inputMEtY
double m_inputMEtY
Definition: MissingMassInput.h:81
DiTauMassTools::MissingMassProb::m_formulaRatio2
TF1 * m_formulaRatio2
Definition: MissingMassProb.h:104
DiTauMassTools::MissingMassInput::m_LFVmode
int m_LFVmode
Definition: MissingMassInput.h:84
DiTauMassTools::MissingMassProb::~MissingMassProb
~MissingMassProb()
Definition: MissingMassProb.cxx:537
TRTCalib_cfilter.p0
p0
Definition: TRTCalib_cfilter.py:129
DiTauMassTools::MissingMassProb::m_paramVectorNuMass
std::vector< TF1 * > m_paramVectorNuMass
Definition: MissingMassProb.h:114
DiTauMassTools::TauTypes::ll
@ ll
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:53
DiTauMassTools::MissingMassProb::mEtAndTauProbability
double mEtAndTauProbability(MissingMassInput &preparedInput)
Definition: MissingMassProb.cxx:580