* SPLFT file SUBROUTINE SPLFT(Y,NDIM) * * Smoothing of poisson distributed data (histograms) * ================================================== * * Usage: CALL SPLFT(Y,NDIM) 4 < NDIM < 257 * * The input data are overwritten by the result. The data are assumed * to be of histogram type with constant weight for all entries, i.e. * "standard deviation of y(i) proportional to sqrt(Y(i))" * In a test sample of distributions the standard deviations of * single data points were reduced by a factor of 0.3, equivalent to * an increase of the number of entries in the histogram by a factor * of more than 10. * The efficiency of smoothing is smaller for histograms with * variable weight. * * Smoothing is done in the following steps: * (1) Transformation of data to stable variance * (2) Smoothing by piecewise fits of orthogonal polynomials, * fourth derivative from result and definition of knot * density [SPOLYN] * (3) Spline fit using knot positions * (4) Backtransformation * (5) Search for peaks and fits of gaussian to peaks, with * improvement of knot density * (6) Spline fit using knot positions repeated * (7) Backtransformation * * Called in SPLFT: * SPOLYN Orthogonal polynomial smoothing * FTORTH fit of orthogonal polynomial * FTSPLN Spline fit with variable knots * SPEAKS Search for peaks in smoothed histograms * FTPEAK Fit of gaussian * * Utility routines: * SPSARP Basic smooth operations * BSPLN4 B-spline calculation * BMDEQU Band-matrix equation * * Author: V. Blobel / October 1994 * REAL Y(*),YS(256) REAL Z(256),RES(16,5),PA(10) * SPLFT - common --------------------------------------------------- PARAMETER (MBSP=100,MKNOT=MBSP+4) COMMON/CMBSP4/NKNOT,T(MKNOT) PARAMETER (MQ=5*MBSP) COMMON/CMQSP4/A(MBSP),B(4),DER(5),FD(256),R(MBSP),Q(MQ) C PARAMETER (MBSP=100,MKNOT=MBSP+4,MQ=5*MBSP) * MBSP = maximum number of spline coefficients C COMMON/CMBSP4/NKNOT,T(MKNOT),A(MBSP),B(4),DER(5),FD(256), C + R(MBSP),Q(MQ) * * T(.) = knot positions (1 ... NKNOT) * A(.) = spline coefficients * B(.) = B-spline values * DER(.) = derivatives * FD(.) = integrated abs of fourth derivative * R(.) = r.h.s. for least squares fit * Q(.) = matrix for least square fit * ------------------------------------------------------------------ * check non-zero content of array and form sum IA=0 IB=0 NOTZ=0 LNEG=0 SUMY=0.0 DO I=1,NDIM IF(Y(I).NE.0.0) THEN IF(Y(I).LT.0.0) LNEG=1 NOTZ=NOTZ+1 IF(IA.EQ.0) IA=I IB=I SUMY=SUMY+ABS(Y(I)) END IF END DO * N= non-zero range of input histogram N=IB-IA+1 IF(IB.EQ.0.OR.N.LT.5.OR.NOTZ.LT.5.OR.N.GT.256) GOTO 100 * * transform non-zero part of array DO I=IA,IB YS(I-IA+1)=Y(I) IF(Y(I).NE.0.0) THEN YI=ABS(Y(I)) Y(I)=SIGN(SQRT(YI)+SQRT(YI+1.0)-1.0,Y(I)) END IF END DO * initial simple smoothing by fits of piecewise orthog. polynomials CALL SPOLYN(Y(IA),N,Z) * loops with spline fits DO LOOP=1,2 * spline fit CALL FTSPLN(Y(IA),N,Z) * inner zero regions of length > 2 are set to zero DO I=IA+1,IB-1 IF(Y(I).EQ.0.0) THEN * three zero values IF(Y(I-1).EQ.0.0.AND.Y(I+1).EQ.0.0) Z(I-IA+1)=0.0 END IF END DO * transform result back to Z(.) SUMZ=0.0 DO I=1,N ZI=ABS(Z(I)) ZI=SIGN(0.25*((ZI+2.0)*ZI/(ZI+1.0))**2,Z(I)) * result must be positive for positive input IF(LNEG.EQ.0.AND.ZI.LT.0.0) ZI=0.0 Z(I)=ZI SUMZ=SUMZ+ABS(Z(I)) END DO IF(LOOP.EQ.1) THEN * search for peaks in smoothed array Z(.) CALL SPEAKS(Z,N,NPEAKS,RES) END IF IF(LOOP.EQ.1.AND.NPEAKS.GT.0) THEN * differentiate knot density FD(.) FDL=0.0 DO I=1,N FDD=FD(I)-FDL FDL=FD(I) FD(I)=FDD END DO * loop over peaks DO J=1,NPEAKS * define fit region ... KA=RES(15,J) KB=RES(16,J) * skip if interval too narrow or peak broad IF(KB-KA.GE.4.AND.RES(5,J).LE.10.0) THEN * ... and define start parameters PA(1)=RES(1,J)-0.5 PA(2)=RES(5,J) PA(3)=RES(2,J)-RES(3,J) PA(4)=RES(12,J) PA(5)=RES(14,J) * fit of gaussian peak to original data CALL FTPEAK(YS ,N,KA,KB,PA) IF(PA(2).GT.0.0) THEN * increase knot density to FD(.) = 1/sigma around peak FDM=0.99/AMAX1(1.00,PA(2)) DO I=1,N FDL=FD(I) X=I ARG =ABS((X-PA(1))/PA(2)) IF(ABS(ARG).LT.3.0.AND.FD(I).LT.FDM) THEN IF(ABS(ARG).LT.2.0) THEN * inside 2 standard deviations FD(I)=AMAX1(FD(I),FDM) ELSE * 2 to 3 standard deviations FD(I)=0.5*(FD(I)+FDM) END IF END IF END DO END IF END IF END DO * sum up knot density CALL SPSARP(FD,N,7) END IF END DO * final normalization to sum of input distribution IF(SUMZ.NE.0.0) SUMZ=1.0/SUMZ DO I=1,N Y(IA+I-1)=Z(I)*SUMY*SUMZ END DO 100 RETURN END SUBROUTINE SPOLYN(Y,N,Z) * * Smooth array Y(.) by fits of orthogonal polynomials. The result is * stored in array Z(.). * * SPLFT - common --------------------------------------------------- PARAMETER (MBSP=100,MKNOT=MBSP+4) COMMON/CMBSP4/NKNOT,T(MKNOT) PARAMETER (MQ=5*MBSP) COMMON/CMQSP4/A(MBSP),B(4),DER(5),FD(256),R(MBSP),Q(MQ) C PARAMETER (MBSP=100,MKNOT=MBSP+4,MQ=5*MBSP) * MBSP = maximum number of spline coefficients C COMMON/CMBSP4/NKNOT,T(MKNOT),A(MBSP),B(4),DER(5),FD(256), C + R(MBSP),Q(MQ) * * T(.) = knot positions (1 ... NKNOT) * A(.) = spline coefficients * B(.) = B-spline values * DER(.) = derivatives * FD(.) = integrated abs of fourth derivative * R(.) = r.h.s. for least squares fit * Q(.) = matrix for least square fit * ------------------------------------------------------------------ REAL Y(*),Z(*) REAL C(5),ZW(13),WW(13),ZW1(13),WW1(13) REAL CUT(5) PARAMETER (EPS=0.01,EPS2=EPS*EPS,EPS4=EPS2*EPS2) * 90 % chi square cuts for 0/2/4/6/8 degrees of freedom DATA CUT/0.0,4.605,7.779,10.645,13.362/ * ... * copy data to Z and ... DO I=1,N FD(I)=0.0 Q(I)=0.0D0 Z(I)=Y(I) END DO * ... smooth using simple 353QH method CALL SPSARP(Z,N,-2) * smooth using fits of orthogonal polynomial M=1 DO I=3,N-2 * maximum value of parameter M MAXM=MIN0(5,N-I-1,I-2) IF(M.GT.MAXM) M=MAXM * fit and derivative around index I using default M value CALL FTORTH(Z(I-M-1),ZW,WW,3+2*M,C) * IF(C(5).GT.CUT(M).AND.M.NE.1) THEN * bad chi square - reduce interval M=M-1 CALL FTORTH(Z(I-M-1),ZW,WW,3+2*M,C) IF(C(5).LE.CUT(M)) THEN * good chi square ELSE * still bad - use nevertheless END IF ELSE IF(C(5).LE.CUT(M)) THEN * good chi square - enlarge interval IF(M.LT.MAXM) THEN * save values DO J=1,2*M+1 ZW1(J)=ZW(J) WW1(J)=WW(J) END DO M=M+1 CALL FTORTH(Z(I-M-1),ZW,WW,3+2*M,C) IF(C(5).LE.CUT(M)) THEN * good chi square ELSE * bad chi square - back to previous M value M=M-1 * restore previous values DO J=1,2*M+1 ZW(J)=ZW1(J) WW(J)=WW1(J) END DO END IF END IF END IF * accumulate values and weights from fit K=0 DO J=I-M-1,I+M+1 K=K+1 Q( J)=Q (J)+WW(K) FD(J)=FD(J)+WW(K)*ZW(K) END DO END DO * determine result as weighted mean DO I=1,N Z(I)=FD(I)/Q(I) END DO * ... DO I=1,N FD(I)=Z(I) END DO * fourth derivative, ... CALL SPSARP(FD,N,9) * ... abs values ... CALL SPSARP(FD,N,8) * .. and hanning (average over three values) CALL SPSARP(FD,N,6) * Transition to knot distances SUM=0.0 DO I=1,N * fourth root H=SQRT(SQRT(FD(I)+EPS4)+3.0*EPS2)-EPS Q(I)=H * transform linear for small arguments, with maximum of 1-eps H=(1.0-EXP(-0.5*H))*(1.0-EPS) FD(I)=H SUM=SUM+H END DO * correct for a minimum of ... knots FMN=AMAX1(3.0,FLOAT(N)/13.0) IF(SUM.LT.FMN) THEN ADD=(FMN-SUM)/FLOAT(N) DO I=1,N FD(I)=FD(I)+ADD END DO END IF * values smoothed by hanning CALL SPSARP(FD,N,6) * sum up CALL SPSARP(FD,N,7) 100 RETURN END SUBROUTINE FTSPLN(Y,N,Z) * * Smooth array Y(.) by spline fits using knot distances, determined * from a previous approximation in array Z(.). The result is stored * in array Z(.). * REAL Y(*),Z(*) * SPLFT - common --------------------------------------------------- PARAMETER (MBSP=100,MKNOT=MBSP+4) COMMON/CMBSP4/NKNOT,T(MKNOT) PARAMETER (MQ=5*MBSP) COMMON/CMQSP4/A(MBSP),B(4),DER(5),FD(256),R(MBSP),Q(MQ) C PARAMETER (MBSP=100,MKNOT=MBSP+4,MQ=5*MBSP) * MBSP = maximum number of spline coefficients C COMMON/CMBSP4/NKNOT,T(MKNOT),A(MBSP),B(4),DER(5),FD(256), C + R(MBSP),Q(MQ) * * T(.) = knot positions (1 ... NKNOT) * A(.) = spline coefficients * B(.) = B-spline values * DER(.) = derivatives * FD(.) = integrated abs of fourth derivative * R(.) = r.h.s. for least squares fit * Q(.) = matrix for least square fit * ------------------------------------------------------------------ * ... * reset optimal NB value and corresponding chi square sum NBOPT=0 CHSM =0.0 * determine number of inner knot positions NSA=AMAX1(1.0,FD(N)+0.5) * three iterations and fourth one with optimal number DO LOOP=1,4 NB=MIN0(57,NSA+LOOP-1,N-4) IF(NB.LT.1) NB=1 * finally use optimal value of NB IF(LOOP.EQ.4) NB=NBOPT * definition of NB+1 different knots ------------------------------ L=0 * define NB+1 different knots DO I=0,NB * fraction of total YF=FLOAT(I)/FLOAT(NB)*FD(N) IF(L.NE.0) THEN IF(YF.LT.FD(L)) L=0 END IF * determine index L 10 IF(YF.GE.FD(L+1).AND.L.LT.N-1) THEN L=L+1 GOTO 10 END IF IF(L.EQ.0) THEN YL=0 ELSE YL=FD(L) END IF T(I+4)=(FLOAT(L)+(YF-YL)/(FD(L+1)-YL))/FLOAT(N) END DO * ... and repeat initial and final knots DO I=1,3 T( I)=T( 4) T(NB+4+I)=T(NB+4) END DO * number of knots ... NKNOT=NB+7 * ... and number of parameters NP=NKNOT-4 IF(LOOP.NE.4) THEN * fit, using integral over bin, in one iteration DO I=1,NP R(I)=0.0 END DO DO I=1,5*NP Q(I)=0.0 END DO * form normal equations INUL=0 DO 60 L=0,2*N CALL BSPLN4(0.5*FLOAT(L)/FLOAT(N),B,I0) AR=R(I0+1)*B(1)+R(I0+2)*B(2)+R(I0+3)*B(3)+R(I0+4)*B(4) IF(L.NE.0) THEN NSUM=I0-INUL+4 C IF(NSUM.GT.5) THEN C WRITE(*,*) NSUM,I0,INUL,L,N C END IF IF(MOD(L,2).EQ.1) THEN * middle point of bin DO I=1,4 DER(I0-INUL+I)=DER(I0-INUL+I)+4.0*B(I) END DO GOTO 60 ELSE * left edge of bin DO I=1,4 DER(I0-INUL+I)=DER(I0-INUL+I)+B(I) END DO END IF DO I=1,5 DER(I)=DER(I)/6.0 END DO * add to normal equations DO I=1,5 R(I+INUL)=R(I+INUL)+DER(I)*Y(L/2) IJ=5*(I+INUL-1) DO J=I,5 IJ=IJ+1 Q(IJ)=Q(IJ)+DER(I)*DER(J) END DO END DO * skip at last edge of last interval IF(L.EQ.2*N) GOTO 60 END IF * left edge of bin INUL=I0 DO I=1,4 DER(I)=B(I) END DO DER(5)=0.0 60 CONTINUE * solve system ... CALL BMDEQU(Q,5,NP,R,DELTA) ELSE * copy best fit DO I=1,NP R(I)=A(I) END DO END IF * reconstruct fitted data with integral over bin CHSQ=0.0 FUN =0.0 DO 70 L=0,2*N CALL BSPLN4(0.5*FLOAT(L)/FLOAT(N),B,I0) AR=R(I0+1)*B(1)+R(I0+2)*B(2)+R(I0+3)*B(3)+R(I0+4)*B(4) IF(L.NE.0) THEN IF(MOD(L,2).EQ.1) THEN * middle point of bin FUN=FUN+4.0*AR GOTO 70 ELSE * left edge of bin FUN=(FUN+AR)/6.0 END IF Z(L/2)=FUN * ... and add to chi square sum CHSQ=CHSQ+(Y(L/2)-FUN)**2 END IF FUN =AR 70 CONTINUE * divide by ndf CHSQ=CHSQ/FLOAT(N-NP) IF(LOOP.LT.4) THEN * compare chi square from fit ... IF(LOOP.EQ.1.OR.CHSQ.LT.CHSM) THEN * ... and save best fit result NBOPT=NB CHSM=CHSQ DO I=1,NP A(I)=R(I) END DO END IF END IF END DO 100 RETURN END SUBROUTINE FTORTH(Y,Z,W,N,C) * * Fit of orthogonal polynomial * * Orthogonal polynomials up to the fourth power are fitted to the * data Y(1)...Y(N), assumed to be equidistant in x (with distance * 1 between adjacent points). * * Result: C(1) = orth. coeff. of first power * C(2) = orth. coeff. of second power * C(3) = orth. coeff. of third power * C(4) = orth. coeff. of fourth power * C(5) = chisquare/ndf after subtraction of fourth order * polynomial from the data * * Restriction: N = 5 or 7 or 9 or 11 or 13 only. * PARAMETER (SNS=1.0) REAL Y(*),Z(*),W(*) REAL C(5),Q(4,20),R(2,5),F(2,5),CUT(5),FC(4) INTEGER INM(5) * 90 % chi cuts for 1 degree of freedom ... DATA CUT1/1.6450/ * ... and chi square for 0/2/4/6/8 degrees of freedom DATA CUT/0.0,4.605,7.779,10.645,13.362/ DATA Q/ 0.632456, 0.534523, 0.316228, 0.119523, + 0.316228,-0.267261,-0.632455,-0.478092, + 0.566947, 0.545545, 0.408248, 0.241747, + 0.377965, 0.000000,-0.408248,-0.564076, + 0.188965,-0.327327,-0.408248, 0.080582, + 0.516398, 0.531816, 0.444950, 0.312894, + 0.387298, 0.132955,-0.222474,-0.469340, + 0.258199,-0.151947,-0.413167,-0.245845, + 0.129099,-0.322888,-0.286039, 0.201145, + 0.476731, 0.512092, 0.458030, 0.354790, + 0.381385, 0.204837,-0.091605,-0.354786, + 0.286039,-0.034139,-0.335887,-0.354787, + 0.190693,-0.204836,-0.351155,-0.059131, + 0.095346,-0.307255,-0.213746, 0.236525, + 0.444750, 0.491690, 0.459934, 0.379458, + 0.370625, 0.245845, 0.000000,-0.252972, + 0.296500, 0.044700,-0.250872,-0.367959, + 0.222375,-0.111747,-0.334496,-0.206977, + 0.148250,-0.223494,-0.292685, 0.042161, + 0.074125,-0.290543,-0.167248, 0.245306/ DATA R/-0.534522, 0.717136,-0.436436, 0.483494, + -0.379868, 0.402291,-0.341394, 0.354787, + -0.312892, 0.321964/ DATA F/0.2635233 ,0.3486083 ,0.6804138E-1,0.4700634E-1, + 0.2648510E-1,0.1303723E-1,0.1272302E-1,0.4927605E-2, + 0.6968684E-2,0.2235864E-2/ DATA INM/0,2,5,9,14/ * ... SN2=1.0 M =N/2-1 IN=INM(M) * determine mean value YQ=0.0 DO I=1,N YQ=YQ+Y(I) END DO YQ=YQ/FLOAT(N) * unit factors FC(1)=1.0 FC(2)=1.0 FC(3)=1.0 FC(4)=1.0 * sum for next coefficients C(1)=0.0 C(2)=(Y(N/2+1)-YQ)*R(1,M) C(3)=0.0 C(4)=(Y(N/2+1)-YQ)*R(2,M) C(5)=(Y(N/2+1)-YQ)**2 DO I=1,N/2 C(1)=C(1)+((Y(N+1-I)-YQ)-(Y(I)-YQ))*Q(1,IN+I) C(2)=C(2)+((Y(N+1-I)-YQ)+(Y(I)-YQ))*Q(2,IN+I) C(3)=C(3)+((Y(N+1-I)-YQ)-(Y(I)-YQ))*Q(3,IN+I) C(4)=C(4)+((Y(N+1-I)-YQ)+(Y(I)-YQ))*Q(4,IN+I) * chi square sum C(5)=C(5)+(Y(N+1-I)-YQ)**2+(Y(I)-YQ)**2 END DO C(5)=C(5)-C(4)**2-C(3)**2-C(2)**2-C(1)**2 IF(M.EQ.1) C(5)=0.0 IF(ABS(C(1)).LT.0.001) C(1)=0.0 IF(ABS(C(2)).LT.0.001) C(2)=0.0 IF(ABS(C(3)).LT.0.001) C(3)=0.0 IF(ABS(C(4)).LT.0.001) C(4)=0.0 * damping of coefficients IF(C(5).LE.CUT(M)) THEN IF(ABS(C(4)).LT.2.588) THEN FC(4)=0.0 C(4) =0.0 IF(ABS(C(3)).LT.CUT1) THEN FC(3)=C(3)**2/(C(3)**2+SN2) C(3) =FC(3)*C(3) IF(ABS(C(2)).LT.CUT1) THEN FC(2)=C(2)**2/(C(2)**2+SN2) C(2) =FC(2)*C(2) END IF END IF END IF END IF * reconstruct function values - center value Z(N/2+1)=YQ+C(2)*R(1,M)+C(4)*R(2,M) WE=1.0/FLOAT(N)+R(1,M)**2+(FC(4)*R(2,M))**2 W(N/2+1)=1.0/WE * ... and left/right values DO I=1,N/2 YE=YQ+FC(2)*C(2)*Q(2,IN+I)+FC(4)*C(4)*Q(4,IN+I) YO= C(1)*Q(1,IN+I)+FC(3)*C(3)*Q(3,IN+I) Z(I )=YE-YO Z(N+1-I)=YE+YO WE=1.0/FLOAT(N)+Q(1,IN+I)**2+(FC(2)*Q(2,IN+I))**2 + +(FC(3)*Q(3,IN+I))**2+(FC(4)*Q(4,IN+I))**2 W(I )=1.0/WE W(N+1-I)=1.0/WE END DO 100 RETURN END SUBROUTINE FTPEAK(HIST,NBINS,IA,IB,PA) * * fit gaussian to bins IA ... IB * Bin content assumed to follow Poisson distribution * PA(1) = mean (in unit of bins) * PA(2) = sigma (in unit of bins) * PA(3) = nr of entries within fit interval * PA(4) = left background level * PA(5) = right background level * PA(6) = mean for correction * PA(7) = sigma for correction * PA(8) = odd correction * PA(9) = even correction REAL HIST(NBINS),FY(256),FMD(256),PA(9) REAL DE(5),DER(5) LOGICAL BGL,BGR REAL Q(5,5),R(5) * ... ARGLIM=ALOG(1.0E-10) XMAX=IB-IA+1 * start assuming background left and right NPA=5 BGL=.TRUE. BGR=.TRUE. CHISUM=0.0 * start with zero modulation MDL=0 01 DO 90 ITER=1,10 * reset normal equations DO J=1,5 R(J)=0.0 DO I=1,5 Q(I,J)=0.0D0 END DO END DO CHLSUM=CHISUM CHISUM=0.0 DELF =0.0 FUN =0.0 FGAU =0.0 FLIN =0.0 II=IA-1 X =II DO 60 I=0,2*(IB-IA+1) IF(MOD(I,2).EQ.1.AND.I.NE.0) THEN * lower bin_edge values from previous index ... II=II+1 FUN=FGAU+FLIN DO J=1,5 DER(J)=DE(J) END DO END IF * gaussian plus linear ARG =-0.5*(X-PA(1))**2/PA(2)**2 IF(ARG.LT.ARGLIM) ARG=ARGLIM FGAU=PA(3)*EXP(ARG) XX=0.5*FLOAT(I) FLIN=(PA(4)*(XMAX-XX)+PA(5)*XX)/XMAX * derivatives of gaussian ... DE(1)=+PA(3)*EXP(ARG)*(X-PA(1))/PA(2)**2 DE(2)= PA(3)*EXP(ARG)*(X-PA(1))**2/PA(2)**3 DE(3)= EXP(ARG) * ... and of linear background DE(4)=(XMAX-XX)/XMAX DE(5)=( XX)/XMAX IF(I.EQ.0) GOTO 60 * add to integrated value IF(MOD(I,2).EQ.1) THEN * bin center FUN=FUN+4.0*(FGAU+FLIN) DO J=1,5 DER(J)=DER(J)+4.0*DE(J) END DO GOTO 60 ELSE * upper bin edge FUN=(FUN+FGAU+FLIN)/6.0 IF(MDL.NE.0) FUN=FMD(II)*FUN DO J=1,5 DER(J)=(DER(J)+DE(J))/6.0 IF(MDL.NE.0) DER(J)=FMD(II)*DER(J) END DO END IF * residuum RES=HIST(II)-FUN FY(II)=FUN * ... and weight for least square WT=1.0 IF(HIST(II).NE.0.0) WT=1.0/ABS(HIST(II)) IF(.NOT.BGL) DER(4)=0.0 IF(.NOT.BGR) DER(5)=0.0 * add to vector and matrix DO J=1,5 R(J)=R(J)+WT*RES*DER(J) DO K=J,5 Q(K-J+1,J)=Q(K-J+1,J)+WT*DER(J)*DER(K) END DO END DO CHISUM=CHISUM+WT*RES*RES 60 X=X+0.5 * solve system of equations CALL BMDEQU(Q,5,5,R,CHSQ) DELP=DELF DO J=1,5 PA(J)=PA(J)+R(J) END DO DELF=CHSQ * tests: peak inside histogram? ... IF(PA(1).LT.0.0.OR.PA(1).GT.FLOAT(NBINS)) GOTO 99 * width and factor positive? IF(PA(2).LE.0.0.OR.PA(3).LE.0.0) GOTO 99 * chisquare per degree of freedom CHINDF=0.0 IF(IB-IA+1-NPA.GT.0) CHINDF=CHISUM/FLOAT(IB-IA+1-NPA) * convergence recognition IF(ITER.EQ.1) GOTO 90 IF(DELF.LT.0.01.AND.AMAX1(DELP,CHLSUM-CHISUM).LT.0.1) GOTO 91 90 CONTINUE GOTO 99 * analyse background 91 L=0 IF(NPA.EQ.5) THEN * background assumed L=4 IF(PA(5).LT.PA(4)) L=5 ELSE IF(NPA.EQ.4) THEN * one-sided background L=4 IF(BGR) L=5 END IF IF(L.EQ.0) GOTO 92 IF(PA(L).GT.0.0) GOTO 92 * assume no background at L_side PA(L)=0.0 NPA=NPA-1 IF(L.EQ.4) BGL=.FALSE. IF(L.EQ.5) BGR=.FALSE. GOTO 01 * analyse peak shape 92 IF(MDL.NE.0) GOTO 100 SUMU2=0.0 SUMG2=0.0 SUMUN=0.0 SUMGN=0.0 DO I=IA,IB X=FLOAT(I)-0.5-PA(1) ARG =-0.5*(X/PA(2))**2 IF(ARG.LT.ARGLIM) ARG=ARGLIM WT=1.0 IF(HIST(I).NE.0.0) WT=1.0/ABS(HIST(I)) WT=WT*FY(I)**2 RAT=HIST(I)/FY(I)-1.0 SUMU2=SUMU2+WT*RAT*EXP(ARG)*X SUMG2=SUMG2+WT*RAT*EXP(ARG)*(X*X-PA(2)*PA(2)) * underflows SUMUN=SUMUN+WT*(EXP(ARG)*X)**2 SUMGN=SUMGN+WT*(EXP(ARG)*(X*X-PA(2)*PA(2)))**2 END DO * test asymmetric modulation FACU=SUMU2/SUMUN TSTU=SUMU2/SQRT(SUMUN) * test asymmetric modulation FACG=SUMG2/SUMGN TSTG=SUMG2/SQRT(SUMGN) IF(FACG.EQ.0.0.AND.FACU.EQ.0.0) GOTO 100 MDL=1 * calculate correction to gaussian PA(6)=PA(1) PA(7)=PA(2) PA(8)=FACU PA(9)=FACG DO I=IA,IB X=FLOAT(I)-0.5-PA(6) ARG =-0.5*(X/PA(7))**2 IF(ARG.LT.ARGLIM) ARG=ARGLIM ADDU=PA(8)*EXP(ARG)*X ADDG=PA(9)*EXP(ARG)*(X*X-PA(7)*PA(7)) FMD(I)=1.0+ADDU+ADDG END DO GOTO 01 99 PA(1)=0.0 PA(2)=0.0 PA(3)=0.0 100 RETURN 107 FORMAT(1X,I5,2F10.2,F9.3,G12.4) 108 FORMAT(1X,I5,3G12.4,' addu addg') END SUBROUTINE SPEAKS(Y,N,NP,R) * ------------------------------------------------------------------ * search for peaks in smoothed histogram * input: * Y(.) = smoothed histogram content * N = number of bins * result: * NP = number of peaks found * R(.) = peak parameters * * peak parameters: j = 1 ... NP * * R( 1,i) peak abscissa * R( 2,j) peak ordinate * R( 3,j) background below peak * R( 4,j) integral * R( 5,j) sigma * R( 6,j) FWHM * R( 7,j) left FWHM abscissa * R( 8,j) left FWHM ordinate * R( 9,j) right FWHM abscissa * R(10,j) right FWHM ordinate * R(11,j) left background abscissa * R(12,j) left background ordinate * R(13,j) right background abscissa * R(14,j) right background ordinate * R(15,j) left background bin * R(16,j) right background bin * * ------------------------------------------------------------------ REAL Y(*),R(16,5) CHARACTER*8 COMT(5) DATA COMT/' X_peak',' height','above bg','integral',' sigma'/ * ... NP=0 SUM=0.0 YMX=0.0 DO I=1,N SUM=SUM+Y(I) YMX=AMAX1(YMX,Y(I)) END DO * ... and analyse smoothed data ==================================== DO 60 I=2,N-1 IF(Y(I).EQ.0.0) GOTO 60 IF(Y(I-1).GT.Y(I)) GOTO 60 IF(Y(I+1).GT.Y(I)) GOTO 60 * first and second derivative DER1=0.5*(Y(I+1)-Y(I-1)) DER2=Y(I+1)-2.0*Y(I)+Y(I-1) IF(DER2.GT.0.0) GOTO 60 XX=0.0 IF(DER2.NE.0.0) XX=-DER1/DER2 IF(ABS(XX).GT.0.5) GOTO 60 * here a maximum (peak) has been found of FM at XM FM=Y(I)+XX*(DER1+0.5*DER2*XX) XM=FLOAT(I)+XX C WRITE(6,102) I,DER1,DER2,XM,FM 102 FORMAT(' ===> DER1 DER2 XM FM ',I4,4G15.5) SND=0. * analyse region below peak ---------------------------------------- RNDMIN=0.0 DER1M =0.0 ISTAT = 0 XL =1.0 J=I 20 IF(J.EQ.2) GOTO 29 J=J-1 DER2=Y(J+1)-2.0*Y(J)+Y(J-1) DER1=0.5*(Y(J+1)-Y(J-1)) IF(ISTAT.EQ.0) THEN * find maximum slope DER1M=AMAX1(DER1M,DER1) * switch to status 1 if second derivative positive IF(DER2.GE.0.0) THEN ISTAT=1 END IF END IF * test slope at 3 sigma IF(ISTAT.EQ.1) THEN * test slope reduction IF(J.GT.2.AND.J.LT.N-1) DER1=0.25*(Y(J+2)-Y(J-2)) IF(DER1.GT.0.2*DER1M) GOTO 20 RN=DER1*(FLOAT(I-J)+XX) RD=FM-Y(J) * constant refers to 3 sigma IF(RD.NE.0.0) THEN RND=RN/RD-0.075 IF(RND.LE.0.0) THEN * more than 2 sigma reached XL=-RND/(SND-RND) GOTO 29 ELSE IF(RNDMIN.GT.0.0) THEN * stop if RND is increasing IF(RND.GT.RNDMIN) GOTO 29 END IF RNDMIN=RND END IF END IF GOTO 20 29 FL=(1.0-XL)*Y(J)+XL*Y(J+1) KL=J XL=FLOAT(J)+XL * analyse region above peak ---------------------------------------- RNDMIN=0.0 DER1M =0.0 ISTAT = 0 XR =1.0 J=I 30 IF(J.EQ.N-1) GOTO 39 J=J+1 DER2=Y(J+1)-2.0*Y(J)+Y(J-1) DER1=0.5*(Y(J+1)-Y(J-1)) IF(ISTAT.EQ.0) THEN * find maximum slope DER1M=AMIN1(DER1M,DER1) IF(DER2.GE.0.0) ISTAT=1 END IF IF(ISTAT.EQ.1) THEN * test slope reduction IF(J.GT.2.AND.J.LT.N-1) DER1=0.25*(Y(J+2)-Y(J-2)) IF(DER1.LT.0.2*DER1M) GOTO 30 RN=-DER1*(FLOAT(J-I)-XX) RD=FM-Y(J) IF(RD.NE.0.0) THEN RND=RN/RD-0.075 IF(RND.LE.0.0) THEN * more than 2 sigma XR=-RND/(SND-RND) GOTO 39 ELSE IF(RNDMIN.GT.0.0) THEN * stop if RND is increasing IF(RND.GT.RNDMIN) GOTO 39 END IF RNDMIN=RND END IF END IF GOTO 30 39 FR=(1.0-XR)*Y(J)+XR*Y(J-1) KR=J XR=FLOAT(J)-XR * test consistency of both sides DO 49 LOOP=1,3 * estimate background level at maximum ----------------------------- FC=0.5*(FL+FR) IF(XR-XL.GT.0.0) FC=(FL*(XR-XM)+FR*(XM-XL))/(XR-XL) FB=FC FS=FM-FB IF(FS.LE.0.0) GOTO 60 IF(FB.LE.0.0.OR.FC.LE.0.0) THEN FL=0.0 FR=0.0 ELSE FL=FL*FB/FC FR=FR*FB/FC END IF * half maximum FD=0.5*(FB+FM) * determine full width at half maximum ----------------------------- DO J=I,KL,-1 IF(Y(J).LE.FD) GOTO 42 END DO GOTO 60 * linear interpolation 42 XLD=FLOAT(J) IF(Y(J).NE.Y(J+1)) XLD=XLD+(FD-Y(J))/(Y(J+1)-Y(J)) DO J=I,KR,+1 IF(Y(J).LE.FD) GOTO 46 END DO GOTO 60 * linear interpolation ... 46 XRD=FLOAT(J-1) IF(Y(J).NE.Y(J-1)) XRD=XRD+(Y(J-1)-FD)/(Y(J-1)-Y(J)) * correct limits XLN=XLD+2.00*(XLD-XM) IF(XL.LT.XLN) THEN XL=XLN KL=XL FL=Y(KL) END IF XRN=XRD+2.00*(XRD-XM) IF(XR.GT.XRN) THEN XR=XRN KR=XR FR=Y(KR) END IF 49 CONTINUE * test symmetry of peak ASYM=AMIN1(2.0+XM-XLD,2.0+XRD-XM)/AMAX1(2.0+XM-XLD,2.0+XRD-XM) IF(ASYM.LT.0.333) GOTO 60 * ... and full width_at-half_max, sigma and integral FWHM= XRD-XLD SIGM=FWHM/2.31 WGHT=2.5066*FS*SIGM BACK=(XR-XL)*FB * reject if signal less than 10 % of (signal +background) ... IF(WGHT/(WGHT+BACK).LT.0.10) GOTO 60 * ... or if signal less than 1 % of total ... IF(WGHT.LT.0.01*SUM) GOTO 60 * ... or if peak less than 2 % of max IF(FS.LT.0.02*YMX) GOTO 60 * store peak parameters in array ----------------------------------- IF(NP.EQ.0) THEN * first peak L =1 NP=1 ELSE * not first peak - determine position of new peak LP=0 DO L=1,NP IF(R(4,L).GT.WGHT) LP=LP+1 END DO IF(LP.LT.NP) THEN DO L=MIN0(5,NP+1),LP+2,-1 DO K=1,16 R(K,L)=R(K,L-1) END DO END DO L=LP+1 NP=MIN0(NP+1,5) ELSE IF(NP.EQ.5) GOTO 60 NP=NP+1 L =NP END IF END IF * store peak parameters under index L R( 1,L)=XM R( 2,L)=FM R( 3,L)=FB R( 4,L)=WGHT R( 5,L)=SIGM R( 6,L)=FWHM R( 7,L)=XLD R( 8,L)=FD R( 9,L)=XRD R(10,L)=FD R(11,L)=XL R(12,L)=FL R(13,L)=XR R(14,L)=FR R(15,L)=KL R(16,L)=KR 60 CONTINUE * ================================================================== 100 RETURN END SUBROUTINE BSPLN4(X,B,I0) * * Calculation of B-spline function values B(1) ... B(4) at X * for cubic B-splines with non-equidistant knots T(1)...T(NKNOT), * multiple knots are allowed * The first non-zero B-spline is the (I0+1)'th B-spline * X-region inside knots: T(4) ... T(NKNOT-3) * Number of coefficients: NKNOT-4 * 4 * S(X) = SUM A(I+I0) * B(I) * I=1 REAL B(4),DL(4),DR(4),SV PARAMETER (MBSP=100,MKNOT=MBSP+4) COMMON/CMBSP4/NKNOT,T(MKNOT) * ... * determination of index L for knot interval for x IL=4 IR=NKNOT-3 XX=AMIN1(AMAX1(X,T(IL)),T(IR)) IF(XX.LT.T(IR)) THEN * binary search 10 L=(IL+IR)/2 IF(XX.LT.T(L)) THEN IR=L GOTO 10 ELSE IF(XX.GT.T(L+1)) THEN IL=L+1 GOTO 10 END IF * case of multiple knots 20 IF(XX.EQ.T(L+1)) THEN L=L+1 GOTO 20 END IF ELSE * case of multiple knots right end L=IR 30 IF(XX.LE.T(L)) THEN L=L-1 GOTO 30 END IF END IF I0=L-4 * compute B(1)...B(4) for given x B(1)=1.0 DO J=1,3 DR(J)=T(L+J)-X DL(J)=X-T(L-J+1) SV=0.0 DO I=1,J TR=B(I)/(DR(I)+DL(J+1-I)) B(I)=SV+DR(I)*TR SV=DL(J+1-I)*TR END DO B(J+1)=SV END DO DO J=1,4 IF(B(J).LT.1.0E-14) B(J)=0.0 END DO 100 RETURN END SUBROUTINE BMDEQU(W,NBAND,NROW,B,DELTA) * * Solution of band matrix equation W * X = B * W(NBAND,NROW) = matrix in "band matrix storage mode * DELTA = scalar product X * F (= function change in fit) * * REAL W(NBAND,NROW),B(NROW),DIAG(257),RATIO,SUM * copy diagonal elements DO N=1,NROW DIAG(N)=W(1,N) END DO * factorisation DO N=1,NROW IF(W(1,N)+DIAG(N).LE.DIAG(N)) THEN * singular case DO J=1,NBAND W(J,N)=0.0 END DO ELSE W(1,N)=1.0/W(1,N) IMAX=MIN0(NBAND-1,NROW-N) JMAX=IMAX DO I=1,IMAX IF(ABS(W(I+1,N)).LT.1.E-18) THEN W(I+1,N)=0.0 END IF RATIO=W(I+1,N)*W(1,N) DO J=1,JMAX W(J,N+I)=W(J,N+I)-W(J+I,N)*RATIO END DO JMAX=JMAX-1 W(I+1,N)=RATIO END DO END IF END DO * copy vector B to auxiliary array DO N=1,NROW DIAG(N)=B(N) END DO * solution by forward ... DO N=1,NROW DO J=1,MIN0(NBAND-1,NROW-N) IF(ABS(W(J+1,N)).LT.1.E-18) THEN W(J+1,N)=0.0 END IF B(J+N)=B(J+N)-W(J+1,N)*B(N) END DO END DO * ... and backward substitution DO NN=1,NROW N=NROW+1-NN B(N)=B(N)*W(1,N) DO J=1,MIN0(NBAND-1,NROW-N) IF(ABS(W(J+1,N)).LT.1.E-18) THEN W(J+1,N)=0.0 END IF B(N)=B(N)-W(J+1,N)*B(J+N) END DO END DO * expected change in chi square SUM=0.0D0 DO N=1,NROW SUM=SUM+DIAG(N)*B(N) END DO DELTA=SUM END SUBROUTINE SPSARP(Z,N,MODE) * * Simple smoothing operations on array (N > 4) * * = -2 (3G53QH) equiv. to 1/2/3/4/5/6 * = -1 (353) equiv. to 1/2/3/ * MODE = 0 (0) set to zero * = 1 (3) replace by median of 3 with extrapolation at edges * = 2 (G) conditional hanning * = 3 (5) replace by median of 5 (of three at edges) * = 4 (3) replace by median of 3 with copy at edges * = 5 (Q) quadratic interpolation of flat 3-areas * = 6 (H) hanning (replace by average of 3) * = 7 (S) sum up (integrate) * = 8 (A) abs * = 9 (A) fouth derivative * REAL Z(*),ZI(3) * Median of values A,B,C DIAN3(A,B,C)=AMAX1(AMIN1(A,B),AMIN1(B,C),AMIN1(C,A)) * ... IF(N.LT.4.AND.MODE.NE.0) GOTO 100 IOPER=MODE IF(MODE+1.EQ.0) IOPER=1 IF(MODE+2.EQ.0) IOPER=2 01 IF(IOPER.EQ.0) THEN * reset to zero DO I=1,N Z(I)=0.0 END DO ELSE IF(IOPER.EQ.1) THEN * (3) replace each value by median of 3 values ZL =Z(1) DO I=1,N-2 ZN =DIAN3(Z(I),Z(I+1),Z(I+2)) Z(I) =ZL ZL =ZN END DO Z(N-1)=ZL * extrapolation at end Z(1) =DIAN3(Z(1),Z( 2),3.0*Z( 2)-2.0*Z( 3)) Z(N) =DIAN3(Z(N),Z(N-1),3.0*Z(N-1)-2.0*Z(N-2)) ELSE IF(IOPER.EQ.2) THEN * (G) conditional hanning ZL =Z(1) DO I=1,N-2 IF((Z(I+1)-Z(I))*(Z(I+1)-Z(I+2)).LT.0) THEN ZN = Z(I+1) ELSE ZN =0.25*(Z(I)+Z(I+1)+Z(I+1)+Z(I+2)) END IF Z(I) =ZL ZL =ZN END DO Z(N-1)=ZL ELSE IF(IOPER.EQ.3) THEN * (5) replace by median of 5 ZL=Z(1) ZN=DIAN3(Z(1),Z(2),Z(3)) ZE=DIAN3(Z(N-2),Z(N-1),Z(N)) ZM=0.0 DO I=1,N-4 MIN=0 IF(Z(I).GT.Z(I+4)) MIN=4 MAX=4-MIN DO K=1,3 IF(Z(I+K).LT.Z(I+MIN)) THEN ZI(K)=Z(I+MIN) MIN=K ELSE IF(Z(I+K).LT.Z(I+MAX)) THEN ZI(K)=Z(I+K) ELSE ZI(K)=Z(I+MAX) MAX=K END IF END IF END DO ZM =DIAN3(ZI(1),ZI(2),ZI(3)) Z(I) =ZL ZL =ZN ZN =ZM END DO Z(N-3)=ZN Z(N-2)=ZM Z(N-1)=ZE ELSE IF(IOPER.EQ.4) THEN * (3C) replace each value by median of 3 values ZL =Z(1) DO I=1,N-2 ZN =DIAN3(Z(I),Z(I+1),Z(I+2)) Z(I) =ZL ZL =ZN END DO Z(N-1)=ZL ELSE IF(IOPER.EQ.5) THEN * (Q) correct flat areas of three identical values by * quadratic interpolation DO I=3,N-2 IF(Z(I).EQ.Z(I-1).AND.Z(I).EQ.Z(I+1)) THEN IF((Z(I)-Z(I-2))*(Z(I)-Z(I+2)).GT.0) THEN J=1 IF(ABS(Z(I)-Z(I-2)).GT.ABS(Z(I)-Z(I+2))) J=-1 Z(I-J) =Z(I+J)-0.5*(Z(I+J+J)-Z(I-J-J)) Z(I) =(8.0*Z(I+J)-3.0*Z(I+J+J)+Z(I-J-J))/6.0 END IF END IF END DO ELSE IF(IOPER.EQ.6) THEN * (H) replace each value by average of 3 ZL =Z(1) DO I=1,N-2 ZN =0.25*(Z(I)+Z(I+1)+Z(I+1)+Z(I+2)) Z(I) =ZL ZL =ZN END DO Z(N-1)=ZL ELSE IF(IOPER.EQ.7) THEN * (S) sum up DO I=2,N Z(I)=Z(I-1)+Z(I) END DO ELSE IF(IOPER.EQ.8) THEN * (A) abs DO I=1,N Z(I)=ABS(Z(I)) END DO ELSE IF(IOPER.EQ.9) THEN ETA=0.5 ETA2=ETA*ETA * (F) fourth derivative SECN=0.5*(3.0*Z(N)-7.0*Z(N-1)+5.0*Z(N-2)-Z(N-3)) DERL=SIGN(SQRT(ETA2+SECN*SECN)-ETA,SECN) SECN=0.5*(3.0*Z(1)-7.0*Z(2)+5.0*Z(3)-Z(4)) DER2=SIGN(SQRT(ETA2+SECN*SECN)-ETA,SECN) DO I=1,N-3 * using 4-point formula SECN=Z(I)-Z(I+1)-Z(I+2)+Z(I+3) * ... with damping of small values Z(I)=DER2 DER2=SIGN(SQRT(ETA2+SECN*SECN)-ETA,SECN) END DO Z(N-2)=DER2 Z(N-1)=DERL * again second derivative SECN=0.5*(3.0*Z(N-1)-7.0*Z(N-2)+5.0*Z(N-3)-Z(N-4)) DERL=SIGN(SQRT(ETA2+SECN*SECN)-ETA,SECN) SECN=0.5*(3.0*Z(1)-7.0*Z(2)+5.0*Z(3)-Z(4)) DER2=SIGN(SQRT(ETA2+SECN*SECN)-ETA,SECN) DO I=1,N-4 * using 4-point formula SECN=Z(I)-Z(I+1)-Z(I+2)+Z(I+3) * ... with damping of small values Z(I)=DER2 DER2=SIGN(SQRT(ETA2+SECN*SECN)-ETA,SECN) END DO Z(N-3)=DER2 Z(N-2)=DERL DO I=N-2,1,-1 Z(I+1)=Z(I) END DO * natural edge values Z(1 )=0.0 Z(N )=0.0 END IF IF(MODE+1.EQ.0) THEN IOPER=IOPER+1 IF(IOPER.LE.3) GOTO 01 ELSE IF(MODE+2.EQ.0) THEN IOPER=IOPER+1 IF(IOPER.LE.6) GOTO 01 END IF 100 RETURN END