OscillatorConstants.h

Go to the documentation of this file.

#ifndef __OSCILLATOR_CONSTANTS__
#define __OSCILLATOR_CONSTANTS__
 
#define DUMMYVAL -999
 
#include <math.h>
 
#include "TFile.h"
#include "TH1.h"
 
#if UseDoubles==1
using FLOAT_T = double;
#else
using FLOAT_T = float;
#endif
 
#include <string>
#include <iostream>
#include <vector>
 
namespace NuOscillator
{
  enum Verbosity{NONE=0,INFO=1,VERBOSE=2};
  
  enum NeutrinoFlavours{kElectron=1,kMuon=2,kTau=3,kSterile1=4,kSterile2=5,kSterile3=6,nNeutrinoFlavours=7};
  
  enum {kNuFlavour_Generated=0,kNuFlavour_Detected=1,nNuFlavours=2};
  
  struct OscillationChannel{
    int GeneratedFlavour;
    int DetectedFlavour;
  };
  
  struct OscillationProbability{
    int NuType;
    OscillationChannel OscChan;
    FLOAT_T Energy;
    FLOAT_T CosineZ;
    FLOAT_T Probability;
  };
  
}
 
inline std::vector<FLOAT_T> ReturnOscParams_Basic() {
  std::vector<FLOAT_T> OscParams_Basic = {3.07e-1,5.28e-1,2.18e-2,7.53e-5,2.509e-3,-1.601};
  return OscParams_Basic;
}
 
inline std::vector<FLOAT_T> ReturnOscParams_Atm() {
  std::vector<FLOAT_T> OscParams_Atm = {3.07e-1,5.28e-1,2.18e-2,7.53e-5,2.509e-3,-1.601,25.0};
  return OscParams_Atm;
}
 
inline std::vector<FLOAT_T> ReturnOscParams_Beam_woYe() {
  std::vector<FLOAT_T> OscParams_Beam_woYe = {3.07e-1,5.28e-1,2.18e-2,7.53e-5,2.509e-3,-1.601,250.0,2.6};
  return OscParams_Beam_woYe;
}
 
inline std::vector<FLOAT_T> ReturnOscParams_Beam_wYe() {
  std::vector<FLOAT_T> OscParams_Beam_wYe = {3.07e-1,5.28e-1,2.18e-2,7.53e-5,2.509e-3,-1.601,250.0,2.6,0.5};
  return OscParams_Beam_wYe;
}
 
inline std::vector<FLOAT_T> ReturnOscParams_Beam_wYe_wDeco() {
  std::vector<FLOAT_T> OscParams_Beam_wYe_wDeco = {3.07e-1,5.28e-1,2.18e-2,7.53e-5,2.509e-3,-1.601,250.0,2.6,0.5,9.48e-18,2.0,1.0}; //9.48e-18
  return OscParams_Beam_wYe_wDeco;
}
 
inline std::vector<FLOAT_T> ReturnOscParams_Beam_wYe_wLIV() {
  std::vector<FLOAT_T> OscParams_Beam_wYe_wLIV = {3.07e-1,5.28e-1,2.18e-2,7.53e-5,2.509e-3,-1.601,250.0,2.6,0.5,1.0e-23,0.0,1.0e-23,0.0,0.0};
  return OscParams_Beam_wYe_wLIV;
}
 
inline std::vector<FLOAT_T> ReturnOscParams_Beam_wYe_wNSI() {
  //6 PMNS params, 6 Epsilons, 3 Deltas, 3 Coupling strengths, baseline, density, A/Z
  std::vector<FLOAT_T> OscParams_Beam_wYe_wNSI = {3.07e-1,5.28e-1,2.18e-2,7.53e-5,2.509e-3,-1.601,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,250,2.6,0.5};
  return OscParams_Beam_wYe_wNSI;
}
 
inline std::vector<std::string> ReturnKnownConfigs() {
  std::vector<std::string> ConfigNames;
 
#if UseCUDAProb3 == 1
  ConfigNames.push_back("./NuOscillatorConfigs/Binned_CUDAProb3.yaml");
#endif
 
#if UseCUDAProb3Linear == 1
  ConfigNames.push_back("./NuOscillatorConfigs/Binned_CUDAProb3Linear.yaml");
#endif
 
#if UseProbGPULinear == 1
  ConfigNames.push_back("./NuOscillatorConfigs/Binned_ProbGPULinear.yaml");
#endif
 
#if UseProb3ppLinear == 1
  ConfigNames.push_back("./NuOscillatorConfigs/Binned_Prob3ppLinear.yaml");
#endif
 
#if UseNuFASTLinear == 1
  ConfigNames.push_back("./NuOscillatorConfigs/Binned_NuFASTLinear.yaml");
#endif
 
#if UseNuSQUIDSLinear == 1
  ConfigNames.push_back("./NuOscillatorConfigs/Binned_NuSQUIDSLinear.yaml");
#endif
  
#if UseOscProb == 1
  ConfigNames.push_back("./NuOscillatorConfigs/Binned_OscProb.yaml");
#endif  
 
  return ConfigNames;
}
 
inline int NeutrinoFlavour_StrToInt(const std::string& NuFlav) {
  if (NuFlav == "Electron" || NuFlav == "electron") {
    return NuOscillator::kElectron;
  } else if (NuFlav == "Muon" || NuFlav == "muon") {
    return NuOscillator::kMuon;
  } else if (NuFlav == "Tau" || NuFlav == "tau") {
    return NuOscillator::kTau;
  } else if (NuFlav == "Sterile1" || NuFlav == "sterile1") {
    return NuOscillator::kSterile1;
  } else if (NuFlav == "Sterile2" || NuFlav == "sterile2") {
    return NuOscillator::kSterile2;
  } else if (NuFlav == "Sterile3" || NuFlav == "sterile3") {
    return NuOscillator::kSterile3;
  } else {
    std::cerr << "Could not convert input string:" << NuFlav << " to known enum value in NeutrinoFlavours" << std::endl;
    throw std::runtime_error("Invalid setup");
  }
  
  return -1;
}
 
inline std::string NeutrinoFlavour_IntToStr(const int NuFlav) {
  switch (NuFlav) {
  case NuOscillator::kElectron: 
    return "Electron";
  case NuOscillator::kMuon:
    return "Muon";
  case NuOscillator::kTau:
    return "Tau";
  case NuOscillator::kSterile1:
    return "Sterile1";
  case NuOscillator::kSterile2:
    return "Sterile2";
  case NuOscillator::kSterile3:
    return "Sterile3";
  default:
    std::cerr << "Recieved unknown NeutrinoFlavour:" << NuFlav << " which is inconsistent with those in enum 'NeutrinoFlavours'" << std::endl;
    throw std::runtime_error("Invalid setup");
  }
  return "";
}
 
inline int Verbosity_StrToInt(const std::string& Verbosity) {
  if (Verbosity == "NONE") {
    return NuOscillator::NONE;
  } else if (Verbosity == "INFO") {
    return NuOscillator::INFO;
  } else if (Verbosity == "VERBOSE") {
    return NuOscillator::VERBOSE;
  } else {
    std::cerr << "Invalid verbosity provided:" << Verbosity << std::endl;
    throw std::runtime_error("Invalid setup");
  }
  
  return -1;
}
 
inline NuOscillator::OscillationChannel ReturnOscillationChannel(const std::string& InputString) {
  std::string GeneratedFlavour = "";
  std::string DetectedFlavour = "";
 
  size_t Delimiter = InputString.find(":");
  if (Delimiter != std::string::npos) {
    GeneratedFlavour = InputString.substr(0,Delimiter);
    DetectedFlavour = InputString.substr(Delimiter+1,InputString.size());
  } else {
    std::cerr << "Expected a string formatted as: 'GeneratedNeutrinoFlavour:DetectedNeutrinoFlavour'" << std::endl;
    std::cerr << "Received:" << InputString << std::endl;
  }
  
  NuOscillator::OscillationChannel OscChannel = {NeutrinoFlavour_StrToInt(GeneratedFlavour),NeutrinoFlavour_StrToInt(DetectedFlavour)};
  return OscChannel;
}
 
inline std::vector<FLOAT_T> logspace(FLOAT_T Emin, FLOAT_T  Emax, int nDiv) {
  if (nDiv==0) {
    std::cerr << "Requested log spacing distribution with 0 divisions" << std::endl;
    throw std::runtime_error("Invalid setup");
  }
 
  std::vector<FLOAT_T> logpoints(nDiv+1, 0.0);
  logpoints[0]=Emin;
 
  if (Emin == 0.) {
    Emin = 0.01;
  }
 
  FLOAT_T Emin_log,Emax_log;
  Emin_log = log10(Emin);
  Emax_log = log10(Emax);
 
  FLOAT_T step_log = (Emax_log - Emin_log)/FLOAT_T(nDiv);
 
  FLOAT_T EE = Emin_log+step_log;
 
  for (int i=1; i<nDiv; i++) {
    logpoints[i] = pow(10.,EE);
    EE += step_log;
  }
 
  logpoints[nDiv]=Emax;
 
  return logpoints;
}
 
inline std::vector<FLOAT_T> linspace(FLOAT_T Emin, FLOAT_T Emax, int nDiv) {
  if (nDiv==0) {
    std::cerr << "Requested linear spacing distribution with 0 divisions" << std::endl;
    throw std::runtime_error("Invalid setup");
  }
 
  std::vector<FLOAT_T> linpoints(nDiv+1, 0.0);
 
  FLOAT_T step_lin = (Emax - Emin)/FLOAT_T(nDiv);
 
  FLOAT_T EE = Emin;
 
  for (int i=0; i<nDiv; i++) {
    if (fabs(EE)<1e-6) {EE = 0.;}
 
    linpoints[i] = EE;
    EE += step_lin;
  }
 
  linpoints[nDiv] = Emax;
 
  return linpoints;
}
 
inline std::vector<FLOAT_T> ReadBinEdgesFromFile(std::string TFileName, std::string HistogramName, int Verbose=NuOscillator::Verbosity::NONE) {
  std::vector<FLOAT_T> BinEdges;
 
  TFile* File = new TFile(TFileName.c_str());
  if (!File || File->IsZombie()) {
    std::cerr << "Could not find file:" << TFileName << std::endl;
    throw std::runtime_error("Invalid input file.");
  }
 
  TH1* Histogram = (TH1*)File->Get(HistogramName.c_str());
  if (!Histogram) {
    std::cerr << "Could not find Histogram:" << HistogramName << " in File:" << TFileName << std::endl;
    throw std::runtime_error("Invalid input file.");
  }
 
  BinEdges.resize(Histogram->GetNbinsX()+1);
  for (int iBin=0;iBin<=Histogram->GetNbinsX();iBin++) {
    BinEdges[iBin] = Histogram->GetBinLowEdge(iBin+1);
  }
 
  File->Close();
  delete File;
 
  if (Verbose >= NuOscillator::INFO) {
    std::cout << "Bin edges successfully read from File:" << TFileName << " , Histogram:" << HistogramName << " :=" << std::endl;
    for (size_t i=0;i<BinEdges.size();i++) {
      std::cout << BinEdges[i] << ", ";
    }
    std::cout << std::endl;
  }
 
  return BinEdges;
}
 
inline std::vector<FLOAT_T> ReturnBinCentersFromBinEdges(std::vector<FLOAT_T> BinEdges) {
  int nBins = BinEdges.size()-1;
  std::vector<FLOAT_T> BinCenters = std::vector<FLOAT_T>(nBins);
 
  for (int iBin=0;iBin<nBins;iBin++) {
    BinCenters[iBin] = (BinEdges[iBin]+BinEdges[iBin+1])/2.0;
  }
 
  return BinCenters;
}
 
#endif