#include <TFile.h>
#include <TTree.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TStyle.h>
#include <TPaveStats.h>
#include <iostream>
#include <vector>
void GLobalTimeVsPhotonCount() {
gStyle->SetOptStat(0);
gStyle->SetOptFit(1111);
gStyle->SetFitFormat("6.3g");
gStyle->SetStatFormat("6.3g");
TFile *file = TFile::Open("Outputtest.root");
if (!file || file->IsZombie()) {
std::cerr << "Error: Could not open file." << std::endl;
return;
}
TTree *tree = dynamic_cast<TTree*>(file->Get("Photon_Hit"));
if (!tree) {
std::cerr << "Error: Could not find tree Photon_Hit." << std::endl;
file->Close();
return;
}
// Canvas setup with improved styling
TCanvas *c1 = new TCanvas("c1", "Time Spectrum", 1400, 800);
c1->SetMargin(0.12, 0.05, 0.12, 0.05); // left, right, bottom, top
gPad->SetGridx();
gPad->SetGridy();
gPad->SetLogy();
double globalTime;
tree->SetBranchAddress("GlobalTime_ns", &globalTime);
TH1F *hist = new TH1F("hist", "Time Spectrum",350, -5, 380);
hist->GetXaxis()->SetTitle("Time [ns]");
hist->GetYaxis()->SetTitle("Number of Counts");
hist->SetMarkerStyle(20);
hist->SetMarkerSize(0.8);
hist->SetLineColor(kBlue);
// Vector to store time values for statistics
std::vector<double> timeValues;
Long64_t nEntries = tree->GetEntries();
for (Long64_t i = 0; i < nEntries; ++i) {
tree->GetEntry(i);
hist->Fill(globalTime);
timeValues.push_back(globalTime);
}
// Calculate statistics for better initial parameters
std::sort(timeValues.begin(), timeValues.end());
double median = timeValues[timeValues.size()/2];
int maxBin = hist->GetMaximumBin();
double peakPos = hist->GetXaxis()->GetBinCenter(maxBin);
double max_val = hist->GetMaximum();
// Improved fit function with better initial parameters
TF1 *fit_func = new TF1("fit_func",
"[0] * ([5] * (exp(-x/[2]) - exp(-x/[1])) / ([2] - [1]) + "
"(1-[5]) * (exp(-x/[4]) - exp(-x/[3])) / ([4] - [3]))",
-4, 400);
// Set parameter names first for clarity
fit_func->SetParName(0, "Amplitude");
fit_func->SetParName(1, "#tau_{rise1}");
fit_func->SetParName(2, "#tau_{decay1}");
fit_func->SetParName(3, "#tau_{rise2}");
fit_func->SetParName(4, "#tau_{decay2}");
fit_func->SetParName(5, "MixingRatio");
fit_func->SetParName(6, "Background");
// Adjusted parameter limits based on physics
fit_func->SetParLimits(0, max_val*0.1, max_val*5); // Amplitude
fit_func->SetParLimits(1, 0.5, 10); // tau_rise1
fit_func->SetParLimits(2, 15, 40); // tau_decay1
fit_func->SetParLimits(3, 0.1, 5); // tau_rise2
fit_func->SetParLimits(4, 80, 150); // tau_decay2
fit_func->SetParLimits(5, 0.4, 0.9); // mixing ratio
fit_func->SetParLimits(6, 0, max_val*0.05); // background
// Improved fit options
fit_func->SetLineColor(kRed);
fit_func->SetLineWidth(2);
hist->Fit(fit_func, "REMN"); // R=fit range, E=error estimation, M=improved fit, N=no draw
// Draw with improved styling
hist->SetTitle(""); // Remove title, it's redundant with axis labels
hist->GetXaxis()->SetTitleSize(0.045);
hist->GetYaxis()->SetTitleSize(0.045);
hist->GetXaxis()->SetLabelSize(0.04);
hist->GetYaxis()->SetLabelSize(0.04);
hist->Draw("E1"); // E1 for error bars
fit_func->Draw("same");
// Add custom stats box
TPaveText *stats = new TPaveText(0.65, 0.45, 0.89, 0.89, "NDC");
stats->SetFillColor(0);
stats->SetBorderSize(1);
stats->SetTextAlign(12);
stats->AddText(Form("#chi^{2}/ndf = %.3f", fit_func->GetChisquare()/fit_func->GetNDF()));
// Add parameters except MixingRatio and Background
for(int i = 0; i < fit_func->GetNpar(); ++i) {
if (i == 5 || i == 6) continue; // Skip "MixingRatio" and "Background"
stats->AddText(Form("%s = %.3f #pm %.3f",
fit_func->GetParName(i),
fit_func->GetParameter(i),
fit_func->GetParError(i)));
}
stats->Draw();
// Print detailed results, excluding MixingRatio and Background
std::cout << "\nDetailed Fit Results:" << std::endl;
std::cout << "Chi2/NDF = " << fit_func->GetChisquare()/fit_func->GetNDF() << std::endl;
std::cout << "Probability = " << fit_func->GetProb() << std::endl;
for(int i = 0; i < fit_func->GetNpar(); ++i) {
if (i == 5 || i == 6) continue; // Skip "MixingRatio" and "Background"
std::cout << fit_func->GetParName(i) << " = "
<< fit_func->GetParameter(i) << " ± "
<< fit_func->GetParError(i) << std::endl;
}
// Save the result
c1->SaveAs("time_spectrum_fit.pdf");
TFile outFile("delaytime.root", "RECREATE");
c1->Write();
outFile.Close();
file->Close();
delete c1;
}
