Logo ROOT  
Reference Guide
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
28#include "RooIntegrator1D.h"
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 <math.h>
37#include "TClass.h"
38
39
41
42// Register this class with RooNumIntConfig
43
44////////////////////////////////////////////////////////////////////////////////
45/// Register RooImproperIntegrator1D, its parameters and capabilities with RooNumIntFactory
46
48{
51}
52
53
54
55////////////////////////////////////////////////////////////////////////////////
56/// Default constructor
57
59 _case(ClosedBothEnds), _xmin(-10), _xmax(10), _useIntegrandLimits(true)
60{
61}
62
63
64////////////////////////////////////////////////////////////////////////////////
65/// Constructor with function binding. The integration range is taken from the
66/// definition in the function binding
67
70 _useIntegrandLimits(true),
71 _origFunc((RooAbsFunc*)&function)
72{
74}
75
76
77
78////////////////////////////////////////////////////////////////////////////////
79/// Constructor with function binding and configuration object. The integration range is taken
80/// from the definition in the function binding
81
84 _useIntegrandLimits(true),
85 _origFunc((RooAbsFunc*)&function),
86 _config(config)
87{
89}
90
91
92
93////////////////////////////////////////////////////////////////////////////////
94/// Constructor with function binding, definition of integration range and configuration object
95
98 _xmin(xmin),
99 _xmax(xmax),
100 _useIntegrandLimits(false),
101 _origFunc((RooAbsFunc*)&function),
102 _config(config)
103{
105}
106
107
108
109////////////////////////////////////////////////////////////////////////////////
110/// Return clone of integrator with given function and configuration. Needed by RooNumIntFactory.
111
113{
114 return new RooImproperIntegrator1D(function,config) ;
115}
116
117
118
119////////////////////////////////////////////////////////////////////////////////
120/// Initialize the integrator, construct and initialize subintegrators
121
123{
124 if(!isValid()) {
125 oocoutE(nullptr,Integration) << "RooImproperIntegrator: cannot integrate invalid function" << std::endl;
126 return;
127 }
128 // Create a new function object that uses the change of vars: x -> 1/x
129 if (function) {
130 _function= std::make_unique<RooInvTransform>(*function);
131 } else {
133 _integrator1.reset();
134 _integrator2.reset();
135 _integrator3.reset();
136 }
137
138 // Helper function to create a new configuration that is just like the one
139 // associated to this integrator, but with a different summation rule.
140 auto makeIntegrator1D = [&](RooAbsFunc const& func,
141 double xmin, double xmax,
143 RooNumIntConfig newConfig{_config}; // copy default configuration
144 newConfig.getConfigSection("RooIntegrator1D").setCatIndex("sumRule", rule);
145 return std::make_unique<RooIntegrator1D>(func, xmin, xmax, newConfig);
146 };
147
148 // partition the integration range into subranges that can each be
149 // handled by RooIntegrator1D
150 switch(_case= limitsCase()) {
151 case ClosedBothEnds:
152 // both limits are finite: use the plain trapezoid integrator
153 _integrator1 = std::make_unique<RooIntegrator1D>(*function,_xmin,_xmax,_config);
154 break;
155 case OpenBothEnds:
156 // both limits are infinite: integrate over (-1,+1) using
157 // the plain trapezoid integrator...
158 _integrator1 = makeIntegrator1D(*function,-1,+1,RooIntegrator1D::Trapezoid);
159 // ...and integrate the infinite tails using the midpoint integrator
160 _integrator2 = makeIntegrator1D(*_function,-1,0,RooIntegrator1D::Midpoint);
161 _integrator3 = makeIntegrator1D(*_function,0,+1,RooIntegrator1D::Midpoint);
162 break;
164 // xmax >= 0 so integrate from (-inf,-1) and (-1,xmax)
165 _integrator1 = makeIntegrator1D(*_function,-1,0,RooIntegrator1D::Midpoint);
167 break;
168 case OpenBelow:
169 // xmax < 0 so integrate from (-inf,xmax)
170 _integrator1 = makeIntegrator1D(*_function,1/_xmax,0,RooIntegrator1D::Midpoint);
171 break;
173 // xmin <= 0 so integrate from (xmin,+1) and (+1,+inf)
174 _integrator1 = makeIntegrator1D(*_function,0,+1,RooIntegrator1D::Midpoint);
176 break;
177 case OpenAbove:
178 // xmin > 0 so integrate from (xmin,+inf)
179 _integrator1 = makeIntegrator1D(*_function,0,1/_xmin,RooIntegrator1D::Midpoint);
180 break;
181 case Invalid:
182 default:
183 _valid= false;
184 }
185}
186
187
188////////////////////////////////////////////////////////////////////////////////
189/// Change our integration limits. Return true if the new limits are
190/// ok, or otherwise false. Always returns false and does nothing
191/// if this object was constructed to always use our integrand's limits.
192
194{
196 oocoutE(nullptr,Integration) << "RooIntegrator1D::setLimits: cannot override integrand's limits" << std::endl;
197 return false;
198 }
199
200 _xmin= *xmin;
201 _xmax= *xmax;
202 return checkLimits();
203}
204
205
206////////////////////////////////////////////////////////////////////////////////
207/// Check if the limits are valid. For this integrator all limit configurations
208/// are valid, but if the limits change between two calculate() calls it
209/// may be necessary to reconfigure (e.g. if an open ended range becomes
210/// a closed range
211
213{
214 // Has either limit changed?
216 if(_xmin == integrand()->getMinLimit(0) &&
217 _xmax == integrand()->getMaxLimit(0)) return true;
218 }
219
220 // The limits have changed: can we use the same strategy?
221 if(limitsCase() != _case) {
222 // Reinitialize embedded integrators, will automatically propagate new limits
223 const_cast<RooImproperIntegrator1D*>(this)->initialize() ;
224 return true ;
225 }
226
227 // Reuse our existing integrators by updating their limits
228 switch(_case) {
229 case ClosedBothEnds:
230 _integrator1->setLimits(_xmin,_xmax);
231 break;
232 case OpenBothEnds:
233 // nothing has changed
234 break;
236 _integrator2->setLimits(-1,_xmax);
237 break;
238 case OpenBelow:
239 _integrator1->setLimits(1/_xmax,0);
240 break;
242 _integrator2->setLimits(_xmin,+1);
243 break;
244 case OpenAbove:
245 _integrator1->setLimits(0,1/_xmin);
246 break;
247 case Invalid:
248 default:
249 return false;
250 }
251 return true;
252}
253
254
255////////////////////////////////////////////////////////////////////////////////
256/// Classify the type of limits we have: OpenBothEnds,ClosedBothEnds,OpenBelow or OpenAbove.
257
259{
260 // Analyze the specified limits to determine which case applies.
261 if(0 == integrand() || !integrand()->isValid()) return Invalid;
262
266 }
267
268 bool inf1= RooNumber::isInfinite(_xmin);
269 bool inf2= RooNumber::isInfinite(_xmax);
270 if(!inf1 && !inf2) {
271 // both limits are finite
272 return ClosedBothEnds;
273 }
274 else if(inf1 && inf2) {
275 // both limits are infinite
276 return OpenBothEnds;
277 }
278 else if(inf1) { // inf2==false
279 if(_xmax >= 0) {
280 return OpenBelowSpansZero;
281 }
282 else {
283 return OpenBelow;
284 }
285 }
286 else { // inf1==false && inf2==true
287 if(_xmin <= 0) {
288 return OpenAboveSpansZero;
289 }
290 else {
291 return OpenAbove;
292 }
293 }
294 // return Invalid; OSF-CC: Statement unreachable
295}
296
297
298////////////////////////////////////////////////////////////////////////////////
299/// Calculate the integral at the given parameter values of the function binding
300
301double RooImproperIntegrator1D::integral(const double* yvec)
302{
303 double result(0);
304 if(_integrator1) result+= _integrator1->integral(yvec);
305 if(_integrator2) result+= _integrator2->integral(yvec);
306 if(_integrator3) result+= _integrator3->integral(yvec);
307 return result;
308}
#define oocoutE(o, a)
Definition: RooMsgService.h:52
#define ClassImp(name)
Definition: Rtypes.h:375
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
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
const char * proto
Definition: civetweb.c:17502
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
RooAbsIntegrator is the abstract interface for integrators of real-valued functions that implement th...
bool isValid() const
Is integrator in valid state.
const RooAbsFunc * integrand() const
Return integrand function binding.
bool _valid
Is integrator in valid state?
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
Special numeric integrator that can handle integrals over open domains.
std::unique_ptr< RooInvTransform > _function
Binding with inverse of function.
std::unique_ptr< RooIntegrator1D > _integrator2
Piece integrator 2.
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()
Default constructor.
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< RooIntegrator1D > _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.
RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const override
Return clone of integrator with given function and configuration. Needed by RooNumIntFactory.
std::unique_ptr< RooIntegrator1D > _integrator1
Piece integrator 1.
RooNumIntConfig _config
Configuration object.
static void registerIntegrator(RooNumIntFactory &fact)
Register RooImproperIntegrator1D, its parameters and capabilities with RooNumIntFactory.
bool _useIntegrandLimits
Use limits in function binding?
static TClass * Class()
RooNumIntConfig holds the configuration parameters of the various numeric integrators used by RooReal...
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
RooNumIntFactory is a factory to instantiate numeric integrators from a given function binding and a ...
bool storeProtoIntegrator(RooAbsIntegrator *proto, const RooArgSet &defConfig, const char *depName="")
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...
static Int_t isInfinite(double x)
Return true if x is infinite by RooNumBer internal specification.
Definition: RooNumber.cxx:57
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:433
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:167
@ Integration
Definition: RooGlobalFunc.h:61