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 RooCmdArg MaxIterations(int nIterations);
0070
0071    static constexpr double OBS = std::numeric_limits<double>::quiet_NaN();
0072
0073    // Helper function for matching precision of a value and its error
0074    static std::pair<double, double> matchPrecision(const std::pair<double, double> &in);
0075
0076    // Static methods that work with the 'first class' object types:
0077    //    Pdfs: RooAbsPdf
0078    //    Datasets: std::pair<RooAbsData,const RooAbsCollection>
0079    //    NLLOptions: RooLinkedList
0080    //    FitOptions: ROOT::Fit::FitConfig
0081
0082    // fit result flags in its constPars list which are global observables with the "global" attribute
0083    static std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>>
0084    generateFrom(RooAbsPdf &pdf, const RooFitResult &fr, bool expected = false, int seed = 0);
0085    static std::shared_ptr<const RooFitResult>
0086    fitTo(RooAbsPdf &pdf, const std::pair<std::shared_ptr<RooAbsData>, std::shared_ptr<const RooAbsCollection>> &data,
0087          const RooLinkedList &nllOpts, const ROOT::Fit::FitConfig &fitConf);
0088    static std::shared_ptr<const RooFitResult> fitTo(RooAbsPdf &pdf,
0089                                                     const std::pair<RooAbsData *, const RooAbsCollection *> &data,
0090                                                     const RooLinkedList &nllOpts, const ROOT::Fit::FitConfig &fitConf);
0091
0092    static xRooNLLVar createNLL(const std::shared_ptr<RooAbsPdf> pdf, const std::shared_ptr<RooAbsData> data,
0093                                const RooLinkedList &nllOpts);
0094    static xRooNLLVar createNLL(RooAbsPdf &pdf, RooAbsData *data, const RooLinkedList &nllOpts);
0095    static xRooNLLVar createNLL(RooAbsPdf &pdf, RooAbsData *data, const RooCmdArg &arg1 = RooCmdArg::none(),
0096                                const RooCmdArg &arg2 = RooCmdArg::none(), const RooCmdArg &arg3 = RooCmdArg::none(),
0097                                const RooCmdArg &arg4 = RooCmdArg::none(), const RooCmdArg &arg5 = RooCmdArg::none(),
0098                                const RooCmdArg &arg6 = RooCmdArg::none(), const RooCmdArg &arg7 = RooCmdArg::none(),
0099                                const RooCmdArg &arg8 = RooCmdArg::none());
0100
0101    static std::shared_ptr<ROOT::Fit::FitConfig> createFitConfig(); // obtain instance of default fit configuration
0102    static std::shared_ptr<RooLinkedList> createNLLOptions();       // obtain instance of default nll options
0103    static std::shared_ptr<RooLinkedList> defaultNLLOptions();      // access default NLL options for modifications
0104    static std::shared_ptr<ROOT::Fit::FitConfig> defaultFitConfig();
0105    static ROOT::Math::IOptions *defaultFitConfigOptions();
0106
0107    static std::shared_ptr<const RooFitResult> minimize(RooAbsReal &nll,
0108                                                        const std::shared_ptr<ROOT::Fit::FitConfig> &fitConfig = nullptr,
0109                                                        const std::shared_ptr<RooLinkedList> &nllOpts = nullptr);
0110    static int minos(RooAbsReal &nll, const RooFitResult &ufit, const char *parName = "",
0111                     const std::shared_ptr<ROOT::Fit::FitConfig> &_fitConfig = nullptr);
0112
0113    // this class is used to store a shared_ptr in a TDirectory's List, so that retrieval of cached fits
0114    // can share the fit result (and avoid re-reading from disk as well)
0115    class StoredFitResult : public TNamed {
0116    public:
0117       StoredFitResult(RooFitResult *_fr);
0118       StoredFitResult(const std::shared_ptr<RooFitResult> &_fr);
0119
0120    public:
0121       std::shared_ptr<RooFitResult> fr; //!
0122       ClassDef(StoredFitResult, 0)
0123    };
0124
0125    // can't use an enum class because pyROOT doesn't handle it
0126    class TestStatistic {
0127    public:
0128       enum Type {
0129          // basic likelihood ratio
0130          tmu = -1,
0131          // upper limit test statistics
0132          qmu = -2,
0133          qmutilde = -3,
0134          // discovery test statistics
0135          q0 = -4,
0136          uncappedq0 = -5,
0137          u0 = -5
0138       };
0139    };
0140
0141    class Asymptotics {
0142
0143    public:
0144       typedef std::vector<std::pair<double, int>> IncompatFunc;
0145
0146       enum PLLType {
0147          TwoSided = 0,
0148          OneSidedPositive, // for exclusions
0149          OneSidedNegative, // for discovery
0150          OneSidedAbsolute, // for exclusions by magnitude
0151          Uncapped,         // for discovery with interest in deficits as well as excesses
0152          Unknown
0153       };
0154
0155       // The incompatibility function (taking mu_hat as an input) is defined by its transitions
0156       // it takes values of -1, 0, or 1 ... when it 0 that means mu_hat is compatible with the hypothesis
0157       // Standard incompatibility functions are parameterized by mu
0158       // Note: the default value is taken to be 1, so an empty vector is function=1
0159       static IncompatFunc IncompatibilityFunction(const PLLType &type, double mu)
0160       {
0161          std::vector<std::pair<double, int>> out;
0162          if (type == TwoSided) {
0163             // standard PLL
0164          } else if (type == OneSidedPositive) {
0165             out.emplace_back(std::make_pair(mu, 0)); // becomes compatible @ mu_hat = mu
0166          } else if (type == OneSidedNegative) {
0167             out.emplace_back(std::make_pair(-std::numeric_limits<double>::infinity(), 0)); // compatible at -inf
0168             out.emplace_back(std::make_pair(mu, 1)); // becomes incompatible at mu_hat = mu
0169          } else if (type == OneSidedAbsolute) {
0170             out.emplace_back(std::make_pair(-std::numeric_limits<double>::infinity(), 0)); // compatible at -inf
0171             out.emplace_back(std::make_pair(-mu, 1)); // incompatible @ mu_hat = -mu
0172             out.emplace_back(std::make_pair(mu, 0));  // compatible again @ mu_hat = mu
0173          } else if (type == Uncapped) {
0174             out.emplace_back(std::make_pair(-std::numeric_limits<double>::infinity(), -1)); // reversed at -inf
0175             out.emplace_back(std::make_pair(mu, 1)); // becomes normal @ mu_hat = mu
0176          } else {
0177             throw std::runtime_error("Unknown PLL Type");
0178          }
0179          return out;
0180       }
0181
0182       // inverse of PValue function
0183       static double k(const IncompatFunc &compatRegions, double pValue, double poiVal, double poiPrimeVal,
0184                       double sigma_mu = 0, double mu_low = -std::numeric_limits<double>::infinity(),
0185                       double mu_high = std::numeric_limits<double>::infinity());
0186
0187       static double k(const PLLType &pllType, double pValue, double mu, double mu_prime, double sigma_mu = 0,
0188                       double mu_low = -std::numeric_limits<double>::infinity(),
0189                       double mu_high = std::numeric_limits<double>::infinity())
0190       {
0191          return k(IncompatibilityFunction(pllType, mu), pValue, mu, mu_prime, sigma_mu, mu_low, mu_high);
0192       }
0193
0194       // Recommend sigma_mu = |mu - mu_prime|/sqrt(pll_mu(asimov_mu_prime))
0195       static double PValue(const IncompatFunc &compatRegions, double k, double mu, double mu_prime, double sigma_mu = 0,
0196                            double mu_low = -std::numeric_limits<double>::infinity(),
0197                            double mu_high = std::numeric_limits<double>::infinity());
0198
0199       static double PValue(const PLLType &pllType, double k, double mu, double mu_prime, double sigma_mu = 0,
0200                            double mu_low = -std::numeric_limits<double>::infinity(),
0201                            double mu_high = std::numeric_limits<double>::infinity())
0202       {
0203          return PValue(IncompatibilityFunction(pllType, mu), k, mu, mu_prime, sigma_mu, mu_low, mu_high);
0204       }
0205
0206       static double Phi_m(double mu, double mu_prime, double a, double sigma, const IncompatFunc &compatRegions);
0207
0208       static int CompatFactor(const IncompatFunc &func, double mu_hat);
0209
0210       static int CompatFactor(int type, double mu, double mu_hat)
0211       {
0212          return CompatFactor(IncompatibilityFunction((PLLType)type, mu), mu_hat);
0213       }
0214
0215       // converts pvalues to significances and finds where they equal the target pvalue
0216       // return is x-axis value with potentially an error on that value if input pVals had errors
0217       // static RooRealVar FindLimit(TGraph *pVals, double target_pVal = 0.05);
0218    };
0219
0220    static std::shared_ptr<RooLinkedList> sDefaultNLLOptions;
0221    static std::shared_ptr<ROOT::Fit::FitConfig> sDefaultFitConfig;
0222
0223    // Run hypothesis test(s) on the given pdf
0224    // Uses hypoPoint binning on model parameters to determine points to scan
0225    // if hypoPoint binning has nBins==0 then will auto-scan (assumes CL=95%, can override with setStringAttribute)
0226    // TODO: specifying number of null and alt toys per point
0227    static TCanvas *
0228    hypoTest(RooWorkspace &w, const xRooFit::Asymptotics::PLLType &pllType = xRooFit::Asymptotics::Unknown)
0229    {
0230       return hypoTest(w, 0, 0, pllType);
0231    }
0232    static TCanvas *hypoTest(RooWorkspace &w, int nToysNull, int nToysAlt,
0233                             const xRooFit::Asymptotics::PLLType &pllType = xRooFit::Asymptotics::Unknown);
0234 };
0235
0236 END_XROOFIT_NAMESPACE
0237
0238 #include "xRooHypoSpace.h"
0239 #include "xRooNLLVar.h"
0240 #include "xRooNode.h"
0241
0242 #endif // include guard