pyoper.f

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

 
C*********************************************************************
 
C...PYOPER
C...Performs operations between histograms.
 
      SUBROUTINE pyoper(ID1,OPER,ID2,ID3,F1,F2)
 
C...Double precision declaration.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
C...Commonblock.
      common/pybins/ihist(4),indx(1000),bin(20000)
      SAVE /pybins/
C...Character variable.
      CHARACTER oper*(*)
 
C...Find initial addresses in memory, and histogram size.
      IF(id1.LE.0.OR.id1.GT.ihist(1)) CALL pyerrm(28,
     &'(PYFACT:) not allowed histogram number')
      is1=indx(id1)
      is2=indx(min(ihist(1),max(1,id2)))
      is3=indx(min(ihist(1),max(1,id3)))
      nx=nint(bin(is3+1))
      IF(oper.EQ.'M'.AND.id3.EQ.0) nx=nint(bin(is2+1))
 
C...Update info on number of histogram entries.
      IF(oper.EQ.'+'.OR.oper.EQ.'-'.OR.oper.EQ.'*'.OR.oper.EQ.'/') THEN
        bin(is3+5)=bin(is1+5)+bin(is2+5)
      ELSEIF(oper.EQ.'A'.OR.oper.EQ.'S'.OR.oper.EQ.'L') THEN
        bin(is3+5)=bin(is1+5)
      ENDIF
 
C...Operations on pair of histograms: addition, subtraction,
C...multiplication, division.
      IF(oper.EQ.'+') THEN
        DO 100 ix=6,8+nx
          bin(is3+ix)=f1*bin(is1+ix)+f2*bin(is2+ix)
  100   CONTINUE
      ELSEIF(oper.EQ.'-') THEN
        DO 110 ix=6,8+nx
          bin(is3+ix)=f1*bin(is1+ix)-f2*bin(is2+ix)
  110   CONTINUE
      ELSEIF(oper.EQ.'*') THEN
        DO 120 ix=6,8+nx
          bin(is3+ix)=f1*bin(is1+ix)*f2*bin(is2+ix)
  120   CONTINUE
      ELSEIF(oper.EQ.'/') THEN
        DO 130 ix=6,8+nx
          fa2=f2*bin(is2+ix)
          IF(abs(fa2).LE.1d-20) THEN
            bin(is3+ix)=0d0
          ELSE
            bin(is3+ix)=f1*bin(is1+ix)/fa2
          ENDIF
  130   CONTINUE
 
C...Operations on single histogram: multiplication+addition,
C...square root+addition, logarithm+addition.
      ELSEIF(oper.EQ.'A') THEN
        DO 140 ix=6,8+nx
          bin(is3+ix)=f1*bin(is1+ix)+f2
  140   CONTINUE
      ELSEIF(oper.EQ.'S') THEN
        DO 150 ix=6,8+nx
          bin(is3+ix)=f1*sqrt(max(0d0,bin(is1+ix)))+f2
  150   CONTINUE
      ELSEIF(oper.EQ.'L') THEN
        zmin=1d20
        DO 160 ix=9,8+nx
          IF(bin(is1+ix).LT.zmin.AND.bin(is1+ix).GT.1d-20)
     &    zmin=0.8d0*bin(is1+ix)
  160   CONTINUE
        DO 170 ix=6,8+nx
          bin(is3+ix)=f1*log10(max(zmin,bin(is1+ix)))+f2
  170   CONTINUE
 
C...Operation on two or three histograms: average and
C...standard deviation.
      ELSEIF(oper.EQ.'M') THEN
        DO 180 ix=6,8+nx
          IF(abs(bin(is1+ix)).LE.1d-20) THEN
            bin(is2+ix)=0d0
          ELSE
            bin(is2+ix)=bin(is2+ix)/bin(is1+ix)
          ENDIF
          IF(id3.NE.0) THEN
            IF(abs(bin(is1+ix)).LE.1d-20) THEN
              bin(is3+ix)=0d0
            ELSE
              bin(is3+ix)=sqrt(max(0d0,bin(is3+ix)/bin(is1+ix)-
     &        bin(is2+ix)**2))
            ENDIF
          ENDIF
          bin(is1+ix)=f1*bin(is1+ix)
  180   CONTINUE
      ENDIF
 
      RETURN
      END