#include "TMath.h"

bool CosTheta_EventRate_MuTauNeutrino_AtmosMuon(){

  //options
 bool tau      = true;//this is true only with 80 strings but can be false with 80 strings
 int nstring   = 9;
 if(nstring==9)tau=false;
 bool integralrate = false;

 bool use_firstguess=false;//if this is false, mc truth cos theta is used
 bool plot_firstguess=false;//if this is false, plotted as mc truth cos theta

  //gSystem->Load("libweighting-module");
  /*style*/
  SetStyle();
  gROOT->ForceStyle();
  //minor option for now
  bool dowrite  = false;
  
  string indir = "/home/aya/IceCube/ehe/March06/RootReducedLeaves/";
  
  if(nstring==80){
    string file_slope2 = "PhysicsWeightedWiscE-2.MergedReduced.root";
    string file_slope1 = "PhysicsWeightedE-1.MergedReduced.root";
    if(tau)string file_slope1_tau = "JulietTauGZKMergedAll.reduced.root";
  }
  else if(nstring==9){
    string file_slope1 = "JulietWeighedWiscE1-S9.MergedReduced.root";
    string file_slope2 = "JulietWeighedWiscE2-S9.MergedReduced.root";
    if(tau)string file_slope1_tau = "";
  }
  /* Read in E^-1 file */
  const char *filenameE1 = (const *char)(indir+file_slope1).c_str();
  TFile *fileE1 = (TFile *) TFile::Open(filenameE1);
  TTree *E1Tree = (TTree *) fileE1->Get("JulietTree");
  
  /* Read in E^-2 file */
  const char *filenameE2 = (const *char)(indir+file_slope2).c_str();
  TFile *fileE2 = (TFile *) TFile::Open(filenameE2);
  TTree *E2Tree = (TTree *) fileE2->Get("JulietTree");
  
  if(tau){
    /* Read in E^-1 Tau file */
    const char *filenameE1Tau = (const *char)(indir+file_slope1_tau).c_str();
    TFile *fileE1Tau = (TFile *) TFile::Open(filenameE1Tau);
    TTree *E1TauTree = (TTree *) fileE1Tau->Get("JulietTree");    
  }

  cout << "*************************************************" << endl;
  cout <<"Energy E-1 file " << filenameE1 << " is " << E1Tree->GetEntries() << endl;
  cout <<"Energy E-2 file " << filenameE2 << " is " << E2Tree->GetEntries() << endl;
  if(tau)cout <<"Energy E-1 Tau file " << filenameE1Tau << " is " << E1TauTree->GetEntries() << endl;
  cout << "Number of string is " << nstring << endl; 
  cout << "*************************************************" << endl;

//Define lowest level cuts --- Our Base Cut
TCut minimum_dom           = "DetectorResponseEvent_.numberOfDomWithRecoHit_>4";
TCut minimum_npe           = "DetectorResponseEvent_.totalBestEstimatedNPE_>200"; 
TCut minimum_fgcostheta    = "firstguessTrack_.cosTheta_>-5"; 

TCut basecut               = (minimum_dom&&minimum_npe);
//Propagation matrix
TCut prop_gzkNeutrino      = "generatorWeight_.julietPropagationMatrixNeutrinoFlux_";
TCut prop_atmMuon          = "generatorWeight_.julietPropagationMatrixAtmMuonFlux_";

//normalization const. for energy we fixed energy bin size to log(Delta E) = 0.1  
TCut energy_bin_norm_slope1 = "(generatorWeight_.maxEnergyLog_-generatorWeight_.minEnergyLog_)/0.1";
TCut energy_bin_norm_slope2 = "(pow(10.0,generatorWeight_.minEnergyLog_*(1.0-generatorWeight_.powerLawIndex_))-pow(10.0,generatorWeight_.maxEnergyLog_*(1.0-generatorWeight_.powerLawIndex_)))*pow(generatorPrimary_.energy_,(generatorWeight_.powerLawIndex_-1.0))/(0.1*log(10.0))";
TCut entry_norm             = "1/Entries$";
TCut entry                  = "Entries$";//E1Tree->Draw(entry);

//normalization const for geometry

TCut area_norm  = "TMath::Pi()*generatorWeight_.injectionSurfaceR_*generatorWeight_.injectionSurfaceR_";// [m^2]
TCut angle_norm = "4.0*TMath::Pi()";// [sr]
TCut sec_year   = "365.0*24.0*60.0*60.0"; // [sec/year]
TCut cm2_m2     = "1.0e4"; //[cm^2/m^2]
TCut m2_km2     = "1.0e6"; //[m^2/km^2]
TCut km2_m2     = "1.0e-6"; //[km^2/m^2]

TCut signal_npe   = "(DetectorResponseEvent_.totalBestEstimatedNPE_ > 5.0e5)"; 

if(use_firstguess)TCut signal_angle = "(firstguessTrack_.cosTheta_<0.1 &&  DetectorResponseEvent_.totalBestEstimatedNPE_ > 1.0e5 )";
else TCut signal_angle = "(generatorPrimary_.cosTheta_<0.1 &&  DetectorResponseEvent_.totalBestEstimatedNPE_ > 1.0e5 )";

//normalized neutrino and atmospheric muonw weights
TCut signalcut = signal_npe||signal_angle; 
TCut effArea_slope1 = energy_bin_norm_slope1*entry_norm*area_norm*km2_m2*angle_norm;
TCut effArea_slope2 = energy_bin_norm_slope2*entry_norm*area_norm*km2_m2*angle_norm;

TCut rate_neutrino_slope1 = effArea_slope1*prop_gzkNeutrino*m2_km2*cm2_m2*sec_year;
TCut rate_neutrino_slope2 = effArea_slope2*prop_gzkNeutrino*m2_km2*cm2_m2*sec_year;
TCut rate_atmMuon_slope1  = effArea_slope1*prop_atmMuon*m2_km2*cm2_m2*sec_year;
TCut rate_atmMuon_slope2  = effArea_slope2*prop_atmMuon*m2_km2*cm2_m2*sec_year;


 //Define atmospherlic muon distribution 
 //y-axis
 double rate_min       =  0.0;
 double rate_max       =  0.5;
 
 //x-axis
 double costheta_min_log    =  -1.;
 double costheta_max_log    =   1.;
 double costheta_bin_size   = 0.05;
 int    costheta_nbin       = (int)(costheta_max_log - costheta_min_log)/costheta_bin_size;


 TF1 *highnpe = new TF1("highnpe","log10(5*pow(10,5))",0,100);highnpe->SetLineWidth(1);highnpe->SetLineStyle(1);
 TF1 *horizon = new TF1("horizon","5",0,100);horizon->SetLineWidth(1);horizon->SetLineStyle(2);

 TH1D* hist_E1_EventRate = new TH1D("hist_E1_EventRate","Event Rate with signal domain cut",  costheta_nbin, costheta_min_log, costheta_max_log);
 TH1D* hist_E2_EventRate = new TH1D("hist_E2_EventRate","atmospheric Muon Event Rate from E^{-2}, func. of Primary Muon Energy", costheta_nbin, costheta_min_log, costheta_max_log);

 if(tau){ TH1D* hist_E1Tau_EventRate = new TH1D("hist_E1Tau_EventRate","Event Rate from E^{-1}, func. of Primary Tau Energy",  costheta_nbin, costheta_min_log, costheta_max_log); hist_E1Tau_EventRate->SetLineColor(kBlack);}
 hist_E1_EventRate->SetLineColor(kBlue);
 hist_E2_EventRate->SetLineColor(kRed);


 TAxis  *xaxis, *yaxis;
 xaxis = hist_E1_EventRate->GetXaxis();
 yaxis = hist_E1_EventRate->GetYaxis();
 yaxis->SetRangeUser(rate_min, rate_max); 
 if(plot_firstguess)xaxis->SetTitle("first-guess cos(#theta)"); 
 else xaxis->SetTitle("Primary MC Truth cos(#theta)"); 
 yaxis->SetTitle("Event Rate [year^{-1}]");

 xaxis = hist_E2_EventRate->GetXaxis();
 yaxis = hist_E2_EventRate->GetYaxis();
 yaxis->SetRangeUser(rate_min, rate_max);
 if(plot_firstguess)xaxis->SetTitle("first-guess cos(#theta)"); 
 else xaxis->SetTitle("Primary MC Truth cos(#theta)"); 
 yaxis->SetTitle("Event Rate [year^{-1}]");

 if(tau){
 xaxis = hist_E1Tau_EventRate->GetXaxis();
 yaxis = hist_E1Tau_EventRate->GetYaxis();
 yaxis->SetRangeUser(rate_min, rate_max); 
 if(plot_firstguess)xaxis->SetTitle("first-guess cos(#theta)"); 
 else xaxis->SetTitle("Primary MC Truth cos(#theta)"); 
 yaxis->SetTitle("Event Rate [year^{-1}]");
 }
 
 //Start Drawing

 TCanvas *canvas_EventRate = new TCanvas("EventRate","EventRate", 900,500);
 canvas_EventRate->Divide(2, 0.001,0.001);  
 canvas_EventRate->SetHighLightColor(10);
 //(basecut&&contained)*neutrino_slope
 /////////////////
 canvas_EventRate->cd(1);
 SetgPad();
 /////////////////
 E1Tree->SetLineColor(kBlue);
 E2Tree->SetLineColor(kRed);
 if(plot_firstguess){
 E2Tree->Draw("firstguessTrack_.cosTheta_", (basecut&&minimum_fgcostheta)*rate_atmMuon_slope2, "same");
 TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp"); 
 h1->Clear();
 h1->SetTitle("Event Rate without signal cut");
 TAxis  *xaxis = htemp->GetXaxis();
 TAxis  *yaxis = htemp->GetYaxis();
 xaxis->SetTitle("First-guess Cos(#theta)");
 yaxis->SetTitle("Event Rate [year^{-1}] without signal cut");
 yaxis->SetRangeUser(0.001, 350000);
 E2Tree->Draw("firstguessTrack_.cosTheta_", (basecut&&minimum_fgcostheta)*rate_atmMuon_slope2, "same");

 E1Tree->Draw("firstguessTrack_.cosTheta_", (basecut&&minimum_fgcostheta)*rate_neutrino_slope1, "same");
 if(tau)E1TauTree->Draw("firstguessTrack_.cosTheta_", (basecut&&minimum_fgcostheta)*rate_neutrino_slope1, "same");


 }else{
 E2Tree->Draw("generatorPrimary_.cosTheta_", (basecut)*rate_atmMuon_slope2, "");
 TH1D   *h1 = (TH1D*)gPad->GetPrimitive("htemp");
 h1->SetTitle("Event Rate without signal cut"); 
 TAxis  *xaxis = htemp->GetXaxis();
 TAxis  *yaxis = htemp->GetYaxis();
 yaxis->SetTitle("Event Rate [year^{-1}] without signal cut");
 xaxis->SetTitle("MC Truth Cos(#theta)");

 E2Tree->Draw("generatorPrimary_.cosTheta_", (basecut)*rate_atmMuon_slope2, "same");
 E1Tree->Draw("generatorPrimary_.cosTheta_", (basecut)*rate_neutrino_slope1, "same");
 if(tau)E1TauTree->Draw("generatorPrimary_.cosTheta_", (basecut)*rate_neutrino_slope1, "same"); }

 ///////////////// 
 canvas_EventRate->cd(2);
 SetgPad();
 /////////////////
 if(plot_firstguess){
 E1Tree->Draw("firstguessTrack_.cosTheta_ >> hist_E1_EventRate", (basecut&&signalcut&&minimum_fgcostheta)*rate_neutrino_slope1, "");
 E2Tree->Draw("firstguessTrack_.cosTheta_ >> hist_E2_EventRate", (basecut&&signalcut&&minimum_fgcostheta)*rate_atmMuon_slope2, "same");
 if(nstring==80&&tau){   
   E1TauTree->Draw("firstguessTrack_.cosTheta_ >> hist_E1Tau_EventRate", (basecut&&signalcut&&minimum_fgcostheta)*rate_neutrino_slope1, "same");
   }
 }else{
 E1Tree->Draw("generatorPrimary_.cosTheta_ >> hist_E1_EventRate", (basecut&&signalcut)*rate_neutrino_slope1, "");
 E2Tree->Draw("generatorPrimary_.cosTheta_ >> hist_E2_EventRate", (basecut&&signalcut)*rate_atmMuon_slope2, "same");
 if(nstring==80&&tau){   
   E1TauTree->Draw("generatorPrimary_.cosTheta_ >> hist_E1Tau_EventRate", (basecut&&signalcut)*rate_neutrino_slope1, "same");
   }
 }
/////////////////////////////////////////////////////////////////////////////////////////
if(integralrate)
   {
  TH1D* hist_E1_Total_EventRate = new TH1D("hist_E1_Total_EventRate","Event Rate, func. of Threshold Energy",  costheta_nbin, costheta_min_log, costheta_max_log);    
  TH1D* hist_E2_Total_EventRate = new TH1D("hist_E2_Total_EventRate","Event Rate, func. of Threshold Energy",  costheta_nbin, costheta_min_log, costheta_max_log);    
  TH1D* hist_E1Tau_Total_EventRate = new TH1D("hist_E1Tau_Total_EventRate","Event Rate, func. of Threshold Energy",  costheta_nbin, costheta_min_log, costheta_max_log);    

 
  xaxis = hist_E1_Total_EventRate->GetXaxis();
  yaxis = hist_E1_Total_EventRate->GetYaxis();
  xaxis->SetTitle("Primary cos(#theta)");  
  yaxis->SetTitle("Event Rate [year^{-1}]");
  xaxis = hist_E2_Total_EventRate->GetXaxis();
  yaxis = hist_E2_Total_EventRate->GetYaxis();
  xaxis->SetTitle("Primary cos(#theta)");  
  yaxis->SetTitle("Event Rate [year^{-1}]");
  xaxis = hist_E1Tau_Total_EventRate->GetXaxis();
  yaxis = hist_E1Tau_Total_EventRate->GetYaxis();
  xaxis->SetTitle("Primary cos(#theta)");  
  yaxis->SetTitle("Event Rate [year^{-1}]");


  double e1_total_rate=0.0, e2_total_rate=0.0, e1tau_total_rate=0.0;

  for(int i=costheta_nbin; i>=0; --i)    
    {
      e1_total_rate += hist_E1_EventRate->GetBinContent(i);
      hist_E1_Total_EventRate->SetBinContent(i, e1_total_rate);

      e2_total_rate += hist_E2_EventRate->GetBinContent(i);
      hist_E2_Total_EventRate->SetBinContent(i, e2_total_rate);

    }

  if(nstring==80&&tau){   

  for(int i=costheta_nbin; i>=0; --i)
    {
     	  e1tau_total_rate += hist_E1Tau_EventRate->GetBinContent(i);
	  hist_E1Tau_Total_EventRate->SetBinContent(i, e1tau_total_rate);
    }
  }

  if(nstring==80&&tau)TCanvas *canvas_EventRate_total = new TCanvas("EventRate_total","EventRate_total", 1200,400);
  else TCanvas *canvas_EventRate_total = new TCanvas("EventRate_total","EventRate_total", 900,500);
  if(nstring==80&&tau)canvas_EventRate_total->Divide(3, 0.001,0.001);
  else canvas_EventRate_total->Divide(2, 0.001,0.001);  
  canvas_EventRate_total->SetHighLightColor(10);


  canvas_EventRate_total->cd(1);
  SetgPad();
  /////////////////
  hist_E1_Total_EventRate->Draw();


  canvas_EventRate_total->cd(2);
  SetgPad();
  /////////////////
  hist_E2_Total_EventRate->Draw();

  if(nstring==80&&tau){
    canvas_EventRate_total->cd(3);
    SetgPad();
    /////////////////
    hist_E1Tau_Total_EventRate->Draw();
  }


   }
/////////////////////////////////////////////////////////////////////////////////////////


 //// output to root file if you want to ///////
 if(dowrite)
   {
     string outfile = "outhist.root";
     TFile *hist_outfile = new TFile(outfile.c_str(),"recreate");

     hist_outfile->Close();
   }
 cout << "exit" << endl;

}
////////////////////////////////
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);
}