source: Sophya/trunk/SophyaExt/CodeMinuit/examples/fcnk0.F@ 3891

Last change on this file since 3891 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: 4.0 KB
RevLine 
[2403]1*
2* $Id: fcnk0.F,v 1.1.1.1 2003-06-11 14:18:31 cmv Exp $
3*
4* $Log: not supported by cvs2svn $
5* Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
6* Minuit
7*
8*
9#include "minuit/pilot.h"
10 SUBROUTINE FCNK0(NPAR,GIN,F,X,IFLAG)
11#include "minuit/d506dp.inc"
12 REAL THPLUI, THMINI
13 DIMENSION X(*),GIN(*)
14 PARAMETER (MXBIN=50)
15 DIMENSION THPLU(MXBIN),THMIN(MXBIN),T(MXBIN),
16 + EVTP(MXBIN),EVTM(MXBIN)
17 DATA NBINS,NEVTOT/ 30,250/
18 DATA (EVTP(IGOD),IGOD=1,30)
19 + /11., 9., 13., 13., 17., 9., 1., 7., 8., 9.,
20 + 6., 4., 6., 3., 7., 4., 7., 3., 8., 4.,
21 + 6., 5., 7., 2., 7., 1., 4., 1., 4., 5./
22 DATA (EVTM(IGOD),IGOD=1,30)
23 + / 0., 0., 0., 0., 0., 0., 0., 0., 1., 1.,
24 + 0., 2., 1., 4., 4., 2., 4., 2., 2., 0.,
25 + 2., 3., 7., 2., 3., 6., 2., 4., 1., 5./
26C
27 XRE = X(1)
28 XIM = X(2)
29 DM = X(5)
30 GAMS = 1.0/X(10)
31 GAML = 1.0/X(11)
32 GAMLS = 0.5*(GAML+GAMS)
33 IF (IFLAG .NE. 1) GO TO 300
34C generate random data
35 STHPLU = 0.
36 STHMIN = 0.
37 DO 200 I= 1, NBINS
38 T(I) = 0.1*REAL(I)
39 TI = T(I)
40 EHALF = EXP(-TI*GAMLS)
41 TH = ((1.0-XRE)**2 + XIM**2) * EXP(-TI*GAML)
42 TH = TH + ((1.0+XRE)**2 + XIM**2) * EXP(-TI*GAMS)
43 TH = TH - 4.0*XIM*SIN(DM*TI) * EHALF
44 STERM = 2.0*(1.0-XRE**2-XIM**2)*COS(DM*TI) * EHALF
45 THPLU(I) = TH + STERM
46 THMIN(I) = TH - STERM
47 STHPLU = STHPLU + THPLU(I)
48 STHMIN = STHMIN + THMIN(I)
49 200 CONTINUE
50 NEVPLU = REAL(NEVTOT)*(STHPLU/(STHPLU+STHMIN))
51 NEVMIN = REAL(NEVTOT)*(STHMIN/(STHPLU+STHMIN))
52 WRITE (6,'(A)') ' LEPTONIC K ZERO DECAYS'
53 WRITE (6,'(A,3I10)') ' PLUS, MINUS, TOTAL=',NEVPLU,NEVMIN,NEVTOT
54 WRITE (6,'(A)')
55 + '0 TIME THEOR+ EXPTL+ THEOR- EXPTL-'
56 SEVTP = 0.
57 SEVTM = 0.
58 DO 250 I= 1, NBINS
59 THPLU(I) = THPLU(I)*REAL(NEVPLU) / STHPLU
60 THMIN(I) = THMIN(I)*REAL(NEVMIN) / STHMIN
61 THPLUI = THPLU(I)
62CCC CALL POISSN(THPLUI,NP,IERROR)
63CCC EVTP(I) = NP
64 SEVTP = SEVTP + EVTP(I)
65 THMINI = THMIN(I)
66CCC CALL POISSN(THMINI,NM,IERROR)
67CCC EVTM(I) = NM
68 SEVTM = SEVTM + EVTM(I)
69 IF (IFLAG .NE. 4)
70 + WRITE (6,'(1X,5G12.4)') T(I),THPLU(I),EVTP(I),THMIN(I),EVTM(I)
71 250 CONTINUE
72 WRITE (6, '(A,2F10.2)') ' DATA EVTS PLUS, MINUS=', SEVTP,SEVTM
73C calculate chisquared
74 300 CONTINUE
75 CHISQ = 0.
76 STHPLU = 0.
77 STHMIN = 0.
78 DO 400 I= 1, NBINS
79 TI = T(I)
80 EHALF = EXP(-TI*GAMLS)
81 TH = ((1.0-XRE)**2 + XIM**2) * EXP(-TI*GAML)
82 TH = TH + ((1.0+XRE)**2 + XIM**2) * EXP(-TI*GAMS)
83 TH = TH - 4.0*XIM*SIN(DM*TI) * EHALF
84 STERM = 2.0*(1.0-XRE**2-XIM**2)*COS(DM*TI) * EHALF
85 THPLU(I) = TH + STERM
86 THMIN(I) = TH - STERM
87 STHPLU = STHPLU + THPLU(I)
88 STHMIN = STHMIN + THMIN(I)
89 400 CONTINUE
90 THP = 0.
91 THM = 0.
92 EVP = 0.
93 EVM = 0.
94 IF (IFLAG .NE. 4) WRITE (6,'(1H0,10X,A,20X,A)')
95 + 'POSITIVE LEPTONS','NEGATIVE LEPTONS'
96 IF (IFLAG .NE. 4) WRITE (6,'(A,3X,A)')
97 + ' TIME THEOR EXPTL CHISQ',
98 + ' TIME THEOR EXPTL CHISQ'
99C
100 DO 450 I= 1, NBINS
101 THPLU(I) = THPLU(I)*SEVTP / STHPLU
102 THMIN(I) = THMIN(I)*SEVTM / STHMIN
103 THP = THP + THPLU(I)
104 THM = THM + THMIN(I)
105 EVP = EVP + EVTP(I)
106 EVM = EVM + EVTM(I)
107C Sum over bins until at least four events found
108 IF (EVP .GT. 3.) THEN
109 CHI1 = (EVP-THP)**2/EVP
110 CHISQ = CHISQ + CHI1
111 IF (IFLAG .NE. 4)
112 + WRITE (6,'(1X,4F9.3)') T(I),THP,EVP,CHI1
113 THP = 0.
114 EVP = 0.
115 ENDIF
116 IF (EVM .GT. 3) THEN
117 CHI2 = (EVM-THM)**2/EVM
118 CHISQ = CHISQ + CHI2
119 IF (IFLAG .NE. 4)
120 + WRITE (6,'(42X,4F9.3)') T(I),THM,EVM,CHI2
121 THM = 0.
122 EVM = 0.
123 ENDIF
124 450 CONTINUE
125 F = CHISQ
126 RETURN
127 END
Note: See TracBrowser for help on using the repository browser.