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

//######################################################################################################
double calculate_shell_edep(TH2D* edepH,double r,int cone_nb)
{
   double edep=0;
   for(int j=0;j<=r;j++)
      for(int i=0;i<cone_nb;i++)
      {
      edep+=edepH->GetBinContent(i,j);
      }
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)/10);
  _f=calculate_shell_edep(h,i,cone_nb)/(5*(1e6-361));
//   cout
//   <<" ----->> r : "<<_r
//   <<" ----->> f : "<<_f
//   <<endl;
  _result->Fill();
}

_result->ResetBranchAddresses(); // disconnect from local variables
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)/10);
         _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))/(5*(1e6-361));
      // cout
      //    << "     r :"<< _r
      //    << "     f :"<< _f
      //    <<endl;
         _result->Fill();
      }

   _result->ResetBranchAddresses(); // disconnect from local variables
   return _result;
}

void Plot_Kernels()
{
   TFile* f =new TFile("result.root");
   TCanvas* c1 = new TCanvas("c1", "",1000, 800);
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);

TMultiGraph *mg = new TMultiGraph("mg", "The fraction of incident energy deposited;sphere radius #[]{cm};#voidb_{#Epsilon}");

_total->SetLineColor(kBlack);
_total->SetLineWidth(2);
_total->Draw("f:r", "");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

_Primary->SetLineColor(kBlack);
_Primary->SetLineWidth(2);
_Primary->Draw("f:r", "");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

_First->SetLineColor(kBlack);
_First->SetLineWidth(2);
_First->Draw("f:r", "");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

_Second->SetLineColor(kBlack);
_Second->SetLineWidth(2);
_Second->Draw("f:r", "");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

_eBremAnnihil->SetLineColor(kBlack);
_eBremAnnihil->SetLineWidth(2);
_eBremAnnihil->Draw("f:r", ""); 
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

_Multiple->SetLineColor(kBlack);
_Multiple->SetLineWidth(2); 
_Multiple->Draw("f:r", "");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

mg->Draw("ALC");
gPad->SetLogy(1);
gPad->SetLogx(1);

   auto *th1 = new TText(5,0.8,"Total");
   th1->SetTextAlign(6); th1->SetTextSize(0.02);
   th1->Draw();

   auto *th2 = new TText(5,0.5,"Primary");
   th2->SetTextAlign(5); th2->SetTextSize(0.02);
   th2->Draw();

   auto *th3 = new TText(5,0.4,"First Scatter");
   th3->SetTextAlign(4); th3->SetTextSize(0.02);
   th3->Draw();

   auto *th4 = new TText(5,0.3,"Second Scatter");
   th4->SetTextAlign(3); th4->SetTextSize(0.02);
   th4->Draw();

   auto *th5 = new TText(5,0.2,"Multiple Scatter");
   th5->SetTextAlign(2); th5->SetTextSize(0.02);
   th5->Draw();

   auto *th6 = new TText(5,0.1,"Bremsstrahlung and Annihilation Scatter");
   th6->SetTextAlign(1); th6->SetTextSize(0.02);
   th6->Draw();
   
gPad->Modified(); gPad->Update();
}
//######################################################################################################



//######################################################################################################
double calculate_voxel_dose_R(TH2D* kernels[5],TH2D* massH,double r,int cone_nb_start,int cone_nb_end)
{
   double edep=0;
   double mass=0;
   double no_photon=1e6-3170;
   // double to_gray=1.609*1e-13;
   // double incident_e=1.25;
   double scale=1/no_photon;

   for(int kid=0;kid<5;kid++)
      for(int i=cone_nb_start;i<cone_nb_end;i++)
         {
         edep+=kernels[kid]->GetBinContent(i,r);
         }
   for(int i=cone_nb_start;i<cone_nb_end;i++)
      {
      mass+=massH->GetBinContent(i,r);
      }
return (edep/mass)*scale;
}

void Plot_r2Dose_R() 
{
   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* Reff = (TH2D*)f->Get("2D_Reffective");
   TH2D* massH = (TH2D*)f->Get("2d_Mass");


// for(int r=0;r<25;r++)
// for(int phi=0;phi<48;phi++)
// cout
// <<" ----->> R = "<<r
// <<" ----->> phi = "<<phi
// <<" ----->> Reff = "<<Reff->GetBinContent(phi,r)
// <<endl;



   double _r;
   double _doser2_angle1;
   double _doser2_angle2;
   double _doser2_angle3;


   TTree* _angle1 = new TTree("Ntuple1", "Edep");
   _angle1->Branch("r", &_r, "r/D");
   _angle1->Branch("f", &_doser2_angle1, "f/D");
   //
      TTree* _angle2 = new TTree("Ntuple1", "Edep");
   _angle2->Branch("r", &_r, "r/D");
   _angle2->Branch("f", &_doser2_angle2, "f/D");

   //
   TTree* _angle3 = new TTree("Ntuple1", "Edep");
   _angle3->Branch("r", &_r, "r/D");
   _angle3->Branch("f", &_doser2_angle3, "f/D");
   //
   //    TTree* _angle4 = new TTree("Ntuple1", "Edep");
   // _angle4->Branch("r", &_r, "r/D");
   // _angle4->Branch("f", &_doser2, "f/D");
   // //
   // TTree* _angle5 = new TTree("Ntuple1", "Edep");
   // _angle5->Branch("r", &_r, "r/D");
   // _angle5->Branch("f", &_doser2, "f/D");
   TH2D* kernels[5];
   kernels[0]=edepPrimaryH;
   kernels[1]=edepFirstH;
   kernels[2]=edepSecondH;
   kernels[3]=edepMultipleH;
   kernels[4]=edepeBremAnnihilH;

   for(int i=0;i<15;i++)
   {
   _r=(0.5*(2*i*i+1))/10;//*(0.5*(2*i*i+1));
   double r1=(0.5*(2*i*i+1))/10;
   double r2=(0.5*(2*(i+1)*(i+1)+1))/10;
   double r=(r2-r1)/2;



   _doser2_angle1=calculate_voxel_dose_R(kernels,massH,i,0,1)*r*r;
   _doser2_angle2=calculate_voxel_dose_R(kernels,massH,i,6,7)*r*r;
   _doser2_angle3=calculate_voxel_dose_R(kernels,massH,i,12,13)*r*r;
   cout
   <<" ----->> R = "<<_r
   <<" ----->> dose = "<<_doser2_angle1
   <<endl;

   //   cout
   //   <<" ----->> r : "<<_r
   //   <<" ----->> f : "<<_doser2
   //   <<endl;
   _angle1->Fill();
   _angle2->Fill();
   _angle3->Fill();

   }
TMultiGraph *mg = new TMultiGraph("mg", ";r #[]{cm};Dose #times r^{2}_{} #[]{#frac{MeV.cm^{2}_{}}{Kg.Photon}} ");

   _angle1->SetLineColor(kBlack);
   _angle1->SetLineWidth(1); 
   _angle1->Draw("f:r","","");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());



   _angle2->SetLineColor(kBlack);
   _angle2->SetLineWidth(2); 
   _angle2->Draw("f:r","","");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());




   _angle3->SetLineColor(kBlack);
   _angle3->SetLineWidth(2); 
   _angle3->Draw("f:r","","");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

   
mg->Draw("ALC");


 
   TLatex latex;
   latex.SetTextSize(0.025);
   latex.SetTextAlign(13);  //align at top
   latex.DrawLatex(.2,.9,"0^{#circ} - 3.75^{#circ}");
   latex.DrawLatex(.3,.9,"22.5^{#circ} - 26.25^{#circ}");
   latex.DrawLatex(.3,.9,"45^{#circ} - 48.75^{#circ}");

   // auto *th1 = new TText(0.80,0.185,"0 - 3.75");
   // th1->SetTextAlign(12); th1->SetTextSize(0.025);
   // th1->Draw();

   // auto *th2 = new TText(0.70,0.165,"22.5 - 26.25");
   // th2->SetTextAlign(11); th2->SetTextSize(0.025);
   // th2->Draw();

   // auto *th3 = new TText(0.60,0.195,"45 - 48.75");
   // th3->SetTextAlign(13); th3->SetTextSize(0.025);
   // th3->Draw();

gPad->SetLogy(1);
// gPad->SetLogx(1);
gPad->Modified(); gPad->Update();
}
//######################################################################################################



//######################################################################################################
double calculate_voxel_dose_Theta(TH2D* kernels,TH2D* massH,double r,int cone_nb)
{
    double edep=0;
   double mass=0;
   double no_photon=1e6-3170;
   // double to_gray=1.609*1e-13;
   // double incident_e=1.25;
   double scale=1/no_photon;
      edep+=kernels->GetBinContent(cone_nb,r);
      mass+=massH->GetBinContent(cone_nb,r);
      cout
       <<"--------- >> R =  "<<r
        <<"--------- >> phi =  "<<3.75*cone_nb
      <<"--------- >> Edep =  "<<edep<<endl;
return (edep/mass)*scale;
}

void Plot_r2Dose_Theta()
{
     TFile* f =new TFile("result.root");
   TCanvas* c1 = new TCanvas("c1", "",8000,10000);
   TH2D* edepPrimaryH = (TH2D*)f->Get("2d_PrimaryedepKernels");
   TH2D* massH = (TH2D*)f->Get("2d_Mass");


// for(int r=0;r<25;r++)
// for(int phi=0;phi<48;phi++)
// cout
// <<" ----->> R = "<<r
// <<" ----->> phi = "<<phi
// <<" ----->> Reff = "<<Reff->GetBinContent(phi,r)
// <<endl;



   double _r;
   double _doser2_angle1;
   double _doser2_angle2;
   double _doser2_angle3;


   TTree* _angle1 = new TTree("Ntuple1", "Edep");
   _angle1->Branch("r", &_r, "r/D");
   _angle1->Branch("f", &_doser2_angle1, "f/D");
   //
      TTree* _angle2 = new TTree("Ntuple1", "Edep");
   _angle2->Branch("r", &_r, "r/D");
   _angle2->Branch("f", &_doser2_angle2, "f/D");

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

   for(int i=0;i<48;i++)
   {
   _r=TMath::Cos((i*3.75*TMath::Pi())/180);//*(0.5*(2*i*i+1));
   //r0
    double r0=(0.5*(2*0*0+1))/10;
    //r1
    double r1=(0.5*(2*1*1+1))/10;
    //r2
    double r2=(0.5*(2*2*2+1))/10;
   _doser2_angle1=  calculate_voxel_dose_Theta(edepPrimaryH,massH,0,i)*r0*r0;
   _doser2_angle2=  calculate_voxel_dose_Theta(edepPrimaryH,massH,1,i)*r1*r1;
   _doser2_angle3=  calculate_voxel_dose_Theta(edepPrimaryH,massH,2,i)*r2*r2;

   _angle1->Fill();
   _angle2->Fill();
   _angle3->Fill();

   }
   // TMultiGraph *mg = new TMultiGraph("mg", ";r (g/cm^{2}_{});Dose x r^{2}_{}");
   TMultiGraph *mg = new TMultiGraph("mg", ";Cos#(){#theta};Dose #times r^{2}_{} #[]{#frac{MeV.cm^{2}_{}}{Kg.Photon}} ");
   // TLatex latex;
   // latex.SetTextSize(0.025);
   // latex.SetTextAlign(13);  //align at top
   // latex.DrawLatex(.2,.9,"r = 2.55 cm");
   // latex.DrawLatex(.3,.9,"r = 0.05 cm");
   // latex.DrawLatex(.3,.9,"r = 10.05 cm");
   _angle1->SetLineColor(kBlack);
   _angle1->SetLineWidth(2); 
   _angle1->Draw("f:r","","L C");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());



   _angle2->SetLineColor(kBlack);
   _angle2->SetLineWidth(2); 
   _angle2->Draw("f:r","","L C");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());




   _angle3->SetLineColor(kBlack);
   _angle3->SetLineWidth(2); 
   _angle3->Draw("f:r","","L C");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

   
mg->Draw("ALC");


 
   // TLatex latex;
   // latex.SetTextSize(0.025);
   // latex.SetTextAlign(13);  //align at top
   // latex.DrawLatex(.2,.9,"0^{#circ} - 3.75^{#circ}");
   // latex.DrawLatex(.3,.9,"22.5^{#circ} - 26.25^{#circ}");
   // latex.DrawLatex(.3,.9,"45^{#circ} - 48.75^{#circ}");

   // auto *th1 = new TText(0.80,0.185,"0 - 3.75");
   // th1->SetTextAlign(12); th1->SetTextSize(0.025);
   // th1->Draw();

   // auto *th2 = new TText(0.70,0.165,"22.5 - 26.25");
   // th2->SetTextAlign(11); th2->SetTextSize(0.025);
   // th2->Draw();

   // auto *th3 = new TText(0.60,0.195,"45 - 48.75");
   // th3->SetTextAlign(13); th3->SetTextSize(0.025);
   // th3->Draw();

gPad->SetLogy(1);
// gPad->SetLogx(1);
gPad->Modified(); gPad->Update();
}


//######################################################################################################
void Plot_DoseR2_CosTheta()
{


   TFile* f =                new TFile("result.root");
TCanvas* c1 =             new TCanvas("c1", "");
// c1->Divide(2,1);

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");
// 0.5*(2*ir*ir+1)

   // const Int_t NBINS = 26;
   // Double_t edges[NBINS + 1] = {
   //    0.0,
   //    // 1e-6,
   //    0.5,
   //    1.5,
   //    4.5,
   //    9.5,
   //    16.5,
   //    25.5,
   //    36.5,
   //    49.5,
   //    64.5,
   //    81.5,
   //    100.5,
   //    121.5,
   //    144.5,
   //    169.5,
   //    196.5,
   //    225.5,
   //    256.5,
   //    289.5,
   //    324.5,
   //    361.5,
   //    400.5,
   //    441.5,
   //    484.5,
   //    529.5,
   //    576.5,
   //    625.5};


   double ttreedose;
   double ttreepol;
   double ttreer;
   TTree* _angle1 = new TTree("Ntuple1", "Edep");
   _angle1->Branch("d", &ttreedose, "d/D");
   _angle1->Branch("p", &ttreepol, "p/D");
   _angle1->Branch("r", &ttreer, "r/D");

   double no_photon=1e6-3170;
   double to_gray=1.609*1e-13;
   double incident_e=1.25;
   double to_centigram=1e2;
   // double scale=(to_gray*to_centigram)/(incident_e*no_photon);
      double scale=1/(no_photon);

// TH2D* isoprimary = new TH2D("Primary Charged Particles","Primary Charged Particles",48,0,48,NBINS,edges);
// TH2D* isofirst = new TH2D("First Scatter Charged Particles","First Scatter Charged Particles",96,0,96,25,0,25);

for(Int_t iphi=0;iphi<48;iphi++)
   for(Int_t ir=0;ir<25;ir++)
   {
      // iso->SetBinContent(iphi,ir,
      // edepFirstH->GetBinContent(iphi,ir)/massH->GetBinContent(iphi,ir)

            double edep=edepPrimaryH->GetBinContent(iphi,ir)
            +edepFirstH->GetBinContent(iphi,ir)
            +edepSecondH->GetBinContent(iphi,ir)
            +edepMultipleH->GetBinContent(iphi,ir)
            +edepeBremAnnihilH->GetBinContent(iphi,ir);
            // double dose = edepPrimaryH->GetBinContent(iphi,ir)/massH->GetBinContent(iphi,ir);
            double dose = edep/massH->GetBinContent(iphi,ir);
            ttreer=0.5*(2*ir*ir+1)/10;
            ttreepol=iphi;
            ttreedose=dose*scale;
            _angle1->Fill();
   }


   _angle1->SetLineColor(kBlack);
   _angle1->SetLineWidth(2); 
TMultiGraph *mg = new TMultiGraph("mg", ";Cos#(){#theta};Dose#timesr^{2}_{} #[]{#frac{MeV.cm^{2}_{}}{Kg.Photon}}");



   _angle1->Draw("d*0.05*0.05:cos((p*3.75*3.14159265358979/180))","r==0.05","");
   mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

   _angle1->Draw("d*0.45*0.45:cos((p*3.75*3.14159265358979/180))","r==0.45","");
   mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

   _angle1->Draw("d*0.15*0.15:cos((p*3.75*3.14159265358979/180))","r==0.15","");
   mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

   _angle1->Draw("d*0.95*0.95:cos((p*3.75*3.14159265358979/180))","r==0.95","");
   mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());
   mg->Draw("ALC");
   gPad->SetLogy(1);




   //  TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
   //  htemp->GetXaxis()->SetTitle("r#[]{cm}");
   //  htemp->GetYaxis()->SetTitle("Dose#timesr^{2}_{} #[]{#frac{MeV.cm^{2}_{}}{Kg.Photon}} "); 





 TLatex latex;
   latex.SetTextSize(0.024);
   latex.SetTextAlign(13);  //align at top
   latex.DrawLatex(.5,5,"r = 0.05 cm");
   latex.DrawLatex(.5,5,"r = 0.45 cm");
   latex.DrawLatex(.5,5,"r = 0.15 cm");
   latex.DrawLatex(.5,5,"r = 0.95 cm");



   // gPad->Update();
   // auto gp = new TGraphPolargram("g",0,0.5, 0., 180.);
   // gp->SetToDegree ();
   // gp->SetNdivPolar(4);
   // gp->SetNdivRadial(4);
   // gp->SetLineColor(1);
   // gp->Draw();
gPad->Modified(); gPad->Update();

// c1->cd();
// c1->Update();
}
//######################################################################################################


//######################################################################################################
void Plot_DoseR2_R(){
   TFile* f =                new TFile("result.root");
TCanvas* c1 =             new TCanvas("c1", "");
// c1->Divide(2,1);

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");
// 0.5*(2*ir*ir+1)

   const Int_t NBINS = 26;
   Double_t edges[NBINS + 1] = {
      0.0,
      // 1e-6,
      0.5,
      1.5,
      4.5,
      9.5,
      16.5,
      25.5,
      36.5,
      49.5,
      64.5,
      81.5,
      100.5,
      121.5,
      144.5,
      169.5,
      196.5,
      225.5,
      256.5,
      289.5,
      324.5,
      361.5,
      400.5,
      441.5,
      484.5,
      529.5,
      576.5,
      625.5};


   double ttreedose;
   double ttreepol;
   double ttreer;
   TTree* _angle1 = new TTree("Ntuple1", "Edep");
   _angle1->Branch("d", &ttreedose, "d/D");
   _angle1->Branch("p", &ttreepol, "p/D");
   _angle1->Branch("r", &ttreer, "r/D");

   double no_photon=1e6-3170;
   double to_gray=1.609*1e-13;
   double incident_e=1.25;
   double to_centigram=1e2;
   // double scale=(to_gray*to_centigram)/(incident_e*no_photon);
      double scale=1/(no_photon);

TH2D* isoprimary = new TH2D("Primary Charged Particles","Primary Charged Particles",48,0,48,NBINS,edges);
TH2D* isofirst = new TH2D("First Scatter Charged Particles","First Scatter Charged Particles",96,0,96,25,0,25);

for(Int_t iphi=0;iphi<48;iphi++)
   for(Int_t ir=0;ir<25;ir++)
   {
      // iso->SetBinContent(iphi,ir,
      // edepFirstH->GetBinContent(iphi,ir)/massH->GetBinContent(iphi,ir)

            double edep=edepPrimaryH->GetBinContent(iphi,ir)
            +edepFirstH->GetBinContent(iphi,ir)
            +edepSecondH->GetBinContent(iphi,ir)
            +edepMultipleH->GetBinContent(iphi,ir)
            +edepeBremAnnihilH->GetBinContent(iphi,ir);
            // double dose = edepPrimaryH->GetBinContent(iphi,ir)/massH->GetBinContent(iphi,ir);
            double dose = edep/massH->GetBinContent(iphi,ir);
            isoprimary->SetBinContent(iphi,0.5*(2*ir*ir+1),dose);
            ttreer=0.5*(2*ir*ir+1)/10;
            ttreepol=iphi;
            ttreedose=dose*scale;
            _angle1->Fill();
   }



   _angle1->SetLineColor(kBlack);
   _angle1->SetLineWidth(2); 



TMultiGraph *mg = new TMultiGraph("mg", ";r#[]{cm};Dose#timesr^{2}_{} #[]{#frac{MeV.cm^{2}_{}}{Kg.Photon}}");






 _angle1->Draw("d*r*r:r","p==0 && r<20","");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

 _angle1->Draw("d*r*r:r","p==6&&r<20","");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

 _angle1->Draw("d*r*r:r","p==12&&r<20","");
mg->Add((TGraph*)gPad->GetPrimitive("Graph")->Clone());

mg->Draw("AL");
gPad->SetLogy(1);

   // TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
   //  htemp->GetXaxis()->SetTitle("r#[]{cm}");
   //  htemp->GetYaxis()->SetTitle("Dose#timesr^{2}_{} #[]{#frac{MeV.cm^{2}_{}}{Kg.Photon}} "); 





 TLatex latex;
   latex.SetTextSize(0.018);
   latex.SetTextAlign(13);  //align at top
   latex.DrawLatex(5,5,"0^{#circ} - 3.75^{#circ}");
   latex.DrawLatex(5,5,"22.5^{#circ} - 26.25^{#circ}");
   latex.DrawLatex(5,5,"45^{#circ} - 48.75^{#circ}");
   // latex.DrawLatex(5,5,"r = 0.15 cm");
   // latex.DrawLatex(.3,.9,"Cos#(){#theta}");
   // latex.DrawLatex(.3,.9,"Dose#timesr^2#[]{}");

   // gPad->Update();
   // auto gp = new TGraphPolargram("g",0,0.5, 0., 180.);
   // gp->SetToDegree ();
   // gp->SetNdivPolar(4);
   // gp->SetNdivRadial(4);
   // gp->SetLineColor(1);
   // gp->Draw();

c1->cd();
c1->Update();
}
//######################################################################################################






//######################################################################################################
double Calculate_BackScatter_Edep_Primary(double E0)
{
   TFile* f =new TFile("result.root");
   TH2D* edepPrimaryH = (TH2D*)f->Get("2d_PrimaryedepKernels");
   double edep=0;
      for(int phi=23;phi<48;phi++)
         for(int r=0;r<25;r++)
         edep+=
         edepPrimaryH->GetBinContent(phi,r);
   return (edep/E0)*100;
}
//######################################################################################################
double Calculate_BackScatter_Edep_Total(double E0)
{
   TFile* f =new TFile("result.root");
   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");
   double edep=0;
      for(int phi=23;phi<48;phi++)
         for(int r=0;r<25;r++)
         edep+=
         edepPrimaryH->GetBinContent(phi,r)
         +edepFirstH->GetBinContent(phi,r)
         +edepSecondH->GetBinContent(phi,r)
         +edepMultipleH->GetBinContent(phi,r)
         +edepeBremAnnihilH->GetBinContent(phi,r);
   return (edep/E0)*100;
}
//######################################################################################################



//######################################################################################################
double Calculate_Ftot(double E0)
{
   double fp;
      TFile* f =new TFile("result.root");
   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");
   double edep=0;
      for(int phi=0;phi<48;phi++)
         for(int r=0;r<25;r++)
         fp+=
         edepPrimaryH->GetBinContent(phi,r)
         +edepFirstH->GetBinContent(phi,r)
         +edepSecondH->GetBinContent(phi,r)
         +edepMultipleH->GetBinContent(phi,r)
         +edepeBremAnnihilH->GetBinContent(phi,r);
   return (fp/E0)*100;
}
//######################################################################################################
double Calculate_Fp(double E0)
{
   double fp;
   TFile* f =new TFile("result.root");
   TH2D* edepPrimaryH = (TH2D*)f->Get("2d_PrimaryedepKernels");
      for(int phi=0;phi<48;phi++)
         for(int r=0;r<25;r++)
         fp+=
         edepPrimaryH->GetBinContent(phi,r);
   return (fp/E0)*100;
}
//######################################################################################################

//######################################################################################################
double Calculate_Fp_without_eioni(double E0)
{
   double fp;
   TFile* f =new TFile("result.root");
   TH2D* edepPrimaryH = (TH2D*)f->Get("2d_PrimaryedepKernels_without_eioni_sec");
      for(int phi=0;phi<48;phi++)
         for(int r=0;r<25;r++)
         fp+=
         edepPrimaryH->GetBinContent(phi,r);
   return (fp/E0)*100;
}
//######################################################################################################


//######################################################################################################
double Calculate_Reff(double E0)
{
   double fp;
   TFile* f =new TFile("result.root");
   TH2D* edepPrimaryH = (TH2D*)f->Get("2d_PrimaryedepKernels_without_eioni_sec");
   TH2D* reff = (TH2D*)f->Get("2D_Reffective");
   double Reff=0;
   double Fp=Calculate_Fp_without_eioni(E0);
      for(int phi=0;phi<48;phi++)
         for(int r=0;r<25;r++)
            {
               double reffij=reff->GetBinContent(phi,r);
               double edepij=edepPrimaryH->GetBinContent(phi,r);
               Reff+=((reffij*edepij)/E0)/Fp;
            }         
   return Reff;
}
//######################################################################################################


//######################################################################################################
double Calculate_Zeff(double E0)
{
   double fp;
   TFile* f =new TFile("result.root");
   TH2D* edepPrimaryH = (TH2D*)f->Get("2d_PrimaryedepKernels_without_eioni_sec");
   TH2D* zeff = (TH2D*)f->Get("2D_Zeffective");
   double Zeff=0;
   double Fp=Calculate_Fp_without_eioni(E0);
      for(int phi=0;phi<48;phi++)
         for(int r=0;r<25;r++)
            {
               double zeffij=zeff->GetBinContent(phi,r);
               double edepij=edepPrimaryH->GetBinContent(phi,r);
               Zeff+=((zeffij*edepij)/E0)/Fp;
            }         
   return Zeff;
}
//######################################################################################################


//######################################################################################################
double Calculate_Yeff(double E0)
{
   double reff=Calculate_Reff(E0);
   double zeff=Calculate_Zeff(E0);
   double Yeff=sqrt(reff*reff + zeff*zeff);
   return Yeff;
}
//######################################################################################################

//######################################################################################################
double Calculate_B11(double E0,double attenuation_coefficient,double density)
{
   double fp;
   TFile* f =new TFile("result.root");
   TH2D* edepPrimaryH = (TH2D*)f->Get("2d_PrimaryedepKernels_without_eioni_sec");
   TH2D* zeff = (TH2D*)f->Get("2D_Zeffective");
   double Zeff=0;
   double Fp=Calculate_Fp_without_eioni(E0);
      for(int phi=0;phi<48;phi++)
         for(int r=0;r<25;r++)
            {
               double zeffij=zeff->GetBinContent(phi,r);
               double edepij=edepPrimaryH->GetBinContent(phi,r);
               Zeff+=((TMath::Exp(zeffij/10*(attenuation_coefficient/density))*edepij)/E0)/Fp;
            }         
   return Zeff;
}
//######################################################################################################

//######################################################################################################
double Calculate_B12(double E0,double attenuation_coefficient,double density) // cm2/g , g/cm3
{

   return 1+(attenuation_coefficient/density)*Calculate_Zeff(E0);
}
//######################################################################################################

//######################################################################################################
double Calculate_B13(double E0,double attenuation_coefficient,double density)
{
   return TMath::Exp((attenuation_coefficient/density)*Calculate_Zeff(E0));
}
//######################################################################################################



void IsoLevelPlot()
{
   TFile* f =new TFile("result.root");
   TCanvas* c1 = new TCanvas("c1", "",8000,10000);
   c1->SetTheta(90);
   c1->SetPhi(0.);

   // double ttreedose;
   // double ttreepol;
   // double ttreer;
   // TTree* _angle1 = new TTree("Ntuple1", "Edep");
   // _angle1->Branch("d", &ttreedose, "d/D");
   // _angle1->Branch("p", &ttreepol, "p/D");
   // _angle1->Branch("r", &ttreer, "r/D");


   TH2D* edepPrimaryH = (TH2D*)f->Get("2d_PrimaryedepKernels");
//    for(int phi=0;phi<96;phi++)
//       for(int r=0;r<25;r++)
//          {
//             ttreer=r;//0.5*(2*r*r+1);
//             ttreepol=phi;
//             ttreedose=edepPrimaryH->GetBinContent(phi,r);
//             _angle1->Fill();
//          }

gStyle->SetOptStat(0.1);
// _angle1->Draw("p:r","","pol surf2 C Z");
edepPrimaryH->Smooth(1);
edepPrimaryH->Draw("pol surf2 C");
edepPrimaryH->GetYaxis()->SetRangeUser(0, 3);
// cout<<edepPrimaryH->GetBinContent(0,0)<<endl;
c1->cd();
c1->Update();
}




void plot() 
{
  double _total_mass_attenuation_coefficient_0_1MeV=0.167;
  double _total_mass_attenuation_coefficient_0_3MeV=0.118;
  double _total_mass_attenuation_coefficient_0_5MeV=0.0966;
  double _total_mass_attenuation_coefficient_1_25MeV=0.0630;
  double _total_mass_attenuation_coefficient_0_8MeV=0.0786;
  double _total_mass_attenuation_coefficient_1_0MeV=0.0706;
  double _total_mass_attenuation_coefficient_1_5MeV=0.0575;
  double _total_mass_attenuation_coefficient_2_0MeV=0.0493;
  double _total_mass_attenuation_coefficient_3_0MeV=0.0396;
  double _total_mass_attenuation_coefficient_4_0MeV=0.0339;
  double _total_mass_attenuation_coefficient_5_0MeV=0.0301;
  double E0=(1e6-361)*5;
// cout
// <<"--->>Primary BackScatter :"<<Calculate_BackScatter_Edep_Primary(E0)<<endl
// <<"--->>Total BackScatter :"  <<Calculate_BackScatter_Edep_Total(E0)<<endl
// <<"--->>Ftot :"  <<Calculate_Ftot(E0)<<endl
// <<"--->>Fp :"    <<Calculate_Fp_without_eioni(E0)<<endl
// <<"--->>Reff :"  <<Calculate_Reff(E0)<<" cm"<<endl
// <<"--->>Zeff :"  <<Calculate_Zeff(E0)<<" cm"<<endl
// <<"--->>Yeff :"  <<Calculate_Yeff(E0)<<" cm"<<endl
// <<"--->>B11 :"   <<Calculate_B11(E0,_total_mass_attenuation_coefficient_1_25MeV,1)<<endl
// <<"--->>B12 :"   <<Calculate_B12(E0,_total_mass_attenuation_coefficient_1_25MeV,1)<<endl
// <<"--->>B13 :"   <<Calculate_B13(E0,_total_mass_attenuation_coefficient_1_25MeV,1)<<endl;

// <<endl;
// Plot_Kernels();
// Plot_DoseR2_CosTheta();
// Plot_DoseR2_R();
IsoLevelPlot();

//  const Int_t NBINS = 26;
//    Double_t edges[NBINS + 1] = {
//       0.0,
//       // 1e-6,
//       0.5,
//       1.5,
//       4.5,
//       9.5,
//       16.5,
//       25.5,
//       36.5,
//       49.5,
//       64.5,
//       81.5,
//       100.5,
//       121.5,
//       144.5,
//       169.5,
//       196.5,
//       225.5,
//       256.5,
//       289.5,
//       324.5,
//       361.5,
//       400.5,
//       441.5,
//       484.5,
//       529.5,
//       576.5,
//       625.5};

// TH2D* _histo2DPrimarydoesKernels_newbining = new TH2D("2d_PrimaryedepKernels", "2d_PrimaryedepKernels", 96, 0, 96/*,25,0,25);*/, NBINS, edges);

// for(int phi=0;phi<48;phi++)
//    for(int r=0;r<25;r++)
//    {
//       if(r==12 && phi==0)
//       {
//          double R=0.5*(2*r*r+1);
//          double value = 100;
//          cout
//          <<"------------------>> R = "<<R
//          <<endl;
//          _histo2DPrimarydoesKernels_newbining->SetBinContent(phi,r,value);
//       }
//    }



// for(int phi=0;phi<48;phi++)
//    for(int r=0;r<25;r++)
//    {
//      double R=(0.5*(2*r*r+1));
//      if(r==0)
//      cout
//      <<"-->> R = "<<R
//      <<"-->> phi = "<<phi
//      <<"-->> value = "<<_histo2DPrimarydoesKernels_newbining->GetBinContent(phi,0)
//      <<endl;
//      cout
//      <<"-->> R = "<<R
//      <<"-->> phi = "<<phi
//      <<"-->> value = "<<_histo2DPrimarydoesKernels_newbining->GetBinContent(phi,r)
//      <<endl;
//    }





}