diff -Naur POWHEG-hvq/POWHEG-hvq.f POWHEG-hvq-new/POWHEG-hvq.f
--- POWHEG-hvq/POWHEG-hvq.f	2007-06-15 12:34:39.000000000 +0200
+++ POWHEG-hvq-new/POWHEG-hvq.f	2007-09-19 15:14:42.000000000 +0200
@@ -1066,6 +1066,9 @@
       endif
 c     generate radiation
       call genptxy
+c     perform random azimuthal rotation of the event
+      call addrnmphi(pq1,pq2,pp)
+c
       call increasecnt('events number')
       if(jflreal.gt.0) call increasecnt('jflreal p')
       if(jflreal.lt.0) call increasecnt('jflreal n')
@@ -1353,6 +1356,38 @@
       p5(5)=sqrt(m2)
       end
 
+
+      subroutine addrnmphi(pq1,pq2,pp)
+c rotate the 2-momenta pq1,pq2 and pp by the same
+c random angle
+      implicit none      
+      real * 8 pq1(2),pq2(2),pp(2)
+      real * 8 sphi,cphi,tmp,tmp1,tmp2
+      real * 8 pi,phi
+      parameter (pi=3.141592653589793d0)
+      real * 8 random
+      external random
+      phi=2*pi*random()
+      sphi=sin(phi)
+      cphi=cos(phi)
+c Rotate the vectors
+      tmp1 = cphi*pq1(1)-sphi*pq1(2)
+      tmp2 = sphi*pq1(1)+cphi*pq1(2)
+      pq1(1)=tmp1
+      pq1(2)=tmp2
+c
+      tmp1 = cphi*pq2(1)-sphi*pq2(2)
+      tmp2 = sphi*pq2(1)+cphi*pq2(2)
+      pq2(1)=tmp1
+      pq2(2)=tmp2
+c
+      tmp1 = cphi*pp(1)-sphi*pp(2)
+      tmp2 = sphi*pp(1)+cphi*pp(2)
+      pp(1)=tmp1
+      pp(2)=tmp2
+c
+      end
+
       subroutine tolab
       implicit none
       real * 8 pwl(2),pzl(2),ppl(2)