#include <IP0DTrackSandMuonControlSample.hxx>
#include <fgdUtils.hxx>
#include <IRawBeamDataSingleton.hxx>
#include <IReconTrack.hxx>
#include <IReconPID.hxx>

IP0DTrackSandMuonControlSample::IP0DTrackSandMuonControlSample() : IControlSampleBase("P0DTrackSandMuon"){

  AddValidTrigger(COMET::Trigger::kBeamSpill);
  fTotPOT = 0;
  fPrevPOT = 0;
  fNSelected = 0;

  // Get Z cut from parameter file for stopping muons
  fZCut = COMET::IOARuntimeParameters::Get().GetParameterD("oaControlSample.P0DTrackSandMuon.ZCut");
}

void IP0DTrackSandMuonControlSample::InitHistograms(){
  gEventsPOT = new TGraph();
  gEventsPOT->SetName("gEventsPOT");
  gDirectory->Append(gEventsPOT);

  hPOTDiff   = new TH1F("hPOTDiff", ";POT to last event (10^{13} POT)", 100, 0.0, 1000.0);
  // From TP0DGeomInfo:
  // fP0DMin (x,y,z)=(-1083.828522,-1107.395000,-3296.476745)
  // fP0DMax (x,y,z)=(1032.953073,1143.495000,-938.753446)
  hNVtx     = new TH1F("hNVtx", ";Vertices;", 20, -0.5, 19.5);
  hNTrk     = new TH1F("hNTrk", ";Tracks;"  , 20, -0.5, 19.5);

  hNTpc1Tracks = new TH1F("hNTpc1Tracks", ";N_{TRACKS};", 20, -0.5, 19.5);
  hNTpc2Tracks = new TH1F("hNTpc2Tracks", ";N_{TRACKS};", 20, -0.5, 19.5);
  hNTpc3Tracks = new TH1F("hNTpc3Tracks", ";N_{TRACKS};", 20, -0.5, 19.5);

  hTrkFront = new TH1F("hTrkFront", ";Track Front (mm)", 237, -3300.0, -930.0);
  hTrkBack  = new TH1F("hTrkBack" , ";Track Front (mm)", 237, -3300.0, -930.0);
  hNSelectedCycles = new TH1F("hNSelectedCycles", ";Cycles;", 49, -0.5, 48.5);
}

bool IP0DTrackSandMuonControlSample::IsEventSelectedInternal(COMET::ICOMETEvent& event) {


  // Get the spill number and time for the event
  int spill = event.GetHeader().GetSpillNumber();
  int mcmTime = event.GetHeader().GetMCMSecond();
  
  // Get the Beam Data
  if(COMET::IRawBeamDataSingleton::Get().GetSpill(spill, mcmTime, 1) ){
    // Get the POT
    double pot = COMET::IRawBeamDataSingleton::Get().GetProtonsPerSpill(4);
    fTotPOT += pot;
  }

  // get outer edges of P0D geometry
  const TVector3& kP0DMinPos=COMET::IGeomInfo::P0D().ActiveMin();
  const TVector3& kP0DMaxPos=COMET::IGeomInfo::P0D().ActiveMax();

  bool select = false;
  int  nSelectedCycles = 0;

  // Loop over Trip-t cycles
  for(int i = 0; i < 24; i++){

    // get the p0d algorithms    
    COMET::IHandle<COMET::IAlgorithmResult> algres = 
      event.GetFit(TString::Format("ReconP0D/p0dCycle%d/IP0DTrackRecon", i));
    if (!algres) continue;


    // Get Vertices
    COMET::IHandle<COMET::IReconObjectContainer> vtxs=algres->GetResultsContainer();
    
    // Require exactly one vertex per cycle
    if( !vtxs || vtxs->size()!=1 ) continue;

    COMET::IHandle<COMET::IReconVertex> vtx=vtxs->front();
    if(!vtx) COMETWarn("Something bad happened - no IReconVertex");

    // Count tracks in vertex

    // N.B. The object here is actually a IReconPID.  This allows us
    // to use the "ReconP0D/p0dCycle%d/IP0DTrackRecon" result, which
    // is a guaranteed path, as oppossed to relying on underlieing
    // algorithms which might change.  As this is referencing the
    // IP0DTrackRecon result, the IReconPID is 100% equivilant to the
    // TReconTracks originally used.
    COMET::IHandle<COMET::IReconObjectContainer> constituents=vtx->GetConstituents();
    int ntrk=0;
    COMET::IHandle<COMET::IReconPID> trk;
    for(int j=0; j < (int)constituents->size(); ++j){
      trk=constituents->at(j);
      if(trk) ++ntrk;
    }

    if( fMakeHistoFile ) hNTrk->Fill(ntrk);

   // Require that the vertex has exactly one track
    if( ntrk!=1 ) continue;
    
    COMET::IHandle<COMET::IHitSelection> trkHits=trk->GetHits();
    
    const TVector3& frontPos=trkHits->front()->GetPosition();
    const TVector3& backPos=trkHits->back()->GetPosition();
    
    // Select tracks entering the face of the P0D and 
    // with Z projection longer than fZCut
    if( fabs(frontPos.Z()- kP0DMinPos.Z()) > 2*unit::cm || 
  //fabs(backPos.Z() - kP0DMaxPos.Z()) > 2*unit::cm)
  backPos.Z() - frontPos.Z() < fZCut )
      continue;

    if( fMakeHistoFile ){
      hTrkFront->Fill(frontPos.z());
      hTrkBack->Fill(backPos.z());
    }

    // If we got this far, we must want to use the whole event (all cycles)
    nSelectedCycles++;
    select = true;
  }

  if( select ){
    fNSelected++;

    if( fMakeHistoFile ){
      gEventsPOT->SetPoint(fNSelected, fNSelected, fTotPOT/1E13);
      hPOTDiff->Fill( (fTotPOT - fPrevPOT)/1E13);
      hNSelectedCycles->Fill(nSelectedCycles);

      int ntpc1 =0, ntpc2 = 0, ntpc3 =0;
      COMET::IHandle<COMET::IAlgorithmResult> TpcResult = event.GetFit("ReconTPC");

      if( TpcResult ){   
  COMET::IHandle<COMET::IReconObjectContainer> TrackCont = TpcResult->GetResultsContainer("TPCPids");

  if( TrackCont ){
    for( COMET::IReconObjectContainer::iterator  tr = TrackCont->begin(); tr < TrackCont->end();tr++) { 
      COMET::IHandle<COMET::IReconPID> trackPiD = *tr;
      if(trackPiD->UsesDetector(COMET::IReconBase::kTPC1)) ntpc1++;
      if(trackPiD->UsesDetector(COMET::IReconBase::kTPC2)) ntpc2++;
      if(trackPiD->UsesDetector(COMET::IReconBase::kTPC3)) ntpc3++;    
    }
  }
      }
      hNTpc1Tracks->Fill(ntpc1);
      hNTpc2Tracks->Fill(ntpc2);
      hNTpc3Tracks->Fill(ntpc3);
    }
    fPrevPOT = fTotPOT;

  }


  return select;
}