NCOscProb.cxx

Go to the documentation of this file.

#include "NCUtils/NCOscProb.h"
#include "NCUtils/NCType.h"

// For three flavour formulae
#include "OscProb/OscCalc.h"

#include "MessageService/MsgService.h"
#include "TComplex.h"
#include "TMath.h"

#include <cassert>
#include <cmath>

CVSID("$Id: NCOscProb.cxx,v 1.44 2014/02/18 03:18:27 rhatcher Exp $");

using namespace NCType;
using std::vector;
using NC::Utility::SQR;

//-------------------------------------------------------
bool NC::OscProb::OscPars::IsEquiv(const NC::OscProb::OscPars* rhs) const
{
  if(!rhs) return false;
  if(rhs->OscillationModel() != OscillationModel()) return false;
  for(int n = 0; n < NCType::kNumParameters; ++n){
    if(rhs->fVals[n].Uninitialized() != fVals[n].Uninitialized())
      return false;
    if(!fVals[n].Uninitialized() && rhs->fVals[n] != fVals[n])
      return false;
  }
  return true;
}

//-------------------------------------------------------
std::ostream& NC::OscProb::OscPars::Print(std::ostream& os) const
{
  os << "[ ";
  for(int n = 0; n < NCType::kNumParameters; ++n)
    if(!fVals[n].Uninitialized()) os << fVals[n] << " ";
  os << "]";
  return os;
}

//-------------------------------------------------------
bool NC::OscProb::OscPars::operator<(const NC::OscProb::OscPars& rhs) const
{
  if(IsEquiv(&rhs)) return false;

  if(fModel < rhs.fModel) return true;
  if(fModel > rhs.fModel) return false;

  for(int n = 0; n < NCType::kNumParameters; ++n){
    if(fVals[n].Uninitialized()) continue;
    assert(!rhs.fVals[n].Uninitialized());
    if(fVals[n] < rhs.fVals[n]) return true;
    if(fVals[n] > rhs.fVals[n]) return false;
  }

  assert(0 && "Not reached");
}

//-------------------------------------------------------
std::ostream& NC::OscProb::operator<<(std::ostream& os,
                                      const NC::OscProb::OscPars& rhs)
{
  return rhs.Print(os);
}

//-------------------------------------------------------
double NC::OscProb::OscProbFs(double ssDeltaMSqr31,
                              double amu,
                              double fs,
                              double ae,
                              int oscMode)
{
  switch(oscMode){
  case NCType::kNuMuToNuMu: return 1. - amu*ssDeltaMSqr31;
  case NCType::kNuMuToNuS: return fs*amu*ssDeltaMSqr31;
  case NCType::kNuMuToNuE: return ae*ssDeltaMSqr31;
  case NCType::kNuMuToNuTau: return (amu*(1.-fs)-ae)*ssDeltaMSqr31;
  case NCType::kNuEToNuE: return 1;
  default:
    assert(0 && "Bad oscMode");
  }
}

//----------------------------------------------------------------------
double NC::OscProb::ThreeFlavor::
TransitionProbability(NCType::EOscMode oscMode,
                      NCType::EEventType interaction,
                      double baseline,
                      double energy) const
{
  if(interaction == NCType::kNC){
    switch(oscMode){
      case NCType::kNuMuToNuTau: return 0.;
      case NCType::kNuMuToNuE:   return 0.;
      case NCType::kNuMuToNuMu:  return 1.;
      case NCType::kNuMuToNuS:
        assert(0 && "TODO - Old code would return 1 here, that's bad");
      case NCType::kNuEToNuE:
        return 1;
    }
  }

  double pars[9];
  pars[OscPar::kTh12] = Theta12();
  pars[OscPar::kTh23] = Theta23();
  pars[OscPar::kTh13] = Theta13();
  pars[OscPar::kDeltaM23] = DeltaMSqr32()*Hierarchy();
  pars[OscPar::kDeltaM12] = DeltaMSqr12();
  pars[OscPar::kDelta] = Delta13();
  pars[OscPar::kDensity] = RockDensity();
  pars[OscPar::kL] = baseline;
  pars[OscPar::kNuAntiNu] = -1; // Treats everything as matter

  // static means it should hold on to derived values from variables that
  // haven't changed.
  static OscCalc oc;
  oc.SetOscParam(pars);

  switch(oscMode){
    case NCType::kNuMuToNuMu:  return oc.Oscillate(14, 14, energy);
    case NCType::kNuMuToNuTau: return oc.Oscillate(16, 14, energy);
    case NCType::kNuMuToNuE:   return oc.Oscillate(12, 14, energy);
    case NCType::kNuEToNuE:    return oc.Oscillate(12, 12, energy);
    case NCType::kNuMuToNuS:
      assert(0 && "TODO - Old code would return 1 here - bad");
  }

  assert(0 && "Not reached");
}

//----------------------------------------------------------------------
double NC::OscProb::SterileFraction::TransitionProbability(
    NCType::EOscMode oscMode,
    NCType::EEventType interaction,
    double baseline,
    double energy) const
{
  double amu = 4.*UMu3Sqr()*(1.-UMu3Sqr());
  double ssDeltaMSqr32 = FindSinSqrDeltaMSqr(energy, baseline, DeltaMSqr32());
  double fs = Fs();
  double ae = 4.*UMu3Sqr()*UE3Sqr();

  if(interaction == NCType::kNC){
    switch(oscMode){
      case NCType::kNuMuToNuTau: return 0.;
      case NCType::kNuMuToNuE:   return 0.;
      case NCType::kNuEToNuE:    return 1.;
      case NCType::kNuMuToNuS:
        return OscProbFs(ssDeltaMSqr32, amu, fs, ae,
                         NCType::kNuMuToNuS);
      case NCType::kNuMuToNuMu:
        return 1.-OscProbFs(ssDeltaMSqr32, amu, fs, ae,
                            NCType::kNuMuToNuS);
    }
  }

  return OscProbFs(ssDeltaMSqr32, amu, fs, ae, oscMode);
}

//----------------------------------------------------------------------
double NC::OscProb::SterileFractionTauNorm::TransitionProbability(
    NCType::EOscMode oscMode,
    NCType::EEventType interaction,
    double baseline,
    double energy) const
{
  double amu = 4.*UMu3Sqr()*(1.-UMu3Sqr());
  double ssDeltaMSqr32 = FindSinSqrDeltaMSqr(energy, baseline, DeltaMSqr32());
  double fs = Fs();
  double ae = 4.*UMu3Sqr()*UE3Sqr();

  if(interaction == NCType::kNC){
    switch(oscMode){
      case NCType::kNuMuToNuTau: return 0.;
      case NCType::kNuMuToNuE:   return 0.;
      case NCType::kNuEToNuE:    return 1.;
      case NCType::kNuMuToNuS:
        return OscProbFs(ssDeltaMSqr32, amu, fs, ae,
                         NCType::kNuMuToNuS);
      case NCType::kNuMuToNuMu:
        return 1.-OscProbFs(ssDeltaMSqr32, amu, fs, ae,
                            NCType::kNuMuToNuS);
    }
  }

  double ret=OscProbFs(ssDeltaMSqr32, amu, fs, ae, oscMode);
  if (oscMode==NCType::kNuMuToNuTau) ret*=TauScale();

  if (ret<0) return 0;
  return ret;
}

//----------------------------------------------------------------------
double NC::OscProb::FindSinSqrDeltaMSqr(double E, double L, double dmsq)
{
  return SQR(TMath::Sin(NCType::k127*L/E*dmsq));
}

//----------------------------------------------------------------------x
NCType::EOscMode NC::OscProb::ToOscMode(int from, int to)
{
  // Use the same oscillation designators for antiparticle oscillations
  if(from < 0) from *= -1;
  if(to < 0) to *= -1;

  // Translate e to nue, mu to numu, tau to nutau
  if(from%2) ++from;
  if(to%2) ++to;

  if(from == 14 && to == 14) return NCType::kNuMuToNuMu;
  if(from == 14 && to == 16) return NCType::kNuMuToNuTau;
  if(from == 14 && to == 12) return NCType::kNuMuToNuE;
  if(from == 14 && to ==  0) return NCType::kNuMuToNuS;
  if(from == 12 && to == 12) return NCType::kNuEToNuE;

  assert(0 && "Don't have an enum for this oscillation mode");
}

//----------------------------------------------------------------------
double NC::OscProb::FourFlavorBase::
TransitionProbability(NCType::EOscMode mode,
                      NCType::EEventType intType,
                      double baseline,
                      double trueE) const
{
  if(intType == NCType::kNC){
    switch(mode){
    case NCType::kNuMuToNuMu:
      return 1-TransitionProbability(NCType::kNuMuToNuS, intType, baseline, trueE);
    case NCType::kNuMuToNuTau:
    case NCType::kNuMuToNuE:
      return 0;
    case NCType::kNuEToNuE:
      return 1;
    case NCType::kNuMuToNuS:
      // Fall through to actually calculate the probability
      break;
    default:
      assert(0 && "Unknown mode");
    }
  }

  fLoverE = baseline/trueE;

  // Precalculate sines and cosines
  const double s13 = TMath::Sin(Theta13());
  const double s23 = TMath::Sin(Theta23());
  const double s14 = TMath::Sin(Theta14());
  const double s24 = TMath::Sin(Theta24());
  const double s34 = TMath::Sin(Theta34());

  const double c13 = TMath::Cos(Theta13());
  const double c23 = TMath::Cos(Theta23());
  const double c14 = TMath::Cos(Theta14());
  const double c24 = TMath::Cos(Theta24());
  const double c34 = TMath::Cos(Theta34());

  // Precalculate complex phases
  const TComplex emd1(1, -Delta1(), true); // 'true' selects polar coordinates
  const TComplex epd1 = TComplex::Conjugate(emd1);
  const TComplex emd2(1, -Delta2(), true);
  const TComplex epd2 = TComplex::Conjugate(emd2);

  // Matrix elements
  const TComplex ue3 = emd1*c14*s13;

  const TComplex ue4 = s14;

  const TComplex umu3 = emd1*emd2*-s13*s14*s24 + c13*s23*c24;

  const TComplex umu4 = emd2*c14*s24;

  const TComplex utau3_1 = emd1*-s13*s14*c24*s34;
  const TComplex utau3_2 = epd2*-c13*s23*s24*s34;
  const TComplex utau3 = c13*c23*c34 + utau3_1 + utau3_2;

  const TComplex utau4 = c14*c24*s34;

  const TComplex us3_1 = emd1*-s13*s14*c24*c34;
  const TComplex us3_2 = epd2*-c13*s23*s24*c34;
  const TComplex us3 = -c13*c23*s34 + us3_1 + us3_2;

  const TComplex us4 = c14*c24*c34;

  // Find the square of the moduli for the elements
  const double ue3Sqr   = ue3.Rho2();
  const double ue4Sqr   = ue4.Rho2();
  const double umu3Sqr  = umu3.Rho2();
  const double umu4Sqr  = umu4.Rho2();
  const double utau3Sqr = utau3.Rho2();
  const double utau4Sqr = utau4.Rho2();
  const double us3Sqr   = us3.Rho2();
  const double us4Sqr   = us4.Rho2();

  // Check unitarity
  const double epsilon = 1e-5;
  const double col3Sum  = ue3Sqr+umu3Sqr+utau3Sqr+us3Sqr;
  const double col4Sum = ue4Sqr+umu4Sqr+utau4Sqr+us4Sqr;
  if(TMath::Abs(col3Sum-1) > epsilon ||
     TMath::Abs(col4Sum-1) > epsilon)
    MSG("NCOscProb", Msg::kWarning) << "Warning: unitarity violated!!!!"
                                    << endl;

  const TComplex umu4conj = TComplex::Conjugate(umu4);

  // Calculate some cross-terms
  const TComplex umuue   = umu4conj*ue4  *umu3*TComplex::Conjugate(ue3);
  const TComplex umuutau = umu4conj*utau4*umu3*TComplex::Conjugate(utau3);
  const TComplex umuus   = umu4conj*us4  *umu3*TComplex::Conjugate(us3);

  double p  = 0;

#ifdef CASEMUTOALPHA
#error Someone already defined CASEMUTOALPHA, have to choose a different name
#endif


#define CASEMUTOALPHA(alpha) \
    p  = umu3Sqr*u##alpha##3Sqr*sinSqrDelta31(); \
    p += umu4Sqr*u##alpha##4Sqr*sinSqrDelta41(); \
    p += (umuu##alpha).Re()*(sinSqrDelta31() - sinSqrDelta43() + sinSqrDelta41()); \
    p *= 4; \
    p += 2.*(umuu##alpha).Im()*(sin2Delta31() - sin2Delta41() + sin2Delta43()); \
    return p

  switch(mode){
  case NCType::kNuMuToNuMu:
    p  = umu3Sqr*(1-umu3Sqr-umu4Sqr)*sinSqrDelta31();
    p += umu4Sqr*(1-umu3Sqr-umu4Sqr)*sinSqrDelta41();
    p += umu4Sqr*umu3Sqr*sinSqrDelta43();
    p *= 4.;
    return 1 - p;

  case NCType::kNuMuToNuTau:
    CASEMUTOALPHA(tau);

  case NCType::kNuMuToNuE:
    CASEMUTOALPHA(e);

  case NCType::kNuMuToNuS:
    CASEMUTOALPHA(s);

  case NCType::kNuEToNuE:
    // OK since there are very few beam nues anyway?
    return 1;

  default:
    assert(0 && "Unknown mode");
  }

#undef CASEMUTOALPHA
}


double NC::OscProb::Decay::TransitionProbability(EOscMode mode,
                                                 EEventType intType,
                                                 double L,
                                                 double E) const
{
  const double ssq = SQR(TMath::Sin(Theta()));
  const double csq = SQR(TMath::Cos(Theta()));
  const double arg = L/(2*E);
  const double exp = TMath::Exp(-Alpha()*arg);

  // Convert to correct units
  const double dmsq = 4*NCType::k127*DeltaMSqr();

  switch(mode){
  case kNuMuToNuMu:
    if(intType == NCType::kNC) return 1-(1-SQR(exp))*ssq;

    return (SQR(csq) + SQR(ssq*exp) +
            2*csq*ssq*exp*TMath::Cos(dmsq*arg));

  case kNuMuToNuTau:
    if(intType == NCType::kNC) return 0;

    return csq*ssq*(1+SQR(exp)-2*TMath::Cos(dmsq*arg)*exp);

  case kNuMuToNuE:
    return 0;

  // Call the thing nu2 decays to "sterile"
  case kNuMuToNuS:
    return (1-SQR(exp))*ssq;

  case kNuEToNuE:
    return 1;
  }

  assert(0 && "Unknown oscillation mode");
}


double NC::OscProb::Decoherence::TransitionProbability(EOscMode mode,
                                                       EEventType intType,
                                                       double baseline,
                                                       double trueEnergy) const
{
  switch(mode){
  case kNuMuToNuMu:
    if(intType == NCType::kNC) return 1;
    return 1-SQR(TMath::Sin(2*Theta()))/2*(1-TMath::Exp(-SQR(Mu())*baseline/(2*trueEnergy)));

  case kNuMuToNuTau:
    return 0;

  case kNuMuToNuE:
    return 0;

  case kNuMuToNuS:
    return 0;

  case kNuEToNuE:
    return 1;
  }

  assert(0 && "Unknown oscillation mode");
}

---