PMNS_Deco.cxx

Go to the documentation of this file.

 
//
// Implementation of oscillations of neutrinos in matter in a
// three-neutrino framework with decoherence.
//
// This  class inherits from the PMNS_DensityMatrix class
//
// jcoelho\@apc.in2p3.fr
 
#include <iostream>
 
#include "PMNS_Deco.h"
#include "exceptions.h"
 
using namespace OscProb;
 
using namespace std;
 
//.............................................................................
PMNS_Deco::PMNS_Deco() : PMNS_DensityMatrix(), fGamma()
{
  SetGamma(2, 0);
  SetGamma(3, 0);
  SetDecoAngle(0);
  SetPower(0);
}
 
//.............................................................................
PMNS_Deco::~PMNS_Deco() {}
 
//.............................................................................
void PMNS_Deco::SetGamma(int j, double val)
{
  THROW_ON_INVALID_ARG(j == 2 || j == 3, j);
 
  if (val < 0) {
    cerr << "WARNING: Gamma_" << j << 1 << " must be positive. "
         << "Setting it to absolute value of input: " << abs(val) << endl;
  }
 
  fGamma[j - 1] = abs(val);
}
 
//.............................................................................
void PMNS_Deco::SetGamma32(double gamma32)
{
  if (gamma32 < 0) {
    cerr << "WARNING: Gamma_32 must be positive. "
         << "Setting it to absolute value of input: " << abs(gamma32) << endl;
  }
 
  double min32 = 0.25 * fGamma[1];
 
  if (fGamma[0] >= 0)
    min32 *= 1 - pow(fGamma[0], 2);
  else
    min32 *= 1 + 3 * pow(fGamma[0], 2);
 
  if (gamma32 < min32) {
    if (fGamma[0] >= 0 || 4 * gamma32 / fGamma[1] < 1) {
      fGamma[0] = sqrt(1 - 4 * gamma32 / fGamma[1]);
      fGamma[2] = 0.25 * fGamma[1] * (1 + 3 * pow(fGamma[0], 2));
    }
    else {
      fGamma[0] = -sqrt((4 * gamma32 / fGamma[1] - 1) / 3);
      fGamma[2] = 0.25 * fGamma[1] * (1 - pow(fGamma[0], 2));
    }
 
    cerr << "WARNING: Impossible to have Gamma32 = " << gamma32
         << " with current Gamma21 and theta parameters." << endl
         << "         Changing the value of cos(theta) to " << fGamma[0]
         << endl;
 
    return;
  }
 
  double arg = fGamma[1] * (4 * gamma32 - fGamma[1] * (1 - pow(fGamma[0], 2)));
 
  if (arg < 0) {
    arg       = 0;
    fGamma[0] = sqrt(1 - 4 * gamma32 / fGamma[1]);
    cerr << "WARNING: Imaginary Gamma31 found. "
         << "Changing the value of cos(theta) to " << fGamma[0] << endl;
  }
 
  double gamma31 = gamma32 + fGamma[0] * (fGamma[0] * fGamma[1] + sqrt(arg));
 
  fGamma[2] = gamma31;
 
  // Sanity check
  double get_gamma32 = GetGamma(3, 2);
  THROW_ON_INVALID_ARG(fabs(gamma32 - get_gamma32) < 1e-6, gamma32, gamma32);
}
 
//.............................................................................
void PMNS_Deco::SetDecoAngle(double th) { fGamma[0] = cos(th); }
 
//.............................................................................
void PMNS_Deco::SetPower(double n) { fPower = n; }
 
//.............................................................................
double PMNS_Deco::GetDecoAngle() { return acos(fGamma[0]); }
 
//.............................................................................
double PMNS_Deco::GetPower() { return fPower; }
 
//.............................................................................
double PMNS_Deco::GetGamma(int i, int j)
{
  if (i < j) {
    cerr << "WARNING: First argument should be larger than second argument"
         << endl
         << "Setting reverse order (Gamma_" << j << i << "). " << endl;
    int temp = i;
    i        = j;
    j        = temp;
  }
  THROW_ON_INVALID_ARG(i == 2 || i == 3, i);
  THROW_ON_INVALID_ARG(j == 1 || j == 2, i);
  THROW_ON_INVALID_ARG(i > j, i, j);
 
  if (j == 1) { return fGamma[i - 1]; }
  else {
    // Combine Gamma31, Gamma21 and an angle (fGamma[0] = cos(th)) to make
    // Gamma32
    double arg =
        fGamma[1] * (4 * fGamma[2] - fGamma[1] * (1 - pow(fGamma[0], 2)));
 
    if (arg < 0) return fGamma[1] - 3 * fGamma[2];
 
    return fGamma[2] + fGamma[1] * pow(fGamma[0], 2) - fGamma[0] * sqrt(arg);
  }
}
 
//.............................................................................
vectorI sort3(const vectorD& x)
{
  vectorI out(3, 0);
 
  // 1st element is smallest
  if (x[0] < x[1] && x[0] < x[2]) {
    // 3rd element is largest
    if (x[1] < x[2]) {
      out[1] = 1;
      out[2] = 2;
    }
    // 2nd element is largest
    else {
      out[1] = 2;
      out[2] = 1;
    }
  }
  // 2nd element is smallest
  else if (x[1] < x[2]) {
    out[0] = 1;
    // 3rd element is largest
    if (x[0] < x[2]) out[2] = 2;
    // 1st element is largest
    else
      out[1] = 2;
  }
  // 3rd element is smallest
  else {
    out[0] = 2;
    // 2nd element is largest
    if (x[0] < x[1]) out[2] = 1;
    // 1st element is largest
    else
      out[1] = 1;
  }
 
  return out;
}
 
//.............................................................................
void PMNS_Deco::PropagatePath(NuPath p)
{
  // Set the neutrino path
  SetCurPath(p);
 
  // Solve for eigensystem
  SolveHam();
 
  // Rotate to effective mass basis
  RotateState(true, fEvec);
 
  // Some ugly way of matching gamma and dmsqr indices
  vectorI dm_idx = sort3(fDm);
  vectorI ev_idx = sort3(fEval);
 
  int idx[3];
  for (int i = 0; i < fNumNus; i++) idx[ev_idx[i]] = dm_idx[i];
 
  // Power law dependency of gamma parameters
  double energyCorr = pow(fEnergy, fPower);
 
  // Convert path length to eV
  double lengthIneV = kKm2eV * p.length;
 
  // Apply phase rotation to off-diagonal elements
  for (int j = 0; j < fNumNus; j++) {
    for (int i = 0; i < j; i++) {
      // Get the appropriate gamma parameters based on sorted indices
      double gamma_ij =
          GetGamma(max(idx[i], idx[j]) + 1, min(idx[i], idx[j]) + 1) *
          energyCorr;
      // Multiply by path length
      gamma_ij *= kGeV2eV * lengthIneV;
      // This is the standard oscillation phase
      double arg = (fEval[i] - fEval[j]) * lengthIneV;
      // apply phase to rho
      fRho[i][j] *= exp(-gamma_ij) * complexD(cos(arg), -sin(arg));
      fRho[j][i] = conj(fRho[i][j]);
    }
  }
 
  // Rotate back to flavour basis
  RotateState(false, fEvec);
}