[FITTING] Added fitting methods for histograms with RooFit

tpc_analysis/FITTING/FittingHelpers.hh

@@ -3,6 +3,7 @@
 
 #include "../HEADERS/service_headers.hh"
 #include "../HEADERS/plotting_headers.hh"
+#include "../HEADERS/physics_headers.hh"
 
 #include "RooRealVar.h"
 #include "RooDataHist.h"
@@ -47,6 +48,7 @@ void SaveFitResultsToCSV(const char* csvFilename,
         csvFile << "FitStart,FitEnd,"
                 << "Mpv,MpvError,SigmaL,SigmaLError,"
                 << "sigmaGauss,sigmaGaussError,"
+                << "kPolya,kPolyaError,thetaPolya,thetaPolyaError,"
                 << "Nsig,NsigError,"
                 << "Chi2,NDF,Chi2NDF,PValue,FitStatus,edm\n";
     }
@@ -57,6 +59,8 @@ void SaveFitResultsToCSV(const char* csvFilename,
     RooRealVar* sigma = (RooRealVar*)params.find("sigmaL");
     RooRealVar* sigmaGauss = (RooRealVar*)params.find("sigmaGauss");
     RooRealVar* nsig = (RooRealVar*)params.find("nsig");
+    RooRealVar* kPolya = (RooRealVar*)params.find("kPolya");
+    RooRealVar* thetaPolya = (RooRealVar*)params.find("thetaPolya");
 
     // Write the results to csv format
     csvFile << std::fixed;
@@ -67,10 +71,14 @@ void SaveFitResultsToCSV(const char* csvFilename,
     csvFile << std::setprecision(8);
     if (mpv) csvFile << mpv->getVal() << "," << mpv->getError() << ",";
     else csvFile << "0,0,";
-    if (sigmaL) csvFile << sigmaL->getVal() << "," << sigmaL->getError() << ",";
+    if (sigma) csvFile << sigma->getVal() << "," << sigma->getError() << ",";
     else csvFile << "0,0,";
     if (sigmaGauss) csvFile << sigmaGauss->getVal() << "," << sigmaGauss->getError() << ",";
     else csvFile << "0,0,";
+    if (kPolya) csvFile << kPolya->getVal() << "," << kPolya->getError() << ",";
+    else csvFile << "0,0,";
+    if (thetaPolya) csvFile << thetaPolya->getVal() << "," << thetaPolya->getError() << ",";
+    else csvFile << "0,0,";
 
     csvFile << std::setprecision(2);
     if (nsig) csvFile << nsig->getVal() << "," << nsig->getError() << ",";

tpc_analysis/FITTING/RooLanGaus.cxx → tpc_analysis/FITTING/RooLanGausFit.cxx

@@ -8,23 +8,23 @@ void LanGausFit(const char* inputRootFile, const char* histogramName,
 {
     // Fit range
     double fit_start = 1;    // a.u.
-    double fit_end = 1000;   // a.u.
+    double fit_end = 300;   // a.u.
 
     // Get data histogram
     TH1D* h_signal = FileReadingUtils::Histograms::getTH1D(inputRootFile, histogramName);
     h_signal->SetDirectory(0);
-    //h_signal->Rebin(2);
+    h_signal->Rebin(4);
 
-    RooRealVar charge("charge", "Charge", 0, 2000, "a.u.");
+    RooRealVar charge("charge", "Charge", 1, 500, "a.u.");
     //charge.setBins(1000, "cache");
     RooDataHist data("data", "Data", charge, h_signal);
 
     // Signal Landau distribution
-    RooRealVar mpv("mpv", "Most Probable Value", 200, 50, 400, "a.u.");
-    RooRealVar sigmaL("sigmaL", "Landau Width", 50, 10, 150, "a.u.");
+    RooRealVar mpv("mpv", "Most Probable Value", 200, 20, 400, "a.u.");
+    RooRealVar sigmaL("sigmaL", "Landau Width", 25, 2, 50, "a.u.");
     RooLandau landau("landau", "Landau", charge, mpv, sigmaL);
 
-    RooRealVar nsig("nsig", "Signal Yield", 2000, 10, 5.e5);
+    RooRealVar nsig("nsig", "Signal Yield", 2000, 10, 5.e8);
 
     // Convolution with Gaussian
     RooRealVar sigmaGauss("sigmaGauss", "Gaussian Width", 20, 1, 100, "a.u.");
@@ -55,21 +55,23 @@ void LanGausFit(const char* inputRootFile, const char* histogramName,
 
     // Plot
     RooPlotMaker plotter(charge, "dlc_gain_fit");
-    plotter.SetLegendPosition(0.58, 0.50, 0.9, 0.88);
+    plotter.SetLegendPosition(0.53, 0.60, 0.9, 0.88);
 
-    plotter.CreatePlot("Collected charge (a.u.)", "Counts");
+    plotter.CreatePlot("Collected primary charge (a.u.)", "Counts");
     plotter.AddData(data, "Data", kBlack, 20);
-    plotter.AddPdf(model, "Fit", kRed, 1);
+    plotter.AddPdf(model, "Landau x Gaussian", kRed, 1);
     //plotter.AddPdf(model, "Signal", kBlue, 2, "conv");
     //plotter.AddPdf(model, "Rot Bg", kGreen, 2, "bg_template_pdf");
-    plotter.AddCustomLegendHeader("Cosmic muon data");
+    plotter.AddCustomLegendHeader("#bf{Data:} Cosmic muons");
+    plotter.AddCustomLegendHeader("#bf{HV:} 300 V");
 
     // Add fit results as text
-    plotter.AddText(Form("mpv (a.u.) = %.4f#kern[0.4]{#pm} %.4f", mpv.getVal(), mpv.getError()), 0.58, 0.45, 0.03);
-    plotter.AddText(Form("#sigma_{Land} (a.u.) = %.4f#kern[0.4]{#pm} %.4f", sigmaL.getVal(), sigmaL.getError()), 0.58, 0.40, 0.03);
-    plotter.AddText(Form("N_{sig} = %.0f#kern[0.4]{#pm} %.0f", nsig.getVal(), nsig.getError()), 0.58, 0.35, 0.03);
-    plotter.AddText(Form("#sigma_{Gauss} (a.u.) = %.4f#kern[0.4]{#pm} %.4f", sigmaGauss.getVal(), sigmaGauss.getError()), 0.58, 0.30, 0.03);
-    plotter.AddText(Form("#chi^{2}/ndf = %.2f / %.0f", chi2, ndf), 0.58, 0.25, 0.03);
+    plotter.AddText(Form("Fit range (a.u.): [%.0f, %.0f]", fit_start, fit_end), 0.53, 0.50, 0.03);
+    plotter.AddText(Form("mpv (a.u.) = %.4f#kern[0.4]{#pm} %.4f", mpv.getVal(), mpv.getError()), 0.53, 0.45, 0.03);
+    plotter.AddText(Form("#sigma_{Land} (a.u.) = %.4f#kern[0.4]{#pm} %.4f", sigmaL.getVal(), sigmaL.getError()), 0.53, 0.40, 0.03);
+    plotter.AddText(Form("N_{sig} = %.0f#kern[0.4]{#pm} %.0f", nsig.getVal(), nsig.getError()), 0.53, 0.35, 0.03);
+    plotter.AddText(Form("#sigma_{Gauss} (a.u.) = %.4f#kern[0.4]{#pm} %.4f", sigmaGauss.getVal(), sigmaGauss.getError()), 0.53, 0.30, 0.03);
+    plotter.AddText(Form("#chi^{2}/ndf = %.2f / %.0f", chi2, ndf), 0.53, 0.25, 0.03);
 
     plotter.Draw();
     plotter.SavePlot(outputCanvasName);

tpc_analysis/FITTING/RooLanPolyaFit.cxx

@@ -0,0 +1,95 @@
+#include "FittingHelpers.hh"
+#include "custom_pdf/RooPolya.h"
+using namespace std;
+using namespace RooFit;
+
+// Important: requires RooPolya custom PDF, you need to compile it first
+
+void LanPolyaFit(const char* inputRootFile, const char* histogramName,
+                const char* outputCanvasName = "canvas.pdf",
+                const char* csvFileName = "data_results.csv") 
+{
+    // Fit range
+    double fit_start = 50;    // a.u.
+    double fit_end = 1000;   // a.u.
+
+    // Get data histogram
+    TH1D* h_signal = FileReadingUtils::Histograms::getTH1D(inputRootFile, histogramName);
+    h_signal->SetDirectory(0);
+    h_signal->Rebin(4);
+
+    RooRealVar charge("charge", "Charge", 5, 1000, "a.u.");
+    //charge.setBins(1000, "cache");
+    RooDataHist data("data", "Data", charge, h_signal);
+
+    // Signal Landau distribution
+    RooRealVar mpv("mpv", "Most Probable Value", 200, 10, 400, "a.u.");
+    RooRealVar sigmaL("sigmaL", "Landau Width", 25, 3, 50, "a.u.");
+    RooRealVar offset("offset", "Landau Offset", 100, 50, 600, "a.u.");
+    RooFormulaVar shifted_charge("shifted_charge", "charge + offset", "charge + offset", RooArgList(charge, offset));
+    RooLandau landau("landau", "Landau", shifted_charge, mpv, sigmaL);
+
+
+    RooRealVar nsig("nsig", "Signal Yield", 2000, 10, 5.e8);
+
+    // Convolution with Polya distribution
+    RooRealVar kPolya("kPolya", "Polya shape", 2, 1.01, 40);
+    RooRealVar thetaPolya("thetaPolya", "Polya scale", 10, 0.01, 50.0, "a.u.");
+    RooPolya polya("polya", "Polya", charge, kPolya, thetaPolya);
+    RooFFTConvPdf conv("conv", "LanPolya", charge, landau, polya);
+    //conv.setConvolutionWindow(RooConst(10), RooConst(200.0));
+
+    // Background (?)
+    RooRealVar nbkg("nbkg", "Background Yield", 100, 10, 1.e5);
+    // (...)
+
+    // Model: Landau * Polya (+ Background)
+    RooAddPdf model("model", "Signal", RooArgList(conv), RooArgList(nsig));
+
+    // Fit 
+    RooFitResult* fitResult = model.fitTo(data, Extended(true), Save(), Range(fit_start, fit_end),
+                                          PrintLevel(-1), Minimizer("Minuit2", "Migrad"), Strategy(2));
+
+    // Calculate fit statistics
+    int nParams = model.getParameters(data)->selectByAttrib("Constant", false)->getSize();
+    double ndf = h_signal->FindBin(fit_end) - h_signal->FindBin(fit_start) + 1 - nParams;
+
+    RooAbsReal* chi2Var = model.createChi2(data, Range(fit_start, fit_end));
+    double chi2 = chi2Var->getVal();
+    double pvalue = (ndf > 0) ? 1 - ROOT::Math::chisquared_cdf(chi2, ndf) : 0.0;
+
+    // Save results
+    SaveFitResultsToCSV(csvFileName, fit_start, fit_end, fitResult, chi2, ndf, pvalue);
+
+    // Plot
+    RooPlotMaker plotter(charge, "dlc_gain_fit");
+    plotter.SetLegendPosition(0.53, 0.60, 0.9, 0.88);
+
+    plotter.CreatePlot("Collected primary charge (a.u.)", "Counts");
+    plotter.AddData(data, "Data", kBlack, 20);
+    plotter.AddPdf(model, "Landau x Polya", kRed, 1);
+    //plotter.AddPdf(model, "Signal", kBlue, 2, "conv");
+    //plotter.AddPdf(model, "Rot Bg", kGreen, 2, "bg_template_pdf");
+    plotter.AddCustomLegendHeader("#bf{Data:} Cosmic muons");
+    plotter.AddCustomLegendHeader("#bf{HV:} 340 V");
+
+    // Add fit results as text
+    plotter.AddText(Form("Fit range (a.u.): [%.0f, %.0f]", fit_start, fit_end), 0.53, 0.50, 0.03);
+    plotter.AddText(Form("mpv (a.u.) = %.4f#kern[0.4]{#pm} %.4f", mpv.getVal(), mpv.getError()), 0.53, 0.45, 0.03);
+    plotter.AddText(Form("#sigma_{Land} (a.u.) = %.4f#kern[0.4]{#pm} %.4f", sigmaL.getVal(), sigmaL.getError()), 0.53, 0.40, 0.03);
+    plotter.AddText(Form("N_{sig} = %.0f#kern[0.4]{#pm} %.0f", nsig.getVal(), nsig.getError()), 0.53, 0.35, 0.03);
+    plotter.AddText(Form("k_{Polya} = %.4f#kern[0.4]{#pm} %.4f", kPolya.getVal(), kPolya.getError()), 0.53, 0.30, 0.03);
+    plotter.AddText(Form("#theta_{Polya} (a.u.) = %.4f#kern[0.4]{#pm} %.4f", thetaPolya.getVal(), thetaPolya.getError()), 0.53, 0.25, 0.03);
+    plotter.AddText(Form("offset (a.u.) = %.4f#kern[0.4]{#pm} %.4f", offset.getVal(), offset.getError()), 0.53, 0.20, 0.03);
+    plotter.AddText(Form("#chi^{2}/ndf = %.2f / %.0f", chi2, ndf), 0.53, 0.15, 0.03);
+
+    plotter.Draw();
+    plotter.SavePlot(outputCanvasName);
+    plotter.Clear();
+
+    // Cleanup
+    delete h_signal;
+    delete chi2Var;
+    delete fitResult;
+}
+

tpc_analysis/FITTING/RooPolyaFit.cxx

@@ -0,0 +1,82 @@
+#include "FittingHelpers.hh"
+#include "custom_pdf/RooPolya.h"
+using namespace std;
+using namespace RooFit;
+
+// Important: requires RooPolya custom PDF, you need to compile it first
+
+void PolyaFit(const char* inputRootFile, const char* histogramName,
+                const char* outputCanvasName = "canvas.pdf",
+                const char* csvFileName = "data_results.csv") 
+{
+    // Fit range
+    double fit_start = 0;    // a.u.
+    double fit_end = 200;   // a.u.
+
+    // Get data histogram
+    TH1D* h_signal = FileReadingUtils::Histograms::getTH1D(inputRootFile, histogramName);
+    h_signal->SetDirectory(0);
+    h_signal->Rebin(2);
+
+    RooRealVar charge("charge", "Charge", 0, 200, "a.u.");
+    //charge.setBins(1000, "cache");
+    RooDataHist data("data", "Data", charge, h_signal);
+
+    // Signal Polya distribution
+    RooRealVar kPolya("kPolya", "Polya shape", 2, 1.05, 40);
+    RooRealVar thetaPolya("thetaPolya", "Polya scale", 10, 0.001, 50.0, "a.u.");
+    RooPolya polya("polya", "Polya", charge, kPolya, thetaPolya);
+
+    RooRealVar nsig("nsig", "Signal Yield", 2000, 10, 5.e8);
+
+    // Background (?)
+    RooRealVar nbkg("nbkg", "Background Yield", 100, 10, 1.e5);
+    // (...)
+
+    // Model: Polya (+ Background)
+    RooAddPdf model("model", "Signal", RooArgList(polya), RooArgList(nsig));
+
+    // Fit 
+    RooFitResult* fitResult = model.fitTo(data, Extended(true), Save(), Range(fit_start, fit_end),
+                                          PrintLevel(-1), Minimizer("Minuit2", "Migrad"), Strategy(2));
+
+    // Calculate fit statistics
+    int nParams = model.getParameters(data)->selectByAttrib("Constant", false)->getSize();
+    double ndf = h_signal->FindBin(fit_end) - h_signal->FindBin(fit_start) + 1 - nParams;
+
+    RooAbsReal* chi2Var = model.createChi2(data, Range(fit_start, fit_end));
+    double chi2 = chi2Var->getVal();
+    double pvalue = (ndf > 0) ? 1 - ROOT::Math::chisquared_cdf(chi2, ndf) : 0.0;
+
+    // Save results
+    SaveFitResultsToCSV(csvFileName, fit_start, fit_end, fitResult, chi2, ndf, pvalue);
+
+    // Plot
+    RooPlotMaker plotter(charge, "dlc_gain_fit");
+    plotter.SetLegendPosition(0.53, 0.60, 0.9, 0.88);
+
+    plotter.CreatePlot("Collected primary charge (a.u.)", "Counts");
+    plotter.AddData(data, "Data", kBlack, 20);
+    plotter.AddPdf(model, "Polya", kRed, 1);
+    //plotter.AddPdf(model, "Signal", kBlue, 2, "conv");
+    //plotter.AddPdf(model, "Rot Bg", kGreen, 2, "bg_template_pdf");
+    plotter.AddCustomLegendHeader("#bf{Data:} Laser");
+    plotter.AddCustomLegendHeader("#bf{HV:} 340 V");
+
+    // Add fit results as text
+    plotter.AddText(Form("Fit range (a.u.): [%.0f, %.0f]", fit_start, fit_end), 0.53, 0.50, 0.03);
+    plotter.AddText(Form("k_{Polya} = %.4f#kern[0.4]{#pm} %.4f", kPolya.getVal(), kPolya.getError()), 0.53, 0.45, 0.03);
+    plotter.AddText(Form("#theta_{Polya} (a.u.) = %.4f#kern[0.4]{#pm} %.4f", thetaPolya.getVal(), thetaPolya.getError()), 0.53, 0.40, 0.03);
+    plotter.AddText(Form("N_{sig} = %.0f#kern[0.4]{#pm} %.0f", nsig.getVal(), nsig.getError()), 0.53, 0.35, 0.03);
+    plotter.AddText(Form("#chi^{2}/ndf = %.2f / %.0f", chi2, ndf), 0.53, 0.30, 0.03);
+
+    plotter.Draw();
+    plotter.SavePlot(outputCanvasName);
+    plotter.Clear();
+
+    // Cleanup
+    delete h_signal;
+    delete chi2Var;
+    delete fitResult;
+}
+

tpc_analysis/FITTING/custom_pdf/RooPolya.cxx

@@ -0,0 +1,46 @@
+#include "RooPolya.h"
+#include "RooMath.h"
+#include "TMath.h"
+
+ClassImp(RooPolya)
+
+RooPolya::RooPolya(const char *name, const char *title,
+                   RooAbsReal& _x,
+                   RooAbsReal& _k,
+                   RooAbsReal& _theta)
+    : RooAbsPdf(name, title),
+      x("x", "Observable", this, _x),
+      k("k", "Polya shape", this, _k),
+      theta("theta", "Polya scale", this, _theta)
+{}
+
+
+RooPolya::RooPolya(const RooPolya& other, const char *name)
+    : RooAbsPdf(other, name),
+      x("x", this, other.x),
+      k("k", this, other.k),
+      theta("theta", this, other.theta)
+{}
+
+
+TObject* RooPolya::clone(const char* newname) const {
+    return new RooPolya(*this, newname);
+}
+
+
+Double_t RooPolya::evaluate() const {
+    double xx = x;
+    double kk = k;
+    double th = theta;
+
+    // Check parameters validity
+    if (xx <= 0.0) return 0.0;
+    if (xx < 1.0e-10) xx = 1.0e-10; // Avoid issues at x=0
+    double val =
+        TMath::Power(xx, kk - 1.0) *
+        TMath::Exp(-xx / th) /
+        (TMath::Gamma(kk) * TMath::Power(th, kk));
+
+    return val;
+}
+

tpc_analysis/FITTING/custom_pdf/RooPolya.h

@@ -0,0 +1,39 @@
+#ifndef ROOPOLYA_H
+#define ROOPOLYA_H
+
+#include "RooAbsPdf.h"
+#include "RooRealProxy.h"
+#include "RooAbsReal.h"
+
+// RooPolya: Probability Density Function implementing the Polya (Gamma) distribution
+// Formula:
+//    f(x; k, theta) = (x^(k-1) * exp(-x/theta)) / (Gamma(k) * theta^k)   for x > 0
+// where k > 0 is the shape parameter and theta > 0 is the scale parameter.
+// Gamma(k) is the gamma function evaluated at k.
+
+class RooPolya : public RooAbsPdf {
+public:
+    RooPolya() {}
+
+    RooPolya(const char *name, const char *title,
+             RooAbsReal& _x,
+             RooAbsReal& _k,
+             RooAbsReal& _theta);
+
+    RooPolya(const RooPolya& other, const char* name = 0);
+    virtual TObject* clone(const char* newname) const override;
+
+    inline virtual ~RooPolya() {}
+
+protected:
+    RooRealProxy x;      // variable
+    RooRealProxy k;      // shape parameter
+    RooRealProxy theta;  // scale parameter
+
+    Double_t evaluate() const override;
+
+public:
+    ClassDefOverride(RooPolya,1)
+};
+
+#endif // ROOPOLYA_H