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

#include "TS800Scintillator.h"
#include "TS800.h"

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

#include <TBranch.h>
#include <TDirectory.h>
#include <TFile.h>
#include <TH1D.h>
#include <TF1.h>
#include <TMath.h>
#include <TSpline.h>

ClassImp(TS800Scintillator);

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// TS800Scintillator                                                    //
//                                                                      //
//                                                                      //
//////////////////////////////////////////////////////////////////////////



//______________________________________________________________________________
TS800Scintillator::TS800Scintillator(const TS800Scintillator & eN) : TObject(eN)
{
  // -- Copy Constructor.

  ((TS800Scintillator&)eN).Copy(*this);
}


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

  de_up     = 0;
  de_down   = 0;
  time_up   = 0;
  time_down = 0;
  tup       = sqrt(-1.0);
  tdown     = sqrt(-1.0);
  de        = sqrt(-1.0);
  time      = sqrt(-1.0);
}


//______________________________________________________________________________
void TS800Scintillator::Calculate(Long64_t entry) 
{
  // -- Calculate an entry.
  //

  Clear();

  // Get the entry.
  fCon->b_time_up->GetEntry(entry);
  fCon->b_time_down->GetEntry(entry);
  fCon->b_de_up->GetEntry(entry);
  fCon->b_de_down->GetEntry(entry);
  fCon->parent->trigger.b_s800->GetEntry(entry);
  // Generate the trigger randomize number.  This maintains the same random number to 
  // convert from UShort_t to Double_t for this entry.  This same number is used in other
  // calculations with the trigger.
  fCon->parent->trigger.fRs800 = fRandom.Rndm();
  Double_t trigRndm = fCon->parent->trigger.fRs800;

  // Calculate quantities.
  if(parent->fAddRndm){
    tup   = (Double_t)(fCon->time_up+fRandom.Rndm()   - fCon->parent->trigger.s800+trigRndm)*timeSlope;
    tdown = (Double_t)(fCon->time_down+fRandom.Rndm() - fCon->parent->trigger.s800+trigRndm)*timeSlope;
    de    = (Double_t)sqrt((Double_t)((fCon->de_up+fRandom.Rndm())*fCon->de_up) + 
			   (Double_t)((fCon->de_down+fRandom.Rndm())*fCon->de_down));
  }else{
    tup   = (Double_t)(fCon->time_up   - fCon->parent->trigger.s800)*timeSlope;
    tdown = (Double_t)(fCon->time_down - fCon->parent->trigger.s800)*timeSlope;
    de    = (Double_t)sqrt((Double_t)(fCon->de_up*fCon->de_up) + (Double_t)(fCon->de_down*fCon->de_down));
  }
  
  #ifdef EXPERIMENT_02023
  time  = (Double_t)(-1*tup);              // Problems with e1.time_down.  EXCLUDING
  #else
  time  = (Double_t)-1*(tup + tdown)/2;
  #endif
  
  if(fE1SplineCorr){
    //     printf("\tCorrecting e1.time using loaded cubic spline.\n");
    //    fCon->parent->crdc2.calc.b_x_gravity->GetEntry(i);
    Correct(e1Spline);
  }
}


//______________________________________________________________________________
void TS800Scintillator::Calibrate(Long64_t entry)
{
  // -- Calibrate an entry.
  //
 
  // Calibration is done in the Calculate method.
  printf("* This method is currently empty.\n");
}


//______________________________________________________________________________
void TS800Scintillator::Copy(TObject &eN) const
{
  // -- Copy this method.
  TObject::Copy((TObject&)eN);
  ((TS800Scintillator&)eN).fRandom.SetSeed(0);
  ((TS800Scintillator&)eN).e1Spline      = e1Spline;
  ((TS800Scintillator&)eN).fE1SplineCorr = fE1SplineCorr;

  ((TS800Scintillator&)eN).name          = name;
  ((TS800Scintillator&)eN).avgTimeOffset = avgTimeOffset;
  ((TS800Scintillator&)eN).timeUpSlope   = timeUpSlope;
  ((TS800Scintillator&)eN).timeDownSlope = timeDownSlope;
  ((TS800Scintillator&)eN).timeSlope     = timeSlope;
  //  Clear();
}

//______________________________________________________________________________
void TS800Scintillator::Correct(TSpline3 *e1Spline)
{
  // -- Correct time for the current entry using a defined cubic spline.
  //

  // Class members now reset to NaNs		
  if(parent->crdc2.calc.x_gravity>fCon->parent->e1.e1Spline->GetXmin() &&
     parent->crdc2.calc.x_gravity<fCon->parent->e1.e1Spline->GetXmax()){
    time = time + (Double_t)(avgTimeOffset - (-1)*fCon->parent->e1.e1Spline->Eval(parent->crdc2.calc.x_gravity));
  }
}


//______________________________________________________________________________
void TS800Scintillator::Correct(Double_t Mean,Double_t Sigma,void *func)
{
  // -- Correct "shadow" peak in time spectrum using a gaussian distribution.
  //    OBSOLETE

  TH1D *h_e1Distro = new TH1D("h_e1Distro","h_e1_Distro",3000,-275,-245);
  TF1  *fnc = new TF1("fnc",func,-30,30,2);
  fnc->SetParameter(0,Mean);
  fnc->SetParameter(1,Sigma);
  h_e1Distro->FillRandom("fnc",100000);

  // First reset branches of parameters to calculate
  b_time->Reset();
  
  // Then loop on entries reading only the required data, calculate, and fill
  for (Long64_t i=0; i<fChain->GetEntriesFast(); i++) {
    Clear();
    // Class members now reset to NaNs		
    
    b_time_up->GetEntry(i);
    parent->trigger.b_s800->GetEntry(i);
    // Class members now loaded with entry i
    tup   = time_up - parent->trigger.s800;
    time  = (Double_t)tup ;              // Problems with e1.time_down.  EXCLUDING

    time  = h_e1Distro->GetRandom();
    
    b_time->Fill();
  }
}


//______________________________________________________________________________
void TS800Scintillator::InitClass(TString iname, TS800* iparent) 
{
  // -- Initialze default data members.
  //

  name   = iname;
  parent = iparent;

  fE1SplineCorr = kFALSE;
  //  fE2SplineCorr = kFALSE;

  avgTimeOffset = 0;
  timeDownSlope = 0;
  timeSlope     = 1.000;   // 100 ps/ch

}


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

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

  b_de_up     = fChain->GetBranch(name+".de_up");
  b_de_down   = fChain->GetBranch(name+".de_down");
  b_time_up   = fChain->GetBranch(name+".time_up");
  b_time_down = fChain->GetBranch(name+".time_down");
  b_tup       = fChain->GetBranch(name+".tup");
  b_tdown     = fChain->GetBranch(name+".tdown");
  b_de        = fChain->GetBranch(name+".de");
  b_time      = fChain->GetBranch(name+".time");
}


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

  //  b_de_up->GetEntry(i);
  //   b_de_down->GetEntry(i);
  //   b_time_up->GetEntry(i);
  //   b_time_down->GetEntry(i);
  //   b_tup->GetEntry(i);
  //   b_tdown->GetEntry(i);
  //   b_de->GetEntry(i);
  //   b_time->GetEntry(i);
}


//______________________________________________________________________________
Bool_t TS800Scintillator::LoadSpline(Char_t *splineFile)
{
  // -- Load the file containing the correcting splines.
  // The E1 scintillator (possibly E2 and E3 also) in the S800 display 
  // reflection effects.  This can be seen in the dependance of time_up vs 
  // crdc2.calc.x_gravity.  To correct this dependance a ProfileX histo can 
  // be made and fit with a cubic spline.  This spline function is should then 
  // be renamed to e1Spline, e2Spline, or e3Spline and saved to a root file. 
  // This method loads the that file and sets the spline functions.
  

  TDirectory *cDir = gDirectory;
  TFile file(splineFile,"READ");

  e1Spline = (TSpline3*)file.FindObjectAny("e1Spline");

  file.Close();
  cDir->cd();
  if(!e1Spline) {printf("Error:  e1Spline function not found.\n"); return kFALSE;}
  else{printf("Loaded e1Spline correcting function (Spline will be used during scintillator calculations).\n");}

  // Use correcting spline during calculation."
  fE1SplineCorr = kTRUE;
  return kTRUE;
}


//______________________________________________________________________________
UShort_t* TS800Scintillator::Unpack(UShort_t *p, UShort_t ID) 
{
  // -- Unpack a subpacket.
  //

  // Even ID: UP channels.  Odd ID: DOWN channels
  if ((ID/2)*2 == ID) {
    de_up     = (UShort_t)((*p++)&0xfff);
    time_up   = (UShort_t)((*p++)&0xfff);
  } else {
    de_down   = (UShort_t)((*p++)&0xfff);
    time_down = (UShort_t)((*p++)&0xfff);
  }
   return p;
}


// Double_t TS800Scintillator::e1GDisto(Double_t *x, Double_t *par)
// {
//   Double_t val = 0.;
//   val = 1/(par[1]*Sqrt(2*Pi()))*Exp(-pow((x[0]-par[0]),2)/(2*par[1]));
//   return val;
// }