*****************************************************************************
*****************************************************************************
The routines MODL2A and MODL2B below has to be replaced in BHLUMI 4.02
in order to correct for the LL bug which I have reprted during today's
meeting of the Bhabha Working Group in CERN.

Note that the "official weight" 
         WTMOD = WTCRU1*WTCRU2*WTSET(142)
is represents now corrected matrix element and 
         WTOLD = WTCRU1*WTCRU2*WTSET(242)
represents old weight with the bug.
This is done in order to alow user to calculate quickly
the correction to any distribution using the relative weight:
         WTCORR = WTCRU1*WTCRU2*(WTSET(142)-WTSET(242))
The same is done with unexponentiated matrix element
where WTSET(272) is included for the same purpose.

Best regards, Staszek Jadach, CERN 3-rd October 1995.


*****************************************************************************
*****************************************************************************

      SUBROUTINE MODL2A(MODE)
!     ***********************
!======================================================================!
!                                                                      !
!    OLD type Matrix element type (A)                                  !
!    similar as in BHLUMI version. 2.01  (Identical in first order)    !
!    but now also in the second order     (January. 94).               !
!    Exponentiated and not exponentiated  (January. 95)                !
!                                                                      !
!======================================================================!
!     ************************
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER( PI = 3.1415926535897932D0, ALFINV=137.03604D0)
      PARAMETER( ALFA= 1D0/ALFINV, ALF1 = 1D0/PI/ALFINV)
      COMMON / BHPAR1 / DEL,EPSCM,THMIN,XMIVIS
      COMMON / BHPAR2 / CMSENE,AMEL
      COMMON / BHPAR3 / KEYRAD,KEYOPT
      COMMON / BHPARZ / AMAZ,GAMMZ,SINW2,GV,GA
      COMMON / TRANSR / TRAN,TRMIN,TRMAX
      COMMON / MOMS1  / TRANP,P1(4),P2(4),PHOT1(100,4),PHSU1(4),NPHOT1
      COMMON / MOMS2  / TRANQ,Q1(4),Q2(4),PHOT2(100,4),PHSU2(4),NPHOT2
      COMMON / MOMZ1  / AL1(50),BE1(50),WTM1(50),MK1(50)
      COMMON / MOMZ2  / AL2(50),BE2(50),WTM2(50),MK2(50)
      COMMON / PSIPHI / TH1,EXT1,EXB1,PSI1,EXW1,
     $                  TH2,EXT2,EXB2,PSI2,EXW2,EXE2,FIF,PHI
      COMMON / WGTALL / WTMOD,WTCRU1,WTCRU2,WTSET(300)
      SAVE   / BHPAR2 /, / BHPAR3 /, / TRANSR /
      SAVE   / MOMS1  /, / MOMS2  /, / MOMZ1  /, / MOMZ2  /
      SAVE   / WGTALL /, / BHPARZ /, / PSIPHI /
      DIMENSION ph1(4),ph2(4)
      DIMENSION sfu(100),sfmu(100),bu10(100),bu11(100),bu11g(100)
      DIMENSION sfl(100),sfml(100),bl10(100),bl11(100),bl11g(100)
! ------------------ Inline functions ------------------
! Elements of single bremss. distribution
      CHI(x)= (1+(1-x)**2)/2
      XDELS(a,b,delt) =
     $ delt*(a**2+b**2)*(a**2+b**2)/((1-a)**2+(1-b)**2)/(a*b)
!==================================================================!
      DATA icont /0/
      icont=icont+1

      s = cmsene**2
      zeta = abs(tran)/s
!*****************************************************
!  Crude MC distribution (S-tilde factors omitted)
      CRUDE  =  s**2/tran**2
      facpq  =  s**2/(tranp*tranq)
      delta  =  amel**2/tran
      deltp  =  amel**2/tranp
      deltq  =  amel**2/tranq
      gam    =  2*alf1 *( dlog(1/delta)-1 )
      gamp   =  2*alf1 *( dlog(1/deltp)-1 )
      gamq   =  2*alf1 *( dlog(1/deltq)-1 )
      dis0   =  facpq *dix1(p1,p2,q1,q2)
!==================================================================!
      bt10u =0
      bt10l =0
      bt11u =0
      bt11l =0
      bt20u =0
      bt20l =0
      bt2ul =0
! ------------------
! beta0
! ------------------
      bt00   = dis0
      bt01   = bt00*(1 + gam )
      bt02   = bt00*(1 + gam +0.5*gam**2)
! ------------------------------------
! Contributions from beta1 UPPER line
! ------------------------------------
      DO 200 i=1,nphot1
! numerically safe variables
      a   = al1(i)
      b   = be1(i)
      v   = a+b-a*b
      y   = a +b*deltp
      z   = b +a*deltp
! soft factor
      sfc  = a*b/(y*z)**2
      sfu(i) = sfc
! one photon bremss. distribution
      wm = a*b/(y*z) +deltp*(y**2+z**2)/((1-y)**2+(1-z)**2)*(y/z+z/y)
      sfmu(i) = wm/(y*z)
      dis10   = sfmu(i)*dis0
      dis11g  = dis10 *(1+gam + gam/2*log((1-b)/(1-a)))
      dis11   = dis10 *(1+gam + gam/4*log((1-b)/(1-a)))
c[[[[[[
c>!#### dis11   = dis10 *(1+gam + gam/2*log((1-b)**2/(1-v)))
c>      dis11   = dis10 *(1+gam + gam/4*log((1-b)**2/(1-v)))
c]]]]]]
! beta1 O(alf1),O(alf2)
      bu10(i) =  dis10   -bt00*sfc
      bu11(i) =  dis11   -bt01*sfc
      bu11g(i)=  dis11g  -bt01*sfc   ! with LL bug
      bt10u   =  bt10u   +bu10(i)/sfc
      bt11u   =  bt11u   +bu11(i)/sfc
  200 CONTINUE
! -------------------------------------
!  Contributions from beta1 LOWER line
! -------------------------------------
      DO 250 i=1,nphot2
      a   = al2(i)
      b   = be2(i)
      v   = a+b-a*b
      y   = a+b*deltq
      z   = b+a*deltq
! soft factor
      sfc  = a*b/(y*z)**2
      sfl(i) = sfc
! one photon bremss. distribution
      wm = a*b/(y*z) +deltq*(y**2+z**2)/((1-y)**2+(1-z)**2)*(y/z+z/y)
      sfml(i) = wm/(y*z)
      dis10   = sfml(i)*dis0
      dis11g  = dis10 *(1+gam + gam/2*log((1-b)/(1-a)))
      dis11   = dis10 *(1+gam + gam/4*log((1-b)/(1-a)))
c[[[[[[
c>!#### dis11   = dis10 *(1+gam + gam/2*log((1-b)**2/(1-v)))
c>      dis11   = dis10 *(1+gam + gam/4*log((1-b)**2/(1-v)))
c]]]]]]
! beta1 O(alf1),O(alf2)
      bl10(i) =  dis10   -bt00*sfc
      bl11(i) =  dis11   -bt01*sfc
      bl11g(i)=  dis11g  -bt01*sfc  ! with LL bug
      bt10l   =  bt10l   +bl10(i)/sfc
      bt11l   =  bt11l   +bl11(i)/sfc
  250 CONTINUE

!------------------
! beta2 upper line
!------------------
      do 253 i=1,nphot1
      do 253 j=i+1,nphot1
! Define photon four-vectors for further local use
      do 252 k=1,4
      ph1(k) = phot1(i,k)
 252  ph2(k) = phot1(j,k)
! Basic variables first photon
      a1  = al1(i)
      b1  = be1(i)
      y1  = a1 +b1*deltp
      z1  = b1 +a1*deltp
      v1  = a1+b1-a1*b1
! Basic variables second photon
      a2  = al1(j)
      b2  = be1(j)
      y2  = a2 +b2*deltp
      z2  = b2 +a2*deltp
      v2  = a2+b2-a2*b2
! Soft factors, starred variables
      sf1  =  sfu(i)
      sf2  =  sfu(j)
      a1st = a1/(1-v2)
      b1st = b1/(1-v2)
      v1st = v1/(1-v2)
      a2st = a2/(1-v1)
      b2st = b2/(1-v1)
      v2st = v2/(1-v1)
      a1st = min(a1st,1-zeta)
      b1st = min(b1st,1-zeta)
      v1st = min(v1st,1-zeta)
      a2st = min(a2st,1-zeta)
      b2st = min(b2st,1-zeta)
      v2st = min(v2st,1-zeta)
! H-functions exact and Leading Log
      gs1    = sf1*chi(v1)
      gs2    = sf2*chi(v2)
      gs1st  = sf1*chi(v1st)
      gs2st  = sf2*chi(v2st)
! The actual two photon distributions
      IF ( (a1-b1)*(a2-b2).le.0d0 ) THEN
!       OPPOSITE directions two photons
        IF ((a1-b1).gt.0d0) THEN
!          FIRST photon in INITIAL state, second in the final
           hs1f   = DIX2(2,ph2,p1,p2,q1,q2) *sfmu(i)
           hs2i   = DIX2(1,ph1,p1,p2,q1,q2) *sfmu(j)
           IF(v1.gt.v2) THEN
               d2uu  = facpq  *hs1f  *gs2
           ELSE
               d2uu  = facpq  *gs1st *hs2i
           ENDIF
        ELSE
!          SECOND photon in INITIAL state, first in the final
           hs1i   = DIX2(1,ph2,p1,p2,q1,q2) *sfmu(i)
           hs2f   = DIX2(2,ph1,p1,p2,q1,q2) *sfmu(j)
           IF(v1.gt.v2) THEN
               d2uu  = facpq  *hs1i *gs2st
           ELSE
               d2uu  = facpq  *gs1  *hs2f
           ENDIF
        ENDIF
      ELSE
!       SAME directions two photons
!         Both in INITIAL  state or Both in FINAL  state
          IF(v1.gt.v2) THEN
            hs1i  = DIX2(1,ph2,p1,p2,q1,q2) *sfmu(i)
            hs1f  = DIX2(2,ph2,p1,p2,q1,q2) *sfmu(i)
            d2uu  = facpq/2 *( hs1f*gs2st + hs1i*gs2)
          ELSE
            hs2i  = DIX2(1,ph1,p1,p2,q1,q2) *sfmu(j)
            hs2f  = DIX2(2,ph1,p1,p2,q1,q2) *sfmu(j)
            d2uu  = facpq/2 *( hs2f*gs1st + hs2i*gs1)
          ENDIF
      ENDIF
      bu20  = d2uu-bu10(i)*sfu(j)-bu10(j)*sfu(i)-bt00*sfu(j)*sfu(i)
      bt20u = bt20u + bu20/sfu(i)/sfu(j)
 253  CONTINUE
! ------------------
! beta2 lower line
! ------------------
      DO 258 i=1,nphot2
      DO 258 j=i+1,nphot2
! Define photon four-vectors for further local use
      do 257 k=1,4
      ph1(k) = phot2(i,k)
 257  ph2(k) = phot2(j,k)
! Basic variables first photon
      a1  = al2(i)
      b1  = be2(i)
      y1  = a1 +b1*deltq
      z1  = b1 +a1*deltq
      v1  = a1+b1-a1*b1
! Basic variables second photon
      a2  = al2(j)
      b2  = be2(j)
      y2  = a2 +b2*deltq
      z2  = b2 +a2*deltq
      v2  = a2+b2-a2*b2
! soft factors, starred variables
      sf1  =  sfl(i)
      sf2  =  sfl(j)
      a1st = a1/(1-v2)
      b1st = b1/(1-v2)
      v1st = v1/(1-v2)
      a2st = a2/(1-v1)
      b2st = b2/(1-v1)
      v2st = v2/(1-v1)
      a1st = min(a1st,1-zeta)
      b1st = min(b1st,1-zeta)
      v1st = min(v1st,1-zeta)
      a2st = min(a2st,1-zeta)
      b2st = min(b2st,1-zeta)
      v2st = min(v2st,1-zeta)
! H-functions exact and Leading Log
      gs1    = sf1*chi(v1)
      gs2    = sf2*chi(v2)
      gs1st  = sf1*chi(v1st)
      gs2st  = sf2*chi(v2st)
! The actual two photon distributions
      IF ( (a1-b1)*(a2-b2).le.0d0 ) THEN
!       OPPOSITE directions two photons
        IF ((a1-b1).gt.0d0) THEN
!          FIRST photon in INITIAL state, second in the final
           hs1f   = DIX2(4,ph2,p1,p2,q1,q2) *sfml(i)
           hs2i   = DIX2(3,ph1,p1,p2,q1,q2) *sfml(j)
           IF(v1.gt.v2) THEN
               d2ll  = facpq  *hs1f  *gs2
           ELSE
               d2ll  = facpq  *gs1st*hs2i
           ENDIF
        ELSE
!          SECOND photon in INITIAL state, first in the final
           hs1i   = DIX2(3,ph2,p1,p2,q1,q2) *sfml(i)
           hs2f   = DIX2(4,ph1,p1,p2,q1,q2) *sfml(j)
           IF(v1.gt.v2) THEN
               d2ll  = facpq  *hs1i *gs2st
           ELSE
               d2ll  = facpq  *gs1 *hs2f
           ENDIF
        ENDIF
      ELSE
!       SAME directions two photons
!         Both in INITIAL  state or Both in FINAL  state
          IF(v1.gt.v2) THEN
            hs1i  = DIX2(3,ph2,p1,p2,q1,q2) *sfml(i)
            hs1f  = DIX2(4,ph2,p1,p2,q1,q2) *sfml(i)
            d2ll  = facpq/2 *( hs1f*gs2st + hs1i*gs2)
          ELSE
            hs2i  = DIX2(3,ph1,p1,p2,q1,q2) *sfml(j)
            hs2f  = DIX2(4,ph1,p1,p2,q1,q2) *sfml(j)
            d2ll  = facpq/2 *( hs2f*gs1st + hs2i*gs1)
          ENDIF
      ENDIF
      bl20  = d2ll-bl10(i)*sfl(j)-bl10(j)*sfl(i)-bt00*sfl(j)*sfl(i)
      bt20l = bt20l + bl20/sfl(i)/sfl(j)
 258  CONTINUE
!=================================================
! [11] One upper and one lower line photons
! Note that bt2ul=0 in the LL.
! WARNING!!! This part is not tested
! It should be completely negligible!!!!
      DO i=1,nphot1
      DO j=1,nphot2
! Basic variables first photon
        a1  = al1(i)
        b1  = be1(i)
        v1  = a1 +b1-a1*b1
! Basic variables second photon
        a2  = al2(j)
        b2  = be2(j)
        v2  = a2+b2-a2*b2
! soft factors
        sf1  =  sfu(i)
        sf2  =  sfl(j)
        IF(v1.gt.v2) THEN
           d2ul= dis0 *sfmu(i) *sfl(j)
         ELSE
           d2ul= dis0 *sfu(i)  *sfml(j)
        ENDIF
        b2ul  = d2ul-bu10(i)*sfl(j)-sfu(i)*bl10(j)-bt00*sfu(i)*sfl(j)
        bt2ul = bt2ul + b2ul/sfu(i)/sfl(j)
      ENDDO
      ENDDO
     

!        **************************************
!        **     Definitions of MC weights    **
!        **************************************
!          ---------------------------------
!          !    UPPER line ONLY,    (A)    !
!          ---------------------------------
! Case of upper line only has separate entries in WTSET because it is 
! analysed in the separate run with separate plotting programs 
! with different entries in vvrho, differerent histogram names etc.
!
! All beta's:  TOTAL O(alf0),O(alf1),O(alf2)
      WTSET( 30) =    bt00                /crude
      WTSET( 31) =   (bt01+bt10u)         /crude
      WTSET( 32) =   (bt02+bt11u +bt20u)  /crude
! Special bet0-type distr. (without chi0) for tests only 
      WTSET( 34) =   ( s**2/(tranp*tranq) )/crude
! Individual beta's in various orders.
! O(alf1)
      WTSET( 35) =   bt01  /crude                  ! beta0
      WTSET( 36) =   bt10u /crude                  ! beta1
! O(alf2)
      WTSET( 37) =   bt02  /crude                  ! beta0
      WTSET( 38) =   bt11u /crude                  ! beta1
      WTSET( 39) =   bt20u /crude                  ! beta2
!          ---------------------------------
!          !  UPPER + LOWER line,    (A)   !
!          ---------------------------------
! All beta's:  O(alf0),O(alf1),O(alf2)
      WTSET( 40) =    bt00  /crude
      WTSET( 41) =   (bt01+bt10u+bt10l) /crude
      WTSET( 42) =   (bt02+bt11u+bt11l+bt20u+bt20l+bt2ul) /crude
! Individual beta's in various orders.
! O(alf1)
      WTSET( 43) =    bt01   /crude                ! beta0
      WTSET( 44) =   (bt10u+bt10l) /crude          ! beta1
      WTSET( 45) =    bt10u  /crude                ! beta1 components
      WTSET( 46) =    bt10l  /crude
! O(alf2)
      WTSET( 47) =    bt02   /crude                ! beta0
      WTSET( 48) =   (bt11u+bt11l)/crude           ! beta1
      WTSET( 49) =   (bt20u+bt20l+bt2ul)/crude     ! beta2
      WTSET( 50) =    bt11u  /crude                ! beta1 components
      WTSET( 51) =    bt11l  /crude
      WTSET( 52) =    bt2ul  /crude
      WTSET( 53) =    bt20u  /crude                ! beta2 components
      WTSET( 54) =    bt20l  /crude

!==================================================================
!==================================================================
!                  Non-exponentiated version                      !
!==================================================================
!==================================================================
      pdel  = del*bcud(p1,p2,phsu1)
      qdel  = del*bcud(q1,q2,phsu2)
      fyfsu = exp(-gamp*dlog(1/pdel)  +gamp/4 -alf1/2)
      fyfsl = exp(-gamq*dlog(1/qdel)  +gamq/4 -alf1/2)
!
      bb0u    =0
      bb1u    =0
      bb2u    =0
      dd0u    =0
      dd1u    =0
      dd2u    =0
      IF(nphot1.eq.0) THEN
        bb0u   = dis0
        bb1u   = bb0u*(1 -gam*dlog(1/pdel))
        bb2u   = bb0u*(1 -gam*dlog(1/pdel) +1./2*(gam*dlog(1/pdel))**2)
        virt1  = 3./4*gam -1./2*alf1
        bb0u   = dis0
        dd1u   = dis0 *(1 -gam*dlog(1/pdel) +virt1)
      ELSEIF(nphot1.eq.1) THEN
        a   = al1(1)
        b   = be1(1)
        y   = a +b*deltp
        z   = b +a*deltp
        sfc = a*b/(y*z)**2
        bb1u = dis0
        wm  = a*b/(y*z) +deltp*(y**2+z**2)/((1-y)**2+(1-z)**2)*(y/z+z/y)
        sfm = wm/(y*z)
        dist1  = dis0 *sfm
        dd1u    = dist1/sfc
      ENDIF
!
      bb1  =0
      bb2  =0
!        **************************************
!        ****   Definitions of MC weights  ****
!        **************************************
! Case of upper-line-only has separate entries in WTSET because it is 
! analysed in the separate run with separate plotting programs 
! with different entries in vvrho, differerent histogram names etc.
!
! YFS formfactor fyfsu has to be removed in the unexponentiated case!
! Baseline
      WTSET(56) =    bb1u               /fyfsu/crude ! 1-photon
      WTSET(57) =    bb2u               /fyfsu/crude ! 2-photons
      WTSET(58) =    bb1                /fyfsu/crude ! 1-photon
      WTSET(59) =    bb2                /fyfsu/crude ! 2-photons
!          ---------------------------------
!          /////    UPPER line ONLY    /////
!          ---------------------------------
! TOTAL O(alf0),O(alf1),O(alf2)
      WTSET(60) =    dd0u                /fyfsu/crude
      WTSET(61) =    dd1u                /fyfsu/crude
      WTSET(62) =    dd2u                /fyfsu/crude
! Baseline
      WTSET(63) =    0
      WTSET(64) =    0
!************************************************************************

!    ////////////////////////////////////////////
!    /////  UNEXP   UPPER + LOWER line      /////
!    ////////////////////////////////////////////

! Entire 0,1,2-photon distributions
      dis0   =0
      dis1   =0
      dis2   =0
!
      virt1y = 6./4*gam -alf1
      virt2y = 0.5*(6./4*gam)**2
      h0     =  dix1(p1,p2,q1,q2)
!=================================================
      IF(nphot1.eq.0.and.nphot2.eq.0) THEN
! [00] No photons
!------
! O(alf0,1,2) Entire distribution 
        dis0= facpq*h0
        dis1= facpq*h0*(1+2*gam*dlog(del))
     $       +facpq*h0*virt1y
        dis2= facpq*h0*(1+2*gam*dlog(del)+0.5*(2*gam*dlog(del))**2)
     $       +facpq*h0*(1+2*gam*dlog(del))*virt1y
     $       +facpq*h0* virt2y
!=================================================
      ELSEIF(nphot1.eq.1.and.nphot2.eq.0) THEN
! [10] One upper line photon 
!  O(alf1,2) Entire distribution
!  Note that 1/2 in the log(1-v) term = {1/2 from two tree diagrams} 
        dis1   = facpq*h0*(1+xdels(al1(1),be1(1),deltp))
        dis2g  = dis1*(1 +(gamp+gamq)*dlog(del) +virt1y
     $        -1/2d0*gam*log(1-al1(1)) -1/2d0*gam*log(1-be1(1)) )
        dis2   = dis1*(1 +(gamp+gamq)*dlog(del) +virt1y
     $        -1/4d0*gam*log(1-al1(1)) -3/4d0*gam*log(1-be1(1)) )
c[[[[[
c>        dis2   = dis1*(1 +(gamp+gamq)*dlog(del) 
c>!####$       +virt1y -gam/2*log(1-v) )
c>     $       +virt1y -gam/4*log(1-v) -gam/2*log(1-be1(1)) )
c]]]]]
!=================================================
      ELSEIF(nphot1.eq.0.and.nphot2.eq.1) THEN
! [01] One lower line photon 
!------
! O(alf1,2) Entire distribution
! Note that 1/2 in the log(1-v) term = {1/2 from two tree diagrams} 
        dis1= facpq*h0*(1+xdels(al2(1),be2(1),deltq))
        dis2g= dis1*(1 +(gamp+gamq)*dlog(del) +virt1y
     $        -1/2d0*gam*log(1-al2(1)) -1/2d0*gam*log(1-be2(1)) )
        dis2 = dis1*(1 +(gamp+gamq)*dlog(del) +virt1y
     $        -1/4d0*gam*log(1-al2(1)) -3/4d0*gam*log(1-be2(1)) )
c[[[[[
c>        dis2= dis1*(1 +(gamp+gamq)*dlog(del) 
c>!####$       +virt1y -gam/2*log(1-v) )
c>     $       +virt1y -gam/4*log(1-v) -gam/2*log(1-be2(1)) )
c]]]]]
!=================================================
      ELSEIF(nphot1.eq.2.and.nphot2.eq.0) THEN
! [20] Two upper line photons 
! Basic distribution
        a1  = al1(1)
        b1  = be1(1)
        a2  = al1(2)
        b2  = be1(2)
        v1  = a1+b1-a1*b1
        v2  = a2+b2-a2*b2
! Define photon four-vectors for further local use
        DO k=1,4
          ph1(k) = phot1(1,k)
          ph2(k) = phot1(2,k)
        ENDDO
        a1st = min(a1/(1-v2),1-zeta)
        b1st = min(b1/(1-v2),1-zeta)
        a2st = min(a2/(1-v1),1-zeta)
        b2st = min(b2/(1-v1),1-zeta)
        v1st = min(v1/(1-v2),1-zeta)
        v2st = min(v2/(1-v1),1-zeta)
!       H-functions exact
        hs1i = DIX2(1,ph2,p1,p2,q1,q2)*(1+xdels(a1,b1,deltp))
        hs1f = DIX2(2,ph2,p1,p2,q1,q2)*(1+xdels(a1,b1,deltp))
        hs2i = DIX2(1,ph1,p1,p2,q1,q2)*(1+xdels(a2,b2,deltp))
        hs2f = DIX2(2,ph1,p1,p2,q1,q2)*(1+xdels(a2,b2,deltp))
!       The actual two photon distributions
        IF ( (a1-b1)*(a2-b2).le.0d0 ) THEN
!       OPPOSITE directions two photons
          IF (a1.gt.b1) THEN
!          FIRST photon in INITIAL state, second in the final
             IF(v1.gt.v2) dis2  = facpq* hs1f      *chi(v2)
             IF(v2.ge.v1) dis2  = facpq* chi(v1st) *hs2i
          ELSE
!          SECOND photon in INITIAL state, first in the final
             IF(v1.gt.v2) dis2  = facpq* hs1i    *chi(v2st)
             IF(v2.ge.v1) dis2  = facpq* chi(v1) *hs2f
          ENDIF
        ELSE
!       SAME directions two photons
          IF(v1.gt.v2) dis2 = facpq/2 *(hs1f*chi(v2st) + hs1i*chi(v2))
          IF(v2.ge.v1) dis2 = facpq/2 *(hs2f*chi(v1st) + hs2i*chi(v1))
        ENDIF
!=================================================
      ELSEIF(nphot1.eq.0.and.nphot2.eq.2) THEN
! [02] Two lower line photons 
! Basic distribution
        a1  = al2(1)
        b1  = be2(1)
        a2  = al2(2)
        b2  = be2(2)
        v1  = a1+b1-a1*b1
        v2  = a2+b2-a2*b2
        a1st = min(a1/(1-v2),1-zeta)
        b1st = min(b1/(1-v2),1-zeta)
        a2st = min(a2/(1-v1),1-zeta)
        b2st = min(b2/(1-v1),1-zeta)
        v1st = min(v1/(1-v2),1-zeta)
        v2st = min(v2/(1-v1),1-zeta)
! Define photon four-vectors for further local use
        DO k=1,4
          ph1(k) = phot2(1,k)
          ph2(k) = phot2(2,k)
        ENDDO
!       H-functions exact
        hs1i = DIX2(3,ph2,p1,p2,q1,q2)*(1+xdels(a1,b1,deltq))
        hs1f = DIX2(4,ph2,p1,p2,q1,q2)*(1+xdels(a1,b1,deltq))
        hs2i = DIX2(3,ph1,p1,p2,q1,q2)*(1+xdels(a2,b2,deltq))
        hs2f = DIX2(4,ph1,p1,p2,q1,q2)*(1+xdels(a2,b2,deltq))
!       The actual two photon distributions
        IF ( (a1-b1)*(a2-b2).le.0d0 ) THEN
!       OPPOSITE directions two photons
          IF (a1.gt.b1) THEN
!          FIRST photon in INITIAL state, second in the final
             IF(v1.gt.v2) dis2  = facpq* hs1f      *chi(v2)
             IF(v2.ge.v1) dis2  = facpq* chi(v1st) *hs2i
          ELSE
!          SECOND photon in INITIAL state, first in the final
             IF(v1.gt.v2) dis2  = facpq* hs1i    *chi(v2st)
             IF(v2.ge.v1) dis2  = facpq* chi(v1) *hs2f
          ENDIF
        ELSE
!       SAME directions two photons
          IF(v1.gt.v2) dis2 = facpq/2 *(hs1f*chi(v2st) + hs1i*chi(v2))
          IF(v2.ge.v1) dis2 = facpq/2 *(hs2f*chi(v1st) + hs2i*chi(v1))
        ENDIF
!=================================================
      ELSEIF(nphot1.eq.1.and.nphot2.eq.1) THEN
! [11] One upper and lower line photons 
        a1  = al1(1)
        b1  = be1(1)
        a2  = al2(1)
        b2  = be2(1)
        v1  = a1+b1-a1*b1
        v2  = a2+b2-a2*b2
        IF(v1.gt.v2) THEN
          dis2 = facpq *dix1(p1,p2,q1,q2)*(1+xdels(a1,b1,deltp))
        ELSE
          dis2 = facpq *dix1(p1,p2,q1,q2)*(1+xdels(a2,b2,deltq))
        ENDIF
      ENDIF
!     =====

!        **************************************
!        ****   Definitions of MC weights  ****
!        ****    Non-exponentiated case    ****
!        **************************************
!          ---------------------------------
!          /////  UPPER + LOWER line   /////
!          ---------------------------------
 
      fyfs = fyfsu*fyfsl
! TOTAL O(alf0),O(alf1),O(alf2)
      WTSET(70) =    dis0                /fyfs/crude
      WTSET(71) =    dis1                /fyfs/crude
      WTSET(72) =    dis2                /fyfs/crude
!
      END

      FUNCTION DIX1(p1,p2,q1,q2)
!     ***************************
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION p1(4),p2(4),q1(4),q2(4)
! ------------------ Inline functions ------------------
      pdot(x1,y1,z1,e1,x2,y2,z2,e2) = 2*(e1*e2-x1*x2-y1*y2-z1*z2)
!
      s  = pdot(p1(1),p1(2),p1(3),p1(4) ,q1(1),q1(2),q1(3),q1(4))
      s1 = pdot(p2(1),p2(2),p2(3),p2(4) ,q2(1),q2(2),q2(3),q2(4))
      u  = pdot(p1(1),p1(2),p1(3),p1(4) ,q2(1),q2(2),q2(3),q2(4))
      u1 = pdot(p2(1),p2(2),p2(3),p2(4) ,q1(1),q1(2),q1(3),q1(4))
      dix1   = (s**2+u**2+s1**2+u1**2)/(4*s**2)
      END

      FUNCTION DIX2(noga,ph,p1,p2,q1,q2)
!     **********************************
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION ph(4),p1(4),p2(4),q1(4),q2(4)
      DIMENSION w1(4),w2(4),v1(4),v2(4)
! ------------------ Inline functions ------------------
      pdot(x1,y1,z1,e1,x2,y2,z2,e2) = 2*(e1*e2-x1*x2-y1*y2-z1*z2)
!
      IF(noga.eq.1) THEN
        DO 1 k=1,4
 1      w1(k) = p1(k) - ph(k)
        s  = pdot(w1(1),w1(2),w1(3),w1(4) ,q1(1),q1(2),q1(3),q1(4))
        s1 = pdot(p2(1),p2(2),p2(3),p2(4) ,q2(1),q2(2),q2(3),q2(4))
        u  = pdot(w1(1),w1(2),w1(3),w1(4) ,q2(1),q2(2),q2(3),q2(4))
        u1 = pdot(p2(1),p2(2),p2(3),p2(4) ,q1(1),q1(2),q1(3),q1(4))
      ELSEIF(noga.eq.2) THEN
        DO 2 k=1,4
 2      w2(k) = p2(k) + ph(k)
        s  = pdot(p1(1),p1(2),p1(3),p1(4) ,q1(1),q1(2),q1(3),q1(4))
        s1 = pdot(w2(1),w2(2),w2(3),w2(4) ,q2(1),q2(2),q2(3),q2(4))
        u  = pdot(p1(1),p1(2),p1(3),p1(4) ,q2(1),q2(2),q2(3),q2(4))
        u1 = pdot(w2(1),w2(2),w2(3),w2(4) ,q1(1),q1(2),q1(3),q1(4))
      ELSEIF(noga.eq.3) THEN
        DO 3 k=1,4
 3      v1(k) = q1(k) - ph(k)
        s  = pdot(p1(1),p1(2),p1(3),p1(4) ,v1(1),v1(2),v1(3),v1(4))
        s1 = pdot(p2(1),p2(2),p2(3),p2(4) ,q2(1),q2(2),q2(3),q2(4))
        u  = pdot(p1(1),p1(2),p1(3),p1(4) ,q2(1),q2(2),q2(3),q2(4))
        u1 = pdot(p2(1),p2(2),p2(3),p2(4) ,v1(1),v1(2),v1(3),v1(4))
      ELSEIF(noga.eq.4) THEN
        DO 4 k=1,4
 4      v2(k) = q2(k) + ph(k)
        s  = pdot(p1(1),p1(2),p1(3),p1(4) ,q1(1),q1(2),q1(3),q1(4))
        s1 = pdot(p2(1),p2(2),p2(3),p2(4) ,v2(1),v2(2),v2(3),v2(4))
        u  = pdot(p1(1),p1(2),p1(3),p1(4) ,v2(1),v2(2),v2(3),v2(4))
        u1 = pdot(p2(1),p2(2),p2(3),p2(4) ,q1(1),q1(2),q1(3),q1(4))
      ELSE
        write(6,*) 'DIX2: Wrong noga = ',noga
      ENDIF
      dix2   = (s**2+u**2+s1**2+u1**2)/(4*s**2)
      END

      SUBROUTINE MODL2B(MODE)
!     ***********************
!======================================================================!
!                                                                      !
!    NEW version  of the EXPONENTIATED matrix element, type (B)        !
!    Started May 92.                                                   !
!    O(alf2) part of the matrix element implemented Oct. 93            !
!                                                                      !
!======================================================================!
!     ************************
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER( PI = 3.1415926535897932D0, ALFINV=137.03604D0)
      PARAMETER( ALFA= 1D0/ALFINV, ALF1 = 1D0/PI/ALFINV)
      COMMON / BHPAR1 / DEL,EPSCM,THMIN,XMIVIS
      COMMON / BHPAR2 / CMSENE,AMEL
      COMMON / BHPAR3 / KEYRAD,KEYOPT
      COMMON / BHPARZ / AMAZ,GAMMZ,SINW2,GV,GA
      COMMON / TRANSR / TRAN,TRMIN,TRMAX
      COMMON / MOMS1  / TRANP,P1(4),P2(4),PHOT1(100,4),PHSU1(4),NPHOT1
      COMMON / MOMS2  / TRANQ,Q1(4),Q2(4),PHOT2(100,4),PHSU2(4),NPHOT2
      COMMON / MOMZ1  / AL1(50),BE1(50),WTM1(50),MK1(50)
      COMMON / MOMZ2  / AL2(50),BE2(50),WTM2(50),MK2(50)
      COMMON / PSIPHI / TH1,EXT1,EXB1,PSI1,EXW1,
     $                  TH2,EXT2,EXB2,PSI2,EXW2,EXE2,FIF,PHI
      COMMON / WGTALL / WTMOD,WTCRU1,WTCRU2,WTSET(300)
      SAVE   / BHPAR2 /, / BHPAR3 /, / TRANSR /
      SAVE   / MOMS1  /, / MOMS2  /, / MOMZ1  /, / MOMZ2  /
      SAVE   / WGTALL /, / BHPARZ /, / PSIPHI /
      DIMENSION sfu(100),bu10(100),bu11(100),hsu(100)
      DIMENSION sfl(100),bl10(100),bl11(100),hsl(100)
      DIMENSION bu11g(100)
      DIMENSION bl11g(100)
! ------------------ Inline functions ------------------
! Dot products
! Elements of single bremss. distribution
      CHI(x)= (1+(1-x)**2)/2
      XDELS(a,b,delt) =
     $ delt*(a**2+b**2)*(a**2+b**2)/((1-a)**2+(1-b)**2)/(a*b)
      XRHOR(a,b,psi) =(1-a)*(1-b)+a*b
     $  +2*sqrt(a*b*max(1-a,0d0)*max(1-b,0d0))*cos(psi)
      XCHIR(a,b,rho,zeta) = ( 1      +(1-zeta/max(1-a,zeta))**2
     $                       +rho**2 +(rho -(1-b)*zeta)**2 )/4
!********** alternative ***********
!*     XSFCM(a,b,aa,bb,y,z,delt) = a*b/(y*z)**2
!*    $ +delt/(y*z)**2 *(y**2+z**2)*(a**2+b**2)/((1-aa)**2+(1-bb)**2)
!*     dis1   = facpq* xchir(a,b,xrhor(a,b,psi1),zeta)
!*    $              * xsfcm(a,b,a,b,y,z,deltp)
!**********************************

!==================================================================!
      s = cmsene**2
      zeta = abs(tran)/s
!*****************************************************
!  Crude MC distribution (S-tilde factors omitted)
      CRUDE  =  s**2/tran**2
      facpq  =  s**2/(tranp*tranq)
      zeta   =  abs(tran)/s
      chi0   =  ( 1+(1-zeta)**2 )/2
      delta  =  amel**2/tran
      deltp  =  amel**2/tranp
      deltq  =  amel**2/tranq
      gam    =  2*alf1 *( dlog(1/delta)-1 )
      gamp   =  2*alf1 *( dlog(1/deltp)-1 )
      gamq   =  2*alf1 *( dlog(1/deltq)-1 )
      psip   = psi1
      psiq   = psi2
!==================================================================!
!    ///////////////////////////////////////
!    /////     UPPER + LOWER line      /////
!    ///////////////////////////////////////
! ------------------
      bt10u =0
      bt10l =0
      bt11u =0
      bt11l =0
      bt20u =0
      bt20l =0
      bt2ul =0
! DEBUG
      bt11ug =0
      bt11lg =0
! ------------------
! beta0
! ------------------
      bt00   =  chi0* s**2 /(tranp*tranq)
      bt01   =  bt00*(1+gam)
      bt02   =  bt00*(1+gam+0.5*gam**2)
! ------------------
! beta1 upper line
! ------------------
      DO i=1,nphot1
      a   = al1(i)
      b   = be1(i)
      y   = a +b*deltp
      z   = b +a*deltp
      v   = a +b -a*b
! soft factor precisely as in mass weight from the generator
      sfc  = a*b/(y*z)**2
      sfu(i) = sfc
! one photon bremss. distribution
      chir    = xchir(a,b,xrhor(a,b,psip),zeta)
      dis10   = facpq * chir *(1+xdels(a,b,deltp)) * sfc
      dis11g  = dis10 *(1+gam + gam/2*log((1-b)/(1-a))) ! with LL bug
      dis11   = dis10 *(1+gam + gam/4*log((1-b)/(1-a)))
c[[[[[
c>!#### dis11   = dis10 *(1+gam + gam/2*log((1-b)**2/(1-v)))
c>      dis11   = dis10 *(1+gam + gam/4*log((1-b)**2/(1-v)))
c]]]]]
      hsu(i) = dis10 / facpq
! beta1 O(alf1),O(alf2)
      bu10(i) =  dis10   -bt00*sfc
      bu11(i) =  dis11   -bt01*sfc
      bt10u   =  bt10u   +bu10(i)/sfc
      bt11u   =  bt11u   +bu11(i)/sfc
! DEBUG
      bu11g(i) =  dis11g   -bt01*sfc ! with LL bug
      bt11ug   =  bt11ug   +bu11g(i)/sfc
      ENDDO
! ------------------
! beta1 lower line
! ------------------
      DO i=1,nphot2
      a   = al2(i)
      b   = be2(i)
      y   = a +b*deltq
      z   = b +a*deltq
      v   = a +b -a*b
! soft factor precisely as in mass weight from the generator
      sfc  = a*b/(y*z)**2
      sfl(i) = sfc
! one photon bremss. distribution
      chir    = xchir(a,b,xrhor(a,b,psiq),zeta)
      dis10   = facpq * chir *(1+xdels(a,b,deltq)) * sfc
      dis11g  = dis10 *(1+gam + gam/2*log((1-b)/(1-a)))
      dis11   = dis10 *(1+gam + gam/4*log((1-b)/(1-a)))
c[[[[[
c>!#### dis11   = dis10 *(1+gam + gam/2*log((1-b)**2/(1-v)))
c>      dis11   = dis10 *(1+gam + gam/4*log((1-b)**2/(1-v)))
c]]]]]
      hsl(i) = dis10 / facpq
! beta1 O(alf1),O(alf2)
      bl10(i) =  dis10   -bt00*sfc
      bl11(i) =  dis11   -bt01*sfc
      bt10l   =  bt10l   +bl10(i)/sfc
      bt11l   =  bt11l   +bl11(i)/sfc
! DEBUG
      bl11g(i) =  dis11g   -bt01*sfc ! with LL bug
      bt11lg   =  bt11lg   +bl11g(i)/sfc
      ENDDO
!------------------
! beta2 upper line
!------------------
      do 153 i=1,nphot1
      do 153 j=i+1,nphot1
! Basic variables first photon
      a1  = al1(i)
      b1  = be1(i)
      y1  = a1 +b1*deltp
      z1  = b1 +a1*deltp
      v1  = a1+b1-a1*b1
! Basic variables second photon
      a2  = al1(j)
      b2  = be1(j)
      y2  = a2 +b2*deltp
      z2  = b2 +a2*deltp
      v2  = a2+b2-a2*b2
! Soft factors, starred variables
      sf1  =  sfu(i)
      sf2  =  sfu(j)
      a1st = a1/(1-v2)
      b1st = b1/(1-v2)
      v1st = v1/(1-v2)
      a2st = a2/(1-v1)
      b2st = b2/(1-v1)
      v2st = v2/(1-v1)
      a1st = min(a1st,1-zeta)
      b1st = min(b1st,1-zeta)
      v1st = min(v1st,1-zeta)
      a2st = min(a2st,1-zeta)
      b2st = min(b2st,1-zeta)
      v2st = min(v2st,1-zeta)
! H-functions exact and Leading Log
      gs1    = sf1*chi(v1)
      gs2    = sf2*chi(v2)
      hs1    = hsu(i)
      hs2    = hsu(j)
      gs1st  = sf1*chi(v1st)
      gs2st  = sf2*chi(v2st)
      hs1st  = xchir(a1st,b1st,xrhor(a1st,b1st,psip),zeta)
     $            *(1+xdels(a1st,b1st,deltp))* sf1
      hs2st  = xchir(a2st,b2st,xrhor(a2st,b2st,psip),zeta)
     $            *(1+xdels(a2st,b2st,deltp))* sf2
! The actual two photon distributions
      IF ( (a1-b1)*(a2-b2).le.0d0 ) THEN
!       OPPOSITE directions two photons
        IF (a1.gt.b1) THEN
!          FIRST photon in INITIAL state, second in the final
           IF(v1.gt.v2) THEN
               d2uu  = facpq  *hs1st*gs2
           ELSE
               d2uu  = facpq  *gs1st*hs2
           ENDIF
        ELSE
!          SECOND photon in INITIAL state, first in the final
           IF(v1.gt.v2) THEN
               d2uu  = facpq  *hs1*gs2st
           ELSE
               d2uu  = facpq  *gs1*hs2st
           ENDIF
        ENDIF
      ELSE
!       SAME directions two photons
        IF(v1.gt.v2) THEN
            d2uu  = facpq/2 *( hs1*gs2st + hs1st*gs2)
        ELSE
            d2uu  = facpq/2 *( hs2*gs1st + hs2st*gs1)
        ENDIF
      ENDIF
      bu20  = d2uu-bu10(i)*sfu(j)-bu10(j)*sfu(i)-bt00*sfu(j)*sfu(i)
      bt20u = bt20u + bu20/sfu(i)/sfu(j)
 153  CONTINUE
! ------------------
! beta2 lower line
! ------------------
      DO 158 i=1,nphot2
      DO 158 j=i+1,nphot2
! Basic variables first photon
      a1  = al2(i)
      b1  = be2(i)
      y1  = a1 +b1*deltq
      z1  = b1 +a1*deltq
      v1  = a1+b1-a1*b1
! Basic variables second photon
      a2  = al2(j)
      b2  = be2(j)
      y2  = a2 +b2*deltq
      z2  = b2 +a2*deltq
      v2  = a2+b2-a2*b2
! soft factors, starred variables
      sf1  =  sfl(i)
      sf2  =  sfl(j)
      a1st = a1/(1-v2)
      b1st = b1/(1-v2)
      v1st = v1/(1-v2)
      a2st = a2/(1-v1)
      b2st = b2/(1-v1)
      v2st = v2/(1-v1)
      a1st = min(a1st,1-zeta)
      b1st = min(b1st,1-zeta)
      v1st = min(v1st,1-zeta)
      a2st = min(a2st,1-zeta)
      b2st = min(b2st,1-zeta)
      v2st = min(v2st,1-zeta)
! H-functions exact and Leading Log
      gs1    = sf1*chi(v1)
      gs2    = sf2*chi(v2)
      hs1    = hsl(i)
      hs2    = hsl(j)
      gs1st  = sf1*chi(v1st)
      gs2st  = sf2*chi(v2st)
      hs1st  = xchir(a1st,b1st,xrhor(a1st,b1st,psiq),zeta)
     $            *(1+xdels(a1st,b1st,deltq))* sf1
      hs2st  = xchir(a2st,b2st,xrhor(a2st,b2st,psiq),zeta)
     $            *(1+xdels(a2st,b2st,deltq))* sf2
! The actual two photon distributions
      IF ( (a1-b1)*(a2-b2).le.0d0 ) THEN
!       OPPOSITE directions two photons
        IF (a1.gt.b1) THEN
!          FIRST photon in INITIAL state, second in the final
           IF(v1.gt.v2) THEN
               d2ll  = facpq  *hs1st*gs2
           ELSE
               d2ll  = facpq  *gs1st*hs2
           ENDIF
        ELSE
!          SECOND photon in INITIAL state, first in the final
           IF(v1.gt.v2) THEN
               d2ll  = facpq  *hs1*gs2st
           ELSE
               d2ll  = facpq  *gs1*hs2st
           ENDIF
        ENDIF
      ELSE
!       SAME directions two photons
        IF(v1.gt.v2) THEN
           d2ll= facpq/2 *( hs1*gs2st  + hs1st*gs2)
       ELSE
           d2ll= facpq/2 *( hs2*gs1st  + hs2st*gs1)
        ENDIF
      ENDIF
      bl20  = d2ll-bl10(i)*sfl(j) -bl10(j)*sfl(i) -bt00*sfl(j)*sfl(i)
      bt20l = bt20l + bl20/sfl(i)/sfl(j)
 158  CONTINUE
!----------------------------
! beta2 upper and lower line
!----------------------------
      DO 160 i=1,nphot1
      DO 160 j=1,nphot2
! Basic variables first photon
      a1  = al1(i)
      b1  = be1(i)
      y1  = a1 +b1*deltp
      z1  = b1 +a1*deltp
      v1  = a1 +b1-a1*b1
! Basic variables second photon
      a2  = al2(j)
      b2  = be2(j)
      y2  = a2 +b2*deltq
      z2  = b2 +a2*deltq
      v2  = a2+b2-a2*b2
! soft factors
      sf1  =  sfu(i)
      sf2  =  sfl(j)
      IF(v1.gt.v2) THEN
         d2ul= facpq  *hsu(i)       *sf2*chi(v2)
       ELSE
         d2ul= facpq  *sf1*chi(v1)  *hsl(j)
      ENDIF
      b2ul  = d2ul-bu10(i)*sfl(j)-sfu(i)*bl10(j)-bt00*sfu(i)*sfl(j)
      bt2ul = bt2ul + b2ul/sfu(i)/sfl(j)
  160 CONTINUE

!        **************************************
!        ****   Definitions of MC weights  ****
!        **************************************
!          ---------------------------------
!          !    UPPER line ONLY,  New      !
!          ---------------------------------
! Case of upper line only has separate entries in WTSET because it is 
! analysed in the separate run with separate plotting programs 
! with different entries in vvrho, differerent histogram names etc.
!
! All beta's:  TOTAL O(alf0),O(alf1),O(alf2)
      WTSET(130) =    bt00                /crude
      WTSET(131) =   (bt01+bt10u)         /crude
      WTSET(132) =   (bt02+bt11u +bt20u)  /crude
! Individual beta's in various orders.
! Special O(alf0) bet0 for tests only (without chi0) called WTZERO
      WTSET(134) =   ( s**2/(tranp*tranq) )/crude
! O(alf1)
      WTSET(135) =   bt01  /crude                  ! beta0
      WTSET(136) =   bt10u /crude                  ! beta1
! O(alf2)
      WTSET(137) =   bt02  /crude                  ! beta0
      WTSET(138) =   bt11u /crude                  ! beta1
      WTSET(139) =   bt20u /crude                  ! beta2
! DEBUG
      WTSET(232) =   (bt02+bt11ug +bt20u)  /crude  ! total (LL bug)
      WTSET(238) =   bt11u /crude                  ! beta1 (LL bug)
!          ---------------------------------
!          !  UPPER + LOWER line, New      !
!          ---------------------------------
! All beta's:  O(alf0),O(alf1),O(alf2)
      WTSET(140) =    bt00  /crude
      WTSET(141) =   (bt01+bt10u+bt10l) /crude
      WTSET(142) =   (bt02+bt11u+bt11l+bt20u+bt20l+bt2ul) /crude
! Individual beta's in various orders.
! O(alf1)
      WTSET(143) =    bt01   /crude                ! beta0
      WTSET(144) =   (bt10u+bt10l) /crude          ! beta1
      WTSET(145) =    bt10u  /crude                ! beta1 components
      WTSET(146) =    bt10l  /crude
! O(alf2)
      WTSET(147) =    bt02   /crude                ! beta0
      WTSET(148) =   (bt11u+bt11l)/crude           ! beta1
      WTSET(149) =   (bt20u+bt20l+bt2ul)/crude     ! beta2
      WTSET(150) =    bt11u  /crude                ! beta1 components
      WTSET(151) =    bt11l  /crude
      WTSET(152) =    bt2ul  /crude
      WTSET(153) =    bt20u  /crude                ! beta2 components
      WTSET(154) =    bt20l  /crude
! DEBUG
      WTSET(242) = (bt02+bt11ug+bt11lg+bt20u+bt20l+bt2ul) /crude
      WTSET(250) =    bt11ug  /crude  ! (LL bug)
      WTSET(251) =    bt11l   /crude  ! (LL bug)
!==================================================================
!==================================================================
!                  Non-exponentiated version                      !
!==================================================================
!==================================================================
      pdel  = del*bcud(p1,p2,phsu1)
      qdel  = del*bcud(q1,q2,phsu2)
      fyfsu = exp(-gamp*dlog(1/pdel)  +gamp/4 -alf1/2)
      fyfsl = exp(-gamq*dlog(1/qdel)  +gamq/4 -alf1/2)
! temporary test CONV distr.
      bcon0u    =0
      bcon1u    =0
      bcon2u    =0
! entire distributions
      dis1uC   =0
      dis0u    =0
      dis1u    =0
      dis2u    =0
! beta-like components
      bt00u    =0
      bt01u    =0
      bt02u    =0
      bt10u    =0
      bt11u    =0
      bt20u    =0
! DEBUG part
      dis2ug   = 0
      bt11ug   = 0
      beta11u  = 0
      beta10u  = 0
      beta11ug = 0
!------------------------------------
! Normal one-line virtual correction
      virt1   =  3d0/4*gam -1d0/2*alf1 
! Tweo lines with one YFS formfactor subtraction
! that is: 5/4gam = 2(3/4)gam -1/4gam
      virt1x = 5d0/4*gam-1d0/2*alf1
      virt2x = 0.5d0*(5d0/4*gam)**2
! Note that in order to shorten our formulas we omitt S-factors,
! They always cancel with S-factors from crude distributions
!    ///////////////////////////////////////
!    /////  UNEXP  UPPER line only     /////
!    ///////////////////////////////////////
      IF(nphot1.eq.0) THEN
!     ====================
! [00] No photons, note here gamp=gam
!------
! Clasical O(alf0,1)
        dis0uC  = facpq*chi0 
        dis1uC  = facpq*chi0 *(1 +virt1 +gamp*dlog(del))
! NCONV  
        bcon0u= facpq*chi0 
        bcon1u= facpq*chi0*(1 +gamp*dlog(del))
        bcon2u= facpq*chi0*(1 +gamp*dlog(del) +1./2*(gamp*dlog(del))**2)
! O(alf0,1,2) Entire distribution 
        dis0u = facpq*chi0
        dis1u = facpq*chi0*(1 +gamp*dlog(del))
     $         +facpq*chi0*virt1x
        dis2u = facpq*chi0*(1+gamp*dlog(del) +1./2*(gamp*dlog(del))**2)
     $         +facpq*chi0*(1+gamp*dlog(del))*virt1x
     $         +facpq*chi0*virt2x
!------
! Beta's for controll
        bt00u= dis0u
        bt01u= dis1u
        bt02u= dis2u
! DEBUG part
        dis2ug = dis2u
      ELSEIF(nphot1.eq.1) THEN
!     ========================
! [10] One upper line photon, note here gamp.ne.gam !!!
!------
! Basic distribution
        a   = al1(1)
        b   = be1(1)
        v   = a+b-a*b
! O(alf1) Entire distributions 
        hs1  = xchir(a,b,xrhor(a,b,psip),zeta)*(1+xdels(a,b,deltp))
!------
! Clasical O(alf1)
        dis1uC = facpq *hs1
! NCONV 
        bcon1u   = facpq*chi0
        bcon2u   = facpq*chi0 *(1 +gamp*dlog(del))
!------
! O(alf1,2) Entire distribution
        dis1u  = facpq*hs1
        dis2ug = facpq*hs1*(1 +gamp*dlog(del) +virt1x 
     $       -1/2d0*gam*log(1-a) -1/2d0*gam*log(1-b) )
        dis2u  = facpq*hs1*(1 +gamp*dlog(del) +virt1x 
     $       -1/4d0*gam*log(1-a) -3/4d0*gam*log(1-b) )
c[[[[[[
c>!####   dis2u  = facpq*hs1*(1 +gamp*dlog(del) +virt1x -gam/2*log(1-v))
c>        dis2u  = facpq*hs1*(1 +gamp*dlog(del) +virt1x 
c>     $       -gam/4*log(1-v)-gam/2*log(1-b) )
c]]]]]]
! Beta's for control
        bt01u  = facpq*chi0
        bt02u  = facpq*chi0 *(1 +gamp*dlog(del) +virt1x)
        bt10u  = dis1u  -bt01u
        bt11u  = dis2u  -bt02u
! [[[[[[[[[[ DEBUG part [[[[[[[[[[[[[[[[[[[[[[[[[[[
        bt11ug = dis2ug -bt02u ! with LL bug
! Standard YFS definition of beta1. Note that bt11u includes also
! virtual gamp*dlog(del) It cancels bt10u contribution from the
! two photon configuration!
        delvir1   =  gamp*dlog(del) +3d0/4*gam -1d0/2*alf1 
        delvirB   =  gamp*dlog(del) +1d0/4*gam -1d0/2*alf1 
        D10u  = facpq*hs1
        D11u  = facpq*hs1*(1 +delvir1
     $                      -1/4d0*gam*log(1-a) -3/4d0*gam*log(1-b) )
        D11ug = facpq*hs1*(1 +delvir1
     $                      -1/2d0*gam*log(1-a) -1/2d0*gam*log(1-b) )
        beta00u  = facpq*chi0
        beta01u  = facpq*chi0*(1+delvir1-delvirB)
        beta10u  = D10u -beta00u
        beta11u  = D11u  -delvirB*D10u -beta01u
        beta11ug = D11ug -delvirB*D10u -beta01u ! with LL bug
! ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]
      ELSEIF(nphot1.eq.2) THEN
!     ========================
! [20] Two upper line photon 
!------
! NCONV 
        bcon2u   = facpq*chi0
!------
! Basic distribution
        a1  = al1(1)
        b1  = be1(1)
        a2  = al1(2)
        b2  = be1(2)
        v1  = a1+b1-a1*b1
        v2  = a2+b2-a2*b2
        a1st = min(a1/(1-v2),1-zeta)
        b1st = min(b1/(1-v2),1-zeta)
        a2st = min(a2/(1-v1),1-zeta)
        b2st = min(b2/(1-v1),1-zeta)
        v1st = min(v1/(1-v2),1-zeta)
        v2st = min(v2/(1-v1),1-zeta)
!       H-functions exact and Leading Log
        hs1 =xchir(a1,b1,xrhor(a1,b1,psip),zeta)*(1+xdels(a1,b1,deltp))
        hs2 =xchir(a2,b2,xrhor(a2,b2,psip),zeta)*(1+xdels(a2,b2,deltp))
        hs1st  = xchir(a1st,b1st,xrhor(a1st,b1st,psip),zeta)
     $            *(1+xdels(a1st,b1st,deltp))
        hs2st  = xchir(a2st,b2st,xrhor(a2st,b2st,psip),zeta)
     $            *(1+xdels(a2st,b2st,deltp))
!       The actual two photon distributions
        IF ( (a1-b1)*(a2-b2).le.0d0 ) THEN
!       OPPOSITE directions two photons
          IF (a1.gt.b1) THEN
!          FIRST photon in INITIAL state, second in the final
             IF(v1.gt.v2) d2u  = hs1st*chi(v2)
             IF(v2.ge.v1) d2u  = chi(v1st)*hs2
          ELSE
!          SECOND photon in INITIAL state, first in the final
             IF(v1.gt.v2) d2u  = hs1*chi(v2st)
             IF(v2.ge.v1) d2u  = chi(v1)*hs2st
          ENDIF
        ELSE
!       SAME directions two photons
          IF(v1.gt.v2) d2u  = 1./2 *( hs1*chi(v2st) + hs1st*chi(v2))
          IF(v2.ge.v1) d2u  = 1./2 *( hs2*chi(v1st) + hs2st*chi(v1))
        ENDIF
!------
! O(alf2) Entire distribution dis2u 
        dis2u = facpq*d2u
!------
! Beta's for control
        bt02u = facpq*chi0
        bt11u = facpq*( (hs1-chi0) +(hs2-chi0) )
        bt20u = dis2u -bt11u  -bt02u
! DEBUG part
        dis2ug = dis2u
        bt11ug = bt11u ! with LL bug
      ENDIF
!     =====
!        **************************************
!        ****   Definitions of MC weights  ****
!        ****    Non-exponentiated case    ****
!        **************************************
! Case of upper-line-only has separate entries in WTSET because it is 
! analysed in the separate run with separate plotting programs 
! with different entries in vvrho, differerent histogram names etc.
!
! YFS formfactor fyfsu has to be removed in the unexponentiated case!
!-----
! NCONV for relative comparisons
      WTSET(156) =   bcon1u      /fyfsu/crude ! NCONV O(alf1)
      WTSET(157) =   bcon2u      /fyfsu/crude ! NCONV O(alf2)
!          ---------------------------------
!          /////    UPPER line ONLY    /////
!          ---------------------------------
! TOTAL O(alf0),O(alf1),O(alf2)
      WTSET(160) =   bt00u       /fyfsu/crude
      WTSET(161) =   dis1u       /fyfsu/crude
      WTSET(162) =   dis2u       /fyfsu/crude
      WTSET(164) =   dis1uC      /fyfsu/crude
! Individual beta's in various orders.
! O(alf1)
      WTSET(165) =   bt01u  /fyfsu/crude            ! bt0
      WTSET(166) =   bt10u  /fyfsu/crude            ! bt1
! O(alf2)
      WTSET(167) =   bt02u  /fyfsu/crude            ! bt0
      WTSET(168) =   bt11u  /fyfsu/crude            ! bt1
      WTSET(169) =   bt20u  /fyfsu/crude            ! bt2
! Debug part
      WTSET(262) =   dis2ug   /fyfsu/crude    ! O(alf2) tot BUG
      WTSET(266) =   beta10u  /fyfsu/crude    ! O(alf1) beta1 true 
      WTSET(268) =   beta11u  /fyfsu/crude    ! O(alf2) beta1 true
      WTSET(298) =   beta11ug /fyfsu/crude    ! O(alf2) beta1 bug
      WTSET(290) =   bt11ug   /fyfsu/crude    ! O(alf2) bt1 BUG

!    ////////////////////////////////////////////
!    /////  UNEXP   UPPER + LOWER line      /////
!    ////////////////////////////////////////////
! CONV for primary tests
      bcon0  =0
      bcon1  =0
      bcon2  =0
! Entire 0,1,2-photon distributions
      dis0   =0
      dis1   =0
      dis2   =0
! beta-like components
      bt00    =0
      bt01    =0
      bt02    =0
      bt10    =0
      bt11    =0
      bt20    =0
! ---------- salami on second order beta1
      bt11rg =0
      bt11uu =0
      bt11ll =0
      bt11ul =0
! ----------salami on second order beta2
      bt20u =0
      bt20l =0
      bt2ul =0
! DEBUG
      dis2g = 0
!
      virt1y = 6./4*gam -alf1
      virt2y = 0.5*(6./4*gam)**2
!=================================================
      IF(nphot1.eq.0.and.nphot2.eq.0) THEN
! [00] No photons
!------CONV 
        bcon0=facpq*chi0
        bcon1=bcon0*(1 -2*gam*dlog(1/del))
        bcon2=bcon0*(1 -2*gam*dlog(1/del)+0.5*(2*gam*dlog(1/del))**2 )
!------
! O(alf0,1,2) Entire distribution 
        dis0= facpq*chi0
        dis1= facpq*chi0*(1+2*gam*dlog(del))
     $       +facpq*chi0*virt1y
        dis2= facpq*chi0*(1+2*gam*dlog(del)+0.5*(2*gam*dlog(del))**2)
     $       +facpq*chi0*(1+2*gam*dlog(del))*virt1y
     $       +facpq*chi0* virt2y
!------
! Beta's for controll
        bt00= dis0
        bt01= dis1
        bt02= dis2
! DEBUG
        dis2g = dis2
!=================================================
      ELSEIF(nphot1.eq.1.and.nphot2.eq.0) THEN
! [10] One upper line photon 
!------CONV
        bcon1  = facpq*chi0
        bcon2  = bcon1  *(1 -(gamp+gamq)*dlog(1/del))
!------Basic distribution
        hs1 = (1+xdels(al1(1),be1(1),deltp))
     $        *xchir(al1(1),be1(1),xrhor(al1(1),be1(1),psip),zeta)
!------
! O(alf1,2) Entire distribution
! Note that 1/2 in the log(1-v) term = {1/2 from two tree diagrams} 
        dis1 = facpq*hs1
        dis2g= facpq*hs1*(1 +(gamp+gamq)*dlog(del) +virt1y
     $       -1/2d0*gam*log(1-al1(1)) -1/2d0*gam*log(1-be1(1)) )
        dis2 = facpq*hs1*(1 +(gamp+gamq)*dlog(del) +virt1y
     $       -1/4d0*gam*log(1-al1(1)) -3/4d0*gam*log(1-be1(1)) )
c[[[[[[
c>        dis2= facpq*hs1*(1 +(gamp+gamq)*dlog(del) 
c>!####$       +virt1y -gam/2*log(1-v) )
c>     $       +virt1y -gam/4*log(1-v) -gam/2*log(1-be1(1)) )
c]]]]]]
! Beta's for control
        bt01 = facpq*chi0
        bt02 = facpq*chi0 *(1 +(gamp+gamq)*dlog(del) +virt1y)
        bt10 = dis1  -bt01
        bt11 = dis2  -bt02
! DEBUG part
        bt11g= dis2g -bt02 ! with LL bug
!/////////  SALAMI on bet11 ///////// 
! bet11 rewritten as bt10 +bt11rg +bt11uu +bt11ll +bt11ul
!        bt11 = facpq*(
!     $  (hs1-chi0)                            ! first order
!     $ +(hs1-chi0)*virt1y -hs1*gam/2*log(1-v) ! regular part
!     $ +(hs1-chi0)*gamp*dlog(del)             ! re_uper-im_uper
!     $ +(hs1-chi0)*gamq*dlog(del)             ! re_uper-im_lower
!     $)
        bt11rg = facpq*((hs1-chi0)*virt1y -hs1*gam/2*log(1-v))
        bt11uu = facpq*(hs1-chi0)*gamp*dlog(del)
        bt11ul = facpq*(hs1-chi0)*gamq*dlog(del)

!=================================================
      ELSEIF(nphot1.eq.0.and.nphot2.eq.1) THEN
! [01] One lower line photon 
!------CONV 
        bcon1  = facpq*chi0
        bcon2  = bcon1  *(1 -(gamp+gamq)*dlog(1/del))
!------Basic distribution
        hs1 = (1+xdels(al2(1),be2(1),deltq))
     $        *xchir(al2(1),be2(1),xrhor(al2(1),be2(1),psiq),zeta)
!------
! O(alf1,2) Entire distribution
! Note that 1/2 in the log(1-v) term = {1/2 from two tree diagrams} 
        dis1 = facpq*hs1
        dis2g= facpq*hs1*(1 +(gamp+gamq)*dlog(del) +virt1y
     $       -1/2d0*gam*log(1-al2(1)) -1/2d0*gam*log(1-be2(1)) )
        dis2 = facpq*hs1*(1 +(gamp+gamq)*dlog(del) +virt1y
     $       -1/4d0*gam*log(1-al2(1)) -3/4d0*gam*log(1-be2(1)) )
![[[[
c>        dis2= facpq*hs1*(1 +(gamp+gamq)*dlog(del) 
c>!####$       +virt1y -gam/2*log(1-v) )
c>     $       +virt1y -gam/4*log(1-v) -gam/2*log(1-be2(1)) )
c]]]]
! Beta's for control
        bt01 = facpq*chi0
        bt02 = facpq*chi0 *(1 +(gamp+gamq)*dlog(del) +virt1y)
        bt10 = dis1  -bt01
        bt11g= dis2g -bt02 ! with LL bug
        bt11 = dis2  -bt02
!/////////  SALAMI on bet11 ///////// 
        bt11rg = facpq*((hs1-chi0)*virt1y -hs1*gam/2*log(1-v))
        bt11ll = facpq*(hs1-chi0)*gamq*dlog(del)
        bt11ul = facpq*(hs1-chi0)*gamp*dlog(del)
!=================================================
      ELSEIF(nphot1.eq.2.and.nphot2.eq.0) THEN
! [20] Two upper line photons 
!------CONV
        bcon2  = facpq*chi0
!------
! Basic distribution
        a1  = al1(1)
        b1  = be1(1)
        a2  = al1(2)
        b2  = be1(2)
        v1  = a1+b1-a1*b1
        v2  = a2+b2-a2*b2
        a1st = min(a1/(1-v2),1-zeta)
        b1st = min(b1/(1-v2),1-zeta)
        a2st = min(a2/(1-v1),1-zeta)
        b2st = min(b2/(1-v1),1-zeta)
        v1st = min(v1/(1-v2),1-zeta)
        v2st = min(v2/(1-v1),1-zeta)
!       H-functions exact and Leading Log
        hs1 =xchir(a1,b1,xrhor(a1,b1,psip),zeta)*(1+xdels(a1,b1,deltp))
        hs2 =xchir(a2,b2,xrhor(a2,b2,psip),zeta)*(1+xdels(a2,b2,deltp))
        hs1st  = xchir(a1st,b1st,xrhor(a1st,b1st,psip),zeta)
     $            *(1+xdels(a1st,b1st,deltp))
        hs2st  = xchir(a2st,b2st,xrhor(a2st,b2st,psip),zeta)
     $            *(1+xdels(a2st,b2st,deltp))
!       The actual two photon distributions
        IF ( (a1-b1)*(a2-b2).le.0d0 ) THEN
!       OPPOSITE directions two photons
          IF (a1.gt.b1) THEN
!          FIRST photon in INITIAL state, second in the final
             IF(v1.gt.v2) d2u  = hs1st*chi(v2)
             IF(v2.ge.v1) d2u  = chi(v1st)*hs2
          ELSE
!          SECOND photon in INITIAL state, first in the final
             IF(v1.gt.v2) d2u  = hs1*chi(v2st)
             IF(v2.ge.v1) d2u  = chi(v1)*hs2st
          ENDIF
        ELSE
!       SAME directions two photons
          IF(v1.gt.v2) d2u  = 1./2 *( hs1*chi(v2st) + hs1st*chi(v2))
          IF(v2.ge.v1) d2u  = 1./2 *( hs2*chi(v1st) + hs2st*chi(v1))
        ENDIF
!------
! O(alf2) Entire distribution dis2u 
        dis2 = facpq*d2u
!------
! Beta's for control
        bt02 = facpq*chi0
        bt11 = facpq*( (hs1-chi0) +(hs2-chi0) )
        bt20 = dis2 -bt11  -bt02
        bt20u= bt20
!/////////  SALAMI on bet11 ///////// 
        bt11uu = bt11
! DEBUG
        dis2g = dis2
!=================================================
      ELSEIF(nphot1.eq.0.and.nphot2.eq.2) THEN
! [02] Two lower line photons 
!------CONV
        bcon2  = facpq*chi0
!------
! Basic distribution
        a1  = al2(1)
        b1  = be2(1)
        a2  = al2(2)
        b2  = be2(2)
        v1  = a1+b1-a1*b1
        v2  = a2+b2-a2*b2
        a1st = min(a1/(1-v2),1-zeta)
        b1st = min(b1/(1-v2),1-zeta)
        a2st = min(a2/(1-v1),1-zeta)
        b2st = min(b2/(1-v1),1-zeta)
        v1st = min(v1/(1-v2),1-zeta)
        v2st = min(v2/(1-v1),1-zeta)
!       H-functions exact and Leading Log
        hs1 =xchir(a1,b1,xrhor(a1,b1,psiq),zeta)*(1+xdels(a1,b1,deltq))
        hs2 =xchir(a2,b2,xrhor(a2,b2,psiq),zeta)*(1+xdels(a2,b2,deltq))
        hs1st  = xchir(a1st,b1st,xrhor(a1st,b1st,psiq),zeta)
     $            *(1+xdels(a1st,b1st,deltq))
        hs2st  = xchir(a2st,b2st,xrhor(a2st,b2st,psiq),zeta)
     $            *(1+xdels(a2st,b2st,deltq))
!       The actual two photon distributions
        IF ( (a1-b1)*(a2-b2).le.0d0 ) THEN
!       OPPOSITE directions two photons
          IF (a1.gt.b1) THEN
!          FIRST photon in INITIAL state, second in the final
             IF(v1.gt.v2) d2u  = hs1st*chi(v2)
             IF(v2.ge.v1) d2u  = chi(v1st)*hs2
          ELSE
!          SECOND photon in INITIAL state, first in the final
             IF(v1.gt.v2) d2u  = hs1*chi(v2st)
             IF(v2.ge.v1) d2u  = chi(v1)*hs2st
          ENDIF
        ELSE
!       SAME directions two photons
          IF(v1.gt.v2) d2u  = 1./2 *( hs1*chi(v2st) + hs1st*chi(v2))
          IF(v2.ge.v1) d2u  = 1./2 *( hs2*chi(v1st) + hs2st*chi(v1))
        ENDIF
!------
! O(alf2) Entire distribution dis2u 
        dis2 = facpq*d2u
!------
! Beta's for control
        bt02 = facpq*chi0
        bt11 = facpq*( (hs1-chi0) +(hs2-chi0) )
        bt20 = dis2 -bt11  -bt02
        bt20l= bt20
c/////////  SALAMI on bet11 ///////// 
        bt11ll = bt11
! DEBUG
        dis2g = dis2
!=================================================
      ELSEIF(nphot1.eq.1.and.nphot2.eq.1) THEN
! [11] One upper and one lower line photons 
!------CONV
        bcon2  = facpq*chi0
!------
! O(alf2) Entire distribution dis2u 
      a1  = al1(1)
      b1  = be1(1)
      a2  = al2(1)
      b2  = be2(1)
      v1  = a1+b1-a1*b1
      v2  = a2+b2-a2*b2
      hs1 = xchir(a1,b1,xrhor(a1,b1,psip),zeta)*(1+xdels(a1,b1,deltp))
      hs2 = xchir(a2,b2,xrhor(a2,b2,psiq),zeta)*(1+xdels(a2,b2,deltq))
      IF(v1.gt.v2) THEN
        d2ul=  hs1     *chi(v2)
      ELSE
        d2ul=  chi(v1) * hs2
      ENDIF
        dis2 = facpq*d2ul
!------
! Beta's for control
        bt02  = facpq*chi0
        bt11  = facpq*( (hs1-chi0) +(hs2-chi0) )
        bt20  = dis2 -bt11 -bt02
        bt2ul = bt20
!/////////  SALAMI on bet11 ///////// 
        bt11ul = bt11
! DEBUG
        dis2g = dis2
      ENDIF
!     =====

!        **************************************
!        ****   Definitions of MC weights  ****
!        ****    Non-exponentiated case    ****
!        **************************************
!          ---------------------------------
!          /////  UPPER + LOWER line   /////
!          ---------------------------------

      fyfs = fyfsu*fyfsl
! TOTAL O(alf0),O(alf1),O(alf2)
      WTSET(170) =    dis0                /fyfs/crude
      WTSET(171) =    dis1                /fyfs/crude
      WTSET(172) =    dis2                /fyfs/crude
!***************************************************************
! Individual beta's in various orders.
! O(alf1)
      WTSET(173) =    bt01         /fyfs/crude       ! beta0
      WTSET(174) =    bt10         /fyfs/crude       ! beta1
! O(alf2)
      WTSET(175) =    bt11uu       /fyfs/crude       ! beta11 up-up
      WTSET(176) =    bt11ll       /fyfs/crude       ! beta11 low-low
      WTSET(177) =    bt02         /fyfs/crude       ! beta0
      WTSET(178) =    bt11         /fyfs/crude       ! beta1
      WTSET(179) =    bt20         /fyfs/crude       ! beta2
      WTSET(180) =    bt11rg       /fyfs/crude       ! beta11 regular 
      WTSET(181) =    bt11ul       /fyfs/crude       ! beta11 up-low
      WTSET(182) =    bt2ul        /fyfs/crude
      WTSET(183) =    bt20u        /fyfs/crude       ! beta2 components
      WTSET(184) =    bt20l        /fyfs/crude
! CONV distribution for tests
      WTSET(185) =   bcon0         /fyfs/crude ! NCONV O(alf0)
      WTSET(186) =   bcon1         /fyfs/crude ! NCONV O(alf1)
      WTSET(187) =   bcon2         /fyfs/crude ! NCONV O(alf2)
! DEBUG part
      WTSET(272) =    dis2g        /fyfs/crude ! Total with LL bug
!***************************************************************
      END
