//_____________________________________________________________________________
// CLC analysis module: below are the examples of what one can do in 
// the interactive session
// 
// note that because of the present limitations of the script compiler this
// module need to reside in the current directory
//_____________________________________________________________________________

/*

  ****** to run:

//  TStnAna x("dataval_ar0192ea_clc.root");
  gInterpreter->ProcessLine(".U /data12/murat/online/./TClcAna2Module.so");
  TStnAna x("dataval_clc.root");
  gInterpreter->ProcessLine(".L /data12/murat/online/./TClcAna2Module.cc+");
  x.AddModule(new TClcAna2Module());
					// using version 1 of the calib const
					// default is 0

  TClcAna2Module* m = (TClcAna2Module*) x.GetModule("ClcAna2");
  m->SetClcCalibVer(2);
  x.Run();
					// to fit residual miscalibration:

  TClcAna2Module* m = (TClcAna2Module*) x.GetModule("ClcAna2"); 
  m->FitW2VsAdc();
  m->Print("dw2_vs_adc");

  gStyle->SetStatW(0.13);
  gStyle->SetStatH(0.12);
  gStyle->SetStatX(0.4);
  gStyle->SetStatY(0.5);

  ****** to display a histogram

  TClcAna2Module::Hist_t* h;
  h = ((TClcAna2Module*) x.GetModule("ClcAna2"))->GetHist();
  h->fWidthVsAdcProf[1]->Draw();

  ****** to process one event (#15 in the ntuple) and to print the CLC data"

  x.ProcessEvent(15)
  TClcAna2Module* m = (TClcAna2Module*) x.GetModule("ClcAna2");
  TClcData* data = m->GetClcBlock();
  data->Print();

*/

#include "TROOT.h"
#include "TF1.h"
#include "TClcAna2Module.hh"

//_____________________________________________________________________________
TClcAna2Module::TClcAna2Module(const char* name, const char* title):
  TStnModule(name,title)
{

  fClcCalibVer  =  0;

  for (int i=0; i<96; i++) {
    fHist.fWidthVsTime    [i] = 0;
    fHist.fWidthVsTimeProf[i] = 0;
    //    fHist.fDtVsT2         [i] = 0;
    fHist.fW2VsTimeProf   [i] = 0;
    fHist.fT3VsTimeProf   [i] = 0;
    fHist.fWidthVsAdc     [i] = 0;
    fHist.fW2VsAdc        [i] = 0;
    fHist.fDw2VsAdc       [i] = 0;
    fHist.fWidthVsAdcProf [i] = 0;
    fHist.fDt3            [i] = 0;
    fHist.fDt4            [i] = 0;
  } 
}

//_____________________________________________________________________________
TClcAna2Module::~TClcAna2Module() {
  DeleteHistograms();
}

//_____________________________________________________________________________
void TClcAna2Module::DeleteHistograms() {
  // delete all the histograms

  if (fHist.fW2VsTimeProf[0]) {
				// delete all the histograms
    for (int i=0; i<96; i++) {
      delete fHist.fWidthVsTime    [i];
      //      delete fHist.fDtVsT2         [i];
      delete fHist.fWidthVsTimeProf[i];
      delete fHist.fW2VsTimeProf   [i];
      delete fHist.fT3VsTimeProf   [i];
      delete fHist.fWidthVsAdc     [i];
      delete fHist.fWidthVsAdcProf [i];
      delete fHist.fW2VsAdc        [i];
      delete fHist.fDw2VsAdc       [i];
      delete fHist.fDt3            [i];
      delete fHist.fDt4            [i];
    }

    fHist.fW2VsTimeProf[0] = 0;
  }
}

//_____________________________________________________________________________
void TClcAna2Module::BookHistograms() {
  char name [200];
  char title[200];
				// book histograms
  DeleteHistograms();

  for (int i=0; i<96; i++) {
     sprintf(name, "%s_wid_vs_time_%02i",GetName(),i);
     sprintf(title,"%s: TDC width vs time channel %02i",GetName(),i);
     fHist.fWidthVsTime[i] = new TH2F(name, title,500 ,900,1400,500,0,1000);

//       sprintf(name, "%s_dt_vs_t2_%02i",GetName(),i);
//       sprintf(title,"%s: TDC T2(i)-T2(ref) vs T2(ref) i=%02i",GetName(),i);
//       fHist.fDtVsT2[i] = new TH2F(name, title,200 ,0,100,200,-50,50);

    sprintf(name, "%s_wid_vs_time_prof%02i",GetName(),i);
    sprintf(title,"%s: TDC width vs time channel %02i (profile)",GetName(),i);
    fHist.fWidthVsTimeProf[i] = new TProfile(name, title,500 ,900,1400,0,1000);

    sprintf(name, "%s_w2_vs_time_prof%02i",GetName(),i);
    sprintf(title,"%s: TDC W2 vs time channel %02i (profile)",GetName(),i);
    fHist.fW2VsTimeProf[i] = new TProfile(name, title,500 ,900,1400,0,1000);

    sprintf(name, "%s_t3_vs_time_prof%02i",GetName(),i);
    sprintf(title,"%s: TDC T3(i) vs Time(i) channel %02i (profile)",GetName(),i);
    fHist.fT3VsTimeProf[i] = new TProfile(name, title,200,1000,1200,0,1000);

    sprintf(name, "%s_wid_vs_adc_prof_%02i",GetName(),i);
    sprintf(title,"%s: TDC width vs ADC channel %02i (profile)",GetName(),i);
    fHist.fWidthVsAdcProf[i] = new TProfile(name, title,40,300,20300,50,450);

    sprintf(name, "%s_wid_vs_adc_%02i",GetName(),i);
    sprintf(title,"%s: Width vs ADC channel %02i ",GetName(),i);
    fHist.fWidthVsAdc[i] = new TH2F(name, title,40,300,20300,200,50,450);

    sprintf(name, "%s_w2_vs_adc_%02i",GetName(),i);
    sprintf(title,"%s: W2 vs ADC channel %02i ",GetName(),i);
    fHist.fW2VsAdc[i] = new TH2F(name, title,40,300,20300,200,50,450);

    sprintf(name, "%s_dw2_vs_adc_%02i",GetName(),i);
    sprintf(title,"%s: Dw2 vs ADC channel %02i ",GetName(),i);
    fHist.fDw2VsAdc[i] = new TH2F(name, title,40,300,20300,300,-300,300);

    sprintf(name, "%s_dt3_%02i",GetName(),i);
    sprintf(title,"%s: Dt3 channel %02i ",GetName(),i);
    fHist.fDt3[i] = new TH1F(name, title,200,-20,20);

    sprintf(name, "%s_dt4_%02i",GetName(),i);
    sprintf(title,"%s: Dt4 channel %02i ",GetName(),i);
    fHist.fDt4[i] = new TH1F(name, title,100,-5,5);
  }

}


//_____________________________________________________________________________
int TClcAna2Module::BeginJob() {

				// register the data block

  fTriggerBlock = (TStnTriggerBlock*) RegisterDataBlock("TriggerBlock",
						       "TStnTriggerBlock");

  fClcBlock     = (TClcDataBlock*)   RegisterDataBlock("ClcDataBlock",
						       "TClcDataBlock");

					// book histograms
  BookHistograms();

  if (! fClcBlock) {
    printf(" >>> branch ** %s ** doesn't exist, disable module \n","ClcBlock");
    fEnabled = 0;
  }

  return 0;
}

//_____________________________________________________________________________
int TClcAna2Module::BeginRun() {
  // initialize run-dependent calibration constants
  
  fClcBlock->SetCalibVersion(fClcCalibVer);
  int rc =  fClcBlock->ReadCalibrations(GetHeaderBlock()->RunNumber());
  return rc;
}

//_____________________________________________________________________________
int TClcAna2Module::Event(int ientry) {
					// in case of a TChain, ientry is the 
					// entry number in the current file
					// read the trigger block to look at
					// the myron flag
  fTriggerBlock->GetEntry(ientry);

  TTsid* tsid = fTriggerBlock->Tsid();
  TTfrd* tfrd = fTriggerBlock->Tfrd();
					// process only requested bucket
  int myron = tsid->MyronFlag();
  if (fMyronFlag >= 0) {
    if (myron != fMyronFlag) return 0;
  }
					// myron is OK, read CLC data
  fClcBlock->GetEntry(ientry);
					// do whatever you want
  if (PrintLevel()) {
    fClcBlock->Print();
  }
					// define reference channels for the 
					// delay calibration, use channels 
					// #1 (west) and #49(east)
  TClcChannel  *ch, *ref_ch;
  float        t0, a0, width, w2, dw2, dt3, dt4;
  int          n_tdc_hits, ich0;

  for (int i=0; i<96; i++) {
    if (i < 48) ich0 = fClcBlock->Module(0)->RefChannel();
    else        ich0 = fClcBlock->Module(1)->RefChannel();

    ref_ch = fClcBlock->Channel(ich0);

    ch         = fClcBlock->Channel(i);
    n_tdc_hits = ch->NTdcHits();

    if (n_tdc_hits == 1) {

      a0     = ch->AdcCounts();
      t0     = ch->TdcWord(0)->LeadingEdge()+132*myron;
      width  = ch->Width(0);
      w2     = ch->W2(0);

      if (ch->AdcCounts() > 300) {
//  	fHist.fWidthVsTime[i]->Fill    (t0,width,1);
	fHist.fWidthVsAdcProf[i]->Fill (a0,width,1);
	fHist.fWidthVsAdc    [i]->Fill (a0,width);
	fHist.fW2VsAdc       [i]->Fill (a0,w2);
	if ( (t0 > 1060) && (t0 < 1100)) {
	  fHist.fWidthVsTimeProf[i]->Fill(t0,width,1);
	  fHist.fW2VsTimeProf[i]->Fill(t0,w2,1);
	  fHist.fT3VsTimeProf[i]->Fill(t0,ch->T3(0),1);
	}
      }
      
      if ((ref_ch->NTdcHits() == 1) && (ref_ch->AdcCounts() > 1500)) {
	dw2    = ref_ch->W2(0)-ch->W2(0);
	fHist.fDw2VsAdc[i]->Fill(a0,dw2);
	if (ch->AdcCounts() > 1500) {
	  dt3    = ref_ch->T3(0)-ch->T3(0);
	  fHist.fDt3[i]->Fill(dt3);
	  dt4    = ref_ch->BestTime(0)-ch->BestTime(0);
	  fHist.fDt4[i]->Fill(dt4);
	}
      }
    }
  }

  return 0;			       
}

//_____________________________________________________________________________
Double_t TClcAna2Module::FitFunction(Double_t *x, Double_t *par) {
  // parametrize threshold dependence of signal time with exponential function

  Double_t arg, f;
  if (par[2]) arg = (x[0] - 300.)/par[2];
  else        arg = 0;

  f = par[0]+par[1]*TMath::Exp(-arg)-par[3]*x[0];
  return f;
}

//_____________________________________________________________________________
Double_t TClcAna2Module::Pol2(Double_t *x, Double_t *par) {
  // 2nd order polinomial

  Double_t f = par[0]+par[1]*x[0]+par[2]*x[0]*x[0];
  return f;
}

//_____________________________________________________________________________
int TClcAna2Module::FitWidthVsAdc() {

  if (fHist.fFitT0) {
    delete fHist.fFitT0;
    delete fHist.fFitT1;
    delete fHist.fFitA0;
    delete fHist.fFitB0;
    delete fHist.fFitChi2;
    delete fHist.fFitDeltaT;
  }

  fHist.fFitT0     = new TH1F("ClcAna2_fit_t0",
			      "ClcAna2: fit T0",100,0,100);
  fHist.fFitT1     = new TH1F("ClcAna2_fit_t1",
			      "ClcAna2: fit T1 ",100,0,100);
  fHist.fFitA0     = new TH1F("ClcAna2_fit_a0",
			      "ClcAna2: fit A0 ",100,0,100);
  fHist.fFitB0     = new TH1F("ClcAna2_fit_b0",
			      "ClcAna2: fit B0 ",100,0,100);
  fHist.fFitChi2   = new TH1F("ClcAna2_fit_chi2",
			      "ClcAna2: fit chi2",100,0,100);

  fHist.fFitDeltaT = new TH1F("ClcAna2_fit_dt",
			      "ClcAna2: relative delays",100,0,100);

  TF1 *func = new TF1("ClcAna2_FitFunction",TClcAna2Module::FitFunction,
		      300,5000,4);

  func->SetParNames("T0","T1","A0","B0");
  func->SetParameters(200,20,200,-0.001);
  func->SetLineWidth(1);

  char name [100];
  for (int i=0; i<96; i++) {
    sprintf(name, "%s_wid_vs_adc_prof_%02i",GetName(),i);
    func->SetParameters(200,20,200,-0.001);
    TProfile* hist = (TProfile*) gROOT->FindObject(name);
    hist->SetMarkerStyle(20);
    hist->SetMarkerSize(0.5);
    hist->Fit("ClcAna2_FitFunction","qr0","");
  }

  TF1* f;
  for (int i=0; i<96; i++) {
    sprintf(name, "%s_wid_vs_adc_prof_%02i",GetName(),i);
    TProfile* hist = (TProfile*) gROOT->FindObject(name);

    f = hist->GetFunction("ClcAna2_FitFunction");
    
    fHist.fFitT0->SetBinContent(i+1,f->GetParameter(0));
    fHist.fFitT0->SetBinError(i+1,f->GetParError(0));

    fHist.fFitT1->SetBinContent(i+1,f->GetParameter(1));
    fHist.fFitT1->SetBinError(i+1,f->GetParError(1));

    fHist.fFitA0->SetBinContent(i+1,f->GetParameter(2));
    fHist.fFitA0->SetBinError(i+1,f->GetParError(2));

    fHist.fFitB0->SetBinContent(i+1,f->GetParameter(3));
    fHist.fFitB0->SetBinError(i+1,f->GetParError(3));

    fHist.fFitChi2->SetBinContent(i+1,f->GetChisquare());

  }
  return 0;
}



//_____________________________________________________________________________
int TClcAna2Module::FitW2VsAdc() {
  // fit residual dependence of W2 on Adc

  TF1 *func = new TF1("ClcAna2_pol2",TClcAna2Module::Pol2,
		      300,6000,3);

  func->SetParNames("A0","A1","A2");
  func->SetLineWidth(1);
  
  TProfile* hprof;

  char name [100];
  for (int i=0; i<96; i++) {
    sprintf(name, "%s_w2_vs_adc_%02i",GetName(),i);
    func->SetParameters(200,0,0);
    TH2F* hist = (TH2F*) gROOT->FindObject(name);
    hprof = hist->ProfileX();
    hprof->SetMarkerStyle(20);
    hprof->SetMarkerSize(0.5);
    hprof->Fit("ClcAna2_pol2","qr0","");
  }
  
  return 0;
}




//_____________________________________________________________________________
int TClcAna2Module::FitGauss(const char* opt) {
  // fit T0 offsets when widths are used

  char name [100], format[100];

  for (int i=0; i<96; i++) {
    if (strcmp(opt,"dw2") == 0) {
      strcpy(format,"%s_dw2_vs_adc_%02i");
      sprintf(name, format,GetName(),i);
      TH2F* hist = (TH2F*) gROOT->FindObject(name);
      TH1D* hist_py = hist->ProjectionY();
      hist_py->SetMarkerStyle(20);
      hist_py->SetMarkerSize(0.5);
      hist_py->Fit("gaus","q0");
      hist_py->GetFunction("gaus")->SetLineWidth(1);
    }
    else if (strcmp(opt,"dt3") == 0) {
      strcpy(format,"%s_dt3_%02i");
      sprintf(name, format,GetName(),i);
      TH1F* hist = (TH1F*) gROOT->FindObject(name);
      hist->SetMarkerStyle(20);
      hist->SetMarkerSize(0.5);
      hist->Fit("gaus","q0");
      hist->GetFunction("gaus")->SetLineWidth(1);
    }
    else if (strcmp(opt,"dt4") == 0) {
      strcpy(format,"%s_dt4_%02i");
      sprintf(name, format,GetName(),i);
      TH1F* hist = (TH1F*) gROOT->FindObject(name);
      hist->SetMarkerStyle(20);
      hist->SetMarkerSize(0.5);
      hist->Fit("gaus","q0");
      hist->GetFunction("gaus")->SetLineWidth(1);
    }
  }
  
  return 0;
}


//_____________________________________________________________________________
int TClcAna2Module::EndRun() {
  return 0;
}

//_____________________________________________________________________________
int TClcAna2Module::EndJob() {
  printf("----- end job: ---- %s\n",GetName());
  return 0;
}

//_____________________________________________________________________________
void TClcAna2Module::Print(Option_t* opt) const {

  char name[200];

  if (opt == 0) return ;

  if (strcmp(opt,"wid_vs_adc") == 0) {
					// print fit results

    printf("----------------------------------------");
    printf("----------------------------------------\n");
    printf("ch#     T0     sig(T0)   T1    sig(T1) ");
    printf("A0    sig(A0)     B1    sig(B1) chi**2\n");
    printf("----------------------------------------");
    printf("----------------------------------------\n");

    for (int i=0; i<96; i++) {
					// print the fit parameters

      TF1* f = fHist.fWidthVsAdcProf[i]->GetFunction("ClcAna2_FitFunction");

      if (f) {
	printf (" %2i  %8.3f %7.3f %7.3f %7.3f  %7.3f %7.3f %9.5f %7.3f %7.3f \n",
		i, 
		(float) f->GetParameter(0), (float) f->GetParError(0),
		(float) f->GetParameter(1), (float) f->GetParError(1),
		(float) f->GetParameter(2), (float) f->GetParError(2),
		(float) f->GetParameter(3), (float) f->GetParError(3),
		(float) f->GetChisquare());
      }
      else {
	printf(" %2i fit error\n",i);
      }
    }
  }
  if (strcmp(opt,"dw2") == 0) {
					// print results of dw2 fit

    printf("---------------------------------------------------------\n");
    printf("ch#  bad    mean  sig(mean)  dw2     sig(dw2)  chi**2    \n");
    printf("---------------------------------------------------------\n");

    for (int i=0; i<96; i++) {
					// print the fit parameters

      sprintf(name,"ClcAna2_dw2_vs_adc_%02i_py",i);
      TH1D* hist = (TH1D*) gROOT->FindObject(name);
      TF1* f = hist->GetFunction("gaus");

      if (f) {
	printf (" %2i  %2i %8.3f %8.3f %8.3f %8.3f %8.3f \n",
		i, 0,
		(float) f->GetParameter(1), (float) f->GetParError(1),
		(float) f->GetParameter(2), (float) f->GetParError(2),
		(float) f->GetChisquare());
      }
      else {
	printf(" %2i fit error\n",i);
      }
    }
  }
  if (strcmp(opt,"dt3") == 0) {
					// print results of dt3 fit

    printf("---------------------------------------------------------\n");
    printf("ch#  bad    T0  sig(T0) width sig(Width)    chi**2       \n");
    printf("---------------------------------------------------------\n");

    for (int i=0; i<96; i++) {
					// print the fit parameters

      sprintf(name,"ClcAna2_dt3_%02i",i);
      TH1D* hist = (TH1D*) gROOT->FindObject(name);
      TF1* f = hist->GetFunction("gaus");

      if (f) {
	printf (" %2i  %2i %8.3f %8.3f %8.3f %8.3f %8.3f \n",
		i, 0,
		(float) f->GetParameter(1), (float) f->GetParError(1),
		(float) f->GetParameter(2), (float) f->GetParError(2),
		(float) f->GetChisquare());
      }
      else {
	printf(" %2i fit error\n",i);
      }
    }
  }
  if (strcmp(opt,"dt4") == 0) {
					// print results of dt3 fit

    printf("---------------------------------------------------------\n");
    printf("ch#  bad    T0  sig(T0) width sig(Width)    chi**2       \n");
    printf("---------------------------------------------------------\n");

    for (int i=0; i<96; i++) {
					// print the fit parameters

      sprintf(name,"ClcAna2_dt4_%02i",i);
      TH1D* hist = (TH1D*) gROOT->FindObject(name);
      TF1* f = hist->GetFunction("gaus");

      if (f) {
	printf (" %2i  %2i %8.3f %8.3f %8.3f %8.3f %8.3f \n",
		i, 0,
		(float) f->GetParameter(1), (float) f->GetParError(1),
		(float) f->GetParameter(2), (float) f->GetParError(2),
		(float) f->GetChisquare());
      }
      else {
	printf(" %2i fit error\n",i);
      }
    }
  }
  else if (strcmp(opt, "calib") == 0) {
					// print calibration constants
    
  printf("---------------------------------");
  printf("----------------------------------------\n");
  printf(" ch  bad      T0     T1W      A0W      B1W        B2W   ");
  printf("    Slope   T0W \n");
  printf("---------------------------------");
  printf("----------------------------------------\n");

    for (int i=0; i<96; i++) {

      TClcChannel* ch = fClcBlock->Channel(i);

      float t0    = ch->T0();
      float t1w   = ch->T1W();
      float a0w   = ch->A0W();
      float b1w   = ch->B1W();
      float b2w   = ch->B2W();
      float t0w   = ch->T0W();
      float slope = ch->Slope();

      sprintf(name,"ClcAna2_w2_vs_adc_%02i_pfx",i);
      TProfile* hist = (TProfile*) gROOT->FindObject(name);
      TF1* f = hist->GetFunction("ClcAna2_pol2");
      if (f) {
				// T0W offsets still need to be redefined

	b1w     -= f->GetParameter(1);
	b2w     -= f->GetParameter(2);
      }
      else {
	printf(" %2i fit error\n",i);
      }
      printf (" %2i  %2i  %8.3f %7.3f  %7.3f %10.3e  %10.3e %7.3f %7.3f \n",
	      i, ch->Bad(), t0, t1w,    a0w,  b1w,   b2w, slope, t0w);
    }
  }
}