1 #include <goofit/PDFs/basic/KinLimitBWPdf.h>
2 #include <goofit/Variable.h>
6 __device__ fptype getMomentum(const fptype &mass, const fptype &pimass, const fptype &d0mass) {
10 double lambda = mass * mass - pimass * pimass - d0mass * d0mass;
12 lambda -= 4 * pimass * pimass * d0mass * d0mass;
17 return sqrt(0.5 * lambda / mass);
20 __device__ fptype bwFactor(const fptype &momentum) {
21 // 2.56 = 1.6^2, comes from radius for spin-1 particle
22 return 1 / sqrt(1.0 + 2.56 * momentum * momentum);
25 __device__ fptype device_KinLimitBW(fptype *evt, fptype *p, unsigned int *indices) {
26 fptype x = evt[RO_CACHE(indices[2 + RO_CACHE(indices[0])])];
27 fptype mean = RO_CACHE(p[RO_CACHE(indices[1])]);
28 fptype width = RO_CACHE(p[RO_CACHE(indices[2])]);
29 fptype d0mass = RO_CACHE(functorConstants[RO_CACHE(indices[3]) + 0]);
30 fptype pimass = RO_CACHE(functorConstants[RO_CACHE(indices[3]) + 1]);
35 fptype pUsingRealMass = getMomentum(mean, pimass, d0mass);
37 if(0 >= pUsingRealMass)
41 fptype pUsingX = getMomentum(x, pimass, d0mass);
42 fptype phspfactor = POW3(pUsingX / pUsingRealMass) * POW2(bwFactor(pUsingX) / bwFactor(pUsingRealMass));
43 fptype phspMassSq = POW2(mean - x * x);
44 fptype phspGammaSq = POW2(width * phspfactor);
46 fptype ret = (phspfactor * mean * width * width) / (phspMassSq + mean * phspGammaSq);
51 __device__ device_function_ptr ptr_to_KinLimitBW = device_KinLimitBW;
53 __host__ KinLimitBWPdf::KinLimitBWPdf(std::string n, Observable _x, Variable mean, Variable width)
55 registerParameter(mean);
56 registerParameter(width);
58 std::vector<unsigned int> pindices;
59 pindices.push_back(mean.getIndex());
60 pindices.push_back(width.getIndex());
61 pindices.push_back(registerConstants(2));
62 setMasses(1.8645, 0.13957);
63 GET_FUNCTION_ADDR(ptr_to_KinLimitBW);
67 __host__ void KinLimitBWPdf::setMasses(fptype bigM, fptype smallM) {
70 constants[1] = smallM;
71 MEMCPY_TO_SYMBOL(functorConstants, constants, 2 * sizeof(fptype), cIndex * sizeof(fptype), cudaMemcpyHostToDevice);