introduce random engine, seed, params struct and randomization function

Author: lubynets

PWGHF/D2H/Macros/HFInvMassFitter.cxx

@@ -39,6 +39,7 @@
 #include <TDatabasePDG.h>
 #include <TLine.h>
 #include <TPaveText.h>
+#include <TRandom3.h>
 #include <TStyle.h>
 #include <TVirtualPad.h>
 
@@ -128,12 +129,18 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit,
                                                    mIntegralSgn(0),
                                                    mHistoTemplateRefl(nullptr),
                                                    mDrawBgPrefit(false),
-                                                   mHighlightPeakRegion(false)
+                                                   mHighlightPeakRegion(false),
+                                                   mRandomSeed(-1),
+                                                   mRandomGen(nullptr)
 {
   // standard constructor
   mHistoInvMass = histoToFit;
   mHistoInvMass->SetName("mHistoInvMass");
   mHistoInvMass->SetDirectory(nullptr);
+  if (mRandomSeed >= 0) {
+    mRandomGen = new TRandom3();
+    mRandomGen->SetSeed(mRandomSeed);
+  }
 }
 
 HFInvMassFitter::~HFInvMassFitter()
@@ -988,3 +995,27 @@ void HFInvMassFitter::setTemplateReflections(TH1* histoRefl)
   mHistoTemplateRefl = histoRefl;
   mHistoTemplateRefl->SetName("mHistoTemplateRefl");
 }
+
+double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& parameterRanges)
+{
+  constexpr double DefaultSigmaFraction{10.};
+  constexpr int MaximalNumberOfIterations{20};
+
+  if (mRandomSeed < 0) {
+    return parameterRanges.initial;
+  }
+  const auto sigma = parameterRanges.sigma < 0 ? (parameterRanges.upper - parameterRanges.lower) / DefaultSigmaFraction : parameterRanges.sigma;
+
+  double result;
+  int nIter{0};
+  do {
+    result = mRandomGen->Gaus(parameterRanges.initial, sigma);
+    ++nIter;
+    if (nIter > MaximalNumberOfIterations) {
+      printf("randomizeInitialFitParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma);
+      throw std::runtime_error("");
+    }
+  } while (result < parameterRanges.lower || result > parameterRanges.upper);
+
+  return result;
+}

PWGHF/D2H/Macros/HFInvMassFitter.h

@@ -28,6 +28,7 @@
 #include <TF1.h>
 #include <TH1.h>
 #include <TNamed.h>
+#include <TRandom3.h>
 #include <TVirtualPad.h>
 
 #include <Rtypes.h>
@@ -39,6 +40,14 @@
 
 class HFInvMassFitter : public TNamed
 {
+ private:
+  struct ParameterRanges {
+    double lower{};
+    double upper{};
+    double initial{};
+    double sigma{-999.};
+  };
+
  public:
   enum TypeOfBkgPdf {
     Expo = 0,
@@ -141,6 +150,7 @@ class HFInvMassFitter : public TNamed
   HFInvMassFitter& operator=(const HFInvMassFitter& source);
   void fillWorkspace(RooWorkspace& w) const;
   void highlightPeakRegion(const RooPlot* plot, Color_t color = kGray + 1, Width_t width = 1, Style_t style = 2) const;
+  double randomizeInitialParameter(const ParameterRanges& parameterRanges);
 
   TH1* mHistoInvMass; // histogram to fit
   std::string mFitOption;
@@ -212,6 +222,8 @@ class HFInvMassFitter : public TNamed
   TH1* mHistoTemplateRefl;   /// reflection histogram
   bool mDrawBgPrefit;        /// draw background after fitting the sidebands
   bool mHighlightPeakRegion; /// draw vertical lines showing the peak region (usually +- 3 sigma)
+  int mRandomSeed;           /// seed for random engine for fit's initial parameters randomization
+  TRandom3* mRandomGen;      /// engine for fit's initial parameters randomization
 
   ClassDefOverride(HFInvMassFitter, 1);
 };