MODULE MyMinuit USE f90_kind PRIVATE ! wrapper for MINUIT ! initialisation and close files ! MINini(IOMIN,CHMIN,IOMOUT,CHMOUT,IOMSAV,CHMSAV) ! MINend() PUBLIC :: MINini,MINend ! minimization ! MINUIT(FCN,CHTITL,NFITP,PCH,PVAL,PERR,PMAT,BND1,BND2,CHI2,IPRT,IERR) ! LinearRegression(X,Y,SIGSQ,A,B,IERR) PUBLIC :: MINUIT, LinearRegression, Qfit ! spline approximation ! YOUT = SPLINE(X,Y,XIN) PUBLIC :: SPLINE PRIVATE :: FPRIME INTEGER,PARAMETER,PUBLIC :: minuitPrec=double INTEGER,SAVE,PRIVATE :: IOIN,IOOUT,IOSAV CONTAINS SUBROUTINE MINini(IOMIN,CHMIN,IOMOUT,CHMOUT,IOMSAV,CHMSAV) INTEGER,INTENT(IN) :: IOMIN,IOMOUT,IOMSAV CHARACTER(len=*),INTENT(IN) :: CHMIN,CHMOUT,CHMSAV INTEGER :: IOS ! save the units in global variables IOIN = IOMIN IOOUT = IOMOUT IOSAV = IOMSAV CALL OPENFILE(IOIN, CHMIN ,IOS) CALL OPENFILE(IOOUT,CHMOUT,IOS) CALL OPENFILE(IOSAV,CHMSAV,IOS) CALL MNINIT(IOIN,IOOUT,IOSAV) END SUBROUTINE MINini SUBROUTINE MINend() INTEGER :: IOS CALL CLOSEFILE(IOin,IOS) CALL CLOSEFILE(IOout,IOS) CALL CLOSEFILE(IOsav,IOS) END SUBROUTINE MINend SUBROUTINE MINUIT & (FCN,CHTITL,NFITP,PCH,PVAL,PERR,PMAT,BND1,BND2,CHI2,IPRT,IERR) !******************************************************************* !*** DOES THE FIT !*** INPUT: FCN : SUBROUTINE that calculates the chisq !*** CHTITL : title of the problem !*** NFITP : number of variable fit parameters !*** PCH : ARRAY: names of fitparameters !*** PVAL : ARRAY: start values of fitparameters !*** PERR : ARRAY: step width for fitparameters !*** if PERR.LE.0 the parameter is fixed !*** BND1,2 : ARRAY: boundaries, =0: no bound !*** IPRT : print level !*** OUTPUT: PVAL : ARRAY: fit values of fitparameters !*** PERR : ARRAY: errors of fitparameters !*** PMAT : MATRIX: covariance matrix !*** CHI2 : chisq !*** IERR : error flag, indicates quality of fit !******************************************************************* CHARACTER(len=*),INTENT(IN) :: CHTITL INTEGER,INTENT(IN) :: NFITP,IPRT CHARACTER(len=10),DIMENSION(NFITP),INTENT(IN) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP),INTENT(IN) :: BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP),INTENT(INOUT) :: PVAL,PERR REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP),INTENT(OUT) :: PMAT REAL(KIND=minuitPrec),INTENT(OUT) :: CHI2 INTEGER,INTENT(OUT) :: IERR INTERFACE SUBROUTINE FCN(NPAR,GRAD,FVAL,XVAL,IFLAG) USE f90_kind INTEGER,INTENT(IN) :: IFLAG , NPAR INTEGER,PARAMETER :: minuitPrec=double REAL(KIND=minuitPrec),DIMENSION(NPAR),INTENT(IN) :: XVAL REAL(KIND=minuitPrec),DIMENSION(NPAR),INTENT(OUT) :: GRAD REAL(KIND=minuitPrec),INTENT(OUT) :: FVAL END SUBROUTINE FCN END INTERFACE INTEGER :: I,IERFLG,NARG,ISTAT,NPARI,NPARX,IVARBL REAL(KIND=minuitPrec) :: FMIN,FEDM,ERRDEF REAL(KIND=minuitPrec),DIMENSION(10) :: ARGLIS LOGICAL,DIMENSION(NFITP) :: LFIX CHARACTER(len=10) :: CHNAM DO I=1,NFITP IF( PERR(I) <= 0.0_minuitPrec )THEN LFIX(I) = .TRUE. PERR(I) = ABS(PERR(I)) ELSE LFIX(I) = .FALSE. ENDIF ENDDO ! printout level: ! -1 : no printout ! 0 : minimum, results only ! 1 : default, input, steps and results ! 2 : + intermediate results ! 3 : maximum output ARGLIS(1) = IPRT CALL MNEXCM(FCN,"SET PRINTOUT",ARGLIS,1,IERFLG,0) CALL MNSETI(CHTITL) IF( IPRT >= 0 )THEN CALL WRITECH(IOout,"+++++++++++++++++++++++++++++++++++++++++") CALL WRITECH(IOout,CHTITL) CALL WRITECH(IOout,"+++++++++++++++++++++++++++++++++++++++++") ENDIF ! set CHI2 to infinity CHI2 = HUGE(CHI2) ! book parameters DO I=1,NFITP CALL MNPARM(I,PCH(I),PVAL(I),PERR(I),BND1(I),BND2(I),IERFLG) IF( IERFLG /= 0 )THEN IERR = -9 GOTO 999 ENDIF ENDDO ! fix parameters DO I=1,NFITP IF( LFIX(I) )THEN NARG = 1 ARGLIS(1) = I CALL MNEXCM(FCN,"FIX",ARGLIS,NARG,IERFLG,0) IF( IERFLG /= 0 )THEN IERR = -8 GOTO 999 ENDIF ENDIF ENDDO ! minimize NARG = 0 CALL MNEXCM(FCN,"MINIMIZE",ARGLIS,NARG,IERFLG,0) ! get the results DO I=1,NFITP CALL MNPOUT(I,CHNAM,PVAL(I),PERR(I),BND1,BND2,IVARBL) ENDDO CALL MNEMAT(PMAT,NFITP) ! quality of fit CALL MNSTAT(FMIN,FEDM,ERRDEF,NPARI,NPARX,ISTAT) CHI2 = FMIN SELECT CASE(ISTAT) CASE (3) ! good fit IERR = 0 CASE (2) ! full covariance matrix, but forced positive defitite IERR = 1 CASE (1) ! covariance matrix approximation only IERR = 2 CASE DEFAULT ! no covariance matrix IERR = 3 END SELECT 999 CONTINUE NARG = 0 CALL MNEXCM(FCN,"CLEAR",ARGLIS,NARG,IERFLG,0) END SUBROUTINE MINUIT SUBROUTINE LinearRegression(X,Y,SIGSQ,A,B,IERR) ! Y(I) = A + B*X(I), SIGSQ = error(Y)**2 ! see Particle Data Book REAL, DIMENSION(:), INTENT(IN) :: X,Y,SIGSQ REAL, INTENT(OUT) :: A,B INTEGER, INTENT(OUT) :: IERR REAL :: A11, A12, A22, g1, g2, D INTEGER :: I IERR = 0 A11 = 0.0 A12 = 0.0 A22 = 0.0 g1 = 0.0 g2 = 0.0 DO I=1,size(SIGSQ) IF( SIGSQ(I) /= 0.0 )THEN A11 = A11 + 1.0/SIGSQ(I) A12 = A12 + X(I)/SIGSQ(I) A22 = A22 + X(I)*X(I)/SIGSQ(I) g1 = g1 + Y(I)/SIGSQ(I) g2 = g2 + X(I)*Y(I)/SIGSQ(I) ENDIF ENDDO D = A11 * A22 - A12 * A12 IF (D /= 0.0) THEN A = (g1 * A22 - g2 * A12) / D B = (g2 * A11 - g1 * A12) / D ELSE A = 0.0 B = 0.0 IERR = 1 ENDIF END SUBROUTINE LinearRegression SUBROUTINE Qfit(Y,X0,SL) ! quartic fit, returns position of maximum/minimum ! via differentiation and Linear regression ! INPUT: Y y-values ! OUTPUT: X0 position of extremum ! SL slope of derivative: SL<0: maximum ! SL>0: minimum ! SL=0: fit failed REAL,DIMENSION(:),INTENT(IN) :: Y REAL,INTENT(OUT) :: X0,SL INTEGER :: I,N,IERR REAL :: A,B REAL,DIMENSION(SIZE(Y)-1) :: X,DY,S N = SIZE(X) X = (/ (I+0.5 , I=1,N) /) DY = Y(2:N+1) - Y(1:N) S = 1.0 CALL LinearRegression(X,DY,S,A,B,IERR) IF( IERR==0 .AND. B/=0.0 )THEN X0 = -A/B SL = B ! no extrapolation IF( X0 < 1.0 ) X0 = 1.0 IF( X0 > N+1.0 ) X0 = N+1.0 ELSE X0 = 0.0 SL = 0.0 ENDIF END SUBROUTINE Qfit !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FUNCTION SPLINE(X,Y,XIN) RESULT(YOUT) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccC !CC INPUT: X x-values with X(i) < X(i+1) !CC Y y-values !CC XIN x-value for which y will be returned !CC -------------------------------------------------------------- !CC WARNING: will crash if X(i) = X(i+1); exeption not caught !CC -------------------------------------------------------------- !CC Function returns y value for a corresponding x value, !CC based on cubic spline. !CC Will never oscillates or overshoot. No need to solve matrix. !CC Also calculate constants for cubic in case needed (for integration). !CC -------------------------------------------------------------- !CC adopted from VBasic routine !CC Developer: C Kruger, Guildford, UK !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccC REAL(KIND=minuitPrec),DIMENSION(:),INTENT(IN) :: X,Y REAL(KIND=minuitPrec),INTENT(IN) :: XIN REAL(KIND=minuitPrec) :: YOUT,A,B,C,D,FPI,FPIM1,F2I,F2IM1 INTEGER :: NPOINTS,LEG,I NPOINTS = SIZE(X) IF( NPOINTS <= 0 )THEN YOUT = 0.0 RETURN ELSEIF( NPOINTS == 1 )THEN YOUT = Y(1) RETURN ENDIF ! 1) Find LineNumber or segment/leg. Linear extrapolate if outside range. IF (XIN >= X(NPOINTS)) THEN ! B=( Y(NPOINTS)-Y(NPOINTS-1) )/( X(NPOINTS)-X(NPOINTS-1) ) ! A= Y(NPOINTS) - B*X(NPOINTS) ! ! YOUT = A + B*XIN YOUT = Y(NPOINTS) RETURN ELSEIF (XIN <= X(1)) THEN ! B=( Y(2)-Y(1) )/( X(2)-X(1) ) ! A= Y(2) - B*X(2) ! YOUT = A + B*XIN YOUT= Y(1) RETURN ELSEIF( NPOINTS == 2 )THEN B=( Y(2)-Y(1) )/( X(2)-X(1) ) A= Y(2) - B*X(2) YOUT = A + B*XIN RETURN ELSE LEG = NPOINTS DO I=2,NPOINTS IF (XIN <= X(I)) THEN LEG=I EXIT ENDIF ENDDO ENDIF ! 2) Calc first derivative (slope) for intermediate points ! 3) Reset first derivative (slope) at first and last point IF (LEG == 2) THEN ! First point has 0 2nd derivative FPI = FPRIME(X(LEG-1:LEG+1),Y(LEG-1:LEG+1)) FPIM1 = (3.0*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) ) - FPI)/2.0 ELSEIF (LEG == NPOINTS) THEN ! Last point has 0 2nd derivative FPIM1 = FPRIME(X(LEG-2:LEG),Y(LEG-2:LEG)) FPI = (3.0*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) ) - FPIM1)/2.0 ELSE FPIM1 = FPRIME(X(LEG-2:LEG ),Y(LEG-2:LEG )) FPI = FPRIME(X(LEG-1:LEG+1),Y(LEG-1:LEG+1)) ENDIF ! 4) Calc second derivative at points F2IM1 = -2.0*( FPI + 2.0*FPIM1 )/( X(LEG)-X(LEG-1) ) & + 6.0*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) )**2 F2I = 2.0*( 2.0*FPI + FPIM1 )/( X(LEG)-X(LEG-1) ) & - 6.0*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) )**2 ! 5) Calc constants for cubic D = ( F2I-F2IM1 )/( 6.0*(X(LEG)-X(LEG-1)) ) C = ( X(LEG)*F2IM1 - X(LEG-1)*F2I )/( 2.0*( X(LEG)-X(LEG-1 )) ) B = ( Y(LEG)-Y(LEG-1) - C*( X(LEG)**2-X(LEG-1)**2 ) - D*( X(LEG)**3-X(LEG-1)**3 ) ) & / ( X(LEG)-X(LEG-1) ) A = Y(LEG-1) - B*X(LEG-1) - C*X(LEG-1)**2 - D*X(LEG-1)**3 YOUT = A + B*XIN + C*XIN**2 + D*XIN**3 END FUNCTION SPLINE FUNCTION FPRIME(X,Y) RESULT(FP) REAL(KIND=minuitPrec),DIMENSION(3),INTENT(IN) :: X,Y REAL(KIND=minuitPrec) :: FP,DX1,DX2,DY1,DY2 DX1 = X(2)-X(1) DX2 = X(3)-X(2) DY1 = Y(2)-Y(1) DY2 = Y(3)-Y(2) IF( DY1==0.0 .OR. DY2==0.0 )THEN ! Prevent divide by 0 FP = 0.0 ELSEIF( DY1*DY2 < 0.0 )THEN ! Pos AND neg slope, assume slope = 0 to prevent overshoot FP = 0.0 ELSE FP = 2.0/( DX1/DY1 + DX2/DY2 ) ENDIF END FUNCTION FPRIME END MODULE MyMinuit