      Subroutine CompAZG2(Ein,Wmin,Wmax,E1,P1,E2,P2)
C
C Generate two photon energies from CompAZ parametrization
C     for given range in gamma-gamma invariant mass W_gg
C
C  Input:  
C  ======
C             Ein   - electron beam energy [GeV]  50.<=Abs(Ein)<=500.
C             Wmin  - minimum value of W_gg
C             Wmax  - maximum value og W_gg
C                     if Wmax<=0 -> up to the limit
C  Output: 
C  =======
C
C  For Ein>0  
C  ---------  E1,E2  - generated photon energies
C             P1,P2  - corresponding photon polarizations
C
C
C  For Ein<0  E1          -  1/L * Int(dL/dW) from Wmin to Wmax for Wmin<Wmax
C  ---------         or   -  1/L * dL/dW  [1/GeV] for Wmin=Wmax,
C                                                         for Ee=ABS(Ein)
C             E2,P1,P2    - not used 
C
C============================================================================= 
C
      Implicit None
C
      Real Ein,Wmin,Wmax,E1,E2,P1,P2
C
      Real Eold,Wold1,Wold2,Method
      Save Eold,Wold1,Wold2,Method
      Data Eold / -1. /
      Data Wold1,Wold2 / -1., -1. /
C
      Real xmax2,Wmax2,S0,Smin,Smax,ymax
      Save Wmax2,Smin,Smax
C
      Real Norm,frac,frac2,S,Ee
      Integer Ngen,Nevt
      Save Norm,Ngen,Nevt,frac2
C
      Real Compaz      
C
C----------------------------------------------------------------------
C
C Check energy
C
      Ee=ABS(Ein)

      If(Ee.LT.50. .OR. Ee.GT.500.)then

        write(6,*)'************************************************* '
        write(6,*)' '
        write(6,*)' '
        write(6,*)'CompAZG2  energy out of range (50.-500.) : ',Ee
        write(6,*)' '
        write(6,*)' '
        write(6,*)'************************************************* '

        STOP
C      ======

        Endif
C
C Check W range
C

      If(Wmin.LT.0. .OR. Wmax.LT.0. .OR.
     +  (Wmax.NE.0. .AND. Wmin.GT.Wmax) )then

        write(6,*)'************************************************* '
        write(6,*)' '
        write(6,*)' '
        write(6,*)'CompAZG2  wrong W range : ',Wmin,Wmax
        write(6,*)' '
        write(6,*)' '
        write(6,*)'************************************************* '

        STOP
C      ======

        Endif

C-----------------------------------------------------------------------
C
C Initialize for new energy - choose method
C      

      If(Ee.NE.Eold .OR. Wmin.NE.Wold1 .OR. Wmax.NE.Wold2)then

        write(6,*)'CompAZG2 initialized for  Ee = ',Ee
        write(6,*)'CompAZG2 requested W range : ',Wmin,Wmax

        xmax2=compaz(0.,Ee,22)
        Wmax2=2.*Ee*xmax2

        If(Wmin.GE.Wmax2)then
          write(6,*)'*********************************************** '
          write(6,*)' '
          write(6,*)' '
          write(6,*)'CompAZG2  requested W range above Wmax : ',Wmax2
          write(6,*)' '
          write(6,*)' '
          write(6,*)'*********************************************** '

          Eold=Ee
          Wold1=Wmin
          Wold2=Wmax
          Norm=0.

          E1=0.
          E2=0.
          P1=0.
          P2=0.

          Return
          EndIf
*
        If(Wmax.EQ.0. .OR. Wmax.GT. Wmax2 )then
           write(6,*)'CompAZG2    W range used   : ',Wmin,Wmax2
        Else
           Wmax2 = Wmax        
           EndIf
*
* Convert to 's_gg' 
*
        s0 = 4.*Ee*Ee
        smin=Wmin*Wmin
        smax=Wmax2*Wmax2
*
* Estimated 2-D space fraction
*

        If(smax.gt.smin)then
          frac2=smax+smax*log(2.*Ee*xmax2 / Wmax2)
          if(wmin.gt.0.)then
              ymax=log(2.*Ee*xmax2 / Wmin)
              frac2=frac2 - smin - smin*ymax
              EndIf
          frac2=frac2/(s0*xmax2*xmax2)
        else
          frac2=0.
          endif

        write(6,*)'Covered fraction of 2-D space : ',frac2

        If(frac2.EQ.0.)then
          Method=1
          write(6,*)'===>>>  Using 1-D generation  method'
        elseif(frac2.LT.0.1)then
          Method=2
          write(6,*)'===>>>  Using 2-D generation in (W,y)'
        else
          Method=3
          write(6,*)'===>>>  Using 2-D generation in (x_1,x_2)'
          call compazgn(-Ee,Norm,P1)
          Norm=Norm*Norm
          endif

        Eold=Ee
        Wold1=Wmin
        Wold2=Wmax

        Ngen=0
        Nevt=0

        EndIf
*
* Generate 
* ========
*
*  1-D
* -----
*
       If(Method.EQ.1)then
          call compazw1(Ein,Wmin,E1,P1,E2,P2)
          Return
          EndIf
*
*  2-D  (W,y)
* -----
*
       If(Method.EQ.2)then
          call compazw2(Ein,Wmin,Wmax2,E1,P1,E2,P2)
          Return
          EndIf
*
*  2-D   (x_1,x_2)
* -----
*
       if(Ein.LT.0.)Then

          If(Ngen.EQ.0)then
             frac=frac2*Norm
          else
             frac=Norm*Nevt/Ngen
             EndIf

          E1=frac
          P1=0.
          E2=0.
          P2=0.

          Return
          EndIf

       Nevt=Nevt+1

 999   continue

         Call CompazGN(Ee,E1,P1)
         Call CompazGN(Ee,E2,P2)
         Ngen=Ngen+1

         S=4.*E1*E2 
         If(S .LT. Smin .OR. S.GT.Smax)goto 999


       Return
       End
C
C==========================================================================
C
      Subroutine CompAZW1(Ein,Wgg,E1,P1,E2,P2)
C
C Generate photon energy from CompAZ parametrization
C
C  Input:  
C  ======
C             Ein   - electron beam energy [GeV]  50.<=Abs(Ein)<=500.
C             Wgg   - requested gamma-gamma cms energy
C
C  Output: 
C  =======
C
C  For Ein>0  E1,E2    - generated photons energies
C             P1,P1    - photon polarizations
C
C  For Ein<0  E1          -  1/L * dL/dW  [1/GeV] for E=ABS(Ein)
C             E2,P1,P2    - not used
C
C
      Implicit None
C
      Real Ein,Wgg,E1,P1,E2,P2
C
C Tables to store binned distribution
C      
      Integer Nbin,Nbin1
      Parameter (Nbin=100, Nbin1=Nbin+1)

      Real Ybin(Nbin1)
      
      Real Prob(Nbin)
      Real Yslope(Nbin)
    
      Real Polar1(Nbin)
      Real Polar2(Nbin)
      Real Pslope1(Nbin)
      Real Pslope2(Nbin)

      Integer Other(Nbin)

      Real Eold,Wold,Norm,x0
      Data Eold / 0. /
      Data Wold / 0. /

      Save Ybin,Prob,Yslope,Polar1,Pslope1,Polar2,Pslope2,Other
      Save Eold,Wold,Norm,x0

      Real Ee,z,z2,xmax,xmax2,Wmax,ymax,ymax2
      Integer i,j,n1,n2,imin,imax
      Real x,y1,y2,x1,x2,dy1,dy2,dy
      Real g1,g2,p11,p12,p21,p22,sum,ysign
      Real pmin,pmax,r,s
      Real rtab(3)

      Real Compaz
C
C Electron energy
C
      Ee=ABS(Ein)

      If(Ee.LT.50. .OR. Ee.GT.500.)then

        write(6,*)'************************************************* '
        write(6,*)' '
        write(6,*)' '
        write(6,*)'CompAZW1  energy out of range (50.-500.) : ',Ee
        write(6,*)' '
        write(6,*)' '
        write(6,*)'************************************************* '

        STOP
C      ======

        Endif
C
C Requested W
C
      If(Wgg .LE. 0.)then

        write(6,*)'************************************************* '
        write(6,*)' '
        write(6,*)' '
        write(6,*)'CompAZW1  requested W not positive : ',Wgg
        write(6,*)' '
        write(6,*)' '
        write(6,*)'************************************************* '

        STOP
C      ======

        Endif

C
C Initialize for new energy
C      

      If(Ee.NE.Eold .OR. Wgg.NE.Wold)then

        write(6,*)'CompAZW1 initialized for  Ee, Wgg = ',Ee,Wgg

        z=compaz(0.,Ee,20)

        xmax=compaz(0.,Ee,21)
        xmax2=compaz(0.,Ee,22)

        Wmax=2.*Ee*xmax2

        If(Wgg.GE.Wmax)then
          write(6,*)'*********************************************** '
          write(6,*)' '
          write(6,*)' '
          write(6,*)'CompAZW1  requested Wgg above Wmax : ',Wmax
          write(6,*)' '
          write(6,*)' '
          write(6,*)'*********************************************** '

          Eold=Ee
          Wold=Wgg
          Norm=0.

          E1=0.
          E2=0.
          P1=0.
          P2=0.

          Return
          EndIf
C        
C Rapidities
C
        ymax2=log(2*Ee*xmax2/Wgg)
C
        ymax=log(2*Ee*xmax/Wgg)
        If(ymax.LT.0.)ymax=abs(ymax)
        If(ymax.gT.ymax2)ymax=ymax2
C
        n1=(nbin-2)*ymax/ymax2 +1
        dy1=ymax/n1

        n2=nbin-n1
        dy2=(ymax2-ymax)/n2

        x0=Wgg/2./Ee

        y1=0.
        Ybin(1)=y1
        dy=dy1

        y2=y1+0.001*dy
        x1=x0*exp(y2)
        x2=x0*exp(-y2)

        g1=compaz(x1,Ee,0)*compaz(x2,Ee,0)
        p11=compaz(x1,Ee,10)
        p21=compaz(x2,Ee,10)

        sum=0.

        Do i=1,nbin

          y2=y1+dy
          Ybin(i+1)=y2

          if(i.eq.n1)y2=y2-0.001*dy

          x1=x0*exp(y2)
          x2=x0*exp(-y2)

          g2=compaz(x1,Ee,0)*compaz(x2,Ee,0)
          p12=compaz(x1,Ee,10)
          p22=compaz(x2,Ee,10)

          Prob(i)=(g1+g2)*dy/2.
          Yslope(i)=(g2-g1)/(g2+g1)

          sum=sum+(g1+g2)*dy/2.

          Polar1(i)=(p11+p12)/2.
          Pslope1(i)=(p12-p11)/(p12+p11)

          Polar2(i)=(p21+p22)/2.
          Pslope2(i)=(p22-p21)/(p22+p21)

          y1=y1+dy
          g1=g2
          p11=p12
          p21=p22

          if(i.eq.n1)then
            dy=dy2
            y2=y1+0.001*dy
            x1=x0*exp(y2)
            x2=x0*exp(-y2)

            g1=compaz(x1,Ee,0)*compaz(x2,Ee,0)
            p11=compaz(x1,Ee,10)
            p21=compaz(x2,Ee,10)
            endif

          EndDo
C
C Normalize
C
          Norm= Wgg/Ee/Ee * sum
          write(6,*)'CompAZW1:  1/L * dL/dW [1/GeV] = ',Norm

          sum=sum/nbin

          Do i=1,Nbin
            Prob(i)=Prob(i)/sum
            Other(i)=-1
            EndDo
C
C Sort table
C           
        Do i=1,Nbin
          pmin=2.
          pmax=0.
          Do j=1,Nbin
             if(Other(j).lt.0.and.Prob(j).lt.pmin)then
               imin=j
               pmin=Prob(j)
               EndIf
             if(Other(j).lt.0.and.Prob(j).ge.pmax)then
               imax=j
               pmax=Prob(j)
               EndIf
             EndDo

          other(imin)=imax
          prob(imax)=prob(imax)-(1.-prob(imin))

          EndDo

        Eold=Ee
        Wold=Wgg

        EndIf

       If(Ein.lt.0.)then
          E1=Norm
          E2=0.
          P1=0.
          P2=0.
          Return
          EndIf

       If(Norm.le.0.)then

          E1=0.
          E2=0.
          P1=0.
          P2=0.

          Return
          EndIf

       call Ranmar(rtab,3)

       if(rtab(1).GE.0.5)then
         i=2*nbin*(rtab(1)-0.5) + 1
         ysign=1.
       else
         i=2*nbin*rtab(1) + 1
         ysign=-1.
         endif

       if(i.gt.nbin)i=nbin

       if(rtab(2).gt.prob(i))i=other(i)

       r=rtab(3)
       s=Yslope(i)

       if(s.eq.0)then
          x=r
       Elseif(s.gt.0)then
          x=(1-s)*(1-s)+4*s*r
          x=(s-1+sqrt(x))/2./s
       Else
          x=(1-s)*(1-s)+4*s*r
          x=(s-1+sqrt(x))/2./s
          EndIf

       y1=Ybin(i)+x*(Ybin(i+1)-Ybin(i))

       If(ysign.GT.0.)then
          E1=Ee*x0*exp(y1)
          P1=Polar1(i)*(1+(2*x-1.)*Pslope1(i))
          E2=Ee*x0*exp(-y1)
          P2=Polar2(i)*(1+(2*x-1.)*Pslope2(i))
       Else
          E2=Ee*x0*exp(y1)
          P2=Polar1(i)*(1+(2*x-1.)*Pslope1(i))
          E1=Ee*x0*exp(-y1)
          P1=Polar2(i)*(1+(2*x-1.)*Pslope2(i))
          EndIf

       Return
       End
*
*=================================================================
*
      Subroutine CompAZW2(Ein,Wmin,Wmax,E1,P1,E2,P2)
C
C Generate photon energy from CompAZ parametrization
C
C  Input:  
C  ======
C             Ein   - electron beam energy [GeV]  50.<=Abs(Ein)<=500.
C             Wmin  - minimum W
C             Wmax  - maximum W [GeV]
C
C  Output: 
C  =======
C
C  For Ein>0  E1,E2    - generated photons energies
C             P1,P1    - photon polarizations
C
C  For Ein<0  E1          -  1/L * Int(dL/dW) from Wmin to Wmax for ABS(Ein)
C             E2,P1,P2    - not used
C
C
      Implicit None
C
      Real Ein,Wmin,Wmax,E1,P1,E2,P2
C
C Tables to store binned distribution
C      
      Integer NbinS,NbinS1,NminS
      Parameter (NbinS=50, NbinS1=NbinS+1, NminS=5)

      Integer NbinY,NbinY1,NminY
      Parameter (NbinY=200, NbinY1=NbinY+1,NminY=20)
      
      Integer Nbin
      Parameter (Nbin=NbinS*NbinY)

      Real Sbin(NbinS1)
      Real Ybin(NbinY1,NbinS1)
      
      Real Prob(Nbin)
      Real Yslope(Nbin)
      Real Sslope(Nbin)
    
      Real Polar1(Nbin)
      Real Polar2(Nbin)
      Real PslopeY1(Nbin)
      Real PslopeS1(Nbin)
      Real PslopeY2(Nbin)
      Real PslopeS2(Nbin)

      Integer Other(Nbin)

      Real Eold,Wold1,Wold2,Norm
      Data Eold / 0. /
      Data Wold1,Wold2 / 0., 0. /

      Save Ybin,Sbin,Prob,Yslope,Sslope,Other
      Save Polar1,PslopeS1,PslopeY1,Polar2,PslopeS2,PslopeY2
      Save Eold,Wold1,Wold2,Norm

      Real Ee,z,xmax,xmax2,Wmax2,Smin,Smax,Smax2,ymax,ymax2
      Integer i,j,is,iy,ny1,ny2,ns1,ns2,imin,imax
      Real S,W,y1,y2,x1,x2,dy1,dy2,ds1,ds2,dy,ds,dynew,dsnew
      Real x0,g2,pl1,pl2,sum,ysign
      Real pmin,pmax,r,sls,sly,shy,shs
      Real rtab(4)

      Real Compaz
C
C Electron energy
C
      Ee=ABS(Ein)

      If(Ee.LT.50. .OR. Ee.GT.500.)then

        write(6,*)'************************************************* '
        write(6,*)' '
        write(6,*)' '
        write(6,*)'CompAZW2  energy out of range (50.-500.) : ',Ee
        write(6,*)' '
        write(6,*)' '
        write(6,*)'************************************************* '

        STOP
C      ======

        Endif
C
C Requested W
C
      If(Wmin .LE. 0. .OR. 
     +   (Wmax.NE.0. .AND. Wmin.GT.Wmax) )then

        write(6,*)'************************************************* '
        write(6,*)' '
        write(6,*)' '
        write(6,*)'CompAZW2  invalid W range : ',Wmin,Wmax
        write(6,*)' '
        write(6,*)' '
        write(6,*)'************************************************* '

        STOP
C      ======

        Endif

C
C Initialize for new energy
C      

      If(Ee.NE.Eold .OR. Wmin.NE.Wold1 .OR. Wmax.NE.Wold2)then

        write(6,*)'CompAZW2 initialized for  Ee = ',Ee
        write(6,*)'CompAZW2 requested W range : ',Wmin,Wmax

        z=compaz(0.,Ee,20)

        xmax=compaz(0.,Ee,21)
        xmax2=compaz(0.,Ee,22)

        Wmax2=2.*Ee*xmax2

        If(Wmin.GE.Wmax2)then
          write(6,*)'*********************************************** '
          write(6,*)' '
          write(6,*)' '
          write(6,*)'CompAZW2  requested W range above Wmax : ',Wmax2
          write(6,*)' '
          write(6,*)' '
          write(6,*)'*********************************************** '

          Eold=Ee
          Wold1=Wmin
          Wold2=Wmax
          Norm=0.

          E1=0.
          E2=0.
          P1=0.
          P2=0.

          Return
          EndIf
C        
        If(Wmax.GT.Wmax2 .OR. Wmax.EQ.0)then
           write(6,*)'CompAZW2    W range used   : ',Wmin,Wmax2
        else
           Wmax2=Wmax
           EndIf
C
C Rapidities for  Wmin
C
        ymax2=log(2*Ee*xmax2/Wmin)
C
        ymax=log(2*Ee*xmax/Wmin)
        If(ymax.LT.0.)ymax=abs(ymax)
        If(ymax.gT.ymax2)ymax=ymax2
C
C y binning
C
        ny1=(NbinY-2*NminY)*ymax/ymax2 + NminY
        ny2=NbinY - ny1
C
C S binning
C
        Smin=Wmin*Wmin
        Smax2=Wmax2*Wmax2
        Smax = 4.*Ee*Ee*xmax*xmax

        If(Smax .GT. Smax2)Smax=Smax2
        If(Smax .LT. Smin)Smax=Smin

        ns1=(NbinS-2*NminS)*(Smax-Smin)/(Smax2-Smin) + NminS
        ns2=NbinS - ns1

        ds1=(Smax-Smin)/ns1
        ds2=(Smax2-Smax)/ns2
C
C Clear tables
C
        Do i=1,Nbin
           Prob(i)=0.
           Yslope(i)=0.
           Sslope(i)=0.
           Polar1(i)=0.
           Polar2(i)=0.
           PslopeY1(i)=0.
           PslopeY2(i)=0.
           PslopeS1(i)=0.
           PslopeS2(i)=0.
           EndDo
C
C Loop over S values
C
        S=Smin-ds1
        ds=ds1
        dsnew=ds

        sum=0.

        do is=0,NbinS

          if(is.eq.ns1)dsnew=ds2

          S=S+ds

          Sbin(is+1)=S
          W=sqrt(S)

          ymax2=log(2*Ee*xmax2/W)

          ymax=log(2*Ee*xmax/W)
          If(ymax.LT.0.)ymax=abs(ymax)
          If(ymax.gT.ymax2)ymax=ymax2

          dy1=ymax/ny1
          dy2=(ymax2-ymax)/ny2

          y1= -dy1
          dy=dy1
          dynew=dy

          Do iy=0,nbiny

            if(iy.eq.ny1)dynew=dy2

            y2=y1+dy
            Ybin(iy+1,is+1)=y2

            x0=W/2./Ee
C
C Low-s low-y bin
C
            If(is.eq.ns1)x0=sqrt(S-0.001*ds)/2./Ee

            if(iy.eq.0)y2=y1+dy+0.001*dynew

            if(iy.eq.ny1)y2=y1+dy-0.001*dy

            x1=x0*exp(y2)
            x2=x0*exp(-y2)

            g2=compaz(x1,Ee,0)*compaz(x2,Ee,0)

            pl1=compaz(x1,Ee,10)
            pl2=compaz(x2,Ee,10)

            i=(is-1)*NbinY+iy

            If(iy.GT.0 .and. is.gt.0 )then
              Prob(i)=Prob(i)+g2*dy*ds/4.
              Yslope(i)=Yslope(i)+g2*dy*ds/4.
              Sslope(i)=Sslope(i)+g2*dy*ds/4.
              Polar1(i)=Polar1(i)+pl1/4.
              Polar2(i)=Polar2(i)+pl2/4.
              PslopeY1(i)=PslopeY1(i)+pl1/4.
              PslopeS1(i)=PslopeS1(i)+pl1/4.
              PslopeY2(i)=PslopeY2(i)+pl2/4.
              PslopeS2(i)=PslopeS2(i)+pl2/4.
              endif
C
C Low-s high-y bin
C
            if(iy.eq.ny1)then
               y2=y1+dy+0.001*dynew

               x1=x0*exp(y2)
               x2=x0*exp(-y2)

               g2=compaz(x1,Ee,0)*compaz(x2,Ee,0)
               pl1=compaz(x1,Ee,10)
               pl2=compaz(x2,Ee,10)
               EndIf

            i=(is-1)*NbinY+iy+1

            If(iy.LT.NbinY .and. is.gt.0 )then
              Prob(i)=Prob(i)+g2*dynew*ds/4.
              Yslope(i)=Yslope(i)-g2*dynew*ds/4.
              Sslope(i)=Sslope(i)+g2*dynew*ds/4.
              Polar1(i)=Polar1(i)+pl1/4.
              Polar2(i)=Polar2(i)+pl2/4.
              PslopeY1(i)=PslopeY1(i)-pl1/4.
              PslopeS1(i)=PslopeS1(i)+pl1/4.
              PslopeY2(i)=PslopeY2(i)-pl2/4.
              PslopeS2(i)=PslopeS2(i)+pl2/4.
              endif
C
C High-s low-y bin
C
            If(is.eq.ns1)x0=sqrt(S+0.001*dsnew)/2./Ee

            if(iy.eq.0)y2=y1+dy+0.001*dy
            if(iy.eq.ny1)y2=y1+dy-0.001*dy

            If(is.eq.ns1 .OR. iy.eq.ny1)then
              x1=x0*exp(y2)
              x2=x0*exp(-y2)

              g2=compaz(x1,Ee,0)*compaz(x2,Ee,0)
              pl1=compaz(x1,Ee,10)
              pl2=compaz(x2,Ee,10)
              EndIf

            i=is*NbinY+iy

            If(iy.GT.0 .and. is.lt.NbinS )then
              Prob(i)=Prob(i)+g2*dy*dsnew/4.
              Yslope(i)=Yslope(i)+g2*dy*dsnew/4.
              Sslope(i)=Sslope(i)-g2*dy*dsnew/4.
              Polar1(i)=Polar1(i)+pl1/4.
              Polar2(i)=Polar2(i)+pl2/4.
              PslopeY1(i)=PslopeY1(i)+pl1/4.
              PslopeS1(i)=PslopeS1(i)-pl1/4.
              PslopeY2(i)=PslopeY2(i)+pl2/4.
              PslopeS2(i)=PslopeS2(i)-pl2/4.
              endif

C
C High-s high-y bin
C
            if(iy.eq.ny1)then
               y2=y1+dy+0.001*dynew

               x1=x0*exp(y2)
               x2=x0*exp(-y2)

               g2=compaz(x1,Ee,0)*compaz(x2,Ee,0)
               pl1=compaz(x1,Ee,10)
               pl2=compaz(x2,Ee,10)
               EndIf

            i=is*NbinY+iy+1

            If(iy.LT.NbinY .and. is.lt.NbinS )then
              Prob(i)=Prob(i)+g2*dynew*dsnew/4.
              Yslope(i)=Yslope(i)-g2*dynew*dsnew/4.
              Sslope(i)=Sslope(i)-g2*dynew*dsnew/4.
              Polar1(i)=Polar1(i)+pl1/4.
              Polar2(i)=Polar2(i)+pl2/4.
              PslopeY1(i)=PslopeY1(i)-pl1/4.
              PslopeS1(i)=PslopeS1(i)-pl1/4.
              PslopeY2(i)=PslopeY2(i)-pl2/4.
              PslopeS2(i)=PslopeS2(i)-pl2/4.
              endif
            
            y1=y1+dy
            dy=dynew
            EndDo

          ds=dsnew
          EndDo
C
C Normalize
C
          sum=0.

          Do i=1,Nbin

            If(Prob(i).GT.0.)Then

              sum=sum+Prob(i)

              Yslope(i)=Yslope(i)/Prob(i)
              Sslope(i)=Sslope(i)/Prob(i)
              PslopeY1(i)=PslopeY1(i)/Polar1(i)
              PslopeS1(i)=PslopeS1(i)/Polar1(i)
              PslopeY2(i)=PslopeY2(i)/Polar2(i)
              PslopeS2(i)=PslopeS2(i)/Polar2(i)

            Else

              Yslope(i)=0.
              Sslope(i)=0.

              EndIf

             EndDo
C
          Norm= sum/Ee/Ee/2.
          write(6,*)'CompAZW2:  1/L * Int(dL/dW)  = ',Norm

          sum=sum/nbin

          Do i=1,Nbin
            Prob(i)=Prob(i)/sum
            Other(i)=-1
            EndDo
C
C Sort table
C           
        Do i=1,Nbin
          pmin=2.
          pmax=0.
          Do j=1,Nbin
             if(Other(j).lt.0.and.Prob(j).lt.pmin)then
               imin=j
               pmin=Prob(j)
               EndIf
             if(Other(j).lt.0.and.Prob(j).ge.pmax)then
               imax=j
               pmax=Prob(j)
               EndIf
             EndDo

          other(imin)=imax
          prob(imax)=prob(imax)-(1.-prob(imin))

          EndDo

        Eold=Ee
        Wold1=Wmin
        Wold2=Wmax

        EndIf

       If(Ein.lt.0.)then
          E1=Norm
          E2=0.
          P1=0.
          P2=0.
          Return
          EndIf

       If(Norm.le.0.)then

          E1=0.
          E2=0.
          P1=0.
          P2=0.

          Return
          EndIf

       call Ranmar(rtab,4)

       if(rtab(1).GE.0.5)then
         i=2*nbin*(rtab(1)-0.5) + 1
         ysign=1.
       else
         i=2*nbin*rtab(1) + 1
         ysign=-1.
         endif

       if(i.gt.nbin)i=nbin

       if(rtab(2).gt.prob(i))i=other(i)

       r=rtab(3)
       sly=Yslope(i)

       if(sly.eq.0)then
          shy=r
       Elseif(sly.gt.0)then
          shy=(1-sly)*(1-sly)+4*sly*r
          shy=(sly-1+sqrt(shy))/2./sly
       Else
          shy=(1-sly)*(1-sly)+4*sly*r
          shy=(sly-1+sqrt(shy))/2./sly
          EndIf

       r=rtab(4)
       sls=Sslope(i)

       if(sls.eq.0)then
          shs=r
       Else
          shs=(1-sls)*(1-sls)+4*sls*r
          shs=(sls-1+sqrt(shs))/2./sls
          EndIf

       iy=mod(i-1,nbiny) + 1
       is=(i-1)/NbinY + 1

       S=Sbin(is)+shs*(Sbin(is+1)-Sbin(is))
       W=Sqrt(S)
       x0=W/2./Ee

       y1=Ybin(iy,is)+shy*(Ybin(iy+1,is)-Ybin(iy,is))
       y2=Ybin(iy,is+1)+shy*(Ybin(iy+1,is+1)-Ybin(iy,is+1))
       y1=shs*y2 + (1.-shs)*y1

       PL1=Polar1(i)*(1+(2*shy-1.)*PslopeY1(i)+(2*shs-1.)*PslopeS1(i))
       PL2=Polar2(i)*(1+(2*shy-1.)*PslopeY2(i)+(2*shs-1.)*PslopeS2(i))

       If(ysign.GT.0.)then
          E1=Ee*x0*exp(y1)
          P1=PL1
          E2=Ee*x0*exp(-y1)
          P2=PL2
       Else
          E2=Ee*x0*exp(y1)
          P2=PL1
          E1=Ee*x0*exp(-y1)
          P1=PL2
          EndIf

       Return
       End
