PMNS_Base.h

Go to the documentation of this file.

 
 
#ifndef PMNS_BASE_H
#define PMNS_BASE_H
 
#include <unordered_set>
 
#include "EigenPoint.h"
 
namespace OscProb {
 
  class PMNS_Base {
    public:
      PMNS_Base(int numNus = 3); 
      virtual ~PMNS_Base();      
 
      // Get the oscillation probability
      virtual double Prob(
          vectorC nu_in,
          int     flvf); 
      virtual double Prob(vectorC nu_in, int flvf,
                          double E); 
      virtual double Prob(vectorC nu_in, int flvf, double E,
                          double L); 
 
      // Get the oscillation probability
      virtual double Prob(
          int flvi,
          int flvf); 
      virtual double Prob(int flvi, int flvf,
                          double E); 
      virtual double Prob(int flvi, int flvf, double E,
                          double L); 
 
      // Get the oscillation probability vector
      virtual vectorD ProbVector(
          vectorC nu_in); 
      virtual vectorD ProbVector(
          vectorC nu_in,
          double  E); 
      virtual vectorD ProbVector(
          vectorC nu_in, double E,
          double L); 
 
      // Get the oscillation probability vector
      virtual vectorD ProbVector(int flvi); 
      virtual vectorD ProbVector(
          int flvi, double E); 
      virtual vectorD ProbVector(
          int flvi, double E,
          double L); 
 
      // Get the oscillation probability matrix
      virtual matrixD ProbMatrix(int nflvi,
                                 int nflvf); 
      virtual matrixD ProbMatrix(
          int nflvi, int nflvf,
          double E); 
      virtual matrixD ProbMatrix(int nflvi, int nflvf, double E,
                                 double L); 
 
      // Get probability averaged over a bin
      virtual double AvgProb(
          vectorC nu_in, int flvf, double E,
          double dE =
              0); 
      virtual double AvgProbLoE(
          vectorC nu_in, int flvf, double LoE,
          double dLoE =
              0); 
 
      // Get probability averaged over a bin
      virtual double AvgProb(
          int flvi, int flvf, double E,
          double dE =
              0); 
      virtual double AvgProbLoE(
          int flvi, int flvf, double LoE,
          double dLoE =
              0); 
 
      // Get probability vector averaged over a bin
      virtual vectorD AvgProbVector(
          vectorC nu_in, double E,
          double dE = 0); 
      virtual vectorD AvgProbVectorLoE(
          vectorC nu_in, double LoE,
          double dLoE =
              0); 
 
      // Get probability vector averaged over a bin
      virtual vectorD AvgProbVector(
          int flvi, double E,
          double dE = 0); 
      virtual vectorD AvgProbVectorLoE(
          int flvi, double LoE,
          double dLoE =
              0); 
 
      // Get probability vector averaged over a bin
      virtual matrixD AvgProbMatrix(
          int nflvi, int nflvf, double E,
          double dE = 0); 
      virtual matrixD AvgProbMatrixLoE(
          int nflvi, int nflvf, double LoE,
          double dLoE =
              0); 
 
      virtual vectorC GetMassEigenstate(
          int mi); 
 
      // Set the oscillation parameters
      virtual void SetAngle(int i, int j,
                            double th); 
      virtual void SetDelta(int i, int j,
                            double delta); 
      virtual void SetDm(int    j,
                         double dm); 
 
      // Get the oscillation parameters
      virtual double GetAngle(int i, int j); 
      virtual double GetDelta(int i, int j); 
      virtual double GetDm(int j); 
 
      // Get the effective oscillation parameters
      virtual double GetDmEff(
          int j); 
 
      // Set default oscillation parameters
      virtual void SetStdPars(); 
 
      // Set energy and anti-neutrino flag
      virtual void SetEnergy(double E);      
      virtual void SetIsNuBar(bool isNuBar); 
 
      // Get energy and anti-neutrino flag
      virtual double GetEnergy();  
      virtual bool   GetIsNuBar(); 
 
      // Set the neutrino path
      virtual void SetPath(NuPath p); 
      virtual void SetPath(double length, double density, double zoa = 0.5,
                           int layer = 0); 
 
      virtual void SetPath(std::vector<NuPath> paths); 
 
      virtual void AddPath(NuPath p); 
      virtual void AddPath(double length, double density, double zoa = 0.5,
                           int layer = 0); 
 
      virtual void ClearPath(); 
 
      virtual void SetLength(double L); 
      virtual void SetDensity(
          double rho);                 
      virtual void SetZoA(double zoa); 
 
      virtual void SetLength(vectorD L);    
      virtual void SetDensity(vectorD rho); 
      virtual void SetZoA(vectorD zoa);     
      virtual void SetLayers(
          std::vector<int> lay); 
 
      // Set a default neutrino path
      virtual void SetStdPath(); 
 
      // Get the neutrino path
      virtual std::vector<NuPath> GetPath(); 
 
      virtual vectorD GetSamplePoints(double LoE, double dLoE);
 
      // Setup the caching system
      virtual void SetUseCache(bool u = true); 
      virtual void ClearCache();               
      virtual void SetMaxCache(int mc = 1e6);  
 
      virtual void SetAvgProbPrec(double prec);
 
    protected:
      // Some useful complex numbers
      static const complexD zero; 
      static const complexD one;  
 
      // Unit conversion constants
      static const double kKm2eV;  
      static const double kK2;     
      static const double kGeV2eV; 
      static const double kNA;     
 
      static const double kGf; 
 
      virtual void InitializeVectors(); 
 
      // Internal caching functions
      virtual bool TryCache();  
      virtual void FillCache(); 
 
      // Auxiliary path functions
      virtual void SetCurPath(
          NuPath p); 
 
      virtual void SetAtt(double att,
                          int    idx); 
      virtual void SetAtt(vectorD att,
                          int     idx); 
 
      // Building and solving
      virtual void RotateH(
          int i, int j,
          matrixC& Ham); 
      virtual void RotateState(
          int i, int j); 
 
      virtual void BuildHms(); 
 
      virtual void SolveHam() = 0; 
 
      // Resetting and propagating
      virtual void ResetToFlavour(
          int flv); 
      virtual void SetPureState(
          vectorC nu_in); 
 
      virtual void PropagatePath(
          NuPath p);            
      virtual void Propagate(); 
 
      virtual double P(
          int flv); 
      virtual vectorD GetProbVector(); 
 
      virtual std::vector<int> GetSortedIndices(
          const vectorD x); 
 
      virtual vectorD ConvertEtoLoE(double E, double dE);
 
      // Attributes
 
      int fNumNus; 
 
      vectorD fDm;    
      matrixD fTheta; 
      matrixD fDelta; 
 
      vectorC fNuState; 
      matrixC fHms;     
 
      vectorC fPhases; 
      vectorC fBuffer; 
 
      vectorD fEval; 
      matrixC fEvec; 
 
      double fEnergy;  
      bool   fIsNuBar; 
 
      std::vector<NuPath> fNuPaths; 
      NuPath              fPath;    
 
      bool fBuiltHms; 
      bool fGotES;    
 
      bool   fUseCache;  
      double fCachePrec; 
      int    fMaxCache;  
 
      double fAvgProbPrec; 
 
      std::unordered_set<EigenPoint> fMixCache; 
      EigenPoint                     fProbe;    
  };
 
  struct IdxCompare {
      IdxCompare(const vectorD& target) : target(target) {}
 
      bool operator()(int a, int b) const { return target[a] < target[b]; }
 
    private:
      const vectorD target;
  };
 
} // namespace OscProb
 
#endif