pythiaeRHIC/pythia/pydisg.f

0001
0002 C*********************************************************************
0003
0004 C...PYDISG
0005 C...Set up a DIS process as gamma* + f -> f, with beam remnant
0006 C...and showering added consecutively. Photon flux by the PYGAGA
0007 C...routine (if at all).
0008
0009       SUBROUTINE PYDISG
0010
0011 C...Double precision and integer declarations.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013       IMPLICIT INTEGER(I-N)
0014       INTEGER PYK,PYCHGE,PYCOMP
0015 C...Parameter statement to help give large particle numbers.
0016       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0017      &KEXCIT=4000000,KDIMEN=5000000)
0018 C...Commonblocks.
0019       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0020       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0021       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0022       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0023       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0024       COMMON/PYINT1/MINT(400),VINT(400)
0025       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYSUBS/,/PYPARS/,/PYINT1/
0026 C...Local arrays.
0027       DIMENSION PMS(4)
0028
0029 C...Choice of subprocess, number of documentation lines
0030       IDOC=7
0031       MINT(3)=IDOC-6
0032       MINT(4)=IDOC
0033       IPU1=MINT(84)+1
0034       IPU2=MINT(84)+2
0035       IPU3=MINT(84)+3
0036       ISIDE=1
0037       IF(MINT(107).EQ.4) ISIDE=2
0038
0039 C...Reset K, P and V vectors. Store incoming particles
0040       DO 110 JT=1,MSTP(126)+20
0041         I=MINT(83)+JT
0042         DO 100 J=1,5
0043           K(I,J)=0
0044           P(I,J)=0D0
0045           V(I,J)=0D0
0046   100   CONTINUE
0047   110 CONTINUE
0048       DO 130 JT=1,2
0049         I=MINT(83)+JT
0050         K(I,1)=21
0051         K(I,2)=MINT(10+JT)
0052         DO 120 J=1,5
0053           P(I,J)=VINT(285+5*JT+J)
0054   120   CONTINUE
0055   130 CONTINUE
0056       MINT(6)=2
0057
0058 C...Store incoming partons in hadronic CM-frame
0059       DO 140 JT=1,2
0060         I=MINT(84)+JT
0061         K(I,1)=14
0062         K(I,2)=MINT(14+JT)
0063         K(I,3)=MINT(83)+2+JT
0064   140 CONTINUE
0065       IF(MINT(15).EQ.22) THEN
0066         P(MINT(84)+1,3)=0.5D0*(VINT(1)+VINT(307)/VINT(1))
0067         P(MINT(84)+1,4)=0.5D0*(VINT(1)-VINT(307)/VINT(1))
0068         P(MINT(84)+1,5)=-SQRT(VINT(307))
0069         P(MINT(84)+2,3)=-0.5D0*VINT(307)/VINT(1)
0070         P(MINT(84)+2,4)=0.5D0*VINT(307)/VINT(1)
0071         KFRES=MINT(16)
0072         ISIDE=2
0073       ELSE
0074         P(MINT(84)+1,3)=0.5D0*VINT(308)/VINT(1)
0075         P(MINT(84)+1,4)=0.5D0*VINT(308)/VINT(1)
0076         P(MINT(84)+2,3)=-0.5D0*(VINT(1)+VINT(308)/VINT(1))
0077         P(MINT(84)+2,4)=0.5D0*(VINT(1)-VINT(308)/VINT(1))
0078         P(MINT(84)+1,5)=-SQRT(VINT(308))
0079         KFRES=MINT(15)
0080         ISIDE=1
0081       ENDIF
0082       SIDESG=(-1D0)**(ISIDE-1)
0083
0084 C...Copy incoming partons to documentation lines.
0085       DO 170 JT=1,2
0086         I1=MINT(83)+4+JT
0087         I2=MINT(84)+JT
0088         K(I1,1)=21
0089         K(I1,2)=K(I2,2)
0090         K(I1,3)=I1-2
0091         DO 150 J=1,5
0092           P(I1,J)=P(I2,J)
0093   150   CONTINUE
0094
0095 C...Second copy for partons before ISR shower, since no such.
0096         I1=MINT(83)+2+JT
0097         K(I1,1)=21
0098         K(I1,2)=K(I2,2)
0099         K(I1,3)=I1-2
0100         DO 160 J=1,5
0101           P(I1,J)=P(I2,J)
0102   160   CONTINUE
0103   170 CONTINUE
0104
0105 C...Define initial partons.
0106       NTRY=0
0107   180 NTRY=NTRY+1
0108       IF(NTRY.GT.100) THEN
0109         MINT(51)=1
0110         RETURN
0111       ENDIF
0112
0113 C...Scattered quark in hadronic CM frame.
0114       I=MINT(83)+7
0115       K(IPU3,1)=3
0116       K(IPU3,2)=KFRES
0117       K(IPU3,3)=I
0118       P(IPU3,5)=PYMASS(KFRES)
0119       P(IPU3,3)=P(IPU1,3)+P(IPU2,3)
0120       P(IPU3,4)=P(IPU1,4)+P(IPU2,4)
0121       P(IPU3,5)=0D0
0122       K(I,1)=21
0123       K(I,2)=KFRES
0124       K(I,3)=MINT(83)+4+ISIDE
0125       P(I,3)=P(IPU3,3)
0126       P(I,4)=P(IPU3,4)
0127       P(I,5)=P(IPU3,5)
0128       N=IPU3
0129       MINT(21)=KFRES
0130       MINT(22)=0
0131
0132 C...No primordial kT, or chosen according to truncated Gaussian or
0133 C...exponential, or (for photon) predetermined or power law.
0134   190 IF(MINT(40+ISIDE).EQ.2.AND.MINT(10+ISIDE).NE.22) THEN
0135         IF(MSTP(91).LE.0) THEN
0136           PT=0D0
0137         ELSEIF(MSTP(91).EQ.1) THEN
0138           PT=PARP(91)*SQRT(-LOG(PYR(0)))
0139         ELSE
0140           RPT1=PYR(0)
0141           RPT2=PYR(0)
0142           PT=-PARP(92)*LOG(RPT1*RPT2)
0143         ENDIF
0144         IF(PT.GT.PARP(93)) GOTO 190
0145       ELSEIF(MINT(106+ISIDE).EQ.3) THEN
0146         PTA=SQRT(VINT(282+ISIDE))
0147         PTB=0D0
0148         IF(MSTP(66).EQ.5.AND.MSTP(93).EQ.1) THEN
0149           PTB=PARP(99)*SQRT(-LOG(PYR(0)))
0150         ELSEIF(MSTP(66).EQ.5.AND.MSTP(93).EQ.2) THEN
0151           RPT1=PYR(0)
0152           RPT2=PYR(0)
0153           PTB=-PARP(99)*LOG(RPT1*RPT2)
0154         ENDIF
0155         IF(PTB.GT.PARP(100)) GOTO 190
0156         PT=SQRT(PTA**2+PTB**2+2D0*PTA*PTB*COS(PARU(2)*PYR(0)))
0157         IF(NTRY.GT.10) PT=PT*0.8D0**(NTRY-10)
0158       ELSEIF(IABS(MINT(14+ISIDE)).LE.8.OR.MINT(14+ISIDE).EQ.21) THEN
0159         IF(MSTP(93).LE.0) THEN
0160           PT=0D0
0161         ELSEIF(MSTP(93).EQ.1) THEN
0162           PT=PARP(99)*SQRT(-LOG(PYR(0)))
0163         ELSEIF(MSTP(93).EQ.2) THEN
0164           RPT1=PYR(0)
0165           RPT2=PYR(0)
0166           PT=-PARP(99)*LOG(RPT1*RPT2)
0167         ELSEIF(MSTP(93).EQ.3) THEN
0168           HA=PARP(99)**2
0169           HB=PARP(100)**2
0170           PT=SQRT(MAX(0D0,HA*(HA+HB)/(HA+HB-PYR(0)*HB)-HA))
0171         ELSE
0172           HA=PARP(99)**2
0173           HB=PARP(100)**2
0174           IF(MSTP(93).EQ.5) HB=MIN(VINT(48),PARP(100)**2)
0175           PT=SQRT(MAX(0D0,HA*((HA+HB)/HA)**PYR(0)-HA))
0176         ENDIF
0177         IF(PT.GT.PARP(100)) GOTO 190
0178       ELSE
0179         PT=0D0
0180       ENDIF
0181       VINT(156+ISIDE)=PT
0182       PHI=PARU(2)*PYR(0)
0183       P(IPU3,1)=PT*COS(PHI)
0184       P(IPU3,2)=PT*SIN(PHI)
0185       P(IPU3,4)=SQRT(P(IPU3,5)**2+PT**2+P(IPU3,3)**2)
0186       PMS(3-ISIDE)=P(IPU3,5)**2+P(IPU3,1)**2+P(IPU3,2)**2
0187       PCP=P(IPU3,4)+ABS(P(IPU3,3))
0188
0189 C...Find one or two beam remnants.
0190       MINT(105)=MINT(102+ISIDE)
0191       MINT(109)=MINT(106+ISIDE)
0192       CALL PYSPLI(MINT(10+ISIDE),MINT(12+ISIDE),KFLCH,KFLSP)
0193       IF(MINT(51).NE.0) THEN
0194         MINT(51)=0
0195         GOTO 180
0196       ENDIF
0197
0198 C...Store first remnant parton, with colour info and kinematics.
0199       I=N+1
0200       K(I,1)=1
0201       K(I,2)=KFLSP
0202       K(I,3)=MINT(83)+ISIDE
0203       P(I,5)=PYMASS(K(I,2))
0204       KCOL=KCHG(PYCOMP(KFLSP),2)
0205       IF(KCOL.NE.0) THEN
0206         K(I,1)=3
0207         KFLS=(3-KCOL*ISIGN(1,KFLSP))/2
0208         K(I,KFLS+3)=MSTU(5)*IPU3
0209         K(IPU3,6-KFLS)=MSTU(5)*I
0210         ICOLR=I
0211       ENDIF
0212       IF(KFLCH.EQ.0) THEN
0213         P(I,1)=-P(IPU3,1)
0214         P(I,2)=-P(IPU3,2)
0215         PMS(ISIDE)=P(I,5)**2+P(I,1)**2+P(I,2)**2
0216         P(I,3)=-P(IPU3,3)
0217         P(I,4)=SQRT(PMS(ISIDE)+P(I,3)**2)
0218         PRP=P(I,4)+ABS(P(I,3))
0219
0220 C...When extra remnant parton or hadron: store extra remnant.
0221       ELSE
0222         I=I+1
0223         K(I,1)=1
0224         K(I,2)=KFLCH
0225         K(I,3)=MINT(83)+ISIDE
0226         P(I,5)=PYMASS(K(I,2))
0227         KCOL=KCHG(PYCOMP(KFLCH),2)
0228         IF(KCOL.NE.0) THEN
0229           K(I,1)=3
0230           KFLS=(3-KCOL*ISIGN(1,KFLCH))/2
0231           K(I,KFLS+3)=MSTU(5)*IPU3
0232           K(IPU3,6-KFLS)=MSTU(5)*I
0233           ICOLR=I
0234         ENDIF
0235
0236 C...Relative transverse momentum when two remnants.
0237         LOOP=0
0238   200   LOOP=LOOP+1
0239         CALL PYPTDI(1,P(I-1,1),P(I-1,2))
0240         P(I-1,1)=P(I-1,1)-0.5D0*P(IPU3,1)
0241         P(I-1,2)=P(I-1,2)-0.5D0*P(IPU3,2)
0242         PMS(3)=P(I-1,5)**2+P(I-1,1)**2+P(I-1,2)**2
0243         P(I,1)=-P(IPU3,1)-P(I-1,1)
0244         P(I,2)=-P(IPU3,2)-P(I-1,2)
0245         PMS(4)=P(I,5)**2+P(I,1)**2+P(I,2)**2
0246
0247 C...Relative distribution of energy for particle into jet plus particle.
0248         IMB=1
0249         IF(MOD(MINT(10+ISIDE)/1000,10).NE.0) IMB=2
0250         IF(MSTP(94).LE.1) THEN
0251           IF(IMB.EQ.1) CHI=PYR(0)
0252           IF(IMB.EQ.2) CHI=1D0-SQRT(PYR(0))
0253           IF(MOD(KFLCH/1000,10).NE.0) CHI=1D0-CHI
0254         ELSEIF(MSTP(94).EQ.2) THEN
0255           CHI=1D0-PYR(0)**(1D0/(1D0+PARP(93+2*IMB)))
0256           IF(MOD(KFLCH/1000,10).NE.0) CHI=1D0-CHI
0257         ELSEIF(MSTP(94).EQ.3) THEN
0258           CALL PYZDIS(1,0,PMS(4),ZZ)
0259           CHI=ZZ
0260         ELSE
0261           CALL PYZDIS(1000,0,PMS(4),ZZ)
0262           CHI=ZZ
0263         ENDIF
0264
0265 C...Construct total transverse mass; reject if too large.
0266         CHI=MAX(1D-8,MIN(1D0-1D-8,CHI))
0267         PMS(ISIDE)=PMS(4)/CHI+PMS(3)/(1D0-CHI)
0268         IF(PMS(ISIDE).GT.P(IPU3,4)**2) THEN
0269           IF(LOOP.LT.10) GOTO 200
0270           GOTO 180
0271         ENDIF
0272         VINT(158+ISIDE)=CHI
0273
0274 C...Subdivide longitudinal momentum according to value selected above.
0275         PRP=SQRT(PMS(ISIDE)+P(IPU3,3)**2)+ABS(P(IPU3,3))
0276         PW1=(1D0-CHI)*PRP
0277         P(I-1,4)=0.5D0*(PW1+PMS(3)/PW1)
0278         P(I-1,3)=0.5D0*(PW1-PMS(3)/PW1)*SIDESG
0279         PW2=CHI*PRP
0280         P(I,4)=0.5D0*(PW2+PMS(4)/PW2)
0281         P(I,3)=0.5D0*(PW2-PMS(4)/PW2)*SIDESG
0282       ENDIF
0283       N=I
0284
0285 C...Boost current and remnant systems to correct frame.
0286       IF(SQRT(PMS(1))+SQRT(PMS(2)).GT.0.99D0*VINT(1)) GOTO 180
0287       DSQLAM=SQRT(MAX(0D0,(VINT(2)-PMS(1)-PMS(2))**2-4D0*PMS(1)*PMS(2)))
0288       DRKC=(VINT(2)+PMS(3-ISIDE)-PMS(ISIDE)+DSQLAM)/
0289      &(2D0*VINT(1)*PCP)
0290       DRKR=(VINT(2)+PMS(ISIDE)-PMS(3-ISIDE)+DSQLAM)/
0291      &(2D0*VINT(1)*PRP)
0292       DBEC=-SIDESG*(DRKC**2-1D0)/(DRKC**2+1D0)
0293       DBER=SIDESG*(DRKR**2-1D0)/(DRKR**2+1D0)
0294       CALL PYROBO(IPU3,IPU3,0D0,0D0,0D0,0D0,DBEC)
0295       CALL PYROBO(IPU3+1,N,0D0,0D0,0D0,0D0,DBER)
0296
0297 C...Let current quark shower; recoil but no showering by colour partner.
0298       QMAX=2D0*SQRT(VINT(309-ISIDE))
0299       MSTJ48=MSTJ(48)
0300       MSTJ(48)=1
0301       PARJ86=PARJ(86)
0302       PARJ(86)=0D0
0303       IF(MSTP(71).EQ.1) CALL PYSHOW(IPU3,ICOLR,QMAX)
0304       MSTJ(48)=MSTJ48
0305       PARJ(86)=PARJ86
0306
0307       RETURN
0308       END