FitterBase Class Reference

Base class for implementing fitting algorithms. More...

#include <Fitters/FitterBase.h>

Inheritance diagram for FitterBase:

Collaboration diagram for FitterBase:

Public Member Functions
	FitterBase (Manager *const fitMan)
	Constructor. More...
virtual	~FitterBase ()
	Destructor for the FitterBase class. More...
void	AddSampleHandler (SampleHandlerInterface *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...
virtual void	RunMCMC ()=0
	The specific fitting algorithm implemented in this function depends on the derived class. It could be Markov Chain Monte Carlo (MCMC), MinuitFit, or another algorithm. 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	RunLLHMap ()
	Perform a general multi-dimensional likelihood scan. More...
void	GetStepScaleBasedOnLLHScan (const std::string &filename="")
	LLH scan is good first estimate of step scale. More...
void	Run2DLLHScan ()
	Perform a 2D likelihood scan. More...
void	RunSigmaVar ()
	Perform a 1D/2D 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...

Protected Member Functions
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	GetScanRange (std::map< std::string, std::vector< double >> &scanRanges) const
	YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D. Barrow. More...
void	GetParameterScanRange (const ParameterHandlerBase *cov, const int i, double &CentralValue, double &lower, double &upper, const int n_points, const std::string &suffix="") const
	Helper function to get parameter scan range, central value. 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) const
	For comparison with other fitting frameworks (like P-Theta) we usually have to apply different parameter values then usual 1, 3 sigma. More...

Protected Attributes
Manager *	fitMan
	The manager for configuration handling. 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< SampleHandlerInterface * >	samples
	Sample holder. More...
unsigned int	TotalNSamples
	Total number of samples used, single SampleHandler can store more than one analysis sample! 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

Base class for implementing fitting algorithms.

This class wraps MaCh3 classes like SampleHandlerInterface and ParameterHandlerBase. It serves as a base for different fitting algorithms and for validation techniques such as LLH scans.

Definition at line 26 of file FitterBase.h.

Constructor & Destructor Documentation

FitterBase()

_MaCh3_Safe_Include_Start_ _MaCh3_Safe_Include_End_ FitterBase::FitterBase

(

Manager *const

fitMan

)

Constructor.

Parameters

fitMan

A pointer to a manager object, which will handle all settings.

Definition at line 15 of file FitterBase.cpp.

                                          : fitMan(man) {
// *************************
  AlgorithmName = "";
  //Get mach3 modes from Manager
  random = std::make_unique<TRandom3>(Get<int>(fitMan->raw()["General"]["Seed"], __FILE__, __LINE__));
 
  // Counter of the accepted # of steps
  accCount = 0;
  step = 0;
  stepStart = 0;
 
  clock = std::make_unique<TStopwatch>();
  stepClock = std::make_unique<TStopwatch>();
  #ifdef MACH3_DEBUG
  // Fit summary and debug info
  debug = GetFromManager<bool>(fitMan->raw()["General"]["Debug"], false, __FILE__ , __LINE__);
  #endif
 
  auto outfile = Get<std::string>(fitMan->raw()["General"]["OutputFile"], __FILE__ , __LINE__);
  // Save output every auto_save steps
  //you don't want this too often https://root.cern/root/html606/TTree_8cxx_source.html#l01229
  auto_save = Get<int>(fitMan->raw()["General"]["MCMC"]["AutoSave"], __FILE__ , __LINE__);
 
  // Set the output file
  outputFile = M3::Open(outfile, "RECREATE", __FILE__, __LINE__);
  outputFile->cd();
  // Set output tree
  outTree = new TTree("posteriors", "Posterior_Distributions");
  // Auto-save every 200MB, the bigger the better https://root.cern/root/html606/TTree_8cxx_source.html#l01229
  outTree->SetAutoSave(-200E6);
 
  FileSaved = false;
  SettingsSaved = false;
  OutputPrepared = false;
 
  //Create TDirectory
  CovFolder = outputFile->mkdir("CovarianceFolder");
  outputFile->cd();
  SampleFolder = outputFile->mkdir("SampleFolder");
  outputFile->cd();
 
  #ifdef MACH3_DEBUG
  // Prepare the output log file
  if (debug) debugFile.open((outfile+".log").c_str());
  #endif
 
  TotalNSamples = 0;
  fTestLikelihood = GetFromManager<bool>(fitMan->raw()["General"]["Fitter"]["FitTestLikelihood"], false, __FILE__ , __LINE__);
}

~FitterBase()

FitterBase::~FitterBase ( )

virtual

Destructor for the FitterBase class.

Definition at line 67 of file FitterBase.cpp.

                        {
// *************************
  SaveOutput();
  if(outputFile != nullptr) delete outputFile;
  outputFile = nullptr;
  MACH3LOG_DEBUG("Closing MaCh3 Fitter Engine");
}

Member Function Documentation

AddSampleHandler()

void FitterBase::AddSampleHandler

(

SampleHandlerInterface *

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.

Parameters

sample

A pointer to a sample PDF object derived from ParameterHandlerBase.

Definition at line 260 of file FitterBase.cpp.

                                                                      {
// *************************
  // Check if any subsample name collides with already-registered subsamples
  for (const auto &s : samples) {
    for (int iExisting = 0; iExisting < s->GetNSamples(); ++iExisting) {
      for (int iNew = 0; iNew < sample->GetNSamples(); ++iNew) {
        if (s->GetSampleTitle(iExisting) == sample->GetSampleTitle(iNew)) {
          MACH3LOG_ERROR(
            "Duplicate sample title '{}' in handler {} detected: "
            "same title exist in handler ", sample->GetSampleTitle(iNew),
            sample->GetName(), s->GetName());
          throw MaCh3Exception(__FILE__, __LINE__);
        }
      }
    }
  }
 
  for (const auto &s : samples) {
    if (s->GetName() == sample->GetName()) {
      MACH3LOG_ERROR("SampleHandler with name '{}' already exists!", sample->GetName());
      MACH3LOG_ERROR("Is it intended?");
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
  }
  // Save additional info from samples
  SampleFolder->cd();
 
  sample->SaveAdditionalInfo(SampleFolder);
  TotalNSamples += sample->GetNSamples();
  MACH3LOG_INFO("Adding {} object, with {} samples", sample->GetName(), sample->GetNSamples());
  samples.push_back(sample);
  outputFile->cd();
}

AddSystObj()

void FitterBase::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.

Parameters

cov	A pointer to a Covariance object derived from ParameterHandlerBase.

Definition at line 296 of file FitterBase.cpp.

                                                            {
// *************************
  MACH3LOG_INFO("Adding systematic object {}, with {} params", cov->GetName(), cov->GetNumParams());
  // KS: Need to make sure we don't have params with same name, otherwise ROOT I/O and parts of MaCh3 will be terribly confused...
  for (size_t s = 0; s < systematics.size(); ++s)
  {
    for (int iPar = 0; iPar < systematics[s]->GetNumParams(); ++iPar)
    {
      for (int i = 0; i < cov->GetNumParams(); ++i)
      {
        if(systematics[s]->GetParName(iPar) == cov->GetParName(i)){
          MACH3LOG_ERROR("ParameterHandler {} has param '{}' which already exists in in {}, with name {}",
                         cov->GetName(), cov->GetParName(i), systematics[s]->GetName(), systematics[s]->GetParName(iPar));
          throw MaCh3Exception(__FILE__ , __LINE__ );
        }
        // Same for fancy name
        if(systematics[s]->GetParFancyName(iPar) == cov->GetParFancyName(i)){
          MACH3LOG_ERROR("ParameterHandler {} has param '{}' which already exists in {}, with name {}",
                         cov->GetName(), cov->GetParFancyName(i), systematics[s]->GetName(), systematics[s]->GetParFancyName(iPar));
          throw MaCh3Exception(__FILE__ , __LINE__ );
        }
      }
    }
  }
 
  systematics.push_back(cov);
 
  CovFolder->cd();
  std::vector<double> n_vec(cov->GetNumParams());
  for (int i = 0; i < cov->GetNumParams(); ++i) {
    n_vec[i] = cov->GetParInit(i);
  }
  cov->GetCovMatrix()->Write(cov->GetName().c_str());
 
  TH2D* CorrMatrix = cov->GetCorrelationMatrix();
  CorrMatrix->Write((cov->GetName() + std::string("_Corr")).c_str());
  delete CorrMatrix;
 
  // If we have yaml config file for covariance let's save it
  YAML::Node Config = cov->GetConfig();
  if(!Config.IsNull())
  {
    TMacro ConfigSave = YAMLtoTMacro(Config, (std::string("Config_") + cov->GetName()));
    ConfigSave.Write();
  }
 
  outputFile->cd();
}

CheckSkipParameter()

bool FitterBase::CheckSkipParameter ( const std::vector< std::string > &  SkipVector, const std::string &  ParamName  ) const

protected

KS: Check whether we want to skip parameter using skip vector.

Definition at line 540 of file FitterBase.cpp.

                                                                                                            {
// *************************
  bool skip = false;
  for(unsigned int is = 0; is < SkipVector.size(); ++is)
  {
    if(ParamName.substr(0, SkipVector[is].length()) == SkipVector[is])
    {
      skip = true;
      break;
    }
  }
  return skip;
}

CustomRange()

void FitterBase::CustomRange ( const std::string &  ParName, const double  sigma, double &  ParamShiftValue  ) const

protected

For comparison with other fitting frameworks (like P-Theta) we usually have to apply different parameter values then usual 1, 3 sigma.

Example YAML format:

SigmaVar:
  Q2_norm_7: { "3": 2.0 }
  SRC_Norm_O: { "-1": 0.5, "1": 1.5, "3": 2.0 }

Definition at line 1301 of file FitterBase.cpp.

                                                                                                        {
// *************************
  if(!fitMan->raw()["SigmaVar"]["CustomRange"]) return;
 
  auto Config = fitMan->raw()["SigmaVar"]["CustomRange"];
 
  const auto sigmaStr = std::to_string(static_cast<int>(std::round(sigma)));
 
  if (Config[ParName] && Config[ParName][sigmaStr]) {
    ParamShiftValue = Config[ParName][sigmaStr].as<double>();
    MACH3LOG_INFO("  ::: setting custom range from config ::: {} -> {}", ParName, ParamShiftValue);
  }
}

DragRace()

void FitterBase::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.

Parameters

NLaps

number of laps, every part of Fitter will be tested with given number of laps and you will get total and average time

Definition at line 458 of file FitterBase.cpp.

                                         {
// *************************
  MACH3LOG_INFO("Let the Race Begin!");
  MACH3LOG_INFO("All tests will be performed with {} threads", M3::GetNThreads());
 
  // Reweight the MC
  for(unsigned int ivs = 0; ivs < samples.size(); ++ivs)
  {
    TStopwatch clockRace;
    clockRace.Start();
    for(int Lap = 0; Lap < NLaps; ++Lap) {
      samples[ivs]->Reweight();
    }
    clockRace.Stop();
    MACH3LOG_INFO("It took {:.4f} s to reweights {} times sample: {}", clockRace.RealTime(), NLaps, samples[ivs]->GetName());
    MACH3LOG_INFO("On average {:.6f}", clockRace.RealTime()/NLaps);
  }
 
  for(unsigned int ivs = 0; ivs < samples.size(); ++ivs)
  {
    TStopwatch clockRace;
    clockRace.Start();
    for(int Lap = 0; Lap < NLaps; ++Lap) {
      samples[ivs]->GetLikelihood();
    }
    clockRace.Stop();
    MACH3LOG_INFO("It took {:.4f} s to calculate  GetLikelihood {} times sample:  {}", clockRace.RealTime(), NLaps, samples[ivs]->GetName());
    MACH3LOG_INFO("On average {:.6f}", clockRace.RealTime()/NLaps);
  }
  // Get vector of proposed steps. If we want to run LLH scan or something else after we need to revert changes after proposing steps multiple times
  std::vector<std::vector<double>> StepsValuesBefore(systematics.size());
  for (size_t s = 0; s < systematics.size(); ++s) {
    StepsValuesBefore[s] = systematics[s]->GetProposed();
  }
  for (size_t s = 0; s < systematics.size(); ++s) {
    TStopwatch clockRace;
    clockRace.Start();
    for(int Lap = 0; Lap < NLaps; ++Lap) {
      systematics[s]->ProposeStep();
    }
    clockRace.Stop();
    MACH3LOG_INFO("It took {:.4f} s to propose step {} times cov:  {}",  clockRace.RealTime(), NLaps, systematics[s]->GetName());
    MACH3LOG_INFO("On average {:.6f}", clockRace.RealTime()/NLaps);
  }
  for (size_t s = 0; s < systematics.size(); ++s) {
    systematics[s]->SetParameters(StepsValuesBefore[s]);
  }
 
  for (size_t s = 0; s < systematics.size(); ++s) {
    TStopwatch clockRace;
    clockRace.Start();
    for(int Lap = 0; Lap < NLaps; ++Lap) {
      systematics[s]->GetLikelihood();
    }
    clockRace.Stop();
    MACH3LOG_INFO("It took {:.4f} s to calculate  get likelihood {} times cov:  {}",  clockRace.RealTime(), NLaps, systematics[s]->GetName());
    MACH3LOG_INFO("On average {:.6f}", clockRace.RealTime()/NLaps);
  }
  MACH3LOG_INFO("End of race");
}

GetName()

std::string FitterBase::GetName ( ) const

inline

Get name of class.

Definition at line 72 of file FitterBase.h.

{return AlgorithmName;};

GetParameterScanRange()

void FitterBase::GetParameterScanRange ( const ParameterHandlerBase *  cov, const int  i, double &  CentralValue, double &  lower, double &  upper, const int  n_points, const std::string &  suffix = ""  ) const

protected

Helper function to get parameter scan range, central value.

Definition at line 556 of file FitterBase.cpp.

                                                                                                                        {
// *************************
  // YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D. Barrow
  std::map<std::string, std::vector<double>> scanRanges;
  const bool isScanRanges = GetScanRange(scanRanges);
 
  double nSigma = GetFromManager<int>(fitMan->raw()["LLHScan"]["LLHScanSigma"], 1., __FILE__, __LINE__);
  bool IsPCA = cov->IsPCA();
 
  // Get the parameter name
  std::string name = cov->GetParFancyName(i);
  if (IsPCA) name += "_PCA";
 
  // Get the parameter priors and bounds
  CentralValue = cov->GetParProp(i);
  if (IsPCA) CentralValue = cov->GetPCAHandler()->GetParPropPCA(i);
 
  double prior = cov->GetParInit(i);
  if (IsPCA) prior = cov->GetPCAHandler()->GetPreFitValuePCA(i);
 
  if (std::abs(CentralValue - prior) > 1e-10) {
    MACH3LOG_INFO("For {} scanning around value {} rather than prior {}", name, CentralValue, prior);
  }
 
  // Get the covariance matrix and do the +/- nSigma
  // Set lower and upper bounds relative the CentralValue
  // Set the parameter ranges between which LLH points are scanned
  lower = CentralValue - nSigma*cov->GetDiagonalError(i);
  upper = CentralValue + nSigma*cov->GetDiagonalError(i);
  // If PCA, transform these parameter values to the PCA basis
  if (IsPCA) {
    lower = CentralValue - nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
    upper = CentralValue + nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
    MACH3LOG_INFO("eval {} = {:.2f}", i, cov->GetPCAHandler()->GetEigenValues()(i));
    MACH3LOG_INFO("CV {} = {:.2f}", i, CentralValue);
    MACH3LOG_INFO("lower {} = {:.2f}", i, lower);
    MACH3LOG_INFO("upper {} = {:.2f}", i, upper);
    MACH3LOG_INFO("nSigma = {:.2f}", nSigma);
  }
 
  // Implementation suggested by D. Barrow
  // If param ranges are specified in scanRanges node, extract it from there
  if(isScanRanges){
    // Find matching entries through std::maps
    auto it = scanRanges.find(name);
    if (it != scanRanges.end() && it->second.size() == 2) { //Making sure the range is has only two entries
      lower = it->second[0];
      upper = it->second[1];
      MACH3LOG_INFO("Found matching param name for setting specified range for {}", name);
      MACH3LOG_INFO("Range for {} = [{:.2f}, {:.2f}]", name, lower, upper);
    }
  }
 
  // Cross-section and flux parameters have boundaries that we scan between, check that these are respected in setting lower and upper variables
  // This also applies for other parameters like osc, etc.
  lower = std::max(lower, cov->GetLowerBound(i));
  upper = std::min(upper, cov->GetUpperBound(i));
  MACH3LOG_INFO("Scanning {} {} with {} steps, from [{:.2f} , {:.2f}], CV = {:.2f}", suffix, name, n_points, lower, upper, CentralValue);
}

GetScanRange()

bool FitterBase::GetScanRange ( std::map< std::string, std::vector< double >> &  scanRanges ) const

protected

YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D. Barrow.

Parameters

scanRanges

A map with user specified parameter ranges

Definition at line 520 of file FitterBase.cpp.

                                                                                    {
// *************************
  bool isScanRanges = false;
  // YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D. Barrow
  if(fitMan->raw()["LLHScan"]["ScanRanges"]){
    YAML::Node scanRangesList = fitMan->raw()["LLHScan"]["ScanRanges"];
    for (auto it = scanRangesList.begin(); it != scanRangesList.end(); ++it) {
      std::string itname = it->first.as<std::string>();
      std::vector<double> itrange = it->second.as<std::vector<double>>();
      // Set the mapping as param_name:param_range
      scanRanges[itname] = itrange;
    }
    isScanRanges = true;
  } else {
    MACH3LOG_INFO("There are no user-defined parameter ranges, so I'll use default param bounds for LLH Scans");
  }
  return isScanRanges;
}

GetStepScaleBasedOnLLHScan()

void FitterBase::GetStepScaleBasedOnLLHScan

(

const std::string &

filename = ""

)

LLH scan is good first estimate of step scale.

Parameters

filename

by default empty, however if specified it will allow to load LLH scan from external file

Definition at line 884 of file FitterBase.cpp.

                                                                           {
// *************************
  TFile* outputFileLLH = nullptr;
  bool ownsfile = false;
  if(outputFileName != ""){
    outputFileLLH = M3::Open(outputFileName, "READ", __FILE__, __LINE__);
    ownsfile = true;
  } else {
    outputFileLLH = outputFile;
  }
  TDirectory *Sample_LLH = outputFileLLH->Get<TDirectory>("Sample_LLH");
  MACH3LOG_INFO("Starting Get Step Scale Based On LLHScan");
 
  if(!Sample_LLH || Sample_LLH->IsZombie())
  {
    MACH3LOG_WARN("Couldn't find Sample_LLH, it looks like LLH scan wasn't run, will do this now");
    RunLLHScan();
    Sample_LLH = outputFileLLH->Get<TDirectory>("Sample_LLH");
  }
  
  for (ParameterHandlerBase *cov : systematics)
  {
    const int npars = cov->GetNumParams();
    std::vector<double> StepScale(npars);
    for (int i = 0; i < npars; ++i)
    {
      std::string name = cov->GetParFancyName(i);
      StepScale[i] = cov->GetIndivStepScale(i);
      TH1D* LLHScan = Sample_LLH->Get<TH1D>((name+"_sam").c_str());
      if(LLHScan == nullptr)
      {
        MACH3LOG_WARN("Couldn't find LLH scan, for {}, skipping", name);
        continue;
      }
      const double LLH_val = std::max(LLHScan->GetBinContent(1), LLHScan->GetBinContent(LLHScan->GetNbinsX()));
      //If there is no sensitivity leave it
      if(LLH_val < 0.001) continue;
 
      // EM: assuming that the likelihood is gaussian, approximate sigma value is given by variation/sqrt(-2LLH)
      // can evaluate this at any point, simple to evaluate it in the first bin of the LLH scan
      // KS: We assume variation is 1 sigma, each dial has different scale so it becomes faff...
      const double Var = 1.;
      const double approxSigma = std::abs(Var)/std::sqrt(LLH_val);
      const double GlobalScale = cov->GetGlobalStepScale();
      // Based on Ewan comment I just took the 1sigma width from the LLH, assuming it was Gaussian, but then had to also scale by 2.38/sqrt(N_params)
      const double TargetStep = approxSigma * 2.38 / std::sqrt(npars);
      // KS: Need to divide by currently used gloalStepScale
      const double NewStepScale = TargetStep / GlobalScale;
 
      StepScale[i] = NewStepScale;
      MACH3LOG_DEBUG("Sigma: {}", approxSigma);
      MACH3LOG_DEBUG("Target Step Size (before accounting for global step size): {}", TargetStep);
      MACH3LOG_DEBUG("Optimal Step Size: {}", NewStepScale);
    }
    cov->SetIndivStepScale(StepScale);
    cov->SaveUpdatedMatrixConfig();
  }
  if(ownsfile && outputFileLLH != nullptr) delete outputFileLLH;
}

PrepareOutput()

void FitterBase::PrepareOutput ( )

protected

Prepare the output file.

Definition at line 151 of file FitterBase.cpp.

                               {
// *******************
  if(OutputPrepared) return;
  //MS: Check if we are fitting the test likelihood, rather than T2K likelihood, and only setup T2K output if not
  if(!fTestLikelihood)
  {
    // Check that we have added samples
    if (!samples.size()) {
      MACH3LOG_CRITICAL("No samples Found! Is this really what you wanted?");
      //throw MaCh3Exception(__FILE__ , __LINE__ );
    }
 
    // Do we want to save proposal? This will break plotting scripts and is heave for disk space and step time. Only use when debugging
    bool SaveProposal = GetFromManager<bool>(fitMan->raw()["General"]["SaveProposal"], false, __FILE__ , __LINE__);
 
    if(SaveProposal) MACH3LOG_INFO("Will save in the chain proposal parameters and LogL");
    // Prepare the output trees
    for (ParameterHandlerBase *cov : systematics) {
      cov->SetBranches(*outTree, SaveProposal);
    }
 
    outTree->Branch("LogL", &logLCurr, "LogL/D");
    if(SaveProposal) outTree->Branch("LogLProp", &logLProp, "LogLProp/D");
    outTree->Branch("accProb", &accProb, "accProb/D");
    outTree->Branch("step", &step, "step/i");
    outTree->Branch("stepTime", &stepTime, "stepTime/D");
 
    // Store individual likelihood components
    // Using double means double as large output file!
    sample_llh.resize(samples.size());
    syst_llh.resize(systematics.size());
 
    for (size_t i = 0; i < samples.size(); ++i) {
      std::stringstream oss, oss2;
      oss << "LogL_sample_" << i;
      oss2 << oss.str() << "/D";
      outTree->Branch(oss.str().c_str(), &sample_llh[i], oss2.str().c_str());
    }
 
    for (size_t i = 0; i < systematics.size(); ++i) {
      std::stringstream oss, oss2;
      oss << "LogL_systematic_" << systematics[i]->GetName();
      oss2 << oss.str() << "/D";
      outTree->Branch(oss.str().c_str(), &syst_llh[i], oss2.str().c_str());
    }
  }
  else
  {
    outTree->Branch("LogL", &logLCurr, "LogL/D");
    outTree->Branch("accProb", &accProb, "accProb/D");
    outTree->Branch("step", &step, "step/i");
    outTree->Branch("stepTime", &stepTime, "stepTime/D");
  }
 
  MACH3LOG_INFO("-------------------- Starting MCMC --------------------");
  #ifdef MACH3_DEBUG
  if (debug) {
    debugFile << "----- Starting MCMC -----" << std::endl;
  }
  #endif
  // Time the progress
 
 
  clock->Start();
 
 
  OutputPrepared = true;
}

ProcessMCMC()

void FitterBase::ProcessMCMC ( )

protected

Process MCMC output.

Definition at line 411 of file FitterBase.cpp.

                             {
// *******************
  if (fitMan == nullptr) return;
 
  // Process the MCMC
  if (GetFromManager<bool>(fitMan->raw()["General"]["ProcessMCMC"], false, __FILE__ , __LINE__)){
    // Make the processor
    MCMCProcessor Processor(std::string(outputFile->GetName()));
 
    Processor.Initialise();
    // Make the TVectorD pointers which hold the processed output
    TVectorD *Central = nullptr;
    TVectorD *Errors = nullptr;
    TVectorD *Central_Gauss = nullptr;
    TVectorD *Errors_Gauss = nullptr;
    TVectorD *Peaks = nullptr;
 
    // Make the postfit
    Processor.GetPostfit(Central, Errors, Central_Gauss, Errors_Gauss, Peaks);
    Processor.DrawPostfit();
 
    // Make the TMatrix pointers which hold the processed output
    TMatrixDSym *Covariance = nullptr;
    TMatrixDSym *Correlation = nullptr;
 
    // Make the covariance matrix
    Processor.GetCovariance(Covariance, Correlation);
    Processor.DrawCovariance();
 
    std::vector<TString> BranchNames = Processor.GetBranchNames();
 
    // Re-open the TFile
    if (!outputFile->IsOpen()) {
      MACH3LOG_INFO("Opening output again to update with means..");
      outputFile = new TFile(Get<std::string>(fitMan->raw()["General"]["OutputFile"], __FILE__, __LINE__).c_str(), "UPDATE");
    }
    Central->Write("PDF_Means");
    Errors->Write("PDF_Errors");
    Central_Gauss->Write("Gauss_Means");
    Errors_Gauss->Write("Errors_Gauss");
    Covariance->Write("Covariance");
    Correlation->Write("Correlation");
  }
}

Run2DLLHScan()

void FitterBase::Run2DLLHScan

(

)

Perform a 2D likelihood scan.

Warning

This operation may take a significant amount of time, especially for complex models.

Definition at line 946 of file FitterBase.cpp.

                              {
// *************************
  // Save the settings into the output file
  SaveSettings();
 
  MACH3LOG_INFO("Starting 2D LLH Scan");
 
  TDirectory *Sample_2DLLH = outputFile->mkdir("Sample_2DLLH");
  auto SkipVector = GetFromManager<std::vector<std::string>>(fitMan->raw()["LLHScan"]["LLHScanSkipVector"], {}, __FILE__ , __LINE__);;
 
  // Number of points we do for each LLH scan
  const int n_points = GetFromManager<int>(fitMan->raw()["LLHScan"]["2DLLHScanPoints"], 20, __FILE__ , __LINE__);
  // We print 5 reweights
  const int countwidth = int(double(n_points)/double(5));
 
  // Loop over the covariance classes
  for (ParameterHandlerBase *cov : systematics)
  {
    // Scan over all the parameters
    // Get the number of parameters
    int npars = cov->GetNumParams();
    bool IsPCA = cov->IsPCA();
    if (IsPCA) npars = cov->GetNParameters();
 
    for (int i = 0; i < npars; ++i)
    {
      std::string name_x = cov->GetParFancyName(i);
      if (IsPCA) name_x += "_PCA";
      // Get the parameter central and bounds
      double central_x, lower_x, upper_x;
      GetParameterScanRange(cov, i, central_x, lower_x, upper_x, n_points, "X");
 
      // KS: Check if we want to skip this parameter
      if(CheckSkipParameter(SkipVector, name_x)) continue;
 
      for (int j = 0; j < i; ++j)
      {
        std::string name_y = cov->GetParFancyName(j);
        if (IsPCA) name_y += "_PCA";
        // KS: Check if we want to skip this parameter
        if(CheckSkipParameter(SkipVector, name_y)) continue;
 
        // Get the parameter central and bounds
        double central_y, lower_y, upper_y;
        GetParameterScanRange(cov, j, central_y, lower_y, upper_y, n_points, "Y");
 
        auto hScanSam = std::make_unique<TH2D>((name_x + "_" + name_y + "_sam").c_str(), (name_x + "_" + name_y + "_sam").c_str(),
                                                n_points, lower_x, upper_x, n_points, lower_y, upper_y);
        hScanSam->SetDirectory(nullptr);
        hScanSam->GetXaxis()->SetTitle(name_x.c_str());
        hScanSam->GetYaxis()->SetTitle(name_y.c_str());
        hScanSam->GetZaxis()->SetTitle("2LLH_sam");
 
        // Scan over the parameter space
        for (int x = 0; x < n_points; ++x)
        {
          if (x % countwidth == 0)
            M3::Utils::PrintProgressBar(x, n_points);
 
          for (int y = 0; y < n_points; ++y)
          {
            // For PCA we have to do it differently
            if (IsPCA) {
              cov->GetPCAHandler()->SetParPropPCA(i, hScanSam->GetXaxis()->GetBinCenter(x+1));
              cov->GetPCAHandler()->SetParPropPCA(j, hScanSam->GetYaxis()->GetBinCenter(y+1));
            } else {
              // Set the parameter
              cov->SetParProp(i, hScanSam->GetXaxis()->GetBinCenter(x+1));
              cov->SetParProp(j, hScanSam->GetYaxis()->GetBinCenter(y+1));
            }
            // Reweight the MC
            for(unsigned int ivs = 0; ivs < samples.size(); ++ivs) {
              samples[ivs]->Reweight();
            }
 
            // Get the -log L likelihoods
            double samplellh = 0;
            for(unsigned int ivs = 0; ivs < samples.size(); ++ivs) {
              samplellh += samples[ivs]->GetLikelihood();
            }
            hScanSam->SetBinContent(x+1, y+1, 2*samplellh);
          }// end loop over y points
        } // end loop over x points
 
        Sample_2DLLH->cd();
        hScanSam->Write();
        // Reset the parameters to their central central values
        if (IsPCA) {
          cov->GetPCAHandler()->SetParPropPCA(i, central_x);
          cov->GetPCAHandler()->SetParPropPCA(j, central_y);
        } else {
          cov->SetParProp(i, central_x);
          cov->SetParProp(j, central_y);
        }
      } //end loop over systematics y
    }//end loop over systematics X
  }//end loop covariance classes
  Sample_2DLLH->Write();
  delete Sample_2DLLH;
}

RunLLHMap()

void FitterBase::RunLLHMap

(

)

Perform a general multi-dimensional likelihood scan.

Author

Tomas Nosek

Definition at line 1049 of file FitterBase.cpp.

                           {
// *************************
  // Save the settings into the output file
  SaveSettings();
 
  MACH3LOG_INFO("Starting {}", __func__);
 
  //KS: Turn it on if you want LLH scan for each ND sample separately, which increase time significantly but can be useful for validating new samples or dials.
  bool PlotLLHScanBySample = GetFromManager<bool>(fitMan->raw()["LLHScan"]["LLHScanBySample"], false, __FILE__ , __LINE__);
  auto ParamsOfInterest = GetFromManager<std::vector<std::string>>(fitMan->raw()["LLHScan"]["LLHParameters"], {}, __FILE__, __LINE__);
 
  if(ParamsOfInterest.empty()) {
    MACH3LOG_WARN("There were no LLH parameters of interest specified to run the LLHMap! LLHMap will not run at all ...");
    return;
  }
 
  // Parameters IDs within the covariance objects
  // ParamsCovIDs = {Name, CovObj, IDinCovObj}
  std::vector<std::tuple<std::string, ParameterHandlerBase*, int>> ParamsCovIDs;
  for(auto& p : ParamsOfInterest) {
    bool found = false;
    for(auto cov : systematics) {
      for(int c = 0; c < cov->GetNumParams(); ++c) {
        if(cov->GetParName(c) == p || cov->GetParFancyName(c) == p) {
          bool add = true;
          for(auto& pc : ParamsCovIDs) {
            if(std::get<1>(pc) == cov && std::get<2>(pc) == c)
            {
              MACH3LOG_WARN("Parameter {} as {}({}) listed multiple times for LLHMap, omitting and using only once!", p, cov->GetName(), c);
              add = false;
              break;
            }
          }
 
          if(add)
            ParamsCovIDs.push_back(std::make_tuple(p, cov, c));
 
          found = true;
          break;
        }
      }
      if(found)
        break;
    }
    if(found)
      MACH3LOG_INFO("Parameter {} found in {} at an index {}.", p, std::get<1>(ParamsCovIDs.back())->GetName(), std::get<2>(ParamsCovIDs.back()));
    else
      MACH3LOG_WARN("Parameter {} not found in any of the systematic covariance objects. Will not scan over this one!", p);
  }
 
  // ParamsRanges["parameter"] = {nPoints, {low, high}}
  std::map<std::string, std::pair<int, std::pair<double, double>>> ParamsRanges;
 
  MACH3LOG_INFO("======================================================================================");
  MACH3LOG_INFO("Performing a general multi-dimensional LogL map scan over following parameters ranges:");
  MACH3LOG_INFO("======================================================================================");
  unsigned long TotalPoints = 1;
 
  double nSigma = GetFromManager<int>(fitMan->raw()["LLHScan"]["LLHScanSigma"], 1., __FILE__, __LINE__);
 
  // TN: Setting up the scan ranges might look like a re-implementation of the
  // FitterBase::GetScanRange, but I guess the use-case here is a bit different.
  // Anyway, just in case, we can discuss and rewrite to everyone's liking!
  for(auto& p : ParamsCovIDs) {
    // Auxiliary vars to help readability
    std::string name = std::get<0>(p);
    int i = std::get<2>(p);
    ParameterHandlerBase* cov = std::get<1>(p);
 
    ParamsRanges[name].first = GetFromManager<int>(fitMan->raw()["LLHScan"]["LLHScanPoints"], 20, __FILE__, __LINE__);
    if(CheckNodeExists(fitMan->raw(),"LLHScan","ScanPoints"))
      ParamsRanges[name].first = GetFromManager<int>(fitMan->raw()["LLHScan"]["ScanPoints"][name], ParamsRanges[name].first, __FILE__, __LINE__);
 
    // Get the parameter priors and bounds
    double CentralValue = cov->GetParProp(i);
 
    bool IsPCA = cov->IsPCA();
    if (IsPCA)
      CentralValue = cov->GetPCAHandler()->GetParPropPCA(i);
 
    double prior = cov->GetParInit(i);
    if (IsPCA)
      prior = cov->GetPCAHandler()->GetPreFitValuePCA(i);
 
    if (std::abs(CentralValue - prior) > 1e-10) {
      MACH3LOG_INFO("For {} scanning around value {} rather than prior {}", name, CentralValue, prior);
    }
    // Get the covariance matrix and do the +/- nSigma
    // Set lower and upper bounds relative the CentralValue
    // Set the parameter ranges between which LLH points are scanned
    double lower = CentralValue - nSigma*cov->GetDiagonalError(i);
    double upper = CentralValue + nSigma*cov->GetDiagonalError(i);
    // If PCA, transform these parameter values to the PCA basis
    if (IsPCA) {
      lower = CentralValue - nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
      upper = CentralValue + nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
      MACH3LOG_INFO("eval {} = {:.2f}", i, cov->GetPCAHandler()->GetEigenValues()(i));
      MACH3LOG_INFO("CV {} = {:.2f}", i, CentralValue);
      MACH3LOG_INFO("lower {} = {:.2f}", i, lower);
      MACH3LOG_INFO("upper {} = {:.2f}", i, upper);
      MACH3LOG_INFO("nSigma = {:.2f}", nSigma);
    }
 
    ParamsRanges[name].second = {lower,upper};
 
    if(CheckNodeExists(fitMan->raw(),"LLHScan","ScanRanges"))
      ParamsRanges[name].second = GetFromManager<std::pair<double,double>>(fitMan->raw()["LLHScan"]["ScanRanges"][name], ParamsRanges[name].second, __FILE__, __LINE__);
 
    MACH3LOG_INFO("{} from {:.4f} to {:.4f} with a {:.5f} step ({} points total)",
                  name, ParamsRanges[name].second.first, ParamsRanges[name].second.second,
                  (ParamsRanges[name].second.second - ParamsRanges[name].second.first)/(ParamsRanges[name].first - 1.),
                  ParamsRanges[name].first);
 
    TotalPoints *= ParamsRanges[name].first;
  }
 
  // TN: Waiting for C++ 20 std::format() function
  MACH3LOG_INFO("In total, looping over {} points, from {} parameters. Estimates for run time:", TotalPoints, ParamsCovIDs.size());
  MACH3LOG_INFO("   1 s per point = {} hours", double(TotalPoints)/3600.);
  MACH3LOG_INFO(" 0.1 s per point = {} hours", double(TotalPoints)/36000.);
  MACH3LOG_INFO("0.01 s per point = {} hours", double(TotalPoints)/360000.);
  MACH3LOG_INFO("==================================================================================");
 
  const int countwidth = int(double(TotalPoints)/double(20));
 
  // Tree to store LogL values
  auto LLHMap = new TTree("llhmap", "LLH Map");
 
  std::vector<double> CovLogL(systematics.size());
  for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc )
  {
    std::string NameTemp = systematics[ivc]->GetName();
    NameTemp = NameTemp.substr(0, NameTemp.find("_cov")) + "_LLH";
    LLHMap->Branch(NameTemp.c_str(), &CovLogL[ivc]);
  }
 
  std::vector<double> SampleClassLogL(samples.size());
  for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
  {
    std::string NameTemp = samples[ivs]->GetName()+"_LLH";
    LLHMap->Branch(NameTemp.c_str(), &SampleClassLogL[ivs]);
  }
 
  double SampleLogL, TotalLogL;
  LLHMap->Branch("Sample_LLH", &SampleLogL);
  LLHMap->Branch("Total_LLH", &TotalLogL);
 
  std::vector<double>SampleSplitLogL;
  if(PlotLLHScanBySample)
  {
    SampleSplitLogL.resize(TotalNSamples);
    int SampleIterator = 0;
    for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
    {
      for(int is = 0; is < samples[ivs]->GetNSamples(); ++is )
      {
        std::string NameTemp =samples[ivs]->GetSampleTitle(is)+"_LLH";
        LLHMap->Branch(NameTemp.c_str(), &SampleSplitLogL[SampleIterator]);
        SampleIterator++;
      }
    }
  }
 
  std::vector<double> ParamsValues(ParamsCovIDs.size());
  for(unsigned int i=0; i < ParamsCovIDs.size(); ++i)
    LLHMap->Branch(std::get<0>(ParamsCovIDs[i]).c_str(), &ParamsValues[i]);
 
  // Setting up the scan
  // Starting at index {0,0,0,...}
  std::vector<unsigned long> idx(ParamsCovIDs.size(), 0);
 
  // loop over scanned points sp
  for(unsigned long sp = 0; sp < TotalPoints; ++sp)
  {
    // At each point need to find the indices and test values to calculate LogL
    for(unsigned int n = 0; n < ParamsCovIDs.size(); ++n)
    {
      // Auxiliaries
      std::string name = std::get<0>(ParamsCovIDs[n]);
      int points  = ParamsRanges[name].first;
      double low  = ParamsRanges[name].second.first;
      double high = ParamsRanges[name].second.second;
 
      // Find the n-th index of the sp-th scanned point
      unsigned long dev = 1;
      for(unsigned int m = 0; m <= n; ++m)
        dev *= ParamsRanges[std::get<0>(ParamsCovIDs[m])].first;
 
      idx[n] = sp % dev;
      if (n > 0)
        idx[n] = idx[n] / ( dev / points );
 
      // Parameter test value = low + ( high - low ) * idx / ( #points - 1 )
      ParamsValues[n] = low + (high-low) * double(idx[n])/double(points-1);
 
      // Now set the covariance objects
      // Auxiliary
      ParameterHandlerBase* cov = std::get<1>(ParamsCovIDs[n]);
      int i = std::get<2>(ParamsCovIDs[n]);
 
      if(cov->IsPCA())
        cov->GetPCAHandler()->SetParPropPCA(i, ParamsValues[n]);
      else
        cov->SetParProp(i, ParamsValues[n]);
    }
 
    // Reweight samples and calculate LogL
    TotalLogL = .0;
    SampleLogL = .0;
 
    for(unsigned int ivs = 0; ivs < samples.size(); ++ivs)
      samples[ivs]->Reweight();
 
    for(unsigned int ivs = 0; ivs < samples.size(); ++ivs)
    {
      SampleClassLogL[ivs] = 2.*samples[ivs]->GetLikelihood();
      SampleLogL += SampleClassLogL[ivs];
    }
    TotalLogL += SampleLogL;
 
    // CovObjs LogL
    for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc )
    {
      CovLogL[ivc] = 2.*systematics[ivc]->GetLikelihood();
      TotalLogL += CovLogL[ivc];
    }
 
    if(PlotLLHScanBySample)
    {
      int SampleIterator = 0;
      for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
      {
        for(int is = 0; is < samples[ivs]->GetNSamples(); ++is)
        {
          SampleSplitLogL[SampleIterator] = 2.*samples[ivs]->GetSampleLikelihood(is);
          SampleIterator++;
        }
      }
    }
 
    LLHMap->Fill();
 
    if (sp % countwidth == 0)
      M3::Utils::PrintProgressBar(sp, TotalPoints);
  }
 
  outputFile->cd();
  LLHMap->Write();
}

RunLLHScan()

void FitterBase::RunLLHScan

(

)

Perform a 1D likelihood scan.

Definition at line 619 of file FitterBase.cpp.

                            {
// *************************
  // Save the settings into the output file
  SaveSettings();
 
  MACH3LOG_INFO("Starting {}", __func__);
 
  //KS: Turn it on if you want LLH scan for each ND sample separately, which increase time significantly but can be useful for validating new samples or dials.
  bool PlotLLHScanBySample = GetFromManager<bool>(fitMan->raw()["LLHScan"]["LLHScanBySample"], false, __FILE__ , __LINE__);
  auto SkipVector = GetFromManager<std::vector<std::string>>(fitMan->raw()["LLHScan"]["LLHScanSkipVector"], {}, __FILE__ , __LINE__);
 
  // Now finally get onto the LLH scan stuff
  // Very similar code to MCMC but never start MCMC; just scan over the parameter space
  std::vector<TDirectory *> Cov_LLH(systematics.size());
  for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc )
  {
    std::string NameTemp = systematics[ivc]->GetName();
    NameTemp = NameTemp.substr(0, NameTemp.find("_cov")) + "_LLH";
    Cov_LLH[ivc] = outputFile->mkdir(NameTemp.c_str());
  }
 
  std::vector<TDirectory *> SampleClass_LLH(samples.size());
  for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
  {
    std::string NameTemp = samples[ivs]->GetName();
    SampleClass_LLH[ivs] = outputFile->mkdir(NameTemp.c_str());
  }
 
  TDirectory *Sample_LLH = outputFile->mkdir("Sample_LLH");
  TDirectory *Total_LLH = outputFile->mkdir("Total_LLH");
 
  std::vector<TDirectory *>SampleSplit_LLH;
  if(PlotLLHScanBySample)
  {
    SampleSplit_LLH.resize(TotalNSamples);
    int SampleIterator = 0;
    for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
    {
      for(int is = 0; is < samples[ivs]->GetNSamples(); ++is )
      {
        SampleSplit_LLH[SampleIterator] = outputFile->mkdir((samples[ivs]->GetSampleTitle(is)+ "_LLH").c_str());
        SampleIterator++;
      }
    }
  }
  // Number of points we do for each LLH scan
  const int n_points = GetFromManager<int>(fitMan->raw()["LLHScan"]["LLHScanPoints"], 100, __FILE__ , __LINE__);
 
  // We print 5 reweights
  const int countwidth = int(double(n_points)/double(5));
   
  // Loop over the covariance classes
  for (ParameterHandlerBase *cov : systematics)
  {
    // Scan over all the parameters
    // Get the number of parameters
    int npars = cov->GetNumParams();
    bool IsPCA = cov->IsPCA();
    if (IsPCA) npars = cov->GetNParameters();
    for (int i = 0; i < npars; ++i)
    {
      // Get the parameter name
      std::string name = cov->GetParFancyName(i);
      if (IsPCA) name += "_PCA";
      // KS: Check if we want to skip this parameter
      if(CheckSkipParameter(SkipVector, name)) continue;
      // Get the parameter central and bounds
      double CentralValue, lower, upper;
      GetParameterScanRange(cov, i, CentralValue, lower, upper, n_points);
      // Make the TH1D
      auto hScan = std::make_unique<TH1D>((name + "_full").c_str(), (name + "_full").c_str(), n_points, lower, upper);
      hScan->SetTitle((std::string("2LLH_full, ") + name + ";" + name + "; -2(ln L_{sample} + ln L_{xsec+flux} + ln L_{det})").c_str());
      hScan->SetDirectory(nullptr);
 
      auto hScanSam = std::make_unique<TH1D>((name + "_sam").c_str(), (name + "_sam").c_str(), n_points, lower, upper);
      hScanSam->SetTitle((std::string("2LLH_sam, ") + name + ";" + name + "; -2(ln L_{sample})").c_str());
      hScanSam->SetDirectory(nullptr);
 
      std::vector<std::unique_ptr<TH1D>> hScanSample(samples.size());
      std::vector<double> nSamLLH(samples.size(), 0.0);
      for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
      {
        std::string NameTemp = samples[ivs]->GetName();
        hScanSample[ivs] = std::make_unique<TH1D>((name+"_"+NameTemp).c_str(), (name+"_" + NameTemp).c_str(), n_points, lower, upper);
        hScanSample[ivs]->SetDirectory(nullptr);
        hScanSample[ivs]->SetTitle(("2LLH_" + NameTemp + ", " + name + ";" + name + "; -2(ln L_{" + NameTemp +"})").c_str());
      }
 
      std::vector<std::unique_ptr<TH1D>> hScanCov(systematics.size());
      std::vector<double> nCovLLH(systematics.size(), 0.0);
      for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc )
      {
        std::string NameTemp = systematics[ivc]->GetName();
        NameTemp = NameTemp.substr(0, NameTemp.find("_cov"));
        hScanCov[ivc] = std::make_unique<TH1D>((name + "_" + NameTemp).c_str(), (name + "_" + NameTemp).c_str(), n_points, lower, upper);
        hScanCov[ivc]->SetDirectory(nullptr);
        hScanCov[ivc]->SetTitle(("2LLH_" + NameTemp + ", " + name + ";" + name + "; -2(ln L_{" + NameTemp +"})").c_str());
      }
 
      std::vector<std::unique_ptr<TH1D>> hScanSamSplit;
      std::vector<double> sampleSplitllh;
      if(PlotLLHScanBySample)
      {
        int SampleIterator = 0;
        hScanSamSplit.resize(TotalNSamples);
        sampleSplitllh.resize(TotalNSamples);
        for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
        {
          for(int is = 0; is < samples[ivs]->GetNSamples(); ++is )
          {
            auto histName = name + samples[ivs]->GetSampleTitle(is);
            auto histTitle = std::string("2LLH_sam, ") + name + ";" + name + "; -2(ln L_{sample})";
            hScanSamSplit[SampleIterator] = std::make_unique<TH1D>(histName.c_str(), histTitle.c_str(), n_points, lower, upper);
             hScanSamSplit[SampleIterator]->SetDirectory(nullptr);
            SampleIterator++;
          }
        }
      }
 
      // Scan over the parameter space
      for (int j = 0; j < n_points; ++j)
      {
        if (j % countwidth == 0)
          M3::Utils::PrintProgressBar(j, n_points);
 
        // For PCA we have to do it differently
        if (IsPCA) {
          cov->GetPCAHandler()->SetParPropPCA(i, hScan->GetBinCenter(j+1));
        } else {
          // Set the parameter
          cov->SetParProp(i, hScan->GetBinCenter(j+1));
        }
 
        // Reweight the MC
        for(unsigned int ivs = 0; ivs < samples.size(); ++ivs ){
          samples[ivs]->Reweight();
        }
        //Total LLH
        double totalllh = 0.;
 
        // Get the -log L likelihoods
        double samplellh = 0.;
 
        for(unsigned int ivs = 0; ivs < samples.size(); ++ivs ) {
          nSamLLH[ivs] = samples[ivs]->GetLikelihood();
          samplellh += nSamLLH[ivs];
        }
 
        for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc ) {
          nCovLLH[ivc] = systematics[ivc]->GetLikelihood();
          totalllh += nCovLLH[ivc];
        }
 
        totalllh += samplellh;
 
        if(PlotLLHScanBySample)
        {
          int SampleIterator = 0;
          for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
          {
            for(int is = 0; is < samples[ivs]->GetNSamples(); ++is)
            {
              sampleSplitllh[SampleIterator] = samples[ivs]->GetSampleLikelihood(is);
              SampleIterator++;
            }
          }
        }
 
        for(unsigned int ivs = 0; ivs < samples.size(); ++ivs ) {
          hScanSample[ivs]->SetBinContent(j+1, 2*nSamLLH[ivs]);
        }
        for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc ) {
          hScanCov[ivc]->SetBinContent(j+1, 2*nCovLLH[ivc]);
        }
 
        hScanSam->SetBinContent(j+1, 2*samplellh);
        hScan->SetBinContent(j+1, 2*totalllh);
 
        if(PlotLLHScanBySample)
        {
          int SampleIterator = 0;
          for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
          {
            for(int is = 0; is < samples[ivs]->GetNSamples(); ++is)
            {
              hScanSamSplit[SampleIterator]->SetBinContent(j+1, 2*sampleSplitllh[SampleIterator]);
              SampleIterator++;
            }
          }
        }
      }
      for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc )
      {
        Cov_LLH[ivc]->cd();
        hScanCov[ivc]->Write();
      }
 
      for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
      {
        SampleClass_LLH[ivs]->cd();
        hScanSample[ivs]->Write();
      }
      Sample_LLH->cd();
      hScanSam->Write();
      Total_LLH->cd();
      hScan->Write();
 
      if(PlotLLHScanBySample)
      {
        int SampleIterator = 0;
        for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
        {
          for(int is = 0; is < samples[ivs]->GetNSamples(); ++is)
          {
            SampleSplit_LLH[SampleIterator]->cd();
            hScanSamSplit[SampleIterator]->Write();
            SampleIterator++;
          }
        }
      }
 
      // Reset the parameters to their CentralValue central values
      if (IsPCA) {
        cov->GetPCAHandler()->SetParPropPCA(i, CentralValue);
      } else {
        cov->SetParProp(i, CentralValue);
      }
    }//end loop over systematics
  }//end loop covariance classes
 
  for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc )
  {
    Cov_LLH[ivc]->Write();
    delete Cov_LLH[ivc];
  }
 
  for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
  {
    SampleClass_LLH[ivs]->Write();
    delete SampleClass_LLH[ivs];
  }
 
  Sample_LLH->Write();
  delete Sample_LLH;
 
  Total_LLH->Write();
  delete Total_LLH;
 
  if(PlotLLHScanBySample)
  {
    int SampleIterator = 0;
    for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
    {
      for(int is = 0; is < samples[ivs]->GetNSamples(); ++is )
      {
        SampleSplit_LLH[SampleIterator]->Write();
        delete SampleSplit_LLH[SampleIterator];
        SampleIterator++;
      }
    }
  }
}

RunMCMC()

virtual void FitterBase::RunMCMC ( )

pure virtual

The specific fitting algorithm implemented in this function depends on the derived class. It could be Markov Chain Monte Carlo (MCMC), MinuitFit, or another algorithm.

Implemented in MCMCBase, PyLikelihoodFit, PyFitterBase, NuDockServerBase, PredictiveThrower, PSO, and MinuitFit.

RunSigmaVar()

void FitterBase::RunSigmaVar

(

)

Perform a 1D/2D sigma var for all samples.

Apply custom range to make easier comparison with p-theta

Definition at line 1395 of file FitterBase.cpp.

                             {
// *************************
  // Save the settings into the output file
  SaveSettings();
 
  bool plot_by_mode = GetFromManager<bool>(fitMan->raw()["SigmaVar"]["PlotByMode"], false);
  bool plot_by_channel = GetFromManager<bool>(fitMan->raw()["SigmaVar"]["PlotByChannel"], false);
  auto SkipVector = GetFromManager<std::vector<std::string>>(fitMan->raw()["SigmaVar"]["SkipVector"], {}, __FILE__ , __LINE__);
 
  if (plot_by_mode) MACH3LOG_INFO("Plotting by sample and mode");
  if (plot_by_channel) MACH3LOG_INFO("Plotting by sample and channel");
  if (!plot_by_mode && !plot_by_channel) MACH3LOG_INFO("Plotting by sample only");
  if (plot_by_mode && plot_by_channel) MACH3LOG_INFO("Plotting by sample, mode and channel");
 
  auto SigmaArray = GetFromManager<std::vector<double>>(fitMan->raw()["SigmaVar"]["SigmaArray"], {-3, -1, 0, 1, 3}, __FILE__ , __LINE__);
  if (std::find(SigmaArray.begin(), SigmaArray.end(), 0.0) == SigmaArray.end()) {
    MACH3LOG_ERROR(":: SigmaArray does not contain 0! Current contents: {} ::", fmt::join(SigmaArray, ", "));
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  TDirectory* SigmaDir = outputFile->mkdir("SigmaVar");
  outputFile->cd();
 
  for (size_t s = 0; s < systematics.size(); ++s)
  {
    for(int i = 0; i < systematics[s]->GetNumParams(); i++)
    {
      std::string ParName = systematics[s]->GetParFancyName(i);
      // KS: Check if we want to skip this parameter
      if(CheckSkipParameter(SkipVector, ParName)) continue;
 
      MACH3LOG_INFO(":: Param {} ::", systematics[s]->GetParFancyName(i));
 
      TDirectory* ParamDir = SigmaDir->mkdir(ParName.c_str());
      ParamDir->cd();
 
      const double ParamCentralValue = systematics[s]->GetParProp(i);
      const double Prior = systematics[s]->GetParInit(i);
      const double ParamLower = systematics[s]->GetLowerBound(i);
      const double ParamUpper = systematics[s]->GetUpperBound(i);
 
      if (std::abs(ParamCentralValue - Prior) > 1e-10) {
        MACH3LOG_INFO("For {} scanning around value {} rather than prior {}", ParName, ParamCentralValue, Prior);
      }
 
      for(unsigned int iSample = 0; iSample < samples.size(); ++iSample)
      {
        auto* MaCh3Sample = samples[iSample];
        std::vector<TDirectory*> SampleDir(MaCh3Sample->GetNSamples());
        for (int SampleIndex = 0; SampleIndex < MaCh3Sample->GetNSamples(); ++SampleIndex) {
          SampleDir[SampleIndex] = ParamDir->mkdir(MaCh3Sample->GetSampleTitle(SampleIndex).c_str());
        }
 
        for (size_t j = 0; j < SigmaArray.size(); ++j) {
          double sigma = SigmaArray[j];
 
          double ParamShiftValue = ParamCentralValue + sigma * std::sqrt((*systematics[s]->GetCovMatrix())(i,i));
          ParamShiftValue = std::max(std::min(ParamShiftValue, ParamUpper), ParamLower);
 
          CustomRange(ParName, sigma, ParamShiftValue);
 
          MACH3LOG_INFO("  - set to {:<5.2f} ({:<2} sigma shift)", ParamShiftValue, sigma);
          systematics[s]->SetParProp(i, ParamShiftValue);
 
          std::ostringstream valStream;
          valStream << std::fixed << std::setprecision(2) << ParamShiftValue;
          std::string valueStr = valStream.str();
 
          std::ostringstream sigmaStream;
          sigmaStream << std::fixed << std::setprecision(2) << std::abs(sigma);
          std::string sigmaStr = sigmaStream.str();
 
          std::string suffix;
          if (sigma == 0) {
            suffix = "_" + ParName + "_nom_val_" + valueStr;
          } else {
            std::string sign = (sigma > 0) ? "p" : "n";
            suffix = "_" + ParName + "_sig_" + sign + sigmaStr + "_val_" + valueStr;
          }
 
          systematics[s]->SetParProp(i, ParamShiftValue);
          MaCh3Sample->Reweight();
 
          WriteHistogramsByMode(MaCh3Sample, suffix, plot_by_mode, plot_by_channel, SampleDir);
        }
        for (int subSampleIndex = 0; subSampleIndex < MaCh3Sample->GetNSamples(); ++subSampleIndex) {
          SampleDir[subSampleIndex]->Close();
          delete SampleDir[subSampleIndex];
        }
        ParamDir->cd();
      }
 
      systematics[s]->SetParProp(i, ParamCentralValue);
      MACH3LOG_INFO("  - set back to CV {:<5.2f}", ParamCentralValue);
      MACH3LOG_INFO("");
      ParamDir->Close();
      delete ParamDir;
      SigmaDir->cd();
    } // end loop over params
  } // end loop over systemics
  SigmaDir->Close();
  delete SigmaDir;
 
  outputFile->cd();
}

SanitiseInputs()

void FitterBase::SanitiseInputs ( )

protected

Remove obsolete memory and make other checks before fit starts.

Todo:

consider expanding into ParmaterHandler and add more sanitisers

Definition at line 221 of file FitterBase.cpp.

                                {
// *******************
  for (size_t i = 0; i < samples.size(); ++i) {
    samples[i]->CleanMemoryBeforeFit();
  }
}

SaveOutput()

void FitterBase::SaveOutput ( )

protected

Save output and close files.

Definition at line 229 of file FitterBase.cpp.

                            {
// *******************
  if(FileSaved) return;
  //Stop Clock
  clock->Stop();
 
  //KS: Some version of ROOT keep spamming about accessing already deleted object which is wrong and not helpful...
  int originalErrorLevel = gErrorIgnoreLevel;
  gErrorIgnoreLevel = kFatal;
 
  outputFile->cd();
  outTree->Write();
 
  MACH3LOG_INFO("{} steps took {:.2e} seconds to complete. ({:.2e}s / step).", step - stepStart, clock->RealTime(), clock->RealTime() / static_cast<double>(step - stepStart));
  MACH3LOG_INFO("{} steps were accepted.", accCount);
  #ifdef MACH3_DEBUG
  if (debug)
  {
    debugFile << "\n\n" << step << " steps took " << clock->RealTime() << " seconds to complete. (" << clock->RealTime() / step << "s / step).\n" << accCount<< " steps were accepted." << std::endl;
    debugFile.close();
  }
  #endif
 
  outputFile->Close();
  FileSaved = true;
 
  gErrorIgnoreLevel = originalErrorLevel;
}

SaveSettings()

void FitterBase::SaveSettings ( )

protected

Save the settings that the MCMC was run with.

Definition at line 77 of file FitterBase.cpp.

                              {
// *******************
  if(SettingsSaved) return;
 
  outputFile->cd();
 
  TDirectory* MaCh3Version = outputFile->mkdir("MaCh3Engine");
  MaCh3Version->cd();
 
  if (std::getenv("MaCh3_ROOT") == nullptr) {
    MACH3LOG_ERROR("Need MaCh3_ROOT environment variable");
    MACH3LOG_ERROR("Please remember about source bin/setup.MaCh3.sh");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  if (std::getenv("MACH3") == nullptr) {
    MACH3LOG_ERROR("Need MACH3 environment variable");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  std::string header_path = std::string(std::getenv("MACH3"));
  header_path += "/version.h";
  FILE* file = fopen(header_path.c_str(), "r");
  //KS: It is better to use experiment specific header file. If given experiment didn't provide it we gonna use one given by Core MaCh3.
  if (!file) {
    header_path = std::string(std::getenv("MaCh3_ROOT"));
    header_path += "/version.h";
  } else {
    fclose(file);
  }
 
  // EM: embed the cmake generated version.h file
  TMacro versionHeader("version_header", "version_header");
  versionHeader.ReadFile(header_path.c_str());
  versionHeader.Write();
 
  if(GetName() == ""){
    MACH3LOG_ERROR("Name of currently used algorithm is {}", GetName());
    MACH3LOG_ERROR("Have you forgotten to modify AlgorithmName?");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  TNamed Engine(GetName(), GetName());
  Engine.Write(GetName().c_str());
 
  MaCh3Version->Write();
  delete MaCh3Version;
 
  outputFile->cd();
 
  fitMan->SaveSettings(outputFile);
 
  MACH3LOG_WARN("\033[0;31mCurrent Total RAM usage is {:.2f} GB\033[0m", M3::Utils::getValue("VmRSS") / 1048576.0);
  MACH3LOG_WARN("\033[0;31mOut of Total available RAM {:.2f} GB\033[0m", M3::Utils::getValue("MemTotal") / 1048576.0);
 
  MACH3LOG_INFO("#####Current Setup#####");
  MACH3LOG_INFO("Number of covariances: {}", systematics.size());
  for(unsigned int i = 0; i < systematics.size(); ++i)
    MACH3LOG_INFO("{}: Cov name: {}, it has {} params", i, systematics[i]->GetName(), systematics[i]->GetNumParams());
  MACH3LOG_INFO("Number of SampleHandlers: {}", samples.size());
  for(unsigned int i = 0; i < samples.size(); ++i) {
    MACH3LOG_INFO("{}: SampleHandler name: {}, it has {} samples",i , samples[i]->GetName(), samples[i]->GetNSamples());
    for(int iSam = 0; iSam < samples[i]->GetNSamples(); ++iSam) {
      MACH3LOG_INFO("   {}: Sample name: {}, with {} osc channels",iSam , samples[i]->GetSampleTitle(iSam), samples[i]->GetNOscChannels(iSam));
    }
  }
  //TN: Have to close the folder in order to write it to disk before SaveOutput is called in the destructor
  CovFolder->Close();
  SampleFolder->Close();
 
  SettingsSaved = true;
}

StartFromPreviousFit()

void FitterBase::StartFromPreviousFit ( const std::string &  FitName )

virtual

Allow to start from previous fit/chain.

Parameters

FitName

Name of previous chain

Reimplemented in MCMCBase.

Definition at line 346 of file FitterBase.cpp.

                                                              {
// *******************
  MACH3LOG_INFO("Getting starting position from {}", FitName);
  TFile *infile = M3::Open(FitName, "READ", __FILE__, __LINE__);
  TTree *posts = infile->Get<TTree>("posteriors");
  unsigned int step_val = 0;
  double log_val = M3::_LARGE_LOGL_;
  posts->SetBranchAddress("step",&step_val);
  posts->SetBranchAddress("LogL",&log_val);
 
  for (size_t s = 0; s < systematics.size(); ++s)
  {
    TDirectory* CovarianceFolder = infile->Get<TDirectory>("CovarianceFolder");
 
    std::string ConfigName = "Config_" + systematics[s]->GetName();
    TMacro *ConfigCov = CovarianceFolder->Get<TMacro>(ConfigName.c_str());
    // KS: Not every covariance uses yaml, if it uses yaml make sure they are identical
    if (ConfigCov != nullptr) {
      // Config which was in MCMC from which we are starting
      YAML::Node CovSettings = TMacroToYAML(*ConfigCov);
      // Config from currently used cov object
      YAML::Node ConfigCurrent = systematics[s]->GetConfig();
 
      if (!compareYAMLNodes(CovSettings, ConfigCurrent))
      {
        MACH3LOG_ERROR("Yaml configs in previous chain (from path {}) and current one are different", FitName);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
      delete ConfigCov;
    }
 
    CovarianceFolder->Close();
    delete CovarianceFolder;
 
    std::vector<double> branch_vals;
    std::vector<std::string> branch_name;
    systematics[s]->MatchMaCh3OutputBranches(posts, branch_vals, branch_name);
    posts->GetEntry(posts->GetEntries()-1);
 
    systematics[s]->SetParameters(branch_vals);
    systematics[s]->AcceptStep();
 
    MACH3LOG_INFO("Printing new starting values for: {}", systematics[s]->GetName());
    systematics[s]->PrintPreFitCurrPropValues();
 
    // Resetting branch addressed to nullptr as we don't want to write into a delected vector out of scope...
    for (int i = 0; i < systematics[s]->GetNumParams(); ++i) {
      posts->SetBranchAddress(systematics[s]->GetParName(i).c_str(), nullptr);
    }
  }
  logLCurr = log_val;
  logLProp = log_val;
  delete posts;
  infile->Close();
  delete infile;
 
  for (size_t s = 0; s < systematics.size(); ++s) {
    if(systematics[s]->GetDoAdaption()){ //Use separate throw matrix for xsec
      systematics[s]->SetNumberOfSteps(step_val);
    }
  }
}

Member Data Documentation

accCount

int FitterBase::accCount

protected

counts accepted steps

Definition at line 121 of file FitterBase.h.

accProb

double FitterBase::accProb

protected

current acceptance prob

Definition at line 119 of file FitterBase.h.

AlgorithmName

std::string FitterBase::AlgorithmName

protected

Name of fitting algorithm that is being used.

Definition at line 170 of file FitterBase.h.

auto_save

int FitterBase::auto_save

protected

auto save every N steps

Definition at line 157 of file FitterBase.h.

clock

std::unique_ptr<TStopwatch> FitterBase::clock

protected

tells global time how long fit took

Definition at line 139 of file FitterBase.h.

CovFolder

TDirectory* FitterBase::CovFolder

protected

Output cov folder.

Definition at line 151 of file FitterBase.h.

FileSaved

bool FitterBase::FileSaved

protected

Checks if file saved not repeat some operations.

Definition at line 163 of file FitterBase.h.

fitMan

Manager* FitterBase::fitMan

protected

The manager for configuration handling.

Definition at line 110 of file FitterBase.h.

fTestLikelihood

bool FitterBase::fTestLikelihood

protected

Necessary for some fitting algorithms like PSO.

Definition at line 160 of file FitterBase.h.

logLCurr

double FitterBase::logLCurr

protected

current likelihood

Definition at line 115 of file FitterBase.h.

logLProp

double FitterBase::logLProp

protected

proposed likelihood

Definition at line 117 of file FitterBase.h.

outputFile

TFile* FitterBase::outputFile

protected

Output.

Definition at line 149 of file FitterBase.h.

OutputPrepared

bool FitterBase::OutputPrepared

protected

Checks if output prepared not repeat some operations.

Definition at line 167 of file FitterBase.h.

outTree

TTree* FitterBase::outTree

protected

Output tree with posteriors.

Definition at line 155 of file FitterBase.h.

random

std::unique_ptr<TRandom3> FitterBase::random

protected

Random number.

Definition at line 146 of file FitterBase.h.

sample_llh

std::vector<double> FitterBase::sample_llh

protected

store the llh breakdowns

Definition at line 126 of file FitterBase.h.

SampleFolder

TDirectory* FitterBase::SampleFolder

protected

Output sample folder.

Definition at line 153 of file FitterBase.h.

samples

std::vector<SampleHandlerInterface*> FitterBase::samples

protected

Sample holder.

Definition at line 131 of file FitterBase.h.

SettingsSaved

bool FitterBase::SettingsSaved

protected

Checks if setting saved not repeat some operations.

Definition at line 165 of file FitterBase.h.

step

unsigned int FitterBase::step

protected

current state

Definition at line 113 of file FitterBase.h.

stepClock

std::unique_ptr<TStopwatch> FitterBase::stepClock

protected

tells how long single step/fit iteration took

Definition at line 141 of file FitterBase.h.

stepStart

unsigned int FitterBase::stepStart

protected

step start, by default 0 if we start from previous chain then it will be different

Definition at line 123 of file FitterBase.h.

stepTime

double FitterBase::stepTime

protected

Time of single step.

Definition at line 143 of file FitterBase.h.

syst_llh

std::vector<double> FitterBase::syst_llh

protected

systematic llh breakdowns

Definition at line 128 of file FitterBase.h.

systematics

std::vector<ParameterHandlerBase*> FitterBase::systematics

protected

Systematic holder.

Definition at line 136 of file FitterBase.h.

TotalNSamples

unsigned int FitterBase::TotalNSamples

protected

Total number of samples used, single SampleHandler can store more than one analysis sample!

Definition at line 133 of file FitterBase.h.

---

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

• Fitters/FitterBase.h
• Fitters/FitterBase.cpp