Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
HistFactoryImpl.cxx
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 * Authors:
4 * Jonas Rembser, CERN 2023
5 *
6 * Copyright (c) 2023, CERN
7 *
8 * Redistribution and use in source and binary forms,
9 * with or without modification, are permitted according to the terms
10 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
11 */
12
14
16
17#include <RooConstVar.h>
18#include <RooGaussian.h>
19#include <RooMsgService.h>
20#include <RooPoisson.h>
21#include <RooProduct.h>
22#include <RooRealVar.h>
23
24namespace RooStats {
25namespace HistFactory {
26namespace Detail {
27
28/**
29 * \brief Configure constrained gamma parameters for fitting.
30 *
31 * This function configures constrained gamma parameters for fitting. If a
32 * given relative sigma is less than or equal to zero or below a threshold, the
33 * gamma parameter is set to be constant. The function also sets reasonable
34 * ranges for the gamma parameter and provides a reasonable starting point for
35 * pre-fit errors.
36 *
37 * @param gammas The gamma parameters to be configured.
38 * @param sigmaRel The relative sigma values to be used for configuring the
39 * limits and errors.
40 * @param minSigma The minimum relative sigma threshold. If a relative sigma is
41 * below this threshold, the gamma parameter is set to be
42 * constant.
43 */
44void configureConstrainedGammas(RooArgList const &gammas, std::span<const double> relSigmas, double minSigma)
45{
46 assert(gammas.size() == relSigmas.size());
47
48 for (std::size_t i = 0; i < gammas.size(); ++i) {
49 auto &gamma = *static_cast<RooRealVar *>(gammas.at(i));
50 double sigmaRel = relSigmas[i];
51
52 // If the sigma is zero, the parameter might as well be constant
53 if (sigmaRel <= 0) {
54 gamma.setConstant(true);
55 continue;
56 }
57
58 // Set reasonable ranges
59 gamma.setMax(1. + 5. * sigmaRel);
60 gamma.setMin(0.);
61 // Set initial error too
62 gamma.setError(sigmaRel);
63
64 // Give reasonable starting point for pre-fit errors by setting it to the
65 // absolute sigma Mostly useful for pre-fit plotting.
66 // Note: in commit 2129c4d920 "[HF] Reduce verbosity of HistFactory."
67 // from 2020, there was a check added to do this only for Gaussian
68 // constrained parameters and for Poisson constrained parameters if they
69 // are stat errors without any justification. In the ROOT 6.30
70 // development cycle, this check got removed again to cause less surprise
71 // to the user.
72 gamma.setError(sigmaRel);
73
74 // If the sigma value is less than a supplied threshold, set the variable to
75 // constant
76 if (sigmaRel < minSigma) {
77 oocxcoutW(nullptr, HistFactory)
78 << "Warning: relative sigma " << sigmaRel << " for \"" << gamma.GetName() << "\" falls below threshold of "
79 << minSigma << ". Setting: " << gamma.GetName() << " to constant" << std::endl;
80 gamma.setConstant(true);
81 }
82 }
83}
84
85// Take a RooArgList of RooAbsReals and create N constraint terms (one for
86// each gamma) whose relative uncertainty is the value of the ith RooAbsReal
88 std::span<const double> relSigmas, double minSigma,
90{
92
93 // Check that there are N elements in the RooArgList
94 if (relSigmas.size() != paramSet.size()) {
95 std::cout << "Error: In createGammaConstraints, encountered bad number of relative sigmas" << std::endl;
96 std::cout << "Given vector with " << relSigmas.size() << " bins,"
97 << " but require exactly " << paramSet.size() << std::endl;
98 throw hf_exc();
99 }
100
101 configureConstrainedGammas(paramSet, relSigmas, minSigma);
102
103 for (std::size_t i = 0; i < paramSet.size(); ++i) {
104
105 RooRealVar &gamma = static_cast<RooRealVar &>(paramSet[i]);
106
107 oocxcoutI(nullptr, HistFactory)
108 << "Creating constraint for: " << gamma.GetName() << ". Type of constraint: " << type << std::endl;
109
110 const double sigmaRel = relSigmas[i];
111
112 // If the sigma is <= 0,
113 // do cont create the term
114 if (sigmaRel <= 0) {
115 oocxcoutI(nullptr, HistFactory)
116 << "Not creating constraint term for " << gamma.GetName() << " because sigma = " << sigmaRel
117 << " (sigma<=0)"
118 << " (bin number = " << i << ")" << std::endl;
119 continue;
120 }
121
122 // Make Constraint Term
123 std::string constrName = std::string(gamma.GetName()) + "_constraint";
124 std::string nomName = std::string("nom_") + gamma.GetName();
125
126 if (type == Constraint::Gaussian) {
127
128 // Type 1 : RooGaussian
129
130 // Make sigma
131 std::string sigmaName = std::string(gamma.GetName()) + "_sigma";
132 auto constrSigma = std::make_unique<RooConstVar>(sigmaName.c_str(), sigmaName.c_str(), sigmaRel);
133
134 // Make "observed" value
135 auto constrNom = std::make_unique<RooRealVar>(nomName.c_str(), nomName.c_str(), 1.0, 0, 10);
136 constrNom->setConstant(true);
137
138 // Make the constraint:
139 auto term = std::make_unique<RooGaussian>(constrName.c_str(), constrName.c_str(), *constrNom, gamma, *constrSigma);
140
141 out.globalObservables.push_back(constrNom.get());
142
143 term->addOwnedComponents(std::move(constrSigma));
144 term->addOwnedComponents(std::move(constrNom));
145
146 out.constraints.emplace_back(std::move(term));
147 } else if (type == Constraint::Poisson) {
148
149 // this is correct Poisson equivalent to a Gaussian with mean 1 and stdev sigma
150 const double tau = 1. / (sigmaRel * sigmaRel);
151
152 // Make nominal "observed" value
153 auto constrNom = std::make_unique<RooRealVar>(nomName.c_str(), nomName.c_str(), tau);
154 constrNom->setMin(0);
155 constrNom->setConstant(true);
156
157 // Make the scaling term
158 std::string scalingName = std::string(gamma.GetName()) + "_tau";
159 auto poissonScaling = std::make_unique<RooConstVar>(scalingName.c_str(), scalingName.c_str(), tau);
160
161 // Make mean for scaled Poisson
162 std::string poisMeanName = std::string(gamma.GetName()) + "_poisMean";
163 auto constrMean = std::make_unique<RooProduct>(poisMeanName.c_str(), poisMeanName.c_str(), gamma, *poissonScaling);
164
165 // Type 2 : RooPoisson
166 auto term = std::make_unique<RooPoisson>(constrName.c_str(), constrName.c_str(), *constrNom, *constrMean);
167 term->setNoRounding(true);
168
169 out.globalObservables.push_back(constrNom.get());
170
171 term->addOwnedComponents(std::move(poissonScaling));
172 term->addOwnedComponents(std::move(constrMean));
173 term->addOwnedComponents(std::move(constrNom));
174
175 out.constraints.emplace_back(std::move(term));
176 } else {
177
178 std::cout << "Error: Did not recognize Stat Error constraint term type: " << type
179 << " for : " << gamma.GetName() << std::endl;
180 throw hf_exc();
181 }
182 } // end loop over parameters
183
184 return out;
185}
186
187} // namespace Detail
188} // namespace HistFactory
189} // namespace RooStats
#define oocxcoutW(o, a)
#define oocxcoutI(o, a)
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
Storage_t::size_type size() const
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition RooArgList.h:110
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void configureConstrainedGammas(RooArgList const &gammas, std::span< const double > relSigmas, double minSigma)
Configure constrained gamma parameters for fitting.
CreateGammaConstraintsOutput createGammaConstraints(RooArgList const &paramList, std::span< const double > relSigmas, double minSigma, Constraint::Type type)
Namespace for the RooStats classes.
Definition CodegenImpl.h:58