From PARTRIDGE@D0BR06.FNAL.GOV Date: Thu, 8 Jan 1998 19:47:53 -0600 (CST) From: PARTRIDGE@D0BR06.FNAL.GOV To: eno@nscpmail.physics.umd.edu, paterno@fnal.gov Cc: PARTRIDGE@D0BR06.FNAL.GOV Subject: Expected significance routine...I will try comparing with Bruce's routine FUNCTION EXPECTED_SIGNIFICANCE(SIGNAL,BACKGROUND) C---------------------------------------------------------------------- C- C- Purpose and Methods : Determine the average statistical significance C- of an experiment given an expected signal and background. Poisson C- statistics are used except when the number of expected background C- becomes large. There are several parameters in this routine, C- with reasonable default values given below. These parameters are: C- C- EPSILON - This parameter is used to truncate the summations of C- Poisson probabilities. It should be small, but significantly C- larger than the relative machine accuracy. The default value C- for this parameter is 1E-10. C- C- MAX_SIGNIFICANCE - The maximum Gaussian significance. This C- parameter is needed because of the limited range of numbers C- that can be represented on a computer. It is suggested that C- this parameter be set close to the Gaussian significance C- corresponding to a probability that is the smallest positive C- number that can be represented by a REAL*8 number. Note that C- the results of this routine are not reliable when the expected C- significance approaches MAX_SIGNIFICANCE. However, in these cases C- the expected significance would presumably not be of much interest. C- The default value for this parameter is 13. C- C- MIN_SIGNIFICANCE - The minimum Gaussian significance. This C- parameter is needed to set a lower bound on the Gaussian C- significance. For experiments where the number of events C- fluctuates below the expected background, the Gaussian C- significance is negative. For the case of an experiment C- fluctuating to no events observed, the Gaussian significance C- is negative infinity. Thus, it is essential to place a lower C- bound on the Gaussian significance for the expected significance C- to be finite. When the expected number of events is small, C- this parameter has a strong influence on the expected C- significance. When the minimum significance is set to zero, C- fluctuations in the number of events below the expected C- background give a null contribution to the expected significance C- and reasonable results are obtained for the cases investigated. C- The default value for this parameter is 0.0. C- C- GAUSSIAN_CUTOFF - A Gaussian approximation for the most probable C- significance is returned rather than the expected significance C- when the expected number of background events is large. The C- The reason for returning the most probable significance C- ( SIGNAL / SQRT ( BACKGROUND ) ) is to protect the routine C- against potentially long calculation times and round-off errors C- when the number of events is large. The default value for this C- parameter is 500. C- C- Inputs : SIGNAL - Expected number of signal events C- BACKGROUND - Expected number of background events C- Outputs : EXPECTED_SIGNIFICANCE - Average statistical significance C- Controls: none C- C- Created 30-DEC-1997 Richard Partridge C- C---------------------------------------------------------------------- IMPLICIT NONE C INTEGER NBACKGROUND ! Number of BG events in an expt INTEGER NTOTAL ! Total number of events in an expt C REAL*8 BACKGROUND ! Expected background REAL*8 BG_PROBABILITY ! Probability for BG giving NBACKGROUND REAL*8 DGAUSN ! Inverse of normal frequency function REAL*8 DLGAMA ! Log of the Gamma function REAL*8 EPSILON ! Truncation parameter REAL*8 EXPECTED_SIGNIFICANCE ! Expected significance REAL*8 FLUCTUATION_PROBABILITY ! Probability for a BG fluctuation REAL*8 GAUSSIAN_CUTOFF ! BG at which Gaussian approx is used REAL*8 MAX_SIGNIFICANCE ! Maximum significance REAL*8 MIN_SIGNIFICANCE ! Minimum significance REAL*8 PROBABILITY ! Probability of seeing NTOTAL evts REAL*8 PROBABILITY_LEFT ! Probability for a larger NTOTAL REAL*8 SIGNAL ! Expected signal REAL*8 SIGNIFICANCE ! Significance of an expt REAL*8 TOTAL ! Expected signal + background C DATA EPSILON / 1.0D-10 / DATA GAUSSIAN_CUTOFF / 500.0D0 / DATA MAX_SIGNIFICANCE / 13.0D0 / DATA MIN_SIGNIFICANCE / 0.0D0 / C C---------------------------------------------------------------------- C C Decide if we should use a Gaussian approximation, S / SQRT(B) C IF ( BACKGROUND .GE. GAUSSIAN_CUTOFF ) THEN EXPECTED_SIGNIFICANCE = SIGNAL / DSQRT ( BACKGROUND ) ELSE C C Poisson calculation. Initialize variables used in loop over C the number of events in an experiment C TOTAL = SIGNAL + BACKGROUND EXPECTED_SIGNIFICANCE = 0.0 PROBABILITY_LEFT = 1.0 NTOTAL = 0 C C Loop over possible number of events in an experiment. Truncate C the loop when the probability of seeing a larger number of C events becomes negligible C DO WHILE ( PROBABILITY_LEFT .GT. EPSILON ) C C Calculate probability for NTOTAL events in an experiment C PROBABILITY = DEXP ( NTOTAL * DLOG ( TOTAL ) - TOTAL - & DLGAMA ( DFLOAT( NTOTAL + 1 ) ) ) PROBABILITY_LEFT = PROBABILITY_LEFT - PROBABILITY C C Check for special case of NTOTAL = 0 and assign minimum significance C IF (NTOTAL.EQ.0) THEN SIGNIFICANCE = MIN_SIGNIFICANCE ELSE C C Normal case. Initialize calculation of the probability that a C background will fluctuate to give at least NTOTAL events C NBACKGROUND = NTOTAL FLUCTUATION_PROBABILITY = 0.0D0 BG_PROBABILITY = 1.0D0 C C Loop over experiments with NTOTAL or more events and calculate the C background fluctuation probability. Terminate the loop when the C probability of additional terms is negligible C DO WHILE & ( BG_PROBABILITY .GT. EPSILON * FLUCTUATION_PROBABILITY ) C C Calculate the probability for the BG fluctuating to NBACKGROUND C BG_PROBABILITY = DEXP ( NBACKGROUND * DLOG ( BACKGROUND ) & - BACKGROUND - DLGAMA ( DFLOAT ( NBACKGROUND + 1 ) ) ) C C Update total fluctuation probability and loop back for the next value C of NBACKGROUND C FLUCTUATION_PROBABILITY = FLUCTUATION_PROBABILITY + & BG_PROBABILITY NBACKGROUND = NBACKGROUND + 1 ENDDO C C Calculate Gaussian significance corresponding to the background C fluctuation probability and ensure that the significance doesn't C exceed specified maximum and minimum values C IF ( FLUCTUATION_PROBABILITY .GE. 1.0D0 ) THEN SIGNIFICANCE = MIN_SIGNIFICANCE ELSEIF ( FLUCTUATION_PROBABILITY .LE. 0.0D0 ) THEN SIGNIFICANCE = MAX_SIGNIFICANCE ELSE SIGNIFICANCE = - DGAUSN ( FLUCTUATION_PROBABILITY ) SIGNIFICANCE = MIN ( SIGNIFICANCE, MAX_SIGNIFICANCE ) SIGNIFICANCE = MAX ( SIGNIFICANCE, MIN_SIGNIFICANCE ) ENDIF ENDIF C C Calculate the contribution from this value of NTOTAL to the expected C significance and loop back for the next value of NTOTAL C EXPECTED_SIGNIFICANCE = EXPECTED_SIGNIFICANCE + & PROBABILITY * SIGNIFICANCE NTOTAL = NTOTAL + 1 ENDDO ENDIF C 999 RETURN C END