sEPDChannelAnalysis/macro/analyze_single_run_FWHM.C

0001 #include <TFile.h>
0002 #include <TH1F.h>
0003 #include <TTree.h>
0004 #include <TSystem.h>
0005 #include <vector>
0006 #include <string>
0007 #include <regex>
0008
0009 // -- RooFit
0010 #include "RooRealVar.h"
0011 #include "RooDataHist.h"
0012 #include "RooLandau.h"
0013 #include "RooGaussian.h"
0014 #include "RooFFTConvPdf.h"
0015 #include "RooPlot.h"
0016 #include "RooFitResult.h"
0017
0018 using namespace RooFit;
0019
0020 const int NUM_CHANNELS = 744;
0021 const double MIN_ENTRIES = 50;
0022
0023 struct ChannelResult {
0024     TH1F* hist = nullptr;
0025     RooPlot* plot = nullptr;
0026     RooRealVar* x = nullptr;
0027     RooRealVar* ml = nullptr;
0028     RooRealVar* sl = nullptr;
0029     RooRealVar* mg = nullptr;
0030     RooRealVar* sg = nullptr;
0031     RooLandau* landau = nullptr;
0032     RooGaussian* gauss = nullptr;
0033     RooFFTConvPdf* pdf = nullptr;
0034     double mpv = -1.0;
0035     double gauss_sigma = -1.0;
0036     double chi2_ndf = -1.0;
0037 };
0038
0039 int extract_run_number(const std::string& filename) {
0040     std::regex pattern("run[_-]?(\\d+).*\\.root");
0041     std::smatch match;
0042     if (std::regex_search(filename, match, pattern)) {
0043         try {
0044             return std::stoi(match[1].str());
0045         } catch (...) {
0046             return -1;
0047         }
0048     }
0049     return -1;
0050 }
0051
0052 ChannelResult process_channel(TH1* hist, int ch) {
0053     ChannelResult result;
0054     if (!hist || hist->GetEntries() < MIN_ENTRIES) return result;
0055
0056     // Find peak bin with ADC > 40 to avoid dark current
0057     int max_bin = -1;
0058     double max_content = -1;
0059     const double ADC_MIN = 40.0;
0060
0061     for(int bin = 1; bin <= hist->GetNbinsX(); bin++) {
0062         double center = hist->GetXaxis()->GetBinCenter(bin);
0063         double content = hist->GetBinContent(bin);
0064         if(center >= ADC_MIN && content > max_content) {
0065             max_content = content;
0066             max_bin = bin;
0067         }
0068     }
0069
0070     if(max_bin == -1) return result;
0071
0072     double x_peak = hist->GetXaxis()->GetBinCenter(max_bin);
0073     double y_peak = hist->GetBinContent(max_bin);
0074
0075     // Calculate FWHM to estimate sigma
0076     double half_max = y_peak/2.0;
0077     double x_left = x_peak, x_right = x_peak;
0078
0079     // Find left half-max point
0080     for(int bin = max_bin; bin >= 1; bin--) {
0081         double center = hist->GetXaxis()->GetBinCenter(bin);
0082         if(center < ADC_MIN) break;
0083
0084         if(hist->GetBinContent(bin) < half_max) {
0085             x_left = center;
0086             break;
0087         }
0088     }
0089
0090
0091     // Find right half-max point
0092     for(int bin = max_bin; bin <= hist->GetNbinsX(); bin++) {
0093         if(hist->GetBinContent(bin) < half_max) {
0094             x_right = hist->GetXaxis()->GetBinCenter(bin);
0095             break;
0096         }
0097     }
0098
0099     // Calculate initial sigma estimate
0100     double fwhm = x_right - x_left;
0101     double sigma_est = fwhm/2.355;  // Convert FWHM to sigma
0102
0103
0104     double fit_lo = x_peak - 2*sigma_est;
0105     double fit_hi = x_peak + 3*sigma_est;
0106
0107     fit_lo = std::max(fit_lo, ADC_MIN);
0108
0109     // Clamp to histogram limits
0110     double h_min = hist->GetXaxis()->GetXmin();
0111     double h_max = hist->GetXaxis()->GetXmax();
0112     fit_lo = std::max(fit_lo, h_min);
0113     fit_hi = std::min(fit_hi, h_max);
0114
0115     // Fallback for bad FWHM estimates
0116     if((fit_hi - fit_lo) < 20) {
0117         fit_lo = x_peak - 40;
0118         fit_hi = x_peak + 40;
0119         fit_lo = std::max(fit_lo, ADC_MIN);
0120         fit_lo = std::max(fit_lo, h_min);
0121         fit_hi = std::min(fit_hi, h_max);
0122     }
0123
0124
0125     result.x = new RooRealVar(Form("x_%d", ch), "ADC", fit_lo, fit_hi);
0126     result.ml = new RooRealVar(Form("ml_%d", ch), "Landau MPV", x_peak, fit_lo, fit_hi);
0127     result.sl = new RooRealVar(Form("sl_%d", ch), "Landau #sigma", 5, 0.1, 50);
0128     result.mg = new RooRealVar(Form("mg_%d", ch), "Gauss #mu", 0, -50, 50);
0129     result.sg = new RooRealVar(Form("sg_%d", ch), "Gauss #sigma", 3, 0.1, 20);
0130
0131     RooDataHist dataHist(Form("data_%d", ch), "Data", *result.x, hist);
0132     result.landau = new RooLandau(Form("landau_%d", ch), "", *result.x, *result.ml, *result.sl);
0133     result.gauss = new RooGaussian(Form("gauss_%d", ch), "", *result.x, *result.mg, *result.sg);
0134
0135     result.x->setBins(4000, "cache");
0136     result.pdf = new RooFFTConvPdf(Form("pdf_%d", ch), "Landau#otimesGauss",
0137                                  *result.x, *result.landau, *result.gauss);
0138
0139     RooFitResult* fit_res = result.pdf->fitTo(dataHist,
0140         Save(true),
0141         PrintLevel(-1),
0142         Range(fit_lo, fit_hi),
0143         NumCPU(4));
0144
0145     if (fit_res && fit_res->status() == 0) {
0146         result.plot = result.x->frame(Title("Pulse Height Spectrum"));
0147         dataHist.plotOn(result.plot);
0148         result.pdf->plotOn(result.plot);
0149
0150         // Calculate MPV using final parameters
0151         const int steps = 500;
0152         double step = (fit_hi - fit_lo)/steps;
0153         double best_x = fit_lo;
0154         double max_val = 0;
0155
0156         for(int i=0; i<steps; i++) {
0157             double x_val = fit_lo + i*step;
0158             result.x->setVal(x_val);
0159             double val = result.pdf->getVal();
0160
0161             if(val > max_val) {
0162                 max_val = val;
0163                 best_x = x_val;
0164             }
0165         }
0166
0167
0168         //chi-squared/ndf
0169
0170         double fit_lo = result.x->getMin();
0171         double fit_hi = result.x->getMax();
0172
0173         int bin_lo = hist->FindBin(fit_lo);
0174         int bin_hi = hist->FindBin(fit_hi);
0175
0176         int nBins = bin_hi - bin_lo + 1;
0177
0178
0179         double chi2 = 0;
0180
0181         int nParams = fit_res->floatParsFinal().getSize();
0182
0183         for(int bin=bin_lo; bin<=bin_hi; bin++) {
0184             double x = hist->GetXaxis()->GetBinCenter(bin);
0185             double data = hist->GetBinContent(bin);
0186             double error = hist->GetBinError(bin);
0187             result.x->setVal(x);
0188             double model = result.pdf->getVal() * hist->GetBinWidth(x) * hist->GetEntries();
0189
0190             if(error > 0) {
0191                 chi2 += pow(data - model, 2) / pow(error, 2);
0192             }
0193         }
0194
0195         int ndf = nBins - nParams;
0196         double chi2_ndf = (ndf > 0) ? chi2/ndf : -1;
0197
0198
0199         result.chi2_ndf = chi2_ndf;
0200
0201         result.mpv = best_x;
0202         result.gauss_sigma = result.sg->getVal();
0203     }
0204
0205     delete fit_res;
0206     return result;
0207 }
0208
0209 void analyze_single_run_FWHM(const char* input_file = "/sphenix/tg/tg01/bulk/ecroft/sEPD/spectra_zvtx_cut/run_47869/run_47869.root_histos.root", const char* output_file = "/sphenix/user/ecroft/sEPDChannelAnalysis/analyze_single_channel_test.root") {
0210     TFile* in_file = TFile::Open(input_file, "READ");
0211     if (!in_file || in_file->IsZombie()) return;
0212
0213     TFile out_file(output_file, "RECREATE");
0214     int run_number = extract_run_number(input_file);
0215
0216     if (run_number >= 0) {
0217         out_file.mkdir(Form("run_%d", run_number));
0218         out_file.cd(Form("run_%d", run_number));
0219
0220         for (int ch = 0; ch < NUM_CHANNELS; ch++) {
0221             TH1F* hist = dynamic_cast<TH1F*>(in_file->Get(Form("channel%d", ch)));
0222             if (!hist) continue;
0223
0224             ChannelResult result = process_channel(hist, ch);
0225             if (!result.plot) continue;
0226
0227             out_file.mkdir(Form("run_%d/ch_%04d", run_number, ch));
0228             out_file.cd(Form("run_%d/ch_%04d", run_number, ch));
0229
0230             hist->Write("original_hist");
0231             result.plot->Write("fit_plot");
0232
0233
0234             TVectorD params(3);
0235             params[0] = result.mpv;
0236             params[1] = result.gauss_sigma;
0237             params[2] = result.chi2_ndf;
0238             params.Write("fit_parameters");
0239
0240             // Cleanup
0241             delete result.plot;
0242             delete result.x;
0243             delete result.ml;
0244             delete result.sl;
0245             delete result.mg;
0246             delete result.sg;
0247             delete result.landau;
0248             delete result.gauss;
0249             delete result.pdf;
0250         }
0251     }
0252
0253     std::cout << "Done!" << std::endl;
0254
0255     in_file->Close();
0256     delete in_file;
0257     out_file.Close();
0258 }