////////////////////////////////////////////////////////////////////
/// \file fitcentroid.cc
///
/// \brief This plots the fitter performance
///
/// \author A Mastbaum <mastbaum@hep.upenn.edu>
///
/// REVISION HISTORY:\n
///     2014-05-29 : P G Jones - Added header information.\n
///
/// \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 <TFile.h>
#include <TMath.h>
#include <TProfile.h>

void make_plots( RAT::DU::DSReader* dsReader, TFile *out_file )
{
  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( "centroid" ) ||  rEV.GetFitResult( "centroid" ).GetVertexCount() == 0 ||
              !rEV.GetFitResult( "centroid" ).GetValid() )
            continue;
          const TVector3 fitPos = rEV.GetFitResult( "centroid" ).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() );
        }
    }

  hCenterRes->Fit("gaus");
  hReconVsRadius->SetMinimum(-600);
  hReconVsRadius->SetMaximum(600);
  hReconVsEnergy->SetMinimum(-600);
  hReconVsEnergy->SetMaximum(600);

  out_file->cd();
  hCenterRes->Write();
  hReconVsRadius->Write();
  hReconVsEnergy->Write();
}