PMNS_Maltoni.cxx

Go to the documentation of this file.

 
//
// Base class for faster oscillations of neutrinos in matter in a
// n-neutrino framework using a Taylor Expansion.
//
//.............................................................................
//
// mloup@apc.in2p3.fr rogan.clark@uclouvain.be
//
 
#include "PMNS_Maltoni.h"
#include <Eigen/Eigenvalues>
 
using namespace std;
 
using namespace OscProb;
 
//.............................................................................
PMNS_Maltoni::PMNS_Maltoni(int numNus) : PMNS_Base(numNus)
{
  fPrem.LoadModel("");
  InitializeTaylorsVectors();
  SetwidthBin(0, 0);
  SetAvgProbPrec(1e-4);
  SetIsOscProbAvg(false);
}
 
//.............................................................................
void PMNS_Maltoni::InitializeTaylorsVectors()
{
  fdensityMatrix = matrixC(fNumNus, vectorC(fNumNus, 0));
 
  flambdaInvE = vectorD(fNumNus, 0);
  fVInvE      = matrixC(fNumNus, vectorC(fNumNus, 0));
  fKInvE      = Eigen::MatrixXcd::Zero(fNumNus, fNumNus);
 
  flambdaCosT = vectorD(fNumNus, 0);
  fVcosT      = matrixC(fNumNus, vectorC(fNumNus, 0));
  fKcosT      = Eigen::MatrixXcd::Zero(fNumNus, fNumNus);
 
  fevolutionMatrixS = matrixC(fNumNus, vectorC(fNumNus, 0));
 
  fSflavor = matrixC(fNumNus, vectorC(fNumNus, 0));
  fKmass   = matrixC(fNumNus, vectorC(fNumNus, 0));
  fKflavor = matrixC(fNumNus, vectorC(fNumNus, 0));
 
  for (int i = 0; i < fNumNus; i++) { fevolutionMatrixS[i][i] = 1; }
 
  fLayer = fPrem.GetPremLayers().size() - 1;
 
  fdl = -1;
 
  fDetRadius = fPrem.GetDetRadius();
}
 
//.............................................................................
void PMNS_Maltoni::SetIsOscProbAvg(bool isOscProbAvg)
{
  fIsOscProbAvg = isOscProbAvg;
}
 
//.............................................................................
void PMNS_Maltoni::SetPremModel(OscProb::PremModel& prem) { fPrem = prem; }
 
//.............................................................................
void PMNS_Maltoni::SetCosT(double cosT) { fcosT = cosT; }
 
//.............................................................................
void PMNS_Maltoni::SetwidthBin(double dE, double dcosT)
{
  fdInvE = dE;
  fdcosT = dcosT;
}
 
//.............................................................................
void PMNS_Maltoni::BuildKcosT()
{
  UpdateHam();
 
  double dL = LnDerivative() * kKm2eV;
 
  for (int j = 0; j < fNumNus; j++) {
    for (int i = 0; i <= j; i++) {
      fKflavor[i][j] = dL * fHam(i, j);
 
      if (i != j) { fKflavor[j][i] = conj(fKflavor[i][j]); }
    }
  }
}
 
//.............................................................................
void PMNS_Maltoni::BuildKE(double L)
{
  double lenghtEV = L * kKm2eV;     // L in eV-1
  double bufK     = lenghtEV * 0.5; // L/2 in eV-1
 
  complexD buffer[fNumNus];
 
  for (int i = 0; i < fNumNus; i++) {
    complexD Hms_kl;
 
    for (int l = 0; l < fNumNus; l++) {
      buffer[l] = 0;
 
      for (int k = 0; k < fNumNus; k++) {
        if (k <= l)
          Hms_kl = fHms[k][l];
        else
          Hms_kl = conj(fHms[l][k]);
 
        if (fIsNuBar && k != l) Hms_kl = conj(Hms_kl);
 
        buffer[l] += conj(fEvec[k][i]) * Hms_kl;
      }
    }
 
    for (int j = 0; j <= i; j++) {
      fKmass[i][j] = 0;
 
      for (int l = 0; l < fNumNus; l++) {
        fKmass[i][j] += buffer[l] * fEvec[l][j];
      }
 
      complexD C;
 
      if (i == j) { C = complexD(1, 0); }
      else {
        double arg = (fEval[i] - fEval[j]) * lenghtEV;
 
        C = (complexD(cos(arg), sin(arg)) - complexD(1, 0.0)) /
            complexD(0.0, arg);
      }
 
      fKmass[i][j] *= bufK * C;
 
      if (i != j) fKmass[j][i] = conj(fKmass[i][j]);
    }
  }
}
 
//.............................................................................
double PMNS_Maltoni::LnDerivative()
{
  double dL = 0;
 
  double L1 = pow(fPrem.GetPremLayers()[fLayer].radius, 2) - fminRsq;
 
  double L2 = -fminRsq;
  if (fLayer > 0) L2 += pow(fPrem.GetPremLayers()[fLayer - 1].radius, 2);
 
  bool ismin = (L2 <= 0 && fcosT < 0);
 
  if (ismin)
    dL = 2 * pow(fDetRadius, 2) * fcosT * pow(L1, -0.5);
  else
    dL = pow(fDetRadius, 2) * fcosT * (pow(L1, -0.5) - pow(L2, -0.5));
 
  if (ismin) fdl = 1;
 
  fLayer += fdl;
 
  return dL;
}
 
//.............................................................................
void PMNS_Maltoni::rotateS()
{
  complexD buffer[fNumNus];
 
  for (int j = 0; j < fNumNus; j++) {
    for (int k = 0; k < fNumNus; k++) {
      buffer[k] = fPhases[k] * conj(fEvec[j][k]);
    }
 
    for (int i = 0; i < fNumNus; i++) {
      fSflavor[i][j] = 0;
      for (int k = 0; k < fNumNus; k++) {
        fSflavor[i][j] += fEvec[i][k] * buffer[k];
      }
    }
  }
}
 
//.............................................................................
void PMNS_Maltoni::rotateK()
{
  complexD buffer[fNumNus];
 
  for (int j = 0; j < fNumNus; j++) {
    for (int k = 0; k < fNumNus; k++) {
      buffer[k] = 0;
 
      for (int l = 0; l < fNumNus; l++) {
        buffer[k] += fKmass[k][l] * conj(fEvec[j][l]);
      }
    }
 
    for (int i = 0; i <= j; i++) {
      fKflavor[i][j] = 0;
 
      for (int k = 0; k < fNumNus; k++) {
        fKflavor[i][j] += fEvec[i][k] * buffer[k];
      }
 
      if (i != j) { fKflavor[j][i] = conj(fKflavor[i][j]); }
    }
  }
}
 
//.............................................................................
void PMNS_Maltoni::MultiplicationRuleS()
{
  complexD save[fNumNus];
 
  for (int j = 0; j < fNumNus; j++) {
    for (int n = 0; n < fNumNus; n++) { save[n] = fevolutionMatrixS[n][j]; }
 
    for (int i = 0; i < fNumNus; i++) {
      fevolutionMatrixS[i][j] = 0;
 
      for (int k = 0; k < fNumNus; k++) {
        fevolutionMatrixS[i][j] += fSflavor[i][k] * save[k];
      }
    }
  }
}
 
//.............................................................................
void PMNS_Maltoni::MultiplicationRuleK(Eigen::MatrixXcd& K)
{
  for (int i = 0; i < fNumNus; i++) {
    complexD buffer[fNumNus];
 
    for (int l = 0; l < fNumNus; l++) {
      for (int k = 0; k < fNumNus; k++) {
        buffer[l] += conj(fevolutionMatrixS[k][i]) * fKflavor[k][l];
      }
    }
 
    for (int j = 0; j <= i; j++) {
      for (int l = 0; l < fNumNus; l++) {
        K(i, j) += buffer[l] * fevolutionMatrixS[l][j];
      }
 
      if (i != j) { K(j, i) = conj(K(i, j)); }
    }
  }
}
 
//.............................................................................
void PMNS_Maltoni::PropagateTaylor()
{
  for (int i = 0; i < int(fNuPaths.size()); i++) {
    PropagatePathTaylor(fNuPaths[i]);
  }
}
 
//.............................................................................
void PMNS_Maltoni::PropagatePathTaylor(NuPath p)
{
  // Set the neutrino path
  SetCurPath(p);
 
  // Solve for eigensystem
  SolveHam();
 
  // Get the evolution matrix in mass basis
  double LengthIneV = kKm2eV * p.length;
  for (int i = 0; i < fNumNus; i++) {
    double arg = fEval[i] * LengthIneV;
    fPhases[i] = complexD(cos(arg), -sin(arg));
  }
 
  // Rotate S in flavor basis
  rotateS();
 
  // if avg on E
  if (fdInvE != 0) {
    // Build KE in mass basis
    BuildKE(p.length);
 
    // Rotate KE in flavor basis
    rotateK();
 
    // Multiply this layer K's with the previous path K's
    MultiplicationRuleK(fKInvE);
  }
  // if avg on cosT
  if (fdcosT != 0) {
    // Build KcosT in mass basis
    BuildKcosT();
 
    // Multiply this layer K's with the previous path K's
    MultiplicationRuleK(fKcosT);
  }
 
  // Multiply this layer S's with the previous path S's
  MultiplicationRuleS();
}
 
//.............................................................................
void PMNS_Maltoni::SolveK(Eigen::MatrixXcd& K, vectorD& lambda, matrixC& V)
{
  if (fNumNus == 4) TemplateSolver<Eigen::Matrix4cd>(K, lambda, V);
  if (fNumNus == 3) TemplateSolver<Eigen::Matrix3cd>(K, lambda, V);
  if (fNumNus == 2)
    TemplateSolver<Eigen::Matrix2cd>(K, lambda, V);
  else
    TemplateSolver<Eigen::MatrixXcd>(K, lambda, V);
}
 
//.............................................................................
template <typename T>
void PMNS_Maltoni::TemplateSolver(Eigen::MatrixXcd& K, vectorD& lambda,
                                  matrixC& V)
{
  Eigen::Ref<T> R(K);
 
  Eigen::SelfAdjointEigenSolver<T> eigensolver(R);
 
  // Fill flambdaInvE and fVInvE vectors from GLoBES arrays
  for (int i = 0; i < fNumNus; i++) {
    lambda[i] = eigensolver.eigenvalues()(i);
    for (int j = 0; j < fNumNus; j++) {
      V[i][j] = eigensolver.eigenvectors()(i, j);
    }
  }
}
 
//.............................................................................
double PMNS_Maltoni::AvgFormula(int flvi, int flvf, double dbin, vectorD lambda,
                                matrixC V)
{
  vectorC SV = vectorC(fNumNus, 0);
 
  for (int i = 0; i < fNumNus; i++) {
    for (int k = 0; k < fNumNus; k++) {
      SV[i] += fevolutionMatrixS[flvf][k] * V[k][i];
    }
  }
 
  complexD buffer[fNumNus];
 
  for (int n = 0; n < fNumNus; n++) { buffer[n] = SV[n] * conj(V[flvi][n]); }
 
  complexD sinc[fNumNus][fNumNus];
  for (int j = 0; j < fNumNus; j++) {
    sinc[j][j] = 1;
    for (int i = 0; i < j; i++) {
      double arg = (lambda[i] - lambda[j]) * dbin / 2;
      sinc[i][j] = sin(arg) / arg;
      sinc[j][i] = sinc[i][j];
    }
  }
 
  complexD P = 0;
 
  for (int j = 0; j < fNumNus; j++) {
    for (int i = 0; i < fNumNus; i++) {
      P += buffer[i] * conj(buffer[j]) * sinc[j][i];
    }
  }
 
  return real(P);
}
 
//.............................................................................
double PMNS_Maltoni::AvgProb(int flvi, int flvf, double E, double dE)
{
  if (fIsOscProbAvg == true) return PMNS_Base::AvgProb(flvi, flvf, E, dE);
 
  // Do nothing if energy not positive
  if (E <= 0) return 0;
 
  if (fNuPaths.empty()) return 0;
 
  // Don't average zero width
  if (dE <= 0) return Prob(flvi, flvf, E);
 
  vectorD LoEbin = ConvertEtoLoE(E, dE);
 
  // Compute average in LoE
  return AvgProbLoE(flvi, flvf, LoEbin[0], LoEbin[1]);
}
 
//.............................................................................
double PMNS_Maltoni::AvgProbLoE(int flvi, int flvf, double LoE, double dLoE)
{
  if (fIsOscProbAvg == true)
    return PMNS_Base::AvgProbLoE(flvi, flvf, LoE, dLoE);
 
  if (LoE <= 0) return 0;
 
  if (fNuPaths.empty()) return 0;
 
  // Don't average zero width
  if (dLoE <= 0) return Prob(flvi, flvf, fPath.length / LoE);
 
  // Get sample points for this bin
  vectorD samples = GetSamplePointsAvgClass(LoE, dLoE);
 
  double avgprob = 0;
  double L       = fPath.length;
  double sumw    = 0;
 
  // Loop over all sample points
  for (int j = 1; j < int(samples.size()); j++) {
    // Set (L/E)^-2 weights
    double w = 1. / pow(samples[j], 2);
 
    avgprob += w * AvgAlgo(flvi, flvf, samples[j], samples[0], L);
 
    // Increment sum of weights
    sumw += w;
  }
 
  // Return weighted average of probabilities
  return avgprob / sumw;
}
 
//.............................................................................
double PMNS_Maltoni::AvgAlgo(int flvi, int flvf, double LoE, double dLoE,
                             double L)
{
  // Set width bin as 1/E
  double d1oE = dLoE / L;
 
  // reset K et S et Ve et lambdaE
  InitializeTaylorsVectors();
 
  SetEnergy(L / LoE);
  SetwidthBin(d1oE, 0);
 
  // Propagate -> get S and K matrix (on the whole path)
  PropagateTaylor();
 
  // DiagolK -> get VE and lambdaE
  SolveK(fKInvE, flambdaInvE, fVInvE);
 
  // return fct avr proba
  return AvgFormula(flvi, flvf, d1oE / kGeV2eV, flambdaInvE, fVInvE);
}
 
//.............................................................................
vectorD PMNS_Maltoni::GetSamplePointsAvgClass(double LoE, double dLoE)
{
  // Set a number of sub-divisions to achieve "good" accuracy
  // This needs to be studied better
  int n_div = ceil(3 * pow(dLoE / LoE, 0.8) * pow(LoE, 0.3) /
                   sqrt(fAvgProbPrec / 1e-4));
 
  // A vector to store sample points
  vectorD Samples;
 
  // Define sub-division width
  Samples.push_back(dLoE / n_div);
 
  // Loop over sub-divisions
  for (int k = 0; k < n_div; k++) {
    // Define sub-division center
    double bctr = LoE - dLoE / 2 + (k + 0.5) * dLoE / n_div;
 
    Samples.push_back(bctr);
 
  } // End of loop over sub-divisions
 
  // Return sample points
  return Samples;
}
 
//.............................................................................
vectorD PMNS_Maltoni::AvgProbVector(int flvi, double E, double dE)
{
  if (fIsOscProbAvg == true) return PMNS_Base::AvgProbVector(flvi, E, dE);
 
  vectorD probs(fNumNus, 0);
 
  // Do nothing if energy is not positive
  if (E <= 0) return probs;
 
  if (fNuPaths.empty()) return probs;
 
  // Don't average zero width
  if (dE <= 0) return ProbVector(flvi, E);
 
  vectorD LoEbin = ConvertEtoLoE(E, dE);
 
  // Compute average in LoE
  return AvgProbVectorLoE(flvi, LoEbin[0], LoEbin[1]);
}
 
//.............................................................................
vectorD PMNS_Maltoni::AvgProbVectorLoE(int flvi, double LoE, double dLoE)
{
  if (fIsOscProbAvg == true)
    return PMNS_Base::AvgProbVectorLoE(flvi, LoE, dLoE);
 
  vectorD probs(fNumNus, 0);
 
  if (LoE <= 0) return probs;
 
  if (fNuPaths.empty()) return probs;
 
  double L = fPath.length;
 
  // Don't average zero width
  if (dLoE <= 0) return ProbVector(flvi, L / LoE);
 
  // Get sample points for this bin
  vectorD samples = GetSamplePointsAvgClass(LoE, dLoE);
 
  double sumw = 0;
 
  // Loop over all sample points
  for (int j = 1; j < int(samples.size()); j++) {
    // Set (L/E)^-2 weights
    double w = 1. / pow(samples[j], 2);
 
    for (int flvf = 0; flvf < fNumNus; flvf++) {
      if (flvf == 0)
        probs[flvf] += w * AvgAlgo(flvi, flvf, samples[j], samples[0], L);
      else
        probs[flvf] +=
            w * AvgFormula(flvi, flvf, fdInvE / kGeV2eV, flambdaInvE, fVInvE);
    }
 
    // Increment sum of weights
    sumw += w;
  }
 
  for (int flvf = 0; flvf < fNumNus; flvf++) {
    // Divide by total sampling weight
    probs[flvf] /= sumw;
  }
 
  // Return weighted average of probabilities
  return probs;
}
 
//.............................................................................
matrixD PMNS_Maltoni::AvgProbMatrix(int nflvi, int nflvf, double E, double dE)
{
  if (fIsOscProbAvg == true)
    return PMNS_Base::AvgProbMatrix(nflvi, nflvf, E, dE);
 
  matrixD probs(nflvi, vectorD(nflvf, 0));
 
  // Do nothing if energy is not positive
  if (E <= 0) return probs;
 
  if (fNuPaths.empty()) return probs;
 
  // Don't average zero width
  if (dE <= 0) return ProbMatrix(nflvi, nflvf, E);
 
  vectorD LoEbin = ConvertEtoLoE(E, dE);
 
  // Compute average in LoE
  return AvgProbMatrixLoE(nflvi, nflvf, LoEbin[0], LoEbin[1]);
}
 
//.............................................................................
matrixD PMNS_Maltoni::AvgProbMatrixLoE(int nflvi, int nflvf, double LoE,
                                       double dLoE)
{
  if (fIsOscProbAvg == true)
    return PMNS_Base::AvgProbMatrixLoE(nflvi, nflvf, LoE, dLoE);
 
  matrixD probs(nflvi, vectorD(nflvf, 0));
 
  if (LoE <= 0) return probs;
 
  if (fNuPaths.empty()) return probs;
 
  double L = fPath.length;
 
  // Don't average zero width
  if (dLoE <= 0) return ProbMatrix(nflvi, nflvf, L / LoE);
 
  // Get sample points for this bin
  vectorD samples = GetSamplePointsAvgClass(LoE, dLoE);
 
  double sumw = 0;
 
  // Loop over all sample points
  for (int j = 1; j < int(samples.size()); j++) {
    // Set (L/E)^-2 weights
    double w = 1. / pow(samples[j], 2);
 
    for (int flvi = 0; flvi < nflvi; flvi++) {
      for (int flvf = 0; flvf < nflvf; flvf++) {
        // if (!TryCacheK())
        if (flvi == 0 && flvf == 0)
          probs[flvi][flvf] +=
              w * AvgAlgo(flvi, flvf, samples[j], samples[0], L);
        else
          probs[flvi][flvf] +=
              w * AvgFormula(flvi, flvf, fdInvE / kGeV2eV, flambdaInvE, fVInvE);
      }
    }
 
    // Increment sum of weights
    sumw += w;
  }
 
  for (int flvi = 0; flvi < nflvi; flvi++) {
    for (int flvf = 0; flvf < nflvf; flvf++) {
      // Divide by total sampling weight
      probs[flvi][flvf] /= sumw;
    }
  }
 
  // Return weighted average of probabilities
  return probs;
}
 
//.............................................................................
double PMNS_Maltoni::AvgProb(int flvi, int flvf, double E, double cosT,
                             double dcosT)
{
  // reset K et S et Ve et lambdaE
  InitializeTaylorsVectors();
 
  SetEnergy(E);
  SetCosT(cosT);
  SetwidthBin(0, dcosT);
 
  fPrem.FillPath(cosT);
  SetPath(fPrem.GetNuPath());
 
  fminRsq = pow(fDetRadius * sqrt(1 - cosT * cosT), 2);
 
  // Propagate -> get S and K matrix (on the whole path)
  PropagateTaylor();
 
  // DiagolK -> get VE and lambdaE
  SolveK(fKcosT, flambdaCosT, fVcosT);
 
  // Compute average
  return AvgFormula(flvi, flvf, fdcosT, flambdaCosT, fVcosT);
}
 
//.............................................................................
double PMNS_Maltoni::AvgAlgoCosT(int flvi, int flvf, double E, double cosT,
                                 double dcosT)
{
  // reset K et S et Ve et lambdaE
  InitializeTaylorsVectors();
 
  SetEnergy(E);
  SetCosT(cosT);
  SetwidthBin(0, dcosT);
 
  fPrem.FillPath(cosT);
  SetPath(fPrem.GetNuPath());
 
  fminRsq = pow(fDetRadius * sqrt(1 - cosT * cosT), 2);
 
  // Propagate -> get S and K matrix (on the whole path)
  PropagateTaylor();
 
  // DiagolK -> get VE and lambdaE
  SolveK(fKcosT, flambdaCosT, fVcosT);
  // Compute average
  return AvgFormula(flvi, flvf, fdcosT, flambdaCosT, fVcosT);
}
 
//.............................................................................
double PMNS_Maltoni::AvgProb(int flvi, int flvf, double E, double dE,
                             double cosT, double dcosT)
{
  if (E <= 0) return 0;
  if (cosT > 0) return 0;
 
  if (fNuPaths.empty()) return 0;
  // Don't average zero width
  if (dE <= 0 && dcosT == 0) return Prob(flvi, flvf, E);
  if (dE <= 0) return AvgProb(flvi, flvf, E, cosT, dcosT);
  if (dcosT == 0) return AvgProb(flvi, flvf, E, dE);
 
  vectorD Ebin = ConvertEtoLoE(E, dE);
 
  // Compute average in LoE
  return AvgProbLoE(flvi, flvf, Ebin[0], Ebin[1], cosT, dcosT);
}
 
//.............................................................................
double PMNS_Maltoni::AvgProbLoE(int flvi, int flvf, double LoE, double dLoE,
                                double cosT, double dcosT)
{
  if (LoE <= 0) return 0;
  if (cosT > 0) return 0;
 
  // if (fNuPaths.empty()) return 0;
 
  // Don't average zero width
  if (dLoE <= 0 && dcosT == 0) return Prob(flvi, flvf, fPath.length / LoE);
  if (dLoE <= 0)
    return AvgProb(flvi, flvf, fPath.length / LoE, cosT, dcosT); 
  if (dcosT == 0) return AvgProbLoE(flvi, flvf, LoE, dLoE);
 
  // Make sample with 1oE and not LoE
  matrixC samples = GetSamplePointsAvgClass(LoE / fPath.length,
                                            dLoE / fPath.length, cosT, dcosT);
 
  int rows = samples.size();
  int cols = samples[0].size();
 
  double avgprob = 0;
  double sumw    = 0;
 
  // Loop over all sample points
  for (int k = 1; k < int(rows); k++) {
    for (int l = 1; l < int(cols); l++) {
      // Set (L/E)^-2 weights
      double w = 1. / pow(real(samples[k][l]), 2);
 
      avgprob +=
          w * AvgAlgo(flvi, flvf, real(samples[k][l]), real(samples[0][0]),
                      imag(samples[k][l]), imag(samples[0][0]));
 
      // Increment sum of weights
      sumw += w;
    }
  }
 
  // Return average of probabilities
  return avgprob / ((sumw));
}
 
//.............................................................................
double PMNS_Maltoni::AvgAlgo(int flvi, int flvf, double InvE, double dInvE,
                             double cosT, double dcosT)
{
  fPrem.FillPath(cosT);
  SetPath(fPrem.GetNuPath());
 
  // reset K et S et Ve et lambdaE
  InitializeTaylorsVectors();
 
  SetEnergy(1 / InvE);
  SetCosT(cosT);
  SetwidthBin(dInvE, dcosT);
 
  fminRsq = pow(fDetRadius * sqrt(1 - cosT * cosT), 2);
 
  // Propagate -> get S and K matrix (on the whole path)
  PropagateTaylor();
 
  // DiagolK -> get VE and lambdaE
  SolveK(fKInvE, flambdaInvE, fVInvE);
  SolveK(fKcosT, flambdaCosT, fVcosT);
 
  return AlgorithmDensityMatrix(flvi, flvf);
}
 
//.............................................................................
double PMNS_Maltoni::AlgorithmDensityMatrix(int flvi, int flvf)
{
  fdensityMatrix[flvi][flvi] = 1;
 
  RotateDensityM(true, fVcosT);
  HadamardProduct(flambdaCosT, fdcosT);
  RotateDensityM(false, fVcosT);
 
  RotateDensityM(true, fVInvE);
  HadamardProduct(flambdaInvE, fdInvE / kGeV2eV);
  RotateDensityM(false, fVInvE);
 
  RotateDensityM(false, fevolutionMatrixS);
 
  return real(fdensityMatrix[flvf][flvf]);
}
 
//.............................................................................
matrixC PMNS_Maltoni::GetSamplePointsAvgClass(double InvE, double dInvE,
                                              double cosT, double dcosT)
{
  // Set a number of sub-divisions to achieve "good" accuracy
  // This needs to be studied better
  int n_divCosT = ceil(380 * pow(InvE, 0.5) * pow(abs(dcosT / cosT), 0.8) /
                       sqrt(fAvgProbPrec / 1e-4));
  int n_divE    = ceil(260 * pow(dInvE / InvE, 0.8) * pow(InvE, 0.6) /
                       sqrt(fAvgProbPrec / 1e-4));
 
  // A matrix to store sample points
  matrixC Samples = matrixC(n_divE + 1, vectorC(n_divCosT + 1, 0));
 
  // Define sub-division width
  Samples[0][0] = complexD(dInvE / n_divE, dcosT / n_divCosT);
 
  // Loop over sub-divisions
  for (int k = 0; k < n_divE; k++) {
    // Define sub-division center for energy
    double bctr_InvE = InvE - dInvE / 2 + (k + 0.5) * dInvE / n_divE;
 
    for (int l = 0; l < n_divCosT; l++) {
      // Define sub-division center for angle
      double bctr_CosT = cosT - dcosT / 2 + (l + 0.5) * dcosT / n_divCosT;
 
      Samples[k + 1][l + 1] = complexD(bctr_InvE, bctr_CosT);
    }
 
  } // End of loop over sub-divisions
 
  // Return sample points
  return Samples;
}
 
//.............................................................................
void PMNS_Maltoni::RotateDensityM(bool to_basis, matrixC V)
{
  matrixC Buffer = matrixC(fNumNus, vectorC(fNumNus, 0));
 
  for (int i = 0; i < fNumNus; i++) {
    for (int j = 0; j < fNumNus; j++) {
      for (int k = 0; k < fNumNus; k++) {
        if (to_basis)
          Buffer[i][j] += fdensityMatrix[i][k] * V[k][j];
        else
          Buffer[i][j] += fdensityMatrix[i][k] * conj(V[j][k]);
      }
    }
  }
 
  for (int i = 0; i < fNumNus; i++) {
    for (int j = i; j < fNumNus; j++) {
      fdensityMatrix[i][j] = 0;
      for (int k = 0; k < fNumNus; k++) {
        if (to_basis)
          fdensityMatrix[i][j] += conj(V[k][i]) * Buffer[k][j];
        else
          fdensityMatrix[i][j] += V[i][k] * Buffer[k][j];
      }
      if (j > i) fdensityMatrix[j][i] = conj(fdensityMatrix[i][j]);
    }
  }
}
 
//.............................................................................
void PMNS_Maltoni::HadamardProduct(vectorD lambda, double dbin)
{
  matrixC sinc = matrixC(fNumNus, vectorC(fNumNus, 0));
  for (int j = 0; j < fNumNus; j++) {
    for (int i = 0; i < j; i++) {
      double arg = (lambda[i] - lambda[j]) * dbin;
      sinc[i][j] = sin(arg) / arg;
 
      sinc[j][i] = sinc[i][j];
    }
  }
 
  for (int i = 0; i < fNumNus; i++) { sinc[i][i] = 1; }
 
  for (int j = 0; j < fNumNus; j++) {
    for (int i = 0; i < fNumNus; i++) {
      fdensityMatrix[i][j] = fdensityMatrix[i][j] * sinc[i][j];
    }
  }
}
 
//.............................................................................
vectorD PMNS_Maltoni::GetSamplePointsAvgClass(double E, double cosT,
                                              double dcosT)
{
  // Set a number of sub-divisions to achieve "good" accuracy
  // This needs to be studied better
  int n_div = ceil(35 * pow(E, -0.4) * pow(abs(dcosT / cosT), 0.8) /
                   sqrt(fAvgProbPrec / 1e-4));
 
  // A vector to store sample points
  vectorD Samples;
 
  // Define sub-division width
  Samples.push_back(dcosT / n_div);
 
  // Loop over sub-divisions
  for (int k = 0; k < n_div; k++) {
    // Define sub-division center
    double bctr = cosT - dcosT / 2 + (k + 0.5) * dcosT / n_div;
 
    Samples.push_back(bctr);
 
  } // End of loop over sub-divisions
 
  // Return sample points
  return Samples;
}
 
//.............................................................................
double PMNS_Maltoni::AvgFormulaExtrapolation(int flvi, int flvf, double dE,
                                             vectorD lambda, matrixC V)
{
  vectorC SV = vectorC(fNumNus, 0);
 
  for (int i = 0; i < fNumNus; i++) {
    for (int k = 0; k < fNumNus; k++) {
      SV[i] += fevolutionMatrixS[flvf][k] * V[k][i];
    }
  }
 
  complexD buffer[fNumNus];
  for (int n = 0; n < fNumNus; n++) { buffer[n] = SV[n] * conj(V[flvi][n]); }
 
  complexD expo[fNumNus][fNumNus];
 
  for (int j = 0; j < fNumNus; j++) {
    for (int i = 0; i < fNumNus; i++) {
      double arg = (lambda[j] - lambda[i]) * dE;
      expo[j][i] = exp(complexD(0.0, arg));
    }
  }
 
  complexD P = 0;
 
  for (int j = 0; j < fNumNus; j++) {
    for (int i = 0; i < fNumNus; i++) {
      P += buffer[i] * conj(buffer[j]) * expo[j][i];
    }
  }
 
  return real(P);
}
 
//.............................................................................
double PMNS_Maltoni::ExtrapolationProb(int flvi, int flvf, double E, double dE)
{
  // reset K et S et Ve et lambdaE
  InitializeTaylorsVectors();
 
  SetEnergy(E);
  SetwidthBin(1 / (E + dE) - 1 / E, 0);
 
  // Propagate -> get S and K matrix (on the whole path)
  PropagateTaylor();
 
  // DiagolK -> get VE and lambdaE
  SolveK(fKInvE, flambdaInvE, fVInvE);
 
  return AvgFormulaExtrapolation(flvi, flvf, fdInvE / kGeV2eV, flambdaInvE,
                                 fVInvE);
}
 
//.............................................................................
double PMNS_Maltoni::ExtrapolationProbLoE(int flvi, int flvf, double LoE,
                                          double dLoE)
{
  // reset K et S et Ve et lambdaE
  InitializeTaylorsVectors();
 
  SetCurPath(AvgPath(fNuPaths));
  double L = fPath.length;
 
  SetEnergy(L / LoE);
  SetwidthBin(dLoE / L, 0);
 
  // Propagate -> get S and K matrix (on the whole path)
  PropagateTaylor();
 
  // DiagolK -> get VE and lambdaE
  SolveK(fKInvE, flambdaInvE, fVInvE);
 
  return AvgFormulaExtrapolation(flvi, flvf, dLoE * kGeV2eV / L, flambdaInvE,
                                 fVInvE);
}
 
//.............................................................................
double PMNS_Maltoni::ExtrapolationProbCosT(int flvi, int flvf, double cosT,
                                           double dcosT)
{
  fPrem.FillPath(cosT);
  SetPath(fPrem.GetNuPath());
 
  // reset K et S et Ve et lambdaE
  InitializeTaylorsVectors();
 
  // SetEnergy(E);
  SetCosT(cosT);
  SetwidthBin(0, dcosT);
 
  fminRsq = pow(fDetRadius * sqrt(1 - cosT * cosT), 2);
 
  // Propagate -> get S and K matrix (on the whole path)
  PropagateTaylor();
 
  // DiagolK -> get VE and lambdaE
  SolveK(fKcosT, flambdaCosT, fVcosT);
 
  return AvgFormulaExtrapolation(flvi, flvf, fdcosT, flambdaCosT, fVcosT);
}