OscCalcPMNSOptEigen.h
Go to the documentation of this file.
1 #ifndef OSC_OSCCALCULATORPMNSOPTEIGEN_H
2 #define OSC_OSCCALCULATORPMNSOPTEIGEN_H
3 
4 #include "OscLib/IOscCalc.h"
5 #include "OscLib/OscParameters.h"
6 #include <complex>
7 #include <vector>
8 
9 
10 #include <Eigen/Dense>
11 #include <Eigen/Eigenvalues>
12 
13 #include <iostream>
14 
15 
16  // Some useful complex numbers
17  static std::complex<double> zero(0.0,0.0);
18  static std::complex<double> one (1.0,0.0);
19 
20  // Unit conversion constants
21  static const double kKm2eV = 5.06773103202e+09; ///< km to eV^-1
22  static const double kK2 = 4.62711492217e-09; ///< mole/GeV^2/cm^3 to eV
23  static const double kGeV2eV = 1.0e+09; ///< GeV to eV
24 
25  //G_F in units of GeV^-2
26  static const double kGf = 1.166371e-5;
27 
28 namespace osc
29 {
30  /// \brief Helper struct for the cache. Might not need this
31 // struct OscParameters {
32 // double Dmsq31() const { return dmsq32 + dmsq21; }
33 // bool operator==(const OscParameters & rhs) const {
34 // return dmsq21 == rhs.dmsq21 &&
35 // dmsq32 == rhs.dmsq32 &&
36 // th12 == rhs.th12 &&
37 // th13 == rhs.th13 &&
38 // th23 == rhs.th23 &&
39 // deltacp == rhs.deltacp &&
40 // L == rhs.L &&
41 // rho == rhs.rho; }
42 // double dmsq21 = 0;
43 // double dmsq32 = 0;
44 // double th12 = 0;
45 // double th13 = 0;
46 // double th23 = 0;
47 // double deltacp = 0;
48 // double L = 0;
49 // double rho = 0;
50 // };
51 //
52 // /// \brief Struct of the cache
53 // struct OscCache {
54 // bool operator==(const OscCache & rhs) const {
55 // if(!(parameters == rhs.parameters)) return false;
56 // for(unsigned int i = 0; i < energies.size(); i++)
57 // if(energies[i] != rhs.energies[i]) return false;
58 // return true;
59 // }
60 //
61 // std::vector<double> energies;
62 // Eigen::ArrayXXd probabilities;
63 // OscParameters parameters;
64 // unsigned int iter = 0;
65 // };
66 
67  /// \brief Wrapper of a system of eigenvectors and eigenvalues
68  struct EigenSystem {
69  const Eigen::Matrix3cd vectors;
70  const Eigen::Vector3d values;
71  };
72 
73  /// \brief A re-optimized version of \ref OscCalcPMNSOpt
74  ///
75  /// Uses a faster caching scheme than OscCalcPMNSOpt
77  {
78  public:
80  OscCalcPMNSOptEigen(std::vector<double> energies) {
81  this->SetCachedEnergies(energies);
82  }
83  ~OscCalcPMNSOptEigen() = default;
84 
85  IOscCalcAdjustable * Copy() const override;
86 
87  double P(int flavBefore, int flavAfter, double E) override;
88  double P(int flavBefore, int flavAfter, double E, bool fast_and_loose);
89  Eigen::ArrayXd P(int flavBefore, int flavAfter,
90  const std::vector<double>& E) override;
91 
92 
93  void SetL (double L ) override {SaveLastParams(); fL = L;}
94  void SetRho (double rho ) override {SaveLastParams(); fRho = rho;}
95  void SetDmsq21(const double& dmsq21) override {SaveLastParams(); fDmsq21 = dmsq21;}
96  void SetDmsq32(const double& dmsq32) override {SaveLastParams(); fDmsq32 = dmsq32;}
97  void SetTh12 (const double& th12 ) override {SaveLastParams(); fTh12 = th12;}
98  void SetTh13 (const double& th13 ) override {SaveLastParams(); fTh13 = th13;}
99  void SetTh23 (const double& th23 ) override {SaveLastParams(); fTh23 = th23;}
100  void SetdCP (const double& dCP ) override {SaveLastParams(); fdCP = dCP;}
101 
102  TMD5* GetParamsHash() const override;
103  std::string Name() const { return name; }
104 
105  private:
106  void FillCache(const std::vector<double> & energies);// move back to private
107  _OscCache<double> fCache; // move back to private
108 
109  std::string name = "OscCalcPMNSOptEigen";
110  int ChannelCacheIdx(int flavBefore, int flavAfter) const;
111 
112  // Fill the cache at the current parameter values
113  void FillCache();
114 
115 
116  // update fLastParams with the current parameters before changing them
117  void SaveLastParams();
118 
119  void SetCachedEnergies(std::vector<double> const & energies);
120 
121  // Build a flavor basis hamiltonian at input parameters. No matter effects
122  Eigen::Matrix3cd BuildHam(const OscParameters & params) const;
123 
124  // Add matter effects onto a flavor basis vacuum hamiltonian
125  Eigen::Matrix3cd AddMatterEffects(const Eigen::Matrix3cd & ham,
126  const double & E,
127  const int & anti,
128  const OscParameters & params) const;
129  Eigen::Matrix3cd AddMatterEffects(const Eigen::Matrix3cd & ham,
130  const double & E,
131  const OscParameters & params) const;
132  Eigen::Matrix3cd AddMatterEffectsAnti(const Eigen::Matrix3cd & ham,
133  const double & E,
134  const OscParameters & params) const;
135 
136  // Build a flavor basis hamiltonian at input parameters. With matter effects
137  Eigen::Matrix3cd MatterHamiltonian(const double & E,
138  const int & anti,
139  const OscParameters & params) const;
140 
141  // Solve the eigenvalue problem for the input hamiltonian
142  EigenSystem Solve(Eigen::Matrix3cd const & ham) const;
143 
144  // rotate eigenvectors through matter
145  Eigen::Matrix3cd PropMatter(Eigen::Matrix3cd const & evec,
146  Eigen::Vector3d const & evals,
147  const OscParameters & params) const;
148 
149  // Calculate a matrix of probabilities corresponding to all
150  // permutations of to and from states.
151  // Note you probably want "to" to be propagated
152  //Eigen::Array33d Probabilities(const Eigen::Matrix3cd & to,
153  //const Eigen::Matrix3cd & from) const;
154 // Eigen::Array33d Probabilities(Eigen::Matrix3cd const & evec,
155 // Eigen::Vector3d const & evals,
156 // const OscParameters & params);
157 
158 
159  void PrintArrayAddresses(const double * data, int size) const;
160  void PrintMatrixAddresses(const double * data,
161  int nrows,
162  int ncols) const;
163  void PrintMatrixAddresses(const Eigen::ArrayXXd & mat) const;
164  bool ParamsAreCached();
165 
167 
168  };
169 
170 
171 
172 }
173 
174 
175 #endif
static std::complex< double > one(1.0, 0.0)
const XML_Char * name
Definition: expat.h:151
static const double kK2
mole/GeV^2/cm^3 to eV
static const double kKm2eV
km to eV^-1
double th23
Definition: runWimpSim.h:98
std::string Name() const
void SetDmsq21(const double &dmsq21) override
double th12
Definition: runWimpSim.h:98
static const double kGf
OscCalcPMNSOptEigen(std::vector< double > energies)
void SetTh13(const double &th13) override
const XML_Char const XML_Char * data
Definition: expat.h:268
#define P(a, b, c, d, e, x)
void SetL(double L) override
void SetdCP(const double &dCP) override
static constexpr double L
Float_t E
Definition: plot.C:20
Helper struct for the cache. Might not need this.
double dCP
const Eigen::Vector3d values
_OscCache< double > fCache
A re-optimized version of OscCalcPMNSOpt.
Oscillation probability calculators.
Definition: Calcs.h:5
void SetTh12(const double &th12) override
double dmsq32
void SetRho(double rho) override
static std::complex< double > zero(0.0, 0.0)
const Eigen::Matrix3cd vectors
General interface to any calculator that lets you set the parameters.
void SetDmsq32(const double &dmsq32) override
Int_t ncols
Definition: plot.C:53
static const double kGeV2eV
GeV to eV.
double th13
Definition: runWimpSim.h:98
void SetTh23(const double &th23) override
Eigen::MatrixXd mat
enum BeamMode string