#include <TMath.h>
#include <TRandom.h>
#include <TRandom3.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH1D.h>
#include <TH2F.h>
#include <TH2.h>
#include <TLegend.h>
#include <TROOT.h>
#include <TStyle.h>
#include <string>
#include <TTree.h>
#include <TFile.h>
#include <cstdlib>
#include <iostream>
#include <TF1.h>
#include <fstream>
#include <ctime>
#include <cstring>
#include <TGraph2D.h>
#include <TF2.h>

/*J. Lerendegui Marco:
TEAFN, Máster Interuniversitario Física Nuclear 2025/26, IFIC (Valencia), 2-6 Febrero 2026.

Root Script to:
- read the individual PET Images (25 coincidences between the 5 crystals of iTED-A and the 5 crystals of iTED-B:  0-5, 0-6, 0-7, 0-8, 0-9, 1-5,  ...,  4-9.)
- sum the images in a single image
- get the projections along the x- and y-axis
- calculate the reconstructed position of the source, the image resolution from the projections and the number of events and fraction of events selected within the energy window and write them in a text output file named XXX_analysis_image.dat
- save the images (2D histograms) and projections (1D histograms) in a root file and the figures in .png format. 
*/

void Plot_PET(const char *fname){
  gROOT->Reset();
  //gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);
  //gStyle->SetPalette(1);
  gStyle->SetOptTitle(0);
  char name[100];
  const int nPETImages {25};
  const Double_t yLim[2] {-60.0,60.0};// {-60.0,60.0};
  const Double_t xLim[2] {-60.0,60.0};//{-60.0,60.0};
  auto *f=new TFile(fname,"READ");
  TH2D *hPET[nPETImages];
  TH1D *hEnergy[nPETImages];
  TH1D *hEnergy_selected[nPETImages];
  for(auto j=0;j<nPETImages;j++){
    sprintf(name,"h%d_Image",j+1);
    f->GetObject(name,hPET[j]);
    sprintf(name,"h%d_Selected_Deposited_Energy",j+1);
    f->GetObject(name,hEnergy_selected[j]);
    hEnergy_selected[j]->SetLineColor(kRed);
    hEnergy_selected[j]->SetLineWidth(3);
    sprintf(name,"h%d_Selected_Total_Energy",j+1);
    f->GetObject(name,hEnergy[j]);
    hEnergy[j]->SetLineColor(kBlue);
    hEnergy[j]->SetLineWidth(3);
  }

  for(auto i=1;i<nPETImages;i++){
    //if (i>50 && i<76)
    //{
    hPET[0]->Add(hPET[i],1.0e0);
    hEnergy[0]->Add(hEnergy[i],1.0e0);
    hEnergy_selected[0]->Add(hEnergy_selected[i],1.0e0);
    cout<<"Adding up "<<i<<" histogram....."<<endl;
    //}
  }
hPET[0]->Rebin2D(1);

  TH1D *hPETx=nullptr,*hPETy=nullptr;
  sprintf(name,"hPETx");
  hPETx=hPET[0]->ProjectionX(name,hPET[0]->GetYaxis()->FindBin(yLim[0]),hPET[0]->GetYaxis()->FindBin(yLim[1]),"[cutg]");
  hPETx->SetLineColor(kRed);
  hPETx->SetLineWidth(3);
  hPETx->GetXaxis()->SetTitle("x [mm]");
  hPETx->GetYaxis()->SetTitle("Counts");
  hPETx->GetYaxis()->CenterTitle();
  hPETx->GetXaxis()->SetTitleSize(0.05);
  hPETx->GetYaxis()->SetTitleSize(0.05);
  hPETx->GetYaxis()->SetTitleOffset(0.85);
  hPETx->GetXaxis()->SetTitleOffset(0.85);
  //hPETx->Scale(1.0/hPETx->GetMaximum());
  sprintf(name,"hPETy");
  hPETy=hPET[0]->ProjectionY(name,hPET[0]->GetXaxis()->FindBin(xLim[0]),hPET[0]->GetXaxis()->FindBin(xLim[1]),"[cutg]");

  hPETy->SetLineColor(kBlue);
  hPETy->SetLineWidth(3);
  hPETy->GetXaxis()->SetTitle("x [mm]");
  hPETy->GetYaxis()->SetTitle("Counts");
  hPETy->GetYaxis()->CenterTitle();
  hPETy->GetXaxis()->SetTitleSize(0.05);
  hPETy->GetYaxis()->SetTitleSize(0.05);
  hPETy->GetYaxis()->SetTitleOffset(0.85);
  hPETy->GetXaxis()->SetTitleOffset(0.85);
  hPETy->Scale(hPETx->GetMaximum()/hPETy->GetMaximum());


// Quantitative analysis:
int bin1_x= hPETx-> FindFirstBinAbove(hPETx->GetMaximum()/2);
int bin2_x= hPETx-> FindLastBinAbove(hPETx->GetMaximum()/2);
double fwhm_x= hPETx->GetBinCenter(bin2_x)-hPETx->GetBinCenter(bin1_x);
double centroid_x=hPETx-> GetBinCenter(hPETx->GetMaximumBin());
cout<< "FWHM [mm] x-projection=  "<<fwhm_x<<endl;
cout<< "Centroid [mm] x-projection=  "<<centroid_x<<endl;
int bin1_y= hPETy-> FindFirstBinAbove(hPETy->GetMaximum()/2);
int bin2_y= hPETy-> FindLastBinAbove(hPETy->GetMaximum()/2);
double fwhm_y= hPETy->GetBinCenter(bin2_y)-hPETy->GetBinCenter(bin1_y);
double centroid_y=hPETy-> GetBinCenter(hPETy->GetMaximumBin());
cout<< "FWHM [mm] y-projection=  "<<fwhm_y<<endl;
cout<< "Centroid [mm] y-projection=  "<<centroid_y<<endl;
double total_counts_x= hEnergy[0]->Integral();
double counts_energy_selected=hEnergy_selected[0]->Integral();
cout<< "Num. Energy counts=  "<<total_counts_x<<endl;
cout<< "Num. Energy selected counts=  "<< counts_energy_selected<< endl;
cout<< "Percent. inside energy window=  " << (counts_energy_selected/total_counts_x)*100<<endl;


ofstream tout(Form("Plots_PET/%s_analysis_image.dat",fname));
// Quantitative analysis:
tout<< "FWHM [mm] x-projection =  "<<fwhm_x<<endl;
tout<< "Centroid [mm] x-projection =  "<<centroid_x<<endl;
tout<< "FWHM [mm] y-projection =  "<<fwhm_y<<endl;
tout<< "Centroid [mm] y-projection =  "<<centroid_y<<endl;
tout<< "No. counts =  "<<total_counts_x<<endl;
tout<< "No. selected counts energy window =  "<< counts_energy_selected<< endl;
tout<< "% inside energy window =  " << (counts_energy_selected/total_counts_x)*100<<endl;

tout.close();


  auto *c2=new TCanvas("c2","c2",700,700);
  hPET[0]->GetXaxis()->SetRangeUser(xLim[0],xLim[1]);
  hPET[0]->GetYaxis()->SetRangeUser(yLim[0],yLim[1]);
  hPET[0]->Draw("colz");
  c2->SaveAs(Form("Plots_PET/%s_SumPETImage.png",fname));

  auto *c1=new TCanvas("c1","c1",700,700);
  hPETx->GetXaxis()->SetRangeUser(xLim[0],xLim[1]);
  hPETx->Draw("histo");
  hPETy->Draw("histo same");
  TLegend* leg = new TLegend(.6, .7, .9, .9);
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);
  leg->AddEntry(hPETx,Form("FWHM-X[mm]: %.1f",fwhm_x),"l");
  leg->AddEntry(hPETy,Form("FWHM-Y[mm]: %.1f",fwhm_y),"l");
  leg->Draw();
  c1->SaveAs(Form("Plots_PET/%s_SumPETImage_X_Y_Projections.png",fname));

  //hRange->Draw("histo same");
  auto *c3=new TCanvas("c3","c3",900,675);
  hEnergy[0]->Draw("histo");
  hEnergy_selected[0]->Draw("histo same");
  TLegend* leg2 = new TLegend(.6, .7, .9, .9);
  leg2->SetFillStyle(0);
  leg2->SetBorderSize(0);
  leg2->AddEntry(hEnergy[0],"E_{A} + E_{B}","l");
  leg2->AddEntry(hEnergy_selected[0],"E_{A} + E_{B} selected PET","l");
  leg2->Draw();
  c3->SaveAs(Form("Plots_PET/%s_SumEnergySpectra_SelectedWindow.png",fname));

TFile *fout = new TFile(Form("%s_Sum.root",fname),"RECREATE");
fout->cd();
hPET[0]->Write("hPET_Sum");
hPETx->Write("hPETx_Sum");
hPETy->Write("hPETy_Sum");
hEnergy[0]->Write("hEnergy_Sum");
hEnergy_selected[0]->Write("hEnergy_selected_Sum");
fout->Close();

}
