pyplot.f

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

 
C*********************************************************************
 
C...PYPLOT
C...Prints a histogram (but does not reset it).
 
      SUBROUTINE pyplot(ID)
 
C...Double precision declaration.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
C...Commonblocks.
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pybins/ihist(4),indx(1000),bin(20000)
      SAVE /pydat1/,/pybins/
C...Local arrays and character variables.
      dimension idati(6), irow(100), ifra(100), dyac(10)
      CHARACTER title*60, out*100, cha(0:11)*1
 
C...Steps in histogram scale. Character sequence.
      DATA dyac/.04,.05,.06,.08,.10,.12,.15,.20,.25,.30/
      DATA cha/'0','1','2','3','4','5','6','7','8','9','X','-'/
 
C...Find initial address in memory; skip if empty histogram.
      IF(id.LE.0.OR.id.GT.ihist(1)) RETURN
      is=indx(id)
      IF(is.EQ.0) RETURN
      IF(nint(bin(is+5)).LE.0) THEN
        WRITE(mstu(11),5000) id
        RETURN
      ENDIF
 
C...Number of histogram lines and x bins.
      lin=ihist(3)-18
      nx=nint(bin(is+1))
 
C...Extract title by conversion from double precision via integer.
      DO 100 it=1,20
        ieq=nint(bin(is+8+nx+it))
        title(3*it-2:3*it)=char(ieq/256**2)//char(mod(ieq,256**2)/256)
     &  //char(mod(ieq,256))
  100 CONTINUE
 
C...Find time; print title.
      CALL pytime(idati)
      IF(idati(1).GT.0) THEN
        WRITE(mstu(11),5100) id, title, (idati(j),j=1,5)
      ELSE
        WRITE(mstu(11),5200) id, title
      ENDIF
 
C...Find minimum and maximum bin content.
      ymin=bin(is+9)
      ymax=bin(is+9)
      DO 110 ix=is+10,is+8+nx
        IF(bin(ix).LT.ymin) ymin=bin(ix)
        IF(bin(ix).GT.ymax) ymax=bin(ix)
  110 CONTINUE
 
C...Determine scale and step size for y axis.
      IF(ymax-ymin.GT.lin*dyac(1)*1d-9) THEN
        IF(ymin.GT.0d0.AND.ymin.LT.0.1d0*ymax) ymin=0d0
        IF(ymax.LT.0d0.AND.ymax.GT.0.1d0*ymin) ymax=0d0
        ipot=int(log10(ymax-ymin)+10d0)-10
        IF(ymax-ymin.LT.lin*dyac(1)*10d0**ipot) ipot=ipot-1
        IF(ymax-ymin.GT.lin*dyac(10)*10d0**ipot) ipot=ipot+1
        dely=dyac(1)
        DO 120 idel=1,9
          IF(ymax-ymin.GE.lin*dyac(idel)*10d0**ipot) dely=dyac(idel+1)
  120   CONTINUE
        dy=dely*10d0**ipot
 
C...Convert bin contents to integer form; fractional fill in top row.
        DO 130 ix=1,nx
          cta=abs(bin(is+8+ix))/dy
          irow(ix)=sign(cta+0.95d0,bin(is+8+ix))
          ifra(ix)=10d0*(cta+1.05d0-dble(int(cta+0.95d0)))
  130   CONTINUE
        irmi=sign(abs(ymin)/dy+0.95d0,ymin)
        irma=sign(abs(ymax)/dy+0.95d0,ymax)
 
C...Print histogram row by row.
        DO 150 ir=irma,irmi,-1
          IF(ir.EQ.0) goto 150
          out=' '
          DO 140 ix=1,nx
            IF(ir.EQ.irow(ix)) out(ix:ix)=cha(ifra(ix))
            IF(ir*(irow(ix)-ir).GT.0) out(ix:ix)=cha(10)
  140     CONTINUE
          WRITE(mstu(11),5300) ir*dely, ipot, out
  150   CONTINUE
 
C...Print sign and value of bin contents.
        ipot=int(log10(max(ymax,-ymin))+10.0001d0)-10
        out=' '
        DO 160 ix=1,nx
          IF(bin(is+8+ix).LT.-10d0**(ipot-4)) out(ix:ix)=cha(11)
          irow(ix)=nint(10d0**(3-ipot)*abs(bin(is+8+ix)))
  160   CONTINUE
        WRITE(mstu(11),5400) out
        DO 180 ir=4,1,-1
          DO 170 ix=1,nx
            out(ix:ix)=cha(mod(irow(ix),10**ir)/10**(ir-1))
  170     CONTINUE
          WRITE(mstu(11),5500) ipot+ir-4, out
  180   CONTINUE
 
C...Print sign and value of lower bin edge.
        ipot=int(log10(max(-bin(is+2),bin(is+3)-bin(is+4)))+
     &  10.0001d0)-10
        out=' '
        DO 190 ix=1,nx
          IF(bin(is+2)+(ix-1)*bin(is+4).LT.-10d0**(ipot-3))
     &    out(ix:ix)=cha(11)
          irow(ix)=nint(10d0**(2-ipot)*abs(bin(is+2)+(ix-1)*bin(is+4)))
  190   CONTINUE
        WRITE(mstu(11),5600) out
        DO 210 ir=3,1,-1
          DO 200 ix=1,nx
            out(ix:ix)=cha(mod(irow(ix),10**ir)/10**(ir-1))
  200     CONTINUE
          WRITE(mstu(11),5500) ipot+ir-3, out
  210   CONTINUE
      ENDIF
 
C...Calculate and print statistics.
      csum=0d0
      cxsum=0d0
      cxxsum=0d0
      DO 220 ix=1,nx
        cta=abs(bin(is+8+ix))
        x=bin(is+2)+(ix-0.5d0)*bin(is+4)
        csum=csum+cta
        cxsum=cxsum+cta*x
        cxxsum=cxxsum+cta*x**2
  220 CONTINUE
      xmean=cxsum/max(csum,1d-20)
      xrms=sqrt(max(0d0,cxxsum/max(csum,1d-20)-xmean**2))
      WRITE(mstu(11),5700) nint(bin(is+5)),xmean,bin(is+6),
     &bin(is+2),bin(is+7),xrms,bin(is+8),bin(is+3)
 
C...Formats for output.
 5000 FORMAT(/5x,'Histogram no',i5,' : no entries')
 5100 FORMAT('1'/5x,'Histogram no',i5,6x,a60,5x,i4,'-',i2,'-',i2,1x,
     &i2,':',i2/)
 5200 FORMAT('1'/5x,'Histogram no',i5,6x,a60/)
 5300 FORMAT(2x,f7.2,'*10**',i2,3x,a100)
 5400 FORMAT(/8x,'Contents',3x,a100)
 5500 FORMAT(9x,'*10**',i2,3x,a100)
 5600 FORMAT(/8x,'Low edge',3x,a100)
 5700 FORMAT(/5x,'Entries  =',i12,1p,6x,'Mean =',d12.4,6x,'Underflow ='
     &,d12.4,6x,'Low edge  =',d12.4/5x,'All chan =',d12.4,6x,
     &'Rms  =',d12.4,6x,'Overflow  =',d12.4,6x,'High edge =',d12.4)
 
      RETURN
      END