////////////////////////////////////////////////////////////////////
/// \file tefitter.cc
///
/// \brief Extracts performance metrics for scintFitter testing.
///
/// \author I Coulter <icoulter@hep.upenn.edu>
///
/// REVISION HISTORY:\n
///     2016-09-29 : I Coulter - Initial commit.\n
//      2018-02-10 : R Lane - change to function arguments necessary for ROOT 6 compatibility
///
/// \details This extracts the reconstructed position and energy
/// versus mc and each other.
///
////////////////////////////////////////////////////////////////////
#include <RAT/DU/DSReader.hh>
#include <RAT/DS/Entry.hh>

#include <TH1D.h>
#include <TFile.h>
#include <TProfile.h>
#include <TVector3.h>

void tefitter(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 Radius (mm);Events per bin", 100, 0, 1000);
  TProfile *hReconVsRadius = new TProfile("hReconVsRadius", "1 MeV electrons (error bars are RMS);(Radius/6005 mm)**3;Reconstructed Pos - True Pos (mm)", 5, 0, 1, 0, 750.0, "s");
  TProfile *hReconVsEnergy = new TProfile("hReconVsEnergy", "Electrons at center (error bars are RMS);Kinetic energy (MeV);Reconstructed Pos (mm)", 5, 0.9, 3.6, 0, 750.0, "s");
  TProfile *hRadiusEnergy = new TProfile("hRadiusEnergy", "Electrons at center (error bars are RMS);(Radius/6005 mm)**3;Reconstructed energy - True energy (MeV)", 5, 0, 1, -1.0, 1.0, "s");
  TProfile *hCenterEnergy = new TProfile("hCenterEnergy", "Electrons at center (error bars are RMS);Kinetic energy (MeV);Reconstructed energy - True energy (MeV)", 5, 0.9, 3.6, -2.0, 2.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( "scintFitter" ) || rEV.GetFitResult( "scintFitter" ).GetVertexCount() == 0 ||
              !rEV.GetFitResult( "scintFitter" ).GetValid() )
            continue;

          const TVector3 fitPos = rEV.GetFitResult( "scintFitter" ).GetVertex(0).GetPosition();

          if( mcEnergy < 1.05 && mcPos.Mag() < 0.5 )
            hCenterRes->Fill( fitPos.Mag() );
          if( mcEnergy > 1.05 && mcPos.Mag() < 0.5 )
            hReconVsEnergy->Fill( mcEnergy, fitPos.Mag() );
          if( mcEnergy < 1.05 && mcPos.Mag() > 0.5 )
            hReconVsRadius->Fill( pow( mcPos.Mag() / 6005.0, 3.0 ), (fitPos - mcPos).Mag() );

          const double fitEnergy = rEV.GetFitResult("scintFitter").GetVertex(0).GetEnergy();
          if( mcPos.Mag() < 0.5 )
            hCenterEnergy->Fill(mcEnergy, fitEnergy-mcEnergy);
          if( mcEnergy < 1.05 && mcPos.Mag() > 0.5 )
            hRadiusEnergy->Fill( pow(mcPos.Mag()/6005.0, 3.0), fitEnergy-mcEnergy );
        }
    }

  hReconVsRadius->SetMaximum(600);
  hReconVsRadius->SetMinimum(0);
  hReconVsEnergy->SetMaximum(600);
  hReconVsEnergy->SetMinimum(0);
  hCenterEnergy->SetMaximum(2);
  hCenterEnergy->SetMinimum(-2);
  hRadiusEnergy->SetMaximum(1);
  hRadiusEnergy->SetMinimum(-1);

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