PYTHIA/pythia/pykmap.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE PYKMAP(IVAR,MVAR,VVAR)
   5
   6 C...Maps a uniform distribution into a distribution of a kinematical
   7 C...variable according to one of the possibilities allowed. It is
   8 C...assumed that kinematical limits have been set by a PYKLIM call.
   9       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  10       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  11       COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
  12       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
  13       COMMON/PYINT1/MINT(400),VINT(400)
  14       COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
  15       SAVE /LUDAT1/,/LUDAT2/
  16       SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/
  17
  18 C...Convert VVAR to tau variable.
  19       ISUB=MINT(1)
  20       ISTSB=ISET(ISUB)
  21       IF(IVAR.EQ.1) THEN
  22         TAUMIN=VINT(11)
  23         TAUMAX=VINT(31)
  24         IF(MVAR.EQ.3.OR.MVAR.EQ.4) THEN
  25           TAURE=VINT(73)
  26           GAMRE=VINT(74)
  27         ELSEIF(MVAR.EQ.5.OR.MVAR.EQ.6) THEN
  28           TAURE=VINT(75)
  29           GAMRE=VINT(76)
  30         ENDIF
  31         IF(MINT(47).EQ.1.AND.(ISTSB.EQ.1.OR.ISTSB.EQ.2.OR.ISTSB.EQ.6))
  32      &  THEN
  33           TAU=1.
  34         ELSEIF(MVAR.EQ.1) THEN
  35           TAU=TAUMIN*(TAUMAX/TAUMIN)**VVAR
  36         ELSEIF(MVAR.EQ.2) THEN
  37           TAU=TAUMAX*TAUMIN/(TAUMIN+(TAUMAX-TAUMIN)*VVAR)
  38         ELSEIF(MVAR.EQ.3.OR.MVAR.EQ.5) THEN
  39           RATGEN=(TAURE+TAUMAX)/(TAURE+TAUMIN)*TAUMIN/TAUMAX
  40           TAU=TAURE*TAUMIN/((TAURE+TAUMIN)*RATGEN**VVAR-TAUMIN)
  41         ELSEIF(MVAR.EQ.4.OR.MVAR.EQ.6) THEN
  42           AUPP=ATAN((TAUMAX-TAURE)/GAMRE)
  43           ALOW=ATAN((TAUMIN-TAURE)/GAMRE)
  44           TAU=TAURE+GAMRE*TAN(ALOW+(AUPP-ALOW)*VVAR)
  45         ELSE
  46           AUPP=LOG(MAX(2E-6,1.-TAUMAX))
  47           ALOW=LOG(MAX(2E-6,1.-TAUMIN))
  48           TAU=1.-EXP(AUPP+VVAR*(ALOW-AUPP))
  49         ENDIF
  50         VINT(21)=MIN(TAUMAX,MAX(TAUMIN,TAU))
  51
  52 C...Convert VVAR to y* variable.
  53       ELSEIF(IVAR.EQ.2) THEN
  54         YSTMIN=VINT(12)
  55         YSTMAX=VINT(32)
  56         TAUE=VINT(21)
  57         IF(ISTSB.GE.3.AND.ISTSB.LE.5) TAUE=VINT(26)
  58         IF(MINT(47).EQ.1) THEN
  59           YST=0.
  60         ELSEIF(MINT(47).EQ.2) THEN
  61           YST=-0.5*LOG(TAUE)
  62         ELSEIF(MINT(47).EQ.3) THEN
  63           YST=0.5*LOG(TAUE)
  64         ELSEIF(MVAR.EQ.1) THEN
  65           YST=YSTMIN+(YSTMAX-YSTMIN)*SQRT(VVAR)
  66         ELSEIF(MVAR.EQ.2) THEN
  67           YST=YSTMAX-(YSTMAX-YSTMIN)*SQRT(1.-VVAR)
  68         ELSEIF(MVAR.EQ.3) THEN
  69           AUPP=ATAN(EXP(YSTMAX))
  70           ALOW=ATAN(EXP(YSTMIN))
  71           YST=LOG(TAN(ALOW+(AUPP-ALOW)*VVAR))
  72         ELSEIF(MVAR.EQ.4) THEN
  73           YST0=-0.5*LOG(TAUE)
  74           AUPP=LOG(MAX(1E-6,EXP(YST0-YSTMIN)-1.))
  75           ALOW=LOG(MAX(1E-6,EXP(YST0-YSTMAX)-1.))
  76           YST=YST0-LOG(1.+EXP(ALOW+VVAR*(AUPP-ALOW)))
  77         ELSE
  78           YST0=-0.5*LOG(TAUE)
  79           AUPP=LOG(MAX(1E-6,EXP(YST0+YSTMIN)-1.))
  80           ALOW=LOG(MAX(1E-6,EXP(YST0+YSTMAX)-1.))
  81           YST=LOG(1.+EXP(AUPP+VVAR*(ALOW-AUPP)))-YST0
  82         ENDIF
  83         VINT(22)=MIN(YSTMAX,MAX(YSTMIN,YST))
  84
  85 C...Convert VVAR to cos(theta-hat) variable.
  86       ELSEIF(IVAR.EQ.3) THEN
  87         RM34=MAX(1E-20,2.*VINT(63)*VINT(64)/(VINT(21)*VINT(2))**2)
  88         RSQM=1.+RM34
  89         IF(2.*VINT(71)**2/(VINT(21)*VINT(2)).LT.0.0001) RM34=MAX(RM34,
  90      &  2.*VINT(71)**2/(VINT(21)*VINT(2)))
  91         CTNMIN=VINT(13)
  92         CTNMAX=VINT(33)
  93         CTPMIN=VINT(14)
  94         CTPMAX=VINT(34)
  95         IF(MVAR.EQ.1) THEN
  96           ANEG=CTNMAX-CTNMIN
  97           APOS=CTPMAX-CTPMIN
  98           IF(ANEG.GT.0..AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
  99             VCTN=VVAR*(ANEG+APOS)/ANEG
 100             CTH=CTNMIN+(CTNMAX-CTNMIN)*VCTN
 101           ELSE
 102             VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
 103             CTH=CTPMIN+(CTPMAX-CTPMIN)*VCTP
 104           ENDIF
 105         ELSEIF(MVAR.EQ.2) THEN
 106           RMNMIN=MAX(RM34,RSQM-CTNMIN)
 107           RMNMAX=MAX(RM34,RSQM-CTNMAX)
 108           RMPMIN=MAX(RM34,RSQM-CTPMIN)
 109           RMPMAX=MAX(RM34,RSQM-CTPMAX)
 110           ANEG=LOG(RMNMIN/RMNMAX)
 111           APOS=LOG(RMPMIN/RMPMAX)
 112           IF(ANEG.GT.0..AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
 113             VCTN=VVAR*(ANEG+APOS)/ANEG
 114             CTH=RSQM-RMNMIN*(RMNMAX/RMNMIN)**VCTN
 115           ELSE
 116             VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
 117             CTH=RSQM-RMPMIN*(RMPMAX/RMPMIN)**VCTP
 118           ENDIF
 119         ELSEIF(MVAR.EQ.3) THEN
 120           RMNMIN=MAX(RM34,RSQM+CTNMIN)
 121           RMNMAX=MAX(RM34,RSQM+CTNMAX)
 122           RMPMIN=MAX(RM34,RSQM+CTPMIN)
 123           RMPMAX=MAX(RM34,RSQM+CTPMAX)
 124           ANEG=LOG(RMNMAX/RMNMIN)
 125           APOS=LOG(RMPMAX/RMPMIN)
 126           IF(ANEG.GT.0..AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
 127             VCTN=VVAR*(ANEG+APOS)/ANEG
 128             CTH=RMNMIN*(RMNMAX/RMNMIN)**VCTN-RSQM
 129           ELSE
 130             VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
 131             CTH=RMPMIN*(RMPMAX/RMPMIN)**VCTP-RSQM
 132           ENDIF
 133         ELSEIF(MVAR.EQ.4) THEN
 134           RMNMIN=MAX(RM34,RSQM-CTNMIN)
 135           RMNMAX=MAX(RM34,RSQM-CTNMAX)
 136           RMPMIN=MAX(RM34,RSQM-CTPMIN)
 137           RMPMAX=MAX(RM34,RSQM-CTPMAX)
 138           ANEG=1./RMNMAX-1./RMNMIN
 139           APOS=1./RMPMAX-1./RMPMIN
 140           IF(ANEG.GT.0..AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
 141             VCTN=VVAR*(ANEG+APOS)/ANEG
 142             CTH=RSQM-1./(1./RMNMIN+ANEG*VCTN)
 143           ELSE
 144             VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
 145             CTH=RSQM-1./(1./RMPMIN+APOS*VCTP)
 146           ENDIF
 147         ELSEIF(MVAR.EQ.5) THEN
 148           RMNMIN=MAX(RM34,RSQM+CTNMIN)
 149           RMNMAX=MAX(RM34,RSQM+CTNMAX)
 150           RMPMIN=MAX(RM34,RSQM+CTPMIN)
 151           RMPMAX=MAX(RM34,RSQM+CTPMAX)
 152           ANEG=1./RMNMIN-1./RMNMAX
 153           APOS=1./RMPMIN-1./RMPMAX
 154           IF(ANEG.GT.0..AND.VVAR*(ANEG+APOS).LE.ANEG) THEN
 155             VCTN=VVAR*(ANEG+APOS)/ANEG
 156             CTH=1./(1./RMNMIN-ANEG*VCTN)-RSQM
 157           ELSE
 158             VCTP=(VVAR*(ANEG+APOS)-ANEG)/APOS
 159             CTH=1./(1./RMPMIN-APOS*VCTP)-RSQM
 160           ENDIF
 161         ENDIF
 162         IF(CTH.LT.0.) CTH=MIN(CTNMAX,MAX(CTNMIN,CTH))
 163         IF(CTH.GT.0.) CTH=MIN(CTPMAX,MAX(CTPMIN,CTH))
 164         VINT(23)=CTH
 165
 166 C...Convert VVAR to tau' variable.
 167       ELSEIF(IVAR.EQ.4) THEN
 168         TAU=VINT(21)
 169         TAUPMN=VINT(16)
 170         TAUPMX=VINT(36)
 171         IF(MINT(47).EQ.1) THEN
 172           TAUP=1.
 173         ELSEIF(MVAR.EQ.1) THEN
 174           TAUP=TAUPMN*(TAUPMX/TAUPMN)**VVAR
 175         ELSEIF(MVAR.EQ.2) THEN
 176           AUPP=(1.-TAU/TAUPMX)**4
 177           ALOW=(1.-TAU/TAUPMN)**4
 178           TAUP=TAU/MAX(1E-7,1.-(ALOW+(AUPP-ALOW)*VVAR)**0.25)
 179         ELSE
 180           AUPP=LOG(MAX(2E-6,1.-TAUPMX))
 181           ALOW=LOG(MAX(2E-6,1.-TAUPMN))
 182           TAUP=1.-EXP(AUPP+VVAR*(ALOW-AUPP))
 183         ENDIF
 184         VINT(26)=MIN(TAUPMX,MAX(TAUPMN,TAUP))
 185
 186 C...Selection of extra variables needed in 2 -> 3 process:
 187 C...pT1, pT2, phi1, phi2, y3 for three outgoing particles.
 188 C...Since no options are available, the functions of PYKLIM
 189 C...and PYKMAP are joint for these choices.
 190       ELSEIF(IVAR.EQ.5) THEN
 191
 192 C...Read out total energy and particle masses.
 193         MINT(51)=0
 194         MPTPK=1
 195         IF(ISUB.EQ.123.OR.ISUB.EQ.124.OR.ISUB.EQ.173.OR.ISUB.EQ.174
 196      &  .OR.ISUB.EQ.178.OR.ISUB.EQ.179) MPTPK=2
 197         SHP=VINT(26)*VINT(2)
 198         SHPR=SQRT(SHP)
 199         PM1=VINT(201)
 200         PM2=VINT(206)
 201         PM3=SQRT(VINT(21))*VINT(1)
 202         IF(PM1+PM2+PM3.GT.0.9999*SHPR) THEN
 203           MINT(51)=1
 204           RETURN
 205         ENDIF
 206         PMRS1=VINT(204)**2
 207         PMRS2=VINT(209)**2
 208
 209 C...Specify coefficients of pT choice; upper and lower limits.
 210         IF(MPTPK.EQ.1) THEN
 211           HWT1=0.4
 212           HWT2=0.4
 213         ELSE
 214           HWT1=0.05
 215           HWT2=0.05
 216         ENDIF
 217         HWT3=1.-HWT1-HWT2
 218         PTSMX1=((SHP-PM1**2-(PM2+PM3)**2)**2-(2.*PM1*(PM2+PM3))**2)/
 219      &  (4.*SHP)
 220         IF(CKIN(52).GT.0.) PTSMX1=MIN(PTSMX1,CKIN(52)**2)
 221         PTSMN1=CKIN(51)**2
 222         PTSMX2=((SHP-PM2**2-(PM1+PM3)**2)**2-(2.*PM2*(PM1+PM3))**2)/
 223      &  (4.*SHP)
 224         IF(CKIN(54).GT.0.) PTSMX2=MIN(PTSMX2,CKIN(54)**2)
 225         PTSMN2=CKIN(53)**2
 226
 227 C...Select transverse momenta according to
 228 C...dp_T^2 * (a + b/(M^2 + p_T^2) + c/(M^2 + p_T^2)^2).
 229         HMX=PMRS1+PTSMX1
 230         HMN=PMRS1+PTSMN1
 231         IF(HMX.LT.1.0001*HMN) THEN
 232           MINT(51)=1
 233           RETURN
 234         ENDIF
 235         HDE=PTSMX1-PTSMN1
 236         RPT=RLU(0)
 237         IF(RPT.LT.HWT1) THEN
 238           PTS1=PTSMN1+RLU(0)*HDE
 239         ELSEIF(RPT.LT.HWT1+HWT2) THEN
 240           PTS1=MAX(PTSMN1,HMN*(HMX/HMN)**RLU(0)-PMRS1)
 241         ELSE
 242           PTS1=MAX(PTSMN1,HMN*HMX/(HMN+RLU(0)*HDE)-PMRS1)
 243         ENDIF
 244         WTPTS1=HDE/(HWT1+HWT2*HDE/(LOG(HMX/HMN)*(PMRS1+PTS1))+
 245      &  HWT3*HMN*HMX/(PMRS1+PTS1)**2)
 246         HMX=PMRS2+PTSMX2
 247         HMN=PMRS2+PTSMN2
 248         IF(HMX.LT.1.0001*HMN) THEN
 249           MINT(51)=1
 250           RETURN
 251         ENDIF
 252         HDE=PTSMX2-PTSMN2
 253         RPT=RLU(0)
 254         IF(RPT.LT.HWT1) THEN
 255           PTS2=PTSMN2+RLU(0)*HDE
 256         ELSEIF(RPT.LT.HWT1+HWT2) THEN
 257           PTS2=MAX(PTSMN2,HMN*(HMX/HMN)**RLU(0)-PMRS2)
 258         ELSE
 259           PTS2=MAX(PTSMN2,HMN*HMX/(HMN+RLU(0)*HDE)-PMRS2)
 260         ENDIF
 261         WTPTS2=HDE/(HWT1+HWT2*HDE/(LOG(HMX/HMN)*(PMRS2+PTS2))+
 262      &  HWT3*HMN*HMX/(PMRS2+PTS2)**2)
 263
 264 C...Select azimuthal angles and check pT choice.
 265         PHI1=PARU(2)*RLU(0)
 266         PHI2=PARU(2)*RLU(0)
 267         PHIR=PHI2-PHI1
 268         PTS3=MAX(0.,PTS1+PTS2+2.*SQRT(PTS1*PTS2)*COS(PHIR))
 269         IF(PTS3.LT.CKIN(55)**2.OR.(CKIN(56).GT.0..AND.PTS3.GT.
 270      &  CKIN(56)**2)) THEN
 271           MINT(51)=1
 272           RETURN
 273         ENDIF
 274
 275 C...Calculate transverse masses and check phase space not closed.
 276         PMS1=PM1**2+PTS1
 277         PMS2=PM2**2+PTS2
 278         PMS3=PM3**2+PTS3
 279         PMT1=SQRT(PMS1)
 280         PMT2=SQRT(PMS2)
 281         PMT3=SQRT(PMS3)
 282         PM12=(PMT1+PMT2)**2
 283         IF(PMT1+PMT2+PMT3.GT.0.9999*SHPR) THEN
 284           MINT(51)=1
 285           RETURN
 286         ENDIF
 287
 288 C...Select rapidity for particle 3 and check phase space not closed.
 289         Y3MAX=LOG((SHP+PMS3-PM12+SQRT(MAX(0.,(SHP-PMS3-PM12)**2-
 290      &  4.*PMS3*PM12)))/(2.*SHPR*PMT3))
 291         IF(Y3MAX.LT.1E-6) THEN
 292           MINT(51)=1
 293           RETURN
 294         ENDIF
 295         Y3=(2.*RLU(0)-1.)*0.999999*Y3MAX
 296         PZ3=PMT3*SINH(Y3)
 297         PE3=PMT3*COSH(Y3)
 298
 299 C...Find momentum transfers in two mirror solutions (in 1-2 frame).
 300         PZ12=-PZ3
 301         PE12=SHPR-PE3
 302         PMS12=PE12**2-PZ12**2
 303         SQL12=SQRT(MAX(0.,(PMS12-PMS1-PMS2)**2-4.*PMS1*PMS2))
 304         IF(SQL12.LT.1E-6*SHP) THEN
 305           MINT(51)=1
 306           RETURN
 307         ENDIF
 308         PMM1=PMS12+PMS1-PMS2
 309         PMM2=PMS12+PMS2-PMS1
 310         TFAC=-SHPR/(2.*PMS12)
 311         T1P=TFAC*(PE12-PZ12)*(PMM1-SQL12)
 312         T1N=TFAC*(PE12-PZ12)*(PMM1+SQL12)
 313         T2P=TFAC*(PE12+PZ12)*(PMM2-SQL12)
 314         T2N=TFAC*(PE12+PZ12)*(PMM2+SQL12)
 315
 316 C...Construct relative mirror weights and make choice.
 317         IF(MPTPK.EQ.1) THEN
 318           WTPU=1.
 319           WTNU=1.
 320         ELSE
 321           WTPU=1./((T1P-PMRS1)*(T2P-PMRS2))**2
 322           WTNU=1./((T1N-PMRS1)*(T2N-PMRS2))**2
 323         ENDIF
 324         WTP=WTPU/(WTPU+WTNU)
 325         WTN=WTNU/(WTPU+WTNU)
 326         EPS=1.
 327         IF(WTN.GT.RLU(0)) EPS=-1.
 328
 329 C...Store result of variable choice and associated weights.
 330         VINT(202)=PTS1
 331         VINT(207)=PTS2
 332         VINT(203)=PHI1
 333         VINT(208)=PHI2
 334         VINT(205)=WTPTS1
 335         VINT(210)=WTPTS2
 336         VINT(211)=Y3
 337         VINT(212)=Y3MAX
 338         VINT(213)=EPS
 339         IF(EPS.GT.0.) THEN
 340           VINT(214)=1./WTP
 341           VINT(215)=T1P
 342           VINT(216)=T2P
 343         ELSE
 344           VINT(214)=1./WTN
 345           VINT(215)=T1N
 346           VINT(216)=T2N
 347         ENDIF
 348         VINT(217)=-0.5*TFAC*(PE12-PZ12)*(PMM2+EPS*SQL12)
 349         VINT(218)=-0.5*TFAC*(PE12+PZ12)*(PMM1+EPS*SQL12)
 350         VINT(219)=0.5*(PMS12-PTS3)
 351         VINT(220)=SQL12
 352       ENDIF
 353
 354       RETURN
 355       END