NuOscProbCalc Namespace Reference

Classes

class  SterileProb

Functions

double AverageOscillationWeight (double energy, double dm2, double sn2, double binwidth=0)
double OscillationWeight (double energy, double dm2, double sn2)
double OscillationWeightMu2Mu (double energy, double dm2, double sn2, double g_mu, double g_tau, double c_mumu, double c_tautau)
double OscillationWeightMu2MuBar (double energy, double dm2, double sn2, double g_mu, double g_tau, double c_mumu, double c_tautau)
double OscillationWeightMu2Tau (double energy, double dm2, double sn2, double g_mu, double g_tau, double c_mumu, double c_tautau)
double OscillationWeightMu2TauBar (double energy, double dm2, double sn2, double g_mu, double g_tau, double c_mumu, double c_tautau)
double DecayWeightCC (double energy, double alpha, double sn2, double dm2=0, double binwidth=0)
double DecayWeightMuToTau (double energy, double alpha, double sn2, double dm2=0, double binwidth=0)
double DecayWeightNC (double energy, double alpha, double sn2)
double DecoherenceWeight (double energy, double mu2, double sn2, double dm2=0, int power=-1, double binwidth=0)
double NSIWeight (double energy, double dm2, double sn2, double epsilon, double sign)
double NSIWeightPrime (double energy, double dm2, double sn2, double epsilon, double epsprime, double sign)
double HighMassSterileWeight (double energy, double dm2, double sn2, double sn224)
double ThreeFlavourNuMuToNuMu (const double energy, double dm232, const double theta23, double dm221, const double theta12, const double theta13, double delta, const double baseline, const int inu)
double ThreeFlavourNuMuToNuTau (const double energy, double dm232, const double theta23, double dm221, const double theta12, const double theta13, double delta, const double baseline, const int inu)
double ThreeFlavourNuMuToNuElectron (const double energy, double dm232, const double theta23, double dm221, const double theta12, const double theta13, double delta, const double baseline, const int inu)
double ThreeFlavourExactNuMuToNuMu (const double energy, double dm232, const double theta23, double dm221, const double theta12, const double theta13, double delta, const double baseline, const int inu)
double ThreeFlavourExactNuMuToNuTau (const double energy, double dm232, const double theta23, double dm221, const double theta12, const double theta13, double delta, const double baseline, const int inu)
double ThreeFlavourExactNuMuToNuElectron (const double energy, double dm232, const double theta23, double dm221, const double theta12, const double theta13, double delta, const double baseline, const int inu)
double FourFlavourDisappearanceWeight (const double energy, double dm232, const double theta23, double dm221, double dm243, const double theta12, const double theta13, const double theta14, const double theta24, const double theta34, double delta1, double delta2, double delta3, const double baseline)
double FourFlavourNuMuToNuSProbability (const double energy, double dm232, const double theta23, double dm221, double dm243, const double theta12, const double theta13, const double theta14, const double theta24, const double theta34, double delta1, double delta2, double delta3, const double baseline)
double FourFlavourNuESurvivalProbability (const double energy, double dm232, const double theta23, double dm221, double dm243, const double theta12, const double theta13, const double theta14, const double theta24, const double theta34, double delta1, double delta2, double delta3, const double baseline)
double FourFlavourNuMuToNuEProbability (const double energy, double dm232, const double theta23, double dm221, double dm243, const double theta12, const double theta13, const double theta14, const double theta24, const double theta34, double delta1, double delta2, double delta3, const double baseline)
double FourFlavourNuMuToNuTauProbability (const double energy, double dm232, const double theta23, double dm221, double dm243, const double theta12, const double theta13, const double theta14, const double theta24, const double theta34, double delta1, double delta2, double delta3, const double baseline)
double LargeExtraDimensionProbability (int initialFlavor, int finalFlavor, double th12, double th13, double th23, double deltacp, double dm221, double dm232, double m0, double a, double energy, double baseLine)
double LargeExtraDimensionNuMuToNuSProbability (double th12, double th13, double th23, double deltacp, double dm221, double dm232, double m0, double a, double energy, double baseLine)
double LargeExtraDimensionDisappearanceWeight (double th12, double th13, double th23, double deltacp, double dm221, double dm232, double m0, double a, double energy, double baseLine)
double LargeExtraDimensionNuESurvivalProbability (double th12, double th13, double th23, double deltacp, double dm221, double dm232, double m0, double a, double energy, double baseLine)
double LargeExtraDimensionNuMuToNuEProbability (double th12, double th13, double th23, double deltacp, double dm221, double dm232, double m0, double a, double energy, double baseLine)
double LargeExtraDimensionNuMuToNuTauProbability (double th12, double th13, double th23, double deltacp, double dm221, double dm232, double m0, double a, double energy, double baseLine)

Function Documentation

double NuOscProbCalc::AverageOscillationWeight ( double  energy,
double  dm2,
double  sn2,
double  binwidth = 0 
)

Definition at line 218 of file NuOscProbCalc.cxx.

References k1267, kFDDistance, and kKmUnits.

Referenced by NuMatrixSpectrum::InverseOscillatePreAveraged(), and NuMatrixSpectrum::OscillatePreAveraged().

00222 {
00223   double arg = 2 * k1267 * dm2 * (kFDDistance / kKmUnits) / energy;
00224   double avg = 0.5 * arg * binwidth / energy;
00225   double sinc = 1;
00226   if(avg!=0) sinc = sin(avg) / avg;
00227 
00228   return 1 - sn2 / 2 * (1 - cos(arg) * sinc );
00229 }

double NuOscProbCalc::DecayWeightCC ( double  energy,
double  alpha,
double  sn2,
double  dm2 = 0,
double  binwidth = 0 
)

The decay weight that should be applied to truly CC events. This takes account of the the additional events that transition to tau.

Definition at line 232 of file NuOscProbCalc.cxx.

References k1267, kFDDistance, kKmUnits, and SQR.

Referenced by NuMatrix1D::DecayCC(), NuDSTAna::Oscillate(), and NuDSTAna::OscillateTrans().

00237 {
00238   double arg1 = - alpha * (kFDDistance/kKmUnits) / (2*energy);
00239   double arg2 = 2 * k1267 * dm2 * (kFDDistance/kKmUnits) / energy;
00240   double avg  = 0.5 * arg2 * binwidth / energy;
00241   double sinc = 1;
00242   if(avg!=0) sinc = sin(avg)/avg;
00243   if(arg1 > 50) arg1 = 50;
00244   const double expPart = exp(arg1);
00245 
00246   return SQR(1-sn2)+SQR(sn2)*SQR(expPart)+2*(1-sn2)*sn2*expPart*cos(arg2)*sinc;
00247 }

double NuOscProbCalc::DecayWeightMuToTau ( double  energy,
double  alpha,
double  sn2,
double  dm2 = 0,
double  binwidth = 0 
)

Probability for muon neutrino to oscillate to tau neutrino, while decay is also occuring.

Definition at line 251 of file NuOscProbCalc.cxx.

References k1267, kFDDistance, kKmUnits, and SQR.

Referenced by NuMatrix1D::DecayMuToTau(), GhostSample::InterpolateOscillatedSpectra(), NuDSTAna::Oscillate(), and NuDSTAna::OscillateTrans().

00256 {
00257   double arg1 = - alpha * (kFDDistance/kKmUnits) / (2*energy);
00258   double arg2 = 2 * k1267 * dm2 * (kFDDistance/kKmUnits) / energy;
00259   double avg  = 0.5 * arg2 * binwidth / energy;
00260   double sinc = 1;
00261   if(avg!=0) sinc = sin(avg)/avg;
00262   if(arg1 > 50) arg1 = 50;
00263   const double expPart = exp(arg1);
00264 
00265   return (1-sn2)*sn2*(1+SQR(expPart)-2*expPart*cos(arg2)*sinc);
00266 }

double NuOscProbCalc::DecayWeightNC ( double  energy,
double  alpha,
double  sn2 
)

The decay weight that should be applied to truly NC events. This includes events that are now nu_tau but still visible in the NC spectrum.

Definition at line 270 of file NuOscProbCalc.cxx.

References kFDDistance, and kKmUnits.

Referenced by ArrayTH1D::Decay(), NuMatrix1D::DecayNC(), NuDSTAna::Oscillate(), and NuDSTAna::OscillateTrans().

00273 {
00274   double arg = -(alpha*kFDDistance/kKmUnits)/energy;
00275   if(arg > 50) arg = 50;
00276   return (1-sn2)+sn2*exp(arg);
00277 }

double NuOscProbCalc::DecoherenceWeight ( double  energy,
double  mu2,
double  sn2,
double  dm2 = 0,
int  power = -1,
double  binwidth = 0 
)

Definition at line 281 of file NuOscProbCalc.cxx.

References k1267, kFDDistance, kKmUnits, and SQR.

Referenced by NuMatrix1D::Decohere(), NuMatrixMethod::DecohereFD(), GhostSample::InterpolateOscillatedSpectra(), NuMatrix1D::InverseDecohere(), NuDSTAna::Oscillate(), NuDSTAna::OscillateTrans(), and NuMatrixMethod::PurityCorrectFDBackgroundDecohereCombined().

00287 {
00288   double arg;
00289   if(power==-1) arg = 2 * k1267 * mu2 * (kFDDistance/kKmUnits) / energy;
00290   else if(power==2) arg = 4 * k1267 * 1e18 * (kFDDistance/kKmUnits) * SQR(energy) / mu2;
00291   else arg = 4 * k1267 * 1e18 * mu2 * (kFDDistance/kKmUnits) * pow(energy,power);
00292 
00293   double arg2 = 2 * k1267 * dm2 * (kFDDistance/kKmUnits) / energy;
00294   double avg = 0.5 * arg2 * binwidth / energy;
00295   if(arg < -50) arg = -50;
00296   double sinc = 1;
00297   if(avg!=0) sinc = sin(avg)/avg;
00298 
00299   return 1 - sn2 * ( 1 - exp(-arg) * cos(arg2) * sinc ) / 2;
00300 }

double NuOscProbCalc::FourFlavourDisappearanceWeight ( const double  energy,
double  dm232,
const double  theta23,
double  dm221,
double  dm243,
const double  theta12,
const double  theta13,
const double  theta14,
const double  theta24,
const double  theta34,
double  delta1,
double  delta2,
double  delta3,
const double  baseline 
)

Definition at line 726 of file NuOscProbCalc.cxx.

References k1267, and kKmUnits.

00733 {
00734 
00735   // Calculate other mass splittings
00736   const double dm231 = dm221 + dm232;
00737   const double dm241 = dm231 + dm243;
00738 
00739   const double c12 = cos(theta12); const double s12 = sin(theta12);
00740   const double c13 = cos(theta13); const double s13 = sin(theta13);
00741   const double c14 = cos(theta14); const double s14 = sin(theta14);
00742   const double c23 = cos(theta23); const double s23 = sin(theta23);
00743   const double c24 = cos(theta24); const double s24 = sin(theta24);
00744 
00745   complex<double> expNegCP13 = complex<double>(cos(delta1), -sin(delta1));
00746   complex<double> expNegCP14 = complex<double>(cos(delta2), -sin(delta2));
00747   complex<double> expNegCP24 = complex<double>(cos(delta3), -sin(delta3));
00748 
00749   complex<double> Umu2  =  c12 * c23 * c24
00750                         -  c24 * s12 * s13 * s23 * conj(expNegCP13)
00751                         -  c13 * s12 * s14 * s24 * expNegCP24 * conj(expNegCP14);
00752 
00753   complex<double> Umu3  =  c13 * c24 * s23
00754                         -  s13 * s14 * s24 * expNegCP13 * expNegCP24 * conj(expNegCP14);
00755   
00756   complex<double> Umu4  =  c14 * s24 * expNegCP24;
00757   
00758   complex<double> i(0.0, 1.0);
00759   
00760   double DeltaM21 = (k1267 * dm221 * baseline / kKmUnits) / energy;
00761   double DeltaM31 = (k1267 * dm231 * baseline / kKmUnits) / energy;
00762   double DeltaM41 = (k1267 * dm241 * baseline / kKmUnits) / energy;
00763   
00764   complex<double> expDeltaM21 = complex<double>(cos(DeltaM21), -sin(DeltaM21));
00765   complex<double> expDeltaM31 = complex<double>(cos(DeltaM31), -sin(DeltaM31));
00766   complex<double> expDeltaM41 = complex<double>(cos(DeltaM41), -sin(DeltaM41));
00767   
00768   double oscProb  = norm(1.0 
00769                          - 2.0 * i * conj(Umu2) * Umu2 * sin(DeltaM21) * expDeltaM21 
00770                          - 2.0 * i * conj(Umu3) * Umu3 * sin(DeltaM31) * expDeltaM31 
00771                          - 2.0 * i * conj(Umu4) * Umu4 * sin(DeltaM41) * expDeltaM41);
00772   
00773   return oscProb;
00774 }

double NuOscProbCalc::FourFlavourNuESurvivalProbability ( const double  energy,
double  dm232,
const double  theta23,
double  dm221,
double  dm243,
const double  theta12,
const double  theta13,
const double  theta14,
const double  theta24,
const double  theta34,
double  delta1,
double  delta2,
double  delta3,
const double  baseline 
)

Definition at line 844 of file NuOscProbCalc.cxx.

References k1267, and kKmUnits.

00851 {
00852 
00853   // Calculate other mass splittings
00854   const double dm231 = dm221 + dm232;
00855   const double dm241 = dm231 + dm243;
00856 
00857   const double s12 = sin(theta12);
00858   const double c13 = cos(theta13); const double s13 = sin(theta13);
00859   const double c14 = cos(theta14); const double s14 = sin(theta14);
00860 
00861   complex<double> expNegCP13 = complex<double>(cos(delta1), -sin(delta1));
00862   complex<double> expNegCP14 = complex<double>(cos(delta2), -sin(delta2));
00863 
00864   complex<double> Ue2   =  c13 * c14 * s12;
00865   complex<double> Ue3   =  c14 * s13 * expNegCP13;
00866   complex<double> Ue4   =  s14 * expNegCP14;
00867   
00868   complex<double> i(0.0, 1.0);
00869   
00870   double DeltaM21 = (k1267 * dm221 * baseline / kKmUnits) / energy;
00871   double DeltaM31 = (k1267 * dm231 * baseline / kKmUnits) / energy;
00872   double DeltaM41 = (k1267 * dm241 * baseline / kKmUnits) / energy;
00873   
00874   complex<double> expDeltaM21 = complex<double>(cos(DeltaM21), -sin(DeltaM21));
00875   complex<double> expDeltaM31 = complex<double>(cos(DeltaM31), -sin(DeltaM31));
00876   complex<double> expDeltaM41 = complex<double>(cos(DeltaM41), -sin(DeltaM41));
00877   
00878   double oscProb = norm(1.0 
00879                         - 2.0 * i * conj(Ue2) * Ue2 * sin(DeltaM21) * expDeltaM21 
00880                         - 2.0 * i * conj(Ue3) * Ue3 * sin(DeltaM31) * expDeltaM31 
00881                         - 2.0 * i * conj(Ue4) * Ue4 * sin(DeltaM41) * expDeltaM41);
00882 
00883   return oscProb;
00884 } // FourFlavourNuESurvivalProbability

double NuOscProbCalc::FourFlavourNuMuToNuEProbability ( const double  energy,
double  dm232,
const double  theta23,
double  dm221,
double  dm243,
const double  theta12,
const double  theta13,
const double  theta14,
const double  theta24,
const double  theta34,
double  delta1,
double  delta2,
double  delta3,
const double  baseline 
)

Definition at line 889 of file NuOscProbCalc.cxx.

References k1267, and kKmUnits.

00896 {
00897 
00898   // Calculate other mass splittings
00899   const double dm231 = dm221 + dm232;
00900   const double dm241 = dm231 + dm243;
00901 
00902   const double c12 = cos(theta12); const double s12 = sin(theta12);
00903   const double c13 = cos(theta13); const double s13 = sin(theta13);
00904   const double c14 = cos(theta14); const double s14 = sin(theta14);
00905   const double c23 = cos(theta23); const double s23 = sin(theta23);
00906   const double c24 = cos(theta24); const double s24 = sin(theta24);
00907 
00908   complex<double> expNegCP13 = complex<double>(cos(delta1), -sin(delta1));
00909   complex<double> expNegCP14 = complex<double>(cos(delta2), -sin(delta2));
00910   complex<double> expNegCP24 = complex<double>(cos(delta3), -sin(delta3));
00911   
00912   complex<double> Umu2  =  c12 * c23 * c24
00913                         -  c24 * s12 * s13 * s23 * conj(expNegCP13)
00914                         -  c13 * s12 * s14 * s24 * expNegCP24 * conj(expNegCP14);
00915 
00916   complex<double> Umu3  =  c13 * c24 * s23
00917                         -  s13 * s14 * s24 * expNegCP13 * expNegCP24 * conj(expNegCP14);
00918   
00919   complex<double> Umu4  =  c14 * s24 * expNegCP24;
00920 
00921 
00922   complex<double> Ue2   =  c13 * c14 * s12;
00923   complex<double> Ue3   =  c14 * s13 * expNegCP13;
00924   complex<double> Ue4   =  s14 * expNegCP14;
00925   
00926   complex<double> i(0.0, 1.0);
00927   
00928   double DeltaM21 = (k1267 * dm221 * baseline / kKmUnits) / energy;
00929   double DeltaM31 = (k1267 * dm231 * baseline / kKmUnits) / energy;
00930   double DeltaM41 = (k1267 * dm241 * baseline / kKmUnits) / energy;
00931   
00932   complex<double> expDeltaM21 = complex<double>(cos(DeltaM21), -sin(DeltaM21));
00933   complex<double> expDeltaM31 = complex<double>(cos(DeltaM31), -sin(DeltaM31));
00934   complex<double> expDeltaM41 = complex<double>(cos(DeltaM41), -sin(DeltaM41));
00935   
00936   double oscProb = norm(-2.0 * i * conj(Umu2) * Ue2 * sin(DeltaM21) * expDeltaM21       
00937                         -2.0 * i * conj(Umu3) * Ue3 * sin(DeltaM31) * expDeltaM31
00938                         -2.0 * i * conj(Umu4) * Ue4 * sin(DeltaM41) * expDeltaM41);
00939 
00940   return oscProb;
00941 } // FourFlavourNuMuToNuEProbability()

double NuOscProbCalc::FourFlavourNuMuToNuSProbability ( const double  energy,
double  dm232,
const double  theta23,
double  dm221,
double  dm243,
const double  theta12,
const double  theta13,
const double  theta14,
const double  theta24,
const double  theta34,
double  delta1,
double  delta2,
double  delta3,
const double  baseline 
)

Definition at line 778 of file NuOscProbCalc.cxx.

References k1267, and kKmUnits.

00785 { 
00786 
00787   // Calculate other mass splittings
00788   const double dm231 = dm221 + dm232;
00789   const double dm241 = dm231 + dm243;
00790 
00791   const double c12 = cos(theta12); const double s12 = sin(theta12);
00792   const double c13 = cos(theta13); const double s13 = sin(theta13);
00793   const double c14 = cos(theta14); const double s14 = sin(theta14);
00794   const double c23 = cos(theta23); const double s23 = sin(theta23);
00795   const double c24 = cos(theta24); const double s24 = sin(theta24);
00796   const double c34 = cos(theta34); const double s34 = sin(theta34);
00797 
00798   complex<double> expNegCP13 = complex<double>(cos(delta1), -sin(delta1));
00799   complex<double> expNegCP14 = complex<double>(cos(delta2), -sin(delta2));
00800   complex<double> expNegCP24 = complex<double>(cos(delta3), -sin(delta3));
00801 
00802   complex<double> Umu2  =  c12 * c23 * c24
00803                         -  c24 * s12 * s13 * s23 * conj(expNegCP13)
00804                         -  c13 * s12 * s14 * s24 * expNegCP24 * conj(expNegCP14);
00805 
00806   complex<double> Umu3  =  c13 * c24 * s23
00807                         -  s13 * s14 * s24 * expNegCP13 * expNegCP24 * conj(expNegCP14);
00808   
00809   complex<double> Umu4  =  c14 * s24 * expNegCP24;
00810   
00811 
00812   complex<double> Us2   =  -c13 * c24 * c34 * s12 * s14 * conj(expNegCP14)
00813                            -c12 * c23 * c34 * s24 * conj(expNegCP24)
00814                            +c34 * s12 * s13 * s23 * s24 * conj(expNegCP13 * expNegCP24)
00815                            +c23 * s12 * s13 * s34 * conj(expNegCP13)
00816                            +c12 * s23 * s34;
00817   
00818   complex<double> Us3   =  -c24 * c34 * s13 * s14 * expNegCP13 * conj(expNegCP14)
00819                            -c13 * c34 * s23 * s24 * conj(expNegCP24)
00820                            -c13 * c23 * s34;
00821   
00822   complex<double> Us4   =  c14 * c24 * c34;
00823   
00824   
00825   complex<double> i(0.0, 1.0);
00826   
00827   double DeltaM21 = (k1267 * dm221 * baseline / kKmUnits) / energy;
00828   double DeltaM31 = (k1267 * dm231 * baseline / kKmUnits) / energy;
00829   double DeltaM41 = (k1267 * dm241 * baseline / kKmUnits) / energy;
00830   
00831   complex<double> expDeltaM21 = complex<double>(cos(DeltaM21), -sin(DeltaM21));
00832   complex<double> expDeltaM31 = complex<double>(cos(DeltaM31), -sin(DeltaM31));
00833   complex<double> expDeltaM41 = complex<double>(cos(DeltaM41), -sin(DeltaM41));
00834   
00835   double oscProb = norm(-2.0 * i * conj(Umu2) * Us2 * sin(DeltaM21) * expDeltaM21
00836                         -2.0 * i * conj(Umu3) * Us3 * sin(DeltaM31) * expDeltaM31
00837                         -2.0 * i * conj(Umu4) * Us4 * sin(DeltaM41) * expDeltaM41);
00838 
00839   return oscProb;
00840 } // FourFlavorNuMuToNuSProbability()

double NuOscProbCalc::FourFlavourNuMuToNuTauProbability ( const double  energy,
double  dm232,
const double  theta23,
double  dm221,
double  dm243,
const double  theta12,
const double  theta13,
const double  theta14,
const double  theta24,
const double  theta34,
double  delta1,
double  delta2,
double  delta3,
const double  baseline 
)

Definition at line 945 of file NuOscProbCalc.cxx.

References k1267, and kKmUnits.

00952 {
00953 
00954   // Calculate other mass splittings
00955   const double dm231 = dm221 + dm232;
00956   const double dm241 = dm231 + dm243;
00957 
00958   const double c12 = cos(theta12); const double s12 = sin(theta12);
00959   const double c13 = cos(theta13); const double s13 = sin(theta13);
00960   const double c14 = cos(theta14); const double s14 = sin(theta14);
00961   const double c23 = cos(theta23); const double s23 = sin(theta23);
00962   const double c24 = cos(theta24); const double s24 = sin(theta24);
00963   const double c34 = cos(theta34); const double s34 = sin(theta34);
00964 
00965   complex<double> expNegCP13 = complex<double>(cos(delta1), -sin(delta1));
00966   complex<double> expNegCP14 = complex<double>(cos(delta2), -sin(delta2));
00967   complex<double> expNegCP24 = complex<double>(cos(delta3), -sin(delta3));
00968 
00969   complex<double> Umu2  =  c12 * c23 * c24
00970                         -  c24 * s12 * s13 * s23 * conj(expNegCP13)
00971                         -  c13 * s12 * s14 * s24 * expNegCP24 * conj(expNegCP14);
00972 
00973   complex<double> Umu3  =  c13 * c24 * s23
00974                         -  s13 * s14 * s24 * expNegCP13 * expNegCP24 * conj(expNegCP14);
00975   
00976   complex<double> Umu4  =  c14 * s24 * expNegCP24;
00977 
00978 
00979   complex<double> Utau2 =  -c12 * c34 * s23
00980                            -c23 * c34 * s12 * s13 * conj(expNegCP13)
00981                            -c13 * c24 * s12 * s14 * s34 * conj(expNegCP14)
00982                            -c12 * c23 * s24 * s34 * conj(expNegCP24)
00983                            +s12 * s13 * s23 * s24 * s34 * conj(expNegCP13 * expNegCP24);
00984 
00985   complex<double> Utau3 =  c13 * c23 * c34
00986                         -  c24 * s13 * s14 * s34 * expNegCP13 * conj(expNegCP14)
00987                         -  c13 * s23 * s24 * s34 * conj(expNegCP24);
00988 
00989   complex<double> Utau4 =  c14 * c24 * s34;
00990 
00991   complex<double> i(0.0, 1.0);
00992   
00993   double DeltaM21 = (k1267 * dm221 * baseline / kKmUnits) / energy;
00994   double DeltaM31 = (k1267 * dm231 * baseline / kKmUnits) / energy;
00995   double DeltaM41 = (k1267 * dm241 * baseline / kKmUnits) / energy;
00996   
00997   complex<double> expDeltaM21 = complex<double>(cos(DeltaM21), -sin(DeltaM21));
00998   complex<double> expDeltaM31 = complex<double>(cos(DeltaM31), -sin(DeltaM31));
00999   complex<double> expDeltaM41 = complex<double>(cos(DeltaM41), -sin(DeltaM41));
01000   
01001   double oscProb =  norm(-2.0 * i * conj(Umu2) * Utau2 * sin(DeltaM21) * expDeltaM21     
01002                          -2.0 * i * conj(Umu3) * Utau3 * sin(DeltaM31) * expDeltaM31
01003                          -2.0 * i * conj(Umu4) * Utau4 * sin(DeltaM41) * expDeltaM41);
01004 
01005   return oscProb;
01006 } // FourFlavourNuMuToNuTauProbability()

double NuOscProbCalc::HighMassSterileWeight ( double  energy,
double  dm2,
double  sn2,
double  sn224 
)

Definition at line 378 of file NuOscProbCalc.cxx.

References MuELoss::e, k1267, kFDDistance, kKmUnits, and SQR.

00380 {
00381   if (energy < 1e-6){return 0.5;}
00382 
00383   //Calculate theta23:
00384   Double_t theta23 = TMath::ASin(sqrt(sn2)) / 2.0;
00385   //Calculate theta24:
00386   Double_t theta24 = TMath::ASin(sqrt(sn224)) / 2.0;
00387 
00388   Double_t s23 = TMath::Sin(theta23);
00389   Double_t c24 = TMath::Cos(theta24);
00390   
00391   Double_t oscProb = 1.0 -
00392     (4.0 * s23*s23 * c24*c24*c24*c24 * (1.0 - s23*s23) *
00393      SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)))
00394     - 0.5 * sn224;
00395 
00396   return oscProb;
00397 }

double NuOscProbCalc::LargeExtraDimensionDisappearanceWeight ( double  th12,
double  th13,
double  th23,
double  deltacp,
double  dm221,
double  dm232,
double  m0,
double  a,
double  energy,
double  baseLine 
)

Definition at line 1702 of file NuOscProbCalc.cxx.

References LargeExtraDimensionProbability().

Referenced by FarOverNearFitLED::MCOscillated().

01706                                                                                             {
01707 
01708   return NuOscProbCalc::LargeExtraDimensionProbability(1, 1, 
01709                                               th12, th13, th23,
01710                                               deltacp,
01711                                               dm221, dm232,
01712                                               m0, a,
01713                                               energy, baseLine);
01714 }

double NuOscProbCalc::LargeExtraDimensionNuESurvivalProbability ( double  th12,
double  th13,
double  th23,
double  deltacp,
double  dm221,
double  dm232,
double  m0,
double  a,
double  energy,
double  baseLine 
)

Definition at line 1717 of file NuOscProbCalc.cxx.

References LargeExtraDimensionProbability().

Referenced by FarOverNearFitLED::MCOscillated().

01721                                                                                                {
01722 
01723   return NuOscProbCalc::LargeExtraDimensionProbability(0, 0, 
01724                                               th12, th13, th23,
01725                                               deltacp,
01726                                               dm221, dm232,
01727                                               m0, a,
01728                                               energy, baseLine);
01729 }

double NuOscProbCalc::LargeExtraDimensionNuMuToNuEProbability ( double  th12,
double  th13,
double  th23,
double  deltacp,
double  dm221,
double  dm232,
double  m0,
double  a,
double  energy,
double  baseLine 
)

Definition at line 1732 of file NuOscProbCalc.cxx.

References LargeExtraDimensionProbability().

Referenced by FarOverNearFitLED::MCOscillated().

01736                                                                                              {
01737 
01738   return NuOscProbCalc::LargeExtraDimensionProbability(1, 0, 
01739                                               th12, th13, th23,
01740                                               deltacp,
01741                                               dm221, dm232,
01742                                               m0, a,
01743                                               energy, baseLine);
01744 }

double NuOscProbCalc::LargeExtraDimensionNuMuToNuSProbability ( double  th12,
double  th13,
double  th23,
double  deltacp,
double  dm221,
double  dm232,
double  m0,
double  a,
double  energy,
double  baseLine 
)

Definition at line 1673 of file NuOscProbCalc.cxx.

References LargeExtraDimensionProbability().

Referenced by FarOverNearFitLED::MCOscillated().

01677                                                                                              {
01678 
01679   return 1. 
01680     - NuOscProbCalc::LargeExtraDimensionProbability(1, 0, 
01681                                            th12, th13, th23,
01682                                            deltacp,
01683                                            dm221, dm232,
01684                                            m0, a,
01685                                            energy, baseLine)
01686     - NuOscProbCalc::LargeExtraDimensionProbability(1, 1, 
01687                                            th12, th13, th23,
01688                                            deltacp,
01689                                            dm221, dm232,
01690                                            m0, a,
01691                                            energy, baseLine)
01692     - NuOscProbCalc::LargeExtraDimensionProbability(1, 2, 
01693                                            th12, th13, th23,
01694                                            deltacp,
01695                                            dm221, dm232,
01696                                            m0, a,
01697                                            energy, baseLine);
01698 
01699 }

double NuOscProbCalc::LargeExtraDimensionNuMuToNuTauProbability ( double  th12,
double  th13,
double  th23,
double  deltacp,
double  dm221,
double  dm232,
double  m0,
double  a,
double  energy,
double  baseLine 
)

Definition at line 1747 of file NuOscProbCalc.cxx.

References LargeExtraDimensionProbability().

Referenced by FarOverNearFitLED::MCOscillated().

01751                                                                                                {
01752 
01753   return NuOscProbCalc::LargeExtraDimensionProbability(1, 2, 
01754                                               th12, th13, th23,
01755                                               deltacp,
01756                                               dm221, dm232,
01757                                               m0, a,
01758                                               energy, baseLine);
01759 }

double NuOscProbCalc::LargeExtraDimensionProbability ( int  initialFlavor,
int  finalFlavor,
double  th12,
double  th13,
double  th23,
double  deltacp,
double  dm221,
double  dm232,
double  m0,
double  a,
double  energy,
double  baseLine 
)

Definition at line 1217 of file NuOscProbCalc.cxx.

References H_GSL, kKmUnits, Munits::m, and n.

Referenced by LargeExtraDimensionDisappearanceWeight(), LargeExtraDimensionNuESurvivalProbability(), LargeExtraDimensionNuMuToNuEProbability(), LargeExtraDimensionNuMuToNuSProbability(), LargeExtraDimensionNuMuToNuTauProbability(), and NuMatrixSpectrum::OscillateLED().

01222                                                                                      {
01223   
01224   /******************************************************************************
01225    * Calculate oscillation probabilities with large extra dimension (LED)       *
01226    ******************************************************************************
01227    * Parameters:                                                                *
01228    * 1   initialFlavor: initial neutrino flavor                                 *
01229    *                    0 for e, 1 for mu and 2 for tau                         *
01230    * 2   finalFlavor:   final neutrino flavor                                   *
01231    * 3   th12;          theta12                                                 *
01232    * 4   th13;          theta13                                                 *
01233    * 5   th23;          theta23                                                 *
01234    * 6   deltacp;       radian, CP phase                                        *
01235    * 7   dm221;         eV2, m_2^2 - m_1^2                                      *
01236    * 8   dm232;         eV2, m_3^2 - m_2^2                                      *
01237    * 9   m0;            eV, smallest mass                                       *
01238    * 10  a;             m, size of the large extra dimension (LED)              *
01239    * 11  energy:        GeV, neutrino energy                                    *
01240    * 12  baseLine:      km, baseline                                            *
01241    ******************************************************************************/    
01242 
01243   //---------------- parameters
01244 
01245   double Prob;
01246   const int NUM_NU_FLAVOURS = 3;            // number of neutrino flavours
01247   bool meff                 = false;        // matter effect
01248   const int modeN           = 5 ;           // max KK mode number used in the computation   
01249 
01250   /*
01251   printf("\n");
01252   printf("th12    = %.6f\n", th12);
01253   printf("th13    = %.6f\n", th13);
01254   printf("th23    = %.6f\n", th23);
01255   printf("deltacp = %.6f\n", deltacp);
01256   printf("dm221   = %.6f\n", dm221);
01257   printf("dm232   = %.6f\n", dm232);
01258   printf("m0      = %.6f\n", m0);
01259   printf("a       = %.6f\n", a);
01260   */
01261 
01262   //---------------- constants
01263 
01264   baseLine = baseLine/kKmUnits;                                 // km, baseline
01265   double rho = 2.8;                                             // g cm-3, earth density
01266   double GF = 1.16637E-5;                                       // GeV-2, Fermi constant
01267   double NA = GSL_CONST_NUM_AVOGADRO;                           // Avogadro number
01268   double GeVtoinvcm = 5.07E13;                                  // 1 GeV = 5.07E13 cm-1
01269   double invcmtoGeV = 1. / GeVtoinvcm;
01270   double cmtoinvGeV = GeVtoinvcm;                               // 1 cm = 5.07E13 GeV-1
01271   double kmtoinvGeV = 1.E5 * GeVtoinvcm;
01272   double eVtoGeV = 1.E-9;                                       // 1 eV = 1E-9 GeV
01273   double mtocm = 100.;                                          // 1 m = 100 cm
01274   double ne = 1.4 * NA;                                         // cm-3, electron number densities
01275   double u = GSL_CONST_MKSA_UNIFIED_ATOMIC_MASS;                // kg, atomic mass unit
01276   double nn = rho / 1000. / u / 2.;                             // cm-3, neutron number density
01277   double VCC =  sqrt(2.) * GF * ne * pow(invcmtoGeV, 3);        // GeV
01278   double VNC = -sqrt(2.) / 2. * GF * nn * pow(invcmtoGeV, 3);   // GeV
01279   const int H_size = NUM_NU_FLAVOURS * (modeN + 1);             // size of Hamiltonian
01280   complex <double> I = sqrt( complex <double> (-1) );           // imaginary unit i
01281   baseLine *= kmtoinvGeV;                                          // GeV-1, Baseline
01282   if(!meff){
01283     VCC = 0.;
01284     VNC = 0.;
01285   }
01286 
01287   /*
01288   printf("\n");
01289   printf("rho       = %.6f\n", rho);
01290   printf("baseLine  = %g\n",   baseLine / GLB_KM_TO_EV(1.) / 1.E9);
01291   printf("ne        = %g\n",   ne);
01292   printf("nn        = %g\n",   nn);
01293   printf("VCC       = %g\n",   VCC);
01294   printf("VNC       = %g\n",   VNC);
01295   printf("VCC + VNC = %g\n", VCC + VNC);
01296   printf("\n");
01297   */
01298 
01299   //---------------- vacuum mixing matrix
01300 
01301   gsl_matrix_complex* U_GSL = NULL; 
01302   U_GSL = gsl_matrix_complex_calloc(NUM_NU_FLAVOURS, NUM_NU_FLAVOURS);
01303   complex <double> (*U_C)[NUM_NU_FLAVOURS] 
01304     = (complex <double> (*)[NUM_NU_FLAVOURS]) gsl_matrix_complex_ptr(U_GSL, 0, 0);
01305   U_C[0][0] = cos(th12)*cos(th13);
01306   U_C[0][1] = sin(th12)*cos(th13);
01307   U_C[0][2] = sin(th13) * exp(-I * deltacp);
01308 
01309   U_C[1][0] = -sin(th12)*cos(th23) - cos(th12)*sin(th23)*sin(th13) * exp(I*deltacp);
01310   U_C[1][1] =  cos(th12)*cos(th23) - sin(th12)*sin(th23)*sin(th13) * exp(I*deltacp);
01311   U_C[1][2] =  sin(th23)*cos(th13);
01312 
01313   U_C[2][0] =  sin(th12)*sin(th23) - cos(th12)*cos(th23)*sin(th13) * exp(I*deltacp);
01314   U_C[2][1] = -cos(th12)*sin(th23) - sin(th12)*cos(th23)*sin(th13) * exp(I*deltacp);
01315   U_C[2][2] =  cos(th23)*cos(th13);
01316 
01317   /*
01318   {
01319     int i, j;
01320     printf("\n");
01321     printf("Mixing Matrix U \n");
01322     for(i = 0; i < NUM_NU_FLAVOURS; i++){
01323       for(j = 0; j < NUM_NU_FLAVOURS; j++){
01324         printf("%f + %fi        ", 
01325                GSL_REAL(gsl_matrix_complex_get(U_GSL, i, j)), 
01326                GSL_IMAG(gsl_matrix_complex_get(U_GSL, i, j)));
01327       }
01328       printf("\n");
01329     }
01330   }
01331   */
01332 
01333   //---------------- vector m (eV, absolute nu mass)
01334 
01335   double m[NUM_NU_FLAVOURS];                                    // eV, mass eigenstate, m[0] = m1, m[1] = m2, m[2] = m3
01336   if(dm232 > 0){                                                // normal hierarchy
01337     m[0] = m0;
01338     m[1] = sqrt(m0 * m0 + dm221);
01339     m[2] = sqrt(dm232 + dm221 + m0 * m0);
01340   }
01341   else{                                                         // inverted hierarchy
01342     m[0] = sqrt(m0 * m0 - dm232 - dm221);
01343     m[1] = sqrt(m0 * m0 - dm232);
01344     m[2] = m0;
01345   }
01346 
01347   /*
01348   {
01349     int i;
01350     printf("\nVector m (eV) \n");
01351     for(i = 0; i < NUM_NU_FLAVOURS; i++)
01352       printf("m[%d] = %g\n", i, m[i]);
01353   }
01354   */
01355 
01356   //---------------- vector xi
01357 
01358   double xi[NUM_NU_FLAVOURS];
01359   {
01360     int i;
01361     for(i = 0; i < NUM_NU_FLAVOURS; i++)
01362       xi[i] = sqrt(2.) * (m[i] * eVtoGeV) * (a * mtocm * cmtoinvGeV);
01363   }
01364 
01365   /*
01366   {
01367     int i;
01368     printf("\nVector xi \n");
01369     for(i = 0; i < NUM_NU_FLAVOURS; i++)
01370       printf("xi[%d] = %g\n", i, xi[i]);
01371   }
01372   */
01373 
01374   //---------------- eta vector
01375 
01376   double eta[NUM_NU_FLAVOURS];
01377   {
01378     int i;
01379     for(i = 0; i < NUM_NU_FLAVOURS; i++ )
01380       eta[i] = (modeN + 0.5) * pow(xi[i], 2);
01381   }
01382 
01383   //---------------- vector V_alpha (GeV)
01384 
01385   double Va[NUM_NU_FLAVOURS];
01386   Va[0] = VCC + VNC;
01387   Va[1] = VNC;
01388   Va[2] = VNC;
01389 
01390   /*
01391   {
01392     int i;
01393     printf("\nVector Va \n");
01394     for(i = 0; i < NUM_NU_FLAVOURS; i++)
01395       printf("Va[%d] = %g\n", i, Va[i]);
01396   }
01397   */
01398 
01399   //---------------- matrix V
01400 
01401   gsl_matrix_complex* V_GSL = NULL; 
01402   V_GSL = gsl_matrix_complex_calloc(NUM_NU_FLAVOURS, NUM_NU_FLAVOURS);
01403   complex <double> (*V_C)[NUM_NU_FLAVOURS] 
01404     = (complex <double> (*)[NUM_NU_FLAVOURS]) gsl_matrix_complex_ptr(V_GSL, 0, 0);
01405   {
01406     int i, j, k;
01407     for(i = 0; i < NUM_NU_FLAVOURS; i++)
01408       for(j = 0; j < NUM_NU_FLAVOURS; j++){
01409         
01410         complex <double> sum = 0.;
01411         for(k = 0; k < NUM_NU_FLAVOURS; k++)
01412           sum += conj(U_C[k][i]) * U_C[k][j] * Va[k];
01413         
01414         V_C[i][j] = 2. * energy * sum * pow(a * mtocm * cmtoinvGeV, 2);
01415       }
01416   }
01417   
01418   /*
01419   {
01420     int i, j;
01421     printf("\n");
01422     printf("Matrix V \n");
01423     for(i = 0; i < NUM_NU_FLAVOURS; i++){
01424       for(j = 0; j < NUM_NU_FLAVOURS; j++){
01425         printf("%.8f + %.8fi        ", 
01426                GSL_REAL(gsl_matrix_complex_get(V_GSL, i, j)), 
01427                GSL_IMAG(gsl_matrix_complex_get(V_GSL, i, j)));
01428       }
01429       printf("\n");
01430     }
01431   }
01432   */
01433 
01434   //---------------- matrix H
01435 
01436   gsl_matrix_complex* H_GSL = NULL;
01437   H_GSL = gsl_matrix_complex_calloc(H_size, H_size);
01438   complex <double> (*H_C)[H_size] 
01439     = (complex <double> (*)[H_size]) gsl_matrix_complex_ptr(H_GSL, 0, 0);
01440   gsl_matrix_complex_set_zero(H_GSL);
01441 
01442   // first 3 by 3 block
01443   {
01444     int i, j;
01445     for(i = 0; i < 3; i++ )
01446       for(j = 0; j < 3; j++ ){
01447       if(i != j)
01448         H_C[i][j] = V_C[i][j];
01449       else 
01450         H_C[i][j] = V_C[i][j] + eta[i];
01451       }
01452   }
01453 
01454   // diagonal blocks
01455   {
01456     int i, j;
01457     for(j = 1; j < modeN + 1; j++)
01458       for(i = 3 * j; i < 3 * j + 3; i++)
01459         H_C[i][i] = j * j;
01460   }
01461 
01462   // first column blocks
01463   {
01464     int i, j;
01465     for(j = 1; j < modeN + 1; j++)
01466       for(i = 3 * j; i < 3 * j + 3; i++)
01467         H_C[i][i - 3 * j] = (double)j * xi[i - 3 * j];
01468   }
01469 
01470   // first row blocks
01471   {
01472     int i, j;
01473     for(j = 1; j < modeN + 1; j++)
01474       for(i = 3 * j; i < 3 * j + 3; i++)
01475         H_C[i - 3 * j][i] = (double)j * xi[i - 3 * j];
01476   }
01477 
01478   /*
01479   {
01480     int i, j;
01481     int row_begin = 1;
01482     int row_end = 3;
01483     int col_begin = 1;
01484     int col_end = 3;
01485     printf("\n");
01486     printf("Matrix H \n");
01487     for(i = row_begin - 1; i < row_end; i++){
01488       for(j = col_begin - 1; j < col_end; j++){
01489         printf("%.8f + %.8fi        ", 
01490                GSL_REAL(gsl_matrix_complex_get(H_GSL, i, j)), 
01491                GSL_IMAG(gsl_matrix_complex_get(H_GSL, i, j)));
01492       }
01493       printf("\n");
01494     }
01495   }
01496   */
01497 
01498   //---------------- eigen values and eigen states of matrix H
01499 
01500   gsl_vector* eval = gsl_vector_alloc (H_size); 
01501   gsl_matrix_complex* evec = gsl_matrix_complex_alloc (H_size, H_size);
01502   gsl_eigen_hermv_workspace* w = gsl_eigen_hermv_alloc (H_size);
01503 
01504   gsl_eigen_hermv (H_GSL, eval, evec, w); 
01505   gsl_eigen_hermv_free (w); 
01506   gsl_eigen_hermv_sort (eval, evec, GSL_EIGEN_SORT_ABS_ASC);
01507 
01508   /*
01509   {
01510     int i, j;
01511     for (i = 0; i < H_size; i++)
01512       {
01513         printf("\n-------------");
01514         double eval_i = gsl_vector_get (eval, i);
01515         gsl_vector_complex_view evec_i = gsl_matrix_complex_column (evec, i);
01516         printf ("eigenvalue = %g \n\n", eval_i);
01517         printf ("eigenvector \n"); 
01518         for (j = 0; j < H_size; ++j){
01519           gsl_complex z = gsl_vector_complex_get(&evec_i.vector, j); 
01520           printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));
01521         } 
01522       }
01523   }
01524   */
01525 
01526   /*
01527   {
01528     int i;
01529     int j = 0;
01530     int k = 1;
01531     gsl_complex dotprod = gsl_complex_rect(0., 0.);
01532     gsl_vector_complex_view evec_1 = gsl_matrix_complex_column (evec, j);
01533     gsl_vector_complex_view evec_2 = gsl_matrix_complex_column (evec, k);
01534     for (i = 0; i < H_size; i++){
01535         gsl_complex z1 = gsl_vector_complex_get(&evec_1.vector, i);
01536         gsl_complex z2 = gsl_vector_complex_get(&evec_2.vector, i);
01537         dotprod = gsl_complex_add(dotprod, gsl_complex_mul(gsl_complex_conjugate(z1), z2));
01538     }
01539     printf("\nevec_%d dot evec_%d = %g + %gi\n", j, k, GSL_REAL(dotprod), GSL_IMAG(dotprod));
01540   }
01541   */
01542 
01543   //---------------- lambda[j][N]
01544 
01545   double lambda[modeN + 1][NUM_NU_FLAVOURS];
01546 
01547   {
01548     int i;
01549     for(i = 0; i < H_size; i++)
01550       if(dm232 > 0)
01551         lambda[ (int)i/NUM_NU_FLAVOURS ][i%NUM_NU_FLAVOURS] = gsl_vector_get (eval, i);
01552       else
01553         lambda[ (int)i/NUM_NU_FLAVOURS ][(i+2)%NUM_NU_FLAVOURS] = gsl_vector_get (eval, i);
01554   }
01555 
01556   /*
01557   {
01558     printf("\nmatrix of lambda[j][N] \n");
01559     int i, j;
01560     for(i = 0; i < modeN + 1; i++){
01561       for(j = 0; j < NUM_NU_FLAVOURS; j++)
01562         printf("%g    ", lambda[i][j]);
01563       printf("\n");
01564     }  
01565   }  
01566   */
01567 
01568   //---------------- cross-check gsl results by solving lambda via pi xi2 / 2 cot (pi lambda) = lambda
01569 
01570   /*
01571   {
01572     int i = 1;
01573     double preFcn;
01574     double fcn;
01575     double epsilon = 1.E-8;
01576     double stepSize = 1.E-4;
01577     double initVal = 0.;
01578     double maxVal = 0.05;
01579     double x;
01580 
01581     x = initVal;
01582     while(x < maxVal){
01583       x += stepSize;
01584       fcn = Pi * pow(xi[i], 2) / 2. * tan(Pi / 2. - Pi * x) - x;
01585       if(preFcn > 0 && fcn < 0){
01586         x -= stepSize;
01587         stepSize /= 2.;
01588       }
01589       else
01590         preFcn = fcn;
01591 
01592       if(abs(fcn) < epsilon){
01593         cout << "pow(x, 2) = " << pow(x, 2) << "  fcn = " << fcn << endl;
01594         break;
01595       }
01596 
01597     }
01598   }
01599   */
01600 
01601   //---------------- vector W_C[i][N]
01602 
01603   complex <double> W_C[modeN + 1][modeN + 1][NUM_NU_FLAVOURS][NUM_NU_FLAVOURS];
01604 
01605   {
01606     int m;
01607     for(m = 0; m < H_size; m++){
01608       gsl_vector_complex_view evec_m = gsl_matrix_complex_column (evec, m);
01609       int n;
01610       for(n = 0; n < H_size; n++){
01611         gsl_complex z = gsl_vector_complex_get(&evec_m.vector, n);
01612         int i, j;
01613         if(dm232 > 0)
01614           i = m % NUM_NU_FLAVOURS;
01615         else
01616           i = (m + 2) % NUM_NU_FLAVOURS;
01617 
01618         j = n % NUM_NU_FLAVOURS;
01619         W_C[(int)m / NUM_NU_FLAVOURS][(int)n / NUM_NU_FLAVOURS][i][j] = GSL_REAL(z) + GSL_IMAG(z) * I;
01620       }
01621     }
01622   }
01623 
01624 
01625   //---------------- transition amplitude A
01626 
01627   int alpha = initialFlavor;
01628   int beta = finalFlavor;
01629   complex <double> A = 0;
01630   {
01631     int i, j, k, n;
01632     for(i = 0; i < NUM_NU_FLAVOURS; i++)
01633       for(j = 0; j < NUM_NU_FLAVOURS; j++)
01634         for(k = 0; k < NUM_NU_FLAVOURS; k++)
01635           for(n = 0; n < modeN + 1; n++){
01636             complex <double> prod = 1.;
01637             prod *= U_C[alpha][i] * conj(U_C[beta][k]);
01638             prod *= conj(W_C[0][n][i][j]) * W_C[0][n][k][j];
01639             prod *= exp(I * lambda[n][j] * baseLine / (2. * energy * pow(a * mtocm * cmtoinvGeV, 2)));
01640             A += prod;
01641           }
01642   }
01643 
01644   Prob = real(A * conj(A));
01645 
01646   //---------------- destroy matrices
01647 
01648   gsl_vector_free(eval); 
01649 
01650   gsl_matrix_complex_free(evec);
01651 
01652   if (U_GSL!=NULL){ 
01653     gsl_matrix_complex_free(U_GSL);
01654     U_GSL = NULL; 
01655   }
01656 
01657   if (V_GSL!=NULL){ 
01658     gsl_matrix_complex_free(V_GSL);
01659     V_GSL = NULL; 
01660   }
01661 
01662   if (H_GSL!=NULL){ 
01663     gsl_matrix_complex_free(H_GSL);
01664     H_GSL = NULL; 
01665   }
01666 
01667   //---------------- return
01668 
01669   return(Prob);
01670 }

double NuOscProbCalc::NSIWeight ( double  energy,
double  dm2,
double  sn2,
double  epsilon,
double  sign 
)

Definition at line 303 of file NuOscProbCalc.cxx.

References k1267, kFDDistance, kKmUnits, and SQR.

Referenced by NuMatrix1D::InverseOscillateNSI(), NuMatrixSpectrum::InverseOscillateNSI(), NuDSTAna::Oscillate(), NuMatrix1D::OscillateNSI(), NuMatrixSpectrum::OscillateNSI(), and NuDSTAna::OscillateTrans().

00308 {
00309   //sign is positive for \pm where plus is for neutrinos
00310 
00311 //   double V = 1.1e-13*Munits::electronvolt;
00312 //   double VE = 1.1e-7*energy; //this is in GeV
00313   double VL = (kFDDistance/kKmUnits)/1900; // V = 1/1900km
00314 
00315   // There are divisions so it fails in the no oscillations case,
00316   // return survival probability =1, don't do epsilon
00317   if((sn2 > 0 && dm2 > 0)){
00318 
00319     double sn = sqrt(sn2);
00320     double B = k1267*dm2*(kFDDistance/kKmUnits)/energy;
00321 
00322     double term1 = sqrt( B*B + sign*2*sn*B*epsilon*VL + SQR(epsilon*VL) );
00323     double term2 = 1 + sign*2*sn*epsilon*VL/B + SQR(epsilon*VL/B);
00324 
00325     double survProb = SQR(TMath::Cos(term1))
00326       + (SQR(TMath::Sin(term1))*(1-sn2))/term2;
00327 
00328 
00329 //     std::cout << "sn2: " << sn2 << ", dm2: " << dm2
00330 //            << ", epsilon: " << epsilon << ", sign: " << sign 
00331 //            << ", survProb= " << survProb  
00332 //            << ", energy= " << energy 
00333 //            << "\nfirst term= " << term1
00334 //            << "\nsecond term= " << term2 
00335 //            << endl;
00336 
00337     return survProb;
00338     
00339   }
00340   else return 1.0;
00341 }

double NuOscProbCalc::NSIWeightPrime ( double  energy,
double  dm2,
double  sn2,
double  epsilon,
double  epsprime,
double  sign 
)

Definition at line 343 of file NuOscProbCalc.cxx.

References k1267, kFDDistance, kKmUnits, and SQR.

00349 {
00350   //sign is positive for \pm where plus is for neutrinos
00351 
00352 //   double V = 1.1e-13*Munits::electronvolt;
00353 //   double VE = 1.1e-7*energy; //this is in GeV
00354   double VL = (kFDDistance/kKmUnits)/1900; // V = 1/1900km
00355 
00356   // There are divisions so it fails in the no oscillations case,
00357   // return survival probability =1, don't do epsilon
00358   if((sn2 > 0 && dm2 > 0)){
00359 
00360     double sn = sqrt(sn2);
00361     double cs = sqrt(1-sn2);
00362     double B = k1267*dm2*(kFDDistance/kKmUnits)/energy;
00363 
00364     double num = B*B*sn2 + SQR(epsilon*VL) + sign*2*B*sn*epsilon*VL;
00365 
00366     double denom = sqrt( B*B + sign*2*B*VL*(sn*epsilon + cs*epsprime) 
00367                          + SQR(epsilon*VL) + SQR(epsprime*VL) );
00368 
00369     double survProb = num/denom;
00370 
00371     return survProb; 
00372   }
00373   else return 1.0;
00374 
00375 }

double NuOscProbCalc::OscillationWeight ( double  energy,
double  dm2,
double  sn2 
)
double NuOscProbCalc::OscillationWeightMu2Mu ( double  energy,
double  dm2,
double  sn2,
double  g_mu,
double  g_tau,
double  c_mumu,
double  c_tautau 
)

Definition at line 51 of file NuOscProbCalc.cxx.

References c_factor, g_factor, k1267, kFDDistance, kKmUnits, and SQR.

Referenced by NuMatrixSpectrum::OscillateMu2Mu(), NuTransSME::OscillationWeight(), and NuDSTAna::OscillationWeightTrans().

00059 {
00060 
00061   //double Theta=0.5*asin(sqrt(sn2));
00062 
00063   double sin2theta = sqrt(sn2);
00064   double cos2theta = cos(asin(sqrt(sn2)));
00065 
00066 
00067   Double_t Theta_Odd = 0.5*(atan2(dm2*sin2theta,((g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00068   Double_t Theta_Evn = 0.5*(atan2(dm2*sin2theta,((-g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00069 
00070   Double_t sin2ThOdd = sin(atan2(dm2*sin2theta,((g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00071   Double_t sin2ThEvn = sin(atan2(dm2*sin2theta,((-g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00072 
00073   double sinSq2ThOdd = sin2ThOdd*sin2ThOdd;
00074   double sinSq2ThEvn = sin2ThEvn*sin2ThEvn;
00075   //double cosSq2ThOdd = 1-sinSq2ThOdd;
00076   //double cosSq2ThEvn = 1-sinSq2ThEvn;
00077 
00078 
00079   Double_t sinThOdd = sin(Theta_Odd);
00080   Double_t sinThEvn = sin(Theta_Evn);
00081 
00082   Double_t sinSqThOdd = sinThOdd*sinThOdd;
00083   Double_t cosSqThOdd = 1-sinSqThOdd;
00084   Double_t sinSqThEvn = sinThEvn*sinThEvn;
00085   Double_t cosSqThEvn = 1-sinSqThEvn;
00086 
00087   //cout<<"energy:"<<energy<<"dm2="<<dm2<<", sn2="<<sn2<<", prob numu2numu="<<0.25*(1-sinSq2ThOdd*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy))+1-sinSq2ThEvn*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)))+0.5*(cosSqThOdd*cosSqThEvn+sinSqThOdd*sinSqThEvn+cosSqThOdd*sinSqThEvn*cos(2*k1267*dm2*(kFDDistance/kKmUnits)/energy)+sinSqThOdd*cosSqThEvn*cos(2*k1267*dm2*(kFDDistance/kKmUnits)/energy))<<endl;
00088 
00089 //nubar to nubar
00090   return 0.25*(1-sinSq2ThOdd*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy))+1-sinSq2ThEvn*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)))+0.5*(cosSqThOdd*cosSqThEvn+sinSqThOdd*sinSqThEvn+cosSqThOdd*sinSqThEvn*cos(2*k1267*dm2*(kFDDistance/kKmUnits)/energy)+sinSqThOdd*cosSqThEvn*cos(2*k1267*dm2*(kFDDistance/kKmUnits)/energy));
00091 //  return 0.5*(1-sinSq2ThOdd*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy))+1-sinSq2ThEvn*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)));
00092 }

double NuOscProbCalc::OscillationWeightMu2MuBar ( double  energy,
double  dm2,
double  sn2,
double  g_mu,
double  g_tau,
double  c_mumu,
double  c_tautau 
)

Definition at line 174 of file NuOscProbCalc.cxx.

References c_factor, g_factor, k1267, kFDDistance, kKmUnits, and SQR.

Referenced by NuMatrixSpectrum::OscillateMu2MuBar(), and NuTransSME::OscillationWeightTrans().

00182 {
00183 //  cout<<"bmu="<<g_mu<<", cmm="<<c_factor*c_mumu<<endl;
00184   //double Theta=0.5*asin(sqrt(sn2));
00185 
00186   double sin2theta = sqrt(sn2);
00187   double cos2theta = cos(asin(sqrt(sn2)));
00188 
00189 
00190   Double_t Theta_Odd = 0.5*(atan2(dm2*sin2theta,((g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00191   Double_t Theta_Evn = 0.5*(atan2(dm2*sin2theta,((-g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00192 
00193   Double_t sin2ThOdd = sin(atan2(dm2*sin2theta,((g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00194   Double_t sin2ThEvn = sin(atan2(dm2*sin2theta,((-g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00195 
00196   double sinSq2ThOdd = sin2ThOdd*sin2ThOdd;
00197   double sinSq2ThEvn = sin2ThEvn*sin2ThEvn;
00198   //double cosSq2ThOdd = 1-sinSq2ThOdd;
00199   //double cosSq2ThEvn = 1-sinSq2ThEvn;
00200 
00201 
00202   Double_t sinThOdd = sin(Theta_Odd);
00203   Double_t sinThEvn = sin(Theta_Evn);
00204 
00205   Double_t sinSqThOdd = sinThOdd*sinThOdd;
00206   Double_t cosSqThOdd = 1-sinSqThOdd;
00207   Double_t sinSqThEvn = sinThEvn*sinThEvn;
00208   Double_t cosSqThEvn = 1-sinSqThEvn;
00209 
00210   //nu to nubar.....
00211 
00212   //cout<<"energy:"<<energy<<"dm2="<<dm2<<", sn2="<<sn2<<", prob numu2numubar="<<0.25*(1-sinSq2ThOdd*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy))+1-sinSq2ThEvn*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)))-0.5*(cosSqThOdd*cosSqThEvn+sinSqThOdd*sinSqThEvn+cosSqThOdd*sinSqThEvn*cos(2*k1267*dm2*(kFDDistance/kKmUnits)/energy)+sinSqThOdd*cosSqThEvn*cos(2*k1267*dm2*(kFDDistance/kKmUnits)/energy))<<endl;
00213                
00214   return 0.25*(1-sinSq2ThOdd*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy))+1-sinSq2ThEvn*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)))-0.5*(cosSqThOdd*cosSqThEvn+sinSqThOdd*sinSqThEvn+cosSqThOdd*sinSqThEvn*cos(2*k1267*dm2*(kFDDistance/kKmUnits)/energy)+sinSqThOdd*cosSqThEvn*cos(2*k1267*dm2*(kFDDistance/kKmUnits)/energy));
00215 
00216 }

double NuOscProbCalc::OscillationWeightMu2Tau ( double  energy,
double  dm2,
double  sn2,
double  g_mu,
double  g_tau,
double  c_mumu,
double  c_tautau 
)

Definition at line 94 of file NuOscProbCalc.cxx.

References c_factor, g_factor, k1267, kFDDistance, kKmUnits, and SQR.

Referenced by NuMatrixSpectrum::OscillateMu2Tau(), NuTransSME::OscillationWeight(), and NuDSTAna::OscillationWeightTrans().

00101 {
00102   //double Theta=0.5*asin(sqrt(sn2));
00103 
00104   double sin2theta = sqrt(sn2);
00105   double cos2theta = cos(asin(sqrt(sn2)));
00106 
00107 
00108   Double_t Theta_Odd = 0.5*(atan2(dm2*sin2theta,((g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00109   Double_t Theta_Evn = 0.5*(atan2(dm2*sin2theta,((-g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00110 
00111   Double_t sin2ThOdd = sin(atan2(dm2*sin2theta,((g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00112   Double_t sin2ThEvn = sin(atan2(dm2*sin2theta,((-g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00113 
00114   double sinSq2ThOdd = sin2ThOdd*sin2ThOdd;
00115   double sinSq2ThEvn = sin2ThEvn*sin2ThEvn;
00116   //double cosSq2ThOdd = 1-sinSq2ThOdd;
00117   //double cosSq2ThEvn = 1-sinSq2ThEvn;
00118 
00119 
00120   //Double_t sinThOdd = sin(Theta_Odd);
00121   //Double_t sinThEvn = sin(Theta_Evn);
00122 
00123   //Double_t sinSqThOdd = sinThOdd*sinThOdd;
00124   //Double_t cosSqThOdd = 1-sinSqThOdd;
00125   //Double_t sinSqThEvn = sinThEvn*sinThEvn;
00126   //Double_t cosSqThEvn = 1-sinSqThEvn;
00127 
00128   //cout<<"energy:"<<energy<<"dm2="<<dm2<<", sn2="<<sn2<<", prob numu2nutau="<<0.25*(sinSq2ThOdd*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy))+sinSq2ThEvn*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)))+0.5*(sin(2*Theta_Odd)*sin(2*Theta_Evn)*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)))<<endl;
00129 
00130   return 0.25*(sinSq2ThOdd*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy))+sinSq2ThEvn*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)))+0.5*(sin(2*Theta_Odd)*sin(2*Theta_Evn)*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)));
00131 
00132  
00133 }

double NuOscProbCalc::OscillationWeightMu2TauBar ( double  energy,
double  dm2,
double  sn2,
double  g_mu,
double  g_tau,
double  c_mumu,
double  c_tautau 
)

Definition at line 136 of file NuOscProbCalc.cxx.

References c_factor, g_factor, k1267, kFDDistance, kKmUnits, and SQR.

Referenced by NuMatrixSpectrum::OscillateMu2TauBar(), and NuTransSME::OscillationWeightTrans().

00143 {
00144   //double Theta=0.5*asin(sqrt(sn2));
00145 
00146   double sin2theta = sqrt(sn2);
00147   double cos2theta = cos(asin(sqrt(sn2)));
00148 
00149 
00150   Double_t Theta_Odd = 0.5*(atan2(dm2*sin2theta,((g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00151   Double_t Theta_Evn = 0.5*(atan2(dm2*sin2theta,((-g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00152 
00153   Double_t sin2ThOdd = sin(atan2(dm2*sin2theta,((g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00154   Double_t sin2ThEvn = sin(atan2(dm2*sin2theta,((-g_factor*(g_mu-g_tau)+c_factor*(c_mumu-c_tautau))*energy*energy*1e18 + dm2*cos2theta)));
00155 
00156   double sinSq2ThOdd = sin2ThOdd*sin2ThOdd;
00157   double sinSq2ThEvn = sin2ThEvn*sin2ThEvn;
00158   //double cosSq2ThOdd = 1-sinSq2ThOdd;
00159   //double cosSq2ThEvn = 1-sinSq2ThEvn;
00160 
00161 
00162   //Double_t sinThOdd = sin(Theta_Odd);
00163   //Double_t sinThEvn = sin(Theta_Evn);
00164   //Double_t sinSqThOdd = sinThOdd*sinThOdd;
00165   //Double_t cosSqThOdd = 1-sinSqThOdd;
00166   //Double_t sinSqThEvn = sinThEvn*sinThEvn;
00167   //Double_t cosSqThEvn = 1-sinSqThEvn;
00168 
00169   //cout<<"energy:"<<energy<<"dm2="<<dm2<<", sn2="<<sn2<<", prob numu2nutaubar="<<0.25*(sinSq2ThOdd*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy))+sinSq2ThEvn*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)))-0.5*(sin(2*Theta_Odd)*sin(2*Theta_Evn)*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)))<<endl;
00170 
00171   return 0.25*(sinSq2ThOdd*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy))+sinSq2ThEvn*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)))-0.5*(sin(2*Theta_Odd)*sin(2*Theta_Evn)*SQR(sin(k1267*dm2*(kFDDistance/kKmUnits)/energy)));
00172  }

double NuOscProbCalc::ThreeFlavourExactNuMuToNuElectron ( const double  energy,
double  dm232,
const double  theta23,
double  dm221,
const double  theta12,
const double  theta13,
double  delta,
const double  baseline,
const int  inu 
)

Definition at line 673 of file NuOscProbCalc.cxx.

References kKmUnits, OscProb::PMNS_Base::P(), OscProb::PMNS_Base::PropMatter(), OscProb::PMNS_Base::ResetToFlavour(), OscProb::PMNS_Fast::SetDeltaMsqrs(), and OscProb::PMNS_Fast::SetMix().

00680 {
00681   //cout << "baseline = " << baseline << ", inu " << inu << endl;
00682   //cout << " " << endl;
00683 
00684   //cout << "energy = " << energy << ", dm232 = " << dm232 << ", theta23 = " 
00685   //<< theta23 << ", dm221 = " << dm221 << ", theta12 = " << theta12 << ", theta13 = " 
00686   //<< theta13 << ", delta = " << delta << endl;
00687 
00688   double density = 2.71/2; 
00689 
00690   double nuAntiNu = 0;
00691   if (12 == inu){nuAntiNu = 1;}
00692   else if (-12 == inu){nuAntiNu = -1;}
00693   else {
00694     cout << "You're using the numu to nutau oscillation probability "
00695          << "on something that isn't a nutau or a numubar"
00696          << endl;
00697     assert(false);
00698   }
00699 
00700   //mark slow PMNS:
00701 //   OscProb::PMNS fPMNS;
00702 //   fPMNS.Reset();
00703 //   fPMNS.SetDeltaMsqrs(dm221, dm232);
00704 //   fPMNS.SetMix(theta12, theta23, theta13,delta);
00705 //   fPMNS.Reset();
00706 //   fPMNS.PropMatter(baseline/kKmUnits, energy, density, nuAntiNu);
00707 //   double prob = fPMNS.P(1,0);
00708 
00709   //joao fast PMNS:
00710   OscProb::PMNS_Fast myPMNS_Fast;
00711   myPMNS_Fast.SetDeltaMsqrs(dm221, dm232);
00712   myPMNS_Fast.SetMix(theta12,theta23,theta13,delta);
00713   myPMNS_Fast.ResetToFlavour(1);
00714   myPMNS_Fast.PropMatter(baseline/kKmUnits, energy, density, nuAntiNu);
00715   double prob = myPMNS_Fast.P(0);
00716 
00717   if(TMath::IsNaN(prob)==1){
00718     prob=0;
00719   }
00720 
00721   return prob;
00722 }

double NuOscProbCalc::ThreeFlavourExactNuMuToNuMu ( const double  energy,
double  dm232,
const double  theta23,
double  dm221,
const double  theta12,
const double  theta13,
double  delta,
const double  baseline,
const int  inu 
)

Definition at line 567 of file NuOscProbCalc.cxx.

References kKmUnits, OscProb::PMNS_Base::P(), OscProb::PMNS_Base::PropMatter(), OscProb::PMNS_Base::ResetToFlavour(), OscProb::PMNS_Fast::SetDeltaMsqrs(), and OscProb::PMNS_Fast::SetMix().

00574 {
00575 
00576   //cout << "baseline = " << baseline << ", inu " << inu << endl;
00577   //cout << " " << endl;
00578 
00579   //cout << "energy = " << energy << ", dm232 = " << dm232 << ", theta23 = " << theta23 << ", dm221 = " << dm221 << ", theta12 = " << theta12 << ", theta13 = " << theta13 << ", delta = " << delta << endl;
00580 
00581   double density = 2.71/2; 
00582 
00583   //cout << "inu = " << inu << endl;
00584   double nuAntiNu = 0;
00585   if (14 == inu){nuAntiNu = 1;}
00586   else if (-14 == inu){nuAntiNu = -1;}
00587   else {
00588     cout << "You're using the numu to numu oscillation probability "
00589          << "on something that isn't a numu or a numubar"
00590          << endl;
00591     assert(false);
00592   }
00593 
00594   //Messier Slow:
00595 //   OscProb::PMNS fPMNS;
00596 //   fPMNS.Reset();
00597 //   fPMNS.SetDeltaMsqrs(dm221, dm232);
00598 //   fPMNS.SetMix(theta12, theta23, theta13,delta);
00599 //   fPMNS.Reset();
00600 //   fPMNS.PropMatter(baseline/kKmUnits, energy, density, nuAntiNu);
00601 //   double prob = fPMNS.P(1,1);
00602   
00603   //joao fast PMNS:
00604   OscProb::PMNS_Fast myPMNS_Fast;
00605   myPMNS_Fast.SetDeltaMsqrs(dm221, dm232);
00606   myPMNS_Fast.SetMix(theta12,theta23,theta13,delta);
00607   myPMNS_Fast.ResetToFlavour(1);
00608   myPMNS_Fast.PropMatter(baseline/kKmUnits, energy, density, nuAntiNu);
00609   double prob = myPMNS_Fast.P(1);
00610 
00611   if(TMath::IsNaN(prob)==1){
00612     prob=0.0;
00613   }
00614 
00615   return prob;
00616 
00617 }

double NuOscProbCalc::ThreeFlavourExactNuMuToNuTau ( const double  energy,
double  dm232,
const double  theta23,
double  dm221,
const double  theta12,
const double  theta13,
double  delta,
const double  baseline,
const int  inu 
)

Definition at line 620 of file NuOscProbCalc.cxx.

References kKmUnits, OscProb::PMNS_Base::P(), OscProb::PMNS_Base::PropMatter(), OscProb::PMNS_Base::ResetToFlavour(), OscProb::PMNS_Fast::SetDeltaMsqrs(), and OscProb::PMNS_Fast::SetMix().

00627 {
00628   //cout << "baseline = " << baseline << ", inu " << inu << endl;
00629   //cout << " " << endl;
00630 
00631   //cout << "energy = " << energy << ", dm232 = " << dm232 << ", theta23 = " 
00632   //<< theta23 << ", dm221 = " << dm221 << ", theta12 = " << theta12 << ", theta13 = " 
00633   //<< theta13 << ", delta = " << delta << endl;
00634 
00635   double density = 2.71/2;
00636 
00637   double nuAntiNu = 0;
00638   if (16 == inu){nuAntiNu = 1;}
00639   else if (-16 == inu){nuAntiNu = -1;}
00640   else {
00641     cout << "You're using the numu to nutau oscillation probability "
00642          << "on something that isn't a nutau or a numubar"
00643          << endl;
00644     assert(false);
00645   }
00646 
00647   //mark slow PMNS:
00648 //   OscProb::PMNS fPMNS;
00649 //   fPMNS.Reset();
00650 //   fPMNS.SetDeltaMsqrs(dm221, dm232);
00651 //   fPMNS.SetMix(theta12, theta23, theta13,delta);
00652 //   fPMNS.Reset();
00653 //   fPMNS.PropMatter(baseline/kKmUnits, energy, density, nuAntiNu);
00654 //   double prob = fPMNS.P(1,2);
00655 
00656   //joao fast PMNS:
00657   OscProb::PMNS_Fast myPMNS_Fast;
00658   myPMNS_Fast.SetDeltaMsqrs(dm221, dm232);
00659   myPMNS_Fast.SetMix(theta12,theta23,theta13,delta);
00660   myPMNS_Fast.ResetToFlavour(1);
00661   myPMNS_Fast.PropMatter(baseline/kKmUnits, energy, density, nuAntiNu);
00662   double prob = myPMNS_Fast.P(2);
00663 
00664 
00665   if(TMath::IsNaN(prob)==1){
00666     prob=0;
00667   }
00668 
00669   return prob;
00670 }

double NuOscProbCalc::ThreeFlavourNuMuToNuElectron ( const double  energy,
double  dm232,
const double  theta23,
double  dm221,
const double  theta12,
const double  theta13,
double  delta,
const double  baseline,
const int  inu 
)

Definition at line 512 of file NuOscProbCalc.cxx.

References MuELoss::e, kKmUnits, OscCalc::Oscillate(), and OscCalc::SetOscParam().

00519 {
00520   if (delta < 1e-8 && delta >= 0){delta = 1e-8;}
00521   if (delta <= 0 && delta > -1e-8){delta = -1e-8;}
00522 
00523   if (dm221 < 1e-8 && dm221 >= 0){dm221 = 1e-8;}
00524   if (dm221 <= 0 && dm221 > -1e-8){dm221 = -1e-8;}
00525 
00526   if (dm232 < 1e-8 && dm232 >= 0){dm232 = 1e-8;}
00527   if (dm232 <= 0 && dm232 > -1e-8){dm232 = -1e-8;}
00528 
00529   //cout << "baseline = " << baseline << ", inu " << inu << endl;
00530   //cout << " " << endl;
00531 
00532   //cout << "energy = " << energy << ", dm232 = " << dm232 << ", theta23 = " 
00533   //<< theta23 << ", dm221 = " << dm221 << ", theta12 = " << theta12 << ", theta13 = " 
00534   //<< theta13 << ", delta = " << delta << endl;
00535 
00536   double density = 2.71; //If you want to turn this off, set it to
00537                          //something small but non-zero
00538   double nuAntiNu = 0;
00539   if (12 == inu){nuAntiNu = 1;}
00540   else if (-12 == inu){nuAntiNu = -1;}
00541   else {
00542     cout << "You're using the numu to nutau oscillation probability "
00543          << "on something that isn't a nutau or a numubar"
00544          << endl;
00545     assert(false);
00546   }
00547 
00548   double par[9] = {theta12, theta23, theta13, dm232,
00549                    dm221, delta, density, baseline/kKmUnits,
00550                    nuAntiNu};
00551 
00552 //   for (Int_t parameter = 0; parameter < 9; ++parameter){
00553 //     cout << "Parameter " << parameter
00554 //       << ": " << par[parameter]
00555 //       << endl;
00556 //   }
00557 //   assert(false);
00558 
00559   static OscCalc oscCalc;
00560   oscCalc.SetOscParam(par);
00561 
00562   return oscCalc.Oscillate(12, 14, energy);
00563 
00564 }

double NuOscProbCalc::ThreeFlavourNuMuToNuMu ( const double  energy,
double  dm232,
const double  theta23,
double  dm221,
const double  theta12,
const double  theta13,
double  delta,
const double  baseline,
const int  inu 
)

Definition at line 401 of file NuOscProbCalc.cxx.

References MuELoss::e, kKmUnits, OscCalc::Oscillate(), and OscCalc::SetOscParam().

00408 {
00409   if (delta < 1e-8 && delta >= 0){delta = 1e-8;}
00410   if (delta <= 0 && delta > -1e-8){delta = -1e-8;}
00411 
00412   if (dm221 < 1e-8 && dm221 >= 0){dm221 = 1e-8;}
00413   if (dm221 <= 0 && dm221 > -1e-8){dm221 = -1e-8;}
00414 
00415   if (dm232 < 1e-8 && dm232 >= 0){dm232 = 1e-8;}
00416   if (dm232 <= 0 && dm232 > -1e-8){dm232 = -1e-8;}
00417 
00418   //cout << "baseline = " << baseline << ", inu " << inu << endl;
00419   //cout << " " << endl;
00420 
00421   //cout << "energy = " << energy << ", dm232 = " << dm232 << ", theta23 = " << theta23 << ", dm221 = " << dm221 << ", theta12 = " << theta12 << ", theta13 = " << theta13 << ", delta = " << delta << endl;
00422 
00423   double density = 2.71; //If you want to turn this off, set it to
00424                          //something small but non-zero
00425 
00426   //cout << "inu = " << inu << endl;
00427   double nuAntiNu = 0;
00428   if (14 == inu){nuAntiNu = 1;}
00429   else if (-14 == inu){nuAntiNu = -1;}
00430   else {
00431     cout << "You're using the numu to numu oscillation probability "
00432          << "on something that isn't a numu or a numubar"
00433          << endl;
00434     assert(false);
00435   }
00436 
00437   double par[9] = {theta12, theta23, theta13, dm232,
00438                    dm221, delta, density, baseline/kKmUnits,
00439                    nuAntiNu};
00440 
00441 //   for (Int_t parameter = 0; parameter < 9; ++parameter){
00442 //     cout << "Parameter " << parameter
00443 //       << ": " << par[parameter]
00444 //       << endl;
00445 //   }
00446 //   assert(false);
00447 
00448   static OscCalc oscCalc;
00449   oscCalc.SetOscParam(par);
00450 
00451   return oscCalc.Oscillate(14, 14, energy);
00452 
00453 }

double NuOscProbCalc::ThreeFlavourNuMuToNuTau ( const double  energy,
double  dm232,
const double  theta23,
double  dm221,
const double  theta12,
const double  theta13,
double  delta,
const double  baseline,
const int  inu 
)

Definition at line 457 of file NuOscProbCalc.cxx.

References MuELoss::e, kKmUnits, OscCalc::Oscillate(), and OscCalc::SetOscParam().

00464 {
00465   if (delta < 1e-8 && delta >= 0){delta = 1e-8;}
00466   if (delta <= 0 && delta > -1e-8){delta = -1e-8;}
00467 
00468   if (dm221 < 1e-8 && dm221 >= 0){dm221 = 1e-8;}
00469   if (dm221 <= 0 && dm221 > -1e-8){dm221 = -1e-8;}
00470 
00471   if (dm232 < 1e-8 && dm232 >= 0){dm232 = 1e-8;}
00472   if (dm232 <= 0 && dm232 > -1e-8){dm232 = -1e-8;}
00473 
00474   //cout << "baseline = " << baseline << ", inu " << inu << endl;
00475   //cout << " " << endl;
00476 
00477   //cout << "energy = " << energy << ", dm232 = " << dm232 << ", theta23 = " 
00478   //<< theta23 << ", dm221 = " << dm221 << ", theta12 = " << theta12 << ", theta13 = " 
00479   //<< theta13 << ", delta = " << delta << endl;
00480 
00481   double density = 2.71; //If you want to turn this off, set it to
00482                          //something small but non-zero
00483   double nuAntiNu = 0;
00484   if (16 == inu){nuAntiNu = 1;}
00485   else if (-16 == inu){nuAntiNu = -1;}
00486   else {
00487     cout << "You're using the numu to nutau oscillation probability "
00488          << "on something that isn't a nutau or a numubar"
00489          << endl;
00490     assert(false);
00491   }
00492 
00493   double par[9] = {theta12, theta23, theta13, dm232,
00494                    dm221, delta, density, baseline/kKmUnits,
00495                    nuAntiNu};
00496 
00497 //   for (Int_t parameter = 0; parameter < 9; ++parameter){
00498 //     cout << "Parameter " << parameter
00499 //       << ": " << par[parameter]
00500 //       << endl;
00501 //   }
00502 //   assert(false);
00503 
00504   static OscCalc oscCalc;
00505   oscCalc.SetOscParam(par);
00506 
00507   return oscCalc.Oscillate(16, 14, energy);
00508 
00509 }


Generated on 19 Jan 2018 for loon by  doxygen 1.6.1