26 #include <TFitResult.h>
27 #include <TFitResultPtr.h>
28 #include <TMatrixDSym.h>
31 #include "Math/VectorUtil.h"
34 using ROOT::Math::PtEtaPhiMVector;
35 using ROOT::Math::PxPyPzMVector;
36 using ROOT::Math::XYVector;
124 float hEmax = 3000.0;
134 m_fMfit_all = std::make_shared<TH1F>(
"MMC_h1",
"M", hNbins, 0.0,
141 std::make_shared<TH1F>(
"MMC_h1NoW",
"M no weight", hNbins, 0.0, hEmax);
143 m_fPXfit1 = std::make_shared<TH1F>(
"MMC_h2",
"Px1", 4 * hNbins, -hEmax,
145 m_fPYfit1 = std::make_shared<TH1F>(
"MMC_h3",
"Py1", 4 * hNbins, -hEmax,
147 m_fPZfit1 = std::make_shared<TH1F>(
"MMC_h4",
"Pz1", 4 * hNbins, -hEmax,
149 m_fPXfit2 = std::make_shared<TH1F>(
"MMC_h5",
"Px2", 4 * hNbins, -hEmax,
151 m_fPYfit2 = std::make_shared<TH1F>(
"MMC_h6",
"Py2", 4 * hNbins, -hEmax,
153 m_fPZfit2 = std::make_shared<TH1F>(
"MMC_h7",
"Pz2", 4 * hNbins, -hEmax,
170 m_fFitting->SetParNames(
"Max",
"Mean",
"InvWidth2");
195 Info(
"DiTauMassTools",
"------------- Raw Input for MissingMassCalculator --------------");
200 Info(
"DiTauMassTools",
"------------- Prepared Input for MissingMassCalculator--------------");
219 double dummy_METres =
240 Info(
"DiTauMassTools",
"Calling DitauMassCalculatorV9lfv");
246 TFile *
outFile = TFile::Open(
"MMC_likelihoods.root",
"UPDATE");
261 TH1D *nosol =
new TH1D(
"nosol",
"nosol", 7, 0, 7);
269 nosol->Write(nosol->GetName(), TObject::kOverwrite);
284 fStuff.
nutau1 = PtEtaPhiMVector(0., 0., 0., 0.);
285 fStuff.
nutau2 = PtEtaPhiMVector(0., 0., 0., 0.);
286 fStuff.
vistau1 = PtEtaPhiMVector(0., 0., 0., 0.);
287 fStuff.
vistau2 = PtEtaPhiMVector(0., 0., 0., 0.);
299 Info(
"DiTauMassTools",
"Retrieving output from fDitauStuffFit");
304 double q1 = (1. - 0.68) / 2.;
331 PtEtaPhiMVector tlvdummy(0., 0., 0., 0.);
332 XYVector metdummy(0., 0.);
373 PtEtaPhiMVector dummy_vec1(0.0, 0.0, 0.0, 0.0);
374 PtEtaPhiMVector dummy_vec2(0.0, 0.0, 0.0, 0.0);
375 for (
int i = 0;
i < 3;
i++) {
397 Info(
"DiTauMassTools",
398 ".........................Other input.....................................");
399 Info(
"DiTauMassTools",
"%s",
403 Info(
"DiTauMassTools",
"%s",
406 Info(
"DiTauMassTools",
"%s",
414 Info(
"DiTauMassTools",
415 "tau1 and tau2 were internally swapped (visible on prepared input printout)");
417 Info(
"DiTauMassTools",
"tau1 and tau2 were NOT internally swapped");
420 Info(
"DiTauMassTools",
"%s",
422 Info(
"DiTauMassTools",
"%s",
424 Info(
"DiTauMassTools",
"%s",
426 Info(
"DiTauMassTools",
"%s",
428 Info(
"DiTauMassTools",
"%s",
430 Info(
"DiTauMassTools",
"%s",
440 const PtEtaPhiMVector *origVisTau1 = 0;
441 const PtEtaPhiMVector *origVisTau2 = 0;
454 Info(
"DiTauMassTools",
455 "------------- Printing Final Results for MissingMassCalculator --------------");
456 Info(
"DiTauMassTools",
457 ".............................................................................");
461 Info(
"DiTauMassTools",
"%s",
464 Info(
"DiTauMassTools",
"%s",
470 Info(
"DiTauMassTools",
" no 4-momentum or MET from this method ");
475 Info(
"DiTauMassTools",
" fit failed ");
484 Info(
"DiTauMassTools",
"%s",
490 Info(
"DiTauMassTools",
"%s",
496 Info(
"DiTauMassTools",
"%s",
502 Info(
"DiTauMassTools",
"%s",
509 Info(
"DiTauMassTools",
"%s",
511 Info(
"DiTauMassTools",
"%s",
514 Info(
"DiTauMassTools",
"%s",
541 const double &phi1,
const double &phi2,
542 int &nsol1,
int &nsol2) {
563 int solution_code = 0;
581 if (pTmiss2 * dPhiSign < 0) {
583 return solution_code;
589 if (pTmiss1 * (-dPhiSign) < 0) {
591 return solution_code;
600 double m4noma1 = m2noma1 * m2noma1;
602 double p2v1proj =
std::pow(pv1proj, 2);
604 double pTmiss2CscDPhi = pTmiss2 / sinDPhi2;
605 double &pTn1 = pTmiss2CscDPhi;
606 double pT2miss2CscDPhi = pTmiss2CscDPhi * pTmiss2CscDPhi;
609 const double discri1 = m4noma1 + 4 * m2noma1 * pTmiss2CscDPhi * pv1proj -
610 4 * (
m_ET2v1 * (
m_m2Nu1 + pT2miss2CscDPhi) - (pT2miss2CscDPhi * p2v1proj));
615 return solution_code;
623 double m4noma2 = m2noma2 * m2noma2;
625 double p2v2proj =
std::pow(pv2proj, 2);
626 double sinDPhi1 = -sinDPhi2;
627 double pTmiss1CscDPhi = pTmiss1 / sinDPhi1;
628 double &pTn2 = pTmiss1CscDPhi;
629 double pT2miss1CscDPhi = pTmiss1CscDPhi * pTmiss1CscDPhi;
631 const double discri2 = m4noma2 + 4 * m2noma2 * pTmiss1CscDPhi * pv2proj -
632 4 * (
m_ET2v2 * (
m_m2Nu2 + pT2miss1CscDPhi) - (pT2miss1CscDPhi * p2v2proj));
637 return solution_code;
643 double sqdiscri1 = sqrt(discri1);
649 double pn1Z = first1 + second1;
651 if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z *
m_tauVec1Pz >
660 pn1Z = first1 - second1;
662 if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z *
m_tauVec1Pz >
674 return solution_code;
679 double sqdiscri2 = sqrt(discri2);
685 double pn2Z = first2 + second2;
687 if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z *
m_tauVec2Pz >
696 pn2Z = first2 - second2;
699 if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z *
m_tauVec2Pz >
710 return solution_code;
726 if (std::abs(pnux - pTmissx) > 0.001 || std::abs(pnuy - pTmissy) > 0.001) {
727 Info(
"DiTauMassTools",
"%s", (
"NuPsolutionV3 ERROR Pnux-Met.X or Pnuy-Met.Y > 0.001 : " +
732 if (std::abs(mtau1plus -
m_mTau) > 0.001 || std::abs(mtau1moins -
m_mTau) > 0.001 ||
733 std::abs(mtau2plus -
m_mTau) > 0.001 || std::abs(mtau2moins -
m_mTau) > 0.001) {
734 Info(
"DiTauMassTools",
"%s", (
"NuPsolutionV3 ERROR tau mass not recovered : " +
741 return solution_code;
746 const PtEtaPhiMVector &tau,
const double &l_nu,
747 std::vector<PtEtaPhiMVector> &nu_vec) {
748 int solution_code = 0;
751 PxPyPzMVector nu(met_vec.X(), met_vec.Y(), 0.0, l_nu);
752 PxPyPzMVector nu2(met_vec.X(), met_vec.Y(), 0.0, l_nu);
754 const double Mtau = 1.777;
756 double msq = (Mtau * Mtau - tau.M() * tau.M() - l_nu * l_nu) /
758 double gamma = nu.Px() * nu.Px() + nu.Py() * nu.Py();
759 double beta = tau.Px() * nu.Px() + tau.Py() * nu.Py() + msq;
760 double a = tau.E() * tau.E() - tau.Pz() * tau.Pz();
761 double b = -2 * tau.Pz() *
beta;
763 if ((
b *
b - 4 *
a *
c) < 0)
764 return solution_code;
767 double pvz1 = (-
b + sqrt(
b *
b - 4 *
a *
c)) / (2 *
a);
768 double pvz2 = (-
b - sqrt(
b *
b - 4 *
a *
c)) / (2 *
a);
770 nu.SetCoordinates(met_vec.X(), met_vec.Y(), pvz1, l_nu);
771 nu2.SetCoordinates(met_vec.X(), met_vec.Y(), pvz2, l_nu);
773 PtEtaPhiMVector return_nu(nu.Pt(), nu.Eta(), nu.Phi(), nu.M());
774 PtEtaPhiMVector return_nu2(nu2.Pt(), nu2.Eta(), nu2.Phi(), nu2.M());
775 nu_vec.push_back(return_nu);
776 nu_vec.push_back(return_nu2);
777 return solution_code;
835 XYVector deltamet_vec;
842 bool paramInsideRange =
false;
861 if (paramInsideRange)
885 if (nsuccesses > 0) {
891 double Px1, Py1, Pz1;
892 double Px2, Py2, Pz2;
893 if (nsuccesses > 0) {
914 if (prob_hist != 0.0)
938 PxPyPzMVector fulltau1, fulltau2;
939 fulltau1.SetCoordinates(Px1, Py1, Pz1, 1.777);
940 fulltau2.SetCoordinates(Px2, Py2, Pz2, 1.777);
960 Info(
"DiTauMassTools",
"Scanning ");
961 Info(
"DiTauMassTools",
" Markov ");
962 Info(
"DiTauMassTools",
"%s",
968 Info(
"DiTauMassTools",
"%s",
975 Info(
"DiTauMassTools",
"%s", (
"!!!----> Warning-3 in "
976 "MissingMassCalculator::DitauMassCalculatorV9Walk() : fit status=" +
979 Info(
"DiTauMassTools",
"%s",
"....... No solution is found. Printing input info .......");
997 Info(
"DiTauMassTools",
" ---------------------------------------------------------- ");
1019 const double Mtau = 1.777;
1039 double METresX_binSize = 2 * N_METsigma * METresX / NiterMET;
1040 double METresY_binSize = 2 * N_METsigma * METresY / NiterMET;
1044 std::vector<PtEtaPhiMVector> nu_vec;
1049 double metprob = 1.0;
1050 double sign_tmp = 0.0;
1051 double tauprob = 1.0;
1052 double totalProb = 0.0;
1056 double met_smear_x = 0.0;
1057 double met_smear_y = 0.0;
1058 double met_smearL = 0.0;
1059 double met_smearP = 0.0;
1061 double angle1 = 0.0;
1101 Info(
"DiTauMassTools",
"Running in dilepton mode");
1106 PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1107 PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1147 double Mlep = tau_tmp.M();
1154 double MnuProb = 1.0;
1156 for (
int i3 = 0; i3 < NiterMnu; i3++)
1158 M_nu = Mnu_binSize * i3;
1159 if (M_nu >= (Mtau - Mlep))
1165 for (
int i4 = 0; i4 < NiterMET + 1; i4++)
1167 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1168 for (
int i5 = 0; i5 < NiterMET + 1; i5++)
1170 met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1171 if (
pow(met_smearL / METresX, 2) +
pow(met_smearP / METresY, 2) >
pow(N_METsigma, 2))
1173 met_smear_x = met_smearL * met_coscovphi - met_smearP * met_sincovphi;
1174 met_smear_y = met_smearL * met_sincovphi + met_smearP * met_coscovphi;
1175 metvec_tmp.SetXY(input_metX + met_smear_x, input_metY + met_smear_y);
1193 for (
unsigned int j1 = 0;
j1 < nu_vec.size();
j1++) {
1196 const double tau1_tmpp = (tau_tmp + nu_vec[
j1]).
P();
1197 angle1 =
Angle(nu_vec[
j1], tau_tmp);
1207 double tauvecprob1j =
1209 if (tauvecprob1j == 0.)
1212 totalProb = tauvecprob1j * metprob * MnuProb * tauprob;
1229 m_fPXfit1->Fill((tau_tmp + nu_vec[
j1]).Px(), totalProb);
1230 m_fPYfit1->Fill((tau_tmp + nu_vec[
j1]).Py(), totalProb);
1231 m_fPZfit1->Fill((tau_tmp + nu_vec[
j1]).Pz(), totalProb);
1235 sign_tmp = -log10(totalProb);
1259 Info(
"DiTauMassTools",
"Running in lepton+tau mode");
1285 PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1286 PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1300 for (
int i4 = 0; i4 < NiterMET + 1; i4++)
1302 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1303 for (
int i5 = 0; i5 < NiterMET + 1; i5++)
1305 met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1306 if (
pow(met_smearL / METresX, 2) +
pow(met_smearP / METresY, 2) >
pow(N_METsigma, 2))
1310 metvec_tmp.SetXY(input_metX + met_smear_x, input_metY + met_smear_y);
1328 for (
unsigned int j1 = 0;
j1 < nu_vec.size();
j1++) {
1331 const double tau1_tmpp = (tau_tmp + nu_vec[
j1]).
P();
1332 angle1 =
Angle(nu_vec[
j1], tau_tmp);
1342 double tauvecprob1j =
1344 if (tauvecprob1j == 0.)
1347 totalProb = tauvecprob1j * metprob * tauprob;
1370 sign_tmp = -log10(totalProb);
1396 Info(
"DiTauMassTools",
"Running in an unknown mode?!?!");
1403 Info(
"DiTauMassTools",
"%s",
1425 if (prob_hist != 0.0)
1443 PxPyPzMVector nu1_tmp(0.0, 0.0, 0.0, 0.0);
1444 PxPyPzMVector nu2_tmp(0.0, 0.0, 0.0, 0.0);
1447 nu2_tmp.SetCoordinates(Px1, Py1, Pz1, 1.777);
1451 nu1_tmp.SetCoordinates(Px1, Py1, Pz1, 1.777);
1464 if (fit_code == 0) {
1466 "DiTauMassTools",
"%s",
1467 (
"!!!----> Warning-3 in MissingMassCalculator::DitauMassCalculatorV9lfv() : fit status=" +
1470 Info(
"DiTauMassTools",
"....... No solution is found. Printing input info .......");
1482 Info(
"DiTauMassTools",
" ---------------------------------------------------------- ");
1493 const double mM =
x[0];
1494 const double mMax =
par[0];
1495 const double mMean =
par[1];
1496 const double mInvWidth2 =
par[2];
1498 const double fitval = mMax * (1 - 4 * mInvWidth2 *
std::pow(mM - mMean, 2));
1511 const int winHalfWidth,
bool debug) {
1528 winHalfWidth == 0)) {
1532 int max_bin = theHist->GetMaximumBin();
1533 maxPos = theHist->GetBinCenter(max_bin);
1536 prob = theHist->GetBinContent(max_bin) /
double(theHist->GetEntries());
1543 int hNbins = theHist->GetNbinsX();
1548 int max_bin = theHist->GetMaximumBin();
1549 int iBinMin = max_bin - winHalfWidth;
1552 int iBinMax = max_bin + winHalfWidth;
1553 if (iBinMax > hNbins)
1554 iBinMax = hNbins - 1;
1557 for (
int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1558 const double weight = theHist->GetBinContent(iBin);
1560 sumx +=
weight * theHist->GetBinCenter(iBin);
1562 maxPos = sumx /
sumw;
1565 prob =
sumw / theHist->GetEntries();
1575 Error(
"DiTauMassTools",
"%s",
1576 (
"ERROR undefined maxHistStrategy:" +
std::to_string(maxHistStrategy)).c_str());
1582 int lastNonZeroBin = -1;
1583 int firstNonZeroBin = -1;
1584 double totalSumw = 0.;
1585 bool firstNullPart =
true;
1586 for (
int iBin = 0; iBin < hNbins; ++iBin) {
1587 const double weight = theHist->GetBinContent(iBin);
1590 lastNonZeroBin = iBin;
1591 if (firstNullPart) {
1592 firstNullPart =
false;
1593 firstNonZeroBin = iBin;
1600 firstNonZeroBin =
std::max(0, firstNonZeroBin - winHalfWidth - 1);
1601 lastNonZeroBin =
std::min(hNbins - 1, lastNonZeroBin + winHalfWidth + 1);
1610 const int nwidth = 2 * winHalfWidth + 1;
1613 for (
int ibin = 0; ibin < nwidth; ++ibin) {
1614 winsum += theHist->GetBinContent(ibin);
1616 double winmax = winsum;
1619 int iBinL = firstNonZeroBin;
1620 int iBinR = iBinL + 2 * winHalfWidth;
1621 bool goingUp =
true;
1626 const double deltawin = theHist->GetBinContent(iBinR) - theHist->GetBinContent(iBinL - 1);
1632 if (winsum > winmax) {
1635 max_bin = (iBinR + iBinL) / 2 - 1;
1646 }
while (iBinR < lastNonZeroBin);
1649 int iBinMin = max_bin - winHalfWidth;
1652 int iBinMax = max_bin + winHalfWidth;
1653 if (iBinMax >= hNbins)
1654 iBinMax = hNbins - 1;
1657 for (
int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1658 const double weight = theHist->GetBinContent(iBin);
1660 sumx +=
weight * theHist->GetBinCenter(iBin);
1663 double maxPosWin = -1.;
1666 maxPosWin = sumx /
sumw;
1673 const double h_rms = theHist->GetRMS(1);
1677 double numerator = 0;
1679 bool nullBin =
false;
1681 for (
int i = iBinMin;
i < iBinMax; ++
i) {
1682 double binError = theHist->GetBinError(
i);
1683 if (binError < 1
e-10) {
1686 double binErrorSquare =
std::pow(binError, 2);
1687 num = theHist->GetBinContent(
i) / (binErrorSquare);
1688 numerator = numerator +
num;
1691 if (numerator < 1
e-10 ||
denominator < 1
e-10 || nullBin ==
true) {
1707 const double binWidth = theHist->GetBinCenter(2) - theHist->GetBinCenter(1);
1708 double fitWidth = (winHalfWidth + 0.5) *
binWidth;
1721 TString fitOption =
debug ?
"QS" :
"QNS";
1727 TFitResultPtr fitRes =
1728 theHist->Fit(
m_fFitting, fitOption,
"", maxPos - fitWidth, maxPos + fitWidth);
1730 double maxPosFit = -1.;
1732 if (
int(fitRes) == 0) {
1735 const double mMax = fitRes->Parameter(0);
1736 const double mMean = fitRes->Parameter(1);
1737 const double mInvWidth2 = fitRes->Parameter(2);
1738 double mMaxError = fitRes->ParError(0);
1740 double mMeanError = fitRes->ParError(1);
1742 double mInvWidth2Error = fitRes->ParError(2);
1745 mInvWidth2Error = 0.;
1746 const double c = mMax * (1 - 4 * mMean * mMean * mInvWidth2);
1747 const double b = 8 * mMax * mMean * mInvWidth2;
1748 const double a = -4 * mMax * mInvWidth2;
1754 const double h_discri =
b *
b - 4 *
a *
c;
1756 const double sqrth_discri = sqrt(h_discri);
1757 const double h_fitLength = sqrth_discri /
a;
1764 maxPosFit = -
b / (2 *
a);
1768 if (maxPosFit >= 0. and std::abs(maxPosFit - maxPosWin) < 0.8 * fitWidth) {
1786 const double &M_nu1,
1787 const double &M_nu2) {
1793 const int solution =
NuPsolutionV3(M_nu1, M_nu2, phi1, phi2, nsol1, nsol2);
1812 const int nsol1,
const int nsol2,
1813 const double &Mvis,
const double &Meff)
1819 Error(
"DiTauMassTools",
"%s",
1823 Error(
"DiTauMassTools",
"%s",
1827 Error(
"DiTauMassTools",
"%s",
1831 Error(
"DiTauMassTools",
"%s",
1835 Error(
"DiTauMassTools",
"%s", (
"refineSolutions ERROR nsol1 " +
std::to_string(nsol1) +
1839 Error(
"DiTauMassTools",
"%s", (
"refineSolutions ERROR nsol1 " +
std::to_string(nsol2) +
1846 Prob->
apply(
preparedInput, -99, -99, PtEtaPhiMVector(0, 0, 0, 0), PtEtaPhiMVector(0, 0, 0, 0),
1847 PtEtaPhiMVector(0, 0, 0, 0), PtEtaPhiMVector(0, 0, 0, 0),
true,
false,
false);
1849 for (
int j1 = 0;
j1 < nsol1; ++
j1) {
1859 const int pickInt = std::abs(10000 *
m_Phi1);
1860 const int pickDigit = pickInt - 10 * (pickInt / 10);
1868 nuvec1_tmpj.SetCoordinates(nuvec1_tmpj.Pt(), nuvec1_tmpj.Eta(), nuvec1_tmpj.Phi(), M_nu1);
1869 tauvecsol1j.SetPxPyPzE(0., 0., 0., 0.);
1870 tauvecsol1j += nuvec1_tmpj;
1875 PtEtaPhiMVector(0, 0, 0, 0), nuvec1_tmpj,
1876 PtEtaPhiMVector(0, 0, 0, 0),
false,
true,
false);
1880 for (
int j2 = 0;
j2 < nsol2; ++
j2) {
1891 const int pickInt = std::abs(10000 *
m_Phi2);
1892 const int pickDigit = pickInt - 10 *
int(pickInt / 10);
1900 nuvec2_tmpj.SetCoordinates(nuvec2_tmpj.Pt(), nuvec2_tmpj.Eta(), nuvec2_tmpj.Phi(), M_nu2);
1901 tauvecsol2j.SetPxPyPzE(0., 0., 0., 0.);
1902 tauvecsol2j += nuvec2_tmpj;
1908 PtEtaPhiMVector(0, 0, 0, 0), nuvec2_tmpj,
false,
true,
false);
1912 if (tauvecprob1j == 0.)
1914 if (tauvecprob2j == 0.)
1917 double totalProb = 1.;
1930 (constProb * tauvecprob1j * tauvecprob2j *
1934 if (totalProb <= 0) {
1936 Warning(
"DiTauMassTools",
"%s",
1937 (
"null proba solution, rejected "+
std::to_string(totalProb)).c_str());
1944 Error(
"DiTauMassTools",
"%s",
1945 (
"refineSolutions ERROR nsol getting larger than nsolfinalmax!!! " +
1948 Error(
"DiTauMassTools",
"%s",
1964 nu1Final.SetPxPyPzE(nuvec1_tmpj.Px(), nuvec1_tmpj.Py(), nuvec1_tmpj.Pz(), nuvec1_tmpj.E());
1965 nu2Final.SetPxPyPzE(nuvec2_tmpj.Px(), nuvec2_tmpj.Py(), nuvec2_tmpj.Pz(), nuvec2_tmpj.E());
1973 if (ngoodsol1 == 0) {
1976 if (ngoodsol2 == 0) {
1983 const PtEtaPhiMVector &nu1,
1984 const PtEtaPhiMVector &vis2,
1985 const PtEtaPhiMVector &nu2,
const double &mmc_mass,
1986 const double &vis_mass,
const double &eff_mass,
1987 const double &dphiTT) {
1998 const double MrecoMvis = mmc_mass / vis_mass;
1999 if (MrecoMvis > 2.6)
2001 const double MrecoMeff = mmc_mass / eff_mass;
2002 if (MrecoMeff > 1.9)
2004 const double e1p1 = nu1.E() / vis1.P();
2005 const double e2p2 = nu2.E() / vis2.P();
2006 if ((e1p1 + e2p2) > 4.5)
2027 const double MrecoMvis = mmc_mass / vis_mass;
2028 const double MrecoMeff = mmc_mass / eff_mass;
2029 const double x = dphiTT > 1.5 ? dphiTT : 1.5;
2030 if ((MrecoMeff + MrecoMvis) > 5.908 - 1.881 *
x + 0.2995 *
x *
x)
2043 double totalProbSumSol = 0.;
2044 double totalProbSumSolOld = 0.;
2045 bool firstPointWithSol =
false;
2047 for (
int isol = 0; isol <
m_nsol; ++isol) {
2052 bool notSureToKeep =
true;
2056 notSureToKeep =
false;
2063 firstPointWithSol =
true;
2075 if (notSureToKeep) {
2078 for (
int isol = 0; isol <
m_nsolOld; ++isol) {
2084 if (!firstPointWithSol && totalProbSumSolOld <= 0.) {
2085 Error(
"DiTauMassTools",
"%s",
2086 (
" ERROR null old probability !!! " +
std::to_string(totalProbSumSolOld) +
" nsolOld " +
2090 }
else if (totalProbSumSol > totalProbSumSolOld) {
2095 }
else if (totalProbSumSol < totalProbSumSolOld * 1
E-6) {
2099 }
else if (
m_nsol <= 0) {
2106 reject = (uMC > totalProbSumSol / totalProbSumSolOld);
2131 bool fillSolution =
true;
2132 bool oldToBeUsed =
false;
2140 fillSolution =
false;
2154 if (!firstPointWithSol) {
2155 fillSolution =
true;
2158 fillSolution =
false;
2165 if (!fillSolution) {
2166 if (firstPointWithSol) {
2169 for (
int isol = 0; isol <
m_nsol; ++isol) {
2181 double solSum2 = 0.;
2183 for (
int isol = 0; isol <
m_nsol; ++isol) {
2187 const PtEtaPhiMVector *pnuvec1_tmpj;
2188 const PtEtaPhiMVector *pnuvec2_tmpj;
2201 const PtEtaPhiMVector &nuvec1_tmpj = *pnuvec1_tmpj;
2202 const PtEtaPhiMVector &nuvec2_tmpj = *pnuvec2_tmpj;
2205 solSum2 += mtautau * mtautau;
2227 if (mtautau != 0. &&
weight != 0.)
2292 for (
int isol = 0; isol <
m_nsol; ++isol) {
2494 Info(
"DiTauMassTools",
" in m_meanbinToBeEvaluated && m_iterNsuc==500 ");
2510 int stopint = stopdouble;
2641 PtEtaPhiMVector Met4vec;
2714 const double &P_tau) {
2716 #ifndef WITHDTHETA3DLIM
2718 if (limit_code == 0)
2720 if (limit_code == 1)
2722 if (limit_code == 2)
2727 if (limit_code == 0)
2729 double par[3] = {0.0, 0.0, 0.0};
2731 if (tau_type == 8) {
2732 if (limit_code == 0)
2738 if (limit_code == 1)
2744 if (limit_code == 2)
2752 if (tau_type >= 0 && tau_type <= 2) {
2753 if (limit_code == 0)
2757 par[2] = -0.0004859;
2759 if (limit_code == 1)
2765 if (limit_code == 2)
2773 if (tau_type >= 3 && tau_type <= 5) {
2774 if (limit_code == 0)
2778 par[2] = -0.0009458;
2780 if (limit_code == 1)
2786 if (limit_code == 2)
2794 if (std::abs(P_tau +
par[1]) > 0.0)
2796 if (limit_code == 0) {
2799 }
else if (
limit > 0.03) {
2803 if (limit < 0.0 || limit > 0.5 * TMath::Pi()) {
2804 limit = 0.5 * TMath::Pi();
2805 }
else if (limit < 0.05 && limit > 0.0) {
2821 const double GEV = 1000.;
2835 if (LFVMode == -1) {
2837 }
else if (LFVMode != -2) {
2849 PtEtaPhiMVector fixedtau1;
2850 fixedtau1.SetCoordinates(tlvTau1.Pt() /
GEV, tlvTau1.Eta(), tlvTau1.Phi(), tlvTau1.M() /
GEV);
2851 PtEtaPhiMVector fixedtau2;
2852 fixedtau2.SetCoordinates(tlvTau2.Pt() /
GEV, tlvTau2.Eta(), tlvTau2.Phi(), tlvTau2.M() /
GEV);
2859 if (mmcType1 == 8 && mmcType2 == 8) {
2861 }
else if (mmcType1 >= 0 && mmcType1 <= 5 && mmcType2 >= 0 && mmcType2 <= 5) {
2877 Error(
"DiTauMassTools",
"MMCCalibrationSet has not been set !. Please use "
2878 "fMMC.SetCalibrationSet(MMCCalibrationSet::MMC2019) or fMMC.SetCalibrationSet(MMCCalibrationSet::MMC2024)"
2920 Info(
"DiTauMassTools",
"correcting sumET");
3018 Info(
"DiTauMassTools",
"Substracting pt1 from sumEt");
3025 Info(
"DiTauMassTools",
"Substracting pt2 from sumEt");
3068 double HtOffset = 0.;
3073 HtOffset = 87.5 - 27.0 *
x;
3092 float hEmax = 3000.0;
3094 m_fMEtP_all = std::make_shared<TH1F>(
"MEtP_h1",
"M", hNbins, -100.0,
3096 m_fMEtL_all = std::make_shared<TH1F>(
"MEtL_h1",
"M", hNbins, -100.0,
3098 m_fMnu1_all = std::make_shared<TH1F>(
"Mnu1_h1",
"M", hNbins, 0.0,
3100 m_fMnu2_all = std::make_shared<TH1F>(
"Mnu2_h1",
"M", hNbins, 0.0,
3102 m_fPhi1_all = std::make_shared<TH1F>(
"Phi1_h1",
"M", hNbins, -10.0,
3104 m_fPhi2_all = std::make_shared<TH1F>(
"Phi2_h1",
"M", hNbins, -10.0,
3127 float hEmax = 3000.0;
3129 m_fMmass_split1 = std::make_shared<TH1F>(
"mass_h1_1",
"M", hNbins, 0.0, hEmax);
3130 m_fMEtP_split1 = std::make_shared<TH1F>(
"MEtP_h1_1",
"M", hNbins, -100.0, 100.0);
3131 m_fMEtL_split1 = std::make_shared<TH1F>(
"MEtL_h1_1",
"M", hNbins, -100.0, 100.0);
3132 m_fMnu1_split1 = std::make_shared<TH1F>(
"Mnu1_h1_1",
"M", hNbins, 0.0, hEmax);
3133 m_fMnu2_split1 = std::make_shared<TH1F>(
"Mnu2_h1_1",
"M", hNbins, 0.0, hEmax);
3134 m_fPhi1_split1 = std::make_shared<TH1F>(
"Phi1_h1_1",
"M", hNbins, -10.0, 10.0);
3135 m_fPhi2_split1 = std::make_shared<TH1F>(
"Phi2_h1_1",
"M", hNbins, -10.0, 10.0);
3136 m_fMmass_split2 = std::make_shared<TH1F>(
"mass_h1_2",
"M", hNbins, 0.0, hEmax);
3137 m_fMEtP_split2 = std::make_shared<TH1F>(
"MEtP_h1_2",
"M", hNbins, -100.0, 100.0);
3138 m_fMEtL_split2 = std::make_shared<TH1F>(
"MEtL_h1_2",
"M", hNbins, -100.0, 100.0);
3139 m_fMnu1_split2 = std::make_shared<TH1F>(
"Mnu1_h1_2",
"M", hNbins, 0.0, hEmax);
3140 m_fMnu2_split2 = std::make_shared<TH1F>(
"Mnu2_h1_2",
"M", hNbins, 0.0, hEmax);
3141 m_fPhi1_split2 = std::make_shared<TH1F>(
"Phi1_h1_2",
"M", hNbins, -10.0, 10.0);
3142 m_fPhi2_split2 = std::make_shared<TH1F>(
"Phi2_h1_2",
"M", hNbins, -10.0, 10.0);