OscCalcGeneral.cxx
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////////
2 // \file OscCalcGeneral.cxx
3 // \brief Source file for general oscillation calculation
4 // \version $Id: OscCalcGeneral.cxx,v 1.5 2012-09-25 04:51:35 gsdavies Exp $
5 // \author
6 /////////////////////////////////////////////////////////////////////////////
8 
9 #include <boost/serialization/array_wrapper.hpp>
10 #pragma GCC diagnostic push
11 #pragma GCC diagnostic ignored "-Wmisleading-indentation"
12 #include <boost/numeric/ublas/vector.hpp>
13 #include <boost/numeric/ublas/matrix.hpp>
14 #pragma GCC diagnostic pop
15 
16 #include <vector>
17 
18 #include <cassert>
19 #include <complex>
20 
21 namespace osc
22 {
23  const unsigned int kNumFlavours = 3;
24 
25  namespace ublas = boost::numeric::ublas;
26  typedef std::complex<long double> val_t;
27  typedef ublas::bounded_array<val_t, kNumFlavours> alloc_t;
28  // Should be a c_matrix, but get ambiguity in triple multiplication
29  typedef ublas::bounded_matrix<val_t, kNumFlavours, kNumFlavours> ComplexMat;
30  typedef ublas::c_vector<val_t, kNumFlavours> ComplexVec;
31  typedef ublas::unit_vector<val_t, alloc_t> UnitVec;
32  const ublas::zero_matrix<val_t, alloc_t> kZeroMat(kNumFlavours, kNumFlavours);
33  const ublas::identity_matrix<val_t, alloc_t> kIdentity(kNumFlavours);
34 
35  // Private data. Doing it this way avoids having to expose complex numbers,
36  // matrices etc. to the user.
38  {
39  Priv() :
40  emutau(0),
44  pmns(kZeroMat),
45  dirty(true)
46  {
47  }
48 
49  long double c13, s13;
50  long double emutau;
51  std::complex<long double> phase;
52 
53  ComplexMat atmos, react, solar;
54  ComplexMat pmns;
55 
56  bool dirty;
57  };
58 
59 
62  {
63  }
64 
66  {
67  delete d;
68  }
69 
71  {
72  OscCalcGeneral* ret = new OscCalcGeneral(*this);
73  ret->d = new Priv;
74  *ret->d = *d;
75  return ret;
76  }
77 
78  void OscCalcGeneral::SetTh23(const double& th23)
79  {
80  fTh23 = th23;
81  d->atmos(2, 2) = d->atmos(1, 1) = cos(th23);
82  d->atmos(1, 2) = sin(th23);
83  d->atmos(2, 1) = -d->atmos(1, 2);
84  d->dirty = true;
85  }
86 
87  void OscCalcGeneral::SetTh13(const double& th13)
88  {
89  fTh13 = th13;
90  d->c13 = cos(th13);
91  d->s13 = sin(th13);
92 
93  d->react(2, 2) = d->react(0, 0) = d->c13;
94  d->react(0, 2) = d->s13*d->phase;
95  d->react(2, 0) = -std::conj(d->react(0, 2));
96  d->dirty = true;
97  }
98 
99  void OscCalcGeneral::SetTh12(const double& th12)
100  {
101  fTh12 = th12;
102  d->solar(1, 1) = d->solar(0, 0) = cos(th12);
103  d->solar(0, 1) = sin(th12);
104  d->solar(1, 0) = -d->solar(0, 1);
105  d->dirty = true;
106  }
107 
108  void OscCalcGeneral::SetdCP(const double& delta)
109  {
110  fdCP = delta;
111 
112  d->phase = std::polar((long double)1, -(long double)delta);
113 
114  d->react(2, 2) = d->react(0, 0) = d->c13;
115  d->react(0, 2) = d->s13*d->phase;
116  d->react(2, 0) = -std::conj(d->react(0, 2));
117  d->dirty = true;
118  }
119 
121  {
122  d->emutau = emutau;
123  }
124 
126  {
127  return d->emutau;
128  }
129 
131  {
132  // Default isn't good enough if we need to consider NSI
133  if(d->emutau) return 0;
135  }
136 
137 
139  {
140  if(d->dirty){
141  ublas::noalias(d->pmns) = ublas::prod(d->atmos, ublas::prod<ComplexMat>(d->react, d->solar));
142  d->dirty = false;
143  }
144  return d->pmns;
145  }
146 
147 
148  // U is the PMNS matrix
149  // mSq are the squared masses. Pick an arbitrary zero and use the dmsqs
150  // E is the neutrino energy
151  ComplexMat VacuumHamiltonian(const ComplexMat& U,
152  const std::vector<long double>& mSq,
153  long double E)
154  {
155  // Conversion factor for units
156  E /= 4*1.267;
157 
158  assert(mSq.size() == kNumFlavours);
159 
160  ComplexMat H = kZeroMat;
161 
162  // Loop over rows and columns of the output
163  for(unsigned int a = 0; a < kNumFlavours; ++a){
164  for(unsigned int b = 0; b < kNumFlavours; ++b){
165  // Form sum over mass states
166  for(unsigned int i = 0; i < kNumFlavours; ++i){
167  H(a, b) += U(a, i)*mSq[i]/(2*E)*std::conj(U(b, i));
168  }
169  }
170  }
171 
172  return H;
173  }
174 
175 
176  // This is calculated assuming matter neutrinos. So far, this can simply be
177  // converted to antineutrinos by taking a minus sign.
178  ComplexMat MatterHamiltonianComponent(long double Ne, long double emutau)
179  {
180  // Need to convert avogadro's constant so that the total term comes out in
181  // units of inverse distance. Note that Ne will be specified in g/cm^-3
182  // I put the following into Wolfram Alpha:
183  // (fermi coupling constant)*((avogadro number)/cm^3)*(reduced planck constant)^2*(speed of light)^2
184  // And then multiplied by 1000 because we specify L in km not m.
185  const long double GF = 1.368e-4;
186 
187  ComplexMat H = kZeroMat;
188 
189  H(0, 0) = sqrt(2)*GF*Ne;
190 
191  // Ignoring conjugates here because we assume e_mutau is real
192  H(1, 2) = H(2, 1) = emutau*sqrt(2)*GF*Ne;
193 
194  return H;
195  }
196 
197 
198  /* TODO: maybe look at
199  http://en.wikipedia.org/wiki/Baker%E2%80%93Campbell%E2%80%93Hausdorff_formula
200  and
201  http://en.wikipedia.org/wiki/Pad%C3%A9_approximant
202  http://en.wikipedia.org/wiki/Pad%C3%A9_table
203  */
204  ComplexMat MatrixExp(const ComplexMat& m2)
205  {
206  // We use the identity e^(a*b) = (e^a)^b
207  // First: ensure the exponent is really small, 65536 = 2^16
208  ComplexMat m = m2/65536;
209 
210  // To first order e^x = 1+x
211  ComplexMat ret = kIdentity+m;
212 
213  // Raise the result to the power 65536, by squaring it 16 times
214  for(int n = 0; n < 16; ++n) ret = ublas::prod(ret, ret);
215 
216  return ret;
217  }
218 
219 
220  ComplexVec EvolveState(ComplexVec A, const ComplexMat& H, long double L)
221  {
222  const std::complex<long double> i = std::complex<long double>(0, 1);
223 
224  return ublas::prod(MatrixExp(-H*L*i), A);
225  }
226 
227  void conjugate_elements(ComplexMat& m)
228  {
229  for(unsigned int i = 0; i < kNumFlavours; ++i)
230  for(unsigned int j = 0; j < kNumFlavours; ++j)
231  m(i, j) = std::conj(m(i, j));
232  }
233 
234  double OscCalcGeneral::P(int from, int to, double E)
235  {
236  const int af = abs(from);
237  const int at = abs(to);
238  assert(af == 12 || af == 14 || af == 16);
239  assert(at == 12 || at == 14 || at == 16);
240 
241  // No matter<->antimatter transitions
242  if(af*at < 0) return 0;
243 
244  ComplexMat U = GetPMNS(d);
245  // P(a->b|U) = P(abar->bbar|U*)
246  if(from < 0) conjugate_elements(U);
247 
248  ComplexVec initState = UnitVec(kNumFlavours, 1); // Display accelerator bias
249  if(af == 12) initState = UnitVec(kNumFlavours, 0);
250  if(af == 16) initState = UnitVec(kNumFlavours, 2);
251 
252  std::vector<long double> mSq;
253  mSq.push_back(0);
254  mSq.push_back(fDmsq21);
255  mSq.push_back(fDmsq21 + fDmsq32);
256 
257  ComplexMat H = VacuumHamiltonian(U, mSq, E);
258  ComplexMat Hm = MatterHamiltonianComponent(fRho, d->emutau);
259  // So far, contribution to the antineutrino Hamiltonian is just the negative
260  // If there were to be any complex stuff here, would have to think about
261  // how it transformed with antineutrinos.
262  if(from < 0) H += -1*Hm; else H += Hm;
263 
264  ComplexVec finalState = EvolveState(initState, H, fL);
265 
266  if(at == 12) return std::norm(finalState(0));
267  if(at == 14) return std::norm(finalState(1));
268  if(at == 16) return std::norm(finalState(2));
269 
270  assert(0 && "Not reached");
271 
272  // should never get to this point, but it quiets a compiler error
273  // a probability of 0 will be so obviously wrong that it should
274  // alert people - ie every probability will be 0
275  return 0.;
276  }
277 
278 } // namespace
ComplexMat MatrixExp(const ComplexMat &m2)
More generic (but probably slower) oscillation calculations.
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
virtual double P(int from, int to, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
T sqrt(T number)
Definition: d0nt_math.hpp:156
ComplexMat MatterHamiltonianComponent(long double Ne, long double emutau)
ComplexMat VacuumHamiltonian(const ComplexMat &U, const std::vector< long double > &mSq, long double E)
void conjugate_elements(ComplexMat &m)
void abs(TH1 *hist)
ublas::unit_vector< val_t, alloc_t > UnitVec
double th12
Definition: runWimpSim.h:98
std::complex< long double > phase
TMD5 * GetParamsHashDefault(const std::string &txt) const
This is only a safe implementation if your calculator only depends on these eight parameters...
virtual void SetTh13(const double &th13) override
virtual void SetTh12(const double &th12) override
static constexpr double L
Float_t E
Definition: plot.C:20
ublas::bounded_matrix< val_t, kNumFlavours, kNumFlavours > ComplexMat
const double a
virtual void SetTh23(const double &th23) override
void SetNSIEpsilonMuTau(double emutau)
std::complex< long double > val_t
const double j
Definition: BetheBloch.cxx:29
T prod(const std::vector< T > &v)
Definition: prod.hpp:17
virtual TMD5 * GetParamsHash() const override
Use to check two calculators are in the same state.
Oscillation probability calculators.
Definition: Calcs.h:5
double GetNSIEpsilonMuTau() const
const ublas::zero_matrix< val_t, alloc_t > kZeroMat(kNumFlavours, kNumFlavours)
const unsigned int kNumFlavours
static constexpr Double_t m2
Definition: Munits.h:145
static const double A
Definition: Units.h:82
T sin(T number)
Definition: d0nt_math.hpp:132
Float_t norm
virtual void SetdCP(const double &dCP) override
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
ublas::c_vector< val_t, kNumFlavours > ComplexVec
const ublas::identity_matrix< val_t, alloc_t > kIdentity(kNumFlavours)
ComplexVec EvolveState(ComplexVec A, const ComplexMat &H, long double L)
double th13
Definition: runWimpSim.h:98
ublas::bounded_array< val_t, kNumFlavours > alloc_t
virtual IOscCalcAdjustable * Copy() const override
ComplexMat GetPMNS(OscCalcGeneral::Priv *d)