// Numeric 1-dimensional convolution operator PDF. This class can convolve any PDF
// with any other PDF
// <p>
// This class should not be used blindly as numeric convolution is computing
// intensive and prone to stability fitting problems. If an analytic convolution
// can be calculated, you should use that or implement it if not available.
// RooNumConvolution implements reasonable defaults that should convolve most
// functions reasonably well, but results strongly depend on the shape of your
// input PDFS so always check your result.
//
// The default integration engine for the numeric convolution is the
// adaptive Gauss-Kronrod method, which empirically seems the most robust
// for this task. You can override the convolution integration settings via
// the RooNumIntConfig object reference returned by the convIntConfig() member
// function
// <p>
// By default the numeric convolution is integrated from -infinity to
// +infinity through a <pre>x -> 1/x</pre> coordinate transformation of the
// tails. For convolution with a very small bandwidth it may be
// advantageous (for both CPU consumption and stability) if the
// integration domain is limited to a finite range. The function
// setConvolutionWindow(mean,width,scale) allows to set a sliding
// window around the x value to be calculated taking a RooAbsReal
// expression for an offset and a width to be taken around the x
// value. These input expression can be RooFormulaVars or other
// function objects although the 3d 'scale' argument 'scale'
// multiplies the width RooAbsReal expression given in the 2nd
// argument, allowing for an appropriate window definition for most
// cases without need for a RooFormulaVar object: e.g. a Gaussian
// resolution PDF do setConvolutionWindow(gaussMean,gaussSigma,5)
// Note that for a 'wide' Gaussian the -inf to +inf integration
// may converge more quickly than that over a finite range!
// <p>
// The default numeric precision is 1e-7, i.e. the global default for
// numeric integration but you should experiment with this value to
// see if it is sufficient for example by studying the number of function
// calls that MINUIT needs to fit your function as function of the
// convolution precision.
// END_HTML
#include "RooFit.h"
#include "Riostream.h"
#include "Riostream.h"
#include "TH2F.h"
#include "RooNumConvolution.h"
#include "RooArgList.h"
#include "RooRealVar.h"
#include "RooFormulaVar.h"
#include "RooCustomizer.h"
#include "RooConvIntegrandBinding.h"
#include "RooNumIntFactory.h"
#include "RooGenContext.h"
#include "RooConvGenContext.h"
#include "RooMsgService.h"
ClassImp(RooNumConvolution)
;
RooNumConvolution::RooNumConvolution() :
_init(kFALSE),
_integrand(0),
_integrator(0),
_callHist(0)
{
}
RooNumConvolution::RooNumConvolution(const char *name, const char *title, RooRealVar& convVar, RooAbsReal& inPdf, RooAbsReal& resmodel, const RooNumConvolution* proto) :
RooAbsReal(name,title),
_init(kFALSE),
_convIntConfig(RooNumIntConfig::defaultConfig()),
_integrand(0),
_integrator(0),
_origVar("origVar","Original Convolution variable",this,convVar),
_origPdf("origPdf","Original Input PDF",this,inPdf),
_origModel("origModel","Original Resolution model",this,resmodel),
_ownedClonedPdfSet("ownedClonePdfSet"),
_ownedClonedModelSet("ownedCloneModelSet"),
_cloneVar(0),
_clonePdf(0),
_cloneModel(0),
_useWindow(kFALSE),
_windowScale(1),
_windowParam("windowParam","Convolution window parameter",this,kFALSE),
_verboseThresh(2000),
_doProf(kFALSE),
_callHist(0)
{
_convIntConfig.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
_convIntConfig.method1DOpen().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
if (proto) {
convIntConfig() = proto->convIntConfig() ;
if (proto->_useWindow) {
setConvolutionWindow((RooAbsReal&)*proto->_windowParam.at(0),(RooAbsReal&)*proto->_windowParam.at(1),proto->_windowScale) ;
}
}
}
RooNumConvolution::RooNumConvolution(const RooNumConvolution& other, const char* name) :
RooAbsReal(other,name),
_init(kFALSE),
_convIntConfig(other._convIntConfig),
_integrand(0),
_integrator(0),
_origVar("origVar",this,other._origVar),
_origPdf("origPdf",this,other._origPdf),
_origModel("origModel",this,other._origModel),
_ownedClonedPdfSet("ownedClonePdfSet"),
_ownedClonedModelSet("ownedCloneModelSet"),
_cloneVar(0),
_clonePdf(0),
_cloneModel(0),
_useWindow(other._useWindow),
_windowScale(other._windowScale),
_windowParam("windowParam",this,other._windowParam),
_verboseThresh(other._verboseThresh),
_doProf(other._doProf),
_callHist(other._callHist)
{
}
void RooNumConvolution::initialize() const
{
_ownedClonedPdfSet.removeAll() ;
_ownedClonedModelSet.removeAll() ;
if (_cloneVar) delete _cloneVar ;
_cloneVar = new RooRealVar(Form("%s_prime",_origVar.arg().GetName()),"Convolution Variable",0) ;
RooCustomizer mgr1(pdf(),"NumConv_PdfClone") ;
mgr1.setCloneBranchSet(_ownedClonedPdfSet) ;
mgr1.replaceArg(var(),*_cloneVar) ;
_clonePdf = (RooAbsReal*) mgr1.build() ;
RooCustomizer mgr2(model(),"NumConv_ModelClone") ;
mgr2.setCloneBranchSet(_ownedClonedModelSet) ;
mgr2.replaceArg(var(),*_cloneVar) ;
_cloneModel = (RooAbsReal*) mgr2.build() ;
_cloneVar->SetName(var().GetName()) ;
_integrand = new RooConvIntegrandBinding(*_clonePdf,*_cloneModel,*_cloneVar,var(),0) ;
_integrator = RooNumIntFactory::instance().createIntegrator(*_integrand,_convIntConfig,1) ;
_integrator->setUseIntegrandLimits(kFALSE) ;
_init = kTRUE ;
}
RooNumConvolution::~RooNumConvolution()
{
}
Double_t RooNumConvolution::evaluate() const
{
if (!_init) initialize() ;
Double_t x = _origVar ;
_integrand->setNormalizationSet(_origVar.nset()) ;
if (_useWindow) {
Double_t center = ((RooAbsReal*)_windowParam.at(0))->getVal() ;
Double_t width = _windowScale * ((RooAbsReal*)_windowParam.at(1))->getVal() ;
_integrator->setLimits(x-center-width,x-center+width) ;
} else {
_integrator->setLimits(-RooNumber::infinity(),RooNumber::infinity()) ;
}
if (_doProf) _integrand->resetNumCall() ;
Double_t ret = _integrator->integral(&x) ;
if (_doProf) {
_callHist->Fill(x,_integrand->numCall()) ;
if (_integrand->numCall()>_verboseThresh) {
coutW(Integration) << "RooNumConvolution::eveluate(" << GetName() << ") WARNING convolution integral at x=" << x
<< " required " << _integrand->numCall() << " function evaluations" << endl ;
}
}
return ret ;
}
Bool_t RooNumConvolution::redirectServersHook(const RooAbsCollection& , Bool_t ,
Bool_t , Bool_t )
{
_init = kFALSE ;
return kFALSE ;
}
void RooNumConvolution::clearConvolutionWindow()
{
_useWindow = kFALSE ;
_windowParam.removeAll() ;
}
void RooNumConvolution::setConvolutionWindow(RooAbsReal& centerParam, RooAbsReal& widthParam, Double_t widthScaleFactor)
{
_useWindow = kTRUE ;
_windowParam.removeAll() ;
_windowParam.add(centerParam) ;
_windowParam.add(widthParam) ;
_windowScale = widthScaleFactor ;
}
void RooNumConvolution::setCallWarning(Int_t threshold)
{
if (threshold<0) {
coutE(InputArguments) << "RooNumConvolution::setCallWarning(" << GetName() << ") ERROR: threshold must be positive, value unchanged" << endl ;
return ;
}
_verboseThresh = threshold ;
}
void RooNumConvolution::setCallProfiling(Bool_t flag, Int_t nbinX, Int_t nbinCall, Int_t nCallHigh)
{
if (flag) {
if (_doProf) {
delete _callHist ;
}
_callHist = new TH2F(Form("callHist_%s",GetName()),Form("Call Profiling of RooNumConvolution %s",GetTitle()),
nbinX,_origVar.min(),_origVar.max(),
nbinCall,0,nCallHigh) ;
_doProf=kTRUE ;
} else if (_doProf) {
delete _callHist ;
_callHist = 0 ;
_doProf = kFALSE ;
}
}
void RooNumConvolution::printCompactTreeHook(ostream& os, const char* indent)
{
os << indent << "RooNumConvolution begin cache" << endl ;
if (_init) {
_cloneVar->printCompactTree(os,Form("%s[Var]",indent)) ;
_clonePdf->printCompactTree(os,Form("%s[Pdf]",indent)) ;
_cloneModel->printCompactTree(os,Form("%s[Mod]",indent)) ;
}
os << indent << "RooNumConvolution end cache" << endl ;
}
RooNumConvolution.cxx:100 RooNumConvolution.cxx:101 RooNumConvolution.cxx:102 RooNumConvolution.cxx:103 RooNumConvolution.cxx:104 RooNumConvolution.cxx:105 RooNumConvolution.cxx:106 RooNumConvolution.cxx:107 RooNumConvolution.cxx:108 RooNumConvolution.cxx:109 RooNumConvolution.cxx:110 RooNumConvolution.cxx:111 RooNumConvolution.cxx:112 RooNumConvolution.cxx:113 RooNumConvolution.cxx:114 RooNumConvolution.cxx:115 RooNumConvolution.cxx:116 RooNumConvolution.cxx:117 RooNumConvolution.cxx:118 RooNumConvolution.cxx:119 RooNumConvolution.cxx:120 RooNumConvolution.cxx:121 RooNumConvolution.cxx:122 RooNumConvolution.cxx:123 RooNumConvolution.cxx:124 RooNumConvolution.cxx:125 RooNumConvolution.cxx:126 RooNumConvolution.cxx:127 RooNumConvolution.cxx:128 RooNumConvolution.cxx:129 RooNumConvolution.cxx:130 RooNumConvolution.cxx:131 RooNumConvolution.cxx:132 RooNumConvolution.cxx:133 RooNumConvolution.cxx:134 RooNumConvolution.cxx:135 RooNumConvolution.cxx:136 RooNumConvolution.cxx:137 RooNumConvolution.cxx:138 RooNumConvolution.cxx:139 RooNumConvolution.cxx:140 RooNumConvolution.cxx:141 RooNumConvolution.cxx:142 RooNumConvolution.cxx:143 RooNumConvolution.cxx:144 RooNumConvolution.cxx:145 RooNumConvolution.cxx:146 RooNumConvolution.cxx:147 RooNumConvolution.cxx:148 RooNumConvolution.cxx:149 RooNumConvolution.cxx:150 RooNumConvolution.cxx:151 RooNumConvolution.cxx:152 RooNumConvolution.cxx:153 RooNumConvolution.cxx:154 RooNumConvolution.cxx:155 RooNumConvolution.cxx:156 RooNumConvolution.cxx:157 RooNumConvolution.cxx:158 RooNumConvolution.cxx:159 RooNumConvolution.cxx:160 RooNumConvolution.cxx:161 RooNumConvolution.cxx:162 RooNumConvolution.cxx:163 RooNumConvolution.cxx:164 RooNumConvolution.cxx:165 RooNumConvolution.cxx:166 RooNumConvolution.cxx:167 RooNumConvolution.cxx:168 RooNumConvolution.cxx:169 RooNumConvolution.cxx:170 RooNumConvolution.cxx:171 RooNumConvolution.cxx:172 RooNumConvolution.cxx:173 RooNumConvolution.cxx:174 RooNumConvolution.cxx:175 RooNumConvolution.cxx:176 RooNumConvolution.cxx:177 RooNumConvolution.cxx:178 RooNumConvolution.cxx:179 RooNumConvolution.cxx:180 RooNumConvolution.cxx:181 RooNumConvolution.cxx:182 RooNumConvolution.cxx:183 RooNumConvolution.cxx:184 RooNumConvolution.cxx:185 RooNumConvolution.cxx:186 RooNumConvolution.cxx:187 RooNumConvolution.cxx:188 RooNumConvolution.cxx:189 RooNumConvolution.cxx:190 RooNumConvolution.cxx:191 RooNumConvolution.cxx:192 RooNumConvolution.cxx:193 RooNumConvolution.cxx:194 RooNumConvolution.cxx:195 RooNumConvolution.cxx:196 RooNumConvolution.cxx:197 RooNumConvolution.cxx:198 RooNumConvolution.cxx:199 RooNumConvolution.cxx:200 RooNumConvolution.cxx:201 RooNumConvolution.cxx:202 RooNumConvolution.cxx:203 RooNumConvolution.cxx:204 RooNumConvolution.cxx:205 RooNumConvolution.cxx:206 RooNumConvolution.cxx:207 RooNumConvolution.cxx:208 RooNumConvolution.cxx:209 RooNumConvolution.cxx:210 RooNumConvolution.cxx:211 RooNumConvolution.cxx:212 RooNumConvolution.cxx:213 RooNumConvolution.cxx:214 RooNumConvolution.cxx:215 RooNumConvolution.cxx:216 RooNumConvolution.cxx:217 RooNumConvolution.cxx:218 RooNumConvolution.cxx:219 RooNumConvolution.cxx:220 RooNumConvolution.cxx:221 RooNumConvolution.cxx:222 RooNumConvolution.cxx:223 RooNumConvolution.cxx:224 RooNumConvolution.cxx:225 RooNumConvolution.cxx:226 RooNumConvolution.cxx:227 RooNumConvolution.cxx:228 RooNumConvolution.cxx:229 RooNumConvolution.cxx:230 RooNumConvolution.cxx:231 RooNumConvolution.cxx:232 RooNumConvolution.cxx:233 RooNumConvolution.cxx:234 RooNumConvolution.cxx:235 RooNumConvolution.cxx:236 RooNumConvolution.cxx:237 RooNumConvolution.cxx:238 RooNumConvolution.cxx:239 RooNumConvolution.cxx:240 RooNumConvolution.cxx:241 RooNumConvolution.cxx:242 RooNumConvolution.cxx:243 RooNumConvolution.cxx:244 RooNumConvolution.cxx:245 RooNumConvolution.cxx:246 RooNumConvolution.cxx:247 RooNumConvolution.cxx:248 RooNumConvolution.cxx:249 RooNumConvolution.cxx:250 RooNumConvolution.cxx:251 RooNumConvolution.cxx:252 RooNumConvolution.cxx:253 RooNumConvolution.cxx:254 RooNumConvolution.cxx:255 RooNumConvolution.cxx:256 RooNumConvolution.cxx:257 RooNumConvolution.cxx:258 RooNumConvolution.cxx:259 RooNumConvolution.cxx:260 RooNumConvolution.cxx:261 RooNumConvolution.cxx:262 RooNumConvolution.cxx:263 RooNumConvolution.cxx:264 RooNumConvolution.cxx:265 RooNumConvolution.cxx:266 RooNumConvolution.cxx:267 RooNumConvolution.cxx:268 RooNumConvolution.cxx:269 RooNumConvolution.cxx:270 RooNumConvolution.cxx:271 RooNumConvolution.cxx:272 RooNumConvolution.cxx:273 RooNumConvolution.cxx:274 RooNumConvolution.cxx:275 RooNumConvolution.cxx:276 RooNumConvolution.cxx:277 RooNumConvolution.cxx:278 RooNumConvolution.cxx:279 RooNumConvolution.cxx:280 RooNumConvolution.cxx:281 RooNumConvolution.cxx:282 RooNumConvolution.cxx:283 RooNumConvolution.cxx:284 RooNumConvolution.cxx:285 RooNumConvolution.cxx:286 RooNumConvolution.cxx:287 RooNumConvolution.cxx:288 RooNumConvolution.cxx:289 RooNumConvolution.cxx:290 RooNumConvolution.cxx:291 RooNumConvolution.cxx:292 RooNumConvolution.cxx:293 RooNumConvolution.cxx:294 RooNumConvolution.cxx:295 RooNumConvolution.cxx:296 RooNumConvolution.cxx:297 RooNumConvolution.cxx:298 RooNumConvolution.cxx:299 RooNumConvolution.cxx:300 RooNumConvolution.cxx:301 RooNumConvolution.cxx:302 RooNumConvolution.cxx:303 RooNumConvolution.cxx:304 RooNumConvolution.cxx:305 RooNumConvolution.cxx:306 RooNumConvolution.cxx:307 RooNumConvolution.cxx:308 RooNumConvolution.cxx:309 RooNumConvolution.cxx:310 RooNumConvolution.cxx:311 RooNumConvolution.cxx:312 RooNumConvolution.cxx:313 RooNumConvolution.cxx:314 RooNumConvolution.cxx:315 RooNumConvolution.cxx:316 RooNumConvolution.cxx:317 RooNumConvolution.cxx:318 RooNumConvolution.cxx:319 RooNumConvolution.cxx:320 RooNumConvolution.cxx:321 RooNumConvolution.cxx:322 RooNumConvolution.cxx:323 RooNumConvolution.cxx:324 RooNumConvolution.cxx:325 RooNumConvolution.cxx:326 RooNumConvolution.cxx:327 RooNumConvolution.cxx:328 RooNumConvolution.cxx:329 RooNumConvolution.cxx:330 RooNumConvolution.cxx:331 RooNumConvolution.cxx:332 RooNumConvolution.cxx:333 RooNumConvolution.cxx:334 RooNumConvolution.cxx:335 RooNumConvolution.cxx:336 RooNumConvolution.cxx:337 RooNumConvolution.cxx:338 RooNumConvolution.cxx:339 RooNumConvolution.cxx:340 RooNumConvolution.cxx:341 RooNumConvolution.cxx:342 RooNumConvolution.cxx:343 RooNumConvolution.cxx:344 RooNumConvolution.cxx:345 RooNumConvolution.cxx:346 RooNumConvolution.cxx:347 RooNumConvolution.cxx:348 RooNumConvolution.cxx:349 RooNumConvolution.cxx:350 RooNumConvolution.cxx:351 RooNumConvolution.cxx:352 RooNumConvolution.cxx:353 RooNumConvolution.cxx:354 RooNumConvolution.cxx:355 RooNumConvolution.cxx:356 RooNumConvolution.cxx:357 RooNumConvolution.cxx:358 RooNumConvolution.cxx:359 RooNumConvolution.cxx:360