RooStats/HistFactory/Systematics.h

0001 // @(#)root/roostats:$Id$
0002 // Author: George Lewis, Kyle Cranmer
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 HISTFACTORY_SYSTEMATICS_H
0012 #define HISTFACTORY_SYSTEMATICS_H
0013
0014 #include <string>
0015 #include <fstream>
0016 #include <iostream>
0017
0018 #include "TH1.h"
0019 #include "TDirectory.h"
0020 #include "RooStats/HistFactory/HistRef.h"
0021
0022 namespace RooStats{
0023 namespace HistFactory {
0024
0025   namespace Constraint {
0026     enum Type{ Gaussian, Poisson };
0027     std::string Name( Type type );
0028     Type GetType( const std::string& Name );
0029   }
0030
0031 /** \class OverallSys
0032  * \ingroup HistFactory
0033  * Configuration for a constrained overall systematic to scale sample normalisations.
0034  */
0035   class OverallSys {
0036
0037   public:
0038
0039     OverallSys() : fLow(0), fHigh(0) {}
0040
0041     void SetName( const std::string& Name ) { fName = Name; }
0042     const std::string& GetName() const { return fName; }
0043
0044     void SetLow( double Low )   { fLow  = Low; }
0045     void SetHigh( double High ) { fHigh = High; }
0046     double GetLow() const { return fLow; }
0047     double GetHigh() const { return fHigh; }
0048
0049     void Print(std::ostream& = std::cout) const;
0050     void PrintXML(std::ostream&) const;
0051
0052   protected:
0053     std::string fName;
0054     double fLow;
0055     double fHigh;
0056
0057   };
0058
0059 /** \class NormFactor
0060  * \ingroup HistFactory
0061  * Configuration for an \a un- constrained overall systematic to scale sample normalisations.
0062  */
0063   class NormFactor {
0064
0065   public:
0066
0067     NormFactor();
0068
0069     void SetName( const std::string& Name ) { fName = Name; }
0070     std::string GetName() const { return fName; }
0071
0072     void SetVal( double Val ) { fVal = Val; }
0073     double GetVal() const { return fVal; }
0074
0075     void SetLow( double Low )   { fLow  = Low; }
0076     void SetHigh( double High ) { fHigh = High; }
0077     double GetLow() const { return fLow; }
0078     double GetHigh() const { return fHigh; }
0079
0080     void Print(std::ostream& = std::cout) const;
0081     void PrintXML(std::ostream&) const;
0082
0083   protected:
0084
0085     std::string fName;
0086     double fVal;
0087     double fLow;
0088     double fHigh;
0089
0090   };
0091
0092
0093   /** ////////////////////////////////////////////////////////////////////////////////////////////
0094    * \class HistogramUncertaintyBase
0095    * \ingroup HistFactory
0096    * Base class to store the up and down variations for histogram uncertainties.
0097    * Use the derived classes for actual models.
0098    */
0099   class HistogramUncertaintyBase {
0100
0101   public:
0102
0103     HistogramUncertaintyBase() : fhLow(nullptr), fhHigh(nullptr) {}
0104     HistogramUncertaintyBase(const std::string& Name) : fName(Name), fhLow(nullptr), fhHigh(nullptr) {}
0105     HistogramUncertaintyBase(const HistogramUncertaintyBase& oth) :
0106     fName{oth.fName},
0107     fInputFileLow{oth.fInputFileLow}, fHistoNameLow{oth.fHistoNameLow}, fHistoPathLow{oth.fHistoPathLow},
0108     fInputFileHigh{oth.fInputFileHigh}, fHistoNameHigh{oth.fHistoNameHigh}, fHistoPathHigh{oth.fHistoPathHigh},
0109     fhLow{oth.fhLow ? static_cast<TH1*>(oth.fhLow->Clone()) : nullptr},
0110     fhHigh{oth.fhHigh ? static_cast<TH1*>(oth.fhHigh->Clone()) : nullptr} {
0111        if (fhLow)
0112           fhLow->SetDirectory(nullptr);
0113        if (fhHigh)
0114           fhHigh->SetDirectory(nullptr);
0115     }
0116     HistogramUncertaintyBase(HistogramUncertaintyBase&&) = default;
0117
0118     virtual ~HistogramUncertaintyBase() {};
0119
0120
0121     // Need deep copies because the class owns its histograms.
0122     HistogramUncertaintyBase& operator=(const HistogramUncertaintyBase& oth) {
0123       fName = oth.fName;
0124       fInputFileLow = oth.fInputFileLow;
0125       fHistoNameLow = oth.fHistoNameLow;
0126       fHistoPathLow = oth.fHistoPathLow;
0127       fInputFileHigh = oth.fInputFileHigh;
0128       fHistoNameHigh = oth.fHistoNameHigh;
0129       fHistoPathHigh = oth.fHistoPathHigh;
0130
0131       TDirectory::TContext ctx{nullptr}; // Don't associate clones to directories
0132       fhLow.reset(oth.fhLow ? static_cast<TH1*>(oth.fhLow->Clone()) : nullptr);
0133       fhHigh.reset(oth.fhHigh ? static_cast<TH1*>(oth.fhHigh->Clone()) : nullptr);
0134
0135       return *this;
0136     }
0137     HistogramUncertaintyBase& operator=(HistogramUncertaintyBase&&) = default;
0138
0139     virtual void Print(std::ostream& = std::cout) const;
0140     virtual void PrintXML(std::ostream&) const = 0;
0141     virtual void writeToFile( const std::string& FileName, const std::string& DirName );
0142
0143     void SetHistoLow(TH1* Low ) {Low->SetDirectory(nullptr); fhLow.reset(Low);}
0144     void SetHistoHigh(TH1* High ) {High->SetDirectory(nullptr); fhHigh.reset(High);}
0145
0146     const TH1* GetHistoLow() const {return fhLow.get();}
0147     const TH1* GetHistoHigh() const {return fhHigh.get();}
0148
0149     void SetName( const std::string& Name ) { fName = Name; }
0150     const std::string& GetName() const { return fName; }
0151
0152     void SetInputFileLow( const std::string& InputFileLow ) { fInputFileLow = InputFileLow; }
0153     void SetInputFileHigh( const std::string& InputFileHigh ) { fInputFileHigh = InputFileHigh; }
0154
0155     const std::string& GetInputFileLow() const { return fInputFileLow; }
0156     const std::string& GetInputFileHigh() const { return fInputFileHigh; }
0157
0158     void SetHistoNameLow( const std::string& HistoNameLow ) { fHistoNameLow = HistoNameLow; }
0159     void SetHistoNameHigh( const std::string& HistoNameHigh ) { fHistoNameHigh = HistoNameHigh; }
0160
0161     const std::string& GetHistoNameLow() const { return fHistoNameLow; }
0162     const std::string& GetHistoNameHigh() const { return fHistoNameHigh; }
0163
0164     void SetHistoPathLow( const std::string& HistoPathLow ) { fHistoPathLow = HistoPathLow; }
0165     void SetHistoPathHigh( const std::string& HistoPathHigh ) { fHistoPathHigh = HistoPathHigh; }
0166
0167     const std::string& GetHistoPathLow() const { return fHistoPathLow; }
0168     const std::string& GetHistoPathHigh() const { return fHistoPathHigh; }
0169
0170   protected:
0171
0172     std::string fName;
0173
0174     std::string fInputFileLow;
0175     std::string fHistoNameLow;
0176     std::string fHistoPathLow;
0177
0178     std::string fInputFileHigh;
0179     std::string fHistoNameHigh;
0180     std::string fHistoPathHigh;
0181
0182     // The Low and High Histograms
0183     std::unique_ptr<TH1> fhLow;
0184     std::unique_ptr<TH1> fhHigh;
0185
0186   };
0187
0188 /** \class HistoSys
0189  * \ingroup HistFactory
0190  * Configuration for a constrained, coherent shape variation of affected samples.
0191  */
0192 class HistoSys final : public HistogramUncertaintyBase {
0193 public:
0194   void PrintXML(std::ostream&) const override;
0195 };
0196
0197 /** \class HistoFactor
0198  * \ingroup HistFactory
0199  * Configuration for an *un*constrained, coherent shape variation of affected samples.
0200  */
0201   class HistoFactor final : public HistogramUncertaintyBase {
0202   public:
0203     void PrintXML(std::ostream&) const override;
0204   };
0205
0206 /** \class ShapeSys
0207  * \ingroup HistFactory
0208  * Constrained bin-by-bin variation of affected histogram.
0209  */
0210   class ShapeSys final : public HistogramUncertaintyBase {
0211
0212   public:
0213     ShapeSys() : fConstraintType(Constraint::Gaussian) {}
0214     ShapeSys(const ShapeSys& other) :
0215       HistogramUncertaintyBase(other),
0216       fConstraintType(other.fConstraintType) {}
0217     ShapeSys& operator=(const ShapeSys& oth) {
0218        if (this == &oth) return *this;
0219        HistogramUncertaintyBase::operator=(oth);
0220        fConstraintType = oth.fConstraintType;
0221        return *this;
0222     }
0223     ShapeSys& operator=(ShapeSys&&) = default;
0224
0225     void SetInputFile( const std::string& InputFile ) { fInputFileHigh = InputFile; }
0226     std::string GetInputFile() const { return fInputFileHigh; }
0227
0228     void SetHistoName( const std::string& HistoName ) { fHistoNameHigh = HistoName; }
0229     std::string GetHistoName() const { return fHistoNameHigh; }
0230
0231     void SetHistoPath( const std::string& HistoPath ) { fHistoPathHigh = HistoPath; }
0232     std::string GetHistoPath() const { return fHistoPathHigh; }
0233
0234     void Print(std::ostream& = std::cout) const override;
0235     void PrintXML(std::ostream&) const override;
0236     void writeToFile( const std::string& FileName, const std::string& DirName ) override;
0237
0238     const TH1* GetErrorHist() const {
0239       return fhHigh.get();
0240     }
0241     void SetErrorHist(TH1* hError) {
0242       fhHigh.reset(hError);
0243     }
0244
0245     void SetConstraintType( Constraint::Type ConstrType ) { fConstraintType = ConstrType; }
0246     Constraint::Type GetConstraintType() const { return fConstraintType; }
0247
0248   protected:
0249     Constraint::Type fConstraintType;
0250   };
0251
0252 /** \class ShapeFactor
0253  * \ingroup HistFactory
0254  * *Un*constrained bin-by-bin variation of affected histogram.
0255  */
0256   class ShapeFactor : public HistogramUncertaintyBase {
0257
0258   public:
0259     void Print(std::ostream& = std::cout) const override;
0260     void PrintXML(std::ostream&) const override;
0261     void writeToFile( const std::string& FileName, const std::string& DirName) override;
0262
0263     void SetInitialShape(TH1* shape) {
0264       fhHigh.reset(shape);
0265     }
0266     const TH1* GetInitialShape() const { return fhHigh.get(); }
0267
0268     void SetConstant(bool constant) { fConstant = constant; }
0269     bool IsConstant() const { return fConstant; }
0270
0271     bool HasInitialShape() const { return fHasInitialShape; }
0272
0273     void SetInputFile( const std::string& InputFile ) {
0274       fInputFileHigh = InputFile;
0275       fHasInitialShape=true;
0276     }
0277     const std::string& GetInputFile() const { return fInputFileHigh; }
0278
0279     void SetHistoName( const std::string& HistoName ) {
0280       fHistoNameHigh = HistoName;
0281       fHasInitialShape=true;
0282     }
0283     const std::string& GetHistoName() const { return fHistoNameHigh; }
0284
0285     void SetHistoPath( const std::string& HistoPath ) {
0286       fHistoPathHigh = HistoPath;
0287       fHasInitialShape=true;
0288     }
0289     const std::string& GetHistoPath() const { return fHistoPathHigh; }
0290
0291   protected:
0292
0293     bool fConstant = false;
0294
0295     // A histogram representing
0296     // the initial shape
0297     bool fHasInitialShape = false;
0298   };
0299
0300 /** \class StatError
0301  * \ingroup HistFactory
0302  * Statistical error of Monte Carlo predictions.
0303  */
0304   class StatError : public HistogramUncertaintyBase {
0305
0306   public:
0307     void Print(std::ostream& = std::cout) const override;
0308     void PrintXML(std::ostream&) const override;
0309     void writeToFile( const std::string& FileName, const std::string& DirName ) override;
0310
0311     void Activate( bool IsActive=true ) { fActivate = IsActive; }
0312     bool GetActivate() const { return fActivate; }
0313
0314     void SetUseHisto( bool UseHisto=true ) { fUseHisto = UseHisto; }
0315     bool GetUseHisto() const { return fUseHisto; }
0316
0317     void SetInputFile( const std::string& InputFile ) { fInputFileHigh = InputFile; }
0318     const std::string& GetInputFile() const { return fInputFileHigh; }
0319
0320     void SetHistoName( const std::string& HistoName ) { fHistoNameHigh = HistoName; }
0321     const std::string& GetHistoName() const { return fHistoNameHigh; }
0322
0323     void SetHistoPath( const std::string& HistoPath ) { fHistoPathHigh = HistoPath; }
0324     const std::string& GetHistoPath() const { return fHistoPathHigh; }
0325
0326
0327     const TH1* GetErrorHist() const {
0328       return fhHigh.get();
0329     }
0330     void SetErrorHist(TH1* Error) {
0331       fhHigh.reset(Error);
0332     }
0333
0334   protected:
0335
0336     bool fActivate = false;
0337     bool fUseHisto = false; // Use an external histogram for the errors
0338   };
0339
0340 /** \class StatErrorConfig
0341  * \ingroup HistFactory
0342  * Configuration to automatically assign nuisance parameters for the statistical
0343  * error of the Monte Carlo simulations.
0344  * The default is to assign a Poisson uncertainty to a bin when its statistical uncertainty
0345  * is larger than 5% of the bin content.
0346  */
0347   class StatErrorConfig {
0348
0349   public:
0350
0351     StatErrorConfig() : fRelErrorThreshold( .05 ), fConstraintType( Constraint::Poisson ) {;}
0352     void Print(std::ostream& = std::cout) const;
0353     void PrintXML(std::ostream&) const;
0354
0355     void SetRelErrorThreshold( double Threshold ) { fRelErrorThreshold = Threshold; }
0356     double GetRelErrorThreshold() const { return fRelErrorThreshold; }
0357
0358     void SetConstraintType( Constraint::Type ConstrType ) { fConstraintType = ConstrType; }
0359     Constraint::Type GetConstraintType() const { return fConstraintType; }
0360
0361   protected:
0362
0363     double fRelErrorThreshold;
0364     Constraint::Type fConstraintType;
0365
0366   };
0367
0368
0369 }
0370 }
0371
0372 #endif