Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooGaussModel.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$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/** \class RooGaussModel
18 \ingroup Roofit
19
20Class RooGaussModel implements a RooResolutionModel that models a Gaussian
21distribution. Object of class RooGaussModel can be used
22for analytical convolutions with classes inheriting from RooAbsAnaConvPdf
23**/
24
25#include "TMath.h"
26#include "Riostream.h"
27#include "RooGaussModel.h"
28#include "RooMath.h"
29#include "RooRealConstant.h"
30#include "RooRandom.h"
31#include "RooBatchCompute.h"
32
33#include "TError.h"
34
36
37namespace {
38
39enum RooGaussBasis {
40 noBasis = 0,
41 expBasisMinus = 1,
42 expBasisSum = 2,
43 expBasisPlus = 3,
44 sinBasisMinus = 11,
45 sinBasisSum = 12,
46 sinBasisPlus = 13,
47 cosBasisMinus = 21,
48 cosBasisSum = 22,
49 cosBasisPlus = 23,
50 linBasisPlus = 33,
51 quadBasisPlus = 43,
52 coshBasisMinus = 51,
53 coshBasisSum = 52,
54 coshBasisPlus = 53,
55 sinhBasisMinus = 61,
56 sinhBasisSum = 62,
57 sinhBasisPlus = 63
58};
59
60enum BasisType {
61 none = 0,
62 expBasis = 1,
63 sinBasis = 2,
64 cosBasis = 3,
65 linBasis = 4,
66 quadBasis = 5,
67 coshBasis = 6,
68 sinhBasis = 7
69};
70
71enum BasisSign { Both = 0, Plus = +1, Minus = -1 };
72
73BasisType getBasisType(int basisCode)
74{
75 return static_cast<BasisType>(basisCode == 0 ? 0 : (basisCode / 10) + 1);
76}
77
78} // namespace
79
82
84
85////////////////////////////////////////////////////////////////////////////////
86
87RooGaussModel::RooGaussModel(const char *name, const char *title, RooAbsRealLValue &xIn, RooAbsReal &_mean,
88 RooAbsReal &_sigma)
89 : RooGaussModel{name, title, xIn, _mean, _sigma, RooRealConstant::value(1), RooRealConstant::value(1)}
90{
91}
92
93////////////////////////////////////////////////////////////////////////////////
94
95RooGaussModel::RooGaussModel(const char *name, const char *title, RooAbsRealLValue &xIn, RooAbsReal &_mean,
96 RooAbsReal &_sigma, RooAbsReal &_msSF)
97 : RooGaussModel{name, title, xIn, _mean, _sigma, _msSF, _msSF}
98{
99}
100
101////////////////////////////////////////////////////////////////////////////////
102
103RooGaussModel::RooGaussModel(const char *name, const char *title, RooAbsRealLValue& xIn,
104 RooAbsReal& _mean, RooAbsReal& _sigma,
105 RooAbsReal& _meanSF, RooAbsReal& _sigmaSF) :
106 RooResolutionModel(name,title,xIn),
107 _flatSFInt(false),
108 _asympInt(false),
109 mean("mean","Mean",this,_mean),
110 sigma("sigma","Width",this,_sigma),
111 msf("msf","Mean Scale Factor",this,_meanSF),
112 ssf("ssf","Sigma Scale Factor",this,_sigmaSF)
113{
114}
115
116////////////////////////////////////////////////////////////////////////////////
117
120 _flatSFInt(other._flatSFInt),
121 _asympInt(other._asympInt),
122 mean("mean",this,other.mean),
123 sigma("sigma",this,other.sigma),
124 msf("msf",this,other.msf),
125 ssf("ssf",this,other.ssf)
126{
127}
128
129////////////////////////////////////////////////////////////////////////////////
130
132{
133 std::string str = name;
134
135 // Remove whitespaces from the input string
136 str.erase(remove(str.begin(),str.end(),' '),str.end());
137
138 if (str == "exp(-@0/@1)") return expBasisPlus ;
139 if (str == "exp(@0/@1)") return expBasisMinus ;
140 if (str == "exp(-abs(@0)/@1)") return expBasisSum ;
141 if (str == "exp(-@0/@1)*sin(@0*@2)") return sinBasisPlus ;
142 if (str == "exp(@0/@1)*sin(@0*@2)") return sinBasisMinus ;
143 if (str == "exp(-abs(@0)/@1)*sin(@0*@2)") return sinBasisSum ;
144 if (str == "exp(-@0/@1)*cos(@0*@2)") return cosBasisPlus ;
145 if (str == "exp(@0/@1)*cos(@0*@2)") return cosBasisMinus ;
146 if (str == "exp(-abs(@0)/@1)*cos(@0*@2)") return cosBasisSum ;
147 if (str == "(@0/@1)*exp(-@0/@1)") return linBasisPlus ;
148 if (str == "(@0/@1)*(@0/@1)*exp(-@0/@1)") return quadBasisPlus ;
149 if (str == "exp(-@0/@1)*cosh(@0*@2/2)") return coshBasisPlus;
150 if (str == "exp(@0/@1)*cosh(@0*@2/2)") return coshBasisMinus;
151 if (str == "exp(-abs(@0)/@1)*cosh(@0*@2/2)") return coshBasisSum;
152 if (str == "exp(-@0/@1)*sinh(@0*@2/2)") return sinhBasisPlus;
153 if (str == "exp(@0/@1)*sinh(@0*@2/2)") return sinhBasisMinus;
154 if (str == "exp(-abs(@0)/@1)*sinh(@0*@2/2)") return sinhBasisSum;
155
156 return 0;
157}
158
159////////////////////////////////////////////////////////////////////////////////
160
162{
163 auto arg1 = static_cast<RooAbsReal*>(basis().getParameter(1));
164 auto arg2 = static_cast<RooAbsReal*>(basis().getParameter(2));
165 double param1 = arg1 ? arg1->getVal() : 0.0;
166 double param2 = arg2 ? arg2->getVal() : 0.0;
167 return evaluate(x, mean * msf, sigma * ssf, param1, param2, _basisCode);
168}
169
171 RooFit::Detail::DataMap const &dataMap) const
172{
173 auto xVals = dataMap.at(x);
174 auto meanVals = dataMap.at(mean);
175 auto meanSfVals = dataMap.at(msf);
176 auto sigmaVals = dataMap.at(sigma);
177 auto sigmaSfVals = dataMap.at(ssf);
178
179 auto param1 = static_cast<RooAbsReal *>(basis().getParameter(1));
180 auto param2 = static_cast<RooAbsReal *>(basis().getParameter(2));
181 const double zeroVal = 0.0;
182 auto param1Vals = param1 ? dataMap.at(param1) : std::span<const double>{&zeroVal, 1};
183 auto param2Vals = param2 ? dataMap.at(param2) : std::span<const double>{&zeroVal, 1};
184
185 BasisType basisType = getBasisType(_basisCode);
186 double basisSign = _basisCode - 10 * (basisType - 1) - 2;
187
188 // We have an implementation also for CUDA right now only for the most used
189 // basis type, which is expBasis. If the need to support other basis types
190 // arises, they can be implemented following this example. Remember to also
191 // adapt RooGaussModel::canComputeBatchWithCuda().
192 if (basisType == expBasis) {
193 RooBatchCompute::ArgVector extraArgs{basisSign};
195 {xVals, meanVals, meanSfVals, sigmaVals, sigmaSfVals, param1Vals}, extraArgs);
196 return;
197 }
198
199 // For now, if the arrays don't have the expected input shape, fall back to the scalar mode
200 if (xVals.size() != size || meanVals.size() != 1 || meanSfVals.size() != 1 || sigmaVals.size() != 1 ||
201 sigmaSfVals.size() != 1 || param1Vals.size() != 1 || param2Vals.size() != 1) {
202 return RooAbsPdf::computeBatch(output, size, dataMap);
203 }
204
205 for (unsigned int i = 0; i < size; ++i) {
206 output[i] = evaluate(xVals[i], meanVals[0] * meanSfVals[0], sigmaVals[0] * sigmaSfVals[0], param1Vals[0],
207 param2Vals[0], _basisCode);
208 }
209}
210
211double RooGaussModel::evaluate(double x, double mean, double sigma, double param1, double param2, int basisCode)
212{
213 // *** 1st form: Straight Gaussian, used for unconvoluted PDF or expBasis with 0 lifetime ***
214 static double root2(std::sqrt(2.)) ;
215 static double root2pi(std::sqrt(2.*std::atan2(0.,-1.))) ;
216 static double rootpi(std::sqrt(std::atan2(0.,-1.))) ;
217
218 BasisType basisType = getBasisType(basisCode);
219 BasisSign basisSign = (BasisSign)( basisCode - 10*(basisType-1) - 2 ) ;
220
221 double tau = (basisCode!=noBasis) ? param1 : 0.0;
222 if (basisType == coshBasis && basisCode!=noBasis ) {
223 double dGamma = param2;
224 if (dGamma==0) basisType = expBasis;
225 }
226
227 if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
228 double xprime = (x-mean)/sigma ;
229 double result = std::exp(-0.5*xprime*xprime)/(sigma*root2pi) ;
230 if (basisCode!=0 && basisSign==Both) result *= 2 ;
231 return result ;
232 }
233
234 // *** 2nd form: 0, used for sinBasis, linBasis, and quadBasis with tau=0 ***
235 if (tau==0) {
236 return 0. ;
237 }
238
239 // *** 3nd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
240 double omega = (basisType==sinBasis || basisType==cosBasis) ? param2 : 0 ;
241 double dgamma = (basisType==sinhBasis || basisType==coshBasis) ? param2 : 0 ;
242 double _x = omega *tau ;
243 double _y = tau*dgamma/2;
244 double xprime = (x-mean)/tau ;
245 double c = sigma/(root2*tau) ;
246 double u = xprime/(2*c) ;
247
248 if (basisType==expBasis || (basisType==cosBasis && _x==0.)) {
249 double result(0) ;
250 if (basisSign!=Minus) result += evalCerf(0,-u,c).real();
251 if (basisSign!=Plus) result += evalCerf(0, u,c).real();
252 return result ;
253 }
254
255 // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
256 if (basisType==sinBasis) {
257 double result(0) ;
258 if (_x==0.) return result ;
259 if (basisSign!=Minus) result += -evalCerf(-_x,-u,c).imag();
260 if (basisSign!=Plus) result += -evalCerf( _x, u,c).imag();
261 return result ;
262 }
263
264 // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
265 if (basisType==cosBasis) {
266 double result(0) ;
267 if (basisSign!=Minus) result += evalCerf(-_x,-u,c).real();
268 if (basisSign!=Plus) result += evalCerf( _x, u,c).real();
269 return result ;
270 }
271
272 // ***8th form: Convolution with exp(-|t|/tau)*cosh(dgamma*t/2), used for coshBasisSum ***
273 if (basisType==coshBasis || basisType ==sinhBasis) {
274 double result(0);
275 int sgn = ( basisType == coshBasis ? +1 : -1 );
276 if (basisSign!=Minus) result += 0.5*( evalCerf(0,-u,c*(1-_y)).real()+sgn*evalCerf(0,-u,c*(1+_y)).real()) ;
277 if (basisSign!=Plus) result += 0.5*(sgn*evalCerf(0, u,c*(1-_y)).real()+ evalCerf(0, u,c*(1+_y)).real()) ;
278 return result ;
279 }
280
281 // *** 6th form: Convolution with (t/tau)*exp(-t/tau), used for linBasis ***
282 if (basisType==linBasis) {
283 R__ASSERT(basisSign==Plus); // This should only be for positive times
284
285 double f0 = std::exp(-xprime+c*c) * RooMath::erfc(-u+c);
286 double f1 = std::exp(-u*u);
287 return (xprime - 2*c*c)*f0 + (2*c/rootpi)*f1 ;
288 }
289
290 // *** 7th form: Convolution with (t/tau)^2*exp(-t/tau), used for quadBasis ***
291 if (basisType==quadBasis) {
292 R__ASSERT(basisSign==Plus); // This should only be for positive times
293
294 double f0 = std::exp(-xprime+c*c) * RooMath::erfc(-u+c);
295 double f1 = std::exp(-u*u);
296 double x2c2 = xprime - 2*c*c;
297 return ( x2c2*x2c2*f0 + (2*c/rootpi)*x2c2*f1 + 2*c*c*f0 );
298 }
299
300 R__ASSERT(0) ;
301 return 0 ;
302}
303
304////////////////////////////////////////////////////////////////////////////////
305
306Int_t RooGaussModel::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
307{
308 switch(_basisCode) {
309
310 // Analytical integration capability of raw PDF
311 case noBasis:
312
313 // Optionally advertise flat integral over sigma scale factor
314 if (_flatSFInt) {
315 if (matchArgs(allVars,analVars,RooArgSet(convVar(),ssf.arg()))) return 2 ;
316 }
317
318 if (matchArgs(allVars,analVars,convVar())) return 1 ;
319 break ;
320
321 // Analytical integration capability of convoluted PDF
322 case expBasisPlus:
323 case expBasisMinus:
324 case expBasisSum:
325 case sinBasisPlus:
326 case sinBasisMinus:
327 case sinBasisSum:
328 case cosBasisPlus:
329 case cosBasisMinus:
330 case cosBasisSum:
331 case linBasisPlus:
332 case quadBasisPlus:
333 case coshBasisMinus:
334 case coshBasisPlus:
335 case coshBasisSum:
336 case sinhBasisMinus:
337 case sinhBasisPlus:
338 case sinhBasisSum:
339
340 // Optionally advertise flat integral over sigma scale factor
341 if (_flatSFInt) {
342
343 if (matchArgs(allVars,analVars,RooArgSet(convVar(),ssf.arg()))) {
344 return 2 ;
345 }
346 }
347
348 if (matchArgs(allVars,analVars,convVar())) return 1 ;
349 break ;
350 }
351
352 return 0 ;
353}
354
355////////////////////////////////////////////////////////////////////////////////
356
357double RooGaussModel::analyticalIntegral(Int_t code, const char* rangeName) const
358{
359 static const double root2 = std::sqrt(2.) ;
360 //static double rootPiBy2 = std::sqrt(std::atan2(0.0,-1.0)/2.0);
361 static const double rootpi = std::sqrt(std::atan2(0.0,-1.0));
362 double ssfInt(1.0) ;
363
364 // Code must be 1 or 2
365 R__ASSERT(code==1||code==2) ;
366 if (code==2) ssfInt = (ssf.max(rangeName)-ssf.min(rangeName)) ;
367
368 BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
369 BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
370
371 // *** 1st form: Straight Gaussian, used for unconvoluted PDF or expBasis with 0 lifetime ***
372 double tau = (_basisCode!=noBasis)?(static_cast<RooAbsReal*>(basis().getParameter(1)))->getVal():0 ;
373 if (basisType == coshBasis && _basisCode!=noBasis ) {
374 double dGamma = (static_cast<RooAbsReal*>(basis().getParameter(2)))->getVal();
375 if (dGamma==0) basisType = expBasis;
376 }
377 if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
378 double xscale = root2*(sigma*ssf);
379 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 1st form" << std::endl ;
380
381 double xpmin = (x.min(rangeName)-(mean*msf))/xscale ;
382 double xpmax = (x.max(rangeName)-(mean*msf))/xscale ;
383
384 double result ;
385 if (_asympInt) { // modified FMV, 07/24/03
386 result = 1.0 ;
387 } else {
388 result = 0.5*(RooMath::erf(xpmax)-RooMath::erf(xpmin)) ;
389 }
390
391 if (_basisCode!=0 && basisSign==Both) result *= 2 ;
392 //cout << "Integral 1st form " << " result= " << result*ssfInt << std::endl;
393 if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 1 " << std::endl; }
394 return result*ssfInt ;
395 }
396
397
398 double omega = ((basisType==sinBasis)||(basisType==cosBasis)) ? (static_cast<RooAbsReal*>(basis().getParameter(2)))->getVal() : 0 ;
399 double dgamma =((basisType==sinhBasis)||(basisType==coshBasis)) ? (static_cast<RooAbsReal*>(basis().getParameter(2)))->getVal() : 0 ;
400
401 // *** 2nd form: unity, used for sinBasis and linBasis with tau=0 (PDF is zero) ***
402 if (tau==0) {
403 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 2nd form" << std::endl ;
404 return 0. ;
405 }
406
407 // *** 3rd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
408 double c = (sigma*ssf)/(root2*tau) ;
409 double xpmin = (x.min(rangeName)-(mean*msf))/tau ;
410 double xpmax = (x.max(rangeName)-(mean*msf))/tau ;
411 double umin = xpmin/(2*c) ;
412 double umax = xpmax/(2*c) ;
413
414 if (basisType==expBasis || (basisType==cosBasis && omega==0.)) {
415 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 3d form tau=" << tau << std::endl ;
416 double result(0) ;
417 if (basisSign!=Minus) result += evalCerfInt(+1,0,tau,-umin,-umax,c).real();
418 if (basisSign!=Plus) result += evalCerfInt(-1,0,tau, umin, umax,c).real();
419 if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 3 " << std::endl; }
420 return result*ssfInt ;
421 }
422
423 // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
424 double _x = omega * tau ;
425 double _y = tau*dgamma/2;
426
427 if (basisType==sinBasis) {
428 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 4th form omega = " << omega << ", tau = " << tau << std::endl ;
429 double result(0) ;
430 if (_x==0) return result*ssfInt ;
431 if (basisSign!=Minus) result += -1*evalCerfInt(+1,-_x,tau,-umin,-umax,c).imag();
432 if (basisSign!=Plus) result += -1*evalCerfInt(-1, _x,tau, umin, umax,c).imag();
433 if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 4 " << std::endl; }
434 return result*ssfInt ;
435 }
436
437 // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
438 if (basisType==cosBasis) {
439 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 5th form omega = " << omega << ", tau = " << tau << std::endl ;
440 double result(0) ;
441 if (basisSign!=Minus) result += evalCerfInt(+1,-_x,tau,-umin,-umax,c).real();
442 if (basisSign!=Plus) result += evalCerfInt(-1, _x,tau, umin, umax,c).real();
443 if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 5 " << std::endl; }
444 return result*ssfInt ;
445 }
446
447 // *** 8th form: Convolution with exp(-|t|/tau)*cosh(dgamma*t/2), used for coshBasis ***
448 // *** 9th form: Convolution with exp(-|t|/tau)*sinh(dgamma*t/2), used for sinhBasis ***
449 if (basisType==coshBasis || basisType == sinhBasis) {
450 if (verboseEval()>0) {std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 8th form tau=" << tau << std::endl ; }
451 double result(0) ;
452 int sgn = ( basisType == coshBasis ? +1 : -1 );
453 if (basisSign!=Minus) result += 0.5*( evalCerfInt(+1,0,tau/(1-_y),-umin,-umax,c*(1-_y)).real()+ sgn*evalCerfInt(+1,0,tau/(1+_y),-umin,-umax,c*(1+_y)).real());
454 if (basisSign!=Plus) result += 0.5*(sgn*evalCerfInt(-1,0,tau/(1-_y), umin, umax,c*(1-_y)).real()+ evalCerfInt(-1,0,tau/(1+_y), umin, umax,c*(1+_y)).real());
455 if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 6 " << std::endl; }
456 return result*ssfInt ;
457 }
458
459 // *** 6th form: Convolution with (t/tau)*exp(-t/tau), used for linBasis ***
460 if (basisType==linBasis) {
461 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 6th form tau=" << tau << std::endl ;
462
463 double f0 = RooMath::erf(-umax) - RooMath::erf(-umin);
464 double f1 = std::exp(-umax*umax) - std::exp(-umin*umin);
465
466 double tmp1 = std::exp(-xpmax)*RooMath::erfc(-umax + c);
467 double tmp2 = std::exp(-xpmin)*RooMath::erfc(-umin + c);
468
469 double f2 = tmp1 - tmp2;
470 double f3 = xpmax*tmp1 - xpmin*tmp2;
471
472 double expc2 = std::exp(c*c);
473
474 return -tau*( f0 +
475 (2*c/rootpi)*f1 +
476 (1 - 2*c*c)*expc2*f2 +
477 expc2*f3
478 )*ssfInt;
479 }
480
481 // *** 7th form: Convolution with (t/tau)*(t/tau)*exp(-t/tau), used for quadBasis ***
482 if (basisType==quadBasis) {
483 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 7th form tau=" << tau << std::endl ;
484
485 double f0 = RooMath::erf(-umax) - RooMath::erf(-umin);
486
487 double tmpA1 = std::exp(-umax*umax);
488 double tmpA2 = std::exp(-umin*umin);
489
490 double f1 = tmpA1 - tmpA2;
491 double f2 = umax*tmpA1 - umin*tmpA2;
492
493 double tmpB1 = std::exp(-xpmax)*RooMath::erfc(-umax + c);
494 double tmpB2 = std::exp(-xpmin)*RooMath::erfc(-umin + c);
495
496 double f3 = tmpB1 - tmpB2;
497 double f4 = xpmax*tmpB1 - xpmin*tmpB2;
498 double f5 = xpmax*xpmax*tmpB1 - xpmin*xpmin*tmpB2;
499
500 double expc2 = std::exp(c*c);
501
502 return -tau*( 2*f0 +
503 (4*c/rootpi)*((1-c*c)*f1 + c*f2) +
504 (2*c*c*(2*c*c-1) + 2)*expc2*f3 - (4*c*c-2)*expc2*f4 + expc2*f5
505 )*ssfInt;
506 }
507 R__ASSERT(0) ;
508 return 0 ;
509}
510
511
512////////////////////////////////////////////////////////////////////////////////
513
514std::complex<double> RooGaussModel::evalCerfInt(double sign, double _x, double tau, double umin, double umax, double c) const
515{
516 std::complex<double> diff(2., 0.);
517 if (!_asympInt) {
518 diff = evalCerf(_x,umin,c);
519 diff -= evalCerf(_x,umax,c);
520 diff += RooMath::erf(umin) - RooMath::erf(umax);
521 diff *= sign;
522 }
523 diff *= std::complex<double>(1., _x);
524 diff *= tau / (1.+_x*_x);
525 return diff;
526}
527
528////////////////////////////////////////////////////////////////////////////////
529
530Int_t RooGaussModel::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
531{
532 return matchArgs(directVars,generateVars,x) ? 1 : 0;
533}
534
535////////////////////////////////////////////////////////////////////////////////
536
538{
539 R__ASSERT(code==1) ;
540 double xmin = x.min();
541 double xmax = x.max();
542 TRandom *generator = RooRandom::randomGenerator();
543 while(true) {
544 double xgen = generator->Gaus(mean*msf,sigma*ssf);
545 if (xgen<xmax && xgen>xmin) {
546 x = xgen ;
547 return ;
548 }
549 }
550}
551
553{
554 return getBasisType(_basisCode) == expBasis;
555}
#define c(i)
Definition RSha256.hxx:101
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define cxcoutE(a)
#define ClassImp(name)
Definition Rtypes.h:377
#define R__ASSERT(e)
Definition TError.h:118
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
static int verboseEval()
Return global level of verbosity for p.d.f. evaluations.
Abstract base class for objects that represent a real value that may appear on the left hand side of ...
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
virtual void computeBatch(double *output, size_t size, RooFit::Detail::DataMap const &) const
Base function for computing multiple values of a RooAbsReal.
bool matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
RooBatchCompute::Config config(RooAbsArg const *arg) const
Definition DataMap.cxx:40
std::span< const double > at(RooAbsArg const *arg, RooAbsArg const *caller=nullptr)
Definition DataMap.cxx:22
RooAbsArg * getParameter(const char *name) const
Return pointer to parameter with given name.
Class RooGaussModel implements a RooResolutionModel that models a Gaussian distribution.
RooRealProxy sigma
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
std::complex< double > evalCerfInt(double sign, double wt, double tau, double umin, double umax, double c) const
RooRealProxy msf
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
void generateEvent(Int_t code) override
Interface for generation of an event using the algorithm corresponding to the specified code.
RooGaussModel()=default
bool canComputeBatchWithCuda() const override
Int_t basisCode(const char *name) const override
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const override
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
void computeBatch(double *output, size_t size, RooFit::Detail::DataMap const &) const override
Base function for computing multiple values of a RooAbsReal.
RooRealProxy mean
double analyticalIntegral(Int_t code, const char *rangeName) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
RooRealProxy ssf
static std::complex< double > erfc(const std::complex< double > z)
complex erfc function
Definition RooMath.cxx:40
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition RooMath.cxx:59
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition RooRandom.cxx:51
Provides static functions to create and keep track of RooRealVar constants.
RooResolutionModel is the base class for PDFs that represent a resolution model that can be convolute...
Int_t _basisCode
Identifier code for selected basis function.
RooAbsRealLValue & convVar() const
Return the convolution variable of the resolution model.
const RooFormulaVar & basis() const
RooTemplateProxy< RooAbsRealLValue > x
Dependent/convolution variable.
double max(const char *rname=nullptr) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
const T & arg() const
Return reference to object held in proxy.
double min(const char *rname=nullptr) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
This is the base class for the ROOT Random number generators.
Definition TRandom.h:27
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition TRandom.cxx:275
const Double_t sigma
Double_t x[n]
Definition legend1.C:17
TF1 * f1
Definition legend1.C:11
std::vector< double > ArgVector
void compute(Config cfg, Computer comp, RestrictArr output, size_t size, const VarVector &vars, ArgVector &extraArgs)
__roohost__ __roodevice__ STD::complex< double > evalCerfApprox(double _x, double u, double c)
use the approximation: erf(z) = exp(-z*z)/(STD::sqrt(pi)*z) to explicitly cancel the divergent exp(y*...
__roohost__ __roodevice__ STD::complex< double > evalCerf(double swt, double u, double c)
Bool_t IsNaN(Double_t x)
Definition TMath.h:892
static void output()