OscCalcPMNSOptEigen.cxx
Go to the documentation of this file.
2 
3 #include <cmath>
4 
5 //#include "OscLib/PMNSOpt.cxx"
6 
7 //#include "zhetrd3.cxx"
8 //#include "zheevc3.cxx"
9 //#include "zheevh3.cxx"
10 //#include "zheevq3.cxx"
11 
12 namespace osc {
13 
14  TMD5*
16  {
17  std::string txt = "PMNSOptEigen";
18  TMD5* ret = new TMD5;
19  ret->Update((unsigned char*)txt.c_str(), txt.size());
20  const int kNumParams = 8;
21  double buf[kNumParams] = {fRho, fL, fDmsq21, fDmsq32,
22  fTh12, fTh13, fTh23, fdCP};
23  ret->Update((unsigned char*)buf, sizeof(double)*kNumParams);
24  ret->Final();
25  return ret;
26  }
27 
28  inline Eigen::Array33d
29  Probabilities(const Eigen::Matrix3cd & to,
30  const Eigen::Matrix3cd & from)
31  {
32  return (to * from.adjoint()).array().abs2();
33  }
34 
35  void
37  {
38  fLastParams.L = fL;
46  }
47 
48  void _sincos(double theta, double & s, double &c)
49  {
50 #ifdef __APPLE_CC__
51  __sincos(theta, &s, &c);
52 #else
53  sincos(theta, &s, &c);
54 #endif
55  }
56 
59  {
61  ret->fCache.clear();
62  return ret;
63  }
64 
65  bool
67  {
68  return
69  fDmsq21 == this->fCache.parameters.dmsq21 &&
70  fDmsq32 == this->fCache.parameters.dmsq32 &&
71  fTh12 == this->fCache.parameters.th12 &&
72  fTh13 == this->fCache.parameters.th13 &&
73  fTh23 == this->fCache.parameters.th23 &&
74  fdCP == this->fCache.parameters.deltacp &&
75  fL == this->fCache.parameters.L &&
76  fRho == this->fCache.parameters.rho;
77  }
78 
79  inline int
80  OscCalcPMNSOptEigen::ChannelCacheIdx(int flavBefore, int flavAfter) const
81  {
82  // rows in the cache are arranged in the following order
83  // 11 21 31 12 22 32 13 23 33 -11 -21 -31 -12 -22 -32 -13 -23 -33
84  // where nue = 1
85  // numu = 2
86  // nutau = 3
87  // and negative is for antineutrino
88  int anti = flavBefore / abs(flavBefore);
89  int i = (abs(flavBefore) - 12) / 2;
90  int j = (abs(flavAfter) - 12) / 2;
91  int idx = (1-anti)/2*9 + (3 * j + i);
92  return idx;
93  }
94 
95 
96  Eigen::VectorXd
97  OscCalcPMNSOptEigen::P(int flavBefore, int flavAfter, const std::vector<double>& E)
98  {
99  // Are there probabilities cached and can we use them?
100  if(fCache.energies.size() != (size_t) fCache.probabilities.cols() &&
101  fCache.energies.size() != 0) { // does a cache exist
102  if(ParamsAreCached()) {
103  if(this->fCache.energies == E)
104  return this->fCache.probabilities.col(ChannelCacheIdx(flavBefore, flavAfter));
105  }
106  }
107  FillCache(E);
108  return this->fCache.probabilities.col(ChannelCacheIdx(flavBefore, flavAfter));
109  }
110 
111  double
112  OscCalcPMNSOptEigen::P(int flavBefore, int flavAfter, double E,
113  bool fast_and_loose)
114  {
115  if(fast_and_loose) {
116  auto e_it = std::find(fCache.energies.begin(),
117  fCache.energies.end(),
118  E);
119  return fCache.probabilities(e_it - fCache.energies.begin(),
120  ChannelCacheIdx(flavBefore, flavAfter));
121  }
122  else {
123  return P(flavBefore, flavAfter, E);
124  }
125  }
126 
127  double
128  OscCalcPMNSOptEigen::P(int flavBefore, int flavAfter, double E)
129  {
130  // Are there probabilities cached and can we use them?
131  if(fCache.energies.size() != (size_t) fCache.probabilities.cols() &&
132  fCache.energies.size() != 0) { // does a cache exist
133  if(ParamsAreCached()) { // do current params match those cached
134  auto e_it = std::find(fCache.energies.begin(),
135  fCache.energies.end(),
136  E);
137  if(e_it != fCache.energies.end()) { // is the given energy cached?
138 
139  return this->fCache.probabilities(e_it - fCache.energies.begin(),
140  ChannelCacheIdx(flavBefore, flavAfter));
141 
142  }
143  }
144  }
145 
146 
147  // If we make it through this logic, just calculate it from scratch
148  const OscParameters params{fDmsq21, fDmsq32,
149  fTh12, fTh13, fTh23,
150  fdCP, fL, fRho};
151 
152  auto anti = 1;
153  if(flavBefore < 0) anti = -1;
154  EigenSystem flavor_states{Solve(MatterHamiltonian(E, anti, params))};
155 
156 
157  // returns Eigen Matrix of probabilities for each transition
158  auto probabilities = Probabilities(PropMatter(flavor_states.vectors,
159  flavor_states.values,
160  params),
161  flavor_states.vectors);
162  return probabilities((abs(flavBefore)-12)/2, (abs(flavAfter)-12)/2);
163 
164  /*
165  auto row = Eigen::Map<const Eigen::RowVectorXd>(probabilities.data(),
166  probabilities.size());
167  return row(ChannelCacheIdx(flavBefore, flavAfter));
168  */
169  }
170 
171 
172 
173  void
175  {
176  // std::cout << "Using OscCalcPMNSOptEigen\n";
177  // std::cout << "Filling cache" << std::endl;
178  const OscParameters params{fDmsq21, fDmsq32,
179  fTh12, fTh13, fTh23,
180  fdCP, fL, fRho};
181  /// probabilities returned are in columns of energy and rows of
182  /// flavors: 0-2 [nue, numu, nutau]
183  /// : 3-5 [antinue, antinumu, antinutau]
184  /// Dimension two is the "from" flavor
185  /// Dimension three is the "to" flavor
186  /// ie. P(numu->nue) is probabilities[energy_idx][0][1]
187  /// P(antinumu-> antinutau) is probabilities[energy_idx][4][5]
188  Eigen::ArrayXXd cache = Eigen::ArrayXXd::Zero(this->fCache.energies.size(), 18);
189 
190  Eigen::Matrix3cd const ham = BuildHam(params);
191  // matter hamiltonian is energy and flavor dependent
192  for(int anti : {-1, 1}) {
193 
194  //#pragma omp parallel for
195  for(unsigned int i = 0; i < this->fCache.energies.size(); i++) {
196  if(fCache.energies[i] <= 0) continue;
197 
198  // solve for flavor eigenstates. This part is expensive
199  EigenSystem const flavor_states =
200  Solve(AddMatterEffects(ham,
201  this->fCache.energies[i],
202  anti,
203  params));
204 
205  // returns Eigen Matrix of probabilities for each transition
206  auto const probabilities = Probabilities(PropMatter(flavor_states.vectors,
207  flavor_states.values,
208  params),
209  flavor_states.vectors);
210 
211  // reshape to a row vector
212  auto const row = Eigen::Map<const Eigen::Matrix<double, 1, 9> >(probabilities.data(), probabilities.size());
213 
214  // insert the matrix of probabilties into cache as a row
215  cache.block(i, (1-anti)/2 * 9, 1, 9) = row;// Eigen::Map<const Eigen::Matrix<double, 1, 9> >(probabilities.data(), probabilities.size());// row;
216  }
217  }
218 
219  this->fCache.probabilities = cache;
220  this->fCache.parameters = params;
221  this->fCache.iter++;
222 
223 // PrintMatrixAddresses(fCache.probabilities.data(),
224 // fCache.probabilities.rows(),
225 // fCache.probabilities.cols());
226 // PrintMatrixAddresses(this->fCache.probabilities);
227  }
228 
229  void
230  OscCalcPMNSOptEigen::PrintMatrixAddresses(const Eigen::ArrayXXd & mat) const
231  {
232 
233  for(int j = 0; j < mat.cols(); j++) {
234  for(int i = 0; i < mat.rows(); i++) {
235  std::cout << &mat(i, j) << " ";
236  }
237  std::cout << std::endl;
238  }
239  }
240 
241  void
243  int size) const
244  {
245  for(int i = 0; i < size; i++)
246  std::cout << data + i << ", (" << data[i] << ") ";
247  std::cout << std::endl;
248  }
249 
250  void
252  int nrows, int ncols) const
253  {
254  for(int i = 0; i < nrows; i++) {
255  for(int j = 0; j < ncols; j++) {
256  std::cout << data + i * ncols + j << " ";
257  }
258  std::cout << std::endl;
259  }
260  }
261 
262  Eigen::Matrix3cd
263  OscCalcPMNSOptEigen::PropMatter(Eigen::Matrix3cd const & evec,
264  Eigen::Vector3d const & evals,
265  const OscParameters & params) const
266  {
267  const Eigen::Array3d temp = evals.array()*kKm2eV*(-1)*params.L;
268  //Eigen::Array3d temp(evals.array());
269 
270  //// propagation phase
271  //temp *= (-1) * kKm2eV * params.L; // NOTE
272  const Eigen::Array3d SIN = temp.sin();
273  const Eigen::Array3d COS = temp.cos();
274  Eigen::Vector3cd PPT;
275  PPT <<
276  //std::complex<double>(std::cos(temp[0]), std::sin(temp[0])),
277  //std::complex<double>(std::cos(temp[1]), std::sin(temp[1])),
278  //std::complex<double>(std::cos(temp[2]), std::sin(temp[2]));
279  std::complex<double>(COS[0], SIN[0]),
280  std::complex<double>(COS[1], SIN[1]),
281  std::complex<double>(COS[2], SIN[2]);
282 
283  // propagated
284  return evec * PPT.asDiagonal();
285  }
286 
287 
288  void
289  OscCalcPMNSOptEigen::FillCache(std::vector<double> const &energies)
290  {
291  this->SetCachedEnergies(energies);
292  // std::cout << "Set Energies " << std::endl;
293  this->FillCache();
294  }
295 
296  void
297  OscCalcPMNSOptEigen::SetCachedEnergies(std::vector<double> const & energies)
298  {
299  this->fCache.energies = energies;
300  }
301 
302 
303  // Eigen::SelfAdjointEigenSolver<Eigen::Matrix3cd>
304  /*
305  EigenSystem
306  OscCalcPMNSOptEigen::SolveZheevh3(const Eigen::Matrix3cd & ham,
307  bool use_eigen)
308  {
309  std::complex<double> A[3][3];
310  std::complex<double> fEvec[3][3];
311  std::complex<double> fEval[3];
312  for(int i = 0; i < 3; i++) {
313  for(int j = 0; j < 3; j++) {
314  A[i][j] = ham(j, i);
315  }
316  }
317  zheevh3(A, fEvec, fEval, use_eigen);
318 
319  Matrix3cd vectors(3,3);
320  Array3d values(3);
321  vectors <<
322  fEvec[0][0], fEvec[1][0], fEvec[2][0],
323  fEvec[0][1], fEvec[1][1], fEvec[2][1],
324  fEvec[0][2], fEvec[1][2], fEvec[2][2];
325  values << fEval[0], fEval[1], fEval[2];
326 
327 
328  return EigenSystem{vectors, values};
329 
330  }
331  */
333  OscCalcPMNSOptEigen::Solve(Eigen::Matrix3cd const & ham) const
334  {
335  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3cd> eig;
336  eig.computeDirect(ham);
337  return EigenSystem{eig.eigenvectors(), eig.eigenvalues()};
338  }
339 
340  Eigen::Matrix3cd
341  OscCalcPMNSOptEigen::AddMatterEffects(const Eigen::Matrix3cd & ham,
342  const double & E,
343  const int & anti,
344  const OscParameters & params) const
345  {
346  double lv = 2 * kGeV2eV*E / params.Dmsq31(); // Osc. length in eV^-1
347  double kr2GNe = kK2*M_SQRT2*kGf * params.rho/2; // Matter potential in eV
348 
349  Eigen::Matrix3cd green_eggs(ham);
350  green_eggs.triangularView<Eigen::Upper>()/=lv;
351 
352  if (anti>0) {
353  // FIXME we only transposehere to be consistent with the original implementation
354  // --- resulting in all that mangling of eigenvectors later
355  green_eggs.transposeInPlace();
356  green_eggs(0,0) += kr2GNe;
357  }
358  else {
359  green_eggs.adjointInPlace();
360  green_eggs(0,0) -= kr2GNe;
361  }
362  return green_eggs;
363  }
364 
365  Eigen::Matrix3cd
366  OscCalcPMNSOptEigen::AddMatterEffectsAnti(const Eigen::Matrix3cd & ham,
367  const double & E,
368  const OscParameters & params) const
369  {
370  double lv = 2 * kGeV2eV*E / params.Dmsq31(); // Osc. length in eV^-1
371  double kr2GNe = kK2*M_SQRT2*kGf * params.rho/2; // Matter potential in eV
372 
373  Eigen::Matrix3cd green_eggs(ham);
374  green_eggs.triangularView<Eigen::Upper>()/=lv;
375 
376  green_eggs.adjointInPlace();
377  green_eggs(0,0) -= kr2GNe;
378 
379  return green_eggs;
380  }
381 
382 
383  Eigen::Matrix3cd
384  OscCalcPMNSOptEigen::AddMatterEffects(const Eigen::Matrix3cd & ham,
385  const double & E,
386  const OscParameters & params) const
387  {
388  double lv = 2 * kGeV2eV*E / params.Dmsq31(); // Osc. length in eV^-1
389  double kr2GNe = kK2*M_SQRT2*kGf * params.rho/2; // Matter potential in eV
390 
391  Eigen::Matrix3cd green_eggs(ham);
392  green_eggs.triangularView<Eigen::Upper>()/=lv;
393 
394  // FIXME we only transposehere to be consistent with the original implementation
395  // --- resulting in all that mangling of eigenvectors later
396  green_eggs.transposeInPlace();
397  green_eggs(0,0) += kr2GNe;
398 
399  return green_eggs;
400  }
401 
402 
403 
404 
405 
406  Eigen::Matrix3cd
408  const int & anti,
409  const OscParameters & params) const
410  {
411  double lv = 2 * kGeV2eV*E / params.Dmsq31(); // Osc. length in eV^-1
412  double kr2GNe = kK2*M_SQRT2*kGf * params.rho/2; // Matter potential in eV
413 
414  auto HLV = BuildHam(params);
415  HLV.triangularView<Eigen::Upper>()/=lv;
416 
417  if (anti>0) {
418  HLV.transposeInPlace(); // FIXME we only transposehere to be consistent with the original implementation --- resulting in all that mangling of eigenvectors later
419  HLV(0,0) += kr2GNe;
420  }
421  else {
422  HLV.adjointInPlace();
423  HLV(0,0) -= kr2GNe;
424  }
425  return HLV;
426  }
427 
428 
429  Eigen::Matrix3cd
431  {
432  // Create temp variables
433  double sij, cij, h00, h11, h01;
434  std::complex<double> expCP, h02, h12;
435 
436  // Hamiltonian in mass base. Only one entry is variable.
437  h11 = params.dmsq21 / params.Dmsq31();
438 
439  // Rotate over theta12
440  _sincos(params.th12, sij, cij);
441 
442  // There are 3 non-zero entries after rephasing so that h22 = 0
443  h00 = h11 * sij * sij - 1;
444  h01 = h11 * sij * cij;
445  h11 = h11 * cij * cij - 1;
446 
447  // Rotate over theta13 with deltaCP
448  _sincos(params.deltacp, sij, cij);
449  expCP = std::complex<double>(cij, -sij);
450 
451  _sincos(params.th13, sij, cij);
452 
453 
454  // There are 5 non-zero entries after rephasing so that h22 = 0
455  h02 = (-h00 * sij * cij) * expCP;
456  h12 = (-h01 * sij) * expCP;
457  h11 -= h00 * sij * sij;
458  h00 *= cij * cij - sij * sij;
459  h01 *= cij;
460 
461  // Finally, rotate over theta23
462  _sincos(params.th23, sij, cij);
463 
464 
465  // Fill the Hamiltonian rephased so that h22 = -h11
466  Eigen::Matrix3cd mHlv;
467  mHlv(0,0) = h00 - 0.5 * h11;
468  mHlv(1,1) = 0.5 * h11 * (cij * cij - sij * sij) + 2 * real(h12) * cij * sij;
469  mHlv(2,2) = -mHlv(1,1);
470 
471  mHlv(0,1) = h02 * sij + h01 * cij;
472  mHlv(0,2) = h02 * cij - h01 * sij;
473  mHlv(1,2) = h12 - (h11 * cij + 2 * real(h12) * sij) * sij;
474 
475  return mHlv;
476  }
477 }
unsigned int iter
Definition: OscParameters.h:67
TMD5 * GetParamsHash() const override
Use to check two calculators are in the same state.
Eigen::Matrix3cd BuildHam(const OscParameters &params) const
Eigen::Matrix3cd AddMatterEffects(const Eigen::Matrix3cd &ham, const double &E, const int &anti, const OscParameters &params) const
static const double kK2
mole/GeV^2/cm^3 to eV
Definition: PMNSOpt.h:42
Eigen::Matrix3cd AddMatterEffectsAnti(const Eigen::Matrix3cd &ham, const double &E, const OscParameters &params) const
IOscCalcAdjustable * Copy() const override
void _sincos(double theta, double &s, double &c)
Definition: OscCalcDMP.cxx:37
void PrintArrayAddresses(const double *data, int size) const
Eigen::Matrix3cd PropMatter(Eigen::Matrix3cd const &evec, Eigen::Vector3d const &evals, const OscParameters &params) const
void abs(TH1 *hist)
double P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
const XML_Char const XML_Char * data
Definition: expat.h:268
void SetCachedEnergies(std::vector< double > const &energies)
EigenSystem Solve(Eigen::Matrix3cd const &ham) const
const XML_Char * s
Definition: expat.h:262
Eigen::Matrix3cd MatterHamiltonian(const double &E, const int &anti, const OscParameters &params) const
static const double kGeV2eV
GeV to eV.
Definition: PMNSOpt.h:43
Float_t E
Definition: plot.C:20
Helper struct for the cache. Might not need this.
int ChannelCacheIdx(int flavBefore, int flavAfter) const
const Eigen::Vector3d values
Eigen::Array33d Probabilities(const Eigen::Matrix3cd &to, const Eigen::Matrix3cd &from)
_OscCache< double > fCache
const double j
Definition: BetheBloch.cxx:29
A re-optimized version of OscCalcPMNSOpt.
Oscillation probability calculators.
Definition: Calcs.h:5
OStream cout
Definition: OStream.cxx:6
std::vector< double > energies
Definition: OscParameters.h:63
TH1F * h11
Definition: plot.C:43
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void Zero()
Array< T, Dynamic, Dynamic > probabilities
Definition: OscParameters.h:65
const Eigen::Matrix3cd vectors
TH1D * h12
Definition: plot_hist.C:27
void PrintMatrixAddresses(const double *data, int nrows, int ncols) const
std::enable_if_t<!std::is_arithmetic_v< T >, void > sincos(const T &x, T *sx, T *cx)
Int_t ncols
Definition: plot.C:53
Eigen::MatrixXd mat
_OscParameters< T > parameters
Definition: OscParameters.h:66
static const double kKm2eV
km to eV^-1
Definition: PMNSOpt.h:41
static const double kGf
Definition: PMNSOpt.h:46