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 rnormlEP(rn1)
 268          YmY = sgmY1*rn1
 269          if (dabs(YmY).gt.Ymed1) go to 13
 270       else if  (r1.lt.Gaus2) then
 271  15      call rnormlEP(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 C
 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 rnormlEP(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 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 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)
 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