[RF] Update xroofit

Author: will-cern

roofit/xroofit/inc/RooFit/xRooFit/xRooNLLVar.h

@@ -83,7 +83,8 @@ class xRooNLLVar : public std::shared_ptr<RooAbsReal> {
 
    void reinitialize();
 
-   void AddOption(const RooCmdArg &opt);
+   void SetOption(const RooCmdArg &opt);
+   [[deprecated("Use SetOption()")]] void AddOption(const RooCmdArg &opt) { SetOption(opt); }
 
    std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>>
    getData() const; // returns pointer to data and snapshot of globs
@@ -240,11 +241,32 @@ class xRooNLLVar : public std::shared_ptr<RooAbsReal> {
       double fNullVal();
       double fAltVal();
 
+      void setNullVal(double val);
+      void setAltVal(double val);
+      void setObsTS(double val, double err)
+      {
+         obs_ts = val;
+         obs_ts_err = err;
+         fPllType = xRooFit::Asymptotics::Unknown;
+      }
+      void addNullToy(double value, double weight = 1., int seed = 0)
+      {
+         fPllType = xRooFit::Asymptotics::Unknown;
+         nullToys.emplace_back(std::make_tuple(seed, value, weight));
+      }
+      void addAltToy(double value, double weight = 1., int seed = 0)
+      {
+         fPllType = xRooFit::Asymptotics::Unknown;
+         altToys.emplace_back(std::make_tuple(seed, value, weight));
+      }
+
       std::shared_ptr<const RooAbsCollection> coords; // pars of the nll that will be held const alongside POI
 
       std::shared_ptr<const RooFitResult> fUfit, fNull_cfit, fAlt_cfit, fLbound_cfit;
       std::shared_ptr<const RooFitResult> fGenFit; // if the data was generated, this is the fit is was generated from
       bool isExpected = false;                     // if genFit, flag says is asimov or not
+      double obs_ts = std::numeric_limits<double>::quiet_NaN(); // only specified for unknown pll types
+      double obs_ts_err = std::numeric_limits<double>::quiet_NaN();
 
       std::shared_ptr<xRooHypoPoint>
          fAsimov; // same as this point but pllType is twosided and data is expected post alt-fit
@@ -289,7 +311,8 @@ class xRooNLLVar : public std::shared_ptr<RooAbsReal> {
 
       bool AddModel(const xRooNode &pdf, const char *validity = "");
 
-      void LoadFits(const char *apath);
+      // the directory where fits are cached from scans
+      TDirectory *fitCache() const { return fFitDb.get(); }
 
       // A points over given parameter, number of points between low and high
       int AddPoints(const char *parName, size_t nPoints, double low, double high);
@@ -368,22 +391,25 @@ class xRooNLLVar : public std::shared_ptr<RooAbsReal> {
 
       std::set<std::pair<std::shared_ptr<RooArgList>, std::shared_ptr<xRooNode>>> fPdfs;
 
-      std::shared_ptr<TFile> fFitDb;
+      std::shared_ptr<TDirectory> fFitDb;
    };
 
    xRooHypoSpace hypoSpace(const char *parName, int nPoints, double low, double high,
                            double alt_value = std::numeric_limits<double>::quiet_NaN(),
-                           const xRooFit::Asymptotics::PLLType &pllType = xRooFit::Asymptotics::Unknown);
+                           const xRooFit::Asymptotics::PLLType &pllType = xRooFit::Asymptotics::Unknown,
+                           int tsType = 0);
    xRooHypoSpace hypoSpace(const char *parName = "",
                            const xRooFit::Asymptotics::PLLType &pllType = xRooFit::Asymptotics::Unknown,
                            double alt_value = std::numeric_limits<double>::quiet_NaN());
    xRooHypoSpace hypoSpace(int nPoints, double low, double high,
                            double alt_value = std::numeric_limits<double>::quiet_NaN(),
                            const xRooFit::Asymptotics::PLLType &pllType = xRooFit::Asymptotics::Unknown);
-   xRooHypoSpace hypoSpace(const char *parName, xRooFit::TestStatistic::Type tsType, int nPoints = 0)
+   xRooHypoSpace hypoSpace(const char *parName, xRooFit::TestStatistic::Type tsType, int nPoints = 0,
+                           double low = -std::numeric_limits<double>::infinity(),
+                           double high = std::numeric_limits<double>::infinity(),
+                           double alt_value = std::numeric_limits<double>::quiet_NaN())
    {
-      return hypoSpace(parName, int(tsType), nPoints, -std::numeric_limits<double>::infinity(),
-                       std::numeric_limits<double>::infinity());
+      return hypoSpace(parName, nPoints, low, high, alt_value, xRooFit::Asymptotics::Unknown, tsType);
    }
 
    std::shared_ptr<RooArgSet> pars(bool stripGlobalObs = true) const;

roofit/xroofit/inc/RooFit/xRooFit/xRooNode.h

@@ -322,7 +322,7 @@ class xRooNode : public TNamed, public std::vector<std::shared_ptr<xRooNode>> {
    xRooNode Vary(const xRooNode &child);
    xRooNode Constrain(const xRooNode &child);
 
-   xRooNode Combine(const xRooNode &rhs); // combine rhs with this node
+   xRooNode Combine(const xRooNode &rhs, bool silent = false); // combine rhs with this node
 
    xRooNode reduced(const std::string &range = "", bool invert = false)
       const; // return a node representing reduced version of this node, will use the SetRange to reduce if blank
@@ -340,6 +340,8 @@ class xRooNode : public TNamed, public std::vector<std::shared_ptr<xRooNode>> {
    void _ShowVars_(bool set = true); // *TOGGLE* *GETTER=_IsShowVars_
    /** @private */
    bool _IsShowVars_() const;
+   /** @private */
+   void _SetAttribute_(const char *name, const char *value = nullptr); // *MENU*
 
    void SetHidden(bool set = true); // *TOGGLE* *GETTER=IsHidden
    bool IsHidden() const;
@@ -392,7 +394,8 @@ class xRooNode : public TNamed, public std::vector<std::shared_ptr<xRooNode>> {
    GetBinError(int bin, const xRooNode &fr = "", int nToys = 0, bool errorsHi = false, bool errorsLo = false) const;
    std::vector<double> GetBinErrors(int binStart = 1, int binEnd = 0, const xRooNode &fr = "", int nToys = 0,
                                     bool errorsHi = false, bool errorsLo = false) const;
-   std::pair<double, double> IntegralAndError(const xRooNode &fr = "", const char *rangeName = nullptr) const;
+   std::pair<double, double> IntegralAndError(const xRooNode &fr = "", const char *rangeName = nullptr, int nToys = 0,
+                                              bool errorsHi = false, bool errorsLo = false) const;
 
    std::vector<double> GetBinErrorsHi(int binStart = 1, int binEnd = 0, const xRooNode &fr = "", int nToys = 0) const
    {
@@ -413,10 +416,13 @@ class xRooNode : public TNamed, public std::vector<std::shared_ptr<xRooNode>> {
 
    // methods to access default content and error
    double GetContent() const { return GetBinContent(fBinNumber); }
-   double GetError(const xRooNode &fr = "") const
+   double GetError(const xRooNode &fr = "", int nToys = 0, bool errorsHi = false, bool errorsLo = false) const
    {
-      return (fBinNumber == -1) ? IntegralAndError(fr).second : GetBinError(fBinNumber, fr);
+      return (fBinNumber == -1) ? IntegralAndError(fr, "", nToys, errorsHi, errorsLo).second
+                                : GetBinError(fBinNumber, fr, nToys, errorsHi, errorsLo);
    }
+   double GetErrorHi(const xRooNode &fr = "", int nToys = 0) const { return GetError(fr, nToys, true, false); }
+   double GetErrorLo(const xRooNode &fr = "", int nToys = 0) const { return GetError(fr, nToys, false, true); }
    double GetData(const xRooNode &data = "obsData") { return GetBinData(fBinNumber, data); }
 
    // methods to access content and covariances of the CHILDREN of a node
@@ -441,13 +447,13 @@ class xRooNode : public TNamed, public std::vector<std::shared_ptr<xRooNode>> {
             int seed = 0); // generate a dataset from a pdf node using given fr - if none given will use current fit
 
    /** @private */
-   void _fit_(const char *constParValues = "", const char* options = ""); // *MENU*
+   void _fit_(const char *constParValues = "", const char *options = "GoF"); // *MENU*
    /** @private */
    void _generate_(const char *name = "", bool expected = false); // *MENU*
    /** @private */
    void _scan_(const char *what = "plr", double nToys = 0, const char *xvar = "", int nPointsX = 0, double lowX = 0,
                double highX = 0 /*, const char* yvar="", int nBinsY=0, double lowY=0, double highY=0*/,
-               const char *constParValues = "", const char* options = ""); // *MENU*
+               const char *constParValues = "", const char *options = ""); // *MENU*
    //    xRooNode fitTo(const char* datasetName) const;
    //    xRooNode fitTo(const xRooNode& _data) const;
    //    xRooNode generate(bool expected=false) const;

roofit/xroofit/src/Asymptotics.cxx

@@ -49,7 +49,7 @@ double xRooFit::Asymptotics::k(const IncompatFunc &compatRegions, double pValue,
       double operator()(double x) const
       {
          double val = PValue(cFunc, x, poiVal, alt_val, sigma_mu, low, high);
-         if (val < 0)
+         if (val < -1e-10)
             kInvalid = true;
          return val - target;
       }
@@ -154,26 +154,22 @@ double xRooFit::Asymptotics::PValue(const IncompatFunc &compatRegions, double k,
    // go through the 4 'regions' ... only two of which will apply
    if (k <= k_high) {
       out2 += ROOT::Math::gaussian_cdf_c(sqrt(k) + Lambda_y);
-      out += 1.0 -
-             Phi_m(poiVal, poi_primeVal, Lambda_y + sqrt(k), sigma, compatRegions);
+      out += 1.0 - Phi_m(poiVal, poi_primeVal, Lambda_y + sqrt(k), sigma, compatRegions);
    } else {
       double Lambda_high = (poiVal - upBound) * (poiVal + upBound - 2. * poi_primeVal) / (sigma * sigma);
       double sigma_high = 2. * (upBound - poiVal) / sigma;
       out2 += ROOT::Math::gaussian_cdf_c((k - Lambda_high) / sigma_high);
-      out += 1.0 -
-             Phi_m(poiVal, poi_primeVal, (k - Lambda_high) / sigma_high, sigma, compatRegions);
+      out += 1.0 - Phi_m(poiVal, poi_primeVal, (k - Lambda_high) / sigma_high, sigma, compatRegions);
    }
 
    if (k <= k_low) {
       out2 += ROOT::Math::gaussian_cdf_c(sqrt(k) - Lambda_y);
-      out += 1.0 +
-             Phi_m(poiVal, poi_primeVal, Lambda_y - sqrt(k), sigma, compatRegions);
+      out += 1.0 + Phi_m(poiVal, poi_primeVal, Lambda_y - sqrt(k), sigma, compatRegions);
    } else {
       double Lambda_low = (poiVal - lowBound) * (poiVal + lowBound - 2. * poi_primeVal) / (sigma * sigma);
       double sigma_low = 2. * (poiVal - lowBound) / sigma;
       out2 += ROOT::Math::gaussian_cdf_c((k - Lambda_low) / sigma_low);
-      out += 1.0 +
-             Phi_m(poiVal, poi_primeVal, (Lambda_low - k) / sigma_low, sigma, compatRegions);
+      out += 1.0 + Phi_m(poiVal, poi_primeVal, (Lambda_low - k) / sigma_low, sigma, compatRegions);
       /*out +=  ROOT::Math::gaussian_cdf((k-Lambda_low)/sigma_low) +
            2*Phi_m(poiVal,poi_primeVal,(Lambda_low - k_low)==0 ? 0 : ((Lambda_low -
          k_low)/sigma_low),sigma,compatRegions)

roofit/xroofit/src/xRooFit.cxx

@@ -1283,10 +1283,13 @@ std::shared_ptr<const RooFitResult> xRooFit::minimize(RooAbsReal &nll,
          }
       }
 
-      if(miniStrat < _minimizer.fitter()->Config().MinimizerOptions().Strategy() && hesse && out->edm() > _minimizer.fitter()->Config().MinimizerOptions().Tolerance()*1e-3 && out->status() != 3) {
+      if (miniStrat < _minimizer.fitter()->Config().MinimizerOptions().Strategy() && hesse &&
+          out->edm() > _minimizer.fitter()->Config().MinimizerOptions().Tolerance() * 1e-3 && out->status() != 3) {
          // hesse may have updated edm by using a better strategy than used in the minimization
          // so print a warning about this
-         std::cerr << "Warning: post-Hesse edm greater than allowed by tolerance. Consider increasing minimization strategy" << std::endl;
+         std::cerr << "Warning: post-Hesse edm " << out->edm() << " greater than allowed by tolerance "
+                   << _minimizer.fitter()->Config().MinimizerOptions().Tolerance() * 1e-3
+                   << ", consider increasing minimization strategy" << std::endl;
          // Dec24: As this is a new warning, will not update status code for now, so edm will be large
          // but in the future we should probably update the code to 3 so that users don't miss this warning.
          // out->setStatus(3); // edm above max
@@ -1400,7 +1403,7 @@ std::shared_ptr<const RooFitResult> xRooFit::minimize(RooAbsReal &nll,
       }
    }
 
-   if(out && out->status() == 0 && minos) {
+   if (out && out->status() == 0 && minos) {
       // call minos if requested on any parameters
       for (auto label : {"xminos", "xMinos"}) {
          std::unique_ptr<RooAbsCollection> pars(floatPars->selectByAttrib(label, true));