Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyhimaxi.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhimaxi.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE pyhimaxi
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  SAVE /ludat1/
11  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
12  SAVE /ludat2/
13  common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
14  SAVE /pyhisubs/
15  common/pyhipars/mstp(200),parp(200),msti(200),pari(200)
16  SAVE /pyhipars/
17  common/pyhiint1/mint(400),vint(400)
18  SAVE /pyhiint1/
19  common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
20  SAVE /pyhiint2/
21  common/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
22  SAVE /pyhiint3/
23  common/pyhiint4/widp(21:40,0:40),wide(21:40,0:40),wids(21:40,3)
24  SAVE /pyhiint4/
25  common/pyhiint5/ngen(0:200,3),xsec(0:200,3)
26  SAVE /pyhiint5/
27  common/pyhiint6/proc(0:200)
28  CHARACTER proc*28
29  SAVE /pyhiint6/
30  CHARACTER cvar(4)*4
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),
33  &sigssm(3)
34  DATA cvar/'tau ','tau''','y* ','cth '/
35 
36 C...Select subprocess to study: skip cases not applicable.
37  vint(143)=1.
38  vint(144)=1.
39  xsec(0,1)=0.
40  DO 350 isub=1,200
41  IF(isub.GE.91.AND.isub.LE.95) THEN
42  xsec(isub,1)=vint(isub+11)
43  IF(msub(isub).NE.1) goto 350
44  goto 340
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
51  ELSE
52  IF(msub(isub).NE.1) goto 350
53  ENDIF
54  mint(1)=isub
55  istsb=iset(isub)
56  IF(isub.EQ.96) istsb=2
57  IF(mstp(122).GE.2) WRITE(mstu(11),1000) isub
58 
59 C...Find resonances (explicit or implicit in cross-section).
60  mint(72)=0
61  kfr1=0
62  IF(istsb.EQ.1.OR.istsb.EQ.3) THEN
63  kfr1=kfpr(isub,1)
64  ELSEIF(isub.GE.71.AND.isub.LE.77) THEN
65  kfr1=25
66  ENDIF
67  IF(kfr1.NE.0) THEN
68  taur1=pmas(kfr1,1)**2/vint(2)
69  gamr1=pmas(kfr1,1)*pmas(kfr1,2)/vint(2)
70  mint(72)=1
71  mint(73)=kfr1
72  vint(73)=taur1
73  vint(74)=gamr1
74  ENDIF
75  IF(isub.EQ.141) THEN
76  kfr2=23
77  taur2=pmas(kfr2,1)**2/vint(2)
78  gamr2=pmas(kfr2,1)*pmas(kfr2,2)/vint(2)
79  mint(72)=2
80  mint(74)=kfr2
81  vint(75)=taur2
82  vint(76)=gamr2
83  ENDIF
84 
85 C...Find product masses and minimum pT of process.
86  sqm3=0.
87  sqm4=0.
88  mint(71)=0
89  vint(71)=ckin(3)
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
94  IF(mint(71).EQ.1) vint(71)=max(ckin(3),ckin(5))
95  IF(isub.EQ.96.AND.mstp(82).LE.1) vint(71)=parp(81)
96  IF(isub.EQ.96.AND.mstp(82).GE.2) vint(71)=0.08*parp(82)
97  ENDIF
98  vint(63)=sqm3
99  vint(64)=sqm4
100 
101 C...Number of points for each variable: tau, tau', y*, cos(theta-hat).
102  npts(1)=2+2*mint(72)
103  IF(mint(43).EQ.1.AND.(istsb.EQ.1.OR.istsb.EQ.2)) npts(1)=1
104  npts(2)=1
105  IF(mint(43).GE.2.AND.(istsb.EQ.3.OR.istsb.EQ.4)) npts(2)=2
106  npts(3)=1
107  IF(mint(43).EQ.4) npts(3)=3
108  npts(4)=1
109  IF(istsb.EQ.2.OR.istsb.EQ.4) npts(4)=5
110  ntry=npts(1)*npts(2)*npts(3)*npts(4)
111 
112 C...Reset coefficients of cross-section weighting.
113  DO 100 j=1,20
114  100 coef(isub,j)=0.
115  coef(isub,1)=1.
116  coef(isub,7)=0.5
117  coef(isub,8)=0.5
118  coef(isub,10)=1.
119  coef(isub,15)=1.
120  mcth=0
121  mtaup=0
122  cth=0.
123  taup=0.
124  sigsam=0.
125 
126 C...Find limits and select tau, y*, cos(theta-hat) and tau' values,
127 C...in grid of phase space points.
128  CALL pyhiklim(1)
129  nacc=0
130  DO 120 itry=1,ntry
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))
133  CALL pyhikmap(1,mtau,0.5)
134  IF(istsb.EQ.3.OR.istsb.EQ.4) CALL pyhiklim(4)
135  ENDIF
136  IF((istsb.EQ.3.OR.istsb.EQ.4).AND.mod(itry-1,npts(3)*npts(4)).
137  &eq.0) THEN
138  mtaup=1+mod((itry-1)/(npts(3)*npts(4)),npts(2))
139  CALL pyhikmap(4,mtaup,0.5)
140  ENDIF
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))
144  CALL pyhikmap(2,myst,0.5)
145  CALL pyhiklim(3)
146  ENDIF
147  IF(istsb.EQ.2.OR.istsb.EQ.4) THEN
148  mcth=1+mod(itry-1,npts(4))
149  CALL pyhikmap(3,mcth,0.5)
150  ENDIF
151  IF(isub.EQ.96) vint(25)=vint(21)*(1.-vint(23)**2)
152 
153 C...Calculate and store cross-section.
154  mint(51)=0
155  CALL pyhiklim(0)
156  IF(mint(51).EQ.1) goto 120
157  nacc=nacc+1
158  mvarpt(nacc,1)=mtau
159  mvarpt(nacc,2)=mtaup
160  mvarpt(nacc,3)=myst
161  mvarpt(nacc,4)=mcth
162  DO 110 j=1,30
163  110 vintpt(nacc,j)=vint(10+j)
164  CALL pyhisigh(nchn,sigs)
165  sigspt(nacc)=sigs
166  IF(sigs.GT.sigsam) sigsam=sigs
167  IF(mstp(122).GE.2) WRITE(mstu(11),1100) mtau,mtaup,myst,mcth,
168  &vint(21),vint(22),vint(23),vint(26),sigs
169  120 CONTINUE
170  IF(sigsam.EQ.0.) THEN
171  WRITE(mstu(11),1200) isub
172  stop
173  ENDIF
174 
175 C...Calculate integrals in tau and y* over maximal phase space limits.
176  taumin=vint(11)
177  taumax=vint(31)
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))/
183  & gamr1
184  ENDIF
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))/
188  & gamr2
189  ENDIF
190  ystmin=0.5*log(taumin)
191  ystmax=-ystmin
192  ayst0=ystmax-ystmin
193  ayst1=0.5*(ystmax-ystmin)**2
194  ayst3=2.*(atan(exp(ystmax))-atan(exp(ystmin)))
195 
196 C...Reset. Sum up cross-sections in points calculated.
197  DO 230 ivar=1,4
198  IF(npts(ivar).EQ.1) goto 230
199  IF(isub.EQ.96.AND.ivar.EQ.4) goto 230
200  nbin=npts(ivar)
201  DO 130 j1=1,nbin
202  narel(j1)=0
203  wtrel(j1)=0.
204  coefu(j1)=0.
205  DO 130 j2=1,nbin
206  130 wtmat(j1,j2)=0.
207  DO 140 iacc=1,nacc
208  ibin=mvarpt(iacc,ivar)
209  narel(ibin)=narel(ibin)+1
210  wtrel(ibin)=wtrel(ibin)+sigspt(iacc)
211 
212 C...Sum up tau cross-section pieces in points used.
213  IF(ivar.EQ.1) THEN
214  tau=vintpt(iacc,11)
215  wtmat(ibin,1)=wtmat(ibin,1)+1.
216  wtmat(ibin,2)=wtmat(ibin,2)+(atau1/atau2)/tau
217  IF(nbin.GE.3) THEN
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)
221  ENDIF
222  IF(nbin.GE.5) THEN
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)
226  ENDIF
227 
228 C...Sum up tau' cross-section pieces in points used.
229  ELSEIF(ivar.EQ.2) THEN
230  tau=vintpt(iacc,11)
231  taup=vintpt(iacc,16)
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/
238  & taup
239 
240 C...Sum up y* and cos(theta-hat) cross-section pieces in points used.
241  ELSEIF(ivar.EQ.3) THEN
242  yst=vintpt(iacc,12)
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)
246  ELSE
247  rm34=2.*sqm3*sqm4/(vintpt(iacc,11)*vint(2))**2
248  rsqm=1.+rm34
249  cthmax=sqrt(1.-4.*vint(71)**2/(taumax*vint(2)))
250  cthmin=-cthmax
251  IF(cthmax.GT.0.9999) rm34=max(rm34,2.*vint(71)**2/
252  & (taumax*vint(2)))
253  acth1=cthmax-cthmin
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)
258  cth=vintpt(iacc,13)
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
264  ENDIF
265  140 CONTINUE
266 
267 C...Check that equation system solvable; else trivial way out.
268  IF(mstp(122).GE.2) WRITE(mstu(11),1300) cvar(ivar)
269  msolv=1
270  DO 150 ibin=1,nbin
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
274  IF(msolv.EQ.0) THEN
275  DO 160 ibin=1,nbin
276  160 coefu(ibin)=1.
277 
278 C...Solve to find relative importance of cross-section pieces.
279  ELSE
280  DO 170 ired=1,nbin-1
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)
290  ENDIF
291 
292 C...Normalize coefficients, with piece shared democratically.
293  coefsu=0.
294  DO 200 ibin=1,nbin
295  coefu(ibin)=max(0.,coefu(ibin))
296  200 coefsu=coefsu+coefu(ibin)
297  IF(ivar.EQ.1) ioff=0
298  IF(ivar.EQ.2) ioff=14
299  IF(ivar.EQ.3) ioff=6
300  IF(ivar.EQ.4) ioff=9
301  IF(coefsu.GT.0.) THEN
302  DO 210 ibin=1,nbin
303  210 coef(isub,ioff+ibin)=parp(121)/nbin+(1.-parp(121))*coefu(ibin)/
304  & coefsu
305  ELSE
306  DO 220 ibin=1,nbin
307  220 coef(isub,ioff+ibin)=1./nbin
308  ENDIF
309  IF(mstp(122).GE.2) WRITE(mstu(11),1500) cvar(ivar),
310  &(coef(isub,ioff+ibin),ibin=1,nbin)
311  230 CONTINUE
312 
313 C...Find two most promising maxima among points previously determined.
314  DO 240 j=1,4
315  iaccmx(j)=0
316  240 sigsmx(j)=0.
317  nmax=0
318  DO 290 iacc=1,nacc
319  DO 250 j=1,30
320  250 vint(10+j)=vintpt(iacc,j)
321  CALL pyhisigh(nchn,sigs)
322  ieq=0
323  DO 260 imv=1,nmax
324  260 IF(abs(sigs-sigsmx(imv)).LT.1e-4*(sigs+sigsmx(imv))) ieq=imv
325  IF(ieq.EQ.0) THEN
326  DO 270 imv=nmax,1,-1
327  iin=imv+1
328  IF(sigs.LE.sigsmx(imv)) goto 280
329  iaccmx(imv+1)=iaccmx(imv)
330  270 sigsmx(imv+1)=sigsmx(imv)
331  iin=1
332  280 iaccmx(iin)=iacc
333  sigsmx(iin)=sigs
334  IF(nmax.LE.1) nmax=nmax+1
335  ENDIF
336  290 CONTINUE
337 
338 C...Read out starting position for search.
339  IF(mstp(122).GE.2) WRITE(mstu(11),1600)
340  sigsam=sigsmx(1)
341  DO 330 imax=1,nmax
342  iacc=iaccmx(imax)
343  mtau=mvarpt(iacc,1)
344  mtaup=mvarpt(iacc,2)
345  myst=mvarpt(iacc,3)
346  mcth=mvarpt(iacc,4)
347  vtau=0.5
348  vyst=0.5
349  vcth=0.5
350  vtaup=0.5
351 
352 C...Starting point and step size in parameter space.
353  DO 320 irpt=1,2
354  DO 310 ivar=1,4
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
368  imov0=1
369  IF(irpt.EQ.1.AND.ivar.EQ.1) imov0=0
370  DO 300 imov=imov0,8
371 
372 C...Define new point in parameter space.
373  IF(imov.EQ.0) THEN
374  inew=2
375  vnew=vvar
376  ELSEIF(imov.EQ.1) THEN
377  inew=3
378  vnew=vvar+vdel
379  ELSEIF(imov.EQ.2) THEN
380  inew=1
381  vnew=vvar-vdel
382  ELSEIF(sigssm(3).GE.max(sigssm(1),sigssm(2)).AND.
383  &vvar+2.*vdel.LT.1.-vmar) THEN
384  vvar=vvar+vdel
385  sigssm(1)=sigssm(2)
386  sigssm(2)=sigssm(3)
387  inew=3
388  vnew=vvar+vdel
389  ELSEIF(sigssm(1).GE.max(sigssm(2),sigssm(3)).AND.
390  &vvar-2.*vdel.GT.vmar) THEN
391  vvar=vvar-vdel
392  sigssm(3)=sigssm(2)
393  sigssm(2)=sigssm(1)
394  inew=1
395  vnew=vvar-vdel
396  ELSEIF(sigssm(3).GE.sigssm(1)) THEN
397  vdel=0.5*vdel
398  vvar=vvar+vdel
399  sigssm(1)=sigssm(2)
400  inew=2
401  vnew=vvar
402  ELSE
403  vdel=0.5*vdel
404  vvar=vvar-vdel
405  sigssm(3)=sigssm(2)
406  inew=2
407  vnew=vvar
408  ENDIF
409 
410 C...Convert to relevant variables and find derived new limits.
411  IF(ivar.EQ.1) THEN
412  vtau=vnew
413  CALL pyhikmap(1,mtau,vtau)
414  IF(istsb.EQ.3.OR.istsb.EQ.4) CALL pyhiklim(4)
415  ENDIF
416  IF(ivar.LE.2.AND.(istsb.EQ.3.OR.istsb.EQ.4)) THEN
417  IF(ivar.EQ.2) vtaup=vnew
418  CALL pyhikmap(4,mtaup,vtaup)
419  ENDIF
420  IF(ivar.LE.2) CALL pyhiklim(2)
421  IF(ivar.LE.3) THEN
422  IF(ivar.EQ.3) vyst=vnew
423  CALL pyhikmap(2,myst,vyst)
424  CALL pyhiklim(3)
425  ENDIF
426  IF(istsb.EQ.2.OR.istsb.EQ.4) THEN
427  IF(ivar.EQ.4) vcth=vnew
428  CALL pyhikmap(3,mcth,vcth)
429  ENDIF
430  IF(isub.EQ.96) vint(25)=vint(21)*(1.-vint(23)**2)
431 
432 C...Evaluate cross-section. Save new maximum. Final maximum.
433  CALL pyhisigh(nchn,sigs)
434  sigssm(inew)=sigs
435  IF(sigs.GT.sigsam) sigsam=sigs
436  IF(mstp(122).GE.2) WRITE(mstu(11),1700) imax,ivar,mvar,imov,
437  &vnew,vint(21),vint(22),vint(23),vint(26),sigs
438  300 CONTINUE
439  310 CONTINUE
440  320 CONTINUE
441  IF(imax.EQ.1) sigs11=sigsam
442  330 CONTINUE
443  xsec(isub,1)=1.05*sigsam
444  340 IF(isub.NE.96) xsec(0,1)=xsec(0,1)+xsec(isub,1)
445  350 CONTINUE
446 
447 C...Print summary table.
448  IF(mstp(122).GE.1) THEN
449  WRITE(mstu(11),1800)
450  WRITE(mstu(11),1900)
451  DO 360 isub=1,200
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)
458  360 CONTINUE
459  WRITE(mstu(11),2100)
460  ENDIF
461 
462 C...Format statements for maximization results.
463  1000 FORMAT(/1x,'Coefficient optimization and maximum search for ',
464  &'subprocess no',i4/1x,'Coefficient modes tau',10x,'y*',9x,
465  &'cth',9x,'tau''',7x,'sigma')
466  1100 FORMAT(1x,4i4,f12.8,f12.6,f12.7,f12.8,1p,e12.4)
467  1200 FORMAT(1x,'Error: requested subprocess ',i3,' has vanishing ',
468  &'cross-section.'/1x,'Execution stopped!')
469  1300 FORMAT(1x,'Coefficients of equation system to be solved for ',a4)
470  1400 FORMAT(1x,1p,7e11.3)
471  1500 FORMAT(1x,'Result for ',a4,':',6f9.4)
472  1600 FORMAT(1x,'Maximum search for given coefficients'/2x,'MAX VAR ',
473  &'MOD MOV VNEW',7x,'tau',7x,'y*',8x,'cth',7x,'tau''',7x,'sigma')
474  1700 FORMAT(1x,4i4,f8.4,f11.7,f9.3,f11.6,f11.7,1p,e12.4)
475  1800 FORMAT(/1x,8('*'),1x,'PYHIMAXI: summary of differential ',
476  &'cross-section maximum search',1x,8('*'))
477  1900 FORMAT(/11x,58('=')/11x,'I',38x,'I',17x,'I'/11x,'I ISUB ',
478  &'Subprocess name',15x,'I Maximum value I'/11x,'I',38x,'I',
479  &17x,'I'/11x,58('=')/11x,'I',38x,'I',17x,'I')
480  2000 FORMAT(11x,'I',2x,i3,3x,a28,2x,'I',2x,1p,e12.4,3x,'I')
481  2100 FORMAT(11x,'I',38x,'I',17x,'I'/11x,58('='))
482 
483  RETURN
484  END