PROGRAM T0FIT C C This program reads in a 2D histogram of time vs wire and fits for C a t0 for each wire. The results are packed into histograms and C written out. C C 20-Jan-1994 (GR) Added the option of fitting a T0 per layer/foil, summing C together all the wires/strips in each layer. C 31-Oct-1994 (GR) Added option not to suppress fitted t0's which are C negative. C 08-Nov-1994 (GR) Don't read in the histograms of those layers we C are not fitting. C Get the histogram offset to the first TDC histogram C from the user. C 13-Nov-1994 (GR) Allow better control of the fit region C 14-Feb-1996 (GR) Some more tweaking of the fit region. C 26-Jun-1998 (GR) Added histogram of ratio of plateau height to "noise" C level C C-- Global symbols $$IMPLICIT $$INCLUDE 'KOFIA$INCLUDE:UTC_GEOM.CMN' $$INCLUDE 'KOFIA$INCLUDE:UTC_QUAL.CMN' C-- PAW stuff INTEGER*4 HBSIZE PARAMETER(HBSIZE=3000000) COMMON/PAWC/HMEMOR(HBSIZE) REAL*4 HMEMOR INTEGER*4 ICYCLE C-- Local variables INTEGER*4 IOFFOUT PARAMETER(IOFFOUT=10000) INTEGER*4 NBUCKET,NWIRE,NSTRIP,NTDC PARAMETER(NBUCKET=250,NWIRE=1152,NSTRIP=768,NTDC=NWIRE+NSTRIP) INTEGER*4 WINMAX PARAMETER(WINMAX=100) CHARACTER*80 INFIL,OUTFIL CHARACTER*4 CNUM INTEGER*4 IOFFINW,IOFFINS INTEGER*4 UNI INTEGER*4 IWIRE,ISTRIP,ITDC INTEGER*4 IBUCKET INTEGER*4 ICTRL INTEGER*4 I INTEGER*4 IDUMMY INTEGER*4 NFILES INTEGER*4 IOFSET INTEGER*4 WWINDOW,SWINDOW,WINDOW INTEGER*4 ISTARTW,ISTARTS,ISTART INTEGER*4 ITLO(18),ITHI(18) INTEGER*4 ILAY,JLAY,ISLAY INTEGER*4 IERR INTEGER*4 NLAYFIT,ILAYFIT(18) INTEGER*4 IFOIL INTEGER*4 ITRIAL,ICONST,MAXTRIAL,ICONST2 REAL*4 TIME(NBUCKET),NHIT(NBUCKET,NTDC),NHITL(NBUCKET,18) REAL*4 SIGMA(NTDC),CHIDOF(NTDC),T0(NTDC),DT0(NTDC) REAL*4 ESIGMA(NTDC),NOENT(NTDC),ARATIO(NTDC) REAL*4 SIGMAW(NWIRE),CHIDOFW(NWIRE),T0W(NWIRE),DT0W(NWIRE) REAL*4 ESIGMAW(NWIRE),NOENTW(NWIRE),ARATIOW(NWIRE) REAL*4 SIGMAS(NSTRIP),CHIDOFS(NSTRIP),T0S(NSTRIP) REAL*4 DT0S(NSTRIP),ESIGMAS(NSTRIP),NOENTS(NSTRIP) REAL*4 ARATIOS(NSTRIP) REAL*4 T0DFLT,CHIPRB,NMIN REAL*4 NSIG1,NSIG2 LOGICAL*4 LWIRE LOGICAL*4 VERBOSE LOGICAL*4 BYLAYER LOGICAL*4 LGOOD LOGICAL*4 NONEGA LOGICAL*4 DOITW(3),DOITS(6) LOGICAL*4 WIDE1,WIDE2 CHARACTER*80 CHTITL INTEGER*4 NX,NY,NWT,LOC REAL*4 XMI,XMA,YMI,YMA EQUIVALENCE (T0W(1),T0(1)) EQUIVALENCE (T0S(1),T0(NWIRE+1)) EQUIVALENCE (DT0W(1),DT0(1)) EQUIVALENCE (DT0S(1),DT0(NWIRE+1)) EQUIVALENCE (SIGMAW(1),SIGMA(1)) EQUIVALENCE (SIGMAS(1),SIGMA(NWIRE+1)) EQUIVALENCE (ESIGMAW(1),ESIGMA(1)) EQUIVALENCE (ESIGMAS(1),ESIGMA(NWIRE+1)) EQUIVALENCE (CHIDOFW(1),CHIDOF(1)) EQUIVALENCE (CHIDOFS(1),CHIDOF(NWIRE+1)) EQUIVALENCE (NOENTW(1),NOENT(1)) EQUIVALENCE (NOENTS(1),NOENT(NWIRE+1)) EQUIVALENCE (ARATIOW(1),ARATIO(1)) EQUIVALENCE (ARATIOS(1),ARATIO(NWIRE+1)) C-- Data statements DATA NMIN/600./ DATA WWINDOW/20/,SWINDOW/40/ DATA ITLO/ 1, 49, 97, 145, & 193, 289, 385, 481, & 577, 721, 865,1009, & 1153,1201,1273,1381,1525,1705/ DATA ITHI/ 48, 96, 144, 192, & 288, 384, 480, 576, & 720, 864,1008,1152, & 1200,1272,1380,1524,1704,1920/ DATA MAXTRIAL/20/ C-- Initialize some parameters WRITE(*,'('' Start of fitting region. Wires and strips '')') ACCEPT *, ISTARTW,ISTARTS WRITE(*,*) ISTARTW,ISTARTS WRITE(*,'('' Length of fitting region. Wires and strips '')') ACCEPT *, WWINDOW,SWINDOW WRITE(*,*) WWINDOW,SWINDOW IF(WWINDOW.GT.WINMAX.OR.SWINDOW.GT.WINMAX) THEN WRITE(*,*) 'Length too long. Increase WINMAX.' STOP ENDIF WRITE(*,'('' Min. no. of entries to make fit '')') ACCEPT *, NMIN WRITE(*,*) NMIN WRITE(*,'('' Type 1 for MINUIT output; 0 otherwise '')') ACCEPT *, IDUMMY IF(IDUMMY.EQ.1) THEN VERBOSE = .TRUE. ELSE VERBOSE = .FALSE. ENDIF WRITE(*,*) IDUMMY IF(VERBOSE) THEN ICTRL = 1000012 ELSE ICTRL = 1000002 ENDIF WRITE(*,'('' Default value for t0 in output file '')') ACCEPT *, T0DFLT WRITE(*,*) T0DFLT WRITE(*,'('' Type 1 to suppress negative t0s; 0 otherwise '')') ACCEPT *, IDUMMY IF(IDUMMY.EQ.1) THEN NONEGA = .TRUE. ELSE NONEGA = .FALSE. ENDIF WRITE(*,'('' Type 1 for t0 fit by layer; 0 otherwise '')') ACCEPT *, IDUMMY IF(IDUMMY.EQ.1) THEN BYLAYER = .TRUE. ELSE BYLAYER = .FALSE. ENDIF WRITE(*,*) IDUMMY WRITE(*,'('' Wire/strip quality file '')') ACCEPT 200, INFIL WRITE(*,201) INFIL WRITE(*,'('' Number of layers to fit (1-18): '')') ACCEPT *, NLAYFIT WRITE(*,*) NLAYFIT WRITE(*,'('' Layer numbers to fit (Wires:1-12, Foils:13-18)'')') ACCEPT *, (ILAYFIT(I),I=1,NLAYFIT) WRITE(*,*) (ILAYFIT(I),I=1,NLAYFIT) DO I=1,3 DOITW(I) = .FALSE. END DO DO I=1,6 DOITS(I) = .FALSE. END DO DO I = 1,NLAYFIT ILAY = ILAYFIT(I) IF(ILAY.GT.12) THEN IFOIL = ILAY-12 DOITS(IFOIL) = .TRUE. ELSE ISLAY = (ILAY-1)/4 + 1 DOITW(ISLAY) = .TRUE. ENDIF END DO C-- Read the quality file CALL READ_QUAL(INFIL,IERR) IF(IERR.NE.0) THEN WRITE(*,*) 'Problem reading wire/strip quality file' STOP ENDIF C-- Initialize HBOOK CALL HLIMIT(HBSIZE) C-- Book output histograms C-- Note that we book some more in the routine FITARRAY C-- Wires CALL HBOOK1(IOFFOUT+1,'T0 for each wire', & NWIRE,0.5,FLOAT(NWIRE)+0.5,0.) CALL HBOOK1(IOFFOUT+2,'Sigma for each wire', & NWIRE,0.5,FLOAT(NWIRE)+0.5,0.) CALL HBOOK1(IOFFOUT+3,'Chisq/dof per wire', & NWIRE,0.5,FLOAT(NWIRE)+0.5,0.) CALL HBOOK2(IOFFOUT+4,'Layer vs Chisq/dof. all wires', & 50,0.,10.,12,0.5,12.5,0.) CALL HBOOK2(IOFFOUT+5,'Layer vs Chisq prob. all wires', & 100,0.,1.,12,0.5,12.5,0.) CALL HBOOK2(IOFFOUT+6,'Layer vs T0 error. all wires', & 100,0.,10.,12,0.5,12.5,0.) CALL HBOOK1(IOFFOUT+7,'No. of entries per wire', & NWIRE,0.5,FLOAT(NWIRE)+0.5,0.) CALL HBOOK2(IOFFOUT+8,'Layer vs Sigma. All wires', & 100,0.,30.,12,0.5,12.5,0.) CALL HBOOK1(IOFFOUT+9,'Plateau ratio vs wire', & NWIRE,0.5,FLOAT(NWIRE)+0.5,0.) CALL HBARX(IOFFOUT+1) CALL HBARX(IOFFOUT+2) CALL HBARX(IOFFOUT+3) C-- Strips CALL HBOOK1(IOFFOUT+11,'T0 for each strip', & NSTRIP,0.5,FLOAT(NSTRIP)+0.5,0.) CALL HBOOK1(IOFFOUT+12,'Sigma for each strip', & NSTRIP,0.5,FLOAT(NSTRIP)+0.5,0.) CALL HBOOK1(IOFFOUT+13,'Chisq/dof per strip', & NSTRIP,0.5,FLOAT(NSTRIP)+0.5,0.) CALL HBOOK2(IOFFOUT+14,'Foil vs Chisq/dof. all strips', & 50,0.,10.,6,0.5,6.5,0.) CALL HBOOK2(IOFFOUT+15,'Foil vs Chisq prob. all strips', & 100,0.,1.,6,0.5,6.5,0.) CALL HBOOK2(IOFFOUT+16,'Foil vs T0 error. all strips', & 100,0.,10.,6,0.5,6.5,0.) CALL HBOOK1(IOFFOUT+17,'No. of entries per strip', & NSTRIP,0.5,FLOAT(NSTRIP)+0.5,0.) CALL HBOOK2(IOFFOUT+18,'Foil vs Sigma. All strips', & 100,0.,30.,6,0.5,6.5,0.) CALL HBOOK1(IOFFOUT+19,'Plateau ratio vs strip', & NSTRIP,0.5,FLOAT(NSTRIP)+0.5,0.) CALL HBARX(IOFFOUT+11) CALL HBARX(IOFFOUT+12) CALL HBARX(IOFFOUT+13) C-- Get the input file name NFILES = 0 7734 CONTINUE WRITE(*,'('' Input histogram file (type none to quit) '')') ACCEPT 200, INFIL 200 FORMAT(A80) WRITE(*,201) INFIL 201 FORMAT(X,A80) IF(INFIL.EQ.'NONE'.OR.INFIL.EQ.'none') GO TO 7735 WRITE(*,'('' Input histogram offset for wires, strips'')') ACCEPT *, IOFFINW,IOFFINS WRITE(*,*) IOFFINW,IOFFINS C-- Open the file NFILES = NFILES+1 UNI = 90 OPEN(UNIT=UNI,FILE=INFIL,FORM='UNFORMATTED', & RECL=1024,ACCESS='DIRECT',STATUS='OLD',READONLY) C-- Read in the histograms CALL HRFILE(UNI,'TIME',' ') IF(NFILES.EQ.1) THEN IOFSET = 0 ELSE IOFSET = 99999 ENDIF IF(DOITW(1)) CALL HRIN(IOFFINW, 999,IOFSET) IF(DOITW(2)) CALL HRIN(IOFFINW+1,999,IOFSET) IF(DOITW(3)) CALL HRIN(IOFFINW+2,999,IOFSET) IF(DOITS(1)) CALL HRIN(IOFFINS ,999,IOFSET) IF(DOITS(2)) CALL HRIN(IOFFINS+1,999,IOFSET) IF(DOITS(3)) CALL HRIN(IOFFINS+2,999,IOFSET) IF(DOITS(4)) CALL HRIN(IOFFINS+3,999,IOFSET) IF(DOITS(5)) CALL HRIN(IOFFINS+4,999,IOFSET) IF(DOITS(6)) CALL HRIN(IOFFINS+5,999,IOFSET) C-- Close input file CALL HREND('TIME') CLOSE(UNIT=UNI) GO TO 7734 C-- Dump the info into arrays 7735 CONTINUE IF(DOITW(1)) CALL HUNPAK(IOFFINW, NHIT(1,1), 'HIST',0) IF(DOITW(2)) CALL HUNPAK(IOFFINW+1,NHIT(1,193),'HIST',0) IF(DOITW(3)) CALL HUNPAK(IOFFINW+2,NHIT(1,577),'HIST',0) IF(DOITS(1)) CALL HUNPAK(IOFFINS ,NHIT(1,1153),'HIST',0) IF(DOITS(2)) CALL HUNPAK(IOFFINS+1,NHIT(1,1201),'HIST',0) IF(DOITS(3)) CALL HUNPAK(IOFFINS+2,NHIT(1,1273),'HIST',0) IF(DOITS(4)) CALL HUNPAK(IOFFINS+3,NHIT(1,1381),'HIST',0) IF(DOITS(5)) CALL HUNPAK(IOFFINS+4,NHIT(1,1525),'HIST',0) IF(DOITS(6)) CALL HUNPAK(IOFFINS+5,NHIT(1,1705),'HIST',0) C-- Get the low edge of the time spectra histogram IF(DOITW(1)) THEN CALL HGIVE(IOFFINW,CHTITL,NX,XMI,XMA,NY,YMI,YMA,NWT,LOC) ELSEIF(DOITW(2)) THEN CALL HGIVE(IOFFINW+1,CHTITL,NX,XMI,XMA,NY,YMI,YMA,NWT,LOC) ELSEIF(DOITW(3)) THEN CALL HGIVE(IOFFINW+2,CHTITL,NX,XMI,XMA,NY,YMI,YMA,NWT,LOC) ELSE WRITE(*,*) 'Unable to get low edge of time spectra' STOP ENDIF C-- Delete the histos from memory IF(DOITW(1)) CALL HDELET(IOFFINW) IF(DOITW(2)) CALL HDELET(IOFFINW+1) IF(DOITW(3)) CALL HDELET(IOFFINW+2) IF(DOITS(1)) CALL HDELET(IOFFINS) IF(DOITS(2)) CALL HDELET(IOFFINS+1) IF(DOITS(3)) CALL HDELET(IOFFINS+2) IF(DOITS(4)) CALL HDELET(IOFFINS+3) IF(DOITS(5)) CALL HDELET(IOFFINS+4) IF(DOITS(6)) CALL HDELET(IOFFINS+5) C-- Get time bin centers DO IBUCKET = 1,NBUCKET TIME(IBUCKET) = XMI+FLOAT(IBUCKET)-0.5 END DO CALL VFILL(T0,NTDC,-1.) CALL VFILL(DT0,NTDC,-1.) CALL VFILL(SIGMA,NTDC,-1.) CALL VFILL(ESIGMA,NTDC,-1.) CALL VFILL(CHIDOF,NTDC,-1.) C-- Loop over layers. Anodes then cathodes DO 110 JLAY = 1,NLAYFIT ILAY = ILAYFIT(JLAY) IF(ILAY.GT.12) THEN IFOIL = ILAY-12 ELSE IFOIL=0 ENDIF C-- Here if we are fitting by layer. IF(BYLAYER) THEN C-- First we have to sum up the entries over all the channels in the layer C-- Ignore the channel if has been flagged as noisy CALL VZERO(NHITL(1,ILAY),NBUCKET) DO ITDC = ITLO(ILAY),ITHI(ILAY) IF(ITDC.LE.NWIRE) THEN LWIRE = .TRUE. IWIRE = ITDC LGOOD = IAND(QWIRE(IWIRE),UTNOIS).EQ.0 ELSE LWIRE = .FALSE. ISTRIP = ITDC-NWIRE LGOOD = IAND(QSTRIP(ISTRIP),UTNOIS).EQ.0 ENDIF IF(LGOOD) THEN DO IBUCKET=1,NBUCKET NHITL(IBUCKET,ILAY) = & NHITL(IBUCKET,ILAY) + NHIT(IBUCKET,ITDC) END DO ENDIF END DO C-- Now make the fit 10 FORMAT(I4) ITRIAL = 1 ICONST = 0 ICONST2 = 0 1001 CONTINUE IF(ILAY.LE.12) THEN WRITE(CNUM,10) ILAY CHTITL = 'T0. Wire layer'//CNUM IOFSET = IOFFOUT+100+ILAY WINDOW = WWINDOW+ICONST-ICONST2 ISTART = ISTARTW+ICONST2 ELSE WRITE(CNUM,10) IFOIL CHTITL = 'T0. Foil'//CNUM IOFSET = IOFFOUT+2100+IFOIL WINDOW = SWINDOW+ICONST-ICONST2 ISTART = ISTARTS+ICONST2 ENDIF CALL FITARRAY(NBUCKET,TIME,NHITL(1,ILAY),ISTART,WINDOW,NMIN, & IOFSET,CHTITL,ICTRL,NOENT(ITLO(ILAY)), & T0(ITLO(ILAY)),DT0(ITLO(ILAY)),SIGMA(ITLO(ILAY)), & ESIGMA(ITLO(ILAY)),CHIDOF(ITLO(ILAY)),CHIPRB, & ARATIO(ITLO(ILAY))) C-- Check to see if the histogram was wide enough to get the full rise. C-- Otherwise, extend the histogram and try fitting again. IF(ILAY.LE.12) THEN NSIG1 = 5.0 NSIG2 = 6.0 ELSE NSIG1 = 4.0 NSIG2 = 5.0 ENDIF IF(T0(ITLO(ILAY))+NSIG1*SIGMA(ITLO(ILAY)).GT. & TIME(ISTART+WINDOW)) THEN ITRIAL = ITRIAL+1 IF(ITRIAL.LE.MAXTRIAL) THEN CALL HDELET(IOFSET) ICONST=ICONST+2 GO TO 1001 ENDIF ENDIF C-- If the histogram is too wide, shrink the range and try fitting again. WIDE1 = .FALSE. WIDE2 = .FALSE. IF(T0(ITLO(ILAY))+NSIG2*SIGMA(ITLO(ILAY)).LT. & TIME(ISTART+WINDOW)) WIDE1 = .TRUE. IF(T0(ITLO(ILAY))-NSIG2*SIGMA(ITLO(ILAY)).GT. & TIME(ISTART)) WIDE2 = .TRUE. IF(WIDE1.OR.WIDE2) THEN ITRIAL = ITRIAL+1 IF(ITRIAL.LE.MAXTRIAL) THEN IF(WIDE1) ICONST=ICONST-2 IF(WIDE2) ICONST2=ICONST2+2 IF(ILAY.LE.12) THEN IF(WWINDOW+ICONST-ICONST2.GE.20) THEN CALL HDELET(IOFSET) GO TO 1001 ENDIF ELSE IF(SWINDOW+ICONST-ICONST2.GE.20) THEN CALL HDELET(IOFSET) GO TO 1001 ENDIF ENDIF ENDIF ENDIF C-- Copy the values for this layer into every channel in the layer DO ITDC = ITLO(ILAY)+1,ITHI(ILAY) NOENT(ITDC) = NOENT(ITLO(ILAY)) T0(ITDC) = T0(ITLO(ILAY)) DT0(ITDC) = DT0(ITLO(ILAY)) SIGMA(ITDC) = SIGMA(ITLO(ILAY)) ESIGMA(ITDC) = ESIGMA(ITLO(ILAY)) CHIDOF(ITDC) = CHIDOF(ITLO(ILAY)) END DO ENDIF C-- Loop over TDC channels in this layer DO 101 ITDC = ITLO(ILAY),ITHI(ILAY) IF(ITDC.LE.NWIRE) THEN LWIRE = .TRUE. IWIRE = ITDC ELSE LWIRE = .FALSE. ISTRIP = ITDC-NWIRE ENDIF C-- Fit the NHIT array for this channel, if we are not doing things by layer. IF(.NOT.BYLAYER) THEN ITRIAL = 1 ICONST = 0 ICONST2 = 0 1002 CONTINUE IF(LWIRE) THEN WRITE(CNUM,10) IWIRE CHTITL = 'T0. Wire '//CNUM IOFSET=IOFFOUT+100+IWIRE LGOOD = IAND(QWIRE(IWIRE),UTNOIS).EQ.0 WINDOW = WWINDOW+ICONST-ICONST2 ISTART = ISTARTW+ICONST2 ELSE WRITE(CNUM,10) ISTRIP CHTITL = 'T0. Strip '//CNUM IOFSET=IOFFOUT+2100+ISTRIP LGOOD = IAND(QSTRIP(ISTRIP),UTNOIS).EQ.0 WINDOW = SWINDOW+ICONST-ICONST2 ISTART = ISTARTS+ICONST2 ENDIF IF(LGOOD) THEN CALL FITARRAY(NBUCKET,TIME,NHIT(1,ITDC),ISTART,WINDOW, & NMIN,IOFSET,CHTITL,ICTRL,NOENT(ITDC), & T0(ITDC),DT0(ITDC),SIGMA(ITDC),ESIGMA(ITDC), & CHIDOF(ITDC),CHIPRB,ARATIO(ITDC)) C-- Check to see if the histogram was wide enough to get the full rise. C-- Otherwise, extend the histogram and try fitting again. IF(LWIRE) THEN NSIG1 = 5.0 NSIG2 = 6.0 ELSE NSIG1 = 4.0 NSIG2 = 5.0 ENDIF IF(T0(ITDC)+NSIG1*SIGMA(ITDC).GT.TIME(ISTART+WINDOW)) THEN ITRIAL = ITRIAL+1 IF(ITRIAL.LE.MAXTRIAL) THEN CALL HDELET(IOFSET) ICONST=ICONST+2 GO TO 1002 ENDIF ENDIF C-- If the histogram is too wide, shrink the range and try fitting again. WIDE1 = .FALSE. WIDE2 = .FALSE. IF(T0(ITDC)+NSIG2*SIGMA(ITDC).LT. & TIME(ISTART+WINDOW)) WIDE1 = .TRUE. IF(T0(ITDC)-NSIG2*SIGMA(ITDC).GT. & TIME(ISTART)) WIDE2 = .TRUE. IF(WIDE1.OR.WIDE2) THEN ITRIAL = ITRIAL+1 IF(ITRIAL.LE.MAXTRIAL) THEN IF(WIDE1) ICONST=ICONST-2 IF(WIDE2) ICONST2=ICONST2+2 IF(LWIRE) THEN IF(WWINDOW+ICONST-ICONST2.GE.20) THEN CALL HDELET(IOFSET) GO TO 1002 ENDIF ELSE IF(SWINDOW+ICONST-ICONST2.GE.20) THEN CALL HDELET(IOFSET) GO TO 1002 ENDIF ENDIF ENDIF ENDIF ENDIF ENDIF IF(LWIRE) THEN CALL HFILL(IOFFOUT+4,CHIDOF(ITDC),FLOAT(ILAY),1.) CALL HFILL(IOFFOUT+5,CHIPRB,FLOAT(ILAY),1.) CALL HFILL(IOFFOUT+6,DT0(ITDC),FLOAT(ILAY),1.) CALL HFILL(IOFFOUT+8,SIGMA(ITDC),FLOAT(ILAY),1.) ELSE CALL HFILL(IOFFOUT+14,CHIDOF(ITDC),FLOAT(IFOIL),1.) CALL HFILL(IOFFOUT+15,CHIPRB,FLOAT(IFOIL),1.) CALL HFILL(IOFFOUT+16,DT0(ITDC),FLOAT(IFOIL),1.) CALL HFILL(IOFFOUT+18,SIGMA(ITDC),FLOAT(IFOIL),1.) ENDIF C-- Check that the numbers can be output without error to the output C-- t0 file IF(ABS(T0(ITDC)).GT.999..OR. & ABS(DT0(ITDC)).GT.999..OR. & ABS(SIGMA(ITDC)).GT.999..OR. & ABS(CHIDOF(ITDC)).GT.999.) THEN T0(ITDC) = -1. DT0(ITDC) = 0. SIGMA(ITDC) = -1. ESIGMA(ITDC) = 0. CHIDOF(ITDC) = -1. CHIPRB = -1. ENDIF 101 CONTINUE 110 CONTINUE C-- Pack the results into output histograms DO ITDC=1,NTDC IF(DT0(ITDC).EQ.-1.) DT0(ITDC) = 0. IF(ESIGMA(ITDC).EQ.-1.) ESIGMA(ITDC) = 0. END DO CALL HPAK(IOFFOUT+1,T0W) CALL HPAKE(IOFFOUT+1,DT0W) CALL HPAK(IOFFOUT+2,SIGMAW) CALL HPAKE(IOFFOUT+2,ESIGMAW) CALL HPAK(IOFFOUT+3,CHIDOFW) CALL HPAK(IOFFOUT+7,NOENTW) CALL HPAK(IOFFOUT+8,ARATIOW) CALL HPAK(IOFFOUT+11,T0S) CALL HPAKE(IOFFOUT+11,DT0S) CALL HPAK(IOFFOUT+12,SIGMAS) CALL HPAKE(IOFFOUT+12,ESIGMAS) CALL HPAK(IOFFOUT+13,CHIDOFS) CALL HPAK(IOFFOUT+17,NOENTS) CALL HPAK(IOFFOUT+18,ARATIOS) C-- Write the t0's to a file C-- For the output file, if the t0 for a channel is -1 at this point, meaning C-- that there was some problem getting the t0, or if the fitted C-- t0 error is large, or if the fitted t0 is negative, then we set the C-- t0 to T0DFLT and the chisquare/dof to -1 DO ITDC=1,NTDC IF(T0(ITDC).EQ.-1..OR.DT0(ITDC).GT.10..OR. & (NONEGA.AND.T0(ITDC).LT.0.)) THEN T0(ITDC) = T0DFLT DT0(ITDC) = 0. SIGMA(ITDC) = -1. CHIDOF(ITDC) = -1. ENDIF END DO WRITE(*,'('' Output t0 file: '')') ACCEPT 200, OUTFIL WRITE(*,201) OUTFIL OPEN(UNIT=66,FILE=OUTFIL,STATUS='NEW') WRITE(66,667) (IWIRE,T0(IWIRE),DT0(IWIRE),SIGMA(IWIRE), & CHIDOF(IWIRE),IWIRE=1,NWIRE) WRITE(66,667) (ISTRIP-NWIRE,T0(ISTRIP),DT0(ISTRIP), & SIGMA(ISTRIP),CHIDOF(ISTRIP),ISTRIP=NWIRE+1,NTDC) 667 FORMAT(X,I6,4F10.5) CLOSE(66) C-- Save the output histograms WRITE(*,'('' Output histogram file name: '')') ACCEPT 200, OUTFIL WRITE(*,201) OUTFIL OPEN(UNIT=66,FILE=OUTFIL,FORM='UNFORMATTED', & RECL=1024,ACCESS='DIRECT',STATUS='NEW') CALL HRFILE(66,'TZERO','N') CALL HROUT(0,ICYCLE,' ') CALL HREND('TZERO') CLOSE(66) STOP END