*
* $Id: hdbini.F,v 1.1.1.1 1996/01/16 17:07:53 mclareni Exp $
*
* $Log: hdbini.F,v $
* Revision 1.1.1.1 1996/01/16 17:07:53 mclareni
* First import
*
*
#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