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

[Olivier Mattelaer]

[lp:~maddevelopers/mg5amcnlo/plugin_plus_py8]

441. By Olivier Mattelaer on 2016-08-25

fix with latest 2.5.0

[Olivier Mattelaer (olivier-mattelaer)]

Here is a first partial run on the code changes:

1) some comment about the pythia8_card.
1.1) Need to explain what HEPMCoutput:file = auto means.
1.2) I do not understand what you mean by
        ! <path> : to select where the HEPMC file must written. It will
        ! therefore *not* be available in the process output.
     Could you rephrase
1.3) The line JetMatching:nJetMax = -1 is present two times. What happens if you modify only one?
1.4) How do we choose which of the MLM/CKKW is running (would be nice to be written)?
1.5) The link to the wiki explaining the MG/PY interface should be written in the cards.

2) somme comment about the run card
2.1) wtf is this:
# 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.

3)some bad change has been done to Template/LO/SubProcesses/cuts.f
forbidding review of the file (100% diff) -> need to be done manually later.

4) in setcuts.f I see:
+c Remember mergeable particles
+ do j=1,pdgs_for_merging_cut(0)
+ if ( pdgs_for_merging_cut(j) .ne. 0
+ $ .and. abs(idup(i,1,iproc)) .eq.pdgs_for_merging_cut(j) ) then
+ is_pdg_for_merging_cut(i)=.true.
+ exit
+ endif
+ enddo
should not we set this on True only if from_decay(i)=false ?

5) I guess that the following is wrong:
amcatnlo_run_interface line 3306:
        # when are we force to change the tag new_run:previous run requiring changes
        upgrade_tag = {'parton': ['parton','pythia','pgs','delphes','shower'],
                       'pythia': ['pythia','pgs','delphes'],
                       'shower': ['shower'],
                       'pgs': ['pgs'],
                       'delphes':['delphes'],
                       'madanalysis5_hadron':['madanalysis5_hadron'],
                       'plot':[]}
your change did not make sense to me (neither why those lines are in AMC@NLO_run_interface).
This will be the source for bug I guess.

need to continue at file
=== modified file 'madgraph/interface/common_run_interface.py'

will continue as soon as possible.

review: Needs Fixing

[Olivier Mattelaer (olivier-mattelaer)]

ok I continue

6) in common_run_interface in the function check_madanalysis5

+ # Maybe the user used something like 'logging.INFO'
+ MA5_options['MA5_stdout_lvl']=eval(lvl)
+ except:
+ try:
+ MA5_options['MA5_stdout_lvl']=eval('logging.%s'%lvl)

would be better to not use eval (this always creates security for the web interface).
Rather use (in both case) getattr(logging, lvl)

7) a bit lower you have the line:
+ logger.warning("LHE event file not found in \n%s\ns"%input_file+
+ "AboringParton level MA5 analysis is disabled.")
what does that mean?

8) a bit below you have
+ file_candidates = glob.glob(pjoin(self.me_dir,'Events',self.run_name,htag))+\
+ glob.glob(pjoin(self.me_dir,'Events',self.run_name,'%s.gz'%htag))

please replace this be call like
 misc.glob('%s.gz'%htag, pjoin(self.me_dir,'Events',self.run_name))
This is done in order to strenghten the security of the function (and limit bug)
Please check that no direct call to glob.glob are remaining in the file. (I have seen many more)

9) need to check for later: the code of ask_madanalysis5_run_configuration seems weird. Need to ensure to check the associate functionality

10) concerning run_madanalysis5:

10.1)remove the
############################# TEMPORARY ########################################
+# Here's how to print the MA5 commands generated by MG5aMC
+# for MA5_runtag, MA5_cmds in MA5_cmds_list:
+# misc.sprint('****************************************')
+# misc.sprint('* Commands for MA5 runtag %s:'%MA5_runtag)
+# misc.sprint('\n'+('\n'.join('* %s'%cmd for cmd in MA5_cmds)))
+# misc.sprint('****************************************')
+############################# TEMPORARY ########################################
(or at least change the temporary)

10.2) Please use the same standard as for other program for the following (one line in bold) [We are not an advertising company]
+ logger.info('-------------------------------------------','$MG:color:GREEN')
+ logger.info('Now starting MadAnalysis5 [arXiv:1206.1599]','$MG:color:GREEN')
+ logger.info('-------------------------------------------','$MG:color:GREEN')

10.3) concerning
+ if MA5_runtag.startswith('_reco_'):
+ # When doing a reconstruction we must first copy the event file
I do not see any cp. so I'm not sure if the comment is accurate.
If it is, please try to avoid to copy event since the filesystem can (often) be quite slow (especially on cluster)

10.4) concerning this:
+ shutil.copy(target, pjoin(self.me_dir,'Events',self.run_name,carboncopy_name))
Is it possible to specify a path directly in MA5 to not have to move the file.
(multiple output are like this actually in that routine --the full HTML for example--)

11) In do_delphes/check_delphes.
you have remove the use of the multicore method to unzip the file, while editing the question.
I can see an advantges on that for condor cluster if cluster_temp_directory is set but in multicore I do not see the point. What is the point?

12) in do_delphes the line
+ if os.path....

[Olivier Mattelaer (olivier-mattelaer)]

Let's go for part III:

14) in madgraph/interface/extended_cmd.py
I see a misc.sprint remaining.
@@ -595,20 +601,22 @@
                 to_rm = len(self.completion_prefix) - 1
                 Nbegidx = len(line.rsplit(os.path.sep, 1)[0]) + 1
                 data = compfunc(Ntext.replace('\ ', ' '), line, Nbegidx, endidx)
+ misc.sprint(data)

15) Note for the functionality review:
Potentially deep change in the auto-completion behavior: Need extensive testing of all auto-completion to check if they still work. (use the opportunity to add unitest on the auto-completion)
- # correct wrong splitting with '-'
- elif line and line[begidx-1] == '-':
+ # correct wrong splitting with '-'/"="
+ elif line and line[begidx-1] in ['-',"=",':']:

16) in extended_cmd, a new parameter self.case is added to the question interfase/
the name is not clear enough if the question is then case sensitive or case insensitive.
Therefore the name has to be modified.

17) in madgraph/interface/launch_ext_program.py @@ -766,7 +766,7 @@
can you check the following change: (the print is kept but not the compilation)
         # Make pythia8
         print "Running make for pythia8 directory"
- misc.compile(cwd=os.path.join(self.running_dir, os.path.pardir), mode='cpp')
+ #misc.compile(cwd=os.path.join(self.running_dir, os.path.pardir), mode='cpp')

18) in madevent_interface for the check_pythia8 command

18.1) you have the line
+ # If not pythia-pgs path

18.2) you have
+ mode = laststep[0].split('=')[-1]
+ if mode not in ['auto', 'pythia', 'delphes']:
should not you also support laststep=pythia8

18.3) The following line are not bugged (output_file not define when untarring:
+ input_file = pjoin(self.me_dir,'Events',self.run_name, 'unweighted_events.lhe')
+ #output_file = pjoin(self.me_dir, 'Events', 'unweighted_events.lhe.gz')
+ if not os.path.exists('%s.gz'%input_file):
+ if os.path.exists(input_file):
+ misc.gzip(input_file, keep=True, stdout=output_file)
+ else:
+ raise self.InvalidCmd('No event file corresponding to %s run. '
+ % self.run_name)

19) complete_shower
This use a useless (and dangerous eval)
+ elif len(args)>1 and args[1] in self._interfaced_showers:
+ return eval("self.complete_%s(text,'%s',%d,%d)"%
+ (args[1],line.replace(args[0]+' ',''),begidx-len(args[0])-1,
+ endidx-len(args[0])-1))

this should be replaced by a nicer (splitting in multiple will also be more readable)
return getattr(self, 'complete_%s' % text)(text, args[1],line.replace(args[0]+' ',''), begidx-len(args[0])-1, endidx-len(args[0])-1))

20) in complete_pythia8:
apply also the replacement for glob.glob in this file/function. (if not done already)

21) in setup_pythia8
Not sure to understand this:
+ # Use defaultSet not to overwrite the current userSet status
+ PY8_Card.defaultSet('HEPMCoutput:file','PY8_hepmc.fifo')
This line never do...

[Olivier Mattelaer (olivier-mattelaer)]

and now the final part for the code review.

35) in banner.py
+ self.add_param('colored_pdgs', '[1,2,3,4,5]')
Why not having it as a list?

36) Valentin are you working on PC? or maybe you like too much xml ;-)
[in banner.py]
# to indicate that he wants to pipe the output. Or \dev\null to turn the

37) in lhe_parser.py: the function do not correspond to the associated comment:
+
+ def unweight(self, outputpath, get_wgt=None, max_wgt=0, trunc_error=0,
+ event_target=0, log_level=logging.INFO, normalization='average'):
         """unweight the current file according to wgt information wgt.
         which can either be a fct of the event or a tag in the rwgt list.
         max_wgt allow to do partial unweighting.
         trunc_error allow for dynamical partial unweighting
         event_target reweight for that many event with maximal trunc_error.
         (stop to write event when target is reached)
+ 'strategy' defines the normalisation mode: 4 -> average. 3-> sum
         """
         if not get_wgt:

38) Why not put this in madevent_interface.py? not really the point of misc for such function.
This can be a function like in misc. madgraph_interface can obviously call it.

+# Return a warning (if applicable) on the consistency of the current Pythia8
+# and MG5_aMC version specified. It is placed here because it should be accessible
+# from both madgraph5_interface and madevent_interface
+#===============================================================================
+def mg5amc_py8_interface_consistency_warning(options):

39)
+def get_MadAnalysis5_interpreter(mg5_path, ma5_path, mg5_interface=None,
+ logstream = sys.stdout, loglevel = logging.INFO, forced = True):

39.1) This function need to work with a with statement.
      To always ensure that the variable are restored at the end.
      or at least restore the variable in case of crash
39.2) the try/except is much too large. many stuff can crash here.
39.3) the import readline might not pass on all machine. (espcecially on VM) so this need to be handle as well (stupid to forbid MA5 on VM right)

To follow the status on those answer, I have put a page on the wiki:
https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/DevelopmentPage/ReviewProgress

[Olivier Mattelaer (olivier-mattelaer)]

Let move forward on the point 3 of the review:
reviewing the change in cuts.f (not done before to the a windows -> linux conversion on end of line)

3.1) The following is at least weird
C COUNT NUMBER OF MASSLESS OUTGOING PARTICLES, SINCE WE DO NOT WANT
C TO APPLY A CUT FOR JUST A SINGLE MASSIVE PARTICLE IN THE FINAL STATE.
        NMASSLESS = 0
        DO I=NINCOMING+1,NEXTERNAL
          if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) .and.
     & DSQRT(DOT( P(0,I), P(0,I))) .LT. 5D0) THEN
            NMASSLESS = NMASSLESS + 1
          ENDIF
        ENDDO
this seems to count in 4F and 5F the b as massles (which sounds bad to me)
should not we use maxjetflavor for this?

3.2) some comment as 3.1 but for ptlund.
C PROCESS WITH EXACTLY TWO MASSLESS OUTGOING PARTICLES IS SPECIAL
C BECAUSE AN ENERGY-SHARING VARIABLE LIKE "Z" DOES NOT MAKE SENSE.
C IN THIS CASE, ONLY APPLY MINIMAL pT W.R.T BEAM CUT.
C THIS CUT WILL ONLY APPLY TO THE TWO-MASSLESS PARTICLE STATE.
        NMASSLESS = 0
        DO I=NINCOMING+1,NEXTERNAL
          if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) .and.
     & DSQRT(DOT( P(0,I), P(0,I))) .LT. 5D0) THEN
            NMASSLESS = NMASSLESS + 1
          ENDIF
        ENDDO

So I marked 3 done in the wiki and add two line in the wiki for those two (identical) point

Otherwise here is two additional point

40) the edition of the card does not work

41) the auto-completion of install MadAnalysis does not work after the first option.
(try to auto-complete install MadAnalysis5 --no_root_in_MA5)

42) the instalation of MA5 fails without the option "--no_root_in_MA5 --no_MA5_further_install"
would be nice to either automatically switch to those options or at least suggest them to the user when it fails

[lp:~maddevelopers/mg5amcnlo/plugin_plus_py8]

442. By Olivier Mattelaer on 2016-08-26

make some of the change found when doing the review

443. By Olivier Mattelaer on 2016-08-27

improve the help command of the card edition interface
- support pythia8 (and madloop)
- parameter can have an optional parameter comment
- this parameter is printed if the user do help on that parameter

[Olivier Mattelaer (olivier-mattelaer)]

43) some acceptance test are failing (related to self.proc_defs @NLO)

44) The question for madanalysis should be changed for the display to the user.
44.1)I would propose:

The following switches determine which programs are run:
 1 Choose the shower/hadronization program: shower = OFF
 2 Choose the detector simulation program: detector = OFF
 3 Choose which the plotting/analysis program: analysis = OFF
 3 Decay particles with the MadSpin module: madspin = OFF
 4 Add weights to the events based on changing model parameters: reweight = OFF
  Either type the switch number (1 to 5) to change its setting,
  Set any switch explicitly (e.g. type 'madspin=ON' at the prompt)
  Type 'help' for the list of all valid option
  Type '0', 'auto', 'done' or just press enter when you are done.

Such that
a) we can add other package easily. (like the fortran analysis present in MC@NLO, or interface to Falstlim/Rivet/...)
b) Such that we can choose to use either MA4 or MA5

44.2) Also, contrary to the other options, the default should be to have the plotting routine activated when installed. (at least for the parton level)

44.3) I would propose to have only MA5 as option (equivalent to parton+hadron if shower is run)
      (potentially have a way in the madanalysis_hadron card to specify that you want to bypass it if needed)

44.4) Neither madanalysis5_parton/hadron have an help message working with this question.

[lp:~maddevelopers/mg5amcnlo/plugin_plus_py8]

444. By Olivier Mattelaer on 2016-08-28

merge with latest 2.5.0

445. By Olivier Mattelaer on 2016-08-28

moving forward in the review

446. By Olivier Mattelaer on 2016-08-28

additional fix related to the review

447. By Valentin Hirschi on 2016-08-29

1. Fixed a path issue.
2. Prevented any potential border effect in the additional_options argument of install and advanced_install

448. By Valentin Hirschi on 2016-08-29

1. Safeguarded against border effects in additional_options of advanced_install
2. Removed the random doubling of '--logging=' which caused troubles.

449. By Valentin Hirschi on 2016-08-29

1. removed left-over print statemnt.

450. By Olivier Mattelaer on 2016-08-29

fix review number 5, 37 and 41

451. By Olivier Mattelaer on 2016-08-29

fix point of the review #12

452. By Valentin Hirschi on 2016-08-30

1. modified the steering of MA5 so that now the switch is of the 'analysis' type.

453. By Valentin Hirschi on 2016-08-30

1. Merged with latest version of 2.5.0 with e+ e- systematics

454. By Valentin Hirschi on 2016-08-30

1. Treated points 1-23 of the review.

455. By Valentin Hirschi on 2016-08-30

1. Minor fix.

[Olivier Mattelaer (olivier-mattelaer)]

OK here is the final part of my review (mainly concerning MA5/interface)

45) install MadAnalysis command should install MA5 and not MA4
    a new command "install MadAnalysis4" should be include for installing MA4

46) the log of MA5 should be written in a file and not on screen.

47) The auto-completion is now broken by the following command:
launch PROC_sm_0 - i
quit
#still has MadEvent completion command

48) the auto-completion of do_xxx returns now multiple times the same function

49) when running only parton level (no shower) an empty line in the html "hadron_MA5" is present...

50) the launch command (and launch -i) is not put in the history of command (so not accessible with up arrow)

51) Having MA5 running (without crash) both before and after the parton-shower. (lot of progress in that directly! but still work requested.)

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/plugin_plus_py8]

456. By Valentin Hirschi on 2016-08-30

1. changed default value of sys_matchscale

457. By Valentin Hirschi on 2016-08-30

1. Esthetics changes

458. By Olivier Mattelaer on 2016-08-30

fix for point #17,24,47,48

459. By Olivier Mattelaer on 2016-08-30

try to fix point 3 and 4. (plus some other small fix not in the review)

460. By Valentin Hirschi on 2016-08-31

1. Several improvements on MA5 steering.

461. By Valentin Hirschi on 2016-08-31

1. Changed install MadAnalysis to MadAnalysis4

462. By Valentin Hirschi on 2016-08-31

1. Removed one consistency_warning_function copy in misc.

463. By Valentin Hirschi on 2016-08-31

1. Attended to the last points of the review

464. By Olivier Mattelaer on 2016-08-31

fixing a unittest

465. By Valentin Hirschi on 2016-08-31

1. Fixed a couple tests

466. By Olivier Mattelaer on 2016-08-31

fix review point #44 + fix ConfigFile/update tests

467. By Olivier Mattelaer on 2016-08-31

fix py8 html display. force sum for py6, improve status of py8

468. By Olivier Mattelaer on 2016-08-31

pass via status_update for reporting that MA5 is running

469. By Valentin Hirschi on 2016-08-31

1. Fixed the unit tests test_w_nlo_gen

470. By Valentin Hirschi on 2016-08-31

1. Small fixes related to what was seen during various tests with Olivier.

471. By Olivier Mattelaer on 2016-08-31

fixing multiple bug found thanks to the acceptance test created_matched_plot

472. By Valentin Hirschi on 2016-08-31

1. Fix MG5 crash when MA5 fail to start when generating cards.

473. By Olivier Mattelaer on 2016-08-31

remove a debug statement

474. By Olivier Mattelaer on 2016-09-01

fixing test related to the automatic consistency of the cards. Checking in detail the new behavior to ensure that no side effect has been introduced

475. By Valentin Hirschi on 2016-09-01

1. Fixed the missing process list when doing nlo parrallel generations.
2. Fixed two parallel test.

476. By Olivier Mattelaer on 2016-09-01

allow question for installation to recognise pre-existing part of installaton

477. By Valentin Hirschi on 2016-09-01

1. Merged witht he bias branch

478. By Olivier Mattelaer on 2016-09-01

modification for the new strategy of installation of MA5

479. By Olivier Mattelaer on 2016-09-01

update IOTEST

480. By Olivier Mattelaer on 2016-09-01

add MA5 banner in MA5 log

481. By Valentin Hirschi on 2016-09-01

1. updated offline tarballs and added two tests to test them.
2. Updated the do_install function features so as to comply with the new HEPToolsInstaller behavior for MA5.

482. By Valentin Hirschi on 2016-09-01

1. Change the default name of the fifo pipe to comply with MadAnalysis need of having the .hepmc extension.

483. By Valentin Hirschi on 2016-09-01

1. Removed a left over test code.

484. By Valentin Hirschi on 2016-09-01

1. Implemented support for the MA5 card tag '@MG5aMC skip_analysis'
2. Improved the printout of the MG5aMC distrib short parrallel tests.

485. By Olivier Mattelaer on 2016-09-01

fixing the short_parralel tests

486. By Valentin Hirschi on 2016-09-01

1. Improved fifo comments in card.
2. Made sure MA5 runs on one fifo only.

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

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

=== modified file 'Template/Common/Cards/delphes_card_default.dat'
--- Template/Common/Cards/delphes_card_default.dat	2015-11-10 11:23:55 +0000
+++ Template/Common/Cards/delphes_card_default.dat	2016-09-01 23:40:50 +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

=== added file 'Template/LO/Cards/pythia8_card_default.dat'
--- Template/LO/Cards/pythia8_card_default.dat	1970-01-01 00:00:00 +0000
+++ Template/LO/Cards/pythia8_card_default.dat	2016-09-01 23:40:50 +0000
@@ -0,0 +1,72 @@
+!
+! 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 
+! ==========================

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

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

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

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

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

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

=== added file 'Template/LO/Source/lhe_event_infos.inc'
--- Template/LO/Source/lhe_event_infos.inc	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/lhe_event_infos.inc	2016-09-01 23:40:50 +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-08-26 13:55:42 +0000
+++ Template/LO/Source/make_opts	2016-09-01 23:40:50 +0000
@@ -2,11 +2,16 @@
 DEFAULT_F_COMPILER=gfortran
 MACFLAG=-mmacosx-version-min=10.7
 DEFAULT_CPP_COMPILER=clang
+MG5AMC_VERSION=SpecifiedByMG5aMCAtRunTime
 STDLIB=-lc++
+PYTHIA8_PATH=NotInstalled
 STDLIB_FLAG=-stdlib=libc++
 #end_of_make_opts_variables
+
+BIASLIBDIR=../../../lib/
+BIASLIBRARY=libbias.$(libext)
+
 # Rest of the makefile
-
 ifeq ($(origin FFLAGS),undefined)
 FFLAGS= -O -w -fbounds-check -fPIC
 #FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
@@ -90,4 +95,9 @@
 else
 alfas_functions=alfas_functions
 llhapdf=
-endif
\ No newline at end of file
+endif
+
+# Helper function to check MG5 version
+define CHECK_MG5AMC_VERSION
+python -c 'import re; from distutils.version import StrictVersion; print StrictVersion("$(MG5AMC_VERSION)") >= StrictVersion("$(1)") if re.match("^[\d\.]+$$","$(MG5AMC_VERSION)") else True;'
+endef
\ No newline at end of file

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

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

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

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

=== modified file 'Template/LO/SubProcesses/setcuts.f'
--- Template/LO/SubProcesses/setcuts.f	2016-06-16 15:03:59 +0000
+++ Template/LO/SubProcesses/setcuts.f	2016-09-01 23:40:50 +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-08-22 13:52:25 +0000
+++ Template/LO/SubProcesses/unwgt.f	2016-09-01 23:40:50 +0000
@@ -240,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)
@@ -441,7 +441,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**************************************************************************
@@ -464,15 +464,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)
@@ -481,19 +482,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
@@ -514,9 +510,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
@@ -533,17 +526,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
@@ -742,6 +742,19 @@
          write(buffclus(nexternal),'(a)')'</clustering>'
       endif
 
+C     If the arguments of write_event have already been set in the
+C     bias module, then the beginning of the routine will directly
+C     jump here.
+
+ 1123 continue
+      if (.not.do_write_events) then
+        return
+      endif
+
+c     Store weight for event
+      nw = nw+1
+      swgt(nw)=wgt
+
       call write_event(lun,pb(0,1),wgt,npart,jpart(1,1),ngroup,
      &   sscale,aaqcd,aaqed,buff,use_syst,s_buff,nclus,buffclus)
       if(btest(mlevel,1))

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

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2016-08-26 13:55:42 +0000
+++ UpdateNotes.txt	2016-09-01 23:40:50 +0000
@@ -4,22 +4,22 @@
 2.5.0 (xx/xx/xx)
      FUNCTIONALITY
      -------------
+       VH+OM: Adding an official interface to PY8 for the parton-shower
        OM: Introduces a new function for LO/NLO interface "systematics"
             This function allows to compute systematics uncertainty from the event sample
             It requires the event sample to have been generated with 
                 - use_syst = T (for LO sample)
                 - store_reweight_info = T (for NLO sample)
             At LO the code is run automatically if use_syst=T (but if SysCalc is installed)
-
+       VH+SP: extend support for CKKWL
+       VH: extend install command to install: lhapdf/pythia8
        VH: adding the possibility to install COLLIER and to use it to reduce the loop-matrix element
        OM+VH: At the first loop/NLO computation, a new question will now be asked to choose which program
            to install to compute the loop. You can still install additional method later via the "install" command
-       OM: add an automatic update of the param_card to write the correct value for all dependent parameter.
-       OM: add the check that the param_card is compatible with the model restriction.
-       OM: NLO/LO Re-weighting works in multi-core
-
-     CODE IMPROVMENT
-     ---------------
+
+
+     CODE IMPROVMENT / small feature
+     -------------------------------
        OM: Modify the structure of the output format such that all the internal format have the same structure
        OM: Adding the Plugin directory. Three kind of plugin are currently supported
            - plugin defining a new type of output format
@@ -27,6 +27,12 @@
            - plugin modifying the main interface of MG5aMCnlo 
            More informations/examples are available here:
            https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Plugin
+       OM: Adding the possiblity of having detailled help at the time of the edition of the cards.
+            help mass / help mt / help nevents provided some information on the parameters.
+       OM: NLO/LO Re-weighting works in multi-core
+       OM: add an automatic update of the param_card to write the correct value for all dependent parameter.
+       OM: add the check that the param_card is compatible with the model restriction.
+       OM: Adding the run_card options "event_norm" for the LO run_card (same meaning as NLO one)
 
      BUG FIXING
      ----------
@@ -34,7 +40,6 @@
        OM: Fix a bug in the reweight_card where relative path was not understood from the local directory 
             where the program was runned by the user.
 
-
 2.4.3 (01/08/16)
         OM: Reduce the amount of log file/output generated for LO run (output can use up to three times less output).
         OM: For the LO combination of events (unweighting) pass to the method previously used for loop-induced.

=== modified file 'VERSION'
--- VERSION	2016-08-26 13:55:42 +0000
+++ VERSION	2016-09-01 23:40:50 +0000
@@ -1,5 +1,6 @@
 version = 2.5.0.alpha
-date = 2016-08-26
+date = 2016-08-28
+
 
 
 

=== modified file 'input/.mg5_configuration_default.txt'
--- input/.mg5_configuration_default.txt	2016-06-16 09:04:34 +0000
+++ input/.mg5_configuration_default.txt	2016-09-01 23:40:50 +0000
@@ -56,10 +56,17 @@
 # timeout = 60
 
 #! Pythia8 path.
-#!  Defines the path to the main pythia8 directory (i.e. that containing
-#!  the pythia8 configure script) .
+#!  Defines the path to the pythia8 installation directory (i.e. the
+#!  on containing the lib, bin and include directories) .
 #!  If using a relative path, that starts from the mg5 directory
-# pythia8_path =
+# pythia8_path = ./HEPTools/pythia8
+
+#! MG5aMC_PY8_interface path
+#!  Defines the path of the C++ driver file that is used by MG5_aMC to
+#!  steer the Pythia8 shower.
+#!  Can be installed directly from within MG5_aMC with the following command:
+#!     MG5_aMC> install mg5amc_py8_interface
+# mg5amc_py8_interface_path = ./HEPTools/MG5aMC_PY8_interface 
 
 #! Herwig++ paths
 #!  specify here the paths also to HepMC ant ThePEG
@@ -132,10 +139,14 @@
 #!  relative path start from main directory
 # delphes_path = ./Delphes
 
-#! MadAnalysis Package [For Drawing output]
+#! MadAnalysis4 fortran-based package [for basic analysis]
 #!  relative path start from main directory
 # madanalysis_path = ./MadAnalysis
 
+#! MadAnalysis5 python-based Package [For advanced analysis]
+#!  relative path start from main directory
+# madanalysis5_path = ./HEPTools/madanalysis5/madanalysis5
+
 #! ExRootAnalysis Package
 #!  relative path start from main directory
 # exrootanalysis_path = ./ExRootAnalysis

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2016-08-15 14:06:30 +0000
+++ madgraph/core/base_objects.py	2016-09-01 23:40:50 +0000
@@ -1270,8 +1270,6 @@
         self.name2part = {}
         for part in self.get("particle_dict").values():
             self.name2part[part.get('name')] = part
-        
-            
 
     def get_lorentz(self, name):
         """return the lorentz object from the associate name"""

=== modified file 'madgraph/interface/amcatnlo_interface.py'
--- madgraph/interface/amcatnlo_interface.py	2016-06-20 11:08:55 +0000
+++ madgraph/interface/amcatnlo_interface.py	2016-09-01 23:40:50 +0000
@@ -323,9 +323,7 @@
 
     def __init__(self, mgme_dir = '', *completekey, **stdin):
         """ Special init tasks for the Loop Interface """
-
         mg_interface.MadGraphCmd.__init__(self, mgme_dir = '', *completekey, **stdin)
-        misc.sprint(type(self.history))
         self.setup()
 
     def setup(self):
@@ -462,6 +460,8 @@
 
         self.proc_validity(myprocdef,'aMCatNLO_%s'%proc_type[1])
 
+        self._curr_proc_defs.append(myprocdef)
+
 #        if myprocdef['perturbation_couplings']!=['QCD']:
 #            message = ""FKS for reals only available in QCD for now, you asked %s" \
 #                        % ', '.join(myprocdef['perturbation_couplings'])"
@@ -519,6 +519,8 @@
         if self._export_format in ['NLO']:
             self._curr_exporter = export_v4.ExportV4Factory(self, noclean, 
                       output_type='amcatnlo',group_subprocesses=group_processes)
+            
+            self._curr_exporter.pass_information_from_cmd(self)
 
         # check if a dir with the same name already exists
         if not force and not noclean and os.path.isdir(self._export_dir)\
@@ -546,6 +548,9 @@
         # Perform export and finalize right away
         self.export(nojpeg, main_file_name, group_processes=group_processes)
 
+        # Pass potential new information generated during the export.
+        self._curr_exporter.pass_information_from_cmd(self)
+
         # Automatically run finalize
         self.finalize(nojpeg)
             
@@ -655,7 +660,9 @@
             global glob_directories_map
             glob_directories_map = []
 
+            # Save processes instances generated
             self.born_processes_for_olp = []
+            self.born_processes = []
             for ime, me in \
                 enumerate(self._curr_matrix_elements.get('matrix_elements')):
                 if not self.options['low_mem_multicore_nlo_generation']:
@@ -667,6 +674,7 @@
                             path,self.options['OLP'])
                     self._fks_directories.extend(self._curr_exporter.fksdirs)
                     self.born_processes_for_olp.append(me.born_matrix_element.get('processes')[0])
+                    self.born_processes.append(me.born_matrix_element.get('processes'))
                 else:
                     glob_directories_map.append(\
                             [self._curr_exporter, me, self._curr_helas_model, 

=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
--- madgraph/interface/amcatnlo_run_interface.py	2016-08-26 16:21:09 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2016-09-01 23:40:50 +0000
@@ -1222,6 +1222,8 @@
         if not mode in ['LO', 'NLO', 'noshower', 'noshowerLO'] \
                                                       and not options['parton']:
             self.run_mcatnlo(evt_file, options)
+            self.exec_cmd('madanalysis5_hadron --no_default', postcmd=False, printcmd=False)
+
         elif mode == 'noshower':
             logger.warning("""You have chosen not to run a parton shower. NLO events without showering are NOT physical.
 Please, shower the Les Houches events before using them for physics analyses.""")
@@ -2962,14 +2964,14 @@
 
         
         # libdl may be needded for pythia 82xx
-        if shower == 'PYTHIA8' and not \
-           os.path.exists(pjoin(self.options['pythia8_path'], 'xmldoc')) and \
-           'dl' not in self.shower_card['extralibs'].split():
-            # 'dl' has to be linked with the extralibs
-            self.shower_card['extralibs'] += ' dl'
-            logger.warning("'dl' was added to extralibs from the shower_card.dat.\n" + \
-                          "It is needed for the correct running of PY8.2xx.\n" + \
-                          "If this library cannot be found on your system, a crash will occur.")
+        #if shower == 'PYTHIA8' and not \
+        #   os.path.exists(pjoin(self.options['pythia8_path'], 'xmldoc')) and \
+        #   'dl' not in self.shower_card['extralibs'].split():
+        #    # 'dl' has to be linked with the extralibs
+        #    self.shower_card['extralibs'] += ' dl'
+        #    logger.warning("'dl' was added to extralibs from the shower_card.dat.\n" + \
+        #                  "It is needed for the correct running of PY8.2xx.\n" + \
+        #                  "If this library cannot be found on your system, a crash will occur.")
 
         misc.call(['./MCatNLO_MadFKS.inputs'], stdout=open(mcatnlo_log, 'w'),
                     stderr=open(mcatnlo_log, 'w'), 
@@ -3298,96 +3300,91 @@
 
     ############################################################################
     def set_run_name(self, name, tag=None, level='parton', reload_card=False):
-         """define the run name, the run_tag, the banner and the results."""
-         
-         # when are we force to change the tag new_run:previous run requiring changes
-         upgrade_tag = {'parton': ['parton','pythia','pgs','delphes','shower'],
-                        'pythia': ['pythia','pgs','delphes'],
-                        'shower': ['shower'],
-                        'pgs': ['pgs'],
-                        'delphes':['delphes'],
-                        'plot':[]}
-         
-         
- 
-         if name == self.run_name:        
-             if reload_card:
-                 run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
-                 self.run_card = banner_mod.RunCardNLO(run_card)
- 
-             #check if we need to change the tag
-             if tag:
-                 self.run_card['run_tag'] = tag
-                 self.run_tag = tag
-                 self.results.add_run(self.run_name, self.run_card)
-             else:
-                 for tag in upgrade_tag[level]:
-                     if getattr(self.results[self.run_name][-1], tag):
-                         tag = self.get_available_tag()
-                         self.run_card['run_tag'] = tag
-                         self.run_tag = tag
-                         self.results.add_run(self.run_name, self.run_card)                        
-                         break
-             return # Nothing to do anymore
-         
-         # save/clean previous run
-         if self.run_name:
-             self.store_result()
-         # store new name
-         self.run_name = name
-         
-         # Read run_card
-         run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
-         self.run_card = banner_mod.RunCardNLO(run_card)
- 
-         new_tag = False
-         # First call for this run -> set the banner
-         self.banner = banner_mod.recover_banner(self.results, level, self.run_name, tag)
-         if tag:
-             self.run_card['run_tag'] = tag
-             new_tag = True
-         elif not self.run_name in self.results and level =='parton':
-             pass # No results yet, so current tag is fine
-         elif not self.run_name in self.results:
-             #This is only for case when you want to trick the interface
-             logger.warning('Trying to run data on unknown run.')
-             self.results.add_run(name, self.run_card)
-             self.results.update('add run %s' % name, 'all', makehtml=True)
-         else:
-             for tag in upgrade_tag[level]:
-                 
-                 if getattr(self.results[self.run_name][-1], tag):
-                     # LEVEL is already define in the last tag -> need to switch tag
-                     tag = self.get_available_tag()
-                     self.run_card['run_tag'] = tag
-                     new_tag = True
-                     break
-             if not new_tag:
-                 # We can add the results to the current run
-                 tag = self.results[self.run_name][-1]['tag']
-                 self.run_card['run_tag'] = tag # ensure that run_tag is correct                
-              
-                     
-         if name in self.results and not new_tag:
-             self.results.def_current(self.run_name)
-         else:
-             self.results.add_run(self.run_name, self.run_card)
- 
-         self.run_tag = self.run_card['run_tag']
- 
-         # Return the tag of the previous run having the required data for this
-         # tag/run to working wel.
-         if level == 'parton':
-             return
-         elif level == 'pythia':
-             return self.results[self.run_name][0]['tag']
-         else:
-             for i in range(-1,-len(self.results[self.run_name])-1,-1):
-                 tagRun = self.results[self.run_name][i]
-                 if tagRun.pythia:
-                     return tagRun['tag']
- 
- 
+        """define the run name, the run_tag, the banner and the results."""
+        
+        # when are we force to change the tag new_run:previous run requiring changes
+        upgrade_tag = {'parton': ['parton','delphes','shower','madanalysis5_hadron'],
+                       'shower': ['shower','delphes','madanalysis5_hadron'],
+                       'delphes':['delphes'],
+                       'madanalysis5_hadron':['madanalysis5_hadron'],
+                       'plot':[]}
+
+        if name == self.run_name:        
+            if reload_card:
+                run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
+                self.run_card = banner_mod.RunCardNLO(run_card)
+
+            #check if we need to change the tag
+            if tag:
+                self.run_card['run_tag'] = tag
+                self.run_tag = tag
+                self.results.add_run(self.run_name, self.run_card)
+            else:
+                for tag in upgrade_tag[level]:
+                    if getattr(self.results[self.run_name][-1], tag):
+                        tag = self.get_available_tag()
+                        self.run_card['run_tag'] = tag
+                        self.run_tag = tag
+                        self.results.add_run(self.run_name, self.run_card)                        
+                        break
+            return # Nothing to do anymore
+        
+        # save/clean previous run
+        if self.run_name:
+            self.store_result()
+        # store new name
+        self.run_name = name
+        
+        # Read run_card
+        run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
+        self.run_card = banner_mod.RunCardNLO(run_card)
+
+        new_tag = False
+        # First call for this run -> set the banner
+        self.banner = banner_mod.recover_banner(self.results, level, self.run_name, tag)
+        if tag:
+            self.run_card['run_tag'] = tag
+            new_tag = True
+        elif not self.run_name in self.results and level =='parton':
+            pass # No results yet, so current tag is fine
+        elif not self.run_name in self.results:
+            #This is only for case when you want to trick the interface
+            logger.warning('Trying to run data on unknown run.')
+            self.results.add_run(name, self.run_card)
+            self.results.update('add run %s' % name, 'all', makehtml=True)
+        else:
+            for tag in upgrade_tag[level]:
+                
+                if getattr(self.results[self.run_name][-1], tag):
+                    # LEVEL is already define in the last tag -> need to switch tag
+                    tag = self.get_available_tag()
+                    self.run_card['run_tag'] = tag
+                    new_tag = True
+                    break
+            if not new_tag:
+                # We can add the results to the current run
+                tag = self.results[self.run_name][-1]['tag']
+                self.run_card['run_tag'] = tag # ensure that run_tag is correct                
+             
+                    
+        if name in self.results and not new_tag:
+            self.results.def_current(self.run_name)
+        else:
+            self.results.add_run(self.run_name, self.run_card)
+
+        self.run_tag = self.run_card['run_tag']
+
+        # Return the tag of the previous run having the required data for this
+        # tag/run to working wel.
+        if level == 'parton':
+            return
+        elif level == 'pythia':
+            return self.results[self.run_name][0]['tag']
+        else:
+            for i in range(-1,-len(self.results[self.run_name])-1,-1):
+                tagRun = self.results[self.run_name][i]
+                if tagRun.pythia:
+                    return tagRun['tag']
 
 
     def store_result(self):
@@ -4457,10 +4454,10 @@
             options['reweightonly'] = False
         
         
-        void = 'NOT INSTALLED'
-        switch_order = ['order', 'fixed_order', 'shower','madspin', 'reweight']
+        void = 'Not installed'
+        switch_order = ['order', 'fixed_order', 'shower','madspin', 'reweight','madanalysis5']
         switch_default = {'order': 'NLO', 'fixed_order': 'OFF', 'shower': void,
-                  'madspin': void,'reweight':'OFF'}
+                  'madspin': void,'reweight':'OFF','madanalysis5':void}
         if not switch:
             switch = switch_default
         else:
@@ -4472,18 +4469,26 @@
                                 'fixed_order': default_switch,
                                 'shower': default_switch,
                                 'madspin': default_switch,
-                                'reweight': default_switch}
+                                'reweight': default_switch,
+                                'madanalysis5':['OFF','HADRON']}
+
+        if not os.path.exists(pjoin(self.me_dir, 'Cards', 
+                                       'madanalysis5_hadron_card_default.dat')):
+            allowed_switch_value['madanalysis5']=[]
 
         description = {'order':  'Perturbative order of the calculation:',
                        'fixed_order': 'Fixed order (no event generation and no MC@[N]LO matching):',
                        'shower': 'Shower the generated events:',
                        'madspin': 'Decay particles with the MadSpin module:',
-                       'reweight': 'Add weights to the events based on changing model parameters:'}
+                       'reweight': 'Add weights to the events based on changing model parameters:',
+                       'madanalysis5':'Run MadAnalysis5 on the events generated:'}
 
         force_switch = {('shower', 'ON'): {'fixed_order': 'OFF'},
                        ('madspin', 'ON'): {'fixed_order':'OFF'},
                        ('reweight', 'ON'): {'fixed_order':'OFF'},
-                       ('fixed_order', 'ON'): {'shower': 'OFF', 'madspin': 'OFF', 'reweight':'OFF'}
+                       ('fixed_order', 'ON'): {'shower': 'OFF', 'madspin': 'OFF', 'reweight':'OFF','madanalysis5':'OFF'},
+                       ('madanalysis5','HADRON'): {'shower': 'ON','fixed_order':'OFF'},
+                       ('shower','OFF'): {'madanalysis5': 'OFF'},   
                        }
         special_values = ['LO', 'NLO', 'aMC@NLO', 'aMC@LO', 'noshower', 'noshowerLO']
 
@@ -4494,10 +4499,12 @@
             switch['shower'] = 'Not available for decay'
             switch['madspin'] = 'Not available for decay'
             switch['reweight'] = 'Not available for decay'
+            switch['madanalysis5'] = 'Not available for decay'
             allowed_switch_value['fixed_order'] = ['ON']