alpgen: new package (#26713)

Co-authored-by: Bernhard Kaindl <43588962+bernhardkaindl@users.noreply.github.com>

7 files changed, 6922 insertions, 0 deletions

diff --git a/var/spack/repos/builtin/packages/alpgen/CMakeLists.txt b/var/spack/repos/builtin/packages/alpgen/CMakeLists.txt
new file mode 100644
index 0000000000..c8bd62928c
--- /dev/null
+++ b/var/spack/repos/builtin/packages/alpgen/CMakeLists.txt
@@ -0,0 +1,69 @@
+cmake_minimum_required(VERSION 2.6 FATAL_ERROR)
+#----------------------------------------------------------------------------
+project(alpgen Fortran)
+#----------------------------------------------------------------------------
+
+#---Compilation and Link options----------------------------------------------
+if(CMAKE_SYSTEM_NAME MATCHES Darwin)
+  set(CMAKE_SHARED_LIBRARY_CREATE_Fortran_FLAGS "${CMAKE_SHARED_LIBRARY_CREATE_Fortran_FLAGS} -flat_namespace -single_module -undefined dynamic_lookup")
+endif()
+
+set(CMAKE_Fortran_FLAGS_DEBUG          "-fPIC -g -fno-automatic -std=legacy")
+set(CMAKE_Fortran_FLAGS_MINSIZEREL     "-fPIC -Os -DNDEBUG -fno-automatic -std=legacy")
+set(CMAKE_Fortran_FLAGS_RELEASE        "-fPIC -O2 -DNDEBUG -fno-automatic -std=legacy")
+set(CMAKE_Fortran_FLAGS_RELWITHDEBINFO "-fPIC -O2 -g -fno-automatic -std=legacy")
+
+
+
+#---Select sources for the various libraries---------------------------------
+
+file(GLOB alpgen_sources ${PROJECT_SOURCE_DIR}/alplib/Acp.f 
+${PROJECT_SOURCE_DIR}/alplib/Acp.f ${PROJECT_SOURCE_DIR}/alplib/Aint.f ${PROJECT_SOURCE_DIR}/alplib/alpgen.f ${PROJECT_SOURCE_DIR}/alplib/alppdf.f 
+${PROJECT_SOURCE_DIR}/alplib/alputi.f ${PROJECT_SOURCE_DIR}/alplib/Asu3.f)
+
+file(GLOB alpgen_alpsho_sources ${PROJECT_SOURCE_DIR}/alplib/alpsho.f)
+
+file(GLOB sources_libs RECURSIVE ${PROJECT_SOURCE_DIR}/*lib)
+file(GLOB sources_work RECURSIVE ${PROJECT_SOURCE_DIR}/*work)
+file(GLOB sources_rest RECURSIVE ${PROJECT_SOURCE_DIR}/DOCS ${PROJECT_SOURCE_DIR}/validation ${PROJECT_SOURCE_DIR}/examples )
+file(GLOB source_readme compare prc.list README)
+
+#---Add libraries------------------------------------------------------------
+
+
+add_library(alpgen SHARED ${alpgen_sources})
+add_library(alpgen_static STATIC ${alpgen_sources})
+set_target_properties(alpgen PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/alplib)
+set_target_properties(alpgen_static PROPERTIES OUTPUT_NAME alpgen ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/alplib)
+
+add_library(alpsho SHARED ${alpgen_alpsho_sources})
+add_library(alpsho_static STATIC ${alpgen_alpsho_sources})
+set_target_properties(alpsho PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/alplib)
+set_target_properties(alpsho_static PROPERTIES OUTPUT_NAME alpsho ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/alplib)
+
+file(GLOB alpgen_herlib_sources ${PROJECT_SOURCE_DIR}/herlib/atoher.f)
+file(GLOB alpgen_pylib_sources ${PROJECT_SOURCE_DIR}/pylib/atopyt.f)
+
+add_library(atoher SHARED ${alpgen_herlib_sources})
+add_library(atoher_static STATIC ${alpgen_herlib_sources})
+set_target_properties(atoher PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/herlib)
+set_target_properties(atoher_static PROPERTIES OUTPUT_NAME atoher ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/herlib)
+
+add_library(atopyt SHARED ${alpgen_pylib_sources})
+add_library(atopyt_static STATIC ${alpgen_pylib_sources})
+set_target_properties(atopyt PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/pylib)
+set_target_properties(atopyt_static PROPERTIES OUTPUT_NAME atopyt ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/herlib)
+
+#---Install build products---------------------------------------------------
+install(TARGETS alpgen alpsho atoher atopyt DESTINATION lib)
+install(TARGETS alpgen_static alpsho_static atoher_static atopyt_static DESTINATION lib)
+
+#---create subdirectory archive with symlinks----------------------------------
+install(DIRECTORY ${sources_libs} ${sources_work} ${sources_rest}  DESTINATION ${CMAKE_INSTALL_PREFIX}/share/sources)
+install(CODE "execute_process(COMMAND cmake -E create_symlink sources ${CMAKE_INSTALL_PREFIX}/share/alpgen-author)") 
+
+install(CODE "file(MAKE_DIRECTORY ${CMAKE_INSTALL_PREFIX}/lib/archive)")
+install(CODE "execute_process(COMMAND cmake -E create_symlink ../libalpgen.a ${CMAKE_INSTALL_PREFIX}/lib/archive/libalpgen.a)")
+install(CODE "execute_process(COMMAND cmake -E create_symlink ../libalpsho.a ${CMAKE_INSTALL_PREFIX}/lib/archive/libalpsho.a)")
+install(CODE "execute_process(COMMAND cmake -E create_symlink ../libatoher.a ${CMAKE_INSTALL_PREFIX}/lib/archive/libatoher.a)")
+install(CODE "execute_process(COMMAND cmake -E create_symlink ../libatopyt.a ${CMAKE_INSTALL_PREFIX}/lib/archive/libatopyt.a)")
diff --git a/var/spack/repos/builtin/packages/alpgen/alpgen-2.1.4-sft.patch b/var/spack/repos/builtin/packages/alpgen/alpgen-2.1.4-sft.patch
new file mode 100644
index 0000000000..d83fda3cf6
--- /dev/null
+++ b/var/spack/repos/builtin/packages/alpgen/alpgen-2.1.4-sft.patch
@@ -0,0 +1,37 @@
+diff -Naur tmp/alplib/alpgen.f tmp1/alplib/alpgen.f
+--- alplib/alpgen.f	2011-03-28 08:19:18.000000000 +0200
++++ alplib/alpgen.f	2014-03-31 15:37:15.000000000 +0200
+@@ -4647,9 +4647,9 @@
+ c     
+       IMPLICIT NONE
+       integer time
+-      character*24 CTIME,now
+-c      now='Day Mon XX hh:mm:ss yyyy'
+-      now=ctime(time())
++      character*25 now
++      now='Day Mon XX hh:mm:ss yyyy'
++c      now=ctime(time())
+       end
+ 
+ 
+diff -Naur tmp/alplib/alputi.f tmp1/alplib/alputi.f
+--- alplib/alputi.f	2010-12-14 09:03:41.000000000 +0100
++++ alplib/alputi.f	2014-03-28 16:37:36.000000000 +0100
+@@ -785,7 +785,7 @@
+      &,NBIN(200),IHIS(200,100),IUSCORE(200),IOSCORE(200)
+      &,IENT(200),HAVG(200),HINT(200),HSIG(200),BOOK(200),TITLE(200)
+      &,NHIST                                                      
+-      CHARACTER TITLE*25,BOOK*3,NOW*24
++      CHARACTER TITLE*25,BOOK*3,NOW*25
+       CHARACTER*(*) LTIT,BTIT,SCALE
+       double precision hist
+       DATA INI/0/
+@@ -841,7 +841,7 @@
+      &,NBIN(200),IHIS(200,100),IUSCORE(200),IOSCORE(200)
+      &,IENT(200),HAVG(200),HINT(200),HSIG(200),BOOK(200),TITLE(200)
+      &,NHIST                                                      
+-      CHARACTER TITLE*25,BOOK*3,SCALE*3,NOW*24
++      CHARACTER TITLE*25,BOOK*3,SCALE*3,NOW*25
+       CHARACTER*(*) LTIT,BTIT,SCA
+       double precision hist
+       CHARACTER*7  PLOT(4)
diff --git a/var/spack/repos/builtin/packages/alpgen/alpgen-214-Darwin-x86_84-gfortran.patch b/var/spack/repos/builtin/packages/alpgen/alpgen-214-Darwin-x86_84-gfortran.patch
new file mode 100644
index 0000000000..015cce4c32
--- /dev/null
+++ b/var/spack/repos/builtin/packages/alpgen/alpgen-214-Darwin-x86_84-gfortran.patch
@@ -0,0 +1,15 @@
+diff --git a/compile.mk b/compile.mk
+index 3d2ce16..7e540ca 100755
+--- a/compile.mk
++++ b/compile.mk
+@@ -17,6 +17,10 @@ ifeq ($(shell uname),Darwin)
+      FFF = g77 -O1 -Wall -fno-automatic -Wno-globals -fno-backslash \
+ 	 -ffast-math
+   endif
++  ifeq ($(shell uname -m), x86_64)
++    FFF = gfortran -O1 -Wall -fno-automatic -fno-backslash -ffast-math
++    FF90 = gfortran -fno-automatic
++  endif
+   ifeq  ($(shell uname -m),i386) 
+      FFF = ifort -V -noautomatic -warn
+      FF90 = ifort -V -noautomatic -warn
diff --git a/var/spack/repos/builtin/packages/alpgen/alpgen-214.patch b/var/spack/repos/builtin/packages/alpgen/alpgen-214.patch
new file mode 100644
index 0000000000..a105488f00
--- /dev/null
+++ b/var/spack/repos/builtin/packages/alpgen/alpgen-214.patch
@@ -0,0 +1,6522 @@
+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
++
++
diff --git a/var/spack/repos/builtin/packages/alpgen/cms_build.sh b/var/spack/repos/builtin/packages/alpgen/cms_build.sh
new file mode 100755
index 0000000000..2a62eb0a68
--- /dev/null
+++ b/var/spack/repos/builtin/packages/alpgen/cms_build.sh
@@ -0,0 +1,104 @@
+#!/bin/bash
+cd 2Qphwork; make gen; cd ..
+cd 2Qwork; make gen; cd ..
+cd 4Qwork; make gen; cd ..
+cd hjetwork; make gen; cd ..
+cd Njetwork; make gen; cd ..
+cd QQhwork; make gen; cd ..
+cd topwork; make gen; cd ..
+cd vbjetwork; make gen; cd ..
+cd wcjetwork; make gen; cd ..
+cd wphjetwork; make gen; cd ..
+cd wphqqwork; make gen; cd ..
+cd wqqwork; make gen; cd ..
+cd zqqwork; make gen; cd ..
+
+cd phjetwork; make gen;
+export USRF=120_180bin
+make gen -f cmsMakefile
+export USRF=180_240bin
+make gen -f cmsMakefile
+export USRF=20_60bin
+make gen -f cmsMakefile
+export USRF=240_300bin
+make gen -f cmsMakefile
+export USRF=300_5000bin  ##this has changed it was up to 7000
+make gen -f cmsMakefile
+export USRF=60_120bin
+make gen -f cmsMakefile
+cd ..
+
+cd wjetwork; make gen;
+export USRF=0ptw100
+make gen -f cmsMakefile
+export USRF=100ptw300
+make gen -f cmsMakefile
+export USRF=300ptw800
+make gen -f cmsMakefile
+export USRF=800ptw1600
+make gen -f cmsMakefile
+export USRF=1600ptw3200 #was not in 213 patch
+make gen -f cmsMakefile
+export USRF=3200ptw5000 #was not in 213 patch
+make gen -f cmsMakefile
+export USRF=VBFHiggsTo2Tau #was not in 213 patch
+make gen -f cmsMakefile
+export USRF=2j_vbf_inv #was not in 213 patch
+make gen -f cmsMakefile
+export USRF=3j_vbf_inv #was not in 213 patch
+make gen -f cmsMakefile
+export USRF=1600ptw #was not in 212 patch
+make gen -f cmsMakefile
+cd ..
+
+cd zjetwork; make gen;
+export USRF=0ptz100
+make gen -f cmsMakefile
+export USRF=100ptz300
+make gen -f cmsMakefile
+export USRF=300ptz800
+make gen -f cmsMakefile
+export USRF=800ptz1600
+make gen -f cmsMakefile
+export USRF=1600ptz3200 #was not in 213 patch
+make gen -f cmsMakefile
+export USRF=3200ptz5000 #was not in 213 patch
+make gen -f cmsMakefile
+export USRF=VBFHiggsTo2Tau
+make gen -f cmsMakefile
+export USRF=2j_vbf_inv #was not in 213 patch
+make gen -f cmsMakefile
+export USRF=3j_vbf_inv #was not in 213 patch
+make gen -f cmsMakefile
+export USRF=1600ptz #was not in 212 patch
+make gen -f cmsMakefile
+cd ..
+
+cd Njetwork; make gen; #whole Njetwork was not in 213 patch
+export USRF=100_160
+make gen -f cmsMakefile
+export USRF=100_180
+make gen -f cmsMakefile
+export USRF=140_180
+make gen -f cmsMakefile
+export USRF=140_5600
+make gen -f cmsMakefile
+export USRF=160_200
+make gen -f cmsMakefile
+export USRF=180_250
+make gen -f cmsMakefile
+export USRF=180_5600
+make gen -f cmsMakefile
+export USRF=200_250
+make gen -f cmsMakefile
+export USRF=20_100
+make gen -f cmsMakefile
+export USRF=20_80
+make gen -f cmsMakefile
+export USRF=250_400
+make gen -f cmsMakefile
+export USRF=400_5600
+make gen -f cmsMakefile
+export USRF=80_140
+make gen -f cmsMakefile
+cd ..
\ No newline at end of file
diff --git a/var/spack/repos/builtin/packages/alpgen/cms_install.sh b/var/spack/repos/builtin/packages/alpgen/cms_install.sh
new file mode 100755
index 0000000000..ed16234a01
--- /dev/null
+++ b/var/spack/repos/builtin/packages/alpgen/cms_install.sh
@@ -0,0 +1,37 @@
+#!/bin/bash
+
+if [ $# -ne 1 ];
+    then echo "Usage: cms_install.sh <prefix>" && exit 1
+fi
+
+i=$1
+
+mkdir -p $i/bin
+mkdir -p $i/alplib
+cp zjetwork/zjet_*gen $i/bin/
+cp wjetwork/wjet_*gen $i/bin/
+cp phjetwork/phjet_*gen $i/bin/
+cp Njetwork/Njet_*gen $i/bin/
+
+cp 2Qphwork/2Qphgen $i/bin/
+cp 2Qwork/2Qgen $i/bin/
+cp 4Qwork/4Qgen $i/bin/
+
+cp hjetwork/hjetgen $i/bin/
+cp Njetwork/Njetgen $i/bin/
+cp phjetwork/phjetgen $i/bin/
+
+cp QQhwork/QQhgen $i/bin/
+cp topwork/topgen $i/bin/
+cp vbjetwork/vbjetgen $i/bin/
+
+cp wcjetwork/wcjetgen $i/bin/
+cp wjetwork/wjetgen $i/bin/
+cp wphjetwork/wphjetgen $i/bin/
+cp wphqqwork/wphqqgen $i/bin/
+cp wqqwork/wqqgen $i/bin/
+
+cp zjetwork/zjetgen $i/bin/
+cp zqqwork/zqqgen $i/bin/
+
+cp -R alplib/* $i/alplib/
\ No newline at end of file
diff --git a/var/spack/repos/builtin/packages/alpgen/package.py b/var/spack/repos/builtin/packages/alpgen/package.py
new file mode 100644
index 0000000000..129d613582
--- /dev/null
+++ b/var/spack/repos/builtin/packages/alpgen/package.py
@@ -0,0 +1,138 @@
+# Copyright 2013-2021 Lawrence Livermore National Security, LLC and other
+# Spack Project Developers. See the top-level COPYRIGHT file for details.
+#
+# SPDX-License-Identifier: (Apache-2.0 OR MIT)
+
+import os
+
+from spack import *
+
+
+class Alpgen(MakefilePackage):
+    """A collection of codes for the generation of
+       multi-parton processes in hadronic collisions."""
+
+    homepage = "http://mlm.home.cern.ch/mlm/alpgen/"
+    url      = "http://mlm.home.cern.ch/mlm/alpgen/V2.1/v214.tgz"
+
+    patch('alpgen-214.patch', when='recipe=cms')
+    patch('alpgen-214-Darwin-x86_84-gfortran.patch', when='platform=darwin recipe=cms')
+    patch('alpgen-2.1.4-sft.patch', when='recipe=sft', level=0)
+
+    depends_on('cmake', type='build', when='recipe=sft')
+
+    variant('recipe', values=('cms', 'sft'), default='sft',
+            description='Select build recipe: CMS for CMS experiment, ' +
+                        'SFT for ATLAS/LHCb/others.')
+
+    version('2.1.4', sha256='2f43f7f526793fe5f81a3a3e1adeffe21b653a7f5851efc599ed69ea13985c5e')
+
+    phases = ['cmake', 'build', 'install']
+
+    # copied from CMakePackage
+    @property
+    def build_dirname(self):
+        """Returns the directory name to use when building the package
+
+        :return: name of the subdirectory for building the package
+        """
+        return 'spack-build-%s' % self.spec.dag_hash(7)
+
+    @property
+    def build_directory(self):
+        """Returns the directory to use when building the package
+
+        :return: directory where to build the package
+        """
+        return os.path.join(self.stage.path, self.build_dirname)
+
+    @property
+    def root_cmakelists_dir(self):
+        """The relative path to the directory containing CMakeLists.txt
+
+        This path is relative to the root of the extracted tarball,
+        not to the ``build_directory``. Defaults to the current directory.
+
+        :return: directory containing CMakeLists.txt
+        """
+        return self.stage.source_path
+
+    def cmake_args(self):
+        """Produces a list containing all the arguments that must be passed to
+        cmake, except:
+
+            * CMAKE_INSTALL_PREFIX
+            * CMAKE_BUILD_TYPE
+
+        which will be set automatically.
+
+        :return: list of arguments for cmake
+        """
+        return []
+
+    @property
+    def std_cmake_args(self):
+        """Standard cmake arguments provided as a property for
+        convenience of package writers
+
+        :return: standard cmake arguments
+        """
+        # standard CMake arguments
+        std_cmake_args = CMakePackage._std_args(self)
+        std_cmake_args += getattr(self, 'cmake_flag_args', [])
+        return std_cmake_args
+
+    # end
+
+    def url_for_version(self, version):
+        root = self.url.rsplit('/', 2)[0]
+        return "{0}/V{1}/v{2}.tgz".format(root, version.up_to(2), version.joined)
+
+    def patch(self):
+        if self.spec.satisfies('recipe=sft'):
+            copy(join_path(os.path.dirname(__file__), 'CMakeLists.txt'),
+                 'CMakeLists.txt')
+
+        if self.spec.satisfies('recipe=cms'):
+            filter_file('-fno-automatic', '-fno-automatic -std=legacy', 'compile.mk')
+
+    @when('recipe=cms')
+    def cmake(self, spec, prefix):
+        return
+
+    @when('recipe=cms')
+    def build(self, spec, prefix):
+        bash = which('bash')
+        bash('./cms_build.sh')
+
+    @when('recipe=cms')
+    def install(self, spec, prefix):
+        bash = which('bash')
+        bash('./cms_install.sh', prefix)
+
+        for root, dirs, files in os.walk(prefix):
+            set_install_permissions(root)
+            for file in files:
+                set_install_permissions(file)
+
+    @when('recipe=sft')
+    def cmake(self, spec, prefix):
+        """Runs ``cmake`` in the build directory"""
+        options = self.std_cmake_args
+        options += self.cmake_args()
+        options.append(os.path.abspath(self.root_cmakelists_dir))
+        with working_dir(self.build_directory, create=True):
+            cmake_x = which('cmake')
+            cmake_x(*options)
+
+    @when('recipe=sft')
+    def build(self, spec, prefix):
+        """Make the build targets"""
+        with working_dir(self.build_directory):
+            make()
+
+    @when('recipe=sft')
+    def install(self, spec, prefix):
+        """Make the install targets"""
+        with working_dir(self.build_directory):
+            make('install')