diff --git a/CMakeLists.txt b/CMakeLists.txt index 0dfc54ab..cf8c8403 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -78,15 +78,15 @@ set( CMAKE_EXPORT_COMPILE_COMMANDS ON ) #Compile and link options #KS: Add -Werror and maybe -Wconversion in the future target_compile_options(MaCh3CompilerOptions INTERFACE - -g - -Wextra - -Wall - -pedantic - -Wshadow - -Wuninitialized - -Wnon-virtual-dtor - -Woverloaded-virtual - ) + -g # Generate debug information + -Wextra # Enable extra warning flags + -Wall # Enable all standard warning flags + -pedantic # Enforce strict ISO compliance + -Wshadow # Warn when a variable declaration shadows one from an outer scope + -Wuninitialized # Warn about uninitialized variables + -Wnon-virtual-dtor # Warn when a class with virtual functions has a non-virtual destructor + -Woverloaded-virtual # Warn when a function declaration hides a virtual function from a base class +) #KS: If Debug is not defined disable it by default DefineEnabledRequiredSwitch(MaCh3_DEBUG_ENABLED FALSE) @@ -94,7 +94,7 @@ DefineEnabledRequiredSwitch(MaCh3_DEBUG_ENABLED FALSE) if(DEFINED DEBUG_LEVEL) SET(MaCh3_DEBUG_ENABLED TRUE) else() - #If MaCh debug was enable but level not, set it to 1. In very rare cases we want to go beyond 1. + #If MaCh3 debug was enable but level not, set it to 1. In very rare cases we want to go beyond 1. if(MaCh3_DEBUG_ENABLED) SET(DEBUG_LEVEL 1) endif() @@ -107,14 +107,14 @@ if(MaCh3_DEBUG_ENABLED) else() #KS: Consider in future __attribute__((always_inline)) see https://indico.cern.ch/event/386232/sessions/159923/attachments/771039/1057534/always_inline_performance.pdf #https://gcc.gnu.org/onlinedocs/gcc-3.3.6/gcc/Optimize-Options.html - target_compile_options(MaCh3CompilerOptions INTERFACE - -O3 - -funroll-loops - --param=max-vartrack-size=100000000 - -finline-limit=100000000 - #-flto # FIXME need more testing - ) - #KS: addomg Link-Time Optimization (LTO) +target_compile_options(MaCh3CompilerOptions INTERFACE + -O3 # Optimize code for maximum speed + -funroll-loops # Unroll loops where possible for performance + --param=max-vartrack-size=100000000 # Set maximum size of variable tracking data to avoid excessive memory usage + -finline-limit=100000000 # Increase the limit for inlining functions to improve performance + #-flto # FIXME need more testing # Enable link-time optimization (commented out for now, needs more testing) +) + #KS: add Link-Time Optimization (LTO) # FIXME previously it wasn't used correctly but would need more testing #target_link_libraries(MaCh3CompilerOptions INTERFACE -flto) endif() diff --git a/covariance/ThrowParms.h b/covariance/ThrowParms.h index a6e6d323..92e8969b 100755 --- a/covariance/ThrowParms.h +++ b/covariance/ThrowParms.h @@ -14,6 +14,8 @@ #include #include +/// @brief No idea really... +/// @warning this is deprecated don't use it class ThrowParms { private: TVectorD *pvals; diff --git a/covariance/covarianceBase.h b/covariance/covarianceBase.h index f1363105..aa131137 100644 --- a/covariance/covarianceBase.h +++ b/covariance/covarianceBase.h @@ -50,6 +50,7 @@ class covarianceBase { /// @brief Constructor For Eigen Value decomp /// @param name Matrix name /// @param file Path to matrix root file + /// @param seed Seed for TRandom3 /// @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. @@ -66,15 +67,19 @@ class covarianceBase { void setName(const char *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 setSingleParameter(const int parNo, const double parVal); /// @brief Set all the covariance matrix parameters to a user-defined value /// @param i Parameter index + /// @param val new value which will be set void setPar(const int i, const double val); /// @brief Set current parameter value /// @param i Parameter index + /// @param val new value which will be set void setParCurrProp(const int i, const double val); /// @brief Set proposed parameter value + /// @param val new value which will be set void setParProp(const int i, const double val) { _fPropVal[i] = val; if (pca) TransferToPCA(); @@ -88,6 +93,8 @@ class covarianceBase { void setFlatPrior(const int i, const bool eL); /// @brief set branches for output file + /// @param tree Tree to which we will save branches + /// @param SaveProposal Normally we only save parameter after is accepted, for debugging purpose it is helpful to see also proposed values. That's what this variable controls void SetBranches(TTree &tree, bool SaveProposal = false); /// @brief Set global step scale for covariance object /// @param scale Value of global step scale @@ -175,6 +182,7 @@ class covarianceBase { TH2D* GetCorrelationMatrix(); /// @brief What parameter Gets reweighted by what amount according to MCMC + /// @param bin simply parameter index inline double calcReWeight(const int bin) { if (bin >= 0 && bin < _fNumPar) { return _fPropVal[bin]; @@ -201,6 +209,7 @@ class covarianceBase { virtual std::vector getNominalArray(); std::vector getPreFitValues(){return _fPreFitValue;} std::vector getGeneratedValues(){return _fGenerated;} + /// @brief Get vector of all proposed parameter values std::vector getProposed() const; /// @brief Get proposed parameter value /// @param i Parameter index @@ -211,8 +220,11 @@ class covarianceBase { /// @brief Get prior parameter value /// @param i Parameter index inline double getParInit(const int i) { return _fPreFitValue[i]; } - + /// @brief Return generated value, although is virtual so class inheriting might actual get nominal not generated. + /// @param i Parameter index virtual double getNominal(const int i) { return getParInit(i); } + /// @brief Return generated value for a given parameter + /// @param i Parameter index inline double GetGenerated(const int i) { return _fGenerated[i];} /// @brief Get upper parameter bound in which it is physically valid /// @param i Parameter index @@ -243,34 +255,38 @@ class covarianceBase { if (fParSigma_PCA[i] < 0) { return true; } else { return false; } } - + /// @brief Get transfer matrix allowing to go from PCA base to normal base inline const TMatrixD getTransferMatrix() { if (!pca) { MACH3LOG_ERROR("Am not running in PCA mode"); throw; } return TransferMat; } - + /// @brief Get eigen vectors of covariance matrix, only works with PCA inline const TMatrixD getEigenVectors() { if (!pca) { MACH3LOG_ERROR("Am not running in PCA mode"); throw; } return eigen_vectors; } - + /// @brief Get eigen values for all parameters, if you want for decomposed only parameters use getEigenValuesMaster inline const TVectorD getEigenValues() { if (!pca) { MACH3LOG_ERROR("Am not running in PCA mode"); throw; } return eigen_values; } - + /// @brief Get eigen value of only decomposed parameters, if you want for all parameters use getEigenValues inline const std::vector getEigenValuesMaster() { if (!pca) { MACH3LOG_ERROR("Am not running in PCA mode"); throw; } return eigen_values_master; } - + /// @brief Set proposed value for parameter in PCA base + /// @param i Parameter index + /// @param value new value inline void setParProp_PCA(const int i, const double value) { if (!pca) { MACH3LOG_ERROR("Am not running in PCA mode"); throw; } fParProp_PCA(i) = value; // And then transfer back to the parameter basis TransferToParam(); } - + /// @brief Set current value for parameter in PCA base + /// @param i Parameter index + /// @param value new value inline void setParCurr_PCA(const int i, const double value) { if (!pca) { MACH3LOG_ERROR("Am not running in PCA mode"); throw; } fParCurr_PCA(i) = value; diff --git a/manager/MaCh3Logger.h b/manager/MaCh3Logger.h index c0836dc9..4082288d 100644 --- a/manager/MaCh3Logger.h +++ b/manager/MaCh3Logger.h @@ -56,15 +56,20 @@ inline void SetMaCh3LoggerFormat() template void LoggerPrint(const std::string& LibName, LogFunc logFunction, Func&& func, Args&&... args) { + // Create a stringstream to capture the output std::stringstream sss; + // Save the original stream buffers std::streambuf* coutBuf = std::cout.rdbuf(); std::streambuf* cerrBuf = std::cerr.rdbuf(); + // Redirect std::cout and std::cerr to the stringstream buffer std::cout.rdbuf(sss.rdbuf()); std::cerr.rdbuf(sss.rdbuf()); + // Call the provided function func(std::forward(args)...); + // Restore the original stream buffers std::cout.rdbuf(coutBuf); std::cerr.rdbuf(cerrBuf); diff --git a/mcmc/PSO.h b/mcmc/PSO.h index c80a0a78..782b0976 100644 --- a/mcmc/PSO.h +++ b/mcmc/PSO.h @@ -2,7 +2,6 @@ // // Created by Emily Ip on 24/2/2023. // -// // Created by Emily Ip on 26/1/2023. // #include @@ -16,9 +15,7 @@ /// @brief Class particle - stores the position, velocity and personal best /// With functions which move particle and update velocity class particle{ - public: - particle(){}; particle(std::vector pos, std::vector vel) : position(pos), velocity(vel){}; /// @brief Destructor @@ -69,16 +66,12 @@ class particle{ double personal_best_value; double curr_value; std::vector personal_best_position; - - }; - /// @brief Class PSO, consist of a vector of object Class Particle and global best /// Takes in the size (number of particle) and number of iteration /// functions includes: finding global best, updating velocity, actual minimisation function class PSO : public LikelihoodFit { - public: /// @brief constructor PSO(manager * const fitMan); diff --git a/samplePDF/FDMCStruct.h b/samplePDF/FDMCStruct.h index a4cea5b5..e951f3ba 100644 --- a/samplePDF/FDMCStruct.h +++ b/samplePDF/FDMCStruct.h @@ -1,7 +1,7 @@ #pragma once -//constructors are same for all three so put in here +/// @brief constructors are same for all three so put in here struct fdmc_base { int nutype; // 2 = numu/signue | -2 = numub | 1 = nue | -1 = nueb int oscnutype; @@ -20,10 +20,10 @@ struct fdmc_base { double** y_var; double **rw_etru; - // xsec bins + /// xsec bins std::list< int > *xsec_norms_bins; - //DB Speedup bits + /// DB Speedup bits double Unity; int* nxsec_norm_pointers; @@ -54,7 +54,7 @@ struct fdmc_base { int **mode; - // DB Atmospheric Parameters + /// DB Atmospheric Parameters const double **osc_w_pointer; double *rw_truecz; diff --git a/samplePDF/Structs.h b/samplePDF/Structs.h index 90eef5f4..3b778cb2 100644 --- a/samplePDF/Structs.h +++ b/samplePDF/Structs.h @@ -75,7 +75,7 @@ constexpr unsigned int str2int(const char* str, int h = 0) { } // ******************* -/// ETA - Normalisations for cross-section parameters +/// @brief ETA - Normalisations for cross-section parameters /// Carrier for whether you want to apply a systematic to an event or not struct XsecNorms4 { // ******************* @@ -97,11 +97,11 @@ struct XsecNorms4 { /// and lower and upper bounds. This can then be passed to IsEventSelected std::vector< std::vector > Selection; - /// Generic vector containing the string of kinematic type - /// This then needs to be converted to a kinematic type enum - /// within a samplePDF daughter class - /// The bounds for each kinematic variable are given in Selection - std::vector< std::string > KinematicVarStr; + /// Generic vector containing the string of kinematic type + /// This then needs to be converted to a kinematic type enum + /// within a samplePDF daughter class + /// The bounds for each kinematic variable are given in Selection + std::vector< std::string > KinematicVarStr; /// Parameter number of this normalisation in current systematic model int index; @@ -117,7 +117,7 @@ enum SplineInterpolation { }; // ************************************************** -/// Convert a LLH type to a string +/// @brief Convert a LLH type to a string inline std::string SplineInterpolation_ToString(const SplineInterpolation i) { // ************************************************** std::string name = ""; @@ -154,7 +154,7 @@ enum SystType { }; // ************************************************** -/// Convert a Syst type type to a string +/// @brief Convert a Syst type type to a string inline std::string SystType_ToString(const SystType i) { // ************************************************** std::string name = ""; @@ -184,7 +184,7 @@ namespace MaCh3Utils { // *************************** // *************************** - // Return mass for given PDG + /// @brief Return mass for given PDG double GetMassFromPDG(int PDG); // *************************** diff --git a/samplePDF/samplePDFBase.cpp b/samplePDF/samplePDFBase.cpp index dad163c1..831aefd2 100644 --- a/samplePDF/samplePDFBase.cpp +++ b/samplePDF/samplePDFBase.cpp @@ -273,7 +273,7 @@ double samplePDFBase::getTestStatLLH(const double data, const double mc, const d return stat+penalty; } break; - //KS: Alterantive calcaution of Barlow-Beeston following Hans Dembinski and Ahmed Abdelmottele arXiv:2206.12346v2 + //KS: Alternative calculation of Barlow-Beeston following Hans Dembinski and Ahmed Abdelmottele arXiv:2206.12346v2 case (kDembinskiAbdelmottele): { //KS: code follows authors implementation from: @@ -334,7 +334,7 @@ double samplePDFBase::getTestStatLLH(const double data, const double mc, const d //KS: Pearson works on assumption that event distribution in each bin is described by a Gaussian which in our case is not fulfilled for all bins, hence use it at your own risk case (kPearson): { - //KS: 2 is beacuese this function returns -LLH not -2LLH + //KS: 2 is because this function returns -LLH not -2LLH const double stat = (data-mc)*(data-mc)/(2*mc); // Return the statistical diff --git a/samplePDF/samplePDFBase.h b/samplePDF/samplePDFBase.h index 8dfe8f0d..0fb0372f 100644 --- a/samplePDF/samplePDFBase.h +++ b/samplePDF/samplePDFBase.h @@ -23,12 +23,13 @@ #include "samplePDF/Structs.h" #include "manager/manager.h" +/// @brief Class responsible for handling implementation of samples used in analysis, reweighting and returning LLH class samplePDFBase { public: samplePDFBase(){}; samplePDFBase(double pot); - + /// @brief destructor virtual ~samplePDFBase(); virtual inline _int_ GetNsamples(){ return nSamples; }; @@ -82,6 +83,10 @@ class samplePDFBase (void) sample; (void) Dimension; throw MaCh3Exception(__FILE__ , __LINE__ , "Not implemented"); }; double getTestStatLLH(double data, double mc); + /// @brief Calculate test statistic for a single bin. Calculation depends on setting of fTestStatistic + /// @param data is data + /// @param mc is mc + /// @param w2 is is Sum(w_{i}^2) (sum of weights squared), which is sigma^2_{MC stats} double getTestStatLLH(const double data, const double mc, const double w2); // Provide a setter for the test-statistic //void SetTestStatistic(TestStatistic test_stat); @@ -111,11 +116,11 @@ class samplePDFBase std::vector* dataSample; std::vector< std::vector >* dataSample2D; - // Contains how many samples we've got + /// Contains how many samples we've got _int_ nSamples; - //KS: number of dimension for this sample + /// KS: number of dimension for this sample int nDims; - //Name of Sample + /// Name of Sample std::vector SampleName; //GetterForModes diff --git a/samplePDF/samplePDFFDBase.h b/samplePDF/samplePDFFDBase.h index 503fc89a..2c7914e3 100644 --- a/samplePDF/samplePDFFDBase.h +++ b/samplePDF/samplePDFFDBase.h @@ -55,7 +55,7 @@ class samplePDFFDBase : virtual public samplePDFBase void addData(TH2D* Data); void addData(std::vector &data); void addData(std::vector< std::vector > &data); - //DB Multi-threaded GetLikelihood + /// DB Multi-threaded GetLikelihood double GetLikelihood(); //=============================================================================== diff --git a/splines/SplineStructs.h b/splines/SplineStructs.h index 5001afdd..a4c20b3e 100644 --- a/splines/SplineStructs.h +++ b/splines/SplineStructs.h @@ -5,7 +5,7 @@ #include "samplePDF/Structs.h" // ******************* -/// CW: Add a struct to hold info about the splinified xsec parameters and help with FindSplineSegment +/// @brief CW: Add a struct to hold info about the splinified xsec parameters and help with FindSplineSegment struct FastSplineInfo { // ******************* /// Number of points in spline @@ -22,7 +22,7 @@ struct FastSplineInfo { }; // ******************************************** -/// CW: Generic xsec class. Can use TF1 or TSpline3 or TSpline5 here, tjoho +/// @brief CW: Generic xsec class. Can use TF1 or TSpline3 or TSpline5 here, tjoho template class XSecStruct{ // ******************************************** @@ -140,12 +140,11 @@ inline void ApplyKnotWeightCap(TGraph* xsecgraph, int splineParsIndex, covarianc } // ************************ -/// CW: A reduced TF1 class only -/// Only saves parameters for each TF1 and how many parameters each parameter set has +/// @brief CW: A reduced TF1 class only. Only saves parameters for each TF1 and how many parameters each parameter set has class TF1_red { // ************************ public: - /// Empty constructor + /// @brief Empty constructor TF1_red() { length = 0; Par = NULL; @@ -194,7 +193,7 @@ class TF1_red { Func = NULL; } - /// Evaluate a variation + /// @brief Evaluate a variation inline double Eval(_float_ var) { /// If we have 5 parameters we're using a fifth order polynomial if (length == 5) { @@ -211,12 +210,12 @@ class TF1_red { } } - /// Set a parameter to a value + /// @brief Set a parameter to a value inline void SetParameter(_int_ Parameter, _float_ Value) { Par[Parameter] = Value; } - /// Get a parameter value + /// @brief Get a parameter value double GetParameter(_int_ Parameter) { if (Parameter > length) { std::cerr << "Error: you requested parameter number " << Parameter << " but length is " << length << " parameters" << std::endl; @@ -226,13 +225,14 @@ class TF1_red { return Par[Parameter]; } - /// Set the size + /// @brief Set the size inline void SetSize(_int_ nSpline) { length = nSpline; Par = new _float_[length]; } - /// Get the size + /// @brief Get the size inline int GetSize() { return length; } + /// @brief Print detailed info inline void Print() { std::cout << "Printing TF1_red: " << std::endl; std::cout << " ParamNo = " << ParamNo << std::endl; @@ -577,7 +577,7 @@ class TSpline3_red { XPos = YResp = NULL; } - /// Find the segment relevant to this variation in x + /// @brief Find the segment relevant to this variation in x /// See root/hist/hist/src/TSpline3::FindX(double) or samplePDFND....::FindSplineSegment inline int FindX(double x) { // The segment we're interested in (klow in ROOT code) @@ -642,7 +642,7 @@ class TSpline3_red { y = YResp[segment]; } - /// CW: Get the number + /// @brief CW: Get the number inline std::string GetName() { std::stringstream ss; ss << ParamNo; @@ -668,24 +668,24 @@ class Truncated_Spline: public TSpline3_red { // ************************ // cubic spline which is flat (returns y_first or y_last) if x outside of knot range public: - /// Empty constructor + /// @brief Empty constructor Truncated_Spline() :TSpline3_red() { } - /// The constructor that takes a TSpline3 pointer and copies in to memory + /// @brief The constructor that takes a TSpline3 pointer and copies in to memory Truncated_Spline(TSpline3* &spline, int Param = -1) :TSpline3_red(spline, Param) { } - /// Empty destructor + /// @brief Empty destructor ~Truncated_Spline() { } - /// Find the segment relevant to this variation in x + /// @brief Find the segment relevant to this variation in x /// See root/hist/hist/src/TSpline3::FindX(double) or samplePDFND....::FindSplineSegment inline int FindX(double x) { // The segment we're interested in (klow in ROOT code) @@ -719,7 +719,7 @@ class Truncated_Spline: public TSpline3_red { return segment; } - /// Evaluate the weight from a variation + /// @brief Evaluate the weight from a variation inline double Eval(double var) { // Get the segment for this variation int segment = FindX(var);