Logo ROOT   6.10/09
Reference Guide
RooAddition.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 RooAddition.cxx
19 \class RooAddition
20 \ingroup Roofitcore
21 
22 RooAddition calculates the sum of a set of RooAbsReal terms, or
23 when constructed with two sets, it sums the product of the terms
24 in the two sets. This class does not (yet) do any smart handling of integrals,
25 i.e. all integrals of the product are handled numerically
26 **/
27 
28 
29 #include "RooFit.h"
30 
31 #include "Riostream.h"
32 #include <math.h>
33 #include <memory>
34 #include <list>
35 #include <algorithm>
36 using namespace std ;
37 
38 #include "RooAddition.h"
39 #include "RooProduct.h"
40 #include "RooAbsReal.h"
41 #include "RooErrorHandler.h"
42 #include "RooArgSet.h"
43 #include "RooNameReg.h"
44 #include "RooNLLVar.h"
45 #include "RooChi2Var.h"
46 #include "RooMsgService.h"
47 
49 ;
50 
51 
52 ////////////////////////////////////////////////////////////////////////////////
53 
55  : _setIter( _set.createIterator() )
56 {
57 }
58 
59 
60 
61 ////////////////////////////////////////////////////////////////////////////////
62 /// Constructor with a single set of RooAbsReals. The value of the function will be
63 /// the sum of the values in sumSet. If takeOwnership is true the RooAddition object
64 /// will take ownership of the arguments in sumSet
65 
66 RooAddition::RooAddition(const char* name, const char* title, const RooArgList& sumSet, Bool_t takeOwnership)
67  : RooAbsReal(name, title)
68  , _set("!set","set of components",this)
69  , _setIter( _set.createIterator() ) // yes, _setIter is defined _after_ _set ;-)
70  , _cacheMgr(this,10)
71 {
72  std::unique_ptr<TIterator> inputIter( sumSet.createIterator() );
73  RooAbsArg* comp ;
74  while((comp = (RooAbsArg*)inputIter->Next())) {
75  if (!dynamic_cast<RooAbsReal*>(comp)) {
76  coutE(InputArguments) << "RooAddition::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
77  << " is not of type RooAbsReal" << endl ;
79  }
80  _set.add(*comp) ;
81  if (takeOwnership) _ownedList.addOwned(*comp) ;
82  }
83 
84 }
85 
86 
87 
88 ////////////////////////////////////////////////////////////////////////////////
89 /// Constructor with two set of RooAbsReals. The value of the function will be
90 ///
91 /// A = sum_i sumSet1(i)*sumSet2(i)
92 ///
93 /// If takeOwnership is true the RooAddition object will take ownership of the arguments in sumSet
94 
95 RooAddition::RooAddition(const char* name, const char* title, const RooArgList& sumSet1, const RooArgList& sumSet2, Bool_t takeOwnership)
96  : RooAbsReal(name, title)
97  , _set("!set","set of components",this)
98  , _setIter( _set.createIterator() ) // yes, _setIter is defined _after_ _set ;-)
99  , _cacheMgr(this,10)
100 {
101  if (sumSet1.getSize() != sumSet2.getSize()) {
102  coutE(InputArguments) << "RooAddition::ctor(" << GetName() << ") ERROR: input lists should be of equal length" << endl ;
104  }
105 
106  std::unique_ptr<TIterator> inputIter1( sumSet1.createIterator() );
107  std::unique_ptr<TIterator> inputIter2( sumSet2.createIterator() );
108  RooAbsArg *comp1(0),*comp2(0) ;
109  while((comp1 = (RooAbsArg*)inputIter1->Next())) {
110  if (!dynamic_cast<RooAbsReal*>(comp1)) {
111  coutE(InputArguments) << "RooAddition::ctor(" << GetName() << ") ERROR: component " << comp1->GetName()
112  << " in first list is not of type RooAbsReal" << endl ;
114  }
115  comp2 = (RooAbsArg*)inputIter2->Next();
116  if (!dynamic_cast<RooAbsReal*>(comp2)) {
117  coutE(InputArguments) << "RooAddition::ctor(" << GetName() << ") ERROR: component " << comp2->GetName()
118  << " in first list is not of type RooAbsReal" << endl ;
120  }
121  // TODO: add flag to RooProduct c'tor to make it assume ownership...
122  TString _name(name);
123  _name.Append( "_[");
124  _name.Append(comp1->GetName());
125  _name.Append( "_x_");
126  _name.Append(comp2->GetName());
127  _name.Append( "]");
128  RooProduct *prod = new RooProduct( _name, _name , RooArgSet(*comp1, *comp2) /*, takeOwnership */ ) ;
129  _set.add(*prod);
130  _ownedList.addOwned(*prod) ;
131  if (takeOwnership) {
132  _ownedList.addOwned(*comp1) ;
133  _ownedList.addOwned(*comp2) ;
134  }
135  }
136 }
137 
138 
139 
140 ////////////////////////////////////////////////////////////////////////////////
141 /// Copy constructor
142 
143 RooAddition::RooAddition(const RooAddition& other, const char* name)
144  : RooAbsReal(other, name)
145  , _set("!set",this,other._set)
146  , _setIter( _set.createIterator() ) // yes, _setIter is defined _after_ _set ;-)
147  , _cacheMgr(other._cacheMgr,this)
148 {
149  // Member _ownedList is intentionally not copy-constructed -- ownership is not transferred
150 }
151 
152 
153 ////////////////////////////////////////////////////////////////////////////////
154 
156 { // Destructor
157  delete _setIter ;
158 }
159 
160 ////////////////////////////////////////////////////////////////////////////////
161 /// Calculate and return current value of self
162 
164 {
165  Double_t sum(0);
166  const RooArgSet* nset = _set.nset() ;
167 
168 // cout << "RooAddition::eval sum = " ;
169 
170  RooFIter setIter = _set.fwdIterator() ;
171  RooAbsReal* comp ;
172  while((comp=(RooAbsReal*)setIter.next())) {
173  Double_t tmp = comp->getVal(nset) ;
174 // cout << tmp << " " ;
175  sum += tmp ;
176  }
177 // cout << " = " << sum << endl ;
178  return sum ;
179 }
180 
181 
182 ////////////////////////////////////////////////////////////////////////////////
183 /// Return the default error level for MINUIT error analysis
184 /// If the addition contains one or more RooNLLVars and
185 /// no RooChi2Vars, return the defaultErrorLevel() of
186 /// RooNLLVar. If the addition contains one ore more RooChi2Vars
187 /// and no RooNLLVars, return the defaultErrorLevel() of
188 /// RooChi2Var. If the addition contains neither or both
189 /// issue a warning message and return a value of 1
190 
192 {
193  RooAbsReal* nllArg(0) ;
194  RooAbsReal* chi2Arg(0) ;
195 
196  RooAbsArg* arg ;
197 
198  RooArgSet* comps = getComponents() ;
199  TIterator* iter = comps->createIterator() ;
200  while((arg=(RooAbsArg*)iter->Next())) {
201  if (dynamic_cast<RooNLLVar*>(arg)) {
202  nllArg = (RooAbsReal*)arg ;
203  }
204  if (dynamic_cast<RooChi2Var*>(arg)) {
205  chi2Arg = (RooAbsReal*)arg ;
206  }
207  }
208  delete iter ;
209  delete comps ;
210 
211  if (nllArg && !chi2Arg) {
212  coutI(Fitting) << "RooAddition::defaultErrorLevel(" << GetName()
213  << ") Summation contains a RooNLLVar, using its error level" << endl ;
214  return nllArg->defaultErrorLevel() ;
215  } else if (chi2Arg && !nllArg) {
216  coutI(Fitting) << "RooAddition::defaultErrorLevel(" << GetName()
217  << ") Summation contains a RooChi2Var, using its error level" << endl ;
218  return chi2Arg->defaultErrorLevel() ;
219  } else if (!nllArg && !chi2Arg) {
220  coutI(Fitting) << "RooAddition::defaultErrorLevel(" << GetName() << ") WARNING: "
221  << "Summation contains neither RooNLLVar nor RooChi2Var server, using default level of 1.0" << endl ;
222  } else {
223  coutI(Fitting) << "RooAddition::defaultErrorLevel(" << GetName() << ") WARNING: "
224  << "Summation contains BOTH RooNLLVar and RooChi2Var server, using default level of 1.0" << endl ;
225  }
226 
227  return 1.0 ;
228 }
229 
230 
231 ////////////////////////////////////////////////////////////////////////////////
232 
234 {
235  _setIter->Reset() ;
236 
237  RooAbsReal* arg;
238  while((arg=(RooAbsReal*)_setIter->Next())) {
239  arg->enableOffsetting(flag) ;
240  }
241 }
242 
243 
244 
245 ////////////////////////////////////////////////////////////////////////////////
246 
248 {
249  _setIter->Reset() ;
250 
251  RooAbsReal* arg;
252  while((arg=(RooAbsReal*)_setIter->Next())) {
253  arg->setData(data,cloneData) ;
254  }
255  return kTRUE ;
256 }
257 
258 
259 
260 ////////////////////////////////////////////////////////////////////////////////
261 
262 void RooAddition::printMetaArgs(ostream& os) const
263 {
264  _setIter->Reset() ;
265 
266  Bool_t first(kTRUE) ;
267  RooAbsArg* arg;
268  while((arg=(RooAbsArg*)_setIter->Next())) {
269  if (!first) { os << " + " ;
270  } else { first = kFALSE ;
271  }
272  os << arg->GetName() ;
273  }
274  os << " " ;
275 }
276 
277 ////////////////////////////////////////////////////////////////////////////////
278 
279 Int_t RooAddition::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
280 {
281  // we always do things ourselves -- actually, always delegate further down the line ;-)
282  analVars.add(allVars);
283 
284  // check if we already have integrals for this combination of factors
285  Int_t sterileIndex(-1);
286  CacheElem* cache = (CacheElem*) _cacheMgr.getObj(&analVars,&analVars,&sterileIndex,RooNameReg::ptr(rangeName));
287  if (cache!=0) {
288  Int_t code = _cacheMgr.lastIndex();
289  return code+1;
290  }
291 
292  // we don't, so we make it right here....
293  cache = new CacheElem;
294  _setIter->Reset();
295  RooAbsReal *arg(0);
296  while( (arg=(RooAbsReal*)_setIter->Next())!=0 ) { // checked in c'tor that this will work...
297  RooAbsReal *I = arg->createIntegral(analVars,rangeName);
298  cache->_I.addOwned(*I);
299  }
300 
301  Int_t code = _cacheMgr.setObj(&analVars,&analVars,(RooAbsCacheElement*)cache,RooNameReg::ptr(rangeName));
302  return 1+code;
303 }
304 
305 ////////////////////////////////////////////////////////////////////////////////
306 /// Calculate integral internally from appropriate integral cache
307 
308 Double_t RooAddition::analyticalIntegral(Int_t code, const char* rangeName) const
309 {
310  // note: rangeName implicit encoded in code: see _cacheMgr.setObj in getPartIntList...
311  CacheElem *cache = (CacheElem*) _cacheMgr.getObjByIndex(code-1);
312  if (cache==0) {
313  // cache got sterilized, trigger repopulation of this slot, then try again...
314  std::unique_ptr<RooArgSet> vars( getParameters(RooArgSet()) );
315  std::unique_ptr<RooArgSet> iset( _cacheMgr.nameSet2ByIndex(code-1)->select(*vars) );
317  Int_t code2 = getAnalyticalIntegral(*iset,dummy,rangeName);
318  assert(code==code2); // must have revived the right (sterilized) slot...
319  return analyticalIntegral(code2,rangeName);
320  }
321  assert(cache!=0);
322 
323  // loop over cache, and sum...
324  std::unique_ptr<TIterator> iter( cache->_I.createIterator() );
325  RooAbsReal *I;
326  double result(0);
327  while ( ( I=(RooAbsReal*)iter->Next() ) != 0 ) result += I->getVal();
328  return result;
329 
330 }
331 
332 
333 
334 ////////////////////////////////////////////////////////////////////////////////
335 
336 std::list<Double_t>* RooAddition::binBoundaries(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const
337 {
338  std::list<Double_t>* sumBinB = 0 ;
339  Bool_t needClean(kFALSE) ;
340 
341  RooFIter iter = _set.fwdIterator() ;
342  RooAbsReal* func ;
343  // Loop over components pdf
344  while((func=(RooAbsReal*)iter.next())) {
345 
346  std::list<Double_t>* funcBinB = func->binBoundaries(obs,xlo,xhi) ;
347 
348  // Process hint
349  if (funcBinB) {
350  if (!sumBinB) {
351  // If this is the first hint, then just save it
352  sumBinB = funcBinB ;
353  } else {
354 
355  std::list<Double_t>* newSumBinB = new std::list<Double_t>(sumBinB->size()+funcBinB->size()) ;
356 
357  // Merge hints into temporary array
358  merge(funcBinB->begin(),funcBinB->end(),sumBinB->begin(),sumBinB->end(),newSumBinB->begin()) ;
359 
360  // Copy merged array without duplicates to new sumBinBArrau
361  delete sumBinB ;
362  delete funcBinB ;
363  sumBinB = newSumBinB ;
364  needClean = kTRUE ;
365  }
366  }
367  }
368 
369  // Remove consecutive duplicates
370  if (needClean) {
371  std::list<Double_t>::iterator new_end = unique(sumBinB->begin(),sumBinB->end()) ;
372  sumBinB->erase(new_end,sumBinB->end()) ;
373  }
374 
375  return sumBinB ;
376 }
377 
378 
379 //_____________________________________________________________________________B
381 {
382  // If all components that depend on obs are binned that so is the product
383 
384  RooFIter iter = _set.fwdIterator() ;
385  RooAbsReal* func ;
386  while((func=(RooAbsReal*)iter.next())) {
387  if (func->dependsOn(obs) && !func->isBinnedDistribution(obs)) {
388  return kFALSE ;
389  }
390  }
391 
392  return kTRUE ;
393 }
394 
395 
396 
397 
398 ////////////////////////////////////////////////////////////////////////////////
399 
400 std::list<Double_t>* RooAddition::plotSamplingHint(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const
401 {
402  std::list<Double_t>* sumHint = 0 ;
403  Bool_t needClean(kFALSE) ;
404 
405  RooFIter iter = _set.fwdIterator() ;
406  RooAbsReal* func ;
407  // Loop over components pdf
408  while((func=(RooAbsReal*)iter.next())) {
409 
410  std::list<Double_t>* funcHint = func->plotSamplingHint(obs,xlo,xhi) ;
411 
412  // Process hint
413  if (funcHint) {
414  if (!sumHint) {
415 
416  // If this is the first hint, then just save it
417  sumHint = funcHint ;
418 
419  } else {
420 
421  std::list<Double_t>* newSumHint = new std::list<Double_t>(sumHint->size()+funcHint->size()) ;
422 
423  // Merge hints into temporary array
424  merge(funcHint->begin(),funcHint->end(),sumHint->begin(),sumHint->end(),newSumHint->begin()) ;
425 
426  // Copy merged array without duplicates to new sumHintArrau
427  delete sumHint ;
428  sumHint = newSumHint ;
429  needClean = kTRUE ;
430  }
431  }
432  }
433 
434  // Remove consecutive duplicates
435  if (needClean) {
436  std::list<Double_t>::iterator new_end = unique(sumHint->begin(),sumHint->end()) ;
437  sumHint->erase(new_end,sumHint->end()) ;
438  }
439 
440  return sumHint ;
441 }
442 
443 
444 
445 ////////////////////////////////////////////////////////////////////////////////
446 /// Return list of all RooAbsArgs in cache element
447 
449 {
450  RooArgList ret(_I) ;
451  return ret ;
452 }
453 
455 {
456  // Destructor
457 }
458 
459 
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TIterator * createIterator(Bool_t dir=kIterForward) const
static long int sum(long int i)
Definition: Factory.cxx:2162
#define coutE(a)
Definition: RooMsgService.h:34
void printMetaArgs(std::ostream &os) const
TIterator * _setIter
Definition: RooAddition.h:63
static void softAbort()
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:86
Bool_t dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=0, Bool_t valueOnly=kFALSE) const
Test whether we depend on (ie, are served by) any object in the specified collection.
Definition: RooAbsArg.cxx:744
virtual void Reset()=0
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
#define coutI(a)
Definition: RooMsgService.h:31
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
virtual std::list< Double_t > * binBoundaries(RooAbsRealLValue &, Double_t, Double_t) const
Definition: RooAbsReal.h:278
virtual RooArgList containedArgs(Action)
Return list of all RooAbsArgs in cache element.
T * getObjByIndex(Int_t index) const
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
STL namespace.
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Calculate integral internally from appropriate integral cache.
Iterator abstract base class.
Definition: TIterator.h:30
virtual std::list< Double_t > * plotSamplingHint(RooAbsRealLValue &, Double_t, Double_t) const
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=0, const TNamed *isetRangeName=0)
friend class RooArgSet
Definition: RooAbsArg.h:469
Double_t evaluate() const
Calculate and return current value of self.
RooArgList _ownedList
Definition: RooAddition.h:61
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
virtual void enableOffsetting(Bool_t)
Definition: RooAbsReal.h:307
virtual Bool_t setData(RooAbsData &, Bool_t=kTRUE)
Definition: RooAbsReal.h:305
Bool_t setData(RooAbsData &data, Bool_t cloneData=kTRUE)
RooListProxy _set
Definition: RooAddition.h:62
Iterator over set.
Definition: RooAddition.h:65
const RooArgSet * nset() const
Definition: RooAbsProxy.h:46
Int_t getSize() const
RooAbsReal * createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create an object that represents the integral of the function over one or more observables listed in ...
Definition: RooAbsReal.cxx:501
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &numVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported...
static const TNamed * ptr(const char *stringPtr)
Return a unique TNamed pointer for given C++ string.
Definition: RooNameReg.cxx:125
RooObjCacheManager _cacheMgr
Definition: RooAddition.h:72
virtual void enableOffsetting(Bool_t)
RooArgSet * getComponents() const
Definition: RooAbsArg.cxx:688
Bool_t isBinnedDistribution(const RooArgSet &obs) const
RooAbsCacheElement is the abstract base class for objects to be stored in RooAbsCache cache manager o...
Int_t lastIndex() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
RooProduct a RooAbsReal implementation that represent the product of a given set of other RooAbsReal ...
Definition: RooProduct.h:32
RooAbsArg * next()
const Bool_t kFALSE
Definition: RtypesCore.h:92
virtual Double_t defaultErrorLevel() const
Return the default error level for MINUIT error analysis If the addition contains one or more RooNLLV...
virtual Bool_t isBinnedDistribution(const RooArgSet &) const
Definition: RooAbsReal.h:277
const RooNameSet * nameSet2ByIndex(Int_t index) const
#define ClassImp(name)
Definition: Rtypes.h:336
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
RooFIter fwdIterator() const
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don&#39;t match any of...
Definition: RooAbsArg.cxx:560
static RooMathCoreReg dummy
double func(double *x, double *p)
Definition: stressTF1.cxx:213
RooArgSet * select(const RooArgSet &list) const
Construct a RooArgSet of objects in input &#39;list&#39; whose names match to those in the internal name list...
Definition: RooNameSet.cxx:189
virtual std::list< Double_t > * plotSamplingHint(RooAbsRealLValue &, Double_t, Double_t) const
Definition: RooAbsReal.h:279
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
virtual Double_t defaultErrorLevel() const
Definition: RooAbsReal.h:189
virtual TObject * Next()=0
virtual ~RooAddition()
double result[121]
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=0)
RooAddition calculates the sum of a set of RooAbsReal terms, or when constructed with two sets...
Definition: RooAddition.h:26
Definition: first.py:1
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
#define I(x, y, z)
virtual std::list< Double_t > * binBoundaries(RooAbsRealLValue &, Double_t, Double_t) const
const Bool_t kTRUE
Definition: RtypesCore.h:91