Skip to content

Commit

Permalink
Merge pull request #183 from mach3-software/feature_RemoveGarbage
Browse files Browse the repository at this point in the history
Remove Garbage
  • Loading branch information
dbarrow257 authored Oct 25, 2024
2 parents df81585 + b1d42b7 commit 210ae41
Show file tree
Hide file tree
Showing 13 changed files with 26 additions and 139 deletions.
12 changes: 6 additions & 6 deletions covariance/covarianceBase.cpp
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
#include "covariance/covarianceBase.h"

// ********************************************
covarianceBase::covarianceBase(const char *name, const char *file) : inputFile(std::string(file)), pca(false) {
covarianceBase::covarianceBase(std::string name, std::string file) : inputFile(std::string(file)), pca(false) {
// ********************************************
MACH3LOG_INFO("Constructing instance of covarianceBase");
init(name, file);
FirstPCAdpar = -999;
LastPCAdpar = -999;
}
// ********************************************
covarianceBase::covarianceBase(const std::vector<std::string>& YAMLFile, const char *name, double threshold, int FirstPCA, int LastPCA) : inputFile(YAMLFile[0].c_str()), matrixName(name), pca(true), eigen_threshold(threshold), FirstPCAdpar(FirstPCA), LastPCAdpar(LastPCA) {
covarianceBase::covarianceBase(const std::vector<std::string>& YAMLFile, std::string name, double threshold, int FirstPCA, int LastPCA) : inputFile(YAMLFile[0].c_str()), matrixName(name), pca(true), eigen_threshold(threshold), FirstPCAdpar(FirstPCA), LastPCAdpar(LastPCA) {
// ********************************************

MACH3LOG_INFO("Constructing instance of covarianceBase using ");
Expand Down Expand Up @@ -105,19 +105,19 @@ void covarianceBase::ConstructPCA() {
}

// ********************************************
void covarianceBase::init(const char *name, const char *file) {
void covarianceBase::init(std::string name, std::string file) {
// ********************************************

// Set the covariance matrix from input ROOT file (e.g. flux, ND280, NIWG)
TFile *infile = new TFile(file, "READ");
TFile *infile = new TFile(file.c_str(), "READ");
if (infile->IsZombie()) {
MACH3LOG_ERROR("Could not open input covariance ROOT file {} !!!", file);
MACH3LOG_ERROR("Was about to retrieve matrix with name {}", name);
throw MaCh3Exception(__FILE__ , __LINE__ );
}

// Should put in a
TMatrixDSym *CovMat = (TMatrixDSym*)(infile->Get(name));
TMatrixDSym *CovMat = (TMatrixDSym*)(infile->Get(name.c_str()));
if (CovMat == nullptr) {
MACH3LOG_ERROR("Could not find covariance matrix name {} in file {}", name, file);
MACH3LOG_ERROR("Are you really sure {} exists in the file?", name);
Expand Down Expand Up @@ -1319,7 +1319,7 @@ std::vector<double> covarianceBase::getNominalArray() {
// KS: Convert covariance matrix to correlation matrix and return TH2D which can be used for fancy plotting
TH2D* covarianceBase::GetCorrelationMatrix() {
// ********************************************
TH2D* hMatrix = new TH2D(getName(), getName(), _fNumPar, 0.0, _fNumPar, _fNumPar, 0.0, _fNumPar);
TH2D* hMatrix = new TH2D(getName().c_str(), getName().c_str(), _fNumPar, 0.0, _fNumPar, _fNumPar, 0.0, _fNumPar);

for(int i = 0; i < _fNumPar; i++)
{
Expand Down
20 changes: 7 additions & 13 deletions covariance/covarianceBase.h
Original file line number Diff line number Diff line change
Expand Up @@ -22,11 +22,11 @@ class covarianceBase {
/// @param threshold PCA threshold from 0 to 1. Default is -1 and means no PCA
/// @param FirstPCAdpar First PCA parameter that will be decomposed.
/// @param LastPCAdpar First PCA parameter that will be decomposed.
covarianceBase(const std::vector<std::string>& YAMLFile, const char *name, double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999);
covarianceBase(const std::vector<std::string>& YAMLFile, std::string name, double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999);
/// @brief "Usual" constructors from root file
/// @param name Matrix name
/// @param file Path to matrix root file
covarianceBase(const char *name, const char *file);
covarianceBase(std::string name, std::string file);

/// @brief Destructor
virtual ~covarianceBase();
Expand All @@ -37,11 +37,11 @@ class covarianceBase {
/// @param cov Covariance matrix which we set and will be used later for evaluation of penalty term
void setCovMatrix(TMatrixDSym *cov);
/// @brief Set matrix name
void setName(const char *name) { matrixName = name; }
void setName(std::string name) { matrixName = name; }
/// @brief change parameter name
/// @param i Parameter index
/// @param name new name which will be set
void setParName(int i, char *name) { _fNames.at(i) = std::string(name); }
void setParName(int i, std::string name) { _fNames.at(i) = name; }
void setSingleParameter(const int parNo, const double parVal);
/// @brief Set all the covariance matrix parameters to a user-defined value
/// @param i Parameter index
Expand Down Expand Up @@ -111,19 +111,13 @@ class covarianceBase {
inline bool getFlatPrior(const int i) { return _fFlatPrior[i]; }

/// @brief Get name of covariance
const char *getName() { return matrixName; }
std::string getName() { return matrixName; }
/// @brief Get name of covariance
/// @param i Parameter index
std::string GetParName(const int i) {return _fNames[i];}
/// @brief Get name of the Parameter
/// @param i Parameter index
const char* GetParName(const int i) const { return _fNames[i].c_str(); }
/// @brief Get fancy name of the Parameter
/// @param i Parameter index
std::string GetParFancyName(const int i) {return _fFancyNames[i];}
/// @brief Get fancy name of the Parameter
/// @param i Parameter index
const char* GetParFancyName(const int i) const { return _fFancyNames[i].c_str(); }
/// @brief Get name of input file
std::string getInputFile() const { return inputFile; }

Expand Down Expand Up @@ -351,7 +345,7 @@ class covarianceBase {
YAML::Node GetConfig() const { return _fYAMLDoc; }
protected:
/// @brief Initialisation of the class using matrix from root file
void init(const char *name, const char *file);
void init(std::string name, std::string file);
/// @brief Initialisation of the class using config
/// @param YAMLFile A vector of strings representing the YAML files used for initialisation of matrix
void init(const std::vector<std::string>& YAMLFile);
Expand Down Expand Up @@ -396,7 +390,7 @@ class covarianceBase {
/// Total number of params, deprecated, please don't use it
int size;
/// Name of cov matrix
const char *matrixName;
std::string matrixName;
/// The covariance matrix
TMatrixDSym *covMatrix;
/// The inverse covariance matrix
Expand Down
9 changes: 2 additions & 7 deletions covariance/covarianceOsc.cpp
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#include "covariance/covarianceOsc.h"

// *************************************
covarianceOsc::covarianceOsc(const std::vector<std::string>& YAMLFile, const char *name, double threshold, int FirstPCA, int LastPCA)
covarianceOsc::covarianceOsc(const std::vector<std::string>& YAMLFile, std::string name, double threshold, int FirstPCA, int LastPCA)
: covarianceBase(YAMLFile, name, threshold, FirstPCA, LastPCA){
// *************************************

Expand All @@ -25,12 +25,7 @@ covarianceOsc::covarianceOsc(const std::vector<std::string>& YAMLFile, const cha
_fPropVal[io] = _fCurrVal[io];
}

/// WARNING HARDCODED
MACH3LOG_CRITICAL("Fixing baseline and density, it is hardcoded sry");
toggleFixParameter("baseline");
toggleFixParameter("density");

/// @todo KS: Technically if we would like to use PCA we have ot initialise parts here...
/// @todo KS: Technically if we would like to use PCA we have to initialise parts here...
flipdelM = false;

randomize();
Expand Down
2 changes: 1 addition & 1 deletion covariance/covarianceOsc.h
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ class covarianceOsc : public covarianceBase
{
public:
/// @brief Constructor
covarianceOsc(const std::vector<std::string>& FileNames, const char *name = "osc_cov", double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999);
covarianceOsc(const std::vector<std::string>& FileNames, std::string name = "osc_cov", double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999);
/// @brief Destructor
virtual ~covarianceOsc();
/// @brief Checks if parameter didn't go outside of physical bounds
Expand Down
2 changes: 1 addition & 1 deletion covariance/covarianceXsec.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
// ETA - YAML constructor
// this will replace the root file constructor but let's keep it in
// to do some validations
covarianceXsec::covarianceXsec(const std::vector<std::string>& YAMLFile, const char *name, double threshold, int FirstPCA, int LastPCA)
covarianceXsec::covarianceXsec(const std::vector<std::string>& YAMLFile, std::string name, double threshold, int FirstPCA, int LastPCA)
: covarianceBase(YAMLFile, name, threshold, FirstPCA, LastPCA){
// ********************************************
InitXsecFromConfig();
Expand Down
17 changes: 1 addition & 16 deletions covariance/covarianceXsec.h
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ class covarianceXsec : public covarianceBase {
/// @param threshold PCA threshold from 0 to 1. Default is -1 and means no PCA
/// @param FirstPCAdpar First PCA parameter that will be decomposed.
/// @param LastPCAdpar First PCA parameter that will be decomposed.
covarianceXsec(const std::vector<std::string>& FileNames, const char *name = "xsec_cov", double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999);
covarianceXsec(const std::vector<std::string>& FileNames, std::string name = "xsec_cov", double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999);
/// @brief Destructor
~covarianceXsec();

Expand Down Expand Up @@ -62,27 +62,12 @@ class covarianceXsec : public covarianceBase {

/// @brief DB Grab the Spline Modes for the relevant DetID
const std::vector< std::vector<int> > GetSplineModeVecFromDetID(const int DetID);
/// @brief DB Grab the Spline Indices for the relevant DetID
const std::vector<int> GetSplineParsIndexFromDetID(const int DetID){return GetParsIndexFromDetID(DetID, SystType::kSpline);}
/// @brief ETA Grab the index of the spline relative to the _fSplineNames vector.
const std::vector<int> GetSplineSystIndexFromDetID(const int DetID){return GetSystIndexFromDetID(DetID, SystType::kSpline);};
/// @brief Grab the index of the syst relative to global numbering.
/// @param Type Type of syst, for example kNorm, kSpline etc
const std::vector<int> GetSystIndexFromDetID(const int DetID, const SystType Type);

/// @brief DB Grab the Number of splines for the relevant DetID
int GetNumSplineParamsFromDetID(const int DetID){return GetNumParamsFromDetID(DetID, SystType::kSpline);}

/// @brief DB Get norm/func parameters depending on given DetID
const std::vector<XsecNorms4> GetNormParsFromDetID(const int DetID);

/// @brief DB Grab the number of Normalisation parameters for the relevant DetID
int GetNumFuncParamsFromDetID(const int DetID){return GetNumParamsFromDetID(DetID, SystType::kFunc);}
/// @brief DB Grab the Functional parameter names for the relevant DetID
const std::vector<std::string> GetFuncParsNamesFromDetID(const int DetID){return GetParsNamesFromDetID(DetID, SystType::kFunc);}
/// @brief DB Grab the Functional parameter indices for the relevant DetID
const std::vector<int> GetFuncParsIndexFromDetID(const int DetID){return GetParsIndexFromDetID(DetID, SystType::kFunc);}

/// @brief KS: For most covariances nominal and fparInit (prior) are the same, however for Xsec those can be different
/// For example Sigma Var are done around nominal in ND280, no idea why though...
std::vector<double> getNominalArray() override
Expand Down
41 changes: 1 addition & 40 deletions mcmc/FitterBase.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,6 @@ FitterBase::FitterBase(manager * const man) : fitMan(man) {
// Clear the samples and systematics
samples.clear();
systematics.clear();
osc = nullptr;

sample_llh = nullptr;
syst_llh = nullptr;
Expand Down Expand Up @@ -182,11 +181,6 @@ void FitterBase::PrepareOutput() {
(*it)->SetBranches(*outTree, SaveProposal);
}

if (osc) {
osc->SetBranches(*outTree, SaveProposal);
outTree->Branch("LogL_osc", &osc_llh, "LogL_osc/D");
}

outTree->Branch("LogL", &logLCurr, "LogL/D");
outTree->Branch("accProb", &accProb, "accProb/D");
outTree->Branch("step", &step, "step/I");
Expand Down Expand Up @@ -287,7 +281,7 @@ void FitterBase::addSystObj(covarianceBase * const cov) {
t_vec.Write((std::string(cov->getName()) + "_prior").c_str());
delete[] n_vec;

cov->getCovMatrix()->Write(cov->getName());
cov->getCovMatrix()->Write(cov->getName().c_str());

TH2D* CorrMatrix = cov->GetCorrelationMatrix();
CorrMatrix->Write((cov->getName() + std::string("_Corr")).c_str());
Expand All @@ -306,39 +300,6 @@ void FitterBase::addSystObj(covarianceBase * const cov) {
return;
}

// *************************
// Add an oscillation handler
// Similar to systematic really, but handles the oscillation weights
void FitterBase::addOscHandler(covarianceOsc * const oscf) {
// *************************
osc = oscf;
if (save_nominal) {
CovFolder->cd();
std::vector<double> vec = oscf->getNominalArray();
size_t n = vec.size();
double *n_vec = new double[n];
for (size_t i = 0; i < n; ++i) {
n_vec[i] = vec[i];
}
TVectorT<double> t_vec(n, n_vec);
TString nameof = TString(oscf->getName());
nameof = nameof.Append("_nom");
t_vec.Write(nameof);
delete[] n_vec;
outputFile->cd();
}

// If we have yaml config file for covariance let's save it
YAML::Node Config = oscf->GetConfig();
if(!Config.IsNull())
{
TMacro ConfigSave = YAMLtoTMacro(Config, (std::string("Config_") + oscf->getName()));
ConfigSave.Write();
}
return;
}


// *******************
// Process the MCMC output to get postfit etc
void FitterBase::ProcessMCMC() {
Expand Down
11 changes: 0 additions & 11 deletions mcmc/FitterBase.h
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@
// MaCh3 Includes
#include "samplePDF/samplePDFBase.h"
#include "covariance/covarianceBase.h"
#include "covariance/covarianceOsc.h"
#include "manager/manager.h"
#include "mcmc/MCMCProcessor.h"

Expand Down Expand Up @@ -37,10 +36,6 @@ class FitterBase {
/// @param cov A pointer to a Covariance object derived from covarianceBase.
void addSystObj(covarianceBase* cov);

/// @brief Adds an oscillation handler for covariance objects.
/// @param oscf A pointer to a covarianceOsc object for forward oscillations.
void addOscHandler(covarianceOsc* oscf);

/// @brief 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.
virtual void runMCMC() = 0;

Expand Down Expand Up @@ -102,9 +97,6 @@ class FitterBase {
/// counts accepted steps
int accCount;

/// LLH for samples/syst objects
/// oscillation covariance llh
double osc_llh;
/// store the llh breakdowns
double *sample_llh;
/// systematic llh breakdowns
Expand All @@ -118,9 +110,6 @@ class FitterBase {
/// Systematic holder
std::vector<covarianceBase*> systematics;

/// handles oscillation parameters
covarianceOsc *osc;

/// tells global time how long fit took
TStopwatch* clock;
/// tells how long single step/fit iteration took
Expand Down
12 changes: 1 addition & 11 deletions mcmc/LikelihoodFit.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -32,10 +32,6 @@ void LikelihoodFit::PrepareFit() {
NParsPCA += systematics[s]->getNpars();
}

if (osc) {
std::cerr<<" Osc not supported "<<std::endl;
throw;
}
//KS: If PCA is note enabled NParsPCA == NPars
MACH3LOG_INFO("Total number of parameters {}", NParsPCA);
}
Expand Down Expand Up @@ -112,13 +108,7 @@ double LikelihoodFit::CalcChi2(const double* x) {
// But since sample reweight is multi-threaded it's probably better to do that
for (size_t i = 0; i < samples.size(); i++)
{
// If we're running with different oscillation parameters for neutrino and anti-neutrino
if (osc) {
samples[i]->reweight();
// If we aren't using any oscillation
} else {
samples[i]->reweight();
}
samples[i]->reweight();
}

//DB for atmospheric event by event sample migration, need to fully reweight all samples to allow event passing prior to likelihood evaluation
Expand Down
19 changes: 0 additions & 19 deletions mcmc/mcmc.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -60,10 +60,6 @@ void mcmc::CheckStep() {
if (accept && !reject) {
logLCurr = logLProp;

if (osc) {
osc->acceptStep();
}

// Loop over systematics and accept
for (size_t s = 0; s < systematics.size(); ++s) {
systematics[s]->acceptStep();
Expand Down Expand Up @@ -136,21 +132,6 @@ void mcmc::ProposeStep() {
// Initiate to false
reject = false;

// Propose steps for the oscillation handlers
if (osc) {
osc->proposeStep();

// Now get the likelihoods for the oscillation
osc_llh = osc->GetLikelihood();

// Add the oscillation likelihoods to the reconfigure likelihoods
llh += osc_llh;

#ifdef DEBUG
if (debug) debugFile << "LLH for oscillation handler: " << llh << std::endl;
#endif
}

// Loop over the systematics and propose the initial step
for (size_t s = 0; s < systematics.size(); ++s) {
// Could throw the initial value here to do MCMC stability studies
Expand Down
10 changes: 1 addition & 9 deletions python/fitter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -128,14 +128,6 @@ void initFitter(py::module &m){
py::arg("cov")
)

.def(
"add_osc_handler",
py::overload_cast<covarianceOsc *>(&FitterBase::addOscHandler),
" Adds an oscillation handler for covariance objects. \n"
" :param oscf: An oscillation handler object for dealing with neutrino oscillation calculations. ",
py::arg("osc")
)

; // End of FitterBase class binding

py::class_<mcmc, FitterBase>(m_fitter, "MCMC")
Expand Down Expand Up @@ -193,4 +185,4 @@ void initFitter(py::module &m){

; // end of PSO class binding

}
}
6 changes: 3 additions & 3 deletions samplePDF/samplePDFFDBase.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -741,9 +741,9 @@ void samplePDFFDBase::SetXsecCov(covarianceXsec *xsec){

//DB Now get this information using the DetID from the config
xsec_norms = XsecCov->GetNormParsFromDetID(SampleDetID);
nFuncParams = XsecCov->GetNumFuncParamsFromDetID(SampleDetID);
funcParsNames = XsecCov->GetFuncParsNamesFromDetID(SampleDetID);
funcParsIndex = XsecCov->GetFuncParsIndexFromDetID(SampleDetID);
nFuncParams = XsecCov->GetNumParamsFromDetID(SampleDetID, SystType::kFunc);
funcParsNames = XsecCov->GetParsNamesFromDetID(SampleDetID, SystType::kFunc);
funcParsIndex = XsecCov->GetParsIndexFromDetID(SampleDetID, SystType::kFunc);

return;
}
Expand Down
Loading

0 comments on commit 210ae41

Please sign in to comment.