Logo ROOT   6.14/05
Reference Guide
RooFFTConvPdf.cxx
Go to the documentation of this file.
1 
2  /*****************************************************************************
3  * Project: RooFit *
4  * *
5  * Copyright (c) 2000-2005, Regents of the University of California *
6  * and Stanford University. All rights reserved. *
7  * *
8  * Redistribution and use in source and binary forms, *
9  * with or without modification, are permitted according to the terms *
10  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
11  *****************************************************************************/
12 
13 //////////////////////////////////////////////////////////////////////////////
14 //
15  //
16  // This class implement a generic one-dimensional numeric convolution of two p.d.f.
17  // and can convolve any two RooAbsPdfs. The class exploits the convolution theorem
18  //
19  // f(x) (*) g(x) --F--> f(k_i) * g(k_i)
20  //
21  // and calculate the convolution by calculate a Real->Complex FFT of both input p.d.fs
22  // multiplying the complex coefficients and performing the reverse Complex->Real FFT
23  // to get the result in the input space. This class using the ROOT FFT Interface to
24  // the (free) FFTW3 package (www.fftw.org) and requires that your ROOT installation is
25  // compiled with the --enable-fftw3 option (instructions for Linux follow)
26  //
27  // Note that the performance in terms of speed and stability of RooFFTConvPdf is
28  // vastly superior to that of RooNumConvPdf
29  //
30  // An important feature of FFT convolutions is that the observable is treated in a
31  // cyclical way. This is correct & desirable behavior for cyclical observables such as angles,
32  // but it may not be for other observables. The effect that is observed is that if
33  // p.d.f is zero at xMin and non-zero at xMax some spillover occurs and
34  // a rising tail may appear at xMin. This effect can be reduced or eliminated by
35  // introducing a buffer zone in the FFT calculation. If this feature is activated
36  // input the sampling array for the FFT calculation is extended in both directions
37  // and filled with repetitions of the lowest bin value and highest bin value
38  // respectively. The buffer bins are stripped again when the FFT output values
39  // are transferred to the p.d.f cache. The default buffer size is 10% of the
40  // observable domain size and can be changed with setBufferFraction() member function.
41  //
42  // This class is a caching p.d.f inheriting from RooAbsCachedPdf. If this p.d.f
43  // is evaluated for a particular value of x, the FFT calculate the values for the
44  // p.d.f at all points in observables space for the given choice of parameters,
45  // which are stored in the cache. Subsequent evaluations of RooFFTConvPdf with
46  // identical parameters will retrieve results from the cache. If one or more
47  // of the parameters change, the cache will be updated.
48  //
49  // The sampling density of the cache is controlled by the binning of the
50  // the convolution observable, which can be changed from RooRealVar::setBins(N)
51  // For good results N should be large (>1000). Additional interpolation of
52  // cache values may improve the result if courser binning are chosen. These can be
53  // set in the constructor or through the setInterpolationOrder() member function.
54  // For N>1000 interpolation will not substantially improve the performance.
55  //
56  // Additionial information on caching activities can be displayed by monitoring
57  // the message stream with topic "Caching" at the INFO level, i.e.
58  // do RooMsgService::instance().addStream(RooMsgService::INFO,Topic("Caching"))
59  // to see these message on stdout
60  //
61  // Multi-dimensional convolutions are not supported yet, but will be in the future
62  // as FFTW can calculate them
63  //
64  // ---
65  //
66  // Installing a copy of FFTW on Linux and compiling ROOT to use it
67  //
68  // 1) Go to www.fftw.org and download the latest stable version (a .tar.gz file)
69  //
70  // If you have root access to your machine and want to make a system installation of FFTW
71  //
72  // 2) Untar fftw-XXX.tar.gz in /tmp, cd into the untarred directory
73  // and type './configure' followed by 'make install'.
74  // This will install fftw in /usr/local/bin,lib etc...
75  //
76  // 3) Start from a source installation of ROOT. If you now have a binary distribution,
77  // first download a source tar ball from root.cern.ch for your ROOT version and untar it.
78  // Run 'configure', following the instruction from 'configure --help' but be sure run 'configure'
79  // with additional flags '--enable-fftw3' and '--enable-roofit', then run 'make'
80  //
81  //
82  // If you do not have root access and want to make a private installation of FFTW
83  //
84  // 2) Make a private install area for FFTW, e.g. /home/myself/fftw
85  //
86  // 3) Untar fftw-XXX.tar.gz in /tmp, cd into the untarred directory
87  // and type './configure --prefix=/home/myself/fftw' followed by 'make install'.
88  // Substitute /home/myself/fftw with a directory of your choice. This
89  // procedure will install FFTW in the location designated by you
90  //
91  // 4) Start from a source installation of ROOT. If you now have a binary distribution,
92  // first download a source tar ball from root.cern.ch for your ROOT version and untar it.
93  // Run 'configure', following the instruction from 'configure --help' but be sure run 'configure'
94  // with additional flags
95  // '--enable-fftw3',
96  // '--with-fftw3-incdir=/home/myself/fftw/include',
97  // '--width-fftw3-libdir=/home/myself/fftw/lib' and
98  // '--enable-roofit'
99  // Then run 'make'
100 //
101 
102 
103 #include "Riostream.h"
104 
105 #include "RooFit.h"
106 #include "RooFFTConvPdf.h"
107 #include "RooAbsReal.h"
108 #include "RooMsgService.h"
109 #include "RooDataHist.h"
110 #include "RooHistPdf.h"
111 #include "RooRealVar.h"
112 #include "TComplex.h"
113 #include "TVirtualFFT.h"
114 #include "RooGenContext.h"
115 #include "RooConvGenContext.h"
116 #include "RooBinning.h"
117 #include "RooLinearVar.h"
118 #include "RooCustomizer.h"
119 #include "RooGlobalFunc.h"
120 #include "RooLinearVar.h"
121 #include "RooConstVar.h"
122 #include "TClass.h"
123 #include "TSystem.h"
124 
125 using namespace std ;
126 
128 
129 
130 
131 ////////////////////////////////////////////////////////////////////////////////
132 /// Constructor for convolution of pdf1 (x) pdf2 in observable convVar. The binning used for the FFT sampling is controlled
133 /// by the binning named "cache" in the convolution observable. The resulting FFT convolved histogram is interpolated at
134 /// order 'ipOrder' A minimum binning of 1000 bins is recommended.
135 
136 RooFFTConvPdf::RooFFTConvPdf(const char *name, const char *title, RooRealVar& convVar, RooAbsPdf& pdf1, RooAbsPdf& pdf2, Int_t ipOrder) :
137  RooAbsCachedPdf(name,title,ipOrder),
138  _x("!x","Convolution Variable",this,convVar),
139  _xprime("!xprime","External Convolution Variable",this,0),
140  _pdf1("!pdf1","pdf1",this,pdf1,kFALSE),
141  _pdf2("!pdf2","pdf2",this,pdf2,kFALSE),
142  _params("!params","effective parameters",this),
143  _bufFrac(0.1),
144  _bufStrat(Extend),
145  _shift1(0),
146  _shift2(0),
147  _cacheObs("!cacheObs","Cached observables",this,kFALSE,kFALSE)
148  {
149  if (!convVar.hasBinning("cache")) {
150  convVar.setBinning(convVar.getBinning(),"cache") ;
151  }
152 
153  _shift2 = (convVar.getMax("cache")+convVar.getMin("cache"))/2 ;
154 
155  calcParams() ;
156 
157  }
158 
159 
160 
161 ////////////////////////////////////////////////////////////////////////////////
162 /// Constructor for convolution of pdf1 (x) pdf2 in observable convVar. The binning used for the FFT sampling is controlled
163 /// by the binning named "cache" in the convolution observable. The resulting FFT convolved histogram is interpolated at
164 /// order 'ipOrder' A minimum binning of 1000 bins is recommended.
165 
166 RooFFTConvPdf::RooFFTConvPdf(const char *name, const char *title, RooAbsReal& pdfConvVar, RooRealVar& convVar, RooAbsPdf& pdf1, RooAbsPdf& pdf2, Int_t ipOrder) :
167  RooAbsCachedPdf(name,title,ipOrder),
168  _x("!x","Convolution Variable",this,convVar,kFALSE,kFALSE),
169  _xprime("!xprime","External Convolution Variable",this,pdfConvVar),
170  _pdf1("!pdf1","pdf1",this,pdf1,kFALSE),
171  _pdf2("!pdf2","pdf2",this,pdf2,kFALSE),
172  _params("!params","effective parameters",this),
173  _bufFrac(0.1),
174  _bufStrat(Extend),
175  _shift1(0),
176  _shift2(0),
177  _cacheObs("!cacheObs","Cached observables",this,kFALSE,kFALSE)
178  {
179  if (!convVar.hasBinning("cache")) {
180  convVar.setBinning(convVar.getBinning(),"cache") ;
181  }
182 
183  _shift2 = (convVar.getMax("cache")+convVar.getMin("cache"))/2 ;
184 
185  calcParams() ;
186  }
187 
188 
189 
190 ////////////////////////////////////////////////////////////////////////////////
191 /// Copy constructor
192 
193 RooFFTConvPdf::RooFFTConvPdf(const RooFFTConvPdf& other, const char* name) :
194  RooAbsCachedPdf(other,name),
195  _x("!x",this,other._x),
196  _xprime("!xprime",this,other._xprime),
197  _pdf1("!pdf1",this,other._pdf1),
198  _pdf2("!pdf2",this,other._pdf2),
199  _params("!params",this,other._params),
200  _bufFrac(other._bufFrac),
201  _bufStrat(other._bufStrat),
202  _shift1(other._shift1),
203  _shift2(other._shift2),
204  _cacheObs("!cacheObs",this,other._cacheObs)
205  {
206  }
207 
208 
209 
210 ////////////////////////////////////////////////////////////////////////////////
211 /// Destructor
212 
214 {
215 }
216 
217 
218 
219 ////////////////////////////////////////////////////////////////////////////////
220 /// Return base name component for cache components in this case 'PDF1_CONV_PDF2'
221 
222 const char* RooFFTConvPdf::inputBaseName() const
223 {
224  static TString name ;
225  name = _pdf1.arg().GetName() ;
226  name.Append("_CONV_") ;
227  name.Append(_pdf2.arg().GetName()) ;
228  return name.Data() ;
229 }
230 
231 
232 
233 
234 ////////////////////////////////////////////////////////////////////////////////
235 /// Return specialized cache subclass for FFT calculations
236 
238 {
239  return new FFTCacheElem(*this,nset) ;
240 }
241 
242 
243 
244 
245 ////////////////////////////////////////////////////////////////////////////////
246 /// Clone input pdf and attach to dataset
247 
249  PdfCacheElem(self,nsetIn),
250  fftr2c1(0),fftr2c2(0),fftc2r(0)
251 {
252  RooAbsPdf* clonePdf1 = (RooAbsPdf*) self._pdf1.arg().cloneTree() ;
253  RooAbsPdf* clonePdf2 = (RooAbsPdf*) self._pdf2.arg().cloneTree() ;
254  clonePdf1->attachDataSet(*hist()) ;
255  clonePdf2->attachDataSet(*hist()) ;
256 
257  // Shift observable
258  RooRealVar* convObs = (RooRealVar*) hist()->get()->find(self._x.arg().GetName()) ;
259 
260  // Install FFT reference range
261  string refName = Form("refrange_fft_%s",self.GetName()) ;
262  convObs->setRange(refName.c_str(),convObs->getMin(),convObs->getMax()) ;
263 
264  if (self._shift1!=0) {
265  RooLinearVar* shiftObs1 = new RooLinearVar(Form("%s_shifted_FFTBuffer1",convObs->GetName()),"shiftObs1",
266  *convObs,RooFit::RooConst(1),RooFit::RooConst(-1*self._shift1)) ;
267 
268  RooArgSet clonedBranches1 ;
269  RooCustomizer cust(*clonePdf1,"fft") ;
270  cust.replaceArg(*convObs,*shiftObs1) ;
271 
272  pdf1Clone = (RooAbsPdf*) cust.build() ;
273 
274  pdf1Clone->addOwnedComponents(*shiftObs1) ;
275  pdf1Clone->addOwnedComponents(*clonePdf1) ;
276 
277  } else {
278  pdf1Clone = clonePdf1 ;
279  }
280 
281  if (self._shift2!=0) {
282  RooLinearVar* shiftObs2 = new RooLinearVar(Form("%s_shifted_FFTBuffer2",convObs->GetName()),"shiftObs2",
283  *convObs,RooFit::RooConst(1),RooFit::RooConst(-1*self._shift2)) ;
284 
285  RooArgSet clonedBranches2 ;
286  RooCustomizer cust(*clonePdf2,"fft") ;
287  cust.replaceArg(*convObs,*shiftObs2) ;
288 
289  pdf1Clone->addOwnedComponents(*shiftObs2) ;
290  pdf1Clone->addOwnedComponents(*clonePdf2) ;
291 
292  pdf2Clone = (RooAbsPdf*) cust.build() ;
293 
294  } else {
295  pdf2Clone = clonePdf2 ;
296  }
297 
298 
299  // Attach cloned pdf to all original parameters of self
300  RooArgSet* fftParams = self.getParameters(*convObs) ;
301 
302  // Remove all cache histogram from fftParams as these
303  // observable need to remain attached to the histogram
304  fftParams->remove(*hist()->get(),kTRUE,kTRUE) ;
305 
306  pdf1Clone->recursiveRedirectServers(*fftParams) ;
307  pdf2Clone->recursiveRedirectServers(*fftParams) ;
308  pdf1Clone->fixAddCoefRange(refName.c_str()) ;
309  pdf2Clone->fixAddCoefRange(refName.c_str()) ;
310 
311  delete fftParams ;
312 
313  // Save copy of original histX binning and make alternate binning
314  // for extended range scanning
315 
316  Int_t N = convObs->numBins() ;
317  Int_t Nbuf = static_cast<Int_t>((N*self.bufferFraction())/2 + 0.5) ;
318  Double_t obw = (convObs->getMax() - convObs->getMin())/N ;
319  Int_t N2 = N+2*Nbuf ;
320 
321  scanBinning = new RooUniformBinning (convObs->getMin()-Nbuf*obw,convObs->getMax()+Nbuf*obw,N2) ;
322  histBinning = convObs->getBinning().clone() ;
323 
324  // Deactivate dirty state propagation on datahist observables
325  // and set all nodes on both pdfs to operMode AlwaysDirty
326  hist()->setDirtyProp(kFALSE) ;
327  convObs->setOperMode(ADirty,kTRUE) ;
328 }
329 
330 
331 ////////////////////////////////////////////////////////////////////////////////
332 /// Suffix for cache histogram (added in addition to suffix for cache name)
333 
335 {
336  return TString(Form("_BufFrac%3.1f_BufStrat%d",_bufFrac,_bufStrat)) ;
337 }
338 
339 
340 
341 ////////////////////////////////////////////////////////////////////////////////
342 /// Returns all RooAbsArg objects contained in the cache element
343 
345 {
347 
348  ret.add(*pdf1Clone) ;
349  ret.add(*pdf2Clone) ;
350  if (pdf1Clone->ownedComponents()) {
351  ret.add(*pdf1Clone->ownedComponents()) ;
352  }
353  if (pdf2Clone->ownedComponents()) {
354  ret.add(*pdf2Clone->ownedComponents()) ;
355  }
356 
357  return ret ;
358 }
359 
360 
361 ////////////////////////////////////////////////////////////////////////////////
362 
364 {
365  delete fftr2c1 ;
366  delete fftr2c2 ;
367  delete fftc2r ;
368 
369  delete pdf1Clone ;
370  delete pdf2Clone ;
371 
372  delete histBinning ;
373  delete scanBinning ;
374 
375 }
376 
377 
378 
379 
380 ////////////////////////////////////////////////////////////////////////////////
381 /// Fill the contents of the cache the FFT convolution output
382 
384 {
385  RooDataHist& cacheHist = *cache.hist() ;
386 
387  ((FFTCacheElem&)cache).pdf1Clone->setOperMode(ADirty,kTRUE) ;
388  ((FFTCacheElem&)cache).pdf2Clone->setOperMode(ADirty,kTRUE) ;
389 
390  // Determine if there other observables than the convolution observable in the cache
391  RooArgSet otherObs ;
392  RooArgSet(*cacheHist.get()).snapshot(otherObs) ;
393 
394  RooAbsArg* histArg = otherObs.find(_x.arg().GetName()) ;
395  if (histArg) {
396  otherObs.remove(*histArg,kTRUE,kTRUE) ;
397  delete histArg ;
398  }
399 
400  //cout << "RooFFTConvPdf::fillCacheObject() otherObs = " << otherObs << endl ;
401 
402  // Handle trivial scenario -- no other observables
403  if (otherObs.getSize()==0) {
405  return ;
406  }
407 
408  // Handle cases where there are other cache slices
409  // Iterator over available slice positions and fill each
410 
411  // Determine number of bins for each slice position observable
412  Int_t n = otherObs.getSize() ;
413  Int_t* binCur = new Int_t[n+1] ;
414  Int_t* binMax = new Int_t[n+1] ;
415  Int_t curObs = 0 ;
416 
417  RooAbsLValue** obsLV = new RooAbsLValue*[n] ;
418  TIterator* iter = otherObs.createIterator() ;
419  RooAbsArg* arg ;
420  Int_t i(0) ;
421  while((arg=(RooAbsArg*)iter->Next())) {
422  RooAbsLValue* lvarg = dynamic_cast<RooAbsLValue*>(arg) ;
423  obsLV[i] = lvarg ;
424  binCur[i] = 0 ;
425  // coverity[FORWARD_NULL]
426  binMax[i] = lvarg->numBins(binningName())-1 ;
427  i++ ;
428  }
429  delete iter ;
430 
431  Bool_t loop(kTRUE) ;
432  while(loop) {
433  // Set current slice position
434  for (Int_t j=0 ; j<n ; j++) { obsLV[j]->setBin(binCur[j],binningName()) ; }
435 
436 // cout << "filling slice: bin of obsLV[0] = " << obsLV[0]->getBin() << endl ;
437 
438  // Fill current slice
439  fillCacheSlice((FFTCacheElem&)cache,otherObs) ;
440 
441  // Determine which iterator to increment
442  while(binCur[curObs]==binMax[curObs]) {
443 
444  // Reset current iterator and consider next iterator ;
445  binCur[curObs]=0 ;
446  curObs++ ;
447 
448  // master termination condition
449  if (curObs==n) {
450  loop=kFALSE ;
451  break ;
452  }
453  }
454 
455  // Increment current iterator
456  binCur[curObs]++ ;
457  curObs=0 ;
458 
459  }
460 
461  delete[] obsLV ;
462  delete[] binMax ;
463  delete[] binCur ;
464 
465 }
466 
467 
468 ////////////////////////////////////////////////////////////////////////////////
469 /// Fill a slice of cachePdf with the output of the FFT convolution calculation
470 
471 void RooFFTConvPdf::fillCacheSlice(FFTCacheElem& aux, const RooArgSet& slicePos) const
472 {
473  // Extract histogram that is the basis of the RooHistPdf
474  RooDataHist& cacheHist = *aux.hist() ;
475 
476  // Sample array of input points from both pdfs
477  // Note that returned arrays have optional buffers zones below and above range ends
478  // to reduce cyclical effects and have been cyclically rotated so that bin containing
479  // zero value is at position zero. Example:
480  //
481  // original: -5 -4 -3 -2 -1 0 +1 +2 +3 +4 +5
482  // add buffer zones: U U -5 -4 -3 -2 -1 0 +1 +2 +3 +4 +5 O O
483  // rotate: 0 +1 +2 +3 +4 +5 O O U U -5 -4 -3 -2 -1
484  //
485  //
486 
487  Int_t N,N2,binShift1,binShift2 ;
488 
489  RooRealVar* histX = (RooRealVar*) cacheHist.get()->find(_x.arg().GetName()) ;
490  if (_bufStrat==Extend) histX->setBinning(*aux.scanBinning) ;
491  Double_t* input1 = scanPdf((RooRealVar&)_x.arg(),*aux.pdf1Clone,cacheHist,slicePos,N,N2,binShift1,_shift1) ;
492  Double_t* input2 = scanPdf((RooRealVar&)_x.arg(),*aux.pdf2Clone,cacheHist,slicePos,N,N2,binShift2,_shift2) ;
493  if (_bufStrat==Extend) histX->setBinning(*aux.histBinning) ;
494 
495 
496 
497 
498  // Retrieve previously defined FFT transformation plans
499  if (!aux.fftr2c1) {
500  aux.fftr2c1 = TVirtualFFT::FFT(1, &N2, "R2CK");
501  aux.fftr2c2 = TVirtualFFT::FFT(1, &N2, "R2CK");
502  aux.fftc2r = TVirtualFFT::FFT(1, &N2, "C2RK");
503  }
504 
505  // Real->Complex FFT Transform on p.d.f. 1 sampling
506  aux.fftr2c1->SetPoints(input1);
507  aux.fftr2c1->Transform();
508 
509  // Real->Complex FFT Transform on p.d.f 2 sampling
510  aux.fftr2c2->SetPoints(input2);
511  aux.fftr2c2->Transform();
512 
513  // Loop over first half +1 of complex output results, multiply
514  // and set as input of reverse transform
515  for (Int_t i=0 ; i<N2/2+1 ; i++) {
516  Double_t re1,re2,im1,im2 ;
517  aux.fftr2c1->GetPointComplex(i,re1,im1) ;
518  aux.fftr2c2->GetPointComplex(i,re2,im2) ;
519  Double_t re = re1*re2 - im1*im2 ;
520  Double_t im = re1*im2 + re2*im1 ;
521  TComplex t(re,im) ;
522  aux.fftc2r->SetPointComplex(i,t) ;
523  }
524 
525  // Reverse Complex->Real FFT transform product
526  aux.fftc2r->Transform() ;
527 
528  Int_t totalShift = binShift1 + (N2-N)/2 ;
529 
530  // Store FFT result in cache
531 
532  TIterator* iter = const_cast<RooDataHist&>(cacheHist).sliceIterator(const_cast<RooAbsReal&>(_x.arg()),slicePos) ;
533  for (Int_t i =0 ; i<N ; i++) {
534 
535  // Cyclically shift array back so that bin containing zero is back in zeroBin
536  Int_t j = i + totalShift ;
537  while (j<0) j+= N2 ;
538  while (j>=N2) j-= N2 ;
539 
540  iter->Next() ;
541  cacheHist.set(aux.fftc2r->GetPointReal(j)) ;
542  }
543  delete iter ;
544 
545  // cacheHist.dump2() ;
546 
547  // Delete input arrays
548  delete[] input1 ;
549  delete[] input2 ;
550 
551 }
552 
553 
554 ////////////////////////////////////////////////////////////////////////////////
555 /// Scan the values of 'pdf' in observable 'obs' using the bin values stored in 'hist' at slice position 'slicePos'
556 /// N is filled with the number of bins defined in hist, N2 is filled with N plus the number of buffer bins
557 /// The return value is an array of doubles of length N2 with the sampled values. The caller takes ownership
558 /// of the array
559 
560 Double_t* RooFFTConvPdf::scanPdf(RooRealVar& obs, RooAbsPdf& pdf, const RooDataHist& hist, const RooArgSet& slicePos,
561  Int_t& N, Int_t& N2, Int_t& zeroBin, Double_t shift) const
562 {
563 
564  RooRealVar* histX = (RooRealVar*) hist.get()->find(obs.GetName()) ;
565 
566  // Calculate number of buffer bins on each size to avoid cyclical flow
567  N = histX->numBins(binningName()) ;
568  Int_t Nbuf = static_cast<Int_t>((N*bufferFraction())/2 + 0.5) ;
569  N2 = N+2*Nbuf ;
570 
571 
572  // Allocate array of sampling size plus optional buffer zones
573  Double_t* array = new Double_t[N2] ;
574 
575  // Set position of non-convolution observable to that of the cache slice that were are processing now
576  hist.get(slicePos) ;
577 
578  // Find bin ID that contains zero value
579  zeroBin = 0 ;
580  if (histX->getMax()>=0 && histX->getMin()<=0) {
581  zeroBin = histX->getBinning().binNumber(0) ;
582  } else if (histX->getMin()>0) {
583  Double_t bw = (histX->getMax() - histX->getMin())/N2 ;
584  zeroBin = Int_t(-histX->getMin()/bw) ;
585  } else {
586  Double_t bw = (histX->getMax() - histX->getMin())/N2 ;
587  zeroBin = Int_t(-1*histX->getMax()/bw) ;
588  }
589 
590  Int_t binShift = Int_t((N2* shift) / (histX->getMax()-histX->getMin())) ;
591 
592  zeroBin += binShift ;
593  while(zeroBin>=N2) zeroBin-= N2 ;
594  while(zeroBin<0) zeroBin+= N2 ;
595 
596  // First scan hist into temp array
597  Double_t *tmp = new Double_t[N2] ;
598  Int_t k(0) ;
599  switch(_bufStrat) {
600 
601  case Extend:
602  // Sample entire extended range (N2 samples)
603  for (k=0 ; k<N2 ; k++) {
604  histX->setBin(k) ;
605  tmp[k] = pdf.getVal(hist.get()) ;
606  }
607  break ;
608 
609  case Flat:
610  // Sample original range (N samples) and fill lower and upper buffer
611  // bins with p.d.f. value at respective boundary
612  {
613  histX->setBin(0) ;
614  Double_t val = pdf.getVal(hist.get()) ;
615  for (k=0 ; k<Nbuf ; k++) {
616  tmp[k] = val ;
617  }
618  for (k=0 ; k<N ; k++) {
619  histX->setBin(k) ;
620  tmp[k+Nbuf] = pdf.getVal(hist.get()) ;
621  }
622  histX->setBin(N-1) ;
623  val = pdf.getVal(hist.get()) ;
624  for (k=0 ; k<Nbuf ; k++) {
625  tmp[N+Nbuf+k] = val ;
626  }
627  }
628  break ;
629 
630  case Mirror:
631  // Sample original range (N samples) and fill lower and upper buffer
632  // bins with mirror image of sampled range
633  for (k=0 ; k<N ; k++) {
634  histX->setBin(k) ;
635  tmp[k+Nbuf] = pdf.getVal(hist.get()) ;
636  }
637  for (k=1 ; k<=Nbuf ; k++) {
638  histX->setBin(k) ;
639  tmp[Nbuf-k] = pdf.getVal(hist.get()) ;
640  histX->setBin(N-k) ;
641  tmp[Nbuf+N+k-1] = pdf.getVal(hist.get()) ;
642  }
643  break ;
644  }
645 
646  // Scan function and store values in array
647  for (Int_t i=0 ; i<N2 ; i++) {
648  // Cyclically shift writing location by zero bin position
649  Int_t j = i - (zeroBin) ;
650  if (j<0) j+= N2 ;
651  if (j>=N2) j-= N2 ;
652  array[i] = tmp[j] ;
653  }
654 
655  // Cleanup
656  delete[] tmp ;
657  return array ;
658 }
659 
660 
661 
662 ////////////////////////////////////////////////////////////////////////////////
663 /// Return the observables to be cached given the normalization set nset
664 ///
665 /// If the cache observables is in nset then this is
666 /// - the convolution observable plus
667 /// - any member of nset that is either a RooCategory,
668 /// - or was previously specified through setCacheObservables().
669 ///
670 /// In case the cache observable is _not_ in nset, then it is
671 /// - the convolution observable plus
672 /// - all member of nset are observables of this p.d.f.
673 ///
674 
676 {
677  // Get complete list of observables
678  RooArgSet* obs1 = _pdf1.arg().getObservables(nset) ;
679  RooArgSet* obs2 = _pdf2.arg().getObservables(nset) ;
680  obs1->add(*obs2,kTRUE) ;
681 
682  // Check if convolution observable is in nset
683  if (nset.contains(_x.arg())) {
684 
685  // Now strip out all non-category observables
686  TIterator* iter = obs1->createIterator() ;
687  RooAbsArg* arg ;
688  RooArgSet killList ;
689  while((arg=(RooAbsArg*)iter->Next())) {
690  if (arg->IsA()->InheritsFrom(RooAbsReal::Class()) && !_cacheObs.find(arg->GetName())) {
691  killList.add(*arg) ;
692  }
693  }
694  delete iter ;
695  obs1->remove(killList) ;
696 
697  // And add back the convolution observables
698  obs1->add(_x.arg(),kTRUE) ;
699 
700  obs1->add(_cacheObs) ;
701 
702  delete obs2 ;
703 
704  } else {
705 
706  // If cacheObs was filled, cache only observables in there
707  if (_cacheObs.getSize()>0) {
708  TIterator* iter = obs1->createIterator() ;
709  RooAbsArg* arg ;
710  RooArgSet killList ;
711  while((arg=(RooAbsArg*)iter->Next())) {
712  if (arg->IsA()->InheritsFrom(RooAbsReal::Class()) && !_cacheObs.find(arg->GetName())) {
713  killList.add(*arg) ;
714  }
715  }
716  delete iter ;
717  obs1->remove(killList) ;
718  }
719 
720 
721  // Make sure convolution observable is always in there
722  obs1->add(_x.arg(),kTRUE) ;
723  delete obs2 ;
724 
725  }
726 
727  return obs1 ;
728 }
729 
730 
731 
732 ////////////////////////////////////////////////////////////////////////////////
733 /// Return the parameters on which the cache depends given normalization
734 /// set nset. For this p.d.f these are the parameters of the input p.d.f.
735 /// but never the convolution variable, it case it is not part of nset
736 
738 {
739  RooArgSet* vars = getVariables() ;
740  RooArgSet* obs = actualObservables(nset) ;
741  vars->remove(*obs) ;
742  delete obs ;
743 
744  return vars ;
745 }
746 
747 
748 
749 ////////////////////////////////////////////////////////////////////////////////
750 /// Return p.d.f. observable (which can be a function) to substitute given
751 /// p.d.f. observable. Substitute x by xprime if xprime is set
752 
754 {
755  if (_xprime.absArg() && string(histObservable.GetName())==_x.absArg()->GetName()) {
756  return (*_xprime.absArg()) ;
757  }
758  return histObservable ;
759 }
760 
761 
762 
763 ////////////////////////////////////////////////////////////////////////////////
764 /// Create appropriate generator context for this convolution. If both input p.d.f.s support
765 /// internal generation, if it is safe to use them and if no observables other than the convolution
766 /// observable are requested for generation, use the specialized convolution generator context
767 /// which implements a smearing strategy in the convolution observable. If not return the
768 /// regular accept/reject generator context
769 
771  const RooArgSet* auxProto, Bool_t verbose) const
772 {
773  RooArgSet vars2(vars) ;
774  vars2.remove(_x.arg(),kTRUE,kTRUE) ;
775  Int_t numAddDep = vars2.getSize() ;
776 
777  RooArgSet dummy ;
778  Bool_t pdfCanDir = (((RooAbsPdf&)_pdf1.arg()).getGenerator(_x.arg(),dummy) != 0 && \
780  Bool_t resCanDir = (((RooAbsPdf&)_pdf2.arg()).getGenerator(_x.arg(),dummy) !=0 &&
782 
783  if (pdfCanDir) {
784  cxcoutI(Generation) << "RooFFTConvPdf::genContext() input p.d.f " << _pdf1.arg().GetName()
785  << " has internal generator that is safe to use in current context" << endl ;
786  }
787  if (resCanDir) {
788  cxcoutI(Generation) << "RooFFTConvPdf::genContext() input p.d.f. " << _pdf2.arg().GetName()
789  << " has internal generator that is safe to use in current context" << endl ;
790  }
791  if (numAddDep>0) {
792  cxcoutI(Generation) << "RooFFTConvPdf::genContext() generation requested for observables other than the convolution observable " << _x.arg().GetName() << endl ;
793  }
794 
795 
796  if (numAddDep>0 || !pdfCanDir || !resCanDir) {
797  // Any resolution model with more dependents than the convolution variable
798  // or pdf or resmodel do not support direct generation
799  cxcoutI(Generation) << "RooFFTConvPdf::genContext() selecting accept/reject generator context because one or both of the input "
800  << "p.d.f.s cannot use internal generator and/or "
801  << "observables other than the convolution variable are requested for generation" << endl ;
802  return new RooGenContext(*this,vars,prototype,auxProto,verbose) ;
803  }
804 
805  // Any other resolution model: use specialized generator context
806  cxcoutI(Generation) << "RooFFTConvPdf::genContext() selecting specialized convolution generator context as both input "
807  << "p.d.fs are safe for internal generator and only "
808  << "the convolution observables is requested for generation" << endl ;
809  return new RooConvGenContext(*this,vars,prototype,auxProto,verbose) ;
810 }
811 
812 
813 
814 ////////////////////////////////////////////////////////////////////////////////
815 /// Change the size of the buffer on either side of the observable range to frac times the
816 /// size of the range of the convolution observable
817 
819 {
820  if (frac<0) {
821  coutE(InputArguments) << "RooFFTConvPdf::setBufferFraction(" << GetName() << ") fraction should be greater than or equal to zero" << endl ;
822  return ;
823  }
824  _bufFrac = frac ;
825 
826  // Sterilize the cache as certain partial results depend on buffer fraction
827  _cacheMgr.sterilize() ;
828 }
829 
830 
831 ////////////////////////////////////////////////////////////////////////////////
832 /// Change strategy to fill the overflow buffer on either side of the convolution observable range.
833 ///
834 /// 'Extend' means is that the input p.d.f convolution observable range is widened to include the buffer range
835 /// 'Flat' means that the buffer is filled with the p.d.f. value at the boundary of the observable range
836 /// 'Mirror' means that the buffer is filled with a ,irror image of the p.d.f. around the convolution observable boundary
837 ///
838 /// The default strategy is extend. If one of the input p.d.f.s is a RooAddPdf, it is configured so that the interpretation
839 /// range of the fraction coefficients is kept at the nominal convolutions observable range (instead of interpreting coefficients
840 /// in the widened range including the buffer)
841 
843 {
844  _bufStrat = bs ;
845 }
846 
847 
848 
849 ////////////////////////////////////////////////////////////////////////////////
850 /// Customized printing of arguments of a RooNumConvPdf to more intuitively reflect the contents of the
851 /// product operator construction
852 
853 void RooFFTConvPdf::printMetaArgs(ostream& os) const
854 {
855  os << _pdf1.arg().GetName() << "(" << _x.arg().GetName() << ") (*) " << _pdf2.arg().GetName() << "(" << _x.arg().GetName() << ") " ;
856 }
857 
858 
859 
860 ////////////////////////////////////////////////////////////////////////////////
861 /// (Re)calculate effective parameters of this p.d.f.
862 
864 {
865  RooArgSet* params1 = _pdf1.arg().getParameters(_x.arg()) ;
866  RooArgSet* params2 = _pdf2.arg().getParameters(_x.arg()) ;
867  _params.removeAll() ;
868  _params.add(*params1) ;
869  _params.add(*params2,kTRUE) ;
870  delete params1 ;
871  delete params2 ;
872 }
873 
874 
875 
876 ////////////////////////////////////////////////////////////////////////////////
877 ///calcParams() ;
878 
879 Bool_t RooFFTConvPdf::redirectServersHook(const RooAbsCollection& /*newServerList*/, Bool_t /*mustReplaceAll*/, Bool_t /*nameChange*/, Bool_t /*isRecursive*/)
880 {
881  return kFALSE ;
882 }
virtual Double_t getMin(const char *name=0) const
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:2073
void calcParams()
(Re)calculate effective parameters of this p.d.f.
RooRealProxy _pdf1
Definition: RooFFTConvPdf.h:69
friend class RooConvGenContext
RooSetProxy _cacheObs
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual void setBin(Int_t ibin, const char *rangeName=0)
Set value to center of bin &#39;ibin&#39; of binning &#39;rangeName&#39; (or of default binning if no range is specif...
virtual RooAbsArg * cloneTree(const char *newname=0) const
Clone tree expression of objects.
Definition: RooAbsArg.cxx:2287
RooAbsBinning * scanBinning
Definition: RooFFTConvPdf.h:93
TIterator * createIterator(Bool_t dir=kIterForward) const
Bool_t hasBinning(const char *name) const
Returns true if variable has a binning with &#39;name&#39;.
Definition: RooRealVar.cxx:253
virtual void fillCacheObject(PdfCacheElem &cache) const
Fill the contents of the cache the FFT convolution output.
Double_t _shift2
#define coutE(a)
Definition: RooMsgService.h:34
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
virtual RooArgSet * actualObservables(const RooArgSet &nset) const
Return the observables to be cached given the normalization set nset.
virtual Double_t getMax(const char *name=0) const
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
virtual const RooArgSet * get() const
Definition: RooDataHist.h:77
BufStrat _bufStrat
#define cxcoutI(a)
Definition: RooMsgService.h:83
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Definition: RooAbsArg.h:194
virtual Int_t numBins(const char *rangeName=0) const
virtual void removeAll()
Remove all argument inset using remove(const RooAbsArg&).
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
const RooArgSet * ownedComponents() const
Definition: RooAbsArg.h:414
RooRealProxy _pdf2
Definition: RooFFTConvPdf.h:70
const RooAbsBinning & getBinning(const char *name=0, Bool_t verbose=kTRUE, Bool_t createOnTheFly=kFALSE) const
Return binning definition with name.
Definition: RooRealVar.cxx:266
void set(Double_t weight, Double_t wgtErr=-1)
Increment the weight of the bin enclosing the coordinates given by &#39;row&#39; by the specified amount...
#define N
void setBufferStrategy(BufStrat bs)
Change strategy to fill the overflow buffer on either side of the convolution observable range...
virtual RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=0, const RooArgSet *auxProto=0, Bool_t verbose=kFALSE) const
Create appropriate generator context for this convolution.
void printMetaArgs(std::ostream &os) const
Customized printing of arguments of a RooNumConvPdf to more intuitively reflect the contents of the p...
virtual RooArgList containedArgs(Action)
Returns all RooAbsArg objects contained in the cache element.
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
FFTCacheElem(const RooFFTConvPdf &self, const RooArgSet *nset)
Clone input pdf and attach to dataset.
virtual PdfCacheElem * createCache(const RooArgSet *nset) const
Return specialized cache subclass for FFT calculations.
virtual void fixAddCoefRange(const char *rangeName=0, Bool_t force=kTRUE)
Fix the interpretation of the coefficient of any RooAddPdf component in the expression tree headed by...
Bool_t addOwnedComponents(const RooArgSet &comps)
Take ownership of the contents of &#39;comps&#39;.
Definition: RooAbsArg.cxx:2273
STL namespace.
virtual RooAbsBinning * clone(const char *name=0) const =0
void attachDataSet(const RooAbsData &set)
Replace server nodes with names matching the dataset variable names with those data set variables...
Definition: RooAbsArg.cxx:1488
Iterator abstract base class.
Definition: TIterator.h:30
Double_t _bufFrac
void setBinning(const RooAbsBinning &binning, const char *name=0)
Add given binning under name &#39;name&#39; with this variable.
Definition: RooRealVar.cxx:345
RooDataSet is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
void Class()
Definition: Class.C:29
virtual TString histNameSuffix() const
Suffix for cache histogram (added in addition to suffix for cache name)
friend class RooArgSet
Definition: RooAbsArg.h:471
virtual RooArgList containedArgs(Action)
Returns all RooAbsArg objects contained in the cache element.
RooAbsGenContext is the abstract base class for generator contexts of RooAbsPdf objects.
friend class FFTCacheElem
Definition: RooFFTConvPdf.h:97
virtual void SetPointComplex(Int_t ipoint, TComplex &c)=0
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
RooRealProxy _x
Definition: RooFFTConvPdf.h:67
virtual void setBin(Int_t ibin, const char *rangeName=0)=0
virtual const char * binningName() const
Double_t bufferFraction() const
Definition: RooFFTConvPdf.h:41
Int_t getSize() const
Double_t * scanPdf(RooRealVar &obs, RooAbsPdf &pdf, const RooDataHist &hist, const RooArgSet &slicePos, Int_t &N, Int_t &N2, Int_t &zeroBin, Double_t shift) const
Scan the values of &#39;pdf&#39; in observable &#39;obs&#39; using the bin values stored in &#39;hist&#39; at slice position ...
void fillCacheSlice(FFTCacheElem &cache, const RooArgSet &slicePosition) const
Fill a slice of cachePdf with the output of the FFT convolution calculation.
virtual Double_t GetPointReal(Int_t ipoint, Bool_t fromInput=kFALSE) const =0
auto * a
Definition: textangle.C:12
static TVirtualFFT * FFT(Int_t ndim, Int_t *n, Option_t *option)
Returns a pointer to the FFT of requested size and type.
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:443
char * Form(const char *fmt,...)
virtual void Transform()=0
RooAbsArg * absArg() const
Definition: RooArgProxy.h:37
virtual RooAbsArg & pdfObservable(RooAbsArg &histObservable) const
Return p.d.f.
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
Class RooGenContext implement a universal generator context for all RooAbsPdf classes that do not hav...
Definition: RooGenContext.h:30
virtual Int_t binNumber(Double_t x) const =0
virtual void SetPoints(const Double_t *data)=0
virtual Bool_t isDirectGenSafe(const RooAbsArg &arg) const
Check if given observable can be safely generated using the pdfs internal generator mechanism (if tha...
Definition: RooAbsPdf.cxx:2159
const Bool_t kFALSE
Definition: RtypesCore.h:88
void replaceArg(const RooAbsArg &orig, const RooAbsArg &subst)
Replace any occurence of arg &#39;orig&#39; with arg &#39;subst&#39;.
virtual void GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput=kFALSE) const =0
RooUniformBinning is an implementation of RooAbsBinning that provides a uniform binning in &#39;n&#39; bins b...
RooAbsCachedPdf is the abstract base class for p.d.f.s that need or want to cache their evaluate() ou...
#define ClassImp(name)
Definition: Rtypes.h:359
RooAbsArg * find(const char *name) const
Find object with given name in list.
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
virtual Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
Definition: RooAbsPdf.cxx:2124
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:532
static RooMathCoreReg dummy
virtual RooArgSet * actualParameters(const RooArgSet &nset) const
Return the parameters on which the cache depends given normalization set nset.
RooSetProxy _params
Definition: RooFFTConvPdf.h:71
Double_t _shift1
void setRange(const char *name, Double_t min, Double_t max)
Set range named &#39;name to [min,max].
Definition: RooRealVar.cxx:448
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects...
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsBinning * histBinning
Definition: RooFFTConvPdf.h:92
RooConstVar & RooConst(Double_t val)
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
void setBufferFraction(Double_t frac)
Change the size of the buffer on either side of the observable range to frac times the size of the ra...
RooAbsArg * build(const char *masterCatState, Bool_t verbose=kFALSE)
Build a clone of the prototype executing all registered &#39;replace&#39; rules and &#39;split&#39; rules for the mas...
virtual TObject * Next()=0
virtual const char * inputBaseName() const
Return base name component for cache components in this case &#39;PDF1_CONV_PDF2&#39;.
void setOperMode(OperMode mode, Bool_t recurseADirty=kTRUE)
Change cache operation mode to given mode.
Definition: RooAbsArg.cxx:1745
virtual ~RooFFTConvPdf()
Destructor.
Abstract base class for objects that are lvalues, i.e.
Definition: RooAbsLValue.h:26
RooObjCacheManager _cacheMgr
Bool_t contains(const RooAbsArg &var) const
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
const RooAbsReal & arg() const
Definition: RooRealProxy.h:43
Bool_t recursiveRedirectServers(const RooAbsCollection &newServerList, Bool_t mustReplaceAll=kFALSE, Bool_t nameChange=kFALSE, Bool_t recurseInNewSet=kTRUE)
Definition: RooAbsArg.cxx:1079
Bool_t redirectServersHook(const RooAbsCollection &newServerList, Bool_t mustReplaceAll, Bool_t nameChange, Bool_t isRecursive)
calcParams() ;
RooRealProxy _xprime
Definition: RooFFTConvPdf.h:68
const Bool_t kTRUE
Definition: RtypesCore.h:87
const Int_t n
Definition: legend1.C:16
char name[80]
Definition: TGX11.cxx:109
virtual Int_t numBins(const char *rangeName=0) const =0
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Overloaded RooArgSet::add() method inserts &#39;var&#39; into set and registers &#39;var&#39; as server to owner with...