]>
Commit | Line | Data |
---|---|---|
c9c6bd3b | 1 | C ********************************************************************** |
2 | C Version 2.0 | |
3 | C Generator of e+e- pairs produced in 1.10.2002 | |
4 | C PbPb collisions at LHC | |
5 | C | |
6 | C Copyright (c) 2002 | |
7 | C Authors: Kai Hencken <k.hencken@unibas.ch> | |
8 | C Yuri Kharlov <Yuri.Kharlov@cern.ch> | |
9 | C Serguei Sadovsky <sadovsky@mx.ihep.su>, | |
10 | C <Serguei.Sadovski@cern.ch> | |
11 | C | |
12 | C Permission to use, copy and distribute this software and its | |
13 | C documentation strictly for non-commercial purposes is hereby granted | |
14 | C without fee, provided that the above copyright notice appears in all | |
15 | C copies and that both the copyright notice and this permission notice | |
16 | C appear in the supporting documentation. The authors make no claims | |
17 | C about the suitability of this software for any purpose. It is | |
18 | C provided "as is" without express or implied warranty. | |
19 | C Any change of the code should be submitted to the authors | |
20 | C | |
21 | C Long write up: K.Hencken,Yu.Kharlov,S.Sadovsky | |
22 | C Internal ALICE Note 2002-27 | |
23 | C | |
24 | C ********************************************************************** | |
25 | subroutine ee_init(Ymin,Ymax,PTmin,PTmax) | |
26 | C----------------------------------------------------------------------- | |
27 | C Generator initialisation | |
28 | C | |
29 | C Input variables: Ymin - minimal value of rapidity \ | |
30 | C Ymax - maximal value of rapidity | of kinematics | |
31 | C PTmin - Pt minimum in MeV/c | range | |
32 | C PTmax - Pt maximum in MeV/c / | |
33 | C----------------------------------------------------------------------- | |
34 | implicit real*8 (A-H,O-Z) | |
35 | external DsdYpY,DsdYmY,DsdXpX,DsdXmX,DsdPhi,DsdXX,DsdYY | |
36 | ||
37 | common/parPh / parPhi(7) | |
38 | common/parYp / parYpY(6) | |
39 | common/parYm / parYmY(6) | |
40 | common/parXp / parXpX(5) | |
41 | common/parXm / parXmX(4) | |
42 | ||
43 | common/eevent/ Xsect2,Dsect2, Xsecttot,Dsecttot, Nevnt | |
44 | common/eepars/ Xmin,Xmax,YpYmin,YpYmax,WYpYmax,YmYmin,YmYmax, | |
45 | & Ymed1,Ymed2,sgmY1,sgmY2,XYsect, | |
46 | & Gaus1,Gaus2,Gauss,Exp1,Exp2,Exp3, | |
47 | & XmXmin,XmXmax,XpXmin,XpXmax,Exmx1,Exmx2,Xmed, | |
48 | & Sgm1,AnorX,Icase | |
49 | ||
50 | real*8 mass | |
51 | ||
52 | data pi / 3.141 592 653 589 793 238 462 643d00 / | |
53 | ||
54 | C | |
55 | C Exact differential Cross section: | |
56 | C - - - - - - - - - - - - - - - - - - - | |
57 | C | |
58 | gm = 2750.0d0 ! Pb gamma factor at the LHC | |
59 | mass = 0.5109991d0 ! electron mass | |
60 | call initdiffcross (gm,mass) | |
61 | C | |
62 | C Cross sections: | |
63 | C - - - - - - - - - - | |
64 | Nevnt = 0 | |
65 | Xsect2 = 0. | |
66 | Dsect2 = 0. | |
67 | ||
68 | Nsd = 1 | |
69 | Npt = 25 ! 7,25,64 | |
70 | Eps = 0.00005 | |
71 | ||
72 | Xmin= Dlog10(Ptmin) | |
73 | Xmax= Dlog10(Ptmax) | |
74 | write(*,*) ' Kinematical limits:' | |
75 | write(*,*) ' Xmin,Xmax=',Xmin,Xmax | |
76 | write(*,*) ' Ymin,Ymax=',Ymin,Ymax | |
77 | write(*,*) | |
78 | ||
79 | Xsec1 = Dtrint(DsdXX,Nsd,Npt,Eps,Xmin,Xmin,Xmin,Xmax,Xmax,Xmin) | |
80 | Xsec2 = Dtrint(DsdXX,Nsd,Npt,Eps,Xmax,Xmax,Xmin,Xmax,Xmax,Xmin) | |
81 | XsecX = Xsec1+Xsec2 | |
82 | ||
83 | Ysec1 = Dtrint(DsdYY,Nsd,Npt,Eps,Ymin,Ymin,Ymin,Ymax,Ymax,Ymin) | |
84 | Ysec2 = Dtrint(DsdYY,Nsd,Npt,Eps,Ymax,Ymax,Ymin,Ymax,Ymax,Ymin) | |
85 | XsecY = Ysec1+Ysec2 | |
86 | ||
87 | IF (Xsec1*Xsec2.le.0.or.Ysec1*Ysec2.le.0.) then | |
88 | write(*,*) ' Error: insufficient accuracy of XY-cross sections' | |
89 | write(*,*) ' Xsec1,Xsec2,Xsec=', Xsec1,Xsec2,XsecX | |
90 | write(*,*) ' Ysec1,Ysec2,Ysec=', Ysec1,Ysec2,XsecY | |
91 | endif | |
92 | ||
93 | XYsect = XsecX*XsecY | |
94 | write(*,*) ' Xsections: Xsec1,Xsec2,XsecX=', Xsec1,Xsec2,XsecX | |
95 | write(*,*) ' Ysec1,Ysec2,XsecY=', Ysec1,Ysec2,XsecY | |
96 | C- write(*,*) ' XsecX*YsecY=', XYsect | |
97 | C | |
98 | XYsect = 2371.5239*(82./137.035)**4*XYsect ! Normalization factor | |
99 | write(*,*) ' Normalized: XsecX*YsecY=', XYsect | |
100 | write(*,*) | |
101 | C | |
102 | C- XsecPhi = Dgauss(DsdPhi,0.0d0, pi, eps) | |
103 | C | |
104 | C Rapidity initialisation: | |
105 | C - - - - - - - - - - - - - - | |
106 | if (Ymin.ge.Ymax) then | |
107 | write(*,*) 'Wrong values of Ymin,Ymax:',Ymin,Ymax | |
108 | stop | |
109 | endif | |
110 | CYpY- | |
111 | C | |
112 | YpYmin = 2.*Ymin | |
113 | YpYmax = 2.*Ymax | |
114 | Yp0 = 0. | |
115 | if (YpYmin*YpYmax.gt.0.) then | |
116 | if(YpYmin.gt.0.) Yp0 = YpYmin | |
117 | if(YpYmin.le.0.) Yp0 = YpYmax | |
118 | endif | |
119 | WYpYmax = DsdYpY(Yp0) | |
120 | C | |
121 | C- write(*,*)' YpY: YpYmin,YpYmax,WYpYmax =',YpYmin,YpYmax,WYpYmax | |
122 | ||
123 | CYmY- | |
124 | C | |
125 | YmYmin = Ymin-Ymax | |
126 | YmYmax = Ymax-Ymin | |
127 | C- write(*,*) ' YmY: YmYmin,YmYmax=',YmYmin,YmYmax | |
128 | C | |
129 | Ymed1 = 0.18 | |
130 | Ymed2 = 4.00 | |
131 | ||
132 | sgmY1 = 1./dsqrt(2.*parYmY(2)) | |
133 | sgmY2 = 1./dsqrt(2.*parYmY(4)) | |
134 | C | |
135 | eps = 0.000001 | |
136 | Gaus1 = Dgauss(DsdYmY ,0.0d0, Ymed1 , eps) | |
137 | Gaus2 = Dgauss(DsdYmY ,Ymed1, Ymed2 , eps) | |
138 | Expon = Dgauss(DsdYmY ,Ymed2, YmYmax, eps) | |
139 | Summ = Gaus1+Gaus2+Expon | |
140 | Gaus2 =(Gaus1+Gaus2)/Summ | |
141 | Gaus1 = Gaus1/Summ | |
142 | Expon = 1. | |
143 | C | |
144 | C Phi initialisation: | |
145 | C - - - - - - - - - - - | |
146 | Exp0 = parPhi(1)*pi | |
147 | Exp1 = parPhi(2)*(1.-dexp(-parPhi(3)*pi))/parPhi(3) | |
148 | Exp2 = parPhi(4)*(1.-dexp(-parPhi(5)*pi))/parPhi(5) | |
149 | Exp3 = parPhi(6)*(1.-dexp(-parPhi(7)*pi))/parPhi(7) | |
150 | Summ = Exp0+Exp1+Exp2+Exp3 | |
151 | C | |
152 | Exp0 = (Exp0+Exp1+Exp2+Exp3)/Summ | |
153 | Exp1 = (Exp1+Exp2+Exp3)/Summ | |
154 | Exp2 = (Exp2+Exp3)/Summ | |
155 | Exp3 = Exp3/Summ | |
156 | C- write(*,*) 'Exp0,Exp1,Exp2,Exp3=',Exp0,Exp1,Exp2,Exp3 | |
157 | C | |
158 | C Pt initialisation: | |
159 | C - - - - - - - - - - - | |
160 | C | |
161 | XmXmin = Xmin-Xmax | |
162 | XmXmax = Xmax-Xmin | |
163 | C | |
164 | XpXmin = 2.*Xmin | |
165 | XpXmax = 2.*Xmax | |
166 | C | |
167 | CXmX- | |
168 | Exmx1 = parXmX(1)*(1.-dexp(-parXmX(2)*dabs(XmXmax)))/parXmX(2) | |
169 | Exmx2 = parXmX(3)*(1.-dexp(-parXmX(4)*dabs(XmXmax)))/parXmX(4) | |
170 | Exmx1 = Exmx1/(Exmx1+Exmx2) | |
171 | Exmx2 = 1. | |
172 | C | |
173 | C- write(*,*) ' XpX: XpXmin,XpXmax=',XpXmin,XpXmax | |
174 | C- write(*,*) ' XmX: XmXmin,XmXmax=',XmXmin,XmXmax | |
175 | C- write(*,*) ' XmX: Exmx1 ,Exmx2 =',Exmx1 ,Exmx2 | |
176 | C | |
177 | CXpX- | |
178 | Icase = 2 ! Gausses and Exponent | |
179 | Xmed = 0.6 | |
180 | if (XpXmax.lt.Xmed) then | |
181 | Icase = 1 ! Gauss only | |
182 | Xmed = XpXmax | |
183 | endif | |
184 | C | |
185 | if (XpXmin.gt.Xmed) then | |
186 | Icase = 3 ! Exponent only | |
187 | Xmed = XpXmin | |
188 | endif | |
189 | C | |
190 | if (Icase.eq.2) then ! Gauss and Exponent | |
191 | eps = 0.000001 | |
192 | Gauss = Dgauss(DsdXpX ,XpXmin, Xmed , eps) | |
193 | Expon = Dgauss(DsdXpX ,Xmed , XpXmax, eps) | |
194 | Summ = Gauss+Expon | |
195 | Gauss = Gauss/Summ | |
196 | Expon = Expon/Summ | |
197 | C- write(*,*) 'XpX, Case 2: Gauss,Expon=',Gauss,Expon | |
198 | endif ! Icase = 2 | |
199 | C | |
200 | C | |
201 | if (Icase.eq.1.or.Icase.eq.2) then ! Gausses only and Icase=2 | |
202 | ||
203 | Sgm1 = 1./dsqrt(2.*parXpX(2)) | |
204 | ||
205 | C- write(*,*) 'XpX, Case 1: Sgm1=',Sgm1 | |
206 | endif | |
207 | C | |
208 | if (Icase.eq.3.or.Icase.eq.2) then ! Exponent only and Icase=2 | |
209 | C | |
210 | AnorX = 1.-dexp(-parXpX(5)*(XpXmax-Xmed)) | |
211 | C | |
212 | C- write(*,*) 'XpX, Case 3: AnorX=',AnorX | |
213 | endif | |
214 | ||
215 | return | |
216 | end | |
217 | * | |
218 | *======================================================================= | |
219 | subroutine ee_event(Ymin,Ymax,PTmin,PTmax, | |
220 | + Ye,Yp,Xe,Xp,Phi,Wtm2) | |
221 | C------------------------------------------------------------------------------ | |
222 | C Produce one event | |
223 | C | |
224 | C Input variables: Ymin - minimal value of rapidity \ | |
225 | C Ymax - maximal value of rapidity | of kinematics | |
226 | C PTmin - Pt minimum in MeV/c | range | |
227 | C PTmax - Pt maximum in MeV/c / | |
228 | C Output variables: Ye,Yp - rapidity of produced e- or e+ | |
229 | C Xe,Xp - log10(pt) of produced e- or e+, pt in MeV/c | |
230 | C Phi - azymuth angle between e- and e+ | |
231 | C Wtm2 - event weight. The sum of these event | |
232 | C weights, divided by the total number | |
233 | C of generated events, gives the integral | |
234 | C cross section of the process of e+e- pair | |
235 | C production in the above mentioned kinematics | |
236 | C range. | |
237 | C Sum of the selected event weights, divided | |
238 | C by the total number of generated events, | |
239 | C gives the integral cross section corresponded | |
240 | C to the set of selected events | |
241 | C------------------------------------------------------------------------------ | |
242 | ||
243 | implicit real*8 (A-H,O-Z) | |
244 | common/parPh / parPhi(7) | |
245 | common/parYp / parYpY(6) | |
246 | common/parYm / parYmY(6) | |
247 | common/parXp / parXpX(5) | |
248 | common/parXm / parXmX(4) | |
249 | common/eevent/ Xsect2,Dsect2, Xsecttot,Dsecttot, Nevnt | |
250 | common/eepars/ Xmin,Xmax,YpYmin,YpYmax,WYpYmax,YmYmin,YmYmax, | |
251 | & Ymed1,Ymed2,sgmY1,sgmY2,XYsect, | |
252 | & Gaus1,Gaus2,Gauss,Exp1,Exp2,Exp3, | |
253 | & XmXmin,XmXmax,XpXmin,XpXmax,Exmx1,Exmx2,Xmed, | |
254 | & Sgm1,AnorX,Icase | |
255 | external DsdYpY,DsdYmY,DsdXpX,DsdXmX,DsdPhi | |
256 | data pi / 3.141 592 653 589 793 238 462 643d00 / | |
257 | ||
258 | C Rapidity distributions: | |
259 | C | |
63e30ba9 | 260 | 10 YpY = YpYmin + (YpYmax-YpYmin)*eernd(0) |
c9c6bd3b | 261 | Wt = DsdYpY(YpY) |
63e30ba9 | 262 | if (Wt.lt.eernd(0)*WYpYmax) go to 10 |
c9c6bd3b | 263 | C |
264 | CYmY- | |
63e30ba9 | 265 | r1 = eernd(0) |
c9c6bd3b | 266 | if (r1.lt.Gaus1) then |
267 | 13 call rnorml(rn1) | |
268 | YmY = sgmY1*rn1 | |
269 | if (dabs(YmY).gt.Ymed1) go to 13 | |
270 | else if (r1.lt.Gaus2) then | |
271 | 15 call rnorml(rn2) | |
272 | YmY = sgmY2*rn2 | |
273 | if (dabs(YmY).lt.Ymed1) go to 15 | |
274 | if (dabs(YmY).gt.Ymed2) go to 15 | |
275 | else | |
63e30ba9 | 276 | 20 r2 = eernd(0) |
c9c6bd3b | 277 | YmY=-Dlog(r2)/parYmY(6) + Ymed2 |
278 | if (YmY.gt.YmYmax) go to 20 | |
63e30ba9 | 279 | r2 = eernd(0) |
c9c6bd3b | 280 | if (r2.lt.0.5) YmY =-YmY |
281 | endif | |
282 | C | |
283 | Ye = 0.5*(YpY+YmY) | |
284 | Yp = 0.5*(YpY-YmY) | |
285 | ||
286 | if (Ye.lt.Ymin.or.Ye.gt.Ymax) go to 10 | |
287 | if (Yp.lt.Ymin.or.Yp.gt.Ymax) go to 10 | |
288 | ||
289 | C | |
290 | C Azimuthal angle: | |
291 | C | |
63e30ba9 | 292 | 30 r1 = eernd(0) |
c9c6bd3b | 293 | if (r1.lt.Exp3) then |
63e30ba9 | 294 | 31 r2 = eernd(0) |
c9c6bd3b | 295 | Phi=-Dlog(r2)/parPhi(7) |
296 | if (Phi.gt.pi) go to 31 | |
297 | go to 50 | |
298 | else if (r1.lt.Exp2) then | |
63e30ba9 | 299 | 32 r2 = eernd(0) |
c9c6bd3b | 300 | Phi=-Dlog(r2)/parPhi(5) |
301 | if (Phi.gt.pi) go to 32 | |
302 | go to 50 | |
303 | else if (r1.lt.Exp1) then | |
63e30ba9 | 304 | 33 r2 = eernd(0) |
c9c6bd3b | 305 | Phi=-Dlog(r2)/parPhi(3) |
306 | if (Phi.gt.pi) go to 33 | |
307 | go to 50 | |
308 | else | |
63e30ba9 | 309 | Phi = pi*eernd(0) |
c9c6bd3b | 310 | endif |
311 | C | |
63e30ba9 | 312 | 50 if (eernd(0).gt.0.5) Phi =-Phi |
c9c6bd3b | 313 | Phi = Phi+pi |
314 | if (Phi.lt.0.03.or.Phi.gt.2.*pi-0.03) go to 30 | |
315 | C | |
316 | C Transverse momentums: | |
317 | C | |
318 | CXpX- | |
319 | ||
320 | 60 Jcase = Icase ! Jcase = Icase, if Icase.ne.2 | |
321 | C | |
322 | if (Icase.eq.2) then ! Gausses and Exponent | |
323 | Jcase = 3 | |
63e30ba9 | 324 | if (eernd(0).lt.Gauss) Jcase = 1 |
c9c6bd3b | 325 | endif ! of Icase 2 |
326 | C | |
327 | if (Jcase.eq.1) then ! Gauss only | |
328 | 65 call rnorml(rn1) | |
329 | XpX=Sgm1*rn1-parXpX(3) | |
330 | if (XpX.lt.XpXmin.or.XpX.gt.Xmed) go to 65 | |
331 | endif ! of Jcase 1 | |
332 | C | |
333 | if (Jcase.eq.3) then ! Exponent only | |
63e30ba9 | 334 | 70 r1 = eernd(0) |
c9c6bd3b | 335 | XpX =-Dlog( 1.- r1*AnorX)/parXpX(5)+Xmed |
336 | endif ! of Jcase 3 | |
337 | ||
338 | CXmX- | |
63e30ba9 | 339 | r1 = eernd(0) |
c9c6bd3b | 340 | if (r1.lt.Exmx1) then |
63e30ba9 | 341 | 81 r2 = eernd(0) |
c9c6bd3b | 342 | XmX=-Dlog(r2)/parXmX(2) |
343 | if (XmX.gt.dabs(XmXmax)) go to 81 | |
344 | else | |
63e30ba9 | 345 | 82 r2 = eernd(0) |
c9c6bd3b | 346 | XmX=-Dlog(r2)/parXmX(4) |
347 | if (XmX.gt.dabs(XmXmax)) go to 82 | |
348 | endif | |
349 | C | |
63e30ba9 | 350 | 85 if (eernd(0).gt.0.5) XmX =-XmX |
c9c6bd3b | 351 | C |
352 | Xe = 0.5*(XpX+XmX) | |
353 | Xp = 0.5*(XpX-XmX) | |
354 | ||
355 | if (Xe.lt.Xmin.or.Xe.gt.Xmax) go to 60 | |
356 | if (Xp.lt.Xmin.or.Xp.gt.Xmax) go to 60 | |
357 | C | |
358 | C --- Good event: | |
359 | C | |
360 | Pte = 10.d0**Xe | |
361 | Ptp = 10.d0**Xp | |
362 | C | |
363 | C Event weight if the events would generate uniformly | |
364 | C | |
365 | Wt = DsdYpY(YpY)*DsdYmY(YmY)*DsdXpX(Xpx)*DsdXmX(XmX)*DsdPhi(Phi) | |
366 | Wt = 2.2706950/Wt | |
367 | C | |
368 | Wtm2=(Pte*Ptp)*Wt ! Transformation factor | |
369 | C | |
370 | C --- Exact diff. cross section (Adrian Alscher 1997, Kai Hencken, May '98): | |
371 | C | |
372 | call Diffcross(Ptp,Yp,Pte,Ye,Phi,Dsigma) | |
373 | Wtm2= XYsect*Dsigma*(Pte*Ptp)*Wtm2 | |
374 | C | |
375 | Nevnt = Nevnt + 1 | |
376 | Xsect2 = Xsect2 + Wtm2 | |
377 | Dsect2 = Dsect2 + Wtm2*Wtm2 | |
378 | C | |
379 | C Calculate cross section and its error accumulated so far | |
380 | C | |
381 | Xsecttot = Xsect2/Nevnt | |
382 | Dsecttot = dsqrt(Dsect2)/Nevnt | |
383 | C | |
384 | return | |
385 | end | |
386 | c | |
387 | C================================================================= | |
388 | Double precision function DsdYpY(Y) | |
389 | c | |
390 | c Y = Yp+Ye | |
391 | c - - - - - - - - - - - - - - - - - - - | |
392 | implicit real*8 (A-H,O-Z) | |
393 | common/parYp/ parYpY(6) | |
394 | C- data parYpY / -48.584, 0.11403E-01, 40.649, 0.38920E-04, | |
395 | C- + 40.283, 0.19238E-01 / | |
396 | data parYpY / -4.8584, 0.11403E-01, 4.0649, 0.38920E-04, | |
397 | + 4.0283, 0.19238E-01 / | |
398 | c | |
399 | DsdYpY = parYpY(1)*dexp(-parYpY(2)*Y**2) | |
400 | + + parYpY(3)*dexp(-parYpY(4)*Y**4) | |
401 | + + parYpY(5)*dexp(-parYpY(6)*Y**2) | |
402 | return | |
403 | end | |
404 | C | |
405 | Double precision function DsdYmY(Y) | |
406 | c | |
407 | c Y = Yp-Ye | |
408 | c - - - - - - - - - - - - - - - - - - - | |
409 | implicit real*8 (A-H,O-Z) | |
410 | common/parYm/ parYmY(6) | |
411 | C- data parYmY / 180., 8., 174.38, 0.17867, 2418.8, 1.3097 / | |
412 | data parYmY / 1.80, 8., 1.7438, 0.17867, 24.188, 1.3097 / | |
413 | C | |
414 | if (abs(Y).lt.0.18) then | |
415 | C1- if (abs(Y).lt.0.00) then | |
416 | DsdYmY = parYmY(1)*dexp(-parYmY(2)*abs(Y)**2) | |
417 | else if (abs(Y).lt.4.00) then | |
418 | DsdYmY = parYmY(3)*dexp(-parYmY(4)*abs(Y)**2) | |
419 | else | |
420 | DsdYmY = parYmY(5)*dexp(-parYmY(6)*abs(Y)) | |
421 | endif | |
422 | ||
423 | return | |
424 | end | |
425 | C | |
426 | Double precision function DsdYY(Yp,Ye) | |
427 | c - - - - - - - - - - - - - - - - - - - - - - | |
428 | implicit real*8 (A-H,O-Z) | |
429 | c | |
430 | YpY = Yp+Ye | |
431 | YmY = Yp-Ye | |
432 | DsdYY = DsdYpY(YpY)*DsdYmY(YmY) | |
433 | return | |
434 | end | |
435 | C | |
436 | C==================================================================== | |
437 | C | |
438 | Double precision function DsdXpX(X) | |
439 | c | |
440 | c X = Xp+Xe | |
441 | c - - - - - - - - - - - - - - - - - - - | |
442 | implicit real*8 (A-H,O-Z) | |
443 | common/parXp/ parXpX(5) | |
444 | ||
445 | C- data parXpX / 83.668, 1.2004, 0.47225, 96.951, 2.3814 / | |
446 | data parXpX / 8.3668, 1.2004, 0.47225, 9.6951, 2.3814 / | |
447 | c | |
448 | if (X.lt.0.6)then | |
449 | DsdXpX = parXpX(1)*dexp(-parXpX(2)*(X+parXpX(3))**2) | |
450 | else | |
451 | DsdXpX = parXpX(4)*dexp(-parXpX(5)*X) | |
452 | endif | |
453 | return | |
454 | end | |
455 | C | |
456 | Double precision function DsdXmX(X) | |
457 | c | |
458 | c X = Xp-Xe | |
459 | c - - - - - - - - - - - - - - - - - - - | |
460 | implicit real*8 (A-H,O-Z) | |
461 | common/parXm/ parXmX(4) | |
462 | c- data parXmX / 950.38, 39.040, 194.92, 3.5660 / | |
463 | data parXmX / 9.5038, 39.040, 1.9492, 3.5660 / | |
464 | c | |
465 | ||
466 | DsdXmX = parXmX(1)*dexp(-parXmX(2)*abs(X)) | |
467 | + + parXmX(3)*dexp(-parXmX(4)*abs(X)) | |
468 | ||
469 | return | |
470 | end | |
471 | C | |
472 | Double precision function DsdXX(Xp,Xe) | |
473 | c - - - - - - - - - - - - - - - - - - - - - - | |
474 | implicit real*8 (A-H,O-Z) | |
475 | c | |
476 | XpX = Xp+Xe | |
477 | XmX = Xp-Xe | |
478 | DsdXX = DsdXpX(XpX)*DsdXmX(XmX) | |
479 | return | |
480 | end | |
481 | C==================================================================== | |
482 | C | |
483 | Double precision function DsdPhi(Phi) ! Phi distribution | |
484 | c | |
485 | c Phi - athimutal angle between e+ and e- (in radians) | |
486 | c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | |
487 | implicit real*8 (A-H,O-Z) | |
488 | common/parPh/ parPhi(7) | |
489 | C ! For 10**7 events | |
490 | data parPhi / 5.4825, 226.00, 11.602, 68.173, ! accuracy= 0.87 % | |
491 | + 1.9509, 0.11864E+04, 109.61 / ! N(Wt>20)= 573 | |
492 | ||
493 | C | |
494 | data pi / 3.141 592 653 589 793 238 462 643d00/ | |
495 | C | |
496 | DsdPhi = parPhi(1) | |
497 | + + parPhi(2)*dexp(-parPhi(3)*dabs(Phi-pi)) | |
498 | + + parPhi(4)*dexp(-parPhi(5)*dabs(Phi-pi)) | |
499 | + + parPhi(6)*dexp(-parPhi(7)*dabs(Phi-pi)) | |
500 | return | |
501 | end | |
502 | C==================================================================== | |
503 | SUBROUTINE rnorml(rnd) | |
504 | * Random generator of normal distribution | |
505 | implicit real*8 (A-H,O-Z) | |
506 | PARAMETER (pi=3.141 592 653 589 793 238 462 643d00) | |
507 | PARAMETER (pi2=2.*pi) | |
508 | C | |
63e30ba9 | 509 | u1 = eernd(0) |
510 | u2 = eernd(0) | |
c9c6bd3b | 511 | p1 = pi2*u1 |
512 | p2 = dsqrt(-2.*dlog(u2)) | |
63e30ba9 | 513 | rnd= dcos(p1)*p2 |
c9c6bd3b | 514 | RETURN |
515 | END |