PMNSOpt.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // $Id: PMNSOpt.cxx,v 1.1 2013/01/19 16:09:57 jcoelho Exp $
3 //
4 // Implementation of oscillations of neutrinos in matter in a
5 // three-neutrino framework. Two optimizations are relevant:
6 // The construction of the Hamiltonian follows DocDB-XXXX (to be posted)
7 // The eigensystem determination is based on the following reference:
8 //
9 //......................................................................
10 //
11 // Int. J. Mod. Phys. C VOLUME 19, NUMBER 03 MARCH 2008
12 //
13 // Efficient numerical diagonalization of hermitian 3x3 matrices
14 //
15 // Joachim Kopp
16 // Max–Planck–Institut für Kernphysik
17 // Postfach 10 39 80, 69029 Heidelberg, Germany
18 // (Received 19 October 2007)
19 //
20 // 523
21 //......................................................................
22 //
23 // Throughout I have taken:
24 // - L to be the neutrino flight distance in km
25 // - E to be the neutrino energy in GeV
26 // - Dmij to be the differences between the mass-squares in eV^2
27 // - Ne to be the electron number density in mole/cm^3
28 // - theta12,theta23,theta13,deltaCP to be in radians
29 //
30 // The code structure follows the implementation written by M. Messier
31 // in the PMNS class.
32 //
33 // joao.coelho@tufts.edu
34 ////////////////////////////////////////////////////////////////////////
35 
36 // Note from templating over stan::math::var
37 // -----------------------------------------
38 // BEWARE. stan::math::var default constructor has a nullptr in its content,
39 // so if you create a std::complex<stan::math::var>, both real & imag components
40 // are initially nullptr.
41 // If you write
42 // std::complex<stan::math::var> zero = 0.0;
43 // you will now have a proper value in the real part, but the imaginary part
44 // will still be null!
45 // That can cause segfaults when we try to perform std::complex overloaded operations
46 // (e.g. arithmetic) that assume that the imaginary part is defined.
47 // Herein we attempt to explicitly call the std::complex constructor with a pair of zeros
48 // std::complex<stan::math::var> zero(0.0, 0.0);
49 // etc. to avoid this happening.
50 //
51 // There is one unavoidable edge case, however.
52 // The overload of std::abs() for std::complex<> internally uses the following test:
53 // if (__s == _Tp())
54 // return __s;
55 // where __s is either the real or imaginary part of the number,
56 // and _Tp is the type std::complex is templated over.
57 // Because the way stan::math::var's operator== is currently written,
58 // it tries to dereference the nullptr inside the anonymous stan::math::var
59 // returned by calling the constructor to compare to __s's value...
60 // causing a segfault.
61 // The workaround is to simply write your own magnitude function and use that.
62 // i.e.
63 // template<typename T>
64 // double Mag(const std::complex<T>& z) { return z.real() * z.real() + z.imag() * z.imag(); }
65 // This was done in the zhe***3 libraries included below.
66 
67 #include "OscLib/PMNSOpt.h"
68 
73 
74 #include <cstdlib>
75 #include <cassert>
76 
77 using namespace osc;
78 
79 //......................................................................
80 
81 template <typename T>
83 {
84  this->SetMix(0.,0.,0.,0.);
85  this->SetDeltaMsqrs(0.,0.);
86  this->ResetToFlavour(1);
87  fCachedNe = 0.0;
88  fCachedE = 1.0;
89  fCachedAnti = 1;
90  fBuiltHlv = false;
91 }
92 
93 
94 //......................................................................
95 
96 template <typename T>
97 void _PMNSOpt<T>::SetMix(const T &th12, const T &th23, const T &th13, const T &deltacp)
98 {
99 
100  fTheta12 = th12;
101  fTheta23 = th23;
102  fTheta13 = th13;
103  fDeltaCP = deltacp;
104 
105  fBuiltHlv = false;
106 
107 }
108 
109 //......................................................................
110 ///
111 /// Set the mass-splittings. These are m_2^2-m_1^2, m_3^2-m_2^2
112 /// and m_3^2-m_1^2 in eV^2
113 ///
114 template <typename T>
115 void _PMNSOpt<T>::SetDeltaMsqrs(const T &dm21, const T &dm32)
116 {
117 
118  fDm21 = dm21;
119  fDm31 = dm32 + dm21;
120 
121  if(fDm31==0) fDm31 = 1.0e-12;
122 
123  fBuiltHlv = false;
124 
125 }
126 
127 //......................................................................
128 ///
129 /// Build H*lv, where H is the Hamiltonian in vacuum on flavour basis
130 /// and lv is the oscillation length
131 ///
132 /// This is a dimentionless hermitian matrix, so only the
133 /// upper triangular part needs to be filled
134 ///
135 /// The construction of the Hamiltonian avoids computing terms that
136 /// are simply zero. This has a big impact in the computation time.
137 /// This construction is described in DocDB-XXXX (to be posted)
138 ///
139 template <typename T>
141 {
142 
143  // Check if anything changed
144  if(fBuiltHlv) return;
145 
146  // Create temp variables
147  T sij, cij, h00, h11, h01;
148  complex expCP(0,0), h02(0,0), h12(0,0);
149 
150  // Hamiltonian in mass base. Only one entry is variable.
151  h11 = fDm21 / fDm31;
152 
153  // Rotate over theta12
154  sij = sin(fTheta12);
155  cij = cos(fTheta12);
156 
157  // There are 3 non-zero entries after rephasing so that h22 = 0
158  h00 = h11 * sij * sij - 1;
159  h01 = h11 * sij * cij;
160  h11 = h11 * cij * cij - 1;
161 
162  // Rotate over theta13 with deltaCP
163  sij = sin(fTheta13);
164  cij = cos(fTheta13);
165  expCP = complex(cos(fDeltaCP), -sin(fDeltaCP));
166 
167  // There are 5 non-zero entries after rephasing so that h22 = 0
168  h02 = (-h00 * sij * cij) * expCP;
169  h12 = (-h01 * sij) * expCP;
170  h11 -= h00 * sij * sij;
171  h00 *= cij * cij - sij * sij;
172  h01 *= cij;
173 
174  // Finally, rotate over theta23
175  sij = sin(fTheta23);
176  cij = cos(fTheta23);
177 
178  // Fill the Hamiltonian rephased so that h22 = -h11
179  // explicit construction of complex vars to avoid problem noted at top of file
180  fHlv[0][0] = complex(h00 - 0.5 * h11, 0);
181  fHlv[1][1] = complex(0.5 * h11 * (cij * cij - sij * sij) + 2 * real(h12) * cij * sij, 0);
182  fHlv[2][2] = -fHlv[1][1];
183 
184  // these are all constructed from complex variables (h02 and h12) so they are ok
185  fHlv[0][1] = h02 * sij + h01 * cij;
186  fHlv[0][2] = h02 * cij - h01 * sij;
187  fHlv[1][2] = h12 - (h11 * cij + 2 * real(h12) * sij) * sij;
188 
189  // Tag as built
190  fBuiltHlv = true;
191 
192 }
193 
194 //......................................................................
195 ///
196 /// Solve the full Hamiltonian for eigenvectors and eigenvalues
197 /// This is using a method from the GLoBES software available at
198 /// http://www.mpi-hd.mpg.de/personalhomes/globes/3x3/
199 /// We should cite them accordingly
200 ///
201 template <typename T>
202 void _PMNSOpt<T>::SolveHam(double E, double Ne, int anti)
203 {
204 
205  // Check if anything has changed before recalculating
206  if(Ne!=fCachedNe || E!=fCachedE || anti!=fCachedAnti || !fBuiltHlv ){
207  fCachedNe = Ne;
208  fCachedE = E;
209  fCachedAnti = anti;
210  this->BuildHlv();
211  }
212  else return;
213 
214  auto lv = 2 * kGeV2eV*E / fDm31; // Osc. length in eV^-1
215  auto kr2GNe = kK2*M_SQRT2*kGf*Ne; // Matter potential in eV
216 
217  // Finish building Hamiltonian in matter with dimension of eV
218  complex A[3][3];
219  for(int i=0;i<3;i++){
220  A[i][i] = fHlv[i][i]/lv;
221  for(int j=i+1;j<3;j++){
222  if(anti>0) A[i][j] = fHlv[i][j]/lv;
223  else A[i][j] = conj(fHlv[i][j])/lv;
224  }
225  }
226  if(anti>0) A[0][0] += kr2GNe;
227  else A[0][0] -= kr2GNe;
228 
229  // Solve Hamiltonian for eigensystem using the GLoBES method
230  zheevh3(A,fEvec,fEval);
231 
232 }
233 
234 ///.....................................................................
235 ///
236 /// Propagate the current neutrino state over a distance L in km
237 /// with an energy E in GeV through constant matter of density
238 /// Ne in mole/cm^3.
239 /// @param anti - +1 = neutrino case, -1 = anti-neutrino case
240 ///
241 template <typename T>
242 void _PMNSOpt<T>::PropMatter(double L, double E, double Ne, int anti)
243 {
244 
245  // Solve Hamiltonian
246  this->SolveHam(E, Ne, anti);
247 
248  // Store coefficients of propagation eigenstates
249  complex nuComp[3];
250 
251  for(int i=0;i<3;i++){
252  nuComp[i] = complex(0,0);
253  for(int j=0;j<3;j++){
254  nuComp[i] += fNuState[j] * conj(fEvec[j][i]);
255  }
256  }
257 
258  for(int i=0;i<3;i++)fNuState[i] = complex(0,0);
259 
260  // Propagate neutrino state
261  for(int j=0;j<3;j++){
262  auto s = sin(-fEval[j] * kKm2eV * L);
263  auto c = cos(-fEval[j] * kKm2eV * L);
264 
265  complex jPart = complex(c, s) * nuComp[j];
266  for(int i=0;i<3;i++){
267  fNuState[i] += jPart * fEvec[i][j];
268  }
269  }
270 
271 }
272 
273 //......................................................................
274 ///
275 /// Do several layers in a row. L and Ne must have the same length
276 ///
277 template <typename T>
278 void _PMNSOpt<T>::PropMatter(const std::list<double>& L,
279  double E,
280  const std::list<double>& Ne,
281  int anti)
282 {
283  if (L.size()!=Ne.size()) abort();
284  auto Li = L.begin();
285  auto Lend = L.end();
286  auto Ni = Ne.begin();
287  for (; Li!=Lend; ++Li, ++Ni) {
288  // For very low densities, use vacumm
289  static const double kRhoCutoff = 1.0E-6; // Ne in moles/cm^3
290  if (*Ni<kRhoCutoff) this->PropVacuum(*Li, E, anti);
291  else this->PropMatter(*Li, E, *Ni, anti);
292  }
293 }
294 
295 
296 ///.....................................................................
297 ///
298 /// We know the vacuum eigensystem, so just write it explicitly
299 /// The eigenvalues depend on energy, so E needs to be provided in GeV
300 /// @param anti - +1 = neutrino case, -1 = anti-neutrino case
301 ///
302 template <typename T>
303 void _PMNSOpt<T>::SetVacuumEigensystem(double E, int anti)
304 {
305 
306  T s12, s23, s13, c12, c23, c13;
307  complex expidelta(cos(fDeltaCP), anti * sin(fDeltaCP));
308 
309  s12 = sin(fTheta12); s23 = sin(fTheta23); s13 = sin(fTheta13);
310  c12 = cos(fTheta12); c23 = cos(fTheta23); c13 = cos(fTheta13);
311 
312  fEvec[0][0] = complex(c12*c13, 0);
313  fEvec[0][1] = complex(s12*c13, 0);
314  fEvec[0][2] = s13*conj(expidelta);
315 
316  fEvec[1][0] = -s12*c23-c12*s23*s13*expidelta;
317  fEvec[1][1] = c12*c23-s12*s23*s13*expidelta;
318  fEvec[1][2] = complex(s23*c13, 0);
319 
320  fEvec[2][0] = s12*s23-c12*c23*s13*expidelta;
321  fEvec[2][1] = -c12*s23-s12*c23*s13*expidelta;
322  fEvec[2][2] = complex(c23*c13, 0);
323 
324  fEval[0] = 0;
325  fEval[1] = fDm21 / (2 * kGeV2eV*E);
326  fEval[2] = fDm31 / (2 * kGeV2eV*E);
327 
328 }
329 
330 ///.....................................................................
331 ///
332 /// Propagate the current neutrino state over a distance L in km
333 /// with an energy E in GeV through vacuum
334 /// @param anti - +1 = neutrino case, -1 = anti-neutrino case
335 ///
336 template <typename T>
337 void _PMNSOpt<T>::PropVacuum(double L, double E, int anti)
338 {
339 
340  this->SetVacuumEigensystem(E, anti);
341 
342  complex nuComp[3];
343 
344  for(int i=0;i<3;i++){
345  nuComp[i] = 0;
346  for(int j=0;j<3;j++){
347  nuComp[i] += fNuState[j] * conj(fEvec[j][i]);
348  }
349  }
350 
351  const T km2EvL = kKm2eV*L; // needed for the templated multiplication below to work
352  for(int i=0;i<3;i++){
353  fNuState[i] = 0;
354  for(int j=0;j<3;j++){
355  complex iEval(0.0,fEval[j]);
356  fNuState[i] += exp(-iEval * km2EvL) * nuComp[j] * fEvec[i][j];
357  }
358  }
359 
360 }
361 
362 //......................................................................
363 ///
364 /// Reset the neutrino state back to a pure flavour where
365 /// it starts
366 ///
367 template <typename T>
369 {
370  int i;
371  for (i=0; i<3; ++i){
372  if (i==flv) fNuState[i] = complex(1, 0);
373  else fNuState[i] = complex(0, 0);
374  }
375 }
376 
377 //......................................................................
378 ///
379 /// Compute oscillation probability of flavour flv
380 ///
381 /// 0 = nue, 1 = numu, 2 = nutau
382 ///
383 template <typename T>
384 T _PMNSOpt<T>::P(int flv) const
385 {
386  assert(flv>=0 && flv<3);
387  return norm(fNuState[flv]);
388 }
389 
390 ////////////////////////////////////////////////////////////////////////
391 // manually instantiate templates for those cases we know about.
392 
393 template class osc::_PMNSOpt<double>;
394 
395 #ifdef OSCLIB_STAN
396 #include "stan/math/rev/scal.hpp"
397 template class osc::_PMNSOpt<stan::math::var>;
398 #endif
virtual void BuildHlv()
Definition: PMNSOpt.cxx:140
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
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
int zheevh3(std::complex< T > A[3][3], std::complex< T > Q[3][3], T w[3])
Definition: zheevh3.h:44
double th12
Definition: runWimpSim.h:98
virtual void SetMix(const T &th12, const T &th23, const T &th13, const T &deltacp)
Definition: PMNSOpt.cxx:97
const XML_Char * s
Definition: expat.h:262
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
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
const double j
Definition: BetheBloch.cxx:29
Oscillation probability calculators.
Definition: Calcs.h:5
TH1F * h11
Definition: plot.C:43
static const double A
Definition: Units.h:82
T sin(T number)
Definition: d0nt_math.hpp:132
Float_t norm
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
TH1D * h12
Definition: plot_hist.C:27
virtual void PropMatter(double L, double E, double Ne, int anti=1)
Definition: PMNSOpt.cxx:242
virtual T P(int flv) const
Definition: PMNSOpt.cxx:384
double T
Definition: Xdiff_gwt.C:5
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