M3 Namespace Reference

Main namespace for MaCh3 software. More...

Namespaces
	detail
	Implementation details only for experts.
	Plotting
	Flexible, experiment-agnostic plotting utilities for MaCh3.
	Utils
	Utility helpers used across MaCh3.

Typedefs
using	float_t = double
using	int_t = int
using	uint_t = unsigned

Enumerations
enum class	BiasFunction { kGaussian , kVonMises , kGeneralisedGaussian }
enum	kInfCrit { kBIC , kDIC , kWAIC , kInfCrits }
	KS: Different Information Criterion tests mostly based Gelman paper. More...
enum	kReweightType {   kGaussian , kGaussian , kTGraph , kTGraph2D ,   kReweightTypes }
	Types of chain reweighting available. More...

Functions
template<typename T >
constexpr T	fmaf_t (T x, T y, T z)
	Function template for fused multiply-add. More...
int	GetNThreads ()
	number of threads which we need for example for TRandom3 More...
void	AddPath (std::string &FilePath)
	Prepends the MACH3 environment path to FilePath if it is not already present. More...
int	GetThreadIndex ()
	thread index inside parallel loop More...
TFile *	Open (const std::string &Name, const std::string &Type, const std::string &File, const int Line)
	Opens a ROOT file with the given name and mode. More...
void	ScaleHistogram (TH1 *Sample_Hist, const double scale)
	Scale histogram to get divided by bin width. More...
void	CheckBinningMatch (const TH1D *Hist1, const TH1D *Hist2, const std::string &File, const int Line)
	KS: Helper function check if data and MC binning matches. More...
void	CheckBinningMatch (const TH2D *Hist1, const TH2D *Hist2, const std::string &File, const int Line)
	KS: Helper function check if data and MC binning matches. More...
void	CheckBinningMatch (TH2Poly *Hist1, TH2Poly *Hist2, const std::string &File, const int Line)
	KS: Helper function check if data and MC binning matches. More...
YAML::Node	PolyToYaml (TH2Poly *Hist, const std::string &YamlName, const std::string &File, const int Line)
	KS: Convert TH2Poly into yaml config accepted by MaCh3. More...
template<typename ObjectType >
std::unique_ptr< ObjectType >	Clone (const ObjectType *obj, const std::string &name="")
	KS: Creates a copy of a ROOT-like object and wraps it in a smart pointer. More...
bool	CaseInsentiveMatch (std::string Text, std::string Pattern)
	Matches a string against a simple wildcard Pattern using regex. Is not case sensitive. More...
bool	CaseInsensitiveMatchAny (std::string Text, const std::vector< std::string > &Patterns)
	Matches a string against a simple wildcard Pattern using regex. Is not case sensitive. More...
double *	MatrixMult (double *A, double *B, int n)
	CW: Multi-threaded matrix multiplication. More...
double **	MatrixMult (double **A, double **B, int n)
	CW: Multi-threaded matrix multiplication. More...
TMatrixD	MatrixMult (TMatrixD A, TMatrixD B)
	CW: Multi-threaded matrix multiplication. More...
void	MatrixVectorMulti (double *_restrict_ VecMulti, double **_restrict_ matrix, const double *_restrict_ vector, const int n)
	KS: Custom function to perform multiplication of matrix and vector with multithreading. More...
double	MatrixVectorMultiSingle (double **_restrict_ matrix, const double *_restrict_ vector, const int Length, const int i)
	KS: Custom function to perform multiplication of matrix and single element which is thread safe. More...
void	FixSampleNamesQuotes (std::string &yamlStr)
	KS: Yaml emitter has problem and drops "", if you have special signs in you like * then there is problem. This bit hacky code adds these "". More...
void	AddTuneValues (YAML::Node &root, const std::vector< double > &Values, const std::string &Tune, const std::vector< std::string > &FancyNames={})
	KS: Add Tune values to YAML covariance matrix. More...
void	MakeCorrelationMatrix (YAML::Node &root, const std::vector< double > &Values, const std::vector< double > &Errors, const std::vector< std::vector< double >> &Correlation, const std::string &OutYAMLName, const std::vector< std::string > &FancyNames={})
	KS: Replace correlation matrix and tune values in YAML covariance matrix. More...
TMacro *	GetConfigMacroFromChain (TDirectory *CovarianceFolder)
	KS: We store configuration macros inside the chain. In the past, multiple configs were stored, which required error-prone hardcoding like "Config_xsec_cov". Therefore, this code maintains backward compatibility by checking the number of macros present and using a hardcoded name only if necessary. More...
TMatrixDSym *	GetCovMatrixFromChain (TDirectory *TempFile)
	KS: Retrieve the cross-section covariance matrix from the given TDirectory. Historically, multiple covariance matrices could be stored, requiring fragile hardcoded paths like "CovarianceFolder/xsec_cov". This function maintains backward compatibility by: More...
std::vector< std::vector< double > >	GetCholeskyDecomposedMatrix (const TMatrixDSym &matrix, const std::string &matrixName)
	Computes Cholesky decomposition of a symmetric positive definite matrix using custom function which can be even 20 times faster. More...
bool	CanDecomposeMatrix (const TMatrixDSym &matrix)
	Checks if a matrix can be Cholesky decomposed. More...
int	MakeMatrixPosDef (TMatrixDSym *cov)
	Makes sure that matrix is positive-definite by adding a small number to on-diagonal elements. More...
void	DumpParamHandlerToFile (const int _fNumPar, const std::vector< double > &_fPreFitValue, const std::vector< double > &_fError, const std::vector< double > &_fLowBound, const std::vector< double > &_fUpBound, const std::vector< double > &_fIndivStepScale, const std::vector< std::string > &_fFancyNames, const std::vector< bool > &_fFlatPrior, const std::vector< SplineParameter > &SplineParams, TMatrixDSym *covMatrix, TH2D *CorrMatrix, const std::string &Name)
	Dump Matrix to ROOT file, useful when we need to pass matrix info to another fitting group. More...
void	DebugPCA (const double sum, const TMatrixD &temp, const TMatrixD &eigen_vectors, const TVectorD &eigen_values, const TMatrixD &TransferMat, const TMatrixD &TransferMatT, const int NumPar, const int FirstPCAdpar, const int LastPCAdpar, const int nKeptPCApars, const double eigen_threshold)
	KS: Let's dump all useful matrices to properly validate PCA. More...

Variables
constexpr static const char *	float_t_str_repr = "D"
constexpr static const double	_BAD_DOUBLE_ = -999.99
	Default value used for double initialisation. More...
constexpr static const int	_BAD_INT_ = -999
	Default value used for int initialisation. More...
constexpr static const double	_DEFAULT_RETURN_VAL_ = -999999.123456
constexpr static const double	Unity_D = 1.
	Some commonly used variables to which we set pointers to. More...
constexpr static const float	Unity_F = 1.
constexpr static const float_t	Unity = Unity_D
constexpr static const double	Zero_D = 0.
constexpr static const float	Zero_F = 0.
constexpr static const float_t	Zero = Zero_D
constexpr static const double	KinematicLowBound = std::numeric_limits<double>::lowest()
	When parameter has no bound this serves as it. Lowest possible value the system. More...
constexpr static const double	KinematicUpBound = std::numeric_limits<double>::max()
	When parameter has no bound this serves as it. Highest possible value the system. More...
constexpr static const double	_LARGE_LOGL_ = 1234567890.0
	Large Likelihood is used it parameter go out of physical boundary, this indicates in MCMC that such step should be removed. More...
constexpr static const double	_LOW_MC_BOUND_ = .00001
	MC prediction lower bound in bin to identify problematic binning definitions and handle LogL calculation. More...
constexpr static const double	DefSplineKnotUpBound = 9999
	Default value for spline knot capping, default mean not capping is being applied. More...
constexpr static const double	DefSplineKnotLowBound = -9999
	Default value for spline knot capping, default mean not capping is being applied. More...
constexpr static const int	UnderOverFlowBin = -1
	Mark bin which is overflow or underflow in MaCh3 binning. More...

Detailed Description

Main namespace for MaCh3 software.

Typedef Documentation

float_t

typedef double M3::float_t

Definition at line 37 of file Core.h.

int_t

using M3::int_t = typedef int

Definition at line 38 of file Core.h.

uint_t

using M3::uint_t = typedef unsigned

Definition at line 39 of file Core.h.

Enumeration Type Documentation

BiasFunction

enum M3::BiasFunction

strong

Enumerator
kGaussian
kVonMises
kGeneralisedGaussian

Definition at line 15 of file MulticanonicalMCMCHandler.h.

                          {
    kGaussian,
    kVonMises,
    kGeneralisedGaussian
  };

kInfCrit

enum M3::kInfCrit

KS: Different Information Criterion tests mostly based Gelman paper.

Enumerator
kBIC	Bayesian Information Criterion.
kDIC	Deviance Information Criterion.
kWAIC	Watanabe-Akaike information criterion.
kInfCrits	This only enumerates.

Definition at line 32 of file StatisticalUtils.h.

                {
    kBIC,      
    kDIC,      
    kWAIC,     
    kInfCrits  
  };

kReweightType

enum M3::kReweightType

Types of chain reweighting available.

Enumerator
kGaussian	Assumes gaussian prior.
kGaussian
kTGraph	Calculates Likelihood based on TGraph.
kTGraph2D	Calculates Likelihood based on TGraph2D.
kReweightTypes	This only enumerates.

Definition at line 33 of file ReweightMCMC.cpp.

                     {
    kGaussian,     
    kTGraph,       
    kTGraph2D,     
    kReweightTypes 
  };

Function Documentation

AddPath()

void M3::AddPath

(

std::string &

FilePath

)

Prepends the MACH3 environment path to FilePath if it is not already present.

Parameters

FilePath

Reference to the file path string to be modified.

Definition at line 382 of file Monitor.cpp.

                                  {
// ***************************************************************************
  //KS:Most inputs are in ${MACH3}/inputs/blarb.root
  if (std::getenv("MACH3") == nullptr) {
    MACH3LOG_ERROR("MACH3 is not defined");
    MACH3LOG_ERROR("Please source your setup script");
    MACH3LOG_ERROR("Read more: https://mach3-software.github.io/MaCh3/FAQ.html");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  std::string MaCh3Path = std::string(std::getenv("MACH3")) + "/";
   // Check if FilePath does NOT start with MaCh3Path
  if (FilePath.find(MaCh3Path) != 0) {
    FilePath.insert(0, MaCh3Path);
  }
}

AddTuneValues()

void M3::AddTuneValues ( YAML::Node &  root, const std::vector< double > &  Values, const std::string &  Tune, const std::vector< std::string > &  FancyNames = {}  )

inline

KS: Add Tune values to YAML covariance matrix.

Parameters

root	The root YAML node to be updated.
Values	The values to add for the specified tune.
Tune	The name of the tune (e.g., "PostFit").
FancyNames	Optional list of fancy names to match systematics (must match Values size if provided).

Definition at line 233 of file ParameterHandlerUtils.h.

                                                                   {}) {
// *************************************
  YAML::Node NodeCopy = YAML::Clone(root);
  YAML::Node systematics = NodeCopy["Systematics"];
 
  if (!systematics || !systematics.IsSequence()) {
    MACH3LOG_ERROR("'Systematics' node is missing or not a sequence in the YAML copy");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  if (!FancyNames.empty() && FancyNames.size() != Values.size()) {
    MACH3LOG_ERROR("Mismatch in sizes: FancyNames has {}, but Values has {}", FancyNames.size(), Values.size());
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  if (FancyNames.empty() && systematics.size() != Values.size()) {
    MACH3LOG_ERROR("Mismatch in sizes: Values has {}, but YAML 'Systematics' has {} entries",
                   Values.size(), systematics.size());
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  if (!FancyNames.empty()) {
    for (std::size_t i = 0; i < FancyNames.size(); ++i) {
      bool matched = false;
      for (std::size_t j = 0; j < systematics.size(); ++j) {
        YAML::Node systematicNode = systematics[j]["Systematic"];
        if (!systematicNode) continue;
        auto nameNode = systematicNode["Names"];
        if (!nameNode || !nameNode["FancyName"]) continue;
        if (nameNode["FancyName"].as<std::string>() == FancyNames[i]) {
          if (!systematicNode["ParameterValues"]) {
            MACH3LOG_ERROR("Missing 'ParameterValues' for matched FancyName '{}'", FancyNames[i]);
            throw MaCh3Exception(__FILE__, __LINE__);
          }
          systematicNode["ParameterValues"][Tune] = M3::Utils::FormatDouble(Values[i], 4);
          matched = true;
          break;
        }
      }
      if (!matched) {
        MACH3LOG_ERROR("Could not find a matching FancyName '{}' in the systematics", FancyNames[i]);
        throw MaCh3Exception(__FILE__, __LINE__);
      }
    }
  } else {
    for (std::size_t i = 0; i < systematics.size(); ++i) {
      YAML::Node systematicNode = systematics[i]["Systematic"];
      if (!systematicNode || !systematicNode["ParameterValues"]) {
        MACH3LOG_ERROR("Missing 'Systematic' or 'ParameterValues' entry at index {}", i);
        throw MaCh3Exception(__FILE__, __LINE__);
      }
      systematicNode["ParameterValues"][Tune] = M3::Utils::FormatDouble(Values[i], 4);
    }
  }
 
  // Convert updated copy to string
  std::string YAMLString = YAMLtoSTRING(NodeCopy);
  FixSampleNamesQuotes(YAMLString);
  // Write to output file
  std::string OutName = "UpdatedMatrixWithTune" + Tune + ".yaml";
  std::ofstream outFile(OutName);
  if (!outFile) {
    MACH3LOG_ERROR("Failed to open file for writing: {}", OutName);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  outFile << YAMLString;
  outFile.close();
}

CanDecomposeMatrix()

bool M3::CanDecomposeMatrix ( const TMatrixDSym &  matrix )

inline

Checks if a matrix can be Cholesky decomposed.

Parameters

matrix

Input symmetric matrix to test

Definition at line 544 of file ParameterHandlerUtils.h.

                                                          {
// *************************************
  TDecompChol chdcmp(matrix);
  return chdcmp.Decompose();
}

CaseInsensitiveMatchAny()

bool M3::CaseInsensitiveMatchAny ( std::string  Text, const std::vector< std::string > &  Patterns  )

inline

Matches a string against a simple wildcard Pattern using regex. Is not case sensitive.

Parameters

Text	Input string to test.
Patterns	Collection wildcard patterns to match against.

Definition at line 66 of file ParameterHandlerUtils.h.

                                                                                          {
  for (size_t i = 0; i < Patterns.size(); i++) {
    if (M3::CaseInsentiveMatch(Text, Patterns[i])) {
      return true;
    }
  }
  return false;
}

CaseInsentiveMatch()

bool M3::CaseInsentiveMatch ( std::string  Text, std::string  Pattern  )

inline

Matches a string against a simple wildcard Pattern using regex. Is not case sensitive.

Parameters

Text	Input string to test.
Pattern	Wildcard pattern to match against.

Definition at line 45 of file ParameterHandlerUtils.h.

                                                                  {
  // Make a copy and to lower case to not be case sensitive
  std::transform(Text.begin(), Text.end(), Text.begin(), ::tolower);
 
  // Convert to low case to not be case sensitive
  std::transform(Pattern.begin(), Pattern.end(), Pattern.begin(), ::tolower);
  try {
    // Replace '*' in the Pattern with '.*' for regex matching
    std::string RegexPattern = "^" + std::regex_replace(Pattern, std::regex("\\*"), ".*") + "$";
    std::regex Regex(RegexPattern);
    return std::regex_match(Text, Regex);
  }
  catch (const std::regex_error& e) {
    MACH3LOG_ERROR("Regex error: {}", e.what());
    return false;
  }
}

CheckBinningMatch() [1/3]

void M3::CheckBinningMatch	(	const TH1D *	Hist1,
		const TH1D *	Hist2,
		const std::string &	File,
		const int	Line
	)

KS: Helper function check if data and MC binning matches.

Definition at line 821 of file HistogramUtils.cpp.

                                                                                                    {
// ***************************************************************************
  if (Hist1->GetNbinsX() != Hist2->GetNbinsX()) {
    MACH3LOG_ERROR("Number of bins does not match for TH1D: {} vs {}", Hist1->GetNbinsX(), Hist2->GetNbinsX());
    throw MaCh3Exception(File, Line);
  }
  for (int i = 1; i <= Hist1->GetNbinsX(); ++i) {
    if (std::fabs(Hist1->GetXaxis()->GetBinLowEdge(i) - Hist2->GetXaxis()->GetBinLowEdge(i)) > 0.001 ||
      std::fabs(Hist1->GetXaxis()->GetBinUpEdge(i) - Hist2->GetXaxis()->GetBinUpEdge(i)) > 0.001) {
      MACH3LOG_ERROR("Bin edges do not match for TH1D at bin {}", i);
      throw MaCh3Exception(File, Line);
    }
  }
}

CheckBinningMatch() [2/3]

void M3::CheckBinningMatch	(	const TH2D *	Hist1,
		const TH2D *	Hist2,
		const std::string &	File,
		const int	Line
	)

KS: Helper function check if data and MC binning matches.

Definition at line 837 of file HistogramUtils.cpp.

                                                                                                    {
// ***************************************************************************
  if (Hist1->GetNbinsX() != Hist2->GetNbinsX() || Hist1->GetNbinsY() != Hist2->GetNbinsY()) {
    MACH3LOG_ERROR("Number of bins does not match for TH2D");
    throw MaCh3Exception(File, Line);
  }
 
  for (int i = 1; i <= Hist1->GetNbinsX(); ++i) {
    if (std::fabs(Hist1->GetXaxis()->GetBinLowEdge(i) - Hist2->GetXaxis()->GetBinLowEdge(i)) > 0.001 ||
      std::fabs(Hist1->GetXaxis()->GetBinUpEdge(i) - Hist2->GetXaxis()->GetBinUpEdge(i)) > 0.001) {
      MACH3LOG_ERROR("X bin edges do not match for TH2D at bin {}", i);
      throw MaCh3Exception(File, Line);
    }
  }
  for (int j = 1; j <= Hist1->GetNbinsY(); ++j) {
    if (std::fabs(Hist1->GetYaxis()->GetBinLowEdge(j) - Hist2->GetYaxis()->GetBinLowEdge(j)) > 0.001 ||
      std::fabs(Hist1->GetYaxis()->GetBinUpEdge(j) - Hist2->GetYaxis()->GetBinUpEdge(j)) > 0.001) {
      MACH3LOG_ERROR("Y bin edges do not match for TH2D at bin {}", j);
      throw MaCh3Exception(File, Line);
    }
  }
}

CheckBinningMatch() [3/3]

void M3::CheckBinningMatch	(	TH2Poly *	Hist1,
		TH2Poly *	Hist2,
		const std::string &	File,
		const int	Line
	)

KS: Helper function check if data and MC binning matches.

Definition at line 862 of file HistogramUtils.cpp.

                                                                                              {
// ***************************************************************************
  int NBins1 = Hist1->GetNumberOfBins();
  int NBins2 = Hist2->GetNumberOfBins();
  if (NBins1 != NBins2) {
    MACH3LOG_ERROR("Number of bins does not match for TH2Poly: {} vs {}", NBins1, NBins2);
    throw MaCh3Exception(File, Line);
  }
  for (int j = 1; j <= NBins1; j++)
  {
    //KS: There is weird offset between bin content and GetBins so this is correct, in spite of looking funny
    TH2PolyBin* polybin1 = static_cast<TH2PolyBin*>(Hist1->GetBins()->At(j - 1));
    TH2PolyBin* polybin2 = static_cast<TH2PolyBin*>(Hist2->GetBins()->At(j - 1));
 
    if( std::fabs(polybin2->GetXMin() - polybin1->GetXMin()) > 0.001 ||
      std::fabs(polybin2->GetXMax() - polybin1->GetXMax()) > 0.001 ||
      std::fabs(polybin2->GetYMin() - polybin1->GetYMin()) > 0.001 ||
      std::fabs(polybin2->GetYMax() - polybin1->GetYMax()) > 0.001  )
    {
      MACH3LOG_ERROR("Binning doesn't match", Hist1->GetTitle());
      MACH3LOG_ERROR("data  x min {} x max {} y min {} y max {}", polybin2->GetXMin(), polybin2->GetXMax(), polybin2->GetYMin(), polybin2->GetYMax());
      MACH3LOG_ERROR("mc    x min {} x max {} y min {} y max {}", polybin1->GetXMin(), polybin1->GetXMax(), polybin1->GetYMin(), polybin1->GetYMax());
      throw MaCh3Exception(File , Line );
    }
  }
}

Clone()

template<typename ObjectType >

std::unique_ptr<ObjectType> M3::Clone	(	const ObjectType *	obj,
		const std::string &	name = ""
	)

KS: Creates a copy of a ROOT-like object and wraps it in a smart pointer.

Template Parameters

ObjectType

The type of the object to clone for example TH1D or TH2Poly.

Parameters

obj	Pointer to the object to clone.
name	Optional argument allowing to set new name of cloned object

Returns

std::unique_ptr<ObjectType> Owning pointer to the cloned object.

Definition at line 155 of file HistogramUtils.h.

                                                                                   {
  std::string cloneName = name.empty() ? obj->GetName() : name;
 
  std::unique_ptr<ObjectType> Hist(static_cast<ObjectType*>(obj->Clone(cloneName.c_str())));
  // Disable ROOT memory management because it causes lot of headache especially as smart pointers are much smarter
  Hist->SetDirectory(nullptr);
 
  return Hist;
}

DebugPCA()

void M3::DebugPCA ( const double  sum, const TMatrixD &  temp, const TMatrixD &  eigen_vectors, const TVectorD &  eigen_values, const TMatrixD &  TransferMat, const TMatrixD &  TransferMatT, const int  NumPar, const int  FirstPCAdpar, const int  LastPCAdpar, const int  nKeptPCApars, const double  eigen_threshold  )

inline

KS: Let's dump all useful matrices to properly validate PCA.

Definition at line 710 of file ParameterHandlerUtils.h.

                                                   {
// ********************************************
  #pragma GCC diagnostic push
  #pragma GCC diagnostic ignored "-Wfloat-conversion"
  int originalErrorWarning = gErrorIgnoreLevel;
  gErrorIgnoreLevel = kFatal;
 
  TDirectory *ogdir = gDirectory;
 
  TFile *PCA_Debug = new TFile("Debug_PCA.root", "RECREATE");
  PCA_Debug->cd();
 
  bool PlotText = true;
  //KS: If we have more than 200 plot becomes unreadable :(
  if(NumPar > 200) PlotText = false;
 
  auto heigen_values     = std::make_unique<TH1D>("eigen_values", "Eigen Values", eigen_values.GetNrows(), 0.0, eigen_values.GetNrows());
  heigen_values->SetDirectory(nullptr);
  auto heigen_cumulative = std::make_unique<TH1D>("heigen_cumulative", "heigen_cumulative", eigen_values.GetNrows(), 0.0, eigen_values.GetNrows());
  heigen_cumulative->SetDirectory(nullptr);
  auto heigen_frac       = std::make_unique<TH1D>("heigen_fractional", "heigen_fractional", eigen_values.GetNrows(), 0.0, eigen_values.GetNrows());
  heigen_frac->SetDirectory(nullptr);
  heigen_values->GetXaxis()->SetTitle("Eigen Vector");
  heigen_values->GetYaxis()->SetTitle("Eigen Value");
 
  double Cumulative = 0;
  for(int i = 0; i < eigen_values.GetNrows(); i++)
  {
    heigen_values->SetBinContent(i+1, (eigen_values)(i));
    heigen_cumulative->SetBinContent(i+1, (eigen_values)(i)/sum + Cumulative);
    heigen_frac->SetBinContent(i+1, (eigen_values)(i)/sum);
    Cumulative += (eigen_values)(i)/sum;
  }
  heigen_values->Write("heigen_values");
  eigen_values.Write("eigen_values");
  heigen_cumulative->Write("heigen_values_cumulative");
  heigen_frac->Write("heigen_values_frac");
 
  TH2D* heigen_vectors = new TH2D(eigen_vectors);
  heigen_vectors->GetXaxis()->SetTitle("Parameter in Normal Base");
  heigen_vectors->GetYaxis()->SetTitle("Parameter in Decomposed Base");
  heigen_vectors->Write("heigen_vectors");
  eigen_vectors.Write("eigen_vectors");
 
  TH2D* SubsetPCA = new TH2D(temp);
  SubsetPCA->GetXaxis()->SetTitle("Parameter in Normal Base");
  SubsetPCA->GetYaxis()->SetTitle("Parameter in Decomposed Base");
 
  SubsetPCA->Write("hSubsetPCA");
  temp.Write("SubsetPCA");
  TH2D* hTransferMat = new TH2D(TransferMat);
  hTransferMat->GetXaxis()->SetTitle("Parameter in Normal Base");
  hTransferMat->GetYaxis()->SetTitle("Parameter in Decomposed Base");
  TH2D* hTransferMatT = new TH2D(TransferMatT);
 
  hTransferMatT->GetXaxis()->SetTitle("Parameter in Decomposed Base");
  hTransferMatT->GetYaxis()->SetTitle("Parameter in Normal Base");
 
  hTransferMat->Write("hTransferMat");
  TransferMat.Write("TransferMat");
  hTransferMatT->Write("hTransferMatT");
  TransferMatT.Write("TransferMatT");
 
  auto c1 = std::make_unique<TCanvas>("c1", " ", 0, 0, 1024, 1024);
  c1->SetBottomMargin(0.1);
  c1->SetTopMargin(0.05);
  c1->SetRightMargin(0.05);
  c1->SetLeftMargin(0.12);
  c1->SetGrid();
 
  gStyle->SetPaintTextFormat("4.1f");
  gStyle->SetOptFit(0);
  gStyle->SetOptStat(0);
  // Make pretty correlation colors (red to blue)
  constexpr int NRGBs = 5;
  TColor::InitializeColors();
  Double_t stops[NRGBs] = { 0.00, 0.25, 0.50, 0.75, 1.00 };
  Double_t red[NRGBs]   = { 0.00, 0.25, 1.00, 1.00, 0.50 };
  Double_t green[NRGBs] = { 0.00, 0.25, 1.00, 0.25, 0.00 };
  Double_t blue[NRGBs]  = { 0.50, 1.00, 1.00, 0.25, 0.00 };
  TColor::CreateGradientColorTable(5, stops, red, green, blue, 255);
  gStyle->SetNumberContours(255);
 
  double maxz = 0;
  double minz = 0;
 
  c1->Print("Debug_PCA.pdf[");
  auto EigenLine = std::make_unique<TLine>(nKeptPCApars, 0, nKeptPCApars, heigen_cumulative->GetMaximum());
  EigenLine->SetLineColor(kPink);
  EigenLine->SetLineWidth(2);
  EigenLine->SetLineStyle(kSolid);
 
  auto text = std::make_unique<TText>(0.5, 0.5, Form("Threshold = %g", eigen_threshold));
  text->SetTextFont (43);
  text->SetTextSize (40);
 
  heigen_values->SetLineColor(kRed);
  heigen_values->SetLineWidth(2);
  heigen_cumulative->SetLineColor(kGreen);
  heigen_cumulative->SetLineWidth(2);
  heigen_frac->SetLineColor(kBlue);
  heigen_frac->SetLineWidth(2);
 
  c1->SetLogy();
  heigen_values->SetMaximum(heigen_cumulative->GetMaximum()+heigen_cumulative->GetMaximum()*0.4);
  heigen_values->Draw();
  heigen_frac->Draw("SAME");
  heigen_cumulative->Draw("SAME");
  EigenLine->Draw("Same");
  text->DrawTextNDC(0.42, 0.84,Form("Threshold = %g", eigen_threshold));
 
  auto leg = std::make_unique<TLegend>(0.2, 0.2, 0.6, 0.5);
  leg->SetTextSize(0.04);
  leg->AddEntry(heigen_values.get(), "Absolute", "l");
  leg->AddEntry(heigen_frac.get(), "Fractional", "l");
  leg->AddEntry(heigen_cumulative.get(), "Cumulative", "l");
 
  leg->SetLineColor(0);
  leg->SetLineStyle(0);
  leg->SetFillColor(0);
  leg->SetFillStyle(0);
  leg->Draw("Same");
 
  c1->Print("Debug_PCA.pdf");
  c1->SetRightMargin(0.15);
  c1->SetLogy(0);
 
  heigen_vectors->SetMarkerSize(0.2);
  minz = heigen_vectors->GetMinimum();
  if (fabs(0-maxz)>fabs(0-minz)) heigen_vectors->GetZaxis()->SetRangeUser(0-fabs(0-maxz),0+fabs(0-maxz));
  else heigen_vectors->GetZaxis()->SetRangeUser(0-fabs(0-minz),0+fabs(0-minz));
  if(PlotText) heigen_vectors->Draw("COLZ TEXT");
  else heigen_vectors->Draw("COLZ");
 
  auto Eigen_Line = std::make_unique<TLine>(0, nKeptPCApars, LastPCAdpar - FirstPCAdpar, nKeptPCApars);
  Eigen_Line->SetLineColor(kGreen);
  Eigen_Line->SetLineWidth(2);
  Eigen_Line->SetLineStyle(kDotted);
  Eigen_Line->Draw("SAME");
  c1->Print("Debug_PCA.pdf");
 
  SubsetPCA->SetMarkerSize(0.2);
  minz = SubsetPCA->GetMinimum();
  if (fabs(0-maxz)>fabs(0-minz)) SubsetPCA->GetZaxis()->SetRangeUser(0-fabs(0-maxz),0+fabs(0-maxz));
  else SubsetPCA->GetZaxis()->SetRangeUser(0-fabs(0-minz),0+fabs(0-minz));
  if(PlotText) SubsetPCA->Draw("COLZ TEXT");
  else SubsetPCA->Draw("COLZ");
  c1->Print("Debug_PCA.pdf");
  delete SubsetPCA;
 
  hTransferMat->SetMarkerSize(0.15);
  minz = hTransferMat->GetMinimum();
  if (fabs(0-maxz)>fabs(0-minz)) hTransferMat->GetZaxis()->SetRangeUser(0-fabs(0-maxz),0+fabs(0-maxz));
  else hTransferMat->GetZaxis()->SetRangeUser(0-fabs(0-minz),0+fabs(0-minz));
  if(PlotText) hTransferMat->Draw("COLZ TEXT");
  else hTransferMat->Draw("COLZ");
  c1->Print("Debug_PCA.pdf");
  delete hTransferMat;
 
  hTransferMatT->SetMarkerSize(0.15);
  minz = hTransferMatT->GetMinimum();
  if (fabs(0-maxz)>fabs(0-minz)) hTransferMatT->GetZaxis()->SetRangeUser(0-fabs(0-maxz),0+fabs(0-maxz));
  else hTransferMatT->GetZaxis()->SetRangeUser(0-fabs(0-minz),0+fabs(0-minz));
  if(PlotText) hTransferMatT->Draw("COLZ TEXT");
  else hTransferMatT->Draw("COLZ");
  c1->Print( "Debug_PCA.pdf");
  delete hTransferMatT;
 
  delete heigen_vectors;
 
  c1->Print("Debug_PCA.pdf]");
  PCA_Debug->Close();
  delete PCA_Debug;
  gErrorIgnoreLevel = originalErrorWarning;
  ogdir->cd();  // go back to original directory
  #pragma GCC diagnostic pop
}

DumpParamHandlerToFile()

void M3::DumpParamHandlerToFile ( const int  _fNumPar, const std::vector< double > &  _fPreFitValue, const std::vector< double > &  _fError, const std::vector< double > &  _fLowBound, const std::vector< double > &  _fUpBound, const std::vector< double > &  _fIndivStepScale, const std::vector< std::string > &  _fFancyNames, const std::vector< bool > &  _fFlatPrior, const std::vector< SplineParameter > &  SplineParams, TMatrixDSym *  covMatrix, TH2D *  CorrMatrix, const std::string &  Name  )

inline

Dump Matrix to ROOT file, useful when we need to pass matrix info to another fitting group.

Warning

This is mostly used for backward compatibility

Definition at line 610 of file ParameterHandlerUtils.h.

                                                          {
// ********************************************
  TFile* outputFile = new TFile(Name.c_str(), "RECREATE");
 
  TObjArray* param_names = new TObjArray();
  TObjArray* spline_interpolation = new TObjArray();
  TObjArray* spline_names = new TObjArray();
 
  TVectorD* param_prior = new TVectorD(_fNumPar);
  TVectorD* flat_prior = new TVectorD(_fNumPar);
  TVectorD* stepscale = new TVectorD(_fNumPar);
  TVectorD* param_lb = new TVectorD(_fNumPar);
  TVectorD* param_ub = new TVectorD(_fNumPar);
 
  TVectorD* param_knot_weight_lb = new TVectorD(_fNumPar);
  TVectorD* param_knot_weight_ub = new TVectorD(_fNumPar);
  TVectorD* error = new TVectorD(_fNumPar);
 
  for(int i = 0; i < _fNumPar; ++i)
  {
    TObjString* nameObj = new TObjString(_fFancyNames[i].c_str());
    param_names->AddLast(nameObj);
 
    TObjString* splineType = new TObjString("TSpline3");
    spline_interpolation->AddLast(splineType);
 
    TObjString* splineName = new TObjString("");
    spline_names->AddLast(splineName);
 
    (*param_prior)[i] = _fPreFitValue[i];
    (*flat_prior)[i] = _fFlatPrior[i];
    (*stepscale)[i] = _fIndivStepScale[i];
    (*error)[i] = _fError[i];
 
    (*param_lb)[i] = _fLowBound[i];
    (*param_ub)[i] = _fUpBound[i];
 
    //Default values
    (*param_knot_weight_lb)[i] = -9999;
    (*param_knot_weight_ub)[i] = +9999;
  }
 
  for (size_t SplineIndex = 0; SplineIndex < SplineParams.size(); SplineIndex++ ) {
    auto SystIndex = SplineParams[SplineIndex].index;
 
    (*param_knot_weight_lb)[SystIndex] = SplineParams.at(SplineIndex)._SplineKnotLowBound;
    (*param_knot_weight_ub)[SystIndex] = SplineParams.at(SplineIndex)._SplineKnotUpBound;
 
    TObjString* splineType = new TObjString(SplineInterpolation_ToString(SplineParams.at(SplineIndex)._SplineInterpolationType).c_str());
    spline_interpolation->AddAt(splineType, SystIndex);
 
    TObjString* splineName = new TObjString(SplineParams[SplineIndex]._fSplineNames.c_str());
    spline_names->AddAt(splineName, SystIndex);
  }
  param_names->Write("xsec_param_names", TObject::kSingleKey);
  delete param_names;
  spline_interpolation->Write("xsec_spline_interpolation", TObject::kSingleKey);
  delete spline_interpolation;
  spline_names->Write("xsec_spline_names", TObject::kSingleKey);
  delete spline_names;
 
  param_prior->Write("xsec_param_prior");
  delete param_prior;
  flat_prior->Write("xsec_flat_prior");
  delete flat_prior;
  stepscale->Write("xsec_stepscale");
  delete stepscale;
  param_lb->Write("xsec_param_lb");
  delete param_lb;
  param_ub->Write("xsec_param_ub");
  delete param_ub;
 
  param_knot_weight_lb->Write("xsec_param_knot_weight_lb");
  delete param_knot_weight_lb;
  param_knot_weight_ub->Write("xsec_param_knot_weight_ub");
  delete param_knot_weight_ub;
  error->Write("xsec_error");
  delete error;
 
  covMatrix->Write("xsec_cov");
  CorrMatrix->Write("hcov");
 
  outputFile->Close();
  delete outputFile;
}

FixSampleNamesQuotes()

void M3::FixSampleNamesQuotes ( std::string &  yamlStr )

inline

KS: Yaml emitter has problem and drops "", if you have special signs in you like * then there is problem. This bit hacky code adds these "".

Parameters

yamlStr

The YAML string to be processed (modified in-place).

Definition at line 198 of file ParameterHandlerUtils.h.

                                                     {
// *************************************
  std::stringstream input(yamlStr);
  std::string line;
  std::string fixedYaml;
  std::regex sampleNamesRegex(R"(SampleNames:\s*\[([^\]]+)\])");
 
  while (std::getline(input, line)) {
    std::smatch match;
    if (std::regex_search(line, match, sampleNamesRegex)) {
      std::string contents = match[1];  // inside the brackets
      std::stringstream ss(contents);
      std::string item;
      std::vector<std::string> quotedItems;
 
      while (std::getline(ss, item, ',')) {
        item = std::regex_replace(item, std::regex(R"(^\s+|\s+$)"), ""); // trim
        quotedItems.push_back("\"" + item + "\"");
      }
 
      std::string replacement = "SampleNames: [" + fmt::format("{}", fmt::join(quotedItems, ", ")) + "]";
      line = std::regex_replace(line, sampleNamesRegex, replacement);
    }
    fixedYaml += line + "\n";
  }
 
  yamlStr = fixedYaml;
}

fmaf_t()

template<typename T >

constexpr T M3::fmaf_t ( T  x, T  y, T  z  )

constexpr

Function template for fused multiply-add.

Definition at line 45 of file Core.h.

                                    {
    #ifdef _LOW_MEMORY_STRUCTS_
    return std::fmaf(x, y, z);
    #else
    return std::fma(x, y, z);
    #endif
  }

GetCholeskyDecomposedMatrix()

std::vector<std::vector<double> > M3::GetCholeskyDecomposedMatrix ( const TMatrixDSym &  matrix, const std::string &  matrixName  )

inline

Computes Cholesky decomposition of a symmetric positive definite matrix using custom function which can be even 20 times faster.

Parameters

matrix	Input symmetric positive definite matrix
matrixName	Identifier for error reporting

Definition at line 508 of file ParameterHandlerUtils.h.

                                                                                                                        {
// *************************************
  const Int_t n = matrix.GetNrows();
  std::vector<std::vector<double>> L(n, std::vector<double>(n, 0.0));
 
  for (Int_t j = 0; j < n; ++j) {
    // Compute diagonal element (must be serial)
    double sum_diag = matrix(j, j);
    for (Int_t k = 0; k < j; ++k) {
      sum_diag -= L[j][k] * L[j][k];
    }
    const double tol = 1e-15;
    if (sum_diag <= tol) {
      MACH3LOG_ERROR("Cholesky decomposition failed for {} (non-positive diagonal)", matrixName);
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    L[j][j] = std::sqrt(sum_diag);
 
    // Compute the rest of the column in parallel
    #ifdef MULTITHREAD
    #pragma omp parallel for
    #endif
    for (Int_t i = j + 1; i < n; ++i) {
      double sum = matrix(i, j);
      for (Int_t k = 0; k < j; ++k) {
        sum -= L[i][k] * L[j][k];
      }
      L[i][j] = sum / L[j][j];
    }
  }
  return L;
}

GetConfigMacroFromChain()

TMacro* M3::GetConfigMacroFromChain ( TDirectory *  CovarianceFolder )

inline

KS: We store configuration macros inside the chain. In the past, multiple configs were stored, which required error-prone hardcoding like "Config_xsec_cov". Therefore, this code maintains backward compatibility by checking the number of macros present and using a hardcoded name only if necessary.

Definition at line 428 of file ParameterHandlerUtils.h.

                                                                     {
// *************************************
  if (!CovarianceFolder) {
    MACH3LOG_ERROR("Null TDirectory passed to {}", __func__);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  TMacro* foundMacro = nullptr;
  int macroCount = 0;
 
  TIter next(CovarianceFolder->GetListOfKeys());
  TKey* key;
  while ((key = dynamic_cast<TKey*>(next()))) {
    if (std::string(key->GetClassName()) == "TMacro") {
      ++macroCount;
      if (macroCount == 1) {
        foundMacro = dynamic_cast<TMacro*>(key->ReadObj());
      }
    }
  }
 
  if (macroCount == 1 && foundMacro) {
    MACH3LOG_INFO("Found single TMacro in directory: using it.");
    return foundMacro;
  } else {
    MACH3LOG_WARN("Found {} TMacro objects. Using hardcoded macro name: Config_xsec_cov.", macroCount);
    TMacro* fallback = CovarianceFolder->Get<TMacro>("Config_xsec_cov");
    if (!fallback) {
      MACH3LOG_WARN("Fallback macro 'Config_xsec_cov' not found in directory.");
    }
    return fallback;
  }
}

GetCovMatrixFromChain()

TMatrixDSym* M3::GetCovMatrixFromChain ( TDirectory *  TempFile )

inline

KS: Retrieve the cross-section covariance matrix from the given TDirectory. Historically, multiple covariance matrices could be stored, requiring fragile hardcoded paths like "CovarianceFolder/xsec_cov". This function maintains backward compatibility by:

• Using the single matrix present if only one exists,
• Otherwise falling back to the hardcoded path. This avoids error-prone assumptions while supporting both old and new formats.

Definition at line 469 of file ParameterHandlerUtils.h.

                                                                {
// *************************************
  if (!TempFile) {
    MACH3LOG_ERROR("Null TDirectory passed to {}.", __func__);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  TMatrixDSym* foundMatrix = nullptr;
  int matrixCount = 0;
 
  TIter next(TempFile->GetListOfKeys());
  TKey* key;
  while ((key = dynamic_cast<TKey*>(next()))) {
    std::string className = key->GetClassName();
    if (className.find("TMatrix") != std::string::npos) {
      ++matrixCount;
      if (matrixCount == 1) {
        foundMatrix = dynamic_cast<TMatrixDSym*>(key->ReadObj());
      }
    }
  }
 
  if (matrixCount == 1 && foundMatrix) {
    MACH3LOG_INFO("Found single TMatrixDSym in directory: using it.");
    return foundMatrix;
  } else {
    MACH3LOG_WARN("Found {} TMatrixDSym objects. Using hardcoded path: xsec_cov.", matrixCount);
    TMatrixDSym* fallback = TempFile->Get<TMatrixDSym>("xsec_cov");
    if (!fallback) {
      MACH3LOG_WARN("Fallback matrix 'xsec_cov' not found.");
    }
    return fallback;
  }
}

GetNThreads()

int M3::GetNThreads

(

)

number of threads which we need for example for TRandom3

Definition at line 372 of file Monitor.cpp.

                  {
// ***************************************************************************
  #ifdef MULTITHREAD
  return omp_get_max_threads();
  #else
  return 1;
  #endif
}

GetThreadIndex()

int M3::GetThreadIndex ( )

inline

thread index inside parallel loop

Definition at line 86 of file Monitor.h.

                              {
    #ifdef MULTITHREAD
    return omp_get_thread_num();
    #else
    return 0;
    #endif
  }

MakeCorrelationMatrix()

void M3::MakeCorrelationMatrix ( YAML::Node &  root, const std::vector< double > &  Values, const std::vector< double > &  Errors, const std::vector< std::vector< double >> &  Correlation, const std::string &  OutYAMLName, const std::vector< std::string > &  FancyNames = {}  )

inline

KS: Replace correlation matrix and tune values in YAML covariance matrix.

Parameters

root	The root YAML node to be updated.
Values	The new values for each systematic.
Errors	The new errors for each systematic.
Correlation	The new correlation matrix (must be square and match Values size).
OutYAMLName	The output filename for the updated YAML.
FancyNames	Optional list of fancy names to match systematics (must match Values size if provided).

Definition at line 314 of file ParameterHandlerUtils.h.

                                                                           {}) {
// *************************************
  if (Values.size() != Errors.size() || Values.size() != Correlation.size()) {
    MACH3LOG_ERROR("Size mismatch between Values, Errors, and Correlation matrix");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  for (const auto& row : Correlation) {
    if (row.size() != Correlation.size()) {
      MACH3LOG_ERROR("Correlation matrix is not square");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
  }
 
  YAML::Node NodeCopy = YAML::Clone(root);
  YAML::Node systematics = NodeCopy["Systematics"];
 
  if (!systematics || !systematics.IsSequence()) {
    MACH3LOG_ERROR("'Systematics' node is missing or not a sequence");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  if (!FancyNames.empty() && FancyNames.size() != Values.size()) {
    MACH3LOG_ERROR("FancyNames size ({}) does not match Values size ({})", FancyNames.size(), Values.size());
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  // Map from FancyName to Systematic node
  std::unordered_map<std::string, YAML::Node> nameToNode;
  for (std::size_t i = 0; i < systematics.size(); ++i) {
    YAML::Node syst = systematics[i]["Systematic"];
    if (!syst || !syst["Names"] || !syst["Names"]["FancyName"]) continue;
    std::string name = syst["Names"]["FancyName"].as<std::string>();
    nameToNode[name] = syst;
  }
 
  if (!FancyNames.empty()) {
    for (std::size_t i = 0; i < FancyNames.size(); ++i) {
      const std::string& name_i = FancyNames[i];
      auto it_i = nameToNode.find(name_i);
      if (it_i == nameToNode.end()) {
        MACH3LOG_ERROR("Could not find FancyName '{}' in YAML", name_i);
        throw MaCh3Exception(__FILE__, __LINE__);
      }
      YAML::Node& syst_i = it_i->second;
 
      syst_i["ParameterValues"]["PreFitValue"] = M3::Utils::FormatDouble(Values[i], 4);
      syst_i["Error"] = M3::Utils::FormatDouble(Errors[i], 4);
 
      YAML::Node correlationsNode = YAML::Node(YAML::NodeType::Sequence);
      for (std::size_t j = 0; j < FancyNames.size(); ++j) {
        if (i == j) continue;
        // KS: Skip if value close to 0
        if (std::abs(Correlation[i][j]) < 1e-8) continue;
        YAML::Node singleEntry;
        singleEntry[FancyNames[j]] = M3::Utils::FormatDouble(Correlation[i][j], 4);
        correlationsNode.push_back(singleEntry);
      }
      syst_i["Correlations"] = correlationsNode;
    }
  } else {
    if (systematics.size() != Values.size()) {
      MACH3LOG_ERROR("Mismatch in sizes: Values has {}, but YAML 'Systematics' has {} entries",
                     Values.size(), systematics.size());
      throw MaCh3Exception(__FILE__, __LINE__);
    }
 
    for (std::size_t i = 0; i < systematics.size(); ++i) {
      YAML::Node syst = systematics[i]["Systematic"];
      if (!syst) {
        MACH3LOG_ERROR("Missing 'Systematic' node at index {}", i);
        throw MaCh3Exception(__FILE__, __LINE__);
      }
 
      syst["ParameterValues"]["PreFitValue"] = M3::Utils::FormatDouble(Values[i], 4);
      syst["Error"] = M3::Utils::FormatDouble(Errors[i], 4);
 
      YAML::Node correlationsNode = YAML::Node(YAML::NodeType::Sequence);
      for (std::size_t j = 0; j < Correlation[i].size(); ++j) {
        if (i == j) continue;
        // KS: Skip if value close to 0
        if (std::abs(Correlation[i][j]) < 1e-8) continue;
        YAML::Node singleEntry;
        const std::string& otherName = systematics[j]["Systematic"]["Names"]["FancyName"].as<std::string>();
        singleEntry[otherName] = M3::Utils::FormatDouble(Correlation[i][j], 4);
        correlationsNode.push_back(singleEntry);
      }
      syst["Correlations"] = correlationsNode;
    }
  }
 
  // Convert and write
  std::string YAMLString = YAMLtoSTRING(NodeCopy);
  FixSampleNamesQuotes(YAMLString);
  std::ofstream outFile(OutYAMLName);
  if (!outFile) {
    MACH3LOG_ERROR("Failed to open file for writing: {}", OutYAMLName);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  outFile << YAMLString;
  outFile.close();
}

MakeMatrixPosDef()

int M3::MakeMatrixPosDef ( TMatrixDSym *  cov )

inline

Makes sure that matrix is positive-definite by adding a small number to on-diagonal elements.

Definition at line 552 of file ParameterHandlerUtils.h.

                                              {
// *************************************
  //DB Save original warning state and then increase it in this function to suppress 'matrix not positive definite' messages
  //Means we no longer need to overload
  int originalErrorWarning = gErrorIgnoreLevel;
  gErrorIgnoreLevel = kFatal;
 
  //DB Loop 1000 times adding 1e-9 which tops out at 1e-6 shift on the diagonal before throwing error
  constexpr int MaxAttempts = 1e5;
  const int matrixSize = cov->GetNrows();
  int iAttempt = 0;
  bool CanDecomp = false;
  
  // HW: We'll store the diagonal first to prevent inflating the matrix too much!
  std::vector<double> original_diagonal(matrixSize, 0.0);
  for (int iVar = 0 ; iVar < matrixSize; iVar++) {
    original_diagonal[iVar] = (*cov)(iVar, iVar);
  }
 
  int attempts = 0;
 
  for (iAttempt = 0; iAttempt < MaxAttempts; iAttempt++) {
    if (CanDecomposeMatrix(*cov)) {
      CanDecomp = true;
      attempts = iAttempt;
      break;
    }
    else {
      #ifdef MULTITHREAD
      #pragma omp parallel for
      #endif 
      for (int iVar = 0 ; iVar < matrixSize; iVar++) {
        if( (*cov)(iVar, iVar)/10 > original_diagonal[iVar]) {
          MACH3LOG_DEBUG("Diagonal element {} has been shifted too much (> original value/10). Stopping further shifts.", iVar);
          original_diagonal[iVar] = 0; // Prevent further shifts for this element
          continue;
        }
        (*cov)(iVar, iVar) += 1e-9;
      }
    }
  }
 
  if (!CanDecomp) {
    MACH3LOG_ERROR("Tried {} times to shift diagonal but still can not decompose the matrix", MaxAttempts);
    MACH3LOG_ERROR("This indicates that something is wrong with the input matrix");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  //DB Resetting warning level
  gErrorIgnoreLevel = originalErrorWarning;
  return attempts;
}

MatrixMult() [1/3]

double** M3::MatrixMult ( double **  A, double **  B, int  n  )

inline

CW: Multi-threaded matrix multiplication.

Definition at line 108 of file ParameterHandlerUtils.h.

                                                          {
  // First make into monolithic array
  double *A_mon = new double[n*n];
  double *B_mon = new double[n*n];
 
  #ifdef MULTITHREAD
  #pragma omp parallel for
  #endif
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j) {
      A_mon[i*n+j] = A[i][j];
      B_mon[i*n+j] = B[i][j];
    }
  }
  //CW: Now call the monolithic calculator
  double *C_mon = MatrixMult(A_mon, B_mon, n);
  delete A_mon;
  delete B_mon;
 
  // Return the double pointer
  double **C = new double*[n];
  #ifdef MULTITHREAD
  #pragma omp parallel for
  #endif
  for (int i = 0; i < n; ++i) {
    C[i] = new double[n];
    for (int j = 0; j < n; ++j) {
      C[i][j] = C_mon[i*n+j];
    }
  }
  delete C_mon;
 
  return C;
}

MatrixMult() [2/3]

double* M3::MatrixMult ( double *  A, double *  B, int  n  )

inline

CW: Multi-threaded matrix multiplication.

Definition at line 76 of file ParameterHandlerUtils.h.

                                                       {
  //CW: First transpose to increse cache hits
  double *BT = new double[n*n];
  #ifdef MULTITHREAD
  #pragma omp parallel for
  #endif
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
      BT[j*n+i] = B[i*n+j];
    }
  }
 
  // Now multiply
  double *C = new double[n*n];
  #ifdef MULTITHREAD
  #pragma omp parallel for
  #endif
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
      double sum = 0;
      for (int k = 0; k < n; k++) {
        sum += A[i*n+k]*BT[j*n+k];
      }
      C[i*n+j] = sum;
    }
  }
  delete BT;
 
  return C;
}

MatrixMult() [3/3]

TMatrixD M3::MatrixMult ( TMatrixD  A, TMatrixD  B  )

inline

CW: Multi-threaded matrix multiplication.

Definition at line 144 of file ParameterHandlerUtils.h.

{
  double *C_mon = MatrixMult(A.GetMatrixArray(), B.GetMatrixArray(), A.GetNcols());
  TMatrixD C;
  C.Use(A.GetNcols(), A.GetNrows(), C_mon);
  return C;
}

MatrixVectorMulti()

void M3::MatrixVectorMulti ( double *_restrict_  VecMulti, double **_restrict_  matrix, const double *_restrict_  vector, const int  n  )

inline

KS: Custom function to perform multiplication of matrix and vector with multithreading.

Parameters

VecMulti	Output Vector, VecMulti = matrix x vector
matrix	This matrix is used for multiplication VecMulti = matrix x vector
vector	This vector is used for multiplication VecMulti = matrix x vector
n	this is size of matrix and vector, we assume matrix is symmetric

Definition at line 158 of file ParameterHandlerUtils.h.

                                                                                                                                     {
// ********************************************
  #ifdef MULTITHREAD
  #pragma omp parallel for
  #endif
  for (int i = 0; i < n; ++i)
  {
    double result = 0.0;
    #ifdef MULTITHREAD
    #pragma omp simd
    #endif
    for (int j = 0; j < n; ++j)
    {
      result += matrix[i][j]*vector[j];
    }
    VecMulti[i] = result;
  }
}

MatrixVectorMultiSingle()

double M3::MatrixVectorMultiSingle ( double **_restrict_  matrix, const double *_restrict_  vector, const int  Length, const int  i  )

inline

KS: Custom function to perform multiplication of matrix and single element which is thread safe.

Parameters

matrix	This matrix is used for multiplication VecMulti = matrix x vector
vector	This vector is used for multiplication VecMulti = matrix x vector
Length	this is size of matrix and vector, we assume matrix is symmetric
i	Element of matrix that we want to multiply

Definition at line 183 of file ParameterHandlerUtils.h.

                                                                                                                                  {
// ********************************************
  double Element = 0.0;
  #ifdef MULTITHREAD
  #pragma omp simd
  #endif
  for (int j = 0; j < Length; ++j) {
    Element += matrix[i][j]*vector[j];
  }
  return Element;
}

Open()

TFile * M3::Open	(	const std::string &	Name,
		const std::string &	Type,
		const std::string &	File,
		const int	Line
	)

Opens a ROOT file with the given name and mode.

This function wraps ROOT’s TFile constructor and checks whether the file was opened successfully and give some useful debugging information

Parameters

Name	The name or path of the file to open.
Type	The file open mode (e.g., "READ", "RECREATE", "UPDATE").
File	The name of the file where the exception occurred.
Line	The line number where the exception occurred.

Returns

Pointer to the opened ROOT file.

Definition at line 782 of file HistogramUtils.cpp.

                                                                                                 {
// **************************************************************************
  TFile* OutFile = new TFile(Name.c_str(), Type.c_str());
 
  // Check if the file is successfully opened and usable
  if (OutFile->IsZombie()) {
    MACH3LOG_ERROR("Failed to open file: {}", Name);
    std::string lowerType = Type;
    std::transform(lowerType.begin(), lowerType.end(), lowerType.begin(), ::tolower);
    if (lowerType == "recreate") {
      MACH3LOG_ERROR("Check if directory exist");
    }
    delete OutFile;
    throw MaCh3Exception(File, Line);
  }
  return OutFile;
}

PolyToYaml()

YAML::Node M3::PolyToYaml	(	TH2Poly *	Hist,
		const std::string &	YamlName,
		const std::string &	File,
		const int	Line
	)

KS: Convert TH2Poly into yaml config accepted by MaCh3.

Definition at line 891 of file HistogramUtils.cpp.

                                                                                                     {
// ***************************************************************************
  if (!Hist) {
    MACH3LOG_ERROR("Null TH2Poly pointer");
    throw MaCh3Exception(File, Line);
  }
 
  YAML::Node bins(YAML::NodeType::Sequence);
  bins.SetStyle(YAML::EmitterStyle::Flow);
  const int NBins = Hist->GetNumberOfBins();
 
  for (int j = 1; j <= NBins; j++)
  {
    TH2PolyBin* polybin = static_cast<TH2PolyBin*>(Hist->GetBins()->At(j - 1));
 
    double xmin = polybin->GetXMin();
    double xmax = polybin->GetXMax();
    double ymin = polybin->GetYMin();
    double ymax = polybin->GetYMax();
 
    YAML::Node xNode(YAML::NodeType::Sequence);
    xNode.SetStyle(YAML::EmitterStyle::Flow);
    xNode.push_back(xmin);
    xNode.push_back(xmax);
 
    YAML::Node yNode(YAML::NodeType::Sequence);
    yNode.SetStyle(YAML::EmitterStyle::Flow);
    yNode.push_back(ymin);
    yNode.push_back(ymax);
 
    YAML::Node bin(YAML::NodeType::Sequence);
    bin.SetStyle(YAML::EmitterStyle::Flow);
    bin.push_back(xNode);
    bin.push_back(yNode);
 
    bins.push_back(bin);
  }
 
  YAML::Node result;
  result[YamlName] = bins;
 
  return result;
}

ScaleHistogram()

void M3::ScaleHistogram	(	TH1 *	Sample_Hist,
		const double	scale
	)

Scale histogram to get divided by bin width.

Definition at line 801 of file HistogramUtils.cpp.

                                                          {
// ***************************************************************************
  for (int j = 0; j <= Sample_Hist->GetNbinsX(); ++j)
  {
    double num = Sample_Hist->GetBinContent(j);
    double numErr = Sample_Hist->GetBinError(j);
    double den = Sample_Hist->GetBinWidth(j);
    double value = 0.;
    double valueErr = 0.;
    if (den != 0)
    {
      value = num/(den/scale);
      valueErr = numErr/(den/scale);
      Sample_Hist->SetBinContent(j,value);
      Sample_Hist->SetBinError(j,valueErr);
    }
  }
}

Variable Documentation

_BAD_DOUBLE_

constexpr static const double M3::_BAD_DOUBLE_ = -999.99

staticconstexpr

Default value used for double initialisation.

Definition at line 53 of file Core.h.

_BAD_INT_

constexpr static const int M3::_BAD_INT_ = -999

staticconstexpr

Default value used for int initialisation.

Definition at line 55 of file Core.h.

_DEFAULT_RETURN_VAL_

constexpr static const double M3::_DEFAULT_RETURN_VAL_ = -999999.123456

staticconstexpr

Definition at line 56 of file Core.h.

_LARGE_LOGL_

constexpr static const double M3::_LARGE_LOGL_ = 1234567890.0

staticconstexpr

Large Likelihood is used it parameter go out of physical boundary, this indicates in MCMC that such step should be removed.

Definition at line 80 of file Core.h.

_LOW_MC_BOUND_

constexpr static const double M3::_LOW_MC_BOUND_ = .00001

staticconstexpr

MC prediction lower bound in bin to identify problematic binning definitions and handle LogL calculation.

Definition at line 83 of file Core.h.

DefSplineKnotLowBound

constexpr static const double M3::DefSplineKnotLowBound = -9999

staticconstexpr

Default value for spline knot capping, default mean not capping is being applied.

Definition at line 88 of file Core.h.

DefSplineKnotUpBound

constexpr static const double M3::DefSplineKnotUpBound = 9999

staticconstexpr

Default value for spline knot capping, default mean not capping is being applied.

Definition at line 86 of file Core.h.

float_t_str_repr

constexpr static const char* M3::float_t_str_repr = "D"

staticconstexpr

Definition at line 40 of file Core.h.

KinematicLowBound

constexpr static const double M3::KinematicLowBound = std::numeric_limits<double>::lowest()

staticconstexpr

When parameter has no bound this serves as it. Lowest possible value the system.

Definition at line 75 of file Core.h.

KinematicUpBound

constexpr static const double M3::KinematicUpBound = std::numeric_limits<double>::max()

staticconstexpr

When parameter has no bound this serves as it. Highest possible value the system.

Definition at line 77 of file Core.h.

UnderOverFlowBin

constexpr static const int M3::UnderOverFlowBin = -1

staticconstexpr

Mark bin which is overflow or underflow in MaCh3 binning.

Definition at line 91 of file Core.h.

Unity

constexpr static const float_t M3::Unity = Unity_D

staticconstexpr

Definition at line 64 of file Core.h.

Unity_D

constexpr static const double M3::Unity_D = 1.

staticconstexpr

Some commonly used variables to which we set pointers to.

Definition at line 59 of file Core.h.

Unity_F

constexpr static const float M3::Unity_F = 1.

staticconstexpr

Definition at line 60 of file Core.h.

Zero

constexpr static const float_t M3::Zero = Zero_D

staticconstexpr

Definition at line 71 of file Core.h.

Zero_D

constexpr static const double M3::Zero_D = 0.

staticconstexpr

Definition at line 66 of file Core.h.

Zero_F

constexpr static const float M3::Zero_F = 0.

staticconstexpr

Definition at line 67 of file Core.h.