RooFit/xRooFit/xRooFit.h

0001 /*
0002  * Project: xRooFit
0003  * Author:
0004  *   Will Buttinger, RAL 2022
0005  *
0006  * Copyright (c) 2022, CERN
0007  *
0008  * Redistribution and use in source and binary forms,
0009  * with or without modification, are permitted according to the terms
0010  * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
0011  */
0012 #include "Config.h"
0013
0014 // when not using the namespace will use the once pragma.
0015 // when using the namespace (as happens in the ROOT build of xRooFit) then
0016 // will effectively use an include guard
0017 #ifdef XROOFIT_USE_PRAGMA_ONCE
0018 #pragma once
0019 #endif
0020 #if (!defined(XROOFIT_USE_PRAGMA_ONCE) && !defined(XROOFIT_XROOFIT_H)) || \
0021    (defined(XROOFIT_USE_PRAGMA_ONCE) && !defined(XROOFIT_XROOFIT_H_XROOFIT))
0022 #ifndef XROOFIT_USE_PRAGMA_ONCE
0023 #define XROOFIT_XROOFIT_H
0024 #else
0025 #define XROOFIT_XROOFIT_H_XROOFIT
0026 // even with using pragma once, need include guard otherwise cannot include this header
0027 // as part of an interpreted file ... the other headers in xRooFit are similarly affected
0028 // however for now users of xRooFit should only need to include the main xRooFit header to use it all
0029 // in future we should try removing the pragma once altogether (undef XROOFIT_USE_PRAGMA_ONCE)
0030 // and see if it has negative consequences anywhere
0031 #endif
0032
0033 /**
0034  * This is the main include for the xRooFit project.
0035  * Including this should give you access to all xRooFit features
0036  */
0037
0038 class RooAbsData;
0039 class RooAbsCollection;
0040 class RooFitResult;
0041 class RooAbsPdf;
0042 class RooAbsReal;
0043 class RooLinkedList;
0044 class RooWorkspace;
0045
0046 #include "Fit/FitConfig.h"
0047
0048 #include "RooCmdArg.h"
0049 #include "TNamed.h"
0050
0051 class TCanvas;
0052
0053 #include <memory>
0054
0055 BEGIN_XROOFIT_NAMESPACE
0056
0057 class xRooNLLVar;
0058
0059 class xRooFit {
0060
0061 public:
0062    static const char *GetVersion();
0063    static const char *GetVersionDate();
0064
0065    // Extra options for NLL creation:
0066    static RooCmdArg ReuseNLL(bool flag); // if should try to reuse the NLL object when it changes dataset
0067    static RooCmdArg Tolerance(double value);
0068    static RooCmdArg StrategySequence(const char *stratSeq); // control minimization strategy sequence
0069    static constexpr double OBS = std::numeric_limits<double>::quiet_NaN();
0070
0071    // Helper function for matching precision of a value and its error
0072    static std::pair<double, double> matchPrecision(const std::pair<double, double> &in);
0073
0074    // Static methods that work with the 'first class' object types:
0075    //    Pdfs: RooAbsPdf
0076    //    Datasets: std::pair<RooAbsData,const RooAbsCollection>
0077    //    NLLOptions: RooLinkedList
0078    //    FitOptions: ROOT::Fit::FitConfig
0079
0080    // fit result flags in its constPars list which are global observables with the "global" attribute
0081    static std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>>
0082    generateFrom(RooAbsPdf &pdf, const RooFitResult &fr, bool expected = false, int seed = 0);
0083    static std::shared_ptr<const RooFitResult>
0084    fitTo(RooAbsPdf &pdf, const std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>> &data,
0085          const RooLinkedList &nllOpts, const ROOT::Fit::FitConfig &fitConf);
0086    static std::shared_ptr<const RooFitResult> fitTo(RooAbsPdf &pdf,
0087                                                     const std::pair<RooAbsData *, const RooAbsCollection *> &data,
0088                                                     const RooLinkedList &nllOpts, const ROOT::Fit::FitConfig &fitConf);
0089
0090    static xRooNLLVar createNLL(const std::shared_ptr<RooAbsPdf> pdf, const std::shared_ptr<RooAbsData> data,
0091                                const RooLinkedList &nllOpts);
0092    static xRooNLLVar createNLL(RooAbsPdf &pdf, RooAbsData *data, const RooLinkedList &nllOpts);
0093    static xRooNLLVar createNLL(RooAbsPdf &pdf, RooAbsData *data, const RooCmdArg &arg1 = RooCmdArg::none(),
0094                                const RooCmdArg &arg2 = RooCmdArg::none(), const RooCmdArg &arg3 = RooCmdArg::none(),
0095                                const RooCmdArg &arg4 = RooCmdArg::none(), const RooCmdArg &arg5 = RooCmdArg::none(),
0096                                const RooCmdArg &arg6 = RooCmdArg::none(), const RooCmdArg &arg7 = RooCmdArg::none(),
0097                                const RooCmdArg &arg8 = RooCmdArg::none());
0098
0099    static std::shared_ptr<ROOT::Fit::FitConfig> createFitConfig(); // obtain instance of default fit configuration
0100    static std::shared_ptr<RooLinkedList> createNLLOptions();       // obtain instance of default nll options
0101    static std::shared_ptr<RooLinkedList> defaultNLLOptions();      // access default NLL options for modifications
0102    static std::shared_ptr<ROOT::Fit::FitConfig> defaultFitConfig();
0103    static ROOT::Math::IOptions *defaultFitConfigOptions();
0104
0105    static std::shared_ptr<const RooFitResult> minimize(RooAbsReal &nll,
0106                                                        const std::shared_ptr<ROOT::Fit::FitConfig> &fitConfig = nullptr,
0107                                                        const std::shared_ptr<RooLinkedList> &nllOpts = nullptr);
0108    static int minos(RooAbsReal &nll, const RooFitResult &ufit, const char *parName = "",
0109                     const std::shared_ptr<ROOT::Fit::FitConfig> &_fitConfig = nullptr);
0110
0111    // this class is used to store a shared_ptr in a TDirectory's List, so that retrieval of cached fits
0112    // can share the fit result (and avoid re-reading from disk as well)
0113    class StoredFitResult : public TNamed {
0114    public:
0115       StoredFitResult(RooFitResult *_fr);
0116       StoredFitResult(const std::shared_ptr<RooFitResult> &_fr);
0117
0118    public:
0119       std::shared_ptr<RooFitResult> fr; //!
0120       ClassDef(StoredFitResult, 0)
0121    };
0122
0123    // can't use an enum class because pyROOT doesn't handle it
0124    class TestStatistic {
0125    public:
0126       enum Type {
0127          // basic likelihood ratio
0128          tmu = -1,
0129          // upper limit test statistics
0130          qmu = -2,
0131          qmutilde = -3,
0132          // discovery test statistics
0133          q0 = -4,
0134          uncappedq0 = -5,
0135          u0 = -5
0136       };
0137    };
0138
0139    class Asymptotics {
0140
0141    public:
0142       typedef std::vector<std::pair<double, int>> IncompatFunc;
0143
0144       enum PLLType {
0145          TwoSided = 0,
0146          OneSidedPositive, // for exclusions
0147          OneSidedNegative, // for discovery
0148          OneSidedAbsolute, // for exclusions by magnitude
0149          Uncapped,         // for discovery with interest in deficits as well as excesses
0150          Unknown
0151       };
0152
0153       // The incompatibility function (taking mu_hat as an input) is defined by its transitions
0154       // it takes values of -1, 0, or 1 ... when it 0 that means mu_hat is compatible with the hypothesis
0155       // Standard incompatibility functions are parameterized by mu
0156       // Note: the default value is taken to be 1, so an empty vector is function=1
0157       static IncompatFunc IncompatibilityFunction(const PLLType &type, double mu)
0158       {
0159          std::vector<std::pair<double, int>> out;
0160          if (type == TwoSided) {
0161             // standard PLL
0162          } else if (type == OneSidedPositive) {
0163             out.emplace_back(std::make_pair(mu, 0)); // becomes compatible @ mu_hat = mu
0164          } else if (type == OneSidedNegative) {
0165             out.emplace_back(std::make_pair(-std::numeric_limits<double>::infinity(), 0)); // compatible at -inf
0166             out.emplace_back(std::make_pair(mu, 1)); // becomes incompatible at mu_hat = mu
0167          } else if (type == OneSidedAbsolute) {
0168             out.emplace_back(std::make_pair(-std::numeric_limits<double>::infinity(), 0)); // compatible at -inf
0169             out.emplace_back(std::make_pair(-mu, 1)); // incompatible @ mu_hat = -mu
0170             out.emplace_back(std::make_pair(mu, 0));  // compatible again @ mu_hat = mu
0171          } else if (type == Uncapped) {
0172             out.emplace_back(std::make_pair(-std::numeric_limits<double>::infinity(), -1)); // reversed at -inf
0173             out.emplace_back(std::make_pair(mu, 1)); // becomes normal @ mu_hat = mu
0174          } else {
0175             throw std::runtime_error("Unknown PLL Type");
0176          }
0177          return out;
0178       }
0179
0180       // inverse of PValue function
0181       static double k(const IncompatFunc &compatRegions, double pValue, double poiVal, double poiPrimeVal,
0182                       double sigma_mu = 0, double mu_low = -std::numeric_limits<double>::infinity(),
0183                       double mu_high = std::numeric_limits<double>::infinity());
0184
0185       static double k(const PLLType &pllType, double pValue, double mu, double mu_prime, double sigma_mu = 0,
0186                       double mu_low = -std::numeric_limits<double>::infinity(),
0187                       double mu_high = std::numeric_limits<double>::infinity())
0188       {
0189          return k(IncompatibilityFunction(pllType, mu), pValue, mu, mu_prime, sigma_mu, mu_low, mu_high);
0190       }
0191
0192       // Recommend sigma_mu = |mu - mu_prime|/sqrt(pll_mu(asimov_mu_prime))
0193       static double PValue(const IncompatFunc &compatRegions, double k, double mu, double mu_prime, double sigma_mu = 0,
0194                            double mu_low = -std::numeric_limits<double>::infinity(),
0195                            double mu_high = std::numeric_limits<double>::infinity());
0196
0197       static double PValue(const PLLType &pllType, double k, double mu, double mu_prime, double sigma_mu = 0,
0198                            double mu_low = -std::numeric_limits<double>::infinity(),
0199                            double mu_high = std::numeric_limits<double>::infinity())
0200       {
0201          return PValue(IncompatibilityFunction(pllType, mu), k, mu, mu_prime, sigma_mu, mu_low, mu_high);
0202       }
0203
0204       static double Phi_m(double mu, double mu_prime, double a, double sigma, const IncompatFunc &compatRegions);
0205
0206       static int CompatFactor(const IncompatFunc &func, double mu_hat);
0207
0208       static int CompatFactor(int type, double mu, double mu_hat)
0209       {
0210          return CompatFactor(IncompatibilityFunction((PLLType)type, mu), mu_hat);
0211       }
0212
0213       // converts pvalues to significances and finds where they equal the target pvalue
0214       // return is x-axis value with potentially an error on that value if input pVals had errors
0215       // static RooRealVar FindLimit(TGraph *pVals, double target_pVal = 0.05);
0216    };
0217
0218    static std::shared_ptr<RooLinkedList> sDefaultNLLOptions;
0219    static std::shared_ptr<ROOT::Fit::FitConfig> sDefaultFitConfig;
0220
0221    // Run hypothesis test(s) on the given pdf
0222    // Uses hypoPoint binning on model parameters to determine points to scan
0223    // if hypoPoint binning has nBins==0 then will auto-scan (assumes CL=95%, can override with setStringAttribute)
0224    // TODO: specifying number of null and alt toys per point
0225    static TCanvas *
0226    hypoTest(RooWorkspace &w, const xRooFit::Asymptotics::PLLType &pllType = xRooFit::Asymptotics::Unknown)
0227    {
0228       return hypoTest(w, 0, 0, pllType);
0229    }
0230    static TCanvas *hypoTest(RooWorkspace &w, int nToysNull, int nToysAlt,
0231                             const xRooFit::Asymptotics::PLLType &pllType = xRooFit::Asymptotics::Unknown);
0232 };
0233
0234 END_XROOFIT_NAMESPACE
0235
0236 #include "xRooHypoSpace.h"
0237 #include "xRooNLLVar.h"
0238 #include "xRooNode.h"
0239
0240 #endif // include guard