--- /dev/null
+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
--- /dev/null
+ 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
--- /dev/null
+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