#include "TMath.h"

bool ChannelWise_NPE_Time(){

  //options
 int nstring   = 80;
 bool doNpe=false; //compare mc truth npe and best estimated npe, atwd npe and fadc npe
 bool doPulseTime=true; //compare first mc hit time and larger charge time, atwd, fadc 
 bool doLETime=true; //compare first mc hit time and larger charge time, atwd, fadc 

  /*style*/
  SetStyle();
  gROOT->ForceStyle();
  //minor option for now
  bool dowrite  = false;
  
  string infile = "/home/aya/IceCube/ehe/March06/MergedRoot/ReWiscE1/JulietE1-0201-80strings-1.added_1.root";
  /* Read in E^-1 file */
  TFile *fileE1 = (TFile *) TFile::Open(infile.c_str());
  TTree *E1Tree = (TTree *) fileE1->Get("JulietTree");
  
  cout << "*************************************************" << endl;
  cout <<"Energy E-1 file " << infile << " is " << E1Tree->GetEntries() << endl;
  cout << "*************************************************" << endl;
  
  //Define lowest level cuts --- Our Base Cut for channels
  TCut minimum_dom           = "DetectorResponseEvent_.numberOfDomWithRecoHit_>4";
  TCut minimum_npe           = "DetectorResponseEvent_.totalBestEstimatedNPE_>200"; 
  TCut minimum_fgcostheta    = "firstguessTrack_.cosTheta_>-5";   

  TCut event_number ="eventNumber>3&&eventNumber<15";
  TCut icecube_string = "AtwdResponseChannels_.string_>0||FadcResponseChannels_.string_>0";
  TCut atwd_time      = "AtwdResponseChannels_.pulseTime_>-2000";
  TCut fadc_time      = "FadcResponseChannels_.pulseTime_>-2000";
  TCut mchit_time     = "McHitChannels_.firstMcHitTime_>-2000";
  
  TCut basecut = (icecube_string&&atwd_time&&fadc_time&&mchit_time&&event_number);
  
  TCut use_fadc = "PulseInfoChannels_.largerFadcCharge_>0";
  TCut use_atwd = "PulseInfoChannels_.largerAtwdCharge_>0";
  
  /*
AtwdResponseChannels_.integratedCharge_
AtwdResponseChannels_.pulseTime_
AtwdResponseChannels_.pulseAmplitude_
AtwdResponseChannels_.pulsePosition_
FadcResponseChannels_.estimatedNPE_
FadcResponseChannels_.pulsePeakTime_
FadcResponseChannels_.pulseTime50_
FadcResponseChannels_.pulseTime80_
FadcResponseChannels_.launchTime_
FadcResponseChannels_.leTime_
*/

  if(doPulseTime){
    double time_min = 0;
    double time_max = 12000; 
    //Start Drawing
    TCanvas *canvas_time = new TCanvas("canvas_time","canvas_time", 800,800);
    canvas_time->Divide(2,2, 0.001,0.001);  
    canvas_time->SetHighLightColor(10);
    /////////////////
    canvas_time->cd(1);
    SetgPad();
    /////////////////
    E1Tree->SetMarkerColor(kBlue);
    E1Tree->Draw("(AtwdResponseChannels_.pulseTime_):(McHitChannels_.firstMcHitTime_)", (basecut&&!use_atwd), "");
    E1Tree->SetMarkerColor(kRed);
    E1Tree->Draw("(AtwdResponseChannels_.pulseTime_):(McHitChannels_.firstMcHitTime_)", (basecut&&use_atwd), "same");

    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Atwd estimated time");
    xaxis->SetTitle("Mc first hit time");
    xaxis->SetRangeUser(time_min, time_max);yaxis->SetRangeUser(time_min, time_max);

    /////////////////
    canvas_time->cd(2);
    SetgPad();
    /////////////////
    E1Tree->SetMarkerColor(kRed);
    E1Tree->Draw("(FadcResponseChannels_.pulseTime_):(McHitChannels_.firstMcHitTime_)", (basecut&&!use_fadc), "");

    E1Tree->SetMarkerColor(kBlue);
    E1Tree->Draw("(FadcResponseChannels_.pulseTime_):(McHitChannels_.firstMcHitTime_)", (basecut&&use_fadc), "same");
    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Fadc estimated time");
    xaxis->SetTitle("Mc first hit time");
    xaxis->SetRangeUser(time_min, time_max);yaxis->SetRangeUser(time_min, time_max);

    /////////////////
    canvas_time->cd(3);
    SetgPad();
    /////////////////
    E1Tree->SetMarkerColor(kBlack);
    E1Tree->Draw("(AtwdResponseChannels_.pulseTime_*PulseInfoChannels_.largerAtwdCharge_+FadcResponseChannels_.pulseTime_*PulseInfoChannels_.largerFadcCharge_):(McHitChannels_.firstMcHitTime_)", (basecut), "");
    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Best estimated time");
    xaxis->SetTitle("Mc first hit time");
    xaxis->SetRangeUser(time_min, time_max);yaxis->SetRangeUser(time_min, time_max);


    canvas_time->cd(4);
    SetgPad();
    E1Tree->SetMarkerColor(kBlue);
    E1Tree->Draw("(FadcResponseChannels_.pulseTime_):(AtwdResponseChannels_.pulseTime_)", (basecut&&use_fadc), "");
    E1Tree->SetMarkerColor(kRed);
    E1Tree->Draw("(FadcResponseChannels_.pulseTime_):(AtwdResponseChannels_.pulseTime_)", (basecut&&use_atwd), "same");

    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Fadc estimated time");
    xaxis->SetTitle("Atwd estimated time");
    xaxis->SetRangeUser(time_min, time_max);yaxis->SetRangeUser(time_min, time_max);
  }//done doTime

  if(doLETime){
    double time_min = 0;
    double time_max = 12000; 
    //Start Drawing
    TCanvas *canvas_letime = new TCanvas("canvas_letime","canvas_letime", 800,800);
    canvas_letime->Divide(2,2, 0.001,0.001);  
    canvas_letime->SetHighLightColor(10);
    /////////////////
    canvas_letime->cd(1);
    SetgPad();
    /////////////////
    E1Tree->SetMarkerColor(kBlue);
    E1Tree->Draw("(AtwdResponseChannels_.leTime_):(McHitChannels_.firstMcHitTime_)", (basecut&&!use_atwd), "");
    E1Tree->SetMarkerColor(kRed);
    E1Tree->Draw("(AtwdResponseChannels_.leTime_):(McHitChannels_.firstMcHitTime_)", (basecut&&use_atwd), "same");

    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Atwd estimated time");
    xaxis->SetTitle("Mc first hit time");
    xaxis->SetRangeUser(time_min, time_max);yaxis->SetRangeUser(time_min, time_max);

    /////////////////
    canvas_letime->cd(2);
    SetgPad();
    /////////////////
    E1Tree->SetMarkerColor(kRed);
    E1Tree->Draw("(FadcResponseChannels_.leTime_):(McHitChannels_.firstMcHitTime_)", (basecut&&!use_fadc), "");

    E1Tree->SetMarkerColor(kBlue);
    E1Tree->Draw("(FadcResponseChannels_.leTime_):(McHitChannels_.firstMcHitTime_)", (basecut&&use_fadc), "same");
    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Fadc estimated time");
    xaxis->SetTitle("Mc first hit time");
    xaxis->SetRangeUser(time_min, time_max);yaxis->SetRangeUser(time_min, time_max);

    /////////////////
    canvas_letime->cd(3);
    SetgPad();
    /////////////////
    E1Tree->SetMarkerColor(kBlack);
    E1Tree->Draw("(AtwdResponseChannels_.leTime_*PulseInfoChannels_.largerAtwdCharge_+FadcResponseChannels_.leTime_*PulseInfoChannels_.largerFadcCharge_):(McHitChannels_.firstMcHitTime_)", (basecut), "");
    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Best estimated time");
    xaxis->SetTitle("Mc first hit time");
    xaxis->SetRangeUser(time_min, time_max);yaxis->SetRangeUser(time_min, time_max);


    canvas_letime->cd(4);
    SetgPad();
    E1Tree->SetMarkerColor(kBlue);
    E1Tree->Draw("(FadcResponseChannels_.leTime_):(AtwdResponseChannels_.leTime_)", (basecut&&use_fadc), "");
    E1Tree->SetMarkerColor(kRed);
    E1Tree->Draw("(FadcResponseChannels_.leTime_):(AtwdResponseChannels_.leTime_)", (basecut&&use_atwd), "same");

    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Fadc estimated time");
    xaxis->SetTitle("Atwd estimated time");
    xaxis->SetRangeUser(time_min, time_max);yaxis->SetRangeUser(time_min, time_max);
  }//done doLeTime


  //////////////////////////////////////////
  if(doNpe){
    double npe_min_log = 0;
    double npe_max_log = 4.5; 
    //Start Drawing
    TCanvas *canvas_npe = new TCanvas("canvas_npe","canvas_npe", 900,900);
    canvas_npe->Divide(2,2, 0.001,0.001);  
    canvas_npe->SetHighLightColor(10);
    /////////////////
    
    canvas_npe->cd(1);
    SetgPad();
    /////////////////
    E1Tree->SetMarkerColor(kBlue);
    E1Tree->Draw("log10(AtwdResponseChannels_.estimatedNPE_):log10(McHitChannels_.mcNpe_)", (basecut&&!use_atwd), "");
    E1Tree->SetMarkerColor(kRed);
    E1Tree->Draw("log10(AtwdResponseChannels_.estimatedNPE_):log10(McHitChannels_.mcNpe_)", (basecut&&use_atwd), "same");

    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Atwd estimated NPE");
    xaxis->SetTitle("Mc NPE");
    xaxis->SetRangeUser(npe_min_log, npe_max_log);yaxis->SetRangeUser(npe_min_log, npe_max_log);

    /////////////////
    canvas_npe->cd(2);
    SetgPad();
    /////////////////
    E1Tree->SetMarkerColor(kRed);
    E1Tree->Draw("log10(FadcResponseChannels_.estimatedNPE_):log10(McHitChannels_.mcNpe_)", (basecut&&!use_fadc), "");
    E1Tree->SetMarkerColor(kBlue);
    E1Tree->Draw("log10(FadcResponseChannels_.estimatedNPE_):log10(McHitChannels_.mcNpe_)", (basecut&&use_fadc), "same");
    E1Tree->SetMarkerColor(kBlack);
    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Fadc estimated NPE");
    xaxis->SetTitle("Mc NPE");
    xaxis->SetRangeUser(npe_min_log, npe_max_log);yaxis->SetRangeUser(npe_min_log, npe_max_log);

    /////////////////
    canvas_npe->cd(3);
    SetgPad();
    /////////////////
    E1Tree->Draw("log10(AtwdResponseChannels_.estimatedNPE_*PulseInfoChannels_.largerAtwdCharge_+FadcResponseChannels_.estimatedNPE_*PulseInfoChannels_.largerFadcCharge_):log10(McHitChannels_.mcNpe_)", (basecut), "");
    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Best estimated NPE");
    xaxis->SetTitle("Mc NPE");
    xaxis->SetRangeUser(npe_min_log, npe_max_log);yaxis->SetRangeUser(npe_min_log, npe_max_log);
    
    
    canvas_npe->cd(4);
    SetgPad();
    E1Tree->SetMarkerColor(kBlue);
    E1Tree->Draw("log10(FadcResponseChannels_.estimatedNPE_):log10(AtwdResponseChannels_.estimatedNPE_)", (basecut&&use_fadc), "");
    E1Tree->SetMarkerColor(kRed);
    E1Tree->Draw("log10(FadcResponseChannels_.estimatedNPE_):log10(AtwdResponseChannels_.estimatedNPE_)", (basecut&&use_atwd), "same");
 
    TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
    TAxis  *xaxis = htemp->GetXaxis();
    TAxis  *yaxis = htemp->GetYaxis();
    yaxis->SetTitle("Fadc estimated NPE");
    xaxis->SetTitle("Atwd estimated NPE");
    xaxis->SetRangeUser(npe_min_log, npe_max_log);yaxis->SetRangeUser(npe_min_log, npe_max_log);
    
  }//done doNpe

}//done
////////////////////////////////

////////////////////////////////
void SetStyle(){

   gStyle->SetOptStat(1111);

   gStyle->SetPadGridX(true);
   gStyle->SetPadGridY(true);

   gStyle->SetStatW(0.4);
   gStyle->SetStatFontSize(0.03);
   gStyle->SetTitleFontSize(18);
   gStyle->SetLabelSize(0.03,"x");
   gStyle->SetLabelSize(0.03,"y");
   gStyle->SetLabelSize(0.03,"z");
   gStyle->SetLabelOffset(0.002,"x");
   gStyle->SetNdivisions(1025,"x");
   gStyle->SetNdivisions(1010,"y");
   gStyle->SetTitleOffset(1.4,"y");
   gStyle->SetTitleOffset(0.8,"x");
   gStyle->SetOptStat(0000);
   
   gStyle->SetTitleBorderSize(0);
   gStyle->SetTitleW(0.6);
   gStyle->SetTitleH(0.06);
   gStyle->SetOptTitle(1);
   gStyle->SetPalette(1);
   
}
void SetgPad(){
  double right  = 0.12;
  //double right  = 0.05;
  double left   = 0.12;
  double top    = 0.1;
  double bottom = 0.08;
  gPad->SetRightMargin(right);
  gPad->SetLeftMargin(left);
  gPad->SetTopMargin(top);
  gPad->SetBottomMargin(bottom);
}