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 // n.b. Stan sets up some type traits that need to be loaded before Eigen is.
68 // Since Eigen gets dragged in via IOscCalc.h we have to get Stan set up before
69 // that is included.
70 #ifdef OSCLIB_STAN
71 #include "OscLib/Stan.h"
72 #endif
73 
74 #include "OscLib/PMNSOpt.h"
75 
80 
81 #include <cstdlib>
82 #include <cassert>
83 
84 using namespace osc;
85 
86 //......................................................................
87 
88 template <typename T>
90 {
91  this->SetMix(0.,0.,0.,0.);
92  this->SetDeltaMsqrs(0.,0.);
93  this->ResetToFlavour(1);
94  fCachedNe = 0.0;
95  fCachedE = 1.0;
96  fCachedAnti = 1;
97  fBuiltHlv = false;
98 }
99 
100 
101 //......................................................................
102 
103 template <typename T>
104 void _PMNSOpt<T>::SetMix(const T &th12, const T &th23, const T &th13, const T &deltacp)
105 {
106 
107  fTheta12 = th12;
108  fTheta23 = th23;
109  fTheta13 = th13;
110  fDeltaCP = deltacp;
111 
112  fBuiltHlv = false;
113 
114 }
115 
116 //......................................................................
117 ///
118 /// Set the mass-splittings. These are m_2^2-m_1^2, m_3^2-m_2^2
119 /// and m_3^2-m_1^2 in eV^2
120 ///
121 template <typename T>
122 void _PMNSOpt<T>::SetDeltaMsqrs(const T &dm21, const T &dm32)
123 {
124 
125  fDm21 = dm21;
126  fDm31 = dm32 + dm21;
127 
128  if(fDm31==0) fDm31 = 1.0e-12;
129 
130  fBuiltHlv = false;
131 
132 }
133 
134 //......................................................................
135 ///
136 /// Build H*lv, where H is the Hamiltonian in vacuum on flavour basis
137 /// and lv is the oscillation length
138 ///
139 /// This is a dimentionless hermitian matrix, so only the
140 /// upper triangular part needs to be filled
141 ///
142 /// The construction of the Hamiltonian avoids computing terms that
143 /// are simply zero. This has a big impact in the computation time.
144 /// This construction is described in DocDB-XXXX (to be posted)
145 ///
146 template <typename T>
148 {
149 
150  // Check if anything changed
151  if(fBuiltHlv) return;
152 
153  // Create temp variables
154  T sij, cij, h00, h11, h01;
155  complex expCP(0,0), h02(0,0), h12(0,0);
156 
157  // Hamiltonian in mass base. Only one entry is variable.
158  h11 = fDm21 / fDm31;
159 
160  // Rotate over theta12
161  sij = sin(fTheta12);
162  cij = cos(fTheta12);
163 
164  // There are 3 non-zero entries after rephasing so that h22 = 0
165  h00 = h11 * sij * sij - 1;
166  h01 = h11 * sij * cij;
167  h11 = h11 * cij * cij - 1;
168 
169  // Rotate over theta13 with deltaCP
170  sij = sin(fTheta13);
171  cij = cos(fTheta13);
172  expCP = complex(cos(fDeltaCP), -sin(fDeltaCP));
173 
174  // There are 5 non-zero entries after rephasing so that h22 = 0
175  h02 = (-h00 * sij * cij) * expCP;
176  h12 = (-h01 * sij) * expCP;
177  h11 -= h00 * sij * sij;
178  h00 *= cij * cij - sij * sij;
179  h01 *= cij;
180 
181  // Finally, rotate over theta23
182  sij = sin(fTheta23);
183  cij = cos(fTheta23);
184 
185  // Fill the Hamiltonian rephased so that h22 = -h11
186  // explicit construction of complex vars to avoid problem noted at top of file
187  fHlv[0][0] = complex(h00 - 0.5 * h11, 0);
188  fHlv[1][1] = complex(0.5 * h11 * (cij * cij - sij * sij) + 2 * real(h12) * cij * sij, 0);
189  fHlv[2][2] = -fHlv[1][1];
190 
191  // these are all constructed from complex variables (h02 and h12) so they are ok
192  fHlv[0][1] = h02 * sij + h01 * cij;
193  fHlv[0][2] = h02 * cij - h01 * sij;
194  fHlv[1][2] = h12 - (h11 * cij + 2 * real(h12) * sij) * sij;
195 
196  // Tag as built
197  fBuiltHlv = true;
198 
199 }
200 
201 //......................................................................
202 ///
203 /// Solve the full Hamiltonian for eigenvectors and eigenvalues
204 /// This is using a method from the GLoBES software available at
205 /// http://www.mpi-hd.mpg.de/personalhomes/globes/3x3/
206 /// We should cite them accordingly
207 ///
208 template <typename T>
209 void _PMNSOpt<T>::SolveHam(double E, double Ne, int anti)
210 {
211 
212  // Check if anything has changed before recalculating
213  if(Ne!=fCachedNe || E!=fCachedE || anti!=fCachedAnti || !fBuiltHlv ){
214  fCachedNe = Ne;
215  fCachedE = E;
216  fCachedAnti = anti;
217  this->BuildHlv();
218  }
219  else return;
220 
221  auto lv = 2 * kGeV2eV*E / fDm31; // Osc. length in eV^-1
222  auto kr2GNe = kK2*M_SQRT2*kGf*Ne; // Matter potential in eV
223 
224  // Finish building Hamiltonian in matter with dimension of eV
225  complex A[3][3];
226  for(int i=0;i<3;i++){
227  A[i][i] = fHlv[i][i]/lv;
228  for(int j=i+1;j<3;j++){
229  if(anti>0) A[i][j] = fHlv[i][j]/lv;
230  else A[i][j] = conj(fHlv[i][j])/lv;
231  }
232  }
233  if(anti>0) A[0][0] += kr2GNe;
234  else A[0][0] -= kr2GNe;
235 
236  // Solve Hamiltonian for eigensystem using the GLoBES method
237  zheevh3(A,fEvec,fEval);
238 
239 }
240 
241 ///.....................................................................
242 ///
243 /// Propagate the current neutrino state over a distance L in km
244 /// with an energy E in GeV through constant matter of density
245 /// Ne in mole/cm^3.
246 /// @param anti - +1 = neutrino case, -1 = anti-neutrino case
247 ///
248 template <typename T>
249 void _PMNSOpt<T>::PropMatter(double L, double E, double Ne, int anti)
250 {
251 
252  // Solve Hamiltonian
253  this->SolveHam(E, Ne, anti);
254 
255  // Store coefficients of propagation eigenstates
256  complex nuComp[3];
257 
258  for(int i=0;i<3;i++){
259  nuComp[i] = complex(0,0);
260  for(int j=0;j<3;j++){
261  nuComp[i] += fNuState[j] * conj(fEvec[j][i]);
262  }
263  }
264 
265  for(int i=0;i<3;i++)fNuState[i] = complex(0,0);
266 
267  // Propagate neutrino state
268  for(int j=0;j<3;j++){
269  auto s = sin(-fEval[j] * kKm2eV * L);
270  auto c = cos(-fEval[j] * kKm2eV * L);
271 
272  complex jPart = complex(c, s) * nuComp[j];
273  for(int i=0;i<3;i++){
274  fNuState[i] += jPart * fEvec[i][j];
275  }
276  }
277 
278 }
279 
280 //......................................................................
281 ///
282 /// Do several layers in a row. L and Ne must have the same length
283 ///
284 template <typename T>
285 void _PMNSOpt<T>::PropMatter(const std::list<double>& L,
286  double E,
287  const std::list<double>& Ne,
288  int anti)
289 {
290  if (L.size()!=Ne.size()) abort();
291  auto Li = L.begin();
292  auto Lend = L.end();
293  auto Ni = Ne.begin();
294  for (; Li!=Lend; ++Li, ++Ni) {
295  // For very low densities, use vacumm
296  static const double kRhoCutoff = 1.0E-6; // Ne in moles/cm^3
297  if (*Ni<kRhoCutoff) this->PropVacuum(*Li, E, anti);
298  else this->PropMatter(*Li, E, *Ni, anti);
299  }
300 }
301 
302 
303 ///.....................................................................
304 ///
305 /// We know the vacuum eigensystem, so just write it explicitly
306 /// The eigenvalues depend on energy, so E needs to be provided in GeV
307 /// @param anti - +1 = neutrino case, -1 = anti-neutrino case
308 ///
309 template <typename T>
310 void _PMNSOpt<T>::SetVacuumEigensystem(double E, int anti)
311 {
312 
313  T s12, s23, s13, c12, c23, c13;
314  complex expidelta(cos(fDeltaCP), anti * sin(fDeltaCP));
315 
316  s12 = sin(fTheta12); s23 = sin(fTheta23); s13 = sin(fTheta13);
317  c12 = cos(fTheta12); c23 = cos(fTheta23); c13 = cos(fTheta13);
318 
319  fEvec[0][0] = complex(c12*c13, 0);
320  fEvec[0][1] = complex(s12*c13, 0);
321  fEvec[0][2] = s13*conj(expidelta);
322 
323  fEvec[1][0] = -s12*c23-c12*s23*s13*expidelta;
324  fEvec[1][1] = c12*c23-s12*s23*s13*expidelta;
325  fEvec[1][2] = complex(s23*c13, 0);
326 
327  fEvec[2][0] = s12*s23-c12*c23*s13*expidelta;
328  fEvec[2][1] = -c12*s23-s12*c23*s13*expidelta;
329  fEvec[2][2] = complex(c23*c13, 0);
330 
331  fEval[0] = 0;
332  fEval[1] = fDm21 / (2 * kGeV2eV*E);
333  fEval[2] = fDm31 / (2 * kGeV2eV*E);
334 
335 }
336 
337 ///.....................................................................
338 ///
339 /// Propagate the current neutrino state over a distance L in km
340 /// with an energy E in GeV through vacuum
341 /// @param anti - +1 = neutrino case, -1 = anti-neutrino case
342 ///
343 template <typename T>
344 void _PMNSOpt<T>::PropVacuum(double L, double E, int anti)
345 {
346 
347  this->SetVacuumEigensystem(E, anti);
348 
349  complex nuComp[3];
350 
351  for(int i=0;i<3;i++){
352  nuComp[i] = 0;
353  for(int j=0;j<3;j++){
354  nuComp[i] += fNuState[j] * conj(fEvec[j][i]);
355  }
356  }
357 
358  const T km2EvL = kKm2eV*L; // needed for the templated multiplication below to work
359  for(int i=0;i<3;i++){
360  fNuState[i] = 0;
361  for(int j=0;j<3;j++){
362  complex iEval(0.0,fEval[j]);
363  fNuState[i] += exp(-iEval * km2EvL) * nuComp[j] * fEvec[i][j];
364  }
365  }
366 
367 }
368 
369 //......................................................................
370 ///
371 /// Reset the neutrino state back to a pure flavour where
372 /// it starts
373 ///
374 template <typename T>
376 {
377  int i;
378  for (i=0; i<3; ++i){
379  if (i==flv) fNuState[i] = complex(1, 0);
380  else fNuState[i] = complex(0, 0);
381  }
382 }
383 
384 //......................................................................
385 ///
386 /// Compute oscillation probability of flavour flv
387 ///
388 /// 0 = nue, 1 = numu, 2 = nutau
389 ///
390 template <typename T>
391 T _PMNSOpt<T>::P(int flv) const
392 {
393  assert(flv>=0 && flv<3);
394  return norm(fNuState[flv]);
395 }
396 
397 ////////////////////////////////////////////////////////////////////////
398 // manually instantiate templates for those cases we know about.
399 
400 template class osc::_PMNSOpt<double>;
401 
402 #ifdef OSCLIB_STAN
403 template class osc::_PMNSOpt<stan::math::var>;
404 #endif
virtual void BuildHlv()
Definition: PMNSOpt.cxx:147
virtual void PropVacuum(double L, double E, int anti=1)
Definition: PMNSOpt.cxx:344
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:375
virtual void SetVacuumEigensystem(double E, int anti)
Definition: PMNSOpt.cxx:310
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:104
const XML_Char * s
Definition: expat.h:262
virtual void SetDeltaMsqrs(const T &dm21, const T &dm32)
Definition: PMNSOpt.cxx:122
static const double kGeV2eV
GeV to eV.
Definition: PMNSOpt.h:43
virtual void SolveHam(double E, double Ne, int anti)
Definition: PMNSOpt.cxx:209
static constexpr double L
Float_t E
Definition: plot.C:20
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:249
virtual T P(int flv) const
Definition: PMNSOpt.cxx:391
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