pythiaeRHIC/pythia/pydiff.f

0001
0002 C*********************************************************************
0003
0004 C...PYDIFF
0005 C...Handles diffractive and elastic scattering.
0006
0007       SUBROUTINE PYDIFF
0008
0009 C...Double precision and integer declarations.
0010       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0011       IMPLICIT INTEGER(I-N)
0012       INTEGER PYK,PYCHGE,PYCOMP
0013 C...Commonblocks.
0014       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0015       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0016       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0017       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0018       COMMON/PYINT1/MINT(400),VINT(400)
0019       SAVE /PYJETS/,/PYDAT1/,/PYPARS/,/PYINT1/
0020       include "mc_set.inc"
0021
0022 C...Reset K, P and V vectors. Store incoming particles.
0023       DO 110 JT=1,MSTP(126)+10
0024         I=MINT(83)+JT
0025         DO 100 J=1,5
0026           K(I,J)=0
0027           P(I,J)=0D0
0028           V(I,J)=0D0
0029   100   CONTINUE
0030   110 CONTINUE
0031       N=MINT(84)
0032       MINT(3)=0
0033       MINT(21)=0
0034       MINT(22)=0
0035       MINT(23)=0
0036       MINT(24)=0
0037       MINT(4)=4
0038       DO 130 JT=1,2
0039         I=MINT(83)+JT
0040         K(I,1)=21
0041         K(I,2)=MINT(10+JT)
0042         DO 120 J=1,5
0043           P(I,J)=VINT(285+5*JT+J)
0044   120   CONTINUE
0045   130 CONTINUE
0046       MINT(6)=2
0047
0048 C...Subprocess; kinematics.
0049       SQLAM=(VINT(2)-VINT(63)-VINT(64))**2-4D0*VINT(63)*VINT(64)
0050       PZ=SQRT(SQLAM)/(2D0*VINT(1))
0051       DO 200 JT=1,2
0052         I=MINT(83)+JT
0053         PE=(VINT(2)+VINT(62+JT)-VINT(65-JT))/(2D0*VINT(1))
0054         KFH=MINT(102+JT)
0055
0056 C...Elastically scattered particle. (Except elastic GVMD states.)
0057         IF(MINT(16+JT).LE.0.AND.(MINT(10+JT).NE.22.OR.
0058      &  MINT(106+JT).NE.3)) THEN
0059           N=N+1
0060           K(N,1)=1
0061           K(N,2)=KFH
0062           K(N,3)=I+2
0063           P(N,3)=PZ*(-1)**(JT+1)
0064           P(N,4)=PE
0065           P(N,5)=SQRT(VINT(62+JT))
0066
0067 C...Decay rho from elastic scattering of gamma with sin**2(theta)
0068 C...distribution of decay products (in rho rest frame).
0069           IF(KFH.EQ.113.AND.MINT(10+JT).EQ.22.AND.MSTP(102).EQ.1) THEN
0070             NSAV=N
0071             DBETAZ=P(N,3)/SQRT(P(N,3)**2+P(N,5)**2)
0072             P(N,3)=0D0
0073             P(N,4)=P(N,5)
0074             CALL PYDECY(NSAV)
0075             IF(N.EQ.NSAV+2.AND.IABS(K(NSAV+1,2)).EQ.211) THEN
0076               PHI=PYANGL(P(NSAV+1,1),P(NSAV+1,2))
0077               CALL PYROBO(NSAV+1,NSAV+2,0D0,-PHI,0D0,0D0,0D0)
0078               THE=PYANGL(P(NSAV+1,3),P(NSAV+1,1))
0079               CALL PYROBO(NSAV+1,NSAV+2,-THE,0D0,0D0,0D0,0D0)
0080   140         CTHE=2D0*PYR(0)-1D0
0081
0082 C... Changing parameters for R_rho with values corresponding to W<7
0083 C... (measured by HERMES     R_rho=1/eps * r0400/(1. - r0400)
0084               PMVIRT=PMAS(PYCOMP(113),1)
0085               R_rho=PARP(165)*(VINT(307)/(PMVIRT**2))**PARP(166)
0086               BEAMAS=PYMASS(11)
0087 C new epsilon (f_L/f_T) as used in pysigh.F with proton mass
0088               eps=1D0/(1D0+(VINT(309)**2*(1D0-2D0*BEAMAS**2/
0089      &            VINT(307)))/(2D0/(1D0+VINT(307)/VINT(309)**2/
0090      &            ebeamEnucl**2)*(1D0-VINT(309)-
0091      &            (VINT(307)/4D0/ebeamEnucl**2))))
0092               r0400=eps*R_rho / ( 1. + eps * R_rho)
0093               w_ang=0.75d0*(1.d0-r0400+(3.d0*r0400-1.d0)*cthe**2.)
0094               if( r0400 .le. 1.d0/3.d0 ) then
0095                  w_ang_max_x = 0.d0
0096               else
0097                  w_ang_max_x = 1.d0
0098               endif
0099               w_ang_max= 0.75d0*(1.d0-r0400+(3.d0*r0400-1.d0)
0100      $             *w_ang_max_x**2.)
0101               IF(PYR(0).gt.w_ang/w_ang_max) GOTO 140
0102 C              IF(1D0-CTHE**2.LT.PYR(0)) GOTO 140
0103               CALL PYROBO(NSAV+1,NSAV+2,ACOS(CTHE),PHI,0D0,0D0,0D0)
0104             ENDIF
0105             CALL PYROBO(NSAV,NSAV+2,0D0,0D0,0D0,0D0,DBETAZ)
0106           ENDIF
0107
0108 C...Diffracted particle: low-mass system to two particles.
0109         ELSEIF(VINT(62+JT).LT.(VINT(66+JT)+PARP(103))**2) THEN
0110           N=N+2
0111           K(N-1,1)=1
0112           K(N,1)=1
0113           K(N-1,3)=I+2
0114           K(N,3)=I+2
0115           PMMAS=SQRT(VINT(62+JT))
0116           NTRY=0
0117   150     NTRY=NTRY+1
0118           IF(NTRY.LT.20) THEN
0119             MINT(105)=MINT(102+JT)
0120             MINT(109)=MINT(106+JT)
0121             CALL PYSPLI(KFH,21,KFL1,KFL2)
0122             CALL PYKFDI(KFL1,0,KFL3,KF1)
0123             IF(KF1.EQ.0) GOTO 150
0124             CALL PYKFDI(KFL2,-KFL3,KFLDUM,KF2)
0125             IF(KF2.EQ.0) GOTO 150
0126           ELSE
0127             KF1=KFH
0128             KF2=111
0129           ENDIF
0130           PM1=PYMASS(KF1)
0131           PM2=PYMASS(KF2)
0132           IF(PM1+PM2+PARJ(64).GT.PMMAS) GOTO 150
0133           K(N-1,2)=KF1
0134           K(N,2)=KF2
0135           P(N-1,5)=PM1
0136           P(N,5)=PM2
0137           PZP=SQRT(MAX(0D0,(PMMAS**2-PM1**2-PM2**2)**2-
0138      &    4D0*PM1**2*PM2**2))/(2D0*PMMAS)
0139           P(N-1,3)=PZP
0140           P(N,3)=-PZP
0141           P(N-1,4)=SQRT(PM1**2+PZP**2)
0142           P(N,4)=SQRT(PM2**2+PZP**2)
0143           CALL PYROBO(N-1,N,ACOS(2D0*PYR(0)-1D0),PARU(2)*PYR(0),
0144      &    0D0,0D0,0D0)
0145           DBETAZ=PZ*(-1)**(JT+1)/SQRT(PZ**2+PMMAS**2)
0146           CALL PYROBO(N-1,N,0D0,0D0,0D0,0D0,DBETAZ)
0147
0148 C...Diffracted particle: valence quark kicked out.
0149         ELSEIF(MSTP(101).EQ.1.OR.(MSTP(101).EQ.3.AND.PYR(0).LT.
0150      &    PARP(101))) THEN
0151           N=N+2
0152           K(N-1,1)=2
0153           K(N,1)=1
0154           K(N-1,3)=I+2
0155           K(N,3)=I+2
0156           MINT(105)=MINT(102+JT)
0157           MINT(109)=MINT(106+JT)
0158           CALL PYSPLI(KFH,21,K(N,2),K(N-1,2))
0159           P(N-1,5)=PYMASS(K(N-1,2))
0160           P(N,5)=PYMASS(K(N,2))
0161           SQLAM=(VINT(62+JT)-P(N-1,5)**2-P(N,5)**2)**2-
0162      &    4D0*P(N-1,5)**2*P(N,5)**2
0163           P(N-1,3)=(PE*SQRT(SQLAM)+PZ*(VINT(62+JT)+P(N-1,5)**2-
0164      &    P(N,5)**2))/(2D0*VINT(62+JT))*(-1)**(JT+1)
0165           P(N-1,4)=SQRT(P(N-1,3)**2+P(N-1,5)**2)
0166           P(N,3)=PZ*(-1)**(JT+1)-P(N-1,3)
0167           P(N,4)=SQRT(P(N,3)**2+P(N,5)**2)
0168
0169 C...Diffracted particle: gluon kicked out.
0170         ELSE
0171           N=N+3
0172           K(N-2,1)=2
0173           K(N-1,1)=2
0174           K(N,1)=1
0175           K(N-2,3)=I+2
0176           K(N-1,3)=I+2
0177           K(N,3)=I+2
0178           MINT(105)=MINT(102+JT)
0179           MINT(109)=MINT(106+JT)
0180           CALL PYSPLI(KFH,21,K(N,2),K(N-2,2))
0181           K(N-1,2)=21
0182           P(N-2,5)=PYMASS(K(N-2,2))
0183           P(N-1,5)=0D0
0184           P(N,5)=PYMASS(K(N,2))
0185 C...Energy distribution for particle into two jets.
0186   160     IMB=1
0187           IF(MOD(KFH/1000,10).NE.0) IMB=2
0188           CHIK=PARP(92+2*IMB)
0189           IF(MSTP(92).LE.1) THEN
0190             IF(IMB.EQ.1) CHI=PYR(0)
0191             IF(IMB.EQ.2) CHI=1D0-SQRT(PYR(0))
0192           ELSEIF(MSTP(92).EQ.2) THEN
0193             CHI=1D0-PYR(0)**(1D0/(1D0+CHIK))
0194           ELSEIF(MSTP(92).EQ.3) THEN
0195             CUT=2D0*0.3D0/VINT(1)
0196   170       CHI=PYR(0)**2
0197             IF((CHI**2/(CHI**2+CUT**2))**0.25D0*(1D0-CHI)**CHIK.LT.
0198      &      PYR(0)) GOTO 170
0199           ELSEIF(MSTP(92).EQ.4) THEN
0200             CUT=2D0*0.3D0/VINT(1)
0201             CUTR=(1D0+SQRT(1D0+CUT**2))/CUT
0202   180       CHIR=CUT*CUTR**PYR(0)
0203             CHI=(CHIR**2-CUT**2)/(2D0*CHIR)
0204             IF((1D0-CHI)**CHIK.LT.PYR(0)) GOTO 180
0205           ELSE
0206             CUT=2D0*0.3D0/VINT(1)
0207             CUTA=CUT**(1D0-PARP(98))
0208             CUTB=(1D0+CUT)**(1D0-PARP(98))
0209   190       CHI=(CUTA+PYR(0)*(CUTB-CUTA))**(1D0/(1D0-PARP(98)))
0210             IF(((CHI+CUT)**2/(2D0*(CHI**2+CUT**2)))**
0211      &      (0.5D0*PARP(98))*(1D0-CHI)**CHIK.LT.PYR(0)) GOTO 190
0212           ENDIF
0213           IF(CHI.LT.P(N,5)**2/VINT(62+JT).OR.CHI.GT.1D0-P(N-2,5)**2/
0214      &    VINT(62+JT)) GOTO 160
0215           SQM=P(N-2,5)**2/(1D0-CHI)+P(N,5)**2/CHI
0216           PZI=(PE*(VINT(62+JT)-SQM)+PZ*(VINT(62+JT)+SQM))/
0217      &    (2D0*VINT(62+JT))
0218           PEI=SQRT(PZI**2+SQM)
0219           PQQP=(1D0-CHI)*(PEI+PZI)
0220           P(N-2,3)=0.5D0*(PQQP-P(N-2,5)**2/PQQP)*(-1)**(JT+1)
0221           P(N-2,4)=SQRT(P(N-2,3)**2+P(N-2,5)**2)
0222           P(N-1,4)=0.5D0*(VINT(62+JT)-SQM)/(PEI+PZI)
0223           P(N-1,3)=P(N-1,4)*(-1)**JT
0224           P(N,3)=PZI*(-1)**(JT+1)-P(N-2,3)
0225           P(N,4)=SQRT(P(N,3)**2+P(N,5)**2)
0226         ENDIF
0227
0228 C...Documentation lines.
0229         K(I+2,1)=21
0230         IF(MINT(16+JT).EQ.0) K(I+2,2)=KFH
0231         IF(MINT(16+JT).NE.0.OR.(MINT(10+JT).EQ.22.AND.
0232      &  MINT(106+JT).EQ.3)) K(I+2,2)=ISIGN(9900000,KFH)+10*(KFH/10)
0233         K(I+2,3)=I
0234         P(I+2,3)=PZ*(-1)**(JT+1)
0235         P(I+2,4)=PE
0236         P(I+2,5)=SQRT(VINT(62+JT))
0237   200 CONTINUE
0238
0239 C...Rotate outgoing partons/particles using cos(theta).
0240       IF(VINT(23).LT.0.9D0) THEN
0241         CALL PYROBO(MINT(83)+3,N,ACOS(VINT(23)),VINT(24),0D0,0D0,0D0)
0242       ELSE
0243         CALL PYROBO(MINT(83)+3,N,ASIN(VINT(59)),VINT(24),0D0,0D0,0D0)
0244       ENDIF
0245
0246       RETURN
0247       END