Merge lp:~maddevelopers/mg5amcnlo/bias into lp:mg5amcnlo/lts

[Valentin Hirschi]

Dear reviewers,

This branch introduces the possibility of using custom modules that can apply a bias weight to MadEvent integrand so as to change the distribution of "unweighted" events to smoothen an observable (for example the leading jet p_t) or to really affect the physics described.

There are several applications to such a module, and the one we have in mind for now is to call Pythia8 directly from MadEvent to compute the CKKW-L merging weights.

I have written the exhaustive details on this implementation on this wiki page

https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/LOEventGenerationBias

I already worked with Olivier to decide on the better structure, so his review should be quick.
But I wanted to add you, Rik and Marco, as well as reviewers because I think there would be much greater applications if we had the same sort of module at NLO!
In particular, we could compute the Pythia8 MC counterterms by calling their trial shower directly from such a module, so that we would have guaranteed consistency with the subsequent showers and lift the current limitations on its parameters!
In the same spirit (i.e. dynamically calling PY8 from the event generator) we might open new approaches to NLO merging, Mixed QCD+EW matching, and even matching in the presence of different resonance structures (say tt~ doubly+singly+non-resonance).

This implementation at LO could provide a basis to start a similar implementation at NLO. Please let me know your opinion on the matter.

Cheers,

[lp:~maddevelopers/mg5amcnlo/bias]

305. By Olivier Mattelaer on 2016-08-23

small code improvment

306. By Valentin Hirschi on 2016-08-24

1. Fixed two tests.

307. By Valentin Hirschi on 2016-08-29

1. Merged

308. By Valentin Hirschi on 2016-08-29

1. Merged with Olivier's change for __type__

309. By Valentin Hirschi on 2016-08-29

1. merged with latest version of plugin_plus_pythia8

310. By Valentin Hirschi on 2016-08-29

1. merged with latest version of plugin_plus_pythia8

311. By Valentin Hirschi on 2016-08-30

1. Added an acceptance test for the bias functionality

312. By Olivier Mattelaer on 2016-09-01

fix problem with bias normalisation

313. By Valentin Hirschi on 2016-09-01

1. merged with plugin_plus_pythia8 and fixed bias-related tests.

314. By Valentin Hirschi on 2016-09-01

1. Removed the PY8_CKKW bias module which should not be part of the release.

315. By Valentin Hirschi on 2016-09-01

1. Merged with latest version of plugin.

Unmerged revisions

=== added directory 'PLUGIN'
=== added file 'PLUGIN/__init__.py'
=== added file 'PythiaMadEventCrossTalk.tar.gz'
Binary files PythiaMadEventCrossTalk.tar.gz	1970-01-01 00:00:00 +0000 and PythiaMadEventCrossTalk.tar.gz	2016-08-23 01:19:18 +0000 differ
=== added file 'Template/LO/Cards/pythia8_card_default.dat'
--- Template/LO/Cards/pythia8_card_default.dat	1970-01-01 00:00:00 +0000
+++ Template/LO/Cards/pythia8_card_default.dat	2016-08-23 01:19:18 +0000
@@ -0,0 +1,62 @@
+!
+! Pythia8 cmd card automatically generated by MadGraph5_aMC@NLO
+!
+! ==================
+! General parameters 
+! ==================
+!
+Main:numberOfEvents      = -1
+!
+! -------------------------------------------------------------------
+! Specify the HEPMC output of the Pythia8 shower. You can set it to:
+!   auto  : to let MG5aMC automatically assign it
+!   fifo  : to have MG5aMC setup the piping of the PY8 output to analysis
+!           tools such as MadAnalysis5. You can also use the syntax 
+!           fifo@<fifo_path> to select where the fifo file should be
+!           created (useful to select a mounted drive that supports fifo)
+!   /dev/null : to turn off the HEPMC output.
+!   <path>    : to select where the HEPMC file must written. It will 
+!               therefore *not* be available in the process output.
+! -------------------------------------------------------------------
+!
+HEPMCoutput:file         = auto
+!
+! --------------------------------------------------------------------
+! Parameters relevant only when performing MLM merging
+! For details, see section 'Jet Matching' on the left-hand menu of 
+!    http://home.thep.lu.se/~torbjorn/pythia81html/Welcome.html
+! --------------------------------------------------------------------
+! If equal to -1.0, MadGraph5_aMC@NLO will set it automatically based 
+! on the parameter 'xqcut' of the run_card.dat 
+JetMatching:qCut         = -1.0
+! Use default kt-MLM to match parton level jets to those produced by the
+! shower. But the other Shower-kt scheme is available too with this option.
+JetMatching:doShowerKt   = off
+! A value of -1 means that it is automatically guessed by MadGraph.
+! It is however always safer to explicitly set it.
+JetMatching:nJetMax      = -1
+!
+! --------------------------------------------------------------------
+! Parameters relevant only when performing CKKWL merging
+! For details, see section 'CKKW-L Merging' on the left-hand menu of 
+!    http://home.thep.lu.se/~torbjorn/pythia81html/Welcome.html
+! --------------------------------------------------------------------
+! Central merging scale values you want to be used.
+! If equal to -1.0, then MadGraph5_aMC@NLO will set this automatically
+! based on the parameter 'ktdurham' of the run_card.dat
+Merging:TMS              = -1.0
+! This must be set manually, according to Pythia8 directives.
+! An example of possible value is 'pp>LEPTONS,NEUTRINOS'
+Merging:Process          = <set_by_user>
+! A value of -1 means that it is automatically guessed by MadGraph.
+! It is however always safer to explicitly set it.
+Merging:nJetMax  		 = -1
+!
+! For all merging schemes, decide wehter you want the merging scale
+! variation computed for only the central weights or all other 
+! PDF and scale variation weights as well
+SysCalc:fullCutVariation = off
+!
+! ==========================
+! User customized parameters 
+! ==========================

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2015-08-11 00:35:53 +0000
+++ Template/LO/Cards/run_card.dat	2016-08-23 01:19:18 +0000
@@ -21,10 +21,6 @@
 #*********************************************************************
   %(run_tag)s     = run_tag ! name of the run 
 #*********************************************************************
-# Run to generate the grid pack                                      *
-#*********************************************************************
-  %(gridpack)s     = gridpack  !True = setting up the grid pack
-#*********************************************************************
 # Number of events and rnd seed                                      *
 # Warning: Do not generate more than 1M events in a single run       *
 # If you want to run Pythia, avoid more than 50k events in a run.    *
@@ -61,51 +57,56 @@
  %(dynamical_scale_choice)s = dynamical_scale_choice ! Choose one of the preselected dynamical choices
  %(scalefact)s  = scalefact        ! scale factor for event-by-event scales
 #*********************************************************************
-# Time of flight information. (-1 means not run)
-#*********************************************************************
- %(time_of_flight)s = time_of_flight ! threshold below which info is not written
-#*********************************************************************
-# Matching - Warning! ickkw > 1 is still beta
-#*********************************************************************
- %(ickkw)s = ickkw            ! 0 no matching, 1 MLM, 2 CKKW matching
- %(highestmult)s = highestmult      ! for ickkw=2, highest mult group
- %(ktscheme)s = ktscheme         ! for ickkw=1, 1 Durham kT, 2 Pythia pTE
+# Type and output format
+#*********************************************************************
+  %(gridpack)s     = gridpack  !True = setting up the grid pack
+  %(time_of_flight)s = time_of_flight ! threshold below which info is not written (-1 means not written)
+  %(lhe_version)s = lhe_version       ! Change the way clustering information pass to shower.        
+  %(clusinfo)s = clusinfo         ! include clustering tag in output
+  %(event_norm)s =  event_norm       ! average/sum. Normalization of the weight in the LHEF
+
+#*********************************************************************
+# Matching parameter (MLM only)
+#*********************************************************************
+ %(ickkw)s = ickkw            ! 0 no matching, 1 MLM
  %(alpsfact)s = alpsfact         ! scale factor for QCD emission vx
  %(chcluster)s = chcluster        ! cluster only according to channel diag
- %(pdfwgt)s = pdfwgt           ! for ickkw=1, perform pdf reweighting
  %(asrwgtflavor)s = asrwgtflavor     ! highest quark flavor for a_s reweight
- %(clusinfo)s = clusinfo         ! include clustering tag in output
- %(lhe_version)s = lhe_version       ! Change the way clustering information pass to shower.        
+ %(auto_ptj_mjj)s  = auto_ptj_mjj  ! Automatic setting of ptj and mjj if xqcut >0
+                                   ! (turn off for VBF and single top processes) 
+ %(xqcut)s   = xqcut   ! minimum kt jet measure between partons
 #*********************************************************************
-#**********************************************************
 #
-#**********************************************************
-# Automatic ptj and mjj cuts if xqcut > 0
-# (turn off for VBF and single top processes)
-#**********************************************************
-   %(auto_ptj_mjj)s  = auto_ptj_mjj  ! Automatic setting of ptj and mjj
-#**********************************************************
-#                                                                    
-#**********************************
-# BW cutoff (M+/-bwcutoff*Gamma)
-#**********************************
-  %(bwcutoff)s  = bwcutoff      ! (M+/-bwcutoff*Gamma)
-#**********************************************************
-# Apply pt/E/eta/dr/mij/kt_durham cuts on decay products or not
-# (note that etmiss/ptll/ptheavy/ht/sorted cuts always apply)
-#*************************************************************
-   %(cut_decays)s  = cut_decays    ! Cut decay products 
-#*************************************************************
+#*********************************************************************
 # Number of helicities to sum per event (0 = all helicities)
 # 0 gives more stable result, but longer run time (needed for
 # long decay chains e.g.).
 # Use >=2 if most helicities contribute, e.g. pure QCD.
-#*************************************************************
+#*********************************************************************
    %(nhel)s  = nhel          ! Number of helicities used per event
-#*******************                                                 
-# Standard Cuts
-#*******************                                                 
+#*********************************************************************
+# Generation bias, check the wiki page below for more information:   *
+#  'cp3.irmp.ucl.ac.be/projects/madgraph/wiki/LOEventGenerationBias' *
+#*********************************************************************
+ %(bias_module)s = bias_module  ! Bias type of bias, [None, ptj_bias, <custom_folder>]
+ %(bias_parameters)s = bias_parameters ! Specifies the parameters of the module.
+#
+#*******************************                                                 
+# Parton level cuts definition *
+#*******************************                                     
 #                                                                    
+#
+#*********************************************************************
+# BW cutoff (M+/-bwcutoff*Gamma) ! Define on/off-shell for "$" and decay  
+#*********************************************************************
+  %(bwcutoff)s  = bwcutoff      ! (M+/-bwcutoff*Gamma)
+#*********************************************************************
+# Apply pt/E/eta/dr/mij/kt_durham cuts on decay products or not
+# (note that etmiss/ptll/ptheavy/ht/sorted cuts always apply)
+#*********************************************************************
+   %(cut_decays)s  = cut_decays    ! Cut decay products 
+#*********************************************************************
+# Standard Cuts                                                      *
 #*********************************************************************
 # Minimum and maximum pt's (for max, -1 means no cut)                *
 #*********************************************************************
@@ -246,21 +247,20 @@
 #*********************************************************************
  %(xetamin)s   = xetamin ! minimum rapidity for two jets in the WBF case  
  %(deltaeta)s   = deltaeta ! minimum rapidity for two jets in the WBF case 
-#*********************************************************************
-# KT DURHAM CUT                                                      *
-#*********************************************************************
- %(ktdurham)s    =  ktdurham        
- %(dparameter)s  =  dparameter 
+#***********************************************************************
+# Turn on either the ktdurham or ptlund cut to activate                *
+# CKKW(L) merging with Pythia8 [arXiv:1410.3012, arXiv:1109.4829]      *
+#***********************************************************************
+ %(ktdurham)s  =  ktdurham        
+ %(dparameter)s   =  dparameter
+ %(ptlund)s  =  ptlund
+ %(pdgs_for_merging_cut)s  =  pdgs_for_merging_cut ! PDGs for two cuts above   
 #*********************************************************************
 # maximal pdg code for quark to be considered as a light jet         *
 # (otherwise b cuts are applied)                                     *
 #*********************************************************************
  %(maxjetflavor)s = maxjetflavor    ! Maximum jet pdg code
 #*********************************************************************
-# Jet measure cuts                                                   *
-#*********************************************************************
- %(xqcut)s   = xqcut   ! minimum kt jet measure between partons
-#*********************************************************************
 #
 #*********************************************************************
 # Store info for systematics studies                                 *
@@ -281,3 +281,4 @@
 # PDF sets and number of members (0 or none for all members).
 %(sys_pdf)s = sys_pdf # matching scales
 # MSTW2008nlo68cl.LHgrid 1  = sys_pdf
+#

=== added directory 'Template/LO/Source/BIAS'
=== added directory 'Template/LO/Source/BIAS/PY8_CKKW'
=== added file 'Template/LO/Source/BIAS/PY8_CKKW/PY8_CKKW.f'
--- Template/LO/Source/BIAS/PY8_CKKW/PY8_CKKW.f	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/PY8_CKKW/PY8_CKKW.f	2016-08-23 01:19:18 +0000
@@ -0,0 +1,180 @@
+C ************************************************************
+C Source for the library implementing a bias function that 
+C implements the Sudakov weight in CKKW directly from Pythia8 
+C ************************************************************
+C
+C The following lines are read by MG5aMC to set what are the 
+C relevant parameters for this bias module.
+C
+C  parameters = {'BeamA': 0.0,
+C                'BeamB': 0.0}
+C
+
+      subroutine bias_wgt(p, original_weight, bias_weight)
+          implicit none
+C
+C Parameters
+C
+          include '../../maxparticles.inc'          
+          include '../../nexternal.inc'
+C
+C Accessingt the details of the event
+C
+          include '../../run_config.inc'
+          include '../../lhe_event_infos.inc'
+C
+C Arguments
+C
+          double precision p(0:3,nexternal)
+          double precision original_weight, bias_weight
+C
+C local variables
+C
+          integer lpp_to_beam(-3:3)
+          data lpp_to_beam/-11,-2212,-2212,0,2212,2212,11/
+c
+c local variables defined in the run_card
+c
+c         Bias module arguments
+          double precision BeamA, BeamB
+
+c         truly local variables
+          integer i,j, n_initial
+          double precision OutputBiasWeight
+          double precision Pythia8eCM
+          integer Pythia8nParticles
+          double precision Pythia8p(5,npart)
+          integer Pythia8BeamA
+          double precision PythiaBeamEnergyA
+          integer Pythia8BeamB
+          double precision PythiaBeamEnergyB
+          character*(maxEventLength) Pythia8EvtRecord
+          integer Pythia8Helicities(npart)
+          integer Pythia8ColorOne(npart)
+          integer Pythia8ColorTwo(npart)          
+          integer Pythia8ID(npart)
+          integer Pythia8Status(npart)
+          integer Pythia8MotherOne(npart)
+          integer Pythia8MotherTwo(npart)          
+          integer Pythia8SubprocessGroup
+          integer Pythia8MurScale
+          integer Pythia8AlphaQCD
+          integer Pythia8AlphaQED
+
+C
+C Global variables
+C
+C
+C Mandatory common block to be defined in bias modules
+C
+          double precision stored_bias_weight
+          data stored_bias_weight/1.0d0/          
+          logical impact_xsec, requires_full_event_info
+C         We only want to bias distributions, but not impact the xsec. 
+          data impact_xsec/.True./
+C         Pythia8 will need the full information for the event
+C          (color, resonances, helicities, etc..)
+          data requires_full_event_info/.True./ 
+          common/bias/stored_bias_weight,impact_xsec,
+     &                requires_full_event_info
+C
+C Access the value of the run parameters in run_card
+C
+          include '../../run.inc'
+          include '../../cuts.inc'
+C
+C Read the definition of the bias parameter from the run_card    
+C
+          include '../bias.inc'
+
+C --------------------
+C BEGIN IMPLEMENTATION
+C --------------------
+
+C        Let's initialize the PY8 variables describing the event
+         Pythia8eCM             = sqrt(4d0*ebeam(1)*ebeam(2))
+         Pythia8EvtRecord       = event_record
+         Pythia8SubprocessGroup = ngroup
+         Pythia8MurScale        = sscale
+         Pythia8AlphaQCD        = aaqcd
+         Pythia8AlphaQED        = aaqed
+         Pythia8nParticles      = npart
+         PythiaBeamEnergyA      = ebeam(1) 
+         PythiaBeamEnergyB      = ebeam(2)
+         do i=1,npart
+           Pythia8ID(i)         = jpart(1,i)
+           Pythia8MotherOne(i)  = jpart(2,i)
+           Pythia8MotherTwo(i)  = jpart(3,i)
+           Pythia8ColorOne(i)   = jpart(4,i)
+           Pythia8ColorTwo(i)   = jpart(5,i)           
+           Pythia8Status(i)     = jpart(6,i)
+           Pythia8Helicities(i) = jpart(7,i)           
+           do j=1,4
+             Pythia8p(j,i)=pb(mod(j,4),i)
+           enddo
+           Pythia8p(5,npart)=pb(4,i)
+         enddo
+         Pythia8BeamA = lpp_to_beam(lpp(1))
+         Pythia8BeamB = lpp_to_beam(lpp(2))
+         n_initial = 0
+         do i=1,npart
+           if (Pythia8Status(i).eq.-1) then
+             n_initial = n_initial + 1
+             if (lpp(n_initial).eq.0) then
+               if (n_initial.eq.1) then
+                 Pythia8BeamA = Pythia8ID(i)
+               elseif (n_initial.eq.2) then
+                 Pythia8BeamB = Pythia8ID(i)
+                 exit
+               endif
+             endif
+             if (n_initial.eq.2) then
+               exit
+             endif
+           endif
+         enddo
+C        If this is a 1 > N decay event, then enforce beamIDs to match
+C        those specified in the event record.
+         if (n_initial.eq.1) then
+           Pythia8BeamB = 0
+           do i=1,npart         
+             if (Pythia8Status(i).eq.-1) then
+               Pythia8BeamA = Pythia8ID(i) 
+               exit
+             endif
+           enddo
+         endif
+C        Make sure to enforce the user-choice of beam if specified.
+         if (idnint(BeamA).ne.0) then
+           Pythia8BeamA = idnint(BeamA)
+         endif
+         if (idnint(BeamB).ne.0) then
+           Pythia8BeamB = idnint(BeamB) 
+         endif
+
+C        Call PY8 to derive the bias weight.
+         call py8_bias_weight( Pythia8eCM,
+     &                         Pythia8BeamA,
+     &                         PythiaBeamEnergyA,
+     &                         Pythia8BeamB,
+     &                         PythiaBeamEnergyB,
+     &                         Pythia8EvtRecord,
+     &                         Pythia8p,
+     &                         Pythia8nParticles,
+     &                         Pythia8MurScale,
+     &                         Pythia8AlphaQCD,
+     &                         Pythia8AlphaQED,
+     &                         Pythia8ID,
+     &                         Pythia8MotherOne,
+     &                         Pythia8MotherTwo,
+     &                         Pythia8ColorOne,
+     &                         Pythia8ColorTwo,
+     &                         Pythia8Status,
+     &                         Pythia8Helicities,
+     &                         OutputBiasWeight    )
+ 
+          bias_weight = OutputBiasWeight
+            
+          return
+
+      end subroutine bias_wgt

=== added file 'Template/LO/Source/BIAS/PY8_CKKW/bias_dependencies'
--- Template/LO/Source/BIAS/PY8_CKKW/bias_dependencies	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/PY8_CKKW/bias_dependencies	2016-08-23 01:19:18 +0000
@@ -0,0 +1,5 @@
+PY8LIBS     = $(shell $(PYTHIA8_PATH)/bin/pythia8-config --libs)
+PY8INCLUDES = $(shell $(PYTHIA8_PATH)/bin/pythia8-config --includedir)
+PY8CXXFLAGS = $(shell $(PYTHIA8_PATH)/bin/pythia8-config --cxxflags)
+
+BIASDEPENDENCIES = $(PY8LIBS) $(PY8CXXFLAGS)

=== added file 'Template/LO/Source/BIAS/PY8_CKKW/makefile'
--- Template/LO/Source/BIAS/PY8_CKKW/makefile	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/PY8_CKKW/makefile	2016-08-23 01:19:18 +0000
@@ -0,0 +1,59 @@
+BIASDIR   = ../
+SOURCEDIR = ../..
+include $(SOURCEDIR)/make_opts
+
+include ./bias_dependencies
+
+all: PY8_CKKW 
+
+clean:
+	$(RM) *.o $(BIASLIBDIR)$(BIASLIBRARY)
+
+INCLUDES = $(BIASDIR)/bias.inc \
+	 	   $(SOURCEDIR)/maxparticles.inc \
+	 	   $(SOURCEDIR)/nexternal.inc \
+		   $(SOURCEDIR)/run_config.inc \
+		   $(SOURCEDIR)/lhe_event_infos.inc \
+		   $(SOURCEDIR)/run.inc \
+		   $(SOURCEDIR)/cuts.inc
+
+#
+# Compilation of the module ptj_bias
+#
+PY8_CKKW.o: PY8_CKKW.f $(INCLUDES) 
+	$(FC) $(FFLAGS) $(LDFLAGS) -c -o PY8_CKKW.o PY8_CKKW.f
+
+#
+# Generate the copy of the relevant cards in the module C++ headers so that they are compiled
+# directly in the executable and do not have to be read at run time.
+#
+param_card.h: ../../../Cards/param_card.dat
+	python -c "out=open('param_card.h','w'); out.write('const std::string ParamCard = \"%s\";'%'\\\n\"\n\"'.join([line for line in open('../../../Cards/param_card.dat').read().replace('\"',r'\"').split('\n')])); out.close();"
+
+run_card.h: ../../../Cards/run_card.dat
+	python -c "out=open('run_card.h','w'); out.write('const std::string RunCard = \"%s\";'%'\\\n\"\n\"'.join([line for line in open('../../../Cards/run_card.dat').read().replace('\"',r'\"').split('\n')])); out.close();"
+
+proc_card.h: ../../../Cards/proc_card_mg5.dat
+	python -c "out=open('proc_card.h','w'); out.write('const std::string ProcCard = \"%s\";'%'\\\n\"\n\"'.join([line for line in open('../../../Cards/proc_card_mg5.dat').read().replace('\"',r'\"').split('\n')])); out.close();"
+
+py8_mg5amc_bias_interface.o: py8_mg5amc_bias_interface.cc param_card.h run_card.h proc_card.h
+	$(DEFAULT_CPP_COMPILER) $(CXXFLAGS) $(PY8CXXFLAGS) $(LDFLAGS) $(PY8LIBS) -c -o py8_mg5amc_bias_interface.o py8_mg5amc_bias_interface.cc
+
+PY8_CKKW: PY8_CKKW.o py8_mg5amc_bias_interface.o
+	$(call CREATELIB, $(BIASLIBDIR)$(BIASLIBRARY), $^)
+
+#
+# List of the requirements for this module.
+# 'VALID' is the keyword that *must* be returned if everything is in order.
+#
+requirements:	
+ifeq ($(shell $(call CHECK_MG5AMC_VERSION,2.4.2)),True)
+ifeq ($(PYTHIA8_PATH),NotInstalled)
+	@echo "Error: This module requires Pythia8. Consider installing it with the command 'install pythia8'."
+else
+	@echo "VALID"
+endif
+else
+	@echo "Error:: MG5aMC is not recent enough (found "$(MG5AMC_VERSION)")"
+	@echo "FAIL"
+endif

=== added file 'Template/LO/Source/BIAS/PY8_CKKW/py8_mg5amc_bias_interface.cc'
--- Template/LO/Source/BIAS/PY8_CKKW/py8_mg5amc_bias_interface.cc	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/PY8_CKKW/py8_mg5amc_bias_interface.cc	2016-08-23 01:19:18 +0000
@@ -0,0 +1,121 @@
+#include "Pythia8/Pythia.h"
+#include <stdlib.h>
+#include <math.h>
+#include <vector>
+#include <sstream>
+#include "param_card.h"
+#include "run_card.h"
+#include "proc_card.h"
+
+bool initialization_done = false;
+
+Pythia8::Pythia pythia;
+// Stringstream of the events in which MG5aMC will write and PY8 read.
+std::stringstream EvtSS;
+
+void generate_header_SS(std::stringstream & headerSS,
+						const int &      Pythia8BeamA,
+						const int &      Pythia8BeamB,
+						const double &   PythiaBeamEnergyA,
+						const double &   PythiaBeamEnergyB) {
+	headerSS<<"<LesHouchesEvents version=\"3.0\">\n";
+	headerSS<<"<header>\n";
+	// Writing proc_card.dat in the banner
+	headerSS<<"<MG5ProcCard>\n";
+	headerSS<<ProcCard<<"\n";	
+	headerSS<<"</MG5ProcCard>\n";
+	// Writing run_card.dat in the banner	
+	headerSS<<"<MGRunCard>\n";
+	headerSS<<RunCard<<"\n";	
+	headerSS<<"</MGRunCard>\n";
+	// Writing param_card.dat in the banner	
+	headerSS<<"<slha>\n";
+	headerSS<<ParamCard<<"\n";	
+	headerSS<<"</slha>\n";
+	headerSS<<"</header>\n";
+	// Writing the init block (only beam related information is relevant here)	
+	headerSS<<"<init>\n";
+	headerSS<<Pythia8BeamA<<" "<<Pythia8BeamB<<" ";
+	headerSS<<std::scientific<<PythiaBeamEnergyA<<" "<<std::scientific<<PythiaBeamEnergyB<<" ";
+	// We fill in the rest of the init information with blanks
+	headerSS<<"0 0 -1 -1 -3 1\n";
+	headerSS<<"-1.0 0.0 1.0 1\n";
+	headerSS<<"<generator name='MadGraph5_aMC@NLO'>please cite 1405.0301</generator>\n";
+	headerSS<<"</init>\n";
+//	headerSS<<"</LesHouchesEvents>\n";
+	return; 
+}
+
+extern "C" {
+//Fortran interface
+  void py8_bias_weight_(const double &   eCM,
+						const int &      Pythia8BeamA,
+						const double &   PythiaBeamEnergyA,
+						const int &      Pythia8BeamB,
+						const double &   PythiaBeamEnergyB,
+		                const char *     Pythia8EvtRecord,
+		                const double *   p, 
+						const int &      nParticles,
+						const double &   MurScale,
+						const double &   AlphaQCD,
+						const double &   AlphaQED,
+						const int *      Pythia8ID,
+						const int *      Pythia8MotherOne,
+						const int *      Pythia8MotherTwo,
+						const int *      Pythia8ColorOne,
+						const int *      Pythia8ColorTwo,
+						const int *      Pythia8Status,
+						const int *      Pythia8Helicities,
+						double &         OutputBiasWeight    ) 
+  {
+	if (!initialization_done) {
+//     /!\ WARNING this example is tailored for p p > z > e+ e- j for now.
+
+		pythia.readString("Merging:doKtMerging=on");
+		pythia.readString("Merging:Process=pp>e+e-");
+		pythia.readString("Merging:TMS=30.0");
+		pythia.readString("Merging:nJetMax=2");
+		pythia.readString("Merging:applyVeto=off");
+		pythia.readString("Beams:readLHEFheaders=on");
+		pythia.readString("Merging:includeWeightInXsection=off");
+
+		pythia.readString("hadronlevel:all=off");
+		pythia.readString("partonlevel:all=off");
+		pythia.readString("SpaceShower:QEDshowerByL=off");
+		pythia.readString("SpaceShower:QEDshowerByQ=off");
+		pythia.readString("ProcessLevel:resonanceDecays=off");
+		pythia.readString("BeamRemnants:primordialKT=off");
+		pythia.readString("TimeShower:QEDshowerByQ=off");
+		pythia.readString("TimeShower:QEDshowerByL=off");
+		pythia.readString("partonlevel:mpi=off");
+		pythia.readString("PartonLevel:FSRinResonances=off");
+		pythia.readString("PartonLevel:Remnants=off");
+		pythia.readString("Check:event=off");
+		
+		pythia.settings.mode("Beams:frameType",4);		
+		std::stringstream headerSS;
+		generate_header_SS(headerSS, Pythia8BeamA,Pythia8BeamB,PythiaBeamEnergyA,PythiaBeamEnergyB);
+		EvtSS.str(headerSS.str());
+		pythia.init(&EvtSS, &headerSS);
+		initialization_done = true;
+	}
+
+	std::string EvtStr(Pythia8EvtRecord);
+	EvtStr = EvtStr.substr(0,EvtStr.find("</event>")+8);
+	// Restart the buffer, clearing EOF state
+	EvtSS.clear();
+	// Set the stringstream buffer to contain the event only
+	EvtSS.str(EvtStr);
+
+//	std::cout<<"---> Event transmitted to Pythiaa8 <---\n"<<EvtSS.str()<<std::endl;
+
+	// Compute the merging weight. All other shower processsing has been disabled in the initialisation above.
+	pythia.next();	
+//	pythia.event.list();
+    // This is the CKKW Sudakov merging weight from the trial shower.
+	std::cout<<pythia.info.mergingWeight()<<std::endl;
+//	exit(0);
+
+	OutputBiasWeight = pythia.info.mergingWeight();
+  }
+}

=== added directory 'Template/LO/Source/BIAS/dummy'
=== added file 'Template/LO/Source/BIAS/dummy/dummy.f'
--- Template/LO/Source/BIAS/dummy/dummy.f	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/dummy/dummy.f	2016-08-23 01:19:18 +0000
@@ -0,0 +1,44 @@
+C ************************************************************
+C Source for the library implementing a dummt bias function 
+C always returns one
+C ************************************************************
+      
+      subroutine bias_wgt(p, original_weight, bias_weight)
+          implicit none
+C
+C Parameters
+C
+          include '../../nexternal.inc'
+C
+C Arguments
+C
+          double precision p(0:3,nexternal)
+          double precision original_weight, bias_weight
+C
+C local variables
+C
+C
+C Global variables
+C
+C Mandatory common block to be defined in bias modules
+C
+          double precision stored_bias_weight
+          data stored_bias_weight/1.0d0/          
+          logical impact_xsec, requires_full_event_info
+C         Not impacting the xsec since the bias is 1.0. Therefore
+C         bias_wgt will not be written in the lhe event file.
+C         Setting it to .True. makes sure that it will not be written.
+          data impact_xsec/.True./
+C         Of course this module does not require the full event
+C         information (color, resonances, helicities, etc..)
+          data requires_full_event_info/.False./ 
+          common/bias/stored_bias_weight,impact_xsec,
+     &                requires_full_event_info
+
+C --------------------
+C BEGIN IMPLEMENTATION
+C --------------------
+
+          bias_weight = 1.0d0
+
+      end subroutine bias_wgt

=== added file 'Template/LO/Source/BIAS/dummy/makefile'
--- Template/LO/Source/BIAS/dummy/makefile	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/dummy/makefile	2016-08-23 01:19:18 +0000
@@ -0,0 +1,21 @@
+
+include ../../make_opts
+
+all: dummy 
+
+clean:
+	$(RM) *.o $(BIASLIBDIR)$(BIASLIBRARY)
+
+#
+# Compilation of the module dummy
+#
+
+dummy: dummy.o
+	$(call CREATELIB, $(BIASLIBDIR)$(BIASLIBRARY), $^)
+
+#
+# List of the requirements for this module.
+# 'VALID' is the keyword that *must* be returned if everything is in order.
+#
+requirements:	
+	@echo "VALID"

=== added directory 'Template/LO/Source/BIAS/ptj_bias'
=== added file 'Template/LO/Source/BIAS/ptj_bias/makefile'
--- Template/LO/Source/BIAS/ptj_bias/makefile	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/ptj_bias/makefile	2016-08-23 01:19:18 +0000
@@ -0,0 +1,27 @@
+include ../../make_opts
+
+all: ptj_bias
+
+clean:
+	$(RM) *.o $(BIASLIBDIR)$(BIASLIBRARY)
+
+#
+# Compilation of the module ptj_bias
+#
+ptj_bias.o: ptj_bias.f ../bias.inc
+	$(FC) $(FFLAGS) $(LDFLAGS) -c -o ptj_bias.o ptj_bias.f
+
+ptj_bias: ptj_bias.o
+	$(call CREATELIB, $(BIASLIBDIR)$(BIASLIBRARY), $^)
+
+#
+# List of the requirements for this module.
+# 'VALID' is the keyword that *must* be returned if everything is in order.
+#
+requirements:	
+ifeq ($(shell $(call CHECK_MG5AMC_VERSION,2.4.2)),True)
+	@echo "VALID"
+else
+	@echo "Error:: MG5aMC is not recent enough (found "$(MG5AMC_VERSION)")"
+	@echo "FAIL"
+endif

=== added file 'Template/LO/Source/BIAS/ptj_bias/ptj_bias.f'
--- Template/LO/Source/BIAS/ptj_bias/ptj_bias.f	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/ptj_bias/ptj_bias.f	2016-08-23 01:19:18 +0000
@@ -0,0 +1,101 @@
+C ************************************************************
+C Source for the library implementing a bias function that 
+C populates the large pt tale of the leading jet. 
+C
+C The two options of this subroutine, that can be set in
+C the run card are:
+C    > (double precision) ptj_bias_target_ptj : target ptj value
+C    > (double precision) ptj_bias_enhancement_power : exponent
+C
+C Schematically, the functional form of the enhancement is
+C    bias_wgt = [ptj(evt)/mean_ptj]^enhancement_power
+C ************************************************************
+C
+C The following lines are read by MG5aMC to set what are the 
+C relevant parameters for this bias module.
+C
+C  parameters = {'ptj_bias_target_ptj': 1000.0,
+C               'ptj_bias_enhancement_power': 4.0}
+C
+
+      subroutine bias_wgt(p, original_weight, bias_weight)
+          implicit none
+C
+C Parameters
+C
+          include '../../maxparticles.inc'          
+          include '../../nexternal.inc'
+
+C
+C Arguments
+C
+          double precision p(0:3,nexternal)
+          double precision original_weight, bias_weight
+C
+C local variables
+C
+          integer i
+          double precision ptj(nexternal)
+          double precision max_ptj
+c
+c local variables defined in the run_card
+c
+          double precision ptj_bias_target_ptj
+          double precision ptj_bias_enhancement_power
+C
+C Global variables
+C
+C
+C Mandatory common block to be defined in bias modules
+C
+          double precision stored_bias_weight
+          data stored_bias_weight/1.0d0/          
+          logical impact_xsec, requires_full_event_info
+C         We only want to bias distributions, but not impact the xsec. 
+          data impact_xsec/.False./
+C         Of course this module does not require the full event
+C         information (color, resonances, helicities, etc..)
+          data requires_full_event_info/.False./ 
+          common/bias/stored_bias_weight,impact_xsec,
+     &                requires_full_event_info
+C
+C Accessingt the details of the event
+C
+          logical is_a_j(nexternal),is_a_l(nexternal),
+     &            is_a_b(nexternal),is_a_a(nexternal),
+     &            is_a_onium(nexternal),is_a_nu(nexternal),
+     &            is_heavy(nexternal),do_cuts(nexternal)
+          common/to_specisa/is_a_j,is_a_a,is_a_l,is_a_b,is_a_nu,
+     &                      is_heavy,is_a_onium,do_cuts
+
+C
+C    Setup the value of the parameters from the run_card    
+C
+      include '../bias.inc'
+
+C --------------------
+C BEGIN IMPLEMENTATION
+C --------------------
+          
+          do i=1,nexternal
+            ptj(i)=-1.0d0
+            if (is_a_j(i)) then
+              ptj(i)=sqrt(p(1,i)**2+p(2,i)**2)
+            endif
+          enddo
+
+          max_ptj=-1.0d0
+          do i=1,nexternal
+            max_ptj = max(max_ptj,ptj(i))
+          enddo
+          if (max_ptj.lt.0.0d0) then
+            bias_weight = 1.0d0
+            return
+          endif
+
+          bias_weight = (max_ptj/ptj_bias_target_ptj)
+     &                                      **ptj_bias_enhancement_power
+
+          return
+
+      end subroutine bias_wgt

=== modified file 'Template/LO/Source/cuts.inc'
--- Template/LO/Source/cuts.inc	2015-04-09 01:31:09 +0000
+++ Template/LO/Source/cuts.inc	2016-08-23 01:19:18 +0000
@@ -88,6 +88,10 @@
       REAL*8 KT_DURHAM, D_PARAMETER
       LOGICAL DO_KT_DURHAM
       COMMON /JET_MEASURE_CUTS/ KT_DURHAM, D_PARAMETER
-      
-      
-      
\ No newline at end of file
+
+C COMMON BLOCK FOR PYTHIA JET MEASURE CUT
+      REAL*8 PT_LUND
+      COMMON /PYTHIA_MEASURE_CUTS/ PT_LUND
+      
+      
+      

=== added file 'Template/LO/Source/lhe_event_infos.inc'
--- Template/LO/Source/lhe_event_infos.inc	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/lhe_event_infos.inc	2016-08-23 01:19:18 +0000
@@ -0,0 +1,16 @@
+      integer jpart(7,-nexternal+3:2*nexternal-3)
+      double precision pb(0:4,-nexternal+3:2*nexternal-3)
+      integer isym(nexternal,99),jsym, npart
+      double precision sscale,aaqcd,aaqed
+      character*1000 buff
+      character*(s_bufflen) s_buff(7)
+      integer nclus
+      character*(clus_bufflen) buffclus(nexternal)
+	  character*(maxEventLength) event_record
+      logical AlreadySetInBiasModule
+
+	  common/to_lhe_event_info/jpart,pb,s_buff,buff,nclus,buffclus,event_record,
+     &  sscale,aaqcd,aaqed,isym,jsym,npart,AlreadySetInBiasModule
+
+      integer ngroup
+      common/to_group/ngroup

=== modified file 'Template/LO/Source/make_opts'
--- Template/LO/Source/make_opts	2016-07-01 22:45:47 +0000
+++ Template/LO/Source/make_opts	2016-08-23 01:19:18 +0000
@@ -1,12 +1,18 @@
+# Defines the variable pertaining to MG5aMC environment
+PYTHIA8_PATH=NotInstalled
+MG5AMC_VERSION=SpecifiedByMG5aMCAtRunTime
 DEFAULT_F2PY_COMPILER=f2py
 DEFAULT_F_COMPILER=gfortran
-MACFLAG=-mmacosx-version-min=10.7
-DEFAULT_CPP_COMPILER=clang
-STDLIB=-lc++
-STDLIB_FLAG=-stdlib=libc++
+MACFLAG=
+DEFAULT_CPP_COMPILER=g++
+STDLIB=-lstdc++
+STDLIB_FLAG=
 #end_of_make_opts_variables
+
+BIASLIBDIR=../../../lib/
+BIASLIBRARY=libbias.$(libext)
+
 # Rest of the makefile
-
 ifeq ($(origin FFLAGS),undefined)
 FFLAGS= -O -w -fbounds-check -fPIC
 #FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
@@ -90,4 +96,9 @@
 else
 alfas_functions=alfas_functions
 llhapdf=
-endif
\ No newline at end of file
+endif
+
+# Helper function to check MG5 version
+define CHECK_MG5AMC_VERSION
+python -c 'import re; from distutils.version import StrictVersion; print StrictVersion("$(MG5AMC_VERSION)") >= StrictVersion("$(1)") if re.match("^[\d\.]+$$","$(MG5AMC_VERSION)") else True;'
+endef

=== modified file 'Template/LO/Source/run.inc'
--- Template/LO/Source/run.inc	2015-04-02 22:56:24 +0000
+++ Template/LO/Source/run.inc	2016-08-23 01:19:18 +0000
@@ -63,3 +63,11 @@
 C
       logical MC_grouped_subproc
       common/to_MC_grouped_subproc/MC_grouped_subproc
+
+C
+C Controls what are the PDGs included in the CKKWl merging procedure, i.e. what
+C are the PDGs subject to the ktdurham cut
+C
+      integer pdgs_for_merging_cut(0:1000)
+      common/TO_MERGE/pdgs_for_merging_cut
+

=== modified file 'Template/LO/Source/rw_events.f'
--- Template/LO/Source/rw_events.f	2014-10-07 13:31:57 +0000
+++ Template/LO/Source/rw_events.f	2016-08-23 01:19:18 +0000
@@ -30,7 +30,6 @@
       character*(s_bufflen) s_buff(*)
       integer nclus
       character*(clus_bufflen) buffclus(*)
-
 c
 c     Local
 c
@@ -46,6 +45,10 @@
 
       data lun_ban/37/
       data banner_open/.false./
+
+      double precision bias_weight
+      logical impact_xsec
+      common/bias/bias_weight,impact_xsec
 c-----
 c  Begin Code
 c-----     
@@ -80,6 +83,21 @@
          buff=''
         endif
       endif
+
+c     Reading the bias weight (if present)
+      read(lun,'(a300)',end=99,err=99) buftmp
+      if(buftmp(1:6).ne.'<rwgt>') then
+         backspace(lun)
+         bias_weight = 1.0d0
+      else
+         do while(buftmp(1:7).ne.'</rwgt>')
+           read(lun,'(a300)',end=99,err=99) buftmp
+           if (buftmp(1:16).eq."<wgt id='bias'> ") then
+              read(buftmp(17:31),'(1e15.7)') bias_weight
+           endif
+         enddo
+      endif
+
 c     Systematics info
       read(lun,'(a)',end=99,err=99) s_buff(1)
       if(s_buff(1).ne.'<mgrwt>') then
@@ -114,9 +132,103 @@
  55   format(i3,5e19.11)         
       end
 
+      subroutine write_event_to_stream(evt_record,P,wgt,nexternal,ic,
+     &     ievent,scale,aqcd, aqed,buff,u_syst,s_buff,nclus,buffclus)
+c********************************************************************
+C     This an *exact* copy of write_event, except that it writes it 
+C     to a character array argument as opposed to an I/O stream.
+c********************************************************************
+      implicit none
+
+      include 'maxparticles.inc'
+      include 'run_config.inc'
+c
+c     parameters
+c
+      double precision pi
+      parameter (pi = 3.1415926d0)
+c
+c     Arguments
+c      
+      character*(maxEventLength) evt_record
+      integer ievent
+      integer nexternal, ic(7,*)
+      double precision P(0:4,*),wgt
+      double precision aqcd, aqed, scale
+      character*1000 buff
+      logical u_syst
+      character*(s_bufflen) s_buff(*)
+      integer nclus
+      character*(clus_bufflen) buffclus(*)
+c
+c     Local
+c
+      integer i,j,k
+      character*(maxEventLength) largeBuff
+c
+c     Global
+c
+      double precision bias_weight
+      logical impact_xsec
+      common/bias/bias_weight,impact_xsec
+
+c-----
+c  Begin Code
+c-----     
+c      aqed= gal(1)*gal(1)/4d0/pi
+c      aqcd = g*g/4d0/pi
+      write(largeBuff,'(a)') '<event>'
+      evt_record=trim(evt_record)//trim(largeBuff)
+      write(largeBuff,'(i2,i4,4e15.7)') nexternal,ievent,wgt,scale,
+     $                                   aqed,aqcd
+      evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      do i=1,nexternal
+         write(largeBuff,51) ic(1,i),ic(6,i),(ic(j,i),j=2,5),
+     $     (p(j,i),j=1,3),p(0,i),p(4,i),0.,real(ic(7,i))
+         evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      enddo
+      if(buff(1:7).eq.'<scales') then
+        write(largeBuff,'(a)') buff(1:len_trim(buff))
+        evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      endif
+      if(buff(1:1).eq.'#') then
+        write(largeBuff,'(a)') buff(1:len_trim(buff))
+        evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      endif
+      if(.not.impact_xsec) then
+          write(largeBuff,'(a)') '<rwgt>'
+          evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+          write(largeBuff,'(a16,1e15.7,a6)') "<wgt id='bias'> ",
+     $                                              bias_weight,"</wgt>"
+          evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+          write(largeBuff,'(a)') '</rwgt>'
+          evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      endif 
+      if(u_syst)then
+         do i=1,7
+            write(largeBuff,'(a)') s_buff(i)(1:len_trim(s_buff(i)))
+            evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+         enddo
+      endif
+      do i=1,nclus
+         write(largeBuff,'(a)') buffclus(i)(1:len_trim(buffclus(i)))
+         evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      enddo
+      write(largeBuff,'(a)') '</event>'
+      evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      return
+ 51   format(i9,5i5,5e19.11,f3.0,f4.0)
+      end
+
+
       subroutine write_event(lun,P,wgt,nexternal,ic,ievent,scale,aqcd,
      $     aqed,buff,u_syst,s_buff,nclus,buffclus)
 c********************************************************************
+c
+c /!\ When making changes to this subroutine, make sure to accordingly
+c     update write_event_to_stream
+c
+c********************************************************************
 c     Writes one event from data file #lun according to LesHouches
 c     ic(1,*) = Particle ID
 c     ic(2.*) = Mothup(1)
@@ -154,6 +266,9 @@
 c
 c     Global
 c
+      double precision bias_weight
+      logical impact_xsec
+      common/bias/bias_weight,impact_xsec
 
 c-----
 c  Begin Code
@@ -168,7 +283,13 @@
      $     (p(j,i),j=1,3),p(0,i),p(4,i),0.,real(ic(7,i))
       enddo
       if(buff(1:7).eq.'<scales') write(lun,'(a)') buff(1:len_trim(buff))
-      if(buff(1:1).eq.'#') write(lun,'(a)') buff(1:len_trim(buff))
+      if(buff(1:1).eq.'#') write(lun,'(a)') buff(1:len_trim(buff))      
+      if(.not.impact_xsec) then
+          write(lun,'(a)') '<rwgt>'
+          write(lun,'(a16,1e15.7,a6)') "<wgt id='bias'> ",bias_weight,
+     $                                                          "</wgt>"
+          write(lun,'(a)') '</rwgt>'
+      endif 
       if(u_syst)then
          do i=1,7
             write(lun,'(a)') s_buff(i)(1:len_trim(s_buff(i)))

=== modified file 'Template/LO/SubProcesses/cuts.f'
--- Template/LO/SubProcesses/cuts.f	2016-06-07 09:33:55 +0000
+++ Template/LO/SubProcesses/cuts.f	2016-08-23 01:19:18 +0000
@@ -1,1319 +1,1750 @@
-      logical function pass_point(p)
-c************************************************************************
-c     This function is called from sample to see if it needs to 
-c     bother calculating the weight from all the different conficurations
-c     You can either just return true, or have it call passcuts
-c************************************************************************
-      implicit none
-c
-c     Arguments
-c
-      double precision p
-c
-c     External
-c
-      logical passcuts
-      external passcuts
-c-----
-c  Begin Code
-c-----
-      pass_point = .true.
-c      pass_point = passcuts(p)
-      end
-C 
-      LOGICAL FUNCTION PASSCUTS(P)
-C**************************************************************************
-C     INPUT:
-C            P(0:3,1)           MOMENTUM OF INCOMING PARTON
-C            P(0:3,2)           MOMENTUM OF INCOMING PARTON
-C            P(0:3,3)           MOMENTUM OF ...
-C            ALL MOMENTA ARE IN THE REST FRAME!!
-C            COMMON/JETCUTS/   CUTS ON JETS
-C     OUTPUT:
-C            TRUE IF EVENTS PASSES ALL CUTS LISTED
-C**************************************************************************
-      IMPLICIT NONE
-c
-c     Constants
-c
-      include 'genps.inc'
-      include 'nexternal.inc'
-C
-C     ARGUMENTS
-C
-      REAL*8 P(0:3,nexternal)
-
-C
-C     LOCAL
-C
-      LOGICAL FIRSTTIME,FIRSTTIME2,pass_bw,notgood,good,foundheavy
-      LOGICAL DEBUG
-      integer i,j,njets,nheavyjets,nleptons,hardj1,hardj2
-      REAL*8 XVAR,ptmax1,ptmax2,htj,tmp,inclht
-      real*8 ptemp(0:3), ptemp2(0:3)
-      character*20 formstr
-C
-C     PARAMETERS
-C
-      real*8 PI
-      parameter( PI = 3.14159265358979323846d0 )
-C
-C     EXTERNAL
-C
-      REAL*8 R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,PtDot
-      logical cut_bw,setclscales
-      external R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,cut_bw,setclscales,PtDot
-C
-C     GLOBAL
-C
-      include 'run.inc'
-      include 'cuts.inc'
-
-      double precision ptjet(nexternal)
-      double precision ptheavyjet(nexternal)
-      double precision ptlepton(nexternal)
-      double precision temp
-
-C VARIABLES TO SPECIFY JETS
-      DOUBLE PRECISION PJET(NEXTERNAL,0:3)
-      DOUBLE PRECISION PTMIN
-      DOUBLE PRECISION PT1,PT2
-      INTEGER K,J1,J2
-
-C VARIABLES FOR KT CUT
-      DOUBLE PRECISION PTNOW,COSTH,PABS1,PABS2
-      DOUBLE PRECISION ETA1,ETA2,COSH_DETA,COS_DPHI,KT1SQ,KT2SQ, DPHI
-
-      double precision etmin(nincoming+1:nexternal),etamax(nincoming+1:nexternal)
-      double precision emin(nincoming+1:nexternal)
-      double precision                    r2min(nincoming+1:nexternal,nincoming+1:nexternal)
-      double precision s_min(nexternal,nexternal)
-      double precision etmax(nincoming+1:nexternal),etamin(nincoming+1:nexternal)
-      double precision emax(nincoming+1:nexternal)
-      double precision r2max(nincoming+1:nexternal,nincoming+1:nexternal)
-      double precision s_max(nexternal,nexternal)
-      double precision ptll_min(nexternal,nexternal),ptll_max(nexternal,nexternal)
-      double precision inclHtmin,inclHtmax
-      common/to_cuts/  etmin, emin, etamax, r2min, s_min,
-     $     etmax, emax, etamin, r2max, s_max, ptll_min, ptll_max, inclHtmin,inclHtmax
-
-      double precision ptjmin4(4),ptjmax4(4),htjmin4(2:4),htjmax4(2:4)
-      logical jetor
-      common/to_jet_cuts/ ptjmin4,ptjmax4,htjmin4,htjmax4,jetor
-
-      double precision ptlmin4(4),ptlmax4(4)
-      common/to_lepton_cuts/ ptlmin4,ptlmax4
-
-c
-c     Special cuts
-c
-
-      integer        lbw(0:nexternal)  !Use of B.W.
-      common /to_BW/ lbw
-C
-C     SPECIAL CUTS
-C
-      LOGICAL  IS_A_J(NEXTERNAL),IS_A_L(NEXTERNAL)
-      LOGICAL  IS_A_B(NEXTERNAL),IS_A_A(NEXTERNAL),IS_A_ONIUM(NEXTERNAL)
-      LOGICAL  IS_A_NU(NEXTERNAL),IS_HEAVY(NEXTERNAL)
-      logical  do_cuts(nexternal)
-      COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,
-     . IS_A_ONIUM, do_cuts
-C
-C     Keep track of whether cuts already calculated for this event
-C
-      LOGICAL CUTSDONE,CUTSPASSED
-      COMMON/TO_CUTSDONE/CUTSDONE,CUTSPASSED
-      DATA CUTSDONE,CUTSPASSED/.FALSE.,.FALSE./
-
-C $B$ MW_NEW_DEF $E$ !this is a tag for MadWeight
-
-      double precision xqcutij(nexternal,nexternal),xqcuti(nexternal)
-      common/to_xqcuts/xqcutij,xqcuti
-
-c jet cluster algorithm
-      integer nQCD !,NJET,JET(nexternal)
-c      double precision plab(0:3, nexternal)
-      double precision pQCD(0:3,nexternal)!,PJET(0:3,nexternal)
-c      double precision rfj,sycut,palg,fastjetdmerge
-c      integer njet_eta
-c     Photon isolation
-      integer nph,nem,nin
-      double precision ptg,chi_gamma_iso,iso_getdrv40
-      double precision Etsum(0:nexternal)
-      real drlist(nexternal)
-      double precision pgamma(0:3,nexternal),pem(0:3,nexternal)
-      logical alliso
-C     Sort array of results: ismode>0 for real, isway=0 for ascending order
-      integer ismode,isway,izero,isorted(nexternal)
-      parameter (ismode=1)
-      parameter (isway=0)
-      parameter (izero=0)
-
-      include 'coupl.inc'
-C
-C
-c
-      DATA FIRSTTIME,FIRSTTIME2/.TRUE.,.TRUE./
-
-c put momenta in common block for couplings.f
-      double precision pp(0:3,max_particles)
-      common /momenta_pp/pp
-
-      DATA DEBUG/.FALSE./
-
-C-----
-C  BEGIN CODE
-C-----
-
-
-
-      PASSCUTS=.TRUE.             !EVENT IS OK UNLESS OTHERWISE CHANGED
-      IF (FIRSTTIME) THEN
-         FIRSTTIME=.FALSE.
-c      Preparation for reweighting by setting up clustering by diagrams
-         call initcluster()
-c
-c
-         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'i8)'
-         write(*,formstr) 'Particle',(i,i=nincoming+1,nexternal)
-         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'f8.1)'
-         write(*,formstr) 'Et >',(etmin(i),i=nincoming+1,nexternal)
-         write(*,formstr) 'E >',(emin(i),i=nincoming+1,nexternal)
-         write(*,formstr) 'Eta <',(etamax(i),i=nincoming+1,nexternal)
-         write(*,formstr) 'xqcut: ',(xqcuti(i),i=nincoming+1,nexternal)
-         write(formstr,'(a,i2.2,a)')'(a,i2,a,',nexternal,'f8.1)'
-         do j=nincoming+1,nexternal-1
-            write(*,formstr) 'd R #',j,'  >',(-0.0,i=nincoming+1,j),
-     &           (r2min(i,j),i=j+1,nexternal)
-            do i=j+1,nexternal
-               r2min(i,j)=r2min(i,j)*dabs(r2min(i,j))    !Since r2 returns distance squared
-               r2max(i,j)=r2max(i,j)*dabs(r2max(i,j))
-            enddo
-         enddo
-         do j=nincoming+1,nexternal-1
-            write(*,formstr) 's min #',j,'>',
-     &           (s_min(i,j),i=nincoming+1,nexternal)
-         enddo
-         do j=nincoming+1,nexternal-1
-            write(*,formstr) 'xqcutij #',j,'>',
-     &           (xqcutij(i,j),i=nincoming+1,nexternal)
-         enddo
-
-cc
-cc     Set the strong coupling
-cc
-c         call set_ren_scale(P,scale)
-c
-cc     Check that the user funtions for setting the scales
-cc     have been edited if the choice of an event-by-event
-cc     scale choice has been made 
-c
-c         if(.not.fixed_ren_scale) then
-c            if(scale.eq.0d0) then
-c               write(6,*) 
-c               write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
-c               write(6,*) ' Dynamical renormalization scale choice '
-c               write(6,*) ' selected but user subroutine' 
-c               write(6,*) ' set_ren_scale not edited in file:setpara.f'
-c               write(6,*) ' Switching to a fixed_ren_scale choice'
-c               write(6,*) ' with scale=zmass'
-c               scale=91.2d0
-c               write(6,*) 'scale=',scale
-c               fixed_ren_scale=.true.
-c               call set_ren_scale(P,scale)
-c            endif
-c         endif
-         
-c         if(.not.fixed_fac_scale) then
-c            call set_fac_scale(P,q2fact)
-c            if(q2fact(1).eq.0d0.or.q2fact(2).eq.0d0) then
-c               write(6,*) 
-c               write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
-c               write(6,*) ' Dynamical renormalization scale choice '
-c               write(6,*) ' selected but user subroutine' 
-c               write(6,*) ' set_fac_scale not edited in file:setpara.f'
-c               write(6,*) ' Switching to a fixed_fac_scale choice'
-c               write(6,*) ' with q2fact(i)=zmass**2'
-c               fixed_fac_scale=.true.
-c               q2fact(1)=91.2d0**2
-c               q2fact(2)=91.2d0**2
-c               write(6,*) 'scales=',q2fact(1),q2fact(2)
-c            endif
-c         endif
-
-         if(fixed_ren_scale) then
-            G = SQRT(4d0*PI*ALPHAS(scale))
-            call update_as_param()
-         endif
-
-c     Put momenta in the common block to zero to start
-         do i=0,3
-            do j=1,max_particles
-               pp(i,j) = 0d0
-            enddo
-         enddo
-         
-      ENDIF ! IF FIRSTTIME
-
-      IF (CUTSDONE) THEN
-         PASSCUTS=CUTSPASSED
-         RETURN
-      ENDIF
-      CUTSDONE=.TRUE.
-c      CUTSPASSED=.FALSE.
-
-c
-c     Make sure have reasonable 4-momenta
-c
-      if (p(0,1) .le. 0d0) then
-         passcuts=.false.
-         return
-      endif
-
-c     Also make sure there's no INF or NAN
-      do i=1,nexternal
-         do j=0,3
-            if(p(j,i).gt.1d32.or.p(j,i).ne.p(j,i))then
-               passcuts=.false.
-               return
-            endif
-         enddo
-      enddo
-      
-c
-c     Limit S_hat
-c
-c      if (x1*x2*stot .gt. 500**2) then
-c         passcuts=.false.
-c         return
-c      endif
-
-C $B$ DESACTIVATE_CUT $E$ !This is a tag for MadWeight
-
-      if(debug) write (*,*) '============================='
-      if(debug) write (*,*) ' EVENT STARTS TO BE CHECKED  '
-      if(debug) write (*,*) '============================='
-c     
-c     p_t min & max cuts
-c     
-      do i=nincoming+1,nexternal
-         if(debug) write (*,*) 'pt(',i,')=',pt(p(0,i)),'   ',etmin(i),
-     $        ':',etmax(i)
-         notgood=(pt(p(0,i)) .lt. etmin(i)).or.
-     &        (etmax(i).ge.0d0.and.pt(p(0,i)) .gt. etmax(i))
-         if (notgood) then
-            if(debug) write (*,*) i,' -> fails'
-            passcuts=.false.
-            return
-         endif
-      enddo
-c
-c    missing ET min & max cut + Invariant mass of leptons and neutrino 
-c    nb: missing Et defined as the vector sum over the neutrino's pt
-c
-c-- reset ptemp(0:3)
-      do j=0,3
-         ptemp(j)=0 ! for the neutrino
-         ptemp2(j)=0 ! for the leptons
-      enddo
-c-  sum over the momenta
-      do i=nincoming+1,nexternal
-         if(is_a_nu(i)) then            
-         if(debug) write (*,*) i,' -> neutrino '
-            do j=0,3
-               ptemp(j)=ptemp(j)+p(j,i)
-            enddo
-         elseif(is_a_l(i)) then            
-         if(debug) write (*,*) i,' -> lepton '
-            do j=0,3
-               ptemp2(j)=ptemp2(j)+p(j,i)
-            enddo
-         endif
-
-      enddo
-c-  check the et
-      if(debug.and.ptemp(0).eq.0d0) write (*,*) 'No et miss in event'
-      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Et miss =',pt(ptemp(0)),'   ',misset,':',missetmax
-      if(debug.and.ptemp2(0).eq.0d0) write (*,*) 'No leptons in event'
-      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Energy of leptons =',pt(ptemp2(0))
-      if(ptemp(0).gt.0d0) then
-         notgood=(pt(ptemp(0)) .lt. misset).or.
-     &        (missetmax.ge.0d0.and.pt(ptemp(0)) .gt. missetmax)
-         if (notgood) then
-            if(debug) write (*,*) ' missing et cut -> fails'
-            passcuts=.false.
-            return
-         endif
-      endif
-      if (mmnl.gt.0d0.or.mmnlmax.ge.0d0)then
-         if(dsqrt(SumDot(ptemp,ptemp2,1d0)).lt.mmnl.or.mmnlmax.ge.0d0.and.dsqrt(SumDot(ptemp, ptemp2,1d0)).gt.mmnlmax) then
-            if(debug) write (*,*) 'lepton invariant mass -> fails'
-            passcuts=.false.
-            return
-         endif
-      endif
-c
-c     pt cut on heavy particles
-c     gives min(pt) for (at least) one heavy particle
-c
-      if(ptheavy.gt.0d0)then
-         passcuts=.false.
-         foundheavy=.false.
-         do i=nincoming+1,nexternal
-            if(is_heavy(i)) then            
-               if(debug) write (*,*) i,' -> heavy '
-               foundheavy=.true.
-               if(pt(p(0,i)).gt.ptheavy) passcuts=.true.
-            endif
-         enddo
-         
-         if(.not.passcuts.and.foundheavy)then
-            if(debug) write (*,*) ' heavy particle cut -> fails'
-            return
-         else
-            passcuts=.true.
-         endif
-      endif
-c     
-c     E min & max cuts
-c     
-      do i=nincoming+1,nexternal
-         if(debug) write (*,*) 'p(0,',i,')=',p(0,i),'   ',emin(i),':',emax(i)
-         notgood=(p(0,i) .le. emin(i)).or.
-     &        (emax(i).ge.0d0 .and. p(0,i) .gt. emax(i))
-         if (notgood) then
-            if(debug) write (*,*) i,' -> fails'
-            passcuts=.false.
-            return
-         endif
-      enddo
-c     
-c     Rapidity  min & max cuts
-c     
-      do i=nincoming+1,nexternal
-         if(debug) write (*,*) 'abs(rap(',i,'))=',abs(rap(p(0,i))),'   ',etamin(i),':',etamax(i)
-         notgood=(etamax(i).ge.0.and.abs(rap(p(0,i))) .gt. etamax(i)).or.
-     &        (abs(rap(p(0,i))) .lt. etamin(i))
-         if (notgood) then
-            if(debug) write (*,*) i,' -> fails'
-            passcuts=.false.
-            return
-         endif
-      enddo
-c     
-c     DeltaR min & max cuts
-c     
-      do i=nincoming+1,nexternal-1
-         do j=i+1,nexternal
-            if(debug) write (*,*) 'r2(',i, ',' ,j,')=',dsqrt(r2(p(0,i),p(0,j)))
-            if(debug) write (*,*) dsqrt(r2min(j,i)),dsqrt(r2max(j,i))
-            if(r2min(j,i).gt.0.or.r2max(j,i).ge.0d0) then
-               tmp=r2(p(0,i),p(0,j))
-               notgood=(tmp .lt. r2min(j,i)).or.
-     $              (r2max(j,i).ge.0d0 .and. tmp .gt. r2max(j,i))
-               if (notgood) then
-                  if(debug) write (*,*) i,j,' -> fails'
-                  passcuts=.false.
-                  return
-               endif
-            endif
-         enddo
-      enddo
-
-
-c     s-channel min & max pt of sum of 4-momenta
-c     
-      do i=nincoming+1,nexternal-1
-         do j=i+1,nexternal
-            if(debug)write (*,*) 'ptll(',i,',',j,')=',PtDot(p(0,i),p(0,j))
-            if(debug)write (*,*) ptll_min(j,i),ptll_max(j,i)
-            if(ptll_min(j,i).gt.0.or.ptll_max(j,i).ge.0d0) then
-               tmp=PtDot(p(0,i),p(0,j))
-               notgood=(tmp .lt. ptll_min(j,i).or.
-     $              ptll_max(j,i).ge.0d0 .and. tmp.gt.ptll_max(j,i))
-               if (notgood) then
-                  if(debug) write (*,*) i,j,' -> fails'
-                  passcuts=.false.
-                  return
-               endif
-            endif
-         enddo
-      enddo
-
-
-
-
-c     
-c     s-channel min & max invariant mass cuts
-c     
-      do i=nincoming+1,nexternal-1
-         do j=i+1,nexternal
-            if(debug) write (*,*) 's(',i,',',j,')=',Sumdot(p(0,i),p(0,j),+1d0)
-            if(debug) write (*,*) s_min(j,i),s_max(j,i)
-            if(s_min(j,i).gt.0.or.s_max(j,i).ge.0d0) then
-               tmp=SumDot(p(0,i),p(0,j),+1d0)
-               if(s_min(j,i).le.s_max(j,i) .or. s_max(j,i).lt.0d0)then
-                  notgood=(tmp .lt. s_min(j,i).or.
-     $                 s_max(j,i).ge.0d0 .and. tmp .gt. s_max(j,i)) 
-                  if (notgood) then
-                     if(debug) write (*,*) i,j,' -> fails'
-                     passcuts=.false.
-                     return
-                  endif
-               else
-                  notgood=(tmp .lt. s_min(j,i).and.tmp .gt. s_max(j,i)) 
-                  if (notgood) then
-                     if(debug) write (*,*) i,j,' -> fails'
-                     passcuts=.false.
-                     return
-                  endif
-               endif
-            endif
-         enddo
-      enddo
-C     $B$DESACTIVATE_BW_CUT$B$ This is a Tag for MadWeight
-c     
-c     B.W. phase space cuts
-c     
-      pass_bw=cut_bw(p)
-c     JA 4/8/11 always check pass_bw
-      if ( pass_bw.and..not.CUTSPASSED) then
-         passcuts=.false.
-         if(debug) write (*,*) ' pass_bw -> fails'
-         return
-      endif
-C     $E$DESACTIVATE_BW_CUT$E$ This is a Tag for MadWeight
-        CUTSPASSED=.FALSE.
-C     
-C     maximal and minimal pt of the jets sorted by pt
-c     
-      njets=0
-      nheavyjets=0
-
-c- fill ptjet with the pt's of the jets.
-      do i=nincoming+1,nexternal
-         if(is_a_j(i)) then
-            njets=njets+1
-            ptjet(njets)=pt(p(0,i))
-         endif
-         if(is_a_b(i)) then
-            nheavyjets=nheavyjets+1
-            ptheavyjet(nheavyjets)=pt(p(0,i))
-         endif
-
-      enddo
-      if(debug) write (*,*) 'not yet ordered ',njets,'   ',ptjet
-
-C----------------------------------------------------------------------------
-C     DURHAM_KT CUT
-C----------------------------------------------------------------------------
-      IF(NJETS.GT.0 .AND.KT_DURHAM.GT.0D0) THEN
-C RESET JET MOMENTA
-      njets=0
-      DO I=1,NEXTERNAL
-        DO J=0,3
-          PJET(I,J) = 0E0
-        ENDDO
-      ENDDO
-
-      do i=nincoming+1,nexternal
-         if(is_a_j(i).and.do_cuts(i)) then
-           njets=njets+1
-           DO J=0,3
-             PJET(NJETS,J) = P(J,I)
-           ENDDO
-         endif
-      enddo
-
-C DURHAM KT SEPARATION CUT
-
-
-        PTMIN = EBEAM(1) + EBEAM(2)
-
-        DO I=1,NJETS
-
-C         PT WITH RESPECT TO Z AXIS FOR HADRONIC COLLISIONS
-          IF ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) THEN
-            PT1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2)
-            PTMIN = MIN( PTMIN, PT1 )
-          ENDIF
-
-          DO J=I+1,NJETS
-C           GET ANGLE BETWEEN JETS
-            PABS1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2 + PJET(I,3)**2)
-            PABS2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2 + PJET(J,3)**2)
-C           CHECK IF 3-MOMENTA DO NOT VANISH
-            IF(PABS1*PABS2 .NE. 0D0) THEN
-              COSTH = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) + PJET(I,3)*PJET(J,3) )/(PABS1*PABS2)
-            ELSE
-C           IF 3-MOMENTA VANISH, MAKE JET COSTH = 1D0 SO THAT JET MEASURE VANISHES
-              COSTH = 1D0
-            ENDIF
-C           GET PT AND ETA OF JETS
-            PT2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2)
-            ETA1 = 0.5D0*LOG( (PJET(I,0) + PJET(I,3)) / (PJET(I,0) - PJET(I,3)) )
-            ETA2 = 0.5D0*LOG( (PJET(J,0) + PJET(J,3)) / (PJET(J,0) - PJET(J,3)) )
-C           GET COSH OF DELTA ETA, COS OF DELTA PHI
-            COSH_DETA = DCOSH( ETA1 - ETA2 )
-            COS_DPHI = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) ) / (PT1*PT2)
-            DPHI = DACOS( COS_DPHI )
-            IF ( (LPP(1).EQ.0) .AND. (LPP(2).EQ.0)) THEN
-C             KT FOR E+E- COLLISION
-              PTNOW = DSQRT( 2D0*MIN(PJET(I,0)**2,PJET(J,0)**2)*( 1D0-COSTH ) )
-             ELSE
-C             HADRONIC KT, FASTJET DEFINITION
-              PTNOW = DSQRT( MIN(PT1**2,PT2**2)*( (ETA1 - ETA2 )**2 + DPHI**2 )/(D_PARAMETER**2) )
-            ENDIF
-
-            PTMIN = MIN( PTMIN, PTNOW )
-
-          ENDDO ! LOOP OVER NJET
-
-        ENDDO ! LOOP OVER NJET
-
-C CHECK COMPATIBILITY WITH CUT
-        IF( (PTMIN .LT. KT_DURHAM)) THEN
-          PASSCUTS = .FALSE.
-          RETURN
-        ENDIF
-      ENDIF ! IF NJETS.GT. 0 .AND. DO_KT_DURHAM
-
-C----------------------------------------------------------------------------
-C----------------------------------------------------------------------------
-
-
-
-c- check existance of jets if jet cuts are on
-      if(njets.lt.1.and.(htjmin.gt.0.or.ptj1min.gt.0).or.
-     $     njets.lt.2.and.ptj2min.gt.0.or.
-     $     njets.lt.3.and.ptj3min.gt.0.or.
-     $     njets.lt.4.and.ptj4min.gt.0)then
-         if(debug) write (*,*) i, ' too few jets -> fails'
-         passcuts=.false.
-         return
-      endif
-
-c - sort jet pts
-      do i=1,njets-1
-         do j=i+1,njets
-            if(ptjet(j).gt.ptjet(i)) then
-               temp=ptjet(i)
-               ptjet(i)=ptjet(j)
-               ptjet(j)=temp
-            endif
-         enddo
-      enddo
-      if(debug) write (*,*) 'ordered ',njets,'   ',ptjet
-c
-c     Use "and" or "or" prescriptions 
-c     
-      inclht=0
-
-      if(njets.gt.0) then
-
-       notgood=.not.jetor
-       if(debug) write (*,*) 'jetor :',jetor  
-       if(debug) write (*,*) '0',notgood   
-      
-      do i=1,min(njets,4)
-            if(debug) write (*,*) i,ptjet(i), '   ',ptjmin4(i),
-     $        ':',ptjmax4(i)
-         if(jetor) then     
-c---  if one of the jets does not pass, the event is rejected
-            notgood=notgood.or.(ptjmax4(i).ge.0d0 .and.
-     $           ptjet(i).gt.ptjmax4(i)) .or.
-     $           (ptjet(i).lt.ptjmin4(i))
-            if(debug) write (*,*) i,' notgood total:', notgood   
-         else
-c---  all cuts must fail to reject the event
-            notgood=notgood.and.(ptjmax4(i).ge.0d0 .and.
-     $              ptjet(i).gt.ptjmax4(i) .or.
-     $              (ptjet(i).lt.ptjmin4(i)))
-            if(debug) write (*,*) i,' notgood total:', notgood   
-         endif
-      enddo
-
-
-      if (notgood) then
-         if(debug) write (*,*) i, ' multiple pt -> fails'
-         passcuts=.false.
-         return
-      endif
-
-c---------------------------
-c      Ht cuts
-C---------------------------
-      htj=ptjet(1)
-
-      do i=2,njets
-         htj=htj+ptjet(i)
-         if(debug) write (*,*) i, 'htj ',htj
-         if(debug.and.i.le.4) write (*,*) 'htmin ',i,' ', htjmin4(i),':',htjmax4(i)
-         if(i.le.4)then
-            if(htj.lt.htjmin4(i) .or.
-     $        htjmax4(i).ge.0d0.and.htj.gt.htjmax4(i)) then
-            if(debug) write (*,*) i, ' ht -> fails'
-            passcuts=.false.
-            return
-            endif
-         endif
-      enddo
-
-      if(htj.lt.htjmin.or.htjmax.ge.0d0.and.htj.gt.htjmax)then
-         if(debug) write (*,*) i, ' htj -> fails'
-         passcuts=.false.
-         return
-      endif
-
-      inclht=htj
-
-      endif !if there are jets 
-
-      if(nheavyjets.gt.0) then
-         do i=1,nheavyjets
-            inclht=inclht+ptheavyjet(i)
-         enddo
-      endif !if there are heavyjets
-
-      if(inclht.lt.inclHtmin.or.
-     $     inclHtmax.ge.0d0.and.inclht.gt.inclHtmax)then
-         if(debug) write (*,*) ' inclhtmin=',inclHtmin,' -> fails'
-         passcuts=.false.
-         return
-      endif
-
-C     
-C     maximal and minimal pt of the leptons sorted by pt
-c     
-      nleptons=0
-
-      if(ptl1min.gt.0.or.ptl2min.gt.0.or.ptl3min.gt.0.or.ptl4min.gt.0.or.
-     $     ptl1max.ge.0d0.or.ptl2max.ge.0d0.or.
-     $     ptl3max.ge.0d0.or.ptl4max.ge.0d0) then
-
-c     - fill ptlepton with the pt's of the leptons.
-         do i=nincoming+1,nexternal
-            if(is_a_l(i)) then
-               nleptons=nleptons+1
-               ptlepton(nleptons)=pt(p(0,i))
-            endif
-         enddo
-         if(debug) write (*,*) 'not yet ordered ',njets,'   ',ptjet
-
-c     - check existance of leptons if lepton cuts are on
-         if(nleptons.lt.1.and.ptl1min.gt.0.or.
-     $        nleptons.lt.2.and.ptl2min.gt.0.or.
-     $        nleptons.lt.3.and.ptl3min.gt.0.or.
-     $        nleptons.lt.4.and.ptl4min.gt.0)then
-            if(debug) write (*,*) i, ' too few leptons -> fails'
-            passcuts=.false.
-            return
-         endif
-
-c     - sort lepton pts
-         do i=1,nleptons-1
-            do j=i+1,nleptons
-               if(ptlepton(j).gt.ptlepton(i)) then
-                  temp=ptlepton(i)
-                  ptlepton(i)=ptlepton(j)
-                  ptlepton(j)=temp
-               endif
-            enddo
-         enddo
-         if(debug) write (*,*) 'ordered ',nleptons,'   ',ptlepton
-
-         if(nleptons.gt.0) then
-
-            notgood = .false.
-            do i=1,min(nleptons,4)
-               if(debug) write (*,*) i,ptlepton(i), '   ',ptlmin4(i),':',ptlmax4(i)
-c---  if one of the leptons does not pass, the event is rejected
-               notgood=notgood.or.
-     $              (ptlmax4(i).ge.0d0.and.ptlepton(i).gt.ptlmax4(i)).or.
-     $              (ptlepton(i).lt.ptlmin4(i))
-               if(debug) write (*,*) i,' notgood total:', notgood   
-            enddo
-
-
-            if (notgood) then
-               if(debug) write (*,*) i, ' multiple pt -> fails'
-               passcuts=.false.
-               return
-            endif
-         endif
-      endif
-C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>     
-C     SPECIAL CUTS
-C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-C     REQUIRE AT LEAST ONE JET WITH PT>XPTJ
-         
-         IF(xptj.gt.0d0) THEN
-            xvar=0
-            do i=nincoming+1,nexternal
-               if(is_a_j(i)) xvar=max(xvar,pt(p(0,i)))
-            enddo
-            if (xvar .lt. xptj) then
-               passcuts=.false.
-               return
-            endif
-         ENDIF
-
-C     REQUIRE AT LEAST ONE PHOTON WITH PT>XPTA
-         
-         IF(xpta.gt.0d0) THEN
-            xvar=0
-            do i=nincoming+1,nexternal
-               if(is_a_a(i)) xvar=max(xvar,pt(p(0,i)))
-            enddo
-            if (xvar .lt. xpta) then
-               passcuts=.false.
-               return
-            endif
-         ENDIF
-
-C     REQUIRE AT LEAST ONE B  WITH PT>XPTB
-         
-         IF(xptb.gt.0d0) THEN
-            xvar=0
-            do i=nincoming+1,nexternal
-               if(is_a_b(i)) xvar=max(xvar,pt(p(0,i)))
-            enddo
-            if (xvar .lt. xptb) then
-               passcuts=.false.
-               return
-            endif
-         ENDIF
-
-C     REQUIRE AT LEAST ONE LEPTON  WITH PT>XPTL
-         
-         IF(xptl.gt.0d0) THEN
-            xvar=0
-            do i=nincoming+1,nexternal
-               if(is_a_l(i)) xvar=max(xvar,pt(p(0,i)))
-            enddo
-            if (xvar .lt. xptl) then
-               passcuts=.false.
-               if(debug) write (*,*) ' xptl -> fails'
-               return
-            endif
-         ENDIF
-C
-C     WBF CUTS: TWO TYPES
-C    
-C     FIRST TYPE:  implemented by FM
-C
-C     1. FIND THE 2 HARDEST JETS
-C     2. REQUIRE |RAP(J)|>XETAMIN
-C     3. REQUIRE RAP(J1)*ETA(J2)<0
-C
-C     SECOND TYPE : added by Simon de Visscher 1-08-2007
-C
-C     1. FIND THE 2 HARDEST JETS
-C     2. REQUIRE |RAP(J1)-RAP(J2)|>DELTAETA
-C     3. REQUIRE RAP(J1)*RAP(J2)<0
-C
-C
-         hardj1=0
-         hardj2=0
-         ptmax1=0d0
-         ptmax2=0d0
-
-C-- START IF AT LEAST ONE OF THE CUTS IS ACTIVATED
-         
-         IF(XETAMIN.GT.0D0.OR.DELTAETA.GT.0D0) THEN
-            
-C-- FIND THE HARDEST JETS
-
-            do i=nincoming+1,nexternal
-               if(is_a_j(i)) then
-c                  write (*,*) i,pt(p(0,i))
-                  if(pt(p(0,i)).gt.ptmax1) then
-                     hardj2=hardj1
-                     ptmax2=ptmax1
-                     hardj1=i
-                     ptmax1=pt(p(0,i))
-                  elseif(pt(p(0,i)).gt.ptmax2) then
-                     hardj2=i
-                     ptmax2=pt(p(0,i))
-                  endif
-c                  write (*,*) hardj1,hardj2,ptmax1,ptmax2
-               endif
-            enddo
-            
-            if (hardj2.eq.0) goto 21 ! bypass vbf cut since not enough jets
-
-C-- NOW APPLY THE CUT I            
-
-            if (abs(rap(p(0,hardj1))) .lt. xetamin
-     &       .or.abs(rap(p(0,hardj2))) .lt. xetamin
-     &       .or.rap(p(0,hardj1))*rap(p(0,hardj2)) .gt.0d0) then
-             passcuts=.false.
-             return
-            endif
-
-            
-C-- NOW APPLY THE CUT II
-            
-            if (abs(rap(p(0,hardj1))-rap(p(0,hardj2))) .lt. deltaeta) then
-             passcuts=.false.
-             return
-            endif
-         
-c            write (*,*) hardj1,hardj2,rap(p(0,hardj1)),rap(p(0,hardj2))
-         
-         ENDIF
-
-c Begin photon isolation
-c NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE
-c     Use is made of parton cm frame momenta. If this must be
-c     changed, pQCD used below must be redefined
-c NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE
-c If we do not require a mimimum jet energy, there's no need to apply
-c jet clustering and all that.
-      if (ptgmin.ne.0d0) then
-
-c Put all (light) QCD partons in momentum array for jet clustering.
-c From the run_card.dat, maxjetflavor defines if b quark should be
-c considered here (via the logical variable 'is_a_jet').  nQCD becomes
-c the number of (light) QCD partons at the real-emission level (i.e. one
-c more than the Born).
-
-      nQCD=0
-      do j=nincoming+1,nexternal
-         if (is_a_j(j)) then
-            nQCD=nQCD+1
-            do i=0,3
-               pQCD(i,nQCD)=p(i,j) ! Use C.o.M. frame momenta
-            enddo
-         endif
-      enddo
-
-        nph=0
-        do j=nincoming+1,nexternal
-          if (is_a_a(j)) then
-            nph=nph+1
-            do i=0,3
-              pgamma(i,nph)=p(i,j) ! Use C.o.M. frame momenta
-            enddo
-          endif
-        enddo
-        if(nph.eq.0) goto 444
-
-        if(isoEM)then
-          nem=nph
-          do k=1,nem
-            do i=0,3
-              pem(i,k)=pgamma(i,k)
-            enddo
-          enddo
-          do j=nincoming+1,nexternal
-            if (is_a_l(j)) then
-              nem=nem+1
-              do i=0,3
-                pem(i,nem)=p(i,j) ! Use C.o.M. frame momenta
-              enddo
-            endif
-          enddo
-        endif
-
-        alliso=.true.
-
-        j=0
-        dowhile(j.lt.nph.and.alliso)
-c Loop over all photons
-          j=j+1
-
-          ptg=pt(pgamma(0,j))
-          if(ptg.lt.ptgmin)then
-            passcuts=.false.
-            return
-          endif
-
-c Isolate from hadronic energy
-          do i=1,nQCD
-            drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pQCD(0,i)))
-          enddo
-          call sortzv(drlist,isorted,nQCD,ismode,isway,izero)
-          Etsum(0)=0.d0
-          nin=0
-          do i=1,nQCD
-            if(dble(drlist(isorted(i))).le.R0gamma)then
-              nin=nin+1
-              Etsum(nin)=Etsum(nin-1)+pt(pQCD(0,isorted(i)))
-            endif
-          enddo
-          do i=1,nin
-            alliso=alliso .and.
-     #        Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
-     #                                  R0gamma,xn,epsgamma,ptg)
-          enddo
-
-c Isolate from EM energy
-          if(isoEM.and.nem.gt.1)then
-            do i=1,nem
-              drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i)))
-            enddo
-            call sortzv(drlist,isorted,nem,ismode,isway,izero)
-c First of list must be the photon: check this, and drop it
-            if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then
-              write(*,*)'Error #1 in photon isolation'
-              write(*,*)j,isorted(1),drlist(isorted(1))
-              stop
-            endif
-            Etsum(0)=0.d0
-            nin=0
-            do i=2,nem
-              if(dble(drlist(isorted(i))).le.R0gamma)then
-                nin=nin+1
-                Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i)))
-              endif
-            enddo
-            do i=1,nin
-              alliso=alliso .and.
-     #          Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
-     #                                    R0gamma,xn,epsgamma,ptg)
-            enddo
-
-          endif
-
-c End of loop over photons
-        enddo
-
-        if(.not.alliso)then
-          passcuts=.false.
-          return
-        endif
-      endif
-
- 444    continue
-c End photon isolation
-
-
-C...Set couplings if event passed cuts
-
- 21   if(.not.fixed_ren_scale) then
-         call set_ren_scale(P,scale)
-         if(scale.gt.0) G = SQRT(4d0*PI*ALPHAS(scale))
-      endif
-
-      if(.not.fixed_fac_scale) then
-         call set_fac_scale(P,q2fact)
-      endif
-
-c
-c     Here we cluster event and reset factorization and renormalization
-c     scales on an event-by-event basis, as well as check xqcut for jets
-c
-c     Note the following condition is the first line of setclscales
-c      if(xqcut.gt.0d0.or.ickkw.gt.0.or.scale.eq.0.or.q2fact(1).eq.0)then
-c     Do not duplicate it since some variable are set for syscalc in the fct
-        if(.not.setclscales(p))then
-           cutsdone=.false.
-           cutspassed=.false.
-           passcuts = .false.
-           if(debug) write (*,*) 'setclscales -> fails'
-           return
-       endif
-c      endif
-
-c     Set couplings in model files
-      if(.not.fixed_ren_scale.or..not.fixed_couplings) then
-         if (.not.fixed_couplings)then
-            do i=0,3
-               do j=1,nexternal
-                  pp(i,j)=p(i,j)
-               enddo
-            enddo
-         endif
-         call update_as_param()
-      endif
-
-      IF (FIRSTTIME2) THEN
-        FIRSTTIME2=.FALSE.
-        write(6,*) 'alpha_s for scale ',scale,' is ', G**2/(16d0*atan(1d0))
-      ENDIF
-
-      if(debug) write (*,*) '============================='
-      if(debug) write (*,*) ' EVENT PASSED THE CUTS       '
-      if(debug) write (*,*) '============================='
-
-      CUTSPASSED=.TRUE.
-
-      RETURN
-      END
-
-
-C
-C     FUNCTION FOR ISOLATION
-C
-
-      function iso_getdrv40(p1,p2)
-      implicit none
-      real*8 iso_getdrv40,p1(0:3),p2(0:3)
-      real*8 iso_getdr
-c
-      iso_getdrv40=iso_getdr(p1(0),p1(1),p1(2),p1(3),
-     #                       p2(0),p2(1),p2(2),p2(3))
-      return
-      end
-
-
-      function iso_getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
-      implicit none
-      real*8 iso_getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
-     # iso_getpseudorap,iso_getdelphi
-c
-      deta=iso_getpseudorap(en1,ptx1,pty1,pl1)-
-     #     iso_getpseudorap(en2,ptx2,pty2,pl2)
-      dphi=iso_getdelphi(ptx1,pty1,ptx2,pty2)
-      iso_getdr=sqrt(dphi**2+deta**2)
-      return
-      end
-
-
-      function iso_getpseudorap(en,ptx,pty,pl)
-      implicit none
-      real*8 iso_getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
-      parameter (tiny=1.d-5)
-c
-      pt=sqrt(ptx**2+pty**2)
-      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
-        eta=sign(1.d0,pl)*1.d8
-      else
-        th=atan2(pt,pl)
-        eta=-log(tan(th/2.d0))
-      endif
-      iso_getpseudorap=eta
-      return
-      end
-
-
-      function iso_getdelphi(ptx1,pty1,ptx2,pty2)
-      implicit none
-      real*8 iso_getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
-      parameter (tiny=1.d-5)
-c
-      pt1=sqrt(ptx1**2+pty1**2)
-      pt2=sqrt(ptx2**2+pty2**2)
-      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
-        tmp=ptx1*ptx2+pty1*pty2
-        tmp=tmp/(pt1*pt2)
-        if(abs(tmp).gt.1.d0+tiny)then
-          write(*,*)'Cosine larger than 1'
-          stop
-        elseif(abs(tmp).ge.1.d0)then
-          tmp=sign(1.d0,tmp)
-        endif
-        tmp=acos(tmp)
-      else
-        tmp=1.d8
-      endif
-      iso_getdelphi=tmp
-      return
-      end
-
-      function chi_gamma_iso(dr,R0,xn,epsgamma,pTgamma)
-c Eq.(3.4) of Phys.Lett. B429 (1998) 369-374 [hep-ph/9801442]
-      implicit none
-      real*8 chi_gamma_iso,dr,R0,xn,epsgamma,pTgamma
-      real*8 tmp,axn
-c
-      axn=abs(xn)
-      tmp=epsgamma*pTgamma
-      if(axn.ne.0.d0)then
-        tmp=tmp*( (1-cos(dr))/(1-cos(R0)) )**axn
-      endif
-      chi_gamma_iso=tmp
-      return
-      end
-
-
-*
-* $Id: sortzv.F,v 1.1.1.1 1996/02/15 17:49:50 mclareni Exp $
-*
-* $Log: sortzv.F,v $
-* Revision 1.1.1.1  1996/02/15 17:49:50  mclareni
-* Kernlib
-*
-*
-c$$$#include "kerngen/pilot.h"
-      SUBROUTINE SORTZV (A,INDEX,N1,MODE,NWAY,NSORT)
-C
-C CERN PROGLIB# M101    SORTZV          .VERSION KERNFOR  3.15  820113
-C ORIG. 02/10/75
-C
-      DIMENSION A(N1),INDEX(N1)
-C
-C
-      N = N1
-      IF (N.LE.0)            RETURN
-      IF (NSORT.NE.0) GO TO 2
-      DO 1 I=1,N
-    1 INDEX(I)=I
-C
-    2 IF (N.EQ.1)            RETURN
-      IF (MODE)    10,20,30
-   10 CALL SORTTI (A,INDEX,N)
-      GO TO 40
-C
-   20 CALL SORTTC(A,INDEX,N)
-      GO TO 40
-C
-   30 CALL SORTTF (A,INDEX,N)
-C
-   40 IF (NWAY.EQ.0) GO TO 50
-      N2 = N/2
-      DO 41 I=1,N2
-      ISWAP = INDEX(I)
-      K = N+1-I
-      INDEX(I) = INDEX(K)
-   41 INDEX(K) = ISWAP
-   50 RETURN
-      END
-*     ========================================
-      SUBROUTINE SORTTF (A,INDEX,N1)
-C
-      DIMENSION A(N1),INDEX(N1)
-C
-      N = N1
-      DO 3 I1=2,N
-      I3 = I1
-      I33 = INDEX(I3)
-      AI = A(I33)
-    1 I2 = I3/2
-      IF (I2) 3,3,2
-    2 I22 = INDEX(I2)
-      IF (AI.LE.A (I22)) GO TO 3
-      INDEX (I3) = I22
-      I3 = I2
-      GO TO 1
-    3 INDEX (I3) = I33
-    4 I3 = INDEX (N)
-      INDEX (N) = INDEX (1)
-      AI = A(I3)
-      N = N-1
-      IF (N-1) 12,12,5
-    5 I1 = 1
-    6 I2 = I1 + I1
-      IF (I2.LE.N) I22= INDEX(I2)
-      IF (I2-N) 7,9,11
-    7 I222 = INDEX (I2+1)
-      IF (A(I22)-A(I222)) 8,9,9
-    8 I2 = I2+1
-      I22 = I222
-    9 IF (AI-A(I22)) 10,11,11
-   10 INDEX(I1) = I22
-      I1 = I2
-      GO TO 6
-   11 INDEX (I1) = I3
-      GO TO 4
-   12 INDEX (1) = I3
-      RETURN
-      END
-*     ========================================
-      SUBROUTINE SORTTI (A,INDEX,N1)
-C
-      INTEGER A,AI
-      DIMENSION A(N1),INDEX(N1)
-C
-      N = N1
-      DO 3 I1=2,N
-      I3 = I1
-      I33 = INDEX(I3)
-      AI = A(I33)
-    1 I2 = I3/2
-      IF (I2) 3,3,2
-    2 I22 = INDEX(I2)
-      IF (AI.LE.A (I22)) GO TO 3
-      INDEX (I3) = I22
-      I3 = I2
-      GO TO 1
-    3 INDEX (I3) = I33
-    4 I3 = INDEX (N)
-      INDEX (N) = INDEX (1)
-      AI = A(I3)
-      N = N-1
-      IF (N-1) 12,12,5
-    5 I1 = 1
-    6 I2 = I1 + I1
-      IF (I2.LE.N) I22= INDEX(I2)
-      IF (I2-N) 7,9,11
-    7 I222 = INDEX (I2+1)
-      IF (A(I22)-A(I222)) 8,9,9
-    8 I2 = I2+1
-      I22 = I222
-    9 IF (AI-A(I22)) 10,11,11
-   10 INDEX(I1) = I22
-      I1 = I2
-      GO TO 6
-   11 INDEX (I1) = I3
-      GO TO 4
-   12 INDEX (1) = I3
-      RETURN
-      END
-*     ========================================
-      SUBROUTINE SORTTC (A,INDEX,N1)
-C
-      INTEGER A,AI
-      DIMENSION A(N1),INDEX(N1)
-C
-      N = N1
-      DO 3 I1=2,N
-      I3 = I1
-      I33 = INDEX(I3)
-      AI = A(I33)
-    1 I2 = I3/2
-      IF (I2) 3,3,2
-    2 I22 = INDEX(I2)
-      IF(ICMPCH(AI,A(I22)))3,3,21
-   21 INDEX (I3) = I22
-      I3 = I2
-      GO TO 1
-    3 INDEX (I3) = I33
-    4 I3 = INDEX (N)
-      INDEX (N) = INDEX (1)
-      AI = A(I3)
-      N = N-1
-      IF (N-1) 12,12,5
-    5 I1 = 1
-    6 I2 = I1 + I1
-      IF (I2.LE.N) I22= INDEX(I2)
-      IF (I2-N) 7,9,11
-    7 I222 = INDEX (I2+1)
-      IF (ICMPCH(A(I22),A(I222))) 8,9,9
-    8 I2 = I2+1
-      I22 = I222
-    9 IF (ICMPCH(AI,A(I22))) 10,11,11
-   10 INDEX(I1) = I22
-      I1 = I2
-      GO TO 6
-   11 INDEX (I1) = I3
-      GO TO 4
-   12 INDEX (1) = I3
-      RETURN
-      END
-*     ========================================
-      FUNCTION ICMPCH(IC1,IC2)
-C     FUNCTION TO COMPARE TWO 4 CHARACTER EBCDIC STRINGS - IC1,IC2
-C     ICMPCH=-1 IF HEX VALUE OF IC1 IS LESS THAN IC2
-C     ICMPCH=0  IF HEX VALUES OF IC1 AND IC2 ARE THE SAME
-C     ICMPCH=+1 IF HEX VALUES OF IC1 IS GREATER THAN IC2
-      I1=IC1
-      I2=IC2
-      IF(I1.GE.0.AND.I2.GE.0)GOTO 40
-      IF(I1.GE.0)GOTO 60
-      IF(I2.GE.0)GOTO 80
-      I1=-I1
-      I2=-I2
-      IF(I1-I2)80,70,60
- 40   IF(I1-I2)60,70,80
- 60   ICMPCH=-1
-      RETURN
- 70   ICMPCH=0
-      RETURN
- 80   ICMPCH=1
-      RETURN
-      END
-
+      logical function pass_point(p)
+c************************************************************************
+c     This function is called from sample to see if it needs to 
+c     bother calculating the weight from all the different conficurations
+c     You can either just return true, or have it call passcuts
+c************************************************************************
+      implicit none
+c
+c     Arguments
+c
+      double precision p
+c
+c     External
+c
+      logical passcuts
+      external passcuts
+c-----
+c  Begin Code
+c-----
+      pass_point = .true.
+c      pass_point = passcuts(p)
+      end
+C 
+      LOGICAL FUNCTION PASSCUTS(P)
+C**************************************************************************
+C     INPUT:
+C            P(0:3,1)           MOMENTUM OF INCOMING PARTON
+C            P(0:3,2)           MOMENTUM OF INCOMING PARTON
+C            P(0:3,3)           MOMENTUM OF ...
+C            ALL MOMENTA ARE IN THE REST FRAME!!
+C            COMMON/JETCUTS/   CUTS ON JETS
+C     OUTPUT:
+C            TRUE IF EVENTS PASSES ALL CUTS LISTED
+C**************************************************************************
+      IMPLICIT NONE
+c
+c     Constants
+c
+      include 'genps.inc'
+      include 'nexternal.inc'
+C
+C     ARGUMENTS
+C
+      REAL*8 P(0:3,nexternal)
+
+C
+C     LOCAL
+C
+      LOGICAL FIRSTTIME,FIRSTTIME2,pass_bw,notgood,good,foundheavy
+      LOGICAL DEBUG
+      integer i,j,njets,nheavyjets,nleptons,hardj1,hardj2
+      REAL*8 XVAR,ptmax1,ptmax2,htj,tmp,inclht
+      real*8 ptemp(0:3), ptemp2(0:3)
+      character*20 formstr
+C
+C     PARAMETERS
+C
+      real*8 PI
+      parameter( PI = 3.14159265358979323846d0 )
+C
+C     EXTERNAL
+C
+      REAL*8 R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,PtDot
+      logical cut_bw,setclscales
+      external R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,cut_bw,setclscales,PtDot
+C
+C     GLOBAL
+C
+      include 'run.inc'
+      include 'cuts.inc'
+
+      double precision ptjet(nexternal)
+      double precision ptheavyjet(nexternal)
+      double precision ptlepton(nexternal)
+      double precision temp
+
+C     VARIABLES TO SPECIFY JETS
+      DOUBLE PRECISION PJET(NEXTERNAL,0:3)
+      DOUBLE PRECISION PTMIN
+      DOUBLE PRECISION PT1,PT2
+
+C     INTEGERS FOR COUNTING.
+      INTEGER K,L,J1,J2
+
+C VARIABLES FOR KT CUT
+      DOUBLE PRECISION PTNOW,COSTH,PABS1,PABS2
+      DOUBLE PRECISION ETA1,ETA2,COSH_DETA,COS_DPHI,KT1SQ,KT2SQ, DPHI
+
+      double precision etmin(nincoming+1:nexternal),etamax(nincoming+1:nexternal)
+      double precision emin(nincoming+1:nexternal)
+      double precision                    r2min(nincoming+1:nexternal,nincoming+1:nexternal)
+      double precision s_min(nexternal,nexternal)
+      double precision etmax(nincoming+1:nexternal),etamin(nincoming+1:nexternal)
+      double precision emax(nincoming+1:nexternal)
+      double precision r2max(nincoming+1:nexternal,nincoming+1:nexternal)
+      double precision s_max(nexternal,nexternal)
+      double precision ptll_min(nexternal,nexternal),ptll_max(nexternal,nexternal)
+      double precision inclHtmin,inclHtmax
+      common/to_cuts/  etmin, emin, etamax, r2min, s_min,
+     $     etmax, emax, etamin, r2max, s_max, ptll_min, ptll_max, inclHtmin,inclHtmax
+
+      double precision ptjmin4(4),ptjmax4(4),htjmin4(2:4),htjmax4(2:4)
+      logical jetor
+      common/to_jet_cuts/ ptjmin4,ptjmax4,htjmin4,htjmax4,jetor
+
+      double precision ptlmin4(4),ptlmax4(4)
+      common/to_lepton_cuts/ ptlmin4,ptlmax4
+
+c
+c     Special cuts
+c
+
+      integer        lbw(0:nexternal)  !Use of B.W.
+      common /to_BW/ lbw
+C
+C     SPECIAL CUTS
+C
+      LOGICAL  IS_A_J(NEXTERNAL),IS_A_L(NEXTERNAL)
+      LOGICAL  IS_A_B(NEXTERNAL),IS_A_A(NEXTERNAL),IS_A_ONIUM(NEXTERNAL)
+      LOGICAL  IS_A_NU(NEXTERNAL),IS_HEAVY(NEXTERNAL)
+      logical  do_cuts(nexternal)
+      COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,
+     . IS_A_ONIUM, do_cuts
+
+C
+C     MERGING SCALE CUT
+C
+C     Retrieve which external particles undergo the ktdurham and ptlund cuts.
+      LOGICAL  is_pdg_for_merging_cut(NEXTERNAL)
+      logical from_decay(-(nexternal+3):nexternal)
+      COMMON /TO_MERGE_CUTS/is_pdg_for_merging_cut, from_decay
+
+C
+C     ADDITIONAL VARIABLES FOR PTLUND CUT
+C
+      INTEGER NMASSLESS
+      DOUBLE PRECISION PINC(NINCOMING,0:3)
+      DOUBLE PRECISION PRADTEMP(0:3), PRECTEMP(0:3), PEMTTEMP(0:3)
+      DOUBLE PRECISION PTMINSAVE, RHOPYTHIA
+      EXTERNAL RHOPYTHIA
+C
+C     FLAVOUR INFORMATION NECESSARY TO RECONSTRUCT PTLUND
+C
+      INTEGER JETFLAVOUR(NEXTERNAL), INCFLAVOUR(NINCOMING)
+      include 'maxamps.inc'
+      integer idup(nexternal,maxproc,maxsproc)
+      integer mothup(2,nexternal)
+      integer icolup(2,nexternal,maxflow,maxsproc)
+      include 'leshouche.inc'
+
+C
+C     Keep track of whether cuts already calculated for this event
+C
+      LOGICAL CUTSDONE,CUTSPASSED
+      COMMON/TO_CUTSDONE/CUTSDONE,CUTSPASSED
+      DATA CUTSDONE,CUTSPASSED/.FALSE.,.FALSE./
+
+C $B$ MW_NEW_DEF $E$ !this is a tag for MadWeight
+
+      double precision xqcutij(nexternal,nexternal),xqcuti(nexternal)
+      common/to_xqcuts/xqcutij,xqcuti
+
+c jet cluster algorithm
+      integer nQCD !,NJET,JET(nexternal)
+c      double precision plab(0:3, nexternal)
+      double precision pQCD(0:3,nexternal)!,PJET(0:3,nexternal)
+c      double precision rfj,sycut,palg,fastjetdmerge
+c      integer njet_eta
+c     Photon isolation
+      integer nph,nem,nin
+      double precision ptg,chi_gamma_iso,iso_getdrv40
+      double precision Etsum(0:nexternal)
+      real drlist(nexternal)
+      double precision pgamma(0:3,nexternal),pem(0:3,nexternal)
+      logical alliso
+C     Sort array of results: ismode>0 for real, isway=0 for ascending order
+      integer ismode,isway,izero,isorted(nexternal)
+      parameter (ismode=1)
+      parameter (isway=0)
+      parameter (izero=0)
+
+      include 'coupl.inc'
+
+C
+C
+c
+      DATA FIRSTTIME,FIRSTTIME2/.TRUE.,.TRUE./
+
+c put momenta in common block for couplings.f
+      double precision pp(0:3,max_particles)
+      common /momenta_pp/pp
+
+      DATA DEBUG/.FALSE./
+
+C-----
+C  BEGIN CODE
+C-----
+
+
+
+      PASSCUTS=.TRUE.             !EVENT IS OK UNLESS OTHERWISE CHANGED
+      IF (FIRSTTIME) THEN
+         FIRSTTIME=.FALSE.
+c      Preparation for reweighting by setting up clustering by diagrams
+         call initcluster()
+c
+c
+         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'i8)'
+         write(*,formstr) 'Particle',(i,i=nincoming+1,nexternal)
+         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'f8.1)'
+         write(*,formstr) 'Et >',(etmin(i),i=nincoming+1,nexternal)
+         write(*,formstr) 'E >',(emin(i),i=nincoming+1,nexternal)
+         write(*,formstr) 'Eta <',(etamax(i),i=nincoming+1,nexternal)
+         write(*,formstr) 'xqcut: ',(xqcuti(i),i=nincoming+1,nexternal)
+         write(formstr,'(a,i2.2,a)')'(a,i2,a,',nexternal,'f8.1)'
+         do j=nincoming+1,nexternal-1
+            write(*,formstr) 'd R #',j,'  >',(-0.0,i=nincoming+1,j),
+     &           (r2min(i,j),i=j+1,nexternal)
+            do i=j+1,nexternal
+               r2min(i,j)=r2min(i,j)*dabs(r2min(i,j))    !Since r2 returns distance squared
+               r2max(i,j)=r2max(i,j)*dabs(r2max(i,j))
+            enddo
+         enddo
+         do j=nincoming+1,nexternal-1
+            write(*,formstr) 's min #',j,'>',
+     &           (s_min(i,j),i=nincoming+1,nexternal)
+         enddo
+         do j=nincoming+1,nexternal-1
+            write(*,formstr) 'xqcutij #',j,'>',
+     &           (xqcutij(i,j),i=nincoming+1,nexternal)
+         enddo
+
+cc
+cc     Set the strong coupling
+cc
+c         call set_ren_scale(P,scale)
+c
+cc     Check that the user funtions for setting the scales
+cc     have been edited if the choice of an event-by-event
+cc     scale choice has been made 
+c
+c         if(.not.fixed_ren_scale) then
+c            if(scale.eq.0d0) then
+c               write(6,*) 
+c               write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
+c               write(6,*) ' Dynamical renormalization scale choice '
+c               write(6,*) ' selected but user subroutine' 
+c               write(6,*) ' set_ren_scale not edited in file:setpara.f'
+c               write(6,*) ' Switching to a fixed_ren_scale choice'
+c               write(6,*) ' with scale=zmass'
+c               scale=91.2d0
+c               write(6,*) 'scale=',scale
+c               fixed_ren_scale=.true.
+c               call set_ren_scale(P,scale)
+c            endif
+c         endif
+         
+c         if(.not.fixed_fac_scale) then
+c            call set_fac_scale(P,q2fact)
+c            if(q2fact(1).eq.0d0.or.q2fact(2).eq.0d0) then
+c               write(6,*) 
+c               write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
+c               write(6,*) ' Dynamical renormalization scale choice '
+c               write(6,*) ' selected but user subroutine' 
+c               write(6,*) ' set_fac_scale not edited in file:setpara.f'
+c               write(6,*) ' Switching to a fixed_fac_scale choice'
+c               write(6,*) ' with q2fact(i)=zmass**2'
+c               fixed_fac_scale=.true.
+c               q2fact(1)=91.2d0**2
+c               q2fact(2)=91.2d0**2
+c               write(6,*) 'scales=',q2fact(1),q2fact(2)
+c            endif
+c         endif
+
+         if(fixed_ren_scale) then
+            G = SQRT(4d0*PI*ALPHAS(scale))
+            call update_as_param()
+         endif
+
+c     Put momenta in the common block to zero to start
+         do i=0,3
+            do j=1,max_particles
+               pp(i,j) = 0d0
+            enddo
+         enddo
+         
+      ENDIF ! IF FIRSTTIME
+
+      IF (CUTSDONE) THEN
+         PASSCUTS=CUTSPASSED
+         RETURN
+      ENDIF
+      CUTSDONE=.TRUE.
+c      CUTSPASSED=.FALSE.
+
+c
+c     Make sure have reasonable 4-momenta
+c
+      if (p(0,1) .le. 0d0) then
+         passcuts=.false.
+         return
+      endif
+
+c     Also make sure there's no INF or NAN
+      do i=1,nexternal
+         do j=0,3
+            if(p(j,i).gt.1d32.or.p(j,i).ne.p(j,i))then
+               passcuts=.false.
+               return
+            endif
+         enddo
+      enddo
+      
+c
+c     Limit S_hat
+c
+c      if (x1*x2*stot .gt. 500**2) then
+c         passcuts=.false.
+c         return
+c      endif
+
+C $B$ DESACTIVATE_CUT $E$ !This is a tag for MadWeight
+
+      if(debug) write (*,*) '============================='
+      if(debug) write (*,*) ' EVENT STARTS TO BE CHECKED  '
+      if(debug) write (*,*) '============================='
+c     
+c     p_t min & max cuts
+c     
+      do i=nincoming+1,nexternal
+         if(debug) write (*,*) 'pt(',i,')=',pt(p(0,i)),'   ',etmin(i),
+     $        ':',etmax(i)
+         notgood=(pt(p(0,i)) .lt. etmin(i)).or.
+     &        (etmax(i).ge.0d0.and.pt(p(0,i)) .gt. etmax(i))
+         if (notgood) then
+            if(debug) write (*,*) i,' -> fails'
+            passcuts=.false.
+            return
+         endif
+      enddo
+c
+c    missing ET min & max cut + Invariant mass of leptons and neutrino 
+c    nb: missing Et defined as the vector sum over the neutrino's pt
+c
+c-- reset ptemp(0:3)
+      do j=0,3
+         ptemp(j)=0 ! for the neutrino
+         ptemp2(j)=0 ! for the leptons
+      enddo
+c-  sum over the momenta
+      do i=nincoming+1,nexternal
+         if(is_a_nu(i)) then            
+         if(debug) write (*,*) i,' -> neutrino '
+            do j=0,3
+               ptemp(j)=ptemp(j)+p(j,i)
+            enddo
+         elseif(is_a_l(i)) then            
+         if(debug) write (*,*) i,' -> lepton '
+            do j=0,3
+               ptemp2(j)=ptemp2(j)+p(j,i)
+            enddo
+         endif
+
+      enddo
+c-  check the et
+      if(debug.and.ptemp(0).eq.0d0) write (*,*) 'No et miss in event'
+      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Et miss =',pt(ptemp(0)),'   ',misset,':',missetmax
+      if(debug.and.ptemp2(0).eq.0d0) write (*,*) 'No leptons in event'
+      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Energy of leptons =',pt(ptemp2(0))
+      if(ptemp(0).gt.0d0) then
+         notgood=(pt(ptemp(0)) .lt. misset).or.
+     &        (missetmax.ge.0d0.and.pt(ptemp(0)) .gt. missetmax)
+         if (notgood) then
+            if(debug) write (*,*) ' missing et cut -> fails'
+            passcuts=.false.
+            return
+         endif
+      endif
+      if (mmnl.gt.0d0.or.mmnlmax.ge.0d0)then
+         if(dsqrt(SumDot(ptemp,ptemp2,1d0)).lt.mmnl.or.mmnlmax.ge.0d0.and.dsqrt(SumDot(ptemp, ptemp2,1d0)).gt.mmnlmax) then
+            if(debug) write (*,*) 'lepton invariant mass -> fails'
+            passcuts=.false.
+            return
+         endif
+      endif
+c
+c     pt cut on heavy particles
+c     gives min(pt) for (at least) one heavy particle
+c
+      if(ptheavy.gt.0d0)then
+         passcuts=.false.
+         foundheavy=.false.
+         do i=nincoming+1,nexternal
+            if(is_heavy(i)) then            
+               if(debug) write (*,*) i,' -> heavy '
+               foundheavy=.true.
+               if(pt(p(0,i)).gt.ptheavy) passcuts=.true.
+            endif
+         enddo
+         
+         if(.not.passcuts.and.foundheavy)then
+            if(debug) write (*,*) ' heavy particle cut -> fails'
+            return
+         else
+            passcuts=.true.
+         endif
+      endif
+c     
+c     E min & max cuts
+c     
+      do i=nincoming+1,nexternal
+         if(debug) write (*,*) 'p(0,',i,')=',p(0,i),'   ',emin(i),':',emax(i)
+         notgood=(p(0,i) .le. emin(i)).or.
+     &        (emax(i).ge.0d0 .and. p(0,i) .gt. emax(i))
+         if (notgood) then
+            if(debug) write (*,*) i,' -> fails'
+            passcuts=.false.
+            return
+         endif
+      enddo
+c     
+c     Rapidity  min & max cuts
+c     
+      do i=nincoming+1,nexternal
+         if(debug) write (*,*) 'abs(rap(',i,'))=',abs(rap(p(0,i))),'   ',etamin(i),':',etamax(i)
+         notgood=(etamax(i).ge.0.and.abs(rap(p(0,i))) .gt. etamax(i)).or.
+     &        (abs(rap(p(0,i))) .lt. etamin(i))
+         if (notgood) then
+            if(debug) write (*,*) i,' -> fails'
+            passcuts=.false.
+            return
+         endif
+      enddo
+c     
+c     DeltaR min & max cuts
+c     
+      do i=nincoming+1,nexternal-1
+         do j=i+1,nexternal
+            if(debug) write (*,*) 'r2(',i, ',' ,j,')=',dsqrt(r2(p(0,i),p(0,j)))
+            if(debug) write (*,*) dsqrt(r2min(j,i)),dsqrt(r2max(j,i))
+            if(r2min(j,i).gt.0.or.r2max(j,i).ge.0d0) then
+               tmp=r2(p(0,i),p(0,j))
+               notgood=(tmp .lt. r2min(j,i)).or.
+     $              (r2max(j,i).ge.0d0 .and. tmp .gt. r2max(j,i))
+               if (notgood) then
+                  if(debug) write (*,*) i,j,' -> fails'
+                  passcuts=.false.
+                  return
+               endif
+            endif
+         enddo
+      enddo
+
+
+c     s-channel min & max pt of sum of 4-momenta
+c     
+      do i=nincoming+1,nexternal-1
+         do j=i+1,nexternal
+            if(debug)write (*,*) 'ptll(',i,',',j,')=',PtDot(p(0,i),p(0,j))
+            if(debug)write (*,*) ptll_min(j,i),ptll_max(j,i)
+            if(ptll_min(j,i).gt.0.or.ptll_max(j,i).ge.0d0) then
+               tmp=PtDot(p(0,i),p(0,j))
+               notgood=(tmp .lt. ptll_min(j,i).or.
+     $              ptll_max(j,i).ge.0d0 .and. tmp.gt.ptll_max(j,i))
+               if (notgood) then
+                  if(debug) write (*,*) i,j,' -> fails'
+                  passcuts=.false.
+                  return
+               endif
+            endif
+         enddo
+      enddo
+
+
+
+
+c     
+c     s-channel min & max invariant mass cuts
+c     
+      do i=nincoming+1,nexternal-1
+         do j=i+1,nexternal
+            if(debug) write (*,*) 's(',i,',',j,')=',Sumdot(p(0,i),p(0,j),+1d0)
+            if(debug) write (*,*) s_min(j,i),s_max(j,i)
+            if(s_min(j,i).gt.0.or.s_max(j,i).ge.0d0) then
+               tmp=SumDot(p(0,i),p(0,j),+1d0)
+               if(s_min(j,i).le.s_max(j,i) .or. s_max(j,i).lt.0d0)then
+                  notgood=(tmp .lt. s_min(j,i).or.
+     $                 s_max(j,i).ge.0d0 .and. tmp .gt. s_max(j,i)) 
+                  if (notgood) then
+                     if(debug) write (*,*) i,j,' -> fails'
+                     passcuts=.false.
+                     return
+                  endif
+               else
+                  notgood=(tmp .lt. s_min(j,i).and.tmp .gt. s_max(j,i)) 
+                  if (notgood) then
+                     if(debug) write (*,*) i,j,' -> fails'
+                     passcuts=.false.
+                     return
+                  endif
+               endif
+            endif
+         enddo
+      enddo
+C     $B$DESACTIVATE_BW_CUT$B$ This is a Tag for MadWeight
+c     
+c     B.W. phase space cuts
+c     
+      pass_bw=cut_bw(p)
+c     JA 4/8/11 always check pass_bw
+      if ( pass_bw.and..not.CUTSPASSED) then
+         passcuts=.false.
+         if(debug) write (*,*) ' pass_bw -> fails'
+         return
+      endif
+C     $E$DESACTIVATE_BW_CUT$E$ This is a Tag for MadWeight
+        CUTSPASSED=.FALSE.
+C     
+C     maximal and minimal pt of the jets sorted by pt
+c     
+      njets=0
+      nheavyjets=0
+
+c- fill ptjet with the pt's of the jets.
+      do i=nincoming+1,nexternal
+         if(is_a_j(i)) then
+            njets=njets+1
+            ptjet(njets)=pt(p(0,i))
+         endif
+         if(is_a_b(i)) then
+            nheavyjets=nheavyjets+1
+            ptheavyjet(nheavyjets)=pt(p(0,i))
+         endif
+
+      enddo
+      if(debug) write (*,*) 'not yet ordered ',njets,'   ',ptjet
+
+C----------------------------------------------------------------------------
+C     DURHAM_KT CUT
+C----------------------------------------------------------------------------
+
+      IF ( KT_DURHAM .GT. 0D0) THEN
+
+C       RESET JET MOMENTA
+        njets=0
+        DO I=1,NEXTERNAL
+          DO J=0,3
+            PJET(I,J) = 0D0
+          ENDDO
+        ENDDO
+
+        do i=nincoming+1,nexternal
+           if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) ) then
+             njets=njets+1
+             DO J=0,3
+               PJET(NJETS,J) = P(J,I)
+             ENDDO
+           endif
+        enddo
+
+C       COUNT NUMBER OF MASSLESS OUTGOING PARTICLES, SINCE WE DO NOT WANT
+C       TO APPLY A CUT FOR JUST A SINGLE MASSIVE PARTICLE IN THE FINAL STATE.
+        NMASSLESS = 0
+        DO I=NINCOMING+1,NEXTERNAL
+          if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) .and. 
+     &                          DSQRT(DOT( P(0,I), P(0,I))) .LT. 5D0) THEN
+            NMASSLESS = NMASSLESS + 1
+          ENDIF
+        ENDDO
+
+C       DURHAM KT SEPARATION CUT
+        IF(NJETS.GT.0 .AND. NMASSLESS .GT. 0) THEN
+
+        PTMIN = EBEAM(1) + EBEAM(2)
+
+        DO I=1,NJETS
+
+C         PT WITH RESPECT TO Z AXIS FOR HADRONIC COLLISIONS
+          IF ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) THEN
+            PT1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2)
+            PTMIN = MIN( PTMIN, PT1 )
+          ENDIF
+
+          DO J=I+1,NJETS
+C           GET ANGLE BETWEEN JETS
+            PABS1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2 + PJET(I,3)**2)
+            PABS2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2 + PJET(J,3)**2)
+C           CHECK IF 3-MOMENTA DO NOT VANISH
+            IF(PABS1*PABS2 .NE. 0D0) THEN
+              COSTH = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) + PJET(I,3)*PJET(J,3) )/(PABS1*PABS2)
+            ELSE
+C           IF 3-MOMENTA VANISH, MAKE JET COSTH = 1D0 SO THAT JET MEASURE VANISHES
+              COSTH = 1D0
+            ENDIF
+C           GET PT AND ETA OF JETS
+            PT2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2)
+            ETA1 = 0.5D0*LOG( (PJET(I,0) + PJET(I,3)) / (PJET(I,0) - PJET(I,3)) )
+            ETA2 = 0.5D0*LOG( (PJET(J,0) + PJET(J,3)) / (PJET(J,0) - PJET(J,3)) )
+C           GET COSH OF DELTA ETA, COS OF DELTA PHI
+            COSH_DETA = DCOSH( ETA1 - ETA2 )
+            COS_DPHI = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) ) / (PT1*PT2)
+            DPHI = DACOS( COS_DPHI )
+            IF ( (LPP(1).EQ.0) .AND. (LPP(2).EQ.0)) THEN
+C             KT FOR E+E- COLLISION
+              PTNOW = DSQRT( 2D0*MIN(PJET(I,0)**2,PJET(J,0)**2)*( 1D0-COSTH ) )
+             ELSE
+C             HADRONIC KT, FASTJET DEFINITION
+              PTNOW = DSQRT( MIN(PT1**2,PT2**2)*( (ETA1 - ETA2 )**2 + DPHI**2 )/(D_PARAMETER**2) )
+            ENDIF
+
+            PTMIN = MIN( PTMIN, PTNOW )
+
+          ENDDO ! LOOP OVER NJET
+
+        ENDDO ! LOOP OVER NJET
+
+C       CHECK COMPATIBILITY WITH CUT
+        IF( (PTMIN .LT. KT_DURHAM)) THEN
+          PASSCUTS = .FALSE.
+          RETURN
+        ENDIF
+
+        ENDIF ! IF NJETS.GT. 0
+
+      ENDIF ! KT_DURHAM .GT. 0D0
+
+C----------------------------------------------------------------------------
+C     PTLUND CUT
+C----------------------------------------------------------------------------
+
+      IF(PT_LUND .GT. 0D0 ) THEN
+
+C       Reset jet momenta
+        NJETS=0
+        DO I=1,NEXTERNAL
+          JETFLAVOUR(I) = 0
+          DO J=0,3
+            PJET(I,J) = 0D0
+          ENDDO
+        ENDDO
+
+C       Fill incoming particle momenta
+        DO I=1,NINCOMING
+          INCFLAVOUR(I) = IDUP(I,1,1)
+          DO J=0,3
+            PINC(I,J) = P(J,I)
+          ENDDO
+        ENDDO
+
+C       Fill final jet momenta
+        DO I=NINCOMING+1,NEXTERNAL
+          if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) ) then
+            NJETS=NJETS+1
+            JETFLAVOUR(NJETS) = IDUP(I,1,1)
+            DO J=0,3
+              PJET(NJETS,J) = P(J,I)
+            ENDDO
+          ENDIF
+        ENDDO
+
+C       PROCESS WITH EXACTLY TWO MASSLESS OUTGOING PARTICLES IS SPECIAL
+C       BECAUSE AN ENERGY-SHARING VARIABLE LIKE "Z" DOES NOT MAKE SENSE.
+C       IN THIS CASE, ONLY APPLY MINIMAL pT W.R.T BEAM CUT.
+C       THIS CUT WILL ONLY APPLY TO THE TWO-MASSLESS PARTICLE STATE.
+        NMASSLESS = 0
+        DO I=NINCOMING+1,NEXTERNAL
+          if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) .and. 
+     &                          DSQRT(DOT( P(0,I), P(0,I))) .LT. 5D0) THEN
+            NMASSLESS = NMASSLESS + 1
+          ENDIF
+        ENDDO
+        IF (NMASSLESS .EQ. 2 .AND. NJETS .EQ. 2 .AND.
+     &                                     NEXTERNAL-NINCOMING .EQ. 2) THEN
+          PTMINSAVE = EBEAM(1) + EBEAM(2)
+          DO I=NINCOMING+1,NEXTERNAL
+            if( .not. from_decay(I) ) then  
+              PTMINSAVE = MIN(PTMINSAVE, PT(p(0,i)))
+            ENDIF
+          ENDDO
+C         CHECK COMPATIBILITY WITH CUT
+          IF ( ((LPP(1).NE.0) .OR. (LPP(2).NE.0)) .AND.
+     &                                         PTMINSAVE .LT. PT_LUND) THEN
+            PASSCUTS = .FALSE.
+            RETURN
+          ENDIF
+C         RESET NJETS TO AVOID FURTHER MERGING SCALE CUT.
+          NJETS=0
+        ENDIF
+
+C       PYTHIA PT SEPARATION CUT
+        IF(NJETS.GT.0 .AND. NMASSLESS .GT. 0) THEN
+
+        PTMINSAVE = EBEAM(1) + EBEAM(2)
+
+        DO I=1,NJETS
+
+          PTMIN = EBEAM(1) + EBEAM(2)
+          PTMINSAVE = MIN(PTMIN,PTMINSAVE)
+
+C         Compute pythia ISR separation between i-jet and incoming.
+C         Only SM-like emissions off the beam are possible.
+          IF ( ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) .AND. 
+     &                               ABS(JETFLAVOUR(I)) .LT. 30 ) THEN
+C           Check separation to first incoming particle
+            DO L=0,3
+              PRADTEMP(L) = PINC(1,L)
+              PEMTTEMP(L) = PJET(I,L)
+              PRECTEMP(L) = PINC(2,L)
+            ENDDO
+            PT1 = RHOPYTHIA(PRADTEMP, PEMTTEMP, PRECTEMP, INCFLAVOUR(1),
+     &                                            JETFLAVOUR(I), -1, -1)
+            PTMIN = MIN( PTMIN, PT1 )
+C           Check separation to second incoming particle
+            DO L=0,3
+              PRADTEMP(L) = PINC(2,L)
+              PEMTTEMP(L) = PJET(I,L)
+              PRECTEMP(L) = PINC(1,L)
+            ENDDO
+            PT2 = RHOPYTHIA(PRADTEMP, PEMTTEMP, PRECTEMP, INCFLAVOUR(2),
+     &                                            JETFLAVOUR(I), -1, -1)
+            PTMIN = MIN( PTMIN, PT2 )
+          ENDIF
+
+C         Compute pythia FSR separation between two jets,
+C         without any knowledge of colour connections
+          DO J=1,NJETS
+            DO K=1,NJETS
+              IF ( I .NE. J .AND. I .NE. K .AND. J .NE. K ) THEN
+
+C               Check separation between final partons i and j, with k as spectator
+                DO L=0,3
+                  PRADTEMP(L) = PJET(J,L)
+                  PEMTTEMP(L) = PJET(I,L)
+                  PRECTEMP(L) = PJET(K,L)
+                ENDDO
+
+                TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP,
+     &                               JETFLAVOUR(J), JETFLAVOUR(I), 1, 1)
+C               Only SM-like emissions off the beam are possible, no additional
+C               BSM particles will be produced as as shower emissions.
+                IF ( ABS(JETFLAVOUR(I)) .LT. 30 ) THEN
+                  PTMIN = MIN(PTMIN, TEMP);
+                ENDIF
+
+                TEMP = RHOPYTHIA( PEMTTEMP, PRADTEMP, PRECTEMP,
+     &                               JETFLAVOUR(I), JETFLAVOUR(J), 1, 1)
+C               Only SM-like emissions off the beam are possible, no additional
+C               BSM particles will be produced as as shower emissions.
+                IF ( ABS(JETFLAVOUR(J)) .LT. 30 ) THEN
+                  PTMIN = MIN(PTMIN, TEMP);
+                ENDIF
+
+              ENDIF
+
+            ENDDO ! LOOP OVER NJET
+          ENDDO ! LOOP OVER NJET
+
+C         Compute pythia FSR separation between two jets, with initial spectator
+          IF ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) THEN
+            DO J=1,NJETS
+
+C             BSM particles can only be radiators, and will not be produced
+C             as shower emissions.
+              IF ( ABS(JETFLAVOUR(I)) .GT. 1000000 ) THEN
+                EXIT
+              ENDIF
+
+C             Allow both initial partons as recoiler
+              IF ( I .NE. J ) THEN
+
+C               Check with first initial as recoiler
+                DO L=0,3
+                  PRADTEMP(L) = PJET(J,L)
+                  PEMTTEMP(L) = PJET(I,L)
+                  PRECTEMP(L) = PINC(1,L)
+                ENDDO
+                TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP, 
+     &                             JETFLAVOUR(J), JETFLAVOUR(I), 1, -1);
+                IF ( LPP(1) .NE. 0 ) THEN
+                  PTMIN = MIN(PTMIN, TEMP);
+                ENDIF
+                DO L=0,3
+                  PRADTEMP(L) = PJET(J,L)
+                  PEMTTEMP(L) = PJET(I,L)
+                  PRECTEMP(L) = PINC(2,L)
+                ENDDO
+                TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP,
+     &                             JETFLAVOUR(J), JETFLAVOUR(I), 1, -1);
+                IF ( LPP(2) .NE. 0 ) THEN
+                  PTMIN = MIN(PTMIN, TEMP);
+                ENDIF
+              ENDIF
+            ENDDO ! LOOP OVER NJET
+          ENDIF
+
+          PTMINSAVE = MIN(PTMIN,PTMINSAVE)
+
+        ENDDO ! LOOP OVER NJET
+
+C       CHECK COMPATIBILITY WITH CUT
+        IF (PTMINSAVE .LT. PT_LUND) THEN
+          PASSCUTS = .FALSE.
+          RETURN
+        ENDIF
+ 
+        ENDIF ! IF NJETS.GT. 0 
+
+      ENDIF ! PT_LUND .GT. 0D0 
+
+C----------------------------------------------------------------------------
+C----------------------------------------------------------------------------
+
+
+
+
+c- check existance of jets if jet cuts are on
+      if(njets.lt.1.and.(htjmin.gt.0.or.ptj1min.gt.0).or.
+     $     njets.lt.2.and.ptj2min.gt.0.or.
+     $     njets.lt.3.and.ptj3min.gt.0.or.
+     $     njets.lt.4.and.ptj4min.gt.0)then
+         if(debug) write (*,*) i, ' too few jets -> fails'
+         passcuts=.false.
+         return
+      endif
+
+c - sort jet pts
+      do i=1,njets-1
+         do j=i+1,njets
+            if(ptjet(j).gt.ptjet(i)) then
+               temp=ptjet(i)
+               ptjet(i)=ptjet(j)
+               ptjet(j)=temp
+            endif
+         enddo
+      enddo
+      if(debug) write (*,*) 'ordered ',njets,'   ',ptjet
+c
+c     Use "and" or "or" prescriptions 
+c     
+      inclht=0
+
+      if(njets.gt.0) then
+
+       notgood=.not.jetor
+       if(debug) write (*,*) 'jetor :',jetor  
+       if(debug) write (*,*) '0',notgood   
+      
+      do i=1,min(njets,4)
+            if(debug) write (*,*) i,ptjet(i), '   ',ptjmin4(i),
+     $        ':',ptjmax4(i)
+         if(jetor) then     
+c---  if one of the jets does not pass, the event is rejected
+            notgood=notgood.or.(ptjmax4(i).ge.0d0 .and.
+     $           ptjet(i).gt.ptjmax4(i)) .or.
+     $           (ptjet(i).lt.ptjmin4(i))
+            if(debug) write (*,*) i,' notgood total:', notgood   
+         else
+c---  all cuts must fail to reject the event
+            notgood=notgood.and.(ptjmax4(i).ge.0d0 .and.
+     $              ptjet(i).gt.ptjmax4(i) .or.
+     $              (ptjet(i).lt.ptjmin4(i)))
+            if(debug) write (*,*) i,' notgood total:', notgood   
+         endif
+      enddo
+
+
+      if (notgood) then
+         if(debug) write (*,*) i, ' multiple pt -> fails'
+         passcuts=.false.
+         return
+      endif
+
+c---------------------------
+c      Ht cuts
+C---------------------------
+      htj=ptjet(1)
+
+      do i=2,njets
+         htj=htj+ptjet(i)
+         if(debug) write (*,*) i, 'htj ',htj
+         if(debug.and.i.le.4) write (*,*) 'htmin ',i,' ', htjmin4(i),':',htjmax4(i)
+         if(i.le.4)then
+            if(htj.lt.htjmin4(i) .or.
+     $        htjmax4(i).ge.0d0.and.htj.gt.htjmax4(i)) then
+            if(debug) write (*,*) i, ' ht -> fails'
+            passcuts=.false.
+            return
+            endif
+         endif
+      enddo
+
+      if(htj.lt.htjmin.or.htjmax.ge.0d0.and.htj.gt.htjmax)then
+         if(debug) write (*,*) i, ' htj -> fails'
+         passcuts=.false.
+         return
+      endif
+
+      inclht=htj
+
+      endif !if there are jets 
+
+      if(nheavyjets.gt.0) then
+         do i=1,nheavyjets
+            inclht=inclht+ptheavyjet(i)
+         enddo
+      endif !if there are heavyjets
+
+      if(inclht.lt.inclHtmin.or.
+     $     inclHtmax.ge.0d0.and.inclht.gt.inclHtmax)then
+         if(debug) write (*,*) ' inclhtmin=',inclHtmin,' -> fails'
+         passcuts=.false.
+         return
+      endif
+
+C     
+C     maximal and minimal pt of the leptons sorted by pt
+c     
+      nleptons=0
+
+      if(ptl1min.gt.0.or.ptl2min.gt.0.or.ptl3min.gt.0.or.ptl4min.gt.0.or.
+     $     ptl1max.ge.0d0.or.ptl2max.ge.0d0.or.
+     $     ptl3max.ge.0d0.or.ptl4max.ge.0d0) then
+
+c     - fill ptlepton with the pt's of the leptons.
+         do i=nincoming+1,nexternal
+            if(is_a_l(i)) then
+               nleptons=nleptons+1
+               ptlepton(nleptons)=pt(p(0,i))
+            endif
+         enddo
+         if(debug) write (*,*) 'not yet ordered ',njets,'   ',ptjet
+
+c     - check existance of leptons if lepton cuts are on
+         if(nleptons.lt.1.and.ptl1min.gt.0.or.
+     $        nleptons.lt.2.and.ptl2min.gt.0.or.
+     $        nleptons.lt.3.and.ptl3min.gt.0.or.
+     $        nleptons.lt.4.and.ptl4min.gt.0)then
+            if(debug) write (*,*) i, ' too few leptons -> fails'
+            passcuts=.false.
+            return
+         endif
+
+c     - sort lepton pts
+         do i=1,nleptons-1
+            do j=i+1,nleptons
+               if(ptlepton(j).gt.ptlepton(i)) then
+                  temp=ptlepton(i)
+                  ptlepton(i)=ptlepton(j)
+                  ptlepton(j)=temp
+               endif
+            enddo
+         enddo
+         if(debug) write (*,*) 'ordered ',nleptons,'   ',ptlepton
+
+         if(nleptons.gt.0) then
+
+            notgood = .false.
+            do i=1,min(nleptons,4)
+               if(debug) write (*,*) i,ptlepton(i), '   ',ptlmin4(i),':',ptlmax4(i)
+c---  if one of the leptons does not pass, the event is rejected
+               notgood=notgood.or.
+     $              (ptlmax4(i).ge.0d0.and.ptlepton(i).gt.ptlmax4(i)).or.
+     $              (ptlepton(i).lt.ptlmin4(i))
+               if(debug) write (*,*) i,' notgood total:', notgood   
+            enddo
+
+
+            if (notgood) then
+               if(debug) write (*,*) i, ' multiple pt -> fails'
+               passcuts=.false.
+               return
+            endif
+         endif
+      endif
+C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>     
+C     SPECIAL CUTS
+C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+C     REQUIRE AT LEAST ONE JET WITH PT>XPTJ
+         
+         IF(xptj.gt.0d0) THEN
+            xvar=0
+            do i=nincoming+1,nexternal
+               if(is_a_j(i)) xvar=max(xvar,pt(p(0,i)))
+            enddo
+            if (xvar .lt. xptj) then
+               passcuts=.false.
+               return
+            endif
+         ENDIF
+
+C     REQUIRE AT LEAST ONE PHOTON WITH PT>XPTA
+         
+         IF(xpta.gt.0d0) THEN
+            xvar=0
+            do i=nincoming+1,nexternal
+               if(is_a_a(i)) xvar=max(xvar,pt(p(0,i)))
+            enddo
+            if (xvar .lt. xpta) then
+               passcuts=.false.
+               return
+            endif
+         ENDIF
+
+C     REQUIRE AT LEAST ONE B  WITH PT>XPTB
+         
+         IF(xptb.gt.0d0) THEN
+            xvar=0
+            do i=nincoming+1,nexternal
+               if(is_a_b(i)) xvar=max(xvar,pt(p(0,i)))
+            enddo
+            if (xvar .lt. xptb) then
+               passcuts=.false.
+               return
+            endif
+         ENDIF
+
+C     REQUIRE AT LEAST ONE LEPTON  WITH PT>XPTL
+         
+         IF(xptl.gt.0d0) THEN
+            xvar=0
+            do i=nincoming+1,nexternal
+               if(is_a_l(i)) xvar=max(xvar,pt(p(0,i)))
+            enddo
+            if (xvar .lt. xptl) then
+               passcuts=.false.
+               if(debug) write (*,*) ' xptl -> fails'
+               return
+            endif
+         ENDIF
+C
+C     WBF CUTS: TWO TYPES
+C    
+C     FIRST TYPE:  implemented by FM
+C
+C     1. FIND THE 2 HARDEST JETS
+C     2. REQUIRE |RAP(J)|>XETAMIN
+C     3. REQUIRE RAP(J1)*ETA(J2)<0
+C
+C     SECOND TYPE : added by Simon de Visscher 1-08-2007
+C
+C     1. FIND THE 2 HARDEST JETS
+C     2. REQUIRE |RAP(J1)-RAP(J2)|>DELTAETA
+C     3. REQUIRE RAP(J1)*RAP(J2)<0
+C
+C
+         hardj1=0
+         hardj2=0
+         ptmax1=0d0
+         ptmax2=0d0
+
+C-- START IF AT LEAST ONE OF THE CUTS IS ACTIVATED
+         
+         IF(XETAMIN.GT.0D0.OR.DELTAETA.GT.0D0) THEN
+            
+C-- FIND THE HARDEST JETS
+
+            do i=nincoming+1,nexternal
+               if(is_a_j(i)) then
+c                  write (*,*) i,pt(p(0,i))
+                  if(pt(p(0,i)).gt.ptmax1) then
+                     hardj2=hardj1
+                     ptmax2=ptmax1
+                     hardj1=i
+                     ptmax1=pt(p(0,i))
+                  elseif(pt(p(0,i)).gt.ptmax2) then
+                     hardj2=i
+                     ptmax2=pt(p(0,i))
+                  endif
+c                  write (*,*) hardj1,hardj2,ptmax1,ptmax2
+               endif
+            enddo
+            
+            if (hardj2.eq.0) goto 21 ! bypass vbf cut since not enough jets
+
+C-- NOW APPLY THE CUT I            
+
+            if (abs(rap(p(0,hardj1))) .lt. xetamin
+     &       .or.abs(rap(p(0,hardj2))) .lt. xetamin
+     &       .or.rap(p(0,hardj1))*rap(p(0,hardj2)) .gt.0d0) then
+             passcuts=.false.
+             return
+            endif
+
+            
+C-- NOW APPLY THE CUT II
+            
+            if (abs(rap(p(0,hardj1))-rap(p(0,hardj2))) .lt. deltaeta) then
+             passcuts=.false.
+             return
+            endif
+         
+c            write (*,*) hardj1,hardj2,rap(p(0,hardj1)),rap(p(0,hardj2))
+         
+         ENDIF
+
+c Begin photon isolation
+c NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE
+c     Use is made of parton cm frame momenta. If this must be
+c     changed, pQCD used below must be redefined
+c NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE
+c If we do not require a mimimum jet energy, there's no need to apply
+c jet clustering and all that.
+      if (ptgmin.ne.0d0) then
+
+c Put all (light) QCD partons in momentum array for jet clustering.
+c From the run_card.dat, maxjetflavor defines if b quark should be
+c considered here (via the logical variable 'is_a_jet').  nQCD becomes
+c the number of (light) QCD partons at the real-emission level (i.e. one
+c more than the Born).
+
+      nQCD=0
+      do j=nincoming+1,nexternal
+         if (is_a_j(j)) then
+            nQCD=nQCD+1
+            do i=0,3
+               pQCD(i,nQCD)=p(i,j) ! Use C.o.M. frame momenta
+            enddo
+         endif
+      enddo
+
+        nph=0
+        do j=nincoming+1,nexternal
+          if (is_a_a(j)) then
+            nph=nph+1
+            do i=0,3
+              pgamma(i,nph)=p(i,j) ! Use C.o.M. frame momenta
+            enddo
+          endif
+        enddo
+        if(nph.eq.0) goto 444
+
+        if(isoEM)then
+          nem=nph
+          do k=1,nem
+            do i=0,3
+              pem(i,k)=pgamma(i,k)
+            enddo
+          enddo
+          do j=nincoming+1,nexternal
+            if (is_a_l(j)) then
+              nem=nem+1
+              do i=0,3
+                pem(i,nem)=p(i,j) ! Use C.o.M. frame momenta
+              enddo
+            endif
+          enddo
+        endif
+
+        alliso=.true.
+
+        j=0
+        dowhile(j.lt.nph.and.alliso)
+c Loop over all photons
+          j=j+1
+
+          ptg=pt(pgamma(0,j))
+          if(ptg.lt.ptgmin)then
+            passcuts=.false.
+            return
+          endif
+
+c Isolate from hadronic energy
+          do i=1,nQCD
+            drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pQCD(0,i)))
+          enddo
+          call sortzv(drlist,isorted,nQCD,ismode,isway,izero)
+          Etsum(0)=0.d0
+          nin=0
+          do i=1,nQCD
+            if(dble(drlist(isorted(i))).le.R0gamma)then
+              nin=nin+1
+              Etsum(nin)=Etsum(nin-1)+pt(pQCD(0,isorted(i)))
+            endif
+          enddo
+          do i=1,nin
+            alliso=alliso .and.
+     #        Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
+     #                                  R0gamma,xn,epsgamma,ptg)
+          enddo
+
+c Isolate from EM energy
+          if(isoEM.and.nem.gt.1)then
+            do i=1,nem
+              drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i)))
+            enddo
+            call sortzv(drlist,isorted,nem,ismode,isway,izero)
+c First of list must be the photon: check this, and drop it
+            if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then
+              write(*,*)'Error #1 in photon isolation'
+              write(*,*)j,isorted(1),drlist(isorted(1))
+              stop
+            endif
+            Etsum(0)=0.d0
+            nin=0
+            do i=2,nem
+              if(dble(drlist(isorted(i))).le.R0gamma)then
+                nin=nin+1
+                Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i)))
+              endif
+            enddo
+            do i=1,nin
+              alliso=alliso .and.
+     #          Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
+     #                                    R0gamma,xn,epsgamma,ptg)
+            enddo
+
+          endif
+
+c End of loop over photons
+        enddo
+
+        if(.not.alliso)then
+          passcuts=.false.
+          return
+        endif
+      endif
+
+ 444    continue
+c End photon isolation
+
+
+C...Set couplings if event passed cuts
+
+ 21   if(.not.fixed_ren_scale) then
+         call set_ren_scale(P,scale)
+         if(scale.gt.0) G = SQRT(4d0*PI*ALPHAS(scale))
+      endif
+
+      if(.not.fixed_fac_scale) then
+         call set_fac_scale(P,q2fact)
+      endif
+
+c
+c     Here we cluster event and reset factorization and renormalization
+c     scales on an event-by-event basis, as well as check xqcut for jets
+c
+c     Note the following condition is the first line of setclscales
+c      if(xqcut.gt.0d0.or.ickkw.gt.0.or.scale.eq.0.or.q2fact(1).eq.0)then
+c     Do not duplicate it since some variable are set for syscalc in the fct
+        if(.not.setclscales(p))then
+           cutsdone=.false.
+           cutspassed=.false.
+           passcuts = .false.
+           if(debug) write (*,*) 'setclscales -> fails'
+           return
+       endif
+c      endif
+
+c     Set couplings in model files
+      if(.not.fixed_ren_scale.or..not.fixed_couplings) then
+         if (.not.fixed_couplings)then
+            do i=0,3
+               do j=1,nexternal
+                  pp(i,j)=p(i,j)
+               enddo
+            enddo
+         endif
+         call update_as_param()
+      endif
+
+      IF (FIRSTTIME2) THEN
+        FIRSTTIME2=.FALSE.
+        write(6,*) 'alpha_s for scale ',scale,' is ', G**2/(16d0*atan(1d0))
+      ENDIF
+
+      if(debug) write (*,*) '============================='
+      if(debug) write (*,*) ' EVENT PASSED THE CUTS       '
+      if(debug) write (*,*) '============================='
+
+      CUTSPASSED=.TRUE.
+
+      RETURN
+      END
+
+
+C
+C     FUNCTION FOR ISOLATION
+C
+
+      function iso_getdrv40(p1,p2)
+      implicit none
+      real*8 iso_getdrv40,p1(0:3),p2(0:3)
+      real*8 iso_getdr
+c
+      iso_getdrv40=iso_getdr(p1(0),p1(1),p1(2),p1(3),
+     #                       p2(0),p2(1),p2(2),p2(3))
+      return
+      end
+
+
+      function iso_getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
+      implicit none
+      real*8 iso_getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
+     # iso_getpseudorap,iso_getdelphi
+c
+      deta=iso_getpseudorap(en1,ptx1,pty1,pl1)-
+     #     iso_getpseudorap(en2,ptx2,pty2,pl2)
+      dphi=iso_getdelphi(ptx1,pty1,ptx2,pty2)
+      iso_getdr=sqrt(dphi**2+deta**2)
+      return
+      end
+
+
+      function iso_getpseudorap(en,ptx,pty,pl)
+      implicit none
+      real*8 iso_getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
+      parameter (tiny=1.d-5)
+c
+      pt=sqrt(ptx**2+pty**2)
+      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
+        eta=sign(1.d0,pl)*1.d8
+      else
+        th=atan2(pt,pl)
+        eta=-log(tan(th/2.d0))
+      endif
+      iso_getpseudorap=eta
+      return
+      end
+
+
+      function iso_getdelphi(ptx1,pty1,ptx2,pty2)
+      implicit none
+      real*8 iso_getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
+      parameter (tiny=1.d-5)
+c
+      pt1=sqrt(ptx1**2+pty1**2)
+      pt2=sqrt(ptx2**2+pty2**2)
+      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
+        tmp=ptx1*ptx2+pty1*pty2
+        tmp=tmp/(pt1*pt2)
+        if(abs(tmp).gt.1.d0+tiny)then
+          write(*,*)'Cosine larger than 1'
+          stop
+        elseif(abs(tmp).ge.1.d0)then
+          tmp=sign(1.d0,tmp)
+        endif
+        tmp=acos(tmp)
+      else
+        tmp=1.d8
+      endif
+      iso_getdelphi=tmp
+      return
+      end
+
+      function chi_gamma_iso(dr,R0,xn,epsgamma,pTgamma)
+c Eq.(3.4) of Phys.Lett. B429 (1998) 369-374 [hep-ph/9801442]
+      implicit none
+      real*8 chi_gamma_iso,dr,R0,xn,epsgamma,pTgamma
+      real*8 tmp,axn
+c
+      axn=abs(xn)
+      tmp=epsgamma*pTgamma
+      if(axn.ne.0.d0)then
+        tmp=tmp*( (1-cos(dr))/(1-cos(R0)) )**axn
+      endif
+      chi_gamma_iso=tmp
+      return
+      end
+
+
+*
+* $Id: sortzv.F,v 1.1.1.1 1996/02/15 17:49:50 mclareni Exp $
+*
+* $Log: sortzv.F,v $
+* Revision 1.1.1.1  1996/02/15 17:49:50  mclareni
+* Kernlib
+*
+*
+c$$$#include "kerngen/pilot.h"
+      SUBROUTINE SORTZV (A,INDEX,N1,MODE,NWAY,NSORT)
+C
+C CERN PROGLIB# M101    SORTZV          .VERSION KERNFOR  3.15  820113
+C ORIG. 02/10/75
+C
+      DIMENSION A(N1),INDEX(N1)
+C
+C
+      N = N1
+      IF (N.LE.0)            RETURN
+      IF (NSORT.NE.0) GO TO 2
+      DO 1 I=1,N
+    1 INDEX(I)=I
+C
+    2 IF (N.EQ.1)            RETURN
+      IF (MODE)    10,20,30
+   10 CALL SORTTI (A,INDEX,N)
+      GO TO 40
+C
+   20 CALL SORTTC(A,INDEX,N)
+      GO TO 40
+C
+   30 CALL SORTTF (A,INDEX,N)
+C
+   40 IF (NWAY.EQ.0) GO TO 50
+      N2 = N/2
+      DO 41 I=1,N2
+      ISWAP = INDEX(I)
+      K = N+1-I
+      INDEX(I) = INDEX(K)
+   41 INDEX(K) = ISWAP
+   50 RETURN
+      END
+*     ========================================
+      SUBROUTINE SORTTF (A,INDEX,N1)
+C
+      DIMENSION A(N1),INDEX(N1)
+C
+      N = N1
+      DO 3 I1=2,N
+      I3 = I1
+      I33 = INDEX(I3)
+      AI = A(I33)
+    1 I2 = I3/2
+      IF (I2) 3,3,2
+    2 I22 = INDEX(I2)
+      IF (AI.LE.A (I22)) GO TO 3
+      INDEX (I3) = I22
+      I3 = I2
+      GO TO 1
+    3 INDEX (I3) = I33
+    4 I3 = INDEX (N)
+      INDEX (N) = INDEX (1)
+      AI = A(I3)
+      N = N-1
+      IF (N-1) 12,12,5
+    5 I1 = 1
+    6 I2 = I1 + I1
+      IF (I2.LE.N) I22= INDEX(I2)
+      IF (I2-N) 7,9,11
+    7 I222 = INDEX (I2+1)
+      IF (A(I22)-A(I222)) 8,9,9
+    8 I2 = I2+1
+      I22 = I222
+    9 IF (AI-A(I22)) 10,11,11
+   10 INDEX(I1) = I22
+      I1 = I2
+      GO TO 6
+   11 INDEX (I1) = I3
+      GO TO 4
+   12 INDEX (1) = I3
+      RETURN
+      END
+*     ========================================
+      SUBROUTINE SORTTI (A,INDEX,N1)
+C
+      INTEGER A,AI
+      DIMENSION A(N1),INDEX(N1)
+C
+      N = N1
+      DO 3 I1=2,N
+      I3 = I1
+      I33 = INDEX(I3)
+      AI = A(I33)
+    1 I2 = I3/2
+      IF (I2) 3,3,2
+    2 I22 = INDEX(I2)
+      IF (AI.LE.A (I22)) GO TO 3
+      INDEX (I3) = I22
+      I3 = I2
+      GO TO 1
+    3 INDEX (I3) = I33
+    4 I3 = INDEX (N)
+      INDEX (N) = INDEX (1)
+      AI = A(I3)
+      N = N-1
+      IF (N-1) 12,12,5
+    5 I1 = 1
+    6 I2 = I1 + I1
+      IF (I2.LE.N) I22= INDEX(I2)
+      IF (I2-N) 7,9,11
+    7 I222 = INDEX (I2+1)
+      IF (A(I22)-A(I222)) 8,9,9
+    8 I2 = I2+1
+      I22 = I222
+    9 IF (AI-A(I22)) 10,11,11
+   10 INDEX(I1) = I22
+      I1 = I2
+      GO TO 6
+   11 INDEX (I1) = I3
+      GO TO 4
+   12 INDEX (1) = I3
+      RETURN
+      END
+*     ========================================
+      SUBROUTINE SORTTC (A,INDEX,N1)
+C
+      INTEGER A,AI
+      DIMENSION A(N1),INDEX(N1)
+C
+      N = N1
+      DO 3 I1=2,N
+      I3 = I1
+      I33 = INDEX(I3)
+      AI = A(I33)
+    1 I2 = I3/2
+      IF (I2) 3,3,2
+    2 I22 = INDEX(I2)
+      IF(ICMPCH(AI,A(I22)))3,3,21
+   21 INDEX (I3) = I22
+      I3 = I2
+      GO TO 1
+    3 INDEX (I3) = I33
+    4 I3 = INDEX (N)
+      INDEX (N) = INDEX (1)
+      AI = A(I3)
+      N = N-1
+      IF (N-1) 12,12,5
+    5 I1 = 1
+    6 I2 = I1 + I1
+      IF (I2.LE.N) I22= INDEX(I2)
+      IF (I2-N) 7,9,11
+    7 I222 = INDEX (I2+1)
+      IF (ICMPCH(A(I22),A(I222))) 8,9,9
+    8 I2 = I2+1
+      I22 = I222
+    9 IF (ICMPCH(AI,A(I22))) 10,11,11
+   10 INDEX(I1) = I22
+      I1 = I2
+      GO TO 6
+   11 INDEX (I1) = I3
+      GO TO 4
+   12 INDEX (1) = I3
+      RETURN
+      END
+*     ========================================
+      FUNCTION ICMPCH(IC1,IC2)
+C     FUNCTION TO COMPARE TWO 4 CHARACTER EBCDIC STRINGS - IC1,IC2
+C     ICMPCH=-1 IF HEX VALUE OF IC1 IS LESS THAN IC2
+C     ICMPCH=0  IF HEX VALUES OF IC1 AND IC2 ARE THE SAME
+C     ICMPCH=+1 IF HEX VALUES OF IC1 IS GREATER THAN IC2
+      I1=IC1
+      I2=IC2
+      IF(I1.GE.0.AND.I2.GE.0)GOTO 40
+      IF(I1.GE.0)GOTO 60
+      IF(I2.GE.0)GOTO 80
+      I1=-I1
+      I2=-I2
+      IF(I1-I2)80,70,60
+ 40   IF(I1-I2)60,70,80
+ 60   ICMPCH=-1
+      RETURN
+ 70   ICMPCH=0
+      RETURN
+ 80   ICMPCH=1
+      RETURN
+      END
+
+c************************************************************************
+c     Returns pTLund (i.e. the Pythia8 evolution variable) between two 
+c     particles with momenta prad and pemt with momentum prec as spectator
+c************************************************************************
+
+      DOUBLE PRECISION FUNCTION RHOPYTHIA(PRAD, PEMT, PREC, FLAVRAD,
+     &                                        FLAVEMT, RADTYPE, RECTYPE)
+
+       IMPLICIT NONE
+c
+c     Arguments
+c
+      DOUBLE PRECISION PRAD(0:3),PEMT(0:3), PREC(0:3)
+      INTEGER FLAVRAD, FLAVEMT, RADTYPE,RECTYPE
+c
+c     Local
+c
+      DOUBLE PRECISION Q(0:3),SUM(0:3), qBR(0:3), qAR(0:3)
+      DOUBLE PRECISION Q2, m2Rad, m2Emt, m2Dip, qBR2, qAR2, x1, x2, z
+      DOUBLE PRECISION m2RadAft, m2EmtAft, m2Rec, m2RadBef, m2ar, rescale
+      DOUBLE PRECISION TEMP, lambda13, k1, k3, m2Final
+      DOUBLE PRECISION m0u, m0d, m0c, m0s, m0t, m0b, m0w, m0z, m0x
+      DOUBLE PRECISION PRECAFT(0:3)
+      INTEGER emtsign
+      INTEGER idRadBef
+      LOGICAL allowed
+
+c-----
+c  Begin Code
+c-----
+
+C     Set masses. Currently no way of getting those?
+      m0u = 0.0
+      m0d = 0.0
+      m0c = 1.5
+      m0s = 0.0
+      m0t = 172.5
+      m0w = 80.4
+      m0z = 91.188
+      m0x = 400.0
+
+C     Store recoiler momentum (since FI splittings require recoiler
+C     rescaling)
+      PRECAFT(0) = PREC(0)
+      PRECAFT(1) = PREC(1)
+      PRECAFT(2) = PREC(2)
+      PRECAFT(3) = PREC(3)
+C     Get sign of emitted momentum
+      emtsign = 1
+      if(radtype .eq. -1) emtsign = -1
+
+C     Get virtuality
+      Q(0) = pRad(0) + emtsign*pEmt(0)
+      Q(1) = pRad(1) + emtsign*pEmt(1)
+      Q(2) = pRad(2) + emtsign*pEmt(2)
+      Q(3) = pRad(3) + emtsign*pEmt(3)
+      Q2   = emtsign * ( Q(0)**2 - Q(1)**2 - Q(2)**2 - Q(3)**2 );
+
+C     Reset allowed
+      allowed = .true.
+
+C     Splitting not possible for negative virtuality.
+      if ( Q2 .lt. 0.0 ) allowed = .false.
+
+C     Try to reconstruct flavour of radiator before emission.
+      idRadBef = 0
+C     gluon radiation: idBef = idAft
+      if (abs(flavEmt) .eq. 21 .or. abs(flavEmt) .eq. 22 ) idRadBef=flavRad
+C     final state gluon splitting: idBef = 21
+      if (radtype .eq.  1 .and. flavEmt .eq. -flavRad) idRadBef=21
+C     final state quark -> gluon conversion
+      if (radtype .eq.  1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. 21) idRadBef=flavEmt
+C     initial state gluon splitting: idBef = -idEmt
+      if (radtype .eq. -1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. 21) idRadBef=-flavEmt
+C     initial state gluon -> quark conversion
+      if (radtype .eq. -1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. flavEmt) idRadBef=21
+C     W-boson radiation
+      if (flavEmt .eq.  24) idRadBef = flavRad+1
+      if (flavEmt .eq. -24) idRadBef = flavRad-1
+
+C     Set particle masses.
+      m2RadAft = 0.0
+      m2EmtAft = 0.0
+      m2Rec    = 0.0
+      m2RadBef = 0.0
+
+      m2RadAft = pRad(0)*pRad(0)-pRad(1)*pRad(1)-pRad(2)*pRad(2)-pRad(3)*pRad(3)
+      m2EmtAft = pEmt(0)*pEmt(0)-pEmt(1)*pEmt(1)-pEmt(2)*pEmt(2)-pEmt(3)*pEmt(3)
+      m2Rec    = pRec(0)*pRec(0)-pRec(1)*pRec(1)-pRec(2)*pRec(2)-pRec(3)*pRec(3)
+
+      if (m2RadAft .lt. 1d-4) m2RadAft = 0.0
+      if (m2EmtAft .lt. 1d-4) m2EmtAft = 0.0
+      if (m2Rec    .lt. 1d-4) m2Rec    = 0.0
+
+      if (abs(flavRad) .ne. 21 .and. abs(flavRad) .ne. 22 .and.
+     &    abs(flavEmt) .ne. 24 .and.
+     &    abs(flavRad) .ne. abs(flavEmt)) then
+        m2RadBef = m2RadAft
+      else if (abs(flavEmt) .eq. 24) then
+        if (idRadBef .ne. 0) then
+          if( abs(idRadBef) .eq. 4 ) m2RadBef       = m0c**2
+          if( abs(idRadBef) .eq. 5 ) m2RadBef       = m0b**2
+          if( abs(idRadBef) .eq. 6 ) m2RadBef       = m0t**2
+          if( abs(idRadBef) .eq. 9000001 ) m2RadBef = m0x**2
+        endif
+      else if (radtype .eq. -1) then
+        if (abs(flavRad) .eq. 21 .and. abs(flavEmt) .eq. 21) m2RadBef = m2EmtAft
+      endif
+
+C     Calculate dipole mass for final-state radiation.
+      m2Final = 0.0
+      m2Final = m2Final + (pRad(0) + pRec(0) + pEmt(0))**2
+      m2Final = m2Final - (pRad(1) + pRec(1) + pEmt(1))**2
+      m2Final = m2Final - (pRad(2) + pRec(2) + pEmt(2))**2
+      m2Final = m2Final - (pRad(3) + pRec(3) + pEmt(3))**2
+
+C     Final state splitting not possible for negative dipole mass.
+      if ( radtype .eq. 1 .and. m2Final .lt. 0.0 ) allowed = .false.
+
+C     Rescale recoiler for final-intial splittings.
+      rescale = 1.0
+      if (radtype .eq. 1 .and. rectype .eq. -1) then
+        m2ar = m2Final - 2.0*Q2 + 2.0*m2RadBef
+        rescale = (1.0 - (Q2 - m2RadBef) / (m2ar-m2RadBef))
+     &           /(1.0 + (Q2 - m2RadBef) / (m2ar-m2RadBef))
+        pRecAft(0) = pRecAft(0)*rescale
+        pRecAft(1) = pRecAft(1)*rescale
+        pRecAft(2) = pRecAft(2)*rescale
+        pRecAft(3) = pRecAft(3)*rescale
+      endif
+
+C     Final-initial splitting not possible for negative rescaling.
+      if ( rescale .lt. 0.0 ) allowed = .false.
+
+C     Construct dipole momentum for FSR.
+      sum(0) = pRad(0) + pRecAft(0) + pEmt(0)
+      sum(1) = pRad(1) + pRecAft(1) + pEmt(1)
+      sum(2) = pRad(2) + pRecAft(2) + pEmt(2)
+      sum(3) = pRad(3) + pRecAft(3) + pEmt(3)
+      m2Dip = sum(0)**2 - sum(1)**2 - sum(2)**2 - sum(3)**2
+
+C     Construct 2->3 variables for FSR
+      x1     = 2. * ( sum(0)*pRad(0) - sum(1)*pRad(1)
+     &             -  sum(2)*pRad(2) - sum(3)*pRad(3) ) / m2Dip
+      x2     = 2. * ( sum(0)*pRecAft(0) - sum(1)*pRecAft(1)
+     &             -  sum(2)*pRecAft(2) - sum(3)*pRecAft(3) ) / m2Dip
+
+C     Final state splitting not possible for ill-defined
+C     3-body-variables.
+      if ( radtype .eq. 1 .and. x1 .lt. 0.0 ) allowed = .false.
+      if ( radtype .eq. 1 .and. x1 .gt. 1.0 ) allowed = .false.
+      if ( radtype .eq. 1 .and. x2 .lt. 0.0 ) allowed = .false.
+      if ( radtype .eq. 1 .and. x2 .gt. 1.0 ) allowed = .false.
+
+C     Auxiliary variables for massive FSR
+      lambda13 = DSQRT( (Q2 - m2RadAft - m2EmtAft )**2 - 4.0 * m2RadAft*m2EmtAft)
+      k1 = ( Q2 - lambda13 + (m2EmtAft - m2RadAft ) ) / ( 2.0 * Q2)
+      k3 = ( Q2 - lambda13 - (m2EmtAft - m2RadAft ) ) / ( 2.0 * Q2)
+
+C     Construct momenta of dipole before/after splitting for ISR
+      qBR(0) = pRad(0) + pRec(0) - pEmt(0)
+      qBR(1) = pRad(1) + pRec(1) - pEmt(1)
+      qBR(2) = pRad(2) + pRec(2) - pEmt(2)
+      qBR(3) = pRad(3) + pRec(3) - pEmt(3)
+      qBR2   = qBR(0)**2 - qBR(1)**2 - qBR(2)**2 - qBR(3)**2
+
+      qAR(0) = pRad(0) + pRec(0)
+      qAR(1) = pRad(1) + pRec(1)
+      qAR(2) = pRad(2) + pRec(2)
+      qAR(3) = pRad(3) + pRec(3)
+      qAR2   = qAR(0)**2 - qAR(1)**2 - qAR(2)**2 - qAR(3)**2
+
+C     Calculate z of splitting, different for FSR and ISR
+      z = 1.0 / (1.0 - k1 -k3) * ( x1 / (2.0-x2) - k3)
+      if(radtype .eq. -1 ) z = qBR2 / qAR2;
+
+C     Splitting not possible for ill-defined energy sharing.
+      if ( z .lt. 0.0 .or. z .gt. 1.0 ) allowed = .false.
+
+C     pT^2 = separation * virtuality (corrected with mass for FSR)
+      if (radtype .eq.  1) temp = z*(1-z)*(Q2 - m2RadBef)
+      if (radtype .eq. -1) temp = (1-z)*Q2
+
+C     Check threshold in ISR
+      if (radtype .ne. 1) then
+        if ((abs(flavRad) .eq. 4 .or. abs(flavEmt) .eq. 4)
+     &    .and. dsqrt(temp) .le. 2.0*m0c**2 ) temp = (1.-z)*(Q2+m0c**2)
+        if ((abs(flavRad) .eq. 5 .or. abs(flavEmt) .eq. 5)
+     &    .and. dsqrt(temp) .le. 2.0*m0b**2 ) temp = (1.-z)*(Q2+m0b**2)
+      endif
+
+C     Kinematically impossible splittings should not be included in the
+C     pT definition!
+      if( .not. allowed) temp = 1d15
+
+      if(temp .lt. 0.0) temp = 0.0
+
+C     Return pT
+      rhoPythia = dsqrt(temp);
+
+      RETURN
+      END

=== modified file 'Template/LO/SubProcesses/makefile'
--- Template/LO/SubProcesses/makefile	2016-05-26 23:36:59 +0000
+++ Template/LO/SubProcesses/makefile	2016-08-23 01:19:18 +0000
@@ -1,6 +1,13 @@
 include ../../Source/make_opts
 FFLAGS+= -w
 
+# Load additional dependencies of the bias module, if present
+ifeq (,$(wildcard ../bias_dependencies))
+BIASDEPENDENCIES =
+else
+include ../bias_dependencies
+endif
+
 # Definitions
 
 LIBDIR = ../../lib/
@@ -17,9 +24,9 @@
     MADLOOP_LIB =
 endif
 
-LINKLIBS = $(LINK_MADLOOP_LIB) $(LINK_LOOP_LIBS) -L../../lib/ -ldhelas -ldsample -lmodel -lgeneric -lpdf -lcernlib $(llhapdf) 
+LINKLIBS = $(LINK_MADLOOP_LIB) $(LINK_LOOP_LIBS) -L../../lib/ -ldhelas -ldsample -lmodel -lgeneric -lpdf -lcernlib $(llhapdf) -lbias 
 
-LIBS = $(LIBDIR)libdhelas.$(libext) $(LIBDIR)libdsample.$(libext) $(LIBDIR)libgeneric.$(libext) $(LIBDIR)libpdf.$(libext) $(LIBDIR)libmodel.$(libext) $(LIBDIR)libcernlib.$(libext) $(MADLOOP_LIB) $(LOOP_LIBS)
+LIBS = $(LIBDIR)libbias.$(libext) $(LIBDIR)libdhelas.$(libext) $(LIBDIR)libdsample.$(libext) $(LIBDIR)libgeneric.$(libext) $(LIBDIR)libpdf.$(libext) $(LIBDIR)libmodel.$(libext) $(LIBDIR)libcernlib.$(libext) $(MADLOOP_LIB) $(LOOP_LIBS)
 
 # Source files
 
@@ -34,7 +41,7 @@
 # Binaries
 
 $(PROG): $(PROCESS) auto_dsig.o $(LIBS)
-	$(FC) -o $(PROG) $(PROCESS) $(LINKLIBS) $(LDFLAGS)
+	$(FC) -o $(PROG) $(PROCESS) $(LINKLIBS) $(LDFLAGS) $(BIASDEPENDENCIES)
 
 gensym: $(SYMMETRY) configs.inc $(LIBDIR)libmodel.$(libext) $(LIBDIR)libgeneric.$(libext)
 	$(FC) -o gensym $(SYMMETRY) -L../../lib/ -lmodel -lgeneric $(LDFLAGS)

=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f	2016-06-07 09:33:55 +0000
+++ Template/LO/SubProcesses/reweight.f	2016-08-23 01:19:18 +0000
@@ -1087,7 +1087,53 @@
       endif
       return
       end
+
+      double precision function custom_bias(p, original_weight, numproc)
+c***********************************************************
+c     Returns a bias weight as instructed by the bias module
+c***********************************************************
+      implicit none
+
+      include 'nexternal.inc'
+      include 'maxparticles.inc'
+      include 'run_config.inc'
+      include 'lhe_event_infos.inc'
+      include 'run.inc'
+
+      DOUBLE PRECISION P(0:3,NEXTERNAL)
+      integer numproc
+
+      double precision original_weight
+
+      double precision bias_weight
+      logical is_bias_dummy, requires_full_event_info
+      common/bias/bias_weight,is_bias_dummy,requires_full_event_info
+
       
+C     If the bias module necessitates the full event information
+C     then we must call write_leshouches here already so as to set it.
+C     The weight specified at this stage is irrelevant since we
+C     use do_write_events set to .False.
+      AlreadySetInBiasModule = .False.      
+      if (requires_full_event_info) then
+        call write_leshouche(p,-1.0d0,numproc,.False.)
+C     Write the event in the string evt_record, part of the
+C     lhe_event_info common block
+        event_record(:) = ''
+        call write_event_to_stream(event_record,pb(0,1),1.0d0,npart,
+     &      jpart(1,1),ngroup,sscale,aaqcd,aaqed,buff,use_syst,
+     &      s_buff, nclus, buffclus)
+        AlreadySetInBiasModule = .True. 
+      else
+        AlreadySetInBiasModule = .False.
+      endif
+C     Apply the bias weight. The default run_card entry 'None' for the 
+c     'bias_weight' option will implement a constant bias_weight of 1.0 below.
+      call bias_wgt(p, original_weight, bias_weight)
+      custom_bias = bias_weight
+
+      end
+
       double precision function rewgt(p)
 c**************************************************
 c   reweight the hard me according to ckkw

=== modified file 'Template/LO/SubProcesses/setcuts.f'
--- Template/LO/SubProcesses/setcuts.f	2016-06-16 15:03:59 +0000
+++ Template/LO/SubProcesses/setcuts.f	2016-08-23 01:19:18 +0000
@@ -81,6 +81,11 @@
       logical  do_cuts(nexternal)
       COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,
      . IS_A_ONIUM, do_cuts 
+
+c     Store which external particles undergo the ktdurham and ptlund cuts.
+      LOGICAL  is_pdg_for_merging_cut(NEXTERNAL)
+      COMMON /TO_MERGE_CUTS/is_pdg_for_merging_cut, from_decay
+
 c
 c
 c     reading parameters
@@ -228,55 +233,16 @@
          if (abs(idup(i,1,iproc)).eq.16) is_a_nu(i)=.true.  ! no cuts on vt vt~
          if (pmass(i).gt.10d0)     is_heavy(i)=.true. ! heavy fs particle
 c-onium
-c         if (idup(i,1,iproc).eq.441)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.10441)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.100441)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.443)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.10443)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.20443)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.100443)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.30443)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.9000443)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.9010443)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.9020443)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.445)   is_a_onium(i)=.true. ! charmonium
-c         if (idup(i,1,iproc).eq.9000445)   is_a_onium(i)=.true. ! charmonium
-
-c         if (idup(i,1,iproc).eq.551)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.10551)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.100551)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.110551)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.200551)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.210551)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.10553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.20553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.30553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.100553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.110553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.120553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.130553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.200553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.210553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.220553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.300553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.9000553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.9010553)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.555)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.10555)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.20555)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.100555)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.110555)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.200555)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.557)   is_a_onium(i)=.true. ! bottomonium
-c         if (idup(i,1,iproc).eq.100557)   is_a_onium(i)=.true. ! bottomonium
-
-c         if (idup(i,1,iproc).eq.541)   is_a_onium(i)=.true. ! Bc
-c         if (idup(i,1,iproc).eq.10541)   is_a_onium(i)=.true. ! Bc
-c         if (idup(i,1,iproc).eq.543)   is_a_onium(i)=.true. ! Bc
-c         if (idup(i,1,iproc).eq.10543)   is_a_onium(i)=.true. ! Bc
-c         if (idup(i,1,iproc).eq.20543)   is_a_onium(i)=.true. ! Bc
-c         if (idup(i,1,iproc).eq.545)   is_a_onium(i)=.true. ! Bc
+
+c        Remember mergeable particles
+         do j=1,pdgs_for_merging_cut(0)
+           if (    pdgs_for_merging_cut(j) .ne. 0
+     $       .and. abs(idup(i,1,iproc)) .eq.pdgs_for_merging_cut(j) ) then
+             is_pdg_for_merging_cut(i)=.true.
+             exit
+            endif
+         enddo
+
       enddo
 
 

=== modified file 'Template/LO/SubProcesses/unwgt.f'
--- Template/LO/SubProcesses/unwgt.f	2016-05-05 11:12:21 +0000
+++ Template/LO/SubProcesses/unwgt.f	2016-08-23 01:19:18 +0000
@@ -236,9 +236,9 @@
 c            call write_event(p,uwgt)
 c            write(29,'(2e15.5)') matrix,wgt
 c $B$ S-COMMENT_C $B$
-            call write_leshouche(p,uwgt,numproc)
+            call write_leshouche(p,uwgt,numproc,.True.)
          elseif (xwgt .gt. 0d0 .and. nw .lt. 5) then
-            call write_leshouche(p,wgt/twgt*1d-6,numproc)
+            call write_leshouche(p,wgt/twgt*1d-6,numproc,.True.)
 c $E$ S-COMMENT_C $E$
          endif
          maxwgt=max(maxwgt,xwgt)
@@ -437,7 +437,7 @@
 c      close(lun)
       end
 
-      SUBROUTINE write_leshouche(p,wgt,numproc)
+      SUBROUTINE write_leshouche(p,wgt,numproc,do_write_events)
 C**************************************************************************
 C     Writes out information for event
 C**************************************************************************
@@ -460,15 +460,16 @@
 c
       double precision p(0:3,nexternal),wgt
       integer numproc
+      logical do_write_events
 c
 c     Local
 c
       integer i,j,k,iini,ifin
       double precision sum_wgt,sum_wgt2, xtarget,targetamp(maxflow)
-      integer ip, np, ic, nc, jpart(7,-nexternal+3:2*nexternal-3)
+      integer ip, np, ic, nc
       integer ida(2),ito(-nexternal+3:nexternal),ns,nres,ires,icloop
       integer iseed
-      double precision pboost(0:3),pb(0:4,-nexternal+3:2*nexternal-3)
+      double precision pboost(0:3)
       double precision beta, get_betaz
       double precision ebi(0:3), ebo(0:3)
       double precision ptcltmp(nexternal), pdum(0:3)
@@ -477,19 +478,14 @@
       integer mothup(2,nexternal)
       integer icolup(2,nexternal,maxflow,maxsproc)
 
-      integer isym(nexternal,99), nsym, jsym
+      integer nsym
 
-      double precision sscale,aaqcd,aaqed
-      integer ievent,npart
+      integer ievent
       logical flip
 
       real ran1
       external ran1
 
-      character*1000 buff
-      character*(s_bufflen) s_buff(7)
-      integer nclus
-      character*(clus_bufflen) buffclus(nexternal)
       character*40 cfmt
 C     
 C     GLOBAL
@@ -510,9 +506,6 @@
       integer           mincfig, maxcfig
       common/to_configs/mincfig, maxcfig
 
-      integer ngroup
-      common/to_group/ngroup
-
       double precision stot,m1,m2
       common/to_stot/stot,m1,m2
 c
@@ -529,17 +522,24 @@
 c      data ncolflow/maxamps*0/
 c      data ncolalt/maxamps*0/
 
+      include 'coupl.inc'
+
+      include 'lhe_event_infos.inc'
+      data AlreadySetInBiasModule/.False./
+
       include 'symswap.inc'
-      include 'coupl.inc'
 C-----
 C  BEGIN CODE
 C-----
       
-      if (nw .ge. maxevents) return
-
-c     Store weight for event
-      nw = nw+1
-      swgt(nw)=wgt
+      if ((nw .ge. maxevents).and.do_write_events) return
+
+C     if all the necessary inputs to write the events have already been
+C     computed in the bias module, then directly jump to write_events
+      if (AlreadySetInBiasModule) then
+        goto 1123
+      endif
+
 c
 c     In case of identical particles symmetry, choose assignment
 c
@@ -746,6 +746,19 @@
          write(buffclus(nexternal),'(a)')'</clustering>'
       endif
 
+C     If the arguments of write_event have already been set in the
+C     bias module, then the beginning of the routine will directly
+C     jump here.
+
+ 1123 continue
+      if (.not.do_write_events) then
+        return
+      endif
+
+c     Store weight for event
+      nw = nw+1
+      swgt(nw)=wgt
+
       call write_event(lun,pb(0,1),wgt,npart,jpart(1,1),ngroup,
      &   sscale,aaqcd,aaqed,buff,use_syst,s_buff,nclus,buffclus)
       if(btest(mlevel,1))

=== modified file 'Template/LO/bin/internal/run_delphes3'
--- Template/LO/bin/internal/run_delphes3	2013-02-25 00:52:53 +0000
+++ Template/LO/bin/internal/run_delphes3	2016-08-23 01:19:18 +0000
@@ -9,6 +9,7 @@
 run=$2
 tag=$3
 cross=$4 
+file=$5
 
 main=`pwd`
 
@@ -22,17 +23,34 @@
   exit
 fi
 
+if [ ${file: -3} == ".gz" ]
+then
+  if [ ${file: -7} == ".hep.gz" ]
+  then
+      gunzip --stdout $file | $delphesdir/DelphesSTDHEP ../Cards/delphes_card.dat  ${run}/${tag}_delphes_events.root      
+  else
+      gunzip --stdout $file | $delphesdir/DelphesHepMC ../Cards/delphes_card.dat  ${run}/${tag}_delphes_events.root      
+  fi
+else
+  if [ ${file: -4} == ".hep" ]
+  then
+      $delphesdir/DelphesSTDHEP ../Cards/delphes_card.dat  ${run}/${tag}_delphes_events.root $file
+  else
+      $delphesdir/DelphesHepMC ../Cards/delphes_card.dat  ${run}/${tag}_delphes_events.root $file      
+  fi
+fi
+
 echo $$ >> ../myprocid
 
-# Set delphes path in delphes_card.dat
-
-$delphesdir/DelphesSTDHEP ../Cards/delphes_card.dat delphes.root pythia_events.hep
-$delphesdir/root2lhco delphes.root delphes_events.lhco
-
-if [ -e delphes_events.lhco ]; then
-# write the delphes banner
-    sed -e "s/^/#/g" ${run}/${run}_${tag}_banner.txt > ${run}/${tag}_delphes_events.lhco
-    echo "##  Integrated weight (pb)  : ${cross}" >> ${run}/${tag}_delphes_events.lhco
-    cat delphes_events.lhco >> ${run}/${tag}_delphes_events.lhco
-    rm -f delphes_events.lhco
-fi
+# Uncomment the following to have the LHCO file:
+#
+#
+#$delphesdir/root2lhco ${run}/${tag}_delphes.root delphes_events.lhco
+#
+#if [ -e delphes_events.lhco ]; then
+## write the delphes banner
+#    sed -e "s/^/#/g" ${run}/${run}_${tag}_banner.txt > ${run}/${tag}_delphes_events.lhco
+#    echo "##  Integrated weight (pb)  : ${cross}" >> ${run}/${tag}_delphes_events.lhco
+#    cat delphes_events.lhco >> ${run}/${tag}_delphes_events.lhco
+#    rm -f delphes_events.lhco
+#fi

=== modified file 'Template/NLO/SubProcesses/BinothLHA.f'
--- Template/NLO/SubProcesses/BinothLHA.f	2016-03-02 09:31:23 +0000
+++ Template/NLO/SubProcesses/BinothLHA.f	2016-08-23 01:19:18 +0000
@@ -85,6 +85,8 @@
 c itself again to perform stability check to make sure no unstable EPS
 c splips unnoticed.
          CALL FORCE_STABILITY_CHECK(.TRUE.)
+         CALL COLLIER_COMPUTE_UV_POLES(.FALSE.)
+         CALL COLLIER_COMPUTE_IR_POLES(.FALSE.)
          firsttime_run = .false.
       endif
       if (firsttime) then
@@ -102,6 +104,7 @@
 c Just set the accuracy found to a positive value as it is not specified
 c once the initial pole check is performed.
          if (mc_hel.eq.0) then
+
             call sloopmatrix_thres(p,virt_wgts,tolerance,accuracies
      $           ,ret_code)
             prec_found = accuracies(0)            
@@ -164,18 +167,19 @@
 c MadLoop initialization PS points.
       cpol=.false.
       if ((firsttime .or. mc_hel.eq.0) .and. mod(ret_code,100)/10.ne.3
-     $     .and. mod(ret_code,100)/10.ne.4) then 
+     $     .and. mod(ret_code,100)/10.ne.4) then
          call getpoles(p,QES2,madfks_double,madfks_single,fksprefact)
          avgPoleRes(1)=(single+madfks_single)/2.0d0
          avgPoleRes(2)=(double+madfks_double)/2.0d0
          PoleDiff(1)=dabs(single - madfks_single)
          PoleDiff(2)=dabs(double - madfks_double)
          if ((dabs(avgPoleRes(1))+dabs(avgPoleRes(2))).ne.0d0) then
-            cpol = .not. (((PoleDiff(1)+PoleDiff(2))/
+            cpol = .not. ((((PoleDiff(1)+PoleDiff(2))/
      $           (dabs(avgPoleRes(1))+dabs(avgPoleRes(2)))) .lt.
-     $           tolerance*10d0)
+     $           tolerance*10d0).or.(mod(ret_code,10).eq.7))
          else
-            cpol = .not.(PoleDiff(1)+PoleDiff(2).lt.tolerance*10d0)
+            cpol = .not.((PoleDiff(1)+PoleDiff(2).lt.tolerance*10d0).or.
+     $                   (mod(ret_code,10).eq.7))
          endif
          if (tolerance.lt.0.0d0) then
             cpol = .false.

=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f	2016-03-02 09:31:23 +0000
+++ Template/NLO/SubProcesses/check_poles.f	2016-08-23 01:19:18 +0000
@@ -116,6 +116,8 @@
 c Make sure that stability checks are always used by MadLoop, even for
 c initialization
       CALL FORCE_STABILITY_CHECK(.TRUE.)
+      CALL COLLIER_COMPUTE_UV_POLES(.TRUE.)
+      CALL COLLIER_COMPUTE_IR_POLES(.TRUE.)
 
 200   continue
           calculatedborn = .false.

=== modified file 'Template/NLO/SubProcesses/makefile_loop.inc'
--- Template/NLO/SubProcesses/makefile_loop.inc	2016-03-02 03:03:20 +0000
+++ Template/NLO/SubProcesses/makefile_loop.inc	2016-08-23 01:19:18 +0000
@@ -24,6 +24,7 @@
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard mp_coef_construction_*.f)) \
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard TIR_interface.f)) \
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard GOLEM_interface.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard COLLIER_interface.f)) \
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard compute_color_flows.f))
 
 # This is the core of madloop computationally wise, so make sure to turn optimizations on and bound checks off.

=== modified file 'Template/loop_material/StandAlone/Cards/MadLoopParams.dat'
--- Template/loop_material/StandAlone/Cards/MadLoopParams.dat	2016-04-06 10:19:13 +0000
+++ Template/loop_material/StandAlone/Cards/MadLoopParams.dat	2016-08-23 01:19:18 +0000
@@ -4,8 +4,8 @@
 
 !     
 #MLReductionLib
-!6|1|2
-! Default :: 6|1|2
+!6|7|1
+! Default :: 6|7|1
 !     The tensor integral reduction library.The current choices are:
 !     1 | CutTools
 !     2 | PJFry++
@@ -13,6 +13,7 @@
 !     4 | Golem95
 !     5 | Samurai
 !     6 | Ninja
+!     7 | COLLIER
 !     One can use the combinations to reduce integral,e.g.
 !     1|2|3 means first use CutTools, if it is not stable, use PJFry++, 
 !     if it is still unstable, use IREGI. If it failed, use QP of CutTools.
@@ -22,15 +23,15 @@
 !     When using quadruple precision with Ninja or CutTools, the reduction will 
 !     always be done in quadruple precision, but the parameters below allow you to
 !     chose if you want to also recompute the *integrand* in quadruple precision.
-!     Doing so is slow and typically does not improve the accuracy for Ninja, so it 
-!     is disabled by default in this case.
+!     Doing so is slow but might improve the accuracy in some situation.
 #UseQPIntegrandForCutTools
 !.TRUE.
 ! Default :: .TRUE.
 #UseQPIntegrandForNinja
-!.FALSE.
-! Default :: .FALSE.
+!.TRUE.
+! Default :: .TRUE.
 !
+
 !     =================================================================================
 !     The parameters below set the parameters for IREGI
 !     =================================================================================
@@ -38,9 +39,9 @@
 #IREGIMODE
 !2
 ! Default :: 2
-!     IMODE=0, IBP reduction
-!     IMODE=1, PaVe reduction
-!     IMODE=2, PaVe reduction with stablility improved by IBP reduction
+!     IREGIMODE=0, IBP reduction
+!     IREGIMODE=1, PaVe reduction
+!     IREGIMODE=2, PaVe reduction with stablility improved by IBP reduction
 
 #IREGIRECY
 !.TRUE.
@@ -59,7 +60,7 @@
 !       2   |  Double precision, loops reduced with propagator with reversed order                               
 !       4   |  Quadruple precision, loops reduced with propagator in original order                               
 !       5   |  Quadruple precision, loops reduced with propagator with reversed order                          
-!      -1   |  Exhaustive automated numerical stability checks. See below for details.
+!      -1   | Exhaustive automated numerical stability checks. See below for details.
 !     
 !     Due to the architecture of the program, you are better off
 !     rerunning the full PS point in quadruple precision than just a single loop 
@@ -84,8 +85,8 @@
 !     Below you can chose if you want to use zero, one or two rotations of the PS point 
 !     in QP.
 #NRotations_DP
-!1
-! Default :: 1
+!0
+! Default :: 0
 #NRotations_QP
 !0
 ! Default :: 0
@@ -114,6 +115,7 @@
 #ImprovePSPoint
 !2
 ! Default :: 2
+
 !     =================================================================================
 !     The parameters below set two CutTools internal parameters accessible to the user. 
 !     =================================================================================
@@ -226,4 +228,71 @@
 !.FALSE.
 ! Default :: .FALSE.
 
+!     =================================================================================
+!     The parameters below set the main parameters for COLLIER
+!     To edit more specific technical COLLIER parameters, modify directly the content
+!     of the subroutine 'INITCOLLIER' in the file 'MadLoopCommons.f'
+!     =================================================================================
+
+!     Decide if COLLIER must be computed multiple times to evaluate the UV pole residues
+!     (Withing a Monte-Carlo performed in MG5aMC, this is automatically disabled internally)
+#COLLIERComputeUVpoles
+!.TRUE.
+! Default :: .TRUE.
+
+!     Decide if COLLIER must be computed multiple times to evaluate the IR pole residues
+!     (Withing a Monte-Carlo performed in MG5aMC, this is automatically disabled internally)
+#COLLIERComputeIRpoles
+!.TRUE.
+! Default :: .TRUE.
+
+!     Decide if COLLIER must be computed multiple times to evaluate the IR pole residues
+#COLLIERRequiredAccuracy
+!1.0d-8
+! Default :: 1.0d-8
+!     A value of -1.0d0 means that it will be automatically set from MLStabThres.
+!     The default value of 1.0d-8 corresponds to the value for which COLLIER's authors
+!     have optimized the library.
+
+!     Decide whether to use COLLIER's internal stability test or the loop-direction
+!     switch test instead.
+#COLLIERUseInternalStabilityTest
+!.TRUE.
+! Default :: .TRUE.
+!     COLLIER's internal stability test is at no extra cost but not as reliable
+!     as the loop-direction switch test, which however doubles the reduction time.
+!     This parameter is only relevant when running MadLoop with CTModeRun=-1.
+!     If you find a large number of unstable points with COLLIER for complicated
+!     processes, set this parameter to .FALSE. to make sure the PS points flagged
+!     as unstable with COLLIER really are so.
+
+!     Set up to which N-loop to use the COLLIER global caching system.
+#COLLIERGlobalCache
+!-1
+! Default :: -1
+!     -1 : Enable the global cache for all loops 
+!      0 : Disable the global cache alltogether
+!      N : Enable the global cache but only for up to N-loops
+
+!     Use the global cache when evaluating the poles as well (more memory consuming)
+!     During a Monte-Carlo it is typically not useful anyway, because the pole  
+!     computation is automatically disabled for COLLIER, irrespectively of the value
+!     of the parameters COLLIERCompute<UV/IR>poles specified above.
+#COLLIERUseCacheForPoles
+!.FALSE.
+! Default :: .FALSE.
+
+!     Choose which branch(es) of COLLIER have to be used
+#COLLIERMode
+!1
+! Default :: 1
+!     COLLIERMode=1 : COLI branch 
+!     COLLIERMode=2 : DD branch
+!     COLLIERMode=3 : Both DD and COLI branch compared
+
+!     Decide if COLLIER can output its information in a log directory.
+#COLLIERCanOutput
+!.FALSE.
+! Default :: .FALSE.
+
 /* End of param file */

=== modified file 'Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.inc'
--- Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.inc	2016-02-08 23:58:00 +0000
+++ Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.inc	2016-08-23 01:19:18 +0000
@@ -30,6 +30,86 @@
 
       end
 
+
+## if(collier_available){
+	  SUBROUTINE INITCOLLIER()
+C 
+C INITIALISATION OF COLLIER
+C
+C  
+C MODULE
+C
+	  USE COLLIER 
+C
+C CONSTANTS
+C
+	  character(len=*) no_folder 
+      parameter (no_folder='')
+	  character(len=*) folderOutput
+      parameter (folderOutput='COLLIER_output')
+
+C     Force COLLIER to completely reset from scratch
+      logical noreset
+      parameter (noreset=.FALSE.)
+C
+C LOCAL VARIABLES 
+C    
+      INTEGER N_CACHES
+C  
+C GLOBAL VARIABLES 
+C
+      include 'MadLoopParams.inc'
+C     Now obtain the overall maximum rank, maximum external lines
+C     and maximal number of processes.
+      include 'global_specs.inc'
+
+C ----------
+C BEGIN CODE
+C ----------
+
+C     Initialize Collier
+      IF (.NOT.COLLIERCanOutput) THEN
+        CALL Init_cll(MAXNEXTERNAL,OVERALLMAXRANK,no_folder,noreset)
+      ELSE
+        CALL Init_cll(MAXNEXTERNAL,OVERALLMAXRANK,folderOutput,noreset)
+	  ENDIF
+
+C     Set target accuracy, be conservative w.r.t user request
+C     Keep in mind that COLLIER has been optimized for 1d-8.
+      IF (COLLIERRequiredAccuracy.eq.-1.0d0) THEN
+        CALL SetReqAcc_cll(MLStabThres*1.0d-3)
+      ELSE
+        CALL SetReqAcc_cll(COLLIERRequiredAccuracy)
+	  ENDIF
+
+C     Set COLLIER mode
+      CALL SetMode_cll(COLLIERMode)
+
+C     Set the global cache strategy
+C     If we use the caches for the poles as well, then it is a total of
+C     4 caches per process.
+      IF (COLLIERUseCacheForPoles) THEN
+	    N_CACHES =4*NPROCS
+	  ELSE
+	    N_CACHES =NPROCS
+	  ENDIF
+      IF (COLLIERGlobalCache.eq.-1) THEN
+	    CALL InitCacheSystem_cll(N_CACHES,MAXNEXTERNAL)
+	  ELSEIF(COLLIERGlobalCache.gt.0) THEN
+	    CALL InitCacheSystem_cll(N_CACHES,COLLIERGlobalCache)
+	  ENDIF
+
+C     Make sure to start by first switching off all cache
+      IF (COLLIERGlobalCache.ne.0) THEN
+        call SwitchOffCacheSystem_cll()
+      ENDIF
+
+C     Specify below your other custom COLLIER parameter settings 
+C     [user_specific_COLLIER_settings]
+	  
+	  END
+## }
+
       SUBROUTINE SET_FORBID_HEL_DOUBLECHECK(ONOFF)
 C
 C Give the possibility to overwrite the value of MadLoopParams.dat
@@ -67,9 +147,10 @@
       DATA ML_INIT/.TRUE./
       common/ML_INIT/ML_INIT
 
-      LOGICAL CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT
-      DATA CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT/.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE./
-      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
+      LOGICAL CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT
+      DATA CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT/.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE./
+      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT, COLLIERINIT
+
 
       character(512) MLPath
       data MLPath/'[[NA]]'/      
@@ -163,6 +244,8 @@
 C         Stable with Samurai.
 C     U == 6
 C         Stable with Ninja in double precision.
+C     U == 7
+C         Stable with COLLIER.
 C     U == 8
 C         Stable with Ninja in quadruple precision.
 C     U == 9
@@ -201,10 +284,10 @@
             STOP 'Only CutTools and Ninja can use quardruple precision'
          ENDIF
       ENDIF
-      IF(MLRed.GE.1.AND.MLRed.LE.6)THEN
+      IF(MLRed.GE.1.AND.MLRed.LE.7)THEN
          SET_RET_CODE_U=MLRed
       ELSE
-         STOP 'Only CutTools, PJFry++, IREGI, Golem95, Samurai and Ninja are available'
+         STOP 'Only CutTools, PJFry++, IREGI, Golem95, Samurai, Ninja and COLLIER are available'
       ENDIF
       END
 
@@ -249,8 +332,11 @@
       ELSEIF(LIBNUM.EQ.6)THEN
 C Ninja 
          CALL DETECT_NINJA(NLOOPLINE,RANK,complex_mass,LPASS)
+      ELSEIF(LIBNUM.EQ.7)THEN
+C Collier 
+         CALL DETECT_COLLIER(NLOOPLINE,RANK,complex_mass,LPASS)
       ELSE
-         STOP 'ONLY CUTTOOLS,PJFry++,IREGI,Golem95 and Samurai are available'
+         STOP 'Only CutTools, PJFry++, IREGI, Golem95, Samurai, Ninja and COLLIER are available'
       ENDIF
       RETURN
       END
@@ -358,6 +444,65 @@
       RETURN
       END
 
+## if(collier_available){
+      SUBROUTINE DETECT_COLLIER(NLOOPLINE,RANK,complex_mass,LPASS)
+C
+C Detect whether Collier can be used or not
+C
+      USE COLLIER
+      IMPLICIT NONE
+C
+C CONSTANTS
+C
+C
+C ARGUMENTS
+C
+      INTEGER NLOOPLINE,RANK
+      LOGICAL complex_mass,LPASS
+C
+C LOCAL VARIABLES
+C
+      INTEGER CURRENT_COLLIERMODE
+C
+C GLOBAL VARIABLES
+C
+      LOGICAL CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT
+      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT, COLLIERINIT
+      include 'MadLoopParams.inc'
+C ----------
+C BEGIN CODE
+C ----------
+      if (.NOT.COLLIERINIT) THEN
+        CALL GETMODE_CLL(CURRENT_COLLIERMODE)
+      ELSE
+        CURRENT_COLLIERMODE = COLLIERMODE
+      ENDIF
+      LPASS=.TRUE.
+	  IF (CURRENT_COLLIERMODE.NE.1) THEN
+C       The DD branch is used and it has limitations
+        IF((NLOOPLINE.GT.6).OR.(RANK.GT.NLOOPLINE)) then
+	      LPASS=.FALSE.
+        endif
+      ELSE
+C       Limitations of the COLI branch are academic.
+        LPASS=.TRUE. 
+	  ENDIF
+      RETURN
+      END
+## }else{
+      SUBROUTINE DETECT_COLLIER(NLOOPLINE,RANK,complex_mass,LPASS)
+C
+C ARGUMENTS
+C
+      INTEGER NLOOPLINE,RANK
+      LOGICAL complex_mass,LPASS
+C
+C     COLLIER is not available in this output. This subroutine is dummy.
+C
+      LPASS=.True.
+      END
+## }
+
       SUBROUTINE DETECT_PJFRY(NLOOPLINE,RANK,complex_mass,LPASS)
 C
 C DETECT whether PJFRY++ CAN BE USED OR NOT

=== modified file 'Template/loop_material/StandAlone/SubProcesses/MadLoopParamReader.f'
--- Template/loop_material/StandAlone/SubProcesses/MadLoopParamReader.f	2016-02-18 11:42:56 +0000
+++ Template/loop_material/StandAlone/SubProcesses/MadLoopParamReader.f	2016-08-23 01:19:18 +0000
@@ -33,6 +33,12 @@
                stop 'CTModeRun must be >= -1 and <=6.'
              endif 
 
+           else if (buff .eq. '#COLLIERGlobalCache') then
+             read(666,*,end=999) COLLIERGlobalCache
+             if (COLLIERGlobalCache .lt. -1) then
+               stop 'COLLIERGlobalCache must be >= -1'
+             endif
+
            else if (buff .eq. '#NRotations_DP') then
              read(666,*,end=999) NRotations_DP
              if (NRotations_DP .lt. 0 .or.
@@ -52,7 +58,14 @@
              if (MLStabThres.lt.0.0d0) then
                  stop 'MLStabThres must be >= 0'
              endif
-          
+
+           else if (buff .eq. '#COLLIERRequiredAccuracy') then
+             read(666,*,end=999) COLLIERRequiredAccuracy
+             if (COLLIERRequiredAccuracy.le.0.0d0.and.
+     &           COLLIERRequiredAccuracy.ne.-1.0d0) then
+                 stop 'COLLIERRequiredAccuracy must be > 0 or = -1.0'
+             endif
+
            else if (buff .eq. '#CTLoopLibrary') then
              read(666,*,end=999) CTLoopLibrary
              if (CTLoopLibrary.lt.2 .or.
@@ -90,6 +103,21 @@
                  stop 'MaxAttempts must be >= 1'
              endif
 
+          else if (buff .eq. '#COLLIERComputeUVpoles') then
+             read(666,*,end=999) COLLIERComputeUVpoles
+
+          else if (buff .eq. '#COLLIERComputeIRpoles') then
+             read(666,*,end=999) COLLIERComputeIRpoles
+
+		  else if (buff .eq. '#COLLIERUseInternalStabilityTest') then
+             read(666,*,end=999) COLLIERUseInternalStabilityTest
+