ATLAS Offline Software
MissingMassCalculator.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 // vim: ts=8 sw=2
6 /*
7  Missing Mass Calculator
8 */
9 //
10 // to be done : tau 4-vect and type should be data member of MMC.
11 
12 // if histogram smoothing
13 //#define SMOOTH
14 
15 #include "DiTauMassTools/MissingMassCalculator.h" // this is for RootCore package
16 #include <fstream>
17 #include <iomanip>
18 #include <iostream>
19 #include <sstream>
20 // #include "MissingMassCalculator.h" // this is for standalone
21 // package
22 
23 #include <TObject.h>
24 // SpeedUp committed from revision 163876
25 #include <TF1.h>
26 #include <TFitResult.h>
27 #include <TFitResultPtr.h>
28 #include <TMatrixDSym.h>
29 #include <TObject.h>
30 #include <TVectorD.h>
31 #include "Math/VectorUtil.h"
32 
34 
35 namespace {
36  constexpr double GEV = 1000.0;
37 }
38 
39 
40 using namespace DiTauMassTools;
41 using ROOT::Math::PtEtaPhiMVector;
42 using ROOT::Math::PxPyPzMVector;
43 using ROOT::Math::XYVector;
46 
47 //______________________________constructor________________________________
49  MMCCalibrationSet::e aset, std::string paramFilePath)
50  : m_randomGen(), Prob(new MissingMassProb(aset, paramFilePath)) {
51  m_mmcCalibrationSet = aset;
53  preparedInput.m_beamEnergy = 6500.0; // for now LHC default is sqrt(S)=7 TeV
54  m_niter_fit1 = 20;
55  m_niter_fit2 = 30;
56  m_niter_fit3 = 10;
57  m_NsucStop = -1;
58  m_NiterRandom = -1; // if the user does not set it to positive value,will be set
59  // in SpaceWalkerInit
60  m_niterRandomLocal = -1; // niterandom which is really used
61  // to be used with RMSSTOP NiterRandom=10000000; // number of random
62  // iterations for lh. Multiplied by 10 for ll, divided by 10 for hh (to be
63  // optimised)
64  // RMSStop=200;// Stop criteria depending of rms of histogram
65  m_reRunWithBestMET = false;
66  m_RMSStop = -1; // disable
67 
68  m_RndmSeedAltering = 0; // can be changed to re-compute with different random seed
69  m_dRmax_tau = 0.4; // changed from 0.2
70  m_nsigma_METscan = -1; // number of sigmas for MET scan
71  m_nsigma_METscan_ll = 3.0; // number of sigmas for MET scan
72  m_nsigma_METscan_lh = 3.0; // number of sigmas for MET scan
73  m_nsigma_METscan_hh = 4.0; // number of sigmas for MET scan (4 for hh 2013)
74  m_nsigma_METscan_lfv_ll = 5.0; // number of sigmas for MET scan (LFV leplep)
75  m_nsigma_METscan_lfv_lh = 5.0; // number of sigmas for MET scan (LFV lephad)
76 
77  m_meanbinStop = -1; // meanbin stopping criterion (-1 if not used)
78  m_proposalTryMEt = -1; // loop on METproposal disable // FIXME should be cleaner
79  m_ProposalTryPhi = -1; // loop on Phiproposal disable
80  m_ProposalTryMnu = -1; // loop on MNuProposal disable
81  m_ProposalTryEtau = -1; // loop on ETauProposal disable
82 
83  Prob->SetUseTauProbability(true); // TauProbability is ON by default DRMERGE comment out for now
84  Prob->SetUseMnuProbability(false); // MnuProbability is OFF by default
85  Prob->SetUseDphiLL(false); // added by Tomas Davidek for lep-lep
86  m_dTheta3d_binMin = 0.0025;
87  m_dTheta3d_binMax = 0.02;
88  preparedInput.m_METresSyst = 0; // no MET resolution systematics by default (+/-1: up/down 1 sigma)
89  preparedInput.m_dataType = 1; // set to "data" by default
90  preparedInput.m_fUseTailCleanup = 1; // cleanup by default for lep-had Moriond 2012 analysis
91  preparedInput.m_fUseDefaults = 0; // use pre-set defaults for various configurations; if set it to 0
92  // if need to study various options
93  m_fUseEfficiencyRecovery = 0; // no re-fit by default
98 
99  preparedInput.m_METScanScheme = 1; // MET-scan scheme: 0- use JER; 1- use simple sumEt & missingHt
100  // for Njet=0 events in (lep-had winter 2012)
101  // MnuScanRange=ParticleConstants::tauMassInMeV / GEV; // range of M(nunu) scan
102  m_MnuScanRange = 1.5; // better value (sacha)
103  preparedInput.m_LFVmode = -1; // by default consider case of H->mu+tau(->ele)
105 
107  m_iterTheta3d = 0;
108  m_debugThisIteration = false;
109  m_lfvLeplepRefit = true;
110  m_SaveLlhHisto = false;
111 
112  m_nsolmax = 4;
114 
115  m_nuvecsol1.resize(m_nsolmax);
116  m_nuvecsol2.resize(m_nsolmax);
117  m_tauvecsol1.resize(m_nsolmax);
118  m_tauvecsol2.resize(m_nsolmax);
119  m_tauvecprob1.resize(m_nsolmax);
120  m_tauvecprob2.resize(m_nsolmax);
121 
122  m_nsol = 0;
127 
128  m_nsolOld = 0;
133 
134  float hEmax = 3000.0; // maximum energy (GeV)
135  // number of bins
136  int hNbins = 1500; // original 2500 for mass, 10000 for P
137  // choice of hNbins also related to size of window for fitting (see
138  // maxFromHist)
139 
140  //--- define histograms for histogram method
141  //--- upper limits need to be revisied in the future!!! It may be not enough
142  // for some analyses
143 
144  m_fMfit_all = std::make_shared<TH1F>("MMC_h1", "M", hNbins, 0.0,
145  hEmax); // all solutions
146  m_fMfit_all->Sumw2(); // allow proper error bin calculation. Slightly slower but
147  // completely negligible
148 
149  // histogram without weight. useful for debugging. negligibly slow until now
151  std::make_shared<TH1F>("MMC_h1NoW", "M no weight", hNbins, 0.0, hEmax); // all solutions
152 
153  m_fPXfit1 = std::make_shared<TH1F>("MMC_h2", "Px1", 4 * hNbins, -hEmax,
154  hEmax); // Px for tau1
155  m_fPYfit1 = std::make_shared<TH1F>("MMC_h3", "Py1", 4 * hNbins, -hEmax,
156  hEmax); // Py for tau1
157  m_fPZfit1 = std::make_shared<TH1F>("MMC_h4", "Pz1", 4 * hNbins, -hEmax,
158  hEmax); // Pz for tau1
159  m_fPXfit2 = std::make_shared<TH1F>("MMC_h5", "Px2", 4 * hNbins, -hEmax,
160  hEmax); // Px for tau2
161  m_fPYfit2 = std::make_shared<TH1F>("MMC_h6", "Py2", 4 * hNbins, -hEmax,
162  hEmax); // Py for tau2
163  m_fPZfit2 = std::make_shared<TH1F>("MMC_h7", "Pz2", 4 * hNbins, -hEmax,
164  hEmax); // Pz for tau2
165 
166  m_fMfit_all->SetDirectory(0);
167 
168  m_fMfit_allNoWeight->SetDirectory(0);
169  m_fPXfit1->SetDirectory(0);
170  m_fPYfit1->SetDirectory(0);
171  m_fPZfit1->SetDirectory(0);
172  m_fPXfit2->SetDirectory(0);
173  m_fPYfit2->SetDirectory(0);
174  m_fPZfit2->SetDirectory(0);
175 
176  // max hist fitting function
177  m_fFitting =
178  new TF1("MMC_maxFitting", this, &MissingMassCalculator::maxFitting, 0., hEmax, 3);
179  // Sets initial parameter names
180  m_fFitting->SetParNames("Max", "Mean", "InvWidth2");
181 
182  if (preparedInput.m_fUseVerbose == 1) {
183  gDirectory->pwd();
184  gDirectory->ls();
185  }
186 
187  if (preparedInput.m_fUseVerbose == 1) {
188  gDirectory->pwd();
189  gDirectory->ls();
190  }
191 }
192 
194 
195 //_____________________________________________________________________________
196 // Main Method to run MissingMassCalculator
198  const xAOD::IParticle *part2,
199  const xAOD::MissingET *met,
200  const int &njets) {
201  m_reRunWithBestMET = false;
202 
204  if (preparedInput.m_fUseVerbose == 1) {
205  Info("DiTauMassTools", "------------- Raw Input for MissingMassCalculator --------------");
206  }
207  FinalizeSettings(part1, part2, met, njets); // rawInput, preparedInput );
209  if (preparedInput.m_fUseVerbose == 1) {
210  Info("DiTauMassTools", "------------- Prepared Input for MissingMassCalculator--------------");
212  }
213 
214  if (preparedInput.m_LFVmode < 0) {
215  // remove argument DiTauMassCalculatorV9Walk work directly on preparedInput
217 
218  // re-running MMC for on failed events
220  // most events where MMC failed happened to have dPhi>2.9. Run re-fit only
221  // on these events
222  if (preparedInput.m_DelPhiTT > 2.9) {
223  // preparedInput.MetVec.Set(-(preparedInput.vistau1+preparedInput.vistau2).Px(),-(preparedInput.vistau1+preparedInput.vistau2).Py());
224  // // replace MET by MPT
225 
226  XYVector dummy_met(-(preparedInput.m_vistau1 + preparedInput.m_vistau2).Px(),
228  preparedInput.m_METcovphi = dummy_met.Phi();
229  double dummy_METres =
232  dummy_METres * std::abs(cos(dummy_met.Phi() - preparedInput.m_MetVec.Phi()));
234  dummy_METres * std::abs(sin(dummy_met.Phi() - preparedInput.m_MetVec.Phi()));
235  if (preparedInput.m_METsigmaP < 5.0)
237  m_nsigma_METscan_lh = 6.0; // increase range of MET scan
238  m_nsigma_METscan_hh = 6.0; // increase range of MET scan
239 
240  OutputInfo.ClearOutput(preparedInput.m_fUseVerbose); // clear output stuff before re-running
241  OutputInfo.m_FitStatus = DitauMassCalculatorV9walk(); // run MMC again
242  }
243  }
244 
245  }
246 
247  // running MMC in LFV mode for reconstructing mass of X->lep+tau
248  else {
249  if (preparedInput.m_fUseVerbose == 1) {
250  Info("DiTauMassTools", "Calling DitauMassCalculatorV9lfv");
251  }
253  }
254 
255  if(m_SaveLlhHisto){
256  TFile *outFile = TFile::Open("MMC_likelihoods.root", "UPDATE");
257  outFile->cd();
259  if (!outFile->GetDirectory(path.c_str()))
260  outFile->mkdir(path.c_str());
261  outFile->cd(path.c_str());
262  m_fMfit_all->Write(m_fMfit_all->GetName(), TObject::kOverwrite);
263  m_fMEtP_all->Write(m_fMEtP_all->GetName(), TObject::kOverwrite);
264  m_fMEtL_all->Write(m_fMEtL_all->GetName(), TObject::kOverwrite);
265  m_fMnu1_all->Write(m_fMnu1_all->GetName(), TObject::kOverwrite);
266  m_fMnu2_all->Write(m_fMnu2_all->GetName(), TObject::kOverwrite);
267  m_fPhi1_all->Write(m_fPhi1_all->GetName(), TObject::kOverwrite);
268  m_fPhi2_all->Write(m_fPhi2_all->GetName(), TObject::kOverwrite);
269  m_fMfit_allNoWeight->Write(m_fMfit_allNoWeight->GetName(), TObject::kOverwrite);
270  m_fMfit_allGraph->Write("Graph", TObject::kOverwrite);
271  TH1D *nosol = new TH1D("nosol", "nosol", 7, 0, 7);
272  nosol->SetBinContent(1, m_testptn1);
273  nosol->SetBinContent(2, m_testptn2);
274  nosol->SetBinContent(3, m_testdiscri1);
275  nosol->SetBinContent(4, m_testdiscri2);
276  nosol->SetBinContent(5, m_nosol1);
277  nosol->SetBinContent(6, m_nosol1);
278  nosol->SetBinContent(7, m_iterNuPV3);
279  nosol->Write(nosol->GetName(), TObject::kOverwrite);
280  outFile->Write();
281  outFile->Close();
282  }
283 
284  DoOutputInfo();
285  PrintResults();
287  return 1;
288 }
289 
290 //-------- clearing ditau container
292  fStuff.Mditau_best = 0.0;
293  fStuff.Sign_best = 1.0E6;
294  fStuff.nutau1 = PtEtaPhiMVector(0., 0., 0., 0.);
295  fStuff.nutau2 = PtEtaPhiMVector(0., 0., 0., 0.);
296  fStuff.vistau1 = PtEtaPhiMVector(0., 0., 0., 0.);
297  fStuff.vistau2 = PtEtaPhiMVector(0., 0., 0., 0.);
298  fStuff.RMSoverMPV = 0.0;
299 
300  return;
301 }
302 
303 //---------------------------- Accessors to output parameters
304 //------------------------
305 // finalizes output information
307  if (OutputInfo.m_FitStatus > 0) {
308  if (preparedInput.m_fUseVerbose == 1) {
309  Info("DiTauMassTools", "Retrieving output from fDitauStuffFit");
310  }
311  // MAXW method : get from fDittauStuffFit
314  double q1 = (1. - 0.68) / 2.;
315  double q2 = 1. - q1;
316  double xq[2], yq[2];
317  xq[0] = q1;
318  xq[1] = q2;
319  m_fMfit_all->GetQuantiles(2, yq, xq);
331  XYVector metmaxw(OutputInfo.m_nuvec1[MMCFitMethod::MAXW].Px() +
336 
340 
341  PtEtaPhiMVector tlvdummy(0., 0., 0., 0.);
342  XYVector metdummy(0., 0.);
350 
351  // MLNU3P method : get from fDittauStuffHisto 4 momentum
365 
366  XYVector metmlnu3p(OutputInfo.m_nuvec1[MMCFitMethod::MLNU3P].Px() +
371 
373  }
374 
377  OutputInfo.m_NSolutions = m_fMfit_all->GetEntries();
378  OutputInfo.m_SumW = m_fMfit_all->GetSumOfWeights();
379 
380  //----------------- Check if input was re-ordered in FinalizeInputStuff() and
381  // restore the original order if needed
382  if (preparedInput.m_InputReorder == 1) {
383  PtEtaPhiMVector dummy_vec1(0.0, 0.0, 0.0, 0.0);
384  PtEtaPhiMVector dummy_vec2(0.0, 0.0, 0.0, 0.0);
385  for (int i = 0; i < 3; i++) {
386  // re-ordering neutrinos
387  dummy_vec1 = OutputInfo.m_nuvec1[i];
388  dummy_vec2 = OutputInfo.m_nuvec2[i];
389  OutputInfo.m_nuvec1[i] = dummy_vec2;
390  OutputInfo.m_nuvec2[i] = dummy_vec1;
391  // re-ordering tau's
392  dummy_vec1 = OutputInfo.m_objvec1[i];
393  dummy_vec2 = OutputInfo.m_objvec2[i];
394  OutputInfo.m_objvec1[i] = dummy_vec2;
395  OutputInfo.m_objvec2[i] = dummy_vec1;
396  }
397  }
398 
399  return;
400 }
401 
402 // Printout of final results
404  if (preparedInput.m_fUseVerbose != 1)
405  return;
406 
407  Info("DiTauMassTools",
408  ".........................Other input.....................................");
409  Info("DiTauMassTools", "%s",
410  ("Beam energy =" + std::to_string(preparedInput.m_beamEnergy) +
411  " sqrt(S) for collisions =" + std::to_string(2.0 * preparedInput.m_beamEnergy))
412  .c_str());
413  Info("DiTauMassTools", "%s",
414  ("CalibrationSet " + MMCCalibrationSet::name[m_mmcCalibrationSet])
415  .c_str());
416  Info("DiTauMassTools", "%s",
417  ("LFV mode " + std::to_string(preparedInput.m_LFVmode) + " seed=" + std::to_string(m_seed))
418  .c_str());
419  Info("DiTauMassTools", "%s", ("usetauProbability=" + std::to_string(Prob->GetUseTauProbability()) +
420  " useTailCleanup=" + std::to_string(preparedInput.m_fUseTailCleanup))
421  .c_str());
422 
423  if (preparedInput.m_InputReorder != 0) {
424  Info("DiTauMassTools",
425  "tau1 and tau2 were internally swapped (visible on prepared input printout)");
426  } else {
427  Info("DiTauMassTools", "tau1 and tau2 were NOT internally swapped");
428  }
429 
430  Info("DiTauMassTools", "%s",
431  (" MEtLMin=" + std::to_string(m_MEtLMin) + " MEtLMax=" + std::to_string(m_MEtLMax)).c_str());
432  Info("DiTauMassTools", "%s",
433  (" MEtPMin=" + std::to_string(m_MEtPMin) + " MEtPMax=" + std::to_string(m_MEtPMax)).c_str());
434  Info("DiTauMassTools", "%s",
435  (" Phi1Min=" + std::to_string(m_Phi1Min) + " Phi1Max=" + std::to_string(m_Phi1Max)).c_str());
436  Info("DiTauMassTools", "%s",
437  (" Phi2Min=" + std::to_string(m_Phi2Min) + " Phi2Max=" + std::to_string(m_Phi2Max)).c_str());
438  Info("DiTauMassTools", "%s",
439  (" Mnu1Min=" + std::to_string(m_Mnu1Min) + " Mnu1Max=" + std::to_string(m_Mnu1Max)).c_str());
440  Info("DiTauMassTools", "%s",
441  (" Mnu2Min=" + std::to_string(m_Mnu2Min) + " Mnu2Max=" + std::to_string(m_Mnu2Max)).c_str());
442 }
443 
444 // Printout of final results
446 
447  if (preparedInput.m_fUseVerbose != 1)
448  return;
449 
450  const PtEtaPhiMVector *origVisTau1 = 0;
451  const PtEtaPhiMVector *origVisTau2 = 0;
452 
453  if (preparedInput.m_InputReorder == 0) {
454  origVisTau1 = &preparedInput.m_vistau1;
455  origVisTau2 = &preparedInput.m_vistau2;
456  } else // input order was flipped
457  {
458  origVisTau1 = &preparedInput.m_vistau2;
459  origVisTau2 = &preparedInput.m_vistau1;
460  }
461 
462  PrintOtherInput();
463 
464  Info("DiTauMassTools",
465  "------------- Printing Final Results for MissingMassCalculator --------------");
466  Info("DiTauMassTools",
467  ".............................................................................");
468  Info("DiTauMassTools", "%s", ("Fit status=" + std::to_string(OutputInfo.m_FitStatus)).c_str());
469 
470  for (int imeth = 0; imeth < MMCFitMethod::MAX; ++imeth) {
471  Info("DiTauMassTools", "%s",
472  ("___ Results for " + MMCFitMethod::name[imeth] + "Method ___")
473  .c_str());
474  Info("DiTauMassTools", "%s",
475  (" signif=" + std::to_string(OutputInfo.m_FitSignificance[imeth])).c_str());
476  Info("DiTauMassTools", "%s", (" mass=" + std::to_string(OutputInfo.m_FittedMass[imeth])).c_str());
477  Info("DiTauMassTools", "%s", (" rms/mpv=" + std::to_string(OutputInfo.m_RMS2MPV)).c_str());
478 
479  if (imeth == MMCFitMethod::MLM) {
480  Info("DiTauMassTools", " no 4-momentum or MET from this method ");
481  continue;
482  }
483 
484  if (OutputInfo.m_FitStatus <= 0) {
485  Info("DiTauMassTools", " fit failed ");
486  }
487 
488  const PtEtaPhiMVector &tlvnu1 = OutputInfo.m_nuvec1[imeth];
489  const PtEtaPhiMVector &tlvnu2 = OutputInfo.m_nuvec2[imeth];
490  const PtEtaPhiMVector &tlvo1 = OutputInfo.m_objvec1[imeth];
491  const PtEtaPhiMVector &tlvo2 = OutputInfo.m_objvec2[imeth];
492  const XYVector &tvmet = OutputInfo.m_FittedMetVec[imeth];
493 
494  Info("DiTauMassTools", "%s",
495  (" Neutrino-1: P=" + std::to_string(tlvnu1.P()) + " Pt=" + std::to_string(tlvnu1.Pt()) +
496  " Eta=" + std::to_string(tlvnu1.Eta()) + " Phi=" + std::to_string(tlvnu1.Phi()) +
497  " M=" + std::to_string(tlvnu1.M()) + " Px=" + std::to_string(tlvnu1.Px()) +
498  " Py=" + std::to_string(tlvnu1.Py()) + " Pz=" + std::to_string(tlvnu1.Pz()))
499  .c_str());
500  Info("DiTauMassTools", "%s",
501  (" Neutrino-2: P=" + std::to_string(tlvnu2.P()) + " Pt=" + std::to_string(tlvnu2.Pt()) +
502  " Eta=" + std::to_string(tlvnu2.Eta()) + " Phi=" + std::to_string(tlvnu2.Phi()) +
503  " M=" + std::to_string(tlvnu2.M()) + " Px=" + std::to_string(tlvnu2.Px()) +
504  " Py=" + std::to_string(tlvnu2.Py()) + " Pz=" + std::to_string(tlvnu2.Pz()))
505  .c_str());
506  Info("DiTauMassTools", "%s",
507  (" Tau-1: P=" + std::to_string(tlvo1.P()) + " Pt=" + std::to_string(tlvo1.Pt()) +
508  " Eta=" + std::to_string(tlvo1.Eta()) + " Phi=" + std::to_string(tlvo1.Phi()) +
509  " M=" + std::to_string(tlvo1.M()) + " Px=" + std::to_string(tlvo1.Px()) +
510  " Py=" + std::to_string(tlvo1.Py()) + " Pz=" + std::to_string(tlvo1.Pz()))
511  .c_str());
512  Info("DiTauMassTools", "%s",
513  (" Tau-2: P=" + std::to_string(tlvo2.P()) + " Pt=" + std::to_string(tlvo2.Pt()) +
514  " Eta=" + std::to_string(tlvo2.Eta()) + " Phi=" + std::to_string(tlvo2.Phi()) +
515  " M=" + std::to_string(tlvo2.M()) + " Px=" + std::to_string(tlvo2.Px()) +
516  " Py=" + std::to_string(tlvo2.Py()) + " Pz=" + std::to_string(tlvo2.Pz()))
517  .c_str());
518 
519  Info("DiTauMassTools", "%s",
520  (" dR(nu1-visTau1)=" + std::to_string(DeltaR(tlvnu1,*origVisTau1))).c_str());
521  Info("DiTauMassTools", "%s",
522  (" dR(nu2-visTau2)=" + std::to_string(DeltaR(tlvnu2,*origVisTau2))).c_str());
523 
524  Info("DiTauMassTools", "%s",
525  (" Fitted MET =" + std::to_string(tvmet.R()) + " Phi=" + std::to_string(tlvnu1.Phi()) +
526  " Px=" + std::to_string(tvmet.X()) + " Py=" + std::to_string(tvmet.Y()))
527  .c_str());
528 
529  Info("DiTauMassTools", "%s", (" Resonance: P=" + std::to_string(OutputInfo.m_totalvec[imeth].P()) +
530  " Pt=" + std::to_string(OutputInfo.m_totalvec[imeth].Pt()) +
531  " Eta=" + std::to_string(OutputInfo.m_totalvec[imeth].Eta()) +
532  " Phi=" + std::to_string(OutputInfo.m_totalvec[imeth].Phi()) +
533  " M=" + std::to_string(OutputInfo.m_totalvec[imeth].M()) +
534  " Px=" + std::to_string(OutputInfo.m_totalvec[imeth].Px()) +
535  " Py=" + std::to_string(OutputInfo.m_totalvec[imeth].Py()) +
536  " Pz=" + std::to_string(OutputInfo.m_totalvec[imeth].Pz()))
537  .c_str());
538  }
539 
540  return;
541 }
542 
543 // returns P1, P2, and theta1 & theta2 solutions
544 // This compute the nu1 nu2 solution in the most efficient way. Wrt to
545 // NuPsolutionV2, the output nu1 nu2 4-vector have non zero mass (if relevant).
546 // It is not optimised for grid running so much less caching is done (which
547 // makes it more readable). Only quantities fixed within an event are cached.
548 // relies on a number of these variables to be initialised before the loop.
549 
550 int MissingMassCalculator::NuPsolutionV3(const double &mNu1, const double &mNu2,
551  const double &phi1, const double &phi2,
552  int &nsol1, int &nsol2) {
553 
554  // Pv1, Pv2 : visible tau decay product momentum
555  // Pn1 Pn2 : neutrino momentum
556  // phi1, phi2 : neutrino azymutal angles
557  // PTmiss2=PTmissy Cos[phi2] - PTmissx Sin[phi2]
558  // PTmiss2cscdphi=PTmiss2/Sin[phi1-phi2]
559  // Pv1proj=Pv1x Cos[phi1] + Pv1y Sin[phi1]
560  // M2noma1=Mtau^2-Mv1^2-Mn1^2
561  // ETv1^2=Ev1^2-Pv1z^2
562 
563  // discriminant : 16 Ev1^2 (M2noma1^2 + 4 M2noma1 PTmiss2cscdphi Pv1proj - 4
564  // (ETv1^2 (Mn1^2 + PTmiss2cscdphi^2) - PTmiss2cscdphi^2 Pv1proj^2))
565  // two solutions for epsilon = +/-1
566  // Pn1z=(1/(2 ETv1^2))(epsilon Ev1 Sqrt[ M2noma1^2 + 4 M2noma1 PTmiss2cscdphi
567  // Pv1proj - 4 (ETv1^2 (Mn1^2 + qPTmiss2cscdphi^2) - PTmiss2cscdphi^2
568  // Pv1proj^2)] + M2noma1 Pv1z + 2 PTmiss2cscdphi Pv1proj Pv1z)
569  // with conditions: M2noma1 + 2 PTmiss2cscdphi Pv1proj + 2 Pn1z Pv1z > 0
570  // PTn1 -> PTmiss2 Csc[phi1 - phi2]
571 
572  // if initialisation precompute some quantities
573  int solution_code = 0; // 0 with no solution, 1 with solution
574  nsol1 = 0;
575  nsol2 = 0;
576 
577  // Variables used to test PTn1 and PTn2 > 0
578 
579  const double &pTmissx = preparedInput.m_MEtX;
580  const double &pTmissy = preparedInput.m_MEtY;
581 
583  double pTmiss2 = pTmissy * m_cosPhi2 - pTmissx * m_sinPhi2;
584 
585  int dPhiSign = 0;
586  dPhiSign = fixPhiRange(phi1 - phi2) > 0 ? +1 : -1;
587 
588  // Test if PTn1 and PTn2 > 0. Then MET vector is between the two neutrino
589  // vector
590 
591  if (pTmiss2 * dPhiSign < 0) {
592  ++m_testptn1;
593  return solution_code;
594  }
595 
597  double pTmiss1 = pTmissy * m_cosPhi1 - pTmissx * m_sinPhi1;
598 
599  if (pTmiss1 * (-dPhiSign) < 0) {
600  ++m_testptn2;
601  return solution_code;
602  }
603 
604  // Variables used to calculate discri1
605 
606  double m2Vis1 = m_tauVec1M * m_tauVec1M;
608  m_m2Nu1 = mNu1 * mNu1;
609  double m2noma1 = m_mTau2 - m_m2Nu1 - m2Vis1;
610  double m4noma1 = m2noma1 * m2noma1;
611  double pv1proj = m_tauVec1Px * m_cosPhi1 + m_tauVec1Py * m_sinPhi1;
612  double p2v1proj = std::pow(pv1proj, 2);
613  double sinDPhi2 = m_cosPhi2 * m_sinPhi1 - m_sinPhi2 * m_cosPhi1; // sin(Phi1-Phi2)
614  double pTmiss2CscDPhi = pTmiss2 / sinDPhi2;
615  double &pTn1 = pTmiss2CscDPhi;
616  double pT2miss2CscDPhi = pTmiss2CscDPhi * pTmiss2CscDPhi;
617 
618  // Test on discri1
619  const double discri1 = m4noma1 + 4 * m2noma1 * pTmiss2CscDPhi * pv1proj -
620  4 * (m_ET2v1 * (m_m2Nu1 + pT2miss2CscDPhi) - (pT2miss2CscDPhi * p2v1proj));
621 
622  if (discri1 < 0) // discriminant negative -> no solution
623  {
624  ++m_testdiscri1;
625  return solution_code;
626  }
627 
628  // Variables used to calculate discri2
629  double m2Vis2 = m_tauVec2M * m_tauVec2M;
631  m_m2Nu2 = mNu2 * mNu2;
632  double m2noma2 = m_mTau2 - m_m2Nu2 - m2Vis2;
633  double m4noma2 = m2noma2 * m2noma2;
634  double pv2proj = m_tauVec2Px * m_cosPhi2 + m_tauVec2Py * m_sinPhi2;
635  double p2v2proj = std::pow(pv2proj, 2);
636  double sinDPhi1 = -sinDPhi2;
637  double pTmiss1CscDPhi = pTmiss1 / sinDPhi1;
638  double &pTn2 = pTmiss1CscDPhi;
639  double pT2miss1CscDPhi = pTmiss1CscDPhi * pTmiss1CscDPhi;
640 
641  const double discri2 = m4noma2 + 4 * m2noma2 * pTmiss1CscDPhi * pv2proj -
642  4 * (m_ET2v2 * (m_m2Nu2 + pT2miss1CscDPhi) - (pT2miss1CscDPhi * p2v2proj));
643 
644  if (discri2 < 0) // discriminant negative -> no solution
645  {
646  ++m_testdiscri2;
647  return solution_code;
648  }
649 
650  // this should be done only once we know there are solutions for nu2
652  m_Ev1 = sqrt(m_E2v1);
653  double sqdiscri1 = sqrt(discri1);
654  double first1 =
655  (m2noma1 * m_tauVec1Pz + 2 * pTmiss2CscDPhi * pv1proj * m_tauVec1Pz) / (2 * m_ET2v1);
656  double second1 = sqdiscri1 * m_Ev1 / (2 * m_ET2v1);
657 
658  // first solution
659  double pn1Z = first1 + second1;
660 
661  if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z * m_tauVec1Pz >
662  0) // Condition for solution to exist
663  {
664  m_nuvecsol1[nsol1].SetPxPyPzE(pTn1 * m_cosPhi1, pTn1 * m_sinPhi1, pn1Z,
665  sqrt(std::pow(pTn1, 2) + std::pow(pn1Z, 2) + m_m2Nu1));
666 
667  ++nsol1;
668  }
669 
670  pn1Z = first1 - second1;
671 
672  if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z * m_tauVec1Pz >
673  0) // Condition for solution to exist
674  {
675 
676  m_nuvecsol1[nsol1].SetPxPyPzE(pTn1 * m_cosPhi1, pTn1 * m_sinPhi1, pn1Z,
677  sqrt(std::pow(pTn1, 2) + std::pow(pn1Z, 2) + m_m2Nu1));
678 
679  ++nsol1;
680  }
681 
682  if (nsol1 == 0) {
683  ++m_nosol1;
684  return solution_code;
685  }
686 
688  m_Ev2 = sqrt(m_E2v2);
689  double sqdiscri2 = sqrt(discri2);
690  double first2 =
691  (m2noma2 * m_tauVec2Pz + 2 * pTmiss1CscDPhi * pv2proj * m_tauVec2Pz) / (2 * m_ET2v2);
692  double second2 = sqdiscri2 * m_Ev2 / (2 * m_ET2v2);
693 
694  // second solution
695  double pn2Z = first2 + second2;
696 
697  if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z * m_tauVec2Pz >
698  0) // Condition for solution to exist
699  {
700  m_nuvecsol2[nsol2].SetPxPyPzE(pTn2 * m_cosPhi2, pTn2 * m_sinPhi2, pn2Z,
701  sqrt(std::pow(pTn2, 2) + std::pow(pn2Z, 2) + m_m2Nu2));
702 
703  ++nsol2;
704  }
705 
706  pn2Z = first2 - second2;
707  ;
708 
709  if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z * m_tauVec2Pz >
710  0) // Condition for solution to exist
711  {
712  m_nuvecsol2[nsol2].SetPxPyPzE(pTn2 * m_cosPhi2, pTn2 * m_sinPhi2, pn2Z,
713  sqrt(std::pow(pTn2, 2) + std::pow(pn2Z, 2) + m_m2Nu2));
714 
715  ++nsol2;
716  }
717 
718  if (nsol2 == 0) {
719  ++m_nosol2;
720  return solution_code;
721  }
722 
723  // Verification if solution exist
724 
725  solution_code = 1;
726  ++m_iterNuPV3;
727 
728  // double check solutions from time to time
729  if (m_iterNuPV3 % 1000 == 1) {
730  double pnux = m_nuvecsol1[0].Px() + m_nuvecsol2[0].Px();
731  double pnuy = m_nuvecsol1[0].Py() + m_nuvecsol2[0].Py();
732  double mtau1plus = (m_nuvecsol1[0] + m_tauVec1).M();
733  double mtau1moins = (m_nuvecsol1[1] + m_tauVec1).M();
734  double mtau2plus = (m_nuvecsol2[0] + m_tauVec2).M();
735  double mtau2moins = (m_nuvecsol2[1] + m_tauVec2).M();
736  if (std::abs(pnux - pTmissx) > 0.001 || std::abs(pnuy - pTmissy) > 0.001) {
737  Info("DiTauMassTools", "%s", ("NuPsolutionV3 ERROR Pnux-Met.X or Pnuy-Met.Y > 0.001 : " +
738  std::to_string(pnux - pTmissx) + " and " +
739  std::to_string(pnuy - pTmissx) + " " + "Invalid solutions")
740  .c_str());
741  }
742  if (std::abs(mtau1plus - m_mTau) > 0.001 || std::abs(mtau1moins - m_mTau) > 0.001 ||
743  std::abs(mtau2plus - m_mTau) > 0.001 || std::abs(mtau2moins - m_mTau) > 0.001) {
744  Info("DiTauMassTools", "%s", ("NuPsolutionV3 ERROR tau mass not recovered : " +
745  std::to_string(mtau1plus) + " " + std::to_string(mtau1moins) + " " +
746  std::to_string(mtau2plus) + " " + std::to_string(mtau2moins))
747  .c_str());
748  }
749  }
750 
751  return solution_code;
752 }
753 
754 // returns solution for Lepton Flavor Violating X->lep+tau studies
755 int MissingMassCalculator::NuPsolutionLFV(const XYVector &met_vec,
756  const PtEtaPhiMVector &tau, const double &l_nu,
757  std::vector<PtEtaPhiMVector> &nu_vec) {
758  int solution_code = 0; // 0 with no solution, 1 with solution
759 
760  nu_vec.clear();
761  PxPyPzMVector nu(met_vec.X(), met_vec.Y(), 0.0, l_nu);
762  PxPyPzMVector nu2(met_vec.X(), met_vec.Y(), 0.0, l_nu);
763 
764  const double Mtau = ParticleConstants::tauMassInMeV / GEV;
765  // double msq = (Mtau*Mtau-tau.M()*tau.M())/2;
766  double msq = (Mtau * Mtau - tau.M() * tau.M() - l_nu * l_nu) /
767  2; // to take into account the fact that 2-nu systema has mass
768  double gamma = nu.Px() * nu.Px() + nu.Py() * nu.Py();
769  double beta = tau.Px() * nu.Px() + tau.Py() * nu.Py() + msq;
770  double a = tau.E() * tau.E() - tau.Pz() * tau.Pz();
771  double b = -2 * tau.Pz() * beta;
772  double c = tau.E() * tau.E() * gamma - beta * beta;
773  if ((b * b - 4 * a * c) < 0)
774  return solution_code; // no solution found
775  else
776  solution_code = 2;
777  double pvz1 = (-b + sqrt(b * b - 4 * a * c)) / (2 * a);
778  double pvz2 = (-b - sqrt(b * b - 4 * a * c)) / (2 * a);
779 
780  nu.SetCoordinates(met_vec.X(), met_vec.Y(), pvz1, l_nu);
781  nu2.SetCoordinates(met_vec.X(), met_vec.Y(), pvz2, l_nu);
782 
783  PtEtaPhiMVector return_nu(nu.Pt(), nu.Eta(), nu.Phi(), nu.M());
784  PtEtaPhiMVector return_nu2(nu2.Pt(), nu2.Eta(), nu2.Phi(), nu2.M());
785  nu_vec.push_back(return_nu);
786  nu_vec.push_back(return_nu2);
787  return solution_code;
788 }
789 
790 // like v9fast, but the parameter space scanning is now factorised out, to allow
791 // flexibility
793 
794  int nsuccesses = 0;
795 
796  int fit_code = 0; // 0==bad, 1==good
799  OutputInfo.m_AveSolRMS = 0.;
800 
801  m_fMfit_all->Reset();
802 
803  if(m_SaveLlhHisto){
804  m_fMEtP_all->Reset();
805  m_fMEtL_all->Reset();
806  m_fMnu1_all->Reset();
807  m_fMnu2_all->Reset();
808  m_fPhi1_all->Reset();
809  m_fPhi2_all->Reset();
810  }
811 
812  m_fMfit_allNoWeight->Reset();
813  m_fPXfit1->Reset();
814  m_fPYfit1->Reset();
815  m_fPZfit1->Reset();
816  m_fPXfit2->Reset();
817  m_fPYfit2->Reset();
818  m_fPZfit2->Reset();
819 
820  // these histograms are used for the floating stopping criterion
821  if (m_fUseFloatStopping) {
822  m_fMmass_split1->Reset();
823  m_fMEtP_split1->Reset();
824  m_fMEtL_split1->Reset();
825  m_fMnu1_split1->Reset();
826  m_fMnu2_split1->Reset();
827  m_fPhi1_split1->Reset();
828  m_fPhi2_split1->Reset();
829  m_fMmass_split2->Reset();
830  m_fMEtP_split2->Reset();
831  m_fMEtL_split2->Reset();
832  m_fMnu1_split2->Reset();
833  m_fMnu2_split2->Reset();
834  m_fPhi1_split2->Reset();
835  m_fPhi2_split2->Reset();
836  }
837 
838  m_prob_tmp = 0.0;
839 
840  m_iter1 = 0;
841 
842  m_totalProbSum = 0;
843  m_mtautauSum = 0;
844 
845  XYVector deltamet_vec;
846 
847  // initialize a spacewalker, which walks the parameter space according to some
848  // algorithm
849  SpaceWalkerInit();
850 
851  while (SpaceWalkerWalk()) {
852  bool paramInsideRange = false;
853  m_nsol = 0;
854 
855  paramInsideRange = checkAllParamInRange();
856 
857  // FIXME if no tau scanning, or symmetric matrices, rotatin is made twice
858  // which is inefficient
859  const double deltaMetx = m_MEtL * m_metCovPhiCos - m_MEtP * m_metCovPhiSin;
860  const double deltaMety = m_MEtL * m_metCovPhiSin + m_MEtP * m_metCovPhiCos;
861 
862  // deltaMetVec.Set(met_smear_x,met_smear_y);
864  preparedInput.m_inputMEtY + deltaMety);
865 
866  // save in global variable for speed sake
870 
871  if (paramInsideRange)
873 
874  // DR for markov chain need to enter handleSolution also when zero solutions
875  handleSolutions();
876  // be careful that with markov, current solution is from now on stored in
877  // XYZOldSolVec
878 
879  if (m_nsol <= 0)
880  continue;
881 
882  // for markov, nsuccess more difficult to define. Decide this is the number
883  // of independent point accepted (hence without weight)
884  nsuccesses = m_markovNAccept;
886 
887  m_iter1 += m_nsol;
888  fit_code = 1;
889 
890  } // while loop
891 
893  OutputInfo.m_NSuccesses = nsuccesses;
894 
895  if (nsuccesses > 0) {
896  OutputInfo.m_AveSolRMS /= nsuccesses;
897  } else {
898  OutputInfo.m_AveSolRMS = -1.;
899  }
900 
901  double Px1, Py1, Pz1;
902  double Px2, Py2, Pz2;
903  if (nsuccesses > 0) {
904 
905  // note that smoothing can slightly change the integral of the histogram
906 
907 #ifdef SMOOTH
908  m_fMfit_all->Smooth();
909  m_fMfit_allNoWeight->Smooth();
910  m_fPXfit1->Smooth();
911  m_fPYfit1->Smooth();
912  m_fPZfit1->Smooth();
913  m_fPXfit2->Smooth();
914  m_fPYfit2->Smooth();
915  m_fPZfit2->Smooth();
916 #endif
917 
918  // default max finding method defined in MissingMassCalculator.h
919  // note that window defined in terms of number of bin, so depend on binning
920  std::vector<double> histInfo(HistInfo::MAXHISTINFO);
922  double prob_hist = histInfo.at(HistInfo::PROB);
923 
924  if (prob_hist != 0.0)
925  m_fDitauStuffHisto.Sign_best = -log10(std::abs(prob_hist));
926  else {
927  // this mean the histogram is empty.
928  // possible but very rare if all entries outside histogram range
929  // fall back to maximum
932  }
933 
936  std::vector<double> histInfoOther(HistInfo::MAXHISTINFO);
937  //---- getting full tau1 momentum
938  Px1 = maxFromHist(m_fPXfit1, histInfoOther);
939  Py1 = maxFromHist(m_fPYfit1, histInfoOther);
940  Pz1 = maxFromHist(m_fPZfit1, histInfoOther);
941 
942  //---- getting full tau2 momentum
943  Px2 = maxFromHist(m_fPXfit2, histInfoOther);
944  Py2 = maxFromHist(m_fPYfit2, histInfoOther);
945  Pz2 = maxFromHist(m_fPZfit2, histInfoOther);
946 
947  //---- setting 4-vecs
948  PxPyPzMVector fulltau1, fulltau2;
949  fulltau1.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
950  fulltau2.SetCoordinates(Px2, Py2, Pz2, ParticleConstants::tauMassInMeV / GEV);
951  // PtEtaPhiMVector fulltau1(_fulltau1.Pt(), _fulltau1.Eta(), _fulltau1.Phi(), _fulltau1.M());
952  //PtEtaPhiMVector fulltau2(_fulltau2.Pt(), _fulltau2.Eta(), _fulltau2.Phi(), _fulltau2.M());
953 
954  if (fulltau1.P() < preparedInput.m_vistau1.P())
955  fulltau1 = 1.01 * preparedInput.m_vistau1; // protection against cases when fitted tau
956  // momentum is smaller than visible tau momentum
957  if (fulltau2.P() < preparedInput.m_vistau2.P())
958  fulltau2 = 1.01 * preparedInput.m_vistau2; // protection against cases when fitted tau
959  // momentum is smaller than visible tau momentum
960  m_fDitauStuffHisto.vistau1 = preparedInput.m_vistau1; // FIXME should also be fitted if tau scan
962  m_fDitauStuffHisto.nutau1 = fulltau1 - preparedInput.m_vistau1; // these are the original tau vis
964  fulltau2 - preparedInput.m_vistau2; // FIXME neutrino mass not necessarily zero
965  }
966 
967  // Note that for v9walk, points outside the METx MEty disk are counted, while
968  // this was not the case for v9
969  if (preparedInput.m_fUseVerbose == 1) {
970  Info("DiTauMassTools", "Scanning ");
971  Info("DiTauMassTools", " Markov ");
972  Info("DiTauMassTools", "%s",
973  (" V9W niters=" + std::to_string(m_iter0) + " " + std::to_string(m_iter1)).c_str());
974  Info("DiTauMassTools", "%s", (" nFullScan " + std::to_string(m_markovNFullScan)).c_str());
975  Info("DiTauMassTools", "%s", (" nRejectNoSol " + std::to_string(m_markovNRejectNoSol)).c_str());
976  Info("DiTauMassTools", "%s", (" nRejectMetro " + std::to_string(m_markovNRejectMetropolis)).c_str());
977  Info("DiTauMassTools", "%s", (" nAccept " + std::to_string(m_markovNAccept)).c_str());
978  Info("DiTauMassTools", "%s",
979  (" probsum " + std::to_string(m_totalProbSum) + " msum " + std::to_string(m_mtautauSum))
980  .c_str());
981  }
982 
983  if (preparedInput.m_fUseVerbose == 1) {
984  if (fit_code == 0) {
985  Info("DiTauMassTools", "%s", ("!!!----> Warning-3 in "
986  "MissingMassCalculator::DitauMassCalculatorV9Walk() : fit status=" +
987  std::to_string(fit_code))
988  .c_str());
989  Info("DiTauMassTools", "%s", "....... No solution is found. Printing input info .......");
990 
991  Info("DiTauMassTools", "%s", (" vis Tau-1: Pt=" + std::to_string(preparedInput.m_vistau1.Pt()) +
992  " M=" + std::to_string(preparedInput.m_vistau1.M()) +
993  " eta=" + std::to_string(preparedInput.m_vistau1.Eta()) +
994  " phi=" + std::to_string(preparedInput.m_vistau1.Phi()) +
996  .c_str());
997  Info("DiTauMassTools", "%s", (" vis Tau-2: Pt=" + std::to_string(preparedInput.m_vistau2.Pt()) +
998  " M=" + std::to_string(preparedInput.m_vistau2.M()) +
999  " eta=" + std::to_string(preparedInput.m_vistau2.Eta()) +
1000  " phi=" + std::to_string(preparedInput.m_vistau2.Phi()) +
1002  .c_str());
1003  Info("DiTauMassTools", "%s", (" MET=" + std::to_string(preparedInput.m_MetVec.R()) +
1004  " Met_X=" + std::to_string(preparedInput.m_MetVec.X()) +
1005  " Met_Y=" + std::to_string(preparedInput.m_MetVec.Y()))
1006  .c_str());
1007  Info("DiTauMassTools", " ---------------------------------------------------------- ");
1008  }
1009  }
1010 
1011  return fit_code;
1012 }
1013 
1015 
1016  // debugThisIteration=false;
1017  m_debugThisIteration = true;
1018 
1019  int fit_code = 0; // 0==bad, 1==good
1022  OutputInfo.m_NTrials = 0;
1024  OutputInfo.m_AveSolRMS = 0.;
1025 
1026  //------- Settings -------------------------------
1027  int NiterMET = m_niter_fit2; // number of iterations for each MET scan loop
1028  int NiterMnu = m_niter_fit3; // number of iterations for Mnu loop
1029  const double Mtau = ParticleConstants::tauMassInMeV / GEV;
1030  double Mnu_binSize = m_MnuScanRange / NiterMnu;
1031 
1032  double METresX = preparedInput.m_METsigmaL; // MET resolution in direction parallel to
1033  // leading jet, for MET scan
1034  double METresY = preparedInput.m_METsigmaP; // MET resolution in direction perpendicular to
1035  // leading jet, for MET scan
1036 
1037  //-------- end of Settings
1038 
1039  // if m_nsigma_METscan was not set by user, set to default values
1040  if(m_nsigma_METscan == -1){
1041  if (preparedInput.m_tauTypes == TauTypes::ll) { // both tau's are leptonic
1043  } else if (preparedInput.m_tauTypes == TauTypes::lh) { // lep had
1045  }
1046  }
1047 
1048  double N_METsigma = m_nsigma_METscan; // number of sigmas for MET scan
1049  double METresX_binSize = 2 * N_METsigma * METresX / NiterMET;
1050  double METresY_binSize = 2 * N_METsigma * METresY / NiterMET;
1051 
1052  int solution = 0;
1053 
1054  std::vector<PtEtaPhiMVector> nu_vec;
1055 
1056  m_totalProbSum = 0;
1057  m_mtautauSum = 0;
1058 
1059  double metprob = 1.0;
1060  double sign_tmp = 0.0;
1061  double tauprob = 1.0;
1062  double totalProb = 0.0;
1063 
1064  m_prob_tmp = 0.0;
1065 
1066  double met_smear_x = 0.0;
1067  double met_smear_y = 0.0;
1068  double met_smearL = 0.0;
1069  double met_smearP = 0.0;
1070 
1071  double angle1 = 0.0;
1072 
1073  if (m_fMfit_all) {
1074  m_fMfit_all->Reset();
1075  }
1076  if (m_fMfit_allNoWeight) {
1077  m_fMfit_allNoWeight->Reset();
1078  }
1079  if (m_fPXfit1) {
1080  m_fPXfit1->Reset();
1081  }
1082  if (m_fPYfit1) {
1083  m_fPYfit1->Reset();
1084  }
1085  if (m_fPZfit1) {
1086  m_fPZfit1->Reset();
1087  }
1088 
1089  int iter0 = 0;
1090  m_iter1 = 0;
1091  m_iter2 = 0;
1092  m_iter3 = 0;
1093  m_iter4 = 0;
1094 
1095  const double met_coscovphi = cos(preparedInput.m_METcovphi);
1096  const double met_sincovphi = sin(preparedInput.m_METcovphi);
1097 
1098  m_iang1low = 0;
1099  m_iang1high = 0;
1100 
1101  // double Mvis=(tau_vec1+tau_vec2).M();
1102  // PtEtaPhiMVector met4vec(0.0,0.0,0.0,0.0);
1103  // met4vec.SetPxPyPzE(met_vec.X(),met_vec.Y(),0.0,met_vec.R());
1104  // double Meff=(tau_vec1+tau_vec2+met4vec).M();
1105  // double met_det=met_vec.R();
1106 
1107  //---------------------------------------------
1108  if (preparedInput.m_tauTypes == TauTypes::ll) // dilepton case
1109  {
1110  if (preparedInput.m_fUseVerbose == 1) {
1111  Info("DiTauMassTools", "Running in dilepton mode");
1112  }
1113  double input_metX = preparedInput.m_MetVec.X();
1114  double input_metY = preparedInput.m_MetVec.Y();
1115 
1116  PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1117  PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1118  int tau_type_tmp;
1119  int tau_ind = 0;
1120 
1121  if (preparedInput.m_LFVmode == 1) // muon case: H->mu+tau(->ele) decays
1122  {
1123  if ((preparedInput.m_vistau1.M() > 0.05 &&
1124  preparedInput.m_vistau2.M() < 0.05) != refit) // choosing lepton from Higgs decay
1125  //When the mass calculator is rerun with refit==true the alternative lepton ordering is used
1126  {
1127  tau_tmp = preparedInput.m_vistau2;
1128  lep_tmp = preparedInput.m_vistau1;
1129  tau_type_tmp = preparedInput.m_type_visTau2;
1130  tau_ind = 2;
1131  } else {
1132  tau_tmp = preparedInput.m_vistau1;
1133  lep_tmp = preparedInput.m_vistau2;
1134  tau_type_tmp = preparedInput.m_type_visTau1;
1135  tau_ind = 1;
1136  }
1137  }
1138  if (preparedInput.m_LFVmode == 0) // electron case: H->ele+tau(->mu) decays
1139  {
1140  if ((preparedInput.m_vistau1.M() < 0.05 &&
1141  preparedInput.m_vistau2.M() > 0.05) != refit) // choosing lepton from Higgs decay
1142  //When the mass calculator is rerun with refit=true the alternative lepton ordering is used
1143  {
1144  tau_tmp = preparedInput.m_vistau2;
1145  lep_tmp = preparedInput.m_vistau1;
1146  tau_type_tmp = preparedInput.m_type_visTau2;
1147  tau_ind = 2;
1148  } else {
1149  tau_tmp = preparedInput.m_vistau1;
1150  lep_tmp = preparedInput.m_vistau2;
1151  tau_type_tmp = preparedInput.m_type_visTau1;
1152  tau_ind = 1;
1153  }
1154  }
1155 
1156  //------- Settings -------------------------------
1157  double Mlep = tau_tmp.M();
1158  // double dMnu_max=m_MnuScanRange-Mlep;
1159  // double Mnu_binSize=dMnu_max/NiterMnu;
1160  //-------- end of Settings
1161 
1162  // double M=Mtau;
1163  double M_nu = 0.0;
1164  double MnuProb = 1.0;
1165  //---------------------------------------------
1166  for (int i3 = 0; i3 < NiterMnu; i3++) //---- loop-3: virtual neutrino mass
1167  {
1168  M_nu = Mnu_binSize * i3;
1169  if (M_nu >= (Mtau - Mlep))
1170  continue;
1171  // M=sqrt(Mtau*Mtau-M_nu*M_nu);
1172  MnuProb = Prob->MnuProbability(preparedInput, M_nu,
1173  Mnu_binSize); // Mnu probability
1174  //---------------------------------------------
1175  for (int i4 = 0; i4 < NiterMET + 1; i4++) // MET_X scan
1176  {
1177  met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1178  for (int i5 = 0; i5 < NiterMET + 1; i5++) // MET_Y scan
1179  {
1180  met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1181  if (pow(met_smearL / METresX, 2) + pow(met_smearP / METresY, 2) > pow(N_METsigma, 2))
1182  continue; // use ellipse instead of square
1183  met_smear_x = met_smearL * met_coscovphi - met_smearP * met_sincovphi;
1184  met_smear_y = met_smearL * met_sincovphi + met_smearP * met_coscovphi;
1185  metvec_tmp.SetXY(input_metX + met_smear_x, input_metY + met_smear_y);
1186 
1187  solution = NuPsolutionLFV(metvec_tmp, tau_tmp, M_nu, nu_vec);
1188 
1189  ++iter0;
1190 
1191  if (solution < 1)
1192  continue;
1193  ++m_iter1;
1194 
1195  // if fast sin cos, result to not match exactly nupsolutionv2, so skip
1196  // test
1197  // SpeedUp no nested loop to compute individual probability
1198  int ngoodsol1 = 0;
1199 
1200  metprob = Prob->MetProbability(preparedInput, met_smearL, met_smearP, METresX, METresY);
1201  if (metprob <= 0)
1202  continue;
1203  for (unsigned int j1 = 0; j1 < nu_vec.size(); j1++) {
1204  if (tau_tmp.E() + nu_vec[j1].E() >= preparedInput.m_beamEnergy)
1205  continue;
1206  const double tau1_tmpp = (tau_tmp + nu_vec[j1]).P();
1207  angle1 = Angle(nu_vec[j1], tau_tmp);
1208 
1209  if (angle1 < dTheta3DLimit(tau_type_tmp, 0, tau1_tmpp)) {
1210  ++m_iang1low;
1211  continue;
1212  } // lower 99% bound
1213  if (angle1 > dTheta3DLimit(tau_type_tmp, 1, tau1_tmpp)) {
1214  ++m_iang1high;
1215  continue;
1216  } // upper 99% bound
1217  double tauvecprob1j =
1218  Prob->dTheta3d_probabilityFast(preparedInput, tau_type_tmp, angle1, tau1_tmpp);
1219  if (tauvecprob1j == 0.)
1220  continue;
1221  tauprob = Prob->TauProbabilityLFV(preparedInput, tau_type_tmp, tau_tmp, nu_vec[j1]);
1222  totalProb = tauvecprob1j * metprob * MnuProb * tauprob;
1223 
1224  m_tautau_tmp.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
1225  m_tautau_tmp += tau_tmp;
1226  m_tautau_tmp += lep_tmp;
1227  m_tautau_tmp += nu_vec[j1];
1228 
1229  const double mtautau = m_tautau_tmp.M();
1230 
1231  m_totalProbSum += totalProb;
1232  m_mtautauSum += mtautau;
1233 
1234  fit_code = 1; // at least one solution is found
1235 
1236  m_fMfit_all->Fill(mtautau, totalProb);
1237  m_fMfit_allNoWeight->Fill(mtautau, 1.);
1238  //----------------- using P*fit to fill Px,y,z_tau
1239  m_fPXfit1->Fill((tau_tmp + nu_vec[j1]).Px(), totalProb);
1240  m_fPYfit1->Fill((tau_tmp + nu_vec[j1]).Py(), totalProb);
1241  m_fPZfit1->Fill((tau_tmp + nu_vec[j1]).Pz(), totalProb);
1242 
1243  if (totalProb > m_prob_tmp) // fill solution with highest probability
1244  {
1245  sign_tmp = -log10(totalProb);
1246  m_prob_tmp = totalProb;
1247  m_fDitauStuffFit.Mditau_best = mtautau;
1248  m_fDitauStuffFit.Sign_best = sign_tmp;
1249  if (tau_ind == 1)
1250  m_fDitauStuffFit.nutau1 = nu_vec[j1];
1251  if (tau_ind == 2)
1252  m_fDitauStuffFit.nutau2 = nu_vec[j1];
1253  }
1254 
1255  ++ngoodsol1;
1256  }
1257 
1258  if (ngoodsol1 == 0)
1259  continue;
1260  m_iter2 += 1;
1261 
1262  m_iter3 += 1;
1263  }
1264  }
1265  }
1266  } else if (preparedInput.m_tauTypes == TauTypes::lh) // lepton+tau case
1267  {
1268  if (preparedInput.m_fUseVerbose == 1) {
1269  Info("DiTauMassTools", "Running in lepton+tau mode");
1270  }
1271  //------- Settings -------------------------------
1272 
1273  //----- Stuff below are for Winter 2012 lep-had analysis only; it has to be
1274  // replaced by a more common scheme once other channels are optimized
1275  // XYVector
1276  // mht_vec((tau_vec1+tau_vec2).Px(),(tau_vec1+tau_vec2).Py()); //
1277  // missing Ht vector for Njet25=0 events const double
1278  // mht=mht_vec.R();
1279  double input_metX = preparedInput.m_MetVec.X();
1280  double input_metY = preparedInput.m_MetVec.Y();
1281 
1282  // double mht_offset=0.0;
1283  // if(InputInfo.UseHT) // use missing Ht (for 0-jet events only for
1284  // now)
1285  // {
1286  // input_metX=-mht_vec.X();
1287  // input_metY=-mht_vec.Y();
1288  // }
1289  // else // use MET (for 0-jet and 1-jet events)
1290  // {
1291  // input_metX=met_vec.X();
1292  // input_metY=met_vec.Y();
1293  // }
1294 
1295  PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1296  PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1297  int tau_type_tmp;
1298  if (preparedInput.m_type_visTau1 == 8) {
1299  tau_tmp = preparedInput.m_vistau2;
1300  lep_tmp = preparedInput.m_vistau1;
1301  tau_type_tmp = preparedInput.m_type_visTau2;
1302  }
1303  if (preparedInput.m_type_visTau2 == 8) {
1304  tau_tmp = preparedInput.m_vistau1;
1305  lep_tmp = preparedInput.m_vistau2;
1306  tau_type_tmp = preparedInput.m_type_visTau1;
1307  }
1308 
1309  //---------------------------------------------
1310  for (int i4 = 0; i4 < NiterMET + 1; i4++) // MET_X scan
1311  {
1312  met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1313  for (int i5 = 0; i5 < NiterMET + 1; i5++) // MET_Y scan
1314  {
1315  met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1316  if (pow(met_smearL / METresX, 2) + pow(met_smearP / METresY, 2) > pow(N_METsigma, 2))
1317  continue; // use ellipse instead of square
1318  met_smear_x = met_smearL * m_metCovPhiCos - met_smearP * m_metCovPhiSin;
1319  met_smear_y = met_smearL * m_metCovPhiSin + met_smearP * m_metCovPhiCos;
1320  metvec_tmp.SetXY(input_metX + met_smear_x, input_metY + met_smear_y);
1321 
1322  solution = NuPsolutionLFV(metvec_tmp, tau_tmp, 0.0, nu_vec);
1323 
1324  ++iter0;
1325 
1326  if (solution < 1)
1327  continue;
1328  ++m_iter1;
1329 
1330  // if fast sin cos, result to not match exactly nupsolutionv2, so skip
1331  // test
1332  // SpeedUp no nested loop to compute individual probability
1333  int ngoodsol1 = 0;
1334 
1335  metprob = Prob->MetProbability(preparedInput, met_smearL, met_smearP, METresX, METresY);
1336  if (metprob <= 0)
1337  continue;
1338  for (unsigned int j1 = 0; j1 < nu_vec.size(); j1++) {
1339  if (tau_tmp.E() + nu_vec[j1].E() >= preparedInput.m_beamEnergy)
1340  continue;
1341  const double tau1_tmpp = (tau_tmp + nu_vec[j1]).P();
1342  angle1 = Angle(nu_vec[j1], tau_tmp);
1343 
1344  if (angle1 < dTheta3DLimit(tau_type_tmp, 0, tau1_tmpp)) {
1345  ++m_iang1low;
1346  continue;
1347  } // lower 99% bound
1348  if (angle1 > dTheta3DLimit(tau_type_tmp, 1, tau1_tmpp)) {
1349  ++m_iang1high;
1350  continue;
1351  } // upper 99% bound
1352  double tauvecprob1j =
1353  Prob->dTheta3d_probabilityFast(preparedInput, tau_type_tmp, angle1, tau1_tmpp);
1354  if (tauvecprob1j == 0.)
1355  continue;
1356  tauprob = Prob->TauProbabilityLFV(preparedInput, tau_type_tmp, tau_tmp, nu_vec[j1]);
1357  totalProb = tauvecprob1j * metprob * tauprob;
1358 
1359  m_tautau_tmp.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
1360  m_tautau_tmp += tau_tmp;
1361  m_tautau_tmp += lep_tmp;
1362  m_tautau_tmp += nu_vec[j1];
1363 
1364  const double mtautau = m_tautau_tmp.M();
1365 
1366  m_totalProbSum += totalProb;
1367  m_mtautauSum += mtautau;
1368 
1369  fit_code = 1; // at least one solution is found
1370 
1371  m_fMfit_all->Fill(mtautau, totalProb);
1372  m_fMfit_allNoWeight->Fill(mtautau, 1.);
1374  // m_fPXfit1->Fill((tau_tmp+nu_vec[j1]).Px(),totalProb);
1375  // m_fPYfit1->Fill((tau_tmp+nu_vec[j1]).Py(),totalProb);
1376  // m_fPZfit1->Fill((tau_tmp+nu_vec[j1]).Pz(),totalProb);
1377 
1378  if (totalProb > m_prob_tmp) // fill solution with highest probability
1379  {
1380  sign_tmp = -log10(totalProb);
1381  m_prob_tmp = totalProb;
1382  m_fDitauStuffFit.Mditau_best = mtautau;
1383  m_fDitauStuffFit.Sign_best = sign_tmp;
1384  if (preparedInput.m_type_visTau1 == 8) {
1385  m_fDitauStuffFit.vistau1 = lep_tmp;
1386  m_fDitauStuffFit.vistau2 = tau_tmp;
1387  m_fDitauStuffFit.nutau2 = nu_vec[j1];
1388  } else if (preparedInput.m_type_visTau2 == 8) {
1389  m_fDitauStuffFit.vistau2 = lep_tmp;
1390  m_fDitauStuffFit.vistau1 = tau_tmp;
1391  m_fDitauStuffFit.nutau1 = nu_vec[j1];
1392  }
1393  }
1394 
1395  ++ngoodsol1;
1396  }
1397 
1398  if (ngoodsol1 == 0)
1399  continue;
1400  m_iter2 += 1;
1401 
1402  m_iter3 += 1;
1403  }
1404  }
1405  } else {
1406  Info("DiTauMassTools", "Running in an unknown mode?!?!");
1407  }
1408 
1409  OutputInfo.m_NTrials = iter0;
1411 
1412  if (preparedInput.m_fUseVerbose == 1) {
1413  Info("DiTauMassTools", "%s",
1414  ("SpeedUp niters=" + std::to_string(iter0) + " " + std::to_string(m_iter1) + " " +
1416  " " + std::to_string(m_iang1high))
1417  .c_str());
1418  }
1419 
1420  if (m_fMfit_all->GetEntries() > 0 && m_iter3 > 0) {
1421 #ifdef SMOOTH
1422  m_fMfit_all->Smooth();
1423  m_fMfit_allNoWeight->Smooth();
1424  m_fPXfit1->Smooth();
1425  m_fPYfit1->Smooth();
1426  m_fPZfit1->Smooth();
1427 #endif
1428 
1429  // default max finding method defined in MissingMassCalculator.h
1430  // note that window defined in terms of number of bin, so depend on binning
1431  std::vector<double> histInfo(HistInfo::MAXHISTINFO);
1433  double prob_hist = histInfo.at(HistInfo::PROB);
1434 
1435  if (prob_hist != 0.0)
1436  m_fDitauStuffHisto.Sign_best = -log10(std::abs(prob_hist));
1437  else {
1438  // this mean the histogram is empty.
1439  // possible but very rare if all entries outside histogram range
1440  // fall back to maximum
1441  m_fDitauStuffHisto.Sign_best = -999.;
1443  }
1444 
1445  if (m_fDitauStuffHisto.Mditau_best > 0.0)
1447  std::vector<double> histInfoOther(HistInfo::MAXHISTINFO);
1448  //---- getting Nu1
1449  double Px1 = maxFromHist(m_fPXfit1, histInfoOther);
1450  double Py1 = maxFromHist(m_fPYfit1, histInfoOther);
1451  double Pz1 = maxFromHist(m_fPZfit1, histInfoOther);
1452  //---- setting 4-vecs
1453  PxPyPzMVector nu1_tmp(0.0, 0.0, 0.0, 0.0);
1454  PxPyPzMVector nu2_tmp(0.0, 0.0, 0.0, 0.0);
1455  if (preparedInput.m_type_visTau1 == 8) {
1456  nu1_tmp = preparedInput.m_vistau1;
1457  nu2_tmp.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
1458  }
1459  if (preparedInput.m_type_visTau2 == 8) {
1460  nu2_tmp = preparedInput.m_vistau2;
1461  nu1_tmp.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
1462  }
1465  }
1466  if (m_lfvLeplepRefit && fit_code==0 && !refit) {
1467  fit_code = DitauMassCalculatorV9lfv(true);
1468  return fit_code;
1469  }
1470 
1471 
1472 
1473  if (preparedInput.m_fUseVerbose == 1) {
1474  if (fit_code == 0) {
1475  Info(
1476  "DiTauMassTools", "%s",
1477  ("!!!----> Warning-3 in MissingMassCalculator::DitauMassCalculatorV9lfv() : fit status=" +
1478  std::to_string(fit_code))
1479  .c_str());
1480  Info("DiTauMassTools", "....... No solution is found. Printing input info .......");
1481 
1482  Info("DiTauMassTools", "%s", (" vis Tau-1: Pt="+std::to_string(preparedInput.m_vistau1.Pt())
1484  +" phi="+std::to_string(preparedInput.m_vistau1.Phi())
1485  +" type="+std::to_string(preparedInput.m_type_visTau1)).c_str());
1486  Info("DiTauMassTools", "%s", (" vis Tau-2: Pt="+std::to_string(preparedInput.m_vistau2.Pt())
1488  +" phi="+std::to_string(preparedInput.m_vistau2.Phi())
1489  +" type="+std::to_string(preparedInput.m_type_visTau2)).c_str());
1490  Info("DiTauMassTools", "%s", (" MET="+std::to_string(preparedInput.m_MetVec.R())+" Met_X="+std::to_string(preparedInput.m_MetVec.X())
1491  +" Met_Y="+std::to_string(preparedInput.m_MetVec.Y())).c_str());
1492  Info("DiTauMassTools", " ---------------------------------------------------------- ");
1493  }
1494  }
1495  return fit_code;
1496 }
1497 
1498 // function to fit maximum
1499 Double_t MissingMassCalculator::maxFitting(Double_t *x, Double_t *par)
1500 // Double_t maxFitting(Double_t *x, Double_t *par)
1501 {
1502  // parabola with parameters max, mean and invwidth
1503  const double mM = x[0];
1504  const double mMax = par[0];
1505  const double mMean = par[1];
1506  const double mInvWidth2 = par[2]; // if param positif distance between intersection of the
1507  // parabola with x axis: 1/Sqrt(mInvWidth2)
1508  const double fitval = mMax * (1 - 4 * mInvWidth2 * std::pow(mM - mMean, 2));
1509  return fitval;
1510 }
1511 
1512 // determine the maximum from the histogram
1513 // if input prob not default , compute also some probability
1514 // MaxHistStrategy : different method to find maximum
1515 // TODO should get the array on work on it
1516 // should also find the effective range of the hist
1517 
1518 double
1519 MissingMassCalculator::maxFromHist(TH1F *theHist, std::vector<double> &histInfo,
1520  const MaxHistStrategy::e maxHistStrategy,
1521  const int winHalfWidth, bool debug) {
1522  // namespace HistInfo
1523  // enum e {
1524  // PROB=0,INTEGRAL,CHI2,DISCRI,TANTHETA,TANTHETAW,FITLENGTH,RMS,RMSVSDISCRI,MAXHISTINFO
1525  // };
1526  double maxPos = 0.;
1527  double prob = 0.;
1528 
1529  for (std::vector<double>::iterator itr = histInfo.begin(); itr != histInfo.end(); ++itr) {
1530  *itr = -1;
1531  }
1532 
1533  histInfo[HistInfo::INTEGRAL] = theHist->Integral();
1534 
1535  if (maxHistStrategy == MaxHistStrategy::MAXBIN ||
1536  ((maxHistStrategy == MaxHistStrategy::MAXBINWINDOW ||
1537  maxHistStrategy == MaxHistStrategy::SLIDINGWINDOW) &&
1538  winHalfWidth == 0)) {
1539 
1540  // simple max search
1541  // original version, simple bin maximum
1542  int max_bin = theHist->GetMaximumBin();
1543  maxPos = theHist->GetBinCenter(max_bin);
1544 
1545  // FIXME GetEntries is unweighted
1546  prob = theHist->GetBinContent(max_bin) / double(theHist->GetEntries());
1547  if (prob > 1.)
1548  prob = 1.;
1549  histInfo[HistInfo::PROB] = prob;
1550  return maxPos;
1551  }
1552 
1553  int hNbins = theHist->GetNbinsX();
1554 
1555  if (maxHistStrategy == MaxHistStrategy::MAXBINWINDOW) {
1556  // average around maximum bin (nearly useless in fact)
1557  // could be faster
1558  int max_bin = theHist->GetMaximumBin();
1559  int iBinMin = max_bin - winHalfWidth;
1560  if (iBinMin < 0)
1561  iBinMin = 0;
1562  int iBinMax = max_bin + winHalfWidth;
1563  if (iBinMax > hNbins)
1564  iBinMax = hNbins - 1;
1565  double sumw = 0;
1566  double sumx = 0;
1567  for (int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1568  const double weight = theHist->GetBinContent(iBin);
1569  sumw += weight;
1570  sumx += weight * theHist->GetBinCenter(iBin);
1571  }
1572  maxPos = sumx / sumw;
1573 
1574  // FIXME GetEntries is unweighted
1575  prob = sumw / theHist->GetEntries();
1576  if (prob > 1.)
1577  prob = 1.;
1578 
1579  return maxPos;
1580  }
1581 
1582  // now compute sliding window anyway
1583  if (maxHistStrategy != MaxHistStrategy::SLIDINGWINDOW &&
1584  maxHistStrategy != MaxHistStrategy::FIT) {
1585  Error("DiTauMassTools", "%s",
1586  ("ERROR undefined maxHistStrategy:" + std::to_string(maxHistStrategy)).c_str());
1587  return -10.;
1588  }
1589 
1590  // first iteration to find the first and last non zero bin, and the histogram
1591  // integral (not same as Entries because of weights)
1592  int lastNonZeroBin = -1;
1593  int firstNonZeroBin = -1;
1594  double totalSumw = 0.;
1595  bool firstNullPart = true;
1596  for (int iBin = 0; iBin < hNbins; ++iBin) {
1597  const double weight = theHist->GetBinContent(iBin);
1598  if (weight > 0) {
1599  totalSumw += weight;
1600  lastNonZeroBin = iBin;
1601  if (firstNullPart) {
1602  firstNullPart = false;
1603  firstNonZeroBin = iBin;
1604  }
1605  }
1606  }
1607 
1608  // enlarge first and last non zero bin with window width to avoid side effect
1609  // (maximum close to the edge)
1610  firstNonZeroBin = std::max(0, firstNonZeroBin - winHalfWidth - 1);
1611  lastNonZeroBin = std::min(hNbins - 1, lastNonZeroBin + winHalfWidth + 1);
1612 
1613  // if null histogram quit
1614  if (firstNullPart)
1615  return maxPos;
1616 
1617  // determine the size of the sliding window in the fit case
1618 
1619  // sliding window
1620  const int nwidth = 2 * winHalfWidth + 1;
1621  double winsum = 0.;
1622 
1623  for (int ibin = 0; ibin < nwidth; ++ibin) {
1624  winsum += theHist->GetBinContent(ibin);
1625  }
1626  double winmax = winsum;
1627 
1628  int max_bin = 0.;
1629  int iBinL = firstNonZeroBin;
1630  int iBinR = iBinL + 2 * winHalfWidth;
1631  bool goingUp = true;
1632 
1633  do {
1634  ++iBinL;
1635  ++iBinR;
1636  const double deltawin = theHist->GetBinContent(iBinR) - theHist->GetBinContent(iBinL - 1);
1637 
1638  if (deltawin < 0) {
1639  if (goingUp) {
1640  // if were climbing and now loose more on the left
1641  // than win on the right. This was a local maxima
1642  if (winsum > winmax) {
1643  // global maximum one so far
1644  winmax = winsum;
1645  max_bin = (iBinR + iBinL) / 2 - 1;
1646  }
1647  goingUp = false; // now going down
1648  }
1649  } else {
1650  // do not care about minima, simply indicate we are going down
1651  goingUp = true;
1652  }
1653 
1654  winsum += deltawin;
1655 
1656  } while (iBinR < lastNonZeroBin);
1657 
1658  // now compute average
1659  int iBinMin = max_bin - winHalfWidth;
1660  if (iBinMin < 0)
1661  iBinMin = 0;
1662  int iBinMax = max_bin + winHalfWidth;
1663  if (iBinMax >= hNbins)
1664  iBinMax = hNbins - 1;
1665  double sumw = 0;
1666  double sumx = 0;
1667  for (int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1668  const double weight = theHist->GetBinContent(iBin);
1669  sumw += weight;
1670  sumx += weight * theHist->GetBinCenter(iBin);
1671  }
1672 
1673  double maxPosWin = -1.;
1674 
1675  if (sumw > 0.) {
1676  maxPosWin = sumx / sumw;
1677  }
1678  // prob if the fraction of events in the window
1679  prob = sumw / totalSumw;
1680 
1681  // Definitions of some useful parameters
1682 
1683  const double h_rms = theHist->GetRMS(1);
1684  histInfo[HistInfo::RMS] = h_rms;
1685 
1686  double num = 0;
1687  double numerator = 0;
1688  double denominator = 0;
1689  bool nullBin = false;
1690 
1691  for (int i = iBinMin; i < iBinMax; ++i) {
1692  double binError = theHist->GetBinError(i);
1693  if (binError < 1e-10) {
1694  nullBin = true;
1695  }
1696  double binErrorSquare = std::pow(binError, 2);
1697  num = theHist->GetBinContent(i) / (binErrorSquare);
1698  numerator = numerator + num;
1699  denominator = denominator + (1 / (binErrorSquare));
1700  }
1701  if (numerator < 1e-10 || denominator < 1e-10 || nullBin == true) {
1702  histInfo[HistInfo::MEANBIN] = -1;
1703  } else {
1704  histInfo[HistInfo::MEANBIN] = sqrt(1 / denominator) / (numerator / denominator);
1705  }
1706 
1707  // stop here if only looking for sliding window
1708  if (maxHistStrategy == MaxHistStrategy::SLIDINGWINDOW) {
1709  return maxPosWin;
1710  }
1711 
1712  maxPos = maxPosWin;
1713  // now FIT maxHistStrategy==MaxHistStrategy::FIT
1714 
1715  // now mass fit in range defined by sliding window
1716  // window will be around maxPos
1717  const double binWidth = theHist->GetBinCenter(2) - theHist->GetBinCenter(1);
1718  double fitWidth = (winHalfWidth + 0.5) * binWidth;
1719  // fit range 2 larger than original window range, 3 if less than 20% of the
1720  // histogram in slinding window
1721 
1722  if (prob > 0.2) {
1723  fitWidth *= 2.;
1724  } else {
1725  fitWidth *= 3.;
1726  }
1727  // fit option : Q == Quiet, no printout S result of the fit returned in
1728  // TFitResultPtr N do not draw the resulting function
1729 
1730  // if debug plot the fitted function
1731  TString fitOption = debug ? "QS" : "QNS";
1732  // root fit
1733  // Sets initial values
1734  m_fFitting->SetParameters(sumw / winHalfWidth, maxPos, 0.0025);
1735  // TFitResultPtr
1736  // fitRes=theHist->Fit("pol2",fitOption,"",maxPos-fitWidth,maxPos+fitWidth);
1737  TFitResultPtr fitRes =
1738  theHist->Fit(m_fFitting, fitOption, "", maxPos - fitWidth, maxPos + fitWidth);
1739 
1740  double maxPosFit = -1.;
1741 
1742  if (int(fitRes) == 0) {
1743  // root fit
1744  histInfo[HistInfo::CHI2] = fitRes->Chi2();
1745  const double mMax = fitRes->Parameter(0);
1746  const double mMean = fitRes->Parameter(1);
1747  const double mInvWidth2 = fitRes->Parameter(2);
1748  double mMaxError = fitRes->ParError(0);
1749  m_PrintmMaxError = mMaxError;
1750  double mMeanError = fitRes->ParError(1);
1751  m_PrintmMeanError = mMeanError;
1752  double mInvWidth2Error = fitRes->ParError(2);
1753  m_PrintmInvWidth2Error = mInvWidth2Error;
1754  mMeanError = 0.; // avoid warning
1755  mInvWidth2Error = 0.; // avoid warning
1756  const double c = mMax * (1 - 4 * mMean * mMean * mInvWidth2);
1757  const double b = 8 * mMax * mMean * mInvWidth2;
1758  const double a = -4 * mMax * mInvWidth2;
1759  // when built in polynomial fit
1760  // const double c=fitRes->Parameter(0);
1761  // const double b=fitRes->Parameter(1);
1762  // const double a=fitRes->Parameter(2);
1763 
1764  const double h_discri = b * b - 4 * a * c;
1765  histInfo[HistInfo::DISCRI] = h_discri;
1766  const double sqrth_discri = sqrt(h_discri);
1767  const double h_fitLength = sqrth_discri / a;
1768  histInfo[HistInfo::FITLENGTH] = h_fitLength;
1769  histInfo[HistInfo::TANTHETA] = 2 * a / sqrth_discri;
1770  histInfo[HistInfo::TANTHETAW] = 2 * a * sumw / sqrth_discri;
1771  histInfo[HistInfo::RMSVSDISCRI] = h_rms / h_fitLength;
1772  // compute maximum position (only if inverted parabola)
1773  if (a < 0)
1774  maxPosFit = -b / (2 * a);
1775  }
1776 
1777  // keep fit result only if within 80% of fit window, and fit succeeded
1778  if (maxPosFit >= 0. and std::abs(maxPosFit - maxPosWin) < 0.8 * fitWidth) {
1779  histInfo[HistInfo::PROB] = prob;
1780  return maxPosFit;
1781  } else {
1782  // otherwise keep the weighted average
1783  // negate prob just to flag such event
1784  prob = -prob;
1785  histInfo[HistInfo::PROB] = prob;
1786  return maxPosWin;
1787  }
1788 }
1789 
1790 // compute probability for any input value,can be called from a pure parameter
1791 // scan
1792 // deltametvec is along phijet
1793 // returns number of solution if positive, return code if negative, vector of
1794 // probability and mass
1795 int MissingMassCalculator::probCalculatorV9fast(const double &phi1, const double &phi2,
1796  const double &M_nu1,
1797  const double &M_nu2) {
1798  // bool debug=true;
1799 
1800  int nsol1;
1801  int nsol2;
1802 
1803  const int solution = NuPsolutionV3(M_nu1, M_nu2, phi1, phi2, nsol1, nsol2);
1804 
1805  if (solution != 1)
1806  return -4;
1807  // refineSolutions ( M_nu1,M_nu2,
1808  // met_smearL,met_smearP,metvec_tmp.R(),
1809  // nsol1, nsol2,m_Mvis,m_Meff);
1810  refineSolutions(M_nu1, M_nu2, nsol1, nsol2, m_Mvis, m_Meff);
1811 
1812  if (m_nsol <= 0)
1813  return 0;
1814 
1815  // success
1816 
1817  return m_nsol; // for backward compatibility
1818 }
1819 
1820 // nuvecsol1 and nuvecsol2 passed by MMC
1821 int MissingMassCalculator::refineSolutions(const double &M_nu1, const double &M_nu2,
1822  const int nsol1, const int nsol2,
1823  const double &Mvis, const double &Meff)
1824 
1825 {
1826  m_nsol = 0;
1827 
1828  if (int(m_probFinalSolVec.size()) < m_nsolfinalmax)
1829  Error("DiTauMassTools", "%s",
1830  ("refineSolutions ERROR probFinalSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1831  .c_str());
1832  if (int(m_mtautauFinalSolVec.size()) < m_nsolfinalmax)
1833  Error("DiTauMassTools", "%s",
1834  ("refineSolutions ERROR mtautauSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1835  .c_str());
1836  if (int(m_nu1FinalSolVec.size()) < m_nsolfinalmax)
1837  Error("DiTauMassTools", "%s",
1838  ("refineSolutions ERROR nu1FinalSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1839  .c_str());
1840  if (int(m_nu2FinalSolVec.size()) < m_nsolfinalmax)
1841  Error("DiTauMassTools", "%s",
1842  ("refineSolutions ERROR nu2FinalSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1843  .c_str());
1844  if (nsol1 > int(m_nsolmax))
1845  Error("DiTauMassTools", "%s", ("refineSolutions ERROR nsol1 " + std::to_string(nsol1) +
1846  " > nsolmax !" + std::to_string(m_nsolmax))
1847  .c_str());
1848  if (nsol2 > int(m_nsolmax))
1849  Error("DiTauMassTools", "%s", ("refineSolutions ERROR nsol1 " + std::to_string(nsol2) +
1850  " > nsolmax !" + std::to_string(m_nsolmax))
1851  .c_str());
1852 
1853  int ngoodsol1 = 0;
1854  int ngoodsol2 = 0;
1855  double constProb =
1856  Prob->apply(preparedInput, -99, -99, PtEtaPhiMVector(0, 0, 0, 0), PtEtaPhiMVector(0, 0, 0, 0),
1857  PtEtaPhiMVector(0, 0, 0, 0), PtEtaPhiMVector(0, 0, 0, 0), true, false, false);
1858 
1859  for (int j1 = 0; j1 < nsol1; ++j1) {
1860  PtEtaPhiMVector &nuvec1_tmpj = m_nuvecsol1[j1];
1861  PtEtaPhiMVector &tauvecsol1j = m_tauvecsol1[j1];
1862  double &tauvecprob1j = m_tauvecprob1[j1];
1863  tauvecprob1j = 0.;
1864  // take first or second solution
1865  // no time to call rndm, switch more or less randomely, according to an
1866  // oscillating switch perturbed by m_phi1
1867  if (nsol1 > 1) {
1868  if (j1 == 0) { // decide at the first solution which one we will take
1869  const int pickInt = std::abs(10000 * m_Phi1);
1870  const int pickDigit = pickInt - 10 * (pickInt / 10);
1871  if (pickDigit < 5)
1872  m_switch1 = !m_switch1;
1873  }
1874  m_switch1 = !m_switch1;
1875  }
1876 
1877  if (!m_switch1) {
1878  nuvec1_tmpj.SetCoordinates(nuvec1_tmpj.Pt(), nuvec1_tmpj.Eta(), nuvec1_tmpj.Phi(), M_nu1);
1879  tauvecsol1j.SetPxPyPzE(0., 0., 0., 0.);
1880  tauvecsol1j += nuvec1_tmpj;
1881  tauvecsol1j += m_tauVec1;
1882  if (tauvecsol1j.E() >= preparedInput.m_beamEnergy)
1883  continue;
1885  PtEtaPhiMVector(0, 0, 0, 0), nuvec1_tmpj,
1886  PtEtaPhiMVector(0, 0, 0, 0), false, true, false);
1887  ++ngoodsol1;
1888  }
1889 
1890  for (int j2 = 0; j2 < nsol2; ++j2) {
1891  PtEtaPhiMVector &nuvec2_tmpj = m_nuvecsol2[j2];
1892  PtEtaPhiMVector &tauvecsol2j = m_tauvecsol2[j2];
1893  double &tauvecprob2j = m_tauvecprob2[j2];
1894  if (j1 == 0) {
1895  tauvecprob2j = 0.;
1896  // take first or second solution
1897  // no time to call rndm, switch more or less randomely, according to an
1898  // oscillating switch perturbed by m_phi2
1899  if (nsol2 > 1) {
1900  if (j2 == 0) { // decide at the first solution which one we will take
1901  const int pickInt = std::abs(10000 * m_Phi2);
1902  const int pickDigit = pickInt - 10 * int(pickInt / 10);
1903  if (pickDigit < 5)
1904  m_switch2 = !m_switch2;
1905  }
1906  m_switch2 = !m_switch2;
1907  }
1908 
1909  if (!m_switch2) {
1910  nuvec2_tmpj.SetCoordinates(nuvec2_tmpj.Pt(), nuvec2_tmpj.Eta(), nuvec2_tmpj.Phi(), M_nu2);
1911  tauvecsol2j.SetPxPyPzE(0., 0., 0., 0.);
1912  tauvecsol2j += nuvec2_tmpj;
1913  tauvecsol2j += m_tauVec2;
1914  if (tauvecsol2j.E() >= preparedInput.m_beamEnergy)
1915  continue;
1916  tauvecprob2j = Prob->apply(preparedInput, -99, preparedInput.m_type_visTau2,
1917  PtEtaPhiMVector(0, 0, 0, 0), m_tauVec2,
1918  PtEtaPhiMVector(0, 0, 0, 0), nuvec2_tmpj, false, true, false);
1919  ++ngoodsol2;
1920  }
1921  }
1922  if (tauvecprob1j == 0.)
1923  continue;
1924  if (tauvecprob2j == 0.)
1925  continue;
1926 
1927  double totalProb = 1.;
1928 
1929  m_tautau_tmp.SetPxPyPzE(0., 0., 0., 0.);
1930  m_tautau_tmp += tauvecsol1j;
1931  m_tautau_tmp += tauvecsol2j;
1932  const double mtautau = m_tautau_tmp.M();
1933 
1934  if (TailCleanUp(m_tauVec1, nuvec1_tmpj, m_tauVec2, nuvec2_tmpj, mtautau, Mvis, Meff,
1935  preparedInput.m_DelPhiTT) == 0) {
1936  continue;
1937  }
1938 
1939  totalProb *=
1940  (constProb * tauvecprob1j * tauvecprob2j *
1942  m_tauVec1, m_tauVec2, nuvec1_tmpj, nuvec2_tmpj, false, false, true));
1943 
1944  if (totalProb <= 0) {
1946  Warning("DiTauMassTools", "%s",
1947  ("null proba solution, rejected "+std::to_string(totalProb)).c_str());
1948  } else {
1949  // only count solution with non zero probability
1950  m_totalProbSum += totalProb;
1951  m_mtautauSum += mtautau;
1952 
1953  if (m_nsol >= int(m_nsolfinalmax)) {
1954  Error("DiTauMassTools", "%s",
1955  ("refineSolutions ERROR nsol getting larger than nsolfinalmax!!! " +
1957  .c_str());
1958  Error("DiTauMassTools", "%s",
1959  (" j1 " + std::to_string(j1) + " j2 " + std::to_string(j2) + " nsol1 " +
1960  std::to_string(nsol1) + " nsol2 " + std::to_string(nsol2))
1961  .c_str());
1962  --m_nsol; // overwrite last solution. However this should really never
1963  // happen
1964  }
1965 
1966  // good solution found, copy in vector
1967  m_mtautauFinalSolVec[m_nsol] = mtautau;
1968  m_probFinalSolVec[m_nsol] = totalProb;
1969 
1970  PtEtaPhiMVector &nu1Final = m_nu1FinalSolVec[m_nsol];
1971  PtEtaPhiMVector &nu2Final = m_nu2FinalSolVec[m_nsol];
1972  // for (int iv=0;iv<4;++iv){
1973 
1974  nu1Final.SetPxPyPzE(nuvec1_tmpj.Px(), nuvec1_tmpj.Py(), nuvec1_tmpj.Pz(), nuvec1_tmpj.E());
1975  nu2Final.SetPxPyPzE(nuvec2_tmpj.Px(), nuvec2_tmpj.Py(), nuvec2_tmpj.Pz(), nuvec2_tmpj.E());
1976  // }
1977 
1978  ++m_nsol;
1979  } // else totalProb<=0
1980 
1981  } // loop j2
1982  } // loop j1
1983  if (ngoodsol1 == 0) {
1984  return -1;
1985  }
1986  if (ngoodsol2 == 0) {
1987  return -2;
1988  }
1989  return m_nsol;
1990 }
1991 
1992 int MissingMassCalculator::TailCleanUp(const PtEtaPhiMVector &vis1,
1993  const PtEtaPhiMVector &nu1,
1994  const PtEtaPhiMVector &vis2,
1995  const PtEtaPhiMVector &nu2, const double &mmc_mass,
1996  const double &vis_mass, const double &eff_mass,
1997  const double &dphiTT) {
1998 
1999  int pass_code = 1;
2001  return pass_code;
2002 
2003  // the Clean-up cuts are specifically for rel16 analyses.
2004  // the will change in rel17 analyses and after the MMC is updated
2005 
2006  if (preparedInput.m_tauTypes == TauTypes::ll) // lepton-lepton channel
2007  {
2008  const double MrecoMvis = mmc_mass / vis_mass;
2009  if (MrecoMvis > 2.6)
2010  return 0;
2011  const double MrecoMeff = mmc_mass / eff_mass;
2012  if (MrecoMeff > 1.9)
2013  return 0;
2014  const double e1p1 = nu1.E() / vis1.P();
2015  const double e2p2 = nu2.E() / vis2.P();
2016  if ((e1p1 + e2p2) > 4.5)
2017  return 0;
2018  if (e2p2 > 4.0)
2019  return 0;
2020  if (e1p1 > 3.0)
2021  return 0;
2022  }
2023 
2024  //-------- these are new cuts for lep-had analysis for Moriond
2025  if (preparedInput.m_tauTypes == TauTypes::lh) // lepton-hadron channel
2026  {
2027 
2033  return pass_code; // don't use TailCleanup for 8 & 13 TeV data
2034 
2035  //--------- leave code uncommented to avoid Compilation warnings
2036  if (Prob->GetUseHT()) {
2037  const double MrecoMvis = mmc_mass / vis_mass;
2038  const double MrecoMeff = mmc_mass / eff_mass;
2039  const double x = dphiTT > 1.5 ? dphiTT : 1.5;
2040  if ((MrecoMeff + MrecoMvis) > 5.908 - 1.881 * x + 0.2995 * x * x)
2041  return 0;
2042  }
2043  }
2044  return pass_code;
2045 }
2046 
2047 // note that if MarkovChain the input solutions can be modified
2049 
2050 {
2051 
2052  bool reject = true;
2053  double totalProbSumSol = 0.;
2054  double totalProbSumSolOld = 0.;
2055  bool firstPointWithSol = false;
2056 
2057  for (int isol = 0; isol < m_nsol; ++isol) {
2058  totalProbSumSol += m_probFinalSolVec[isol];
2059  }
2060 
2061  double uMC = -1.;
2062  bool notSureToKeep = true;
2063  // note : if no solution, the point is treated as having a zero probability
2064  if (m_fullParamSpaceScan) {
2065  reject = false; // accept anyway in this mode
2066  notSureToKeep = false; // do not need to test on prob
2067  if (m_nsol <= 0) {
2068  // if initial full scaning and no sol : continue
2069  m_markovNFullScan += 1;
2070  } else {
2071  // if we were in in full scan mode and we have a solution, switch it off
2072  m_fullParamSpaceScan = false;
2073  firstPointWithSol = true; // as this is the first point without a solution
2074  // there is no old sol
2075  m_iter0 = 0; // reset the counter so that separately the full scan pphase
2076  // and the markov phase use m_niterRandomLocal points
2077  // hack for hh : allow 10 times less iteration for markov than for the
2078  // fullscan phase
2080  m_niterRandomLocal /= 10;
2081  }
2082  }
2083  }
2084 
2085  if (notSureToKeep) {
2086  // apply Metropolis algorithm to decide to keep this point.
2087  // compute the probability of the previous point and the current one
2088  for (int isol = 0; isol < m_nsolOld; ++isol) {
2089  totalProbSumSolOld += m_probFinalSolOldVec[isol];
2090  }
2091 
2092  // accept anyway if null old probability (should only happen for the very
2093  // first point with a solution)
2094  if (!firstPointWithSol && totalProbSumSolOld <= 0.) {
2095  Error("DiTauMassTools", "%s",
2096  (" ERROR null old probability !!! " + std::to_string(totalProbSumSolOld) + " nsolOld " +
2098  .c_str());
2099  reject = false;
2100  } else if (totalProbSumSol > totalProbSumSolOld) {
2101  // if going up, accept anyway
2102  reject = false;
2103  // else if (totalProbSumSol < 1E-16) { // if null target probability,
2104  // reject anyway
2105  } else if (totalProbSumSol < totalProbSumSolOld * 1E-6) { // if ratio of probability <1e6, point
2106  // will be accepted only every 1E6
2107  // iteration, so can reject anyway
2108  reject = true;
2109  } else if (m_nsol <= 0) { // new parametrisation give prob too small to
2110  // trigger above condition if no solution is found
2111  reject = true;
2112  } else {
2113  // if going down, reject with a probability
2114  // 1-totalProbSum/totalProbSumOld)
2115  uMC = m_randomGen.Rndm();
2116  reject = (uMC > totalProbSumSol / totalProbSumSolOld);
2117  }
2118  } // if reject
2119 
2120  // proceed with the handling of the solutions wether the old or the new ones
2121 
2122  // optionally fill the vectors with the complete list of points (for all
2123  // walkstrategy)
2124 
2125  if (reject) {
2126  // current point reset to the previous one
2127  // Note : only place where m_MEtP etc... are modified outside spacewalkerXYZ
2128  m_MEtP = m_MEtP0;
2129  m_MEtL = m_MEtL0;
2130  m_Phi1 = m_Phi10;
2131  m_Phi2 = m_Phi20;
2132  m_eTau1 = m_eTau10;
2133  m_eTau2 = m_eTau20;
2134  if (m_scanMnu1)
2135  m_Mnu1 = m_Mnu10;
2136  if (m_scanMnu2)
2137  m_Mnu2 = m_Mnu20;
2138  }
2139 
2140  // default case : fill the histogram with solution, using current point
2141  bool fillSolution = true;
2142  bool oldToBeUsed = false;
2143 
2144  // now handle the reject or accept cases
2145  // the tricky thing is that for markov, we accept the old point as soon as a
2146  // new accepted point is found with a weight equal to one plus the number of
2147  // rejected point inbetween
2148 
2149  if (reject) {
2150  fillSolution = false; // do not fill solution, just count number of replication
2152  if (m_nsol <= 0) {
2153  m_markovNRejectNoSol += 1;
2154  } else {
2156  }
2157 
2158  } else {
2159  // if accept, will fill solution (except for very first point) but taking
2160  // the values from the previous point
2161  if (!m_fullParamSpaceScan) {
2162  m_markovNAccept += 1;
2163  }
2164  if (!firstPointWithSol) {
2165  fillSolution = true;
2166  oldToBeUsed = true;
2167  } else {
2168  fillSolution = false;
2169  }
2170  } // else reject
2171 
2172  // if do not fill solution exit now
2173  // for the first point with solution we need to copy the new sol into the old
2174  // one before leaving
2175  if (!fillSolution) {
2176  if (firstPointWithSol) {
2177  // current point is the future previous one
2178  m_nsolOld = m_nsol;
2179  for (int isol = 0; isol < m_nsol; ++isol) {
2182  m_nu1FinalSolOldVec[isol] = m_nu1FinalSolVec[isol];
2183  m_nu2FinalSolOldVec[isol] = m_nu2FinalSolVec[isol];
2184  }
2185  }
2186  return;
2187  }
2188 
2189  // compute RMS of the different solutions
2190  double solSum = 0.;
2191  double solSum2 = 0.;
2192 
2193  for (int isol = 0; isol < m_nsol; ++isol) {
2194  ++m_iter5;
2195  double totalProb;
2196  double mtautau;
2197  const PtEtaPhiMVector *pnuvec1_tmpj;
2198  const PtEtaPhiMVector *pnuvec2_tmpj;
2199 
2200  if (oldToBeUsed) {
2201  totalProb = m_probFinalSolOldVec[isol];
2202  mtautau = m_mtautauFinalSolOldVec[isol];
2203  pnuvec1_tmpj = &m_nu1FinalSolOldVec[isol];
2204  pnuvec2_tmpj = &m_nu2FinalSolOldVec[isol];
2205  } else {
2206  totalProb = m_probFinalSolVec[isol];
2207  mtautau = m_mtautauFinalSolVec[isol];
2208  pnuvec1_tmpj = &m_nu1FinalSolVec[isol];
2209  pnuvec2_tmpj = &m_nu2FinalSolVec[isol];
2210  }
2211  const PtEtaPhiMVector &nuvec1_tmpj = *pnuvec1_tmpj;
2212  const PtEtaPhiMVector &nuvec2_tmpj = *pnuvec2_tmpj;
2213 
2214  solSum += mtautau;
2215  solSum2 += mtautau * mtautau;
2216 
2217  double weight;
2218  // MarkovChain : accepted events already distributed according to
2219  // probability distribution, so weight is 1. acutally to have a proper
2220  // estimate of per bin error, instead of putting several time the same point
2221  // when metropolis alg reject one (or no solution), rather put it with the
2222  // multiplicity weight. Should only change the error bars might change if
2223  // weighted markov chain are used there is also an issue with the 4 very
2224  // close nearly identical solution
2226  1; // incremented only when a point is rejected, hence need to add 1
2227 
2228  m_fMfit_all->Fill(mtautau, weight);
2229 
2230  if(m_SaveLlhHisto){
2231  m_fMEtP_all->Fill(m_MEtP, weight);
2232  m_fMEtL_all->Fill(m_MEtL, weight);
2233  m_fMnu1_all->Fill(m_Mnu1, weight);
2234  m_fMnu2_all->Fill(m_Mnu2, weight);
2235  m_fPhi1_all->Fill(m_Phi1, weight);
2236  m_fPhi2_all->Fill(m_Phi2, weight);
2237  if (mtautau != 0. && weight != 0.)
2238  m_fMfit_allGraph->SetPoint(m_iter0, mtautau, -TMath::Log(weight));
2239  }
2240 
2241  m_fMfit_allNoWeight->Fill(mtautau, 1.);
2242 
2243  // m_fPXfit1->Fill(nuvec1_tmpj.Px(),weight);
2244  // m_fPYfit1->Fill(nuvec1_tmpj.Py(),weight);
2245  // m_fPZfit1->Fill(nuvec1_tmpj.Pz(),weight);
2246  // m_fPXfit2->Fill(nuvec2_tmpj.Px(),weight);
2247  // m_fPYfit2->Fill(nuvec2_tmpj.Py(),weight);
2248  // m_fPZfit2->Fill(nuvec2_tmpj.Pz(),weight);
2249 
2250  //----------------- using P*fit to fill Px,y,z_tau
2251  // Note that the original vistau are used there deliberately,
2252  // since they will be subtracted after histogram fitting
2253  // DR, kudos Antony Lesage : do not create temporary TLV within each Fill,
2254  // saves 10% CPU
2255  m_fPXfit1->Fill(preparedInput.m_vistau1.Px() + nuvec1_tmpj.Px(), totalProb);
2256  m_fPYfit1->Fill(preparedInput.m_vistau1.Py() + nuvec1_tmpj.Py(), totalProb);
2257  m_fPZfit1->Fill(preparedInput.m_vistau1.Pz() + nuvec1_tmpj.Pz(), totalProb);
2258  m_fPXfit2->Fill(preparedInput.m_vistau2.Px() + nuvec2_tmpj.Px(), totalProb);
2259  m_fPYfit2->Fill(preparedInput.m_vistau2.Py() + nuvec2_tmpj.Py(), totalProb);
2260  m_fPZfit2->Fill(preparedInput.m_vistau2.Pz() + nuvec2_tmpj.Pz(), totalProb);
2261 
2262  // fill histograms for floating stopping criterion, split randomly
2263  if (m_fUseFloatStopping) {
2264  if (m_randomGen.Rndm() <= 0.5) {
2265  m_fMmass_split1->Fill(mtautau, weight);
2266  m_fMEtP_split1->Fill(m_MEtP, weight);
2267  m_fMEtL_split1->Fill(m_MEtL, weight);
2268  m_fMnu1_split1->Fill(m_Mnu1, weight);
2269  m_fMnu2_split1->Fill(m_Mnu2, weight);
2270  m_fPhi1_split1->Fill(m_Phi1, weight);
2271  m_fPhi2_split1->Fill(m_Phi2, weight);
2272  } else {
2273  m_fMmass_split2->Fill(mtautau, weight);
2274  m_fMEtP_split2->Fill(m_MEtP, weight);
2275  m_fMEtL_split2->Fill(m_MEtL, weight);
2276  m_fMnu1_split2->Fill(m_Mnu1, weight);
2277  m_fMnu2_split2->Fill(m_Mnu2, weight);
2278  m_fPhi1_split2->Fill(m_Phi1, weight);
2279  m_fPhi2_split2->Fill(m_Phi2, weight);
2280  }
2281  }
2282 
2283  if (totalProb > m_prob_tmp) // fill solution with highest probability
2284  {
2285  m_prob_tmp = totalProb;
2286  m_fDitauStuffFit.Mditau_best = mtautau;
2287  m_fDitauStuffFit.Sign_best = -log10(totalProb);
2288  ;
2289  m_fDitauStuffFit.nutau1 = nuvec1_tmpj;
2290  m_fDitauStuffFit.nutau2 = nuvec2_tmpj;
2293  }
2294  } // loop on solutions
2295 
2296  m_markovCountDuplicate = 0; // now can reset the duplicate count
2297 
2298  if (oldToBeUsed) {
2299  // current point is the future previous one
2300  // TLV copy not super efficient but not dramatic
2301  m_nsolOld = m_nsol;
2302  for (int isol = 0; isol < m_nsol; ++isol) {
2305  m_nu1FinalSolOldVec[isol] = m_nu1FinalSolVec[isol];
2306  m_nu2FinalSolOldVec[isol] = m_nu2FinalSolVec[isol];
2307  }
2308  }
2309 
2310  // compute rms of solutions
2311  const double solRMS = sqrt(solSum2 / m_nsol - std::pow(solSum / m_nsol, 2));
2312  OutputInfo.m_AveSolRMS += solRMS;
2313 
2314  return;
2315 }
2316 
2318  // FIXME could use function pointer to switch between functions
2319  m_nsolOld = 0;
2320 
2321  double METresX = preparedInput.m_METsigmaL; // MET resolution in direction parallel to MET
2322  // resolution major axis, for MET scan
2323  double METresY = preparedInput.m_METsigmaP; // MET resolution in direction perpendicular to
2324  // to MET resolution major axis, for MET scan
2325 
2326  // precompute some quantities and store in m_ data members
2327  precomputeCache();
2334  }
2335 
2336  // if m_nsigma_METscan was not set by user, set to default values
2337  if(m_nsigma_METscan == -1){
2338  if (preparedInput.m_tauTypes == TauTypes::ll) // both tau's are leptonic
2339  {
2341  } else if (preparedInput.m_tauTypes == TauTypes::lh) // lep had
2342  {
2344  } else // hh
2345  {
2347  }
2348  }
2349 
2351 
2354 
2355  m_walkWeight = 1.;
2356 
2357  // dummy initial value to avoid printout with random values
2358  m_Phi10 = 0.;
2359  m_Phi20 = 0.;
2360  m_MEtL0 = 0.;
2361  m_MEtP0 = 0.;
2362  m_Mnu10 = 0.;
2363  m_Mnu20 = 0.;
2364 
2366 
2367  // seeds the random generator in a reproducible way from the phi of both tau;
2368  double aux = std::abs(m_tauVec1Phi + double(m_tauVec2Phi) / 100. / TMath::Pi()) * 100;
2369  m_seed = (aux - floor(aux)) * 1E6 * (1 + m_RndmSeedAltering) + 13;
2370 
2371  m_randomGen.SetSeed(m_seed);
2372  // int Niter=Niter_fit1; // number of points for each dR loop
2373  // int NiterMET=Niter_fit2; // number of iterations for each MET scan loop
2374  // int NiterMnu=Niter_fit3; // number of iterations for Mnu loop
2375 
2376  // approximately compute the number of points from the grid scanning
2377  // divide by abritry number to recover timing with still better results
2378  // m_NiterRandom=(NiterMET+1)*(NiterMET+1)*4*Niter*Niter/10;
2379 
2383 
2387 
2388  m_Mnu1Min = 0.;
2389  m_scanMnu1 = false;
2390  m_Mnu1 = m_Mnu1Min;
2391 
2392  // for markov chain use factor 2
2394 
2395  // NiterRandom set by user (default is -1). If negative, defines the default
2396  // here. no more automatic scaling for ll hl hh
2397  if (m_NiterRandom <= 0) {
2398  m_niterRandomLocal = 100000; // number of iterations for Markov for lh
2400  m_niterRandomLocal *= 2; // multiplied for ll , unchecked
2402  m_niterRandomLocal *= 5; // divided for hh ,checked
2403  } else {
2405  }
2406 
2407  if (preparedInput.m_type_visTau1 == 8) {
2408  // m_Mnu1Max=m_mTau-m_tauVec1M;
2411  m_scanMnu1 = true;
2412  }
2413 
2414  m_Mnu2Min = 0.;
2415  m_scanMnu2 = false;
2416  m_Mnu2 = m_Mnu2Min;
2417  if (preparedInput.m_type_visTau2 == 8) {
2418  // m_Mnu2Max=m_mTau-m_tauVec2M;
2421  m_scanMnu2 = true;
2422  }
2423 
2424  m_MEtLMin = -m_nsigma_METscan * METresX;
2425  m_MEtLMax = +m_nsigma_METscan * METresX;
2427 
2428  m_MEtPMin = -m_nsigma_METscan * METresY;
2429  m_MEtPMax = +m_nsigma_METscan * METresY;
2431 
2432  m_eTau1Min = -1;
2433  m_eTau1Max = -1;
2434  m_eTau2Min = -1;
2435  m_eTau2Max = -1;
2436 
2437  m_switch1 = true;
2438  m_switch2 = true;
2439 
2441  m_rmsStop = m_RMSStop;
2442 
2443  m_iter0 = -1;
2444  m_iterNuPV3 = 0;
2445  m_testptn1 = 0;
2446  m_testptn2 = 0;
2447  m_testdiscri1 = 0;
2448  m_testdiscri2 = 0;
2449  m_nosol1 = 0;
2450  m_nosol2 = 0;
2451  m_iterNsuc = 0;
2452  if (m_meanbinStop > 0) {
2453  m_meanbinToBeEvaluated = true;
2454  } else {
2455  m_meanbinToBeEvaluated = false;
2456  }
2457 
2461  m_markovNAccept = 0;
2462  m_markovNFullScan = 0;
2463  // set full parameter space scannning for the first steps, until a solution is
2464  // found
2465  m_fullParamSpaceScan = true;
2466  // size of step. Needs to be tune. Start with simple heuristic.
2467  if (m_proposalTryMEt < 0) {
2468  m_MEtProposal = m_MEtPRange / 30.;
2469  } else {
2471  }
2472  if (m_ProposalTryPhi < 0) {
2473  m_PhiProposal = 0.04;
2474  } else {
2476  }
2477  // FIXME if m_Mnu1Range !ne m_Mnu2Range same proposal will be done
2478  if (m_scanMnu1) {
2479  if (m_ProposalTryMnu < 0) {
2480  m_MnuProposal = m_Mnu1Range / 10.;
2481  } else {
2483  }
2484  }
2485  if (m_scanMnu2) {
2486  if (m_ProposalTryMnu < 0) {
2487  m_MnuProposal = m_Mnu2Range / 10.;
2488  } else {
2490  }
2491  }
2492 }
2493 
2494 // iterator. walk has internal counters, should only be used in a while loop
2495 // so far only implement grid strategy
2496 // act on MMC data member to be fast
2498  preparedInput.m_MEtX = -999.;
2499  preparedInput.m_MEtY = -999.;
2500 
2501  ++m_iter0;
2502 
2503  if (m_meanbinToBeEvaluated && m_iterNsuc == 500) {
2504  Info("DiTauMassTools", " in m_meanbinToBeEvaluated && m_iterNsuc==500 ");
2505  // for markov chain m_iterNsuc is the number of *accepted* points, so there
2506  // can be several iterations without any increment of m_iterNsuc. Hence need
2507  // to make sure meanbin is evaluated only once
2508  m_meanbinToBeEvaluated = false;
2509 
2510  // Meanbin stopping criterion
2511  std::vector<double> histInfo(HistInfo::MAXHISTINFO);
2512  // SLIDINGWINDOW strategy to avoid doing the parabola fit now given it will
2513  // not be use
2515  double meanbin = histInfo.at(HistInfo::MEANBIN);
2516  if (meanbin < 0) {
2517  m_nsucStop = -1; // no meaningful meanbin switch back to niter criterion
2518  } else {
2519  double stopdouble = 500 * std::pow((meanbin / m_meanbinStop), 2);
2520  int stopint = stopdouble;
2521  m_nsucStop = stopint;
2522  }
2523  if (m_nsucStop < 500)
2524  return false;
2525  }
2526  // should be outside m_meanbinStop test
2527  if (m_iterNsuc == m_nsucStop)
2528  return false; // Critere d'arret pour nombre de succes
2529 
2530  if (m_iter0 == m_niterRandomLocal)
2531  return false; // for now simple stopping criterion on number of iteration
2532 
2533  // floating stopping criterion, reduces run-time for lh, hh by a factor ~2 and ll by roughly
2534  // factor ~3 check if every scanned variable and resulting mass thermalised after N (default 10k) iterations
2535  // and then every M (default 1k) iterations do this by checking that the means of the split distributions is
2536  // comparable within X% (default 5%) of their sigma
2538  if (std::abs(m_fMEtP_split1->GetMean() - m_fMEtP_split2->GetMean()) <= m_fUseFloatStoppingComp * m_fMEtP_split1->GetRMS()) {
2539  if (std::abs(m_fMEtL_split1->GetMean() - m_fMEtL_split2->GetMean()) <=
2541  if (std::abs(m_fMnu1_split1->GetMean() - m_fMnu1_split2->GetMean()) <=
2543  if (std::abs(m_fMnu2_split1->GetMean() - m_fMnu2_split2->GetMean()) <=
2545  if (std::abs(m_fPhi1_split1->GetMean() - m_fPhi1_split2->GetMean()) <=
2547  if (std::abs(m_fPhi2_split1->GetMean() - m_fPhi2_split2->GetMean()) <=
2549  if (std::abs(m_fMmass_split1->GetMean() - m_fMmass_split2->GetMean()) <=
2551  return false;
2552  }
2553  }
2554  }
2555  }
2556  }
2557  }
2558  }
2559  }
2560 
2561  if (m_fullParamSpaceScan) {
2562  // as long as no solution found need to randomise on the full parameter
2563  // space
2564 
2565  // cut the corners in MissingET (not optimised at all)
2566  // not needed if distribution is already gaussian
2567  do {
2568  m_MEtP = m_MEtPMin + m_MEtPRange * m_randomGen.Rndm();
2569  m_MEtL = m_MEtLMin + m_MEtLRange * m_randomGen.Rndm();
2570  } while (!checkMEtInRange());
2571 
2572  if (m_scanMnu1) {
2573  m_Mnu1 = m_Mnu1Min + m_Mnu1Range * m_randomGen.Rndm();
2574  }
2575 
2576  if (m_scanMnu2) {
2577  m_Mnu2 = m_Mnu2Min + m_Mnu2Range * m_randomGen.Rndm();
2578  }
2579 
2580  m_Phi1 = m_Phi1Min + m_Phi1Range * m_randomGen.Rndm();
2581  m_Phi2 = m_Phi2Min + m_Phi2Range * m_randomGen.Rndm();
2582 
2583  return true;
2584  }
2585 
2586  // here the real markov chain takes place : "propose" the new point
2587  // note that if one parameter goes outside range, this should not be fixed
2588  // here but later in handleSolution, otherwise would cause a bias
2589 
2590  // m_MEtP0 etc... also store the position of the previous Markov Chain step,
2591  // which is needed by the algorithm
2592  m_MEtP0 = m_MEtP;
2593  m_MEtL0 = m_MEtL;
2594 
2596 
2598 
2599  if (m_scanMnu1) {
2600  m_Mnu10 = m_Mnu1;
2602  }
2603 
2604  if (m_scanMnu2) {
2605  m_Mnu20 = m_Mnu2;
2607  }
2608 
2609  m_Phi10 = m_Phi1;
2611 
2612  m_Phi20 = m_Phi2;
2613 
2615 
2616  return true;
2617 }
2618 
2619 // compute cached values (this value do not change within one call of MMC,
2620 // except for tau e scanning) return true if cache was already uptodatexs
2622 
2623  // copy tau 4 vect. If tau E scanning, these vectors will be modified
2626 
2627  const XYVector &metVec = preparedInput.m_MetVec;
2628 
2629  bool same = true;
2644 
2650 
2651  PtEtaPhiMVector Met4vec;
2652  Met4vec.SetPxPyPzE(preparedInput.m_MetVec.X(), preparedInput.m_MetVec.Y(), 0.0,
2653  preparedInput.m_MetVec.R());
2654  same = updateDouble((m_tauVec1 + m_tauVec2 + Met4vec).M(), m_Meff) && same;
2655 
2657  // note that if useHT met_vec is actually -HT
2658  same = updateDouble(metVec.X(), preparedInput.m_inputMEtX) && same;
2659  same = updateDouble(metVec.Y(), preparedInput.m_inputMEtY) && same;
2660  same = updateDouble(metVec.R(), preparedInput.m_inputMEtT) && same;
2661 
2662  return same;
2663 }
2664 
2665 // return true if all parameters are within their domain
2667 
2668  if (m_scanMnu1) {
2669  if (m_Mnu1 < m_Mnu1Min)
2670  return false;
2671  if (m_Mnu1 > m_Mnu1Max)
2672  return false;
2673  if (m_Mnu1 > m_mTau - m_tauVec1M)
2674  return false;
2675  }
2676 
2677  if (m_scanMnu2) {
2678  if (m_Mnu2 < m_Mnu2Min)
2679  return false;
2680  if (m_Mnu2 > m_Mnu2Max)
2681  return false;
2682  if (m_Mnu2 > m_mTau - m_tauVec2M)
2683  return false;
2684  }
2685 
2686  // FIXME note that since there is a coupling between Met and tau, should
2687  // rigorously test both together however since the 3 sigma range is just a
2688  // hack, it is probably OK
2689 
2690  if (m_Phi1 < m_Phi1Min)
2691  return false;
2692  if (m_Phi1 > m_Phi1Max)
2693  return false;
2694 
2695  if (m_Phi2 < m_Phi2Min)
2696  return false;
2697  if (m_Phi2 > m_Phi2Max)
2698  return false;
2699 
2700  if (!checkMEtInRange())
2701  return false;
2702 
2703  return true;
2704 }
2705 
2706 // return true if Met is within disk instead of withing square (cut the corners)
2708  // check MEt is in allowed range
2709  // range is 3sigma disk ("cutting the corners")
2713  return false;
2714  } else {
2715  return true;
2716  }
2717 }
2718 
2719 // ----- returns dTheta3D lower and upper boundaries:
2720 // limit_code=0: 99% lower limit
2721 // limit_code=1; 99% upper limit
2722 // limit_code=2; 95% upper limit
2723 double MissingMassCalculator::dTheta3DLimit(const int &tau_type, const int &limit_code,
2724  const double &P_tau) {
2725 
2726 #ifndef WITHDTHETA3DLIM
2727  // make the test ineffective if desired
2728  if (limit_code == 0)
2729  return 0.;
2730  if (limit_code == 1)
2731  return 10.;
2732  if (limit_code == 2)
2733  return 10.;
2734 #endif
2735 
2736  double limit = 1.0;
2737  if (limit_code == 0)
2738  limit = 0.0;
2739  double par[3] = {0.0, 0.0, 0.0};
2740  // ---- leptonic tau's
2741  if (tau_type == 8) {
2742  if (limit_code == 0) // lower 99% limit
2743  {
2744  par[0] = 0.3342;
2745  par[1] = -0.3376;
2746  par[2] = -0.001377;
2747  }
2748  if (limit_code == 1) // upper 99% limit
2749  {
2750  par[0] = 3.243;
2751  par[1] = -12.87;
2752  par[2] = 0.009656;
2753  }
2754  if (limit_code == 2) // upper 95% limit
2755  {
2756  par[0] = 2.927;
2757  par[1] = -7.911;
2758  par[2] = 0.007783;
2759  }
2760  }
2761  // ---- 1-prong tau's
2762  if (tau_type >= 0 && tau_type <= 2) {
2763  if (limit_code == 0) // lower 99% limit
2764  {
2765  par[0] = 0.2673;
2766  par[1] = -14.8;
2767  par[2] = -0.0004859;
2768  }
2769  if (limit_code == 1) // upper 99% limit
2770  {
2771  par[0] = 9.341;
2772  par[1] = -15.88;
2773  par[2] = 0.0333;
2774  }
2775  if (limit_code == 2) // upper 95% limit
2776  {
2777  par[0] = 6.535;
2778  par[1] = -8.649;
2779  par[2] = 0.00277;
2780  }
2781  }
2782  // ---- 3-prong tau's
2783  if (tau_type >= 3 && tau_type <= 5) {
2784  if (limit_code == 0) // lower 99% limit
2785  {
2786  par[0] = 0.2308;
2787  par[1] = -15.24;
2788  par[2] = -0.0009458;
2789  }
2790  if (limit_code == 1) // upper 99% limit
2791  {
2792  par[0] = 14.58;
2793  par[1] = -6.043;
2794  par[2] = -0.00928;
2795  }
2796  if (limit_code == 2) // upper 95% limit
2797  {
2798  par[0] = 8.233;
2799  par[1] = -0.3018;
2800  par[2] = -0.009399;
2801  }
2802  }
2803 
2804  if (std::abs(P_tau + par[1]) > 0.0)
2805  limit = par[0] / (P_tau + par[1]) + par[2];
2806  if (limit_code == 0) {
2807  if (limit < 0.0) {
2808  limit = 0.0;
2809  } else if (limit > 0.03) {
2810  limit = 0.03;
2811  }
2812  } else {
2813  if (limit < 0.0 || limit > 0.5 * TMath::Pi()) {
2814  limit = 0.5 * TMath::Pi();
2815  } else if (limit < 0.05 && limit > 0.0) {
2816  limit = 0.05; // parameterization only runs up to P~220 GeV in this regime
2817  // will set an upper bound of 0.05
2818  }
2819  }
2820 
2821  return limit;
2822 }
2823 
2824 // checks units of input variables, converts into [GeV] if needed, make all
2825 // possible corrections DR new : now a second structure preparedInput is derived
2826 // from the input one which only has direct user input
2828  const xAOD::IParticle *part2,
2829  const xAOD::MissingET *met,
2830  const int &njets) {
2831  int mmcType1 = mmcType(part1);
2832  if (mmcType1 < 0)
2833  return; // return CP::CorrectionCode::Error;
2834 
2835  int mmcType2 = mmcType(part2);
2836  if (mmcType2 < 0)
2837  return; // return CP::CorrectionCode::Error;
2838 
2839  preparedInput.SetLFVmode(-2); // initialise LFV mode value for this event with being *not* LFV
2840  // if(getLFVMode(part1, part2, mmcType1, mmcType2) ==
2841  // CP::CorrectionCode::Error) {
2843  int LFVMode = getLFVMode(part1, part2, mmcType1, mmcType2);
2844  if (LFVMode == -1) {
2845  return; // return CP::CorrectionCode::Error;
2846  } else if (LFVMode != -2) {
2847  preparedInput.SetLFVmode(LFVMode);
2848  }
2849  }
2850 
2851  // this will be in MeV but MMC allows MeV
2852  // assume the mass is correct as well
2853  PtEtaPhiMVector tlvTau1(part1->pt(), part1->eta(), part1->phi(), part1->m());
2854  PtEtaPhiMVector tlvTau2(part2->pt(), part2->eta(), part2->phi(), part2->m());
2855 
2856  // Convert to GeV. In principle, MMC should cope with MeV but should check
2857  // thoroughly
2858  PtEtaPhiMVector fixedtau1;
2859  fixedtau1.SetCoordinates(tlvTau1.Pt() / GEV, tlvTau1.Eta(), tlvTau1.Phi(), tlvTau1.M() / GEV);
2860  PtEtaPhiMVector fixedtau2;
2861  fixedtau2.SetCoordinates(tlvTau2.Pt() / GEV, tlvTau2.Eta(), tlvTau2.Phi(), tlvTau2.M() / GEV);
2862 
2863  preparedInput.SetVisTauType(0, mmcType1);
2864  preparedInput.SetVisTauType(1, mmcType2);
2865  preparedInput.SetVisTauVec(0, fixedtau1);
2866  preparedInput.SetVisTauVec(1, fixedtau2);
2867 
2868  if (mmcType1 == 8 && mmcType2 == 8) {
2870  } else if (mmcType1 >= 0 && mmcType1 <= 5 && mmcType2 >= 0 && mmcType2 <= 5) {
2872  } else {
2874  }
2876  Info("DiTauMassTools", "%s", ("running for tau types "+std::to_string(preparedInput.m_type_visTau1)+" "+std::to_string(preparedInput.m_type_visTau2)).c_str());
2877  XYVector met_vec(met->mpx() / GEV, met->mpy() / GEV);
2878  preparedInput.SetMetVec(met_vec);
2880  Info("DiTauMassTools", "%s", ("passing SumEt="+std::to_string(met->sumet() / GEV)).c_str());
2881  preparedInput.SetSumEt(met->sumet() / GEV);
2882  preparedInput.SetNjet25(njets);
2883 
2884  // check that the calibration set has been chosen explicitly, otherwise abort
2886  Error("DiTauMassTools", "MMCCalibrationSet has not been set !. Please use "
2887  "fMMC.SetCalibrationSet(MMCCalibrationSet::MMC2019) or fMMC.SetCalibrationSet(MMCCalibrationSet::MMC2024)"
2888  ". Abort now. ");
2889  std::abort();
2890  }
2891  //----------- Re-ordering input info, to make sure there is no dependence of
2892  // results on input order
2893  // this might be needed because a random scan is used
2894  // highest pT tau is always first
2895  preparedInput.m_InputReorder = 0; // set flag to 0 by default, i.e. no re-ordering
2897  preparedInput.m_type_visTau2 == 8) // if hadron-lepton, reorder to have lepton first
2898  {
2900  1; // re-order to be done, this flag is to be checked in DoOutputInfo()
2901  } else if (!((preparedInput.m_type_visTau2 >= 0 && preparedInput.m_type_visTau2 <= 5) &&
2902  preparedInput.m_type_visTau1 == 8)) // if not lep-had nor had lep, reorder if tau1 is
2903  // after tau2 clockwise
2904  {
2905  if (fixPhiRange(preparedInput.m_vistau1.Phi() - preparedInput.m_vistau2.Phi()) > 0) {
2906  preparedInput.m_InputReorder = 1; // re-order to be done, this flag is to be
2907  // checked in DoOutputInfo()
2908  }
2909  }
2910 
2911  if (preparedInput.m_InputReorder == 1) // copy and re-order
2912  {
2916  }
2917  //--------- re-ordering is done ---------------------------------------
2918 
2920  std::abs(Phi_mpi_pi(preparedInput.m_vistau1.Phi() - preparedInput.m_vistau2.Phi()));
2921 
2922  for (unsigned int i = 0; i < preparedInput.m_jet4vecs.size(); i++) {
2923  // correcting sumEt, give priority to SetMetScanParamsUE()
2924  if (preparedInput.m_METScanScheme == 0) {
2925  if ((preparedInput.m_METsigmaP < 0.1 || preparedInput.m_METsigmaL < 0.1) &&
2927  preparedInput.m_jet4vecs[i].Pt() > 20.0) {
2928  if (preparedInput.m_fUseVerbose == 1) {
2929  Info("DiTauMassTools", "correcting sumET");
2930  }
2932  }
2933  }
2934  }
2935 
2936  // give priority to SetVisTauType, only do this if type_visTau1 and
2937  // type_visTau2 are not set
2938  /*if(type_visTau1<0 && type_visTau2<0 && Nprong_tau1>-1 && Nprong_tau2>-1)
2939  {
2940  if(Nprong_tau1==0) type_visTau1 = 8; // leptonic tau
2941  else if( Nprong_tau1==1) type_visTau1 = 0; // set to 1p0n for now, may use
2942 different solution later like explicit integer for this case that pantau info is
2943 not available? else if( Nprong_tau1==3) type_visTau1 = 3; // set to 3p0n for
2944 now, see above if(Nprong_tau2==0) type_visTau2 = 8; // leptonic tau else if(
2945 Nprong_tau2==1) type_visTau2 = 0; // set to 1p0n for now, see above else if(
2946 Nprong_tau2==3) type_visTau2=3; // set to 3p0n for now, see above
2947  }
2948  */
2949  // checking input mass of hadronic tau-1
2950  // one prong
2951  // // checking input mass of hadronic tau-1
2952  // DRMERGE LFV addition
2955  preparedInput.m_vistau1.M() != 1.1) {
2957  preparedInput.m_vistau1.Phi(), 1.1);
2958  }
2960  preparedInput.m_vistau1.M() != 1.35) {
2962  preparedInput.m_vistau1.Phi(), 1.35);
2963  }
2964  // checking input mass of hadronic tau-2
2966  preparedInput.m_vistau2.M() != 1.1) {
2968  preparedInput.m_vistau2.Phi(), 1.1);
2969  }
2971  preparedInput.m_vistau2.M() != 1.35) {
2973  preparedInput.m_vistau2.Phi(), 1.35);
2974  }
2975  } else {
2976  // DRMERGE end LFV addition
2978  preparedInput.m_vistau1.M() != 0.8) {
2980  preparedInput.m_vistau1.Phi(), 0.8);
2981  }
2982  // 3 prong
2984  preparedInput.m_vistau1.M() != 1.2) {
2986  preparedInput.m_vistau1.Phi(), 1.2);
2987  }
2988  // checking input mass of hadronic tau-2
2989  // one prong
2991  preparedInput.m_vistau2.M() != 0.8) {
2993  preparedInput.m_vistau2.Phi(), 0.8);
2994  }
2995  // 3 prong
2997  preparedInput.m_vistau2.M() != 1.2) {
2999  preparedInput.m_vistau2.Phi(), 1.2);
3000  }
3001  } // DRDRMERGE LFV else closing
3002 
3003  // correcting sumEt for electron pt, give priority to SetMetScanParamsUE()
3004  // DR20150615 in tag 00-00-11 and before. The following was done before the
3005  // mass of the hadronic tau was set which mean that sumEt was wrongly
3006  // corrected for the hadronic tau pt if the hadronic tau mass was set to zero
3007  // Sasha 08/12/15: don't do electron Pt subtraction for high mass studies; in
3008  // the future, need to check if lepton Pt needs to be subtracted for both ele
3009  // and muon
3010  if (preparedInput.m_METsigmaP < 0.1 || preparedInput.m_METsigmaL < 0.1) {
3011 
3012  // T. Davidek: hack for lep-lep -- subtract lepton pT both for muon and
3013  // electron
3017  preparedInput.m_vistau1.M() < 0.12 && preparedInput.m_vistau2.M() < 0.12) { // lep-lep channel
3022  } else {
3023  // continue with the original code
3026  if (preparedInput.m_fUseVerbose == 1) {
3027  Info("DiTauMassTools", "Substracting pt1 from sumEt");
3028  }
3030  }
3033  if (preparedInput.m_fUseVerbose == 1) {
3034  Info("DiTauMassTools", "Substracting pt2 from sumEt");
3035  }
3037  }
3038  }
3039  }
3040 
3041  // controling TauProbability settings for UPGRADE studies
3044  if ((preparedInput.m_vistau1.M() < 0.12 && preparedInput.m_vistau2.M() > 0.12) ||
3045  (preparedInput.m_vistau2.M() < 0.12 && preparedInput.m_vistau1.M() > 0.12)) {
3046  Prob->SetUseTauProbability(true); // lep-had case
3047  }
3048  if (preparedInput.m_vistau1.M() > 0.12 && preparedInput.m_vistau2.M() > 0.12) {
3049  Prob->SetUseTauProbability(false); // had-had case
3050  }
3051  }
3052 
3053  // change Beam Energy for different running conditions
3055 
3056  //--------------------- pre-set defaults for Run-2. To disable pre-set
3057  // defaults set fUseDefaults=0
3058  if (preparedInput.m_fUseDefaults == 1) {
3060  SetNsigmaMETscan_ll(4.0);
3061  SetNsigmaMETscan_lh(4.0);
3062  SetNsigmaMETscan_hh(4.0);
3064  if ((preparedInput.m_vistau1.M() < 0.12 && preparedInput.m_vistau2.M() > 0.12) ||
3065  (preparedInput.m_vistau2.M() < 0.12 && preparedInput.m_vistau1.M() > 0.12))
3066  Prob->SetUseTauProbability(false); // lep-had
3068  Prob->SetUseTauProbability(true); // had-had
3069  Prob->SetUseMnuProbability(false);
3070  }
3071  }
3072 
3073  // compute HTOffset if relevant
3074  if (Prob->GetUseHT()) // use missing Ht for Njet25=0 events
3075  {
3076  // dPhi(l-t) dependence of misHt-trueMET
3077  double HtOffset = 0.;
3078  // proper for hh
3080  // hh
3081  double x = preparedInput.m_DelPhiTT;
3082  HtOffset = 87.5 - 27.0 * x;
3083  }
3084 
3085  preparedInput.m_HtOffset = HtOffset;
3086 
3087  // if use HT, replace MET with HT
3089  preparedInput.m_MHtSigma2; // sigma of 2nd Gaussian for missing Ht resolution
3091 
3092  PtEtaPhiMVector tauSum = preparedInput.m_vistau1 + preparedInput.m_vistau2;
3093  preparedInput.m_MetVec.SetXY(-tauSum.Px(), -tauSum.Py()); // WARNING this replace metvec by -mht
3094  }
3095 }
3096 
3099  if(!m_SaveLlhHisto) return;
3100 
3101  float hEmax = 3000.0; // maximum energy (GeV)
3102  int hNbins = 1500;
3103  m_fMEtP_all = std::make_shared<TH1F>("MEtP_h1", "M", hNbins, -100.0,
3104  100.); // all solutions
3105  m_fMEtL_all = std::make_shared<TH1F>("MEtL_h1", "M", hNbins, -100.0,
3106  100.); // all solutions
3107  m_fMnu1_all = std::make_shared<TH1F>("Mnu1_h1", "M", hNbins, 0.0,
3108  hEmax); // all solutions
3109  m_fMnu2_all = std::make_shared<TH1F>("Mnu2_h1", "M", hNbins, 0.0,
3110  hEmax); // all solutions
3111  m_fPhi1_all = std::make_shared<TH1F>("Phi1_h1", "M", hNbins, -10.0,
3112  10.); // all solutions
3113  m_fPhi2_all = std::make_shared<TH1F>("Phi2_h1", "M", hNbins, -10.0,
3114  10.); // all solutions
3115  m_fMfit_allGraph = std::make_shared<TGraph>(); // all solutions
3116 
3117  m_fMEtP_all->Sumw2();
3118  m_fMEtL_all->Sumw2();
3119  m_fMnu1_all->Sumw2();
3120  m_fMnu2_all->Sumw2();
3121  m_fPhi1_all->Sumw2();
3122  m_fPhi2_all->Sumw2();
3123 
3124  m_fMEtP_all->SetDirectory(0);
3125  m_fMEtL_all->SetDirectory(0);
3126  m_fMnu1_all->SetDirectory(0);
3127  m_fMnu2_all->SetDirectory(0);
3128  m_fPhi1_all->SetDirectory(0);
3129  m_fPhi2_all->SetDirectory(0);
3130 }
3131 
3134  if(!m_fUseFloatStopping) return;
3135 
3136  float hEmax = 3000.0; // maximum energy (GeV)
3137  int hNbins = 1500;
3138  m_fMmass_split1 = std::make_shared<TH1F>("mass_h1_1", "M", hNbins, 0.0, hEmax);
3139  m_fMEtP_split1 = std::make_shared<TH1F>("MEtP_h1_1", "M", hNbins, -100.0, 100.0);
3140  m_fMEtL_split1 = std::make_shared<TH1F>("MEtL_h1_1", "M", hNbins, -100.0, 100.0);
3141  m_fMnu1_split1 = std::make_shared<TH1F>("Mnu1_h1_1", "M", hNbins, 0.0, hEmax);
3142  m_fMnu2_split1 = std::make_shared<TH1F>("Mnu2_h1_1", "M", hNbins, 0.0, hEmax);
3143  m_fPhi1_split1 = std::make_shared<TH1F>("Phi1_h1_1", "M", hNbins, -10.0, 10.0);
3144  m_fPhi2_split1 = std::make_shared<TH1F>("Phi2_h1_1", "M", hNbins, -10.0, 10.0);
3145  m_fMmass_split2 = std::make_shared<TH1F>("mass_h1_2", "M", hNbins, 0.0, hEmax);
3146  m_fMEtP_split2 = std::make_shared<TH1F>("MEtP_h1_2", "M", hNbins, -100.0, 100.0);
3147  m_fMEtL_split2 = std::make_shared<TH1F>("MEtL_h1_2", "M", hNbins, -100.0, 100.0);
3148  m_fMnu1_split2 = std::make_shared<TH1F>("Mnu1_h1_2", "M", hNbins, 0.0, hEmax);
3149  m_fMnu2_split2 = std::make_shared<TH1F>("Mnu2_h1_2", "M", hNbins, 0.0, hEmax);
3150  m_fPhi1_split2 = std::make_shared<TH1F>("Phi1_h1_2", "M", hNbins, -10.0, 10.0);
3151  m_fPhi2_split2 = std::make_shared<TH1F>("Phi2_h1_2", "M", hNbins, -10.0, 10.0);
3152 
3153  m_fMmass_split1->Sumw2();
3154  m_fMEtP_split1->Sumw2();
3155  m_fMEtL_split1->Sumw2();
3156  m_fMnu1_split1->Sumw2();
3157  m_fMnu2_split1->Sumw2();
3158  m_fPhi1_split1->Sumw2();
3159  m_fPhi2_split1->Sumw2();
3160  m_fMmass_split2->Sumw2();
3161  m_fMEtP_split2->Sumw2();
3162  m_fMEtL_split2->Sumw2();
3163  m_fMnu1_split2->Sumw2();
3164  m_fMnu2_split2->Sumw2();
3165  m_fPhi1_split2->Sumw2();
3166  m_fPhi2_split2->Sumw2();
3167 
3168  m_fMmass_split1->SetDirectory(0);
3169  m_fMEtP_split1->SetDirectory(0);
3170  m_fMEtL_split1->SetDirectory(0);
3171  m_fMnu1_split1->SetDirectory(0);
3172  m_fMnu2_split1->SetDirectory(0);
3173  m_fPhi1_split1->SetDirectory(0);
3174  m_fPhi2_split1->SetDirectory(0);
3175  m_fMmass_split2->SetDirectory(0);
3176  m_fMEtP_split2->SetDirectory(0);
3177  m_fMEtL_split2->SetDirectory(0);
3178  m_fMnu1_split2->SetDirectory(0);
3179  m_fMnu2_split2->SetDirectory(0);
3180  m_fPhi1_split2->SetDirectory(0);
3181  m_fPhi2_split2->SetDirectory(0);
3182 }
DiTauMassTools::MMCCalibrationSet::MMC2024
@ MMC2024
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
DiTauMassTools::MissingMassCalculator::DitauStuff::Mditau_best
double Mditau_best
Definition: MissingMassCalculator.h:54
DiTauMassTools::MissingMassCalculator::m_Phi1
double m_Phi1
Definition: MissingMassCalculator.h:144
DiTauMassTools::MissingMassProb::SetUseDphiLL
void SetUseDphiLL(bool val)
Definition: MissingMassProb.h:57
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
DiTauMassTools::MissingMassCalculator::m_fFitting
TF1 * m_fFitting
Definition: MissingMassCalculator.h:234
DiTauMassTools::MissingMassCalculator::NuPsolutionLFV
int NuPsolutionLFV(const XYVector &met_vec, const PtEtaPhiMVector &tau, const double &m_nu, std::vector< PtEtaPhiMVector > &nu_vec)
Definition: MissingMassCalculator.cxx:755
DiTauMassTools::MissingMassCalculator::m_fMnu1_split2
std::shared_ptr< TH1F > m_fMnu1_split2
Definition: MissingMassCalculator.h:229
DiTauMassTools::MissingMassCalculator::m_testdiscri2
int m_testdiscri2
Definition: MissingMassCalculator.h:118
ParticleConstants::PDG2011::tauMassInMeV
constexpr double tauMassInMeV
the mass of the tau (in MeV)
Definition: ParticleConstants.h:32
DiTauMassTools::MissingMassCalculator::m_eTau2
double m_eTau2
Definition: MissingMassCalculator.h:145
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::MissingMassCalculator::m_tautau_tmp
PtEtaPhiMVector m_tautau_tmp
Definition: MissingMassCalculator.h:91
DiTauMassTools::MissingMassCalculator::m_switch2
bool m_switch2
Definition: MissingMassCalculator.h:123
DiTauMassTools::MissingMassCalculator::SaveLlhHisto
void SaveLlhHisto(const bool val)
Definition: MissingMassCalculator.cxx:3097
DiTauMassTools::MissingMassCalculator::m_nuvecsol1
std::vector< PtEtaPhiMVector > m_nuvecsol1
Definition: MissingMassCalculator.h:80
DiTauMassTools::MaxHistStrategy::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:30
DiTauMassTools::MissingMassCalculator::SpaceWalkerInit
void SpaceWalkerInit()
Definition: MissingMassCalculator.cxx:2317
DiTauMassTools::MissingMassCalculator::m_MEtPMax
double m_MEtPMax
Definition: MissingMassCalculator.h:149
DiTauMassTools::MissingMassCalculator::m_fMnu2_split2
std::shared_ptr< TH1F > m_fMnu2_split2
Definition: MissingMassCalculator.h:230
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:128
DiTauMassTools::MissingMassCalculator::precomputeCache
bool precomputeCache()
Definition: MissingMassCalculator.cxx:2621
DiTauMassTools::MissingMassCalculator::m_MEtL0
double m_MEtL0
Definition: MissingMassCalculator.h:147
DiTauMassTools::MissingMassCalculator::NuPsolutionV3
int NuPsolutionV3(const double &mNu1, const double &mNu2, const double &phi1, const double &phi2, int &nsol1, int &nsol2)
Definition: MissingMassCalculator.cxx:550
DiTauMassTools::MissingMassCalculator::DitauStuff::nutau1
PtEtaPhiMVector nutau1
Definition: MissingMassCalculator.h:56
DiTauMassTools::MMCFitMethod::MAX
@ MAX
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:46
DiTauMassTools::MMCCalibrationSet::MMC2015HIGHMASS
@ MMC2015HIGHMASS
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::MissingMassCalculator::m_ProposalTryPhi
double m_ProposalTryPhi
Definition: MissingMassCalculator.h:138
DiTauMassTools::MissingMassProb::GetUseTauProbability
bool GetUseTauProbability()
Definition: MissingMassProb.h:52
DiTauMassTools::MissingMassCalculator::m_nosol2
int m_nosol2
Definition: MissingMassCalculator.h:120
DiTauMassTools::TauTypes::hh
@ hh
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:53
DiTauMassTools::MissingMassCalculator::m_E2v1
double m_E2v1
Definition: MissingMassCalculator.h:192
DiTauMassTools::MissingMassCalculator::m_fMnu2_all
std::shared_ptr< TH1F > m_fMnu2_all
Definition: MissingMassCalculator.h:205
DiTauMassTools::MissingMassCalculator::m_fDitauStuffHisto
DitauStuff m_fDitauStuffHisto
Definition: MissingMassCalculator.h:250
DiTauMassTools::MissingMassCalculator::m_iter5
int m_iter5
Definition: MissingMassCalculator.h:101
DiTauMassTools::MissingMassCalculator::m_tauvecprob1
std::vector< double > m_tauvecprob1
Definition: MissingMassCalculator.h:85
DiTauMassTools::MissingMassCalculator::m_PrintmMaxError
double m_PrintmMaxError
Definition: MissingMassCalculator.h:133
DiTauMassTools::MissingMassCalculator::m_m2Nu1
double m_m2Nu1
Definition: MissingMassCalculator.h:188
DiTauMassTools::MMCCalibrationSet::name
const std::string name[MAXMMCCALIBRATIONSET]
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:41
DiTauMassTools::MissingMassInput::m_METsigmaP
double m_METsigmaP
Definition: MissingMassInput.h:62
DiTauMassTools::MissingMassCalculator::m_Ev2
double m_Ev2
Definition: MissingMassCalculator.h:194
DiTauMassTools::MissingMassCalculator::m_Phi2Max
double m_Phi2Max
Definition: MissingMassCalculator.h:149
DiTauMassTools::MissingMassCalculator::m_tauVec2Phi
double m_tauVec2Phi
Definition: MissingMassCalculator.h:181
DiTauMassTools::MissingMassCalculator::m_fMEtL_split1
std::shared_ptr< TH1F > m_fMEtL_split1
Definition: MissingMassCalculator.h:221
DiTauMassTools::MissingMassInput::m_METScanScheme
int m_METScanScheme
Definition: MissingMassInput.h:74
DiTauMassTools::MissingMassCalculator::m_fPhi2_split1
std::shared_ptr< TH1F > m_fPhi2_split1
Definition: MissingMassCalculator.h:225
DiTauMassTools::MissingMassCalculator::m_lfvLeplepRefit
bool m_lfvLeplepRefit
Definition: MissingMassCalculator.h:93
DiTauMassTools::MissingMassCalculator::refineSolutions
int refineSolutions(const double &M_nu1, const double &M_nu2, const int nsol1, const int nsol2, const double &Mvis, const double &Meff)
Definition: MissingMassCalculator.cxx:1821
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
DiTauMassTools::MMCFitMethod::MLM
@ MLM
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:46
DMTest::P
P_v1 P
Definition: P.h:23
DiTauMassTools::MissingMassCalculator::m_metCovPhiCos
double m_metCovPhiCos
Definition: MissingMassCalculator.h:153
DiTauMassTools::MissingMassCalculator::m_eventNumber
int m_eventNumber
Definition: MissingMassCalculator.h:109
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
met::DeltaR
@ DeltaR
Definition: METRecoCommon.h:11
DiTauMassTools::MissingMassInput::m_metVec
XYVector m_metVec
Definition: MissingMassInput.h:79
DiTauMassTools::MissingMassCalculator::DitauMassCalculatorV9lfv
int DitauMassCalculatorV9lfv(bool refit)
Definition: MissingMassCalculator.cxx:1014
DiTauMassTools::MissingMassProb::MET
void MET(MissingMassInput &preparedInput)
Definition: MissingMassProb.cxx:1166
DiTauMassTools::MissingMassCalculator::m_niter_fit2
int m_niter_fit2
Definition: MissingMassCalculator.h:253
DiTauMassTools::MissingMassCalculator::m_tauVec1
PtEtaPhiMVector m_tauVec1
Definition: MissingMassCalculator.h:180
DiTauMassTools::MissingMassCalculator::DitauStuff::vistau2
PtEtaPhiMVector vistau2
Definition: MissingMassCalculator.h:59
DiTauMassTools::MissingMassInput::m_inputMEtT
double m_inputMEtT
Definition: MissingMassInput.h:81
CSV_InDetExporter.new
new
Definition: CSV_InDetExporter.py:145
DiTauMassTools::MissingMassCalculator::m_nsolOld
int m_nsolOld
Definition: MissingMassCalculator.h:160
DiTauMassTools::MissingMassCalculator::m_Mvis
double m_Mvis
Definition: MissingMassCalculator.h:196
DiTauMassTools::MissingMassProb::SetUseMnuProbability
void SetUseMnuProbability(bool val)
Definition: MissingMassProb.h:54
DiTauMassTools::MissingMassCalculator::PrintOtherInput
void PrintOtherInput()
Definition: MissingMassCalculator.cxx:403
DiTauMassTools::MissingMassInput::m_METsigmaL
double m_METsigmaL
Definition: MissingMassInput.h:63
DiTauMassTools::MissingMassCalculator::m_PhiProposal
double m_PhiProposal
Definition: MissingMassCalculator.h:152
DiTauMassTools::MissingMassCalculator::m_fPYfit2
std::shared_ptr< TH1F > m_fPYfit2
Definition: MissingMassCalculator.h:215
DiTauMassTools::MissingMassCalculator::m_fPhi1_all
std::shared_ptr< TH1F > m_fPhi1_all
Definition: MissingMassCalculator.h:206
DiTauMassTools::MissingMassCalculator::DitauStuff::vistau1
PtEtaPhiMVector vistau1
Definition: MissingMassCalculator.h:58
DiTauMassTools::MissingMassCalculator::m_Phi2
double m_Phi2
Definition: MissingMassCalculator.h:144
DiTauMassTools::MissingMassOutput::m_FittedMass
double m_FittedMass[MMCFitMethod::MAX]
Definition: MissingMassOutput.h:55
DiTauMassTools::MissingMassCalculator::m_meanbinToBeEvaluated
bool m_meanbinToBeEvaluated
Definition: MissingMassCalculator.h:125
DiTauMassTools::MissingMassCalculator::m_tauVec2
PtEtaPhiMVector m_tauVec2
Definition: MissingMassCalculator.h:180
DiTauMassTools::MissingMassCalculator::m_fPXfit2
std::shared_ptr< TH1F > m_fPXfit2
Definition: MissingMassCalculator.h:214
DiTauMassTools::MissingMassCalculator::m_nsigma_METscan_lfv_ll
double m_nsigma_METscan_lfv_ll
Definition: MissingMassCalculator.h:98
DiTauMassTools::MissingMassInput::m_vistau1
PtEtaPhiMVector m_vistau1
Definition: MissingMassInput.h:54
DiTauMassTools::MissingMassCalculator::m_ET2v2
double m_ET2v2
Definition: MissingMassCalculator.h:191
DiTauMassTools::MissingMassOutput::m_FitSignificance
double m_FitSignificance[MMCFitMethod::MAX]
Definition: MissingMassOutput.h:54
DiTauMassTools::MissingMassCalculator::m_MEtLMin
double m_MEtLMin
Definition: MissingMassCalculator.h:148
DiTauMassTools::MissingMassCalculator::m_mmcCalibrationSet
MMCCalibrationSet::e m_mmcCalibrationSet
Definition: MissingMassCalculator.h:65
DiTauMassTools::MissingMassCalculator::m_niterRandomLocal
int m_niterRandomLocal
Definition: MissingMassCalculator.h:74
DiTauMassTools::MissingMassCalculator::m_nsigma_METscan_hh
double m_nsigma_METscan_hh
Definition: MissingMassCalculator.h:98
DiTauMassTools::MissingMassCalculator::m_iang1low
int m_iang1low
Definition: MissingMassCalculator.h:101
DiTauMassTools::MissingMassCalculator::m_MEtP0
double m_MEtP0
Definition: MissingMassCalculator.h:147
DiTauMassTools::MissingMassInput::SetLFVmode
void SetLFVmode(int val)
Definition: MissingMassInput.h:47
DiTauMassTools::MissingMassCalculator::m_fullParamSpaceScan
bool m_fullParamSpaceScan
Definition: MissingMassCalculator.h:158
DiTauMassTools::MissingMassOutput::m_FittedMassLowerError
double m_FittedMassLowerError[MMCFitMethod::MAX]
Definition: MissingMassOutput.h:57
DiTauMassTools::MissingMassCalculator::m_Mnu2Max
double m_Mnu2Max
Definition: MissingMassCalculator.h:149
DiTauMassTools::MissingMassCalculator::m_Phi20
double m_Phi20
Definition: MissingMassCalculator.h:147
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
DiTauMassTools::MissingMassCalculator::m_tauVec1M
double m_tauVec1M
Definition: MissingMassCalculator.h:182
DiTauMassTools::HistInfo::TANTHETA
@ TANTHETA
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:35
DiTauMassTools::MissingMassCalculator::m_cosPhi2
double m_cosPhi2
Definition: MissingMassCalculator.h:176
plotBeamSpotVxVal.sumw
int sumw
Definition: plotBeamSpotVxVal.py:235
DiTauMassTools::MissingMassCalculator::m_nsigma_METscan2
double m_nsigma_METscan2
Definition: MissingMassCalculator.h:97
DiTauMassTools::MissingMassOutput::m_hMfit_all
std::shared_ptr< TH1F > m_hMfit_all
Definition: MissingMassOutput.h:65
DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingComp
double m_fUseFloatStoppingComp
Definition: MissingMassCalculator.h:71
DiTauMassTools::MissingMassCalculator::m_MEtP
double m_MEtP
Definition: MissingMassCalculator.h:144
DiTauMassTools::MaxHistStrategy::MAXBINWINDOW
@ MAXBINWINDOW
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:30
DiTauMassTools::MissingMassCalculator::m_tauVec2Py
double m_tauVec2Py
Definition: MissingMassCalculator.h:184
DiTauMassTools::HistInfo::MEANBIN
@ MEANBIN
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:35
DiTauMassTools::MMCFitMethod::name
const std::string name[MAX]
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:47
xAOD::IParticle
Class providing the definition of the 4-vector interface.
Definition: Event/xAOD/xAODBase/xAODBase/IParticle.h:41
DiTauMassTools::MissingMassCalculator::m_iterNuPV3
int m_iterNuPV3
Definition: MissingMassCalculator.h:114
covarianceTool.prob
prob
Definition: covarianceTool.py:678
DiTauMassTools::MissingMassCalculator::m_nsucStop
int m_nsucStop
Definition: MissingMassCalculator.h:75
Phi_mpi_pi
__HOSTDEV__ double Phi_mpi_pi(double)
Definition: GeoRegion.cxx:7
python.TrigEgammaFastCaloHypoTool.same
def same(val, tool)
Definition: TrigEgammaFastCaloHypoTool.py:10
x
#define x
DiTauMassTools::MissingMassCalculator::m_fMEtP_split1
std::shared_ptr< TH1F > m_fMEtP_split1
Definition: MissingMassCalculator.h:220
DiTauMassTools::MissingMassCalculator::m_mTau2
double m_mTau2
Definition: MissingMassCalculator.h:143
DiTauMassTools::MissingMassCalculator::m_fMnu2_split1
std::shared_ptr< TH1F > m_fMnu2_split1
Definition: MissingMassCalculator.h:223
DiTauMassTools::MMCCalibrationSet::LFVMMC2012
@ LFVMMC2012
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::MissingMassCalculator::m_fMEtP_all
std::shared_ptr< TH1F > m_fMEtP_all
Definition: MissingMassCalculator.h:202
DiTauMassTools::MissingMassCalculator::maxFitting
Double_t maxFitting(Double_t *x, Double_t *par)
Definition: MissingMassCalculator.cxx:1499
DiTauMassTools::MissingMassCalculator::SetNsigmaMETscan_ll
void SetNsigmaMETscan_ll(const double val)
Definition: MissingMassCalculator.h:390
DiTauMassTools::MissingMassCalculator::m_fUseFloatStopping
bool m_fUseFloatStopping
Definition: MissingMassCalculator.h:68
binWidth
void binWidth(TH1 *h)
Definition: listroot.cxx:80
DiTauMassTools::MissingMassCalculator::m_Mnu2
double m_Mnu2
Definition: MissingMassCalculator.h:144
DiTauMassTools::MissingMassCalculator::m_nsigma_METscan_lh
double m_nsigma_METscan_lh
Definition: MissingMassCalculator.h:98
DiTauMassTools::MissingMassCalculator::m_walkWeight
double m_walkWeight
Definition: MissingMassCalculator.h:175
DiTauMassTools::MissingMassCalculator::m_MEtPMin
double m_MEtPMin
Definition: MissingMassCalculator.h:148
DiTauMassTools::MissingMassCalculator::m_ET2v1
double m_ET2v1
Definition: MissingMassCalculator.h:190
DiTauMassTools::MissingMassCalculator::m_Mnu20
double m_Mnu20
Definition: MissingMassCalculator.h:147
DiTauMassTools::MissingMassCalculator::m_debugThisIteration
bool m_debugThisIteration
Definition: MissingMassCalculator.h:93
DiTauMassTools::fastSinCos
void fastSinCos(const double &phi, double &sinPhi, double &cosPhi)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:42
DiTauMassTools::MissingMassOutput::m_hMfit_allNoWeight
std::shared_ptr< TH1F > m_hMfit_allNoWeight
Definition: MissingMassOutput.h:66
DiTauMassTools::MissingMassCalculator::m_fPZfit1
std::shared_ptr< TH1F > m_fPZfit1
Definition: MissingMassCalculator.h:213
DiTauMassTools::updateDouble
bool updateDouble(const double in, double &out)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:75
DiTauMassTools::MissingMassInput::m_METcovphi
double m_METcovphi
Definition: MissingMassInput.h:61
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:190
xAOD::EgammaParameters::deltaPhi1
@ deltaPhi1
difference between the cluster eta (1st sampling) and the eta of the track extrapolated to the 1st sa...
Definition: EgammaEnums.h:196
DiTauMassTools::MissingMassCalculator::m_markovNRejectNoSol
int m_markovNRejectNoSol
Definition: MissingMassCalculator.h:129
DiTauMassTools::MissingMassCalculator::m_fMEtL_split2
std::shared_ptr< TH1F > m_fMEtL_split2
Definition: MissingMassCalculator.h:228
DiTauMassTools::MissingMassInput::SetVisTauVec
void SetVisTauVec(int i, const PtEtaPhiMVector &vec)
Definition: MissingMassInput.cxx:123
DiTauMassTools::MissingMassInput::m_tauTypes
TauTypes::e m_tauTypes
Definition: MissingMassInput.h:85
DiTauMassTools::MissingMassCalculator::m_Mnu2Min
double m_Mnu2Min
Definition: MissingMassCalculator.h:148
DiTauMassTools::MissingMassCalculator::SpaceWalkerWalk
bool SpaceWalkerWalk()
Definition: MissingMassCalculator.cxx:2497
DiTauMassTools::MissingMassCalculator::m_Mnu1Range
double m_Mnu1Range
Definition: MissingMassCalculator.h:151
DiTauMassTools::HistInfo::MAXHISTINFO
@ MAXHISTINFO
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:35
DiTauMassTools::HistInfo::PROB
@ PROB
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:35
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::MissingMassCalculator::m_ProposalTryMnu
double m_ProposalTryMnu
Definition: MissingMassCalculator.h:139
part1
Definition: part1.py:1
DiTauMassTools
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:24
DiTauMassTools::MissingMassCalculator::m_nsigma_METscan_ll
double m_nsigma_METscan_ll
Definition: MissingMassCalculator.h:97
DiTauMassTools::MissingMassCalculator::m_fPXfit1
std::shared_ptr< TH1F > m_fPXfit1
Definition: MissingMassCalculator.h:211
DiTauMassTools::MissingMassCalculator::m_fMnu1_all
std::shared_ptr< TH1F > m_fMnu1_all
Definition: MissingMassCalculator.h:204
DiTauMassTools::MissingMassCalculator::m_fPhi1_split2
std::shared_ptr< TH1F > m_fPhi1_split2
Definition: MissingMassCalculator.h:231
DiTauMassTools::MissingMassCalculator::m_Mnu10
double m_Mnu10
Definition: MissingMassCalculator.h:147
DiTauMassTools::MaxDelPhi
double MaxDelPhi(int tau_type, double Pvis, double dRmax_tau)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:15
DiTauMassTools::MissingMassCalculator::m_switch1
bool m_switch1
Definition: MissingMassCalculator.h:122
met
Definition: IMETSignificance.h:24
DiTauMassTools::MissingMassCalculator::m_fMEtP_split2
std::shared_ptr< TH1F > m_fMEtP_split2
Definition: MissingMassCalculator.h:227
DiTauMassTools::MissingMassCalculator::m_Phi10
double m_Phi10
Definition: MissingMassCalculator.h:147
DiTauMassTools::MissingMassCalculator::m_nuvecsol2
std::vector< PtEtaPhiMVector > m_nuvecsol2
Definition: MissingMassCalculator.h:81
DiTauMassTools::MissingMassCalculator::m_nCallprobCalculatorV9fast
int m_nCallprobCalculatorV9fast
Definition: MissingMassCalculator.h:95
DiTauMassTools::MissingMassCalculator::m_Mnu1
double m_Mnu1
Definition: MissingMassCalculator.h:144
DiTauMassTools::MissingMassCalculator::preparedInput
MissingMassInput preparedInput
Definition: MissingMassCalculator.h:340
DiTauMassTools::MissingMassInput::m_htOffset
double m_htOffset
Definition: MissingMassInput.h:83
DiTauMassTools::MissingMassCalculator::m_Phi2Min
double m_Phi2Min
Definition: MissingMassCalculator.h:148
DiTauMassTools::MissingMassInput::m_SumEt
double m_SumEt
Definition: MissingMassInput.h:64
doubleTestComp.j1
j1
Definition: doubleTestComp.py:21
DiTauMassTools::MissingMassProb::MnuProbability
double MnuProbability(MissingMassInput &preparedInput, double mnu, double binsize)
Definition: MissingMassProb.cxx:945
DiTauMassTools::MissingMassCalculator::FinalizeSettings
void FinalizeSettings(const xAOD::IParticle *part1, const xAOD::IParticle *part2, const xAOD::MissingET *met, const int &njets)
Definition: MissingMassCalculator.cxx:2827
DiTauMassTools::MissingMassCalculator::m_fMnu1_split1
std::shared_ptr< TH1F > m_fMnu1_split1
Definition: MissingMassCalculator.h:222
TrigVtx::gamma
@ gamma
Definition: TrigParticleTable.h:27
DiTauMassTools::MissingMassInput::m_Nprong_tau2
int m_Nprong_tau2
Definition: MissingMassInput.h:59
DiTauMassTools::MissingMassCalculator::m_nu1FinalSolOldVec
std::vector< PtEtaPhiMVector > m_nu1FinalSolOldVec
Definition: MissingMassCalculator.h:163
DiTauMassTools::MissingMassCalculator::m_Ev1
double m_Ev1
Definition: MissingMassCalculator.h:195
ParseInputs.gDirectory
gDirectory
Definition: Final2012/ParseInputs.py:133
ParticleConstants.h
lumiFormat.i
int i
Definition: lumiFormat.py:85
DiTauMassTools::MaxHistStrategy::MAXBIN
@ MAXBIN
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:30
DiTauMassTools::MissingMassInput::m_MHtSigma2
double m_MHtSigma2
Definition: MissingMassInput.h:69
DiTauMassTools::MissingMassCalculator::m_MnuScanRange
double m_MnuScanRange
Definition: MissingMassCalculator.h:264
DiTauMassTools::MissingMassCalculator::m_tauVec1Pz
double m_tauVec1Pz
Definition: MissingMassCalculator.h:183
DiTauMassTools::MissingMassCalculator::m_ProposalTryEtau
double m_ProposalTryEtau
Definition: MissingMassCalculator.h:140
DiTauMassTools::MissingMassOutput::ClearOutput
void ClearOutput(bool fUseVerbose)
Definition: MissingMassOutput.cxx:24
DiTauMassTools::MissingMassCalculator::DitauStuff::nutau2
PtEtaPhiMVector nutau2
Definition: MissingMassCalculator.h:57
DiTauMassTools::MissingMassCalculator::~MissingMassCalculator
~MissingMassCalculator()
Definition: MissingMassCalculator.cxx:193
DiTauMassTools::MissingMassCalculator::m_PrintmInvWidth2Error
double m_PrintmInvWidth2Error
Definition: MissingMassCalculator.h:135
DiTauMassTools::MissingMassCalculator::m_niter_fit3
int m_niter_fit3
Definition: MissingMassCalculator.h:254
DiTauMassTools::MissingMassInput::m_MEtT
double m_MEtT
Definition: MissingMassInput.h:82
DiTauMassTools::MissingMassCalculator::m_fMfit_all
std::shared_ptr< TH1F > m_fMfit_all
Definition: MissingMassCalculator.h:201
DiTauMassTools::getLFVMode
int getLFVMode(const xAOD::IParticle *p1, const xAOD::IParticle *p2, int mmcType1, int mmcType2)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:83
DiTauMassTools::MissingMassCalculator::m_seed
int m_seed
Definition: MissingMassCalculator.h:110
DiTauMassTools::MissingMassInput::m_Nprong_tau1
int m_Nprong_tau1
Definition: MissingMassInput.h:58
DiTauMassTools::MissingMassCalculator::m_tauVec2P
double m_tauVec2P
Definition: MissingMassCalculator.h:185
DiTauMassTools::MissingMassCalculator::probCalculatorV9fast
int probCalculatorV9fast(const double &phi1, const double &phi2, const double &M_nu1, const double &M_nu2)
Definition: MissingMassCalculator.cxx:1795
DiTauMassTools::MissingMassCalculator::m_Phi1Max
double m_Phi1Max
Definition: MissingMassCalculator.h:149
DiTauMassTools::MMCCalibrationSet::MMC2019
@ MMC2019
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::MissingMassCalculator::m_tauVec1Phi
double m_tauVec1Phi
Definition: MissingMassCalculator.h:181
xAOD::MissingET_v1
Principal data object for Missing ET.
Definition: MissingET_v1.h:25
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::MissingMassCalculator::m_tauVec1E
double m_tauVec1E
Definition: MissingMassCalculator.h:186
DiTauMassTools::MissingMassProb::GetUseHT
bool GetUseHT()
Definition: MissingMassProb.h:49
DiTauMassTools::MaxHistStrategy::FIT
@ FIT
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:30
DiTauMassTools::MissingMassCalculator::m_randomGen
TRandom2 m_randomGen
Definition: MissingMassCalculator.h:63
DiTauMassTools::MissingMassCalculator::m_scanMnu2
bool m_scanMnu2
Definition: MissingMassCalculator.h:178
GEV
#define GEV
Definition: PrintPhotonSF.cxx:25
MissingMassCalculator.h
DiTauMassTools::MissingMassInput::SetNjet25
void SetNjet25(int val)
Definition: MissingMassInput.cxx:107
DiTauMassTools::MissingMassCalculator::SetNsigmaMETscan_hh
void SetNsigmaMETscan_hh(const double val)
Definition: MissingMassCalculator.h:392
DiTauMassTools::MissingMassCalculator::m_nsigma_METscan
double m_nsigma_METscan
Definition: MissingMassCalculator.h:97
DiTauMassTools::MMCCalibrationSet::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::MissingMassProb::GetUseMnuProbability
bool GetUseMnuProbability()
Definition: MissingMassProb.h:55
DiTauMassTools::MissingMassCalculator::DitauMassCalculatorV9walk
int DitauMassCalculatorV9walk()
Definition: MissingMassCalculator.cxx:792
DiTauMassTools::MissingMassOutput::m_NSolutions
int m_NSolutions
Definition: MissingMassOutput.h:69
DiTauMassTools::MaxHistStrategy::SLIDINGWINDOW
@ SLIDINGWINDOW
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:30
DiTauMassTools::MissingMassCalculator::m_tauVec2M
double m_tauVec2M
Definition: MissingMassCalculator.h:182
DiTauMassTools::MissingMassInput::m_dataType
int m_dataType
Definition: MissingMassInput.h:60
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
DiTauMassTools::mmcType
int mmcType(const xAOD::IParticle *part)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:136
DiTauMassTools::MissingMassCalculator::m_nsol
int m_nsol
Definition: MissingMassCalculator.h:166
DiTauMassTools::MissingMassInput::m_fUseVerbose
bool m_fUseVerbose
Definition: MissingMassInput.h:78
DiTauMassTools::MissingMassCalculator::m_reRunWithBestMET
bool m_reRunWithBestMET
Definition: MissingMassCalculator.h:197
DiTauMassTools::MissingMassCalculator::m_sinPhi2
double m_sinPhi2
Definition: MissingMassCalculator.h:176
DiTauMassTools::MissingMassCalculator::m_iter4
int m_iter4
Definition: MissingMassCalculator.h:101
DiTauMassTools::MissingMassCalculator::m_scanMnu1
bool m_scanMnu1
Definition: MissingMassCalculator.h:178
ReadTripsProbsFromCool.denominator
denominator
Definition: ReadTripsProbsFromCool.py:96
DiTauMassTools::MissingMassCalculator::m_beamEnergy
double m_beamEnergy
Definition: MissingMassCalculator.h:99
WriteCalibToCool.swap
swap
Definition: WriteCalibToCool.py:94
DiTauMassTools::MissingMassCalculator::checkAllParamInRange
bool checkAllParamInRange()
Definition: MissingMassCalculator.cxx:2666
DiTauMassTools::MissingMassOutput::m_objvec2
PtEtaPhiMVector m_objvec2[MMCFitMethod::MAX]
Definition: MissingMassOutput.h:61
DiTauMassTools::MissingMassCalculator::m_fMfit_allNoWeight
std::shared_ptr< TH1F > m_fMfit_allNoWeight
Definition: MissingMassCalculator.h:209
DiTauMassTools::MissingMassCalculator::m_mtautauSum
double m_mtautauSum
Definition: MissingMassCalculator.h:107
DiTauMassTools::MissingMassOutput::m_RMS2MPV
double m_RMS2MPV
Definition: MissingMassOutput.h:64
DiTauMassTools::MissingMassCalculator::m_markovNRejectMetropolis
int m_markovNRejectMetropolis
Definition: MissingMassCalculator.h:130
DiTauMassTools::MissingMassCalculator::m_eTau1Max
double m_eTau1Max
Definition: MissingMassCalculator.h:156
DiTauMassTools::MissingMassCalculator::checkMEtInRange
bool checkMEtInRange()
Definition: MissingMassCalculator.cxx:2707
DiTauMassTools::MissingMassCalculator::m_fPhi2_all
std::shared_ptr< TH1F > m_fPhi2_all
Definition: MissingMassCalculator.h:207
DiTauMassTools::MissingMassOutput::m_objvec1
PtEtaPhiMVector m_objvec1[MMCFitMethod::MAX]
Definition: MissingMassOutput.h:59
DiTauMassTools::MissingMassCalculator::m_fMmass_split1
std::shared_ptr< TH1F > m_fMmass_split1
Definition: MissingMassCalculator.h:219
DiTauMassTools::MissingMassCalculator::m_iterTheta3d
int m_iterTheta3d
Definition: MissingMassCalculator.h:102
DiTauMassTools::MissingMassInput::m_HtOffset
double m_HtOffset
Definition: MissingMassInput.h:71
part2
Definition: part2.py:1
DiTauMassTools::MissingMassCalculator::m_dRmax_tau
double m_dRmax_tau
Definition: MissingMassCalculator.h:262
DiTauMassTools::MissingMassCalculator::m_nsigma_METscan_lfv_lh
double m_nsigma_METscan_lfv_lh
Definition: MissingMassCalculator.h:99
DiTauMassTools::MissingMassProb::setParamRatio
void setParamRatio(int tau, int tautype)
Definition: MissingMassProb.cxx:190
DiTauMassTools::MissingMassCalculator::m_nsolmax
int m_nsolmax
Definition: MissingMassCalculator.h:73
DiTauMassTools::MissingMassCalculator::m_markovNFullScan
int m_markovNFullScan
Definition: MissingMassCalculator.h:128
DiTauMassTools::MissingMassCalculator::m_markovNAccept
int m_markovNAccept
Definition: MissingMassCalculator.h:131
DQPostProcessTest.outFile
outFile
Comment Out Those You do not wish to run.
Definition: DQPostProcessTest.py:36
DiTauMassTools::MissingMassCalculator::m_testdiscri1
int m_testdiscri1
Definition: MissingMassCalculator.h:117
DiTauMassTools::MissingMassCalculator::m_iter3
int m_iter3
Definition: MissingMassCalculator.h:101
DiTauMassTools::MissingMassCalculator::DitauStuff::RMSoverMPV
double RMSoverMPV
Definition: MissingMassCalculator.h:60
DiTauMassTools::MissingMassCalculator::m_Phi2Range
double m_Phi2Range
Definition: MissingMassCalculator.h:151
trigbs_pickEvents.num
num
Definition: trigbs_pickEvents.py:76
DiTauMassTools::MissingMassOutput::m_NSuccesses
int m_NSuccesses
Definition: MissingMassOutput.h:68
DiTauMassTools::MissingMassCalculator::m_Phi1Range
double m_Phi1Range
Definition: MissingMassCalculator.h:151
DiTauMassTools::MissingMassCalculator::RunMissingMassCalculator
int RunMissingMassCalculator(const xAOD::IParticle *part1, const xAOD::IParticle *part2, const xAOD::MissingET *met, const int &njets)
Definition: MissingMassCalculator.cxx:197
DiTauMassTools::MissingMassCalculator::m_tauVec2Px
double m_tauVec2Px
Definition: MissingMassCalculator.h:184
DiTauMassTools::HistInfo::RMS
@ RMS
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:35
DiTauMassTools::MissingMassInput::m_MetVec
XYVector m_MetVec
Definition: MissingMassInput.h:53
DiTauMassTools::MissingMassInput::m_type_visTau2
int m_type_visTau2
Definition: MissingMassInput.h:57
DiTauMassTools::MissingMassCalculator::m_nu1FinalSolVec
std::vector< PtEtaPhiMVector > m_nu1FinalSolVec
Definition: MissingMassCalculator.h:169
DiTauMassTools::MissingMassCalculator::m_fPZfit2
std::shared_ptr< TH1F > m_fPZfit2
Definition: MissingMassCalculator.h:216
DiTauMassTools::MissingMassCalculator::m_niter_fit1
int m_niter_fit1
Definition: MissingMassCalculator.h:252
DiTauMassTools::MissingMassCalculator::DitauStuff
Definition: MissingMassCalculator.h:53
DiTauMassTools::MissingMassProb::MetProbability
double MetProbability(MissingMassInput &preparedInput, const double &met1, const double &met2, const double &MetSigma1, const double &MetSigma2)
Definition: MissingMassProb.cxx:559
DiTauMassTools::MissingMassOutput::m_FittedMassUpperError
double m_FittedMassUpperError[MMCFitMethod::MAX]
Definition: MissingMassOutput.h:56
debug
const bool debug
Definition: MakeUncertaintyPlots.cxx:53
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
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
DiTauMassTools::MissingMassCalculator::m_eTau2Max
double m_eTau2Max
Definition: MissingMassCalculator.h:157
DiTauMassTools::MissingMassCalculator::m_proposalTryMEt
double m_proposalTryMEt
Definition: MissingMassCalculator.h:137
DiTauMassTools::MissingMassCalculator::m_probFinalSolOldVec
std::vector< double > m_probFinalSolOldVec
Definition: MissingMassCalculator.h:161
DiTauMassTools::MissingMassCalculator::m_tauvecsol1
std::vector< PtEtaPhiMVector > m_tauvecsol1
Definition: MissingMassCalculator.h:83
DiTauMassTools::MissingMassInput::PrintInputInfo
void PrintInputInfo()
Definition: MissingMassInput.cxx:58
DiTauMassTools::MissingMassCalculator::m_mtautauFinalSolOldVec
std::vector< double > m_mtautauFinalSolOldVec
Definition: MissingMassCalculator.h:162
DiTauMassTools::MissingMassProb::SetUseTauProbability
void SetUseTauProbability(bool val)
Definition: MissingMassProb.h:51
DiTauMassTools::MissingMassOutput::m_AveSolRMS
double m_AveSolRMS
Definition: MissingMassOutput.h:71
DiTauMassTools::MissingMassCalculator::m_tauVec1P
double m_tauVec1P
Definition: MissingMassCalculator.h:185
DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingCheckFreq
int m_fUseFloatStoppingCheckFreq
Definition: MissingMassCalculator.h:70
DiTauMassTools::MissingMassCalculator::m_eTau2Min
double m_eTau2Min
Definition: MissingMassCalculator.h:157
DiTauMassTools::MissingMassCalculator::m_fMfit_allGraph
std::shared_ptr< TGraph > m_fMfit_allGraph
Definition: MissingMassCalculator.h:208
xAOD::EgammaParameters::deltaPhi2
@ deltaPhi2
difference between the cluster phi (second sampling) and the phi of the track extrapolated to the sec...
Definition: EgammaEnums.h:204
DiTauMassTools::MissingMassCalculator::m_iter0
int m_iter0
Definition: MissingMassCalculator.h:113
DiTauMassTools::MissingMassCalculator::m_fPhi1_split1
std::shared_ptr< TH1F > m_fPhi1_split1
Definition: MissingMassCalculator.h:224
DiTauMassTools::MissingMassCalculator::TailCleanUp
int TailCleanUp(const PtEtaPhiMVector &vis1, const PtEtaPhiMVector &nu1, const PtEtaPhiMVector &vis2, const PtEtaPhiMVector &nu2, const double &mmc_mass, const double &vis_mass, const double &eff_mass, const double &dphiTT)
Definition: MissingMassCalculator.cxx:1992
DiTauMassTools::MissingMassCalculator::m_iang1high
int m_iang1high
Definition: MissingMassCalculator.h:101
DiTauMassTools::fixPhiRange
double fixPhiRange(const double &phi)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:21
DiTauMassTools::MissingMassCalculator::m_tauVec2Pz
double m_tauVec2Pz
Definition: MissingMassCalculator.h:184
DiTauMassTools::MissingMassCalculator::DitauStuff::Sign_best
double Sign_best
Definition: MissingMassCalculator.h:55
DiTauMassTools::MissingMassInput::m_fUseTailCleanup
bool m_fUseTailCleanup
Definition: MissingMassInput.h:77
DiTauMassTools::HistInfo::INTEGRAL
@ INTEGRAL
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:35
DiTauMassTools::MissingMassCalculator::m_MnuProposal
double m_MnuProposal
Definition: MissingMassCalculator.h:152
DiTauMassTools::MissingMassCalculator::m_fPhi2_split2
std::shared_ptr< TH1F > m_fPhi2_split2
Definition: MissingMassCalculator.h:232
DiTauMassTools::MissingMassOutput::m_nuvec1
PtEtaPhiMVector m_nuvec1[MMCFitMethod::MAX]
Definition: MissingMassOutput.h:58
DiTauMassTools::MissingMassCalculator::maxFromHist
double maxFromHist(TH1F *theHist, std::vector< double > &histInfo, const MaxHistStrategy::e maxHistStrategy=MaxHistStrategy::FIT, const int winHalfWidth=2, bool debug=false)
Definition: MissingMassCalculator.cxx:1519
DiTauMassTools::MissingMassCalculator::metvec_tmp
XYVector metvec_tmp
Definition: MissingMassCalculator.h:421
DiTauMassTools::MissingMassOutput::m_FittedMetVec
XYVector m_FittedMetVec[MMCFitMethod::MAX]
Definition: MissingMassOutput.h:63
DiTauMassTools::MissingMassCalculator::DoOutputInfo
void DoOutputInfo()
Definition: MissingMassCalculator.cxx:306
DiTauMassTools::MissingMassCalculator::m_RMSStop
int m_RMSStop
Definition: MissingMassCalculator.h:257
DiTauMassTools::MissingMassCalculator::m_NsucStop
int m_NsucStop
Definition: MissingMassCalculator.h:256
DiTauMassTools::MissingMassOutput::m_nuvec2
PtEtaPhiMVector m_nuvec2[MMCFitMethod::MAX]
Definition: MissingMassOutput.h:60
DiTauMassTools::MissingMassCalculator::m_eTau1Min
double m_eTau1Min
Definition: MissingMassCalculator.h:156
DiTauMassTools::MissingMassCalculator::m_testptn2
int m_testptn2
Definition: MissingMassCalculator.h:116
DiTauMassTools::HistInfo::DISCRI
@ DISCRI
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:35
DiTauMassTools::MissingMassCalculator::SetNsigmaMETscan_lh
void SetNsigmaMETscan_lh(const double val)
Definition: MissingMassCalculator.h:391
DiTauMassTools::MissingMassInput::m_MEtX
double m_MEtX
Definition: MissingMassInput.h:82
DiTauMassTools::MissingMassCalculator::MissingMassCalculator
MissingMassCalculator(MMCCalibrationSet::e aset, std::string paramFilePath)
Definition: MissingMassCalculator.cxx:48
DiTauMassTools::MissingMassInput::m_vistau2
PtEtaPhiMVector m_vistau2
Definition: MissingMassInput.h:55
DiTauMassTools::MissingMassInput::SetVisTauType
void SetVisTauType(int i, int tautype)
Definition: MissingMassInput.cxx:117
DiTauMassTools::MissingMassCalculator::m_iter1
int m_iter1
Definition: MissingMassCalculator.h:101
DiTauMassTools::MissingMassCalculator::m_fUseEfficiencyRecovery
bool m_fUseEfficiencyRecovery
Definition: MissingMassCalculator.h:67
DiTauMassTools::MissingMassCalculator::m_MEtL
double m_MEtL
Definition: MissingMassCalculator.h:144
DiTauMassTools::MissingMassCalculator::m_NiterRandom
int m_NiterRandom
Definition: MissingMassCalculator.h:255
DiTauMassTools::MissingMassProb
Definition: MissingMassProb.h:28
DiTauMassTools::MissingMassCalculator::m_MEtPRange
double m_MEtPRange
Definition: MissingMassCalculator.h:151
DiTauMassTools::MissingMassOutput::m_NTrials
int m_NTrials
Definition: MissingMassOutput.h:67
DiTauMassTools::MissingMassCalculator::m_eTau1
double m_eTau1
Definition: MissingMassCalculator.h:145
DiTauMassTools::MissingMassInput::m_MEtY
double m_MEtY
Definition: MissingMassInput.h:82
DiTauMassTools::MissingMassCalculator::m_dTheta3d_binMin
double m_dTheta3d_binMin
Definition: MissingMassCalculator.h:260
DiTauMassTools::MissingMassCalculator::m_MEtProposal
double m_MEtProposal
Definition: MissingMassCalculator.h:152
DiTauMassTools::MissingMassCalculator::m_MEtLRange
double m_MEtLRange
Definition: MissingMassCalculator.h:151
DiTauMassTools::MissingMassCalculator::m_iterNsuc
int m_iterNsuc
Definition: MissingMassCalculator.h:121
a
TList * a
Definition: liststreamerinfos.cxx:10
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
DiTauMassTools::MissingMassInput::ClearInput
void ClearInput()
Definition: MissingMassInput.cxx:24
DiTauMassTools::MissingMassProb::setParamNuMass
void setParamNuMass()
Definition: MissingMassProb.cxx:153
DiTauMassTools::MMCFitMethod::MAXW
@ MAXW
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:46
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
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
DiTauMassTools::MissingMassOutput::m_FitStatus
int m_FitStatus
Definition: MissingMassOutput.h:53
DiTauMassTools::MissingMassCalculator::m_eTau20
double m_eTau20
Definition: MissingMassCalculator.h:146
DiTauMassTools::MissingMassCalculator::m_fPYfit1
std::shared_ptr< TH1F > m_fPYfit1
Definition: MissingMassCalculator.h:212
DiTauMassTools::MissingMassCalculator::m_markovCountDuplicate
int m_markovCountDuplicate
Definition: MissingMassCalculator.h:127
DiTauMassTools::MMCFitMethod::MLNU3P
@ MLNU3P
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:46
DiTauMassTools::MMCCalibrationSet::MAXMMCCALIBRATIONSET
@ MAXMMCCALIBRATIONSET
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::HistInfo::FITLENGTH
@ FITLENGTH
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:35
DiTauMassTools::HistInfo::RMSVSDISCRI
@ RMSVSDISCRI
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:35
DiTauMassTools::MissingMassCalculator::m_Mnu1Min
double m_Mnu1Min
Definition: MissingMassCalculator.h:148
DiTauMassTools::MissingMassInput::m_jet4vecs
std::vector< PtEtaPhiMVector > m_jet4vecs
Definition: MissingMassInput.h:65
DiTauMassTools::MissingMassInput::m_type_visTau1
int m_type_visTau1
Definition: MissingMassInput.h:56
DiTauMassTools::MMCCalibrationSet::MMC2016MC15C
@ MMC2016MC15C
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::MissingMassCalculator::m_rmsStop
int m_rmsStop
Definition: MissingMassCalculator.h:76
DiTauMassTools::MissingMassInput::m_METresSyst
int m_METresSyst
Definition: MissingMassInput.h:75
DiTauMassTools::MissingMassCalculator::m_tauVec1Px
double m_tauVec1Px
Definition: MissingMassCalculator.h:183
DiTauMassTools::MissingMassCalculator::handleSolutions
void handleSolutions()
Definition: MissingMassCalculator.cxx:2048
DiTauMassTools::MissingMassCalculator::m_SaveLlhHisto
bool m_SaveLlhHisto
Definition: MissingMassCalculator.h:93
DiTauMassTools::MissingMassCalculator::m_m2Nu2
double m_m2Nu2
Definition: MissingMassCalculator.h:189
DiTauMassTools::MissingMassCalculator::SetUseFloatStopping
void SetUseFloatStopping(const bool val)
Definition: MissingMassCalculator.cxx:3132
DiTauMassTools::MissingMassCalculator::m_RndmSeedAltering
int m_RndmSeedAltering
Definition: MissingMassCalculator.h:258
DiTauMassTools::MissingMassInput::m_InputReorder
int m_InputReorder
Definition: MissingMassInput.h:73
DiTauMassTools::MissingMassOutput::m_totalvec
PtEtaPhiMVector m_totalvec[MMCFitMethod::MAX]
Definition: MissingMassOutput.h:62
DiTauMassTools::MissingMassCalculator::m_iter2
int m_iter2
Definition: MissingMassCalculator.h:101
DiTauMassTools::HistInfo::TANTHETAW
@ TANTHETAW
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:35
DiTauMassTools::MissingMassOutput::m_SumW
double m_SumW
Definition: MissingMassOutput.h:70
DiTauMassTools::MissingMassCalculator::m_nu2FinalSolVec
std::vector< PtEtaPhiMVector > m_nu2FinalSolVec
Definition: MissingMassCalculator.h:170
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
DiTauMassTools::MissingMassCalculator::m_probFinalSolVec
std::vector< double > m_probFinalSolVec
Definition: MissingMassCalculator.h:167
DiTauMassTools::MissingMassCalculator::m_Phi1Min
double m_Phi1Min
Definition: MissingMassCalculator.h:148
DiTauMassTools::MissingMassInput::m_beamEnergy
double m_beamEnergy
Definition: MissingMassInput.h:72
DiTauMassTools::MissingMassCalculator::m_meanbinStop
double m_meanbinStop
Definition: MissingMassCalculator.h:77
DiTauMassTools::MissingMassCalculator::m_mtautauFinalSolVec
std::vector< double > m_mtautauFinalSolVec
Definition: MissingMassCalculator.h:168
DiTauMassTools::MMCCalibrationSet::UPGRADE
@ UPGRADE
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::MissingMassCalculator::m_fMmass_split2
std::shared_ptr< TH1F > m_fMmass_split2
Definition: MissingMassCalculator.h:226
DiTauMassTools::MissingMassCalculator::m_PrintmMeanError
double m_PrintmMeanError
Definition: MissingMassCalculator.h:134
DiTauMassTools::MissingMassCalculator::m_totalProbSum
double m_totalProbSum
Definition: MissingMassCalculator.h:106
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingMinIter
int m_fUseFloatStoppingMinIter
Definition: MissingMassCalculator.h:69
DiTauMassTools::MissingMassCalculator::m_MEtLMax
double m_MEtLMax
Definition: MissingMassCalculator.h:149
DiTauMassTools::MissingMassCalculator::m_dTheta3d_binMax
double m_dTheta3d_binMax
Definition: MissingMassCalculator.h:261
DiTauMassTools::MissingMassInput::SetMetVec
void SetMetVec(const XYVector &vec)
Definition: MissingMassInput.cxx:130
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
DiTauMassTools::MissingMassInput::m_fUseDefaults
bool m_fUseDefaults
Definition: MissingMassInput.h:76
DiTauMassTools::MissingMassCalculator::m_tauvecprob2
std::vector< double > m_tauvecprob2
Definition: MissingMassCalculator.h:86
DiTauMassTools::MissingMassInput::SetSumEt
void SetSumEt(double sumEt)
Definition: MissingMassInput.cxx:112
DiTauMassTools::MissingMassCalculator::m_testptn1
int m_testptn1
Definition: MissingMassCalculator.h:115
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
DiTauMassTools::HistInfo::CHI2
@ CHI2
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:35
DiTauMassTools::MissingMassCalculator::m_cosPhi1
double m_cosPhi1
Definition: MissingMassCalculator.h:176
DiTauMassTools::MissingMassCalculator::m_nsolfinalmax
int m_nsolfinalmax
Definition: MissingMassCalculator.h:73
DiTauMassTools::MissingMassCalculator::PrintResults
void PrintResults()
Definition: MissingMassCalculator.cxx:445
python.compressB64.c
def c
Definition: compressB64.py:93
updateCoolNtuple.limit
int limit
Definition: updateCoolNtuple.py:44
DiTauMassTools::MissingMassCalculator::m_prob_tmp
double m_prob_tmp
Definition: MissingMassCalculator.h:104
DiTauMassTools::MissingMassInput::m_inputMEtY
double m_inputMEtY
Definition: MissingMassInput.h:81
DiTauMassTools::MissingMassCalculator::m_sinPhi1
double m_sinPhi1
Definition: MissingMassCalculator.h:176
DiTauMassTools::MissingMassCalculator::m_Mnu2Range
double m_Mnu2Range
Definition: MissingMassCalculator.h:151
DiTauMassTools::MissingMassCalculator::m_tauvecsol2
std::vector< PtEtaPhiMVector > m_tauvecsol2
Definition: MissingMassCalculator.h:84
DiTauMassTools::MissingMassCalculator::dTheta3DLimit
double dTheta3DLimit(const int &tau_type, const int &limit_code, const double &P_tau)
Definition: MissingMassCalculator.cxx:2723
DiTauMassTools::MissingMassCalculator::m_Mnu1Max
double m_Mnu1Max
Definition: MissingMassCalculator.h:149
DiTauMassTools::MissingMassCalculator::m_metCovPhiSin
double m_metCovPhiSin
Definition: MissingMassCalculator.h:153
doubleTestComp.j2
j2
Definition: doubleTestComp.py:22
DiTauMassTools::MissingMassInput::m_LFVmode
int m_LFVmode
Definition: MissingMassInput.h:84
DiTauMassTools::MissingMassCalculator::OutputInfo
MissingMassOutput OutputInfo
Definition: MissingMassCalculator.h:341
DiTauMassTools::MissingMassCalculator::m_Meff
double m_Meff
Definition: MissingMassCalculator.h:196
DiTauMassTools::MissingMassCalculator::m_nu2FinalSolOldVec
std::vector< PtEtaPhiMVector > m_nu2FinalSolOldVec
Definition: MissingMassCalculator.h:164
DiTauMassTools::MissingMassCalculator::Prob
MissingMassProb * Prob
Definition: MissingMassCalculator.h:342
DiTauMassTools::MissingMassCalculator::m_tauVec1Py
double m_tauVec1Py
Definition: MissingMassCalculator.h:183
DiTauMassTools::MissingMassCalculator::m_mTau
double m_mTau
Definition: MissingMassCalculator.h:143
DiTauMassTools::MissingMassCalculator::m_eTau10
double m_eTau10
Definition: MissingMassCalculator.h:146
DiTauMassTools::MissingMassCalculator::m_fDitauStuffFit
DitauStuff m_fDitauStuffFit
Definition: MissingMassCalculator.h:249
DiTauMassTools::MissingMassCalculator::m_tauVec2E
double m_tauVec2E
Definition: MissingMassCalculator.h:187
DiTauMassTools::TauTypes::ll
@ ll
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:53
DiTauMassTools::MissingMassCalculator::m_E2v2
double m_E2v2
Definition: MissingMassCalculator.h:193
DiTauMassTools::MissingMassCalculator::ClearDitauStuff
void ClearDitauStuff(DitauStuff &fStuff)
Definition: MissingMassCalculator.cxx:291
DiTauMassTools::MissingMassCalculator::m_nosol1
int m_nosol1
Definition: MissingMassCalculator.h:119
DiTauMassTools::MissingMassCalculator::m_fMEtL_all
std::shared_ptr< TH1F > m_fMEtL_all
Definition: MissingMassCalculator.h:203