| 3820ca8e |
1 | |
| 2 | CDECK ID>, HWHBGF. |
| 3 | |
| 4 | *CMZ :- -26/04/91 11.11.55 by Bryan Webber |
| 5 | |
| 6 | *-- Author : Giovanni Abbiendi & Luca Stanco |
| 7 | |
| 8 | C----------------------------------------------------------------------- |
| 9 | |
| 10 | SUBROUTINE HWHBGF |
| 11 | |
| 12 | C----------------------------------------------------------------------- |
| 13 | |
| 14 | C Order Alpha_s processes in charged lepton-hadron collisions |
| 15 | |
| 16 | C |
| 17 | |
| 18 | C Process code IPROC has to be set in the Main Program |
| 19 | |
| 20 | C the following codes IPROC may be selected |
| 21 | |
| 22 | C |
| 23 | |
| 24 | C 9100 : NC BOSON-GLUON FUSION |
| 25 | |
| 26 | C 9100+IQK (IQK=1,...,6) : produced flavour is IQK |
| 27 | |
| 28 | C 9107 : produced J/psi + gluon |
| 29 | |
| 30 | C |
| 31 | |
| 32 | C 9110 : NC QCD COMPTON |
| 33 | |
| 34 | C 9110+IQK (IQK=1,...,12) : struck parton is IQK |
| 35 | |
| 36 | C |
| 37 | |
| 38 | C 9130 : NC order alpha_s processes (9100+9110) |
| 39 | |
| 40 | C |
| 41 | |
| 42 | C Select maximum and minimum generated flavour when IQK=0 |
| 43 | |
| 44 | C setting IFLMIN and IFLMAX in the Main Program |
| 45 | |
| 46 | C (allowed values from 1 to 6), default are 1 and 5 |
| 47 | |
| 48 | C allowing d,u,s,c,b,dbar,ubar,sbar,cbar,bbar |
| 49 | |
| 50 | C |
| 51 | |
| 52 | C CHARGED CURRENT Boson-Gluon Fusion processes |
| 53 | |
| 54 | C 9141 : CC s cbar (c sbar) |
| 55 | |
| 56 | C 9142 : CC b cbar (c bbar) |
| 57 | |
| 58 | C 9143 : CC s tbar (t cbar) |
| 59 | |
| 60 | C 9144 : CC b tbar (t bbar) |
| 61 | |
| 62 | C |
| 63 | |
| 64 | C other inputs : Q2MIN,Q2MAX,YBMIN,YBMAX,PTMIN,EMMIN,EMMAX |
| 65 | |
| 66 | C when IPROC=(1)9107 : as above but Q2WWMN, Q2WWMX substitute |
| 67 | |
| 68 | C Q2MIN and Q2MAX (EPA is used); ZJMAX cut |
| 69 | |
| 70 | C |
| 71 | |
| 72 | C Add 10000 to suppress soft remnant fragmentation |
| 73 | |
| 74 | C |
| 75 | |
| 76 | C Mean EVWGT = cross section in nanoBarn |
| 77 | |
| 78 | C |
| 79 | |
| 80 | C----------------------------------------------------------------------- |
| 81 | |
| 82 | INCLUDE 'HERWIG61.INC' |
| 83 | |
| 84 | DOUBLE PRECISION HWR,LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP, |
| 85 | |
| 86 | & ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT,FSIGMA(18), |
| 87 | |
| 88 | & SIGSUM,PROB,PRAN,PVRT(4),X |
| 89 | |
| 90 | INTEGER IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,IPROO,LEPFIN,ID1,ID2,I,IDD |
| 91 | |
| 92 | LOGICAL CHARGD,INCLUD(18),INSIDE(18) |
| 93 | |
| 94 | EXTERNAL HWR |
| 95 | |
| 96 | SAVE LEPFIN,ID1,ID2,FSIGMA,SIGSUM |
| 97 | |
| 98 | COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF, |
| 99 | |
| 100 | & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL, |
| 101 | |
| 102 | & IPROO,CHARGD,INCLUD,INSIDE |
| 103 | |
| 104 | C---Initialization |
| 105 | |
| 106 | IF (FSTWGT) THEN |
| 107 | |
| 108 | C---LEP = 1 FOR LEPTONS, -1 FOR ANTILEPTONS |
| 109 | |
| 110 | LEP=ZERO |
| 111 | |
| 112 | IF (IDHW(1).GE.121.AND.IDHW(1).LE.126) THEN |
| 113 | |
| 114 | LEP=ONE |
| 115 | |
| 116 | ELSEIF (IDHW(1).GE.127.AND.IDHW(1).LE.132) THEN |
| 117 | |
| 118 | LEP=-ONE |
| 119 | |
| 120 | ENDIF |
| 121 | |
| 122 | IF (LEP.EQ.ZERO) CALL HWWARN('HWHBGF',500,*999) |
| 123 | |
| 124 | IPROO=MOD(IPROC,100)/10 |
| 125 | |
| 126 | IF (IPROO.EQ.0.OR.IPROO.EQ.4) THEN |
| 127 | |
| 128 | IQK=MOD(IPROC,10) |
| 129 | |
| 130 | IFL=IQK |
| 131 | |
| 132 | IF (IQK.EQ.7) IFL=164 |
| 133 | |
| 134 | CHARGD=IPROO.EQ.4 |
| 135 | |
| 136 | ELSEIF (IPROO.EQ.1.OR.IPROO.EQ.2) THEN |
| 137 | |
| 138 | IQK=MOD(IPROC,100)-10 |
| 139 | |
| 140 | IFL=IQK+6 |
| 141 | |
| 142 | CHARGD=.FALSE. |
| 143 | |
| 144 | ELSEIF (IPROO.EQ.3) THEN |
| 145 | |
| 146 | IQK=0 |
| 147 | |
| 148 | IFL=0 |
| 149 | |
| 150 | CHARGD=.FALSE. |
| 151 | |
| 152 | ELSE |
| 153 | |
| 154 | CALL HWWARN('HWHBGF',501,*999) |
| 155 | |
| 156 | ENDIF |
| 157 | |
| 158 | C |
| 159 | |
| 160 | LEPFIN = IDHW(1) |
| 161 | |
| 162 | IF(CHARGD) THEN |
| 163 | |
| 164 | LEPFIN = IDHW(1)+1 |
| 165 | |
| 166 | IF (IQK.EQ.1) THEN |
| 167 | |
| 168 | IFLAVU=4 |
| 169 | |
| 170 | IFLAVD=3 |
| 171 | |
| 172 | ID1 = 3 |
| 173 | |
| 174 | ID2 = 10 |
| 175 | |
| 176 | ELSEIF (IQK.EQ.2) THEN |
| 177 | |
| 178 | IFLAVU=4 |
| 179 | |
| 180 | IFLAVD=5 |
| 181 | |
| 182 | ID1 = 5 |
| 183 | |
| 184 | ID2 = 10 |
| 185 | |
| 186 | ELSEIF (IQK.EQ.3) THEN |
| 187 | |
| 188 | IFLAVU=6 |
| 189 | |
| 190 | IFLAVD=3 |
| 191 | |
| 192 | ID1 = 3 |
| 193 | |
| 194 | ID2 =12 |
| 195 | |
| 196 | ELSE |
| 197 | |
| 198 | IFLAVU=6 |
| 199 | |
| 200 | IFLAVD=5 |
| 201 | |
| 202 | ID1 = 5 |
| 203 | |
| 204 | ID2 =12 |
| 205 | |
| 206 | ENDIF |
| 207 | |
| 208 | IF (LEP.EQ.-1.) THEN |
| 209 | |
| 210 | IDD=ID1 |
| 211 | |
| 212 | ID1=ID2-6 |
| 213 | |
| 214 | ID2=IDD+6 |
| 215 | |
| 216 | ENDIF |
| 217 | |
| 218 | ENDIF |
| 219 | |
| 220 | C |
| 221 | |
| 222 | IF (IQK.EQ.0) THEN |
| 223 | |
| 224 | DO I=1,18 |
| 225 | |
| 226 | INCLUD(I)=.TRUE. |
| 227 | |
| 228 | ENDDO |
| 229 | |
| 230 | IMIN=1 |
| 231 | |
| 232 | IMAX=18 |
| 233 | |
| 234 | DO I=1,6 |
| 235 | |
| 236 | IF (I.LT.IFLMIN.OR.I.GT.IFLMAX) INCLUD(I)=.FALSE. |
| 237 | |
| 238 | ENDDO |
| 239 | |
| 240 | DO I=7,18 |
| 241 | |
| 242 | IF (I.LE.12) THEN |
| 243 | |
| 244 | IF (I-6.LT.IFLMIN.OR.I-6.GT.IFLMAX) INCLUD(I)=.FALSE. |
| 245 | |
| 246 | ELSE |
| 247 | |
| 248 | IF (I-12.LT.IFLMIN.OR.I-12.GT.IFLMAX) INCLUD(I)=.FALSE. |
| 249 | |
| 250 | ENDIF |
| 251 | |
| 252 | ENDDO |
| 253 | |
| 254 | IF (IPROO.EQ.0) THEN |
| 255 | |
| 256 | DO I=7,18 |
| 257 | |
| 258 | INCLUD(I)=.FALSE. |
| 259 | |
| 260 | ENDDO |
| 261 | |
| 262 | IMIN=IFLMIN |
| 263 | |
| 264 | IMAX=IFLMAX |
| 265 | |
| 266 | ELSEIF (IPROO.EQ.1.OR.IPROO.EQ.2) THEN |
| 267 | |
| 268 | DO I=1,6 |
| 269 | |
| 270 | INCLUD(I)=.FALSE. |
| 271 | |
| 272 | ENDDO |
| 273 | |
| 274 | IMIN=IFLMIN+6 |
| 275 | |
| 276 | IMAX=IFLMAX+12 |
| 277 | |
| 278 | ELSEIF (IPROO.EQ.3) THEN |
| 279 | |
| 280 | IMIN=IFLMIN |
| 281 | |
| 282 | IMAX=IFLMAX+12 |
| 283 | |
| 284 | ENDIF |
| 285 | |
| 286 | ELSEIF (IQK.NE.0 .AND. (.NOT.CHARGD)) THEN |
| 287 | |
| 288 | DO I=1,18 |
| 289 | |
| 290 | INCLUD(I)=.FALSE. |
| 291 | |
| 292 | ENDDO |
| 293 | |
| 294 | IF (IFL.LE.18) THEN |
| 295 | |
| 296 | INCLUD(IFL)=.TRUE. |
| 297 | |
| 298 | IMIN=IFL |
| 299 | |
| 300 | IMAX=IFL |
| 301 | |
| 302 | ELSEIF (IFL.EQ.164) THEN |
| 303 | |
| 304 | INCLUD(7)=.TRUE. |
| 305 | |
| 306 | IMIN=7 |
| 307 | |
| 308 | IMAX=7 |
| 309 | |
| 310 | ENDIF |
| 311 | |
| 312 | ENDIF |
| 313 | |
| 314 | ENDIF |
| 315 | |
| 316 | C---End of initialization |
| 317 | |
| 318 | IF(GENEV) THEN |
| 319 | |
| 320 | IF (.NOT.CHARGD) THEN |
| 321 | |
| 322 | IF (IQK.EQ.0) THEN |
| 323 | |
| 324 | PRAN= SIGSUM * HWR() |
| 325 | |
| 326 | PROB=ZERO |
| 327 | |
| 328 | DO 10 IFL=IMIN,IMAX |
| 329 | |
| 330 | IF (.NOT.INSIDE(IFL)) GOTO 10 |
| 331 | |
| 332 | PROB=PROB+FSIGMA(IFL) |
| 333 | |
| 334 | IF (PROB.GE.PRAN) GOTO 20 |
| 335 | |
| 336 | 10 CONTINUE |
| 337 | |
| 338 | ENDIF |
| 339 | |
| 340 | C---at this point the subprocess has been selected (IFL) |
| 341 | |
| 342 | 20 CONTINUE |
| 343 | |
| 344 | IF (IFL.LE.6) THEN |
| 345 | |
| 346 | C---Boson-Gluon Fusion event |
| 347 | |
| 348 | IDHW(NHEP+1)=IDHW(1) |
| 349 | |
| 350 | IDHW(NHEP+2)=13 |
| 351 | |
| 352 | IDHW(NHEP+3)=15 |
| 353 | |
| 354 | IDHW(NHEP+4)=LEPFIN |
| 355 | |
| 356 | IDHW(NHEP+5)=IFL |
| 357 | |
| 358 | IDHW(NHEP+6)=IFL+6 |
| 359 | |
| 360 | ELSEIF (IFL.GE.7.AND.IFL.LE.18) THEN |
| 361 | |
| 362 | C---QCD_Compton event |
| 363 | |
| 364 | IDHW(NHEP+1)=IDHW(1) |
| 365 | |
| 366 | IDHW(NHEP+2)=IFL-6 |
| 367 | |
| 368 | IDHW(NHEP+3)=15 |
| 369 | |
| 370 | IDHW(NHEP+4)=LEPFIN |
| 371 | |
| 372 | IDHW(NHEP+5)=IFL-6 |
| 373 | |
| 374 | IDHW(NHEP+6)=13 |
| 375 | |
| 376 | ELSEIF (IFL.EQ.164) THEN |
| 377 | |
| 378 | C---gamma+gluon-->J/Psi+gluon |
| 379 | |
| 380 | IDHW(NHEP+1)=IDHW(1) |
| 381 | |
| 382 | IDHW(NHEP+2)=13 |
| 383 | |
| 384 | IDHW(NHEP+3)=15 |
| 385 | |
| 386 | IDHW(NHEP+4)=LEPFIN |
| 387 | |
| 388 | IDHW(NHEP+5)=164 |
| 389 | |
| 390 | IDHW(NHEP+6)=13 |
| 391 | |
| 392 | ELSE |
| 393 | |
| 394 | CALL HWWARN('HWHBGF',503,*999) |
| 395 | |
| 396 | ENDIF |
| 397 | |
| 398 | ELSE |
| 399 | |
| 400 | C---Charged current event of specified flavours |
| 401 | |
| 402 | IDHW(NHEP+1)=IDHW(1) |
| 403 | |
| 404 | IDHW(NHEP+2)=13 |
| 405 | |
| 406 | IDHW(NHEP+3)=15 |
| 407 | |
| 408 | IDHW(NHEP+4)=LEPFIN |
| 409 | |
| 410 | IDHW(NHEP+5)=ID1 |
| 411 | |
| 412 | IDHW(NHEP+6)=ID2 |
| 413 | |
| 414 | ENDIF |
| 415 | |
| 416 | C |
| 417 | |
| 418 | DO 1 I=NHEP+1,NHEP+6 |
| 419 | |
| 420 | 1 IDHEP(I)=IDPDG(IDHW(I)) |
| 421 | |
| 422 | C |
| 423 | |
| 424 | C---Codes common for all processes |
| 425 | |
| 426 | ISTHEP(NHEP+1)=111 |
| 427 | |
| 428 | ISTHEP(NHEP+2)=112 |
| 429 | |
| 430 | ISTHEP(NHEP+3)=110 |
| 431 | |
| 432 | ISTHEP(NHEP+4)=113 |
| 433 | |
| 434 | ISTHEP(NHEP+5)=114 |
| 435 | |
| 436 | ISTHEP(NHEP+6)=114 |
| 437 | |
| 438 | C |
| 439 | |
| 440 | DO I=NHEP+1,NHEP+6 |
| 441 | |
| 442 | JMOHEP(1,I)=NHEP+3 |
| 443 | |
| 444 | JDAHEP(1,I)=0 |
| 445 | |
| 446 | ENDDO |
| 447 | |
| 448 | C---Incoming lepton |
| 449 | |
| 450 | JMOHEP(2,NHEP+1)=NHEP+4 |
| 451 | |
| 452 | JDAHEP(2,NHEP+1)=NHEP+4 |
| 453 | |
| 454 | C---Hard Process C.M. |
| 455 | |
| 456 | JMOHEP(1,NHEP+3)=NHEP+1 |
| 457 | |
| 458 | JMOHEP(2,NHEP+3)=NHEP+2 |
| 459 | |
| 460 | JDAHEP(1,NHEP+3)=NHEP+4 |
| 461 | |
| 462 | JDAHEP(2,NHEP+3)=NHEP+6 |
| 463 | |
| 464 | C---Outgoing lepton |
| 465 | |
| 466 | JMOHEP(2,NHEP+4)=NHEP+1 |
| 467 | |
| 468 | JDAHEP(2,NHEP+4)=NHEP+1 |
| 469 | |
| 470 | C |
| 471 | |
| 472 | IF (IFL.LE.6 .OR. CHARGD) THEN |
| 473 | |
| 474 | C---Codes for boson-gluon fusion processes |
| 475 | |
| 476 | C--- Incoming gluon |
| 477 | |
| 478 | JMOHEP(2,NHEP+2)=NHEP+6 |
| 479 | |
| 480 | JDAHEP(2,NHEP+2)=NHEP+5 |
| 481 | |
| 482 | C--- Outgoing quark |
| 483 | |
| 484 | JMOHEP(2,NHEP+5)=NHEP+2 |
| 485 | |
| 486 | JDAHEP(2,NHEP+5)=NHEP+6 |
| 487 | |
| 488 | C--- Outgoing antiquark |
| 489 | |
| 490 | JMOHEP(2,NHEP+6)=NHEP+5 |
| 491 | |
| 492 | JDAHEP(2,NHEP+6)=NHEP+2 |
| 493 | |
| 494 | ELSEIF (IFL.GE.7 .AND. IFL.LE.12) THEN |
| 495 | |
| 496 | C---Codes for V+q --> q+g |
| 497 | |
| 498 | C--- Incoming quark |
| 499 | |
| 500 | JMOHEP(2,NHEP+2)=NHEP+5 |
| 501 | |
| 502 | JDAHEP(2,NHEP+2)=NHEP+6 |
| 503 | |
| 504 | C--- Outgoing quark |
| 505 | |
| 506 | JMOHEP(2,NHEP+5)=NHEP+6 |
| 507 | |
| 508 | JDAHEP(2,NHEP+5)=NHEP+2 |
| 509 | |
| 510 | C--- Outgoing gluon |
| 511 | |
| 512 | JMOHEP(2,NHEP+6)=NHEP+2 |
| 513 | |
| 514 | JDAHEP(2,NHEP+6)=NHEP+5 |
| 515 | |
| 516 | ELSEIF (IFL.GE.13 .AND. IFL.LE.18) THEN |
| 517 | |
| 518 | C---Codes for V+qbar --> qbar+g |
| 519 | |
| 520 | C--- Incoming antiquark |
| 521 | |
| 522 | JMOHEP(2,NHEP+2)=NHEP+6 |
| 523 | |
| 524 | JDAHEP(2,NHEP+2)=NHEP+5 |
| 525 | |
| 526 | C--- Outgoing antiquark |
| 527 | |
| 528 | JMOHEP(2,NHEP+5)=NHEP+2 |
| 529 | |
| 530 | JDAHEP(2,NHEP+5)=NHEP+6 |
| 531 | |
| 532 | C--- Outgoing gluon |
| 533 | |
| 534 | JMOHEP(2,NHEP+6)=NHEP+5 |
| 535 | |
| 536 | JDAHEP(2,NHEP+6)=NHEP+2 |
| 537 | |
| 538 | ELSEIF (IFL.EQ.164) THEN |
| 539 | |
| 540 | C---Codes for Gamma+gluon --> J/Psi+gluon |
| 541 | |
| 542 | C--- Incoming gluon |
| 543 | |
| 544 | JMOHEP(2,NHEP+2)=NHEP+6 |
| 545 | |
| 546 | JDAHEP(2,NHEP+2)=NHEP+6 |
| 547 | |
| 548 | C--- Outgoing J/Psi |
| 549 | |
| 550 | JMOHEP(2,NHEP+5)=NHEP+1 |
| 551 | |
| 552 | JDAHEP(2,NHEP+5)=NHEP+1 |
| 553 | |
| 554 | C--- Outgoing gluon |
| 555 | |
| 556 | JMOHEP(2,NHEP+6)=NHEP+2 |
| 557 | |
| 558 | JDAHEP(2,NHEP+6)=NHEP+2 |
| 559 | |
| 560 | ENDIF |
| 561 | |
| 562 | C---Computation of momenta in Laboratory frame of reference |
| 563 | |
| 564 | CALL HWHBKI |
| 565 | |
| 566 | NHEP=NHEP+6 |
| 567 | |
| 568 | C Decide which quark radiated and assign production vertices |
| 569 | |
| 570 | IF (IFL.LE.6) THEN |
| 571 | |
| 572 | C Boson-Gluon fusion case |
| 573 | |
| 574 | IF (1-Z.LT.HWR()) THEN |
| 575 | |
| 576 | C Gluon splitting to quark |
| 577 | |
| 578 | CALL HWVZRO(4,VHEP(1,NHEP-1)) |
| 579 | |
| 580 | CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP),PVRT) |
| 581 | |
| 582 | CALL HWUDKL(IFL,PVRT,VHEP(1,NHEP)) |
| 583 | |
| 584 | CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-4)) |
| 585 | |
| 586 | ELSE |
| 587 | |
| 588 | C Gluon splitting to antiquark |
| 589 | |
| 590 | CALL HWVZRO(4,VHEP(1,NHEP)) |
| 591 | |
| 592 | CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP-1),PVRT) |
| 593 | |
| 594 | CALL HWUDKL(IFL,PVRT,VHEP(1,NHEP-1)) |
| 595 | |
| 596 | CALL HWVEQU(4,VHEP(1,NHEP-1),VHEP(1,NHEP-4)) |
| 597 | |
| 598 | ENDIF |
| 599 | |
| 600 | ELSEIF (IFL.GE.7.AND.IFL.LE.18) THEN |
| 601 | |
| 602 | C QCD Compton case |
| 603 | |
| 604 | X=1/(1+SHAT/Q2) |
| 605 | |
| 606 | IF (1.LT.HWR()*(1+(1-X-Z)**2+6*X*(1-X)*Z*(1-Z))) THEN |
| 607 | |
| 608 | C Incoming quark radiated the gluon |
| 609 | |
| 610 | CALL HWVZRO(4,VHEP(1,NHEP-1)) |
| 611 | |
| 612 | CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP),PVRT) |
| 613 | |
| 614 | CALL HWUDKL(IFL-6,PVRT,VHEP(1,NHEP)) |
| 615 | |
| 616 | CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-4)) |
| 617 | |
| 618 | ELSE |
| 619 | |
| 620 | C Outgoing quark radiated the gluon |
| 621 | |
| 622 | CALL HWVZRO(4,VHEP(1,NHEP-4)) |
| 623 | |
| 624 | CALL HWVSUM(4,PHEP(1,NHEP-1),PHEP(1,NHEP),PVRT) |
| 625 | |
| 626 | CALL HWUDKL(IFL-6,PVRT,VHEP(1,NHEP)) |
| 627 | |
| 628 | CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-1)) |
| 629 | |
| 630 | ENDIF |
| 631 | |
| 632 | ENDIF |
| 633 | |
| 634 | C---HERWIG gets confused if lepton momentum is different from beam |
| 635 | |
| 636 | C momentum, which it can be if incoming hadron has negative virtuality |
| 637 | |
| 638 | C As a temporary fix, simply copy the momentum. |
| 639 | |
| 640 | C Momentum conservation somehow gets taken care of HWBGEN! |
| 641 | |
| 642 | call hwvequ(5,phep(1,1),phep(1,nhep-5)) |
| 643 | |
| 644 | ELSE |
| 645 | |
| 646 | EVWGT=ZERO |
| 647 | |
| 648 | C---generation of the 5 variables Y,Q2,SHAT,Z,PHI and Jacobian computation |
| 649 | |
| 650 | C---in the largest phase space avalaible for selected processes and |
| 651 | |
| 652 | C---filling of logical vector INSIDE to tag contributing ones |
| 653 | |
| 654 | CALL HWHBRN (*999) |
| 655 | |
| 656 | C---calculate differential cross section corresponding to the chosen |
| 657 | |
| 658 | C---variables and the weight for MC generation |
| 659 | |
| 660 | IF (IQK.EQ.0) THEN |
| 661 | |
| 662 | C---many subprocesses included |
| 663 | |
| 664 | DO I=1,18 |
| 665 | |
| 666 | FSIGMA(I)=ZERO |
| 667 | |
| 668 | ENDDO |
| 669 | |
| 670 | SIGSUM=ZERO |
| 671 | |
| 672 | DO I=IMIN,IMAX |
| 673 | |
| 674 | IF (INSIDE(I)) THEN |
| 675 | |
| 676 | IFL=I |
| 677 | |
| 678 | DSIGMA=ZERO |
| 679 | |
| 680 | CALL HWHBSG |
| 681 | |
| 682 | FSIGMA(I)=DSIGMA |
| 683 | |
| 684 | SIGSUM=SIGSUM+DSIGMA |
| 685 | |
| 686 | ENDIF |
| 687 | |
| 688 | ENDDO |
| 689 | |
| 690 | EVWGT=SIGSUM * AJACOB |
| 691 | |
| 692 | ELSE |
| 693 | |
| 694 | C---only one subprocess included |
| 695 | |
| 696 | CALL HWHBSG |
| 697 | |
| 698 | EVWGT= DSIGMA * AJACOB |
| 699 | |
| 700 | ENDIF |
| 701 | |
| 702 | IF (EVWGT.LT.ZERO) EVWGT=ZERO |
| 703 | |
| 704 | ENDIF |
| 705 | |
| 706 | 999 END |
| 707 | |
| 708 | CDECK ID>, HWHBKI. |
| 709 | |
| 710 | *CMZ :- -26/04/91 13.19.32 by Federico Carminati |
| 711 | |
| 712 | *-- Author : Giovanni Abbiendi & Luca Stanco |
| 713 | |
| 714 | C---------------------------------------------------------------------- |
| 715 | |
| 716 | SUBROUTINE HWHBKI |
| 717 | |
| 718 | C---------------------------------------------------------------------- |
| 719 | |
| 720 | C gives the fourmomenta in the laboratory system for the particles |
| 721 | |
| 722 | C of the hard 2-->3 subprocess, to match with HERWIG routines of |
| 723 | |
| 724 | C jet evolution. |
| 725 | |
| 726 | C---------------------------------------------------------------------- |
| 727 | |
| 728 | INCLUDE 'HERWIG61.INC' |
| 729 | |
| 730 | DOUBLE PRECISION HWUECM,HWUPCM,HWUSQR,LEP,Y,Q2,SHAT,Z,PHI,AJACOB, |
| 731 | |
| 732 | & DSIGMA,ME,MP,ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT, |
| 733 | |
| 734 | & PGAMMA(5),SG,MF1,MF2,EP,PP,EL,PL,E1,E2,Q1,COSBET,SINBET,COSTHE, |
| 735 | |
| 736 | & SINTHE,SINAZI,COSAZI,ROTAZI(3,3),EGAM,A,PPROT,MREMIN,PGAM,PEP(5), |
| 737 | |
| 738 | & COSPHI,SINPHI,ROT(3,3),EPROT,PROTON(5),MPART |
| 739 | |
| 740 | INTEGER IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,IPROO,I,IHAD,J,IS,ICMF |
| 741 | |
| 742 | LOGICAL CHARGD,INCLUD(18),INSIDE(18) |
| 743 | |
| 744 | EXTERNAL HWUECM,HWUPCM,HWUSQR |
| 745 | |
| 746 | COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF, |
| 747 | |
| 748 | & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL, |
| 749 | |
| 750 | & IPROO,CHARGD,INCLUD,INSIDE |
| 751 | |
| 752 | C |
| 753 | |
| 754 | IHAD=2 |
| 755 | |
| 756 | IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD) |
| 757 | |
| 758 | C---Set masses |
| 759 | |
| 760 | IF (CHARGD) THEN |
| 761 | |
| 762 | MPART=ZERO |
| 763 | |
| 764 | MF1=RMASS(IDHW(NHEP+5)) |
| 765 | |
| 766 | MF2=RMASS(IDHW(NHEP+6)) |
| 767 | |
| 768 | MREMIN=MP |
| 769 | |
| 770 | ELSE |
| 771 | |
| 772 | IS = IFL |
| 773 | |
| 774 | IF (IFL.EQ.164) IS=IQK |
| 775 | |
| 776 | MPART=ZERO |
| 777 | |
| 778 | IF (IFL.GE.7.AND.IFL.LE.18) MPART=RMASS(IFL-6) |
| 779 | |
| 780 | MF1=MFIN1(IS) |
| 781 | |
| 782 | MF2=MFIN2(IS) |
| 783 | |
| 784 | MREMIN = MREMIF(IS) |
| 785 | |
| 786 | ENDIF |
| 787 | |
| 788 | C---Calculation of kinematical variables for the generated event |
| 789 | |
| 790 | C in the center of mass frame of the incoming boson and parton |
| 791 | |
| 792 | C with parton along +z |
| 793 | |
| 794 | EGAM = HWUECM (SHAT, -Q2, MPART**2) |
| 795 | |
| 796 | PGAM = SQRT( EGAM**2 + Q2 ) |
| 797 | |
| 798 | EP = RSHAT-EGAM |
| 799 | |
| 800 | PP = PGAM |
| 801 | |
| 802 | A = (W2+Q2-MP**2)/TWO |
| 803 | |
| 804 | PPROT = (A*PGAM-EGAM*SQRT(A**2+MP**2*Q2))/Q2 |
| 805 | |
| 806 | IF (PPROT.LT.ZERO) CALL HWWARN('HWHBKI',101,*999) |
| 807 | |
| 808 | EPROT = SQRT(PPROT**2+MP**2) |
| 809 | |
| 810 | IF ((EPROT+PPROT).LT.(EP+PP)) CALL HWWARN('HWHBKI',102,*999) |
| 811 | |
| 812 | EL = ( PGAM / PPROT * SMA - Q2 ) / TWO |
| 813 | |
| 814 | + / (EGAM + PGAM / PPROT * EPROT) |
| 815 | |
| 816 | IF (EL.GT.ME) THEN |
| 817 | |
| 818 | PL = SQRT ( EL**2 - ME**2 ) |
| 819 | |
| 820 | ELSE |
| 821 | |
| 822 | CALL HWWARN ('HWHBKI',103,*999) |
| 823 | |
| 824 | ENDIF |
| 825 | |
| 826 | COSBET = (TWO * EPROT * EL - SMA) / (TWO * PPROT * PL) |
| 827 | |
| 828 | IF ( ABS(COSBET) .GE. ONE ) THEN |
| 829 | |
| 830 | COSBET = SIGN (ONE,COSBET) |
| 831 | |
| 832 | SINBET = ZERO |
| 833 | |
| 834 | ELSE |
| 835 | |
| 836 | SINBET = SQRT (ONE - COSBET**2) |
| 837 | |
| 838 | ENDIF |
| 839 | |
| 840 | SG = ME**2 + MPART**2 + Q2 + TWO * RSHAT * EL |
| 841 | |
| 842 | IF (SG.LE.(RSHAT+ML)**2 .OR. SG.GE.(RS-MREMIN)**2) |
| 843 | |
| 844 | + CALL HWWARN ('HWHBKI',104,*999) |
| 845 | |
| 846 | Q1 = HWUPCM( RSHAT, MF1, MF2) |
| 847 | |
| 848 | E1 = SQRT(Q1**2+MF1**2) |
| 849 | |
| 850 | E2 = SQRT(Q1**2+MF2**2) |
| 851 | |
| 852 | IF (Q1 .GT. ZERO) THEN |
| 853 | |
| 854 | COSTHE=(TWO*EP*E1 - Z*(SHAT+Q2))/(TWO*PP*Q1) |
| 855 | |
| 856 | IF (ABS(COSTHE) .GT. ONE) THEN |
| 857 | |
| 858 | COSTHE=SIGN(ONE,COSTHE) |
| 859 | |
| 860 | SINTHE=ZERO |
| 861 | |
| 862 | ELSE |
| 863 | |
| 864 | SINTHE=SQRT(ONE-COSTHE**2) |
| 865 | |
| 866 | ENDIF |
| 867 | |
| 868 | ELSE |
| 869 | |
| 870 | COSTHE=ZERO |
| 871 | |
| 872 | SINTHE=ONE |
| 873 | |
| 874 | ENDIF |
| 875 | |
| 876 | C---Initial lepton |
| 877 | |
| 878 | PHEP(1,NHEP+1)=PL*SINBET |
| 879 | |
| 880 | PHEP(2,NHEP+1)=ZERO |
| 881 | |
| 882 | PHEP(3,NHEP+1)=PL*COSBET |
| 883 | |
| 884 | PHEP(4,NHEP+1)=EL |
| 885 | |
| 886 | PHEP(5,NHEP+1)=RMASS(IDHW(1)) |
| 887 | |
| 888 | C---Initial Hadron |
| 889 | |
| 890 | PROTON(1)=ZERO |
| 891 | |
| 892 | PROTON(2)=ZERO |
| 893 | |
| 894 | PROTON(3)=PPROT |
| 895 | |
| 896 | PROTON(4)=EPROT |
| 897 | |
| 898 | CALL HWUMAS (PROTON) |
| 899 | |
| 900 | C---Initial parton |
| 901 | |
| 902 | PHEP(1,NHEP+2)=ZERO |
| 903 | |
| 904 | PHEP(2,NHEP+2)=ZERO |
| 905 | |
| 906 | PHEP(3,NHEP+2)=PP |
| 907 | |
| 908 | PHEP(4,NHEP+2)=EP |
| 909 | |
| 910 | PHEP(5,NHEP+2)=MPART |
| 911 | |
| 912 | C---HARD SUBPROCESS 2-->3 CENTRE OF MASS |
| 913 | |
| 914 | PHEP(1,NHEP+3)=PHEP(1,NHEP+1)+PHEP(1,NHEP+2) |
| 915 | |
| 916 | PHEP(2,NHEP+3)=PHEP(2,NHEP+1)+PHEP(2,NHEP+2) |
| 917 | |
| 918 | PHEP(3,NHEP+3)=PHEP(3,NHEP+1)+PHEP(3,NHEP+2) |
| 919 | |
| 920 | PHEP(4,NHEP+3)=PHEP(4,NHEP+1)+PHEP(4,NHEP+2) |
| 921 | |
| 922 | CALL HWUMAS ( PHEP(1,NHEP+3) ) |
| 923 | |
| 924 | C---Virtual boson |
| 925 | |
| 926 | PGAMMA(1)=ZERO |
| 927 | |
| 928 | PGAMMA(2)=ZERO |
| 929 | |
| 930 | PGAMMA(3)=-PGAM |
| 931 | |
| 932 | PGAMMA(4)=EGAM |
| 933 | |
| 934 | PGAMMA(5)=HWUSQR(Q2) |
| 935 | |
| 936 | C---Scattered lepton |
| 937 | |
| 938 | PHEP(1,NHEP+4)=PHEP(1,NHEP+1)-PGAMMA(1) |
| 939 | |
| 940 | PHEP(2,NHEP+4)=PHEP(2,NHEP+1)-PGAMMA(2) |
| 941 | |
| 942 | PHEP(3,NHEP+4)=PHEP(3,NHEP+1)-PGAMMA(3) |
| 943 | |
| 944 | PHEP(4,NHEP+4)=PHEP(4,NHEP+1)-PGAMMA(4) |
| 945 | |
| 946 | PHEP(5,NHEP+4)=RMASS(IDHW(1)) |
| 947 | |
| 948 | IF (CHARGD) PHEP(5,NHEP+4)=ZERO |
| 949 | |
| 950 | C---First Final parton: quark (or J/psi) in Boson-Gluon Fusion |
| 951 | |
| 952 | C--- quark or antiquark in QCD Compton |
| 953 | |
| 954 | PHEP(1,NHEP+5)=Q1*SINTHE*COS(PHI) |
| 955 | |
| 956 | PHEP(2,NHEP+5)=Q1*SINTHE*SIN(PHI) |
| 957 | |
| 958 | PHEP(3,NHEP+5)=Q1*COSTHE |
| 959 | |
| 960 | PHEP(4,NHEP+5)=E1 |
| 961 | |
| 962 | PHEP(5,NHEP+5)=MF1 |
| 963 | |
| 964 | C---Second Final parton: antiquark in Boson-Gluon Fusion |
| 965 | |
| 966 | C--- gluon in QCD Compton |
| 967 | |
| 968 | PHEP(1,NHEP+6)=-PHEP(1,NHEP+5) |
| 969 | |
| 970 | PHEP(2,NHEP+6)=-PHEP(2,NHEP+5) |
| 971 | |
| 972 | PHEP(3,NHEP+6)=-PHEP(3,NHEP+5) |
| 973 | |
| 974 | PHEP(4,NHEP+6)=E2 |
| 975 | |
| 976 | PHEP(5,NHEP+6)=MF2 |
| 977 | |
| 978 | C---Boost to lepton-hadron CM frame |
| 979 | |
| 980 | PEP(1) = PHEP(1,NHEP+1) |
| 981 | |
| 982 | PEP(2) = PHEP(2,NHEP+1) |
| 983 | |
| 984 | PEP(3) = PHEP(3,NHEP+1) + PPROT |
| 985 | |
| 986 | PEP(4) = PHEP(4,NHEP+1) + EPROT |
| 987 | |
| 988 | CALL HWUMAS (PEP) |
| 989 | |
| 990 | DO I=1,6 |
| 991 | |
| 992 | CALL HWULOF (PEP,PHEP(1,NHEP+I),PHEP(1,NHEP+I)) |
| 993 | |
| 994 | ENDDO |
| 995 | |
| 996 | CALL HWULOF (PEP,PROTON,PROTON) |
| 997 | |
| 998 | CALL HWULOF (PEP,PGAMMA,PGAMMA) |
| 999 | |
| 1000 | C---Rotation around y-axis to align lepton beam with z-axis |
| 1001 | |
| 1002 | COSPHI = PHEP(3,NHEP+1) / |
| 1003 | |
| 1004 | & SQRT( PHEP(1,NHEP+1)**2 + PHEP(3,NHEP+1)**2 ) |
| 1005 | |
| 1006 | SINPHI = PHEP(1,NHEP+1) / |
| 1007 | |
| 1008 | & SQRT( PHEP(1,NHEP+1)**2 + PHEP(3,NHEP+1)**2 ) |
| 1009 | |
| 1010 | DO I=1,3 |
| 1011 | |
| 1012 | DO J=1,3 |
| 1013 | |
| 1014 | ROT(I,J)=ZERO |
| 1015 | |
| 1016 | ENDDO |
| 1017 | |
| 1018 | ENDDO |
| 1019 | |
| 1020 | ROT(1,1) = COSPHI |
| 1021 | |
| 1022 | ROT(1,3) = -SINPHI |
| 1023 | |
| 1024 | ROT(2,2) = ONE |
| 1025 | |
| 1026 | ROT(3,1) = SINPHI |
| 1027 | |
| 1028 | ROT(3,3) = COSPHI |
| 1029 | |
| 1030 | DO I=1,6 |
| 1031 | |
| 1032 | CALL HWUROF (ROT,PHEP(1,NHEP+I),PHEP(1,NHEP+I)) |
| 1033 | |
| 1034 | ENDDO |
| 1035 | |
| 1036 | CALL HWUROF (ROT,PROTON,PROTON) |
| 1037 | |
| 1038 | CALL HWUROF (ROT,PGAMMA,PGAMMA) |
| 1039 | |
| 1040 | C---Boost to the LAB frame |
| 1041 | |
| 1042 | ICMF=3 |
| 1043 | |
| 1044 | DO I=1,6 |
| 1045 | |
| 1046 | CALL HWULOB (PHEP(1,ICMF),PHEP(1,NHEP+I),PHEP(1,NHEP+I)) |
| 1047 | |
| 1048 | ENDDO |
| 1049 | |
| 1050 | CALL HWULOB (PHEP(1,ICMF),PROTON,PROTON) |
| 1051 | |
| 1052 | CALL HWULOB (PHEP(1,ICMF),PGAMMA,PGAMMA) |
| 1053 | |
| 1054 | C---Random azimuthal rotation |
| 1055 | |
| 1056 | CALL HWRAZM (ONE,COSAZI,SINAZI) |
| 1057 | |
| 1058 | DO I=1,3 |
| 1059 | |
| 1060 | DO J=1,3 |
| 1061 | |
| 1062 | ROTAZI(I,J)=ZERO |
| 1063 | |
| 1064 | ENDDO |
| 1065 | |
| 1066 | ENDDO |
| 1067 | |
| 1068 | ROTAZI(1,1) = COSAZI |
| 1069 | |
| 1070 | ROTAZI(1,2) = SINAZI |
| 1071 | |
| 1072 | ROTAZI(2,1) = -SINAZI |
| 1073 | |
| 1074 | ROTAZI(2,2) = COSAZI |
| 1075 | |
| 1076 | ROTAZI(3,3) = ONE |
| 1077 | |
| 1078 | DO I=1,6 |
| 1079 | |
| 1080 | CALL HWUROF (ROTAZI,PHEP(1,NHEP+I),PHEP(1,NHEP+I)) |
| 1081 | |
| 1082 | ENDDO |
| 1083 | |
| 1084 | CALL HWUROF (ROTAZI,PROTON,PROTON) |
| 1085 | |
| 1086 | CALL HWUROF (ROTAZI,PGAMMA,PGAMMA) |
| 1087 | |
| 1088 | 999 END |
| 1089 | |
| 1090 | CDECK ID>, HWHBRN. |
| 1091 | |
| 1092 | *CMZ :- -03/07/95 19.02.12 by Giovanni Abbiendi |
| 1093 | |
| 1094 | *-- Author : Giovanni Abbiendi & Luca Stanco |
| 1095 | |
| 1096 | C----------------------------------------------------------------------- |
| 1097 | |
| 1098 | SUBROUTINE HWHBRN (*) |
| 1099 | |
| 1100 | C---------------------------------------------------------------------- |
| 1101 | |
| 1102 | C Returns a point in the phase space (Y,Q2,SHAT,Z,PHI) and the |
| 1103 | |
| 1104 | C corresponding Jacobian factor AJACOB |
| 1105 | |
| 1106 | C Fill the logical vector INSIDE to tag contributing subprocesses |
| 1107 | |
| 1108 | C to the cross-section |
| 1109 | |
| 1110 | C----------------------------------------------------------------------- |
| 1111 | |
| 1112 | INCLUDE 'HERWIG61.INC' |
| 1113 | |
| 1114 | DOUBLE PRECISION HWRUNI,HWR,HWUPCM,LEP,Y,Q2,SHAT,Z,PHI,AJACOB, |
| 1115 | |
| 1116 | & DSIGMA,ME,MP,ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT, |
| 1117 | |
| 1118 | & MF1,MF2,YMIN,YMAX,YJAC,Q2INF,Q2SUP,Q2JAC,EMW2,ZMIN,ZMAX,ZJAC, |
| 1119 | |
| 1120 | & GAMMA2,LAMBDA,PHIJAC,ZINT,ZLMIN,ZL,EMW,TMIN,TMAX,EMLMIN,EMLMAX, |
| 1121 | |
| 1122 | & SHMIN,EMMIF(18),EMMAF(18),WMIF(18),WMIN,MREMIN,YMIF(18),Q1CM(18), |
| 1123 | |
| 1124 | & Q2MAF(18),EMMAWF(18),ZMIF(18),ZMAF(18),PLMAX,PINC,SHINF,SHSUP, |
| 1125 | |
| 1126 | & SHJAC,CTHLIM,Q1,DETDSH,SRY,SRY0,SRY1 |
| 1127 | |
| 1128 | INTEGER IQK,IFLAVU,IFLAVD,I,IMIN,IMAX,IFL,IPROO,IHAD,NTRY,DEBUG |
| 1129 | |
| 1130 | LOGICAL CHARGD,INCLUD(18),INSIDE(18) |
| 1131 | |
| 1132 | EXTERNAL HWRUNI,HWR,HWUPCM |
| 1133 | |
| 1134 | SAVE EMLMIN,EMLMAX,EMMIF,EMMAF,MREMIN,MF1,MF2,YMIF, |
| 1135 | |
| 1136 | & YMIN,YMAX,WMIN,WMIF |
| 1137 | |
| 1138 | COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF, |
| 1139 | |
| 1140 | & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL, |
| 1141 | |
| 1142 | & IPROO,CHARGD,INCLUD,INSIDE |
| 1143 | |
| 1144 | EQUIVALENCE (EMW,RMASS(198)) |
| 1145 | |
| 1146 | C |
| 1147 | |
| 1148 | IHAD=2 |
| 1149 | |
| 1150 | IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD) |
| 1151 | |
| 1152 | C---Initialization |
| 1153 | |
| 1154 | IF (FSTWGT.OR.IHAD.NE.2) THEN |
| 1155 | |
| 1156 | ME = RMASS(IDHW(1)) |
| 1157 | |
| 1158 | MP = RMASS(IDHW(IHAD)) |
| 1159 | |
| 1160 | RS = PHEP(5,3) |
| 1161 | |
| 1162 | SMA = RS**2-ME**2-MP**2 |
| 1163 | |
| 1164 | PINC = HWUPCM(RS,ME,MP) |
| 1165 | |
| 1166 | C---Charged current |
| 1167 | |
| 1168 | IF (CHARGD) THEN |
| 1169 | |
| 1170 | ML=RMASS(IDHW(1)+1) |
| 1171 | |
| 1172 | YMAX = ONE - TWO*ML*MP / SMA |
| 1173 | |
| 1174 | YMAX = MIN(YMAX,YBMAX) |
| 1175 | |
| 1176 | MREMIN=MP |
| 1177 | |
| 1178 | IF (LEP.EQ.ONE) THEN |
| 1179 | |
| 1180 | MF1=RMASS(IFLAVD) |
| 1181 | |
| 1182 | MF2=RMASS(IFLAVU) |
| 1183 | |
| 1184 | ELSE |
| 1185 | |
| 1186 | MF1=RMASS(IFLAVU) |
| 1187 | |
| 1188 | MF2=RMASS(IFLAVD) |
| 1189 | |
| 1190 | ENDIF |
| 1191 | |
| 1192 | SHMIN = MF1**2+MF2**2 + TWO * PTMIN**2 + |
| 1193 | |
| 1194 | + TWO * SQRT(PTMIN**2+MF1**2) * SQRT(PTMIN**2+MF2**2) |
| 1195 | |
| 1196 | EMLMIN=MAX(EMMIN,SQRT(SHMIN)) |
| 1197 | |
| 1198 | EMLMAX=MIN(EMMAX,RS-ML-MREMIN) |
| 1199 | |
| 1200 | DEBUG=1 |
| 1201 | |
| 1202 | IF (EMLMIN.GT.EMLMAX) GOTO 888 |
| 1203 | |
| 1204 | WMIN=EMLMIN+MREMIN |
| 1205 | |
| 1206 | PLMAX=HWUPCM(RS,ML,WMIN) |
| 1207 | |
| 1208 | YMIN = ONE-TWO*(SQRT(PINC**2+MP**2)*SQRT(PLMAX**2+ML**2)+ |
| 1209 | |
| 1210 | + PINC*PLMAX)/SMA |
| 1211 | |
| 1212 | YMIN = MAX(YMIN,YBMIN) |
| 1213 | |
| 1214 | DEBUG=2 |
| 1215 | |
| 1216 | IF (YMIN.GT.YMAX) GOTO 888 |
| 1217 | |
| 1218 | ELSE |
| 1219 | |
| 1220 | C---Neutral current |
| 1221 | |
| 1222 | ML = ME |
| 1223 | |
| 1224 | YMAX = ONE - TWO*ML*MP / SMA |
| 1225 | |
| 1226 | YMAX = MIN(YMAX,YBMAX) |
| 1227 | |
| 1228 | DO I=1,18 |
| 1229 | |
| 1230 | YMIF(I)=ZERO |
| 1231 | |
| 1232 | EMMIF(I)=ZERO |
| 1233 | |
| 1234 | EMMAF(I)=ZERO |
| 1235 | |
| 1236 | WMIF(I)=ZERO |
| 1237 | |
| 1238 | IF (I.LE.8) THEN |
| 1239 | |
| 1240 | C---Boson-Gluon Fusion (also J/Psi) and QCD Compton with struck u or d |
| 1241 | |
| 1242 | MREMIF(I)=MP |
| 1243 | |
| 1244 | IF (I.LE.6) THEN |
| 1245 | |
| 1246 | MFIN1(I)=RMASS(I) |
| 1247 | |
| 1248 | MFIN2(I)=RMASS(I+6) |
| 1249 | |
| 1250 | ELSE |
| 1251 | |
| 1252 | MFIN1(I)=RMASS(I-6) |
| 1253 | |
| 1254 | MFIN2(I)=ZERO |
| 1255 | |
| 1256 | ENDIF |
| 1257 | |
| 1258 | ELSE |
| 1259 | |
| 1260 | C---QCD Compton with struck non-valence parton |
| 1261 | |
| 1262 | MREMIF(I)=MP+RMASS(I-6) |
| 1263 | |
| 1264 | MFIN1(I)=RMASS(I-6) |
| 1265 | |
| 1266 | MFIN2(I)=ZERO |
| 1267 | |
| 1268 | ENDIF |
| 1269 | |
| 1270 | ENDDO |
| 1271 | |
| 1272 | IF (IFL.EQ.164) THEN |
| 1273 | |
| 1274 | C---J/Psi |
| 1275 | |
| 1276 | MFIN1(7)=RMASS(164) |
| 1277 | |
| 1278 | MFIN2(7)=ZERO |
| 1279 | |
| 1280 | ENDIF |
| 1281 | |
| 1282 | C---y boundaries for different flavours and processes |
| 1283 | |
| 1284 | DO 100 I=IMIN,IMAX |
| 1285 | |
| 1286 | IF (INCLUD(I)) THEN |
| 1287 | |
| 1288 | MF1=MFIN1(I) |
| 1289 | |
| 1290 | MF2=MFIN2(I) |
| 1291 | |
| 1292 | MREMIN=MREMIF(I) |
| 1293 | |
| 1294 | SHMIN = MF1**2+MF2**2 + TWO * PTMIN**2 + |
| 1295 | |
| 1296 | + TWO * SQRT(PTMIN**2+MF1**2) * SQRT(PTMIN**2+MF2**2) |
| 1297 | |
| 1298 | EMMIF(I) = MAX(EMMIN,SQRT(SHMIN)) |
| 1299 | |
| 1300 | EMMAF(I) = MIN(EMMAX,RS-ML-MREMIN) |
| 1301 | |
| 1302 | IF (EMMIF(I).GT.EMMAF(I)) THEN |
| 1303 | |
| 1304 | INCLUD(I)=.FALSE. |
| 1305 | |
| 1306 | CALL HWWARN('HWHBRN',3,*999) |
| 1307 | |
| 1308 | GOTO 100 |
| 1309 | |
| 1310 | ENDIF |
| 1311 | |
| 1312 | WMIF(I) = EMMIF(I)+MREMIF(I) |
| 1313 | |
| 1314 | WMIN = WMIF(I) |
| 1315 | |
| 1316 | PLMAX = HWUPCM(RS,ML,WMIN) |
| 1317 | |
| 1318 | YMIF(I)=ONE-TWO*(SQRT(PINC**2+MP**2)*SQRT(PLMAX**2+ML**2)+ |
| 1319 | |
| 1320 | + PINC*PLMAX)/SMA |
| 1321 | |
| 1322 | IF (YMIF(I).GT.YMAX) THEN |
| 1323 | |
| 1324 | INCLUD(I)=.FALSE. |
| 1325 | |
| 1326 | CALL HWWARN('HWHBRN',4,*999) |
| 1327 | |
| 1328 | GOTO 100 |
| 1329 | |
| 1330 | ENDIF |
| 1331 | |
| 1332 | ENDIF |
| 1333 | |
| 1334 | 100 CONTINUE |
| 1335 | |
| 1336 | C---considering the largest boundaries |
| 1337 | |
| 1338 | EMLMIN=EMMIF(IMIN) |
| 1339 | |
| 1340 | EMLMAX=EMMAF(IMIN) |
| 1341 | |
| 1342 | IF (IPROO.EQ.3) THEN |
| 1343 | |
| 1344 | EMLMIN=MIN(EMMIF(IMIN),EMMIF(IMIN+6)) |
| 1345 | |
| 1346 | EMLMAX=MAX(EMMAF(IMIN),EMMAF(IMIN+6)) |
| 1347 | |
| 1348 | ENDIF |
| 1349 | |
| 1350 | DEBUG=3 |
| 1351 | |
| 1352 | IF (EMLMIN.GT.EMLMAX) GOTO 888 |
| 1353 | |
| 1354 | YMIN=YMIF(IMIN) |
| 1355 | |
| 1356 | IF (IPROO.EQ.3) YMIN=MIN(YMIF(IMIN),YMIF(IMIN+6)) |
| 1357 | |
| 1358 | YMIN = MAX(YMIN,YBMIN) |
| 1359 | |
| 1360 | DEBUG=4 |
| 1361 | |
| 1362 | IF (YMIN.GT.YMAX) GOTO 888 |
| 1363 | |
| 1364 | WMIN = WMIF(IMIN) |
| 1365 | |
| 1366 | MREMIN = MREMIF(IMIN) |
| 1367 | |
| 1368 | MF1=MFIN1(IMIN) |
| 1369 | |
| 1370 | MF2=MFIN2(IMIN) |
| 1371 | |
| 1372 | IF (IPROO.EQ.3) THEN |
| 1373 | |
| 1374 | WMIN = MIN(WMIF(IMIN),WMIF(IMIN+6)) |
| 1375 | |
| 1376 | MREMIN = MIN(MREMIF(IMIN),MREMIF(IMIN+6)) |
| 1377 | |
| 1378 | ENDIF |
| 1379 | |
| 1380 | ENDIF |
| 1381 | |
| 1382 | ENDIF |
| 1383 | |
| 1384 | C---Random generation in largest phase space |
| 1385 | |
| 1386 | Y=ZERO |
| 1387 | |
| 1388 | Q2=ZERO |
| 1389 | |
| 1390 | SHAT=ZERO |
| 1391 | |
| 1392 | Z=ZERO |
| 1393 | |
| 1394 | PHI=ZERO |
| 1395 | |
| 1396 | AJACOB=ZERO |
| 1397 | |
| 1398 | C---y generation |
| 1399 | |
| 1400 | IF (.NOT.CHARGD) THEN |
| 1401 | |
| 1402 | IF (IFL.LE.5.OR.(IFL.GE.7.AND.IFL.LE.18)) THEN |
| 1403 | |
| 1404 | SRY0 = SQRT(YMIN) |
| 1405 | |
| 1406 | SRY1 = SQRT(YMAX) |
| 1407 | |
| 1408 | SRY = HWRUNI(0,SRY0,SRY1) |
| 1409 | |
| 1410 | Y = SRY**2 |
| 1411 | |
| 1412 | YJAC = TWO*SRY*(SRY1-SRY0) |
| 1413 | |
| 1414 | ELSEIF (IFL.EQ.6) THEN |
| 1415 | |
| 1416 | Y = SQRT(HWRUNI(0,YMIN**2,YMAX**2)) |
| 1417 | |
| 1418 | YJAC = HALF * (YMAX**2-YMIN**2) / Y |
| 1419 | |
| 1420 | ELSEIF (IFL.EQ.164) THEN |
| 1421 | |
| 1422 | C---in J/psi photoproduction Y and Q2 are given by the Equivalent Photon |
| 1423 | |
| 1424 | C Approximation |
| 1425 | |
| 1426 | 10 NTRY=0 |
| 1427 | |
| 1428 | 20 NTRY=NTRY+1 |
| 1429 | |
| 1430 | IF (NTRY.GT.NETRY) CALL HWWARN('HWHBRN',50,*10) |
| 1431 | |
| 1432 | Y = (YMIN/YMAX)**HWR()*YMAX |
| 1433 | |
| 1434 | IF (ONE+(ONE-Y)**2.LT.TWO*HWR()) GOTO 20 |
| 1435 | |
| 1436 | YJAC=(TWO*LOG(YMAX/YMIN)-TWO*(YMAX-YMIN) |
| 1437 | |
| 1438 | & +HALF*(YMAX**2-YMIN**2)) |
| 1439 | |
| 1440 | ENDIF |
| 1441 | |
| 1442 | ELSE |
| 1443 | |
| 1444 | IF (IPRO.EQ.5) THEN |
| 1445 | |
| 1446 | Y = EXP(HWRUNI(0,LOG(YMIN),LOG(YMAX))) |
| 1447 | |
| 1448 | YJAC = Y * LOG(YMAX/YMIN) |
| 1449 | |
| 1450 | ELSE |
| 1451 | |
| 1452 | Y = HWRUNI(0,YMIN,YMAX) |
| 1453 | |
| 1454 | YJAC = YMAX - YMIN |
| 1455 | |
| 1456 | ENDIF |
| 1457 | |
| 1458 | ENDIF |
| 1459 | |
| 1460 | C---Q**2 generation |
| 1461 | |
| 1462 | Q2INF = ME**2*Y**2 / (ONE-Y) |
| 1463 | |
| 1464 | Q2SUP = MP**2 + SMA*Y - WMIN**2 |
| 1465 | |
| 1466 | IF (IFL.EQ.164) THEN |
| 1467 | |
| 1468 | Q2INF = MAX(Q2INF,Q2WWMN) |
| 1469 | |
| 1470 | Q2SUP = MIN(Q2SUP,Q2WWMX) |
| 1471 | |
| 1472 | ELSE |
| 1473 | |
| 1474 | Q2INF = MAX(Q2INF,Q2MIN) |
| 1475 | |
| 1476 | Q2SUP = MIN(Q2SUP,Q2MAX) |
| 1477 | |
| 1478 | ENDIF |
| 1479 | |
| 1480 | DEBUG=5 |
| 1481 | |
| 1482 | IF (Q2INF .GT. Q2SUP) GOTO 888 |
| 1483 | |
| 1484 | C |
| 1485 | |
| 1486 | IF (.NOT.CHARGD) THEN |
| 1487 | |
| 1488 | IF (IFL.EQ.164) THEN |
| 1489 | |
| 1490 | Q2 = EXP(HWRUNI(0,LOG(Q2INF),LOG(Q2SUP))) |
| 1491 | |
| 1492 | Q2JAC = LOG(Q2SUP/Q2INF) |
| 1493 | |
| 1494 | ELSEIF (Q2INF.LT.RMASS(4)**2) THEN |
| 1495 | |
| 1496 | Q2 = EXP(HWRUNI(0,LOG(Q2INF),LOG(Q2SUP))) |
| 1497 | |
| 1498 | Q2JAC = Q2 * LOG(Q2SUP/Q2INF) |
| 1499 | |
| 1500 | ELSE |
| 1501 | |
| 1502 | Q2 = Q2INF*Q2SUP/HWRUNI(0,Q2INF,Q2SUP) |
| 1503 | |
| 1504 | Q2JAC = Q2**2 * (Q2SUP-Q2INF)/(Q2SUP*Q2INF) |
| 1505 | |
| 1506 | ENDIF |
| 1507 | |
| 1508 | ELSE |
| 1509 | |
| 1510 | EMW2=EMW**2 |
| 1511 | |
| 1512 | Q2=(Q2INF+EMW2)*(Q2SUP+EMW2)/(HWRUNI(0,Q2INF,Q2SUP)+EMW2)-EMW2 |
| 1513 | |
| 1514 | Q2JAC=(Q2+EMW2)**2*(Q2SUP-Q2INF)/((Q2SUP+EMW2)*(Q2INF+EMW2)) |
| 1515 | |
| 1516 | ENDIF |
| 1517 | |
| 1518 | W2 = MP**2 + SMA*Y - Q2 |
| 1519 | |
| 1520 | C---s_hat generation |
| 1521 | |
| 1522 | SHINF = EMLMIN **2 |
| 1523 | |
| 1524 | SHSUP = (MIN(SQRT(W2)-MREMIN,EMLMAX))**2 |
| 1525 | |
| 1526 | DEBUG=6 |
| 1527 | |
| 1528 | IF (SHINF .GT. SHSUP) GOTO 888 |
| 1529 | |
| 1530 | C |
| 1531 | |
| 1532 | IF (IPRO.EQ.91) THEN |
| 1533 | |
| 1534 | IF (.NOT.CHARGD) THEN |
| 1535 | |
| 1536 | SHAT = SHINF*SHSUP/HWRUNI(0,SHINF,SHSUP) |
| 1537 | |
| 1538 | SHJAC = SHAT**2 * (SHSUP-SHINF)/(SHSUP*SHINF) |
| 1539 | |
| 1540 | ELSE |
| 1541 | |
| 1542 | SHAT = EXP(HWRUNI(0,LOG(SHINF),LOG(SHSUP))) |
| 1543 | |
| 1544 | SHJAC = SHAT*(LOG(SHSUP/SHINF)) |
| 1545 | |
| 1546 | ENDIF |
| 1547 | |
| 1548 | ELSE |
| 1549 | |
| 1550 | EMW2=EMW**2 |
| 1551 | |
| 1552 | IF (SHINF.GT.EMW2+10*GAMW*EMW) THEN |
| 1553 | |
| 1554 | SHAT = SHINF*SHSUP/HWRUNI(0,SHINF,SHSUP) |
| 1555 | |
| 1556 | SHJAC = SHAT**2 * (SHSUP-SHINF)/(SHSUP*SHINF) |
| 1557 | |
| 1558 | ELSEIF (SHSUP.LT.EMW2-10*EMW*GAMW) THEN |
| 1559 | |
| 1560 | SHAT = HWRUNI(0,SHINF,SHSUP) |
| 1561 | |
| 1562 | SHJAC = SHSUP-SHINF |
| 1563 | |
| 1564 | ELSE |
| 1565 | |
| 1566 | TMIN=ATAN((SHINF-EMW2)/(GAMW*EMW)) |
| 1567 | |
| 1568 | TMAX=ATAN((SHSUP-EMW2)/(GAMW*EMW)) |
| 1569 | |
| 1570 | SHAT = GAMW*EMW*TAN(HWRUNI(0,TMIN,TMAX))+EMW2 |
| 1571 | |
| 1572 | SHJAC=((SHAT-EMW2)**2+(GAMW*EMW)**2)/(GAMW*EMW)*(TMAX-TMIN) |
| 1573 | |
| 1574 | ENDIF |
| 1575 | |
| 1576 | ENDIF |
| 1577 | |
| 1578 | DETDSH = ONE/SMA/Y |
| 1579 | |
| 1580 | SHJAC=SHJAC*DETDSH |
| 1581 | |
| 1582 | RSHAT = SQRT (SHAT) |
| 1583 | |
| 1584 | C--- z generation |
| 1585 | |
| 1586 | ZMIN = 10E10 |
| 1587 | |
| 1588 | ZMAX = -ONE |
| 1589 | |
| 1590 | IF (.NOT.CHARGD) THEN |
| 1591 | |
| 1592 | DO I=1,18 |
| 1593 | |
| 1594 | Q1CM(I) = ZERO |
| 1595 | |
| 1596 | ZMIF(I) = ZERO |
| 1597 | |
| 1598 | ZMAF(I) = ZERO |
| 1599 | |
| 1600 | ENDDO |
| 1601 | |
| 1602 | DO 150 I=IMIN,IMAX |
| 1603 | |
| 1604 | IF (INCLUD(I)) THEN |
| 1605 | |
| 1606 | Q1CM(I) = HWUPCM( RSHAT, MFIN1(I), MFIN2(I) ) |
| 1607 | |
| 1608 | IF (Q1CM(I) .LT. PTMIN) THEN |
| 1609 | |
| 1610 | ZMAF(I)=-ONE |
| 1611 | |
| 1612 | GOTO 150 |
| 1613 | |
| 1614 | ENDIF |
| 1615 | |
| 1616 | CTHLIM = SQRT(ONE - (PTMIN / Q1CM(I))**2) |
| 1617 | |
| 1618 | GAMMA2 = SHAT + MFIN1(I)**2 - MFIN2(I)**2 |
| 1619 | |
| 1620 | LAMBDA = (SHAT-MFIN1(I)**2-MFIN2(I)**2)**2 - |
| 1621 | |
| 1622 | + 4.D0*MFIN1(I)**2*MFIN2(I)**2 |
| 1623 | |
| 1624 | ZMIF(I) = (GAMMA2 - SQRT(LAMBDA)*CTHLIM)/TWO/SHAT |
| 1625 | |
| 1626 | ZMIF(I) = MAX(ZMIF(I),ZERO) |
| 1627 | |
| 1628 | ZMAF(I) = (GAMMA2 + SQRT(LAMBDA)*CTHLIM)/TWO/SHAT |
| 1629 | |
| 1630 | ZMAF(I) = MIN(ZMAF(I),ONE) |
| 1631 | |
| 1632 | ZMIN = MIN( ZMIN, ZMIF(I) ) |
| 1633 | |
| 1634 | ZMAX = MAX( ZMAX, ZMAF(I) ) |
| 1635 | |
| 1636 | ENDIF |
| 1637 | |
| 1638 | 150 CONTINUE |
| 1639 | |
| 1640 | IF (IFL.EQ.164) ZMAX=MIN(ZMAX,ZJMAX) |
| 1641 | |
| 1642 | ELSE |
| 1643 | |
| 1644 | Q1 = HWUPCM(RSHAT,MF1,MF2) |
| 1645 | |
| 1646 | DEBUG=7 |
| 1647 | |
| 1648 | IF (Q1.LT.PTMIN) GOTO 888 |
| 1649 | |
| 1650 | CTHLIM = SQRT(ONE-(PTMIN/Q1)**2) |
| 1651 | |
| 1652 | GAMMA2 = SHAT+MF1**2-MF2**2 |
| 1653 | |
| 1654 | LAMBDA = (SHAT-MF1**2-MF2**2)**2-4.D0*MF1**2*MF2**2 |
| 1655 | |
| 1656 | ZMIN = (GAMMA2-SQRT(LAMBDA)*CTHLIM)/TWO/SHAT |
| 1657 | |
| 1658 | ZMIN = MAX(ZMIN,1D-6) |
| 1659 | |
| 1660 | ZMAX = (GAMMA2+SQRT(LAMBDA)*CTHLIM)/TWO/SHAT |
| 1661 | |
| 1662 | ZMAX = MIN(ZMAX,ONE-1D-6) |
| 1663 | |
| 1664 | ENDIF |
| 1665 | |
| 1666 | DEBUG=8 |
| 1667 | |
| 1668 | IF (ZMIN .GT. ZMAX) GOTO 888 |
| 1669 | |
| 1670 | ZLMIN = LOG(ZMIN/(ONE-ZMIN)) |
| 1671 | |
| 1672 | ZINT = LOG(ZMAX/(ONE-ZMAX)) - LOG(ZMIN/(ONE-ZMIN)) |
| 1673 | |
| 1674 | ZL = ZLMIN+HWR()*ZINT |
| 1675 | |
| 1676 | Z = EXP(ZL)/(ONE+EXP(ZL)) |
| 1677 | |
| 1678 | ZJAC = Z*(ONE-Z)*ZINT |
| 1679 | |
| 1680 | C |
| 1681 | |
| 1682 | DEBUG=9 |
| 1683 | |
| 1684 | IF ((Y.LT.YMIN.OR.Y.GT.YMAX).OR.(Q2.LT.Q2INF.OR.Q2.GT.Q2SUP).OR. |
| 1685 | |
| 1686 | + (SHAT.LT.SHINF.OR.SHAT.GT.SHSUP).OR.(Z.LT.ZMIN.OR.Z.GT.ZMAX)) |
| 1687 | |
| 1688 | + GOTO 888 |
| 1689 | |
| 1690 | C---Phi generation |
| 1691 | |
| 1692 | PHI = HWRUNI(0,ZERO,2*PIFAC) |
| 1693 | |
| 1694 | PHIJAC = 2 * PIFAC |
| 1695 | |
| 1696 | IF (IFL.EQ.164) PHIJAC=ONE |
| 1697 | |
| 1698 | C |
| 1699 | |
| 1700 | AJACOB = YJAC * Q2JAC * SHJAC * ZJAC * PHIJAC |
| 1701 | |
| 1702 | C |
| 1703 | |
| 1704 | IF (IQK.NE.0.OR.IPRO.EQ.5) GOTO 999 |
| 1705 | |
| 1706 | C---contributing subprocesses: filling of logical vector INSIDE |
| 1707 | |
| 1708 | DO I=1,18 |
| 1709 | |
| 1710 | INSIDE(I)=.FALSE. |
| 1711 | |
| 1712 | Q2MAF(I)=ZERO |
| 1713 | |
| 1714 | EMMAWF(I)=ZERO |
| 1715 | |
| 1716 | ENDDO |
| 1717 | |
| 1718 | DO 200 I=IMIN,IMAX |
| 1719 | |
| 1720 | IF (INCLUD(I)) THEN |
| 1721 | |
| 1722 | IF ( Y.LT.YMIF(I) ) GOTO 200 |
| 1723 | |
| 1724 | C |
| 1725 | |
| 1726 | Q2MAF(I) = MP**2 + SMA*Y - WMIF(I)**2 |
| 1727 | |
| 1728 | Q2MAF(I) = MIN( Q2MAF(I), Q2MAX) |
| 1729 | |
| 1730 | IF (Q2INF .GT. Q2MAF(I)) GOTO 200 |
| 1731 | |
| 1732 | IF (Q2.LT.Q2INF .OR. Q2.GT.Q2MAF(I)) GOTO 200 |
| 1733 | |
| 1734 | C |
| 1735 | |
| 1736 | EMMAWF(I) = SQRT(W2) - MREMIF(I) |
| 1737 | |
| 1738 | EMMAWF(I) = MIN( EMMAWF(I), EMLMAX ) |
| 1739 | |
| 1740 | C |
| 1741 | |
| 1742 | IF (EMMIF(I) .GT. EMMAWF(I)) GOTO 200 |
| 1743 | |
| 1744 | IF (SHAT.LT.EMMIF(I)**2.OR.SHAT.GT.EMMAWF(I)**2) GOTO 200 |
| 1745 | |
| 1746 | C |
| 1747 | |
| 1748 | IF (ZMIF(I) .GT. ZMAF(I)) GOTO 200 |
| 1749 | |
| 1750 | IF (Z.LT.ZMIF(I) .OR. Z.GT.ZMAF(I)) GOTO 200 |
| 1751 | |
| 1752 | INSIDE(I)=.TRUE. |
| 1753 | |
| 1754 | ENDIF |
| 1755 | |
| 1756 | 200 CONTINUE |
| 1757 | |
| 1758 | 999 RETURN |
| 1759 | |
| 1760 | 888 EVWGT=ZERO |
| 1761 | |
| 1762 | C---UNCOMMENT THIS LINE TO GET A DEBUGGING WARNING FOR NO PHASE-SPACE |
| 1763 | |
| 1764 | C CALL HWWARN('HWHBRN',DEBUG,*777) |
| 1765 | |
| 1766 | 777 RETURN 1 |
| 1767 | |
| 1768 | END |
git://git.uio.no / u / mrichter / AliRoot.git / blame