28 #include <TFitResult.h>
29 #include <TFitResultPtr.h>
30 #include <TMatrixDSym.h>
121 float hEmax = 3000.0;
131 m_fMfit_all = std::make_shared<TH1F>(
"MMC_h1",
"M", hNbins, 0.0,
133 #ifdef SAVELIKELIHOODHISTO
134 m_fMEtP_all = std::make_shared<TH1F>(
"MEtP_h1",
"M", hNbins, -100.0,
136 m_fMEtL_all = std::make_shared<TH1F>(
"MEtL_h1",
"M", hNbins, -100.0,
138 m_fMnu1_all = std::make_shared<TH1F>(
"Mnu1_h1",
"M", hNbins, 0.0,
140 m_fMnu2_all = std::make_shared<TH1F>(
"Mnu2_h1",
"M", hNbins, 0.0,
142 m_fPhi1_all = std::make_shared<TH1F>(
"Phi1_h1",
"M", hNbins, -10.0,
144 m_fPhi2_all = std::make_shared<TH1F>(
"Phi2_h1",
"M", hNbins, -10.0,
152 m_fMmass_split1 = std::make_shared<TH1F>(
"mass_h1_1",
"M", hNbins, 0.0, hEmax);
153 m_fMEtP_split1 = std::make_shared<TH1F>(
"MEtP_h1_1",
"M", hNbins, -100.0, 100.0);
154 m_fMEtL_split1 = std::make_shared<TH1F>(
"MEtL_h1_1",
"M", hNbins, -100.0, 100.0);
155 m_fMnu1_split1 = std::make_shared<TH1F>(
"Mnu1_h1_1",
"M", hNbins, 0.0, hEmax);
156 m_fMnu2_split1 = std::make_shared<TH1F>(
"Mnu2_h1_1",
"M", hNbins, 0.0, hEmax);
157 m_fPhi1_split1 = std::make_shared<TH1F>(
"Phi1_h1_1",
"M", hNbins, -10.0, 10.0);
158 m_fPhi2_split1 = std::make_shared<TH1F>(
"Phi2_h1_1",
"M", hNbins, -10.0, 10.0);
159 m_fMmass_split2 = std::make_shared<TH1F>(
"mass_h1_2",
"M", hNbins, 0.0, hEmax);
160 m_fMEtP_split2 = std::make_shared<TH1F>(
"MEtP_h1_2",
"M", hNbins, -100.0, 100.0);
161 m_fMEtL_split2 = std::make_shared<TH1F>(
"MEtL_h1_2",
"M", hNbins, -100.0, 100.0);
162 m_fMnu1_split2 = std::make_shared<TH1F>(
"Mnu1_h1_2",
"M", hNbins, 0.0, hEmax);
163 m_fMnu2_split2 = std::make_shared<TH1F>(
"Mnu2_h1_2",
"M", hNbins, 0.0, hEmax);
164 m_fPhi1_split2 = std::make_shared<TH1F>(
"Phi1_h1_2",
"M", hNbins, -10.0, 10.0);
165 m_fPhi2_split2 = std::make_shared<TH1F>(
"Phi2_h1_2",
"M", hNbins, -10.0, 10.0);
167 #ifdef SAVELIKELIHOODHISTO
195 std::make_shared<TH1F>(
"MMC_h1NoW",
"M no weight", hNbins, 0.0, hEmax);
197 m_fPXfit1 = std::make_shared<TH1F>(
"MMC_h2",
"Px1", 4 * hNbins, -hEmax,
199 m_fPYfit1 = std::make_shared<TH1F>(
"MMC_h3",
"Py1", 4 * hNbins, -hEmax,
201 m_fPZfit1 = std::make_shared<TH1F>(
"MMC_h4",
"Pz1", 4 * hNbins, -hEmax,
203 m_fPXfit2 = std::make_shared<TH1F>(
"MMC_h5",
"Px2", 4 * hNbins, -hEmax,
205 m_fPYfit2 = std::make_shared<TH1F>(
"MMC_h6",
"Py2", 4 * hNbins, -hEmax,
207 m_fPZfit2 = std::make_shared<TH1F>(
"MMC_h7",
"Pz2", 4 * hNbins, -hEmax,
211 #ifdef SAVELIKELIHOODHISTO
249 m_fFitting->SetParNames(
"Max",
"Mean",
"InvWidth2");
274 Info(
"DiTauMassTools",
"------------- Raw Input for MissingMassCalculator --------------");
279 Info(
"DiTauMassTools",
"------------- Prepared Input for MissingMassCalculator--------------");
298 double dummy_METres =
319 Info(
"DiTauMassTools",
"Calling DitauMassCalculatorV9lfv");
323 #ifdef SAVELIKELIHOODHISTO
324 TFile *
outFile = TFile::Open(
"MMC_likelihoods.root",
"UPDATE");
339 TH1D *nosol =
new TH1D(
"nosol",
"nosol", 7, 0, 7);
347 nosol->Write(nosol->GetName(), TObject::kOverwrite);
361 fStuff.
nutau1 = TLorentzVector(0., 0., 0., 0.);
362 fStuff.
nutau2 = TLorentzVector(0., 0., 0., 0.);
363 fStuff.
vistau1 = TLorentzVector(0., 0., 0., 0.);
364 fStuff.
vistau2 = TLorentzVector(0., 0., 0., 0.);
376 Info(
"DiTauMassTools",
"Retrieving output from fDitauStuffFit");
381 double q1 = (1. - 0.68) / 2.;
414 TLorentzVector tlvdummy(0., 0., 0., 0.);
415 TVector2 metdummy(0., 0.);
456 TLorentzVector dummy_vec1(0.0, 0.0, 0.0, 0.0);
457 TLorentzVector dummy_vec2(0.0, 0.0, 0.0, 0.0);
458 for (
int i = 0;
i < 3;
i++) {
480 Info(
"DiTauMassTools",
481 ".........................Other input.....................................");
482 Info(
"DiTauMassTools",
"%s",
486 Info(
"DiTauMassTools",
"%s",
489 Info(
"DiTauMassTools",
"%s",
497 Info(
"DiTauMassTools",
498 "tau1 and tau2 were internally swapped (visible on prepared input printout)");
500 Info(
"DiTauMassTools",
"tau1 and tau2 were NOT internally swapped");
503 Info(
"DiTauMassTools",
"%s",
505 Info(
"DiTauMassTools",
"%s",
507 Info(
"DiTauMassTools",
"%s",
509 Info(
"DiTauMassTools",
"%s",
511 Info(
"DiTauMassTools",
"%s",
513 Info(
"DiTauMassTools",
"%s",
523 const TLorentzVector *origVisTau1 = 0;
524 const TLorentzVector *origVisTau2 = 0;
537 Info(
"DiTauMassTools",
538 "------------- Printing Final Results for MissingMassCalculator --------------");
539 Info(
"DiTauMassTools",
540 ".............................................................................");
544 Info(
"DiTauMassTools",
"%s",
547 Info(
"DiTauMassTools",
"%s",
553 Info(
"DiTauMassTools",
" no 4-momentum or MET from this method ");
558 Info(
"DiTauMassTools",
" fit failed ");
567 Info(
"DiTauMassTools",
"%s",
573 Info(
"DiTauMassTools",
"%s",
579 Info(
"DiTauMassTools",
"%s",
585 Info(
"DiTauMassTools",
"%s",
592 Info(
"DiTauMassTools",
"%s",
593 (
" dR(nu1-visTau1)=" +
std::to_string(tlvnu1.DeltaR(*origVisTau1))).c_str());
594 Info(
"DiTauMassTools",
"%s",
595 (
" dR(nu2-visTau2)=" +
std::to_string(tlvnu2.DeltaR(*origVisTau2))).c_str());
597 Info(
"DiTauMassTools",
"%s",
624 const double &phi1,
const double &phi2,
625 int &nsol1,
int &nsol2) {
646 int solution_code = 0;
664 if (pTmiss2 * dPhiSign < 0) {
666 return solution_code;
672 if (pTmiss1 * (-dPhiSign) < 0) {
674 return solution_code;
683 double m4noma1 = m2noma1 * m2noma1;
685 double p2v1proj =
std::pow(pv1proj, 2);
687 double pTmiss2CscDPhi = pTmiss2 / sinDPhi2;
688 double &pTn1 = pTmiss2CscDPhi;
689 double pT2miss2CscDPhi = pTmiss2CscDPhi * pTmiss2CscDPhi;
692 const double discri1 = m4noma1 + 4 * m2noma1 * pTmiss2CscDPhi * pv1proj -
693 4 * (
m_ET2v1 * (
m_m2Nu1 + pT2miss2CscDPhi) - (pT2miss2CscDPhi * p2v1proj));
698 return solution_code;
706 double m4noma2 = m2noma2 * m2noma2;
708 double p2v2proj =
std::pow(pv2proj, 2);
709 double sinDPhi1 = -sinDPhi2;
710 double pTmiss1CscDPhi = pTmiss1 / sinDPhi1;
711 double &pTn2 = pTmiss1CscDPhi;
712 double pT2miss1CscDPhi = pTmiss1CscDPhi * pTmiss1CscDPhi;
714 const double discri2 = m4noma2 + 4 * m2noma2 * pTmiss1CscDPhi * pv2proj -
715 4 * (
m_ET2v2 * (
m_m2Nu2 + pT2miss1CscDPhi) - (pT2miss1CscDPhi * p2v2proj));
720 return solution_code;
726 double sqdiscri1 = sqrt(discri1);
732 double pn1Z = first1 + second1;
734 if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z *
m_tauVec1Pz >
743 pn1Z = first1 - second1;
745 if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z *
m_tauVec1Pz >
757 return solution_code;
762 double sqdiscri2 = sqrt(discri2);
768 double pn2Z = first2 + second2;
770 if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z *
m_tauVec2Pz >
779 pn2Z = first2 - second2;
782 if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z *
m_tauVec2Pz >
793 return solution_code;
809 if (std::abs(pnux - pTmissx) > 0.001 || std::abs(pnuy - pTmissy) > 0.001) {
810 Info(
"DiTauMassTools",
"%s", (
"NuPsolutionV3 ERROR Pnux-Met.X or Pnuy-Met.Y > 0.001 : " +
815 if (std::abs(mtau1plus -
m_mTau) > 0.001 || std::abs(mtau1moins -
m_mTau) > 0.001 ||
816 std::abs(mtau2plus -
m_mTau) > 0.001 || std::abs(mtau2moins -
m_mTau) > 0.001) {
817 Info(
"DiTauMassTools",
"%s", (
"NuPsolutionV3 ERROR tau mass not recovered : " +
824 return solution_code;
829 const TLorentzVector &tau,
const double &l_nu,
830 std::vector<TLorentzVector> &nu_vec) {
831 int solution_code = 0;
834 TLorentzVector nu(0.0, 0.0, 0.0, 0.0);
835 TLorentzVector nu2(0.0, 0.0, 0.0, 0.0);
836 nu.SetXYZM(met_vec.Px(), met_vec.Py(), 0.0, l_nu);
837 nu2.SetXYZM(met_vec.Px(), met_vec.Py(), 0.0, l_nu);
839 const double Mtau = 1.777;
841 double msq = (Mtau * Mtau - tau.M() * tau.M() - l_nu * l_nu) /
843 double gamma = nu.Px() * nu.Px() + nu.Py() * nu.Py();
844 double beta = tau.Px() * nu.Px() + tau.Py() * nu.Py() + msq;
845 double a = tau.E() * tau.E() - tau.Pz() * tau.Pz();
846 double b = -2 * tau.Pz() *
beta;
848 if ((
b *
b - 4 *
a *
c) < 0)
849 return solution_code;
852 double pvz1 = (-
b + sqrt(
b *
b - 4 *
a *
c)) / (2 *
a);
853 double pvz2 = (-
b - sqrt(
b *
b - 4 *
a *
c)) / (2 *
a);
854 nu.SetXYZM(met_vec.Px(), met_vec.Py(), pvz1, l_nu);
855 nu2.SetXYZM(met_vec.Px(), met_vec.Py(), pvz2, l_nu);
856 nu_vec.push_back(nu);
857 nu_vec.push_back(nu2);
858 return solution_code;
873 #ifdef SAVELIKELIHOODHISTO
914 TVector2 deltamet_vec;
921 bool paramInsideRange =
false;
940 if (paramInsideRange)
964 if (nsuccesses > 0) {
970 double Px1, Py1, Pz1;
971 double Px2, Py2, Pz2;
972 if (nsuccesses > 0) {
993 if (prob_hist != 0.0)
1017 TLorentzVector fulltau1, fulltau2;
1019 fulltau1.SetXYZM(Px1, Py1, Pz1, 1.777);
1020 fulltau2.SetXYZM(Px2, Py2, Pz2, 1.777);
1037 Info(
"DiTauMassTools",
"Scanning ");
1038 Info(
"DiTauMassTools",
" Markov ");
1039 Info(
"DiTauMassTools",
"%s",
1045 Info(
"DiTauMassTools",
"%s",
1051 if (fit_code == 0) {
1052 Info(
"DiTauMassTools",
"%s", (
"!!!----> Warning-3 in "
1053 "MissingMassCalculator::DitauMassCalculatorV9Walk() : fit status=" +
1056 Info(
"DiTauMassTools",
"%s",
"....... No solution is found. Printing input info .......");
1074 Info(
"DiTauMassTools",
" ---------------------------------------------------------- ");
1096 const double Mtau = 1.777;
1112 double METresX_binSize = 2 * N_METsigma * METresX / NiterMET;
1113 double METresY_binSize = 2 * N_METsigma * METresY / NiterMET;
1117 std::vector<TLorentzVector> nu_vec;
1122 double metprob = 1.0;
1123 double sign_tmp = 0.0;
1124 double tauprob = 1.0;
1125 double totalProb = 0.0;
1129 double met_smear_x = 0.0;
1130 double met_smear_y = 0.0;
1131 double met_smearL = 0.0;
1132 double met_smearP = 0.0;
1134 double angle1 = 0.0;
1174 Info(
"DiTauMassTools",
"Running in dilepton mode");
1179 TLorentzVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1180 TLorentzVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1220 double Mlep = tau_tmp.M();
1227 double MnuProb = 1.0;
1229 for (
int i3 = 0; i3 < NiterMnu; i3++)
1231 M_nu = Mnu_binSize * i3;
1232 if (M_nu >= (Mtau - Mlep))
1238 for (
int i4 = 0; i4 < NiterMET + 1; i4++)
1240 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1241 for (
int i5 = 0; i5 < NiterMET + 1; i5++)
1243 met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1244 if (
pow(met_smearL / METresX, 2) +
pow(met_smearP / METresY, 2) >
pow(N_METsigma, 2))
1246 met_smear_x = met_smearL * met_coscovphi - met_smearP * met_sincovphi;
1247 met_smear_y = met_smearL * met_sincovphi + met_smearP * met_coscovphi;
1248 metvec_tmp.Set(input_metX + met_smear_x, input_metY + met_smear_y);
1266 for (
unsigned int j1 = 0;
j1 < nu_vec.size();
j1++) {
1269 const double tau1_tmpp = (tau_tmp + nu_vec[
j1]).P();
1270 angle1 =
Angle(nu_vec[
j1], tau_tmp);
1280 double tauvecprob1j =
1282 if (tauvecprob1j == 0.)
1285 totalProb = tauvecprob1j * metprob * MnuProb * tauprob;
1302 m_fPXfit1->Fill((tau_tmp + nu_vec[
j1]).Px(), totalProb);
1304 m_fPZfit1->Fill((tau_tmp + nu_vec[
j1]).Pz(), totalProb);
1308 sign_tmp = -log10(totalProb);
1332 Info(
"DiTauMassTools",
"Running in lepton+tau mode");
1358 TLorentzVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1359 TLorentzVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1373 for (
int i4 = 0; i4 < NiterMET + 1; i4++)
1375 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1376 for (
int i5 = 0; i5 < NiterMET + 1; i5++)
1378 met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1379 if (
pow(met_smearL / METresX, 2) +
pow(met_smearP / METresY, 2) >
pow(N_METsigma, 2))
1383 metvec_tmp.Set(input_metX + met_smear_x, input_metY + met_smear_y);
1401 for (
unsigned int j1 = 0;
j1 < nu_vec.size();
j1++) {
1404 const double tau1_tmpp = (tau_tmp + nu_vec[
j1]).P();
1405 angle1 =
Angle(nu_vec[
j1], tau_tmp);
1415 double tauvecprob1j =
1417 if (tauvecprob1j == 0.)
1420 totalProb = tauvecprob1j * metprob * tauprob;
1443 sign_tmp = -log10(totalProb);
1469 Info(
"DiTauMassTools",
"Running in an unknown mode?!?!");
1476 Info(
"DiTauMassTools",
"%s",
1498 if (prob_hist != 0.0)
1516 TLorentzVector nu1_tmp(0.0, 0.0, 0.0, 0.0);
1517 TLorentzVector nu2_tmp(0.0, 0.0, 0.0, 0.0);
1520 nu2_tmp.SetXYZM(Px1, Py1, Pz1, 1.777);
1524 nu1_tmp.SetXYZM(Px1, Py1, Pz1, 1.777);
1537 if (fit_code == 0) {
1539 "DiTauMassTools",
"%s",
1540 (
"!!!----> Warning-3 in MissingMassCalculator::DitauMassCalculatorV9lfv() : fit status=" +
1543 Info(
"DiTauMassTools",
"....... No solution is found. Printing input info .......");
1555 Info(
"DiTauMassTools",
" ---------------------------------------------------------- ");
1566 const double mM =
x[0];
1567 const double mMax =
par[0];
1568 const double mMean =
par[1];
1569 const double mInvWidth2 =
par[2];
1571 const double fitval = mMax * (1 - 4 * mInvWidth2 *
std::pow(mM - mMean, 2));
1584 const int winHalfWidth,
bool debug) {
1601 winHalfWidth == 0)) {
1605 int max_bin = theHist->GetMaximumBin();
1606 maxPos = theHist->GetBinCenter(max_bin);
1609 prob = theHist->GetBinContent(max_bin) /
double(theHist->GetEntries());
1616 int hNbins = theHist->GetNbinsX();
1621 int max_bin = theHist->GetMaximumBin();
1622 int iBinMin = max_bin - winHalfWidth;
1625 int iBinMax = max_bin + winHalfWidth;
1626 if (iBinMax > hNbins)
1627 iBinMax = hNbins - 1;
1630 for (
int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1631 const double weight = theHist->GetBinContent(iBin);
1633 sumx +=
weight * theHist->GetBinCenter(iBin);
1635 maxPos = sumx /
sumw;
1638 prob =
sumw / theHist->GetEntries();
1648 Error(
"DiTauMassTools",
"%s",
1649 (
"ERROR undefined maxHistStrategy:" +
std::to_string(maxHistStrategy)).c_str());
1655 int lastNonZeroBin = -1;
1656 int firstNonZeroBin = -1;
1657 double totalSumw = 0.;
1658 bool firstNullPart =
true;
1659 for (
int iBin = 0; iBin < hNbins; ++iBin) {
1660 const double weight = theHist->GetBinContent(iBin);
1663 lastNonZeroBin = iBin;
1664 if (firstNullPart) {
1665 firstNullPart =
false;
1666 firstNonZeroBin = iBin;
1673 firstNonZeroBin =
std::max(0, firstNonZeroBin - winHalfWidth - 1);
1674 lastNonZeroBin =
std::min(hNbins - 1, lastNonZeroBin + winHalfWidth + 1);
1683 const int nwidth = 2 * winHalfWidth + 1;
1686 for (
int ibin = 0; ibin < nwidth; ++ibin) {
1687 winsum += theHist->GetBinContent(ibin);
1689 double winmax = winsum;
1692 int iBinL = firstNonZeroBin;
1693 int iBinR = iBinL + 2 * winHalfWidth;
1694 bool goingUp =
true;
1699 const double deltawin = theHist->GetBinContent(iBinR) - theHist->GetBinContent(iBinL - 1);
1705 if (winsum > winmax) {
1708 max_bin = (iBinR + iBinL) / 2 - 1;
1719 }
while (iBinR < lastNonZeroBin);
1722 int iBinMin = max_bin - winHalfWidth;
1725 int iBinMax = max_bin + winHalfWidth;
1726 if (iBinMax >= hNbins)
1727 iBinMax = hNbins - 1;
1730 for (
int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1731 const double weight = theHist->GetBinContent(iBin);
1733 sumx +=
weight * theHist->GetBinCenter(iBin);
1736 double maxPosWin = -1.;
1739 maxPosWin = sumx /
sumw;
1746 const double h_rms = theHist->GetRMS(1);
1750 double numerator = 0;
1752 bool nullBin =
false;
1754 for (
int i = iBinMin;
i < iBinMax; ++
i) {
1755 double binError = theHist->GetBinError(
i);
1756 if (binError < 1
e-10) {
1759 double binErrorSquare =
std::pow(binError, 2);
1760 num = theHist->GetBinContent(
i) / (binErrorSquare);
1761 numerator = numerator +
num;
1764 if (numerator < 1
e-10 ||
denominator < 1
e-10 || nullBin ==
true) {
1780 const double binWidth = theHist->GetBinCenter(2) - theHist->GetBinCenter(1);
1781 double fitWidth = (winHalfWidth + 0.5) *
binWidth;
1794 TString fitOption =
debug ?
"QS" :
"QNS";
1800 TFitResultPtr fitRes =
1801 theHist->Fit(
m_fFitting, fitOption,
"", maxPos - fitWidth, maxPos + fitWidth);
1803 double maxPosFit = -1.;
1805 if (
int(fitRes) == 0) {
1808 const double mMax = fitRes->Parameter(0);
1809 const double mMean = fitRes->Parameter(1);
1810 const double mInvWidth2 = fitRes->Parameter(2);
1811 double mMaxError = fitRes->ParError(0);
1813 double mMeanError = fitRes->ParError(1);
1815 double mInvWidth2Error = fitRes->ParError(2);
1818 mInvWidth2Error = 0.;
1819 const double c = mMax * (1 - 4 * mMean * mMean * mInvWidth2);
1820 const double b = 8 * mMax * mMean * mInvWidth2;
1821 const double a = -4 * mMax * mInvWidth2;
1827 const double h_discri =
b *
b - 4 *
a *
c;
1829 const double sqrth_discri = sqrt(h_discri);
1830 const double h_fitLength = sqrth_discri /
a;
1837 maxPosFit = -
b / (2 *
a);
1841 if (maxPosFit >= 0. and std::abs(maxPosFit - maxPosWin) < 0.8 * fitWidth) {
1859 const double &M_nu1,
1860 const double &M_nu2) {
1866 const int solution =
NuPsolutionV3(M_nu1, M_nu2, phi1, phi2, nsol1, nsol2);
1885 const int nsol1,
const int nsol2,
1886 const double &Mvis,
const double &Meff)
1892 Error(
"DiTauMassTools",
"%s",
1896 Error(
"DiTauMassTools",
"%s",
1900 Error(
"DiTauMassTools",
"%s",
1904 Error(
"DiTauMassTools",
"%s",
1908 Error(
"DiTauMassTools",
"%s", (
"refineSolutions ERROR nsol1 " +
std::to_string(nsol1) +
1912 Error(
"DiTauMassTools",
"%s", (
"refineSolutions ERROR nsol1 " +
std::to_string(nsol2) +
1919 Prob->
apply(
preparedInput, -99, -99, TLorentzVector(0, 0, 0, 0), TLorentzVector(0, 0, 0, 0),
1920 TLorentzVector(0, 0, 0, 0), TLorentzVector(0, 0, 0, 0),
true,
false,
false);
1922 for (
int j1 = 0;
j1 < nsol1; ++
j1) {
1932 const int pickInt = std::abs(10000 *
m_Phi1);
1933 const int pickDigit = pickInt - 10 * (pickInt / 10);
1941 nuvec1_tmpj.SetXYZM(nuvec1_tmpj.Px(), nuvec1_tmpj.Py(), nuvec1_tmpj.Pz(), M_nu1);
1942 tauvecsol1j.SetPxPyPzE(0., 0., 0., 0.);
1943 tauvecsol1j += nuvec1_tmpj;
1948 TLorentzVector(0, 0, 0, 0), nuvec1_tmpj,
1949 TLorentzVector(0, 0, 0, 0),
false,
true,
false);
1953 for (
int j2 = 0;
j2 < nsol2; ++
j2) {
1964 const int pickInt = std::abs(10000 *
m_Phi2);
1965 const int pickDigit = pickInt - 10 *
int(pickInt / 10);
1973 nuvec2_tmpj.SetXYZM(nuvec2_tmpj.Px(), nuvec2_tmpj.Py(), nuvec2_tmpj.Pz(), M_nu2);
1974 tauvecsol2j.SetPxPyPzE(0., 0., 0., 0.);
1975 tauvecsol2j += nuvec2_tmpj;
1981 TLorentzVector(0, 0, 0, 0), nuvec2_tmpj,
false,
true,
false);
1985 if (tauvecprob1j == 0.)
1987 if (tauvecprob2j == 0.)
1990 double totalProb = 1.;
2003 (constProb * tauvecprob1j * tauvecprob2j *
2007 if (totalProb <= 0) {
2009 Warning(
"DiTauMassTools",
"%s",
2010 (
"null proba solution, rejected "+
std::to_string(totalProb)).c_str());
2017 Error(
"DiTauMassTools",
"%s",
2018 (
"refineSolutions ERROR nsol getting larger than nsolfinalmax!!! " +
2021 Error(
"DiTauMassTools",
"%s",
2037 nu1Final.SetPxPyPzE(nuvec1_tmpj.Px(), nuvec1_tmpj.Py(), nuvec1_tmpj.Pz(), nuvec1_tmpj.E());
2038 nu2Final.SetPxPyPzE(nuvec2_tmpj.Px(), nuvec2_tmpj.Py(), nuvec2_tmpj.Pz(), nuvec2_tmpj.E());
2046 if (ngoodsol1 == 0) {
2049 if (ngoodsol2 == 0) {
2056 const TLorentzVector &nu1,
2057 const TLorentzVector &vis2,
2058 const TLorentzVector &nu2,
const double &mmc_mass,
2059 const double &vis_mass,
const double &eff_mass,
2060 const double &dphiTT) {
2071 const double MrecoMvis = mmc_mass / vis_mass;
2072 if (MrecoMvis > 2.6)
2074 const double MrecoMeff = mmc_mass / eff_mass;
2075 if (MrecoMeff > 1.9)
2077 const double e1p1 = nu1.E() / vis1.P();
2078 const double e2p2 = nu2.E() / vis2.P();
2079 if ((e1p1 + e2p2) > 4.5)
2099 const double MrecoMvis = mmc_mass / vis_mass;
2100 const double MrecoMeff = mmc_mass / eff_mass;
2101 const double x = dphiTT > 1.5 ? dphiTT : 1.5;
2102 if ((MrecoMeff + MrecoMvis) > 5.908 - 1.881 *
x + 0.2995 *
x *
x)
2111 const int &tau_type2) {
2112 double Fscale = 1.0;
2134 float p0, p1, p2,
p3, p4, p5, p6, p7;
2135 if ((tau_type1 >= 0 && tau_type1 <= 2) || (tau_type2 >= 0 && tau_type2 <= 2))
2137 if ((tau_type1 >= 3 && tau_type1 <= 5) || (tau_type2 >= 3 && tau_type2 <= 5))
2148 double scale3 = p4 / (p5 + p6 *
mass) + p7;
2150 Fscale = scale3 /
scale1;
2152 scale1 = p0 / (p1 + p2 * 91.2) +
p3;
2153 scale3 = p4 / (p5 + p6 * 91.2) + p7;
2154 Fscale = scale3 /
scale1;
2163 return 1.0 / Fscale;
2172 double totalProbSumSol = 0.;
2173 double totalProbSumSolOld = 0.;
2174 bool firstPointWithSol =
false;
2176 for (
int isol = 0; isol <
m_nsol; ++isol) {
2181 bool notSureToKeep =
true;
2185 notSureToKeep =
false;
2192 firstPointWithSol =
true;
2204 if (notSureToKeep) {
2207 for (
int isol = 0; isol <
m_nsolOld; ++isol) {
2213 if (!firstPointWithSol && totalProbSumSolOld <= 0.) {
2214 Error(
"DiTauMassTools",
"%s",
2215 (
" ERROR null old probability !!! " +
std::to_string(totalProbSumSolOld) +
" nsolOld " +
2219 }
else if (totalProbSumSol > totalProbSumSolOld) {
2224 }
else if (totalProbSumSol < totalProbSumSolOld * 1
E-6) {
2228 }
else if (
m_nsol <= 0) {
2235 reject = (uMC > totalProbSumSol / totalProbSumSolOld);
2260 bool fillSolution =
true;
2261 bool oldToBeUsed =
false;
2269 fillSolution =
false;
2283 if (!firstPointWithSol) {
2284 fillSolution =
true;
2287 fillSolution =
false;
2294 if (!fillSolution) {
2295 if (firstPointWithSol) {
2298 for (
int isol = 0; isol <
m_nsol; ++isol) {
2310 double solSum2 = 0.;
2312 for (
int isol = 0; isol <
m_nsol; ++isol) {
2316 const TLorentzVector *pnuvec1_tmpj;
2317 const TLorentzVector *pnuvec2_tmpj;
2330 const TLorentzVector &nuvec1_tmpj = *pnuvec1_tmpj;
2331 const TLorentzVector &nuvec2_tmpj = *pnuvec2_tmpj;
2334 solSum2 += mtautau * mtautau;
2348 #ifdef SAVELIKELIHOODHISTO
2355 if (mtautau != 0. &&
weight != 0.)
2419 for (
int isol = 0; isol <
m_nsol; ++isol) {
2618 Info(
"DiTauMassTools",
" in m_meanbinToBeEvaluated && m_iterNsuc==500 ");
2634 int stopint = stopdouble;
2765 TLorentzVector Met4vec;
2838 const double &P_tau) {
2840 #ifndef WITHDTHETA3DLIM
2842 if (limit_code == 0)
2844 if (limit_code == 1)
2846 if (limit_code == 2)
2851 if (limit_code == 0)
2853 double par[3] = {0.0, 0.0, 0.0};
2855 if (tau_type == 8) {
2856 if (limit_code == 0)
2862 if (limit_code == 1)
2868 if (limit_code == 2)
2876 if (tau_type >= 0 && tau_type <= 2) {
2877 if (limit_code == 0)
2881 par[2] = -0.0004859;
2883 if (limit_code == 1)
2889 if (limit_code == 2)
2897 if (tau_type >= 3 && tau_type <= 5) {
2898 if (limit_code == 0)
2902 par[2] = -0.0009458;
2904 if (limit_code == 1)
2910 if (limit_code == 2)
2918 if (std::abs(P_tau +
par[1]) > 0.0)
2920 if (limit_code == 0) {
2923 }
else if (
limit > 0.03) {
2927 if (limit < 0.0 || limit > 0.5 * TMath::Pi()) {
2928 limit = 0.5 * TMath::Pi();
2929 }
else if (limit < 0.05 && limit > 0.0) {
2940 double theta2,
double phi2,
double &P1,
2942 int solution_code = 0;
2945 double D =
sin(theta1) *
sin(theta2) *
sin(phi2 - phi1);
2946 if (std::abs(D) > 0.0)
2948 P1 = (met_vec.Px() *
sin(phi2) - met_vec.Py() *
cos(phi2)) *
sin(theta2) / D;
2949 P2 = (met_vec.Py() *
cos(phi1) - met_vec.Px() *
sin(phi1)) *
sin(theta1) / D;
2950 if (P1 > 0.0 && P2 > 0.0)
2953 return solution_code;
2960 const TLorentzVector &tau_vec2,
2961 const TVector2 &met_vec,
double &Mrec) {
2966 int coll_code =
NuPsolution(met_vec, tau_vec1.Theta(), tau_vec1.Phi(), tau_vec2.Theta(),
2967 tau_vec2.Phi(), P_nu1, P_nu2);
2968 if (coll_code == 1) {
2970 TLorentzVector nu1(P_nu1 *
sin(tau_vec1.Theta()) *
cos(tau_vec1.Phi()),
2971 P_nu1 *
sin(tau_vec1.Theta()) *
sin(tau_vec1.Phi()),
2972 P_nu1 *
cos(tau_vec1.Theta()), P_nu1);
2973 TLorentzVector nu2(P_nu2 *
sin(tau_vec2.Theta()) *
cos(tau_vec2.Phi()),
2974 P_nu2 *
sin(tau_vec2.Theta()) *
sin(tau_vec2.Phi()),
2975 P_nu2 *
cos(tau_vec2.Theta()), P_nu2);
2976 Mrec = (nu1 + nu2 + tau_vec1 + tau_vec2).M();
2988 const double GEV = 1000.;
3002 if (LFVMode == -1) {
3004 }
else if (LFVMode != -2) {
3011 TLorentzVector tlvTau1 =
part1->p4();
3012 TLorentzVector tlvTau2 =
part2->p4();
3016 TLorentzVector fixedtau1;
3017 fixedtau1.SetPtEtaPhiM(tlvTau1.Pt() /
GEV, tlvTau1.Eta(), tlvTau1.Phi(), tlvTau1.M() /
GEV);
3018 TLorentzVector fixedtau2;
3019 fixedtau2.SetPtEtaPhiM(tlvTau2.Pt() /
GEV, tlvTau2.Eta(), tlvTau2.Phi(), tlvTau2.M() /
GEV);
3026 if (mmcType1 == 8 && mmcType2 == 8) {
3028 }
else if (mmcType1 >= 0 && mmcType1 <= 5 && mmcType2 >= 0 && mmcType2 <= 5) {
3044 Error(
"DiTauMassTools",
"MMCCalibrationSet has not been set !. Please use "
3045 "fMMC.SetCalibrationSet(MMCCalibrationSetV2::MMC2019)"
3087 Info(
"DiTauMassTools",
"correcting sumET");
3184 Info(
"DiTauMassTools",
"Substracting pt1 from sumEt");
3191 Info(
"DiTauMassTools",
"Substracting pt2 from sumEt");
3234 double HtOffset = 0.;
3239 HtOffset = 87.5 - 27.0 *
x;