OscProb::PMNS_Base Class Reference

#include <PMNS_Base.h>

Inheritance diagram for OscProb::PMNS_Base:
OscProb::PMNS_Fast OscProb::PMNS_Sterile OscProb::PMNS_Transitions OscProb::PMNS_NSI

List of all members.

Public Member Functions

 PMNS_Base (int numNus=3)
virtual ~PMNS_Base ()
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.
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 void SetLength (double L)
virtual void SetEnergy (double E)
virtual void SetDensity (double Ne)
virtual void SetIsNuBar (bool isNuBar)
virtual double Prob (int flvi, int flvf)
virtual double Prob (int flvi, int flvf, double E)
virtual double Prob (int flvi, int flvf, double E, double L)
virtual double AvgProb (int flvi, int flvf, double E, double dE=0)
virtual double AvgProbLoE (int flvi, int flvf, double LoE, double dLoE=0)

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.
virtual void BuildHms ()
virtual void SolveHam ()
 Solve the full Hamiltonian for eigenvectors and eigenvalues.
virtual bool SamePars (double E, double Ne, bool isNuBar)

Protected Attributes

std::complex< double > zero
std::complex< double > one
const double kKm2eV
 km to eV^-1
const double kK2
 mole/GeV^2/cm^3 to eV
const double kGeV2eV
 GeV to eV.
const double kGf
int fNumNus
std::vector< double > fDm
 m^2_i - m^2_1 in vacuum
std::vector< std::vector
< double > > 
fTheta
 theta[i][j] mixing angle
std::vector< std::vector
< double > > 
fDelta
 delta[i][j] CP violating phase
std::vector< complexfNuState
 The neutrino current state.
std::vector< std::vector
< complex > > 
fHms
 matrix H*2E in eV^2
std::vector< double > fEval
 Eigenvalues of the Hamiltonian.
std::vector< std::vector
< complex > > 
fEvec
 Eigenvectors of the Hamiltonian.
double fLength
 Path length in km.
double fNe
 Path density in mol/cm^3.
double fEnergy
 Neutrino energy in GeV.
bool fIsNuBar
 Nu-NuBar flag.
bool fBuiltHms
 Tag to avoid rebuilding Hms.
bool fGotES
 Tag to avoid recalculating eigensystem.

Detailed Description

Definition at line 18 of file PMNS_Base.h.


Member Typedef Documentation

typedef std::complex<double> OscProb::PMNS_Base::complex [protected]

Definition at line 88 of file PMNS_Base.h.


Constructor & Destructor Documentation

PMNS_Base::PMNS_Base ( int  numNus = 3  ) 

Definition at line 31 of file PMNS_Base.cxx.

References fNumNus, ResetToFlavour(), and SetStdPars().

00031                                : 
00032 zero(0,0), one(1,0),
00033 kKm2eV( 5.06773103202e+09 ), 
00034 kK2( 4.62711492217e-09 ),
00035 kGeV2eV( 1.0e+09 ),
00036 kGf( 1.166371E-5 )
00037 {
00038 
00039   fNumNus = numNus;
00040   this->SetStdPars();
00041   this->ResetToFlavour(1);
00042 
00043 }

PMNS_Base::~PMNS_Base (  )  [virtual]

Definition at line 45 of file PMNS_Base.cxx.

00045 {}


Member Function Documentation

double PMNS_Base::AvgProb ( int  flvi,
int  flvf,
double  E,
double  dE = 0 
) [virtual]

Return the probability of flvi going to flvf averaged over a bin of width dE

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

Get average oscillation probability over an energy range. This gets transformed into L/E, since the oscillation terms have arguments linear in L/E and not E.

Parameters:
flvi - initial flavor (0,1,2) for (nue, numu, nutau)
flvf - final flavor (0,1,2) for (nue, numu, nutau)
E - Energy center in GeV
dE - Energy width in GeV

Definition at line 513 of file PMNS_Base.cxx.

References AvgProbLoE(), fLength, and Prob().

00514 {
00515 
00516   // Do nothing if energy is not positive
00517   if(E<=0) return 0;
00518 
00519   // Don't average zero width
00520   if(dE<=0) return Prob(flvi, flvf, E);
00521 
00522   // Define L/E variables
00523   double LoE = 0;
00524   double dLoE = 0;
00525   
00526   // Set a minimum energy
00527   double minE = 0.1 * E;
00528   
00529   // Transform range to L/E
00530   // Full range if low edge > minE
00531   if(E-dE/2 > minE){
00532     LoE = 0.5 * (fLength/(E-dE/2) + fLength/(E+dE/2));
00533     dLoE = fLength/(E-dE/2) - fLength/(E+dE/2);
00534   }
00535   // Else start at minE
00536   else{
00537     LoE = 0.5 * (fLength/minE + fLength/(E+dE/2));
00538     dLoE = fLength/minE - fLength/(E+dE/2);
00539   }
00540 
00541   // Compute average in LoE
00542   return AvgProbLoE(flvi, flvf, LoE, dLoE);
00543 
00544 } 

double PMNS_Base::AvgProbLoE ( int  flvi,
int  flvf,
double  LoE,
double  dLoE = 0 
) [virtual]

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

Get average oscillation probability over an L/E range. The probabilities are weighted by (L/E)^-2 so that event density is flat in energy. This avoids giving too much weight to low energies. Better approximations would be achieved if we used an interpolated event density.

Parameters:
flvi - initial flavor (0,1,2) for (nue, numu, nutau)
flvf - final flavor (0,1,2) for (nue, numu, nutau)
LoE - L/E center in km/GeV
dLoE - L/E width in km/GeV

Definition at line 559 of file PMNS_Base.cxx.

References fEnergy, fEval, fGotES, fIsNuBar, fLength, fNe, fNumNus, k1267, kGeV2eV, kKm2eV, Prob(), SamePars(), SetEnergy(), SolveHam(), and TMath::Sort().

Referenced by AvgProb().

00560 {
00561 
00562   // Don't average zero width
00563   if(dLoE<=0) return Prob(flvi, flvf, fLength / LoE);
00564 
00565   // Check if energy changed from last computation
00566   if(!SamePars(fLength / LoE, fNe, fIsNuBar)){
00567     SetEnergy(fLength / LoE);
00568   }
00569 
00570   // If Hamiltonian changed, solve for eigensystem
00571   if(!fGotES) this->SolveHam();
00572 
00573   // Define conversion factor [km/GeV -> 1/(4 eV^2)]
00574   const double k1267 = kKm2eV / (4 * kGeV2eV);
00575   
00576   // Get list of all effective Dm^2
00577   vector<double> effDm;
00578   
00579   for(int i=0; i<fNumNus-1; i++){
00580     for(int j=i+1; j<fNumNus; j++){
00581       effDm.push_back( 2 * kGeV2eV * fEnergy * TMath::Abs(fEval[j] - fEval[i]) );
00582     }
00583   }
00584   
00585   int numDm = effDm.size();
00586   
00587   // Set a number of sub-divisions to achieve "good" accuracy
00588   // This needs to be studied better
00589   int n_div = TMath::Ceil( 20 * pow(dLoE/LoE,0.8) );
00590 
00591   double sumw = 0;
00592   double prob = 0;
00593   
00594   // Loop over sub-divisions
00595   for(int k=0; k<n_div; k++){
00596 
00597     // Define sub-division center and width
00598     double bctr = LoE - dLoE/2 + (k+0.5)*dLoE/n_div;
00599     double bwdt = dLoE/n_div;
00600     
00601     // Make a list of indices sorted by Dm^2 value
00602     vector<int> Idx(numDm, 0);
00603     TMath::Sort(numDm, &effDm[0], &Idx[0]);
00604 
00605     // Make a vector of L/E sample values
00606     // Initialized in the sub-division center
00607     vector<double> samples;
00608     samples.push_back(bctr);
00609 
00610     // Loop over all Dm^2 to average each frequency
00611     // This will recursively sample points in smaller
00612     // bins so that all relevant frequencies are used
00613     for(int i=0; i<numDm; i++){
00614       
00615       // Copy the list of sample L/E values
00616       vector<double> prev = samples;
00617 
00618       // Redefine bin width to lie within full sub-division
00619       double Width = 2*TMath::Min(prev[0] - (bctr - bwdt/2), (bctr + bwdt/2) - prev[0]);
00620       
00621       // Compute oscillation argument sorted from lowest  to highest
00622       const double arg = k1267 * effDm[Idx[i]] * Width;
00623     
00624       // Skip small oscillation values unless it's the last one
00625       if(arg<0.9 && i!=numDm-1) continue;
00626     
00627       // Reset samples to redefine them
00628       samples.clear();
00629 
00630       // Loop over previous samples
00631       for(int j=0; j<int(prev.size()); j++){
00632       
00633         // Compute new sample points around old samples
00634         // This is based on a oscillatory quadrature rule 
00635         double sample = (1/sqrt(3)) * (Width/2);
00636         if(arg!=0) sample = TMath::ACos(TMath::Sin(arg)/arg)/arg * (Width/2);
00637       
00638         // Add samples above and below center
00639         samples.push_back(prev[j]-sample);
00640         samples.push_back(prev[j]+sample);
00641 
00642       }
00643 
00644     }// End of loop over Dm^2
00645 
00646     // Loop over all sample points
00647     for(int j=0; j<int(samples.size()); j++){
00648 
00649       // Set (L/E)^-2 weights
00650       double w = 1./pow(samples[j],2);
00651 
00652       // Add weighted probability
00653       prob += w * Prob(flvi, flvf, fLength / samples[j]);
00654       
00655       // Increment sum of weights
00656       sumw += w;
00657       
00658     }
00659     
00660   }// End of loop over sub-divisions
00661 
00662   // Return weighted average of probabilities
00663   return prob / sumw;
00664 
00665 }

void PMNS_Base::BuildHms (  )  [protected, virtual]

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.

Definition at line 313 of file PMNS_Base.cxx.

References fBuiltHms, fDm, fHms, fNumNus, and RotateH().

Referenced by OscProb::PMNS_Sterile::SolveHam(), OscProb::PMNS_Fast::SolveHam(), OscProb::PMNS_NSI::SolveHam(), and OscProb::PMNS_Transitions::SolveHam().

00314 {
00315 
00316   // Check if anything changed
00317   if(fBuiltHms) return;
00318 
00319   for(int j=0; j<fNumNus; j++){
00320     // Set mass splitting
00321     fHms[j][j] = fDm[j];
00322     // Reset off-diagonal elements
00323     for(int i=0; i<j; i++){
00324       fHms[i][j] = 0;
00325     }
00326     // Rotate j neutrinos
00327     for(int i=0; i<j; i++){
00328       RotateH(i,j);
00329     }
00330   }
00331   
00332   // Tag as built
00333   fBuiltHms = true;
00334 
00335 }

void PMNS_Base::InitializeVectors (  )  [protected, virtual]

Set vector sizes and initialize elements to zero

Definition at line 51 of file PMNS_Base.cxx.

References fDelta, fDm, fEval, fEvec, fHms, fNumNus, fNuState, fTheta, and zero.

Referenced by SetStdPars().

00052 {
00053 
00054   fDm    = vector<double>(fNumNus, 0);
00055   fTheta = vector< vector<double> >(fNumNus, vector<double>(fNumNus,0));
00056   fDelta = vector< vector<double> >(fNumNus, vector<double>(fNumNus,0));
00057 
00058   fNuState = vector<complex>(fNumNus, zero);
00059   fHms     = vector< vector<complex> >(fNumNus, vector<complex>(fNumNus,zero));
00060 
00061   fEval = vector<double>(fNumNus, 0);
00062   fEvec = vector< vector<complex> >(fNumNus, vector<complex>(fNumNus,zero));
00063 
00064 }

double PMNS_Base::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 446 of file PMNS_Base.cxx.

References fNumNus, and fNuState.

Referenced by OscCalc::OscillateLSND(), OscCalc::OscillateNSI(), OscFit::OscillationCalculator::OscProb(), Prob(), NuOscProbCalc::ThreeFlavourExactNuMuToNuElectron(), NuOscProbCalc::ThreeFlavourExactNuMuToNuMu(), and NuOscProbCalc::ThreeFlavourExactNuMuToNuTau().

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

double PMNS_Base::Prob ( int  flvi,
int  flvf,
double  E,
double  L 
) [virtual]

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

Get oscillation probability

Parameters:
flvi - initial flavor (0,1,2) for (nue, numu, nutau)
flvf - final flavor (0,1,2) for (nue, numu, nutau)
E - Energy in GeV
L - Path length in km

Definition at line 492 of file PMNS_Base.cxx.

References Prob(), SetEnergy(), and SetLength().

00493 {
00494 
00495   SetEnergy(E);
00496   SetLength(L);
00497 
00498   return Prob(flvi, flvf);
00499 
00500 }

double PMNS_Base::Prob ( int  flvi,
int  flvf,
double  E 
) [virtual]

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

Get oscillation probability

Parameters:
flvi - initial flavor (0,1,2) for (nue, numu, nutau)
flvf - final flavor (0,1,2) for (nue, numu, nutau)
E - Energy in GeV

Definition at line 475 of file PMNS_Base.cxx.

References Prob(), and SetEnergy().

00476 {
00477 
00478   SetEnergy(E);
00479 
00480   return Prob(flvi, flvf);
00481 
00482 }

double PMNS_Base::Prob ( int  flvi,
int  flvf 
) [virtual]

Give a list of path lengths and densities Not implemented yet Return the probability of flvi going to flvf

Parameters:
flvi - initial flavor (0,1,2) = (nue,numu,nutau)
flvf - final flavor (0,1,2) = (nue,numu,nutau)
E - Neutrino energy in GeV
L - Path lentgh in km

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

Get oscillation probability

Parameters:
flvi - initial flavor (0,1,2) for (nue, numu, nutau)
flvf - final flavor (0,1,2) for (nue, numu, nutau)

Definition at line 458 of file PMNS_Base.cxx.

References fEnergy, fIsNuBar, fLength, fNe, P(), PropMatter(), and ResetToFlavour().

Referenced by AvgProb(), AvgProbLoE(), and Prob().

00459 {
00460 
00461   ResetToFlavour(flvi);
00462   PropMatter(fLength, fEnergy, fNe, fIsNuBar ? -1 : 1);
00463 
00464   return P(flvf);
00465 
00466 }

void PMNS_Base::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 412 of file PMNS_Base.cxx.

References PropMatter().

00416 {
00417   if (L.size()!=Ne.size()) abort();
00418   std::list<double>::const_iterator Li  (L.begin());
00419   std::list<double>::const_iterator Lend(L.end());
00420   std::list<double>::const_iterator Ni  (Ne.begin());
00421   for (; Li!=Lend; ++Li, ++Ni) {
00422     PropMatter(*Li, E, *Ni, anti);
00423   }
00424 }

void PMNS_Base::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 - positive = neutrino case, negative = 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 mol/cm^3.

Parameters:
anti - positive = neutrino case, negative = anti-neutrino case

Definition at line 374 of file PMNS_Base.cxx.

References fEval, fEvec, fGotES, fNumNus, fNuState, kKm2eV, SamePars(), SetDensity(), SetEnergy(), SetIsNuBar(), SolveHam(), and zero.

Referenced by OscCalc::OscillateLSND(), OscCalc::OscillateNSI(), Prob(), OscFit::OscillationCalculator::PropMatter(), PropMatter(), OscProb::PMNS_Fast::PropVacuum(), NuOscProbCalc::ThreeFlavourExactNuMuToNuElectron(), NuOscProbCalc::ThreeFlavourExactNuMuToNuMu(), and NuOscProbCalc::ThreeFlavourExactNuMuToNuTau().

00375 {
00376 
00377   // Check if anything changed from last computation
00378   // If it did, set parameters to new values
00379   if(!SamePars(E, Ne, anti <= 0)){
00380     SetEnergy(E);
00381     SetDensity(Ne);
00382     SetIsNuBar(anti <= 0);
00383   }
00384 
00385   // If Hamiltonian changed, solve for eigensystem
00386   if(!fGotES) this->SolveHam();
00387 
00388   // Store coefficients of propagation eigenstates
00389   vector<complex> nuComp(fNumNus, zero);
00390   for(int i=0;i<fNumNus;i++){
00391     nuComp[i] = 0;
00392     for(int j=0;j<fNumNus;j++){
00393       nuComp[i] += fNuState[j] * conj(fEvec[j][i]);
00394     }
00395   }
00396 
00397   // Propagate neutrino state
00398   for(int i=0;i<fNumNus;i++){
00399     fNuState[i] = 0;
00400     for(int j=0;j<fNumNus;j++){
00401       double arg = fEval[j] * kKm2eV*L;
00402       fNuState[i] +=  complex(cos(arg),-sin(arg)) * nuComp[j] * fEvec[i][j];
00403     }
00404   }
00405 
00406 }

void PMNS_Base::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 431 of file PMNS_Base.cxx.

References fNumNus, fNuState, one, and zero.

Referenced by OscCalc::OscillateLSND(), OscCalc::OscillateNSI(), PMNS_Base(), Prob(), OscFit::OscillationCalculator::ResetWithNewFlavour(), OscFit::OscillationCalculator::SetPMNS(), NuOscProbCalc::ThreeFlavourExactNuMuToNuElectron(), NuOscProbCalc::ThreeFlavourExactNuMuToNuMu(), and NuOscProbCalc::ThreeFlavourExactNuMuToNuTau().

00432 {
00433   register int i;
00434   for (i=0; i<fNumNus; ++i){
00435     if (i==flv) fNuState[i] = one;
00436     else        fNuState[i] = zero;
00437   }
00438 }

void PMNS_Base::RotateH ( int  i,
int  j 
) [protected, virtual]

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 179 of file PMNS_Base.cxx.

References fDelta, fHms, and fTheta.

Referenced by BuildHms().

00179                                   {
00180 
00181   // Do nothing if angle is zero
00182   if(fTheta[i][j]==0) return;
00183 
00184   double fSinBuffer = sin(fTheta[i][j]);
00185   double fCosBuffer = cos(fTheta[i][j]);
00186 
00187   double  fHmsBufferD;
00188   complex fHmsBufferC;
00189 
00190   // With Delta
00191   if(i+1<j){
00192     complex fExpBuffer = complex(cos(fDelta[i][j]), -sin(fDelta[i][j]));
00193 
00194     // General case
00195     if(i>0){
00196       // Top columns
00197       for(int k=0; k<i; k++){
00198         fHmsBufferC = fHms[k][i];
00199 
00200         fHms[k][i] *= fCosBuffer;
00201         fHms[k][i] += fHms[k][j] * fSinBuffer * conj(fExpBuffer);
00202 
00203         fHms[k][j] *= fCosBuffer;
00204         fHms[k][j] -= fHmsBufferC * fSinBuffer * fExpBuffer;
00205       }
00206 
00207       // Middle row and column
00208       for(int k=i+1; k<j; k++){
00209         fHmsBufferC = fHms[k][j];
00210     
00211         fHms[k][j] *= fCosBuffer;
00212         fHms[k][j] -= conj(fHms[i][k]) * fSinBuffer * fExpBuffer;
00213 
00214         fHms[i][k] *= fCosBuffer;
00215         fHms[i][k] += fSinBuffer * fExpBuffer * conj(fHmsBufferC);
00216       }
00217 
00218       // Nodes ij
00219       fHmsBufferC = fHms[i][i];
00220       fHmsBufferD = real(fHms[j][j]);
00221 
00222       fHms[i][i] *= fCosBuffer * fCosBuffer;
00223       fHms[i][i] += 2 * fSinBuffer * fCosBuffer * real(fHms[i][j] * conj(fExpBuffer));
00224       fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;
00225 
00226       fHms[j][j] *= fCosBuffer * fCosBuffer;
00227       fHms[j][j] += fSinBuffer * fHmsBufferC * fSinBuffer;
00228       fHms[j][j] -= 2 * fSinBuffer * fCosBuffer * real(fHms[i][j] * conj(fExpBuffer));
00229   
00230       fHms[i][j] -= 2 * fSinBuffer * real(fHms[i][j] * conj(fExpBuffer)) * fSinBuffer * fExpBuffer;
00231       fHms[i][j] += fSinBuffer * fCosBuffer * (fHmsBufferD - fHmsBufferC) * fExpBuffer;
00232 
00233     }
00234     // First rotation on j (No top columns)
00235     else{
00236       // Middle rows and columns
00237       for(int k=i+1; k<j; k++){
00238         fHms[k][j] = -conj(fHms[i][k]) * fSinBuffer * fExpBuffer;
00239 
00240         fHms[i][k] *= fCosBuffer;
00241       }
00242 
00243       // Nodes ij
00244       fHmsBufferD = real(fHms[i][i]);
00245 
00246       fHms[i][j] = fSinBuffer * fCosBuffer * (fHms[j][j] - fHmsBufferD) * fExpBuffer;
00247 
00248       fHms[i][i] *= fCosBuffer * fCosBuffer;
00249       fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;
00250 
00251       fHms[j][j] *= fCosBuffer * fCosBuffer;
00252       fHms[j][j] += fSinBuffer * fHmsBufferD * fSinBuffer;
00253     }
00254   
00255   }
00256   // Without Delta (No middle rows or columns: j = i+1)
00257   else{
00258     // General case
00259     if(i>0){
00260       // Top columns
00261       for(int k=0; k<i; k++){
00262         fHmsBufferC = fHms[k][i];
00263 
00264         fHms[k][i] *= fCosBuffer;
00265         fHms[k][i] += fHms[k][j] * fSinBuffer;
00266 
00267         fHms[k][j] *= fCosBuffer;
00268         fHms[k][j] -= fHmsBufferC * fSinBuffer;
00269       }
00270 
00271       // Nodes ij
00272       fHmsBufferC = fHms[i][i];
00273       fHmsBufferD = real(fHms[j][j]);
00274 
00275       fHms[i][i] *= fCosBuffer * fCosBuffer;
00276       fHms[i][i] += 2 * fSinBuffer * fCosBuffer * real(fHms[i][j]);
00277       fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;
00278 
00279       fHms[j][j] *= fCosBuffer * fCosBuffer;
00280       fHms[j][j] += fSinBuffer * fHmsBufferC * fSinBuffer;
00281       fHms[j][j] -= 2 * fSinBuffer * fCosBuffer * real(fHms[i][j]);
00282   
00283       fHms[i][j] -= 2 * fSinBuffer * real(fHms[i][j]) * fSinBuffer;
00284       fHms[i][j] += fSinBuffer * fCosBuffer * (fHmsBufferD - fHmsBufferC);
00285 
00286     }
00287     // First rotation (theta12)
00288     else{
00289 
00290       fHms[i][j] = fSinBuffer * fCosBuffer * fHms[j][j];
00291 
00292       fHms[i][i] = fSinBuffer * fHms[j][j] * fSinBuffer;
00293 
00294       fHms[j][j] *= fCosBuffer * fCosBuffer;
00295   
00296     }
00297   }
00298 
00299 }

bool PMNS_Base::SamePars ( double  E,
double  Ne,
bool  isNuBar 
) [protected, virtual]

Check if parameters are the same as stored

Definition at line 341 of file PMNS_Base.cxx.

References fBuiltHms, fEnergy, fIsNuBar, and fNe.

Referenced by AvgProbLoE(), and PropMatter().

00342 {
00343 
00344   // Check if energy is the same
00345   bool isSame = (E == fEnergy);
00346   
00347   // Check if nu-nubar tag is the same
00348   isSame *= (isNuBar == fIsNuBar);
00349   
00350   // Check if density is the same
00351   isSame *= (Ne == fNe);
00352   
00353   // Check if mass and mixing are the same
00354   isSame *= fBuiltHms;
00355   
00356   return isSame;
00357 
00358 }

void PMNS_Base::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 101 of file PMNS_Base.cxx.

References fBuiltHms, fNumNus, and fTheta.

Referenced by OscProb::PMNS_Fast::SetMix(), SetStdPars(), and OscCalc::UpdateLSND().

00102 {
00103 
00104   if(i>j){
00105     cout << "First argument should be smaller than second argument" << endl;
00106     cout << "Setting reverse order (Theta" << j << i << "). " << endl;
00107     int temp = i;
00108     i = j;
00109     j = temp;
00110   }
00111   if(i<1 || i>fNumNus-1 || j<2 || j>fNumNus){
00112     cout << "Theta" << i << j << " not valid for " << fNumNus;
00113     cout << " neutrinos. Doing nothing." << endl;
00114     return;
00115   }
00116 
00117   fTheta[i-1][j-1] = th;
00118 
00119   fBuiltHms = false;
00120 
00121 }

void PMNS_Base::SetDelta ( int  i,
int  j,
double  delta 
) [virtual]

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

Definition at line 127 of file PMNS_Base.cxx.

References fBuiltHms, fDelta, and fNumNus.

Referenced by OscProb::PMNS_Fast::SetMix(), and OscCalc::UpdateLSND().

00128 {
00129 
00130   if(i>j){
00131     cout << "First argument should be smaller than second argument" << endl;
00132     cout << "Setting reverse order (Delta" << j << i << "). " << endl;
00133     int temp = i;
00134     i = j;
00135     j = temp;
00136   }
00137   if(i<1 || i>fNumNus-1 || j<2 || j>fNumNus){
00138     cout << "Delta" << i << j << " not valid for " << fNumNus;
00139     cout << " neutrinos. Doing nothing." << endl;
00140     return;
00141   }
00142   if(i+1==j){
00143     cout << "Rotation " << i << j << " is real. Doing nothing." << endl;
00144     return;
00145   }
00146 
00147   fDelta[i-1][j-1] = delta;
00148 
00149   fBuiltHms = false;
00150 
00151 }

virtual void OscProb::PMNS_Base::SetDensity ( double  Ne  )  [inline, virtual]

Definition at line 63 of file PMNS_Base.h.

References fGotES, and fNe.

Referenced by PropMatter(), and SetStdPars().

00063 { fNe = Ne; fGotES = false; }

void PMNS_Base::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 157 of file PMNS_Base.cxx.

References fBuiltHms, fDm, and fNumNus.

Referenced by OscProb::PMNS_Fast::SetDeltaMsqrs(), SetStdPars(), and OscCalc::UpdateLSND().

00158 {
00159 
00160   if(j<2 || j>fNumNus){
00161     cout << "Dm" << j << "1 not valid for " << fNumNus;
00162     cout << " neutrinos. Doing nothing." << endl;
00163     return;
00164   }
00165 
00166   fDm[j-1] = dm;
00167 
00168   fBuiltHms = false;
00169 
00170 }

virtual void OscProb::PMNS_Base::SetEnergy ( double  E  )  [inline, virtual]

Definition at line 62 of file PMNS_Base.h.

References fEnergy, and fGotES.

Referenced by AvgProbLoE(), Prob(), PropMatter(), and SetStdPars().

00062 { fEnergy = E; fGotES = false; }

virtual void OscProb::PMNS_Base::SetIsNuBar ( bool  isNuBar  )  [inline, virtual]

Definition at line 64 of file PMNS_Base.h.

References fGotES, and fIsNuBar.

Referenced by PropMatter(), and SetStdPars().

00064 { fIsNuBar = isNuBar; fGotES = false; }

virtual void OscProb::PMNS_Base::SetLength ( double  L  )  [inline, virtual]

Set oscillation

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

Definition at line 61 of file PMNS_Base.h.

References fLength.

Referenced by Prob(), and SetStdPars().

00061 { fLength = L; }

void PMNS_Base::SetStdPars (  )  [virtual]

Set standard 3-flavor parameters.

Set standard oscillation parameters from PDG. Length from MINOS FD with approximate peak energy Density is approximate from CRUST2.0 (~2.8 g/cm^3 = ~1.4 mol/cm^3)

Definition at line 73 of file PMNS_Base.cxx.

References MuELoss::e, fNumNus, InitializeVectors(), SetAngle(), SetDensity(), SetDm(), SetEnergy(), SetIsNuBar(), and SetLength().

Referenced by PMNS_Base().

00074 {
00075 
00076   SetDensity(1.4);
00077   SetLength(735);
00078   SetEnergy(1.5);
00079   SetIsNuBar(false);
00080 
00081   InitializeVectors();
00082 
00083   if(fNumNus>2){
00084     SetAngle(1,2, 0.588);
00085     SetAngle(1,3, 0.154);                                          
00086     SetAngle(2,3, 0.722);
00087     SetDm(2, 7.54e-5);
00088     SetDm(3, 2.47e-3);
00089   }
00090   else if(fNumNus==2){
00091     SetAngle(1,2, 0.712);
00092     SetDm(2, 2.43e-3);
00093   }
00094   
00095 }

void PMNS_Base::SolveHam (  )  [protected, virtual]

Solve the full Hamiltonian for eigenvectors and eigenvalues.

Solve the full Hamiltonian Not implemented in base class.

Reimplemented in OscProb::PMNS_Fast, OscProb::PMNS_NSI, OscProb::PMNS_Sterile, and OscProb::PMNS_Transitions.

Definition at line 365 of file PMNS_Base.cxx.

Referenced by AvgProbLoE(), and PropMatter().

00365 {}


Member Data Documentation

std::vector< std::vector<double> > OscProb::PMNS_Base::fDelta [protected]

delta[i][j] CP violating phase

Definition at line 120 of file PMNS_Base.h.

Referenced by InitializeVectors(), RotateH(), SetDelta(), and OscProb::PMNS_Fast::SetVacuumEigensystem().

std::vector<double> OscProb::PMNS_Base::fDm [protected]

m^2_i - m^2_1 in vacuum

Definition at line 118 of file PMNS_Base.h.

Referenced by BuildHms(), InitializeVectors(), SetDm(), and OscProb::PMNS_Fast::SetVacuumEigensystem().

double OscProb::PMNS_Base::fEnergy [protected]
std::vector<double> OscProb::PMNS_Base::fEval [protected]
std::vector< std::vector<complex> > OscProb::PMNS_Base::fEvec [protected]
bool OscProb::PMNS_Base::fGotES [protected]
std::vector< std::vector<complex> > OscProb::PMNS_Base::fHms [protected]
bool OscProb::PMNS_Base::fIsNuBar [protected]
double OscProb::PMNS_Base::fLength [protected]

Path length in km.

Definition at line 128 of file PMNS_Base.h.

Referenced by AvgProb(), AvgProbLoE(), Prob(), and SetLength().

double OscProb::PMNS_Base::fNe [protected]
int OscProb::PMNS_Base::fNumNus [protected]
std::vector<complex> OscProb::PMNS_Base::fNuState [protected]

The neutrino current state.

Definition at line 122 of file PMNS_Base.h.

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

std::vector< std::vector<double> > OscProb::PMNS_Base::fTheta [protected]

theta[i][j] mixing angle

Definition at line 119 of file PMNS_Base.h.

Referenced by InitializeVectors(), RotateH(), SetAngle(), and OscProb::PMNS_Fast::SetVacuumEigensystem().

const double OscProb::PMNS_Base::kGeV2eV [protected]
const double OscProb::PMNS_Base::kGf [protected]
const double OscProb::PMNS_Base::kK2 [protected]

mole/GeV^2/cm^3 to eV

Definition at line 96 of file PMNS_Base.h.

Referenced by OscProb::PMNS_Sterile::SolveHam(), OscProb::PMNS_NSI::SolveHam(), and OscProb::PMNS_Fast::SolveHam().

const double OscProb::PMNS_Base::kKm2eV [protected]

km to eV^-1

Definition at line 95 of file PMNS_Base.h.

Referenced by AvgProbLoE(), and PropMatter().

std::complex<double> OscProb::PMNS_Base::one [protected]

Definition at line 92 of file PMNS_Base.h.

Referenced by ResetToFlavour().

std::complex<double> OscProb::PMNS_Base::zero [protected]

Definition at line 91 of file PMNS_Base.h.

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


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

Generated on 3 Oct 2018 for loon by  doxygen 1.6.1