#include <RAT/MedianCut.hh>
#include <RAT/FitterPMT.hh>
#include <RAT/Processor.hh>
using namespace RAT;
using namespace RAT::PMTSelectors;

#include <vector>
#include <string>
#include <algorithm> // std::sort

using std::string;
using std::vector;
using std::sort;

void MedianCut::Initialise( const std::string& )
{
    // Default values
    fLowCut = -50.0;
    fHighCut = 50.0;
}

void MedianCut::SetD( const std::string& param, double value )
{
    if( param == string( "lowLimit" ) )
        fLowCut = value;
    else if( param == string( "highLimit" ) )
        fHighCut = value;
    else
        throw Processor::ParamUnknown( param );
}

vector<FitterPMT> MedianCut::GetSelectedPMTs( const std::vector<FitterPMT>& data,
                                              const DS::FitVertex& )
{

    // First find the median
    double medianTime = GetMedianTime( data );

    // Now select the PMTs
    vector<FitterPMT> selectedPMTs;

    for( vector<FitterPMT>::const_iterator iPMT = data.begin(); iPMT != data.end(); iPMT++ ){
        const double hitTime = iPMT->GetTime();
        if( hitTime > (medianTime+fLowCut) && hitTime < (medianTime+fHighCut) ){
            selectedPMTs.push_back( *iPMT );
        }
    }
    return selectedPMTs;

}

double MedianCut::GetMedianTime( const std::vector<FitterPMT>& data ) const
{

    vector<double> times(data.size(), 0);
    vector<FitterPMT>::const_iterator iPMT;
    int i;
    for( iPMT=data.begin(), i=0; iPMT != data.end(); iPMT++, i++ )
    {
        times[i] = iPMT->GetTime();
    }

    sort( times.begin(), times.end() );

    if( data.size()%2==0 )
        return ( times[data.size()/2] + times[data.size()/2 + 1] ) / 2;
    return times[data.size()/2];

}