MinuitFit Class Reference

Implementation of Minuit fitting algorithm [17]. More...

#include <Fitters/MinuitFit.h>

Inheritance diagram for MinuitFit:

Collaboration diagram for MinuitFit:

Public Member Functions
	MinuitFit (manager *const fitMan)
	Constructor. More...
virtual	~MinuitFit ()
	Destructor. More...
void	RunMCMC () override
	Actual implementation of Minuit Fit algorithm. More...
Public Member Functions inherited from LikelihoodFit
	LikelihoodFit (manager *const fitMan)
	Constructor. More...
virtual	~LikelihoodFit ()
	Destructor. More...
virtual double	CalcChi2 (const double *x)
	Chi2 calculation over all included samples and syst objects. More...
int	GetNPars ()
	Get total number of params, this sums over all covariance objects. More...
Public Member Functions inherited from FitterBase
	FitterBase (manager *const fitMan)
	Constructor. More...
virtual	~FitterBase ()
	Destructor for the FitterBase class. More...
void	AddSampleHandler (SampleHandlerBase *sample)
	This function adds a sample PDF object to the analysis framework. The sample PDF object will be utilized in fitting procedures or likelihood scans. More...
void	AddSystObj (ParameterHandlerBase *cov)
	This function adds a Covariance object to the analysis framework. The Covariance object will be utilized in fitting procedures or likelihood scans. More...
void	DragRace (const int NLaps=100)
	Calculates the required time for each sample or covariance object in a drag race simulation. Inspired by Dan's feature. More...
void	RunLLHScan ()
	Perform a 1D likelihood scan. More...
void	GetStepScaleBasedOnLLHScan ()
	LLH scan is good first estimate of step scale. More...
void	Run2DLLHScan ()
	Perform a 2D likelihood scan. More...
void	RunSigmaVar ()
	Perform a 2D and 1D sigma var for all samples. More...
void	RunSigmaVarFD ()
	Perform a 1D sigma var for all samples. More...
virtual void	StartFromPreviousFit (const std::string &FitName)
	Allow to start from previous fit/chain. More...
std::string	GetName () const
	Get name of class. More...

Private Attributes
std::unique_ptr< ROOT::Math::Minimizer >	minuit
	Pointer to minimizer, which most often is Minuit. More...

Additional Inherited Members
Protected Member Functions inherited from LikelihoodFit
void	PrepareFit ()
	prepare output and perform sanity checks More...
Protected Member Functions inherited from FitterBase
void	ProcessMCMC ()
	Process MCMC output. More...
void	PrepareOutput ()
	Prepare the output file. More...
void	SaveOutput ()
	Save output and close files. More...
void	SanitiseInputs ()
	Remove obsolete memory and make other checks before fit starts. More...
void	SaveSettings ()
	Save the settings that the MCMC was run with. More...
bool	GetScaneRange (std::map< std::string, std::vector< double >> &scanRanges)
	YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D. Barrow. More...
bool	CheckSkipParameter (const std::vector< std::string > &SkipVector, const std::string &ParamName) const
	KS: Check whether we want to skip parameter using skip vector. More...
void	CustomRange (const std::string &ParName, const double sigma, double &ParamShiftValue)
	For comparison with P-Theta we usually have to apply different parameter values then usual 1, 3 sigma. More...
Protected Attributes inherited from LikelihoodFit
int	NPars
	Number of all parameters from all covariances. More...
int	NParsPCA
	Number of all parameters from all covariances in PCA base. More...
bool	fMirroring
	Flag telling if mirroring is used or not. More...
Protected Attributes inherited from FitterBase
manager *	fitMan
	The manager. More...
unsigned int	step
	current state More...
double	logLCurr
	current likelihood More...
double	logLProp
	proposed likelihood More...
double	accProb
	current acceptance prob More...
int	accCount
	counts accepted steps More...
unsigned int	stepStart
	step start, by default 0 if we start from previous chain then it will be different More...
std::vector< double >	sample_llh
	store the llh breakdowns More...
std::vector< double >	syst_llh
	systematic llh breakdowns More...
std::vector< SampleHandlerBase * >	samples
	Sample holder. More...
unsigned int	TotalNSamples
	Total number of samples used. More...
std::vector< ParameterHandlerBase * >	systematics
	Systematic holder. More...
std::unique_ptr< TStopwatch >	clock
	tells global time how long fit took More...
std::unique_ptr< TStopwatch >	stepClock
	tells how long single step/fit iteration took More...
double	stepTime
	Time of single step. More...
std::unique_ptr< TRandom3 >	random
	Random number. More...
TFile *	outputFile
	Output. More...
TDirectory *	CovFolder
	Output cov folder. More...
TDirectory *	SampleFolder
	Output sample folder. More...
TTree *	outTree
	Output tree with posteriors. More...
int	auto_save
	auto save every N steps More...
bool	fTestLikelihood
	Necessary for some fitting algorithms like PSO. More...
bool	FileSaved
	Checks if file saved not repeat some operations. More...
bool	SettingsSaved
	Checks if setting saved not repeat some operations. More...
bool	OutputPrepared
	Checks if output prepared not repeat some operations. More...
std::string	AlgorithmName
	Name of fitting algorithm that is being used. More...

Detailed Description

Implementation of Minuit fitting algorithm [17].

Author

Kamil Skwarczynski

Definition at line 14 of file MinuitFit.h.

Constructor & Destructor Documentation

MinuitFit()

MinuitFit::MinuitFit

(

manager *const

fitMan

)

Constructor.

Todo:

KS: Make this in future configurable, for more see: https://root.cern.ch/doc/master/classROOT_1_1Math_1_1Minimizer.html

Definition at line 5 of file MinuitFit.cpp.

                                 : LikelihoodFit(man) {
// *******************
  AlgorithmName = "MinuitFit";
  // Minimizer type: determines the underlying implementation.
  // Available types include:
  //   - "Minuit2" (recommended modern option)
  //   - "Minuit"  (legacy)
  //   - "Fumili"
  //   - "GSLMultiMin" (for gradient-free minimization)
  //   - "GSLMultiFit"
  //   - "GSLSimAn" (Simulated Annealing)
  const std::string MinimizerType = "Minuit2";
  // Minimizer algorithm (specific to the selected type).
  // For Minuit2, the following algorithms are available:
  //   - "Migrad"   : gradient-based minimization (default)
  //   - "Simplex"  : Nelder-Mead simplex method (derivative-free)
  //   - "Combined" : combination of Simplex and Migrad
  //   - "Scan"     : parameter grid scan
  const std::string MinimizerAlgo = "Migrad";
 
  MACH3LOG_INFO("Creating instance of Minimizer with {} and {}", MinimizerType, MinimizerAlgo);
 
  minuit = std::unique_ptr<ROOT::Math::Minimizer>(
    ROOT::Math::Factory::CreateMinimizer(MinimizerType.c_str(), MinimizerAlgo.c_str()));
}

~MinuitFit()

MinuitFit::~MinuitFit ( )

virtual

Destructor.

Definition at line 34 of file MinuitFit.cpp.

                      {
// *************************
}

Member Function Documentation

RunMCMC()

void MinuitFit::RunMCMC ( )

overridevirtual

Actual implementation of Minuit Fit algorithm.

Implements FitterBase.

Definition at line 41 of file MinuitFit.cpp.

                        {
// *******************
  PrepareFit();
 
  // Remove obsolete memory and make other checks before fit starts
  SanitiseInputs();
 
  //KS: For none PCA this will be equal to normal parameters
  const int NparsMinuitFull = NPars;
  const int NparsMinuit = NParsPCA;
 
  //KS: Set SetFunction we will Minimize
  ROOT::Math::Functor fChi2(this, &MinuitFit::CalcChi2, NparsMinuit);
  minuit->SetFunction(fChi2);
 
  //KS: add config or something
  minuit->SetPrintLevel(2);
  minuit->SetTolerance(0.01);
  minuit->SetMaxFunctionCalls(fitMan->raw()["General"]["Minuit2"]["NSteps"].as<unsigned>());
  minuit->SetMaxIterations(10000);
 
  MACH3LOG_INFO("Preparing Minuit");
  int ParCounter = 0;
 
  for (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
  {
    if(!(*it)->IsPCA())
    {
      for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
      {
        //KS: Index, name, prior, step scale [different to MCMC],
        minuit->SetVariable(ParCounter, ((*it)->GetParName(i)), (*it)->GetParInit(i), (*it)->GetDiagonalError(i)/10);
        minuit->SetVariableValue(ParCounter, (*it)->GetParInit(i));
        //KS: lower bound, upper bound, if Mirroring enabled then ignore
        if(!fMirroring) minuit->SetVariableLimits(ParCounter, (*it)->GetLowerBound(i), (*it)->GetUpperBound(i));
        if((*it)->IsParameterFixed(i))
        {
          minuit->FixVariable(ParCounter);
        }
      }
    }
    else
    {
      for(int i = 0; i < (*it)->GetNParameters(); ++i, ++ParCounter)
      {
        minuit->SetVariable(ParCounter, Form("%i_PCA", i), (*it)->GetPCAHandler()->GetParPropPCA(i), (*it)->GetPCAHandler()->GetEigenValuesMaster()[i]/10);
        if((*it)->GetPCAHandler()->IsParameterFixedPCA(i))
        {
          minuit->FixVariable(ParCounter);
        }
      }
    }
  }
 
  minuit->SetPrintLevel(2);
 
  MACH3LOG_INFO("Starting MIGRAD");
  minuit->Minimize();
 
  MACH3LOG_INFO("Starting HESSE");
  minuit->Hesse();
  outputFile->cd();
  
  TVectorD* MinuitParValue = new TVectorD(NparsMinuitFull);
  TVectorD* MinuitParError = new TVectorD(NparsMinuitFull);
  TMatrixDSym* Postmatrix = new TMatrixDSym(NparsMinuitFull);
 
  for(int i = 0; i < NparsMinuitFull; ++i)
  {
    (*MinuitParValue)(i) = 0;
    (*MinuitParError)(i) = 0;
    for(int j = 0; j < NparsMinuitFull; ++j)
    {
      (*Postmatrix)(i,j) = 0;
      (*Postmatrix)(i,j) = minuit->CovMatrix(i,j);
    }
  }
 
  ParCounter = 0;
  const double *X = minuit->X();
  const double *err = minuit->Errors();
  for (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
  {
    if(!(*it)->IsPCA())
    {
      for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
      {
        double ParVal = X[ParCounter];
        //KS: Basically apply mirroring for parameters out of bounds
        if(fMirroring)
        {
          if(ParVal < (*it)->GetLowerBound(i))
          {
            ParVal = (*it)->GetLowerBound(i) + ((*it)->GetLowerBound(i) - ParVal);
          }
          else if (ParVal > (*it)->GetUpperBound(i))
          {
            ParVal = (*it)->GetUpperBound(i) - ( ParVal - (*it)->GetUpperBound(i));
          }
        }
        (*MinuitParValue)(ParCounter) = ParVal;
        (*MinuitParError)(ParCounter) = err[ParCounter];
        //KS: For fixed params HESS will not calculate error so we need to pass prior error
        if((*it)->IsParameterFixed(i))
        {
          (*MinuitParError)(ParCounter) = (*it)->GetDiagonalError(i);
          (*Postmatrix)(ParCounter,ParCounter) = (*MinuitParError)(ParCounter) * (*MinuitParError)(ParCounter);
        }
      }
    }
    else
    {
      //KS: We need to convert parameters from PCA to normal base
      TVectorD ParVals((*it)->GetNumParams());
      TVectorD ParVals_PCA((*it)->GetNParameters());
 
      TVectorD ErrorVals((*it)->GetNumParams());
      TVectorD ErrorVals_PCA((*it)->GetNParameters());
 
      TMatrixD MatrixVals((*it)->GetNumParams(), (*it)->GetNumParams());
      TMatrixD MatrixVals_PCA((*it)->GetNParameters(), (*it)->GetNParameters());
 
      //First save them
      //KS: This code is super convoluted as MaCh3 can store separate matrices while Minuit has one matrix. In future this will be simplified, keep it like this for now.
      const int StartVal = ParCounter;
      for(int i = 0; i < (*it)->GetNParameters(); ++i, ++ParCounter)
      {
        ParVals_PCA(i) = X[ParCounter];
        ErrorVals_PCA(i) = err[ParCounter];
        int ParCounterMatrix = StartVal;
        for(int j = 0; j < (*it)->GetNParameters(); ++j, ++ParCounterMatrix)
        {
          MatrixVals_PCA(i,j) = minuit->CovMatrix(ParCounter,ParCounterMatrix);
        }
      }
      ParVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ParVals_PCA;
      ErrorVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ErrorVals_PCA;
      MatrixVals.Mult(((*it)->GetPCAHandler()->GetTransferMatrix()),MatrixVals_PCA);
 
      ParCounter = StartVal;
      //KS: Now after going from PCA to normal let';s save it
      for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
      {
        (*MinuitParValue)(ParCounter) = ParVals(i);
        (*MinuitParError)(ParCounter) = std::fabs(ErrorVals(i));
        int ParCounterMatrix = StartVal;
        for(int j = 0; j < (*it)->GetNumParams(); ++j, ++ParCounterMatrix)
        {
          (*Postmatrix)(ParCounter,ParCounterMatrix)  = MatrixVals(i,j);
        }
        //If fixed take prior
        if((*it)->GetPCAHandler()->IsParameterFixedPCA(i))
        {
          (*MinuitParError)(ParCounter) = (*it)->GetDiagonalError(i);
          (*Postmatrix)(ParCounter,ParCounter) = (*MinuitParError)(ParCounter) * (*MinuitParError)(ParCounter);
        }
      }
    }
  }
 
  MinuitParValue->Write("MinuitParValue");
  MinuitParError->Write("MinuitParError");
  Postmatrix->Write("Postmatrix");
  delete MinuitParValue;
  delete MinuitParError;
  delete Postmatrix;
  // Save all the output
  SaveOutput();
}

Member Data Documentation

minuit

std::unique_ptr<ROOT::Math::Minimizer> MinuitFit::minuit

private

Pointer to minimizer, which most often is Minuit.

Definition at line 25 of file MinuitFit.h.

---

The documentation for this class was generated from the following files:

• Fitters/MinuitFit.h
• Fitters/MinuitFit.cpp