PMNSOpt.h
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \class _PMNSOpt
3 ///
4 /// \brief Implementation of oscillations of neutrinos in matter in a
5 /// three-neutrino framework.
6 ///
7 /// Two optimizations are relevant:
8 /// The construction of the Hamiltonian follows DocDB-XXXX (to be posted)
9 /// The eigensystem determination is based on the following reference:
10 ///
11 ///......................................................................
12 ///
13 /// Int. J. Mod. Phys. C VOLUME 19, NUMBER 03 MARCH 2008
14 ///
15 /// Efficient numerical diagonalization of hermitian 3x3 matrices
16 ///
17 /// Joachim Kopp
18 /// Max–Planck–Institut für Kernphysik
19 /// Postfach 10 39 80, 69029 Heidelberg, Germany
20 /// (Received 19 October 2007)
21 ///
22 /// 523
23 ///......................................................................
24 ///
25 /// The code structure follows the implementation written by M. Messier
26 /// in the PMNS class.
27 ///
28 /// \version $Id: PMNSOpt.h,v 1.1 2013/01/19 16:09:57 jcoelho Exp $
29 ///
30 /// @author joao.coelho@tufts.edu
31 ////////////////////////////////////////////////////////////////////////
32 #ifndef PMNSOPT_H
33 #define PMNSOPT_H
34 #include <list>
35 #include <complex>
36 
37 namespace osc
38 {
39 
40  // Unit conversion constants
41  static const double kKm2eV = 5.06773103202e+09; ///< km to eV^-1
42  static const double kK2 = 4.62711492217e-09; ///< mole/GeV^2/cm^3 to eV
43  static const double kGeV2eV = 1.0e+09; ///< GeV to eV
44 
45  //G_F in units of GeV^-2
46  static const double kGf = 1.166371e-5;
47 
48  /// Optimized version of \ref PMNS
49  template <typename T>
50  class _PMNSOpt
51  {
52  public:
53  _PMNSOpt();
54  virtual ~_PMNSOpt() = default;
55 
56  /// Set the parameters of the PMNS matrix
57  /// @param th12 - The angle theta_12 in radians
58  /// @param th23 - The angle theta_23 in radians
59  /// @param th13 - The angle theta_13 in radians
60  /// @param deltacp - The CPV phase delta_CP in radians
61  virtual void SetMix(const T &th12, const T &th23, const T &th13, const T &deltacp);
62 
63  /// Set the mass-splittings
64  /// @param dm21 - m2^2-m1^2 in eV^2
65  /// @param dm32 - m3^2-m2^2 in eV^2
66  virtual void SetDeltaMsqrs(const T &dm21, const T &dm32);
67 
68  /// Propagate a neutrino through a slab of matter
69  /// @param L - length of slab (km)
70  /// @param E - neutrino energy in GeV
71  /// @param Ne - electron number density of matter in mole/cm^3
72  /// @param anti - +1 = neutrino case, -1 = anti-neutrino case
73  virtual void PropMatter(double L, double E, double Ne, int anti=1);
74  virtual void PropMatter(const std::list<double>& L,
75  double E,
76  const std::list<double>& Ne,
77  int anti);
78 
79  /// Propagate a neutrino through vacuum
80  /// @param L - length of slab (km)
81  /// @param E - neutrino energy in GeV
82  /// @param anti - +1 = neutrino case, -1 = anti-neutrino case
83  virtual void PropVacuum(double L, double E, int anti=1);
84 
85  /// Return the probability of final state in flavour flv
86  /// @param flv - final flavor (0,1,2) = (nue,numu,nutau)
87  virtual T P(int flv) const;
88 
89  /// Erase memory of neutrino propagate and reset neutrino
90  /// to pure flavour flv. Preserves values of mixing and mass-splittings
91  /// @param flv - final flavor (0,1,2) = (nue,numu,nutau)
92  virtual void ResetToFlavour(int flv=1);
93 
94  protected:
95  // A shorthand...
96  typedef std::complex<T> complex;
97 
98  /// Build H*lv, where H is the Hamiltonian in vacuum on flavour basis
99  /// and lv is the oscillation length
100  virtual void BuildHlv();
101 
102  /// Solve the full Hamiltonian for eigenvectors and eigenvalues
103  /// @param E - neutrino energy in GeV
104  /// @param Ne - electron number density of matter in mole/cm^3
105  /// @param anti - +1 = neutrino case, -1 = anti-neutrino case
106  virtual void SolveHam(double E, double Ne, int anti);
107 
108  /// Set the eigensystem to the analytic solution of the vacuum Hamiltonian
109  /// @param E - neutrino energy in GeV
110  /// @param anti - +1 = neutrino case, -1 = anti-neutrino case
111  virtual void SetVacuumEigensystem(double E, int anti);
112 
113  T fDm21; ///< m^2_2 - m^2_1 in vacuum
114  T fDm31; ///< m^2_3 - m^2_1 in vacuum
115  T fTheta12; ///< theta12 mixing angle
116  T fTheta23; ///< theta23 mixing angle
117  T fTheta13; ///< theta13 mixing angle
118  T fDeltaCP; ///< CP violating phase
119  complex fHlv[3][3]; ///< dimensionless matrix H*lv
120  complex fEvec[3][3]; ///< Eigenvectors of the Hamiltonian
121  T fEval[3]; ///< Eigenvalues of the Hamiltonian
122  complex fNuState[3]; ///< The neutrino current state
123  double fCachedNe; ///< Cached electron density
124  double fCachedE; ///< Cached neutrino energy
125  int fCachedAnti; ///< Cached nu/nubar selector
126  bool fBuiltHlv; ///< Tag to avoid rebuilding Hlv
127  };
128 
129  // preserve older behavior
131 }
132 #endif
133 ////////////////////////////////////////////////////////////////////////
virtual ~_PMNSOpt()=default
T fTheta12
theta12 mixing angle
Definition: PMNSOpt.h:115
virtual void BuildHlv()
Definition: PMNSOpt.cxx:140
T fTheta23
theta23 mixing angle
Definition: PMNSOpt.h:116
virtual void PropVacuum(double L, double E, int anti=1)
Definition: PMNSOpt.cxx:337
static const double kK2
mole/GeV^2/cm^3 to eV
Definition: PMNSOpt.h:42
double th23
Definition: runWimpSim.h:98
double fCachedNe
Cached electron density.
Definition: PMNSOpt.h:123
virtual void ResetToFlavour(int flv=1)
Definition: PMNSOpt.cxx:368
virtual void SetVacuumEigensystem(double E, int anti)
Definition: PMNSOpt.cxx:303
Optimized version of PMNS.
Definition: PMNSOpt.h:50
T fDm21
m^2_2 - m^2_1 in vacuum
Definition: PMNSOpt.h:113
double th12
Definition: runWimpSim.h:98
_PMNSOpt< double > PMNSOpt
Definition: PMNSOpt.h:130
int fCachedAnti
Cached nu/nubar selector.
Definition: PMNSOpt.h:125
virtual void SetMix(const T &th12, const T &th23, const T &th13, const T &deltacp)
Definition: PMNSOpt.cxx:97
virtual void SetDeltaMsqrs(const T &dm21, const T &dm32)
Definition: PMNSOpt.cxx:115
static const double kGeV2eV
GeV to eV.
Definition: PMNSOpt.h:43
virtual void SolveHam(double E, double Ne, int anti)
Definition: PMNSOpt.cxx:202
static constexpr double L
Float_t E
Definition: plot.C:20
std::complex< T > complex
Definition: PMNSOpt.h:96
T fDm31
m^2_3 - m^2_1 in vacuum
Definition: PMNSOpt.h:114
T fDeltaCP
CP violating phase.
Definition: PMNSOpt.h:118
double fCachedE
Cached neutrino energy.
Definition: PMNSOpt.h:124
complex fNuState[3]
The neutrino current state.
Definition: PMNSOpt.h:122
Oscillation probability calculators.
Definition: Calcs.h:5
T fEval[3]
Eigenvalues of the Hamiltonian.
Definition: PMNSOpt.h:121
bool fBuiltHlv
Tag to avoid rebuilding Hlv.
Definition: PMNSOpt.h:126
T fTheta13
theta13 mixing angle
Definition: PMNSOpt.h:117
virtual void PropMatter(double L, double E, double Ne, int anti=1)
Definition: PMNSOpt.cxx:242
complex fEvec[3][3]
Eigenvectors of the Hamiltonian.
Definition: PMNSOpt.h:120
virtual T P(int flv) const
Definition: PMNSOpt.cxx:384
double T
Definition: Xdiff_gwt.C:5
complex fHlv[3][3]
dimensionless matrix H*lv
Definition: PMNSOpt.h:119
double th13
Definition: runWimpSim.h:98
static const double kKm2eV
km to eV^-1
Definition: PMNSOpt.h:41
static const double kGf
Definition: PMNSOpt.h:46