// constants
const double syncTimeBase = 47.163878; // [ns]

// ntuple variables
int scrod = 0;                                        
int channel = 0;
int window = 0;
int pixel = 0;
float tdc = 0;
float tdc0 = 0;
float adc = 0;
float width = 0;
float integral = 0;
int flag = 0;

// function prototypes
double corTime(double tdc, int window, TH1F* tcor);
string toString(int n);


// -- code ---------------------------------------------------------

/**
 * @param scrodID scrod ID
 * @param fileName file containing the ntuple
 * @param dir directory with time base corrections (output of tbcFit.C)
 */
void tbcCheck(int scrodID, std::string fileName, std::string dir = "tbc")
{

  std::vector<TH1F*> corrections(128, 0);

  int nn = 0;
  std::string fName = dir + "/tbcScrod" + toString(scrodID) + ".root";
  TFile* fcor = TFile::Open(fName.c_str());
  if(!fcor) return;
  for(int k = 0; k < 16; k++) {
    int chan = k * 8 + 7;
    std::string hname = "tcor_ch" + toString(chan);
    TH1F* h = (TH1F*) fcor->Get(hname.c_str());
    if(!h) continue;
    for(int i = 0; i < 8; i++) corrections[chan - i] = h;
    nn++;
  }
  if(nn == 0) {
    cout << "No correctios available for scrod "<< scrodID << endl;
    return;
  }

  double mcTime[128];
  for(int i = 0; i < 128; i++) mcTime[i] = 0;
  TFile* fmc = TFile::Open("laserChannelTimesMC.root");
  if(fmc) {
    TProfile* profTime = (TProfile*) fmc->Get("profTime");
    for(int i = 0; i < 128; i++) mcTime[i] = profTime->GetBinContent(i+1);
  } else {
    cout << "No T0 from MC - file not found" << endl;
  }
  TFile* f = TFile::Open(fileName.c_str());
  if(!f) return;
  TTree* tree = (TTree*) f->Get("tree");

  tree->SetBranchAddress("scrod", &scrod);
  tree->SetBranchAddress("channel", &channel);
  tree->SetBranchAddress("window", &window);
  tree->SetBranchAddress("pixel", &pixel);
  tree->SetBranchAddress("tdc", &tdc);
  tree->SetBranchAddress("tdc0", &tdc0);
  tree->SetBranchAddress("adc", &adc);
  tree->SetBranchAddress("width", &width);
  tree->SetBranchAddress("integral", &integral);
  tree->SetBranchAddress("flag", &flag);

  TH2F* h0 = new TH2F("h0", "laser signal (uncorrected time base)",
		      128, 0, 128, 200, -5, 5);
  TH2F* h1 = new TH2F("h1", "laser signal (time base corrected)",
		      128, 0, 128, 200, -5, 5);
  TH2F* h2 = new TH2F("h2", "laser signal (mean corrected)",
		      128, 0, 128, 200, -5, 5);
  TH2F* h3 = new TH2F("h3", "laser signal (T0 corrected)",
		      128, 0, 128, 200, -5, 5);
  TH2F* h4 = new TH2F("h4", "laser signal (T0 corrected)",
		      128, 0, 128, 1000, -5, 100);
  h0->GetXaxis()->SetTitle("channel number");
  h0->GetYaxis()->SetTitle("time [ns]");
  h1->GetXaxis()->SetTitle("channel number");
  h1->GetYaxis()->SetTitle("time [ns]");
  h2->GetXaxis()->SetTitle("channel number");
  h2->GetYaxis()->SetTitle("time [ns]");
  h3->GetXaxis()->SetTitle("channel number");
  h3->GetYaxis()->SetTitle("time [ns]");
  h4->GetXaxis()->SetTitle("channel number");
  h4->GetYaxis()->SetTitle("time [ns]");
  
  int nev = (int)tree->GetEntries();

  for(int iev = 0; iev < nev; iev++) {
    tree->GetEntry(iev);
    if(flag == 2000 && scrod == scrodID) {
      TH1F* tcor = corrections[channel];
      if(!tcor) continue;
      double t = (tdc - tdc0) * syncTimeBase / 128 + 100;
      h0->Fill(channel, t);
      double time =  corTime(tdc, window, tcor) - corTime(tdc0, window, tcor) + 100;
      h1->Fill(channel, time);
      h2->Fill(channel, time);
    }
  }    

  TProfile* h2_pfx = h2->ProfileX();
  double t0[128];
  for(int i = 0; i < 128; i++) t0[i] = h2_pfx->GetBinContent(i + 1);

  for(int iter = 0; iter < 4; iter++) {
    h2->Reset();
    for(int iev = 0; iev < nev; iev++) {
      tree->GetEntry(iev);
      if(flag == 2000 && scrod == scrodID) {
	TH1F* tcor = corrections[channel];
	double time =  corTime(tdc, window, tcor) - corTime(tdc0, window, tcor) + 100;
	h2->Fill(channel, time - t0[channel]);
      }
    }    
    delete h2_pfx;
    h2_pfx = h2->ProfileX();
    for(int i = 0; i < 128; i++) t0[i] += h2_pfx->GetBinContent(i + 1);
  }  

  TProfile* prof = new TProfile("prof", "laser signal average", 128, 0, 128, -10, 10);
  prof->GetXaxis()->SetTitle("channel number");
  prof->GetYaxis()->SetTitle("time [ns]");
  for(int iev = 0; iev < nev; iev++) {
    tree->GetEntry(iev);
    if(flag == 2000 && scrod == scrodID) {
      TH1F* tcor = corrections[channel];
      double time =  corTime(tdc, window, tcor) - corTime(tdc0, window, tcor) + 100;
      if(time - t0[channel] < 5) prof->Fill(channel, time);
      h3->Fill(channel, time - t0[channel] + mcTime[channel]);
      h4->Fill(channel, time - t0[channel] + mcTime[channel]);
    }
  }    
  
  
}


string toString(int n)
{
  stringstream ss;
  ss << n;
  string str;
  ss >> str;
  return str;
}


double corTime(double tdc, int window, TH1F* tcor)
{

  if(!tcor) return 0;
  
  int sampleNum = int(tdc);
  double frac = tdc - sampleNum;

  sampleNum += window * 64; // counted from window 0
  int n = sampleNum / 256;
  int k = sampleNum % 256;
  double time = n * syncTimeBase * 2 + tcor->GetBinContent(k+1); // from sample 0 window 0
  time += (tcor->GetBinContent(k+2) - tcor->GetBinContent(k+1)) * frac; // add fraction
  
  return time;
  
}