Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix FastMTT memory leak #269

Merged
merged 2 commits into from
Jun 5, 2024
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 6 additions & 73 deletions include/SVFit/FastMTT.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include "Math/LorentzVector.h"
#include "TBenchmark.h"
#include "TMatrixD.h"
#include <Math/Vector4D.h>
#include <bitset>
#include <string>
#include <tuple>
Expand Down Expand Up @@ -72,7 +73,7 @@ class Likelihood {
std::vector<double> 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;
Expand All @@ -94,81 +95,19 @@ class 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);
ROOT::Math::PtEtaPhiMVector
run(const std::vector<fastmtt::MeasuredTauLepton> &, 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<double> minimumPosition;
Expand All @@ -189,14 +128,8 @@ class FastMTT {
/// Step sizes for each minimized variable
std::vector<double> stepSizes;

ROOT::Math::Functor *likelihoodFunctor;

Likelihood myLikelihood;

LorentzVector tau1P4, tau2P4, bestP4;

TBenchmark clock;

int verbosity;
};

Expand Down
159 changes: 28 additions & 131 deletions src/SVFit/FastMTT.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -7,25 +7,27 @@
#include "TMath.h"
#include "TMatrixD.h"
#include "TVector2.h"
#include <Math/Vector4D.h>

#include "Math/BasicMinimizer.h"
#include "TF1.h"

#include "../../include/SVFit/FastMTT.hxx"
#include "../../include/SVFit/MeasuredTauLepton.hxx"
#include "../../include/utility/Logger.hxx"

///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
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);
harrypuuter marked this conversation as resolved.
Show resolved Hide resolved
}
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -107,13 +109,13 @@ void Likelihood::setParameters(const std::vector<double> &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);
}
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -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<std::string> varNames = {"x1", "x2"};
nVariables = varNames.size();
std::vector<double> initialValues(nVariables, 0.5);
std::vector<double> 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<double> 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<double> &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) {
Expand All @@ -397,19 +355,22 @@ bool FastMTT::compareLeptons(
}
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
void FastMTT::run(
const std::vector<fastmtt::MeasuredTauLepton> &measuredTauLeptons,
const double &measuredMETx, const double &measuredMETy,
const TMatrixD &covMET) {

ROOT::Math::PtEtaPhiMVector
FastMTT::run(const std::vector<fastmtt::MeasuredTauLepton> &measuredTauLeptons,
const double &measuredMETx, const double &measuredMETy,
const TMatrixD &covMET) {
Likelihood myLikelihood;
std::vector<double> shapeParams = {6, 1.0 / 1.15};
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;
return ROOT::Math::PtEtaPhiMVector();
}

std::vector<fastmtt::MeasuredTauLepton> sortedMeasuredTauLeptons =
Expand All @@ -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;

Expand All @@ -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<double, double> 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);
}
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
}
Loading
Loading