73enum BasisSign { Both = 0, Plus = +1, Minus = -1 };
75BasisType getBasisType(
int basisCode)
77 return static_cast<BasisType
>(basisCode == 0 ? 0 : (basisCode / 10) + 1);
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)
122 _flatSFInt(other._flatSFInt),
123 _asympInt(other._asympInt),
124 mean(
"mean",this,other.mean),
126 msf(
"msf",this,other.msf),
127 ssf(
"ssf",this,other.ssf)
135 std::string str =
name;
138 str.erase(remove(str.begin(),str.end(),
' '),str.end());
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;
167 double param1 = arg1 ? arg1->
getVal() : 0.0;
168 double param2 = arg2 ? arg2->getVal() : 0.0;
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);
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};
189 BasisType basisType = getBasisType(
_basisCode);
190 double basisSign =
_basisCode - 10 * (basisType - 1) - 2;
196 if (basisType == expBasis) {
197 std::array<double, 1> extraArgs{basisSign};
199 {xVals, meanVals, meanSfVals, sigmaVals, sigmaSfVals, param1Vals}, extraArgs);
204 if (xVals.size() !=
size || meanVals.size() != 1 || meanSfVals.size() != 1 || sigmaVals.size() != 1 ||
205 sigmaSfVals.size() != 1 || param1Vals.size() != 1 || param2Vals.size() != 1) {
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],
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.))) ;
222 BasisType basisType = getBasisType(
basisCode);
223 BasisSign basisSign = (BasisSign)(
basisCode - 10*(basisType-1) - 2 ) ;
225 double tau = (
basisCode!=noBasis) ? param1 : 0.0;
226 if (basisType == coshBasis &&
basisCode!=noBasis ) {
227 double dGamma = param2;
228 if (dGamma==0) basisType = expBasis;
231 if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
233 double result = std::exp(-0.5*xprime*xprime)/(
sigma*root2pi) ;
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) ;
252 if (basisType==expBasis || (basisType==cosBasis && _x==0.)) {
254 if (basisSign!=Minus)
result += evalCerf(0,-u,
c).real();
255 if (basisSign!=Plus)
result += evalCerf(0, u,
c).real();
260 if (basisType==sinBasis) {
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();
269 if (basisType==cosBasis) {
271 if (basisSign!=Minus)
result += evalCerf(-_x,-u,
c).real();
272 if (basisSign!=Plus)
result += evalCerf( _x, u,
c).real();
277 if (basisType==coshBasis || basisType ==sinhBasis) {
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()) ;
286 if (basisType==linBasis) {
290 double f1 = std::exp(-u*u);
291 return (xprime - 2*
c*
c)*f0 + (2*
c/rootpi)*
f1 ;
295 if (basisType==quadBasis) {
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 );
363 static const double root2 = std::sqrt(2.) ;
365 static const double rootpi = std::sqrt(std::atan2(0.0,-1.0));
370 if (code==2) ssfInt = (
ssf.max(rangeName)-
ssf.min(rangeName)) ;
373 BasisSign basisSign = (BasisSign)(
_basisCode - 10*(basisType-1) - 2 ) ;
377 if (basisType == coshBasis &&
_basisCode!=noBasis ) {
379 if (dGamma==0) basisType = expBasis;
381 if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
383 if (
verboseEval()>0) std::cout <<
"RooGaussModel::analyticalIntegral(" <<
GetName() <<
") 1st form" << std::endl ;
385 double xpmin = (
x.min(rangeName)-(
mean*
msf))/xscale ;
386 double xpmax = (
x.max(rangeName)-(
mean*
msf))/xscale ;
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 ;
407 if (
verboseEval()>0) std::cout <<
"RooGaussModel::analyticalIntegral(" <<
GetName() <<
") 2nd form" << std::endl ;
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) ;
418 if (basisType==expBasis || (basisType==cosBasis && omega==0.)) {
419 if (
verboseEval()>0) std::cout <<
"RooGaussModel::analyticalIntegral(" <<
GetName() <<
") 3d form tau=" << tau << std::endl ;
428 double _x = omega * tau ;
429 double _y = tau*dgamma/2;
431 if (basisType==sinBasis) {
432 if (
verboseEval()>0) std::cout <<
"RooGaussModel::analyticalIntegral(" <<
GetName() <<
") 4th form omega = " << omega <<
", tau = " << tau << std::endl ;
434 if (_x==0)
return result*ssfInt ;
442 if (basisType==cosBasis) {
443 if (
verboseEval()>0) std::cout <<
"RooGaussModel::analyticalIntegral(" <<
GetName() <<
") 5th form omega = " << omega <<
", tau = " << tau << std::endl ;
453 if (basisType==coshBasis || basisType == sinhBasis) {
454 if (
verboseEval()>0) {std::cout <<
"RooGaussModel::analyticalIntegral(" <<
GetName() <<
") 8th form tau=" << tau << std::endl ; }
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());
464 if (basisType==linBasis) {
465 if (
verboseEval()>0) std::cout <<
"RooGaussModel::analyticalIntegral(" <<
GetName() <<
") 6th form tau=" << tau << std::endl ;
468 double f1 = std::exp(-umax*umax) - std::exp(-umin*umin);
473 double f2 = tmp1 - tmp2;
474 double f3 = xpmax*tmp1 - xpmin*tmp2;
476 double expc2 = std::exp(
c*
c);
480 (1 - 2*
c*
c)*expc2*f2 +
486 if (basisType==quadBasis) {
487 if (
verboseEval()>0) std::cout <<
"RooGaussModel::analyticalIntegral(" <<
GetName() <<
") 7th form tau=" << tau << std::endl ;
491 double tmpA1 = std::exp(-umax*umax);
492 double tmpA2 = std::exp(-umin*umin);
494 double f1 = tmpA1 - tmpA2;
495 double f2 = umax*tmpA1 - umin*tmpA2;
500 double f3 = tmpB1 - tmpB2;
501 double f4 = xpmax*tmpB1 - xpmin*tmpB2;
502 double f5 = xpmax*xpmax*tmpB1 - xpmin*xpmin*tmpB2;
504 double expc2 = std::exp(
c*
c);
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
520 std::complex<double> diff(2., 0.);
522 diff = evalCerf(_x,umin,
c);
523 diff -= evalCerf(_x,umax,
c);
527 diff *= std::complex<double>(1., _x);
528 diff *= tau / (1.+_x*_x);
536 return matchArgs(directVars,generateVars,
x) ? 1 : 0;
544 double xmin =
x.min();
545 double xmax =
x.max();
549 if (xgen<xmax && xgen>
xmin) {
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
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
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...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
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.
std::span< const double > at(RooAbsArg const *arg, RooAbsArg const *caller=nullptr)
std::span< double > output()
RooBatchCompute::Config config(RooAbsArg const *arg) const
Class RooGaussModel implements a RooResolutionModel that models a Gaussian distribution.
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
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.
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...
double analyticalIntegral(Int_t code, const char *rangeName) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
static std::complex< double > erfc(const std::complex< double > z)
complex erfc function
static std::complex< double > erf(const std::complex< double > z)
complex erf function
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
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.
This is the base class for the ROOT Random number generators.
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...
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)