add SVFit FastMTT package

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})

data/custom_top_sf/btag_eff/convert_hist_to_corrlib.py

@@ -51,7 +51,6 @@
 
 
 for y in corrs.keys():
-
     cset[y] = correctionlib.schemav2.CorrectionSet(
         schema_version=2,
         description="{} ({})".format(file_description, y),

data/custom_top_sf/electron/convert_hist_to_corrlib.py

@@ -52,7 +52,6 @@
 
 
 for y in corrs.keys():
-
     cset[y] = correctionlib.schemav2.CorrectionSet(
         schema_version=2,
         description="{} ({})".format(file_description, y),

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 <TGraphErrors.h>
+
+#include <string>
+#include <vector>
+
+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<ROOT::Math::Cartesian3D<double>>
+    Vector;
+/**
+   \typedef SVfitStandalone::LorentzVector
+   \brief   lorentz vector (equivalent to reco::Candidate::LorentzVector)
+*/
+typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>> 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<GraphPoint> &);
+
+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

include/SVFit/FastMTT.hxx

@@ -0,0 +1,204 @@
+#ifndef FastMTT_FastMTT_H
+#define FastMTT_FastMTT_H
+
+#include <bitset>
+#include <string>
+#include <tuple>
+#include <vector>
+
+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<ROOT::Math::PxPyPzE4D<double>> 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<double> &parameters);
+
+    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<double> 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<double, 5u>; // first powers of a number
+    // array with dimensions 2x3x5 used to store the powers of the coefficients
+    // of two vectors
+    std::array<std::array<PowTable, 3u>, 2u> allpTpows;
+};
+//////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////
+class FastMTT {
+
+  public:
+    FastMTT();
+
+    ~FastMTT();
+
+    void initialize();
+
+    /// Run fastMTT algorithm for given input
+    void run(const std::vector<fastmtt::MeasuredTauLepton> &, 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<double> &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<double, double> 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<double> minimumPosition;
+
+    /// Mimimum value
+    double minimumValue;
+
+    /// Dimension of minimalization space
+    unsigned int nVariables;
+
+    /// Names of variables to be minimized
+    std::vector<std::string> varNames;
+
+    /// Values of variables to be minimized
+
+    std::vector<double> variables;
+
+    /// Step sizes for each minimized variable
+    std::vector<double> stepSizes;
+
+    ROOT::Math::Functor *likelihoodFunctor;
+
+    Likelihood myLikelihood;
+
+    LorentzVector tau1P4, tau2P4, bestP4;
+
+    TBenchmark clock;
+
+    int verbosity;
+};
+
+#endif