unbinned Maximum Likelihood Method - decay curve fitting

#include "RooRealVar.h" // For variables
#include "RooDataSet.h" // For unbinned datasets
#include "RooExponential.h" // For exponential PDF
#include "RooAddPdf.h" // For adding PDFs (used for extended term)
#include "RooPlot.h" // For plotting
#include "TCanvas.h"
#include "TPad.h"
#include "TRandom3.h"

using namespace RooFit;     // Use RooFit namespace for convenience

// Execute the following in a single Jupyter cell:

// For reproducibility (optional)
gRandom->SetSeed(42);

// Define true yield for generation
int N_true = 100;

// Step 1: Define the observable 't' (time, range 0 to 5)
RooRealVar t("t", "time", 0, 5);  // Variable name, title, min, max

// Step 2: Define the parameter 'tau' (lifetime, initial guess -1, range -5 to -0.1)
RooRealVar tau("tau", "lifetime", -1, -5, -0.1);  // Name, title, initial value, min, max

// Step 3: Create the exponential PDF: f(t|tau) = (1/|tau|) * exp(-t/|tau|) for decay (tau negative)
RooExponential expo("expo", "exponential PDF", t, tau);  // PDF name, title, observable, parameter

// ===================== Data generation =====================
RooRealVar nsig_gen("nsig_gen","signal yield(gen)", N_true);
RooAddPdf  model_gen("model_gen","expo(gen)", RooArgList(expo), RooArgList(nsig_gen));

// Generate toy data with extended term (Poisson-fluctuated total events)
// (In practice, replace with your real data)
RooDataSet* dataset = model_gen.generate(t, Extended(true));  // Generate unbinned dataset

// ===================== EML fit (extended: directly fit yield nsig) =====================
RooRealVar nsig("nsig","signal yield", N_true*0.8, 0.0, N_true*3);
RooAddPdf  model("model","expo (EML)", RooArgList(expo), RooArgList(nsig));

// Extended fit: include Extended(true)
// To switch to unbinned non-extended maximum likelihood fit, remove Extended(true); yield will be treated as normalized to data total
model.fitTo(*dataset, PrintLevel(-1), Extended(true));  // Fit to data; quiet mode

// Step 6: Print the fitted parameter values
tau.Print();  // Outputs the estimated tau with error
nsig.Print(); // Outputs the fitted number of events with error

// Step 7: Plot the data and fitted PDF with a residual (pull) plot below
TCanvas* canvas = new TCanvas("canvas", "Exponential Fit with Pull Plot", 600, 600);  // Larger canvas for two pads

// Split canvas: top pad for main plot, bottom for pulls
TPad* pad1 = new TPad("pad1", "Main plot", 0, 0.3, 1, 1.0);
TPad* pad2 = new TPad("pad2", "Pull plot", 0, 0.0, 1, 0.3);
pad1->Draw();
pad2->Draw();

pad1->cd();  // Switch to top pad
RooPlot* frame = t.frame(Title("Exponential Decay Fit"));  // Create plot frame

// Histogram overlay (with binning); Use DataError(RooAbsData::Poisson) to handle Poisson statistics properly, 
//providing asymmetric error bars for low-count bins and avoiding plotting error bars for empty bins. 
//The default SumW2 assumes Gaussian errors (sqrt(sum w^2)), which can lead to issues like division by zero 
//in pull calculations for empty bins or underestimated uncertainties in low-statistics bins.
dataset->plotOn(frame, Binning(20), DataError(RooAbsData::Poisson));  
model.plotOn(frame, LineColor(kRed));  // Fitted curve
frame->Draw();  // Draw the frame

pad2->cd();  // Switch to bottom pad
RooPlot* frame_pull = t.frame(Title("Pull Distribution"));  // Create pull frame
RooHist* hpull = frame->pullHist();  // Compute Residual: (data - fit)/sigma per bin
frame_pull->addPlotable(hpull, "P");  // Add as points
frame_pull->SetYTitle("Pull");  // Label y-axis
frame_pull->Draw();  // Draw the pull plot

canvas->Draw();

[#1] INFO:Fitting -- RooAbsPdf::fitTo(model) fixing normalization set for coefficient determination to observables in data
[#1] INFO:Fitting -- using generic CPU library compiled with no vectorizations
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_model_model_genData) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
RooRealVar::tau = -1.16743 +/- 0.123584  L(-5 - -0.1) 
RooRealVar::nsig = 98.9814 +/- 9.94  L(0 - 300)