#include <RAT/Ext0NuTimeRes.hh>
#include <RAT/DB.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DU/LightPathCalculator.hh>
#include <RAT/DU/Utility.hh>

#include <math.h>

using namespace RAT;
using namespace RAT::Classifiers;
using namespace RAT::DS;
using namespace RAT::DU;
using namespace std;


std::string Ext0NuTimeRes::GetName() const {
    return Name();
}

void Ext0NuTimeRes::Initialise(const std::string& param) {
    fParam = param;
}

void Ext0NuTimeRes::BeginOfRun(DS::Run&) {
    string tableName;
    string materialIndex;

    if (fParam.empty()) {
        tableName = "0NUEXT_TIMERESID_Tl208AV"; //default background
        try {
            materialIndex = DB::Get()->GetLink("GEO","inner_av")->GetS("material");
            DB::Get()->GetLink(tableName, materialIndex)->GetD("time_first");
        } catch (...) {
            materialIndex = ""; //fallback scintillator
        }
    }else {
        const size_t hyphenPos = fParam.find("-");
        if (hyphenPos != string::npos) {
            tableName = "0NUEXT_TIMERESID_"+fParam.substr(0, hyphenPos);
            materialIndex = fParam.substr(hyphenPos+1);
        } else {
            tableName = "0NUEXT_TIMERESID_"+fParam;
            try {
                materialIndex = DB::Get()->GetLink("GEO","inner_av")->GetS("material");
                DB::Get()->GetLink(tableName, materialIndex)->GetD("time_first");
            } catch (...) {
                materialIndex = ""; //fallback scintillator
            }
        }
    }
    info << "Ext0NuTimeRes: Using " << tableName << "[" << materialIndex << "] weights" << endl;
    DBLinkPtr pdftable = DB::Get()->GetLink(tableName, materialIndex);

    fTimeFirst = pdftable->GetD("time_first");
    fTimeLast = pdftable->GetD("time_last");
    fTimeStep = pdftable->GetD("time_step");
    fNormalisationFlag = pdftable->GetI("normalisation_flag");
    fOffset = pdftable->GetI("offset");

    fWeights = pdftable->GetDArray("weights");
}


void Ext0NuTimeRes::EndOfRun(DS::Run&) {
}

void Ext0NuTimeRes::DefaultSeed() {
    fEventPos = TVector3();
    fEventTime = 0.0;
}

void Ext0NuTimeRes::SetSeed(const DS::FitResult& seed) {
    fEventPos = seed.GetVertex(0).GetPosition();
    fEventTime = seed.GetVertex(0).GetTime();
}

DS::ClassifierResult Ext0NuTimeRes::GetClassification() {
    const PMTInfo &pmtInfo = Utility::Get()->GetPMTInfo();
    const EffectiveVelocity &effectiveVelocity = Utility::Get()->GetEffectiveVelocity();
    LightPathCalculator lightPathCalc = Utility::Get()->GetLightPathCalculator();

    //Reset from previous event
    fClassifierResult.Reset();

    //No PMTs -> do nothing
    if (fPMTData.empty())
        return fClassifierResult;

    //Don't classify events less hits then offset
    if (static_cast<int>(fPMTData.size())<fOffset){
        return fClassifierResult;
    }

    //Reset discriminant
    fDiscriminant = 0.0;

    std::vector<double> timeResiduals;
    timeResiduals.reserve(fPMTData.size());

    //Loop over PMTs and calculate time residuals
    for (vector<FitterPMT>::const_iterator pmt = fPMTData.begin(); pmt != fPMTData.end(); ++pmt) {
        // Calculate the hit time residuals
        const TVector3 pmtPos = pmtInfo.GetPosition(pmt->GetID());
        lightPathCalc.CalcByPosition(fEventPos, pmtPos);
        const double transitTime = effectiveVelocity.CalcByDistance(lightPathCalc.GetDistInInnerAV(), lightPathCalc.GetDistInAV(), lightPathCalc.GetDistInWater());
        timeResiduals.push_back(pmt->GetTime() - transitTime - fEventTime);

    }

    //Sort timeresiduals and get the offset time
    double timeOffset =0.0;
    if (fOffset!=-1){
        std::sort(timeResiduals.begin(), timeResiduals.end());
        timeOffset = timeResiduals.at(fOffset-1);
    }

  //Loop over time residuals and get the time bin
    for (size_t i = 0; i< timeResiduals.size(); i++){
        double offsetResidual = timeResiduals.at(i) - timeOffset;
        size_t bin;
        if (offsetResidual>fTimeLast){
            bin = static_cast<size_t>((fTimeLast-fTimeFirst)/fTimeStep) -1; //overflow bin
        }else if (offsetResidual<fTimeFirst){
            bin = 0; //underflow bin
        }else{
            bin = static_cast<size_t>((offsetResidual-fTimeFirst)/fTimeStep);
        }

        //Add the weights
        fDiscriminant += fWeights[bin];
    }

    //Normalise to one hit
  if (fNormalisationFlag)
    fDiscriminant /= static_cast<double>(fPMTData.size());

    //Store & return results
    fClassifierResult.SetClassification("discriminant", fDiscriminant);
    fClassifierResult.SetValid(true);
    return fClassifierResult;
}