SUBROUTINE YPHD(DATA,N,TH,CENTER,SIGMA,NITER,TMIN,TMAX, * ORDER,FWIDTH,HMIN,INTPHD,DIFPHD, * FITPHD,CENT,SIG,PEAK,A0,A1,A2,RMS,STATUS) C C Module Number: 14.4.1 C C Module Name: YPHD C C Keyphrase: C ---------- C Pulse height distribution fit routine C C Description: C ------------ C This routine computes the differential pulse height distribution C and fits a Gaussian with a polynomial baseline to the PHD. C CENTER and SIGMA are used as initial guesses of the gaussian C center and sigma. An area of plus or minus 25 data points C are searched for a peak location used as the initial guess to C the nonlinear fit routine. C C FORTRAN name: YPHD C C Keywords of accessed files and tables: C ------------------------------------- C none C Subroutines called: C ------------------- C SDAS C USMPUT C CDBS C MRQMIN, YDIFPH C C HISTORY: C ------- C version date Author description C 1 DEC 88 D. Lindler Designed and coded C-------------------------------------------------------------------------- C C UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87 C INTEGER STDOUT PARAMETER (STDOUT = 1) INTEGER STDERR PARAMETER (STDERR = 2) INTEGER USRLOG PARAMETER (USRLOG = 4) C C input parameters C C -- pha data array REAL DATA(512,256) C -- NUMBER OF LEVELS IN DATA INTEGER N C -- THRESHOLD LEVELS DOUBLE PRECISION TH(N) C -- INITIAL GUESS FOR CENTER DOUBLE PRECISION CENTER C -- INITIAL GUESS FOR SIGMA DOUBLE PRECISION SIGMA C -- NUMBER OF ITERATIONS INTEGER NITER DOUBLE PRECISION TMIN,TMAX C -- RANGE TO PERFORM FIT INTEGER ORDER C -- ORDER OF POLY. BASELINE INTEGER FWIDTH C -- DIFF. FILTER WIDTH DOUBLE PRECISION HMIN C -- MIN. ALLOWED PEAK HEIGHT C C OUTPUT PARAMETERS C C INTEGERAL, DIFFERENTIAL, AND FITTED PHD C DOUBLE PRECISION INTPHD(N,512),DIFPHD(N,512),FITPHD(N,512) C C FIT PARAMETERS C DOUBLE PRECISION CENT(512),SIG(512),PEAK(512),A0(512),A1(512), * A2(512),RMS(512) C C STATUS C INTEGER STATUS C C LOCAL VARIABLES C C -- LOOP INDICES INTEGER I,J,ITER C -- DATA MINS/MAXS DOUBLE PRECISION MIN,MAX C -- RAW DATA SIGMAS DOUBLE PRECISION SIGMAS(256) C -- FIT PARAMETERS DOUBLE PRECISION PAR(6) C -- NUMBER OF LEVELS IN DIFPHD INTEGER I1,I2 C -- ERROR PROCESSING CHARACTER*130 CONTXT C -- SCRATCH SCALAR DOUBLE PRECISION VAL C -- INDEX IN TH AT CENTER INTEGER ICENT C -- APPROX. PEAK LOCATION DOUBLE PRECISION ACENT C C FIT ROUTINE JUNK C INTEGER LISTA(6),FIRST,LAST,NFIT,NTERMS DOUBLE PRECISION COVAR(6,6),ALPHA(6,6),ALAMDA,CHISQ,DERIV(6), * SUMSQ C C DATA DECLARATIONS C EXTERNAL YPHAFN DATA LISTA/1,2,3,4,5,6/ DATA SIGMAS/256*1.0/ C-------------------------------------------------------------------- C C COMPUTE INTEGRAL PHD BY TRANSPOSING DATA C DO 10 I=1,512 DO 10 J=1,N INTPHD(J,I)=DATA(I,J) 10 CONTINUE C C COMPUTE DIFFERENTIAL PHD C CALL YDIFPH(N,INTPHD,TH,TMIN,TMAX,FWIDTH,DIFPHD,FIRST,LAST) IF((LAST-FIRST).LT.10)THEN CONTXT='TOO FEW POINTS TO PERFORM FIT' GO TO 999 ENDIF C C DETERMINE SEARCH AREA FOR APPROXIMATE PEAK LOCATION C DO 25 I=1,N IF(TH(I).LT.CENTER)ICENT=I 25 CONTINUE I1=ICENT-25 I2=ICENT+25 IF(I1.LT.FIRST) I2=FIRST IF(I2.GT.LAST) I2=LAST IF(I2.LE.I1)THEN CONTXT='INVALID SPECIFICATION OF CENTER' GO TO 999 ENDIF C C LOOP ON EACH DIODE AND PERFORM FIT C DO 100 I=1,512 C C INITIALIZE OUTPUT FITTED VALUES C CENT(I)=0.0 SIG(I)=0.0 PEAK(I)=0.0 A0(I)=0.0 A1(I)=0.0 A2(I)=0.0 C C FIND INITIAL GUESS BY SEARCHING REGION BETWEEN I1 AND I2 C MIN=1.0E10 MAX=0.0 DO 30 J=I1,I2 VAL=DIFPHD(J,I)+DIFPHD(J+1,I)+DIFPHD(J-1,I) VAL=VAL/3.0 IF(VAL.GT.MAX)THEN MAX=VAL ACENT=TH(J) ENDIF IF(VAL.LT.MIN)MIN=VAL 30 CONTINUE C C CHECK IF MAX IS LARGE ENOUGH C IF (MAX.LT.HMIN) GO TO 95 C C SET UP INITIAL SEARCH PARAMETERS C C PAR(1) = PEAK OF GAUSSIAN C PAR(2) = CENTER OF GAUSSIAN C PAR(3) = SIGMA OF GAUSSIAN C PAR(4) = CONSTANT TERM OF BASELINE C PAR(5) = SLOPE OF BASELINE C PAR(6) = QUADRATIC TERM OF BASELINE PAR(1)=MAX PAR(2)=ACENT PAR(3)=SIGMA PAR(4)=0.0 PAR(5)=0.0 PAR(6)=0.0 C C PERFORM FIT -------------------------------------------------------------- C C INITIALIZE ROUTINE C NFIT=LAST-FIRST+1 NTERMS=4+ORDER ALAMDA=-1.0 CALL MRQMIN(TH(FIRST),DIFPHD(FIRST,I),SIGMAS,NFIT,PAR, * NTERMS,LISTA,NTERMS, * COVAR,ALPHA,NTERMS,CHISQ,YPHAFN,ALAMDA,STATUS) IF(STATUS.NE.0)GO TO 95 C C ITERATE NITER TIMES C DO 80 ITER=1,NITER CALL MRQMIN(TH(FIRST),DIFPHD(FIRST,I),SIGMAS,NFIT, * PAR,NTERMS,LISTA,NTERMS,COVAR, * ALPHA,NTERMS,CHISQ,YPHAFN,ALAMDA,STATUS) IF(STATUS.NE.0)GO TO 95 80 CONTINUE C C CLEAN UP C ALAMDA=0.0 CALL MRQMIN(TH(FIRST),DIFPHD(FIRST,I),SIGMAS,NFIT,PAR, * NTERMS,LISTA,NTERMS, * COVAR,ALPHA,NTERMS,CHISQ,YPHAFN,ALAMDA,STATUS) IF(STATUS.NE.0)GO TO 95 C C COMPUTE FITTED PHD AND RMS C DO 90 J=1,N CALL YPHAFN(TH(J),PAR,FITPHD(J,I),DERIV,NTERMS) 90 CONTINUE SUMSQ=0.0 DO 91 J=I1,I2 SUMSQ=SUMSQ+(FITPHD(J,I)-DIFPHD(J,I))**2 91 CONTINUE RMS(I)=DSQRT(SUMSQ/NFIT) GO TO 98 C C FAILURE C 95 WRITE(CONTXT,99)I-1 99 FORMAT('No fit computed for diode ',I4) CALL UMSPUT(CONTXT,STDOUT,0,STATUS) RMS(I)=0.0 DO 96 J=1,6 96 PAR(J)=0.0 DO 97 J=1,512 97 FITPHD(J,I)=0.0 98 PEAK(I)=PAR(1) CENT(I)=PAR(2) SIG(I)=PAR(3) A0(I)=PAR(4) A1(I)=PAR(5) A2(I)=PAR(6) 100 CONTINUE STATUS=0 GO TO 1000 999 CALL UMSPUT(CONTXT,STDOUT+STDERR,0,STATUS) STATUS=1 1000 RETURN END