SUBROUTINE YPHA C C Module number: 14.4 C C Module name: YPHA C C Keyphrase: C ---------- C FOS discriminator thresholds and pulse height analysis C C Description: C ------------ C Using FOS PHA observationS, an analysis of each channel is done C by fitting a Gaussian with a polynomial baseline to the differential C pulse distribution. The differential pulse height distribution C is computed from the integral pulse height distribution by C using a smoothing differential filter. C C Fortran Name: YPHA.FOR C C Keywords of accessed files and tables: C -------------------------------------- C INPUT input input PHA observation file name C INTPHD output integral phd file name C DIFPHD output differential pdh file name C FITPHD output fitted phd file name C TABLE output Theshold table file name C T0 input first threshold in INFILE C DELT input threshold stepsize in INFILE C CENTER input initial guess for peak centers (in threshold units) C SIGMA input initial guess for sigma of Gaussians fitted. C NITER input Number of least sqaures iterations C TMIN input Minimum threshold to fit C TMAX input Maximum threshold to fit C ORDER input Order of the polynomial baseline C FWIDTH input Differential filter width C HMIN input Min. allowed peak height to attempt fit C C Subroutines called: C ------------------- C SDAS: * uclgs* , umsput * uttinn, utppti, utcdef, utrpt*, utcpt*, uthad*, uttclo, utccre * uttopn, utpgti, utcfnd, utrgt*, utcgt*, uthgt*, uttclo * uimotp, uimxtp, uimctp, uimopn, uimgid, uimclo, uhdgs* uimclo C OTHERS: C yphd C C HISTORY: C -------- C version date Author Description C 1 Dec 88 D. Lindler Designed and coded C------------------------------------------------------------------------------- C C FILE I/O ACCESS MODES C INTEGER RDONLY PARAMETER (RDONLY = 1) C C CODES FOR DATA TYPES C INTEGER TYINT PARAMETER (TYINT = 4) INTEGER TYREAL PARAMETER (TYREAL = 6) INTEGER TYDOUB PARAMETER (TYDOUB = 7) 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 UHDAS HEADER PARM TYPES -- CB, DAO, 5-SEP-87 C INTEGER GENHDR PARAMETER (GENHDR = 0) C C THIS SECTION IS FOR PARAMETERS RELEVANT TO TABLE I/O. C C THESE MAY BE SET BY UTPPTI AND/OR READ BY UTPGTI: C C NUMBER OF ROWS TO ALLOCATE INTEGER TBALLR PARAMETER (TBALLR = 3) C INCREASE ALLOC NUM OF ROWS INTEGER TBWTYP PARAMETER (TBWTYP = 5) C MAXIMUM NUMBER OF COLUMNS INTEGER TBMXCL PARAMETER (TBMXCL = 7) C TYPE = COLUMN-ORDERED TABLE INTEGER TBTYPC PARAMETER (TBTYPC = 12) C C ERROR PROCESSING PARAMETERS C INTEGER STATUS,ISTAT,ISTATS(20) CHARACTER*130 CONTXT C C KEYWORD PARAMETERS C DOUBLE PRECISION T0,DELT,CENTER,SIGMA,TMIN,TMAX,HMIN INTEGER NITER,FWIDTH,ORDER C C OUTPUT TABLE COLUMN NAMES C CHARACTER*64 TABLE CHARACTER*8 COLNAM(7),CUNIT(7),CFORM(7) INTEGER CTYPE(7),IDOUT,COLIDS(7) DOUBLE PRECISION CENT(512),PEAK(512),SIG(512),A0(512),A1(512) DOUBLE PRECISION A2(512),RMS(512) C C INPUT DATA ARRARY FOR 512 DIODES BY GCOUNT C REAL DATA(512,256) INTEGER GCOUNT C C DATA ARRAYS FOR INTEGRAL, DIFFERENTIAL, AND FITTED PULSE HEIGHT DISTRIBUTIONS C DOUBLE PRECISION INTPHD(256,512),DIFPHD(256,512),FITPHD(256,512) C C ARRAY OF INPUT THRESHOLD LEVELS C DOUBLE PRECISION TH(256) C C FILE PARAMETERS C C C INPUT FILE I/O C CHARACTER*130 INPUT C --->TEMPLATE CHARACTER*64 NAME C --->FILE NAME INTEGER IDTEMP C --->TEMPLATE ID INTEGER IDIN C --->FILE IDS INTEGER NAXIS,DTYPE,DIMEN(8) C --->DATA DESCRIPTIONS C C OUTPUT FILES C CHARACTER*64 INTPH,DIFPH,FITPH C C DETECTOR NUMBER C CHARACTER*5 DET C C LOOP INDEX C INTEGER I C C DATA INITIALIZATION C DATA COLNAM/'CENTER','PEAK','SIGMA','A0','A1','A2','RMS'/ DATA CTYPE/7*TYREAL/ DATA CFORM/'F6.1','F6.1','F6.2','F8.1','F8.3',' ','F8.2'/ DATA CUNIT/7*' '/ C C----------------------------------------------------------------------------- C C C GET KEYWORD PARAMETERS C CALL UCLGST('input',INPUT,ISTATS(1)) CALL UCLGST('intphd',INTPH,ISTATS(2)) CALL UCLGST('difphd',DIFPH,ISTATS(3)) CALL UCLGST('fitphd',FITPH,ISTATS(4)) CALL UCLGST('table',TABLE,ISTATS(5)) CALL UCLGSD('t0',T0,ISTATS(6)) CALL UCLGSD('delt',DELT,ISTATS(7)) CALL UCLGSD('center',CENTER,ISTATS(8)) CALL UCLGSD('sigma',SIGMA,ISTATS(9)) CALL UCLGSI('niter',NITER,ISTATS(10)) CALL UCLGSD('tmin',TMIN,ISTATS(11)) CALL UCLGSD('tmax',TMAX,ISTATS(12)) CALL UCLGSD('hmin',HMIN,ISTATS(13)) CALL UCLGSI('order',ORDER,ISTATS(14)) CALL UCLGSI('fwidth',FWIDTH,ISTATS(15)) DO 5 I=1,15 IF(ISTATS(I).NE.0)THEN CONTXT='Error getting CL parameter' GO TO 999 ENDIF 5 CONTINUE C--------------------------------------------------------------------------- C READ INPUT OBSERVATIONS AND EXTRACTED ENG. DATA FILES C GCOUNT=0 C C OPEN TEMPLATE C CALL UIMOTP(INPUT,IDTEMP,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error opening input filename template '//INPUT GO TO 999 ENDIF C C GET NEXT FILE NAME C 20 CALL UIMXTP(IDTEMP,NAME,ISTAT) IF(ISTAT.LT.0)GO TO 100 IF(ISTAT.NE.0)THEN CONTXT='Error getting filename from template '//INPUT GO TO 999 ENDIF C C OPEN INPUT FILE C CALL UIMOPN(NAME,RDONLY,IDIN,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error opening input file '//NAME GO TO 999 ENDIF C C READ IMAGE INFO C CALL UIMGID(IDIN,DTYPE,NAXIS,DIMEN,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error reading input file '//NAME GO TO 999 ENDIF C C CHECK FOR VALID DATA C IF ( DIMEN(1) .NE. 512 ) THEN CONTXT='Input data must be 512 points' GO TO 999 ELSE IF ( (NAXIS .EQ. 2) .AND. (DIMEN(2) .NE. 1) ) THEN CONTXT='2-D input image must be 512x1' GOTO 999 ELSE IF (NAXIS .GT. 2) THEN CONTXT='Only 1- or 2-D images supported' GOTO 999 ENDIF C C READ DATA C GCOUNT=GCOUNT+1 IF(GCOUNT.GT.256)THEN CONTXT='MAX of 256 input data vectors allowed' GO TO 999 ENDIF CALL UIGL1R(IDIN,DATA(1,GCOUNT),ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error reading input file '//NAME GO TO 999 ENDIF C C GET DETECTOR NUMBER IF IT IS FIRST GROUP C IF(GCOUNT.EQ.1)THEN CALL UHDGST(IDIN,'DETECTOR',DET,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error getting DETECTOR number from '//NAME GO TO 999 ENDIF ENDIF C C CLOSE IMAGE C CALL UIMCLO(IDIN,ISTAT) C C GO GET NEXT IMAGE C GO TO 20 100 CONTINUE C C DONE READING INPUT DATA ------------------------------------------------ C WRITE(CONTXT,99)GCOUNT 99 FORMAT(I5,' Threshold levels read') CALL UMSPUT(CONTXT,STDOUT,0,ISTAT) IF(GCOUNT.LT.15)THEN CONTXT='Insufficient input data files specified' GO TO 999 ENDIF C C WE ARE NOW READY TO DO SOME REAL WORK ----------------------------------- C C C COMPUTE THRESHOLD ARRAY C DO 30 I=1,GCOUNT TH(I)=T0+(I-1)*DELT 30 CONTINUE C C COMPUTE INTEGRAL AND DIFFERENTRIAL PHD C AND PERFORM FITS C CALL YPHD(DATA,GCOUNT,TH,CENTER,SIGMA,NITER,TMIN,TMAX, * ORDER,FWIDTH,HMIN,INTPHD,DIFPHD, * FITPHD,CENT,SIG,PEAK,A0,A1,A2,RMS,STATUS) IF(STATUS.NE.0)THEN CONTXT='ERROR FITTING PULSE HEIGHT DIST.' GO TO 999 ENDIF C C WRITE INTEGRAL PHA --------------------------------------------------------- C DIMEN(1)=GCOUNT DIMEN(2)=512 CALL UIMCRE(INTPH,TYREAL,2,DIMEN,IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error opening output file '//INTPH GO TO 999 ENDIF CALL UIPS2D(IDOUT,1,GCOUNT,1,512,INTPHD,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error writing to output file '//INTPH GO TO 999 ENDIF CALL UHDAST(IDOUT,'DETECTOR',DET,' ',GENHDR,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error writing to output file '//INTPH GO TO 999 ENDIF CALL UIMCLO(IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error closing output file '//INTPH GO TO 999 ENDIF C C WRITE DIFFERENTIAL PHA C CALL UIMCRE(DIFPH,TYREAL,2,DIMEN,IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error opening output file '//DIFPH GO TO 999 ENDIF CALL UIPS2D(IDOUT,1,GCOUNT,1,512,DIFPHD,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error writing to output file '//DIFPH GO TO 999 ENDIF CALL UHDAST(IDOUT,'DETECTOR',DET,' ',GENHDR,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error writing to output file '//DIFPH GO TO 999 ENDIF CALL UIMCLO(IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error closing output file '//DIFPH GO TO 999 ENDIF C C WRITE FITTED PHD C CALL UIMCRE(FITPH,TYREAL,2,DIMEN,IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error opening output file '//FITPH GO TO 999 ENDIF CALL UIPS2D(IDOUT,1,GCOUNT,1,512,FITPHD,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error writing to output file '//FITPH GO TO 999 ENDIF CALL UHDAST(IDOUT,'DETECTOR',DET,' ',GENHDR,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error writing to output file '//FITPH GO TO 999 ENDIF CALL UIMCLO(IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error closing output file '//FITPH GO TO 999 ENDIF C C WRITE OUTPUT TABLE FILE ---------------------------------------------- C C CALL UTTINN(TABLE,IDOUT,ISTATS(1)) CALL UTCDEF(IDOUT,COLNAM,CUNIT,CFORM,CTYPE,7,COLIDS,ISTATS(2)) CALL UTTCRE(IDOUT,ISTATS(3)) DO 200 I=1,3 IF(ISTATS(I).NE.0)THEN CONTXT='Error creating output table '//TABLE GO TO 999 ENDIF 200 CONTINUE C C COPY RESULTS TO TABLE C CALL UTCPTD(IDOUT,COLIDS(1),1,512,CENT,ISTATS(1)) CALL UTCPTD(IDOUT,COLIDS(2),1,512,PEAK,ISTATS(2)) CALL UTCPTD(IDOUT,COLIDS(3),1,512,SIG,ISTATS(3)) CALL UTCPTD(IDOUT,COLIDS(4),1,512,A0,ISTATS(4)) CALL UTCPTD(IDOUT,COLIDS(5),1,512,A1,ISTATS(5)) CALL UTCPTD(IDOUT,COLIDS(6),1,512,A2,ISTATS(6)) CALL UTCPTD(IDOUT,COLIDS(7),1,512,RMS,ISTATS(7)) DO 210 I=1,6 IF(ISTATS(I).NE.0)THEN CONTXT='Error writing to output table' GO TO 999 ENDIF 210 CONTINUE CALL UTHADT(IDOUT,'DETECTOR',DET,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error writing to output table header' GO TO 999 ENDIF CALL UTTCLO(IDOUT,ISTAT) IF(ISTAT.NE.0)THEN CONTXT='Error closing output table '//TABLE GO TO 999 ENDIF C C DONE ------------------------------------------------------------------- C GO TO 1000 999 CALL UMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT) 1000 RETURN END