From 9aa3f747ac4b0bad6d007f6cb0dd8653eee09286 Mon Sep 17 00:00:00 2001 From: Sebastian Brommer Date: Wed, 5 Jun 2024 15:51:19 +0200 Subject: [PATCH 1/2] fix FastMTT memory leak --- include/SVFit/FastMTT.hxx | 79 ++----------------- src/SVFit/FastMTT.cxx | 159 +++++++------------------------------- src/quantities.cxx | 11 ++- 3 files changed, 41 insertions(+), 208 deletions(-) diff --git a/include/SVFit/FastMTT.hxx b/include/SVFit/FastMTT.hxx index cfb0eff7..fbd8e0f3 100644 --- a/include/SVFit/FastMTT.hxx +++ b/include/SVFit/FastMTT.hxx @@ -4,6 +4,7 @@ #include "Math/LorentzVector.h" #include "TBenchmark.h" #include "TMatrixD.h" +#include #include #include #include @@ -72,7 +73,7 @@ class Likelihood { std::vector parameters; /// Bit word coding enabled likelihood components - std::bitset<128> compnentsBitWord; + std::bitset<128> componentsBitWord; // precomputed values used to reduce the number of redundant calculations double mVis1OverTauSquare; @@ -94,81 +95,19 @@ class 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); + ROOT::Math::PtEtaPhiMVector + run(const std::vector &, const double &, + const double &, const TMatrixD &); 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; + void scan(Likelihood &myLikelihood); /// Minimum location std::vector minimumPosition; @@ -189,14 +128,8 @@ class FastMTT { /// Step sizes for each minimized variable std::vector stepSizes; - ROOT::Math::Functor *likelihoodFunctor; - - Likelihood myLikelihood; - LorentzVector tau1P4, tau2P4, bestP4; - TBenchmark clock; - int verbosity; }; diff --git a/src/SVFit/FastMTT.cxx b/src/SVFit/FastMTT.cxx index 02d5f8a4..eb8b6638 100644 --- a/src/SVFit/FastMTT.cxx +++ b/src/SVFit/FastMTT.cxx @@ -7,12 +7,14 @@ #include "TMath.h" #include "TMatrixD.h" #include "TVector2.h" +#include #include "Math/BasicMinimizer.h" #include "TF1.h" #include "../../include/SVFit/FastMTT.hxx" #include "../../include/SVFit/MeasuredTauLepton.hxx" +#include "../../include/utility/Logger.hxx" /////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// @@ -20,12 +22,12 @@ Likelihood::Likelihood() { covMET.ResizeTo(2, 2); - compnentsBitWord.reset(); + componentsBitWord.reset(); enableComponent(fastMTT::MASS); enableComponent(fastMTT::MET); - /// experimental components. Disabled by default. - disableComponent(fastMTT::PX); - disableComponent(fastMTT::PY); + // /// experimental components. Disabled by default. + // disableComponent(fastMTT::PX); + // disableComponent(fastMTT::PY); } /////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// @@ -107,13 +109,13 @@ void Likelihood::setParameters(const std::vector &aPars) { /////////////////////////////////////////////////////////////////// void Likelihood::enableComponent(fastMTT::likelihoodComponent aCompIndex) { - compnentsBitWord.set(aCompIndex); + componentsBitWord.set(aCompIndex); } /////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// void Likelihood::disableComponent(fastMTT::likelihoodComponent aCompIndex) { - compnentsBitWord.reset(aCompIndex); + componentsBitWord.reset(aCompIndex); } /////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// @@ -308,77 +310,33 @@ double Likelihood::value(const double *x) const { const auto testMET = testP4 - leg1P4 - leg2P4; double value = -1.0; - if (compnentsBitWord.test(fastMTT::MET)) + if (componentsBitWord.test(fastMTT::MET)) value *= metTF(recoMET, testMET, covMET); - if (compnentsBitWord.test(fastMTT::MASS)) + if (componentsBitWord.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() {} /////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// -FastMTT::~FastMTT() { delete minimizer; } +FastMTT::~FastMTT() {} /////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// void FastMTT::initialize() { - minimizer = - ROOT::Math::Factory::CreateMinimizer(minimizerName, minimizerAlgorithm); - minimizer->SetMaxFunctionCalls(100000); - minimizer->SetMaxIterations(100000); - minimizer->SetTolerance(0.01); - + // Set the start parameters 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) { @@ -397,11 +355,14 @@ bool FastMTT::compareLeptons( } /////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// -void FastMTT::run( - const std::vector &measuredTauLeptons, - const double &measuredMETx, const double &measuredMETy, - const TMatrixD &covMET) { - +ROOT::Math::PtEtaPhiMVector +FastMTT::run(const std::vector &measuredTauLeptons, + const double &measuredMETx, const double &measuredMETy, + const TMatrixD &covMET) { + Likelihood myLikelihood; + std::vector shapeParams = {6, 1.0 / 1.15}; + myLikelihood.setParameters(shapeParams); + FastMTT::initialize(); bestP4 = LorentzVector(); //////////////////////////////////////////// @@ -409,7 +370,7 @@ void FastMTT::run( std::cout << "Number of MeasuredTauLepton is " << measuredTauLeptons.size() << " a user shouls pass exactly two leptons." << std::endl; - return; + return ROOT::Math::PtEtaPhiMVector(); } std::vector sortedMeasuredTauLeptons = @@ -430,51 +391,17 @@ void FastMTT::run( aLepton2.decayMode()); myLikelihood.setMETInputs(aMET, covMET); - scan(); - // minimize(); + scan(myLikelihood); tau1P4 = aLepton1.p4() * (1.0 / minimumPosition[0]); tau2P4 = aLepton2.p4() * (1.0 / minimumPosition[1]); bestP4 = tau1P4 + tau2P4; + // Logger::get("FastMTT::run")->warn("bestP4: {}", bestP4); + ROOT::Math::PtEtaPhiMVector result = (ROOT::Math::PtEtaPhiMVector)bestP4; + return result; } -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// -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"); +void FastMTT::scan(Likelihood &myLikelihood) { double lh = 0.0; double bestLH = 0.0; @@ -501,34 +428,4 @@ void FastMTT::scan() { 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); -} -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// +} \ No newline at end of file diff --git a/src/quantities.cxx b/src/quantities.cxx index 81d8c9f7..3db8abaa 100644 --- a/src/quantities.cxx +++ b/src/quantities.cxx @@ -284,6 +284,7 @@ p4_fastmtt(ROOT::RDF::RNode df, const std::string &outputname, 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) { + // initialize the FastMTT algorithm 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, @@ -332,12 +333,12 @@ p4_fastmtt(ROOT::RDF::RNode df, const std::string &outputname, 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 = + FastMTTAlgo.run(measuredTauLeptons, met.X(), met.Y(), covMET); // 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 result; }; return df.Define(outputname, calculate_fast_mtt, {pt_1, pt_2, eta_1, eta_2, phi_1, phi_2, mass_1, mass_2, @@ -492,6 +493,8 @@ ROOT::RDF::RNode deltaR(ROOT::RDF::RNode df, const std::string &outputname, const std::string &p_1_p4, const std::string &p_2_p4) { auto calculate_deltaR = [](ROOT::Math::PtEtaPhiMVector &p_1_p4, ROOT::Math::PtEtaPhiMVector &p_2_p4) { + if (p_1_p4.pt() < 0.0 || p_2_p4.pt() < 0.0) + return default_float; return (float)ROOT::Math::VectorUtil::DeltaR(p_1_p4, p_2_p4); }; return df.Define(outputname, calculate_deltaR, {p_1_p4, p_2_p4}); @@ -575,7 +578,7 @@ ROOT::RDF::RNode pt_ttjj(ROOT::RDF::RNode df, const std::string &outputname, /** * @brief function used to calculate the pt two leading jets - If the number of jets is less than 2, the quantity is set to 10 + If the number of jets is less than 2, the quantity is set to -10 * instead. * * @param df name of the dataframe From 860d2a08fe0d959cb1bc89c45fb138886c060834 Mon Sep 17 00:00:00 2001 From: Sebastian Brommer Date: Wed, 5 Jun 2024 16:55:54 +0200 Subject: [PATCH 2/2] format file --- src/SVFit/FastMTT.cxx | 77 ++++++++++++++----------------------------- 1 file changed, 24 insertions(+), 53 deletions(-) diff --git a/src/SVFit/FastMTT.cxx b/src/SVFit/FastMTT.cxx index eb8b6638..213d2d7f 100644 --- a/src/SVFit/FastMTT.cxx +++ b/src/SVFit/FastMTT.cxx @@ -16,24 +16,16 @@ #include "../../include/SVFit/MeasuredTauLepton.hxx" #include "../../include/utility/Logger.hxx" -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// -Likelihood::Likelihood() { +Likelihood::Likelihood() { covMET.ResizeTo(2, 2); - componentsBitWord.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{}; @@ -43,8 +35,7 @@ std::array getPowTable(double x) { } return powerTable; }; -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// + void Likelihood::setLeptonInputs(const LorentzVector &aLeg1P4, const LorentzVector &aLeg2P4, int aLeg1DecayType, int aLeg2DecayType, @@ -90,37 +81,33 @@ void Likelihood::setLeptonInputs(const LorentzVector &aLeg1P4, leg1DecayMode = aLeg1DecayMode; leg2DecayMode = aLeg2DecayMode; } -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// + + void Likelihood::setMETInputs(const LorentzVector &aMET, const TMatrixD &aCovMET) { recoMET = aMET; covMET = aCovMET; } -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// -void Likelihood::setParameters(const std::vector &aPars) { + +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) { + +void Likelihood::enableComponent(fastMTT::likelihoodComponent aCompIndex) { componentsBitWord.set(aCompIndex); } -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// -void Likelihood::disableComponent(fastMTT::likelihoodComponent aCompIndex) { + +void Likelihood::disableComponent(fastMTT::likelihoodComponent aCompIndex) { componentsBitWord.reset(aCompIndex); } -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// -double Likelihood::massLikelihood(const double &m) const { + +double Likelihood::massLikelihood(const double &m) const { double coeff1 = parameters[0]; double coeff2 = parameters[1]; double mScaled = m * coeff2; @@ -157,14 +144,12 @@ double Likelihood::massLikelihood(const double &m) const { 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]; @@ -266,14 +251,12 @@ double Likelihood::ptLikelihood(const double &pTTauTau, int type) const { 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); @@ -296,8 +279,7 @@ double Likelihood::metTF(const LorentzVector &metP4, const LorentzVector &nuP4, return const_MET * TMath::Exp(-0.5 * pull2); } -////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////// + double Likelihood::value(const double *x) const { double x1Min = std::min(1.0, mVis1OverTauSquare); @@ -316,14 +298,11 @@ double Likelihood::value(const double *x) const { value *= massLikelihood(testP4.M()); return value; } -////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////// + FastMTT::FastMTT() {} -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// + FastMTT::~FastMTT() {} -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// + void FastMTT::initialize() { // Set the start parameters @@ -335,14 +314,12 @@ void FastMTT::initialize() { minimumValue = 999.0; verbosity = 0; } -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// + + 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) @@ -353,8 +330,8 @@ bool FastMTT::compareLeptons( return false; return (measuredTauLepton1.pt() > measuredTauLepton2.pt()); } -/////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////// + + ROOT::Math::PtEtaPhiMVector FastMTT::run(const std::vector &measuredTauLeptons, const double &measuredMETx, const double &measuredMETy, @@ -364,15 +341,12 @@ FastMTT::run(const std::vector &measuredTauLeptons, myLikelihood.setParameters(shapeParams); FastMTT::initialize(); bestP4 = LorentzVector(); - //////////////////////////////////////////// - if (measuredTauLeptons.size() != 2) { std::cout << "Number of MeasuredTauLepton is " << measuredTauLeptons.size() << " a user shouls pass exactly two leptons." << std::endl; return ROOT::Math::PtEtaPhiMVector(); } - std::vector sortedMeasuredTauLeptons = measuredTauLeptons; std::sort(sortedMeasuredTauLeptons.begin(), sortedMeasuredTauLeptons.end(), @@ -404,9 +378,7 @@ FastMTT::run(const std::vector &measuredTauLeptons, void FastMTT::scan(Likelihood &myLikelihood) { 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; @@ -423,7 +395,6 @@ void FastMTT::scan(Likelihood &myLikelihood) { theMinimum[1] = x[1]; } } - minimumPosition[0] = theMinimum[0]; minimumPosition[1] = theMinimum[1]; minimumValue = bestLH;