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

Last change on this file since 3632 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.