source: Sophya/trunk/SophyaExt/CodeMinuit/examples/fcnk0.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: 4.0 KB
Line 
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.