////////////////////////////////////////////////////////////////////
/// \file ELLIE.cc
///
/// \brief Extracts fibre information from file for testing.
///
/// \author J Wilson <j.r.wilson@qmul.ac.uk>
///
/// REVISION HISTORY:\n
///     2014-05-29 : P G Jones - Added header information.\n
//      2018-02-10 : R Lane - change to function arguments necessary for ROOT 6 compatibility
///
/// \details This extracts the photon angle from beam, wavelength and
/// times.
///
////////////////////////////////////////////////////////////////////

#include <RAT/DU/DSReader.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/MC.hh>
#include <RAT/DS/MCParticle.hh>

#include <TH1D.h>
#include <TFile.h>

#include <math.h>

void ELLIE(std::string event_file, std::string outfile)
{
  RAT::DU::DSReader* dsReader = new RAT::DU::DSReader(event_file);
  TFile *outtfile = new TFile(outfile.c_str(),"RECREATE");

  TH1D *hang = new TH1D("hang","FT055B, TELLIE503 Angle",180,0,60);
  hang->SetXTitle("Photon Angle from beam (degrees)");
  hang->SetYTitle("Amplitude");

  TH1D *hwave = new TH1D("hwave","TELLIE503 wavelength",200,400,600);
  hwave->SetXTitle("Photon Wavelength (nm)");
  hwave->SetYTitle("Amplitude");

  TH1D *htime = new TH1D("htime","TELLIE503 time",100,-20.,20.);
  htime->SetXTitle("Photon Time (ns)");
  htime->SetYTitle("Amplitude");

  const TVector3 fibredir(-1.3e-05,0.894427,0.447213); // Load from FIBREInfo?? PHIL FIXME

  for( size_t iEntry = 0; iEntry < dsReader->GetEntryCount(); iEntry++ )
    {
      const RAT::DS::Entry& rDS = dsReader->GetEntry( iEntry );
      const RAT::DS::MC& rMC = rDS.GetMC();
      for( size_t iParticle = 0; iParticle < rMC.GetMCParticleCount(); iParticle++ )
        {
          const RAT::DS::MCParticle& rMCParticle = rMC.GetMCParticle( iParticle );
          const TVector3 mcMom = rMCParticle.GetMomentum();
          const double mcEnergy = rMCParticle.GetKineticEnergy();
          const double mcTime = rMCParticle.GetTime();
          hwave->Fill(1.23984193e-3/mcEnergy);
          htime->Fill(mcTime);
          hang->Fill(180/M_PI*mcMom.Angle(fibredir));
        }
    }
  outtfile->cd();
  hang->Write();
  hwave->Write();
  htime->Write();
  outtfile->Close();
  delete outtfile;
  delete dsReader;
}