*     RUNPAW file with copy of commons
 
 
 
      SUBROUTINE RHBOOK
*
*     interface between RUN and HBOOK for graphical display
*
*     UNCOMA: ntuple data accesible to user ----------------------------
      COMMON/UNCOMA/IDMB,NY,WT,X,Y(100)
      REAL REC(102),XY(0:100)
      EQUIVALENCE (REC(1),WT),(XY(0),X)
*     IDMB = flag = 1 for Data, =2 for Monte Carlo, =3 for Background
*     NY   = number of measured variables
*     WT   = weight
*     X    = true variable (defined only in Monte Carlo) = XY(0)
*     Y(j) = measured variables, j = 1 ... NY            = XY(j)
*     REC(.) = array with: WT, X, Y(.) (by equivalence)
*     ==================================================================
*     UNCOMB: internal variables and arrays ----------------------------
      PARAMETER (LENGDP=50000)
      PARAMETER (NKNMAX=63,NHU=4*(NKNMAX-3))
      COMMON/UNCOMB/HUH(NHU),HUS(NHU),HUF(NHU),HUW(NHU),HUA(NHU),
     +     HUP(NHU),ASPMC(NKNMAX),
     +     NHA,ITAB,NTAB,JTAB(12,33),INWT,INIX,LUNP,PERC,LPFLAG,
     +     XYMIN(0:32),XYMAX(0:32),JXY(0:32),LENGD,
     +     NKN,NDF,XLOW,XHIG,NCNS,FCNS(2,10),NBDEF,BINS(31),NIKN,
     +    WLIN,LOPT(10),NOPT,NV,NDV(3),NBV(3),XYFL(3),XYFH(3),
     +    NPROD,NBN,NSPACE,NDFA,FACLUM,LUSFUN,GAMMA,DATLUM,XMCLUM,
     +    HIST(120,3,32),HSM(120),XYL(0:32),XYH(0:32)
      REAL STAB(12,33)
      EQUIVALENCE (JTAB(1,1),STAB(1,1))
*
*     ITAB         = index of last ntuple
*     NTAB         = number of different ntuple types
*     J/STAB( 1,.) = 1 data  =2 moca  =3 bg
*     J/STAB( 2,.) = hollerith text
*     J/STAB( 3,.) =
*     J/STAB( 4,.) = Lumi  Lumi  evtnr    (first argument
*     J/STAB( 5,.) =
*     J/STAB( 6,.) =
*     J/STAB( 7,.) =
*     J/STAB( 8,.) =
*     J/STAB( 9,.) = 1.0   wfac2           (second argument)
*     J/STAB(10,.) =
*     J/STAB(11,.) = ifg for LOOK ntuples
*     J/STAB(12,.) = inr for LOOK ntuples
*
*     INWT = index weight
*     INIX = index X
*
*
*     NV = number of fit variables                (from text input)
*     NDV(j) = index of fit variable                   - " -
*     NBV(j) = number of bins                          - " -
*
*     XLOW, XHIG = lower and upper limit of x     (from text input)
*        mandatory, not modified by program
*
*     XMIN(j),XMAX(j) = min and max value of Y(j), determined from
*        ntuples with positive weight (except background)
*     JXY(j) = 0 if XMIN/XMAX undefined
*     JXY(j) = 1 if XMIN/XMAX defined
*
*     XYL(j),XYH(j) = limits of histogram, rounded for 120 bins,
*                     for final histograms
*
*     XYFL(K),XYFH(K) = limits as used by unfolding fit
*
*
*     Where are bin limits used:
*
*     for unfolding fit, usually with small number of bins: XYFL/XYFH
*
*     for final histograms, with 120 bins: XYL/XYH
*
*
*
*
*     HIST(.,.,.) = histogram of meaured variables, filled in PASS 3
*     ==================================================================
*     UNCOMC: character strings-----------------------------------------
      CHARACTER*50 TITLE,VTEXT
      CHARACTER*8 SOPTE
      COMMON/UNCOMC/TITLE,VTEXT(0:32),SOPTE(10)
*     TITLE    =  text as comment for problem
*     VTEXT(j) =  text as comment for variable j
*     SOPTE(i) =  selected options
*     ------------------------------------------------------------------
      SAVE /UNCOMA/,/UNCOMB/,/UNCOMC/
*     ==================================================================
      COMMON/RUNRES/XW(527),RSP(2082)
*     ------------------------------------------------------------------
*     Internal common for matrices, vectors etc for up to 63 knots
*     and up to 30 data points
      PARAMETER (NKMAX=63,NDAMAX=30)
      COMMON/UINTER/UM(NKMAX*NKMAX),HE((NKMAX+NKMAX*NKMAX)/2),
     +              GR(NKMAX),XX(NKMAX),AM(NKMAX),DE(NKMAX),
     +              SCV(NKMAX),
     +              SCM(NKMAX*NKMAX),
     +              FUN,SF,SD,TAU
*     M=NKN-1
*     UM = tranformation matrix
*     HE = hessian, replaced by covariance matrix
*     GR = gradient, replaced by step vector
*     XX = coefficient
*     AM = amplitudes     M
*     DE = eigenvalues    M
*     SC = scratch storage
*     ------------------------------------------------------------------
      SAVE /RUNRES/,/UINTER/
 
      CHARACTER*2 TNUM
      EXTERNAL HBFUN
      COMMON/CHBFUN/NXFUN,XAFUN,XBFUN,HCONT(1000)
      COMMON/UNCOMD/H(LENGDP)
      INTEGER NH(10000)
      EQUIVALENCE (H(1),NH(1))
      REAL RH(11,2),SH(2),S(4),RSREL(2),RSNUM(2)
      REAL BX(2),BY(2)
      REAL HUT(120),HUQ(120),HISTO(200)
      CHARACTER*7  COLOR(7), CHNUM*3
      DOUBLE PRECISION SUMDB,SUMMC
      DATA COLOR/'WHITE','RED','GREEN','BLUE','YELLOW','MAGENTA','CYAN'/
*     ...
      CALL SMTEXT(' ')
      CALL SMTEXT('-> Subroutine RHBOOK starting ...')
*
*     ------------------------------------------------------------------
*     plots of group 1: Result of unfolding and related distribution
*
*
*     Unfolding result as curve
      NKN=RSP(1)+0.5
      XLOW=RSP(2)
      XHIG=RSP(3)
 
      K=XW(1)+1.5
*     Data points from unfolding
      NX=200
      CALL HBOOK1(101,'Fig.101  Unfolding result',NX,XLOW,XHIG,0.0)
      CALL HBOOK1(102,'Fig.102  Unfolding result',NX,XLOW,XHIG,0.0)
 
      DO I=1,NX
       HISTO(I)=0.0
      END DO
      NN=XW(1)+0.5
      II=0
      DO I=1,NN
       II=II+I
       DY=SQRT(XW(NN+NN+2+II))
       YY=XW(NN+2+I)
       XE=XW(1+I)
       XR=XW(2+I)
*      calculate bins
       ZE=1.0+(XE-XLOW)/(XHIG-XLOW)*FLOAT(NX)
       IE=ZE
       IF(ZE-FLOAT(IE).GT.0.5) IE=IE+1
*
       ZR=1.0+(XR-XLOW)/(XHIG-XLOW)*FLOAT(NX)
       IR=ZR
       IF(ZR-FLOAT(IR).GT.0.5) IR=IR+1
*      all bins to same value
       DO J=IE,IR
        HISTO(J)=YY
       END DO
      END DO
      CALL HPAK(101,HISTO)
      CALL HPAK(102,HISTO)
*     curves
      NXFUN=NX
      XAFUN=XLOW
      XBFUN=XHIG
      DO I=1,NX
       XF=(FLOAT(NX+NX-I-I+1)*XLOW+FLOAT(I+I-1)*XHIG)/FLOAT(NX+NX)
       FU=ESPF(XF,RSP(4),NKN,XLOW,XHIG)
       TX=FU*USFUN(XF)
       HCONT(I)=FACLUM*TX
      END DO
      CALL HFUNC(102,HBFUN)
      DO I=1,NX
       XF=(FLOAT(NX+NX-I-I+1)*XLOW+FLOAT(I+I-1)*XHIG)/FLOAT(NX+NX)
       TT=ESPF(XF,UM(1+(K-1)*NKN),NKN,XLOW,XHIG)*USFUN(XF)
       HCONT(I)=FACLUM*TT
      END DO
      CALL HBFUN1(103,'Fig.103  (N+1).th orthogonal function',
     +            NX,XLOW,XHIG,HBFUN)
*     figure 104: function USFUN
      NX=200
      DO I=1,NX
*      X in center of bin
       XF=(FLOAT(NX+NX-I-I+1)*XLOW+FLOAT(I+I-1)*XHIG)/FLOAT(NX+NX)
       TX=USFUN(XF)
       HCONT(I)=TX
      END DO
      NXFUN=NX
      XAFUN=XLOW
      XBFUN=XHIG
      CALL HBFUN1(104,'Fig.104  User function USFUN(X)',
     +     200,XLOW,XHIG,HBFUN)
*     figures 105 and 106
      NX=NHA
      NXFUN=NX
*     book ...
      CALL HBOOK1(105,
     +     'Fig.105  MOCA input distribution and spline fit',
     +     NX,XLOW,XHIG,0.0)
*     ... and global fill
      CALL HPAK(105,HUH)
      DO I=1,NX
*      X in center of bin
       XF=(FLOAT(NX+NX-I-I+1)*XLOW+FLOAT(I+I-1)*XHIG)/FLOAT(NX+NX)
       TX=ESPF(XF,ASPMC,NKN,XLOW,XHIG)
       HCONT(I)=TX/GAMMA
      END DO
      CALL HFUNC(105,HBFUN)
 
      NX=200
      CALL HBOOK1(106,
     +    'Fig.106  Unfolding distribution f mult(x)',
     +     NX,XLOW,XHIG,0.0)
      DO I=1,NX
*      X in center of bin
       XF=(FLOAT(NX+NX-I-I+1)*XLOW+FLOAT(I+I-1)*XHIG)/FLOAT(NX+NX)
       FU=ESPF(XF,RSP(4),NKN,XLOW,XHIG)
       HCONT(I)=FACLUM*FU
      END DO
      NXFUN=NX
      CALL HFUNC(106,HBFUN)
 
*     figure 107: Weight as a function of x
 
      CALL HBOOK1(107,'Fig.107  Weight as a function of x',
     +            NHA,XLOW,XHIG,0.0)
      CALL HPAK(107,HUW)
*
*     ------------------------------------------------------------------
*     Plots of group 2: Orthogonal functions
*
      NX=NKN*4
*     ... is number of bins for function plot
      DO K=1,NIKN
*      calculate function first
       DO I=1,NX
*       X in center of bin
        X=(FLOAT(NX+NX-I-I+1)*XLOW+FLOAT(I+I-1)*XHIG)/FLOAT(NX+NX)
        F=ESPF(X,UM(1+(K-1)*NKN),NKN,XLOW,XHIG)
*       fill into common array for HBOOK function
        HCONT(I)=F
       END DO
       TNUM='  '
       IF(K.GE.10) TNUM(1:1)=CHAR(ICHAR('0')+K/10)
       TNUM(2:2)=CHAR(ICHAR('0')+MOD(K,10))
*      information for function
       NXFUN=NX
       XAFUN=XLOW
       XBFUN=XHIG
       CALL HBFUN1(200+K,TNUM//'. orthogonal function ',
     +             NX,XLOW,XHIG,HBFUN)
      END DO
*
*     ------------------------------------------------------------------
*     Plots of group 3: Fit variables and B-splines
*
*     amplitudes
      NX=NKN-1
      NX=NIKN
      CALL HBOOK1(301,
     + 'Fig.301  Amplitudes (Cut at 12) and 10*damping factor',
     +            NX,0.5,FLOAT(NX)+0.5,0.0)
      DO I=1,NX
       HCONT(I)=MIN(12.0,ABS(AM(I)))
      END DO
      CALL HPAK(301,HCONT)
*     curve of damping factor *10.0
      DO I=1,NX
C      HCONT(I)=1.0
C      IF(DE(I).NE.0.0) THEN
          HCONT(I)=10.0/(1.0+TAU*DE(I))
C      END IF
      END DO
      NXFUN=NX
      XAFUN=0.5
      XBFUN=0.5+FLOAT(NX)
      CALL HFUNC(301,HBFUN)
 
*     ------------------------------------------------------------------
*     Plots of group 4: MOCA distributions and B-spline contributions
*
*     not implemented for HBOOK
*
*     Plots of groups 5 to 8 -------------------------------------------
*
*     Groups 6 and 8 not implemented for HBOOK
      DO J=1,32
       IF(JXY(J).NE.0) THEN
          TNUM='  '
          IF(J/10.NE.0) TNUM(1:1)=CHAR(ICHAR('0')+J/10)
          TNUM(2:2)=CHAR(ICHAR('0')+MOD(J,10))
*         length of text for variable
          DO L=50,1,-1
          IF(VTEXT(J)(L:L).NE.' ') GOTO 66
          END DO
          L=1
 66       LV=L
*         event numbers in bins of constant width
*         transfer data histogram to LOOK
          NX=120
          NXFUN=NX
          XAFUN=XYL(J)
          XBFUN=XYH(J)
          TNUM(1:1)=CHAR(ICHAR('0')+J/10)
*         histogram and ...
          CALL HBOOK1(500+J,
     +        'Fig.5'//TNUM//' '//VTEXT(J)(1:LV),NX,
     +         XYL(J),XYH(J),0.0)
          CALL HPAK(500+J,HIST(1,1,J))
*         copy sum of MOCA and background ...
          DO I=1,120
           HCONT(I)=HIST(I,2,J)
          END DO
*         ... smooth it ...
          CALL SPLFT(HCONT,120)
*         ... curve
          CALL HFUNC(500+J,HBFUN)
*         calculate difference
          CALL HBOOK1(700+J,
     +     'FIG.7'//TNUM//' Measured - expected for y('//TNUM//')',
     +      NX,XYL(J),XYH(J),0.0)
          DO I=1,120
           HSM(I)=HIST(I,1,J)-HIST(I,2,J)
          END DO
          CALL HPAK(700+J,HSM)
       END IF
      END DO
      CALL SMTEXT('-- Subroutine RHBOOK ending.')
      CALL SMTEXT(' ')
  100 RETURN
      END
 
 
 
      SUBROUTINE THBOOK
*
*     HBOOK test histograms
*
      COMMON/CHBFUN/NXFUN,XAFUN,XBFUN,HCONT(1000)
      EXTERNAL HBFUN
*     common for TEST histograms ---------------------------------------
      PARAMETER (NBINS=80,NDIM=10+NBINS)
      COMMON/CRTEST/FTI(101),
     + DATAXT(NDIM),DATAXA(NDIM),DATAYT(NDIM),DATAYA(NDIM),DATAYB(NDIM),
     + RANDXT(NDIM),RANDXA(NDIM),RANDYT(NDIM),RANDYA(NDIM),BACKYT(NDIM),
     + NDATA,NMOCA,NDACK,NBACK
*     end of common for TEST histograms --------------------------------
     +
*     ...
*     data histograms X
      CALL AHISTH(901,DATAXT,
     +    'Fig.901  True data x, all')
      CALL AHISTH(902,DATAXA,'Fig.902  True data x, accepted')
*     data histograms Y
      CALL AHISTH(903,DATAYT,
     +    'Fig.903  Data y, incl. background')
      CALL AHISTH(904,DATAYA,'Fig.904 Data y, aacepted ' )
      CALL AHISTH(905,DATAYB,'Fig.905  Background in data ')
*     moca histograms X
      CALL AHISTH(906,RANDXT,
     +    'Fig.906  True moca x, all')
      CALL AHISTH(907,RANDXA,'Fig.907  True MOCA x, accepted')
*     moca histograms Y
      CALL AHISTH(908,RANDYT,
     +    'Fig.908  Moca y  ')
      CALL AHISTH(909,RANDYA,'Fig. 909  MOCA y, accepted ')
*     back histograms Y
      CALL AHISTH(910,BACKYT,'Fig. 910  Background ')
*     plot histograms
 
*     integrated luminosities for data ...
      DLUMI=FLOAT(NDATA)        ! assuming total cross section of 1
*
*     plot true function as figure 911
      XA=0.0-0.01
      XB=2.0+0.01
      NX=101
      FACTOR=DLUMI*2.0/FLOAT(NBINS)
      DO I=1,NX
*      X in center of bin
       X=(FLOAT(NX+NX-I-I+1)*XA+FLOAT(I+I-1)*XB)/FLOAT(NX+NX)
       HCONT(I)=FACTOR*FTI(I)
      END DO
      NXFUN=NX
      XAFUN=XA
      XBFUN=XB
      CALL HBFUN1(911,'Fig. 911  True function ',NX,XA,XB,HBFUN)
*
  100 RETURN
      END
 
      FUNCTION HBFUN(X)
*     function to allow function plots with HBOOK
      COMMON/CHBFUN/NXFUN,XAFUN,XBFUN,HCONT(1000)
*     ...
      J=1.0+(X-XAFUN)/(XBFUN-XAFUN)*FLOAT(NXFUN)
      IF(J.LE.0.OR.J.GT.NXFUN) THEN
         HBFUN=0.0
      ELSE
         HBFUN=HCONT(J)
      END IF
      END
 
      SUBROUTINE AHISTH(IFG,A,TEXT)
*     histogram to graphic (HBOOK)
      REAL A(*)
      CHARACTER*(*) TEXT
*     ...
      NB=A(4)
      IF(A(5).EQ.0.0.OR.NB.EQ.0) RETURN
      CALL HBOOK1(IFG,TEXT,NB,A(2),A(3),0.0)
      CALL HPAK(IFG,A(11))
      END