Classes | Public Member Functions | Protected Types | Protected Member Functions | Protected Attributes | List of all members
osc::PMNS_Sterile Class Reference

#include "/cvmfs/nova.opensciencegrid.org/externals/osclib/v00.07/src/OscLib/PMNS_Sterile.h"

Classes

struct  Priv
 

Public Member Functions

 PMNS_Sterile (int NumNus)
 
virtual ~PMNS_Sterile ()
 
virtual void SetAngle (int i, int j, double th)
 
virtual void SetDelta (int i, int j, double delta)
 
virtual void SetDm (int i, double dm)
 
virtual void SetStdPars ()
 Set standard 3-flavor parameters. More...
 
virtual void PropMatter (double L, double E, double Ne, int anti=1)
 
virtual void PropMatter (const std::list< double > &L, double E, const std::list< double > &Ne, int anti)
 
virtual double P (int flv) const
 
virtual void ResetToFlavour (int flv=1)
 
virtual int GetNFlavors () const
 Getters. More...
 
virtual double GetDm (int i) const
 
virtual double GetAngle (int i, int j) const
 
virtual double GetDelta (int i, int j) const
 

Protected Types

typedef std::complex< double > complex
 

Protected Member Functions

virtual void InitializeVectors ()
 
virtual void RotateH (int i, int j)
 Rotate the Hamiltonian by theta_ij and delta_ij. More...
 
virtual void BuildHms ()
 
virtual void SolveHam (double E, double Ne, int anti)
 

Protected Attributes

int fNumNus
 
std::vector< double > fDm
 m^2_i - m^2_1 in vacuum More...
 
std::vector< std::vector< double > > fTheta
 theta[i][j] mixing angle More...
 
std::vector< std::vector< double > > fDelta
 delta[i][j] CP violating phase More...
 
std::vector< complexfNuState
 The neutrino current state. More...
 
std::vector< std::vector< complex > > fHms
 matrix H*2E in eV^2 More...
 
double fCachedNe
 Cached electron density. More...
 
double fCachedE
 Cached neutrino energy. More...
 
int fCachedAnti
 Cached nu/nubar selector. More...
 
bool fBuiltHms
 Tag to avoid rebuilding Hms. More...
 
Privd
 

Detailed Description

Definition at line 18 of file PMNS_Sterile.h.

Member Typedef Documentation

typedef std::complex<double> osc::PMNS_Sterile::complex
protected

Definition at line 65 of file PMNS_Sterile.h.

Constructor & Destructor Documentation

PMNS_Sterile::PMNS_Sterile ( int  NumNus)

Definition at line 96 of file PMNS_Sterile.cxx.

References fCachedAnti, fCachedE, fCachedNe, fNumNus, ResetToFlavour(), and SetStdPars().

97  : d(new PMNS_Sterile::Priv(NumNus))
98 {
99 
100  fNumNus = NumNus;
101  this->SetStdPars();
102  this->ResetToFlavour(1);
103  fCachedNe = 0.0;
104  fCachedE = 1.0;
105  fCachedAnti = 1;
106 }
double fCachedNe
Cached electron density.
Definition: PMNS_Sterile.h:91
virtual void SetStdPars()
Set standard 3-flavor parameters.
virtual void ResetToFlavour(int flv=1)
int fCachedAnti
Cached nu/nubar selector.
Definition: PMNS_Sterile.h:93
double fCachedE
Cached neutrino energy.
Definition: PMNS_Sterile.h:92
PMNS_Sterile::~PMNS_Sterile ( )
virtual

Definition at line 108 of file PMNS_Sterile.cxx.

References d.

109 {
110  delete d;
111 }

Member Function Documentation

void PMNS_Sterile::BuildHms ( )
protectedvirtual

Build Hms = H*2E, where H is the Hamiltonian in vacuum on flavour basis and E is the neutrino energy. This is effectively the matrix of masses squared.

Build Hms = H*2E, where H is the Hamiltonian in vacuum on flavour basis and E is the neutrino energy in eV. Hms is effectively the matrix of masses squared.

This is a hermitian matrix, so only the upper triangular part needs to be filled

The construction of the Hamiltonian avoids computing terms that are simply zero. This has a big impact in the computation time. This construction is described in DocDB-XXXX (to be posted)

Definition at line 361 of file PMNS_Sterile.cxx.

References fBuiltHms, fDm, fHms, fNumNus, MECModelEnuComparisons::i, calib::j, and RotateH().

Referenced by SolveHam().

362 {
363 
364  // Check if anything changed
365  if(fBuiltHms) return;
366 
367  for(int j=0; j<fNumNus; j++){
368  // Set mass splitting
369  fHms[j][j] = fDm[j];
370  // Reset off-diagonal elements
371  for(int i=0; i<j; i++){
372  fHms[i][j] = 0;
373  }
374  // Rotate j neutrinos
375  for(int i=0; i<j; i++){
376  this->RotateH(i,j);
377  }
378  }
379 
380  // Tag as built
381  fBuiltHms = true;
382 
383 }
std::vector< double > fDm
m^2_i - m^2_1 in vacuum
Definition: PMNS_Sterile.h:84
virtual void RotateH(int i, int j)
Rotate the Hamiltonian by theta_ij and delta_ij.
const double j
Definition: BetheBloch.cxx:29
bool fBuiltHms
Tag to avoid rebuilding Hms.
Definition: PMNS_Sterile.h:94
std::vector< std::vector< complex > > fHms
matrix H*2E in eV^2
Definition: PMNS_Sterile.h:89
virtual double osc::PMNS_Sterile::GetAngle ( int  i,
int  j 
) const
inlinevirtual

Definition at line 59 of file PMNS_Sterile.h.

References fTheta.

Referenced by osc::OscCalcSterile::GetAngle().

59 { return fTheta[i-1][j-1]; }
const double j
Definition: BetheBloch.cxx:29
std::vector< std::vector< double > > fTheta
theta[i][j] mixing angle
Definition: PMNS_Sterile.h:85
virtual double osc::PMNS_Sterile::GetDelta ( int  i,
int  j 
) const
inlinevirtual

Definition at line 60 of file PMNS_Sterile.h.

References fDelta.

Referenced by osc::OscCalcSterile::GetDelta().

60 { return fDelta[i-1][j-1]; }
const double j
Definition: BetheBloch.cxx:29
std::vector< std::vector< double > > fDelta
delta[i][j] CP violating phase
Definition: PMNS_Sterile.h:86
virtual double osc::PMNS_Sterile::GetDm ( int  i) const
inlinevirtual

Definition at line 58 of file PMNS_Sterile.h.

References fDm.

Referenced by osc::OscCalcSterile::GetDm().

58 { return fDm[i-1]; }
std::vector< double > fDm
m^2_i - m^2_1 in vacuum
Definition: PMNS_Sterile.h:84
virtual int osc::PMNS_Sterile::GetNFlavors ( ) const
inlinevirtual

Getters.

Definition at line 57 of file PMNS_Sterile.h.

References fNumNus.

Referenced by osc::OscCalcSterile::GetNFlavors().

57 { return fNumNus; }
void PMNS_Sterile::InitializeVectors ( )
protectedvirtual

Definition at line 113 of file PMNS_Sterile.cxx.

References fDelta, fDm, fHms, fNumNus, fNuState, fTheta, and zero().

Referenced by SetStdPars().

114 {
115 
116  fDm = vector<double>(fNumNus, 0);
117  fTheta = vector< vector<double> >(fNumNus, vector<double>(fNumNus,0));
118  fDelta = vector< vector<double> >(fNumNus, vector<double>(fNumNus,0));
119 
120  fNuState = vector<complex>(fNumNus, zero);
121  fHms = vector< vector<complex> >(fNumNus, vector<complex>(fNumNus,zero));
122 
123 }
std::vector< double > fDm
m^2_i - m^2_1 in vacuum
Definition: PMNS_Sterile.h:84
std::vector< complex > fNuState
The neutrino current state.
Definition: PMNS_Sterile.h:88
std::vector< std::vector< double > > fTheta
theta[i][j] mixing angle
Definition: PMNS_Sterile.h:85
std::vector< std::vector< complex > > fHms
matrix H*2E in eV^2
Definition: PMNS_Sterile.h:89
std::vector< std::vector< double > > fDelta
delta[i][j] CP violating phase
Definition: PMNS_Sterile.h:86
static std::complex< double > zero(0.0, 0.0)
double PMNS_Sterile::P ( int  flv) const
virtual

Return the probability of final state in flavour flv

Parameters
flv- final flavor (0,1,2) = (nue,numu,nutau)

Compute oscillation probability of flavour flv

0 = nue, 1 = numu, 2 = nutau

Definition at line 509 of file PMNS_Sterile.cxx.

References ana::assert(), fNumNus, fNuState, and norm.

Referenced by osc::OscCalcSterile::P().

510 {
511  assert(flv>=0 && flv<fNumNus);
512  return norm(fNuState[flv]);
513 }
std::vector< complex > fNuState
The neutrino current state.
Definition: PMNS_Sterile.h:88
Float_t norm
assert(nhit_max >=nhit_nbins)
void PMNS_Sterile::PropMatter ( double  L,
double  E,
double  Ne,
int  anti = 1 
)
virtual

Propagate a neutrino through a slab of matter

Parameters
L- length of slab (km)
E- neutrino energy in GeV
Ne- electron number density of matter in mole/cm^3
anti- +1 = neutrino case, -1 = anti-neutrino case

.....................................................................

Propagate the current neutrino state over a distance L in km with an energy E in GeV through constant matter of density Ne in mole/cm^3.

Parameters
anti- +1 = neutrino case, -1 = anti-neutrino case

Definition at line 439 of file PMNS_Sterile.cxx.

References d, stan::math::exp(), osc::PMNS_Sterile::Priv::fEval, osc::PMNS_Sterile::Priv::fEvec, fNumNus, fNuState, MECModelEnuComparisons::i, calib::j, osc::kKm2eV, SolveHam(), and zero().

Referenced by osc::OscCalcSterile::P(), and PropMatter().

440 {
441 
442  // Solve Hamiltonian
443  this->SolveHam(E, Ne, anti);
444 
445  // Store coefficients of propagation eigenstates
446  vector<complex> nuComp(fNumNus, zero);
447  for(int i=0;i<fNumNus;i++){
448  nuComp[i] = 0;
449  for(int j=0;j<fNumNus;j++){
450  gsl_complex buf = gsl_matrix_complex_get(d->fEvec,j,i);
451  complex evecji = complex( GSL_REAL(buf), GSL_IMAG(buf) );
452  nuComp[i] += fNuState[j] * conj( evecji );
453  }
454  }
455 
456  // Propagate neutrino state
457  for(int i=0;i<fNumNus;i++){
458  fNuState[i] = 0;
459  for(int j=0;j<fNumNus;j++){
460  gsl_complex buf = gsl_matrix_complex_get(d->fEvec,i,j);
461  complex evecij = complex( GSL_REAL(buf), GSL_IMAG(buf) );
462  complex iEval(0.0,gsl_vector_get(d->fEval,j));
463  fNuState[i] += exp(-iEval * kKm2eV*L) * nuComp[j] * evecij;
464  }
465  }
466 
467 }
gsl_matrix_complex * fEvec
static constexpr double L
Float_t E
Definition: plot.C:20
std::vector< complex > fNuState
The neutrino current state.
Definition: PMNS_Sterile.h:88
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
std::complex< double > complex
Definition: PMNS_Sterile.h:65
const double j
Definition: BetheBloch.cxx:29
virtual void SolveHam(double E, double Ne, int anti)
static const double kKm2eV
km to eV^-1
Definition: PMNSOpt.h:41
static std::complex< double > zero(0.0, 0.0)
void PMNS_Sterile::PropMatter ( const std::list< double > &  L,
double  E,
const std::list< double > &  Ne,
int  anti 
)
virtual

Do several layers in a row. L and Ne must have the same length

Definition at line 473 of file PMNS_Sterile.cxx.

References PropMatter().

477 {
478  if (L.size()!=Ne.size()) abort();
479  std::list<double>::const_iterator Li (L.begin());
480  std::list<double>::const_iterator Lend(L.end());
481  std::list<double>::const_iterator Ni (Ne.begin());
482  for (; Li!=Lend; ++Li, ++Ni) {
483  //cout << *Li << " km and " << *Ni << " mol/cm^3" << endl;
484  this->PropMatter(*Li, E, *Ni, anti);
485  }
486 }
virtual void PropMatter(double L, double E, double Ne, int anti=1)
static constexpr double L
Float_t E
Definition: plot.C:20
void PMNS_Sterile::ResetToFlavour ( int  flv = 1)
virtual

Erase memory of neutrino propagate and reset neutrino to pure flavour flv. Preserves values of mixing and mass-splittings

Parameters
flv- final flavor (0,1,2) = (nue,numu,nutau)

Reset the neutrino state back to a pure flavour where it starts

Definition at line 494 of file PMNS_Sterile.cxx.

References fNumNus, fNuState, MECModelEnuComparisons::i, one(), and zero().

Referenced by osc::OscCalcSterile::P(), and PMNS_Sterile().

495 {
496  int i;
497  for (i=0; i<fNumNus; ++i){
498  if (i==flv) fNuState[i] = one;
499  else fNuState[i] = zero;
500  }
501 }
std::vector< complex > fNuState
The neutrino current state.
Definition: PMNS_Sterile.h:88
static std::complex< double > one(1.0, 0.0)
static std::complex< double > zero(0.0, 0.0)
void PMNS_Sterile::RotateH ( int  i,
int  j 
)
protectedvirtual

Rotate the Hamiltonian by theta_ij and delta_ij.

Rotate the Hamiltonian by the angle theta_ij and phase delta_ij. The rotations assume all off-diagonal elements with i > j are zero. This is correct if the order of rotations is chosen appropriately and it speeds up computation by skipping null terms

Definition at line 226 of file PMNS_Sterile.cxx.

References std::cos(), fDelta, fHms, fTheta, MECModelEnuComparisons::i, calib::j, and std::sin().

Referenced by BuildHms().

226  {
227 
228  // Do nothing if angle is zero
229  if(fTheta[i][j]==0) return;
230 
231  double fSinBuffer = sin(fTheta[i][j]);
232  double fCosBuffer = cos(fTheta[i][j]);
233 
234  double fHmsBufferD;
235  complex fHmsBufferC;
236 
237  // With Delta
238  if(i+1<j){
239  complex fExpBuffer = complex(cos(fDelta[i][j]), -sin(fDelta[i][j]));
240 
241  // General case
242  if(i>0){
243  // Top columns
244  for(int k=0; k<i; k++){
245  fHmsBufferC = fHms[k][i];
246 
247  fHms[k][i] *= fCosBuffer;
248  fHms[k][i] += fHms[k][j] * fSinBuffer * conj(fExpBuffer);
249 
250  fHms[k][j] *= fCosBuffer;
251  fHms[k][j] -= fHmsBufferC * fSinBuffer * fExpBuffer;
252  }
253 
254  // Middle row and column
255  for(int k=i+1; k<j; k++){
256  fHmsBufferC = fHms[k][j];
257 
258  fHms[k][j] *= fCosBuffer;
259  fHms[k][j] -= conj(fHms[i][k]) * fSinBuffer * fExpBuffer;
260 
261  fHms[i][k] *= fCosBuffer;
262  fHms[i][k] += fSinBuffer * fExpBuffer * conj(fHmsBufferC);
263  }
264 
265  // Nodes ij
266  fHmsBufferC = fHms[i][i];
267  fHmsBufferD = real(fHms[j][j]);
268 
269  fHms[i][i] *= fCosBuffer * fCosBuffer;
270  fHms[i][i] += 2 * fSinBuffer * fCosBuffer * real(fHms[i][j] * conj(fExpBuffer));
271  fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;
272 
273  fHms[j][j] *= fCosBuffer * fCosBuffer;
274  fHms[j][j] += fSinBuffer * fHmsBufferC * fSinBuffer;
275  fHms[j][j] -= 2 * fSinBuffer * fCosBuffer * real(fHms[i][j] * conj(fExpBuffer));
276 
277  fHms[i][j] -= 2 * fSinBuffer * real(fHms[i][j] * conj(fExpBuffer)) * fSinBuffer * fExpBuffer;
278  fHms[i][j] += fSinBuffer * fCosBuffer * (fHmsBufferD - fHmsBufferC) * fExpBuffer;
279 
280  }
281  // First rotation on j (No top columns)
282  else{
283  // Middle rows and columns
284  for(int k=i+1; k<j; k++){
285  fHms[k][j] = -conj(fHms[i][k]) * fSinBuffer * fExpBuffer;
286 
287  fHms[i][k] *= fCosBuffer;
288  }
289 
290  // Nodes ij
291  fHmsBufferD = real(fHms[i][i]);
292 
293  fHms[i][j] = fSinBuffer * fCosBuffer * (fHms[j][j] - fHmsBufferD) * fExpBuffer;
294 
295  fHms[i][i] *= fCosBuffer * fCosBuffer;
296  fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;
297 
298  fHms[j][j] *= fCosBuffer * fCosBuffer;
299  fHms[j][j] += fSinBuffer * fHmsBufferD * fSinBuffer;
300  }
301 
302  }
303  // Without Delta (No middle rows or columns: j = i+1)
304  else{
305  // General case
306  if(i>0){
307  // Top columns
308  for(int k=0; k<i; k++){
309  fHmsBufferC = fHms[k][i];
310 
311  fHms[k][i] *= fCosBuffer;
312  fHms[k][i] += fHms[k][j] * fSinBuffer;
313 
314  fHms[k][j] *= fCosBuffer;
315  fHms[k][j] -= fHmsBufferC * fSinBuffer;
316  }
317 
318  // Nodes ij
319  fHmsBufferC = fHms[i][i];
320  fHmsBufferD = real(fHms[j][j]);
321 
322  fHms[i][i] *= fCosBuffer * fCosBuffer;
323  fHms[i][i] += 2 * fSinBuffer * fCosBuffer * real(fHms[i][j]);
324  fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;
325 
326  fHms[j][j] *= fCosBuffer * fCosBuffer;
327  fHms[j][j] += fSinBuffer * fHmsBufferC * fSinBuffer;
328  fHms[j][j] -= 2 * fSinBuffer * fCosBuffer * real(fHms[i][j]);
329 
330  fHms[i][j] -= 2 * fSinBuffer * real(fHms[i][j]) * fSinBuffer;
331  fHms[i][j] += fSinBuffer * fCosBuffer * (fHmsBufferD - fHmsBufferC);
332 
333  }
334  // First rotation (theta12)
335  else{
336 
337  fHms[i][j] = fSinBuffer * fCosBuffer * fHms[j][j];
338 
339  fHms[i][i] = fSinBuffer * fHms[j][j] * fSinBuffer;
340 
341  fHms[j][j] *= fCosBuffer * fCosBuffer;
342 
343  }
344  }
345 
346 }
std::complex< double > complex
Definition: PMNS_Sterile.h:65
const double j
Definition: BetheBloch.cxx:29
T sin(T number)
Definition: d0nt_math.hpp:132
std::vector< std::vector< double > > fTheta
theta[i][j] mixing angle
Definition: PMNS_Sterile.h:85
T cos(T number)
Definition: d0nt_math.hpp:78
std::vector< std::vector< complex > > fHms
matrix H*2E in eV^2
Definition: PMNS_Sterile.h:89
std::vector< std::vector< double > > fDelta
delta[i][j] CP violating phase
Definition: PMNS_Sterile.h:86
void PMNS_Sterile::SetAngle ( int  i,
int  j,
double  th 
)
virtual

Set the parameters of the PMNS matrix

Parameters
th- The angle theta_ij in radians

Set mixing angles in radians. Requires i < j.

Definition at line 148 of file PMNS_Sterile.cxx.

References om::cout, allTimeWatchdog::endl, fBuiltHms, fNumNus, fTheta, MECModelEnuComparisons::i, calib::j, and APDHVSetting::temp.

Referenced by osc::OscCalcSterile::SetAngle(), and SetStdPars().

149 {
150 
151  if(i>j){
152  cout << "First argument should be smaller than second argument" << endl;
153  cout << "Setting reverse order (Theta" << j << i << "). " << endl;
154  int temp = i;
155  i = j;
156  j = temp;
157  }
158  if(i<1 || i>fNumNus-1 || j<2 || j>fNumNus){
159  cout << "Theta" << i << j << " not valid for " << fNumNus;
160  cout << " neutrinos. Doing nothing." << endl;
161  return;
162  }
163 
164  fTheta[i-1][j-1] = th;
165 
166  fBuiltHms = false;
167 
168 }
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
bool fBuiltHms
Tag to avoid rebuilding Hms.
Definition: PMNS_Sterile.h:94
std::vector< std::vector< double > > fTheta
theta[i][j] mixing angle
Definition: PMNS_Sterile.h:85
void PMNS_Sterile::SetDelta ( int  i,
int  j,
double  delta 
)
virtual

Set CP violating phases in radians. Requires i+1 < j.

Definition at line 174 of file PMNS_Sterile.cxx.

References om::cout, delta, allTimeWatchdog::endl, fBuiltHms, fDelta, fNumNus, MECModelEnuComparisons::i, calib::j, and APDHVSetting::temp.

Referenced by osc::OscCalcSterile::SetDelta().

175 {
176 
177  if(i>j){
178  cout << "First argument should be smaller than second argument" << endl;
179  cout << "Setting reverse order (Delta" << j << i << "). " << endl;
180  int temp = i;
181  i = j;
182  j = temp;
183  }
184  if(i<1 || i>fNumNus-1 || j<2 || j>fNumNus){
185  cout << "Delta" << i << j << " not valid for " << fNumNus;
186  cout << " neutrinos. Doing nothing." << endl;
187  return;
188  }
189  if(i+1==j){
190  cout << "Rotation " << i << j << " is real. Doing nothing." << endl;
191  return;
192  }
193 
194  fDelta[i-1][j-1] = delta;
195 
196  fBuiltHms = false;
197 
198 }
double delta
Definition: runWimpSim.h:98
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
bool fBuiltHms
Tag to avoid rebuilding Hms.
Definition: PMNS_Sterile.h:94
std::vector< std::vector< double > > fDelta
delta[i][j] CP violating phase
Definition: PMNS_Sterile.h:86
void PMNS_Sterile::SetDm ( int  j,
double  dm 
)
virtual

Set the mass-splittings

Parameters
dmi1- mi^2-m1^2 in eV^2

Set the mass-splittings. These are m_j^2-m_1^2

Definition at line 204 of file PMNS_Sterile.cxx.

References om::cout, allTimeWatchdog::endl, fBuiltHms, fDm, and fNumNus.

Referenced by osc::OscCalcSterile::SetDm(), and SetStdPars().

205 {
206 
207  if(j<2 || j>fNumNus){
208  cout << "Dm" << j << "1 not valid for " << fNumNus;
209  cout << " neutrinos. Doing nothing." << endl;
210  return;
211  }
212 
213  fDm[j-1] = dm;
214 
215  fBuiltHms = false;
216 
217 }
std::vector< double > fDm
m^2_i - m^2_1 in vacuum
Definition: PMNS_Sterile.h:84
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
bool fBuiltHms
Tag to avoid rebuilding Hms.
Definition: PMNS_Sterile.h:94
void PMNS_Sterile::SetStdPars ( )
virtual

Set standard 3-flavor parameters.

Definition at line 125 of file PMNS_Sterile.cxx.

References e, fNumNus, InitializeVectors(), SetAngle(), and SetDm().

Referenced by PMNS_Sterile().

126 {
127 
128  this->InitializeVectors();
129 
130  if(fNumNus>2){
131  this->SetAngle(1,2,0.6);
132  this->SetAngle(1,3,0.16);
133  this->SetAngle(2,3,0.7);
134  this->SetDm(2,7.6e-5);
135  this->SetDm(3,2.4e-3);
136  }
137  else if(fNumNus==2){
138  this->SetAngle(1,2,0.7);
139  this->SetDm(2,2.4e-3);
140  }
141 
142 }
virtual void SetAngle(int i, int j, double th)
virtual void SetDm(int i, double dm)
Float_t e
Definition: plot.C:35
virtual void InitializeVectors()
void PMNS_Sterile::SolveHam ( double  E,
double  Ne,
int  anti 
)
protectedvirtual

Solve the full Hamiltonian for eigenvectors and eigenvalues

Parameters
E- neutrino energy in GeV
Ne- electron number density of matter in mole/cm^3
anti- +1 = neutrino case, -1 = anti-neutrino case

Solve the full Hamiltonian for eigenvectors and eigenvalues This uses GSL to solve the eigensystem.

The matter effect assumes # of neutrons = # of electrons.

Definition at line 392 of file PMNS_Sterile.cxx.

References BuildHms(), d, E, fBuiltHms, fCachedAnti, fCachedE, fCachedNe, osc::PMNS_Sterile::Priv::fEval, osc::PMNS_Sterile::Priv::fEvec, fHms, fNumNus, osc::PMNS_Sterile::Priv::H_GSL, MECModelEnuComparisons::i, calib::j, osc::kGeV2eV, osc::kGf, osc::kK2, and osc::PMNS_Sterile::Priv::W_GSL.

Referenced by PropMatter().

393 {
394 
395  // Check if anything has changed before recalculating
396  if(Ne!=fCachedNe || E!=fCachedE || anti!=fCachedAnti || !fBuiltHms ){
397  fCachedNe = Ne;
398  fCachedE = E;
399  fCachedAnti = anti;
400  this->BuildHms();
401  }
402  else return;
403 
404  double lv = 2 * kGeV2eV*E; // 2E in eV
405  double kr2GNe = kK2*M_SQRT2*kGf*Ne; // Matter potential in eV
406 
407  // Finish building Hamiltonian in matter with dimension of eV
408 
409  for(size_t i=0;i<size_t(fNumNus);i++){
410  complex buf = fHms[i][i]/lv;
411  *gsl_matrix_complex_ptr(d->H_GSL, i, i) = gsl_complex_rect(real(buf),0);
412  for(size_t j=i+1;j<size_t(fNumNus);j++){
413  buf = fHms[i][j]/lv;
414  if(anti>0) *gsl_matrix_complex_ptr(d->H_GSL, j, i) = gsl_complex_rect(real(buf),-imag(buf));
415  else *gsl_matrix_complex_ptr(d->H_GSL, j, i) = gsl_complex_rect(real(buf), imag(buf));
416  }
417  if(i>2){
418  // Subtract NC coherent forward scattering from sterile neutrinos. See arXiv:hep-ph/0606054v3, eq. 3.15, for example.
419  if(anti>0) *gsl_matrix_complex_ptr(d->H_GSL, i, i) = gsl_complex_add_real(gsl_matrix_complex_get(d->H_GSL,i,i) , kr2GNe/2);
420  else *gsl_matrix_complex_ptr(d->H_GSL, i, i) = gsl_complex_add_real(gsl_matrix_complex_get(d->H_GSL,i,i) , -kr2GNe/2);;
421  }
422  }
423  // Add nue CC coherent forward scattering from sterile neutrinos.
424  if(anti>0) *gsl_matrix_complex_ptr(d->H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(d->H_GSL,0,0) , kr2GNe);
425  else *gsl_matrix_complex_ptr(d->H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(d->H_GSL,0,0) , -kr2GNe);
426 
427  // Solve Hamiltonian for eigensystem
428  gsl_eigen_hermv(d->H_GSL, d->fEval, d->fEvec, d->W_GSL);
429 
430 }
static const double kK2
mole/GeV^2/cm^3 to eV
Definition: PMNSOpt.h:42
gsl_matrix_complex * fEvec
double fCachedNe
Cached electron density.
Definition: PMNS_Sterile.h:91
virtual void BuildHms()
gsl_matrix_complex * H_GSL
static const double kGeV2eV
GeV to eV.
Definition: PMNSOpt.h:43
gsl_eigen_hermv_workspace * W_GSL
Float_t E
Definition: plot.C:20
const double j
Definition: BetheBloch.cxx:29
int fCachedAnti
Cached nu/nubar selector.
Definition: PMNS_Sterile.h:93
bool fBuiltHms
Tag to avoid rebuilding Hms.
Definition: PMNS_Sterile.h:94
double fCachedE
Cached neutrino energy.
Definition: PMNS_Sterile.h:92
std::vector< std::vector< complex > > fHms
matrix H*2E in eV^2
Definition: PMNS_Sterile.h:89
static const double kGf
Definition: PMNSOpt.h:46

Member Data Documentation

Priv* osc::PMNS_Sterile::d
protected

Definition at line 96 of file PMNS_Sterile.h.

Referenced by PropMatter(), SolveHam(), and ~PMNS_Sterile().

bool osc::PMNS_Sterile::fBuiltHms
protected

Tag to avoid rebuilding Hms.

Definition at line 94 of file PMNS_Sterile.h.

Referenced by BuildHms(), SetAngle(), SetDelta(), SetDm(), and SolveHam().

int osc::PMNS_Sterile::fCachedAnti
protected

Cached nu/nubar selector.

Definition at line 93 of file PMNS_Sterile.h.

Referenced by PMNS_Sterile(), and SolveHam().

double osc::PMNS_Sterile::fCachedE
protected

Cached neutrino energy.

Definition at line 92 of file PMNS_Sterile.h.

Referenced by PMNS_Sterile(), and SolveHam().

double osc::PMNS_Sterile::fCachedNe
protected

Cached electron density.

Definition at line 91 of file PMNS_Sterile.h.

Referenced by PMNS_Sterile(), and SolveHam().

std::vector< std::vector<double> > osc::PMNS_Sterile::fDelta
protected

delta[i][j] CP violating phase

Definition at line 86 of file PMNS_Sterile.h.

Referenced by GetDelta(), InitializeVectors(), RotateH(), and SetDelta().

std::vector<double> osc::PMNS_Sterile::fDm
protected

m^2_i - m^2_1 in vacuum

Definition at line 84 of file PMNS_Sterile.h.

Referenced by BuildHms(), GetDm(), InitializeVectors(), and SetDm().

std::vector< std::vector<complex> > osc::PMNS_Sterile::fHms
protected

matrix H*2E in eV^2

Definition at line 89 of file PMNS_Sterile.h.

Referenced by BuildHms(), InitializeVectors(), RotateH(), and SolveHam().

int osc::PMNS_Sterile::fNumNus
protected
std::vector<complex> osc::PMNS_Sterile::fNuState
protected

The neutrino current state.

Definition at line 88 of file PMNS_Sterile.h.

Referenced by InitializeVectors(), P(), PropMatter(), and ResetToFlavour().

std::vector< std::vector<double> > osc::PMNS_Sterile::fTheta
protected

theta[i][j] mixing angle

Definition at line 85 of file PMNS_Sterile.h.

Referenced by GetAngle(), InitializeVectors(), RotateH(), and SetAngle().


The documentation for this class was generated from the following files: