Implement NuFAST

Apps/DragRace.cpp

@@ -54,6 +54,11 @@ int main() {
   ConfigNames.push_back("./Configs/Unbinned_CUDAProb3.yaml");
 #endif
 
+#if UseNuFASTLinear == 1
+  ConfigNames.push_back("./Configs/Binned_NuFASTLinear.yaml");
+  ConfigNames.push_back("./Configs/Unbinned_NuFASTLinear.yaml");
+#endif
+
 #if UseCUDAProb3Linear == 1
   ConfigNames.push_back("./Configs/Binned_CUDAProb3Linear.yaml");
   ConfigNames.push_back("./Configs/Unbinned_CUDAProb3Linear.yaml");
@@ -114,7 +119,7 @@ int main() {
     auto t1 = high_resolution_clock::now();
     for (int iThrow=0;iThrow<nThrows;iThrow++) {      
       //Throw dcp to some new value
-      FLOAT_T RandVal = rand();
+      FLOAT_T RandVal = static_cast <FLOAT_T> (rand()) / static_cast <FLOAT_T> (RAND_MAX);
       OscParams_Atm[5] = RandVal;
       OscParams_Beam[5] = RandVal;
 

Apps/SingleOscProbCalcer.cpp

@@ -17,6 +17,10 @@
 #include "OscProbCalcer/OscProbCalcer_ProbGPULinear.h"
 #endif
 
+#if UseNuFASTLinear==1
+#include "OscProbCalcer/OscProbCalcer_NuFASTLinear.h"
+#endif
+
 #include <iostream>
 #include <math.h>
 #include <chrono>
@@ -109,6 +113,16 @@ int main(int argc, char **argv) {
     throw;
 #endif
   }
+
+  else if (OscProbCalcerImplementationToCreate == "NuFASTLinear") {
+#if UseNuFASTLinear==1
+    OscProbCalcerNuFASTLinear* NuFASTLinear = new OscProbCalcerNuFASTLinear(OscProbCalcerConfigname);
+    Calcer = (OscProbCalcerBase*)NuFASTLinear;
+#else
+    std::cerr << "Oscillator was requsted to create " << OscProbCalcerImplementationToCreate << " OscProbCalcer but Use" << OscProbCalcerImplementationToCreate << " is undefined. Indicates problem in setup" << std::endl;
+    throw;
+#endif
+  }
 
   else {
     std::cerr << "Executable was requsted to create " << OscProbCalcerImplementationToCreate << " OscProbCalcer but this is not implemented within " << __FILE__ << std::endl;

CMakeLists.txt

@@ -38,7 +38,7 @@ endif()
 add_library(NuOscillatorCompilerOptions INTERFACE)
 
 message(STATUS "NuOscillator Features: ")
-foreach(VAR UseGPU UseMultithreading UseDoubles UseCUDAProb3 UseCUDAProb3Linear UseProb3ppLinear UseProbGPULinear)
+foreach(VAR UseGPU UseMultithreading UseDoubles UseCUDAProb3 UseCUDAProb3Linear UseProb3ppLinear UseProbGPULinear UseNuFASTLinear)
 	    if(NOT DEFINED ${VAR})
 	    	   message(FATAL_ERROR "${VAR} not defined")
 	    endif()
@@ -204,7 +204,7 @@ if(${UseGPU} EQUAL 1)
 endif()
 
 if(${UseProb3ppLinear} EQUAL 1)
-    CPMFindPackage(
+    CPMAddPackage(
         NAME Prob3plusplus
         VERSION 3.10.4
         GITHUB_REPOSITORY rogerwendell/Prob3plusplus
@@ -221,6 +221,16 @@ if(${UseProb3ppLinear} EQUAL 1)
         LIBRARY DESTINATION lib/)
 endif()
 
+if(${UseNuFASTLinear} EQUAL 1)
+    CPMAddPackage(
+      	NAME NuFAST
+      	GIT_TAG main
+      	GITHUB_REPOSITORY PeterDenton/NuFast
+      	DOWNLOAD_ONLY
+      )
+      target_compile_definitions(NuOscillatorCompilerOptions INTERFACE UseNuFASTLinear=1)
+endif()
+
 if(${UseGPU} EQUAL 1)
     target_compile_definitions(NuOscillatorCompilerOptions INTERFACE UseGPU=1)
 else()

Configs/Binned_NuFASTLinear.yaml

@@ -0,0 +1,12 @@
+General:
+  Verbosity: "NONE"
+  CosineZIgnored: true
+  CalculationType: "Binned"
+
+OscProbCalcer:
+  - Config: "NuFASTLinear:./build/Linux/Configs/OscProbCalcerConfigs/NuFASTLinear.yaml"
+
+Binned:
+  FileName: "./Inputs/ExampleAtmosphericBinning.root"
+  EnergyAxisHistName: "EnergyAxisBinning"
+  CosineZAxisHistName: "CosineZAxisBinning"

Configs/OscProbCalcerConfigs/NuFASTLinear.yaml

@@ -0,0 +1,15 @@
+General:
+  Verbosity: "NONE"
+
+OscProbCalcerSetup:
+  - NuFASTLinear:
+    OscChannelMapping:
+      - Entry: "Electron:Electron"
+      - Entry: "Electron:Muon"
+      - Entry: "Electron:Tau"
+      - Entry: "Muon:Electron"
+      - Entry: "Muon:Muon"
+      - Entry: "Muon:Tau"
+      - Entry: "Tau:Electron"
+      - Entry: "Tau:Muon"
+      - Entry: "Tau:Tau"

Configs/Unbinned_NuFASTLinear.yaml

@@ -0,0 +1,7 @@
+General:
+  Verbosity: "NONE"
+  CosineZIgnored: true
+  CalculationType: "Unbinned"
+
+OscProbCalcer:
+  - Config: "NuFASTLinear:./build/Linux/Configs/OscProbCalcerConfigs/NuFASTLinear.yaml"

OscProbCalcer/CMakeLists.txt

@@ -51,6 +51,12 @@ if(${UseProb3ppLinear} EQUAL 1)
 	   list(APPEND HEADERS OscProbCalcer_Prob3ppLinear.h)
 endif()
 
+if(${UseNuFASTLinear} EQUAL 1)
+	   include_directories(${NuFAST_SOURCE_DIR})
+	   target_sources(OscProbCalcer PRIVATE OscProbCalcer_NuFASTLinear.cpp)
+	   list(APPEND HEADERS OscProbCalcer_NuFASTLinear.h)
+endif()
+
 if (${UseGPU} EQUAL 1)
    set_target_properties(OscProbCalcer PROPERTIES
 				    CUDA_SEPARABLE_COMPILATION ON

OscProbCalcer/OscProbCalcer_NuFASTLinear.cpp

@@ -0,0 +1,106 @@
+#include "OscProbCalcer_NuFASTLinear.h"
+
+#include <iostream>
+
+#include "c++/NuFast.cpp"
+
+OscProbCalcerNuFASTLinear::OscProbCalcerNuFASTLinear(std::string ConfigName_, int Instance_) : OscProbCalcerBase(ConfigName_,"NuFASTLinear",Instance_)
+{
+  //=======
+  //Grab information from the config
+
+  //=======
+
+  fNOscParams = kNOscParams;
+
+  fNNeutrinoTypes = 2;
+  InitialiseNeutrinoTypesArray(fNNeutrinoTypes);
+  fNeutrinoTypes[0] = Nu;
+  fNeutrinoTypes[1] = Nubar;
+
+  // This implementation only considers linear propagation, thus no requirement to set cosineZ array
+  IgnoreCosineZBinning(true);
+}
+
+void OscProbCalcerNuFASTLinear::SetupPropagator() {
+}
+
+void OscProbCalcerNuFASTLinear::CalculateProbabilities(std::vector<FLOAT_T> OscParams) {
+  double L, E, rho, Ye, probs_returned[3][3];
+  double s12sq, s13sq, s23sq, delta, Dmsq21, Dmsq31;
+  int N_Newton;
+
+  // ------------------------------- //
+  // Set the experimental parameters //
+  // ------------------------------- //
+  L = OscParams[kPATHL]; // km
+  rho = OscParams[kDENS]; // g/cc
+
+  //Ye = OscParams[kELECDENS];
+  Ye = 0.5;
+
+  // --------------------------------------------------------------------- //
+  // Set the number of Newton-Raphson iterations which sets the precision. //
+  // 0 is close to the single precision limit and is better than DUNE/HK   //
+  // in the high statistics regime. Increasig N_Newton to 1,2,... rapidly  //
+  // improves the precision at a modest computational cost                 //
+  // --------------------------------------------------------------------- //
+  N_Newton = 1;
+
+  // ------------------------------------- //
+  // Set the vacuum oscillation parameters //
+  // ------------------------------------- //
+  s12sq = OscParams[kTH12];
+  s13sq = OscParams[kTH13];
+  s23sq = OscParams[kTH23];
+  delta = OscParams[kDCP];
+  Dmsq21 = OscParams[kDM12];
+
+  //Need to convert OscParams[kDM23] to kDM31
+  Dmsq31 = 2.5e-3; // eV^2
+
+  // ------------------------------------------ //
+  // Calculate all 9 oscillationa probabilities //
+  // ------------------------------------------ //
+
+  for (int iOscProb=0;iOscProb<fNEnergyPoints;iOscProb++) {    
+    for (int iNuType=0;iNuType<fNNeutrinoTypes;iNuType++) {
+
+      //+ve energy for neutrinos, -ve energy for antineutrinos
+      E = fEnergyArray[iOscProb] * fNeutrinoTypes[iNuType];
+
+      Probability_Matter_LBL(s12sq, s13sq, s23sq, delta, Dmsq21, Dmsq31, L, E, rho, Ye, N_Newton, &probs_returned);
+
+      for (int iOscChannel=0;iOscChannel<fNOscillationChannels;iOscChannel++) {
+	// Mapping which links the oscillation channel, neutrino type and energy index to the fWeightArray index
+	int IndexToFill = iNuType*fNOscillationChannels*fNEnergyPoints + iOscChannel*fNEnergyPoints;
+
+
+	double Weight = probs_returned[fOscillationChannels[iOscChannel].GeneratedFlavour-1][fOscillationChannels[iOscChannel].DetectedFlavour-1];
+
+	//Cancel floating point precision
+	if (Weight<0. && Weight>-1e-6) {Weight = 0.;}
+
+	if (Weight<0. || Weight > 1.) {
+	  std::cout << "s12sq:" << s12sq << " s13sq:" << s13sq << " s23sq:" << s23sq << " delta:" << delta << " Dmsq21:" << Dmsq21 << " Dmsq31:" << Dmsq31 << " L:" << L << " E:" << E << " rho:" << rho << " Ye:" << Ye << " N_Newton:" << N_Newton << std::endl;
+	  std::cout << "iOscProb:" << iOscProb << " iNuType:" << iNuType << " iOscChannel:" << iOscChannel << " IndexToFill:" << IndexToFill << " fWeightArray[IndexToFill+iOscProb]:" << Weight << std::endl;
+	  throw;
+	}
+
+	fWeightArray[IndexToFill+iOscProb] = Weight;
+      }
+
+    }
+  }
+
+}
+
+int OscProbCalcerNuFASTLinear::ReturnWeightArrayIndex(int NuTypeIndex, int OscChanIndex, int EnergyIndex, int CosineZIndex) {
+  int IndexToReturn = NuTypeIndex*fNOscillationChannels*fNEnergyPoints + OscChanIndex*fNEnergyPoints + EnergyIndex;
+  return IndexToReturn;
+}
+
+long OscProbCalcerNuFASTLinear::DefineWeightArraySize() {
+  long nCalculationPoints = fNEnergyPoints * fNOscillationChannels * fNNeutrinoTypes;
+  return nCalculationPoints;
+}

OscProbCalcer/OscProbCalcer_NuFASTLinear.h

@@ -0,0 +1,82 @@
+#ifndef __OSCILLATOR_NUFASTLINEAR_H__
+#define __OSCILLATOR_NUFASTLINEAR_H__
+
+#include "OscProbCalcerBase.h"
+
+/**
+ * @file OscProbCalcer_NuFASTLinear.h
+ *
+ * @class OscProbCalcerNuFASTLinear
+ *
+ * @brief Oscillation calculation engine for linear propagation in NuFAST.
+ */
+class OscProbCalcerNuFASTLinear : public OscProbCalcerBase {
+ public:
+
+  /**
+   * @brief Default constructor
+   *
+   * @param ConfigName_ Name of config used to setup the OscProbCalcerNuFASTLinear() instance
+   * @param Instance_ Which entry of the OscProbCalcerSetup config block should be read in the case where there are multiple OscProbCalcers to be initialised
+   */
+  OscProbCalcerNuFASTLinear(std::string ConfigName_="", int Instance_=0);
+
+  // ========================================================================================================================================================================
+  // Functions which need implementation specific code
+
+  /**
+   * @brief Setup NuFAST specific variables
+   */  
+  void SetupPropagator();
+
+  /**
+   * @brief Calculate some oscillation probabilities for a particular oscillation parameter set
+   *
+   * Calculator oscillation probabilities in NuFAST. This function both calculates and stores
+   * the oscillation probabilities in #fWeightArray.
+   *
+   * @param OscParams The parameter set to calculate oscillation probabilities at
+   */
+  void CalculateProbabilities(std::vector<FLOAT_T> OscParams);
+
+  /**
+   * @brief Return implementation specific index in the weight array for a specific combination of neutrino oscillation channel, energy and cosine zenith
+   * 
+   * @param NuTypeIndex The index in #fNeutrinoTypes (neutrino/antinuetrino) to return the pointer for 
+   * @param OscChanIndex The index in #fOscillationChannels to return the pointer for 
+   * @param EnergyIndex The index in #fEnergyArray to return the pointer for 
+   * @param CosineZIndex The index in #fCosineZArray to return the pointer for 
+   *
+   * @return Index in #fWeightArray which corresponds to the given inputs
+   */
+  int ReturnWeightArrayIndex(int NuTypeIndex, int OscNuIndex, int EnergyIndex, int CosineZIndex=-1);
+
+  /**
+   * @brief Define the size of fWeightArray
+   *
+   * This is implementation specific because because NuFAST is setup to calculate all oscillation channels together, whilst others calculate only a single oscillation channel.
+   *
+   * @return Length that #fWeightArray should be initialised to
+   */
+  long DefineWeightArraySize();
+
+  // ========================================================================================================================================================================
+  // Functions which help setup implementation specific code
+
+  // ========================================================================================================================================================================
+  // Variables which are needed for implementation specific code
+
+  /**
+   * @brief Definition of oscillation parameters which are expected in this ProbGPU implementation
+   */
+  //enum OscParams{kTH12, kTH23, kTH13, kDM12, kDM23, kDCP, kPATHL, kDENS, kELECDENS, kNOscParams};
+  enum OscParams{kTH12, kTH23, kTH13, kDM12, kDM23, kDCP, kPATHL, kDENS, kNOscParams};
+
+  /**
+   * @brief Define the neutrino and antineutrino values expected by this implementation
+   */
+  enum NuType{Nu=1,Nubar=-1};
+
+};
+
+#endif

Oscillator/OscillatorBase.cpp

@@ -16,6 +16,10 @@
 #include "OscProbCalcer/OscProbCalcer_ProbGPULinear.h"
 #endif
 
+#if UseNuFASTLinear==1
+#include "OscProbCalcer/OscProbCalcer_NuFASTLinear.h"
+#endif
+
 #include <iostream>
 
 OscillatorBase::OscillatorBase(std::string ConfigName_) {
@@ -163,6 +167,17 @@ OscProbCalcerBase* OscillatorBase::InitialiseOscProbCalcer(std::string OscProbCa
     throw;
 #endif
   }
+
+  else if (OscProbCalcerImplementationToCreate == "NuFASTLinear") {
+#if UseNuFASTLinear==1
+    OscProbCalcerNuFASTLinear* NuFASTLinear = new OscProbCalcerNuFASTLinear(OscProbCalcerConfigname,Instance);
+    Calcer = (OscProbCalcerBase*)NuFASTLinear;
+    if (fVerbose >= INFO) {std::cout << "Initalised OscProbCalcer Implementation:" << Calcer->ReturnImplementationName() << " in OscillatorBase object" << std::endl;}
+#else
+    std::cerr << "Oscillator was requsted to create " << OscProbCalcerImplementationToCreate << " OscProbCalcer but Use" << OscProbCalcerImplementationToCreate << " is undefined. Indicates problem in setup" << std::endl;
+    throw;
+#endif
+  }
 
   else {
     std::cerr << "Oscillator was requsted to create " << OscProbCalcerImplementationToCreate << " OscProbCalcer but this is not implemented within " << __FILE__ << std::endl;