9 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
11 common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
13 common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
19 common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
21 common/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
23 common/pyhiint4/widp(21:40,0:40),wide(21:40,0:40),wids(21:40,3)
25 common/pyhiint5/
ngen(0:200,3),
xsec(0:200,3)
27 common/pyhiint6/proc(0:200)
31 dimension npts(4),mvarpt(200,4),vintpt(200,30),sigspt(200),
32 &narel(6),wtrel(6),wtmat(6,6),coefu(6),iaccmx(4),sigsmx(4),
34 DATA cvar/
'tau ',
'tau''',
'y* ',
'cth '/
41 IF(isub.GE.91.AND.isub.LE.95)
THEN
43 IF(msub(isub).NE.1) goto 350
45 ELSEIF(isub.EQ.96)
THEN
46 IF(
mint(43).NE.4) goto 350
47 IF(msub(95).NE.1.AND.
mstp(81).LE.0.AND.
mstp(131).LE.0) goto 350
48 ELSEIF(isub.EQ.11.OR.isub.EQ.12.OR.isub.EQ.13.OR.isub.EQ.28.OR.
49 &isub.EQ.53.OR.isub.EQ.68)
THEN
50 IF(msub(isub).NE.1.OR.msub(95).EQ.1) goto 350
52 IF(msub(isub).NE.1) goto 350
56 IF(isub.EQ.96) istsb=2
57 IF(
mstp(122).GE.2)
WRITE(mstu(11),1000) isub
62 IF(istsb.EQ.1.OR.istsb.EQ.3)
THEN
64 ELSEIF(isub.GE.71.AND.isub.LE.77)
THEN
68 taur1=pmas(kfr1,1)**2/
vint(2)
69 gamr1=pmas(kfr1,1)*pmas(kfr1,2)/
vint(2)
77 taur2=pmas(kfr2,1)**2/
vint(2)
78 gamr2=pmas(kfr2,1)*pmas(kfr2,2)/
vint(2)
90 IF(istsb.EQ.2.OR.istsb.EQ.4)
THEN
91 IF(kfpr(isub,1).NE.0) sqm3=pmas(kfpr(isub,1),1)**2
92 IF(kfpr(isub,2).NE.0) sqm4=pmas(kfpr(isub,2),1)**2
93 IF(min(sqm3,sqm4).LT.ckin(6)**2)
mint(71)=1
103 IF(
mint(43).EQ.1.AND.(istsb.EQ.1.OR.istsb.EQ.2)) npts(1)=1
105 IF(
mint(43).GE.2.AND.(istsb.EQ.3.OR.istsb.EQ.4)) npts(2)=2
107 IF(
mint(43).EQ.4) npts(3)=3
109 IF(istsb.EQ.2.OR.istsb.EQ.4) npts(4)=5
110 ntry=npts(1)*npts(2)*npts(3)*npts(4)
131 IF(mod(itry-1,npts(2)*npts(3)*npts(4)).EQ.0)
THEN
132 mtau=1+(itry-1)/(npts(2)*npts(3)*npts(4))
134 IF(istsb.EQ.3.OR.istsb.EQ.4) CALL
pyhiklim(4)
136 IF((istsb.EQ.3.OR.istsb.EQ.4).AND.mod(itry-1,npts(3)*npts(4)).
138 mtaup=1+mod((itry-1)/(npts(3)*npts(4)),npts(2))
141 IF(mod(itry-1,npts(3)*npts(4)).EQ.0) CALL
pyhiklim(2)
142 IF(mod(itry-1,npts(4)).EQ.0)
THEN
143 myst=1+mod((itry-1)/npts(4),npts(3))
147 IF(istsb.EQ.2.OR.istsb.EQ.4)
THEN
148 mcth=1+mod(itry-1,npts(4))
156 IF(
mint(51).EQ.1) goto 120
163 110 vintpt(nacc,
j)=
vint(10+
j)
166 IF(sigs.GT.sigsam) sigsam=sigs
167 IF(
mstp(122).GE.2)
WRITE(mstu(11),1100) mtau,mtaup,myst,mcth,
170 IF(sigsam.EQ.0.)
THEN
171 WRITE(mstu(11),1200) isub
178 atau1=
log(taumax/taumin)
179 atau2=(taumax-taumin)/(taumax*taumin)
180 IF(npts(1).GE.3)
THEN
181 atau3=
log(taumax/taumin*(taumin+taur1)/(taumax+taur1))/taur1
182 atau4=(atan((taumax-taur1)/gamr1)-atan((taumin-taur1)/gamr1))/
185 IF(npts(1).GE.5)
THEN
186 atau5=
log(taumax/taumin*(taumin+taur2)/(taumax+taur2))/taur2
187 atau6=(atan((taumax-taur2)/gamr2)-atan((taumin-taur2)/gamr2))/
190 ystmin=0.5*
log(taumin)
193 ayst1=0.5*(ystmax-ystmin)**2
194 ayst3=2.*(atan(exp(ystmax))-atan(exp(ystmin)))
198 IF(npts(ivar).EQ.1) goto 230
199 IF(isub.EQ.96.AND.ivar.EQ.4) goto 230
208 ibin=mvarpt(iacc,ivar)
209 narel(ibin)=narel(ibin)+1
210 wtrel(ibin)=wtrel(ibin)+sigspt(iacc)
215 wtmat(ibin,1)=wtmat(ibin,1)+1.
216 wtmat(ibin,2)=wtmat(ibin,2)+(atau1/atau2)/
tau
218 wtmat(ibin,3)=wtmat(ibin,3)+(atau1/atau3)/(
tau+taur1)
219 wtmat(ibin,4)=wtmat(ibin,4)+(atau1/atau4)*
tau/
220 & ((
tau-taur1)**2+gamr1**2)
223 wtmat(ibin,5)=wtmat(ibin,5)+(atau1/atau5)/(
tau+taur2)
224 wtmat(ibin,6)=wtmat(ibin,6)+(atau1/atau6)*
tau/
225 & ((
tau-taur2)**2+gamr2**2)
229 ELSEIF(ivar.EQ.2)
THEN
232 taupmn=vintpt(iacc,6)
233 taupmx=vintpt(iacc,26)
234 ataup1=
log(taupmx/taupmn)
235 ataup2=((1.-
tau/taupmx)**4-(1.-
tau/taupmn)**4)/(4.*
tau)
236 wtmat(ibin,1)=wtmat(ibin,1)+1.
237 wtmat(ibin,2)=wtmat(ibin,2)+(ataup1/ataup2)*(1.-
tau/taup)**3/
241 ELSEIF(ivar.EQ.3)
THEN
243 wtmat(ibin,1)=wtmat(ibin,1)+(ayst0/ayst1)*(yst-ystmin)
244 wtmat(ibin,2)=wtmat(ibin,2)+(ayst0/ayst1)*(ystmax-yst)
245 wtmat(ibin,3)=wtmat(ibin,3)+(ayst0/ayst3)/cosh(yst)
247 rm34=2.*sqm3*sqm4/(vintpt(iacc,11)*
vint(2))**2
249 cthmax=sqrt(1.-4.*
vint(71)**2/(taumax*
vint(2)))
251 IF(cthmax.GT.0.9999) rm34=
max(rm34,2.*
vint(71)**2/
254 acth2=
log(
max(rm34,rsqm-cthmin)/
max(rm34,rsqm-cthmax))
255 acth3=
log(
max(rm34,rsqm+cthmax)/
max(rm34,rsqm+cthmin))
256 acth4=1./
max(rm34,rsqm-cthmax)-1./
max(rm34,rsqm-cthmin)
257 acth5=1./
max(rm34,rsqm+cthmin)-1./
max(rm34,rsqm+cthmax)
259 wtmat(ibin,1)=wtmat(ibin,1)+1.
260 wtmat(ibin,2)=wtmat(ibin,2)+(acth1/acth2)/
max(rm34,rsqm-cth)
261 wtmat(ibin,3)=wtmat(ibin,3)+(acth1/acth3)/
max(rm34,rsqm+cth)
262 wtmat(ibin,4)=wtmat(ibin,4)+(acth1/acth4)/
max(rm34,rsqm-cth)**2
263 wtmat(ibin,5)=wtmat(ibin,5)+(acth1/acth5)/
max(rm34,rsqm+cth)**2
268 IF(
mstp(122).GE.2)
WRITE(mstu(11),1300) cvar(ivar)
271 IF(
mstp(122).GE.2)
WRITE(mstu(11),1400) (wtmat(ibin,ired),
272 &ired=1,
nbin),wtrel(ibin)
273 150
IF(narel(ibin).EQ.0) msolv=0
281 DO 170 ibin=ired+1,
nbin
282 rqt=wtmat(ibin,ired)/wtmat(ired,ired)
283 wtrel(ibin)=wtrel(ibin)-rqt*wtrel(ired)
284 DO 170 icoe=ired,
nbin
285 170 wtmat(ibin,icoe)=wtmat(ibin,icoe)-rqt*wtmat(ired,icoe)
286 DO 190 ired=
nbin,1,-1
287 DO 180 icoe=ired+1,
nbin
288 180 wtrel(ired)=wtrel(ired)-wtmat(ired,icoe)*coefu(icoe)
289 190 coefu(ired)=wtrel(ired)/wtmat(ired,ired)
295 coefu(ibin)=
max(0.,coefu(ibin))
296 200 coefsu=coefsu+coefu(ibin)
298 IF(ivar.EQ.2) ioff=14
301 IF(coefsu.GT.0.)
THEN
303 210 coef(isub,ioff+ibin)=
parp(121)/
nbin+(1.-
parp(121))*coefu(ibin)/
307 220 coef(isub,ioff+ibin)=1./
nbin
309 IF(
mstp(122).GE.2)
WRITE(mstu(11),1500) cvar(ivar),
310 &(coef(isub,ioff+ibin),ibin=1,
nbin)
320 250
vint(10+
j)=vintpt(iacc,
j)
324 260
IF(abs(sigs-sigsmx(imv)).LT.1
e-4*(sigs+sigsmx(imv))) ieq=imv
328 IF(sigs.LE.sigsmx(imv)) goto 280
329 iaccmx(imv+1)=iaccmx(imv)
330 270 sigsmx(imv+1)=sigsmx(imv)
334 IF(nmax.LE.1) nmax=nmax+1
339 IF(
mstp(122).GE.2)
WRITE(mstu(11),1600)
355 IF(npts(ivar).EQ.1) goto 310
356 IF(ivar.EQ.1) vvar=vtau
357 IF(ivar.EQ.2) vvar=vtaup
358 IF(ivar.EQ.3) vvar=vyst
359 IF(ivar.EQ.4) vvar=vcth
360 IF(ivar.EQ.1) mvar=mtau
361 IF(ivar.EQ.2) mvar=mtaup
362 IF(ivar.EQ.3) mvar=myst
363 IF(ivar.EQ.4) mvar=mcth
364 IF(irpt.EQ.1) vdel=0.1
365 IF(irpt.EQ.2) vdel=
max(0.01,min(0.05,vvar-0.02,0.98-vvar))
366 IF(irpt.EQ.1) vmar=0.02
367 IF(irpt.EQ.2) vmar=0.002
369 IF(irpt.EQ.1.AND.ivar.EQ.1) imov0=0
376 ELSEIF(imov.EQ.1)
THEN
379 ELSEIF(imov.EQ.2)
THEN
382 ELSEIF(sigssm(3).GE.
max(sigssm(1),sigssm(2)).AND.
383 &vvar+2.*vdel.LT.1.-vmar)
THEN
389 ELSEIF(sigssm(1).GE.
max(sigssm(2),sigssm(3)).AND.
390 &vvar-2.*vdel.GT.vmar)
THEN
396 ELSEIF(sigssm(3).GE.sigssm(1))
THEN
414 IF(istsb.EQ.3.OR.istsb.EQ.4) CALL
pyhiklim(4)
416 IF(ivar.LE.2.AND.(istsb.EQ.3.OR.istsb.EQ.4))
THEN
417 IF(ivar.EQ.2) vtaup=vnew
422 IF(ivar.EQ.3) vyst=vnew
426 IF(istsb.EQ.2.OR.istsb.EQ.4)
THEN
427 IF(ivar.EQ.4) vcth=vnew
435 IF(sigs.GT.sigsam) sigsam=sigs
436 IF(
mstp(122).GE.2)
WRITE(mstu(11),1700) imax,ivar,mvar,imov,
441 IF(imax.EQ.1) sigs11=sigsam
443 xsec(isub,1)=1.05*sigsam
448 IF(
mstp(122).GE.1)
THEN
452 IF(msub(isub).NE.1.AND.isub.NE.96) goto 360
453 IF(isub.EQ.96.AND.
mint(43).NE.4) goto 360
454 IF(isub.EQ.96.AND.msub(95).NE.1.AND.
mstp(81).LE.0) goto 360
455 IF(msub(95).EQ.1.AND.(isub.EQ.11.OR.isub.EQ.12.OR.isub.EQ.13.OR.
456 & isub.EQ.28.OR.isub.EQ.53.OR.isub.EQ.68)) goto 360
457 WRITE(mstu(11),2000) isub,proc(isub),
xsec(isub,1)
463 1000
FORMAT(/1
x,
'Coefficient optimization and maximum search for ',
464 &
'subprocess no',
i4/1
x,
'Coefficient modes tau',10
x,
'y*',9
x,
465 &
'cth',9
x,
'tau''',7
x,
'sigma')
466 1100
FORMAT(1
x,4
i4,f12.8,f12.6,f12.7,f12.8,1
p,e12.4)
467 1200
FORMAT(1
x,
'Error: requested subprocess ',
i3,
' has vanishing ',
468 &
'cross-section.'/1
x,
'Execution stopped!')
469 1300
FORMAT(1
x,
'Coefficients of equation system to be solved for ',a4)
470 1400
FORMAT(1
x,1
p,7e11.3)
471 1500
FORMAT(1
x,
'Result for ',a4,
':',6f9.4)
472 1600
FORMAT(1
x,
'Maximum search for given coefficients'/2
x,
'MAX VAR ',
473 &
'MOD MOV VNEW',7
x,
'tau',7
x,
'y*',8
x,
'cth',7
x,
'tau''',7
x,
'sigma')
474 1700
FORMAT(1
x,4
i4,f8.4,f11.7,f9.3,f11.6,f11.7,1
p,e12.4)
475 1800
FORMAT(/1
x,8(
'*'),1
x,
'PYHIMAXI: summary of differential ',
476 &
'cross-section maximum search',1
x,8(
'*'))
477 1900
FORMAT(/11
x,58(
'=')/11
x,
'I',38
x,
'I',17
x,
'I'/11
x,
'I ISUB ',
478 &
'Subprocess name',15
x,
'I Maximum value I'/11
x,
'I',38
x,
'I',
479 &17
x,
'I'/11
x,58(
'=')/11
x,
'I',38
x,
'I',17
x,
'I')
480 2000
FORMAT(11
x,
'I',2
x,
i3,3
x,a28,2
x,
'I',2
x,1
p,e12.4,3
x,
'I')
481 2100
FORMAT(11
x,
'I',38
x,
'I',17
x,
'I'/11
x,58(
'='))