program eventanal CHARACTER*40 FILES(17)*40, TEXT40*40,TEXT20*20,COMMENT*300 CHARACTER chtop*8 REAL SCOPEDATA(500,8), AMPLITUDE(3),XPOS(3),XFIN REAL CGCOEFF(5,2),time,XRAL,XRAR REAL COEFF(6,3),WIN(2,3),SCALE,AVG,SIGMA,CHISQ,DEV(3) INTEGER NCHANNEL,NWORD,IODATA, IODATA2,IODATA3, BADCOUNT LOGICAL GOOD INTEGER NWPAWC common / pawpoly / coeff common /bad/ badcount common /counter/ j C Paw memory stuff PARAMETER (NWPAWC=1000000) REAL HMEM COMMON/PAWC/HMEM(NWPAWC) CALL HLIMIT(NWPAWC) C Paw memory stuff ends XPOS(1) = 1.595 XPOS(2) = 0. XPOS(3) = -1.595 IODATA=31 IODATA2=32 IODATA3=33 OPEN(IODATA2,FILE="coefficients.txt") OPEN(IODATA3,FILE="ratiofunctions.txt") OPEN(45,FILE="xmea_diff_differr_sig_sigerr.txt") OPEN(46,FILE="xmea_diff_differr_sig_sigerr_forratio.txt") READ(IODATA3,*) CGCOEFF(5,1),CGCOEFF(4,1),CGCOEFF(3,1), + CGCOEFF(2,1),CGCOEFF(1,1) READ(IODATA3,*) CGCOEFF(5,2),CGCOEFF(4,2),CGCOEFF(3,2), + CGCOEFF(2,2),CGCOEFF(1,2) CLOSE(IODATA3) CALL GETDATAFILES(FILES) chtop="scope" CALL HROPEN(33,chtop,"finalplots2.rzhist", "N", 1024,ierr) call hbook1(900, 'Residual vs. Run Number', 17, 0.5, 17.5,0.) call hbook1(901, 'Sigma vs. Run Number', 17, 0.5, 17.5,0.) CALL HBOOK1(902, "INTERPOLATED FUNCTION RR",150,0.,1.5,0.) CALL HBOOK1(903, "INTERPOLATED FUNCTION RL",150,0.,1.5,0.) CALL HBOOK1(904, "Sigma vs. Run Number for Ratio",17,0.5,17.5,0.) CALL HBOOK1(905, "Residual vs. Run Number for Ratio",17,0.5,17.5, + 0.) R1 = 0. DO I=1,150 CALL XRATIO(R1,R1,O1,O2) R1 = R1 + 0.01 CALL HFILL(902,R1,0.,O1) CALL HFILL(903,R1,0.,O2) ENDDO DO I=1,17 OPEN(IODATA,FILE=FILES(I),FORM="unformatted") READ(IODATA) TEXT20 READ(TEXT20,*) NCHANNEL, NWORD, TIME READ(IODATA) COMMENT READ(COMMENT(13:19),*) XMEASURED XMEASURED = XMEASURED-21.983 WRITE(TEXT40,*) "Run",I," Centroid Distribution" CALL HBOOK1(I,TEXT40,80, XMEASURED-1.,XMEASURED+1.,0.) WRITE(TEXT40,*) "Run",I," RL X Distribution" CALL HBOOK1(20+I,TEXT40,80,XMEASURED-1.,XMEASURED+1.,0.) WRITE(TEXT40,*) "Run",I," RR X Distribution" CALL HBOOK1(40+I,TEXT40,80,XMEASURED-1.,XMEASURED+1.,0.) WRITE(TEXT40,*) "Run",I," Left Ratio" CALL HBOOK1(60+I,TEXT40,100, 0.,1.5,0.) WRITE(TEXT40,*) "Run",I," Right Ratio" CALL HBOOK1(80+I,TEXT40,100, 0.,1.5,0.) CALL GETPARAMS(COEFF,WIN,IODATA2) do J=1,500 CALL READEVENT(IODATA,SCOPEDATA) enddo do J=501,1000 CALL READEVENT(IODATA,SCOPEDATA) CALL FILTER(SCOPEDATA, GOOD) IF (GOOD) THEN CALL FITTER(AMPLITUDE,SCOPEDATA,WIN,COEFF) CALL XCALC(XFIN,XPOS,AMPLITUDE) CALL XRAT(AMPLITUDE,XRAL,XRAR,CGCOEFF) CALL HFILL(20+I,XRAL,0.,1.) CALL HFILL(40+I,XRAR,0.,1.) CALL HFILL(60+I,AMPLITUDE(1)/AMPLITUDE(2),0.,1.) CALL HFILL(80+I,AMPLITUDE(3)/AMPLITUDE(2),0.,1.) CALL HFILL(I,XFIN,0.,1.) ENDIF ENDDO CALL HFITGA(I,SCALE,AVG,SIGMA,CHISQ,2,DEV) if (chisq .gt. 1.) then do z=1,3 dev(z) = dev(z)*sqrt(chisq) enddo endif write(45,*) XMEASURED,(AVG-XMEASURED)*1000.,dev(2)*1000. + ,SIGMA*1000.,dev(3)*1000. CALL HFILL(900,REAL(I),0.,AVG-xmeasured) CALL HFILL(901,REAL(I),0.,SIGMA) CALL HFITGA(20+I,SCALE,AVG,SIGMA,CHISQ,2,DEV) IF (I .LE. 4) then write(46,*) XMEASURED,(AVG-XMEASURED)*1000.,dev(2)*1000. + ,SIGMA*1000.,dev(3)*1000. call HFILL(904,REAL(I),0.,SIGMA) call HFILL(905,REAL(I),0.,AVG-XMEASURED) ENDIF CALL HFITGA(40+I,SCALE,AVG,SIGMA,CHISQ,2,DEV) IF (I .ge. 11) then write(46,*) XMEASURED,(AVG-XMEASURED)*1000.,dev(2)*1000. + ,SIGMA*1000.,dev(3)*1000. call HFILL(904,REAL(I),0.,SIGMA) call HFILL(905,REAL(I),0.,AVG-XMEASURED) ENDIF CLOSE(IODATA) ENDDO close(iodata2) close(45) icycle=0 call hrout(0,icycle,"N") call hrend(chtop) RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE GETDATAFILES(FILEN) IMPLICIT NONE CHARACTER*40 FILEN(17) FILEN(1) = '../outfiles/1_11_02_edg18_19.out' FILEN(2) = '../outfiles/1_11_02_strip19_1.out' FILEN(3) = '../outfiles/1_11_02_strip19_2.out' FILEN(4) = '../outfiles/1_11_02_strip19_3.out' FILEN(5) = '../outfiles/1_11_02_strip19_4.out' FILEN(6) = '../outfiles/1_11_02_strip19_5.out' FILEN(7) = '../outfiles/1_11_02_strip19_6.out' FILEN(8) = '../outfiles/1_11_02_strip19_7.out' FILEN(9) = '../outfiles/1_11_02_cntr19.out' FILEN(10) = '../outfiles/1_11_02_strip19_8.out' FILEN(11) = '../outfiles/1_11_02_strip19_9.out' FILEN(12) = '../outfiles/1_11_02_strip19_10.out' FILEN(13) = '../outfiles/1_11_02_strip19_11.out' FILEN(14) = '../outfiles/1_11_02_strip19_12.out' FILEN(15) = '../outfiles/1_11_02_strip19_13.out' FILEN(16) = '../outfiles/1_11_02_strip19_14.out' FILEN(17) = '../outfiles/1_11_02_edg19_20.out' RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE READEVENT(IODATA,SCOPEDATA) IMPLICIT NONE REAL SCOPEDATA(500,8), SCDATA(500) CHARACTER*20 TEXTDATA INTEGER IODATA, EVENT, PAD,I,J READ(IODATA) TEXTDATA READ(TEXTDATA,90) EVENT 90 FORMAT(6X,I7) DO I=1,8 READ(IODATA) TEXTDATA READ(TEXTDATA, 91) PAD 91 FORMAT(4X, I7) READ(IODATA) SCDATA DO J=1,500 SCOPEDATA(J,I) = SCDATA(J) ENDDO ENDDO RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE FILTER(SCOPEDATA,GOOD) IMPLICIT NONE LOGICAL GOOD REAL SCOPEDATA(500,8),MIN,TEMP,DIFF INTEGER POS,SCOPENUM, I,M,Q, BADCOUNT common /bad/ badcount GOOD = .TRUE. MIN = 0. DO I=1,8 DO M=50,150 TEMP = SCOPEDATA(M,I) IF (TEMP .LE. MIN) THEN MIN = TEMP POS = M SCOPENUM = I ENDIF ENDDO ENDDO IF ((MIN .LE. -0.12) .AND. (MIN .GE. -0.35)) THEN DO Q=POS,150 DIFF = (SCOPEDATA(Q,SCOPENUM)-SCOPEDATA(Q+1,SCOPENUM))/MIN IF (ABS(DIFF) .GE. 0.2) THEN GOOD = .FALSE. badcount = badcount +1 ENDIF ENDDO ELSE badcount = badcount+1 GOOD = .FALSE. ENDIF RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE FITTER(AMPLITUDE,SCOPEDATA,WIN,COEFF) IMPLICIT NONE EXTERNAL UFUNC REAL AMPLITUDE(3),SCOPEDATA(500,8),DEV(2),chi,CHI2 REAL COEFF(6,3),win(2,3),XVAL,HMIN,HSUM, BESTCHI2 DOUBLE PRECISION UFUNC,p(1) INTEGER IODATA2,I,j,k,l,lbest CHARACTER TEXT20*20 common /counter/ j common /comp/ p COMMON /pawpoly2/ I p(1) = 1.d0 DO I=4,6 CALL HBOOK1(300,"FIT",INT(WIN(2,I-3)-WIN(1,I-3)), + WIN(1,I-3) + ,WIN(2,I-3),0.) CALL HPAK(300,SCOPEDATA(INT(WIN(1,I-3)),I)) CALL HFITS(300,UFUNC,1,P,CHI,2,DEV) CALL HDELET(300) CHI2=0. XVAL = WIN(1,I-3)+0.5 do K=INT(WIN(1,I-3)),INT(WIN(2,I-3)) XVAL = XVAL+1 CHI2 = CHI2+(SCOPEDATA(K,I)*10000.-UFUNC(XVAL)*10000.)**2 + /ABS(SCOPEDATA(K,I)*10000.) ENDDO CHI2 = CHI2 / (WIN(2,I-3)-WIN(1,I-3)-1) if (chi2 .ge. 2.) then CALL HBOOK1(300,"SMOOTH",INT(WIN(2,I-3)-WIN(1,I-3))+10, + WIN(1,I-3)-10.,WIN(2,I-3),0.) CALL HPAK(300,SCOPEDATA(INT(WIN(1,I-3))-10,I)) CALL HSMOOF(300,1,CHI) IF (WIN(2,I-3)-WIN(1,I-3) .LT. 60) THEN AMPLITUDE(I-3)= HMIN(300) ELSE AMPLITUDE(I-3) = HSUM(300)/(WIN(2,I-3)-WIN(1,I-3)) ENDIF CALL HDELET(300) ELSE CALL LOCALMAXMIN(AMPLITUDE(I-3),WIN(1,I-3),WIN(2,I-3) + ,COEFF,P,I-3) ENDIF ENDDO RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE LOCALMAXMIN(EXT, LSTART, RSTART,COEFF,P,I) IMPLICIT NONE REAL M,L,R, LFUN, RFUN, MFUN, LSTART, RSTART, COEFF(6,3),EXT DOUBLE PRECISION UFUNC,p(1) INTEGER I L=LSTART R=RSTART DO WHILE(ABS(R-L) .GE. 0.01) M=(L+R)/2. LFUN=real(P(1))*(5*COEFF(6,I)*L**4+4*COEFF(5,I)*L**3+ + 3*COEFF(4,I)*L**2 + +2*COEFF(3,I)*L+COEFF(2,I)) RFUN=real(P(1))*(5*COEFF(6,I)*R**4+4*COEFF(5,I)*R**3+ + 3*COEFF(4,I)*R**2 + +2*COEFF(3,I)*R+COEFF(2,I)) MFUN=real(P(1))*(5*COEFF(6,I)*M**4+4*COEFF(5,I)*M**3+ + 3*COEFF(4,I)*M**2 + +2*COEFF(3,I)*M+COEFF(2,I)) IF(RFUN*MFUN .LE. 0) THEN L=M ELSE R=M ENDIF ENDDO EXT = real(UFUNC(m)) RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE XCALC(XMIC,XPOS,AMPLITUDE) IMPLICIT NONE REAL XTEMP,XPOS(3),AMPLITUDE(3),XMIC, XMICRO INTEGER M XTEMP = 0. DO M=1,3 XTEMP = XTEMP + AMPLITUDE(M)*XPOS(M) ENDDO XTEMP = XTEMP / (AMPLITUDE(1)+AMPLITUDE(2)+AMPLITUDE(3)) XMIC = xmicro(xtemp) RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE XRAT(AMP,XL,XR,CG) IMPLICIT NONE REAL AMP(3),X,XR,XL,ERRL,ERRR,CG(5,2) C CALL xratio(AMP(3)/AMP(2),AMP(1)/AMP(2),XR,XL) CALL FINDX(AMP(1)/AMP(2),CG,2,XL) CALL FINDX(AMP(3)/AMP(2),CG,1,XR) C ERRL = (AMP(2)/AMP(1))**2*(AMP(1)*AMP(2)/(AMP(1) C + +AMP(2))) C ERRR = (AMP(2)/AMP(3))**2*(AMP(3)*AMP(2)/(AMP(3) C + +AMP(2))) C X=SOMETHIHNG HERE EVENTUALLY...MAYBE return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE FINDX(RATIO,CG,I,X) REAL RATIO,CG(5,2),X,M,R,L,MFUN,RFUN,LUN INTEGER I IF (I.EQ.1) THEN L=-1.5 R=0. ELSE L=0. R=1.5 ENDIF DO WHILE (ABS(R-L) .GE. 0.001) M=(R+L)/2. LFUN=CG(5,I)*L**4+CG(4,I)*L**3+CG(3,I)*L**2+CG(2,I)*L+ + CG(1,I)-RATIO RFUN=CG(5,I)*R**4+CG(4,I)*R**3+CG(3,I)*R**2+CG(2,I)*R+ + CG(1,I)-RATIO MFUN=CG(5,I)*M**4+CG(4,I)*M**3+CG(3,I)*M**2+CG(2,I)*M+ + CG(1,I)-RATIO IF(MFUN*RFUN .LT. 0.) THEN L = M ELSE R = M ENDIF ENDDO X = M RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC double precision FUNCTION UFUNC(X) REAL coeff(6,3),TMP,x,xval INTEGER I DOUBLE PRECISION p(1) common /pawpoly/ coeff common /comp/ p COMMON /pawpoly2/ I TMP = real(p(1))*(COEFF(6,I-3)*X**5+COEFF(5,I-3)*X**4+ + COEFF(4,I-3)*X**3 + +COEFF(3,I-3)*X**2+COEFF(2,I-3)*X+COEFF(1,I-3)) UFUNC = TMP RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE GETPARAMS(COEFF,WIN,IODATA2) REAL WIN(2,3),COEFF(6,3) INTEGER IODATA2 READ(IODATA2,*) WIN(1, 1),WIN(2, 1),COEFF(6, 1) + ,COEFF(5,1),COEFF(4,1), + COEFF(3,1),COEFF(2,1),COEFF(1,1) READ(IODATA2,*) WIN(1, 2),WIN(2, 2),COEFF(6, 2) + ,COEFF(5,2),COEFF(4,2), + COEFF(3,2),COEFF(2,2),COEFF(1,2) READ(IODATA2,*) WIN(1, 3),WIN(2, 3),COEFF(6, 3) + ,COEFF(5,3),COEFF(4,3), + COEFF(3,3),COEFF(2,3),COEFF(1,3) RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine xratio(rr,rl,rrfun,rlfun) dimension xrr(17),xrl(17),ymicro(17) data xrl/ 0.1923, 0.1897, 0.1876,0.1853, 0.1966, 0.1971, 0.2191, > 0.2308, 0.2444, 0.2648, 0.2918, 0.3386, 0.3942, 0.4762, > 0.5827, 0.7467, 0.9724/ data xrr/0.9787, 0.7669, 0.5795, 0.4610, 0.3655, 0.2828, 0.2440, > 0.2083, 0.1878, 0.1731, 0.1614, 0.1536, 0.1517, 0.1511, > 0.1559, 0.1600, 0.1726/ data ymicro/-0.7980, -0.7030, -0.6030, -0.5030, -0.4030, > -0.3030, -0.2030, -0.1030, 0.0000, 0.0970, > 0.1970, 0.2970, 0.3970, 0.4970, 0.5970, > 0.6970, 0.7970/ rrfun=XINTR4(ymicro,xrr,17,rr) rlfun=XINTR4(ymicro,xrl,17,rl) C write (6,*) rrfun, rr, rlfun, rl return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC function xmicro(centroid) dimension xcg(17),ymicro(17) data xcg/ -0.5778, -0.4705, -0.3537, -0.2671, -0.1725, > -0.0923, -0.0272, 0.0249, 0.0631, 0.1017, > 0.1431, 0.1978, 0.2501, 0.3186, 0.3916, > 0.4908, 0.5948/ data ymicro/-0.7980, -0.7030, -0.6030, -0.5030, -0.4030, > -0.3030, -0.2030, -0.1030, 0.0000, 0.0970, > 0.1970, 0.2970, 0.3970, 0.4970, 0.5970, > 0.6970, 0.7970/ xmicro=XINTR4(ymicro,xcg,17,centroid) * write(*,*)'Centroid =',centroid,' Compute xmicro=',xmicro return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C XINTR4 FUNCTION XINTR4(Y,X,N,X0) C C 4 POINT INTERPOLATION FROM A TABLE WITH N ENTRIES C GIVEN Y=F(X) IN TABULAR FORM, C FOR A GIVEN VALUE OF X0, C THE FUNCTION XINTR4 RETURNS THE VALUE OF Y0 DIMENSION Y(1),X(1),VAL(4),ARG(4) C DO 15 J=1,N IF(X(J).GE.X0)GO TO 16 15 CONTINUE J=N-1 16 CONTINUE IF(X(J).EQ.X0)GO TO 17 IF(J.EQ.N)J=N-1 IF(J.LT.3)J=3 F=Y(J)/100. J3=J-3 DO 70 M=1,4 J3M=J3+M VAL(M)=Y(J3M) ARG(M)=X(J3M) 70 CONTINUE CALL CFI4(Y0,X0,VAL,ARG,F,IER) XINTR4=Y0 RETURN 17 CONTINUE XINTR4=Y(J) RETURN END C CFI4 C 4 POINT CONTINUED FRACTION INTERPOLATION C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE CFI4(Y,X,VAL,ARG,EPS,IER) C C Y=INTERPOLATED VALUE FOR THE ARGUMENT X . C VAL=VALUE ARRAY . C ARG=ARGUMENT ARRAY . C EPS=DESIRED ACCURACY . C OUTPUT ERROR PARAMETER IER = 0 ACCURACY O.K. C = 1 ACCURACY CAN NOT BE TESTED C IN 4TH ORDER INTERPOLATION . C = 2 TWO IDENTICAL ELEMENTS IN THE C ARGUMENT ARRAY . C REAL*4 ARG(4),VAL(4),X,Y,Z,P1,P2,P3,Q1,Q2,Q3,AUX,H IER=1 Y=VAL(1) P2=1. P3=Y Q2=0. Q3=1. DO 16 I=2,4 II=0 P1=P2 P2=P3 Q1=Q2 Q2=Q3 Z=Y JEND=I-1 3 AUX=VAL(I) DO 10 J=1,JEND H=VAL(I)-VAL(J) IF(ABS(H).GT.1.E-6*ABS(VAL(I)))GO TO 9 IF(ARG(I).EQ.ARG(J))GO TO 17 IF(J.LT.JEND)GO TO 8 II=II+1 III=I+II IF(III.GT.4)GO TO 19 VAL(I)=VAL(III) VAL(III)=AUX AUX=ARG(I) ARG(I)=ARG(III) ARG(III)=AUX GO TO 3 8 VAL(I)=1.E36 GO TO 10 9 VAL(I)=(ARG(I)-ARG(J))/H 10 CONTINUE P3=VAL(I)*P2+(X-ARG(I-1))*P1 Q3=VAL(I)*Q2+(X-ARG(I-1))*Q1 IF(Q3.EQ.0.)GO TO 12 Y=P3/Q3 GO TO 13 12 Y=1.E36 13 DELT=ABS(Z-Y) IF(DELT.LE.EPS)GO TO 19 16 CONTINUE RETURN 17 IER=2 RETURN 19 IER=0 RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc