#include <RAT/EnergyPromptHits.hh>
#include <RAT/DB.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/FitVertex.hh>

#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/DU/ChanHWStatus.hh>

#include <RAT/PMTSelector.hh>
#include <RAT/PMTSelectorFactory.hh>
#include <RAT/FitterPMT.hh>
#include <RAT/SelectorMethod.hh>

#include <cmath>

using namespace RAT;
using namespace RAT::Methods;
using namespace std;

void
EnergyPromptHits::Initialise( const string& param )
{
  fIndex = param;
}

void
EnergyPromptHits::BeginOfRun( DS::Run& run )
{
  DB *db = DB::Get();
  string index = fIndex;
  if( fIndex.empty() )
    {
      try{ index = db->GetLink("GEO","inner_av")->GetS("material"); }
      catch( DBNotFoundError& e ) { /* Will use defaults instead */ }
    }
  // Get default parameters
  DBLinkPtr dbEPHLink = db->GetLink("FIT_ENERGY_PROMPT_HITS");
  DBLinkGroup grp = db->GetLinkGroup("FIT_ENERGY_PROMPT_HITS");
  DBLinkGroup::iterator it;
  // Check for material specific parameters
  // If not, keep default values
  for(it=grp.begin(); it != grp.end(); ++it){
    if(index == it->first){
      dbEPHLink = db->GetLink("FIT_ENERGY_PROMPT_HITS", index);
      break;
    }
  }

  fPromptWindow = dbEPHLink->GetDArray( "prompt_window" );
  fWorkingPMTs = dbEPHLink->GetD( "working_pmts" );
  fNpromptPerMeV = dbEPHLink->GetD("nprompt_per_mev");

  fTimeResidualCut = PMTSelectors::PMTSelectorFactory::Get()->GetPMTSelector( "timeResidualCut" );
  fTimeResidualCut->BeginOfRun( run );
}

void
EnergyPromptHits::EndOfRun( DS::Run& )
{
  fPromptWindow.clear();
}

void
EnergyPromptHits::DefaultSeed()
{
  // Radius = 0.0, direction = -1z, time = 0.0
  DS::FitVertex vertex;
  vertex.SetPosition( TVector3( 0.0, 0.0, 0.0 ), false, true );
  vertex.SetDirection( TVector3( 0.0, 0.0, -1.0 ), false, true );
  vertex.SetTime( 0.0, false, true );
  fSeedResult.SetVertex( 0, vertex );
}

DS::FitResult
EnergyPromptHits::GetBestFit()
{
  fFitResult.Reset();
  if( fPMTData.empty() )
    return fFitResult;
  CopySeedToResult();

  DS::FitVertex fitVertex = fFitResult.GetVertex( 0 );

  // Get number of working PMTs
  const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
  unsigned int numChannels = pmtInfo.GetCount();
  const DU::ChanHWStatus& channelHardwareStatus = DU::Utility::Get()->GetChanHWStatus();
  unsigned int totalActiveChannels = 0;
  for( size_t lcn = 0; lcn < numChannels; lcn++ ){
    if( pmtInfo.GetType( lcn ) != DU::PMTInfo::NORMAL &&
        pmtInfo.GetType( lcn ) != DU::PMTInfo::HQE)
      continue;
    if( channelHardwareStatus.IsEnabled() ){
      if( channelHardwareStatus.IsTubeOnline( lcn ) ){
        totalActiveChannels++;
      }
    }
  }

  // Get prompt PMTs
  fTimeResidualCut->SetD( "lowLimit", fPromptWindow[0] );
  fTimeResidualCut->SetD( "highLimit", fPromptWindow[1] );
  fTimeResidualCut->SetI( "useGroup", 1 );
  vector<FitterPMT> promptPMTs = fTimeResidualCut->GetSelectedPMTs( fPMTData, fitVertex );
  size_t Nprompt = promptPMTs.size();

  // Set fit result
  fitVertex.SetEnergy( ( static_cast<double>( Nprompt ) / fNpromptPerMeV ) * ( fWorkingPMTs / static_cast<double>( totalActiveChannels ) ) );
  fitVertex.SetEnergyErrors( 1.0 );
  fFitResult.SetVertex( 0, fitVertex );
  return fFitResult;
}