/******************************************************************************
 --- PHD.C ---
 
 Calculate pulse height distribution (PHD) of output of Aoki shaping amp.
 *****************************************************************************/

#include "TTree.h"
#include "TH1D.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TText.h"

void PHD()
{
  const int kChN = 4;
  const int kCh[kChN] = {2, 3, 5, 8};
  const double kScale[kChN] = {.5, .5, .5, .5};
  /*
  const int kChN = 7;
  const int kCh[kChN] = {2, 3, 5, 6, 7, 8, 11};
  // Modify binning mistake
  const double kScale[kChN] = {.5, .5, 1, .5, .5, 1., 1};
  */
  TFile* file[kChN];
  TTree* tree[kChN];

  TH1F* PHD[kChN];

  for(int i=0; i<kChN; i++){
    file[i] = new TFile(Form("ch%d.root", kCh[i]), "read");
    tree[i] = new TTree(Form("ch%d", kCh[i]), Form("ch%d", kCh[i]));
    
    double pedestal, height, charge;
    int time;
    tree[i]->Branch("pedestal", &pedestal, "pedestal/D");
    tree[i]->Branch("height",   &height,   "height/D");
    tree[i]->Branch("charge",   &charge,   "charge/D");
    tree[i]->Branch("time",     &time, "time/I");

    for(int j=0; j<10000; j++){
      TTree* t;
      file[i]->GetObject("tree", t);

      TH1D* hist;
      t->SetBranchAddress("curve", &hist);
      t->SetBranchAddress("tv_sec", &time);
      t->GetEntry(j);

      pedestal = hist->Integral(1, 1000)/1000.;
      height = hist->GetMinimum()*kScale[i] - pedestal;
      charge = hist->Integral(1001, 1500)/50.*hist->GetBinWidth(1);
           
      delete hist;
      hist = 0;
      
      tree[i]->Fill();
    } // j
    tree[i]->Draw(Form("-height*1000.>>PHD%d", i));
    PHD[i] = (TH1F*)gDirectory->Get(Form("PHD%d", i));
  } // i
  
  TCanvas* can1 = new TCanvas("can1", "can1", 808, 828);
  can1->Divide(2, 4, 0.005, 0.01);
  
  TCanvas* can2 = new TCanvas("can2", "can2", 808, 828);
  can2->Divide(2, 4, 0.005, 0.01);

  TF1* f[kChN];
  for(int i=0; i<kChN; i++){
    f[i] = new TF1(Form("f%d", i), "[0]*exp(-x/[1])", 2., 20.);
    f[i]->SetParameters(4e2, 5);

    can1->cd(i+1);
    gPad->SetLogy();
    PHD[i]->SetTitle(Form("ch%d;Voltage [mV];Entries", kCh[i]));
    PHD[i]->GetXaxis()->SetRangeUser(0, 30);
    PHD[i]->SetMaximum(1e4);
    PHD[i]->Draw();
    PHD[i]->Fit(Form("f%d", i), "0", "", 2, 20);
    f[i]->SetLineColor(2);
    f[i]->Draw("same");
    double mean = f[i]->GetParameter(1);

    TText* txt = new TText(10, 1e3,
			   Form("Single Pulse Height = %1.2lf [mV]", mean));
    txt->Draw();

    can2->cd(i+1);
    tree[i]->Draw("-height:charge");
  } // i

  return;
}