alpgen-214.patch

diff --git a/Njetlib/Njet.inc b/Njetlib/Njet.inc
index 32dc414..67b57f5 100644
--- a/Njetlib/Njet.inc
+++ b/Njetlib/Njet.inc
@@ -8,3 +8,7 @@ c event kinematics, explicit labelling of flavours, leptons, etc
       double precision pin,pout,pjet,ptj,etaj,drjj
       common/usrevt/pin(4,2),pout(4,maxpar-2),pjet(4,maxpar),
      + ptj(maxpar),etaj(maxpar),drjj(maxpar,maxpar)
+c
+c my counters
+      integer ggclsn,gqclsn,qqclsn
+      common/mycounters/ggclsn,gqclsn,qqclsn
diff --git a/Njetwork/Njetusr_100_160.f b/Njetwork/Njetusr_100_160.f
new file mode 100644
index 0000000..5c1476d
--- /dev/null
+++ b/Njetwork/Njetusr_100_160.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.100 .or. ptj(ord(njets)).ge.160) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_100_180.f b/Njetwork/Njetusr_100_180.f
new file mode 100644
index 0000000..3ae31ff
--- /dev/null
+++ b/Njetwork/Njetusr_100_180.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.100 .or. ptj(ord(njets)).ge.180) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_140_180.f b/Njetwork/Njetusr_140_180.f
new file mode 100644
index 0000000..93d7173
--- /dev/null
+++ b/Njetwork/Njetusr_140_180.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.140 .or. ptj(ord(njets)).ge.180) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_140_5600.f b/Njetwork/Njetusr_140_5600.f
new file mode 100644
index 0000000..54c8f29
--- /dev/null
+++ b/Njetwork/Njetusr_140_5600.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.140 .or. ptj(ord(njets)).ge.ptjmax) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_160_200.f b/Njetwork/Njetusr_160_200.f
new file mode 100644
index 0000000..ef9eb20
--- /dev/null
+++ b/Njetwork/Njetusr_160_200.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.160 .or. ptj(ord(njets)).ge.200) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_180_250.f b/Njetwork/Njetusr_180_250.f
new file mode 100644
index 0000000..f80e0bb
--- /dev/null
+++ b/Njetwork/Njetusr_180_250.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.180 .or. ptj(ord(njets)).ge.250) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_180_5600.f b/Njetwork/Njetusr_180_5600.f
new file mode 100644
index 0000000..8feb15d
--- /dev/null
+++ b/Njetwork/Njetusr_180_5600.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.180 .or. ptj(ord(njets)).ge.ptjmax) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_200_250.f b/Njetwork/Njetusr_200_250.f
new file mode 100644
index 0000000..75ce624
--- /dev/null
+++ b/Njetwork/Njetusr_200_250.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.200 .or. ptj(ord(njets)).ge.250) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_20_100.f b/Njetwork/Njetusr_20_100.f
new file mode 100644
index 0000000..c08dda0
--- /dev/null
+++ b/Njetwork/Njetusr_20_100.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.20 .or. ptj(ord(njets)).ge.100) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_20_80.f b/Njetwork/Njetusr_20_80.f
new file mode 100644
index 0000000..08452b3
--- /dev/null
+++ b/Njetwork/Njetusr_20_80.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.20 .or. ptj(ord(njets)).ge.80) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_250_400.f b/Njetwork/Njetusr_250_400.f
new file mode 100644
index 0000000..e31a24f
--- /dev/null
+++ b/Njetwork/Njetusr_250_400.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.250 .or. ptj(ord(njets)).ge.400) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_400_5600.f b/Njetwork/Njetusr_400_5600.f
new file mode 100644
index 0000000..48ef581
--- /dev/null
+++ b/Njetwork/Njetusr_400_5600.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.400 .or. ptj(ord(njets)).ge.ptjmax) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/Njetusr_80_140.f b/Njetwork/Njetusr_80_140.f
new file mode 100644
index 0000000..a7c5b3c
--- /dev/null
+++ b/Njetwork/Njetusr_80_140.f
@@ -0,0 +1,143 @@
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      ptbin=10e0
+      ptmax=400e0
+      xmbin=4e0
+      xmmax=400e0
+      do i=1,njets
+        call mbook(i,'pt_jet',ptbin,0e0,ptmax)
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mbook(i*10+j,'dR',0.1,0e0,5e0)
+        enddo
+      enddo
+      end
+
+      subroutine usrcut(lnot,weight)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c                                                                c
+c     Applies kinematical cuts to the final state during the phase
+c     -space generation                                          c
+c                                                                c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer lnot,ord(10)
+      double precision cutkin(10)
+      common/loccut/cutkin
+      double precision weight
+
+      weight=1.d0
+      lnot=0
+c
+      call alusor(ptj,njets,ord,2)
+      if(ptj(ord(njets)).lt.80 .or. ptj(ord(njets)).ge.140) goto 10
+c
+ 20   return
+c
+ 10   lnot=1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer i,j
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+      do i=1,200
+         call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,njets-1
+        do j=i+1,njets
+          call mtop(i*10+j,99,'dR['//jet(i)//jet(j)//']',' ','LIN')
+        enddo
+      enddo
+c
+      close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      integer n
+      character *15 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'Njet.inc'
+      real*8 mQQ,wgt
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      do i=1,njets-1
+        do j=i+1,njets
+          call mfill(i*10+j,real(drjj(ord(njets+1-i),ord(njets+1-j)))
+     $         ,rwgt)
+        enddo
+      enddo
+c     
+      end
+
+
+
diff --git a/Njetwork/cmsMakefile b/Njetwork/cmsMakefile
new file mode 100644
index 0000000..5a72d96
--- /dev/null
+++ b/Njetwork/cmsMakefile
@@ -0,0 +1,4 @@
+include ../compile.mk
+prc=Njet
+usrfun=$(USRF)
+include ../alplib/cms_alpgen.mk
\ No newline at end of file
diff --git a/alplib/cms_alpgen.mk b/alplib/cms_alpgen.mk
new file mode 100644
index 0000000..f62763b
--- /dev/null
+++ b/alplib/cms_alpgen.mk
@@ -0,0 +1,75 @@
+# DO NOT EDIT FROM HERE ON:
+#
+# DEFINE DIRECTORY AND FILE ALIASES
+alp= ../alplib
+her= ../herlib
+prclib= ../$(prc)lib
+prcusr=.
+prcfile=$(prclib)/$(prc)
+execfile=$(prc)_$(usrfun)gen
+usrfile=$(prcusr)/$(prc)usr_$(usrfun)
+
+VF90= ../VF90
+
+# DEFINE FILE GROUPS:
+# Files used for the parton-level event genertaion:
+ALPGEN= $(alp)/alpgen.f $(alp)/Acp.f $(alp)/Aint.f $(alp)/alputi.f \
+	$(alp)/alppdf.f $(alp)/Asu3.f
+PARTON= $(prcfile).f $(usrfile).f $(ALPGEN)
+
+# Include files
+INC=  $(prclib)/$(prc).inc $(alp)/alpgen.inc
+
+# include files' dependencies
+$(PARTON): $(INC)
+
+# object files
+OBJ=$(PARTON:.f=.o) $(PARTON90:.f90=.o)
+
+
+# compilation
+%.o: %.f $(PARTON) $(INC)
+	$(FFF) -c -o $*.o $*.f 
+$(prclib)/XXX.o90 : $(alp)/A90.f90 $(prclib)/ini_$(prc).f90
+	cd $(prclib); cp $(alp)/A90.f90 XXX.f90; \
+	cat $(prclib)/ini_$(prc).f90 >> XXX.f90;\
+	$(FF90)  -c XXX.f90; cp XXX.o XXX.o90
+
+$(prclib)/XXX.o90V : $(alp)/A90.f90 $(prclib)/ini_$(prc).f90
+	cd $(VF90); cp $(alp)/A90.f90 XXX.f90;\
+        cat $(prclib)/ini_$(prc).f90 >> XXX.f90;\
+	$(FF90V) -c XXX.f90; cp XXX.o $(prclib)/XXX.o90V; mv XXX.o $(prclib); \
+	rm -f *.vo; rm -f V*.inc;
+
+# fortran77 version
+gen: $(OBJ)
+	$(FFF) -o $(execfile) $(usrfile).o $(prcfile).o \
+	$(alp)/alpgen.o $(alp)/alputi.o $(alp)/alppdf.o \
+	$(alp)/Acp.o $(alp)/Asu3.o $(alp)/Aint.o
+# fortran90 version
+gen90: $(usrfile).o $(prcfile).o $(prclib)/$(prc).inc\
+	$(alp)/alpgen.o $(alp)/alputi.o $(alp)/alppdf.o \
+	$(alp)/Aint90.o $(prclib)/XXX.o90 $(alp)/alpgen.inc 
+	$(FF90) -o $(execfile)90 $(usrfile).o $(prcfile).o \
+	$(alp)/alpgen.o $(alp)/alputi.o $(alp)/alppdf.o \
+	$(alp)/Aint90.o $(prclib)/XXX.o
+# fortran90 version, Vast/Veridian compyler
+gen90V: $(usrfile).o $(prcfile).o $(prclib)/$(prc).inc\
+	$(alp)/alpgen.o $(alp)/alputi.o $(alp)/alppdf.o \
+	$(alp)/Aint90.o $(prclib)/XXX.o90V $(alp)/alpgen.inc 
+	$(FFF) -o $(execfile)90V $(usrfile).o $(prcfile).o \
+	$(alp)/alpgen.o $(alp)/alputi.o $(alp)/alppdf.o \
+	$(prclib)/XXX.o $(alp)/Aint90.o $(VF90)/libvast90.a
+
+# DIRECTORY CLEANUP UTILITIES:
+#
+# remove object files only
+cleanobj:
+	-rm $(PARTON:.f=.o) $(PARTON90:.f90=.o) $(prcusr)/../*/*.o90*
+
+# remove object files, etc
+cleanall:
+	-rm $(OBJ) $(prcusr)/fort.* $(prcusr)/*.top $(prcusr)/*.par \
+	$(prcusr)/*.wgt $(prcusr)/*.unw $(prcusr)/*.mon \
+	$(prcusr)/*.stat $(prcusr)/../*/*.o90*
+
diff --git a/phjetwork/cmsMakefile b/phjetwork/cmsMakefile
new file mode 100644
index 0000000..9c416bd
--- /dev/null
+++ b/phjetwork/cmsMakefile
@@ -0,0 +1,4 @@
+include ../compile.mk
+prc=phjet
+usrfun=$(USRF)
+include ../alplib/cms_alpgen.mk
diff --git a/phjetwork/phjetusr_120_180bin.f b/phjetwork/phjetusr_120_180bin.f
new file mode 100644
index 0000000..fa1f4bd
--- /dev/null
+++ b/phjetwork/phjetusr_120_180bin.f
@@ -0,0 +1,137 @@
+c      data resc/1d3/
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      double precision wusr
+      integer lnot,i
+      double precision PTGMIN, PTGMAX, ptg
+      integer init
+      data init /0/
+
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+      PTGMIN = 120d0
+      PTGMAX = 180d0
+       ptg=sqrt(pphot(1,1)**2+pphot(2,1)**2) 
+c      if(cuts not passed) goto 10
+C      write (*,*)  ptg , ptp(1)
+       if(ptg.le.PTGMIN) goto 10
+       if(ptg.gt.PTGMAX) goto 10 
+ 5    return
+ 10     lnot= 1
+       return 
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+      call mbook(4,'pt_4',ptbin,0e0,ptmax)
+      call mbook(5,'pt_5',ptbin,0e0,ptmax)
+      call mbook(6,'pt_6',ptbin,0e0,ptmax)
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
diff --git a/phjetwork/phjetusr_180_240bin.f b/phjetwork/phjetusr_180_240bin.f
new file mode 100644
index 0000000..dd21cb6
--- /dev/null
+++ b/phjetwork/phjetusr_180_240bin.f
@@ -0,0 +1,137 @@
+c      data resc/1d3/
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      double precision wusr
+      integer lnot,i
+      double precision PTGMIN, PTGMAX, ptg
+      integer init
+      data init /0/
+
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+      PTGMIN = 180d0
+      PTGMAX = 240d0
+       ptg=sqrt(pphot(1,1)**2+pphot(2,1)**2) 
+c      if(cuts not passed) goto 10
+C      write (*,*)  ptg , ptp(1)
+       if(ptg.le.PTGMIN) goto 10
+       if(ptg.gt.PTGMAX) goto 10 
+ 5    return
+ 10     lnot= 1
+       return 
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+      call mbook(4,'pt_4',ptbin,0e0,ptmax)
+      call mbook(5,'pt_5',ptbin,0e0,ptmax)
+      call mbook(6,'pt_6',ptbin,0e0,ptmax)
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
diff --git a/phjetwork/phjetusr_20_60bin.f b/phjetwork/phjetusr_20_60bin.f
new file mode 100644
index 0000000..1bdd6af
--- /dev/null
+++ b/phjetwork/phjetusr_20_60bin.f
@@ -0,0 +1,137 @@
+c      data resc/1d3/
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      double precision wusr
+      integer lnot,i
+      double precision PTGMIN, PTGMAX, ptg
+      integer init
+      data init /0/
+
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+      PTGMIN = 20d0
+      PTGMAX = 60d0
+       ptg=sqrt(pphot(1,1)**2+pphot(2,1)**2) 
+c      if(cuts not passed) goto 10
+C      write (*,*)  ptg , ptp(1)
+       if(ptg.le.PTGMIN) goto 10
+       if(ptg.gt.PTGMAX) goto 10 
+ 5    return
+ 10     lnot= 1
+       return 
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+      call mbook(4,'pt_4',ptbin,0e0,ptmax)
+      call mbook(5,'pt_5',ptbin,0e0,ptmax)
+      call mbook(6,'pt_6',ptbin,0e0,ptmax)
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
diff --git a/phjetwork/phjetusr_240_300bin.f b/phjetwork/phjetusr_240_300bin.f
new file mode 100644
index 0000000..d72191a
--- /dev/null
+++ b/phjetwork/phjetusr_240_300bin.f
@@ -0,0 +1,137 @@
+c      data resc/1d3/
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      double precision wusr
+      integer lnot,i
+      double precision PTGMIN, PTGMAX, ptg
+      integer init
+      data init /0/
+
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+      PTGMIN = 240d0
+      PTGMAX = 300d0
+       ptg=sqrt(pphot(1,1)**2+pphot(2,1)**2) 
+c      if(cuts not passed) goto 10
+C      write (*,*)  ptg , ptp(1)
+       if(ptg.le.PTGMIN) goto 10
+       if(ptg.gt.PTGMAX) goto 10 
+ 5    return
+ 10     lnot= 1
+       return 
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+      call mbook(4,'pt_4',ptbin,0e0,ptmax)
+      call mbook(5,'pt_5',ptbin,0e0,ptmax)
+      call mbook(6,'pt_6',ptbin,0e0,ptmax)
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
diff --git a/phjetwork/phjetusr_300_5000bin.f b/phjetwork/phjetusr_300_5000bin.f
new file mode 100644
index 0000000..5fc57ab
--- /dev/null
+++ b/phjetwork/phjetusr_300_5000bin.f
@@ -0,0 +1,137 @@
+c      data resc/1d3/
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      double precision wusr
+      integer lnot,i
+      double precision PTGMIN, PTGMAX, ptg
+      integer init
+      data init /0/
+
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+      PTGMIN = 300d0
+      PTGMAX = 5000d0
+       ptg=sqrt(pphot(1,1)**2+pphot(2,1)**2) 
+c      if(cuts not passed) goto 10
+C      write (*,*)  ptg , ptp(1)
+       if(ptg.le.PTGMIN) goto 10
+       if(ptg.gt.PTGMAX) goto 10 
+ 5    return
+ 10     lnot= 1
+       return 
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+      call mbook(4,'pt_4',ptbin,0e0,ptmax)
+      call mbook(5,'pt_5',ptbin,0e0,ptmax)
+      call mbook(6,'pt_6',ptbin,0e0,ptmax)
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
diff --git a/phjetwork/phjetusr_60_120bin.f b/phjetwork/phjetusr_60_120bin.f
new file mode 100644
index 0000000..56cce90
--- /dev/null
+++ b/phjetwork/phjetusr_60_120bin.f
@@ -0,0 +1,137 @@
+c      data resc/1d3/
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      double precision wusr
+      integer lnot,i
+      double precision PTGMIN, PTGMAX, ptg
+      integer init
+      data init /0/
+
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+      PTGMIN = 60d0
+      PTGMAX = 120d0
+       ptg=sqrt(pphot(1,1)**2+pphot(2,1)**2) 
+c      if(cuts not passed) goto 10
+C      write (*,*)  ptg , ptp(1)
+       if(ptg.le.PTGMIN) goto 10
+       if(ptg.gt.PTGMAX) goto 10 
+ 5    return
+ 10     lnot= 1
+       return 
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+      call mbook(4,'pt_4',ptbin,0e0,ptmax)
+      call mbook(5,'pt_5',ptbin,0e0,ptmax)
+      call mbook(6,'pt_6',ptbin,0e0,ptmax)
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,njets
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'phjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,njets
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
diff --git a/wjetwork/cmsMakefile b/wjetwork/cmsMakefile
new file mode 100644
index 0000000..c493b7d
--- /dev/null
+++ b/wjetwork/cmsMakefile
@@ -0,0 +1,4 @@
+include ../compile.mk
+prc=wjet
+usrfun=$(USRF)
+include ../alplib/cms_alpgen.mk
\ No newline at end of file
diff --git a/wjetwork/wjetusr_0ptw100.f b/wjetwork/wjetusr_0ptw100.f
new file mode 100644
index 0000000..34cadb8
--- /dev/null
+++ b/wjetwork/wjetusr_0ptw100.f
@@ -0,0 +1,157 @@
+c-------------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-------------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     pt(lept)>ptlmin  etmiss > minetmiss 
+c     abs(eta(lept)) < etalmax 
+c     lepton/jet isolation 
+c
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer lnot
+      double precision wusr
+
+c     smaria@cern.ch sep.2005 for tails cut
+      real ptw
+c
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point. 
+c     if(cut-not-passed) goto 10
+c     smaria@cern.ch sep.2005 for tails cut
+      ptw=sqrt(pw(1)**2+pw(2)**2)
+      if(ptw.gt.100) goto 10
+
+      return
+ 10      lnot= 1
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      character*1 ijet(6)
+      integer i
+      data ijet/'1','2','3','4','5','6'/
+      ptbin=2.5e0
+      ptmax=100*ptbin
+      xmbin=4e0
+      xmmax=400e0
+c
+      do i=1,min(5,njets)
+        call mbook(i,'pt j'//ijet(i),ptbin,0e0,ptmax)
+        call mbook(5+i,'eta j'//ijet(i),0.1,-3e0,3e0)
+      enddo
+      call mbook(12,'ptlept',2.,0e0,200.)
+      call mbook(13,'mW',0.5,70.,110.)
+      call mbook(14,'etal',0.2,-5.,5.)
+      end
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c debug
+      integer idbg
+      double precision fcount
+      common/fldbg/fcount(16),idbg
+      data idbg/0/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(5,njets)
+        call mtop(i,99,'pt j'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(5,njets)
+        call mtop(5+i,99,'eta j'//jet(i),' ','LIN')
+      enddo
+c
+      call mtop(12,99,'ptl',' ','LIN')
+      call mtop(13,99,'mW',' ','LIN')
+      call mtop(14,99,'etal',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,1000000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real*8 wgt,xmw
+      real rwgt
+      integer i,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+      call mfill(12,real(ptlep),rwgt)
+      if(njets.eq.0) return
+      call alusor(ptj,njets,ord,2)
+      do i=1,min(5,njets)
+        call mfill(i,real(ptj(ord(njets-i+1))),rwgt)
+        call mfill(5+i,real(etaj(ord(njets-i+1))),rwgt)
+      enddo
+      xmw=sqrt(pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2)
+      call mfill(13,real(xmw),rwgt)
+      call mfill(14,real(etalep),rwgt)
+      end
+      
diff --git a/wjetwork/wjetusr_100ptw300.f b/wjetwork/wjetusr_100ptw300.f
new file mode 100644
index 0000000..5659736
--- /dev/null
+++ b/wjetwork/wjetusr_100ptw300.f
@@ -0,0 +1,158 @@
+c-------------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-------------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     pt(lept)>ptlmin  etmiss > minetmiss 
+c     abs(eta(lept)) < etalmax 
+c     lepton/jet isolation 
+c
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer lnot
+      double precision wusr
+
+c     smaria@cern.ch sep.2005 for tails cut
+      real ptw
+c
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point. 
+c     if(cut-not-passed) goto 10
+c     smaria@cern.ch sep.2005 for tails cut
+      ptw=sqrt(pw(1)**2+pw(2)**2)
+      if(ptw.le.100) goto 10
+      if(ptw.gt.300) goto 10
+
+      return
+ 10      lnot= 1
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      character*1 ijet(6)
+      integer i
+      data ijet/'1','2','3','4','5','6'/
+      ptbin=2.5e0
+      ptmax=100*ptbin
+      xmbin=4e0
+      xmmax=400e0
+c
+      do i=1,min(5,njets)
+        call mbook(i,'pt j'//ijet(i),ptbin,0e0,ptmax)
+        call mbook(5+i,'eta j'//ijet(i),0.1,-3e0,3e0)
+      enddo
+      call mbook(12,'ptlept',2.,0e0,200.)
+      call mbook(13,'mW',0.5,70.,110.)
+      call mbook(14,'etal',0.2,-5.,5.)
+      end
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c debug
+      integer idbg
+      double precision fcount
+      common/fldbg/fcount(16),idbg
+      data idbg/0/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(5,njets)
+        call mtop(i,99,'pt j'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(5,njets)
+        call mtop(5+i,99,'eta j'//jet(i),' ','LIN')
+      enddo
+c
+      call mtop(12,99,'ptl',' ','LIN')
+      call mtop(13,99,'mW',' ','LIN')
+      call mtop(14,99,'etal',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,1000000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real*8 wgt,xmw
+      real rwgt
+      integer i,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+      call mfill(12,real(ptlep),rwgt)
+      if(njets.eq.0) return
+      call alusor(ptj,njets,ord,2)
+      do i=1,min(5,njets)
+        call mfill(i,real(ptj(ord(njets-i+1))),rwgt)
+        call mfill(5+i,real(etaj(ord(njets-i+1))),rwgt)
+      enddo
+      xmw=sqrt(pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2)
+      call mfill(13,real(xmw),rwgt)
+      call mfill(14,real(etalep),rwgt)
+      end
+      
diff --git a/wjetwork/wjetusr_1600ptw.f b/wjetwork/wjetusr_1600ptw.f
new file mode 100644
index 0000000..7f7b972
--- /dev/null
+++ b/wjetwork/wjetusr_1600ptw.f
@@ -0,0 +1,157 @@
+c-------------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-------------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     pt(lept)>ptlmin  etmiss > minetmiss 
+c     abs(eta(lept)) < etalmax 
+c     lepton/jet isolation 
+c
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer lnot
+      double precision wusr
+
+c     smaria@cern.ch sep.2005 for tails cut
+      real ptw
+c
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point. 
+c     if(cut-not-passed) goto 10
+c     smaria@cern.ch sep.2005 for tails cut
+      ptw=sqrt(pw(1)**2+pw(2)**2)
+      if(ptw.le.1600) goto 10
+
+      return
+ 10      lnot= 1
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      character*1 ijet(6)
+      integer i
+      data ijet/'1','2','3','4','5','6'/
+      ptbin=2.5e0
+      ptmax=100*ptbin
+      xmbin=4e0
+      xmmax=400e0
+c
+      do i=1,min(5,njets)
+        call mbook(i,'pt j'//ijet(i),ptbin,0e0,ptmax)
+        call mbook(5+i,'eta j'//ijet(i),0.1,-3e0,3e0)
+      enddo
+      call mbook(12,'ptlept',2.,0e0,200.)
+      call mbook(13,'mW',0.5,70.,110.)
+      call mbook(14,'etal',0.2,-5.,5.)
+      end
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c debug
+      integer idbg
+      double precision fcount
+      common/fldbg/fcount(16),idbg
+      data idbg/0/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(5,njets)
+        call mtop(i,99,'pt j'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(5,njets)
+        call mtop(5+i,99,'eta j'//jet(i),' ','LIN')
+      enddo
+c
+      call mtop(12,99,'ptl',' ','LIN')
+      call mtop(13,99,'mW',' ','LIN')
+      call mtop(14,99,'etal',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,1000000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real*8 wgt,xmw
+      real rwgt
+      integer i,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+      call mfill(12,real(ptlep),rwgt)
+      if(njets.eq.0) return
+      call alusor(ptj,njets,ord,2)
+      do i=1,min(5,njets)
+        call mfill(i,real(ptj(ord(njets-i+1))),rwgt)
+        call mfill(5+i,real(etaj(ord(njets-i+1))),rwgt)
+      enddo
+      xmw=sqrt(pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2)
+      call mfill(13,real(xmw),rwgt)
+      call mfill(14,real(etalep),rwgt)
+      end
+      
diff --git a/wjetwork/wjetusr_1600ptw3200.f b/wjetwork/wjetusr_1600ptw3200.f
new file mode 100644
index 0000000..395b816
--- /dev/null
+++ b/wjetwork/wjetusr_1600ptw3200.f
@@ -0,0 +1,158 @@
+c-------------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-------------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     pt(lept)>ptlmin  etmiss > minetmiss 
+c     abs(eta(lept)) < etalmax 
+c     lepton/jet isolation 
+c
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer lnot
+      double precision wusr
+
+c     smaria@cern.ch sep.2005 for tails cut
+      real ptw
+c
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point. 
+c     if(cut-not-passed) goto 10
+c     smaria@cern.ch sep.2005 for tails cut
+      ptw=sqrt(pw(1)**2+pw(2)**2)
+      if(ptw.le.1600) goto 10
+      if(ptw.gt.3200) goto 10
+
+      return
+ 10      lnot= 1
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      character*1 ijet(6)
+      integer i
+      data ijet/'1','2','3','4','5','6'/
+      ptbin=2.5e0
+      ptmax=100*ptbin
+      xmbin=4e0
+      xmmax=400e0
+c
+      do i=1,min(5,njets)
+        call mbook(i,'pt j'//ijet(i),ptbin,0e0,ptmax)
+        call mbook(5+i,'eta j'//ijet(i),0.1,-3e0,3e0)
+      enddo
+      call mbook(12,'ptlept',2.,0e0,200.)
+      call mbook(13,'mW',0.5,70.,110.)
+      call mbook(14,'etal',0.2,-5.,5.)
+      end
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c debug
+      integer idbg
+      double precision fcount
+      common/fldbg/fcount(16),idbg
+      data idbg/0/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(5,njets)
+        call mtop(i,99,'pt j'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(5,njets)
+        call mtop(5+i,99,'eta j'//jet(i),' ','LIN')
+      enddo
+c
+      call mtop(12,99,'ptl',' ','LIN')
+      call mtop(13,99,'mW',' ','LIN')
+      call mtop(14,99,'etal',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,1000000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real*8 wgt,xmw
+      real rwgt
+      integer i,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+      call mfill(12,real(ptlep),rwgt)
+      if(njets.eq.0) return
+      call alusor(ptj,njets,ord,2)
+      do i=1,min(5,njets)
+        call mfill(i,real(ptj(ord(njets-i+1))),rwgt)
+        call mfill(5+i,real(etaj(ord(njets-i+1))),rwgt)
+      enddo
+      xmw=sqrt(pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2)
+      call mfill(13,real(xmw),rwgt)
+      call mfill(14,real(etalep),rwgt)
+      end
+      
diff --git a/wjetwork/wjetusr_2j_vbf_inv.f b/wjetwork/wjetusr_2j_vbf_inv.f
new file mode 100644
index 0000000..430eebe
--- /dev/null
+++ b/wjetwork/wjetusr_2j_vbf_inv.f
@@ -0,0 +1,219 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-)
+c     lepton/jet isolation (if l+l-)
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      double precision wusr
+      integer lnot
+c local params
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj
+      integer i,imax,imin,j1,j2
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point.
+c     if(cut-not-passed) goto 10
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+ 
+      if((etaj(j1)-etaj(j2)).ge.0.and.(etaj(j1)-etaj(j2)).le. 2) goto 10
+      if((etaj(j1)-etaj(j2)).le.0.and.(etaj(j1)-etaj(j2)).ge.-2) goto 10
+ 
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2 -
+     &             (pjet(1,j1)+pjet(1,j2))**2 -
+     &             (pjet(2,j1)+pjet(2,j2))**2 -
+     &             (pjet(3,j1)+pjet(3,j2))**2 )
+ 
+      if(xmjj.le.300) goto 10
+c
+ 5    return
+c     if(cut-not-passed) goto 10
+ 10   lnot= 1
+      return
+      end
+ 
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+      call mbook(15,'detajj',0.2,0.,10.)
+      call mbook(16,'mjj',100.,0.,5000.)
+ 
+      end
+ 
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+      call mtop(15,99,'detajj',' ','LIN')
+      call mtop(16,99,'mjj',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+ 
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif
+      endif
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real*8 wgt,tmp,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c ===================================================
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj, deta
+      integer j1,j2
+c ===================================================
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+c ===================================================
+c  USR will add possible extra cuts at this point.
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+c
+      deta = abs(etaj(j2)-etaj(j1))
+      if(abs(etaj(j1)-etaj(j2)).ge.abs(etaj(j2)-etaj(j1))) then
+         deta = abs(etaj(j1)-etaj(j2))
+      endif
+c
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2
+     &           - (pjet(1,j1)+pjet(1,j2))**2
+     &           - (pjet(2,j1)+pjet(2,j2))**2
+     &           - (pjet(3,j1)+pjet(3,j2))**2 )
+c
+      call mfill(15,real(deta),rwgt)
+      call mfill(16,real(xmjj),rwgt)
+c ===================================================
+      end
+ 
+
diff --git a/wjetwork/wjetusr_300ptw800.f b/wjetwork/wjetusr_300ptw800.f
new file mode 100644
index 0000000..84917d2
--- /dev/null
+++ b/wjetwork/wjetusr_300ptw800.f
@@ -0,0 +1,158 @@
+c-------------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-------------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     pt(lept)>ptlmin  etmiss > minetmiss 
+c     abs(eta(lept)) < etalmax 
+c     lepton/jet isolation 
+c
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer lnot
+      double precision wusr
+
+c     smaria@cern.ch sep.2005 for tails cut
+      real ptw
+c
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point. 
+c     if(cut-not-passed) goto 10
+c     smaria@cern.ch sep.2005 for tails cut
+      ptw=sqrt(pw(1)**2+pw(2)**2)
+      if(ptw.le.300) goto 10
+      if(ptw.gt.800) goto 10
+
+      return
+ 10      lnot= 1
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      character*1 ijet(6)
+      integer i
+      data ijet/'1','2','3','4','5','6'/
+      ptbin=2.5e0
+      ptmax=100*ptbin
+      xmbin=4e0
+      xmmax=400e0
+c
+      do i=1,min(5,njets)
+        call mbook(i,'pt j'//ijet(i),ptbin,0e0,ptmax)
+        call mbook(5+i,'eta j'//ijet(i),0.1,-3e0,3e0)
+      enddo
+      call mbook(12,'ptlept',2.,0e0,200.)
+      call mbook(13,'mW',0.5,70.,110.)
+      call mbook(14,'etal',0.2,-5.,5.)
+      end
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c debug
+      integer idbg
+      double precision fcount
+      common/fldbg/fcount(16),idbg
+      data idbg/0/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(5,njets)
+        call mtop(i,99,'pt j'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(5,njets)
+        call mtop(5+i,99,'eta j'//jet(i),' ','LIN')
+      enddo
+c
+      call mtop(12,99,'ptl',' ','LIN')
+      call mtop(13,99,'mW',' ','LIN')
+      call mtop(14,99,'etal',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,1000000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real*8 wgt,xmw
+      real rwgt
+      integer i,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+      call mfill(12,real(ptlep),rwgt)
+      if(njets.eq.0) return
+      call alusor(ptj,njets,ord,2)
+      do i=1,min(5,njets)
+        call mfill(i,real(ptj(ord(njets-i+1))),rwgt)
+        call mfill(5+i,real(etaj(ord(njets-i+1))),rwgt)
+      enddo
+      xmw=sqrt(pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2)
+      call mfill(13,real(xmw),rwgt)
+      call mfill(14,real(etalep),rwgt)
+      end
+ 
diff --git a/wjetwork/wjetusr_3200ptw5000.f b/wjetwork/wjetusr_3200ptw5000.f
new file mode 100644
index 0000000..f4ce29b
--- /dev/null
+++ b/wjetwork/wjetusr_3200ptw5000.f
@@ -0,0 +1,158 @@
+c-------------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-------------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     pt(lept)>ptlmin  etmiss > minetmiss 
+c     abs(eta(lept)) < etalmax 
+c     lepton/jet isolation 
+c
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer lnot
+      double precision wusr
+
+c     smaria@cern.ch sep.2005 for tails cut
+      real ptw
+c
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point. 
+c     if(cut-not-passed) goto 10
+c     smaria@cern.ch sep.2005 for tails cut
+      ptw=sqrt(pw(1)**2+pw(2)**2)
+      if(ptw.le.3200) goto 10
+      if(ptw.gt.5000) goto 10
+
+      return
+ 10      lnot= 1
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      character*1 ijet(6)
+      integer i
+      data ijet/'1','2','3','4','5','6'/
+      ptbin=2.5e0
+      ptmax=100*ptbin
+      xmbin=4e0
+      xmmax=400e0
+c
+      do i=1,min(5,njets)
+        call mbook(i,'pt j'//ijet(i),ptbin,0e0,ptmax)
+        call mbook(5+i,'eta j'//ijet(i),0.1,-3e0,3e0)
+      enddo
+      call mbook(12,'ptlept',2.,0e0,200.)
+      call mbook(13,'mW',0.5,70.,110.)
+      call mbook(14,'etal',0.2,-5.,5.)
+      end
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c debug
+      integer idbg
+      double precision fcount
+      common/fldbg/fcount(16),idbg
+      data idbg/0/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(5,njets)
+        call mtop(i,99,'pt j'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(5,njets)
+        call mtop(5+i,99,'eta j'//jet(i),' ','LIN')
+      enddo
+c
+      call mtop(12,99,'ptl',' ','LIN')
+      call mtop(13,99,'mW',' ','LIN')
+      call mtop(14,99,'etal',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,1000000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real*8 wgt,xmw
+      real rwgt
+      integer i,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+      call mfill(12,real(ptlep),rwgt)
+      if(njets.eq.0) return
+      call alusor(ptj,njets,ord,2)
+      do i=1,min(5,njets)
+        call mfill(i,real(ptj(ord(njets-i+1))),rwgt)
+        call mfill(5+i,real(etaj(ord(njets-i+1))),rwgt)
+      enddo
+      xmw=sqrt(pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2)
+      call mfill(13,real(xmw),rwgt)
+      call mfill(14,real(etalep),rwgt)
+      end
+      
diff --git a/wjetwork/wjetusr_3j_vbf_inv.f b/wjetwork/wjetusr_3j_vbf_inv.f
new file mode 100644
index 0000000..430eebe
--- /dev/null
+++ b/wjetwork/wjetusr_3j_vbf_inv.f
@@ -0,0 +1,219 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-)
+c     lepton/jet isolation (if l+l-)
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      double precision wusr
+      integer lnot
+c local params
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj
+      integer i,imax,imin,j1,j2
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point.
+c     if(cut-not-passed) goto 10
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+ 
+      if((etaj(j1)-etaj(j2)).ge.0.and.(etaj(j1)-etaj(j2)).le. 2) goto 10
+      if((etaj(j1)-etaj(j2)).le.0.and.(etaj(j1)-etaj(j2)).ge.-2) goto 10
+ 
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2 -
+     &             (pjet(1,j1)+pjet(1,j2))**2 -
+     &             (pjet(2,j1)+pjet(2,j2))**2 -
+     &             (pjet(3,j1)+pjet(3,j2))**2 )
+ 
+      if(xmjj.le.300) goto 10
+c
+ 5    return
+c     if(cut-not-passed) goto 10
+ 10   lnot= 1
+      return
+      end
+ 
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+      call mbook(15,'detajj',0.2,0.,10.)
+      call mbook(16,'mjj',100.,0.,5000.)
+ 
+      end
+ 
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+      call mtop(15,99,'detajj',' ','LIN')
+      call mtop(16,99,'mjj',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+ 
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif
+      endif
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real*8 wgt,tmp,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c ===================================================
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj, deta
+      integer j1,j2
+c ===================================================
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+c ===================================================
+c  USR will add possible extra cuts at this point.
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+c
+      deta = abs(etaj(j2)-etaj(j1))
+      if(abs(etaj(j1)-etaj(j2)).ge.abs(etaj(j2)-etaj(j1))) then
+         deta = abs(etaj(j1)-etaj(j2))
+      endif
+c
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2
+     &           - (pjet(1,j1)+pjet(1,j2))**2
+     &           - (pjet(2,j1)+pjet(2,j2))**2
+     &           - (pjet(3,j1)+pjet(3,j2))**2 )
+c
+      call mfill(15,real(deta),rwgt)
+      call mfill(16,real(xmjj),rwgt)
+c ===================================================
+      end
+ 
+
diff --git a/wjetwork/wjetusr_800ptw1600.f b/wjetwork/wjetusr_800ptw1600.f
new file mode 100644
index 0000000..5f47b04
--- /dev/null
+++ b/wjetwork/wjetusr_800ptw1600.f
@@ -0,0 +1,158 @@
+c-------------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-------------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     pt(lept)>ptlmin  etmiss > minetmiss 
+c     abs(eta(lept)) < etalmax 
+c     lepton/jet isolation 
+c
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer lnot
+      double precision wusr
+
+c     smaria@cern.ch sep.2005 for tails cut
+      real ptw
+c
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point. 
+c     if(cut-not-passed) goto 10
+c     smaria@cern.ch sep.2005 for tails cut
+      ptw=sqrt(pw(1)**2+pw(2)**2)
+      if(ptw.le.800) goto 10
+      if(ptw.gt.1600) goto 10
+
+      return
+ 10      lnot= 1
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      character*1 ijet(6)
+      integer i
+      data ijet/'1','2','3','4','5','6'/
+      ptbin=2.5e0
+      ptmax=100*ptbin
+      xmbin=4e0
+      xmmax=400e0
+c
+      do i=1,min(5,njets)
+        call mbook(i,'pt j'//ijet(i),ptbin,0e0,ptmax)
+        call mbook(5+i,'eta j'//ijet(i),0.1,-3e0,3e0)
+      enddo
+      call mbook(12,'ptlept',2.,0e0,200.)
+      call mbook(13,'mW',0.5,70.,110.)
+      call mbook(14,'etal',0.2,-5.,5.)
+      end
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c debug
+      integer idbg
+      double precision fcount
+      common/fldbg/fcount(16),idbg
+      data idbg/0/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(5,njets)
+        call mtop(i,99,'pt j'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(5,njets)
+        call mtop(5+i,99,'eta j'//jet(i),' ','LIN')
+      enddo
+c
+      call mtop(12,99,'ptl',' ','LIN')
+      call mtop(13,99,'mW',' ','LIN')
+      call mtop(14,99,'etal',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,1000000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real*8 wgt,xmw
+      real rwgt
+      integer i,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+      call mfill(12,real(ptlep),rwgt)
+      if(njets.eq.0) return
+      call alusor(ptj,njets,ord,2)
+      do i=1,min(5,njets)
+        call mfill(i,real(ptj(ord(njets-i+1))),rwgt)
+        call mfill(5+i,real(etaj(ord(njets-i+1))),rwgt)
+      enddo
+      xmw=sqrt(pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2)
+      call mfill(13,real(xmw),rwgt)
+      call mfill(14,real(etalep),rwgt)
+      end
+
diff --git a/wjetwork/wjetusr_VBFHiggsTo2Tau.f b/wjetwork/wjetusr_VBFHiggsTo2Tau.f
new file mode 100644
index 0000000..c8d65f8
--- /dev/null
+++ b/wjetwork/wjetusr_VBFHiggsTo2Tau.f
@@ -0,0 +1,231 @@
+c-------------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-------------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     pt(lept)>ptlmin  etmiss > minetmiss 
+c     abs(eta(lept)) < etalmax 
+c     lepton/jet isolation 
+c
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer lnot
+      double precision wusr
+c local params
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj
+      integer i,j1,j2
+c
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point. 
+c     if(cut-not-passed) goto 10
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection      
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+
+      if((etaj(j1)-etaj(j2)).ge.0.and.(etaj(j1)-etaj(j2)).le. 2) goto 10
+      if((etaj(j1)-etaj(j2)).le.0.and.(etaj(j1)-etaj(j2)).ge.-2) goto 10
+
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2 -
+     &             (pjet(1,j1)+pjet(1,j2))**2 -      
+     &             (pjet(2,j1)+pjet(2,j2))**2 -      
+     &             (pjet(3,j1)+pjet(3,j2))**2 )
+
+      if(xmjj.le.300) goto 10
+
+ 5    return
+c     if(cut-not-passed) goto 10
+
+      return
+ 10   lnot= 1
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      character*1 ijet(6)
+      integer i
+      data ijet/'1','2','3','4','5','6'/
+      ptbin=2.5e0
+      ptmax=100*ptbin
+      xmbin=4e0
+      xmmax=400e0
+c
+      do i=1,min(5,njets)
+        call mbook(i,'pt j'//ijet(i),ptbin,0e0,ptmax)
+        call mbook(5+i,'eta j'//ijet(i),0.1,-3e0,3e0)
+      enddo
+      call mbook(12,'ptlept',2.,0e0,200.)
+      call mbook(13,'mW',0.5,70.,110.)
+      call mbook(14,'etal',0.2,-5.,5.)
+      call mbook(15,'detajj',0.2,0.,10.)
+      call mbook(16,'mjj',100.,0.,5000.)
+      end
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c debug
+      integer idbg
+      double precision fcount
+      common/fldbg/fcount(16),idbg
+      data idbg/0/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(5,njets)
+        call mtop(i,99,'pt j'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(5,njets)
+        call mtop(5+i,99,'eta j'//jet(i),' ','LIN')
+      enddo
+c
+      call mtop(12,99,'ptl',' ','LIN')
+      call mtop(13,99,'mW',' ','LIN')
+      call mtop(14,99,'etal',' ','LIN')
+      call mtop(15,99,'detajj',' ','LIN')
+      call mtop(16,99,'mjj',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,1000000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'wjet.inc'
+      real*8 wgt,xmw
+      real rwgt
+      integer i,jproc,ord(10)
+c ===================================================
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj, deta
+      integer j1,j2
+c ===================================================
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+      call mfill(12,real(ptlep),rwgt)
+      if(njets.eq.0) return
+      call alusor(ptj,njets,ord,2)
+      do i=1,min(5,njets)
+        call mfill(i,real(ptj(ord(njets-i+1))),rwgt)
+        call mfill(5+i,real(etaj(ord(njets-i+1))),rwgt)
+      enddo
+      xmw=sqrt(pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2)
+      call mfill(13,real(xmw),rwgt)
+      call mfill(14,real(etalep),rwgt)
+c ===================================================
+c  USR will add possible extra cuts at this point. 
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection      
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+c
+      deta = abs(etaj(j2)-etaj(j1))
+      if(abs(etaj(j1)-etaj(j2)).ge.abs(etaj(j2)-etaj(j1))) then
+         deta = abs(etaj(j1)-etaj(j2))
+      endif
+c
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2 
+     &           - (pjet(1,j1)+pjet(1,j2))**2       
+     &           - (pjet(2,j1)+pjet(2,j2))**2       
+     &           - (pjet(3,j1)+pjet(3,j2))**2 )
+c
+      call mfill(15,real(deta),rwgt)
+      call mfill(16,real(xmjj),rwgt)
+c ===================================================
+      end
+      
diff --git a/zjetwork/cmsMakefile b/zjetwork/cmsMakefile
new file mode 100644
index 0000000..eb67e38
--- /dev/null
+++ b/zjetwork/cmsMakefile
@@ -0,0 +1,4 @@
+include ../compile.mk
+prc=zjet
+usrfun=$(USRF)
+include ../alplib/cms_alpgen.mk
diff --git a/zjetwork/zjetusr_0ptz100.f b/zjetwork/zjetusr_0ptz100.f
new file mode 100644
index 0000000..d2a8fd1
--- /dev/null
+++ b/zjetwork/zjetusr_0ptz100.f
@@ -0,0 +1,149 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-) 
+c     lepton/jet isolation (if l+l-) 
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      double precision wusr
+      integer lnot
+
+      real ptz
+
+c local params
+      double precision emax,emin
+      integer i,imax,imin
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+
+c     if(cut-not-passed) goto 10
+c      for tails cut
+      ptz=sqrt(pz(1)**2+pz(2)**2)
+      if(ptz.lt.0) goto 10
+      if(ptz.gt.100) goto 10
+      return
+ 10   lnot= 1
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
+
+
diff --git a/zjetwork/zjetusr_100ptz300.f b/zjetwork/zjetusr_100ptz300.f
new file mode 100644
index 0000000..8270cc3
--- /dev/null
+++ b/zjetwork/zjetusr_100ptz300.f
@@ -0,0 +1,149 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-) 
+c     lepton/jet isolation (if l+l-) 
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      double precision wusr
+      integer lnot
+
+      real ptz
+
+c local params
+      double precision emax,emin
+      integer i,imax,imin
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+
+c     if(cut-not-passed) goto 10
+c      for tails cut
+      ptz=sqrt(pz(1)**2+pz(2)**2)
+      if(ptz.lt.100) goto 10
+      if(ptz.gt.300) goto 10
+      return
+ 10   lnot= 1
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
+
+
diff --git a/zjetwork/zjetusr_1600ptz.f b/zjetwork/zjetusr_1600ptz.f
new file mode 100644
index 0000000..ed6b814
--- /dev/null
+++ b/zjetwork/zjetusr_1600ptz.f
@@ -0,0 +1,148 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-) 
+c     lepton/jet isolation (if l+l-) 
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      double precision wusr
+      integer lnot
+
+      real ptz
+
+c local params
+      double precision emax,emin
+      integer i,imax,imin
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+
+c     if(cut-not-passed) goto 10
+c      for tails cut
+      ptz=sqrt(pz(1)**2+pz(2)**2)
+      if(ptz.lt.1600) goto 10
+      return
+ 10   lnot= 1
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
+
+
diff --git a/zjetwork/zjetusr_1600ptz3200.f b/zjetwork/zjetusr_1600ptz3200.f
new file mode 100644
index 0000000..1f30941
--- /dev/null
+++ b/zjetwork/zjetusr_1600ptz3200.f
@@ -0,0 +1,149 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-) 
+c     lepton/jet isolation (if l+l-) 
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      double precision wusr
+      integer lnot
+
+      real ptz
+
+c local params
+      double precision emax,emin
+      integer i,imax,imin
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+
+c     if(cut-not-passed) goto 10
+c      for tails cut
+      ptz=sqrt(pz(1)**2+pz(2)**2)
+      if(ptz.lt.1600) goto 10
+      if(ptz.gt.3200) goto 10
+      return
+ 10   lnot= 1
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
+
+
diff --git a/zjetwork/zjetusr_2j_vbf_inv.f b/zjetwork/zjetusr_2j_vbf_inv.f
new file mode 100644
index 0000000..7e33f64
--- /dev/null
+++ b/zjetwork/zjetusr_2j_vbf_inv.f
@@ -0,0 +1,219 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-) 
+c     lepton/jet isolation (if l+l-) 
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      double precision wusr
+      integer lnot
+c local params
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj
+      integer i,imax,imin,j1,j2
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point. 
+c     if(cut-not-passed) goto 10
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection      
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+
+      if((etaj(j1)-etaj(j2)).ge.0.and.(etaj(j1)-etaj(j2)).le. 2) goto 10
+      if((etaj(j1)-etaj(j2)).le.0.and.(etaj(j1)-etaj(j2)).ge.-2) goto 10
+
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2 -
+     &             (pjet(1,j1)+pjet(1,j2))**2 -      
+     &             (pjet(2,j1)+pjet(2,j2))**2 -      
+     &             (pjet(3,j1)+pjet(3,j2))**2 )
+
+      if(xmjj.le.300) goto 10
+c
+ 5    return
+c     if(cut-not-passed) goto 10
+ 10   lnot= 1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+      call mbook(15,'detajj',0.2,0.,10.)
+      call mbook(16,'mjj',100.,0.,5000.)
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+      call mtop(15,99,'detajj',' ','LIN')
+      call mtop(16,99,'mjj',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c ===================================================
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj, deta
+      integer j1,j2
+c ===================================================
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+c ===================================================
+c  USR will add possible extra cuts at this point. 
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection      
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+c
+      deta = abs(etaj(j2)-etaj(j1))
+      if(abs(etaj(j1)-etaj(j2)).ge.abs(etaj(j2)-etaj(j1))) then
+         deta = abs(etaj(j1)-etaj(j2))
+      endif
+c
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2 
+     &           - (pjet(1,j1)+pjet(1,j2))**2       
+     &           - (pjet(2,j1)+pjet(2,j2))**2       
+     &           - (pjet(3,j1)+pjet(3,j2))**2 )
+c
+      call mfill(15,real(deta),rwgt)
+      call mfill(16,real(xmjj),rwgt)
+c ===================================================
+      end
+
+
diff --git a/zjetwork/zjetusr_300ptz800.f b/zjetwork/zjetusr_300ptz800.f
new file mode 100644
index 0000000..fdf4439
--- /dev/null
+++ b/zjetwork/zjetusr_300ptz800.f
@@ -0,0 +1,149 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-) 
+c     lepton/jet isolation (if l+l-) 
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      double precision wusr
+      integer lnot
+
+      real ptz
+
+c local params
+      double precision emax,emin
+      integer i,imax,imin
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+
+c     if(cut-not-passed) goto 10
+c      for tails cut
+      ptz=sqrt(pz(1)**2+pz(2)**2)
+      if(ptz.lt.300) goto 10
+      if(ptz.gt.800) goto 10
+      return
+ 10   lnot= 1
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
+
+
diff --git a/zjetwork/zjetusr_3200ptz5000.f b/zjetwork/zjetusr_3200ptz5000.f
new file mode 100644
index 0000000..a57b645
--- /dev/null
+++ b/zjetwork/zjetusr_3200ptz5000.f
@@ -0,0 +1,149 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-) 
+c     lepton/jet isolation (if l+l-) 
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      double precision wusr
+      integer lnot
+
+      real ptz
+
+c local params
+      double precision emax,emin
+      integer i,imax,imin
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+
+c     if(cut-not-passed) goto 10
+c      for tails cut
+      ptz=sqrt(pz(1)**2+pz(2)**2)
+      if(ptz.lt.3200) goto 10
+      if(ptz.gt.5000) goto 10
+      return
+ 10   lnot= 1
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
+
+
diff --git a/zjetwork/zjetusr_3j_vbf_inv.f b/zjetwork/zjetusr_3j_vbf_inv.f
new file mode 100644
index 0000000..7e33f64
--- /dev/null
+++ b/zjetwork/zjetusr_3j_vbf_inv.f
@@ -0,0 +1,219 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-) 
+c     lepton/jet isolation (if l+l-) 
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      double precision wusr
+      integer lnot
+c local params
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj
+      integer i,imax,imin,j1,j2
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point. 
+c     if(cut-not-passed) goto 10
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection      
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+
+      if((etaj(j1)-etaj(j2)).ge.0.and.(etaj(j1)-etaj(j2)).le. 2) goto 10
+      if((etaj(j1)-etaj(j2)).le.0.and.(etaj(j1)-etaj(j2)).ge.-2) goto 10
+
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2 -
+     &             (pjet(1,j1)+pjet(1,j2))**2 -      
+     &             (pjet(2,j1)+pjet(2,j2))**2 -      
+     &             (pjet(3,j1)+pjet(3,j2))**2 )
+
+      if(xmjj.le.300) goto 10
+c
+ 5    return
+c     if(cut-not-passed) goto 10
+ 10   lnot= 1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+      call mbook(15,'detajj',0.2,0.,10.)
+      call mbook(16,'mjj',100.,0.,5000.)
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+      call mtop(15,99,'detajj',' ','LIN')
+      call mtop(16,99,'mjj',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c ===================================================
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj, deta
+      integer j1,j2
+c ===================================================
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+c ===================================================
+c  USR will add possible extra cuts at this point. 
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection      
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+c
+      deta = abs(etaj(j2)-etaj(j1))
+      if(abs(etaj(j1)-etaj(j2)).ge.abs(etaj(j2)-etaj(j1))) then
+         deta = abs(etaj(j1)-etaj(j2))
+      endif
+c
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2 
+     &           - (pjet(1,j1)+pjet(1,j2))**2       
+     &           - (pjet(2,j1)+pjet(2,j2))**2       
+     &           - (pjet(3,j1)+pjet(3,j2))**2 )
+c
+      call mfill(15,real(deta),rwgt)
+      call mfill(16,real(xmjj),rwgt)
+c ===================================================
+      end
+
+
diff --git a/zjetwork/zjetusr_800ptz1600.f b/zjetwork/zjetusr_800ptz1600.f
new file mode 100644
index 0000000..e334085
--- /dev/null
+++ b/zjetwork/zjetusr_800ptz1600.f
@@ -0,0 +1,149 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-) 
+c     lepton/jet isolation (if l+l-) 
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      double precision wusr
+      integer lnot
+
+      real ptz
+
+c local params
+      double precision emax,emin
+      integer i,imax,imin
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+
+c     if(cut-not-passed) goto 10
+c      for tails cut
+      ptz=sqrt(pz(1)**2+pz(2)**2)
+      if(ptz.lt.800) goto 10
+      if(ptz.gt.1600) goto 10
+      return
+ 10   lnot= 1
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+      end
+
+
diff --git a/zjetwork/zjetusr_VBFHiggsTo2Tau.f b/zjetwork/zjetusr_VBFHiggsTo2Tau.f
new file mode 100644
index 0000000..7e33f64
--- /dev/null
+++ b/zjetwork/zjetusr_VBFHiggsTo2Tau.f
@@ -0,0 +1,219 @@
+c-----------------------------------------------------------------
+      subroutine usrcut(lnot,wusr)
+c-----------------------------------------------------------------
+c PRIMARY CUTS ALREADY APPLIED TO PHASE-SPACE GENERATION:
+c     ptjmin < pt(jet) < ptjmax for all light jets
+c     -etajmax < eta(jet) < etajmax for all light jets
+c     delta R(jj) > drjmin for all (light jet, light jet) pairs
+c     mllmin < m(l+l-) < mllmax
+c     pt(lept)>ptlmin (if l+l-) or  etmiss > minetmiss (if nu nubar)
+c     abs(eta(lept)) < etalmax (if l+l-) 
+c     lepton/jet isolation (if l+l-) 
+c USE THIS ROUTINE TO ENFORCE OTHER CUTS    
+c
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      double precision wusr
+      integer lnot
+c local params
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj
+      integer i,imax,imin,j1,j2
+c initialize output parameters
+      lnot=0
+      wusr=1d0
+c
+c  USR will add possible extra cuts at this point. 
+c     if(cut-not-passed) goto 10
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection      
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+
+      if((etaj(j1)-etaj(j2)).ge.0.and.(etaj(j1)-etaj(j2)).le. 2) goto 10
+      if((etaj(j1)-etaj(j2)).le.0.and.(etaj(j1)-etaj(j2)).ge.-2) goto 10
+
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2 -
+     &             (pjet(1,j1)+pjet(1,j2))**2 -      
+     &             (pjet(2,j1)+pjet(2,j2))**2 -      
+     &             (pjet(3,j1)+pjet(3,j2))**2 )
+
+      if(xmjj.le.300) goto 10
+c
+ 5    return
+c     if(cut-not-passed) goto 10
+ 10   lnot= 1
+      return
+      end
+
+c-------------------------------------------------------------------
+      subroutine alshis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real ptbin,ptmax,xmbin,xmmax
+      ptbin=2e0
+      ptmax=200e0
+      xmbin=4e0
+      xmmax=400e0
+c
+      call mbook(1,'pt_1',2.*ptbin,0e0,2.*ptmax)
+      call mbook(2,'pt_2',ptbin,0e0,ptmax)
+      call mbook(3,'pt_3',ptbin,0e0,ptmax)
+c
+      call mbook(11,'eta_1',0.2,-5.0,5.0)
+      call mbook(12,'eta_2',0.2,-5.0,5.0)
+      call mbook(13,'eta_3',0.2,-5.0,5.0)
+      call mbook(15,'detajj',0.2,0.,10.)
+      call mbook(16,'mjj',100.,0.,5000.)
+
+      end
+
+c-------------------------------------------------------------------
+      subroutine alfhis
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer i
+      real  xnorm
+      character *1 jet(9)
+      data jet/'1','2','3','4','5','6','7','8','9'/
+c
+      open(unit=99,file=topfile,err=101,status='unknown')
+      if(imode.le.1) then
+         xnorm=sngl(avgwgt/totwgt)
+      elseif(imode.eq.2) then
+         xnorm=1e0/real(unwev)
+      else
+         write(*,*) 'imode type not allowed, stop'
+         stop
+      endif
+c
+      do i=1,200
+         if(i.ne.61) call mopera(i,'F',i,i,xnorm,1.)
+         call mfinal(i)
+      enddo 
+c
+      do i=1,min(3,njets)
+         call mtop(i,99,'pt'//jet(i),' ','LOG')
+      enddo
+      do i=1,min(3,njets)
+         call mtop(10+i,99,'eta'//jet(i),' ','LIN')
+      enddo
+      call mtop(15,99,'detajj',' ','LIN')
+      call mtop(16,99,'mjj',' ','LIN')
+c
+ 100  close(99)
+ 101  return
+      end
+
+      subroutine monitor(n,mon_fname)
+c     This routine is called by default every 100K events.
+c     The user can use it to get regular updates on the run
+c     while this is progressing. Textual output can be written to file
+c     fname, where partial cross-sections and and generation
+c     efficiencies have already been printed by default
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      integer n
+      character *50 mon_fname
+c
+      if(evgen) then
+         if(mod(n,100000).eq.0) then
+c     save histograms' contents
+            call msave
+c     print out histograms
+            call alfhis
+c     restore original contents, to proceed with analysis
+            call mrestore
+         endif 
+      endif 
+      end
+c-------------------------------------------------------------------
+      subroutine aleana(jproc,wgt)
+c     analyse event, fill histograms
+c-------------------------------------------------------------------
+      implicit none
+      include 'alpgen.inc'
+      include 'zjet.inc'
+      real*8 wgt,tmp,ptlep,etmiss,xmz,mll
+      real rwgt
+      integer i,j,jproc,ord(10)
+c ===================================================
+      double precision emax,emin, ptjbig, ptjets(maxpar-2), xmjj, deta
+      integer j1,j2
+c ===================================================
+c
+      rwgt=real(wgt)
+      if(rwgt.lt.0e0) then
+         write(*,*) 'negative wgt=',wgt
+         return
+      elseif (rwgt.eq.0e0) then
+         return
+      endif 
+c
+c     reordering according to pt
+      call alusor(ptj,njets,ord,2)              
+      do i=1,min(3,njets)
+         call mfill(i,real(ptj(ord(njets+1-i))),rwgt)
+         call mfill(10+i,real(etaj(ord(njets+1-i))),rwgt)
+      enddo
+c
+c ===================================================
+c  USR will add possible extra cuts at this point. 
+      do i = 1, njets
+         ptjets(i) = ptj(i)
+      enddo
+c  first max pt jet selection
+      ptjbig = -10.
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j1 = i
+         endif
+      enddo
+      ptjets(j1) = -10.
+      ptjbig = -10.
+c  second max pt jet selection      
+      do i = 1, njets
+         if(ptjets(i).gt.ptjbig) then
+            ptjbig = ptjets(i)
+            j2 = i
+         endif
+      enddo
+c
+      deta = abs(etaj(j2)-etaj(j1))
+      if(abs(etaj(j1)-etaj(j2)).ge.abs(etaj(j2)-etaj(j1))) then
+         deta = abs(etaj(j1)-etaj(j2))
+      endif
+c
+      xmjj = sqrt( (pjet(4,j1)+pjet(4,j2))**2 
+     &           - (pjet(1,j1)+pjet(1,j2))**2       
+     &           - (pjet(2,j1)+pjet(2,j2))**2       
+     &           - (pjet(3,j1)+pjet(3,j2))**2 )
+c
+      call mfill(15,real(deta),rwgt)
+      call mfill(16,real(xmjj),rwgt)
+c ===================================================
+      end
+
+