///////////////////////////////////////////////////////////////////////////////
/// \file fresnel.cc
///
/// \brief Extracts the transmission factor versus angle
///
/// \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 transmission factor versus angle is 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 fresnel(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* temp = new TH1D( "temp", "Temp", 180, -90.0, 90.0 );
  TH1D* transmission = new TH1D( "transmission", "Transmission factor", 180, -90.0, 90.0 );
  transmission->SetXTitle( "Incident angle [degrees]" );
  transmission->SetYTitle( "Transmission factor" );

  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();
      for( size_t iTrack = 0; iTrack < mcTrackIDs.size(); iTrack++ )
        {
          const RAT::DS::MCTrack& rMCTrack = rMC.GetMCTrack( mcTrackIDs.at(iTrack) );
          if( rMCTrack.GetPDGCode() != 0 )
            continue;

          const double startDir = acos( rMCTrack.GetMCTrackStep( 0 ).GetMomentum().Unit().Dot( TVector3( 1.0, 0.0, 0.0 ) ) ) / 3.14 * 180.0;
          temp->Fill( startDir );
          const RAT::DS::MCTrackStep& lastStep = rMCTrack.GetMCTrackStep( rMCTrack.GetMCTrackStepCount() - 1 );
          if( lastStep.GetPosition().x() > 0.0 )
            {
              const double endDir = acos( lastStep.GetMomentum().Unit().Dot( TVector3( 1.0, 0.0, 0.0 ) ) ) / 3.14 * 180.0;
              transmission->Fill( startDir );
            }
        }
    }
  transmission->Divide( temp );
  outtfile->cd();
  transmission->Write();
  outtfile->Close();
  delete outtfile;
  delete dsReader;
}