////////////////////////////////////////////////////////////////////
/// \file PlotSunDirection.cc
///
/// \brief Function to show how to use the SunDirection function
///
/// \author N Barros <nfbarros@hep.upenn.edu>
///
/// REVISION HISTORY:\n
///     2015-02-10 : N Barros - First Revision.\n
///
/// \details This file contains functions which show how to use the different utils
/// 		 related to solar data.
///
/// Contains: PlotCosThetaSun( fileName, fitName )
///
////////////////////////////////////////////////////////////////////

#include <RAT/DU/DSReader.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/EV.hh>
#include <RAT/DS/PMT.hh>
#include <RAT/DS/FitResult.hh>
#include <RAT/DS/UniversalTime.hh>
#include <RAT/SunUtil.hh>
#include <RAT/DB.hh>

#include <TH1D.h>
#include <TMath.h>
#include <TCanvas.h>
#include <TVector3.h>

#include <string>

/// Plot cosine of the angle between the reconstructed direction and the direction to the Sun.
///
/// @param[in] fileName of the RAT::DS root file to analyse
/// @param[in] fitName of the fitter whose results we want to use. By default it uses waterFitter
void PlotCosThetaSun( const std::string& fileName, const std::string& fitName )
{
  TH1D* hCosThetaSun = new TH1D( "hCosThetaSun", "Event cos(#theta_{#odot})", 100, -1.0, 1.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 );
  // Loop through entrys in rootfile
  for( size_t iEntry = 0; iEntry < dsReader.GetEntryCount(); iEntry++ )
    {
      const RAT::DS::Entry& rDS = dsReader.GetEntry( iEntry );

      // Look at each triggered event in the entry
      for( size_t iEV = 0; iEV < rDS.GetEVCount(); iEV++ )
        {
          // Get the EV
          const RAT::DS::EV& rEV = rDS.GetEV( iEV );

          // Check for non-existant or poor fits
          if( !rEV.FitResultExists(fitName) ) continue;                                // No fit vertex exists
          if( !rEV.GetFitResult(fitName).GetVertex(0).ContainsDirection() ) continue;  // Fit vertex doesn't contain a direction fit
          if( !rEV.GetFitResult(fitName).GetVertex(0).ValidDirection() ) continue;     // Direction fit was marked as an invalid fit

          // Should be left with only the events with successful direction fits from the "fitName" fitter

          RAT::DS::FitResult fResult = rEV.GetFitResult(fitName);               // Get fit result
          RAT::DS::FitVertex fVertex = fResult.GetVertex(0);                    // Get first fit vertex
          TVector3 fDirection = fVertex.GetDirection();                         // Get reconstructed direction

          // Get the direction
          RAT::DS::UniversalTime fTime = rEV.GetUniversalTime();

          TVector3 fSunDir = RAT::SunDirection(fTime.GetDays(),fTime.GetSeconds(),fTime.GetNanoSeconds());
          hCosThetaSun->Fill( fDirection.Dot(-1.0*fSunDir));
        }
    }
  TCanvas* can = new TCanvas("cDir","Direction Plots", 600, 600);
  can->cd();
  hCosThetaSun->GetYaxis()->SetTitle( "Counts" );
  hCosThetaSun->GetXaxis()->SetTitle( "#cos(#theta_{#odot})" );
  hCosThetaSun->Draw();
}