pyhimaxi.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhimaxi.f

    
C*********************************************************************  
    
      SUBROUTINE pyhimaxi 
    
C...Finds optimal set of coefficients for kinematical variable selection    
C...and the maximum of the part of the differential cross-section used  
C...in the event weighting. 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /ludat2/ 
      common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200) 
      SAVE /pyhisubs/ 
      common/pyhipars/mstp(200),parp(200),msti(200),pari(200) 
      SAVE /pyhipars/ 
      common/pyhiint1/mint(400),vint(400) 
      SAVE /pyhiint1/ 
      common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2) 
      SAVE /pyhiint2/ 
      common/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)  
      SAVE /pyhiint3/ 
      common/pyhiint4/widp(21:40,0:40),wide(21:40,0:40),wids(21:40,3) 
      SAVE /pyhiint4/ 
      common/pyhiint5/ngen(0:200,3),xsec(0:200,3) 
      SAVE /pyhiint5/ 
      common/pyhiint6/proc(0:200) 
      CHARACTER proc*28 
      SAVE /pyhiint6/ 
      CHARACTER cvar(4)*4   
      dimension npts(4),mvarpt(200,4),vintpt(200,30),sigspt(200),   
     &narel(6),wtrel(6),wtmat(6,6),coefu(6),iaccmx(4),sigsmx(4),    
     &sigssm(3) 
      DATA cvar/'tau ','tau''','y*  ','cth '/   
    
C...Select subprocess to study: skip cases not applicable.  
      vint(143)=1.  
      vint(144)=1.  
      xsec(0,1)=0.  
      DO 350 isub=1,200 
      IF(isub.GE.91.AND.isub.LE.95) THEN    
        xsec(isub,1)=vint(isub+11)  
        IF(msub(isub).NE.1) goto 350    
        goto 340    
      ELSEIF(isub.EQ.96) THEN   
        IF(mint(43).NE.4) goto 350  
        IF(msub(95).NE.1.AND.mstp(81).LE.0.AND.mstp(131).LE.0) goto 350 
      ELSEIF(isub.EQ.11.OR.isub.EQ.12.OR.isub.EQ.13.OR.isub.EQ.28.OR.   
     &isub.EQ.53.OR.isub.EQ.68) THEN    
        IF(msub(isub).NE.1.OR.msub(95).EQ.1) goto 350   
      ELSE  
        IF(msub(isub).NE.1) goto 350    
      ENDIF 
      mint(1)=isub  
      istsb=iset(isub)  
      IF(isub.EQ.96) istsb=2    
      IF(mstp(122).GE.2) WRITE(mstu(11),1000) isub  
    
C...Find resonances (explicit or implicit in cross-section).    
      mint(72)=0    
      kfr1=0    
      IF(istsb.EQ.1.OR.istsb.EQ.3) THEN 
        kfr1=kfpr(isub,1)   
      ELSEIF(isub.GE.71.AND.isub.LE.77) THEN    
        kfr1=25 
      ENDIF 
      IF(kfr1.NE.0) THEN    
        taur1=pmas(kfr1,1)**2/vint(2)   
        gamr1=pmas(kfr1,1)*pmas(kfr1,2)/vint(2) 
        mint(72)=1  
        mint(73)=kfr1   
        vint(73)=taur1  
        vint(74)=gamr1  
      ENDIF 
      IF(isub.EQ.141) THEN  
        kfr2=23 
        taur2=pmas(kfr2,1)**2/vint(2)   
        gamr2=pmas(kfr2,1)*pmas(kfr2,2)/vint(2) 
        mint(72)=2  
        mint(74)=kfr2   
        vint(75)=taur2  
        vint(76)=gamr2  
      ENDIF 
    
C...Find product masses and minimum pT of process.  
      sqm3=0.   
      sqm4=0.   
      mint(71)=0    
      vint(71)=ckin(3)  
      IF(istsb.EQ.2.OR.istsb.EQ.4) THEN 
        IF(kfpr(isub,1).NE.0) sqm3=pmas(kfpr(isub,1),1)**2  
        IF(kfpr(isub,2).NE.0) sqm4=pmas(kfpr(isub,2),1)**2  
        IF(min(sqm3,sqm4).LT.ckin(6)**2) mint(71)=1 
        IF(mint(71).EQ.1) vint(71)=max(ckin(3),ckin(5)) 
        IF(isub.EQ.96.AND.mstp(82).LE.1) vint(71)=parp(81)  
        IF(isub.EQ.96.AND.mstp(82).GE.2) vint(71)=0.08*parp(82) 
      ENDIF 
      vint(63)=sqm3 
      vint(64)=sqm4 
    
C...Number of points for each variable: tau, tau', y*, cos(theta-hat).  
      npts(1)=2+2*mint(72)  
      IF(mint(43).EQ.1.AND.(istsb.EQ.1.OR.istsb.EQ.2)) npts(1)=1    
      npts(2)=1 
      IF(mint(43).GE.2.AND.(istsb.EQ.3.OR.istsb.EQ.4)) npts(2)=2    
      npts(3)=1 
      IF(mint(43).EQ.4) npts(3)=3   
      npts(4)=1 
      IF(istsb.EQ.2.OR.istsb.EQ.4) npts(4)=5    
      ntry=npts(1)*npts(2)*npts(3)*npts(4)  
    
C...Reset coefficients of cross-section weighting.  
      DO 100 j=1,20 
  100 coef(isub,j)=0.   
      coef(isub,1)=1.   
      coef(isub,7)=0.5  
      coef(isub,8)=0.5  
      coef(isub,10)=1.  
      coef(isub,15)=1.  
      mcth=0    
      mtaup=0   
      cth=0.    
      taup=0.   
      sigsam=0. 
    
C...Find limits and select tau, y*, cos(theta-hat) and tau' values, 
C...in grid of phase space points.  
      CALL pyhiklim(1)    
      nacc=0    
      DO 120 itry=1,ntry    
      IF(mod(itry-1,npts(2)*npts(3)*npts(4)).EQ.0) THEN 
        mtau=1+(itry-1)/(npts(2)*npts(3)*npts(4))   
        CALL pyhikmap(1,mtau,0.5) 
        IF(istsb.EQ.3.OR.istsb.EQ.4) CALL pyhiklim(4) 
      ENDIF 
      IF((istsb.EQ.3.OR.istsb.EQ.4).AND.mod(itry-1,npts(3)*npts(4)).    
     &eq.0) THEN    
        mtaup=1+mod((itry-1)/(npts(3)*npts(4)),npts(2)) 
        CALL pyhikmap(4,mtaup,0.5)    
      ENDIF 
      IF(mod(itry-1,npts(3)*npts(4)).EQ.0) CALL pyhiklim(2)   
      IF(mod(itry-1,npts(4)).EQ.0) THEN 
        myst=1+mod((itry-1)/npts(4),npts(3))    
        CALL pyhikmap(2,myst,0.5) 
        CALL pyhiklim(3)  
      ENDIF 
      IF(istsb.EQ.2.OR.istsb.EQ.4) THEN 
        mcth=1+mod(itry-1,npts(4))  
        CALL pyhikmap(3,mcth,0.5) 
      ENDIF 
      IF(isub.EQ.96) vint(25)=vint(21)*(1.-vint(23)**2) 
    
C...Calculate and store cross-section.  
      mint(51)=0    
      CALL pyhiklim(0)    
      IF(mint(51).EQ.1) goto 120    
      nacc=nacc+1   
      mvarpt(nacc,1)=mtau   
      mvarpt(nacc,2)=mtaup  
      mvarpt(nacc,3)=myst   
      mvarpt(nacc,4)=mcth   
      DO 110 j=1,30 
  110 vintpt(nacc,j)=vint(10+j) 
      CALL pyhisigh(nchn,sigs)    
      sigspt(nacc)=sigs 
      IF(sigs.GT.sigsam) sigsam=sigs    
      IF(mstp(122).GE.2) WRITE(mstu(11),1100) mtau,mtaup,myst,mcth, 
     &vint(21),vint(22),vint(23),vint(26),sigs  
  120 CONTINUE  
      IF(sigsam.EQ.0.) THEN 
        WRITE(mstu(11),1200) isub   
        stop    
      ENDIF 
    
C...Calculate integrals in tau and y* over maximal phase space limits.  
      taumin=vint(11)   
      taumax=vint(31)   
      atau1=log(taumax/taumin)  
      atau2=(taumax-taumin)/(taumax*taumin) 
      IF(npts(1).GE.3) THEN 
        atau3=log(taumax/taumin*(taumin+taur1)/(taumax+taur1))/taur1    
        atau4=(atan((taumax-taur1)/gamr1)-atan((taumin-taur1)/gamr1))/  
     &  gamr1   
      ENDIF 
      IF(npts(1).GE.5) THEN 
        atau5=log(taumax/taumin*(taumin+taur2)/(taumax+taur2))/taur2    
        atau6=(atan((taumax-taur2)/gamr2)-atan((taumin-taur2)/gamr2))/  
     &  gamr2   
      ENDIF 
      ystmin=0.5*log(taumin)    
      ystmax=-ystmin    
      ayst0=ystmax-ystmin   
      ayst1=0.5*(ystmax-ystmin)**2  
      ayst3=2.*(atan(exp(ystmax))-atan(exp(ystmin)))    
    
C...Reset. Sum up cross-sections in points calculated.  
      DO 230 ivar=1,4   
      IF(npts(ivar).EQ.1) goto 230  
      IF(isub.EQ.96.AND.ivar.EQ.4) goto 230 
      nbin=npts(ivar)   
      DO 130 j1=1,nbin  
      narel(j1)=0   
      wtrel(j1)=0.  
      coefu(j1)=0.  
      DO 130 j2=1,nbin  
  130 wtmat(j1,j2)=0.   
      DO 140 iacc=1,nacc    
      ibin=mvarpt(iacc,ivar)    
      narel(ibin)=narel(ibin)+1 
      wtrel(ibin)=wtrel(ibin)+sigspt(iacc)  
    
C...Sum up tau cross-section pieces in points used. 
      IF(ivar.EQ.1) THEN    
        tau=vintpt(iacc,11) 
        wtmat(ibin,1)=wtmat(ibin,1)+1.  
        wtmat(ibin,2)=wtmat(ibin,2)+(atau1/atau2)/tau   
        IF(nbin.GE.3) THEN  
          wtmat(ibin,3)=wtmat(ibin,3)+(atau1/atau3)/(tau+taur1) 
          wtmat(ibin,4)=wtmat(ibin,4)+(atau1/atau4)*tau/    
     &    ((tau-taur1)**2+gamr1**2) 
        ENDIF   
        IF(nbin.GE.5) THEN  
          wtmat(ibin,5)=wtmat(ibin,5)+(atau1/atau5)/(tau+taur2) 
          wtmat(ibin,6)=wtmat(ibin,6)+(atau1/atau6)*tau/    
     &    ((tau-taur2)**2+gamr2**2) 
        ENDIF   
    
C...Sum up tau' cross-section pieces in points used.    
      ELSEIF(ivar.EQ.2) THEN    
        tau=vintpt(iacc,11) 
        taup=vintpt(iacc,16)    
        taupmn=vintpt(iacc,6)   
        taupmx=vintpt(iacc,26)  
        ataup1=log(taupmx/taupmn)   
        ataup2=((1.-tau/taupmx)**4-(1.-tau/taupmn)**4)/(4.*tau) 
        wtmat(ibin,1)=wtmat(ibin,1)+1.  
        wtmat(ibin,2)=wtmat(ibin,2)+(ataup1/ataup2)*(1.-tau/taup)**3/   
     &  taup    
    
C...Sum up y* and cos(theta-hat) cross-section pieces in points used.   
      ELSEIF(ivar.EQ.3) THEN    
        yst=vintpt(iacc,12) 
        wtmat(ibin,1)=wtmat(ibin,1)+(ayst0/ayst1)*(yst-ystmin)  
        wtmat(ibin,2)=wtmat(ibin,2)+(ayst0/ayst1)*(ystmax-yst)  
        wtmat(ibin,3)=wtmat(ibin,3)+(ayst0/ayst3)/cosh(yst) 
      ELSE  
        rm34=2.*sqm3*sqm4/(vintpt(iacc,11)*vint(2))**2  
        rsqm=1.+rm34    
        cthmax=sqrt(1.-4.*vint(71)**2/(taumax*vint(2))) 
        cthmin=-cthmax  
        IF(cthmax.GT.0.9999) rm34=max(rm34,2.*vint(71)**2/  
     &  (taumax*vint(2)))   
        acth1=cthmax-cthmin 
        acth2=log(max(rm34,rsqm-cthmin)/max(rm34,rsqm-cthmax))  
        acth3=log(max(rm34,rsqm+cthmax)/max(rm34,rsqm+cthmin))  
        acth4=1./max(rm34,rsqm-cthmax)-1./max(rm34,rsqm-cthmin) 
        acth5=1./max(rm34,rsqm+cthmin)-1./max(rm34,rsqm+cthmax) 
        cth=vintpt(iacc,13) 
        wtmat(ibin,1)=wtmat(ibin,1)+1.  
        wtmat(ibin,2)=wtmat(ibin,2)+(acth1/acth2)/max(rm34,rsqm-cth)    
        wtmat(ibin,3)=wtmat(ibin,3)+(acth1/acth3)/max(rm34,rsqm+cth)    
        wtmat(ibin,4)=wtmat(ibin,4)+(acth1/acth4)/max(rm34,rsqm-cth)**2 
        wtmat(ibin,5)=wtmat(ibin,5)+(acth1/acth5)/max(rm34,rsqm+cth)**2 
      ENDIF 
  140 CONTINUE  
    
C...Check that equation system solvable; else trivial way out.  
      IF(mstp(122).GE.2) WRITE(mstu(11),1300) cvar(ivar)    
      msolv=1   
      DO 150 ibin=1,nbin    
      IF(mstp(122).GE.2) WRITE(mstu(11),1400) (wtmat(ibin,ired),    
     &ired=1,nbin),wtrel(ibin)  
  150 IF(narel(ibin).EQ.0) msolv=0  
      IF(msolv.EQ.0) THEN   
        DO 160 ibin=1,nbin  
  160   coefu(ibin)=1.  
    
C...Solve to find relative importance of cross-section pieces.  
      ELSE  
        DO 170 ired=1,nbin-1    
        DO 170 ibin=ired+1,nbin 
        rqt=wtmat(ibin,ired)/wtmat(ired,ired)   
        wtrel(ibin)=wtrel(ibin)-rqt*wtrel(ired) 
        DO 170 icoe=ired,nbin   
  170   wtmat(ibin,icoe)=wtmat(ibin,icoe)-rqt*wtmat(ired,icoe)  
        DO 190 ired=nbin,1,-1   
        DO 180 icoe=ired+1,nbin 
  180   wtrel(ired)=wtrel(ired)-wtmat(ired,icoe)*coefu(icoe)    
  190   coefu(ired)=wtrel(ired)/wtmat(ired,ired)    
      ENDIF 
    
C...Normalize coefficients, with piece shared democratically.   
      coefsu=0. 
      DO 200 ibin=1,nbin    
      coefu(ibin)=max(0.,coefu(ibin))   
  200 coefsu=coefsu+coefu(ibin) 
      IF(ivar.EQ.1) ioff=0  
      IF(ivar.EQ.2) ioff=14 
      IF(ivar.EQ.3) ioff=6  
      IF(ivar.EQ.4) ioff=9  
      IF(coefsu.GT.0.) THEN 
        DO 210 ibin=1,nbin  
  210   coef(isub,ioff+ibin)=parp(121)/nbin+(1.-parp(121))*coefu(ibin)/ 
     &  coefsu  
      ELSE  
        DO 220 ibin=1,nbin  
  220   coef(isub,ioff+ibin)=1./nbin    
      ENDIF 
      IF(mstp(122).GE.2) WRITE(mstu(11),1500) cvar(ivar),   
     &(coef(isub,ioff+ibin),ibin=1,nbin)    
  230 CONTINUE  
    
C...Find two most promising maxima among points previously determined.  
      DO 240 j=1,4  
      iaccmx(j)=0   
  240 sigsmx(j)=0.  
      nmax=0    
      DO 290 iacc=1,nacc    
      DO 250 j=1,30 
  250 vint(10+j)=vintpt(iacc,j) 
      CALL pyhisigh(nchn,sigs)    
      ieq=0 
      DO 260 imv=1,nmax 
  260 IF(abs(sigs-sigsmx(imv)).LT.1e-4*(sigs+sigsmx(imv))) ieq=imv  
      IF(ieq.EQ.0) THEN 
        DO 270 imv=nmax,1,-1    
        iin=imv+1   
        IF(sigs.LE.sigsmx(imv)) goto 280    
        iaccmx(imv+1)=iaccmx(imv)   
  270   sigsmx(imv+1)=sigsmx(imv)   
        iin=1   
  280   iaccmx(iin)=iacc    
        sigsmx(iin)=sigs    
        IF(nmax.LE.1) nmax=nmax+1   
      ENDIF 
  290 CONTINUE  
    
C...Read out starting position for search.  
      IF(mstp(122).GE.2) WRITE(mstu(11),1600)   
      sigsam=sigsmx(1)  
      DO 330 imax=1,nmax    
      iacc=iaccmx(imax) 
      mtau=mvarpt(iacc,1)   
      mtaup=mvarpt(iacc,2)  
      myst=mvarpt(iacc,3)   
      mcth=mvarpt(iacc,4)   
      vtau=0.5  
      vyst=0.5  
      vcth=0.5  
      vtaup=0.5 
    
C...Starting point and step size in parameter space.    
      DO 320 irpt=1,2   
      DO 310 ivar=1,4   
      IF(npts(ivar).EQ.1) goto 310  
      IF(ivar.EQ.1) vvar=vtau   
      IF(ivar.EQ.2) vvar=vtaup  
      IF(ivar.EQ.3) vvar=vyst   
      IF(ivar.EQ.4) vvar=vcth   
      IF(ivar.EQ.1) mvar=mtau   
      IF(ivar.EQ.2) mvar=mtaup  
      IF(ivar.EQ.3) mvar=myst   
      IF(ivar.EQ.4) mvar=mcth   
      IF(irpt.EQ.1) vdel=0.1    
      IF(irpt.EQ.2) vdel=max(0.01,min(0.05,vvar-0.02,0.98-vvar))    
      IF(irpt.EQ.1) vmar=0.02   
      IF(irpt.EQ.2) vmar=0.002  
      imov0=1   
      IF(irpt.EQ.1.AND.ivar.EQ.1) imov0=0   
      DO 300 imov=imov0,8   
    
C...Define new point in parameter space.    
      IF(imov.EQ.0) THEN    
        inew=2  
        vnew=vvar   
      ELSEIF(imov.EQ.1) THEN    
        inew=3  
        vnew=vvar+vdel  
      ELSEIF(imov.EQ.2) THEN    
        inew=1  
        vnew=vvar-vdel  
      ELSEIF(sigssm(3).GE.max(sigssm(1),sigssm(2)).AND. 
     &vvar+2.*vdel.LT.1.-vmar) THEN 
        vvar=vvar+vdel  
        sigssm(1)=sigssm(2) 
        sigssm(2)=sigssm(3) 
        inew=3  
        vnew=vvar+vdel  
      ELSEIF(sigssm(1).GE.max(sigssm(2),sigssm(3)).AND. 
     &vvar-2.*vdel.GT.vmar) THEN    
        vvar=vvar-vdel  
        sigssm(3)=sigssm(2) 
        sigssm(2)=sigssm(1) 
        inew=1  
        vnew=vvar-vdel  
      ELSEIF(sigssm(3).GE.sigssm(1)) THEN   
        vdel=0.5*vdel   
        vvar=vvar+vdel  
        sigssm(1)=sigssm(2) 
        inew=2  
        vnew=vvar   
      ELSE  
        vdel=0.5*vdel   
        vvar=vvar-vdel  
        sigssm(3)=sigssm(2) 
        inew=2  
        vnew=vvar   
      ENDIF 
    
C...Convert to relevant variables and find derived new limits.  
      IF(ivar.EQ.1) THEN    
        vtau=vnew   
        CALL pyhikmap(1,mtau,vtau)    
        IF(istsb.EQ.3.OR.istsb.EQ.4) CALL pyhiklim(4) 
      ENDIF 
      IF(ivar.LE.2.AND.(istsb.EQ.3.OR.istsb.EQ.4)) THEN 
        IF(ivar.EQ.2) vtaup=vnew    
        CALL pyhikmap(4,mtaup,vtaup)  
      ENDIF 
      IF(ivar.LE.2) CALL pyhiklim(2)  
      IF(ivar.LE.3) THEN    
        IF(ivar.EQ.3) vyst=vnew 
        CALL pyhikmap(2,myst,vyst)    
        CALL pyhiklim(3)  
      ENDIF 
      IF(istsb.EQ.2.OR.istsb.EQ.4) THEN 
        IF(ivar.EQ.4) vcth=vnew 
        CALL pyhikmap(3,mcth,vcth)    
      ENDIF 
      IF(isub.EQ.96) vint(25)=vint(21)*(1.-vint(23)**2) 
    
C...Evaluate cross-section. Save new maximum. Final maximum.    
      CALL pyhisigh(nchn,sigs)    
      sigssm(inew)=sigs 
      IF(sigs.GT.sigsam) sigsam=sigs    
      IF(mstp(122).GE.2) WRITE(mstu(11),1700) imax,ivar,mvar,imov,  
     &vnew,vint(21),vint(22),vint(23),vint(26),sigs 
  300 CONTINUE  
  310 CONTINUE  
  320 CONTINUE  
      IF(imax.EQ.1) sigs11=sigsam   
  330 CONTINUE  
      xsec(isub,1)=1.05*sigsam  
  340 IF(isub.NE.96) xsec(0,1)=xsec(0,1)+xsec(isub,1)   
  350 CONTINUE  
    
C...Print summary table.    
      IF(mstp(122).GE.1) THEN   
        WRITE(mstu(11),1800)    
        WRITE(mstu(11),1900)    
        DO 360 isub=1,200   
        IF(msub(isub).NE.1.AND.isub.NE.96) goto 360 
        IF(isub.EQ.96.AND.mint(43).NE.4) goto 360   
        IF(isub.EQ.96.AND.msub(95).NE.1.AND.mstp(81).LE.0) goto 360 
        IF(msub(95).EQ.1.AND.(isub.EQ.11.OR.isub.EQ.12.OR.isub.EQ.13.OR.    
     &  isub.EQ.28.OR.isub.EQ.53.OR.isub.EQ.68)) goto 360   
        WRITE(mstu(11),2000) isub,proc(isub),xsec(isub,1)   
  360   CONTINUE    
        WRITE(mstu(11),2100)    
      ENDIF 
    
C...Format statements for maximization results. 
 1000 FORMAT(/1x,'Coefficient optimization and maximum search for ',    
     &'subprocess no',i4/1x,'Coefficient modes     tau',10x,'y*',9x,    
     &'cth',9x,'tau''',7x,'sigma')  
 1100 FORMAT(1x,4i4,f12.8,f12.6,f12.7,f12.8,1p,e12.4)   
 1200 FORMAT(1x,'Error: requested subprocess ',i3,' has vanishing ',    
     &'cross-section.'/1x,'Execution stopped!')
 1300 FORMAT(1x,'Coefficients of equation system to be solved for ',a4) 
 1400 FORMAT(1x,1p,7e11.3)  
 1500 FORMAT(1x,'Result for ',a4,':',6f9.4) 
 1600 FORMAT(1x,'Maximum search for given coefficients'/2x,'MAX VAR ',  
     &'MOD MOV   VNEW',7x,'tau',7x,'y*',8x,'cth',7x,'tau''',7x,'sigma') 
 1700 FORMAT(1x,4i4,f8.4,f11.7,f9.3,f11.6,f11.7,1p,e12.4)   
 1800 FORMAT(/1x,8('*'),1x,'PYHIMAXI: summary of differential ',  
     &'cross-section maximum search',1x,8('*')) 
 1900 FORMAT(/11x,58('=')/11x,'I',38x,'I',17x,'I'/11x,'I  ISUB  ',  
     &'Subprocess name',15x,'I  Maximum value  I'/11x,'I',38x,'I',  
     &17x,'I'/11x,58('=')/11x,'I',38x,'I',17x,'I')  
 2000 FORMAT(11x,'I',2x,i3,3x,a28,2x,'I',2x,1p,e12.4,3x,'I')    
 2100 FORMAT(11x,'I',38x,'I',17x,'I'/11x,58('='))   
    
      RETURN    
      END