Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooImproperIntegrator1D.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/**
18\file RooImproperIntegrator1D.cxx
19\class RooImproperIntegrator1D
20\ingroup Roofitcore
21
22Special numeric integrator that can handle integrals over open domains.
23To this end the range is cut in up three pieces: [-inf,-1],[-1,+1] and [+1,inf]
24and the outer two pieces, if required are calculated using a 1/x transform
25**/
26
29#include "RooInvTransform.h"
30#include "RooNumber.h"
31#include "RooNumIntFactory.h"
32#include "RooArgSet.h"
33#include "RooMsgService.h"
34
35#include "Riostream.h"
36#include <cmath>
37#include "TClass.h"
38
39
40
41// Register this class with RooNumIntConfig
42
43////////////////////////////////////////////////////////////////////////////////
44/// Register RooImproperIntegrator1D, its parameters and capabilities with RooNumIntFactory
45
47{
48 auto creator = [](const RooAbsFunc &function, const RooNumIntConfig &config) {
49 return std::make_unique<RooImproperIntegrator1D>(function, config);
50 };
51
52 fact.registerPlugin("RooImproperIntegrator1D", creator, {},
53 /*canIntegrate1D=*/true,
54 /*canIntegrate2D=*/false,
55 /*canIntegrateND=*/false,
56 /*canIntegrateOpenEnded=*/true,
57 /*depName=*/"RooIntegrator1D");
58}
59
60
61////////////////////////////////////////////////////////////////////////////////
62/// Constructor with function binding. The integration range is taken from the
63/// definition in the function binding
64
67 _useIntegrandLimits(true),
68 _origFunc(const_cast<RooAbsFunc*>(&function))
69{
71}
72
73
74
75////////////////////////////////////////////////////////////////////////////////
76/// Constructor with function binding and configuration object. The integration range is taken
77/// from the definition in the function binding
78
81 _useIntegrandLimits(true),
82 _origFunc(const_cast<RooAbsFunc*>(&function)),
83 _config(config)
84{
86}
87
88
89
90////////////////////////////////////////////////////////////////////////////////
91/// Constructor with function binding, definition of integration range and configuration object
92
95 _xmin(xmin),
96 _xmax(xmax),
97 _useIntegrandLimits(false),
98 _origFunc(const_cast<RooAbsFunc*>(&function)),
99 _config(config)
100{
102}
103
104
105
106////////////////////////////////////////////////////////////////////////////////
107/// Initialize the integrator, construct and initialize subintegrators
108
110{
111 if(!isValid()) {
112 oocoutE(nullptr,Integration) << "RooImproperIntegrator: cannot integrate invalid function" << std::endl;
113 return;
114 }
115 // Create a new function object that uses the change of vars: x -> 1/x
116 if (function) {
117 _function= std::make_unique<RooInvTransform>(*function);
118 } else {
120 _integrator1.reset();
121 _integrator2.reset();
122 _integrator3.reset();
123 }
124
125 // Helper function to create a new configuration that is just like the one
126 // associated to this integrator, but with a different summation rule.
127 auto makeIntegrator1D = [&](RooAbsFunc const& func,
128 double xmin, double xmax,
130 RooNumIntConfig newConfig{_config}; // copy default configuration
131 newConfig.getConfigSection("RooIntegrator1D").setCatIndex("sumRule", rule);
132 return std::make_unique<RooRombergIntegrator>(func, xmin, xmax, newConfig);
133 };
134
135 // partition the integration range into subranges that can each be
136 // handled by RooIntegrator1D
137 switch(_case= limitsCase()) {
138 case ClosedBothEnds:
139 // both limits are finite: use the plain trapezoid integrator
140 _integrator1 = std::make_unique<RooRombergIntegrator>(*function,_xmin,_xmax,_config);
141 break;
142 case OpenBothEnds:
143 // both limits are infinite: integrate over (-1,+1) using
144 // the plain trapezoid integrator...
145 _integrator1 = makeIntegrator1D(*function,-1,+1,RooRombergIntegrator::Trapezoid);
146 // ...and integrate the infinite tails using the midpoint integrator
147 _integrator2 = makeIntegrator1D(*_function,-1,0,RooRombergIntegrator::Midpoint);
148 _integrator3 = makeIntegrator1D(*_function,0,+1,RooRombergIntegrator::Midpoint);
149 break;
151 // xmax >= 0 so integrate from (-inf,-1) and (-1,xmax)
152 _integrator1 = makeIntegrator1D(*_function,-1,0,RooRombergIntegrator::Midpoint);
154 break;
155 case OpenBelow:
156 // xmax < 0 so integrate from (-inf,xmax)
158 break;
160 // xmin <= 0 so integrate from (xmin,+1) and (+1,+inf)
161 _integrator1 = makeIntegrator1D(*_function,0,+1,RooRombergIntegrator::Midpoint);
163 break;
164 case OpenAbove:
165 // xmin > 0 so integrate from (xmin,+inf)
167 break;
168 case Invalid:
169 default:
170 _valid= false;
171 }
172}
173
174
175////////////////////////////////////////////////////////////////////////////////
176/// Change our integration limits. Return true if the new limits are
177/// ok, or otherwise false. Always returns false and does nothing
178/// if this object was constructed to always use our integrand's limits.
179
181{
183 oocoutE(nullptr,Integration) << "RooImproperIntegrator1D::setLimits: cannot override integrand's limits" << std::endl;
184 return false;
185 }
186
187 _xmin= *xmin;
188 _xmax= *xmax;
189 return checkLimits();
190}
191
192
193////////////////////////////////////////////////////////////////////////////////
194/// Check if the limits are valid. For this integrator all limit configurations
195/// are valid, but if the limits change between two calculate() calls it
196/// may be necessary to reconfigure (e.g. if an open ended range becomes
197/// a closed range
198
200{
201 // Has either limit changed?
203 if(_xmin == integrand()->getMinLimit(0) &&
204 _xmax == integrand()->getMaxLimit(0)) return true;
205 }
206
207 // The limits have changed: can we use the same strategy?
208 if(limitsCase() != _case) {
209 // Reinitialize embedded integrators, will automatically propagate new limits
210 const_cast<RooImproperIntegrator1D*>(this)->initialize() ;
211 return true ;
212 }
213
214 // Reuse our existing integrators by updating their limits
215 switch(_case) {
216 case ClosedBothEnds:
217 _integrator1->setLimits(_xmin,_xmax);
218 break;
219 case OpenBothEnds:
220 // nothing has changed
221 break;
223 _integrator2->setLimits(-1,_xmax);
224 break;
225 case OpenBelow:
226 _integrator1->setLimits(1/_xmax,0);
227 break;
229 _integrator2->setLimits(_xmin,+1);
230 break;
231 case OpenAbove:
232 _integrator1->setLimits(0,1/_xmin);
233 break;
234 case Invalid:
235 default:
236 return false;
237 }
238 return true;
239}
240
241
242////////////////////////////////////////////////////////////////////////////////
243/// Classify the type of limits we have: OpenBothEnds,ClosedBothEnds,OpenBelow or OpenAbove.
244
246{
247 // Analyze the specified limits to determine which case applies.
248 if(nullptr == integrand() || !integrand()->isValid()) return Invalid;
249
253 }
254
255 bool inf1= RooNumber::isInfinite(_xmin);
256 bool inf2= RooNumber::isInfinite(_xmax);
257 if(!inf1 && !inf2) {
258 // both limits are finite
259 return ClosedBothEnds;
260 }
261 else if(inf1 && inf2) {
262 // both limits are infinite
263 return OpenBothEnds;
264 }
265 else if(inf1) { // inf2==false
266 if(_xmax >= 0) {
267 return OpenBelowSpansZero;
268 }
269 else {
270 return OpenBelow;
271 }
272 }
273 else { // inf1==false && inf2==true
274 if(_xmin <= 0) {
275 return OpenAboveSpansZero;
276 }
277 else {
278 return OpenAbove;
279 }
280 }
281 // return Invalid; OSF-CC: Statement unreachable
282}
283
284
285////////////////////////////////////////////////////////////////////////////////
286/// Calculate the integral at the given parameter values of the function binding
287
288double RooImproperIntegrator1D::integral(const double* yvec)
289{
290 double result(0);
291 if(_integrator1) result+= _integrator1->integral(yvec);
292 if(_integrator2) result+= _integrator2->integral(yvec);
293 if(_integrator3) result+= _integrator3->integral(yvec);
294 return result;
295}
RooAbsReal & function()
RooAbsReal * _origFunc
Original function.
bool initialize()
#define oocoutE(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 result
float xmin
float xmax
bool setCatIndex(const char *name, Int_t newVal=0, bool verbose=false)
Set index value of a RooAbsCategoryLValue stored in set with given name to newVal.
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition RooAbsFunc.h:27
virtual double getMaxLimit(UInt_t dimension) const =0
virtual double getMinLimit(UInt_t dimension) const =0
Abstract interface for integrators of real-valued functions that implement the RooAbsFunc interface.
bool isValid() const
Is integrator in valid state.
const RooAbsFunc * integrand() const
Return integrand function binding.
bool _valid
Is integrator in valid state?
Special numeric integrator that can handle integrals over open domains.
std::unique_ptr< RooRombergIntegrator > _integrator1
Piece integrator 1.
std::unique_ptr< RooInvTransform > _function
Binding with inverse of function.
LimitsCase limitsCase() const
Classify the type of limits we have: OpenBothEnds,ClosedBothEnds,OpenBelow or OpenAbove.
double integral(const double *yvec=nullptr) override
Calculate the integral at the given parameter values of the function binding.
RooImproperIntegrator1D(const RooAbsFunc &function)
Constructor with function binding.
LimitsCase _case
Configuration of limits.
double _xmax
Value of limits.
bool checkLimits() const override
Check if the limits are valid.
RooAbsFunc * _origFunc
Original function binding.
std::unique_ptr< RooRombergIntegrator > _integrator3
Piece integrator 3.
void initialize(const RooAbsFunc *function=nullptr)
Initialize the integrator, construct and initialize subintegrators.
bool setLimits(double *xmin, double *xmax) override
Change our integration limits.
std::unique_ptr< RooRombergIntegrator > _integrator2
Piece integrator 2.
RooNumIntConfig _config
Configuration object.
static void registerIntegrator(RooNumIntFactory &fact)
Register RooImproperIntegrator1D, its parameters and capabilities with RooNumIntFactory.
bool _useIntegrandLimits
Use limits in function binding?
Holds the configuration parameters of the various numeric integrators used by RooRealIntegral.
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
Factory to instantiate numeric integrators from a given function binding and a given configuration.
bool registerPlugin(std::string const &name, Creator const &creator, const RooArgSet &defConfig, bool canIntegrate1D, bool canIntegrate2D, bool canIntegrateND, bool canIntegrateOpenEnded, const char *depName="")
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...
static constexpr int isInfinite(double x)
Return true if x is infinite by RooNumber internal specification.
Definition RooNumber.h:27