///////////////////////////////////////////////////////////////////////////////
/// \file cerenkov.cc
///
/// \brief Extracts Cerenkov photon properties from tracking
///
/// \author P G Jones <p.g.jones@qmul.ac.uk>
///
/// REVISION HISTORY:\n
///     2014-07-21 : P G Jones - Moved to rat-tools validation.\n
///     2016-10-19 : M Stringer - Changes to how MCTracks are accessed (PR #1508).
//      2018-02-10 : R Lane - change to function arguments necessary for ROOT 6 compatibility
///
/// \details The Number, wavelength direction and time of Cerenkov photons are
/// extracted.
///
///////////////////////////////////////////////////////////////////////////////
#include <RAT/DU/DSReader.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/MC.hh>
#include <RAT/DS/MCTrack.hh>
#include <RAT/DS/MCTrackStep.hh>

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

void cerenkov(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* hNum = new TH1D( "Num", "Number of emitted photons", 200, 0.0, 4000.0 );
  hNum->SetXTitle( "Number of photons per event." );
  hNum->SetYTitle( "Count per 20 photon bin." );
  TH1D* hWave = new TH1D( "Wave", "Wavelength of emitted photons", 80, 100.0, 900.0 );
  hWave->SetXTitle( "Wavelength nm." );
  hWave->SetYTitle( "Count per 10nm bin." );
  TH1D* hDir = new TH1D( "Dir", "Direction of emitted photons relative to initial electron", 10, -1.0, 1.0 );
  hDir->SetXTitle( "cos(#theta)." );
  hDir->SetYTitle( "Count per 0.1 cos(#theta) bin." );
  TH1D* hTime = new TH1D( "Time", "Time of photon emission", 10, 0.0, 10.0 );
  hTime->SetXTitle( "Emission time ns." );
  hTime->SetYTitle( "Count per 1 ns bin." );

  for( size_t iEntry = 0; iEntry < dsReader->GetEntryCount(); iEntry++ )
    {
      const RAT::DS::Entry& rDS = dsReader->GetEntry( iEntry );
      const RAT::DS::MC& rMC = rDS.GetMC();
      std::vector<unsigned int> mcTrackIDs = rMC.GetMCTrackIDs();
      const TVector3 initialDir = rMC.GetMCParticle( 0 ).GetMomentum().Unit();

      unsigned int numPhotons = 0;
      for( size_t iMCTrack = 0; iMCTrack < mcTrackIDs.size(); iMCTrack++ )
        {
          const RAT::DS::MCTrack& rMCTrack = rMC.GetMCTrack( mcTrackIDs.at(iMCTrack) );
          if( rMCTrack.GetPDGCode() != 0 ) // Optical photon is 0
            continue; // Only want optical photon tracks
          const RAT::DS::MCTrackStep& initialStep = rMCTrack.GetMCTrackStep( 0 );
          if( initialStep.GetProcess() != string( "Cerenkov" ) )
            continue; // Only want Cerenkov photons
          numPhotons++;
          hWave->Fill( 0.001242375 / initialStep.GetKineticEnergy() ); // MeV to nm
          hDir->Fill( initialStep.GetMomentum().Unit().Dot( initialDir ) );
          hTime->Fill( initialStep.GetGlobalTime() );
        }
      hNum->Fill( numPhotons );
    }

  outtfile->cd();
  hNum->Write();
  hWave->Write();
  hDir->Write();
  hTime->Write();
  outtfile->Close();
  delete outtfile;
  delete dsReader;
}