#include "TMath.h"

bool PrimaryEnergy_NPE_MuTauNeutrino_AtmosMuon(){

  //options

  bool doLog    = true;//if npe is in log or not
  bool doWeight = false;//
  bool tau      = true;//this is true only with 80 strings but can be false with 80 strings
  bool use_cob  = true; //either one only can be true //cob 500
  if(!use_cob)bool use_cscd_pos  = true; //cscd 600

  int nstring = 80;

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

 TCut cscd_r600="sqrt(pow(generatorCascade_.positionX_,2)+pow(generatorCascade_.positionY_,2)+pow(generatorCascade_.positionZ_,2))<600";
 TCut cob_r500 ="sqrt(pow(firstguessTrack_.cobX_,2)+pow(firstguessTrack_.cobY_,2)+pow(firstguessTrack_.cobZ_,2))<500";

//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;
 
 if(use_cob)TCut contained = cob_r500;
 else if (use_cscd_pos) TCut contained = cscd_r600;
 else TCut contained = "";
 //TCut contained = "sqrt(generatorCascade_.positionX_*generatorCascade_.positionX_+generatorCascade_.positionY_*generatorCascade_.positionY_+generatorCascade_.positionZ_*generatorCascade_.positionZ_)<600";

 //Define atmospherlic muon distribution 
 //y-axis
 int    npe_nbin        =  100;
 if(doLog){
   double npe_min_log   =  1;
   double npe_max_log   =  8;
   npe_nbin        = (int)(npe_max_log - npe_min_log)/0.1;
 }else{
   double npe_min       =  0;
   double npe_max       =  1e8;
 }
 //x-axis
 double energy_min_log    = 4.;
 double energy_max_log    = 12.;
 int    energy_nbin       = (int)(energy_max_log - energy_min_log)/0.1;

 //z-axis
 double atmMuon_min     =  1.0e-22;
 double atmMuon_max     =  1.0e-14;
 double neutrino_min    =  1.0e-30;
 double neutrino_max    =  1.0e-21;
 double tauneutrino_min    = 1.0e-40;
 double tauneutrino_max    =  1.0e-21;

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

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

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

 TAxis  *xaxis, *yaxis;
 xaxis = hist_AtmMuon_NPE->GetXaxis();
 yaxis = hist_AtmMuon_NPE->GetYaxis();
 zaxis = hist_AtmMuon_NPE->GetZaxis();
 hist_AtmMuon_NPE->SetContour(40);
 zaxis->SetRangeUser(atmMuon_min, atmMuon_max); 

 xaxis->SetTitle("log_{10} (Primary Energy [GeV]) "); 
 if(doLog)yaxis->SetTitle("log_{10} Estimated event-sum NPE");
 else yaxis->SetTitle("Estimated event-sum NPE");


 
 xaxis = hist_Neutrino_NPE->GetXaxis();
 yaxis = hist_Neutrino_NPE->GetYaxis();
 zaxis = hist_Neutrino_NPE->GetZaxis();
 hist_Neutrino_NPE->SetContour(40);
 zaxis->SetRangeUser(neutrino_min, neutrino_max); 

 xaxis->SetTitle("log_{10} (Primary Energy [GeV]) "); 
 if(doLog)yaxis->SetTitle("log_{10} Estimated event-sum NPE");
 else yaxis->SetTitle("Estimated event-sum NPE");


 
 if(tau){
   xaxis = hist_TauNeutrino_NPE->GetXaxis();
   yaxis = hist_TauNeutrino_NPE->GetYaxis();
   zaxis = hist_TauNeutrino_NPE->GetZaxis();
   hist_TauNeutrino_NPE->SetContour(40);
   if(doLog)yaxis->SetTitle("log_{10} Estimated event-sum NPE");
   else yaxis->SetTitle("Estimated event-sum NPE");
   xaxis->SetTitle("log_{10} (Primary Energy [GeV]) "); 
   zaxis->SetRangeUser(tauneutrino_min, tauneutrino_max); 
 }
 
 xaxis = hist_E1_NPE->GetXaxis();
 yaxis = hist_E1_NPE->GetYaxis();
 zaxis = hist_E1_NPE->GetZaxis();
 hist_E1_NPE->SetContour(40);
 if(doLog)yaxis->SetTitle("log_{10} Estimated event-sum NPE");
 else yaxis->SetTitle("Estimated event-sum NPE");
 xaxis->SetTitle("log_{10} (Primary Energy [GeV]) "); 
 
 xaxis = hist_E2_NPE->GetXaxis();
 yaxis = hist_E2_NPE->GetYaxis();
 zaxis = hist_E2_NPE->GetZaxis();
 hist_E2_NPE->SetContour(40);
 if(doLog)yaxis->SetTitle("log_{10} Estimated event-sum NPE");
 else yaxis->SetTitle("Estimated event-sum NPE");
 xaxis->SetTitle("log_{10} (Primary Energy [GeV]) "); 

 if(tau){
   xaxis = hist_E1Tau_NPE->GetXaxis();
   yaxis = hist_E1Tau_NPE->GetYaxis();
   zaxis = hist_E1Tau_NPE->GetZaxis();
   hist_E1Tau_NPE->SetContour(40);
   if(doLog)yaxis->SetTitle("log_{10} Estimated event-sum NPE");
   else yaxis->SetTitle("Estimated event-sum NPE");
   xaxis->SetTitle("log_{10} (Primary Energy [GeV]) "); 
 }
 
 //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_):(generatorPrimary_.energyLog_) >> hist_Neutrino_NPE",  
			 (basecut&&contained)*neutrino_slope1, "coltz");
   else E1Tree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_Neutrino_NPE", 
		     (basecut&&contained)*neutrino_slope1, "coltz");
   highnpe->Draw("same");
   horizon->Draw("same");

 } else if(!doWeight) {


   if(doLog)E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_E1_NPE",
			 basecut&&(!contained), "");
   else E1Tree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_E1_NPE", 
		     basecut&&(!contained), "");

   
   if(doLog)E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_E1_NPE_c", 
			 basecut&&contained, "");
   else E1Tree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_E1_NPE_c", 
		     basecut&&contained, "");

   if(!doWeight){
     hist_E1_NPE->SetMarkerColor(9);
     hist_E1_NPE->SetMarkerStyle(2);
     hist_E1_NPE->DrawCopy();
     hist_E1_NPE_c->SetMarkerColor(8);
     hist_E1_NPE_c->SetMarkerStyle(2);
     hist_E1_NPE_c->DrawCopy("same");
     //highnpe->Draw("same");
     //horizon->Draw("same");
   }
 }
 
 
 canvas_NPE->cd(2);
 SetgPad();
 if(doWeight){
   if(doLog) E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_AtmMuon_NPE", 
			  (basecut&&contained)*atmMuon_slope2, "coltz");
   else E2Tree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_AtmMuon_NPE", 
		     (basecut&&contained)*atmMuon_slope2, "coltz");
   highnpe->Draw("same");
   horizon->Draw("same");
 }else{
   if(doLog) E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_E2_NPE",
			  basecut&&(!contained), "");

   else E2Tree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_E2_NPE",
		     basecut&&(!contained), "");

   if(doLog) E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_E2_NPE_c",
			  basecut&&contained, "");

   else E2Tree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_E2_NPE_c",
		     basecut&&contained, "");

    if(!doWeight){

      hist_E2_NPE->SetMarkerColor(9);
      hist_E2_NPE->SetMarkerStyle(2);
      hist_E2_NPE_c->SetMarkerColor(8);
      hist_E2_NPE_c->SetMarkerStyle(2);
      hist_E2_NPE->DrawCopy();
      hist_E2_NPE_c->DrawCopy("same");
      //highnpe->Draw("same");
      //horizon->Draw("same");
    }
 }


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

     if(doLog)E1TauTree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_E1Tau_NPE_c", 
			      basecut&&contained, "");
     else E1TauTree->Draw("(DetectorResponseEvent_.totalBestEstimatedNPE_):(generatorPrimary_.energyLog_) >> hist_E1Tau_NPE_c", 
			  basecut&&contained, "");
     
     if(!doWeight){

       hist_E1Tau_NPE->SetMarkerColor(10);
       hist_E1Tau_NPE->SetMarkerStyle(2);
       //hist_E1Tau_NPE->Clear();
       hist_E1Tau_NPE->Draw("");
       hist_E1Tau_NPE_c->SetMarkerColor(8);
       hist_E1Tau_NPE_c->SetMarkerStyle(2);
       hist_E1Tau_NPE_c->DrawCopy("samecoltz");
     }
   }
   
 }
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////



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