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