#include <IFgdStoppingCosmicControlSample.hxx>

IFgdStoppingCosmicControlSample::IFgdStoppingCosmicControlSample() : IControlSampleBase("FgdStoppingCosmic"){
     fgdTimeBinner = new IFGDTimeBinner();
     fInitialized = false;

     AddValidTrigger(COMET::Trigger::kFGDCosmic);
     AddValidTrigger(COMET::Trigger::kTFBCosmic);

     // Get edge cut (in mm ) from parameters file
     fEdgeCut = COMET::IOARuntimeParameters::Get().GetParameterD("oaControlSample.FgdStoppingCosmic.EdgeCut");
}


void IFgdStoppingCosmicControlSample::InitHistograms(){

  hNTimeBins  = new TH1F("hNTimeBins", ";Time Bins; Events",  10, -0.5,  9.5);
  hNFgdHits   = new TH1F("hNFgdHits" , "; Hits; Events"    , 100, -0.5, 99.5);
  hNFgd1Hits  = new TH1F("hNFgd1Hits", "; Hits; Events"    , 100, -0.5, 99.5);
  hNFgd2Hits  = new TH1F("hNFgd2Hits", "; Hits; Events"    , 100, -0.5, 99.5);
  hT0         = new TH1F("hT0",        "; T_{0} (ns)"      , 100, -10000.0, 10000.0);
  hT1         = new TH1F("hT1",        "; T_{1} (ns)"      , 100, -10000.0, 10000.0);
  hDT10       = new TH1F("hDT10"     , "; #Delta T (ns); Events/100 ns", 100, 0.0, 10000.0);

  return;
}

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


  if(!fInitialized){
    // get plane numbers corresponding to outer planes
    usFGD1XModule = COMET::IGeomInfo::FGD().GetFirstFGD1Plane("X");  // Upstreammost FGD1 X module	
    usFGD1YModule = COMET::IGeomInfo::FGD().GetFirstFGD1Plane("Y");  // Upstreammost FGD1 Y module	
    dsFGD1XModule = COMET::IGeomInfo::FGD().GetLastFGD1Plane("X");   // Downstreammost FGD1 X module
    dsFGD1YModule = COMET::IGeomInfo::FGD().GetLastFGD1Plane("Y");   // Downstreammost FGD1 Y module
    
    usFGD2XModule = COMET::IGeomInfo::FGD().GetFirstFGD2Plane("X");  // Upstreammost FGD2 X module	
    usFGD2YModule = COMET::IGeomInfo::FGD().GetFirstFGD2Plane("Y");  // Upstreammost FGD2 Y module	
    dsFGD2XModule = COMET::IGeomInfo::FGD().GetLastFGD2Plane("X");	  // Downstreammost FGD2 X module
    dsFGD2YModule = COMET::IGeomInfo::FGD().GetLastFGD2Plane("Y");	  // Downstreammost FGD2 Y module
    fInitialized = true;
  }

  COMET::IHandle<COMET::IHitSelection> fgdHitSelect = event.GetHitSelection("fgd");

  // abort if no FGD hits.
  if(!fgdHitSelect)return false;

  // run the time binner
  TFGDHitBins theTimeBinHits = fgdTimeBinner->BinHits(fgdHitSelect);
  Int_t nTimeBin = theTimeBinHits.size();

  // abort if no time bins
  if(nTimeBin == 0)return false;

  // hits on outer FGD modules
  Bool_t usFGD1X = false;
  Bool_t usFGD1Y = false;
  Bool_t dsFGD1X = false;
  Bool_t dsFGD1Y = false;

  Bool_t usFGD2X = false;
  Bool_t usFGD2Y = false;
  Bool_t dsFGD2X = false;
  Bool_t dsFGD2Y = false;

  // hits on peripheral bars
  Bool_t outFGD1Hit = false;
  Bool_t outFGD2Hit = false;

  Int_t nFGD1 = 0;
  Int_t nFGD2 = 0;
  Double_t avgTime0     = 0;
  Double_t totalCharge0 = 0;
  // loop over hits in first time bin
  for (COMET::IHitSelection::const_iterator timeBinHit = theTimeBinHits[0]->begin();
             timeBinHit != theTimeBinHits[0]->end(); timeBinHit++) {  

    Int_t plane = COMET::IGeomInfo::FGD().ActivePlane((*timeBinHit)->GetPosition().Z());
    if     ( plane == usFGD1XModule )usFGD1X = true;
    else if( plane == usFGD1YModule )usFGD1Y = true;
    else if( plane == dsFGD1XModule )dsFGD1X = true;
    else if( plane == dsFGD1YModule )dsFGD1Y = true;

    if     ( plane == usFGD2XModule )usFGD2X = true;
    else if( plane == usFGD2YModule )usFGD2Y = true;
    else if( plane == dsFGD2XModule )dsFGD2X = true;
    else if( plane == dsFGD2YModule )dsFGD2Y = true;
    
    totalCharge0 += (*timeBinHit)->GetCharge();
    avgTime0     += (*timeBinHit)->GetCharge() * (*timeBinHit)->GetTime();
 
    // record if we have hits on the outer edge of the detector (10 cm)
    COMET::IGeometryId theID =  (*timeBinHit)->GetGeomId();
    Int_t FGD = COMET::IGeomInfo::FGD().GetFGD(theID);
    TVector3 hitPos = (*timeBinHit)->GetPosition();
    if(FGD == 0){
      nFGD1++;
      if(COMET::IGeomInfo::FGD().EdgeDistanceFGD1(hitPos.X(), hitPos.Y()) < fEdgeCut){
  outFGD1Hit = true;
      }
    } else if(FGD==1){
      nFGD2++;
      if(COMET::IGeomInfo::FGD().EdgeDistanceFGD2(hitPos.X(), hitPos.Y()) < fEdgeCut){
  outFGD2Hit = true;
      }
    }
    
  }
  avgTime0 /= totalCharge0;
  
  if(fMakeHistoFile) hT0->Fill(avgTime0);
  
  // track comes from downstream through FGD2 and stops in FGD1
  Bool_t stopFGD1 = (dsFGD2X && dsFGD2Y) &&  (usFGD2X && usFGD2Y) && 
                    (dsFGD1X && dsFGD1Y) && !(usFGD1X || usFGD1Y);

  Bool_t stopFGD2 = (usFGD1X && usFGD1Y) &&  (dsFGD1X && dsFGD1Y) && 
                    (usFGD2X && usFGD2Y) && !(dsFGD2X || dsFGD2Y);

  // either is fine, but we don't want edge hits in the stopping FGD
  Bool_t accept   =  (stopFGD1 && !outFGD1Hit) || (stopFGD2 && !outFGD2Hit) ;

  // track comes from downstream through FGD2 and stops in FGD1
  if(accept && fMakeHistoFile){
    hNTimeBins->Fill(nTimeBin);
    hNFgdHits->Fill(theTimeBinHits[0]->size());
    hNFgd1Hits->Fill(nFGD1);
    hNFgd2Hits->Fill(nFGD2);
    
    if(nTimeBin == 2){
      Double_t avgTime1     = 0;
      Double_t totalCharge1 = 0;
      for (COMET::IHitSelection::const_iterator timeBinHit = theTimeBinHits[1]->begin();
     timeBinHit != theTimeBinHits[1]->end(); timeBinHit++) {  
  totalCharge1 += (*timeBinHit)->GetCharge();
  avgTime1     += (*timeBinHit)->GetCharge() * (*timeBinHit)->GetTime();
      }
      avgTime1 /= totalCharge1;
      hT1->Fill(avgTime1);
      hDT10->Fill(avgTime1 - avgTime0);
    }
  }

  return accept;

}