Feature core 1.2.0

covariance/covarianceXsec.cpp

@@ -1,4 +1,5 @@
 #include "covariance/covarianceXsec.h"
+#include "samplePDF/Structs.h"
 
 // ********************************************
 // ETA - YAML constructor
@@ -114,10 +115,29 @@ const std::vector<std::string> covarianceXsec::GetSplineParsNamesFromDetID(const
   for (auto &pair : _fSystToGlobalSystIndexMap[SystType::kSpline]) {
     auto &SplineIndex = pair.first;
     auto &SystIndex = pair.second;
-    if ((GetParDetID(SystIndex) & DetID )){
+    if (AppliesToDetID(SystIndex, DetID)) { //If parameter applies to required DetID
       returnVec.push_back(_fSplineNames.at(SplineIndex));
     }
+
+  }
+  return returnVec;
+}
+
+const std::vector<SplineInterpolation> covarianceXsec::GetSplineInterpolationFromDetID(const int DetID) {
+  std::vector<SplineInterpolation> returnVec;
+  for (auto &pair : _fSystToGlobalSystIndexMap[SystType::kSpline]) {
+    auto &SplineIndex = pair.first;
+    auto &SystIndex = pair.second;
+
+    if (AppliesToDetID(SystIndex, DetID)) { //If parameter applies to required DetID
+      returnVec.push_back(SplineParams.at(SplineIndex)._SplineInterpolationType);
+    }
+
+    // if ((GetParDetID(SystIndex) & DetID )){
+    //   returnVec.push_back(SplineParams.at(SplineIndex)._SplineInterpolationType);
+    // }
   }
+
   return returnVec;
 }
 

covariance/covarianceXsec.h

@@ -32,6 +32,8 @@ class covarianceXsec : public covarianceBase {
     /// @brief Get interpolation type for a given parameter
     /// @param i spline parameter index, not confuse with global index
     inline SplineInterpolation GetParSplineInterpolation(const int i) {return SplineParams.at(i)._SplineInterpolationType;}
+    /// @brief Get the interpolation types for splines affecting a particular DetID
+    const std::vector<SplineInterpolation> GetSplineInterpolationFromDetID(int DetID);
     /// @brief Get the name of the spline associated with the spline at index i
     /// @param i spline parameter index, not to be confused with global index
     std::string GetParSplineName(const int i) {return _fSplineNames[i];}

mcmc/MaCh3Factory.h

@@ -124,9 +124,7 @@ std::vector<SampleType*> MaCh3SamplePDFFactory(const std::vector<std::string>& S
   for (size_t i = 0; i < SampleConfig.size(); ++i)
   {
     // Instantiate the sample using the specified class type
-    SampleType* Sample = new SampleType(SampleConfig[i], xsec);
-    Sample->SetXsecCov(xsec);
-    if (osc != nullptr) Sample->SetOscCov(osc);
+    SampleType* Sample = new SampleType(SampleConfig[i], xsec, osc);
     Sample->reweight();
 
     // Obtain sample name and create a TString version for histogram naming

python/samplePDF.cpp

@@ -274,19 +274,6 @@ void initSamplePDF(py::module &m){
             py::arg("mc_version"),
             py::arg("xsec_cov")
         )
-
-        .def(
-            "set_xsec_cov", 
-            &samplePDFFDBase::SetXsecCov, 
-            "Set the cross section covariance matrix object."
-        )
-
-        .def(
-            "set_osc_cov", 
-            &samplePDFFDBase::SetOscCov, 
-            "Set the oscillation parameter covariance matrix object."
-        )
-
     ;
 
     /* Not sure if this will be needed in future versions of MaCh3 so leaving commented for now

samplePDF/FarDetectorCoreInfoStruct.h

@@ -10,16 +10,14 @@ struct FarDetectorCoreInfo {
 
   ~FarDetectorCoreInfo(){ delete [] isNC; }
 
-  int nutype; // 2 = numu/signue | -2 = numub | 1 = nue | -1 = nueb           
-  int oscnutype;    
-  int nupdg;
-  int nupdgUnosc;
   bool signal; // true if signue                                              
   int nEvents; // how many MC events are there
   double ChannelIndex;
   std::string flavourName;
 
   std::vector<int*> Target; // target the interaction was on
+  std::vector<const int*> nupdg;
+  std::vector<const int*> nupdgUnosc;
 
   int SampleDetID;
 

samplePDF/Structs.cpp

@@ -1,7 +1,13 @@
 #include "samplePDF/Structs.h"
 
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wfloat-conversion"
+#include "Constants/OscillatorConstants.h"
+#pragma GCC diagnostic pop
+
 #include "TList.h"
 #include "TObjArray.h"
+#include "manager/MaCh3Logger.h"
 
 namespace MaCh3Utils {
 
@@ -71,6 +77,31 @@ namespace MaCh3Utils {
     return 0;
   }
 
+  int PDGToNuOscillatorFlavour(int NuPdg){
+    int NuOscillatorFlavour = _BAD_INT_;
+    switch(std::abs(NuPdg)){
+      case NuPDG::kNue:
+        NuOscillatorFlavour = NuOscillator::kElectron;
+        break;
+      case NuPDG::kNumu:
+        NuOscillatorFlavour = NuOscillator::kMuon;
+        break;
+      case NuPDG::kNutau:
+        NuOscillatorFlavour = NuOscillator::kTau;
+        break;
+      default:
+        MACH3LOG_ERROR("Unknown Nuetrino PDG {}, cannot convert to NuOscillator type", NuPdg);
+        break;
+    }
+
+    //This is very cheeky but if the PDG is negative then multiply the PDG by -1
+    // This is consistent with the treatment that NuOscillator expects as enums only
+    // exist for the generic matter flavour and not the anti-matter version
+    if(NuPdg < 0){NuOscillatorFlavour *= -1;}
+
+    return NuOscillatorFlavour;
+  }
+
 
   //DB Anything added here must be of the form 2^X, where X is an integer
   //

samplePDF/Structs.h

@@ -270,6 +270,13 @@ namespace MaCh3Utils {
   /// @brief Return mass for given PDG
   double GetMassFromPDG(int PDG);
   // ***************************
+
+  // ***************************
+  /// @brief Convert from PDG flavour to NuOscillator type
+  /// beware that in the case of anti-neutrinos the NuOscillator
+  /// type simply gets multiplied by -1
+  int PDGToNuOscillatorFlavour(int PDG);
+  // ***************************
 
   extern std::unordered_map<int,int> KnownDetIDsMap;
   extern int nKnownDetIDs;