PMNS_Avg.cxx

Go to the documentation of this file.

 
//
// Implementation of oscillations of neutrinos in matter in a
// three-neutrino framework with a first order Taylor expansion.
//
// This  class inherits from the PMNS_Fast class
//
 
#include "PMNS_Avg.h"
#include "MatrixDecomp/zheevh3.h"
#include "PremModel.h"
#include <algorithm>
#include <iostream>
 
using namespace OscProb;
 
using namespace std;
 
//.............................................................................
PMNS_Avg::PMNS_Avg() : PMNS_Fast()
{
  fPrem.LoadModel("");
 
  InitializeTaylorsVectors();
 
  SetwidthBin(0, 0);
 
  SetAvgProbPrec(1e-4);
}
 
//.............................................................................
PMNS_Avg::~PMNS_Avg() {}
 
//.............................................................................
void PMNS_Avg::InitializeTaylorsVectors()
{
  fdensityMatrix = matrixC(fNumNus, vectorC(fNumNus, 0));
 
  flambdaInvE = vectorD(fNumNus, 0);
  fVInvE      = matrixC(fNumNus, vectorC(fNumNus, 0));
 
  flambdaCosT = vectorD(fNumNus, 0);
  fVcosT      = matrixC(fNumNus, vectorC(fNumNus, 0));
 
  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;
 
    for (int j = 0; j < fNumNus; j++) {
      fKInvE[i][j] = 0;
      fKcosT[i][j] = 0;
    }
  }
 
  fLayer = fPrem.GetPremLayers().size() - 1;
 
  fdl = -1;
 
  fDetRadius = fPrem.GetDetRadius();
}
 
//.............................................................................
void PMNS_Avg::SetPremModel(OscProb::PremModel& prem) { fPrem = prem; }
 
//.............................................................................
void PMNS_Avg::SetCosT(double cosT) { fcosT = cosT; }
 
//.............................................................................
void PMNS_Avg::SetwidthBin(double dE, double dcosT)
{
  fdInvE = dE;
  fdcosT = dcosT;
}
 
//.............................................................................
void PMNS_Avg::BuildKE(double L)
{
  double lenghtEV = L * kKm2eV;     // L in eV-1
  double bufK     = lenghtEV * 0.5; // L/2 in eV-1
 
  complexD buffer[3];
 
  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]);
    }
  }
}
 
//.............................................................................
void PMNS_Avg::BuildKcosT(double L)
{
  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]); }
    }
  }
}
 
//.............................................................................
double PMNS_Avg::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_Avg::rotateS()
{
  complexD buffer[3];
 
  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_Avg::rotateK()
{
  complexD buffer[3];
 
  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_Avg::MultiplicationRuleS()
{
  complexD save[3];
 
  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_Avg::MultiplicationRuleK(complexD K[3][3])
{
  for (int i = 0; i < fNumNus; i++) {
    complexD buffer[3];
 
    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_Avg::PropagateTaylor()
{
  for (int i = 0; i < int(fNuPaths.size()); i++) {
    PropagatePathTaylor(fNuPaths[i]);
  }
}
 
//.............................................................................
void PMNS_Avg::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(p.length);
 
    // 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_Avg::SolveK(complexD K[3][3], vectorD& lambda, matrixC& V)
{
  double   fEvalGLoBES[3];
  complexD fEvecGLoBES[3][3];
 
  // Solve K for eigensystem using the GLoBES method
  zheevh3(K, fEvecGLoBES, fEvalGLoBES);
 
  // Fill flambdaInvE and fVInvE vectors from GLoBES arrays
  for (int i = 0; i < fNumNus; i++) {
    lambda[i] = fEvalGLoBES[i];
    for (int j = 0; j < fNumNus; j++) { V[i][j] = fEvecGLoBES[i][j]; }
  }
}
 
//.............................................................................
double PMNS_Avg::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[3];
 
  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_Avg::AvgProb(int flvi, int flvf, double E, double dE)
{
  if (E <= 0) return 0;
 
  if (fNuPaths.empty()) return 0;
 
  // Don't average zero width
  if (dE <= 0) return Prob(flvi, flvf, E);
 
  vectorD Ebin = ConvertEtoLoE(E, dE);
 
  // return fct avr proba
  return AvgProbLoE(flvi, flvf, Ebin[0], Ebin[1]);
}
 
//.............................................................................
double PMNS_Avg::AvgProbLoE(int flvi, int flvf, double LoE, double 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 = GetSamplePoints(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++) {
    double w = 1. / pow(samples[j], 2);
 
    avgprob += w * AvgAlgo(flvi, flvf, samples[j], samples[0], L);
 
    sumw += w;
  }
 
  // Return average of probabilities
  return avgprob / sumw;
}
 
//.............................................................................
double PMNS_Avg::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);
}
 
//.............................................................................
double PMNS_Avg::AvgProb(int flvi, int flvf, double E, double cosT,
                         double dcosT)
{
  if (cosT > 0) return 0;
 
  // if (fNuPaths.empty()) return 0; //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 
  // Don't average zero width
  if (dcosT == 0) return Prob(flvi, flvf, E);
 
  vectorD samples = GetSamplePoints(E, cosT, dcosT);
 
  double avgprob = 0;
 
  // Loop over all sample points
  for (int j = 1; j < int(samples.size()); j++) {
    avgprob += AvgAlgoCosT(flvi, flvf, E, samples[j], samples[0]);
  }
 
  // Return average of probabilities
  return avgprob / (samples.size() - 1);
}
 
//.............................................................................
double PMNS_Avg::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);
 
  // return fct avr proba
  return AvgFormula(flvi, flvf, fdcosT, flambdaCosT, fVcosT);
}
 
//.............................................................................
double PMNS_Avg::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);
 
  return AvgProbLoE(flvi, flvf, Ebin[0], Ebin[1], cosT, dcosT);
}
 
//.............................................................................
double PMNS_Avg::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 =
      GetSamplePoints(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++) {
      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]));
 
      sumw += w;
    }
  }
 
  // Return average of probabilities
  return avgprob / ((sumw));
}
 
//.............................................................................
double PMNS_Avg::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_Avg::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]);
}
 
//.............................................................................
void PMNS_Avg::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_Avg::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_Avg::GetSamplePoints(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_Avg::GetSamplePoints(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;
}
 
//.............................................................................
matrixC PMNS_Avg::GetSamplePoints(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;
}
 
//.............................................................................
double PMNS_Avg::AvgFormulaExtrapolation(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[3];
 
  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]) * dbin;
      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_Avg::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_Avg::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_Avg::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);
}