Feature flexible binning

CMakeLists.txt

@@ -49,7 +49,7 @@ include(${CMAKE_CURRENT_LIST_DIR}/cmake/Modules/MaCh3Dependencies.cmake)
 
 ############################  C++ Compiler  ####################################
 if (NOT DEFINED CMAKE_CXX_STANDARD OR "${CMAKE_CXX_STANDARD} " STREQUAL " ")
-  SET(CMAKE_CXX_STANDARD 11)
+  SET(CMAKE_CXX_STANDARD 17)
 endif()
 
 # KS: If C++ standard is lower than C++ standard used for ROOT compilation things will go terribly wrong

OscClass/OscClass_CUDAProb3.cpp

@@ -911,11 +911,11 @@ void Oscillator::Reset(int NeutrinoSignIndex, int InitialNeutrinoIndex, int Fina
 
 }
 
-bool Oscillator::isAlreadyCalculated(double* oscpar, double prodH, double Yp_Val) {
+bool Oscillator::isAlreadyCalculated(const double** oscpar, double prodH, double Yp_Val) {
 
   bool fAlreadyCalculated = true;
   for (int i=0;i<nOscpars;i++) {
-    if (oscpar[i]!=foscpar[i]) {
+    if (*oscpar[i]!=foscpar[i]) {
       fAlreadyCalculated = false;
     }
   }
@@ -940,9 +940,9 @@ bool Oscillator::isAlreadyCalculated(double* oscpar, double prodH, double Yp_Val
   return fAlreadyCalculated;
 }
 
-void Oscillator::SaveParams(double* oscpar, double prodH, double Yp_Val) {
+void Oscillator::SaveParams(const double** oscpar, double prodH, double Yp_Val) {
   for (int i=0;i<nOscpars;i++) {
-    foscpar[i] = oscpar[i];
+    foscpar[i] = *oscpar[i];
   }
   fprodH = prodH;
   fYp_Val = Yp_Val;
@@ -952,7 +952,7 @@ void Oscillator::SaveParams(double* oscpar, double prodH, double Yp_Val) {
   fBinning[3] = nFineEnergy;
 }
 
-void  Oscillator::FillOscillogram(double* oscpar, double prodH, double Yp_Val) {
+void  Oscillator::FillOscillogram(const double** oscpar, double prodH, double Yp_Val) {
   if (isAlreadyCalculated(oscpar,prodH,Yp_Val)) {
     return;
   }
@@ -961,8 +961,8 @@ void  Oscillator::FillOscillogram(double* oscpar, double prodH, double Yp_Val) {
 
   //DB oscpar, as given from MaCh3, expresses the mixing angles in sin^2(theta). The propagator expects them in theta
   for (int iOscPar=0;iOscPar<3;iOscPar++) {
-    if (oscpar[iOscPar] < 0) {
-      std::cerr << "Invalid oscillation parameter (Can not sqrt this value)!:" << oscpar[iOscPar] << std::endl;
+    if (*oscpar[iOscPar] < 0) {
+      std::cerr << "Invalid oscillation parameter (Can not sqrt this value)!:" << *oscpar[iOscPar] << std::endl;
       throw;
     }
   }

OscClass/OscClass_CUDAProb3.h

@@ -41,7 +41,7 @@ class Oscillator {
 
   int GetOscillogramNBins(int Switcher);
 
-  void FillOscillogram(double* oscpar, double prodH, double Yp_Val=0.468);
+  void FillOscillogram(const double** oscpar, double prodH, double Yp_Val=0.468);
   double ReturnProb(double NeutrinoEnergy, double Cosz, int InitialFlavour, int FinalFlavour);
 
   void RebinOscillogram(int Switcher, std::vector<double> NewBinning);
@@ -103,8 +103,8 @@ class Oscillator {
   void FillSecondaryHistograms(int NeutrinoSignIndex, int InitialNeutrinoIndex, int FinalNeutrinoIndex);
   void FillPrimaryHistograms(int NeutrinoSignIndex, int InitialNeutrinoIndex, int FinalNeutrinoIndex);
 
-  bool isAlreadyCalculated(double* oscpar, double prodH, double Yp_Val);
-  void SaveParams(double* oscpar, double prodH, double Yp_Val);
+  bool isAlreadyCalculated(const double** oscpar, double prodH, double Yp_Val);
+  void SaveParams(const double** oscpar, double prodH, double Yp_Val);
   void ResetSavedParams();
 
   void PrintBox(Box Box1);

covariance/covarianceBase.cpp

@@ -681,17 +681,15 @@ void covarianceBase::setSingleParameter(const int parNo, const double parVal) {
 // *************************************
   _fPropVal[parNo] = parVal;
   _fCurrVal[parNo] = parVal;
-  MACH3LOG_DEBUG("Setting {}(parameter {}) to {}", GetParName(parNo), parNo, parVal);
-
+  MACH3LOG_DEBUG("Setting {} (parameter {}) to {})" GetParName(parNo),  parNo, parVal);
   if (pca) TransferToPCA();
 }
 // ********************************************
 void covarianceBase::setParCurrProp(const int parNo, const double parVal) {
 // ********************************************
   _fPropVal[parNo] = parVal;
   _fCurrVal[parNo] = parVal;
-  MACH3LOG_DEBUG("Setting {}(parameter {}) to {}", GetParName(parNo), parNo, parVal);
-
+  MACH3LOG_DEBUG("Setting {} (parameter {}) to {})" GetParName(parNo),  parNo, parVal);
   if (pca) TransferToPCA();
 }
 
@@ -997,16 +995,22 @@ void covarianceBase::setParameters(const std::vector<double>& 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(std::isnan(pars[i])) {
-        std::cerr << "Error: trying to set parameter value to a nan for parameter " << GetParName(i) << " in matrix " << matrixName << ". This will not go well!" << std::endl;
+	  //Make sure that you are actually passing a number to set the parameter to
+	  if(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 {
 		_fPropVal[i] = pars[i];
 	  }
     }
   }
 
+  MACH3LOG_DEBUG("Set parameters to:");
+  std::cout << "Set parameters to: " << std::endl;
+  for(int par_i = 0 ; par_i < pars.size() ; ++par_i){
+	MACH3LOG_DEBUG("{} : {}", par_i, _fPropVal[par_i]);
+  }
+
   // And if pca make the transfer
   if (pca) {
     TransferToPCA();
@@ -1029,6 +1033,7 @@ void covarianceBase::SetBranches(TTree &tree, bool SaveProposal) {
       tree.Branch(Form("%s_PCA", _fNames[i].c_str()), (double*)&(fParCurr_PCA.GetMatrixArray()[i]), Form("%s_PCA/D", _fNames[i].c_str()));
     }
   }
+
   if(SaveProposal)
   {
     // loop over parameters and set a branch
@@ -1043,6 +1048,7 @@ void covarianceBase::SetBranches(TTree &tree, bool SaveProposal) {
     }
   }
 }
+
 // ********************************************
 void covarianceBase::setStepScale(const double scale) {
 // ********************************************

covariance/covarianceBase.h

@@ -147,6 +147,8 @@ class covarianceBase {
   /// @brief Get diagonal error for ith parameter
   /// @param i Parameter index
   inline double getDiagonalError(const int i) { return std::sqrt((*covMatrix)(i,i)); }
+  // Get the error for the ith parameter
+  const double GetError(const int i) {return _fError[i];}
 
   /// @brief Adaptive Step Tuning Stuff
   void resetIndivStepScale();
@@ -175,19 +177,6 @@ class covarianceBase {
   /// @brief KS: Convert covariance matrix to correlation matrix and return TH2D which can be used for fancy plotting
   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];
-    } else {
-      MACH3LOG_WARN("Specified bin is <= 0 OR bin > npar!");
-      MACH3LOG_WARN("bin = {}, npar = {}", bin, _fNumPar);
-      MACH3LOG_WARN("This won't ruin much that this step in the MCMC, but does indicate something wrong in memory!");
-      return 1.0;
-    }
-    return 1.0;
-  }
   //========
   //ETA - This might be a bit squiffy? If the vector gots moved from say a
   //push_back then the pointer is no longer valid... maybe need a better