OscCalcDMP.cxx
Go to the documentation of this file.
1 #include "OscLib/OscCalcDMP.h"
2 #include "OscLib/PMNS_DMP.h"
3 
4 #include <cmath>
5 
6 namespace osc
7 {
8  template<typename T>
9  TMD5*
11  {
12  std::string txt = "DMP";
13  TMD5* ret = new TMD5;
14  ret->Update((unsigned char*)txt.c_str(), txt.size());
15  const int kNumParams = 8;
16  T buf[kNumParams] = {this->fRho, this->fL, this->fDmsq21, this->fDmsq32,
17  this->fTh12, this->fTh13, this->fTh23, this->fdCP};
18  ret->Update((unsigned char*)buf, sizeof(T)*kNumParams);
19  ret->Final();
20  return ret;
21  }
22 
23 
24  template<typename T>
26  {
27  this->fLastParams.L = this->fL;
28  this->fLastParams.rho = this->fRho;
29  this->fLastParams.dmsq21 = this->fDmsq21;
30  this->fLastParams.dmsq32 = this->fDmsq32;
31  this->fLastParams.th12 = this->fTh12;
32  this->fLastParams.th13 = this->fTh13;
33  this->fLastParams.th23 = this->fTh23;
34  this->fLastParams.deltacp = this->fdCP;
35  }
36 
37  inline void _sincos(double theta, double &s, double &c)
38  {
39 #ifdef __APPLE_CC__
40  __sincos(theta, &s, &c);
41 #else
42  sincos(theta, &s, &c);
43 #endif
44  }
45 
46  template<typename T>
49  {
50  return new _OscCalcDMP<T>(*this);
51  }
52 
53  template<typename T>
55  {
56  return this->fDmsq21 == this->fCache.parameters.dmsq21 &&
57  this->fDmsq32 == this->fCache.parameters.dmsq32 &&
58  this->fTh12 == this->fCache.parameters.th12 &&
59  this->fTh13 == this->fCache.parameters.th13 &&
60  this->fTh23 == this->fCache.parameters.th23 &&
61  this->fdCP == this->fCache.parameters.deltacp &&
62  this->fL == this->fCache.parameters.L &&
63  this->fRho == this->fCache.parameters.rho;
64  }
65 
66  template<typename T>
67  inline int _OscCalcDMP<T>::ChannelCacheIdx(int flavBefore, int flavAfter) const
68  {
69  // rows in the cache are arranged in the following order
70  // 11 21 31 12 22 32 13 23 33 -11 -21 -31 -12 -22 -32 -13 -23 -33
71  // where nue = 1
72  // numu = 2
73  // nutau = 3
74  // and negative is for antineutrino
75  int anti = flavBefore / abs(flavBefore);
76  int i = (abs(flavBefore) - 12) / 2;
77  int j = (abs(flavAfter) - 12) / 2;
78  int idx = (1 - anti) / 2 * 9 + (3 * j + i);
79  return idx;
80  }
81 
82 
83  template<typename T>
84  Array<T, Dynamic, 1> _OscCalcDMP<T>::P(int flavBefore, int flavAfter, const std::vector<double> &E)
85  {
86  if (fCache.energies.size() != (size_t) fCache.probabilities.cols() &&
87  fCache.energies.size() != 0)
88  { // does a cache exist
89  if (ParamsAreCached())
90  {
91  if (this->fCache.energies == E)
92  return this->fCache.probabilities.col(ChannelCacheIdx(flavBefore, flavAfter));
93  }
94  }
95  FillCache(E);
96  return this->fCache.probabilities.col(ChannelCacheIdx(flavBefore, flavAfter));
97  }
98 
99  template<typename T>
100  T _OscCalcDMP<T>::P(int flavBefore, int flavAfter, double E,
101  bool fast_and_loose)
102  {
103  if (fast_and_loose)
104  {
105  auto e_it = std::find(fCache.energies.begin(),
106  fCache.energies.end(),
107  E);
108  return fCache.probabilities(e_it - fCache.energies.begin(),
109  ChannelCacheIdx(flavBefore, flavAfter));
110  } else
111  {
112  return P(flavBefore, flavAfter, E);
113  }
114  }
115 
116  template<typename T>
117  T _OscCalcDMP<T>::P(int flavBefore, int flavAfter, double E)
118  {
119  if(std::is_arithmetic_v<T> &&
120  fCache.energies.size() != (size_t) fCache.probabilities.cols() &&
121  fCache.energies.size() != 0)
122  { // does a cache exist
123  if (ParamsAreCached())
124  { // do current params match those cached
125  auto e_it = std::find(fCache.energies.begin(),
126  fCache.energies.end(),
127  E);
128  if (e_it != fCache.energies.end())
129  { // is the given energy cached?
130 
131  return this->fCache.probabilities(e_it - fCache.energies.begin(),
132  ChannelCacheIdx(flavBefore, flavAfter));
133 
134  }
135  }
136  }
137  // if we get here, the cache is stale.
138  // TODO. write function ValidateCache
139  this->FillCache({E});
140  return this->fCache.probabilities(0, ChannelCacheIdx(flavBefore, flavAfter));
141  }
142 
143 
144  template<typename T>
146  {
147  using ArrayXXT = Array<T, Dynamic, Dynamic>;
148 
149  const _OscParameters<T> params{this->fDmsq21, this->fDmsq32, this->fTh12, this->fTh13, this->fTh23, this->fdCP,
150  this->fL, this->fRho};
151  _PMNS_DMP<T> OO(this->fCache.ENERGIES, this->fRho, this->fL);
152  ArrayXXT _cache = OO.P(params);
153 
154  int nbins = this->fCache.energies.size();
155  ArrayXXT cache(nbins, 18);
156  cache.block(0, 0, nbins, 9) = _cache.block(0, 0, nbins, 9);
157  cache.block(0, 9, nbins, 9) = _cache.block(nbins, 0, nbins, 9);
158 
159  this->fCache.probabilities = cache;
160  this->fCache.parameters = params;
161  this->fCache.iter++;
162  }
163 
164  template<typename T>
165  void _OscCalcDMP<T>::FillCache(std::vector<double> const &energies)
166  {
167  this->SetCachedEnergies(energies);
168  this->FillCache();
169  }
170 
171  template<typename T>
172  void _OscCalcDMP<T>::SetCachedEnergies(std::vector<double> const &energies)
173  {
174  this->fCache.energies = energies;
175  int nbins = energies.size();
176  ArrayXd EE(2 * nbins);
177  EE.segment(0, nbins) = Map<const ArrayXd>(energies.data(), energies.size());
178  EE.segment(nbins, nbins) = -1 * Map<const ArrayXd>(energies.data(), energies.size());
179  this->fCache.ENERGIES = EE;
180  }
181 
182 } // namespace
183 
184 //---------------------------------------------------------------------------
185 template class osc::_OscCalcDMP<double>;
186 
187 #ifdef OSCLIB_STAN
188 #include "stan/math/rev/scal.hpp"
189 template class osc::_OscCalcDMP<stan::math::var>;
190 #endif
Array< T, Dynamic, Dynamic > P()
Definition: OscCalcDMP.h:49
void _sincos(double theta, double &s, double &c)
Definition: OscCalcDMP.cxx:37
ProbArray P(_OscParameters< T > const &params)
Definition: PMNS_DMP.cxx:65
Helper struct for the cache. Might not need this.
Definition: StanTypedefs.h:25
void abs(TH1 *hist)
void sincos(T &x, Eigen::ArrayX< U > *sx, Eigen::ArrayX< U > *cx)
virtual void FillCache()
Definition: OscCalcDMP.cxx:145
#define P(a, b, c, d, e, x)
const XML_Char * s
Definition: expat.h:262
const int nbins
Definition: cellShifts.C:15
TMD5 * GetParamsHash() const override
Use to check two calculators are in the same state.
Definition: OscCalcDMP.cxx:10
bool ParamsAreCached()
Definition: OscCalcDMP.cxx:54
Float_t E
Definition: plot.C:20
const double j
Definition: BetheBloch.cxx:29
_IOscCalcAdjustable< T > * Copy() const override
Definition: OscCalcDMP.cxx:48
Oscillation probability calculators.
Definition: Calcs.h:5
void SetCachedEnergies(std::vector< double > const &energies)
Definition: OscCalcDMP.cxx:172
General interface to any calculator that lets you set the parameters.
void SaveLastParams()
Definition: OscCalcDMP.cxx:25
double T
Definition: Xdiff_gwt.C:5
int ChannelCacheIdx(int flavBefore, int flavAfter) const
Definition: OscCalcDMP.cxx:67
enum BeamMode string