#include "fQEventManager.h"

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <string>
#include <bitset>
#include <stdio.h>

#include <TFile.h>
#include <TH1F.h>

#include "fiTQun_shared.h"

#ifndef NOSKLIBRARIES

#include "skheadC.h"
#include "skparmC.h"
#include "skbadcC.h"
#include "sktqC.h"

#else

#include "WCSimWrap.h"

#endif

// In g77, "real*4" functions return double
//#ifdef FQ_G77
#if defined(f2cFortran)&&!defined(gFortran)
#define FORTRETURN_FLOAT double
#else
#define FORTRETURN_FLOAT float
#endif

extern"C" {
#ifndef NOSKLIBRARIES
  void fortinit_(char*, char*, int, int);
  int fortread_();
  void skclosef_(int*);
  
  void skread_(int *);
  void set_timing_gate_nsec2_(float *, float *);
  
  void kzmnum_(int *, int *);
  void skrdsub_(int *, int *, int *);
  FORTRETURN_FLOAT diftevsk_(int *);
  
  void savetqbanks_();
  void restoretqbanks_();
  
  int getfakepi0_(int *, int *, int *);
  
#endif
}

fQEventManager* fQEventManager::staticthis = NULL;

char fQEventManager::ofname[1024];

fQEventManager::fQEventManager(){
  
  fQuiet=1;
  
  flgHist=0;// Don't output tHit histograms
//  flgHist=1;
  
  if (!(fiTQun_shared::Get()->AreYouReady())) {
    std::cout << "fQEventManager::fQEventManager: fiTQun_shared not initialized!" << std::endl;
    exit(-1);
  }
  
  flgData = false;
  
  PMTGainCorr = 1.;
  WAttenLCorr = 1.;
  
#ifndef NOSKLIBRARIES
  SKConsts_tvar = new SKTVarConsts(fiTQun_shared::Get()->GetSKVer());
#endif
  
  fNPMT = fiTQun_shared::Get()->MAXNPMT();
  
  PMT_iHit = new int[fNPMT];
  PMT_q = new double[fNPMT];
  PMT_t = new double[fNPMT];
  
  PMTall_nhit=0;
  PMTall_icab = new int[30*fNPMT];
  PMTall_q = new double[30*fNPMT];
  PMTall_t = new double[30*fNPMT];
  MinMax_PMTall_t[0]=0.;
  MinMax_PMTall_t[1]=0.;
  
  nEvt=0;
}

fQEventManager::~fQEventManager(){
  
  delete[] PMT_iHit;
  delete[] PMT_q;
  delete[] PMT_t;
  
  delete[] PMTall_icab;
  delete[] PMTall_q;
  delete[] PMTall_t;
  
  std::cout << "Destructing fQEventManager!!" << std::endl;
}


int fQEventManager::SetGate(float t_Start, float t_End, bool flgShiftGate, bool flgMaskDeadTimePMT){
  
  for (int icab=0; icab<fNPMT; icab++) {
    PMT_iHit[icab]=0;
    PMT_q[icab]=0.;
    PMT_t[icab]=0.;
    if (fiTQun_shared::Get()->PMTtype[icab]<0) PMT_iHit[icab]=-1;//ignore missing PMTs
  }
  
  if (flgMaskDeadTimePMT) {
    std::cout << "fQEventManager::SetGate: Masking dead-time PMTs!" << std::endl;
    for (int ih=0; ih<PMTall_nhit; ih++) {
      int icab = PMTall_icab[ih];
      if (PMTall_t[ih]>t_Start-1000. && PMTall_t[ih]<t_Start) PMT_iHit[icab]=-1;// Mask dead-time PMT
    }// put A/B channel masking for ATM?
  }
  
  int nhitPMT=0;
  
#ifndef NOSKLIBRARIES
  
  if (fiTQun_shared::Get()->GetSKVer()<=3) {// ATM
    
    for (int icab=0; icab<fNPMT; icab++) {
      if (combad_.ibad[icab]!=0) PMT_iHit[icab]=-1;//ignore bad channels
    }
    
    for (int ih=0; ih<PMTall_nhit; ih++) {
      if (!(PMTall_t[ih]>=t_Start && PMTall_t[ih]<t_End)) continue;
      // In window!
      int icab = PMTall_icab[ih];
      
      if (PMT_iHit[icab]==0) nhitPMT++;
      else if (PMT_iHit[icab]==1) {
        if (PMT_t[icab]<PMTall_t[ih]) continue;// take earlier hit within time window!
      }
      else continue;
      
      PMT_q[icab]=PMTall_q[ih];
      PMT_t[icab]=PMTall_t[ih];
      
      PMT_iHit[icab]=1;//hit!
      
    }
    
  }
  else {// QBEE
    double tOfst=0.;
    if (flgShiftGate) {
      savetqbanks_();// store TQREAL banks before modification
      set_timing_gate_nsec2_(&t_Start,&t_End);
      
      int LUN=-10;
      skread_(&LUN);
      
      tOfst=t_Start;
    }
    
    /// charge and time info
    for (int icab=0; icab<fNPMT; icab++) {
      PMT_q[icab]=skq_.qisk[icab];
      PMT_t[icab]=skt_.tisk[icab]+tOfst;// time is in the original frame, before shift.
      if (combad_.ibad[icab]!=0) PMT_iHit[icab]=-1;//ignore bad channels
    }
    
    /// following fills PMT hit information
    for (int ih=0; ih<skq_.nqisk; ih++) {
      int icab = skchnl_.ihcab[ih] - 1;//cable # in C++
      if (!(PMT_t[icab]>=t_Start && PMT_t[icab]<t_End)) continue;
      if (PMT_iHit[icab]!=0) continue;
      PMT_iHit[icab]=1;//hit!
      nhitPMT++;
    }
    
    if (flgShiftGate) {
      // Shift TQ(A)REAL time window back to original for output to ZBS bank
      restoretqbanks_();
    }
  }
  
#else
  
  // fill charge and time info from wcsim root file
  WCSimWrap * wc  = WCSimWrap::Get();
  
  // now fill with hits info
  int nsubev = wc->NSubEvt();
  if (nsubev==0)nsubev=1;
  for (int isub=0; isub<nsubev; isub++){
    WCSimRootTrigger * evt = wc->SubEvt(isub);
    double aSubToffs = wc->GetTOffset(isub);
    for (int ich = 0; ich<evt->GetNcherenkovdigihits(); ich++){
      WCSimRootCherenkovDigiHit *cDigiHit = dynamic_cast<WCSimRootCherenkovDigiHit*>( evt->GetCherenkovDigiHits()->At(ich) );
      int icab = cDigiHit->GetTubeId() - 1; //convert to C++ index
      
      double tmp_t = cDigiHit->GetT()+aSubToffs;
      
      if (!(tmp_t>=t_Start && tmp_t<t_End)) continue;
      if (PMT_iHit[icab]==-1) continue;
      
      PMT_q[icab] = cDigiHit->GetQ();
      PMT_t[icab] = tmp_t;
      
      PMT_iHit[icab]=1;//hit!
      nhitPMT++;
    }
  }
#endif
  
  return nhitPMT;
  
}

void fQEventManager::LoadEvent() {
  
  flgData = false;
  
  PMTGainCorr = 1.;
  WAttenLCorr = 1.;
  
#ifndef NOSKLIBRARIES
  
  int iRun=skhead_.nrunsk;//Run number
  int iSub=skhead_.nsubsk;//Run number
  int iEv=skhead_.nevsk;//Run number
  
  int fpi0nrun, fpi0nsub, fpi0nev;
  int iret = getfakepi0_(&fpi0nrun,&fpi0nsub,&fpi0nev);
  if (iret) {// contains FAKEPI0 bank
    std::cout << "Detected FAKEPI0 bank!" << std::endl;
    iRun = fpi0nrun;
    iSub = fpi0nsub;
    iEv  = fpi0nev;
  }
  
  std::cout << Form("Loading event: Run %d, Sub %d, Event %d",iRun,iSub,iEv) << std::endl;
  
  if (skhead_.mdrnsk!=0 && skhead_.mdrnsk!=999999) {// Data
    
    flgData = true;
    
    std::cout << "  Obtaining detector parametors for Run" << iRun << std::endl;
    
    PMTGainCorr = SKConsts_tvar->GetGainFact(iRun);
    std::cout << "      PMT gain correction factor: " << PMTGainCorr << std::endl;
    
    WAttenLCorr = SKConsts_tvar->GetAttenLFact(iRun);
    std::cout << "      Attenuation length correction factor: " << WAttenLCorr << std::endl;
    
  }
  else flgData = false;// MC
  
#else
  
#endif
  
  GetAllHits();
  
  nEvt++;
}

void fQEventManager::GetAllHits() {
  
  TFile *fHist=NULL;
  TH1F *h_thit_all, *h_thit_sub;
  
  if (flgHist) {
    
    if (nEvt>0) {
      fHist = new TFile(Form("%s_tHit.root",ofname),"UPDATE");
    }
    else {
      fHist = new TFile(Form("%s_tHit.root",ofname),"RECREATE");
    }
    
    h_thit_all = new TH1F(Form("h_thit_%d_all",nEvt),Form("Raw hit time, iEvt=%d",nEvt),10000,-50000,50000);
    
  }
  
  PMTall_nhit=0;
  
  n_Gate=0;
  
  MinMax_PMTall_t[0]=1e300;
  MinMax_PMTall_t[1]=-1e300;
  
#ifndef NOSKLIBRARIES
  
  if (fiTQun_shared::Get()->GetSKVer()<=3) {// ATM
    
    int parent_48bitClock[3];
    
    int LUN=10;
    int ierr;
    
    int nsub_ATM,isub_ATM;
    kzmnum_(&nsub_ATM,&isub_ATM);// Get total number of subevents
    if (isub_ATM!=1) std::cout << Form("fQEventManager::GetAllHits: isub_ATM!=1, isub=%d",isub_ATM) << std::endl;
    
    std::cout << Form("ATM: Reading %d subevents...",nsub_ATM) << std::endl;
    for (isub_ATM=0; isub_ATM<nsub_ATM; isub_ATM++) {
      
      skrdsub_(&LUN,&isub_ATM,&ierr);
      
      if (isub_ATM==0) {
        // Save the time of the parent event
        for (int i=0; i<3; i++) parent_48bitClock[i] = skhead_.nt48sk[i];
      }
      
      double tdiff_sub = diftevsk_(parent_48bitClock);// In ns
      std::cout << Form(" Subevent %d : tdiff = %.0fns",isub_ATM,tdiff_sub);
      
      if (flgHist) {
        h_thit_sub = new TH1F(Form("h_thit_%d_sub%d",nEvt,isub_ATM),Form("Raw hit time, iEvt=%d, iSub=%d",nEvt,isub_ATM),10000,-50000,50000);
      }
      
      if (n_Gate>=100) {
        std::cout << "n_Gate>=100!" << std::endl;
        continue;
//        exit(-1);
      }
      
      t_Gate[n_Gate][0]=1e100;
      t_Gate[n_Gate][1]=-1e100;
      for (int ih=0; ih<skq_.nqisk; ih++) {
        int icab = skchnl_.ihcab[ih] - 1;//cable # in C++
        
        double q_tmp=skq_.qisk[icab];
        double t_tmp=skt_.tisk[icab]+tdiff_sub;
        
        // Ignore hits that occur prior to the HE software trigger
        if (t_tmp<800.) continue;
        
        if (!(t_tmp<16000.)) continue; // make it consistent with SK4
        
//        if (!(q_tmp>qThresh)) continue;
        if (fiTQun_shared::Get()->PMTtype[icab]==-1) continue;//ignore missing PMTs
        if (combad_.ibad[icab]!=0) continue;//ignore bad channels
        
        PMTall_icab[PMTall_nhit]=icab;
        PMTall_q[PMTall_nhit]=q_tmp;
        PMTall_t[PMTall_nhit]=t_tmp;
        
        if (flgHist) h_thit_sub->Fill(PMTall_t[PMTall_nhit]);
        
        if (PMTall_t[PMTall_nhit]<MinMax_PMTall_t[0]) MinMax_PMTall_t[0]=PMTall_t[PMTall_nhit];
        if (PMTall_t[PMTall_nhit]>MinMax_PMTall_t[1]) MinMax_PMTall_t[1]=PMTall_t[PMTall_nhit];
        
        if (PMTall_t[PMTall_nhit]<t_Gate[n_Gate][0]) t_Gate[n_Gate][0]=PMTall_t[PMTall_nhit];
        if (PMTall_t[PMTall_nhit]>t_Gate[n_Gate][1]) t_Gate[n_Gate][1]=PMTall_t[PMTall_nhit];
        
        PMTall_nhit++;
        
      }
      
      
      if (n_Gate>0) {
        if (!(t_Gate[n_Gate][0]>t_Gate[n_Gate-1][1])) {
          t_Gate[n_Gate][0]=t_Gate[n_Gate-1][1];
        }
      }
      
      std::cout << Form(", Gate=[%.1f,%.1f]",t_Gate[n_Gate][0],t_Gate[n_Gate][1]);
      
      if (t_Gate[n_Gate][0]<t_Gate[n_Gate][1]) {// gate properly defined!
        n_Gate++;
      }
      else {
        std::cout << " - skipping!";
      }
      std::cout << std::endl;
      
      if (flgHist) {
        h_thit_sub->Write();
      }
    }
    
    std::cout << Form("n_Gate=%d",n_Gate) << std::endl;
    
    isub_ATM=0;
    skrdsub_(&LUN,&isub_ATM,&ierr);//restore
    
  }
  else {// QBEE
    
    for (int i=0; i<sktqz_.nqiskz; i++) {
      
      // Ignore hits that occur prior to the HE software trigger
      if (sktqz_.tiskz[i]<800.) continue;
      
      if (!(sktqz_.tiskz[i]<16000.)) continue; // dark noise is only simulated ~+17us in NEUTRON skdetsim
      
      std::bitset<64>  bihtiflz(sktqz_.ihtiflz[i]);
      if (!bihtiflz.test(1)) continue;// Only consider hits within 40us gate
      
      int icab = sktqz_.icabiz[i]-1;// C++ indexing
      if (!(icab<fNPMT)) continue;// ID PMT only
      
//      if (!(sktqz_.qiskz[i]>qThresh)) continue;
      if (fiTQun_shared::Get()->PMTtype[icab]==-1) continue;//ignore missing PMTs
      if (combad_.ibad[icab]!=0) continue;//ignore bad channels
      
      PMTall_icab[PMTall_nhit]=icab;
      PMTall_q[PMTall_nhit]=sktqz_.qiskz[i];
      PMTall_t[PMTall_nhit]=sktqz_.tiskz[i];
      
      if (PMTall_t[PMTall_nhit]<MinMax_PMTall_t[0]) MinMax_PMTall_t[0]=PMTall_t[PMTall_nhit];
      if (PMTall_t[PMTall_nhit]>MinMax_PMTall_t[1]) MinMax_PMTall_t[1]=PMTall_t[PMTall_nhit];
      
      PMTall_nhit++;
    }
  }
  
#else
  WCSimWrap * wc  = WCSimWrap::Get(); // get wcsim information
  
  if (!fQuiet) std::cout<<"<fQEventManager::GetAllHits> nsubev="<<wc->NSubEvt()<<std::endl;
  for (int isubev=0; isubev<wc->NSubEvt(); isubev++){
    WCSimRootTrigger * evt = wc->SubEvt(isubev);
    double aSubToffs = wc->GetTOffset(isubev);
    if (!fQuiet) std::cout<<"<fQEventManager::GetAllHits> SubEv"<<isubev<<" nDigits="<<evt->GetNcherenkovdigihits()<<std::endl;
//    if (!fQuiet) std::cout<<"<fQEventManager::GetAllHits> qThresh="<<qThresh<<std::endl;
    for (int i=0; i<evt->GetNcherenkovdigihits(); i++) {
      WCSimRootCherenkovDigiHit *cDigiHit = dynamic_cast<WCSimRootCherenkovDigiHit*> ( evt->GetCherenkovDigiHits()->At(i) );
      if (!fQuiet) std::cout<<"<fQEventManager::GetAllHits> digit"<<i
        <<" tube="<<cDigiHit->GetTubeId()
        <<" T, Q="<<cDigiHit->GetT()+aSubToffs<<" ,"<<cDigiHit->GetQ()<<std::endl;
      // Ignore hits that occur prior to the HE software trigger
//      if ( cDigiHit->GetT()+aSubToffs<800.){
//        if (!fQuiet) std::cout<<"<fQEventManager::GetAllHits> WARNING! T-offs="
//          <<cDigiHit->GetT()+aSubToffs<<std::endl;
//        continue;
//      }
      
      int icab = cDigiHit->GetTubeId()-1;//convert to c++ index
      if (icab>=fNPMT){//ID PMT
        if (!fQuiet) std::cout<<"<fQEventManager::GetAllHits> WARNING! icab>fNPMT ="
          <<icab<<std::endl;
        continue;
      }
      
//      if (!(cDigiHit->GetQ()>qThresh)){
//        if (!fQuiet) std::cout<<"<fQEventManager::GetAllHits> WARNING! Q<=qThresh  ="
//          <<cDigiHit->GetQ()<<std::endl;
//        continue;
//      }
      if (fiTQun_shared::Get()->PMTtype[icab]==-1){
        if (!fQuiet) std::cout<<"<fQEventManager::GetAllHits> WARNING! PMTtype["<<icab
          <<"] = "<<fiTQun_shared::Get()->PMTtype[icab]
          <<std::endl;
        continue;//ignore missing PMTs
      }
      
      PMTall_icab[PMTall_nhit]=icab;
      PMTall_q[PMTall_nhit]=cDigiHit->GetQ();
      PMTall_t[PMTall_nhit]=cDigiHit->GetT()+aSubToffs;
      
      if (PMTall_t[PMTall_nhit]<MinMax_PMTall_t[0]) MinMax_PMTall_t[0]=PMTall_t[PMTall_nhit];
      if (PMTall_t[PMTall_nhit]>MinMax_PMTall_t[1]) MinMax_PMTall_t[1]=PMTall_t[PMTall_nhit];
      
      PMTall_nhit++;
    }
  }
#endif
  
  if (flgHist) {
    
    for (int ih=0; ih<PMTall_nhit; ih++) {
      h_thit_all->Fill(PMTall_t[ih]);
    }
    
    h_thit_all->Write();
    
    delete fHist;
    
  }
  
}