diff --git a/CMakeLists.txt b/CMakeLists.txt index 29b712df..6fc3f8b8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -287,7 +287,8 @@ include_directories(${CMAKE_SOURCE_DIR}/include) file(GLOB SOURCES_1 ${CMAKE_SOURCE_DIR}/src/*.cxx) file(GLOB SOURCES_2 ${CMAKE_SOURCE_DIR}/src/utility/*.cxx - ${CMAKE_SOURCE_DIR}/src/RecoilCorrections/*.cxx) + ${CMAKE_SOURCE_DIR}/src/RecoilCorrections/*.cxx + ${CMAKE_SOURCE_DIR}/src/SVFit/*.cxx) set(SOURCES ${SOURCES_1} ${SOURCES_2}) add_library(CROWNLIB SHARED ${SOURCES}) target_include_directories(CROWNLIB PRIVATE ${CMAKE_SOURCE_DIR} ${ROOT_INCLUDE_DIRS}) diff --git a/include/SVFit/AuxFunctions.hxx b/include/SVFit/AuxFunctions.hxx new file mode 100644 index 00000000..437059ef --- /dev/null +++ b/include/SVFit/AuxFunctions.hxx @@ -0,0 +1,130 @@ +#ifndef TauAnalysis_ClassicSVfit_svFitAuxFunctions_h +#define TauAnalysis_ClassicSVfit_svFitAuxFunctions_h + +#include "Math/LorentzVector.h" +#include "Math/Vector3D.h" +#include + +#include +#include + +namespace fastmtt { +inline double square(double x) { return x * x; } + +inline double cube(double x) { return x * x * x; } + +inline double fourth(double x) { return x * x * x * x; } + +inline double fifth(double x) { return x * x * x * x * x; } + +inline double sixth(double x) { return x * x * x * x * x * x; } + +inline double seventh(double x) { return x * x * x * x * x * x * x; } + +inline double eigth(double x) { return x * x * x * x * x * x * x * x; } + +const double epsilon = 1E-3; +//----------------------------------------------------------------------------- +// define masses, widths and lifetimes of particles +// relevant for computing values of likelihood functions in SVfit algorithm +// +// NOTE: the values are taken from +// K. Nakamura et al. (Particle Data Group), +// J. Phys. G 37, 075021 (2010) +// +const double electronMass = 0.51100e-3; // GeV +const double electronMass2 = electronMass * electronMass; +const double muonMass = 0.10566; // GeV +const double muonMass2 = muonMass * muonMass; + +const double chargedPionMass = 0.13957; // GeV +const double chargedPionMass2 = chargedPionMass * chargedPionMass; +const double neutralPionMass = 0.13498; // GeV +const double neutralPionMass2 = neutralPionMass * neutralPionMass; + +const double rhoMesonMass = 0.77526; // GeV +const double rhoMesonMass2 = rhoMesonMass * rhoMesonMass; +const double a1MesonMass = 1.230; // GeV +const double a1MesonMass2 = a1MesonMass * a1MesonMass; + +const double tauLeptonMass = 1.77685; // GeV +const double tauLeptonMass2 = tauLeptonMass * tauLeptonMass; +const double tauLeptonMass3 = tauLeptonMass2 * tauLeptonMass; +const double tauLeptonMass4 = tauLeptonMass3 * tauLeptonMass; +const double cTauLifetime = 8.711e-3; // centimeters + +const double hbar_c = 0.1973; // GeV fm +const double ctau = 87.e+9; // tau lifetime = 87 microns, converted to fm +const double GammaTau = hbar_c / ctau; +const double GammaTauToElec = GammaTau * 0.178; // BR(tau -> e) = 17.8% +const double GammaTauToMu = GammaTau * 0.174; // BR(tau -> mu) = 17.4% +const double GammaTauToHad = GammaTau * 0.648; // BR(tau -> hadrons) = 64.8% + +#ifdef XSECTION_NORMALIZATION +const double GF = 1.166e-5; // in units of GeV^-2 +const double GFfactor = square(GF) / square(M_PI); +const double M2 = 16. * M_PI * cube(tauLeptonMass) * GammaTauToHad; +#endif + +const double conversionFactor = + 1.e+10 * square(hbar_c); // conversion factor from GeV^-2 to picobarn = + // 10^-40m//FIX ME store this +const double constFactor = 2. * conversionFactor / eigth(2. * M_PI); +const double matrixElementNorm = + square(M_PI / (tauLeptonMass * + GammaTau)); // CV: multiply matrix element by factor + // (Pi/(mTau GammaTau))^2 from Luca's write-up + +// const double v2 = square(246.22); // GeV^2 +//----------------------------------------------------------------------------- + +/** + \typedef SVfitStandalone::Vector + \brief spacial momentum vector (equivalent to reco::Candidate::Vector) +*/ +typedef ROOT::Math::DisplacementVector3D> + Vector; +/** + \typedef SVfitStandalone::LorentzVector + \brief lorentz vector (equivalent to reco::Candidate::LorentzVector) +*/ +typedef ROOT::Math::LorentzVector> LorentzVector; + +double roundToNdigits(double, int = 3); + +struct GraphPoint { + double x_; + double xErr_; + double y_; + double yErr_; + double mTest_step_; +}; +TGraphErrors *makeGraph(const std::string &, const std::vector &); + +void extractResult(TGraphErrors *, double &, double &, double &, int = 0); + +Vector normalize(const Vector &); +double compScalarProduct(const Vector &, const Vector &); +Vector compCrossProduct(const Vector &, const Vector &); + +double compCosThetaNuNu(double, double, double, double, double, double); +double compPSfactor_tauToLepDecay(double, double, double, double, double, + double, double); +double compPSfactor_tauToHadDecay(double, double, double, double, double, + double); + +const int maxNumberOfDimensions = 6; +const int numberOfLegs = 2; + +struct integrationParameters { + + int idx_X_ = -1; + int idx_phi_ = -1; + int idx_VisPtShift_ = -1; + int idx_mNuNu_ = -1; + + void reset(); +}; + +} // namespace fastmtt +#endif diff --git a/include/SVFit/FastMTT.hxx b/include/SVFit/FastMTT.hxx new file mode 100644 index 00000000..61d2f0df --- /dev/null +++ b/include/SVFit/FastMTT.hxx @@ -0,0 +1,204 @@ +#ifndef FastMTT_FastMTT_H +#define FastMTT_FastMTT_H + +#include +#include +#include +#include + +namespace fastmtt { +class MeasuredTauLepton; +} + +namespace ROOT { +namespace Math { +class Minimizer; +class Functor; +} // namespace Math +} // namespace ROOT +class TVector2; + +#include "Math/LorentzVector.h" +#include "TBenchmark.h" +#include "TMatrixD.h" + +typedef ROOT::Math::LorentzVector> LorentzVector; + +namespace fastMTT { +double likelihoodFunc(double *x, double *par); + +enum likelihoodComponent { MET, MASS, PX, PY, ENERGY, IP }; +} // namespace fastMTT + +class Likelihood { + + public: + Likelihood(); + + ~Likelihood(); + + double value(const double *x) const; + + void setLeptonInputs(const LorentzVector &aLeg1P4, + const LorentzVector &aLeg2P4, int aLeg1DecayType, + int aLeg2DecayType, int aLeg1DecayMode, + int aLeg2DecayMode); + + void setMETInputs(const LorentzVector &aMET, const TMatrixD &aCovMET); + + void setParameters(const std::vector ¶meters); + + void enableComponent(fastMTT::likelihoodComponent aCompIndex); + + void disableComponent(fastMTT::likelihoodComponent aCompIndex); + + double massLikelihood(const double &m) const; + + double ptLikelihood(const double &pTTauTau, int type) const; + + double metTF(const LorentzVector &metP4, const LorentzVector &nuP4, + const TMatrixD &covMET) const; + + private: + LorentzVector leg1P4, leg2P4; + LorentzVector recoMET; + + TMatrixD covMET; + + double mVis, mVisLeg1, mVisLeg2; + + int leg1DecayType, leg2DecayType; + int leg1DecayMode, leg2DecayMode; + + std::vector parameters; + + /// Bit word coding enabled likelihood components + std::bitset<128> compnentsBitWord; + + // precomputed values used to reduce the number of redundant calculations + double mVis1OverTauSquare; + double mVis2OverTauSquare; + using PowTable = std::array; // first powers of a number + // array with dimensions 2x3x5 used to store the powers of the coefficients + // of two vectors + std::array, 2u> allpTpows; +}; +////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////// +class FastMTT { + + public: + FastMTT(); + + ~FastMTT(); + + void initialize(); + + /// Run fastMTT algorithm for given input + void run(const std::vector &, const double &, + const double &, const TMatrixD &); + + /// Set likelihood shape parameters. Two parameters are expected: + /// power of 1/mVis, and scaling factor of mTest + /// in case of less than two parameters the default values are used: + ///{6, 1.0/1.15}; + void setLikelihoodParams(const std::vector &aPars); + + /// Set the bit word for switching on particular likelihood components + void enableComponent(fastMTT::likelihoodComponent aCompIndex); + + /// Set the bit word for switching off particular likelihood components + void disableComponent(fastMTT::likelihoodComponent aCompIndex); + + /// Retrieve the four momentum corresponding to the likelihood maximum + const LorentzVector &getBestP4() const { return bestP4; } + + /// Retrieve the tau1 four momentum corresponding to the likelihood maximum + const LorentzVector &getTau1P4() const { return tau1P4; } + + /// Retrieve the tau1 four momentum corresponding to the likelihood maximum + const LorentzVector &getTau2P4() const { return tau2P4; } + + /// Retrieve x values corresponding to the likelihood maximim + std::tuple getBestX() const; + + /// Retrieve the best likelihood value + double getBestLikelihood() const; + + /// Retrieve the likelihood value for given x values + double getLikelihoodForX(double *x) const; + + /// Retrieve the CPU timing for given methods + /// Possible values: + /// scan + /// minimize + double getCpuTime(const std::string &method); + + /// Retrieve the CPU timing for given methods + /// Possible values: + /// scan + /// minimize + double getRealTime(const std::string &method); + + private: + static bool + compareLeptons(const fastmtt::MeasuredTauLepton &measuredTauLepton1, + const fastmtt::MeasuredTauLepton &measuredTauLepton2); + + /// Run the minimalization procedure for given inputs. + /// Results are stored in internal variables accesed by + /// relevant get methods. + void minimize(); + + /// Run a scan over x1 and x2 [0,1] rectangle for given inputs. + /// Results are stored in internal variables accesed by + /// relevant get methods. + void scan(); + + // Minimizer types and algorithms. + // minimizerName minimizerAlgorithm + // Minuit /Minuit2 Migrad, Simplex,Combined,Scan (default is + // Migrad) + // Minuit2 Fumili2 + // Fumili + // GSLMultiMin ConjugateFR, ConjugatePR, BFGS, + // BFGS2, SteepestDescent + // GSLMultiFit + // GSLSimAn + // Genetic + std::string minimizerName; + std::string minimizerAlgorithm; + + ROOT::Math::Minimizer *minimizer; + + /// Minimum location + std::vector minimumPosition; + + /// Mimimum value + double minimumValue; + + /// Dimension of minimalization space + unsigned int nVariables; + + /// Names of variables to be minimized + std::vector varNames; + + /// Values of variables to be minimized + + std::vector variables; + + /// Step sizes for each minimized variable + std::vector stepSizes; + + ROOT::Math::Functor *likelihoodFunctor; + + Likelihood myLikelihood; + + LorentzVector tau1P4, tau2P4, bestP4; + + TBenchmark clock; + + int verbosity; +}; + +#endif diff --git a/include/SVFit/MeasuredTauLepton.hxx b/include/SVFit/MeasuredTauLepton.hxx new file mode 100644 index 00000000..6255d706 --- /dev/null +++ b/include/SVFit/MeasuredTauLepton.hxx @@ -0,0 +1,108 @@ +#ifndef TauAnalysis_ClassicSVfit_MeasuredTauLepton_h +#define TauAnalysis_ClassicSVfit_MeasuredTauLepton_h + +#include "AuxFunctions.hxx" + +namespace fastmtt { +class MeasuredTauLepton { + public: + /** + \enum MeasuredTauLepton::kDecayType + \brief enumeration of all tau decay types + */ + enum kDecayType { + kUndefinedDecayType, + kTauToHadDecay, /* < hadronic tau lepton decay */ + kTauToElecDecay, /* < tau lepton decay to electron */ + kTauToMuDecay /* < tau lepton decay to muon */ + }; + + MeasuredTauLepton(); + MeasuredTauLepton(int, double, double, double, double, int = -1); + MeasuredTauLepton(const MeasuredTauLepton &); + ~MeasuredTauLepton(); + + /// return decay type of the tau lepton + int type() const; + + /// return pt of the measured tau lepton in labframe + double pt() const; + /// return pseudo-rapidity of the measured tau lepton in labframe + double eta() const; + /// return azimuthal angle of the measured tau lepton in labframe + double phi() const; + /// return visible mass in labframe + double mass() const; + + /// return visible energy in labframe + double energy() const; + /// return px of the measured tau lepton in labframe + double px() const; + /// return py of the measured tau lepton in labframe + double py() const; + /// return pz of the measured tau lepton in labframe + double pz() const; + + /// return the measured tau lepton momentum in labframe + double p() const; + + /// return decay mode of the reconstructed hadronic tau decay + int decayMode() const; + + /// return the lorentz vector in the labframe + LorentzVector p4() const; + + /// return the momentum vector in the labframe + Vector p3() const; + + /// return auxiliary data-members to speed-up numerical computations + double cosPhi_sinTheta() const; + double sinPhi_sinTheta() const; + double cosTheta() const; + + void roundToNdigits(unsigned int nDigis = 3); + + protected: + /// set visible momentum in all coordinates systems + void initialize(); + + private: + /// decay type + int type_; + + /// visible momentum in labframe (in polar coordinates) + double pt_; + double eta_; + double phi_; + double mass_; + + /// visible momentum in labframe (in cartesian coordinates) + double energy_; + double px_; + double py_; + double pz_; + + /// visible momentum in labframe (magnitude); + double p_; + + /// decay mode (hadronic tau decays only) + int decayMode_; + + /// visible momentum in labframe (four-vector) + LorentzVector p4_; + + /// visible momentum in labframe + Vector p3_; + + /// mass of visible tau decay products (recomputed to reduce rounding + /// errors) + double preciseVisMass_; + + /// auxiliary data-members to speed-up numerical computations + double cosPhi_sinTheta_; + double sinPhi_sinTheta_; + double cosTheta_; +}; +} // namespace fastmtt + +#endif diff --git a/include/quantities.hxx b/include/quantities.hxx index 4ff59416..0f9a1daa 100644 --- a/include/quantities.hxx +++ b/include/quantities.hxx @@ -26,6 +26,16 @@ ROOT::RDF::RNode charge(ROOT::RDF::RNode df, const std::string &outputname, const std::string &chargecolumn); ROOT::RDF::RNode m_vis(ROOT::RDF::RNode df, const std::string &outputname, const std::vector &inputvectors); +ROOT::RDF::RNode +p4_fastmtt(ROOT::RDF::RNode df, const std::string &outputname, + const std::string &pt_1, const std::string &pt_2, + const std::string &eta_1, const std::string &eta_2, + const std::string &phi_1, const std::string &phi_2, + const std::string &mass_1, const std::string &mass_2, + const std::string &met_pt, const std::string &met_phi, + const std::string &met_cov_xx, const std::string &met_cov_xy, + const std::string &met_cov_yy, const std::string &decay_mode_1, + const std::string &decay_mode_2, const std::string &finalstate); ROOT::RDF::RNode pt_vis(ROOT::RDF::RNode df, const std::string &outputname, const std::vector &inputvectors); ROOT::RDF::RNode pzetamissvis(ROOT::RDF::RNode df, diff --git a/src/SVFit/AuxFunctions.cxx b/src/SVFit/AuxFunctions.cxx new file mode 100644 index 00000000..ac548140 --- /dev/null +++ b/src/SVFit/AuxFunctions.cxx @@ -0,0 +1,280 @@ +#include "../../include/SVFit/AuxFunctions.hxx" + +#include +#include +#include + +namespace fastmtt { + +double roundToNdigits(double x, int n) { + double tmp = TMath::Power(10., n); + if (x != 0.) { + tmp /= TMath::Power(10., TMath::Floor(TMath::Log10(TMath::Abs(x)))); + } + double x_rounded = TMath::Nint(x * tmp) / tmp; + // std::cout << ": x = " << x << ", x_rounded = " << + // x_rounded << std::endl; + return x_rounded; +} + +TGraphErrors *makeGraph(const std::string &graphName, + const std::vector &graphPoints) { + // std::cout << ":" << std::endl; + size_t numPoints = graphPoints.size(); + // std::cout << " numPoints = " << numPoints << std::endl; + TGraphErrors *graph = new TGraphErrors(numPoints); + graph->SetName(graphName.data()); + for (size_t iPoint = 0; iPoint < numPoints; ++iPoint) { + const GraphPoint &graphPoint = graphPoints[iPoint]; + graph->SetPoint(iPoint, graphPoint.x_, graphPoint.y_); + graph->SetPointError(iPoint, graphPoint.xErr_, graphPoint.yErr_); + } + return graph; +} + +void extractResult(TGraphErrors *graph, double &mass, double &massErr, + double &Lmax, int verbosity) { + // determine range of mTest values that are within ~2 sigma interval within + // maximum of likelihood function + double x_Lmax = 0.; + double y_Lmax = 0.; + double idxPoint_Lmax = -1; + for (int iPoint = 0; iPoint < graph->GetN(); ++iPoint) { + double x, y; + graph->GetPoint(iPoint, x, y); + if (y > y_Lmax) { + x_Lmax = x; + y_Lmax = y; + idxPoint_Lmax = iPoint; + } + } + + double xMin = 1.e+6; + double xMax = 0.; + for (int iPoint = 0; iPoint < graph->GetN(); ++iPoint) { + double x, y; + graph->GetPoint(iPoint, x, y); + if (x < xMin) + xMin = x; + if (x > xMax) + xMax = x; + } + + // fit log-likelihood function within ~2 sigma interval within maximum + // with parabola + std::vector graphPoints_forFit; + double xMin_fit = 1.e+6; + double xMax_fit = 0.; + for (int iPoint = 0; iPoint < graph->GetN(); ++iPoint) { + double x, y; + graph->GetPoint(iPoint, x, y); + double xErr = graph->GetErrorX(iPoint); + double yErr = graph->GetErrorY(iPoint); + // std::cout << "point #" << iPoint << ": x = " << x << " +/- " << xErr + // << ", y = " << y << " +/- " << yErr << std::endl; + if (y > (1.e-1 * y_Lmax) && TMath::Abs(iPoint - idxPoint_Lmax) <= 5) { + GraphPoint graphPoint; + graphPoint.x_ = x; + graphPoint.xErr_ = xErr; + if ((x - xErr) < xMin_fit) + xMin_fit = x - xErr; + if ((x + xErr) > xMax_fit) + xMax_fit = x + xErr; + graphPoint.y_ = -TMath::Log(y); + graphPoint.yErr_ = yErr / y; + graphPoints_forFit.push_back(graphPoint); + } + } + + TGraphErrors *likelihoodGraph_forFit = + makeGraph("svFitLikelihoodGraph_forFit", graphPoints_forFit); + int numPoints = likelihoodGraph_forFit->GetN(); + bool useFit = false; + if (numPoints >= 3) { + TF1 *fitFunction = + new TF1("fitFunction", "TMath::Power((x - [0])/[1], 2.) + [2]", + xMin_fit, xMax_fit); + fitFunction->SetParameter(0, x_Lmax); + fitFunction->SetParameter(1, 0.20 * x_Lmax); + fitFunction->SetParameter(2, -TMath::Log(y_Lmax)); + + std::string fitOptions = "NSQ"; + // if ( !verbosity ) fitOptions.append("Q"); + TFitResultPtr fitResult = + likelihoodGraph_forFit->Fit(fitFunction, fitOptions.data()); + if (fitResult.Get()) { + if (verbosity >= 1) { + std::cout << "fitting graph of p versus M(test) in range " + << xMin_fit << ".." << xMax_fit + << ", result:" << std::endl; + std::cout << " parameter #0 = " << fitFunction->GetParameter(0) + << " +/- " << fitFunction->GetParError(0) + << std::endl; + std::cout << " parameter #1 = " << fitFunction->GetParameter(1) + << " +/- " << fitFunction->GetParError(1) + << std::endl; + std::cout << " parameter #2 = " << fitFunction->GetParameter(2) + << " +/- " << fitFunction->GetParError(2) + << std::endl; + std::cout << "chi^2 = " << fitResult->Chi2() << std::endl; + } + if (fitResult->Chi2() < (10. * numPoints) && + fitFunction->GetParameter(0) > xMin && + fitFunction->GetParameter(0) < xMax && + TMath::Abs(fitFunction->GetParameter(0) - x_Lmax) < + (0.10 * x_Lmax)) { + mass = fitFunction->GetParameter(0); + massErr = TMath::Sqrt(square(fitFunction->GetParameter(1)) + + square(fitFunction->GetParError(0))); + Lmax = TMath::Exp(-fitFunction->GetParameter(2)); + // std::cout << "fit: mass = " << mass << " +/- " << massErr << + // " (Lmax = " << Lmax << ")" << std::endl; + useFit = true; + } + } else { + std::cerr << "Warning in : Fit did not converge !!" + << std::endl; + } + delete fitFunction; + } + if (!useFit) { + mass = x_Lmax; + massErr = TMath::Sqrt(0.5 * (square(x_Lmax - xMin_fit) + + square(xMax_fit - x_Lmax))) / + TMath::Sqrt(2. * TMath::Log(10.)); + Lmax = y_Lmax; + // std::cout << "graph: mass = " << mass << " +/- " << massErr << " + // (Lmax = " << Lmax << ")" << std::endl; + } + + delete likelihoodGraph_forFit; +} + +Vector normalize(const Vector &p) { + double p_x = p.x(); + double p_y = p.y(); + double p_z = p.z(); + double mag2 = square(p_x) + square(p_y) + square(p_z); + if (mag2 <= 0.) + return p; + double mag = TMath::Sqrt(mag2); + return Vector(p_x / mag, p_y / mag, p_z / mag); +} + +double compScalarProduct(const Vector &p1, const Vector &p2) { + return (p1.x() * p2.x() + p1.y() * p2.y() + p1.z() * p2.z()); +} + +Vector compCrossProduct(const Vector &p1, const Vector &p2) { + double p3_x = p1.y() * p2.z() - p1.z() * p2.y(); + double p3_y = p1.z() * p2.x() - p1.x() * p2.z(); + double p3_z = p1.x() * p2.y() - p1.y() * p2.x(); + return Vector(p3_x, p3_y, p3_z); +} + +double compCosThetaNuNu(double visEn, double visP, double visMass2, + double nunuEn, double nunuP, double nunuMass2) { + double cosThetaNuNu = + (visEn * nunuEn - 0.5 * (tauLeptonMass2 - (visMass2 + nunuMass2))) / + (visP * nunuP); + return cosThetaNuNu; +} + +double compPSfactor_tauToLepDecay(double x, double visEn, double visP, + double visMass, double nunuEn, double nunuP, + double nunuMass) { + // std::cout << ":" << std::endl; + // std::cout << " x = " << x << std::endl; + // std::cout << " visEn = " << visEn << std::endl; + // std::cout << " visP = " << visP << std::endl; + // std::cout << " visMass = " << visMass << std::endl; + // std::cout << " nunuEn = " << nunuEn << std::endl; + // std::cout << " nunuP = " << nunuP << std::endl; + // std::cout << " nunuMass = " << nunuMass << std::endl; + double visMass2 = square(visMass); + double nunuMass2 = square(nunuMass); + if (x >= (visMass2 / tauLeptonMass2) && x <= 1. && + nunuMass2 < ((1. - x) * tauLeptonMass2)) { // physical solution + double tauEn_rf = + (tauLeptonMass2 + nunuMass2 - visMass2) / (2. * nunuMass); + double visEn_rf = tauEn_rf - nunuMass; + if (!(tauEn_rf >= tauLeptonMass && visEn_rf >= visMass)) + return 0.; + double I = + nunuMass2 * + (2. * tauEn_rf * visEn_rf - + (2. / 3.) * TMath::Sqrt((square(tauEn_rf) - tauLeptonMass2) * + (square(visEn_rf) - visMass2))); +#ifdef XSECTION_NORMALIZATION + I *= GFfactor; +#endif + double cosThetaNuNu = + compCosThetaNuNu(visEn, visP, visMass2, nunuEn, nunuP, nunuMass2); + if (!(cosThetaNuNu >= (-1. + epsilon) && cosThetaNuNu <= +1.)) + return 0.; + double PSfactor = + (visEn + nunuEn) * I / + (8. * visP * square(x) * + TMath::Sqrt(square(visP) + square(nunuP) + + 2. * visP * nunuP * cosThetaNuNu + tauLeptonMass2)); +//------------------------------------------------------------------------- +// CV: fudge factor to reproduce literature value for cross-section times +// branching fraction +#ifdef XSECTION_NORMALIZATION + PSfactor *= 2.; +#endif + //------------------------------------------------------------------------- + return PSfactor; + } else { + return 0.; + } +} + +double compPSfactor_tauToHadDecay(double x, double visEn, double visP, + double visMass, double nuEn, double nuP) { + // std::cout << ":" << std::endl; + // std::cout << " x = " << x << std::endl; + // std::cout << " visEn = " << visEn << std::endl; + // std::cout << " visP = " << visP << std::endl; + // std::cout << " visMass = " << visMass << std::endl; + // std::cout << " nuEn = " << nuEn << std::endl; + // std::cout << " nuP = " << nuP << std::endl; + double visMass2 = square(visMass); + if (x >= (visMass2 / tauLeptonMass2) && x <= 1.) { // physical solution + double cosThetaNu = + compCosThetaNuNu(visEn, visP, visMass2, nuEn, nuP, 0.); + // std::cout << "cosThetaNu = " << cosThetaNu << std::endl; + if (!(cosThetaNu >= (-1. + epsilon) && cosThetaNu <= +1.)) + return 0.; + double PSfactor = + (visEn + nuEn) / + (8. * visP * square(x) * + TMath::Sqrt(square(visP) + square(nuP) + + 2. * visP * nuP * cosThetaNu + tauLeptonMass2)); + PSfactor *= 1.0 / (tauLeptonMass2 - visMass2); +//------------------------------------------------------------------------- +// CV: multiply by constant matrix element, +// chosen such that the branching fraction of the tau to decay into hadrons +// is reproduced +// const double M2 +// = 16.*TMath::Pi()*cube(tauLeptonMass)*GammaTauToHad/(tauLeptonMass2 - +// visMass2); Remove multiplication as it add to execution time, and does not +// alter the result. +#ifdef XSECTION_NORMALIZATION + PSfactor *= M2; +#endif + //------------------------------------------------------------------------- + return PSfactor; + } else { + return 0.; + } +} + +void integrationParameters::reset() { + idx_X_ = -1; + idx_phi_ = -1; + idx_VisPtShift_ = -1; + idx_mNuNu_ = -1; +} + +} // namespace fastmtt diff --git a/src/SVFit/FastMTT.cxx b/src/SVFit/FastMTT.cxx new file mode 100644 index 00000000..02d5f8a4 --- /dev/null +++ b/src/SVFit/FastMTT.cxx @@ -0,0 +1,534 @@ +#include + +#include "Math/Factory.h" +#include "Math/Functor.h" +#include "Math/Minimizer.h" + +#include "TMath.h" +#include "TMatrixD.h" +#include "TVector2.h" + +#include "Math/BasicMinimizer.h" +#include "TF1.h" + +#include "../../include/SVFit/FastMTT.hxx" +#include "../../include/SVFit/MeasuredTauLepton.hxx" + +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +Likelihood::Likelihood() { + + covMET.ResizeTo(2, 2); + + compnentsBitWord.reset(); + enableComponent(fastMTT::MASS); + enableComponent(fastMTT::MET); + /// experimental components. Disabled by default. + disableComponent(fastMTT::PX); + disableComponent(fastMTT::PY); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +Likelihood::~Likelihood() {} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +template +std::array getPowTable(double x) { + std::array powerTable{}; + powerTable[0] = 1.; + for (unsigned int i = 1; i <= max_order; ++i) { + powerTable[i] = powerTable[i - 1] * x; + } + return powerTable; +}; +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void Likelihood::setLeptonInputs(const LorentzVector &aLeg1P4, + const LorentzVector &aLeg2P4, + int aLeg1DecayType, int aLeg2DecayType, + int aLeg1DecayMode, int aLeg2DecayMode) { + leg1P4 = aLeg1P4; + leg2P4 = aLeg2P4; + using PowTable = std::array; + auto getPowTableLam = [](double x) { + PowTable powerTable{}; + powerTable[0] = 1.; + for (unsigned int i = 1; i <= 4; ++i) { + powerTable[i] = powerTable[i - 1] * x; + } + return powerTable; + }; + allpTpows = std::array, 2u>{ + {{{getPowTableLam(leg1P4.Px()), getPowTableLam(leg1P4.Py()), + getPowTableLam(leg1P4.Pz())}}, + {{getPowTableLam(leg2P4.Px()), getPowTableLam(leg2P4.Py()), + getPowTableLam(leg2P4.Pz())}}}}; + + mVis = (leg1P4 + leg2P4).M(); + + mVisLeg1 = leg1P4.M(); + mVisLeg2 = leg2P4.M(); + + if (aLeg1DecayType == fastmtt::MeasuredTauLepton::kTauToHadDecay && + mVisLeg1 > 1.5) { + mVisLeg1 = 0.3; + } + if (aLeg2DecayType == fastmtt::MeasuredTauLepton::kTauToHadDecay && + mVisLeg2 > 1.5) { + mVisLeg2 = 0.3; + } + + const double &mTau = fastmtt::tauLeptonMass; + mVis1OverTauSquare = std::pow(mVisLeg1 / mTau, 2); + mVis2OverTauSquare = std::pow(mVisLeg2 / mTau, 2); + + leg1DecayType = aLeg1DecayType; + leg2DecayType = aLeg2DecayType; + + leg1DecayMode = aLeg1DecayMode; + leg2DecayMode = aLeg2DecayMode; +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void Likelihood::setMETInputs(const LorentzVector &aMET, + const TMatrixD &aCovMET) { + recoMET = aMET; + covMET = aCovMET; +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void Likelihood::setParameters(const std::vector &aPars) { + + parameters = aPars; + if (parameters.size() < 2) + parameters = std::vector{6, 1.0 / 1.15}; +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void Likelihood::enableComponent(fastMTT::likelihoodComponent aCompIndex) { + + compnentsBitWord.set(aCompIndex); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void Likelihood::disableComponent(fastMTT::likelihoodComponent aCompIndex) { + + compnentsBitWord.reset(aCompIndex); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +double Likelihood::massLikelihood(const double &m) const { + + double coeff1 = parameters[0]; + double coeff2 = parameters[1]; + double mScaled = m * coeff2; + + if (mScaled < mVis) + return 0.0; + + const double mVS2 = std::pow(mVis / mScaled, 2); + + double x1Min = std::min(1.0, mVis1OverTauSquare); + double x2Min = std::max(mVis2OverTauSquare, mVS2); + double x2Max = std::min(1.0, mVS2 / x1Min); + + if (x2Max < x2Min) + return 0.0; + + double jacobiFactor = 2.0 * std::pow(mVis, 2) * std::pow(mScaled, -coeff1); + double x2IntegralTerm = log(x2Max) - log(x2Min); + + double value = x2IntegralTerm; + if (leg1DecayType != fastmtt::MeasuredTauLepton::kTauToHadDecay) { + + double mNuNuIntegralTermLeg1 = + mVS2 * (std::pow(x2Max, -1) - std::pow(x2Min, -1)); + value += mNuNuIntegralTermLeg1; + } + if (leg2DecayType != fastmtt::MeasuredTauLepton::kTauToHadDecay) { + double mNuNuIntegralTermLeg2 = mVS2 * x2IntegralTerm - (x2Max - x2Min); + value += mNuNuIntegralTermLeg2; + } + + /// The E9 factor to get values around 1.0 + value *= 1E9 * jacobiFactor; + + return value; +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +double Likelihood::ptLikelihood(const double &pTTauTau, int type) const { + + /// Protection against numerical singularity in phase space volume. + if (std::abs(pTTauTau) < 0.5) + return 0.0; + + const auto pT1pow = allpTpows[0][type]; + const auto pT2pow = allpTpows[1][type]; + const auto pT1 = pT1pow[1]; + const auto pT2 = pT2pow[1]; + const auto pTTauTauPow = getPowTable<5>(pTTauTau); + + Double_t x1Min = std::min(1.0, mVis1OverTauSquare); + Double_t x2Min = std::min(1.0, mVis2OverTauSquare); + + Double_t x2Max = 1.0; + Double_t x1Max = 1.0; + + Double_t a_x2 = x1Min * pT2 / (x1Min * pTTauTau - pT1); + Double_t b_x2 = x1Max * pT2 / (x1Max * pTTauTau - pT1); + + const bool is_x1_vs_x2_falling = (-pT2 * pT1) < 0; + const double x1_singularity = pT1 / pTTauTau; + const bool x2_vs_x1_hasSingularity = + x1_singularity > 0.0 && x1_singularity < 1.0; + if (x2_vs_x1_hasSingularity && x1_singularity < x1Min) + return 0.0; + + if (is_x1_vs_x2_falling) { + x2Min = std::max(x2Min, b_x2); + x2Max = std::min(x2Max, a_x2); + if (x2_vs_x1_hasSingularity && x2Max < 0) + x2Max = 1.0; + } else { + x2Min = std::max(x2Min, a_x2); + x2Max = std::min(x2Max, b_x2); + if (x2_vs_x1_hasSingularity && x2Max < 0) + x2Max = 1.0; + } + + if (x2Min < 0) + x2Min = 0.0; + if (x2Min > x2Max) + return 0.0; + + Double_t mNuNuIntegral = 0.0; + Double_t x2 = std::min(1.0, x2Max); + + const Double_t term1 = pT2 - pTTauTau * x2; + const Double_t log_term1 = log(std::abs(term1)); + const Double_t term1Square = std::pow(term1, 2); + + Double_t integralMax = + (pT1 * (pTTauTau * x2 + pT2pow[2] / term1 + 2 * pT2 * log_term1)) / + pTTauTauPow[3]; + if (leg1DecayType != fastmtt::MeasuredTauLepton::kTauToHadDecay) { + mNuNuIntegral = + -pT1pow[2] * + (2 * pTTauTau * x2 + + (pT2pow[2] * (5 * pT2 - 6 * pTTauTau * x2)) / term1Square + + 6 * pT2 * log_term1) / + (2 * pTTauTauPow[4]); + } + if (leg2DecayType != fastmtt::MeasuredTauLepton::kTauToHadDecay) { + mNuNuIntegral += -pT1 / (2 * pTTauTauPow[5]) * + (2 * pT2 * pTTauTau * (-3 * pT1 + 2 * pTTauTau) * x2 + + pTTauTauPow[2] * (-pT1 + pTTauTau) * pow(x2, 2) + + (pT2pow[4] * pT1) / term1Square + + (2 * pT2pow[3] * (-4 * pT1 + pTTauTau)) / term1 + + 6 * pT2pow[2] * (-2 * pT1 + pTTauTau) * log_term1); + } + integralMax += mNuNuIntegral; + + x2 = x2Min; + const Double_t term2 = pT2 - pTTauTau * x2; + const Double_t log_term2 = log(std::abs(term2)); + const Double_t term2Square = std::pow(term2, 2); + + Double_t integralMin = + (pT1 * (pTTauTau * x2 + pT2pow[2] / term2 + 2 * pT2 * log_term2)) / + pTTauTauPow[3]; + if (leg1DecayType != fastmtt::MeasuredTauLepton::kTauToHadDecay) { + mNuNuIntegral = + -pT1pow[2] * + (2 * pTTauTau * x2 + + (pT2pow[2] * (5 * pT2 - 6 * pTTauTau * x2)) / term2Square + + 6 * pT2 * log_term2) / + (2 * pTTauTauPow[4]); + } + if (leg2DecayType != fastmtt::MeasuredTauLepton::kTauToHadDecay) { + mNuNuIntegral += -pT1 / (2 * pTTauTauPow[5]) * + (2 * pT2 * pTTauTau * (-3 * pT1 + 2 * pTTauTau) * x2 + + pTTauTauPow[2] * (-pT1 + pTTauTau) * pow(x2, 2) + + (pT2pow[4] * pT1) / term2Square + + (2 * pT2pow[3] * (-4 * pT1 + pTTauTau)) / term2 + + 6 * pT2pow[2] * (-2 * pT1 + pTTauTau) * log_term2); + } + + integralMin += mNuNuIntegral; + + Double_t value = integralMax - integralMin; + + /// The 1E4 factor to get values around 1.0 + value *= 1E4; + + return std::abs(value); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +double Likelihood::metTF(const LorentzVector &metP4, const LorentzVector &nuP4, + const TMatrixD &covMET) const { + + const double aMETx = metP4.X(); + const double aMETy = metP4.Y(); + + double invCovMETxx = covMET(1, 1); + double invCovMETxy = -covMET(0, 1); + double invCovMETyx = -covMET(1, 0); + double invCovMETyy = covMET(0, 0); + double covDet = invCovMETxx * invCovMETyy - invCovMETxy * invCovMETyx; + + if (std::abs(covDet) < 1E-10) { + std::cerr << "Error: Cannot invert MET covariance Matrix (det=0) !!" + << "METx: " << aMETy << " METy: " << aMETy << std::endl; + return 0; + } + double const_MET = 1. / (2. * M_PI * TMath::Sqrt(covDet)); + double residualX = aMETx - (nuP4.X()); + double residualY = aMETy - (nuP4.Y()); + + double pull2 = + residualX * (invCovMETxx * residualX + invCovMETxy * residualY) + + residualY * (invCovMETyx * residualX + invCovMETyy * residualY); + pull2 /= covDet; + + return const_MET * TMath::Exp(-0.5 * pull2); +} +////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////// +double Likelihood::value(const double *x) const { + + double x1Min = std::min(1.0, mVis1OverTauSquare); + double x2Min = std::min(1.0, mVis2OverTauSquare); + + if (x[0] < x1Min || x[1] < x2Min) + return 0.0; + + const auto testP4 = leg1P4 * (1.0 / x[0]) + leg2P4 * (1.0 / x[1]); + const auto testMET = testP4 - leg1P4 - leg2P4; + + double value = -1.0; + if (compnentsBitWord.test(fastMTT::MET)) + value *= metTF(recoMET, testMET, covMET); + if (compnentsBitWord.test(fastMTT::MASS)) + value *= massLikelihood(testP4.M()); + if (compnentsBitWord.test(fastMTT::PX)) + value *= ptLikelihood(testP4.Px(), 0); + if (compnentsBitWord.test(fastMTT::PY)) + value *= ptLikelihood(testP4.Py(), 1); + return value; +} +////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////// +FastMTT::FastMTT() { + + minimizerName = "Minuit2"; + minimizerAlgorithm = "Migrad"; + initialize(); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +FastMTT::~FastMTT() { delete minimizer; } +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void FastMTT::initialize() { + + minimizer = + ROOT::Math::Factory::CreateMinimizer(minimizerName, minimizerAlgorithm); + minimizer->SetMaxFunctionCalls(100000); + minimizer->SetMaxIterations(100000); + minimizer->SetTolerance(0.01); + + std::vector varNames = {"x1", "x2"}; + nVariables = varNames.size(); + std::vector initialValues(nVariables, 0.5); + std::vector stepSizes(nVariables, 0.01); + minimumPosition = initialValues; + minimumValue = 999.0; + + for (unsigned int iVar = 0; iVar < nVariables; ++iVar) { + minimizer->SetVariable(iVar, varNames[iVar].c_str(), + initialValues[iVar], stepSizes[iVar]); + } + + std::vector shapeParams = {6, 1.0 / 1.15}; + setLikelihoodParams(shapeParams); + likelihoodFunctor = + new ROOT::Math::Functor(&myLikelihood, &Likelihood::value, nVariables); + minimizer->SetFunction(*likelihoodFunctor); + + verbosity = 0; +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void FastMTT::setLikelihoodParams(const std::vector &aPars) { + + myLikelihood.setParameters(aPars); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void FastMTT::enableComponent(fastMTT::likelihoodComponent aCompIndex) { + + myLikelihood.enableComponent(aCompIndex); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void FastMTT::disableComponent(fastMTT::likelihoodComponent aCompIndex) { + + myLikelihood.disableComponent(aCompIndex); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +bool FastMTT::compareLeptons( + const fastmtt::MeasuredTauLepton &measuredTauLepton1, + const fastmtt::MeasuredTauLepton &measuredTauLepton2) { + + using namespace fastmtt; + + if ((measuredTauLepton1.type() == MeasuredTauLepton::kTauToElecDecay || + measuredTauLepton1.type() == MeasuredTauLepton::kTauToMuDecay) && + measuredTauLepton2.type() == MeasuredTauLepton::kTauToHadDecay) + return true; + if ((measuredTauLepton2.type() == MeasuredTauLepton::kTauToElecDecay || + measuredTauLepton2.type() == MeasuredTauLepton::kTauToMuDecay) && + measuredTauLepton1.type() == MeasuredTauLepton::kTauToHadDecay) + return false; + return (measuredTauLepton1.pt() > measuredTauLepton2.pt()); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void FastMTT::run( + const std::vector &measuredTauLeptons, + const double &measuredMETx, const double &measuredMETy, + const TMatrixD &covMET) { + + bestP4 = LorentzVector(); + //////////////////////////////////////////// + + if (measuredTauLeptons.size() != 2) { + std::cout << "Number of MeasuredTauLepton is " + << measuredTauLeptons.size() + << " a user shouls pass exactly two leptons." << std::endl; + return; + } + + std::vector sortedMeasuredTauLeptons = + measuredTauLeptons; + std::sort(sortedMeasuredTauLeptons.begin(), sortedMeasuredTauLeptons.end(), + compareLeptons); + + double metLength = + sqrt(std::pow(measuredMETx, 2) + std::pow(measuredMETy, 2)); + LorentzVector aMET = + LorentzVector(measuredMETx, measuredMETy, 0, metLength); + + const fastmtt::MeasuredTauLepton &aLepton1 = measuredTauLeptons[0]; + const fastmtt::MeasuredTauLepton &aLepton2 = measuredTauLeptons[1]; + + myLikelihood.setLeptonInputs(aLepton1.p4(), aLepton2.p4(), aLepton1.type(), + aLepton2.type(), aLepton1.decayMode(), + aLepton2.decayMode()); + myLikelihood.setMETInputs(aMET, covMET); + + scan(); + // minimize(); + + tau1P4 = aLepton1.p4() * (1.0 / minimumPosition[0]); + tau2P4 = aLepton2.p4() * (1.0 / minimumPosition[1]); + bestP4 = tau1P4 + tau2P4; +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void FastMTT::minimize() { + + clock.Reset(); + clock.Start("minimize"); + + minimizer->SetVariableLimits(0, 0.01, 1.0); + minimizer->SetVariableLimits(1, 0.01, 1.0); + + minimizer->SetVariable(0, "x1", 0.5, 0.1); + minimizer->SetVariable(1, "x2", 0.5, 0.1); + + minimizer->Minimize(); + + const double *theMinimum = minimizer->X(); + minimumPosition[0] = theMinimum[0]; + minimumPosition[1] = theMinimum[1]; + minimumValue = minimizer->MinValue(); + + if (true || minimizer->Status() != 0) { + std::cout << " minimizer " + << " Status: " << minimizer->Status() + << " nCalls: " << minimizer->NCalls() + << " nIterations: " << minimizer->NIterations() + << " x1Max: " << theMinimum[0] << " x2Max: " << theMinimum[1] + << " max LLH: " << minimizer->MinValue() + << " m: " << bestP4.M() << std::endl; + } + clock.Stop("minimize"); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +void FastMTT::scan() { + + clock.Reset(); + clock.Start("scan"); + + double lh = 0.0; + double bestLH = 0.0; + + double x[2] = {0.5, 0.5}; + + double theMinimum[2] = {0.75, 0.75}; + const int nGridPoints = 100; + const double gridFactor = 1. / nGridPoints; + int nCalls = 0; + for (int iX2 = 1; iX2 < nGridPoints; ++iX2) { + x[1] = iX2 * gridFactor; + for (int iX1 = 1; iX1 < nGridPoints; ++iX1) { + x[0] = iX1 * gridFactor; + lh = myLikelihood.value(x); + ++nCalls; + if (lh < bestLH) { + bestLH = lh; + theMinimum[0] = x[0]; + theMinimum[1] = x[1]; + } + } + + minimumPosition[0] = theMinimum[0]; + minimumPosition[1] = theMinimum[1]; + minimumValue = bestLH; + } + clock.Stop("scan"); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +double FastMTT::getCpuTime(const std::string &method) { + + return clock.GetCpuTime(method.c_str()); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +double FastMTT::getRealTime(const std::string &method) { + + return clock.GetRealTime(method.c_str()); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +std::tuple FastMTT::getBestX() const { + + return std::make_tuple(minimumPosition[0], minimumPosition[1]); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +double FastMTT::getBestLikelihood() const { return minimumValue; } +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// +double FastMTT::getLikelihoodForX(double *x) const { + + return myLikelihood.value(x); +} +/////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////// diff --git a/src/SVFit/MeasuredTauLepton.cxx b/src/SVFit/MeasuredTauLepton.cxx new file mode 100644 index 00000000..d2de8cfc --- /dev/null +++ b/src/SVFit/MeasuredTauLepton.cxx @@ -0,0 +1,104 @@ +#include "../../include/SVFit/MeasuredTauLepton.hxx" +#include "../../include/SVFit/AuxFunctions.hxx" + +#include + +using namespace fastmtt; + +MeasuredTauLepton::MeasuredTauLepton() + : type_(kUndefinedDecayType), pt_(0.), eta_(0.), phi_(0.), mass_(0.), + decayMode_(-1) { + initialize(); +} + +MeasuredTauLepton::MeasuredTauLepton(int type, double pt, double eta, + double phi, double mass, int decayMode) + : type_(type), pt_(pt), eta_(eta), phi_(phi), mass_(mass), + decayMode_(decayMode) { + double minVisMass = electronMass; + double maxVisMass = tauLeptonMass; + if (type_ == kTauToElecDecay) { + minVisMass = electronMass; + maxVisMass = minVisMass; + } else if (type_ == kTauToMuDecay) { + minVisMass = muonMass; + maxVisMass = minVisMass; + } else if (type_ == kTauToHadDecay) { + if (decayMode_ == -1) { + minVisMass = chargedPionMass; + maxVisMass = 1.5; + } else if (decayMode_ == 0) { + minVisMass = chargedPionMass; + maxVisMass = minVisMass; + } else { + minVisMass = 0.3; + maxVisMass = 1.5; + } + } + preciseVisMass_ = mass_; + if (preciseVisMass_ < minVisMass) + preciseVisMass_ = minVisMass; + if (preciseVisMass_ > maxVisMass) + preciseVisMass_ = maxVisMass; + initialize(); +} + +MeasuredTauLepton::MeasuredTauLepton(const MeasuredTauLepton &measuredTauLepton) + : type_(measuredTauLepton.type()), pt_(measuredTauLepton.pt()), + eta_(measuredTauLepton.eta()), phi_(measuredTauLepton.phi()), + mass_(measuredTauLepton.mass()), + decayMode_(measuredTauLepton.decayMode()) { + preciseVisMass_ = measuredTauLepton.mass(); + + initialize(); +} + +MeasuredTauLepton::~MeasuredTauLepton() {} + +int MeasuredTauLepton::type() const { return type_; } + +double MeasuredTauLepton::pt() const { return pt_; } +double MeasuredTauLepton::eta() const { return eta_; } +double MeasuredTauLepton::phi() const { return phi_; } +double MeasuredTauLepton::mass() const { return preciseVisMass_; } + +double MeasuredTauLepton::energy() const { return energy_; } +double MeasuredTauLepton::px() const { return px_; } +double MeasuredTauLepton::py() const { return py_; } +double MeasuredTauLepton::pz() const { return pz_; } + +double MeasuredTauLepton::p() const { return p_; } + +int MeasuredTauLepton::decayMode() const { return decayMode_; } + +LorentzVector MeasuredTauLepton::p4() const { return p4_; } + +Vector MeasuredTauLepton::p3() const { return p3_; } + +double MeasuredTauLepton::cosPhi_sinTheta() const { return cosPhi_sinTheta_; } +double MeasuredTauLepton::sinPhi_sinTheta() const { return sinPhi_sinTheta_; } +double MeasuredTauLepton::cosTheta() const { return cosTheta_; } + +void MeasuredTauLepton::roundToNdigits(unsigned int nDigis) { + pt_ = fastmtt::roundToNdigits(pt_, nDigis); + eta_ = fastmtt::roundToNdigits(eta_, nDigis); + phi_ = fastmtt::roundToNdigits(phi_, nDigis); + mass_ = fastmtt::roundToNdigits(mass_, nDigis); + initialize(); +} + +void MeasuredTauLepton::initialize() { + // CV: relations between pT and p, energy taken from + // http://en.wikipedia.org/wiki/Pseudorapidity + p_ = pt_ * TMath::CosH(eta_); + px_ = pt_ * TMath::Cos(phi_); + py_ = pt_ * TMath::Sin(phi_); + pz_ = pt_ * TMath::SinH(eta_); + energy_ = TMath::Sqrt(p_ * p_ + preciseVisMass_ * preciseVisMass_); + p4_ = LorentzVector(px_, py_, pz_, energy_); + p3_ = Vector(px_, py_, pz_); + double theta = p4_.theta(); + cosPhi_sinTheta_ = TMath::Cos(phi_) * TMath::Sin(theta); + sinPhi_sinTheta_ = TMath::Sin(phi_) * TMath::Sin(theta); + cosTheta_ = TMath::Cos(theta); +} diff --git a/src/quantities.cxx b/src/quantities.cxx index 7c238c1a..7cd0d9bb 100644 --- a/src/quantities.cxx +++ b/src/quantities.cxx @@ -1,6 +1,8 @@ #ifndef GUARD_QUANTITIES_H #define GUARD_QUANTITIES_H +#include "../include/SVFit/FastMTT.hxx" +#include "../include/SVFit/MeasuredTauLepton.hxx" #include "../include/basefunctions.hxx" #include "../include/defaults.hxx" #include "../include/utility/Logger.hxx" @@ -176,6 +178,107 @@ ROOT::RDF::RNode m_vis(ROOT::RDF::RNode df, const std::string &outputname, }, inputvectors); } + +/** + * @brief Function used to calculate the FastMTT p4 from the given inputs. The + * implementation is based on + * https://github.com/SVfit/ClassicSVfit/tree/fastMTT_19_02_2019 + * + * @param df The dataframe to add the quantity to + * @param outputname name of the new column containing the lorentz vector value + * @param pt_1 the name of the column containing the pt of the first particle + * @param pt_2 the name of the column containing the pt of the second particle + * @param eta_1 the name of the column containing the eta of the first particle + * @param eta_2 the name of the column containing the eta of the second particle + * @param phi_1 the name of the column containing the phi of the first particle + * @param phi_2 the name of the column containing the phi of the second particle + * @param mass_1 the name of the column containing the mass of the first + * particle + * @param mass_2 the name of the column containing the mass of the second + * particle + * @param met_pt the name of the column containing the met pt + * @param met_phi the name of the column containing the met phi + * @param met_cov_xx the name of the column containing the met covariance xx + * @param met_cov_xy the name of the column containing the met covariance xy + * @param met_cov_yy the name of the column containing the met covariance yy + * @param decay_mode_1 the name of the column containing the decay mode of the + * first particle + * @param decay_mode_2 the name of the column containing the decay mode of the + * second particle + * @param finalstate the final state of the ditaudecay. Supported are "mt", + * "et", "tt", "em" + * @return ROOT::RDF::RNode + */ +ROOT::RDF::RNode +p4_fastmtt(ROOT::RDF::RNode df, const std::string &outputname, + const std::string &pt_1, const std::string &pt_2, + const std::string &eta_1, const std::string &eta_2, + const std::string &phi_1, const std::string &phi_2, + const std::string &mass_1, const std::string &mass_2, + const std::string &met_pt, const std::string &met_phi, + const std::string &met_cov_xx, const std::string &met_cov_xy, + const std::string &met_cov_yy, const std::string &decay_mode_1, + const std::string &decay_mode_2, const std::string &finalstate) { + auto calculate_fast_mtt = + [finalstate](const float &pt_1, const float &pt_2, const float &eta_1, + const float &eta_2, const float &phi_1, const float &phi_2, + const float &mass_1, const float &mass_2, + const float &met_pt, const float &met_phi, + const float &met_cov_xx, const float &met_cov_xy, + const float &met_cov_yy, const int &decay_mode_1, + const int &decay_mode_2) { + std::vector measuredTauLeptons; + TMatrixD covMET(2, 2); + covMET[0][0] = met_cov_xx; + covMET[1][0] = met_cov_xy; + covMET[0][1] = met_cov_xy; + covMET[1][1] = met_cov_yy; + // build the met lorentz vector + ROOT::Math::PtEtaPhiMVector met(met_pt, 0.0, met_phi, 0.0); + // set the decay modes according to the final state + auto decay_obj_1 = fastmtt::MeasuredTauLepton::kTauToHadDecay; + auto decay_obj_2 = fastmtt::MeasuredTauLepton::kTauToHadDecay; + int dm_1, dm_2; + if (finalstate == "mt") { + dm_1 = -1; + dm_2 = decay_mode_2; + auto decay_obj_1 = fastmtt::MeasuredTauLepton::kTauToMuDecay; + auto decay_obj_2 = fastmtt::MeasuredTauLepton::kTauToHadDecay; + } else if (finalstate == "et") { + dm_1 = -1; + dm_2 = decay_mode_2; + auto decay_obj_1 = fastmtt::MeasuredTauLepton::kTauToElecDecay; + auto decay_obj_2 = fastmtt::MeasuredTauLepton::kTauToHadDecay; + } else if (finalstate == "tt") { + dm_1 = decay_mode_1; + dm_2 = decay_mode_2; + } else if (finalstate == "em") { + dm_1 = -1; + dm_2 = -1; + auto decay_obj_1 = fastmtt::MeasuredTauLepton::kTauToElecDecay; + auto decay_obj_2 = fastmtt::MeasuredTauLepton::kTauToMuDecay; + } else { + Logger::get("FastMTT")->error( + "Final state {} not supported by FastMTT", finalstate); + return (ROOT::Math::PtEtaPhiMVector)LorentzVector(); + } + measuredTauLeptons.push_back(fastmtt::MeasuredTauLepton( + decay_obj_1, pt_1, eta_1, phi_1, mass_1, dm_1)); + measuredTauLeptons.push_back(fastmtt::MeasuredTauLepton( + decay_obj_2, pt_2, eta_2, phi_2, mass_2, dm_2)); + FastMTT FastMTTAlgo; + FastMTTAlgo.run(measuredTauLeptons, met.X(), met.Y(), covMET); + LorentzVector result = FastMTTAlgo.getBestP4(); + // ROOT::Math::PtEtaPhiMVector result(_result.Pt(), _result.Eta(), + // _result.Phi(), _result.M()); + Logger::get("FastMTT")->debug("FastMTT result: {}", result.M()); + return (ROOT::Math::PtEtaPhiMVector)result; + }; + return df.Define(outputname, calculate_fast_mtt, + {pt_1, pt_2, eta_1, eta_2, phi_1, phi_2, mass_1, mass_2, + met_pt, met_phi, met_cov_xx, met_cov_xy, met_cov_yy, + decay_mode_1, decay_mode_2}); +} /// Function to calculate the visible pt from a pair of lorentz vectors and /// add it to the dataframe. The visible pt is calculated as the pt of the /// lorentz vector of the dilepton system.