~frederix/+junk/matrix

module parameters
! User input:
  ! next     = number of external gluons (initial+final)
  ! npoints  = number of PS points to compute (negative to keep going)
  ! npermtry = number of colour flow permutations to include per PS
  !            point (needs to be greater than 2 and smaller than the
  !            maximal, (next-1)!)
  ! alphas   = the value of the strong coupling
  ! sqrtshat = collision energy
  integer, parameter   :: next=5
  integer, parameter   :: npermtry=24
  integer*8, parameter :: npoints=10000
  real*8, parameter    :: alphas=0.12d0
  real*8, parameter    :: sqrtshat=1000d0
! constants:  
  integer, parameter :: ninc=2,nfin=next-ninc
  real*8, parameter :: pi=3.14159265358979323846d0
! External momenta (generated by 'call get_momenta()'): 
  real*8, dimension(0:3,next) :: p
! helicity (polarisation) and permutation used for amplitude
! computation:
  integer, dimension(next) :: ip
! final or initial state particle:
  integer, dimension(next) :: ifinal
  data ifinal(1:next) / -1,-1,nfin*1 /
  complex*16, dimension(npermtry,0:maskr(next)) :: amp
end module parameters

program matrix
  use parameters
  implicit none
  integer :: ih,iperm,jperm,i,ib1,ib2,ib,ih1,ih2,towrite=1
  integer*8 :: iden,ipoint=0,inonzero=0,nperm
  integer, dimension(next,npermtry) :: ips
  real*8, dimension(next) :: pmass=0d0
  real*8 :: amp2,integral=0d0,uncertainty=0d0,flux,matrix_element,wgt, &
       perm_vol
  real*8, parameter :: conv=389379660d0
  ! number of colour permutations
  nperm=factorial(next-1)
  ! colour, polarisation incoming gluons: 8*8, 2*2
  ! identical final state particle factor: nfin!
  iden=8*8*2*2 * factorial(nfin)
  ! flux factor
  flux=1d0/(2d0*sqrtshat**2)
  do
     if (npoints.gt.0 .and. inonzero.eq.npoints) exit
     ipoint=ipoint+1
     call get_momenta(sqrtshat,pmass,wgt)
     do i=1,2
        p(0,i)=-p(0,i) ! treat all momenta as outgoing
     enddo
     ! weight from Ramdo is logarithmic and does not
     ! include 1/(2pi)^(3n-4). Compensate...     
     wgt=exp(wgt)/((2*pi)**(3*nfin-4))

     ! Check if PS point passes generation cuts
     if (fail_cuts()) cycle
     inonzero=inonzero+1
     amp2=0d0
     ! Determine which npermtry colour ordered amplitudes to
     ! compute. Do this randomly.
     do iperm=1,npermtry
        if (npermtry.eq.nperm) then
           ! Summing over all permutations. 
           call set_iperm(iperm,ips(1,iperm))
        elseif (npermtry.lt.nperm) then
           if (iperm.eq.1) then
              call get_iperm_rn(ips(1,iperm))
           else
111           call get_iperm_rn(ips(1,iperm))
              ! second (or third, etc.) should always be different
              ! from previous ones. Simply regenerate. This is not
              ! optimal, but sufficient if npermtry<<nperm, which
              ! should be the case for large number of gluons
              do i=1,iperm-1
                 if (all(ips(1:next,iperm).eq.ips(1:next,i))) goto 111
              enddo
           endif
        else
           write (*,*) 'npermtry should be smaller than',nperm
           stop 1
        endif
     enddo
     do iperm=1,npermtry
        ip(1:next)=ips(1:next,iperm)
        ! Compute the matrix colour ordered amplitude, first computing
        ! the wavefunctions correspoding to the sum of the first
        ! next-1 particles, dotting that into the gluon polarisation
        ! for the next particle.
        call compute_amplitude(iperm)
     enddo

!!$     do iperm=1,1
!!$        write (*,*) ips(1:next,iperm) ,'=>',iperm
!!$        do i=0,maskr(next)
!!$           write (*,*) i,amp(iperm,i)
!!$        enddo
!!$     enddo

     ! Loop over all colour ordered amplitudes computed, and multiply
     ! by the colour factor to get the matrix element squared. If
     ! npermtry<nperm, also include the weight from the MC-ing over
     ! the colour-ordered amplitudes. The colour factor for the
     ! off-diagonal contributions includes a factor 2, because only
     ! upper triangle of colour matrix is included.
     do ih=0,maskr(next)
        do iperm=1,npermtry
           ! do some bit moving to align the helicity with the
           ! permutation iperm
           ih1=0
           do ib=1,next
              call mvbits(ih,ips(ib,iperm)-1,1,ih1,ib-1)
           enddo
           do jperm=iperm,npermtry
              ! do some bit moving to align the helicity with the
              ! permutation jperm
              ih2=0
              do ib=1,next
                 call mvbits(ih,ips(ib,jperm)-1,1,ih2,ib-1)
              enddo
              ! Compute the weight from the MC-ing over permutations
              if (iperm.eq.jperm) then
                 perm_vol=dble(nperm)/dble(npermtry)
              else
                 perm_vol=dble(nperm*(nperm-1))/dble(npermtry*(npermtry-1))
              endif
              ! Compute the amplitude squared
              amp2=amp2+dble(amp(iperm,ih1)*dconjg(amp(jperm,ih2)))* &
                   color_factor_list(ips(1,iperm),ips(1,jperm))*perm_vol
           enddo
        enddo
     enddo
     ! Include the strong coupling and the identical particle factor.
     amp2=amp2*(4*pi*alphas)**nfin/dble(iden)
     ! Include the flux factor, the phase-space weight, and the
     ! conversion from GeV->pb.
     matrix_element=amp2*flux*wgt*conv

     ! Add the current matrix element to the total cross section
     integral=integral+matrix_element
     uncertainty=uncertainty+matrix_element**2
     if (inonzero.gt.10*towrite .and. towrite.lt.10000) towrite=towrite*10
     if (mod(inonzero,towrite).eq.0) then
        write (88,*) ipoint,inonzero,matrix_element,integral/dble(ipoint), &
             sqrt(abs(uncertainty/dble(ipoint)-(integral/dble(ipoint))**2)/dble(ipoint)), &
             100d0*sqrt(abs(uncertainty/dble(ipoint)-(integral/dble(ipoint))**2)/dble(ipoint))/(integral/dble(ipoint)) 
        write (*,*) ipoint,inonzero,matrix_element,integral/dble(ipoint), &
             sqrt(abs(uncertainty/dble(ipoint)-(integral/dble(ipoint))**2)/dble(ipoint)), &
             100d0*sqrt(abs(uncertainty/dble(ipoint)-(integral/dble(ipoint))**2)/dble(ipoint))/(integral/dble(ipoint)) 
        call flush(88)
     endif
  enddo
  ! Print the final result to the screen
  integral=integral/dble(npoints)
  uncertainty=uncertainty/dble(npoints)
  uncertainty=sqrt(abs(uncertainty-integral**2)/dble(npoints))
  if (npoints.gt.10) write (*,*) 'sigma:',integral,'+/-',uncertainty,'(',100d0*uncertainty/integral,'%)'

contains
  logical function fail_cuts()
    ! Cuts on the phase-space point.
    use parameters
    implicit none
    integer i,j
    fail_cuts=.false.
    do i=3,next
       if (pt(p(0,i)).lt.60d0) then
          fail_cuts=.true.
          return
       endif
       if (abs(eta(p(0,i))).gt.2d0) then
          fail_cuts=.true.
          return
       endif
       if (i.ne.next) then
          do j=i+1,next
             if (DeltaR(p(0,i),p(0,j)).lt.0.7d0) then
                fail_cuts=.true.
                return
             endif
          enddo
       endif
    enddo
  end function fail_cuts
  
  real*8 function pt(p)
    ! transverse momentum of 'p'
    implicit none
    real*8, dimension(0:3) :: p
    pt=sqrt(p(1)**2+p(2)**2)
  end function pt
  
  real*8 function eta(p)
    ! pseudo-rapidity of 'p'
    implicit none
    real*8, dimension(0:3) :: p
    real*8 :: theta
    theta=acos(p(3)/sqrt(p(1)**2+p(2)**2+p(3)**2))
    eta=-log(dtan(theta/2d0))
  end function eta

  real*8 function delta_phi(p1,p2)
    ! azimuthal difference of 'p1' and 'p2'
    implicit none
    real*8, dimension(0:3) :: p1,p2
    real*8 :: denom
    denom=pt(p1)*pt(p2)
    delta_phi=acos((p1(1)*p2(1)+p1(2)*p2(2))/denom)
  end function delta_phi

  real*8 function deltaR(p1,p2)
    ! Distance (Delta-R) between 'p1' and 'p2'
    implicit none
    real*8, dimension(0:3) :: p1,p2
    deltaR=sqrt(delta_phi(p1,p2)**2+(eta(p1)-eta(p2))**2)
  end function deltaR

  real*8 function color_factor(iperm,jperm)
    ! Compute colour factor for permutation numbers 'iperm' and
    ! 'jperm'
    use parameters
    implicit none
    integer :: iperm,jperm
    integer, dimension(next) :: iper,jper,ips
    call set_iperm(iperm,ips)
    iper(1:next)=ips(:)
    call set_iperm(jperm,ips)
    jper(1:next)=ips(:)
    color_factor=color_factor_list(iper,jper)
  end function color_factor

  real*8 function color_factor_list(iper,jper)
    ! Given two permutations (iper and jper) compute corresponding
    ! colour factor. If the two configurations are identical, include
    ! an additional factor 2.
    use parameters
    implicit none
    integer :: i,j,k,id0,jd0,id1,jd1,colfac,nperms
    integer, dimension(next) :: iper,jper
    integer, dimension(0:1,next) :: idel,jdel
    logical, dimension(next) :: lper
    ! add the factor 2 if need be
    if (all(iper.eq.jper)) then
       colfac=1
    else
       colfac=2
    endif
    ! Determine the delta's to contract
    do i=1,next
       ! delta's of the amplitude
       idel(0,i)=iper(i)
       jdel(0,i)=iper(mod(i,next)+1)
       ! delta's of the conjugate amplitude
       idel(1,i)=jper(i)
       jdel(1,i)=jper(mod(i,next)+1)
    enddo
    ! No closed string of delta's found so far: set all the labels to
    ! false
    lper(1:next)=.false.
    ! Loop over the delta's.
    do k=1,next
       ! If this delta already belonged to a previously closed
       ! string. Skip it
       if (lper(k)) cycle
       ! The two labels of the delta in the amplitude
       id0=idel(0,k)
       jd0=jdel(0,k)
       ! Go searching to see to which delta's it is connected
       do
          do j=1,next
             if (jdel(1,j).eq.jd0) then
                ! Found the delta in the conjugate amplitude to which
                ! the delta in the amplitude is connected to
                id1=idel(1,j)
                exit
             endif
          enddo
          if (id0.eq.id1) then
             ! The delta in the conjugate amplitude, links back to the
             ! starting delta. Hence, one closed string found: add a
             ! factor 3 to the colour factor and start with the next
             ! delta in the amplitude.
             colfac=colfac*3
             exit
          endif
          do i=1,next
             if (idel(0,i).eq.id1) then
                ! Found delta in the amplitude to which the delta in
                ! the conjugate amplitude connects to. Hence, it's
                ! part of a closed string, so set corresponding lper
                ! to true.
                jd0=jdel(0,i)
                lper(i)=.true.
                exit
             endif
          enddo
       enddo
    enddo
    color_factor_list=dble(colfac)
  end function color_factor_list
  
  integer*8 function factorial(ifact)
    ! computes the factorial of 'ifact'. Uses long integers to deal
    ! with large numbers
    implicit none
    integer, value :: ifact
    factorial=1
    do while (ifact.gt.1)
       factorial=factorial*ifact
       ifact=ifact-1
    enddo
  end function factorial
end program matrix

subroutine set_iperm(iperm,ipss)
  ! Given permutation number 'iperm', fills the corresponding
  ! permutation. This is slow when the number of permutations is
  ! large, and easily becomes the most time-consuming part of the
  ! calculation if used too often.
  use parameters
  implicit none
  integer :: iperm,i
  integer, dimension(next-1) :: ips
  integer, dimension(next) :: ipss
  integer :: iflag
  do i=1,next-1
     ips(i)=i
  enddo
  do i=1,iperm-1
     call ipnext(ips,next-1,iflag)
  enddo
  ipss(1:next-1)=ips(:)
  ipss(next)=next
end subroutine set_iperm

subroutine get_iperm_rn(ips)
  ! Computes a random permutation 'ips'.
  use parameters
  implicit none
  integer :: i,j,k,itemp
  integer,dimension(1:next-1) :: ileft
  integer,dimension(next) :: ips
  real*8, external :: rn
  ! Create a pool of numbers to pick from
  do i=1,next
     ileft(i)=i
  enddo
  ! Randomly remove a number from the pool 'ileft' and add it to the
  ! 'ips' permutation list.
  do i=1,next-1
     k=0
     itemp=int(rn(1)*(next-i))+1
     do j=1,itemp
        k=k+1
        do while (ileft(k).eq.0)
           k=k+1
        enddo
     enddo
     ips(i)=ileft(k)
     ileft(k)=0
  enddo
  ! Last number in permutation is always 'next': this removes all cyclic
  ! permutations.
  ips(next)=next
end subroutine get_iperm_rn

subroutine compute_amplitude(iperm)
  use parameters
  implicit none
  integer,parameter :: mw=maskr(next-1),mm=maskr(next-1)+1,zero=0
  integer :: imom,i,ih,level,j,isplit,jsplit,iwf,nnext,iperm, &
       a3,b3,a4,b4,c4,ia,ib,ic,ext,nh,icount
  complex*16, dimension(4,mm,mm) :: wf
  integer,dimension(mm) :: nwf
  real*8, dimension(0:3,mm) :: pp
  complex*16 :: propagator,contribution
  complex*16,dimension(4) :: wfout
  complex*16,parameter :: ci=(0d0,1d0)
  complex*16,dimension(4,0:1) :: wffinal

  nwf=0
  nnext=next-1
  do i=1,next
     ext=0
     ext=ibset(ext,i-1)
     pp(0:3,ext)=p(0:3,ip(i))
     do ih=0,1
        nwf(ext)=nwf(ext)+1
        call v_ext(pp(0,ext),ih,1,wf(1,nwf(ext),ext))
     enddo
  enddo

  do level=1,nnext-1
     do i=1,nnext-level
        ext=0
        nh=1
        do j=0,level
           ext=ibset(ext,i+j-1)
           nh=nh*2
        enddo
        nwf(ext)=nh
        wf(1:4,1:nwf(ext),ext)=(0d0,0d0)

        pp(0:3,ext)=0d0
        do j=1,nnext
           if (btest(ext,j-1)) pp(0:3,ext)=pp(0:3,ext)+pp(0:3,ibset(zero,j-1))
        enddo

        do isplit=1+trailz(ext),popcnt(ext)-1+trailz(ext)
           a3=merge_bits(ext,zero,maskr(isplit))
           b3=merge_bits(zero,ext,maskr(isplit))
           do ib=1,nwf(b3)
              do ia=1,nwf(a3)
                 call gluon3(wf(1,ia,a3),pp(0,a3), &
                             wf(1,ib,b3),pp(0,b3), &
                             wfout)
                 icount=(ib-1)*nwf(a3)+ia
                 wf(1:4,icount,ext)=wf(1:4,icount,ext)+wfout(1:4)
              enddo
           enddo
           do jsplit=isplit+1,popcnt(ext)-1+trailz(ext)
              a4=a3
              b4=merge_bits(b3,zero,maskr(jsplit))
              c4=merge_bits(zero,b3,maskr(jsplit))
              do ic=1,nwf(c4)
                 do ib=1,nwf(b4)
                    do ia=1,nwf(a4)
                       call gluon4(wf(1,ia,a4),wf(1,ib,b4),wf(1,ic,c4),wfout)
                       icount=(ic-1)*nwf(a4)*nwf(b4)+(ib-1)*nwf(a4)+ia
                       wf(1:4,icount,ext)=wf(1:4,icount,ext)+wfout(1:4)
                    enddo
                 enddo
              enddo
           enddo
        enddo
        if (level.ne.nnext-1) then
           propagator=-ci/(pp(0,ext)**2-pp(1,ext)**2-pp(2,ext)**2-pp(3,ext)**2)
           wf(1:4,1:nwf(ext),ext)=wf(1:4,1:nwf(ext),ext)*propagator
        endif
     enddo
  enddo

  amp(iperm,:)=0d0
  do iwf=1,nwf(ext)
     do i=1,2
        contribution= &
                    (wf(1,iwf,ext)*wf(1,i,ibset(zero,nnext))- &
                     wf(2,iwf,ext)*wf(2,i,ibset(zero,nnext))- &
                     wf(3,iwf,ext)*wf(3,i,ibset(zero,nnext))- &
                     wf(4,iwf,ext)*wf(4,i,ibset(zero,nnext)))
        ih=iwf-1
        if (i.eq.2) ih=ibset(ih,next-1)
        amp(iperm,ih)=amp(iperm,ih)+contribution
     enddo
  enddo

end subroutine compute_amplitude
   
subroutine v_ext(p,ihel,ifinal,wf)
  ! External gluon wavefunction. From HELAS.
  implicit none
  integer :: ihel,ifinal
  real*8, dimension(0:3) :: p
  complex*16, dimension(4) :: wf
  real*8, parameter :: rzero=0d0,rhalf=0.5d0,sqh=sqrt(0.5d0)
  real*8 :: hel,pt2,pp,pt,pzpt
  hel = dble(2*ihel-1)
  pt2 = p(1)**2+p(2)**2
  pp = min(p(0),sqrt(pt2+p(3)**2))
  pt = min(pp,sqrt(pt2))
  pp = p(0)
  pt = sqrt(p(1)**2+p(2)**2)
  wf(1) = dcmplx( rZero )
  wf(4) = dcmplx( hel*pt/pp*sqh )
  if ( pt.ne.rZero ) then
     pzpt = p(3)/(pp*pt)*sqh*hel
     wf(2) = dcmplx( -p(1)*pzpt , -ifinal*p(2)/pt*sqh )
     wf(3) = dcmplx( -p(2)*pzpt ,  ifinal*p(1)/pt*sqh )
  else
     wf(2) = dcmplx( -hel*sqh )
     wf(3) = dcmplx( rZero , ifinal*sign(sqh,p(3)) )
  endif
end subroutine v_ext

subroutine gluon3(wf1,pwf1,wf2,pwf2,wf)
  ! Colour-ordered three gluon interaction
  implicit none
  complex*16,dimension(4) :: wf1,wf2,wf
  real*8,dimension(0:3) :: pwf1,pwf2
  complex*16, parameter :: prefact=(0d0,1d0)/sqrt(2d0)
  complex*16 :: TMP1,TMP2,TMP3
  TMP1 = (wf1(1)*wf2(1)-wf1(2)*wf2(2)-wf1(3)*wf2(3)-wf1(4)*wf2(4))
  TMP2 = (wf1(1)*pwf2(0)-wf1(2)*pwf2(1)-wf1(3)*pwf2(2)-wf1(4)*pwf2(3))
  TMP3 = (wf2(1)*pwf1(0)-wf2(2)*pwf1(1)-wf2(3)*pwf1(2)-wf2(4)*pwf1(3))
  wf(1) = prefact*(TMP1*(pwf1(0)-pwf2(0))+2d0*TMP2*wf2(1)-2d0*TMP3*wf1(1))
  wf(2) = prefact*(TMP1*(pwf1(1)-pwf2(1))+2d0*TMP2*wf2(2)-2d0*TMP3*wf1(2))
  wf(3) = prefact*(TMP1*(pwf1(2)-pwf2(2))+2d0*TMP2*wf2(3)-2d0*TMP3*wf1(3))
  wf(4) = prefact*(TMP1*(pwf1(3)-pwf2(3))+2d0*TMP2*wf2(4)-2d0*TMP3*wf1(4))
end subroutine gluon3

subroutine gluon4(wf1,wf2,wf3,wf)
  ! Colour-ordered four gluon interaction
  implicit none
  complex*16,dimension(4) :: wf1,wf2,wf3,wf
  complex*16, parameter :: prefact=(0d0,0.5d0)
  complex*16 :: TMP1,TMP2,TMP3
  TMP1 = (wf1(1)*wf2(1)-wf1(2)*wf2(2)-wf1(3)*wf2(3)-wf1(4)*wf2(4))
  TMP2 = (wf1(1)*wf3(1)-wf1(2)*wf3(2)-wf1(3)*wf3(3)-wf1(4)*wf3(4))
  TMP3 = (wf2(1)*wf3(1)-wf2(2)*wf3(2)-wf2(3)*wf3(3)-wf2(4)*wf3(4))
  wf(1) = prefact*(2d0*wf2(1)*TMP2-wf1(1)*TMP3-wf3(1)*TMP1)
  wf(2) = prefact*(2d0*wf2(2)*TMP2-wf1(2)*TMP3-wf3(2)*TMP1)
  wf(3) = prefact*(2d0*wf2(3)*TMP2-wf1(3)*TMP3-wf3(3)*TMP1)
  wf(4) = prefact*(2d0*wf2(4)*TMP2-wf1(4)*TMP3-wf3(4)*TMP1)
end subroutine gluon4

subroutine ipnext(ia, n, flag)
  ! Compute next permutation starting from the list 'ia' (of length
  ! 'n').
  implicit none
  integer ::n,flag,i,j,itemp
  integer, dimension(*) :: ia
  i=n-1
  do while (i.ge.1 .and. ia(i).gt.ia(i+1))
     i=i-1
  enddo
  if (i.lt.1) then
     flag = -1
     return
  endif
  j = n
  do while (ia(i).gt.ia(j))
     j = j-1
  enddo
  itemp = ia(i)
  ia(i) = ia(j)
  ia(j) = itemp
  i = i+1
  j = n
  do while (i.lt.j)
     itemp = ia(i)
     ia(i) = ia(j)
     ia(j) = itemp
     i = i+1
     j = j-1
  enddo
  flag = 1
  return
end subroutine ipnext