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

[Olivier Mattelaer]

Hi guys,

As you know, I plan to release this version as beta (then only accessible via launchpad) next week.

They are still a couple of merge in progress. It would be great to have them include.
The most important (which are not yet approve) are:

lp:~maddevelopers/mg5amcnlo/sm_width
lp:~maddevelopers/mg5amcnlo/2.3.0_nopdftransfer (extremelly important for our LLN cluster. Marco any news on the review?)
lp:~maddevelopers/mg5amcnlo/LoopInduced_splitted_refine

And obviously we (i.e. Hua-Sheng) need to fix the huge crash linked to IREGI. Since without that we can not run anything in practise.

Do you see any point missing? Do you want to have additional stuff in that beta version?

On my side, I might try to implement a "bridge" version of MadSpin (Since MS can not handle loop induce) but we will see if I can do it in time.

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/2.3]

279. By Olivier Mattelaer on 2015-03-06

Fix a bug in presence of large number of resonances/conflict of BW

280. By Paolo Torrielli on 2015-03-06

merged with 2.2.0_shower_set

281. By Rikkert Frederix on 2015-03-06

merge with veto_xsec

282. By Rikkert Frederix on 2015-03-06

converted some parameters in the run_card from integer to float.

283. By Olivier Mattelaer on 2015-03-06

better fix of the acceptance test (all file are in temporary directory

284. By Olivier Mattelaer on 2015-03-07

Fix an error in the if statement

285. By Olivier Mattelaer on 2015-03-07

fix a stupid bug

286. By Olivier Mattelaer on 2015-03-09

Fixing a problem in the formatting of string (not preserving case) when using the set command

[Rikkert Frederix (frederix)]

Hi Olivier,

I don't think it makes sense to rush this too much and release a beta version. We can wait 2-3 weeks more to have a slightly more stable version out, I think. In any case, I would like to have the HwU_shower branch to be part of it as well (it's currently under review).

Best,
Rikkert

[lp:~maddevelopers/mg5amcnlo/2.3]

287. By Olivier Mattelaer on 2015-03-09

fix a bug in the addmodel function introduces in 2.2.3

288. By Valentin Hirschi on 2015-03-10

1. Merge the fix for the I/O crash in local multicore when reading MadLoop5_resources/HelFilter.dat

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rik,

My feeling is that the version is stable enough for a beta.
I have two reason to rush a bit this version
1) We have a lot of pressure to release the loop-induced code
2) This version is going to include the no_pdf_transfer utilities which is needed to reduce the pressure of our code on the LLN cluster. I promise the LLN IT team to have that
option available as soon as possible. They (=Jerome) have work a lot to have that functionality available on the ingrid cluster in the hope that our tool will stop (or at least decrease)
to slow down the full filesystem.

Now I’m fine to postpone it if this can help.

> In any case, I would like to have the HwU_shower branch to be part of it as well (it's currently under review).

looks like well in progress, I will also review that branch later today (china time).

Cheers,

Olivier

On 09 Mar 2015, at 21:35, Rikkert Frederix <email address hidden> wrote:

> Hi Olivier,
>
> I don't think it makes sense to rush this too much and release a beta version. We can wait 2-3 weeks more to have a slightly more stable version out, I think. In any case, I would like to have the HwU_shower branch to be part of it as well (it's currently under review).
>
> Best,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3/+merge/252090
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/2.3.

[lp:~maddevelopers/mg5amcnlo/2.3]

289. By marco zaro on 2015-03-10

two small fixes in export_fks relevant to BSM cases when only diagrams with
>=4 particle-vertices exist

290. By Valentin Hirschi on 2015-03-10

1. Removed the double-checking of the vanishing helicity configurations for the actual runs. This solves the issue of the multiple write access from different jobs of a given channel when running in local multicore mode without cluster_tmp_path.

291. By Valentin Hirschi on 2015-03-10

1. Small fix that prevented 'check stability' to work properly.

292. By Rikkert Frederix on 2015-03-11

merge with the HwU_shower branch

293. By Rikkert Frederix on 2015-03-11

Updated the UpdateNotes

294. By Rikkert Frederix on 2015-03-11

some trivial fixes in NLO/Source/setrun

295. By Olivier Mattelaer on 2015-03-11

allow faster identification of jobs which are pointless to compute during survey time. Improve error message with MultiCore class

296. By Olivier Mattelaer on 2015-03-11

allow survey.sh/refine_splitted.sh to remove log output to reduce the charge on the cluster.

297. By marco zaro on 2015-03-12

merged with 2.2.3-LOfromNLO

298. By Olivier Mattelaer on 2015-03-12

improve logging information in LoopInduced

299. By Valentin Hirschi on 2015-03-13

1. Fix in template a bug introduced in the merge of LOfromNLO.

300. By Olivier Mattelaer on 2015-03-13

fix bug in LoopInduced if some channel are in 0result

301. By Olivier Mattelaer on 2015-03-13

Fix problem with the re-weighting method

302. By Paolo Torrielli on 2015-03-13

Bug fix for Pythia8 in conjunction with LHAPDF6.
Thanks to Federico Demartin for spotting this bug.

303. By Paolo Torrielli on 2015-03-13

Bug fix for Herwig++ in conjunction with LHAPDF6

304. By Olivier Mattelaer on 2015-03-14

allow -l option to take digit input for more refine level of print

305. By Valentin Hirschi on 2015-03-17

1. merged HwUNewFunction

306. By Rikkert Frederix on 2015-03-17

Added support for gnuplot v5 in the HwU histograms.py.
Also, added a string to the gnuplot output that contains the command
that was executed to get this output (standalone mode only). Useful to
keep track which curve is which when many are combined.

307. By Valentin Hirschi on 2015-03-17

1. Removed unnecessary print statements in histograms.py

308. By Valentin Hirschi on 2015-03-17

1. Small modification of the plotting of multiple .HwU file so that the
   '#i' prefixes given to the histogram types do not stack when re-plotting
   an output HwU.

[Olivier Mattelaer (olivier-mattelaer)]

Ok I will release this in a beta during the weekend.
This version was lating since too long.

Likely that the bridge review will not be ontime, but it can wait.
On the otherhand the 2.3.0_nopdftransfer is very important and should be merge in any case.

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/2.3]

309. By Valentin Hirschi on 2015-03-19

1. Merged with branch for the new Madloop monitoring for loop-induced simulation.

310. By Rikkert Frederix on 2015-03-19

Fixed a bug in the Pythia8 driver files: the call to pythia.init()
should be AFTER the UserHooks (in particular the JetMatching for FxFx)
has been setup. This bug was introduced in v.2.2.3 and gave the wrong
results without any error message.

311. By Rikkert Frederix on 2015-03-19

also updated the pythia8.1 driver files (see previous commit)

312. By Olivier Mattelaer on 2015-03-19

improve the print_result command for easy parsing option

313. By Olivier Mattelaer on 2015-03-19

merge the no pdf transfer branch

314. By Olivier Mattelaer on 2015-03-20

nicer implementation of the check of unweight event in LoopInduce mode and nicer printout

315. By Olivier Mattelaer on 2015-03-20

Fix a bug in the reweighting interface

316. By Valentin Hirschi on 2015-03-20

1. Small esthetical change in histograms.py

317. By Paolo Torrielli on 2015-03-20

added some HwU shower analyses

318. By Olivier Mattelaer on 2015-03-20

merge with the bridge branch --with bridge mode deactivate-- but with none mode + fix a cluster problem

319. By Olivier Mattelaer on 2015-03-20

remove a debug statement

320. By Olivier Mattelaer on 2015-03-23

Fix a problem in the matching/merging leading to wrong QCD/QED assignment (leading to crash)

321. By Olivier Mattelaer on 2015-03-23

Fix a problem of MadSpin in gridpack mode for the mssm model

322. By Paolo Torrielli on 2015-03-23

other HwU shower analyses

323. By marco zaro on 2015-03-24

 fix in born_fks.inc: bug introduced with the LI branch related to scalar octets

324. By Olivier Mattelaer on 2015-03-25

allow Auto in the reweight interface/fix small bug in recreating html db/ add new function for the FourMomenta class

325. By Olivier Mattelaer on 2015-03-25

correct a problem on condor for submitting a second run in the same directory

326. By marco zaro on 2015-03-25

the logs are collected in a more memory-efficient way

327. By Rikkert Frederix on 2015-03-25

small improvements for the NNLL+NLO veto cross section runs.

328. By Marco Zaro <mzaro@ingrid-ui1> on 2015-03-26

crucial fix in copy_lhapdf_set

329. By Paolo Torrielli on 2015-03-30

Rik's fix to the bug in the MG5_aMC + PY8 interface, concerning
the theoretical-uncertainty bands at (N)LO+PS, as plotted by
topdrawer analyses.

330. By Olivier Mattelaer on 2015-03-31

Fix a tricky bug in the cmd interface when multiple file where simultaneously used as script for various part of the code (i.e. typical use of MadSpin/reweight)

331. By Olivier Mattelaer on 2015-03-31

fix the problem of missing initialisation of the text_editor

332. By Olivier Mattelaer on 2015-03-31

Fix problem in single core mode

333. By Paolo Torrielli on 2015-04-02

changed the default PDF set for PY8, from CTEQ6L to NNPDF 2.3 QCD+QED LO

334. By Olivier Mattelaer on 2015-04-05

merge the fix of MadSpin done by Pierre

335. By Olivier Mattelaer on 2015-04-05

fix a problem with the <init> block

336. By Valentin Hirschi on 2015-04-08

1. Merged the branch implementing loop-induced simulation features in MadGraph5_aMCatNLO using the syntax [noborn=<pert_orders>] in the process definition.

[Valentin Hirschi (valentin-hirschi)]

Hi all,

We have finalized our changes in the loop-induced branch, merged it and fixed the tests.
This revision 2.3 is now to be considered frozen and we will release the code tomorrow.
If you have any opinion against this, please let us know now.

Olivier & Valentin

[lp:~maddevelopers/mg5amcnlo/2.3]

337. By Valentin Hirschi on 2015-04-08

1. Fixed a diagram_generation test which has incorrectly fixed before.

338. By Valentin Hirschi on 2015-04-08

1. Fixed MadLoop compilation command in the parallel tests.
2. Reverted two personal modification for testing purposes which were left in the code.

339. By Valentin Hirschi on 2015-04-08

1. Fixed a minor bug in the helicity filter double-checking procedure of MadLoop.
2. Fixed the way the loop parallel test adjust MadLoopParams.dat by now using the new python banner handling system.

340. By Valentin Hirschi on 2015-04-08

1. Improved the debug printouts on the MadLoop Initialization stage.

341. By Valentin Hirschi on 2015-04-08

1. Updated IOTest with the minor fix on MadLoop helicity doublecheck procedure.

342. By Valentin Hirschi on 2015-04-08

1. Small fix brought to the ML5EW steering file as well.

343. By marco zaro on 2015-04-08

typos fixed in .mg5_configuration_default.txt

344. By Olivier Mattelaer on 2015-04-08

Update Version/UpdateNotes

345. By Olivier Mattelaer on 2015-04-09

Fixing two of the parralel tests. One is still failing

346. By Olivier Mattelaer on 2015-04-09

Fix the parralel test which was crashing for interference computation

347. By Valentin Hirschi on 2015-04-09

1. Improved the loop instability security to also apply on the survey, but not its first iteration, and to automatically not apply for channels which seem to be analytically zero (i.e. with a very small relative across contribution).
2. Fixed an issue with the sqrts definition in check_sa.f for loop induced processes.

348. By Valentin Hirschi on 2015-04-10

1. Synchronized some IOTests.

349. By Olivier Mattelaer on 2015-04-10

Fix a problem of multiple run with pythia on for loop-induced where all samples were identical

350. By Valentin Hirschi on 2015-04-10

1. Fixed a comment.
2. Removed another hardcoded path in the MODEL makefile.

351. By Valentin Hirschi on 2015-04-10

1. Updated the parrallel test as a result of having doubled the center of mass energy in certain cases.

[Olivier Mattelaer (olivier-mattelaer)]

Ok Version released as Golden Master.

So please do not push any modication to this branch.

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/2.3]

352. By Valentin Hirschi on 2015-04-10

1. Fixed a bug in histograms.py pointed out by Rik.
2. Vetoed from the multichanneling loop-induced diagram above boxes.

353. By Valentin Hirschi on 2015-04-10

1. Updated version and notes date.

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2015-02-03 15:18:06 +0000
+++ MadSpin/decay.py	2015-04-10 09:55:25 +0000
@@ -63,7 +63,7 @@
 from madgraph import MG5DIR, MadGraph5Error
 import madgraph.various.misc as misc
 #import time
-import tests.unit_tests.various.test_aloha as test_aloha
+import tests.parallel_tests.test_aloha as test_aloha
 
 class MadSpinError(MadGraph5Error):
     pass
@@ -441,6 +441,7 @@
             self['tree'][propa_id] = {'nbody': len(proc.get('legs'))-1,\
                                       'label':proc.get('legs')[0].get('id')}
             # loop over the child
+            child_propa_id = propa_id
             for c_nb,leg in enumerate(proc.get('legs')):
                 if c_nb == 0:
                     continue
@@ -449,12 +450,14 @@
                 self["tree"][propa_id]["d%s" % c_nb]["label"] = c_pid
                 self["tree"][propa_id]["d%s" % c_nb]["labels"] = [c_pid]
                 if c_pid in to_decay:
-                    self["tree"][propa_id]["d%s" % c_nb]["index"] = propa_id-1
+                    child_propa_id -= 1
+                    self["tree"][propa_id]["d%s" % c_nb]["index"] = child_propa_id
                     self.nb_decays += 1
-                    add_decay(to_decay[c_pid].pop(), propa_id-1)
+                    child_propa_id = add_decay(to_decay[c_pid].pop(), child_propa_id)
                 else:
                     self.nexternal += 1
                     self["tree"][propa_id]["d%s" % c_nb]["index"] = self.nexternal
+            return child_propa_id
         
         # launch the recursive loop
         add_decay(process)
@@ -1063,9 +1066,11 @@
 
         s_and_t_channels = []
 
-        minvert = min([max([d for d in config if d][0].get_vertex_leg_numbers()) \
-                                             for config in configs])
-
+        vert_list = [max([d for d in config if d][0].get_vertex_leg_numbers()) \
+          for config in configs if [d for d in config if d][0].\
+                                             get_vertex_leg_numbers()!=[]]
+        minvert = min(vert_list) if vert_list!=[] else 0
+    
         # Number of subprocesses
         #    nsubprocs = len(configs[0])
 
@@ -1075,7 +1080,7 @@
 
         for iconfig, helas_diags in enumerate(configs):
             if any([vert > minvert for vert in
-                            [d for d in helas_diags if d][0].get_vertex_leg_numbers()]):
+                    [d for d in helas_diags if d][0].get_vertex_leg_numbers()]):
                     # Only 3-vertices allowed in configs.inc
                     continue
             nconfigs += 1
@@ -1315,12 +1320,12 @@
 
         # Maybe the branching fractions are already given in the banner:
         self.extract_br_from_banner(self.banner)
-        to_decay = [p for p in to_decay if not p in self.br]
+        to_decay = list(to_decay)
         for part in to_decay[:]:
             if part in mgcmd._multiparticles:
                 to_decay += [self.pid2label[id] for id in mgcmd._multiparticles[part]]
                 to_decay.remove(part)
-        to_decay = list(set(to_decay))
+        to_decay = list(set([p for p in to_decay if not p in self.br]))
         
         if to_decay:
             logger.info('We need to recalculate the branching fractions for %s' % ','.join(to_decay))
@@ -2855,6 +2860,7 @@
         #to remove potential pointless decay in the diagram generation.
         resonances = decay_misc.get_all_resonances(self.banner, 
                          self.mgcmd, self.mscmd.list_branches.keys())
+
         logger.debug('List of resonances:%s' % resonances)
         path_me = os.path.realpath(self.path_me) 
         width = width_estimate(resonances, path_me, self.banner, self.model,
@@ -3526,7 +3532,8 @@
                 else:
                     # now we need to write the decay products in the event
                     # follow the decay chain order, so that we can easily keep track of the mother index
-                        
+                       
+                    map_to_part_number={}
                     for res in range(-1,-len(decay_struct[part]["tree"].keys())-1,-1):
                         index_res_for_mom=decay_struct[part]['mg_tree'][-res-1][0]
                         if (res==-1):
@@ -3548,8 +3555,8 @@
                                 "mass":mass,"helicity":helicity}
                             decayed_event.event2mg[part_number]=part_number
 
-                        mothup1=part_number
-                        mothup2=part_number
+                            map_to_part_number[res]=part_number
+   
 #
 #             Extract color information so that we can write the color flow
 #
@@ -3654,7 +3661,10 @@
                             decay_struct[part]["tree"][indexd1]["colup2"]=d1colup2
                             istup=2                    
                             mass=mom.m
-                        
+                            map_to_part_number[indexd1]=part_number 
+ 
+                        mothup1=map_to_part_number[res]
+                        mothup2=map_to_part_number[res]
                         decayed_event.particle[part_number]={"pid":pid,\
                                 "istup":istup,"mothup1":mothup1,"mothup2":mothup2,\
                                 "colup1":d1colup1,"colup2":d1colup2,"momentum":mom,\
@@ -3685,10 +3695,10 @@
                             decay_struct[part]["tree"][indexd2]["colup1"]=d2colup1
                             decay_struct[part]["tree"][indexd2]["colup2"]=d2colup2
                             mass=mom.m
-
-
-                        mothup1=part_number-2
-                        mothup2=part_number-2
+                            map_to_part_number[indexd2]=part_number
+
+                        mothup1=map_to_part_number[res]
+                        mothup2=map_to_part_number[res]
                         decayed_event.particle[part_number]={"pid":pid,"istup":istup,\
                            "mothup1":mothup1,"mothup2":mothup2,"colup1":d2colup1,\
                            "colup2":d2colup2,\
@@ -3776,8 +3786,7 @@
                 assert production['total_br'] - min(total_br) < 1e-4
         
         self.branching_ratio = max(total_br) * eff
-        
-        
+
         #self.banner['madspin'] += ms_banner
         # Update cross-section in the banner
         if 'mggenerationinfo' in self.banner:
@@ -3785,9 +3794,15 @@
             for i,line in enumerate(mg_info):
                 if 'Events' in line:
                     if eff == 1:
+                        self.err_branching_ratio = 0
                         continue
-                    nb_event =  int(int(mg_info[i].split()[-1]) * eff) 
+                    initial_event = int(mg_info[i].split()[-1])
+                    nb_event =  int(initial_event * eff) 
                     mg_info[i] = '#  Number of Events        :       %i' % nb_event
+                    if eff >0.5:
+                        self.err_branching_ratio = max(total_br) * math.sqrt(initial_event - eff * initial_event)/initial_event
+                    else:
+                        self.err_branching_ratio = max(total_br) * math.sqrt(eff * initial_event)/initial_event
                     continue
                 if ':' not in line:
                     continue
@@ -3801,6 +3816,9 @@
                 else:
                     mg_info[i] = '%s : %s' % (info, value * self.branching_ratio)
                 self.banner['mggenerationinfo'] = '\n'.join(mg_info)
+                
+            
+        
         if 'init' in self.banner:
             new_init =''
             for line in self.banner['init'].split('\n'):

=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2015-02-10 02:53:57 +0000
+++ MadSpin/interface_madspin.py	2015-04-10 09:55:25 +0000
@@ -14,9 +14,11 @@
 ################################################################################
 """ Command interface for MadSpin """
 from __future__ import division
+import collections
 import logging
 import math
 import os
+import random
 import re
 import shutil
 import glob
@@ -33,11 +35,13 @@
 import madgraph.various.misc as misc
 import madgraph.iolibs.files as files
 import madgraph.various.banner as banner
+import madgraph.various.lhe_parser as lhe_parser
 
 import models.import_ufo as import_ufo
 import models.check_param_card as check_param_card
 import MadSpin.decay as madspin
 
+
 logger = logging.getLogger('decay.stdout') # -> stdout
 logger_stderr = logging.getLogger('decay.stderr') # ->stderr
 cmd_logger = logging.getLogger('cmdprint2') # -> print
@@ -64,6 +68,8 @@
         
         self.decay = madspin.decay_misc()
         self.model = None
+        self.mode = "madspin" # can be flat/bridge change the way the decay is done.
+                              # note amc@nlo does not support bridge.
         
         self.options = {'max_weight': -1, 
                         'curr_dir': os.path.realpath(os.getcwd()),
@@ -73,7 +79,8 @@
                         'nb_sigma':0,
                         'ms_dir':None,
                         'max_running_process':100,
-                        'onlyhelicity': False}
+                        'onlyhelicity': False,
+                        'spinmode': "madspin"}
         
 
         
@@ -83,6 +90,7 @@
         self.to_decay={}
         self.mg5cmd = master_interface.MasterCmd()
         self.seed = None
+        self.err_branching_ratio = 0
         
         
         if event_path:
@@ -293,9 +301,17 @@
         """checking the validity of the set command"""
         
         if len(args) < 2:
-            raise self.InvalidCmd('set command requires at least two argument.')
-        
-        valid = ['max_weight','seed','curr_dir']
+            if args and '=' in args[0]:
+                name, value = args[0].split('=')
+                args[0]= name
+                args.append(value)
+            else:
+                raise self.InvalidCmd('set command requires at least two argument.')
+        
+        if args[1].strip() == '=':
+            args.pop(1)
+        
+        valid = ['max_weight','seed','curr_dir', 'spinmode']
         if args[0] not in self.options and args[0] not in valid:
             raise self.InvalidCmd('Unknown options %s' % args[0])        
     
@@ -316,6 +332,11 @@
         elif args[0] == 'curr_dir':
             if not os.path.isdir(args[1]):
                 raise self.InvalidCmd('second argument should be a path to a existing directory')
+        
+        elif args[0] == "spinmode":
+            if args[1].lower() not in ["full", "bridge", "none"]:
+                raise self.InvalidCmd("spinmode can only take one of those 3 value: full/bridge/none") 
+        
 
         
     def do_set(self, line):
@@ -324,7 +345,7 @@
         args = self.split_arg(line)
         self.check_set(args)
         
-        if args[0] in  ['max_weight', 'BW_effect','ms_dir']:
+        if args[0] in  ['max_weight', 'BW_effect','ms_dir', 'spinmode']:
             self.options[args[0]] = args[1]
             if args[0] == 'ms_dir':
                 self.options['curr_dir'] = self.options['ms_dir']
@@ -346,7 +367,7 @@
 
         # Format
         if len(args) == 1:
-            opts = self.options.keys() + ['seed']
+            opts = self.options.keys() + ['seed', "spinmode"]
             return self.list_completion(text, opts) 
         elif len(args) == 2:
             if args[1] == 'BW_effect':
@@ -357,7 +378,9 @@
             curr_path = pjoin(*[a for a in args \
                                                    if a.endswith(os.path.sep)])
             return self.path_completion(text, curr_path, only_dirs = True)
-        
+        elif args[1] == "spinmode":
+            return self.list_completion(text, ["full", "none"], line)
+         
     def help_set(self):
         """help the set command"""
         
@@ -373,6 +396,8 @@
         print '       It has a period of 2**19937-1.'
         print '   - set max_running_process VALUE: allow to limit the number of open file used by the code'
         print '       The number of running is raising like 2*VALUE'
+        print '   - set spinmode=none: mode with simple file merging. No spin correlation attempt.'
+        print '       This mode allows 3 (and more) body decay.'
     
     def do_define(self, line):
         """ """
@@ -419,6 +444,11 @@
     def do_launch(self, line):
         """end of the configuration launched the code"""
         
+        if self.options["spinmode"] in ["none"]:
+            return self.run_bridge(line)
+        elif self.options["spinmode"] == "bridge":
+            raise Exception, "Bridge mode not yet available due to lack of validation."
+        
         if self.options['ms_dir'] and os.path.exists(pjoin(self.options['ms_dir'], 'madspin.pkl')):
             return self.run_from_pickle()
         
@@ -463,6 +493,11 @@
         generate_all.run()
                         
         self.branching_ratio = generate_all.branching_ratio
+        try:
+            self.err_branching_ratio = generate_all.err_branching_ratio
+        except Exception:
+            self.err_branching_ratio = 0
+            
         evt_path = self.events_file.name
         try:
             self.events_file.close()
@@ -489,7 +524,7 @@
             ms_card_to_copy = pjoin(base_path,'madspin_card_for_%s.dat'%evt_name)
             count = 0    
             while os.path.exists(ms_card_to_copy):
-                count +=1
+                count += 1
                 ms_card_to_copy = pjoin(base_path,'madspin_card_for_%s_%d.dat'%\
                                                                (evt_name,count))
             files.cp(str(ms_card_path),str(ms_card_to_copy))
@@ -520,11 +555,19 @@
         
         if not hasattr(self.banner, 'param_card'):
             self.banner.charge_card('slha')
+        
+        # Special treatment for the mssm. Convert the param_card to the correct
+        # format
+        if self.banner.get('model').startswith('mssm-') or self.banner.get('model')=='mssm':
+            self.banner.param_card = check_param_card.convert_to_mg5card(\
+                    self.banner.param_card, writting=False)
+            
         for name, block in self.banner.param_card.items():
             if name.startswith('decay'):
                 continue
+                        
             orig_block = generate_all.banner.param_card[name]
-            if block != orig_block:
+            if block != orig_block:                
                 raise Exception, """The directory %s is specific to a mass spectrum. 
                 Your event file is not compatible with this one. (Different param_card: %s different)
                 orig block:
@@ -547,6 +590,7 @@
         
         generate_all.ending_run()
         self.branching_ratio = generate_all.branching_ratio
+        self.err_branching_ratio = generate_all.err_branching_ratio 
         evt_path = self.events_file.name
         try:
             self.events_file.close()
@@ -559,6 +603,266 @@
         if not self.mother:
             logger.info("Decayed events have been written in %s.gz" % decayed_evt_file)    
     
+    def run_bridge(self, line):
+        """Run the Bridge Algorithm"""
+        
+        # 1. Read the event file to check which decay to perform and the number
+        #   of event to generate for each type of particle.
+        # 2. Generate the events requested
+        # 3. perform the merge of the events.
+        #    if not enough events. re-generate the missing one.
+                
+        # First define an utility function for generating events when needed
+        def generate_events(pdg, nb_event, mg5, restrict_file=None, cumul=False):
+            """generate new events for this particle
+               restrict_file allow to only generate a subset of the definition
+               cumul allow to merge all the definition in one run (add process)
+                     to generate events according to cross-section
+            """
+            
+            part = self.model.get_particle(pdg)
+            name = part.get_name()
+            out = {}
+            logger.info("generate %s decay event for particle %s" % (nb_event, name))
+            if name not in self.list_branches:
+                return out
+            for i,proc in enumerate(self.list_branches[name]):
+                if restrict_file and i not in restrict_file:
+                    continue
+                decay_dir = pjoin(self.path_me, "decay_%s_%s" %(str(pdg).replace("-","x"),i))
+                if not os.path.exists(decay_dir):
+                    if cumul:
+                        mg5.exec_cmd("generate %s" % proc)
+                        for j,proc2 in enumerate(self.list_branches[name][1:]):
+                            if restrict_file and j not in restrict_file:
+                                raise Exception # Do not see how this can happen
+                            mg5.exec_cmd("add process %s" % proc2)
+                        mg5.exec_cmd("output %s -f" % decay_dir)
+                    else:
+                        mg5.exec_cmd("generate %s" % proc)
+                        mg5.exec_cmd("output %s -f" % decay_dir)
+                options = dict(mg5.options)
+                import madgraph.interface.madevent_interface as madevent_interface
+                me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
+                                                decay_dir), options=options)
+                me5_cmd.options["automatic_html_opening"] = False
+                run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))
+                run_card["nevents"] = int(1.2*nb_event)
+                run_card["iseed"] = self.seed
+                run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
+                self.seed += 1
+                me5_cmd.exec_cmd("generate_events run_01 -f")
+                me5_cmd.exec_cmd("exit")     
+                out[i] = lhe_parser.EventFile(pjoin(decay_dir, "Events", 'run_01', 'unweighted_events.lhe.gz'))            
+                if cumul:
+                    break
+            
+            return out
+
+        
+        args = self.split_arg(line)
+
+        #0. Define the path where to write the file
+        self.path_me = os.path.realpath(self.options['curr_dir']) 
+        if self.options['ms_dir']:
+            self.path_me = os.path.realpath(self.options['ms_dir'])
+            if not os.path.exists(self.path_me):
+                os.mkdir(self.path_me) 
+        else:
+            # cleaning
+            for name in glob.glob(pjoin(self.path_me, "decay_*_*")):
+                shutil.rmtree(name)
+
+        self.events_file.close()
+        orig_lhe = lhe_parser.EventFile(self.events_file.name)
+        
+        to_decay = collections.defaultdict(int)
+        nb_event = 0
+        for event in orig_lhe:
+            nb_event +=1
+            for particle in event:
+                if particle.status == 1 and particle.pdg in self.final_state:
+                    # final state and tag as to decay
+                    to_decay[particle.pdg] += 1
+
+        # Handle the banner of the output file
+        if not self.seed:
+            import random
+            self.seed = random.randint(0, int(30081*30081))
+            self.do_set('seed %s' % self.seed)
+            logger.info('Will use seed %s' % self.seed)
+            self.history.insert(0, 'set seed %s' % self.seed)
+
+        if self.seed > 30081*30081: # can't use too big random number
+            msg = 'Random seed too large ' + str(self.seed) + ' > 30081*30081'
+            raise Exception, msg
+
+        self.options['seed'] = self.seed
+        
+        text = '%s\n' % '\n'.join([ line for line in self.history if line])
+        self.banner.add_text('madspin' , text)
+
+
+        # 2. Generate the events requested
+        with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
+            mg5 = self.mg5cmd
+            modelpath = self.model.get('modelpath')
+            mg5.exec_cmd("import model %s" % modelpath)      
+            to_event = {}
+            for pdg, nb_needed in to_decay.items():
+                #check if a splitting is needed
+                if nb_needed == nb_event:
+                    to_event[pdg] = generate_events(pdg, nb_needed, mg5)
+                elif nb_needed %  nb_event == 0:
+                    nb_mult = nb_needed // nb_event
+                    part = self.model.get_particle(pdg)
+                    name = part.get_name()
+                    if name not in self.list_branches:
+                        continue
+                    elif len(self.list_branches[name]) == nb_mult:
+                        to_event[pdg] = generate_events(pdg, nb_event, mg5)
+                    else:
+                        to_event[pdg] = generate_events(pdg, nb_needed, mg5, cumul=True)
+                else:
+                    part = self.model.get_particle(pdg)
+                    name = part.get_name()
+                    if name not in self.list_branches or len(self.list_branches[name]) == 0:
+                        continue
+                    raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")
+                    
+                     
+                
+        
+        # Compute the branching ratio.
+        br = 1
+        for (pdg, event_files) in to_event.items():
+            if not event_files:
+                continue
+            totwidth = float(self.banner.get('param', 'decay', abs(pdg)).value)
+            if to_decay[pdg] == nb_event:
+                # Exactly one particle of this type to decay by event
+                pwidth = sum([event_files[k].cross for k in event_files])
+                if pwidth > 1.01 * totwidth:
+                    logger.critical("Branching ratio larger than one for %s " % pdg) 
+                br *= pwidth / totwidth
+            elif to_decay[pdg] % nb_event == 0:
+                # More than one particle of this type to decay by event
+                # Need to check the number of event file to check if we have to 
+                # make separate type of decay or not.
+                nb_mult = to_decay[pdg] // nb_event
+                if nb_mult == len(event_files):
+                    for k in event_files:
+                        pwidth = event_files[k].cross
+                        if pwidth > 1.01 * totwidth:
+                            logger.critical("Branching ratio larger than one for %s " % pdg)                       
+                        br *= pwidth / totwidth
+                    br *= math.factorial(nb_mult)
+                else:
+                    pwidth = sum(event_files[k].cross for k in event_files)
+                    if pwidth > 1.01 * totwidth:
+                        logger.critical("Branching ratio larger than one for %s " % pdg) 
+                    br *= (pwidth / totwidth)**nb_mult
+            else:
+                raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")
+        self.branching_ratio = br
+        # modify the cross-section in the init block of the banner
+        self.banner.scale_init_cross(self.branching_ratio)
+                    
+        
+        # 3. Merge the various file together.
+        output_lhe = lhe_parser.EventFile(orig_lhe.name.replace('.lhe', '_decayed.lhe.gz'), 'w')
+        self.banner.write(output_lhe, close_tag=False)
+        
+        # initialise object which store not use event due to wrong helicity
+        bufferedEvents_decay = {}
+        for pdg in to_event:
+            bufferedEvents_decay[pdg] = [{}] * len(to_event[pdg])
+        
+        import time
+        start = time.time()
+        counter = 0
+        orig_lhe.seek(0)
+        for event in orig_lhe:
+            if counter and counter % 1000 == 0 and float(str(counter)[1:]) ==0:
+                print "decaying event number %s [%s s]" % (counter, time.time()-start)
+            counter +=1
+            
+            # use random order for particles to avoid systematics when more than 
+            # one type of decay is asked.
+            particles = [p for p in event if int(p.status) == 1.0]
+            random.shuffle(particles)
+            ids = [particle.pid for particle in particles]
+            for i,particle in enumerate(particles):
+                # check if we need to decay the particle 
+                if particle.pdg not in self.final_state or particle.pdg not in to_event:
+                    continue # nothing to do for this particle
+                # check how the decay need to be done
+                nb_decay = len(to_event[particle.pdg])
+                if nb_decay == 0:
+                    continue #nothing to do for this particle
+                if nb_decay == 1:
+                    decay_file = to_event[particle.pdg][0]
+                    decay_file_nb = 0
+                elif ids.count(particle.pdg) == nb_decay:
+                    decay_file = to_event[particle.pdg][ids[:i].count(particle.pdg)]
+                    decay_file_nb = ids[:i].count(particle.pdg)
+                else:
+                    #need to select the file according to the associate cross-section
+                    r = random.random()
+                    tot = sum(events.cross for events in to_event[particle.pdg])
+                    r = r * tot
+                    cumul = 0
+                    for j,events in enumerate(to_event[particle.pdg]):
+                        cumul += events.cross
+                        if r < cumul:
+                            decay_file = events
+                            decay_file_nb = j
+                        else:
+                            break
+                        
+                # ok start the procedure
+                helicity = particle.helicity
+                bufferedEvents = bufferedEvents_decay[particle.pdg][decay_file_nb]
+                
+                # now that we have the file to read. find the associate event
+                # checks if we have one event in memory
+                if helicity in bufferedEvents and bufferedEvents[helicity]:
+                    decay = bufferedEvents[helicity].pop()
+                else:
+                    # read the event file up to completion
+                    while 1:
+                        try:
+                            decay = decay_file.next()
+                        except StopIteration:
+                            # check how far we are
+                            ratio = counter / nb_event 
+                            needed = 1.05 * to_decay[particle.pdg]/ratio - counter
+                            needed = min(1000, max(needed, 1000))
+                            with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
+                                new_file = generate_events(particle.pdg, needed, mg5, [decay_file_nb])
+                            to_event[particle.pdg].update(new_file)
+                            decay_file = to_event[particle.pdg][decay_file_nb]
+                            continue
+                        if helicity == decay[0].helicity or helicity==9 or \
+                                            self.options["spinmode"] == "none":
+                            break # use that event
+                        # not valid event store it for later
+                        if helicity not in bufferedEvents:
+                            bufferedEvents[helicity] = [decay]
+                        elif len(bufferedEvents[helicity]) < 200:
+                            # only add to the buffering if the buffer is not too large
+                            bufferedEvents[helicity].append(decay)
+                # now that we have the event make the merge
+                particle.add_decay(decay)
+            # change the weight associate to the event
+            event.wgt *= self.branching_ratio
+            wgts = event.parse_reweight()
+            for key in wgts:
+                wgts[key] *= self.branching_ratio
+            # all particle have been decay if needed
+            output_lhe.write(str(event))
+        output_lhe.write('</LesHouchesEvent>\n')        
+                    
     
     def load_model(self, name, use_mg_default, complex_mass=False):
         """load the model"""

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2014-07-17 13:34:15 +0000
+++ Template/LO/Cards/run_card.dat	2015-04-10 09:55:25 +0000
@@ -19,58 +19,59 @@
 #*********************************************************************
 # Tag name for the run (one word)                                    *
 #*********************************************************************
-  tag_1     = run_tag ! name of the run 
+  %(run_tag)s     = run_tag ! name of the run 
 #*********************************************************************
 # Run to generate the grid pack                                      *
 #*********************************************************************
-  .false.     = gridpack  !True = setting up 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.    *
 #*********************************************************************
-  10000 = nevents ! Number of unweighted events requested 
-      0       = iseed   ! rnd seed (0=assigned automatically=default))
+  %(nevents)s = nevents ! Number of unweighted events requested 
+  %(iseed)s   = iseed   ! rnd seed (0=assigned automatically=default))
 #*********************************************************************
 # Collider type and energy                                           *
 # lpp: 0=No PDF, 1=proton, -1=antiproton, 2=photon from proton,      *
 #                                         3=photon from electron     *
 #*********************************************************************
-        1     = lpp1    ! beam 1 type 
-        1     = lpp2    ! beam 2 type
-     6500     = ebeam1  ! beam 1 total energy in GeV
-     6500     = ebeam2  ! beam 2 total energy in GeV
+     %(lpp1)s        = lpp1    ! beam 1 type 
+     %(lpp2)s        = lpp2    ! beam 2 type
+     %(ebeam1)s     = ebeam1  ! beam 1 total energy in GeV
+     %(ebeam2)s     = ebeam2  ! beam 2 total energy in GeV
 #*********************************************************************
 # Beam polarization from -100 (left-handed) to 100 (right-handed)    *
 #*********************************************************************
-        0     = polbeam1 ! beam polarization for beam 1
-        0     = polbeam2 ! beam polarization for beam 2
+     %(polbeam1)s     = polbeam1 ! beam polarization for beam 1
+     %(polbeam2)s     = polbeam2 ! beam polarization for beam 2
 #*********************************************************************
 # PDF CHOICE: this automatically fixes also alpha_s and its evol.    *
 #*********************************************************************
- 'nn23lo1'    = pdlabel     ! PDF set                                     
- 230000    = lhaid     ! if pdlabel=lhapdf, this is the lhapdf number
+     %(pdlabel)s    = pdlabel     ! PDF set                                     
+     %(lhaid)s    = lhaid     ! if pdlabel=lhapdf, this is the lhapdf number
 #*********************************************************************
 # Renormalization and factorization scales                           *
 #*********************************************************************
- F        = fixed_ren_scale  ! if .true. use fixed ren scale
- F        = fixed_fac_scale  ! if .true. use fixed fac scale
- 91.1880  = scale            ! fixed ren scale
- 91.1880  = dsqrt_q2fact1    ! fixed fact scale for pdf1
- 91.1880  = dsqrt_q2fact2    ! fixed fact scale for pdf2
- 1        = scalefact        ! scale factor for event-by-event scales
+ %(fixed_ren_scale)s = fixed_ren_scale  ! if .true. use fixed ren scale
+ %(fixed_fac_scale)s        = fixed_fac_scale  ! if .true. use fixed fac scale
+ %(scale)s  = scale            ! fixed ren scale
+ %(dsqrt_q2fact1)s  = dsqrt_q2fact1    ! fixed fact scale for pdf1
+ %(dsqrt_q2fact2)s  = dsqrt_q2fact2    ! fixed fact scale for pdf2
+ %(dynamical_scale_choice)s = dynamical_scale_choice ! Select one of the preselect dynamical choice
+ %(scalefact)s  = scalefact        ! scale factor for event-by-event scales
 #*********************************************************************
 # Matching - Warning! ickkw > 1 is still beta
 #*********************************************************************
- 0        = ickkw            ! 0 no matching, 1 MLM, 2 CKKW matching
- 1        = highestmult      ! for ickkw=2, highest mult group
- 1        = ktscheme         ! for ickkw=1, 1 Durham kT, 2 Pythia pTE
- 1        = alpsfact         ! scale factor for QCD emission vx
- F        = chcluster        ! cluster only according to channel diag
- T        = pdfwgt           ! for ickkw=1, perform pdf reweighting
- 5        = asrwgtflavor     ! highest quark flavor for a_s reweight
- T        = clusinfo         ! include clustering tag in output
- 3.0      = lhe_version       ! Change the way clustering information pass to shower.        
+ %(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
+ %(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.        
 #*********************************************************************
 #**********************************************************
 #
@@ -78,25 +79,25 @@
 # Automatic ptj and mjj cuts if xqcut > 0
 # (turn off for VBF and single top processes)
 #**********************************************************
-   T  = auto_ptj_mjj  ! Automatic setting of ptj and mjj
+   %(auto_ptj_mjj)s  = auto_ptj_mjj  ! Automatic setting of ptj and mjj
 #**********************************************************
 #                                                                    
 #**********************************
 # BW cutoff (M+/-bwcutoff*Gamma)
 #**********************************
-  15  = bwcutoff      ! (M+/-bwcutoff*Gamma)
+  %(bwcutoff)s  = bwcutoff      ! (M+/-bwcutoff*Gamma)
 #**********************************************************
 # Apply pt/E/eta/dr/mij cuts on decay products or not
 # (note that etmiss/ptll/ptheavy/ht/sorted cuts always apply)
 #**********************************************************
-   T  = cut_decays    ! Cut decay products 
+   %(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.
 #*************************************************************
-   0  = nhel          ! Number of helicities used per event
+   %(nhel)s  = nhel          ! Number of helicities used per event
 #*******************                                                 
 # Standard Cuts
 #*******************                                                 
@@ -104,159 +105,157 @@
 #*********************************************************************
 # Minimum and maximum pt's (for max, -1 means no cut)                *
 #*********************************************************************
- 20  = ptj       ! minimum pt for the jets 
-  0  = ptb       ! minimum pt for the b 
- 10  = pta       ! minimum pt for the photons 
- 10  = ptl       ! minimum pt for the charged leptons 
-  0  = misset    ! minimum missing Et (sum of neutrino's momenta)
-  0  = ptheavy   ! minimum pt for one heavy final state
- 1.0 = ptonium   ! minimum pt for the quarkonium states
- -1  = ptjmax    ! maximum pt for the jets
- -1  = ptbmax    ! maximum pt for the b
- -1  = ptamax    ! maximum pt for the photons
- -1  = ptlmax    ! maximum pt for the charged leptons
- -1  = missetmax ! maximum missing Et (sum of neutrino's momenta)
+ %(ptj)s  = ptj       ! minimum pt for the jets 
+ %(ptb)s  = ptb       ! minimum pt for the b 
+ %(pta)s  = pta       ! minimum pt for the photons 
+ %(ptl)s  = ptl       ! minimum pt for the charged leptons 
+ %(misset)s  = misset    ! minimum missing Et (sum of neutrino's momenta)
+ %(ptheavy)s  = ptheavy   ! minimum pt for one heavy final state
+ %(ptjmax)s  = ptjmax    ! maximum pt for the jets
+ %(ptbmax)s  = ptbmax    ! maximum pt for the b
+ %(ptamax)s  = ptamax    ! maximum pt for the photons
+ %(ptlmax)s  = ptlmax    ! maximum pt for the charged leptons
+ %(missetmax)s  = missetmax ! maximum missing Et (sum of neutrino's momenta)
 #*********************************************************************
 # Minimum and maximum E's (in the center of mass frame)              *
 #*********************************************************************
-  0  = ej     ! minimum E for the jets 
-  0  = eb     ! minimum E for the b 
-  0  = ea     ! minimum E for the photons 
-  0  = el     ! minimum E for the charged leptons 
- -1   = ejmax ! maximum E for the jets
- -1   = ebmax ! maximum E for the b
- -1   = eamax ! maximum E for the photons
- -1   = elmax ! maximum E for the charged leptons
+  %(ej)s  = ej     ! minimum E for the jets 
+  %(eb)s  = eb     ! minimum E for the b 
+  %(ea)s  = ea     ! minimum E for the photons 
+  %(el)s  = el     ! minimum E for the charged leptons 
+  %(ejmax)s   = ejmax ! maximum E for the jets
+ %(ebmax)s   = ebmax ! maximum E for the b
+ %(eamax)s   = eamax ! maximum E for the photons
+ %(elmax)s   = elmax ! maximum E for the charged leptons
 #*********************************************************************
 # Maximum and minimum absolute rapidity (for max, -1 means no cut)   *
 #*********************************************************************
-   5  = etaj    ! max rap for the jets 
-  -1  = etab    ! max rap for the b
- 2.5  = etaa    ! max rap for the photons 
- 2.5  = etal    ! max rap for the charged leptons 
- 0.6  = etaonium ! max rap for the quarkonium states
-   0  = etajmin ! min rap for the jets
-   0  = etabmin ! min rap for the b
-   0  = etaamin ! min rap for the photons
-   0  = etalmin ! main rap for the charged leptons
+  %(etaj)s = etaj    ! max rap for the jets 
+  %(etab)s  = etab    ! max rap for the b
+ %(etaa)s  = etaa    ! max rap for the photons 
+ %(etal)s  = etal    ! max rap for the charged leptons 
+ %(etajmin)s  = etajmin ! min rap for the jets
+ %(etabmin)s  = etabmin ! min rap for the b
+ %(etaamin)s  = etaamin ! min rap for the photons
+ %(etalmin)s  = etalmin ! main rap for the charged leptons
 #*********************************************************************
 # Minimum and maximum DeltaR distance                                *
 #*********************************************************************
- 0.4 = drjj    ! min distance between jets 
- 0   = drbb    ! min distance between b's 
- 0.4 = drll    ! min distance between leptons 
- 0.4 = draa    ! min distance between gammas 
- 0   = drbj    ! min distance between b and jet 
- 0.4 = draj    ! min distance between gamma and jet 
- 0.4 = drjl    ! min distance between jet and lepton 
- 0   = drab    ! min distance between gamma and b 
- 0   = drbl    ! min distance between b and lepton 
- 0.4 = dral    ! min distance between gamma and lepton 
- -1  = drjjmax ! max distance between jets
- -1  = drbbmax ! max distance between b's
- -1  = drllmax ! max distance between leptons
- -1  = draamax ! max distance between gammas
- -1  = drbjmax ! max distance between b and jet
- -1  = drajmax ! max distance between gamma and jet
- -1  = drjlmax ! max distance between jet and lepton
- -1  = drabmax ! max distance between gamma and b
- -1  = drblmax ! max distance between b and lepton
- -1  = dralmax ! maxdistance between gamma and lepton
+ %(drjj)s = drjj    ! min distance between jets 
+ %(drbb)s   = drbb    ! min distance between b's 
+ %(drll)s = drll    ! min distance between leptons 
+ %(draa)s = draa    ! min distance between gammas 
+ %(drbj)s   = drbj    ! min distance between b and jet 
+ %(draj)s = draj    ! min distance between gamma and jet 
+ %(drjl)s = drjl    ! min distance between jet and lepton 
+ %(drab)s   = drab    ! min distance between gamma and b 
+ %(drbl)s   = drbl    ! min distance between b and lepton 
+ %(dral)s = dral    ! min distance between gamma and lepton 
+ %(drjjmax)s  = drjjmax ! max distance between jets
+ %(drbbmax)s  = drbbmax ! max distance between b's
+ %(drllmax)s  = drllmax ! max distance between leptons
+ %(draamax)s  = draamax ! max distance between gammas
+ %(drbjmax)s  = drbjmax ! max distance between b and jet
+ %(drajmax)s  = drajmax ! max distance between gamma and jet
+ %(drjlmax)s  = drjlmax ! max distance between jet and lepton
+ %(drabmax)s  = drabmax ! max distance between gamma and b
+ %(drblmax)s  = drblmax ! max distance between b and lepton
+ %(dralmax)s  = dralmax ! maxdistance between gamma and lepton
 #*********************************************************************
 # Minimum and maximum invariant mass for pairs                       *
 # WARNING: for four lepton final state mmll cut require to have      *
 #          different lepton masses for each flavor!                  *           
 #*********************************************************************
- 0   = mmjj    ! min invariant mass of a jet pair 
- 0   = mmbb    ! min invariant mass of a b pair 
- 0   = mmaa    ! min invariant mass of gamma gamma pair
- 0   = mmll    ! min invariant mass of l+l- (same flavour) lepton pair
- -1  = mmjjmax ! max invariant mass of a jet pair
- -1  = mmbbmax ! max invariant mass of a b pair
- -1  = mmaamax ! max invariant mass of gamma gamma pair
- -1  = mmllmax ! max invariant mass of l+l- (same flavour) lepton pair
+ %(mmjj)s   = mmjj    ! min invariant mass of a jet pair 
+ %(mmbb)s   = mmbb    ! min invariant mass of a b pair 
+ %(mmaa)s   = mmaa    ! min invariant mass of gamma gamma pair
+ %(mmll)s   = mmll    ! min invariant mass of l+l- (same flavour) lepton pair
+ %(mmjjmax)s  = mmjjmax ! max invariant mass of a jet pair
+ %(mmbbmax)s  = mmbbmax ! max invariant mass of a b pair
+ %(mmaamax)s  = mmaamax ! max invariant mass of gamma gamma pair
+ %(mmllmax)s  = mmllmax ! max invariant mass of l+l- (same flavour) lepton pair
 #*********************************************************************
 # Minimum and maximum invariant mass for all letpons                 *
 #*********************************************************************
- 0   = mmnl    ! min invariant mass for all letpons (l+- and vl) 
- -1  = mmnlmax ! max invariant mass for all letpons (l+- and vl) 
+ %(mmnl)s   = mmnl    ! min invariant mass for all letpons (l+- and vl) 
+ %(mmnlmax)s  = mmnlmax ! max invariant mass for all letpons (l+- and vl) 
 #*********************************************************************
 # Minimum and maximum pt for 4-momenta sum of leptons                *
 #*********************************************************************
- 0   = ptllmin  ! Minimum pt for 4-momenta sum of leptons(l and vl)
- -1  = ptllmax  ! Maximum pt for 4-momenta sum of leptons(l and vl)
+ %(ptllmin)s   = ptllmin  ! Minimum pt for 4-momenta sum of leptons(l and vl)
+ %(ptllmax)s  = ptllmax  ! Maximum pt for 4-momenta sum of leptons(l and vl)
 #*********************************************************************
 # Inclusive cuts                                                     *
 #*********************************************************************
- 0  = xptj ! minimum pt for at least one jet  
- 0  = xptb ! minimum pt for at least one b 
- 0  = xpta ! minimum pt for at least one photon 
- 0  = xptl ! minimum pt for at least one charged lepton 
+ %(xptj)s  = xptj ! minimum pt for at least one jet  
+ %(xptb)s  = xptb ! minimum pt for at least one b 
+ %(xpta)s  = xpta ! minimum pt for at least one photon 
+ %(xptl)s  = xptl ! minimum pt for at least one charged lepton 
 #*********************************************************************
 # Control the pt's of the jets sorted by pt                          *
 #*********************************************************************
- 0   = ptj1min ! minimum pt for the leading jet in pt
- 0   = ptj2min ! minimum pt for the second jet in pt
- 0   = ptj3min ! minimum pt for the third jet in pt
- 0   = ptj4min ! minimum pt for the fourth jet in pt
- -1  = ptj1max ! maximum pt for the leading jet in pt 
- -1  = ptj2max ! maximum pt for the second jet in pt
- -1  = ptj3max ! maximum pt for the third jet in pt
- -1  = ptj4max ! maximum pt for the fourth jet in pt
- 0   = cutuse  ! reject event if fails any (0) / all (1) jet pt cuts
+ %(ptj1min)s   = ptj1min ! minimum pt for the leading jet in pt
+ %(ptj2min)s   = ptj2min ! minimum pt for the second jet in pt
+ %(ptj3min)s   = ptj3min ! minimum pt for the third jet in pt
+ %(ptj4min)s   = ptj4min ! minimum pt for the fourth jet in pt
+ %(ptj1max)s  = ptj1max ! maximum pt for the leading jet in pt 
+ %(ptj2max)s  = ptj2max ! maximum pt for the second jet in pt
+ %(ptj3max)s  = ptj3max ! maximum pt for the third jet in pt
+ %(ptj4max)s  = ptj4max ! maximum pt for the fourth jet in pt
+ %(cutuse)s   = cutuse  ! reject event if fails any (0) / all (1) jet pt cuts
 #*********************************************************************
 # Control the pt's of leptons sorted by pt                           *
 #*********************************************************************
- 0   = ptl1min ! minimum pt for the leading lepton in pt
- 0   = ptl2min ! minimum pt for the second lepton in pt
- 0   = ptl3min ! minimum pt for the third lepton in pt
- 0   = ptl4min ! minimum pt for the fourth lepton in pt
- -1  = ptl1max ! maximum pt for the leading lepton in pt 
- -1  = ptl2max ! maximum pt for the second lepton in pt
- -1  = ptl3max ! maximum pt for the third lepton in pt
- -1  = ptl4max ! maximum pt for the fourth lepton in pt
+ %(ptl1min)s   = ptl1min ! minimum pt for the leading lepton in pt
+ %(ptl2min)s   = ptl2min ! minimum pt for the second lepton in pt
+ %(ptl3min)s   = ptl3min ! minimum pt for the third lepton in pt
+ %(ptl4min)s   = ptl4min ! minimum pt for the fourth lepton in pt
+ %(ptl1max)s  = ptl1max ! maximum pt for the leading lepton in pt 
+ %(ptl2max)s  = ptl2max ! maximum pt for the second lepton in pt
+ %(ptl3max)s  = ptl3max ! maximum pt for the third lepton in pt
+ %(ptl4max)s  = ptl4max ! maximum pt for the fourth lepton in pt
 #*********************************************************************
 # Control the Ht(k)=Sum of k leading jets                            *
 #*********************************************************************
- 0   = htjmin ! minimum jet HT=Sum(jet pt)
- -1  = htjmax ! maximum jet HT=Sum(jet pt)
- 0   = ihtmin  !inclusive Ht for all partons (including b)
- -1  = ihtmax  !inclusive Ht for all partons (including b)
- 0   = ht2min ! minimum Ht for the two leading jets
- 0   = ht3min ! minimum Ht for the three leading jets
- 0   = ht4min ! minimum Ht for the four leading jets
- -1  = ht2max ! maximum Ht for the two leading jets
- -1  = ht3max ! maximum Ht for the three leading jets
- -1  = ht4max ! maximum Ht for the four leading jets
+ %(htjmin)s   = htjmin ! minimum jet HT=Sum(jet pt)
+ %(htjmax)s  = htjmax ! maximum jet HT=Sum(jet pt)
+ %(ihtmin)s   = ihtmin  !inclusive Ht for all partons (including b)
+ %(ihtmax)s  = ihtmax  !inclusive Ht for all partons (including b)
+ %(ht2min)s   = ht2min ! minimum Ht for the two leading jets
+ %(ht3min)s   = ht3min ! minimum Ht for the three leading jets
+ %(ht4min)s   = ht4min ! minimum Ht for the four leading jets
+ %(ht2max)s  = ht2max ! maximum Ht for the two leading jets
+ %(ht3max)s  = ht3max ! maximum Ht for the three leading jets
+ %(ht4max)s  = ht4max ! maximum Ht for the four leading jets
 #***********************************************************************
 # Photon-isolation cuts, according to hep-ph/9801442                   *
 # When ptgmin=0, all the other parameters are ignored                  *
 # When ptgmin>0, pta and draj are not going to be used                 *
 #***********************************************************************
-   0 = ptgmin ! Min photon transverse momentum
- 0.4 = R0gamma ! Radius of isolation code
- 1.0 = xn ! n parameter of eq.(3.4) in hep-ph/9801442
- 1.0 = epsgamma ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
- .true. = isoEM ! isolate photons from EM energy (photons and leptons)
+ %(ptgmin)s = ptgmin ! Min photon transverse momentum
+ %(r0gamma)s = R0gamma ! Radius of isolation code
+ %(xn)s = xn ! n parameter of eq.(3.4) in hep-ph/9801442
+ %(epsgamma)s = epsgamma ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
+ %(isoem)s = isoEM ! isolate photons from EM energy (photons and leptons)
 #*********************************************************************
 # WBF cuts                                                           *
 #*********************************************************************
- 0   = xetamin ! minimum rapidity for two jets in the WBF case  
- 0   = deltaeta ! minimum rapidity for two jets in the WBF case 
+ %(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                                                      *
 #*********************************************************************
- -1    =  ktdurham        
- 0.4  =  dparameter 
+ %(ktdurham)s    =  ktdurham        
+ %(dparameter)s  =  dparameter 
 #*********************************************************************
 # maximal pdg code for quark to be considered as a light jet         *
 # (otherwise b cuts are applied)                                     *
 #*********************************************************************
- 4 = maxjetflavor    ! Maximum jet pdg code
+ %(maxjetflavor)s = maxjetflavor    ! Maximum jet pdg code
 #*********************************************************************
 # Jet measure cuts                                                   *
 #*********************************************************************
- 0   = xqcut   ! minimum kt jet measure between partons
+ %(xqcut)s   = xqcut   ! minimum kt jet measure between partons
 #*********************************************************************
 #
 #*********************************************************************
@@ -265,16 +264,16 @@
 #          meaningful ONLY if plotted taking matchscale              *
 #          reweighting into account!                                 *
 #*********************************************************************
-   F  = use_syst      ! Enable systematics studies
+   %(use_syst)s  = use_syst      ! Enable systematics studies
 #
 #**************************************
 # Parameter of the systematics study
 #  will be used by SysCalc (if installed)
 #**************************************                                  
 #
-0.5 1 2 = sys_scalefact  # factorization/renormalization scale factor
-0.5 1 2 = sys_alpsfact  # \alpha_s emission scale factors
-30 50 = sys_matchscale # variation of merging scale
+%(sys_scalefact)s = sys_scalefact  # factorization/renormalization scale factor
+%(sys_alpsfact)s = sys_alpsfact  # \alpha_s emission scale factors
+%(sys_matchscale)s = sys_matchscale # variation of merging scale
 # PDF sets and number of members (0 or none for all members).
-CT10nlo.LHgrid = sys_pdf # matching scales
+%(sys_pdf)s = sys_pdf # matching scales
 # MSTW2008nlo68cl.LHgrid 1  = sys_pdf

=== removed file 'Template/LO/Source/PDF/pdfwrap_lhapdf.f'
--- Template/LO/Source/PDF/pdfwrap_lhapdf.f	2012-11-07 05:57:53 +0000
+++ Template/LO/Source/PDF/pdfwrap_lhapdf.f	1970-01-01 00:00:00 +0000
@@ -1,68 +0,0 @@
-      subroutine pdfwrap
-      implicit none
-C
-C     INCLUDE
-C
-      include 'pdf.inc'
-      include '../alfas.inc'
-      real*8 zmass
-      data zmass/91.188d0/
-      Character*150 LHAPath
-      character*20 parm(20)
-      double precision value(20)
-      real*8 alphasPDF
-      external alphasPDF
-
-
-c-------------------
-c     START THE CODE
-c-------------------      
-
-c     initialize the pdf set
-      call FindPDFPath(LHAPath)
-      CALL SetPDFPath(LHAPath)
-      value(1)=lhaid
-      parm(1)='DEFAULT'
-      call pdfset(parm,value)
-      call GetOrderAs(nloop)
-      nloop=nloop+1  
-      asmz=alphasPDF(zmass)
-      
-      return
-      end
- 
-
-      subroutine FindPDFPath(LHAPath)
-c********************************************************************
-c generic subroutine to open the table files in the right directories
-c********************************************************************
-      implicit none
-c
-      Character LHAPath*150,up*3
-      data up/'../'/
-      logical exists
-      integer i
-
-c     first try in the current directory
-      LHAPath='PDFsets'
-      Inquire(File=LHAPath, exist=exists)
-      if(exists)return
-c     then try one directory up
-      LHAPath=up//LHAPath
-      Inquire(File=LHAPath, exist=exists)
-      if(exists)return
-c     finally try in the lib directory
-      LHAPath='lib/PDFsets'
-      Inquire(File=LHAPath, exist=exists)
-      if(exists)return
-      do i=1,6
-         LHAPath=up//LHAPath
-         Inquire(File=LHAPath, exist=exists)
-         if(exists)return
-      enddo
-      print*,'Could not find PDFsets directory, quitting'
-      stop 1
-      
-      return
-      end
-

=== modified file 'Template/LO/Source/cuts.inc'
--- Template/LO/Source/cuts.inc	2014-11-22 14:28:19 +0000
+++ Template/LO/Source/cuts.inc	2015-04-10 09:55:25 +0000
@@ -26,7 +26,7 @@
       REAL*8 DRBJmax,DRAJmax,DRJLmax,DRABmax,DRBLmax,DRALmax
       REAL*8 MMJJmax,MMLLmax,MMAAmax,MMBBmax !max inv mass 
       REAL*8 MMNL,MMNLMAX ! invariant mass of all leptons
-      REAL*8 cutuse
+      integer cutuse
       REAL*8 ptj1min,ptj2min,ptj3min,ptj4min
       REAL*8 ptj1max,ptj2max,ptj3max,ptj4max
       REAL*8 ptl1min,ptl2min,ptl3min,ptl4min

=== modified file 'Template/LO/Source/dsample.f'
--- Template/LO/Source/dsample.f	2014-11-05 17:34:19 +0000
+++ Template/LO/Source/dsample.f	2015-04-10 09:55:25 +0000
@@ -76,6 +76,10 @@
       integer                   neventswritten
       common /to_eventswritten/ neventswritten
 
+      integer th_nunwgt
+      double precision th_maxwgt
+      common/theoretical_unwgt_max/th_maxwgt, th_nunwgt
+
 c
 c     External
 c
@@ -105,6 +109,9 @@
       if (nsteps .lt. 1) nsteps=1
       nwrite = itmax*ncall/nsteps
 
+C     Fix for 2>1 process where ndim is 2 and not 1
+      ninvar = max(2,ninvar)
+
       call sample_init(ndim,ncall,itmax,ninvar,nconfigs)
       call graph_init
       do i=1,itmax
@@ -200,24 +207,32 @@
       if (chi2 .gt. 1) tsigma=tsigma*sqrt(chi2)
 c     JA 02/2011 Added twgt to results.dat to allow event generation in
 c     first iteration for gridpack runs
+C     OM 02/2015 Added maxwgt (target of the secondary unweight) to allow splitted
+C        generation of event.
       if (icor .eq. 0) then
-         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5,e13.5)')tmean,tsigma,0.0,
-     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt, trmean
+         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5,3e13.5,i9)')tmean,tsigma, 0.0,
+     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt, trmean, 
+     &     maxwgt, th_maxwgt, th_nunwgt
       else
-         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5,e13.5)')tmean,0.0,tsigma,
-     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt, trmean
+         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5,3e13.5,i9)')tmean,0.0,tsigma,
+     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt, trmean, 
+     &     maxwgt, th_maxwgt, th_nunwgt
       endif
 c      do i=1,cur_it-1
       do i=cur_it-itsum,cur_it-1
          write(66,'(i4,5e15.5)') i,xmean(i),xsigma(i),xeff(i),xwmax(i),xrmean(i)
       enddo
+c     Write out MadLoop statistics, if any
+      call output_run_statistics(66)
       flush(66)
       close(66, status='KEEP')
       else
          open(unit=66,file='results.dat',status='unknown')
-         write(66,'(3e12.5,2i9,i5,i9,3e10.3)')0.,0.,0.,kevent,nw,
-     &     1,0,0.,0.,0.
+         write(66,'(3e12.5,2i9,i5,i9,5e10.3,i9)')0.,0.,0.,kevent,nw,
+     &     1,0,0.,0.,0.,0.,0.,0
          write(66,'(i4,5e15.5)') 1,0.,0.,0.,0.,0.
+c        Write out MadLoop statistics, if any
+         call output_run_statistics(66)
          flush(66)
          close(66, status='KEEP')
 
@@ -361,25 +376,31 @@
       if (nun .lt. 0) nun=-nun   !Case when wrote maximun number allowed
       if (chi2 .gt. 1) tsigma=tsigma*sqrt(chi2)
 c     JA 02/2011 Added twgt to results.dat to allow event generation in
-c     first iteration for gridpack runs
+c     first iteration for gridpack runs +02/2015 maxwgt 
       if (icor .eq. 0) then
-         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5,e13.5)')tmean,tsigma,0.0,
-     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt,trmean
+         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5,3e13.5, i9)')tmean,tsigma,0.0,
+     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt,trmean, 
+     &    maxwgt, th_maxwgt, th_nunwgt
       else
-         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5,e13.5)')tmean,0.0,tsigma,
-     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt,trmean
+         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5,3e13.5,i9)')tmean,0.0,tsigma,
+     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt,trmean,
+     &    maxwgt, th_maxwgt, th_nunwgt
       endif
 c      do i=1,cur_it-1
       do i=cur_it-itsum,cur_it-1
          write(66,'(i4,5e15.5)') i,xmean(i),xsigma(i),xeff(i),xwmax(i),xrmean(i)
       enddo
+c     Write out MadLoop statistics, if any
+      call output_run_statistics(66)      
       flush(66)
       close(66, status='KEEP')
       else
          open(unit=66,file='results.dat',status='unknown')
-         write(66,'(3e12.5,2i9,i5,i9,3e10.3)')0.,0.,0.,kevent,nw,
-     &     1,0,0.,0.,0.
+         write(66,'(3e12.5,2i9,i5,i9,5e10.3,i9)')0.,0.,0.,kevent,nw,
+     &     1,0,0.,0.,0.,0.,0.,0
          write(66,'(i4,5e15.5)') 1,0.,0.,0.,0.,0.
+c        Write out MadLoop statistics, if any
+         call output_run_statistics(66)
          flush(66)
          close(66, status='KEEP')
 
@@ -387,6 +408,71 @@
 
       end
 
+      subroutine output_run_statistics(outUnit)
+c***********************************************************************
+c     Writes out the madloop runtime statistics to the unit in argument
+c***********************************************************************
+      use StringCast
+      implicit none
+c
+c     Arguments
+c
+      integer outUnit
+C
+C     Local
+C
+      double precision t_after
+c
+c     Global
+c
+      INTEGER U_RETURN_CODES(0:9)
+      INTEGER T_RETURN_CODES(0:9)
+      INTEGER H_RETURN_CODES(0:9)
+      DOUBLE PRECISION AVG_TIMING
+      DOUBLE PRECISION MAX_PREC, MIN_PREC
+      INTEGER N_EVALS
+      DATA U_RETURN_CODES/10*0/
+      DATA T_RETURN_CODES/10*0/
+      DATA H_RETURN_CODES/10*0/
+      DATA MAX_PREC /-1.0d0/
+      DATA MIN_PREC /1.0d99/
+      DATA AVG_TIMING/0.0d0/
+      DATA N_EVALS/0/
+      COMMON/MADLOOPSTATS/AVG_TIMING,MAX_PREC,MIN_PREC,N_EVALS,
+     &       U_RETURN_CODES,T_RETURN_CODES,H_RETURN_CODES
+
+      DOUBLE PRECISION CUMULATED_TIMING
+      DATA CUMULATED_TIMING/0.0d0/
+      COMMON/GENERAL_STATS/CUMULATED_TIMING
+
+c-----
+c  Begin Code
+c-----
+      call cpu_time(t_after)
+      CUMULATED_TIMING = t_after - CUMULATED_TIMING
+
+      if (N_EVALS.eq.0) then
+        return
+      endif
+      
+      write(outUnit,*) '<run_statistics> '
+      write(outUnit,33) '<u_return_code>',U_RETURN_CODES,'</u_return_code>'
+      write(outUnit,33) '<t_return_code>',T_RETURN_CODES,'</t_return_code>'
+      write(outUnit,33) '<h_return_code>',H_RETURN_CODES,'</h_return_code>'
+      write(outUnit,*) '<average_time>'//trim(toStr_real(AVG_TIMING))
+     & //'</average_time>'
+      write(outUnit,*) '<cumulated_time>'//trim(toStr_real(CUMULATED_TIMING))
+     & //'</cumulated_time>'
+      write(outUnit,*) '<max_prec>'//trim(toStr_real(MAX_PREC))//'</max_prec>'
+      write(outUnit,*) '<min_prec>'//trim(toStr_real(MIN_PREC))//'</min_prec>'
+      write(outUnit,*) '<n_evals>'//trim(toStr_int(N_EVALS))//'</n_evals>'   
+      write(outUnit,*) '</run_statistics>'
+      
+33    FORMAT( a15,i12,',',i12',',i12',',i12',',i12',
+     &        ',i12',',i12',',i12',',i12',',i12,a16)
+
+      end subroutine
+
       subroutine sample_writehtm()
 c***********************************************************************
 c     Writes out results of run in html format
@@ -473,6 +559,7 @@
 c
       include 'genps.inc'
       include 'maxconfigs.inc'
+      include 'run.inc'
 c
 c     Arguments
 c
@@ -482,9 +569,13 @@
 c     Local
 c
       integer i, j
+      integer get_maxsproc
 c
 c     Global
 c
+      double precision force_max_wgt
+      common/unwgt_secondary_max/force_max_wgt
+
       integer                                      nsteps
       character*40          result_file,where_file
       common /sample_status/result_file,where_file,nsteps
@@ -510,12 +601,20 @@
       logical            flat_grid
       common/to_readgrid/flat_grid                !Tells if grid read from file
 
+      double precision twgt, maxwgt,swgt(maxevents)
+      integer                             lun, nw, itminx
+      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itminx
+      
       integer              icor
       common/to_correlated/icor
 
       logical               zooming
       common /to_zoomchoice/zooming
 
+      logical read_grid_file
+      data read_grid_file/.False./
+      common/read_grid_file/read_grid_file
+
       data use_cut/2/            !Grid: 0=fixed , 1=standard, 2=non-zero
       data ituple/1/             !1=htuple, 2=sobel 
       data Minvar(1,1)/-1/       !No special variable mapping
@@ -615,11 +714,12 @@
 c
       flat_grid=.true.
       open(unit=25,file='ftn25',status='unknown',err=102)
-      read(25,fmt='(4f20.17)', err=101, end=101)
-     .     ((grid(2,i,j),i=1,ng),j=1,maxinvar)
-c      write(*,*) 'Got the grid OK, now getting alpha'
-c      read(25,fmt='(4f20.17)',err=101,end=101)(alpha(i),i=1,maxconfigs)
+      read(25,*, err=1011, end=1012)
+     .     ((grid(2,i,j),i=1,ng),j=1,invar)
+      read(25,*) twgt, force_max_wgt
+      call read_discrete_grids(25)
       write(*,*) 'Grid read from file'
+      read_grid_file=.true.
       flat_grid=.false.
       close(25)
 c
@@ -650,14 +750,21 @@
          enddo
       endif
       goto 103
+ 1011  write(*,*) 'fail to open file'
+       goto 101
+ 1012  write(*,*) 'fail to read data'
+       goto 101
  101  close(25)
 c      write(*,*) 'Tried reading it',i,j
  102  write(*,*) 'Error opening grid'
+
 c
 c     Unable to read grid, using uniform grid and equal points in
 c     each configuration
 c
-      write(*,*) 'Using Uniform Grid!'
+      read_grid_file=.false.
+      write(*,*) 'Using Uniform Grid!', maxinvar
+      force_max_wgt = -1d0
       do j = 1, maxinvar
          do i = 1, ng
             grid(2, i, j) = xgmin+ (xgmax-xgmin)*(i / dble(ng))**1
@@ -674,6 +781,14 @@
 c      write(*,*) 'Forwarding random number generator'
 
  103  write(*,*) 'Grid defined OK'
+
+C     sanity check that we have a minimal number of event      
+      if ( MC_GROUPED_SUBPROC )then
+         events = max(events, maxtries)
+      else 
+         events = max(events, 2*maxtries*get_maxsproc())
+      endif
+
       end
 
       subroutine setgrid(j,xo,a,itype)
@@ -848,6 +963,97 @@
       endif
       end
 
+      subroutine write_discrete_grids(stream_id, grid_type)
+c************************************************************************
+c     Write out the grid using the DiscreteSampler module
+c************************************************************************
+      use DiscreteSampler          
+      implicit none
+      integer, intent(in)                           :: stream_id
+      character(len=*)                              :: grid_type
+      logical MC_grouped_subproc
+      common/to_MC_grouped_subproc/MC_grouped_subproc
+      INTEGER                    ISUM_HEL
+      LOGICAL                    MULTI_CHANNEL
+      COMMON/TO_MATRIX/ISUM_HEL, MULTI_CHANNEL
+c
+c     Begin code
+c
+
+      if (ISUM_HEL.ne.0.and.DS_get_dim_status('Helicity').ge.1) then
+        call DS_write_grid(stream_id, dim_name='Helicity', 
+     &                                              grid_type=grid_type)
+      elseif(ISUM_HEL.eq.0)then
+        call write_good_hel(stream_id)  
+      endif
+      
+
+
+      if(MC_grouped_subproc.and.
+     &             DS_get_dim_status('grouped_processes').ge.1) then
+        call DS_write_grid(stream_id, dim_name='grouped_processes', 
+     &                                              grid_type=grid_type)
+      endif
+
+      end subroutine write_discrete_grids
+
+      subroutine read_discrete_grids(stream_id)
+c************************************************************************
+c     Write out the grid using the DiscreteSampler module
+c************************************************************************
+      use DiscreteSampler          
+      implicit none
+      integer, intent(in)                           :: stream_id
+      INTEGER                    ISUM_HEL
+      LOGICAL                    MULTI_CHANNEL
+      COMMON/TO_MATRIX/ISUM_HEL, MULTI_CHANNEL
+      
+      if (ISUM_HEL.eq.0)then
+         call read_good_hel(stream_id)
+      endif
+      call DS_load_grid(stream_id)
+
+      end subroutine read_discrete_grids
+
+      subroutine sample_get_discrete_x(wgt,picked_bin,iconfig,dim_name)
+c************************************************************************
+c     Returns maxdim random numbers between 0 and 1, and the wgt
+c     associated with this set of points, and the iteration number
+c     This routine chooses the point within the range specified by
+c     xmin and xmax for dimension j in configuration ipole
+c************************************************************************
+      use DiscreteSampler
+
+      implicit none
+      include 'genps.inc'
+C     Subroutine arguments
+      integer picked_bin
+      character(len=*) dim_name
+      real*8 wgt
+C     This variable iconfig is what corresponds to ipole in sample_get_x
+C     and is used for random number generation
+      integer iconfig
+C     Local variables
+      real*8 jacobian
+      real*8 rdm
+      integer dummy
+c     
+c      Begin code
+c
+C     Fetch a random number bewteen 0.0 and 1.0
+c     The fourth argument is not used and therefore a dummy
+      dummy = 0
+      call ntuple(rdm,0.0d0,1.0d0,dummy,iconfig)
+C     Pick a point using the DiscreteSampler module
+      CALL DS_get_point(dim_name, rdm, picked_bin, jacobian, 'norm') 
+C     Store the helicity sampling jacobian so that it can be divided out
+c     of wgt later when adding an entry to the DiscreteSampler helicity
+c      grid. Also we don't want to multiply wgt by it yet since this is
+c     taken care of at the level of matrix<i> already.
+      hel_jacobian = jacobian
+      
+      end subroutine sample_get_discrete_x
+
       subroutine sample_get_x(wgt, x, j, ipole, xmin, xmax)
 c************************************************************************
 c     Returns maxdim random numbers between 0 and 1, and the wgt
@@ -1168,6 +1374,93 @@
 
       end
 
+C
+C     Subroutine to take care of the update of the discrete grids
+C     (used for helicity and the matrix<i> choice in the grouped case
+C     as implented in the DiscreteSampler module.
+C
+      subroutine add_entry_to_discrete_dimensions(wgt)
+      use DiscreteSampler
+      implicit none
+c
+c     Constants
+c
+      include 'genps.inc'
+c
+c     Arguments
+c
+      double precision wgt
+c
+c     Local
+c
+c
+c     Global
+c
+      INTEGER                    ISUM_HEL
+      LOGICAL                    MULTI_CHANNEL
+      COMMON/TO_MATRIX/ISUM_HEL, MULTI_CHANNEL
+      logical cutsdone, cutspassed
+      COMMON/TO_CUTSDONE/CUTSDONE,CUTSPASSED
+c
+c     Begin code
+c
+c       It is important to divide the wgt stored in the grid by the 
+c       corresponding jacobian otherwise it flattens the sampled
+c       distribution.
+C       Also, if HEL_PICKED is equal to -1, it means that MadEvent
+C       is in the initialization stage where all helicity were probed
+c       and added individually to the grid directly by matrix<i>.f so
+c       that they shouldn't be added here.
+        if(ISUM_HEL.ne.0.and.HEL_PICKED.ne.-1.and.
+     &                            (.NOT.CUTSDONE.or.CUTSPASSED)) then
+          call DS_add_entry('Helicity',HEL_PICKED,(wgt/hel_jacobian))
+        endif
+
+      end subroutine add_entry_to_discrete_dimensions
+
+C
+C     Subroutine to take care of the update of the discrete grids
+C     (used for helicity and the matrix<i> choice in the grouped case
+C     as implented in the DiscreteSampler module.
+C
+      subroutine update_discrete_dimensions()
+      use DiscreteSampler
+      implicit none
+c
+c     Constants
+c
+      include 'genps.inc'
+c
+c     Arguments
+c
+c
+c     Local
+c
+      type(SampledDimension) tmp_dim      
+c
+c     Global
+c
+      INTEGER                    ISUM_HEL
+      LOGICAL                    MULTI_CHANNEL
+      COMMON/TO_MATRIX/ISUM_HEL, MULTI_CHANNEL
+      logical MC_grouped_subproc
+      common/to_MC_grouped_subproc/MC_grouped_subproc
+c
+c     Begin code
+c
+      if(ISUM_HEL.ne.0) then
+        call DS_update_grid('Helicity', filterZeros=.True.)
+        tmp_dim = DS_get_dimension(ref_grid,'Helicity')
+C       Security in case of all helicity vanishing (G1 of gg > qq )
+        if (size(tmp_dim%bins).eq.0) then
+          call none_pass(-1)
+        endif
+      endif
+      if(MC_grouped_subproc.and.DS_get_dim_status('grouped_processes').ne.-1) then
+        call DS_update_grid('grouped_processes', filterZeros=.True.)
+      endif
+
+      end subroutine update_discrete_dimensions
 
       subroutine sample_put_point(wgt, point, iteration,ipole)
 c**************************************************************************
@@ -1190,12 +1483,12 @@
 c
 c     Local
 c
-      integer i, j, k, knt, non_zero, nun,itsum
+      integer i, j, k, knt, nun,itsum
       double precision vol,xnmin,xnmax,tot,xdum,tmp1,chi2tmp
       double precision rc, dr, xo, xn, x(maxinvar), dum(ng-1)
       save vol,knt
       double precision  chi2
-      save chi2, non_zero
+      save chi2
       double precision wmax1,ddumb
       save wmax1
       double precision twgt1,xchi2,xxmean,tmeant,tsigmat
@@ -1211,6 +1504,13 @@
 c
 c     Global
 c
+      integer th_nunwgt
+      double precision th_maxwgt
+      common/theoretical_unwgt_max/th_maxwgt, th_nunwgt
+
+      double precision force_max_wgt
+      common/unwgt_secondary_max/force_max_wgt
+
       double precision    accur
       common /to_accuracy/accur
 
@@ -1221,8 +1521,8 @@
       common/to_result/mean,rmean,sigma
 
       double precision grid2(0:ng,maxinvar)
-      integer               inon_zero(ng,maxinvar)
-      common/to_grid2/grid2,inon_zero
+      integer               inon_zero(ng,maxinvar), non_zero
+      common/to_grid2/grid2,inon_zero,non_zero
 
       double precision tmean, trmean, tsigma
       integer             dim, events, itm, kn, cur_it, invar, configs
@@ -1275,6 +1575,7 @@
 c-----
 c  Begin Code
 c-----
+
       if (first_time) then
          first_time = .false.
          twgt1 = 0d0       !
@@ -1305,6 +1606,27 @@
       endif
 
       if (iteration .eq. cur_it) then
+c        Add the current point to the DiscreteSamplerGrid
+         call add_entry_to_discrete_dimensions(wgt)
+         if (kn.eq.0) then
+            ! ensure that all cumulative variable are at zero (usefull for reset)
+            twgt1 = 0d0         !
+            iavg = 0            !Vars for averging to increase err estimate
+            navg = 1            !
+            wmax1= 99d99
+            wmax = -1d0
+            mean = 0d0
+            rmean = 0d0
+            sigma = 0d0
+            chi2 = 0d0
+            non_zero = 0
+            vol = 1d0 / dble(events * itm)
+            knt = events
+            do i=1,maxconfigs
+               psect(i)=0d0
+            enddo
+         endif
+
          kn = kn + 1
          if (.true.) then       !Average points to increase error estimate
             twgt1=twgt1+dabs(wgt)     !This doesn't change anything should remove
@@ -1389,11 +1711,37 @@
 c     Now if done with an iteration, print out stats, rebin, reset
 c         
 c         if (kn .eq. events) then
-         if (kn .ge. max_events .and. non_zero .lt. 5) then
+         if (kn .ge. max_events .and. non_zero .le. 5) then
             call none_pass(max_events)
          endif
          if (non_zero .eq. events .or. (kn .gt. 200*events .and.
      $        non_zero .gt. 5)) then
+
+c          # special mode where we store information to combine them
+           if(use_cut.eq.-2)then
+                open(unit=22, file="grid_information")
+                write(22,*) non_zero, ng, invar
+                write(22,*) ((grid(1,i,j),i=1,ng),j=1,invar)
+                write(22,*) ((grid(2,i,j),i=1,ng),j=1,invar)
+                write(22,*) ((inon_zero(i,j),i=1,ng),j=1,invar)
+                write(22,*) (xmin(j), j=1,invar)
+                write(22,*) (xmax(j), j=1,invar)
+                write(22,*) mean, rmean, sigma, wmax, kn,events, force_max_wgt
+c               In order not to write out the reference grid but just
+c               the points which were added for this last iteration,
+c               we write out the discrete 'running' grids before the
+c               update of the reference grid.
+                call write_discrete_grids(22,'all')
+                close(22)
+           endif
+
+C
+C    Now updated the discrete dimensions of the DiscreteSampler module
+C    used for sampling helicity configurations and matrix<i> config
+C    choice in the grouped case.
+C
+           call update_discrete_dimensions()
+
             mean=mean*dble(events)/dble(non_zero)
             rmean=rmean*dble(events)/dble(non_zero)
             twgt1=twgt1*dble(events)/dble(non_zero)
@@ -1420,6 +1768,7 @@
                sigma = (sigma/vol/vol-non_zero*mean*mean*navg)  !knt replaced by non_zero
      .              / dble(knt-1) / dble(knt)
             else
+
                sigma = (sigma/vol/vol - knt*mean*mean)
      .              / dble(knt-1) / dble(knt)
             endif
@@ -1511,13 +1860,16 @@
 c     &           wmax*(dble(non_zero)/dble(kn)),
 c     &           dble(non_zero)/dble(kn)*100.,'%'
 c            close(22)
+
 c------
 c    Here we will double the number of events requested for the next run
 c-----
  23         events = 2 * events
             vol = 1d0/dble(events*itm)
             knt = events
-            twgt = mean / (dble(itm)*dble(events))
+            if (use_cut.ne.-2) then
+              twgt = mean / (dble(itm)*dble(events))
+            endif
 c            write(*,*) 'New number of events',events,twgt
 
             mean = 0d0
@@ -1526,12 +1878,13 @@
             cur_it = cur_it + 1
             kn = 0
             wmax = -1d0
+
 c
 c     do so adjusting of weights according to number of events in bin
 c
             do j=1,invar
                do i = 1, ng
-                  if (use_cut .ne. 2 .and.
+                  if (abs(use_cut) .ne. 2 .and.
      &                use_cut .ne. 3 .and. use_cut .ne. 5)
      $                 inon_zero(i,j) = 0
                   if (use_cut .eq. 3) grid(1,i,j)=grid2(i,j)
@@ -1698,7 +2051,11 @@
             if (tsigma .gt. 0d0 .and. cur_it .gt. itmin .and. accur .gt. 0d0) then
 
                xxmean = tmean/tsigma
-               xchi2 = (chi2/xxmean/xxmean-tsigma)/dble(cur_it-2)               
+               if (cur_it.ne.2)then
+                  xchi2 = (chi2/xxmean/xxmean-tsigma)/dble(cur_it-2)               
+               else
+                  xchi2 = 0d0
+               endif
                write(*,'(a,4f8.3)') ' Accuracy: ',sqrt(xchi2/tsigma),
      &              accur,1/sqrt(tsigma),xchi2
 c               write(*,*) 'We got it',1d0/sqrt(tsigma), accur
@@ -1716,10 +2073,12 @@
                   write(*, 80) real(tmean), real(tsigma), real(trmean), real(chi2)
                   if (use_cut .ne. 0) then
                   open(26, file='ftn26',status='unknown')
-                  write(26,fmt='(4f20.17)')
-     $                 ((grid(2,i,j),i=1,ng),j=1,maxinvar)
-                  write(26,fmt='(4f20.17)') (alpha(i),i=1,maxconfigs)
-                  close(26)
+                  write(26,fmt='(4f21.17)')
+     $                 ((grid(2,i,j),i=1,ng),j=1,invar)
+                  write(26,*) twgt, force_max_wgt
+c                  write(26,fmt='(4f21.16)') (alpha(i),i=1,maxconfigs)
+                  call write_discrete_grids(26,'ref')
+                  close(26)                  
                   endif
                   call sample_writehtm()
 c                  open(unit=22,file=result_file,status='old',
@@ -1727,7 +2086,11 @@
 c                  write(22, 80) real(tmean), real(tsigma), real(chi2)
 c 122              close(22)
                   tsigma = tsigma*sqrt(chi2)  !This gives the 68% confidence cross section
-                  call store_events
+                  if (use_cut.eq.-2)then
+                    call store_events(force_max_wgt, .False.)
+                  else
+                     call store_events(-1d0, .True.)
+                  endif
                   cur_it = itm+2
                   return
                endif
@@ -1742,7 +2105,12 @@
 c
 c               nun = n_unwgted()
 c               write(*,*) 'Estimated events',nun, accur
-               call store_events
+               if (use_cut.eq.-2) then
+                  call store_events(force_max_wgt, .False.)
+               else
+                  call store_events(-1d0, .True.)
+               endif
+
                nun = neventswritten
 c               tmp1 = tmean / tsigma
 c               chi2tmp = (chi2/tmp1/tmp1-tsigma)/dble(cur_it-2)
@@ -1775,9 +2143,11 @@
                   write(*, 80) real(tmean), real(tsigma), real(chi2)
                   if (use_cut .ne. 0) then
                   open(26, file='ftn26',status='unknown')
-                  write(26,fmt='(4f20.17)')
-     $                 ((grid(2,i,j),i=1,ng),j=1,maxinvar)
-                  write(26,fmt='(4f20.17)') (alpha(i),i=1,maxconfigs)
+                  write(26,fmt='(4f21.17)')
+     $                 ((grid(2,i,j),i=1,ng),j=1,invar)
+                  write(26,*) twgt, force_max_wgt
+c                  write(26,fmt='(4f21.17)') (alpha(i),i=1,maxconfigs)
+                  call write_discrete_grids(26,'ref')
                   close(26)
                   endif
                   call sample_writehtm()
@@ -1794,7 +2164,11 @@
 
 
             if (cur_it .gt. itm) then               
-               call store_events
+               if (use_cut.eq.-2)then
+                  call store_events(force_max_wgt, .False.)
+               else
+                  call store_events(-1d0, .True.)
+               endif
                tmean = tmean / tsigma
                trmean = trmean / tsigma
                chi2 = (chi2 / tmean / tmean - tsigma) / dble(itm - 1)
@@ -1806,9 +2180,11 @@
      .              13X,21HChi**2 per DoF.     =,f12.4/1X,79(1H-))
                if (use_cut .ne. 0) then
                open(26, file='ftn26',status='unknown')
-               write(26,fmt='(4f20.17)')
-     $              ((grid(2,i,j),i=1,ng),j=1,maxinvar)
-               write(26,fmt='(4f20.17)') (alpha(i),i=1,maxconfigs)
+               write(26,fmt='(4f21.17)')
+     $              ((grid(2,i,j),i=1,ng),j=1,invar)
+               write(26,*) twgt, force_max_wgt
+               call write_discrete_grids(26,'ref')
+c               write(26,fmt='(4f21.17)') (alpha(i),i=1,maxconfigs)
                close(26)
                endif
                call sample_writehtm()
@@ -1867,9 +2243,11 @@
  222  format(a10,I3,3x,a6,e10.4,a16,e10.3,e12.3,3x,f5.1,a1)
 
       open(unit=66,file='results.dat',status='unknown')
-      write(66,'(3e12.5,2i9,i5,i9,3e10.3)')0.,0.,0.,0,0,
-     &     0,1,0.,0.,0.
+      write(66,'(3e12.5,2i9,i5,i9,5e10.3,i9)')0.,0.,0.,0,0,
+     &     0,1,0.,0.,0.,0.,0.,0
       write(66,'(i4,5e15.5)') 1,0.,0.,0.,0.,0.
+c     Write out MadLoop statistics, if any
+      call output_run_statistics(66)
       flush(66)
       close(66, status='KEEP')
 
@@ -1879,6 +2257,10 @@
       write(67,*)
       close(67)
 
+      open(unit=67, file='grid_information')
+      write(67,*) ''
+      close(67)
+
       stop
       end
             
@@ -1912,7 +2294,7 @@
 c  Begin Code
 c-----
       kmax=k
-      do i=k+1,ng-1
+      do i=k+1,ng
          if (grid(1,i,j) .gt. 0d0) kmax=i
       enddo
       xo = grid(1,k,j)
@@ -1929,7 +2311,7 @@
       end do
 c      grid(1, ng, j) = (xn + xo) / 2d0  !Original without kmax stuff
       grid(1, kmax, j) = (xn + xo) / 2d0
-      x(j) = x(j) + grid(1, ng, j)
+      x(j) = x(j) + grid(1, kmax, j)
       end
 
       double precision function xbin(y,j)
@@ -2097,5 +2479,44 @@
       return
       end
 
+      subroutine reset_cumulative_variable()
+C       Reset to zero all the variable which evaluates the cross-section.
+C       grid information for the current-grid/non-zero entry/...
+C       This is used to avoid the (small) bias introduce in the first iteration
+C       Due to the initialization of the helicity sum.
+      implicit none
+      include 'genps.inc'
+
+      double precision grid2(0:ng,maxinvar)
+      integer               inon_zero(ng,maxinvar), non_zero
+      common/to_grid2/grid2,inon_zero, non_zero
+      double precision    grid(2, ng, 0:maxinvar)
+      common /data_grid/ grid
+
+      double precision tmean, trmean, tsigma
+      integer             dim, events, itm, kn, cur_it, invar, configs
+      common /sample_common/
+     .     tmean, trmean, tsigma, dim, events, itm, kn, cur_it, invar, configs
+
+
+C     LOCAL
+      integer i,j
+
+      write(*,*) "RESET CUMULATIVE VARIABLE"
+      non_zero = 0
+      do j=1,maxinvar
+         do i=1,ng -1
+            inon_zero = 0
+            grid2(i,j) = 0
+            grid(1,i,j) = 0
+         enddo
+      enddo
+      tmean = 0.0
+      trmean = 0.0
+      tsigma = 0.0
+      kn = 0
+      return 
+      end
+
 
 

=== modified file 'Template/LO/Source/gen_ximprove.f'
--- Template/LO/Source/gen_ximprove.f	2013-12-09 02:09:49 +0000
+++ Template/LO/Source/gen_ximprove.f	2015-04-10 09:55:25 +0000
@@ -262,6 +262,7 @@
       double precision xt(lmaxconfigs),elimit
       double precision yerr,ysec,rerr
       logical fopened
+
 c-----
 c  Begin Code
 c-----
@@ -424,6 +425,10 @@
 c      write(26,15) 'if [[ "$PBS_O_WORKDIR" != "" ]]; then' 
 c      write(26,15) '    cd $PBS_O_WORKDIR'
 c      write(26,15) 'fi'
+      write(26,15) 'if [[ -e MadLoop5_resources.tar.gz && ! -e MadLoop5_resources ]]; then'
+      write(26,15) 'tar -xzf MadLoop5_resources.tar.gz'
+      write(26,15) 'fi'
+
       write(26,15) 'k=run1_app.log'
       write(lun,15) 'script=' // fname
 c      write(lun,15) 'rm -f wait.$script >& /dev/null'
@@ -613,9 +618,9 @@
 c
          mjobs = (goal_lum*xsec(io(np))*1000 / MaxEventsPerJob + 0.9)
 c         write(*,*) "Working on Channel ",i,io(np),xt(np), goal_lum*xsec(io(np))*1000 /maxeventsperjob
-         if (mjobs .gt. 26)  then
+         if (mjobs .gt. 130)  then
             write(*,*) 'Error in gen_ximprove.f, too many events requested ',mjobs*maxeventsperjob
-            mjobs=26
+            mjobs=130
          endif
          if (mjobs .lt. 1 .or. .not. split_channels)  mjobs=1
 c
@@ -651,9 +656,24 @@
 
          ip = index(gn(io(np)),'/')
          if (mjobs .gt. 1) then
-            write(26,'(3a)') 'j=',gn(io(np))(1:ip-1),cjobs(ijob:ijob)
+
+            if (ip.eq.3) then
+                write(26,'(a2,a2,a,i1)') 'j=',gn(io(np))(1:ip-1),cjobs(MODULO(ijob-1,26)+1:MODULO(ijob-1,26)+1),
+     &                                              ijob/26
+            else if(ip.eq.4) then
+                write(26,'(a2,a3,a,i1)') 'j=',gn(io(np))(1:ip-1),cjobs(MODULO(ijob-1,26)+1:MODULO(ijob-1,26)+1),
+     &                                              ijob/26
+            else if(ip.eq.5) then
+               write(26,'(a2,a4,a,i1)') 'j=',gn(io(np))(1:ip-1),cjobs(MODULO(ijob-1,26)+1:MODULO(ijob-1,26)+1),
+     &                                              ijob/26
+            else if(ip.eq.6) then
+               write(26,'(a2,a5,a,i1)') 'j=',gn(io(np))(1:ip-1),cjobs(MODULO(ijob-1,26)+1:MODULO(ijob-1,26)+1),
+     &                                              ijob/26
+           else
+               stop 1
+           endif
          else
-            write(26,'(3a)') 'j=',gn(io(np))(1:ip-1)
+             write(26,'(3a)') 'j=',gn(io(np))(1:ip-1)
          endif
 c
 c     Now write the commands

=== modified file 'Template/LO/Source/genps.inc'
--- Template/LO/Source/genps.inc	2012-04-16 08:48:53 +0000
+++ Template/LO/Source/genps.inc	2015-04-10 09:55:25 +0000
@@ -32,4 +32,9 @@
       REAL*8 LIMHEL
       PARAMETER(LIMHEL=1e-6)
       INTEGER MAXTRIES
-      PARAMETER(MAXTRIES=100)
+      PARAMETER(MAXTRIES=10)
+C     To pass the helicity configuration chosen by the DiscreteSampler to 
+C     matrix<i>.f
+      double precision hel_jacobian
+	  INTEGER HEL_PICKED
+	  COMMON/HEL_PICKED/HEL_PICKED,hel_jacobian

=== modified file 'Template/LO/Source/run.inc'
--- Template/LO/Source/run.inc	2014-06-27 13:04:42 +0000
+++ Template/LO/Source/run.inc	2015-04-10 09:55:25 +0000
@@ -6,9 +6,10 @@
 c
       real*8          scale,scalefact,alpsfact
       logical         fixed_ren_scale,fixed_fac_scale,fixed_couplings,hmult
-      integer         ickkw,nhmult,asrwgtflavor
+      integer         ickkw,nhmult,asrwgtflavor, dynamical_scale_choice
       common/to_scale/scale,scalefact,alpsfact,fixed_ren_scale,fixed_fac_scale,
-     $                fixed_couplings,ickkw,nhmult,hmult,asrwgtflavor
+     $                fixed_couplings,ickkw,nhmult,hmult,asrwgtflavor,
+     $                dynamical_scale_choice
 c
 c     Collider
 c
@@ -56,3 +57,9 @@
       logical clusinfo
       double precision lhe_version
       COMMON/TO_LHEFORMAT/lhe_version,clusinfo
+
+c     
+C     Controls wheter to perform Monte-Carlo sampling over grouped subprocesses
+C
+      logical MC_grouped_subproc
+      common/to_MC_grouped_subproc/MC_grouped_subproc

=== modified file 'Template/LO/Source/setrun.f'
--- Template/LO/Source/setrun.f	2015-01-29 12:27:51 +0000
+++ Template/LO/Source/setrun.f	2015-04-10 09:55:25 +0000
@@ -77,6 +77,9 @@
 c
       include 'run_card.inc'
 
+c     if no matching ensure that no pdfreweight are done
+      if (ickkw.eq.0) pdfwgt = .false.
+
       q2fact(1) = sf1**2      ! fact scale**2 for pdf1
       q2fact(2) = sf2**2      ! fact scale**2 for pdf2     
 

=== modified file 'Template/LO/SubProcesses/cuts.f'
--- Template/LO/SubProcesses/cuts.f	2014-03-10 17:54:39 +0000
+++ Template/LO/SubProcesses/cuts.f	2015-04-10 09:55:25 +0000
@@ -1005,13 +1005,15 @@
 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
-      if(xqcut.gt.0d0.or.ickkw.gt.0.or.scale.eq.0.or.q2fact(1).eq.0)then
+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
          if(debug) write (*,*) ' setclscales -> fails'
          passcuts=.false.
          return
        endif
-      endif
+c      endif
 
 c     Set couplings in model files
       if(.not.fixed_ren_scale.or..not.fixed_couplings) then

=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2013-11-05 18:08:13 +0000
+++ Template/LO/SubProcesses/genps.f	2015-04-10 09:55:25 +0000
@@ -23,6 +23,9 @@
 c     OUTPUTS: wgt      == updated weight after choosing points
 c              x        == points choosen from sample grid
 c              p        == transformed points call is f(p(x))
+c     COMMON:
+c           hel_picked  == Modified integer in gen_ps.inc to pass the
+c                          chosen helicity configuration to matrix<i>.f
 c
 c**************************************************************************
       implicit none
@@ -48,10 +51,18 @@
 c
 c     Global
 c
+      INTEGER                    ISUM_HEL
+      LOGICAL                    MULTI_CHANNEL
+      COMMON/TO_MATRIX/ISUM_HEL, MULTI_CHANNEL
 c-----
 c  Begin Code
 c-----
       call gen_mom(iconfig,mincfig,maxcfig,invar,wgt,x,p)
+C     Pick the helicity configuration from the DiscreteSampler if user
+C     decided to perform MC over helicity configurations.
+      if(ISUM_HEL.ne.0) then
+        call sample_get_discrete_x(wgt,hel_picked,iconfig,'Helicity')
+      endif
       end
 
       subroutine gen_mom(iconfig,mincfig,maxcfig,invar,wgt,x,p1)

=== modified file 'Template/LO/SubProcesses/makefile'
--- Template/LO/SubProcesses/makefile	2013-11-19 14:23:47 +0000
+++ Template/LO/SubProcesses/makefile	2015-04-10 09:55:25 +0000
@@ -7,9 +7,19 @@
 BINDIR = ../../bin/
 PROG   = madevent
 
-LINKLIBS =  -L../../lib/ -ldhelas -ldsample -lmodel -lgeneric -lpdf -lcernlib $(lhapdf)
-
-LIBS = $(LIBDIR)libdhelas.$(libext) $(LIBDIR)libdsample.$(libext) $(LIBDIR)libgeneric.$(libext) $(LIBDIR)libpdf.$(libext) $(LIBDIR)libmodel.$(libext) $(LIBDIR)libcernlib.$(libext)
+ifneq ("$(wildcard ../MadLoop_makefile_definitions)","")
+    include ../MadLoop_makefile_definitions
+else
+    LINK_LOOP_LIBS =
+    LOOP_LIBS =
+    LOOP_INCLUDE =
+    LINK_MADLOOP_LIB =
+    MADLOOP_LIB =
+endif
+
+LINKLIBS =  -L../../lib/ -ldsample -lmodel -lgeneric -lpdf -lcernlib $(lhapdf) $(LINK_MADLOOP_LIB) $(LINK_LOOP_LIBS) -L../../lib/ -ldhelas
+
+LIBS = $(LIBDIR)libdhelas.$(libext) $(LIBDIR)libdsample.$(libext) $(LIBDIR)libgeneric.$(libext) $(LIBDIR)libpdf.$(libext) $(LIBDIR)libmodel.$(libext) $(LIBDIR)libcernlib.$(libext) $(MADLOOP_LIB) $(LOOP_LIBS)
 
 # Source files
 
@@ -18,17 +28,16 @@
 	 idenparts.o \
          $(patsubst %.f,%.o,$(wildcard auto_dsig*.f)) \
          $(patsubst %.f,%.o,$(wildcard matrix*.f))
-SYMMETRY = symmetry.o setcuts.o cuts.o cluster.o myamp.o genps.o \
-	   initcluster.o setscales.o reweight.o get_color.o idenparts.o \
-           $(patsubst %.f,%.o,$(wildcard matrix*.f))
+
+SYMMETRY = symmetry.o idenparts.o 
 
 # Binaries
 
 $(PROG): $(PROCESS) auto_dsig.o $(LIBS)
 	$(FC) $(FFLAGS) -o $(PROG) $(PROCESS) $(LINKLIBS)
 
-gensym: $(SYMMETRY) configs.inc $(LIBS)
-	$(FC) $(FFLAGS) -o gensym $(SYMMETRY) $(LINKLIBS)
+gensym: $(SYMMETRY) configs.inc $(LIBDIR)libmodel.$(libext) $(LIBDIR)libgeneric.$(libext)
+	$(FC) $(FFLAGS) -o gensym $(SYMMETRY) -L../../lib/ -lmodel -lgeneric
 
 $(LIBDIR)libmodel.$(libext): ../../Cards/param_card.dat
 	cd ../../Source/MODEL; make
@@ -39,12 +48,13 @@
 $(LIBDIR)libpdf.$(libext): 
 	cd ../../Source/PDF; make
 
+# Add source so that the compiler finds the DiscreteSampler module.
+%.o: %.f
+	$(FC) $(FFLAGS) -c $< -I../../Source/
+
 # Dependencies
 
-driver.f: genps.inc qmass.inc
-setcuts.f: qmass.inc
-qmass.inc: nexternal.inc
-	   touch qmass.inc
+driver.f: genps.inc
 symmetry.o: genps.inc nexternal.inc configs.inc run_config.inc
 genps.o: genps.inc nexternal.inc configs.inc
 cuts.o: genps.inc nexternal.inc pmass.inc

=== modified file 'Template/LO/SubProcesses/myamp.f'
--- Template/LO/SubProcesses/myamp.f	2015-01-05 14:10:20 +0000
+++ Template/LO/SubProcesses/myamp.f	2015-04-10 09:55:25 +0000
@@ -651,7 +651,7 @@
 
       write(*,*),'Impossible BW configuration'
       open(unit=66,file='results.dat',status='unknown')
-      write(66,'(3e12.5,2i9,i5,i9,3e10.3)')0.,0.,0.,0,0,1,0,0.,0.,0.
+      write(66,'(3e12.5,2i9,i5,i9,4e10.3)')0.,0.,0.,0,0,1,0,0.,0.,0.,0.
       write(66,'(i4,5e15.5)') 1,0.,0.,0.,0.,0.
       close(66)
       end

=== added file 'Template/LO/SubProcesses/refine.sh'
--- Template/LO/SubProcesses/refine.sh	1970-01-01 00:00:00 +0000
+++ Template/LO/SubProcesses/refine.sh	2015-04-10 09:55:25 +0000
@@ -0,0 +1,68 @@
+#!/bin/bash
+
+# For support of LHAPATH in cluster mode
+if [ $CLUSTER_LHAPATH ]; then
+  export LHAPATH=$CLUSTER_LHAPATH;
+fi
+
+if [[ -e MadLoop5_resources.tar.gz && ! -e MadLoop5_resources ]]; then
+tar -xzf MadLoop5_resources.tar.gz
+fi
+k=%(name)s_app.log
+script=%(script_name)s                         
+
+grid_directory=%(base_directory)s
+j=%(directory)s
+     if [[ ! -e $j ]]; then
+          mkdir $j
+          if [[ -e $grid_directory/ftn26 ]];then
+             cp $grid_directory/ftn26 $j/ftn25
+          fi 
+          if [[ ! -e ../../SubProcesses ]];then
+          	 if [[ -e ftn26 ]]; then
+          	 	cp ./ftn26 $j/ftn25
+          	 fi
+          fi	
+     fi
+     
+     cd $j
+     rm -f $k
+     rm -f moffset.dat >& /dev/null
+      echo   %(offset)s  > moffset.dat
+     if  [[ -e ftn26 ]]; then
+          cp ftn26 ftn25
+     fi
+     # create the input file
+         echo "    %(nevents)s       %(maxiter)s       %(miniter)s" >& input_sg.txt
+         echo "    %(precision)s" >> input_sg.txt
+     if [[ ! -e ftn25 ]]; then
+         echo "2" >> input_sg.txt   # grid refinement
+         echo "1" >> input_sg.txt   # suppress amplitude
+
+     else
+         echo "%(grid_refinment)s" >> input_sg.txt
+         echo "1" >> input_sg.txt
+     fi
+     echo "%(nhel)s" >> input_sg.txt
+     echo "%(channel)s" >> input_sg.txt
+
+     # run the executable. The loop is design to avoid
+     # filesystem problem (executable not found)
+     for((try=1;try<=16;try+=1)); 
+     do
+         ../madevent 2>&1 >> $k <input_sg.txt | tee -a $k;
+         if [ -s $k ]
+         then
+             break
+         else
+             echo $try > fail.log 
+         fi
+     done
+     echo "" >> $k; echo "ls status:" >> $k; ls >> $k
+     cat $k >> log.txt
+     
+     if [[ -e ftn26 ]]; then
+         cp ftn26 ftn25
+     fi
+     cd ../
+

=== added file 'Template/LO/SubProcesses/refine_splitted.sh'
--- Template/LO/SubProcesses/refine_splitted.sh	1970-01-01 00:00:00 +0000
+++ Template/LO/SubProcesses/refine_splitted.sh	2015-04-10 09:55:25 +0000
@@ -0,0 +1,72 @@
+#!/bin/bash
+
+# For support of LHAPATH in cluster mode
+if [ $CLUSTER_LHAPATH ]; then
+  export LHAPATH=$CLUSTER_LHAPATH;
+fi
+if [[ -e MadLoop5_resources.tar.gz && ! -e MadLoop5_resources ]]; then
+tar -xzf MadLoop5_resources.tar.gz
+fi
+k=run1_app.log
+script=refine_splitted.sh
+# Argument
+# 1st argument the Directory name
+grid_directory=$1;
+# 2st argument Directory where to find the grid input
+base_directory=$2; 
+# 3st argument the offset
+offset=$3;
+
+
+# prepare the directory where to run
+if [[ ! -e $grid_directory ]]; then
+    # Do not exists     
+    mkdir $grid_directory;
+else
+   rm -rf $grid_directory/$k;
+   rm -rf $grid_directory/input_app.txt;
+   rm -rf $grid_directory/ftn25;
+   rm -rf $grid_directory/ftn26;
+fi
+# handle input file
+if [[ -e $base_directory ]]; then
+    cp $base_directory/ftn26 $grid_directory/ftn25;
+    cp $base_directory/input_app.txt $grid_directory/input_app.txt;
+elif [[ -e ./ftn26 ]]; then
+    cp ./ftn26 $grid_directory/ftn25;
+    cp ./input_app.txt $grid_directory/input_app.txt;
+else
+    exit 1;
+fi
+
+# Move to the running directory
+cd $grid_directory;
+
+# Put the correct offset
+rm -f moffset.dat >& /dev/null;
+echo   $offset > moffset.dat;
+
+# run the executable. The loop is design to avoid
+# filesystem problem (executable not found)
+for((try=1;try<=16;try+=1)); 
+do
+    ../madevent 2>&1 >> $k <input_app.txt | tee -a $k;
+    if [ -s $k ]
+    then
+        break
+    else
+        echo $try > fail.log 
+    fi
+done
+echo "" >> $k; echo "ls status:" >> $k; ls >> $k
+# Perform some cleaning to keep less file on disk/transfer less file.
+subdir=${grid_directory##*_}
+if [[ $subdir -ne 1 &&  -s results.dat && $MG5DEBUG != true ]]; then
+	 rm -f ftn25 &> /dev/null
+         rm -f ftn26 &> /dev/null
+	 rm -f log.txt &> /dev/null
+         rm -f *.log &> /dev/null
+	 rm -f moffset.dat &> /dev/null
+	 rm -f fail.log &> /dev/null
+fi
+cd ../
\ No newline at end of file

=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f	2014-09-09 10:27:43 +0000
+++ Template/LO/SubProcesses/reweight.f	2015-04-10 09:55:25 +0000
@@ -538,9 +538,18 @@
       external isqcd, isjet, isparton, cluster, isjetvx, alphas, ifsno
       setclscales=.true.
 
-      if(ickkw.le.0.and.xqcut.le.0d0.and.q2fact(1).gt.0.and.scale.gt.0)
-     $     return
-
+      if(ickkw.le.0.and.xqcut.le.0d0.and.q2fact(1).gt.0.and.scale.gt.0) then
+         if(use_syst)then
+            s_scale=scale
+            n_qcd=nqcd(iconfig)
+            n_alpsem=0
+            do i=1,2
+               n_pdfrw(i)=0
+            enddo
+            s_rwfact=1d0
+         endif
+      return
+      endif
 c   
 c   Cluster the configuration
 c   
@@ -729,7 +738,7 @@
            endif
            if(ipart(2,imocl(n)).gt.2)then ! ipart(1) set and not IS line
 c             The ishft gives the FS particle corresponding to imocl
-              if(ipdgcl(ishft(1,ipart(2,imocl(n))-1),igraphs(1),iproc).ne.21.and.
+              if(ipdgcl(ishft(1,ipart(2,imocl(n))-1),igraphs(1),iproc).ne.21.or.
      $                                   ipdgcl(imocl(n),igraphs(1),iproc).ne.21) then
 c                 The second condition is to prevent the case of ggh where the gluon split in quark later.
 c                 The first quark is already remove so we shouldn't remove this one.      
@@ -993,7 +1002,7 @@
             q2fact(1) = pt2ijcl(jfirst(1))
       elseif(jcentral(2).eq.0)then
             q2fact(2) = pt2ijcl(jfirst(2))
-      elseif(ickkw.eq.2.or.pdfwgt)then
+      elseif(ickkw.eq.2.or.(pdfwgt.and.ickkw.gt.0))then
 c     Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
 c     f1(x1,pt2E) is given by DSIG, just need to set scale.
 c     Use the minimum scale found for fact scale in ME
@@ -1175,7 +1184,7 @@
          write(*,*) 'Set process number ',ipsel
       endif
 
-      if (use_syst.and.igraphs(1).eq.0) igraphs(1) = 1 ! happens if use_syst=T BUT fix scale
+      if (use_syst.and.igraphs(1).eq.0) igraphs(1) = iconfig ! happens if use_syst=T BUT fix scale
 c     Set incoming particle identities
       ipdgcl(1,igraphs(1),iproc)=idup(1,ipsel,iproc)
       ipdgcl(2,igraphs(1),iproc)=idup(2,ipsel,iproc)
@@ -1197,8 +1206,10 @@
          enddo
       endif
 
-      if(ickkw.le.0) goto 100
-
+      if(ickkw.le.0)then
+         asref=0 ! usefull for syscalc
+         goto 100
+      endif
 c   Preparing graph particle information (ipart, needed to keep track of
 c   external particle clustering scales)
       do i=1,nexternal
@@ -1210,7 +1221,7 @@
      $              max(pt2min,p(0,i)**2-p(1,i)**2-p(2,i)**2-p(3,i)**2)
             endif
             pt2pdf(ishft(1,i-1))=pt2prev(ishft(1,i-1))
-         else if(pdfwgt) then
+         else if(pdfwgt.and.ickkw.gt.0) then
             pt2pdf(ishft(1,i-1))=0d0
          endif
          ipart(1,ishft(1,i-1))=i
@@ -1323,7 +1334,7 @@
               endif
            endif
         endif
-        if(ickkw.eq.2.or.pdfwgt) then
+        if(ickkw.eq.2.or.(pdfwgt.and.ickkw.gt.0)) then
 c       Perform PDF and, if ickkw=2, Sudakov reweighting
           isvx=.false.
           do i=1,2
@@ -1495,7 +1506,7 @@
       if(ickkw.eq.2.and.lpp(1).eq.0.and.lpp(2).eq.0)then
          q2fact(1)=pt2min
          q2fact(2)=q2fact(1)
-      else if (ickkw.eq.1.and.pdfwgt) then
+      else if (ickkw.gt.0.and.pdfwgt) then
          q2fact(1)=q2bck(1)
          q2fact(2)=q2bck(2)         
          if (btest(mlevel,3))
@@ -1512,14 +1523,21 @@
 c     Set reweight factor for systematics studies
       if(use_syst)then
          s_rwfact = rewgt
+         
 c     Need to multiply by: initial PDF, alpha_s^n_qcd to get
 c     factor in front of matrix element
          do i=1,2
-            s_rwfact=s_rwfact*pdg2pdf(abs(lpp(IB(i))),
+            if (lpp(IB(i)).ne.0) then
+                s_rwfact=s_rwfact*pdg2pdf(abs(lpp(IB(i))),
      $           i_pdgpdf(1,i)*sign(1,lpp(IB(i))),
      $           s_xpdf(1,i),s_qpdf(1,i))
+            endif
          enddo
-         s_rwfact=s_rwfact*asref**n_qcd
+         if (asref.gt.0d0.and.n_qcd.le.nexternal)then
+            s_rwfact=s_rwfact*asref**n_qcd
+         else
+            s_rwfact=0d0
+         endif
       endif
 
       return

=== modified file 'Template/LO/SubProcesses/setcuts.f'
--- Template/LO/SubProcesses/setcuts.f	2014-10-07 13:31:57 +0000
+++ Template/LO/SubProcesses/setcuts.f	2015-04-10 09:55:25 +0000
@@ -110,7 +110,6 @@
 c     setup masses for the final-state particles
 c
       include 'pmass.inc'
-      include 'qmass.inc'
 
 C-----
 C  BEGIN CODE
@@ -185,7 +184,6 @@
 c
       do i=nincoming+1,nexternal
          do_cuts(i)=.true.
-         if(nincoming.eq.1) do_cuts(i)=.false.
          if(.not.cut_decays.and.from_decay(i)) do_cuts(i)=.false.
          is_a_j(i)=.false.
          is_a_l(i)=.false.
@@ -322,7 +320,7 @@
             endif
 c        PHOTON
             if(is_a_a(i))then
-                 etmin(i) = max(pta, ptgmin, eamax)
+                 etmin(i) = max(pta, ptgmin, ea)
                  SMIN = SMIN + etmin(i)
                  etmax(i)=ptamax
                  emin(i)=ea
@@ -458,7 +456,7 @@
       inclHtmin=ihtmin
       inclHtmax=ihtmax
 
-      jetor = cutuse.eq.0d0
+      jetor = cutuse.eq.0
 
 c
 c   EXTRA LEPTON CUTS

=== modified file 'Template/LO/SubProcesses/setscales.f'
--- Template/LO/SubProcesses/setscales.f	2014-05-26 01:30:36 +0000
+++ Template/LO/SubProcesses/setscales.f	2015-04-10 09:55:25 +0000
@@ -42,8 +42,43 @@
 c     start
 c----------
 
-      rscale=0d0
-
+      if (dynamical_scale_choice.eq.-1) then
+c         Cluster external states until reducing the system to a 2->2 topology whose transverse mass is used for setting the scale.
+c         This is not done in this file due to the clustering.
+         rscale=0d0
+      elseif(dynamical_scale_choice.eq.1) then
+c         Total transverse energy of the event.         
+          rscale=0d0
+          do i=3,nexternal
+             rscale=rscale+et(P(0,i))
+          enddo      
+      elseif(dynamical_scale_choice.eq.2) then
+c         sum of the transverse mass
+c         m^2+pt^2=p(0)^2-p(3)^2=(p(0)+p(3))*(p(0)-p(3))
+          rscale=0d0
+          do i=3,nexternal
+            rscale=rscale+dsqrt(max(0d0,(P(0,i)+P(3,i))*(P(0,i)-P(3,i))))
+          enddo
+          rscale=rscale
+      elseif(dynamical_scale_choice.eq.3) then
+c         sum of the transverse mass divide by 2
+c         m^2+pt^2=p(0)^2-p(3)^2=(p(0)+p(3))*(p(0)-p(3))
+          rscale=0d0
+          do i=3,nexternal
+            rscale=rscale+dsqrt(max(0d0,(P(0,i)+P(3,i))*(P(0,i)-P(3,i))))
+          enddo
+          rscale=rscale/2d0
+      elseif(dynamical_scale_choice.eq.4) then
+c         \sqrt(s), partonic energy
+          rscale=dsqrt(max(0d0,2d0*dot(P(0,1),P(0,2))))
+      elseif(dynamical_scale_choice.eq.0) then
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+cc      USER DEFINE SCALE: ENTER YOUR CODE HERE                                  cc
+cc      to use this code you need to set                                         cc
+cc                 dymamical_scale_choice to 0 in the run_card                   cc
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+         write(*,*) "scale not define by the user"
+         stop 21
 c
 c-some examples of dynamical scales
 c
@@ -64,36 +99,17 @@
 
 c      rscale=sqrt(rscale)
 
-c      rscale=rscale*scalefact
-
-
-c---------------------------------------
-c-- total transverse energy of the event 
-c---------------------------------------
-c     rscale=0d0
-c     do i=3,nexternal
-c      rscale=rscale+et(P(0,i))
-c     enddo
-
-c--------------------------------------
-c-- scale^2 = \sum_i  (pt_i^2+m_i^2)  
-c--------------------------------------
-c     rscale=0d0
-c     do i=3,nexternal
-c      rscale=rscale+pt(P(0,i))**2+dot(p(0,i),p(0,i))
-c     enddo
-c     rscale=dsqrt(rscale)
-
-c--------------------------------------
-c-- \sqrt(s): partonic energy
-c--------------------------------------
-c     rscale=dsqrt(2d0*dot(P(0,1),P(0,2)))
+
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+cc      USER DEFINE SCALE: END of USER CODE                                      cc
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      endif
 
       return
       end
 
 
-      subroutine set_fac_scale(P,q2fact)
+      subroutine set_fac_scale(P,q2factorization)
 c----------------------------------------------------------------------
 c     This is the USER-FUNCTION to calculate the factorization 
 c     scales^2 on event-by-event basis.
@@ -102,15 +118,16 @@
 
 c     INCLUDE and COMMON
 c
-c      include 'genps.inc'
+      include 'genps.inc'
       include 'nexternal.inc'
       include 'coupl.inc'
+      include 'run.inc'
 c--masses and poles
 c
 c     ARGUMENTS
 c      
       REAL*8 P(0:3,nexternal)
-      real*8 q2fact(2)
+      real*8 q2factorization(2)
 c
 c     EXTERNAL
 c
@@ -126,16 +143,21 @@
 c     start
 c----------
       
-      
-      q2fact(1)=0d0             !factorization scale**2 for pdf1
-      q2fact(2)=0d0             !factorization scale**2 for pdf2
-
-c      call set_ren_scale(P,q2fact(1))
-c
-c      q2fact(1)=q2fact(1)**2
-c
-c      q2fact(2)=q2fact(1)       !factorization scale**2 for pdf2
-      
+      if (dynamical_scale_choice.eq.-1) then
+c         Cluster external states until reducing the system to a 2->2 topology whose transverse mass is used for setting the scale.
+c         This is not done in this file due to the clustering.
+         q2factorization(1)=0d0          !factorization scale**2 for pdf1
+         q2factorization(2)=0d0          !factorization scale**2 for pdf2
+      elseif(dynamical_scale_choice.eq.0) then
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+cc      USER DEFINE SCALE: ENTER YOUR CODE HERE                                  cc
+cc      to use this code you need to set                                         cc
+cc                 dymamical_scale_choice to 0 in the run_card                   cc
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c         default: use the renormalization scale
+          call set_ren_scale(P,q2factorization(1))
+          q2factorization(1)=q2factorization(1)**2
+          q2factorization(2)=q2factorization(1)   !factorization scale**2 for pdf2
 
 c
 c-some examples of dynamical scales
@@ -144,26 +166,40 @@
 c---------------------------------------
 c-- total transverse energy of the event 
 c---------------------------------------
-c     q2fact(1)=0d0
+c     q2factorization(1)=0d0
 c     do i=3,nexternal
-c      q2fact(1)= q2fact(1)+et(P(0,i))**2
+c      q2factorization(1)= q2factorization(1)+et(P(0,i))**2
 c     enddo
-c     q2fact(2)=q2fact(1)  
+c     q2factorization(2)=q2factorization(1)  
 
 c--------------------------------------
 c-- scale^2 = \sum_i  (pt_i^2+m_i^2)  
 c--------------------------------------
-c     q2fact(1)=0d0
+c     q2factorization(1)=0d0
 c     do i=3,nexternal
-c      q2fact(1)=q2fact(1)+pt(P(0,i))**2+dot(p(0,i),p(0,i))
+c      q2factorization(1)=q2factorization(1)+pt(P(0,i))**2+dot(p(0,i),p(0,i))
 c     enddo
-c     q2fact(2)=q2fact(1)  
+c     q2factorization(2)=q2factorization(1)  
 
 c--------------------------------------
 c-- \sqrt(s): partonic energy
 c--------------------------------------
-c     q2fact(1)=2d0*dot(P(0,1),P(0,2))
-c     q2fact(2)=q2fact(1)  
+c     q2factorization(1)=2d0*dot(P(0,1),P(0,2))
+c     q2factorization(2)=q2factorization(1)  
+
+
+
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+cc      USER DEFINE SCALE: END of USER CODE                                      cc
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      else
+          call set_ren_scale(P,q2factorization(1))
+          q2factorization(1)=q2factorization(1)**2
+          q2factorization(2)=q2factorization(1)   !factorization scale**2 for pdf2
+      endif
+
+
+      
 
       
       return

=== added file 'Template/LO/SubProcesses/survey.sh'
--- Template/LO/SubProcesses/survey.sh	1970-01-01 00:00:00 +0000
+++ Template/LO/SubProcesses/survey.sh	2015-04-10 09:55:25 +0000
@@ -0,0 +1,78 @@
+#!/bin/bash
+
+# For support of LHAPATH in cluster mode
+if [ $CLUSTER_LHAPATH ]; then 
+  export LHAPATH=$CLUSTER_LHAPATH;
+fi
+
+if [[ -e MadLoop5_resources.tar.gz && ! -e MadLoop5_resources ]]; then
+tar -xzf MadLoop5_resources.tar.gz;
+fi
+
+k=run1_app.log;
+script=ajob1;                      
+offset=$1;
+shift;
+subdir=$offset;
+for i in $@ ; do
+     j=G$i;
+     if [[ $offset == *.* ]];then
+	 subdir=${offset%.*};
+	 offset=${offset##*.};
+	 j=G${i}_${subdir};     
+     elif [[ $offset -gt 0 ]]; then
+	 j=G${i}_${subdir};
+     fi
+     if [[ ! -e $j ]]; then
+         mkdir $j;
+     fi
+     cd $j;
+     if [[ $offset -eq 0 ]]; then
+	 rm -f ftn25 ftn26 ftn99;
+	 rm -f $k;
+     else
+	 echo   "$offset"  > moffset.dat;
+     fi
+     if [[ $offset -eq $subdir ]]; then
+	 rm -f ftn25 ftn26 ftn99;
+	 rm -f $k;
+     else
+        if [[ -e ../ftn25 ]]; then
+	    cp ../ftn25 .;
+	fi
+     fi
+     if [[ ! -e input_app.txt  ]]; then
+	 cat ../input_app.txt >& input_app.txt;
+     fi
+     echo $i >> input_app.txt;
+
+     for((try=1;try<=10;try+=1)); 
+     do
+     ../madevent 2>&1 >> $k <input_app.txt | tee -a $k;
+     if [ -s $k ]
+     then
+         break;
+     else
+     sleep 1;
+     fi
+     done
+#     rm -f ftn25 ftn99
+#     rm -f ftn26
+     echo "" >> $k; echo "ls status:" >> $k; ls >> $k;
+     cp $k log.txt;
+# Perform some cleaning to keep less file on disk/transfer less file.
+     if [[ $subdir -ne 1 &&  -s results.dat && $MG5DEBUG -ne true ]]; then
+	 rm -f ftn25 &> /dev/null
+         rm -f ftn26 &> /dev/null
+	 rm -f log.txt &> /dev/null
+         rm -f *.log &> /dev/null
+	 rm -f moffset.dat &> /dev/null
+     fi
+     cd ../;
+ 
+done;
+
+# Cleaning 
+
+
+

=== modified file 'Template/LO/SubProcesses/unwgt.f'
--- Template/LO/SubProcesses/unwgt.f	2014-10-07 13:31:57 +0000
+++ Template/LO/SubProcesses/unwgt.f	2015-04-10 09:55:25 +0000
@@ -133,7 +133,7 @@
 c  Begin Code
 c-----
 c      write(*,*) 'storing Events'
-      call store_events
+      call store_events(-1d0, .True.)
       rewind(lun)
       nw = 0
       maxwgt = 0d0
@@ -213,10 +213,14 @@
       endif
       end
 
-      subroutine store_events()
+      subroutine store_events(force_max_wgt, scale_to_xsec)
 C**************************************************************************
 C     Takes events from scratch file (lun) and writes them to a permanent
 c     file  events.dat
+c     if force_max_weight =-1, then get it automatically (for a given truncation)
+c     if xscale=0 then the sum of the weight will be reweighted to the cross-section.
+c     computed from the last 3 iteration. otherwise the weight of each event
+c     will be multiply by that value.
 C**************************************************************************
       IMPLICIT NONE
 c
@@ -228,18 +232,22 @@
 c
 c     Arguments
 c
+      double precision force_max_wgt 
+      logical scale_to_xsec
 c
 c     Local
 c
       integer i, lunw, ic(7,2*nexternal-3), n, j
       logical done
       double precision wgt,p(0:4,2*nexternal-3)
-      double precision xsec,xsecabs,xerr,xscale,xtot
-      double precision xsum, xover
-      double precision target_wgt,orig_Wgt(maxevents)
+      double precision xsec,xsecabs,xerr,xtot
+      double precision xsum, xover, target_wgt
+      double precision orig_Wgt(maxevents)
+      double precision xscale
       logical store_event(maxevents)
       integer iseed, nover, nstore
       double precision scale,aqcd,aqed
+      double precision random
       integer ievent
       character*1000 buff
       logical u_syst
@@ -255,6 +263,11 @@
 
       integer                   neventswritten
       common /to_eventswritten/ neventswritten
+      
+      integer th_nunwgt
+      double precision th_maxwgt
+      common/theoretical_unwgt_max/th_maxwgt, th_nunwgt
+
 c      save neventswritten
 
       integer ngroup
@@ -273,9 +286,14 @@
 c
 c     First scale all of the events to the total cross section
 c
+
       if (nw .le. 0) return
-      call sample_result(xsecabs,xsec,xerr,itmin)
-      if (xsecabs .le. 0) return   !Fix by TS 12/3/2010
+      if (scale_to_xsec) then
+         call sample_result(xsecabs,xsec,xerr,itmin)
+         if (xsecabs .le. 0) return !Fix by TS 12/3/2010
+      else
+         xscale = nw*twgt
+      endif
       xtot=0
       call dsort(nw, swgt)
       do i=1,nw
@@ -291,12 +309,20 @@
          i = i-1
       enddo
       if (i .lt. nw) i=i+1
-      target_wgt = dabs(swgt(i))
+      th_maxwgt = dabs(swgt(i))
+      if ( force_max_wgt.lt.0)then
+         target_wgt = dabs(swgt(i))
+      else if (.not.scale_to_xsec) then
+         target_wgt = force_max_wgt / xscale
+      else
+         stop 1
+      endif
 c
 c     Select events which will be written
 c
       xsum = 0d0
       nstore = 0
+      th_nunwgt = 0
       rewind(lun)
       done = .false. 
       do i=1,nw
@@ -306,16 +332,25 @@
          else
             wgt = 0d0
          endif
-         if (dabs(wgt) .gt. target_wgt*xran1(iseed)) then
+         random = xran1(iseed)
+         if (dabs(wgt) .gt. target_wgt*random) then
             xsum=xsum+max(dabs(wgt),target_Wgt)
             store_event(i)=.true.
             nstore=nstore+1
          else
             store_event(i) = .false.
          endif
+c           we use the same seed for the two evaluation of the unweighting efficiency
+         if (dabs(wgt) .gt. th_maxwgt*random) then
+            th_nunwgt = th_nunwgt +1
+         endif
       enddo
-      xscale = xsecabs/xsum
+      if (scale_to_xsec)then
+         xscale = xsecabs/xsum
+      endif
       target_wgt = target_wgt*xscale
+      th_maxwgt = th_maxwgt*xscale
+
       rewind(lun)
 c     JA 8/17/2011 Don't check for previously stored events
 c      if (nstore .le. neventswritten) then
@@ -351,13 +386,22 @@
       enddo
       write(*,*) 'Found ',nw,' events.'
       write(*,*) 'Wrote ',i ,' events.'
-      write(*,*) 'Actual xsec ',xsec
-      write(*,*) 'Correct abs xsec ',xsecabs
-      write(*,*) 'Event xsec ', xtot
+      if (scale_to_xsec)then
+         write(*,*) 'Actual xsec ',xsec
+         write(*,*) 'Correct abs xsec ',xsecabs
+         write(*,*) 'Event xsec ', xtot
+      endif
       write(*,*) 'Events wgts > 1: ', nover
       write(*,*) '% Cross section > 1: ',xover, xover/xtot*100.
       neventswritten = i
+      maxwgt = target_wgt
+      if (force_max_wgt.lt.0)then
+         th_maxwgt = target_wgt
+         th_nunwgt = neventswritten
+      endif
+
  99   close(lunw)
+      
 c      close(lun)
       end
 

=== modified file 'Template/LO/bin/madevent'
--- Template/LO/bin/madevent	2014-06-29 23:09:00 +0000
+++ Template/LO/bin/madevent	2015-04-10 09:55:25 +0000
@@ -132,9 +132,14 @@
 try:
     if __debug__ and options.logging == 'INFO':
         options.logging = 'DEBUG'    
+    if options.logging.isdigit():
+        level = int(options.logging)
+    else:
+        level = eval('logging.' + options.logging)
+
     logging.config.fileConfig(os.path.join(root_path, 'internal', 'me5_logging.conf'))
-    logging.root.setLevel(eval('logging.' + options.logging))
-    logging.getLogger('madgraph').setLevel(eval('logging.' + options.logging))
+    logging.root.setLevel(level)
+    logging.getLogger('madgraph').setLevel(level)
 except:
     pass
 import internal.madevent_interface as cmd_interface

=== modified file 'Template/MadWeight/Cards/run_card.dat'
--- Template/MadWeight/Cards/run_card.dat	2012-03-09 09:18:24 +0000
+++ Template/MadWeight/Cards/run_card.dat	2015-04-10 09:55:25 +0000
@@ -1,8 +1,7 @@
 #*********************************************************************
-#                       MadGraph/MadEvent                            *
-#                  http://madgraph.hep.uiuc.edu                      *
+#                       MadGraph5_aMC@NLO                            *
 #                                                                    *
-#                        run_card.dat                                *
+#                     run_card.dat MadEvent                          *
 #                                                                    *
 #  This file is used to set the parameters of the run.               *
 #                                                                    *
@@ -15,183 +14,273 @@
 #
 #*******************                                                 
 # Running parameters
-#*******************                                                 
-#                                                                    
-#*********************************************************************
-# Tag name for the run (one word)                                    *
-#*********************************************************************
-  ''     = run_tag ! name of the run. overwritten by the MW card 
-#*********************************************************************
-# Run to generate the grid pack                                      *
-#*********************************************************************
-  .false.     = gridpack  !True = setting up the grid pack
-#*********************************************************************
-# Number of events and rnd seed                                      *
-#*********************************************************************
-      0       = iseed   ! rnd seed 
+#*******************
+#*********************************************************************
+# 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.    *
+#*********************************************************************
+  %(iseed)s   = iseed   ! rnd seed (0=assigned automatically=default))
 #*********************************************************************
 # Collider type and energy                                           *
+# lpp: 0=No PDF, 1=proton, -1=antiproton, 2=photon from proton,      *
+#                                         3=photon from electron     *
 #*********************************************************************
-        1     = lpp1  ! beam 1 type (0=NO PDF)
-        1     = lpp2  ! beam 2 type (0=NO PDF)
-     7000     = ebeam1  ! beam 1 energy in GeV
-     7000     = ebeam2  ! beam 2 energy in GeV
+     %(lpp1)s        = lpp1    ! beam 1 type 
+     %(lpp2)s        = lpp2    ! beam 2 type
+     %(ebeam1)s     = ebeam1  ! beam 1 total energy in GeV
+     %(ebeam2)s     = ebeam2  ! beam 2 total energy in GeV
 #*********************************************************************
 # Beam polarization from -100 (left-handed) to 100 (right-handed)    *
 #*********************************************************************
-        0     = polbeam1 ! beam polarization for beam 1
-        0     = polbeam2 ! beam polarization for beam 2
+     %(polbeam1)s     = polbeam1 ! beam polarization for beam 1
+     %(polbeam2)s     = polbeam2 ! beam polarization for beam 2
 #*********************************************************************
 # PDF CHOICE: this automatically fixes also alpha_s and its evol.    *
 #*********************************************************************
- 'cteq6l1'    = pdlabel     ! PDF set                                     
+     %(pdlabel)s    = pdlabel     ! PDF set                                     
+     %(lhaid)s    = lhaid     ! if pdlabel=lhapdf, this is the lhapdf number
 #*********************************************************************
 # Renormalization and factorization scales                           *
 #*********************************************************************
- T        = fixed_ren_scale  ! if .true. use fixed ren scale (false is beta)
- T        = fixed_fac_scale  ! if .true. use fixed fac scale (false is beta)
- 91.1880  = scale            ! fixed ren scale
- 91.1880  = dsqrt_q2fact1    ! fixed fact scale for pdf1
- 91.1880  = dsqrt_q2fact2    ! fixed fact scale for pdf2
- 1        = scalefact        ! scale factor for event-by-event scales
+ %(fixed_ren_scale)s = fixed_ren_scale  ! if .true. use fixed ren scale
+ %(fixed_fac_scale)s        = fixed_fac_scale  ! if .true. use fixed fac scale
+ %(scale)s  = scale            ! fixed ren scale
+ %(dsqrt_q2fact1)s  = dsqrt_q2fact1    ! fixed fact scale for pdf1
+ %(dsqrt_q2fact2)s  = dsqrt_q2fact2    ! fixed fact scale for pdf2
+ %(dynamical_scale_choice)s = dynamical_scale_choice ! Select one of the preselect dynamical choice 
+ %(scalefact)s  = scalefact        ! scale factor for event-by-event scales
+
+
+
+##################################################################################
+##################################################################################
+##        INFORMATION USEFULL FOR MADWEIGHT STOPS HERE                          ##
+##################################################################################
+##################################################################################
+#**********************************
+# BW cutoff (M+/-bwcutoff*Gamma)
+#**********************************
+  %(bwcutoff)s  = bwcutoff      ! (M+/-bwcutoff*Gamma)
+#*********************************************************************
+# Tag name for the run (one word)                                    *
+#*********************************************************************
+   ' '    = run_tag ! name of the run 
+  %(nevents)s = nevents ! Number of unweighted events requested 
+#*********************************************************************
+# Run to generate the grid pack                                      *
+#*********************************************************************
+  %(gridpack)s     = gridpack  !True = setting up the grid pack
 #*********************************************************************
 # Matching - Warning! ickkw > 1 is still beta
 #*********************************************************************
- 0        = ickkw            ! 0 no matching, 1 MLM, 2 CKKW matching
- 1        = highestmult      ! for ickkw=2, highest mult group
- 1        = ktscheme         ! for ickkw=1, 1 Durham kT, 2 Pythia pTE
- 1        = alpsfact         ! scale factor for QCD emission vx
- F        = chcluster        ! cluster only according to channel diag
- F        = pdfwgt           ! for ickkw=1, perform pdf reweighting
-#*********************************************************************
-#                                                                    
-#**********************************
-# BW cutoff (M+/-bwcutoff*Gamma)
-#**********************************
-  40  = bwcutoff    ! desactivate by default in the MadWeight_card
-#*******************                                                 
-# Standard Cuts    ! desactivate by default in the MadWeight_card
-#*******************                                                 
-#                                                                    
-#*********************************************************************
-# Minimum and maximum pt's                                           *
-#*********************************************************************
-  0  = ptj       ! minimum pt for the jets 
-  0  = ptb       ! minimum pt for the b 
-  0  = pta       ! minimum pt for the photons 
-  0  = ptl       ! minimum pt for the charged leptons 
-  0  = misset    ! minimum missing Et (sum of neutrino's momenta)
-  0  = ptheavy   ! minimum pt for one heavy final state
- 1.0 = ptonium   ! minimum pt for the quarkonium states
- 1d5 = ptjmax    ! maximum pt for the jets
- 1d5 = ptbmax    ! maximum pt for the b
- 1d5 = ptamax    ! maximum pt for the photons
- 1d5 = ptlmax    ! maximum pt for the charged leptons
- 1d5 = missetmax ! maximum missing Et (sum of neutrino's momenta)
-#*********************************************************************
-# Minimum and maximum E's (in the lab frame)                         *
-#*********************************************************************
-  0  = ej     ! minimum E for the jets 
-  0  = eb     ! minimum E for the b 
-  0  = ea     ! minimum E for the photons 
-  0  = el     ! minimum E for the charged leptons 
- 1d5  = ejmax ! maximum E for the jets
- 1d5  = ebmax ! maximum E for the b
- 1d5  = eamax ! maximum E for the photons
- 1d5  = elmax ! maximum E for the charged leptons
-#*********************************************************************
-# Maximum and minimum rapidity                                       *
-#*********************************************************************
- 1d2  = etaj    ! max rap for the jets 
- 1d2  = etab    ! max rap for the b 
- 1d2  = etaa    ! max rap for the photons 
- 1d2  = etal    ! max rap for the charged leptons 
- 1d2  = etaonium ! max rap for the quarkonium states
- 0d0  = etajmin ! min rap for the jets
- 0d0  = etabmin ! min rap for the b
- 0d0  = etaamin ! min rap for the photons
- 0d0  = etalmin ! main rap for the charged leptons
+ %(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
+ %(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
+#**********************************************************
+#                                                                    
+#**********************************************************
+# Apply pt/E/eta/dr/mij 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
+#*******************                                                 
+#                                                                    
+#*********************************************************************
+# Minimum and maximum pt's (for max, -1 means no cut)                *
+#*********************************************************************
+ %(ptj)s  = ptj       ! minimum pt for the jets 
+ %(ptb)s  = ptb       ! minimum pt for the b 
+ %(pta)s  = pta       ! minimum pt for the photons 
+ %(ptl)s  = ptl       ! minimum pt for the charged leptons 
+ %(misset)s  = misset    ! minimum missing Et (sum of neutrino's momenta)
+ %(ptheavy)s  = ptheavy   ! minimum pt for one heavy final state
+ %(ptjmax)s  = ptjmax    ! maximum pt for the jets
+ %(ptbmax)s  = ptbmax    ! maximum pt for the b
+ %(ptamax)s  = ptamax    ! maximum pt for the photons
+ %(ptlmax)s  = ptlmax    ! maximum pt for the charged leptons
+ %(missetmax)s  = missetmax ! maximum missing Et (sum of neutrino's momenta)
+#*********************************************************************
+# Minimum and maximum E's (in the center of mass frame)              *
+#*********************************************************************
+  %(ej)s  = ej     ! minimum E for the jets 
+  %(eb)s  = eb     ! minimum E for the b 
+  %(ea)s  = ea     ! minimum E for the photons 
+  %(el)s  = el     ! minimum E for the charged leptons 
+  %(ejmax)s   = ejmax ! maximum E for the jets
+ %(ebmax)s   = ebmax ! maximum E for the b
+ %(eamax)s   = eamax ! maximum E for the photons
+ %(elmax)s   = elmax ! maximum E for the charged leptons
+#*********************************************************************
+# Maximum and minimum absolute rapidity (for max, -1 means no cut)   *
+#*********************************************************************
+  %(etaj)s = etaj    ! max rap for the jets 
+  %(etab)s  = etab    ! max rap for the b
+ %(etaa)s  = etaa    ! max rap for the photons 
+ %(etal)s  = etal    ! max rap for the charged leptons 
+ %(etajmin)s  = etajmin ! min rap for the jets
+ %(etabmin)s  = etabmin ! min rap for the b
+ %(etaamin)s  = etaamin ! min rap for the photons
+ %(etalmin)s  = etalmin ! main rap for the charged leptons
 #*********************************************************************
 # Minimum and maximum DeltaR distance                                *
 #*********************************************************************
- 0.  = drjj    ! min distance between jets 
- 0   = drbb    ! min distance between b's 
- 0.  = drll    ! min distance between leptons 
- 0.  = draa    ! min distance between gammas 
- 0   = drbj    ! min distance between b and jet 
- 0.  = draj    ! min distance between gamma and jet 
- 0.  = drjl    ! min distance between jet and lepton 
- 0   = drab    ! min distance between gamma and b 
- 0   = drbl    ! min distance between b and lepton 
- 0.  = dral    ! min distance between gamma and lepton 
- 1d2 = drjjmax ! max distance between jets
- 1d2 = drbbmax ! max distance between b's
- 1d2 = drllmax ! max distance between leptons
- 1d2 = draamax ! max distance between gammas
- 1d2 = drbjmax ! max distance between b and jet
- 1d2 = drajmax ! max distance between gamma and jet
- 1d2 = drjlmax ! max distance between jet and lepton
- 1d2 = drabmax ! max distance between gamma and b
- 1d2 = drblmax ! max distance between b and lepton
- 1d2 = dralmax ! maxdistance between gamma and lepton
+ %(drjj)s = drjj    ! min distance between jets 
+ %(drbb)s   = drbb    ! min distance between b's 
+ %(drll)s = drll    ! min distance between leptons 
+ %(draa)s = draa    ! min distance between gammas 
+ %(drbj)s   = drbj    ! min distance between b and jet 
+ %(draj)s = draj    ! min distance between gamma and jet 
+ %(drjl)s = drjl    ! min distance between jet and lepton 
+ %(drab)s   = drab    ! min distance between gamma and b 
+ %(drbl)s   = drbl    ! min distance between b and lepton 
+ %(dral)s = dral    ! min distance between gamma and lepton 
+ %(drjjmax)s  = drjjmax ! max distance between jets
+ %(drbbmax)s  = drbbmax ! max distance between b's
+ %(drllmax)s  = drllmax ! max distance between leptons
+ %(draamax)s  = draamax ! max distance between gammas
+ %(drbjmax)s  = drbjmax ! max distance between b and jet
+ %(drajmax)s  = drajmax ! max distance between gamma and jet
+ %(drjlmax)s  = drjlmax ! max distance between jet and lepton
+ %(drabmax)s  = drabmax ! max distance between gamma and b
+ %(drblmax)s  = drblmax ! max distance between b and lepton
+ %(dralmax)s  = dralmax ! maxdistance between gamma and lepton
 #*********************************************************************
 # Minimum and maximum invariant mass for pairs                       *
+# WARNING: for four lepton final state mmll cut require to have      *
+#          different lepton masses for each flavor!                  *           
 #*********************************************************************
- 0   = mmjj    ! min invariant mass of a jet pair 
- 0   = mmbb    ! min invariant mass of a b pair 
- 0   = mmaa    ! min invariant mass of gamma gamma pair
- 0   = mmll    ! min invariant mass of l+l- (same flavour) lepton pair
- 1d5 = mmjjmax ! max invariant mass of a jet pair
- 1d5 = mmbbmax ! max invariant mass of a b pair
- 1d5 = mmaamax ! max invariant mass of gamma gamma pair
- 1d5 = mmllmax ! max invariant mass of l+l- (same flavour) lepton pair
+ %(mmjj)s   = mmjj    ! min invariant mass of a jet pair 
+ %(mmbb)s   = mmbb    ! min invariant mass of a b pair 
+ %(mmaa)s   = mmaa    ! min invariant mass of gamma gamma pair
+ %(mmll)s   = mmll    ! min invariant mass of l+l- (same flavour) lepton pair
+ %(mmjjmax)s  = mmjjmax ! max invariant mass of a jet pair
+ %(mmbbmax)s  = mmbbmax ! max invariant mass of a b pair
+ %(mmaamax)s  = mmaamax ! max invariant mass of gamma gamma pair
+ %(mmllmax)s  = mmllmax ! max invariant mass of l+l- (same flavour) lepton pair
 #*********************************************************************
 # Minimum and maximum invariant mass for all letpons                 *
 #*********************************************************************
- 0   = mmnl    ! min invariant mass for all letpons (l+- and vl) 
- 1d5 = mmnlmax ! max invariant mass for all letpons (l+- and vl) 
+ %(mmnl)s   = mmnl    ! min invariant mass for all letpons (l+- and vl) 
+ %(mmnlmax)s  = mmnlmax ! max invariant mass for all letpons (l+- and vl) 
+#*********************************************************************
+# Minimum and maximum pt for 4-momenta sum of leptons                *
+#*********************************************************************
+ %(ptllmin)s   = ptllmin  ! Minimum pt for 4-momenta sum of leptons(l and vl)
+ %(ptllmax)s  = ptllmax  ! Maximum pt for 4-momenta sum of leptons(l and vl)
 #*********************************************************************
 # Inclusive cuts                                                     *
 #*********************************************************************
- 0  = xptj ! minimum pt for at least one jet  
- 0  = xptb ! minimum pt for at least one b 
- 0  = xpta ! minimum pt for at least one photon 
- 0  = xptl ! minimum pt for at least one charged lepton 
+ %(xptj)s  = xptj ! minimum pt for at least one jet  
+ %(xptb)s  = xptb ! minimum pt for at least one b 
+ %(xpta)s  = xpta ! minimum pt for at least one photon 
+ %(xptl)s  = xptl ! minimum pt for at least one charged lepton 
 #*********************************************************************
 # Control the pt's of the jets sorted by pt                          *
 #*********************************************************************
- 0   = ptj1min ! minimum pt for the leading jet in pt
- 0   = ptj2min ! minimum pt for the second jet in pt
- 0   = ptj3min ! minimum pt for the third jet in pt
- 0   = ptj4min ! minimum pt for the fourth jet in pt
- 1d5 = ptj1max ! maximum pt for the leading jet in pt 
- 1d5 = ptj2max ! maximum pt for the second jet in pt
- 1d5 = ptj3max ! maximum pt for the third jet in pt
- 1d5 = ptj4max ! maximum pt for the fourth jet in pt
- 0   = cutuse  ! reject event if fails any (0) / all (1) jet pt cuts
+ %(ptj1min)s   = ptj1min ! minimum pt for the leading jet in pt
+ %(ptj2min)s   = ptj2min ! minimum pt for the second jet in pt
+ %(ptj3min)s   = ptj3min ! minimum pt for the third jet in pt
+ %(ptj4min)s   = ptj4min ! minimum pt for the fourth jet in pt
+ %(ptj1max)s  = ptj1max ! maximum pt for the leading jet in pt 
+ %(ptj2max)s  = ptj2max ! maximum pt for the second jet in pt
+ %(ptj3max)s  = ptj3max ! maximum pt for the third jet in pt
+ %(ptj4max)s  = ptj4max ! maximum pt for the fourth jet in pt
+ %(cutuse)s   = cutuse  ! reject event if fails any (0) / all (1) jet pt cuts
+#*********************************************************************
+# Control the pt's of leptons sorted by pt                           *
+#*********************************************************************
+ %(ptl1min)s   = ptl1min ! minimum pt for the leading lepton in pt
+ %(ptl2min)s   = ptl2min ! minimum pt for the second lepton in pt
+ %(ptl3min)s   = ptl3min ! minimum pt for the third lepton in pt
+ %(ptl4min)s   = ptl4min ! minimum pt for the fourth lepton in pt
+ %(ptl1max)s  = ptl1max ! maximum pt for the leading lepton in pt 
+ %(ptl2max)s  = ptl2max ! maximum pt for the second lepton in pt
+ %(ptl3max)s  = ptl3max ! maximum pt for the third lepton in pt
+ %(ptl4max)s  = ptl4max ! maximum pt for the fourth lepton in pt
 #*********************************************************************
 # Control the Ht(k)=Sum of k leading jets                            *
 #*********************************************************************
- 0   = htjmin ! minimum jet HT=Sum(jet pt)
- 1d5 = htjmax ! maximum jet HT=Sum(jet pt)
- 0   = ht2min ! minimum Ht for the two leading jets
- 0   = ht3min ! minimum Ht for the three leading jets
- 0   = ht4min ! minimum Ht for the four leading jets
- 1d5 = ht2max ! maximum Ht for the two leading jets
- 1d5 = ht3max ! maximum Ht for the three leading jets
- 1d5 = ht4max ! maximum Ht for the four leading jets
+ %(htjmin)s   = htjmin ! minimum jet HT=Sum(jet pt)
+ %(htjmax)s  = htjmax ! maximum jet HT=Sum(jet pt)
+ %(ihtmin)s   = ihtmin  !inclusive Ht for all partons (including b)
+ %(ihtmax)s  = ihtmax  !inclusive Ht for all partons (including b)
+ %(ht2min)s   = ht2min ! minimum Ht for the two leading jets
+ %(ht3min)s   = ht3min ! minimum Ht for the three leading jets
+ %(ht4min)s   = ht4min ! minimum Ht for the four leading jets
+ %(ht2max)s  = ht2max ! maximum Ht for the two leading jets
+ %(ht3max)s  = ht3max ! maximum Ht for the three leading jets
+ %(ht4max)s  = ht4max ! maximum Ht for the four leading jets
+#***********************************************************************
+# Photon-isolation cuts, according to hep-ph/9801442                   *
+# When ptgmin=0, all the other parameters are ignored                  *
+# When ptgmin>0, pta and draj are not going to be used                 *
+#***********************************************************************
+ %(ptgmin)s = ptgmin ! Min photon transverse momentum
+ %(r0gamma)s = R0gamma ! Radius of isolation code
+ %(xn)s = xn ! n parameter of eq.(3.4) in hep-ph/9801442
+ %(epsgamma)s = epsgamma ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
+ %(isoem)s = isoEM ! isolate photons from EM energy (photons and leptons)
 #*********************************************************************
 # WBF cuts                                                           *
 #*********************************************************************
- 0   = xetamin ! minimum rapidity for two jets in the WBF case  
- 0   = deltaeta ! minimum rapidity for two jets in the WBF case 
-#*********************************************************************
-# maximal pdg code for quark to be considered as a jet               *
-# otherwise b cuts are applied                                       *
-#*********************************************************************
- 4 = maxjetflavor
+ %(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 
+#*********************************************************************
+# 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                                                   *
 #*********************************************************************
- 0   = xqcut   ! minimum kt jet measure between partons
-#*********************************************************************
+ %(xqcut)s   = xqcut   ! minimum kt jet measure between partons
+#*********************************************************************
+#
+#*********************************************************************
+# Store info for systematics studies                                 *
+# WARNING: If use_syst is T, matched Pythia output is                *
+#          meaningful ONLY if plotted taking matchscale              *
+#          reweighting into account!                                 *
+#*********************************************************************
+   F  = use_syst      ! Enable systematics studies
+#
+#**************************************
+# Parameter of the systematics study
+#  will be used by SysCalc (if installed)
+#**************************************                                  
+#
+%(sys_scalefact)s = sys_scalefact  # factorization/renormalization scale factor
+%(sys_alpsfact)s = sys_alpsfact  # \alpha_s emission scale factors
+%(sys_matchscale)s = sys_matchscale # variation of merging scale
+# PDF sets and number of members (0 or none for all members).
+%(sys_pdf)s = sys_pdf # matching scales
+# MSTW2008nlo68cl.LHgrid 1  = sys_pdf

=== modified file 'Template/MadWeight/src/gen_ps.f'
--- Template/MadWeight/src/gen_ps.f	2014-01-06 16:54:40 +0000
+++ Template/MadWeight/src/gen_ps.f	2015-04-10 09:55:25 +0000
@@ -254,17 +254,22 @@
 c     var_num=1 means that we generate a theta
        if (var_num.eq.1) then
 
+      if (gam.gt.0) then
       point_max=c_point +5d0*gam
       if (point_max.gt.pi) point_max=pi
       point_min=c_point -5d0*gam
       if (point_min.lt.0d0) point_min=0d0
+      else
+        point_max=pi
+        point_min=0d0
+      endif
 
        gen_point=(point_max-point_min)*x+point_min
 
 c     var_num=2 means that we generate a phi (note that phi is a cyclic variable) 
       elseif(var_num.eq.2) then
 
-      if(gam.lt.(2d0*pi/10)) then
+      if(gam.lt.(2d0*pi/10).and.gam.gt.0d0) then
       point_max=c_point +5d0*gam
       point_min=c_point -5d0*gam
       else

=== modified file 'Template/MadWeight/src/run.inc'
--- Template/MadWeight/src/run.inc	2014-01-07 10:50:42 +0000
+++ Template/MadWeight/src/run.inc	2015-04-10 09:55:25 +0000
@@ -6,9 +6,9 @@
 c
       real*8          scale,scalefact,alpsfact
       logical         fixed_ren_scale,fixed_fac_scale,fixed_couplings,hmult
-      integer         ickkw,nhmult,asrwgtflavor
+      integer         ickkw,nhmult,asrwgtflavor,dynamical_scale_choice
       common/to_scale/scale,scalefact,alpsfact,fixed_ren_scale,fixed_fac_scale,
-     $                fixed_couplings,ickkw,nhmult,hmult,asrwgtflavor
+     $                fixed_couplings,ickkw,nhmult,hmult,asrwgtflavor,dynamical_scale_choice
 c
 c     Collider
 c

=== modified file 'Template/MadWeight/src/setrun.f'
--- Template/MadWeight/src/setrun.f	2014-01-07 10:50:42 +0000
+++ Template/MadWeight/src/setrun.f	2015-04-10 09:55:25 +0000
@@ -262,7 +262,7 @@
  	call get_real   (npara,param,value,"ptj3max",ptj3max,1d5)
  	call get_real   (npara,param,value,"ptj4min",ptj4min,0d0)
  	call get_real   (npara,param,value,"ptj4max",ptj4max,1d5)
-	call get_real   (npara,param,value,"cutuse",cutuse,0d0)
+	call get_integer   (npara,param,value,"cutuse",cutuse,0)
 
 c*********************************************************************
 c Check  Ht                                                          *

=== modified file 'Template/NLO/Cards/FO_analyse_card.dat'
--- Template/NLO/Cards/FO_analyse_card.dat	2013-12-02 16:53:51 +0000
+++ Template/NLO/Cards/FO_analyse_card.dat	2015-04-10 09:55:25 +0000
@@ -1,21 +1,21 @@
 #######################################################################
 #                             
-# This file contains the settings for analyses to be linked to aMC@NLO
-# fixed order runs. Analyse files are meant to be put (or linked)
-# inside <PROCDIR>/FixedOrderAnalysis/ (<PROCDIR> is the name of the
-# exported process directory). See the
-# <PROCDIR>/FixedOrderAnalysis/analysis_template.f file for details on
-# how to write your own analysis.
+# This file contains the settings for analyses to be linked to fixed
+# order runs. Analysis files are meant to be put (or linked) inside
+# <PROCDIR>/FixedOrderAnalysis/ (<PROCDIR> is the name of the exported
+# process directory). See the
+# <PROCDIR>/FixedOrderAnalysis/analysis_*_template.f file for details
+# on how to write your own analysis.
 #                                                                               
 #######################################################################
 #
-# Analysis format. Can either be 'topdrawer', 'root', or 'none'.
-# Topdrawer is human-readable text format, which allows for easy
-# conversion to other formats. When choosing topdrawer, the
-# histogramming package 'dbook.f' is included in the code, while when
-# choosing root the 'rbook_fe8.f' and 'rbook_be8.cc' are included. If
-# 'none' is chosen, all the other entries below have to be set empty.
-FO_ANALYSIS_FORMAT = topdrawer
+# Analysis format. Can either be 'topdrawer', 'root', 'HwU' or 'none'.
+# When choosing HwU, it comes with a GnuPlot wrapper. When choosing
+# topdrawer, the histogramming package 'dbook.f' is included in the
+# code, while when choosing root the 'rbook_fe8.f' and 'rbook_be8.cc'
+# are included. If 'none' is chosen, all the other entries below have
+# to be set empty.
+FO_ANALYSIS_FORMAT = HwU
 #
 # Needed extra-libraries (FastJet is already linked):
 FO_EXTRALIBS = 
@@ -31,12 +31,13 @@
 # User's analysis (to be put in the <PROCDIR>/FixedOrderAnalysis/
 # directory). Please use .o as extension and white spaces to separate
 # files.
-FO_ANALYSE = analysis_td_template.o
+FO_ANALYSE = analysis_HwU_template.o
 #
 #
 ## When linking with root, the following settings are a working
-## example on lxplus (CERN). When using this, comment out the lines
-## above and replace <PATH_TO_ROOT> with the physical path to root,
+## example on lxplus (CERN) as of July 2014. When using this, comment
+## out the lines above and replace <PATH_TO_ROOT> with the physical
+## path to root,
 ## e.g. /afs/cern.ch/sw/lcg/app/releases/ROOT/5.34.11/x86_64-slc6-gcc46-dbg/root/
 #FO_ANALYSIS_FORMAT = root
 #FO_EXTRALIBS = Core Cint Hist Matrix MathCore RIO dl Thread

=== modified file 'Template/NLO/Cards/run_card.dat'
--- Template/NLO/Cards/run_card.dat	2014-09-11 05:33:30 +0000
+++ Template/NLO/Cards/run_card.dat	2015-04-10 09:55:25 +0000
@@ -19,98 +19,104 @@
 #***********************************************************************
 # Tag name for the run (one word)                                      *
 #***********************************************************************
-  tag_1     = run_tag ! name of the run 
+  %(run_tag)s     = run_tag ! name of the run 
 #***********************************************************************
 # Number of LHE events (and their normalization) and the required      *
 # (relative) accuracy on the Xsec.                                     *
 # These values are ignored for fixed order runs                        *
 #***********************************************************************
- 10000 = nevents ! Number of unweighted events requested 
-    -1 = req_acc ! Required accuracy (-1=auto determined from nevents)
-    -1 = nevt_job! Max number of events per job in event generation. 
+ %(nevents)s = nevents ! Number of unweighted events requested 
+ %(req_acc)s = req_acc ! Required accuracy (-1=auto determined from nevents)
+ %(nevt_job)s = nevt_job! Max number of events per job in event generation. 
                  !  (-1= no split).
 #***********************************************************************
 # Normalize the weights of LHE events such that they sum or average to *
 # the total cross section                                              *
 #***********************************************************************
- average = event_norm      ! average or sum
+ %(event_norm)s = event_norm      ! average or sum
 #***********************************************************************
 # Number of points per itegration channel (ignored for aMC@NLO runs)   *
 #***********************************************************************
- 0.01   = req_acc_FO       ! Required accuracy (-1=ignored, and use the 
+ %(req_acc_fo)s   = req_acc_FO       ! Required accuracy (-1=ignored, and use the 
  	                   ! number of points and iter. below)
 # These numbers are ignored except if req_acc_FO is equal to -1
- 5000   = npoints_FO_grid  ! number of points to setup grids
- 4      = niters_FO_grid   ! number of iter. to setup grids
- 10000  = npoints_FO       ! number of points to compute Xsec
- 6      = niters_FO        ! number of iter. to compute Xsec
+ %(npoints_fo_grid)s   = npoints_FO_grid  ! number of points to setup grids
+ %(niters_fo_grid)s      = niters_FO_grid   ! number of iter. to setup grids
+ %(npoints_fo)s  = npoints_FO       ! number of points to compute Xsec
+ %(niters_fo)s      = niters_FO        ! number of iter. to compute Xsec
 #***********************************************************************
 # Random number seed                                                   *
 #***********************************************************************
-     0    = iseed       ! rnd seed (0=assigned automatically=default))
+ %(iseed)s    = iseed       ! rnd seed (0=assigned automatically=default))
 #***********************************************************************
 # Collider type and energy                                             *
 #***********************************************************************
-    1   = lpp1    ! beam 1 type (0 = no PDF)
-    1   = lpp2    ! beam 2 type (0 = no PDF)
- 6500   = ebeam1  ! beam 1 energy in GeV
- 6500   = ebeam2  ! beam 2 energy in GeV
+ %(lpp1)s   = lpp1    ! beam 1 type (0 = no PDF)
+ %(lpp2)s   = lpp2    ! beam 2 type (0 = no PDF)
+ %(ebeam1)s   = ebeam1  ! beam 1 energy in GeV
+ %(ebeam2)s   = ebeam2  ! beam 2 energy in GeV
 #***********************************************************************
 # PDF choice: this automatically fixes also alpha_s(MZ) and its evol.  *
 #***********************************************************************
- nn23nlo   = pdlabel   ! PDF set
- 244600    = lhaid     ! if pdlabel=lhapdf, this is the lhapdf number
+ %(pdlabel)s   = pdlabel   ! PDF set
+ %(lhaid)s    = lhaid     ! if pdlabel=lhapdf, this is the lhapdf number
 #***********************************************************************
 # Include the NLO Monte Carlo subtr. terms for the following parton    *
 # shower (HERWIG6 | HERWIGPP | PYTHIA6Q | PYTHIA6PT | PYTHIA8)         *
 # WARNING: PYTHIA6PT works only for processes without FSR!!!!          *
 #***********************************************************************
-  HERWIG6   = parton_shower
+  %(parton_shower)s   = parton_shower
 #***********************************************************************
 # Renormalization and factorization scales                             *
 # (Default functional form for the non-fixed scales is the sum of      *
 # the transverse masses of all final state particles and partons. This *
 # can be changed in SubProcesses/set_scales.f)                         *
 #***********************************************************************
- F        = fixed_ren_scale  ! if .true. use fixed ren scale
- F        = fixed_fac_scale  ! if .true. use fixed fac scale
- 91.188   = muR_ref_fixed    ! fixed ren reference scale 
- 91.188   = muF1_ref_fixed   ! fixed fact reference scale for pdf1
- 91.188   = muF2_ref_fixed   ! fixed fact reference scale for pdf2
+ %(fixed_ren_scale)s        = fixed_ren_scale  ! if .true. use fixed ren scale
+ %(fixed_fac_scale)s        = fixed_fac_scale  ! if .true. use fixed fac scale
+ %(mur_ref_fixed)s   = muR_ref_fixed    ! fixed ren reference scale 
+ %(muf1_ref_fixed)s   = muF1_ref_fixed   ! fixed fact reference scale for pdf1
+ %(muf2_ref_fixed)s   = muF2_ref_fixed   ! fixed fact reference scale for pdf2
+ %(dynamical_scale_choice)s = dynamical_scale_choice ! Select one of the preselect dynamical choice
 #***********************************************************************
 # Renormalization and factorization scales (advanced and NLO options)  *
 #***********************************************************************
- F        = fixed_QES_scale  ! if .true. use fixed Ellis-Sexton scale
- 91.188   = QES_ref_fixed    ! fixed Ellis-Sexton reference scale
- 1        = muR_over_ref     ! ratio of current muR over reference muR
- 1        = muF1_over_ref    ! ratio of current muF1 over reference muF1
- 1        = muF2_over_ref    ! ratio of current muF2 over reference muF2
- 1        = QES_over_ref     ! ratio of current QES over reference QES
+ %(fixed_qes_scale)s = fixed_QES_scale  ! if .true. use fixed Ellis-Sexton scale
+ %(qes_ref_fixed)s = QES_ref_fixed    ! fixed Ellis-Sexton reference scale
+ %(mur_over_ref)s  = muR_over_ref     ! ratio of current muR over reference muR
+ %(muf1_over_ref)s = muF1_over_ref    ! ratio of current muF1 over reference muF1
+ %(muf2_over_ref)s = muF2_over_ref    ! ratio of current muF2 over reference muF2
+ %(qes_over_ref)s  = QES_over_ref     ! ratio of current QES over reference QES
 #*********************************************************************** 
 # Reweight flags to get scale dependence and PDF uncertainty           *
 # For scale dependence: factor rw_scale_up/down around central scale   *
 # For PDF uncertainty: use LHAPDF with supported set                   *
 #***********************************************************************
- .true.   = reweight_scale   ! reweight to get scale dependence
-  0.5     = rw_Rscale_down   ! lower bound for ren scale variations
-  2.0     = rw_Rscale_up     ! upper bound for ren scale variations
-  0.5     = rw_Fscale_down   ! lower bound for fact scale variations
-  2.0     = rw_Fscale_up     ! upper bound for fact scale variations
- .false.  = reweight_PDF     ! reweight to get PDF uncertainty
- 244601   = PDF_set_min      ! First of the error PDF sets 
- 244700   = PDF_set_max      ! Last of the error PDF sets
-#***********************************************************************
-# Merging - WARNING! Applies merging only at the hard-event level.     *
-# After showering an MLM-type merging should be applied as well.       *
-# See http://amcatnlo.cern.ch/FxFx_merging.htm for more details.       *
-#***********************************************************************
- 0        = ickkw            ! 0 no merging, 3 FxFx merging, 4 UNLOPS
+ %(reweight_scale)s   = reweight_scale   ! reweight to get scale dependence
+ %(rw_rscale_down)s   = rw_Rscale_down   ! lower bound for ren scale variations
+ %(rw_rscale_up)s     = rw_Rscale_up     ! upper bound for ren scale variations
+ %(rw_fscale_down)s     = rw_Fscale_down   ! lower bound for fact scale variations
+ %(rw_fscale_up)s     = rw_Fscale_up     ! upper bound for fact scale variations
+ %(reweight_pdf)s  = reweight_PDF     ! reweight to get PDF uncertainty
+ %(pdf_set_min)s   = PDF_set_min      ! First of the error PDF sets 
+ %(pdf_set_max)s   = PDF_set_max      ! Last of the error PDF sets
+#***********************************************************************
+# ickkw parameter:                                                     *
+#   0: No merging                                                      *
+#   3: FxFx Merging - WARNING! Applies merging only at the hard-event  *
+#      level. After showering an MLM-type merging should be applied as *
+#      well. See http://amcatnlo.cern.ch/FxFx_merging.htm for details. *
+#   4: UNLOPS merging (with pythia8 only). No interface from within    *
+#      MG5_aMC available, but available in Pythia8.                    *
+#  -1: NNLL+NLO jet-veto computation. See arxiv:1412.8408 [hep-ph].    *
+#***********************************************************************
+ %(ickkw)s        = ickkw
 #***********************************************************************
 #
 #***********************************************************************
 # BW cutoff (M+/-bwcutoff*Gamma)                                       *
 #***********************************************************************
- 15  = bwcutoff
+ %(bwcutoff)s  = bwcutoff
 #***********************************************************************
 # Cuts on the jets                                                     *
 # Jet clustering is performed by FastJet.
@@ -118,37 +124,37 @@
 # considerably softer than the analysis cuts.                          *
 # (more specific cuts can be specified in SubProcesses/cuts.f)         *
 #***********************************************************************
-   1  = jetalgo   ! FastJet jet algorithm (1=kT, 0=C/A, -1=anti-kT)
- 0.7  = jetradius ! The radius parameter for the jet algorithm
-  10  = ptj       ! Min jet transverse momentum
-  -1  = etaj      ! Max jet abs(pseudo-rap) (a value .lt.0 means no cut)
+ %(jetalgo)s  = jetalgo   ! FastJet jet algorithm (1=kT, 0=C/A, -1=anti-kT)
+ %(jetradius)s  = jetradius ! The radius parameter for the jet algorithm
+ %(ptj)s  = ptj       ! Min jet transverse momentum
+ %(etaj)s  = etaj      ! Max jet abs(pseudo-rap) (a value .lt.0 means no cut)
 #***********************************************************************
 # Cuts on the charged leptons (e+, e-, mu+, mu-, tau+ and tau-)        *
 # (more specific gen cuts can be specified in SubProcesses/cuts.f)     *
 #***********************************************************************
-   0  = ptl     ! Min lepton transverse momentum
-  -1  = etal    ! Max lepton abs(pseudo-rap) (a value .lt.0 means no cut)
-   0  = drll    ! Min distance between opposite sign lepton pairs
-   0  = drll_sf ! Min distance between opp. sign same-flavor lepton pairs
-   0  = mll     ! Min inv. mass of all opposite sign lepton pairs
-  30  = mll_sf  ! Min inv. mass of all opp. sign same-flavor lepton pairs
+  %(ptl)s  = ptl     ! Min lepton transverse momentum
+  %(etal)s  = etal    ! Max lepton abs(pseudo-rap) (a value .lt.0 means no cut)
+  %(drll)s  = drll    ! Min distance between opposite sign lepton pairs
+  %(drll_sf)s  = drll_sf ! Min distance between opp. sign same-flavor lepton pairs
+  %(mll)s  = mll     ! Min inv. mass of all opposite sign lepton pairs
+  %(mll_sf)s  = mll_sf  ! Min inv. mass of all opp. sign same-flavor lepton pairs
 #***********************************************************************
 # Photon-isolation cuts, according to hep-ph/9801442                   *
 # When ptgmin=0, all the other parameters are ignored                  *
 #***********************************************************************
-  20  = ptgmin    ! Min photon transverse momentum
-  -1  = etagamma  ! Max photon abs(pseudo-rap)
- 0.4  = R0gamma   ! Radius of isolation code
- 1.0  = xn        ! n parameter of eq.(3.4) in hep-ph/9801442
- 1.0  = epsgamma  ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
- .true.  = isoEM  ! isolate photons from EM energy (photons and leptons)
+  %(ptgmin)s  = ptgmin    ! Min photon transverse momentum
+  %(etagamma)s  = etagamma  ! Max photon abs(pseudo-rap)
+  %(r0gamma)s  = R0gamma   ! Radius of isolation code
+  %(xn)s  = xn        ! n parameter of eq.(3.4) in hep-ph/9801442
+  %(epsgamma)s  = epsgamma  ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
+  %(isoem)s  = isoEM  ! isolate photons from EM energy (photons and leptons)
 #***********************************************************************
 # Maximal PDG code for quark to be considered a jet when applying cuts.*
 # At least all massless quarks of the model should be included here.   *
 #***********************************************************************
- 4 = maxjetflavor
+ %(maxjetflavor)s = maxjetflavor
 #***********************************************************************
 # For aMCfast+APPLGRID use in PDF fitting (http://amcfast.hepforge.org)*
 #***********************************************************************
- 0 = iappl   ! aMCfast switch (0=OFF, 1=prepare APPLgrids, 2=fill grids)
+ %(iappl)s = iappl   ! aMCfast switch (0=OFF, 1=prepare APPLgrids, 2=fill grids)
 #***********************************************************************

=== modified file 'Template/NLO/Cards/shower_card.dat'
--- Template/NLO/Cards/shower_card.dat	2015-01-13 15:36:47 +0000
+++ Template/NLO/Cards/shower_card.dat	2015-04-10 09:55:25 +0000
@@ -21,7 +21,7 @@
 #***********************************************************************
 nevents      = -1     # N evts to shower (< 0 = all)
 nsplit_jobs  = 1      # N jobs to run in parallel (< 100!!)
-combine_td   = T      # combine the topdrawer files if nsplit_jobs > 1
+combine_td   = T      # combine the topdrawer/HwU files if nsplit_jobs>1
 maxprint     = 2      # N evts to print in the log
 maxerrs      = 0.1    # max fraction of errors
 rnd_seed     = 0      # 1st random seed (0 = default)
@@ -100,6 +100,8 @@
 # create a StdHEP/HepMC file, but to directly run an own analysis (to  *
 # be placed in HWAnalyzer or analogous MCatNLO subfolders).            *
 # Please use files in those folders as examples.                       *
+# Remember that if your analysis uses hbook or is in the HwU format,   *
+# you must also add to hbook.o or HwU.o to the ANALYSE list as well.   *
 #***********************************************************************
 EXTRALIBS    = stdhep Fmcfio     # Extra-libraries (not LHAPDF) 
                                  # Default: "stdhep Fmcfio"

=== added file 'Template/NLO/FixedOrderAnalysis/HwU.f'
--- Template/NLO/FixedOrderAnalysis/HwU.f	1970-01-01 00:00:00 +0000
+++ Template/NLO/FixedOrderAnalysis/HwU.f	2015-04-10 09:55:25 +0000
@@ -0,0 +1,410 @@
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+C                                                                          C
+C                  HwU: Histograms with Uncertainties                      C
+C                   By Rikkert Frederix, Dec. 2014                         C
+C                                                                          C
+C     Book, fill and write out histograms. Particularly suited for NLO     C
+C     computations with correlations between points (ie. event and         C
+C     counter-event) and multiple weights for given points (e.g. scale     C
+C     and PDF uncertainties through reweighting).                          C
+C                                                                          C
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+      
+c To be called once at the start of each run. Initialises the packages
+c and sets the number of weights that need to be included for each point.
+      subroutine HwU_inithist(nweights,wgt_info)
+      implicit none
+      include "HwU.inc"
+      integer i,nweights
+      character*(*) wgt_info(*)
+      do i=1,max_plots
+         booked(i)=.false.
+      enddo
+c     Number of weights associated to each point. Note that the first
+c     weight should always be the 'central value' and it should not be
+c     zero if any of the other weights are non-zero.
+      nwgts=nweights
+      if (nwgts.gt.max_wgts) then
+         write (*,*) 'ERROR: increase max_wgts in HwU histogramming'
+     $        ,max_wgts,nwgts
+         stop 1
+      endif
+      do i=1,nwgts
+         wgts_info(i)=wgt_info(i)
+      enddo
+      return
+      end
+
+      subroutine set_error_estimation(input)
+c Error estimation
+c     input 0: Simply sum all PS points; error is estimated by assuming
+c        Binomial statistics: use for unweighted events. (Assumes that
+c        all events have equal weight (up to a sign)).
+c     input 1: Simply sum all PS points; error is estimated by the
+c        variance between the separate iterations. This assumes that
+c        each iteration has the same number of PS points (This is the
+c        topdrawer default in MG5_aMC).
+c     input 2: Sum PS points for a given iteration and error estimate by
+c       square root of the sum of the squares. Perform a weighted average
+c       iteration-by-iteration
+      implicit none
+      integer input
+      integer error_estimation
+      common /HwU_common2/ error_estimation
+      if (input.ge.0 .and. input.le.2) then
+         error_estimation=input
+      else
+         write (*,*) 'unknown error estimation',input
+         stop 1
+      endif
+      return
+      end
+      
+      block data HwU
+c set the default for the error estimation method. To reduce the size of
+c the executable, put the error_estimation variable in a separate common
+c block. If we would have included the 'HwU.inc' file here, that
+c complete common block seems to be included in the size executable
+c (approx. 110 MB).
+      integer error_estimation
+      common /HwU_common2/ error_estimation
+      data error_estimation /1/
+      end
+      
+c Book the histograms at the start of the run. Give a 'label' (an
+c integer) that identifies the plot when filling it and a title
+c ('title_l') for each plot. Also the number of bins ('nbin_l') and the
+c plot range (from 'xmin' to 'xmax') should be given.
+      subroutine HwU_book(label,title_l,nbin_l,xmin,xmax)
+      implicit none
+      include "HwU.inc"
+      integer label,nbin_l,i,j
+      character*(*) title_l
+      double precision xmin,xmax
+c     Check that label and number of bins are reasonable
+      if (label.gt.max_plots) then
+         write (*,*) 'ERROR: increase max_plots in HwU histogramming'
+     $        ,max_plots, label
+         stop 1
+      endif
+      if (nbin_l.gt.max_bins) then
+         write (*,*) 'ERROR: increase max_bins in HwU histogramming'
+     $        ,max_bins,nbin_l
+         stop 1
+      endif
+      booked(label)=.true.
+      title(label)=title_l
+      nbin(label)=nbin_l
+c     Compute the bin width
+      step(label)=(xmax-xmin)/dble(nbin(label))
+      do i=1,nbin(label)
+c     Compute the lower and upper bin edges
+         histxl(label,i)=xmin+step(label)*dble(i-1)
+         histxm(label,i)=xmin+step(label)*dble(i)
+c     Set all the bins to zero.
+         do j=1,nwgts
+            histy(j,label,i)=0d0
+            histy_acc(j,label,i)=0d0
+         enddo
+         histi(label,i)=0
+         histy2(label,i)=0d0
+         histy_err(label,i)=0d0
+      enddo
+      np=0
+      return
+      end
+      
+c Fill the histograms identified with 'label' with a point at x with
+c weights giving by 'wgts'. The dimension of the 'wgts' array should
+c always be as large as the number of weights ('nweights') specified in
+c the 'HwU_inithist' subroutine. That means that each point should have
+c the same number of weights.
+      subroutine HwU_fill(label,x,wgts)
+      implicit none
+      include "HwU.inc"
+      integer label,i,j,bin
+      double precision x, wgts(*)
+c     If central weight is zero do not add this point.
+      if (wgts(1).eq.0d0) return
+c     Check if point is within plotting range
+      if (x.lt.histxl(label,1) .or.
+     $     x.gt.histxm(label,nbin(label))) return
+c     Compute the bin to fill
+      bin=int((x-histxl(label,1))/step(label)) +1
+c     To prevent numerical inaccuracy, double check that bin is
+c     reasonable
+      if (bin.lt.1 .or. bin.gt.nbin(label)) return
+c     Check if we already had another (correlated) point in this bin. If
+c     so, add the current weight to that point.
+      do i=1,np
+         if ( p_label(i).eq.label .and. p_bin(i).eq.bin) then
+            do j=1,nwgts
+               p_wgts(j,i)=p_wgts(j,i)+wgts(j)
+            enddo
+            return
+         endif
+      enddo
+c     If a new bin, add it to the list of points
+      np=np+1
+      if (np.gt.max_points) then
+         write (*,*) 'ERROR: increase max_points in HwU histogramming'
+     $        ,max_points
+         stop 1
+      endif
+      p_label(np)=label
+      p_bin(np)=bin
+      do j=1,nwgts
+         p_wgts(j,np)=wgts(j)
+      enddo
+      return
+      end
+
+c Call after all correlated contributions for a give phase-space
+c point. I.e., every time you get a new set of random numbers from
+c MINT/VEGAS. It adds the current list of points to the histograms. Add
+c the squares to compute the statistical uncertainty on the bin. Do the
+c second only for the weight corresponding to the 'central value'. In
+c this way, correlations between events and counter-events can be
+c correctly taken into account.
+      subroutine HwU_add_points
+      implicit none
+      include "HwU.inc"
+      integer i,j
+      do i=1,np
+         do j=1,nwgts
+            histy(j,p_label(i),p_bin(i))=
+     $           histy(j,p_label(i),p_bin(i))+p_wgts(j,i)
+         enddo
+         histi(p_label(i),p_bin(i))=histi(p_label(i),p_bin(i))+1
+         histy2(p_label(i),p_bin(i))=
+     $        histy2(p_label(i),p_bin(i))+p_wgts(1,i)**2
+      enddo
+      np=0
+      return
+      end
+
+
+c *****For plotting during fixed order f(N)LO runs*****
+c Call after every iteration. This adds the histograms of the current
+c iteration ('histy') to the accumulated results ('histy_acc'), with the
+c uncertainty estimate given in 'histy_err' and empties the arrays for
+c the current iteration so that they can be filled with the next
+c iteration.
+      subroutine HwU_accum_iter(inclde,nPSpoints)
+      implicit none
+      include "HwU.inc"
+      logical inclde
+      integer nPSpoints,label,i,j
+      double precision nPSinv,etot,vtot(max_wgts),niter,y_squared
+      data niter /0d0/
+      nPSinv = 1d0/dble(nPSpoints)
+      if (inclde) niter = niter+1d0
+      do label=1,max_plots
+         if (.not. booked(label)) cycle
+         if (inclde) then
+            call accumulate_results(label,nPSinv,niter)
+         endif
+c     Reset the histo of the current iteration to zero so that they are
+c     ready for the next iteration.
+         do i=1,nbin(label)
+            do j=1,nwgts
+               histy(j,label,i)=0d0
+            enddo
+            histy2(label,i)=0d0
+            histi(label,i)=0
+         enddo
+      enddo
+      return
+      end
+
+c *****For plotting unweighted events****
+c Overwrites the accumulated results ('histy_acc') with the current
+c results ('histy') and provides an uncertainty estimate in
+c 'histy_err'. This means that during the plotting of events, at
+c intermediate stages this function can be called (together with
+c HwU_output) to write intermediate plots to disk.
+      subroutine finalize_histograms(nPSpoints)
+      implicit none
+      include "HwU.inc"
+      integer label,nPSpoints,i,j
+      double precision nPSinv,niter
+      nPSinv=1d0/dble(nPSpoints)
+      niter=1d0
+      do label=1,max_plots
+         if (.not. booked(label)) cycle
+         do i=1,nbin(label)
+c     Set all the bins to zero.
+            do j=1,nwgts
+               histy_acc(j,label,i)=0d0
+            enddo
+            histy_err(label,i)=0d0
+         enddo
+         call accumulate_results(label,nPSinv,niter)
+      enddo
+      return
+      end
+
+c This adds the histograms of the current iteration ('histy') to the
+c accumulated results ('histy_acc'), with the uncertainty estimate given
+c in 'histy_err'. When error_estimation==2 it adds using a weight
+c corresponding to the statistical uncertainties on the 'central value'
+c weights ('etot' (from 'histy2') and 'histy_err', respectively), if
+c error_estimation==1 it simply averages over the iterations, if
+c error_estimation==0 it uses Poisson statistics to compute the
+c uncertainty. Note that this means that during the filling of the
+c histograms the central value should not be zero if any of the other
+c weights are non-zero.
+      subroutine accumulate_results(label,nPSinv,niter)
+      implicit none
+      include "HwU.inc"
+      integer label,i,j
+      double precision nPSinv,etot,vtot(max_wgts),niter,y_squared
+      integer error_estimation
+      common /HwU_common2/ error_estimation
+      if (error_estimation.eq.2) then
+c     Use the weighted average bin-by-bin. This is not really justified
+c     for fNLO computations, because for bins with low statistics, the
+c     variance computation cannot be trusted.
+         do i=1,nbin(label)
+c     Skip bin if no entries
+            if (histi(label,i).eq.0) cycle
+c     Divide weights by the number of PS points. This means that this is
+c     now normalised to the total cross section in that bin
+            do j=1,nwgts
+               vtot(j)=histy(j,label,i)*nPSinv
+            enddo
+c     Error estimation of the current bin
+            etot=sqrt(abs(histy2(label,i)*nPSinv-vtot(1)**2)*nPSinv)
+c     Include "Bessel's correction" to have a corrected (even though
+c     still biased) estimator of the standard deviation.
+            if (histi(label,i).gt.1) then
+               etot=etot* sqrt(dble(histi(label,i))
+     &                    /(dble(histi(label,i))-1.5d0))
+            else
+               etot=0.999d49
+            endif
+c     If the error estimation of the accumulated results is still zero
+c     (i.e. no points were added yet, e.g. because it is the first
+c     iteration) simply copy the results of this iteration over the
+c     accumulated results.
+            if (histy_err(label,i).eq.0d0) then
+               do j=1,nwgts
+                  histy_acc(j,label,i)=vtot(j)
+               enddo
+               histy_err(label,i)=etot
+            else
+c     Add the results of the current iteration to the accumulated results
+               do j=1,nwgts
+                  histy_acc(j,label,i)=(histy_acc(j,label,i)
+     &                 /histy_err(label,i)+vtot(j)/etot)/(1d0
+     &                 /histy_err(label,i) + 1d0/etot)
+               enddo
+               histy_err(label,i)=
+     $              1d0/sqrt(1d0/histy_err(label,i)**2+1d0/etot**2)
+            endif
+         enddo
+      elseif(error_estimation.eq.1) then
+c     simply sum the weights in the bins
+         do i=1,nbin(label)
+            if (histi(label,i).eq.0 .and.
+     &           histy_acc(1,label,i).eq.0d0) cycle
+            if (niter.ne.1d0) y_squared=((niter-1)*histy_err(label,i))
+     &           **2+(niter-1)*histy_acc(1,label,i)**2
+            do j=1,nwgts
+               vtot(j)=histy(j,label,i)*nPSinv
+               histy_acc(j,label,i)=(histy_acc(j,label,i)*(niter-1d0)
+     &              +vtot(j))/niter
+            enddo
+c     base the error on the variance in the results per iteration. For a
+c     small number of iterations, this underestimates the actual
+c     uncertainty.
+            if (niter.eq.1d0) then
+               histy_err(label,i)=0d0
+            else
+               histy_err(label,i)= sqrt(((y_squared+vtot(1)**2)/niter
+     &              -histy_acc(1,label,i)**2)/niter)
+            endif
+         enddo
+      elseif(error_estimation.eq.0) then
+c     simply sum the weights in the bins
+         do i=1,nbin(label)
+            if (histi(label,i).eq.0 .and.
+     &           histy_acc(1,label,i).eq.0d0) cycle
+            do j=1,nwgts
+               vtot(j)=histy(j,label,i)*nPSinv
+               histy_acc(j,label,i)=(histy_acc(j,label,i)*(niter-1d0)
+     &              +vtot(j))/niter
+            enddo
+c     Error estimation of the current bin using Poisson statistics,
+c     making sure we normalise with the number of iterations.
+            if (histi(label,i).ne.0) then
+               etot=sqrt(histy2(label,i))*nPSinv
+               histy_err(label,i)= sqrt(((niter-1d0)*histy_err(label,i))
+     &              **2+etot**2)/niter
+            else
+               histy_err(label,i)=(niter-1d0)/niter*histy_err(label,i)
+            endif
+         enddo
+      endif
+      end
+      
+c Write the histograms to disk at the end of the run, multiplying the
+c output by 'xnorm'
+      subroutine HwU_output(unit,xnorm)
+      implicit none
+      include "HwU.inc"
+      integer unit,i,j,label
+      integer max_length
+      parameter (max_length=(max_wgts+3)*17)
+      character*(max_length) buffer
+      character*4 str_nbin
+      double precision xnorm
+c     column info: x_min, x_max, y (central value), dy, {extra
+c     weights}. Use columns with a width of 17 characters.
+      write (buffer( 1:17),'(a)')'##&       xmin  &'
+      write (buffer(18:34),'(a)')'            xmax '
+      write (buffer(35:51),'(a2,a15)') ' &',wgts_info(1)(1:15)
+      write (buffer(52:68),'(a)')' &            dy '
+      do j=2,nwgts
+         write (buffer((j+2)*17+1:(j+3)*17),'(a2,a15)')
+     $        ' &',wgts_info(j)(1:15)
+      enddo
+      write (unit,'(a)') buffer(1:(nwgts+3)*17)
+      write (unit,'(a)') ''
+      do label=1,max_plots
+         if (.not. booked(label)) cycle
+c     title
+c        For some weird reason, it is no possible to include directly
+c        the integer in two line below with the format '(12a,i4.4,3a,a,2a)'
+c        this is why I have to do it in two steps instead.
+         write (str_nbin,'(i4.4)') nbin(label)
+         write (unit,'(12a,4a,3a,a,2a)') '<histogram> ',str_nbin,
+     &                                       ' " ',title(label),' "'
+c     data
+         do i=1,nbin(label)
+           write (buffer( 1:16),'(2x,e14.7)') histxl(label,i)
+           write (buffer(17:32),'(2x,e14.7)') histxm(label,i)
+           write (buffer(33:48),'(2x,e14.7)') histy_acc(1,label,i)*xnorm
+           write (buffer(49:64),'(2x,e14.7)') histy_err(label,i)*xnorm
+           do j=2,nwgts
+              write (buffer((j+2)*16+1:(j+3)*16),'(2x,e14.7)')
+     $             histy_acc(j,label,i)*xnorm
+           enddo
+           write (unit,'(a)') buffer(1:(nwgts+3)*16)
+         enddo
+c     2 empty lines after each plot
+         write (unit,'(12a)') '<\histogram>'
+         write (unit,'(a)') ''
+         write (unit,'(a)') ''         
+      enddo
+      return
+      end
+
+c dummy subroutine
+      subroutine accum(idummy)
+      integer idummy
+      end
+c dummy subroutine
+      subroutine addfil(string)
+      character*(*) string
+      end

=== added file 'Template/NLO/FixedOrderAnalysis/HwU.inc'
--- Template/NLO/FixedOrderAnalysis/HwU.inc	1970-01-01 00:00:00 +0000
+++ Template/NLO/FixedOrderAnalysis/HwU.inc	2015-04-10 09:55:25 +0000
@@ -0,0 +1,22 @@
+*     -*-fortran-*-
+
+      integer max_plots,max_bins,max_wgts,max_points
+      parameter (max_plots=200)
+      parameter (max_bins=100)
+      parameter (max_wgts=300)
+      parameter (max_points=max_plots*40)
+
+      logical booked(max_plots)
+      integer nbin(max_plots),nwgts,np,p_bin(max_points)
+     &     ,p_label(max_points),histi(max_plots,max_bins)
+      character*50 title(max_plots)
+      character*15 wgts_info(max_wgts)
+      double precision histy(max_wgts,max_plots,max_bins)
+     $     ,histy_acc(max_wgts,max_plots,max_bins),histy2(max_plots
+     $     ,max_bins),histy_err(max_plots,max_bins),histxl(max_plots
+     $     ,max_bins),histxm(max_plots,max_bins),step(max_plots)
+     $     ,p_wgts(max_wgts,max_points)
+
+      common/HwU_common/histy,histy_acc,histy2,histy_err,histxl,histxm
+     &     ,p_wgts,step,histi,nbin,p_bin,p_label,np,nwgts
+     &     ,booked,title,wgts_info

=== added file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_V.f'
--- Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_V.f	1970-01-01 00:00:00 +0000
+++ Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_V.f	2015-04-10 09:55:25 +0000
@@ -0,0 +1,139 @@
+c
+c Example analysis for "p p > w+ [QCD]" process.
+c Example analysis for "p p > w- [QCD]" process.
+c Example analysis for "p p > z [QCD]" process.
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_begin(nwgt,weights_info)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      integer nwgt
+      character*(*) weights_info(*)
+      integer j,kk,l
+      character*6 cc(2)
+      data cc/'|T@NLO','|T@LO '/
+      real * 8 xmi,xms
+      call HwU_inithist(nwgt,weights_info)
+      xmi=40.d0
+      xms=120.d0
+      do j=1,2
+         l=(j-1)*5
+         call HwU_book(l+ 1,'V pt     '//cc(j),100,0.d0,200.d0)
+         call HwU_book(l+ 2,'V log pt '//cc(j),100,0.d0,5.d0)
+         call HwU_book(l+ 3,'V y      '//cc(j), 78,-9.d0,9.d0)
+         call HwU_book(l+ 4,'V eta    '//cc(j), 78,-9.d0,9.d0)
+         call HwU_book(l+ 5,'mV       '//cc(j), 80,xmi,xms)
+      enddo
+      return
+      end
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_end(dummy)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      double precision dummy
+      call HwU_write_file
+      return                
+      end
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'nexternal.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision p(0:4,nexternal)
+      double precision wgts(*)
+      integer ibody
+      double precision wgt,var
+      integer i,kk,l
+      double precision www,pv(0:3),xmv,ptv,yv,etav
+      double precision getrapidity,getpseudorap,dot
+      external getrapidity,getpseudorap,dot
+      if (nexternal.ne.4) then
+         write (*,*) 'error #1 in analysis_fill: '/
+     &        /'only for process "p p > V [QCD]"'
+         stop 1
+      endif
+      if (.not. (abs(ipdg(1)).le.5 .or. ipdg(1).eq.21)) then
+         write (*,*) 'error #2 in analysis_fill: '/
+     &        /'only for process "p p > V [QCD]"'
+         stop 1
+      endif
+      if (.not. (abs(ipdg(2)).le.5 .or. ipdg(2).eq.21)) then
+         write (*,*) 'error #3 in analysis_fill: '/
+     &        /'only for process "p p > V [QCD]"'
+         stop 1
+      endif
+      if (.not. (abs(ipdg(4)).le.5 .or. ipdg(4).eq.21)) then
+         write (*,*) 'error #4 in analysis_fill: '/
+     &        /'only for process "p p > V [QCD]"'
+         stop 1
+      endif
+      if (abs(ipdg(3)).ne.24.and.ipdg(3).ne.23) then
+         write (*,*) 'error #5 in analysis_fill: '/
+     &        /'only for process "p p > V [QCD]"'
+         stop 1
+      endif
+C
+      do i=0,3
+        pv(i)=p(i,3)
+      enddo
+      xmv=sqrt(max(dot(pv,pv),0d0))
+      ptv=sqrt(max(pv(1)**2+pv(2)**2,0d0))
+      yv=getrapidity(pv(0),pv(3))
+      etav=getpseudorap(pv(0),pv(1),pv(2),pv(3))
+C
+      do i=1,2
+         l=(i-1)*5
+         if (ibody.ne.3 .and.i.eq.2) cycle
+         call HwU_fill(l+1,ptv,WGTS)
+         if(ptv.gt.0) call HwU_fill(l+2,log10(ptv),WGTS)
+         call HwU_fill(l+3,yv,WGTS)
+         call HwU_fill(l+4,etav,WGTS)
+         call HwU_fill(l+5,xmv,WGTS)
+      enddo
+C
+ 999  return      
+      end
+
+
+
+      function getrapidity(en,pl)
+      implicit none
+      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
+      parameter (tiny=1.d-8)
+      xplus=en+pl
+      xminus=en-pl
+      if(xplus.gt.tiny.and.xminus.gt.tiny)then
+         if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny)then
+            y=0.5d0*log( xplus/xminus  )
+         else
+            y=sign(1.d0,pl)*1.d8
+         endif
+      else 
+         y=sign(1.d0,pl)*1.d8
+      endif
+      getrapidity=y
+      return
+      end
+
+
+      function getpseudorap(en,ptx,pty,pl)
+      implicit none
+      real*8 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
+      getpseudorap=eta
+      return
+      end

=== added file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_h.f'
--- Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_h.f	1970-01-01 00:00:00 +0000
+++ Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_h.f	2015-04-10 09:55:25 +0000
@@ -0,0 +1,234 @@
+c
+c Example analysis for "p p > h [QCD]" process.
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_begin(nwgt,weights_info)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      integer nwgt
+      character*(*) weights_info(*)
+      integer kk,l
+      call HwU_inithist(nwgt,weights_info)
+      call HwU_book(l+1,'Higgs pT ' ,100,0.d0,200.d0)
+      call HwU_book(l+2,'Higgs pT ' ,100,0.d0,500.d0)
+      call HwU_book(l+3,'Higgs log[pT] ' ,98,0.1d0,5.d0)
+      call HwU_book(l+4,'Higgs pT,|y_H|<2 ' ,100,0.d0,200.d0)
+      call HwU_book(l+5,'Higgs pT,|y_H|<2 ' ,100,0.d0,500.d0)
+      call HwU_book(l+6,'Higgs log[pT],|y_H|<2 ' ,98,0.1d0,5.d0)
+
+      call HwU_book(l+7,'j1 pT ' ,100,0.d0,200.d0)
+      call HwU_book(l+8,'j1 pT ' ,100,0.d0,500.d0)
+      call HwU_book(l+9,'j1 log[pT] ' ,98,0.1d0,5.d0)
+      call HwU_book(l+10,'j1 pT,|y_j1|<2 ' ,100,0.d0,200.d0)
+      call HwU_book(l+11,'j1 pT,|y_j1|<2 ' ,100,0.d0,500.d0)
+      call HwU_book(l+12,'j1 log[pT],|y_j1|<2 ' ,98,0.1d0,5.d0)
+
+      call HwU_book(l+13,'Inc j pT ' ,100,0.d0,200.d0)
+      call HwU_book(l+14,'Inc j pT ' ,100,0.d0,500.d0)
+      call HwU_book(l+15,'Inc j log[pT] ' ,98,0.1d0,5.d0)
+      call HwU_book(l+16,'Inc j pT,|y_Ij|<2 ' ,100,0.d0,2.d2)
+      call HwU_book(l+17,'Inc j pT,|y_Ij|<2 ' ,100,0.d0,5.d2)
+      call HwU_book(l+18,'Inc j log[pT],|y_Ij|<2' ,98,0.1d0,5.d0)
+
+      call HwU_book(l+19,'Higgs y ' ,60,-6.d0,6.d0)
+      call HwU_book(l+20,'Higgs y,pT_H>10GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+21,'Higgs y,pT_H>30GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+22,'Higgs y,pT_H>50GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+23,'Higgs y,pT_H>70GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+24,'Higgs y,pt_H>90GeV ' ,100,-6.d0,6.d0)
+
+      call HwU_book(l+25,'j1 y ' ,60,-6.d0,6.d0)
+      call HwU_book(l+26,'j1 y,pT_j1>10GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+27,'j1 y,pT_j1>30GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+28,'j1 y,pT_j1>50GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+29,'j1 y,pT_j1>70GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+30,'j1 y,pT_j1>90GeV ' ,100,-6.d0,6.d0)
+
+      call HwU_book(l+31,'H-j1 y ' ,60,-6.d0,6.d0)
+      call HwU_book(l+32,'H-j1 y,pT_j1>10GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+33,'H-j1 y,pT_j1>30GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+34,'H-j1 y,pT_j1>50GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+35,'H-j1 y,pT_j1>70GeV ' ,100,-6.d0,6.d0)
+      call HwU_book(l+36,'H-j1 y,pT_j1>90GeV ' ,100,-6.d0,6.d0)
+      
+      call HwU_book(l+37,'njets ' ,11,-0.5d0,10.5d0)
+      call HwU_book(l+38,'njets,|y_j|<2.5 ' ,11,-0.5d0,10.5d0)
+      call HwU_book(l+39,'xsec ' ,3,-0.5d0,2.5d0)
+      return
+      end
+
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_end(dummy)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      double precision dummy
+      call HwU_write_file
+      return                
+      end
+
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'nexternal.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision p(0:4,nexternal)
+      double precision wgts(*)
+      integer ibody
+      double precision wgt,var
+      integer i,kk,l
+      double precision www,pph(0:3),ppj1(0:3),pth,yh,ptj1,yj1,y_central
+     $     ,ptcut,njdble,njcdble
+      integer njet,njet_central,nj
+      double precision getrapidity
+      external getrapidity
+      if (nexternal.ne.4) then
+         write (*,*) 'error #1 in analysis_fill: '/
+     &        /'only for process "p p > h [QCD]"'
+         stop 1
+      endif
+      if (.not. (abs(ipdg(1)).le.5 .or. ipdg(1).eq.21)) then
+         write (*,*) 'error #2 in analysis_fill: '/
+     &        /'only for process "p p > h [QCD]"'
+         stop 1
+      endif
+      if (.not. (abs(ipdg(2)).le.5 .or. ipdg(2).eq.21)) then
+         write (*,*) 'error #3 in analysis_fill: '/
+     &        /'only for process "p p > h [QCD]"'
+         stop 1
+      endif
+      if (.not. (abs(ipdg(4)).le.5 .or. ipdg(4).eq.21)) then
+         write (*,*) 'error #4 in analysis_fill: '/
+     &        /'only for process "p p > h [QCD]"'
+         stop 1
+      endif
+      if (ipdg(3).ne.25) then
+         write (*,*) 'error #5 in analysis_fill: '/
+     &        /'only for process "p p > h [QCD]"'
+         stop 1
+      endif
+c
+      do i=0,3
+         pph(i)=p(i,3)
+         ppj1(i)=p(i,4)
+      enddo
+c Higgs variables
+      pth=sqrt(pph(1)**2+pph(2)**2)
+      yh=getrapidity(pph(0),pph(3))
+c hardest jet variables
+      ptj1=sqrt(ppj1(1)**2+ppj1(2)**2)
+      yj1=getrapidity(ppj1(0),ppj1(3))
+      njet=0
+      njet_central=0
+ 8    y_central=2.5d0
+      ptcut=1d1
+      if(ptj1.gt.ptcut)njet=1
+      if(ptj1.gt.ptcut.and.abs(yj1).le.y_central)njet_central=1
+      njdble=dble(njet)
+      njcdble=dble(njet_central)
+C
+      l=0
+      call HwU_fill(l+1,pth,WGTS)
+      call HwU_fill(l+2,pth,WGTS)
+      if(pth.gt.0.d0)call HwU_fill(l+3,log10(pth),WGTS)
+      if(abs(yh).le.2.d0)then
+         call HwU_fill(l+4,pth,WGTS)
+         call HwU_fill(l+5,pth,WGTS)
+         if(pth.gt.0.d0)call HwU_fill(l+6,log10(pth),WGTS)
+      endif
+c
+      if(njet.ge.1)then
+         call HwU_fill(l+7,ptj1,WGTS)
+         call HwU_fill(l+8,ptj1,WGTS)
+         if(ptj1.gt.0.d0)call HwU_fill(l+9,log10(ptj1),WGTS)
+         if(abs(yj1).le.2.d0)then
+            call HwU_fill(l+10,ptj1,WGTS)
+            call HwU_fill(l+11,ptj1,WGTS)
+            if(ptj1.gt.0.d0)call HwU_fill(l+12,log10(ptj1),WGTS)
+         endif
+c
+         do nj=1,njet
+            call HwU_fill(l+13,ptj1,WGTS)
+            call HwU_fill(l+14,ptj1,WGTS)
+            if(ptj1.gt.0.d0)call HwU_fill(l+15,log10(ptj1),WGTS)
+            if(abs(yj1).le.2.d0)then
+               call HwU_fill(l+16,ptj1,WGTS)
+               call HwU_fill(l+17,ptj1,WGTS)
+               if(ptj1.gt.0d0)call HwU_fill(l+18,log10(ptj1),WGTS)
+            endif
+         enddo
+      endif
+c
+      call HwU_fill(l+19,yh,WGTS)
+      if(pth.ge.10.d0) call HwU_fill(l+20,yh,WGTS)
+      if(pth.ge.30.d0) call HwU_fill(l+21,yh,WGTS)
+      if(pth.ge.50.d0) call HwU_fill(l+22,yh,WGTS)
+      if(pth.ge.70.d0) call HwU_fill(l+23,yh,WGTS)
+      if(pth.ge.90.d0) call HwU_fill(l+24,yh,WGTS)  
+c     
+      if(njet.ge.1)then
+         call HwU_fill(l+25,yj1,WGTS)
+         if(ptj1.ge.10.d0) call HwU_fill(l+26,yj1,WGTS)
+         if(ptj1.ge.30.d0) call HwU_fill(l+27,yj1,WGTS)
+         if(ptj1.ge.50.d0) call HwU_fill(l+28,yj1,WGTS)
+         if(ptj1.ge.70.d0) call HwU_fill(l+29,yj1,WGTS)
+         if(ptj1.ge.90.d0) call HwU_fill(l+30,yj1,WGTS)
+c     
+         call HwU_fill(l+31,yh-yj1,WGTS)
+         if(ptj1.ge.10.d0) call HwU_fill(l+32,yh-yj1,WGTS)
+         if(ptj1.ge.30.d0) call HwU_fill(l+33,yh-yj1,WGTS)
+         if(ptj1.ge.50.d0) call HwU_fill(l+34,yh-yj1,WGTS)
+         if(ptj1.ge.70.d0) call HwU_fill(l+35,yh-yj1,WGTS)
+         if(ptj1.ge.90.d0) call HwU_fill(l+36,yh-yj1,WGTS)
+      endif
+c
+      call HwU_fill(l+37,njdble,WGTS)
+      call HwU_fill(l+38,njcdble,WGTS)
+      call HwU_fill(l+39,1d0,WGTS)
+c      
+C
+ 999  return      
+      end
+
+
+
+      function getrapidity(en,pl)
+      implicit none
+      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
+      parameter (tiny=1.d-8)
+      xplus=en+pl
+      xminus=en-pl
+      if(xplus.gt.tiny.and.xminus.gt.tiny)then
+         if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny)then
+            y=0.5d0*log( xplus/xminus  )
+         else
+            y=sign(1.d0,pl)*1.d8
+         endif
+      else 
+         y=sign(1.d0,pl)*1.d8
+      endif
+      getrapidity=y
+      return
+      end
+
+
+      function getpseudorap(en,ptx,pty,pl)
+      implicit none
+      real*8 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
+      getpseudorap=eta
+      return
+      end

=== added file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_hjj.f'
--- Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_hjj.f	1970-01-01 00:00:00 +0000
+++ Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_hjj.f	2015-04-10 09:55:25 +0000
@@ -0,0 +1,777 @@
+c
+c Example analysis for "p p > h j j [QCD]" process.
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_begin(nwgt,weights_info)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      integer nwgt
+      character*(*) weights_info(*)
+      integer i,kk,l
+      character*8 cc(2)
+      data cc/'        ','vbfcuts '/
+      double precision pi
+      PARAMETER (PI=3.14159265358979312D0)
+      include 'dbook.inc'
+      real*8 vetomin, vetomax
+      integer nbinveto
+      common /to_veto_hist/vetomin,vetomax,nbinveto
+c      
+      call HwU_inithist(nwgt,weights_info)
+      vetomin = 0d0
+      vetomax = 100d0
+      nbinveto = 50
+      do i=1,2
+         l=(i-1)*54
+         call HwU_book(l+  1,'total ' //cc(i),5,0.5d0,5.5d0)
+
+         call HwU_book(l+  2,'Higgs pT ' //cc(i),50,0.d0,400.d0)
+         call HwU_book(l+  3,'Higgs pT ' //cc(i),50,0.d0,800.d0)
+         call HwU_book(l+  4,'Higgs logpT ' //cc(i),50,0.d0,4.d0)
+         call HwU_book(l+  5,'Higgs eta ' //cc(i),50,-6.d0,6.d0)
+         call HwU_book(l+  6,'Higgs y ' //cc(i),50,-6.d0,6.d0)
+
+         call HwU_book(l+  7,'j1 pT ' //cc(i),50,0.d0,400.d0)
+         call HwU_book(l+  8,'j1 pT ' //cc(i),50,0.d0,800.d0)
+         call HwU_book(l+  9,'j1 logpT ' //cc(i),50,0.d0,4.d0)
+         call HwU_book(l+ 10,'j1 eta ' //cc(i),50,-6.d0,6.d0)
+         call HwU_book(l+ 11,'j1 y ' //cc(i),50,-6.d0,6.d0)
+
+         call HwU_book(l+ 12,'j2 pT ' //cc(i),50,0.d0,400.d0)
+         call HwU_book(l+ 13,'j2 pT ' //cc(i),50,0.d0,800.d0)
+         call HwU_book(l+ 14,'j2 logpT ' //cc(i),50,0.d0,4.d0)
+         call HwU_book(l+ 15,'j2 eta ' //cc(i),50,-6.d0,6.d0)
+         call HwU_book(l+ 16,'j2 y ' //cc(i),50,-6.d0,6.d0)
+
+         call HwU_book(l+ 17,'j3 pT ' //cc(i),50,0.d0,400.d0)
+         call HwU_book(l+ 18,'j3 pT ' //cc(i),50,0.d0,800.d0)
+         call HwU_book(l+ 19,'j3 logpT ' //cc(i),50,0.d0,4.d0)
+         call HwU_book(l+ 20,'j3 eta ' //cc(i),50,-6.d0,6.d0)
+         call HwU_book(l+ 21,'j3 y ' //cc(i),50,-6.d0,6.d0)
+
+         call HwU_book(l+ 22,'H+j1 pT ' //cc(i),50,0.d0,400.d0)
+         call HwU_book(l+ 23,'H+j1 pT ' //cc(i),50,0.d0,800.d0)
+         call HwU_book(l+ 24,'H+j1 logpT ' //cc(i),50,0.d0,4.d0)
+         call HwU_book(l+ 25,'H+j1 eta ' //cc(i),50,-6.d0,6.d0)
+         call HwU_book(l+ 26,'H+j1 y ' //cc(i),50,-6.d0,6.d0)
+
+         call HwU_book(l+ 27,'j1+j2 pT ' //cc(i),50,0.d0,400.d0)
+         call HwU_book(l+ 28,'j1+j2 pT ' //cc(i),50,0.d0,800.d0)
+         call HwU_book(l+ 29,'j1+j2 logpT ' //cc(i),50,0.d0,4.d0)
+         call HwU_book(l+ 30,'j1+j2 eta ' //cc(i),50,-6.d0,6.d0)
+         call HwU_book(l+ 31,'j1+j2 y ' //cc(i),50,-6.d0,6.d0)
+
+         call HwU_book(l+ 32,'syst pT ' //cc(i),50,0.d0,400.d0)
+         call HwU_book(l+ 33,'syst pT ' //cc(i),50,0.d0,800.d0)
+         call HwU_book(l+ 34,'syst logpT ' //cc(i),50,0.d0,4.d0)
+         call HwU_book(l+ 35,'syst eta ' //cc(i),50,-10.d0,10.d0)
+         call HwU_book(l+ 36,'syst y ' //cc(i),50,-6.d0,6.d0)
+
+         call HwU_book(l+ 37,'Dphi H-j1 ' //cc(i),50,0d0,pi)
+         call HwU_book(l+ 38,'Dphi H-j2 ' //cc(i),50,0d0,pi)
+         call HwU_book(l+ 39,'Dphi j1-j2 ' //cc(i),50,0d0,pi)
+         
+         call HwU_book(l+ 40,'DR H-j1 ' //cc(i),50,0d0,10.d0)
+         call HwU_book(l+ 41,'DR H-j2 ' //cc(i),50,0d0,10.d0)
+         call HwU_book(l+ 42,'DR j1-j2 ' //cc(i),50,0d0,10.d0)
+
+         call HwU_book(l+ 43,'mj1j2 ' //cc(i),50,0d0,3000.d0)
+
+c Nason-Oleari plots (hep-ph/0911.5299)
+         call HwU_book(l+ 44,'|yj1-yj2| ' //cc(i),25,0.d0,10.d0)
+         call HwU_book(l+ 45,'yj3_rel ' //cc(i),50,-6.d0,6.d0)
+         call HwU_book(l+ 46,'njets ' //cc(i),100,-0.5d0,9.5d0)
+         call HwU_book(l+ 47,'ptrel_j1 ' //cc(i),50,0.d0,200.d0)
+         call HwU_book(l+ 48,'ptrel_j2 ' //cc(i),50,0.d0,200.d0)
+         call HwU_book(l+ 49,'P-veto ' //cc(i), dble(nbinveto),vetomin
+     &        ,vetomax)
+         call HwU_book(l+ 50,'jveto pT ' //cc(i), dble(nbinveto),vetomin
+     &        ,vetomax)
+         call HwU_book(l+ 51,'jveto pT ' //cc(i), dble(nbinveto),
+     &        vetomin,2d0*vetomax)
+         call HwU_book(l+ 52,'jveto logpT ' //cc(i),50,0.d0,4.d0)
+         call HwU_book(l+ 53,'jveto eta ' //cc(i),50,-6.d0,6.d0)
+         call HwU_book(l+ 54,'jveto y ' //cc(i),50,-6.d0,6.d0)
+
+      enddo
+      return
+      end
+
+
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_end(dummy)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      double precision dummy
+      call HwU_write_file
+      return                
+      end
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'nexternal.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision p(0:4,nexternal)
+      double precision wgts(*)
+      integer ibody
+      double precision wgt,var
+      integer i,j,k,kk,l
+      double precision www,pQCD(0:3,nexternal),palg,rfj,sycut,yjmax
+     $     ,pjet(0:3,nexternal),ptjet(nexternal),yjet(nexternal)
+     $     ,etajet(nexternal),ptj_tag,deltay12,mj1j2min,ph(0:3),pj1(0:3)
+     $     ,pj2(0:3),pj3(0:3),pjj(0:3),pHj(0:3),psyst(0:3),pjveto(0:3)
+     $     ,ptH,etaH,yH,njdble,ptj1,etaj1,yj1,ptHj,etaHj,yHj,DphiHj1
+     $     ,DRHj1,ptj2,etaj2,yj2,ptjj,etajj,yjj,ptsyst,etasyst,ysyst
+     $     ,DphiHj2,Dphij1j2,DRHj2,DRj1j2,mj1j2,Dyj1j2,ptj3,etaj3,yj3
+     $     ,yj3rel,chy1,shy1,chy1mo,chy2,shy2,chy2mo,ptrel_j1,ptrel_j2
+     $     ,ppboost(0:3,nexternal),prel_j1(0:3),prel_j2(0:3)
+     $     ,pj1boost(0:3),pj2boost(0:3),pt_veto,xsecup2
+      logical pass_tag_cuts,flag
+      integer nQCD,jet(nexternal),ij1y,ij2y,ij3y,njet,njety,ijveto
+     $     ,ijvetoy,ij1,ij2,ij3
+      double precision getptv4,getrapidityv4,getpseudorapv4,getdelphiv4
+     $     ,getdrv4,getinvmv4,getmod
+      external getptv4,getrapidityv4,getpseudorapv4,getdelphiv4,getdrv4
+     $     ,getinvmv4,getmod
+      real*8 vetomin, vetomax
+      integer nbinveto
+      common /to_veto_hist/vetomin,vetomax,nbinveto
+      real*8 xd(1:3)
+      data (xd(i),i=1,3)/0d0,0d0,1d0/
+      if (nexternal.ne.6) then
+         write (*,*) 'error #1 in analysis_fill: '/
+     &        /'only for process "p p > h j j [QCD]"'
+         stop 1
+      endif
+      if (.not. (abs(ipdg(1)).le.5 .or. ipdg(1).eq.21)) then
+         write (*,*) 'error #2 in analysis_fill: '/
+     &        /'only for process "p p > h j j [QCD]"'
+         stop 1
+      endif
+      if (.not. (abs(ipdg(2)).le.5 .or. ipdg(2).eq.21)) then
+         write (*,*) 'error #3 in analysis_fill: '/
+     &        /'only for process "p p > h j j [QCD]"'
+         stop 1
+      endif
+      if (.not. (abs(ipdg(4)).le.5 .or. ipdg(4).eq.21)) then
+         write (*,*) 'error #4 in analysis_fill: '/
+     &        /'only for process "p p > h j j [QCD]"'
+         stop 1
+      endif
+      if (.not. (abs(ipdg(5)).le.5 .or. ipdg(5).eq.21)) then
+         write (*,*) 'error #5 in analysis_fill: '/
+     &        /'only for process "p p > h j j [QCD]"'
+         stop 1
+      endif
+      if (.not. (abs(ipdg(6)).le.5 .or. ipdg(6).eq.21)) then
+         write (*,*) 'error #6 in analysis_fill: '/
+     &        /'only for process "p p > h j j [QCD]"'
+         stop 1
+      endif
+      if (ipdg(3).ne.25) then
+         write (*,*) 'error #7 in analysis_fill: '/
+     &        /'only for process "p p > h j j [QCD]"'
+         stop 1
+      endif
+c
+
+c Put all (light) QCD partons in momentum array for jet clustering.
+      nQCD=0
+      do j=nincoming+1,nexternal
+         if (abs(ipdg(j)).le.5 .or. ipdg(j).eq.21) then
+            nQCD=nQCD+1
+            do i=0,3
+               pQCD(i,nQCD)=p(i,j) 
+            enddo
+         endif
+      enddo
+      
+C---CLUSTER THE EVENT
+      palg  = -1.d0
+      rfj   = 0.4d0
+      sycut = 20d0
+      yjmax = 4.5d0
+      do i=1,nexternal
+         do j=0,3
+            pjet(j,i)=0d0
+         enddo
+         ptjet(i)=0d0
+         yjet(i)=0d0
+         etajet(i)=0d0
+         jet(i)=0
+      enddo
+      ij1y=0
+      ij2y=0
+      ij3y=0
+      njet=0
+      njety=0
+      ijveto = 0
+      ijvetoy = 0
+
+
+c******************************************************************************
+c     call FASTJET to get all the jets
+c
+c     INPUT:
+c     input momenta:               pQCD(0:3,nexternal), energy is 0th component
+c     number of input momenta:     nQCD
+c     radius parameter:            rfj
+c     minumum jet pt:              sycut
+c     jet algorithm:               palg, 1.0=kt, 0.0=C/A, -1.0 = anti-kt
+c
+c     OUTPUT:
+c     jet momenta:                             pjet(0:3,nexternal), E is 0th cmpnt
+c     the number of jets (with pt > SYCUT):    njet
+c     the jet for a given particle 'i':        jet(i),   note that this is
+c     the particle in pQCD, which doesn't necessarily correspond to the particle
+c     label in the process
+c
+      call amcatnlo_fastjetppgenkt(pQCD,nQCD,rfj,sycut,palg,pjet,njet
+     $     ,jet)
+c
+c******************************************************************************
+      do i=1,njet
+         ptjet(i)=getptv4(pjet(0,i))
+         if(i.gt.1)then
+            if (ptjet(i).gt.ptjet(i-1)) then
+               write (*,*) "Error 1: jets should be ordered in pt"
+               stop
+            endif
+         endif
+         yjet(i)=getrapidityv4(pjet(0,i))
+         etajet(i)=getpseudorapv4(pjet(0,i))
+c look for veto jet without y cuts
+         if (i.gt.2.and.yjet(i).gt.min(yjet(1),yjet(2)).and.
+     &        yjet(i).lt.max(yjet(1),yjet(2)).and.ijveto.eq.0) ijveto=i
+
+C now look for jets within the rapidity cuts
+         if (dabs(yjet(i)).lt.yjmax) then
+            njety=njety+1
+            if (ij1y.eq.0) then
+               ij1y = i
+            else if (ij2y.eq.0) then
+               ij2y = i
+            else if (ij3y.eq.0) then
+               ij3y = i
+            endif
+c look for veto jet with y cuts
+            if (ij3y.gt.0.and.
+     &           yjet(i).gt.min(yjet(ij1y),yjet(ij2y)).and.
+     &           yjet(i).lt.max(yjet(ij1y),yjet(ij2y)).and.ijvetoy.eq.0) 
+     &           ijvetoy = i
+         endif
+      enddo
+
+c Nason-Oleari cuts (hep-ph/0911.5299)
+      ptj_tag  = 20d0
+      deltay12 = 4.d0
+      mj1j2min = 600d0
+
+c this is the loop for w-o / w vbf cuts
+      do i=1,2
+      if(i.eq.1) then 
+         ij1 = 1
+         ij2 = 2
+         ij3 = 3
+      endif
+      if(i.eq.2) then
+         njet = njety
+         ijveto = ijvetoy
+         ij1 = ij1y
+         ij2 = ij2y
+         ij3 = ij3y
+      endif
+
+c Load momenta
+      xsecup2=1d0
+      do k=0,3
+         pH(k)   =p(k,3)
+         pj1(k)  =pjet(k,ij1)
+         pj2(k)  =pjet(k,ij2)
+         pj3(k)  =pjet(k,ij3)
+         pjj(k)  =pjet(k,ij1)+pjet(k,ij2)
+         pHj(k)  =pjet(k,ij1)+pH(k)
+         psyst(k)=pjet(k,ij1)+pjet(k,ij2)+pH(k)
+         pjveto(k)=pjet(k,ijveto)
+      enddo
+
+c Define observables
+c Higgs
+         ptH     = getptv4(pH)
+         etaH    = getpseudorapv4(pH)
+         yH      = getrapidityv4(pH)
+         njdble  = dble(njet)
+c At least one jet
+      if(njet.ge.1)then
+        ptj1    = getptv4(pj1)
+        etaj1   = getpseudorapv4(pj1)
+        yj1     = getrapidityv4(pj1)
+        ptHj    = getptv4(pHj)
+        etaHj   = getpseudorapv4(pHj)
+        yHj     = getrapidityv4(pHj)
+        DphiHj1 = getdelphiv4(pH,pj1)
+        DRHj1   = getdrv4(pH,pj1)
+      endif
+c At least two jets
+      if(njet.ge.2)then
+        ptj2    = getptv4(pj2)
+        etaj2   = getpseudorapv4(pj2)
+        yj2     = getrapidityv4(pj2)
+        ptjj    = getptv4(pjj)
+        etajj   = getpseudorapv4(pjj)
+        yjj     = getrapidityv4(pjj)
+        ptsyst  = getptv4(psyst)
+        etasyst = getpseudorapv4(psyst)
+        ysyst   = getrapidityv4(psyst)
+        DphiHj2 = getdelphiv4(pH,pj2)
+        Dphij1j2= getdelphiv4(pj1,pj2)
+        DRHj2   = getdrv4(pH,pj2)
+        DRj1j2  = getdrv4(pj1,pj2)
+        mj1j2   = getinvmv4(pjj)
+        Dyj1j2  = abs(yj1-yj2)
+      endif
+c At least three jets
+      if(njet.ge.3)then
+        ptj3    = getptv4(pj3)
+        etaj3   = getpseudorapv4(pj3)
+        yj3     = getrapidityv4(pj3)
+        yj3rel  = yj3-(yj1+yj2)/2d0
+      endif
+c
+      chy1=cosh(yj1)
+      shy1=sinh(yj1)
+      chy1mo=chy1-1.d0
+      chy2=cosh(yj2)
+      shy2=sinh(yj2)
+      chy2mo=chy2-1.d0
+
+      call boostwdir2(chy1,shy1,chy1mo,xd,pj1,pj1boost)
+      call boostwdir2(chy2,shy2,chy2mo,xd,pj2,pj2boost)
+      ptrel_j1=0d0
+      ptrel_j2=0d0
+
+      pass_tag_cuts = njety.ge.2 .and.
+     &                ptj1.ge.ptj_tag .and.
+     &                ptj2.ge.ptj_tag .and.
+     &                abs(yj1-yj2).ge.deltay12 .and.
+     &                yj1*yj2.le.0d0 .and.
+     &                mj1j2.ge.mj1j2min 
+
+      if(i.eq.1) then 
+         flag=.true.
+      endif
+
+      if(i.eq.2) then
+         flag=pass_tag_cuts
+      endif
+
+
+      do j=1,nQCD
+         if(njet.ge.1.and.jet(j).eq.1)then
+           call boostwdir2(chy1,shy1,chy1mo,xd,pQCD(0,j),ppboost(0,j))
+           call getwedge(ppboost(0,j),pj1boost,prel_j1)
+           ptrel_j1=ptrel_j1+getmod(prel_j1)/getmod(pj1boost)
+         elseif(njet.ge.2.and.jet(j).eq.2)then
+           call boostwdir2(chy2,shy2,chy2mo,xd,pQCD(0,j),ppboost(0,j))
+           call getwedge(ppboost(0,j),pj2boost,prel_j2)
+           ptrel_j2=ptrel_j2+getmod(prel_j2)/getmod(pj2boost)
+         endif
+      enddo
+
+      l=(i-1)*54
+      if(flag)then
+         call HwU_fill(l+  1,1d0,wgts)
+         call HwU_fill(l+  2,ptH,wgts)
+         call HwU_fill(l+  3,ptH,wgts)
+         if(ptH.gt.0d0) call HwU_fill(l+  4,log10(ptH),wgts)
+         call HwU_fill(l+  5,etaH,wgts)
+         call HwU_fill(l+  6,yH,wgts)
+         call HwU_fill(l+ 46,njdble,wgts)
+
+         if(njet.ge.1)then
+            call HwU_fill(l+  7,ptj1,wgts)
+            call HwU_fill(l+  8,ptj1,wgts)
+            if (ptj1.gt.0d0) call HwU_fill(l+  9,log10(ptj1),wgts)
+            call HwU_fill(l+ 10,etaj1,wgts)
+            call HwU_fill(l+ 11,yj1,wgts)
+            call HwU_fill(l+ 22,ptHj,wgts)
+            call HwU_fill(l+ 23,ptHj,wgts)
+            if(ptHj.gt.0d0) call HwU_fill(l+ 24,log10(ptHj),wgts)
+            call HwU_fill(l+ 25,etaHj,wgts)
+            call HwU_fill(l+ 26,yHj,wgts)
+            call HwU_fill(l+ 37,DphiHj1,wgts)
+            call HwU_fill(l+ 40,DRHj1,wgts)
+            call HwU_fill(l+ 47,ptrel_j1,wgts)
+         endif 
+         
+         if(njet.ge.2)then
+            call HwU_fill(l+ 12,ptj2,wgts)
+            call HwU_fill(l+ 13,ptj2,wgts)
+            if (ptj2.gt.0d0) call HwU_fill(l+ 14,log10(ptj2),wgts)
+            call HwU_fill(l+ 15,etaj2,wgts)
+            call HwU_fill(l+ 16,yj2,wgts)
+            call HwU_fill(l+ 27,ptjj,wgts)
+            call HwU_fill(l+ 28,ptjj,wgts)
+            if(ptjj.gt.0d0) call HwU_fill(l+ 29,log10(ptjj),wgts)
+            call HwU_fill(l+ 30,etajj,wgts)
+            call HwU_fill(l+ 31,yjj,wgts)
+            call HwU_fill(l+ 32,ptsyst,wgts)
+            call HwU_fill(l+ 33,ptsyst,wgts)
+            if(ptsyst.gt.0d0) call HwU_fill(l+ 34,log10(ptsyst),wgts)
+            call HwU_fill(l+ 35,etasyst,wgts)
+            call HwU_fill(l+ 36,ysyst,wgts)
+            call HwU_fill(l+ 38,DphiHj2,wgts)
+            call HwU_fill(l+ 39,Dphij1j2,wgts)
+            call HwU_fill(l+ 41,DRHj2,wgts)
+            call HwU_fill(l+ 42,DRj1j2,wgts)
+            call HwU_fill(l+ 43,mj1j2,wgts)
+            call HwU_fill(l+ 44,Dyj1j2,wgts)
+            call HwU_fill(l+ 48,ptrel_j2,wgts)
+         endif
+         
+         if(njet.ge.3)then
+            call HwU_fill(l+ 17,ptj3,wgts)
+            call HwU_fill(l+ 18,ptj3,wgts)
+            if(ptj3.gt.0d0) call HwU_fill(l+ 19,log10(ptj3),wgts)
+            call HwU_fill(l+ 20,etaj3,wgts)
+            call HwU_fill(l+ 21,yj3,wgts)
+            call HwU_fill(l+ 45,yj3rel,wgts)
+         endif
+         if (ijveto.gt.0) then
+            pt_veto = getptv4(pjveto)
+            do k=1,nbinveto
+               if (pt_veto.gt. (vetomin+(vetomax-vetomin)*dble(k-1)
+     $              /dble(nbinveto))) then
+                  call HwU_fill(l+49, (vetomax-vetomin)* dble(k)
+     $                 /dble(nbinveto)*0.99, wgts/xsecup2)
+               endif
+            enddo
+            call HwU_fill(l+50,pt_veto,wgts)