Logo ROOT   6.10/09
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 
22 Special numeric integrator that can handle integrals over open domains.
23 To this end the range is cut in up three pieces: [-inf,-1],[-1,+1] and [+1,inf]
24 and the outer two pieces, if required are calculated using a 1/x transform
25 **/
26 
27 
28 #include "RooFit.h"
29 
32 #include "RooIntegrator1D.h"
33 #include "RooInvTransform.h"
34 #include "RooNumber.h"
35 #include "RooNumIntFactory.h"
36 #include "RooArgSet.h"
37 #include "RooMsgService.h"
38 
39 #include "Riostream.h"
40 #include <math.h>
41 #include "TClass.h"
42 
43 
44 
45 using namespace std;
46 
48 ;
49 
50 // Register this class with RooNumIntConfig
51 
52 ////////////////////////////////////////////////////////////////////////////////
53 /// Register RooImproperIntegrator1D, its parameters and capabilities with RooNumIntFactory
54 
56 {
59 }
60 
61 
62 
63 ////////////////////////////////////////////////////////////////////////////////
64 /// Default constructor
65 
67  _case(ClosedBothEnds), _xmin(-10), _xmax(10), _useIntegrandLimits(kTRUE),
68  _origFunc(0), _function(0), _integrator1(0), _integrator2(0), _integrator3(0)
69 {
70 }
71 
72 
73 ////////////////////////////////////////////////////////////////////////////////
74 /// Constructor with function binding. The integration range is taken from the
75 /// definition in the function binding
76 
81  _function(0),
82  _integrator1(0),
83  _integrator2(0),
84  _integrator3(0)
85 {
86  initialize(&function) ;
87 }
88 
89 
90 
91 ////////////////////////////////////////////////////////////////////////////////
92 /// Constructor with function binding and configuration object. The integration range is taken
93 /// from the definition in the function binding
94 
99  _function(0),
100  _config(config),
101  _integrator1(0),
102  _integrator2(0),
103  _integrator3(0)
104 {
105  initialize(&function) ;
106 }
107 
108 
109 
110 ////////////////////////////////////////////////////////////////////////////////
111 /// Constructor with function binding, definition of integration range and configuration object
112 
115  _xmin(xmin),
116  _xmax(xmax),
119  _function(0),
120  _config(config),
121  _integrator1(0),
122  _integrator2(0),
123  _integrator3(0)
124 {
125  initialize(&function) ;
126 }
127 
128 
129 
130 ////////////////////////////////////////////////////////////////////////////////
131 /// Return clone of integrator with given function and configuration. Needed by RooNumIntFactory.
132 
134 {
135  return new RooImproperIntegrator1D(function,config) ;
136 }
137 
138 
139 
140 ////////////////////////////////////////////////////////////////////////////////
141 /// Initialize the integrator, construct and initialize subintegrators
142 
144 {
145  if(!isValid()) {
146  oocoutE((TObject*)0,Integration) << "RooImproperIntegrator: cannot integrate invalid function" << endl;
147  return;
148  }
149  // Create a new function object that uses the change of vars: x -> 1/x
150  if (function) {
151  _function= new RooInvTransform(*function);
152  } else {
153  function = _origFunc ;
154  if (_integrator1) {
155  delete _integrator1 ;
156  _integrator1 = 0 ;
157  }
158  if (_integrator2) {
159  delete _integrator2 ;
160  _integrator2 = 0 ;
161  }
162  if (_integrator3) {
163  delete _integrator3 ;
164  _integrator3 = 0 ;
165  }
166  }
167 
168  // partition the integration range into subranges that can each be
169  // handled by RooIntegrator1D
170  switch(_case= limitsCase()) {
171  case ClosedBothEnds:
172  // both limits are finite: use the plain trapezoid integrator
174  break;
175  case OpenBothEnds:
176  // both limits are infinite: integrate over (-1,+1) using
177  // the plain trapezoid integrator...
179  // ...and integrate the infinite tails using the midpoint integrator
182  break;
183  case OpenBelowSpansZero:
184  // xmax >= 0 so integrate from (-inf,-1) and (-1,xmax)
187  break;
188  case OpenBelow:
189  // xmax < 0 so integrate from (-inf,xmax)
191  break;
192  case OpenAboveSpansZero:
193  // xmin <= 0 so integrate from (xmin,+1) and (+1,+inf)
196  break;
197  case OpenAbove:
198  // xmin > 0 so integrate from (xmin,+inf)
200  break;
201  case Invalid:
202  default:
203  _valid= kFALSE;
204  }
205 }
206 
207 
208 ////////////////////////////////////////////////////////////////////////////////
209 /// Destructor
210 
212 {
213  if(0 != _integrator1) delete _integrator1;
214  if(0 != _integrator2) delete _integrator2;
215  if(0 != _integrator3) delete _integrator3;
216  if(0 != _function) delete _function;
217 }
218 
219 
220 ////////////////////////////////////////////////////////////////////////////////
221 /// Change our integration limits. Return kTRUE if the new limits are
222 /// ok, or otherwise kFALSE. Always returns kFALSE and does nothing
223 /// if this object was constructed to always use our integrand's limits.
224 
226 {
227  if(_useIntegrandLimits) {
228  oocoutE((TObject*)0,Integration) << "RooIntegrator1D::setLimits: cannot override integrand's limits" << endl;
229  return kFALSE;
230  }
231 
232  _xmin= *xmin;
233  _xmax= *xmax;
234  return checkLimits();
235 }
236 
237 
238 ////////////////////////////////////////////////////////////////////////////////
239 /// Check if the limits are valid. For this integrator all limit configurations
240 /// are valid, but if the limits change between two calculate() calls it
241 /// may be necessary to reconfigure (e.g. if an open ended range becomes
242 /// a closed range
243 
245 {
246  // Has either limit changed?
247  if (_useIntegrandLimits) {
248  if(_xmin == integrand()->getMinLimit(0) &&
249  _xmax == integrand()->getMaxLimit(0)) return kTRUE;
250  }
251 
252  // The limits have changed: can we use the same strategy?
253  if(limitsCase() != _case) {
254  // Reinitialize embedded integrators, will automatically propagate new limits
255  const_cast<RooImproperIntegrator1D*>(this)->initialize() ;
256  return kTRUE ;
257  }
258 
259  // Reuse our existing integrators by updating their limits
260  switch(_case) {
261  case ClosedBothEnds:
263  break;
264  case OpenBothEnds:
265  // nothing has changed
266  break;
267  case OpenBelowSpansZero:
269  break;
270  case OpenBelow:
272  break;
273  case OpenAboveSpansZero:
275  break;
276  case OpenAbove:
278  break;
279  case Invalid:
280  default:
281  return kFALSE;
282  }
283  return kTRUE;
284 }
285 
286 
287 ////////////////////////////////////////////////////////////////////////////////
288 /// Classify the type of limits we have: OpenBothEnds,ClosedBothEnds,OpenBelow or OpenAbove.
289 
291 {
292  // Analyze the specified limits to determine which case applies.
293  if(0 == integrand() || !integrand()->isValid()) return Invalid;
294 
295  if (_useIntegrandLimits) {
296  _xmin= integrand()->getMinLimit(0);
297  _xmax= integrand()->getMaxLimit(0);
298  }
299 
302  if(!inf1 && !inf2) {
303  // both limits are finite
304  return ClosedBothEnds;
305  }
306  else if(inf1 && inf2) {
307  // both limits are infinite
308  return OpenBothEnds;
309  }
310  else if(inf1) { // inf2==false
311  if(_xmax >= 0) {
312  return OpenBelowSpansZero;
313  }
314  else {
315  return OpenBelow;
316  }
317  }
318  else { // inf1==false && inf2==true
319  if(_xmin <= 0) {
320  return OpenAboveSpansZero;
321  }
322  else {
323  return OpenAbove;
324  }
325  }
326  // return Invalid; OSF-CC: Statement unreachable
327 }
328 
329 
330 ////////////////////////////////////////////////////////////////////////////////
331 /// Calculate the integral at the given parameter values of the function binding
332 
334 {
335  Double_t result(0);
336  if(0 != _integrator1) result+= _integrator1->integral(yvec);
337  if(0 != _integrator2) result+= _integrator2->integral(yvec);
338  if(0 != _integrator3) result+= _integrator3->integral(yvec);
339  return result;
340 }
std::string GetName(const std::string &scope_name)
Definition: Cppyy.cxx:145
RooAbsIntegrator is the abstract interface for integrators of real-valued functions that implement th...
float xmin
Definition: THbookFile.cxx:93
LimitsCase limitsCase() const
Classify the type of limits we have: OpenBothEnds,ClosedBothEnds,OpenBelow or OpenAbove.
RooNumIntConfig holds the configuration parameters of the various numeric integrators used by RooReal...
RooNumIntFactory is a factory to instantiate numeric integrators from a given function binding and a ...
bool Bool_t
Definition: RtypesCore.h:59
static Int_t isInfinite(Double_t x)
Return true if x is infinite by RooNumBer internal specification.
Definition: RooNumber.cxx:58
static void registerIntegrator(RooNumIntFactory &fact)
Register RooImproperIntegrator1D, its parameters and capabilities with RooNumIntFactory.
STL namespace.
Bool_t isValid() const
void initialize(const RooAbsFunc *function=0)
Initialize the integrator, construct and initialize subintegrators.
void Class()
Definition: Class.C:29
#define oocoutE(o, a)
Definition: RooMsgService.h:47
Bool_t setLimits(Double_t *xmin, Double_t *xmax)
Change our integration limits.
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:146
Lightweight function binding that returns the inverse of an input function binding Apply the change o...
virtual Double_t integral(const Double_t *yvec=0)
Calculate numeric integral at given set of function binding parameters.
RooImproperIntegrator1D()
Default constructor.
virtual Bool_t checkLimits() const
Check if the limits are valid.
virtual Double_t getMinLimit(UInt_t dimension) const =0
Special numeric integrator that can handle integrals over open domains.
float xmax
Definition: THbookFile.cxx:93
Bool_t setLimits(Double_t *xmin, Double_t *xmax)
Change our integration limits.
const Bool_t kFALSE
Definition: RtypesCore.h:92
const RooAbsFunc * integrand() const
RooIntegrator1D implements an adaptive one-dimensional numerical integration algorithm.
#define ClassImp(name)
Definition: Rtypes.h:336
virtual Double_t getMaxLimit(UInt_t dimension) const =0
double Double_t
Definition: RtypesCore.h:55
Mother of all ROOT objects.
Definition: TObject.h:37
virtual Double_t integral(const Double_t *yvec=0)
Calculate the integral at the given parameter values of the function binding.
const char * proto
Definition: civetweb.c:11652
double result[121]
virtual RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const
Return clone of integrator with given function and configuration. Needed by RooNumIntFactory.
virtual ~RooImproperIntegrator1D()
Destructor.
const Bool_t kTRUE
Definition: RtypesCore.h:91
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23
Bool_t storeProtoIntegrator(RooAbsIntegrator *proto, const RooArgSet &defConfig, const char *depName="")
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...