include/root/RooFFTConvPdf.h

0001 /*****************************************************************************
0002  * Project: RooFit                                                           *
0003  *                                                                           *
0004  * Copyright (c) 2000-2005, Regents of the University of California          *
0005  *                          and Stanford University. All rights reserved.    *
0006  *                                                                           *
0007  * Redistribution and use in source and binary forms,                        *
0008  * with or without modification, are permitted according to the terms        *
0009  * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
0010  *****************************************************************************/
0011
0012 #ifndef ROOFFTCONVPDF
0013 #define ROOFFTCONVPDF
0014
0015 #include "RooAbsCachedPdf.h"
0016 #include "RooRealProxy.h"
0017 #include "RooSetProxy.h"
0018 #include "RooAbsReal.h"
0019 #include "RooHistPdf.h"
0020 #include "TVirtualFFT.h"
0021
0022 class RooRealVar;
0023
0024 ///PDF for the numerical (FFT) convolution of two PDFs.
0025 class RooFFTConvPdf : public RooAbsCachedPdf {
0026 public:
0027
0028   RooFFTConvPdf() {
0029     // coverity[UNINIT_CTOR]
0030   } ;
0031   RooFFTConvPdf(const char *name, const char *title, RooRealVar& convVar, RooAbsPdf& pdf1, RooAbsPdf& pdf2, Int_t ipOrder=2);
0032   RooFFTConvPdf(const char *name, const char *title, RooAbsReal& pdfConvVar, RooRealVar& convVar, RooAbsPdf& pdf1, RooAbsPdf& pdf2, Int_t ipOrder=2);
0033   RooFFTConvPdf(const RooFFTConvPdf& other, const char* name=nullptr) ;
0034   TObject* clone(const char* newname=nullptr) const override { return new RooFFTConvPdf(*this,newname); }
0035   ~RooFFTConvPdf() override ;
0036
0037   void setShift(double val1, double val2) { _shift1 = val1 ; _shift2 = val2 ; }
0038   void setCacheObservables(const RooArgSet& obs) { _cacheObs.removeAll() ; _cacheObs.add(obs) ; }
0039   const RooArgSet& cacheObservables() const { return _cacheObs ; }
0040
0041   /// Return value of buffer fraction applied in FFT calculation array beyond either
0042   /// end of the observable domain to reduce cyclical effects
0043   double bufferFraction() const {
0044     return _bufFrac ;
0045   }
0046
0047   enum BufStrat { Extend=0, Mirror=1, Flat=2 } ;
0048   /// Return the strategy currently used to fill the buffer:
0049   /// 'Extend' means is that the input p.d.f convolution observable range is widened to include the buffer range
0050   /// 'Flat' means that the buffer is filled with the p.d.f. value at the boundary of the observable range
0051   /// 'Mirror' means that the buffer is filled with a mirror image of the p.d.f. around the convolution observable boundary
0052   BufStrat bufferStrategy() const {
0053     return _bufStrat ;
0054   }
0055   void setBufferStrategy(BufStrat bs) ;
0056   void setBufferFraction(double frac) ;
0057
0058   void printMetaArgs(std::ostream& os) const override ;
0059
0060   // Propagate maximum value estimate of pdf1 as convolution can only result in lower max values
0061   Int_t getMaxVal(const RooArgSet& vars) const override { return _pdf1.arg().getMaxVal(vars) ; }
0062   double maxVal(Int_t code) const override { return _pdf1.arg().maxVal(code) ; }
0063
0064   RooAbsReal const &getConvVar() const { return *_x; }
0065
0066   RooAbsReal const &getPdf1() const { return *_pdf1; }
0067
0068   RooAbsReal const &getPdf2() const { return *_pdf2; }
0069
0070   RooAbsReal const *getPdfConvVar() const { return dynamic_cast<RooAbsReal *>(_xprime.absArg()); }
0071
0072   protected:
0073
0074   RooRealProxy _x ;      ///< Convolution observable
0075   RooRealProxy _xprime ; ///< Input function representing value of convolution observable
0076   RooRealProxy _pdf1 ;   ///< First input p.d.f
0077   RooRealProxy _pdf2 ;   ///< Second input p.d.f
0078   RooSetProxy _params ;  ///< Effective parameters of this p.d.f.
0079
0080   void calcParams() ;
0081
0082   std::vector<double>  scanPdf(RooRealVar& obs, RooAbsPdf& pdf, const RooDataHist& hist, const RooArgSet& slicePos, Int_t& N, Int_t& N2, Int_t& zeroBin, double shift) const ;
0083
0084   class FFTCacheElem : public PdfCacheElem {
0085   public:
0086     FFTCacheElem(const RooFFTConvPdf& self, const RooArgSet* nset) ;
0087
0088     RooArgList containedArgs(Action) override ;
0089
0090     std::unique_ptr<TVirtualFFT> fftr2c1;
0091     std::unique_ptr<TVirtualFFT> fftr2c2;
0092     std::unique_ptr<TVirtualFFT> fftc2r;
0093
0094     std::unique_ptr<RooAbsPdf> pdf1Clone;
0095     std::unique_ptr<RooAbsPdf> pdf2Clone;
0096
0097     std::unique_ptr<RooAbsBinning> histBinning;
0098     std::unique_ptr<RooAbsBinning> scanBinning;
0099   };
0100
0101   friend class FFTCacheElem ;
0102
0103   double evaluate() const override { RooArgSet dummy(_x.arg()) ; return getVal(&dummy) ; } ; // dummy
0104   const char* inputBaseName() const override ;
0105   RooFit::OwningPtr<RooArgSet> actualObservables(const RooArgSet& nset) const override ;
0106   RooFit::OwningPtr<RooArgSet> actualParameters(const RooArgSet& nset) const override ;
0107   RooAbsArg& pdfObservable(RooAbsArg& histObservable) const override ;
0108   void fillCacheObject(PdfCacheElem& cache) const override ;
0109   void fillCacheSlice(FFTCacheElem& cache, const RooArgSet& slicePosition) const ;
0110
0111   PdfCacheElem* createCache(const RooArgSet* nset) const override ;
0112   TString histNameSuffix() const override ;
0113
0114   // mutable std::map<const RooHistPdf*,CacheAuxInfo*> _cacheAuxInfo ; //! Auxiliary Cache information (do not persist)
0115   double _bufFrac ; // Sampling buffer size as fraction of domain size
0116   BufStrat _bufStrat ; // Strategy to fill the buffer
0117
0118   double  _shift1 ;
0119   double  _shift2 ;
0120
0121   RooAbsGenContext* genContext(const RooArgSet &vars, const RooDataSet *prototype=nullptr,
0122                                        const RooArgSet* auxProto=nullptr, bool verbose= false) const override ;
0123
0124   friend class RooConvGenContext ;
0125   RooSetProxy  _cacheObs ; ///< Non-convolution observables that are also cached
0126
0127 private:
0128
0129   void prepareFFTBinning(RooRealVar& convVar) const;
0130
0131   ClassDefOverride(RooFFTConvPdf,1) // Convolution operator p.d.f based on numeric Fourier transforms
0132 };
0133
0134 #endif