Event generator for e+e- pair production (Yu.Kharlov)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 6 Nov 2002 10:15:32 +0000 (10:15 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 6 Nov 2002 10:15:32 +0000 (10:15 +0000)
EPEMGEN/diffcross.f [new file with mode: 0644]
EPEMGEN/dtrint.f [new file with mode: 0644]
EPEMGEN/epemgen.f [new file with mode: 0644]
EPEMGEN/libepemgen.pkg [new file with mode: 0644]
TEPEMGEN/TEPEMGENLinkDef.h [new file with mode: 0644]
TEPEMGEN/TEcommon.h [new file with mode: 0644]
TEPEMGEN/TEpEmGen.cxx [new file with mode: 0644]
TEPEMGEN/TEpEmGen.h [new file with mode: 0644]
TEPEMGEN/libTEPEMGEN.pkg [new file with mode: 0644]

diff --git a/EPEMGEN/diffcross.f b/EPEMGEN/diffcross.f
new file mode 100644 (file)
index 0000000..09213b9
--- /dev/null
@@ -0,0 +1,934 @@
+C**********************************************************************
+C*                                                                    *
+C* Exact calculation of the total differential e+ e- -pair production *
+C* in Relativistic Heavy Ion Collisions for a point particle in an    *
+C* external field approach.                                           *
+C*                                                                    *
+C* For details see the publication:                                   *
+C* "Multiple electromagnetic electron positron pair production in     *
+C*  relativistic heavy ion collisions"                                *
+C*  Adrian Alscher, Kai Hencken, Dirk Trautmann, and Gerhard Baur     *
+C*  Phys. Rev. A55 (1997) 396.                                        *
+C*                                                                    *
+C* Copyright (c) 1997, 1998, 2002, Adrian Alscher and Kai Hencken     *
+C*                                                                    *
+C* Permission to use, copy and distribute this software and its       *
+C* documentation strictly for non-commercial purposes is hereby       *
+C* granted without fee, provided that the above copyright notice      *
+C* appears in all copies, that both the copyright notice and this     *
+C* permission notice appear in the supporting documentation and that  *
+C* the use of this program is acknowledged in scientific              *
+C* publications (see reference above). The authors make no claims     *
+C* about the suitability of this software for any purpose. It is      *
+C* provided "as is" without express or implied warranty. Any change   *
+C* of the code should be submitted to the authors                     *
+C*                                                                    *
+C* To use this program at LHC energies, please make sure that         *
+C* "double precision" variables should better be real*16              *
+C**********************************************************************
+
+C======================================================================
+C
+C call this routine first to initialize some parameter needed in the
+C function. 
+C
+C gm = Gamma_cm, that is, gm of each ion (~100 for RHIC ~3000 for LHC)
+C mass = mass of the produced particle in MeV (~0.511 for e,~100 for mu)
+C======================================================================
+
+      SUBROUTINE Initdiffcross (gm,mass)
+      IMPLICIT NONE
+      DOUBLE PRECISION gm,mass
+
+      DOUBLE PRECISION gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
+      COMMON/PHYSPARAM/gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
+
+      DOUBLE PRECISION ARCOSH,x
+      ARCOSH(x)=LOG(x+SQRT(x**2-1D0))
+
+      gamma=gm
+      beta=sqrt(1D0-1D0/gamma**2)
+      m=mass
+      
+      w1xw1 = 1D0/gamma**2
+      w2xw2 = 1D0/gamma**2
+      w1xw2 = 2D0-1D0/gamma**2
+      wl=SQRT (w1xw1)
+      wy=ARCOSH (gamma)
+      
+      RETURN 
+      END
+
+C ============================================================================
+C
+C Diffcross calculates the fivefold differential cross section
+C
+C dsigma/ dp+t dp-t dy+ dy- ddeltaphi
+C
+C the trivial integration over the total phi-dependence is not included
+C
+C ppvt= absolute value of the positron transverse momentum (MeV)
+C pmvt=  - " -                electron - " -
+C dphi= phi-angle between the electron and positron transverse momentum
+C yp = rapidity of the positron
+C ym = rapidity of the electron
+C      defined to make E = sqrt(pt^2 + m^2) cosh (y)
+C      and             Pz= sqrt(pt^2 + m^2) sinh (y)
+C
+C dsigma = differential cross section in kbarn/MeV^4/ (Zalpha)**4
+C
+C to get the total cross section, you have to integrate over
+C
+C Integral dsigma dyp dym ddphi 2 pi ppt dppt pmt dpmt
+C
+C (see also source for sigma.f)
+C
+C======================================================================
+
+      SUBROUTINE Diffcross(ppvt,yp,pmvt,ym,dphi,dsigma)
+      
+      IMPLICIT NONE
+      
+      DOUBLE PRECISION ppvt,pmvt,dphi,dsigma,yp,ym
+      DOUBLE PRECISION ppvl,pmvl
+      
+      DOUBLE PRECISION k1(2)
+      DOUBLE PRECISION kd(2)
+      DOUBLE PRECISION kx(2)
+      DOUBLE PRECISION pmt(2)
+      DOUBLE PRECISION ppt(2)
+      DOUBLE PRECISION qb
+      DOUBLE PRECISION mk1(2),mk1d(2),mk1x(2)
+      DOUBLE PRECISION mkd(2),mkd1(2)
+      DOUBLE PRECISION mkx(2),mkx1(2)
+      
+      DOUBLE PRECISION Iz2,Id1,Id2,Id3,Iv1,Iv2
+      
+      DOUBLE PRECISION PI
+      PARAMETER (PI=3.141592653589793238462643D0)
+      
+      DOUBLE PRECISION gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
+      COMMON/PHYSPARAM/gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
+      
+      DOUBLE PRECISION m0,m1,md,mx
+      DOUBLE PRECISION N1,N2,N3,N4,N5
+      DOUBLE PRECISION N6,N7,N8,N9,N10
+      DOUBLE PRECISION N11,N12,N13,N14,N15
+      DOUBLE PRECISION N16,N17,N18,NT
+      
+      DOUBLE PRECISION pmlxpml,pmlxppl,pmlxql
+      DOUBLE PRECISION pmtxpmt,pmtxppt,pptxppt
+      DOUBLE PRECISION pplxppl,pplxql,qlxql
+      DOUBLE PRECISION w1xpml,w1xppl,w1xql
+      DOUBLE PRECISION w2xpml,w2xppl,w2xql
+      
+      DOUBLE PRECISION r1,rd,rx
+
+      INTEGER setzero
+      COMMON /SETZPARAM/ setzero
+      INTEGER badcount
+      COMMON/badpar/badcount
+
+      setzero=0
+      
+      ppt(1)=ppvt
+      ppt(2)=0D0
+      ppvl=SQRT(ppvt**2+m**2)
+      
+      pmt(1)=pmvt*cos(dphi)
+      pmt(2)=pmvt*sin(dphi)
+      pmvl=SQRT(pmvt**2+m**2)
+      
+      pmlxpml = pmvl**2
+      pplxppl = ppvl**2
+      pmlxppl = pmvl*ppvl*COSH(yp-ym)
+      w2xpml = wl*pmvl*COSH(ym+wy)
+      w2xppl = wl*ppvl*COSH(yp+wy)
+      w1xpml = wl*pmvl*COSH(ym-wy)
+      w1xppl = wl*ppvl*COSH(yp-wy)
+      w1xql=0D0
+      w2xql=w2xpml+w2xppl
+      qb=gamma*(w2xppl+w2xpml)/SINH(2D0*wy)
+      qlxql=-qb**2
+      pmlxql = qb*pmvl*SINH(wy-ym)
+      pplxql = qb*ppvl*SINH(wy-yp)
+      
+      pmtxpmt =-pmvt**2
+      pmtxppt =-pmvt*ppvt*cos(dphi)
+      pptxppt =-ppvt**2
+      
+C======================================================================
+C
+C Definition of propagatorterms
+C
+      m0 = -qlxql
+      
+      k1(1) = -ppt(1) -pmt(1)
+      k1(2) = -pmt(2)
+      m1=-qlxql-pplxppl-pmlxpml+2D0*(pplxql+pmlxql-pmlxppl)
+      
+      kd(1)= -pmt(1)
+      kd(2)= -pmt(2)
+      md= m**2 - (qlxql+pmlxpml-2D0*pmlxql)
+      
+      kx(1)= -ppt(1)   
+      kx(2)=0D0
+      mx= m**2 - (qlxql+pplxppl-2D0*pplxql)
+      
+      r1 = m1 - m0 + k1(1)*k1(1)+k1(2)*k1(2)
+      rd = md - m0 + kd(1)*kd(1)+kd(2)*kd(2)
+      rx = mx - m0 + kx(1)*kx(1)+kx(2)*kx(2)
+      
+      mk1(1)= -k1(1)
+      mk1(2)= -k1(2)
+      mkd(1)= -kd(1)
+      mkd(2)= -kd(2)
+      mkx(1)= -kx(1)
+      mkx(2)= -kx(2)
+      mk1d(1)= k1(1)-kd(1)
+      mk1d(2)= k1(2)-kd(2)
+      mk1x(1)= k1(1)-kx(1)
+      mk1x(2)= k1(2)-kx(2)
+      mkd1(1)= kd(1)-k1(1)
+      mkd1(2)= kd(2)-k1(2)
+      mkx1(1)= kx(1)-k1(1)
+      mkx1(2)= kx(2)-k1(2)
+
+C
+C calculate all different integrals
+C
+
+      N1 = Id1(mkd,mk1d,md,m0,m1) * (  - 2*w1xw1*w2xw2 )
+      
+      N2 = Id1(mkx,mk1x,mx,m0,m1) * (  - 2*w1xw1*w2xw2 )
+      
+      N3 = Iv1(mk1,mkd1,mkx1,m1,m0,md,mx) * ( 4*w1xw1*w2xw2*
+     &  pmlxpml + 8*w1xw1*w2xw2*pmlxppl - 8*w1xw1*w2xw2*pmlxql + 4*
+     &  w1xw1*w2xw2*pmtxpmt + 8*w1xw1*w2xw2*pmtxppt + 4*w1xw1*w2xw2*
+     &  pplxppl - 8*w1xw1*w2xw2*pplxql + 4*w1xw1*w2xw2*pptxppt + 4*
+     &  w1xw1*w2xw2*rd + 4*w1xw1*w2xw2*rx - 16*w1xw1*w2xpml*w2xppl + 
+     &  8*w1xw1*w2xpml*w2xql + 8*w1xw1*w2xppl*w2xql - 8*w1xw2**2*
+     &  pmlxpml - 16*w1xw2**2*pmlxppl + 16*w1xw2**2*pmlxql - 8*
+     &  w1xw2**2*pmtxpmt - 16*w1xw2**2*pmtxppt - 8*w1xw2**2*pplxppl
+     &  + 16*w1xw2**2*pplxql - 8*w1xw2**2*pptxppt - 8*w1xw2**2*rd - 
+     &  8*w1xw2**2*rx + 16*w1xw2*w1xpml*w2xpml + 16*w1xw2*w1xpml*
+     &  w2xppl - 16*w1xw2*w1xpml*w2xql + 16*w1xw2*w1xppl*w2xpml + 16*
+     &  w1xw2*w1xppl*w2xppl - 16*w1xw2*w1xppl*w2xql - 16*w1xw2*w1xql*
+     &  w2xpml - 16*w1xw2*w1xql*w2xppl - 8*w1xpml**2*w2xw2 + 8*w1xpml
+     &  *w1xql*w2xw2 - 8*w1xppl**2*w2xw2 + 8*w1xppl*w1xql*w2xw2 )
+      
+      N4 = Id1(mk1,mkd1,m1,m0,md) * (  - 2*w1xw1*w2xw2 + 8*
+     &  w1xw2**2 )
+      
+      N5 = Id2(mk1d,mkd,md,m1,m0) * ( 4*w1xw1*w2xw2*pmlxppl + 4*
+     &  w1xw1*w2xw2*pmtxppt - 4*w1xw1*w2xw2*pplxql + 4*w1xw1*w2xw2*
+     &  m**2 + 2*w1xw1*w2xw2*r1 - 2*w1xw1*w2xw2*rd - 8*w1xw1*w2xpml*
+     &  w2xppl + 8*w1xw1*w2xppl*w2xql )
+      
+      N6 = Id1(mk1,mkx1,m1,m0,mx) * (  - 2*w1xw1*w2xw2 + 8*
+     &  w1xw2**2 )
+      
+      N7 = Id2(mk1x,mkx,mx,m1,m0) * ( 4*w1xw1*w2xw2*pmlxppl - 4*
+     &  w1xw1*w2xw2*pmlxql + 4*w1xw1*w2xw2*pmtxppt + 4*w1xw1*w2xw2*
+     &  m**2 + 2*w1xw1*w2xw2*r1 - 2*w1xw1*w2xw2*rx - 8*w1xw1*w2xpml*
+     &  w2xppl + 8*w1xw1*w2xpml*w2xql )
+      
+      N8 = Id1(k1,kd,m0,m1,md) * ( 2*w1xw1*w2xw2 )
+      
+      N9 = Id2(kd,k1,m0,md,m1) * (  - 4*w1xw1*w2xw2*pmlxpml
+     &  + 4*w1xw1*w2xw2*pmlxql - 4*w1xw1*w2xw2*pmtxpmt + 4*w1xw1*
+     &  w2xw2*m**2 - 2*w1xw1*w2xw2*rd + 8*w1xpml**2*w2xw2 - 8*w1xpml*
+     &  w1xql*w2xw2 )
+      
+      N10 = Id1(k1,kx,m0,m1,mx) * ( 2*w1xw1*w2xw2 )
+      
+      N11 = Id2(kx,k1,m0,mx,m1) * (  - 4*w1xw1*w2xw2*pplxppl
+     &  + 4*w1xw1*w2xw2*pplxql - 4*w1xw1*w2xw2*pptxppt + 4*w1xw1*
+     &  w2xw2*m**2 - 2*w1xw1*w2xw2*rx + 8*w1xppl**2*w2xw2 - 8*w1xppl*
+     &  w1xql*w2xw2 )
+      
+      N12 = Iv2(k1,kd,kx,m0,m1,md,mx) * ( 8*w1xw1*w2xw2*
+     &  pmlxpml*pplxppl - 8*w1xw1*w2xw2*pmlxpml*pplxql + 8*w1xw1*
+     &  w2xw2*pmlxpml*pptxppt - 8*w1xw1*w2xw2*pmlxpml*m**2 + 4*w1xw1*
+     &  w2xw2*pmlxpml*rx - 8*w1xw1*w2xw2*pmlxppl*qlxql - 8*w1xw1*
+     &  w2xw2*pmlxppl*m0 - 8*w1xw1*w2xw2*pmlxql*pplxppl + 16*w1xw1
+     &  *w2xw2*pmlxql*pplxql - 8*w1xw1*w2xw2*pmlxql*pptxppt + 8*w1xw1
+     &  *w2xw2*pmlxql*m**2 - 8*w1xw1*w2xw2*pmlxql*rx + 8*w1xw1*w2xw2*
+     &  pmtxpmt*pplxppl - 8*w1xw1*w2xw2*pmtxpmt*pplxql + 8*w1xw1*
+     &  w2xw2*pmtxpmt*pptxppt - 8*w1xw1*w2xw2*pmtxpmt*m**2 + 4*w1xw1*
+     &  w2xw2*pmtxpmt*rx - 8*w1xw1*w2xw2*pmtxppt*qlxql - 8*w1xw1*
+     &  w2xw2*pmtxppt*m0 - 8*w1xw1*w2xw2*pplxppl*m**2 + 4*w1xw1*
+     &  w2xw2*pplxppl*rd + 8*w1xw1*w2xw2*pplxql*m**2 - 8*w1xw1*w2xw2*
+     &  pplxql*rd - 8*w1xw1*w2xw2*pptxppt*m**2 + 4*w1xw1*w2xw2*
+     &  pptxppt*rd - 8*w1xw1*w2xw2*qlxql*m**2 + 8*w1xw1*w2xw2*m**4 - 
+     &  8*w1xw1*w2xw2*m**2*m0 - 4*w1xw1*w2xw2*m**2*rd - 4*w1xw1*
+     &  w2xw2*m**2*rx + 4*w1xw1*w2xw2*rd*rx + 16*w1xw1*w2xpml*w2xppl*
+     &  qlxql +  16*w1xw1*w2xpml*
+     &  w2xppl*m0 - 16*w1xw1*w2xpml*w2xql*pplxql + 8*w1xw1*w2xpml*
+     &  w2xql*rx - 16*w1xw1*w2xppl*w2xql*pmlxql + 8*w1xw1*w2xppl*
+     &  w2xql*rd + 16*w1xw1*w2xql**2*pmlxppl + 16*w1xw1*w2xql**2*
+     &  pmtxppt + 16*w1xw1*w2xql**2*m**2 - 16*w1xw2**2*pmlxpml*
+     &  pplxppl + 16*w1xw2**2*pmlxpml*pplxql - 16*w1xw2**2*pmlxpml*
+     &  pptxppt + 16*w1xw2**2*pmlxpml*m**2 - 8*w1xw2**2*pmlxpml*rx + 
+     &  16*w1xw2**2*pmlxppl*qlxql + 16*w1xw2**2*pmlxppl*m0 + 16*
+     &  w1xw2**2*pmlxql*pplxppl - 32*w1xw2**2*pmlxql*pplxql + 16*
+     &  w1xw2**2*pmlxql*pptxppt - 16*w1xw2**2*pmlxql*m**2 + 16*
+     &  w1xw2**2*pmlxql*rx - 16*w1xw2**2*pmtxpmt*pplxppl + 16*
+     &  w1xw2**2*pmtxpmt*pplxql - 16*w1xw2**2*pmtxpmt*pptxppt + 16*
+     &  w1xw2**2*pmtxpmt*m**2 - 8*w1xw2**2*pmtxpmt*rx + 16*w1xw2**2*
+     &  pmtxppt*qlxql + 16*w1xw2**2*pmtxppt*m0 + 16*w1xw2**2*
+     &  pplxppl*m**2 - 8*w1xw2**2*pplxppl*rd - 16*w1xw2**2*pplxql*
+     &  m**2 + 16*w1xw2**2*pplxql*rd + 16*w1xw2**2*pptxppt*m**2 - 8*
+     &  w1xw2**2*pptxppt*rd + 16*w1xw2**2*qlxql
+     &  *m**2 - 16*w1xw2**2*m**4 + 16*w1xw2**2*m**2*m0 + 8*
+     &  w1xw2**2*m**2*rd + 8*w1xw2**2*m**2*rx - 8*w1xw2**2*rd*rx + 32
+     &  *w1xw2*w1xpml*w2xpml*pplxppl - 32*w1xw2*w1xpml*w2xpml*pplxql
+     &  + 32*w1xw2*w1xpml*w2xpml*pptxppt - 32*w1xw2*w1xpml*w2xpml*
+     &  m**2 + 16*w1xw2*w1xpml*w2xpml*rx - 16*w1xw2*w1xpml*w2xppl*
+     &  qlxql - 16*w1xw2*w1xpml*w2xppl*m0 - 16*w1xw2*w1xpml*w2xql*
+     &  pplxppl + 32*w1xw2*w1xpml*w2xql*pplxql - 16*w1xw2*w1xpml*
+     &  w2xql*pptxppt + 16*w1xw2*w1xpml*w2xql*m**2 - 16*w1xw2*w1xpml*
+     &  w2xql*rx - 16*w1xw2*w1xppl*w2xpml*qlxql - 16*w1xw2*w1xppl*
+     &  w2xpml*m0 + 32*w1xw2*w1xppl*w2xppl*pmlxpml - 32*w1xw2*
+     &  w1xppl*w2xppl*pmlxql + 32*w1xw2*w1xppl*w2xppl*pmtxpmt - 32*
+     &  w1xw2*w1xppl*w2xppl*m**2 + 16*w1xw2*w1xppl*w2xppl*rd - 16*
+     &  w1xw2*w1xppl*w2xql*pmlxpml + 32*w1xw2*w1xppl*w2xql*pmlxql - 
+     &  16*w1xw2*w1xppl*w2xql*pmtxpmt + 16*w1xw2*w1xppl*w2xql*m**2 - 
+     &  16*w1xw2*w1xppl*w2xql*rd - 16*w1xw2*w1xql*w2xpml*pplxppl + 32
+     &  *w1xw2*w1xql*w2xpml*pplxql - 16*w1xw2*w1xql
+     &  *w2xpml*pptxppt + 16*w1xw2*w1xql*w2xpml*m**2 - 16*w1xw2*w1xql
+     &  *w2xpml*rx - 16*w1xw2*w1xql*w2xppl*pmlxpml + 32*w1xw2*w1xql*
+     &  w2xppl*pmlxql - 16*w1xw2*w1xql*w2xppl*pmtxpmt + 16*w1xw2*
+     &  w1xql*w2xppl*m**2 - 16*w1xw2*w1xql*w2xppl*rd - 32*w1xw2*w1xql
+     &  *w2xql*pmlxppl - 32*w1xw2*w1xql*w2xql*pmtxppt - 32*w1xw2*
+     &  w1xql*w2xql*m**2 - 16*w1xpml**2*w2xw2*pplxppl + 16*w1xpml**2*
+     &  w2xw2*pplxql - 16*w1xpml**2*w2xw2*pptxppt + 16*w1xpml**2*
+     &  w2xw2*m**2 - 8*w1xpml**2*w2xw2*rx + 32*w1xpml*w1xppl*w2xw2*
+     &  pmlxppl - 16*w1xpml*w1xppl*w2xw2*pmlxql + 32*w1xpml*w1xppl*
+     &  w2xw2*pmtxppt - 16*w1xpml*w1xppl*w2xw2*pplxql + 16*w1xpml*
+     &  w1xppl*w2xw2*qlxql + 32*w1xpml*w1xppl*w2xw2*m**2 + 16*w1xpml*
+     &  w1xppl*w2xw2*m0 + 8*w1xpml*w1xppl*w2xw2*rd + 8*w1xpml*
+     &  w1xppl*w2xw2*rx - 64*w1xpml*w1xppl*w2xpml*w2xppl + 32*w1xpml*
+     &  w1xppl*w2xpml*w2xql + 32*w1xpml*w1xppl*w2xppl*w2xql - 32*
+     &  w1xpml*w1xppl*w2xql**2 - 16*w1xpml*w1xql*w2xw2*pmlxppl - 16*
+     &  w1xpml*w1xql*w2xw2*pmtxppt + 16*w1xpml*w1xql*
+     &  w2xw2*pplxppl - 16*w1xpml*w1xql*w2xw2*pplxql + 16*w1xpml*
+     &  w1xql*w2xw2*pptxppt - 32*w1xpml*w1xql*w2xw2*m**2 + 8*w1xpml*
+     &  w1xql*w2xw2*rx + 32*w1xpml*w1xql*w2xpml*w2xppl - 16*w1xppl**2
+     &  *w2xw2*pmlxpml + 16*w1xppl**2*w2xw2*pmlxql - 16*w1xppl**2*
+     &  w2xw2*pmtxpmt + 16*w1xppl**2*w2xw2*m**2 - 8*w1xppl**2*w2xw2*
+     &  rd + 16*w1xppl*w1xql*w2xw2*pmlxpml - 16*w1xppl*w1xql*w2xw2*
+     &  pmlxppl - 16*w1xppl*w1xql*w2xw2*pmlxql + 16*w1xppl*w1xql*
+     &  w2xw2*pmtxpmt - 16*w1xppl*w1xql*w2xw2*pmtxppt - 32*w1xppl*
+     &  w1xql*w2xw2*m**2 + 8*w1xppl*w1xql*w2xw2*rd + 32*w1xppl*w1xql*
+     &  w2xpml*w2xppl + 16*w1xql**2*w2xw2*pmlxppl + 16*w1xql**2*w2xw2
+     &  *pmtxppt + 16*w1xql**2*w2xw2*m**2 - 32*w1xql**2*w2xpml*w2xppl
+     &  )
+      
+      N13 = Id2(k1,kd,m0,m1,md) * ( 4*w1xw1*w2xw2*pmlxql + 4*
+     &  w1xw1*w2xw2*pplxql - 2*w1xw1*w2xw2*r1 - 8*w1xw1*w2xpml*w2xql
+     &  - 8*w1xw1*w2xppl*w2xql + 8*w1xw2**2*pmlxpml - 16*w1xw2**2*
+     &  pmlxql + 8*w1xw2**2*pmtxpmt - 8*w1xw2**2*m**2 + 8*w1xw2**2*rd
+     &  - 16*w1xw2*w1xpml*w2xpml + 16*w1xw2*w1xpml*w2xppl + 16*w1xw2
+     &  *w1xpml*w2xql + 16*w1xw2*w1xql*w2xpml - 16*w1xpml*w1xppl*
+     &  w2xw2 )
+      
+      N14 = Id3(k1,kd,m0,m1,md) * ( 4*w1xw1*w2xw2*pmlxpml*
+     &  pmlxppl + 4*w1xw1*w2xw2*pmlxpml*pmtxppt - 8*w1xw1*w2xw2*
+     &  pmlxpml*pplxql + 4*w1xw1*w2xw2*pmlxpml*m**2 + 4*w1xw1*w2xw2*
+     &  pmlxpml*r1 - 4*w1xw1*w2xw2*pmlxpml*rd + 4*w1xw1*w2xw2*pmlxppl
+     &  *pmtxpmt - 4*w1xw1*w2xw2*pmlxppl*qlxql - 4*w1xw1*w2xw2*
+     &  pmlxppl*m**2 - 4*w1xw1*w2xw2*pmlxppl*m0 + 8*w1xw1*w2xw2*
+     &  pmlxql*pplxql - 4*w1xw1*w2xw2*pmlxql*r1 + 4*w1xw1*w2xw2*
+     &  pmlxql*rd + 4*w1xw1*w2xw2*pmtxpmt*pmtxppt - 8*w1xw1*w2xw2*
+     &  pmtxpmt*pplxql + 4*w1xw1*w2xw2*pmtxpmt*m**2 + 4*w1xw1*w2xw2*
+     &  pmtxpmt*r1 - 4*w1xw1*w2xw2*pmtxpmt*rd - 4*w1xw1*w2xw2*pmtxppt
+     &  *qlxql - 4*w1xw1*w2xw2*pmtxppt*m**2 - 4*w1xw1*w2xw2*pmtxppt*
+     &  m0 + 8*w1xw1*w2xw2*pplxql*m**2 - 4*w1xw1*w2xw2*pplxql*rd
+     &  - 4*w1xw1*w2xw2*qlxql*m**2 - 4*w1xw1*w2xw2*m**4 - 4*w1xw1*
+     &  w2xw2*m**2*m0 - 4*w1xw1*w2xw2*m**2*r1 + 4*w1xw1*w2xw2*m**2
+     &  *rd + 2*w1xw1*w2xw2*r1*rd - 2*w1xw1*w2xw2*rd**2 - 8*w1xw1*
+     &  w2xpml*w2xppl*pmlxpml - 8*w1xw1*w2xpml*w2xppl*pmtxpmt + 8*
+     &  w1xw1*w2xpml*w2xppl*qlxql + 8*w1xw1*w2xpml*w2xppl*m**2
+     &  + 8*w1xw1*w2xpml*w2xppl*m0 + 16*w1xw1*w2xppl*w2xql*
+     &  pmlxpml - 16*w1xw1*w2xppl*w2xql*pmlxql + 16*w1xw1*w2xppl*
+     &  w2xql*pmtxpmt - 16*w1xw1*w2xppl*w2xql*m**2 + 8*w1xw1*w2xppl*
+     &  w2xql*rd - 16*w1xw2*w1xpml*w2xppl*pmlxpml + 32*w1xw2*w1xpml*
+     &  w2xppl*pmlxql - 16*w1xw2*w1xpml*w2xppl*pmtxpmt - 16*w1xw2*
+     &  w1xpml*w2xppl*qlxql + 16*w1xw2*w1xpml*w2xppl*m**2 - 16*w1xw2*
+     &  w1xpml*w2xppl*m0 - 16*w1xw2*w1xpml*w2xppl*rd - 16*
+     &  w1xpml**2*w2xw2*pmlxppl - 16*w1xpml**2*w2xw2*pmtxppt + 16*
+     &  w1xpml**2*w2xw2*pplxql - 16*w1xpml**2*w2xw2*m**2 - 8*
+     &  w1xpml**2*w2xw2*r1 + 8*w1xpml**2*w2xw2*rd + 32*w1xpml**2*
+     &  w2xpml*w2xppl - 32*w1xpml**2*w2xppl*w2xql + 8*w1xpml*w1xppl*
+     &  w2xw2*pmlxpml - 16*w1xpml*w1xppl*w2xw2*pmlxql + 8*w1xpml*
+     &  w1xppl*w2xw2*pmtxpmt + 8*w1xpml*w1xppl*w2xw2*qlxql - 8*w1xpml
+     &  *w1xppl*w2xw2*m**2 + 8*w1xpml*w1xppl*w2xw2*m0 + 8*w1xpml*
+     &  w1xppl*w2xw2*rd + 16*w1xpml*w1xql*w2xw2*pmlxppl + 16*w1xpml*
+     &  w1xql*w2xw2*pmtxppt - 16*w1xpml*w1xql*w2xw2*
+     &  pplxql + 16*w1xpml*w1xql*w2xw2*m**2 + 8*w1xpml*w1xql*w2xw2*r1
+     &  - 8*w1xpml*w1xql*w2xw2*rd - 32*w1xpml*w1xql*w2xpml*w2xppl + 
+     &  32*w1xpml*w1xql*w2xppl*w2xql )
+      
+      N15 = Id2(k1,kx,m0,m1,mx) * ( 4*w1xw1*w2xw2*pmlxql + 4*
+     &  w1xw1*w2xw2*pplxql - 2*w1xw1*w2xw2*r1 - 8*w1xw1*w2xpml*w2xql
+     &  - 8*w1xw1*w2xppl*w2xql + 8*w1xw2**2*pplxppl - 16*w1xw2**2*
+     &  pplxql + 8*w1xw2**2*pptxppt - 8*w1xw2**2*m**2 + 8*w1xw2**2*rx
+     &  + 16*w1xw2*w1xppl*w2xpml - 16*w1xw2*w1xppl*w2xppl + 16*w1xw2
+     &  *w1xppl*w2xql + 16*w1xw2*w1xql*w2xppl - 16*w1xpml*w1xppl*
+     &  w2xw2 )
+      
+      N16 = Id3(k1,kx,m0,m1,mx) * ( 4*w1xw1*w2xw2*pmlxppl*
+     &  pplxppl + 4*w1xw1*w2xw2*pmlxppl*pptxppt - 4*w1xw1*w2xw2*
+     &  pmlxppl*qlxql - 4*w1xw1*w2xw2*pmlxppl*m**2 - 4*w1xw1*w2xw2*
+     &  pmlxppl*m0 - 8*w1xw1*w2xw2*pmlxql*pplxppl + 8*w1xw1*w2xw2*
+     &  pmlxql*pplxql - 8*w1xw1*w2xw2*pmlxql*pptxppt + 8*w1xw1*w2xw2*
+     &  pmlxql*m**2 - 4*w1xw1*w2xw2*pmlxql*rx + 4*w1xw1*w2xw2*pmtxppt
+     &  *pplxppl + 4*w1xw1*w2xw2*pmtxppt*pptxppt - 4*w1xw1*w2xw2*
+     &  pmtxppt*qlxql - 4*w1xw1*w2xw2*pmtxppt*m**2 - 4*w1xw1*w2xw2*
+     &  pmtxppt*m0 + 4*w1xw1*w2xw2*pplxppl*m**2 + 4*w1xw1*w2xw2*
+     &  pplxppl*r1 - 4*w1xw1*w2xw2*pplxppl*rx - 4*w1xw1*w2xw2*pplxql*
+     &  r1 + 4*w1xw1*w2xw2*pplxql*rx + 4*w1xw1*w2xw2*pptxppt*m**2 + 4
+     &  *w1xw1*w2xw2*pptxppt*r1 - 4*w1xw1*w2xw2*pptxppt*rx - 4*w1xw1*
+     &  w2xw2*qlxql*m**2 - 4*w1xw1*w2xw2*m**4 - 4*w1xw1*w2xw2*m**2*
+     &  m0 - 4*w1xw1*w2xw2*m**2*r1 + 4*w1xw1*w2xw2*m**2*rx + 2*
+     &  w1xw1*w2xw2*r1*rx - 2*w1xw1*w2xw2*rx**2 - 8*w1xw1*w2xpml*
+     &  w2xppl*pplxppl - 8*w1xw1*w2xpml*w2xppl*pptxppt + 8*w1xw1*
+     &  w2xpml*w2xppl*qlxql + 8*w1xw1*w2xpml*w2xppl*m**2
+     &  + 8*w1xw1*w2xpml*w2xppl*m0 + 16*w1xw1*w2xpml*w2xql*
+     &  pplxppl - 16*w1xw1*w2xpml*w2xql*pplxql + 16*w1xw1*w2xpml*
+     &  w2xql*pptxppt - 16*w1xw1*w2xpml*w2xql*m**2 + 8*w1xw1*w2xpml*
+     &  w2xql*rx - 16*w1xw2*w1xppl*w2xpml*pplxppl + 32*w1xw2*w1xppl*
+     &  w2xpml*pplxql - 16*w1xw2*w1xppl*w2xpml*pptxppt - 16*w1xw2*
+     &  w1xppl*w2xpml*qlxql + 16*w1xw2*w1xppl*w2xpml*m**2 - 16*w1xw2*
+     &  w1xppl*w2xpml*m0 - 16*w1xw2*w1xppl*w2xpml*rx + 8*w1xpml*
+     &  w1xppl*w2xw2*pplxppl - 16*w1xpml*w1xppl*w2xw2*pplxql + 8*
+     &  w1xpml*w1xppl*w2xw2*pptxppt + 8*w1xpml*w1xppl*w2xw2*qlxql - 8
+     &  *w1xpml*w1xppl*w2xw2*m**2 + 8*w1xpml*w1xppl*w2xw2*m0 + 8*
+     &  w1xpml*w1xppl*w2xw2*rx - 16*w1xppl**2*w2xw2*pmlxppl + 16*
+     &  w1xppl**2*w2xw2*pmlxql - 16*w1xppl**2*w2xw2*pmtxppt - 16*
+     &  w1xppl**2*w2xw2*m**2 - 8*w1xppl**2*w2xw2*r1 + 8*w1xppl**2*
+     &  w2xw2*rx + 32*w1xppl**2*w2xpml*w2xppl - 32*w1xppl**2*w2xpml*
+     &  w2xql + 16*w1xppl*w1xql*w2xw2*pmlxppl - 16*w1xppl*w1xql*w2xw2
+     &  *pmlxql + 16*w1xppl*w1xql*w2xw2*
+     &  pmtxppt + 16*w1xppl*w1xql*w2xw2*m**2 + 8*w1xppl*w1xql*w2xw2*
+     &  r1 - 8*w1xppl*w1xql*w2xw2*rx - 32*w1xppl*w1xql*w2xpml*w2xppl
+     &  + 32*w1xppl*w1xql*w2xpml*w2xql )
+      
+      N17 = Iz2(k1,m0,m1) * (  - 8*w1xw2**2 )
+      
+      N18 = Id1(mkd1,mkx1,m1,md,mx) * ( 4*w1xw1*w2xw2 - 8*w1xw2**2 )
+      
+C
+C dsigma is summ of all terms
+C
+      NT=N1+N2+N3+N4+N5+N6+N7+N8+N9+N10+N11+N12+N13+N14+
+     &  N15+N16+N17+N18
+
+C
+C correction from w/u
+      NT=NT*4D0/beta**2
+C 1/(2pi)**6 from d3p, (2*pi)**2 from F.T.
+      NT=NT/(2*pi)**6*(2*pi)**2
+C from 1/2E+ 1/2E-
+      NT=NT/4D0
+C transform from MeV^-2 to kbarn
+      NT=NT*(1.9733D0)**2/10D0
+      
+      dsigma=NT
+      IF((setzero.EQ.1).OR.(dsigma.LT.0)) THEN
+      dsigma=0D0
+      badcount=badcount+1
+      ENDIF 
+      END
+      
+C========================================================================
+C All the differential integral forms are calculated here
+C
+      DOUBLE PRECISION FUNCTION Iz0(x,u,v)
+      IMPLICIT NONE
+      
+      DOUBLE PRECISION x(2)
+      DOUBLE PRECISION u,v
+      DOUBLE PRECISION s,arg
+      
+      INTEGER setzero
+      COMMON /SETZPARAM/ setzero
+      
+      DOUBLE PRECISION tepxx
+      
+      DOUBLE PRECISION PI
+      PARAMETER (PI=3.141592653589793238462643D0)
+      
+      tepxx=x(1)*x(1)+x(2)*x(2)
+      s=SQRT((tepxx+u+v)**2 - 4*u*v)
+      arg=(tepxx+u+v+s)**2/(4*u*v)
+      IF(arg.lt.0D0) THEN
+      setzero=1
+      Iz0=0D0
+      ELSE 
+      Iz0=pi * LOG((tepxx+u+v+s)**2/(4*u*v)) / s        
+      ENDIF 
+      END
+      
+C     ---------------------------------------------------------------     
+
+      DOUBLE PRECISION FUNCTION Iz1(x,u,v)
+      IMPLICIT NONE
+      DOUBLE PRECISION Iz0
+      EXTERNAL         Iz0
+
+      DOUBLE PRECISION x(2)
+      DOUBLE PRECISION u,v
+      DOUBLE PRECISION a,b,s
+      
+      DOUBLE PRECISION tepxx
+      DOUBLE PRECISION PI
+      PARAMETER (PI=3.141592653589793238462643D0)
+      
+      tepxx=x(1)*x(1)+x(2)*x(2)
+      
+      s=SQRT((tepxx+u+v)**2 - 4*u*v)
+      A= (2*(tepxx+u-v+s)) / (s**2*(tepxx+u+v+s)) 
+     &  - 1/(s*u)
+      B= -Iz0(x,u,v)/pi * (tepxx+u-v)/s**2     
+      Iz1=(-1)*pi*(A + B)
+      
+      END
+      
+C     -----------------------------------------------------------
+      
+      DOUBLE PRECISION FUNCTION Iz2(x,u,v)
+      IMPLICIT NONE
+      
+      DOUBLE PRECISION Iz0,Iz1
+      EXTERNAL         Iz0,Iz1
+
+      DOUBLE PRECISION x(2)
+      DOUBLE PRECISION u,v
+      DOUBLE PRECISION s,A,B,C,D,E
+      
+      DOUBLE PRECISION tepxx
+      
+      DOUBLE PRECISION PI
+      PARAMETER (PI=3.141592653589793238462643D0)
+
+      tepxx=x(1)*x(1)+x(2)*x(2)
+      
+      s=SQRT((tepxx+u+v)**2 - 4*u*v)
+      A= (2*(tepxx-u+v-s))/(s**3*(tepxx+u+v+s))
+      B= -2*(tepxx+u-v+s)*( 
+     &  2*(tepxx-u+v)*(tepxx+u+v+s) +
+     &  s*(tepxx-u+v+s)  ) / 
+     &  (s**4*(tepxx+u+v+s)**2)
+      C= (tepxx-u+v)/(u*s**3)
+      D= (Iz1(x,v,u)/pi)*(tepxx+u-v)/(s**2)
+      E= (Iz0(x,u,v)/pi)*( 
+     &  1/s**2 + (tepxx+u-v)*2*(tepxx-u+v) /
+     &  (s**4) )
+      Iz2= pi*(A + B + C +D + E)
+      
+      END
+
+C     -------------------------------------------------------------
+      
+      DOUBLE PRECISION FUNCTION Id0(x,y,u,v,w)
+      IMPLICIT NONE
+      DOUBLE PRECISION Iz0,Iz1
+      EXTERNAL         Iz0,Iz1
+      DOUBLE PRECISION x(2),y(2),dyx(2)
+      DOUBLE PRECISION u,v,w
+      DOUBLE PRECISION A,AXY,B,C,D,E,RX,RY
+      DOUBLE PRECISION tepxx,tepxy,tepyy
+
+      tepxx=x(1)*x(1)+x(2)*x(2)
+      tepyy=y(1)*y(1)+y(2)*y(2)
+      tepxy=x(1)*y(1)+x(2)*y(2)
+      
+      rx= v-u+tepxx
+      ry= w-u+tepyy
+      axy= x(1)*y(2)-x(2)*y(1)
+      
+      A= rx**2*tepyy - 2*rx*ry*tepxy + ry**2*tepxx
+      B= 4*u*axy**2 + A
+      C= 2*axy**2 + (rx+ry)*tepxy - rx*tepyy - ry*tepxx
+      D= rx*tepyy - ry*tepxy
+      E= ry*tepxx - rx*tepxy
+      dyx(1)= y(1)-x(1)
+      dyx(2)= y(2)-x(2)
+      Id0=1/B*(C*Iz0(dyx,v,w) + D*Iz0(y,u,w) + E*Iz0(x,u,v))
+      
+      END
+      
+C     -------------------------------------------------------------
+      
+      DOUBLE PRECISION FUNCTION Id1(x,y,u,v,w)
+      IMPLICIT NONE
+      DOUBLE PRECISION Id0,Iz0,Iz1
+      EXTERNAL         Id0,Iz0,Iz1
+      
+      DOUBLE PRECISION x(2),y(2),dyx(2)
+      DOUBLE PRECISION u,v,w
+      
+      DOUBLE PRECISION A,AXY,B,C,D,E,F,G,H,RX,RY
+      DOUBLE PRECISION tepxx,tepxy,tepyy
+      
+      tepxx=x(1)*x(1)+x(2)*x(2)
+      tepxy=x(1)*y(1)+x(2)*y(2)
+      tepyy=y(1)*y(1)+y(2)*y(2)
+      
+      rx= v-u+tepxx
+      ry= w-u+tepyy
+      axy= x(1)*y(2)-x(2)*y(1)
+      A= rx**2*tepyy - 2*rx*ry*tepxy + ry**2*tepxx
+      B= 4*u*axy**2 + A
+      C= -4*axy**2 -2*(-rx*tepyy+(rx+ry)*tepxy-ry*tepxx)
+      D= -2*tepxy + tepxx + tepyy
+      E= -tepyy + tepxy
+      F= rx*tepyy - ry*tepxy
+      G= -tepxx + tepxy
+      H= ry*tepxx - rx*tepxy
+      
+      dyx(1)= y(1)-x(1)
+      dyx(2)= y(2)-x(2)
+      
+      Id1= -( (C/B)*Id0(x,y,u,v,w) + (1/B)*(
+     &  D*Iz0(dyx,v,w) +
+     &  E*Iz0(y,u,w) - F*Iz1(y,u,w) +
+     &  G*Iz0(x,u,v) - H*Iz1(x,u,v)   ) )
+      
+      END
+      
+C     --------------------------------------------------------
+      
+      DOUBLE PRECISION FUNCTION Id2(x,y,u,v,w)
+      IMPLICIT NONE
+      DOUBLE PRECISION Id1,Id0,Iz0,Iz1,Iz2
+      EXTERNAL         Id1,Id0,Iz0,Iz1,Iz2
+      
+      DOUBLE PRECISION x(2),y(2),mx(2),dyx(2)
+      DOUBLE PRECISION u,v,w
+      
+      DOUBLE PRECISION A0,A1,A2,A3,A4,AXY
+      DOUBLE PRECISION B,C1,C2,C3,D1,D2,D3
+      DOUBLE PRECISION E1,E2,E3,F1,F2,F3
+      DOUBLE PRECISION G1,G2,G3,RX,RY
+      DOUBLE PRECISION tepxx,tepxy,tepyy
+      
+      tepxx=x(1)*x(1)+x(2)*x(2)
+      tepxy=x(1)*y(1)+x(2)*y(2)
+      tepyy=y(1)*y(1)+y(2)*y(2)
+      
+      mx(1)=-x(1)
+      mx(2)=-x(2)
+      dyx(1)= y(1)-x(1)
+      dyx(2)= y(2)-x(2)
+      rx= v-u+tepxx
+      ry= w-u+tepyy
+      axy= x(1)*y(2)-x(2)*y(1)
+      
+      A0= rx**2*tepyy - 2*rx*ry*tepxy + ry**2*tepxx
+      B= 4*u*axy**2 + A0
+      A1= 2*(tepyy-tepxy)*B
+      A2= (-4*axy**2 +2*(rx*tepyy-(rx+ry)*tepxy+ry*tepxx))*
+     &  2*2*(rx*tepyy-ry*tepxy)
+      A3= -4*axy**2 +2*(rx*tepyy-(rx+ry)*tepxy+ry*tepxx)
+      A4= -2*(rx*tepyy-ry*tepxy)
+      C1= tepxy-tepyy
+      D1= 2*axy**2+(rx+ry)*tepxy-rx*tepyy-ry*tepxx
+      E1= tepyy
+      F1= -tepxy
+      G1= ry*tepxx-rx*tepxy
+      C2= -2*tepxy+tepxx+tepyy
+      D2= -tepyy + tepxy
+      E2= rx*tepyy - ry*tepxy
+      F2= -tepxx + tepxy
+      G2= ry*tepxx - rx*tepxy
+      C3= -2*tepxy+tepxx+tepyy
+      D3= -tepyy
+      E3= -tepxx+tepxy
+      F3= tepxy
+      G3= -(ry*tepxx-rx*tepxy)
+      
+      Id2= ((A1-A2)/B**2)*Id0(x,y,u,v,w) + 
+     &  (A3/B**2)*( 
+     &  C1*Iz0(dyx,v,w) - D1*Iz1(dyx,v,w) + 
+     &  E1*Iz0(y,u,w) +
+     &  F1*Iz0(x,u,v) - G1*Iz1(mx,v,u) )  +      
+     &  (A4/B**2)*(                                        
+     &  C2*Iz0(dyx,v,w) +
+     &  D2*Iz0(y,u,w) - E2*Iz1(y,u,w) +
+     &  F2*Iz0(x,u,v) - G2*Iz1(x,u,v)  ) +
+     &  (1/B)*(
+     &  - C3*Iz1(dyx,v,w) 
+     &  + D3*Iz1(y,u,w)
+     &  - E3*Iz1(mx,v,u) + F3*Iz1(x,u,v) - G3*Iz2(mx,v,u) )
+      
+      END
+      
+C     --------------------------------------------------------
+      
+      DOUBLE PRECISION FUNCTION Id3(x,y,u,v,w)
+      IMPLICIT NONE
+      DOUBLE PRECISION Id2,Id1,Id0,Iz0,Iz1,Iz2
+      EXTERNAL         Id2,Id1,Id0,Iz0,Iz1,Iz2
+      
+      DOUBLE PRECISION x(2),y(2),dyx(2),mdyx(2),mx(2),my(2)
+      DOUBLE PRECISION u,v,w
+      
+      DOUBLE PRECISION A0,A1,A2,A3A,A3B,A3C,A3D
+      DOUBLE PRECISION A4,A5A,A5B,A5C,A5D,A5E,A5F
+      DOUBLE PRECISION A6,A8A,A8B,A8C,A9,A10,AXY
+      DOUBLE PRECISION B,C2,C5,C6,C7,C8,C9,C10,C11
+      DOUBLE PRECISION D2,D5,D6,D7,D8,D9,D10,D11
+      DOUBLE PRECISION E2,E5,E6,E7,E8,E9,E10,E11
+      DOUBLE PRECISION F2,F5,F6,F7,F8,F9,F10
+      DOUBLE PRECISION G2,G5,G6,G7,G8,G9,G10
+      DOUBLE PRECISION RX,RY
+      DOUBLE PRECISION tepxx,tepxy,tepyy
+
+      tepxx=x(1)*x(1)+x(2)*x(2)
+      tepxy=x(1)*y(1)+x(2)*y(2)
+      tepyy=y(1)*y(1)+y(2)*y(2)
+      
+      dyx(1)= y(1)-x(1)
+      dyx(2)= y(2)-x(2)
+      mdyx(1)=-dyx(1)
+      mdyx(2)=-dyx(2)
+      mx(1)=-x(1)
+      mx(2)=-x(2)
+      my(1)=-y(1)
+      my(2)=-y(2)
+      
+      rx= v-u+tepxx
+      ry= w-u+tepyy
+      axy= x(1)*y(2)-x(2)*y(1)
+      
+      A0= rx**2*tepyy - 2*rx*ry*tepxy + ry**2*tepxx
+      B= 4*u*axy**2 + A0
+      
+      A1= 2*(tepyy-tepxy)*(-2)*2*(-rx*tepxy+ry*tepxx)
+      A2= 2*(tepyy-tepxy)
+      C2= tepxy-tepxx
+      D2= 2*axy**2+(rx+ry)*tepxy-rx*tepyy-ry*tepxx
+      E2= -tepxy
+      F2= rx*tepyy-ry*tepxy
+      G2= tepxx
+      
+      A3a= 2*(-tepxy+tepxx) * 2*(rx*tepyy-ry*tepxy)     
+      A3b= (-4*axy**2 +2*(rx*tepyy-(rx+ry)*tepxy+ry*tepxx))
+     +  * 2*(-tepxy)
+      A3c= (-4*axy**2 +2*(rx*tepyy-(rx+ry)*tepxy+ry*tepxx))
+     &  * 2*(rx*tepyy-ry*tepxy)                  
+      A3d= 3*2*(-rx*tepxy+ry*tepxx)
+      
+      A4= A3c
+      A5a= 2*(-tepxy+tepxx)
+      A5b= -4*axy**2 +2*(rx*tepyy-(rx+ry)*tepxy+ry*tepxx)
+      A5c= 2*2*(-rx*tepxy+ ry*tepxx)
+      A5d= -2*(rx*tepyy-ry*tepxy)
+      A5e= -2*(-tepxy)
+      A5f= -2*(rx*tepyy-ry*tepxy)*(2)*(-rx*tepxy+ry*tepxx)
+      C5= tepxy-tepyy
+      D5= 2*axy**2+(rx+ry)*tepxy-rx*tepyy-ry*tepxx
+      E5= tepyy
+      F5= -tepxy
+      G5= ry*tepxx-rx*tepxy
+      
+      A6= -2*(rx*tepyy- ry*tepxy)
+      C6= tepxy-tepxx
+      D6= 2*axy**2+(rx+ry)*tepxy-rx*tepyy-ry*tepxx
+      E6= -tepxy
+      F6= rx*tepyy-ry*tepxy
+      G6= tepxx
+      
+      C7= tepxy-tepyy
+      D7= -(tepxy-tepxx)
+      E7= 2*axy**2+(rx+ry)*tepxy-rx*tepyy-ry*tepxx
+      F7= tepyy
+      G7= -tepxx
+      
+      A8a= -2*(-tepxy)
+      A8b= -2*(rx*tepyy-ry*tepxy)
+      A8c=  2*2*(-rx*tepxy+ry*tepxx)
+      C8= -2*tepxy+tepxx+tepyy
+      D8= -tepyy + tepxy
+      E8= rx*tepyy - ry*tepxy
+      F8= -tepxx + tepxy
+      G8= ry*tepxx - rx*tepxy
+      
+      A9= -2*(rx*tepyy-ry*tepxy)
+      C9= -2*tepxy+tepxx+tepyy
+      D9= -tepyy+tepxy  
+      E9= tepxy
+      F9= -(rx*tepyy-ry*tepxy)
+      G9= -tepxx
+      
+      A10= -2*(-rx*tepxy+ry*tepxx)
+      C10= -2*tepxy+tepxx+tepyy
+      D10= -tepyy
+      E10= -tepxx+tepxy
+      F10= tepxy
+      G10= -(ry*tepxx-rx*tepxy)
+      
+      C11= -(-2*tepxy+tepyy+tepxx)
+      D11= -tepyy
+      E11= tepxx
+      
+       Id3= -( (A1/B**2)*Id0(x,y,u,v,w) + (A2/B**2)*(
+     &  C2*Iz0(dyx,v,w) - D2*Iz1(mdyx,w,v)+
+     &  E2*Iz0(y,u,w) - F2*Iz1(my,w,u) + 
+     &  G2*Iz0(x,u,v) )      
+     &  -( ((A3a+A3b)*B - A3c*A3d)/B**3 * Id0(x,y,u,v,w) +
+     &  (A4/B**3)*(
+     &  C2*Iz0(dyx,v,w) - D2*Iz1(mdyx,w,v)+
+     &  E2*Iz0(y,u,w) - F2*Iz1(my,w,u)+
+     &  G2*Iz0(x,u,v) ) ) +
+     &  (A5a*B - A5b*A5c)/B**3*( A5d*Id0(x,y,u,v,w) + 
+     &  C5*Iz0(dyx,v,w) - D5*Iz1(dyx,v,w) + 
+     &  E5*Iz0(y,u,w) +
+     &  F5*Iz0(x,u,v) - G5*Iz1(mx,v,u) ) +
+     &  A5b/B**2 * ( (A5e - A5f/B) * Id0(x,y,u,v,w) + 
+     &  (A6/B)*(  
+     &  C6*Iz0(dyx,v,w) - D6*Iz1(mdyx,w,v) + 
+     &  E6*Iz0(y,u,w) - F6*Iz1(my,w,u)+
+     &  G6*Iz0(x,u,v)  ) 
+     &  - C7*Iz1(mdyx,w,v) + D7*Iz1(dyx,v,w) 
+     &  + E7*Iz2(mdyx,w,v)
+     &  - F7*Iz1(my,w,u) 
+     &  + G7*Iz1(mx,v,u) )+
+     &  (A8a*B - A8b*A8c)/B**3*(
+     &  C8*Iz0(dyx,v,w) +
+     &  D8*Iz0(y,u,w) - E8*Iz1(y,u,w) + 
+     &  F8*Iz0(x,u,v) - G8*Iz1(x,u,v)  ) +
+     &  A9/B**2*(  
+     &  - C9*Iz1(mdyx,w,v) 
+     &  - D9*Iz1(my,w,u) + E9*Iz1(y,u,w) - F9*Iz2(y,u,w) 
+     &  + G9*Iz1(x,u,v)  ) +
+     &  A10/B**2*(                  
+     &  - C10*Iz1(dyx,v,w) 
+     &  + D10*Iz1(y,u,w)
+     &  - E10*Iz1(mx,v,u) + F10*Iz1(x,u,v) - G10*Iz2(mx,v,u) )+
+     &  1/B*(
+     &  - C11*Iz2(dyx,v,w) 
+     &  - D11*Iz2(y,u,w)
+     &  + E11*Iz2(x,u,v)  ))
+      
+      END
+      
+C     --------------------------------------------------------
+      
+      DOUBLE PRECISION FUNCTION Iv0(x,y,z,u,v,w,w2)
+      IMPLICIT NONE
+      DOUBLE PRECISION Id0
+      EXTERNAL         Id0
+      DOUBLE PRECISION x(2),y(2),z(2),dyx(2),dzx(2)
+      DOUBLE PRECISION u,v,w,w2
+      DOUBLE PRECISION rx,ry,rz
+      DOUBLE PRECISION axy,ayz,azx
+      DOUBLE PRECISION A
+
+      rx= v-u+x(1)*x(1)+x(2)*x(2)
+      ry= w-u+y(1)*y(1)+y(2)*y(2)
+      rz= w2-u+z(1)*z(1)+z(2)*z(2)
+      
+      axy= x(1)*y(2)-x(2)*y(1)
+      ayz= y(1)*z(2)-y(2)*z(1)
+      azx= z(1)*x(2)-z(2)*x(1)
+      
+      dyx(1)= y(1)-x(1)
+      dyx(2)= y(2)-x(2)
+      dzx(1)= z(1)-x(1)
+      dzx(2)= z(2)-x(2)
+      
+      A= 1/(axy*rz + ayz*rx + azx*ry)
+      
+      Iv0= A*(-(axy+ayz+azx)*Id0(dyx,dzx,v,w,w2)
+     &  +axy*Id0(x,y,u,v,w)
+     &  +ayz*Id0(y,z,u,w,w2)
+     &  +azx*Id0(x,z,u,v,w2) )
+      
+      END
+     
+C     --------------------------------------------------------
+
+      DOUBLE PRECISION FUNCTION Iv1(x,y,z,u,v,w,w2)
+      IMPLICIT NONE
+      DOUBLE PRECISION Iv0,Id0,Id1
+      EXTERNAL         Iv0,Id0,Id1
+      
+      DOUBLE PRECISION x(2),y(2),z(2),dyx(2),dzx(2)
+      DOUBLE PRECISION u,v,w,w2
+      DOUBLE PRECISION A,B,RX,RY,RZ
+      DOUBLE PRECISION AXY,AYZ,AZX
+      
+      dyx(1)= y(1)-x(1)
+      dyx(2)= y(2)-x(2)
+      dzx(1)= z(1)-x(1)
+      dzx(2)= z(2)-x(2)
+      
+      rx= v-u+x(1)*x(1)+x(2)*x(2)
+      ry= w-u+y(1)*y(1)+y(2)*y(2)
+      rz= w2-u+z(1)*z(1)+z(2)*z(2)
+      
+      axy= x(1)*y(2)-x(2)*y(1)
+      ayz= y(1)*z(2)-y(2)*z(1)
+      azx= z(1)*x(2)-z(2)*x(1)
+      
+      A= 1/(axy*rz + ayz*rx + azx*ry)
+      B= axy+ayz+azx
+      Iv1=-((B*A**2)*(-(axy+ayz+azx)*Id0(dyx,dzx,v,w,w2)
+     &  +axy*Id0(x,y,u,v,w)
+     &  +ayz*Id0(y,z,u,w,w2)
+     &  +azx*Id0(x,z,u,v,w2) )+
+     &  A*( -axy*Id1(x,y,u,v,w)
+     &  -ayz*Id1(y,z,u,w,w2)
+     &  -azx*Id1(x,z,u,v,w2) ))
+      
+      END
+
+C     --------------------------------------------------------
+      
+      DOUBLE PRECISION FUNCTION Iv2(x,y,z,u,v,w,w2)
+      IMPLICIT NONE
+      DOUBLE PRECISION Iv1,Iv0,Id0,Id1,Id2
+      EXTERNAL         Iv1,Iv0,Id0,Id1,Id2
+      DOUBLE PRECISION x(2),mx(2),y(2),z(2),dyx(2)
+      DOUBLE PRECISION dzx(2)
+      DOUBLE PRECISION u,v,w,w2
+      DOUBLE PRECISION axy,azx,ayz
+      DOUBLE PRECISION rx,ry,rz,A,B
+      
+      mx(1)=-x(1)
+      mx(2)=-x(2)
+      dyx(1)= y(1)-x(1)
+      dyx(2)= y(2)-x(2)
+      dzx(1)= z(1)-x(1)
+      dzx(2)= z(2)-x(2)
+      
+      rx= v-u+x(1)*x(1)+x(2)*x(2)
+      ry= w-u+y(1)*y(1)+y(2)*y(2)
+      rz= w2-u+z(1)*z(1)+z(2)*z(2)
+      
+      axy= x(1)*y(2)-x(2)*y(1)
+      azx= z(1)*x(2)-z(2)*x(1)
+      ayz= y(1)*z(2)-y(2)*z(1)
+      
+      A= 1/(axy*rz + ayz*rx + azx*ry)
+      B= axy+ayz+azx
+      
+      Iv2=(-B*2D0*ayz*A**3)*(
+     &  -B*Id0(dyx,dzx,v,w,w2)
+     &  +axy*Id0(x,y,u,v,w)
+     &  +ayz*Id0(y,z,u,w,w2)
+     &  +azx*Id0(x,z,u,v,w2) )+
+     &  (B*A**2)*(B*Id1(dyx,dzx,v,w,w2)
+     &  -axy*Id1(mx,dyx,v,u,w)
+     &  -azx*Id1(mx,dzx,v,u,w2) )+
+     &  (-ayz*A**2)*(-axy*Id1(x,y,u,v,w)
+     &  -ayz*Id1(y,z,u,w,w2)
+     &  -azx*Id1(x,z,u,v,w2) )+
+     &  A*(axy*Id2(x,y,u,v,w)
+     &  +azx*Id2(x,z,u,v,w2) )
+      
+      END
diff --git a/EPEMGEN/dtrint.f b/EPEMGEN/dtrint.f
new file mode 100644 (file)
index 0000000..438e72b
--- /dev/null
@@ -0,0 +1,229 @@
+      FUNCTION DTRINT(F,NSD,NPT,EPS,X1,Y1,X2,Y2,X3,Y3)
+C     INTEGRATION OVER A TRIANGLE USING A 7-, 25- OR 64-POINT FORMULA,
+C     WITH OR WITHOUT SUBDIVISION OF THE TRIANGLE.
+C     CERNLIB program D105.
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      CHARACTER NAME*(*)
+      CHARACTER*80 ERRTXT
+      PARAMETER (NAME = 'DTRINT')
+      PARAMETER (Z1 = 1, HALF = Z1/2, KMX = 35)
+
+      DIMENSION R(KMX),XP1(KMX),XP2(KMX),XP3(KMX),YP1(KMX),YP2(KMX)
+      DIMENSION YP3(KMX)
+      DIMENSION U(8,3),V(32,3),W(32,3),UU(5,2),WW(5,2)
+      DIMENSION G1(32),G2(32)
+      DIMENSION JZ(3),JZ0(3),JZ1(3),IZ(3,3)
+
+      DATA JZ /2,5,8/, JZ0 /2,10,32/, JZ1 /3,5,0/
+      DATA ((IZ(I,M),I=1,3),M=1,3) /1,0,1, 5,5,1, 8,24,8/
+
+      DATA (U(J,1),J=1,2)
+     1/0.89871 34926 76543 661D0, 0.52985 79358 94884 910D0/
+
+      DATA (UU(J,1),J=1,3)
+     1/0.20257 30146 46912 678D0, 0.94028 41282 10230 180D0,
+     2 0.66666 66666 66666 667D0/
+
+      DATA (V(J,1),J=1,2)
+     1/0.69614 04780 29630 984D0,-0.41042 61923 15345 269D0/
+
+      DATA (W(J,1),J=1,2)
+     1/0.12593 91805 44827 153D0, 0.13239 41527 88506 181D0/
+
+      DATA (WW(J,1),J=1,3)
+     1/0.12593 91805 44827 153D0, 0.13239 41527 88506 181D0,
+     2 0.22500 00000 00000 000D0/
+
+      DATA (U(J,2),J=1,5)
+     1/0.09853 50857 98826 426D0, 0.30453 57266 46363 905D0,
+     2 0.56202 51897 52613 856D0, 0.80198 65821 26391 827D0,
+     3 0.96019 01429 48531 258D0/
+
+      DATA (UU(J,2),J=1,5)
+     1/0.09853 50857 98826 426D0, 0.30453 57266 46363 905D0,
+     2 0.56202 51897 52613 856D0, 0.80198 65821 26391 827D0,
+     3 0.96019 01429 48531 258D0/
+
+      DATA (V(J,2),J=1,10)
+     1/0.08929 05088 68733 569D0, 0.27596 41378 55221 135D0,
+     2 0.50929 58998 63672 021D0, 0.72674 40774 36169 444D0,
+     3 0.87010 49558 08923 811D0, 0.05305 81196 71298 357D0,
+     4 0.16398 31426 29800 463D0, 0.30263 33161 88105 613D0,
+     5 0.43184 51615 91612 961D0, 0.51703 29238 43772 854D0/
+
+      DATA (W(J,2),J=1,10)
+     1/0.00373 11043 33755 67687D0, 0.01751 09983 64327 66347D0,
+     2 0.03468 30128 62731 40026D0, 0.03960 81662 64094 70756D0,
+     3 0.02293 01607 03185 09559D0, 0.00753 74033 90655 24076D0,
+     4 0.03537 49042 20966 93175D0, 0.07006 50090 06743 44063D0,
+     5 0.08001 45747 72320 84819D0, 0.04632 24438 58996 77269D0/
+
+      DATA (WW(J,2),J=1,5)
+     1/0.00895 88135 94562 71712D0, 0.04204 59349 74644 15024D0,
+     2 0.08327 79304 30389 93562D0, 0.09510 37941 15908 01948D0,
+     3 0.05505 79713 28939 62198D0/
+
+      DATA (U(J,3),J=1,8)
+     1/0.04463 39552 89969 851D0, 0.14436 62570 42145 571D0,
+     2 0.28682 47571 44430 519D0, 0.45481 33151 96573 351D0,
+     3 0.62806 78354 16727 698D0, 0.78569 15206 04369 242D0,
+     4 0.90867 63921 00206 044D0, 0.98222 00848 52636 548D0/
+
+      DATA (V(J,3),J=1,32)
+     1/0.04286 15345 20322 596D0, 0.13863 34522 58088 400D0,
+     2 0.27543 49048 78165 863D0, 0.43675 26131 83286 138D0,
+     3 0.60312 71715 43047 645D0, 0.75449 15975 72500 770D0,
+     4 0.87259 27221 72605 828D0, 0.94321 59843 32136 212D0,
+     5 0.03555 83759 33897 592D0, 0.11501 17574 55156 297D0,
+     6 0.22850 36689 09272 431D0, 0.36233 45216 98467 603D0,
+     7 0.50036 05900 18245 933D0, 0.62593 40960 53638 777D0,
+     8 0.72391 20204 03394 633D0, 0.78250 18150 44463 514D0,
+     9 0.02345 65900 87635 536D0, 0.07586 91469 73958 963D0,
+     A 0.15073 57058 45778 370D0, 0.23901 91375 97290 126D0,
+     B 0.33007 00031 37485 188D0, 0.41290 63582 74039 218D0,
+     C 0.47753 88941 74496 369D0, 0.51618 84882 60827 229D0,
+     D 0.00818 74136 31782 437D0, 0.02648 17727 48961 059D0,
+     E 0.05261 35967 85690 189D0, 0.08342 85178 75344 723D0,
+     F 0.11520 93988 52684 066D0, 0.14412 30431 93925 944D0,
+     G 0.16668 27291 29138 200D0, 0.18017 31901 16990 201D0/
+
+      DATA (W(J,3),J=1,32)
+     1/0.00033 35674 06495 41982D0, 0.00180 62109 19037 15084D0,
+     2 0.00459 97558 03491 41419D0, 0.00801 72595 31391 49525D0,
+     3 0.01073 50189 73158 61631D0, 0.01138 87974 04616 51588D0,
+     4 0.00922 38453 90918 29977D0, 0.00450 98127 16079 21752D0,
+     5 0.00073 27880 81649 19485D0, 0.00396 79231 50289 07586D0,
+     6 0.01010 48428 76312 33624D0, 0.01761 24888 63394 88637D0,
+     7 0.02358 29214 92410 93797D0, 0.02501 91560 68339 84265D0,
+     8 0.02026 31427 34638 24614D0, 0.00990 72539 59652 71520D0,
+     9 0.00103 37234 54873 38862D0, 0.00559 74371 44935 08001D0,
+     A 0.01425 46165 12792 75399D0, 0.02484 54407 11607 46855D0,
+     B 0.03326 77614 32852 32482D0, 0.03529 38169 93822 26192D0,
+     C 0.02858 46432 80634 70277D0, 0.01397 58834 07425 66299D0,
+     D 0.00119 51124 99230 79556D0, 0.00647 13314 41724 90169D0,
+     E 0.01648 01043 12102 39366D0, 0.02872 44103 85925 30995D0,
+     F 0.03846 16575 37508 13161D0, 0.04080 40290 04108 74571D0,
+     G 0.03304 73922 30182 37768D0, 0.01615 78542 78398 33562D0/
+
+      IF(NPT .EQ. 7) THEN
+       M=1
+      ELSE IF(NPT .EQ. 25) THEN
+       M=2
+      ELSE IF(NPT .EQ. 64) THEN
+       M=3
+      ELSE
+       H=0
+       WRITE(ERRTXT,101) NPT
+       CALL MTLPRT(NAME,'D105.1',ERRTXT)
+       GO TO 99
+      END IF
+      K=0
+      C1=X1
+      D1=Y1
+      C2=X2
+      D2=Y2
+      C3=X3
+      D3=Y3
+      A11=HALF*(C2+C3)-C1
+      A12=HALF*(C3-C2)
+      A21=HALF*(D2+D3)-D1
+      A22=HALF*(D3-D2)
+      DO 1 J = 1,JZ(M)
+      G1(J)=A11*U(J,M)+C1
+      G2(J)=A21*U(J,M)+D1
+      DO 1 I = IZ(1,M),IZ(2,M),IZ(3,M)
+      G1(J+I)=G1(J)
+    1 G2(J+I)=G2(J)
+      S=0
+      DO 2 J = 1,JZ0(M)
+      H1=A12*V(J,M)
+      H2=A22*V(J,M)
+    2 S=S+W(J,M)*(F(G1(J)+H1,G2(J)+H2)+F(G1(J)-H1,G2(J)-H2))
+      DO 3 J = 1,JZ1(M)
+    3 S=S+WW(J,M)*F(A11*UU(J,M)+C1,A21*UU(J,M)+D1)
+      S=HALF*ABS(C1*(D2-D3)+C2*(D3-D1)+C3*(D1-D2))*S
+      H=S
+      IF(NSD .EQ. 0) GO TO 99
+      H=0
+
+   10 K=K+1
+      SUM0=S
+      U1=C1
+      V1=D1
+      U2=C2
+      V2=D2
+      U3=C3
+      V3=D3
+   11 C1=HALF*(U2+U3)
+      D1=HALF*(V2+V3)
+      C2=U1
+      D2=V1
+      C3=U2
+      D3=V2
+      XP1(K)=C1
+      YP1(K)=D1
+      XP2(K)=C2
+      YP2(K)=D2
+      XP3(K)=U3
+      YP3(K)=V3
+      A11=HALF*(C2+U3)-C1
+      A12=HALF*(U3-C2)
+      A21=HALF*(D2+V3)-D1
+      A22=HALF*(V3-D2)
+      DO 4 J = 1,JZ(M)
+      G1(J)=A11*U(J,M)+C1
+      G2(J)=A21*U(J,M)+D1
+      DO 4 I = IZ(1,M),IZ(2,M),IZ(3,M)
+      G1(J+I)=G1(J)
+    4 G2(J+I)=G2(J)
+      S=0
+      DO 5 J = 1,JZ0(M)
+      H1=A12*V(J,M)
+      H2=A22*V(J,M)
+    5 S=S+W(J,M)*(F(G1(J)+H1,G2(J)+H2)+F(G1(J)-H1,G2(J)-H2))
+      DO 6 J = 1,JZ1(M)
+    6 S=S+WW(J,M)*F(A11*UU(J,M)+C1,A21*UU(J,M)+D1)
+      S=HALF*ABS(C1*(D2-V3)+C2*(V3-D1)+U3*(D1-D2))*S
+      R(K)=S
+      A11=HALF*(C2+C3)-C1
+      A12=HALF*(C3-C2)
+      A21=HALF*(D2+D3)-D1
+      A22=HALF*(D3-D2)
+      DO 7 J = 1,JZ(M)
+      G1(J)=A11*U(J,M)+C1
+      G2(J)=A21*U(J,M)+D1
+      DO 7 I = IZ(1,M),IZ(2,M),IZ(3,M)
+      G1(J+I)=G1(J)
+    7 G2(J+I)=G2(J)
+      S=0
+      DO 8 J = 1,JZ0(M)
+      H1=A12*V(J,M)
+      H2=A22*V(J,M)
+    8 S=S+W(J,M)*(F(G1(J)+H1,G2(J)+H2)+F(G1(J)-H1,G2(J)-H2))
+      DO 9 J = 1,JZ1(M)
+    9 S=S+WW(J,M)*F(A11*UU(J,M)+C1,A21*UU(J,M)+D1)
+      S=HALF*ABS(C1*(D2-D3)+C2*(D3-D1)+C3*(D1-D2))*S
+      SUM=S+R(K)
+      IF(ABS(SUM0-SUM) .GT. EPS*(1+ABS(SUM))) THEN
+       IF(K .LT. KMX) GO TO 10
+       H=0
+       CALL MTLPRT(NAME,'D105.2','TOO HIGH ACCURACY REQUIRED')
+       GO TO 99
+      ELSE
+       H=H+SUM
+       K=K-1
+       IF(K .LE. 0)  GO TO 99
+       U1=XP1(K)
+       V1=YP1(K)
+       U2=XP2(K)
+       V2=YP2(K)
+       U3=XP3(K)
+       V3=YP3(K)
+       SUM0=R(K)
+       GO TO 11
+      END IF
+   99 DTRINT=H
+      RETURN
+  101 FORMAT('INCORRECT NUMBER OF POINTS =',I5)
+      END
diff --git a/EPEMGEN/epemgen.f b/EPEMGEN/epemgen.f
new file mode 100644 (file)
index 0000000..4608424
--- /dev/null
@@ -0,0 +1,515 @@
+C **********************************************************************
+C                                                  Version 2.0
+C     Generator of e+e- pairs produced in            1.10.2002
+C                  PbPb collisions at LHC             
+C                                            
+C     Copyright (c) 2002
+C     Authors:          Kai Hencken         <k.hencken@unibas.ch>
+C                       Yuri Kharlov        <Yuri.Kharlov@cern.ch>
+C                       Serguei Sadovsky     <sadovsky@mx.ihep.su>,
+C                                            <Serguei.Sadovski@cern.ch>
+C
+C Permission to use, copy and distribute this software and its
+C documentation strictly for non-commercial purposes is hereby granted
+C without fee, provided that the above copyright notice appears in all
+C copies and that both the copyright notice and this permission notice
+C appear in the supporting documentation. The authors make no claims
+C about the suitability of this software for any purpose. It is
+C provided "as is" without express or implied warranty.
+C Any change of the code should be submitted to the authors
+C
+C     Long write up:   K.Hencken,Yu.Kharlov,S.Sadovsky
+C                      Internal ALICE Note 2002-27 
+C                      
+C **********************************************************************
+      subroutine ee_init(Ymin,Ymax,PTmin,PTmax)
+C-----------------------------------------------------------------------
+C    Generator initialisation
+C
+C    Input variables:  Ymin      - minimal value of rapidity  \
+C                      Ymax      - maximal value of rapidity   | of kinematics 
+C                      PTmin     - Pt minimum in MeV/c        | range
+C                      PTmax     - Pt maximum in MeV/c        /
+C-----------------------------------------------------------------------
+      implicit real*8 (A-H,O-Z)
+      external DsdYpY,DsdYmY,DsdXpX,DsdXmX,DsdPhi,DsdXX,DsdYY
+
+      common/parPh / parPhi(7)
+      common/parYp / parYpY(6)
+      common/parYm / parYmY(6)
+      common/parXp / parXpX(5)
+      common/parXm / parXmX(4)
+
+      common/eevent/ Xsect2,Dsect2, Xsecttot,Dsecttot, Nevnt
+      common/eepars/ Xmin,Xmax,YpYmin,YpYmax,WYpYmax,YmYmin,YmYmax,
+     &               Ymed1,Ymed2,sgmY1,sgmY2,XYsect,
+     &               Gaus1,Gaus2,Gauss,Exp1,Exp2,Exp3,
+     &               XmXmin,XmXmax,XpXmin,XpXmax,Exmx1,Exmx2,Xmed,
+     &               Sgm1,AnorX,Icase
+
+      real*8 mass  
+
+      data pi     / 3.141 592 653 589 793 238 462 643d00 /
+
+C
+C  Exact differential Cross section:
+C - - - - - - - - - - - - - - - - - - -
+C
+      gm = 2750.0d0             ! Pb gamma factor at the LHC
+      mass = 0.5109991d0        ! electron mass 
+      call initdiffcross (gm,mass)
+C     
+C  Cross sections:
+C - - - - - - - - - -         
+      Nevnt  = 0
+      Xsect2 = 0.
+      Dsect2 = 0.
+      
+      Nsd = 1
+      Npt = 25                  !  7,25,64
+      Eps = 0.00005
+
+      Xmin=    Dlog10(Ptmin)
+      Xmax=    Dlog10(Ptmax)
+      write(*,*) ' Kinematical limits:'
+      write(*,*) '          Xmin,Xmax=',Xmin,Xmax 
+      write(*,*) '          Ymin,Ymax=',Ymin,Ymax
+      write(*,*) 
+
+      Xsec1  = Dtrint(DsdXX,Nsd,Npt,Eps,Xmin,Xmin,Xmin,Xmax,Xmax,Xmin)
+      Xsec2  = Dtrint(DsdXX,Nsd,Npt,Eps,Xmax,Xmax,Xmin,Xmax,Xmax,Xmin)
+      XsecX  = Xsec1+Xsec2
+      
+      Ysec1  = Dtrint(DsdYY,Nsd,Npt,Eps,Ymin,Ymin,Ymin,Ymax,Ymax,Ymin)
+      Ysec2  = Dtrint(DsdYY,Nsd,Npt,Eps,Ymax,Ymax,Ymin,Ymax,Ymax,Ymin)
+      XsecY  = Ysec1+Ysec2
+
+      IF (Xsec1*Xsec2.le.0.or.Ysec1*Ysec2.le.0.) then
+         write(*,*) ' Error: insufficient accuracy of XY-cross sections'  
+         write(*,*) ' Xsec1,Xsec2,Xsec=', Xsec1,Xsec2,XsecX
+         write(*,*) ' Ysec1,Ysec2,Ysec=', Ysec1,Ysec2,XsecY
+      endif
+      
+      XYsect = XsecX*XsecY
+      write(*,*) ' Xsections: Xsec1,Xsec2,XsecX=', Xsec1,Xsec2,XsecX
+      write(*,*) '            Ysec1,Ysec2,XsecY=', Ysec1,Ysec2,XsecY
+C-    write(*,*) '                  XsecX*YsecY=', XYsect
+C
+      XYsect = 2371.5239*(82./137.035)**4*XYsect ! Normalization factor 
+      write(*,*) ' Normalized:      XsecX*YsecY=', XYsect
+      write(*,*) 
+C      
+C-                XsecPhi  = Dgauss(DsdPhi,0.0d0,  pi, eps)
+C
+C  Rapidity initialisation:
+C - - - - - - - - - - - - - -  
+      if (Ymin.ge.Ymax) then
+         write(*,*) 'Wrong values of Ymin,Ymax:',Ymin,Ymax
+         stop
+      endif      
+CYpY-
+C
+      YpYmin = 2.*Ymin
+      YpYmax = 2.*Ymax
+      Yp0    = 0.
+      if (YpYmin*YpYmax.gt.0.) then 
+         if(YpYmin.gt.0.) Yp0 = YpYmin            
+         if(YpYmin.le.0.) Yp0 = YpYmax
+      endif
+      WYpYmax = DsdYpY(Yp0)
+C     
+C-    write(*,*)' YpY: YpYmin,YpYmax,WYpYmax =',YpYmin,YpYmax,WYpYmax
+
+CYmY-
+C
+      YmYmin = Ymin-Ymax
+      YmYmax = Ymax-Ymin
+C-    write(*,*) ' YmY:  YmYmin,YmYmax=',YmYmin,YmYmax
+C     
+      Ymed1  = 0.18
+      Ymed2  = 4.00
+      
+      sgmY1  = 1./dsqrt(2.*parYmY(2))
+      sgmY2  = 1./dsqrt(2.*parYmY(4))
+C     
+      eps = 0.000001
+      Gaus1 = Dgauss(DsdYmY ,0.0d0, Ymed1 , eps)
+      Gaus2 = Dgauss(DsdYmY ,Ymed1, Ymed2 , eps)
+      Expon = Dgauss(DsdYmY ,Ymed2, YmYmax, eps)
+      Summ  = Gaus1+Gaus2+Expon
+      Gaus2 =(Gaus1+Gaus2)/Summ
+      Gaus1 = Gaus1/Summ
+      Expon = 1.
+C     
+C     Phi initialisation:
+C   - - - - - - - - - - - 
+      Exp0 = parPhi(1)*pi
+      Exp1 = parPhi(2)*(1.-dexp(-parPhi(3)*pi))/parPhi(3)
+      Exp2 = parPhi(4)*(1.-dexp(-parPhi(5)*pi))/parPhi(5)
+      Exp3 = parPhi(6)*(1.-dexp(-parPhi(7)*pi))/parPhi(7)
+      Summ = Exp0+Exp1+Exp2+Exp3
+C     
+      Exp0 = (Exp0+Exp1+Exp2+Exp3)/Summ
+      Exp1 = (Exp1+Exp2+Exp3)/Summ   
+      Exp2 = (Exp2+Exp3)/Summ 
+      Exp3 =  Exp3/Summ
+C-    write(*,*) 'Exp0,Exp1,Exp2,Exp3=',Exp0,Exp1,Exp2,Exp3
+C
+C  Pt initialisation:
+C - - - - - - - - - - -        
+C
+      XmXmin = Xmin-Xmax
+      XmXmax = Xmax-Xmin
+C     
+      XpXmin = 2.*Xmin
+      XpXmax = 2.*Xmax
+C     
+CXmX-
+      Exmx1 = parXmX(1)*(1.-dexp(-parXmX(2)*dabs(XmXmax)))/parXmX(2)
+      Exmx2 = parXmX(3)*(1.-dexp(-parXmX(4)*dabs(XmXmax)))/parXmX(4)
+      Exmx1 = Exmx1/(Exmx1+Exmx2)
+      Exmx2 = 1. 
+C     
+C-    write(*,*) ' XpX:  XpXmin,XpXmax=',XpXmin,XpXmax
+C-    write(*,*) ' XmX:  XmXmin,XmXmax=',XmXmin,XmXmax
+C-    write(*,*) ' XmX:  Exmx1 ,Exmx2 =',Exmx1 ,Exmx2
+C     
+CXpX-
+      Icase = 2                 !   Gausses and Exponent 
+      Xmed  = 0.6
+      if (XpXmax.lt.Xmed) then
+         Icase = 1              !   Gauss only
+         Xmed  = XpXmax
+      endif            
+C
+      if (XpXmin.gt.Xmed) then
+         Icase = 3              !   Exponent only
+         Xmed  = XpXmin
+      endif
+C      
+      if (Icase.eq.2)     then  !   Gauss and Exponent
+         eps = 0.000001
+         Gauss = Dgauss(DsdXpX ,XpXmin, Xmed  , eps)
+         Expon = Dgauss(DsdXpX ,Xmed  , XpXmax, eps)
+         Summ  = Gauss+Expon
+         Gauss = Gauss/Summ 
+         Expon = Expon/Summ
+C-       write(*,*) 'XpX, Case 2: Gauss,Expon=',Gauss,Expon     
+      endif                     ! Icase = 2
+C
+C
+      if (Icase.eq.1.or.Icase.eq.2) then !   Gausses only and Icase=2
+         
+         Sgm1 =       1./dsqrt(2.*parXpX(2))
+         
+C-       write(*,*) 'XpX, Case 1:  Sgm1=',Sgm1        
+      endif
+C      
+      if (Icase.eq.3.or.Icase.eq.2) then !    Exponent only and Icase=2   
+C
+         AnorX = 1.-dexp(-parXpX(5)*(XpXmax-Xmed))
+C     
+C-       write(*,*) 'XpX, Case 3:  AnorX=',AnorX
+      endif
+         
+      return
+      end
+*
+*=======================================================================
+      subroutine ee_event(Ymin,Ymax,PTmin,PTmax,
+     +                    Ye,Yp,Xe,Xp,Phi,Wtm2)
+C------------------------------------------------------------------------------
+C    Produce one event 
+C
+C    Input variables:  Ymin     - minimal value of rapidity  \
+C                      Ymax     - maximal value of rapidity   | of kinematics 
+C                      PTmin    - Pt minimum in MeV/c        | range
+C                      PTmax    - Pt maximum in MeV/c       / 
+C    Output variables: Ye,Yp    - rapidity of produced e- or e+
+C                      Xe,Xp    - log10(pt) of produced e- or e+, pt in MeV/c
+C                      Phi      - azymuth angle between e- and e+
+C                      Wtm2     - event weight. The sum of these event 
+C                                 weights, divided by the total number 
+C                                 of generated events, gives the integral 
+C                                 cross section of the process of e+e- pair 
+C                                production in the above mentioned kinematics
+C                                 range.
+C                                Sum of the selected event weights, divided 
+C                                 by the total number of generated events, 
+C                                 gives the integral cross section corresponded 
+C                                 to the set of selected events
+C------------------------------------------------------------------------------
+
+      implicit real*8 (A-H,O-Z)
+      common/parPh / parPhi(7)
+      common/parYp / parYpY(6)
+      common/parYm / parYmY(6)
+      common/parXp / parXpX(5)
+      common/parXm / parXmX(4)
+      common/eevent/ Xsect2,Dsect2, Xsecttot,Dsecttot, Nevnt
+      common/eepars/ Xmin,Xmax,YpYmin,YpYmax,WYpYmax,YmYmin,YmYmax,
+     &               Ymed1,Ymed2,sgmY1,sgmY2,XYsect,
+     &               Gaus1,Gaus2,Gauss,Exp1,Exp2,Exp3,
+     &               XmXmin,XmXmax,XpXmin,XpXmax,Exmx1,Exmx2,Xmed,
+     &               Sgm1,AnorX,Icase
+      external DsdYpY,DsdYmY,DsdXpX,DsdXmX,DsdPhi 
+      data pi     / 3.141 592 653 589 793 238 462 643d00 /
+
+C  Rapidity distributions: 
+C
+ 10   YpY = YpYmin + (YpYmax-YpYmin)*pyr(0)
+      Wt  = DsdYpY(YpY)
+      if (Wt.lt.pyr(0)*WYpYmax) go to 10
+C
+CYmY-
+      r1 = pyr(0)
+      if  (r1.lt.Gaus1) then
+ 13      call rnorml(rn1)
+         YmY = sgmY1*rn1
+         if (dabs(YmY).gt.Ymed1) go to 13
+      else if  (r1.lt.Gaus2) then
+ 15      call rnorml(rn2)
+         YmY = sgmY2*rn2
+         if (dabs(YmY).lt.Ymed1) go to 15
+         if (dabs(YmY).gt.Ymed2) go to 15
+      else
+ 20      r2 = pyr(0)
+         YmY=-Dlog(r2)/parYmY(6) + Ymed2
+         if (YmY.gt.YmYmax) go to 20    
+         r2 = pyr(0)
+         if (r2.lt.0.5) YmY =-YmY
+      endif
+C
+      Ye = 0.5*(YpY+YmY)
+      Yp = 0.5*(YpY-YmY)
+      
+      if (Ye.lt.Ymin.or.Ye.gt.Ymax) go to 10
+      if (Yp.lt.Ymin.or.Yp.gt.Ymax) go to 10      
+
+C      
+C Azimuthal angle:       
+C 
+ 30   r1 = pyr(0)
+      if (r1.lt.Exp3) then
+ 31      r2 = pyr(0)
+         Phi=-Dlog(r2)/parPhi(7)
+         if (Phi.gt.pi) go to 31
+         go to 50            
+      else if (r1.lt.Exp2) then
+ 32      r2 = pyr(0)
+         Phi=-Dlog(r2)/parPhi(5)
+         if (Phi.gt.pi) go to 32
+         go to 50  
+      else if (r1.lt.Exp1) then 
+ 33      r2 = pyr(0)
+         Phi=-Dlog(r2)/parPhi(3)
+         if (Phi.gt.pi) go to 33
+         go to 50
+      else
+         Phi = pi*pyr(0)
+      endif
+C     
+ 50   if (pyr(0).gt.0.5) Phi =-Phi
+      Phi = Phi+pi
+      if (Phi.lt.0.03.or.Phi.gt.2.*pi-0.03) go to 30 
+C
+C Transverse momentums: 
+C     
+CXpX-
+
+  60  Jcase = Icase             !       Jcase = Icase, if Icase.ne.2  
+C
+      if (Icase.eq.2) then      !       Gausses and Exponent       
+         Jcase = 3
+         if (pyr(0).lt.Gauss) Jcase = 1
+      endif                     ! of Icase 2
+C           
+      if (Jcase.eq.1) then      !       Gauss only
+ 65      call rnorml(rn1)
+         XpX=Sgm1*rn1-parXpX(3)
+         if (XpX.lt.XpXmin.or.XpX.gt.Xmed) go to 65
+      endif                     ! of Jcase 1         
+C
+      if (Jcase.eq.3) then      !       Exponent only
+ 70      r1  = pyr(0)
+         XpX =-Dlog( 1.- r1*AnorX)/parXpX(5)+Xmed            
+      endif                     ! of Jcase 3 
+
+CXmX-
+      r1 = pyr(0)
+      if (r1.lt.Exmx1) then
+ 81      r2 = pyr(0)
+         XmX=-Dlog(r2)/parXmX(2)
+         if (XmX.gt.dabs(XmXmax)) go to 81           
+      else 
+ 82      r2 = pyr(0)
+         XmX=-Dlog(r2)/parXmX(4)
+         if (XmX.gt.dabs(XmXmax)) go to 82  
+      endif
+C
+ 85   if (pyr(0).gt.0.5) XmX =-XmX
+C
+      Xe = 0.5*(XpX+XmX)
+      Xp = 0.5*(XpX-XmX)
+      
+      if (Xe.lt.Xmin.or.Xe.gt.Xmax) go to 60
+      if (Xp.lt.Xmin.or.Xp.gt.Xmax) go to 60 
+C
+C --- Good event: 
+C
+      Pte = 10.d0**Xe
+      Ptp = 10.d0**Xp
+C
+C Event weight if the events would generate uniformly
+C      
+      Wt = DsdYpY(YpY)*DsdYmY(YmY)*DsdXpX(Xpx)*DsdXmX(XmX)*DsdPhi(Phi)
+      Wt = 2.2706950/Wt
+C
+      Wtm2=(Pte*Ptp)*Wt          ! Transformation factor
+C
+C --- Exact diff. cross section (Adrian Alscher 1997, Kai Hencken, May '98):  
+C
+      call Diffcross(Ptp,Yp,Pte,Ye,Phi,Dsigma)
+      Wtm2= XYsect*Dsigma*(Pte*Ptp)*Wtm2
+C   
+      Nevnt  = Nevnt  + 1
+      Xsect2 = Xsect2 + Wtm2
+      Dsect2 = Dsect2 + Wtm2*Wtm2
+C
+C     Calculate cross section and its error accumulated so far
+C
+      Xsecttot = Xsect2/Nevnt
+      Dsecttot = dsqrt(Dsect2)/Nevnt
+C
+      return
+      end
+c      
+C=================================================================
+      Double precision function DsdYpY(Y)    
+c
+c     Y = Yp+Ye
+c - - - - - - - - - - - - - - - - - - -
+      implicit real*8 (A-H,O-Z)
+      common/parYp/   parYpY(6)
+C-    data parYpY /  -48.584, 0.11403E-01, 40.649, 0.38920E-04,
+C-   +                                     40.283, 0.19238E-01 / 
+      data parYpY /  -4.8584, 0.11403E-01, 4.0649, 0.38920E-04,
+     +                                     4.0283, 0.19238E-01 / 
+c       
+      DsdYpY = parYpY(1)*dexp(-parYpY(2)*Y**2)
+     +       + parYpY(3)*dexp(-parYpY(4)*Y**4)
+     +       + parYpY(5)*dexp(-parYpY(6)*Y**2)
+      return
+      end
+C 
+      Double precision function DsdYmY(Y)    
+c
+c     Y = Yp-Ye
+c - - - - - - - - - - - - - - - - - - -
+      implicit real*8 (A-H,O-Z)
+      common/parYm/   parYmY(6)
+C-    data parYmY /   180.,  8., 174.38, 0.17867, 2418.8, 1.3097 /
+      data parYmY /   1.80,  8., 1.7438, 0.17867, 24.188, 1.3097 /
+C      
+      if (abs(Y).lt.0.18) then
+C1-      if (abs(Y).lt.0.00) then
+         DsdYmY  = parYmY(1)*dexp(-parYmY(2)*abs(Y)**2)      
+      else if (abs(Y).lt.4.00) then
+         DsdYmY  = parYmY(3)*dexp(-parYmY(4)*abs(Y)**2)
+      else
+         DsdYmY  = parYmY(5)*dexp(-parYmY(6)*abs(Y))
+      endif
+         
+      return
+      end
+C      
+      Double precision function DsdYY(Yp,Ye)    
+c - - - - - - - - - - - - - - - - - - - - - - 
+      implicit real*8 (A-H,O-Z)
+c
+      YpY   = Yp+Ye
+      YmY   = Yp-Ye
+      DsdYY = DsdYpY(YpY)*DsdYmY(YmY)
+      return
+      end
+C 
+C====================================================================
+C 
+      Double precision function DsdXpX(X)    
+c
+c     X = Xp+Xe
+c - - - - - - - - - - - - - - - - - - -
+      implicit real*8 (A-H,O-Z)
+      common/parXp/   parXpX(5)
+
+C-    data parXpX /   83.668, 1.2004, 0.47225, 96.951, 2.3814 /
+      data parXpX /   8.3668, 1.2004, 0.47225, 9.6951, 2.3814 /
+c       
+      if (X.lt.0.6)then 
+         DsdXpX = parXpX(1)*dexp(-parXpX(2)*(X+parXpX(3))**2)         
+      else
+         DsdXpX = parXpX(4)*dexp(-parXpX(5)*X)
+      endif
+      return
+      end
+C 
+      Double precision function DsdXmX(X)    
+c
+c     X = Xp-Xe
+c - - - - - - - - - - - - - - - - - - -
+      implicit real*8 (A-H,O-Z)
+      common/parXm/   parXmX(4)
+c-    data parXmX /   950.38, 39.040, 194.92, 3.5660 /
+      data parXmX /   9.5038, 39.040, 1.9492, 3.5660 /
+c      
+      DsdXmX = parXmX(1)*dexp(-parXmX(2)*abs(X))   
+     +       + parXmX(3)*dexp(-parXmX(4)*abs(X))
+         
+      return
+      end        
+C      
+      Double precision function DsdXX(Xp,Xe)    
+c - - - - - - - - - - - - - - - - - - - - - - 
+      implicit real*8 (A-H,O-Z)
+c
+      XpX   = Xp+Xe
+      XmX   = Xp-Xe
+      DsdXX = DsdXpX(XpX)*DsdXmX(XmX)
+      return
+      end
+C====================================================================
+C
+      Double precision function DsdPhi(Phi)  ! Phi distribution 
+c
+c     Phi - athimutal angle between e+ and e- (in radians)
+c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
+      implicit real*8 (A-H,O-Z)
+      common/parPh/ parPhi(7)
+C                                                         !  For 10**7 events
+      data parPhi / 5.4825, 226.00, 11.602,      68.173,  !  accuracy= 0.87 % 
+     +                     1.9509, 0.11864E+04, 109.61 / !  N(Wt>20)= 573  
+C
+      data pi     / 3.141 592 653 589 793 238 462 643d00/
+C
+      DsdPhi  = parPhi(1)
+     +        + parPhi(2)*dexp(-parPhi(3)*dabs(Phi-pi)) 
+     +        + parPhi(4)*dexp(-parPhi(5)*dabs(Phi-pi))
+     +        + parPhi(6)*dexp(-parPhi(7)*dabs(Phi-pi))            
+      return
+      end
+C====================================================================
+      SUBROUTINE rnorml(rnd)
+*     Random generator of normal distribution
+      implicit real*8 (A-H,O-Z)
+      PARAMETER (pi=3.141 592 653 589 793 238 462 643d00)
+      PARAMETER (pi2=2.*pi)
+C      
+      u1 = pyr(0)
+      u2 = pyr(0)
+      p1 = pi2*u1
+      p2 = dsqrt(-2.*dlog(u2))
+      rnd=  dcos(p1)*p2
+      RETURN
+      END
diff --git a/EPEMGEN/libepemgen.pkg b/EPEMGEN/libepemgen.pkg
new file mode 100644 (file)
index 0000000..f29b771
--- /dev/null
@@ -0,0 +1,6 @@
+CSRCS:=
+
+FSRCS:= \
+epemgen.f diffcross.f dtrint.f
+
+PACKFFLAGS      := $(filter-out -O%,$(FFLAGS))
diff --git a/TEPEMGEN/TEPEMGENLinkDef.h b/TEPEMGEN/TEPEMGENLinkDef.h
new file mode 100644 (file)
index 0000000..1498d48
--- /dev/null
@@ -0,0 +1,9 @@
+#ifdef __CINT__
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+#pragma link C++ class TEpEmGen+;
+
+#endif
diff --git a/TEPEMGEN/TEcommon.h b/TEPEMGEN/TEcommon.h
new file mode 100644 (file)
index 0000000..b85dbcf
--- /dev/null
@@ -0,0 +1,32 @@
+#ifndef ROOT_TEcommon
+#define ROOT_TEcommon
+//------------------------------------------------------------------------
+// TEcommon is an interface COMMON blocks of the fortran event generator of
+// single e+e- pair production in ultraperipheral PbPb collisions
+// at 5.5 GeV/c
+//%
+// Yuri.Kharlov@cern.ch
+// 9 October 2002
+//------------------------------------------------------------------------
+
+#ifndef __CFORTRAN_LOADED
+#include "cfortran.h"
+#endif
+
+// c++ interface to the f77 program - event generator of
+// e+e- pair production in ultraperipheral ion collisions
+// Author: Yuri Kharlov, 20 September 2002
+
+extern "C" {
+  typedef struct {
+    Double_t Xsect2;
+    Double_t Dsect2;
+    Double_t Xsecttot;
+    Double_t Dsecttot;
+    Int_t    Nevnt;
+  } EeventCommon;
+#define EEVENT COMMON_BLOCK(EEVENT,eevent)
+  COMMON_BLOCK_DEF(EeventCommon,EEVENT);
+}
+
+#endif
diff --git a/TEPEMGEN/TEpEmGen.cxx b/TEPEMGEN/TEpEmGen.cxx
new file mode 100644 (file)
index 0000000..3b28593
--- /dev/null
@@ -0,0 +1,111 @@
+/**************************************************************************
+ * Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ *                                                                        *
+ *                                                                        *
+ * Copyright(c) 1997, 1998, 2002, Adrian Alscher and Kai Hencken          *
+ * See $ALICE_ROOT/EpEmGen/diffcross.f for full Copyright notice          *
+ *                                                                        *
+ *                                                                        *
+ * Copyright(c) 2002 Kai Hencken, Yuri Kharlov, Serguei Sadovsky          *
+ * See $ALICE_ROOT/EpEmGen/epemgen.f for full Copyright notice            *
+ *                                                                        *
+ **************************************************************************/
+//------------------------------------------------------------------------
+// TEpEmGen is an interface class to fortran event generator of
+// single e+e- pair production in ultraperipheral PbPb collisions
+// at 5.5 GeV/c
+//%
+// Yuri.Kharlov@cern.ch
+// 9 October 2002
+//------------------------------------------------------------------------
+
+#include "TEpEmGen.h"
+#include "TClonesArray.h"
+#include "TParticle.h"
+#include "TEcommon.h"
+
+#ifndef WIN32
+# define ee_init  ee_init_
+# define ee_event ee_event_
+#else
+# define ee_init  EE_INIT
+# define ee_event EE_EVENT
+#endif
+
+extern "C" {
+  void ee_init  (Double_t &ymin, Double_t &ymax, Double_t &xmin, Double_t &xmax);
+  void ee_event (Double_t &ymin, Double_t &ymax, Double_t &xmin, Double_t &xmax,
+                Double_t &yE,   Double_t &yP,   Double_t &xE,   Double_t &xP, 
+                Double_t &phi,  Double_t &w);
+}
+
+ClassImp(TEpEmGen)
+
+//------------------------------------------------------------------------------
+TEpEmGen::TEpEmGen() : TGenerator("TEpEmGen","TEpEmGen")
+{
+// TEpEmGen constructor: creates a TClonesArray in which it will store all
+// particles. Note that there may be only one functional TEpEmGen object
+// at a time, so it's not use to create more than one instance of it.
+
+  delete fParticles; // was allocated as TObjArray in TGenerator
+  fParticles = new TClonesArray("TParticle",2);
+
+}
+
+//------------------------------------------------------------------------------
+TEpEmGen::~TEpEmGen()
+{
+  // Destroys the object, deletes and disposes all TParticles currently on list.
+  if (fParticles) {
+    fParticles->Delete();
+    delete fParticles;
+    fParticles = 0;
+  }
+}
+
+//______________________________________________________________________________
+void TEpEmGen::GenerateEvent(Double_t ymin, Double_t ymax, Double_t ptmin, Double_t ptmax,
+                            Double_t &yElectron, Double_t &yPositron,
+                            Double_t &xElectron, Double_t &xPositron,
+                            Double_t &phi12,     Double_t &weight)
+{
+  //produce one event
+  ee_event(ymin,ymax,ptmin,ptmax,
+          yElectron,yPositron,xElectron,xPositron,
+          phi12,weight);
+}
+
+//______________________________________________________________________________
+void TEpEmGen::Initialize(Double_t ymin, Double_t ymax, Double_t ptmin, Double_t ptmax)
+{
+  // Initialize EpEmGen
+  Double_t ptminMeV = ptmin*1000;
+  Double_t ptmaxMeV = ptmax*1000;
+  ee_init(ymin,ymax,ptminMeV,ptmaxMeV);
+}
+
+//______________________________________________________________________________
+Double_t TEpEmGen::GetXsection()
+{
+  // Return cross section accumulated so far
+  return EEVENT.Xsecttot;
+}
+
+//______________________________________________________________________________
+Double_t TEpEmGen::GetDsection()
+{
+  // Return cross section error accumulated so far
+  return EEVENT.Dsecttot;
+}
diff --git a/TEPEMGEN/TEpEmGen.h b/TEPEMGEN/TEpEmGen.h
new file mode 100644 (file)
index 0000000..803d132
--- /dev/null
@@ -0,0 +1,39 @@
+#ifndef ROOT_TEpEmGen
+#define ROOT_TEpEmGen
+/* Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+//------------------------------------------------------------------------
+// TEpEmGen is an interface class to fortran event generator of
+// single e+e- pair production in ultraperipheral PbPb collisions
+// at 5.5 GeV/c
+//%
+// Yuri.Kharlov@cern.ch
+// 9 October 2002
+//------------------------------------------------------------------------
+
+#include "TGenerator.h"
+
+// c++ interface to the f77 program - event generator of
+// e+e- pair production in ultraperipheral ion collisions
+// Author: Yuri Kharlov, 20 September 2002
+
+class TEpEmGen : public TGenerator {
+
+public:
+  TEpEmGen();
+  virtual ~TEpEmGen();
+
+  void Initialize    (Double_t ymin, Double_t ymax, Double_t ptmin, Double_t ptmax);
+  void GenerateEvent (Double_t ymin, Double_t ymax, Double_t ptmin, Double_t ptmax,
+                     Double_t &yElectron, Double_t &yPositron,
+                     Double_t &xElectron, Double_t &xPositron,
+                     Double_t &phi12,     Double_t &weight);
+  Double_t GetXsection();
+  Double_t GetDsection();
+
+  ClassDef(TEpEmGen,1)  //Interface to EpEmGen Event Generator
+};
+
+#endif
diff --git a/TEPEMGEN/libTEPEMGEN.pkg b/TEPEMGEN/libTEPEMGEN.pkg
new file mode 100644 (file)
index 0000000..7dc7089
--- /dev/null
@@ -0,0 +1,7 @@
+SRCS= TEpEmGen.cxx
+
+HDRS= TEpEmGen.h
+
+DHDR:=TEPEMGENLinkDef.h
+
+EXPORT:=TEpEmGen.h