parametric_fit/src/ConstrainedFitter.cc

0001 #include "CustomHistogramFit/ConstrainedFitter.h"
0002
0003 #include "Math/Minimizer.h"
0004 #include "Math/Functor.h"
0005 #include "Math/Factory.h"
0006
0007 #include <algorithm>
0008
0009 ConstrainedFitter::ConstrainedFitter(TH1* hist, std::vector<FitFunctions::FitModel> models)
0010   : hist_(hist) {
0011
0012   modelName_ = "";
0013   singleModelFuns_.resize(models.size());
0014   singleNParams_.resize(models.size());
0015   nParams_ = 0;
0016   int iparGlobal = 0;
0017   for (int im = 0; im < (int) models.size(); im++) {
0018     FitFunctions::FitModel model = models[im];
0019     FitFunctions::modelStructure modelInfo = FitFunctions::GetModel(model);
0020     modelName_ += (im > 0 ? " and " : "") + modelInfo.name;
0021     singleModelFuns_[im] = modelInfo.fun;
0022     singleNParams_[im] = modelInfo.nParams;
0023     nParams_ += modelInfo.nParams;
0024     int ipar = 0;
0025     //std::cout << "singleNParams_[" << im << "] = " << modelInfo.nParams << std::endl;
0026     for (; ipar < modelInfo.nParams; ipar++, iparGlobal++) {
0027       /*std::cout << "paramNames_ = [";
0028       for (int i = 0; i < (int) paramNames_.size() - 1; i++) {
0029         std::cout << paramNames_[i] << ", ";
0030       }
0031       std::cout << paramNames_[(int)paramNames_.size() - 1] << "]\n";*/
0032       if (std::find(paramNames_.begin(), paramNames_.end(), modelInfo.paramNames[ipar]) == paramNames_.end()) {
0033         paramNames_.push_back(modelInfo.paramNames[ipar]);
0034       } else {
0035         paramNames_.push_back(modelInfo.paramNames[ipar] + "_2");
0036       }
0037       //std::cout << "iparGlobal = " << iparGlobal << ", parameter: '" << paramNames_.back() << "'\n";
0038     }
0039     //std::cout << "model '" << modelInfo.name << "' " << modelInfo.nParams << " parameters\n";
0040   }
0041   //std::cout << "word 0\n";
0042   CombinedModelFun_ = [this] (const double *x, const double *par) { return this->CombineModelFunctions(x, par); };
0043   //std::cout << "word 1\n";
0044   //std::cout << "nParams_ = " << nParams_ << std::endl;
0045   initialParams_.resize(nParams_, 0.0);
0046   paramBounds_.resize(nParams_, {0.0, 0.0});
0047   paramStep_ = 0.001;
0048   std::cout << "combined model '" << modelName_ << "' " << nParams_ << " parameters\n";
0049
0050   fitFunction_ = new TF1("fitFunction", CombinedModelFun_, hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax(), nParams_);
0051   for (int ipar = 0; ipar < nParams_; ipar++) fitFunction_->SetParName(ipar, paramNames_[ipar].c_str());
0052 }
0053
0054 ConstrainedFitter::~ConstrainedFitter()
0055 {
0056   delete fitFunction_;
0057 }
0058
0059 void ConstrainedFitter::SetInitialParams(const std::vector<double>& params) {
0060   if (params.size() != static_cast<size_t>(nParams_)) {
0061     std::cerr << "Error in ConstrainedFitter::SetInitialParams: Incorrect number of parameters provided! Expected "
0062               << nParams_ << " parameters for the model '" << modelName_ << "'. "
0063               << params.size() << " given instead." << std::endl;
0064     return;
0065   }
0066   initialParams_ = params;
0067
0068   for (int ipar = 0; ipar < nParams_; ipar++) {
0069     fitFunction_->SetParameter(ipar, initialParams_[ipar]);
0070   }
0071 }
0072
0073 void ConstrainedFitter::SetParamBounds(const std::vector<std::pair<double, double>>& bounds) {
0074   if (bounds.size() != static_cast<size_t>(nParams_)) {
0075     std::cerr << "Error in ConstrainedFitter::SetParamBounds: Incorrect number of bounds provided! Expected "
0076               << nParams_ << " parameters for the model '" << modelName_ << "'. "
0077               << bounds.size() << " given instead." << std::endl;
0078     return;
0079   }
0080   paramBounds_ = bounds;
0081 }
0082
0083 void ConstrainedFitter::PerformFit() {
0084   if (!hist_) {
0085     std::cerr << "Error in ConstrainedFitter::PerformFit: No histogram provided!" << std::endl;
0086     return;
0087   }
0088
0089   ROOT::Math::Minimizer *minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
0090
0091   // Make a lambda function
0092   // Chi2Function cannot be made static
0093   std::function<double(double const*)> const Chi2Lambda = [this] (const double *params) { return this->Chi2Function(params); };
0094
0095   ROOT::Math::Functor functor(Chi2Lambda, nParams_);
0096   minimizer->SetFunction(functor);
0097
0098   // Set:
0099   // - parameter name
0100   // - initial parameter
0101   // - parameter step
0102   // - parameter bounds
0103   for (int i = 0; i < nParams_; i++) {
0104     // std::cout << "minimizer->SetLimitedVariable(" << i << ", \"" << paramNames_[i] << "\", " << initialParams_[i] << ", " << paramStep_ << ", " << paramBounds_[i].first << ", " << paramBounds_[i].second << ");\n";
0105     minimizer->SetLimitedVariable(i, paramNames_[i], initialParams_[i], paramStep_, paramBounds_[i].first, paramBounds_[i].second);
0106   }
0107
0108   bool success = minimizer->Minimize();
0109
0110   if (!success) {
0111     std::cout << "Minimization failed!" << std::endl;
0112   } else {
0113     std::cout << "Minimization successful!" << std::endl;
0114   }
0115   std::cout << "Final Chi2: " << Chi2Function(minimizer->X()) << std::endl;
0116
0117   const double *optimized_params = minimizer->X();
0118   const double *param_errors = minimizer->Errors();
0119
0120   for (int i = 0; i < nParams_; i++)
0121   {
0122     fitFunction_->SetParameter(i, optimized_params[i]);
0123     fitFunction_->SetParError(i, param_errors[i]);
0124   }
0125 }
0126
0127 void ConstrainedFitter::PrintResults() {
0128   if (!fitFunction_) {
0129     std::cerr << "Error: Fit function not initialized! Cannot print results" << std::endl;
0130     return;
0131   }
0132
0133   std::cout << "Fit Results:" << std::endl;
0134   for (int i = 0; i < nParams_; ++i) {
0135     double value, error;
0136     value = fitFunction_->GetParameter(i);
0137     error = fitFunction_->GetParError(i);
0138     std::cout << "Parameter " << i << ": " << value << " +- " << error << std::endl;
0139   }
0140 }
0141
0142 double ConstrainedFitter::CombineModelFunctions(const double* x, const double* par) {
0143   double sum = 0;
0144   int nPars = 0;
0145   int iparGlobal = 0; // parameter index in the combined model
0146   FitFunctions::modelFunction modelFun;
0147   for (int ifun = 0; ifun < (int) singleModelFuns_.size(); ifun++) {
0148     // Iterate over model functions
0149     modelFun = singleModelFuns_[ifun];
0150     nPars = singleNParams_[ifun];
0151     double *params = new double[nPars];
0152     int ipar = 0; // parameter index in the single model
0153     for (; ipar < nPars; ipar++, iparGlobal++) {
0154       params[ipar] = par[iparGlobal];
0155     }
0156     sum += modelFun(x, params);
0157     delete params;
0158   }
0159   return sum;
0160 }
0161
0162 double ConstrainedFitter::Chi2Function(const double* par) {
0163   if (hist_ == nullptr) {
0164     std::cerr << "Error: Histogram not found in minimization function!" << std::endl;
0165     exit(1);
0166   }
0167   fitFunction_->SetParameters(par);
0168
0169   double chi2 = 0;
0170   for (int bin = 2; bin <= hist_->GetNbinsX(); ++bin) {
0171     double x = hist_->GetBinCenter(bin);
0172     double dx = hist_->GetBinWidth(bin);
0173     double y_data_error = hist_->GetBinError(bin);
0174     double y_data = hist_->GetBinContent(bin);
0175     double y_fit = fitFunction_->Eval(x);
0176
0177     // Skip bins outside of the fit range
0178     int nbWindows = fitRange_.size();
0179     // If no fit range is given, consider that the fit is performed over the full histogram
0180     if (nbWindows >= 1) {
0181       bool inFitRange = false;
0182       for (int iw = 0; iw < nbWindows; iw++) {
0183         if ((fitRange_[iw].first < x) && (x < fitRange_[iw].second)) {
0184           inFitRange = true;
0185           break;
0186         }
0187       }
0188       if (!inFitRange) continue;
0189     }
0190
0191     // Decrease error in certain regions to constrain the fit
0192     /*if (x > 0.08 && x < 0.10) // Improve fit on left tail
0193       y_data_error /= 1000;
0194     if (0.18 < x) // Improve fit on right tail
0195       y_data_error /= 1000;
0196       if (y_data > 0.5 * peak_height) // Improve fit on the peak
0197       y_data_error /= 0.01;*/
0198
0199     // Skip empty_bins
0200     if (y_data == 0 || y_data_error == 0) continue;
0201
0202
0203     // Check for reasonable values in fit evaluation
0204     if (std::isnan(y_fit) || std::isinf(y_fit))
0205     {
0206       std::cout << "NaN or Inf detected in fit evaluation at bin " << bin << std::endl;
0207       std::cout << "Parameters were the following:\n";
0208       for (int i = 0; i < nParams_; i++)
0209       {
0210         std::cout << "parameter " << i << ": " << par[i] << std::endl;
0211       }
0212       //exit(1);
0213     }
0214
0215     double increment = std::pow((y_fit - y_data) / y_data_error, 2);
0216     if (std::isnan(increment) || std::isinf(increment))
0217       chi2 += 1e7;
0218     else
0219       chi2 += increment;
0220   }
0221   return chi2;
0222 }
0223