////////////////////////////////////////////////////////////////////////////////
/// \file fitquad.cc
///
/// \brief This plots the fitter performance
///
/// \author P Jones <p.g.jones@qmul.ac.uk>
///
/// REVISION HISTORY:\n
///     2014-06-16 : P G Jones - new file.\n
//      2018-02-10 : R Lane - change to function arguments necessary for ROOT 6 compatibility
///
/// \details The performance is defined as the resolution in both
/// position and energy verus radius, energy etc...
///
////////////////////////////////////////////////////////////////////////////////
#include <RAT/DU/DSReader.hh>
#include <RAT/DS/Entry.hh>

#include <TH1D.h>
#include <TF1.h>
#include <TFile.h>
#include <TMath.h>
#include <TProfile.h>

void fitquad(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 *hCenterRes = new TH1D("hCenterRes", "1 MeV electrons at the center;Reconstructed X (mm);Events per bin", 100, -6000, 6000);
  TProfile *hReconVsRadius = new TProfile("hReconVsRadius", "1 MeV electrons (error bars are RMS);(Radius/6005 mm)**3;Reconstructed X - True X (mm)", 5, 0, 1, -750.0, 750.0, "s");
  TProfile *hReconVsEnergy = new TProfile("hReconVsEnergy", "Electrons at center (error bars are RMS);Kinetic energy (MeV);Reconstructed X (mm)", 5, 0.9, 3.6, -750.0, 750.0, "s");

  for( size_t iEntry = 0; iEntry < dsReader->GetEntryCount(); iEntry++ )
    {
      const RAT::DS::Entry& rDS = dsReader->GetEntry( iEntry );
      const TVector3 mcPos = rDS.GetMC().GetMCParticle(0).GetPosition();
      const double mcEnergy = rDS.GetMC().GetMCParticle(0).GetKineticEnergy();

      for( size_t iEV = 0; iEV < rDS.GetEVCount(); iEV++ )
        {
          const RAT::DS::EV& rEV = rDS.GetEV( iEV );
          if( !rEV.FitResultExists( "quad" ) ||  rEV.GetFitResult( "quad" ).GetVertexCount() == 0 ||
              !rEV.GetFitResult( "quad" ).GetValid() )
            continue;
          const TVector3 fitPos = rEV.GetFitResult( "quad" ).GetVertex(0).GetPosition();

          if( mcEnergy < 1.05 && mcPos.X() < 0.5 )
            hCenterRes->Fill( fitPos.X() );
          if( mcEnergy > 1.05 && mcPos.X() < 0.5 )
            hReconVsEnergy->Fill( mcEnergy, fitPos.X() );
          if( mcEnergy < 1.05 && mcPos.X() > 0.5 )
            hReconVsRadius->Fill( TMath::Power( mcPos.Mag() / 6005, 3.0 ), fitPos.X() - mcPos.X() );
        }
    }

  TF1 *g1    = new TF1("g1","gaus",-2000,2000);
  hCenterRes->Fit(g1,"R");
  hReconVsRadius->SetMinimum(-600);
  hReconVsRadius->SetMaximum(600);
  hReconVsEnergy->SetMinimum(-600);
  hReconVsEnergy->SetMaximum(600);

  outtfile->cd();
  hCenterRes->Write();
  hReconVsRadius->Write();
  hReconVsEnergy->Write();
  outtfile->Close();
  delete outtfile;
  delete dsReader;
}