hijing/src/pyhidiff.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE PYHIDIFF
0005
0006 C...Handles diffractive and elastic scattering.
0007       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0008       SAVE /LUJETS/
0009       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0010       SAVE /LUDAT1/
0011       COMMON/PYHIPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0012       SAVE /PYHIPARS/
0013       COMMON/PYHIINT1/MINT(400),VINT(400)
0014       SAVE /PYHIINT1/
0015
0016 C...Reset K, P and V vectors. Store incoming particles.
0017       DO 100 JT=1,MSTP(126)+10
0018       I=MINT(83)+JT
0019       DO 100 J=1,5
0020       K(I,J)=0
0021       P(I,J)=0.
0022   100 V(I,J)=0.
0023       N=MINT(84)
0024       MINT(3)=0
0025       MINT(21)=0
0026       MINT(22)=0
0027       MINT(23)=0
0028       MINT(24)=0
0029       MINT(4)=4
0030       DO 110 JT=1,2
0031       I=MINT(83)+JT
0032       K(I,1)=21
0033       K(I,2)=MINT(10+JT)
0034       P(I,5)=VINT(2+JT)
0035       P(I,3)=VINT(5)*(-1)**(JT+1)
0036   110 P(I,4)=SQRT(P(I,3)**2+P(I,5)**2)
0037       MINT(6)=2
0038
0039 C...Subprocess; kinematics.
0040       ISUB=MINT(1)
0041       SQLAM=(VINT(2)-VINT(63)-VINT(64))**2-4.*VINT(63)*VINT(64)
0042       PZ=SQRT(SQLAM)/(2.*VINT(1))
0043       DO 150 JT=1,2
0044       I=MINT(83)+JT
0045       PE=(VINT(2)+VINT(62+JT)-VINT(65-JT))/(2.*VINT(1))
0046
0047 C...Elastically scattered particle.
0048       IF(MINT(16+JT).LE.0) THEN
0049         N=N+1
0050         K(N,1)=1
0051         K(N,2)=K(I,2)
0052         K(N,3)=I+2
0053         P(N,3)=PZ*(-1)**(JT+1)
0054         P(N,4)=PE
0055         P(N,5)=P(I,5)
0056
0057 C...Diffracted particle: valence quark kicked out.
0058       ELSEIF(MSTP(101).EQ.1) THEN
0059         N=N+2
0060         K(N-1,1)=2
0061         K(N,1)=1
0062         K(N-1,3)=I+2
0063         K(N,3)=I+2
0064         CALL PYHISPLI(K(I,2),21,K(N,2),K(N-1,2))
0065         P(N-1,5)=ULMASS(K(N-1,2))
0066         P(N,5)=ULMASS(K(N,2))
0067         SQLAM=(VINT(62+JT)-P(N-1,5)**2-P(N,5)**2)**2-
0068      &  4.*P(N-1,5)**2*P(N,5)**2
0069         P(N-1,3)=(PE*SQRT(SQLAM)+PZ*(VINT(62+JT)+P(N-1,5)**2-
0070      &  P(N,5)**2))/(2.*VINT(62+JT))*(-1)**(JT+1)
0071         P(N-1,4)=SQRT(P(N-1,3)**2+P(N-1,5)**2)
0072         P(N,3)=PZ*(-1)**(JT+1)-P(N-1,3)
0073         P(N,4)=SQRT(P(N,3)**2+P(N,5)**2)
0074
0075 C...Diffracted particle: gluon kicked out.
0076       ELSE
0077         N=N+3
0078         K(N-2,1)=2
0079         K(N-1,1)=2
0080         K(N,1)=1
0081         K(N-2,3)=I+2
0082         K(N-1,3)=I+2
0083         K(N,3)=I+2
0084         CALL PYHISPLI(K(I,2),21,K(N,2),K(N-2,2))
0085         K(N-1,2)=21
0086         P(N-2,5)=ULMASS(K(N-2,2))
0087         P(N-1,5)=0.
0088         P(N,5)=ULMASS(K(N,2))
0089 C...Energy distribution for particle into two jets.
0090   120   IMB=1
0091         IF(MOD(K(I,2)/1000,10).NE.0) IMB=2
0092         CHIK=PARP(92+2*IMB)
0093         IF(MSTP(92).LE.1) THEN
0094           IF(IMB.EQ.1) CHI=RLU(0)
0095           IF(IMB.EQ.2) CHI=1.-SQRT(RLU(0))
0096         ELSEIF(MSTP(92).EQ.2) THEN
0097           CHI=1.-RLU(0)**(1./(1.+CHIK))
0098         ELSEIF(MSTP(92).EQ.3) THEN
0099           CUT=2.*0.3/VINT(1)
0100   130     CHI=RLU(0)**2
0101           IF((CHI**2/(CHI**2+CUT**2))**0.25*(1.-CHI)**CHIK.LT.
0102      &    RLU(0)) GOTO 130
0103         ELSE
0104           CUT=2.*0.3/VINT(1)
0105           CUTR=(1.+SQRT(1.+CUT**2))/CUT
0106   140     CHIR=CUT*CUTR**RLU(0)
0107           CHI=(CHIR**2-CUT**2)/(2.*CHIR)
0108           IF((1.-CHI)**CHIK.LT.RLU(0)) GOTO 140
0109         ENDIF
0110         IF(CHI.LT.P(N,5)**2/VINT(62+JT).OR.CHI.GT.1.-P(N-2,5)**2/
0111      &  VINT(62+JT)) GOTO 120
0112         SQM=P(N-2,5)**2/(1.-CHI)+P(N,5)**2/CHI
0113         IF((SQRT(SQM)+PARJ(32))**2.GE.VINT(62+JT)) GOTO 120
0114         PZI=(PE*(VINT(62+JT)-SQM)+PZ*(VINT(62+JT)+SQM))/
0115      &  (2.*VINT(62+JT))
0116         PEI=SQRT(PZI**2+SQM)
0117         PQQP=(1.-CHI)*(PEI+PZI)
0118         P(N-2,3)=0.5*(PQQP-P(N-2,5)**2/PQQP)*(-1)**(JT+1)
0119         P(N-2,4)=SQRT(P(N-2,3)**2+P(N-2,5)**2)
0120         P(N-1,3)=(PZ-PZI)*(-1)**(JT+1)
0121         P(N-1,4)=ABS(P(N-1,3))
0122         P(N,3)=PZI*(-1)**(JT+1)-P(N-2,3)
0123         P(N,4)=SQRT(P(N,3)**2+P(N,5)**2)
0124       ENDIF
0125
0126 C...Documentation lines.
0127       K(I+2,1)=21
0128       IF(MINT(16+JT).EQ.0) K(I+2,2)=MINT(10+JT)
0129       IF(MINT(16+JT).NE.0) K(I+2,2)=10*(MINT(10+JT)/10)
0130       K(I+2,3)=I
0131       P(I+2,3)=PZ*(-1)**(JT+1)
0132       P(I+2,4)=PE
0133       P(I+2,5)=SQRT(VINT(62+JT))
0134   150 CONTINUE
0135
0136 C...Rotate outgoing partons/particles using cos(theta).
0137       CALL LUDBRB(MINT(83)+3,N,ACOS(VINT(23)),VINT(24),0D0,0D0,0D0)
0138
0139       RETURN
0140       END