+++ /dev/null
-
-C*********************************************************************
-
- SUBROUTINE PYRAND
-
-C...Generates quantities characterizing the high-pT scattering at the
-C...parton level according to the matrix elements. Chooses incoming,
-C...reacting partons, their momentum fractions and one of the possible
-C...subprocesses.
- COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
- COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
- COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
- COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
- COMMON/PYINT1/MINT(400),VINT(400)
- COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
- COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
- COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
- COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
- COMMON/PYINT7/SIGT(0:6,0:6,0:5)
- COMMON/PYINT9/DXSEC(0:200)
- DOUBLE PRECISION DXSEC
- COMMON/PYUPPR/NUP,KUP(20,7),PUP(20,5),NFUP,IFUP(10,2),Q2UP(0:10)
- SAVE /LUDAT1/,/LUDAT2/
- SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/,/PYINT4/,
- &/PYINT5/,/PYINT7/,/PYINT9/,/PYUPPR/
- DIMENSION XPQ(-25:25),PMM(2),PDIF(4),BHAD(4)
-
-C...Parameters and data used in elastic/diffractive treatment.
- DATA EPS/0.0808/, ALP/0.25/, CRES/2./, PMRC/1.062/, SMP/0.880/
- DATA BHAD/2.3,1.4,1.4,0.23/
-
-C...Initial values, specifically for (first) semihard interaction.
- MINT(10)=0
- MINT(17)=0
- MINT(18)=0
- VINT(143)=1.
- VINT(144)=1.
- MFAIL=0
- IF(MSTP(171).EQ.1.AND.MSTP(172).EQ.2) MFAIL=1
- ISUB=0
- LOOP=0
- 100 LOOP=LOOP+1
- MINT(51)=0
-
-C...Choice of process type - first event of pileup.
- IF(MINT(82).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96)) THEN
-
-C...For gamma-p or gamma-gamma first pick between alternatives.
- IF(MINT(121).GT.1) CALL PYSAVE(4,IGA)
- MINT(122)=IGA
-
-C...For gamma + gamma with different nature, flip at random.
- IF(MINT(11).EQ.22.AND.MINT(12).EQ.22.AND.MINT(123).GE.4.AND.
- & RLU(0).GT.0.5) THEN
- MINTSV=MINT(41)
- MINT(41)=MINT(42)
- MINT(42)=MINTSV
- MINTSV=MINT(45)
- MINT(45)=MINT(46)
- MINT(46)=MINTSV
- MINTSV=MINT(107)
- MINT(107)=MINT(108)
- MINT(108)=MINTSV
- IF(MINT(47).EQ.2.OR.MINT(47).EQ.3) MINT(47)=5-MINT(47)
- ENDIF
-
-C...Pick process type.
- RSUB=XSEC(0,1)*RLU(0)
- DO 110 I=1,200
- IF(MSUB(I).NE.1) GOTO 110
- ISUB=I
- RSUB=RSUB-XSEC(I,1)
- IF(RSUB.LE.0.) GOTO 120
- 110 CONTINUE
- 120 IF(ISUB.EQ.95) ISUB=96
- IF(ISUB.EQ.96) CALL PYMULT(2)
-
-C...Choice of inclusive process type - pileup events.
- ELSEIF(MINT(82).GE.2.AND.ISUB.EQ.0) THEN
- RSUB=VINT(131)*RLU(0)
- ISUB=96
- IF(RSUB.GT.SIGT(0,0,5)) ISUB=94
- IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)) ISUB=93
- IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)+SIGT(0,0,3)) ISUB=92
- IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)+SIGT(0,0,3)+SIGT(0,0,2))
- & ISUB=91
- IF(ISUB.EQ.96) CALL PYMULT(2)
- ENDIF
- IF(MINT(82).EQ.1) NGEN(0,1)=NGEN(0,1)+1
- IF(MINT(82).EQ.1) NGEN(ISUB,1)=NGEN(ISUB,1)+1
- IF(ISUB.EQ.96.AND.LOOP.EQ.1.AND.MINT(82).EQ.1)
- &NGEN(97,1)=NGEN(97,1)+1
- MINT(1)=ISUB
- ISTSB=ISET(ISUB)
-
-C...Find resonances (explicit or implicit in cross-section).
- MINT(72)=0
- KFR1=0
- IF(ISTSB.EQ.1.OR.ISTSB.EQ.3.OR.ISTSB.EQ.5) THEN
- KFR1=KFPR(ISUB,1)
- ELSEIF(ISUB.EQ.24.OR.ISUB.EQ.25.OR.ISUB.EQ.110.OR.ISUB.EQ.165.OR.
- &ISUB.EQ.171.OR.ISUB.EQ.176) THEN
- KFR1=23
- ELSEIF(ISUB.EQ.23.OR.ISUB.EQ.26.OR.ISUB.EQ.166.OR.ISUB.EQ.172.OR.
- &ISUB.EQ.177) THEN
- KFR1=24
- ELSEIF(ISUB.GE.71.AND.ISUB.LE.77) THEN
- KFR1=25
- IF(MSTP(46).EQ.5) THEN
- KFR1=30
- PMAS(30,1)=PARP(45)
- PMAS(30,2)=PARP(45)**3/(96.*PARU(1)*PARP(47)**2)
- ENDIF
- ENDIF
- CKMX=CKIN(2)
- IF(CKMX.LE.0.) CKMX=VINT(1)
- IF(KFR1.NE.0) THEN
- IF(CKIN(1).GT.PMAS(KFR1,1)+20.*PMAS(KFR1,2).OR.
- & CKMX.LT.PMAS(KFR1,1)-20.*PMAS(KFR1,2)) KFR1=0
- ENDIF
- IF(KFR1.NE.0) THEN
- TAUR1=PMAS(KFR1,1)**2/VINT(2)
- GAMR1=PMAS(KFR1,1)*PMAS(KFR1,2)/VINT(2)
- MINT(72)=1
- MINT(73)=KFR1
- VINT(73)=TAUR1
- VINT(74)=GAMR1
- ENDIF
- IF(ISUB.EQ.141) THEN
- KFR2=23
- TAUR2=PMAS(KFR2,1)**2/VINT(2)
- GAMR2=PMAS(KFR2,1)*PMAS(KFR2,2)/VINT(2)
- IF(CKIN(1).GT.PMAS(KFR2,1)+20.*PMAS(KFR2,2).OR.
- & CKMX.LT.PMAS(KFR2,1)-20.*PMAS(KFR2,2)) KFR2=0
- IF(KFR2.NE.0.AND.KFR1.NE.0) THEN
- MINT(72)=2
- MINT(74)=KFR2
- VINT(75)=TAUR2
- VINT(76)=GAMR2
- ELSEIF(KFR2.NE.0) THEN
- KFR1=KFR2
- TAUR1=TAUR2
- GAMR1=GAMR2
- MINT(72)=1
- MINT(73)=KFR1
- VINT(73)=TAUR1
- VINT(74)=GAMR1
- ENDIF
- ENDIF
-
-C...Find product masses and minimum pT of process,
-C...optionally with broadening according to a truncated Breit-Wigner.
- VINT(63)=0.
- VINT(64)=0.
- MINT(71)=0
- VINT(71)=CKIN(3)
- IF(MINT(82).GE.2) VINT(71)=0.
- VINT(80)=1.
- IF(ISTSB.EQ.2.OR.ISTSB.EQ.4) THEN
- NBW=0
- DO 130 I=1,2
- IF(KFPR(ISUB,I).EQ.0) THEN
- ELSEIF(MSTP(42).LE.0.OR.PMAS(LUCOMP(KFPR(ISUB,I)),2).LT.
- & PARP(41)) THEN
- VINT(62+I)=PMAS(LUCOMP(KFPR(ISUB,I)),1)**2
- ELSE
- NBW=NBW+1
- ENDIF
- 130 CONTINUE
- IF(NBW.GE.1) THEN
- CALL PYOFSH(4,0,KFPR(ISUB,1),KFPR(ISUB,2),0.,PQM3,PQM4)
- IF(MINT(51).EQ.1) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- IF(MFAIL.EQ.1) THEN
- MSTI(61)=1
- RETURN
- ENDIF
- GOTO 100
- ENDIF
- VINT(63)=PQM3**2
- VINT(64)=PQM4**2
- ENDIF
- IF(MIN(VINT(63),VINT(64)).LT.CKIN(6)**2) MINT(71)=1
- IF(MINT(71).EQ.1) VINT(71)=MAX(CKIN(3),CKIN(5))
- ELSEIF(ISTSB.EQ.6) THEN
- CALL PYOFSH(6,0,KFPR(ISUB,1),KFPR(ISUB,2),0.,PQM3,PQM4)
- IF(MINT(51).EQ.1) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- IF(MFAIL.EQ.1) THEN
- MSTI(61)=1
- RETURN
- ENDIF
- GOTO 100
- ENDIF
- VINT(63)=PQM3**2
- VINT(64)=PQM4**2
- ENDIF
-
-C...Prepare for additional variable choices in 2 -> 3.
- IF(ISTSB.EQ.5) THEN
- VINT(201)=0.
- IF(KFPR(ISUB,2).GT.0) VINT(201)=PMAS(KFPR(ISUB,2),1)
- VINT(206)=VINT(201)
- VINT(204)=PMAS(23,1)
- IF(ISUB.EQ.124) VINT(204)=PMAS(24,1)
- IF(ISUB.EQ.121.OR.ISUB.EQ.122.OR.ISUB.EQ.181.OR.ISUB.EQ.182.OR.
- & ISUB.EQ.186.OR.ISUB.EQ.187) VINT(204)=VINT(201)
- VINT(209)=VINT(204)
- ENDIF
-
-C...Select incoming VDM particle (rho/omega/phi/J/psi).
- IF(ISTSB.NE.0.AND.(MINT(101).GE.2.OR.MINT(102).GE.2).AND.
- &(MINT(123).EQ.2.OR.MINT(123).EQ.5.OR.MINT(123).EQ.7)) THEN
- VRN=RLU(0)*SIGT(0,0,5)
- IF(MINT(101).LE.1) THEN
- I1MN=0
- I1MX=0
- ELSE
- I1MN=1
- I1MX=MINT(101)
- ENDIF
- IF(MINT(102).LE.1) THEN
- I2MN=0
- I2MX=0
- ELSE
- I2MN=1
- I2MX=MINT(102)
- ENDIF
- DO 150 I1=I1MN,I1MX
- KFV1=110*I1+3
- DO 140 I2=I2MN,I2MX
- KFV2=110*I2+3
- VRN=VRN-SIGT(I1,I2,5)
- IF(VRN.LE.0.) GOTO 160
- 140 CONTINUE
- 150 CONTINUE
- 160 IF(MINT(101).GE.2) MINT(103)=KFV1
- IF(MINT(102).GE.2) MINT(104)=KFV2
- ENDIF
-
- IF(ISTSB.EQ.0) THEN
-C...Elastic scattering or single or double diffractive scattering.
-
-C...Select incoming particle (rho/omega/phi/J/psi for VDM) and mass.
- MINT(103)=MINT(11)
- MINT(104)=MINT(12)
- PMM(1)=VINT(3)
- PMM(2)=VINT(4)
- IF(MINT(101).GE.2.OR.MINT(102).GE.2) THEN
- JJ=ISUB-90
- VRN=RLU(0)*SIGT(0,0,JJ)
- IF(MINT(101).LE.1) THEN
- I1MN=0
- I1MX=0
- ELSE
- I1MN=1
- I1MX=MINT(101)
- ENDIF
- IF(MINT(102).LE.1) THEN
- I2MN=0
- I2MX=0
- ELSE
- I2MN=1
- I2MX=MINT(102)
- ENDIF
- DO 180 I1=I1MN,I1MX
- KFV1=110*I1+3
- DO 170 I2=I2MN,I2MX
- KFV2=110*I2+3
- VRN=VRN-SIGT(I1,I2,JJ)
- IF(VRN.LE.0.) GOTO 190
- 170 CONTINUE
- 180 CONTINUE
- 190 IF(MINT(101).GE.2) THEN
- MINT(103)=KFV1
- PMM(1)=ULMASS(KFV1)
- ENDIF
- IF(MINT(102).GE.2) THEN
- MINT(104)=KFV2
- PMM(2)=ULMASS(KFV2)
- ENDIF
- ENDIF
-
-C...Side/sides of diffractive system.
- MINT(17)=0
- MINT(18)=0
- IF(ISUB.EQ.92.OR.ISUB.EQ.94) MINT(17)=1
- IF(ISUB.EQ.93.OR.ISUB.EQ.94) MINT(18)=1
-
-C...Find masses of particles and minimal masses of diffractive states.
- DO 200 JT=1,2
- PDIF(JT)=PMM(JT)
- VINT(66+JT)=PDIF(JT)
- IF(MINT(16+JT).EQ.1) PDIF(JT)=PDIF(JT)+PARP(102)
- 200 CONTINUE
- SH=VINT(2)
- SQM1=PMM(1)**2
- SQM2=PMM(2)**2
- SQM3=PDIF(1)**2
- SQM4=PDIF(2)**2
- SMRES1=(PMM(1)+PMRC)**2
- SMRES2=(PMM(2)+PMRC)**2
-
-C...Find elastic slope and lower limit diffractive slope.
- IHA=MAX(2,IABS(MINT(103))/110)
- IF(IHA.GE.5) IHA=1
- IHB=MAX(2,IABS(MINT(104))/110)
- IF(IHB.GE.5) IHB=1
- IF(ISUB.EQ.91) THEN
- BMN=2.*BHAD(IHA)+2.*BHAD(IHB)+4.*SH**EPS-4.2
- ELSEIF(ISUB.EQ.92) THEN
- BMN=MAX(2.,2.*BHAD(IHB))
- ELSEIF(ISUB.EQ.93) THEN
- BMN=MAX(2.,2.*BHAD(IHA))
- ELSEIF(ISUB.EQ.94) THEN
- BMN=2.*ALP*4.
- ENDIF
-
-C...Determine maximum possible t range and coefficient of generation.
- SQLA12=(SH-SQM1-SQM2)**2-4.*SQM1*SQM2
- SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4
- THA=SH-(SQM1+SQM2+SQM3+SQM4)+(SQM1-SQM2)*(SQM3-SQM4)/SH
- THB=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH
- THC=(SQM3-SQM1)*(SQM4-SQM2)+(SQM1+SQM4-SQM2-SQM3)*
- & (SQM1*SQM4-SQM2*SQM3)/SH
- THL=-0.5*(THA+THB)
- THU=THC/THL
- THRND=EXP(MAX(-50.,BMN*(THL-THU)))-1.
-
-C...Select diffractive mass/masses according to dm^2/m^2.
- 210 DO 220 JT=1,2
- IF(MINT(16+JT).EQ.0) THEN
- PDIF(2+JT)=PDIF(JT)
- ELSE
- PMMIN=PDIF(JT)
- PMMAX=MAX(VINT(2+JT),VINT(1)-PDIF(3-JT))
- PDIF(2+JT)=PMMIN*(PMMAX/PMMIN)**RLU(0)
- ENDIF
- 220 CONTINUE
- SQM3=PDIF(3)**2
- SQM4=PDIF(4)**2
-
-C..Additional mass factors, including resonance enhancement.
- IF(PDIF(3)+PDIF(4).GE.VINT(1)) GOTO 210
- IF(ISUB.EQ.92) THEN
- FSD=(1.-SQM3/SH)*(1.+CRES*SMRES1/(SMRES1+SQM3))
- IF(FSD.LT.RLU(0)*(1.+CRES)) GOTO 210
- ELSEIF(ISUB.EQ.93) THEN
- FSD=(1.-SQM4/SH)*(1.+CRES*SMRES2/(SMRES2+SQM4))
- IF(FSD.LT.RLU(0)*(1.+CRES)) GOTO 210
- ELSEIF(ISUB.EQ.94) THEN
- FDD=(1.-(PDIF(3)+PDIF(4))**2/SH)*(SH*SMP/(SH*SMP+SQM3*SQM4))*
- & (1.+CRES*SMRES1/(SMRES1+SQM3))*(1.+CRES*SMRES2/(SMRES2+SQM4))
- IF(FDD.LT.RLU(0)*(1.+CRES)**2) GOTO 210
- ENDIF
-
-C...Select t according to exp(Bmn*t) and correct to right slope.
- TH=THU+LOG(1.+THRND*RLU(0))/BMN
- IF(ISUB.GE.92) THEN
- IF(ISUB.EQ.92) THEN
- BADD=2.*ALP*LOG(SH/SQM3)
- IF(BHAD(IHB).LT.1.) BADD=MAX(0.,BADD+2.*BHAD(IHB)-2.)
- ELSEIF(ISUB.EQ.93) THEN
- BADD=2.*ALP*LOG(SH/SQM4)
- IF(BHAD(IHA).LT.1.) BADD=MAX(0.,BADD+2.*BHAD(IHA)-2.)
- ELSEIF(ISUB.EQ.94) THEN
- BADD=2.*ALP*(LOG(EXP(4.)+SH/(ALP*SQM3*SQM4))-4.)
- ENDIF
- IF(EXP(MAX(-50.,BADD*(TH-THU))).LT.RLU(0)) GOTO 210
- ENDIF
-
-C...Check whether m^2 and t choices are consistent.
- SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4
- THA=SH-(SQM1+SQM2+SQM3+SQM4)+(SQM1-SQM2)*(SQM3-SQM4)/SH
- THB=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH
- IF(THB.LE.1E-8) GOTO 210
- THC=(SQM3-SQM1)*(SQM4-SQM2)+(SQM1+SQM4-SQM2-SQM3)*
- & (SQM1*SQM4-SQM2*SQM3)/SH
- THLM=-0.5*(THA+THB)
- THUM=THC/THLM
- IF(TH.LT.THLM.OR.TH.GT.THUM) GOTO 210
-
-C...Information to output.
- VINT(21)=1.
- VINT(22)=0.
- VINT(23)=MIN(1.,MAX(-1.,(THA+2.*TH)/THB))
- VINT(45)=TH
- VINT(59)=2.*SQRT(MAX(0.,-(THC+THA*TH+TH**2)))/THB
- VINT(63)=PDIF(3)**2
- VINT(64)=PDIF(4)**2
-
-C...Note: in the following, by In is meant the integral over the
-C...quantity multiplying coefficient cn.
-C...Choose tau according to h1(tau)/tau, where
-C...h1(tau) = c1 + I1/I2*c2*1/tau + I1/I3*c3*1/(tau+tau_R) +
-C...I1/I4*c4*tau/((s*tau-m^2)^2+(m*Gamma)^2) +
-C...I1/I5*c5*1/(tau+tau_R') +
-C...I1/I6*c6*tau/((s*tau-m'^2)^2+(m'*Gamma')^2) +
-C...I1/I7*c7*tau/(1.-tau), and
-C...c1 + c2 + c3 + c4 + c5 + c6 + c7 = 1.
- ELSEIF(ISTSB.GE.1.AND.ISTSB.LE.6) THEN
- CALL PYKLIM(1)
- IF(MINT(51).NE.0) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- IF(MFAIL.EQ.1) THEN
- MSTI(61)=1
- RETURN
- ENDIF
- GOTO 100
- ENDIF
- RTAU=RLU(0)
- MTAU=1
- IF(RTAU.GT.COEF(ISUB,1)) MTAU=2
- IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)) MTAU=3
- IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)) MTAU=4
- IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4))
- & MTAU=5
- IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
- & COEF(ISUB,5)) MTAU=6
- IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
- & COEF(ISUB,5)+COEF(ISUB,6)) MTAU=7
- CALL PYKMAP(1,MTAU,RLU(0))
-
-C...2 -> 3, 4 processes:
-C...Choose tau' according to h4(tau,tau')/tau', where
-C...h4(tau,tau') = c1 + I1/I2*c2*(1 - tau/tau')^3/tau' +
-C...I1/I3*c3*1/(1 - tau'), and c1 + c2 + c3 = 1.
- IF(ISTSB.GE.3.AND.ISTSB.LE.5) THEN
- CALL PYKLIM(4)
- IF(MINT(51).NE.0) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- IF(MFAIL.EQ.1) THEN
- MSTI(61)=1
- RETURN
- ENDIF
- GOTO 100
- ENDIF
- RTAUP=RLU(0)
- MTAUP=1
- IF(RTAUP.GT.COEF(ISUB,18)) MTAUP=2
- IF(RTAUP.GT.COEF(ISUB,18)+COEF(ISUB,19)) MTAUP=3
- CALL PYKMAP(4,MTAUP,RLU(0))
- ENDIF
-
-C...Choose y* according to h2(y*), where
-C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) +
-C...I0/I3*c3*1/cosh(y*) + I0/I4*c4*1/(1-exp(y*-y*max)) +
-C...I0/I5*c5*1/(1-exp(-y*-y*min)), I0 = y*max-y*min,
-C...and c1 + c2 + c3 + c4 + c5 = 1.
- CALL PYKLIM(2)
- IF(MINT(51).NE.0) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- IF(MFAIL.EQ.1) THEN
- MSTI(61)=1
- RETURN
- ENDIF
- GOTO 100
- ENDIF
- RYST=RLU(0)
- MYST=1
- IF(RYST.GT.COEF(ISUB,8)) MYST=2
- IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
- IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)+COEF(ISUB,10)) MYST=4
- IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)+COEF(ISUB,10)+
- & COEF(ISUB,11)) MYST=5
- CALL PYKMAP(2,MYST,RLU(0))
-
-C...2 -> 2 processes:
-C...Choose cos(theta-hat) (cth) according to h3(cth), where
-C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) +
-C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2,
-C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products),
-C...and c0 + c1 + c2 + c3 + c4 = 1.
- CALL PYKLIM(3)
- IF(MINT(51).NE.0) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- IF(MFAIL.EQ.1) THEN
- MSTI(61)=1
- RETURN
- ENDIF
- GOTO 100
- ENDIF
- IF(ISTSB.EQ.2.OR.ISTSB.EQ.4.OR.ISTSB.EQ.6) THEN
- RCTH=RLU(0)
- MCTH=1
- IF(RCTH.GT.COEF(ISUB,13)) MCTH=2
- IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)) MCTH=3
- IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)+COEF(ISUB,15)) MCTH=4
- IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)+COEF(ISUB,15)+
- & COEF(ISUB,16)) MCTH=5
- CALL PYKMAP(3,MCTH,RLU(0))
- ENDIF
-
-C...2 -> 3 : select pT1, phi1, pT2, phi2, y3 for 3 outgoing.
- IF(ISTSB.EQ.5) THEN
- CALL PYKMAP(5,0,0.)
- IF(MINT(51).NE.0) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- IF(MFAIL.EQ.1) THEN
- MSTI(61)=1
- RETURN
- ENDIF
- GOTO 100
- ENDIF
- ENDIF
-
-C...Low-pT or multiple interactions (first semihard interaction).
- ELSEIF(ISTSB.EQ.9) THEN
- CALL PYMULT(3)
- ISUB=MINT(1)
-
-C...Generate user-defined process: kinematics plus weight.
- ELSEIF(ISTSB.EQ.11) THEN
- MSTI(51)=0
- CALL PYUPEV(ISUB,SIGS)
- IF(NUP.LE.0) THEN
- MINT(51)=2
- MSTI(51)=1
- IF(MINT(82).EQ.1) THEN
- NGEN(0,1)=NGEN(0,1)-1
- NGEN(0,2)=NGEN(0,2)-1
- NGEN(ISUB,1)=NGEN(ISUB,1)-1
- ENDIF
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- RETURN
- ENDIF
-
-C...Construct 'trivial' kinematical variables needed.
- KFL1=KUP(1,2)
- KFL2=KUP(2,2)
- VINT(41)=2.*PUP(1,4)/VINT(1)
- VINT(42)=2.*PUP(2,4)/VINT(1)
- VINT(21)=VINT(41)*VINT(42)
- VINT(22)=0.5*LOG(VINT(41)/VINT(42))
- VINT(44)=VINT(21)*VINT(2)
- VINT(43)=SQRT(MAX(0.,VINT(44)))
- VINT(56)=Q2UP(0)
- VINT(55)=SQRT(MAX(0.,VINT(56)))
-
-C...Construct other kinematical variables needed (approximately).
- VINT(23)=0.
- VINT(26)=VINT(21)
- VINT(45)=-0.5*VINT(44)
- VINT(46)=-0.5*VINT(44)
- VINT(49)=VINT(43)
- VINT(50)=VINT(44)
- VINT(51)=VINT(55)
- VINT(52)=VINT(56)
- VINT(53)=VINT(55)
- VINT(54)=VINT(56)
- VINT(25)=0.
- VINT(48)=0.
- DO 230 IUP=3,NUP
- IF(KUP(IUP,1).EQ.1) VINT(25)=VINT(25)+2.*(PUP(IUP,5)**2+
- & PUP(IUP,1)**2+PUP(IUP,2)**2)/VINT(1)
- IF(KUP(IUP,1).EQ.1) VINT(48)=VINT(48)+0.5*(PUP(IUP,1)**2+
- & PUP(IUP,2)**2)
- 230 CONTINUE
- VINT(47)=SQRT(VINT(48))
-
-C...Calculate structure function weights.
- IF(MINT(47).GE.2) THEN
- DO 250 I=3-MIN(2,MINT(45)),MIN(2,MINT(46))
- MINT(105)=MINT(102+I)
- MINT(109)=MINT(106+I)
- IF(MSTP(57).LE.1) THEN
- CALL PYSTFU(MINT(10+I),VINT(40+I),Q2UP(0),XPQ)
- ELSE
- CALL PYSTFL(MINT(10+I),VINT(40+I),Q2UP(0),XPQ)
- ENDIF
- DO 240 KFL=-25,25
- XSFX(I,KFL)=XPQ(KFL)
- 240 CONTINUE
- 250 CONTINUE
- ENDIF
- ENDIF
-
-C...Choose azimuthal angle.
- VINT(24)=PARU(2)*RLU(0)
-
-C...Check against user cuts on kinematics at parton level.
- MINT(51)=0
- IF((ISUB.LE.90.OR.ISUB.GT.100).AND.ISTSB.LE.10) CALL PYKLIM(0)
- IF(MINT(51).NE.0) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- IF(MFAIL.EQ.1) THEN
- MSTI(61)=1
- RETURN
- ENDIF
- GOTO 100
- ENDIF
- IF(MINT(82).EQ.1.AND.MSTP(141).GE.1.AND.ISTSB.LE.10) THEN
- MCUT=0
- IF(MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+MSUB(95).EQ.0)
- & CALL PYKCUT(MCUT)
- IF(MCUT.NE.0) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- IF(MFAIL.EQ.1) THEN
- MSTI(61)=1
- RETURN
- ENDIF
- GOTO 100
- ENDIF
- ENDIF
-
-C...Calculate differential cross-section for different subprocesses.
- IF(ISTSB.LE.10) CALL PYSIGH(NCHN,SIGS)
- SIGSOR=SIGS
- SIGLPT=SIGT(0,0,5)
-
-C...Multiply cross-section by user-defined weights.
- IF(MSTP(173).EQ.1) THEN
- SIGS=PARP(173)*SIGS
- DO 260 ICHN=1,NCHN
- SIGH(ICHN)=PARP(173)*SIGH(ICHN)
- 260 CONTINUE
- SIGLPT=PARP(173)*SIGLPT
- ENDIF
- WTXS=1.
- SIGSWT=SIGS
- VINT(99)=1.
- VINT(100)=1.
- IF(MINT(82).EQ.1.AND.MSTP(142).GE.1) THEN
- IF(ISUB.NE.96.AND.MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+
- & MSUB(95).EQ.0) CALL PYEVWT(WTXS)
- SIGSWT=WTXS*SIGS
- VINT(99)=WTXS
- IF(MSTP(142).EQ.1) VINT(100)=1./WTXS
- ENDIF
-
-C...Calculations for Monte Carlo estimate of all cross-sections.
- IF(MINT(82).EQ.1.AND.ISUB.LE.90.OR.ISUB.GE.96) THEN
- IF(MSTP(142).LE.1) THEN
- XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
- DXSEC(ISUB)=DXSEC(ISUB)+SIGS
- ELSE
- XSEC(ISUB,2)=XSEC(ISUB,2)+SIGSWT
- DXSEC(ISUB)=DXSEC(ISUB)+SIGSWT
- ENDIF
- ELSEIF(MINT(82).EQ.1) THEN
- XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
- DXSEC(ISUB)=DXSEC(ISUB)+SIGS
- ENDIF
- IF((ISUB.EQ.95.OR.ISUB.EQ.96).AND.LOOP.EQ.1.AND.MINT(82).EQ.1)
- &THEN
- XSEC(97,2)=XSEC(97,2)+SIGLPT
- DXSEC(97)=DXSEC(97)+SIGLPT
- ENDIF
-
-C...Multiple interactions: store results of cross-section calculation.
- IF(MINT(50).EQ.1.AND.MSTP(82).GE.3) THEN
- VINT(153)=SIGSOR
- CALL PYMULT(4)
- ENDIF
-
-C...Check that weight not negative.
- VIOL=SIGSWT/XSEC(ISUB,1)
- IF(ISUB.EQ.96.AND.MSTP(173).EQ.1) VIOL=VIOL/PARP(174)
- IF(MSTP(123).LE.0) THEN
- IF(VIOL.LT.-1E-3) THEN
- WRITE(MSTU(11),5000) VIOL,NGEN(0,3)+1
- WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
- STOP
- ENDIF
- ELSE
- IF(VIOL.LT.MIN(-1E-3,VINT(109))) THEN
- VINT(109)=VIOL
- WRITE(MSTU(11),5200) VIOL,NGEN(0,3)+1
- WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
- ENDIF
- ENDIF
-
-C...Weighting using estimate of maximum of differential cross-section.
- IF(MFAIL.EQ.0) THEN
- IF(VIOL.LT.RLU(0)) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- GOTO 100
- ENDIF
- ELSEIF(ISUB.NE.95.AND.ISUB.NE.96) THEN
- IF(VIOL.LT.RLU(0)) THEN
- MSTI(61)=1
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- RETURN
- ENDIF
- ELSE
- RATND=SIGLPT/XSEC(95,1)
- IF(LOOP.EQ.1.AND.RATND.LT.RLU(0)) THEN
- MSTI(61)=1
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- RETURN
- ENDIF
- VIOL=VIOL/RATND
- IF(VIOL.LT.RLU(0)) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- GOTO 100
- ENDIF
- ENDIF
-
-C...Check for possible violation of estimated maximum of differential
-C...cross-section used in weighting.
- IF(MSTP(123).LE.0) THEN
- IF(VIOL.GT.1.) THEN
- WRITE(MSTU(11),5300) VIOL,NGEN(0,3)+1
- WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
- STOP
- ENDIF
- ELSEIF(MSTP(123).EQ.1) THEN
- IF(VIOL.GT.VINT(108)) THEN
- VINT(108)=VIOL
- IF(VIOL.GT.1.) THEN
- MINT(10)=1
- WRITE(MSTU(11),5400) VIOL,NGEN(0,3)+1
- WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),
- & VINT(26)
- ENDIF
- ENDIF
- ELSEIF(VIOL.GT.VINT(108)) THEN
- VINT(108)=VIOL
- IF(VIOL.GT.1.) THEN
- MINT(10)=1
- XDIF=XSEC(ISUB,1)*(VIOL-1.)
- XSEC(ISUB,1)=XSEC(ISUB,1)+XDIF
- IF(MSUB(ISUB).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96))
- & XSEC(0,1)=XSEC(0,1)+XDIF
- WRITE(MSTU(11),5400) VIOL,NGEN(0,3)+1
- WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
- IF(ISUB.LE.9) THEN
- WRITE(MSTU(11),5500) ISUB,XSEC(ISUB,1)
- ELSEIF(ISUB.LE.99) THEN
- WRITE(MSTU(11),5600) ISUB,XSEC(ISUB,1)
- ELSE
- WRITE(MSTU(11),5700) ISUB,XSEC(ISUB,1)
- ENDIF
- VINT(108)=1.
- ENDIF
- ENDIF
-
-C...Multiple interactions: choose impact parameter.
- VINT(148)=1.
- IF(MINT(50).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GE.96).AND.MSTP(82).GE.3)
- &THEN
- CALL PYMULT(5)
- IF(VINT(150).LT.RLU(0)) THEN
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
- IF(MFAIL.EQ.1) THEN
- MSTI(61)=1
- RETURN
- ENDIF
- GOTO 100
- ENDIF
- ENDIF
- IF(MINT(82).EQ.1) NGEN(0,2)=NGEN(0,2)+1
- IF(MINT(82).EQ.1.AND.MSUB(95).EQ.1) THEN
- IF(ISUB.LE.90.OR.ISUB.GE.95) NGEN(95,1)=NGEN(95,1)+1
- IF(ISUB.LE.90.OR.ISUB.GE.96) NGEN(96,2)=NGEN(96,2)+1
- ENDIF
- IF(ISUB.LE.90.OR.ISUB.GE.96) MINT(31)=MINT(31)+1
-
-C...Choose flavour of reacting partons (and subprocess).
- IF(ISTSB.GE.11) GOTO 280
- RSIGS=SIGS*RLU(0)
- QT2=VINT(48)
- RQQBAR=PARP(87)*(1.-(QT2/(QT2+(PARP(88)*PARP(82))**2))**2)
- IF(ISUB.NE.95.AND.(ISUB.NE.96.OR.MSTP(82).LE.1.OR.
- &RLU(0).GT.RQQBAR)) THEN
- DO 270 ICHN=1,NCHN
- KFL1=ISIG(ICHN,1)
- KFL2=ISIG(ICHN,2)
- MINT(2)=ISIG(ICHN,3)
- RSIGS=RSIGS-SIGH(ICHN)
- IF(RSIGS.LE.0.) GOTO 280
- 270 CONTINUE
-
-C...Multiple interactions: choose qq~ preferentially at small pT.
- ELSEIF(ISUB.EQ.96) THEN
- MINT(105)=MINT(103)
- MINT(109)=MINT(107)
- CALL PYSPLI(MINT(11),21,KFL1,KFLDUM)
- MINT(105)=MINT(104)
- MINT(109)=MINT(108)
- CALL PYSPLI(MINT(12),21,KFL2,KFLDUM)
- MINT(1)=11
- MINT(2)=1
- IF(KFL1.EQ.KFL2.AND.RLU(0).LT.0.5) MINT(2)=2
-
-C...Low-pT: choose string drawing configuration.
- ELSE
- KFL1=21
- KFL2=21
- RSIGS=6.*RLU(0)
- MINT(2)=1
- IF(RSIGS.GT.1.) MINT(2)=2
- IF(RSIGS.GT.2.) MINT(2)=3
- ENDIF
-
-C...Reassign QCD process. Partons before initial state radiation.
- 280 IF(MINT(2).GT.10) THEN
- MINT(1)=MINT(2)/10
- MINT(2)=MOD(MINT(2),10)
- ENDIF
- IF(MINT(82).EQ.1.AND.MSTP(111).GE.0) NGEN(MINT(1),2)=
- &NGEN(MINT(1),2)+1
- MINT(15)=KFL1
- MINT(16)=KFL2
- MINT(13)=MINT(15)
- MINT(14)=MINT(16)
- VINT(141)=VINT(41)
- VINT(142)=VINT(42)
- VINT(151)=0.
- VINT(152)=0.
-
-C...Calculate x value of photon for parton inside photon inside e.
- DO 310 JT=1,2
- MINT(18+JT)=0
- VINT(154+JT)=0.
- MSPLI=0
- IF(JT.EQ.1.AND.MINT(43).LE.2) MSPLI=1
- IF(JT.EQ.2.AND.MOD(MINT(43),2).EQ.1) MSPLI=1
- IF(IABS(MINT(14+JT)).LE.8.OR.MINT(14+JT).EQ.21) MSPLI=MSPLI+1
- IF(MSPLI.EQ.2) THEN
- KFLH=MINT(14+JT)
- XHRD=VINT(140+JT)
- Q2HRD=VINT(54)
- MINT(105)=MINT(102+JT)
- MINT(109)=MINT(106+JT)
- IF(MSTP(57).LE.1) THEN
- CALL PYSTFU(22,XHRD,Q2HRD,XPQ)
- ELSE
- CALL PYSTFL(22,XHRD,Q2HRD,XPQ)
- ENDIF
- WTMX=4.*XPQ(KFLH)
- IF(MSTP(13).EQ.2) THEN
- Q2PMS=Q2HRD/PMAS(11,1)**2
- WTMX=WTMX*LOG(MAX(2.,Q2PMS*(1.-XHRD)/XHRD**2))
- ENDIF
- 290 XE=XHRD**RLU(0)
- XG=MIN(0.999999,XHRD/XE)
- IF(MSTP(57).LE.1) THEN
- CALL PYSTFU(22,XG,Q2HRD,XPQ)
- ELSE
- CALL PYSTFL(22,XG,Q2HRD,XPQ)
- ENDIF
- WT=(1.+(1.-XE)**2)*XPQ(KFLH)
- IF(MSTP(13).EQ.2) WT=WT*LOG(MAX(2.,Q2PMS*(1.-XE)/XE**2))
- IF(WT.LT.RLU(0)*WTMX) GOTO 290
- MINT(18+JT)=1
- VINT(154+JT)=XE
- DO 300 KFLS=-25,25
- XSFX(JT,KFLS)=XPQ(KFLS)
- 300 CONTINUE
- ENDIF
- 310 CONTINUE
-
-C...Pick scale where photon is resolved.
- IF(MINT(107).EQ.3) VINT(283)=PARP(15)**2*
- &(VINT(54)/PARP(15)**2)**RLU(0)
- IF(MINT(108).EQ.3) VINT(284)=PARP(15)**2*
- &(VINT(54)/PARP(15)**2)**RLU(0)
- IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
-
-C...Format statements for differential cross-section maximum violations.
- 5000 FORMAT(1X,'Error: negative cross-section fraction',1P,E11.3,1X,
- &'in event',1X,I7,'.'/1X,'Execution stopped!')
- 5100 FORMAT(1X,'ISUB = ',I3,'; Point of violation:'/1X,'tau =',1P,
- &E11.3,', y* =',E11.3,', cthe = ',0P,F11.7,', tau'' =',1P,E11.3)
- 5200 FORMAT(1X,'Warning: negative cross-section fraction',1P,E11.3,1X,
- &'in event',1X,I7)
- 5300 FORMAT(1X,'Error: maximum violated by',1P,E11.3,1X,
- &'in event',1X,I7,'.'/1X,'Execution stopped!')
- 5400 FORMAT(1X,'Warning: maximum violated by',1P,E11.3,1X,
- &'in event',1X,I7)
- 5500 FORMAT(1X,'XSEC(',I1,',1) increased to',1P,E11.3)
- 5600 FORMAT(1X,'XSEC(',I2,',1) increased to',1P,E11.3)
- 5700 FORMAT(1X,'XSEC(',I3,',1) increased to',1P,E11.3)
-
- RETURN
- END
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff