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< complex >	fNuState
	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...
Priv *	d

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().

  : d(new PMNS_Sterile::Priv(NumNus))
{

  fNumNus = NumNus;
  this->SetStdPars();
  this->ResetToFlavour(1);
  fCachedNe = 0.0;
  fCachedE =  1.0;
  fCachedAnti = 1;
}

PMNS_Sterile::~PMNS_Sterile ( )

virtual

Definition at line 108 of file PMNS_Sterile.cxx.

References d.

{
  delete d;
}

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().

{

  // Check if anything changed
  if(fBuiltHms) return;

  for(int j=0; j<fNumNus; j++){
    // Set mass splitting
    fHms[j][j] = fDm[j];
    // Reset off-diagonal elements
    for(int i=0; i<j; i++){
      fHms[i][j] = 0;
    }
    // Rotate j neutrinos
    for(int i=0; i<j; i++){
      this->RotateH(i,j);
    }
  }
  
  // Tag as built
  fBuiltHms = true;

}

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().

{ return fTheta[i-1][j-1]; }

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().

{ return fDelta[i-1][j-1]; }

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().

{ return fDm[i-1]; }

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().

{ 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().

{

  fDm    = vector<double>(fNumNus, 0);
  fTheta = vector< vector<double> >(fNumNus, vector<double>(fNumNus,0));
  fDelta = vector< vector<double> >(fNumNus, vector<double>(fNumNus,0));

  fNuState = vector<complex>(fNumNus, zero);
  fHms     = vector< vector<complex> >(fNumNus, vector<complex>(fNumNus,zero));

}

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().

{
  assert(flv>=0 && flv<fNumNus);
  return norm(fNuState[flv]);
}

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().

{

  // Solve Hamiltonian
  this->SolveHam(E, Ne, anti);

  // Store coefficients of propagation eigenstates
  vector<complex> nuComp(fNumNus, zero);
  for(int i=0;i<fNumNus;i++){
    nuComp[i] = 0;
    for(int j=0;j<fNumNus;j++){
      gsl_complex buf = gsl_matrix_complex_get(d->fEvec,j,i);
      complex evecji = complex( GSL_REAL(buf), GSL_IMAG(buf) );
      nuComp[i] += fNuState[j] * conj( evecji );
    }
  }

  // Propagate neutrino state
  for(int i=0;i<fNumNus;i++){
    fNuState[i] = 0;
    for(int j=0;j<fNumNus;j++){
      gsl_complex buf = gsl_matrix_complex_get(d->fEvec,i,j);
      complex evecij = complex( GSL_REAL(buf), GSL_IMAG(buf) );
      complex iEval(0.0,gsl_vector_get(d->fEval,j));
      fNuState[i] +=  exp(-iEval * kKm2eV*L) * nuComp[j] * evecij;
    }
  }

}

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().

{
  if (L.size()!=Ne.size()) abort();
  std::list<double>::const_iterator Li  (L.begin());
  std::list<double>::const_iterator Lend(L.end());
  std::list<double>::const_iterator Ni  (Ne.begin());
  for (; Li!=Lend; ++Li, ++Ni) {
    //cout << *Li << " km and " << *Ni << " mol/cm^3" << endl;
    this->PropMatter(*Li, E, *Ni, anti);
  }
}

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().

{
  int i;
  for (i=0; i<fNumNus; ++i){
    if (i==flv) fNuState[i] = one;
    else        fNuState[i] = zero;
  }
}

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().

                                     {

  // Do nothing if angle is zero
  if(fTheta[i][j]==0) return;

  double fSinBuffer = sin(fTheta[i][j]);
  double fCosBuffer = cos(fTheta[i][j]);

  double  fHmsBufferD;
  complex fHmsBufferC;

  // With Delta
  if(i+1<j){
    complex fExpBuffer = complex(cos(fDelta[i][j]), -sin(fDelta[i][j]));

    // General case
    if(i>0){
      // Top columns
      for(int k=0; k<i; k++){
        fHmsBufferC = fHms[k][i];

        fHms[k][i] *= fCosBuffer;
        fHms[k][i] += fHms[k][j] * fSinBuffer * conj(fExpBuffer);

        fHms[k][j] *= fCosBuffer;
        fHms[k][j] -= fHmsBufferC * fSinBuffer * fExpBuffer;
      }

      // Middle row and column
      for(int k=i+1; k<j; k++){
        fHmsBufferC = fHms[k][j];
    
        fHms[k][j] *= fCosBuffer;
        fHms[k][j] -= conj(fHms[i][k]) * fSinBuffer * fExpBuffer;

        fHms[i][k] *= fCosBuffer;
        fHms[i][k] += fSinBuffer * fExpBuffer * conj(fHmsBufferC);
      }

      // Nodes ij
      fHmsBufferC = fHms[i][i];
      fHmsBufferD = real(fHms[j][j]);

      fHms[i][i] *= fCosBuffer * fCosBuffer;
      fHms[i][i] += 2 * fSinBuffer * fCosBuffer * real(fHms[i][j] * conj(fExpBuffer));
      fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;

      fHms[j][j] *= fCosBuffer * fCosBuffer;
      fHms[j][j] += fSinBuffer * fHmsBufferC * fSinBuffer;
      fHms[j][j] -= 2 * fSinBuffer * fCosBuffer * real(fHms[i][j] * conj(fExpBuffer));
  
      fHms[i][j] -= 2 * fSinBuffer * real(fHms[i][j] * conj(fExpBuffer)) * fSinBuffer * fExpBuffer;
      fHms[i][j] += fSinBuffer * fCosBuffer * (fHmsBufferD - fHmsBufferC) * fExpBuffer;

    }
    // First rotation on j (No top columns)
    else{
      // Middle rows and columns
      for(int k=i+1; k<j; k++){
        fHms[k][j] = -conj(fHms[i][k]) * fSinBuffer * fExpBuffer;

        fHms[i][k] *= fCosBuffer;
      }

      // Nodes ij
      fHmsBufferD = real(fHms[i][i]);

      fHms[i][j] = fSinBuffer * fCosBuffer * (fHms[j][j] - fHmsBufferD) * fExpBuffer;

      fHms[i][i] *= fCosBuffer * fCosBuffer;
      fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;

      fHms[j][j] *= fCosBuffer * fCosBuffer;
      fHms[j][j] += fSinBuffer * fHmsBufferD * fSinBuffer;
    }
  
  }
  // Without Delta (No middle rows or columns: j = i+1)
  else{
    // General case
    if(i>0){
      // Top columns
      for(int k=0; k<i; k++){
        fHmsBufferC = fHms[k][i];

        fHms[k][i] *= fCosBuffer;
        fHms[k][i] += fHms[k][j] * fSinBuffer;

        fHms[k][j] *= fCosBuffer;
        fHms[k][j] -= fHmsBufferC * fSinBuffer;
      }

      // Nodes ij
      fHmsBufferC = fHms[i][i];
      fHmsBufferD = real(fHms[j][j]);

      fHms[i][i] *= fCosBuffer * fCosBuffer;
      fHms[i][i] += 2 * fSinBuffer * fCosBuffer * real(fHms[i][j]);
      fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;

      fHms[j][j] *= fCosBuffer * fCosBuffer;
      fHms[j][j] += fSinBuffer * fHmsBufferC * fSinBuffer;
      fHms[j][j] -= 2 * fSinBuffer * fCosBuffer * real(fHms[i][j]);
  
      fHms[i][j] -= 2 * fSinBuffer * real(fHms[i][j]) * fSinBuffer;
      fHms[i][j] += fSinBuffer * fCosBuffer * (fHmsBufferD - fHmsBufferC);

    }
    // First rotation (theta12)
    else{

      fHms[i][j] = fSinBuffer * fCosBuffer * fHms[j][j];

      fHms[i][i] = fSinBuffer * fHms[j][j] * fSinBuffer;

      fHms[j][j] *= fCosBuffer * fCosBuffer;
  
    }
  }

}

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().

{

  if(i>j){
    cout << "First argument should be smaller than second argument" << endl;
    cout << "Setting reverse order (Theta" << j << i << "). " << endl;
    int temp = i;
    i = j;
    j = temp;
  }
  if(i<1 || i>fNumNus-1 || j<2 || j>fNumNus){
    cout << "Theta" << i << j << " not valid for " << fNumNus;
    cout << " neutrinos. Doing nothing." << endl;
    return;
  }

  fTheta[i-1][j-1] = th;

  fBuiltHms = false;

}

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().

{

  if(i>j){
    cout << "First argument should be smaller than second argument" << endl;
    cout << "Setting reverse order (Delta" << j << i << "). " << endl;
    int temp = i;
    i = j;
    j = temp;
  }
  if(i<1 || i>fNumNus-1 || j<2 || j>fNumNus){
    cout << "Delta" << i << j << " not valid for " << fNumNus;
    cout << " neutrinos. Doing nothing." << endl;
    return;
  }
  if(i+1==j){
    cout << "Rotation " << i << j << " is real. Doing nothing." << endl;
    return;
  }

  fDelta[i-1][j-1] = delta;

  fBuiltHms = false;

}

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().

{

  if(j<2 || j>fNumNus){
    cout << "Dm" << j << "1 not valid for " << fNumNus;
    cout << " neutrinos. Doing nothing." << endl;
    return;
  }

  fDm[j-1] = dm;

  fBuiltHms = false;

}

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().

{

  this->InitializeVectors();

  if(fNumNus>2){
    this->SetAngle(1,2,0.6);
    this->SetAngle(1,3,0.16);
    this->SetAngle(2,3,0.7);
    this->SetDm(2,7.6e-5);
    this->SetDm(3,2.4e-3);
  }
  else if(fNumNus==2){
    this->SetAngle(1,2,0.7);
    this->SetDm(2,2.4e-3);
  }
  
}

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().

{

  // Check if anything has changed before recalculating
  if(Ne!=fCachedNe || E!=fCachedE || anti!=fCachedAnti || !fBuiltHms ){
    fCachedNe = Ne;
    fCachedE = E;
    fCachedAnti = anti;
    this->BuildHms();
  }
  else return;

  double lv = 2 * kGeV2eV*E;          // 2E in eV 
  double kr2GNe = kK2*M_SQRT2*kGf*Ne; // Matter potential in eV

  // Finish building Hamiltonian in matter with dimension of eV

  for(size_t i=0;i<size_t(fNumNus);i++){
    complex buf = fHms[i][i]/lv;
    *gsl_matrix_complex_ptr(d->H_GSL, i, i) = gsl_complex_rect(real(buf),0);
    for(size_t j=i+1;j<size_t(fNumNus);j++){
      buf = fHms[i][j]/lv;
      if(anti>0) *gsl_matrix_complex_ptr(d->H_GSL, j, i) = gsl_complex_rect(real(buf),-imag(buf));
      else       *gsl_matrix_complex_ptr(d->H_GSL, j, i) = gsl_complex_rect(real(buf), imag(buf));
    }
    if(i>2){
      // Subtract NC coherent forward scattering from sterile neutrinos. See arXiv:hep-ph/0606054v3, eq. 3.15, for example. 
      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);
      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);;
    }
  }
  // Add nue CC coherent forward scattering from sterile neutrinos. 
  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);
  else       *gsl_matrix_complex_ptr(d->H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(d->H_GSL,0,0) , -kr2GNe);

  // Solve Hamiltonian for eigensystem
  gsl_eigen_hermv(d->H_GSL, d->fEval, d->fEvec, d->W_GSL);
  
}

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

Definition at line 82 of file PMNS_Sterile.h.

Referenced by BuildHms(), GetNFlavors(), InitializeVectors(), P(), PMNS_Sterile(), PropMatter(), ResetToFlavour(), SetAngle(), SetDelta(), SetDm(), SetStdPars(), and SolveHam().

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:

• OscLib/PMNS_Sterile.h
• OscLib/PMNS_Sterile.cxx