PMNS_DMP.cxx
Go to the documentation of this file.
1 // std_isnan needs to precede the header
2 // because it injects an important Stan overload of std::isnan
3 // into the std namespace that needs to be there
4 // BEFORE the Eigen headers are seen (due to the '#pragma once').
5 // PMNS_DMP.h #includes Eigen/Eigen, so /shrug
6 #ifdef OSCLIB_STAN
8 #endif
9 
10 #include "OscLib/PMNS_DMP.h"
11 
12 namespace osc
13 {
14 
15 
16  template <typename T>
18 
19  _s212 = sin(2 * params.th12);
20  _s213 = sin(2 * params.th13);
21  _c212 = cos(2 * params.th12);
22  _c213 = cos(2 * params.th13);
23 
24  _Dmsq21 = params.dmsq21;
25  _Dmsq31 = params.Dmsq31();
26  _cd = cos(params.deltacp);
27  _sd = sin(params.deltacp);
28 
29  T s12 = sin(params.th12);
30  _Dmsqee = _Dmsq21 + params.dmsq32 - s12*s12 * _Dmsq21;
31 
32  _s13 = sin(params.th13);
33  _s23 = sin(params.th23);
34 
35  _s13sq = _s13*_s13;
36  _s23sq = _s23*_s23;
37 
38  _c13sq = 1 - _s13sq;
39  _c23sq = 1 - _s23sq;
40 
41  _c13 = sqrt(_c13sq);
42  _c23 = sqrt(_c23sq);
43 
44  }
45 
46  template <typename T>
48  _OscParameters<T> p = {params[0], params[1], params[2], params[3], params[4], params[5]};
49  this->SetParameters(p);
50  }
51 
52  template <typename T>
54  _OscParameters<T> p = {params[0], params[1], params[2], params[3], params[4], params[5]};
55  return this->P(p);
56  }
57 
58  template <typename T>
60  _OscParameters<T> p = { dmsq21, dmsq32, th12, th13, th23, delta};
61  return this->P(p);
62  }
63 
64  template <typename T>
66 
67  this->SetParameters(params);
68 
69  const ParamArray DMSQEEA = _Dmsqee * sqrt(square(_c213*_ONE - 1./_Dmsqee * _AA) + _s213*_s213*_ONE);
70 
71  const ParamArray C2PHI = (_Dmsqee * _c213 * _ONE - _AA) / DMSQEEA;
72  const ParamArray A12 = 0.5 * (_AA + _Dmsqee*_ONE - DMSQEEA);
73 
74  const ParamArray CPHISQ = 0.5 * (_ONE + C2PHI);
75  const ParamArray SPHISQ = 0.5 * (_ONE - C2PHI);
76  const ParamArray S2PHI = _s213 * _Dmsqee * _ONE / DMSQEEA;
77  const ParamArray CPHI13SQ = CPHISQ * _c13sq + SPHISQ * _s13sq + S2PHI * _c13 * _s13;
78 
79  const ParamArray DL21 = _Dmsq21 * sqrt(square(_c212*_ONE - A12 / _Dmsq21) + CPHI13SQ * _s212 * _s212);
80  const ParamArray C2PSI = (_Dmsq21 * _c212 * _ONE - A12) / DL21;
81  const ParamArray CPSISQ = 0.5 * (_ONE + C2PSI);
82  const ParamArray SPSISQ = 0.5 * (_ONE - C2PSI);
83 
84  const ParamArray DL31 = _Dmsq31*_ONE + 0.25 * _AA + 0.5 * (DL21 - _Dmsq21*_ONE) + 0.75 * (DMSQEEA - _Dmsqee*_ONE);
85 
86  const ParamArray JRRM = _s23 * _c23 * sqrt(SPHISQ * CPSISQ * SPSISQ);
87  const ParamArray JRM = JRRM * CPHISQ;
88  const ParamArray D = -JRM * _sd;
89 
90  ProbArray DD(_ENERGIES.rows(), 3);
91  DD.col(0) = DL21 * _L4E;
92  DD.col(1) = DL31 * _L4E;
93  DD.col(2) = DD.col(1)-DD.col(0);
94  const ProbArray SDD = DD.sin();
95 
96  ProbArray PP(_ENERGIES.rows(), 9);
97 
98  // e -> e
99  ParamArray C31 = -1 * CPHISQ * SPHISQ * CPSISQ;
100  ParamArray C32 = -1 * CPHISQ * SPHISQ * SPSISQ;
101  ParamArray C21 = -1 * CPHISQ.square() * SPSISQ * CPSISQ;
102  PP.col(0) = _ONE + 4 * (C31 * SDD.col(1).square() + C32 * SDD.col(2).square() + C21 * SDD.col(0).square());
103 
104  // mu -> mu
105  C31 = -1 * CPHISQ * _s23sq * (_c23sq * SPSISQ + _s23sq * SPHISQ * CPSISQ) \
106  - 2 * _s23sq * JRM * _cd;
107  C32 = -1 * CPHISQ * _s23sq * (_c23sq * CPSISQ + _s23sq * SPHISQ * SPSISQ) \
108  + 2 * _s23sq * JRM * _cd;
109  C21 = -1 * (_c23sq * CPSISQ + _s23sq * SPHISQ * SPSISQ) * (_c23sq * SPSISQ + _s23sq * SPHISQ * CPSISQ) \
110  - 2 * (_c23sq - SPHISQ * _s23sq) * C2PSI * JRRM * _cd + square(2 * JRRM * _cd);
111  PP.col(4) = _ONE + 4 * (C31 * SDD.col(1).square() + C32 * SDD.col(2).square() + C21 * SDD.col(0).square());
112 
113  C31 = _s23sq * SPHISQ * CPHISQ * CPSISQ + JRM * _cd;
114  C32 = _s23sq * SPHISQ * CPHISQ * SPSISQ - JRM * _cd;
115  C21 = CPHISQ * SPSISQ * CPSISQ * (_c23sq* _ONE - SPHISQ * _s23sq) + JRM * _cd * C2PSI;
116  // Pemu
117  PP.col(3) = 4 * (C31 * SDD.col(1).square() + C32 * SDD.col(2).square() + C21 * SDD.col(0).square()) - 8 * D * SDD.col(0) * SDD.col(1) * SDD.col(2);
118  // Pmue
119  PP.col(1) = 4 * (C31 * SDD.col(1).square() + C32 * SDD.col(2).square() + C21 * SDD.col(0).square()) + 8 * D * SDD.col(0) * SDD.col(1) * SDD.col(2);
120 
121  // e->tau = 1 - ee - emu
122  PP.col(6) = _ONE - PP.col(0) - PP.col(3);
123  // tau->e = 1 - ee - mue
124  PP.col(2) = _ONE - PP.col(0) - PP.col(1);
125  // mu->tau = 1 - mue - mumu
126  PP.col(7) = _ONE - PP.col(1) - PP.col(4);
127  // tau->mu = 1 - emu - mumu
128  PP.col(5) = _ONE - PP.col(3) - PP.col(4);
129  // tau->tau = 1 - etau - mutau
130  PP.col(8) = _ONE - PP.col(6) - PP.col(7);
131 
132  return PP.isNaN().select(0,PP);
133  }
134 }
135 
136 
137 ////////////////////////////////////////////////////////////////////////
138 // manually instantiate templates for those cases we know about.
139 
140 template class osc::_PMNS_DMP<double>;
141 
142 #ifdef OSCLIB_STAN
143 #include "stan/math/rev/scal.hpp"
144 template class osc::_PMNS_DMP<stan::math::var>;
145 #endif
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
ProbArray P(_OscParameters< T > const &params)
Definition: PMNS_DMP.cxx:65
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SetParameters(_OscParameters< T > const &params)
Definition: PMNS_DMP.cxx:17
double th12
Definition: runWimpSim.h:98
constexpr T square(T x)
Definition: pow.h:21
Array< T, Dynamic, Dynamic > ProbArray
Definition: PMNS_DMP.h:26
Array< T, Dynamic, 1 > ParamArray
Definition: PMNS_DMP.h:25
#define P(a, b, c, d, e, x)
Oscillation probability calculators.
Definition: Calcs.h:5
T sin(T number)
Definition: d0nt_math.hpp:132
double dmsq32
T cos(T number)
Definition: d0nt_math.hpp:78
double T
Definition: Xdiff_gwt.C:5
double th13
Definition: runWimpSim.h:98