1 #include <goofit/PDFs/physics/TruthResolution_Aux.h>
5 __device__ fptype device_truth_resolution(fptype coshterm,
15 unsigned int * /*indices*/) {
18 ret += coshterm * cosh(ymixing * dtime);
19 ret += costerm * cos(xmixing * dtime);
20 ret -= 2 * sinhterm * sinh(ymixing * dtime);
22 * sin(xmixing * dtime); // Notice sign difference wrt to Mikhail's code, because I have AB* and he has A*B.
25 // printf("device_truth_resolution %f %f %f %f %f\n", coshterm, costerm, sinhterm, sinterm, dtime);
29 __device__ fptype device_truth_resolution_average_tau(
30 fptype A2, fptype B2, fptype ABr, fptype ABi, fptype xmixing, fptype ymixing, fptype tau) {
35 fptype averagetau = ((xmixing * xmixing + 1) * (ymixing * ymixing - 1)
36 * (((a * tau * (xmixing * xmixing - 1.) + 2. * b * tau * xmixing))
37 / ((xmixing * xmixing + 1.) * (xmixing * xmixing + 1.))
38 + (c * (-(tau * ymixing * ymixing) - tau) + d * (2. * tau) * ymixing)
39 / ((ymixing * ymixing - 1.) * (ymixing * ymixing - 1.))))
40 / ((ymixing * ymixing - 1) * (b * xmixing - a) + (xmixing * xmixing + 1) * (c - d * ymixing));
41 // printf("device avg tau: %.5g with A2: %.5g, B2: %.5g, ABr:%.5g, ABi:%.5g, x:%.5g, y:%.5g, tau:%.5g \n",
42 // averagetau, A2, B2, ABr, ABi, xmixing, ymixing, tau);
46 __device__ device_resfunction_ptr ptr_to_truth = device_truth_resolution;
47 __device__ device_calc_tau_fcn_ptr ptr_to_calc_tau = device_truth_resolution_average_tau;
49 TruthResolution::TruthResolution()
50 : MixingTimeResolution() {
51 GET_FUNCTION_ADDR(ptr_to_truth);
53 GET_FUNCTION_ADDR(ptr_to_calc_tau);
54 setCalcTauIdx(GooPdf::findFunctionIdx(host_fcn_ptr));
56 TruthResolution::~TruthResolution() = default;
58 fptype TruthResolution::normalisation(
59 fptype di1, fptype di2, fptype di3, fptype di4, fptype tau, fptype xmixing, fptype ymixing) const {
60 fptype timeIntegralOne = tau / (1 - ymixing * ymixing);
61 fptype timeIntegralTwo = tau / (1 + xmixing * xmixing);
62 fptype timeIntegralThr = ymixing * timeIntegralOne;
63 fptype timeIntegralFou = xmixing * timeIntegralTwo;
65 fptype ret = timeIntegralOne * (di1 + di2);
66 ret += timeIntegralTwo * (di1 - di2);
67 ret -= 2 * timeIntegralThr * di3;
68 ret -= 2 * timeIntegralFou * di4;