| File: | generators/evtgen/models/besiii/src/EvtDsToKKpi.cc |
| Warning: | line 192, column 7 Value stored to 'tmp1' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | // Model: EvtDsToKKpi |
| 2 | // This file is an amplitude model for Ds+ -> K- K+ pi+. |
| 3 | // The model is from the BESIII Collaboration in PRD 104, 012016 (2021). DOI: https://doi.org/10.1103/PhysRevD.104.012016 |
| 4 | // |
| 5 | // Permission to use these files in basf2 was generously granted by the BESIII Collaboration. |
| 6 | // |
| 7 | // Please cite the original reference for any public/published results where this model was used. |
| 8 | |
| 9 | #include <iomanip> |
| 10 | #include <cmath> |
| 11 | #include <string> |
| 12 | #include <EvtGenBase/EvtCPUtil.hh> |
| 13 | #include <EvtGenBase/EvtTensor4C.hh> |
| 14 | |
| 15 | #include <EvtGenBase/EvtPatches.hh> |
| 16 | #include <fstream> |
| 17 | #include <stdlib.h> |
| 18 | #include <EvtGenBase/EvtParticle.hh> |
| 19 | #include <EvtGenBase/EvtGenKine.hh> |
| 20 | #include <EvtGenBase/EvtPDL.hh> |
| 21 | #include <EvtGenBase/EvtReport.hh> |
| 22 | #include <EvtGenBase/EvtResonance.hh> |
| 23 | #include <EvtGenBase/EvtResonance2.hh> |
| 24 | #include <string> |
| 25 | #include <EvtGenBase/EvtConst.hh> |
| 26 | #include <EvtGenBase/EvtFlatte.hh> |
| 27 | #include <EvtGenBase/EvtDecayTable.hh> |
| 28 | |
| 29 | #include <generators/evtgen/EvtGenModelRegister.h> |
| 30 | #include <generators/evtgen/models/besiii/EvtDsToKKpi.h> |
| 31 | |
| 32 | namespace Belle2 { |
| 33 | |
| 34 | /** register the model in EvtGen */ |
| 35 | B2_EVTGEN_REGISTER_MODEL(EvtDsToKKpi)namespace { Belle2::EvtGenModelRegister::Factory<EvtDsToKKpi > EvtGenModelFactory_EvtDsToKKpi; }; |
| 36 | |
| 37 | EvtDsToKKpi::~EvtDsToKKpi() {} |
| 38 | |
| 39 | std::string EvtDsToKKpi::getName() |
| 40 | { |
| 41 | return "DsToKKpi"; |
| 42 | } |
| 43 | |
| 44 | EvtDecayBase* EvtDsToKKpi::clone() |
| 45 | { |
| 46 | return new EvtDsToKKpi; |
| 47 | } |
| 48 | |
| 49 | void EvtDsToKKpi::init() |
| 50 | { |
| 51 | checkNArg(0); |
| 52 | checkNDaug(3); |
| 53 | checkSpinParent(EvtSpinType::SCALAR); |
| 54 | checkSpinDaughter(0, EvtSpinType::SCALAR); |
| 55 | checkSpinDaughter(1, EvtSpinType::SCALAR); |
| 56 | checkSpinDaughter(2, EvtSpinType::SCALAR); |
| 57 | |
| 58 | |
| 59 | phi[0] = 0; |
| 60 | phi[1] = 6.1944e+00; |
| 61 | phi[2] = 4.7398e+00; |
| 62 | phi[3] = 2.9047e+00; |
| 63 | phi[4] = 1.0068e+00; |
| 64 | phi[5] = 5.8035e-01; |
| 65 | |
| 66 | rho[0] = 1; |
| 67 | rho[1] = 1.0963e+00; |
| 68 | rho[2] = 2.7818e+00; |
| 69 | rho[3] = 1.2570e+00; |
| 70 | rho[4] = 7.7351e-01; |
| 71 | rho[5] = 5.6277e-01; |
| 72 | |
| 73 | mass[0] = 8.9581e-01; |
| 74 | mass[1] = 1.019461e+00; |
| 75 | mass[2] = 0.919; |
| 76 | mass[3] = 1.4712e+00; |
| 77 | mass[4] = 1.7220e+00; |
| 78 | mass[5] = 1.3500e+00; |
| 79 | |
| 80 | width[0] = 4.7400e-02; |
| 81 | width[1] = 4.2660e-03; |
| 82 | width[2] = 0.272; |
| 83 | width[3] = 2.7000e-01; |
| 84 | width[4] = 1.3500e-01; |
| 85 | width[5] = 2.6500e-01; |
| 86 | |
| 87 | modetype[0] = 0; |
| 88 | modetype[1] = 1; |
| 89 | modetype[2] = 13; |
| 90 | modetype[3] = 3; |
| 91 | modetype[4] = 4; |
| 92 | modetype[5] = 5; |
| 93 | |
| 94 | |
| 95 | mD = 1.86484; |
| 96 | mDs = 1.96835; |
| 97 | rRes = 1.5; |
| 98 | rD = 3.0; |
| 99 | mkstr = 0.89555; |
| 100 | mk0 = 0.497611; |
| 101 | mass_Kaon = 0.49368; |
| 102 | mass_Pion = 0.13957; |
| 103 | mass_Pi0 = 0.1349768; |
| 104 | mass_EtaP = 0.95778; |
| 105 | mass_Eta = 0.547862; |
| 106 | math_pi = 3.1415926; |
| 107 | afRatio = 2.04835; //BABAR // afRatio= 1.23202; //BES2 |
| 108 | int GG[4][4] = { {1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1} }; |
| 109 | for (int i = 0; i < 4; i++) { |
| 110 | for (int j = 0; j < 4; j++) { |
| 111 | G[i][j] = GG[i][j]; |
| 112 | } |
| 113 | } |
| 114 | } |
| 115 | |
| 116 | void EvtDsToKKpi::initProbMax() |
| 117 | { |
| 118 | setProbMax(35638.0); |
| 119 | } |
| 120 | |
| 121 | void EvtDsToKKpi::decay(EvtParticle* p) |
| 122 | { |
| 123 | p->initializePhaseSpace(getNDaug(), getDaugs()); |
| 124 | EvtVector4R D1 = p->getDaug(0)->getP4(); |
| 125 | EvtVector4R D2 = p->getDaug(1)->getP4(); |
| 126 | EvtVector4R D3 = p->getDaug(2)->getP4(); |
| 127 | double P1[4], P2[4], P3[4]; |
| 128 | P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3); |
| 129 | P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3); |
| 130 | P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3); |
| 131 | |
| 132 | double value; |
| 133 | int g0[6] = {1, 1, 1, 1, 1, 1}; |
| 134 | int nstates = 6; |
| 135 | calEvaMy(P1, P2, P3, mass, width, rho, phi, g0, modetype, nstates, value); |
| 136 | setProb(value); |
| 137 | return ; |
| 138 | } |
| 139 | |
| 140 | void EvtDsToKKpi::MIP_LineShape(double sa, double pro[2]) |
| 141 | { |
| 142 | double m0 = mass[2], cKK = width[2]; |
| 143 | pro[0] = sqrt(1 / (((m0 * m0) - sa) * ((m0 * m0) - sa) + (cKK * m0 * sqrt(1 - 4 * mass_Kaon * mass_Kaon / sa)) * (cKK * m0 * sqrt( |
| 144 | 1 - 4 * mass_Kaon * mass_Kaon / sa)))); |
| 145 | pro[1] = 0; |
| 146 | } |
| 147 | void EvtDsToKKpi::calEvaMy(double* pKm, double* pKp, double* pPi, double* mass1, double* width1, double* amp, double* phase, |
| 148 | int* g0, int* modetype_param, int nstates, double& Result) // Renamed modetype to modetype_param to avoid shadowing member variable |
| 149 | { |
| 150 | double pF0_980[4], pPhi1020[4], pKstr[4], pD[4], p23[4]; |
| 151 | pKp[0] = sqrt(mass_Kaon * mass_Kaon + pKp[1] * pKp[1] + pKp[2] * pKp[2] + pKp[3] * pKp[3]); |
| 152 | pKm[0] = sqrt(mass_Kaon * mass_Kaon + pKm[1] * pKm[1] + pKm[2] * pKm[2] + pKm[3] * pKm[3]); |
| 153 | pPi[0] = sqrt(mass_Pion * mass_Pion + pPi[1] * pPi[1] + pPi[2] * pPi[2] + pPi[3] * pPi[3]); |
| 154 | for (int u = 0; u < 4; u++) { |
| 155 | pF0_980[u] = pKp[u] + pKm[u]; |
| 156 | pPhi1020[u] = pKp[u] + pKm[u]; |
| 157 | pKstr[u] = pKm[u] + pPi[u]; |
| 158 | p23[u] = pKp[u] + pPi[u]; |
| 159 | pD[u] = pKp[u] + pKm[u] + pPi[u]; |
| 160 | } |
| 161 | double cof[2], amp_tmp1[2], amp_tmp[2], amp_PDF[2], PDF[2]; |
| 162 | amp_PDF[0] = 0; |
| 163 | amp_PDF[1] = 0; |
| 164 | PDF[0] = 0; |
| 165 | PDF[1] = 0; |
| 166 | double temp_PDF, tmp1; |
| 167 | double pro[2], B[3]; |
| 168 | double t1kstr1[4], t1phi1020[4], t1D1[4], t1D2[4]; |
| 169 | |
| 170 | double sD, sF0_980, sF0_1710, sF0_1370, sKp, sKm, sPi, sKstr, sPhi1020, sKstr1430; |
| 171 | sF0_980 = SCADot(pF0_980, pF0_980); |
| 172 | sF0_1710 = sF0_980; |
| 173 | sF0_1370 = sF0_980; |
| 174 | sKstr = SCADot(pKstr, pKstr); |
| 175 | sD = SCADot(pD, pD); |
| 176 | sKstr1430 = sKstr; |
| 177 | sPhi1020 = SCADot(pPhi1020, pPhi1020); |
| 178 | sKp = SCADot(pKp, pKp); |
| 179 | sKm = SCADot(pKm, pKm); |
| 180 | sPi = SCADot(pPi, pPi); |
| 181 | |
| 182 | calt1(pKm, pPi, t1kstr1); |
| 183 | calt1(pKm, pKp, t1phi1020); |
| 184 | calt1(pKstr, pKp, t1D1); |
| 185 | calt1(pPhi1020, pPi, t1D2); |
| 186 | |
| 187 | for (int i = 0; i < nstates; i++) { |
| 188 | amp_tmp[0] = 0; |
| 189 | amp_tmp[1] = 0; |
| 190 | amp_tmp1[0] = 0; |
| 191 | amp_tmp1[1] = 0; |
| 192 | tmp1 = 0; |
Value stored to 'tmp1' is never read | |
| 193 | temp_PDF = 0; |
| 194 | cof[0] = amp[i] * cos(phase[i]); |
| 195 | cof[1] = amp[i] * sin(phase[i]); |
| 196 | if (modetype_param[i] == 13) { |
| 197 | //a0(980) and f0(980) mixture |
| 198 | double amp_a0[2] = {0}; |
| 199 | double sa0_980 = sF0_980; |
| 200 | MIP_LineShape(sa0_980, pro); |
| 201 | B[0] = barrier(0, sa0_980, sKp, sKm, rRes); |
| 202 | temp_PDF = 1; |
| 203 | tmp1 = temp_PDF * B[0]; |
| 204 | amp_a0[0] = tmp1 * pro[0]; |
| 205 | amp_a0[1] = tmp1 * pro[1]; |
| 206 | amp_tmp1[0] = amp_a0[0] ; |
| 207 | amp_tmp1[1] = amp_a0[1] ; |
| 208 | } else if (modetype_param[i] == 0) { |
| 209 | //K*(892) K+ |
| 210 | if (g0[i] == 0) { |
| 211 | pro[0] = 1; |
| 212 | pro[1] = 0; |
| 213 | } |
| 214 | bool neo = true; |
| 215 | if (neo) { |
| 216 | double sBC = SCADot(p23, p23); |
| 217 | if (g0[i] == 1) propagatorRBWNeoKstr892(mass1[i], width1[i], sKstr, sPi, sKm, rRes, 1, pro); |
| 218 | B[0] = barrierNeo(1, sKstr, sPi, sKm, rRes, mass1[i]); |
| 219 | B[1] = barrierNeoDs(1, sD, sKstr, sKp, rD, mDs, mass1[i]); |
| 220 | temp_PDF = (sBC - sPhi1020 + ((sD - sKp) * (sKm - sPi) / (sKstr))); |
| 221 | } else { |
| 222 | if (g0[i] == 1) propagatorRBW(mass1[i], width1[i], sKstr, sKm, sPi, rRes, 1, pro); |
| 223 | B[0] = barrier(1, sKstr, sKm, sPi, rRes); |
| 224 | B[1] = barrier(1, sD, sKstr, sKp, rD); |
| 225 | for (int m = 0; m < 4; m++) { |
| 226 | for (int j = 0; j < 4; j++) { |
| 227 | temp_PDF += t1D1[m] * t1kstr1[j] * G[m][j]; |
| 228 | } |
| 229 | } |
| 230 | } |
| 231 | tmp1 = temp_PDF * B[0] * B[1]; |
| 232 | amp_tmp1[0] = tmp1 * pro[0]; |
| 233 | amp_tmp1[1] = tmp1 * pro[1]; |
| 234 | } else if (modetype_param[i] == 1) { |
| 235 | if (g0[i] == 0) { |
| 236 | pro[0] = 1; |
| 237 | pro[1] = 0; |
| 238 | } |
| 239 | bool neo = true; |
| 240 | if (neo) { |
| 241 | if (g0[i] == 1) propagatorRBWNeo(mass1[i], width1[i], sPhi1020, sKm, sKp, rRes, 1, pro); |
| 242 | B[0] = barrierNeo(1, sPhi1020, sKp, sKm, rRes, mass1[i]); |
| 243 | B[1] = barrierNeoDs(1, sD, sPhi1020, sPi, rD, mDs, mass1[i]); |
| 244 | double sBC = SCADot(p23, p23); |
| 245 | temp_PDF = (sKstr - sBC + ((sD - sPi) * (sKp - sKm) / (sKstr))); |
| 246 | } else { |
| 247 | if (g0[i] == 1) propagatorRBW(mass1[i], width1[i], sPhi1020, sKm, sKp, rRes, 1, pro); |
| 248 | B[0] = barrier(1, sPhi1020, sKp, sKm, rRes); |
| 249 | B[1] = barrier(1, sD, sPhi1020, sPi, rD); |
| 250 | for (int m = 0; m < 4; m++) { |
| 251 | for (int j = 0; j < 4; j++) { |
| 252 | temp_PDF += t1D2[m] * t1phi1020[j] * G[m][j]; |
| 253 | } |
| 254 | } |
| 255 | } |
| 256 | |
| 257 | tmp1 = temp_PDF * B[0] * B[1]; |
| 258 | amp_tmp1[0] = tmp1 * pro[0]; |
| 259 | amp_tmp1[1] = tmp1 * pro[1]; |
| 260 | } else if (modetype_param[i] == 3) { |
| 261 | //Kstr1430 K S |
| 262 | double sKm2[2] = {sKm, mass_EtaP * mass_EtaP}; |
| 263 | double sPi2[2] = {sPi, mass_Kaon * mass_Kaon}; |
| 264 | propagatorKstr1430(mass1[i], sKstr1430, sKm2, sPi2, pro); |
| 265 | B[0] = barrier(0, sPhi1020, sKp, sKm, rRes); |
| 266 | tmp1 = 1 * B[0]; |
| 267 | amp_tmp1[0] = tmp1 * pro[0]; |
| 268 | amp_tmp1[1] = tmp1 * pro[1]; |
| 269 | |
| 270 | } else if (modetype_param[i] == 4) { |
| 271 | if (g0[i] == 1) propagatorRBWNeo(mass1[i], width1[i], sF0_1710, sKp, sKm, rRes, 0, pro); |
| 272 | if (g0[i] == 0) { |
| 273 | pro[0] = 1; |
| 274 | pro[1] = 0; |
| 275 | } |
| 276 | B[0] = barrier(0, sF0_1710, sKp, sKm, rRes); |
| 277 | temp_PDF = 1; |
| 278 | tmp1 = temp_PDF * B[0]; |
| 279 | amp_tmp1[0] = tmp1 * pro[0]; |
| 280 | amp_tmp1[1] = tmp1 * pro[1]; |
| 281 | } else if (modetype_param[i] == 5) { |
| 282 | //f0(1370) Pi+ S |
| 283 | if (g0[i] == 1) propagatorRBWNeo(mass1[i], width1[i], sF0_1370, sKp, sKm, rRes, 0, pro); |
| 284 | if (g0[i] == 0) { |
| 285 | pro[0] = 1; |
| 286 | pro[1] = 0; |
| 287 | } |
| 288 | B[0] = barrier(0, sF0_1370, sKp, sKm, rRes); |
| 289 | tmp1 = 1 * B[0]; |
| 290 | amp_tmp1[0] = tmp1 * pro[0]; |
| 291 | amp_tmp1[1] = tmp1 * pro[1]; |
| 292 | } |
| 293 | amp_tmp[0] = amp_tmp1[0]; |
| 294 | amp_tmp[1] = amp_tmp1[1]; |
| 295 | Com_Multi(amp_tmp, cof, amp_PDF); |
| 296 | PDF[0] += amp_PDF[0]; |
| 297 | PDF[1] += amp_PDF[1]; |
| 298 | } |
| 299 | |
| 300 | double value = PDF[0] * PDF[0] + PDF[1] * PDF[1]; |
| 301 | if (value <= 0) { |
| 302 | value = 1e-20; |
| 303 | } |
| 304 | Result = value; |
| 305 | } |
| 306 | |
| 307 | void EvtDsToKKpi::Com_Multi(double a1[2], double a2[2], double res[2]) |
| 308 | { |
| 309 | res[0] = a1[0] * a2[0] - a1[1] * a2[1]; |
| 310 | res[1] = a1[1] * a2[0] + a1[0] * a2[1]; |
| 311 | } |
| 312 | void EvtDsToKKpi::Com_Divide(double a1[2], double a2[2], double res[2]) |
| 313 | { |
| 314 | res[0] = (a1[0] * a2[0] + a1[1] * a2[1]) / (a2[0] * a2[0] + a2[1] * a2[1]); |
| 315 | res[1] = (a1[1] * a2[0] - a1[0] * a2[1]) / (a2[0] * a2[0] + a2[1] * a2[1]); |
| 316 | } |
| 317 | |
| 318 | double EvtDsToKKpi::SCADot(double a1[4], double a2[4]) |
| 319 | { |
| 320 | double _cal = 0.; |
| 321 | _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3]; |
| 322 | return _cal; |
| 323 | } |
| 324 | double EvtDsToKKpi::barrier(int l, double sa, double sb, double sc, double r) |
| 325 | { |
| 326 | double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb; |
| 327 | if (q < 0) q = 1e-16; |
| 328 | double z; |
| 329 | z = q * r * r; |
| 330 | double F = 0.0; |
| 331 | if (l == 0) F = 1; |
| 332 | if (l == 1) F = sqrt(2 * z / (1 + z)); |
| 333 | if (l == 2) F = sqrt(13 * z * z / (9 + 3 * z + z * z)); |
| 334 | return F; |
| 335 | } |
| 336 | double EvtDsToKKpi::barrierNeo(int l, double sa, double sb, double sc, double r, double mR) |
| 337 | { |
| 338 | double pAB = ((sa - sb - sc) * (sa - sb - sc) / 4.0 - (sb * sc)) / sa; |
| 339 | double pR = ((mR * mR - sb - sc) * (mR * mR - sb - sc) / 4.0 - (sb * sc)) / (mR * mR); |
| 340 | double zAB = pAB * r * r; |
| 341 | double zR = pR * r * r; |
| 342 | double F = 0.0; |
| 343 | if (l == 0) F = 1; |
| 344 | if (l == 1) F = sqrt((1 + zR) / (1 + zAB)); |
| 345 | if (l == 2) F = sqrt((9 + 3 * zAB + zAB * zAB) / (9 + 3 * zAB + zAB * zAB)); |
| 346 | return F; |
| 347 | } |
| 348 | double EvtDsToKKpi::barrierNeoDs(int l, double sa, double sb, double sc, double r, double mR, double mb) |
| 349 | { |
| 350 | double pAB = ((sa - sb - sc) * (sa - sb - sc) / 4.0 - (sb * sc)) / sa; |
| 351 | double pR = ((sa - mb * mb - sc) * (sa - mb * mb - sc) / 4.0 - (mb * mb * sc)) / (mR * mR); |
| 352 | double zAB = pAB * r * r; |
| 353 | double zR = pR * r * r; |
| 354 | double F = 0.0; |
| 355 | if (l == 0) F = 1; |
| 356 | if (l == 1) F = sqrt((1 + zR) / (1 + zAB)); |
| 357 | if (l == 2) F = sqrt((9 + 3 * zAB + zAB * zAB) / (9 + 3 * zAB + zAB * zAB)); |
| 358 | return F; |
| 359 | } |
| 360 | //------------------spin------------------------------------------- |
| 361 | void EvtDsToKKpi::calt1(double daug1[4], double daug2[4], double t1[4]) |
| 362 | { |
| 363 | double p, pq; |
| 364 | double pa[4], qa[4]; |
| 365 | for (int i = 0; i < 4; i++) { |
| 366 | pa[i] = daug1[i] + daug2[i]; |
| 367 | qa[i] = daug1[i] - daug2[i]; |
| 368 | } |
| 369 | p = SCADot(pa, pa); |
| 370 | pq = SCADot(pa, qa); |
| 371 | for (int i = 0; i < 4; i++) { |
| 372 | t1[i] = qa[i] - pq / p * pa[i]; |
| 373 | } |
| 374 | } |
| 375 | void EvtDsToKKpi::calt2(double daug1[4], double daug2[4], double t2[4][4]) |
| 376 | { |
| 377 | double p, r; |
| 378 | double pa[4], t1[4]; |
| 379 | calt1(daug1, daug2, t1); |
| 380 | r = SCADot(t1, t1); |
| 381 | for (int i = 0; i < 4; i++) { |
| 382 | pa[i] = daug1[i] + daug2[i]; |
| 383 | } |
| 384 | p = SCADot(pa, pa); |
| 385 | for (int i = 0; i < 4; i++) { |
| 386 | for (int j = 0; j < 4; j++) { |
| 387 | t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * (G[i][j] - pa[i] * pa[j] / p); |
| 388 | } |
| 389 | } |
| 390 | } |
| 391 | //-------------------prop-------------------------------------------- |
| 392 | void EvtDsToKKpi::propagator(double mass_param, double width_param, double sx, double prop[2]) |
| 393 | { |
| 394 | double a[2], b[2]; |
| 395 | a[0] = 1; |
| 396 | a[1] = 0; |
| 397 | b[0] = mass_param * mass_param - sx; |
| 398 | b[1] = -mass_param * width_param; |
| 399 | Com_Divide(a, b, prop); |
| 400 | } |
| 401 | double EvtDsToKKpi::wid(double mass_param, double sa, double sb, double sc, double r, int l) |
| 402 | { |
| 403 | double widm = 0.; |
| 404 | double q = 0.; |
| 405 | double q0 = 0.; |
| 406 | double sa0 = mass_param * mass_param; |
| 407 | double m = sqrt(sa); |
| 408 | q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb; |
| 409 | if (q < 0) q = 1e-16; |
| 410 | q0 = (sa0 + sb - sc) * (sa0 + sb - sc) / (4 * sa0) - sb; |
| 411 | if (q0 < 0) q0 = 1e-16; |
| 412 | double z = q * r * r; |
| 413 | double z0 = q0 * r * r; |
| 414 | double F = 0; |
| 415 | if (l == 0) F = 1; |
| 416 | if (l == 1) F = sqrt((1 + z0) / (1 + z)); |
| 417 | if (l == 2) F = sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z)); |
| 418 | double t = q / q0; |
| 419 | widm = pow(t, l + 0.5) * mass_param / m * F * F; |
| 420 | return widm; |
| 421 | } |
| 422 | |
| 423 | void EvtDsToKKpi::Flatte_rhoab(double sa, double sb, double sc, double rho_param[2]) |
| 424 | { |
| 425 | double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb; |
| 426 | if (q > 0) { |
| 427 | rho_param[0] = 2 * sqrt(q / sa); |
| 428 | rho_param[1] = 0; |
| 429 | } else if (q < 0) { |
| 430 | rho_param[0] = 0; |
| 431 | rho_param[1] = 2 * sqrt(-q / sa); |
| 432 | } |
| 433 | } |
| 434 | |
| 435 | void EvtDsToKKpi::propagatorFlatte(double mass_param, double width_param __attribute__((unused)), double sx, double* sb, double* sc, |
| 436 | double prop[2]) |
| 437 | { |
| 438 | double unit[2] = {1.0}; |
| 439 | double ci[2] = {0, 1}; |
| 440 | double rho1[2]; |
| 441 | Flatte_rhoab(sx, sb[0], sc[0], rho1); |
| 442 | double rho2[2]; |
| 443 | Flatte_rhoab(sx, sb[1], sc[1], rho2); |
| 444 | double g1_f980 = 0.165, g2_f980 = 0.69465; |
| 445 | double tmp1[2] = {g1_f980, 0}; |
| 446 | double tmp11[2]; |
| 447 | double tmp2[2] = {g2_f980, 0}; |
| 448 | double tmp22[2]; |
| 449 | Com_Multi(tmp1, rho1, tmp11); |
| 450 | Com_Multi(tmp2, rho2, tmp22); |
| 451 | double tmp3[2] = {tmp11[0] + tmp22[0], tmp11[1] + tmp22[1]}; |
| 452 | double tmp31[2]; |
| 453 | Com_Multi(tmp3, ci, tmp31); |
| 454 | double tmp4[2] = {mass_param* mass_param - sx - tmp31[0], -1.0 * tmp3[1]}; |
| 455 | Com_Divide(unit, tmp4, prop); |
| 456 | } |
| 457 | void EvtDsToKKpi::propagator980(double mass_param, double sx, double* sb, double* sc, double prop[2]) |
| 458 | { |
| 459 | double unit[2] = {1.0}; |
| 460 | double ci[2] = {0, 1}; |
| 461 | double rho1[2]; |
| 462 | Flatte_rhoab(sx, sb[0], sc[0], rho1); |
| 463 | double rho2[2]; |
| 464 | Flatte_rhoab(sx, sb[1], sc[1], rho2); |
| 465 | double gK_f980 = 0.69466, gPi_f980 = 0.165; |
| 466 | double tmp1[2] = {gK_f980, 0}; |
| 467 | double tmp11[2]; |
| 468 | double tmp2[2] = {gPi_f980, 0}; |
| 469 | double tmp22[2]; |
| 470 | Com_Multi(tmp1, rho1, tmp11); |
| 471 | Com_Multi(tmp2, rho2, tmp22); |
| 472 | double tmp3[2] = {tmp11[0] + tmp22[0], tmp11[1] + tmp22[1]}; |
| 473 | double tmp31[2]; |
| 474 | Com_Multi(tmp3, ci, tmp31); |
| 475 | double tmp4[2] = {mass_param* mass_param - sx - tmp31[0], -1.0 * tmp31[1]}; |
| 476 | Com_Divide(unit, tmp4, prop); |
| 477 | } |
| 478 | |
| 479 | void EvtDsToKKpi::propagatora0980(double mass_param, double sx, double* sb, double* sc, double prop[2]) |
| 480 | { |
| 481 | double unit[2] = {1.0}; |
| 482 | double ci[2] = {0, 1}; |
| 483 | double rho1[2]; |
| 484 | Flatte_rhoab(sx, sb[0], sc[0], rho1); |
| 485 | double rho2[2]; |
| 486 | Flatte_rhoab(sx, sb[1], sc[1], rho2); |
| 487 | double gKK_a980 = 0.892 * 0.341, gPiEta_a980 = 0.341; |
| 488 | double tmp1[2] = {gKK_a980, 0}; |
| 489 | double tmp11[2]; |
| 490 | double tmp2[2] = {gPiEta_a980, 0}; |
| 491 | double tmp22[2]; |
| 492 | Com_Multi(tmp1, rho1, tmp11); |
| 493 | Com_Multi(tmp2, rho2, tmp22); |
| 494 | double tmp3[2] = {tmp11[0] + tmp22[0], tmp11[1] + tmp22[1]}; |
| 495 | double tmp31[2]; |
| 496 | Com_Multi(tmp3, ci, tmp31); |
| 497 | double tmp4[2] = {mass_param* mass_param - sx - tmp31[0], -1.0 * tmp31[1]}; |
| 498 | Com_Divide(unit, tmp4, prop); |
| 499 | } |
| 500 | |
| 501 | void EvtDsToKKpi::propagatorKstr1430(double mass_param, double sx, double* sb, double* sc, double prop[2]) |
| 502 | { |
| 503 | double unit[2] = {1.0}; |
| 504 | double ci[2] = {0, 1}; |
| 505 | double rho1[2]; |
| 506 | Flatte_rhoab(sx, sb[0], sc[0], rho1); |
| 507 | double rho2[2]; |
| 508 | Flatte_rhoab(sx, sb[1], sc[1], rho2); |
| 509 | double gKPi_Kstr1430 = 0.2990, gEtaPK_Kstr1430 = 0.0529; |
| 510 | double tmp1[2] = {gKPi_Kstr1430, 0}; |
| 511 | double tmp11[2]; |
| 512 | double tmp2[2] = {gEtaPK_Kstr1430, 0}; |
| 513 | double tmp22[2]; |
| 514 | Com_Multi(tmp1, rho1, tmp11); |
| 515 | Com_Multi(tmp2, rho2, tmp22); |
| 516 | double tmp3[2] = {tmp11[0] + tmp22[0], tmp11[1] + tmp22[1]}; |
| 517 | double tmp31[2]; |
| 518 | Com_Multi(tmp3, ci, tmp31); |
| 519 | double tmp4[2] = {mass_param* mass_param - sx - tmp31[0], -1.0 * tmp31[1]}; |
| 520 | Com_Divide(unit, tmp4, prop); |
| 521 | } |
| 522 | |
| 523 | void EvtDsToKKpi::propagatorRBW(double mass_param, double width_param, double sa, double sb, double sc, double r, int l, |
| 524 | double prop[2]) |
| 525 | { |
| 526 | double a[2], b[2]; |
| 527 | a[0] = 1; |
| 528 | a[1] = 0; |
| 529 | b[0] = mass_param * mass_param - sa; |
| 530 | b[1] = -mass_param * width_param * wid(mass_param, sa, sb, sc, r, l); |
| 531 | Com_Divide(a, b, prop); |
| 532 | } |
| 533 | |
| 534 | void EvtDsToKKpi::propagatorRBWNeo(double mass_param, double width_param, double sa, double sb, double sc, |
| 535 | double r __attribute__((unused)), int l, double prop[2]) |
| 536 | { |
| 537 | double a[2], b[2]; |
| 538 | a[0] = 1; |
| 539 | a[1] = 0; |
| 540 | b[0] = mass_param * mass_param - sa; |
| 541 | double tmp1 = (sa - sb - sc); |
| 542 | double tmp2 = sb * sc; |
| 543 | double pAB = sqrt((tmp1 * tmp1 / 4.0 - tmp2) / sa); |
| 544 | double pR = sqrt(((mass_param * mass_param - sb - sc) * (mass_param * mass_param - sb - sc) / 4.0 - (sb * sc)) / |
| 545 | (mass_param * mass_param)); |
| 546 | double fR = sqrt(1.0 + 1.5 * 1.5 * pR * pR) / sqrt(1.0 + 1.5 * 1.5 * pAB * pAB); |
| 547 | double power = 1; |
| 548 | if (!l) { |
| 549 | power = 1; |
| 550 | fR = 1; |
| 551 | } else if (l == 1) { |
| 552 | power = 3; |
| 553 | } |
| 554 | double gammaAB = width_param * pow(pAB / pR, power) * (mass_param / sqrt(sa)) * fR * fR; |
| 555 | b[1] = -mass_param * gammaAB; |
| 556 | Com_Divide(a, b, prop); |
| 557 | } |
| 558 | |
| 559 | void EvtDsToKKpi::propagatorRBWNeoKstr892(double mass_param, double width_param, double sa, double sb, double sc, |
| 560 | double r __attribute__((unused)), int l, double prop[2]) |
| 561 | { |
| 562 | double a[2], b[2]; |
| 563 | a[0] = 1; |
| 564 | a[1] = 0; |
| 565 | b[0] = mass_param * mass_param - sa; |
| 566 | double tmp1 = (sa - sb - sc); |
| 567 | double tmp2 = sb * sc; |
| 568 | double pAB = sqrt((tmp1 * tmp1 / 4.0 - tmp2) / sa); |
| 569 | double pR = sqrt(((mass_param * mass_param - sb - sc) * (mass_param * mass_param - sb - sc) / 4.0 - (sb * sc)) / |
| 570 | (mass_param * mass_param)); |
| 571 | double fR = sqrt(1.0 + 1.5 * 1.5 * pR * pR) / sqrt(1.0 + 1.5 * 1.5 * pAB * pAB); |
| 572 | double power = 1; |
| 573 | if (!l) { |
| 574 | power = 1; |
| 575 | } else if (l == 1) { |
| 576 | power = 3; |
| 577 | } |
| 578 | double gammaAB = width_param * pow(pAB / pR, power) * (mass_param / sqrt(sa)) * fR * fR; |
| 579 | b[1] = -mass_param * gammaAB; |
| 580 | Com_Divide(a, b, prop); |
| 581 | } |
| 582 | |
| 583 | //------------GS---used by rho---------------------------- |
| 584 | double EvtDsToKKpi::h(double m, double q) |
| 585 | { |
| 586 | double h = 2 / math_pi * q / m * log((m + 2 * q) / (2 * mass_Pion)); |
| 587 | return h; |
| 588 | } |
| 589 | double EvtDsToKKpi::dh(double mass_param, double q0) |
| 590 | { |
| 591 | double dh = h(mass_param, q0) * (1.0 / (8 * q0 * q0) - 1.0 / (2 * mass_param * mass_param)) + 1.0 / |
| 592 | (2 * math_pi * mass_param * mass_param); |
| 593 | return dh; |
| 594 | } |
| 595 | double EvtDsToKKpi::f(double mass_param, double sx, double q0, double q) |
| 596 | { |
| 597 | double m = sqrt(sx); |
| 598 | double f = mass_param * mass_param / (pow(q0, 3)) * (q * q * (h(m, q) - h(mass_param, |
| 599 | q0)) + (mass_param * mass_param - sx) * q0 * q0 * dh(mass_param, q0)); |
| 600 | return f; |
| 601 | } |
| 602 | double EvtDsToKKpi::d(double mass_param, double q0) |
| 603 | { |
| 604 | double d = 3.0 / math_pi * mass_Pion * mass_Pion / (q0 * q0) * log((mass_param + 2 * q0) / (2 * mass_Pion)) + mass_param / |
| 605 | (2 * math_pi * q0) |
| 606 | - (mass_Pion * mass_Pion * mass_param) / (math_pi * pow(q0, 3)); |
| 607 | return d; |
| 608 | } |
| 609 | |
| 610 | //rho |
| 611 | void EvtDsToKKpi::propagatorGS(double mass_param, double width_param, double sa, double sb, double sc, double r, int l, |
| 612 | double prop[2]) |
| 613 | { |
| 614 | double a[2], b[2]; |
| 615 | double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb; |
| 616 | double sa0 = mass_param * mass_param; |
| 617 | double q0 = (sa0 + sb - sc) * (sa0 + sb - sc) / (4 * sa0) - sb; |
| 618 | if (q < 0) q = 1e-16; |
| 619 | if (q0 < 0) q0 = 1e-16; |
| 620 | q = sqrt(q); |
| 621 | q0 = sqrt(q0); |
| 622 | a[0] = 1 + d(mass_param, q0) * width_param / mass_param; |
| 623 | a[1] = 0; |
| 624 | b[0] = mass_param * mass_param - sa + width_param * f(mass_param, sa, q0, q); |
| 625 | b[1] = -mass_param * width_param * wid(mass_param, sa, sb, sc, r, l); |
| 626 | Com_Divide(a, b, prop); |
| 627 | } |
| 628 | |
| 629 | |
| 630 | |
| 631 | } // Belle2 namespace |