#include "TMath.h"

bool NPE_MuTauNeutrino_AtmosMuon(){

  bool dowrite  = false;
  bool doLog    = true;
  bool doWeight = true;//this is only shape observation not truly normalized
  bool tau      = false;//this is true only with 80 strings but can be false with 80 strings

  int nstring=9;

  //gSystem->Load("libweighting-module");
  /*style*/
  SetStyle();
  gROOT->ForceStyle();
  
  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 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);

//normalized neutrino and atmospheric muonw weights
TCut neutrino_slope1 = prop_gzkNeutrino*energy_bin_norm_slope1*entry_norm;
TCut neutrino_slope2 = prop_gzkNeutrino*energy_bin_norm_slope2*entry_norm;
TCut atmMuon_slope1  = prop_atmMuon*energy_bin_norm_slope1*entry_norm;
TCut atmMuon_slope2  = prop_atmMuon*energy_bin_norm_slope2*entry_norm;

//Define atmospherlic muon distribution 
//y-axis
 int    npe_nbin        =  100;
 if(doLog){
   double npe_min_log   =  1;
   double npe_max_log   =  8;
 }else{
   double npe_min       =  0;
   double npe_max       =  1e8;
 }
 //x-axis
 int    costheta_nbin   =  40;
 double costheta_min    = -1.0;
 double costheta_max    =  1.0;
 //z-axis
 double atmMuon_min     =  1.0e-22;
 double atmMuon_max     =  1.0e-15;
 double neutrino_min    =  1.0e-25;
 double neutrino_max    =  1.0e-21;

 if(doLog)TH1D* hist_E2_NPE = new TH1D("hist_E2_NPE","Estimated NPE from E^{-2} Mu sample", npe_nbin, npe_min_log, npe_max_log);
 else TH1D* hist_E2_NPE = new TH1D("hist_E2_NPE","Estimated NPE from E^{-2} Mu sample",  npe_nbin, npe_min, npe_max);
 if(doLog)TH1D* hist_E1_NPE = new TH1D("hist_E1_NPE","Estimated NPE from E^{-1} Mu sample",  npe_nbin, npe_min_log, npe_max_log);
 else TH1D* hist_E1_NPE = new TH1D("hist_E1_NPE","Estimated NPE from E^{-1} Mu sample",  npe_nbin, npe_min, npe_max);
 if(tau){
 if(doLog)TH1D* hist_E1Tau_NPE = new TH1D("hist_E1Tau_NPE","Estimated NPE from E^{-1} Tau sample",  npe_nbin, npe_min_log, npe_max_log);
 else TH1D* hist_E1Tau_NPE = new TH1D("hist_E1Tau_NPE","Estimated NPE from E^{-1} Tau sample",  npe_nbin, npe_min, npe_max);
 }

 if(doLog)TH1D* hist_AtmMuon_NPE = new TH1D("hist_AtmMuon_NPE","Atmospheric Muon weighed estimated NPE",npe_nbin, npe_min_log,npe_max_log);
 else TH1D* hist_AtmMuon_NPE = new TH1D("hist_AtmMuon_NPE","Atmospheric Muon weighed estimated NPE",npe_nbin, npe_min, npe_max);
 if(doLog)TH1D* hist_Neutrino_NPE = new TH1D("hist_Neutrino_NPE","GZK Mu Neutrino weighed estimated NPE",npe_nbin, npe_min_log, npe_max_log);
 else TH1D* hist_Neutrino_NPE = new TH1D("hist_Neutrino_NPE","GZK Mu Neutrino weighed estimated NPE",npe_nbin, npe_min, npe_max); 
 if(tau){
 if(doLog)TH1D* hist_TauNeutrino_NPE = new TH1D("hist_TauNeutrino_NPE","GZK Tau Neutrino weighed estimated NPE",npe_nbin, npe_min_log, npe_max_log);
 else TH1D* hist_TauNeutrino_NPE = new TH1D("hist_TauNeutrino_NPE","GZK Tau Neutrino weighed estimated NPE",npe_nbin, npe_min, npe_max);
 }

 TAxis  *xaxis, *yaxis;
 xaxis = hist_AtmMuon_NPE->GetXaxis();
 yaxis = hist_AtmMuon_NPE->GetYaxis();
 if(doLog)xaxis->SetTitle("log_{10} Estimated event-sum NPE");
 else xaxis->SetTitle("Estimated event-sum NPE");
 yaxis->SetTitle("Atmospheric Muon Flux"); 

 xaxis = hist_Neutrino_NPE->GetXaxis();
 yaxis = hist_Neutrino_NPE->GetYaxis();
 if(doLog)xaxis->SetTitle("log_{10} Estimated event-sum NPE");
 else xaxis->SetTitle("Estimated event-sum NPE");
 yaxis->SetTitle("GZK Mu Neutrino Flux"); 
 if(tau){
 xaxis = hist_TauNeutrino_NPE->GetXaxis();
 yaxis = hist_TauNeutrino_NPE->GetYaxis();
 if(doLog)xaxis->SetTitle("log_{10} Estimated event-sum NPE");
 else xaxis->SetTitle("Estimated event-sum NPE");
 yaxis->SetTitle("GZK Tau Neturino Flux"); 
 }
 xaxis = hist_E1_NPE->GetXaxis();
 yaxis = hist_E1_NPE->GetYaxis();
 if(doLog)xaxis->SetTitle("log_{10} Estimated event-sum NPE");
 else xaxis->SetTitle("Estimated event-sum NPE");
 yaxis->SetTitle("number of events"); 

 xaxis = hist_E2_NPE->GetXaxis();
 yaxis = hist_E2_NPE->GetYaxis();
 if(doLog)xaxis->SetTitle("log_{10} Estimated event-sum NPE");
 else xaxis->SetTitle("Estimated event-sum NPE");
 yaxis->SetTitle("number of events "); 
 if(tau){
 xaxis = hist_E1Tau_NPE->GetXaxis();
 yaxis = hist_E1Tau_NPE->GetYaxis();
 if(doLog)xaxis->SetTitle("log_{10} Estimated event-sum NPE");
 else xaxis->SetTitle("Estimated event-sum NPE");
 yaxis->SetTitle("number of events "); 
 }
 //Start Drawing
 if(nstring==80&&tau)TCanvas *canvas_NPE = new TCanvas("NPE","Npe", 1200,400);
 else TCanvas *canvas_NPE = new TCanvas("NPE","Npe", 900,500);
 if(nstring==80&&tau)canvas_NPE->Divide(3, 0.001,0.001);
 else canvas_NPE->Divide(2, 0.001,0.001);  
 canvas_NPE->SetHighLightColor(10);

 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
 canvas_NPE->cd(1);
 SetgPad();
 if(doWeight){
   if(doLog)E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_Neutrino_NPE",  basecut*neutrino_slope1, "");
   else E1Tree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_Neutrino_NPE", basecut*neutrino_slope1, "");
 }else{
      if(doLog)E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E1_NPE", basecut, "");
      else E1Tree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E1_NPE", basecut, "");
 }

 canvas_NPE->cd(2);
 SetgPad();
 if(doWeight){
   if(doLog) E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_AtmMuon_NPE", basecut*atmMuon_slope2, "");
   else E2Tree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_AtmMuon_NPE", basecut*atmMuon_slope2, "");
 }else{
   if(doLog) E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E2_NPE",basecut, "");

   else E2Tree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E2_NPE",basecut, "");
 }

 if(nstring==80&&tau){
   
   canvas_NPE->cd(3);
   SetgPad();
   if(doWeight){
     if(doLog)E1TauTree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_TauNeutrino_NPE",basecut*neutrino_slope1,"");
     else E1TauTree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_TauNeutrino_NPE", basecut*neutrino_slope1, "");
   }else{
     if(doLog)E1TauTree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E1Tau_NPE", basecut, "");
     else E1TauTree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E1Tau_NPE", basecut, "");
   }
   
 }
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////



 //// output to root file if you want to ///////
 if(dowrite)
   {
     string outfile = "outhist.root";
     TFile *hist_outfile = new TFile(outfile.c_str(),"recreate");
     hist_Neutrino_NPE->Write();
     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.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);
}