DataUtil/OscCalc.cxx

Go to the documentation of this file.

#include "OscCalc.h"
#include "TMath.h"
#include <cmath>
#include <cstdlib>
#include <iostream>
#include "MessageService/MsgService.h"

using namespace std; 

CVSID("$Id: OscCalc.cxx,v 1.1 2009/04/28 20:13:11 boehm Exp $");

OscCalc::OscCalc()
{
  for(int i = 0; i < 9; i++) fOscPar[i] = 0;  
  double par[9] = {0};
  SetOscParam(par, true);
}

double OscCalc::Oscillate(int nuFlavor, int nonOscNuFlavor, double Energy, double * par)
{
  SetOscParam(par);
  return Oscillate(nuFlavor, nonOscNuFlavor, Energy); 
}


double OscCalc::Oscillate(int nuFlavor, int nonOscNuFlavor,  double Energy)
{
  double x = Energy;

  if(nonOscNuFlavor > 0) SetOscParam(OscPar::kNuAntiNu, 1);
  else SetOscParam(OscPar::kNuAntiNu, -1);

  double prob = 0.0;
  if(abs(nonOscNuFlavor)==14) {
    if(abs(nuFlavor)==12)      prob = OscCalc::MuToElec(x);
    else if(abs(nuFlavor)==14) prob = OscCalc::MuToMu(x);
    else if(abs(nuFlavor)==16) prob = OscCalc::MuToTau(x);
  }
  if(abs(nonOscNuFlavor)==12) {
    if(abs(nuFlavor)==14) prob = OscCalc::ElecToMu(x);
    else if(abs(nuFlavor)==12) prob = OscCalc::ElecToElec(x);
    else if(abs(nuFlavor)==16) prob = OscCalc::ElecToTau(x);
  }
  if(abs(nonOscNuFlavor)==16) {
    if(abs(nuFlavor)==14) prob = OscCalc::TauToMu(x);
    else if(abs(nuFlavor)==12) prob = OscCalc::TauToElec(x);
    else if(abs(nuFlavor)==16) prob = OscCalc::TauToTau(x);
  } 

  if(prob < -1e-4 && Energy > 0.05){  //I'm sorry but the oscillations are just noisy at low E 
      MSG("OscCalc",Msg::kInfo)<<"P(<"<<nonOscNuFlavor<<"->"<<nuFlavor
                               <<") less than zero at E = "<<x<<" prob = "<<prob<<endl;
      prob = 0;
  } 

  return prob;
}

void OscCalc::SetOscParam(OscPar::OscPar_t pos, double val, bool force)
{
  if(fabs(fOscPar[pos] - val) > 1e-9 || force){
    fOscPar[pos] = val;
    if(pos < 3){
       fSinTh[pos] = sin(val);
       fCosTh[pos] = cos(val);
       fSin2Th[pos] = sin(2*val);
       fCos2Th[pos] = cos(2*val);
    }
    if(pos == OscPar::kDensity){
      double ne = OscPar::z_a*OscPar::A_av*val; //electron density #/cm^{3}
      fElecDensity = ne*OscPar::invCmToeV*OscPar::invCmToGeV*OscPar::invCmToGeV;
      //electron density with units Gev^{2} eV
      //Gev^{2} to cancel with GeV^{-2} in Gf
                                                                                
       fV = OscPar::root2*OscPar::Gf*fElecDensity; //eV
    }
    if(pos == OscPar::kDelta){
       fCosDCP = cos(val); 
       fSinDCP = sin(val);
    }
  }
}

void OscCalc::SetOscParam(double *par, bool force)
{
   for(int i = 0; i < int(OscPar::kUnknown); i++){
     SetOscParam(OscPar::OscPar_t(i), par[i], force);
   }
}

void OscCalc::GetOscParam(double *par)
{
   for(int i = 0; i < int(OscPar::kUnknown); i++){
     par[i] = fOscPar[i];
   }
}

double OscCalc::GetOscParam(OscPar::OscPar_t pos)
{
   if(pos >=0 && pos < OscPar::kUnknown)
     return fOscPar[pos];
   else 
     return -9999;
}


double OscCalc::MuToElec(double x)
{
  double E = x; //energy

  //Building the standard terms
  double L = fOscPar[OscPar::kL];  //baseline
  double dmsq_23 = fOscPar[OscPar::kDeltaM23];
  double dmsq_12 = fOscPar[OscPar::kDeltaM12];
  double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}

  double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
  double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
                                                                                
  double cos_th23 = fCosTh[OscPar::kTh23];
  double cos_th12 = fCosTh[OscPar::kTh12];
  double sin_th13 = fSinTh[OscPar::kTh13];
  double cos_th13 = fCosTh[OscPar::kTh13];  
                                                                              
  double sin_2th23 = fSin2Th[OscPar::kTh23];
  double sin_2th12 = fSin2Th[OscPar::kTh12];
  double cos_2th13 = fCos2Th[OscPar::kTh13];
  double cos_2th12 = fCos2Th[OscPar::kTh12];

  double sinsq_th23 = fSinTh[OscPar::kTh23]*fSinTh[OscPar::kTh23];
  double sinsq_th12 = fSinTh[OscPar::kTh12]*fSinTh[OscPar::kTh12];

  double d_cp = fOscPar[OscPar::kDelta];
  double cos_dcp = fCosDCP;
  double sin_dcp = fSinDCP;

  //Building the more complicated terms
  double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
  double A = 2*fV*E*1e9/(dmsq_13);
  double alpha = dmsq_12/dmsq_13;

  // A and d_cp both change sign for antineutrinos
  double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
  A *= plusminus;
  d_cp *= plusminus;
  sin_dcp *= plusminus;

  //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)

  double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));

  double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
  if(fabs(alpha) > 1e-10){  //want to be careful here
    double temp = cos_2th12 - A/alpha;
    C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
  }

  //More complicated sin and cosine terms
  double cosC13Delta = cos(C13*Delta);
  double sinC13Delta = sin(C13*Delta);
  double sin1pADelta = sin((A+1)*Delta);
  double cos1pADelta = cos((A+1)*Delta);
  double sinADelta = sin(A*Delta);
  double sinAm1Delta = sin((A-1)*Delta);
  double cosdpDelta = cos(d_cp+Delta);
  double sinApam2Delta = sin((A+alpha-2)*Delta);
  double cosApam2Delta = cos((A+alpha-2)*Delta);

  double cosaC12Delta = 0;
  double sinaC12Delta = 0; 
  
  if(fabs(alpha) > 1e-10){
    cosaC12Delta = cos(alpha*C12*Delta);
    sinaC12Delta = sin(alpha*C12*Delta);
  } // otherwise not relevant

  //First we calculate the terms for the alpha expansion (good to all orders in th13)
  // this is the equivalent of Eq 47 & 48 corrected for Mu to E instead of E to Mu

  // Leading order term 
  double p1 = sinsq_th23*sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);

  // terms that appear at order alpha
  //first work out the vacuum case since we get 0/0 otherwise.......
  double p2Inner = Delta*cosC13Delta;

  if(fabs(A) > 1e-9)
    p2Inner = Delta*cosC13Delta*(1-A*cos_2th13)/C13 -
                      A*sinC13Delta*(cos_2th13-A)/(C13*C13);

  double p2 = -2*sinsq_th12*sinsq_th23*sinsq_2th13*sinC13Delta/(C13*C13)*p2Inner*alpha;


  //again working out vacuum first.....
  double p3Inner = Delta* cos_th13* cos_th13*(-2*sin_dcp*sinC13Delta*sinC13Delta+2*cos_dcp*sinC13Delta*cosC13Delta);

  if(fabs(A) > 1e-9)
    p3Inner = (sinC13Delta/(A*C13*C13))*(- sin_dcp*(cosC13Delta - cos1pADelta)*C13
                         + cos_dcp*(C13*sin1pADelta - (1-A*cos_2th13)*sinC13Delta));

  double p3 = sin_2th12*sin_2th23*sin_th13*p3Inner*alpha;

  //  p1 + p2 + p3 is the complete contribution for this expansion
  
  // Now for the expansion in orders of sin(th13) (good to all order alpha) 
  //  this is the equivalent of Eq 65 and 66

  // leading order term
  double pa1 = 0.0, pa2 = 0.0;

  // no problems here when A -> 0
  if(fabs(alpha) > 1e-10){
    // leading order term
    pa1 = cos_th23*cos_th23*sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12);

    // and now to calculate the first order in s13 term
    double t1 = (cos_2th12 - A/alpha)/C12 
                  - alpha*A*C12*sinsq_2th12/(2*(1-alpha)*C12*C12);
    double t2 = -cos_dcp*(sinApam2Delta-sinaC12Delta*t1);
  
    double t3 = -(cosaC12Delta-cosApam2Delta)*sin_dcp;
 
    double denom = (1-A-alpha+A*alpha*cos_th12*cos_th12)*C12;
    double t4 = sin_2th12*sin_2th23*(1-alpha)*sinaC12Delta/denom;

    pa2 = t4*(t3+t2)*sin_th13;
  }
  //pa1+pa2 is the complete contribution from this expansion

  // In order to combine the information correctly we need to add the two
  //  expansions and subtract off the terms that are in both (alpha^1, s13^1) 
  //  these may be taken from the expansion to second order in both parameters
  //  Equation 31 

  double t1 = Delta*sinC13Delta*cosdpDelta;
  if(fabs(A) > 1e-9) t1 = sinADelta*cosdpDelta*sinAm1Delta/(A*(A-1));

  double repeated = 2*alpha*sin_2th12*sin_2th23*sin_th13*t1;

  //  Calculate the total probability
  double totalP = p1+p2+p3 + (pa1+pa2) - repeated;

  return totalP;
}

double OscCalc::MuToTau(double x)
{
  double E = x; //energy

  double L = fOscPar[OscPar::kL];  //baseline
  double dmsq_23 = fOscPar[OscPar::kDeltaM23];
  double dmsq_12 = fOscPar[OscPar::kDeltaM12];
  double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}

  double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
  double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
  double sinsq_2th23 = fSin2Th[OscPar::kTh23]*fSin2Th[OscPar::kTh23];
                                                                                
  double cos_th13 = fCosTh[OscPar::kTh13];
  double cos_th12 = fCosTh[OscPar::kTh12];
  double sin_th12 = fSinTh[OscPar::kTh12];
  double sin_th13 = fSinTh[OscPar::kTh13];
                                                                                
  double sin_2th23 = fSin2Th[OscPar::kTh23];
//  double sin_2th13 = fSin2Th[OscPar::kTh13];
  double sin_2th12 = fSin2Th[OscPar::kTh12];
                                                                                
  double cos_2th23 = fCos2Th[OscPar::kTh23];
  double cos_2th13 = fCos2Th[OscPar::kTh13];
  double cos_2th12 = fCos2Th[OscPar::kTh12];

  double d_cp = fOscPar[OscPar::kDelta];
  double cos_dcp = fCosDCP;
  double sin_dcp = fSinDCP;

  //Building the more complicated terms                                                                              
  double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
  double A = 2*fV*E*1e9/(dmsq_13);
  double alpha = dmsq_12/dmsq_13;
  
  // A and d_cp both change sign for antineutrinos
  double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
  A *= plusminus;
  d_cp *= plusminus;
  sin_dcp *= plusminus;

  //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)
                                                                                                                       
  double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));
                                                                                                                       
  double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
  if(fabs(alpha) > 1e-10){  //want to be careful here
    double temp = cos_2th12 - A/alpha;
    C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
  }
                                                                                                                       
  //More complicated sin and cosine terms
  double cosC13Delta = cos(C13*Delta);
  double sinC13Delta = sin(C13*Delta);
  double sin1pADelta = sin((A+1)*Delta);
  double cos1pADelta = cos((A+1)*Delta);
  double sinADelta = sin((A)*Delta);

  double sin1pAmCDelta = sin(0.5*(A+1-C13)*Delta);
  double sin1pApCDelta = sin(0.5*(A+1+C13)*Delta);

  double cosaC12Delta = 0;
  double sinaC12Delta = 0;
                                                                                                                       
  if(fabs(alpha) > 1e-10){
    cosaC12Delta = cos(alpha*C12*Delta);
    sinaC12Delta = sin(alpha*C12*Delta);
  } // otherwise not relevant

  double sinApam2Delta = sin((A+alpha-2)*Delta);
  double cosApam2Delta = cos((A+alpha-2)*Delta);
  double sinAm1Delta = sin((A-1)*Delta);
  double cosAm1Delta = cos((A-1)*Delta);
  double sinDelta = sin(Delta);
  double sin2Delta = sin(2*Delta);

  double cosaC12pApam2Delta = 0;                                                                                                                       
  if(fabs(alpha) > 1e-10){
    cosaC12pApam2Delta = cos((alpha*C12+A+alpha-2)*Delta);
  }

  //First we calculate the terms for the alpha expansion (good to all orders in th13)
  // this is the equivalent of Eq 49 & 50 corrected for Mu to E instead of E to Mu

  // Leading order term
  double pmt_0 = 0.5*sinsq_2th23;
  pmt_0 *= (1 - (cos_2th13-A)/C13)*sin1pAmCDelta*sin1pAmCDelta 
             +  (1 + (cos_2th13-A)/C13)*sin1pApCDelta*sin1pApCDelta
             - 0.5*sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);

  // terms that appear at order alpha
  double t0, t1, t2, t3;
  t0 = (cos_th12*cos_th12-sin_th12*sin_th12*sin_th13*sin_th13
          *(1+2*sin_th13*sin_th13*A+A*A)/(C13*C13))*cosC13Delta*sin1pADelta*2;
  t1 = 2*(cos_th12*cos_th12*cos_th13*cos_th13-cos_th12*cos_th12*sin_th13*sin_th13
         +sin_th12*sin_th12*sin_th13*sin_th13
         +(sin_th12*sin_th12*sin_th13*sin_th13-cos_th12*cos_th12)*A);
  t1 *= sinC13Delta*cos1pADelta/C13;

  t2 =  sin_th12*sin_th12*sinsq_2th13*sinC13Delta/(C13*C13*C13);
  t2 *= A/Delta*sin1pADelta+A/Delta*(cos_2th13-A)/C13*sinC13Delta
          - (1-A*cos_2th13)*cosC13Delta;

  double pmt_1 = -0.5*sinsq_2th23*Delta*(t0+t1+t2);   

  t0 = t1 = t2 = t3 = 0.0;

  t0 = cosC13Delta-cos1pADelta;
  t1 = 2*cos_th13*cos_th13*sin(d_cp)*sinC13Delta/C13*t0;
  t2 = -cos_2th23*cos_dcp*(1+A)*t0*t0;

  t3  = cos_2th23*cos_dcp*(sin1pADelta+(cos_2th13-A)/C13*sinC13Delta);
  t3 *= (1+2*sin_th13*sin_th13*A + A*A)*sinC13Delta/C13 - (1+A)*sin1pADelta;

//  cout<<t1<<"  "<<t2<<"  "<<t3<<endl;

  if(fabs(A) > 1e-9) 
    pmt_1 = pmt_1 + (t1+t2+t3)*sin_th13*sin_2th12*sin_2th23/(2*A*cos_th13*cos_th13);
  else
    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);

  pmt_1 *= alpha;

  //  pmt_0 + pmt_1 is the complete contribution for this expansion
                                                                                                                       
  // Now for the expansion in orders of sin(th13) (good to all order alpha)
  //  this is the equivalent of Eq 67 and 68
                                                                                                                       
  // leading order term
  double pmt_a0 =  0.5*sinsq_2th23;

  pmt_a0 *= 1 - 0.5*sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12)
              - cosaC12pApam2Delta
              - (1 - (cos_2th12 - A/alpha)/C12)*sinaC12Delta*sinApam2Delta;
            
  double denom = (1-A-alpha+A*alpha*cos_th12*cos_th12)*C12;

  t0 = (cosaC12Delta-cosApam2Delta)*(cosaC12Delta-cosApam2Delta);
           
  t1 = (cos_2th12 - A/alpha)/C12*sinaC12Delta+sinApam2Delta;
            
  t2 = ((cos_2th12 - A/alpha)/C12+2*(1-alpha)/(alpha*A*C12))*sinaC12Delta
             + sinApam2Delta;

  t3 = (alpha*A*C12)/2.0*cos_2th23*cos_dcp*(t0 + t1*t2);
  
  t3 += sin_dcp*(1-alpha)*(cosaC12Delta-cosApam2Delta)*sinaC12Delta;

  double pmt_a1 = sin_th13*sin_2th12*sin_2th23/denom*t3;

  // pmt_a1+pmt_a2 is the complete contribution from this expansion
                                                                                                                       
  // In order to combine the information correctly we need to add the two
  //  expansions and subtract off the terms that are in both (alpha^1, s13^1)
  //  and lower order terms
  //  these may be taken from the expansion to second order in both parameters
  //  Equation 34


  // Now for the term of order alpha * s13 or lower order!
  t0 = t1 = t2 = t3 = 0.0;

  t1 = +sin_dcp*sinDelta*sinADelta*sinAm1Delta/(A*(A-1));
  t2 = -1/(A-1)*cos_dcp*sinDelta*(A*sinDelta-sinADelta*cosAm1Delta/A)*cos_2th23;
  t0 =  2*alpha*sin_2th12*sin_2th23*sin_th13*(t1+t2);

  t1 = sinsq_2th23*sinDelta*sinDelta 
       - alpha*sinsq_2th23*cos_th12*cos_th12*Delta*sin2Delta;

  double repeated = t0+t1;

  //  Calculate the total probability
  double totalP = pmt_0 + pmt_1 + pmt_a0 + pmt_a1 - repeated;

  return totalP;
}


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

double OscCalc::ElecToTau(double x)
{
//  EtoTau is the same as E->Mu wich sinsq_th23 <-> cossq_th23, sin(2th23) <->-sin(2th23)
  double origCos = fCosTh[OscPar::kTh23];
  double origSin = fSinTh[OscPar::kTh23];
  double orig2Sin = fSin2Th[OscPar::kTh23];

  fCosTh[OscPar::kTh23] = origSin;
  fSinTh[OscPar::kTh23] = origCos;
  fSin2Th[OscPar::kTh23] = -orig2Sin;

  double prob = OscCalc::ElecToMu(x);

  //restore the world
  fCosTh[OscPar::kTh23] = origCos;
  fSinTh[OscPar::kTh23] = origSin;
  fSin2Th[OscPar::kTh23] = orig2Sin;

  return prob;

/*  This is an option as well but I think results in more cycles
  double p1 = 1. - OscCalc::ElecToMu(x) - OscCalc::ElecToElec(x);
  if(p1 < 0){ cout<<"Damnation! "<<x[0]<<" "<<p1<<endl; p1 = 0;}
  return p1;
*/ 
}

double OscCalc::ElecToMu(double x)
{
  // Flip delta to reverse direction
  double oldSinDelta = fSinDCP;
  double oldDelta =  fOscPar[OscPar::kDelta];

  fOscPar[OscPar::kDelta] = -oldDelta;
  fSinDCP = -oldSinDelta;
                                                                                                                       
  double prob = OscCalc::MuToElec(x);
                                                                                                                       
  //restore the world
  fOscPar[OscPar::kDelta] = oldDelta;
  fSinDCP = oldSinDelta;
 
  return prob;   
}

double OscCalc::ElecToElec(double x)
{
  double E = x; //energy
                                                                                                                       
  //Building the standard terms
  double L = fOscPar[OscPar::kL];  //baseline
  double dmsq_23 = fOscPar[OscPar::kDeltaM23];
  double dmsq_12 = fOscPar[OscPar::kDeltaM12];
  double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}
                                                                                                                       
  double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
  double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
                                                                                                                       
  double sin_th12 = fSinTh[OscPar::kTh12];
 
//  double cos_2th23 = fCos2Th[OscPar::kTh23];
  double cos_2th13 = fCos2Th[OscPar::kTh13];
  double cos_2th12 = fCos2Th[OscPar::kTh12];
                                                                                                                       
  double d_cp = fOscPar[OscPar::kDelta];
  double sin_dcp = fSinDCP;
                                                                                                                       
  //Building the more complicated terms
  double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
  double A = 2*fV*E*1e9/(dmsq_13);
  double alpha = dmsq_12/dmsq_13;
                                                                                                                       
  // A and d_cp both change sign for antineutrinos
  double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
  A *= plusminus;
  d_cp *= plusminus;
  sin_dcp *= plusminus;
                                                                                                                       
  //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)
  double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));
                                                                                                                       
  double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
  if(fabs(alpha) > 1e-10){  //want to be careful here
    double temp = cos_2th12 - A/alpha;
    C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
  }
                                                                                                                       
  //More complicated sin and cosine terms
  double cosC13Delta = cos(C13*Delta);
  double sinC13Delta = sin(C13*Delta);
                                                                                                                       
  double cosaC12Delta = 0;
  double sinaC12Delta = 0;
                                                                                                                       
  if(fabs(alpha) > 1e-10){
    cosaC12Delta = cos(alpha*C12*Delta);
    sinaC12Delta = sin(alpha*C12*Delta);
  } // otherwise not relevant
                                                                                                                       
  //First we calculate the terms for the alpha expansion (good to all orders in th13)
  // this is the equivalent of Eq 45 & 46 corrected for Mu to E instead of E to Mu
                                                                                                                       
  // Leading order term
  double p1 = 1 - sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);
                                                                                                                       
  // terms that appear at order alpha
  double p2Inner = Delta*cosC13Delta*(1-A*cos_2th13)/C13 -
                      A*sinC13Delta*(cos_2th13-A)/(C13*C13);
                                                                                                                       
  double p2 = +2*sin_th12*sin_th12*sinsq_2th13*sinC13Delta/(C13*C13)*p2Inner*alpha;
  //  p1 + p2 is the complete contribution for this expansion
                                                                                                                       
  // Now for the expansion in orders of sin(th13) (good to all order alpha)
  //  this is the equivalent of Eq 63 and 64
                                                                                                                       
  // leading order term
  double pa1 = 1.0, pa2 = 0.0;
                                                                                                                       
  if(fabs(alpha) > 1e-10){
    // leading order term
    pa1 = 1.0 - sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12);
  }
  //pa1 is the complete contribution from this expansion, there is no order s13^1 term
                                                                                                                       
  // In order to combine the information correctly we need to add the two
  //  expansions and subtract off the terms that are in both (alpha^1, s13^1)
  //  these may be taken from the expansion to second order in both parameters
  //  Equation 30
                                                                                                                       
  double repeated = 1;
                                                                                                                       
  //  Calculate the total probability
  double totalP = p1+p2 + (pa1+pa2) - repeated;
                                                                                                                       
  return totalP;
}

double OscCalc::TauToTau(double x)
{
//  TautoTau is the same as Mu->Mu wich sinsq_th23 <-> cossq_th23, sin(2th23) <->-sin(2th23)
  double origCos = fCosTh[OscPar::kTh23];
  double origSin = fSinTh[OscPar::kTh23];
  double orig2Sin = fSin2Th[OscPar::kTh23];

  fCosTh[OscPar::kTh23] = origSin;
  fSinTh[OscPar::kTh23] = origCos;
  fSin2Th[OscPar::kTh23] = -orig2Sin;

  double prob = OscCalc::MuToMu(x);
  
  //restore the world
  fCosTh[OscPar::kTh23] = origCos;
  fSinTh[OscPar::kTh23] = origSin;
  fSin2Th[OscPar::kTh23] = orig2Sin;
  
  return prob;
} 
  
double OscCalc::TauToMu(double x)
{
  // Flip delta to reverse direction
  double oldSinDelta = fSinDCP;
  double oldDelta =  fOscPar[OscPar::kDelta];
  
  fOscPar[OscPar::kDelta] = -oldDelta;
  fSinDCP = -oldSinDelta;
                              
  double prob = OscCalc::MuToTau(x);
                                                                                      
  //restore the world
  fOscPar[OscPar::kDelta] = oldDelta;
  fSinDCP = oldSinDelta;
  
  return prob;
} 
  
double OscCalc::TauToElec(double x)
{
  // Flip delta to reverse direction
  double oldSinDelta = fSinDCP;
  double oldDelta =  fOscPar[OscPar::kDelta];
   
  fOscPar[OscPar::kDelta] = -oldDelta;
  fSinDCP = -oldSinDelta; 
                               
  double prob = OscCalc::ElecToTau(x);
                                 
  //restore the world 
  fOscPar[OscPar::kDelta] = oldDelta;
  fSinDCP = oldSinDelta; 
   
  return prob; 
}  

---