SUBROUTINE VPHA * * Module Number: 15.3.1 * * Module Name: phav * * Keyphrase: * ---------- * Analyze count rate as a function of the discriminator setting and determine * the best discriminator setting. * * Description: * ------------ * Fit the count rate as a function of the discriminator setting. * The measured count rate will be the INTEGRATED pulse height distribution * (IPHD) which is assumed to have the form: * * IPHD = 0.5 * P * SQRT(2*pi) * [1 - erf((x-Q)/(R*SQRT(2)))] + * S * (exp(-x/T)) + U * * where erf is the error function. This is from the assumption that the * pulse height distribution (PHD) has this form: * * PHD = P/R * exp-((x-Q)**2/2*R**2) + S/T * exp(-x/T) * * where P = coeff(1), Q = coeff(2)..., U = coeff(6). * A non-linear least square fitting algorithm is used to determine the fitting * function. * for the assumed pulse height distribution function (PHD), i. e. a combination * of gaussian and exponential, find out the discriminator setting at which: * (1) the setting is 50% of the setting of the gaussian peak, (scheme = * 'HALF_PEAK', = default) * (2) the setting is three HWHM lower than the setting of the gaussian peak, * (scheme = '3HWHM_BELOW_PEAK') * (3) the difference between the area under the gaussian and the area under * the exponential is maximum (scheme = 'MAX_AREA_DIFF'). * Maximum number of input data points is 2000. * Input data are in single precision but are double precision during the * least square calculation. * * FORTRAN Name: VPHA.FOR * * Keywords of Accessed Files and Tables: * -------------------------------------- * Name I/O Description / Comments * * Subroutines Called: * ------------------- * CDBS: * VPHAGT, VPHAIN, VPHALS, VPHAPT * SDAS: * UMSPUT * OTHERS: * None * * History: * -------- * Version Date Author Description * 1 11-15-85 J.-C. Hsu Design and coding * 2 09-20-87 J.-C. Hsu F77 standard * *------------------------------------------------------------------------------- * * === local: * --temperature REAL TEMP(2000), * --lower and upper limits of temperature : TMPMIN, TMPMAX, * --mean temperature : TMPAVE, * --tolerance of sigma-square convergence * --in least square calculation : TOLERN, * --fraction of coefficients modification : FRAC, * --weighting flag : WTFLAG(2), * --areas : EXAREA, GAAREA, TOTAL, DIFF, MAXDIF, * --optimum discriminator setting : OPTIM, DISMIN, DISMAX, * --attributes : HIVOLT * --epoch of the observation DOUBLE PRECISION EPOCH(2000), * --lower and upper limit of epoch : EPMIN, EPMAX, * --mean epoch : EPAVE, * --x and y arrays in double precision : X(2000), Y(2000), : VAR(2), DUMMY(10), * --x and y standard deviation array in double precision : XERR(2000), YERR(2000), * --covariance matrix : MATRIX(6,6), * --coefficients of the fitting equation : COEFF(6), * --chi-squared of the fit : CHISQ * --number of input data points INTEGER NPTS, * --dimension of the covariance matrix : DIM, * --number of iteration in least square : ITER, * --number of variablesand coefficients : NVAR, NCOEFF, * --detector ID : DETID, * --loop indices : I, * --return status : STATUS, STATOK * --scheme of determining best * --discriminator setting CHARACTER*16 SCHEME * --target name CHARACTER*20 TARGET * --VMS file name of output table CHARACTER*128 OFILE * --error message context CHARACTER*130 CONTXT *=========================begin hsp.inc========================================= * --status return code INTEGER OK, ERRNUM(20) * --message destination and priority INTEGER DEST, PRIO DATA OK /0/ DATA ERRNUM /701, 702, 703, 704, 705, 706, 707, 708, 709, 710, : 711, 712, 713, 714, 715, 716, 717, 718, 719, 720/ DATA DEST, PRIO /1, 0/ *=========================end hsp.inc=========================================== *------------------------------------------------------------------------------- * * (maximum) dimension of covariance matrix * DIM = 6 * * input parameters and data * CALL VPHAGT (X, Y, YERR, TEMP, EPOCH, NPTS, : WTFLAG, ITER, TOLERN, FRAC, : HIVOLT, DETID, TARGET, SCHEME, : OFILE, STATUS) IF (STATUS .NE. OK .OR. NPTS .GT. 2000) THEN CONTXT = 'cannot get parameters/data or have more than ' : // '2000 input data points' GO TO 999 END IF * * get initial coefficients * CALL VPHAIN (X, Y, NPTS, : COEFF, NCOEFF, NVAR, STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'error getting initial coefficients' GO TO 999 END IF * CALL VPHALS (X, Y, XERR, YERR, NPTS, WTFLAG, FRAC, : NCOEFF, DIM, NVAR, : COEFF, TOLERN, ITER, : CHISQ, MATRIX, STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'error in least square fitting' GO TO 999 END IF * * find the range of discriminator setting, temperature, and epoch, and * calculate their average * TMPAVE = 0. EPAVE = 0.D0 DISMIN = X(1) DISMAX = X(1) TMPMIN = TEMP(1) TMPMAX = TEMP(1) EPMIN = EPOCH(1) EPMAX = EPOCH(1) * DO 10 I = 1, NPTS TMPAVE = TMPAVE + TEMP(I) EPAVE = EPAVE + EPOCH(I) * IF (X(I) .LT. DISMIN) THEN DISMIN = X(I) ELSE IF (X(I) .GT. DISMAX) THEN DISMAX = X(I) END IF * IF (TEMP(I) .LT. TMPMIN) THEN TMPMIN = TEMP(I) ELSE IF (TEMP(I) .GT. TMPMAX) THEN TMPMAX = TEMP(I) END IF * IF (EPOCH(I) .LT. EPMIN) THEN EPMIN = EPOCH(I) ELSE IF (EPOCH(I) .GT. EPMAX) THEN EPMAX = EPOCH(I) END IF 10 CONTINUE * TMPAVE = TMPAVE / REAL(NPTS) EPAVE = EPAVE / DBLE(NPTS) * * find optimum discriminator threshold setting * MAXDIF = 0 OPTIM = 0 * * case 1: optimum setting = half way between 0 and the setting where the * Gaussian peak is * IF (SCHEME(1:9) .EQ. 'HALF_PEAK') THEN OPTIM = 0.5 * COEFF(2) * * case 2: optimum setting = three Gaussian HWHM less than the setting where * the Gaussian is * ELSE IF (SCHEME(1:16) .EQ. '3HWHM_BELOW_PEAK') THEN OPTIM = COEFF(2) - 3. * 1.17741 * COEFF(3) * * case 3: optimum setting is where the difference between the area under * the Gaussian component and the area under the exponential * component is maximum * ELSE IF (SCHEME(1:13) .EQ. 'MAX_AREA_DIFF') THEN * * determine the area difference of each discriminator setting * DO 20 I = 1, NPTS * * exclude points to the right of the gaussian peak * IF (X(I) .LT. COEFF(2)) THEN EXAREA = COEFF(4) * EXP(-X(I)/COEFF(5)) / : COEFF(5) VAR(1) = X(I) VAR(2) = 0.D0 CALL VPHAFN (VAR, COEFF, : TOTAL, DUMMY, DUMMY) GAAREA = TOTAL - EXAREA - COEFF(6) DIFF = GAAREA - EXAREA * * find maximum area difference * IF (DIFF .GT. MAXDIF) THEN OPTIM = X(I) MAXDIF = DIFF END IF END IF 20 CONTINUE ELSE STATUS = ERRNUM(1) CONTXT = 'illegal name of discriminator determination ' : // 'scheme' GO TO 999 END IF * * round off optimum setting to the nearest integer * OPTIM = REAL(NINT(OPTIM)) * * check if the determined optimum setting is positive * IF (OPTIM .LE. 0.) THEN CONTXT = 'determined optimum discriminator setting ' : // 'is negative' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) END IF * * write result to output table * CALL VPHAPT (OFILE, SCHEME, OPTIM, NPTS, WTFLAG, COEFF, CHISQ, : ITER, TOLERN, FRAC, HIVOLT, DETID, TARGET, : TMPMIN, TMPMAX, TMPAVE, EPMIN, EPMAX, EPAVE, : STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'cannot write to output table' GO TO 999 END IF * STATUS = OK GO TO 1000 * * write error message * 999 CALL UMSPUT ('VPHA: ' // CONTXT, DEST, PRIO, STATOK) * 1000 RETURN END