source: Sophya/trunk/SophyaExt/CodeMinuit/code/mnmnot.F@ 3369

Last change on this file since 3369 was 2403, checked in by cmv, 22 years ago

Creation du module de code source de MINUIT (CERNLIB) extrait par CMV

cmv 11/06/2003
File size: 5.4 KB
RevLine 
[2403]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"
12CC Performs a MINOS error analysis on one parameter.
13CC The parameter ILAX is varied, and the minimum of the
14CC function with respect to the other parameters is followed
15CC until it crosses the value FMIN+UP.
16CC
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')
22C . . 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
59C 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
66C . . . . . Nota Bene: from here on, NPAR=MPAR-1
67C Remember: MNFIXP squeezes IT out of X, XT, WERR, and VHMAT,
68C 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
77C . 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
89C 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)
96C loop to hit AMIN+UP
97 KE1CR = ILAX
98 KE2CR = 0
99 XMIDCR = U(ILAX)
100 XDIRCR = DELU
101C
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
109C . 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
115C . . . . . . . . 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.
129C
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
139C . . parameter finished. reset v
140C 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
156C 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
167C in any case
168 700 CONTINUE
169 ITAUR = 0
170 NFCNMX = NFMXIN
171 ISTRAT = ISTRAV
172 RETURN
173 END
Note: See TracBrowser for help on using the repository browser.