////////////////////////////////////////////////////////////////////
/// \file PlotPMTResponse.cc
///
/// \brief Functions to select and plot MCPhoton features
///
/// \author P G Jones <p.g.jones@qmul.ac.uk>
///
/// REVISION HISTORY:\n
///     2014-05-26: P. Jones - First Revision.\n
///
/// \details Functions to plot the PMT response.
///
////////////////////////////////////////////////////////////////////

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

#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 <RAT/DB.hh>

#include <TH1D.h>
#include <TStyle.h>
#include <TMath.h>
#include <TLegend.h>

using namespace std;


/// Plot the photoelectron and reflection responses
///
/// The photoelectron response is the number of photoelectrons/number of hits per bin.
/// The reflection response is the number of reflected photons/number of hits per bin.
///
/// @param[in] fileName of the RAT::DS root file to analyse
/// @return the histogram plot
TH1D* PlotAngularResponse( const string& fileName )
{
  gStyle->SetOptStat(111111); // to show overflow / underflow stats

  TH1D* hResponse = new TH1D( "hResponse", "PMT Response as function of angle", 90, 0.0, 90.0 );
  TH1D* hReflected = new TH1D( "hReflected", "PMT Response as function of angle", 90, 0.0, 90.0 );
  TH1D* hHits = new TH1D( "hHits", "PMT hits as function of angle", 90, 0.0, 90.0 );
  // If this is being done on data that does not require remote database connection
  // eg.: a simple simulation with default run number (0)
  // We can disable the remote connections:
  //
  // NOTE: Don't do this if you are using real data!!!
  RAT::DB::Get()->SetAirplaneModeStatus(true);

  RAT::DU::DSReader dsReader( fileName );
  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 );
              const double xyRadius = sqrt( mcPhoton.GetInPosition().X() * mcPhoton.GetInPosition().X() + mcPhoton.GetInPosition().Y() * mcPhoton.GetInPosition().Y() );
              if( mcPhoton.GetInPosition().Z() < 132.0 || xyRadius > 137.0 )
                continue;
              const double angle = TMath::ACos( mcPhoton.GetInDirection().Dot( pmtDirection ) ) * TMath::RadToDeg();
              hHits->Fill( angle );
              if( mcPhoton.GetFate() == RAT::DS::MCPhoton::ePhotoelectron )
                hResponse->Fill( angle );
              else if( mcPhoton.GetFate() == RAT::DS::MCPhoton::eReflected )
                hReflected->Fill( angle );
            }
        }
    }
  hResponse->Sumw2();
  hHits->Sumw2();
  hReflected->Sumw2();
  hResponse->Divide( hHits );
  hResponse->Draw("E");
  hReflected->Divide( hHits );
  hReflected->SetLineColor( kRed );
  //hReflected->Draw("SAMEE");
  TLegend* t1 = new TLegend( 0.7, 0.7, 0.88, 0.88 );
  t1->AddEntry( hResponse, "Signal", "l" );
  //t1->AddEntry( hReflected, "Reflection", "l" );
  t1->SetLineColor( kWhite );
  t1->SetFillColor( kWhite );
  t1->Draw();
  return hResponse;
}


/// Plot the photoelectron and reflection responses
///
/// The photoelectron response is the number of photoelectrons/number of hits per bin.
/// The reflection response is the number of reflected photons/number of hits per bin.
///
/// @param[in] fileName of the RAT::DS root file to analyse
/// @return the histogram plot
TH1D* PlotPositionResponse( const string& fileName )
{
  gStyle->SetOptStat(111111); // to show overflow / underflow stats

  TH1D* hResponse = new TH1D( "hResponse", "PMT Response as function of position", 200, 0.0, 200.0 );
  TH1D* hReflected = new TH1D( "hReflected", "PMT Response as function of position", 200, 0.0, 200.0 );
  TH1D* hHits = new TH1D( "hHits", "PMT hits as function of position", 200, 0.0, 200.0 );
  // If this is being done on data that does not require remote database connection
  // eg.: a simple simulation with default run number (0)
  // We can disable the remote connections:
  //
  // NOTE: Don't do this if you are using real data!!!
  RAT::DB::Get()->SetAirplaneModeStatus(true);

  RAT::DU::DSReader dsReader( fileName );
  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() < 132.0 )
                continue;
              const double xyRadius = sqrt( mcPhoton.GetInPosition().X() * mcPhoton.GetInPosition().X() + mcPhoton.GetInPosition().Y() * mcPhoton.GetInPosition().Y() );
              hHits->Fill( xyRadius );
              if( mcPhoton.GetFate() == RAT::DS::MCPhoton::ePhotoelectron )
                hResponse->Fill( xyRadius );
              else if( mcPhoton.GetFate() == RAT::DS::MCPhoton::eReflected )
                hReflected->Fill( xyRadius );
            }
        }
    }
  hResponse->Sumw2();
  hHits->Sumw2();
  hReflected->Sumw2();
  hResponse->Divide( hHits );
  hResponse->Draw("E");
  hReflected->Divide( hHits );
  hReflected->SetLineColor( kRed );
  hReflected->Draw("SAMEE");
  TLegend* t1 = new TLegend( 0.7, 0.7, 0.88, 0.88 );
  t1->AddEntry( hResponse, "Signal", "l" );
  t1->AddEntry( hReflected, "Reflection", "l" );
  t1->SetLineColor( kWhite );
  t1->SetFillColor( kWhite );
  t1->Draw();
  return hResponse;
}

/// Plot the photoelectron and reflection responses
///
/// The photoelectron response is the number of photoelectrons/number of hits per bin.
/// The reflection response is the number of reflected photons/number of hits per bin.
///
/// @param[in] fileName of the RAT::DS root file to analyse
/// @return the histogram plot
TH1D* PlotWavelengthResponse( const string& fileName )
{
  gStyle->SetOptStat(111111); // to show overflow / underflow stats

  TH1D* hResponse = new TH1D( "hResponse", "PMT Response as function of wavelength", 600, 200.0, 800.0 );
  TH1D* hReflected = new TH1D( "hReflected", "PMT Response as function of wavelength", 600, 200.0, 800.0 );
  TH1D* hHits = new TH1D( "hHits", "PMT hits as function of wavelength", 600, 200.0, 800.0 );
  // If this is being done on data that does not require remote database connection
  // eg.: a simple simulation with default run number (0)
  // We can disable the remote connections:
  //
  // NOTE: Don't do this if you are using real data!!!
  RAT::DB::Get()->SetAirplaneModeStatus(true);

  RAT::DU::DSReader dsReader( fileName );
  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() < 132.0 )
                continue;
              const double wavelength = twopi * 197.32705e-6 / mcPhoton.GetEnergy(); // In nm
              hHits->Fill( wavelength );
              if( mcPhoton.GetFate() == RAT::DS::MCPhoton::ePhotoelectron )
                hResponse->Fill( wavelength );
              else if( mcPhoton.GetFate() == RAT::DS::MCPhoton::eReflected )
                hReflected->Fill( wavelength );
            }
        }
    }
  hResponse->Sumw2();
  hHits->Sumw2();
  hReflected->Sumw2();
  hResponse->Divide( hHits );
  hResponse->Draw("E");
  hReflected->Divide( hHits );
  hReflected->SetLineColor( kRed );
  hReflected->Draw("SAMEE");
  TLegend* t1 = new TLegend( 0.7, 0.7, 0.88, 0.88 );
  t1->AddEntry( hResponse, "Signal", "l" );
  t1->AddEntry( hReflected, "Reflection", "l" );
  t1->SetLineColor( kWhite );
  t1->SetFillColor( kWhite );
  t1->Draw();
  return hResponse;
}