///////////////////////////////////////////////////////////////////////////////
/// \file pmtresponse.cc
///
/// \brief Extracts the PMT response plots
///
/// \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
//      2018-02-10 : R Lane - change to function arguments necessary for ROOT 6 compatibility
///
/// \details Extracts the PMT response for testing
///
///////////////////////////////////////////////////////////////////////////////

#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>
#include <TStyle.h>
#include <TMath.h>

#include <math.h>

void pmtresponse(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");

  gStyle->SetOptStat(111111); // to show overflow / underflow stats

  const int thetaBins = 90; const double thetaLow = 0.0; const double thetaHigh = 89.0;
  const int lambdaBins = 49; const double lambdaLow = 220.0; const double lambdaHigh = 710.0;
  TH1D* hWaveResponse = new TH1D( "hWaveResponse", "PMT Response as function of wavelength", lambdaBins, lambdaLow, lambdaHigh );
  TH1D* hResponse = new TH1D( "hResponse", "PMT Response as function of angle", thetaBins, thetaLow, thetaHigh );
  TH1D* hReflected = new TH1D( "hReflected", "PMT Response as function of angle", thetaBins, thetaLow, thetaHigh );
  TH1D* hHits = new TH1D( "hHits", "PMT hits as function of angle", thetaBins, thetaLow, thetaHigh );

  for( size_t iEntry = 0; iEntry < dsReader->GetEntryCount(); iEntry++ )
    {
      const RAT::DS::Entry& rDS = dsReader->GetEntry( iEntry );
      const RAT::DS::MC& mc = rDS.GetMC();
      for( size_t iPMT = 0; iPMT < mc.GetMCPMTCount(); iPMT++ )
        {
          const RAT::DS::MCPMT& mcPMT = mc.GetMCPMT( iPMT );
          const TVector3 pmtDirection = TVector3( 0.0, 0.0, -1.0 );
          for( size_t iPhoton = 0; iPhoton < mcPMT.GetMCPhotonCount(); iPhoton++ )
            {
              const RAT::DS::MCPhoton& mcPhoton = mcPMT.GetMCPhoton( iPhoton );
              if( mcPhoton.GetInPosition().Z() < 122.0 )
                continue;
              const double angle = TMath::ACos( mcPhoton.GetInDirection().Dot( pmtDirection ) ) * TMath::RadToDeg();
              const double wavelength = 2*M_PI * 197.32705e-6 / mcPhoton.GetEnergy(); // In nm
              hHits->Fill( angle );
              if( mcPhoton.GetFate() == RAT::DS::MCPhoton::ePhotoelectron )
                {
                  hResponse->Fill( angle );
                  hWaveResponse->Fill( wavelength );
                }
              else if( mcPhoton.GetFate() == RAT::DS::MCPhoton::eReflected )
                hReflected->Fill( angle );
            }
        }
    }
  outtfile->cd();
  hResponse->Sumw2();
  hHits->Sumw2();
  hReflected->Sumw2();
  hResponse->Divide( hHits );
  hReflected->Divide( hHits );
  hResponse->Write();
  hReflected->Write();
  hWaveResponse->Write();
  outtfile->Close();
  delete outtfile;
  delete dsReader;
}