1 C **********************************************************************
3 C Generator of e+e- pairs produced in 1.10.2002
4 C PbPb collisions at LHC
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>
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
21 C Long write up: K.Hencken,Yu.Kharlov,S.Sadovsky
22 C Internal ALICE Note 2002-27
24 C **********************************************************************
25 subroutine ee_init(Ymin,Ymax,PTmin,PTmax)
26 C-----------------------------------------------------------------------
27 C Generator initialisation
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
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)
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,
52 data pi / 3.141 592 653 589 793 238 462 643d00 /
55 C Exact differential Cross section:
56 C - - - - - - - - - - - - - - - - - - -
58 gm = 2750.0d0 ! Pb gamma factor at the LHC
59 mass = 0.5109991d0 ! electron mass
60 call initdiffcross (gm,mass)
74 write(*,*) ' Kinematical limits:'
75 write(*,*) ' Xmin,Xmax=',Xmin,Xmax
76 write(*,*) ' Ymin,Ymax=',Ymin,Ymax
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)
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)
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
94 write(*,*) ' Xsections: Xsec1,Xsec2,XsecX=', Xsec1,Xsec2,XsecX
95 write(*,*) ' Ysec1,Ysec2,XsecY=', Ysec1,Ysec2,XsecY
96 C- write(*,*) ' XsecX*YsecY=', XYsect
98 XYsect = 2371.5239*(82./137.035)**4*XYsect ! Normalization factor
99 write(*,*) ' Normalized: XsecX*YsecY=', XYsect
102 C- XsecPhi = Dgauss(DsdPhi,0.0d0, pi, eps)
104 C Rapidity initialisation:
105 C - - - - - - - - - - - - - -
106 if (Ymin.ge.Ymax) then
107 write(*,*) 'Wrong values of Ymin,Ymax:',Ymin,Ymax
115 if (YpYmin*YpYmax.gt.0.) then
116 if(YpYmin.gt.0.) Yp0 = YpYmin
117 if(YpYmin.le.0.) Yp0 = YpYmax
119 WYpYmax = DsdYpY(Yp0)
121 C- write(*,*)' YpY: YpYmin,YpYmax,WYpYmax =',YpYmin,YpYmax,WYpYmax
127 C- write(*,*) ' YmY: YmYmin,YmYmax=',YmYmin,YmYmax
132 sgmY1 = 1./dsqrt(2.*parYmY(2))
133 sgmY2 = 1./dsqrt(2.*parYmY(4))
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
144 C Phi initialisation:
145 C - - - - - - - - - - -
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
152 Exp0 = (Exp0+Exp1+Exp2+Exp3)/Summ
153 Exp1 = (Exp1+Exp2+Exp3)/Summ
154 Exp2 = (Exp2+Exp3)/Summ
156 C- write(*,*) 'Exp0,Exp1,Exp2,Exp3=',Exp0,Exp1,Exp2,Exp3
159 C - - - - - - - - - - -
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)
173 C- write(*,*) ' XpX: XpXmin,XpXmax=',XpXmin,XpXmax
174 C- write(*,*) ' XmX: XmXmin,XmXmax=',XmXmin,XmXmax
175 C- write(*,*) ' XmX: Exmx1 ,Exmx2 =',Exmx1 ,Exmx2
178 Icase = 2 ! Gausses and Exponent
180 if (XpXmax.lt.Xmed) then
181 Icase = 1 ! Gauss only
185 if (XpXmin.gt.Xmed) then
186 Icase = 3 ! Exponent only
190 if (Icase.eq.2) then ! Gauss and Exponent
192 Gauss = Dgauss(DsdXpX ,XpXmin, Xmed , eps)
193 Expon = Dgauss(DsdXpX ,Xmed , XpXmax, eps)
197 C- write(*,*) 'XpX, Case 2: Gauss,Expon=',Gauss,Expon
201 if (Icase.eq.1.or.Icase.eq.2) then ! Gausses only and Icase=2
203 Sgm1 = 1./dsqrt(2.*parXpX(2))
205 C- write(*,*) 'XpX, Case 1: Sgm1=',Sgm1
208 if (Icase.eq.3.or.Icase.eq.2) then ! Exponent only and Icase=2
210 AnorX = 1.-dexp(-parXpX(5)*(XpXmax-Xmed))
212 C- write(*,*) 'XpX, Case 3: AnorX=',AnorX
218 *=======================================================================
219 subroutine ee_event(Ymin,Ymax,PTmin,PTmax,
220 + Ye,Yp,Xe,Xp,Phi,Wtm2)
221 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
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------------------------------------------------------------------------------
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,
255 external DsdYpY,DsdYmY,DsdXpX,DsdXmX,DsdPhi
256 data pi / 3.141 592 653 589 793 238 462 643d00 /
258 C Rapidity distributions:
260 10 YpY = YpYmin + (YpYmax-YpYmin)*eernd(0)
262 if (Wt.lt.eernd(0)*WYpYmax) go to 10
266 if (r1.lt.Gaus1) then
267 13 call rnormlEP(rn1)
269 if (dabs(YmY).gt.Ymed1) go to 13
270 else if (r1.lt.Gaus2) then
271 15 call rnormlEP(rn2)
273 if (dabs(YmY).lt.Ymed1) go to 15
274 if (dabs(YmY).gt.Ymed2) go to 15
277 YmY=-Dlog(r2)/parYmY(6) + Ymed2
278 if (YmY.gt.YmYmax) go to 20
280 if (r2.lt.0.5) YmY =-YmY
286 if (Ye.lt.Ymin.or.Ye.gt.Ymax) go to 10
287 if (Yp.lt.Ymin.or.Yp.gt.Ymax) go to 10
295 Phi=-Dlog(r2)/parPhi(7)
296 if (Phi.gt.pi) go to 31
298 else if (r1.lt.Exp2) then
300 Phi=-Dlog(r2)/parPhi(5)
301 if (Phi.gt.pi) go to 32
303 else if (r1.lt.Exp1) then
305 Phi=-Dlog(r2)/parPhi(3)
306 if (Phi.gt.pi) go to 33
312 50 if (eernd(0).gt.0.5) Phi =-Phi
314 if (Phi.lt.0.03.or.Phi.gt.2.*pi-0.03) go to 30
316 C Transverse momentums:
320 60 Jcase = Icase ! Jcase = Icase, if Icase.ne.2
322 if (Icase.eq.2) then ! Gausses and Exponent
324 if (eernd(0).lt.Gauss) Jcase = 1
327 if (Jcase.eq.1) then ! Gauss only
328 65 call rnormlEP(rn1)
329 XpX=Sgm1*rn1-parXpX(3)
330 if (XpX.lt.XpXmin.or.XpX.gt.Xmed) go to 65
333 if (Jcase.eq.3) then ! Exponent only
335 XpX =-Dlog( 1.- r1*AnorX)/parXpX(5)+Xmed
340 if (r1.lt.Exmx1) then
342 XmX=-Dlog(r2)/parXmX(2)
343 if (XmX.gt.dabs(XmXmax)) go to 81
346 XmX=-Dlog(r2)/parXmX(4)
347 if (XmX.gt.dabs(XmXmax)) go to 82
350 85 if (eernd(0).gt.0.5) XmX =-XmX
355 if (Xe.lt.Xmin.or.Xe.gt.Xmax) go to 60
356 if (Xp.lt.Xmin.or.Xp.gt.Xmax) go to 60
363 C Event weight if the events would generate uniformly
365 Wt = DsdYpY(YpY)*DsdYmY(YmY)*DsdXpX(Xpx)*DsdXmX(XmX)*DsdPhi(Phi)
368 Wtm2=(Pte*Ptp)*Wt ! Transformation factor
370 C --- Exact diff. cross section (Adrian Alscher 1997, Kai Hencken, May '98):
372 call Diffcross(Ptp,Yp,Pte,Ye,Phi,Dsigma)
373 Wtm2= XYsect*Dsigma*(Pte*Ptp)*Wtm2
376 Xsect2 = Xsect2 + Wtm2
377 Dsect2 = Dsect2 + Wtm2*Wtm2
379 C Calculate cross section and its error accumulated so far
381 Xsecttot = Xsect2/Nevnt
382 Dsecttot = dsqrt(Dsect2)/Nevnt
387 C=================================================================
388 Double precision function DsdYpY(Y)
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 /
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)
405 Double precision function DsdYmY(Y)
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 /
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)
420 DsdYmY = parYmY(5)*dexp(-parYmY(6)*abs(Y))
426 Double precision function DsdYY(Yp,Ye)
427 c - - - - - - - - - - - - - - - - - - - - - -
428 implicit real*8 (A-H,O-Z)
432 DsdYY = DsdYpY(YpY)*DsdYmY(YmY)
436 C====================================================================
438 Double precision function DsdXpX(X)
441 c - - - - - - - - - - - - - - - - - - -
442 implicit real*8 (A-H,O-Z)
443 common/parXp/ parXpX(5)
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 /
449 DsdXpX = parXpX(1)*dexp(-parXpX(2)*(X+parXpX(3))**2)
451 DsdXpX = parXpX(4)*dexp(-parXpX(5)*X)
456 Double precision function DsdXmX(X)
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 /
466 DsdXmX = parXmX(1)*dexp(-parXmX(2)*abs(X))
467 + + parXmX(3)*dexp(-parXmX(4)*abs(X))
472 Double precision function DsdXX(Xp,Xe)
473 c - - - - - - - - - - - - - - - - - - - - - -
474 implicit real*8 (A-H,O-Z)
478 DsdXX = DsdXpX(XpX)*DsdXmX(XmX)
481 C====================================================================
483 Double precision function DsdPhi(Phi) ! Phi distribution
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)
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
494 data pi / 3.141 592 653 589 793 238 462 643d00/
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))
502 C====================================================================
503 SUBROUTINE rnormlEP(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)
512 p2 = dsqrt(-2.*dlog(u2))