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

[Olivier Mattelaer]

Hi,

I think that it is now time to release this branch.
Any objections? Know bugs that need to be fixed? Any important merge that we have to wait for?

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/2.5.1]

330. By Olivier Mattelaer on 2016-10-25

allow PYTHONPATH to not have valid path inside it

331. By Olivier Mattelaer on 2016-10-25

fix problem with plugin for ./bin/generate_events script

332. By Olivier Mattelaer on 2016-10-26

avoid problem with xml in the run_card

333. By Paolo Torrielli on 2016-10-27

Merging with 2.3.4_new_HWPP: interface to HW7 available from this merge

[Rikkert Frederix (frederix)]

Hi Olivier,

Go ahead!

Cheers,
Rik

review: Approve

[lp:~maddevelopers/mg5amcnlo/2.5.1]

334. By Olivier Mattelaer on 2016-10-31

include antiname in name2part

[Olivier Mattelaer (olivier-mattelaer)]

I'm waiting the green light of Valentin.
He asked me to wait for one bug in the PY8 interface to be fixed before releasing this.
In principle, this should be fixed in less than 24h

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/2.5.1]

335. By Olivier Mattelaer on 2016-11-01

fix systematics for e+ e- initial state + change web adress for uiuc

336. By Olivier Mattelaer on 2016-11-02

1) allow py8 to run in single core
2) add MA5 card to the detect_card_type algorithm (to allow to specify a path to a file)

337. By Valentin Hirschi on 2016-11-03

1. Fixed the presence of 'Merging:Process' when there is no merging.
2. Fixed an issue of nJetMax = -1 because of the new behaviour of user_set.
3. added the new hidden parameter HEPMCOutput:scaling, automatically set to 1.0e9 so as to have it in picobarn.

338. By Valentin Hirschi on 2016-11-03

1. Fixed an issue of HEPMC weight normalization when using PY8 parallelization.

339. By Olivier Mattelaer on 2016-11-03

1) Fix a problem with plugin directory specified via PYTHONPATH
2) add two hidden paramereter in the run_card (both LO and NLO)
systematics_program [none, syscalc, systematics, auto]
systematics_arguments
to allow to run systematics on the flight

[Olivier Mattelaer (olivier-mattelaer)]

Ok I'm starting the process of merging the branch.
The release should be done by tomorrow.

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/2.5.1]

340. By Olivier Mattelaer on 2016-11-03

put systematics to OFF for LO interference processes

341. By Valentin Hirschi on 2016-11-04

1. Fixed parallelization of PY8 when enforcing a number of events to be processed while keeping the excessive events so that they can be used if they are skipped.

342. By Olivier Mattelaer on 2016-11-04

VERSION/UpdateNotes + small fix from Valentin

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2016-07-06 08:54:04 +0000
+++ MadSpin/decay.py	2016-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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 (in mm) below which the invariant livetime 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,27 +246,24 @@
 #*********************************************************************
  %(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                                 *
-# WARNING: If use_syst is T, matched Pythia output is                *
-#          meaningful ONLY if plotted taking matchscale              *
-#          reweighting into account!                                 *
+# WARNING: Do not use for interference type of computation           *
 #*********************************************************************
    %(use_syst)s  = use_syst      ! Enable systematics studies
 #
@@ -279,5 +276,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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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/myamp.f'
--- Template/LO/SubProcesses/myamp.f	2016-06-16 15:03:59 +0000
+++ Template/LO/SubProcesses/myamp.f	2016-11-04 09:43:15 +0000
@@ -212,7 +212,8 @@
 c
             onshell = (abs(xmass - prmass(i,iconfig)) .lt.
      $           bwcutoff*prwidth(i,iconfig).and.
-     $           prwidth(i,iconfig)/prmass(i,iconfig).lt.0.1d0)
+     $           (prwidth(i,iconfig)/prmass(i,iconfig).lt.0.1d0.or.
+     $            gForceBW(i,iconfig).eq.1))
             if(onshell)then
 c     Remove on-shell forbidden s-channels (gForceBW=2) (JA 2/10/11)
               if(gForceBW(i,iconfig).eq.2) then

=== 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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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-11-04 09:43:15 +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

=== removed file 'Template/NLO/MCatNLO/HWPPAnalyzer/HepMCFortran.h'
--- Template/NLO/MCatNLO/HWPPAnalyzer/HepMCFortran.h	2013-09-26 09:34:42 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/HepMCFortran.h	1970-01-01 00:00:00 +0000
@@ -1,154 +0,0 @@
-// -*- C++ -*-
-#ifndef MCATNLO_HepMCFortran_H
-#define MCATNLO_HepMCFortran_H
-//
-// This is the declaration of the HepMCFortran class.
-//
-
-#include "ThePEG/Handlers/AnalysisHandler.h"
-#include "Herwig++/Utilities/Histogram.h"
-
-#include "HepMC/IO_BaseClass.h"
-
-#include "HepMCFortran.fh"
-
-namespace MCatNLO {
-
-using namespace ThePEG;
-using namespace Herwig;
-
-/**
- * Here is the documentation of the HepMCFortran class.
- *
- * @see \ref HepMCFortranInterfaces "The interfaces"
- * defined for HepMCFortran.
- */
-class HepMCFortran: public AnalysisHandler {
-
-public:
-
-  /** @name Virtual functions required by the AnalysisHandler class. */
-  //@{
-  /**
-   * Analyze a given Event. Note that a fully generated event
-   * may be presented several times, if it has been manipulated in
-   * between. The default version of this function will call transform
-   * to make a lorentz transformation of the whole event, then extract
-   * all final state particles and call analyze(tPVector) of this
-   * analysis object and those of all associated analysis objects. The
-   * default version will not, however, do anything on events which
-   * have not been fully generated, or have been manipulated in any
-   * way.
-   * @param event pointer to the Event to be analyzed.
-   * @param ieve the event number.
-   * @param loop the number of times this event has been presented.
-   * If negative the event is now fully generated.
-   * @param state a number different from zero if the event has been
-   * manipulated in some way since it was last presented.
-   */
-  virtual void analyze(tEventPtr event, long ieve, int loop, int state);
-  //@}
-
-public:
-
-  /**
-   * The standard Init function used to initialize the interfaces.
-   * Called exactly once for each class by the class description system
-   * before the main function starts or
-   * when this class is dynamically loaded.
-   */
-  static void Init();
-
-protected:
-
-  /** @name Standard Interfaced functions. */
-  //@{
-  /**
-   * Initialize this object. Called in the run phase just before
-   * a run begins.
-   */
-  virtual void doinitrun();
-
-  /**
-   * Finalize this object. Called in the run phase just after a
-   * run has ended. Used eg. to write out statistics.
-   */
-  virtual void dofinish();
-  //@}
-
-protected:
-
-  /** @name Clone Methods. */
-  //@{
-  /**
-   * Make a simple clone of this object.
-   * @return a pointer to the new object.
-   */
-  virtual IBPtr clone() const;
-
-  /** Make a clone of this object, possibly modifying the cloned object
-   * to make it sane.
-   * @return a pointer to the new object.
-   */
-  virtual IBPtr fullclone() const;
-  //@}
-
-private:
-
-  /**
-   * The static object used to initialize the description of this class.
-   * Indicates that this is an concrete class without persistent data.
-   */
-  static NoPIOClassDescription<HepMCFortran> initHepMCFortran;
-
-  /**
-   * The assignment operator is private and must never be called.
-   * In fact, it should not even be implemented.
-   */
-  HepMCFortran & operator=(const HepMCFortran &);
-
-
-private:
-
-  HepMC::IO_BaseClass *_hepevtio;
-
-};
-
-}
-
-#include "ThePEG/Utilities/ClassTraits.h"
-
-namespace ThePEG {
-
-/** @cond TRAITSPECIALIZATIONS */
-
-/** This template specialization informs ThePEG about the
- *  base classes of HepMCFortran. */
-template <>
-struct BaseClassTrait<MCatNLO::HepMCFortran,1> {
-  /** Typedef of the first base class of HepMCFortran. */
-  typedef AnalysisHandler NthBase;
-};
-
-/** This template specialization informs ThePEG about the name of
- *  the HepMCFortran class and the shared object where it is defined. */
-template <>
-struct ClassTraits<MCatNLO::HepMCFortran>
-  : public ClassTraitsBase<MCatNLO::HepMCFortran> {
-  /** Return a platform-independent class name */
-  static string className() { return "MCatNLO::HepMCFortran"; }
-  /**
-   * The name of a file containing the dynamic library where the class
-   * HepMCFortran is implemented. It may also include several, space-separated,
-   * libraries if the class HepMCFortran depends on other classes (base classes
-   * excepted). In this case the listed libraries will be dynamically
-   * linked in the order they are specified.
-   */
-  static string library() { return "HepMCFortran.so"; }
-};
-
-/** @endcond */
-
-}
-
-#endif /* MCATNLO_HepMCFortran_H */

=== added file 'Template/NLO/MCatNLO/HWPPAnalyzer/HepMCFortran2.h'
--- Template/NLO/MCatNLO/HWPPAnalyzer/HepMCFortran2.h	1970-01-01 00:00:00 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/HepMCFortran2.h	2016-11-04 09:43:15 +0000
@@ -0,0 +1,154 @@
+// -*- C++ -*-
+#ifndef MCATNLO_HepMCFortran_H
+#define MCATNLO_HepMCFortran_H
+//
+// This is the declaration of the HepMCFortran class.
+//
+
+#include "ThePEG/Handlers/AnalysisHandler.h"
+#include "Herwig++/Utilities/Histogram.h"
+
+#include "HepMC/IO_BaseClass.h"
+
+#include "HepMCFortran.fh"
+
+namespace MCatNLO {
+
+using namespace ThePEG;
+using namespace Herwig;
+
+/**
+ * Here is the documentation of the HepMCFortran class.
+ *
+ * @see \ref HepMCFortranInterfaces "The interfaces"
+ * defined for HepMCFortran.
+ */
+class HepMCFortran: public AnalysisHandler {
+
+public:
+