////////////////////////////////////////////////////////////////////////
//
// Implements functions to average neutrino paths
//
// jcoelho\@apc.in2p3.fr
////////////////////////////////////////////////////////////////////////

#include "NuPath.h"

using namespace std;

using namespace OscProb;

//......................................................................
///
/// Get the merged average of two paths
///
/// This method will merge two paths and take their average density
/// weighted by Z/A and path length.
///
/// The Z/A will be the average weighted by path length
///
/// @param p1 - The first path to merge
/// @param p2 - The second path to merge
/// @return The merged path
///
NuPath OscProb::AvgPath(NuPath& p1, NuPath& p2){

  // Start with the first path
  NuPath mergedPath = p1;

  // Add the second length
  mergedPath.length += p2.length;

  // Compute weighted average of Z/A
  mergedPath.zoa = (p1.zoa*p1.length + p2.zoa*p2.length) / (p1.length + p2.length);

  // Compute weighted average of density
  mergedPath.density = (p1.density*p1.zoa*p1.length + p2.density*p2.zoa*p2.length) / (p1.zoa*p1.length + p2.zoa*p2.length);

  // return merged path
  return mergedPath;

}

//......................................................................
///
/// Get the merged average of a vector of paths
///
/// This method will merge a set of paths and take their average density
/// weighted by Z/A and path length.
///
/// The Z/A will be the average weighted by path length
///
/// @param pv - vector of paths to merge
/// @return The merged path
///
NuPath OscProb::AvgPath(std::vector<NuPath>& pv){

  // Get size of vector
  int np = pv.size();
  
  // Start with the first path
  NuPath mergedPath;
  
  // If vector is not empty, start on first path
  if(np>0) mergedPath = pv[0];
  else return mergedPath;

  // Merge each of the following paths
  for(int i=1; i<np; i++){
    mergedPath = AvgPath(mergedPath, pv[i]);
  }

  // return merged path
  return mergedPath;

}

//......................................................................
///
/// Merge two specific paths by their indices in a path vector
///
/// @param inputPath - The original vector of paths to merge
/// @param j,k - The indices of the two paths to merge
/// @return The merged vector of paths
///
vector<NuPath> OscProb::MergePaths(std::vector<NuPath>& inputPath, int j, int k){

  // Output vector
  vector<NuPath> mergedPath;

  // Loop over input paths
  for(int i=0; i<inputPath.size(); i++){

    // If first index, merge the paths j and k
    if(i==j) mergedPath.push_back(AvgPath(inputPath[j], inputPath[k]));
    // If not second index add the path as is
    else if(i!=k) mergedPath.push_back(inputPath[i]);

  }// End of loop

  // return merged vector
  return mergedPath;

}