// ---------------------------------------------------------------------------------
// -----                           PndSttRawAna                                -----
// -----                        Author: S.Costanza                             -----
// ---------------------------------------------------------------------------------

#include <iostream>
#include <fstream>
using namespace std;

//ROOT headers
#include "TClonesArray.h"
#include "TH1F.h"
#include "TF1.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TMatrixT.h"
#include "TGraph.h"
#include "TMultiGraph.h"

//Fair heaeders
#include "FairRootManager.h"
#include "TFolder.h"
#include "TList.h"

//Stt headers
#include "PndSttRawHit.h"
#include "PndSttRawAna.h"

PndSttRawAna::PndSttRawAna()
{ 
  fPersistence = kTRUE;
  fevtn = 0;
  fVerbose = 0;
}


PndSttRawAna::~PndSttRawAna()
{
}


InitStatus PndSttRawAna::Init()
{
  FairRootManager *ioman = FairRootManager::Instance();
  if ( ! ioman ) {
    cout << "-E- PndSttRawAna::Init: "
	 << "RootManager not instantiated!" << endl;
    return kFATAL;
  }

  fRawData = (TClonesArray*) ioman->GetObject("SttRawHit");
  if (!fRawData) {
    cout << "-E- PndSttRawAnat::Init: "
	 << "No input raw array!" << endl;
    return kERROR;
  }

 
  fInFile = ioman->GetOutFile();
  fInFile->cd();
  fInFile->mkdir("Histos");

  CreateHistos();

 return kSUCCESS;

}

void PndSttRawAna::Exec(Option_t *option)
{
  if(fVerbose==0 && fevtn%1000==0) 
    cout << "Event Number "<<fevtn<<endl;
  else if(fVerbose >= 3) cout << "Event Number "<<fevtn<<endl;

  fevtn++;
  
  Int_t layer = 0;
  Int_t tubeID = 0;
  Double_t rawled = 0., rawzct = 0., rawcfd = 0.;
  Double_t fpgaled = 0., fpgazct = 0., fpgacfd = 0.;

  PndSttRawHit *hit;
  for (Int_t i = 0; i<fRawData->GetEntriesFast(); i++) {
    hit = (PndSttRawHit*) fRawData->At(i);
    
    layer = hit->GetLayer();
    tubeID = hit->GetTubeID();

    rawled = hit->GetRawLed();
    rawzct = hit->GetRawZct();
    rawcfd = hit->GetRawCfd();
    hrawled0[tubeID]->Fill(rawled);
    hrawzct0[tubeID]->Fill(rawzct);
    hrawcfd0[tubeID]->Fill(rawcfd);

    fpgaled = hit->GetFpgaLed();
    fpgazct = hit->GetFpgaZct();
    fpgacfd = hit->GetFpgaCfd();
    hfpgaled0[tubeID]->Fill(fpgaled);
    hfpgazct0[tubeID]->Fill(fpgazct);
    hfpgacfd0[tubeID]->Fill(fpgacfd);
  }
}


void PndSttRawAna::FinishTask()
{
  ReverseTDCspectra();
  WriteHistos();
  FindDeltaT();
  WriteSumHisto();
  WriteOffset();
  Calibration();
  ClearHistos();
}


void PndSttRawAna::CreateHistos()
{
  for (int i=1; i<fNTubes; i++) {

    // raw drift time plots
    TString hrawled0name = "rawLEDrev_"; hrawled0name += i;
    TString hrawzct0name = "rawZCTrev_"; hrawzct0name += i;
    TString hrawcfd0name = "rawCFDrev_"; hrawcfd0name += i;
    TString hrawledname  = "rawLED_";    hrawledname  += i;
    TString hrawzctname  = "rawZCT_";    hrawzctname  += i;
    TString hrawcfdname  = "rawCFD_";    hrawcfdname  += i;

//     hrawled0[i] = new TH1F(hrawled0name, hrawled0name, 100, -1000, -400);
//     hrawzct0[i] = new TH1F(hrawzct0name, hrawzct0name, 100, -1000, -400);
//     hrawcfd0[i] = new TH1F(hrawcfd0name, hrawcfd0name, 100, -1000, -400);
//     hrawled[i]  = new TH1F(hrawledname,  hrawledname,  100, 400, 1000);
//     hrawzct[i]  = new TH1F(hrawzctname,  hrawzctname,  100, 400, 1000);
//     hrawcfd[i]  = new TH1F(hrawcfdname,  hrawcfdname,  100, 400, 1000);
    hrawled0[i] = new TH1F(hrawled0name, hrawled0name, 100, 0, 700);
    hrawzct0[i] = new TH1F(hrawzct0name, hrawzct0name, 100, 0, 700);
    hrawcfd0[i] = new TH1F(hrawcfd0name, hrawcfd0name, 100, 0, 700);
    hrawled[i]  = new TH1F(hrawledname,  hrawledname,  100, 0, 700);
    hrawzct[i]  = new TH1F(hrawzctname,  hrawzctname,  100, 0, 700);
    hrawcfd[i]  = new TH1F(hrawcfdname,  hrawcfdname,  100, 0, 700);

    // FPGA drift time plots
    TString hfpgaled0name = "fpgaLEDrev_"; hfpgaled0name += i;
    TString hfpgazct0name = "fpgaZCTrev_"; hfpgazct0name += i;
    TString hfpgacfd0name = "fpgaCFDrev_"; hfpgacfd0name += i;
    TString hfpgaledname  = "fpgaLED_";    hfpgaledname  += i;
    TString hfpgazctname  = "fpgaZCT_";    hfpgazctname  += i;
    TString hfpgacfdname  = "fpgaCFD_";    hfpgacfdname  += i;

    hfpgaled0[i] = new TH1F(hfpgaled0name, hfpgaled0name, 100, 0, 700);
    hfpgazct0[i] = new TH1F(hfpgazct0name, hfpgazct0name, 100, 10, 700);
    hfpgacfd0[i] = new TH1F(hfpgacfd0name, hfpgacfd0name, 100, 0, 700);
    hfpgaled[i]  = new TH1F(hfpgaledname,  hfpgaledname,  100, 0, 700);
    hfpgazct[i]  = new TH1F(hfpgazctname,  hfpgazctname,  100, 10, 700);
    hfpgacfd[i]  = new TH1F(hfpgacfdname,  hfpgacfdname,  100, 0, 700);
  }

  hsumled = new TH1F("sumLED", "sumLED", 100, 0, 600);
}


void PndSttRawAna::ReverseTDCspectra() 
{
  for (int i=1; i<fNTubes; i++) {
    for (int bin=1; bin<101; bin++) {
//       hrawled[i]->SetBinContent(bin, hrawled0[i]->GetBinContent(101-bin));
//       hrawzct[i]->SetBinContent(bin, hrawzct0[i]->GetBinContent(101-bin));
//       hrawcfd[i]->SetBinContent(bin, hrawcfd0[i]->GetBinContent(101-bin));
      hrawled[i]->SetBinContent(bin, hrawled0[i]->GetBinContent(bin));
      hrawzct[i]->SetBinContent(bin, hrawzct0[i]->GetBinContent(bin));
      hrawcfd[i]->SetBinContent(bin, hrawcfd0[i]->GetBinContent(bin));

      hfpgaled[i]->SetBinContent(bin, hfpgaled0[i]->GetBinContent(101-bin));
      hfpgazct[i]->SetBinContent(bin, hfpgazct0[i]->GetBinContent(101-bin));
      hfpgacfd[i]->SetBinContent(bin, hfpgacfd0[i]->GetBinContent(101-bin));
    }
  }
}


void PndSttRawAna::WriteHistos()
{
  fInFile->cd("Histos");
  for (int i=0; i<8; i++) {
    TString cnrled = "RawLED_layer_"; cnrled += i+1;
    TString cnrzct = "RawZCT_layer_"; cnrzct += i+1;
    TString cnrcfd = "RawCFD_layer_"; cnrcfd += i+1;
    crled[i] = new TCanvas(cnrled, cnrled, 900, 900);
    crzct[i] = new TCanvas(cnrzct, cnrzct, 900, 900);
    crcfd[i] = new TCanvas(cnrcfd, cnrcfd, 900, 900);
    crled[i]->Divide(4,4);
    crzct[i]->Divide(4,4);
    crcfd[i]->Divide(4,4);

    TString cnfled = "FpgaLED_layer_"; cnfled += i+1;
    TString cnfzct = "FpgaZCT_layer_"; cnfzct += i+1;
    TString cnfcfd = "FpgaCFD_layer_"; cnfcfd += i+1;
    cfled[i] = new TCanvas(cnfled, cnfled, 900, 900);
    cfzct[i] = new TCanvas(cnfzct, cnfzct, 900, 900);
    cfcfd[i] = new TCanvas(cnfcfd, cnfcfd, 900, 900);
    cfled[i]->Divide(4,4);
    cfzct[i]->Divide(4,4);
    cfcfd[i]->Divide(4,4);

    for (int j=1; j<=16; j++) {
      crled[i]->cd(j); hrawled[i*16+j]->Draw();
      crzct[i]->cd(j); hrawzct[i*16+j]->Draw();
      crcfd[i]->cd(j); hrawcfd[i*16+j]->Draw();
      cfled[i]->cd(j); hfpgaled[i*16+j]->Draw();
      cfzct[i]->cd(j); hfpgazct[i*16+j]->Draw();
      cfcfd[i]->cd(j); hfpgacfd[i*16+j]->Draw();
    }
    crled[i]->Write();
    crzct[i]->Write();
    crcfd[i]->Write();
    cfled[i]->Write();
    cfzct[i]->Write();
    cfcfd[i]->Write();
  }
  for (int i=1; i<fNTubes; i++) {
    hrawled[i]->Write();
    hrawzct[i]->Write();
    hrawcfd[i]->Write();
    hfpgaled[i]->Write();
    hfpgazct[i]->Write();
    hfpgacfd[i]->Write();
  }
}
	

void PndSttRawAna::WriteSumHisto()
{
  fInFile->cd("Histos");
  hsumled->Write();
}


void PndSttRawAna::FindDeltaT() 
{

  fTInfo.ResizeTo(fNTubes, 5);
  // [0] = tpeak; [1] = t0; [2] = tmax; [3] = DeltaT; [4] = noise
  for (int i=0; i<fNTubes; i++) {
    for (int j=0; j<3; j++)
      fTInfoBin[i][j] = 0.;
    for (int j=0; j<5; j++)
      fTInfo[i][j] = 0.;
  }

  Double_t maximum = 0;
  Int_t NBins = 0;

  // for each tube
  for (int i=1; i<fNTubes; i++) 
    {
      maximum         = hrawled[i]->GetMaximum();
      if (maximum!=0) {
	fTInfoBin[i][0] = hrawled[i]->GetMaximumBin();
	fTInfo[i][0]    = hrawled[i]->GetBinCenter(fTInfoBin[i][0]);
	Double_t thresh = 0.1;
	for (int j=fTInfoBin[i][0]; j>0; j--) {
	  if (hrawled[i]->GetBinContent(j)<thresh*fTInfo[i][0]) {
	    fTInfoBin[i][1] = j;
	    fTInfo[i][1] = hrawled[i]->GetBinCenter(j);
	    break;
	  }
	}
	for (int j=fTInfoBin[i][0]; j<101; j++) {
	  if (hrawled[i]->GetBinContent(j)<thresh*fTInfo[i][0]) {
	    fTInfoBin[i][2] = j;
	    fTInfo[i][2] = hrawled[i]->GetBinCenter(j);
	    break;
	  }
	}
      }
      fTInfo[i][3] = fTInfo[i][2] - fTInfo[i][1];
       if (fVerbose>0) {
	cout << "Tube " << i << ": tpeak= " << fTInfo[i][0] << ", t0 = " 
	     << fTInfo[i][1] << ", tmax = " << fTInfo[i][2] << ", Deltat = " 
	     << fTInfo[i][3] << endl;
      }

       // find noise level
       TF1 *f = new TF1("f", "pol0", fTInfoBin[i][2], 1000);
       hrawled[i]->Fit(f, "Rq");
       Double_t par;
       fTInfo[i][4] = f->GetParameter(par);

       NBins = fTInfoBin[i][2] - fTInfoBin[i][1] + 1;
       for (int j=1; j<=NBins; j++)
	 //	 hsumled->SetBinContent(j, hrawled[i]->GetBinContent(j+fTInfoBin[i][1])-fTInfo[i][4]);
	 hsumled->AddBinContent(j, hrawled[i]->GetBinContent(j+fTInfoBin[i][1])-fTInfo[i][4]);
    }
}


void PndSttRawAna::Calibration() 
{

  Double_t Ntot = hsumled->Integral();
  Int_t tmaxbin = 0;
  for (int i=100; i>0; i--) {
    if (hsumled->GetBinContent(i) != 0) {
      tmaxbin = i;
      break;
    }
  }
  Double_t t[tmaxbin+1], rt[tmaxbin+1];
  for (int bin=0; bin<=tmaxbin; bin++) {
    if (bin==0) t[bin] = 0;
    else 
      t[bin] = hsumled->GetBinCenter(bin);
    rt[bin] = hsumled->Integral(0, bin)/Ntot*(fRadiusS-fRadiusW) + fRadiusW;
  }

  TMultiGraph *mg = new TMultiGraph("calib", "calib");
  TGraph *hcalib = new TGraph(tmaxbin+1, t, rt);
  hcalib->SetMarkerStyle(5);
  hcalib->SetLineWidth(2);
  mg->Add(hcalib, "P");

//   TF1 *fitcal = new TF1("fitcal", "pol5", 0, t[tmaxbin]);
//   mg->Fit(fitcal, "Rq");
//   Double_t params[6];
//   fitcal->GetParameters(params);
//   cout << "Fit parameters are: " << params[0] << ", " << params[1] << ", " 
//        << params[2] << ", " << params[3] << ", " << params[4] << ", " 
//        << params[5] << endl;
  TF1 *fitcal = new TF1("fitcal", "pol4", 0, t[tmaxbin]);
  mg->Fit(fitcal, "Rq");
  Double_t params[5];
  fitcal->GetParameters(params);

  TString fname = "rt_fitparams.dat";
  ofstream fout;
  fout.open(fname, ios::out);
  fout << params[0] << "\t" << params[1] << "\t" << params[2] 
       << "\t" << params[3] << "\t" << params[4] << endl;
  fout.close();

  if (fVerbose>1)
    cout << "Fit parameters are: " << params[0] << ", " << params[1] << ", " 
	 << params[2] << ", " << params[3] << ", " << params[4] << endl;

  fInFile->cd("Histos");
  mg->Write();
}


void PndSttRawAna::WriteOffset() 
{

  TString fname = "toffset.dat";
  ofstream fout;
  fout.open(fname, ios::out);
  
  // output: tubeID   t0 offset [ns]   Deltat [ns]   noise
  for (int i=1; i<fNTubes; i++) 
    fout << i << "\t" << fTInfo[i][1] << "\t" << fTInfo[i][3] 
	 << "\t" << fTInfo[i][4] << endl;

  fout.close();

}


void PndSttRawAna::ClearHistos() 
{
  delete *hrawled0;
  delete *hrawzct0;
  delete *hrawcfd0;
  delete *hrawled;
  delete *hrawzct;
  delete *hrawcfd;
  delete *hfpgaled0;
  delete *hfpgazct0;
  delete *hfpgacfd0;
  delete *hfpgaled;
  delete *hfpgazct;
  delete *hfpgacfd;
  delete *crled;
  delete *crzct;
  delete *crcfd;
  delete *cfled;
  delete *cfzct;
  delete *cfcfd;

}

	
ClassImp(PndSttRawAna)