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
269 if (dabs(YmY).gt.Ymed1) go to 13
270 else if (r1.lt.Gaus2) then
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
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 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)
512 p2 = dsqrt(-2.*dlog(u2))