Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pymaxi.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pymaxi.f
1 
2 C*********************************************************************
3 
4 C...PYMAXI
5 C...Finds optimal set of coefficients for kinematical variable selection
6 C...and the maximum of the part of the differential cross-section used
7 C...in the event weighting.
8 
9  SUBROUTINE pymaxi
10 
11 C...Double precision and integer declarations.
12  IMPLICIT DOUBLE PRECISION(a-h, o-z)
13  IMPLICIT INTEGER(i-n)
14  INTEGER pyk,pychge,pycomp
15 C...Parameter statement to help give large particle numbers.
16  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
17  &kexcit=4000000,kdimen=5000000)
18 
19 C...User process initialization commonblock.
20  INTEGER maxpup
21  parameter(maxpup=100)
22  INTEGER idbmup,pdfgup,pdfsup,idwtup,nprup,lprup
23  DOUBLE PRECISION ebmup,xsecup,xerrup,xmaxup
24  common/heprup/idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
25  &idwtup,nprup,xsecup(maxpup),xerrup(maxpup),xmaxup(maxpup),
26  &lprup(maxpup)
27  SAVE /heprup/
28 
29 C...Commonblocks.
30  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
31  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
32  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
33  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
34  common/pypars/mstp(200),parp(200),msti(200),pari(200)
35  common/pyint1/mint(400),vint(400)
36  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
37  common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
38  common/pyint4/mwid(500),wids(500,5)
39  common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
40  common/pyint6/proc(0:500)
41  CHARACTER proc*28
42  common/pyint7/sigt(0:6,0:6,0:5)
43  SAVE /pydat1/,/pydat2/,/pydat3/,/pysubs/,/pypars/,/pyint1/,
44  &/pyint2/,/pyint3/,/pyint4/,/pyint5/,/pyint6/,/pyint7/
45 C...Local arrays, character variables and data.
46  CHARACTER cvar(4)*4
47  dimension npts(4),mvarpt(500,4),vintpt(500,30),sigspt(500),
48  &narel(7),wtrel(7),wtmat(7,7),wtreln(7),coefu(7),coefo(7),
49  &iaccmx(4),sigsmx(4),sigssm(3),pmmn(2)
50  DATA cvar/'tau ','tau''','y* ','cth '/
51  DATA sigssm/3*0d0/
52 
53 C...Initial values and loop over subprocesses.
54  nposi=0
55  vint(143)=1d0
56  vint(144)=1d0
57  xsec(0,1)=0d0
58  DO 460 isub=1,500
59  mint(1)=isub
60  mint(51)=0
61 
62 C...Find maximum weight factors for photon flux.
63  IF(msub(isub).EQ.1.OR.(isub.GE.91.AND.isub.LE.100)) THEN
64  IF(mint(141).NE.0.OR.mint(142).NE.0) CALL pygaga(2,wtgaga)
65  ENDIF
66 
67 C...Select subprocess to study: skip cases not applicable.
68  IF(iset(isub).EQ.11) THEN
69  IF(msub(isub).NE.1) goto 460
70 C...User process intialization: cross section model dependent.
71  IF(iabs(idwtup).EQ.1) THEN
72  IF(idwtup.GT.0.AND.xmaxup(kfpr(isub,1)).LT.0d0) CALL
73  & pyerrm(26,'(PYMAXI:) Negative XMAXUP for user process')
74  xsec(isub,1)=1.00000001d-9*abs(xmaxup(kfpr(isub,1)))
75  ELSE
76  IF((idwtup.EQ.2.OR.idwtup.EQ.3).AND.
77  & xsecup(kfpr(isub,1)).LT.0d0) CALL
78  & pyerrm(26,'(PYMAXI:) Negative XSECUP for user process')
79  IF(idwtup.EQ.2.AND.xmaxup(kfpr(isub,1)).LT.0d0) CALL
80  & pyerrm(26,'(PYMAXI:) Negative XMAXUP for user process')
81  xsec(isub,1)=1.00000001d-9*abs(xsecup(kfpr(isub,1)))
82  ENDIF
83  IF(mint(141).NE.0.OR.mint(142).NE.0) xsec(isub,1)=
84  & wtgaga*xsec(isub,1)
85  nposi=nposi+1
86  goto 450
87  ELSEIF(isub.GE.91.AND.isub.LE.95) THEN
88  CALL pysigh(nchn,sigs)
89  xsec(isub,1)=sigs
90  IF(mint(141).NE.0.OR.mint(142).NE.0) xsec(isub,1)=
91  & wtgaga*xsec(isub,1)
92  IF(msub(isub).NE.1) goto 460
93  nposi=nposi+1
94  goto 450
95  ELSEIF(isub.EQ.99.AND.msub(isub).EQ.1) THEN
96  CALL pysigh(nchn,sigs)
97  xsec(isub,1)=sigs
98  IF(mint(141).NE.0.OR.mint(142).NE.0) xsec(isub,1)=
99  & wtgaga*xsec(isub,1)
100  IF(xsec(isub,1).EQ.0d0) THEN
101  msub(isub)=0
102  ELSE
103  nposi=nposi+1
104  ENDIF
105  goto 450
106  ELSEIF(isub.EQ.96) THEN
107  IF(mint(50).EQ.0) goto 460
108  IF(msub(95).NE.1.AND.mod(mstp(81),10).LE.0.AND.mstp(131).LE.0)
109  & goto 460
110  IF(mint(49).EQ.0.AND.mstp(131).EQ.0) goto 460
111  ELSEIF(isub.EQ.11.OR.isub.EQ.12.OR.isub.EQ.13.OR.isub.EQ.28.OR.
112  & isub.EQ.53.OR.isub.EQ.68) THEN
113  IF(msub(isub).NE.1.OR.msub(95).EQ.1) goto 460
114  ELSEIF(isub.GE.381.AND.isub.LE.386) THEN
115  IF(msub(isub).NE.1.OR.msub(95).EQ.1) goto 460
116  ELSE
117  IF(msub(isub).NE.1) goto 460
118  ENDIF
119  istsb=iset(isub)
120  IF(isub.EQ.96) istsb=2
121  IF(mstp(122).GE.2) WRITE(mstu(11),5000) isub
122  mwtxs=0
123  IF(mstp(142).GE.1.AND.isub.NE.96.AND.msub(91)+msub(92)+msub(93)+
124  & msub(94)+msub(95).EQ.0) mwtxs=1
125 
126 C...Find resonances (explicit or implicit in cross-section).
127  mint(72)=0
128  kfr1=0
129  IF(istsb.EQ.1.OR.istsb.EQ.3.OR.istsb.EQ.5) THEN
130  kfr1=kfpr(isub,1)
131  ELSEIF(isub.EQ.24.OR.isub.EQ.25.OR.isub.EQ.110.OR.isub.EQ.165
132  & .OR.isub.EQ.171.OR.isub.EQ.176) THEN
133  kfr1=23
134  ELSEIF(isub.EQ.23.OR.isub.EQ.26.OR.isub.EQ.166.OR.isub.EQ.172
135  & .OR.isub.EQ.177) THEN
136  kfr1=24
137  ELSEIF(isub.GE.71.AND.isub.LE.77) THEN
138  kfr1=25
139  IF(mstp(46).EQ.5) THEN
140  kfr1=89
141  pmas(89,1)=parp(45)
142  pmas(89,2)=parp(45)**3/(96d0*paru(1)*parp(47)**2)
143  ENDIF
144  ELSEIF(isub.EQ.194) THEN
145  kfr1=ktechn+113
146  ELSEIF(isub.EQ.195) THEN
147  kfr1=ktechn+213
148  ELSEIF(isub.GE.361.AND.isub.LE.368) THEN
149  kfr1=ktechn+113
150  ELSEIF(isub.GE.370.AND.isub.LE.377) THEN
151  kfr1=ktechn+213
152  ENDIF
153  ckmx=ckin(2)
154  IF(ckmx.LE.0d0) ckmx=vint(1)
155  kcr1=pycomp(kfr1)
156  IF(kfr1.NE.0) THEN
157  IF(ckin(1).GT.pmas(kcr1,1)+20d0*pmas(kcr1,2).OR.
158  & ckmx.LT.pmas(kcr1,1)-20d0*pmas(kcr1,2)) kfr1=0
159  ENDIF
160  IF(kfr1.NE.0) THEN
161  taur1=pmas(kcr1,1)**2/vint(2)
162  IF(kfr1.EQ.ktechn+113) THEN
163  CALL pytecm(s1,s2)
164  taur1=s1/vint(2)
165  ENDIF
166  gamr1=pmas(kcr1,1)*pmas(kcr1,2)/vint(2)
167  mint(72)=1
168  mint(73)=kfr1
169  vint(73)=taur1
170  vint(74)=gamr1
171  ENDIF
172  kfr2=0
173  IF(isub.EQ.141.OR.isub.EQ.194.OR.(isub.GE.364.AND.isub.LE.368))
174  $ THEN
175  kfr2=23
176  IF(isub.EQ.194) THEN
177  kfr2=ktechn+223
178  ELSEIF(isub.GE.364.AND.isub.LE.368) THEN
179  kfr2=ktechn+223
180  ENDIF
181  kcr2=pycomp(kfr2)
182  taur2=pmas(kcr2,1)**2/vint(2)
183  IF(kfr2.EQ.ktechn+223) THEN
184  CALL pytecm(s1,s2)
185  taur2=s2/vint(2)
186  ENDIF
187  gamr2=pmas(kcr2,1)*pmas(kcr2,2)/vint(2)
188  IF(ckin(1).GT.pmas(kcr2,1)+20d0*pmas(kcr2,2).OR.
189  & ckmx.LT.pmas(kcr2,1)-20d0*pmas(kcr2,2)) kfr2=0
190  IF(kfr2.NE.0.AND.kfr1.NE.0) THEN
191  mint(72)=2
192  mint(74)=kfr2
193  vint(75)=taur2
194  vint(76)=gamr2
195  ELSEIF(kfr2.NE.0) THEN
196  kfr1=kfr2
197  taur1=taur2
198  gamr1=gamr2
199  mint(72)=1
200  mint(73)=kfr1
201  vint(73)=taur1
202  vint(74)=gamr1
203  kfr2=0
204  ENDIF
205  ENDIF
206 
207 C...Find product masses and minimum pT of process.
208  sqm3=0d0
209  sqm4=0d0
210  mint(71)=0
211  vint(71)=ckin(3)
212  vint(80)=1d0
213  IF(istsb.EQ.2.OR.istsb.EQ.4) THEN
214  nbw=0
215  DO 110 i=1,2
216  pmmn(i)=0d0
217  IF(kfpr(isub,i).EQ.0) THEN
218  ELSEIF(mstp(42).LE.0.OR.pmas(pycomp(kfpr(isub,i)),2).LT.
219  & parp(41)) THEN
220  IF(i.EQ.1) sqm3=pmas(pycomp(kfpr(isub,i)),1)**2
221  IF(i.EQ.2) sqm4=pmas(pycomp(kfpr(isub,i)),1)**2
222  ELSE
223  nbw=nbw+1
224 C...This prevents SUSY/t particles from becoming too light.
225  kflw=kfpr(isub,i)
226  IF(kflw/ksusy1.EQ.1.OR.kflw/ksusy1.EQ.2) THEN
227  kcw=pycomp(kflw)
228  pmmn(i)=pmas(kcw,1)
229  DO 100 idc=mdcy(kcw,2),mdcy(kcw,2)+mdcy(kcw,3)-1
230  IF(mdme(idc,1).GT.0.AND.brat(idc).GT.1e-4) THEN
231  pmsum=pmas(pycomp(kfdp(idc,1)),1)+
232  & pmas(pycomp(kfdp(idc,2)),1)
233  IF(kfdp(idc,3).NE.0) pmsum=pmsum+
234  & pmas(pycomp(kfdp(idc,3)),1)
235  pmmn(i)=min(pmmn(i),pmsum)
236  ENDIF
237  100 CONTINUE
238  ELSEIF(kflw.EQ.6) THEN
239  pmmn(i)=pmas(24,1)+pmas(5,1)
240  ENDIF
241  ENDIF
242  110 CONTINUE
243  IF(nbw.GE.1) THEN
244  ckin41=ckin(41)
245  ckin43=ckin(43)
246  ckin(41)=max(pmmn(1),ckin(41))
247  ckin(43)=max(pmmn(2),ckin(43))
248  CALL pyofsh(3,0,kfpr(isub,1),kfpr(isub,2),0d0,pqm3,pqm4)
249  ckin(41)=ckin41
250  ckin(43)=ckin43
251  IF(mint(51).EQ.1) THEN
252  WRITE(mstu(11),5100) isub
253  msub(isub)=0
254  goto 460
255  ENDIF
256  sqm3=pqm3**2
257  sqm4=pqm4**2
258  ENDIF
259  IF(min(sqm3,sqm4).LT.ckin(6)**2) mint(71)=1
260  IF(mint(71).EQ.1) vint(71)=max(ckin(3),ckin(5))
261  IF(isub.EQ.96.AND.mstp(82).LE.1) THEN
262  vint(71)=parp(81)*(vint(1)/parp(89))**parp(90)
263  ELSEIF(isub.EQ.96) THEN
264  vint(71)=0.08d0*parp(82)*(vint(1)/parp(89))**parp(90)
265  ENDIF
266  ENDIF
267  vint(63)=sqm3
268  vint(64)=sqm4
269 
270 C...Prepare for additional variable choices in 2 -> 3.
271  IF(istsb.EQ.5) THEN
272  vint(201)=0d0
273  IF(kfpr(isub,2).GT.0) vint(201)=pmas(pycomp(kfpr(isub,2)),1)
274  vint(206)=vint(201)
275  IF(isub.EQ.401.OR.isub.EQ.402) vint(206)=pmas(5,1)
276  vint(204)=pmas(23,1)
277  IF(isub.EQ.124.OR.isub.EQ.174.OR.isub.EQ.179.OR.isub.EQ.351)
278  & vint(204)=pmas(24,1)
279  IF(isub.EQ.352) vint(204)=pmas(pycomp(9900024),1)
280  IF(isub.EQ.121.OR.isub.EQ.122.OR.isub.EQ.181.OR.isub.EQ.182
281  & .OR.isub.EQ.186.OR.isub.EQ.187.OR.isub.EQ.401.OR.isub.EQ.402)
282  & vint(204)=vint(201)
283  vint(209)=vint(204)
284  IF(isub.EQ.401.OR.isub.EQ.402) vint(209)=vint(206)
285  ENDIF
286 
287 C...Number of points for each variable: tau, tau', y*, cos(theta-hat).
288  npts(1)=2+2*mint(72)
289  IF(mint(47).EQ.1) THEN
290  IF(istsb.EQ.1.OR.istsb.EQ.2) npts(1)=1
291  ELSEIF(mint(47).GE.5) THEN
292  IF(istsb.LE.2.OR.istsb.GT.5) npts(1)=npts(1)+1
293  ENDIF
294  npts(2)=1
295  IF(istsb.GE.3.AND.istsb.LE.5) THEN
296  IF(mint(47).GE.2) npts(2)=2
297  IF(mint(47).GE.5) npts(2)=3
298  ENDIF
299  npts(3)=1
300  IF(mint(47).EQ.4.OR.mint(47).EQ.5) THEN
301  npts(3)=3
302  IF(mint(45).EQ.3) npts(3)=npts(3)+1
303  IF(mint(46).EQ.3) npts(3)=npts(3)+1
304  ENDIF
305  npts(4)=1
306  IF(istsb.EQ.2.OR.istsb.EQ.4) npts(4)=5
307  ntry=npts(1)*npts(2)*npts(3)*npts(4)
308 
309 C...Reset coefficients of cross-section weighting.
310  DO 120 j=1,20
311  coef(isub,j)=0d0
312  120 CONTINUE
313  coef(isub,1)=1d0
314  coef(isub,8)=0.5d0
315  coef(isub,9)=0.5d0
316  coef(isub,13)=1d0
317  coef(isub,18)=1d0
318  mcth=0
319  mtaup=0
320  metaup=0
321  vint(23)=0d0
322  vint(26)=0d0
323  sigsam=0d0
324 
325 C...Find limits and select tau, y*, cos(theta-hat) and tau' values,
326 C...in grid of phase space points.
327  CALL pyklim(1)
328  metau=mint(51)
329  nacc=0
330  DO 150 itry=1,ntry
331  mint(51)=0
332  IF(metau.EQ.1) goto 150
333  IF(mod(itry-1,npts(2)*npts(3)*npts(4)).EQ.0) THEN
334  mtau=1+(itry-1)/(npts(2)*npts(3)*npts(4))
335  IF(mtau.GT.2+2*mint(72)) mtau=7
336  rtau=0.5d0
337 C...Special case when both resonances have same mass,
338 C...as is often the case in process 194.
339  IF(mint(72).EQ.2) THEN
340  IF(abs(pmas(kcr2,1)-pmas(kcr1,1)).LT.
341  & 0.01d0*(pmas(kcr2,1)+pmas(kcr1,1))) THEN
342  IF(mtau.EQ.3.OR.mtau.EQ.4) THEN
343  rtau=0.4d0
344  ELSEIF(mtau.EQ.5.OR.mtau.EQ.6) THEN
345  rtau=0.6d0
346  ENDIF
347  ENDIF
348  ENDIF
349  CALL pykmap(1,mtau,rtau)
350  IF(istsb.GE.3.AND.istsb.LE.5) CALL pyklim(4)
351  metaup=mint(51)
352  ENDIF
353  IF(metaup.EQ.1) goto 150
354  IF(istsb.GE.3.AND.istsb.LE.5.AND.mod(itry-1,npts(3)*npts(4))
355  & .EQ.0) THEN
356  mtaup=1+mod((itry-1)/(npts(3)*npts(4)),npts(2))
357  CALL pykmap(4,mtaup,0.5d0)
358  ENDIF
359  IF(mod(itry-1,npts(3)*npts(4)).EQ.0) THEN
360  CALL pyklim(2)
361  meyst=mint(51)
362  ENDIF
363  IF(meyst.EQ.1) goto 150
364  IF(mod(itry-1,npts(4)).EQ.0) THEN
365  myst=1+mod((itry-1)/npts(4),npts(3))
366  IF(myst.EQ.4.AND.mint(45).NE.3) myst=5
367  CALL pykmap(2,myst,0.5d0)
368  CALL pyklim(3)
369  mecth=mint(51)
370  ENDIF
371  IF(mecth.EQ.1) goto 150
372  IF(istsb.EQ.2.OR.istsb.EQ.4) THEN
373  mcth=1+mod(itry-1,npts(4))
374  CALL pykmap(3,mcth,0.5d0)
375  ENDIF
376  IF(isub.EQ.96) vint(25)=vint(21)*(1d0-vint(23)**2)
377 
378 C...Store position and limits.
379  mint(51)=0
380  CALL pyklim(0)
381  IF(mint(51).EQ.1) goto 150
382  nacc=nacc+1
383  mvarpt(nacc,1)=mtau
384  mvarpt(nacc,2)=mtaup
385  mvarpt(nacc,3)=myst
386  mvarpt(nacc,4)=mcth
387  DO 130 j=1,30
388  vintpt(nacc,j)=vint(10+j)
389  130 CONTINUE
390 
391 C...Normal case: calculate cross-section.
392  IF(istsb.NE.5) THEN
393  CALL pysigh(nchn,sigs)
394  IF(mwtxs.EQ.1) THEN
395  CALL pyevwt(wtxs)
396  sigs=wtxs*sigs
397  ENDIF
398 
399 C..2 -> 3: find highest value out of a number of tries.
400  ELSE
401  sigs=0d0
402  DO 140 ikin3=1,mstp(129)
403  CALL pykmap(5,0,0d0)
404  IF(mint(51).EQ.1) goto 140
405  CALL pysigh(nchn,sigtmp)
406  IF(mwtxs.EQ.1) THEN
407  CALL pyevwt(wtxs)
408  sigtmp=wtxs*sigtmp
409  ENDIF
410  IF(sigtmp.GT.sigs) sigs=sigtmp
411  140 CONTINUE
412  ENDIF
413 
414 C...Store cross-section.
415  sigspt(nacc)=sigs
416  IF(sigs.GT.sigsam) sigsam=sigs
417  IF(mstp(122).GE.2) WRITE(mstu(11),5200) mtau,myst,mcth,mtaup,
418  & vint(21),vint(22),vint(23),vint(26),sigs
419  150 CONTINUE
420  IF(nacc.EQ.0) THEN
421  WRITE(mstu(11),5100) isub
422  msub(isub)=0
423  goto 460
424  ELSEIF(sigsam.EQ.0d0) THEN
425  WRITE(mstu(11),5300) isub
426  msub(isub)=0
427  goto 460
428  ENDIF
429  IF(isub.NE.96) nposi=nposi+1
430 
431 C...Calculate integrals in tau over maximal phase space limits.
432  taumin=vint(11)
433  taumax=vint(31)
434  atau1=log(taumax/taumin)
435  IF(npts(1).GE.2) THEN
436  atau2=(taumax-taumin)/(taumax*taumin)
437  ENDIF
438  IF(npts(1).GE.4) THEN
439  atau3=log(taumax/taumin*(taumin+taur1)/(taumax+taur1))/taur1
440  atau4=(atan((taumax-taur1)/gamr1)-atan((taumin-taur1)/gamr1))/
441  & gamr1
442  ENDIF
443  IF(npts(1).GE.6) THEN
444  atau5=log(taumax/taumin*(taumin+taur2)/(taumax+taur2))/taur2
445  atau6=(atan((taumax-taur2)/gamr2)-atan((taumin-taur2)/gamr2))/
446  & gamr2
447  ENDIF
448  IF(npts(1).GT.2+2*mint(72)) THEN
449  atau7=log(max(2d-10,1d0-taumin)/max(2d-10,1d0-taumax))
450  ENDIF
451 
452 C...Reset. Sum up cross-sections in points calculated.
453  DO 320 ivar=1,4
454  IF(npts(ivar).EQ.1) goto 320
455  IF(isub.EQ.96.AND.ivar.EQ.4) goto 320
456  nbin=npts(ivar)
457  DO 170 j1=1,nbin
458  narel(j1)=0
459  wtrel(j1)=0d0
460  coefu(j1)=0d0
461  DO 160 j2=1,nbin
462  wtmat(j1,j2)=0d0
463  160 CONTINUE
464  170 CONTINUE
465  DO 180 iacc=1,nacc
466  ibin=mvarpt(iacc,ivar)
467  IF(ivar.EQ.1.AND.ibin.EQ.7) ibin=3+2*mint(72)
468  IF(ivar.EQ.3.AND.ibin.EQ.5.AND.mint(45).NE.3) ibin=4
469  narel(ibin)=narel(ibin)+1
470  wtrel(ibin)=wtrel(ibin)+sigspt(iacc)
471 
472 C...Sum up tau cross-section pieces in points used.
473  IF(ivar.EQ.1) THEN
474  tau=vintpt(iacc,11)
475  wtmat(ibin,1)=wtmat(ibin,1)+1d0
476  wtmat(ibin,2)=wtmat(ibin,2)+(atau1/atau2)/tau
477  IF(nbin.GE.4) THEN
478  wtmat(ibin,3)=wtmat(ibin,3)+(atau1/atau3)/(tau+taur1)
479  wtmat(ibin,4)=wtmat(ibin,4)+(atau1/atau4)*tau/
480  & ((tau-taur1)**2+gamr1**2)
481  ENDIF
482  IF(nbin.GE.6) THEN
483  wtmat(ibin,5)=wtmat(ibin,5)+(atau1/atau5)/(tau+taur2)
484  wtmat(ibin,6)=wtmat(ibin,6)+(atau1/atau6)*tau/
485  & ((tau-taur2)**2+gamr2**2)
486  ENDIF
487  IF(nbin.GT.2+2*mint(72)) THEN
488  wtmat(ibin,nbin)=wtmat(ibin,nbin)+(atau1/atau7)*
489  & tau/max(2d-10,1d0-tau)
490  ENDIF
491 
492 C...Sum up tau' cross-section pieces in points used.
493  ELSEIF(ivar.EQ.2) THEN
494  tau=vintpt(iacc,11)
495  taup=vintpt(iacc,16)
496  taupmn=vintpt(iacc,6)
497  taupmx=vintpt(iacc,26)
498  ataup1=log(taupmx/taupmn)
499  ataup2=((1d0-tau/taupmx)**4-(1d0-tau/taupmn)**4)/(4d0*tau)
500  wtmat(ibin,1)=wtmat(ibin,1)+1d0
501  wtmat(ibin,2)=wtmat(ibin,2)+(ataup1/ataup2)*
502  & (1d0-tau/taup)**3/taup
503  IF(nbin.GE.3) THEN
504  ataup3=log(max(2d-10,1d0-taupmn)/max(2d-10,1d0-taupmx))
505  wtmat(ibin,3)=wtmat(ibin,3)+(ataup1/ataup3)*
506  & taup/max(2d-10,1d0-taup)
507  ENDIF
508 
509 C...Sum up y* cross-section pieces in points used.
510  ELSEIF(ivar.EQ.3) THEN
511  yst=vintpt(iacc,12)
512  ystmin=vintpt(iacc,2)
513  ystmax=vintpt(iacc,22)
514  ayst0=ystmax-ystmin
515  ayst1=0.5d0*(ystmax-ystmin)**2
516  ayst2=ayst1
517  ayst3=2d0*(atan(exp(ystmax))-atan(exp(ystmin)))
518  wtmat(ibin,1)=wtmat(ibin,1)+(ayst0/ayst1)*(yst-ystmin)
519  wtmat(ibin,2)=wtmat(ibin,2)+(ayst0/ayst2)*(ystmax-yst)
520  wtmat(ibin,3)=wtmat(ibin,3)+(ayst0/ayst3)/cosh(yst)
521  IF(mint(45).EQ.3) THEN
522  taue=vintpt(iacc,11)
523  IF(istsb.GE.3.AND.istsb.LE.5) taue=vintpt(iacc,16)
524  yst0=-0.5d0*log(taue)
525  ayst4=log(max(1d-10,exp(yst0-ystmin)-1d0)/
526  & max(1d-10,exp(yst0-ystmax)-1d0))
527  wtmat(ibin,4)=wtmat(ibin,4)+(ayst0/ayst4)/
528  & max(1d-10,1d0-exp(yst-yst0))
529  ENDIF
530  IF(mint(46).EQ.3) THEN
531  taue=vintpt(iacc,11)
532  IF(istsb.GE.3.AND.istsb.LE.5) taue=vintpt(iacc,16)
533  yst0=-0.5d0*log(taue)
534  ayst5=log(max(1d-10,exp(yst0+ystmax)-1d0)/
535  & max(1d-10,exp(yst0+ystmin)-1d0))
536  wtmat(ibin,nbin)=wtmat(ibin,nbin)+(ayst0/ayst5)/
537  & max(1d-10,1d0-exp(-yst-yst0))
538  ENDIF
539 
540 C...Sum up cos(theta-hat) cross-section pieces in points used.
541  ELSE
542  rm34=max(1d-20,2d0*sqm3*sqm4/(vintpt(iacc,11)*vint(2))**2)
543  rsqm=1d0+rm34
544  cthmax=sqrt(1d0-4d0*vint(71)**2/(taumax*vint(2)))
545  cthmin=-cthmax
546  IF(cthmax.GT.0.9999d0) rm34=max(rm34,2d0*vint(71)**2/
547  & (taumax*vint(2)))
548  acth1=cthmax-cthmin
549  acth2=log(max(rm34,rsqm-cthmin)/max(rm34,rsqm-cthmax))
550  acth3=log(max(rm34,rsqm+cthmax)/max(rm34,rsqm+cthmin))
551  acth4=1d0/max(rm34,rsqm-cthmax)-1d0/max(rm34,rsqm-cthmin)
552  acth5=1d0/max(rm34,rsqm+cthmin)-1d0/max(rm34,rsqm+cthmax)
553  cth=vintpt(iacc,13)
554  wtmat(ibin,1)=wtmat(ibin,1)+1d0
555  wtmat(ibin,2)=wtmat(ibin,2)+(acth1/acth2)/
556  & max(rm34,rsqm-cth)
557  wtmat(ibin,3)=wtmat(ibin,3)+(acth1/acth3)/
558  & max(rm34,rsqm+cth)
559  wtmat(ibin,4)=wtmat(ibin,4)+(acth1/acth4)/
560  & max(rm34,rsqm-cth)**2
561  wtmat(ibin,5)=wtmat(ibin,5)+(acth1/acth5)/
562  & max(rm34,rsqm+cth)**2
563  ENDIF
564  180 CONTINUE
565 
566 C...Check that equation system solvable.
567  IF(mstp(122).GE.2) WRITE(mstu(11),5400) cvar(ivar)
568  msolv=1
569  wtrels=0d0
570  DO 190 ibin=1,nbin
571  IF(mstp(122).GE.2) WRITE(mstu(11),5500) (wtmat(ibin,ired),
572  & ired=1,nbin),wtrel(ibin)
573  IF(narel(ibin).EQ.0) msolv=0
574  wtrels=wtrels+wtrel(ibin)
575  190 CONTINUE
576  IF(abs(wtrels).LT.1d-20) msolv=0
577 
578 C...Solve to find relative importance of cross-section pieces.
579  IF(msolv.EQ.1) THEN
580  DO 200 ibin=1,nbin
581  wtreln(ibin)=max(0.1d0,wtrel(ibin)/wtrels)
582  200 CONTINUE
583  DO 230 ired=1,nbin-1
584  DO 220 ibin=ired+1,nbin
585  IF(abs(wtmat(ired,ired)).LT.1d-20) THEN
586  msolv=0
587  goto 260
588  ENDIF
589  rqt=wtmat(ibin,ired)/wtmat(ired,ired)
590  wtrel(ibin)=wtrel(ibin)-rqt*wtrel(ired)
591  DO 210 icoe=ired,nbin
592  wtmat(ibin,icoe)=wtmat(ibin,icoe)-rqt*wtmat(ired,icoe)
593  210 CONTINUE
594  220 CONTINUE
595  230 CONTINUE
596  DO 250 ired=nbin,1,-1
597  DO 240 icoe=ired+1,nbin
598  wtrel(ired)=wtrel(ired)-wtmat(ired,icoe)*coefu(icoe)
599  240 CONTINUE
600  coefu(ired)=wtrel(ired)/wtmat(ired,ired)
601  250 CONTINUE
602  ENDIF
603 
604 C...Share evenly if failure.
605  260 IF(msolv.EQ.0) THEN
606  DO 270 ibin=1,nbin
607  coefu(ibin)=1d0
608  wtreln(ibin)=0.1d0
609  IF(wtrels.GT.0d0) wtreln(ibin)=max(0.1d0,
610  & wtrel(ibin)/wtrels)
611  270 CONTINUE
612  ENDIF
613 
614 C...Normalize coefficients, with piece shared democratically.
615  coefsu=0d0
616  wtrels=0d0
617  DO 280 ibin=1,nbin
618  coefu(ibin)=max(0d0,coefu(ibin))
619  coefsu=coefsu+coefu(ibin)
620  wtrels=wtrels+wtreln(ibin)
621  280 CONTINUE
622  IF(coefsu.GT.0d0) THEN
623  DO 290 ibin=1,nbin
624  coefo(ibin)=parp(122)/nbin+(1d0-parp(122))*0.5d0*
625  & (coefu(ibin)/coefsu+wtreln(ibin)/wtrels)
626  290 CONTINUE
627  ELSE
628  DO 300 ibin=1,nbin
629  coefo(ibin)=1d0/nbin
630  300 CONTINUE
631  ENDIF
632  IF(ivar.EQ.1) ioff=0
633  IF(ivar.EQ.2) ioff=17
634  IF(ivar.EQ.3) ioff=7
635  IF(ivar.EQ.4) ioff=12
636  DO 310 ibin=1,nbin
637  icof=ioff+ibin
638  IF(ivar.EQ.1.AND.ibin.GT.2+2*mint(72)) icof=7
639  IF(ivar.EQ.3.AND.ibin.EQ.4.AND.mint(45).NE.3) icof=icof+1
640  coef(isub,icof)=coefo(ibin)
641  310 CONTINUE
642  IF(mstp(122).GE.2) WRITE(mstu(11),5600) cvar(ivar),
643  & (coefo(ibin),ibin=1,nbin)
644  320 CONTINUE
645 
646 C...Find two most promising maxima among points previously determined.
647  DO 330 j=1,4
648  iaccmx(j)=0
649  sigsmx(j)=0d0
650  330 CONTINUE
651  nmax=0
652  DO 390 iacc=1,nacc
653  DO 340 j=1,30
654  vint(10+j)=vintpt(iacc,j)
655  340 CONTINUE
656  IF(istsb.NE.5) THEN
657  CALL pysigh(nchn,sigs)
658  IF(mwtxs.EQ.1) THEN
659  CALL pyevwt(wtxs)
660  sigs=wtxs*sigs
661  ENDIF
662  ELSE
663  sigs=0d0
664  DO 350 ikin3=1,mstp(129)
665  CALL pykmap(5,0,0d0)
666  IF(mint(51).EQ.1) goto 350
667  CALL pysigh(nchn,sigtmp)
668  IF(mwtxs.EQ.1) THEN
669  CALL pyevwt(wtxs)
670  sigtmp=wtxs*sigtmp
671  ENDIF
672  IF(sigtmp.GT.sigs) sigs=sigtmp
673  350 CONTINUE
674  ENDIF
675  ieq=0
676  DO 360 imv=1,nmax
677  IF(abs(sigs-sigsmx(imv)).LT.1d-4*(sigs+sigsmx(imv))) ieq=imv
678  360 CONTINUE
679  IF(ieq.EQ.0) THEN
680  DO 370 imv=nmax,1,-1
681  iin=imv+1
682  IF(sigs.LE.sigsmx(imv)) goto 380
683  iaccmx(imv+1)=iaccmx(imv)
684  sigsmx(imv+1)=sigsmx(imv)
685  370 CONTINUE
686  iin=1
687  380 iaccmx(iin)=iacc
688  sigsmx(iin)=sigs
689  IF(nmax.LE.1) nmax=nmax+1
690  ENDIF
691  390 CONTINUE
692 
693 C...Read out starting position for search.
694  IF(mstp(122).GE.2) WRITE(mstu(11),5700)
695  sigsam=sigsmx(1)
696  DO 440 imax=1,nmax
697  iacc=iaccmx(imax)
698  mtau=mvarpt(iacc,1)
699  mtaup=mvarpt(iacc,2)
700  myst=mvarpt(iacc,3)
701  mcth=mvarpt(iacc,4)
702  vtau=0.5d0
703  vyst=0.5d0
704  vcth=0.5d0
705  vtaup=0.5d0
706 
707 C...Starting point and step size in parameter space.
708  DO 430 irpt=1,2
709  DO 420 ivar=1,4
710  IF(npts(ivar).EQ.1) goto 420
711  IF(ivar.EQ.1) vvar=vtau
712  IF(ivar.EQ.2) vvar=vtaup
713  IF(ivar.EQ.3) vvar=vyst
714  IF(ivar.EQ.4) vvar=vcth
715  IF(ivar.EQ.1) mvar=mtau
716  IF(ivar.EQ.2) mvar=mtaup
717  IF(ivar.EQ.3) mvar=myst
718  IF(ivar.EQ.4) mvar=mcth
719  IF(irpt.EQ.1) vdel=0.1d0
720  IF(irpt.EQ.2) vdel=max(0.01d0,min(0.05d0,vvar-0.02d0,
721  & 0.98d0-vvar))
722  IF(irpt.EQ.1) vmar=0.02d0
723  IF(irpt.EQ.2) vmar=0.002d0
724  imov0=1
725  IF(irpt.EQ.1.AND.ivar.EQ.1) imov0=0
726  DO 410 imov=imov0,8
727 
728 C...Define new point in parameter space.
729  IF(imov.EQ.0) THEN
730  inew=2
731  vnew=vvar
732  ELSEIF(imov.EQ.1) THEN
733  inew=3
734  vnew=vvar+vdel
735  ELSEIF(imov.EQ.2) THEN
736  inew=1
737  vnew=vvar-vdel
738  ELSEIF(sigssm(3).GE.max(sigssm(1),sigssm(2)).AND.
739  & vvar+2d0*vdel.LT.1d0-vmar) THEN
740  vvar=vvar+vdel
741  sigssm(1)=sigssm(2)
742  sigssm(2)=sigssm(3)
743  inew=3
744  vnew=vvar+vdel
745  ELSEIF(sigssm(1).GE.max(sigssm(2),sigssm(3)).AND.
746  & vvar-2d0*vdel.GT.vmar) THEN
747  vvar=vvar-vdel
748  sigssm(3)=sigssm(2)
749  sigssm(2)=sigssm(1)
750  inew=1
751  vnew=vvar-vdel
752  ELSEIF(sigssm(3).GE.sigssm(1)) THEN
753  vdel=0.5d0*vdel
754  vvar=vvar+vdel
755  sigssm(1)=sigssm(2)
756  inew=2
757  vnew=vvar
758  ELSE
759  vdel=0.5d0*vdel
760  vvar=vvar-vdel
761  sigssm(3)=sigssm(2)
762  inew=2
763  vnew=vvar
764  ENDIF
765 
766 C...Convert to relevant variables and find derived new limits.
767  ilerr=0
768  IF(ivar.EQ.1) THEN
769  vtau=vnew
770  CALL pykmap(1,mtau,vtau)
771  IF(istsb.GE.3.AND.istsb.LE.5) THEN
772  CALL pyklim(4)
773  IF(mint(51).EQ.1) ilerr=1
774  ENDIF
775  ENDIF
776  IF(ivar.LE.2.AND.istsb.GE.3.AND.istsb.LE.5.AND.
777  & ilerr.EQ.0) THEN
778  IF(ivar.EQ.2) vtaup=vnew
779  CALL pykmap(4,mtaup,vtaup)
780  ENDIF
781  IF(ivar.LE.2.AND.ilerr.EQ.0) THEN
782  CALL pyklim(2)
783  IF(mint(51).EQ.1) ilerr=1
784  ENDIF
785  IF(ivar.LE.3.AND.ilerr.EQ.0) THEN
786  IF(ivar.EQ.3) vyst=vnew
787  CALL pykmap(2,myst,vyst)
788  CALL pyklim(3)
789  IF(mint(51).EQ.1) ilerr=1
790  ENDIF
791  IF((istsb.EQ.2.OR.istsb.EQ.4.OR.istsb.EQ.6).AND.
792  & ilerr.EQ.0) THEN
793  IF(ivar.EQ.4) vcth=vnew
794  CALL pykmap(3,mcth,vcth)
795  ENDIF
796  IF(isub.EQ.96) vint(25)=vint(21)*(1.-vint(23)**2)
797 
798 C...Evaluate cross-section. Save new maximum. Final maximum.
799  IF(ilerr.NE.0) THEN
800  sigs=0.
801  ELSEIF(istsb.NE.5) THEN
802  CALL pysigh(nchn,sigs)
803  IF(mwtxs.EQ.1) THEN
804  CALL pyevwt(wtxs)
805  sigs=wtxs*sigs
806  ENDIF
807  ELSE
808  sigs=0d0
809  DO 400 ikin3=1,mstp(129)
810  CALL pykmap(5,0,0d0)
811  IF(mint(51).EQ.1) goto 400
812  CALL pysigh(nchn,sigtmp)
813  IF(mwtxs.EQ.1) THEN
814  CALL pyevwt(wtxs)
815  sigtmp=wtxs*sigtmp
816  ENDIF
817  IF(sigtmp.GT.sigs) sigs=sigtmp
818  400 CONTINUE
819  ENDIF
820  sigssm(inew)=sigs
821  IF(sigs.GT.sigsam) sigsam=sigs
822  IF(mstp(122).GE.2) WRITE(mstu(11),5800) imax,ivar,mvar,
823  & imov,vnew,vint(21),vint(22),vint(23),vint(26),sigs
824  410 CONTINUE
825  420 CONTINUE
826  430 CONTINUE
827  440 CONTINUE
828  IF(mstp(121).EQ.1) sigsam=parp(121)*sigsam
829  xsec(isub,1)=1.05d0*sigsam
830  IF(mint(141).NE.0.OR.mint(142).NE.0) xsec(isub,1)=
831  & wtgaga*xsec(isub,1)
832  450 CONTINUE
833  IF(mstp(173).EQ.1.AND.isub.NE.96) xsec(isub,1)=
834  & parp(174)*xsec(isub,1)
835  IF(isub.NE.96) xsec(0,1)=xsec(0,1)+xsec(isub,1)
836  460 CONTINUE
837  mint(51)=0
838 
839 C...Print summary table.
840  IF(mint(121).EQ.1.AND.nposi.EQ.0) THEN
841  IF(mstp(127).NE.1) THEN
842  WRITE(mstu(11),5900)
843  CALL pystop(1)
844  ELSE
845  WRITE(mstu(11),6400)
846  msti(53)=1
847  ENDIF
848  ENDIF
849  IF(mstp(122).GE.1) THEN
850  WRITE(mstu(11),6000)
851  WRITE(mstu(11),6100)
852  DO 470 isub=1,500
853  IF(msub(isub).NE.1.AND.isub.NE.96) goto 470
854  IF(isub.EQ.96.AND.mint(50).EQ.0) goto 470
855  IF(isub.EQ.96.AND.msub(95).NE.1.AND.mod(mstp(81),10).LE.0)
856  & goto 470
857  IF(isub.EQ.96.AND.mint(49).EQ.0.AND.mstp(131).EQ.0) goto 470
858  IF(msub(95).EQ.1.AND.(isub.EQ.11.OR.isub.EQ.12.OR.isub.EQ.13
859  & .OR.isub.EQ.28.OR.isub.EQ.53.OR.isub.EQ.68)) goto 470
860  IF(msub(95).EQ.1.AND.isub.GE.381.AND.isub.LE.386) goto 470
861  WRITE(mstu(11),6200) isub,proc(isub),xsec(isub,1)
862  470 CONTINUE
863  WRITE(mstu(11),6300)
864  ENDIF
865 
866 C...Format statements for maximization results.
867  5000 FORMAT(/1x,'Coefficient optimization and maximum search for ',
868  &'subprocess no',i4/1x,'Coefficient modes tau',10x,'y*',9x,
869  &'cth',9x,'tau''',7x,'sigma')
870  5100 FORMAT(1x,'Warning: requested subprocess ',i3,' has no allowed ',
871  &'phase space.'/1x,'Process switched off!')
872  5200 FORMAT(1x,4i4,f12.8,f12.6,f12.7,f12.8,1p,d12.4)
873  5300 FORMAT(1x,'Warning: requested subprocess ',i3,' has vanishing ',
874  &'cross-section.'/1x,'Process switched off!')
875  5400 FORMAT(1x,'Coefficients of equation system to be solved for ',a4)
876  5500 FORMAT(1x,1p,8d11.3)
877  5600 FORMAT(1x,'Result for ',a4,':',7f9.4)
878  5700 FORMAT(1x,'Maximum search for given coefficients'/2x,'MAX VAR ',
879  &'MOD MOV VNEW',7x,'tau',7x,'y*',8x,'cth',7x,'tau''',7x,'sigma')
880  5800 FORMAT(1x,4i4,f8.4,f11.7,f9.3,f11.6,f11.7,1p,d12.4)
881  5900 FORMAT(1x,'Error: no requested process has non-vanishing ',
882  &'cross-section.'/1x,'Execution stopped!')
883  6000 FORMAT(/1x,8('*'),1x,'PYMAXI: summary of differential ',
884  &'cross-section maximum search',1x,8('*'))
885  6100 FORMAT(/11x,58('=')/11x,'I',38x,'I',17x,'I'/11x,'I ISUB ',
886  &'Subprocess name',15x,'I Maximum value I'/11x,'I',38x,'I',
887  &17x,'I'/11x,58('=')/11x,'I',38x,'I',17x,'I')
888  6200 FORMAT(11x,'I',2x,i3,3x,a28,2x,'I',2x,1p,d12.4,3x,'I')
889  6300 FORMAT(11x,'I',38x,'I',17x,'I'/11x,58('='))
890  6400 FORMAT(1x,'Error: no requested process has non-vanishing ',
891  &'cross-section.'/
892  &1x,'Execution will stop if you try to generate events.')
893 
894  RETURN
895  END