#include "TMath.h"

bool Nch_NPE_MuTauNeutrino_AtmosMuon(){

  bool dowrite=false;
  bool doLog=false;
  bool doAtm=true;
  bool tau = false;
  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 = "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 = "";
  }
/* 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");

cout << "*************************************************" << endl;
cout <<"Energy E-1 file " << filenameE1 << " is " << E1Tree->GetEntries() << endl;
cout <<"Energy E-2 file " << filenameE2 << " is " << E2Tree->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    nch_nbin        =  100;
 if(doLog){
   double nch_min_log   =  0;
   double nch_max_log   =  5;
 }else{
   double nch_min       =  0;
   if(nstring==80)double nch_max       =  80*60;
   if(nstring==9)double nch_max       =  10*60;
 }
 //x-axis
 int    npe_nbin   =  100;
 double npe_min    =  2;
 double npe_max    =  8;
 //z-axis
 double atmMuon_min     =  1.0e-22;
 double atmMuon_max     =  1.0e-15;
 double neutrino_min    =  1.0e-28;
 double neutrino_max    =  1.0e-21;

 if(doLog)TH2D* hist_AtmMuon_Nch_Npe = new TH2D("hist_AtmMuon_Nch_Npe","Atmospheric Muon from E^{-2} sample", 
					    nch_nbin, nch_min_log, nch_max_log, npe_nbin, npe_min, npe_max);
 else TH2D* hist_AtmMuon_Nch_Npe = new TH2D("hist_AtmMuon_Nch_Npe","Atmospheric Muon from E^{-2} sample", 
					    nch_nbin, nch_min, nch_max, npe_nbin, npe_min, npe_max);
 if(doLog)TH2D* hist_Neutrino_Nch_Npe = new TH2D("hist_Neutrino_Nch_Npe","GZK Neutrino from E^{-1} sample", 
					    nch_nbin, nch_min_log, nch_max_log, npe_nbin, npe_min, npe_max);
 else TH2D* hist_Neutrino_Nch_Npe = new TH2D("hist_Neutrino_Nch_Npe","GZK Neutrino from E^{-1} sample", 
						  nch_nbin, nch_min, nch_max, npe_nbin, npe_min, npe_max);
 TAxis  *xaxis, *yaxis, *zaxis;
 xaxis = hist_AtmMuon_Nch_Npe->GetXaxis();
 yaxis = hist_AtmMuon_Nch_Npe->GetYaxis();
 zaxis = hist_AtmMuon_Nch_Npe->GetZaxis();
 hist_AtmMuon_Nch_Npe->SetContour(40);
 if(doLog)xaxis->SetTitle("log_{10} Number of DOMs with Reco Pulse");
 else xaxis->SetTitle("Number of DOMs with Reco Pulse");
 yaxis->SetTitle("log_{10} NPE"); 
 zaxis->SetRangeUser(atmMuon_min, atmMuon_max); 

 xaxis = hist_Neutrino_Nch_Npe->GetXaxis();
 yaxis = hist_Neutrino_Nch_Npe->GetYaxis();
 zaxis = hist_Neutrino_Nch_Npe->GetZaxis();
 hist_Neutrino_Nch_Npe->SetContour(40);
 if(doLog)xaxis->SetTitle("log_{10} Number of DOMs with Reco Pulse");
 else xaxis->SetTitle("Number of DOMs with Reco Pulse");
 yaxis->SetTitle("log_{10} NPE"); 
 zaxis->SetRangeUser(neutrino_min, neutrino_max); 

 //Start Drawing
 if(doAtm)TCanvas *canvas_Nch_Npe = new TCanvas("Nch_Npe","Nch : NPE", 900,500);
 else TCanvas *canvas_Nch_Npe = new TCanvas("Nch_Npe","Nch : NPE", 500,500);
 if(doAtm)canvas_Nch_Npe->Divide(2, 0.001,0.001);  
 canvas_Nch_Npe->SetHighLightColor(10);
 TLine *nch53 = new TLine(53,2,53,8);nch53->SetLineWidth(3);nch53->SetLineStyle(1);nch53->SetLineColor(kBlack);
 TLine *npe4 = new TLine(0,4,600,4);npe4->SetLineWidth(3);npe4->SetLineStyle(1);npe4->SetLineColor(kRed);

 canvas_Nch_Npe->cd(1);
 SetgPad();
 if(doLog)E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_):log10(DetectorResponseEvent_.numberOfDomWithRecoHit_) >> hist_Neutrino_Nch_Npe", 
	      basecut*neutrino_slope1, "coltz");
 else E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_):(DetectorResponseEvent_.numberOfDomWithRecoHit_) >> hist_Neutrino_Nch_Npe", 
	      basecut*neutrino_slope1, "coltz");
 nch53->Draw("same");
 npe4->Draw("same");
 if(doAtm){
 canvas_Nch_Npe->cd(2);
 SetgPad();
 if(doLog) E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_):log10(DetectorResponseEvent_.numberOfDomWithRecoHit_) >> hist_AtmMuon_Nch_Npe", basecut*atmMuon_slope2, "coltz");
 else E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_):(DetectorResponseEvent_.numberOfDomWithRecoHit_) >> hist_AtmMuon_Nch_Npe",	      basecut*atmMuon_slope2, "coltz");
 nch53->Draw("same");
 npe4->Draw("same");
 }


 //TF2 *npe4 = new TF2("npe4","4", 0,900,0,900);npe4->SetLineWidth(1);npe4->SetLineStyle(2);


 //npe4->Draw("same");

 //// 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_Nch_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.16;
  double left   = 0.12;
  double top    = 0.1;
  double bottom = 0.08;
  gPad->SetRightMargin(right);
  gPad->SetLeftMargin(left);
  gPad->SetTopMargin(top);
  gPad->SetBottomMargin(bottom);
}