| File: | generators/kkmc/photospp/photosC.cc |
| Warning: | line 279, column 7 Value stored to 'PHINT' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | #include "Photos.h" |
| 2 | #include "forW-MEc.h" |
| 3 | #include "forZ-MEc.h" |
| 4 | #include "pairs.h" |
| 5 | #include "Log.h" |
| 6 | #include <cstdio> |
| 7 | #include <cmath> |
| 8 | #include <iostream> |
| 9 | #include "photosC.h" |
| 10 | #include "HEPEVT_struct.h" |
| 11 | #include "PhotosUtilities.h" |
| 12 | using std::cout; |
| 13 | using std::endl; |
| 14 | using std::max; |
| 15 | using namespace Photospp; |
| 16 | using namespace PhotosUtilities; |
| 17 | |
| 18 | namespace Photospp { |
| 19 | |
| 20 | // Instantiating structs declared in photosC.h |
| 21 | |
| 22 | struct HEPEVT hep; |
| 23 | struct HEPEVT pho; |
| 24 | |
| 25 | struct PHOCOP phocop; |
| 26 | struct PHNUM phnum; |
| 27 | struct PHOKEY phokey; |
| 28 | struct PHOSTA phosta; |
| 29 | struct PHOLUN pholun; |
| 30 | struct PHOPHS phophs; |
| 31 | struct TOFROM tofrom; |
| 32 | struct PHOPRO phopro; |
| 33 | struct PHOREST phorest; |
| 34 | struct PHWT phwt; |
| 35 | struct PHOCORWT phocorwt; |
| 36 | struct PHOMOM phomom; |
| 37 | struct PHOCMS phocms; |
| 38 | struct PHOEXP phoexp; |
| 39 | struct PHLUPY phlupy; |
| 40 | struct DARKR darkr; |
| 41 | /** Logical function used deep inside algorithm to check if emitted |
| 42 | particles are to emit. For mother it blocks the vertex, |
| 43 | but for daughters individually: bad sisters will not prevent electron to emit. |
| 44 | top quark has further exception method. */ |
| 45 | bool F(int m, int i) |
| 46 | { |
| 47 | return Photos::IPHQRK_setQarknoEmission(0, i) && (i <= 41 || i > 100) |
| 48 | && i != 21 |
| 49 | && i != 2101 && i != 3101 && i != 3201 |
| 50 | && i != 1103 && i != 2103 && i != 2203 |
| 51 | && i != 3103 && i != 3203 && i != 3303; |
| 52 | } |
| 53 | |
| 54 | |
| 55 | // --- can be used with VARIANT A. For B use PHINT1 or 2 -------------- |
| 56 | //---------------------------------------------------------------------- |
| 57 | // |
| 58 | // PHINT: PHotos universal INTerference correction weight |
| 59 | // |
| 60 | // Purpose: calculates correction weight as expressed by |
| 61 | // formula (17) from CPC 79 (1994), 291. |
| 62 | // |
| 63 | // Input Parameters: Common /PHOEVT/, with photon added. |
| 64 | // |
| 65 | // Output Parameters: correction weight |
| 66 | // |
| 67 | // Author(s): Z. Was, P.Golonka Created at: 19/01/05 |
| 68 | // Last Update: 23/06/13 |
| 69 | // |
| 70 | //---------------------------------------------------------------------- |
| 71 | |
| 72 | double PHINT(int IDUM) |
| 73 | { |
| 74 | |
| 75 | double PHINT2; |
| 76 | double EPS1[4], EPS2[4], PH[4], PL[4]; |
| 77 | static int i = 1; |
| 78 | int K, L; |
| 79 | // DOUBLE PRECISION EMU,MCHREN,BETA,COSTHG,MPASQR,XPH, XC1, XC2 |
| 80 | double XNUM1, XNUM2, XDENO, XC1, XC2; |
| 81 | |
| 82 | // REAL*8 PHOCHA |
| 83 | //-- |
| 84 | |
| 85 | // Calculate polarimetric vector: ph, eps1, eps2 are orthogonal |
| 86 | |
| 87 | for (K = 1; K <= 4; K++) { |
| 88 | PH[K - i] = pho.phep[pho.nhep - i][K - i]; |
| 89 | EPS2[K - i] = 1.0; |
| 90 | } |
| 91 | |
| 92 | |
| 93 | PHOEPS(PH, EPS2, EPS1); |
| 94 | PHOEPS(PH, EPS1, EPS2); |
| 95 | |
| 96 | |
| 97 | XNUM1 = 0.0; |
| 98 | XNUM2 = 0.0; |
| 99 | XDENO = 0.0; |
| 100 | |
| 101 | for (K = pho.jdahep[1 - i][1 - i]; K <= pho.nhep - i; K++) { //! or jdahep[1-i][2-i] |
| 102 | |
| 103 | // momenta of charged particle in PL |
| 104 | |
| 105 | for (L = 1; L <= 4; L++) PL[L - i] = pho.phep[K - i][L - i]; |
| 106 | |
| 107 | // scalar products: epsilon*p/k*p |
| 108 | |
| 109 | XC1 = - PHOCHA(pho.idhep[K - i]) * |
| 110 | (PL[1 - i] * EPS1[1 - i] + PL[2 - i] * EPS1[2 - i] + PL[3 - i] * EPS1[3 - i]) / |
| 111 | (PH[4 - i] * PL[4 - i] - PH[1 - i] * PL[1 - i] - PH[2 - i] * PL[2 - i] - PH[3 - i] * PL[3 - i]); |
| 112 | |
| 113 | XC2 = - PHOCHA(pho.idhep[K - i]) * |
| 114 | (PL[1 - i] * EPS2[1 - i] + PL[2 - i] * EPS2[2 - i] + PL[3 - i] * EPS2[3 - i]) / |
| 115 | (PH[4 - i] * PL[4 - i] - PH[1 - i] * PL[1 - i] - PH[2 - i] * PL[2 - i] - PH[3 - i] * PL[3 - i]); |
| 116 | |
| 117 | |
| 118 | // accumulate the currents |
| 119 | XNUM1 = XNUM1 + XC1; |
| 120 | XNUM2 = XNUM2 + XC2; |
| 121 | |
| 122 | XDENO = XDENO + XC1 * XC1 + XC2 * XC2; |
| 123 | } |
| 124 | |
| 125 | PHINT2 = (XNUM1 * XNUM1 + XNUM2 * XNUM2) / XDENO; |
| 126 | return PHINT2; |
| 127 | |
| 128 | } |
| 129 | |
| 130 | |
| 131 | |
| 132 | //---------------------------------------------------------------------- |
| 133 | // |
| 134 | // PHINT: PHotos INTerference (Old version kept for tests only. |
| 135 | // |
| 136 | // Purpose: Calculates interference between emission of photons from |
| 137 | // different possible chaged daughters stored in |
| 138 | // the HEP common /PHOEVT/. |
| 139 | // |
| 140 | // Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/ |
| 141 | // |
| 142 | // |
| 143 | // Output Parameters: |
| 144 | // |
| 145 | // |
| 146 | // Author(s): Z. Was, Created at: 10/08/93 |
| 147 | // Last Update: 15/03/99 |
| 148 | // |
| 149 | //---------------------------------------------------------------------- |
| 150 | |
| 151 | double PHINT1(int IDUM) |
| 152 | { |
| 153 | |
| 154 | double PHINT; |
| 155 | |
| 156 | /* |
| 157 | DOUBLE PRECISION MCHSQR,MNESQR |
| 158 | REAL*8 PNEUTR |
| 159 | COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5) |
| 160 | DOUBLE PRECISION COSTHG,SINTHG |
| 161 | REAL*8 XPHMAX,XPHOTO |
| 162 | COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG |
| 163 | |
| 164 | */ |
| 165 | double MPASQR, XX, BETA; |
| 166 | bool IFINT; |
| 167 | int K, IDENT; |
| 168 | double& COSTHG = phophs.costhg; |
| 169 | double& XPHOTO = phophs.xphoto; |
| 170 | //double *PNEUTR = phomom.pneutr; // unused variable |
| 171 | double& MCHSQR = phomom.mchsqr; |
| 172 | double& MNESQR = phomom.mnesqr; |
| 173 | |
| 174 | static int i = 1; |
| 175 | IDENT = pho.nhep; |
| 176 | // |
| 177 | for (K = pho.jdahep[1 - i][2 - i]; K >= pho.jdahep[1 - i][1 - i]; K--) { |
| 178 | if (pho.idhep[K - i] != 22) { |
| 179 | IDENT = K; |
| 180 | break; |
| 181 | } |
| 182 | } |
| 183 | |
| 184 | // check if there is a photon |
| 185 | IFINT = pho.nhep > IDENT; |
| 186 | // check if it is two body + gammas reaction |
| 187 | IFINT = IFINT && (IDENT - pho.jdahep[1 - i][1 - i]) == 1; |
| 188 | // check if two body was particle antiparticle |
| 189 | IFINT = IFINT && pho.idhep[pho.jdahep[1 - i][1 - i] - i] == -pho.idhep[IDENT - i]; |
| 190 | // check if particles were charged |
| 191 | IFINT = IFINT && PHOCHA(pho.idhep[IDENT - i]) != 0; |
| 192 | // calculates interference weight contribution |
| 193 | if (IFINT) { |
| 194 | MPASQR = pho.phep[1 - i][5 - i] * pho.phep[1 - i][5 - i]; |
| 195 | XX = 4.0 * MCHSQR / MPASQR * (1.0 - XPHOTO) / (1.0 - XPHOTO + (MCHSQR - MNESQR) / MPASQR) / (1.0 - XPHOTO + |
| 196 | (MCHSQR - MNESQR) / MPASQR); |
| 197 | BETA = sqrt(1.0 - XX); |
| 198 | PHINT = 2.0 / (1.0 + COSTHG * COSTHG * BETA * BETA); |
| 199 | } else { |
| 200 | PHINT = 1.0; |
| 201 | } |
| 202 | |
| 203 | return PHINT; |
| 204 | } |
| 205 | |
| 206 | |
| 207 | //---------------------------------------------------------------------- |
| 208 | // |
| 209 | // PHINT: PHotos INTerference |
| 210 | // |
| 211 | // Purpose: Calculates interference between emission of photons from |
| 212 | // different possible chaged daughters stored in |
| 213 | // the HEP common /PHOEVT/. |
| 214 | // |
| 215 | // Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/ |
| 216 | // |
| 217 | // |
| 218 | // Output Parameters: |
| 219 | // |
| 220 | // |
| 221 | // Author(s): Z. Was, Created at: 10/08/93 |
| 222 | // Last Update: |
| 223 | // |
| 224 | //---------------------------------------------------------------------- |
| 225 | |
| 226 | double PHINT2(int IDUM) |
| 227 | { |
| 228 | |
| 229 | |
| 230 | /* |
| 231 | DOUBLE PRECISION MCHSQR,MNESQR |
| 232 | REAL*8 PNEUTR |
| 233 | COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5) |
| 234 | DOUBLE PRECISION COSTHG,SINTHG |
| 235 | REAL*8 XPHMAX,XPHOTO |
| 236 | COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG |
| 237 | */ |
| 238 | double& COSTHG = phophs.costhg; |
| 239 | double& XPHOTO = phophs.xphoto; |
| 240 | double& MCHSQR = phomom.mchsqr; |
| 241 | double& MNESQR = phomom.mnesqr; |
| 242 | |
| 243 | |
| 244 | double MPASQR, XX, BETA, pq1[4], pq2[4], pphot[4]; |
| 245 | double SS, PP2, PP, E1, E2, q1, q2, costhe, PHINT; |
| 246 | bool IFINT; |
| 247 | int K, k, IDENT; |
| 248 | static int i = 1; |
| 249 | IDENT = pho.nhep; |
| 250 | // |
| 251 | for (K = pho.jdahep[1 - i][2 - i]; K >= pho.jdahep[1 - i][1 - i]; K--) { |
| 252 | if (pho.idhep[K - i] != 22) { |
| 253 | IDENT = K; |
| 254 | break; |
| 255 | } |
| 256 | } |
| 257 | |
| 258 | // check if there is a photon |
| 259 | IFINT = pho.nhep > IDENT; |
| 260 | // check if it is two body + gammas reaction |
| 261 | IFINT = IFINT && (IDENT - pho.jdahep[1 - i][1 - i]) == 1; |
| 262 | // check if two body was particle antiparticle (we improve on it ! |
| 263 | // IFINT= IFINT.AND.pho.idhep(JDAPHO(1,1)).EQ.-pho.idhep(IDENT) |
| 264 | // check if particles were charged |
| 265 | IFINT = IFINT && fabs(PHOCHA(pho.idhep[IDENT - i])) > 0.01; |
| 266 | // check if they have both charge |
| 267 | IFINT = IFINT && fabs(PHOCHA(pho.idhep[pho.jdahep[1 - i][1 - i] - i])) > 0.01; |
| 268 | // calculates interference weight contribution |
| 269 | if (IFINT) { |
| 270 | MPASQR = pho.phep[1 - i][5 - i] * pho.phep[1 - i][5 - i]; |
| 271 | XX = 4.0 * MCHSQR / MPASQR * (1.0 - XPHOTO) / pow(1. - XPHOTO + (MCHSQR - MNESQR) / MPASQR, 2); |
| 272 | BETA = sqrt(1.0 - XX); |
| 273 | PHINT = 2.0 / (1.0 + COSTHG * COSTHG * BETA * BETA); |
| 274 | SS = MPASQR * (1.0 - XPHOTO); |
| 275 | PP2 = ((SS - MCHSQR - MNESQR) * (SS - MCHSQR - MNESQR) - 4 * MCHSQR * MNESQR) / SS / 4; |
| 276 | PP = sqrt(PP2); |
| 277 | E1 = sqrt(PP2 + MCHSQR); |
| 278 | E2 = sqrt(PP2 + MNESQR); |
| 279 | PHINT = (E1 + E2) * (E1 + E2) / ((E2 + COSTHG * PP) * (E2 + COSTHG * PP) + (E1 - COSTHG * PP) * (E1 - COSTHG * PP)); |
Value stored to 'PHINT' is never read | |
| 280 | // return PHINT; |
| 281 | // |
| 282 | q1 = PHOCHA(pho.idhep[pho.jdahep[1 - i][1 - i] - i]); |
| 283 | q2 = PHOCHA(pho.idhep[IDENT - i]); |
| 284 | for (k = 1; k <= 4; k++) { |
| 285 | pq1[k - i] = pho.phep[pho.jdahep[1 - i][1 - i] - i][k - i]; |
| 286 | pq2[k - i] = pho.phep[pho.jdahep[1 - i][1 - i] + 1 - i][k - i]; |
| 287 | pphot[k - i] = pho.phep[pho.nhep - i][k - i]; |
| 288 | } |
| 289 | costhe = (pphot[1 - i] * pq1[1 - i] + pphot[2 - i] * pq1[2 - i] + pphot[3 - i] * pq1[3 - i]); |
| 290 | costhe = costhe / sqrt(pq1[1 - i] * pq1[1 - i] + pq1[2 - i] * pq1[2 - i] + pq1[3 - i] * pq1[3 - i]); |
| 291 | costhe = costhe / sqrt(pphot[1 - i] * pphot[1 - i] + pphot[2 - i] * pphot[2 - i] + pphot[3 - i] * pphot[3 - i]); |
| 292 | // |
| 293 | // --- this IF checks whether JDAPHO(1,1) was MCH or MNE. |
| 294 | // --- COSTHG angle (and in-generation variables) may be better choice |
| 295 | // --- than costhe. note that in the formulae below amplitudes were |
| 296 | // --- multiplied by (E2+COSTHG*PP)*(E1-COSTHG*PP). |
| 297 | if (COSTHG * costhe > 0) { |
| 298 | |
| 299 | PHINT = pow(q1 * (E2 + COSTHG * PP) - q2 * (E1 - COSTHG * PP), |
| 300 | 2) / (q1 * q1 * (E2 + COSTHG * PP) * (E2 + COSTHG * PP) + q2 * q2 * (E1 - COSTHG * PP) * (E1 - COSTHG * PP)); |
| 301 | } else { |
| 302 | |
| 303 | PHINT = pow(q1 * (E1 - COSTHG * PP) - q2 * (E2 + COSTHG * PP), |
| 304 | 2) / (q1 * q1 * (E1 - COSTHG * PP) * (E1 - COSTHG * PP) + q2 * q2 * (E2 + COSTHG * PP) * (E2 + COSTHG * PP)); |
| 305 | } |
| 306 | } else { |
| 307 | PHINT = 1.0; |
| 308 | } |
| 309 | return PHINT; |
| 310 | } |
| 311 | |
| 312 | |
| 313 | //***************************************************************** |
| 314 | //***************************************************************** |
| 315 | //***************************************************************** |
| 316 | // beginning of the class of methods reading from PH_HEPEVT |
| 317 | //***************************************************************** |
| 318 | //***************************************************************** |
| 319 | //***************************************************************** |
| 320 | |
| 321 | |
| 322 | //---------------------------------------------------------------------- |
| 323 | // |
| 324 | // PHOTOS: PHOton radiation in decays event DuMP routine |
| 325 | // |
| 326 | // Purpose: Print event record. |
| 327 | // |
| 328 | // Input Parameters: Common /PH_HEPEVT/ |
| 329 | // |
| 330 | // Output Parameters: None |
| 331 | // |
| 332 | // Author(s): B. van Eijk Created at: 05/06/90 |
| 333 | // Last Update: 20/06/13 |
| 334 | // |
| 335 | //---------------------------------------------------------------------- |
| 336 | void PHODMP() |
| 337 | { |
| 338 | |
| 339 | double SUMVEC[5]; |
| 340 | int I, J; |
| 341 | static int i = 1; |
| 342 | const char eq80[81] = "================================================================================"; |
| 343 | const char X29[30] = " "; |
| 344 | const char X23[24 ] = " "; |
| 345 | const char X1[2] = " "; |
| 346 | const char X2[3] = " "; |
| 347 | const char X3[4] = " "; |
| 348 | const char X4[5] = " "; |
| 349 | const char X6[7] = " "; |
| 350 | const char X7[8] = " "; |
| 351 | FILE* PHLUN = stdoutstdout; |
| 352 | |
| 353 | for (I = 0; I < 5; I++) SUMVEC[I] = 0.0; |
| 354 | //-- |
| 355 | //-- Print event number... |
| 356 | fprintf(PHLUN, "%s", eq80); |
| 357 | fprintf(PHLUN, "%s Event No.: %10i\n", X29, hep.nevhep); |
| 358 | fprintf(PHLUN, "%s Particle Parameters\n", X6); |
| 359 | fprintf(PHLUN, "%s Nr %s Type %s Parent(s) %s Daughter(s) %s Px %s Py %s Pz %s E %s Inv. M.\n", X1, X3, X3, X2, X6, X7, X7, X7, X4); |
| 360 | for (I = 1; I <= hep.nhep; I++) { |
| 361 | //-- |
| 362 | //-- For 'stable particle' calculate vector momentum sum |
| 363 | if (hep.jdahep[I - i][1 - i] == 0) { |
| 364 | for (J = 1; J <= 4; J++) { |
| 365 | SUMVEC[J - i] = SUMVEC[J - i] + hep.phep[I - i][J - i]; |
| 366 | } |
| 367 | if (hep.jmohep[I - i][2 - i] == 0) { |
| 368 | fprintf(PHLUN, "%4i %7i %s %4i %s Stable %9.2f %9.2f %9.2f %9.2f %9.2f\n", I, hep.idhep[I - i], X3, hep.jmohep[I - i][1 - i], X7, |
| 369 | hep.phep[I - i][1 - i], hep.phep[I - i][2 - i], hep.phep[I - i][3 - i], hep.phep[I - i][4 - i], hep.phep[I - i][5 - i]); |
| 370 | } else { |
| 371 | fprintf(PHLUN, "%4i %7i %4i - %4i %s Stable %9.2f %9.2f %9.2f %9.2f %9.2f\n", I, hep.idhep[I - i], hep.jmohep[I - i][1 - i], |
| 372 | hep.jmohep[I - i][2 - i], X4, hep.phep[I - i][1 - i], hep.phep[I - i][2 - i], hep.phep[I - i][3 - i], hep.phep[I - i][4 - i], |
| 373 | hep.phep[I - i][5 - i]); |
| 374 | } |
| 375 | } else { |
| 376 | if (hep.jmohep[I - i][2 - i] == 0) { |
| 377 | fprintf(PHLUN, "%4i %7i %s %4i %s %4i - %4i %9.2f %9.2f %9.2f %9.2f %9.2f\n", I, hep.idhep[I - i], X3, hep.jmohep[I - i][1 - i], |
| 378 | X2, hep.jdahep[I - i][1 - i], hep.jdahep[I - i][2 - i], hep.phep[I - i][1 - i], hep.phep[I - i][2 - i], hep.phep[I - i][3 - i], |
| 379 | hep.phep[I - i][4 - i], hep.phep[I - i][5 - i]); |
| 380 | } else { |
| 381 | fprintf(PHLUN, "%4i %7i %4i - %4i %4i - %4i %9.2f %9.2f %9.2f %9.2f %9.2f\n", I, hep.idhep[I - i], hep.jmohep[I - i][1 - i], |
| 382 | hep.jmohep[I - i][2 - i], hep.jdahep[I - i][1 - i], hep.jdahep[I - i][2 - i], hep.phep[I - i][1 - i], hep.phep[I - i][2 - i], |
| 383 | hep.phep[I - i][3 - i], hep.phep[I - i][4 - i], hep.phep[I - i][5 - i]); |
| 384 | } |
| 385 | } |
| 386 | } |
| 387 | SUMVEC[5 - i] = sqrt(SUMVEC[4 - i] * SUMVEC[4 - i] - SUMVEC[1 - i] * SUMVEC[1 - i] - SUMVEC[2 - i] * SUMVEC[2 - i] - SUMVEC[3 - i] * |
| 388 | SUMVEC[3 - i]); |
| 389 | fprintf(PHLUN, "%s Vector Sum: %9.2f %9.2f %9.2f %9.2f %9.2f\n", X23, SUMVEC[1 - i], SUMVEC[2 - i], SUMVEC[3 - i], SUMVEC[4 - i], |
| 390 | SUMVEC[5 - i]); |
| 391 | |
| 392 | fprintf(PHLUN, " Vertices:\n"); |
| 393 | fprintf(PHLUN, " Nr Type Vx Vy Vz T\n"); |
| 394 | for (I = 1; I <= hep.nhep; I++) { |
| 395 | fprintf(PHLUN, " %2i %7.2f %7.2f %7.2f %7.2f\n", I, hep.vhep[I - i][1 - i], hep.vhep[I - i][2 - i], hep.vhep[I - i][3 - i], |
| 396 | hep.vhep[I - i][4 - i]); |
| 397 | } |
| 398 | |
| 399 | |
| 400 | |
| 401 | |
| 402 | |
| 403 | // 9030 FORMAT(1H ,I4,I7,3X,I4,9X,'Stable',2X,5F9.2) |
| 404 | //"%4i %7i %s %4i %s Stable %s %9.2f %9.2f %9.2f %9.2f %9.2f " X3,9X,X2 |
| 405 | |
| 406 | // 9050 FORMAT(1H ,I4,I7,3X,I4,6X,I4,' - ',I4,5F9.2) |
| 407 | //"%4i %7i %s %4i %s %4i - %4i %9.2f %9.2f %9.2f %9.2f %9.2f " X3,X6 |
| 408 | |
| 409 | |
| 410 | |
| 411 | |
| 412 | //"%4i %7i %4i - %4i %s Stable %s %9.2f %9.2f %9.2f %9.2f %9.2f " X5,X2 |
| 413 | |
| 414 | |
| 415 | //9060 FORMAT(1H ,I4,I7,I4,' - ',I4,2X,I4,' - ',I4,5F9.2) |
| 416 | //"%4i %7i %4i - %4i %s %4i - %4i %9.2f %9.2f %9.2f %9.2f %9.2f " X2, |
| 417 | } |
| 418 | |
| 419 | |
| 420 | |
| 421 | //---------------------------------------------------------------------- |
| 422 | // |
| 423 | // PHLUPAB: debugging tool |
| 424 | // |
| 425 | // Purpose: NONE, eventually may printout content of the |
| 426 | // /PH_HEPEVT/ common |
| 427 | // |
| 428 | // Input Parameters: Common /PH_HEPEVT/ and /PHNUM/ |
| 429 | // latter may have number of the event. |
| 430 | // |
| 431 | // Output Parameters: None |
| 432 | // |
| 433 | // Author(s): Z. Was Created at: 30/05/93 |
| 434 | // Last Update: 20/06/13 |
| 435 | // |
| 436 | //---------------------------------------------------------------------- |
| 437 | |
| 438 | void PHLUPAB(int IPOINT) |
| 439 | { |
| 440 | char name[12] = "/PH_HEPEVT/"; |
| 441 | int I, J; |
| 442 | static int IPOIN0 = -5; |
| 443 | static int i = 1; |
| 444 | double SUM[5]; |
| 445 | FILE* PHLUN = stdoutstdout; |
| 446 | |
| 447 | if (IPOIN0 < 0) { |
| 448 | IPOIN0 = 400000; // ! maximal no-print point |
| 449 | phlupy.ipoin = IPOIN0; |
| 450 | phlupy.ipoinm = 400001; // ! minimal no-print point |
| 451 | } |
| 452 | |
| 453 | if (IPOINT <= phlupy.ipoinm || IPOINT >= phlupy.ipoin) return; |
| 454 | if ((int)phnum.iev < 1000) { |
| 455 | for (I = 1; I <= 5; I++) SUM[I - i] = 0.0; |
| 456 | |
| 457 | fprintf(PHLUN, "EVENT NR= %i WE ARE TESTING %s at IPOINT=%i \n", (int)phnum.iev, name, IPOINT); |
| 458 | fprintf(PHLUN, " ID p_x p_y p_z E m ID-MO_DA1 ID-MO_DA2\n"); |
| 459 | I = 1; |
| 460 | fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", hep.idhep[I - i], hep.phep[I - i][1 - i], hep.phep[I - i][2 - i], |
| 461 | hep.phep[I - i][3 - i], hep.phep[I - i][4 - i], hep.phep[I - i][5 - i], hep.jdahep[I - i][1 - i], hep.jdahep[I - i][2 - i]); |
| 462 | I = 2; |
| 463 | fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", hep.idhep[I - i], hep.phep[I - i][1 - i], hep.phep[I - i][2 - i], |
| 464 | hep.phep[I - i][3 - i], hep.phep[I - i][4 - i], hep.phep[I - i][5 - i], hep.jdahep[I - i][1 - i], hep.jdahep[I - i][2 - i]); |
| 465 | fprintf(PHLUN, " \n"); |
| 466 | for (I = 3; I <= hep.nhep; I++) { |
| 467 | fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", hep.idhep[I - i], hep.phep[I - i][1 - i], hep.phep[I - i][2 - i], |
| 468 | hep.phep[I - i][3 - i], hep.phep[I - i][4 - i], hep.phep[I - i][5 - i], hep.jmohep[I - i][1 - i], hep.jmohep[I - i][2 - i]); |
| 469 | for (J = 1; J <= 4; J++) SUM[J - i] = SUM[J - i] + hep.phep[I - i][J - i]; |
| 470 | } |
| 471 | |
| 472 | |
| 473 | SUM[5 - i] = sqrt(fabs(SUM[4 - i] * SUM[4 - i] - SUM[1 - i] * SUM[1 - i] - SUM[2 - i] * SUM[2 - i] - SUM[3 - i] * SUM[3 - i])); |
| 474 | fprintf(PHLUN, " SUM %14.9f %14.9f %14.9f %14.9f %14.9f\n", SUM[1 - i], SUM[2 - i], SUM[3 - i], SUM[4 - i], SUM[5 - i]); |
| 475 | |
| 476 | } |
| 477 | |
| 478 | |
| 479 | // 10 FORMAT(1X,' ID ','p_x ','p_y ','p_z ', |
| 480 | //$ 'E ','m ', |
| 481 | //$ 'ID-MO_DA1','ID-MO DA2' ) |
| 482 | // 20 FORMAT(1X,I4,5(F14.9),2I9) |
| 483 | //"%i4 %14.9f %14.9f %14.9f %14.9f %i9 i9" |
| 484 | // 30 FORMAT(1X,' SUM',5(F14.9)) |
| 485 | } |
| 486 | |
| 487 | |
| 488 | |
| 489 | |
| 490 | |
| 491 | |
| 492 | |
| 493 | |
| 494 | |
| 495 | //---------------------------------------------------------------------- |
| 496 | // |
| 497 | // PHLUPA: debugging tool |
| 498 | // |
| 499 | // Purpose: NONE, eventually may printout content of the |
| 500 | // /PHOEVT/ common |
| 501 | // |
| 502 | // Input Parameters: Common /PHOEVT/ and /PHNUM/ |
| 503 | // latter may have number of the event. |
| 504 | // |
| 505 | // Output Parameters: None |
| 506 | // |
| 507 | // Author(s): Z. Was Created at: 30/05/93 |
| 508 | // Last Update: 21/06/13 |
| 509 | // |
| 510 | //---------------------------------------------------------------------- |
| 511 | |
| 512 | void PHLUPA(int IPOINT) |
| 513 | { |
| 514 | char name[9] = "/PHOEVT/"; |
| 515 | int I, J; |
| 516 | static int IPOIN0 = -5; |
| 517 | static int i = 1; |
| 518 | double SUM[5]; |
| 519 | FILE* PHLUN = stdoutstdout; |
| 520 | |
| 521 | if (IPOIN0 < 0) { |
| 522 | IPOIN0 = 400000; // ! maximal no-print point |
| 523 | phlupy.ipoin = IPOIN0; |
| 524 | phlupy.ipoinm = 400001; // ! minimal no-print point |
| 525 | } |
| 526 | |
| 527 | if (IPOINT <= phlupy.ipoinm || IPOINT >= phlupy.ipoin) return; |
| 528 | if ((int)phnum.iev < 1000) { |
| 529 | for (I = 1; I <= 5; I++) SUM[I - i] = 0.0; |
| 530 | |
| 531 | fprintf(PHLUN, "EVENT NR= %i WE ARE TESTING %s at IPOINT=%i \n", (int)phnum.iev, name, IPOINT); |
| 532 | fprintf(PHLUN, " ID p_x p_y p_z E m ID-MO_DA1 ID-MO_DA2\n"); |
| 533 | I = 1; |
| 534 | fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", pho.idhep[I - i], pho.phep[I - i][1 - i], pho.phep[I - i][2 - i], |
| 535 | pho.phep[I - i][3 - i], pho.phep[I - i][4 - i], pho.phep[I - i][5 - i], pho.jdahep[I - i][1 - i], pho.jdahep[I - i][2 - i]); |
| 536 | I = 2; |
| 537 | fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", pho.idhep[I - i], pho.phep[I - i][1 - i], pho.phep[I - i][2 - i], |
| 538 | pho.phep[I - i][3 - i], pho.phep[I - i][4 - i], pho.phep[I - i][5 - i], pho.jdahep[I - i][1 - i], pho.jdahep[I - i][2 - i]); |
| 539 | fprintf(PHLUN, " \n"); |
| 540 | for (I = 3; I <= pho.nhep; I++) { |
| 541 | fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", pho.idhep[I - i], pho.phep[I - i][1 - i], pho.phep[I - i][2 - i], |
| 542 | pho.phep[I - i][3 - i], pho.phep[I - i][4 - i], pho.phep[I - i][5 - i], pho.jmohep[I - i][1 - i], pho.jmohep[I - i][2 - i]); |
| 543 | for (J = 1; J <= 4; J++) SUM[J - i] = SUM[J - i] + pho.phep[I - i][J - i]; |
| 544 | } |
| 545 | |
| 546 | |
| 547 | SUM[5 - i] = sqrt(fabs(SUM[4 - i] * SUM[4 - i] - SUM[1 - i] * SUM[1 - i] - SUM[2 - i] * SUM[2 - i] - SUM[3 - i] * SUM[3 - i])); |
| 548 | fprintf(PHLUN, " SUM %14.9f %14.9f %14.9f %14.9f %14.9f\n", SUM[1 - i], SUM[2 - i], SUM[3 - i], SUM[4 - i], SUM[5 - i]); |
| 549 | |
| 550 | } |
| 551 | |
| 552 | |
| 553 | // 10 FORMAT(1X,' ID ','p_x ','p_y ','p_z ', |
| 554 | //$ 'E ','m ', |
| 555 | //$ 'ID-MO_DA1','ID-MO DA2' ) |
| 556 | // 20 FORMAT(1X,I4,5(F14.9),2I9) |
| 557 | //"%4i %14.9f %14.9f %14.9f %14.9f %9i %9i" |
| 558 | // 30 FORMAT(1X,' SUM',5(F14.9)) |
| 559 | } |
| 560 | |
| 561 | |
| 562 | void PHOtoRF() |
| 563 | { |
| 564 | |
| 565 | |
| 566 | // COMMON /PH_TOFROM/ QQ[4],XM,th1,fi1 |
| 567 | double PP[4], RR[4]; |
| 568 | |
| 569 | int K, L; |
| 570 | static int i = 1; |
| 571 | |
| 572 | for (K = 1; K <= 4; K++) { |
| 573 | tofrom.QQ[K - i] = 0.0; |
| 574 | } |
| 575 | for (L = hep.jdahep[hep.jmohep[hep.nhep - i][1 - i] - i][1 - i]; L <= hep.jdahep[hep.jmohep[hep.nhep - i][1 - i] - i][2 - i]; L++) { |
| 576 | for (K = 1; K <= 4; K++) { |
| 577 | tofrom.QQ[K - i] = tofrom.QQ[K - i] + hep.phep[L - i][K - i]; |
| 578 | } |
| 579 | } |
| 580 | tofrom.XM = tofrom.QQ[4 - i] * tofrom.QQ[4 - i] - tofrom.QQ[3 - i] * tofrom.QQ[3 - i] - tofrom.QQ[2 - i] * tofrom.QQ[2 - i] - |
| 581 | tofrom.QQ[1 - i] * tofrom.QQ[1 - i]; |
| 582 | if (tofrom.XM > 0.0) tofrom.XM = sqrt(tofrom.XM); |
| 583 | if (tofrom.XM <= 0.0) return; |
| 584 | |
| 585 | for (L = 1; L <= hep.nhep; L++) { |
| 586 | for (K = 1; K <= 4; K++) { |
| 587 | PP[K - i] = hep.phep[L - i][K - i]; |
| 588 | } |
| 589 | bostdq(1, tofrom.QQ, PP, RR); |
| 590 | for (K = 1; K <= 4; K++) { |
| 591 | hep.phep[L - i][K - i] = RR[K - i]; |
| 592 | } |
| 593 | } |
| 594 | |
| 595 | tofrom.fi1 = 0.0; |
| 596 | tofrom.th1 = 0.0; |
| 597 | if (fabs(hep.phep[1 - i][1 - i]) + fabs(hep.phep[1 - i][2 - i]) > 0.0) tofrom.fi1 = PHOAN1(hep.phep[1 - i][1 - i], |
| 598 | hep.phep[1 - i][2 - i]); |
| 599 | if (fabs(hep.phep[1 - i][1 - i]) + fabs(hep.phep[1 - i][2 - i]) + fabs(hep.phep[1 - i][3 - i]) > 0.0) |
| 600 | tofrom.th1 = PHOAN2(hep.phep[1 - i][3 - i], |
| 601 | sqrt(hep.phep[1 - i][1 - i] * hep.phep[1 - i][1 - i] + hep.phep[1 - i][2 - i] * hep.phep[1 - i][2 - i])); |
| 602 | |
| 603 | for (L = 1; L <= hep.nhep; L++) { |
| 604 | for (K = 1; K <= 4; K++) { |
| 605 | RR[K - i] = hep.phep[L - i][K - i]; |
| 606 | } |
| 607 | |
| 608 | PHORO3(-tofrom.fi1, RR); |
| 609 | PHORO2(-tofrom.th1, RR); |
| 610 | for (K = 1; K <= 4; K++) { |
| 611 | hep.phep[L - i][K - i] = RR[K - i]; |
| 612 | } |
| 613 | } |
| 614 | |
| 615 | return; |
| 616 | } |
| 617 | |
| 618 | void PHOtoLAB() |
| 619 | { |
| 620 | |
| 621 | // // REAL*8 QQ(4),XM,th1,fi1 |
| 622 | // COMMON /PH_TOFROM/ QQ,XM,th1,fi1 |
| 623 | double PP[4], RR[4]; |
| 624 | int K, L; |
| 625 | static int i = 1; |
| 626 | |
| 627 | if (tofrom.XM <= 0.0) return; |
| 628 | |
| 629 | |
| 630 | for (L = 1; L <= hep.nhep; L++) { |
| 631 | for (K = 1; K <= 4; K++) { |
| 632 | PP[K - i] = hep.phep[L - i][K - i]; |
| 633 | } |
| 634 | |
| 635 | PHORO2(tofrom.th1, PP); |
| 636 | PHORO3(tofrom.fi1, PP); |
| 637 | bostdq(-1, tofrom.QQ, PP, RR); |
| 638 | |
| 639 | for (K = 1; K <= 4; K++) { |
| 640 | hep.phep[L - i][K - i] = RR[K - i]; |
| 641 | } |
| 642 | |
| 643 | // same for vertex |
| 644 | |
| 645 | for (K = 1; K <= 4; K++) { |
| 646 | PP[K - i] = hep.vhep[L - i][K - i]; |
| 647 | } |
| 648 | |
| 649 | PHORO2(tofrom.th1, PP); |
| 650 | PHORO3(tofrom.fi1, PP); |
| 651 | bostdq(-1, tofrom.QQ, PP, RR); |
| 652 | |
| 653 | for (K = 1; K <= 4; K++) { |
| 654 | hep.vhep[L - i][K - i] = RR[K - i]; |
| 655 | } |
| 656 | } |
| 657 | return; |
| 658 | } |
| 659 | void PHOcorrectDARK(int IDaction) |
| 660 | { |
| 661 | int K, L; |
| 662 | static int i = 1; |
| 663 | // if(darkr.ifspecial==1) |
| 664 | { |
| 665 | for (L = 1; L <= hep.nhep; L++) { |
| 666 | if (hep.idhep[L - i] == darkr.IDspecial) { |
| 667 | //PHODMP(); |
| 668 | hep.jdahep[L - i][1 - i] = L - i + 1 + 1; |
| 669 | hep.jdahep[L - i][2 - i] = L - i + 2 + 1; |
| 670 | hep.jmohep[L - i + 1][1 - i] = L - i + 1; |
| 671 | hep.jmohep[L - i + 2][1 - i] = L - i + 1; |
| 672 | |
| 673 | hep.jmohep[L - i - 2][2 - i] = 0; |
| 674 | hep.jmohep[L - i - 1][2 - i] = 0; |
| 675 | hep.jmohep[L - i][2 - i] = 0; |
| 676 | |
| 677 | hep.jmohep[L - i + 1][2 - i] = 0; |
| 678 | hep.jmohep[L - i + 2][2 - i] = 0; |
| 679 | // cout << "adres prosze= " << hep.jmohep[L-i][1-i] << endl; |
| 680 | hep.jdahep[ hep.jmohep[L - i][1 - i] - i ][2 - i] = hep.jdahep[ hep.jmohep[L - i][1 - i] - i ][2 - i] - 2; |
| 681 | //PHODMP(); |
| 682 | } |
| 683 | } |
| 684 | } |
| 685 | } |
| 686 | |
| 687 | |
| 688 | |
| 689 | |
| 690 | |
| 691 | |
| 692 | // 2) GENERAL INTERFACE: |
| 693 | // PHOTOS_GET |
| 694 | // PHOTOS_MAKE |
| 695 | |
| 696 | |
| 697 | // COMMONS: |
| 698 | // NAME USED IN SECT. # OF OC// Comment |
| 699 | // PHOQED 1) 2) 3 Flags whether emisson to be gen. |
| 700 | // PHOLUN 1) 4) 6 Output device number |
| 701 | // PHOCOP 1) 3) 4 photon coupling & min energy |
| 702 | // PHPICO 1) 3) 4) 5 PI & 2*PI |
| 703 | // PHSEED 1) 4) 3 RN seed |
| 704 | // PHOSTA 1) 4) 3 Status information |
| 705 | // PHOKEY 1) 2) 3) 7 Keys for nonstandard application |
| 706 | // PHOVER 1) 1 Version info for outside |
| 707 | // HEPEVT 2) 2 PDG common |
| 708 | // PH_HEPEVT2) 8 PDG common internal |
| 709 | // PHOEVT 2) 3) 10 PDG branch |
| 710 | // PHOIF 2) 3) 2 emission flags for PDG branch |
| 711 | // PHOMOM 3) 5 param of char-neutr system |
| 712 | // PHOPHS 3) 5 photon momentum parameters |
| 713 | // PHOPRO 3) 4 var. for photon rep. (in branch) |
| 714 | // PHOCMS 2) 3 parameters of boost to branch CMS |
| 715 | // PHNUM 4) 1 event number from outside |
| 716 | //---------------------------------------------------------------------- |
| 717 | |
| 718 | |
| 719 | //---------------------------------------------------------------------- |
| 720 | // |
| 721 | // PHOTOS_MAKE: General search routine |
| 722 | // |
| 723 | // Purpose: Search through the /PH_HEPEVT/ standard HEP common, sta- |
| 724 | // rting from the IPPAR-th particle. Whenevr branching |
| 725 | // point is found routine PHTYPE(IP) is called. |
| 726 | // Finally if calls on PHTYPE(IP) modified entries, common |
| 727 | // /PH_HEPEVT/ is ordered. |
| 728 | // |
| 729 | // Input Parameter: IPPAR: Pointer to decaying particle in |
| 730 | // /PH_HEPEVT/ and the common itself, |
| 731 | // |
| 732 | // Output Parameters: Common /PH_HEPEVT/, either with or without |
| 733 | // new particles added. |
| 734 | // |
| 735 | // Author(s): Z. Was, B. van Eijk Created at: 26/11/89 |
| 736 | // Last Update: 30/08/93 |
| 737 | // |
| 738 | //---------------------------------------------------------------------- |
| 739 | |
| 740 | void PHOTOS_MAKE_C(int IPARR) |
| 741 | { |
| 742 | static int i = 1; |
| 743 | int IPPAR, I, J, NLAST, MOTHER; |
| 744 | //-- |
| 745 | PHLUPAB(3); |
| 746 | |
| 747 | // write(*,*) 'at poczatek' |
| 748 | // PHODMP(); |
| 749 | IPPAR = abs(IPARR); |
| 750 | //-- Store pointers for cascade treatement... |
| 751 | NLAST = hep.nhep; |
| 752 | |
| 753 | |
| 754 | //-- |
| 755 | //-- Check decay multiplicity and minimum of correctness.. |
| 756 | if ((hep.jdahep[IPPAR - i][1 - i] == 0) || (hep.jmohep[hep.jdahep[IPPAR - i][1 - i] - i][1 - i] != IPPAR)) return; |
| 757 | |
| 758 | PHOtoRF(); |
| 759 | |
| 760 | // write(*,*) 'at przygotowany' |
| 761 | // PHODMP(); |
| 762 | |
| 763 | //-- |
| 764 | //-- single branch mode |
| 765 | //-- IPPAR is original position where the program was called |
| 766 | |
| 767 | //-- let-s do generation |
| 768 | PHTYPE(IPPAR); |
| 769 | |
| 770 | |
| 771 | //-- rearrange /PH_HEPEVT/ for added particles. |
| 772 | //-- at present this may be not needed as information |
| 773 | //-- is set at HepMC level. |
| 774 | if (hep.nhep > NLAST) { |
| 775 | for (I = NLAST + 1; I <= hep.nhep; I++) { |
| 776 | //-- |
| 777 | //-- Photon mother and vertex... |
| 778 | MOTHER = hep.jmohep[I - i][1 - i]; |
| 779 | hep.jdahep[MOTHER - i][2 - i] = I; |
| 780 | |
| 781 | // TP: This copies vertex position from mothers to daughters |
| 782 | // overwriting daughters position (if set) |
| 783 | //for( J=1;J<=4;J++){ |
| 784 | // hep.vhep[I-i][J-i]=hep.vhep[I-1-i][J-i]; |
| 785 | //} |
| 786 | |
| 787 | } |
| 788 | } |
| 789 | // write(*,*) 'at po dzialaniu ' |
| 790 | // PHODMP(); |
| 791 | PHOtoLAB(); |
| 792 | PHOcorrectDARK(1); |
| 793 | // write(*,*) 'at koniec' |
| 794 | // PHODMP(); |
| 795 | return; |
| 796 | } |
| 797 | |
| 798 | |
| 799 | |
| 800 | |
| 801 | //---------------------------------------------------------------------- |
| 802 | // |
| 803 | // PHCORK: corrects kinmatics of subbranch needed if host program |
| 804 | // produces events with the shaky momentum conservation |
| 805 | // |
| 806 | // Input Parameters: Common /PHOEVT/, MODCOR |
| 807 | // MODCOR >0 type of action |
| 808 | // =1 no action |
| 809 | // =2 corrects energy from mass |
| 810 | // =3 corrects mass from energy |
| 811 | // =4 corrects energy from mass for |
| 812 | // particles up to .4 GeV mass, |
| 813 | // for heavier ones corrects mass, |
| 814 | // =5 most complete correct also of mother |
| 815 | // often necessary for exponentiation. |
| 816 | // =0 execution mode |
| 817 | // |
| 818 | // Output Parameters: corrected /PHOEVT/ |
| 819 | // |
| 820 | // Author(s): P.Golonka, Z. Was Created at: 01/02/99 |
| 821 | // Modified : 07/07/13 |
| 822 | //---------------------------------------------------------------------- |
| 823 | |
| 824 | void PHCORK(int MODCOR) |
| 825 | { |
| 826 | |
| 827 | double M, P2, PX, PY, PZ, E, EN, XMS; |
| 828 | int I, K; |
| 829 | FILE* PHLUN = stdoutstdout; |
| 830 | |
| 831 | |
| 832 | static int MODOP = 0; |
| 833 | static int IPRINT = 0; |
| 834 | static double MCUT = 0.4; |
| 835 | static int i = 1; |
| 836 | |
| 837 | if (MODCOR != 0) { |
| 838 | // INITIALIZATION |
| 839 | MODOP = MODCOR; |
| 840 | |
| 841 | fprintf(PHLUN, "Message from PHCORK(MODCOR):: initialization\n"); |
| 842 | if (MODOP == 1) fprintf(PHLUN, "MODOP=1 -- no corrections on event: DEFAULT\n"); |
| 843 | else if (MODOP == 2) fprintf(PHLUN, "MODOP=2 -- corrects Energy from mass\n"); |
| 844 | else if (MODOP == 3) fprintf(PHLUN, "MODOP=3 -- corrects mass from Energy\n"); |
| 845 | else if (MODOP == 4) { |
| 846 | fprintf(PHLUN, "MODOP=4 -- corrects Energy from mass to Mcut\n"); |
| 847 | fprintf(PHLUN, " and mass from energy above Mcut\n"); |
| 848 | fprintf(PHLUN, " Mcut=%6.3f GeV", MCUT); |
| 849 | } else if (MODOP == 5) fprintf(PHLUN, "MODOP=5 -- corrects Energy from mass+flow\n"); |
| 850 | |
| 851 | else { |
| 852 | fprintf(PHLUN, "PHCORK wrong MODCOR=%4i\n", MODCOR); |
| 853 | exit(-1); |
| 854 | } |
| 855 | return; |
| 856 | } |
| 857 | |
| 858 | if (MODOP == 0 && MODCOR == 0) { |
| 859 | fprintf(PHLUN, "PHCORK lack of initialization\n"); |
| 860 | exit(-1); |
| 861 | } |
| 862 | |
| 863 | // execution mode |
| 864 | // ============== |
| 865 | // ============== |
| 866 | |
| 867 | |
| 868 | PX = 0.0; |
| 869 | PY = 0.0; |
| 870 | PZ = 0.0; |
| 871 | E = 0.0; |
| 872 | |
| 873 | if (MODOP == 1) { |
| 874 | // ----------------------- |
| 875 | // In this case we do nothing |
| 876 | return; |
| 877 | } else if (MODOP == 2) { |
| 878 | // ----------------------- |
| 879 | // lets loop thru all daughters and correct their energies |
| 880 | // according to E^2=p^2+m^2 |
| 881 | |
| 882 | for (I = 3; I <= pho.nhep; I++) { |
| 883 | |
| 884 | PX = PX + pho.phep[I - i][1 - i]; |
| 885 | PY = PY + pho.phep[I - i][2 - i]; |
| 886 | PZ = PZ + pho.phep[I - i][3 - i]; |
| 887 | |
| 888 | P2 = pho.phep[I - i][1 - i] * pho.phep[I - i][1 - i] + pho.phep[I - i][2 - i] * pho.phep[I - i][2 - i] + pho.phep[I - i][3 - i] * |
| 889 | pho.phep[I - i][3 - i]; |
| 890 | |
| 891 | EN = sqrt(pho.phep[I - i][5 - i] * pho.phep[I - i][5 - i] + P2); |
| 892 | |
| 893 | if (IPRINT == 1)fprintf(PHLUN, "CORRECTING ENERGY OF %6i: %14.9f => %14.9f\n", I, pho.phep[I - i][4 - i], EN); |
| 894 | |
| 895 | pho.phep[I - i][4 - i] = EN; |
| 896 | E = E + pho.phep[I - i][4 - i]; |
| 897 | |
| 898 | } |
| 899 | } |
| 900 | |
| 901 | else if (MODOP == 5) { |
| 902 | // ----------------------- |
| 903 | //C lets loop thru all daughters and correct their energies |
| 904 | //C according to E^2=p^2+m^2 |
| 905 | |
| 906 | for (I = 3; I <= pho.nhep; I++) { |
| 907 | PX = PX + pho.phep[I - i][1 - i]; |
| 908 | PY = PY + pho.phep[I - i][2 - i]; |
| 909 | PZ = PZ + pho.phep[I - i][3 - i]; |
| 910 | |
| 911 | P2 = pho.phep[I - i][1 - i] * pho.phep[I - i][1 - i] + pho.phep[I - i][2 - i] * pho.phep[I - i][2 - i] + pho.phep[I - i][3 - i] * |
| 912 | pho.phep[I - i][3 - i]; |
| 913 | |
| 914 | EN = sqrt(pho.phep[I - i][5 - i] * pho.phep[I - i][5 - i] + P2); |
| 915 | |
| 916 | if (IPRINT == 1)fprintf(PHLUN, "CORRECTING ENERGY OF %6i: %14.9f => %14.9f\n", I, pho.phep[I - i][4 - i], EN); |
| 917 | |
| 918 | pho.phep[I - i][4 - i] = EN; |
| 919 | E = E + pho.phep[I - i][4 - i]; |
| 920 | |
| 921 | } |
| 922 | for (K = 1; K <= 4; K++) { |
| 923 | pho.phep[1 - i][K - i] = 0.0; |
| 924 | for (I = 3; I <= pho.nhep; I++) { |
| 925 | pho.phep[1 - i][K - i] = pho.phep[1 - i][K - i] + pho.phep[I - i][K - i]; |
| 926 | } |
| 927 | } |
| 928 | XMS = sqrt(pho.phep[1 - i][4 - i] * pho.phep[1 - i][4 - i] - pho.phep[1 - i][3 - i] * pho.phep[1 - i][3 - i] - pho.phep[1 - i][2 - |
| 929 | i] * pho.phep[1 - i][2 - i] - pho.phep[1 - i][1 - i] * pho.phep[1 - i][1 - i]); |
| 930 | pho.phep[1 - i][5 - i] = XMS; |
| 931 | } else if (MODOP == 3) { |
| 932 | // ----------------------- |
| 933 | |
| 934 | // lets loop thru all daughters and correct their masses |
| 935 | // according to E^2=p^2+m^2 |
| 936 | |
| 937 | for (I = 3; I <= pho.nhep; I++) { |
| 938 | |
| 939 | PX = PX + pho.phep[I - i][1 - i]; |
| 940 | PY = PY + pho.phep[I - i][2 - i]; |
| 941 | PZ = PZ + pho.phep[I - i][3 - i]; |
| 942 | E = E + pho.phep[I - i][4 - i]; |
| 943 | |
| 944 | P2 = pho.phep[I - i][1 - i] * pho.phep[I - i][1 - i] + pho.phep[I - i][2 - i] * pho.phep[I - i][2 - i] + pho.phep[I - i][3 - i] * |
| 945 | pho.phep[I - i][3 - i]; |
| 946 | |
| 947 | M = sqrt(fabs(pho.phep[I - i][4 - i] * pho.phep[I - i][4 - i] - P2)); |
| 948 | |
| 949 | if (IPRINT == 1) fprintf(PHLUN, "CORRECTING MASS OF %6i: %14.9f => %14.9f\n", I, pho.phep[I - i][5 - i], M); |
| 950 | |
| 951 | pho.phep[I - i][5 - i] = M; |
| 952 | |
| 953 | } |
| 954 | |
| 955 | } else if (MODOP == 4) { |
| 956 | // ----------------------- |
| 957 | |
| 958 | // lets loop thru all daughters and correct their masses |
| 959 | // or energies according to E^2=p^2+m^2 |
| 960 | |
| 961 | for (I = 3; I <= pho.nhep; I++) { |
| 962 | |
| 963 | PX = PX + pho.phep[I - i][1 - i]; |
| 964 | PY = PY + pho.phep[I - i][2 - i]; |
| 965 | PZ = PZ + pho.phep[I - i][3 - i]; |
| 966 | P2 = pho.phep[I - i][1 - i] * pho.phep[I - i][1 - i] + pho.phep[I - i][2 - i] * pho.phep[I - i][2 - i] + pho.phep[I - i][3 - i] * |
| 967 | pho.phep[I - i][3 - i]; |
| 968 | M = sqrt(fabs(pho.phep[I - i][4 - i] * pho.phep[I - i][4 - i] - P2)); |
| 969 | |
| 970 | |
| 971 | if (M > MCUT) { |
| 972 | if (IPRINT == 1) fprintf(PHLUN, "CORRECTING MASS OF %6i: %14.9f => %14.9f\n", I, pho.phep[I - i][5 - i], M); |
| 973 | pho.phep[I - i][5 - i] = M; |
| 974 | E = E + pho.phep[I - i][4 - i]; |
| 975 | } else { |
| 976 | |
| 977 | EN = sqrt(pho.phep[I - i][5 - i] * pho.phep[I - i][5 - i] + P2); |
| 978 | if (IPRINT == 1) fprintf(PHLUN, "CORRECTING ENERGY OF %6i: %14.9f =>% 14.9f\n", I, pho.phep[I - i][4 - i], EN); |
| 979 | |
| 980 | pho.phep[I - i][4 - i] = EN; |
| 981 | E = E + pho.phep[I - i][4 - i]; |
| 982 | } |
| 983 | |
| 984 | |
| 985 | } |
| 986 | } |
| 987 | |
| 988 | // ----- |
| 989 | |
| 990 | if (IPRINT == 1) { |
| 991 | fprintf(PHLUN, "CORRECTING MOTHER"); |
| 992 | fprintf(PHLUN, "PX:%14.9f =>%14.9f", pho.phep[1 - i][1 - i], PX - pho.phep[2 - i][1 - i]); |
| 993 | fprintf(PHLUN, "PY:%14.9f =>%14.9f", pho.phep[1 - i][2 - i], PY - pho.phep[2 - i][2 - i]); |
| 994 | fprintf(PHLUN, "PZ:%14.9f =>%14.9f", pho.phep[1 - i][3 - i], PZ - pho.phep[2 - i][3 - i]); |
| 995 | fprintf(PHLUN, " E:%14.9f =>%14.9f", pho.phep[1 - i][4 - i], E - pho.phep[2 - i][4 - i]); |
| 996 | } |
| 997 | |
| 998 | pho.phep[1 - i][1 - i] = PX - pho.phep[2 - i][1 - i]; |
| 999 | pho.phep[1 - i][2 - i] = PY - pho.phep[2 - i][2 - i]; |
| 1000 | pho.phep[1 - i][3 - i] = PZ - pho.phep[2 - i][3 - i]; |
| 1001 | pho.phep[1 - i][4 - i] = E - pho.phep[2 - i][4 - i]; |
| 1002 | |
| 1003 | |
| 1004 | P2 = pho.phep[1 - i][1 - i] * pho.phep[1 - i][1 - i] + pho.phep[1 - i][2 - i] * pho.phep[1 - i][2 - i] + pho.phep[1 - i][3 - i] * |
| 1005 | pho.phep[1 - i][3 - i]; |
| 1006 | if (pho.phep[1 - i][4 - i]*pho.phep[1 - i][4 - i] > P2) { |
| 1007 | M = sqrt(pho.phep[1 - i][4 - i] * pho.phep[1 - i][4 - i] - P2); |
| 1008 | if (IPRINT == 1)fprintf(PHLUN, " M: %14.9f => %14.9f\n", pho.phep[1 - i][5 - i], M); |
| 1009 | pho.phep[1 - i][5 - i] = M; |
| 1010 | } |
| 1011 | |
| 1012 | PHLUPA(25); |
| 1013 | |
| 1014 | } |
| 1015 | |
| 1016 | |
| 1017 | |
| 1018 | |
| 1019 | |
| 1020 | |
| 1021 | //---------------------------------------------------------------------- |
| 1022 | // |
| 1023 | // PHOTOS: PHOton radiation in decays DOing of KINematics |
| 1024 | // |
| 1025 | // Purpose: Starting from the charged particle energy/momentum, |
| 1026 | // PNEUTR, photon energy fraction and photon angle with |
| 1027 | // respect to the axis formed by charged particle energy/ |
| 1028 | // momentum vector and PNEUTR, scale the energy/momentum, |
| 1029 | // keeping the original direction of the neutral system in |
| 1030 | // the lab. frame untouched. |
| 1031 | // |
| 1032 | // Input Parameters: IP: Pointer to decaying particle in |
| 1033 | // /PHOEVT/ and the common itself |
| 1034 | // NCHARB: pointer to the charged radiating |
| 1035 | // daughter in /PHOEVT/. |
| 1036 | // NEUDAU: pointer to the first neutral daughter |
| 1037 | // Output Parameters: Common /PHOEVT/, with photon added. |
| 1038 | // |
| 1039 | // Author(s): Z. Was, B. van Eijk Created at: 26/11/89 |
| 1040 | // Last Update: 27/05/93 |
| 1041 | // |
| 1042 | //---------------------------------------------------------------------- |
| 1043 | |
| 1044 | void PHODO(int IP, int NCHARB, int NEUDAU) |
| 1045 | { |
| 1046 | static int i = 1; |
| 1047 | double QNEW, QOLD, EPHOTO, PMAVIR; |
| 1048 | double GNEUT, DATA; |
| 1049 | double CCOSTH, SSINTH, PVEC[4], PARNE; |
| 1050 | double TH3, FI3, TH4, FI4, FI5, ANGLE; |
| 1051 | int I, J, FIRST, LAST; |
| 1052 | double& COSTHG = phophs.costhg; |
| 1053 | double& SINTHG = phophs.sinthg; |
| 1054 | double& XPHOTO = phophs.xphoto; |
| 1055 | double* PNEUTR = phomom.pneutr; |
| 1056 | double& MNESQR = phomom.mnesqr; |
| 1057 | |
| 1058 | //-- |
| 1059 | EPHOTO = XPHOTO * pho.phep[IP - i][5 - i] / 2.0; |
| 1060 | PMAVIR = sqrt(pho.phep[IP - i][5 - i] * (pho.phep[IP - i][5 - i] - 2.0 * EPHOTO)); |
| 1061 | //-- |
| 1062 | //-- Reconstruct kinematics of charged particle and neutral system |
| 1063 | phorest.fi1 = PHOAN1(PNEUTR[1 - i], PNEUTR[2 - i]); |
| 1064 | //-- |
| 1065 | //-- Choose axis along z of PNEUTR, calculate angle between x and y |
| 1066 | //-- components and z and x-y plane and perform Lorentz transform... |
| 1067 | phorest.th1 = PHOAN2(PNEUTR[3 - i], sqrt(PNEUTR[1 - i] * PNEUTR[1 - i] + PNEUTR[2 - i] * PNEUTR[2 - i])); |
| 1068 | PHORO3(-phorest.fi1, PNEUTR); |
| 1069 | PHORO2(-phorest.th1, PNEUTR); |
| 1070 | //-- |
| 1071 | //-- Take away photon energy from charged particle and PNEUTR ! Thus |
| 1072 | //-- the onshell charged particle decays into virtual charged particle |
| 1073 | //-- and photon. The virtual charged particle mass becomes: |
| 1074 | //-- SQRT(pho.phep[5,IP)*(pho.phep[5,IP)-2*EPHOTO)). Construct new PNEUTR mo- |
| 1075 | //-- mentum in the rest frame of the parent: |
| 1076 | //-- 1) Scaling parameters... |
| 1077 | QNEW = PHOTRI(PMAVIR, PNEUTR[5 - i], pho.phep[NCHARB - i][5 - i]); |
| 1078 | QOLD = PNEUTR[3 - i]; |
| 1079 | GNEUT = (QNEW * QNEW + QOLD * QOLD + MNESQR) / (QNEW * QOLD + sqrt((QNEW * QNEW + MNESQR) * (QOLD * QOLD + MNESQR))); |
| 1080 | if (GNEUT < 1.0) { |
| 1081 | DATA = 0.0; |
| 1082 | PHOERR(4, "PHOKIN", DATA); |
| 1083 | } |
| 1084 | PARNE = GNEUT - sqrt(max(GNEUT * GNEUT - 1.0, 0.0)); |
| 1085 | //-- |
| 1086 | //-- 2) ...reductive boost... |
| 1087 | PHOBO3(PARNE, PNEUTR); |
| 1088 | //-- |
| 1089 | //-- ...calculate photon energy in the reduced system... |
| 1090 | pho.nhep = pho.nhep + 1; |
| 1091 | pho.isthep[pho.nhep - i] = 1; |
| 1092 | pho.idhep[pho.nhep - i] = 22; |
| 1093 | //-- Photon mother and daughter pointers ! |
| 1094 | pho.jmohep[pho.nhep - i][1 - i] = IP; |
| 1095 | pho.jmohep[pho.nhep - i][2 - i] = 0; |
| 1096 | pho.jdahep[pho.nhep - i][1 - i] = 0; |
| 1097 | pho.jdahep[pho.nhep - i][2 - i] = -1; |
| 1098 | pho.phep[pho.nhep - i][4 - i] = EPHOTO * pho.phep[IP - i][5 - i] / PMAVIR; |
| 1099 | //-- |
| 1100 | //-- ...and photon momenta |
| 1101 | CCOSTH = -COSTHG; |
| 1102 | SSINTH = SINTHG; |
| 1103 | TH3 = PHOAN2(CCOSTH, SSINTH); |
| 1104 | FI3 = TWOPI * Photos::randomDouble(); |
| 1105 | pho.phep[pho.nhep - i][1 - i] = pho.phep[pho.nhep - i][4 - i] * SINTHG * cos(FI3); |
| 1106 | pho.phep[pho.nhep - i][2 - i] = pho.phep[pho.nhep - i][4 - i] * SINTHG * sin(FI3); |
| 1107 | //-- |
| 1108 | //-- Minus sign because axis opposite direction of charged particle ! |
| 1109 | pho.phep[pho.nhep - i][3 - i] = -pho.phep[pho.nhep - i][4 - i] * COSTHG; |
| 1110 | pho.phep[pho.nhep - i][5 - i] = 0.0; |
| 1111 | //-- |
| 1112 | //-- Rotate in order to get photon along z-axis |
| 1113 | PHORO3(-FI3, PNEUTR); |
| 1114 | PHORO3(-FI3, pho.phep[pho.nhep - i]); |
| 1115 | PHORO2(-TH3, PNEUTR); |
| 1116 | PHORO2(-TH3, pho.phep[pho.nhep - i]); |
| 1117 | ANGLE = EPHOTO / pho.phep[pho.nhep - i][4 - i]; |
| 1118 | //-- |
| 1119 | //-- Boost to the rest frame of decaying particle |
| 1120 | PHOBO3(ANGLE, PNEUTR); |
| 1121 | PHOBO3(ANGLE, pho.phep[pho.nhep - i]); |
| 1122 | //-- |
| 1123 | //-- Back in the parent rest frame but PNEUTR not yet oriented ! |
| 1124 | FI4 = PHOAN1(PNEUTR[1 - i], PNEUTR[2 - i]); |
| 1125 | TH4 = PHOAN2(PNEUTR[3 - i], sqrt(PNEUTR[1 - i] * PNEUTR[1 - i] + PNEUTR[2 - i] * PNEUTR[2 - i])); |
| 1126 | PHORO3(FI4, PNEUTR); |
| 1127 | PHORO3(FI4, pho.phep[pho.nhep - i]); |
| 1128 | //-- |
| 1129 | for (I = 2; I <= 4; I++) PVEC[I - i] = 0.0; |
| 1130 | PVEC[1 - i] = 1.0; |
| 1131 | |
| 1132 | PHORO3(-FI3, PVEC); |
| 1133 | PHORO2(-TH3, PVEC); |
| 1134 | PHOBO3(ANGLE, PVEC); |
| 1135 | PHORO3(FI4, PVEC); |
| 1136 | PHORO2(-TH4, PNEUTR); |
| 1137 | PHORO2(-TH4, pho.phep[pho.nhep - i]); |
| 1138 | PHORO2(-TH4, PVEC); |
| 1139 | FI5 = PHOAN1(PVEC[1 - i], PVEC[2 - i]); |
| 1140 | //-- |
| 1141 | //-- Charged particle restores original direction |
| 1142 | PHORO3(-FI5, PNEUTR); |
| 1143 | PHORO3(-FI5, pho.phep[pho.nhep - i]); |
| 1144 | PHORO2(phorest.th1, PNEUTR); |
| 1145 | PHORO2(phorest.th1, pho.phep[pho.nhep - i]); |
| 1146 | PHORO3(phorest.fi1, PNEUTR); |
| 1147 | PHORO3(phorest.fi1, pho.phep[pho.nhep - i]); |
| 1148 | //-- See whether neutral system has multiplicity larger than 1... |
| 1149 | |
| 1150 | if ((pho.jdahep[IP - i][2 - i] - pho.jdahep[IP - i][1 - i]) > 1) { |
| 1151 | //-- Find pointers to components of 'neutral' system |
| 1152 | //-- |
| 1153 | FIRST = NEUDAU; |
| 1154 | LAST = pho.jdahep[IP - i][2 - i]; |
| 1155 | for (I = FIRST; I <= LAST; I++) { |
| 1156 | if (I != NCHARB && (pho.jmohep[I - i][1 - i] == IP)) { |
| 1157 | //-- |
| 1158 | //-- Reconstruct kinematics... |
| 1159 | PHORO3(-phorest.fi1, pho.phep[I - i]); |
| 1160 | PHORO2(-phorest.th1, pho.phep[I - i]); |
| 1161 | //-- |
| 1162 | //-- ...reductive boost |
| 1163 | PHOBO3(PARNE, pho.phep[I - i]); |
| 1164 | //-- |
| 1165 | //-- Rotate in order to get photon along z-axis |
| 1166 | PHORO3(-FI3, pho.phep[I - i]); |
| 1167 | PHORO2(-TH3, pho.phep[I - i]); |
| 1168 | //-- |
| 1169 | //-- Boost to the rest frame of decaying particle |
| 1170 | PHOBO3(ANGLE, pho.phep[I - i]); |
| 1171 | //-- |
| 1172 | //-- Back in the parent rest-frame but PNEUTR not yet oriented. |
| 1173 | PHORO3(FI4, pho.phep[I - i]); |
| 1174 | PHORO2(-TH4, pho.phep[I - i]); |
| 1175 | //-- |
| 1176 | //-- Charged particle restores original direction |
| 1177 | PHORO3(-FI5, pho.phep[I - i]); |
| 1178 | PHORO2(phorest.th1, pho.phep[I - i]); |
| 1179 | PHORO3(phorest.fi1, pho.phep[I - i]); |
| 1180 | } |
| 1181 | } |
| 1182 | } else { |
| 1183 | //-- |
| 1184 | // ...only one 'neutral' particle in addition to photon! |
| 1185 | for (J = 1; J <= 4; J++) pho.phep[NEUDAU - i][J - i] = PNEUTR[J - i]; |
| 1186 | } |
| 1187 | //-- |
| 1188 | //-- All 'neutrals' treated, fill /PHOEVT/ for charged particle... |
| 1189 | for (J = 1; J <= 3; J++) pho.phep[NCHARB - i][J - i] = -(pho.phep[pho.nhep - i][J - i] + PNEUTR[J - i]); |
| 1190 | pho.phep[NCHARB - i][4 - i] = pho.phep[IP - i][5 - i] - (pho.phep[pho.nhep - i][4 - i] + PNEUTR[4 - i]); |
| 1191 | //-- |
| 1192 | } |
| 1193 | |
| 1194 | |
| 1195 | //---------------------------------------------------------------------- |
| 1196 | // |
| 1197 | // PHOTOS: PHOtos Boson W correction weight |
| 1198 | // |
| 1199 | // Purpose: calculates correction weight due to amplitudes of |
| 1200 | // emission from W boson. |
| 1201 | // |
| 1202 | // |
| 1203 | // |
| 1204 | // |
| 1205 | // |
| 1206 | // Input Parameters: Common /PHOEVT/, with photon added. |
| 1207 | // wt to be corrected |
| 1208 | // |
| 1209 | // |
| 1210 | // |
| 1211 | // Output Parameters: wt |
| 1212 | // |
| 1213 | // Author(s): G. Nanava, Z. Was Created at: 13/03/03 |
| 1214 | // Last Update: 08/07/13 |
| 1215 | // |
| 1216 | //---------------------------------------------------------------------- |
| 1217 | |
| 1218 | void PHOBW(double* WT) |
| 1219 | { |
| 1220 | static int i = 1; |
| 1221 | int I; |
| 1222 | double EMU, MCHREN, BETA, COSTHG, MPASQR, XPH; |
| 1223 | //-- |
| 1224 | if (abs(pho.idhep[1 - i]) == 24 && |
| 1225 | abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i]) >= 11 && |
| 1226 | abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i]) <= 16 && |
| 1227 | abs(pho.idhep[pho.jdahep[1 - i][1 - i] + 1 - i]) >= 11 && |
| 1228 | abs(pho.idhep[pho.jdahep[1 - i][1 - i] + 1 - i]) <= 16) { |
| 1229 | |
| 1230 | if ( |
| 1231 | abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i]) == 11 || |
| 1232 | abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i]) == 13 || |
| 1233 | abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i]) == 15) { |
| 1234 | I = pho.jdahep[1 - i][1 - i]; |
| 1235 | } else { |
| 1236 | I = pho.jdahep[1 - i][1 - i] + 1; |
| 1237 | } |
| 1238 | |
| 1239 | EMU = pho.phep[I - i][4 - i]; |
| 1240 | MCHREN = fabs(pow(pho.phep[I - i][4 - i], 2) - pow(pho.phep[I - i][3 - i], 2) |
| 1241 | - pow(pho.phep[I - i][2 - i], 2) - pow(pho.phep[I - i][1 - i], 2)); |
| 1242 | BETA = sqrt(1.0 - MCHREN / pho.phep[I - i][4 - i] / pho.phep[I - i][4 - i]); |
| 1243 | COSTHG = (pho.phep[I - i][3 - i] * pho.phep[pho.nhep - i][3 - i] + pho.phep[I - i][2 - i] * pho.phep[pho.nhep - i][2 - i] |
| 1244 | + pho.phep[I - i][1 - i] * pho.phep[pho.nhep - i][1 - i]) / |
| 1245 | sqrt(pho.phep[I - i][3 - i] * pho.phep[I - i][3 - i] + pho.phep[I - i][2 - i] * pho.phep[I - i][2 - i] + pho.phep[I - i][1 - i] * |
| 1246 | pho.phep[I - i][1 - i]) / |
| 1247 | sqrt(pho.phep[pho.nhep - i][3 - i] * pho.phep[pho.nhep - i][3 - i] + pho.phep[pho.nhep - i][2 - i] * pho.phep[pho.nhep - i][2 - i] + |
| 1248 | pho.phep[pho.nhep - i][1 - i] * pho.phep[pho.nhep - i][1 - i]); |
| 1249 | MPASQR = pho.phep[1 - i][4 - i] * pho.phep[1 - i][4 - i]; |
| 1250 | XPH = pho.phep[pho.nhep - i][4 - i]; |
| 1251 | *WT = (*WT) * (1 - 8 * EMU * XPH * (1 - COSTHG * BETA) * |
| 1252 | (MCHREN + 2 * XPH * sqrt(MPASQR)) / |
| 1253 | (MPASQR * MPASQR) / (1 - MCHREN / MPASQR) / (4 - MCHREN / MPASQR)); |
| 1254 | } |
| 1255 | // write(*,*) pho.idhep[1),pho.idhep[pho.jdahep[1,1)),pho.idhep[pho.jdahep[1,1)+1) |
| 1256 | // write(*,*) emu,xph,costhg,beta,mpasqr,mchren |
| 1257 | |
| 1258 | } |
| 1259 | |
| 1260 | |
| 1261 | |
| 1262 | //---------------------------------------------------------------------- |
| 1263 | // |
| 1264 | // PHOTOS: PHOton radiation in decays control FACtor |
| 1265 | // |
| 1266 | // Purpose: This is the control function for the photon spectrum and |
| 1267 | // final weighting. It is called from PHOENE for genera- |
| 1268 | // ting the raw photon energy spectrum (MODE=0) and in PHO- |
| 1269 | // COR to scale the final weight (MODE=1). The factor con- |
| 1270 | // sists of 3 terms. Addition of the factor FF which mul- |
| 1271 | // tiplies PHOFAC for MODE=0 and divides PHOFAC for MODE=1, |
| 1272 | // does not affect the results for the MC generation. An |
| 1273 | // appropriate choice for FF can speed up the calculation. |
| 1274 | // Note that a too small value of FF may cause weight over- |
| 1275 | // flow in PHOCOR and will generate a warning, halting the |
| 1276 | // execution. PRX should be included for repeated calls |
| 1277 | // for the same event, allowing more particles to radiate |
| 1278 | // photons. At the first call IREP=0, for more than 1 |
| 1279 | // charged decay products, IREP >= 1. Thus, PRSOFT (no |
| 1280 | // photon radiation probability in the previous calls) |
| 1281 | // appropriately scales the strength of the bremsstrahlung. |
| 1282 | // |
| 1283 | // Input Parameters: MODE, PROBH, XF |
| 1284 | // |
| 1285 | // Output Parameter: Function value |
| 1286 | // |
| 1287 | // Author(s): S. Jadach, Z. Was Created at: 01/01/89 |
| 1288 | // B. van Eijk, P.Golonka Last Update: 09/07/13 |
| 1289 | // |
| 1290 | //---------------------------------------------------------------------- |
| 1291 | |
| 1292 | double PHOFAC(int MODE) |
| 1293 | { |
| 1294 | static double FF = 0.0, PRX = 0.0; |
| 1295 | |
| 1296 | if (phokey.iexp) return 1.0; // In case of exponentiation this routine is useles |
| 1297 | |
| 1298 | if (MODE == -1) { |
| 1299 | PRX = 1.0; |
| 1300 | FF = 1.0; |
| 1301 | phopro.probh = 0.0; |
| 1302 | } else if (MODE == 0) { |
| 1303 | if (phopro.irep == 0) PRX = 1.0; |
| 1304 | PRX = PRX / (1.0 - phopro.probh); |
| 1305 | FF = 1.0; |
| 1306 | //-- |
| 1307 | //-- Following options are not considered for the time being... |
| 1308 | //-- (1) Good choice, but does not save very much time: |
| 1309 | //-- FF=(1.0-sqrt(phopro.xf)/2.0)/(1.0+sqrt(phopro.xf)/2.0) |
| 1310 | //-- (2) Taken from the blue, but works without weight overflows... |
| 1311 | //-- FF=(1.0-phopro.xf/(1-pow((1-sqrt(phopro.xf)),2)))*(1+(1-sqrt(phopro.xf))/sqrt(1-phopro.xf))/2.0 |
| 1312 | return FF * PRX; |
| 1313 | } else { |
| 1314 | return 1.0 / FF; |
| 1315 | } |
| 1316 | |
| 1317 | return NAN(__builtin_nanf ("")); |
| 1318 | } |
| 1319 | |
| 1320 | |
| 1321 | |
| 1322 | // ###### |
| 1323 | // replace with, |
| 1324 | // ###### |
| 1325 | |
| 1326 | //---------------------------------------------------------------------- |
| 1327 | // |
| 1328 | // PHOTOS: PHOton radiation in decays CORrection weight from |
| 1329 | // matrix elements This version for spin 1/2 is verified for |
| 1330 | // W decay only |
| 1331 | // Purpose: Calculate photon angle. The reshaping functions will |
| 1332 | // have to depend on the spin S of the charged particle. |
| 1333 | // We define: ME = 2 * S + 1 ! |
| 1334 | // THIS IS POSSIBLY ALWAYS BETTER THAN PHOCOR() |
| 1335 | // |
| 1336 | // Input Parameters: MPASQR: Parent mass squared, |
| 1337 | // MCHREN: Renormalised mass of charged system, |
| 1338 | // ME: 2 * spin + 1 determines matrix element |
| 1339 | // |
| 1340 | // Output Parameter: Function value. |
| 1341 | // |
| 1342 | // Author(s): Z. Was, B. van Eijk, G. Nanava Created at: 26/11/89 |
| 1343 | // Last Update: 01/11/12 |
| 1344 | // |
| 1345 | //---------------------------------------------------------------------- |
| 1346 | |
| 1347 | double PHOCORN(double MPASQR, double MCHREN, int ME) |
| 1348 | { |
| 1349 | double wt1, wt2, wt3; |
| 1350 | double beta0, beta1, XX, YY, DATA; |
| 1351 | double S1, PHOCOR; |
| 1352 | double& COSTHG = phophs.costhg; |
| 1353 | double& XPHMAX = phophs.xphmax; |
| 1354 | double& XPHOTO = phophs.xphoto; |
| 1355 | double& MCHSQR = phomom.mchsqr; |
| 1356 | double& MNESQR = phomom.mnesqr; |
| 1357 | |
| 1358 | |
| 1359 | |
| 1360 | //-- |
| 1361 | //-- Shaping (modified by ZW)... |
| 1362 | XX = 4.0 * MCHSQR / MPASQR * (1.0 - XPHOTO) / pow(1.0 - XPHOTO + (MCHSQR - MNESQR) / MPASQR, 2); |
| 1363 | if (ME == 1) { |
| 1364 | S1 = MPASQR * (1.0 - XPHOTO); |
| 1365 | beta0 = 2 * PHOTRI(1.0, sqrt(MCHSQR / MPASQR), sqrt(MNESQR / MPASQR)); |
| 1366 | beta1 = 2 * PHOTRI(1.0, sqrt(MCHSQR / S1), sqrt(MNESQR / S1)); |
| 1367 | wt1 = (1.0 - COSTHG * sqrt(1.0 - MCHREN)) |
| 1368 | / ((1.0 + pow(1.0 - XPHOTO / XPHMAX, 2)) / 2.0) * XPHOTO; // de-presampler |
| 1369 | |
| 1370 | wt2 = beta1 / beta0 * XPHOTO; //phase space jacobians |
| 1371 | wt3 = beta1 * beta1 * (1.0 - COSTHG * COSTHG) * (1.0 - XPHOTO) / XPHOTO / XPHOTO |
| 1372 | / pow(1.0 + MCHSQR / S1 - MNESQR / S1 - beta1 * COSTHG, 2) / 2.0; // matrix element |
| 1373 | } else if (ME == 2) { |
| 1374 | S1 = MPASQR * (1.0 - XPHOTO); |
| 1375 | beta0 = 2 * PHOTRI(1.0, sqrt(MCHSQR / MPASQR), sqrt(MNESQR / MPASQR)); |
| 1376 | beta1 = 2 * PHOTRI(1.0, sqrt(MCHSQR / S1), sqrt(MNESQR / S1)); |
| 1377 | wt1 = (1.0 - COSTHG * sqrt(1.0 - MCHREN)) |
| 1378 | / ((1.0 + pow(1.0 - XPHOTO / XPHMAX, 2)) / 2.0) * XPHOTO; // de-presampler |
| 1379 | |
| 1380 | wt2 = beta1 / beta0 * XPHOTO; // phase space jacobians |
| 1381 | |
| 1382 | wt3 = beta1 * beta1 * (1.0 - COSTHG * COSTHG) * (1.0 - XPHOTO) / XPHOTO / XPHOTO // matrix element |
| 1383 | / pow(1.0 + MCHSQR / S1 - MNESQR / S1 - beta1 * COSTHG, 2) / 2.0 ; |
| 1384 | wt3 = wt3 * (1 - XPHOTO / XPHMAX + 0.5 * pow(XPHOTO / XPHMAX, 2)) / (1 - XPHOTO / XPHMAX); |
| 1385 | // print*,"wt3=",wt3 |
| 1386 | phocorwt.phocorwt3 = wt3; |
| 1387 | phocorwt.phocorwt2 = wt2; |
| 1388 | phocorwt.phocorwt1 = wt1; |
| 1389 | |
| 1390 | // YY=0.5D0*(1.D0-XPHOTO/XPHMAX+1.D0/(1.D0-XPHOTO/XPHMAX)) |
| 1391 | // phwt.beta=SQRT(1.D0-XX) |
| 1392 | // wt1=(1.D0-COSTHG*SQRT(1.D0-MCHREN))/(1.D0-COSTHG*phwt.beta) |
| 1393 | // wt2=(1.D0-XX/YY/(1.D0-phwt.beta**2*COSTHG**2))*(1.D0+COSTHG*phwt.beta)/2.D0 |
| 1394 | // wt3=1.D0 |
| 1395 | } else if ((ME == 3) || (ME == 4) || (ME == 5)) { |
| 1396 | YY = 1.0; |
| 1397 | phwt.beta = sqrt(1.0 - XX); |
| 1398 | wt1 = (1.0 - COSTHG * sqrt(1.0 - MCHREN)) / (1.0 - COSTHG * phwt.beta); |
| 1399 | wt2 = (1.0 - XX / YY / (1.0 - phwt.beta * phwt.beta * COSTHG * COSTHG)) * (1.0 + COSTHG * phwt.beta) / 2.0; |
| 1400 | wt3 = (1.0 + pow(1.0 - XPHOTO / XPHMAX, 2) - pow(XPHOTO / XPHMAX, 3)) / |
| 1401 | (1.0 + pow(1.0 - XPHOTO / XPHMAX, 2)); |
| 1402 | } else { |
| 1403 | DATA = (ME - 1.0) / 2.0; |
| 1404 | PHOERR(6, "PHOCORN", DATA); |
| 1405 | YY = 1.0; |
| 1406 | phwt.beta = sqrt(1.0 - XX); |
| 1407 | wt1 = (1.0 - COSTHG * sqrt(1.0 - MCHREN)) / (1.0 - COSTHG * phwt.beta); |
| 1408 | wt2 = (1.0 - XX / YY / (1.0 - phwt.beta * phwt.beta * COSTHG * COSTHG)) * (1.0 + COSTHG * phwt.beta) / 2.0; |
| 1409 | wt3 = 1.0; |
| 1410 | } |
| 1411 | wt2 = wt2 * PHOFAC(1); |
| 1412 | PHOCOR = wt1 * wt2 * wt3; |
| 1413 | |
| 1414 | phopro.corwt = PHOCOR; |
| 1415 | if (PHOCOR > 1.0) { |
| 1416 | DATA = PHOCOR; |
| 1417 | PHOERR(3, "PHOCOR", DATA); |
| 1418 | } |
| 1419 | return PHOCOR; |
| 1420 | } |
| 1421 | |
| 1422 | |
| 1423 | |
| 1424 | |
| 1425 | |
| 1426 | //---------------------------------------------------------------------- |
| 1427 | // |
| 1428 | // PHOTOS: PHOton radiation in decays CORrection weight from |
| 1429 | // matrix elements |
| 1430 | // |
| 1431 | // Purpose: Calculate photon angle. The reshaping functions will |
| 1432 | // have to depend on the spin S of the charged particle. |
| 1433 | // We define: ME = 2 * S + 1 ! |
| 1434 | // |
| 1435 | // Input Parameters: MPASQR: Parent mass squared, |
| 1436 | // MCHREN: Renormalised mass of charged system, |
| 1437 | // ME: 2 * spin + 1 determines matrix element |
| 1438 | // |
| 1439 | // Output Parameter: Function value. |
| 1440 | // |
| 1441 | // Author(s): Z. Was, B. van Eijk Created at: 26/11/89 |
| 1442 | // Last Update: 21/03/93 |
| 1443 | // |
| 1444 | //---------------------------------------------------------------------- |
| 1445 | |
| 1446 | double PHOCOR(double MPASQR, double MCHREN, int ME) |
| 1447 | { |
| 1448 | double XX, YY, DATA; |
| 1449 | double PHOC; |
| 1450 | int IscaNLO; |
| 1451 | double& COSTHG = phophs.costhg; |
| 1452 | double& XPHMAX = phophs.xphmax; |
| 1453 | double& XPHOTO = phophs.xphoto; |
| 1454 | double& MCHSQR = phomom.mchsqr; |
| 1455 | double& MNESQR = phomom.mnesqr; |
| 1456 | |
| 1457 | |
| 1458 | //-- |
| 1459 | //-- Shaping (modified by ZW)... |
| 1460 | XX = 4.0 * MCHSQR / MPASQR * (1.0 - XPHOTO) / pow((1.0 - XPHOTO + (MCHSQR - MNESQR) / MPASQR), 2); |
| 1461 | if (ME == 1) { |
| 1462 | YY = 1.0; |
| 1463 | phwt.wt3 = (1.0 - XPHOTO / XPHMAX) / ((1.0 + pow((1.0 - XPHOTO / XPHMAX), 2)) / 2.0); |
| 1464 | } else if (ME == 2) { |
| 1465 | YY = 0.5 * (1.0 - XPHOTO / XPHMAX + 1.0 / (1.0 - XPHOTO / XPHMAX)); |
| 1466 | phwt.wt3 = 1.0; |
| 1467 | } else if ((ME == 3) || (ME == 4) || (ME == 5)) { |
| 1468 | YY = 1.0; |
| 1469 | phwt.wt3 = (1.0 + pow(1.0 - XPHOTO / XPHMAX, 2) - pow(XPHOTO / XPHMAX, 3)) / |
| 1470 | (1.0 + pow(1.0 - XPHOTO / XPHMAX, 2)); |
| 1471 | } else { |
| 1472 | DATA = (ME - 1.0) / 2.0; |
| 1473 | PHOERR(6, "PHOCOR", DATA); |
| 1474 | YY = 1.0; |
| 1475 | phwt.wt3 = 1.0; |
| 1476 | } |
| 1477 | |
| 1478 | |
| 1479 | phwt.beta = sqrt(1.0 - XX); |
| 1480 | phwt.wt1 = (1.0 - COSTHG * sqrt(1.0 - MCHREN)) / (1.0 - COSTHG * phwt.beta); |
| 1481 | phwt.wt2 = (1.0 - XX / YY / (1.0 - phwt.beta * phwt.beta * COSTHG * COSTHG)) * (1.0 + COSTHG * phwt.beta) / 2.0; |
| 1482 | |
| 1483 | |
| 1484 | IscaNLO = Photos::meCorrectionWtForScalar; |
| 1485 | if (ME == 1 && IscaNLO == 1) { // this switch NLO in scalar decays. |
| 1486 | // overrules default calculation. |
| 1487 | // Need tests including basic ones |
| 1488 | PHOC = PHOCORN(MPASQR, MCHREN, ME); |
| 1489 | phwt.wt1 = 1.0; |
| 1490 | phwt.wt2 = 1.0; |
| 1491 | phwt.wt3 = PHOC; |
| 1492 | } else { |
| 1493 | phwt.wt2 = phwt.wt2 * PHOFAC(1); |
| 1494 | } |
| 1495 | PHOC = phwt.wt1 * phwt.wt2 * phwt.wt3; |
| 1496 | |
| 1497 | phopro.corwt = PHOC; |
| 1498 | if (PHOC > 1.0) { |
| 1499 | DATA = PHOC; |
| 1500 | PHOERR(3, "PHOCOR", DATA); |
| 1501 | } |
| 1502 | return PHOC; |
| 1503 | } |
| 1504 | |
| 1505 | |
| 1506 | //---------------------------------------------------------------------- |
| 1507 | // |
| 1508 | // PHOTWO: PHOtos but TWO mothers allowed |
| 1509 | // |
| 1510 | // Purpose: Combines two mothers into one in /PHOEVT/ |
| 1511 | // necessary eg in case of g g (q qbar) --> t tbar |
| 1512 | // |
| 1513 | // Input Parameters: Common /PHOEVT/ (/PHOCMS/) |
| 1514 | // |
| 1515 | // Output Parameters: Common /PHOEVT/, (stored mothers) |
| 1516 | // |
| 1517 | // Author(s): Z. Was Created at: 5/08/93 |
| 1518 | // Last Update:10/08/93 |
| 1519 | // |
| 1520 | //---------------------------------------------------------------------- |
| 1521 | |
| 1522 | void PHOTWO(int MODE) |
| 1523 | { |
| 1524 | |
| 1525 | int I; |
| 1526 | static int i = 1; |
| 1527 | double MPASQR; |
| 1528 | bool IFRAD; |
| 1529 | // logical IFRAD is used to tag cases when two mothers may be |
| 1530 | // merged to the sole one. |
| 1531 | // So far used in case: |
| 1532 | // 1) of t tbar production |
| 1533 | // |
| 1534 | // t tbar case |
| 1535 | if (MODE == 0) { |
| 1536 | IFRAD = (pho.idhep[1 - i] == 21) && (pho.idhep[2 - i] == 21); |
| 1537 | IFRAD = IFRAD || (pho.idhep[1 - i] == -pho.idhep[2 - i] && abs(pho.idhep[1 - i]) <= 6); |
| 1538 | IFRAD = IFRAD && (abs(pho.idhep[3 - i]) == 6) && (abs(pho.idhep[4 - i]) == 6); |
| 1539 | MPASQR = pow(pho.phep[1 - i][4 - i] + pho.phep[2 - i][4 - i], 2) - pow(pho.phep[1 - i][3 - i] + pho.phep[2 - i][3 - i], 2) |
| 1540 | - pow(pho.phep[1 - i][2 - i] + pho.phep[2 - i][2 - i], 2) - pow(pho.phep[1 - i][1 - i] + pho.phep[2 - i][1 - i], 2); |
| 1541 | IFRAD = IFRAD && (MPASQR > 0.0); |
| 1542 | if (IFRAD) { |
| 1543 | //.....combining first and second mother |
| 1544 | for (I = 1; I <= 4; I++) { |
| 1545 | pho.phep[1 - i][I - i] = pho.phep[1 - i][I - i] + pho.phep[2 - i][I - i]; |
| 1546 | } |
| 1547 | pho.phep[1 - i][5 - i] = sqrt(MPASQR); |
| 1548 | //.....removing second mother, |
| 1549 | for (I = 1; I <= 5; I++) { |
| 1550 | pho.phep[2 - i][I - i] = 0.0; |
| 1551 | } |
| 1552 | } |
| 1553 | } else { |
| 1554 | // boosting of the mothers to the reaction frame not implemented yet. |
| 1555 | // to do it in mode 0 original mothers have to be stored in new comon (?) |
| 1556 | // and in mode 1 boosted to cms. |
| 1557 | } |
| 1558 | } |
| 1559 | |
| 1560 | |
| 1561 | |
| 1562 | //---------------------------------------------------------------------- |
| 1563 | // |
| 1564 | // PHOTOS: PHOtos CDE-s |
| 1565 | // |
| 1566 | // Purpose: Keep definitions for PHOTOS QED correction Monte Carlo. |
| 1567 | // |
| 1568 | // Input Parameters: None |
| 1569 | // |
| 1570 | // Output Parameters: None |
| 1571 | // |
| 1572 | // Author(s): Z. Was, B. van Eijk Created at: 29/11/89 |
| 1573 | // Last Update: 10/08/93 |
| 1574 | // |
| 1575 | // ========================================================= |
| 1576 | // General Structure Information: = |
| 1577 | // ========================================================= |
| 1578 | //: ROUTINES: |
| 1579 | // 1) INITIALIZATION (all in C++ now) |
| 1580 | // 2) GENERAL INTERFACE: |
| 1581 | // PHOBOS |
| 1582 | // PHOIN |
| 1583 | // PHOTWO (specific interface |
| 1584 | // PHOOUT |
| 1585 | // PHOCHK |
| 1586 | // PHTYPE (specific interface |
| 1587 | // PHOMAK (specific interface |
| 1588 | // 3) QED PHOTON GENERATION: |
| 1589 | // PHINT |
| 1590 | // PHOBW |
| 1591 | // PHOPRE |
| 1592 | // PHOOMA |
| 1593 | // PHOENE |
| 1594 | // PHOCOR |
| 1595 | // PHOFAC |
| 1596 | // PHODO |
| 1597 | // 4) UTILITIES: |
| 1598 | // PHOTRI |
| 1599 | // PHOAN1 |
| 1600 | // PHOAN2 |
| 1601 | // PHOBO3 |
| 1602 | // PHORO2 |
| 1603 | // PHORO3 |
| 1604 | // PHOCHA |
| 1605 | // PHOSPI |
| 1606 | // PHOERR |
| 1607 | // PHOREP |
| 1608 | // PHLUPA |
| 1609 | // PHCORK |
| 1610 | // IPHQRK |
| 1611 | // IPHEKL |
| 1612 | // COMMONS: |
| 1613 | // NAME USED IN SECT. # OF OC// Comment |
| 1614 | // PHOQED 1) 2) 3 Flags whether emisson to be gen. |
| 1615 | // PHOLUN 1) 4) 6 Output device number |
| 1616 | // PHOCOP 1) 3) 4 photon coupling & min energy |
| 1617 | // PHPICO 1) 3) 4) 5 PI & 2*PI |
| 1618 | // PHOSTA 1) 4) 3 Status information |
| 1619 | // PHOKEY 1) 2) 3) 7 Keys for nonstandard application |
| 1620 | // PHOVER 1) 1 Version info for outside |
| 1621 | // HEPEVT 2) 2 PDG common |
| 1622 | // PH_HEPEVT2) 8 PDG common internal |
| 1623 | // PHOEVT 2) 3) 10 PDG branch |
| 1624 | // PHOIF 2) 3) 2 emission flags for PDG branch |
| 1625 | // PHOMOM 3) 5 param of char-neutr system |
| 1626 | // PHOPHS 3) 5 photon momentum parameters |
| 1627 | // PHOPRO 3) 4 var. for photon rep. (in branch) |
| 1628 | // PHOCMS 2) 3 parameters of boost to branch CMS |
| 1629 | // PHNUM 4) 1 event number from outside |
| 1630 | //---------------------------------------------------------------------- |
| 1631 | |
| 1632 | |
| 1633 | //---------------------------------------------------------------------- |
| 1634 | // |
| 1635 | // PHOIN: PHOtos INput |
| 1636 | // |
| 1637 | // Purpose: copies IP branch of the common /PH_HEPEVT/ into /PHOEVT/ |
| 1638 | // moves branch into its CMS system. |
| 1639 | // |
| 1640 | // Input Parameters: IP: pointer of particle starting branch |
| 1641 | // to be copied |
| 1642 | // BOOST: Flag whether boost to CMS was or was |
| 1643 | // . replace stri not performed. |
| 1644 | // |
| 1645 | // Output Parameters: Commons: /PHOEVT/, /PHOCMS/ |
| 1646 | // |
| 1647 | // Author(s): Z. Was Created at: 24/05/93 |
| 1648 | // Last Update: 16/11/93 |
| 1649 | // |
| 1650 | //---------------------------------------------------------------------- |
| 1651 | void PHOIN(int IP, bool* BOOST, int* NHEP0) |
| 1652 | { |
| 1653 | int FIRST, LAST, I, LL, IP2, J, NA; |
| 1654 | double PB; |
| 1655 | static int i = 1; |
| 1656 | int& nhep0 = *NHEP0; |
| 1657 | // double &BET[3]=BET; |
| 1658 | double& GAM = phocms.gam; |
| 1659 | double* BET = phocms.bet; |
| 1660 | |
| 1661 | //-- |
| 1662 | // let-s calculate size of the little common entry |
| 1663 | FIRST = hep.jdahep[IP - i][1 - i]; |
| 1664 | LAST = hep.jdahep[IP - i][2 - i]; |
| 1665 | pho.nhep = 3 + LAST - FIRST + hep.nhep - nhep0; |
| 1666 | pho.nevhep = pho.nhep; |
| 1667 | |
| 1668 | // let-s take in decaying particle |
| 1669 | pho.idhep[1 - i] = hep.idhep[IP - i]; |
| 1670 | pho.jdahep[1 - i][1 - i] = 3; |
| 1671 | pho.jdahep[1 - i][2 - i] = 3 + LAST - FIRST; |
| 1672 | for (I = 1; I <= 5; I++) pho.phep[1 - i][I - i] = hep.phep[IP - i][I - i]; |
| 1673 | for (I = 1; I <= 4; I++) pho.vhep[1 - i][I - i] = hep.vhep[IP - i][I - i]; |
| 1674 | |
| 1675 | // let-s take in eventual second mother |
| 1676 | IP2 = hep.jmohep[hep.jdahep[IP - i][1 - i] - i][2 - i]; |
| 1677 | if ((IP2 != 0) && (IP2 != IP)) { |
| 1678 | pho.idhep[2 - i] = hep.idhep[IP2 - i]; |
| 1679 | pho.jdahep[2 - i][1 - i] = 3; |
| 1680 | pho.jdahep[2 - i][2 - i] = 3 + LAST - FIRST; |
| 1681 | for (I = 1; I <= 5; I++) |
| 1682 | pho.phep[2 - i][I - i] = hep.phep[IP2 - i][I - i]; |
| 1683 | for (I = 1; I <= 4; I++) |
| 1684 | pho.vhep[2 - i][I - i] = hep.vhep[IP2 - i][I - i]; |
| 1685 | } else { |
| 1686 | pho.idhep[2 - i] = 0; |
| 1687 | for (I = 1; I <= 5; I++) pho.phep[2 - i][I - i] = 0.0; |
| 1688 | for (I = 1; I <= 4; I++) pho.vhep[2 - i][I - i] = 0.0; |
| 1689 | } |
| 1690 | |
| 1691 | // let-s take in daughters |
| 1692 | for (LL = 0; LL <= LAST - FIRST; LL++) { |
| 1693 | pho.idhep[3 + LL - i] = hep.idhep[FIRST + LL - i]; |
| 1694 | pho.jmohep[3 + LL - i][1 - i] = hep.jmohep[FIRST + LL - i][1 - i]; |
| 1695 | if (hep.jmohep[FIRST + LL - i][1 - i] == IP) pho.jmohep[3 + LL - i][1 - i] = 1; |
| 1696 | for (I = 1; I <= 5; I++) pho.phep[3 + LL - i][I - i] = hep.phep[FIRST + LL - i][I - i]; |
| 1697 | for (I = 1; I <= 4; I++) pho.vhep[3 + LL - i][I - i] = hep.vhep[FIRST + LL - i][I - i]; |
| 1698 | } |
| 1699 | if (hep.nhep > nhep0) { |
| 1700 | // let-s take in illegitimate daughters |
| 1701 | NA = 3 + LAST - FIRST; |
| 1702 | for (LL = 1; LL <= hep.nhep - nhep0; LL++) { |
| 1703 | pho.idhep[NA + LL - i] = hep.idhep[nhep0 + LL - i]; |
| 1704 | pho.jmohep[NA + LL - i][1 - i] = hep.jmohep[nhep0 + LL - i][1 - i]; |
| 1705 | if (hep.jmohep[nhep0 + LL - i][1 - i] == IP) pho.jmohep[NA + LL - i][1 - i] = 1; |
| 1706 | for (I = 1; I <= 5; I++) pho.phep[NA + LL - i][I - i] = hep.phep[nhep0 + LL - i][I - i]; |
| 1707 | for (I = 1; I <= 4; I++) pho.vhep[NA + LL - i][I - i] = hep.vhep[nhep0 + LL - i][I - i]; |
| 1708 | } |
| 1709 | //-- there is hep.nhep-nhep0 daugters more. |
| 1710 | pho.jdahep[1 - i][2 - i] = 3 + LAST - FIRST + hep.nhep - nhep0; |
| 1711 | } |
| 1712 | if (pho.idhep[pho.nhep - i] == 22) PHLUPA(100001); |
| 1713 | // if (pho.idhep[pho.nhep-i]==22) exit(-1); |
| 1714 | PHCORK(0); |
| 1715 | if (pho.idhep[pho.nhep - i] == 22) PHLUPA(100002); |
| 1716 | |
| 1717 | // special case of t tbar production process |
| 1718 | if (phokey.iftop) PHOTWO(0); |
| 1719 | *BOOST = false; |
| 1720 | |
| 1721 | //-- Check whether parent is in its rest frame... |
| 1722 | // ZBW ND 27.07.2009: |
| 1723 | // bug reported by Vladimir Savinov localized and fixed. |
| 1724 | // protection against rounding error was back-firing if soft |
| 1725 | // momentum of mother was physical. Consequence was that PHCORK was |
| 1726 | // messing up masses of final state particles in vertex of the decay. |
| 1727 | // Only configurations with previously generated photons of energy fraction |
| 1728 | // smaller than 0.0001 were affected. Effect was numerically insignificant. |
| 1729 | |
| 1730 | // IF ( (ABS(pho.phep[4,1)-pho.phep[5,1)).GT.pho.phep[5,1)*1.D-8) |
| 1731 | // $ .AND.(pho.phep[5,1).NE.0)) THEN |
| 1732 | |
| 1733 | if ((fabs(pho.phep[1 - i][1 - i] + fabs(pho.phep[1 - i][2 - i]) + fabs(pho.phep[1 - i][3 - i])) > |
| 1734 | pho.phep[1 - i][5 - i] * 1.E-8) && (pho.phep[1 - i][5 - i] != 0)) { |
| 1735 | |
| 1736 | *BOOST = true; |
| 1737 | //PHOERR(404,"PHOIN",1.0); // we need to improve this warning: program should never |
| 1738 | // enter this place |
| 1739 | // may be exit(-1); |
| 1740 | //-- |
| 1741 | //-- Boost daughter particles to rest frame of parent... |
| 1742 | //-- Resultant neutral system already calculated in rest frame ! |
| 1743 | for (J = 1; J <= 3; J++) BET[J - i] = -pho.phep[1 - i][J - i] / pho.phep[1 - i][5 - i]; |
| 1744 | GAM = pho.phep[1 - i][4 - i] / pho.phep[1 - i][5 - i]; |
| 1745 | for (I = pho.jdahep[1 - i][1 - i]; I <= pho.jdahep[1 - i][2 - i]; I++) { |
| 1746 | PB = BET[1 - i] * pho.phep[I - i][1 - i] + BET[2 - i] * pho.phep[I - i][2 - i] + BET[3 - i] * pho.phep[I - i][3 - i]; |
| 1747 | for (J = 1; J <= 3; |
| 1748 | J++) pho.phep[I - i][J - i] = pho.phep[I - i][J - i] + BET[J - i] * (pho.phep[I - i][4 - i] + PB / (GAM + 1.0)); |
| 1749 | pho.phep[I - i][4 - i] = GAM * pho.phep[I - i][4 - i] + PB; |
| 1750 | } |
| 1751 | //-- Finally boost mother as well |
| 1752 | I = 1; |
| 1753 | PB = BET[1 - i] * pho.phep[I - i][1 - i] + BET[2 - i] * pho.phep[I - i][2 - i] + BET[3 - i] * pho.phep[I - i][3 - i]; |
| 1754 | for (J = 1; J <= 3; J++) pho.phep[I - i][J - i] = pho.phep[I - i][J - i] + BET[J - i] * (pho.phep[I - i][4 - i] + PB / (GAM + 1.0)); |
| 1755 | |
| 1756 | pho.phep[I - i][4 - i] = GAM * pho.phep[I - i][4 - i] + PB; |
| 1757 | |
| 1758 | for (J = 1; J <= 3; J++) BET[J - i] = -pho.phep[1 - i][J - i] / pho.phep[1 - i][5 - i]; |
| 1759 | GAM = pho.phep[1 - i][4 - i] / pho.phep[1 - i][5 - i]; |
| 1760 | for (I = pho.jdahep[1 - i][1 - i]; I <= pho.jdahep[1 - i][2 - i]; I++) { |
| 1761 | PB = BET[1 - i] * pho.vhep[I - i][1 - i] + BET[2 - i] * pho.vhep[I - i][2 - i] + BET[3 - i] * pho.vhep[I - i][3 - i]; |
| 1762 | for (J = 1; J <= 3; J++) pho.vhep[I - i][J - i] = pho.vhep[I - i][J - i] + BET[J - i] * (pho.vhep[I - i][4 - i] + PB / (GAM + 1.0)); |
| 1763 | pho.vhep[I - i][4 - i] = GAM * pho.vhep[I - i][4 - i] + PB; |
| 1764 | } |
| 1765 | //-- Finally boost mother as well |
| 1766 | I = 1; |
| 1767 | PB = BET[1 - i] * pho.vhep[I - i][1 - i] + BET[2 - i] * pho.vhep[I - i][2 - i] + BET[3 - i] * pho.vhep[I - i][3 - i]; |
| 1768 | for (J = 1; J <= 3; J++) pho.vhep[I - i][J - i] = pho.vhep[I - i][J - i] + BET[J - i] * (pho.vhep[I - i][4 - i] + PB / (GAM + 1.0)); |
| 1769 | |
| 1770 | pho.vhep[I - i][4 - i] = GAM * pho.vhep[I - i][4 - i] + PB; |
| 1771 | } |
| 1772 | |
| 1773 | |
| 1774 | // special case of t tbar production process |
| 1775 | if (phokey.iftop) PHOTWO(1); |
| 1776 | PHLUPA(2); |
| 1777 | if (pho.idhep[pho.nhep - i] == 22) PHLUPA(10000); |
| 1778 | //if (pho.idhep[pho.nhep-1-i]==22) exit(-1); // this is probably form very old times ... |
| 1779 | return; |
| 1780 | } |
| 1781 | |
| 1782 | |
| 1783 | //---------------------------------------------------------------------- |
| 1784 | // |
| 1785 | // PHOOUT: PHOtos OUTput |
| 1786 | // |
| 1787 | // Purpose: copies back IP branch of the common /PH_HEPEVT/ from |
| 1788 | // /PHOEVT/ moves branch back from its CMS system. |
| 1789 | // |
| 1790 | // Input Parameters: IP: pointer of particle starting branch |
| 1791 | // to be given back. |
| 1792 | // BOOST: Flag whether boost to CMS was or was |
| 1793 | // . not performed. |
| 1794 | // |
| 1795 | // Output Parameters: Common /PHOEVT/, |
| 1796 | // |
| 1797 | // Author(s): Z. Was Created at: 24/05/93 |
| 1798 | // Last Update: |
| 1799 | // |
| 1800 | //---------------------------------------------------------------------- |
| 1801 | void PHOOUT(int IP, bool BOOST, int nhep0) |
| 1802 | { |
| 1803 | int LL, FIRST, LAST, I; |
| 1804 | int NN, J, K, NA; |
| 1805 | double PB; |
| 1806 | double& GAM = phocms.gam; |
| 1807 | double* BET = phocms.bet; |
| 1808 | |
| 1809 | static int i = 1; |
| 1810 | if (pho.nhep == pho.nevhep) return; |
| 1811 | //-- When parent was not in its rest-frame, boost back... |
| 1812 | PHLUPA(10); |
| 1813 | if (BOOST) { |
| 1814 | //PHOERR(404,"PHOOUT",1.0); // we need to improve this warning: program should never |
| 1815 | // enter this place |
| 1816 | |
| 1817 | double phocms_check = fabs(1 - GAM) + fabs(BET[1 - i]) + fabs(BET[2 - i]) + fabs(BET[3 - i]); |
| 1818 | if (phocms_check > 0.001) { |
| 1819 | Log::Error() << "Msg. from PHOOUT: possible problems with boosting due to the rounding errors." << endl |
| 1820 | << "Boost parameters: (" << GAM << "," |
| 1821 | << BET[1 - i] << "," << BET[2 - i] << "," << BET[3 - i] << ")" << endl |
| 1822 | << "should be equal to: (1,0,0,0) up to at least several digits." << endl; |
| 1823 | } else { |
| 1824 | Log::Warning() << "Msg. from PHOOUT: possible problems with boosting due to the rounding errors." << endl |
| 1825 | << "Boost parameters: (" << GAM << "," |
| 1826 | << BET[1 - i] << "," << BET[2 - i] << "," << BET[3 - i] << ")" << endl |
| 1827 | << "should be equal to: (1,0,0,0) up to at least several digits." << endl; |
| 1828 | } |
| 1829 | |
| 1830 | for (J = pho.jdahep[1 - i][1 - i]; J <= pho.jdahep[1 - i][2 - i]; J++) { |
| 1831 | PB = -BET[1 - i] * pho.phep[J - i][1 - i] - BET[2 - i] * pho.phep[J - i][2 - i] - BET[3 - i] * pho.phep[J - i][3 - i]; |
| 1832 | for (K = 1; K <= 3; K++) pho.phep[J - i][K - i] = pho.phep[J - i][K - i] - BET[K - i] * (pho.phep[J - i][4 - i] + PB / (GAM + 1.0)); |
| 1833 | pho.phep[J - i][4 - i] = GAM * pho.phep[J - i][4 - i] + PB; |
| 1834 | } |
| 1835 | |
| 1836 | //-- ...boost photon, or whatever else has shown up |
| 1837 | for (NN = pho.nevhep + 1; NN <= pho.nhep; NN++) { |
| 1838 | PB = -BET[1 - i] * pho.phep[NN - i][1 - i] - BET[2 - i] * pho.phep[NN - i][2 - i] - BET[3 - i] * pho.phep[NN - i][3 - i]; |
| 1839 | for (K = 1; K <= 3; |
| 1840 | K++) pho.phep[NN - i][K - i] = pho.phep[NN - i][K - i] - BET[K - i] * (pho.phep[NN - i][4 - i] + PB / (GAM + 1.0)); |
| 1841 | pho.phep[NN - i][4 - i] = GAM * pho.phep[NN][4 - i] + PB; |
| 1842 | } |
| 1843 | |
| 1844 | // same for vertex |
| 1845 | for (J = pho.jdahep[1 - i][1 - i]; J <= pho.jdahep[1 - i][2 - i]; J++) { |
| 1846 | PB = -BET[1 - i] * pho.vhep[J - i][1 - i] - BET[2 - i] * pho.vhep[J - i][2 - i] - BET[3 - i] * pho.vhep[J - i][3 - i]; |
| 1847 | for (K = 1; K <= 3; K++) pho.vhep[J - i][K - i] = pho.vhep[J - i][K - i] - BET[K - i] * (pho.vhep[J - i][4 - i] + PB / (GAM + 1.0)); |
| 1848 | pho.vhep[J - i][4 - i] = GAM * pho.vhep[J - i][4 - i] + PB; |
| 1849 | } |
| 1850 | |
| 1851 | //-- ...boost photon, or whatever else has shown up |
| 1852 | for (NN = pho.nevhep + 1; NN <= pho.nhep; NN++) { |
| 1853 | PB = -BET[1 - i] * pho.vhep[NN - i][1 - i] - BET[2 - i] * pho.vhep[NN - i][2 - i] - BET[3 - i] * pho.vhep[NN - i][3 - i]; |
| 1854 | for (K = 1; K <= 3; |
| 1855 | K++) pho.vhep[NN - i][K - i] = pho.vhep[NN - i][K - i] - BET[K - i] * (pho.vhep[NN - i][4 - i] + PB / (GAM + 1.0)); |
| 1856 | pho.vhep[NN - i][4 - i] = GAM * pho.vhep[NN][4 - i] + PB; |
| 1857 | } |
| 1858 | } |
| 1859 | PHCORK(0); // we have to use it because it clears input |
| 1860 | // for grandaughters modified in C++ |
| 1861 | FIRST = hep.jdahep[IP - i][1 - i]; |
| 1862 | LAST = hep.jdahep[IP - i][2 - i]; |
| 1863 | // let-s take in original daughters |
| 1864 | for (LL = 0; LL <= LAST - FIRST; LL++) { |
| 1865 | hep.idhep[FIRST + LL - i] = pho.idhep[3 + LL - i]; |
| 1866 | for (I = 1; I <= 5; I++) hep.phep[FIRST + LL - i][I - i] = pho.phep[3 + LL - i][I - i]; |
| 1867 | for (I = 1; I <= 4; I++) hep.vhep[FIRST + LL - i][I - i] = pho.vhep[3 + LL - i][I - i]; |
| 1868 | } |
| 1869 | |
| 1870 | // let-s take newcomers to the end of HEPEVT. |
| 1871 | NA = 3 + LAST - FIRST; |
| 1872 | for (LL = 1; LL <= pho.nhep - NA; LL++) { |
| 1873 | hep.idhep[nhep0 + LL - i] = pho.idhep[NA + LL - i]; |
| 1874 | hep.isthep[nhep0 + LL - i] = pho.isthep[NA + LL - i]; |
| 1875 | hep.jmohep[nhep0 + LL - i][1 - i] = IP; |
| 1876 | hep.jmohep[nhep0 + LL - i][2 - i] = hep.jmohep[hep.jdahep[IP - i][1 - i] - i][2 - i]; |
| 1877 | hep.jdahep[nhep0 + LL - i][1 - i] = 0; |
| 1878 | hep.jdahep[nhep0 + LL - i][2 - i] = 0; |
| 1879 | for (I = 1; I <= 5; I++) hep.phep[nhep0 + LL - i][I - i] = pho.phep[NA + LL - i][I - i]; |
| 1880 | for (I = 1; I <= 4; I++) hep.vhep[nhep0 + LL - i][I - i] = pho.vhep[NA + LL - i][I - i]; |
| 1881 | } |
| 1882 | hep.nhep = hep.nhep + pho.nhep - pho.nevhep; |
| 1883 | PHLUPA(20); |
| 1884 | return; |
| 1885 | } |
| 1886 | |
| 1887 | //---------------------------------------------------------------------- |
| 1888 | // |
| 1889 | // PHOCHK: checking branch. |
| 1890 | // |
| 1891 | // Purpose: checks whether particles in the common block /PHOEVT/ |
| 1892 | // can be served by PHOMAK. |
| 1893 | // JFIRST is the position in /PH_HEPEVT/ (!) of the first |
| 1894 | // daughter of sub-branch under action. |
| 1895 | // |
| 1896 | // |
| 1897 | // Author(s): Z. Was Created at: 22/10/92 |
| 1898 | // Last Update: 11/12/00 |
| 1899 | // |
| 1900 | //---------------------------------------------------------------------- |
| 1901 | // ******************** |
| 1902 | |
| 1903 | void PHOCHK(int JFIRST) |
| 1904 | { |
| 1905 | |
| 1906 | int IDABS, NLAST, I; |
| 1907 | bool IFRAD; |
| 1908 | int IDENT, K; |
| 1909 | static int i = 1, IPPAR = 1; |
| 1910 | |
| 1911 | NLAST = pho.nhep; |
| 1912 | // |
| 1913 | |
| 1914 | for (I = IPPAR; I <= NLAST; I++) { |
| 1915 | IDABS = abs(pho.idhep[I - i]); |
| 1916 | // possibly call on PHZODE is a dead (to be omitted) code. |
| 1917 | pho.qedrad[I - i] = pho.qedrad[I - i] && F(0, IDABS) && F(0, abs(pho.idhep[1 - i])) |
| 1918 | && (pho.idhep[2 - i] == 0); |
| 1919 | |
| 1920 | if (I > 2) pho.qedrad[I - i] = pho.qedrad[I - i] && hep.qedrad[JFIRST + I - IPPAR - 2 - i]; |
| 1921 | } |
| 1922 | |
| 1923 | //-- |
| 1924 | // now we go to special cases, where pho.qedrad[I) will be overwritten |
| 1925 | //-- |
| 1926 | IDENT = pho.nhep; |
| 1927 | if (phokey.iftop) { |
| 1928 | // special case of top pair production |
| 1929 | for (K = pho.jdahep[1 - i][2 - i]; K >= pho.jdahep[1 - i][1 - i]; K--) { |
| 1930 | if (pho.idhep[K - i] != 22) { |
| 1931 | IDENT = K; |
| 1932 | break; // from loop over K |
| 1933 | } |
| 1934 | } |
| 1935 | |
| 1936 | IFRAD = ((pho.idhep[1 - i] == 21) && (pho.idhep[2 - i] == 21)) |
| 1937 | || ((abs(pho.idhep[1 - i]) <= 6) && (pho.idhep[2 - i] == (-pho.idhep[1 - i]))); |
| 1938 | IFRAD = IFRAD |
| 1939 | && (abs(pho.idhep[3 - i]) == 6) && (pho.idhep[4 - i] == (-pho.idhep[3 - i])) |
| 1940 | && (IDENT == 4); |
| 1941 | if (IFRAD) { |
| 1942 | for (I = IPPAR; I <= NLAST; I++) { |
| 1943 | pho.qedrad[I - i] = true; |
| 1944 | if (I > 2) pho.qedrad[I - i] = pho.qedrad[I - i] && hep.qedrad[JFIRST + I - IPPAR - 2 - i]; |
| 1945 | } |
| 1946 | } |
| 1947 | } |
| 1948 | //-- |
| 1949 | //-- |
| 1950 | if (phokey.iftop) { |
| 1951 | // special case of top decay |
| 1952 | for (K = pho.jdahep[1 - i][2 - i]; K >= pho.jdahep[1 - i][1 - i]; K--) { |
| 1953 | if (pho.idhep[K - i] != 22) { |
| 1954 | IDENT = K; |
| 1955 | break; |
| 1956 | } |
| 1957 | } |
| 1958 | IFRAD = ((abs(pho.idhep[1 - i]) == 6) && (pho.idhep[2 - i] == 0)); |
| 1959 | IFRAD = IFRAD |
| 1960 | && (((abs(pho.idhep[3 - i]) == 24) && (abs(pho.idhep[4 - i]) == 5)) |
| 1961 | || ((abs(pho.idhep[3 - i]) == 5) && (abs(pho.idhep[4 - i]) == 24))) |
| 1962 | && (IDENT == 4); |
| 1963 | |
| 1964 | if (IFRAD) { |
| 1965 | for (I = IPPAR; I <= NLAST; I++) { |
| 1966 | pho.qedrad[I - i] = true; |
| 1967 | if (I > 2) pho.qedrad[I - i] = (pho.qedrad[I - i] && hep.qedrad[JFIRST + I - IPPAR - 2 - i]); |
| 1968 | } |
| 1969 | } |
| 1970 | } |
| 1971 | //-- |
| 1972 | //-- |
| 1973 | return; |
| 1974 | } |
| 1975 | |
| 1976 | |
| 1977 | |
| 1978 | //---------------------------------------------------------------------- |
| 1979 | // |
| 1980 | // PHOTOS: PHOton radiation in decays calculation of photon ENErgy |
| 1981 | // fraction |
| 1982 | // |
| 1983 | // Purpose: Subroutine returns photon energy fraction (in (parent |
| 1984 | // mass)/2 units) for the decay bremsstrahlung. |
| 1985 | // |
| 1986 | // Input Parameters: MPASQR: Mass of decaying system squared, |
| 1987 | // XPHCUT: Minimum energy fraction of photon, |
| 1988 | // XPHMAX: Maximum energy fraction of photon. |
| 1989 | // |
| 1990 | // Output Parameter: MCHREN: Renormalised mass squared, |
| 1991 | // BETA: Beta factor due to renormalisation, |
| 1992 | // XPHOTO: Photon energy fraction, |
| 1993 | // XF: Correction factor for PHOFA// |
| 1994 | // |
| 1995 | // Author(s): S. Jadach, Z. Was Created at: 01/01/89 |
| 1996 | // B. van Eijk, P.Golonka Last Update: 11/07/13 |
| 1997 | // |
| 1998 | //---------------------------------------------------------------------- |
| 1999 | |
| 2000 | void PHOENE(double MPASQR, double* pMCHREN, double* pBETA, double* pBIGLOG, int IDENT) |
| 2001 | { |
| 2002 | double DATA; |
| 2003 | double PRSOFT, PRHARD; |
| 2004 | double PRKILL, RRR; |
| 2005 | int K, IDME; |
| 2006 | double PRSUM; |
| 2007 | static int i = 1; |
| 2008 | double& MCHREN = *pMCHREN; |
| 2009 | double& BETA = *pBETA; |
| 2010 | double& BIGLOG = *pBIGLOG; |
| 2011 | int& NCHAN = phoexp.nchan; |
| 2012 | double& XPHMAX = phophs.xphmax; |
| 2013 | double& XPHOTO = phophs.xphoto; |
| 2014 | double& MCHSQR = phomom.mchsqr; |
| 2015 | double& MNESQR = phomom.mnesqr; |
| 2016 | |
| 2017 | //-- |
| 2018 | if (XPHMAX <= phocop.xphcut) { |
| 2019 | BETA = PHOFAC(-1); // to zero counter, here beta is dummy |
| 2020 | XPHOTO = 0.0; |
| 2021 | return; |
| 2022 | } |
| 2023 | //-- Probabilities for hard and soft bremstrahlung... |
| 2024 | MCHREN = 4.0 * MCHSQR / MPASQR / pow(1.0 + MCHSQR / MPASQR, 2); |
| 2025 | BETA = sqrt(1.0 - MCHREN); |
| 2026 | |
| 2027 | #ifdef VARIANTB |
| 2028 | // ----------- VARIANT B ------------------ |
| 2029 | // we replace 1D0/BETA*BIGLOG with (1.0/BETA*BIGLOG+2*phokey.fint) |
| 2030 | // for integral of new crude |
| 2031 | BIGLOG = log(MPASQR / MCHSQR * (1.0 + BETA) * (1.0 + BETA) / 4.0 * |
| 2032 | pow(1.0 + MCHSQR / MPASQR, 2)); |
| 2033 | PRHARD = phocop.alpha / PI * (1.0 / BETA * BIGLOG + 2 * phokey.fint) |
| 2034 | * (log(XPHMAX / phocop.xphcut) - .75 + phocop.xphcut / XPHMAX - .25 * phocop.xphcut * phocop.xphcut / XPHMAX / XPHMAX); |
| 2035 | PRHARD = PRHARD * PHOCHA(IDENT) * PHOCHA(IDENT) * phokey.fsec; |
| 2036 | // ----------- END OF VARIANT B ------------------ |
| 2037 | #else |
| 2038 | // ----------- VARIANT A ------------------ |
| 2039 | BIGLOG = log(MPASQR / MCHSQR * (1.0 + BETA) * (1.0 + BETA) / 4.0 * |
| 2040 | pow(1.0 + MCHSQR / MPASQR, 2)); |
| 2041 | PRHARD = phocop.alpha / PI * (1.0 / BETA * BIGLOG) * |
| 2042 | (log(XPHMAX / phocop.xphcut) - .75 + phocop.xphcut / XPHMAX - .25 * phocop.xphcut * phocop.xphcut / XPHMAX / XPHMAX); |
| 2043 | PRHARD = PRHARD * PHOCHA(IDENT) * PHOCHA(IDENT) * phokey.fsec * phokey.fint; |
| 2044 | //me_channel_(&IDME); |
| 2045 | IDME = HEPEVT_struct::ME_channel; |
| 2046 | // write(*,*) 'KANALIK IDME=',IDME |
| 2047 | if (IDME == 0) { |
| 2048 | // do nothing |
| 2049 | } |
| 2050 | |
| 2051 | else if (IDME == 1) { |
| 2052 | PRHARD = PRHARD / (1.0 + 0.75 * phocop.alpha / PI); // NLO |
| 2053 | } else if (IDME == 2) { |
| 2054 | // work on virtual crrections in W decay to be done. |
| 2055 | } else { |
| 2056 | cout << "problem with ME_CHANNEL IDME= " << IDME << endl; |
| 2057 | exit(-1); |
| 2058 | } |
| 2059 | |
| 2060 | //----------- END OF VARIANT A ------------------ |
| 2061 | #endif |
| 2062 | if (phopro.irep == 0) phopro.probh = 0.0; |
| 2063 | PRKILL = 0.0; |
| 2064 | if (phokey.iexp) { // IEXP |
| 2065 | NCHAN = NCHAN + 1; |
| 2066 | if (phoexp.expini) { // EXPINI |
| 2067 | phoexp.pro[NCHAN - i] = PRHARD + 0.25 * (1.0 + phokey.fint); // we store hard photon emission prob |
| 2068 | //for leg NCHAN |
| 2069 | PRHARD = 0.0; // to kill emission at initialization call |
| 2070 | phopro.probh = PRHARD; |
| 2071 | } else { // EXPINI |
| 2072 | PRSUM = 0.0; |
| 2073 | for (K = NCHAN; K <= phoexp.NX; K++) PRSUM = PRSUM + phoexp.pro[K - i]; |
| 2074 | PRHARD = PRHARD / PRSUM; // note that PRHARD may be smaller than |
| 2075 | //phoexp.pro[NCHAN) because it is calculated |
| 2076 | // for kinematical configuartion as is |
| 2077 | // (with effects of previous photons) |
| 2078 | PRKILL = phoexp.pro[NCHAN - i] / PRSUM - PRHARD; |
| 2079 | |
| 2080 | } // EXPINI |
| 2081 | PRSOFT = 1.0 - PRHARD; |
| 2082 | } else { // IEXP |
| 2083 | PRHARD = PRHARD * PHOFAC(0); // PHOFAC is used to control eikonal |
| 2084 | // formfactors for non exp version only |
| 2085 | // here PHOFAC(0)=1 at least now. |
| 2086 | phopro.probh = PRHARD; |
| 2087 | } // IEXP |
| 2088 | PRSOFT = 1.0 - PRHARD; |
| 2089 | //-- |
| 2090 | //-- Check on kinematical bounds |
| 2091 | if (phokey.iexp) { |
| 2092 | if (PRSOFT < -5.0E-8) { |
| 2093 | DATA = PRSOFT; |
| 2094 | PHOERR(2, "PHOENE", DATA); |
| 2095 | } |
| 2096 | } else { |
| 2097 | if (PRSOFT < 0.1) { |
| 2098 | DATA = PRSOFT; |
| 2099 | PHOERR(2, "PHOENE", DATA); |
| 2100 | } |
| 2101 | } |
| 2102 | |
| 2103 | RRR = Photos::randomDouble(); |
| 2104 | if (RRR < PRSOFT) { |
| 2105 | //-- |
| 2106 | //-- No photon... (ie. photon too soft) |
| 2107 | XPHOTO = 0.0; |
| 2108 | if (RRR < PRKILL) XPHOTO = -5.0; //No photon...no further trials |
| 2109 | } else { |
| 2110 | //-- |
| 2111 | //-- Hard photon... (ie. photon hard enough). |
| 2112 | //-- Calculate Altarelli-Parisi Kernel |
| 2113 | do { |
| 2114 | XPHOTO = exp(Photos::randomDouble() * log(phocop.xphcut / XPHMAX)); |
| 2115 | XPHOTO = XPHOTO * XPHMAX; |
| 2116 | } while (Photos::randomDouble() > ((1.0 + pow(1.0 - XPHOTO / XPHMAX, 2)) / 2.0)); |
| 2117 | } |
| 2118 | |
| 2119 | //-- |
| 2120 | //-- Calculate parameter for PHOFAC function |
| 2121 | phopro.xf = 4.0 * MCHSQR * MPASQR / pow(MPASQR + MCHSQR - MNESQR, 2); |
| 2122 | return; |
| 2123 | } |
| 2124 | |
| 2125 | |
| 2126 | //---------------------------------------------------------------------- |
| 2127 | // |
| 2128 | // PHOTOS: Photon radiation in decays |
| 2129 | // |
| 2130 | // Purpose: Order (alpha) radiative corrections are generated in |
| 2131 | // the decay of the IPPAR-th particle in the HEP-like |
| 2132 | // common /PHOEVT/. Photon radiation takes place from one |
| 2133 | // of the charged daughters of the decaying particle IPPAR |
| 2134 | // WT is calculated, eventual rejection will be performed |
| 2135 | // later after inclusion of interference weight. |
| 2136 | // |
| 2137 | // Input Parameter: IPPAR: Pointer to decaying particle in |
| 2138 | // /PHOEVT/ and the common itself, |
| 2139 | // |
| 2140 | // Output Parameters: Common /PHOEVT/, either with or without a |
| 2141 | // photon(s) added. |
| 2142 | // WT weight of the configuration |
| 2143 | // |
| 2144 | // Author(s): Z. Was, B. van Eijk Created at: 26/11/89 |
| 2145 | // Last Update: 12/07/13 |
| 2146 | // |
| 2147 | //---------------------------------------------------------------------- |
| 2148 | |
| 2149 | void PHOPRE(int IPARR, double* pWT, int* pNEUDAU, int* pNCHARB) |
| 2150 | { |
| 2151 | int CHAPOI[pho.nmxhep]; |
| 2152 | double MINMAS, MPASQR, MCHREN; |
| 2153 | double EPS, DEL1, DEL2, DATA, BIGLOG; |
| 2154 | double MASSUM; |
| 2155 | int IP, IPPAR, I, J, ME, NCHARG, NEUPOI, NLAST; |
| 2156 | int IDABS; |
| 2157 | double WGT; |
| 2158 | int IDME; |
| 2159 | double a, b; |
| 2160 | double& WT = *pWT; |
| 2161 | int& NEUDAU = *pNEUDAU; |
| 2162 | int& NCHARB = *pNCHARB; |
| 2163 | double& COSTHG = phophs.costhg; |
| 2164 | double& SINTHG = phophs.sinthg; |
| 2165 | double& XPHOTO = phophs.xphoto; |
| 2166 | double& XPHMAX = phophs.xphmax; |
| 2167 | double* PNEUTR = phomom.pneutr; |
| 2168 | double& MCHSQR = phomom.mchsqr; |
| 2169 | double& MNESQR = phomom.mnesqr; |
| 2170 | |
| 2171 | static int i = 1; |
| 2172 | |
| 2173 | //-- |
| 2174 | IPPAR = IPARR; |
| 2175 | //-- Store pointers for cascade treatement... |
| 2176 | IP = IPPAR; |
| 2177 | NLAST = pho.nhep; |
| 2178 | |
| 2179 | //-- |
| 2180 | //-- Check decay multiplicity.. |
| 2181 | if (pho.jdahep[IP - i][1 - i] == 0) return; |
| 2182 | |
| 2183 | //-- |
| 2184 | //-- Loop over daughters, determine charge multiplicity |
| 2185 | |
| 2186 | NCHARG = 0; |
| 2187 | phopro.irep = 0; |
| 2188 | MINMAS = 0.0; |
| 2189 | MASSUM = 0.0; |
| 2190 | for (I = pho.jdahep[IP - i][1 - i]; I <= pho.jdahep[IP - i][2 - i]; I++) { |
| 2191 | //-- |
| 2192 | //-- |
| 2193 | //-- Exclude marked particles, quarks and gluons etc... |
| 2194 | IDABS = abs(pho.idhep[I - i]); |
| 2195 | if (pho.qedrad[I - pho.jdahep[IP - i][1 - i] + 3 - i]) { |
| 2196 | if (PHOCHA(pho.idhep[I - i]) != 0) { |
| 2197 | NCHARG = NCHARG + 1; |
| 2198 | if (NCHARG > pho.nmxhep) { |
| 2199 | DATA = NCHARG; |
| 2200 | PHOERR(1, "PHOTOS", DATA); |
| 2201 | } |
| 2202 | CHAPOI[NCHARG - i] = I; |
| 2203 | } |
| 2204 | MINMAS = MINMAS + pho.phep[I - i][5 - i] * pho.phep[I - i][5 - i]; |
| 2205 | } |
| 2206 | MASSUM = MASSUM + pho.phep[I - i][5 - i]; |
| 2207 | } |
| 2208 | |
| 2209 | if (NCHARG != 0) { |
| 2210 | //-- |
| 2211 | //-- Check that sum of daughter masses does not exceed parent mass |
| 2212 | if ((pho.phep[IP - i][5 - i] - MASSUM) / pho.phep[IP - i][5 - i] > 2.0 * phocop.xphcut) { |
| 2213 | //-- |
| 2214 | label30: |
| 2215 | |
| 2216 | // do{ |
| 2217 | |
| 2218 | for (J = 1; J <= 3; J++) PNEUTR[J - i] = -pho.phep[CHAPOI[NCHARG - i] - i][J - i]; |
| 2219 | PNEUTR[4 - i] = pho.phep[IP - i][5 - i] - pho.phep[CHAPOI[NCHARG - i] - i][4 - i]; |
| 2220 | //-- |
| 2221 | //-- Calculate invariant mass of 'neutral' etc. systems |
| 2222 | MPASQR = pho.phep[IP - i][5 - i] * pho.phep[IP - i][5 - i]; |
| 2223 | MCHSQR = pow(pho.phep[CHAPOI[NCHARG - i] - i][5 - i], 2); |
| 2224 | if ((pho.jdahep[IP - i][2 - i] - pho.jdahep[IP - i][1 - i]) == 1) { |
| 2225 | NEUPOI = pho.jdahep[IP - i][1 - i]; |
| 2226 | if (NEUPOI == CHAPOI[NCHARG - i]) NEUPOI = pho.jdahep[IP - i][2 - i]; |
| 2227 | MNESQR = pho.phep[NEUPOI - i][5 - i] * pho.phep[NEUPOI - i][5 - i]; |
| 2228 | PNEUTR[5 - i] = pho.phep[NEUPOI - i][5 - i]; |
| 2229 | } else { |
| 2230 | MNESQR = pow(PNEUTR[4 - i], 2) - pow(PNEUTR[1 - i], 2) - pow(PNEUTR[2 - i], 2) - pow(PNEUTR[3 - i], 2); |
| 2231 | MNESQR = max(MNESQR, MINMAS - MCHSQR); |
| 2232 | PNEUTR[5 - i] = sqrt(MNESQR); |
| 2233 | } |
| 2234 | |
| 2235 | //-- |
| 2236 | //-- Determine kinematical limit... |
| 2237 | XPHMAX = (MPASQR - pow(PNEUTR[5 - i] + pho.phep[CHAPOI[NCHARG - i] - i][5 - i], 2)) / MPASQR; |
| 2238 | |
| 2239 | //-- |
| 2240 | //-- Photon energy fraction... |
| 2241 | PHOENE(MPASQR, &MCHREN, &phwt.beta, &BIGLOG, pho.idhep[CHAPOI[NCHARG - i] - i]); |
| 2242 | //-- |
| 2243 | |
| 2244 | if (XPHOTO < -4.0) { |
| 2245 | NCHARG = 0; // we really stop trials |
| 2246 | XPHOTO = 0.0; // in this case !! |
| 2247 | //-- Energy fraction not too large (very seldom) ? Define angle. |
| 2248 | } else if ((XPHOTO < phocop.xphcut) || (XPHOTO > XPHMAX)) { |
| 2249 | //-- |
| 2250 | //-- No radiation was accepted, check for more daughters that may ra- |
| 2251 | //-- diate and correct radiation probability... |
| 2252 | NCHARG = NCHARG - 1; |
| 2253 | if (NCHARG > 0) phopro.irep = phopro.irep + 1; |
| 2254 | if (NCHARG > 0) goto label30; |
| 2255 | } else { |
| 2256 | //-- |
| 2257 | //-- Angle is generated in the frame defined by charged vector and |
| 2258 | //-- PNEUTR, distribution is taken in the infrared limit... |
| 2259 | EPS = MCHREN / (1.0 + phwt.beta); |
| 2260 | //-- |
| 2261 | //-- Calculate sin(theta) and cos(theta) from interval variables |
| 2262 | DEL1 = (2.0 - EPS) * pow(EPS / (2.0 - EPS), Photos::randomDouble()); |
| 2263 | DEL2 = 2.0 - DEL1; |
| 2264 | |
| 2265 | #ifdef VARIANTB |
| 2266 | // ----------- VARIANT B ------------------ |
| 2267 | // corrections for more efiicient interference correction, |
| 2268 | // instead of doubling crude distribution, we add flat parallel channel |
| 2269 | if (Photos::randomDouble() < BIGLOG / phwt.beta / (BIGLOG / phwt.beta + 2 * phokey.fint)) { |
| 2270 | COSTHG = (1.0 - DEL1) / phwt.beta; |
| 2271 | SINTHG = sqrt(DEL1 * DEL2 - MCHREN) / phwt.beta; |
| 2272 | } else { |
| 2273 | COSTHG = -1.0 + 2 * Photos::randomDouble(); |
| 2274 | SINTHG = sqrt(1.0 - COSTHG * COSTHG); |
| 2275 | } |
| 2276 | |
| 2277 | if (phokey.fint > 1.0) { |
| 2278 | |
| 2279 | WGT = 1.0 / (1.0 - phwt.beta * COSTHG); |
| 2280 | WGT = WGT / (WGT + phokey.fint); |
| 2281 | // WGT=1.0 // ?? |
| 2282 | } else { |
| 2283 | WGT = 1.0; |
| 2284 | } |
| 2285 | // |
| 2286 | // ----------- END OF VARIANT B ------------------ |
| 2287 | #else |
| 2288 | // ----------- VARIANT A ------------------ |
| 2289 | COSTHG = (1.0 - DEL1) / phwt.beta; |
| 2290 | SINTHG = sqrt(DEL1 * DEL2 - MCHREN) / phwt.beta; |
| 2291 | WGT = 1.0; |
| 2292 | // ----------- END OF VARIANT A ------------------ |
| 2293 | #endif |
| 2294 | //-- |
| 2295 | //-- Determine spin of particle and construct code for matrix element |
| 2296 | ME = (int)(2.0 * PHOSPI(pho.idhep[CHAPOI[NCHARG - i] - i]) + 1.0); |
| 2297 | //-- |
| 2298 | //-- Weighting procedure with 'exact' matrix element, reconstruct kine- |
| 2299 | //-- matics for photon, neutral and charged system and update /PHOEVT/. |
| 2300 | //-- Find pointer to the first component of 'neutral' system |
| 2301 | for (I = pho.jdahep[IP - i][1 - i]; I <= pho.jdahep[IP - i][2 - i]; I++) { |
| 2302 | if (I != CHAPOI[NCHARG - i]) { |
| 2303 | NEUDAU = I; |
| 2304 | goto label51; //break; // to 51 |
| 2305 | } |
| 2306 | } |
| 2307 | //-- |
| 2308 | //-- Pointer not found... |
| 2309 | DATA = NCHARG; |
| 2310 | PHOERR(5, "PHOKIN", DATA); |
| 2311 | label51: |
| 2312 | |
| 2313 | NCHARB = CHAPOI[NCHARG - i]; |
| 2314 | NCHARB = NCHARB - pho.jdahep[IP - i][1 - i] + 3; |
| 2315 | NEUDAU = NEUDAU - pho.jdahep[IP - i][1 - i] + 3; |
| 2316 | |
| 2317 | IDME = HEPEVT_struct::ME_channel; |
| 2318 | // two options introduced temporarily. |
| 2319 | // In future always PHOCOR-->PHOCORN |
| 2320 | // Tests and adjustment of wts for Znlo needed. |
| 2321 | // otherwise simple change. PHOCORN implements |
| 2322 | // exact ME for scalar to 2 scalar decays. |
| 2323 | if (IDME == 2) { |
| 2324 | b = PHOCORN(MPASQR, MCHREN, ME); |
| 2325 | WT = b * WGT; |
| 2326 | WT = WT / (1 - XPHOTO / XPHMAX + 0.5 * pow(XPHOTO / XPHMAX, 2)) * (1 - XPHOTO / XPHMAX) / 2; // factor to go to WnloWT |
| 2327 | } else if (IDME == 1) { |
| 2328 | |
| 2329 | a = PHOCOR(MPASQR, MCHREN, ME); |
| 2330 | b = PHOCORN(MPASQR, MCHREN, ME); |
| 2331 | WT = b * WGT ; |
| 2332 | WT = WT * phwt.wt1 * phwt.wt2 * phwt.wt3 / phocorwt.phocorwt1 / phocorwt.phocorwt2 / phocorwt.phocorwt3; // factor to go to ZnloWT |
| 2333 | // write(*,*) ' -----------' |
| 2334 | // write(*,*) phwt.wt1,' ',phwt.wt2,' ',phwt.wt3 |
| 2335 | // write(*,*) phocorwt.phocorwt1,' ',phocorwt.phocorwt2,' ',phocorwt.phocorwt3 |
| 2336 | } else { |
| 2337 | a = PHOCOR(MPASQR, MCHREN, ME); |
| 2338 | WT = a * WGT; |
| 2339 | // WT=b*WGT; // /(1-XPHOTO/XPHMAX+0.5*pow(XPHOTO/XPHMAX,2))*(1-XPHOTO/XPHMAX)/2; |
| 2340 | } |
| 2341 | |
| 2342 | |
| 2343 | |
| 2344 | } |
| 2345 | } else { |
| 2346 | DATA = pho.phep[IP - i][5 - i] - MASSUM; |
| 2347 | PHOERR(10, "PHOTOS", DATA); |
| 2348 | } |
| 2349 | } |
| 2350 | |
| 2351 | //-- |
| 2352 | return; |
| 2353 | } |
| 2354 | |
| 2355 | |
| 2356 | //---------------------------------------------------------------------- |
| 2357 | // |
| 2358 | // PHOMAK: PHOtos MAKe |
| 2359 | // |
| 2360 | // Purpose: Single or double bremstrahlung radiative corrections |
| 2361 | // are generated in the decay of the IPPAR-th particle in |
| 2362 | // the HEP common /PH_HEPEVT/. Example of the use of |
| 2363 | // general tools. |
| 2364 | // |
| 2365 | // Input Parameter: IPPAR: Pointer to decaying particle in |
| 2366 | // /PH_HEPEVT/ and the common itself |
| 2367 | // |
| 2368 | // Output Parameters: Common /PH_HEPEVT/, either with or without |
| 2369 | // particles added. |
| 2370 | // |
| 2371 | // Author(s): Z. Was, Created at: 26/05/93 |
| 2372 | // Last Update: 29/01/05 |
| 2373 | // |
| 2374 | //---------------------------------------------------------------------- |
| 2375 | |
| 2376 | void PHOMAK(int IPPAR, int NHEP0) |
| 2377 | { |
| 2378 | |
| 2379 | double DATA; |
| 2380 | int IP, NCHARG, IDME; |
| 2381 | int IDUM; |
| 2382 | int NCHARB, NEUDAU; |
| 2383 | double RN, WT; |
| 2384 | bool BOOST; |
| 2385 | static int i = 1; |
| 2386 | //-- |
| 2387 | IP = IPPAR; |
| 2388 | IDUM = 1; |
| 2389 | NCHARG = 0; |
| 2390 | //-- |
| 2391 | PHOIN(IP, &BOOST, &NHEP0); |
| 2392 | PHOCHK(hep.jdahep[IP - i][1 - i]); |
| 2393 | WT = 0.0; |
| 2394 | PHOPRE(1, &WT, &NEUDAU, &NCHARB); |
| 2395 | |
| 2396 | if (WT == 0.0) return; |
| 2397 | RN = Photos::randomDouble(); |
| 2398 | // PHODO is caling randomDouble(), thus change of series if it is moved before if |
| 2399 | PHODO(1, NCHARB, NEUDAU); |
| 2400 | |
| 2401 | #ifdef VARIANTB |
| 2402 | // we eliminate divisions /phokey.fint in variant B. ??? |
| 2403 | #endif |
| 2404 | // get ID of channel dependent ME, ID=0 means no |
| 2405 | |
| 2406 | IDME = HEPEVT_struct::ME_channel; |
| 2407 | // corrections for matrix elements |
| 2408 | // controlled by IDME |
| 2409 | // write(*,*) 'KANALIK IDME=',IDME |
| 2410 | |
| 2411 | if (IDME == 0) { // default |
| 2412 | |
| 2413 | if (phokey.interf) WT = WT * PHINT(IDUM); |
| 2414 | if (phokey.ifw) PHOBW(&WT); // extra weight for leptonic W decay |
| 2415 | } else if (IDME == 2) { // ME weight for leptonic W decay |
| 2416 | |
| 2417 | PhotosMEforW::PHOBWnlo(&WT); |
| 2418 | WT = WT * 2.0; |
| 2419 | } else if (IDME == 1) { // ME weight for leptonic Z decay |
| 2420 | |
| 2421 | WT = WT * PhotosMEforZ::phwtnlo(); |
| 2422 | } else { |
| 2423 | cout << "problem with ME_CHANNEL IDME= " << IDME << endl; |
| 2424 | exit(-1); |
| 2425 | } |
| 2426 | |
| 2427 | #ifndef VARIANTB |
| 2428 | WT = WT / phokey.fint; // FINT must be in variant A |
| 2429 | #endif |
| 2430 | |
| 2431 | DATA = WT; |
| 2432 | if (WT > 1.0) PHOERR(3, "WT_INT", DATA); |
| 2433 | // weighting |
| 2434 | if (RN <= WT) { |
| 2435 | PHOOUT(IP, BOOST, NHEP0); |
| 2436 | } |
| 2437 | return; |
| 2438 | } |
| 2439 | |
| 2440 | //---------------------------------------------------------------------- |
| 2441 | // |
| 2442 | // PHTYPE: Central manadgement routine. |
| 2443 | // |
| 2444 | // Purpose: defines what kind of the |
| 2445 | // actions will be performed at point ID. |
| 2446 | // |
| 2447 | // Input Parameters: ID: pointer of particle starting branch |
| 2448 | // in /PH_HEPEVT/ to be treated. |
| 2449 | // |
| 2450 | // Output Parameters: Common /PH_HEPEVT/. |
| 2451 | // |
| 2452 | // Author(s): Z. Was Created at: 24/05/93 |
| 2453 | // P. Golonka Last Update: 27/06/04 |
| 2454 | // |
| 2455 | //---------------------------------------------------------------------- |
| 2456 | void PHTYPE(int ID) |
| 2457 | { |
| 2458 | |
| 2459 | int K; |
| 2460 | double PRSUM, ESU; |
| 2461 | int NHEP0; |
| 2462 | bool IPAIR; |
| 2463 | bool IPHOT; |
| 2464 | double RN, SUM; |
| 2465 | bool IFOUR; |
| 2466 | int& NCHAN = phoexp.nchan; |
| 2467 | |
| 2468 | static int i = 1; |
| 2469 | |
| 2470 | |
| 2471 | //-- |
| 2472 | IFOUR = phokey.itre; // we can make internal choice whether |
| 2473 | // we want 3 or four photons at most. |
| 2474 | IPAIR = false; |
| 2475 | IPAIR = Photos::IfPair; |
| 2476 | IPHOT = true; |
| 2477 | IPHOT = Photos::IfPhot; |
| 2478 | //-- Check decay multiplicity.. |
| 2479 | if (hep.jdahep[ID - i][1 - i] == 0) return; |
| 2480 | // if (hep.jdahep[ID-i][1-i]==hep.jdahep[ID-i][2-i]) return; |
| 2481 | //-- |
| 2482 | NHEP0 = hep.nhep; |
| 2483 | |
| 2484 | // initialization of pho.qedrad for new event. |
| 2485 | // some of (old style and doubling) restrictions introduced with PHOCHK, |
| 2486 | // also new pairs have emissions blocked with pho.qedrad[] |
| 2487 | // most of the restrictions are introduced prior decay vertex is copied |
| 2488 | // to struct pho. |
| 2489 | |
| 2490 | // Establish size for the struct pho: number of daughters + 2 places for mothers (no grandmothers) |
| 2491 | // This solution is `hep.nhep hardy'. Use of hep.nhep was perilous |
| 2492 | // if decaying particle (ID-i) was the first in the event. That was the case of EvtGen |
| 2493 | // interface. We adopt to such non-standard HepMC fill. |
| 2494 | // NOTE: here 'max' is used as a safety for future changes to hep or pho content. |
| 2495 | // TP ZW (26.09.15): Thanks to Michal Kreps and John Back |
| 2496 | |
| 2497 | int pho_size = max(NHEP0, (hep.jdahep[ID - i][2 - i] - hep.jdahep[ID - i][1 - i] + 1) + 2); |
| 2498 | |
| 2499 | for (int I = 0; I < pho_size; ++I) { |
| 2500 | pho.qedrad[I] = true; |
| 2501 | } |
| 2502 | |
| 2503 | double elMass = 0.000511; |
| 2504 | double muMass = 0.1057; |
| 2505 | double STRENG = 0.5; |
| 2506 | |
| 2507 | if (IPAIR) { |
| 2508 | |
| 2509 | switch (Photos::momentumUnit) { |
| 2510 | case Photos::GEV: |
| 2511 | elMass = 0.000511; |
| 2512 | muMass = 0.1057; |
| 2513 | break; |
| 2514 | case Photos::MEV: |
| 2515 | elMass = 0.511; |
| 2516 | muMass = 105.7; |
| 2517 | break; |
| 2518 | default: |
| 2519 | Log::Error() << "GEV or MEV unit must be set for pair emission" << endl; |
| 2520 | break; |
| 2521 | }; |
| 2522 | PHOPAR(ID, NHEP0, 11, elMass, &STRENG); |
| 2523 | PHOPAR(ID, NHEP0, 13, muMass, &STRENG); |
| 2524 | } |
| 2525 | //-- |
| 2526 | if (IPHOT) { |
| 2527 | if (phokey.iexp) { |
| 2528 | phoexp.expini = true; // Initialization/cleaning |
| 2529 | for (NCHAN = 1; NCHAN <= phoexp.NX; NCHAN++) |
| 2530 | phoexp.pro[NCHAN - i] = 0.0; |
| 2531 | NCHAN = 0; |
| 2532 | |
| 2533 | phokey.fsec = 1.0; |
| 2534 | PHOMAK(ID, NHEP0); // Initialization/crude formfactors into |
| 2535 | // phoexp.pro[NCHAN) |
| 2536 | phoexp.expini = false; |
| 2537 | RN = Photos::randomDouble(); |
| 2538 | PRSUM = 0.0; |
| 2539 | for (K = 1; K <= phoexp.NX; K++)PRSUM = PRSUM + phoexp.pro[K - i]; |
| 2540 | |
| 2541 | ESU = exp(-PRSUM); |
| 2542 | // exponent for crude Poissonian multiplicity |
| 2543 | // distribution, will be later overwritten |
| 2544 | // to give probability for k |
| 2545 | SUM = ESU; |
| 2546 | // distribuant for the crude Poissonian |
| 2547 | // at first for k=0 |
| 2548 | for (K = 1; K <= 100; K++) { // hard coded max (photon) multiplicity is 100 |
| 2549 | if (RN < SUM) break; |
| 2550 | ESU = ESU * PRSUM / K; // we get at K ESU=EXP(-PRSUM)*PRSUM**K/K! |
| 2551 | SUM = SUM + ESU; // thus we get distribuant at K. |
| 2552 | NCHAN = 0; |
| 2553 | PHOMAK(ID, NHEP0); // LOOPING |
| 2554 | if (SUM > 1.0 - phokey.expeps) break; |
| 2555 | } |
| 2556 | |
| 2557 | } else if (IFOUR) { |
| 2558 | //-- quatro photon emission |
| 2559 | phokey.fsec = 1.0; |
| 2560 | RN = Photos::randomDouble(); |
| 2561 | if (RN >= 23.0 / 24.0) { |
| 2562 | PHOMAK(ID, NHEP0); |
| 2563 | PHOMAK(ID, NHEP0); |
| 2564 | PHOMAK(ID, NHEP0); |
| 2565 | PHOMAK(ID, NHEP0); |
| 2566 | } else if (RN >= 17.0 / 24.0) { |
| 2567 | PHOMAK(ID, NHEP0); |
| 2568 | PHOMAK(ID, NHEP0); |
| 2569 | } else if (RN >= 9.0 / 24.0) { |
| 2570 | PHOMAK(ID, NHEP0); |
| 2571 | } else { |
| 2572 | } |
| 2573 | } else if (phokey.itre) { |
| 2574 | //-- triple photon emission |
| 2575 | phokey.fsec = 1.0; |
| 2576 | RN = Photos::randomDouble(); |
| 2577 | if (RN >= 5.0 / 6.0) { |
| 2578 | PHOMAK(ID, NHEP0); |
| 2579 | PHOMAK(ID, NHEP0); |
| 2580 | PHOMAK(ID, NHEP0); |
| 2581 | } else if (RN >= 2.0 / 6.0) { |
| 2582 | PHOMAK(ID, NHEP0); |
| 2583 | } |
| 2584 | } else if (phokey.isec) { |
| 2585 | //-- double photon emission |
| 2586 | phokey.fsec = 1.0; |
| 2587 | RN = Photos::randomDouble(); |
| 2588 | if (RN >= 0.5) { |
| 2589 | PHOMAK(ID, NHEP0); |
| 2590 | PHOMAK(ID, NHEP0); |
| 2591 | } |
| 2592 | } else { |
| 2593 | //-- single photon emission |
| 2594 | phokey.fsec = 1.0; |
| 2595 | PHOMAK(ID, NHEP0); |
| 2596 | } |
| 2597 | } |
| 2598 | //-- |
| 2599 | //-- lepton anti-lepton pair(s) |
| 2600 | // we prepare to migrate half of tries to before photons accordingly to LL |
| 2601 | // pho.qedrad is not yet used by PHOPAR |
| 2602 | if (IPAIR) { |
| 2603 | PHOPAR(ID, NHEP0, 11, elMass, &STRENG); |
| 2604 | PHOPAR(ID, NHEP0, 13, muMass, &STRENG); |
| 2605 | } |
| 2606 | |
| 2607 | // Fill Photos::EventNo in user main program to provide |
| 2608 | // debug input for specific events, e.g.: |
| 2609 | // if(Photos::EventNo==1331094) printf("PHOTOS: event no: %10i finished\n",Photos::EventNo); |
| 2610 | } |
| 2611 | |
| 2612 | /*---------------------------------------------------------------------- |
| 2613 | |
| 2614 | PHOTOS: Photon radiation in decays |
| 2615 | |
| 2616 | Purpose: e+e- pairs are generated in |
| 2617 | the decay of the IPPAR-th particle in the HEP-like |
| 2618 | common /PHOEVT/. Radiation takes place from one |
| 2619 | of the charged daughters of the decaying particle IPPAR |
| 2620 | |
| 2621 | |
| 2622 | |
| 2623 | Input Parameter: IPPAR: Pointer to decaying particle in |
| 2624 | /PHOEVT/ and the common itself, |
| 2625 | NHEP0 length of the /HEPEVT/ entry |
| 2626 | before starting any activity on this |
| 2627 | IPPAR decay. |
| 2628 | Output Parameters: Common /HEPEVT/, either with or without a |
| 2629 | e+e-(s) added. |
| 2630 | |
| 2631 | |
| 2632 | Author(s): Z. Was, Created at: 01/06/93 |
| 2633 | Last Update: |
| 2634 | |
| 2635 | ----------------------------------------------------------------------*/ |
| 2636 | void PHOPAR(int IPARR, int NHEP0, int idlep, double masslep, double* pSTRENG) |
| 2637 | { |
| 2638 | double PCHAR[4], PNEU[4], PELE[4], PPOZ[4], BUF[4]; |
| 2639 | int IP, IPPAR, NLAST; |
| 2640 | bool BOOST, JESLI; |
| 2641 | static int i = 1; |
| 2642 | IPPAR = IPARR; |
| 2643 | |
| 2644 | double& STRENG = *pSTRENG; |
| 2645 | |
| 2646 | // Store pointers for cascade treatment... |
| 2647 | IP = 0; |
| 2648 | NLAST = pho.nhep; |
| 2649 | // Check decay multiplicity.. |
| 2650 | PHOIN(IPPAR, &BOOST, &NHEP0); |
| 2651 | PHOCHK(pho.jdahep[IP][0]); // should be loop over all mothers? |
| 2652 | PHLUPA(100); |
| 2653 | |
| 2654 | if (pho.jdahep[IP][0] == 0) return; |
| 2655 | if (pho.jdahep[IP][0] == pho.jdahep[IP][1]) return; |
| 2656 | |
| 2657 | // Loop over charged daughters |
| 2658 | for (int I = pho.jdahep[IP][0] - i; I <= pho.jdahep[IP][1] - i; ++I) { |
| 2659 | |
| 2660 | |
| 2661 | // Skip this particle if it has no charge |
| 2662 | if (PHOCHA(pho.idhep[I]) == 0) continue; |
| 2663 | |
| 2664 | int IDABS = abs(pho.idhep[I]); |
| 2665 | // at the moment the following re-definition make not much sense as constraints |
| 2666 | // were already checked before for photons tries. |
| 2667 | // we have to come back to this when we will have pairs emitted before photons. |
| 2668 | |
| 2669 | pho.qedrad[I] = pho.qedrad[I] && F(0, IDABS) && F(0, abs(pho.idhep[1 - i])) |
| 2670 | && (pho.idhep[2 - i] == 0); |
| 2671 | |
| 2672 | if (!pho.qedrad[I]) continue; // |
| 2673 | |
| 2674 | |
| 2675 | // Set 3-vectors |
| 2676 | for (int J = 0; J < 3; ++J) { |
| 2677 | PCHAR[J] = pho.phep[I][J]; |
| 2678 | PNEU [J] = -pho.phep[I][J]; |
| 2679 | } |
| 2680 | |
| 2681 | // Set energy |
| 2682 | PNEU[3] = pho.phep[IP][3] - pho.phep[I][3]; |
| 2683 | PCHAR[3] = pho.phep[I][3]; |
| 2684 | // Set mass |
| 2685 | double AMCH = pho.phep[I][4]; |
| 2686 | //here we attempt generating pair from PCHAR. One of the charged |
| 2687 | //decay products; that is why algorithm works in a loop. |
| 2688 | //PNEU is four vector of all decay products except PCHAR |
| 2689 | //we do not care on rare cases when two pairs could be generated |
| 2690 | //we assume it is negligibly rare and fourth order in alpha anyway |
| 2691 | //TRYPAR should take as an input electron mass. |
| 2692 | //then it can be used for muons. |
| 2693 | // printf ("wrotki %10.7f\n",STRENG); |
| 2694 | /* |
| 2695 | double PCH[4]={0}; |
| 2696 | double PNEu[4]={0}; |
| 2697 | double CC1; |
| 2698 | double CC2; |
| 2699 | |
| 2700 | for(int K = 0; K<4; ++K) { |
| 2701 | PCH[K]=PCHAR[K]; |
| 2702 | PNEu[K]=PNEU[K]; |
| 2703 | } |
| 2704 | */ |
| 2705 | // printf ("idlep,pdgidid= %10i %10i\n",idlep,pho.idhep[I]); |
| 2706 | |
| 2707 | // arrangements for the case when emitted lept6ons have |
| 2708 | // the same flavour as emitters |
| 2709 | bool sameflav = abs(idlep) == abs(pho.idhep[I]); |
| 2710 | int idsign = 1; |
| 2711 | if (pho.idhep[I] < 0) idsign = -1; // this is to ensure |
| 2712 | //that first lepton has the same PDGID as emitter |
| 2713 | |
| 2714 | trypar(&JESLI, &STRENG, AMCH, masslep, PCHAR, PNEU, PELE, PPOZ, &sameflav); |
| 2715 | // printf ("rowerek %10.7f\n",STRENG); |
| 2716 | //emitted pair four momenta are stored in PELE PPOZ |
| 2717 | //then JESLI=.true. |
| 2718 | /* |
| 2719 | if (JESLI) { |
| 2720 | // printf ("PCHAR %10.7f %10.7f %10.7f %10.7f\n",PCHAR[0],PCHAR[1],PCHAR[2],PCHAR[3]); |
| 2721 | //printf ("PNEU %10.7f %10.7f %10.7f %10.7f\n",PNEU[0],PNEU[1],PNEU[2],PNEU[3]); |
| 2722 | |
| 2723 | // printf ("PNEu %10.7f %10.7f %10.7f %10.7f\n",PNEu[0],PNEu[1],PNEu[2],PNEu[3]); |
| 2724 | |
| 2725 | printf ("PELE %10.7f %10.7f %10.7f %10.7f\n",PELE[0],PELE[1],PELE[2],PELE[3]); |
| 2726 | printf ("PPOZ %10.7f %10.7f %10.7f %10.7f\n",PPOZ[0],PPOZ[1],PPOZ[2],PPOZ[3]); |
| 2727 | printf ("-----------------\n"); |
| 2728 | printf ("PCH %10.7f %10.7f %10.7f %10.7f\n",PCH[0],PCH[1],PCH[2],PCH[3]); |
| 2729 | CC1=(PELE[0]*PCHAR[0]+PELE[1]*PCHAR[1]+PELE[2]*PCHAR[2])/sqrt(PELE[0]*PELE[0]+PELE[1]*PELE[1]+PELE[2]*PELE[2])/sqrt(PCHAR[0]*PCHAR[0]+PCHAR[1]*PCHAR[1]+PCHAR[2]*PCHAR[2]); |
| 2730 | CC2=(PPOZ[0]*PCHAR[0]+PPOZ[1]*PCHAR[1]+PPOZ[2]*PCHAR[2])/sqrt(PPOZ[0]*PPOZ[0]+PPOZ[1]*PPOZ[1]+PPOZ[2]*PPOZ[2])/sqrt(PCHAR[0]*PCHAR[0]+PCHAR[1]*PCHAR[1]+PCHAR[2]*PCHAR[2]); |
| 2731 | |
| 2732 | printf ("-=================-\n"); |
| 2733 | |
| 2734 | } |
| 2735 | */ |
| 2736 | // If JESLI = true, we modify old particles of the vertex |
| 2737 | if (JESLI) { |
| 2738 | PHLUPA(1010); |
| 2739 | |
| 2740 | // we have to correct 4-momenta |
| 2741 | // of all decay products |
| 2742 | // we use PARTRA for that |
| 2743 | // PELE PPOZ are in right frame |
| 2744 | for (int J = pho.jdahep[IP][0] - i; J <= pho.jdahep[IP][1] - i; ++J) { |
| 2745 | for (int K = 0; K < 4; ++K) { |
| 2746 | BUF[K] = pho.phep[J][K]; |
| 2747 | } |
| 2748 | if (J == I) partra(1, BUF); |
| 2749 | else partra(-1, BUF); |
| 2750 | |
| 2751 | for (int K = 0; K < 4; ++K) { |
| 2752 | pho.phep[J][K] = BUF[K]; |
| 2753 | } |
| 2754 | /* |
| 2755 | if (J == I){ |
| 2756 | printf ("PCHar %10.7f %10.7f %10.7f %10.7f\n",pho.phep[J][0],pho.phep[J][1],pho.phep[J][2],pho.phep[J][3]); |
| 2757 | printf ("c1= %10.7f\n",CC1); |
| 2758 | printf ("c2= %10.7f\n",CC2); |
| 2759 | printf ("-=#####################################====-\n"); |
| 2760 | } |
| 2761 | */ |
| 2762 | } |
| 2763 | |
| 2764 | PHLUPA(1011); |
| 2765 | |
| 2766 | if (darkr.ifspecial == 1) { |
| 2767 | // virtual adding to vertex |
| 2768 | pho.nhep = pho.nhep + 1; |
| 2769 | pho.isthep[pho.nhep - i] = 2; |
| 2770 | pho.idhep [pho.nhep - i] = darkr.IDspecial; |
| 2771 | pho.jmohep[pho.nhep - i][0] = IP; |
| 2772 | pho.jmohep[pho.nhep - i][1] = 0; |
| 2773 | pho.jdahep[pho.nhep - i][0] = 0; |
| 2774 | pho.jdahep[pho.nhep - i][1] = 0; |
| 2775 | pho.qedrad[pho.nhep - i] = false; |
| 2776 | |
| 2777 | pho.phep[pho.nhep - i][4] = sqrt(-(PELE[0] + PPOZ[0]) * (PELE[0] + PPOZ[0]) |
| 2778 | - (PELE[1] + PPOZ[1]) * (PELE[1] + PPOZ[1]) |
| 2779 | - (PELE[2] + PPOZ[2]) * (PELE[2] + PPOZ[2]) |
| 2780 | + (PELE[3] + PPOZ[3]) * (PELE[3] + PPOZ[3])); |
| 2781 | } |
| 2782 | //double life=darkr.SpecialLife; |
| 2783 | double life = darkr.SpecialLife * (-log(Photos::randomDouble())); |
| 2784 | // life = pho.phep[pho.nhep-i][4]; // this is helpful for test. Instead of vertex we get momentum |
| 2785 | // here was missing if(darkr.ifspecial==1) zbw:16.09.2021 |
| 2786 | if (darkr.ifspecial == 1) { |
| 2787 | for (int K = 1; K <= 4; ++K) { |
| 2788 | pho.phep[pho.nhep - i][K - i] = PELE[K - i] + PPOZ[K - i]; |
| 2789 | pho.vhep[pho.nhep - i][K - i] = pho.vhep[IP][K - i] |
| 2790 | + (PELE[K - i] + PPOZ[K - i]) / pho.phep[pho.nhep - i][4] * life; |
| 2791 | } |
| 2792 | } |
| 2793 | |
| 2794 | // electron: adding to vertex |
| 2795 | pho.nhep = pho.nhep + 1; |
| 2796 | pho.isthep[pho.nhep - i] = 1; |
| 2797 | pho.idhep [pho.nhep - i] = idlep * idsign; |
| 2798 | pho.jmohep[pho.nhep - i][0] = IP; |
| 2799 | pho.jmohep[pho.nhep - i][1] = 0; |
| 2800 | pho.jdahep[pho.nhep - i][0] = 0; |
| 2801 | pho.jdahep[pho.nhep - i][1] = 0; |
| 2802 | pho.qedrad[pho.nhep - i] = false; |
| 2803 | |
| 2804 | |
| 2805 | for (int K = 1; K <= 4; ++K) { |
| 2806 | pho.phep[pho.nhep - i][K - i] = PELE[K - i]; |
| 2807 | if (darkr.ifspecial == 1) |
| 2808 | pho.vhep[pho.nhep - i][K - i] = pho.vhep[pho.nhep - i - 1][K - i]; |
| 2809 | } |
| 2810 | |
| 2811 | pho.phep[pho.nhep - i][4] = masslep; |
| 2812 | |
| 2813 | // positron: adding |
| 2814 | pho.nhep = pho.nhep + 1; |
| 2815 | pho.isthep[pho.nhep - i] = 1; |
| 2816 | pho.idhep [pho.nhep - i] = -idlep * idsign; |
| 2817 | pho.jmohep[pho.nhep - i][0] = IP; |
| 2818 | pho.jmohep[pho.nhep - i][1] = 0; |
| 2819 | pho.jdahep[pho.nhep - i][0] = 0; |
| 2820 | pho.jdahep[pho.nhep - i][1] = 0; |
| 2821 | pho.qedrad[pho.nhep - i] = false; |
| 2822 | |
| 2823 | for (int K = 1; K <= 4; ++K) { |
| 2824 | pho.phep[pho.nhep - i][K - i] = PPOZ[K - i]; |
| 2825 | if (darkr.ifspecial == 1) |
| 2826 | pho.vhep[pho.nhep - i][K - i] = pho.vhep[pho.nhep - i - 1][K - i]; |
| 2827 | } |
| 2828 | |
| 2829 | // for mc-test with KORALW, mumu from mu mu emissions: BEGIN |
| 2830 | /* |
| 2831 | double RRX[2]; |
| 2832 | for( int k=0;k<=1;k++) RRX[k]=Photos::randomDouble(); |
| 2833 | |
| 2834 | for(int KK=0;KK<=pho.nhep-i;KK++){ |
| 2835 | for(int KJ=KK+1;KJ<=pho.nhep-i;KJ++){ |
| 2836 | // 1 <-> 3 |
| 2837 | if(RRX[0]>.5&&pho.idhep[KK]==13&&pho.idhep[KJ]==13){ |
| 2838 | for( int k=0;k<=3;k++){ |
| 2839 | double stored=pho.phep[KK][k]; |
| 2840 | pho.phep[KK][k]=pho.phep[KJ][k]; |
| 2841 | pho.phep[KJ][k]=stored; |
| 2842 | } |
| 2843 | } |
| 2844 | // 2 <-> 4 |
| 2845 | |
| 2846 | if(RRX[1]>.5&&pho.idhep[KK]==-13&&pho.idhep[KJ]==-13){ |
| 2847 | for( int k=0;k<=3;k++){ |
| 2848 | double stored=pho.phep[KK][k]; |
| 2849 | pho.phep[KK][k]=pho.phep[KJ][k]; |
| 2850 | pho.phep[KJ][k]=stored; |
| 2851 | |
| 2852 | } |
| 2853 | } |
| 2854 | |
| 2855 | } |
| 2856 | } |
| 2857 | |
| 2858 | // for mc-test with KORALW, mumu from mu mu emissions: END |
| 2859 | */ |
| 2860 | |
| 2861 | pho.phep[pho.nhep - i][4] = masslep; |
| 2862 | PHCORK(0); |
| 2863 | // write in |
| 2864 | PHLUPA(1012); |
| 2865 | PHOOUT(IPPAR, BOOST, NHEP0); |
| 2866 | PHOIN(IPPAR, &BOOST, &NHEP0); |
| 2867 | PHLUPA(1013); |
| 2868 | } // end of if (JESLI) |
| 2869 | } // end of loop over charged particles |
| 2870 | } |
| 2871 | |
| 2872 | |
| 2873 | } // namespace Photospp |
| 2874 |