parametric_fit/src/FitFunctions.cc

0001 #include "CustomHistogramFit/FitFunctions.h"
0002 #include <iostream>
0003 #include <cmath>
0004
0005 namespace FitFunctions {
0006
0007   // Implementation of the fitting functions
0008
0009   double gauss(const double* x, const double* par) {
0010     double shape =
0011       par[0] * std::exp(-0.5 * std::pow((x[0] - par[1]) / par[2], 2));
0012     return shape;
0013   }
0014
0015   double poly2(const double* x, const double* par) {
0016     // Second order polynomial parametrized in terms of Bernstein polynomials for the background
0017     double poly =
0018       par[0] * std::pow(1 - x[0], 2) +
0019       par[1] * 2 * x[0] * (1 - x[0]) +
0020       par[2] * std::pow(x[0], 2);
0021     return poly;
0022   }
0023
0024   double poly3(const double* x, const double* par) {
0025     // Third order polynomial parametrized in terms of Bernstein polynomials for the background
0026     double poly =
0027       par[0] * std::pow(1 - x[0], 3) +
0028       par[1] * 3 * x[0] * std::pow(1 - x[0], 2) +
0029       par[2] * 3 * std::pow(x[0], 2) * (1 - x[0]) +
0030       par[3] * std::pow(x[0], 3);
0031     return poly;
0032   }
0033
0034   double poly5(const double* x, const double* par) {
0035     // Fifth order polynomial parametrized in terms of Bernstein polynomials for the background
0036     double poly =
0037       par[0] * std::pow(1 - x[0], 5) +
0038       par[1] * 5 * x[0] * std::pow(1 - x[0], 4) +
0039       par[2] * 10 * std::pow(x[0], 2) * std::pow(1 - x[0], 3) +
0040       par[3] * 10 * std::pow(x[0], 3) * std::pow(1 - x[0], 2) +
0041       par[4] * 5 * std::pow(x[0], 4) * (1 - x[0]) +
0042       par[5] * std::pow(x[0], 5);
0043     return poly;
0044   }
0045
0046   double argusRoofit(const double* x, const double* par)   {
0047     // Argus Background just as it is defined in Roofit
0048     double argusBG =
0049       par[0] * x[0] *
0050     std::pow(1 - std::pow(x[0] / par[1], 2), par[3]) *
0051       std::exp(par[2] * std::pow(1 - x[0] / par[1], 2));
0052     return argusBG;
0053   }
0054
0055   double argusModified(const double* x, const double* par) {
0056     // Revisited formula for the Argus background
0057     // Observable
0058     double m = x[0]; // Input invariant mass
0059
0060     // Parameters
0061     double M = par[0]; // Maximum of the background
0062     double m_offset = par[1]; // start of the shape
0063     double m_max = par[2]; // x-coord for which the maximum is reached
0064     double p = par[3]; // Defines how much the background varies as x changes. The smaller p the flatter.
0065     double q = par[4]; // Defines how the background increases at low x. The smaller q, the steeper increase
0066
0067     // Shift invariant mass wrt the offset
0068     double m_shift = std::abs(m - m_offset);
0069     double m_max_shift = std::max(std::abs(m_max - m_offset), 1e-6);
0070
0071     // Compute the background
0072     double argusBG =
0073       M * std::exp(q * std::log(m_shift / m_max_shift) +
0074                    q / p * (1 - std::pow(m_shift / m_max_shift, p)));
0075     return argusBG;
0076   }
0077
0078   // Modified exponential function
0079   double chat1(const double* x, const double* par) {
0080     // Observable
0081     double M = x[0];
0082
0083     // Parameters
0084     double A = par[0]; // Normalization constant
0085     double Mth = par[1]; // Mass threshold
0086     double B = par[2]; // Control for the steepness of the initial increase
0087     double C = par[3]; // Control for the steepness of the exponential decrease
0088
0089     double result =
0090       A * std::exp(B * std::log(M - Mth) - C * ( M - Mth));
0091     return result;
0092   }
0093
0094   // Threshold function
0095   double chat2(const double* x, const double* par) {
0096     // Observable
0097     double M = x[0];
0098
0099     // Parameters
0100     double A = par[0]; // Normalization constant
0101     double Mth = par[1]; // Mass threshold
0102     double B = par[2]; // Control for the steepness of the initial increase
0103     double C = par[3]; // Control for the steepness of the exponential decrease
0104
0105     double result =
0106       A * (1 - std::exp(- (M - Mth) / B )) * std::exp(- C * (M - Mth));
0107     return result;
0108   }
0109
0110   // Polynomial function with threshold
0111   double chat3(const double* x, const double* par) {
0112     // Observable
0113     double M = x[0];
0114
0115     // Parameters
0116     double A = par[0]; // Normalization constant
0117     double Mth = par[1]; // Mass threshold
0118     double B = par[2]; // First shape parameter
0119     double C = par[3]; // Second shape parameter
0120     double D = par[4]; // Third shape parameter
0121
0122     double result =
0123       A * std::pow(M - Mth, B) * (1 + C * (M - Mth) + D * std::pow(M - Mth, 2));
0124     return result;
0125   }
0126
0127   // Table of models
0128   const std::map<FitModel, modelStructure> models = {
0129     {FitModel::GAUSS, {gauss, "Gaussian function", 3, {"Norm", "mean", "sigma"}}},
0130     {FitModel::POLY2, {poly2, "Second order polynomial (Bernstein)", 3, {"b0", "b1", "b2"}}},
0131     {FitModel::POLY3, {poly3, "Third order polynomial (Bernstein)", 4, {"b0", "b1", "b2", "b3"}}},
0132     {FitModel::POLY5, {poly5, "Fifth order polynomial (Bernstein)", 6, {"b0", "b1", "b2", "b3", "b4", "b5"}}},
0133     {FitModel::ARGUSROOFIT, {argusRoofit, "Argus Background (same as RooFit)", 4, {"N", "m0", "c", "p"}}},
0134     {FitModel::ARGUSMODIFIED, {argusModified, "Argus Background (revisited)", 5, {"M", "m_offset", "m_max", "p", "q"}}},
0135     {FitModel::CHAT1, {chat1, "Modified exponential function", 4, {"A", "Mth", "B", "C"}}},
0136     {FitModel::CHAT2, {chat2, "Threshold function", 4, {"A", "Mth", "B", "C"}}},
0137     {FitModel::CHAT3, {chat3, "Polynomial function with threshold", 5, {"A", "Mth", "B", "C", "D"}}}};
0138
0139   // Function to get model information
0140   const modelStructure& GetModel(FitModel model) {
0141     auto it = models.find(model);
0142     if (it == models.end()) {
0143       std::cerr << "Error in FitFunctions::GetModel: Requested fit model does not exist.\n";
0144       exit(1);
0145     }
0146     return it->second;
0147   }
0148
0149 } // namespace FitFunctions
0150
0151
0152