////////////////////////////////////////////////////////////////////
/// \file PlotFlatMap.cc
///
/// \brief Example how to plot a detector flatmap
///
/// \author F Descamps <fbdescamps@lbl.gov>
///
/// REVISION HISTORY:\n
///     2015-07-16 : F Descamps - First Revision.\n
///
/// \details This example plots a flatmap of the detector filled
/// with number of hits seen by each PMT.
///
////////////////////////////////////////////////////////////////////

#include <RAT/DU/DSReader.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/MC.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/FlatMap.hh>

#include <TGraph2D.h>
#include <TCanvas.h>
#include <TVector3.h>
#include <TVector2.h>

#include <string>
using namespace std;

TGraph2D* PlotHitPMTs( const string& fileName )
{
   // 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);
 // This plots the detector flatmap with the number of hits per PMT
  RAT::DU::DSReader dsReader( fileName );
  const RAT::DU::PMTInfo& pmtInfo = RAT::DU::Utility::Get()->GetPMTInfo();
  TGraph2D* PMTStatus = new TGraph2D;
  std::vector<int> countHits;
  countHits.resize(pmtInfo.GetCount());
  for( size_t i = 0; i<pmtInfo.GetCount(); i++) countHits[i] = 0;

  for( size_t iEntry = 0; iEntry < dsReader.GetEntryCount(); iEntry++ )
    {
      const RAT::DS::Entry& rDS = dsReader.GetEntry( iEntry );
      for( size_t iEV = 0; iEV < rDS.GetEVCount(); iEV++ )
        {
          const RAT::DS::EV& rEV = rDS.GetEV( iEV );
          const RAT::DS::UncalPMTs& PMTs = rEV.GetUncalPMTs();
          for( size_t iPMT = 0; iPMT < PMTs.GetCount(); iPMT++ )
            {
              const RAT::DS::PMTUncal& pmtUncal = PMTs.GetPMT( iPMT );
              countHits[pmtUncal.GetID()]++;
            }
        }
    }
  // Now fill the flatmap
  for(size_t i = 0; i<pmtInfo.GetCount(); i++)
    {
      // This is the vector to hold the 2D flat position:
      TVector2 PMTPosFlat;
      TVector3 PMTPos = pmtInfo.GetPosition( i );
      const double size = 1.0;
      const double rotation = 2.12;
      // Transform the 3D position of the pmt (in the TVector3 PMTPos)
      // to a 2D position (fill PMTPosFlat)
      RAT::SphereToIcosahedron(PMTPos,PMTPosFlat, size, rotation);
      if (PMTPosFlat.X() != PMTPosFlat.X()) continue; // This checks for NAN
      // Fill the 2D graph, in my case, the z-axis is the number of hits
      PMTStatus->SetPoint(i,PMTPosFlat.X(),PMTPosFlat.Y(),countHits[i]);
    }
  PMTStatus->SetMarkerStyle(21);
  PMTStatus->SetMarkerSize(0.5);
  PMTStatus->Draw("ZCOLPCOL");
  return PMTStatus;
}