pythiaeRHIC/pythia/pykmap.f

0001
0002 C*********************************************************************
0003
0004 C...PYKMAP
0005 C...Maps a uniform distribution into a distribution of a kinematical
0006 C...variable according to one of the possibilities allowed. It is
0007 C...assumed that kinematical limits have been set by a PYKLIM call.
0008
0009       SUBROUTINE PYKMAP(IVAR,MVAR,VVAR)
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...Commonblocks.
0016       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0017       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0018       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0019       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0020       COMMON/PYINT1/MINT(400),VINT(400)
0021       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0022       SAVE /PYDAT1/,/PYDAT2/,/PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/
0023
0024 C...Convert VVAR to tau variable.
0025       ISUB=MINT(1)
0026       ISTSB=ISET(ISUB)
0027       IF(IVAR.EQ.1) THEN
0028         TAUMIN=VINT(11)
0029         TAUMAX=VINT(31)
0030         IF(MVAR.EQ.3.OR.MVAR.EQ.4) THEN
0031           TAURE=VINT(73)
0032           GAMRE=VINT(74)
0033         ELSEIF(MVAR.EQ.5.OR.MVAR.EQ.6) THEN
0034           TAURE=VINT(75)
0035           GAMRE=VINT(76)
0036         ENDIF
0037         IF(MINT(47).EQ.1.AND.(ISTSB.EQ.1.OR.ISTSB.EQ.2)) THEN
0038           TAU=1D0
0039         ELSEIF(MVAR.EQ.1) THEN
0040           TAU=TAUMIN*(TAUMAX/TAUMIN)**VVAR
0041         ELSEIF(MVAR.EQ.2) THEN
0042           TAU=TAUMAX*TAUMIN/(TAUMIN+(TAUMAX-TAUMIN)*VVAR)
0043         ELSEIF(MVAR.EQ.3.OR.MVAR.EQ.5) THEN
0044           RATGEN=(TAURE+TAUMAX)/(TAURE+TAUMIN)*TAUMIN/TAUMAX
0045           TAU=TAURE*TAUMIN/((TAURE+TAUMIN)*RATGEN**VVAR-TAUMIN)
0046         ELSEIF(MVAR.EQ.4.OR.MVAR.EQ.6) THEN
0047           AUPP=ATAN((TAUMAX-TAURE)/GAMRE)
0048           ALOW=ATAN((TAUMIN-TAURE)/GAMRE)
0049           TAU=TAURE+GAMRE*TAN(ALOW+(AUPP-ALOW)*VVAR)
0050         ELSEIF(MINT(47).EQ.5) THEN
0051           AUPP=LOG(MAX(2D-10,1D0-TAUMAX))
0052           ALOW=LOG(MAX(2D-10,1D0-TAUMIN))
0053           TAU=1D0-EXP(AUPP+VVAR*(ALOW-AUPP))
0054         ELSE
0055           AUPP=LOG(MAX(1D-10,1D0-TAUMAX))
0056           ALOW=LOG(MAX(1D-10,1D0-TAUMIN))
0057           TAU=1D0-EXP(AUPP+VVAR*(ALOW-AUPP))
0058         ENDIF
0059         VINT(21)=MIN(TAUMAX,MAX(TAUMIN,TAU))
0060
0061 C...Convert VVAR to y* variable.
0062       ELSEIF(IVAR.EQ.2) THEN
0063         YSTMIN=VINT(12)
0064         YSTMAX=VINT(32)
0065         TAUE=VINT(21)
0066         IF(ISTSB.GE.3.AND.ISTSB.LE.5) TAUE=VINT(26)
0067         IF(MINT(47).EQ.1) THEN
0068           YST=0D0
0069         ELSEIF(MINT(47).EQ.2.OR.MINT(47).EQ.6) THEN
0070           YST=-0.5D0*LOG(TAUE)
0071         ELSEIF(MINT(47).EQ.3.OR.MINT(47).EQ.7) THEN
0072           YST=0.5D0*LOG(TAUE)
0073         ELSEIF(MVAR.EQ.1) THEN
0074           YST=YSTMIN+(YSTMAX-YSTMIN)*SQRT(VVAR)
0075         ELSEIF(MVAR.EQ.2) THEN
0076           YST=YSTMAX-(YSTMAX-YSTMIN)*SQRT(1D0-VVAR)
0077         ELSEIF(MVAR.EQ.3) THEN
0078           AUPP=ATAN(EXP(YSTMAX))
0079           ALOW=ATAN(EXP(YSTMIN))
0080           YST=LOG(TAN(ALOW+(AUPP-ALOW)*VVAR))
0081         ELSEIF(MVAR.EQ.4) THEN
0082           YST0=-0.5D0*LOG(TAUE)
0083           AUPP=LOG(MAX(1D-10,EXP(YST0-YSTMIN)-1D0))
0084           ALOW=LOG(MAX(1D-10,EXP(YST0-YSTMAX)-1D0))
0085           YST=YST0-LOG(1D0+EXP(ALOW+VVAR*(AUPP-ALOW)))
0086         ELSE
0087           YST0=-0.5D0*LOG(TAUE)
0088           AUPP=LOG(MAX(1D-10,EXP(YST0+YSTMIN)-1D0))
0089           ALOW=LOG(MAX(1D-10,EXP(YST0+YSTMAX)-1D0))
0090           YST=LOG(1D0+EXP(AUPP+VVAR*(ALOW-AUPP)))-YST0
0091         ENDIF
0092         VINT(22)=MIN(YSTMAX,MAX(YSTMIN,YST))
0093
0094 C...Convert VVAR to cos(theta-hat) variable.
0095       ELSEIF(IVAR.EQ.3) THEN
0096         RM34=MAX(1D-20,2D0*VINT(63)*VINT(64)/(VINT(21)*VINT(2))**2)
0097         RSQM=1D0+RM34
0098         IF(2D0*VINT(71)**2/(VINT(21)*VINT(2)).LT.0.0001D0)
0099      &  RM34=MAX(RM34,2D0*VINT(71)**2/(VINT(21)*VINT(2)))
0100         CTNMIN=VINT(13)
0101         CTNMAX=VINT(33)
0102         CTPMIN=VINT(14)
0103         CTPMAX=VINT(34)
0104         IF(MVAR.EQ.1) THEN
0105           ANEG=CTNMAX-CTNMIN
0106           APOS=CTPMAX-CTPMIN
0107           IF(ANEG.GT.0D0.AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
0108             VCTN=VVAR*(ANEG+APOS)/ANEG
0109             CTH=CTNMIN+(CTNMAX-CTNMIN)*VCTN
0110           ELSE
0111             VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
0112             CTH=CTPMIN+(CTPMAX-CTPMIN)*VCTP
0113           ENDIF
0114         ELSEIF(MVAR.EQ.2) THEN
0115           RMNMIN=MAX(RM34,RSQM-CTNMIN)
0116           RMNMAX=MAX(RM34,RSQM-CTNMAX)
0117           RMPMIN=MAX(RM34,RSQM-CTPMIN)
0118           RMPMAX=MAX(RM34,RSQM-CTPMAX)
0119           ANEG=LOG(RMNMIN/RMNMAX)
0120           APOS=LOG(RMPMIN/RMPMAX)
0121           IF(ANEG.GT.0D0.AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
0122             VCTN=VVAR*(ANEG+APOS)/ANEG
0123             CTH=RSQM-RMNMIN*(RMNMAX/RMNMIN)**VCTN
0124           ELSE
0125             VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
0126             CTH=RSQM-RMPMIN*(RMPMAX/RMPMIN)**VCTP
0127           ENDIF
0128         ELSEIF(MVAR.EQ.3) THEN
0129           RMNMIN=MAX(RM34,RSQM+CTNMIN)
0130           RMNMAX=MAX(RM34,RSQM+CTNMAX)
0131           RMPMIN=MAX(RM34,RSQM+CTPMIN)
0132           RMPMAX=MAX(RM34,RSQM+CTPMAX)
0133           ANEG=LOG(RMNMAX/RMNMIN)
0134           APOS=LOG(RMPMAX/RMPMIN)
0135           IF(ANEG.GT.0D0.AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
0136             VCTN=VVAR*(ANEG+APOS)/ANEG
0137             CTH=RMNMIN*(RMNMAX/RMNMIN)**VCTN-RSQM
0138           ELSE
0139             VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
0140             CTH=RMPMIN*(RMPMAX/RMPMIN)**VCTP-RSQM
0141           ENDIF
0142         ELSEIF(MVAR.EQ.4) THEN
0143           RMNMIN=MAX(RM34,RSQM-CTNMIN)
0144           RMNMAX=MAX(RM34,RSQM-CTNMAX)
0145           RMPMIN=MAX(RM34,RSQM-CTPMIN)
0146           RMPMAX=MAX(RM34,RSQM-CTPMAX)
0147           ANEG=1D0/RMNMAX-1D0/RMNMIN
0148           APOS=1D0/RMPMAX-1D0/RMPMIN
0149           IF(ANEG.GT.0D0.AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
0150             VCTN=VVAR*(ANEG+APOS)/ANEG
0151             CTH=RSQM-1D0/(1D0/RMNMIN+ANEG*VCTN)
0152           ELSE
0153             VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
0154             CTH=RSQM-1D0/(1D0/RMPMIN+APOS*VCTP)
0155           ENDIF
0156         ELSEIF(MVAR.EQ.5) THEN
0157           RMNMIN=MAX(RM34,RSQM+CTNMIN)
0158           RMNMAX=MAX(RM34,RSQM+CTNMAX)
0159           RMPMIN=MAX(RM34,RSQM+CTPMIN)
0160           RMPMAX=MAX(RM34,RSQM+CTPMAX)
0161           ANEG=1D0/RMNMIN-1D0/RMNMAX
0162           APOS=1D0/RMPMIN-1D0/RMPMAX
0163           IF(ANEG.GT.0D0.AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
0164             VCTN=VVAR*(ANEG+APOS)/ANEG
0165             CTH=1D0/(1D0/RMNMIN-ANEG*VCTN)-RSQM
0166           ELSE
0167             VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
0168             CTH=1D0/(1D0/RMPMIN-APOS*VCTP)-RSQM
0169           ENDIF
0170         ENDIF
0171         IF(CTH.LT.0D0) CTH=MIN(CTNMAX,MAX(CTNMIN,CTH))
0172         IF(CTH.GT.0D0) CTH=MIN(CTPMAX,MAX(CTPMIN,CTH))
0173         VINT(23)=CTH
0174
0175 C...Convert VVAR to tau' variable.
0176       ELSEIF(IVAR.EQ.4) THEN
0177         TAU=VINT(21)
0178         TAUPMN=VINT(16)
0179         TAUPMX=VINT(36)
0180         IF(MINT(47).EQ.1) THEN
0181           TAUP=1D0
0182         ELSEIF(MVAR.EQ.1) THEN
0183           TAUP=TAUPMN*(TAUPMX/TAUPMN)**VVAR
0184         ELSEIF(MVAR.EQ.2) THEN
0185           AUPP=(1D0-TAU/TAUPMX)**4
0186           ALOW=(1D0-TAU/TAUPMN)**4
0187           TAUP=TAU/MAX(1D-10,1D0-(ALOW+(AUPP-ALOW)*VVAR)**0.25D0)
0188         ELSEIF(MINT(47).EQ.5) THEN
0189           AUPP=LOG(MAX(2D-10,1D0-TAUPMX))
0190           ALOW=LOG(MAX(2D-10,1D0-TAUPMN))
0191           TAUP=1D0-EXP(AUPP+VVAR*(ALOW-AUPP))
0192         ELSE
0193           AUPP=LOG(MAX(1D-10,1D0-TAUPMX))
0194           ALOW=LOG(MAX(1D-10,1D0-TAUPMN))
0195           TAUP=1D0-EXP(AUPP+VVAR*(ALOW-AUPP))
0196         ENDIF
0197         VINT(26)=MIN(TAUPMX,MAX(TAUPMN,TAUP))
0198
0199 C...Selection of extra variables needed in 2 -> 3 process:
0200 C...pT1, pT2, phi1, phi2, y3 for three outgoing particles.
0201 C...Since no options are available, the functions of PYKLIM
0202 C...and PYKMAP are joint for these choices.
0203       ELSEIF(IVAR.EQ.5) THEN
0204
0205 C...Read out total energy and particle masses.
0206         MINT(51)=0
0207         MPTPK=1
0208         IF(ISUB.EQ.123.OR.ISUB.EQ.124.OR.ISUB.EQ.173.OR.ISUB.EQ.174
0209      &  .OR.ISUB.EQ.178.OR.ISUB.EQ.179.OR.ISUB.EQ.351.OR.ISUB.EQ.352)
0210      &  MPTPK=2
0211         SHP=VINT(26)*VINT(2)
0212         SHPR=SQRT(SHP)
0213         PM1=VINT(201)
0214         PM2=VINT(206)
0215         PM3=SQRT(VINT(21))*VINT(1)
0216         IF(PM1+PM2+PM3.GT.0.9999D0*SHPR) THEN
0217           MINT(51)=1
0218           RETURN
0219         ENDIF
0220         PMRS1=VINT(204)**2
0221         PMRS2=VINT(209)**2
0222
0223 C...Specify coefficients of pT choice; upper and lower limits.
0224         IF(MPTPK.EQ.1) THEN
0225           HWT1=0.4D0
0226           HWT2=0.4D0
0227         ELSE
0228           HWT1=0.05D0
0229           HWT2=0.05D0
0230         ENDIF
0231         HWT3=1D0-HWT1-HWT2
0232         PTSMX1=((SHP-PM1**2-(PM2+PM3)**2)**2-(2D0*PM1*(PM2+PM3))**2)/
0233      &  (4D0*SHP)
0234         IF(CKIN(52).GT.0D0) PTSMX1=MIN(PTSMX1,CKIN(52)**2)
0235         PTSMN1=CKIN(51)**2
0236         PTSMX2=((SHP-PM2**2-(PM1+PM3)**2)**2-(2D0*PM2*(PM1+PM3))**2)/
0237      &  (4D0*SHP)
0238         IF(CKIN(54).GT.0D0) PTSMX2=MIN(PTSMX2,CKIN(54)**2)
0239         PTSMN2=CKIN(53)**2
0240
0241 C...Select transverse momenta according to
0242 C...dp_T^2 * (a + b/(M^2 + p_T^2) + c/(M^2 + p_T^2)^2).
0243         HMX=PMRS1+PTSMX1
0244         HMN=PMRS1+PTSMN1
0245         IF(HMX.LT.1.0001D0*HMN) THEN
0246           MINT(51)=1
0247           RETURN
0248         ENDIF
0249         HDE=PTSMX1-PTSMN1
0250         RPT=PYR(0)
0251         IF(RPT.LT.HWT1) THEN
0252           PTS1=PTSMN1+PYR(0)*HDE
0253         ELSEIF(RPT.LT.HWT1+HWT2) THEN
0254           PTS1=MAX(PTSMN1,HMN*(HMX/HMN)**PYR(0)-PMRS1)
0255         ELSE
0256           PTS1=MAX(PTSMN1,HMN*HMX/(HMN+PYR(0)*HDE)-PMRS1)
0257         ENDIF
0258         WTPTS1=HDE/(HWT1+HWT2*HDE/(LOG(HMX/HMN)*(PMRS1+PTS1))+
0259      &  HWT3*HMN*HMX/(PMRS1+PTS1)**2)
0260         HMX=PMRS2+PTSMX2
0261         HMN=PMRS2+PTSMN2
0262         IF(HMX.LT.1.0001D0*HMN) THEN
0263           MINT(51)=1
0264           RETURN
0265         ENDIF
0266         HDE=PTSMX2-PTSMN2
0267         RPT=PYR(0)
0268         IF(RPT.LT.HWT1) THEN
0269           PTS2=PTSMN2+PYR(0)*HDE
0270         ELSEIF(RPT.LT.HWT1+HWT2) THEN
0271           PTS2=MAX(PTSMN2,HMN*(HMX/HMN)**PYR(0)-PMRS2)
0272         ELSE
0273           PTS2=MAX(PTSMN2,HMN*HMX/(HMN+PYR(0)*HDE)-PMRS2)
0274         ENDIF
0275         WTPTS2=HDE/(HWT1+HWT2*HDE/(LOG(HMX/HMN)*(PMRS2+PTS2))+
0276      &  HWT3*HMN*HMX/(PMRS2+PTS2)**2)
0277
0278 C...Select azimuthal angles and check pT choice.
0279         PHI1=PARU(2)*PYR(0)
0280         PHI2=PARU(2)*PYR(0)
0281         PHIR=PHI2-PHI1
0282         PTS3=MAX(0D0,PTS1+PTS2+2D0*SQRT(PTS1*PTS2)*COS(PHIR))
0283         IF(PTS3.LT.CKIN(55)**2.OR.(CKIN(56).GT.0D0.AND.PTS3.GT.
0284      &  CKIN(56)**2)) THEN
0285           MINT(51)=1
0286           RETURN
0287         ENDIF
0288
0289 C...Calculate transverse masses and check phase space not closed.
0290         PMS1=PM1**2+PTS1
0291         PMS2=PM2**2+PTS2
0292         PMS3=PM3**2+PTS3
0293         PMT1=SQRT(PMS1)
0294         PMT2=SQRT(PMS2)
0295         PMT3=SQRT(PMS3)
0296         PM12=(PMT1+PMT2)**2
0297         IF(PMT1+PMT2+PMT3.GT.0.9999D0*SHPR) THEN
0298           MINT(51)=1
0299           RETURN
0300         ENDIF
0301
0302 C...Select rapidity for particle 3 and check phase space not closed.
0303         Y3MAX=LOG((SHP+PMS3-PM12+SQRT(MAX(0D0,(SHP-PMS3-PM12)**2-
0304      &  4D0*PMS3*PM12)))/(2D0*SHPR*PMT3))
0305         IF(Y3MAX.LT.1D-6) THEN
0306           MINT(51)=1
0307           RETURN
0308         ENDIF
0309         Y3=(2D0*PYR(0)-1D0)*0.999999D0*Y3MAX
0310         PZ3=PMT3*SINH(Y3)
0311         PE3=PMT3*COSH(Y3)
0312
0313 C...Find momentum transfers in two mirror solutions (in 1-2 frame).
0314         PZ12=-PZ3
0315         PE12=SHPR-PE3
0316         PMS12=PE12**2-PZ12**2
0317         SQL12=SQRT(MAX(0D0,(PMS12-PMS1-PMS2)**2-4D0*PMS1*PMS2))
0318         IF(SQL12.LT.1D-6*SHP) THEN
0319           MINT(51)=1
0320           RETURN
0321         ENDIF
0322         PMM1=PMS12+PMS1-PMS2
0323         PMM2=PMS12+PMS2-PMS1
0324         TFAC=-SHPR/(2D0*PMS12)
0325         T1P=TFAC*(PE12-PZ12)*(PMM1-SQL12)
0326         T1N=TFAC*(PE12-PZ12)*(PMM1+SQL12)
0327         T2P=TFAC*(PE12+PZ12)*(PMM2-SQL12)
0328         T2N=TFAC*(PE12+PZ12)*(PMM2+SQL12)
0329
0330 C...Construct relative mirror weights and make choice.
0331         IF(MPTPK.EQ.1.OR.ISUB.EQ.351.OR.ISUB.EQ.352) THEN
0332           WTPU=1D0
0333           WTNU=1D0
0334         ELSE
0335           WTPU=1D0/((T1P-PMRS1)*(T2P-PMRS2))**2
0336           WTNU=1D0/((T1N-PMRS1)*(T2N-PMRS2))**2
0337         ENDIF
0338         WTP=WTPU/(WTPU+WTNU)
0339         WTN=WTNU/(WTPU+WTNU)
0340         EPS=1D0
0341         IF(WTN.GT.PYR(0)) EPS=-1D0
0342
0343 C...Store result of variable choice and associated weights.
0344         VINT(202)=PTS1
0345         VINT(207)=PTS2
0346         VINT(203)=PHI1
0347         VINT(208)=PHI2
0348         VINT(205)=WTPTS1
0349         VINT(210)=WTPTS2
0350         VINT(211)=Y3
0351         VINT(212)=Y3MAX
0352         VINT(213)=EPS
0353         IF(EPS.GT.0D0) THEN
0354           VINT(214)=1D0/WTP
0355           VINT(215)=T1P
0356           VINT(216)=T2P
0357         ELSE
0358           VINT(214)=1D0/WTN
0359           VINT(215)=T1N
0360           VINT(216)=T2N
0361         ENDIF
0362         VINT(217)=-0.5D0*TFAC*(PE12-PZ12)*(PMM2+EPS*SQL12)
0363         VINT(218)=-0.5D0*TFAC*(PE12+PZ12)*(PMM1+EPS*SQL12)
0364         VINT(219)=0.5D0*(PMS12-PTS3)
0365         VINT(220)=SQL12
0366       ENDIF
0367
0368       RETURN
0369       END