program bbbrem

      implicit double precision(a-h,o-z)
      common / experc / roots,rk0
      common / labmom / p1(4),p2(4),q1(4),q2(4),qk(4),d1,d2,t,weight
      common / ranchn / nran
* this is the original file received from R. Kleiss
 
* read input parameters
      read(*,*) roots
      read(*,*) rk0
      read(*,*) nevent
      read(*,*) nran
 
* iitialize bookkeeping
      wmax = 0
      wnil = 0
      wneg = 0
      w0 = 0
      w1 = 0
      w2 = 0
 
* start of event generation loop
      do 9999 k = 1,nevent
 
* generate an event
        call bcube
 
* bookkeeping on the weight
        if(weight.gt.wmax) wmax = weight
        if(weight.eq.0d0) wnil = wnil+1
        if(weight.lt.0d0) wneg = wneg+1
        w0 = w0+1
        w1 = w1+weight
        w2 = w2+weight*weight
 
* end of event generation loop
 9999 continue
 
* analyze weight distribution, and print results
      smc = w1/w0
      ser = dsqrt(w2-w1*w1/w0)/w0
      ranchk = random(0d0)
      write(*,1) ranchk,w0,w1,w2,wmax,wnil,wneg,smc,ser
    1 format(' ******** cross section evaluation ********'/,
     .       ' random number check       ',f15.13/,
     .       ' total # of events         ',f15.0/,
     .       ' sum of weights            ',d15.6/,
     .       ' sum of (weight**2)s       ',d15.6/,
     .       ' max weight occurred       ',d15.6/,
     .       ' # of zero weights         ',f15.0/,
     .       ' # of negative weights     ',f15.0/,
     .       ' computed  cross section   ',d15.6,' millibarn'/,
     .       '                   +/-     ',d15.6,' millibarn'/,
     .       ' ******************************************')
      end
 
      subroutine bcube
      implicit logical(a-z)
      double precision
     . roots,rk0,p1,p2,q1,q2,qk,alpha,rme,tomb,pi,twopi,
     . s,rme2,rme2s,rls,z0,a1,a2,ac,sigapp,eb,pb,rin2pb,random,
     . z,y,q0,temp1,tmin,tmax,sy,w2,rlamx,b,t,rlam,eps,rl,vgam,
     . cgam,sgam,phi,phig,ql,qt,q(4),r0,w,rin2w,rinr0w,eta,phat1,
     . phat2,phat3,phatt,phatl,sfhat,cfhat,sthat,vthat,cthat,sfg,
     . cfg,temp2,veg,qkhat(4),c1,vnumx,vdenx,vnumn,vdenn,c2,rlabl,
     . rhat4,etar,rhat1,rhat2,rhat3,zz,s1,d1,d2,rind1,rind12,
     . rind2,rind22,aa0,aa1,aa2,aa3,aa4,rmex,rmap,weight
      integer init
      save
 
* experimental constants
      common / experc / roots,rk0
* momenta in the lab frame
      common / labmom / p1(4),p2(4),q1(4),q2(4),qk(4),d1,d2,t,weight
 
      data init/0/
* start of initialization
      if(init.eq.0) then
        init = 1
        write(*,901)
  901   format(' ******************************************'/,
     .         ' *** bbbrem : beam-beam bremsstrahlung  ***'/,
     .         ' *** authors r. kleiss and h. burkhardt ***'/,
     .         ' ******************************************')
 
* physical constants
        alpha = 1d0/137.036d0
        rme = 0.51099906d-03
        tomb = 3.8937966d05 /1d6
 
* mathematical constants
        pi = 4*datan(1d0)
        twopi = 2*pi
 
* derived constants
        write(*,902) roots,rk0
  902   format(' total energy              ',f15.6,' gev'/,
     .         ' minimum photon energy     ',f15.6,' * beam energy')
        if(rk0.le.0d0.or.rk0.ge.1d0) then
          write(*,*) 'wrong value for rk0'
          stop
        endif
        s = roots**2
        rme2 = rme**2
        rme2s = rme2/s
        rls = -dlog(rme2s)
        z0 = rk0/(1-rk0)
 
* approximate total cross section
        a1 = dlog((1+z0)/z0)
        a2 = (dlog(1+z0))/z0
        ac = a1/(a1+a2)
        sigapp = 8*alpha**3/rme2*(-dlog(rme2s))*(a1+a2)*tomb
        write(*,903) sigapp
  903   format(' approximate cross section ',d15.6,' millibarn')
 
* the initial-state momenta
        eb = roots*0.5d0
        pb = dsqrt(eb*eb-rme2)
        rin2pb = 0.5d0/pb
        p1(1) = 0
        p1(2) = 0
        p1(3) = -pb
        p1(4) = eb
        q1(1) = 0
        q1(2) = 0
        q1(3) = pb
        q1(4) = eb
 
* end of initialization
      endif
 
* generate z
      if(random(1).lt.ac) then
         temp1=random(2)
         z = 1d0/(temp1*(dexp(a1*random(3))-1))
      else
        z = z0/random(4)
      endif
 
* bounds on t
      y = rme2s*z
      q0 = eb*y
      temp1 = pb*pb-eb*q0
      temp2 = temp1*temp1-rme2-q0*q0
* exit if temp2<0 (very very rare): put weight to 0
* the `else' clause extends to the end of the routine
      if(temp2.lt.0d0) then
        write(*,904) temp2
  904   format(' y too large: delta_t^2 = ',d15.6)
        weight = 0d0
      else
        tmin = -2*(temp1+dsqrt(temp2))
        tmax = rme2*s*y*y/tmin
 
* generate t
        sy = s*y
        w2 = sy+rme2
        temp1 = sy+tmax
        rlamx = dsqrt(temp1*temp1-4*w2*tmax)
        if(temp1.le.0d0) then
          temp1 = rlamx-temp1
        else
          temp1 = -4*w2*tmax/(rlamx+temp1)
        endif
    1   continue
          b = dexp(random(5)*dlog(1+2*sy/temp1))
          t = -b*z*z*rme2/((b-1)*(b*z+b-1))
          if(t.lt.tmin) then
            write(*,905) t,tmin
  905       format(' t = ',d15.6,'   < t_min =',d15.6)
            goto 1
          endif
        continue
 
* generate cgam
        rlam = dsqrt((sy-t)*(sy-t)-4*rme2*t)
        eps = 4*rme2*w2/(rlam*(rlam+w2+rme2-t))
        rl = dlog((2+eps)/eps)
        vgam = eps*(dexp(random(6)*rl)-1)
        cgam = 1-vgam
        sgam = dsqrt(vgam*(2-vgam))
 
* generate azimuthal angles
        phi = twopi*random(7)
        phig = twopi*random(8)
 
* construct momentum transfer q(mu)
        ql = (2*eb*q0-t)*rin2pb
        qt = dsqrt((tmax-t)*(t-tmin))*rin2pb
        q(1) = qt*dsin(phi)
        q(2) = qt*dcos(phi)
        q(3) = ql
        q(4) = q0
 
* construct momentum of outgoing positron in lab frame
        q2(1) = q1(1)-q(1)
        q2(2) = q1(2)-q(2)
        q2(3) = q1(3)-q(3)
        q2(4) = q1(4)-q(4)
 
* find euler angles of p1(mu) in cm frame
        r0 = eb+q0
        w = dsqrt(w2)
        rin2w = 0.5d0/w
        rinr0w = 1d0/(r0+w)
        eta = -(sy+2*w*q0+t)*rin2w*rinr0w
        phat1 = -q(1)*(1+eta)
        phat2 = -q(2)*(1+eta)
        phat3 = pb*eta-ql*(1+eta)
        phatl = rlam*rin2w
        phatt = dsqrt(phat1*phat1+phat2*phat2)
        sfhat = phat1/phatt
        cfhat = phat2/phatt
        sthat = phatt/phatl
        if(phat3.gt.0d0) then
          vthat = sthat*sthat/(1-dsqrt(1-sthat*sthat))
        else
          vthat = sthat*sthat/(1+dsqrt(1-sthat*sthat))
        endif
        cthat = vthat-1
 
* rotate using these euler angles to get the qk direction in the cm
        sfg = dsin(phig)
        cfg = dcos(phig)
        temp1 = sgam*sfg
        temp2 = cthat*sgam*cfg+sthat*cgam
        veg = vthat+vgam-vthat*vgam-sthat*sgam*cfg
        qkhat(4) = sy*rin2w
        qkhat(1) = qkhat(4)*( cfhat*temp1+sfhat*temp2)
        qkhat(2) = qkhat(4)*(-sfhat*temp1+cfhat*temp2)
        qkhat(3) = qkhat(4)*(veg-1)
 
* boost the photon momentum to the lab frame
        temp1 = pb*qkhat(3)
        if(temp1.gt.0d0) then
          temp2 = (rme2*qkhat(4)*qkhat(4)
     .            +pb*pb*(qkhat(1)*qkhat(1)+qkhat(2)*qkhat(2)))
     .            /(eb*qkhat(4)+temp1)
        else
          temp2 = eb*qkhat(4)-temp1
        endif
        qk(4) = (temp2+qkhat(4)*q(4)+qkhat(1)*q(1)
     .                +qkhat(2)*q(2)+qkhat(3)*q(3))/w
        temp1 = (qk(4)+qkhat(4))*rinr0w
        qk(1) = qkhat(1)+temp1*q(1)
        qk(2) = qkhat(2)+temp1*q(2)
        qk(3) = qkhat(3)+temp1*(-pb+ql)
 
* construct p2 by momentum conservation
        p2(1) = -q2(1)-qk(1)
        p2(2) = -q2(2)-qk(2)
        p2(3) = -q2(3)-qk(3)
        p2(4) = -q2(4)-qk(4)+roots
 
* impose cut on the photon energy: qk(4)>eb*rk0
        if(qk(4).lt.eb*rk0) then
          weight = 0d0
        else
 
* the event is now accepted: compute matrix element and weight
* compute fudge factor c1
          c1 = dlog(1+z)/dlog((2+eps)/eps)
 
* compute fudge factor c2
          temp1 = sy-tmax
          vnumx = dsqrt(temp1*temp1-4*rme2*tmax)+temp1
          temp1 = sy+tmax
          if(temp1.lt.0d0) then
            vdenx = dsqrt(temp1*temp1-4*w2*tmax)-temp1
          else
            vdenx = -4*w2*tmax/(dsqrt(temp1*temp1-4*w2*tmax)+temp1)
          endif
          temp1 = sy-tmin
          vnumn = dsqrt(temp1*temp1-4*rme2*tmin)+temp1
          temp1 = sy+tmin
          if(temp1.lt.0d0) then
            vdenn = dsqrt(temp1*temp1-4*w2*tmin)-temp1
          else
            vdenn = -4*w2*tmin/(dsqrt(temp1*temp1-4*w2*tmin)+temp1)
          endif
          c2 = 2*rls/dlog((vnumx*vdenn)/(vdenx*vnumn))
 
* compute vector (small) r in cm frame, and (big) z
          rlabl = (t-2*rme2*y)*rin2pb
          rhat4 = -(2*rme2*y+(1-y)*t)*rin2w
          etar = rhat4*rinr0w
          rhat1 = -q(1)*(1+etar)
          rhat2 = -q(2)*(1+etar)
          rhat3 = rlabl+(pb-ql)*etar
          zz = s*(rhat4*qkhat(4)-rhat1*qkhat(1)
     .         -rhat2*qkhat(2)-rhat3*qkhat(3))
 
* the other invariants
          s1 = 4*eb*(eb-qk(4))
          d1 = sy*rlam*(eps+vgam)*rin2w*rin2w
          d2 = 0.5d0*sy
 
* the exact matrix element
          rind1 = 1d0/d1
          rind12 = rind1*rind1
          rind2 = 1d0/d2
          rind22 = rind2*rind2
c* kleiss-burkhardt cross section multiplied by t**2
          temp1 = s+t-2*d2
          temp2 = s1+t+2*d1
          aa0 = (s*s+s1*s1+temp1*temp1+temp2*temp2)*rind1*rind2*(-t)
          aa1 = -4*rme2*zz*zz*rind12*rind22
          aa2 = -8*rme2*(d1*d1+d2*d2)*rind1*rind2
          aa3 = 16*rme2*rme2*(d1-d2)*zz*rind12*rind22
          aa4 = -16*rme2*rme2*rme2*(d1-d2)*(d1-d2)*rind12*rind22
          rmex = aa0+aa1+aa2+aa3+aa4
 
* the approximate matrix element without c1,2, multiplied by t**2
          rmap = 4*s*s*rind1*rind2*(-t)*c1*c2
 
* the weight
          weight = rmex/rmap*sigapp
 
* the weight is now defined for both accepted and rejected events
        endif
      endif
 1001 format(3(a10,d15.6))
      end
 
      function random(r)
      implicit double precision(a-h,o-z)
      common / ranchn / nran
      save
      data init/0/
      if(init.eq.0) then
        init = 1
        if(nran.eq.1) then
          write(*,*) ' random: very crude random numbers chosen'
        elseif(nran.eq.2) then
          write(*,*) ' random: quasi-random numbers chosen'
        else
          write(*,*) ' random: error: nran = ',nran,' not allowed'
          stop
        endif
      endif
      if(nran.eq.1) then
        random = rando1(r)
      elseif(nran.eq.2) then
        random = rando2(r)
      endif
      end
 
      function rando1(r)
      implicit double precision(a-h,o-z)
      save
      data init/0/
      if(init.eq.0) then
        init = 1
        m = 2**20
        rm = 1d0/(1d0*m)
        ia = 2**9+5
        ic = 1
        iseed = 100461
      endif
    1 iseed = mod(ia*iseed+ic,m)
      if(iseed.eq.0) goto 1
      rando1 = iseed*rm
      end
 
      function rando2(r)
      implicit double precision(a-h,o-z)
      save
      parameter(n = 3)
      dimension nprime(10),s(10)
      data init/0/
      data nprime/2,3,5,7,11,13,17,19,23,29/
      if(init.eq.0) then
        init = 1
        do 1 k = 1,n
          s(k) = dsqrt(nprime(k)*1d0)
    1   continue
      endif
      do 3 k = 1,n
        ss = 0
        do 2 j=1,n
          ss = ss+s(j)
    2   continue
        s(k) = dmod(ss,1d0)
    3 continue
      rando2 = s(1)
      end