#include <TTree.h>
#include <TFile.h>
#include <TH2D.h>


double calculate_shell_edep(TH2D* edepH,double r,int cone_nb)
{
   double edep=0;
   double mass=0;
   for(int j=0;j<=r;j++)
      for(int i=0;i<cone_nb;i++)
      {
      edep+=edepH->GetBinContent(i,j);
      // mass+=massH->GetBinContent(i,r);

      }
return edep;
}
TTree* Calculate_Edep_Fractions(TH2D* h,int cone_nb)
{

double _r;
double _f;


TTree* _result = new TTree("Ntuple1", "Edep");
_result->Branch("r", &_r, "r/D");
_result->Branch("f", &_f, "f/D");



for(int i=0;i<=25;i++)
{
  _r=(0.5*(2*i*i+1)/100);
  _f=calculate_shell_edep(h,i,cone_nb)/(1.25*1e6-3186);
//   cout
//   <<" ----->> r : "<<_r
//   <<" ----->> f : "<<_f
//   <<endl;
  _result->Fill();
}

return _result;
}
TTree* Calculate_Total_edep_frachtion(TH2D* h1,TH2D* h2,TH2D* h3,TH2D*h4,TH2D*h5)
{
   double _r;
   double _f;


   TTree* _result = new TTree("Ntuple1", "Edep");
   _result->Branch("r", &_r, "r/D");
   _result->Branch("f", &_f, "f/D");
      for(int i=0;i<25;i++)
      {
         _r=(0.5*(2*i*i+1)/100);
         _f=(calculate_shell_edep(h1,i,48)
         +calculate_shell_edep(h2,i,48)
         +calculate_shell_edep(h3,i,48)
         +calculate_shell_edep(h4,i,48)
         +calculate_shell_edep(h5,i,48))/(1.25*1e6-3186);
      // cout
      //    << "     r :"<< _r
      //    << "     f :"<< _f
      //    <<endl;
         _result->Fill();
      }
_result->SetLineColor(kBlack);
_result->SetLineWidth(2); 
_result->Draw("f:r","","L C");
TH1F *htemp1 = (TH1F*)gPad->GetPrimitive("htemp");
htemp1->SetTitle("The fraction of incident energy deposited ; r (g/cm2) ; fraction of energy deposited within shell");
      return _result;

}
void Plot_Kernels()
{
   TFile* f =new TFile("result.root");
TCanvas* c1 = new TCanvas("c1", "",8000,10000);
TH2D* edepPrimaryH = (TH2D*)f->Get("2d_PrimaryedepKernels");
TH2D* edepFirstH = (TH2D*)f->Get("2D_FirstedepKernels");
TH2D* edepSecondH = (TH2D*)f->Get("2D_SecondedepKernels");
TH2D* edepMultipleH = (TH2D*)f->Get("2D_MultipleedepKernels");
TH2D* edepeBremAnnihilH = (TH2D*)f->Get("2D_eBremAnnihiledepKernels");
TH2D* massH = (TH2D*)f->Get("2d_Mass");

TTree* _Primary=Calculate_Edep_Fractions(edepPrimaryH,48);
TTree* _First=Calculate_Edep_Fractions(edepFirstH,48);
TTree* _Multiple=Calculate_Edep_Fractions(edepMultipleH,48);
TTree* _Second=Calculate_Edep_Fractions(edepSecondH,48);
TTree* _eBremAnnihil=Calculate_Edep_Fractions(edepeBremAnnihilH,48);
TTree* _total = Calculate_Total_edep_frachtion(edepPrimaryH,edepFirstH,edepSecondH,edepMultipleH,edepeBremAnnihilH);


// _total->SetLineColor(kBlack);
// _total->SetLineWidth(2); 
// _total->Draw("f:r","","L same");

_Primary->SetLineColor(kBlack);
_Primary->SetLineWidth(2); 
_Primary->Draw("f:r","","L same c");

_First->SetLineColor(kBlack);
_First->SetLineWidth(2); 
_First->Draw("f:r","","L same C");

_Second->SetLineColor(kBlack);
_Second->SetLineWidth(2); 
_Second->Draw("f:r","","L same C");

_eBremAnnihil->SetLineColor(kBlack);
_eBremAnnihil->SetLineWidth(2); 
_eBremAnnihil->Draw("f:r","","L same C"); 

_Multiple->SetLineColor(kBlack);
_Multiple->SetLineWidth(2); 
_Multiple->Draw("f:r","","L same C");

   gPad->SetLogy(1);
   gPad->SetLogx(1);
   c1->cd();
   c1->Update();
}


void plot() 
{ 
Plot_Kernels();
}