source: Sophya/trunk/SophyaExt/CodeMinuit/code/mnpfit.F

Last change on this file 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: 1.6 KB
RevLine 
[2403]1*
2* $Id: mnpfit.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:31 mclareni
6* Minuit
7*
8*
9#include "minuit/pilot.h"
10 SUBROUTINE MNPFIT(PARX2P,PARY2P,NPAR2P,COEF2P,SDEV2P)
11#include "minuit/d506dp.inc"
12C
13C to fit a parabola to npar2p points
14C
15C npar2p no. of points
16C parx2p(i) x value of point i
17C pary2p(i) y value of point i
18C
19C coef2p(1...3) coefficients of the fitted parabola
20C y=coef2p(1) + coef2p(2)*x + coef2p(3)*x**2
21C sdev2p= variance
22C method : chi**2 = min equation solved explicitly
23 DIMENSION PARX2P(NPAR2P),PARY2P(NPAR2P),COEF2P(NPAR2P)
24 DIMENSION CZ(3)
25C
26 DO 3 I=1,3
27 3 CZ(I)=0.
28 SDEV2P=0.
29 IF(NPAR2P.LT.3) GO TO 10
30 F=NPAR2P
31C--- center x values for reasons of machine precision
32 XM=0.
33 DO 2 I=1,NPAR2P
34 2 XM=XM+PARX2P(I)
35 XM=XM/F
36 X2=0.
37 X3=0.
38 X4=0.
39 Y=0.
40 Y2=0.
41 XY=0.
42 X2Y=0.
43 DO 1 I=1,NPAR2P
44 S=PARX2P(I)-XM
45 T=PARY2P(I)
46 S2=S*S
47 X2=X2+S2
48 X3=X3+S*S2
49 X4=X4+S2*S2
50 Y=Y+T
51 Y2=Y2+T*T
52 XY=XY+S*T
53 X2Y=X2Y+S2*T
54 1 CONTINUE
55 A=(F*X4-X2**2)*X2-F*X3**2
56 IF(A.EQ.0.) GOTO 10
57 CZ(3)=(X2*(F*X2Y-X2*Y)-F*X3*XY)/A
58 CZ(2)=(XY-X3*CZ(3))/X2
59 CZ(1)=(Y-X2*CZ(3))/F
60 IF(NPAR2P.EQ.3) GOTO 6
61 SDEV2P=Y2-(CZ(1)*Y+CZ(2)*XY+CZ(3)*X2Y)
62 IF(SDEV2P.LT.0.) SDEV2P=0.
63 SDEV2P=SDEV2P/(F-3.)
64 6 CZ(1)=CZ(1)+XM*(XM*CZ(3)-CZ(2))
65 CZ(2)=CZ(2)-2.*XM*CZ(3)
66 10 CONTINUE
67 DO 11 I=1,3
68 11 COEF2P(I)=CZ(I)
69 RETURN
70 END
Note: See TracBrowser for help on using the repository browser.