PYTHIA/pythia/pymaxi.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE PYMAXI
   5
   6 C...Finds optimal set of coefficients for kinematical variable selection
   7 C...and the maximum of the part of the differential cross-section used
   8 C...in the event weighting.
   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       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
  16       COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
  17       COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
  18       COMMON/PYINT6/PROC(0:200)
  19       CHARACTER PROC*28
  20       COMMON/PYINT7/SIGT(0:6,0:6,0:5)
  21       SAVE /LUDAT1/,/LUDAT2/
  22       SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/,/PYINT4/,
  23      &/PYINT5/,/PYINT6/,/PYINT7/
  24       CHARACTER CVAR(4)*4
  25       DIMENSION NPTS(4),MVARPT(500,4),VINTPT(500,30),SIGSPT(500),
  26      &NAREL(7),WTREL(7),WTMAT(7,7),WTRELN(7),COEFU(7),COEFO(7),
  27      &IACCMX(4),SIGSMX(4),SIGSSM(3)
  28       DATA CVAR/'tau ','tau''','y*  ','cth '/
  29       DATA SIGSSM/3*0./
  30
  31 C...Select subprocess to study: skip cases not applicable.
  32       NPOSI=0
  33       VINT(143)=1.
  34       VINT(144)=1.
  35       XSEC(0,1)=0.
  36       DO 440 ISUB=1,200
  37       MINT(51)=0
  38       IF(ISET(ISUB).EQ.11) THEN
  39         XSEC(ISUB,1)=1.00001*COEF(ISUB,1)
  40         NPOSI=NPOSI+1
  41         GOTO 430
  42       ELSEIF(ISUB.GE.91.AND.ISUB.LE.95) THEN
  43         XSEC(ISUB,1)=SIGT(0,0,ISUB-90)
  44         IF(MSUB(ISUB).NE.1) GOTO 440
  45         NPOSI=NPOSI+1
  46         GOTO 430
  47       ELSEIF(ISUB.EQ.96) THEN
  48         IF(MINT(50).EQ.0) GOTO 440
  49         IF(MSUB(95).NE.1.AND.MSTP(81).LE.0.AND.MSTP(131).LE.0) GOTO 440
  50         IF(MINT(49).EQ.0.AND.MSTP(131).EQ.0) GOTO 440
  51       ELSEIF(ISUB.EQ.11.OR.ISUB.EQ.12.OR.ISUB.EQ.13.OR.ISUB.EQ.28.OR.
  52      &ISUB.EQ.53.OR.ISUB.EQ.68) THEN
  53         IF(MSUB(ISUB).NE.1.OR.MSUB(95).EQ.1) GOTO 440
  54       ELSE
  55         IF(MSUB(ISUB).NE.1) GOTO 440
  56       ENDIF
  57       MINT(1)=ISUB
  58       ISTSB=ISET(ISUB)
  59       IF(ISUB.EQ.96) ISTSB=2
  60       IF(MSTP(122).GE.2) WRITE(MSTU(11),5000) ISUB
  61       MWTXS=0
  62       IF(MSTP(142).GE.1.AND.ISUB.NE.96.AND.MSUB(91)+MSUB(92)+MSUB(93)+
  63      &MSUB(94)+MSUB(95).EQ.0) MWTXS=1
  64
  65 C...Find resonances (explicit or implicit in cross-section).
  66       MINT(72)=0
  67       KFR1=0
  68       IF(ISTSB.EQ.1.OR.ISTSB.EQ.3.OR.ISTSB.EQ.5) THEN
  69         KFR1=KFPR(ISUB,1)
  70       ELSEIF(ISUB.EQ.24.OR.ISUB.EQ.25.OR.ISUB.EQ.110.OR.ISUB.EQ.165.OR.
  71      &ISUB.EQ.171.OR.ISUB.EQ.176) THEN
  72         KFR1=23
  73       ELSEIF(ISUB.EQ.23.OR.ISUB.EQ.26.OR.ISUB.EQ.166.OR.ISUB.EQ.172.OR.
  74      &ISUB.EQ.177) THEN
  75         KFR1=24
  76       ELSEIF(ISUB.GE.71.AND.ISUB.LE.77) THEN
  77         KFR1=25
  78         IF(MSTP(46).EQ.5) THEN
  79           KFR1=30
  80           PMAS(30,1)=PARP(45)
  81           PMAS(30,2)=PARP(45)**3/(96.*PARU(1)*PARP(47)**2)
  82         ENDIF
  83       ENDIF
  84       CKMX=CKIN(2)
  85       IF(CKMX.LE.0.) CKMX=VINT(1)
  86       IF(KFR1.NE.0) THEN
  87         IF(CKIN(1).GT.PMAS(KFR1,1)+20.*PMAS(KFR1,2).OR.
  88      &  CKMX.LT.PMAS(KFR1,1)-20.*PMAS(KFR1,2)) KFR1=0
  89       ENDIF
  90       IF(KFR1.NE.0) THEN
  91         TAUR1=PMAS(KFR1,1)**2/VINT(2)
  92         GAMR1=PMAS(KFR1,1)*PMAS(KFR1,2)/VINT(2)
  93         MINT(72)=1
  94         MINT(73)=KFR1
  95         VINT(73)=TAUR1
  96         VINT(74)=GAMR1
  97       ENDIF
  98       IF(ISUB.EQ.141) THEN
  99         KFR2=23
 100         TAUR2=PMAS(KFR2,1)**2/VINT(2)
 101         GAMR2=PMAS(KFR2,1)*PMAS(KFR2,2)/VINT(2)
 102         IF(CKIN(1).GT.PMAS(KFR2,1)+20.*PMAS(KFR2,2).OR.
 103      &  CKMX.LT.PMAS(KFR2,1)-20.*PMAS(KFR2,2)) KFR2=0
 104         IF(KFR2.NE.0.AND.KFR1.NE.0) THEN
 105           MINT(72)=2
 106           MINT(74)=KFR2
 107           VINT(75)=TAUR2
 108           VINT(76)=GAMR2
 109         ELSEIF(KFR2.NE.0) THEN
 110           KFR1=KFR2
 111           TAUR1=TAUR2
 112           GAMR1=GAMR2
 113           MINT(72)=1
 114           MINT(73)=KFR1
 115           VINT(73)=TAUR1
 116           VINT(74)=GAMR1
 117         ENDIF
 118       ENDIF
 119
 120 C...Find product masses and minimum pT of process.
 121       SQM3=0.
 122       SQM4=0.
 123       MINT(71)=0
 124       VINT(71)=CKIN(3)
 125       VINT(80)=1.
 126       IF(ISTSB.EQ.2.OR.ISTSB.EQ.4) THEN
 127         NBW=0
 128         DO 100 I=1,2
 129         IF(KFPR(ISUB,I).EQ.0) THEN
 130         ELSEIF(MSTP(42).LE.0.OR.PMAS(LUCOMP(KFPR(ISUB,I)),2).LT.
 131      &  PARP(41)) THEN
 132           IF(I.EQ.1) SQM3=PMAS(LUCOMP(KFPR(ISUB,I)),1)**2
 133           IF(I.EQ.2) SQM4=PMAS(LUCOMP(KFPR(ISUB,I)),1)**2
 134         ELSE
 135           NBW=NBW+1
 136         ENDIF
 137   100   CONTINUE
 138         IF(NBW.GE.1) THEN
 139           CALL PYOFSH(3,0,KFPR(ISUB,1),KFPR(ISUB,2),0.,PQM3,PQM4)
 140           IF(MINT(51).EQ.1) THEN
 141             WRITE(MSTU(11),5100) ISUB
 142             MSUB(ISUB)=0
 143             GOTO 440
 144           ENDIF
 145           SQM3=PQM3**2
 146           SQM4=PQM4**2
 147         ENDIF
 148         IF(MIN(SQM3,SQM4).LT.CKIN(6)**2) MINT(71)=1
 149         IF(MINT(71).EQ.1) VINT(71)=MAX(CKIN(3),CKIN(5))
 150         IF(ISUB.EQ.96.AND.MSTP(82).LE.1) VINT(71)=PARP(81)
 151         IF(ISUB.EQ.96.AND.MSTP(82).GE.2) VINT(71)=0.08*PARP(82)
 152       ELSEIF(ISTSB.EQ.6) THEN
 153         CALL PYOFSH(5,0,KFPR(ISUB,1),KFPR(ISUB,2),0.,PQM3,PQM4)
 154         IF(MINT(51).EQ.1) THEN
 155           WRITE(MSTU(11),5100) ISUB
 156           MSUB(ISUB)=0
 157           GOTO 440
 158         ENDIF
 159         SQM3=PQM3**2
 160         SQM4=PQM4**2
 161       ENDIF
 162       VINT(63)=SQM3
 163       VINT(64)=SQM4
 164
 165 C...Prepare for additional variable choices in 2 -> 3.
 166       IF(ISTSB.EQ.5) THEN
 167         VINT(201)=0.
 168         IF(KFPR(ISUB,2).GT.0) VINT(201)=PMAS(KFPR(ISUB,2),1)
 169         VINT(206)=VINT(201)
 170         VINT(204)=PMAS(23,1)
 171         IF(ISUB.EQ.124) VINT(204)=PMAS(24,1)
 172         IF(ISUB.EQ.121.OR.ISUB.EQ.122.OR.ISUB.EQ.181.OR.ISUB.EQ.182.OR.
 173      &  ISUB.EQ.186.OR.ISUB.EQ.187) VINT(204)=VINT(201)
 174         VINT(209)=VINT(204)
 175       ENDIF
 176
 177 C...Number of points for each variable: tau, tau', y*, cos(theta-hat).
 178       NPTS(1)=2+2*MINT(72)
 179       IF(MINT(47).EQ.1) THEN
 180         IF(ISTSB.EQ.1.OR.ISTSB.EQ.2.OR.ISTSB.EQ.6) NPTS(1)=1
 181       ELSEIF(MINT(47).EQ.5) THEN
 182         IF(ISTSB.LE.2.OR.ISTSB.GE.6) NPTS(1)=NPTS(1)+1
 183       ENDIF
 184       NPTS(2)=1
 185       IF(ISTSB.GE.3.AND.ISTSB.LE.5) THEN
 186         IF(MINT(47).GE.2) NPTS(2)=2
 187         IF(MINT(47).EQ.5) NPTS(2)=3
 188       ENDIF
 189       NPTS(3)=1
 190       IF(MINT(47).GE.4) NPTS(3)=3
 191       IF(MINT(45).EQ.3) NPTS(3)=NPTS(3)+1
 192       IF(MINT(46).EQ.3) NPTS(3)=NPTS(3)+1
 193       NPTS(4)=1
 194       IF(ISTSB.EQ.2.OR.ISTSB.EQ.4.OR.ISTSB.EQ.6) NPTS(4)=5
 195       NTRY=NPTS(1)*NPTS(2)*NPTS(3)*NPTS(4)
 196
 197 C...Reset coefficients of cross-section weighting.
 198       DO 110 J=1,20
 199       COEF(ISUB,J)=0.
 200   110 CONTINUE
 201       COEF(ISUB,1)=1.
 202       COEF(ISUB,8)=0.5
 203       COEF(ISUB,9)=0.5
 204       COEF(ISUB,13)=1.
 205       COEF(ISUB,18)=1.
 206       MCTH=0
 207       MTAUP=0
 208       METAUP=0
 209       VINT(23)=0.
 210       VINT(26)=0.
 211       SIGSAM=0.
 212
 213 C...Find limits and select tau, y*, cos(theta-hat) and tau' values,
 214 C...in grid of phase space points.
 215       CALL PYKLIM(1)
 216       METAU=MINT(51)
 217       NACC=0
 218       DO 140 ITRY=1,NTRY
 219       MINT(51)=0
 220       IF(METAU.EQ.1) GOTO 140
 221       IF(MOD(ITRY-1,NPTS(2)*NPTS(3)*NPTS(4)).EQ.0) THEN
 222         MTAU=1+(ITRY-1)/(NPTS(2)*NPTS(3)*NPTS(4))
 223         IF(MTAU.GT.2+2*MINT(72)) MTAU=7
 224         CALL PYKMAP(1,MTAU,0.5)
 225         IF(ISTSB.GE.3.AND.ISTSB.LE.5) CALL PYKLIM(4)
 226         METAUP=MINT(51)
 227       ENDIF
 228       IF(METAUP.EQ.1) GOTO 140
 229       IF(ISTSB.GE.3.AND.ISTSB.LE.5.AND.MOD(ITRY-1,NPTS(3)*NPTS(4))
 230      &.EQ.0) THEN
 231         MTAUP=1+MOD((ITRY-1)/(NPTS(3)*NPTS(4)),NPTS(2))
 232         CALL PYKMAP(4,MTAUP,0.5)
 233       ENDIF
 234       IF(MOD(ITRY-1,NPTS(3)*NPTS(4)).EQ.0) THEN
 235         CALL PYKLIM(2)
 236         MEYST=MINT(51)
 237       ENDIF
 238       IF(MEYST.EQ.1) GOTO 140
 239       IF(MOD(ITRY-1,NPTS(4)).EQ.0) THEN
 240         MYST=1+MOD((ITRY-1)/NPTS(4),NPTS(3))
 241         IF(MYST.EQ.4.AND.MINT(45).NE.3) MYST=5
 242         CALL PYKMAP(2,MYST,0.5)
 243         CALL PYKLIM(3)
 244         MECTH=MINT(51)
 245       ENDIF
 246       IF(MECTH.EQ.1) GOTO 140
 247       IF(ISTSB.EQ.2.OR.ISTSB.EQ.4.OR.ISTSB.EQ.6) THEN
 248         MCTH=1+MOD(ITRY-1,NPTS(4))
 249         CALL PYKMAP(3,MCTH,0.5)
 250       ENDIF
 251       IF(ISUB.EQ.96) VINT(25)=VINT(21)*(1.-VINT(23)**2)
 252
 253 C...Store position and limits.
 254       MINT(51)=0
 255       CALL PYKLIM(0)
 256       IF(MINT(51).EQ.1) GOTO 140
 257       NACC=NACC+1
 258       MVARPT(NACC,1)=MTAU
 259       MVARPT(NACC,2)=MTAUP
 260       MVARPT(NACC,3)=MYST
 261       MVARPT(NACC,4)=MCTH
 262       DO 120 J=1,30
 263       VINTPT(NACC,J)=VINT(10+J)
 264   120 CONTINUE
 265
 266 C...Normal case: calculate cross-section.
 267       IF(ISTSB.NE.5) THEN
 268         CALL PYSIGH(NCHN,SIGS)
 269         IF(MWTXS.EQ.1) THEN
 270           CALL PYEVWT(WTXS)
 271           SIGS=WTXS*SIGS
 272         ENDIF
 273
 274 C..2 -> 3: find highest value out of a number of tries.
 275       ELSE
 276         SIGS=0.
 277         DO 130 IKIN3=1,MSTP(129)
 278         CALL PYKMAP(5,0,0.)
 279         IF(MINT(51).EQ.1) GOTO 130
 280         CALL PYSIGH(NCHN,SIGTMP)
 281         IF(MWTXS.EQ.1) THEN
 282           CALL PYEVWT(WTXS)
 283           SIGTMP=WTXS*SIGTMP
 284         ENDIF
 285         IF(SIGTMP.GT.SIGS) SIGS=SIGTMP
 286   130   CONTINUE
 287       ENDIF
 288
 289 C...Store cross-section.
 290       SIGSPT(NACC)=SIGS
 291       IF(SIGS.GT.SIGSAM) SIGSAM=SIGS
 292       IF(MSTP(122).GE.2) WRITE(MSTU(11),5200) MTAU,MYST,MCTH,MTAUP,
 293      &VINT(21),VINT(22),VINT(23),VINT(26),SIGS
 294   140 CONTINUE
 295       IF(NACC.EQ.0) THEN
 296         WRITE(MSTU(11),5100) ISUB
 297         MSUB(ISUB)=0
 298         GOTO 440
 299       ELSEIF(SIGSAM.EQ.0.) THEN
 300         WRITE(MSTU(11),5300) ISUB
 301         MSUB(ISUB)=0
 302         GOTO 440
 303       ENDIF
 304       IF(ISUB.NE.96) NPOSI=NPOSI+1
 305
 306 C...Calculate integrals in tau over maximal phase space limits.
 307       TAUMIN=VINT(11)
 308       TAUMAX=VINT(31)
 309       ATAU1=LOG(TAUMAX/TAUMIN)
 310       IF(NPTS(1).GE.2) THEN
 311         ATAU2=(TAUMAX-TAUMIN)/(TAUMAX*TAUMIN)
 312       ENDIF
 313       IF(NPTS(1).GE.4) THEN
 314         ATAU3=LOG(TAUMAX/TAUMIN*(TAUMIN+TAUR1)/(TAUMAX+TAUR1))/TAUR1
 315         ATAU4=(ATAN((TAUMAX-TAUR1)/GAMR1)-ATAN((TAUMIN-TAUR1)/GAMR1))/
 316      &  GAMR1
 317       ENDIF
 318       IF(NPTS(1).GE.6) THEN
 319         ATAU5=LOG(TAUMAX/TAUMIN*(TAUMIN+TAUR2)/(TAUMAX+TAUR2))/TAUR2
 320         ATAU6=(ATAN((TAUMAX-TAUR2)/GAMR2)-ATAN((TAUMIN-TAUR2)/GAMR2))/
 321      &  GAMR2
 322       ENDIF
 323       IF(NPTS(1).GT.2+2*MINT(72)) THEN
 324         ATAU7=LOG(MAX(2E-6,1.-TAUMIN)/MAX(2E-6,1.-TAUMAX))
 325       ENDIF
 326
 327 C...Reset. Sum up cross-sections in points calculated.
 328       DO 300 IVAR=1,4
 329       IF(NPTS(IVAR).EQ.1) GOTO 300
 330       IF(ISUB.EQ.96.AND.IVAR.EQ.4) GOTO 300
 331       NBIN=NPTS(IVAR)
 332       DO 160 J1=1,NBIN
 333       NAREL(J1)=0
 334       WTREL(J1)=0.
 335       COEFU(J1)=0.
 336       DO 150 J2=1,NBIN
 337       WTMAT(J1,J2)=0.
 338   150 CONTINUE
 339   160 CONTINUE
 340       DO 170 IACC=1,NACC
 341       IBIN=MVARPT(IACC,IVAR)
 342       IF(IVAR.EQ.1.AND.IBIN.EQ.7) IBIN=3+2*MINT(72)
 343       IF(IVAR.EQ.3.AND.IBIN.EQ.5.AND.MINT(45).NE.3) IBIN=4
 344       NAREL(IBIN)=NAREL(IBIN)+1
 345       WTREL(IBIN)=WTREL(IBIN)+SIGSPT(IACC)
 346
 347 C...Sum up tau cross-section pieces in points used.
 348       IF(IVAR.EQ.1) THEN
 349         TAU=VINTPT(IACC,11)
 350         WTMAT(IBIN,1)=WTMAT(IBIN,1)+1.
 351         WTMAT(IBIN,2)=WTMAT(IBIN,2)+(ATAU1/ATAU2)/TAU
 352         IF(NBIN.GE.4) THEN
 353           WTMAT(IBIN,3)=WTMAT(IBIN,3)+(ATAU1/ATAU3)/(TAU+TAUR1)
 354           WTMAT(IBIN,4)=WTMAT(IBIN,4)+(ATAU1/ATAU4)*TAU/
 355      &    ((TAU-TAUR1)**2+GAMR1**2)
 356         ENDIF
 357         IF(NBIN.GE.6) THEN
 358           WTMAT(IBIN,5)=WTMAT(IBIN,5)+(ATAU1/ATAU5)/(TAU+TAUR2)
 359           WTMAT(IBIN,6)=WTMAT(IBIN,6)+(ATAU1/ATAU6)*TAU/
 360      &    ((TAU-TAUR2)**2+GAMR2**2)
 361         ENDIF
 362         IF(NBIN.GT.2+2*MINT(72)) THEN
 363           WTMAT(IBIN,NBIN)=WTMAT(IBIN,NBIN)+(ATAU1/ATAU7)*
 364      &    TAU/MAX(2E-6,1.-TAU)
 365         ENDIF
 366
 367 C...Sum up tau' cross-section pieces in points used.
 368       ELSEIF(IVAR.EQ.2) THEN
 369         TAU=VINTPT(IACC,11)
 370         TAUP=VINTPT(IACC,16)
 371         TAUPMN=VINTPT(IACC,6)
 372         TAUPMX=VINTPT(IACC,26)
 373         ATAUP1=LOG(TAUPMX/TAUPMN)
 374         ATAUP2=((1.-TAU/TAUPMX)**4-(1.-TAU/TAUPMN)**4)/(4.*TAU)
 375         WTMAT(IBIN,1)=WTMAT(IBIN,1)+1.
 376         WTMAT(IBIN,2)=WTMAT(IBIN,2)+(ATAUP1/ATAUP2)*(1.-TAU/TAUP)**3/
 377      &  TAUP
 378         IF(NBIN.GE.3) THEN
 379           ATAUP3=LOG(MAX(2E-6,1.-TAUPMN)/MAX(2E-6,1.-TAUPMX))
 380           WTMAT(IBIN,3)=WTMAT(IBIN,3)+(ATAUP1/ATAUP3)*
 381      &    TAUP/MAX(2E-6,1.-TAUP)
 382         ENDIF
 383
 384 C...Sum up y* cross-section pieces in points used.
 385       ELSEIF(IVAR.EQ.3) THEN
 386         YST=VINTPT(IACC,12)
 387         YSTMIN=VINTPT(IACC,2)
 388         YSTMAX=VINTPT(IACC,22)
 389         AYST0=YSTMAX-YSTMIN
 390         AYST1=0.5*(YSTMAX-YSTMIN)**2
 391         AYST2=AYST1
 392         AYST3=2.*(ATAN(EXP(YSTMAX))-ATAN(EXP(YSTMIN)))
 393         WTMAT(IBIN,1)=WTMAT(IBIN,1)+(AYST0/AYST1)*(YST-YSTMIN)
 394         WTMAT(IBIN,2)=WTMAT(IBIN,2)+(AYST0/AYST2)*(YSTMAX-YST)
 395         WTMAT(IBIN,3)=WTMAT(IBIN,3)+(AYST0/AYST3)/COSH(YST)
 396         IF(MINT(45).EQ.3) THEN
 397           TAUE=VINTPT(IACC,11)
 398           IF(ISTSB.GE.3.AND.ISTSB.LE.5) TAUE=VINTPT(IACC,16)
 399           YST0=-0.5*LOG(TAUE)
 400           AYST4=LOG(MAX(1E-6,EXP(YST0-YSTMIN)-1.)/
 401      &    MAX(1E-6,EXP(YST0-YSTMAX)-1.))
 402           WTMAT(IBIN,4)=WTMAT(IBIN,4)+(AYST0/AYST4)/
 403      &    MAX(1E-6,1.-EXP(YST-YST0))
 404         ENDIF
 405         IF(MINT(46).EQ.3) THEN
 406           TAUE=VINTPT(IACC,11)
 407           IF(ISTSB.GE.3.AND.ISTSB.LE.5) TAUE=VINTPT(IACC,16)
 408           YST0=-0.5*LOG(TAUE)
 409           AYST5=LOG(MAX(1E-6,EXP(YST0+YSTMAX)-1.)/
 410      &    MAX(1E-6,EXP(YST0+YSTMIN)-1.))
 411           WTMAT(IBIN,NBIN)=WTMAT(IBIN,NBIN)+(AYST0/AYST5)/
 412      &    MAX(1E-6,1.-EXP(-YST-YST0))
 413         ENDIF
 414
 415 C...Sum up cos(theta-hat) cross-section pieces in points used.
 416       ELSE
 417         RM34=MAX(1E-20,2.*SQM3*SQM4/(VINTPT(IACC,11)*VINT(2))**2)
 418         RSQM=1.+RM34
 419         CTHMAX=SQRT(1.-4.*VINT(71)**2/(TAUMAX*VINT(2)))
 420         CTHMIN=-CTHMAX
 421         IF(CTHMAX.GT.0.9999) RM34=MAX(RM34,2.*VINT(71)**2/
 422      &  (TAUMAX*VINT(2)))
 423         ACTH1=CTHMAX-CTHMIN
 424         ACTH2=LOG(MAX(RM34,RSQM-CTHMIN)/MAX(RM34,RSQM-CTHMAX))
 425         ACTH3=LOG(MAX(RM34,RSQM+CTHMAX)/MAX(RM34,RSQM+CTHMIN))
 426         ACTH4=1./MAX(RM34,RSQM-CTHMAX)-1./MAX(RM34,RSQM-CTHMIN)
 427         ACTH5=1./MAX(RM34,RSQM+CTHMIN)-1./MAX(RM34,RSQM+CTHMAX)
 428         CTH=VINTPT(IACC,13)
 429         WTMAT(IBIN,1)=WTMAT(IBIN,1)+1.
 430         WTMAT(IBIN,2)=WTMAT(IBIN,2)+(ACTH1/ACTH2)/MAX(RM34,RSQM-CTH)
 431         WTMAT(IBIN,3)=WTMAT(IBIN,3)+(ACTH1/ACTH3)/MAX(RM34,RSQM+CTH)
 432         WTMAT(IBIN,4)=WTMAT(IBIN,4)+(ACTH1/ACTH4)/MAX(RM34,RSQM-CTH)**2
 433         WTMAT(IBIN,5)=WTMAT(IBIN,5)+(ACTH1/ACTH5)/MAX(RM34,RSQM+CTH)**2
 434       ENDIF
 435   170 CONTINUE
 436
 437 C...Check that equation system solvable; else trivial way out.
 438       IF(MSTP(122).GE.2) WRITE(MSTU(11),5400) CVAR(IVAR)
 439       MSOLV=1
 440       WTRELS=0.
 441       DO 180 IBIN=1,NBIN
 442       IF(MSTP(122).GE.2) WRITE(MSTU(11),5500) (WTMAT(IBIN,IRED),
 443      &IRED=1,NBIN),WTREL(IBIN)
 444       IF(NAREL(IBIN).EQ.0) MSOLV=0
 445       WTRELS=WTRELS+WTREL(IBIN)
 446   180 CONTINUE
 447       IF(MSOLV.EQ.0) THEN
 448         DO 190 IBIN=1,NBIN
 449         COEFU(IBIN)=1.
 450         WTRELN(IBIN)=0.1
 451         IF(WTRELS.GT.0.) WTRELN(IBIN)=MAX(0.1,WTREL(IBIN)/WTRELS)
 452   190   CONTINUE
 453
 454 C...Solve to find relative importance of cross-section pieces.
 455       ELSE
 456         DO 200 IBIN=1,NBIN
 457         WTRELN(IBIN)=MAX(0.1,WTREL(IBIN)/WTRELS)
 458   200   CONTINUE
 459         DO 230 IRED=1,NBIN-1
 460         DO 220 IBIN=IRED+1,NBIN
 461         RQT=WTMAT(IBIN,IRED)/WTMAT(IRED,IRED)
 462         WTREL(IBIN)=WTREL(IBIN)-RQT*WTREL(IRED)
 463         DO 210 ICOE=IRED,NBIN
 464         WTMAT(IBIN,ICOE)=WTMAT(IBIN,ICOE)-RQT*WTMAT(IRED,ICOE)
 465   210   CONTINUE
 466   220   CONTINUE
 467   230   CONTINUE
 468         DO 250 IRED=NBIN,1,-1
 469         DO 240 ICOE=IRED+1,NBIN
 470         WTREL(IRED)=WTREL(IRED)-WTMAT(IRED,ICOE)*COEFU(ICOE)
 471   240   CONTINUE
 472         COEFU(IRED)=WTREL(IRED)/WTMAT(IRED,IRED)
 473   250   CONTINUE
 474       ENDIF
 475
 476 C...Normalize coefficients, with piece shared democratically.
 477       COEFSU=0.
 478       WTRELS=0.
 479       DO 260 IBIN=1,NBIN
 480       COEFU(IBIN)=MAX(0.,COEFU(IBIN))
 481       COEFSU=COEFSU+COEFU(IBIN)
 482       WTRELS=WTRELS+WTRELN(IBIN)
 483   260 CONTINUE
 484       IF(COEFSU.GT.0.) THEN
 485         DO 270 IBIN=1,NBIN
 486         COEFO(IBIN)=PARP(122)/NBIN+(1.-PARP(122))*0.5*
 487      &  (COEFU(IBIN)/COEFSU+WTRELN(IBIN)/WTRELS)
 488   270   CONTINUE
 489       ELSE
 490         DO 280 IBIN=1,NBIN
 491         COEFO(IBIN)=1./NBIN
 492   280   CONTINUE
 493       ENDIF
 494       IF(IVAR.EQ.1) IOFF=0
 495       IF(IVAR.EQ.2) IOFF=17
 496       IF(IVAR.EQ.3) IOFF=7
 497       IF(IVAR.EQ.4) IOFF=12
 498       DO 290 IBIN=1,NBIN
 499       ICOF=IOFF+IBIN
 500       IF(IVAR.EQ.1.AND.IBIN.GT.2+2*MINT(72)) ICOF=7
 501       IF(IVAR.EQ.3.AND.IBIN.EQ.4.AND.MINT(45).NE.3) ICOF=ICOF+1
 502       COEF(ISUB,ICOF)=COEFO(IBIN)
 503   290 CONTINUE
 504       IF(MSTP(122).GE.2) WRITE(MSTU(11),5600) CVAR(IVAR),
 505      &(COEFO(IBIN),IBIN=1,NBIN)
 506   300 CONTINUE
 507
 508 C...Find two most promising maxima among points previously determined.
 509       DO 310 J=1,4
 510       IACCMX(J)=0
 511       SIGSMX(J)=0.
 512   310 CONTINUE
 513       NMAX=0
 514       DO 370 IACC=1,NACC
 515       DO 320 J=1,30
 516       VINT(10+J)=VINTPT(IACC,J)
 517   320 CONTINUE
 518       IF(ISTSB.NE.5) THEN
 519         CALL PYSIGH(NCHN,SIGS)
 520         IF(MWTXS.EQ.1) THEN
 521           CALL PYEVWT(WTXS)
 522           SIGS=WTXS*SIGS
 523         ENDIF
 524       ELSE
 525         SIGS=0.
 526         DO 330 IKIN3=1,MSTP(129)
 527         CALL PYKMAP(5,0,0.)
 528         IF(MINT(51).EQ.1) GOTO 330
 529         CALL PYSIGH(NCHN,SIGTMP)
 530         IF(MWTXS.EQ.1) THEN
 531           CALL PYEVWT(WTXS)
 532           SIGTMP=WTXS*SIGTMP
 533         ENDIF
 534         IF(SIGTMP.GT.SIGS) SIGS=SIGTMP
 535   330   CONTINUE
 536       ENDIF
 537       IEQ=0
 538       DO 340 IMV=1,NMAX
 539       IF(ABS(SIGS-SIGSMX(IMV)).LT.1E-4*(SIGS+SIGSMX(IMV))) IEQ=IMV
 540   340 CONTINUE
 541       IF(IEQ.EQ.0) THEN
 542         DO 350 IMV=NMAX,1,-1
 543         IIN=IMV+1
 544         IF(SIGS.LE.SIGSMX(IMV)) GOTO 360
 545         IACCMX(IMV+1)=IACCMX(IMV)
 546         SIGSMX(IMV+1)=SIGSMX(IMV)
 547   350   CONTINUE
 548         IIN=1
 549   360   IACCMX(IIN)=IACC
 550         SIGSMX(IIN)=SIGS
 551         IF(NMAX.LE.1) NMAX=NMAX+1
 552       ENDIF
 553   370 CONTINUE
 554
 555 C...Read out starting position for search.
 556       IF(MSTP(122).GE.2) WRITE(MSTU(11),5700)
 557       SIGSAM=SIGSMX(1)
 558       DO 420 IMAX=1,NMAX
 559       IACC=IACCMX(IMAX)
 560       MTAU=MVARPT(IACC,1)
 561       MTAUP=MVARPT(IACC,2)
 562       MYST=MVARPT(IACC,3)
 563       MCTH=MVARPT(IACC,4)
 564       VTAU=0.5
 565       VYST=0.5
 566       VCTH=0.5
 567       VTAUP=0.5
 568
 569 C...Starting point and step size in parameter space.
 570       DO 410 IRPT=1,2
 571       DO 400 IVAR=1,4
 572       IF(NPTS(IVAR).EQ.1) GOTO 400
 573       IF(IVAR.EQ.1) VVAR=VTAU
 574       IF(IVAR.EQ.2) VVAR=VTAUP
 575       IF(IVAR.EQ.3) VVAR=VYST
 576       IF(IVAR.EQ.4) VVAR=VCTH
 577       IF(IVAR.EQ.1) MVAR=MTAU
 578       IF(IVAR.EQ.2) MVAR=MTAUP
 579       IF(IVAR.EQ.3) MVAR=MYST
 580       IF(IVAR.EQ.4) MVAR=MCTH
 581       IF(IRPT.EQ.1) VDEL=0.1
 582       IF(IRPT.EQ.2) VDEL=MAX(0.01,MIN(0.05,VVAR-0.02,0.98-VVAR))
 583       IF(IRPT.EQ.1) VMAR=0.02
 584       IF(IRPT.EQ.2) VMAR=0.002
 585       IMOV0=1
 586       IF(IRPT.EQ.1.AND.IVAR.EQ.1) IMOV0=0
 587       DO 390 IMOV=IMOV0,8
 588
 589 C...Define new point in parameter space.
 590       IF(IMOV.EQ.0) THEN
 591         INEW=2
 592         VNEW=VVAR
 593       ELSEIF(IMOV.EQ.1) THEN
 594         INEW=3
 595         VNEW=VVAR+VDEL
 596       ELSEIF(IMOV.EQ.2) THEN
 597         INEW=1
 598         VNEW=VVAR-VDEL
 599       ELSEIF(SIGSSM(3).GE.MAX(SIGSSM(1),SIGSSM(2)).AND.
 600      &VVAR+2.*VDEL.LT.1.-VMAR) THEN
 601         VVAR=VVAR+VDEL
 602         SIGSSM(1)=SIGSSM(2)
 603         SIGSSM(2)=SIGSSM(3)
 604         INEW=3
 605         VNEW=VVAR+VDEL
 606       ELSEIF(SIGSSM(1).GE.MAX(SIGSSM(2),SIGSSM(3)).AND.
 607      &VVAR-2.*VDEL.GT.VMAR) THEN
 608         VVAR=VVAR-VDEL
 609         SIGSSM(3)=SIGSSM(2)
 610         SIGSSM(2)=SIGSSM(1)
 611         INEW=1
 612         VNEW=VVAR-VDEL
 613       ELSEIF(SIGSSM(3).GE.SIGSSM(1)) THEN
 614         VDEL=0.5*VDEL
 615         VVAR=VVAR+VDEL
 616         SIGSSM(1)=SIGSSM(2)
 617         INEW=2
 618         VNEW=VVAR
 619       ELSE
 620         VDEL=0.5*VDEL
 621         VVAR=VVAR-VDEL
 622         SIGSSM(3)=SIGSSM(2)
 623         INEW=2
 624         VNEW=VVAR
 625       ENDIF
 626
 627 C...Convert to relevant variables and find derived new limits.
 628       IF(IVAR.EQ.1) THEN
 629         VTAU=VNEW
 630         CALL PYKMAP(1,MTAU,VTAU)
 631         IF(ISTSB.GE.3.AND.ISTSB.LE.5) CALL PYKLIM(4)
 632       ENDIF
 633       IF(IVAR.LE.2.AND.ISTSB.GE.3.AND.ISTSB.LE.5) THEN
 634         IF(IVAR.EQ.2) VTAUP=VNEW
 635         CALL PYKMAP(4,MTAUP,VTAUP)
 636       ENDIF
 637       IF(IVAR.LE.2) CALL PYKLIM(2)
 638       IF(IVAR.LE.3) THEN
 639         IF(IVAR.EQ.3) VYST=VNEW
 640         CALL PYKMAP(2,MYST,VYST)
 641         CALL PYKLIM(3)
 642       ENDIF
 643       IF(ISTSB.EQ.2.OR.ISTSB.EQ.4.OR.ISTSB.EQ.6) THEN
 644         IF(IVAR.EQ.4) VCTH=VNEW
 645         CALL PYKMAP(3,MCTH,VCTH)
 646       ENDIF
 647       IF(ISUB.EQ.96) VINT(25)=VINT(21)*(1.-VINT(23)**2)
 648
 649 C...Evaluate cross-section. Save new maximum. Final maximum.
 650       IF(ISTSB.NE.5) THEN
 651         CALL PYSIGH(NCHN,SIGS)
 652         IF(MWTXS.EQ.1) THEN
 653           CALL PYEVWT(WTXS)
 654           SIGS=WTXS*SIGS
 655         ENDIF
 656       ELSE
 657         SIGS=0.
 658         DO 380 IKIN3=1,MSTP(129)
 659         CALL PYKMAP(5,0,0.)
 660         IF(MINT(51).EQ.1) GOTO 380
 661         CALL PYSIGH(NCHN,SIGTMP)
 662         IF(MWTXS.EQ.1) THEN
 663           CALL PYEVWT(WTXS)
 664           SIGTMP=WTXS*SIGTMP
 665         ENDIF
 666         IF(SIGTMP.GT.SIGS) SIGS=SIGTMP
 667   380   CONTINUE
 668       ENDIF
 669       SIGSSM(INEW)=SIGS
 670       IF(SIGS.GT.SIGSAM) SIGSAM=SIGS
 671       IF(MSTP(122).GE.2) WRITE(MSTU(11),5800) IMAX,IVAR,MVAR,IMOV,
 672      &VNEW,VINT(21),VINT(22),VINT(23),VINT(26),SIGS
 673   390 CONTINUE
 674   400 CONTINUE
 675   410 CONTINUE
 676   420 CONTINUE
 677       IF(MSTP(121).EQ.1) SIGSAM=PARP(121)*SIGSAM
 678       XSEC(ISUB,1)=1.05*SIGSAM
 679   430 CONTINUE
 680       IF(MSTP(173).EQ.1.AND.ISUB.NE.96) XSEC(ISUB,1)=
 681      &PARP(174)*XSEC(ISUB,1)
 682       IF(ISUB.NE.96) XSEC(0,1)=XSEC(0,1)+XSEC(ISUB,1)
 683   440 CONTINUE
 684       MINT(51)=0
 685
 686 C...Print summary table.
 687       IF(NPOSI.EQ.0) THEN
 688         WRITE(MSTU(11),5900)
 689         STOP
 690       ENDIF
 691       IF(MSTP(122).GE.1) THEN
 692         WRITE(MSTU(11),6000)
 693         WRITE(MSTU(11),6100)
 694         DO 450 ISUB=1,200
 695         IF(MSUB(ISUB).NE.1.AND.ISUB.NE.96) GOTO 450
 696         IF(ISUB.EQ.96.AND.MINT(50).EQ.0) GOTO 450
 697         IF(ISUB.EQ.96.AND.MSUB(95).NE.1.AND.MSTP(81).LE.0) GOTO 450
 698         IF(ISUB.EQ.96.AND.MINT(49).EQ.0.AND.MSTP(131).EQ.0) GOTO 450
 699         IF(MSUB(95).EQ.1.AND.(ISUB.EQ.11.OR.ISUB.EQ.12.OR.ISUB.EQ.13.OR.
 700      &  ISUB.EQ.28.OR.ISUB.EQ.53.OR.ISUB.EQ.68)) GOTO 450
 701         WRITE(MSTU(11),6200) ISUB,PROC(ISUB),XSEC(ISUB,1)
 702   450   CONTINUE
 703         WRITE(MSTU(11),6300)
 704       ENDIF
 705
 706 C...Format statements for maximization results.
 707  5000 FORMAT(/1X,'Coefficient optimization and maximum search for ',
 708      &'subprocess no',I4/1X,'Coefficient modes     tau',10X,'y*',9X,
 709      &'cth',9X,'tau''',7X,'sigma')
 710  5100 FORMAT(1X,'Warning: requested subprocess ',I3,' has no allowed ',
 711      &'phase space.'/1X,'Process switched off!')
 712  5200 FORMAT(1X,4I4,F12.8,F12.6,F12.7,F12.8,1P,E12.4)
 713  5300 FORMAT(1X,'Warning: requested subprocess ',I3,' has vanishing ',
 714      &'cross-section.'/1X,'Process switched off!')
 715  5400 FORMAT(1X,'Coefficients of equation system to be solved for ',A4)
 716  5500 FORMAT(1X,1P,8E11.3)
 717  5600 FORMAT(1X,'Result for ',A4,':',7F9.4)
 718  5700 FORMAT(1X,'Maximum search for given coefficients'/2X,'MAX VAR ',
 719      &'MOD MOV   VNEW',7X,'tau',7X,'y*',8X,'cth',7X,'tau''',7X,'sigma')
 720  5800 FORMAT(1X,4I4,F8.4,F11.7,F9.3,F11.6,F11.7,1P,E12.4)
 721  5900 FORMAT(1X,'Error: no requested process has non-vanishing ',
 722      &'cross-section.'/1X,'Execution stopped!')
 723  6000 FORMAT(/1X,8('*'),1X,'PYMAXI: summary of differential ',
 724      &'cross-section maximum search',1X,8('*'))
 725  6100 FORMAT(/11X,58('=')/11X,'I',38X,'I',17X,'I'/11X,'I  ISUB  ',
 726      &'Subprocess name',15X,'I  Maximum value  I'/11X,'I',38X,'I',
 727      &17X,'I'/11X,58('=')/11X,'I',38X,'I',17X,'I')
 728  6200 FORMAT(11X,'I',2X,I3,3X,A28,2X,'I',2X,1P,E12.4,3X,'I')
 729  6300 FORMAT(11X,'I',38X,'I',17X,'I'/11X,58('='))
 730
 731       RETURN
 732       END