#include <vector>
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TFrame.h"
#include "TH1F.h"
#include "TBenchmark.h"
#include "TRandom.h"
#include "TSystem.h"
#include "TMath.h"
#include "Math/Interpolator.h"

using namespace std;

void code_pfva2025(){
    
    // Root dosyasını aç
    TFile *infile = TFile::Open("PFVA2025_PlastikScintillator.root","READ");
    
    // Histogramları çağır
    TH1F *h_EJ_Cs137 = (TH1F*) infile->Get("hEnergy_EJ_Cs137");
    TH1F *h_EJ_Co60 = (TH1F*) infile->Get("hEnergy_EJ_Co60");
    TH1F *h_EJ_Na22 = (TH1F*) infile->Get("hEnergy_EJ_Na22");
    TH1F *h_EJ_Mn54 = (TH1F*) infile->Get("hEnergy_EJ_Mn54");
    TH1F *h_IstinyePS_Cs137 = (TH1F*) infile->Get("hEnergy_IstinyePS_Cs137");

    TLegend* legend1 = new TLegend(0.1,0.7,0.48,0.9);
    legend1->AddEntry(h_EJ_Na22,"^{22}Na","L");
    legend1->AddEntry(h_EJ_Cs137,"^{137}Cs","L");
    legend1->AddEntry(h_EJ_Mn54,"^{54}Mn","L");
    legend1->AddEntry(h_EJ_Co60,"^{60}Co","L");
    
    TLegend* legend2 = new TLegend(0.1,0.7,0.48,0.9);
    legend2->AddEntry(h_EJ_Cs137,"EJ","L");
    legend2->AddEntry(h_IstinyePS_Cs137,"IstinyePS","L");

    TF1 *f_Cs137 = new TF1("fit_Cs137","gaus",270,385); // 4 Par. Fit
    gStyle->SetOptFit(1111);
    f_Cs137->SetParameter(0,4403.41);
    f_Cs137->SetParameter(1,114.793);
    f_Cs137->SetParameter(2, 21.5557);

    TF1 *f_Co60 = new TF1("fit_Co60","gaus",540,850); // 4 Par. Fit
    gStyle->SetOptFit(1111);
    f_Co60->SetParameter(0,3675.66);
    f_Co60->SetParameter(1,240.137);
    f_Co60->SetParameter(2, 48.8852);
    
    TF1 *f_Na22_1 = new TF1("fit_Na22_1","gaus",185,250); // 4 Par. Fit
    gStyle->SetOptFit(1111);
    f_Na22_1->SetParameter(0,8.61365e+03 );
    f_Na22_1->SetParameter(1,7.90463e+01);
    f_Na22_1->SetParameter(2, 2.02349e+01 );
     
    TF1 *f_Na22_2 = new TF1("fit_Na22_2","gaus",620,800); // 4 Par. Fit
    gStyle->SetOptFit(1111);
    f_Na22_2->SetParameter(0,1.16243e+03  );
    f_Na22_2->SetParameter(1,2.55871e+02);
    f_Na22_2->SetParameter(2, 3.82774e+01 );

    TF1 *f_Mn54 = new TF1("fit_Mn54","gaus",355,480); // 4 Par. Fit
    gStyle->SetOptFit(1111);
    f_Mn54->SetParameter(0,4.93233e+0);
    f_Mn54->SetParameter(1,1.52077e+02);
    f_Mn54->SetParameter(2, 2.66871e+01 );
    
    h_EJ_Co60->Fit("fit_Co60","R");
    h_EJ_Mn54->Fit("fit_Mn54","R");
    h_EJ_Na22->Fit("fit_Na22_1","R");
    h_EJ_Na22->Fit("fit_Na22_2","R+");
    h_EJ_Cs137->Fit("fit_Cs137","R");
    
    TF1 *f_Cs137_2 = new TF1("fit_Cs137_2","gaus",170,300); // 4 Par. Fit
    gStyle->SetOptFit(1111);
    f_Cs137_2->SetParameter(0,4403.41);
    f_Cs137_2->SetParameter(1,114.793);
    f_Cs137_2->SetParameter(2, 21.5557);
    h_IstinyePS_Cs137->Fit("fit_Cs137_2","R");
    
    TCanvas *c1 = new TCanvas("c1","EJ All RS");
    h_EJ_Co60->SetLineColor(8);
    h_EJ_Co60->SetLineWidth(3);
    h_EJ_Cs137->SetLineColor(2);
    h_EJ_Cs137->SetLineWidth(3);
    h_EJ_Na22->SetLineColor(4);
    h_EJ_Na22->SetLineWidth(3);
    h_EJ_Mn54->SetLineWidth(3);
    h_EJ_Na22->GetXaxis()->SetTitle("ADC Channel");
    h_EJ_Na22->GetYaxis()->SetTitle("Event");
    h_EJ_Na22->Draw();
    h_EJ_Mn54->Draw("same");
    h_EJ_Co60->Draw("same");
    h_EJ_Cs137->Draw("same");
    legend1->Draw("same");
    
    TCanvas *c2 = new TCanvas("c2","EJ vs IstinyePS Cs137");
    h_IstinyePS_Cs137->SetLineColor(4);
    h_IstinyePS_Cs137->SetLineWidth(3);
    h_IstinyePS_Cs137->Draw();
    h_EJ_Cs137->Draw("same");
    legend2->Draw("same");
    
    float peak_Cs137 = f_Cs137->GetParameter(1);
    float peak_Co60 = f_Co60->GetParameter(1);
    float peak_Mn54 = f_Mn54->GetParameter(1);
    float peak_Na22_1= f_Na22_1->GetParameter(1);
    float peak_Na22_2= f_Na22_2->GetParameter(1);
    float peak_Cs137_2 = f_Cs137_2->GetParameter(1);
    
    float ymax_Cs137 = f_Cs137->Eval(peak_Cs137);
    float ymax_Co60 = f_Co60->Eval(peak_Co60);
    float ymax_Mn54 = f_Mn54->Eval(peak_Mn54);
    float ymax_Na22_1 = f_Na22_1->Eval(peak_Na22_1);
    float ymax_Na22_2 = f_Na22_2->Eval(peak_Na22_2);
    float ymax_Cs137_2 = f_Cs137_2->Eval(peak_Cs137_2);

    
    float  i_Cs137, i_Co60, i_Na22_1, i_Na22_2, i_Mn54, i_Cs137_2 ;
    
    float ratio;
    float per = 0.80;
    
    cout<<"Gaus peak Cs137 = "<<peak_Cs137<<endl;
    for(float i=peak_Cs137; i<1000; i=i+0.5){
        float y = f_Cs137->Eval(i);
        ratio = y / ymax_Cs137;
        if(ratio<per) {
            i_Cs137 = i;
            cout<<"Ratio Cs137 = "<<ratio<<" at "<<i_Cs137<<endl;
            break;}
        }

    cout<<"Gaus peak Co60 = "<<peak_Co60<<endl;
    for(float i=peak_Co60; i<1000; i=i+0.5){
        float y = f_Co60->Eval(i);
        ratio = y / ymax_Co60;
        if(ratio<per) {
            i_Co60 = i;
            cout<<"Ratio Co60 = "<<ratio<<" at "<<i_Co60<<endl;
            break;}
        }
    
    cout<<"Gaus peak Mn54 = "<<peak_Mn54<<endl;
    for(float i=peak_Mn54; i<1000; i=i+0.5){
        float y = f_Mn54->Eval(i);
        ratio = y / ymax_Mn54;
        if(ratio<per) {
            i_Mn54 = i;
            cout<<"Ratio Mn54 = "<<ratio<<" at "<<i_Mn54<<endl;
            break;}
        }
    
    cout<<"Gaus peak Na22_1 = "<<peak_Na22_1<<endl;
    for(float i=peak_Na22_1; i<1000; i=i+0.5){
        float y = f_Na22_1->Eval(i);
        ratio = y / ymax_Na22_1;
        if(ratio<per) {
            i_Na22_1 = i;
            cout<<"Ratio Na22_1 = "<<ratio<<" at "<<i_Na22_1<<endl;
            break;}
        }
     
     cout<<"Gaus peak Na22_2 = "<<peak_Na22_2<<endl;
     for(float i=peak_Na22_2; i<1000; i=i+0.5){
         float y = f_Na22_2->Eval(i);
         ratio = y / ymax_Na22_2;
         if(ratio<per) {
             i_Na22_2 = i;
             cout<<"Ratio Na22_2 = "<<ratio<<" at "<<i_Na22_2<<endl;
             break;}
         }
     
    cout<<"Gaus peak Cs137 IstiyePS = "<<peak_Cs137_2<<endl;
    for(float i=peak_Cs137_2; i<1000; i=i+0.5){
        float y = f_Cs137_2->Eval(i);
        ratio = y / ymax_Cs137_2;
        if(ratio<per) {
            i_Cs137_2 = i;
            cout<<"Ratio Cs137 IstinyePS = "<<ratio<<" at "<<i_Cs137_2<<endl;
            break;}
        }
    
    double compton[5] = {341., 478., 640., 1038., 1062.}; //Na22, Cs137, Mn54, Co60, Na22;

    double lin[5] = {i_Na22_1, i_Cs137, i_Mn54, i_Co60, i_Na22_2};
    
    TGraph *gr = new TGraph(5, compton, lin);
    TCanvas *c3 = new TCanvas("c3","Linearity");
    gr->SetMarkerStyle(5);
    gr->Draw("AP");
    
    cout<<"{"<<i_Na22_1<<","<<i_Cs137<<","<<i_Mn54<<","<<i_Co60<<","<<i_Na22_2<<"}"<<endl;
    
    float RLY_IstinyePS = 100.*(i_Cs137_2/i_Cs137);
    
    cout<<"Istinye RLY = "<<RLY_IstinyePS<<endl;
    
}