OscCalcPMNSOpt.cxx
Go to the documentation of this file.
2 
3 namespace osc
4 {
5  //---------------------------------------------------------------------------
6  template <typename T>
8  : fMixIdx(0), fDmIdx(0), fLRIdx(0)
9  {
10  }
11 
12  //---------------------------------------------------------------------------
13  template <typename T>
15  {
16  for(int i = 0; i < 2; ++i)
17  for(auto it: fPMNSOpt[i])
18  delete it.second.pmns;
19  }
20 
21  //---------------------------------------------------------------------------
22  template <typename T>
24  {
26  // Raw pointers were blindly copied, so we'd be in trouble when the
27  // destructors are called. More importantly, having two calculators sharing
28  // one PMNS object will lead to both of them getting confused as to whether
29  // they're up-to-date or not. Just clear out the cache entirely in the copy
30  // and let it be repopulated.
31  ret->fPMNSOpt[0].clear();
32  ret->fPMNSOpt[1].clear();
33  return ret;
34  }
35 
36  //---------------------------------------------------------------------------
37  template <typename T>
38  T _OscCalcPMNSOpt<T>::P(int flavBefore, int flavAfter, double E)
39  {
40  // Normal usage of a calculator in a fit is to configure one set of
41  // oscillation parameters and then calculate probabilities for all flavours
42  // and for a range of energies. The underlying PMNSOpt object has to redo
43  // its expensive diagonalization every time the energy changes, but not
44  // when simply switching to a different flavour. Unfortunately callers in
45  // CAFAna wind up making the calls in the wrong order to take advantage of
46  // this. So, we maintain a map of calculators for different energies (also
47  // for neutrinos and antineutrinos) and select the right one if available,
48  // on the hope that it's already configured to do what we're asked, since
49  // it might already have done a different flavour. This works because
50  // callers often pass the same set of energies (bin centers) over and over
51  // again.
52 
53  // If the caches get too large clear them. The caller is probably giving
54  // the precise energy of each event?
55  for(int i = 0; i < 2; ++i){
56  if(fPMNSOpt[i].size() > 10000){
57  for(auto it: fPMNSOpt[i]) delete it.second.pmns;
58  fPMNSOpt[i].clear();
59  }
60  }
61 
62  const int anti = (flavBefore > 0) ? +1 : -1;
63  assert(flavAfter/anti > 0);
64 
65  Val_t& calc = fPMNSOpt[(1+anti)/2][E];
66  if(!calc.pmns) calc.pmns = new _PMNSOpt<T>;
67 
68  int i = -1, j = -1;
69  if(abs(flavBefore) == 12) i = 0;
70  if(abs(flavBefore) == 14) i = 1;
71  if(abs(flavBefore) == 16) i = 2;
72  if(abs(flavAfter) == 12) j = 0;
73  if(abs(flavAfter) == 14) j = 1;
74  if(abs(flavAfter) == 16) j = 2;
75  assert(i >= 0 && j >= 0);
76 
77  const bool dirty = (fMixIdx > calc.mixIdx ||
78  fDmIdx > calc.dmIdx ||
79  fLRIdx > calc.lrIdx);
80 
81  // If the calculator has out of date parameters update them
82  if(fMixIdx > calc.mixIdx){
83  calc.pmns->SetMix(this->fTh12, this->fTh23, this->fTh13, this->fdCP);
84  calc.mixIdx = fMixIdx;
85  }
86  if(fDmIdx > calc.dmIdx){
87  calc.pmns->SetDeltaMsqrs(this->fDmsq21, this->fDmsq32);
88  calc.dmIdx = fDmIdx;
89  }
90 
91  if(dirty){
92  // Cache results for all nine flavour combinations
93  for(int ii = 0; ii < 3; ++ii){
94  calc.pmns->ResetToFlavour(ii);
95  // Assume Z/A=0.5
96  const double Ne = this->fRho/2;
97  calc.pmns->PropMatter(this->fL, E, Ne, anti);
98  for(int jj = 0; jj < 3; ++jj){
99  calc.P[ii][jj] = calc.pmns->P(jj);
100  }
101  }
102 
103  calc.lrIdx = fLRIdx;
104  }
105 
106  return calc.P[i][j];
107  }
108 }
109 
110 //---------------------------------------------------------------------------
111 template class osc::_OscCalcPMNSOpt<double>;
112 
113 #ifdef OSCLIB_STAN
114 #include "stan/math/rev/scal.hpp"
116 #endif
set< int >::iterator it
Optimized version of PMNS.
Definition: PMNSOpt.h:50
void abs(TH1 *hist)
osc::OscCalcDumb calc
Float_t E
Definition: plot.C:20
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
T P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
const double j
Definition: BetheBloch.cxx:29
_IOscCalcAdjustable< T > * Copy() const override
Oscillation probability calculators.
Definition: Calcs.h:5
T P[3][3]
Cache of oscillation probabilities.
General interface to any calculator that lets you set the parameters.
assert(nhit_max >=nhit_nbins)
std::unordered_map< double, Val_t > fPMNSOpt[2]
double T
Definition: Xdiff_gwt.C:5
_PMNSOpt< T > * pmns
The calculator itself.