#include "TMath.h"
bool NPE_CosTheta_Tau(){

  bool dowrite=false;
  gSystem->Load("libweighting-module");
  /*style*/
  SetStyle();
  gROOT->ForceStyle();

  string indir       = "/home/aya/IceCube/ehe/March06/RootReducedLeaves/";
  string file_slope1 = "JulietTauGZKMergedAll.reduced.root";
  
// 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");

 cout <<"********************************************************************"<< endl;
 cout <<"Energy E-1 file Tau " << filenameE1 << " is " << E1Tree->GetEntries()<< 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 consts 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;

 //y-axis 
 int    npe_nbin      =  100;
 double npe_min       =  0.0e0;
 double npe_max       =  1.0e5;
 //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-14;

 double neutrino_min   = 1.0e-25;
 double neutrino_max   = 1.0e-21;

 TH2D* hist_TauNeutrino_NPE_CosTheta = new TH2D("hist_TauNeutrino_NPE_CosTheta","GZK Tau Neutrino from E^{-1} sample", 
					    npe_nbin, npe_min, npe_max, costheta_nbin, costheta_min, costheta_max);
 TAxis  *xaxis, *yaxis, *zaxis;
 xaxis = hist_TauNeutrino_NPE_CosTheta->GetXaxis();
 yaxis = hist_TauNeutrino_NPE_CosTheta->GetYaxis();
 zaxis = hist_TauNeutrino_NPE_CosTheta->GetZaxis();
 hist_TauNeutrino_NPE_CosTheta->SetContour(40);
 xaxis->SetTitle("Estimated NPE");
 yaxis->SetTitle("MC Truth cos #theta"); 
 zaxis->SetRangeUser(neutrino_min, neutrino_max); 

 //Start Drawing
 TCanvas *canvas_NPE_CosTheta = new TCanvas("NPE_CosTheta","Npe : cos(theta)", 500, 500);
 canvas_NPE_CosTheta->Divide(1, 0.001);  
 canvas_NPE_CosTheta->SetHighLightColor(10);

 canvas_NPE_CosTheta->cd(1);
 SetgPad();
 E1Tree->Draw("generatorPrimary_.cosTheta_:DetectorResponseEvent_.totalBestEstimatedNPE_ >> hist_TauNeutrino_NPE_CosTheta", 
	      basecut*neutrino_slope1, "coltz");// cut*weight
 canvas_NPE_CosTheta->Update();

 if(dowrite)
   {
     string outfile = "outhist.root";
     TFile *hist_outfile = new TFile(outfile.c_str(),"recreate");
     hist_TauNeutrino_NPE_CosTheta->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(0.02);
   gStyle->SetLabelSize(0.03,"x");
   gStyle->SetLabelSize(0.03,"y");
   gStyle->SetLabelSize(0.03,"z");
   gStyle->SetLabelOffset(0.002,"x");
   gStyle->SetNdivisions(1010,"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);
}
//////////////////////////////////