OscCalc Class Reference

#include <OscCalc.h>

List of all members.

Public Member Functions

 OscCalc ()
void GetOscParam (double *par)
double GetOscParam (OscPar::OscPar_t pos)
void SetOscParam (double *par, bool force=false)
void SetOscParam (OscPar::OscPar_t pos, double val, bool force=false)
double Oscillate (int nuFlavor, int nonOscNuFlavor, double Energy)
double Oscillate (int nuFlavor, int nonOscNuFlavor, double Energy, double *par)
double MuToElec (double x)
double MuToTau (double x)
double MuToMu (double x)
double ElecToElec (double x)
double ElecToMu (double x)
double ElecToTau (double x)
double TauToElec (double x)
double TauToTau (double x)
double TauToMu (double x)
double OscillateNSI (int nuFlavor, int nonOscNuFlavor, double E)

Private Member Functions

void UpdateNSI ()

Private Attributes

OscProb::PMNS_NSI myPMNS_NSI
bool fMixDirty
bool fDmSqrDirty
bool fNSIDirty
double fOscPar [OscPar::kNumParameters]
double fSinTh [3]
double fCosTh [3]
double fSin2Th [3]
double fCos2Th [3]
double fElecDensity
double fV
double fSinDCP
double fCosDCP


Detailed Description

Definition at line 61 of file OscCalc.h.


Constructor & Destructor Documentation

OscCalc::OscCalc (  ) 

Definition at line 14 of file OscCalc.cxx.

References OscPar::kNumParameters, and SetOscParam().

00015 {
00016   double par[OscPar::kNumParameters] = {0}; 
00017   SetOscParam(par, true);
00018 }


Member Function Documentation

double OscCalc::ElecToElec ( double  x  ) 

Definition at line 525 of file OscCalc.cxx.

References MuELoss::e, fCos2Th, fOscPar, fSin2Th, fSinDCP, fSinTh, fV, OscPar::invKmToeV, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kL, OscPar::kNuAntiNu, OscPar::kTh12, and OscPar::kTh13.

Referenced by Oscillate(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), and Extrapolate2D::SetNominalOscProb().

00526 {
00527   double E = x; //energy
00528                                                                                                                        
00529   //Building the standard terms
00530   double L = fOscPar[OscPar::kL];  //baseline
00531   double dmsq_23 = fOscPar[OscPar::kDeltaM23];
00532   double dmsq_12 = fOscPar[OscPar::kDeltaM12];
00533   double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}
00534                                                                                                                        
00535   double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
00536   double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
00537                                                                                                                        
00538   double sin_th12 = fSinTh[OscPar::kTh12];
00539  
00540 //  double cos_2th23 = fCos2Th[OscPar::kTh23];
00541   double cos_2th13 = fCos2Th[OscPar::kTh13];
00542   double cos_2th12 = fCos2Th[OscPar::kTh12];
00543                                                                                                                        
00544   double d_cp = fOscPar[OscPar::kDelta];
00545   double sin_dcp = fSinDCP;
00546                                                                                                                        
00547   //Building the more complicated terms
00548   double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
00549   double A = 2*fV*E*1e9/(dmsq_13);
00550   double alpha = dmsq_12/dmsq_13;
00551                                                                                                                        
00552   // A and d_cp both change sign for antineutrinos
00553   double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
00554   A *= plusminus;
00555   d_cp *= plusminus;
00556   sin_dcp *= plusminus;
00557                                                                                                                        
00558   //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)
00559   double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));
00560                                                                                                                        
00561   double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
00562   if(fabs(alpha) > 1e-10){  //want to be careful here
00563     double temp = cos_2th12 - A/alpha;
00564     C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
00565   }
00566                                                                                                                        
00567   //More complicated sin and cosine terms
00568   double cosC13Delta = cos(C13*Delta);
00569   double sinC13Delta = sin(C13*Delta);
00570                                                                                                                        
00571   double cosaC12Delta = 0;
00572   double sinaC12Delta = 0;
00573                                                                                                                        
00574   if(fabs(alpha) > 1e-10){
00575     cosaC12Delta = cos(alpha*C12*Delta);
00576     sinaC12Delta = sin(alpha*C12*Delta);
00577   } // otherwise not relevant
00578                                                                                                                        
00579   //First we calculate the terms for the alpha expansion (good to all orders in th13)
00580   // this is the equivalent of Eq 45 & 46 corrected for Mu to E instead of E to Mu
00581                                                                                                                        
00582   // Leading order term
00583   double p1 = 1 - sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);
00584                                                                                                                        
00585   // terms that appear at order alpha
00586   double p2Inner = Delta*cosC13Delta*(1-A*cos_2th13)/C13 -
00587                       A*sinC13Delta*(cos_2th13-A)/(C13*C13);
00588                                                                                                                        
00589   double p2 = +2*sin_th12*sin_th12*sinsq_2th13*sinC13Delta/(C13*C13)*p2Inner*alpha;
00590   //  p1 + p2 is the complete contribution for this expansion
00591                                                                                                                        
00592   // Now for the expansion in orders of sin(th13) (good to all order alpha)
00593   //  this is the equivalent of Eq 63 and 64
00594                                                                                                                        
00595   // leading order term
00596   double pa1 = 1.0, pa2 = 0.0;
00597                                                                                                                        
00598   if(fabs(alpha) > 1e-10){
00599     // leading order term
00600     pa1 = 1.0 - sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12);
00601   }
00602   //pa1 is the complete contribution from this expansion, there is no order s13^1 term
00603                                                                                                                        
00604   // In order to combine the information correctly we need to add the two
00605   //  expansions and subtract off the terms that are in both (alpha^1, s13^1)
00606   //  these may be taken from the expansion to second order in both parameters
00607   //  Equation 30
00608                                                                                                                        
00609   double repeated = 1;
00610                                                                                                                        
00611   //  Calculate the total probability
00612   double totalP = p1+p2 + (pa1+pa2) - repeated;
00613                                                                                                                        
00614   return totalP;
00615 }

double OscCalc::ElecToMu ( double  x  ) 

Definition at line 507 of file OscCalc.cxx.

References fOscPar, fSinDCP, OscPar::kDelta, and MuToElec().

Referenced by ElecToTau(), Oscillate(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), and Extrapolate2D::SetNominalOscProb().

00508 {
00509   // Flip delta to reverse direction
00510   double oldSinDelta = fSinDCP;
00511   double oldDelta =  fOscPar[OscPar::kDelta];
00512 
00513   fOscPar[OscPar::kDelta] = -oldDelta;
00514   fSinDCP = -oldSinDelta;
00515                                                                                                                        
00516   double prob = OscCalc::MuToElec(x);
00517                                                                                                                        
00518   //restore the world
00519   fOscPar[OscPar::kDelta] = oldDelta;
00520   fSinDCP = oldSinDelta;
00521  
00522   return prob;   
00523 }

double OscCalc::ElecToTau ( double  x  ) 

Definition at line 480 of file OscCalc.cxx.

References ElecToMu(), fCosTh, fSin2Th, fSinTh, and OscPar::kTh23.

Referenced by Oscillate(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), Extrapolate2D::SetNominalOscProb(), and TauToElec().

00481 {
00482 //  EtoTau is the same as E->Mu wich sinsq_th23 <-> cossq_th23, sin(2th23) <->-sin(2th23)
00483   double origCos = fCosTh[OscPar::kTh23];
00484   double origSin = fSinTh[OscPar::kTh23];
00485   double orig2Sin = fSin2Th[OscPar::kTh23];
00486 
00487   fCosTh[OscPar::kTh23] = origSin;
00488   fSinTh[OscPar::kTh23] = origCos;
00489   fSin2Th[OscPar::kTh23] = -orig2Sin;
00490 
00491   double prob = OscCalc::ElecToMu(x);
00492 
00493   //restore the world
00494   fCosTh[OscPar::kTh23] = origCos;
00495   fSinTh[OscPar::kTh23] = origSin;
00496   fSin2Th[OscPar::kTh23] = orig2Sin;
00497 
00498   return prob;
00499 
00500 /*  This is an option as well but I think results in more cycles
00501   double p1 = 1. - OscCalc::ElecToMu(x) - OscCalc::ElecToElec(x);
00502   if(p1 < 0){ cout<<"Damnation! "<<x[0]<<" "<<p1<<endl; p1 = 0;}
00503   return p1;
00504 */ 
00505 }

double OscCalc::GetOscParam ( OscPar::OscPar_t  pos  ) 

Definition at line 147 of file OscCalc.cxx.

References fOscPar, and OscPar::kNumParameters.

00148 {
00149    if(pos >=0 && pos < OscPar::kNumParameters)
00150      return fOscPar[pos];
00151    else 
00152      return -9999;
00153 }

void OscCalc::GetOscParam ( double *  par  ) 

Definition at line 140 of file OscCalc.cxx.

References fOscPar, and OscPar::kNumParameters.

Referenced by NueStandard::GetOscParam(), NueSystematic::GetOscParams(), Extrapolate2D::GetPrediction(), Extrapolate2D::InvertMassHierarchy(), Extrapolate2D_Simple::OscillatePrediction(), Extrapolate2D::OscillatePrediction(), and Extrapolate2D::ReadExtrap().

00141 {
00142    for(int i = 0; i < int(OscPar::kNumParameters); i++){
00143      par[i] = fOscPar[i];
00144    }
00145 }

double OscCalc::MuToElec ( double  x  ) 

Definition at line 156 of file OscCalc.cxx.

References MuELoss::e, fCos2Th, fCosDCP, fCosTh, fOscPar, fSin2Th, fSinDCP, fSinTh, fV, OscPar::invKmToeV, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kL, OscPar::kNuAntiNu, OscPar::kTh12, OscPar::kTh13, and OscPar::kTh23.

Referenced by ElecToMu(), MuToMu(), Oscillate(), Extrapolate2D_Simple::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D_Simple::OscillatePrediction_SepNuNuBar(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), and Extrapolate2D::SetNominalOscProb().

00157 {
00158   double E = x; //energy
00159 
00160   //Building the standard terms
00161   double L = fOscPar[OscPar::kL];  //baseline
00162   double dmsq_23 = fOscPar[OscPar::kDeltaM23];
00163   double dmsq_12 = fOscPar[OscPar::kDeltaM12];
00164   double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}
00165 
00166   double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
00167   double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
00168                                                                                 
00169   double cos_th23 = fCosTh[OscPar::kTh23];
00170   double cos_th12 = fCosTh[OscPar::kTh12];
00171   double sin_th13 = fSinTh[OscPar::kTh13];
00172   double cos_th13 = fCosTh[OscPar::kTh13];  
00173                                                                               
00174   double sin_2th23 = fSin2Th[OscPar::kTh23];
00175   double sin_2th12 = fSin2Th[OscPar::kTh12];
00176   double cos_2th13 = fCos2Th[OscPar::kTh13];
00177   double cos_2th12 = fCos2Th[OscPar::kTh12];
00178 
00179   double sinsq_th23 = fSinTh[OscPar::kTh23]*fSinTh[OscPar::kTh23];
00180   double sinsq_th12 = fSinTh[OscPar::kTh12]*fSinTh[OscPar::kTh12];
00181 
00182   double d_cp = fOscPar[OscPar::kDelta];
00183   double cos_dcp = fCosDCP;
00184   double sin_dcp = fSinDCP;
00185 
00186   //Building the more complicated terms
00187   double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
00188   double A = 2*fV*E*1e9/(dmsq_13);
00189   double alpha = dmsq_12/dmsq_13;
00190 
00191   // A and d_cp both change sign for antineutrinos
00192   double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
00193   A *= plusminus;
00194   d_cp *= plusminus;
00195   sin_dcp *= plusminus;
00196 
00197   //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)
00198 
00199   double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));
00200 
00201   double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
00202   if(fabs(alpha) > 1e-10){  //want to be careful here
00203     double temp = cos_2th12 - A/alpha;
00204     C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
00205   }
00206 
00207   //More complicated sin and cosine terms
00208   double cosC13Delta = cos(C13*Delta);
00209   double sinC13Delta = sin(C13*Delta);
00210   double sin1pADelta = sin((A+1)*Delta);
00211   double cos1pADelta = cos((A+1)*Delta);
00212   double sinADelta = sin(A*Delta);
00213   double sinAm1Delta = sin((A-1)*Delta);
00214   double cosdpDelta = cos(d_cp+Delta);
00215   double sinApam2Delta = sin((A+alpha-2)*Delta);
00216   double cosApam2Delta = cos((A+alpha-2)*Delta);
00217 
00218   double cosaC12Delta = 0;
00219   double sinaC12Delta = 0; 
00220   
00221   if(fabs(alpha) > 1e-10){
00222     cosaC12Delta = cos(alpha*C12*Delta);
00223     sinaC12Delta = sin(alpha*C12*Delta);
00224   } // otherwise not relevant
00225 
00226   //First we calculate the terms for the alpha expansion (good to all orders in th13)
00227   // this is the equivalent of Eq 47 & 48 corrected for Mu to E instead of E to Mu
00228 
00229   // Leading order term 
00230   double p1 = sinsq_th23*sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);
00231 
00232   // terms that appear at order alpha
00233   //first work out the vacuum case since we get 0/0 otherwise.......
00234   double p2Inner = Delta*cosC13Delta;
00235 
00236   if(fabs(A) > 1e-9)
00237     p2Inner = Delta*cosC13Delta*(1-A*cos_2th13)/C13 -
00238                       A*sinC13Delta*(cos_2th13-A)/(C13*C13);
00239 
00240   double p2 = -2*sinsq_th12*sinsq_th23*sinsq_2th13*sinC13Delta/(C13*C13)*p2Inner*alpha;
00241 
00242 
00243   //again working out vacuum first.....
00244   double p3Inner = Delta* cos_th13* cos_th13*(-2*sin_dcp*sinC13Delta*sinC13Delta+2*cos_dcp*sinC13Delta*cosC13Delta);
00245 
00246   if(fabs(A) > 1e-9)
00247     p3Inner = (sinC13Delta/(A*C13*C13))*(- sin_dcp*(cosC13Delta - cos1pADelta)*C13
00248                          + cos_dcp*(C13*sin1pADelta - (1-A*cos_2th13)*sinC13Delta));
00249 
00250   double p3 = sin_2th12*sin_2th23*sin_th13*p3Inner*alpha;
00251 
00252   //  p1 + p2 + p3 is the complete contribution for this expansion
00253   
00254   // Now for the expansion in orders of sin(th13) (good to all order alpha) 
00255   //  this is the equivalent of Eq 65 and 66
00256 
00257   // leading order term
00258   double pa1 = 0.0, pa2 = 0.0;
00259 
00260   // no problems here when A -> 0
00261   if(fabs(alpha) > 1e-10){
00262     // leading order term
00263     pa1 = cos_th23*cos_th23*sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12);
00264 
00265     // and now to calculate the first order in s13 term
00266     double t1 = (cos_2th12 - A/alpha)/C12 
00267                   - alpha*A*C12*sinsq_2th12/(2*(1-alpha)*C12*C12);
00268     double t2 = -cos_dcp*(sinApam2Delta-sinaC12Delta*t1);
00269   
00270     double t3 = -(cosaC12Delta-cosApam2Delta)*sin_dcp;
00271  
00272     double denom = (1-A-alpha+A*alpha*cos_th12*cos_th12)*C12;
00273     double t4 = sin_2th12*sin_2th23*(1-alpha)*sinaC12Delta/denom;
00274 
00275     pa2 = t4*(t3+t2)*sin_th13;
00276   }
00277   //pa1+pa2 is the complete contribution from this expansion
00278 
00279   // In order to combine the information correctly we need to add the two
00280   //  expansions and subtract off the terms that are in both (alpha^1, s13^1) 
00281   //  these may be taken from the expansion to second order in both parameters
00282   //  Equation 31 
00283 
00284   double t1 = Delta*sinC13Delta*cosdpDelta;
00285   if(fabs(A) > 1e-9) t1 = sinADelta*cosdpDelta*sinAm1Delta/(A*(A-1));
00286 
00287   double repeated = 2*alpha*sin_2th12*sin_2th23*sin_th13*t1;
00288 
00289   //  Calculate the total probability
00290   double totalP = p1+p2+p3 + (pa1+pa2) - repeated;
00291 
00292   return totalP;
00293 }

double OscCalc::MuToMu ( double  x  ) 

Definition at line 473 of file OscCalc.cxx.

References MuELoss::e, MuToElec(), and MuToTau().

Referenced by Oscillate(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), Extrapolate2D::SetNominalOscProb(), and TauToTau().

00474 {
00475   double p1 = 1. - OscCalc::MuToTau(x) - OscCalc::MuToElec(x);
00476   if(p1 < 1e-4){ cout<<"P(mu->mu) less than zero Damnation! "<<x<<" "<<p1<<endl; p1 = 0;}
00477   return p1;
00478 }

double OscCalc::MuToTau ( double  x  ) 

Definition at line 295 of file OscCalc.cxx.

References MuELoss::e, fCos2Th, fCosDCP, fCosTh, fOscPar, fSin2Th, fSinDCP, fSinTh, fV, OscPar::invKmToeV, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kL, OscPar::kNuAntiNu, OscPar::kTh12, OscPar::kTh13, and OscPar::kTh23.

Referenced by MuToMu(), Oscillate(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), Extrapolate2D::SetNominalOscProb(), and TauToMu().

00296 {
00297   double E = x; //energy
00298 
00299   double L = fOscPar[OscPar::kL];  //baseline
00300   double dmsq_23 = fOscPar[OscPar::kDeltaM23];
00301   double dmsq_12 = fOscPar[OscPar::kDeltaM12];
00302   double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}
00303 
00304   double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
00305   double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
00306   double sinsq_2th23 = fSin2Th[OscPar::kTh23]*fSin2Th[OscPar::kTh23];
00307                                                                                 
00308   double cos_th13 = fCosTh[OscPar::kTh13];
00309   double cos_th12 = fCosTh[OscPar::kTh12];
00310   double sin_th12 = fSinTh[OscPar::kTh12];
00311   double sin_th13 = fSinTh[OscPar::kTh13];
00312                                                                                 
00313   double sin_2th23 = fSin2Th[OscPar::kTh23];
00314 //  double sin_2th13 = fSin2Th[OscPar::kTh13];
00315   double sin_2th12 = fSin2Th[OscPar::kTh12];
00316                                                                                 
00317   double cos_2th23 = fCos2Th[OscPar::kTh23];
00318   double cos_2th13 = fCos2Th[OscPar::kTh13];
00319   double cos_2th12 = fCos2Th[OscPar::kTh12];
00320 
00321   double d_cp = fOscPar[OscPar::kDelta];
00322   double cos_dcp = fCosDCP;
00323   double sin_dcp = fSinDCP;
00324 
00325   //Building the more complicated terms                                                                              
00326   double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
00327   double A = 2*fV*E*1e9/(dmsq_13);
00328   double alpha = dmsq_12/dmsq_13;
00329   
00330   // A and d_cp both change sign for antineutrinos
00331   double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
00332   A *= plusminus;
00333   d_cp *= plusminus;
00334   sin_dcp *= plusminus;
00335 
00336   //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)
00337                                                                                                                        
00338   double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));
00339                                                                                                                        
00340   double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
00341   if(fabs(alpha) > 1e-10){  //want to be careful here
00342     double temp = cos_2th12 - A/alpha;
00343     C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
00344   }
00345                                                                                                                        
00346   //More complicated sin and cosine terms
00347   double cosC13Delta = cos(C13*Delta);
00348   double sinC13Delta = sin(C13*Delta);
00349   double sin1pADelta = sin((A+1)*Delta);
00350   double cos1pADelta = cos((A+1)*Delta);
00351   double sinADelta = sin((A)*Delta);
00352 
00353   double sin1pAmCDelta = sin(0.5*(A+1-C13)*Delta);
00354   double sin1pApCDelta = sin(0.5*(A+1+C13)*Delta);
00355 
00356   double cosaC12Delta = 0;
00357   double sinaC12Delta = 0;
00358                                                                                                                        
00359   if(fabs(alpha) > 1e-10){
00360     cosaC12Delta = cos(alpha*C12*Delta);
00361     sinaC12Delta = sin(alpha*C12*Delta);
00362   } // otherwise not relevant
00363 
00364   double sinApam2Delta = sin((A+alpha-2)*Delta);
00365   double cosApam2Delta = cos((A+alpha-2)*Delta);
00366   double sinAm1Delta = sin((A-1)*Delta);
00367   double cosAm1Delta = cos((A-1)*Delta);
00368   double sinDelta = sin(Delta);
00369   double sin2Delta = sin(2*Delta);
00370 
00371   double cosaC12pApam2Delta = 0;                                                                                                                       
00372   if(fabs(alpha) > 1e-10){
00373     cosaC12pApam2Delta = cos((alpha*C12+A+alpha-2)*Delta);
00374   }
00375 
00376   //First we calculate the terms for the alpha expansion (good to all orders in th13)
00377   // this is the equivalent of Eq 49 & 50 corrected for Mu to E instead of E to Mu
00378 
00379   // Leading order term
00380   double pmt_0 = 0.5*sinsq_2th23;
00381   pmt_0 *= (1 - (cos_2th13-A)/C13)*sin1pAmCDelta*sin1pAmCDelta 
00382              +  (1 + (cos_2th13-A)/C13)*sin1pApCDelta*sin1pApCDelta
00383              - 0.5*sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);
00384 
00385   // terms that appear at order alpha
00386   double t0, t1, t2, t3;
00387   t0 = (cos_th12*cos_th12-sin_th12*sin_th12*sin_th13*sin_th13
00388           *(1+2*sin_th13*sin_th13*A+A*A)/(C13*C13))*cosC13Delta*sin1pADelta*2;
00389   t1 = 2*(cos_th12*cos_th12*cos_th13*cos_th13-cos_th12*cos_th12*sin_th13*sin_th13
00390          +sin_th12*sin_th12*sin_th13*sin_th13
00391          +(sin_th12*sin_th12*sin_th13*sin_th13-cos_th12*cos_th12)*A);
00392   t1 *= sinC13Delta*cos1pADelta/C13;
00393 
00394   t2 =  sin_th12*sin_th12*sinsq_2th13*sinC13Delta/(C13*C13*C13);
00395   t2 *= A/Delta*sin1pADelta+A/Delta*(cos_2th13-A)/C13*sinC13Delta
00396           - (1-A*cos_2th13)*cosC13Delta;
00397 
00398   double pmt_1 = -0.5*sinsq_2th23*Delta*(t0+t1+t2);   
00399 
00400   t0 = t1 = t2 = t3 = 0.0;
00401 
00402   t0 = cosC13Delta-cos1pADelta;
00403   t1 = 2*cos_th13*cos_th13*sin(d_cp)*sinC13Delta/C13*t0;
00404   t2 = -cos_2th23*cos_dcp*(1+A)*t0*t0;
00405 
00406   t3  = cos_2th23*cos_dcp*(sin1pADelta+(cos_2th13-A)/C13*sinC13Delta);
00407   t3 *= (1+2*sin_th13*sin_th13*A + A*A)*sinC13Delta/C13 - (1+A)*sin1pADelta;
00408 
00409 //  cout<<t1<<"  "<<t2<<"  "<<t3<<endl;
00410 
00411   if(fabs(A) > 1e-9) 
00412     pmt_1 = pmt_1 + (t1+t2+t3)*sin_th13*sin_2th12*sin_2th23/(2*A*cos_th13*cos_th13);
00413   else
00414     pmt_1 = pmt_1 + sin_th13*sin_2th12*sin_2th23*cos_th13*cos_th13*Delta*(2*sin_dcp*sinC13Delta*sinC13Delta+cos_dcp*cos_2th23*2*sinC13Delta*cosC13Delta);
00415 
00416   pmt_1 *= alpha;
00417 
00418   //  pmt_0 + pmt_1 is the complete contribution for this expansion
00419                                                                                                                        
00420   // Now for the expansion in orders of sin(th13) (good to all order alpha)
00421   //  this is the equivalent of Eq 67 and 68
00422                                                                                                                        
00423   // leading order term
00424   double pmt_a0 =  0.5*sinsq_2th23;
00425 
00426   pmt_a0 *= 1 - 0.5*sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12)
00427               - cosaC12pApam2Delta
00428               - (1 - (cos_2th12 - A/alpha)/C12)*sinaC12Delta*sinApam2Delta;
00429             
00430   double denom = (1-A-alpha+A*alpha*cos_th12*cos_th12)*C12;
00431 
00432   t0 = (cosaC12Delta-cosApam2Delta)*(cosaC12Delta-cosApam2Delta);
00433            
00434   t1 = (cos_2th12 - A/alpha)/C12*sinaC12Delta+sinApam2Delta;
00435             
00436   t2 = ((cos_2th12 - A/alpha)/C12+2*(1-alpha)/(alpha*A*C12))*sinaC12Delta
00437              + sinApam2Delta;
00438 
00439   t3 = (alpha*A*C12)/2.0*cos_2th23*cos_dcp*(t0 + t1*t2);
00440   
00441   t3 += sin_dcp*(1-alpha)*(cosaC12Delta-cosApam2Delta)*sinaC12Delta;
00442 
00443   double pmt_a1 = sin_th13*sin_2th12*sin_2th23/denom*t3;
00444 
00445   // pmt_a1+pmt_a2 is the complete contribution from this expansion
00446                                                                                                                        
00447   // In order to combine the information correctly we need to add the two
00448   //  expansions and subtract off the terms that are in both (alpha^1, s13^1)
00449   //  and lower order terms
00450   //  these may be taken from the expansion to second order in both parameters
00451   //  Equation 34
00452 
00453 
00454   // Now for the term of order alpha * s13 or lower order!
00455   t0 = t1 = t2 = t3 = 0.0;
00456 
00457   t1 = +sin_dcp*sinDelta*sinADelta*sinAm1Delta/(A*(A-1));
00458   t2 = -1/(A-1)*cos_dcp*sinDelta*(A*sinDelta-sinADelta*cosAm1Delta/A)*cos_2th23;
00459   t0 =  2*alpha*sin_2th12*sin_2th23*sin_th13*(t1+t2);
00460 
00461   t1 = sinsq_2th23*sinDelta*sinDelta 
00462        - alpha*sinsq_2th23*cos_th12*cos_th12*Delta*sin2Delta;
00463 
00464   double repeated = t0+t1;
00465 
00466   //  Calculate the total probability
00467   double totalP = pmt_0 + pmt_1 + pmt_a0 + pmt_a1 - repeated;
00468 
00469   return totalP;
00470 }

double OscCalc::Oscillate ( int  nuFlavor,
int  nonOscNuFlavor,
double  Energy,
double *  par 
)

Definition at line 20 of file OscCalc.cxx.

References Oscillate(), and SetOscParam().

00021 {
00022   SetOscParam(par);
00023   return Oscillate(nuFlavor, nonOscNuFlavor, Energy); 
00024 }

double OscCalc::Oscillate ( int  nuFlavor,
int  nonOscNuFlavor,
double  Energy 
)

Definition at line 27 of file OscCalc.cxx.

References ElecToElec(), ElecToMu(), ElecToTau(), Msg::kInfo, OscPar::kNuAntiNu, MSG, MuToElec(), MuToMu(), MuToTau(), SetOscParam(), TauToElec(), TauToMu(), and TauToTau().

Referenced by NueSystematic::DoOscCalc(), ANtpTruthInfoBeamAna::GetOscProb(), NueStandard::GetOscWeight(), Oscillate(), NueSensitivity::OscillateMatter(), NueFCSensitivity::OscillateMatter(), NuOscProbCalc::ThreeFlavourNuMuToNuElectron(), NuOscProbCalc::ThreeFlavourNuMuToNuMu(), NuOscProbCalc::ThreeFlavourNuMuToNuTau(), and NC::OscProb::ThreeFlavor::TransitionProbability().

00028 {
00029 
00030 
00031   double x = Energy;
00032 
00033   if(nonOscNuFlavor > 0) SetOscParam(OscPar::kNuAntiNu, 1);
00034   else SetOscParam(OscPar::kNuAntiNu, -1);
00035 
00036   double prob = 0.0;
00037   if(abs(nonOscNuFlavor)==14) {
00038     if(abs(nuFlavor)==12)      prob = OscCalc::MuToElec(x);
00039     else if(abs(nuFlavor)==14) prob = OscCalc::MuToMu(x);
00040     else if(abs(nuFlavor)==16) prob = OscCalc::MuToTau(x);
00041   }
00042   if(abs(nonOscNuFlavor)==12) {
00043     if(abs(nuFlavor)==14) prob = OscCalc::ElecToMu(x);
00044     else if(abs(nuFlavor)==12) prob = OscCalc::ElecToElec(x);
00045     else if(abs(nuFlavor)==16) prob = OscCalc::ElecToTau(x);
00046   }
00047   if(abs(nonOscNuFlavor)==16) {
00048     if(abs(nuFlavor)==14) prob = OscCalc::TauToMu(x);
00049     else if(abs(nuFlavor)==12) prob = OscCalc::TauToElec(x);
00050     else if(abs(nuFlavor)==16) prob = OscCalc::TauToTau(x);
00051   } 
00052 
00053   if(prob < -1e-4 && Energy > 0.05){  //I'm sorry but the oscillations are just noisy at low E 
00054       MSG("OscCalc",Msg::kInfo)<<"P(<"<<nonOscNuFlavor<<"->"<<nuFlavor
00055                                <<") less than zero at E = "<<x<<" prob = "<<prob<<endl;
00056       prob = 0;
00057   } 
00058 
00059   return prob;
00060 }

double OscCalc::OscillateNSI ( int  nuFlavor,
int  nonOscNuFlavor,
double  E 
)

Definition at line 674 of file OscCalc.cxx.

References fOscPar, OscPar::kDensity, OscPar::kL, OscPar::kNuAntiNu, myPMNS_NSI, OscProb::PMNS_Fast::P(), OscProb::PMNS_Fast::PropMatter(), OscProb::PMNS_Fast::ResetToFlavour(), UpdateNSI(), and OscPar::z_a.

Referenced by NueStandard::GetOscWeight().

00675 {     
00676   UpdateNSI();
00677 
00678   int NuAntiNu = int(fOscPar[OscPar::kNuAntiNu]);
00679   double L = fOscPar[OscPar::kL];  //baseline
00680   
00681   // electron density in mol/cm^3
00682   double Ne = fOscPar[OscPar::kDensity]*OscPar::z_a;
00683 
00684   // Flavors
00685   int flvi = (abs(nonOscNuFlavor)-12)/2;           // Initial flavour: {0,1,2} for {e,mu,tau}
00686   int flvf = (abs(NuFlavor)-12)/2;                 // Final flavour: {0,1,2} for {e,mu,tau} 
00687 
00688   if(NuFlavor * NuAntiNu < 0){
00689     cout << "Warning: NuFlavor is " << (NuFlavor > 0 ? "positive" : "negative") 
00690          << ", but NuAntiNu is " << (NuAntiNu > 0 ? "positive" : "negative") << endl;
00691     cout << "Setting NuAntiNu = " << - NuAntiNu << " to agree with sign of NuFlavor..." << endl;
00692     NuAntiNu *= -1;
00693   }//Induce Flavor state if NuAntiNu is not set correctly
00694 
00695   // Reset the state to flvi
00696   myPMNS_NSI.ResetToFlavour(flvi);
00697 
00698   // Propagate through matter
00699   myPMNS_NSI.PropMatter(L, E, Ne, NuAntiNu);
00700 
00701   // Get the probability of flvf
00702   return myPMNS_NSI.P(flvf);
00703   
00704 
00705 }

void OscCalc::SetOscParam ( OscPar::OscPar_t  pos,
double  val,
bool  force = false 
)

Definition at line 62 of file OscCalc.cxx.

References OscPar::A_av, MuELoss::e, fCos2Th, fCosDCP, fCosTh, fDmSqrDirty, fElecDensity, fMixDirty, fNSIDirty, fOscPar, fSin2Th, fSinDCP, fSinTh, fV, OscPar::Gf, OscPar::invCmToeV, OscPar::invCmToGeV, OscPar::kDelta, OscPar::kDelta_mutau, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kDensity, OscPar::kEps_ee, OscPar::root2, and OscPar::z_a.

00063 {
00064   if(fabs(fOscPar[pos] - val) > 1e-9 || force){
00065     fOscPar[pos] = val;
00066     if(pos < 3){
00067        fSinTh[pos] = sin(val);
00068        fCosTh[pos] = cos(val);
00069        fSin2Th[pos] = sin(2*val);
00070        fCos2Th[pos] = cos(2*val);
00071        fMixDirty = true;
00072     }
00073     if(pos == OscPar::kDensity){
00074       double ne = OscPar::z_a*OscPar::A_av*val; //electron density #/cm^{3}
00075       fElecDensity = ne*OscPar::invCmToeV*OscPar::invCmToGeV*OscPar::invCmToGeV;
00076       //electron density with units Gev^{2} eV
00077       //Gev^{2} to cancel with GeV^{-2} in Gf
00078                                                                                 
00079        fV = OscPar::root2*OscPar::Gf*fElecDensity; //eV
00080     }
00081     if(pos == OscPar::kDelta){
00082        fCosDCP = cos(val); 
00083        fSinDCP = sin(val);
00084        fMixDirty = true;
00085     }
00086     if(pos == OscPar::kDeltaM23 || pos == OscPar::kDeltaM12){
00087        fDmSqrDirty = true;
00088     }
00089     if(pos >= OscPar::kEps_ee && pos <= OscPar::kDelta_mutau){
00090        fNSIDirty = true;
00091     }
00092   }
00093 }

void OscCalc::SetOscParam ( double *  par,
bool  force = false 
)

Definition at line 95 of file OscCalc.cxx.

References OscPar::kNumParameters.

Referenced by NueSystematic::DoOscCalc(), ANtpTruthInfoBeamAna::GetOscProb(), Extrapolate2D::InitializeOscCalc(), Extrapolate2D::InvertMassHierarchy(), OscCalc(), Oscillate(), NueSensitivity::OscillateMatter(), NueFCSensitivity::OscillateMatter(), Extrapolate2D_Simple::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D_Simple::OscillatePrediction_SepNuNuBar(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), Extrapolate2D::ReadExtrap(), Extrapolate2D_Simple::ReadFNFile(), Extrapolate2D::ReadFNFile(), Extrapolate2D_Simple::ReadFNFileNo3D(), NueSensitivity::RunStandardApproach(), NueFCSensitivity::RunStandardApproach(), NueStandard::SetDefaultOscParam(), Extrapolate2D::SetDeltaCP(), Extrapolate2D::SetNominalOscProb(), NueStandard::SetOscNoMatter(), Extrapolate2D::SetOscPar(), NueStandard::SetOscParam(), NueSensitivity::SetOscParamBase(), NueFCSensitivity::SetOscParamBase(), NNTrain::SetOscParamBase(), NueSystematic::SetOscParams(), Extrapolate2D::SetSinSq2Th13(), NuOscProbCalc::ThreeFlavourNuMuToNuElectron(), NuOscProbCalc::ThreeFlavourNuMuToNuMu(), NuOscProbCalc::ThreeFlavourNuMuToNuTau(), and NC::OscProb::ThreeFlavor::TransitionProbability().

00096 {
00097    for(int i = 0; i < int(OscPar::kNumParameters); i++){
00098      SetOscParam(OscPar::OscPar_t(i), par[i], force);
00099    }
00100 }

double OscCalc::TauToElec ( double  x  ) 

Definition at line 656 of file OscCalc.cxx.

References ElecToTau(), fOscPar, fSinDCP, and OscPar::kDelta.

Referenced by Oscillate().

00657 {
00658   // Flip delta to reverse direction
00659   double oldSinDelta = fSinDCP;
00660   double oldDelta =  fOscPar[OscPar::kDelta];
00661    
00662   fOscPar[OscPar::kDelta] = -oldDelta;
00663   fSinDCP = -oldSinDelta; 
00664                                
00665   double prob = OscCalc::ElecToTau(x);
00666                                  
00667   //restore the world 
00668   fOscPar[OscPar::kDelta] = oldDelta;
00669   fSinDCP = oldSinDelta; 
00670    
00671   return prob; 
00672 }  

double OscCalc::TauToMu ( double  x  ) 

Definition at line 638 of file OscCalc.cxx.

References fOscPar, fSinDCP, OscPar::kDelta, and MuToTau().

Referenced by Oscillate().

00639 {
00640   // Flip delta to reverse direction
00641   double oldSinDelta = fSinDCP;
00642   double oldDelta =  fOscPar[OscPar::kDelta];
00643   
00644   fOscPar[OscPar::kDelta] = -oldDelta;
00645   fSinDCP = -oldSinDelta;
00646                               
00647   double prob = OscCalc::MuToTau(x);
00648                                                                                       
00649   //restore the world
00650   fOscPar[OscPar::kDelta] = oldDelta;
00651   fSinDCP = oldSinDelta;
00652   
00653   return prob;
00654 } 

double OscCalc::TauToTau ( double  x  ) 

Definition at line 617 of file OscCalc.cxx.

References fCosTh, fSin2Th, fSinTh, OscPar::kTh23, and MuToMu().

Referenced by Oscillate().

00618 {
00619 //  TautoTau is the same as Mu->Mu wich sinsq_th23 <-> cossq_th23, sin(2th23) <->-sin(2th23)
00620   double origCos = fCosTh[OscPar::kTh23];
00621   double origSin = fSinTh[OscPar::kTh23];
00622   double orig2Sin = fSin2Th[OscPar::kTh23];
00623 
00624   fCosTh[OscPar::kTh23] = origSin;
00625   fSinTh[OscPar::kTh23] = origCos;
00626   fSin2Th[OscPar::kTh23] = -orig2Sin;
00627 
00628   double prob = OscCalc::MuToMu(x);
00629   
00630   //restore the world
00631   fCosTh[OscPar::kTh23] = origCos;
00632   fSinTh[OscPar::kTh23] = origSin;
00633   fSin2Th[OscPar::kTh23] = orig2Sin;
00634   
00635   return prob;
00636 } 

void OscCalc::UpdateNSI (  )  [private]

Definition at line 102 of file OscCalc.cxx.

References fDmSqrDirty, fMixDirty, fNSIDirty, fOscPar, OscPar::kDelta, OscPar::kDelta_emu, OscPar::kDelta_etau, OscPar::kDelta_mutau, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kEps_ee, OscPar::kEps_emu, OscPar::kEps_etau, OscPar::kEps_mumu, OscPar::kEps_mutau, OscPar::kEps_tautau, OscPar::kTh12, OscPar::kTh13, OscPar::kTh23, myPMNS_NSI, OscProb::PMNS_Fast::SetDeltaMsqrs(), OscProb::PMNS_Fast::SetMix(), and OscProb::PMNS_NSI::SetNSI().

Referenced by OscillateNSI().

00102                        {   
00103 
00104    // Set the mass splittings
00105    if(fDmSqrDirty){
00106       double dmsq_32 = fOscPar[OscPar::kDeltaM23];
00107       double dmsq_21 = fOscPar[OscPar::kDeltaM12];
00108    
00109       myPMNS_NSI.SetDeltaMsqrs(dmsq_21, dmsq_32);
00110    }
00111    
00112    // Set the mixings
00113    if(fMixDirty){
00114       double theta12 = fOscPar[OscPar::kTh12];
00115       double theta13 = fOscPar[OscPar::kTh13];
00116       double theta23 = fOscPar[OscPar::kTh23];
00117    
00118       double Delta_cp = fOscPar[OscPar::kDelta];
00119    
00120       myPMNS_NSI.SetMix(theta12,theta23,theta13,Delta_cp);
00121    }
00122    
00123    // Set the NSI parameters
00124    if(fNSIDirty){
00125       double Eps_ee = fOscPar[OscPar::kEps_ee];
00126       double Eps_emu = fOscPar[OscPar::kEps_emu];
00127       double Eps_etau = fOscPar[OscPar::kEps_etau];
00128       double Eps_mumu = fOscPar[OscPar::kEps_mumu];
00129       double Eps_mutau = fOscPar[OscPar::kEps_mutau];
00130       double Eps_tautau = fOscPar[OscPar::kEps_tautau];
00131    
00132       double Delta_emu = fOscPar[OscPar::kDelta_emu];
00133       double Delta_etau = fOscPar[OscPar::kDelta_etau];
00134       double Delta_mutau = fOscPar[OscPar::kDelta_mutau];
00135    
00136       myPMNS_NSI.SetNSI(Eps_ee,Eps_emu,Eps_etau,Eps_mumu,Eps_mutau,Eps_tautau,Delta_emu,Delta_etau,Delta_mutau);   
00137    }
00138 }


Member Data Documentation

double OscCalc::fCos2Th[3] [private]

Definition at line 106 of file OscCalc.h.

Referenced by ElecToElec(), MuToElec(), MuToTau(), and SetOscParam().

double OscCalc::fCosDCP [private]

Definition at line 110 of file OscCalc.h.

Referenced by MuToElec(), MuToTau(), and SetOscParam().

double OscCalc::fCosTh[3] [private]

Definition at line 104 of file OscCalc.h.

Referenced by ElecToTau(), MuToElec(), MuToTau(), SetOscParam(), and TauToTau().

bool OscCalc::fDmSqrDirty [private]

Definition at line 98 of file OscCalc.h.

Referenced by SetOscParam(), and UpdateNSI().

double OscCalc::fElecDensity [private]

Definition at line 107 of file OscCalc.h.

Referenced by SetOscParam().

bool OscCalc::fMixDirty [private]

Definition at line 97 of file OscCalc.h.

Referenced by SetOscParam(), and UpdateNSI().

bool OscCalc::fNSIDirty [private]

Definition at line 99 of file OscCalc.h.

Referenced by SetOscParam(), and UpdateNSI().

double OscCalc::fOscPar[OscPar::kNumParameters] [private]

Definition at line 101 of file OscCalc.h.

Referenced by ElecToElec(), ElecToMu(), GetOscParam(), MuToElec(), MuToTau(), OscillateNSI(), SetOscParam(), TauToElec(), TauToMu(), and UpdateNSI().

double OscCalc::fSin2Th[3] [private]

Definition at line 105 of file OscCalc.h.

Referenced by ElecToElec(), ElecToTau(), MuToElec(), MuToTau(), SetOscParam(), and TauToTau().

double OscCalc::fSinDCP [private]

Definition at line 109 of file OscCalc.h.

Referenced by ElecToElec(), ElecToMu(), MuToElec(), MuToTau(), SetOscParam(), TauToElec(), and TauToMu().

double OscCalc::fSinTh[3] [private]

Definition at line 103 of file OscCalc.h.

Referenced by ElecToElec(), ElecToTau(), MuToElec(), MuToTau(), SetOscParam(), and TauToTau().

double OscCalc::fV [private]

Definition at line 108 of file OscCalc.h.

Referenced by ElecToElec(), MuToElec(), MuToTau(), and SetOscParam().

OscProb::PMNS_NSI OscCalc::myPMNS_NSI [private]

Definition at line 96 of file OscCalc.h.

Referenced by OscillateNSI(), and UpdateNSI().


The documentation for this class was generated from the following files:
Generated on Mon Aug 11 01:06:39 2014 for loon by  doxygen 1.4.7