SUBROUTINE ZPHD(DATA,N,TH,CENTER,SIGMA,NITER,BOX, * INTPHD,DIFPHD, * FITPHD,CENT,SIG,PEAK,A0,A1,STATUS) C C Module Number: 13.4.1 C C Module Name: ZPHD 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 linear 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: ZPHD C C Keywords of accessed files and tables: C ------------------------------------- C none C Subroutines called: C ------------------- C SDAS C USMPUT C OTHERS: C MRQMIN - NUMERICAL RECIPE FIT ROUTINE C C HISTORY: C ------- C version date Author description C 1 Jan, 87 D. Lindler Designed and coded C 2 MAR 88 D. Lindler New I/O and standards C-------------------------------------------------------------------------- C C UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87 C INTEGER STDOUT PARAMETER (STDOUT = 1) INTEGER STDERR PARAMETER (STDERR = 2) C C input parameters C C -- pha data array REAL DATA(512,N) 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 C -- SIZE OF SMOOTHING BOX. INTEGER BOX 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) 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(200) C -- FIT PARAMETERS DOUBLE PRECISION PAR(5) C -- NUMBER OF LEVELS IN DIFPHD INTEGER N1 C -- SUBSCRIPT RANGE INTEGER I1,I2 C -- ERROR PROCESSING CHARACTER*130 CONTXT C -- SCRATCH SCALAR DOUBLE PRECISION VAL C -- THRESHOLD STEP SIZE DOUBLE PRECISION DELT C -- INDEX IN TH AT CENTER INTEGER ICENT C -- APPROX. PEAK LOCATION DOUBLE PRECISION ACENT DOUBLE PRECISION TMP(150) C C C FIT ROUTINE JUNK C INTEGER LISTA(5) DOUBLE PRECISION COVAR(5,5),ALPHA(5,5),ALAMDA,CHISQ,DERIV(5) C C DATA DECLARATIONS C EXTERNAL ZPHAFN DATA LISTA/1,2,3,4,5/ DATA SIGMAS/200*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 DELT=TH(2)-TH(1) N1=N-1 DO 20 I=1,512 DO 20 J=1,N1 DIFPHD(J,I)=(INTPHD(J,I)-INTPHD(J+1,I))/DELT 20 CONTINUE C C SMOOTH THE DIFFERENTIAL PHD C DO 27 I=1,512 DO 24 J=1,N IDELT1=BOX/2 IDELT2=BOX/2 IF (I.EQ.1)IDELT1=0 IF (I.EQ.N)IDELT2=0 SUM=0.0 NSUM=0 DO 26 K=J-IDELT1,J+IDELT2 SUM=SUM+DIFPHD(K,I) NSUM=NSUM+1 26 CONTINUE IF (NSUM.NE.0)THEN TMP(J)=SUM/NSUM ELSE write (16,*)i,j TMP(J)=0. ENDIF 24 CONTINUE DO 28 J=1,N DIFPHD(J,I)=TMP(J) 28 CONTINUE 27 CONTINUE C C DETERMINE SEARCH AREA FOR APPROXIMATE PEAK LOCATION C DO 25 I=1,N1 IF(TH(I).LT.CENTER)ICENT=I 25 CONTINUE I1=ICENT-25 I2=ICENT+25 IF(I1.LT.2) I1=1 IF(I2.GT.(N1-1)) I2=N1-1 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 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.20) .OR. ( (MAX*0.5).LT.MIN ))GO TO 100 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(1)=MAX PAR(2)=ACENT PAR(3)=SIGMA PAR(4)=0.0 PAR(5)=0.0 C C PERFORM FIT -------------------------------------------------------------- C C INITIALIZE ROUTINE C ALAMDA=-1.0 CALL MRQMIN(TH,DIFPHD(1,I),SIGMAS,N1,PAR,5,LISTA,5, * COVAR,ALPHA,5,CHISQ,ZPHAFN,ALAMDA,STATUS) IF(STATUS.NE.0)GO TO 100 C C ITERATE NITER TIMES C DO 80 ITER=1,NITER CALL MRQMIN(TH,DIFPHD(1,I),SIGMAS,N1,PAR,5,LISTA,5, * COVAR,ALPHA,5,CHISQ,ZPHAFN,ALAMDA,STATUS) IF(STATUS.NE.0)GO TO 100 80 CONTINUE C C CLEAN UP C ALAMDA=0.0 CALL MRQMIN(TH,DIFPHD(1,I),SIGMAS,N1,PAR,5,LISTA,5, * COVAR,ALPHA,5,CHISQ,ZPHAFN,ALAMDA,STATUS) IF(STATUS.NE.0)GO TO 100 C C COMPUTE FITTED PHD C DO 90 J=1,N1 CALL ZPHAFN(TH(J),PAR,FITPHD(J,I),DERIV,5) 90 CONTINUE PEAK(I)=PAR(1) CENT(I)=PAR(2) SIG(I)=PAR(3) A0(I)=PAR(4) A1(I)=PAR(5) 100 CONTINUE STATUS=0 GO TO 1000 999 CALL UMSPUT(CONTXT,STDOUT+STDERR,0,STATUS) STATUS=1 1000 RETURN END