From 757734410c972032a93396a61dc698936783fe2d Mon Sep 17 00:00:00 2001
From: Antonio De Maria <antonio.de.maria@cern.ch>
Date: Sat, 20 Jul 2024 07:38:44 +0200
Subject: [PATCH 1/2] Remove unused function in MMC

---
 .../DiTauMassTools/MissingMassCalculatorV2.h  |  11 --
 .../Root/MissingMassCalculatorV2.cxx          | 110 +-----------------
 2 files changed, 1 insertion(+), 120 deletions(-)

diff --git a/PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/MissingMassCalculatorV2.h b/PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/MissingMassCalculatorV2.h
index a3dd64f6856d..e5981213c6b8 100644
--- a/PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/MissingMassCalculatorV2.h
+++ b/PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/MissingMassCalculatorV2.h
@@ -243,8 +243,6 @@ class MissingMassCalculatorV2 {
   DitauStuff m_fDitauStuffFit; // results based on fit method
   DitauStuff m_fDitauStuffHisto; // results based on histo method
 
-  int m_fApplyMassScale; // switch to apply mass scale correction
-
   int m_niter_fit1; // number of iterations for dR-dPhi scan 
   int m_niter_fit2; // number of iterations for MET-scan 
   int m_niter_fit3; // number of iterations for Mnu-scan 
@@ -253,8 +251,6 @@ class MissingMassCalculatorV2 {
   int m_RMSStop;
   int m_RndmSeedAltering; // reset seed (not necessary by default)
 
-
-  int m_fJERsyst; // switch for JER systematics
   double m_dTheta3d_binMin; // minimal step size for dTheta3D
   double m_dTheta3d_binMax; // maximum step size for dTheta3D
   double m_dRmax_tau; // maximum dR(nu-visTau)
@@ -267,9 +263,6 @@ class MissingMassCalculatorV2 {
   void DoOutputInfo();
   void PrintOtherInput();
   void PrintResults();
-  int NuPsolution(TVector2 met_vec, double theta1, double phi1, 
-		  double theta2, double phi2, double &P1, double &P2); // keep this version for simple tests
-
 
   inline int NuPsolutionV3(const double & mNu1, const double & mNu2, const double & phi1, const double & phi2, 
 			   int & nsol1, int & nsol2);
@@ -358,9 +351,6 @@ public:
   void SetdTheta3d_binMin(const double val) { m_dTheta3d_binMin=val; } // minimal step size for dTheta3D
   void SetEventNumber(const int eventNumber) { m_eventNumber = eventNumber; }
 
-  void SetJERsyst(const int val) { m_fJERsyst=val; }
-  void SetApplyMassScale(const int val) { m_fApplyMassScale=val; } 
-
   void SetMnuScanRange(const double val) { m_MnuScanRange=val; }
 
   void SetProposalTryMEt(const double val) {m_proposalTryMEt=val; }
@@ -410,7 +400,6 @@ public:
 
 
   //-------- Get results;
-  int StandardCollApprox(const TLorentzVector & tau_vec1, const TLorentzVector & tau_vec2, const TVector2 & met_vec, double &Mrec); // standard collinear approximation
   Double_t maxFitting(Double_t *x, Double_t *par);
   
   // compute maximum from histo
diff --git a/PhysicsAnalysis/TauID/DiTauMassTools/Root/MissingMassCalculatorV2.cxx b/PhysicsAnalysis/TauID/DiTauMassTools/Root/MissingMassCalculatorV2.cxx
index b8b90eb62939..29cca5bc756b 100644
--- a/PhysicsAnalysis/TauID/DiTauMassTools/Root/MissingMassCalculatorV2.cxx
+++ b/PhysicsAnalysis/TauID/DiTauMassTools/Root/MissingMassCalculatorV2.cxx
@@ -74,9 +74,7 @@ MissingMassCalculatorV2::MissingMassCalculatorV2(
   Prob->SetUseDphiLL(false);         // added by Tomas Davidek for lep-lep
   m_dTheta3d_binMin = 0.0025;
   m_dTheta3d_binMax = 0.02;
-  m_fJERsyst = 0;                 // no JER systematics by default (+/-1: up/down 1 sigma)
   preparedInput.m_METresSyst = 0; // no MET resolution systematics by default (+/-1: up/down 1 sigma)
-  m_fApplyMassScale = 0;          // don't apply mass scale correction by default
   preparedInput.m_dataType = 1;   // set to "data" by default
   preparedInput.m_fUseTailCleanup = 1; // cleanup by default for lep-had Moriond 2012 analysis
   preparedInput.m_fUseDefaults = 0; // use pre-set defaults for various configurations; if set it to 0
@@ -401,13 +399,7 @@ void MissingMassCalculatorV2::DoOutputInfo() {
                          OutputInfo.m_nuvec2[MMCFitMethodV2::MAXW].Py());
     OutputInfo.m_FittedMetVec[MMCFitMethodV2::MAXW] = metmaxw;
 
-    // MLM method : can only get MMC, rest is dummy
-    double scale = MassScale(MMCFitMethodV2::MAXW, m_fDitauStuffHisto.Mditau_best,
-                             preparedInput.m_type_visTau1,
-                             preparedInput.m_type_visTau2); // only for histo method for now. In
-                                                          // practice disabled by default
-
-    OutputInfo.m_FittedMass[MMCFitMethodV2::MLM] = scale * m_fDitauStuffHisto.Mditau_best;
+    OutputInfo.m_FittedMass[MMCFitMethodV2::MLM] = m_fDitauStuffHisto.Mditau_best;
     OutputInfo.m_FittedMassLowerError[MMCFitMethodV2::MLM] = yq[0];
     OutputInfo.m_FittedMassUpperError[MMCFitMethodV2::MLM] = yq[1];
 
@@ -2106,63 +2098,6 @@ int MissingMassCalculatorV2::TailCleanUp(const TLorentzVector &vis1,
   return pass_code;
 }
 
-//-------- This function applies correction to compensate for the off-set
-double MissingMassCalculatorV2::MassScale(int method, double mass, const int &tau_type1,
-                                                   const int &tau_type2) {
-  double Fscale = 1.0;
-  // calibration for rel16 lep-had analysis only
-  if (m_fApplyMassScale == 1) {
-    if (preparedInput.m_tauTypes == TauTypes::lh) {
-      if (method != 1)
-        return 1.0;
-      //      float p0, p1, p2, p3;
-      //      if(tau_type1==1 || tau_type2==1) // 1-prong tau's
-      //        {
-      //          p0=3.014; p1=-71.86; p2=1.018; p3=0.8912;
-      //          if(mass>91.2) Fscale=p0/(p1+p2*mass)+p3;
-      //          else Fscale=p0/(p1+p2*91.2)+p3;
-      //        }
-      //      if(tau_type1==3 || tau_type2==3) // 3-prong tau's
-      //        {
-      //          p0=0.4576; p1=-84.22; p2=0.9783; p3=0.9136;
-      //          if(mass>91.2) Fscale=p0/(p1+p2*mass)+p3;
-      //          else Fscale=p0/(p1+p2*91.2)+p3;
-      //        }
-      //      if(Fscale>1.0) Fscale=1.0;
-      //      if(Fscale<0.89) Fscale=0.89;
-
-      float p0, p1, p2, p3, p4, p5, p6, p7;
-      if ((tau_type1 >= 0 && tau_type1 <= 2) || (tau_type2 >= 0 && tau_type2 <= 2))
-        return 1.0;                                                                 // 1-prong tau's
-      if ((tau_type1 >= 3 && tau_type1 <= 5) || (tau_type2 >= 3 && tau_type2 <= 5)) // 3-prong tau's
-      {
-        p0 = 3.014;
-        p1 = -71.86;
-        p2 = 1.018;
-        p3 = 0.8912;
-        p4 = 0.4576;
-        p5 = -84.22;
-        p6 = 0.9783;
-        p7 = 0.9136;
-        double scale1 = p0 / (p1 + p2 * mass) + p3;
-        double scale3 = p4 / (p5 + p6 * mass) + p7;
-        if (mass > 91.2)
-          Fscale = scale3 / scale1;
-        else {
-          scale1 = p0 / (p1 + p2 * 91.2) + p3;
-          scale3 = p4 / (p5 + p6 * 91.2) + p7;
-          Fscale = scale3 / scale1;
-        }
-      }
-      if (Fscale > 1.0)
-        Fscale = 1.0;
-      if (Fscale < 0.95)
-        Fscale = 0.95;
-    }
-  }
-  return 1.0 / Fscale;
-}
-
 // note that if MarkovChain the input solutions can be modified
 void MissingMassCalculatorV2::handleSolutions()
 
@@ -2935,49 +2870,6 @@ double MissingMassCalculatorV2::dTheta3DLimit(const int &tau_type, const int &li
   return limit;
 }
 
-// returns P(nu1) & P(nu2)
-int MissingMassCalculatorV2::NuPsolution(TVector2 met_vec, double theta1, double phi1,
-                                                  double theta2, double phi2, double &P1,
-                                                  double &P2) {
-  int solution_code = 0; // 0== no solution, 1==with solution
-  P1 = 0.0;
-  P2 = 0.0;
-  double D = sin(theta1) * sin(theta2) * sin(phi2 - phi1);
-  if (std::abs(D) > 0.0) // matrix deteriminant is non-zero
-  {
-    P1 = (met_vec.Px() * sin(phi2) - met_vec.Py() * cos(phi2)) * sin(theta2) / D;
-    P2 = (met_vec.Py() * cos(phi1) - met_vec.Px() * sin(phi1)) * sin(theta1) / D;
-    if (P1 > 0.0 && P2 > 0.0)
-      solution_code = 1;
-  }
-  return solution_code;
-}
-
-// standard collinear approximation
-// it returns code=0 if collinear approximation can't be applied
-// and code=1 and Mrec if collinear approximation was applied
-int MissingMassCalculatorV2::StandardCollApprox(const TLorentzVector &tau_vec1,
-                                                         const TLorentzVector &tau_vec2,
-                                                         const TVector2 &met_vec, double &Mrec) {
-  int code = 0;
-  Mrec = 0.0;
-  double P_nu1 = 0.0;
-  double P_nu2 = 0.0;
-  int coll_code = NuPsolution(met_vec, tau_vec1.Theta(), tau_vec1.Phi(), tau_vec2.Theta(),
-                              tau_vec2.Phi(), P_nu1, P_nu2);
-  if (coll_code == 1) {
-    code = 1;
-    TLorentzVector nu1(P_nu1 * sin(tau_vec1.Theta()) * cos(tau_vec1.Phi()),
-                       P_nu1 * sin(tau_vec1.Theta()) * sin(tau_vec1.Phi()),
-                       P_nu1 * cos(tau_vec1.Theta()), P_nu1);
-    TLorentzVector nu2(P_nu2 * sin(tau_vec2.Theta()) * cos(tau_vec2.Phi()),
-                       P_nu2 * sin(tau_vec2.Theta()) * sin(tau_vec2.Phi()),
-                       P_nu2 * cos(tau_vec2.Theta()), P_nu2);
-    Mrec = (nu1 + nu2 + tau_vec1 + tau_vec2).M();
-  }
-  return code;
-}
-
 // checks units of input variables, converts into [GeV] if needed, make all
 // possible corrections DR new : now a second structure preparedInput is derived
 // from the input one which only has direct user input
-- 
GitLab


From 76a8f5c23cc0d48abad6082f9a66f89457abc8b9 Mon Sep 17 00:00:00 2001
From: Antonio De Maria <antonio.de.maria@cern.ch>
Date: Thu, 25 Jul 2024 16:04:15 +0200
Subject: [PATCH 2/2] further cleaning

---
 .../DiTauMassTools/HelperFunctions.h          |  2 -
 .../DiTauMassTools/Root/HelperFunctions.cxx   | 26 ------
 .../DiTauMassTools/Root/MissingMassProb.cxx   | 90 -------------------
 3 files changed, 118 deletions(-)

diff --git a/PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h b/PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h
index 3b69f2bc36fc..b8ffd41197b6 100644
--- a/PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h
+++ b/PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h
@@ -53,8 +53,6 @@ namespace TauTypes
 // see source file of MissingMassProb for further reasoning
 template <typename T> void ignore(T &&){}
 
-int getFirstBinBelowMax(const std::shared_ptr<TH1F>& hist, double max, double targetVal);
-int getFirstBinAboveMax(const std::shared_ptr<TH1F>& hist, double max, double targetVal);
 double Angle(const TLorentzVector & vec1, const TLorentzVector & vec2);
 bool updateDouble  (const double in, double & out) ;
 void fastSinCos (const double & phi, double & sinPhi, double & cosPhi);
diff --git a/PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx b/PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx
index 2a0f187adc78..a4cdf27c87c4 100644
--- a/PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx
+++ b/PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx
@@ -12,32 +12,6 @@
 
 using namespace DiTauMassTools;
 
-int DiTauMassTools::getFirstBinAboveMax(const std::shared_ptr<TH1F>& hist, double max, double targetVal) 
-{
-  int maxBin = hist->FindBin(max);
-  int targetBin = maxBin;
-  for(int i=maxBin; i<=hist->GetNbinsX(); i++){
-    if(hist->GetBinContent(i)<targetVal){
-      targetBin = i-1;
-      break;
-    }
-  }
-  return targetBin;
-}
-
-int DiTauMassTools::getFirstBinBelowMax(const std::shared_ptr<TH1F>& hist, double max, double targetVal) 
-{
-  int maxBin = hist->FindBin(max);
-  int targetBin = maxBin;
-  for(int i=1; i<=maxBin; i++){
-    if(hist->GetBinContent(i)>=targetVal){
-      targetBin = i;
-      break;
-    }
-  }
-  return targetBin;
-}
-
 double DiTauMassTools::MaxDelPhi(int tau_type, double Pvis, double dRmax_tau)
 {
   return dRmax_tau+0*Pvis*tau_type; // hack to avoid warning
diff --git a/PhysicsAnalysis/TauID/DiTauMassTools/Root/MissingMassProb.cxx b/PhysicsAnalysis/TauID/DiTauMassTools/Root/MissingMassProb.cxx
index 142a220d7ca8..9978cebb868f 100644
--- a/PhysicsAnalysis/TauID/DiTauMassTools/Root/MissingMassProb.cxx
+++ b/PhysicsAnalysis/TauID/DiTauMassTools/Root/MissingMassProb.cxx
@@ -1153,96 +1153,6 @@ double MissingMassProb::dTheta3Dparam(const int & parInd, const int & tau_type,
   return 0.;
 }
 
-/*
-// Sasha: keep this for now, may need in the future
-// returns analytical P(theta)-probability for given tau-topology
-// decayType==1 for leptonic decays and 0 for hadronic decays
-// uses product of probabilities
-double MissingMassProb::AngularProbability(TLorentzVector nu_vec, TLorentzVector vis_vec, int decayType) {
-  double prob=0.0;
-  double M=1.777;
-  double angl=0.0;
-  double P=(vis_vec+nu_vec).P();
-  double V=P/sqrt(P*P+M*M); // tau speed
-  double dA=dRmax_tau/(2.0*Niter_fit1);
-
-  if(decayType==1) // leptonic tau for now
-    {
-      // exact formular for energy probability is sqrt(1-V^2)/(2*V*p_0)*dE
-      double m_1=nu_vec.M();
-      double m_2=vis_vec.M();
-      double E_nu=(M*M+m_1*m_1-m_2*m_2)/(2.0*M);
-      if(E_nu<=m_1) return 0.0;
-      double P_nu=sqrt(pow(E_nu,2)-pow(m_1,2));
-      double prob1=0.5*sqrt(1-V*V)/(P_nu*V); // energy probability
-
-      angl=Angle(vis_vec,vis_vec+nu_vec); // using lepton direction
-      double det1=1.0-V*cos(angl+dA);
-      double det2= (angl-dA)>0.0 ? 1.0-V*cos(angl-dA) : 1.0-V;
-      double prob2=std::abs(1.0/det1-1.0/det2)*(1.0-V*V)/(2.0*V); // using massless approximation for leptons
-      prob=prob1*prob2;
-    }
-  if(decayType==0) // hadronic tau
-    {
-      // exact formula for energy probability is sqrt(1-V^2)/(2*V*p_0)*dE
-      // drop p_0 because it's a contstant for a given hadronic tau
-      double prob1=0.5*sqrt(1-V*V)/V; // "energy" probability
-
-      angl=Angle(nu_vec,vis_vec+nu_vec); // using neutrino direction
-      double det1=1.0-V*cos(angl+dA);
-      double det2= (angl-dA)>0.0 ? 1.0-V*cos(angl-dA) : 1.0-V;
-      double prob2=std::abs(1.0/det1-1.0/det2)*(1.0-V*V)/(2.0*V);
-      prob=prob1*prob2;
-    }
-  return prob;
-}
-
-// returns probability of angle between two tau's
-// assuming massless tau's for now, should be small effect for M>M_Z
-double MissingMassProb::ResonanceProbability(TLorentzVector vec1, TLorentzVector vec2) {
-
-  double prob=1.0;
-  double boson_P=(vec1+vec2).P();
-  if(boson_P==0.0) return 1.0;
-  double boson_E=(vec1+vec2).E();
-  double boson_V=0.0;
-  if(boson_E>0.0) boson_V=boson_P/boson_E;
-  else return 1.0E-10;
-
-  double testMass=(vec1+vec2).M();
-  double m=1.777; // tau mass
-  double E_tau=testMass/2.0;
-  double P_tau=sqrt(pow(E_tau,2)-m*m);
-  double Evis_lim[2];
-  Evis_lim[0]=(E_tau-boson_V*P_tau)/sqrt(1.0-boson_V*boson_V);
-  Evis_lim[1]=(E_tau+boson_V*P_tau)/sqrt(1.0-boson_V*boson_V);
-  if(vec1.E()<Evis_lim[0] || vec1.E()>Evis_lim[1]) return 1.0E-20;
-  if(vec2.E()<Evis_lim[0] || vec2.E()>Evis_lim[1]) return 1.0E-20;
-
-  double prob1=0.5*sqrt(1-boson_V*boson_V)/(P_tau*boson_V);
-
-  if(vec1.P()*vec2.P()>0)
-    {
-      double theta=acos((vec1.Px()*vec2.Px()+vec1.Py()*vec2.Py()+vec1.Pz()*vec2.Pz())/(vec1.P()*vec2.P()));
-      if(boson_V>0.0 && boson_V<1.0)
-        {
-          if(boson_V<cos(theta/2)) return 1.0E-10;
-          double num=(1.0-boson_V*boson_V)*cos(theta/2);
-          double denom=4*boson_V*sin(theta/2)*sin(theta/2)*sqrt(boson_V*boson_V-cos(theta/2)*cos(theta/2));
-          prob=num/denom;
-        }
-      else
-        {
-          if(std::abs(theta-TMath::Pi())>0.0001) return 1.0E-10;
-        }
-    }
-  else return 1.0E-10;
-  prob=prob*prob1;
-  if(prob<1.0E-20) prob=1.0E-20;
-  return prob;
-}
-*/
-
 void MissingMassProb::MET(MissingMassInput& preparedInput){
   // compute MET resolution (eventually use HT)
   if(preparedInput.m_METsigmaP<0.1 || preparedInput.m_METsigmaL<0.1)
-- 
GitLab