26 #include <TFitResult.h>
27 #include <TFitResultPtr.h>
28 #include <TMatrixDSym.h>
31 #include "Math/VectorUtil.h"
36 constexpr
double GEV = 1000.0;
41 using ROOT::Math::PtEtaPhiMVector;
42 using ROOT::Math::PxPyPzMVector;
43 using ROOT::Math::XYVector;
134 float hEmax = 3000.0;
144 m_fMfit_all = std::make_shared<TH1F>(
"MMC_h1",
"M", hNbins, 0.0,
151 std::make_shared<TH1F>(
"MMC_h1NoW",
"M no weight", hNbins, 0.0, hEmax);
153 m_fPXfit1 = std::make_shared<TH1F>(
"MMC_h2",
"Px1", 4 * hNbins, -hEmax,
155 m_fPYfit1 = std::make_shared<TH1F>(
"MMC_h3",
"Py1", 4 * hNbins, -hEmax,
157 m_fPZfit1 = std::make_shared<TH1F>(
"MMC_h4",
"Pz1", 4 * hNbins, -hEmax,
159 m_fPXfit2 = std::make_shared<TH1F>(
"MMC_h5",
"Px2", 4 * hNbins, -hEmax,
161 m_fPYfit2 = std::make_shared<TH1F>(
"MMC_h6",
"Py2", 4 * hNbins, -hEmax,
163 m_fPZfit2 = std::make_shared<TH1F>(
"MMC_h7",
"Pz2", 4 * hNbins, -hEmax,
180 m_fFitting->SetParNames(
"Max",
"Mean",
"InvWidth2");
205 Info(
"DiTauMassTools",
"------------- Raw Input for MissingMassCalculator --------------");
210 Info(
"DiTauMassTools",
"------------- Prepared Input for MissingMassCalculator--------------");
229 double dummy_METres =
250 Info(
"DiTauMassTools",
"Calling DitauMassCalculatorV9lfv");
256 TFile *
outFile = TFile::Open(
"MMC_likelihoods.root",
"UPDATE");
271 TH1D *nosol =
new TH1D(
"nosol",
"nosol", 7, 0, 7);
279 nosol->Write(nosol->GetName(), TObject::kOverwrite);
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.);
309 Info(
"DiTauMassTools",
"Retrieving output from fDitauStuffFit");
314 double q1 = (1. - 0.68) / 2.;
341 PtEtaPhiMVector tlvdummy(0., 0., 0., 0.);
342 XYVector metdummy(0., 0.);
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++) {
407 Info(
"DiTauMassTools",
408 ".........................Other input.....................................");
409 Info(
"DiTauMassTools",
"%s",
413 Info(
"DiTauMassTools",
"%s",
416 Info(
"DiTauMassTools",
"%s",
424 Info(
"DiTauMassTools",
425 "tau1 and tau2 were internally swapped (visible on prepared input printout)");
427 Info(
"DiTauMassTools",
"tau1 and tau2 were NOT internally swapped");
430 Info(
"DiTauMassTools",
"%s",
432 Info(
"DiTauMassTools",
"%s",
434 Info(
"DiTauMassTools",
"%s",
436 Info(
"DiTauMassTools",
"%s",
438 Info(
"DiTauMassTools",
"%s",
440 Info(
"DiTauMassTools",
"%s",
450 const PtEtaPhiMVector *origVisTau1 = 0;
451 const PtEtaPhiMVector *origVisTau2 = 0;
464 Info(
"DiTauMassTools",
465 "------------- Printing Final Results for MissingMassCalculator --------------");
466 Info(
"DiTauMassTools",
467 ".............................................................................");
471 Info(
"DiTauMassTools",
"%s",
474 Info(
"DiTauMassTools",
"%s",
480 Info(
"DiTauMassTools",
" no 4-momentum or MET from this method ");
485 Info(
"DiTauMassTools",
" fit failed ");
494 Info(
"DiTauMassTools",
"%s",
500 Info(
"DiTauMassTools",
"%s",
506 Info(
"DiTauMassTools",
"%s",
512 Info(
"DiTauMassTools",
"%s",
519 Info(
"DiTauMassTools",
"%s",
521 Info(
"DiTauMassTools",
"%s",
524 Info(
"DiTauMassTools",
"%s",
551 const double &phi1,
const double &phi2,
552 int &nsol1,
int &nsol2) {
573 int solution_code = 0;
591 if (pTmiss2 * dPhiSign < 0) {
593 return solution_code;
599 if (pTmiss1 * (-dPhiSign) < 0) {
601 return solution_code;
610 double m4noma1 = m2noma1 * m2noma1;
612 double p2v1proj =
std::pow(pv1proj, 2);
614 double pTmiss2CscDPhi = pTmiss2 / sinDPhi2;
615 double &pTn1 = pTmiss2CscDPhi;
616 double pT2miss2CscDPhi = pTmiss2CscDPhi * pTmiss2CscDPhi;
619 const double discri1 = m4noma1 + 4 * m2noma1 * pTmiss2CscDPhi * pv1proj -
620 4 * (
m_ET2v1 * (
m_m2Nu1 + pT2miss2CscDPhi) - (pT2miss2CscDPhi * p2v1proj));
625 return solution_code;
633 double m4noma2 = m2noma2 * m2noma2;
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;
641 const double discri2 = m4noma2 + 4 * m2noma2 * pTmiss1CscDPhi * pv2proj -
642 4 * (
m_ET2v2 * (
m_m2Nu2 + pT2miss1CscDPhi) - (pT2miss1CscDPhi * p2v2proj));
647 return solution_code;
653 double sqdiscri1 = sqrt(discri1);
659 double pn1Z = first1 + second1;
661 if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z *
m_tauVec1Pz >
670 pn1Z = first1 - second1;
672 if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z *
m_tauVec1Pz >
684 return solution_code;
689 double sqdiscri2 = sqrt(discri2);
695 double pn2Z = first2 + second2;
697 if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z *
m_tauVec2Pz >
706 pn2Z = first2 - second2;
709 if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z *
m_tauVec2Pz >
720 return solution_code;
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 : " +
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 : " +
751 return solution_code;
756 const PtEtaPhiMVector &tau,
const double &l_nu,
757 std::vector<PtEtaPhiMVector> &nu_vec) {
758 int solution_code = 0;
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);
766 double msq = (Mtau * Mtau - tau.M() * tau.M() - l_nu * l_nu) /
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;
773 if ((
b *
b - 4 *
a *
c) < 0)
774 return solution_code;
777 double pvz1 = (-
b + sqrt(
b *
b - 4 *
a *
c)) / (2 *
a);
778 double pvz2 = (-
b - sqrt(
b *
b - 4 *
a *
c)) / (2 *
a);
780 nu.SetCoordinates(met_vec.X(), met_vec.Y(), pvz1, l_nu);
781 nu2.SetCoordinates(met_vec.X(), met_vec.Y(), pvz2, l_nu);
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;
845 XYVector deltamet_vec;
852 bool paramInsideRange =
false;
871 if (paramInsideRange)
895 if (nsuccesses > 0) {
901 double Px1, Py1, Pz1;
902 double Px2, Py2, Pz2;
903 if (nsuccesses > 0) {
924 if (prob_hist != 0.0)
948 PxPyPzMVector fulltau1, fulltau2;
970 Info(
"DiTauMassTools",
"Scanning ");
971 Info(
"DiTauMassTools",
" Markov ");
972 Info(
"DiTauMassTools",
"%s",
978 Info(
"DiTauMassTools",
"%s",
985 Info(
"DiTauMassTools",
"%s", (
"!!!----> Warning-3 in "
986 "MissingMassCalculator::DitauMassCalculatorV9Walk() : fit status=" +
989 Info(
"DiTauMassTools",
"%s",
"....... No solution is found. Printing input info .......");
1007 Info(
"DiTauMassTools",
" ---------------------------------------------------------- ");
1049 double METresX_binSize = 2 * N_METsigma * METresX / NiterMET;
1050 double METresY_binSize = 2 * N_METsigma * METresY / NiterMET;
1054 std::vector<PtEtaPhiMVector> nu_vec;
1059 double metprob = 1.0;
1060 double sign_tmp = 0.0;
1061 double tauprob = 1.0;
1062 double totalProb = 0.0;
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;
1071 double angle1 = 0.0;
1111 Info(
"DiTauMassTools",
"Running in dilepton mode");
1116 PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1117 PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1157 double Mlep = tau_tmp.M();
1164 double MnuProb = 1.0;
1166 for (
int i3 = 0; i3 < NiterMnu; i3++)
1168 M_nu = Mnu_binSize * i3;
1169 if (M_nu >= (Mtau - Mlep))
1175 for (
int i4 = 0; i4 < NiterMET + 1; i4++)
1177 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1178 for (
int i5 = 0; i5 < NiterMET + 1; i5++)
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))
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);
1203 for (
unsigned int j1 = 0;
j1 < nu_vec.size();
j1++) {
1206 const double tau1_tmpp = (tau_tmp + nu_vec[
j1]).
P();
1207 angle1 =
Angle(nu_vec[
j1], tau_tmp);
1217 double tauvecprob1j =
1219 if (tauvecprob1j == 0.)
1222 totalProb = tauvecprob1j * metprob * MnuProb * tauprob;
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);
1245 sign_tmp = -log10(totalProb);
1269 Info(
"DiTauMassTools",
"Running in lepton+tau mode");
1295 PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1296 PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1310 for (
int i4 = 0; i4 < NiterMET + 1; i4++)
1312 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1313 for (
int i5 = 0; i5 < NiterMET + 1; i5++)
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))
1320 metvec_tmp.SetXY(input_metX + met_smear_x, input_metY + met_smear_y);
1338 for (
unsigned int j1 = 0;
j1 < nu_vec.size();
j1++) {
1341 const double tau1_tmpp = (tau_tmp + nu_vec[
j1]).
P();
1342 angle1 =
Angle(nu_vec[
j1], tau_tmp);
1352 double tauvecprob1j =
1354 if (tauvecprob1j == 0.)
1357 totalProb = tauvecprob1j * metprob * tauprob;
1380 sign_tmp = -log10(totalProb);
1406 Info(
"DiTauMassTools",
"Running in an unknown mode?!?!");
1413 Info(
"DiTauMassTools",
"%s",
1435 if (prob_hist != 0.0)
1453 PxPyPzMVector nu1_tmp(0.0, 0.0, 0.0, 0.0);
1454 PxPyPzMVector nu2_tmp(0.0, 0.0, 0.0, 0.0);
1474 if (fit_code == 0) {
1476 "DiTauMassTools",
"%s",
1477 (
"!!!----> Warning-3 in MissingMassCalculator::DitauMassCalculatorV9lfv() : fit status=" +
1480 Info(
"DiTauMassTools",
"....... No solution is found. Printing input info .......");
1492 Info(
"DiTauMassTools",
" ---------------------------------------------------------- ");
1503 const double mM =
x[0];
1504 const double mMax =
par[0];
1505 const double mMean =
par[1];
1506 const double mInvWidth2 =
par[2];
1508 const double fitval = mMax * (1 - 4 * mInvWidth2 *
std::pow(mM - mMean, 2));
1521 const int winHalfWidth,
bool debug) {
1538 winHalfWidth == 0)) {
1542 int max_bin = theHist->GetMaximumBin();
1543 maxPos = theHist->GetBinCenter(max_bin);
1546 prob = theHist->GetBinContent(max_bin) /
double(theHist->GetEntries());
1553 int hNbins = theHist->GetNbinsX();
1558 int max_bin = theHist->GetMaximumBin();
1559 int iBinMin = max_bin - winHalfWidth;
1562 int iBinMax = max_bin + winHalfWidth;
1563 if (iBinMax > hNbins)
1564 iBinMax = hNbins - 1;
1567 for (
int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1568 const double weight = theHist->GetBinContent(iBin);
1570 sumx +=
weight * theHist->GetBinCenter(iBin);
1572 maxPos = sumx /
sumw;
1575 prob =
sumw / theHist->GetEntries();
1585 Error(
"DiTauMassTools",
"%s",
1586 (
"ERROR undefined maxHistStrategy:" +
std::to_string(maxHistStrategy)).c_str());
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);
1600 lastNonZeroBin = iBin;
1601 if (firstNullPart) {
1602 firstNullPart =
false;
1603 firstNonZeroBin = iBin;
1610 firstNonZeroBin =
std::max(0, firstNonZeroBin - winHalfWidth - 1);
1611 lastNonZeroBin =
std::min(hNbins - 1, lastNonZeroBin + winHalfWidth + 1);
1620 const int nwidth = 2 * winHalfWidth + 1;
1623 for (
int ibin = 0; ibin < nwidth; ++ibin) {
1624 winsum += theHist->GetBinContent(ibin);
1626 double winmax = winsum;
1629 int iBinL = firstNonZeroBin;
1630 int iBinR = iBinL + 2 * winHalfWidth;
1631 bool goingUp =
true;
1636 const double deltawin = theHist->GetBinContent(iBinR) - theHist->GetBinContent(iBinL - 1);
1642 if (winsum > winmax) {
1645 max_bin = (iBinR + iBinL) / 2 - 1;
1656 }
while (iBinR < lastNonZeroBin);
1659 int iBinMin = max_bin - winHalfWidth;
1662 int iBinMax = max_bin + winHalfWidth;
1663 if (iBinMax >= hNbins)
1664 iBinMax = hNbins - 1;
1667 for (
int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1668 const double weight = theHist->GetBinContent(iBin);
1670 sumx +=
weight * theHist->GetBinCenter(iBin);
1673 double maxPosWin = -1.;
1676 maxPosWin = sumx /
sumw;
1683 const double h_rms = theHist->GetRMS(1);
1687 double numerator = 0;
1689 bool nullBin =
false;
1691 for (
int i = iBinMin;
i < iBinMax; ++
i) {
1692 double binError = theHist->GetBinError(
i);
1693 if (binError < 1
e-10) {
1696 double binErrorSquare =
std::pow(binError, 2);
1697 num = theHist->GetBinContent(
i) / (binErrorSquare);
1698 numerator = numerator +
num;
1701 if (numerator < 1
e-10 ||
denominator < 1
e-10 || nullBin ==
true) {
1717 const double binWidth = theHist->GetBinCenter(2) - theHist->GetBinCenter(1);
1718 double fitWidth = (winHalfWidth + 0.5) *
binWidth;
1731 TString fitOption =
debug ?
"QS" :
"QNS";
1737 TFitResultPtr fitRes =
1738 theHist->Fit(
m_fFitting, fitOption,
"", maxPos - fitWidth, maxPos + fitWidth);
1740 double maxPosFit = -1.;
1742 if (
int(fitRes) == 0) {
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);
1750 double mMeanError = fitRes->ParError(1);
1752 double mInvWidth2Error = fitRes->ParError(2);
1755 mInvWidth2Error = 0.;
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;
1764 const double h_discri =
b *
b - 4 *
a *
c;
1766 const double sqrth_discri = sqrt(h_discri);
1767 const double h_fitLength = sqrth_discri /
a;
1774 maxPosFit = -
b / (2 *
a);
1778 if (maxPosFit >= 0. and std::abs(maxPosFit - maxPosWin) < 0.8 * fitWidth) {
1796 const double &M_nu1,
1797 const double &M_nu2) {
1803 const int solution =
NuPsolutionV3(M_nu1, M_nu2, phi1, phi2, nsol1, nsol2);
1822 const int nsol1,
const int nsol2,
1823 const double &Mvis,
const double &Meff)
1829 Error(
"DiTauMassTools",
"%s",
1833 Error(
"DiTauMassTools",
"%s",
1837 Error(
"DiTauMassTools",
"%s",
1841 Error(
"DiTauMassTools",
"%s",
1845 Error(
"DiTauMassTools",
"%s", (
"refineSolutions ERROR nsol1 " +
std::to_string(nsol1) +
1849 Error(
"DiTauMassTools",
"%s", (
"refineSolutions ERROR nsol1 " +
std::to_string(nsol2) +
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);
1859 for (
int j1 = 0;
j1 < nsol1; ++
j1) {
1869 const int pickInt = std::abs(10000 *
m_Phi1);
1870 const int pickDigit = pickInt - 10 * (pickInt / 10);
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;
1885 PtEtaPhiMVector(0, 0, 0, 0), nuvec1_tmpj,
1886 PtEtaPhiMVector(0, 0, 0, 0),
false,
true,
false);
1890 for (
int j2 = 0;
j2 < nsol2; ++
j2) {
1901 const int pickInt = std::abs(10000 *
m_Phi2);
1902 const int pickDigit = pickInt - 10 *
int(pickInt / 10);
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;
1918 PtEtaPhiMVector(0, 0, 0, 0), nuvec2_tmpj,
false,
true,
false);
1922 if (tauvecprob1j == 0.)
1924 if (tauvecprob2j == 0.)
1927 double totalProb = 1.;
1940 (constProb * tauvecprob1j * tauvecprob2j *
1944 if (totalProb <= 0) {
1946 Warning(
"DiTauMassTools",
"%s",
1947 (
"null proba solution, rejected "+
std::to_string(totalProb)).c_str());
1954 Error(
"DiTauMassTools",
"%s",
1955 (
"refineSolutions ERROR nsol getting larger than nsolfinalmax!!! " +
1958 Error(
"DiTauMassTools",
"%s",
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());
1983 if (ngoodsol1 == 0) {
1986 if (ngoodsol2 == 0) {
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) {
2008 const double MrecoMvis = mmc_mass / vis_mass;
2009 if (MrecoMvis > 2.6)
2011 const double MrecoMeff = mmc_mass / eff_mass;
2012 if (MrecoMeff > 1.9)
2014 const double e1p1 = nu1.E() / vis1.P();
2015 const double e2p2 = nu2.E() / vis2.P();
2016 if ((e1p1 + e2p2) > 4.5)
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)
2053 double totalProbSumSol = 0.;
2054 double totalProbSumSolOld = 0.;
2055 bool firstPointWithSol =
false;
2057 for (
int isol = 0; isol <
m_nsol; ++isol) {
2062 bool notSureToKeep =
true;
2066 notSureToKeep =
false;
2073 firstPointWithSol =
true;
2085 if (notSureToKeep) {
2088 for (
int isol = 0; isol <
m_nsolOld; ++isol) {
2094 if (!firstPointWithSol && totalProbSumSolOld <= 0.) {
2095 Error(
"DiTauMassTools",
"%s",
2096 (
" ERROR null old probability !!! " +
std::to_string(totalProbSumSolOld) +
" nsolOld " +
2100 }
else if (totalProbSumSol > totalProbSumSolOld) {
2105 }
else if (totalProbSumSol < totalProbSumSolOld * 1
E-6) {
2109 }
else if (
m_nsol <= 0) {
2116 reject = (uMC > totalProbSumSol / totalProbSumSolOld);
2141 bool fillSolution =
true;
2142 bool oldToBeUsed =
false;
2150 fillSolution =
false;
2164 if (!firstPointWithSol) {
2165 fillSolution =
true;
2168 fillSolution =
false;
2175 if (!fillSolution) {
2176 if (firstPointWithSol) {
2179 for (
int isol = 0; isol <
m_nsol; ++isol) {
2191 double solSum2 = 0.;
2193 for (
int isol = 0; isol <
m_nsol; ++isol) {
2197 const PtEtaPhiMVector *pnuvec1_tmpj;
2198 const PtEtaPhiMVector *pnuvec2_tmpj;
2211 const PtEtaPhiMVector &nuvec1_tmpj = *pnuvec1_tmpj;
2212 const PtEtaPhiMVector &nuvec2_tmpj = *pnuvec2_tmpj;
2215 solSum2 += mtautau * mtautau;
2237 if (mtautau != 0. &&
weight != 0.)
2302 for (
int isol = 0; isol <
m_nsol; ++isol) {
2504 Info(
"DiTauMassTools",
" in m_meanbinToBeEvaluated && m_iterNsuc==500 ");
2520 int stopint = stopdouble;
2651 PtEtaPhiMVector Met4vec;
2724 const double &P_tau) {
2726 #ifndef WITHDTHETA3DLIM
2728 if (limit_code == 0)
2730 if (limit_code == 1)
2732 if (limit_code == 2)
2737 if (limit_code == 0)
2739 double par[3] = {0.0, 0.0, 0.0};
2741 if (tau_type == 8) {
2742 if (limit_code == 0)
2748 if (limit_code == 1)
2754 if (limit_code == 2)
2762 if (tau_type >= 0 && tau_type <= 2) {
2763 if (limit_code == 0)
2767 par[2] = -0.0004859;
2769 if (limit_code == 1)
2775 if (limit_code == 2)
2783 if (tau_type >= 3 && tau_type <= 5) {
2784 if (limit_code == 0)
2788 par[2] = -0.0009458;
2790 if (limit_code == 1)
2796 if (limit_code == 2)
2804 if (std::abs(P_tau +
par[1]) > 0.0)
2806 if (limit_code == 0) {
2809 }
else if (
limit > 0.03) {
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) {
2844 if (LFVMode == -1) {
2846 }
else if (LFVMode != -2) {
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);
2868 if (mmcType1 == 8 && mmcType2 == 8) {
2870 }
else if (mmcType1 >= 0 && mmcType1 <= 5 && mmcType2 >= 0 && mmcType2 <= 5) {
2886 Error(
"DiTauMassTools",
"MMCCalibrationSet has not been set !. Please use "
2887 "fMMC.SetCalibrationSet(MMCCalibrationSet::MMC2019) or fMMC.SetCalibrationSet(MMCCalibrationSet::MMC2024)"
2929 Info(
"DiTauMassTools",
"correcting sumET");
3027 Info(
"DiTauMassTools",
"Substracting pt1 from sumEt");
3034 Info(
"DiTauMassTools",
"Substracting pt2 from sumEt");
3077 double HtOffset = 0.;
3082 HtOffset = 87.5 - 27.0 *
x;
3101 float hEmax = 3000.0;
3103 m_fMEtP_all = std::make_shared<TH1F>(
"MEtP_h1",
"M", hNbins, -100.0,
3105 m_fMEtL_all = std::make_shared<TH1F>(
"MEtL_h1",
"M", hNbins, -100.0,
3107 m_fMnu1_all = std::make_shared<TH1F>(
"Mnu1_h1",
"M", hNbins, 0.0,
3109 m_fMnu2_all = std::make_shared<TH1F>(
"Mnu2_h1",
"M", hNbins, 0.0,
3111 m_fPhi1_all = std::make_shared<TH1F>(
"Phi1_h1",
"M", hNbins, -10.0,
3113 m_fPhi2_all = std::make_shared<TH1F>(
"Phi2_h1",
"M", hNbins, -10.0,
3136 float hEmax = 3000.0;
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);