diff --git a/Diagnostics/README.md b/Diagnostics/README.md index 8e6c36689..07bdfb233 100644 --- a/Diagnostics/README.md +++ b/Diagnostics/README.md @@ -1,4 +1,4 @@ -### Plotting and Diagnostic +### Diagnostic **ProcessMCMC** - The main app you want to use for analysing the ND280 chain. It prints posterior distribution after burn-in the cut. Moreover, you can compare two/three different chains. There are a few options you can modify easily inside the app like selection, burn-in cut, and whether to plot xse+flux or only flux. Other functionality
  1. Produce a covariance matrix with multithreading (which requires lots of RAM due to caching)
    2. diff --git a/covariance/AdaptiveMCMCHandler.h b/covariance/AdaptiveMCMCHandler.h index 3706d56a8..400c71460 100644 --- a/covariance/AdaptiveMCMCHandler.h +++ b/covariance/AdaptiveMCMCHandler.h @@ -1,7 +1,6 @@ #pragma once // MaCh3 Includes -#include "manager/MaCh3Logger.h" #include "manager/manager.h" #include "covariance/CovarianceUtils.h" diff --git a/covariance/CMakeLists.txt b/covariance/CMakeLists.txt index b9af2d5a1..3239a17f6 100644 --- a/covariance/CMakeLists.txt +++ b/covariance/CMakeLists.txt @@ -1,5 +1,4 @@ set(HEADERS - ThrowParms.h CovarianceUtils.h covarianceBase.h covarianceOsc.h @@ -9,7 +8,6 @@ set(HEADERS ) add_library(Covariance SHARED - ThrowParms.cpp covarianceBase.cpp covarianceOsc.cpp covarianceXsec.cpp diff --git a/covariance/PCAHandler.h b/covariance/PCAHandler.h index 80a475952..add35f483 100644 --- a/covariance/PCAHandler.h +++ b/covariance/PCAHandler.h @@ -1,7 +1,6 @@ #pragma once // MaCh3 Includes -#include "manager/MaCh3Logger.h" #include "manager/manager.h" #include "covariance/CovarianceUtils.h" diff --git a/covariance/ThrowParms.cpp b/covariance/ThrowParms.cpp deleted file mode 100755 index 170e3b4ba..000000000 --- a/covariance/ThrowParms.cpp +++ /dev/null @@ -1,104 +0,0 @@ -#include "ThrowParms.h" - -// ************************ -// The constructor -// Taking a TVectorD of central value parameter -// Taking a TMatrixDSym of covariance matrix -ThrowParms::ThrowParms(TVectorD &parms, TMatrixDSym &covm) { -// ************************ - - // Copy over the central value parameters and covariance matrix - npars = parms.GetNrows(); - pvals = new TVectorD(npars); - covar = new TMatrixDSym(npars); - - // Copy the central values and covariance matrix - (*pvals) = parms; - (*covar) = covm; - - // Check that the covariance matrix Cholesky decomposes - TDecompChol chdcmp(*covar); - if (!chdcmp.Decompose()) { - std::cerr << "Cholesky decomposition failed for " << std::endl; - throw; - } - - // Running - chel_dec = new TMatrixD(chdcmp.GetU()); - CheloskyDecomp((*chel_dec)); -} - -ThrowParms::~ThrowParms() -{ - delete pvals; - delete covar; - delete chel_dec; -} - -void ThrowParms::ThrowSet(std::vector &parms) { - // Empty the vector of parameter values that has been passed - if(!parms.empty()) parms.clear(); - - parms.resize(npars); - - int half_pars = npars/2+npars%2; - TVectorD std_rand(npars); - - for(int j=0; j=1.){ - u = 2.*rand.Rndm()-1.; - v = 2.*rand.Rndm()-1.; - s = u*u+v*v; - } - - z[0] = u*sqrt(-2.*TMath::Log(s)/s); - z[1] = v*sqrt(-2.*TMath::Log(s)/s); -} - - -// ************************************ -void ThrowParms::CheloskyDecomp(TMatrixD &chel_mat) { -// ************************************ - - for (int i = 0; i < npars; i++) { - for (int j = 0; j < npars; j++) { - //if diagonal element - if (i == j){ - chel_mat(i, i) = (*covar)(i, i); - for (int k = 0; k <= i-1; k++) { - chel_mat(i,i) = chel_mat(i,i)-pow(chel_mat(i,k),2); - } - chel_mat(i,i) = sqrt(chel_mat(i,i)); - //if lower half - } else if (j < i) { - chel_mat(i, j) = (*covar)(i, j); - for(int k=0; k<=j-1; k++) { - chel_mat(i,j) = chel_mat(i,j)-chel_mat(i,k)*chel_mat(j,k); - } - chel_mat(i,j) = chel_mat(i,j)/chel_mat(j,j); - } else chel_mat(i,j) = 0.; - } // end inner for loop - } // end outer for loop - -} - diff --git a/covariance/ThrowParms.h b/covariance/ThrowParms.h deleted file mode 100755 index 92e8969bf..000000000 --- a/covariance/ThrowParms.h +++ /dev/null @@ -1,38 +0,0 @@ -#pragma once - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -/// @brief No idea really... -/// @warning this is deprecated don't use it -class ThrowParms { - private: - TVectorD *pvals; - TMatrixDSym *covar; - TMatrixD *chel_dec; - int npars; - TRandom3 rand; - - void CheloskyDecomp(TMatrixD &chel_mat); - void StdNormRand(double *z); - - public: - ThrowParms(TVectorD &parms, TMatrixDSym &covm); - ~ThrowParms(); - - void SetSeed (int seed = 0) {rand.SetSeed(seed);}; - int GetSize() {return npars;}; - void ThrowSet(std::vector &parms); - -}; diff --git a/covariance/covarianceBase.cpp b/covariance/covarianceBase.cpp index ed43b8f42..e937f5547 100644 --- a/covariance/covarianceBase.cpp +++ b/covariance/covarianceBase.cpp @@ -531,57 +531,6 @@ std::vector covarianceBase::getProposed() const { return props; } -// Throw nominal values -void covarianceBase::throwNominal(bool nomValues, int seed) { - TVectorD* vec = new TVectorD(_fNumPar); - for (int i = 0; i < _fNumPar; i++) { - (*vec)(i) = 1.0; - } - - ThrowParms* nom_throws = new ThrowParms(*vec, (*covMatrix)); - nom_throws->SetSeed(seed); - std::vector nominal = getNominalArray(); - nominal.clear(); - nominal.resize(_fNumPar); - - // If we want to put the nominals somewhere else than user specified - // Don't fully understand this though: won't we have to reweight the MC somehow? - // nominal[i] is used in GetLikelihood() as the penalty term, so we're essentially setting a random parameter penalty term? - if (!nomValues) - { - bool throw_again = true; - - while(throw_again == true) - { - throw_again = false; - MACH3LOG_INFO("Setting {} nominal values to random throws.", getName()); - nom_throws->ThrowSet(nominal); - - for (int i = 0; i < _fNumPar; i++) - { - // if parameter is fixed, dont throw - if (_fError[i] < 0) { - nominal[i] = 1.0; - continue; - } - - if (nominal[i] < 0) { - nominal[i] = 0.0; - throw_again = true; - } - } - } - } else { - // If we want nominal values, set all entries to 1 (defined as nominal in MaCh3) - for (int i = 0; i < int(nominal.size()); i++) { - nominal[i] = 1.0; - } - } - - delete nom_throws; - delete vec; -} - // ************************************* // Throw the parameters according to the covariance matrix // This shouldn't be used in MCMC code ase it can break Detailed Balance; @@ -995,8 +944,8 @@ void covarianceBase::setParameters(const std::vector& pars) { unsigned int parsSize = pars.size(); for (unsigned int i = 0; i < parsSize; i++) { - //Make sure that you are actually passing a number to set the parameter to - if(isnan(pars[i])) { + //Make sure that you are actually passing a number to set the parameter to + if(std::isnan(pars[i])) { MACH3LOG_ERROR("Error: trying to set parameter value to a nan for parameter {} in matrix {}. This will not go well!", GetParName(i), matrixName); throw; } else { @@ -1004,7 +953,6 @@ void covarianceBase::setParameters(const std::vector& pars) { } } } - // And if pca make the transfer if (pca) { TransferToPCA(); @@ -1199,8 +1147,7 @@ void covarianceBase::printIndivStepScale() { std::cout << "============================================================" << std::endl; } // ******************************************** -//Makes sure that matrix is positive-definite (so no error is thrown when -//throwNominal() is called) by adding a small number to on-diagonal elements +//Makes sure that matrix is positive-definite by adding a small number to on-diagonal elements void covarianceBase::MakePosDef(TMatrixDSym *cov) { // ******************************************** //DB Save original warning state and then increase it in this function to suppress 'matrix not positive definite' messages diff --git a/covariance/covarianceBase.h b/covariance/covarianceBase.h index 346b9c216..756bb6ef3 100644 --- a/covariance/covarianceBase.h +++ b/covariance/covarianceBase.h @@ -4,11 +4,9 @@ #include "manager/manager.h" #include "samplePDF/Structs.h" #include "covariance/CovarianceUtils.h" -#include "covariance/ThrowParms.h" #include "covariance/AdaptiveMCMCHandler.h" #include "covariance/PCAHandler.h" - #ifndef _LARGE_LOGL_ /// Large Likelihood is used it parameter go out of physical boundary, this indicates in MCMC that such step should eb removed #define _LARGE_LOGL_ 1234567890.0 diff --git a/covariance/covarianceXsec.h b/covariance/covarianceXsec.h index 99e0bea82..849319306 100644 --- a/covariance/covarianceXsec.h +++ b/covariance/covarianceXsec.h @@ -9,9 +9,6 @@ // MaCh3 includes #include "covariance/covarianceBase.h" #include "samplePDF/Structs.h" -#include "manager/YamlHelper.h" - -#include "yaml-cpp/yaml.h" /// @brief Class responsible for handling of systematic error parameters with different types defined in the config. Like spline, normalisation parameters etc. /// @see For more details, visit the [Wiki](https://github.com/mach3-software/MaCh3/wiki/02.-Implementation-of-Systematic). diff --git a/manager/manager.h b/manager/manager.h index 68ccdb53d..53dd472af 100644 --- a/manager/manager.h +++ b/manager/manager.h @@ -1,17 +1,5 @@ #pragma once -// C++ includes -#include -#include -#include -#include -#include - -// ROOT include -#include "TTree.h" -#include "TBranch.h" -#include "TMacro.h" - // MaCh3 Includes #include "manager/MaCh3Logger.h" #include "samplePDF/Structs.h" diff --git a/plotting/plottingUtils/plottingUtils.cpp b/plotting/plottingUtils/plottingUtils.cpp index 93cb2e94c..74f8e27a0 100644 --- a/plotting/plottingUtils/plottingUtils.cpp +++ b/plotting/plottingUtils/plottingUtils.cpp @@ -36,7 +36,7 @@ TH1D TGraphToTH1D(TGraph graph, std::string newName, std::string newTitle) { } // get the bin edges - Double_t binEdges[nPoints + 1]; + std::vector binEdges(nPoints + 1); binEdges[0] = pointsX[0] - (pointsX[1] - pointsX[0]) / 2.0; binEdges[nPoints] = pointsX[nPoints - 1] + (pointsX[nPoints - 1] - pointsX[nPoints - 2]) / 2.0; @@ -46,7 +46,7 @@ TH1D TGraphToTH1D(TGraph graph, std::string newName, std::string newTitle) { binEdges[pointId] = (pointsX[pointId] + pointsX[pointId - 1]) / 2.0; } - TH1D retHist = TH1D(name.c_str(), title.c_str(), nPoints, binEdges); + TH1D retHist = TH1D(name.c_str(), title.c_str(), nPoints, binEdges.data()); for (int binId = 0; binId < nPoints; binId++) { @@ -80,4 +80,4 @@ std::vector> TGraphToVector(TGraph graph) { return ret; } -} // namespace MaCh3Plotting \ No newline at end of file +} // namespace MaCh3Plotting diff --git a/samplePDF/ShiftFunctors.h b/samplePDF/ShiftFunctors.h index 5b28f4c9f..d4b4cb2a9 100644 --- a/samplePDF/ShiftFunctors.h +++ b/samplePDF/ShiftFunctors.h @@ -1,6 +1,6 @@ #pragma once -#include "samplePDFFDBase.h" +#include "samplePDF/Structs.h" class BaseFuncPar{ public: diff --git a/samplePDF/Structs.cpp b/samplePDF/Structs.cpp index 84c90e268..3eac171ae 100644 --- a/samplePDF/Structs.cpp +++ b/samplePDF/Structs.cpp @@ -1,4 +1,4 @@ -#include "Structs.h" +#include "samplePDF/Structs.h" #include "TList.h" #include "TObjArray.h" diff --git a/samplePDF/samplePDFBase.h b/samplePDF/samplePDFBase.h index 3cb43b38f..3030ff469 100644 --- a/samplePDF/samplePDFBase.h +++ b/samplePDF/samplePDFBase.h @@ -1,10 +1,7 @@ #pragma once //C++ includes -#include -#include #include -#include //ROOT includes #include "TTree.h" @@ -13,11 +10,8 @@ #include "TMath.h" #include "TFile.h" #include "TROOT.h" -#include "TRandom.h" -#include "TSpline.h" #include "TRandom3.h" #include "TString.h" -#include "TMath.h" //MaCh3 includes #include "samplePDF/Structs.h" diff --git a/samplePDF/samplePDFFDBase.h b/samplePDF/samplePDFFDBase.h index 4fece720f..fb9e72383 100644 --- a/samplePDF/samplePDFFDBase.h +++ b/samplePDF/samplePDFFDBase.h @@ -1,24 +1,11 @@ #pragma once //C++ includes -#include -#include -#include #include -#include //ROOT includes -#include "TTree.h" -#include "TH1D.h" -#include "TH2D.h" #include "THStack.h" #include "TLegend.h" -#include "TMath.h" -#include "TFile.h" -#include "TGraph2DErrors.h" -#include "TROOT.h" -#include "TRandom.h" -#include "TString.h" //MaCh3 includes #include "OscClass/OscClass_CUDAProb3.h" @@ -32,9 +19,6 @@ #include "samplePDF/FDMCStruct.h" #include "samplePDF/ShiftFunctors.h" - -#define USEBETA 0 - class samplePDFFDBase : virtual public samplePDFBase { public: