Major update of the CMake compilation:
[u/mrichter/AliRoot.git] / TEPEMGEN / epemgen.f
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
260  10   YpY = YpYmin + (YpYmax-YpYmin)*eernd(0)
261       Wt  = DsdYpY(YpY)
262       if (Wt.lt.eernd(0)*WYpYmax) go to 10
263 C
264 CYmY-
265       r1 = eernd(0)
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
276  20      r2 = eernd(0)
277          YmY=-Dlog(r2)/parYmY(6) + Ymed2
278          if (YmY.gt.YmYmax) go to 20    
279          r2 = eernd(0)
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
292  30   r1 = eernd(0)
293       if (r1.lt.Exp3) then
294  31      r2 = eernd(0)
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
299  32      r2 = eernd(0)
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 
304  33      r2 = eernd(0)
305          Phi=-Dlog(r2)/parPhi(3)
306          if (Phi.gt.pi) go to 33
307          go to 50
308       else
309          Phi = pi*eernd(0)
310       endif
311 C     
312  50   if (eernd(0).gt.0.5) Phi =-Phi
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
324          if (eernd(0).lt.Gauss) Jcase = 1
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
334  70      r1  = eernd(0)
335          XpX =-Dlog( 1.- r1*AnorX)/parXpX(5)+Xmed             
336       endif                     ! of Jcase 3 
337
338 CXmX-
339       r1 = eernd(0)
340       if (r1.lt.Exmx1) then
341  81      r2 = eernd(0)
342          XmX=-Dlog(r2)/parXmX(2)
343          if (XmX.gt.dabs(XmXmax)) go to 81           
344       else 
345  82      r2 = eernd(0)
346          XmX=-Dlog(r2)/parXmX(4)
347          if (XmX.gt.dabs(XmXmax)) go to 82  
348       endif
349 C
350  85   if (eernd(0).gt.0.5) XmX =-XmX
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
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
436 C====================================================================
437
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
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      
509       u1 = eernd(0)
510       u2 = eernd(0)
511       p1 = pi2*u1
512       p2 = dsqrt(-2.*dlog(u2))
513       rnd= dcos(p1)*p2
514       RETURN
515       END