ATLAS Offline Software
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Public Attributes | Private Attributes | Static Private Attributes | List of all members
DiTauMassTools::MissingMassProb Class Reference

#include <MissingMassProb.h>

Collaboration diagram for DiTauMassTools::MissingMassProb:

Public Member Functions

 MissingMassProb (MMCCalibrationSet::e aset, const std::string &paramFilePath)
 ~MissingMassProb ()
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)
void setParamAngle (const PtEtaPhiMVector &tauvec, int tau, int tautype)
void setParamRatio (int tau, int tautype)
void setParamNuMass ()
TF1 * GetFormulaAngle1 ()
TF1 * GetFormulaAngle2 ()
TF1 * GetFormulaRatio1 ()
TF1 * GetFormulaRatio2 ()
TF1 * GetFormulaNuMass ()
void SetAllowUseHT (bool allowUseHT)
bool GetAllowUseHT ()
void SetUseHT (bool val)
bool GetUseHT ()
void SetUseTauProbability (bool val)
bool GetUseTauProbability ()
void SetUseMnuProbability (bool val)
bool GetUseMnuProbability ()
void SetUseDphiLL (bool val)
bool GetUseDphiLL ()
double MetProbability (MissingMassInput &preparedInput, const double &met1, const double &met2, const double &MetSigma1, const double &MetSigma2)
double dTheta3Dparam (const int &parInd, const int &tau_type, const double &P_tau, const double *par)
double dTheta3d_probabilityFast (MissingMassInput &preparedInput, const int &tau_type, const double &dTheta3d, const double &P_tau)
double myDelThetaHadFunc (double *x, double *par)
double myDelThetaLepFunc (double *x, double *par)
double MHtProbability (MissingMassInput &preparedInput, const double &d_mhtX, const double &d_mhtY, const double &mht, const double &trueMetGuess, const double &mht_offset)
double MHtProbabilityHH (MissingMassInput &preparedInput, const double &d_mhtX, const double &d_mhtY, const double &mht, const double &trueMetGuess, const double &mht_offset)
void MET (MissingMassInput &preparedInput)
double mEtAndTauProbability (MissingMassInput &preparedInput)
double MnuProbability (MissingMassInput &preparedInput, double mnu, double binsize)
double MnuProbability (MissingMassInput &preparedInput, double mnu)
double TauProbability (MissingMassInput &preparedInput, const int &type1, const PtEtaPhiMVector &vis1, const PtEtaPhiMVector &nu1, const int &type2, const PtEtaPhiMVector &vis2, const PtEtaPhiMVector &nu2)
double TauProbability (MissingMassInput &preparedInput, const int &type1, const PtEtaPhiMVector &vis1, const PtEtaPhiMVector &nu1, const int &type2, const PtEtaPhiMVector &vis2, const PtEtaPhiMVector &nu2, const double &detmet)
double TauProbabilityLFV (MissingMassInput &preparedInput, const int &type1, const PtEtaPhiMVector &vis1, const PtEtaPhiMVector &nu1)

Static Public Member Functions

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)
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)
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)
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)
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)
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)
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)
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)

Public Attributes

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
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
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

Private Attributes

TF1 * m_formulaAngle1 = new TF1("formulaAngle1", "[0]*exp(-[2]*(log((x+[3])/[1]))**2)")
TF1 * m_formulaAngle2 = new TF1("formulaAngle2", "[0]*exp(-[2]*(log((x+[3])/[1]))**2)")
TF1 * m_formulaRatio1
TF1 * m_formulaRatio2
TF1 * m_formulaRatioLep1 = new TF1("formulaRatio1", "gaus(0)+expo(3)")
TF1 * m_formulaRatioLep2 = new TF1("formulaRatio2", "gaus(0)+expo(3)")
TF1 * m_formulaRatioHad1 = new TF1("formulaRatio1", "gaus(0)")
TF1 * m_formulaRatioHad2 = new TF1("formulaRatio2", "gaus(0)")
TF1 * m_formulaNuMass = new TF1("formulaNuMass", "pol6")
std::vector< TF1 * > m_paramVectorAngle
std::vector< TF1 * > m_paramVectorAngleLep
std::vector< TF1 * > m_paramVectorRatio
std::vector< TF1 * > m_paramVectorRatioLep
std::vector< TF1 * > m_paramVectorNuMass
std::string m_paramFilePath
TFile * m_fParams
MMCCalibrationSet::e m_mmcCalibrationSet
bool m_allowUseHT
bool m_UseHT
bool m_fUseTauProbability
bool m_fUseMnuProbability
bool m_fUseDphiLL

Static Private Attributes

static thread_local double s_fit_param [2][3][6][5]
static thread_local double s_ter_sigma_par [2][10][3]

Detailed Description

Definition at line 28 of file MissingMassProb.h.

Constructor & Destructor Documentation

◆ MissingMassProb()

MissingMassProb::MissingMassProb ( MMCCalibrationSet::e aset,
const std::string & paramFilePath )

MMC2011 parameterisation

MMC2012 parameterisation

Definition at line 221 of file MissingMassProb.cxx.

221 {
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 }
231 if (aset == MMCCalibrationSet::MMC2019 || aset == MMCCalibrationSet::MMC2024) {
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){
237 readInParams(m_fParams, aset, m_paramVectorNuMass, m_paramVectorAngleLep, m_paramVectorRatioLep, m_paramVectorAngle, m_paramVectorRatio);
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}
PathResolverFindCalibFile
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
Definition PathResolver.cxx:325
DiTauMassTools::MissingMassProb::m_paramVectorRatioLep
std::vector< TF1 * > m_paramVectorRatioLep
Definition MissingMassProb.h:113
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_paramVectorNuMass
std::vector< TF1 * > m_paramVectorNuMass
Definition MissingMassProb.h:114
DiTauMassTools::MissingMassProb::s_ter_sigma_par
static thread_local double s_ter_sigma_par[2][10][3]
Definition MissingMassProb.h:98
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_paramVectorAngle
std::vector< TF1 * > m_paramVectorAngle
Definition MissingMassProb.h:110
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::s_fit_param
static thread_local double s_fit_param[2][3][6][5]
Definition MissingMassProb.h:97
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
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::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::MissingMassProb::m_paramVectorRatio
std::vector< TF1 * > m_paramVectorRatio
Definition MissingMassProb.h:112
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::MissingMassProb::m_paramVectorAngleLep
std::vector< TF1 * > m_paramVectorAngleLep
Definition MissingMassProb.h:111
DiTauMassTools::MissingMassProb::m_mmcCalibrationSet
MMCCalibrationSet::e m_mmcCalibrationSet
Definition MissingMassProb.h:118
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::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::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::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

◆ ~MissingMassProb()

MissingMassProb::~MissingMassProb ( )

Definition at line 537 of file MissingMassProb.cxx.

537 {
538}

Member Function Documentation

◆ apply()

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 = false,
bool oneTau = false,
bool twoTau = false )

Definition at line 540 of file MissingMassProb.cxx.

540 {
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}
hist_file_dump.f
f
Definition hist_file_dump.py:140

◆ dTheta3d_probabilityFast()

double MissingMassProb::dTheta3d_probabilityFast ( MissingMassInput & preparedInput,
const int & tau_type,
const double & dTheta3d,
const double & P_tau )

Definition at line 1042 of file MissingMassProb.cxx.

1042 {
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 {
1063 if(m_mmcCalibrationSet==MMCCalibrationSet::UPGRADE
1064 || m_mmcCalibrationSet==MMCCalibrationSet::LFVMMC2012
1065 || m_mmcCalibrationSet==MMCCalibrationSet::MMC2015HIGHMASS)
1066 myDelThetaParam[i]=dTheta3Dparam(i,tau_code,P_tau,s_fit_param[1][tau_code][i]);
1067 }
1068 double dTheta3dVal=dTheta3d;
1069
1070 if (tau_type==8) prob=myDelThetaLepFunc(&dTheta3dVal, myDelThetaParam);
1071 else prob=myDelThetaHadFunc(&dTheta3dVal, myDelThetaParam);
1072
1073 if (false)
1074 {
1075
1076 if( preparedInput.m_fUseVerbose==1 && (prob>1.0 || prob<0.0))
1077 {
1078 Warning("DiTauMassTools", "---- WARNING in MissingMassCalculator::dTheta3d_probabilityFast() ----");
1079 Warning("DiTauMassTools", "%s", ("..... wrong probability="+std::to_string(prob)).c_str());
1080 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());
1081 Warning("DiTauMassTools", "____________________________________________________________");
1082 prob=1.0E-10;
1083 }
1084 }
1085
1086 return prob;
1087}
DiTauMassTools::MissingMassProb::myDelThetaLepFunc
double myDelThetaLepFunc(double *x, double *par)
Definition MissingMassProb.cxx:1115
DiTauMassTools::MissingMassProb::dTheta3Dparam
double dTheta3Dparam(const int &parInd, const int &tau_type, const double &P_tau, const double *par)
Definition MissingMassProb.cxx:1139
DiTauMassTools::MissingMassProb::myDelThetaHadFunc
double myDelThetaHadFunc(double *x, double *par)
Definition MissingMassProb.cxx:1090

◆ dTheta3d_probabilityFastWrapper()

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 )
static

Definition at line 49 of file MissingMassProb.cxx.

49 {
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}
DiTauMassTools::Angle
double Angle(const VectorType1 &vec1, const VectorType2 &vec2)
Definition PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:61

◆ dTheta3d_probabilityNewWrapper()

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 )
static

Definition at line 94 of file MissingMassProb.cxx.

94 {
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}

◆ dTheta3Dparam()

double MissingMassProb::dTheta3Dparam ( const int & parInd,
const int & tau_type,
const double & P_tau,
const double * par )

Definition at line 1139 of file MissingMassProb.cxx.

1139 {
1140 //tau_type 0: l, 1:1-prong, 3:3-prong
1141 if(P_tau<0.0) return 0.0;
1142
1143
1144 if(parInd==0) {
1145 if (m_mmcCalibrationSet==MMCCalibrationSet::MMC2015HIGHMASS
1146 || m_mmcCalibrationSet==MMCCalibrationSet::UPGRADE
1147 || m_mmcCalibrationSet==MMCCalibrationSet::LFVMMC2012){
1148 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;
1149 }
1150 }
1151 else { // parInd==0
1152 if (m_mmcCalibrationSet==MMCCalibrationSet::MMC2015HIGHMASS
1153 || m_mmcCalibrationSet==MMCCalibrationSet::UPGRADE
1154 || m_mmcCalibrationSet==MMCCalibrationSet::LFVMMC2012){
1155 if(tau_type==0) return par[0]*(exp(-par[1]*P_tau)+par[2]/P_tau)+par[3]+par[4]*P_tau;
1156 else return par[0]*(exp(-par[1]*sqrt(P_tau))+par[2]/P_tau)+par[3]+par[4]*P_tau;
1157 }
1158 }
1159 return 0.;
1160}
pow
constexpr int pow(int base, int exp) noexcept
Definition ap_fixedTest.cxx:15

◆ GetAllowUseHT()

bool DiTauMassTools::MissingMassProb::GetAllowUseHT ( )
inline

Definition at line 46 of file MissingMassProb.h.

46{ return m_allowUseHT;}

◆ GetFormulaAngle1()

TF1 * DiTauMassTools::MissingMassProb::GetFormulaAngle1 ( )
inline

Definition at line 39 of file MissingMassProb.h.

39{return m_formulaAngle1;}

◆ GetFormulaAngle2()

TF1 * DiTauMassTools::MissingMassProb::GetFormulaAngle2 ( )
inline

Definition at line 40 of file MissingMassProb.h.

40{return m_formulaAngle2;}

◆ GetFormulaNuMass()

TF1 * DiTauMassTools::MissingMassProb::GetFormulaNuMass ( )
inline

Definition at line 43 of file MissingMassProb.h.

43{return m_formulaNuMass;}

◆ GetFormulaRatio1()

TF1 * DiTauMassTools::MissingMassProb::GetFormulaRatio1 ( )
inline

Definition at line 41 of file MissingMassProb.h.

41{return m_formulaRatio1;}

◆ GetFormulaRatio2()

TF1 * DiTauMassTools::MissingMassProb::GetFormulaRatio2 ( )
inline

Definition at line 42 of file MissingMassProb.h.

42{return m_formulaRatio2;}

◆ GetUseDphiLL()

bool DiTauMassTools::MissingMassProb::GetUseDphiLL ( )
inline

Definition at line 58 of file MissingMassProb.h.

58{ return m_fUseDphiLL; }

◆ GetUseHT()

bool DiTauMassTools::MissingMassProb::GetUseHT ( )
inline

Definition at line 49 of file MissingMassProb.h.

49{ return m_UseHT; } // if use HT

◆ GetUseMnuProbability()

bool DiTauMassTools::MissingMassProb::GetUseMnuProbability ( )
inline

Definition at line 55 of file MissingMassProb.h.

55{ return m_fUseMnuProbability; }

◆ GetUseTauProbability()

bool DiTauMassTools::MissingMassProb::GetUseTauProbability ( )
inline

Definition at line 52 of file MissingMassProb.h.

52{ return m_fUseTauProbability; }

◆ MET()

void MissingMassProb::MET ( MissingMassInput & preparedInput)

Definition at line 1162 of file MissingMassProb.cxx.

1162 {
1163 // compute MET resolution (eventually use HT)
1164 if(preparedInput.m_METsigmaP<0.1 || preparedInput.m_METsigmaL<0.1)
1165 {
1166 if(m_mmcCalibrationSet==MMCCalibrationSet::LFVMMC2012)
1167 {
1168 if(preparedInput.m_fUseVerbose==1) { Info("DiTauMassTools", "Attempting to set LFV MMC settings"); }
1169 double mT1 = mT(preparedInput.m_vistau1,preparedInput.m_MetVec);
1170 double mT2 = mT(preparedInput.m_vistau2,preparedInput.m_MetVec);
1171 int sr_switch = 0;
1172 if(preparedInput.m_vistau1.M()<0.12 && preparedInput.m_vistau2.M()<0.12) // lep-lep case
1173 {
1174 if(preparedInput.m_LFVmode==0) // e+tau(->mu) case
1175 {
1176 if(preparedInput.m_vistau1.M()<0.05 && preparedInput.m_vistau2.M()>0.05)
1177 {
1178 if(mT1>mT2) sr_switch = 0; // SR1
1179 else sr_switch = 1; // SR2
1180 }
1181 else
1182 {
1183 if(mT1>mT2) sr_switch = 1; // SR2
1184 else sr_switch = 0; // SR1
1185 }
1186 }
1187 if(preparedInput.m_LFVmode==1) // mu+tau(->e) case
1188 {
1189 if(preparedInput.m_vistau1.M()>0.05 && preparedInput.m_vistau2.M()<0.05)
1190 {
1191 if(mT1>mT2) sr_switch = 0; // SR1
1192 else sr_switch = 1; // SR2
1193 }
1194 else
1195 {
1196 if(mT1>mT2) sr_switch = 1; // SR2
1197 else sr_switch = 0; // SR1
1198 }
1199 }
1200 }
1201 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
1202 {
1203 if(preparedInput.m_vistau1.M()<0.12 && preparedInput.m_vistau2.M()>0.12)
1204 {
1205 if(mT1>mT2) sr_switch = 0; // SR1
1206 else sr_switch = 1; // SR2
1207 }
1208 else
1209 {
1210 if(mT1>mT2) sr_switch = 1; // SR2
1211 else sr_switch = 0; // SR1
1212 }
1213 }
1214
1215 m_UseHT = false;
1216 if(preparedInput.m_Njet25==0) // 0-jet
1217 {
1218 double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1219 double METresScale;
1220 double METoffset;
1221 if(sr_switch==0)
1222 {
1223 METresScale = 0.41*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1224 METoffset = 7.36*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1225 }
1226 else
1227 {
1228 METresScale = 0.34*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1229 METoffset = 6.61*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1230 }
1231 if(preparedInput.m_fUseVerbose==1) {
1232 Info("DiTauMassTools", "%s", ("SumEt = "+std::to_string(preparedInput.m_SumEt)).c_str());
1233 Info("DiTauMassTools", "%s", ("METoffset = "+std::to_string(METoffset)).c_str());
1234 Info("DiTauMassTools", "%s", ("METresScale = "+std::to_string(METresScale)).c_str());
1235 }
1236
1237 double sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1238 preparedInput.m_METsigmaP = sigma;
1239 preparedInput.m_METsigmaL = sigma;
1240 if(preparedInput.m_fUseVerbose==1) {
1241 Info("DiTauMassTools", "%s", ("=> METsigmaP = "+std::to_string(preparedInput.m_METsigmaP)).c_str());
1242 Info("DiTauMassTools", "%s", ("=> METsigmaL = "+std::to_string(preparedInput.m_METsigmaL)).c_str());
1243 }
1244 }
1245 if(preparedInput.m_Njet25>0) // Inclusive 1-jet
1246 {
1247 double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1248 double sigma = 0.;
1249 double METresScale;
1250 double METoffset;
1251 if(sr_switch==0)
1252 {
1253 METresScale = 0.38*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1254 METoffset = 7.96*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1255 }
1256 else
1257 {
1258 METresScale = 0.39*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1259 METoffset = 6.61*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1260 }
1261
1262 // MET resolution can't be perfect in presence of other objects (i.e., electrons, jets, taus), so assume minSumEt = 5.0 for now
1263 sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1264 preparedInput.m_METsigmaP = sigma;
1265 preparedInput.m_METsigmaL = sigma;
1266 } // Njet25>0
1267 }
1268 else //LFV
1269 {
1270 //DRDRMERGE end addition
1271
1272 if(preparedInput.m_METScanScheme==1) // default for Winter 2012 and further
1273 {
1274 //LEP-HAD
1275 if ( preparedInput.m_tauTypes==TauTypes::lh ) // lephad case
1276 {
1277 //0-jet
1278 if(preparedInput.m_Njet25==0)//0-jet
1279 {
1280 // placeholder for 2019 tune
1281 if (m_mmcCalibrationSet==MMCCalibrationSet::MMC2019 ||
1282 m_mmcCalibrationSet==MMCCalibrationSet::MMC2024){
1283 if(preparedInput.m_MetVec.R()<20.0) // 0-jet low MET case
1284 {
1285 if(std::abs(preparedInput.m_DelPhiTT)>2.95 && m_allowUseHT) // use mHt only if dPhi(lep-tau)>2.95
1286 {
1287 m_UseHT = true;
1288 // giving priority to external settings
1289 if(preparedInput.m_MHtSigma1<0.0) preparedInput.m_MHtSigma1 = 4.822;
1290 if(preparedInput.m_MHtSigma2<0.0) preparedInput.m_MHtSigma2 = 10.31;
1291 if(preparedInput.m_MHtGaussFr<0.0) preparedInput.m_MHtGaussFr = 6.34E-5;
1292 }
1293 else
1294 {
1295 m_UseHT = false;
1296 double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1297 double METresScale = 0.32*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1298 double METoffset = 5.38*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1299 double sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1300 preparedInput.m_METsigmaP = sigma;
1301 preparedInput.m_METsigmaL = sigma;
1302 }
1303 }
1304 else // 0-jet high MET case
1305 {
1306 if(std::abs(preparedInput.m_DelPhiTT)>2.8 && m_allowUseHT) // use mHt only if dPhi(lep-tau)>2.8
1307 {
1308 m_UseHT = true;
1309 // giving priority to external settings
1310 if(preparedInput.m_MHtSigma1<0.0) preparedInput.m_MHtSigma1 = 7.5;
1311 if(preparedInput.m_MHtSigma2<0.0) preparedInput.m_MHtSigma2 = 13.51;
1312 if(preparedInput.m_MHtGaussFr<0.0) preparedInput.m_MHtGaussFr = 6.81E-4;
1313 preparedInput.m_METsigmaP = preparedInput.m_MHtSigma2; // sigma of 2nd Gaussian for missing Ht resolution
1314 preparedInput.m_METsigmaL = preparedInput.m_MHtSigma2;
1315 }
1316 else
1317 {
1318 m_UseHT = false;
1319 double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1320 double METresScale = 0.87*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1321 double METoffset = 4.16*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1322 double sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1323 preparedInput.m_METsigmaP = sigma;
1324 preparedInput.m_METsigmaL = sigma;
1325 }
1326 } // high MET
1327 } // MMC2019
1328 // 2015 high-mass tune; avergare MET resolution for Mh=600,1000 mass points
1329 else if (m_mmcCalibrationSet==MMCCalibrationSet::MMC2015HIGHMASS)
1330 {
1331 m_UseHT = false;
1332 double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1333 double METresScale = 0.65*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1334 double METoffset = 5.0*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1335 double sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1336 preparedInput.m_METsigmaP = sigma;
1337 preparedInput.m_METsigmaL = sigma;
1338 } // MMC2015HIGHMASS
1339 } // 0 jet
1340 //1-jet
1341 else if(preparedInput.m_Njet25>0) // Inclusive 1-jet and VBF lep-had categories for Winter 2012
1342 {
1343 double sigmaSyst=0.10; // 10% systematics for now (be conservative)
1344 double sigma=0.;
1345 // 2015 high-mass tune; average MET resolution for Mh=400,600 mass points (they look consistent);
1346 if (m_mmcCalibrationSet==MMCCalibrationSet::MMC2015HIGHMASS)
1347 {
1348 double METresScale=0.86*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1349 double METoffset=3.0*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1350 // MET resolution can't be perfect in presence of other objects (i.e., electrons, jets, taus), so assume minSumEt=5.0 for now
1351 sigma= preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1352 }
1353 //2019
1354 else if (m_mmcCalibrationSet==MMCCalibrationSet::MMC2019 ||
1355 m_mmcCalibrationSet==MMCCalibrationSet::MMC2024)
1356 {
1357 double x = preparedInput.m_DelPhiTT;
1358 double dphi_scale = x > 0.3 ? 0.9429 - 0.059*x + 0.054*x*x : 0.728;
1359 double METoffset = 1.875*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1360 double METresScale1 = 8.914*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1361 double METresScale2 = -8.53*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1362
1363
1364 sigma = preparedInput.m_SumEt > 80.0 ? METoffset + METresScale1*TMath::Log(sqrt(preparedInput.m_SumEt)+METresScale2) : 5.0;
1365 sigma = sigma * dphi_scale;
1366 }
1367 //
1368
1369 preparedInput.m_METsigmaP=sigma;
1370 preparedInput.m_METsigmaL=sigma;
1371 } // Njet25>0
1372
1373 } // lep-had
1374
1375 //HAD-HAD
1376 else if(preparedInput.m_tauTypes==TauTypes::hh) // had-had events
1377 {
1378 if(preparedInput.m_Njet25==0 && m_mmcCalibrationSet==MMCCalibrationSet::MMC2015HIGHMASS) //0-jet high mass hadhad
1379 {
1380 // 2015 high-mass tune; average of all mass points
1381 // double METresScale=-1.;
1382 // double METoffset=-1.;
1383 double sigmaSyst=0.10; // 10% systematics for now (be conservative)
1384
1385 double METresScale=0.9*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1386 double METoffset=-1.8*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1387 double sigma= preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : std::abs(METoffset);
1388 preparedInput.m_METsigmaP=sigma;
1389 preparedInput.m_METsigmaL=sigma;
1390
1391 }
1392 else if(preparedInput.m_Njet25==0 &&
1393 (m_mmcCalibrationSet==MMCCalibrationSet::MMC2019 ||
1394 m_mmcCalibrationSet==MMCCalibrationSet::MMC2024))
1395 {
1396 double sigmaSyst=0.10; // 10% systematics for now (be conservative)
1397 double x = preparedInput.m_DelPhiTT;
1398 double dphi_scale = x > 2.5 ? 11.0796 - 4.61236*x + 0.423617*x*x : 2.;
1399 double METoffset = -8.51013*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1400 double METresScale1 = 8.54378*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1401 double METresScale2 = -3.97146*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1402 double sigma= preparedInput.m_SumEt>80.0 ? METoffset+METresScale1*TMath::Log(sqrt(preparedInput.m_SumEt)+METresScale2) : 5.;
1403 sigma = sigma*dphi_scale;
1404
1405 preparedInput.m_METsigmaP=sigma;
1406 preparedInput.m_METsigmaL=sigma;
1407
1408 }
1409 else if(preparedInput.m_Njet25==0 && m_allowUseHT) // 0-jet high MET had-had category for Winter 2012
1410 {
1411
1412 m_UseHT=true; // uncomment this line to enable HT also for HH (crucial)
1413 // updated for final cuts, may 21 2012
1414 if(preparedInput.m_MHtSigma1<0.0) preparedInput.m_MHtSigma1=5.972;
1415 if(preparedInput.m_MHtSigma2<0.0) preparedInput.m_MHtSigma2=13.85;
1416 // if(MHtGaussFr<0.0) MHtGaussFr=0.4767; // don't directly use 2nd fraction
1417 }
1418 //1-jet
1419 else // Inclusive 1-jet and VBF categories
1420 {
1421 double METresScale=-1.;
1422 double METoffset=-1.;
1423 double sigmaSyst=0.10; // 10% systematics for now (be conservative)
1424
1425 // previous value in trunk
1426 if(m_mmcCalibrationSet==MMCCalibrationSet::MMC2015HIGHMASS) {
1427 // 2015 high-mass tune; average of all mass points
1428 METresScale = 1.1*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1429 METoffset = -5.0*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1430 }
1431 // MET resolution can't be perfect in presence of other objects (i.e., electrons, jets, taus), so assume minSumEt=5.0 for now
1432 double sigma = preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : std::abs(METoffset);
1433
1434 if(m_mmcCalibrationSet==MMCCalibrationSet::MMC2019 ||
1435 m_mmcCalibrationSet==MMCCalibrationSet::MMC2024) {
1436 double x = preparedInput.m_DelPhiTT;
1437 double dphi_scale = x > 0.6 ? 1.42047 - 0.666644*x + 0.199986*x*x : 1.02;
1438 METoffset = 1.19769*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1439 double METresScale1 = 5.61687*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1440 double METresScale2 = -4.2076*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1441 sigma= preparedInput.m_SumEt>115.0 ? METoffset+METresScale1*TMath::Log(sqrt(preparedInput.m_SumEt)+METresScale2) : 12.1;
1442 sigma = sigma*dphi_scale;
1443 } //for hh 2016 mc15c
1444
1445 preparedInput.m_METsigmaP = sigma;
1446 preparedInput.m_METsigmaL = sigma;
1447
1448 }// 1 jet
1449 } // had-had
1450 //LEP-LEP
1451 else if(preparedInput.m_tauTypes==TauTypes::ll) // setup for LEP-LEP channel
1452 {
1453 if(m_mmcCalibrationSet==MMCCalibrationSet::MMC2015HIGHMASS) // placeholder for 2015 high-mass tune; for now it is the same as 2012
1454 {
1455 m_UseHT = false;
1456 double sigmaSyst = 0.10; // 10% systematics for now (be conservative)
1457 double METresScale = -1.0;
1458 double METoffset = -1.0;
1459 double sigma = 5.0;
1460 // tune is based on cuts for Run-1 paper analysis
1461 if(preparedInput.m_Njet25==0)
1462 {
1463 // use tune for emebedding
1464 METresScale=-0.4307*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1465 METoffset=7.06*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1466 double METresScale2=0.07693*(1.0+preparedInput.m_METresSyst*sigmaSyst); // quadratic term
1467 // this is a tune for Higgs125
1468 // METresScale=-0.5355*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1469 // METoffset=11.5*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1470 // double METresScale2=0.07196*(1.0+preparedInput.m_METresSyst*sigmaSyst); // quadratic term
1471 sigma= preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt)+METresScale2*preparedInput.m_SumEt : METoffset;
1472 }
1473 if(preparedInput.m_Njet25>0)
1474 {
1475 // use tune for embedding
1476 METresScale=0.8149*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1477 METoffset=5.343*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1478 // this is a tune for Higgs125
1479 // METresScale=0.599*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1480 // METoffset=8.223*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1481 sigma= preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt) : METoffset;
1482 }
1483 preparedInput.m_METsigmaP = sigma;
1484 preparedInput.m_METsigmaL = sigma;
1485 } // end of MMC2015HIGHMASS
1486
1487 if(m_mmcCalibrationSet==MMCCalibrationSet::MMC2019 ||
1488 m_mmcCalibrationSet==MMCCalibrationSet::MMC2024)
1489 // 2019 leplep
1490 {
1491 m_UseHT=false;
1492 double sigmaSyst=0.10; // 10% systematics for now (be conservative)
1493 double METresScale=-1.0;
1494 double METoffset=-1.0;
1495 double sigma=5.0;
1496 double min_sigma = 2.0;
1497 // tune is based on cuts for Run-1 paper analysis
1498 if(preparedInput.m_Njet25==0)
1499 {
1500 // Madgraph Ztautau MET param
1501 METoffset = 4.22581*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1502 METresScale = 0.03818*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1503 double METresScale2= 0.12623;
1504 sigma= preparedInput.m_SumEt>0.0 ? METoffset+METresScale*sqrt(preparedInput.m_SumEt)+METresScale2*preparedInput.m_SumEt : min_sigma;
1505 if (m_fUseDphiLL) {
1506 double p0 = 2.60131;
1507 double p1const = 1.22427;
1508 double p2quad = -1.71261;
1509 double DphiLL = std::abs(Phi_mpi_pi(preparedInput.m_vistau1.Phi()-preparedInput.m_vistau2.Phi()));
1510 sigma *= (DphiLL < p0) ? p1const : p1const+
1511 p2quad*p0*p0 - 2*p2quad*p0*DphiLL+p2quad*DphiLL*DphiLL;
1512 }
1513 if (sigma < min_sigma) sigma = min_sigma;
1514 }
1515 if(preparedInput.m_Njet25>0)
1516 {
1517 // Madgraph Ztautau MET param
1518 METoffset = 5.42506*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1519 METresScale = 5.36760*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1520 double METoffset2 = -4.86808*(1.0+preparedInput.m_METresSyst*sigmaSyst);
1521 if (preparedInput.m_SumEt > 0.0) {
1522 double x = sqrt(preparedInput.m_SumEt);
1523 sigma = (x+METoffset2 > 1) ? METoffset+METresScale*log(x+METoffset2) : METoffset;
1524 } else {
1525 sigma = METoffset;
1526 }
1527 if (m_fUseDphiLL) {
1528 double p0 = 2.24786;
1529 double p1const = 0.908597;
1530 double p2quad = 0.544577;
1531 double DphiLL = std::abs(Phi_mpi_pi(preparedInput.m_vistau1.Phi()-preparedInput.m_vistau2.Phi()));
1532 sigma *= (DphiLL < p0) ? p1const : p1const+
1533 p2quad*p0*p0 - 2*p2quad*p0*DphiLL+p2quad*DphiLL*DphiLL;
1534 }
1535 if (sigma < min_sigma) sigma = min_sigma;
1536 }
1537 preparedInput.m_METsigmaP=sigma;
1538 preparedInput.m_METsigmaL=sigma;
1539 }//2016 mc15c
1540
1541 } // lep-lep
1542
1543 } //preparedInput.METScanScheme
1544
1545 if(preparedInput.m_METScanScheme==0) // old scheme with JER
1546 {
1547 if(preparedInput.m_dataType==0 || preparedInput.m_dataType==1) preparedInput.SetMetScanParamsUE(preparedInput.m_SumEt,preparedInput.m_METcovphi,preparedInput.m_dataType);
1548 else preparedInput.SetMetScanParamsUE(preparedInput.m_SumEt,preparedInput.m_METcovphi,0);
1549 }
1550 }
1551 } // end else LFV
1552
1553 return;
1554}
Phi_mpi_pi
__HOSTDEV__ double Phi_mpi_pi(double)
Definition GeoRegion.cxx:10
x
#define x
DiTauMassTools::MissingMassInput::SetMetScanParamsUE
void SetMetScanParamsUE(double sumEt, double phi_scan=0.0, int data_code=0)
Definition MissingMassInput.cxx:89
DiTauMassTools::mT
double mT(const VectorType &vec, const XYVector &met_vec)
Definition PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:72
TRTCalib_cfilter.p0
p0
Definition TRTCalib_cfilter.py:129
ZDCMsg::Info
@ Info
Definition ZDCMsg.h:20
pdg_comparison.sigma
sigma
Definition pdg_comparison.py:324

◆ mEtAndTauProbability()

double MissingMassProb::mEtAndTauProbability ( MissingMassInput & preparedInput)

Definition at line 580 of file MissingMassProb.cxx.

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}
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
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::MissingMassProb::MetProbability
double MetProbability(MissingMassInput &preparedInput, const double &met1, const double &met2, const double &MetSigma1, const double &MetSigma2)
Definition MissingMassProb.cxx:559

◆ mEtAndTauProbabilityWrapper()

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 )
static

Definition at line 39 of file MissingMassProb.cxx.

39 {
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}

◆ MetProbability()

double MissingMassProb::MetProbability ( MissingMassInput & preparedInput,
const double & met1,
const double & met2,
const double & MetSigma1,
const double & MetSigma2 )

Definition at line 559 of file MissingMassProb.cxx.

559 {
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}

◆ MetProbabilityWrapper()

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 )
static

Definition at line 29 of file MissingMassProb.cxx.

29 {
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}

◆ MHtProbability()

double MissingMassProb::MHtProbability ( MissingMassInput & preparedInput,
const double & d_mhtX,
const double & d_mhtY,
const double & mht,
const double & trueMetGuess,
const double & mht_offset )

Definition at line 1003 of file MissingMassProb.cxx.

1004 {
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}

◆ MHtProbabilityHH()

double MissingMassProb::MHtProbabilityHH ( MissingMassInput & preparedInput,
const double & d_mhtX,
const double & d_mhtY,
const double & mht,
const double & trueMetGuess,
const double & mht_offset )

Definition at line 1020 of file MissingMassProb.cxx.

1021 {
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}

◆ MnuProbability() [1/2]

double MissingMassProb::MnuProbability ( MissingMassInput & preparedInput,
double mnu )

Definition at line 976 of file MissingMassProb.cxx.

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}

◆ MnuProbability() [2/2]

double MissingMassProb::MnuProbability ( MissingMassInput & preparedInput,
double mnu,
double binsize )

Definition at line 945 of file MissingMassProb.cxx.

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}

◆ MnuProbabilityNewWrapper()

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 )
static

Definition at line 136 of file MissingMassProb.cxx.

136 {
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}

◆ MnuProbabilityWrapper()

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 )
static

Definition at line 78 of file MissingMassProb.cxx.

78 {
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}

◆ myDelThetaHadFunc()

double MissingMassProb::myDelThetaHadFunc ( double * x,
double * par )

Definition at line 1090 of file MissingMassProb.cxx.

1091{
1092 double fitval=1.0E-10;
1093 if(x[0]>TMath::Pi() || x[0]<0.0) return fitval;
1094 const double arg=x[0];
1095 const double arg_L=arg;
1096 const double mean=par[1];
1097 const double sigmaG=par[2];
1098 const double mpv=par[3];
1099 const double sigmaL=par[4];
1100
1101 if (m_mmcCalibrationSet==MMCCalibrationSet::MMC2015HIGHMASS
1102 || m_mmcCalibrationSet==MMCCalibrationSet::UPGRADE
1103 || m_mmcCalibrationSet==MMCCalibrationSet::LFVMMC2012){
1104 const double norm=sqrt(2.0*TMath::Pi());
1105 const double g1=TMath::Gaus(arg,mean,sigmaG)/norm;
1106 const double g2=TMath::Landau(arg_L,mpv,sigmaL)/norm;
1107 fitval=par[0]*g1/sigmaG+par[5]*g2/sigmaL;
1108 }
1109
1110 if(fitval<0.0) return 0.0;
1111 return fitval;
1112}
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
create_dcsc_inputs_sqlite.arg
list arg
Definition create_dcsc_inputs_sqlite.py:48

◆ myDelThetaLepFunc()

double MissingMassProb::myDelThetaLepFunc ( double * x,
double * par )

Definition at line 1115 of file MissingMassProb.cxx.

1116{
1117 double fitval=1.0E-10;
1118 if(x[0]>TMath::Pi() || x[0]<0.0) return fitval;
1119 double arg=x[0]/par[1];
1120
1121 double normL=par[5];
1122 if(normL<0.0) normL=0.0;
1123
1124 if(arg<1) arg=sqrt(std::abs(arg));
1125 else arg=arg*arg;
1126 const double arg_L=x[0];
1127 const double mean=1.0;
1128 const double sigmaG=par[2];
1129 const double mpv=par[3];
1130 const double sigmaL=par[4];
1131 const double g1=normL*TMath::Gaus(arg,mean,sigmaG);
1132 const double g2=TMath::Landau(arg_L,mpv,sigmaL);
1133 fitval=par[0]*(g1+g2)/(1.0+normL);
1134 if(fitval<0.0) return 0.0;
1135 return fitval;
1136}

◆ SetAllowUseHT()

void DiTauMassTools::MissingMassProb::SetAllowUseHT ( bool allowUseHT)
inline

Definition at line 45 of file MissingMassProb.h.

45{ m_allowUseHT=allowUseHT;}

◆ setParamAngle()

void MissingMassProb::setParamAngle ( const PtEtaPhiMVector & tauvec,
int tau,
int tautype )

Definition at line 161 of file MissingMassProb.cxx.

161 {
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}

◆ setParamNuMass()

void MissingMassProb::setParamNuMass ( )

Definition at line 153 of file MissingMassProb.cxx.

153 {
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}

◆ setParamRatio()

void MissingMassProb::setParamRatio ( int tau,
int tautype )

Definition at line 190 of file MissingMassProb.cxx.

190 {
191 int type = tautype;
192 if(tautype > 4 && tautype < 8) type = 4;
193 if (tautype <= 5){
194 if(tau==1){
195 m_formulaRatio1 = m_formulaRatioHad1;
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 {
200 m_formulaRatio2 = m_formulaRatioHad2;
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){
207 m_formulaRatio1 = m_formulaRatioLep1;
208 for(int i=0; i<m_formulaRatio1->GetNpar(); i++){
209 m_formulaRatio1->SetParameter(i, m_paramVectorRatioLep[0]->GetParameter(i));
210 }
211 } else {
212 m_formulaRatio2 = m_formulaRatioLep2;
213 for(int i=0; i<m_formulaRatio2->GetNpar(); i++){
214 m_formulaRatio2->SetParameter(i, m_paramVectorRatioLep[0]->GetParameter(i));
215 }
216 }
217 }
218}

◆ SetUseDphiLL()

void DiTauMassTools::MissingMassProb::SetUseDphiLL ( bool val)
inline

Definition at line 57 of file MissingMassProb.h.

57{m_fUseDphiLL = val;}
Pythia8_RapidityOrderMPI.val
val
Definition Pythia8_RapidityOrderMPI.py:14

◆ SetUseHT()

void DiTauMassTools::MissingMassProb::SetUseHT ( bool val)
inline

Definition at line 48 of file MissingMassProb.h.

48{ m_UseHT=val; }

◆ SetUseMnuProbability()

void DiTauMassTools::MissingMassProb::SetUseMnuProbability ( bool val)
inline

Definition at line 54 of file MissingMassProb.h.

54{ m_fUseMnuProbability=val; }

◆ SetUseTauProbability()

void DiTauMassTools::MissingMassProb::SetUseTauProbability ( bool val)
inline

Definition at line 51 of file MissingMassProb.h.

51{ m_fUseTauProbability=val; }

◆ TauProbability() [1/2]

double MissingMassProb::TauProbability ( MissingMassInput & preparedInput,
const int & type1,
const PtEtaPhiMVector & vis1,
const PtEtaPhiMVector & nu1,
const int & type2,
const PtEtaPhiMVector & vis2,
const PtEtaPhiMVector & nu2 )

Definition at line 657 of file MissingMassProb.cxx.

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}
GEV
#define GEV
Definition PrintPhotonSF.cxx:25
ParticleConstants::PDG2011::tauMassInMeV
constexpr double tauMassInMeV
the mass of the tau (in MeV)
Definition ParticleConstants.h:32

◆ TauProbability() [2/2]

double MissingMassProb::TauProbability ( MissingMassInput & preparedInput,
const int & type1,
const PtEtaPhiMVector & vis1,
const PtEtaPhiMVector & nu1,
const int & type2,
const PtEtaPhiMVector & vis2,
const PtEtaPhiMVector & nu2,
const double & detmet )

Definition at line 730 of file MissingMassProb.cxx.

731 {
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}
scale2
#define scale2
Definition JetAttributeHisto.cxx:42
scale1
#define scale1
Definition JetAttributeHisto.cxx:41
P
static Double_t P(Double_t *tt, Double_t *par)
Definition LArPhysWaveHECTool.cxx:280
AnalysisUtils::Delta::R
double R(const INavigable4Momentum *p1, const double v_eta, const double v_phi)
Definition AnalysisMisc.h:49

◆ TauProbabilityLFV()

double MissingMassProb::TauProbabilityLFV ( MissingMassInput & preparedInput,
const int & type1,
const PtEtaPhiMVector & vis1,
const PtEtaPhiMVector & nu1 )

Definition at line 620 of file MissingMassProb.cxx.

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}

◆ TauProbabilityNewWrapper()

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 )
static

Definition at line 122 of file MissingMassProb.cxx.

122 {
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}

◆ TauProbabilityWrapper()

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 )
static

Definition at line 70 of file MissingMassProb.cxx.

70 {
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}

Member Data Documentation

◆ m_allowUseHT

bool DiTauMassTools::MissingMassProb::m_allowUseHT
private

Definition at line 120 of file MissingMassProb.h.

◆ m_formulaAngle1

TF1* DiTauMassTools::MissingMassProb::m_formulaAngle1 = new TF1("formulaAngle1", "[0]*exp(-[2]*(log((x+[3])/[1]))**2)")
private

Definition at line 101 of file MissingMassProb.h.

◆ m_formulaAngle2

TF1* DiTauMassTools::MissingMassProb::m_formulaAngle2 = new TF1("formulaAngle2", "[0]*exp(-[2]*(log((x+[3])/[1]))**2)")
private

Definition at line 102 of file MissingMassProb.h.

◆ m_formulaNuMass

TF1* DiTauMassTools::MissingMassProb::m_formulaNuMass = new TF1("formulaNuMass", "pol6")
private

Definition at line 109 of file MissingMassProb.h.

◆ m_formulaRatio1

TF1* DiTauMassTools::MissingMassProb::m_formulaRatio1
private

Definition at line 103 of file MissingMassProb.h.

◆ m_formulaRatio2

TF1* DiTauMassTools::MissingMassProb::m_formulaRatio2
private

Definition at line 104 of file MissingMassProb.h.

◆ m_formulaRatioHad1

TF1* DiTauMassTools::MissingMassProb::m_formulaRatioHad1 = new TF1("formulaRatio1", "gaus(0)")
private

Definition at line 107 of file MissingMassProb.h.

◆ m_formulaRatioHad2

TF1* DiTauMassTools::MissingMassProb::m_formulaRatioHad2 = new TF1("formulaRatio2", "gaus(0)")
private

Definition at line 108 of file MissingMassProb.h.

◆ m_formulaRatioLep1

TF1* DiTauMassTools::MissingMassProb::m_formulaRatioLep1 = new TF1("formulaRatio1", "gaus(0)+expo(3)")
private

Definition at line 105 of file MissingMassProb.h.

◆ m_formulaRatioLep2

TF1* DiTauMassTools::MissingMassProb::m_formulaRatioLep2 = new TF1("formulaRatio2", "gaus(0)+expo(3)")
private

Definition at line 106 of file MissingMassProb.h.

◆ m_fParams

TFile* DiTauMassTools::MissingMassProb::m_fParams
private

Definition at line 117 of file MissingMassProb.h.

◆ m_fUseDphiLL

bool DiTauMassTools::MissingMassProb::m_fUseDphiLL
private

Definition at line 124 of file MissingMassProb.h.

◆ m_fUseMnuProbability

bool DiTauMassTools::MissingMassProb::m_fUseMnuProbability
private

Definition at line 123 of file MissingMassProb.h.

◆ m_fUseTauProbability

bool DiTauMassTools::MissingMassProb::m_fUseTauProbability
private

Definition at line 122 of file MissingMassProb.h.

◆ m_mmcCalibrationSet

MMCCalibrationSet::e DiTauMassTools::MissingMassProb::m_mmcCalibrationSet
private

Definition at line 118 of file MissingMassProb.h.

◆ m_paramFilePath

std::string DiTauMassTools::MissingMassProb::m_paramFilePath
private

Definition at line 116 of file MissingMassProb.h.

◆ m_paramVectorAngle

std::vector<TF1*> DiTauMassTools::MissingMassProb::m_paramVectorAngle
private

Definition at line 110 of file MissingMassProb.h.

◆ m_paramVectorAngleLep

std::vector<TF1*> DiTauMassTools::MissingMassProb::m_paramVectorAngleLep
private

Definition at line 111 of file MissingMassProb.h.

◆ m_paramVectorNuMass

std::vector<TF1*> DiTauMassTools::MissingMassProb::m_paramVectorNuMass
private

Definition at line 114 of file MissingMassProb.h.

◆ m_paramVectorRatio

std::vector<TF1*> DiTauMassTools::MissingMassProb::m_paramVectorRatio
private

Definition at line 112 of file MissingMassProb.h.

◆ m_paramVectorRatioLep

std::vector<TF1*> DiTauMassTools::MissingMassProb::m_paramVectorRatioLep
private

Definition at line 113 of file MissingMassProb.h.

◆ 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)> > DiTauMassTools::MissingMassProb::m_probListConstant

Definition at line 92 of file MissingMassProb.h.

◆ 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)> > DiTauMassTools::MissingMassProb::m_probListOneTau

Definition at line 93 of file MissingMassProb.h.

◆ 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)> > DiTauMassTools::MissingMassProb::m_probListTwoTau

Definition at line 94 of file MissingMassProb.h.

◆ m_UseHT

bool DiTauMassTools::MissingMassProb::m_UseHT
private

Definition at line 121 of file MissingMassProb.h.

◆ s_fit_param

thread_local double MissingMassProb::s_fit_param
staticprivate

Definition at line 97 of file MissingMassProb.h.

◆ s_ter_sigma_par

thread_local double MissingMassProb::s_ter_sigma_par
staticprivate

Definition at line 98 of file MissingMassProb.h.

The documentation for this class was generated from the following files:
Generated on for ATLAS Offline Software by doxygen 1.14.0