PROGRAM MINU02 C C Example of program using MINUIT (cernlib D506) in Fortran-callable mode C (see Reference manual version 94.1 section 3.3 and example in section 6.2) C C INPUTS: unit 5 (minu02.input) for COMMAND input ONLY in data-driven mode C histograms (minu02.hrin) for DATA input C C OUTPUTS: unit 6 (minu02.output) for MINUIT and user output C C (C) Luciano Ramello 1999, 2009 C PARAMETER (NWHB=1000000) COMMON /PAWC/ PIPPO(NWHB) C.. Declarations for the MINUIT fits EXTERNAL FCN CHARACTER*50 TITLE DOUBLE PRECISION ARGLIS(10) PARAMETER (MAXPAR=4) DIMENSION NPRM(MAXPAR) CHARACTER*10 PNAM(MAXPAR) DOUBLE PRECISION STRT(MAXPAR), STEP(MAXPAR), BLOW(MAXPAR), & BHIG(MAXPAR) DOUBLE PRECISION ZERO,ONE,TWO,THREE,FOUR,FIVE ! numeric arguments for MNEXCM C.. Common with some additional user control parameters COMMON/USRPAR/NTOTPR,XMIN,XMAX ! number of param., fit range REAL XMIN,XMAX DATA TITLE / 'Fit of muon lifetime ' / DATA NPRM / 1, 2, 3, 4 / DATA PNAM / 'Const1', 'Slope1', 'Const2', 'Slope2' / DATA STRT / 4.0, -0.5, 10., -10. / DATA STEP / 1., 0.1, 1., 1. / DATA BLOW / 0.1, -5., 1., -100. / DATA BHIG / 50., -0.1, 100., -1. / DATA ZERO,ONE,TWO,THREE,FOUR,FIVE / 0.,1.,2.,3.,4.,5. / C.. Initialize HBOOK working space CALL HLIMIT(NWHB) C.. Initialize I/O C.. use the default logical unit numbers (5, 6 and 7) CALL MNINIT(5,6,7) ! define unit no's for INPUT, OUTPUT and SAVE-command C.. open unit 5 only for data-driven mode (not for Fortran-callable mode) C OPEN(UNIT=5,FILE='minu02.input',STATUS='OLD') OPEN(UNIT=6,FILE='minu02.output',STATUS='NEW',FORM='FORMATTED') C.. Specify a title CALL MNSETI(TITLE) C.. Define the number of parameters and range of the fit NTOTPR=4 ! number of fit param.: must not exceed MAXPAR XMIN=0.4 ! lower limit of fit (microseconds) XMAX=4.5 ! upper limit of fit (microseconds) C.. Define the MINUIT parameters DO I=1,NTOTPR CALL MNPARM(NPRM(I),PNAM(I),STRT(I),STEP(I) &, BLOW(I),BHIG(I),IERFLG) IF(IERFLG.NE.0) THEN WRITE(6,*) ' MNPARM: unable to set param. ',I STOP ENDIF ENDDO ! I=1,NTOTPR WRITE(6,*) ' Main: MINUIT parameters have been set' C.. Execute MINUIT commands (third arg. "ARGLIS" must be dimensioned when command needs it) CALL MNEXCM(FCN,'CALL FCN', ONE, 1,IERFLG,0) C CALL MNEXCM(FCN,'SIMPLEX', ZERO, 0,IERFLG,0) C CALL MNEXCM(FCN,'MIGRAD', ZERO, 0,IERFLG,0) CALL MNEXCM(FCN,'MINIMIZE', ZERO, 0,IERFLG,0) CALL MNEXCM(FCN,'CALL FCN',THREE, 1,IERFLG,0) ARGLIS(1)=2. ! first param. ARGLIS(2)=4. ! second param. ARGLIS(3)=2. ! number of standard dev. ARGLIS(4)=50. ! number of grid points CALL MNEXCM(FCN,'CONTOUR',ARGLIS,4,IERFLG,0) C.. Close output file(s) CLOSE(UNIT=6) STOP END SUBROUTINE FCN(NPAR,DERIV,FCHISQ,XPAR,IFLAG) C-- MINUIT user subroutine C-- The function to be minimized in this case is a CHI-square IMPLICIT DOUBLE PRECISION(A-H, O-Z) DIMENSION DERIV(*),XPAR(*) C-- Declare explicitly REAL variables if double precision is unwanted REAL X1,X2,X3,X4,CHICHI REAL XAR(200),YAR(200),DXAR(200),DYAR(200) SAVE XAR,YAR,DXAR,DYAR ! save values across calls REAL XMI,XMA,YMI,YMA CHARACTER*80 CHTITL C-- Save integer variable NDATA between calls: INTEGER NDATA SAVE NDATA COMMON/USRPAR/NTOTPR,XMIN,XMAX ! number of param., fit range REAL XMIN,XMAX IF(IFLAG.EQ.1) THEN C-- Read input data: WRITE(6,*)' Reading input data from histogram file' C-- Get data from histogram file --> arrays XAR,YAR (errors DXAR,DYAR) ID=2 CALL HRGET(ID,'minu02.hrin',' ') CALL HGIVE(ID,CHTITL,NX,XMI,XMA,NY,YMI,YMA,NWT,LOC) CALL HREBIN(ID,XAR,YAR,DXAR,DYAR,NX,1,NX) CALL HDELET(ID) C-- NDATA=0 DO KK=1,NX C-- Exclude bad points (empty bin or too big relative error) C-- and points outside fit range: IF(DYAR(KK).GT.0..AND.DYAR(KK).LE.0.50*YAR(KK)) THEN IF(XAR(KK).GE.XMIN.AND.XAR(KK).LE.XMAX) THEN NDATA=NDATA+1 XAR(NDATA)=XAR(KK) YAR(NDATA)=YAR(KK) DXAR(NDATA)=DXAR(KK) DYAR(NDATA)=DYAR(KK) ENDIF ! XAR(KK).GE.XMIN ... ENDIF ! DYAR(KK).GT.0 ... ENDDO ! KK=1,NX WRITE(6,*) ' FCN: retained ',NDATA,' data points' ELSE IF(IFLAG.EQ.2) THEN C-- Calculate the first derivatives (optional) ENDIF C-- Always calculate the value of FCHISQ: FCHISQ=0.0 NDF=0 DO I=1,NDATA C C-- Function to be fitted: sum of two exponentials (or simple exponential) C FIT=EXP(XPAR(1)+XPAR(2)*XAR(I)) ! use PAW convention for expon. fits IF(NTOTPR.EQ.4) THEN FIT=FIT+EXP(XPAR(3)+XPAR(4)*XAR(I)) ENDIF NDF=NDF+1 FCHISQ=FCHISQ+(YAR(I)-FIT)*(YAR(I)-FIT)/DYAR(I)/DYAR(I) ENDDO ! i=1,ndata NDF=NDF-NTOTPR ! compute n. of degrees of freedom CHICHI=FCHISQ ! save the chi**2 value IF(IFLAG.EQ.3) THEN C-- Perform final calculations, output fitted data, etc.: X1=XPAR(1) X2=XPAR(2) X3=XPAR(3) X4=XPAR(4) C-- Print a resume' of the fit CALL RESUME(X1,X2,X3,X4,CHICHI,NDF) ENDIF RETURN END SUBROUTINE RESUME(X1,X2,X3,X4,CHICHI,NDF) C C-- Print short "resume'" of fit: C C-- X1 = param. no. 1 C-- X2 = param. no. 2 C-- X3 = param. no. 3 C-- X4 = param. no. 4 C COMMON/USRPAR/NTOTPR,XMIN,XMAX ! number of param., fit range REAL XMIN,XMAX CHINDF=CHICHI/NDF WRITE(6,1) CHICHI,CHINDF,NDF 1 FORMAT(/' 2' / & ' CHI ',E12.5 / & ' ---- = ------------ = ',F10.3 / & ' NDF ',I6 ) WRITE(6,2) X1,X2,X3,X4 2 FORMAT(/' P1= ',E12.5,' P2= ',E12.5, & ' P3= ',E12.5,' P4= ',E12.5) RETURN END