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
37#include <array>
38
39namespace {
40
41enum RooGaussBasis {
42 noBasis = 0,
43 expBasisMinus = 1,
44 expBasisSum = 2,
45 expBasisPlus = 3,
46 sinBasisMinus = 11,
47 sinBasisSum = 12,
48 sinBasisPlus = 13,
49 cosBasisMinus = 21,
50 cosBasisSum = 22,
51 cosBasisPlus = 23,
52 linBasisPlus = 33,
53 quadBasisPlus = 43,
54 coshBasisMinus = 51,
55 coshBasisSum = 52,
56 coshBasisPlus = 53,
57 sinhBasisMinus = 61,
58 sinhBasisSum = 62,
59 sinhBasisPlus = 63
60};
61
62enum BasisType {
63 none = 0,
64 expBasis = 1,
65 sinBasis = 2,
66 cosBasis = 3,
67 linBasis = 4,
68 quadBasis = 5,
69 coshBasis = 6,
70 sinhBasis = 7
71};
72
73enum BasisSign { Both = 0, Plus = +1, Minus = -1 };
74
75BasisType getBasisType(int basisCode)
76{
77 return static_cast<BasisType>(basisCode == 0 ? 0 : (basisCode / 10) + 1);
78}
79
80} // namespace
81
84
86
87////////////////////////////////////////////////////////////////////////////////
88
89RooGaussModel::RooGaussModel(const char *name, const char *title, RooAbsRealLValue &xIn, RooAbsReal &_mean,
90 RooAbsReal &_sigma)
91 : RooGaussModel{name, title, xIn, _mean, _sigma, RooRealConstant::value(1), RooRealConstant::value(1)}
92{
93}
94
95////////////////////////////////////////////////////////////////////////////////
96
97RooGaussModel::RooGaussModel(const char *name, const char *title, RooAbsRealLValue &xIn, RooAbsReal &_mean,
98 RooAbsReal &_sigma, RooAbsReal &_msSF)
99 : RooGaussModel{name, title, xIn, _mean, _sigma, _msSF, _msSF}
100{
101}
102
103////////////////////////////////////////////////////////////////////////////////
104
105RooGaussModel::RooGaussModel(const char *name, const char *title, RooAbsRealLValue& xIn,
106 RooAbsReal& _mean, RooAbsReal& _sigma,
107 RooAbsReal& _meanSF, RooAbsReal& _sigmaSF) :
108 RooResolutionModel(name,title,xIn),
109 _flatSFInt(false),
110 _asympInt(false),
111 mean("mean","Mean",this,_mean),
112 sigma("sigma","Width",this,_sigma),
113 msf("msf","Mean Scale Factor",this,_meanSF),
114 ssf("ssf","Sigma Scale Factor",this,_sigmaSF)
115{
116}
117
118////////////////////////////////////////////////////////////////////////////////
119
122 _flatSFInt(other._flatSFInt),
123 _asympInt(other._asympInt),
124 mean("mean",this,other.mean),
125 sigma("sigma",this,other.sigma),
126 msf("msf",this,other.msf),
127 ssf("ssf",this,other.ssf)
128{
129}
130
131////////////////////////////////////////////////////////////////////////////////
132
134{
135 std::string str = name;
136
137 // Remove whitespaces from the input string
138 str.erase(remove(str.begin(),str.end(),' '),str.end());
139
140 if (str == "exp(-@0/@1)") return expBasisPlus ;
141 if (str == "exp(@0/@1)") return expBasisMinus ;
142 if (str == "exp(-abs(@0)/@1)") return expBasisSum ;
143 if (str == "exp(-@0/@1)*sin(@0*@2)") return sinBasisPlus ;
144 if (str == "exp(@0/@1)*sin(@0*@2)") return sinBasisMinus ;
145 if (str == "exp(-abs(@0)/@1)*sin(@0*@2)") return sinBasisSum ;
146 if (str == "exp(-@0/@1)*cos(@0*@2)") return cosBasisPlus ;
147 if (str == "exp(@0/@1)*cos(@0*@2)") return cosBasisMinus ;
148 if (str == "exp(-abs(@0)/@1)*cos(@0*@2)") return cosBasisSum ;
149 if (str == "(@0/@1)*exp(-@0/@1)") return linBasisPlus ;
150 if (str == "(@0/@1)*(@0/@1)*exp(-@0/@1)") return quadBasisPlus ;
151 if (str == "exp(-@0/@1)*cosh(@0*@2/2)") return coshBasisPlus;
152 if (str == "exp(@0/@1)*cosh(@0*@2/2)") return coshBasisMinus;
153 if (str == "exp(-abs(@0)/@1)*cosh(@0*@2/2)") return coshBasisSum;
154 if (str == "exp(-@0/@1)*sinh(@0*@2/2)") return sinhBasisPlus;
155 if (str == "exp(@0/@1)*sinh(@0*@2/2)") return sinhBasisMinus;
156 if (str == "exp(-abs(@0)/@1)*sinh(@0*@2/2)") return sinhBasisSum;
157
158 return 0;
159}
160
161////////////////////////////////////////////////////////////////////////////////
162
164{
165 auto arg1 = static_cast<RooAbsReal*>(basis().getParameter(1));
166 auto arg2 = static_cast<RooAbsReal*>(basis().getParameter(2));
167 double param1 = arg1 ? arg1->getVal() : 0.0;
168 double param2 = arg2 ? arg2->getVal() : 0.0;
169 return evaluate(x, mean * msf, sigma * ssf, param1, param2, _basisCode);
170}
171
173{
174 std::span<double> output = ctx.output();
175 std::size_t size = output.size();
176
177 auto xVals = ctx.at(x);
178 auto meanVals = ctx.at(mean);
179 auto meanSfVals = ctx.at(msf);
180 auto sigmaVals = ctx.at(sigma);
181 auto sigmaSfVals = ctx.at(ssf);
182
183 auto param1 = static_cast<RooAbsReal *>(basis().getParameter(1));
184 auto param2 = static_cast<RooAbsReal *>(basis().getParameter(2));
185 const double zeroVal = 0.0;
186 auto param1Vals = param1 ? ctx.at(param1) : std::span<const double>{&zeroVal, 1};
187 auto param2Vals = param2 ? ctx.at(param2) : std::span<const double>{&zeroVal, 1};
188
189 BasisType basisType = getBasisType(_basisCode);
190 double basisSign = _basisCode - 10 * (basisType - 1) - 2;
191
192 // We have an implementation also for CUDA right now only for the most used
193 // basis type, which is expBasis. If the need to support other basis types
194 // arises, they can be implemented following this example. Remember to also
195 // adapt RooGaussModel::canComputeBatchWithCuda().
196 if (basisType == expBasis) {
197 std::array<double, 1> extraArgs{basisSign};
199 {xVals, meanVals, meanSfVals, sigmaVals, sigmaSfVals, param1Vals}, extraArgs);
200 return;
201 }
202
203 // For now, if the arrays don't have the expected input shape, fall back to the scalar mode
204 if (xVals.size() != size || meanVals.size() != 1 || meanSfVals.size() != 1 || sigmaVals.size() != 1 ||
205 sigmaSfVals.size() != 1 || param1Vals.size() != 1 || param2Vals.size() != 1) {
206 return RooAbsPdf::doEval(ctx);
207 }
208
209 for (unsigned int i = 0; i < size; ++i) {
210 output[i] = evaluate(xVals[i], meanVals[0] * meanSfVals[0], sigmaVals[0] * sigmaSfVals[0], param1Vals[0],
211 param2Vals[0], _basisCode);
212 }
213}
214
215double RooGaussModel::evaluate(double x, double mean, double sigma, double param1, double param2, int basisCode)
216{
217 // *** 1st form: Straight Gaussian, used for unconvoluted PDF or expBasis with 0 lifetime ***
218 static double root2(std::sqrt(2.)) ;
219 static double root2pi(std::sqrt(2.*std::atan2(0.,-1.))) ;
220 static double rootpi(std::sqrt(std::atan2(0.,-1.))) ;
221
222 BasisType basisType = getBasisType(basisCode);
223 BasisSign basisSign = (BasisSign)( basisCode - 10*(basisType-1) - 2 ) ;
224
225 double tau = (basisCode!=noBasis) ? param1 : 0.0;
226 if (basisType == coshBasis && basisCode!=noBasis ) {
227 double dGamma = param2;
228 if (dGamma==0) basisType = expBasis;
229 }
230
231 if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
232 double xprime = (x-mean)/sigma ;
233 double result = std::exp(-0.5*xprime*xprime)/(sigma*root2pi) ;
234 if (basisCode!=0 && basisSign==Both) result *= 2 ;
235 return result ;
236 }
237
238 // *** 2nd form: 0, used for sinBasis, linBasis, and quadBasis with tau=0 ***
239 if (tau==0) {
240 return 0. ;
241 }
242
243 // *** 3nd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
244 double omega = (basisType==sinBasis || basisType==cosBasis) ? param2 : 0 ;
245 double dgamma = (basisType==sinhBasis || basisType==coshBasis) ? param2 : 0 ;
246 double _x = omega *tau ;
247 double _y = tau*dgamma/2;
248 double xprime = (x-mean)/tau ;
249 double c = sigma/(root2*tau) ;
250 double u = xprime/(2*c) ;
251
252 if (basisType==expBasis || (basisType==cosBasis && _x==0.)) {
253 double result(0) ;
254 if (basisSign!=Minus) result += evalCerf(0,-u,c).real();
255 if (basisSign!=Plus) result += evalCerf(0, u,c).real();
256 return result ;
257 }
258
259 // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
260 if (basisType==sinBasis) {
261 double result(0) ;
262 if (_x==0.) return result ;
263 if (basisSign!=Minus) result += -evalCerf(-_x,-u,c).imag();
264 if (basisSign!=Plus) result += -evalCerf( _x, u,c).imag();
265 return result ;
266 }
267
268 // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
269 if (basisType==cosBasis) {
270 double result(0) ;
271 if (basisSign!=Minus) result += evalCerf(-_x,-u,c).real();
272 if (basisSign!=Plus) result += evalCerf( _x, u,c).real();
273 return result ;
274 }
275
276 // ***8th form: Convolution with exp(-|t|/tau)*cosh(dgamma*t/2), used for coshBasisSum ***
277 if (basisType==coshBasis || basisType ==sinhBasis) {
278 double result(0);
279 int sgn = ( basisType == coshBasis ? +1 : -1 );
280 if (basisSign!=Minus) result += 0.5*( evalCerf(0,-u,c*(1-_y)).real()+sgn*evalCerf(0,-u,c*(1+_y)).real()) ;
281 if (basisSign!=Plus) result += 0.5*(sgn*evalCerf(0, u,c*(1-_y)).real()+ evalCerf(0, u,c*(1+_y)).real()) ;
282 return result ;
283 }
284
285 // *** 6th form: Convolution with (t/tau)*exp(-t/tau), used for linBasis ***
286 if (basisType==linBasis) {
287 R__ASSERT(basisSign==Plus); // This should only be for positive times
288
289 double f0 = std::exp(-xprime+c*c) * RooMath::erfc(-u+c);
290 double f1 = std::exp(-u*u);
291 return (xprime - 2*c*c)*f0 + (2*c/rootpi)*f1 ;
292 }
293
294 // *** 7th form: Convolution with (t/tau)^2*exp(-t/tau), used for quadBasis ***
295 if (basisType==quadBasis) {
296 R__ASSERT(basisSign==Plus); // This should only be for positive times
297
298 double f0 = std::exp(-xprime+c*c) * RooMath::erfc(-u+c);
299 double f1 = std::exp(-u*u);
300 double x2c2 = xprime - 2*c*c;
301 return ( x2c2*x2c2*f0 + (2*c/rootpi)*x2c2*f1 + 2*c*c*f0 );
302 }
303
304 R__ASSERT(0) ;
305 return 0 ;
306}
307
308////////////////////////////////////////////////////////////////////////////////
309
310Int_t RooGaussModel::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
311{
312 switch(_basisCode) {
313
314 // Analytical integration capability of raw PDF
315 case noBasis:
316
317 // Optionally advertise flat integral over sigma scale factor
318 if (_flatSFInt) {
319 if (matchArgs(allVars,analVars,RooArgSet(convVar(),ssf.arg()))) return 2 ;
320 }
321
322 if (matchArgs(allVars,analVars,convVar())) return 1 ;
323 break ;
324
325 // Analytical integration capability of convoluted PDF
326 case expBasisPlus:
327 case expBasisMinus:
328 case expBasisSum:
329 case sinBasisPlus:
330 case sinBasisMinus:
331 case sinBasisSum:
332 case cosBasisPlus:
333 case cosBasisMinus:
334 case cosBasisSum:
335 case linBasisPlus:
336 case quadBasisPlus:
337 case coshBasisMinus:
338 case coshBasisPlus:
339 case coshBasisSum:
340 case sinhBasisMinus:
341 case sinhBasisPlus:
342 case sinhBasisSum:
343
344 // Optionally advertise flat integral over sigma scale factor
345 if (_flatSFInt) {
346
347 if (matchArgs(allVars,analVars,RooArgSet(convVar(),ssf.arg()))) {
348 return 2 ;
349 }
350 }
351
352 if (matchArgs(allVars,analVars,convVar())) return 1 ;
353 break ;
354 }
355
356 return 0 ;
357}
358
359////////////////////////////////////////////////////////////////////////////////
360
361double RooGaussModel::analyticalIntegral(Int_t code, const char* rangeName) const
362{
363 static const double root2 = std::sqrt(2.) ;
364 //static double rootPiBy2 = std::sqrt(std::atan2(0.0,-1.0)/2.0);
365 static const double rootpi = std::sqrt(std::atan2(0.0,-1.0));
366 double ssfInt(1.0) ;
367
368 // Code must be 1 or 2
369 R__ASSERT(code==1||code==2) ;
370 if (code==2) ssfInt = (ssf.max(rangeName)-ssf.min(rangeName)) ;
371
372 BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
373 BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
374
375 // *** 1st form: Straight Gaussian, used for unconvoluted PDF or expBasis with 0 lifetime ***
376 double tau = (_basisCode!=noBasis)?(static_cast<RooAbsReal*>(basis().getParameter(1)))->getVal():0 ;
377 if (basisType == coshBasis && _basisCode!=noBasis ) {
378 double dGamma = (static_cast<RooAbsReal*>(basis().getParameter(2)))->getVal();
379 if (dGamma==0) basisType = expBasis;
380 }
381 if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
382 double xscale = root2*(sigma*ssf);
383 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 1st form" << std::endl ;
384
385 double xpmin = (x.min(rangeName)-(mean*msf))/xscale ;
386 double xpmax = (x.max(rangeName)-(mean*msf))/xscale ;
387
388 double result ;
389 if (_asympInt) { // modified FMV, 07/24/03
390 result = 1.0 ;
391 } else {
392 result = 0.5*(RooMath::erf(xpmax)-RooMath::erf(xpmin)) ;
393 }
394
395 if (_basisCode!=0 && basisSign==Both) result *= 2 ;
396 //cout << "Integral 1st form " << " result= " << result*ssfInt << std::endl;
397 if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 1 " << std::endl; }
398 return result*ssfInt ;
399 }
400
401
402 double omega = ((basisType==sinBasis)||(basisType==cosBasis)) ? (static_cast<RooAbsReal*>(basis().getParameter(2)))->getVal() : 0 ;
403 double dgamma =((basisType==sinhBasis)||(basisType==coshBasis)) ? (static_cast<RooAbsReal*>(basis().getParameter(2)))->getVal() : 0 ;
404
405 // *** 2nd form: unity, used for sinBasis and linBasis with tau=0 (PDF is zero) ***
406 if (tau==0) {
407 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 2nd form" << std::endl ;
408 return 0. ;
409 }
410
411 // *** 3rd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
412 double c = (sigma*ssf)/(root2*tau) ;
413 double xpmin = (x.min(rangeName)-(mean*msf))/tau ;
414 double xpmax = (x.max(rangeName)-(mean*msf))/tau ;
415 double umin = xpmin/(2*c) ;
416 double umax = xpmax/(2*c) ;
417
418 if (basisType==expBasis || (basisType==cosBasis && omega==0.)) {
419 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 3d form tau=" << tau << std::endl ;
420 double result(0) ;
421 if (basisSign!=Minus) result += evalCerfInt(+1,0,tau,-umin,-umax,c).real();
422 if (basisSign!=Plus) result += evalCerfInt(-1,0,tau, umin, umax,c).real();
423 if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 3 " << std::endl; }
424 return result*ssfInt ;
425 }
426
427 // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
428 double _x = omega * tau ;
429 double _y = tau*dgamma/2;
430
431 if (basisType==sinBasis) {
432 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 4th form omega = " << omega << ", tau = " << tau << std::endl ;
433 double result(0) ;
434 if (_x==0) return result*ssfInt ;
435 if (basisSign!=Minus) result += -1*evalCerfInt(+1,-_x,tau,-umin,-umax,c).imag();
436 if (basisSign!=Plus) result += -1*evalCerfInt(-1, _x,tau, umin, umax,c).imag();
437 if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 4 " << std::endl; }
438 return result*ssfInt ;
439 }
440
441 // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
442 if (basisType==cosBasis) {
443 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 5th form omega = " << omega << ", tau = " << tau << std::endl ;
444 double result(0) ;
445 if (basisSign!=Minus) result += evalCerfInt(+1,-_x,tau,-umin,-umax,c).real();
446 if (basisSign!=Plus) result += evalCerfInt(-1, _x,tau, umin, umax,c).real();
447 if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 5 " << std::endl; }
448 return result*ssfInt ;
449 }
450
451 // *** 8th form: Convolution with exp(-|t|/tau)*cosh(dgamma*t/2), used for coshBasis ***
452 // *** 9th form: Convolution with exp(-|t|/tau)*sinh(dgamma*t/2), used for sinhBasis ***
453 if (basisType==coshBasis || basisType == sinhBasis) {
454 if (verboseEval()>0) {std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 8th form tau=" << tau << std::endl ; }
455 double result(0) ;
456 int sgn = ( basisType == coshBasis ? +1 : -1 );
457 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());
458 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());
459 if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 6 " << std::endl; }
460 return result*ssfInt ;
461 }
462
463 // *** 6th form: Convolution with (t/tau)*exp(-t/tau), used for linBasis ***
464 if (basisType==linBasis) {
465 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 6th form tau=" << tau << std::endl ;
466
467 double f0 = RooMath::erf(-umax) - RooMath::erf(-umin);
468 double f1 = std::exp(-umax*umax) - std::exp(-umin*umin);
469
470 double tmp1 = std::exp(-xpmax)*RooMath::erfc(-umax + c);
471 double tmp2 = std::exp(-xpmin)*RooMath::erfc(-umin + c);
472
473 double f2 = tmp1 - tmp2;
474 double f3 = xpmax*tmp1 - xpmin*tmp2;
475
476 double expc2 = std::exp(c*c);
477
478 return -tau*( f0 +
479 (2*c/rootpi)*f1 +
480 (1 - 2*c*c)*expc2*f2 +
481 expc2*f3
482 )*ssfInt;
483 }
484
485 // *** 7th form: Convolution with (t/tau)*(t/tau)*exp(-t/tau), used for quadBasis ***
486 if (basisType==quadBasis) {
487 if (verboseEval()>0) std::cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 7th form tau=" << tau << std::endl ;
488
489 double f0 = RooMath::erf(-umax) - RooMath::erf(-umin);
490
491 double tmpA1 = std::exp(-umax*umax);
492 double tmpA2 = std::exp(-umin*umin);
493
494 double f1 = tmpA1 - tmpA2;
495 double f2 = umax*tmpA1 - umin*tmpA2;
496
497 double tmpB1 = std::exp(-xpmax)*RooMath::erfc(-umax + c);
498 double tmpB2 = std::exp(-xpmin)*RooMath::erfc(-umin + c);
499
500 double f3 = tmpB1 - tmpB2;
501 double f4 = xpmax*tmpB1 - xpmin*tmpB2;
502 double f5 = xpmax*xpmax*tmpB1 - xpmin*xpmin*tmpB2;
503
504 double expc2 = std::exp(c*c);
505
506 return -tau*( 2*f0 +
507 (4*c/rootpi)*((1-c*c)*f1 + c*f2) +
508 (2*c*c*(2*c*c-1) + 2)*expc2*f3 - (4*c*c-2)*expc2*f4 + expc2*f5
509 )*ssfInt;
510 }
511 R__ASSERT(0) ;
512 return 0 ;
513}
514
515
516////////////////////////////////////////////////////////////////////////////////
517
518std::complex<double> RooGaussModel::evalCerfInt(double sign, double _x, double tau, double umin, double umax, double c) const
519{
520 std::complex<double> diff(2., 0.);
521 if (!_asympInt) {
522 diff = evalCerf(_x,umin,c);
523 diff -= evalCerf(_x,umax,c);
524 diff += RooMath::erf(umin) - RooMath::erf(umax);
525 diff *= sign;
526 }
527 diff *= std::complex<double>(1., _x);
528 diff *= tau / (1.+_x*_x);
529 return diff;
530}
531
532////////////////////////////////////////////////////////////////////////////////
533
534Int_t RooGaussModel::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
535{
536 return matchArgs(directVars,generateVars,x) ? 1 : 0;
537}
538
539////////////////////////////////////////////////////////////////////////////////
540
542{
543 R__ASSERT(code==1) ;
544 double xmin = x.min();
545 double xmax = x.max();
546 TRandom *generator = RooRandom::randomGenerator();
547 while(true) {
548 double xgen = generator->Gaus(mean*msf,sigma*ssf);
549 if (xgen<xmax && xgen>xmin) {
550 x = xgen ;
551 return ;
552 }
553 }
554}
555
557{
558 return getBasisType(_basisCode) == expBasis;
559}
#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:382
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
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
bool matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
virtual void doEval(RooFit::EvalContext &) const
Base function for computing multiple values of a RooAbsReal.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
std::span< const double > at(RooAbsArg const *arg, RooAbsArg const *caller=nullptr)
std::span< double > output()
RooBatchCompute::Config config(RooAbsArg const *arg) const
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
void doEval(RooFit::EvalContext &) const override
Base function for computing multiple values of a RooAbsReal.
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...
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:48
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.
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
void compute(Config cfg, Computer comp, std::span< double > output, VarSpan vars, ArgSpan 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()