OscCalcAnalytic.cxx
Go to the documentation of this file.
2 
3 #include <cassert>
4 
5 // Stan doesn't provide sincos(). Make sure not to hide the fast one for double
6 template<class T> std::enable_if_t<!std::is_arithmetic_v<T>, void> sincos(const T& x, T* sx, T* cx)
7 {
8  *sx = sin(x);
9  *cx = cos(x);
10 }
11 
12 namespace osc
13 {
14  template<class T> inline __attribute__((always_inline)) cmplx<T>
15  operator*(const cmplx<T>& x, const cmplx<T>& y)
16  {
17  return cmplx<T>(x.re*y.re - x.im*y.im, x.re*y.im + x.im*y.re);
18  }
19 
20  template<class T, class U> inline __attribute__((always_inline)) cmplx<T>
21  operator*(const U& x, const cmplx<T>& y)
22  {
23  return cmplx<T>(x*y.re, x*y.im);
24  }
25 
26  template<class T, class U> inline __attribute__((always_inline)) cmplx<T>
27  operator*(const cmplx<T>& x, const U& y)
28  {
29  return cmplx<T>(x.re*y, x.im*y);
30  }
31 
32  template<class T, class U> inline __attribute__((always_inline)) cmplx<T>
33  operator/(const cmplx<T>& x, const U& y)
34  {
35  return cmplx<T>(x.re/y, x.im/y);
36  }
37 
38  template<class T> inline __attribute__((always_inline)) cmplx<T>
39  operator+(const cmplx<T>& x, const cmplx<T>& y)
40  {
41  return cmplx<T>(x.re + y.re, x.im + y.im);
42  }
43 
44  template<class T, class U> inline __attribute__((always_inline)) cmplx<T>
45  operator+(const cmplx<T>& x, const U& y)
46  {
47  return cmplx<T>(x.re + y, x.im);
48  }
49 
50  template<class T, class U> inline __attribute__((always_inline)) cmplx<T>
51  operator+(const U& x, const cmplx<T>& y)
52  {
53  return cmplx<T>(x + y.re, y.im);
54  }
55 
56  template<class T> inline __attribute__((always_inline)) cmplx<T>
57  operator-(const cmplx<T>& x, const cmplx<T>& y)
58  {
59  return cmplx<T>(x.re - y.re, x.im - y.im);
60  }
61 
62  template<class T, class U> inline __attribute__((always_inline)) cmplx<T>
63  operator-(const cmplx<T>& x, const U& y)
64  {
65  return cmplx<T>(x.re - y, x.im);
66  }
67 
68  template<class T, class U> inline __attribute__((always_inline)) cmplx<T>
69  operator-(const U& x, const cmplx<T>& y)
70  {
71  return cmplx<T>(x - y.re, -y.im);
72  }
73 
74  template<class T> inline __attribute__((always_inline)) cmplx<T>
75  operator-(const cmplx<T>& x)
76  {
77  return cmplx<T>(-x.re, -x.im);
78  }
79 
80  //---------------------------------------------------------------------------
82  : fDirty12(true), fDirty13(true), fDirty23(true), fDirtyCP(true), fDirtyMasses(true),
83  Ue3(0, 0), Um2(0, 0), Ut2(0, 0),
84  Hem(0, 0), Het(0, 0), Hmt(0, 0),
85  Mem(0, 0), Met(0, 0), Mmt(0, 0)
86  {
87  }
88 
89  //---------------------------------------------------------------------------
91  {
92  }
93 
94  //---------------------------------------------------------------------------
96  Copy() const
97  {
98  return new _OscCalcAnalytic<T>(*this);
99  }
100 
101  //---------------------------------------------------------------------------
102  template<class T> void _OscCalcAnalytic<T>::SetL(double L)
103  {
104  if(L == this->fL) return;
105 
106  this->fL = L;
107  fProbCache.clear();
108  }
109 
110  //---------------------------------------------------------------------------
111  template<class T> void _OscCalcAnalytic<T>::SetRho(double rho)
112  {
113  if(rho == this->fRho) return;
114 
115  this->fRho = rho;
116  fProbCache.clear();
117  }
118 
119  //---------------------------------------------------------------------------
120  template<class T> void _OscCalcAnalytic<T>::SetDmsq21(const T& dmsq21)
121  {
122  if constexpr(std::is_arithmetic_v<T>) if(dmsq21 == this->fDmsq21) return;
123 
124  this->fDmsq21 = dmsq21;
125  fDirtyMasses = true;
126  }
127 
128  //---------------------------------------------------------------------------
129  template<class T> void _OscCalcAnalytic<T>::SetDmsq32(const T& dmsq32)
130  {
131  if constexpr(std::is_arithmetic_v<T>) if(dmsq32 == this->fDmsq32) return;
132 
133  this->fDmsq32 = dmsq32;
134  fDirtyMasses = true;
135  }
136 
137  //---------------------------------------------------------------------------
138  template<class T> void _OscCalcAnalytic<T>::SetTh23(const T& th23)
139  {
140  if constexpr(std::is_arithmetic_v<T>) if(th23 == this->fTh23) return;
141 
142  this->fTh23 = th23;
143  fDirty23 = true;
144  }
145 
146  //---------------------------------------------------------------------------
147  template<class T> void _OscCalcAnalytic<T>::SetTh13(const T& th13)
148  {
149  if constexpr(std::is_arithmetic_v<T>) if(th13 == this->fTh13) return;
150 
151  this->fTh13 = th13;
152  fDirty13 = true;
153  }
154 
155  //---------------------------------------------------------------------------
156  template<class T> void _OscCalcAnalytic<T>::SetTh12(const T& th12)
157  {
158  if constexpr(std::is_arithmetic_v<T>) if(th12 == this->fTh12) return;
159 
160  this->fTh12 = th12;
161  fDirty12 = true;
162  }
163 
164  //---------------------------------------------------------------------------
165  template<class T> void _OscCalcAnalytic<T>::SetdCP(const T& delta)
166  {
167  if constexpr(std::is_arithmetic_v<T>) if(delta == this->fdCP) return;
168 
169  this->fdCP = delta;
170  fDirtyCP = true;
171  }
172 
173  //---------------------------------------------------------------------------
174  template<class T> TMD5* _OscCalcAnalytic<T>::GetParamsHash() const
175  {
177  }
178 
179  //---------------------------------------------------------------------------
180  template<class T> double _OscCalcAnalytic<T>::Hmat()
181  {
182  // Need to convert avogadro's constant so that the total term comes out in
183  // units of inverse distance. Note that Ne will be specified in g/cm^-3
184  // I put the following into Wolfram Alpha:
185  // (fermi coupling constant)*((avogadro number)/cm^3)*(reduced planck constant)^2*(speed of light)^2
186  // And then multiplied by 1000 because we specify L in km not m.
187  const double GF = 1.368e-4;
188  const double Ne = this->fRho;
189  return sqrt(2)*GF*Ne;
190  }
191 
192  //---------------------------------------------------------------------------
193  template<class T> T sqr(T x){return x*x;}
194  template<class T> T cube(T x){return x*x*x;}
195 
196  const double sqrt3 = sqrt(3);
197 
198  //---------------------------------------------------------------------------
199  /// Solve x^3 + b*x^2 + c*x + d = 0
200  template<class T> std::array<T, 3> SolveCubic(T b, T c, T d)
201  {
202  b /= 3;
203  // Now solving x^3 + 3*b*x^2 + c*x + d = 0
204 
205  // https://en.wikipedia.org/wiki/Cubic_equation#Depressed_cubic
206  const T p = c/3 - sqr(b);
207  const T q = 2*cube(b) - b*c + d;
208  // Now solving t^3 + 3*p*t + q = 0
209 
210  // https://en.wikipedia.org/wiki/Cubic_equation#Trigonometric_solution_for_three_real_roots
211  //
212  // Cardano's formula looks simpler, but you still have to involve trig
213  // operations anyway for the cube root, and it's harder to keep everything
214  // explicitly real.
215 
216  const T s = sqrt(-p);
217 
218  const T r = acos(q/(2*p*s)) / 3;
219 
220  // t_k = 2*s * cos(r/3 + 2pi*k/3)
221 
222  // Use cos(a+b) = cosa*cosb - sina*sinb to save one trig operation
223  T sinr, cosr;
224  sincos(r, &sinr, &cosr);
225 
226  const T t0 = 2*s*cosr;
227  const T t1 = s*(sqrt3*sinr - cosr);
228  const T t2 = -t0-t1;
229 
230  return {t0 - b, t1 - b, t2 - b};
231  }
232 
233  //---------------------------------------------------------------------------
235  {
236  // const T a = 1; // cube term
237 
238  const T b = -Mee-Mmm-Mtt; // square term
239 
240  const T Lem = Mem.norm();
241  const T Let = Met.norm();
242  const T Lmt = Mmt.norm();
243 
244  // Matrix is Hermitian
245  const T c = Mee*Mmm + Mee*Mtt + Mmm*Mtt - Lem - Let - Lmt; // linear term
246 
247  const T d = (Mee*Lmt + Mmm*Let + Mtt*Lem
248  - Mee*Mmm*Mtt
249  -2*(Mem.re * Mmt.re * Met.re + Mem.re * Mmt.im * Met.im + Mem.im * Mmt.re * Met.im - Mem.im * Mmt.im * Met.re)); // const term
250 
251  const std::array<T, 3> xs = SolveCubic(b, c, d);
252 
253  // Overall phase doesn't matter, which allows us to save one exponentiation
254  T c10, s10, c20, s20;
255  sincos(xs[1]-xs[0], &s10, &c10);
256  sincos(xs[2]-xs[0], &s20, &c20);
257 
258  const T ei0 = 1/(3*sqr(xs[0]) + 2*b*xs[0] + c);
259  const cmplx<T> ei1 = cmplx(c10, s10)/(3*sqr(xs[1]) + 2*b*xs[1] + c);
260  const cmplx<T> ei2 = cmplx(c20, s20)/(3*sqr(xs[2]) + 2*b*xs[2] + c);
261 
262  // TODO, is there a way to more cheaply calculate these combinations?
263  return {
264  ei0 + ei1 + ei2,
265  xs[0]*ei0 + xs[1]*ei1 + xs[2]*ei2,
266  sqr(xs[0])*ei0 + sqr(xs[1])*ei1 + sqr(xs[2])*ei2};
267  }
268 
269  //---------------------------------------------------------------------------
270  template<class T> void _OscCalcAnalytic<T>::UpdatePMNS()
271  {
272  Ue2 = s12*c13;
273  Um3 = s23*c13;
274  Ut3 = c23*c13;
275 
276  const cmplx<T> phase(cCP, sCP);
277  Ue3 = s13*phase.conj();
278  Um2 = c12*c23-(s12*s23*s13)*phase;
279  Ut2 = -c12*s23-(s12*c23*s13)*phase;
280  }
281 
282  //---------------------------------------------------------------------------
283  template<class T> void _OscCalcAnalytic<T>::UpdateHamiltonian()
284  {
285  // d1 would be zero, so all those terms drop out
286  const T d2 = this->fDmsq21;
287  const T d3 = this->fDmsq21 + this->fDmsq32;
288 
289  Hee = d2 * sqr(Ue2) + d3 * Ue3.norm();
290  Hmm = d2 * Um2.norm() + d3 * sqr(Um3);
291  Htt = d2 * Ut2.norm() + d3 * sqr(Ut3);
292 
293  Hem = d2 * Ue2 * Um2.conj() + d3 * Ue3 * Um3;
294  Het = d2 * Ue2 * Ut2.conj() + d3 * Ue3 * Ut3;
295  Hmt = d2 * Um2 * Ut2.conj() + d3 * Um3 * Ut3;
296 
297  fProbCache.clear();
298  }
299 
300  //---------------------------------------------------------------------------
301  template<class T> T _OscCalcAnalytic<T>::Probs::P(int from, int to) const
302  {
303  // convert flavours to indices into matrix
304  const int i0 = (from-12)/2;
305  const int i1 = (to-12)/2;
306 
307  // Exploit unitarity
308  switch(i0*3+i1){
309  case 0: return Pee;
310  case 1: return Pme;
311  case 2: return 1-Pee-Pme; // Pte
312  case 3: return Pem;
313  case 4: return Pmm;
314  case 5: return 1-Pem-Pmm; // Ptm
315  case 6: return 1-Pee-Pem; // Pet
316  case 7: return 1-Pme-Pmm; // Pmt
317  case 8: return Pee+Pem+Pme+Pmm-1; // Ptt
318  default: abort();
319  }
320  }
321 
322  //---------------------------------------------------------------------------
323  template<class T> T _OscCalcAnalytic<T>::P(int from, int to, double E)
324  {
325  // -E effectively flips rho and conjugates H
326  if(from < 0) return P(-from, -to, -E);
327 
328  assert(from > 0 && to > 0);
329 
330  assert(from == 12 || from == 14 || from == 16);
331  assert(to == 12 || to == 14 || to == 16);
332 
333  const bool dirtyAngles = fDirty12 || fDirty13 || fDirty23 || fDirtyCP;
334 
335  if(dirtyAngles){
336  if(fDirty12) sincos(this->fTh12, &s12, &c12);
337  if(fDirty13) sincos(this->fTh13, &s13, &c13);
338  if(fDirty23) sincos(this->fTh23, &s23, &c23);
339  if(fDirtyCP) sincos(this->fdCP, &sCP, &cCP);
340  UpdatePMNS();
341  UpdateHamiltonian();
342  }
343  else{
344  if(fDirtyMasses){
345  UpdateHamiltonian();
346  }
347  else{
348  auto it = fProbCache.find(E);
349  if(it != fProbCache.end()) return it->second.P(from, to);
350  }
351  }
352 
354 
355  const double k = -this->fL * 2*1.267 / E;
356  Mee = Hee * k - this->fL * Hmat();
357  Mem = Hem * k;
358  Mmm = Hmm * k;
359  Met = Het * k;
360  Mmt = Hmt * k;
361  Mtt = Htt * k;
362 
363  // Matrix exponent is based on https://www.wolframalpha.com/input/?i=matrixExp+%5B%5Br%2Cs%2Ct%5D%2C%5Bu%2Cv%2Cw%5D%2C%5Bx%2Cy%2Cz%5D%5D
364 
365  const Eigenvalues es = GetEigenvalues(E);
366 
367  const T Aee = Mmm*Mtt - Mmt.norm();
368  const T Amm = Mee*Mtt - Met.norm();
369  const cmplx<T> Aem = Met*Mmt.conj() - Mem*Mtt;
370 
371  const Probs ps((Aee *es.sume - (Mmm+Mtt) *es.sumxe + es.sumxxe).norm(),
372  (Aem.conj()*es.sume + Mem.conj()*es.sumxe ).norm(),
373  (Aem *es.sume + Mem *es.sumxe ).norm(),
374  (Amm *es.sume - (Mee+Mtt) *es.sumxe + es.sumxxe).norm());
375 
376  fProbCache.emplace(E, ps);
377 
378  return ps.P(from, to);
379  }
380 
381 } // namespace
382 
383 
384 // Instantiate
385 template class osc::_OscCalcAnalytic<double>;
386 
387 #ifdef OSCLIB_STAN
388 #include "stan/math/rev/scal.hpp"
390 #endif
virtual TMD5 * GetParamsHash() const override
Use to check two calculators are in the same state.
virtual void SetRho(double rho) override
TDC::value_type operator-(TDC lhs, TDC rhs)
Definition: BaseProducts.h:47
virtual void SetL(double L) override
__attribute__((always_inline)) cmplx< T > operator*(const cmplx< T > &x
set< int >::iterator it
virtual void SetdCP(const T &dCP) override
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
HepLorentzVector operator/(const HepLorentzVector &, double a)
TDC::value_type operator+(TDC lhs, TDC rhs)
Definition: BaseProducts.h:49
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual void SetTh13(const T &th13) override
virtual void SetDmsq21(const T &dmsq21) override
double th12
Definition: runWimpSim.h:98
virtual T P(int from, int to, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
T acos(T number)
Definition: d0nt_math.hpp:54
virtual void SetTh23(const T &th23) override
TMD5 * GetParamsHashDefault(const std::string &txt) const
This is only a safe implementation if your calculator only depends on these eight parameters...
Definition: IOscCalc.cxx:44
#define P(a, b, c, d, e, x)
virtual _IOscCalcAdjustable< T > * Copy() const override
virtual void SetDmsq32(const T &dmsq32) override
const XML_Char * s
Definition: expat.h:262
float_mat operator*(const float_mat &a, const float_mat &b)
matrix multiplication.
T cube(T x)
static constexpr double L
Float_t E
Definition: plot.C:20
std::complex takes a lot of care with inf/nan which we don&#39;t want
T sqr(T x)
Float_t d
Definition: plot.C:236
double t2
Oscillation probability calculators.
Definition: Calcs.h:5
const double sqrt3
T sin(T number)
Definition: d0nt_math.hpp:132
Float_t norm
double dmsq32
virtual void SetTh12(const T &th12) override
General interface to any calculator that lets you set the parameters.
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
double T
Definition: Xdiff_gwt.C:5
const cmplx< T > & y
std::enable_if_t<!std::is_arithmetic_v< T >, void > sincos(const T &x, T *sx, T *cx)
std::unordered_map< double, Probs > fProbCache
double th13
Definition: runWimpSim.h:98
std::array< T, 3 > SolveCubic(T b, T c, T d)
Solve x^3 + b*x^2 + c*x + d = 0.