*
* $Id$
*
* $Log$
* Revision 1.1  1996/01/16 17:07:52  mclareni
* Initial revision
*
*
#include "hbook/pilot.h"
*CMZ :  4.20/03 23/07/93  17.46.20  by  Rene Brun
*-- Author :    R. J. Genik II   23/10/92
      SUBROUTINE HDBINI(ID1,ID2,TOL,NBINS,CHOPT,ERRORS)
C----------------------------------------------------------------------
C-
C-   Purpose and Methods : Initialization routine for HDIFFB
C-                         1) set option flags and check for validity
C-                         2) set overall normalization
C-                         3) fill the local common block
C-
C-   Inputs  : IDR,IDD,TOL,NBINS,CHOPT
C-   Outputs : ERRORS, If DEBUG option on, various messages
C-              fills /HDBCOM/ the HDIFFB common block
C-   Controls: NONE
C-
C-   Created :  3-DEC-1990   James T. McKinley
C-                           Michigan State University, USA
C-
C-
C-   MODIFIED:  24-SEP-1992  R. J. Genik II
C-                           Michigan State University, USA
C-
C-
C---------------------------------------------------------------------
C
C
C
C----------------------------------------------------------------------
C     Local Variable Declaration
C----------------------------------------------------------------------
C
      INTEGER NBINS,ID1,ID2,NUMZSR,NUMZSD
      INTEGER LOCR,NCXR,NCYR,LEXCR,LEYCR,UEXCR,UEYCR,NWTR
      INTEGER LOCD,NCXD,NCYD,LEXCD,LEYCD,UEXCD,UEYCD,NWTD
      INTEGER LCIDR,LCIDD,LCONTR,LCONTD,LWR,LWD,JBIT
      INTEGER NENTR,NENTD,I,J,K,L,PROFIR,PROFID
      REAL    HGCONT,TOTALR,TOTALD,SUMR,SUMD,TOL
      LOGICAL ERRORS,NGCONR,NGCOND,ERBARR,ERBARD
      CHARACTER*80 TITLR,TITLD
      CHARACTER*(*) CHOPT
C
C
#include "hbook/hcbook.inc"
#include "hbook/hcbits.inc"
#include "hbook/hcunit.inc"
#include "hbook/hcprin.inc"
#include "hbook/hcdifb.inc"
C
C=====================================================================C
C SECTION 1  -  Initialize variables, error check, set flags          C
C=====================================================================C
C   Initialize the variables
C
      ERRORS = .FALSE.
C                                     ! No errors yet
      WEIGHR = .FALSE.
C                                     ! Histogram 1 not weighted
      WEIGHD = .FALSE.
C                                     ! 2 is not either
      TWODIM = .FALSE.
C                                     ! Assume 1-D for now
      PROFIL = .FALSE.
C                                     ! Assume not profile for now
      PSDMR  = .TRUE.
C                                     ! Assume std deviation of mean for error
      PSDMD  = .TRUE.
C                                     ! bar, if it is a profile histogram
      NGCONR = .FALSE.
C                                     ! Negative contents flag
      NGCOND = .FALSE.
C                                     ! Negative contents flag
C
C----------------------------------------------------------------------
C  Change ID1 and ID2 to IDR and IDD to keep reference and data
C  straight internally.
C----------------------------------------------------------------------
C
      IDR=ID1
      IDD=ID2
C
C----------------------------------------------------------------------
C   Set the default output device for the debug dump (DUMPDV) to the
C   common default device (LOUT)
C----------------------------------------------------------------------
C
      DUMPDV = LOUT
C
C----------------------------------------------------------------------
C   Set the default accuracy for small S stat to 2 digits
C----------------------------------------------------------------------
C
      ACDIGT = 2.
C----------------------------------------------------------------------
C   Make sure the tolerance is not 0  (Warning)
C----------------------------------------------------------------------
C
      IF (TOL .LE. 0.) CALL HBUG('+Zero tolerance ',
     +  'HDIFFB',IDR)
C
C----------------------------------------------------------------------
C   Get booking information for IDR and IDD
C----------------------------------------------------------------------
C
      CALL HGIVE(IDR, TITLR, NCXR, LEXCR, UEXCR, NCYR, LEYCR, UEYCR,
     +    NWTR, LOCR)
      CALL HGIVE(IDD, TITLD, NCXD, LEXCD, UEXCD, NCYD, LEYCD, UEYCD,
     +    NWTD, LOCD)
C
C----------------------------------------------------------------------
C   Check that the histograms are of the same dimension.  If Y channels
C   are present, set flag for two dimensions.
C----------------------------------------------------------------------
C
      IF ((NCYR.EQ.0 .OR. NCYD.EQ.0) .AND. NCYR+NCYD.NE.0) THEN
        CALL HBUG('Both histograms must be the same dimension.',
     +      'HDIFFB', IDR)
        ERRORS=.TRUE.
      ENDIF
      IF (NCXR .NE. NCXD .OR. NCYR .NE. NCYD) THEN
        CALL HBUG('Number of channels is different.', 'HDIFFB', IDR)
        ERRORS=.TRUE.
      ENDIF
      IF (NCYR .GT. 0) TWODIM=.TRUE.
C
C----------------------------------------------------------------------
C   Decode option string
C----------------------------------------------------------------------
C
      CALL HUOPTC(CHOPT,OPTST,OPTS)
C
C----------------------------------------------------------------------
C  If A option is selected require that errors bars exist.
C----------------------------------------------------------------------
C
      CALL HFIND(IDR,'HDIFFB')
      LCIDR=LCID
      ERBARR = LQ(LQ(LCIDR-1)).NE.0
      IF((OPTS(AOPTN).EQ.1).AND..NOT.ERBARR) THEN
        CALL HBUG('A option with no error bars on reference histogram',
     +    'HDIFFB',IDR)
        ERRORS = .TRUE.
      ENDIF
C
C----------------------------------------------------------------------
C   Determine whether the histograms are of profile type
C   and if so, set profile flag.
C----------------------------------------------------------------------
C
C  Get address of IDR (REF) and decode flags I1-I230 (in HCBITS)
C
      CALL HDCOFL
      PROFIR=I8
      PROFIL = (PROFIR.EQ.1)
C   Find kind of error bar used for reference profile histograms
      IF (PROFIL) THEN
        LCONTR=LQ(LCIDR-1)
        LWR=LQ(LCONTR)
        PSDMR = JBIT(IQ(LWR),1).EQ.0
C                             !True if using spread instead of SD mean
C
      ENDIF
C
C  Get address of IDD (DAT) and decode flags I1-I230 (in HCBITS)
C
      CALL HFIND(IDD,'HDIFFB')
      LCIDD=LCID
      CALL HDCOFL
      PROFID=I8
C
C
      ERBARD = LQ(LQ(LCIDD-1)).NE.0
C   Find kind of error bar used for data profile histograms
      IF (PROFIL) THEN
        LCONTD=LQ(LCIDD-1)
        LWD=LQ(LCONTD)
        PSDMD = JBIT(IQ(LWD),1).EQ.0
      ENDIF
C
C----------------------------------------------------------------------
C   Check that both histograms are profile type
C----------------------------------------------------------------------
C
C
      IF (PROFIR.NE.PROFID) THEN
        CALL HBUG('Both histograms must be standard or profile type',
     +      'HDIFFB',IDR)
        ERRORS = .TRUE.
      ENDIF
C
C----------------------------------------------------------------------
C   Display options string and TOL to user for Debug Option
C----------------------------------------------------------------------
C
C
      IF (OPTS(DEBUG) .EQ. 1) THEN
        WRITE(DUMPDV, FMT=200) IDR, IDD
        WRITE(DUMPDV, FMT=210) CHOPT, TOL
      ENDIF
C
C----------------------------------------------------------------------
C   Options N,O,U,L,R,T,B not allowed for profile histograms
C   Check parameterization for OPTS index values (numbers to save space)
C----------------------------------------------------------------------
C
      IF((PROFIR+PROFID.NE.0).AND.((OPTS(1)+OPTS(3)+OPTS(4)+OPTS(9)+
     +    OPTS(10)+OPTS(11)+OPTS(12)).NE.0))THEN
        CALL HBUG('options N,O,U,L,R,T, or B used with profile hist.',
     +             'HDIFFB',IDR)
        ERRORS=.TRUE.
      ENDIF
C
C----------------------------------------------------------------------
C   Screen out incompatible options.  S,C, or A. (Warning)
C   Set S option as default.
C----------------------------------------------------------------------
C
      IF (OPTS(SOPTN) + OPTS(COPTN) + OPTS(AOPTN) .GT. 1) THEN
        CALL HBUG('Only one comparison at a time, please','HDIFFB',
     +      IDR)
        OPTS(SOPTN) = 1
C                                     ! Default to S case
        OPTS(COPTN) = 0
C                                     ! Turn off both the
        OPTS(AOPTN) = 0
C                                     ! C and A modes.
      ENDIF
C
C
C----------------------------------------------------------------------
C   If no test options are selected, default to S
C----------------------------------------------------------------------
C
      IF (OPTS(SOPTN)+OPTS(COPTN)+OPTS(AOPTN).EQ.0) OPTS(SOPTN)=1
C
C
C----------------------------------------------------------------------
C   Make sure both histograms have the same binning.  The routines that
C   return the content and error bars don't worry about the different
C   binning, but it may cause bad results, so warn user.
C----------------------------------------------------------------------
C
      IF (LEXCR .NE. LEXCD) THEN
        CALL HBUG('+Different binning ','HDIFFB',IDR)
      ENDIF
      IF (NCYR.NE.0.AND.LEYCR .NE. LEYCD) THEN
        CALL HBUG('+Different binning ','HDIFFB',IDR)
      ENDIF
C
C
C----------------------------------------------------------------------
C   Prepare for calculation of indices in DIFFS
C----------------------------------------------------------------------
C
      BEGINI = 1
C                                     ! First bin, no underflow, X
      ENDI = NCXR
C                                     ! Last bin, no overflow, X
C
      IF (TWODIM) THEN
        BEGINJ = 1
C                                     ! First bin, no underflow, Y
        ENDJ = NCYR
C                                     ! Last bin, no overflow, Y
      ELSE
        BEGINJ = 0
C                                     ! Zero bins
        ENDJ = 0
C                                     ! Zero is also the last bin
      ENDIF
C
C
C----------------------------------------------------------------------
C   Find the starting and ending points with O/U flow bins
C----------------------------------------------------------------------
C
      IF (OPTS(UFLOW).EQ.1.OR.OPTS(XUNDR).EQ.1) BEGINI = 0
      IF (OPTS(OFLOW).EQ.1.OR.OPTS(XOVER).EQ.1) ENDI = NCXR+1
      IF ((OPTS(UFLOW).EQ.1.OR.OPTS(YUNDR).EQ.1).AND.TWODIM) BEGINJ=0
      IF((OPTS(OFLOW).EQ.1.OR.OPTS(YOVER).EQ.1).AND.TWODIM)ENDJ=NCYR+1
C
C
C----------------------------------------------------------------------
C   Find the XSIZ value (total length across the X-Axis)
C----------------------------------------------------------------------
C
      XSIZ = ENDI - BEGINI + 1
C                                   ! Used in INDEX calculation
C
C
C----------------------------------------------------------------------
C   Verify that DIFFS has enough room for the output.
C----------------------------------------------------------------------
C
      K = XSIZ*(ENDJ-BEGINJ+1)
      IF (K.GT.NBINS)THEN
        CALL HBUG('Not enough bins in DIFFS to hold result','HDIFFB',
     +      IDR)
        WRITE(DUMPDV, FMT=910) K
        ERRORS = .TRUE.
      ENDIF
C
C
C----------------------------------------------------------------------
C   Jump ship if any errors have been detected
C----------------------------------------------------------------------
C
      IF (ERRORS) GO TO 999
C
C======================================================================
C   Section 2: Continue with the second stage of initialization.
C======================================================================
C
C
C----------------------------------------------------------------------
C   Find the sum of the weights for each histogram
C----------------------------------------------------------------------
C
C   Clear sums of number of bins with zero contents
C
      NUMZSR = 0
      NUMZSD = 0
C
C
      SUMR = 0.
C                                     ! Clear sums in "good" regions
      SUMD = 0.
      TOTALR =0.
C                                     ! Clear total sums
      TOTALD =0.
      L=0
C                                     ! Initialize L for 1-D
      IF (TWODIM) L=1
C                                     ! Add for Y overflow only on 2-D
C
      DO 41 J=0, NCYR+L
C                                     ! including O/U flows
        DO 40 I=0, NCXR+1
C                                     ! Sum each bin contents
          R = HGCONT(IDR,I,J,1)
C                                     ! Save content
          IF(R.LT.0.)NGCONR=.TRUE.
          TOTALR=TOTALR+R
C                                     ! Add to total sum
          IF(I.GE.BEGINI .AND. I.LE.ENDI .AND. J.GE.BEGINJ .AND. J.LE.
     +        ENDJ) THEN
            SUMR=SUMR+R
C                                     ! Add to sum if in requested region
C
C
            IF (R.EQ.0.)  NUMZSR = NUMZSR + 1
C
C                                     ! Add to number of zeroes
C
          ENDIF
   40   CONTINUE
   41 CONTINUE
C
      DO 46 J=0, NCYD+L
C                                     ! including O/U flows
        DO 45 I=0, NCXD+1
C                                     ! Sum each bin contents
          D = HGCONT(IDD,I,J,1)
C                                     ! Save content
          IF(D.LT.0.)NGCOND=.TRUE.
          TOTALD=TOTALD+D
C                                     ! Add to total sum
          IF(I.GE.BEGINI .AND. I.LE.ENDI .AND. J.GE.BEGINJ .AND. J.LE.
     +        ENDJ) THEN
            SUMD=SUMD+D
C                                     ! Add to sum if in requested region
C
C
            IF (D.EQ.0.)  NUMZSD = NUMZSD + 1
C
C                                     ! Add to number of zeroes
          ENDIF
   45   CONTINUE
   46 CONTINUE
C
C
C----------------------------------------------------------------------
C   If sum of contents is 0, exit
C----------------------------------------------------------------------
C
      IF ((SUMR .EQ. 0).AND.(.NOT.NGCONR)) THEN
        CALL HBUG('Sum of hitsogram contents is zero!', 'HDIFFB', IDR)
        ERRORS = .TRUE.
      ENDIF
      IF ((SUMD .EQ. 0).AND.(.NOT.NGCOND)) THEN
        CALL HBUG('Sum of hitsogram contents is zero!', 'HDIFFB', IDD)
        ERRORS = .TRUE.
      ENDIF
C
C
C----------------------------------------------------------------------
C   Find out if reference histogram is weighted.
C   Note: We flag as weighted if the number of entries is not equal
C         any of the below:
C
C         - Sum of contents (filled via HFILL)
C         - NBINS (filled via HPAK)
C         - NBINS - #zeroes (filled via HPAK, but had zeores in filling
C           array. In this case HNOENT returns number of non-zero bins)
C
C
C----------------------------------------------------------------------
C
      CALL HNOENT(IDR, NENTR)
      WEIGHR = ((REAL(NENTR) .NE. TOTALR).AND.(NENTR .NE. NBINS)
     +  .AND.(NENTR.NE.(NBINS-NUMZSR)) .AND. .NOT.PROFIL)
      IF (WEIGHR) THEN
C
        IF ((OPTS(3)+OPTS(4)+OPTS(9)+OPTS(10)+
     +      OPTS(11)+OPTS(12)).NE.0) THEN
          CALL HBUG('U/O/R/L/T/B Option with weighted events',
     +        'HDIFFB', IDR)
          ERRORS = .TRUE.
        ENDIF
C
        IF (TWODIM) THEN
          CALL HBUG(
     + 'Weighted or saturated 2-D histogram, results are unreliable!'
     + ,'HDIFFB',IDR)
        ELSEIF (.NOT.ERBARR) THEN
          CALL HBUG('Weighted events and no HBARX','HDIFFB',IDR)
        ENDIF
      ENDIF
C
C
C----------------------------------------------------------------------
C   Check if data histogram is weighted.
C----------------------------------------------------------------------
C
      CALL HNOENT(IDD, NENTD)
      WEIGHD = ((REAL(NENTD).NE.TOTALD).AND.(NENTD.NE.NBINS)
     +  .AND.(NENTD.NE.(NBINS-NUMZSD)).AND.(.NOT.PROFIL))
      IF (WEIGHD) THEN
C
        IF ((OPTS(3)+OPTS(4)+OPTS(9)+OPTS(10)+
     +      OPTS(11)+OPTS(12)).NE.0) THEN
          CALL HBUG('U/O/R/L/T/B Option with weighted events',
     +        'HDIFFB', IDD)
          ERRORS = .TRUE.
        ENDIF
C
        IF (TWODIM) THEN
          CALL HBUG(
     +'Weighted or saturated 2-D histogram, results are unreliable!'
     +,'HDIFFB',IDD)
        ELSEIF (.NOT.ERBARD) THEN
          CALL HBUG('Weighted events and no HBARX','HDIFFB',IDD)
        ENDIF
      ENDIF
C
C
C----------------------------------------------------------------------
C   Beam up if any errors have occured
C----------------------------------------------------------------------
C
      IF (ERRORS) GO TO 999
C
C
C
C----------------------------------------------------------------------
C  Find scaling factor for overall normalization difference between R+D
C----------------------------------------------------------------------
C
C
      IF (OPTS(NORMD).EQ.1) THEN
        LAMBDA = 1
C                                 ! Turn off scaling for N option
      ELSE
        LAMBDA = 1
        IF (SUMR.NE.0) LAMBDA = SUMD/SUMR
C                                 ! Scale to contents in desired region
      ENDIF
C
      IF (OPTS(DEBUG) .EQ. 1) THEN
C                                 ! Tell user weighted status for Debug
        WRITE(DUMPDV, FMT=220) WEIGHR,WEIGHD
        WRITE(DUMPDV, FMT=240) LAMBDA, SUMR, SUMD
C                                 ! Debugging dump
      ENDIF
 
C---------------------------------------------------------------------
C     Put LOG(bigp) - 1. into /HDBCOM/ used in routines to
C     avoid overflow. We use -1. because usually we will get a log
C     approxitmation to test for possible overflow before we
C     actually calculate the number.
C---------------------------------------------------------------------
C
      LNBIGP = LOG(ABS(REAL(BIGP))) - 1.
C
C======================================================================
C  Formats
C======================================================================
C
C
C
C----------------------------------------------------------------------
C     Initialization messages
C----------------------------------------------------------------------
C
  200 FORMAT('1','HDIFFB debugging dump is now on for histograms R=',
     +    I10,' + D=',I10,'.')
  210 FORMAT(1X,'Option string: ',A,'  TOL=',E10.4)
  220 FORMAT(1X,'Ref, Data histograms are of weighted type? ',2L2)
  240 FORMAT(1X,'Ratio',E10.4,' REF Cont',E10.4,' DAT Cont',E10.4)
C
C
C----------------------------------------------------------------------
C     Error message for not enough room in input arrays
C----------------------------------------------------------------------
C
  910 FORMAT(1X,'This histograms requires',I6,' bins for the result.')
C
C
  999 RETURN
      END
