| fe4da5cc |
1 | |
| 2 | C********************************************************************* |
| 3 | |
| 4 | SUBROUTINE PYDIFF |
| 5 | |
| 6 | C...Handles diffractive and elastic scattering. |
| 7 | COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) |
| 8 | COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) |
| 9 | COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) |
| 10 | COMMON/PYINT1/MINT(400),VINT(400) |
| 11 | SAVE /LUJETS/,/LUDAT1/ |
| 12 | SAVE /PYPARS/,/PYINT1/ |
| 13 | DOUBLE PRECISION DBETAZ |
| 14 | |
| 15 | C...Reset K, P and V vectors. Store incoming particles. |
| 16 | DO 110 JT=1,MSTP(126)+10 |
| 17 | I=MINT(83)+JT |
| 18 | DO 100 J=1,5 |
| 19 | K(I,J)=0 |
| 20 | P(I,J)=0. |
| 21 | V(I,J)=0. |
| 22 | 100 CONTINUE |
| 23 | 110 CONTINUE |
| 24 | N=MINT(84) |
| 25 | MINT(3)=0 |
| 26 | MINT(21)=0 |
| 27 | MINT(22)=0 |
| 28 | MINT(23)=0 |
| 29 | MINT(24)=0 |
| 30 | MINT(4)=4 |
| 31 | DO 130 JT=1,2 |
| 32 | I=MINT(83)+JT |
| 33 | K(I,1)=21 |
| 34 | K(I,2)=MINT(10+JT) |
| 35 | DO 120 J=1,5 |
| 36 | P(I,J)=VINT(285+5*JT+J) |
| 37 | 120 CONTINUE |
| 38 | 130 CONTINUE |
| 39 | MINT(6)=2 |
| 40 | |
| 41 | C...Subprocess; kinematics. |
| 42 | SQLAM=(VINT(2)-VINT(63)-VINT(64))**2-4.*VINT(63)*VINT(64) |
| 43 | PZ=SQRT(SQLAM)/(2.*VINT(1)) |
| 44 | DO 200 JT=1,2 |
| 45 | I=MINT(83)+JT |
| 46 | PE=(VINT(2)+VINT(62+JT)-VINT(65-JT))/(2.*VINT(1)) |
| 47 | KFH=MINT(102+JT) |
| 48 | |
| 49 | C...Elastically scattered particle. |
| 50 | IF(MINT(16+JT).LE.0) THEN |
| 51 | N=N+1 |
| 52 | K(N,1)=1 |
| 53 | K(N,2)=KFH |
| 54 | K(N,3)=I+2 |
| 55 | P(N,3)=PZ*(-1)**(JT+1) |
| 56 | P(N,4)=PE |
| 57 | P(N,5)=SQRT(VINT(62+JT)) |
| 58 | |
| 59 | C...Decay rho from elastic scattering of gamma with sin**2(theta) |
| 60 | C...distribution of decay products (in rho rest frame). |
| 61 | IF(KFH.EQ.113.AND.MINT(10+JT).EQ.22.AND.MSTP(102).EQ.1) THEN |
| 62 | NSAV=N |
| 63 | DBETAZ=DBLE(P(N,3))/SQRT(DBLE(P(N,3))**2+DBLE(P(N,5))**2) |
| 64 | P(N,3)=0. |
| 65 | P(N,4)=P(N,5) |
| 66 | CALL LUDECY(NSAV) |
| 67 | IF(N.EQ.NSAV+2.AND.IABS(K(NSAV+1,2)).EQ.211) THEN |
| 68 | PHI=ULANGL(P(NSAV+1,1),P(NSAV+1,2)) |
| 69 | CALL LUDBRB(NSAV+1,NSAV+2,0.,-PHI,0D0,0D0,0D0) |
| 70 | THE=ULANGL(P(NSAV+1,3),P(NSAV+1,1)) |
| 71 | CALL LUDBRB(NSAV+1,NSAV+2,-THE,0.,0D0,0D0,0D0) |
| 72 | 140 CTHE=2.*RLU(0)-1. |
| 73 | IF(1.-CTHE**2.LT.RLU(0)) GOTO 140 |
| 74 | CALL LUDBRB(NSAV+1,NSAV+2,ACOS(CTHE),PHI,0D0,0D0,0D0) |
| 75 | ENDIF |
| 76 | CALL LUDBRB(NSAV,NSAV+2,0.,0.,0D0,0D0,DBETAZ) |
| 77 | ENDIF |
| 78 | |
| 79 | C...Diffracted particle: low-mass system to two particles. |
| 80 | ELSEIF(VINT(62+JT).LT.(VINT(66+JT)+PARP(103))**2) THEN |
| 81 | N=N+2 |
| 82 | K(N-1,1)=1 |
| 83 | K(N,1)=1 |
| 84 | K(N-1,3)=I+2 |
| 85 | K(N,3)=I+2 |
| 86 | PMMAS=SQRT(VINT(62+JT)) |
| 87 | NTRY=0 |
| 88 | 150 NTRY=NTRY+1 |
| 89 | IF(NTRY.LT.20) THEN |
| 90 | MINT(105)=MINT(102+JT) |
| 91 | MINT(109)=MINT(106+JT) |
| 92 | CALL PYSPLI(KFH,21,KFL1,KFL2) |
| 93 | CALL LUKFDI(KFL1,0,KFL3,KF1) |
| 94 | IF(KF1.EQ.0) GOTO 150 |
| 95 | CALL LUKFDI(KFL2,-KFL3,KFLDUM,KF2) |
| 96 | IF(KF2.EQ.0) GOTO 150 |
| 97 | ELSE |
| 98 | KF1=KFH |
| 99 | KF2=111 |
| 100 | ENDIF |
| 101 | PM1=ULMASS(KF1) |
| 102 | PM2=ULMASS(KF2) |
| 103 | IF(PM1+PM2+PARJ(64).GT.PMMAS) GOTO 150 |
| 104 | K(N-1,2)=KF1 |
| 105 | K(N,2)=KF2 |
| 106 | P(N-1,5)=PM1 |
| 107 | P(N,5)=PM2 |
| 108 | PZP=SQRT(MAX(0.,(PMMAS**2-PM1**2-PM2**2)**2-4.*PM1**2*PM2**2))/ |
| 109 | & (2.*PMMAS) |
| 110 | P(N-1,3)=PZP |
| 111 | P(N,3)=-PZP |
| 112 | P(N-1,4)=SQRT(PM1**2+PZP**2) |
| 113 | P(N,4)=SQRT(PM2**2+PZP**2) |
| 114 | CALL LUDBRB(N-1,N,ACOS(2.*RLU(0)-1.),PARU(2)*RLU(0),0D0,0D0,0D0) |
| 115 | DBETAZ=DBLE(PZ)*(-1)**(JT+1)/SQRT(DBLE(PZ)**2+DBLE(PMMAS)**2) |
| 116 | CALL LUDBRB(N-1,N,0.,0.,0D0,0D0,DBETAZ) |
| 117 | |
| 118 | C...Diffracted particle: valence quark kicked out. |
| 119 | ELSEIF(MSTP(101).EQ.1.OR.(MSTP(101).EQ.3.AND.RLU(0).LT.PARP(101))) |
| 120 | &THEN |
| 121 | N=N+2 |
| 122 | K(N-1,1)=2 |
| 123 | K(N,1)=1 |
| 124 | K(N-1,3)=I+2 |
| 125 | K(N,3)=I+2 |
| 126 | MINT(105)=MINT(102+JT) |
| 127 | MINT(109)=MINT(106+JT) |
| 128 | CALL PYSPLI(KFH,21,K(N,2),K(N-1,2)) |
| 129 | P(N-1,5)=ULMASS(K(N-1,2)) |
| 130 | P(N,5)=ULMASS(K(N,2)) |
| 131 | SQLAM=(VINT(62+JT)-P(N-1,5)**2-P(N,5)**2)**2- |
| 132 | & 4.*P(N-1,5)**2*P(N,5)**2 |
| 133 | P(N-1,3)=(PE*SQRT(SQLAM)+PZ*(VINT(62+JT)+P(N-1,5)**2- |
| 134 | & P(N,5)**2))/(2.*VINT(62+JT))*(-1)**(JT+1) |
| 135 | P(N-1,4)=SQRT(P(N-1,3)**2+P(N-1,5)**2) |
| 136 | P(N,3)=PZ*(-1)**(JT+1)-P(N-1,3) |
| 137 | P(N,4)=SQRT(P(N,3)**2+P(N,5)**2) |
| 138 | |
| 139 | C...Diffracted particle: gluon kicked out. |
| 140 | ELSE |
| 141 | N=N+3 |
| 142 | K(N-2,1)=2 |
| 143 | K(N-1,1)=2 |
| 144 | K(N,1)=1 |
| 145 | K(N-2,3)=I+2 |
| 146 | K(N-1,3)=I+2 |
| 147 | K(N,3)=I+2 |
| 148 | MINT(105)=MINT(102+JT) |
| 149 | MINT(109)=MINT(106+JT) |
| 150 | CALL PYSPLI(KFH,21,K(N,2),K(N-2,2)) |
| 151 | K(N-1,2)=21 |
| 152 | P(N-2,5)=ULMASS(K(N-2,2)) |
| 153 | P(N-1,5)=0. |
| 154 | P(N,5)=ULMASS(K(N,2)) |
| 155 | C...Energy distribution for particle into two jets. |
| 156 | 160 IMB=1 |
| 157 | IF(MOD(KFH/1000,10).NE.0) IMB=2 |
| 158 | CHIK=PARP(92+2*IMB) |
| 159 | IF(MSTP(92).LE.1) THEN |
| 160 | IF(IMB.EQ.1) CHI=RLU(0) |
| 161 | IF(IMB.EQ.2) CHI=1.-SQRT(RLU(0)) |
| 162 | ELSEIF(MSTP(92).EQ.2) THEN |
| 163 | CHI=1.-RLU(0)**(1./(1.+CHIK)) |
| 164 | ELSEIF(MSTP(92).EQ.3) THEN |
| 165 | CUT=2.*0.3/VINT(1) |
| 166 | 170 CHI=RLU(0)**2 |
| 167 | IF((CHI**2/(CHI**2+CUT**2))**0.25*(1.-CHI)**CHIK.LT. |
| 168 | & RLU(0)) GOTO 170 |
| 169 | ELSEIF(MSTP(92).EQ.4) THEN |
| 170 | CUT=2.*0.3/VINT(1) |
| 171 | CUTR=(1.+SQRT(1.+CUT**2))/CUT |
| 172 | 180 CHIR=CUT*CUTR**RLU(0) |
| 173 | CHI=(CHIR**2-CUT**2)/(2.*CHIR) |
| 174 | IF((1.-CHI)**CHIK.LT.RLU(0)) GOTO 180 |
| 175 | ELSE |
| 176 | CUT=2.*0.3/VINT(1) |
| 177 | CUTA=CUT**(1.-PARP(98)) |
| 178 | CUTB=(1.+CUT)**(1.-PARP(98)) |
| 179 | 190 CHI=(CUTA+RLU(0)*(CUTB-CUTA))**(1./(1.-PARP(98))) |
| 180 | IF(((CHI+CUT)**2/(2.*(CHI**2+CUT**2)))** |
| 181 | & (0.5*PARP(98))*(1.-CHI)**CHIK.LT.RLU(0)) GOTO 190 |
| 182 | ENDIF |
| 183 | IF(CHI.LT.P(N,5)**2/VINT(62+JT).OR.CHI.GT.1.-P(N-2,5)**2/ |
| 184 | & VINT(62+JT)) GOTO 160 |
| 185 | SQM=P(N-2,5)**2/(1.-CHI)+P(N,5)**2/CHI |
| 186 | IF((SQRT(SQM)+PARJ(32))**2.GE.VINT(62+JT)) GOTO 160 |
| 187 | PZI=(PE*(VINT(62+JT)-SQM)+PZ*(VINT(62+JT)+SQM))/ |
| 188 | & (2.*VINT(62+JT)) |
| 189 | PEI=SQRT(PZI**2+SQM) |
| 190 | PQQP=(1.-CHI)*(PEI+PZI) |
| 191 | P(N-2,3)=0.5*(PQQP-P(N-2,5)**2/PQQP)*(-1)**(JT+1) |
| 192 | P(N-2,4)=SQRT(P(N-2,3)**2+P(N-2,5)**2) |
| 193 | P(N-1,4)=0.5*(VINT(62+JT)-SQM)/(PEI+PZI) |
| 194 | P(N-1,3)=P(N-1,4)*(-1)**JT |
| 195 | P(N,3)=PZI*(-1)**(JT+1)-P(N-2,3) |
| 196 | P(N,4)=SQRT(P(N,3)**2+P(N,5)**2) |
| 197 | ENDIF |
| 198 | |
| 199 | C...Documentation lines. |
| 200 | K(I+2,1)=21 |
| 201 | IF(MINT(16+JT).EQ.0) K(I+2,2)=KFH |
| 202 | IF(MINT(16+JT).NE.0) K(I+2,2)=10*(KFH/10) |
| 203 | K(I+2,3)=I |
| 204 | P(I+2,3)=PZ*(-1)**(JT+1) |
| 205 | P(I+2,4)=PE |
| 206 | P(I+2,5)=SQRT(VINT(62+JT)) |
| 207 | 200 CONTINUE |
| 208 | |
| 209 | C...Rotate outgoing partons/particles using cos(theta). |
| 210 | IF(VINT(23).LT.0.9) THEN |
| 211 | CALL LUDBRB(MINT(83)+3,N,ACOS(VINT(23)),VINT(24),0D0,0D0,0D0) |
| 212 | ELSE |
| 213 | CALL LUDBRB(MINT(83)+3,N,ASIN(VINT(59)),VINT(24),0D0,0D0,0D0) |
| 214 | ENDIF |
| 215 | |
| 216 | RETURN |
| 217 | END |
git://git.uio.no / u / mrichter / AliRoot.git / blame