#include <IOADatabase.hxx>
#include <cometEventLoop.hxx>
#include <ITFBChannelId.hxx>
#include <ITripTDigitBank.hxx>
#include <ITripTXTimeBank.hxx>
#include <IMidasTripTDigitItr.hxx>
#include <IMidasTripTDigit.hxx>
#include <IRawDataHeader.hxx>
#include <ITripTInfoBank.hxx>
#include <IGeometryDatabase.hxx>
#include <IMCMBank.hxx>
#include <TTree.h>

#include "P0D_GainAnalysis.hxx"


using namespace COMET;


struct pbs_t { 
              Int_t startUnixTime;    //unix time of first event
              Int_t events;           //all events
              Int_t p0dEvents;        //p0d events
              Int_t gt300ADC;         //p0d events with 15 hits > 300ADC
              Int_t fidcut_gt300ADC;  //p0d events within fidicial cut with 15 hits > 300ADC
 }; 
pbs_t pbs;
TTree *treePBS = new TTree("T","p0d beam statistics");
   
bool HasGoodTDC(const COMET::IMidasTripTDigit& mDigit) {
    UInt_t time_offset = mDigit.GetTimeOffset();
    Bool_t err = ((   time_offset == RAWDPT_DVX_FMT0_TIMEOVERFLOW)
                  || (time_offset == RAWDPT_DVX_FMT0_TIMENOHIT)
                  || (time_offset == RAWDPT_DVX_FMT0_TIMEMAX) );  
    return !err;
}

class IP0dHealth: public COMET::ICOMETEventLoopFunction {
public:
  IP0dHealth() :fDB(COMET::IGeometryDatabase::Get()),
    fTFBs(300, -1),
    fReadEventCount(0),
          fBarTrig(0),
          fBarTrig300(0),
          hRawHits_x(0),
          hRawHits_y(0),
          hRawHits_z(0),
          hRawHits_xz(0),
          hRawHits_yz(0),
          fTrigger(0),
          fOutputFile(0){ }

  virtual ~IP0dHealth() {};

  void Initialize();
  void Finalize(ICOMETOutput * const output);
  using COMET::ICOMETEventLoopFunction::Finalize;
  
  virtual bool SetOption(std::string option,std::string value="") {
    if (option == "outputfile" ){
      fOutputFileName = value;
      return true;
    }
    return false;
  }

  bool operator () (COMET::ICOMETEvent& event);

private:
  COMET::IGeometryDatabase& fDB;
  std::vector<int> fTFBs;
  int fReadEventCount;
  bool fAllocated;
  

  //Timing Histograms
  TH1I* fBarTrig;
  TH1I* fBarTrig300;
  TH1I *hRawHits_x,*hRawHits_y,*hRawHits_z;
  TH2I *hRawHits_xz,*hRawHits_yz;

  //Histogram of trigger types
  TH1I* fTrigger;

  bool fStart;
  int fRunId;
  int fEventCounter;

  TFile* fOutputFile;
  std::string fOutputFileName;

};

void IP0dHealth::Initialize() {
  
  TH1::AddDirectory(false);

  fTrigger = new TH1I("Triggers","Triggers",256,0,256);

  //treePBS->Branch("startUnixTime",&pbs.startUnixTime,"startUnixTime/I:events:p0dEvents:gt300ADC:fidcut_gt300ADC");
  treePBS->Branch("startUnixTime",&pbs.startUnixTime,"startUnixTime/I");
  treePBS->Branch("events",&pbs.events,"events/I");
  treePBS->Branch("p0dEvents",&pbs.p0dEvents,"p0dEvents/I");
  treePBS->Branch("gt300ADC",&pbs.gt300ADC,"gt300ADC/I");
  treePBS->Branch("fidcut_gt300ADC",&pbs.fidcut_gt300ADC,"fidcut_gt300ADC/I");
  
  pbs.events=pbs.p0dEvents=pbs.gt300ADC=pbs.fidcut_gt300ADC=0;
  fEventCounter=0;
  fStart=true;
}

void IP0dHealth::Finalize(ICOMETOutput* const output){
  COMETLog("Data had " << fEventCounter << " events in it.");
  fOutputFile->cd();
  if ( fBarTrig ) {
    fBarTrig->Write();
  }
  if ( fBarTrig300 ) {
    fBarTrig300->Write();
  }
  if ( hRawHits_x ) {
    hRawHits_x->Write();
  }
  if ( hRawHits_y ) {
    hRawHits_y->Write();
  }
  if ( hRawHits_z ) {
    hRawHits_z->Write();
  }
  if ( hRawHits_xz ) {
    hRawHits_xz->Write();
  }
  if ( hRawHits_yz ) {
    hRawHits_yz->Write();
  }
  if ( fTrigger ) {
    fTrigger->Write();
  }

  pbs.events = fReadEventCount;
  treePBS->Fill();
  treePBS->Write();
  
  fOutputFile->Close();
}

bool IP0dHealth::operator () (COMET::ICOMETEvent& event) {

  unsigned int icurTime;
  int hitcnt300, hitcnt300fid;
  const COMET::IP0DChannelMap& P0DChanMap = fDB.P0D();

  hitcnt300=hitcnt300fid=0;
  fRunId = event.GetContext().GetRun();
  const int eventId = event.GetContext().GetEvent();
  
  if(eventId%50 == 0){
    COMETLog("\t Run " << fRunId
       << " event "<<eventId);
  }
  if(fStart==true && fRunId!=eventId){
    fOutputFile = new TFile( fOutputFileName != "" ?
           fOutputFileName.c_str() :
           Form("p0d-health-beam-%d.root", fRunId),
           "RECREATE");
    fStart=false;
  } 

  fTrigger->Fill(event.GetHeader().GetTriggerBits());
  if(event.GetHeader().GetTriggerBits()==1){ //beam events only
    
    if(fEventCounter==0){
      icurTime=event.GetContext().GetTimeStamp();
      pbs.startUnixTime=icurTime;
      
      fBarTrig300 = new TH1I(Form("300ADC Hit Timing-%d",icurTime),Form("300ADC Hit Timing-%d",icurTime),12000,-6000,6000);
      fBarTrig300->GetXaxis()->SetTitle("Time (2.5ns per bin)");
      fBarTrig300->GetYaxis()->SetTitle("N");
      
      fBarTrig = new TH1I(Form("Hit Timing-%d",icurTime),Form("Hit Timing-%d",icurTime),12000,-6000,6000);
      fBarTrig->GetXaxis()->SetTitle("Time (2.5ns per bin)");
      fBarTrig->GetYaxis()->SetTitle("N");
      
      //create 1d hists for x, y, z
      hRawHits_x = new TH1I(Form("Hit Dist x-%d",icurTime),Form("Hit Dist x-%d",icurTime),126,0,126);
      hRawHits_x->GetXaxis()->SetTitle("x layer");
      hRawHits_y = new TH1I(Form("Hit Dist y-%d",icurTime),Form("Hit Dist y-%d",icurTime),134,0,134);
      hRawHits_y->GetXaxis()->SetTitle("y layer");
      hRawHits_z = new TH1I(Form("Hit Dist z-%d",icurTime),Form("Hit Dist z-%d",icurTime),40,0,40);
      hRawHits_z->GetXaxis()->SetTitle("z p0dule");
      
      //create 2d hists for xz,yz
      hRawHits_xz = new TH2I(Form("Hit Dist xz-%d",icurTime),Form("Hit Dist xz-%d",icurTime),40,0,40,126,0,126);
      hRawHits_xz->GetYaxis()->SetTitle("x layer");
      hRawHits_xz->GetXaxis()->SetTitle("z p0dule");
      hRawHits_yz = new TH2I(Form("Hit Dist yz-%d",icurTime),Form("Hit Dist yz-%d",icurTime),40,0,40,134,0,134);
      hRawHits_yz->GetYaxis()->SetTitle("y layer");
      hRawHits_yz->GetXaxis()->SetTitle("z p0dule");
      
      fEventCounter=0;
    }
    
    fEventCounter += 1;

    IHandle<ICOMETRawEvent> rawEvt(event.GetRaw());
    assert(rawEvt);

    //const int first = 0;
    //const int last = 22;
    
    int trigtime = 0XDEADBEEF;
    //int timesincebusy_mcm,evttime;
    
    IHandle<IMCMBank> mcmInfo;
    while(mcmInfo = rawEvt->GetMidasBank<COMET::IMCMBank>("", mcmInfo)){
      std::stringstream rr;
      rr << mcmInfo->GetName();
      std::string name(rr.str());
      std::string p0d="P";
      if(name[0]!=p0d) continue;
      trigtime = mcmInfo->GetUnixTimeSSecTrig();
      //timesincebusy_mcm = mcmInfo->GetTimeSinceBusy();
      //evttime = mcmInfo->GetEventTime();
      
    }
    
    IHandle<ITripTDigitBank> triptBank;
    while (triptBank = rawEvt->GetMidasBank<ITripTDigitBank>("",triptBank) ) {
      if (triptBank->GetSubDetector() != 1){
        // Not a P0D Bank, ignore it
        COMETVerbose("Not a P0D Bank, ignore it");
        continue;
      }
      //reset event category counters
      hitcnt300=0;
      hitcnt300fid=0;
      
      IMidasTripTDigitItr itr(triptBank->GetMidasTripTDigitItr());
      while (!itr.EOD() ) {
        IMidasTripTDigit digit(itr.Get());
        ITFBChannelId tfbChannelId(digit.GetChannelId());
        //const int rmm = digit.GetRMMNum();
        //const int port = digit.GetTFBNum();
        //const int tripT = digit.GetTripTNum();
        //const int channel = digit.GetChannelNum(); //Unused
        //const int capacitor = digit.GetIntegrationNum();
        const int time_tfb = digit.GetTime();
        //const int time_offset = digit.GetTimeOffset(); //Unused
        const int charge = digit.GetHighGainADC();

        //ignore first and last integration cycle as they tend to shift the pedestal values (in error)
        //if (capacitor == first || capacitor == last) continue;
      
        //const int id = (rmm*48 + port)*64 + tripT*16 + channel;
        //const int id = tfbChannelId.GetCableId();
      
        COMET::IGeometryId geometryId;
        if (!fDB.GetGeometryId(geometryId, tfbChannelId)) {
          COMETVerbose("Skip: !GetGeometryId()");
          continue;
        }
      
        int p0dule, layer, bar;
      
        if(!P0DChanMap.GeomIDToLocation(p0dule, layer, bar, geometryId)) {
          COMETVerbose("Skip: !GeomIDToLocation()");
          continue;
        }
         
        //int position = p0dule*140+50+bar+layer*6000;
        int timeRelToTrig_tfb = (time_tfb)%400000000-trigtime*4;
        //Fill timing histograms
        if(charge>300){
          fBarTrig300->Fill(timeRelToTrig_tfb);
          hitcnt300++;
          //fiducial counts
          if(bar>2 && ((layer==0 && bar>3 && bar<122) || (layer=1 && bar>3 && bar<130))){
            hitcnt300fid++;
          }
          
          hRawHits_z->Fill(p0dule);
          if(layer==0){
            hRawHits_x->Fill(bar);
            hRawHits_xz->Fill(p0dule,bar);
          }else{
            hRawHits_y->Fill(bar);
            hRawHits_yz->Fill(p0dule,bar);
          }
        }
        fBarTrig->Fill(timeRelToTrig_tfb);
      }
    }
  }
  
  if(hitcnt300>15) pbs.gt300ADC++;
  if(hitcnt300fid>15) pbs.fidcut_gt300ADC++;
 
  ++fReadEventCount;
  return true;
  
  
}

int main(int argc, char **argv) {
  IP0dHealth userCode;
  cometEventLoop(argc,argv,userCode,1);
}