root/RooStats/HypoTestInverterResult.h

0001 // @(#)root/roostats:$Id$
0002 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
0003 /*************************************************************************
0004  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
0005  * All rights reserved.                                                  *
0006  *                                                                       *
0007  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0008  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0009  *************************************************************************/
0010
0011 #ifndef ROOSTATS_HypoTestInverterResult
0012 #define ROOSTATS_HypoTestInverterResult
0013
0014 #include "RooStats/SimpleInterval.h"
0015
0016 #include "RooStats/HypoTestResult.h"
0017
0018 #include <vector>
0019
0020 class RooRealVar;
0021
0022 namespace RooStats {
0023
0024 class SamplingDistribution;
0025
0026 class HypoTestInverterResult : public SimpleInterval {
0027
0028 public:
0029
0030    /// default constructor
0031    explicit HypoTestInverterResult(const char *name = nullptr);
0032
0033    /// constructor
0034    HypoTestInverterResult( const char* name,
0035                            const RooRealVar& scannedVariable,
0036                            double cl ) ;
0037
0038    HypoTestInverterResult( const HypoTestInverterResult& other, const char* name );
0039
0040    /// destructor
0041    ~HypoTestInverterResult() override;
0042
0043    /// operator =
0044    HypoTestInverterResult& operator = (const HypoTestInverterResult& other);
0045
0046    /// remove points that appear to have failed.
0047    int ExclusionCleanup();
0048
0049    /// merge with the content of another HypoTestInverterResult object
0050    bool Add( const HypoTestInverterResult& otherResult );
0051
0052    ///add the result of a single point (an HypoTestRsult)
0053    bool Add( double x, const HypoTestResult & result );
0054
0055    /// function to return the value of the parameter of interest for the i^th entry in the results
0056    double GetXValue( int index ) const ;
0057
0058    /// function to return the value of the confidence level for the i^th entry in the results
0059    double GetYValue( int index ) const ;
0060
0061    /// function to return the estimated error on the value of the confidence level for the i^th entry in the results
0062    double GetYError( int index ) const ;
0063
0064    /// return the observed CLsplusb value  for the i-th entry
0065    double CLsplusb( int index) const;
0066
0067    /// return the observed CLb value  for the i-th entry
0068    double CLb( int index) const;
0069
0070    /// return the observed CLb value  for the i-th entry
0071    double CLs( int index) const;
0072
0073    /// return the observed CLsplusb value  for the i-th entry
0074    double CLsplusbError( int index) const;
0075
0076    /// return the observed CLb value  for the i-th entry
0077    double CLbError( int index) const;
0078
0079    /// return the observed CLb value  for the i-th entry
0080    double CLsError( int index) const;
0081
0082    /// return a pointer to the i^th result object
0083    HypoTestResult* GetResult( int index ) const ;
0084
0085    double GetLastYValue( ) const  { return GetYValue(  fXValues.size()-1); }
0086
0087    double GetLastXValue( ) const  { return GetXValue(  fXValues.size()-1); }
0088
0089    double GetLastYError( ) const  { return GetYError(  fXValues.size()-1); }
0090
0091    HypoTestResult * GetLastResult( ) const  { return GetResult(  fXValues.size()-1); }
0092
0093    /// number of entries in the results array
0094    int ArraySize() const { return fXValues.size(); };
0095
0096    int FindIndex(double xvalue) const;
0097
0098    /// set the size of the test (rate of Type I error) (eg. 0.05 for a 95% Confidence Interval)
0099    virtual void SetTestSize( double size ) { fConfidenceLevel = 1.-size; }
0100
0101    /// set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
0102    void SetConfidenceLevel( double cl ) override { fConfidenceLevel = cl; }
0103
0104    /// set CLs threshold for exclusion cleanup function
0105    inline void SetCLsCleanupThreshold( double th ) { fCLsCleanupThreshold = th; }
0106
0107    /// flag to switch between using CLsb (default) or CLs as confidence level
0108    void UseCLs( bool on = true ) { fUseCLs = on; }
0109
0110    /// query if one sided result
0111    bool IsOneSided() const { return !fIsTwoSided; }
0112    /// query if two sided result
0113    bool IsTwoSided() const { return fIsTwoSided; }
0114
0115    /// lower and upper bound of the confidence interval (to get upper/lower limits, multiply the size( = 1-confidence level ) by 2
0116    double LowerLimit() override;
0117    double UpperLimit() override;
0118
0119    /// rough estimation of the error on the computed bound of the confidence interval
0120    /// Estimate of lower limit error
0121    ///function evaluates only a rough error on the lower limit. Be careful when using this estimation
0122    double LowerLimitEstimatedError();
0123
0124    /// Estimate of lower limit error
0125    ///function evaluates only a rough error on the lower limit. Be careful when using this estimation
0126    double UpperLimitEstimatedError();
0127
0128    /// return expected distribution of p-values (Cls or Clsplusb)
0129
0130    SamplingDistribution * GetExpectedPValueDist(int index) const;
0131
0132    SamplingDistribution * GetBackgroundTestStatDist(int index ) const;
0133
0134    SamplingDistribution * GetSignalAndBackgroundTestStatDist(int index) const;
0135
0136    /// same in terms of alt and null
0137    SamplingDistribution * GetNullTestStatDist(int index) const {
0138       return  GetSignalAndBackgroundTestStatDist(index);
0139    }
0140    SamplingDistribution * GetAltTestStatDist(int index) const {
0141       return  GetBackgroundTestStatDist(index);
0142    }
0143
0144    /// get expected lower limit distributions
0145    /// implemented using interpolation
0146    ///  The size for the sampling distribution is given (by default is given by the average number of toy/point)
0147    SamplingDistribution* GetLowerLimitDistribution() const { return GetLimitDistribution(true); }
0148
0149    /// get expected upper limit distributions
0150    /// implemented using interpolation
0151    SamplingDistribution* GetUpperLimitDistribution() const { return GetLimitDistribution(false); }
0152
0153    /// get Limit value corresponding at the desired nsigma level (0) is median -1 sigma is 1 sigma
0154    double GetExpectedLowerLimit(double nsig = 0, const char * opt = "" ) const ;
0155
0156    /// get Limit value corresponding at the desired nsigma level (0) is median -1 sigma is 1 sigma
0157    double GetExpectedUpperLimit(double nsig = 0, const char * opt = "") const ;
0158
0159
0160    double FindInterpolatedLimit(double target, bool lowSearch = false, double xmin=1, double xmax=0.0);
0161
0162    enum InterpolOption_t { kLinear, kSpline };
0163
0164    /// set the interpolation option, linear (kLinear ) or spline  (kSpline)
0165    void SetInterpolationOption( InterpolOption_t opt) { fInterpolOption = opt; }
0166
0167    InterpolOption_t GetInterpolationOption() const { return fInterpolOption; }
0168
0169 private:
0170
0171
0172    double CalculateEstimatedError(double target, bool lower = true, double xmin = 1, double xmax = 0.0);
0173
0174    int FindClosestPointIndex(double target, int mode = 0, double xtarget = 0);
0175
0176    SamplingDistribution* GetLimitDistribution(bool lower ) const;
0177
0178    double GetExpectedLimit(double nsig, bool lower, const char * opt = "" ) const ;
0179
0180    double GetGraphX(const TGraph & g, double y0, bool lowSearch, double &xmin, double &xmax) const;
0181    double GetGraphX(const TGraph & g, double y0, bool lowSearch = true) const {
0182       double xmin=1; double xmax = 0;
0183       return GetGraphX(g,y0,lowSearch,xmin,xmax);
0184    }
0185
0186
0187 protected:
0188
0189    bool fUseCLs;
0190    bool fIsTwoSided;                  ///< two sided scan (look for lower/upper limit)
0191    bool fInterpolateLowerLimit;
0192    bool fInterpolateUpperLimit;
0193    bool fFittedLowerLimit;
0194    bool fFittedUpperLimit;
0195    InterpolOption_t fInterpolOption;  ///< interpolation option (linear or spline)
0196
0197    double fLowerLimitError;
0198    double fUpperLimitError;
0199
0200    double fCLsCleanupThreshold;
0201
0202    static double fgAsymptoticMaxSigma;  ///< max sigma value used to scan asymptotic expected p values
0203    static int fgAsymptoticNumPoints;    ///< number of points used to build expected p-values
0204
0205    std::vector<double> fXValues;
0206
0207    TList fYObjects;       ///< list of HypoTestResult for each point
0208    TList fExpPValues;     ///< list of expected sampling distribution for each point
0209
0210    friend class HypoTestInverter;
0211    friend class HypoTestInverterPlot;
0212
0213    ClassDefOverride(HypoTestInverterResult,5)  // HypoTestInverterResult class
0214 };
0215 }
0216
0217 #endif