| File: | generators/kkmc/photospp/forZ-MEc.cc |
| Warning: | line 102, column 5 Value stored to 'RaZ' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | #include "forZ-MEc.h" |
| 2 | #include "Photos.h" |
| 3 | #include "PhotosUtilities.h" |
| 4 | #include "photosC.h" |
| 5 | #include <cmath> |
| 6 | #include <cstdio> |
| 7 | #include <cstdlib> |
| 8 | #include <iostream> |
| 9 | using std::cout; |
| 10 | using std::endl; |
| 11 | using namespace Photospp; |
| 12 | using namespace PhotosUtilities; |
| 13 | |
| 14 | namespace Photospp { |
| 15 | |
| 16 | double (*PhotosMEforZ::currentVakPol)(double[4], double[4], double[4], double[4], double[4], int, int) = default_VakPol; |
| 17 | |
| 18 | // from photosC.cxx |
| 19 | |
| 20 | void PHODMP(); |
| 21 | double PHINT(int idumm); |
| 22 | // ---------------------------------------------------------------------- |
| 23 | // PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION |
| 24 | // IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK |
| 25 | // NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE |
| 26 | // IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY) |
| 27 | // SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF) |
| 28 | // AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE |
| 29 | // KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS |
| 30 | // |
| 31 | // called by : EVENTE, EVENTM, FUNTIH, ..... |
| 32 | // ---------------------------------------------------------------------- |
| 33 | |
| 34 | void PhotosMEforZ::GIVIZO(int IDFERM, int IHELIC, double* SIZO3, double* CHARGE, int* KOLOR) |
| 35 | { |
| 36 | // |
| 37 | int IH, IDTYPE, IC, LEPQUA, IUPDOW; |
| 38 | if (IDFERM == 0 || abs(IDFERM) > 4 || abs(IHELIC) != 1) { |
| 39 | cout << "STOP IN GIVIZO: WRONG PARAMS" << endl; |
| 40 | exit(-1); |
| 41 | } |
| 42 | |
| 43 | IH = IHELIC; |
| 44 | IDTYPE = abs(IDFERM); |
| 45 | IC = IDFERM / IDTYPE; |
| 46 | LEPQUA = (int)(IDTYPE * 0.4999999); |
| 47 | IUPDOW = IDTYPE - 2 * LEPQUA - 1; |
| 48 | *CHARGE = (-IUPDOW + 2.0 / 3.0 * LEPQUA) * IC; |
| 49 | *SIZO3 = 0.25 * (IC - IH) * (1 - 2 * IUPDOW); |
| 50 | *KOLOR = 1 + 2 * LEPQUA; |
| 51 | //** NOTE THAT CONVENTIONALY Z0 COUPLING IS |
| 52 | //** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ)) |
| 53 | return; |
| 54 | } |
| 55 | |
| 56 | |
| 57 | //////////////////////////////////////////////////////////////////////////// |
| 58 | /// // |
| 59 | /// This routine provides unsophisticated Born differential cross section // |
| 60 | /// at the crude x-section level, with Z and gamma s-chanel exchange. // |
| 61 | /////////////////////////////////////////////////////////////////////////// |
| 62 | double PhotosMEforZ::PHBORNM(double svar, double costhe, double T3e, double qe, double T3f, double qf, int NCf) |
| 63 | { |
| 64 | |
| 65 | double s, Sw2, MZ, MZ2, GammZ, AlfInv, GFermi; // t,MW,MW2, |
| 66 | double Ve, Ae, thresh; // sum,deno, |
| 67 | double xe, yf, xf, ye, ff0, ff1, amx2, amfin, Vf, Af; |
| 68 | double ReChiZ, SqChiZ, RaZ; //,RaW,ReChiW,SqChiW; |
| 69 | double Born; //, BornS; |
| 70 | // int KeyZet,HadMin,KFbeam; |
| 71 | // int i,ke,KFfin,kf,IsGenerated,iKF; |
| 72 | int KeyWidFix; |
| 73 | |
| 74 | AlfInv = 137.0359895; |
| 75 | GFermi = 1.16639e-5; |
| 76 | |
| 77 | //-------------------------------------------------------------------- |
| 78 | s = svar; |
| 79 | //------------------------------ |
| 80 | // EW paratemetrs taken from BornV |
| 81 | MZ = 91.187; |
| 82 | GammZ = 2.50072032; |
| 83 | Sw2 = .22276773; |
| 84 | //------------------------------ |
| 85 | // Z and gamma couplings to beams (electrons) |
| 86 | // Z and gamma couplings to final fermions |
| 87 | // Loop over all flavours defined in m_xpar(400+i) |
| 88 | |
| 89 | |
| 90 | //------ incoming fermion |
| 91 | Ve = 2 * T3e - 4 * qe * Sw2; |
| 92 | Ae = 2 * T3e; |
| 93 | //------ final fermion couplings |
| 94 | amfin = 0.000511; // m_xpar(kf+6) |
| 95 | Vf = 2 * T3f - 4 * qf * Sw2; |
| 96 | Af = 2 * T3f; |
| 97 | if (fabs(costhe) > 1.0) { |
| 98 | cout << "+++++STOP in PHBORN: costhe>0 =" << costhe << endl; |
| 99 | exit(-1); |
| 100 | } |
| 101 | MZ2 = MZ * MZ; |
| 102 | RaZ = (GFermi * MZ2 * AlfInv) / (sqrt(2.0) * 8.0 * PI); // |
Value stored to 'RaZ' is never read | |
| 103 | RaZ = 1 / (16.0 * Sw2 * (1.0 - Sw2)); |
| 104 | KeyWidFix = 1; // fixed width |
| 105 | KeyWidFix = 0; // variable width |
| 106 | if (KeyWidFix == 0) { |
| 107 | ReChiZ = (s - MZ2) * s / ((s - MZ2) * (s - MZ2) + (GammZ * s / MZ) * (GammZ * s / MZ)) * RaZ; // variable width |
| 108 | SqChiZ = s * s / ((s - MZ2) * (s - MZ2) + (GammZ * s / MZ) * (GammZ * s / MZ)) * RaZ * RaZ; // variable width |
| 109 | } else { |
| 110 | ReChiZ = (s - MZ2) * s / ((s - MZ2) * (s - MZ2) + (GammZ * MZ) * (GammZ * MZ)) * RaZ; // fixed width |
| 111 | SqChiZ = s * s / ((s - MZ2) * (s - MZ2) + (GammZ * MZ) * (GammZ * MZ)) * RaZ * RaZ; // fixed width |
| 112 | } |
| 113 | xe = Ve * Ve + Ae * Ae; |
| 114 | xf = Vf * Vf + Af * Af; |
| 115 | ye = 2 * Ve * Ae; |
| 116 | yf = 2 * Vf * Af; |
| 117 | ff0 = qe * qe * qf * qf + 2 * ReChiZ * qe * qf * Ve * Vf + SqChiZ * xe * xf; |
| 118 | ff1 = +2 * ReChiZ * qe * qf * Ae * Af + SqChiZ * ye * yf; |
| 119 | Born = (1.0 + costhe * costhe) * ff0 + 2.0 * costhe * ff1; |
| 120 | // Colour factor |
| 121 | Born = NCf * Born; |
| 122 | // Crude method of correcting threshold, cos(theta) depencence incorrect!!! |
| 123 | if (svar < 4.0 * amfin * amfin) { |
| 124 | thresh = 0.0; |
| 125 | } else if (svar < 16.0 * amfin * amfin) { |
| 126 | amx2 = 4.0 * amfin * amfin / svar; |
| 127 | thresh = sqrt(1.0 - amx2) * (1.0 + amx2 / 2.0); |
| 128 | } else { |
| 129 | thresh = 1.0; |
| 130 | } |
| 131 | |
| 132 | Born = Born * thresh; |
| 133 | return Born; |
| 134 | } |
| 135 | |
| 136 | |
| 137 | // ---------------------------------------------------------------------- |
| 138 | // THIS ROUTINE CALCULATES BORN ASYMMETRY. |
| 139 | // IT EXPLOITS THE FACT THAT BORN X. SECTION = A + B*C + D*C**2 |
| 140 | // |
| 141 | // called by : EVENTM |
| 142 | // ---------------------------------------------------------------------- |
| 143 | // |
| 144 | double PhotosMEforZ::AFBCALC(double SVAR, int IDEE, int IDFF) |
| 145 | { |
| 146 | int KOLOR, KOLOR1; |
| 147 | double T3e, qe, T3f, qf, A, B; |
| 148 | GIVIZO(IDEE, -1, &T3e, &qe, &KOLOR); |
| 149 | GIVIZO(IDFF, -1, &T3f, &qf, &KOLOR1); |
| 150 | |
| 151 | A = PHBORNM(SVAR, 0.5, T3e, qe, T3f, qf, KOLOR * KOLOR1); |
| 152 | B = PHBORNM(SVAR, -0.5, T3e, qe, T3f, qf, KOLOR * KOLOR1); |
| 153 | return (A - B) / (A + B) * 5.0 / 2.0 * 3.0 / 8.0; |
| 154 | } |
| 155 | |
| 156 | |
| 157 | int PhotosMEforZ::GETIDEE(int IDE) |
| 158 | { |
| 159 | |
| 160 | int IDEE; |
| 161 | IDEE = -555; |
| 162 | if ((IDE == 11) || (IDE == 13) || (IDE == 15)) { |
| 163 | IDEE = 2; |
| 164 | } else if ((IDE == -11) || (IDE == -13) || (IDE == -15)) { |
| 165 | IDEE = -2; |
| 166 | } else if ((IDE == 12) || (IDE == 14) || (IDE == 16)) { |
| 167 | IDEE = 1; |
| 168 | } else if ((IDE == -12) || (IDE == -14) || (IDE == -16)) { |
| 169 | IDEE = -1; |
| 170 | } else if ((IDE == 1) || (IDE == 3) || (IDE == 5)) { |
| 171 | IDEE = 4; |
| 172 | } else if ((IDE == -1) || (IDE == -3) || (IDE == -5)) { |
| 173 | IDEE = -4; |
| 174 | } else if ((IDE == 2) || (IDE == 4) || (IDE == 6)) { |
| 175 | IDEE = 3; |
| 176 | } else if ((IDE == - 2) || (IDE == -4) || (IDE == -6)) { |
| 177 | IDEE = -3; |
| 178 | } |
| 179 | if (IDEE == -555) {cout << " ERROR IN GETIDEE of PHOTS Z-ME: I3= &4i" << IDEE << endl;} |
| 180 | return IDEE; |
| 181 | } |
| 182 | |
| 183 | |
| 184 | |
| 185 | |
| 186 | //---------------------------------------------------------------------- |
| 187 | // |
| 188 | // PHASYZ: PHotosASYmmetry of Z |
| 189 | // |
| 190 | // Purpose: Calculates born level asymmetry for Z |
| 191 | // between distributions (1-C)**2 and (1+C)**2 |
| 192 | // At present dummy, requrires effective Z and gamma |
| 193 | // Couplings and also spin polarization states |
| 194 | // For initial and final states. |
| 195 | // To be correct this function need to be tuned |
| 196 | // to host generator. Axes orientation polarisation |
| 197 | // conventions etc etc. |
| 198 | // Modularity of PHOTOS would break. |
| 199 | // |
| 200 | // Input Parameters: SVAR |
| 201 | // |
| 202 | // Output Parameters: Function value |
| 203 | // |
| 204 | // Author(s): Z. Was Created at: 10/12/05 |
| 205 | // Last Update: 19/06/13 |
| 206 | // |
| 207 | //---------------------------------------------------------------------- |
| 208 | double PhotosMEforZ::PHASYZ(double SVAR, int IDE, int IDF) |
| 209 | { |
| 210 | |
| 211 | double AFB; |
| 212 | int IDEE, IDFF; |
| 213 | |
| 214 | IDEE = abs(GETIDEE(IDE)); |
| 215 | IDFF = abs(GETIDEE(IDF)); |
| 216 | AFB = -AFBCALC(SVAR, IDEE, IDFF); |
| 217 | // AFB=0 |
| 218 | return 4.0 / 3.0 * AFB; |
| 219 | // write(*,*) 'IDE=',IDE,' IDF=',IDF,' SVAR=',SVAR,'AFB=',AFB |
| 220 | } |
| 221 | |
| 222 | //---------------------------------------------------------------------- |
| 223 | // |
| 224 | // PHWTNLO: PHotosWTatNLO |
| 225 | // |
| 226 | // Purpose: calculates instead of interference weight |
| 227 | // complete NLO weight for vector boson decays |
| 228 | // of pure vector (or pseudovector) couplings |
| 229 | // Proper orientation of beams required. |
| 230 | // This is not standard in PHOTOS. |
| 231 | // At NLO more tuning than in standard is needed. |
| 232 | // |
| 233 | // |
| 234 | // |
| 235 | // Input Parameters: as in function declaration |
| 236 | // |
| 237 | // Output Parameters: Function value |
| 238 | // |
| 239 | // Author(s): Z. Was Created at: 08/12/05 |
| 240 | // Last Update: 20/06/13 |
| 241 | // |
| 242 | //---------------------------------------------------------------------- |
| 243 | double PhotosMEforZ::Zphwtnlo(double svar, double xk, int IDHEP3, int IREP, double qp[4], double qm[4], double ph[4], double pp[4], |
| 244 | double pm[4], double COSTHG, double BETA, double th1, int IDE, int IDF) |
| 245 | { |
| 246 | double C, s, xkaM, xkaP, t, u, t1, u1, BT, BU; |
| 247 | double waga, wagan2; |
| 248 | static int i = 1; |
| 249 | int IBREM; |
| 250 | |
| 251 | |
| 252 | // IBREM is spurious but it numbers branches of MUSTRAAL |
| 253 | IBREM = 1; |
| 254 | if (IREP == 1) IBREM = -1; |
| 255 | |
| 256 | // we calculate C and S, note that TH1 exists in MUSTRAAL as well. |
| 257 | |
| 258 | C = cos(th1); // this parameter is calculated outside of the class |
| 259 | |
| 260 | // from off line application we had: |
| 261 | if (IBREM == -1) C = -C; |
| 262 | // ... we may want to re-check it. |
| 263 | s = sqrt(1.0 - C * C); |
| 264 | |
| 265 | if (IBREM == 1) { |
| 266 | xkaM = (qp[4 - i] * ph[4 - i] - qp[3 - i] * ph[3 - i] - qp[2 - i] * ph[2 - i] - qp[1 - i] * ph[1 - i]) / xk; |
| 267 | xkaP = (qm[4 - i] * ph[4 - i] - qm[3 - i] * ph[3 - i] - qm[2 - i] * ph[2 - i] - qm[1 - i] * ph[1 - i]) / xk; |
| 268 | } else { |
| 269 | xkaP = (qp[4 - i] * ph[4 - i] - qp[3 - i] * ph[3 - i] - qp[2 - i] * ph[2 - i] - qp[1 - i] * ph[1 - i]) / xk; |
| 270 | xkaM = (qm[4 - i] * ph[4 - i] - qm[3 - i] * ph[3 - i] - qm[2 - i] * ph[2 - i] - qm[1 - i] * ph[1 - i]) / xk; |
| 271 | } |
| 272 | |
| 273 | // XK=2*PHEP(4,nhep)/PHEP(4,1)/xphmax ! it is not used becuse here |
| 274 | // ! order of emissions is meaningless |
| 275 | // |
| 276 | // DELTA=2*PHEP(5,4)**2/svar/(1+(1-XK)**2)*(xKAP/xKAM+xKAM/xKAP) |
| 277 | // waga=svar/4./xkap |
| 278 | // waga=waga*(1.D0-COSTHG*BETA) ! sprawdzone 1= svar/xKAp/4 * (1.D0-COSTHG*BETA) |
| 279 | // waga=waga*(1-delta) /wt2 ! sprawdzone ze to jest =2/(1.D0+COSTHG*BETA) |
| 280 | // ! czyli ubija de-interferencje |
| 281 | |
| 282 | |
| 283 | // this is true only for intermediate resonances with afb=0! |
| 284 | t = 2 * (qp[4 - i] * pp[4 - i] - qp[3 - i] * pp[3 - i] - qp[2 - i] * pp[2 - i] - qp[1 - i] * pp[1 - i]); |
| 285 | u = 2 * (qm[4 - i] * pp[4 - i] - qm[3 - i] * pp[3 - i] - qm[2 - i] * pp[2 - i] - qm[1 - i] * pp[1 - i]); |
| 286 | u1 = 2 * (qp[4 - i] * pm[4 - i] - qp[3 - i] * pm[3 - i] - qp[2 - i] * pm[2 - i] - qp[1 - i] * pm[1 - i]); |
| 287 | t1 = 2 * (qm[4 - i] * pm[4 - i] - qm[3 - i] * pm[3 - i] - qm[2 - i] * pm[2 - i] - qm[1 - i] * pm[1 - i]); |
| 288 | |
| 289 | // basically irrelevant lines ... |
| 290 | t = t - (qp[4 - i] * qp[4 - i] - qp[3 - i] * qp[3 - i] - qp[2 - i] * qp[2 - i] - qp[1 - i] * qp[1 - i]); |
| 291 | u = u - (qm[4 - i] * qm[4 - i] - qm[3 - i] * qm[3 - i] - qm[2 - i] * qm[2 - i] - qm[1 - i] * qm[1 - i]); |
| 292 | u1 = u1 - (qp[4 - i] * qp[4 - i] - qp[3 - i] * qp[3 - i] - qp[2 - i] * qp[2 - i] - qp[1 - i] * qp[1 - i]); |
| 293 | t1 = t1 - (qm[4 - i] * qm[4 - i] - qm[3 - i] * qm[3 - i] - qm[2 - i] * qm[2 - i] - qm[1 - i] * qm[1 - i]); |
| 294 | |
| 295 | // we adjust to what is f-st,s-nd beam flavour |
| 296 | if (IDE * IDHEP3 > 0) { |
| 297 | BT = 1.0 + PHASYZ(svar, IDE, IDF); |
| 298 | BU = 1.0 - PHASYZ(svar, IDE, IDF); |
| 299 | } else { |
| 300 | BT = 1.0 - PHASYZ(svar, IDE, IDF); |
| 301 | BU = 1.0 + PHASYZ(svar, IDE, IDF); |
| 302 | } |
| 303 | wagan2 = 2 * (BT * t * t + BU * u * u + BT * t1 * t1 + BU * u1 * u1) |
| 304 | / (1 + (1 - xk) * (1 - xk)) * 2.0 / (BT * (1 - C) * (1 - C) + BU * (1 + C) * (1 + C)) / svar / svar; |
| 305 | |
| 306 | wagan2 = wagan2 * VakPol(qp, qm, ph, pp, pm, IDE, IDF); // default VakWpol returns 1. Hook for the user function |
| 307 | //! waga=waga*wagan2 |
| 308 | //! waga=waga*(1-delta) /wt2 ! sprawdzone ze to jest =2/(1.D0+COSTHG*BETA) |
| 309 | waga = 2 / (1.0 + COSTHG * BETA) * wagan2; |
| 310 | //! % * svar/4./xkap*(1.D0-COSTHG*BETA)*sqrt(1.0-xk) |
| 311 | |
| 312 | |
| 313 | if (wagan2 <= 3.8) return waga; |
| 314 | |
| 315 | // |
| 316 | // exceptional case wagan2>3.8 |
| 317 | // it should correspond to extremely high bremssthahlung in multiphot conf. |
| 318 | // |
| 319 | FILE* PHLUN = stdoutstdout; |
| 320 | |
| 321 | |
| 322 | // fprintf(PHLUN,"") 'phwtnlo= ',phwtnlo |
| 323 | // fprintf(PHLUN,"") 'idhepy= ',IDHEP[1-i],IDHEP[2-i],IDHEP[3-i],IDHEP[4-i],IDHEP(5) |
| 324 | fprintf(PHLUN, " IDE= %i IDF= %i", IDE, IDF); |
| 325 | fprintf(PHLUN, "bt,bu,bt+bu= %f %f %f", BT, BU, BT + BU); |
| 326 | PHODMP(); // we will activate this once PHODMP(); is re-written |
| 327 | |
| 328 | fprintf(PHLUN, " "); |
| 329 | fprintf(PHLUN, "%i %i <-- IREP,IBREM", IREP, IBREM); |
| 330 | //! fprintf(PHLUN,"%f %f %f %f %f") 'pneutr= ',phomom.pneutr[0],phomom.pneutr[1],phomom.pneutr[2],phomom.pneutr[3],phomom.pneutr[4]; |
| 331 | fprintf(PHLUN, "%f %f %f %f qp = ", qp[0], qp[1], qp[2], qp[3]); |
| 332 | fprintf(PHLUN, "%f %f %f %f qm = ", qm[0], qm[1], qm[2], qm[3]); |
| 333 | fprintf(PHLUN, " "); |
| 334 | fprintf(PHLUN, "%f %f %f %f ph = ", ph[0], ph[1], ph[2], ph[3]); |
| 335 | // fprintf(PHLUN,"") 'p1= ',PHEP(1,1),PHEP(2,1),PHEP(3,1),PHEP(4,1) |
| 336 | // fprintf(PHLUN,"") 'p2= ',PHEP(1,2),PHEP(2,2),PHEP(3,2),PHEP(4,2) |
| 337 | // fprintf(PHLUN,"") 'p3= ',PHEP(1,3),PHEP(2,3),PHEP(3,3),PHEP(4,3) |
| 338 | // fprintf(PHLUN,"") 'p4= ',PHEP(1,4),PHEP(2,4),PHEP(3,4),PHEP(4,4) |
| 339 | // fprintf(PHLUN,"") 'p5= ',PHEP(1,5),PHEP(2,5),PHEP(3,5),PHEP(4,5) |
| 340 | |
| 341 | fprintf(PHLUN, " c= %f theta= %f", C, th1); |
| 342 | // fprintf(PHLUN,"") 'photos waga daje ... IBREM=',IBREM,' waga=',waga |
| 343 | // fprintf(PHLUN,"") 'xk,COSTHG,c',xk,COSTHG,c |
| 344 | // fprintf(PHLUN,"") svar/4./xkap*(1.D0-COSTHG*BETA), |
| 345 | // $ (1-delta) /wt2 *(1.D0+COSTHG*BETA)/2, wagan2 |
| 346 | // fprintf(PHLUN,"") ' delta, wt2,beta', delta, wt2,beta |
| 347 | fprintf(PHLUN, " - "); |
| 348 | fprintf(PHLUN, "t,u = %f %f", t, u); |
| 349 | fprintf(PHLUN, "t1,u1 = %f %f", t1, u1); |
| 350 | fprintf(PHLUN, "sredniaki = %f %f", svar * (1 - C) / 2, svar * (1 + C) / 2); |
| 351 | // ! fprintf(PHLUN,"") 'xk= %f c= %f COSTHG= %f' ,xk,c,COSTHG |
| 352 | fprintf(PHLUN, "PHASYZ(svar)=',%f,' svar= %f',' waga= %f", PHASYZ(svar, IDE, IDF), svar, waga); |
| 353 | fprintf(PHLUN, " - "); |
| 354 | fprintf(PHLUN, "BT-part= %f BU-part= %f", |
| 355 | 2 * (BT * t * t + BT * t1 * t1) |
| 356 | / (1 + (1 - xk) * (1 - xk)) * 2.0 / (BT * (1 - C) * (1 - C)) / svar / svar, |
| 357 | 2 * (BU * u * u + BU * u1 * u1) |
| 358 | / (1 + (1 - xk) * (1 - xk)) * 2.0 / (BU * (1 + C) * (1 + C)) / svar / svar); |
| 359 | fprintf(PHLUN, "BT-part*BU-part= %f wagan2= %f", |
| 360 | 2 * (BT * t * t + BT * t1 * t1) |
| 361 | / (1 + (1 - xk) * (1 - xk)) * 2.0 / (BT * (1 - C) * (1 - C)) / svar / svar |
| 362 | * 2 * (BU * u * u + BU * u1 * u1) |
| 363 | / (1 + (1 - xk) * (1 - xk)) * 2.0 / (BU * (1 + C) * (1 + C)) / svar / svar, wagan2); |
| 364 | |
| 365 | fprintf(PHLUN, "wagan2= %f", wagan2); |
| 366 | fprintf(PHLUN, " ################### "); |
| 367 | |
| 368 | // wagan2=wagan2*VakPol(qp,qm,ph,pp,pm,IDE,IDF); // default VakWpol returns 1. Hook for the user function |
| 369 | wagan2 = 3.8; // ! overwrite |
| 370 | waga = 2 / (1.0 + COSTHG * BETA) * wagan2 ; |
| 371 | // % * svar/4./xkap*(1.D0-COSTHG*BETA)*sqrt(1.0-xk) |
| 372 | |
| 373 | |
| 374 | |
| 375 | |
| 376 | return waga; |
| 377 | |
| 378 | } |
| 379 | |
| 380 | double PhotosMEforZ::VakPol(double qp[4], double qm[4], double ph[4], double pp[4], double pm[4], int IDE, int IDF) |
| 381 | { |
| 382 | return PhotosMEforZ::currentVakPol(qp, qm, ph, pp, pm, IDE, IDF); |
| 383 | } |
| 384 | |
| 385 | double PhotosMEforZ::default_VakPol(double qp[4], double qm[4], double ph[4], double pp[4], double pm[4], int IDE, int IDF) |
| 386 | { |
| 387 | return 1; // that is default. |
| 388 | } |
| 389 | |
| 390 | void PhotosMEforZ::set_VakPol(double (*fun)(double[4], double[4], double[4], double[4], double[4], int, int)) |
| 391 | { |
| 392 | if (fun == NULL__null) PhotosMEforZ::currentVakPol = default_VakPol; |
| 393 | else PhotosMEforZ::currentVakPol = fun; |
| 394 | } |
| 395 | |
| 396 | //---------------------------------------------------------------------- |
| 397 | // |
| 398 | // PHWTNLO: PHotosWTatNLO |
| 399 | // |
| 400 | // Purpose: calculates instead of interference weight |
| 401 | // complete NLO weight for vector boson decays |
| 402 | // of pure vector (or pseudovector) couplings |
| 403 | // Proper orientation of beams required. |
| 404 | // Uses Zphwtnlo encapsulating actual matrix element |
| 405 | // At NLO more tuning on kinematical conf. |
| 406 | // than in standard is needed. |
| 407 | // Works with KORALZ and KKM// |
| 408 | // Note some commented out commons from MUSTAAL, KORALZ |
| 409 | // |
| 410 | // Input Parameters: Common /PHOEVT/ /PHOPS/ /PHOREST/ /PHOPRO/ |
| 411 | // |
| 412 | // Output Parameters: Function value |
| 413 | // |
| 414 | // Author(s): Z. Was Created at: 08/12/05 |
| 415 | // Last Update: 23/06/13 |
| 416 | // |
| 417 | //---------------------------------------------------------------------- |
| 418 | |
| 419 | double PhotosMEforZ::phwtnlo() |
| 420 | { |
| 421 | // fi3 orientation of photon, fi1,th1 orientation of neutral |
| 422 | |
| 423 | // COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG |
| 424 | |
| 425 | // COMMON /PHOREST/ FI3,fi1,th1 |
| 426 | // COMMON /PHWT/ BETA,WT1,WT2,WT3 |
| 427 | // COMMON/PHOPRO/PROBH,CORWT,XF,IREP |
| 428 | // COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5) |
| 429 | // static double PI=3.141592653589793238462643; |
| 430 | static int i = 1; |
| 431 | int K, L, IDHEP3, IDUM = 0; |
| 432 | int IDE, IDF; |
| 433 | double QP[4], QM[4], PH[4], QQ[4], PP[4], PM[4], QQS[4]; |
| 434 | double XK, ENE, svar; |
| 435 | |
| 436 | // REAL*8 s,c,svar,xkaM,xkaP,xk,phwtnlo,xdumm,PHINT |
| 437 | // REAL*8 ENE,a,t,u,t1,u1,wagan2,waga,PHASYZ,BT,BU,ENEB |
| 438 | // INTEGER IBREM,K,L,IREP,IDUM,IDHEP3 |
| 439 | // integer icont,ide,idf |
| 440 | // REAL*8 delta |
| 441 | |
| 442 | ///////////////////// |
| 443 | // phlupa(299500); |
| 444 | |
| 445 | |
| 446 | ///////////////////// |
| 447 | // phlupa(299500); |
| 448 | |
| 449 | XK = 2.0 * pho.phep[pho.nhep - i][4 - i] / pho.phep[1 - i][4 - i]; |
| 450 | |
| 451 | // XK=2.0*pho.phep[pho.nhep-i][4-i]/pho.phep[1-i][4-i]/phophs.xphmax; // it is not used becuse here |
| 452 | //order of emissions is meaningless |
| 453 | if (pho.nhep <= 4) XK = 0.0; |
| 454 | // the mother must be Z or gamma* !!!! |
| 455 | |
| 456 | if (XK > 1.0e-10 && (pho.idhep[1 - i] == 22 || pho.idhep[1 - i] == 23)) { |
| 457 | |
| 458 | // write(*,*) 'nhep=',nhep |
| 459 | // DO K=1,3 ENDDO |
| 460 | // IF (K.EQ.1) IBREM= 1 |
| 461 | // IF (K.EQ.2) IBREM=-1 |
| 462 | // ICONT=ICONT+1 |
| 463 | // IBREM=IBREX ! that will be input parameter. |
| 464 | // IBREM=IBREY ! that IS now input parameter. |
| 465 | |
| 466 | // We initialize twice 4-vectors, here and again later after boost |
| 467 | // must be the same way. Important is how the reduction procedure will work. |
| 468 | // It seems at present that the beams must be translated to be back to back. |
| 469 | // this may be done after initialising, thus on 4-vectors. |
| 470 | |
| 471 | for (K = 1; K < 5; K++) { |
| 472 | PP[K - i] = pho.phep[1 - i][K - i]; |
| 473 | PM[K - i] = pho.phep[2 - i][K - i]; |
| 474 | QP[K - i] = pho.phep[3 - i][K - i]; |
| 475 | QM[K - i] = pho.phep[4 - i][K - i]; |
| 476 | PH[K - i] = pho.phep[pho.nhep - i][K - i]; |
| 477 | QQ[K - i] = 0.0; |
| 478 | QQS[K - i] = QP[K - i] + QM[K - i]; |
| 479 | } |
| 480 | |
| 481 | |
| 482 | PP[4 - i] = (pho.phep[1 - i][4 - i] + pho.phep[2 - i][4 - i]) / 2.0; |
| 483 | PM[4 - i] = (pho.phep[1 - i][4 - i] + pho.phep[2 - i][4 - i]) / 2.0; |
| 484 | PP[3 - i] = PP[4 - i]; |
| 485 | PM[3 - i] = -PP[4 - i]; |
| 486 | |
| 487 | for (L = 5; L <= pho.nhep - 1; L++) { |
| 488 | for (K = 1; K < 5; K++) { |
| 489 | QQ [K - i] = QQ [K - i] + pho.phep[L - i][K - i]; |
| 490 | QQS[K - i] = QQS[K - i] + pho.phep[L - i][K - i]; |
| 491 | } |
| 492 | } |
| 493 | |
| 494 | // go to the restframe of 3 |
| 495 | PHOB(1, QQS, QP); |
| 496 | PHOB(1, QQS, QM); |
| 497 | PHOB(1, QQS, QQ); |
| 498 | ENE = (QP[4 - i] + QM[4 - i] + QQ[4 - i]) / 2; |
| 499 | |
| 500 | // preserve direction of emitting particle and wipeout QQ |
| 501 | if (phopro.irep == 1) { |
| 502 | double a = sqrt(ENE * ENE - pho.phep[3 - i][5 - i] * pho.phep[3 - i][5 - i]) / sqrt(QM[4 - i] * QM[4 - i] - pho.phep[3 - i][5 - i] |
| 503 | * pho.phep[3 - i][5 - i]); |
| 504 | QM[1 - i] = QM[1 - i] * a; |
| 505 | QM[2 - i] = QM[2 - i] * a; |
| 506 | QM[3 - i] = QM[3 - i] * a; |
| 507 | QP[1 - i] = -QM[1 - i]; |
| 508 | QP[2 - i] = -QM[2 - i]; |
| 509 | QP[3 - i] = -QM[3 - i]; |
| 510 | } else { |
| 511 | double a = sqrt(ENE * ENE - pho.phep[3 - i][5 - i] * pho.phep[3 - i][5 - i]) / sqrt(QP[4 - i] * QP[4 - i] - pho.phep[3 - i][5 - i] |
| 512 | * pho.phep[3 - i][5 - i]); |
| 513 | QP[1 - i] = QP[1 - i] * a; |
| 514 | QP[2 - i] = QP[2 - i] * a; |
| 515 | QP[3 - i] = QP[3 - i] * a; |
| 516 | QM[1 - i] = -QP[1 - i]; |
| 517 | QM[2 - i] = -QP[2 - i]; |
| 518 | QM[3 - i] = -QP[3 - i]; |
| 519 | } |
| 520 | QP[4 - i] = ENE; |
| 521 | QM[4 - i] = ENE; |
| 522 | // go back to reaction frame (QQ eliminated) |
| 523 | PHOB(-1, QQS, QP); |
| 524 | PHOB(-1, QQS, QM); |
| 525 | PHOB(-1, QQS, QQ); |
| 526 | |
| 527 | svar = pho.phep[1 - i][4 - i] * pho.phep[1 - i][4 - i]; |
| 528 | |
| 529 | IDE = hep.idhep[1 - i]; |
| 530 | IDF = hep.idhep[4 - i]; |
| 531 | if (abs(hep.idhep[4 - i]) == abs(hep.idhep[3 - i])) IDF = hep.idhep[3 - i]; |
| 532 | |
| 533 | IDHEP3 = pho.idhep[3 - i]; |
| 534 | return Zphwtnlo(svar, XK, IDHEP3, phopro.irep, QP, QM, PH, PP, PM, phophs.costhg, phwt.beta, phorest.th1, IDE, IDF); |
| 535 | } else { |
| 536 | // in other cases we just use default setups. |
| 537 | return PHINT(IDUM); |
| 538 | } |
| 539 | } |
| 540 | |
| 541 | } // namespace Photospp |
| 542 |