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.

virtual	~FitterBase ()
	Destructor for the FitterBase class.

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.

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.

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.

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.

void	RunLLHScan ()
	Perform a 1D likelihood scan.

void	GetStepScaleBasedOnLLHScan ()
	LLH scan is good first estimate of step scale.

void	Run2DLLHScan ()
	Perform a 2D likelihood scan.

void	RunSigmaVar ()
	Perform a 2D and 1D sigma var for all samples.

virtual void	StartFromPreviousFit (const std::string &FitName)
	Allow to start from previous fit/chain.

virtual std::string	GetName () const
	Get name of class.

Protected Member Functions
void	ProcessMCMC ()
	Process MCMC output.

void	PrepareOutput ()
	Prepare the output file.

void	SaveOutput ()
	Save output and close files.

void	SanitiseInputs ()
	Remove obsolete memory and make other checks before fit starts.

void	SaveSettings ()
	Save the settings that the MCMC was run with.

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.

bool	CheckSkipParameter (const std::vector< std::string > &SkipVector, const std::string &ParamName) const
	KS: Check whether we want to skip parameter using skip vector.

Protected Attributes
manager *	fitMan
	The manager.

unsigned int	step
	current state

double	logLCurr
	current likelihood

double	logLProp
	proposed likelihood

double	accProb
	current acceptance prob

int	accCount
	counts accepted steps

int	stepStart
	step start if restarting

std::vector< double >	sample_llh
	store the llh breakdowns

std::vector< double >	syst_llh
	systematic llh breakdowns

std::vector< SampleHandlerBase * >	samples
	Sample holder.

unsigned int	TotalNSamples
	Total number of samples used.

std::vector< ParameterHandlerBase * >	systematics
	Systematic holder.

std::unique_ptr< TStopwatch >	clock
	tells global time how long fit took

std::unique_ptr< TStopwatch >	stepClock
	tells how long single step/fit iteration took

double	stepTime
	Time of single step.

std::unique_ptr< TRandom3 >	random
	Random number.

TFile *	outputFile
	Output.

TDirectory *	CovFolder
	Output cov folder.

TDirectory *	SampleFolder
	Output sample folder.

TTree *	outTree
	Output tree with posteriors.

int	auto_save
	auto save every N steps

bool	fTestLikelihood
	Necessary for some fitting algorithms like PSO.

bool	FileSaved
	Checks if file saved not repeat some operations.

bool	SettingsSaved
	Checks if setting saved not repeat some operations.

bool	OutputPrepared
	Checks if output prepared not repeat some operations.

Detailed Description

Base class for implementing fitting algorithms.

See also: For more details, visit the Wiki.

Definition at line 23 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) {
// *************************
  //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 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__);
 
  #ifdef MULTITHREAD
  //KS: TODO This should help with performance when saving entries to ROOT file. I didn't have time to validate hence commented out
  //Based on other tests it is really helpful
  //ROOT::EnableImplicitMT();
  #endif
  // 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 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 71 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 ( 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.

Parameters

sample A pointer to a sample PDF object derived from ParameterHandlerBase.

Definition at line 246 of file FitterBase.cpp.

                                                                  {
// *************************
  //Check if the sample has a unique name
  for (const auto &s : samples) {
    if (s->GetTitle() == sample->GetTitle()) {
      MACH3LOG_ERROR("SampleHandler with name '{}' already exists!", sample->GetTitle());
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
  }
  // Save additional info from samples
  SampleFolder->cd();
 
  sample->SaveAdditionalInfo(SampleFolder);
  TotalNSamples += sample->GetNsamples();
  MACH3LOG_INFO("Adding {} object, with {} samples", sample->GetTitle(), 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 267 of file FitterBase.cpp.

                                                            {
// *************************
  MACH3LOG_INFO("Adding systematic object {}, with {} params", cov->GetName(), cov->GetNumParams());
  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 500 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;
}

◆ 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 418 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]->GetTitle());
    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]->GetTitle());
    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()

virtual std::string FitterBase::GetName ( ) const

inlinevirtual

Get name of class.

Reimplemented in LikelihoodFit, mcmc, MinuitFit, PSO, and PyFitterBase.

Definition at line 66 of file FitterBase.h.

66{return "FitterBase";};

◆ GetScaneRange()

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

protected

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

Definition at line 480 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 ( )

LLH scan is good first estimate of step scale.

Definition at line 826 of file FitterBase.cpp.

                                            {
// *************************
  TDirectory *Sample_LLH = outputFile->Get<TDirectory>("Sample_LLH");
  MACH3LOG_INFO("Starting Get Step Scale Based On LLHScan");
 
  if(Sample_LLH->IsZombie())
  {
    MACH3LOG_WARN("Couldn't find Sample_LLH, it looks like LLH scan wasn't run, will do this now");
    RunLLHScan();
  }
 
  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 = TMath::Abs(Var)/std::sqrt(LLH_val);
 
      // 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 NewStepScale = approxSigma * 2.38/std::sqrt(npars);
      StepScale[i] = NewStepScale;
      MACH3LOG_DEBUG("Sigma: {}", approxSigma);
      MACH3LOG_DEBUG("optimal Step Size: {}", NewStepScale);
    }
    cov->SetIndivStepScale(StepScale);
    cov->SaveUpdatedMatrixConfig();
  }
}

◆ PrepareOutput()

void FitterBase::PrepareOutput ( )

protected

Prepare the output file.

Definition at line 148 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__);
 
    // Prepare the output trees
    for (ParameterHandlerBase *cov : systematics) {
      cov->SetBranches(*outTree, SaveProposal);
    }
 
    outTree->Branch("LogL", &logLCurr, "LogL/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 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 371 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"]["Output"]["Filename"], __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 875 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));
 
  std::map<std::string, std::vector<double>> scanRanges;
  const bool isScanRanges = GetScaneRange(scanRanges);
 
  const double nSigma = GetFromManager<int>(fitMan->raw()["LLHScan"]["LLHScanSigma"], 1., __FILE__, __LINE__);
 
  // 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 priors and bounds
      double prior_x = cov->GetParInit(i);
      if (IsPCA) prior_x = cov->GetPCAHandler()->GetParCurrPCA(i);
 
      // Get the covariance matrix and do the +/- nSigma
      // Set lower and upper bounds relative the prior
      double lower_x = prior_x - nSigma*cov->GetDiagonalError(i);
      double upper_x = prior_x + nSigma*cov->GetDiagonalError(i);
      // If PCA, transform these parameter values to the PCA basis
      if (IsPCA) {
        lower_x = prior_x - nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
        upper_x = prior_x + nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
        MACH3LOG_INFO("eval {} = {:.2f}", i, cov->GetPCAHandler()->GetEigenValues()(i));
        MACH3LOG_INFO("prior {} = {:.2f}", i, prior_x);
        MACH3LOG_INFO("lower {} = {:.2f}", i, lower_x);
        MACH3LOG_INFO("upper {} = {:.2f}", i, upper_x);
        MACH3LOG_INFO("nSigma = {:.2f}", nSigma);
      }
      // 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_x);
        if (it != scanRanges.end() && it->second.size() == 2) { //Making sure the range is has only two entries
          lower_x = it->second[0];
          upper_x = it->second[1];
          MACH3LOG_INFO("Found matching param name for setting specified range for {}", name_x);
          MACH3LOG_INFO("Range for {} = [{:.2f}, {:.2f}]", name_x, lower_x, upper_x);
        }
      }
 
      // Cross-section and flux parameters have boundaries that we scan between, check that these are respected in setting lower and upper variables
      lower_x = std::max(lower_x, cov->GetLowerBound(i));
      upper_x = std::min(upper_x, cov->GetUpperBound(i));
      // 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 priors and bounds
        double prior_y = cov->GetParInit(j);
        if (IsPCA) prior_y = cov->GetPCAHandler()->GetParCurrPCA(j);
 
        // Set lower and upper bounds relative the prior
        double lower_y = prior_y - nSigma*cov->GetDiagonalError(j);
        double upper_y = prior_y + nSigma*cov->GetDiagonalError(j);
        // If PCA, transform these parameter values to the PCA basis
        if (IsPCA) {
          lower_y = prior_y - nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(j));
          upper_y = prior_y + nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(j));
          MACH3LOG_INFO("eval {} = {:.2f}", i, cov->GetPCAHandler()->GetEigenValues()(j));
          MACH3LOG_INFO("prior {} = {:.2f}", i, prior_y);
          MACH3LOG_INFO("lower {} = {:.2f}", i, lower_y);
          MACH3LOG_INFO("upper {} = {:.2f}", i, upper_y);
          MACH3LOG_INFO("nSigma = {:.2f}", nSigma);
        }
        // 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_y);
          if (it != scanRanges.end() && it->second.size() == 2) { //Making sure the range is has only two entries
            lower_y = it->second[0];
            upper_y = it->second[1];
            MACH3LOG_INFO("Found matching param name for setting specified range for {}", name_y);
            MACH3LOG_INFO("Range for {} = [{:.2f}, {:.2f}]", name_y, lower_y, upper_y);
          }
        }
 
        // Cross-section and flux parameters have boundaries that we scan between, check that these are respected in setting lower and upper variables
        lower_y = std::max(lower_y, cov->GetLowerBound(j));
        upper_y = std::min(upper_y, cov->GetUpperBound(j));
        MACH3LOG_INFO("Scanning X {} with {} steps, from {:.2f} - {:.2f}, prior = {}", name_x, n_points, lower_x, upper_x, prior_x);
        MACH3LOG_INFO("Scanning Y {} with {} steps, from {:.2f} - {:.2f}, prior = {}", name_y, n_points, lower_y, upper_y, prior_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)
            MaCh3Utils::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 prior central values
        if (IsPCA) {
          cov->GetPCAHandler()->SetParPropPCA(i, prior_x);
          cov->GetPCAHandler()->SetParPropPCA(j, prior_y);
        } else {
          cov->SetParProp(i, prior_x);
          cov->SetParProp(j, prior_y);
        }
      } //end loop over systematics y
    }//end loop over systematics X
  }//end loop covariance classes
  Sample_2DLLH->Write();
  delete Sample_2DLLH;
}

◆ RunLLHScan()

void FitterBase::RunLLHScan ( )

Perform a 1D likelihood scan.

Definition at line 516 of file FitterBase.cpp.

                            {
// *************************
  // Save the settings into the output file
  SaveSettings();
 
  MACH3LOG_INFO("Starting LLH Scan");
 
  //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]->GetTitle();
    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]->GetSampleName(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__);
  double nSigma = GetFromManager<int>(fitMan->raw()["LLHScan"]["LLHScanSigma"], 1., __FILE__, __LINE__);
 
  // We print 5 reweights
  const int countwidth = int(double(n_points)/double(5));
 
  // 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 = GetScaneRange(scanRanges);
   
  // 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 priors and bounds
      double prior = cov->GetParInit(i);
      if (IsPCA) prior = cov->GetPCAHandler()->GetParCurrPCA(i);
 
      // Get the covariance matrix and do the +/- nSigma
      // Set lower and upper bounds relative the prior
      // Set the parameter ranges between which LLH points are scanned
      double lower = prior - nSigma*cov->GetDiagonalError(i);
      double upper = prior + nSigma*cov->GetDiagonalError(i);
      // If PCA, transform these parameter values to the PCA basis
      if (IsPCA) {
        lower = prior - nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
        upper = prior + nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
        MACH3LOG_INFO("eval {} = {:.2f}", i, cov->GetPCAHandler()->GetEigenValues()(i));
        MACH3LOG_INFO("prior {} = {:.2f}", i, prior);
        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}], prior = {:.2f}", name, n_points, lower, upper, prior);
 
      // 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());
      for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
      {
        std::string NameTemp = samples[ivs]->GetTitle();
        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());
        nSamLLH[ivs] = 0.;
      }
 
      std::vector<std::unique_ptr<TH1D>> hScanCov(systematics.size());
      std::vector<double> nCovLLH(systematics.size());
      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());
        nCovLLH[ivc] = 0.;
      }
 
      std::vector<TH1D*> hScanSamSplit;
      std::vector<double> sampleSplitllh;
      if(PlotLLHScanBySample)
      {
        int SampleIterator = 0;
        for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
        {
          hScanSamSplit.resize(TotalNSamples);
          sampleSplitllh.resize(TotalNSamples);
          for(int is = 0; is < samples[ivs]->GetNsamples(); ++is )
          {
            hScanSamSplit[SampleIterator] = new TH1D((name+samples[ivs]->GetSampleName(is)).c_str(), (name+samples[ivs]->GetSampleName(is)).c_str(), n_points, lower, upper);
            hScanSamSplit[SampleIterator]->SetTitle((std::string("2LLH_sam, ") + name + ";" + name + "; -2(ln L_{sample})").c_str());
            SampleIterator++;
          }
        }
      }
 
      // Scan over the parameter space
      for (int j = 0; j < n_points; ++j)
      {
        if (j % countwidth == 0)
          MaCh3Utils::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();
            delete hScanSamSplit[SampleIterator];
            SampleIterator++;
          }
        }
      }
 
      // Reset the parameters to their prior central values
      if (IsPCA) {
        cov->GetPCAHandler()->SetParPropPCA(i, prior);
      } else {
        cov->SetParProp(i, prior);
      }
    }//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 mcmc, MinuitFit, PSO, PyFitterBase, PyLikelihoodFit, and LikelihoodFit.

◆ RunSigmaVar()

void FitterBase::RunSigmaVar ( )

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

Warning: Code uses TH2Poly

Definition at line 1043 of file FitterBase.cpp.

                             {
// *************************
  // Save the settings into the output file
  SaveSettings();
 
  MACH3LOG_INFO("Starting Sigma Variation");
 
  // Number of variations we want
  constexpr int numVar = 5;
  //-3 -1 0 +1 +3 sigma variation
  constexpr int sigmaArray[numVar] = {-3, -1, 0, 1, 3};
 
  outputFile->cd();
 
  //KS: this is only relevant if PlotByMode is turned on
  //Checking each mode is time consuming so we only consider one which are relevant for particular analysis
  constexpr int nRelevantModes = 2;
 
  bool DoByMode = GetFromManager<int>(fitMan->raw()["General"]["DoByMode"], false, __FILE__ , __LINE__);
 
  //KS: If true it will make additional plots with LLH sample contribution in each bin, should make it via config file...
  bool PlotLLHperBin = false;
  auto SkipVector = GetFromManager<std::vector<std::string>>(fitMan->raw()["LLHScan"]["LLHScanSkipVector"], {}, __FILE__ , __LINE__);;
 
  for (ParameterHandlerBase *cov : systematics)
  {
    TMatrixDSym *Cov = cov->GetCovMatrix();
 
    if(cov->IsPCA())
    {
      MACH3LOG_ERROR("Using PCAed matrix not implemented within sigma var code, I am sorry :(");
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
 
    // Loop over xsec parameters
    for (int i = 0; i < cov->GetNumParams(); ++i)
    {
      // Get the parameter name
      std::string name = cov->GetParFancyName(i);
      // KS: Check if we want to skip this parameter
      if(CheckSkipParameter(SkipVector, name)) continue;
 
      outputFile->cd();
      TDirectory* dirArryDial = outputFile->mkdir(name.c_str());
      std::vector<TDirectory*> dirArrySample(TotalNSamples);
 
      int SampleIterator = 0;
      // Get each sample and how it's responded to our reweighted parameter
      for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
      {
        for(int k = 0; k < samples[ivs]->GetNsamples(); k++ )
        {
          std::string title = std::string(samples[ivs]->GetPDF(k)->GetName());
          dirArryDial->cd();
          dirArrySample[SampleIterator] = dirArryDial->mkdir(title.c_str());
          SampleIterator++;
        }
      }
 
      // Get the initial value of ith parameter
      double init = cov->GetParInit(i);
 
      std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_x(numVar);
      std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_y(numVar);
      std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_x_norm(numVar);
      std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_y_norm(numVar);
 
      // Set up for single mode
      TH1D ****sigmaArray_mode_x = nullptr;
      TH1D ****sigmaArray_mode_y = nullptr;
      if (DoByMode)
      {
        sigmaArray_mode_x = new TH1D***[numVar]();
        sigmaArray_mode_y = new TH1D***[numVar]();
      }
 
      MACH3LOG_INFO("{:<20}{}", name, "");
 
      // Make asymmetric errors just in case
      for (int j = 0; j < numVar; ++j)
      {
        // New value = prior + variation*1sigma uncertainty
        double paramVal = cov->GetParInit(i)+sigmaArray[j]*std::sqrt((*Cov)(i,i));
 
        // Check the bounds on the parameter
        paramVal = std::max(cov->GetLowerBound(i), std::min(paramVal, cov->GetUpperBound(i)));
 
        // Set the parameter
        cov->SetParProp(i, paramVal);
        // And reweight the sample
        for(unsigned int ivs = 0; ivs < samples.size(); ivs++) {
          samples[ivs]->Reweight();
        }
 
        sigmaArray_x[j].resize(TotalNSamples);
        sigmaArray_y[j].resize(TotalNSamples);
        sigmaArray_x_norm[j].resize(TotalNSamples);
        sigmaArray_y_norm[j].resize(TotalNSamples);
 
        if (DoByMode)
        {
          sigmaArray_mode_x[j] = new TH1D**[TotalNSamples]();
          sigmaArray_mode_y[j] = new TH1D**[TotalNSamples]();
        }
 
        SampleIterator = 0;
        // Get each sample and how it's responded to our reweighted parameter
        for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
        {
          for (int k = 0; k < samples[ivs]->GetNsamples(); ++k)
          {
            // Make a string of the double
            std::ostringstream ss;
            ss << paramVal;
            std::string parVarTitle = name + "_" + ss.str();
 
            auto currSamp = M3::Clone<TH2Poly>(static_cast<TH2Poly*>(samples[ivs]->GetPDF(k)));
            // Set a descriptiv-ish title
            std::string title_long = std::string(currSamp->GetName())+"_"+parVarTitle;
 
            // Enable the mode histograms AFTER addSelection is called
            //Get the 1d binning we want. Let's just use SetupBinning to get this as it already exists
            std::vector<double> xbins;
            std::vector<double> ybins;
            samples[ivs]->SetupBinning(M3::int_t(k), xbins, ybins);
 
            //KS:here we loop over all reaction modes defined in "RelevantModes[nRelevantModes]"
            if (DoByMode)
            {
              //KS: this is only relevant if PlotByMode is turned on
              //Checking each mode is time consuming so we only consider one which are relevant for particular analysis
              MaCh3Modes_t RelevantModes[nRelevantModes] = {samples[ivs]->GetMaCh3Modes()->GetMode("CCQE"), samples[ivs]->GetMaCh3Modes()->GetMode("2p2h")};
 
              sigmaArray_mode_x[j][SampleIterator] = new TH1D*[nRelevantModes]();
              sigmaArray_mode_y[j][SampleIterator] = new TH1D*[nRelevantModes]();
              // Now get the TH2D mode variations
              std::string mode_title_long;
 
              for(int ir = 0; ir < nRelevantModes; ir++)
              {
                auto currSampMode = M3::Clone<TH2Poly>(static_cast<TH2Poly*>(samples[ivs]->GetPDFMode(k, RelevantModes[ir])));
 
                mode_title_long = title_long + "_" + samples[ivs]->GetMaCh3Modes()->GetMaCh3ModeName(RelevantModes[ir]);
                currSampMode->SetNameTitle(mode_title_long.c_str(), mode_title_long.c_str());
                dirArrySample[SampleIterator]->cd();
                currSampMode->Write();
 
                sigmaArray_mode_x[j][SampleIterator][ir] = PolyProjectionX(currSampMode.get(), (mode_title_long+"_xProj").c_str(), xbins);
                sigmaArray_mode_x[j][SampleIterator][ir]->SetDirectory(nullptr);
                sigmaArray_mode_y[j][SampleIterator][ir] = PolyProjectionY(currSampMode.get(), (mode_title_long+"_yProj").c_str(), ybins);
                sigmaArray_mode_y[j][SampleIterator][ir]->SetDirectory(nullptr);
              }
            }
 
            //KS: This will give different results depending if data or Asimov, both have their uses
            if (PlotLLHperBin)
            {
              auto currLLHSamp = M3::Clone<TH2Poly>(static_cast<TH2Poly*>(samples[ivs]->GetPDF(k)));
              currLLHSamp->Reset("");
              currLLHSamp->Fill(0.0, 0.0, 0.0);
 
              TH2Poly* MCpdf = static_cast<TH2Poly*>(samples[ivs]->GetPDF(k));
              TH2Poly* Datapdf = static_cast<TH2Poly*>(samples[ivs]->GetData(k));
              TH2Poly* W2pdf = samples[ivs]->GetW2(k);
 
              for(int bin = 1; bin < currLLHSamp->GetNumberOfBins()+1; bin++)
              {
                const double mc = MCpdf->GetBinContent(bin);
                const double dat = Datapdf->GetBinContent(bin);
                const double w2 = W2pdf->GetBinContent(bin);
                currLLHSamp->SetBinContent(bin, samples[ivs]->GetTestStatLLH(dat, mc, w2));
              }
              currLLHSamp->SetNameTitle((title_long+"_LLH").c_str() ,(title_long+"_LLH").c_str());
              dirArrySample[SampleIterator]->cd();
              currLLHSamp->Write();
            }
 
            // Project down onto x axis
            sigmaArray_x[j][SampleIterator] = std::unique_ptr<TH1D>(PolyProjectionX(currSamp.get(), (title_long+"_xProj").c_str(), xbins));
            sigmaArray_x[j][SampleIterator]->SetDirectory(nullptr);
            sigmaArray_x[j][SampleIterator]->GetXaxis()->SetTitle(currSamp->GetXaxis()->GetTitle());
            sigmaArray_y[j][SampleIterator] = std::unique_ptr<TH1D>(PolyProjectionY(currSamp.get(), (title_long+"_yProj").c_str(), ybins));
            sigmaArray_y[j][SampleIterator]->SetDirectory(nullptr);
            sigmaArray_y[j][SampleIterator]->GetXaxis()->SetTitle(currSamp->GetYaxis()->GetTitle());
 
            sigmaArray_x_norm[j][SampleIterator] = M3::Clone<TH1D>(sigmaArray_x[j][SampleIterator].get());
            sigmaArray_x_norm[j][SampleIterator]->Scale(1., "width");
            sigmaArray_y_norm[j][SampleIterator] = M3::Clone<TH1D>(sigmaArray_y[j][SampleIterator].get());
            sigmaArray_y_norm[j][SampleIterator]->SetDirectory(nullptr);
            sigmaArray_y_norm[j][SampleIterator]->Scale(1., "width");
 
            currSamp->SetNameTitle(title_long.c_str(), title_long.c_str());
            dirArrySample[k]->cd();
            currSamp->Write();
 
            sigmaArray_x[j][k]->Write();
            sigmaArray_y[j][k]->Write();
            SampleIterator++;
          }//End loop over samples
        }
      } // End looping over variation
 
      // Restore the parameter to prior value
      cov->SetParProp(i, init);
 
      SampleIterator = 0;
      // Get each sample and how it's responded to our reweighted parameter
      for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
      {
        for (int k = 0; k < samples[ivs]->GetNsamples(); ++k)
        {
          std::string title = std::string(samples[ivs]->GetPDF(k)->GetName()) + "_" + name;
          auto var_x = MakeAsymGraph(sigmaArray_x[1][SampleIterator].get(), sigmaArray_x[2][SampleIterator].get(), sigmaArray_x[3][SampleIterator].get(), (title+"_X").c_str());
          auto var_y = MakeAsymGraph(sigmaArray_y[1][SampleIterator].get(), sigmaArray_y[2][SampleIterator].get(), sigmaArray_y[3][SampleIterator].get(), (title+"_Y").c_str());
 
          auto var_x_norm = MakeAsymGraph(sigmaArray_x_norm[1][SampleIterator].get(), sigmaArray_x_norm[2][SampleIterator].get(), sigmaArray_x_norm[3][SampleIterator].get(), (title+"_X_norm").c_str());
          auto var_y_norm = MakeAsymGraph(sigmaArray_y_norm[1][SampleIterator].get(), sigmaArray_y_norm[2][SampleIterator].get(), sigmaArray_y_norm[3][SampleIterator].get(), (title+"_Y_norm").c_str());
 
          dirArrySample[SampleIterator]->cd();
          var_x->Write();
          var_y->Write();
          var_x_norm->Write();
          var_y_norm->Write();
 
          //KS: here we loop over all reaction modes defined in "RelevantModes[nRelevantModes]"
          if (DoByMode)
          {
            //KS: this is only relevant if PlotByMode is turned on
            //Checking each mode is time consuming so we only consider one which are relevant for particular analysis
            MaCh3Modes_t RelevantModes[nRelevantModes] = {samples[ivs]->GetMaCh3Modes()->GetMode("CCQE"), samples[ivs]->GetMaCh3Modes()->GetMode("2p2h")};
 
            for(int ir = 0; ir < nRelevantModes;ir++)
            {
              auto var_mode_x = MakeAsymGraph(sigmaArray_mode_x[1][SampleIterator][ir], sigmaArray_mode_x[2][SampleIterator][ir], sigmaArray_mode_x[3][SampleIterator][ir], (title+"_"+samples[ivs]->GetMaCh3Modes()->GetMaCh3ModeName(RelevantModes[ir])+"_X").c_str());
              auto var_mode_y = MakeAsymGraph(sigmaArray_mode_y[1][SampleIterator][ir], sigmaArray_mode_y[2][SampleIterator][ir], sigmaArray_mode_y[3][SampleIterator][ir], (title+"_"+samples[ivs]->GetMaCh3Modes()->GetMaCh3ModeName(RelevantModes[ir])+"_Y").c_str());
 
              dirArrySample[SampleIterator]->cd();
              var_mode_x->Write();
              var_mode_y->Write();
            } // end for nRelevantModes
          } // end if mode
 
          SampleIterator++;
        }//End loop over samples(k)
      }//end looping over sample object
      SampleIterator = 0;
      for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
      {
        for (int k = 0; k < samples[ivs]->GetNsamples(); ++k)
        {
          dirArrySample[SampleIterator]->Close();
          delete dirArrySample[SampleIterator];
          SampleIterator++;
        }
      }
 
      dirArryDial->Close();
      delete dirArryDial;
 
      if (DoByMode)
      {
        for (int j = 0; j < numVar; ++j)
        {
          SampleIterator = 0;
          for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
          {
            for (int k = 0; k < samples[ivs]->GetNsamples(); ++k)
            {
              for(int ir = 0; ir < nRelevantModes;ir++)
              {
                delete sigmaArray_mode_x[j][SampleIterator][ir];
                delete sigmaArray_mode_y[j][SampleIterator][ir];
              }// end for nRelevantModes
              delete[] sigmaArray_mode_x[j][SampleIterator];
              delete[] sigmaArray_mode_y[j][SampleIterator];
              SampleIterator++;
            }// end for samples
          }
        }
        delete[] sigmaArray_mode_x;
        delete[] sigmaArray_mode_y;
      }
    } // end looping over xsec parameters (i)
  } // end looping over covarianceBase objects
}

◆ 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 213 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 221 of file FitterBase.cpp.

                            {
// *******************
  if(FileSaved) return;
  //Stop Clock
  clock->Stop();
 
  outputFile->cd();
  outTree->Write();
 
  MACH3LOG_INFO("{} steps took {:.2f} seconds to complete. ({:.2f}s / step).", step - stepStart, clock->RealTime(), clock->RealTime() / static_cast<double>(step - stepStart));
  MACH3LOG_INFO("{} steps were accepted.", accCount);
  #ifdef 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;
}

◆ SaveSettings()

void FitterBase::SaveSettings ( )

protected

Save the settings that the MCMC was run with.

Definition at line 82 of file FitterBase.cpp.

                              {
// *******************
  if(SettingsSaved) return;
 
  outputFile->cd();
 
  TDirectory* MaCh3Version = outputFile->mkdir("MaCh3Engine");
  MaCh3Version->cd();
 
  if (std::getenv("MaCh3_ROOT") == NULL) {
    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") == NULL) {
    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();
 
  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", MaCh3Utils::getValue("VmRSS") / 1048576.0);
  MACH3LOG_WARN("\033[0;31mOut of Total available RAM {:.2f} GB\033[0m", MaCh3Utils::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, {} OscChannels",i , samples[i]->GetTitle(), samples[i]->GetNsamples(), samples[i]->GetNOscChannels());
 
  //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

Todo:: implement some check that number of params matches etc

Reimplemented in mcmc.

Definition at line 295 of file FitterBase.cpp.

                                                              {
// *******************
  MACH3LOG_INFO("Getting starting position from {}", FitName);
 
  TFile *infile = new TFile(FitName.c_str(), "READ");
  TTree *posts = infile->Get<TTree>("posteriors");
  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 and current one are different", FitName);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
 
      delete ConfigCov;
    }
 
    CovarianceFolder->Close();
    delete CovarianceFolder;
 
    std::vector<double> branch_vals(systematics[s]->GetNumParams(), M3::_BAD_DOUBLE_);
    for (int i = 0; i < systematics[s]->GetNumParams(); ++i) {
      posts->SetBranchAddress(systematics[s]->GetParName(i).c_str(), &branch_vals[i]);
    }
    posts->GetEntry(posts->GetEntries()-1);
 
    for (int i = 0; i < systematics[s]->GetNumParams(); ++i) {
      if(branch_vals[i] == M3::_BAD_DOUBLE_)
      {
        MACH3LOG_ERROR("Parameter {} is unvitalised with value {}", i, branch_vals[i]);
        MACH3LOG_ERROR("Please check more precisely chain you passed {}", FitName);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
    }
    systematics[s]->SetParameters(branch_vals);
    systematics[s]->AcceptStep();
 
    MACH3LOG_INFO("Printing new starting values for: {}", systematics[s]->GetName());
    systematics[s]->PrintNominalCurrProp();
 
    // Resetting branch adressed 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;
 
  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 103 of file FitterBase.h.

◆ accProb

double FitterBase::accProb

protected

current acceptance prob

Definition at line 101 of file FitterBase.h.

◆ auto_save

int FitterBase::auto_save

protected

auto save every N steps

Definition at line 139 of file FitterBase.h.

◆ clock

std::unique_ptr<TStopwatch> FitterBase::clock

protected

tells global time how long fit took

Definition at line 121 of file FitterBase.h.

◆ CovFolder

TDirectory* FitterBase::CovFolder

protected

Output cov folder.

Definition at line 133 of file FitterBase.h.

◆ FileSaved

bool FitterBase::FileSaved

protected

Checks if file saved not repeat some operations.

Definition at line 145 of file FitterBase.h.

◆ fitMan

manager* FitterBase::fitMan

protected

The manager.

Definition at line 92 of file FitterBase.h.

◆ fTestLikelihood

bool FitterBase::fTestLikelihood

protected

Necessary for some fitting algorithms like PSO.

Definition at line 142 of file FitterBase.h.

◆ logLCurr

double FitterBase::logLCurr

protected

current likelihood

Definition at line 97 of file FitterBase.h.

◆ logLProp

double FitterBase::logLProp

protected

proposed likelihood

Definition at line 99 of file FitterBase.h.

◆ outputFile

TFile* FitterBase::outputFile

protected

Output.

Definition at line 131 of file FitterBase.h.

◆ OutputPrepared

bool FitterBase::OutputPrepared

protected

Checks if output prepared not repeat some operations.

Definition at line 149 of file FitterBase.h.

◆ outTree

TTree* FitterBase::outTree

protected

Output tree with posteriors.

Definition at line 137 of file FitterBase.h.

◆ random

std::unique_ptr<TRandom3> FitterBase::random

protected

Random number.

Definition at line 128 of file FitterBase.h.

◆ sample_llh

std::vector<double> FitterBase::sample_llh

protected

store the llh breakdowns

Definition at line 108 of file FitterBase.h.

◆ SampleFolder

TDirectory* FitterBase::SampleFolder

protected

Output sample folder.

Definition at line 135 of file FitterBase.h.

◆ samples

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

protected

Sample holder.

Definition at line 113 of file FitterBase.h.

◆ SettingsSaved

bool FitterBase::SettingsSaved

protected

Checks if setting saved not repeat some operations.

Definition at line 147 of file FitterBase.h.

◆ step

unsigned int FitterBase::step

protected

current state

Definition at line 95 of file FitterBase.h.

◆ stepClock

std::unique_ptr<TStopwatch> FitterBase::stepClock

protected

tells how long single step/fit iteration took

Definition at line 123 of file FitterBase.h.

◆ stepStart

int FitterBase::stepStart

protected

step start if restarting

Definition at line 105 of file FitterBase.h.

◆ stepTime

double FitterBase::stepTime

protected

Time of single step.

Definition at line 125 of file FitterBase.h.

◆ syst_llh

std::vector<double> FitterBase::syst_llh

protected

systematic llh breakdowns

Definition at line 110 of file FitterBase.h.

◆ systematics

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

protected

Systematic holder.

Definition at line 118 of file FitterBase.h.

◆ TotalNSamples

unsigned int FitterBase::TotalNSamples

protected

Total number of samples used.

Definition at line 115 of file FitterBase.h.

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

Fitters/FitterBase.h
Fitters/FitterBase.cpp

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ FitterBase()

◆ ~FitterBase()

Member Function Documentation

◆ AddSampleHandler()

◆ AddSystObj()

◆ CheckSkipParameter()

◆ DragRace()

◆ GetName()

◆ GetScaneRange()

◆ GetStepScaleBasedOnLLHScan()

◆ PrepareOutput()

◆ ProcessMCMC()

◆ Run2DLLHScan()

◆ RunLLHScan()

◆ RunMCMC()

◆ RunSigmaVar()

◆ SanitiseInputs()

◆ SaveOutput()

◆ SaveSettings()

◆ StartFromPreviousFit()

Member Data Documentation

◆ accCount

◆ accProb

◆ auto_save

◆ clock

◆ CovFolder

◆ FileSaved

◆ fitMan

◆ fTestLikelihood

◆ logLCurr

◆ logLProp

◆ outputFile

◆ OutputPrepared

◆ outTree

◆ random

◆ sample_llh

◆ SampleFolder

◆ samples

◆ SettingsSaved

◆ step

◆ stepClock

◆ stepStart

◆ stepTime

◆ syst_llh

◆ systematics

◆ TotalNSamples