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

[Olivier Mattelaer]

Hi guys,

As announced, I freeze now this version.
Please do not apply any change to this version now.
All modification should be done in 2.5.1 (not yet created).
In case of critical bug, please contact me first.

I will put this version available as beta, as soon as I have run all the parallel tests
(as usual count at least a week).

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/2.5.0]

297. By Olivier Mattelaer on 2016-09-03

small fix/new tests

298. By Olivier Mattelaer on 2016-09-03

fix an mg5py8interface installation problem + fix tests

299. By Olivier Mattelaer on 2016-09-05

fix html for py8 with merging + fix decay_comparator parralel tests

300. By Olivier Mattelaer on 2016-09-05

add default card for ma5 to avoid crash if those card are not created by ma5

301. By Olivier Mattelaer on 2016-09-05

fix a series of small bug

302. By Olivier Mattelaer on 2016-09-06

change mssm model to slha2 convention all the way

303. By Olivier Mattelaer on 2016-09-06

update tests
 - fix some linked to mssm
 - pass some acceptance test to parralel tests
 - improve speed of some tests
include some of Val fix of 2.5.1

304. By Olivier Mattelaer on 2016-09-06

fixing short parralel test related to the mssm modification

305. By Olivier Mattelaer on 2016-09-07

fixing the automatic installation of MA4 for the acceptance tests + make mg5 to handle correctly crash of MA5 at the time of creation of the cards

306. By Olivier Mattelaer on 2016-09-07

1) allow to filter function for the IOTest
2) remove the ML print_banner fct for the IOTest
3) fix a bug in the tests/test_manager related to IOTest
4) secure the interface to MA5

307. By Olivier Mattelaer on 2016-09-07

fix create_release

308. By Valentin Hirschi on 2016-09-07

1. Added the shortcut 'simplepy8' to turnoff advanced slow features of PY8.

309. By Valentin Hirschi on 2016-09-07

1. Added a comment in PY8 card concerning the possibility of turning off mpi.

310. By Olivier Mattelaer on 2016-09-07

edit interface modification
- add a new shortcut mpi
- improve shortcut treatment (possibility to define type of argument/number of arguement)
- add suggestions of parameter if the set command is unknow

311. By Valentin Hirschi on 2016-09-07

1. Fixed typo in mpi comment in PY8 card.

312. By Valentin Hirschi on 2016-09-07

1. modified compeltion for PY8 so as to show all parameters, not only visible ones.

313. By Olivier Mattelaer on 2016-09-07

remove a useless import triggering depreciationWarning

314. By Olivier Mattelaer on 2016-09-07

secure combine_events and systematics if problem with the zipping. + fixing combine events in presence of partials unweight

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2016-07-06 08:54:04 +0000
+++ MadSpin/decay.py	2016-09-07 23:55:24 +0000
@@ -64,7 +64,7 @@
 from madgraph import MG5DIR, MadGraph5Error
 import madgraph.various.misc as misc
 #import time
-import tests.parallel_tests.test_aloha as test_aloha
+
 
 class MadSpinError(MadGraph5Error):
     pass
@@ -2615,7 +2615,7 @@
     
 
     @misc.mute_logger()
-    @test_aloha.set_global()
+    @misc.set_global()
     def generate_all_matrix_element(self):
         """generate the full series of matrix element needed by Madspin.
         i.e. the undecayed and the decay one. And associate those to the 

=== added directory 'PLUGIN'
=== added file 'PLUGIN/__init__.py'
=== modified file 'Template/Common/Cards/delphes_card_ATLAS.dat'
--- Template/Common/Cards/delphes_card_ATLAS.dat	2015-11-10 11:23:55 +0000
+++ Template/Common/Cards/delphes_card_ATLAS.dat	2016-09-07 23:55:24 +0000
@@ -31,6 +31,8 @@
 
   NeutrinoFilter
   GenJetFinder
+  GenMissingET
+  
   FastJetFinder
 
   JetEnergyScale
@@ -138,9 +140,9 @@
   # set ResolutionFormula {resolution formula as a function of eta and pt}
 
   # resolution formula for charged hadrons
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
 }
 
 ###################################
@@ -154,9 +156,9 @@
   # set ResolutionFormula {resolution formula as a function of eta and energy}
 
   # resolution formula for electrons
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
 }
 
 ###############################
@@ -168,11 +170,10 @@
   set OutputArray muons
 
   # set ResolutionFormula {resolution formula as a function of eta and pt}
-
   # resolution formula for muons
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.02^2 + pt^2*2.0e-4^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*3.0e-4^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*6.0e-4^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
 }
 
 ##############
@@ -312,7 +313,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
 }
 
 #################
@@ -358,7 +359,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
 }
 
 #################
@@ -392,7 +393,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.25
 }
 
 ###################
@@ -456,6 +457,18 @@
 }
 
 
+#########################
+# Gen Missing ET merger
+########################
+
+module Merger GenMissingET {
+# add InputArray InputArray
+  add InputArray NeutrinoFilter/filteredParticles
+  set MomentumOutputArray momentum
+}
+
+
+
 ############
 # Jet finder
 ############
@@ -530,23 +543,36 @@
 # tau-tagging
 #############
 
-module TauTagging TauTagging {
+module TrackCountingTauTagging TauTagging {
+ 
   set ParticleInputArray Delphes/allParticles
   set PartonInputArray Delphes/partons
+  set TrackInputArray TrackMerger/tracks
   set JetInputArray JetEnergyScale/jets
 
-  set DeltaR 0.5
+  set DeltaR 0.2
+  set DeltaRTrack 0.2
 
+  set TrackPTMin 1.0
+ 
   set TauPTMin 1.0
-
   set TauEtaMax 2.5
 
-  # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt}
+  # instructions: {n-prongs} {eff} 
+  
+  # 1 - one prong efficiency
+  # 2 - two or more efficiency
+  # -1 - one prong mistag rate
+  # -2 - two or more mistag rate
+ 
+  set BitNumber 0
+ 
+  # taken from ATL-PHYS-PUB-2015-045 (medium working point)
+  add EfficiencyFormula {1} {0.70}
+  add EfficiencyFormula {2} {0.60}
+  add EfficiencyFormula {-1} {0.02}
+  add EfficiencyFormula {-2} {0.01}
 
-  # default efficiency formula (misidentification rate)
-  add EfficiencyFormula {0} {0.01}
-  # efficiency formula for tau-jets
-  add EfficiencyFormula {15} {0.6}
 }
 
 #####################################################
@@ -582,6 +608,8 @@
   add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
 
   add Branch GenJetFinder/jets GenJet Jet
+  add Branch GenMissingET/momentum GenMissingET MissingET
+
   add Branch UniqueObjectFinder/jets Jet Jet
   add Branch UniqueObjectFinder/electrons Electron Electron
   add Branch UniqueObjectFinder/photons Photon Photon

=== modified file 'Template/Common/Cards/delphes_card_CMS.dat'
--- Template/Common/Cards/delphes_card_CMS.dat	2015-11-10 11:23:55 +0000
+++ Template/Common/Cards/delphes_card_CMS.dat	2016-09-07 23:55:24 +0000
@@ -31,6 +31,8 @@
 
   NeutrinoFilter
   GenJetFinder
+  GenMissingET
+  
   FastJetFinder
 
   JetEnergyScale
@@ -139,9 +141,9 @@
 
   # resolution formula for charged hadrons
   # based on arXiv:1405.6569
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
 }
 
 ###################################
@@ -156,9 +158,9 @@
 
   # resolution formula for electrons
   # based on arXiv:1405.6569
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
 }
 
 ###############################
@@ -172,9 +174,9 @@
   # set ResolutionFormula {resolution formula as a function of eta and pt}
 
   # resolution formula for muons
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*2.0e-4^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.02^2 + pt^2*3.0e-4^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*6.0e-4^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
 }
 
 ##############
@@ -266,8 +268,10 @@
   add EnergyFraction {3122} {0.3 0.7}
 
   # set ECalResolutionFormula {resolution formula as a function of eta and energy}
-  set ECalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.007^2 + energy*0.07^2 + 0.35^2)  +
-                             (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
+  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
+  set ECalResolutionFormula {                  (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
 
   # set HCalResolutionFormula {resolution formula as a function of eta and energy}
   set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
@@ -317,7 +321,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
 }
 
 #################
@@ -363,7 +367,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
 }
 
 #################
@@ -399,7 +403,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.25
 }
 
 ###################
@@ -463,6 +467,17 @@
   set JetPTMin 20.0
 }
 
+#########################
+# Gen Missing ET merger
+########################
+
+module Merger GenMissingET {
+# add InputArray InputArray
+  add InputArray NeutrinoFilter/filteredParticles
+  set MomentumOutputArray momentum
+}
+
+
 
 ############
 # Jet finder
@@ -526,7 +541,7 @@
   # based on arXiv:1211.4462
   
   # default efficiency formula (misidentification rate)
-  add EfficiencyFormula {0} {0.01+0.00038*pt}
+  add EfficiencyFormula {0} {0.01+0.000038*pt}
 
   # efficiency formula for c-jets (misidentification rate)
   add EfficiencyFormula {4} {0.25*tanh(0.018*pt)*(1/(1+ 0.0013*pt))}
@@ -591,6 +606,8 @@
   add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
 
   add Branch GenJetFinder/jets GenJet Jet
+  add Branch GenMissingET/momentum GenMissingET MissingET
+
   add Branch UniqueObjectFinder/jets Jet Jet
   add Branch UniqueObjectFinder/electrons Electron Electron
   add Branch UniqueObjectFinder/photons Photon Photon

=== modified file 'Template/Common/Cards/delphes_card_default.dat'
--- Template/Common/Cards/delphes_card_default.dat	2015-11-10 11:23:55 +0000
+++ Template/Common/Cards/delphes_card_default.dat	2016-09-07 23:55:24 +0000
@@ -31,6 +31,8 @@
 
   NeutrinoFilter
   GenJetFinder
+  GenMissingET
+  
   FastJetFinder
 
   JetEnergyScale
@@ -139,9 +141,9 @@
 
   # resolution formula for charged hadrons
   # based on arXiv:1405.6569
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
 }
 
 ###################################
@@ -156,9 +158,9 @@
 
   # resolution formula for electrons
   # based on arXiv:1405.6569
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
 }
 
 ###############################
@@ -172,9 +174,9 @@
   # set ResolutionFormula {resolution formula as a function of eta and pt}
 
   # resolution formula for muons
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*2.0e-4^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.02^2 + pt^2*3.0e-4^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*6.0e-4^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
 }
 
 ##############
@@ -266,8 +268,10 @@
   add EnergyFraction {3122} {0.3 0.7}
 
   # set ECalResolutionFormula {resolution formula as a function of eta and energy}
-  set ECalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.007^2 + energy*0.07^2 + 0.35^2)  +
-                             (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
+  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
+  set ECalResolutionFormula {                  (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
 
   # set HCalResolutionFormula {resolution formula as a function of eta and energy}
   set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
@@ -317,7 +321,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
 }
 
 #################
@@ -363,7 +367,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
 }
 
 #################
@@ -399,7 +403,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.25
 }
 
 ###################
@@ -463,6 +467,17 @@
   set JetPTMin 20.0
 }
 
+#########################
+# Gen Missing ET merger
+########################
+
+module Merger GenMissingET {
+# add InputArray InputArray
+  add InputArray NeutrinoFilter/filteredParticles
+  set MomentumOutputArray momentum
+}
+
+
 
 ############
 # Jet finder
@@ -526,7 +541,7 @@
   # based on arXiv:1211.4462
   
   # default efficiency formula (misidentification rate)
-  add EfficiencyFormula {0} {0.01+0.00038*pt}
+  add EfficiencyFormula {0} {0.01+0.000038*pt}
 
   # efficiency formula for c-jets (misidentification rate)
   add EfficiencyFormula {4} {0.25*tanh(0.018*pt)*(1/(1+ 0.0013*pt))}
@@ -591,6 +606,8 @@
   add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
 
   add Branch GenJetFinder/jets GenJet Jet
+  add Branch GenMissingET/momentum GenMissingET MissingET
+
   add Branch UniqueObjectFinder/jets Jet Jet
   add Branch UniqueObjectFinder/electrons Electron Electron
   add Branch UniqueObjectFinder/photons Photon Photon

=== removed file 'Template/Common/bin/internal/run_pythia'
--- Template/Common/bin/internal/run_pythia	2015-10-26 10:18:49 +0000
+++ Template/Common/bin/internal/run_pythia	1970-01-01 00:00:00 +0000
@@ -1,23 +0,0 @@
-#!/bin/bash
-#
-#  This runs pythia on the unweighted_events.dat
-#
-#  Usage: run_pythia [pydir]
-#  where pydir is the path to the pythia executable
-
-pydir=$1
-main=`pwd`
-
-if [ ! -e ../Cards/pythia_card.dat ]; then
-  echo "No pythia_card.dat found. Quitting..."
-  exit
-fi
-
-echo $$ >> ../myprocid
-# shower and hadronize event through Pythia
-#echo "" >> ../Cards/pythia_card.dat
-#echo "      LHAPATH=$pydir/PDFsets" >> ../Cards/pythia_card.dat
-export PDG_MASS_TBL=$pydir/mass_width_2004.mc
-rm -rf Events/hepmv_event.hepmc &> /dev/null
-$pydir/pythia && touch pythia.done
-

=== added file 'Template/LO/Cards/madanalysis5_hadron_card_default.dat'
--- Template/LO/Cards/madanalysis5_hadron_card_default.dat	1970-01-01 00:00:00 +0000
+++ Template/LO/Cards/madanalysis5_hadron_card_default.dat	2016-09-07 23:55:24 +0000
@@ -0,0 +1,3 @@
+# This card is used only if MA5 failed to create a default for this run
+# We therefore use as default: do nothing
+@MG5aMC skip_analysis
\ No newline at end of file

=== added file 'Template/LO/Cards/madanalysis5_parton_card_default.dat'
--- Template/LO/Cards/madanalysis5_parton_card_default.dat	1970-01-01 00:00:00 +0000
+++ Template/LO/Cards/madanalysis5_parton_card_default.dat	2016-09-07 23:55:24 +0000
@@ -0,0 +1,3 @@
+# This card is used only if MA5 failed to create a default for this run
+# We therefore use as default: do nothing
+@MG5aMC skip_analysis
\ No newline at end of file

=== 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-09-07 23:55:24 +0000
@@ -0,0 +1,77 @@
+!
+! Pythia8 cmd card automatically generated by MadGraph5_aMC@NLO
+! For more information on the use of the MG5aMC / Pythia8 interface, visit
+!    https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/LOPY8Merging
+!
+! ==================
+! General parameters 
+! ==================
+!
+Main:numberOfEvents      = -1
+!
+! -------------------------------------------------------------------
+! Specify the HEPMC output of the Pythia8 shower. You can set it to:
+!   auto      : MG5aMC will automatically place it the run_<i> directory
+!   /dev/null : to turn off the HEPMC output.
+!   <path>    : to select where the HEPMC file must written. It will 
+!               therefore not be placed in the run_<i> directory. The
+!               specified path, if not absolute, will be relative to 
+!               the Event/run_<i> directory of the process output.
+!   fifo      : to have MG5aMC setup the piping of the PY8 output to 
+!               analysis tools such as MadAnalysis5.  
+!   fifo@<fifo_path> :
+!               Same as 'fifo', but selecting a custom path to create the
+!               fifo pipe. (useful to select a mounted drive that supports 
+!               fifo). Note that the fifo file extension *must* be '.hepmc.fifo'.
+! -------------------------------------------------------------------
+!
+HEPMCoutput:file         = auto
+!
+! --------------------------------------------------------------------
+! Parameters relevant only when performing MLM merging, which can be
+! turned on by setting ickkw to '1' in the run_card and chosing a 
+! positive value for the parameter xqcut.
+! 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 CKKW-L merging, which can
+! be turned on by setting the parameter 'ptlund' *or* 'ktdurham' to
+! a positive value. 
+! 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 
+! ==========================
+!
+! By default, Pythia8 generates multi-parton interaction events. This is
+! often irrelevant for phenomenology and very slow. You can turn this 
+! feature off by uncommenting the line below if so desired.
+!partonlevel:mpi = off

=== 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-09-07 23:55:24 +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,55 @@
  %(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.        
-#*********************************************************************
-#**********************************************************
-#
-#**********************************************************
-# 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
-#**********************************************************
+ %(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
+#*********************************************************************
+#
+#*********************************************************************
+# handling of the helicities:
+#  0: sum over all helicities
+#  1: importance sampling over helicities
+#*********************************************************************
+   %(nhel)s  = nhel          ! using helicities importance sampling or not.
+#*********************************************************************
+# 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)
-#**********************************
+#
+#*********************************************************************
+# 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 
-#*************************************************************
-# 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
-#*******************                                                 
-#                                                                    
+#*********************************************************************
+# Standard Cuts                                                      *
 #*********************************************************************
 # Minimum and maximum pt's (for max, -1 means no cut)                *
 #*********************************************************************
@@ -246,21 +246,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                                 *
@@ -279,5 +278,6 @@
 %(sys_alpsfact)s = sys_alpsfact  # \alpha_s emission scale factors
 %(sys_matchscale)s = sys_matchscale # variation of merging scale
 # PDF sets and number of members (0 or none for all members).
-%(sys_pdf)s = sys_pdf # matching scales
+%(sys_pdf)s = sys_pdf # separate by && if more than one set.
 # MSTW2008nlo68cl.LHgrid 1  = sys_pdf
+#

=== added directory 'Template/LO/Source/BIAS'
=== 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-09-07 23:55:24 +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-09-07 23:55:24 +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-09-07 23:55:24 +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-09-07 23:55:24 +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-09-07 23:55:24 +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-09-07 23:55:24 +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-09-07 23:55:24 +0000
@@ -2,11 +2,16 @@
 DEFAULT_F_COMPILER=gfortran
 MACFLAG=-mmacosx-version-min=10.7
 DEFAULT_CPP_COMPILER=clang
+MG5AMC_VERSION=SpecifiedByMG5aMCAtRunTime
 STDLIB=-lc++
+PYTHIA8_PATH=NotInstalled
 STDLIB_FLAG=-stdlib=libc++
 #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 +95,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
\ No newline at end of file

=== 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-09-07 23:55:24 +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-09-07 23:55:24 +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,e16.7e3,3e15.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
@@ -162,13 +277,20 @@
 c      aqcd = g*g/4d0/pi
 
       write(lun,'(a)') '<event>'
-      write(lun,'(i2,i4,4e15.7)') nexternal,ievent,wgt,scale,aqed,aqcd
+      write(lun,'(i2,i4,e16.7e3,3e15.7)') nexternal,ievent,wgt,scale,aqed,aqcd
+      write(*,'(i2,i4,e15.7e3,3e15.7)') nexternal,ievent,wgt,scale,aqed,aqcd
       do i=1,nexternal
          write(lun,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))
       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-09-07 23:55:24 +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.
+     &                        is_a_j(i).or.is_a_b(i)) 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. 
+     &                          is_a_j(i).or.is_a_b(i)) 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/genps.f'
--- Template/LO/SubProcesses/genps.f	2016-06-16 15:03:59 +0000
+++ Template/LO/SubProcesses/genps.f	2016-09-07 23:55:24 +0000
@@ -216,9 +216,9 @@
             if(ebeam(1).lt.m1) ebeam(1)=m1
             if(ebeam(2).lt.m2) ebeam(2)=m2
             pi1(0)=ebeam(1)
-            pi1(3)=sqrt(ebeam(1)**2-m1**2)
+            pi1(3)=sqrt(max(ebeam(1)**2-m1**2, 0d0))
             pi2(0)=ebeam(2)
-            pi2(3)=-sqrt(ebeam(2)**2-m2**2)
+            pi2(3)=-sqrt(max(ebeam(2)**2-m2**2, 0d0))
             stot=m1**2+m2**2+2*(pi1(0)*pi2(0)-pi1(3)*pi2(3))
          endif
          write(*,'(x,a,f13.2)') 'Set CM energy to ',sqrt(stot)

=== modified file 'Template/LO/SubProcesses/makefile'
--- Template/LO/SubProcesses/makefile	2016-05-26 23:36:59 +0000
+++ Template/LO/SubProcesses/makefile	2016-09-07 23:55:24 +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-09-07 23:55:24 +0000
@@ -994,21 +994,21 @@
             pt2ijcl(nexternal-2)=q2fact(1)
             if(nexternal.gt.3) pt2ijcl(nexternal-3)=q2fact(1)
          else
-            q2fact(1)=pt2ijcl(nexternal-2)
-            q2fact(2)=q2fact(1)
+            q2fact(1)=scalefact**2*pt2ijcl(nexternal-2)
+            q2fact(2)=scalefact**2*q2fact(1)
          endif
       elseif(jcentral(1).eq.0)then
-            q2fact(1) = pt2ijcl(jfirst(1))
+            q2fact(1) = scalefact**2*pt2ijcl(jfirst(1))
       elseif(jcentral(2).eq.0)then
-            q2fact(2) = pt2ijcl(jfirst(2))
+            q2fact(2) = scalefact**2*pt2ijcl(jfirst(2))
       elseif(ickkw.eq.2.or.(pdfwgt.and.ickkw.gt.0))then
 c     Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
 c     f1(x1,pt2E) is given by DSIG, just need to set scale.
 c     Use the minimum scale found for fact scale in ME
          if(jlast(1).gt.0.and.jfirst(1).le.jlast(1))
-     $        q2fact(1)=min(pt2ijcl(jfirst(1)),q2fact(1))
+     $        q2fact(1)=scalefact**2*min(pt2ijcl(jfirst(1)),q2fact(1))
          if(jlast(2).gt.0.and.jfirst(2).le.jlast(2))
-     $        q2fact(2)=min(pt2ijcl(jfirst(2)),q2fact(2))
+     $        q2fact(2)=scalefact**2*min(pt2ijcl(jfirst(2)),q2fact(2))
       endif
 
 c     Check that factorization scale is >= 2 GeV
@@ -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-09-07 23:55:24 +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,17 @@
          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)
+     $       .and..not.from_decay(i) ) 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-09-07 23:55:24 +0000
@@ -150,6 +150,10 @@
       denom = pin(0)*pout(0) - pin(3)*pout(3)
       if (denom.ne.0d0) then
          get_betaz = (pin(3) * pout(0) - pout(3) * pin(0)) / denom
+      else if (pin(0).eq.pin(3)) then
+         get_betaz = (pin(0)**2 - pout(0)**2)/(pin(0)**2 + pout(0)**2)
+      else if (pin(0).eq.abs(pin(3))) then
+         get_betaz = (pout(0)**2 - pin(0)**2)/(pin(0)**2 + pout(0)**2)
       else
          get_betaz = 0d0
       endif
@@ -236,9 +240,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)
@@ -336,7 +340,8 @@
 c
       xsum = 0d0
       i = nw
-      do while (xsum-dabs(swgt(i))*(nw-i) .lt. xtot*trunc_max .and. i .gt. 2)
+      do while (xsum-dabs(swgt(i))*(nw-i) .lt. xtot*trunc_max
+     $ .and. i .gt. 2)
          xsum = xsum + dabs(swgt(i))
          i = i-1
       enddo
@@ -437,7 +442,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 +465,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 +483,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 +511,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 +527,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
@@ -596,53 +601,45 @@
       pboost(2)=0d0
       pboost(3)=0d0
       if (nincoming.eq.2) then
-         if (xbk(1) .gt. 0d0 .and. xbk(1) .lt. 1d0 .and.
-     $       xbk(2) .gt. 0d0 .and. xbk(2) .lt. 1d0) then
-            pboost(0) = -1d0 ! special flag
-            ! construct the beam momenta in each frame and compute the related (z)boost
-            ebi(0) = p(0,1)/xbk(1)
-            ebi(1) = 0
-            ebi(2) = 0
-            ebi(3) = DSQRT(ebi(0)**2-m1**2)
-            ebo(0) = ebeam(1)
-            ebo(1) = 0
-            ebo(2) = 0
-            ebo(3) = DSQRT(ebo(0)**2-m1**2)
-            beta = get_betaz(ebi, ebo)
-c           eta = sqrt(xbk(1)*ebeam(1)/
-c     $                 (xbk(2)*ebeam(2)))
-c            pboost(0)=p(0,1)*(eta + 1d0/eta)
-c            pboost(3)=p(0,1)*(eta - 1d0/eta)
-         else if (xbk(1) .gt. 0d0 .and. xbk(1) .lt. 1d0 .and.
-     $       xbk(2) .eq. 1d0) then
-            pboost(0)=xbk(1)*ebeam(1)+ebeam(2)
-            pboost(3)=xbk(1)*ebeam(1)-
-     $           sqrt(ebeam(2)**2-m2**2)
-         else if (xbk(1) .eq. 1d0 .and.
-     $       xbk(2) .gt. 0d0 .and. xbk(2) .lt. 1d0) then
-            pboost(0)=ebeam(1)+xbk(2)*ebeam(2)
-            pboost(3)=sqrt(ebeam(1)**2-m1**2)-
-     $           xbk(2)*ebeam(2)
-         else if (xbk(1) .eq. 1d0 .and. xbk(2) .eq. 1d0) then
-            pboost(0)=ebeam(1)+ebeam(2)
-            pboost(3)=sqrt(ebeam(1)**2-m1**2)-
-     $           sqrt(ebeam(2)**2-m2**2)
+         if (xbk(1) .gt. 0d0 .and. xbk(1) .le. 1d0 .and.
+     $       xbk(2) .gt. 0d0 .and. xbk(2) .le. 1d0) then
+            if(xbk(1).eq.1d0.or.pmass(1).eq.0d0) then
+               ! construct the beam momenta in each frame and compute the related (z)boost
+               ebi(0) = p(0,1)/xbk(1) ! this assumes that particle 1 is massless or mass equal to beam
+               ebi(1) = 0
+               ebi(2) = 0
+               ebi(3) = DSQRT(ebi(0)**2-m1**2)
+               ebo(0) = ebeam(1)
+               ebo(1) = 0
+               ebo(2) = 0
+               ebo(3) = DSQRT(ebo(0)**2-m1**2)
+               beta = get_betaz(ebi, ebo)
+            else
+               ebi(0) = p(0,2)/xbk(2) ! this assumes that particle 2 is massless or mass equal to beam
+               ebi(1) = 0
+               ebi(2) = 0
+               ebi(3) = -1d0*DSQRT(ebi(0)**2-m2**2)
+               ebo(0) = ebeam(2)
+               ebo(1) = 0
+               ebo(2) = 0
+               ebo(3) = -1d0*DSQRT(ebo(0)**2-m2**2)
+               beta = get_betaz(ebi, ebo)
+               ! wrong boost if both parton are massive!
+            endif
          else
             write(*,*) 'Warning bad x1 or x2 in write_leshouche',
      $           xbk(1),xbk(2)
          endif
-      endif
-      if (pboost(0).ne.-1d0) then
+         do j=1,nexternal
+            call zboost_with_beta(p(0,j),beta,pb(0,isym(j,jsym)))
+            pb(4,isym(j,jsym))=pmass(j)
+         enddo
+      else
          do j=1,nexternal
             call boostx(p(0,j),pboost,pb(0,isym(j,jsym)))
             ! Add mass information in pb(4)
             pb(4,isym(j,jsym))=pmass(j)
          enddo
-      else
-         do j=1,nexternal
-            call zboost_with_beta(p(0,j),beta,pb(0,isym(j,jsym)))
-            pb(4,isym(j,jsym))=pmass(j)
-         enddo
       endif
 c
 c     Add info on resonant mothers
@@ -746,6 +743,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-09-07 23:55:24 +0000
@@ -9,6 +9,7 @@
 run=$2
 tag=$3
 cross=$4 
+file=$5
 
 main=`pwd`
 
@@ -22,17 +23,35 @@
   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
+#    gzip ${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-09-07 23:55:24 +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-09-07 23:55:24 +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-09-07 23:55:24 +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-09-07 23:55:24 +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