//  
//  trackQuality.C
//  <Enter discription here>
//
//  Created by Christian Elsasser on 31.05.13.
//  University of Zurich, elsasser@cern.ch
//  Copyright only for commercial use
//

//  General include
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <stdlib.h>
#include <limits>

#include <time.h>


//  ROOT include
#include "incBasicROOT.h"
#include "incDrawROOT.h"
#include "incHistROOT.h"
#include "incIOROOT.h"
#include "incRooFit.h"


#include "basicROOTTools.h"
#include "basicRooFitTools.h"
#include "TSpline.h"


#include "lhcbstyle.C"

#include "deplTool.h"




int trackQuality(int argc, char* argv[]);


int main(int argc, char* argv[]){
  
  TStyle* style = lhcbStyle();
  
  int argc_copy = argc;
  
  TApplication theApp("Analysis", &argc, argv);
  argc = theApp.Argc();
  argv = theApp.Argv();
  
  style->SetPadLeftMargin(0.12);
  
  int exit = trackQuality(argc,argv);
  
  Printf("End");
  theApp.Run(kTRUE); 
  
  return exit;
}


//  Executable method
int trackQuality(int argc, char* argv[]){
  
  
  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetTitleOffset(1.2,"Y");
  gStyle->SetTitleOffset(1.0,"X");

  const char* diskVar = std::getenv ("DISK");
  
  int fill      = 1944;
  TString det("IT");
  TString lay("T3X2");
  
  int caliStep  = 0;
  
  int val       = 7;
  
  
  Printf("=========================================");
  Info("trackQuality","Analyze: ");
  Printf("               Detector:  %s",det.Data());
  Printf("               Layer:     %s",lay.Data());
  Printf("               Fill:      %d",fill);
  
  TString dn_data(Form("%s/data/ST/Aging_Tuples/%s/%s"
                        ,diskVar,det.Data(),lay.Data()));
  TString fn_data(Form("%s/%d.root",
                       dn_data.Data(),fill));
  TString tn_sig(Form("STADCTrackMonitor/HitInfo/%s",
                      lay.Data()));
  TString tn_bkg(Form("STADCTrackMonitor/NoiseInfo/%s",
                      lay.Data()));
  
  

  TFile* f_data       = TFile::Open(fn_data.Data());
  
  if (!f_data) {
    Error("trackQuality","Data file does not exist!");
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  TTree* t_sig = (TTree*)f_data->Get(tn_sig.Data());
  if (!t_sig) {
    Error("trackQuality","Signal tree not found!");
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  TTree* t_bkg = (TTree*)f_data->Get(tn_bkg.Data());
  if (!t_bkg) {
    Error("trackQuality","Background tree not found!");
    Printf("=========================================\n");
    return EXIT_FAILURE;
  }
  
  
  int snCut = 12;
  
  TH2F* h_s = new TH2F("h_s","Signal histo",
                        100,0.0,5.0,
                        100,0.0,1.0);
  h_s->SetXTitle("#chi^{2}_{track}/ndf");
  h_s->SetYTitle("Ghost Prob.");
  h_s->SetZTitle("#it{S/B}");
  h_s->Sumw2();
  
  TH2F* h_n = (TH2F*)h_s->Clone("h_n");
  h_n->SetTitle("Noise histo");
  TH2F* h_sn = (TH2F*)h_s->Clone("h_sn");
  h_sn->SetTitle("Signal/Noise histo");
  
  t_sig->Draw("GhostP:TrChi2/TrNDoF>>h_s",Form("val%d>%d",
                                               val,snCut),"goff");
  t_sig->Draw("GhostP:TrChi2/TrNDoF>>h_n",Form("val%d<%d",
                                               val,snCut),"goff");
  
  h_sn->Divide(h_s,h_n,1.0,1.0);
  h_sn->SetMaximum(10.0);
  h_sn->SetMinimum(-1e-9);
  
  
  
  
  TH1F* h_tot_s = new TH1F("h_tot_s","All Signal Histo",
                           40,0.0,100.0);
  h_tot_s->SetXTitle("#it{r} [cm]");
  h_tot_s->SetYTitle("(Rejection) Efficiency");
  h_tot_s->SetLineColor(kRed);
  h_tot_s->SetMarkerColor(kRed);
  h_tot_s->SetFillColor(kRed-7);
  h_tot_s->Sumw2();
  
  
  TH1F* h_sel_s = (TH1F*)h_tot_s->Clone("h_sel_s");
  h_sel_s->SetTitle("Sel Signal Histo");
  TH1F* h_eff_s = (TH1F*)h_tot_s->Clone("h_eff_s");
  h_eff_s->SetTitle("Eff Signal Histo");
  
  TH1F* h_tot_n = (TH1F*)h_tot_s->Clone("h_tot_n");
  h_tot_n->SetTitle("Sel Noise Histo");
  h_tot_n->SetLineColor(kBlue);
  h_tot_n->SetMarkerColor(kBlue);
  h_tot_n->SetFillColor(kBlue-7);
  
  TH1F* h_sel_n = (TH1F*)h_tot_n->Clone("h_sel_n");
  h_sel_n->SetTitle("Sel. Noise Histo");
  TH1F* h_pur_n = (TH1F*)h_tot_n->Clone("h_eff_n");
  h_pur_n->SetTitle("Pur. Noise Histo");
  
  TCut c_sel = STTool::cu_track_TT;
  if (strcmp(det.Data(),"IT")==0) {
    c_sel = STTool::cu_track_IT;
  }
  
  t_sig->Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_tot_s",
              Form("val%d>%d",val,snCut),"goff");
  t_sig->Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_tot_n",
              Form("val%d<%d",val,snCut),"goff");
  
  t_sig->Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_sel_s",
              Form("(val%d>%d) && (%s)",val,snCut,c_sel.GetTitle()),"goff");
  t_sig->Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_sel_n",
              Form("(val%d<%d) && !(%s)",val,snCut,c_sel.GetTitle()),"goff");
  
  h_eff_s->Divide(h_sel_s,h_tot_s,1.0,1.0,"b");
  h_pur_n->Divide(h_sel_n,h_tot_n,1.0,1.0,"b");
  
  TCanvas* c_eff = new TCanvas("c_eff","Eff./Pur. Canvas",800,600);
  h_eff_s->SetMaximum(1.1);
  h_eff_s->SetMinimum(-0.1);
  h_pur_n->SetMaximum(1.1);
  h_pur_n->SetMinimum(-0.1);
  h_eff_s->Draw();
  h_pur_n->Draw("same");
  
  Info("trackQuality","Eff/Pur");
  Printf("         Radius      : Eff.    /Pur.");
  for (int i=1; i<=h_eff_s->GetNbinsX() && i<=h_pur_n->GetNbinsX(); i++) {
    Printf("         r = %8.4f cm: %8.4f/%8.4f",
           h_eff_s->GetBinCenter(i),h_eff_s->GetBinContent(i),h_pur_n->GetBinContent(i));
  }
  
  
  gStyle->SetPadLeftMargin(0.20);
  gStyle->SetPadRightMargin(0.25);
  gStyle->SetPadBottomMargin(0.15);
  gStyle->SetPadTopMargin(0.05);
  
  TCanvas* c_sn = new TCanvas("c_sn","Signal/Noise Canvas",800,600);
  h_sn->Draw("colz");
  
  TCanvas* c_sn_cut = new TCanvas("c_sn_cut","Signal/Noise Canvas with cuts",800,600);
  h_sn->Draw("colz");
  
  TLine* line = new TLine();
  line->SetLineWidth(4);
  line->SetLineColor(kWhite);
  if (strcmp(det.Data(),"TT")==0) {
    line->DrawLine(0.0,0.01,1.6,0.16);
    line->DrawLine(3.0,0.07,1.6,0.16);
    line->DrawLine(3.0,0.07,3.0,0.0);
  }else{
    line->DrawLine(0.0,0.0,3.5,0.35);
    line->DrawLine(3.5,0.35,5.0,0.20);
  }
  h_sn->Draw("axissame");
  
  
  Printf("=========================================\n");
  return EXIT_SUCCESS;
}