38 x(
"x",
"Dependent",this,_x),
39 mean(
"mean",
"Mean",this,_mean),
40 sigma(
"sigma",
"Width",this,_sigma)
50 mean(
"mean",this,other.mean),
73constexpr double p1[5] = {0.4259894875,-0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253};
74constexpr double q1[5] = {1.0 ,-0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063};
76constexpr double p2[5] = {0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211};
77constexpr double q2[5] = {1.0 , 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714};
79constexpr double p3[5] = {0.1788544503, 0.09359161662,0.006325387654, 0.00006611667319,-0.000002031049101};
80constexpr double q3[5] = {1.0 , 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675};
82constexpr double p4[5] = {0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186};
83constexpr double q4[5] = {1.0 , 106.8615961, 337.6496214, 2016.712389, 1597.063511};
85constexpr double p5[5] = {1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910};
86constexpr double q5[5] = {1.0 , 156.9424537, 3745.310488, 9834.698876, 66924.28357};
88constexpr double p6[5] = {1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109};
89constexpr double q6[5] = {1.0 , 651.4101098, 56974.73333, 165917.4725, -2815759.939};
91constexpr double a1[3] = {0.04166666667,-0.01996527778, 0.02709538966};
92constexpr double a2[2] = {-1.845568670,-4.284640743};
94template<
class Tx,
class Tmean,
class Tsigma>
95void compute(
size_t batchSize,
96 double * __restrict
output,
97 Tx X, Tmean M, Tsigma
S)
99 const double NaN = std::nan(
"");
100 constexpr size_t block=256;
103 for (
size_t i=0; i<batchSize; i+=
block) {
104 const size_t stop = (i+
block < batchSize) ?
block : batchSize-i ;
106 for (
size_t j=0; j<stop; j++) {
107 v[j] = (X[i+j]-M[i+j]) /
S[i+j];
108 output[i+j] = (p2[0]+(p2[1]+(p2[2]+(p2[3]+p2[4]*
v[j])*
v[j])*
v[j])*
v[j]) /
109 (q2[0]+(q2[1]+(q2[2]+(q2[3]+q2[4]*
v[j])*
v[j])*
v[j])*
v[j]);
112 for (
size_t j=0; j<stop; j++) {
113 const bool mask =
S[i+j] > 0;
117 if (!mask)
v[j] =
NaN;
122 for (
size_t j=0; j<stop; j++) {
126 output[i+j] = (p3[0]+(p3[1]+(p3[2]+(p3[3]+p3[4]*
v[j])*
v[j])*
v[j])*
v[j]) /
127 (q3[0]+(q3[1]+(q3[2]+(q3[3]+q3[4]*
v[j])*
v[j])*
v[j])*
v[j]);
128 }
else if (
v[j] < 12) {
130 output[i+j] = u*u*(p4[0]+(p4[1]+(p4[2]+(p4[3]+p4[4]*u)*u)*u)*u) /
131 (q4[0]+(q4[1]+(q4[2]+(q4[3]+q4[4]*u)*u)*u)*u);
132 }
else if (
v[j] < 50) {
134 output[i+j] = u*u*(p5[0]+(p5[1]+(p5[2]+(p5[3]+p5[4]*u)*u)*u)*u) /
135 (q5[0]+(q5[1]+(q5[2]+(q5[3]+q5[4]*u)*u)*u)*u);
136 }
else if (
v[j] < 300) {
138 output[i+j] = u*u*(p6[0]+(p6[1]+(p6[2]+(p6[3]+p6[4]*u)*u)*u)*u) /
139 (q6[0]+(q6[1]+(q6[2]+(q6[3]+q6[4]*u)*u)*u)*u);
142 output[i+j] = u*u*(1 +(a2[0] +a2[1]*u)*u );
144 }
else if (
v[j] < -1) {
148 (p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*
v[j])*
v[j])*
v[j])*
v[j])/
149 (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*
v[j])*
v[j])*
v[j])*
v[j]);
152 if (u < 1
e-10)
output[i+j] = 0.0;
156 output[i+j] = 0.3989422803*(ue/
us)*(1+(a1[0]+(a1[1]+a1[2]*u)*u)*u);
181 const bool batchX = !xData.empty();
182 const bool batchMean = !meanData.empty();
183 const bool batchSigma = !sigmaData.empty();
185 if (!batchX && !batchMean && !batchSigma) {
188 batchSize =
findSize({ xData, meanData, sigmaData });
191 if (batchX && !batchMean && !batchSigma ) {
194 else if (!batchX && batchMean && !batchSigma ) {
197 else if (batchX && batchMean && !batchSigma ) {
200 else if (!batchX && !batchMean && batchSigma ) {
203 else if (batchX && !batchMean && batchSigma ) {
206 else if (!batchX && batchMean && batchSigma ) {
209 else if (batchX && batchMean && batchSigma ) {
210 compute(batchSize,
output.data(), xData, meanData, sigmaData);
219 if (
matchArgs(directVars,generateVars,
x))
return 1 ;
227 assert(1 == code); (
void)code;
typedef void((*Func_t)())
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize, const RooArgSet *const normSet=nullptr, Tag_t ownerTag=kUnspecified)
Make a batch and return a span pointing to the pdf-local memory.
Little adapter that gives a bracket operator to types that don't have one.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
BatchHelpers::BatchData _batchData
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Landau distribution p.d.f.
void generateEvent(Int_t code)
Interface for generation of an event using the algorithm corresponding to the specified code.
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Compute in batches.
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
A simple container to hold a batch of data values.
double min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
double max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
RooSpan< const typename T::value_type > getValBatch(std::size_t begin, std::size_t batchSize) const
Retrieve a batch of real or category data.
virtual Double_t Landau(Double_t mean=0, Double_t sigma=1)
Generate a random number following a Landau distribution with location parameter mu and scale paramet...
size_t findSize(std::vector< RooSpan< const double > > parameters)
This function returns the minimum size of the non-zero-sized batches.
RooArgSet S(const RooAbsArg &v1)
void checkRangeOfParameters(const RooAbsReal *callingClass, std::initializer_list< const RooAbsReal * > pars, double min=-std::numeric_limits< double >::max(), double max=std::numeric_limits< double >::max(), bool limitsInAllowedRange=false, std::string extraMessage="")
Check if the parameters have a range, and warn if the range extends below / above the set limits.
static constexpr double us
Double_t Landau(Double_t x, Double_t mpv=0, Double_t sigma=1, Bool_t norm=kFALSE)
The LANDAU function.
static void output(int code)