#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 <cmath>
#include <ctime>
#include <cstring>
#include <TGraph2D.h>
#include <TF2.h>
#include <TH3.h>
#include <TH3D.h>

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

Root Script to:
- Get the 2D Compton image from the file and plot it.
- Get the 1D projections along the x- and y-axis from the file and plot thme.
- Calculate the reconstructed position of the source and the image resolution from the projections and write them in a text output file named XXX_analysis_image.dat
- Save the plot of the images (2D histograms) and projections (1D histograms) in .png format. 
*/

double FWHM(TH1D *h);

void Plot_Compton(const char *fname){
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);
  //gStyle->SetPalette(1);
  gStyle->SetOptTitle(0);
  const int nIteds {1};
  TH2D *hImage[nIteds];
  TH1D *hProjectionX[nIteds];
  TH1D *hProjectionY[nIteds];
  const Double_t yLim[2]={-20.0,20.0};
  const Int_t nColor[nIteds]={2};
  TFile *f[nIteds];
  //char fname[300];
  for(auto i=0;i<nIteds;i++){
    //if(i==0){sprintf(fname,"%s_Analytical_iTED_s_%d_nDegree_%d_zp0.root",name,0,nDegree);}
    //if(i==1){sprintf(fname,"%s_Analytical_iTED_s_%d_nDegree_%d_zp0.root",name,6,nDegree);}
   // if(i==0){sprintf(fname,"%s",fname);}
  //  if(i==1){sprintf(fname,"%s",fname);}

    int xBinImage, yBinImage, zBinImage;
    float xImage, yImage;

    f[i]=new TFile(fname,"READ");
    f[i]->GetObject("hAnalytical",hImage[i]);
    hImage[i]->SetTitle("Compton Image iTED-A");
    
   hImage[i]->GetMaximumBin(xBinImage, yBinImage, zBinImage);

   xImage = hImage[i]->GetXaxis()->GetBinCenter(xBinImage);
   yImage = hImage[i]->GetYaxis()->GetBinCenter(yBinImage);
   
   // xBinImage=hImage[i]->GetXaxis()->FindBin(0.);
   // yBinImage=hImage[i]->GetYaxis()->FindBin(Ypos);

  //  xBinImage=hImage[i]->GetXaxis()->FindBin(xImage);
   // yBinImage=hImage[i]->GetYaxis()->FindBin(yImage);


    int nBinsCut = 1; // +-nBinsCut will be taken for the cutted projection.

    hProjectionX[i] = hImage[i]->ProjectionX("hAnalyticalx", yBinImage - nBinsCut, yBinImage + nBinsCut, "");
    hProjectionX[i]->GetYaxis()->SetTitleSize(0.05);
    hProjectionX[i]->GetXaxis()->SetTitleSize(0.05);
    hProjectionX[i]->GetYaxis()->SetTitleOffset(0.85);
    hProjectionX[i]->GetXaxis()->SetTitleOffset(0.85);
    hProjectionX[i]->GetXaxis()->SetTitle("x [mm]");
    hProjectionX[i]->GetYaxis()->SetTitle("Arb. Units");
    hProjectionX[i]->GetYaxis()->CenterTitle();
    hProjectionX[i]->SetLineWidth(3);
    hProjectionX[i]->SetLineColor(kRed);

    //  l->AddEntry(hProjectionX[i], "x-central-Projection", "l");
   // l->AddEntry(hProjectionX[i], Form("x-Proj: FWHM = %.1f mm", FWHM(hProjectionX[i])), "l");

    hProjectionY[i] = hImage[i]->ProjectionY("hAnalyticaly", xBinImage - nBinsCut, xBinImage + nBinsCut, "");
    hProjectionY[i]->GetYaxis()->SetTitleSize(0.05);
    hProjectionY[i]->GetXaxis()->SetTitleSize(0.05);
    hProjectionY[i]->GetYaxis()->SetTitleOffset(0.85);
    hProjectionY[i]->GetXaxis()->SetTitleOffset(0.85);
    hProjectionY[i]->GetXaxis()->SetTitle("y [mm]");
    hProjectionY[i]->GetYaxis()->SetTitle("Arb. Units");
    hProjectionY[i]->GetYaxis()->CenterTitle();
    hProjectionY[i]->SetLineWidth(3);
    hProjectionY[i]->SetLineColor(kBlue);
    //   l->AddEntry(hAnalyticaly, "y-central-Projection", "l");
   // l->AddEntry(hAnalyticaly, Form("x-Proj: FWHM = %.1f mm", FWHM(hAnalyticaly)), "l");tf(fname,"Projection_%d",i);

   //f[i]->GetObject("hAnalyticalx",hProjectionX[i]);
   // f[i]->GetObject("hAnalyticaly",hProjectionY[i]);

    /*hProjectionX[i]->SetLineColor(nColor[i]);
    hProjectionX[i]->SetLineWidth(3);
    hProjectionX[i]->GetXaxis()->SetTitle("X [mm]");
    hProjectionX[i]->GetYaxis()->SetTitle("Arb. units");
    hProjectionX[i]->GetYaxis()->CenterTitle();
    hProjectionX[i]->GetXaxis()->SetTitleSize(0.05);
    hProjectionX[i]->GetYaxis()->SetTitleSize(0.05);
    hProjectionX[i]->GetYaxis()->SetTitleOffset(0.85);
    hProjectionX[i]->GetXaxis()->SetTitleOffset(0.85);*/
    //sprintf(fname,"X-Projection %3.1f < y [mm] < %3.1f",yLim[0],yLim[1]);
   // hProjectionX[i]->SetTitle(fname);
    //hProjectionX[i]->GetYaxis()->SetMaxDigits(3);


  }

int bin1_x= hProjectionX[0]-> FindFirstBinAbove(hProjectionX[0]->GetMaximum()/2);
cout<<" First Bin above half-maximum X-projection:"<< bin1_x<<endl;
cout<<" X-coordinate first bin above half-maximum X-projection:"<< hProjectionX[0]->GetBinCenter(bin1_x)<<endl;
int bin2_x= hProjectionX[0]-> FindLastBinAbove(hProjectionX[0]->GetMaximum()/2);
cout<<" Last Bin above half-maximum X-projection:"<< bin2_x<<endl;
cout<<" X-coordinate last bin above half-maximum X-projection:"<< hProjectionX[0]->GetBinCenter(bin2_x)<<endl;

//double fwhm_x= hProjectionX[0]->GetBinCenter(bin2_x)-hProjectionX[0]->GetBinCenter(bin1_x);
double fwhm_x=FWHM(hProjectionX[0]);
double centroid_x=hProjectionX[0]-> GetBinCenter(hProjectionX[0]->GetMaximumBin());
cout<< "FWHM [mm] x-projection=  "<<fwhm_x<<endl;
cout<< "Centroid [mm] x-projection=  "<<centroid_x<<endl;

int bin1_y= hProjectionY[0]-> FindFirstBinAbove(hProjectionY[0]->GetMaximum()/2);
cout<<" First Bin above half-maximum Y-projection:"<< bin1_y<<endl;
cout<<" Y-coordinate first bin above half-maximum Y-projection:"<< hProjectionY[0]->GetBinCenter(bin1_y)<<endl;
int bin2_y= hProjectionY[0]-> FindLastBinAbove(hProjectionY[0]->GetMaximum()/2);
cout<<" Last Bin above half-maximum Y-projection:"<< bin1_y<<endl;
cout<<" Y-coordinate last bin above half-maximum Y-projection:"<< hProjectionY[0]->GetBinCenter(bin1_y)<<endl;
//double fwhm_y= hProjectionY[0]->GetBinCenter(bin2_y)-hProjectionY[0]->GetBinCenter(bin1_y);
double fwhm_y=FWHM(hProjectionY[0]);
double centroid_y=hProjectionY[0]-> GetBinCenter(hProjectionY[0]->GetMaximumBin());
cout<< "FWHM [mm] y-projection=  "<<fwhm_y<<endl;
cout<< "Centroid [mm] y-projection=  "<<centroid_y<<endl;

ofstream tout(Form("Plots_Compton/%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.close();

  auto *c1=new TCanvas("c1","c1",1400,700);
  c1->Divide(2);
  c1->cd(1);
  hImage[0]->Draw("cont4");
  hImage[0]->GetXaxis()->SetRangeUser(-175.0,175.0);
  hImage[0]->GetYaxis()->SetRangeUser(-175.0,175.0);
  c1->cd(2);
  hProjectionX[0]->GetXaxis()->SetRangeUser(-175.0,175.0);
  hProjectionX[0]->Draw("histo");
  hProjectionY[0]->Draw("histosame");
  TLegend* leg = new TLegend(.55, .75, .9, .9);
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);
  leg->AddEntry(hProjectionX[0],Form("FWHM-X (mm): %.1f",fwhm_x),"l");
  leg->AddEntry(hProjectionY[0],Form("FWHM-Y (mm): %.1f",fwhm_y),"l");

  leg->Draw();

c1->SaveAs(Form("Plots_Compton/%s_ComptonImage_and_X_Y_Projections.root",fname));
c1->SaveAs(Form("Plots_Compton/%s_ComptonImage_and_X_Y_Projections.png",fname));
/*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 FWHM(TH1D *h) {
    bool find1 = false;

    double x1 = 0;
    double x2 = 0;

    for (int i = 0; i < h->GetNbinsX(); i++) {
        if ((h->GetBinContent(i + 1) >= h->GetMaximum() / 2) && (!find1)) {
            x1 = h->GetBinCenter(i + 1);
           // if (verbose == 1) cout << " Previous x1: " << x1 << endl;
            for (float x = h->GetBinCenter(i); x <= h->GetBinCenter(i + 1); x += 0.1) {
            //    if (verbose == 1)cout << " Half Max Value :" << 0.5 * h->GetMaximum() << " --> Value for x: " << x << " --> " << h->Interpolate(x) << endl;
                if (h->Interpolate(x) >= h->GetMaximum() / 2) x1 = x - 0.05;
            }

            find1 = true;
        }
        if ((i + 1 > h->GetMaximumBin())) {
            if (h->GetBinContent(i + 1) <= h->GetMaximum() / 2) {
                x2 = h->GetBinCenter(i + 1);
          //      if (verbose == 1)cout << " Previous x2: " << x2 << endl;
                for (float x = h->GetBinCenter(i); x <= h->GetBinCenter(i + 1); x += 0.1) {
                //    if (verbose == 1)cout << " Half Max Value :" << 0.5 * h->GetMaximum() << " --> Value for x: " << x << " --> " << h->Interpolate(x) << endl;
                    if (h->Interpolate(x) <= h->GetMaximum() / 2) x2 = x - 0.05;
                }
                break;
            } else if (i + 1 == h->GetNbinsX()) {
                x2 = h->GetBinCenter(h->GetNbinsX());
                break;
            }
        }
    }
    //getchar();
   // cout << "\nx1: " << x1 << endl;
  //  cout << " x2: " << x2 << endl;
   // cout << " FWHM: " << x2 - x1 << endl;
    return x2 - x1;
}




