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