1 | *
|
---|
2 | * $Id: mnmnot.F,v 1.1.1.1 2003-06-11 14:18:28 cmv Exp $
|
---|
3 | *
|
---|
4 | * $Log: not supported by cvs2svn $
|
---|
5 | * Revision 1.1.1.1 1996/03/07 14:31:30 mclareni
|
---|
6 | * Minuit
|
---|
7 | *
|
---|
8 | *
|
---|
9 | #include "minuit/pilot.h"
|
---|
10 | SUBROUTINE MNMNOT(FCN,ILAX,ILAX2,VAL2PL,VAL2MI,FUTIL)
|
---|
11 | #include "minuit/d506dp.inc"
|
---|
12 | CC Performs a MINOS error analysis on one parameter.
|
---|
13 | CC The parameter ILAX is varied, and the minimum of the
|
---|
14 | CC function with respect to the other parameters is followed
|
---|
15 | CC until it crosses the value FMIN+UP.
|
---|
16 | CC
|
---|
17 | #include "minuit/d506cm.inc"
|
---|
18 | EXTERNAL FCN,FUTIL
|
---|
19 | DIMENSION XDEV(MNI),W(MNI),GCC(MNI)
|
---|
20 | CHARACTER*4 CPOS,CNEG,CSIG
|
---|
21 | PARAMETER (CPOS='POSI',CNEG='NEGA')
|
---|
22 | C . . save and prepare start vals
|
---|
23 | ISW2 = ISW(2)
|
---|
24 | ISW4 = ISW(4)
|
---|
25 | SIGSAV = EDM
|
---|
26 | ISTRAV = ISTRAT
|
---|
27 | DC = DCOVAR
|
---|
28 | LNEWMN = .FALSE.
|
---|
29 | APSI = EPSI*0.5
|
---|
30 | ABEST=AMIN
|
---|
31 | MPAR=NPAR
|
---|
32 | NFMXIN = NFCNMX
|
---|
33 | DO 125 I= 1, MPAR
|
---|
34 | 125 XT(I) = X(I)
|
---|
35 | DO 130 J= 1, MPAR*(MPAR+1)/2
|
---|
36 | 130 VTHMAT(J) = VHMAT(J)
|
---|
37 | DO 135 I= 1, MPAR
|
---|
38 | GCC(I) = GLOBCC(I)
|
---|
39 | 135 W(I) = WERR(I)
|
---|
40 | IT = NIOFEX(ILAX)
|
---|
41 | ERP(IT) = 0.
|
---|
42 | ERN(IT) = 0.
|
---|
43 | CALL MNINEX(XT)
|
---|
44 | UT = U(ILAX)
|
---|
45 | IF (NVARL(ILAX) .EQ. 1) THEN
|
---|
46 | ALIM(ILAX) = UT -100.*W(IT)
|
---|
47 | BLIM(ILAX) = UT +100.*W(IT)
|
---|
48 | ENDIF
|
---|
49 | NDEX = IT*(IT+1)/2
|
---|
50 | XUNIT = SQRT(UP/VTHMAT(NDEX))
|
---|
51 | MARC = 0
|
---|
52 | DO 162 I= 1, MPAR
|
---|
53 | IF (I .EQ. IT) GO TO 162
|
---|
54 | MARC = MARC + 1
|
---|
55 | IMAX = MAX(IT,I)
|
---|
56 | INDX = IMAX*(IMAX-1)/2 + MIN(IT,I)
|
---|
57 | XDEV(MARC) = XUNIT*VTHMAT(INDX)
|
---|
58 | 162 CONTINUE
|
---|
59 | C fix the parameter in question
|
---|
60 | CALL MNFIXP (IT,IERR)
|
---|
61 | IF (IERR .GT. 0) THEN
|
---|
62 | WRITE (ISYSWR,'(A,I5,A,I5)')
|
---|
63 | + ' MINUIT ERROR. CANNOT FIX PARAMETER',ILAX,' INTERNAL',IT
|
---|
64 | GO TO 700
|
---|
65 | ENDIF
|
---|
66 | C . . . . . Nota Bene: from here on, NPAR=MPAR-1
|
---|
67 | C Remember: MNFIXP squeezes IT out of X, XT, WERR, and VHMAT,
|
---|
68 | C not W, VTHMAT
|
---|
69 | DO 500 ISIG= 1,2
|
---|
70 | IF (ISIG .EQ. 1) THEN
|
---|
71 | SIG = 1.0
|
---|
72 | CSIG = CPOS
|
---|
73 | ELSE
|
---|
74 | SIG = -1.0
|
---|
75 | CSIG = CNEG
|
---|
76 | ENDIF
|
---|
77 | C . sig=sign of error being calcd
|
---|
78 | IF (ISW(5) .GT. 1) WRITE (ISYSWR,806) CSIG,ILAX,CPNAM(ILAX)
|
---|
79 | 806 FORMAT (/' DETERMINATION OF ',A4,'TIVE MINOS ERROR FOR PARAMETER',
|
---|
80 | + I3, 2X ,A)
|
---|
81 | IF (ISW(2).LE.0) CALL MNWARN('D','MINOS','NO COVARIANCE MATRIX.')
|
---|
82 | NLIMIT = NFCN + NFMXIN
|
---|
83 | ISTRAT = MAX(ISTRAV-1,0)
|
---|
84 | DU1 = W(IT)
|
---|
85 | U(ILAX) = UT + SIG*DU1
|
---|
86 | U(ILAX) = MIN(U(ILAX),BLIM(ILAX))
|
---|
87 | U(ILAX) = MAX(U(ILAX),ALIM(ILAX))
|
---|
88 | DELU = U(ILAX) - UT
|
---|
89 | C stop if already at limit with negligible step size
|
---|
90 | IF (ABS(DELU)/(ABS(UT)+ABS(U(ILAX))) .LT. EPSMAC) GO TO 440
|
---|
91 | FAC = DELU/W(IT)
|
---|
92 | DO 185 I= 1, NPAR
|
---|
93 | 185 X(I) = XT(I) + FAC*XDEV(I)
|
---|
94 | IF (ISW(5) .GT. 1) WRITE (ISYSWR,801) ILAX,UT,DELU,U(ILAX)
|
---|
95 | 801 FORMAT (/' PARAMETER',I4,' SET TO',E11.3,' + ',E10.3,' = ',E12.3)
|
---|
96 | C loop to hit AMIN+UP
|
---|
97 | KE1CR = ILAX
|
---|
98 | KE2CR = 0
|
---|
99 | XMIDCR = U(ILAX)
|
---|
100 | XDIRCR = DELU
|
---|
101 | C
|
---|
102 | AMIN = ABEST
|
---|
103 | NFCNMX = NLIMIT - NFCN
|
---|
104 | CALL MNCROS(FCN,AOPT,IERCR,FUTIL)
|
---|
105 | IF (ABEST-AMIN .GT. 0.01*UP) GO TO 650
|
---|
106 | IF (IERCR .EQ. 1) GO TO 440
|
---|
107 | IF (IERCR .EQ. 2) GO TO 450
|
---|
108 | IF (IERCR .EQ. 3) GO TO 460
|
---|
109 | C . error successfully calculated
|
---|
110 | EROS = XMIDCR-UT + AOPT*XDIRCR
|
---|
111 | IF (ISW(5) .GT. 1) WRITE (ISYSWR,808) CSIG,ILAX,CPNAM(ILAX),EROS
|
---|
112 | 808 FORMAT (/9X,4HTHE ,A4, 29HTIVE MINOS ERROR OF PARAMETER,I3, 2H
|
---|
113 | +, ,A10, 4H, IS ,E12.4)
|
---|
114 | GO TO 480
|
---|
115 | C . . . . . . . . failure returns
|
---|
116 | 440 IF (ISW(5) .GE. 1) WRITE(ISYSWR,807) CSIG,ILAX,CPNAM(ILAX)
|
---|
117 | 807 FORMAT (5X,'THE ',A4,'TIVE MINOS ERROR OF PARAMETER',I3,', ',A,
|
---|
118 | +', EXCEEDS ITS LIMIT.'/)
|
---|
119 | EROS = UNDEFI
|
---|
120 | GO TO 480
|
---|
121 | 450 IF (ISW(5) .GE. 1) WRITE (ISYSWR, 802) CSIG,ILAX,NFMXIN
|
---|
122 | 802 FORMAT (9X,'THE ',A,'TIVE MINOS ERROR',I4,' REQUIRES MORE THAN',
|
---|
123 | + I5,' FUNCTION CALLS.'/)
|
---|
124 | EROS = 0.
|
---|
125 | GO TO 480
|
---|
126 | 460 IF (ISW(5) .GE. 1) WRITE (ISYSWR, 805) CSIG,ILAX
|
---|
127 | 805 FORMAT (25X,A,'TIVE MINOS ERROR NOT CALCULATED FOR PARAMETER',I4/)
|
---|
128 | EROS = 0.
|
---|
129 | C
|
---|
130 | 480 IF (ISW(5) .GT. 1) WRITE (ISYSWR,'(5X, 74(1H*))')
|
---|
131 | IF (SIG .LT. ZERO) THEN
|
---|
132 | ERN(IT) = EROS
|
---|
133 | IF (ILAX2.GT.0 .AND. ILAX2.LE.NU) VAL2MI = U(ILAX2)
|
---|
134 | ELSE
|
---|
135 | ERP(IT) = EROS
|
---|
136 | IF (ILAX2.GT.0 .AND. ILAX2.LE.NU) VAL2PL = U(ILAX2)
|
---|
137 | ENDIF
|
---|
138 | 500 CONTINUE
|
---|
139 | C . . parameter finished. reset v
|
---|
140 | C normal termination
|
---|
141 | ITAUR = 1
|
---|
142 | CALL MNFREE(1)
|
---|
143 | DO 550 J= 1, MPAR*(MPAR+1)/2
|
---|
144 | 550 VHMAT(J) = VTHMAT(J)
|
---|
145 | DO 595 I= 1, MPAR
|
---|
146 | WERR(I) = W(I)
|
---|
147 | GLOBCC(I) = GCC(I)
|
---|
148 | 595 X(I) = XT(I)
|
---|
149 | CALL MNINEX (X)
|
---|
150 | EDM = SIGSAV
|
---|
151 | AMIN = ABEST
|
---|
152 | ISW(2) = ISW2
|
---|
153 | ISW(4) = ISW4
|
---|
154 | DCOVAR = DC
|
---|
155 | GO TO 700
|
---|
156 | C new minimum
|
---|
157 | 650 LNEWMN = .TRUE.
|
---|
158 | ISW(2) = 0
|
---|
159 | DCOVAR = 1.
|
---|
160 | ISW(4) = 0
|
---|
161 | SAV = U(ILAX)
|
---|
162 | ITAUR = 1
|
---|
163 | CALL MNFREE(1)
|
---|
164 | U(ILAX) = SAV
|
---|
165 | CALL MNEXIN(X)
|
---|
166 | EDM = BIGEDM
|
---|
167 | C in any case
|
---|
168 | 700 CONTINUE
|
---|
169 | ITAUR = 0
|
---|
170 | NFCNMX = NFMXIN
|
---|
171 | ISTRAT = ISTRAV
|
---|
172 | RETURN
|
---|
173 | END
|
---|