#include"TpcdEdx.h"
#include<iostream>
#include<cassert>
#include<cmath>
#include<algorithm>

ClassImp(TpcdEdx)

bool dedxSort(PAIR p1,PAIR p2){
  //TODO: improve!
  if(fabs(p1.second)<1.E-10) return false;
  if(fabs(p2.second)<1.E-10) return true;
  return p1.first/p1.second < p2.first/p2.second;
}

TpcdEdx::TpcdEdx()
  :  _mom(TVector3(0.,0.,0.))
{
  ;
}

PAIR 
TpcdEdx::getdEdx(unsigned int i) const {
  if(i < nEntries())
    return _data[i];
  else
    ;  //TODO Exception
}


double TpcdEdx::simpleMean(){
  double sumdE=0.;
  double sumdx=0.;
  for(unsigned int i=0;i<_data.size();++i){
	sumdE+=_data[i].first;
	sumdx+=_data[i].second;
  }
  if(std::fabs(sumdx)>1.E-10) return sumdE/sumdx;
  else return 0;
}

double TpcdEdx::truncMean(double trashLow, double trashHigh){
  int N = _data.size();

  //std::cout<<"N="<<N<<std::endl;
  
  sort(_data.begin(),_data.end(),dedxSort);

  double cutLow = trashLow;
  double cutHigh = trashHigh;  //throw away this percentage of entries (high tail)

  int intCutLow = (int)(cutLow*N);
  int intCutHigh = (int)(cutHigh*N);

  double sum=0.;
  int NinSum=0;
  for(int i=intCutLow;i<(N-intCutHigh);++i){
	if(fabs(_data[i].second)>1.E-10){
	  sum+=_data[i].first/_data[i].second;
	  ++NinSum;
	}
  }
  //std::cout << "NinSum " << NinSum  << std::endl;
  if(NinSum>0)return sum/NinSum;
  return -1.;
}

std::vector<double>
TpcdEdx::getESpectrum(bool sorted) const{
  std::vector<double> ret;
  ret.reserve(_data.size());
  for(unsigned int i=0; i<_data.size(); i++)
    ret.push_back(_data[i].first);
  if(sorted)
    std::sort(ret.begin(), ret.end());
  return ret;
}