// UnCal.cc
// Contact person: Gabriel Orebi Gann <gorebigann@lbl.gov> - Author
// See UnCal.hh for more details
//———————————————————————//

#include <CLHEP/Random/RandGauss.h>

#include <G4UImessenger.hh>
#include <Randomize.hh>

#include <RAT/DS/PMT.hh>
#include <RAT/DS/MCEV.hh>
#include <RAT/UnCal.hh>
#include <RAT/DB.hh>
#include <RAT/Log.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>

#include <iostream>
#include <sstream>
#include <cmath>
#include <algorithm>
#include <vector>
using namespace std;

namespace RAT {

UnCal::UnCal() {

}

void UnCal::BeginOfRun(DS::Run &run) {
  fNcrate = 19;
  fNcard = 16;
  fNchannel = 32;
  fNcell = 16;
  fNTotCell = fNcrate * fNcard *fNchannel * fNcell;
  fTotChannel = fNcrate * fNcard *fNchannel;

  fUnCalbank = DB::Get()->GetLink("DAQ");
  int UnCalType = fUnCalbank->GetI("uncal_type");
  int ECAType = fUnCalbank->GetI("eca_type");
  int ECAPDSTpattern = fUnCalbank->GetI("eca_pdst_pattern");
  int ECAPDSTrate = fUnCalbank->GetI("eca_pdst_rate");
  int ECATSLPpattern = fUnCalbank->GetI("eca_tslp_pattern");
  int ECATSLPrate = fUnCalbank->GetI("eca_tslp_rate");
  //int PCAType = fUnCalbank->GetI("pca_type");
  UnCal::UnCalibStatus(UnCalType, fDoECA, fDoPCA, fDoMC);
  bool mc;
  UnCal::UnCalibStatus(ECAType, fDoPed, fDoTslp, mc);
  fDoUncal = true;
  if(!fDoECA && !fDoPCA && !fDoMC)fDoUncal = false;


  //Get the pattern and rate
  std::stringstream pdst_pattern_rate, tslp_pattern_rate;
  pdst_pattern_rate << ECAPDSTpattern <<"_"<<ECAPDSTrate;
  tslp_pattern_rate << ECATSLPpattern <<"_"<<ECATSLPrate;

  info << " (ECACal::BeginOfRun) Using pattern_rate " << pdst_pattern_rate.str() << " for PDST " <<newline;
  info << " (ECACal::BeginOfRun) Using pattern_rate " << tslp_pattern_rate.str() << " for TSLP " <<newline;
  fPedsbank = DB::Get()->GetLink("ECA_PDST", pdst_pattern_rate.str());
  fTslpbank = DB::Get()->GetLink("ECA_TSLP", tslp_pattern_rate.str());

  if(fDoPed){
    try{
      fPedsbank->GetS("index");
    }
    catch(DBNotFoundError){
      UnCal::Die("",5,ECAPDSTpattern,ECAPDSTrate);
    }
    fPedQHS = fPedsbank->GetFArrayFromD("pdst_qhs");
    fPedQHL = fPedsbank->GetFArrayFromD("pdst_qhl");
    fPedQLX = fPedsbank->GetFArrayFromD("pdst_qlx");
    //fPedTAC = fPedsbank->GetDArray("pdst_tac");
    fPedSigQHS = fPedsbank->GetFArrayFromD("pdst_qhswidth");
    fPedSigQHL = fPedsbank->GetFArrayFromD("pdst_qhlwidth");
    fPedSigQLX = fPedsbank->GetFArrayFromD("pdst_qlxwidth");
    fPedSigTAC = fPedsbank->GetFArrayFromD("pdst_tacwidth");
  }

  if(fDoTslp){
    try{
      fTslpbank->GetS("index");
    }
    catch(DBNotFoundError){
      UnCal::Die("",5,ECAPDSTpattern,ECAPDSTrate);
    }
    fTslope = fTslpbank->GetFArrayFromD("tslp_points");

    fNSlopePoints = fTslpbank->GetI("tslp_ntimes");
    fNdim = static_cast<int>(fNSlopePoints / 32)+1;

    fTslopeX = fTslpbank->GetFArrayFromD("tslp_times");

    fTslopeBadTemp = fTslpbank->GetIArray("tslp_bad");
    // Sort status words by cell
    fTslopeBad.resize(fNTotCell);
    if((int)fTslopeBadTemp.size() != fNTotCell*fNdim){
      warn<<"=================================================="<<newline;
      warn<<"UnCal Proc: PROBLEM WITH number of TSlope points (a)"<<newline;
      warn<<fTslopeBadTemp.size()<<" status words in ratdb"<<newline;
      warn<<"But we expect "<<fNTotCell<<" cells * "<<fNdim<<" words ";
      warn<<" = "<<fNTotCell*fNdim<<" words"<<newline;
      warn<<"=================================================="<<newline;
    }
    int counter = 0;
    for(int n=0; n<fNTotCell; ++n){
      fTslopeBad[n].resize(fNdim);
      for(int j=0; j<fNdim; ++j){
        fTslopeBad[n][j] = fTslopeBadTemp[counter];
	++counter;
      }
    }
    if(counter != fNTotCell*fNdim){
      warn<<"=================================================="<<newline;
      warn<<"UnCal Proc: PROBLEM WITH number of TSlope points (b)"<<newline;
      warn<<counter<<" status words read in"<<newline;
      warn<<"But we expect "<<fNTotCell<<" cells * "<<fNdim<<" words ";
      warn<<" = "<<fNTotCell*fNdim<<" words"<<newline;
      warn<<"=================================================="<<newline;
    }
  }

  fQHS = fQHL = fQLX = 0.0;
  fTAC = 0.0;
  fsQHS = fsQHL = fsQLX = fsTAC = 0;
  fBits = new BitManip();

}

void UnCal::EndOfRun(DS::Run &run) {
  (void)run;
}

UnCal::~UnCal()
{
}

void UnCal::UnCalibStatus(int flag, bool &eca, bool &pca, bool &mc)
{
  eca = false;
  pca = false;
  mc = false;

  if(flag==9){mc = true;}
  else{
    if(flag%10==1)eca = true;
    if(int((flag%100)/10)==1)pca = true;
  }

  return;
}

void UnCal::UnCalibrateEvent(DS::Entry& ds, int iev, double GTtime)
{
  const DU::PMTInfo& pmtList = DU::Utility::Get()->GetPMTInfo();
  DS::EV& ev = ds.GetEV(iev);
  DS::MCEV& mcev = ds.GetMCEV(iev);
  size_t npmttrue = mcev.GetMCHits().GetAllCount();
  size_t nhit = ev.GetUncalPMTs().GetCount(); // NOTE normal only
  size_t npmtuncal = ev.GetUncalPMTs().GetAllCount();
  ULong64_t totalQ = 0;

  // If no uncalibration requested, and no PMTUncal to be digitised
  if(!fDoUncal && (npmtuncal==0))return;

  // For a regular run we expect:
  // Data: SHOULD NOT BE UNCALIBRATING DATA!
  // MC:   NHit == 0 until we create some PMTUncal
  //      NPMTUncal == 0 since we haven't uncalibrated yet

  // Die, if the numbers of Cal / UnCal PMTs / NHit are not as expected
  if(nhit != 0)UnCal::Die("MC",1,nhit,npmttrue);
  if(npmtuncal!=0)UnCal::Die("MC",2,npmtuncal,0);

  for(size_t ih = 0; ih<npmttrue; ++ih){
    // Get simulated values
    const DS::MCHit &truth = mcev.GetMCHits().GetAllPMT(ih);
    fQHS = truth.GetQHS();
    fQHL = truth.GetQHL();
    fQLX = truth.GetQLX();
    fTAC = GTtime - truth.GetTime();
    if(fTAC<0)UnCal::Die("",4,(int)truth.GetTime(),(int)GTtime);
    const unsigned int cell = truth.GetCellID();
    const unsigned int lcn = truth.GetID();
    const unsigned int cccc = truth.GetCrateCardChannelCell();
    // Uncalibrate in requested ways
    if(fDoECA){
      // Convert tac to ADC counts before adding the ped
      if(fDoTslp)UnCal::TSlope(cccc);
      //if(!(temp>200.&& temp<3700.))G4cout<<"TAC = "<<temp<<" from ns val "<<tac<<G4endl;

      // Add pedestals + smearing
      if(fDoPed)UnCal::Pedestals(cccc);
    }

    // Digitise
    fsQHS = UnCal::DigitiseQ(fQHS);
    fsQHL = UnCal::DigitiseQ(fQHL);
    fsQLX = UnCal::DigitiseQ(fQLX);
    fsTAC = UnCal::DigitiseT(fTAC);
    totalQ += fsQHS;

    // Create new PMTUncal and fill it
    DS::PMTUncal pmt;
    pmt.SetPMTParameters(lcn,fsQHS,fsQHL,fsQLX,fsTAC,cell,truth.GetChanFlags());
    ev.GetUncalPMTs().AddPMT( pmt, pmtList.GetType( lcn ) );
  }
  ev.SetNhits( ev.GetUncalPMTs().GetCount() );
  ev.SetTotalCharge(totalQ);

  return;

}

// Add pedestals, plus intrinsic pedestal smearing
void UnCal::Pedestals(int cccc)
{
  // Add pedestal median, plus smearing, to Q values
  fQHS += fPedQHS[cccc] + fPedSigQHS[cccc]*CLHEP::RandGauss::shoot();
  fQHL += fPedQHL[cccc] + fPedSigQHL[cccc]*CLHEP::RandGauss::shoot();
  fQLX += fPedQLX[cccc] + fPedSigQLX[cccc]*CLHEP::RandGauss::shoot();
  // Add pedestal smearing to T value
  // IF tac=0 already, it was a bad TSlope UnCal, so leave at 0
  double throwme = fPedSigTAC[cccc]*CLHEP::RandGauss::shoot();
  if(fTAC>0) fTAC += throwme;
}

// Convert from ns back to ADC counts
void UnCal::TSlope(int cccc)
{
  // Assume time reversal already handled in the event builder
  double result = 0.;

  // Get the TAC slope points for this cell
  // And find where in time the hit fell
  // ONLY USE GOOD POINTS ON TAC SLOPE (i.e. not bad/susp)
  // but using the curl region is fine
  std::vector<int> bad = fTslopeBad[cccc];
  float *adc = new float[fNSlopePoints];
  int ilow = -1;
  int ihigh = 99;
  //bool status = 0;
  int lastgood = 0;
  int penultimate = 0;

  for(int i=0;i<fNSlopePoints;++i){
    int arrayentry = cccc*fNSlopePoints + i;
    adc[i] = fTslope[arrayentry];

    // only want good tac points
    if(fBits->TestBit(bad,i))continue;

    penultimate = lastgood;
    lastgood = i;

    // bracket value
    if((fTAC>=fTslopeX[i])&&(i>ilow))ilow = i;
    if((fTAC<=fTslopeX[i])&&(i<ihigh))ihigh=i;//status=1;}
  }

  // If the hit was < time[0], it was too late, lose it
  // NO: temp fix, if: 0<hitT<time[0] then extrap from points 0,1
  if(ilow<0){
    //UnCal::DontDie("late",1,time,fTslopeX[0],fTslopeX[fNSlopePoints-1]);
    //result = 0.;
    //return result;
    // Set brackets to the first 2 points and extrap from that slope
    ilow = 0;
    ihigh = 1;
  }

  // If the hit was >time[fNSlopePoints-1], it was too early(!), lose it
  // (Latter should never happen, since max TSlopeX > GTValid, unless the points were bad)
  // NO: might happen for bad points so again, extrapolate based on last two good points
  if(ihigh>=fNSlopePoints){
//  UnCal::DontDie("early",1,(int)fTAC,(int)fTslopeX[0],(int)fTslopeX[fNSlopePoints-1]);
//  result = 0;
//  fTAC = result;
//  delete[] adc;
//  return;
    ihigh = lastgood;
    ilow = penultimate;
    debug<<"UnCalib: "<<cccc<<" cell at "<<fTAC<<" after "<<lastgood<<" so using points "<<ilow<<" and "<<ihigh<<newline;
  }

  if(ilow==ihigh){ // on a TAC point
    result = adc[ilow];
    fTAC = result;
    delete[] adc;
    return;
  }

  // Otherwise, interpolate as expected between points [ilow, ihigh]
  double grad = (adc[ihigh] - adc[ilow]) / (fTslopeX[ihigh] - fTslopeX[ilow]);
  double icpt = adc[ihigh] - grad*fTslopeX[ihigh];
  result = grad*fTAC + icpt;
  fTAC = result;
  delete[] adc;
}

unsigned short UnCal::DigitiseT(double RelHitTime)
{
   //eventually this will have to return a 12-bit integer,
   //in ADC or TDC counts
   unsigned short time = (unsigned short) round(RelHitTime);
   if(time>4095)time=4095;
   //if(!(time>200.&& time<3700.))G4cout<<"TAC = "<<time<<" (ADC) from "<<RelHitTime<<G4endl;
   return time;
}

unsigned short UnCal::DigitiseQ(float HitCharge)
{
   //Eventually this will have to return a 12-bit integer,
   //in ADC counts
   unsigned short charge = (unsigned short) round(HitCharge);
   if(charge>4095)charge=4095;
   return charge;
}

bool UnCal::FileExists(std::string filename)
{
  struct stat stFileInfo;
  bool open;
  int intStat;

  // Attempt to get the file attributes
  intStat = stat(filename.c_str(),&stFileInfo);
  if(intStat == 0) {
    open = true;
  } else {
    open = false;
  }

  return(open);
}

// Handle errors
void UnCal::DontDie(std::string message, int flag_err, int info1, int info2, int info3)
{
  warn << YELLOW << " UnCal SETBACK:" << CLR << newline;

  if(flag_err==1){ // PMT hit is outside TAC region
    warn<<"This PMT, hit at "<<info1<<"ns before the GT, was too "<<message<<" in time."<<newline;
    warn<<"The TAC slope measurements extend from "<<info2<<" to "<<info3<<newline;
    warn<<"Something is wrong, setting hit time to 0."<<newline;
    warn<<"(Don't worry about this message at the moment: the PMT hit times need fixing anyway)"<<newline;
    warn<<newline;
  }

  if(flag_err==3){ // PMT hit is outside TAC region
    warn<<"This PMT hit occurred at "<<info1<<newline;
    warn<<"But the TAC slope measurements only extend to "<<info2<<newline;
    warn<<"Something is wrong, setting hit time to 0."<<newline;
    warn<<"(Don't worry about this message at the moment: the PMT hit times need fixing anyway)"<<newline;
    warn<<newline;
  }

}

// Handle FATAL errors
void UnCal::Die(std::string message, int flag_err, int info1, int info2, int return_code)
{
  warn << RED << " UnCal PROBLEM:" << CLR << newline;
  if(message==""){/**/}

  if(flag_err==1){ // NHit != 0
    warn<<"Discrepancy in # events:"<<newline;
    warn<<"We have "<<info1<<" NHit but "<<info2<<" PMTTruth in this event"<<newline;
    warn<<"Something is wrong."<<newline;
    warn<<newline;
  }

  if(flag_err==2){ // NPMTUncal != 0
    warn<<"Discrepancy in # events:"<<newline;
    warn<<"We have "<<info1<<" PMTUncal in this event"<<newline;
    warn<<"But we have not uncalibrated yet..."<<newline;
    warn<<"Something is wrong."<<newline;
    warn<<newline;
  }

  if(flag_err==3){ // PMT hit is outside TAC region
    warn<<"This PMT hit occurred at "<<info1<<newline;
    warn<<"But the TAC slope measurements only extend to "<<info2<<newline;
    warn<<"Something is wrong."<<newline;
    warn<<newline;
  }

  if(flag_err==4){ // Hit was after GT!!
    warn<<"This PMT hit occurred at "<<info1<<newline;
    warn<<"But the GT was at "<<info2<<newline;
    warn<<"Something is wrong."<<newline;
    warn<<newline;
  }
  if(flag_err==5){
    warn<<"Requested pattern ("<<info1<<") or rate ("<<info2<<") do not exists in the DB."<<newline;
  }

  exit(return_code);
}


} // namespace RAT