diff --git a/covariance/covarianceBase.cpp b/covariance/covarianceBase.cpp index e6033ac1..d0b7adb1 100644 --- a/covariance/covarianceBase.cpp +++ b/covariance/covarianceBase.cpp @@ -1,7 +1,7 @@ #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); @@ -9,7 +9,7 @@ covarianceBase::covarianceBase(const char *name, const char *file) : inputFile(s LastPCAdpar = -999; } // ******************************************** -covarianceBase::covarianceBase(const std::vector& 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& 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 "); @@ -105,11 +105,11 @@ 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); @@ -117,7 +117,7 @@ void covarianceBase::init(const char *name, const char *file) { } // 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); @@ -1319,7 +1319,7 @@ std::vector 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++) { diff --git a/covariance/covarianceBase.h b/covariance/covarianceBase.h index ab4f207f..73ec3523 100644 --- a/covariance/covarianceBase.h +++ b/covariance/covarianceBase.h @@ -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& YAMLFile, const char *name, double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999); + covarianceBase(const std::vector& 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(); @@ -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 @@ -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; } @@ -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& YAMLFile); @@ -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 diff --git a/covariance/covarianceOsc.cpp b/covariance/covarianceOsc.cpp index b4b20f03..f5aeb009 100644 --- a/covariance/covarianceOsc.cpp +++ b/covariance/covarianceOsc.cpp @@ -1,7 +1,7 @@ #include "covariance/covarianceOsc.h" // ************************************* -covarianceOsc::covarianceOsc(const std::vector& YAMLFile, const char *name, double threshold, int FirstPCA, int LastPCA) +covarianceOsc::covarianceOsc(const std::vector& YAMLFile, std::string name, double threshold, int FirstPCA, int LastPCA) : covarianceBase(YAMLFile, name, threshold, FirstPCA, LastPCA){ // ************************************* @@ -25,12 +25,7 @@ covarianceOsc::covarianceOsc(const std::vector& 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(); diff --git a/covariance/covarianceOsc.h b/covariance/covarianceOsc.h index d13381ca..43e98f16 100644 --- a/covariance/covarianceOsc.h +++ b/covariance/covarianceOsc.h @@ -8,7 +8,7 @@ class covarianceOsc : public covarianceBase { public: /// @brief Constructor - covarianceOsc(const std::vector& FileNames, const char *name = "osc_cov", double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999); + covarianceOsc(const std::vector& 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 diff --git a/covariance/covarianceXsec.cpp b/covariance/covarianceXsec.cpp index 547645fe..fefae18c 100644 --- a/covariance/covarianceXsec.cpp +++ b/covariance/covarianceXsec.cpp @@ -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& YAMLFile, const char *name, double threshold, int FirstPCA, int LastPCA) +covarianceXsec::covarianceXsec(const std::vector& YAMLFile, std::string name, double threshold, int FirstPCA, int LastPCA) : covarianceBase(YAMLFile, name, threshold, FirstPCA, LastPCA){ // ******************************************** InitXsecFromConfig(); diff --git a/covariance/covarianceXsec.h b/covariance/covarianceXsec.h index 487d5b97..0936cc36 100644 --- a/covariance/covarianceXsec.h +++ b/covariance/covarianceXsec.h @@ -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& FileNames, const char *name = "xsec_cov", double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999); + covarianceXsec(const std::vector& FileNames, std::string name = "xsec_cov", double threshold = -1, int FirstPCAdpar = -999, int LastPCAdpar = -999); /// @brief Destructor ~covarianceXsec(); @@ -62,27 +62,12 @@ class covarianceXsec : public covarianceBase { /// @brief DB Grab the Spline Modes for the relevant DetID const std::vector< std::vector > GetSplineModeVecFromDetID(const int DetID); - /// @brief DB Grab the Spline Indices for the relevant DetID - const std::vector 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 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 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 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 GetFuncParsNamesFromDetID(const int DetID){return GetParsNamesFromDetID(DetID, SystType::kFunc);} - /// @brief DB Grab the Functional parameter indices for the relevant DetID - const std::vector 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 getNominalArray() override diff --git a/mcmc/FitterBase.cpp b/mcmc/FitterBase.cpp index 229ea786..96751fa3 100644 --- a/mcmc/FitterBase.cpp +++ b/mcmc/FitterBase.cpp @@ -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; @@ -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"); @@ -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()); @@ -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 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 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() { diff --git a/mcmc/FitterBase.h b/mcmc/FitterBase.h index de875f7d..14d9d365 100644 --- a/mcmc/FitterBase.h +++ b/mcmc/FitterBase.h @@ -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" @@ -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; @@ -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 @@ -118,9 +110,6 @@ class FitterBase { /// Systematic holder std::vector systematics; - /// handles oscillation parameters - covarianceOsc *osc; - /// tells global time how long fit took TStopwatch* clock; /// tells how long single step/fit iteration took diff --git a/mcmc/LikelihoodFit.cpp b/mcmc/LikelihoodFit.cpp index 9b8b34ef..cdc8696b 100644 --- a/mcmc/LikelihoodFit.cpp +++ b/mcmc/LikelihoodFit.cpp @@ -32,10 +32,6 @@ void LikelihoodFit::PrepareFit() { NParsPCA += systematics[s]->getNpars(); } - if (osc) { - std::cerr<<" Osc not supported "<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 diff --git a/mcmc/mcmc.cpp b/mcmc/mcmc.cpp index 1be703e5..bc3c582b 100644 --- a/mcmc/mcmc.cpp +++ b/mcmc/mcmc.cpp @@ -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(); @@ -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 diff --git a/python/fitter.cpp b/python/fitter.cpp index 91d31796..403ee308 100644 --- a/python/fitter.cpp +++ b/python/fitter.cpp @@ -128,14 +128,6 @@ void initFitter(py::module &m){ py::arg("cov") ) - .def( - "add_osc_handler", - py::overload_cast(&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_(m_fitter, "MCMC") @@ -193,4 +185,4 @@ void initFitter(py::module &m){ ; // end of PSO class binding -} \ No newline at end of file +} diff --git a/samplePDF/samplePDFFDBase.cpp b/samplePDF/samplePDFFDBase.cpp index 0548ad33..2f69dca5 100644 --- a/samplePDF/samplePDFFDBase.cpp +++ b/samplePDF/samplePDFFDBase.cpp @@ -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; } diff --git a/splines/splineFDBase.cpp b/splines/splineFDBase.cpp index 99784a32..a8075393 100644 --- a/splines/splineFDBase.cpp +++ b/splines/splineFDBase.cpp @@ -54,11 +54,11 @@ bool splineFDBase::AddSample(std::string SampleName, int DetID, std::vectorGetNumSplineParamsFromDetID(DetID); + int nSplineParam = xsec->GetNumParamsFromDetID(DetID, SystType::kSpline); nSplineParams.push_back(nSplineParam); //This holds the global index of the spline i.e. 0 -> _fNumPar - std::vector GlobalSystIndex_Sample = xsec->GetGlobalSystIndexFromDetID(DetID, kSpline); + std::vector GlobalSystIndex_Sample = xsec->GetGlobalSystIndexFromDetID(DetID, SystType::kSpline); //Keep track of this for all the samples GlobalSystIndex.push_back(GlobalSystIndex_Sample);