#include "TH1F.h"
#include <iostream>
#include "TChain.h"
#include "TStyle.h"
#include "TFile.h"
#include "TRandom3.h"
#include "TF1.h"

using namespace::std;

void aligncalo_60hr(std::string inputFileName){

  gStyle->SetOptFit(1111);
  gStyle->SetOptStat(10);
  TFile *file  = new TFile(inputFileName.c_str());
  TTree *tree = (TTree*)file->Get("FastRotation/frTree");

  int64_t nEntries = tree->GetEntries();
  cout << "\n\tEntries = " << nEntries << ";\n" << endl;

  float time;
  float energy;
  int caloNum;
  tree->SetBranchAddress("time",&time);
  tree->SetBranchAddress("caloNum",&caloNum);
  tree->SetBranchAddress("energy",&energy);

  double timeShift = 0; // unit: time ticks, 1  = 1.25ns
  double tSignalLo = 0.0; // [us]
  double tSignalHi = 400 + tSignalLo; // [us]
	
  double magicTime = 0.149143; // [us]
  double magic = magicTime/24;
 
  TH1D *hSignalSum = new TH1D("hSignalSum","Decay Positron on all Calos;Time [#mus];Entries",400050,0.0,400.05);
  TH1D *hSignalRatioNum = new TH1D("hSignalRatioNum","Numerator (Half of decays);Time [#mus];Entries",400050,0.0,400.05);
  TH1D *hSignalRatioDen = new TH1D("hSignalRatioDen","Denominator (Smeared Times);Time [#mus];Entries",400050,0.0,400.05);
  TH1D *hSignalRatio = new TH1D("hSignalRatio","FR Signal - Ratio;Time [#mus];Entries",400050,0.0,400.05);
  TH1D *hSignalWeighted = new TH1D("hSignalWeighted","FR Signal - Background Removed;Time [#mus];Entries",400050,0.0,400.05);

  double targetPerc = 0;
  TRandom3* rng = new TRandom3(12345);

  int count = 0;
  for(int64_t entry = 0; entry < nEntries; entry++){

    if(100*float(entry) / nEntries > targetPerc){
      cout << "Processed " << 100*float(entry) / nEntries  << "%" << endl;
      targetPerc += 1;
    }

    tree->GetEntry(entry);
    double shiftedTime = (time - timeShift)*1.25 / 1000 - caloNum*magic;

    if(energy > 1500 and shiftedTime < tSignalHi and shiftedTime > tSignalLo){

      hSignalSum->Fill(shiftedTime);

      // Ratio histograms
      if (rng->Uniform() < 0.5){
	hSignalRatioNum->Fill(shiftedTime);
      } else {
	hSignalRatioDen->Fill(shiftedTime + (rng->Uniform()-0.5)*magicTime);
      }
    }
  }

  //////////////////////////////////////////////////////////////////////////////////////////////
  // Fit wiggle and remove from weighted histogram
  //////////////////////////////////////////////////////////////////////////////////////////////

  // Bin out fast rotation (incoming bins were 1 ns, so rebin by 150)
  TH1D* hSignalWiggle = (TH1D*)hSignalSum->Clone("hSignalWiggle");
  hSignalWiggle->Rebin(150);

  // Fit function (normal 5 parameter wiggle)
  TF1 *fitWiggle = new TF1("fitWiggle", "[0]*exp(-x/[1])*(1 + [2]*cos([3]*x + [4]))", 30.0, 95.0);
  fitWiggle->SetParameters(1, 64.4, 0.4, 1.4, 0);

  // Choose good starting guess by comparing integral of function and histogram
  double normGuess = hSignalWiggle->Integral(hSignalWiggle->GetXaxis()->FindBin(30), hSignalWiggle->GetXaxis()->FindBin(95), "WIDTH") / fitWiggle->Integral(30,95);
  fitWiggle->SetParameter(0,normGuess);

  // Fit function
  hSignalWiggle->Fit(fitWiggle,"NREML");

  // Set normalisation to 1 to just keep function without normalisation
  TF1* fitWiggleNorm = (TF1*)fitWiggle->Clone("fitWiggleNorm");
  fitWiggleNorm->SetParameter(0,1);

  // Make histogram divided by wiggle and lifetime
  for(int j = 1; j <= hSignalSum->GetNbinsX(); j++) {
    hSignalWeighted->SetBinContent(j, hSignalSum->GetBinContent(j) / fitWiggleNorm->Eval(hSignalSum->GetBinCenter(j)) );
    hSignalWeighted->SetBinError(j, hSignalSum->GetBinError(j) / fitWiggleNorm->Eval(hSignalSum->GetBinCenter(j)) );
  }
  hSignalWeighted->ResetStats();

  //////////////////////////////////////////////////////////////////////////////////////////////
  // Make Ratio Histogram
  //////////////////////////////////////////////////////////////////////////////////////////////

  // Scale denominator to 1 at t = 0 (same as fit function);
  hSignalRatioDen->Scale(fitWiggleNorm->Integral(hSignalRatioDen->GetXaxis()->GetXmin(), hSignalRatioDen->GetXaxis()->GetXmax()) / hSignalRatioDen->Integral("WIDTH")); 

  for(int iBin = 1; iBin <= hSignalRatioNum->GetNbinsX(); iBin++){
    if(hSignalRatioDen->GetBinContent(iBin) > 0){
      double U = hSignalRatioNum->GetBinContent(iBin);
      double V = hSignalRatioDen->GetBinContent(iBin);
      double UErr = hSignalRatioNum->GetBinError(iBin);
      double VErr = hSignalRatioDen->GetBinError(iBin);
      double ratio = U/V;
      double ratioErr = U > 0 ? ratio * sqrt(pow(UErr/U,2) + pow(VErr/V,2)) : 0;
      hSignalRatio->SetBinContent(iBin,ratio);
      hSignalRatio->SetBinError(iBin,ratioErr);
    }
  }
  hSignalRatio->ResetStats();


  TFile *fOutput = new TFile("AfterAlignCalo_60hr_Ratio.root", "RECREATE");
                               
  hSignalSum->Write();                  
  hSignalRatio->Write();
  hSignalRatioNum->Write();
  hSignalRatioDen->Write();
  hSignalWeighted->Write();

  fOutput->Close();

}