#include <CLHEP/Units/PhysicalConstants.h>
using CLHEP::pi; using CLHEP::twopi; using CLHEP::halfpi;

#include <RAT/TellieMonitorProc.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/Run.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>
using namespace RAT;

#include <TMath.h>
#include <TFile.h>
#include <TH2F.h>

#include <math.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
using namespace std;




TellieMonitorProc::TellieMonitorProc(): Processor("TellieMonitorProc")
{
  fOutputName = "TELLIE_PROC";
  fNoPlot = false;
}


TellieMonitorProc::~TellieMonitorProc()
{
  if(fHitMap)
    delete fHitMap;
}


void TellieMonitorProc::SetS( const std::string& param,
                              const std::string& value )
{
  if( param == "file" )
    fOutputName = value;
  else
    throw ParamUnknown( param );
}


void TellieMonitorProc::SetI( const std::string& param,
                              const int value )
{
  if( param == string("noPlot") )
    {
      if( value > 0 )
        fNoPlot = true;
      else
        fNoPlot = false;
    }
  else
    throw ParamUnknown( param );
}


void TellieMonitorProc::BeginOfRun(DS::Run&)
{
  fTrigBit = 0x8000;
  // Set any old elements to zero
  fill(fHitCount.begin(), fHitCount.end(), 0);
  // Resize for the given run
  fHitCount.resize(DU::Utility::Get()->GetPMTInfo().GetCount(), 0);
  fNEvents = 0;
  fHitMap = new TH2F("TellieHitMap", ";Phi;Theta", 360, -pi, pi, 180, 0, pi);
  fHitMap->SetDirectory(NULL);
}

void TellieMonitorProc::EndOfRun(DS::Run& run)
{
  // Get the most frequently hit PMT
  int maxPMT = distance(fHitCount.begin(), max_element(fHitCount.begin(), fHitCount.end()));
  const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
  TVector3 maxPMTPos = pmtInfo.GetPosition(maxPMT);
  int numDirect = 0;
  int numReflected = 0;
  double posDirect[3] = {0.0, 0.0, 0.0};
  double posReflected[3] = {0.0, 0.0, 0.0};

  for(unsigned int i=0;i<fHitCount.size();i++)
    {
      if(pmtInfo.GetType(i) != DU::PMTInfo::NORMAL || fHitCount[i]==0)
        continue;

      TVector3 pmtPos = pmtInfo.GetPosition(i);
      if(pmtPos.Dot(maxPMTPos)>0)
        {
          // direct
          numDirect += fHitCount[i];
          posDirect[0] += fHitCount[i] * pmtPos.X() / pmtPos.Mag();
          posDirect[1] += fHitCount[i] * pmtPos.Y() / pmtPos.Mag();
          posDirect[2] += fHitCount[i] * pmtPos.Z() / pmtPos.Mag();
        }
      else
        {
          // reflected
          numReflected += fHitCount[i];
          posReflected[0] += fHitCount[i] * pmtPos.X() / pmtPos.Mag();
          posReflected[1] += fHitCount[i] * pmtPos.Y() / pmtPos.Mag();
          posReflected[2] += fHitCount[i] * pmtPos.Z() / pmtPos.Mag();
        }
    }

  double radReflected = 0;
  double radDirect = 0;
  for(int i=0;i<3;i++)
    {
      posDirect[i] /= static_cast<double>(numDirect);
      posReflected[i] /= static_cast<double>(numReflected);
      radDirect += (posDirect[i] * posDirect[i]);
      radReflected += (posReflected[i] * posReflected[i]);
    }
  radDirect = sqrt(radDirect);
  radReflected = sqrt(radReflected);

  TVector3 direct;
  TVector3 reflected;

  direct.SetX(posDirect[0]);
  direct.SetY(posDirect[1]);
  direct.SetZ(posDirect[2]);

  reflected.SetX(posReflected[0]);
  reflected.SetY(posReflected[1]);
  reflected.SetZ(posReflected[2]);

  stringstream ss;
  ss << std::setfill('0') << std::setw(10) << run.GetRunID();

  ofstream tellieOut((fOutputName + "_" + ss.str() + ".dat").c_str());
  tellieOut << "Tellie triggers " << fNEvents << endl;
  tellieOut << "Direct\t" << direct.Phi() << "\t" << direct.Theta() << endl;
  tellieOut << "Reflected\t" << reflected.Phi() << "\t" << reflected.Theta() << endl;
  tellieOut.close();

  if(!fNoPlot)
    {
      TFile tellieFile((fOutputName + "_" + ss.str() + ".root").c_str(), "recreate");
      tellieFile.cd();
      fHitMap->Write();
      tellieFile.Close();
    }

  // finally, delete newed
  delete fHitMap;
  fHitMap = NULL;

}

Processor::Result TellieMonitorProc::DSEvent(DS::Run&, DS::Entry& ds)
{
  const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
  for( size_t iEV = 0; iEV < ds.GetEVCount(); iEV++ )
    {
      DS::EV& ev = ds.GetEV( iEV );
      // Only run for events with the correct trigger bit
      if(!(ev.GetTrigType() & fTrigBit))
        return OK;

      for(size_t ipmt = 0; ipmt < ev.GetUncalPMTs().GetCount(); ipmt++ )
        {
          DS::PMTUncal& pmt = ev.GetUncalPMTs().GetPMT( ipmt );
          fHitCount[pmt.GetID()]++;
          if(!fNoPlot)
            fHitMap->Fill(pmtInfo.GetPosition(pmt.GetID()).Phi(), pmtInfo.GetPosition(pmt.GetID()).Theta());
        }
      fNEvents++;
    }
  return OK;
}