* * $Id: hdbcop.F,v 1.1.1.1 1996/01/16 17:07:53 mclareni Exp $ * * $Log: hdbcop.F,v $ * Revision 1.1.1.1 1996/01/16 17:07:53 mclareni * First import * * #include "hbook/pilot.h" *CMZ : 4.21/01 28/10/93 08.26.07 by R. J. Genik II *-- Author : R. J. Genik II 23/10/92 SUBROUTINE HDBCOP(TOL,NBINS,NBAD,DIFFS) C---------------------------------------------------------------------- C- C- Purpose and Methods: Statistical compatibility of data histogram C- bin by bin with a reference histogram C- (HDIFFB C option) C- C- Inputs : TOL,NBINS C- Outputs : NBAD,DIFFS, If DEBUG option on, various messages C- Controls: None C- C- Created 3-DEC-1990 James T. McKinley, Michigan State University, USA C- C- C- MODIFIED 1-OCT-1992 R. J. Genik II, Michigan State University, USA C- MODIFIED 27-OCT-1993 R. J. Genik II, Michigan State University, USA C- Divide by zero protection increased. See C- comment below. C- C--------------------------------------------------------------------- C Local variable declarations for HDIFFB C option C--------------------------------------------------------------------- C INTEGER I,J,INDEX,ND #if defined(CERNLIB_DOUBLE) DOUBLE PRECISION DTEMPU,DTEMPL,DGAGNC,DGAPNC #endif #if !defined(CERNLIB_DOUBLE) REAL GAGNC,GAPNC #endif REAL ZVAL,MEAND,L,U,HGCONT,PROB,FREQ,GAMDIS C C--------------------------------------------------------------------- C... passed varaibles, input and output C--------------------------------------------------------------------- C INTEGER NBINS,NBAD REAL TOL,DIFFS(NBINS) #include "hbook/hcdifb.inc" C C C NBAD=0 ZVAL=0. DO 110 J = BEGINJ, ENDJ DO 100 I=BEGINI, ENDI R = HGCONT(IDR, I, J, 1) C ! Get value from Ref HG D = HGCONT(IDD, I, J, 1) C ! Value from Dat HG INDEX = I - BEGINI + 1 C ! Compute position in DIFFS IF (TWODIM) INDEX = INDEX + XSIZ*(J - BEGINJ) C C--------------------------------------------------------------------- C default is to pass C--------------------------------------------------------------------- C DIFFS(INDEX) = 1.0 C C--------------------------------------------------------------------- C C Do actual comparisons. NOTE: R = ID1 = IDR = REFERENCE HISTOGRAM C D = ID2 = IDD = DATA HISTOGRAM C C--------------------------------------------------------------------- C Check for negative contents (and fail) C--------------------------------------------------------------------- C IF(R.LT.0.)THEN DIFFS(INDEX) = 0.0 C !absolute fail WRITE(DUMPDV,FMT=900) I,J,IDR GOTO 90 ENDIF IF(D.LT.0.)THEN DIFFS(INDEX) = 0.0 C !absolute fail WRITE(DUMPDV,FMT=900) I,J,IDD GOTO 90 ENDIF C C--------------------------------------------------------------------- C If option Z has been selected and Ref bin = 0, skip bin C--------------------------------------------------------------------- C IF((OPTS(ZEROS).EQ.1) .AND. (R .EQ. 0)) THEN IF (OPTS(DEBUG) .EQ. 1) WRITE(DUMPDV,FMT=400) I,J C ! Debug opt GO TO 90 ENDIF C C--------------------------------------------------------------------- C If R=D=0 then skip the comp. C--------------------------------------------------------------------- C IF((R+D).EQ.0.) THEN IF (OPTS(DEBUG).EQ.1) WRITE(DUMPDV,FMT=410) 'C', I,J GOTO 90 ENDIF C C--------------------------------------------------------------------- C Get expected value for D, =lambda* (=R) C or =lambda*sum(Wr) C--------------------------------------------------------------------- C MEAND = LAMBDA*R C C--------------------------------------------------------------------- C This is the update referred to above C 27 October 1993 - R. J. Genik II C The calculation of ZVAL will be done here, now, after we make C sure no divide by zero error will occur. C The check for negative contents assures that R and D are equal C to or greater than zero. Then, we are assured that the case C R=D=0 has been handled. If D=0.NE.R, the calculations and C functions called properly return the probability that 0 C was seen when we expected a non-zero result. (The case C lambda = 0 is handled in HDBINI) The case R=0.NE.D is handled C below: we issue an absolute fail because the specification of C the routine states that the error bar is taken from the C reference histogram. We do not allow the user to specify an C error bar in this option, that is what the A-option is for. C e.g. the question: I expected zero and got 3 is not answered. C--------------------------------------------------------------------- C C--------------------------------------------------------------------- C Calculate the SIGD value for Poisson C--------------------------------------------------------------------- C C C SIGD=SQRT(MEAND) C IF((SIGD.EQ.0.).AND.(D.NE.R))THEN DIFFS(INDEX)=0.0 C ! Absolute fail IF(OPTS(DEBUG).EQ.1) WRITE(DUMPDV,FMT=470) I GOTO 90 ENDIF C C---------------------------------------------------------------------- C Compute ZVAL (how many sigmas), used to determine C level of accuracy required and in Large Stat case C---------------------------------------------------------------------- C ZVAL=(D-MEAND)/SIGD C C--------------------------------------------------------------------- C C====================================================================== C Large statistics or weighted C====================================================================== C IF((MEAND.GT.1000000).AND.(TOL.GE.0.01) + .OR.(WEIGHD))THEN C ZVAL=ABS(ZVAL) C ! for Large stat, we remove the sign C C---------------------------------------------------------------------- C Here is where the approximation is: C C We assume that Poisson ==> Gaussian C Hence, we have symmetric distribution with C standard deviation of sqrt(mean), or sqrt(sum w**2), C and then use z value as the test statistic. C C C--------------------------------------------------------------------- C Compute the probabilities C--------------------------------------------------------------------- C DIFFS(INDEX) = 2. - 2.*FREQ(ZVAL) C C---------------------------------------------------------------------- C Do a debugging dump if requested C---------------------------------------------------------------------- C IF (OPTS(DEBUG).EQ.1) THEN WRITE(DUMPDV,FMT=600) R, D, MEAND, SIGD WRITE(DUMPDV,FMT=800)I,J,ZVAL**2,DIFFS(INDEX) ENDIF ELSE C C====================================================================== C Small statistics C====================================================================== C ND = INT(D + 0.5) C !avoid roundoff problems C C---------------------------------------------------------------------- C If we have small tol, are far away from the mean, and are within C the range of called function accuracy, do more accurate L, U C calculation. C---------------------------------------------------------------------- C IF (ND.EQ.0) THEN L = EXP(-MEAND) ELSE IF ((TOL.LT.0.0001).AND.(ZVAL.LT.-3.)) THEN #if defined(CERNLIB_DOUBLE) DTEMPL = DGAGNC(REAL(ND)+1.,MEAND) #endif #if !defined(CERNLIB_DOUBLE) DTEMPL = GAGNC(REAL(ND)+1.,MEAND) #endif L = REAL(DTEMPL) ELSE L = PROB(2.*MEAND, (2*ND+ 2)) ENDIF ENDIF C ! Get lower tail probability IF ((TOL.LT.0.0001).AND.(ZVAL.GT.3.).AND. + (MEAND).LT.10000.) THEN C C---------------------------------------------------------------------- C - The above avoids convergence error in DGAPNC C---------------------------------------------------------------------- C #if defined(CERNLIB_DOUBLE) DTEMPU = DGAPNC(REAL(ND),MEAND) #endif #if !defined(CERNLIB_DOUBLE) DTEMPU = GAPNC(REAL(ND),MEAND) #endif U = REAL(DTEMPU) ELSE IF (ND.EQ.0) THEN U = 1. C C---------------------------------------------------------------------- C - The below avoids convergence error in GAMDIS C---------------------------------------------------------------------- C ELSE IF ((ZVAL.GT.2.0).AND.(MEAND.LT.5000.)) THEN U = GAMDIS(MEAND,REAL(ND)) ELSE C C---------------------------------------------------------------------- C - PROB uses 26.4.14 from Abramowitz, et al, for large C values of mean and degrees of freedom. Accurate to C 3 digits for dof.gt.100, and it doesn't generate errors C for our guarenteed range. C---------------------------------------------------------------------- C U = 1. - PROB(2.*MEAND,2*ND) ENDIF ENDIF C C---------------------------------------------------------------------- C Calculate DIFFS. Use smallest tail * 2 C---------------------------------------------------------------------- C DIFFS(INDEX)=2.*MIN(L,U) C C C--------------------------------------------------------------------- C Display debugging dump if requested C--------------------------------------------------------------------- C IF (OPTS(DEBUG).EQ.1) THEN WRITE(DUMPDV, FMT=610) R,D,L,U ENDIF ENDIF 90 CONTINUE C C---------------------------------------------------------------------- C Check if bin is below user's tol level C---------------------------------------------------------------------- C IF(DIFFS(INDEX) .LT. TOL) NBAD=NBAD+1 C C--------------------------------------------------------------------- C Display debugging dump if requested C--------------------------------------------------------------------- C IF (OPTS(DEBUG).EQ.1) THEN WRITE(DUMPDV, FMT=810) I,J,DIFFS(INDEX),NBAD ENDIF 100 CONTINUE 110 CONTINUE C C C--------------------------------------------------------------------- C Formats C--------------------------------------------------------------------- C C Special case indicators and debug dump info... C--------------------------------------------------------------------- 400 FORMAT('0','Reference bin ',I6,',',I6, + '=0, Z opt, so bin passed') 410 FORMAT('0','Ref=Dat=0 with opts ',A,' in bin ',I6,',',I6) 470 FORMAT('0','C opt and Ref. error bar= 0, thus bin ',I6, + I6,' fails.') C C C option data 600 FORMAT('0','REF',E10.4,' DAT',E10.4,' Expect',E10.4,' EBar',E10. + 4) 610 FORMAT('0','REF',E10.4,' DAT',E10.4,' L',E10.4,' U',E10.4) C C Result for each bin 800 FORMAT(1X,'Bin ',I6,',',I6,': CHISQ',E10.4,' Diffs',E10.4) 810 FORMAT(1X,'Bin ',I6,',',I6,': DIFFS',E10.4,' No. Bad',I6) 900 FORMAT(1X,'Negative bin contents for BIN=',I6,',',I6,'ID=',I6, + 5X, 'Negative bin contents not allowed for S and C options') C END