CMake: Retrieve Git information
[u/mrichter/AliRoot.git] / TEPEMGEN / epemgen.f
CommitLineData
c9c6bd3b 1C **********************************************************************
2C Version 2.0
3C Generator of e+e- pairs produced in 1.10.2002
4C PbPb collisions at LHC
5C
6C Copyright (c) 2002
7C Authors: Kai Hencken <k.hencken@unibas.ch>
8C Yuri Kharlov <Yuri.Kharlov@cern.ch>
9C Serguei Sadovsky <sadovsky@mx.ihep.su>,
10C <Serguei.Sadovski@cern.ch>
11C
12C Permission to use, copy and distribute this software and its
13C documentation strictly for non-commercial purposes is hereby granted
14C without fee, provided that the above copyright notice appears in all
15C copies and that both the copyright notice and this permission notice
16C appear in the supporting documentation. The authors make no claims
17C about the suitability of this software for any purpose. It is
18C provided "as is" without express or implied warranty.
19C Any change of the code should be submitted to the authors
20C
21C Long write up: K.Hencken,Yu.Kharlov,S.Sadovsky
22C Internal ALICE Note 2002-27
23C
24C **********************************************************************
25 subroutine ee_init(Ymin,Ymax,PTmin,PTmax)
26C-----------------------------------------------------------------------
27C Generator initialisation
28C
29C Input variables: Ymin - minimal value of rapidity \
30C Ymax - maximal value of rapidity | of kinematics
31C PTmin - Pt minimum in MeV/c | range
32C PTmax - Pt maximum in MeV/c /
33C-----------------------------------------------------------------------
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
54C
55C Exact differential Cross section:
56C - - - - - - - - - - - - - - - - - - -
57C
58 gm = 2750.0d0 ! Pb gamma factor at the LHC
59 mass = 0.5109991d0 ! electron mass
60 call initdiffcross (gm,mass)
61C
62C Cross sections:
63C - - - - - - - - - -
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
96C- write(*,*) ' XsecX*YsecY=', XYsect
97C
98 XYsect = 2371.5239*(82./137.035)**4*XYsect ! Normalization factor
99 write(*,*) ' Normalized: XsecX*YsecY=', XYsect
100 write(*,*)
101C
102C- XsecPhi = Dgauss(DsdPhi,0.0d0, pi, eps)
103C
104C Rapidity initialisation:
105C - - - - - - - - - - - - - -
106 if (Ymin.ge.Ymax) then
107 write(*,*) 'Wrong values of Ymin,Ymax:',Ymin,Ymax
108 stop
109 endif
110CYpY-
111C
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)
120C
121C- write(*,*)' YpY: YpYmin,YpYmax,WYpYmax =',YpYmin,YpYmax,WYpYmax
122
123CYmY-
124C
125 YmYmin = Ymin-Ymax
126 YmYmax = Ymax-Ymin
127C- write(*,*) ' YmY: YmYmin,YmYmax=',YmYmin,YmYmax
128C
129 Ymed1 = 0.18
130 Ymed2 = 4.00
131
132 sgmY1 = 1./dsqrt(2.*parYmY(2))
133 sgmY2 = 1./dsqrt(2.*parYmY(4))
134C
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.
143C
144C Phi initialisation:
145C - - - - - - - - - - -
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
151C
152 Exp0 = (Exp0+Exp1+Exp2+Exp3)/Summ
153 Exp1 = (Exp1+Exp2+Exp3)/Summ
154 Exp2 = (Exp2+Exp3)/Summ
155 Exp3 = Exp3/Summ
156C- write(*,*) 'Exp0,Exp1,Exp2,Exp3=',Exp0,Exp1,Exp2,Exp3
157C
158C Pt initialisation:
159C - - - - - - - - - - -
160C
161 XmXmin = Xmin-Xmax
162 XmXmax = Xmax-Xmin
163C
164 XpXmin = 2.*Xmin
165 XpXmax = 2.*Xmax
166C
167CXmX-
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.
172C
173C- write(*,*) ' XpX: XpXmin,XpXmax=',XpXmin,XpXmax
174C- write(*,*) ' XmX: XmXmin,XmXmax=',XmXmin,XmXmax
175C- write(*,*) ' XmX: Exmx1 ,Exmx2 =',Exmx1 ,Exmx2
176C
177CXpX-
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
184C
185 if (XpXmin.gt.Xmed) then
186 Icase = 3 ! Exponent only
187 Xmed = XpXmin
188 endif
189C
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
197C- write(*,*) 'XpX, Case 2: Gauss,Expon=',Gauss,Expon
198 endif ! Icase = 2
199C
200C
201 if (Icase.eq.1.or.Icase.eq.2) then ! Gausses only and Icase=2
202
203 Sgm1 = 1./dsqrt(2.*parXpX(2))
204
205C- write(*,*) 'XpX, Case 1: Sgm1=',Sgm1
206 endif
207C
208 if (Icase.eq.3.or.Icase.eq.2) then ! Exponent only and Icase=2
209C
210 AnorX = 1.-dexp(-parXpX(5)*(XpXmax-Xmed))
211C
212C- 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)
221C------------------------------------------------------------------------------
222C Produce one event
223C
224C Input variables: Ymin - minimal value of rapidity \
225C Ymax - maximal value of rapidity | of kinematics
226C PTmin - Pt minimum in MeV/c | range
227C PTmax - Pt maximum in MeV/c /
228C Output variables: Ye,Yp - rapidity of produced e- or e+
229C Xe,Xp - log10(pt) of produced e- or e+, pt in MeV/c
230C Phi - azymuth angle between e- and e+
231C Wtm2 - event weight. The sum of these event
232C weights, divided by the total number
233C of generated events, gives the integral
234C cross section of the process of e+e- pair
235C production in the above mentioned kinematics
236C range.
237C Sum of the selected event weights, divided
238C by the total number of generated events,
239C gives the integral cross section corresponded
240C to the set of selected events
241C------------------------------------------------------------------------------
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
258C Rapidity distributions:
259C
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 263C
264CYmY-
63e30ba9 265 r1 = eernd(0)
c9c6bd3b 266 if (r1.lt.Gaus1) then
c9e3082d 267 13 call rnormlEP(rn1)
c9c6bd3b 268 YmY = sgmY1*rn1
269 if (dabs(YmY).gt.Ymed1) go to 13
270 else if (r1.lt.Gaus2) then
c9e3082d 271 15 call rnormlEP(rn2)
c9c6bd3b 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
282C
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
289C
290C Azimuthal angle:
291C
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
311C
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
315C
316C Transverse momentums:
317C
318CXpX-
319
320 60 Jcase = Icase ! Jcase = Icase, if Icase.ne.2
321C
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
326C
327 if (Jcase.eq.1) then ! Gauss only
c9e3082d 328 65 call rnormlEP(rn1)
c9c6bd3b 329 XpX=Sgm1*rn1-parXpX(3)
330 if (XpX.lt.XpXmin.or.XpX.gt.Xmed) go to 65
331 endif ! of Jcase 1
332C
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
338CXmX-
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
349C
63e30ba9 350 85 if (eernd(0).gt.0.5) XmX =-XmX
c9c6bd3b 351C
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
357C
358C --- Good event:
359C
360 Pte = 10.d0**Xe
361 Ptp = 10.d0**Xp
362C
363C Event weight if the events would generate uniformly
364C
365 Wt = DsdYpY(YpY)*DsdYmY(YmY)*DsdXpX(Xpx)*DsdXmX(XmX)*DsdPhi(Phi)
366 Wt = 2.2706950/Wt
367C
368 Wtm2=(Pte*Ptp)*Wt ! Transformation factor
369C
370C --- Exact diff. cross section (Adrian Alscher 1997, Kai Hencken, May '98):
371C
372 call Diffcross(Ptp,Yp,Pte,Ye,Phi,Dsigma)
373 Wtm2= XYsect*Dsigma*(Pte*Ptp)*Wtm2
374C
375 Nevnt = Nevnt + 1
376 Xsect2 = Xsect2 + Wtm2
377 Dsect2 = Dsect2 + Wtm2*Wtm2
378C
379C Calculate cross section and its error accumulated so far
380C
381 Xsecttot = Xsect2/Nevnt
382 Dsecttot = dsqrt(Dsect2)/Nevnt
383C
384 return
385 end
386c
387C=================================================================
388 Double precision function DsdYpY(Y)
389c
390c Y = Yp+Ye
391c - - - - - - - - - - - - - - - - - - -
392 implicit real*8 (A-H,O-Z)
393 common/parYp/ parYpY(6)
394C- data parYpY / -48.584, 0.11403E-01, 40.649, 0.38920E-04,
395C- + 40.283, 0.19238E-01 /
396 data parYpY / -4.8584, 0.11403E-01, 4.0649, 0.38920E-04,
397 + 4.0283, 0.19238E-01 /
398c
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
404C
405 Double precision function DsdYmY(Y)
406c
407c Y = Yp-Ye
408c - - - - - - - - - - - - - - - - - - -
409 implicit real*8 (A-H,O-Z)
410 common/parYm/ parYmY(6)
411C- 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 /
413C
414 if (abs(Y).lt.0.18) then
415C1- 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
425C
426 Double precision function DsdYY(Yp,Ye)
427c - - - - - - - - - - - - - - - - - - - - - -
428 implicit real*8 (A-H,O-Z)
429c
430 YpY = Yp+Ye
431 YmY = Yp-Ye
432 DsdYY = DsdYpY(YpY)*DsdYmY(YmY)
433 return
434 end
435C
436C====================================================================
437C
438 Double precision function DsdXpX(X)
439c
440c X = Xp+Xe
441c - - - - - - - - - - - - - - - - - - -
442 implicit real*8 (A-H,O-Z)
443 common/parXp/ parXpX(5)
444
445C- 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 /
447c
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
455C
456 Double precision function DsdXmX(X)
457c
458c X = Xp-Xe
459c - - - - - - - - - - - - - - - - - - -
460 implicit real*8 (A-H,O-Z)
461 common/parXm/ parXmX(4)
462c- data parXmX / 950.38, 39.040, 194.92, 3.5660 /
463 data parXmX / 9.5038, 39.040, 1.9492, 3.5660 /
464c
465
466 DsdXmX = parXmX(1)*dexp(-parXmX(2)*abs(X))
467 + + parXmX(3)*dexp(-parXmX(4)*abs(X))
468
469 return
470 end
471C
472 Double precision function DsdXX(Xp,Xe)
473c - - - - - - - - - - - - - - - - - - - - - -
474 implicit real*8 (A-H,O-Z)
475c
476 XpX = Xp+Xe
477 XmX = Xp-Xe
478 DsdXX = DsdXpX(XpX)*DsdXmX(XmX)
479 return
480 end
481C====================================================================
482C
483 Double precision function DsdPhi(Phi) ! Phi distribution
484c
485c Phi - athimutal angle between e+ and e- (in radians)
486c - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
487 implicit real*8 (A-H,O-Z)
488 common/parPh/ parPhi(7)
489C ! 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
493C
494 data pi / 3.141 592 653 589 793 238 462 643d00/
495C
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
502C====================================================================
c9e3082d 503 SUBROUTINE rnormlEP(rnd)
c9c6bd3b 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)
508C
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