// Author: Andrew M. Rogers, NSCL 07/01/2006
// Original Design: Daniel Bazin, NSCL
//* Copyright (C) 2006-2008 Andrew M. Rogers

#include "TS800Crdc.h"
#include "TS800.h"
#include "nrutil.h"

#include <fstream>
#include <iostream>
#include <math.h>

#include <TDirectory.h>
#include <TFile.h>
#include <TNamed.h>
#include <TList.h>

ClassImp(TS800Crdc);

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// TS800Crdc                                                            //
//                                                                      //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

extern short gNumrecError, gNumrecMessage;


//______________________________________________________________________________
TS800Crdc::TS800Crdc(const TS800Crdc &crdc) : TNamed(crdc)
{
  // -- Copy Constructor.

  ((TS800Crdc&)crdc).Copy(*this);
}


//______________________________________________________________________________
void TS800Crdc::InitFix()
{
  pad.parent  = this;
  calc.parent = this;
  calc.pad    = &pad;
  channel     = new Short_t[4096];
  sample      = new Short_t[4096];
  data        = new Short_t[4096];
}


//______________________________________________________________________________
Int_t TS800Crdc::Calculate(Long64_t entry) 
{
  // -- Calculate an entry.
  //

  if(!fChain){printf("TTree or TChain has not been initialized!!!\n"); return -1;}

  Long64_t ibad  = 0;
  Double_t yInt  = sqrt(-1.0);
  short restore  = gNumrecMessage;
  gNumrecMessage = 0;

  // Notify of any corrections being applied.
  if(entry==0){
  if(calc.fInterpPoly==kTRUE)        { printf("Using POLYNOMIAL INTERPOLATION for bad pads.\n");}
  else if(calc.fInterpSpline3==kTRUE){ printf("Using CUBIC SPLINE INTERPOLATION for bad pads.\n");}
  else if(calc.fSatFlat==kTRUE)      { printf("Using FLATTENING for saturated pads.\n");}
  else if(fYSplineCorr==kTRUE)       { printf("Using Y Correcting Cubic Spline.\n");
    if(!strcmp(fYSplineMatchPnt,"begin") || !strcmp(fYSplineMatchPnt,"end")){
      printf("* Invalid CRDC Y Spline Matching Point.\n");
    }
  }
  else{                                printf("NO CRDC PAD CORRECTIONS APPLIED.\n");}
  }

  Clear();
  fCon->b_filled->GetEntry(entry);
  if (fCon->filled > 0) {
    // Choose the method to use for the pad calculation.
    if (calc.method == 1) {
      calc.CalculateGravity();
      //      x = x_offset[0] + x_slope[0] * calc.x_gravity;
      x = x_offset[1] 
	+ x_slope[1]*x_offset[0] 
	+ x_slope[0]*x_slope[1]*calc.x_gravity;
      if (isnan(calc.x_gravity)) ibad++;
    }
    if (calc.method == 2) {
      calc.CalculateFit();
      //      x = x_offset[0] + x_slope[0] * calc.x_fit;
      x = x_offset[1] 
	+ x_slope[1]*x_offset[0] 
	+ x_slope[0]*x_slope[1]*calc.x_fit;
      if (isnan(calc.x_fit)) ibad++;
    }
  }
  fCon->b_tac->GetEntry(entry);

  // Uncorrected y calibration with radomization.
  if(parent->fAddRndm && !fTACSplineCorr && !fAnodeSplineCorr){
    if(entry==0) printf("* Using uncorrected y calibration with radomization.\n");
    y    = y_offset[1] + y_slope[1]*y_offset[0] + y_slope[0]*y_slope[1]*(fCon->tac+fRandom.Rndm());
    yInt = y_offset[1] + y_slope[1]*y_offset[0];

  // TAC corrected y calibration without radomization.
  }else if(!parent->fAddRndm && fTACSplineCorr) {
    if(entry==0) printf("* Using TAC corrected y calibration without radomization.\n");
    // Apply TAC Correcting Cubic Spline.
    if(x>fTACSpline->GetXmin() && x<fTACSpline->GetXmax() && !isnan(x)){
      fCon->tac -= (UShort_t)(fTACSpline->Eval(x)-fTACSplineLevel);
    }
    tac = fCon->tac;
    y    = y_offset[1] + y_slope[1]*y_offset[0] + y_slope[0]*y_slope[1]*(Double_t)(fCon->tac);
    yInt = y_offset[1] + y_slope[1]*y_offset[0];
    if(fYSplineCorr){
      Double_t xkb,ykb,xke,yke;
      fYSpline->GetKnot(0,xkb,ykb);
      fYSpline->GetKnot(fYSpline->GetNp()-1,xke,yke);
      if     (strcmp(fYSplineMatchPnt,"begin")==0){
	if(y>xkb && y<=xke) y = fYSpline->Eval(y);
	if(y>xke)           y = fYSpline->Derivative(xke)*(y-xke)+yke; // point slope form for a line.
      }
      else if(strcmp(fYSplineMatchPnt,"end")==0){
	if(y<xke && y>=xkb) y = fYSpline->Eval(y);
	if(y<xkb)           y = fYSpline->Derivative(xkb)*(y-xkb)+ykb; // point slope form for a line.
      }
    }

  // TAC corrected y calibration with radomization.
  }else if(parent->fAddRndm && fTACSplineCorr) {
    if(entry==0) printf("* Using TAC corrected y calibration with radomization.\n");
    // Apply TAC Correcting Cubic Spline.
    if(x>fTACSpline->GetXmin() && x<fTACSpline->GetXmax() && !isnan(x)){
      fCon->tac -= (UShort_t)(fTACSpline->Eval(x)-fTACSplineLevel);
    }
    tac = fCon->tac;
    y    = y_offset[1] + y_slope[1]*y_offset[0] + y_slope[0]*y_slope[1]*(fCon->tac+fRandom.Rndm());
    yInt = y_offset[1] + y_slope[1]*y_offset[0];
    if(fYSplineCorr){
      Double_t xkb,ykb,xke,yke;
      fYSpline->GetKnot(0,xkb,ykb);
      fYSpline->GetKnot(fYSpline->GetNp()-1,xke,yke);
      if     (strcmp(fYSplineMatchPnt,"begin")==0){
	if(y>xkb && y<=xke) y = fYSpline->Eval(y);
	if(y>xke)           y = fYSpline->Derivative(xke)*(y-xke)+yke; // point slope form for a line.
      }
      else if(strcmp(fYSplineMatchPnt,"end")==0){
	if(y<xke && y>=xkb) y = fYSpline->Eval(y);
	if(y<xkb)           y = fYSpline->Derivative(xkb)*(y-xkb)+ykb; // point slope form for a line.
      }
    }

  // Anode corrected y calibration with radomization.
  }else if(parent->fAddRndm && fAnodeSplineCorr) {
    if(entry==0) printf("* Using Anode corrected y calibration with radomization.\n");
    fCon->b_anode->GetEntry(entry);
    if(calc.x_gravity>fAnodeSpline->GetXmin() && calc.x_gravity<fAnodeSpline->GetXmax() && !isnan(calc.x_gravity)){
      anode = (UShort_t)(fCon->anode - (UShort_t)(fAnodeSpline->Eval(calc.x_gravity)-fAnodeSplineLevel));
    }
    tac  = (UShort_t)(fCon->tac + (Double_t)(anode)*fAnodeCorrSlope + fAnodeCorrOffset);
    y    = y_offset[1] + y_slope[1]*y_offset[0] + y_slope[0]*y_slope[1]*(tac+fRandom.Rndm());
    yInt = y_offset[1] + y_slope[1]*y_offset[0];

  // Uncorrected y calibration.
  }else{
    if(entry==0) printf("* Using uncorrected y calibration without radomization.\n");
    y    = y_offset[1] + y_slope[1]*y_offset[0] + y_slope[0]*y_slope[1]*(Double_t)fCon->tac;
    yInt = y_offset[1] + y_slope[1]*y_offset[0];
  }

  // Apply any CRDC Shifts.
  // To apply y gain shift first wubtract off y-intersept, multiply by gain correction, then add back the y-intersept.
  //  y -= yInt;  y *= fYGainShift;  y += yInt;
  //  y *= fYGainShift;
  y += fYShift;  
  x += fXShift;
  
  //  cout << "* Failed calculations on " << name << ": " << 100.*ibad/i << " %\n";
  gNumrecMessage = restore;

  return 1;
}


//______________________________________________________________________________
void TS800Crdc::InitClass(TString iname,TS800 *theParent) 
{
  // -- Initialize the default class data.
  //

  parent = theParent;
  fName   = iname;
  
  fTACSplineCorr    = kFALSE;
  fTACSplineLevel   = 0;
  fAnodeSplineCorr  = kFALSE;
  fAnodeSplineLevel = 0;
  fYSplineCorr      = kFALSE;
  //  fYSplineLevel     = 0;
  debug    = false; 
  maxwidth = 32;
  filled   = 0;
  fAnodeCorrSlope  = 1;
  fAnodeCorrOffset = 0;
  x_offset[0] = -288;
  x_slope[0]  = 2.54;    // Constant value set by the finite pad size.
  y_offset[0] = 0;
  y_slope[0]  = 1;
  x_offset[1] = 0;
  x_slope[1]  = 1;       // Constant value set by the finite pad size.
  y_offset[1] = 0;
  y_slope[1]  = 1;
  fXShift     = 0;
  fYShift     = 0;
  fYGainShift = 0;
  
  pad.parent  = this;
  calc.parent = this;
  calc.pad    = &pad;
  channel     = new Short_t[4096];
  sample      = new Short_t[4096];
  data        = new Short_t[4096];
  pad.InitClass(fName+".pad");
  calc.InitClass(fName+".calc");


  // Prepare histos.
  Char_t   name[1000],title[1000];
  sprintf(name,"fHYShifts_%s",GetName());
  sprintf(title,"%s Y Shifts",GetName());
  fHYShifts.SetNameTitle(name,title);
  fHYShifts.SetBins(2000,1,2001);

  sprintf(name,"fHYGainShifts_%s",GetName());
  sprintf(title,"%s Y Gain Shifts",GetName());
  fHYGainShifts.SetNameTitle(name,title);
  fHYGainShifts.SetBins(2000,1,2001);
  for(Int_t i=0; i<fHYGainShifts.GetNbinsX(); i++) fHYGainShifts.SetBinContent(i,1.0);
  
}


//______________________________________________________________________________
void TS800Crdc::InitTree(TTree *tree) 
{
  // -- Initialize the class to tree.
  //

  fChain   = tree;
  fCurrent = -1;
  //  TNamed *aState = (TNamed*)fChain->GetUserInfo()->FindObject("Analysis State");

  b_data    = fChain->GetBranch(fName+".data");
  b_sample  = fChain->GetBranch(fName+".sample");
  b_channel = fChain->GetBranch(fName+".channel");
  b_filled  = fChain->GetBranch(fName+".filled");
  b_x       = fChain->GetBranch(fName+".x");
  b_y       = fChain->GetBranch(fName+".y");
  b_anode   = fChain->GetBranch(fName+".anode");
  b_tac     = fChain->GetBranch(fName+".tac");

  pad.InitTree(tree);
  calc.InitTree(tree);
}


//______________________________________________________________________________
void TS800Crdc::Clear(Option_t*) 
{
  // -- Clear the data members.
  //

  x     = sqrt(-1.0);
  y     = sqrt(-1.0);
  anode = 0;
  tac   = 0;
  //  pad.Clear();
  calc.Clear();
}


//______________________________________________________________________________
void TS800Crdc::Copy(TObject &crdc) const
{
  // -- Copy this method.

  TObject::Copy((TObject&)crdc);
  ((TS800Crdc&)crdc).channels       = channels;
  ((TS800Crdc&)crdc).tag            = tag;
  ((TS800Crdc&)crdc).maxwidth       = maxwidth;
  ((TS800Crdc&)crdc).threshold      = threshold;
  ((TS800Crdc&)crdc).sampleBegin    = sampleBegin;
  ((TS800Crdc&)crdc).sampleWidth    = sampleWidth;
  ((TS800Crdc&)crdc).debug          = debug;
  ((TS800Crdc&)crdc).fAnodeSplineCorr  = fAnodeSplineCorr;
  ((TS800Crdc&)crdc).fAnodeSpline      = fAnodeSpline;
  ((TS800Crdc&)crdc).fAnodeSplineLevel = fAnodeSplineLevel;
  ((TS800Crdc&)crdc).fTACSplineCorr    = fTACSplineCorr;
  ((TS800Crdc&)crdc).fTACSpline        = fTACSpline;
  ((TS800Crdc&)crdc).fTACSplineLevel   = fTACSplineLevel;
  ((TS800Crdc&)crdc).fYSplineCorr      = fYSplineCorr;
  ((TS800Crdc&)crdc).fYSpline          = fYSpline;
  ((TS800Crdc&)crdc).fYSplineMatchPnt  = fYSplineMatchPnt;
//   ((TS800Crdc&)crdc).fYSplineLevel     = fYSplineLevel;

  
  pad.Copy(((TS800Crdc&)crdc).pad);
  calc.Copy(((TS800Crdc&)crdc).calc);
  //  TS800CrdcCalc calc(calc);
  ((TS800Crdc&)crdc).pad.parent  = &((TS800Crdc&)crdc);
  ((TS800Crdc&)crdc).calc.parent = &((TS800Crdc&)crdc);
  //  ((TS800Crdc&)crdc).calc.pad    = ((TS800Crdc&)crdc).calc.pad;

  ((TS800Crdc&)crdc).fAnodeCorrOffset = fAnodeCorrOffset;
  ((TS800Crdc&)crdc).fAnodeCorrSlope  = fAnodeCorrSlope;
  ((TS800Crdc&)crdc).x_offset[0] = x_offset[0];
  ((TS800Crdc&)crdc).x_slope[0]  = x_slope[0];
  ((TS800Crdc&)crdc).y_offset[0] = y_offset[0];
  ((TS800Crdc&)crdc).y_slope[0]  = y_slope[0];
  ((TS800Crdc&)crdc).x_offset[1] = x_offset[1];
  ((TS800Crdc&)crdc).x_slope[1]  = x_slope[1];
  ((TS800Crdc&)crdc).y_offset[1] = y_offset[1];
  ((TS800Crdc&)crdc).y_slope[1]  = y_slope[1];
  ((TS800Crdc&)crdc).fXShift     = fXShift;
  ((TS800Crdc&)crdc).fYShift     = fYShift;
  ((TS800Crdc&)crdc).fYGainShift = fYGainShift;
  fHYShifts.Copy(((TS800Crdc&)crdc).fHYShifts);
  fHYGainShifts.Copy(((TS800Crdc&)crdc).fHYGainShifts);
  ((TS800Crdc&)crdc).fName       = fName;

  //  Clear();
}


//______________________________________________________________________________
void TS800Crdc::GetEntry(int i) 
{
  // --
  //

  //	pad.GetEntry(i);
  calc.GetEntry(i);
  //	b_data->GetEntry(i);
  //	b_sample->GetEntry(i);
  //	b_channel->GetEntry(i);
  //	b_filled->GetEntry(i);
  b_x->GetEntry(i);
  b_y->GetEntry(i);
  b_anode->GetEntry(i);
  b_tac->GetEntry(i);
}


//______________________________________________________________________________
Bool_t TS800Crdc::LoadYShifts(Char_t *yShiftFile)
{
  // -- Load the file containing the CRDC Y shifts.
  // 
  // Commented lines are denoted by the '#' character.  Each line has the following
  // format,
  //
  // <RUN NUMBER>     <MEAN>     <SHIFT>    <GAINSHIFT>
  //
  // NOTE:  Currently only the shift is used for the drift correction!!!
  //
  
  ifstream iFile(yShiftFile,ios::in);
  Char_t   line[1000];
  Int_t    runNum;
  Double_t shift;
  
  
  if(!iFile.is_open()) {printf("* File could not be opened!\n"); return kFALSE;}

  // Ignore commented lines.
  iFile.getline(line,1000);
  while(line[0]=='#' || line[0]==' ') iFile.getline(line,1000);

  // Read in file.
  Bool_t firstLine = kTRUE;
  while (!iFile.eof()){
    if(!firstLine) iFile.getline(line,1000);
    firstLine = kFALSE;
    //    sscanf(line,"%i %lf %le",&runNum,&shift,&gshift);
    sscanf(line,"%i %*f %lf %*e",&runNum,&shift);
    fHYShifts.SetBinContent(runNum,shift);
    //    fHYGainShifts.SetBinContent(runNum,gshift);
  }

  iFile.close();
  printf("* Loading %s Y shifts . . .  \t\t\t\t [OK]\n",GetName());
  
  return kTRUE;
}


//______________________________________________________________________________
Bool_t TS800Crdc::LoadSpline(Char_t *splineFile,Int_t crdcID, Int_t type)
{
  // -- Load the file containing the correcting splines.

  TDirectory *cDir = gDirectory;
  TFile spFile(splineFile,"READ");
  //  TSpline3 *sptmp;
  //  if(!file.IsOpen()){ printf("Error:  File %s could not be opened.\n",splineFile); return kFALSE;}

  Char_t spName[100];
  TSpline3 *tmpSp;
  if(type==0){
    sprintf(spName,"crdc%iAnodeSpline",crdcID);
    tmpSp = (TSpline3*)spFile.FindObjectAny(spName);
    fAnodeSpline = (TSpline3*)tmpSp->Clone(spName);
  }else if(type==1){
    sprintf(spName,"crdc%iTACSpline",crdcID);
    fTACSpline = (TSpline3*)spFile.FindObjectAny(spName);
    //    sptmp      = (TSpline3*)spFile.FindObjectAny(spName);
    //    fTACSpline = (TSpline3*)sptmp->Clone("sp_TACSpline");    
  }else if(type==2){
    sprintf(spName,"crdc%iYSpline",crdcID);
    fYSpline = (TSpline3*)spFile.FindObjectAny(spName);
    //    sptmp    = (TSpline3*)spFile.FindObjectAny(spName);
    //    fYSpline = (TSpline3*)sptmp->Clone("sp_YSpline");    
  }else{
    printf("* ERROR:  Invalid type.\n");
  }

  spFile.Close();
  cDir->cd();
  if(type==0 && !fAnodeSpline) {printf("Error:  %s function not found.\n",spName); return kFALSE;}
  if(type==1 && !fTACSpline  ) {printf("Error:  %s function not found.\n",spName); return kFALSE;}
  if(type==2 && !fYSpline)     {printf("Error:  %s function not found.\n",spName); return kFALSE;}
  else{printf("Loaded %s correcting function (Spline will be used during CRDC%i calculations).\n",spName,crdcID);}

  // Use correcting spline during calculation."
  if(type==0) fAnodeSplineCorr = kTRUE;
  if(type==1) fTACSplineCorr   = kTRUE;
  if(type==2) fYSplineCorr     = kTRUE;
  return kTRUE;
}


//______________________________________________________________________________
void TS800Crdc::SetYShift(Int_t runNum) 
{
  // -- 
  
  //  Int_t cbin;
  Double_t shift = fHYShifts.GetBinContent(runNum);
  fYShift = shift;
  //  if(cbin==runNum) fYShift = shift;
  //  else printf("* WARNING:  Bin %i and run number %i does not match!\n",cbin,runNum);
}


//______________________________________________________________________________
void TS800Crdc::SetYGainShift(Int_t runNum) 
{
  // -- 
  
  //  Int_t cbin;
  Double_t shift = fHYGainShifts.GetBinContent(runNum);
  fYGainShift = shift;
  //  if(cbin==runNum) fYShift = shift;
  //  else printf("* WARNING:  Bin %i and run number %i does not match!\n",cbin,runNum);
}


//______________________________________________________________________________
UShort_t* TS800Crdc::Unpack(UShort_t *p) 
{
  // -- Unpack the subpacket.
  //

  if(*(p+1) == tag+1){ // we are reading raw data
    p = Unpackrawdata(p);
  }
  if(*(p+1) == tag+5){ // we are reading anodes
    anode = (UShort_t)*(p+2);
    tac   = (UShort_t)*(p+3);
    p     += 4;
  }
  return p;
}


//______________________________________________________________________________
UShort_t* TS800Crdc::Unpackrawdata(UShort_t *p) 
{
  // -- Unpack the subpacket RAW data.
  //

  static ULong_t total=0, failed=0;
  Short_t isample, ichannel, cdata[4], connector, previous_sample=0, ch;
  
  memset(data, 0, sizeof(data));
  memset(channel, 0, sizeof(channel));
  memset(sample, 0, sizeof(sample));
  filled = 0;
  
  UShort_t *pStore = p;
  bool mes1=true, mes2=true, mes3=true, mes4=true;
  UShort_t length = *p++;
  Short_t  i      = length-3;
  p++;	// skip packet id
  threshold = *p++;
  while(i>0) {
    if((*p)>>15 != 1) {
      cout << "CRDC data is corrupted!" << endl;
      p++; i--;
      continue;
    }else{
      isample  = (Short_t)((*p)&0x7FC0)>>6;
      ichannel = (Short_t)(*p)&0x3F;
      if (i == length-3) {sampleBegin = isample; previous_sample = isample;}
    }
    p++; i--;
    memset(cdata, 0, sizeof(cdata));
    while((*p)>>15 == 0) {
      connector        = (Short_t)((*p)&0x0C00)>>10;
      cdata[connector] = (Short_t)(*p)&0x3FF;
      p++; i--;
      if(i==0) break;
    }
    if(isample < sampleBegin || isample > sampleBegin+maxwidth) {
      if(debug)
	printf("Warning in Crdc Unpack: inconsistent sample number: %d (first: %d)\n",
	       isample, sampleBegin);
      mes1 = false;
      continue;
    }
    if(isample < previous_sample) {
      if (debug)
	printf("Warning in Crdc Unpack: sample number lower than previous: %d (previous: %d)\n",
	       isample, previous_sample);
      mes2 = false;
      continue;
    }
    previous_sample = isample;
    for(Int_t j=0; j<4; j++) {
      ch = ichannel+j*64;
      if(cdata[j]!=0 && ch<channels) {
	if(cdata[j]<threshold) {
	  if(debug)
	    printf ("Warning in Crdc Unpack: data lower than threshold: %d (threshold: %d)\n",
		    cdata[j], threshold);
	  mes3 = false;
	}else{
	  data[filled]    = cdata[j];
	  channel[filled] = ch;
	  sample[filled]  = isample;
	  filled++;
	}
      }else if(cdata[j]!=0 && ch>=channels) {
	if(debug)
	  printf ("Warning in Crdc Unpack: channel greater than limit: %d (limit: %d)\n",
		  ch, channels);
	mes4 = false;
      }
    }
    sampleWidth = isample - sampleBegin + 1;
  }
  // integrate samples into pads
  //	if (filled > 0) pad.Integrate();
  if(!mes1 || !mes2 || !mes3 || !mes4) failed++;
  total++;
  if(failed==1000) {
    if (debug)
      printf ("Errors in Crdc Unpackings: %g%%\n", 1.0*failed/total*100);
    total = 0;
    failed = 0;
  }
  return (pStore+length);
}


//______________________________________________________________________________
void TS800Crdc::Snapshot1D(TH1F *histogram, Int_t entry) 
{
  // -- Draw the 1D snapshot.
  //

  histogram->Reset();
  pad.b_cal->GetEntry(entry);
  for (Int_t i=0; i<channels; i++) {
    if (!isnan(pad.cal[i])) histogram->Fill(i, pad.cal[i]);
  }
  histogram->Draw();
}


//______________________________________________________________________________
void TS800Crdc::Snapshot2D(TH2F *histogram, Int_t entry) 
{
  // -- Draw the 2D snapshot.
  //

  histogram->Reset();
  b_sample->GetEntry(entry);
  b_channel->GetEntry(entry);
  b_data->GetEntry(entry);
  b_filled->GetEntry(entry);
  for (int i=0; i<filled; i++) {
    histogram->Fill(channel[i], sample[i], data[i]);
  }
  histogram->Draw("COLZ");
}