#include <ITpcTrackSandMuonControlSample.hxx>
#include <IReconPID.hxx>

ITpcTrackSandMuonControlSample::ITpcTrackSandMuonControlSample() : IControlSampleBase("TpcTrackSandMuon"){

  AddValidTrigger(COMET::Trigger::kBeamSpill); 

  fTotPOT    = 0;
  fPrevPOT   = 0;
  fNSelected = 0;

}

void ITpcTrackSandMuonControlSample::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);
  hNFgd1Tracks = new TH1F("hNFgd1Tracks", ";N_{TRACKS};", 20, -0.5, 19.5);
  hNFgd2Tracks = new TH1F("hNFgd2Tracks", ";N_{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);
  hNTpc3_12    = new TH1F("hNTpc3_12"   , ";N_{TRACKS};", 20, -0.5, 19.5);
  hNTpc2_13    = new TH1F("hNTpc2_13"   , ";N_{TRACKS};", 20, -0.5, 19.5);
  hNTpc1_23    = new TH1F("hNTpc1_23"   , ";N_{TRACKS};", 20, -0.5, 19.5);

}

bool ITpcTrackSandMuonControlSample::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;
  }

  COMET::IHandle<COMET::IAlgorithmResult> TpcResult = event.GetFit("ReconTPC");
  if(!TpcResult) return false;
    
  COMET::IHandle<COMET::IReconObjectContainer> TrackCont = TpcResult->GetResultsContainer("TPCPids");
  if(!TrackCont) return false;
  
  
  int ntpc1 =0, ntpc2 = 0, ntpc3 =0;
  //  int total_tpc = TrackCont->size();
  for( COMET::IReconObjectContainer::iterator  tr = TrackCont->begin(); tr < TrackCont->end();tr++) { 
    
    COMET::IHandle<COMET::IReconPID> trackPiD = *tr;    
    // Don't count the delta rays
    if ( trackPiD->GetMomentum() < 50. )
      continue;

    if(trackPiD->UsesDetector(COMET::IReconBase::kTPC1)) ntpc1++;
    if(trackPiD->UsesDetector(COMET::IReconBase::kTPC2)) ntpc2++;
    if(trackPiD->UsesDetector(COMET::IReconBase::kTPC3)) ntpc3++;    
    
  }

  if(fMakeHistoFile){
    hNTpc1Tracks->Fill(ntpc1);
    hNTpc2Tracks->Fill(ntpc2);
    hNTpc3Tracks->Fill(ntpc3);
  }

  //if(ntpc1 <= 1 && ntpc2 <= 1 && ntpc3 <= 1){
  //  if(total_tpc >= 2) return true;
  //}

  // require at least two of the tpcs to have exactly one track
  bool pass = (ntpc1 == 1 && ntpc2 == 1) || (ntpc1 == 1 && ntpc3 == 1) || (ntpc2 == 1 && ntpc3 == 1);

  if(pass){
    fNSelected++;

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

      if(ntpc1 == 1 && ntpc2 == 1)hNTpc3_12->Fill(ntpc3);
      if(ntpc1 == 1 && ntpc3 == 1)hNTpc2_13->Fill(ntpc2);
      if(ntpc2 == 1 && ntpc3 == 1)hNTpc1_23->Fill(ntpc1);

      int nfgd1 =0, nfgd2 = 0;
      COMET::IHandle<COMET::IAlgorithmResult> FgdResult = event.GetFit("ReconFGD");
      if(FgdResult){
  COMET::IHandle<COMET::IReconObjectContainer> TrackCont = FgdResult->GetResultsContainer("fittedFgdReconTracks_allFGD");
      if(TrackCont){
    for( COMET::IReconObjectContainer::iterator  tr = TrackCont->begin(); tr < TrackCont->end();tr++) {     
      COMET::IHandle<COMET::IReconTrack> track = *tr;
      if(track->UsesDetector(COMET::IReconBase::kFGD1)) nfgd1++;
      if(track->UsesDetector(COMET::IReconBase::kFGD2)) nfgd2++;    
    }
  }
      }
      hNFgd1Tracks->Fill(nfgd1);
      hNFgd2Tracks->Fill(nfgd2);
    }
    fPrevPOT = fTotPOT;
  }

  return pass;

}