Merge lp:~olivier-mattelaer/mg5amcnlo/2.5.3 into lp:mg5amcnlo/lts

[Olivier Mattelaer]

Hi,

Any objection for releasing 2.5.3?

Cheers,

Olivier

[lp:~olivier-mattelaer/mg5amcnlo/2.5.3]

302. By olivier-mattelaer on 2017-03-04

fixing one small bug

303. By olivier-mattelaer on 2017-03-06

fixing a tricky error related to multi_run with pgs

304. By olivier-mattelaer on 2017-03-06

ensure that the group number is correctly propagated even on distributted system

305. By Valentin Hirschi on 2017-03-07

1. Fixed an issues found by Huasheng on the computation of the loop symmetry factor in very special cases.

306. By Hua-Sheng Shao on 2017-03-07

1) add suffice for avh_olo.f90 and avh_olo.o in create.py of IREGI/src/oneloop; 2) Correct one lhaid of PDF4LHC15 in histograms.py

307. By Valentin Hirschi on 2017-03-07

1. Enforced sequential build of CutTools in our compilation steered by MG5aMC.
   Nowadays almost all version of make support the '-j' options and an increasing
   number of environment are now automatically turning on parallel build, so this must
   be done.
   Ideally one would either fix CutTools target dependencies to support paralell build
   or enforce sequential build directly in its makefile, but unfortunately none of these
   two options are easy to implement.

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

I guess this means that everyone approves. So I freeze the version here.
If you want to have something in 2.5.3. Please contact me since apriori any version after this one will not be part of 2.5.3 release. (I might still make one or two push depending if some of the tests fails or not)

Cheers,

Olivier

review: Approve

[lp:~olivier-mattelaer/mg5amcnlo/2.5.3]

308. By Rikkert Frederix on 2017-03-08

small improvement in the length of the buffer. Needed for processes
with more than 100 integration channels at the Born level.

[Rikkert Frederix (frederix)]

Hi Olivier,

Just submitted and pushed a final improvement. So, now it's good to go.

Cheers,
Rik

review: Approve

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2016-12-09 00:39:21 +0000
+++ MadSpin/decay.py	2017-03-08 10:27:27 +0000
@@ -1919,7 +1919,6 @@
         sd=sd/(N-1.0)
         
         return mean, sd
-     
 
 class decay_all_events(object):
     
@@ -2026,7 +2025,6 @@
         self.all_decay = {}
 
 
- 
         # generate BR and all the square matrix element based on the banner.
         pickle_info = pjoin(self.path_me,"production_me", "all_ME.pkl")
         if not options["use_old_dir"] or not os.path.exists(pickle_info):
@@ -2042,7 +2040,7 @@
                 save_load_object.save_to_file(pickle_info,
                                           (self.all_ME,self.all_decay,self.width_estimator))                
         
-        if not self.options["onlyhelicity"]:
+        if not self.options["onlyhelicity"] and self.options['spinmode'] != 'onshell':
             resonances = self.width_estimator.resonances
             logger.debug('List of resonances: %s' % resonances)
             self.extract_resonances_mass_width(resonances) 
@@ -2669,9 +2667,8 @@
                 #self.banner.get('proc_card').get('multiparticles'):
                 mgcmd.do_define("%s = %s" % (name, ' '.join(`i` for i in pdgs)))
             
-
+        
         mgcmd.exec_cmd("set group_subprocesses False")
-
         logger.info('generating the production square matrix element')
         start = time.time()
         commandline=''
@@ -2762,7 +2759,7 @@
 
         # 4. compute the full matrix element -----------------------------------
         if not self.options["onlyhelicity"]:
-            logger.info('generating the full square matrix element (with decay)')
+            logger.info('generating the full matrix element squared (with decay)')
             start = time.time()
             to_decay = self.mscmd.list_branches.keys()
             decay_text = []
@@ -4009,3 +4006,283 @@
                 del external
 
         self.calculator = {}
+
+
+
+
+class decay_all_events_onshell(decay_all_events):
+    """special mode for onshell production"""
+
+    @misc.mute_logger()
+    @misc.set_global()
+    def generate_all_matrix_element(self):
+        """generate the full series of matrix element needed by Madspin.
+        i.e. the undecayed and the decay one. And associate those to the 
+        madspin production_topo object"""
+
+        # 1. compute the partial width        
+        # 2. compute the production matrix element
+        # 3. create the all_topology object 
+        # 4. compute the full matrix element (use the partial to throw away 
+        #     pointless decay.
+        # 5. add the decay information to the all_topology object (with branching
+        #     ratio)  
+        
+        
+        # 0. clean previous run ------------------------------------------------
+        path_me = self.path_me
+        try:
+            shutil.rmtree(pjoin(path_me,'madspin_me'))
+        except Exception: 
+            pass       
+        
+        # 1. compute the partial width------------------------------------------
+        #self.get_branching_ratio()
+        
+        # 2. compute the production matrix element -----------------------------
+        processes = [line[9:].strip() for line in self.banner.proc_card
+                     if line.startswith('generate')]
+        processes += [' '.join(line.split()[2:]) for line in self.banner.proc_card
+                      if re.search('^\s*add\s+process', line)]
+        
+        mgcmd = self.mgcmd
+        modelpath = self.model.get('modelpath+restriction')
+
+        commandline="import model %s" % modelpath
+        if not self.model.mg5_name:
+            commandline += ' --modelname'
+            
+        mgcmd.exec_cmd(commandline)
+        # Handle the multiparticle of the banner        
+        #for name, definition in self.mscmd.multiparticles:
+        if hasattr(self.mscmd, 'multiparticles_ms'):
+            for name, pdgs in  self.mscmd.multiparticles_ms.items():
+                if name == 'all':
+                    continue
+                #self.banner.get('proc_card').get('multiparticles'):
+                mgcmd.do_define("%s = %s" % (name, ' '.join(`i` for i in pdgs)))
+            
+        
+        mgcmd.exec_cmd("set group_subprocesses False")
+        logger.info('generating the production square matrix element')
+        start = time.time()
+        commandline=''
+        for proc in processes:
+            # deal with @ syntax need to move it after the decay specification
+            if '@' in proc:
+                proc, proc_nb = proc.split('@')
+                try:
+                    proc_nb = int(proc_nb)
+                except ValueError:
+                    raise MadSpinError, 'MadSpin didn\'t allow order restriction after the @ comment: \"%s\" not valid' % proc_nb
+                proc_nb = '@ %i' % proc_nb 
+            else:
+                proc_nb = ''
+                        
+            if ',' in proc:
+                raise MadSpinError, 'MadSpin can not decay event which comes from a decay chain.'+\
+                        '\n  The full decay chain should either be handle by MadGraph or by Masdspin.'
+            
+            if '[' not in proc:
+                commandline+="add process %s %s  --no_warning=duplicate;" % (proc, proc_nb)
+            else:
+                process, order, final = re.split('\[\s*(.*)\s*\]', proc)
+                commandline+="add process %s %s --no_warning=duplicate;" % (process, proc_nb)
+                if not order:
+                    continue
+                elif not order.startswith('virt='):
+                    if '=' in order:
+                        order = order.split('=',1)[1]
+                    # define the list of particles that are needed for the radiateion
+                    pert = fks_common.find_pert_particles_interactions(
+                         mgcmd._curr_model,pert_order = order)['soft_particles']
+                    commandline += "define pert_%s = %s;" % (order, ' '.join(map(str,pert)) )
+                    
+                    # check if we have to increase by one the born order
+                    if '%s=' % order in process:
+                        result=re.split(' ',process)
+                        process=''
+                        for r in result:
+                            if '%s=' % order in r:
+                                ior=re.split('=',r)
+                                r='QCD=%i' % (int(ior[1])+1)
+                            process=process+r+' '
+                    #handle special tag $ | / @
+                    result = re.split('([/$@]|\w+=\w+)', process, 1)                    
+                    if len(result) ==3:
+                        process, split, rest = result
+                        commandline+="add process %s pert_%s %s%s %s --no_warning=duplicate;" % (process, order ,split, rest, proc_nb)
+                    else:
+                        commandline +='add process %s pert_%s  %s --no_warning=duplicate;' % (process,order, proc_nb)                                       
+#        commandline = commandline.replace('add process', 'generate',1)
+#        logger.info(commandline)
+#        
+#        mgcmd.exec_cmd(commandline, precmd=True)
+#        commandline = 'output standalone_msP %s %s' % \
+#        (pjoin(path_me,'production_me'), ' '.join(self.list_branches.keys()))        
+#        mgcmd.exec_cmd(commandline, precmd=True)        
+#        logger.info('Done %.4g' % (time.time()-start))
+
+        # 3. Create all_ME + topology objects ----------------------------------
+#        matrix_elements = mgcmd._curr_matrix_elements.get_matrix_elements()
+#        self.all_ME.adding_me(matrix_elements, pjoin(path_me,'production_me'))
+        
+        # 4. compute the full matrix element -----------------------------------
+        logger.info('generating the full matrix element squared (with decay)')
+#        start = time.time()
+        to_decay = self.mscmd.list_branches.keys()
+        decay_text = []
+        for decays in self.mscmd.list_branches.values():
+            for decay in  decays:
+                if '=' not in decay:
+                    decay += ' QCD=99'
+                if ',' in decay:
+                    decay_text.append('(%s)' % decay)
+                else:
+                    decay_text.append(decay)
+        decay_text = ', '.join(decay_text)
+#        commandline = ''
+        
+        
+        for proc in processes:
+            # deal with @ syntax need to move it after the decay specification
+            if '@' in proc:
+                proc, proc_nb = proc.split('@')
+                try:
+                    proc_nb = int(proc_nb)
+                except ValueError:
+                    raise MadSpinError, 'MadSpin didn\'t allow order restriction after the @ comment: \"%s\" not valid' % proc_nb
+                proc_nb = '@ %i' % proc_nb 
+            else:
+                proc_nb = '' 
+            
+            if '[' not in proc:
+                nb_comma = proc.count(',')
+                if nb_comma == 0:
+                    commandline+="add process %s, %s %s  --no_warning=duplicate;" % (proc, decay_text, proc_nb)
+                elif nb_comma == 1:
+                    before, after = proc.split(',')
+                    commandline+="add process %s, %s, (%s, %s) %s  --no_warning=duplicate;" % (before, decay_text, after, decay_text, proc_nb)
+                else:
+                    raise Exception, 'too much decay at MG level. this can not be done for the moment)'
+            else:
+                process, order, final = re.split('\[\s*(.*)\s*\]', proc)
+                commandline+="add process %s, %s %s  --no_warning=duplicate;" % (process, decay_text, proc_nb)
+                if not order:
+                    continue
+                elif not order.startswith('virt='):
+                    if '=' in order:
+                        order = order.split('=',1)[1]
+                    # define the list of particles that are needed for the radiateion
+                    pert = fks_common.find_pert_particles_interactions(
+                         mgcmd._curr_model,pert_order = order)['soft_particles']
+                    commandline += "define pert_%s = %s;" % (order, ' '.join(map(str,pert)) )
+                    
+                    # check if we have to increase by one the born order
+                    if '%s=' % order in process:
+                        result=re.split(' ',process)
+                        process=''
+                        for r in result:
+                            if '%s=' % order in r:
+                                ior=re.split('=',r)
+                                r='QCD=%i' % (int(ior[1])+1)
+                            process=process+r+' '
+                    #handle special tag $ | / @
+                    result = re.split('([/$@]|\w+=\w+)', process, 1)                    
+                    if len(result) ==3:
+                        process, split, rest = result
+                        commandline+="add process %s pert_%s %s%s , %s %s --no_warning=duplicate;" % \
+                              (process, order, split, rest, decay_text, proc_nb)
+                    else:
+                        commandline +='add process %s pert_%s, %s %s  --no_warning=duplicate;' % \
+                                           (process, order, decay_text, proc_nb)
+#        commandline = commandline.replace('add process', 'generate',1)
+#        logger.info(commandline)
+#        mgcmd.exec_cmd(commandline, precmd=True)
+        # remove decay with 0 branching ratio.
+#        mgcmd.remove_pointless_decay(self.banner.param_card)
+#        commandline = 'output standalone_msF %s %s' % (pjoin(path_me,'full_me'),
+#                                                      ' '.join(self.list_branches.keys()))
+#        mgcmd.exec_cmd(commandline, precmd=True)
+#        logger.info('Done %.4g' % (time.time()-start))
+                 
+
+
+
+        # 5. add the decay information to the all_topology object --------------                        
+#        for matrix_element in mgcmd._curr_matrix_elements.get_matrix_elements():
+#            me_path = pjoin(path_me,'full_me', 'SubProcesses', \
+#                       "P%s" % matrix_element.get('processes')[0].shell_string())
+#            self.all_ME.add_decay(matrix_element, me_path)
+
+        # 5.b import production matrix elements (+ related info) in the full process directory
+#        list_prodfiles=['matrix_prod.f','configs_production.inc','props_production.inc','nexternal_prod.inc']
+#        for tag in self.all_ME:
+#            prod_path=self.all_ME[tag]['path']
+#            nfinal=len(self.all_ME[tag]['base_order'][1])
+#            for dico in self.all_ME[tag]['decays']:
+#                full_path=dico['path']
+                #print prod_path
+                #print full_path
+                #print ' '
+#                for item in list_prodfiles:
+                     #print full_path
+#                     prodfile=pjoin(prod_path,item)
+#                     destination=pjoin(full_path,item)
+#                     shutil.copyfile(prodfile, destination)  
+                # we need to write the file config_decays.inc
+#                self.generate_configs_file(nfinal,dico,full_path)
+                
+
+#        if self.options["onlyhelicity"]:
+#            return
+
+        # 6. generate decay only part ------------------------------------------
+        logger.info('generate matrix element for decay only (1 - > N).')
+#        start = time.time()
+#        commandline = ''
+        i=0
+        for processes in self.list_branches.values():
+            for proc in processes:
+                commandline+="add process %s @%i --no_warning=duplicate --standalone;" % (proc,i)
+                i+=1        
+        commandline = commandline.replace('add process', 'generate',1)
+        mgcmd.exec_cmd(commandline, precmd=True)
+        # remove decay with 0 branching ratio.
+        #mgcmd.remove_pointless_decay(self.banner.param_card)
+        #
+        commandline = 'output standalone %s' % pjoin(path_me,'madspin_me')
+        logger.info(commandline)
+        mgcmd.exec_cmd(commandline, precmd=True)
+        logger.info('Done %.4g' % (time.time()-start))  
+        self.all_me = {}
+        
+        # store information about matrix element
+        for matrix_element in mgcmd._curr_matrix_elements.get_matrix_elements():
+            me_string = matrix_element.get('processes')[0].shell_string()
+            for me in matrix_element.get('processes'):   
+                dirpath = pjoin(path_me,'madspin_me', 'SubProcesses', "P%s" % me_string)
+                # get the orignal order:
+                initial = []
+                final = [l.get('id') for l in me.get_legs_with_decays()\
+                          if l.get('state') or initial.append(l.get('id'))]
+                order = (tuple(initial), tuple(final))
+                initial.sort(), final.sort()
+                tag = (tuple(initial), tuple(final))
+                self.all_me[tag] = {'pdir': "P%s" % me_string, 'order': order}
+
+
+        return self.all_me
+    
+
+
+    def compile(self):
+        logger.info('Compiling code')
+        misc.compile(cwd=pjoin(self.path_me,'madspin_me', 'Source'),
+                     nb_core=self.mgcmd.options['nb_core'])        
+        misc.compile(['all'],cwd=pjoin(self.path_me,'madspin_me', 'SubProcesses'),
+                     nb_core=self.mgcmd.options['nb_core'])
+
+
+    
+    

=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2016-11-25 15:38:07 +0000
+++ MadSpin/interface_madspin.py	2017-03-08 10:27:27 +0000
@@ -21,6 +21,8 @@
 import random
 import re
 import shutil
+import sys
+import time
 import glob
 
 
@@ -33,6 +35,7 @@
 import madgraph.interface.madgraph_interface as mg_interface
 import madgraph.interface.master_interface as master_interface
 import madgraph.interface.madevent_interface as madevent_interface
+import madgraph.interface.reweight_interface as rwgt_interface
 import madgraph.various.misc as misc
 import madgraph.iolibs.files as files
 import madgraph.iolibs.export_v4 as export_v4
@@ -44,6 +47,7 @@
 import MadSpin.decay as madspin
 
 
+
 logger = logging.getLogger('decay.stdout') # -> stdout
 logger_stderr = logging.getLogger('decay.stderr') # ->stderr
 cmd_logger = logging.getLogger('cmdprint2') # -> print
@@ -128,7 +132,7 @@
                 inputfile = inputfile + '.gz'
             else: 
                 raise self.InvalidCmd('No such file or directory : %s' % inputfile)
-        
+
         if inputfile.endswith('.gz'):
             misc.gunzip(inputfile)
             inputfile = inputfile[:-3]
@@ -371,8 +375,8 @@
                 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")
+            if args[1].lower() not in ["full", "onshell", "none"]:
+                raise self.InvalidCmd("spinmode can only take one of those 3 value: full/onshell/none")
              
         elif args[0] == "run_card":
             if self.options['spinmode'] == "madspin":
@@ -438,7 +442,7 @@
                                                    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)
+            return self.list_completion(text, ["full","onshell", "none"], line)
          
     def help_set(self):
         """help the set command"""
@@ -516,14 +520,16 @@
         launch
         '''
 
-    #@misc.mute_logger()
+    @misc.mute_logger()
     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"] == "onshell":
+            return self.run_onshell(line)
         elif self.options["spinmode"] == "bridge":
-            raise Exception, "Bridge mode not yet available due to lack of validation."
+            raise Exception, "Bridge mode not available."
         
         if self.options['ms_dir'] and os.path.exists(pjoin(self.options['ms_dir'], 'madspin.pkl')):
             return self.run_from_pickle()
@@ -693,6 +699,8 @@
         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"""
         
@@ -700,111 +708,7 @@
         #   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)
-            if not part:
-                return {}# this particle is not defined in the current model so ignore it
-            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)
-                    if self.options['ms_dir']:
-                        misc.sprint("start gridpack!")
-                        # we are in gridpack mode -> create it
-                        me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
-                                                decay_dir), options=options)
-                        me5_cmd.options["automatic_html_opening"] = False
-                        me5_cmd.options["madanalysis5_path"] = None
-                        me5_cmd.options["madanalysis_path"] = None
-                        try:
-                            os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
-                            os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
-                        except Exception,error:
-                            logger.debug(error)
-                            pass
-                        misc.sprint(os.listdir(pjoin(decay_dir,'Cards')))
-                        if self.options["run_card"]:
-                            run_card = self.options["run_card"]
-                        else:
-                            run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))                        
-                        
-                        run_card["iseed"] = self.seed
-                        run_card['gridpack'] = True
-                        run_card['systematics_program'] = False
-                        run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
-                        param_card = self.banner['slha']
-                        open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
-                        self.seed += 1
-                        # actually creation
-                        me5_cmd.exec_cmd("generate_events run_01 -f")
-                        me5_cmd.exec_cmd("exit")                        
-                        #remove pointless informat
-                        misc.call(["rm", "Cards", "bin", 'Source', 'SubProcesses'], cwd=decay_dir)
-                        misc.call(['tar', '-xzpvf', 'run_01_gridpack.tar.gz'], cwd=decay_dir)
-                
-                # Now generate the events
-
-                if not self.options['ms_dir']:
-                    me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
-                                                    decay_dir), options=mg5.options)
-                    me5_cmd.options["automatic_html_opening"] = False
-                    me5_cmd.options["madanalysis5_path"] = None
-                    me5_cmd.options["madanalysis_path"] = None
-                    try:
-                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
-                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
-                    except Exception,error:
-                        logger.debug(error)
-                        pass                    
-                    
-                    if self.options["run_card"]:
-                        run_card = self.options["run_card"]
-                    else:
-                        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["systematics_program"] = 'None'
-                    run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
-                    param_card = self.banner['slha']
-                    open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
-                    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'))            
-                else:
-                    misc.call(['run.sh', str(int(1.2*nb_event)), str(self.seed)], cwd=decay_dir)     
-                    out[i] = lhe_parser.EventFile(pjoin(decay_dir, 'events.lhe.gz'))            
-                if cumul:
-                    break
-            
-            return out
-
+        #    if not enough events. re-generate the missing one.
         
         args = self.split_arg(line)
 
@@ -854,11 +758,11 @@
             if not self.model:
                 modelpath = self.model.get('modelpath+restriction')
                 mg5.exec_cmd("import model %s" % modelpath)      
-            to_event = {}
+            evt_decayfile = {} 
             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)
+                    evt_decayfile[pdg] = self.generate_events(pdg, nb_needed, mg5)
                 elif nb_needed %  nb_event == 0:
                     nb_mult = nb_needed // nb_event
                     part = self.model.get_particle(pdg)
@@ -866,9 +770,9 @@
                     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)
+                        evt_decayfile[pdg] = self.generate_events(pdg, nb_event, mg5)
                     else:
-                        to_event[pdg] = generate_events(pdg, nb_needed, mg5, cumul=True)
+                        evt_decayfile[pdg] = self.generate_events(pdg, nb_needed, mg5, cumul=True)
                 else:
                     part = self.model.get_particle(pdg)
                     name = part.get_name()
@@ -881,7 +785,7 @@
         
         # Compute the branching ratio.
         br = 1
-        for (pdg, event_files) in to_event.items():
+        for (pdg, event_files) in evt_decayfile.items():
             if not event_files:
                 continue
             totwidth = float(self.banner.get('param', 'decay', abs(pdg)).value)
@@ -927,8 +831,8 @@
         
         # 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])
+        for pdg in evt_decayfile:
+            bufferedEvents_decay[pdg] = [{}] * len(evt_decayfile[pdg])
         
         import time
         start = time.time()
@@ -946,25 +850,25 @@
             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:
+                if particle.pdg not in self.final_state or particle.pdg not in evt_decayfile:
                     continue # nothing to do for this particle
                 # check how the decay need to be done
-                nb_decay = len(to_event[particle.pdg])
+                nb_decay = len(evt_decayfile[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 = evt_decayfile[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 = evt_decayfile[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(to_event[particle.pdg][key].cross for key in to_event[particle.pdg])
+                    tot = sum(evt_decayfile[particle.pdg][key].cross for key in evt_decayfile[particle.pdg])
                     r = r * tot
                     cumul = 0
-                    for j,events in to_event[particle.pdg].items():
+                    for j,events in evt_decayfile[particle.pdg].items():
                         cumul += events.cross
                         if r < cumul:
                             decay_file = events
@@ -991,9 +895,9 @@
                             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]
+                                new_file = self.generate_events(particle.pdg, needed, mg5, [decay_file_nb])
+                            evt_decayfile[particle.pdg].update(new_file)
+                            decay_file = evt_decayfile[particle.pdg][decay_file_nb]
                             continue
                         if helicity == decay[0].helicity or helicity==9 or \
                                             self.options["spinmode"] == "none":
@@ -1013,7 +917,7 @@
                 wgts[key] *= self.branching_ratio
             # all particle have been decay if needed
             output_lhe.write(str(event))
-        output_lhe.write('</LesHouchesEvent>\n')        
+        output_lhe.write('</LesHouchesEvents>\n')        
                     
     
     def load_model(self, name, use_mg_default, complex_mass=False):
@@ -1042,6 +946,520 @@
         self.mg5cmd._curr_model = self.model
         self.mg5cmd.process_model()
 
+    def generate_events(self, pdg, nb_event, mg5, restrict_file=None, cumul=False,
+                        output_width=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
+        """
+        
+        if cumul:
+            width = 0.
+        else:   
+            width = 1.
+        part = self.model.get_particle(pdg)
+        if not part:
+            return {}# this particle is not defined in the current model so ignore it
+        name = part.get_name()
+        out = {}
+        logger.info("generate %s decay event for particle %s" % (int(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)
+                if self.options['ms_dir']:
+                    # we are in gridpack mode -> create it
+                    me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
+                                            decay_dir), options=options)
+                    me5_cmd.options["automatic_html_opening"] = False
+                    me5_cmd.options["madanalysis5_path"] = None
+                    me5_cmd.options["madanalysis_path"] = None
+                    me5_cmd.allow_notification_center = False
+                    try:
+                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
+                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
+                    except Exception,error:
+                        logger.debug(error)
+                        pass            
+                    if self.options["run_card"]:
+                        run_card = self.options["run_card"]
+                    else:
+                        run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))                        
+                    
+                    run_card["iseed"] = self.seed
+                    run_card['gridpack'] = True
+                    run_card['systematics_program'] = False
+                    run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
+                    param_card = self.banner['slha']
+                    open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
+                    self.seed += 1
+                    # actually creation
+                    me5_cmd.exec_cmd("generate_events run_01 -f")
+                    if output_width:
+                        if cumul:
+                            width += me5_cmd.results.current['cross']
+                        else:
+                            width *= me5_cmd.results.current['cross']
+                    me5_cmd.exec_cmd("exit")                        
+                    #remove pointless informat
+                    misc.call(["rm", "Cards", "bin", 'Source', 'SubProcesses'], cwd=decay_dir)
+                    misc.call(['tar', '-xzpvf', 'run_01_gridpack.tar.gz'], cwd=decay_dir)
+            
+            # Now generate the events
+
+            if not self.options['ms_dir']:
+                me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
+                                                decay_dir), options=mg5.options)
+                me5_cmd.options["automatic_html_opening"] = False
+                me5_cmd.options["automatic_html_opening"] = False
+                me5_cmd.options["madanalysis5_path"] = None
+                me5_cmd.options["madanalysis_path"] = None
+                me5_cmd.allow_notification_center = False
+                try:
+                    os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
+                    os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
+                except Exception,error:
+                    logger.debug(error)
+                    pass                    
+                if self.options["run_card"]:
+                    run_card = self.options["run_card"]
+                else:
+                    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["systematics_program"] = 'None'
+                run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
+                param_card = self.banner['slha']
+                open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
+                self.seed += 1
+                me5_cmd.exec_cmd("generate_events run_01 -f")
+                if output_width:
+                    if cumul:    
+                        width += me5_cmd.results.current['cross']
+                    else:
+                        width *= me5_cmd.results.current['cross']
+                if run_card["nevents"] > 1.01 * me5_cmd.results.current['nb_event']:
+                    logger.critical('The number of event generated is only %s/%s. This typically indicates that you need specify cut on the decay process.',me5_cmd.results.current['nb_event'], run_card["nevents"])
+                    logger.critical('We strongly suggest that you cancel/discard this run.')
+                me5_cmd.exec_cmd("exit")
+                out[i] = lhe_parser.EventFile(pjoin(decay_dir, "Events", 'run_01', 'unweighted_events.lhe.gz'))            
+            else:
+                misc.call(['run.sh', str(int(1.2*nb_event)), str(self.seed)], cwd=decay_dir)     
+                out[i] = lhe_parser.EventFile(pjoin(decay_dir, 'events.lhe.gz'))            
+            if cumul:
+                break
+        
+        if not output_width:
+            return out
+        else:
+            return out, width
+
+    def run_onshell(self, line):
+        """Run the onshell Algorithm"""
+        
+        # 1. Read the event file to check which decay to perform and the number
+        #   of event to generate for each type of particle. (assume efficiency=1 for spin 0
+        #   otherwise efficiency=2
+        # 2. Generate the associated events
+        # 3. generate the various matrix-element (production/decay/production+decay) 
+        
+        # 3. comput
+        # 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
+
+        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 misc.glob("decay_*_*", self.path_me):
+                shutil.rmtree(name)
+
+        self.events_file.close()
+        orig_lhe = lhe_parser.EventFile(self.events_file.name)
+
+        #count the number of particle need to be decayed.
+        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
+        #misc.sprint(to_decay)
+        #misc.sprint("import the mode -> temporary with logging")   
+        with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
+            mg5 = self.mg5cmd
+            if not self.model:
+                modelpath = self.model.get('modelpath+restriction')
+                mg5.exec_cmd("import model %s" % modelpath)  
+                self.model = mg5._curr_model
+
+        # Handle the banner of the output file
+        if not self.seed:
+            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
+        nevents_for_max = self.options['Nevents_for_max_weigth']
+        if nevents_for_max == 0 :
+            nevents_for_max = 75
+        nevents_for_max *= self.options['max_weight_ps_point']
+        
+        with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
+            mg5 = self.mg5cmd
+            if not self.model:
+                modelpath = self.model.get('modelpath+restriction')
+                mg5.exec_cmd("import model %s" % modelpath)      
+            evt_decayfile = {}
+            br = 1.
+            for pdg, nb_needed in to_decay.items():
+                # muliply by expected effeciency of generation
+                spin = self.model.get_particle(pdg).get('spin')
+                if spin ==1:
+                    efficiency = 1.1
+                else:
+                    efficiency = 2.0
+                    
+                totwidth = self.banner.get('param_card', 'decay', abs(pdg)).value
+                #check if a splitting is needed
+                if nb_needed == nb_event:
+                    nb_needed = int(efficiency*nb_needed) + nevents_for_max   
+                    evt_decayfile[pdg], pwidth = self.generate_events(pdg, nb_needed, mg5, output_width=True)
+                    if pwidth > 1.01*totwidth:
+                        logger.warning('partial width (%s) larger than total width (%s) --from param_card--')
+                    elif pwidth > totwidth:
+                        pwidth = totwidth
+                    br = pwidth / totwidth
+                elif nb_needed %  nb_event == 0:
+                    nb_mult = nb_needed // nb_event
+                    nb_needed = int(efficiency*nb_needed) +nevents_for_max *nb_mult
+                    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:
+                        evt_decayfile[pdg], pwidth = self.generate_events(pdg, nb_event, mg5, output_width=True)
+                        if pwidth > 1.01*totwidth:
+                            logger.warning('partial width (%s) larger than total width (%s) --from param_card--')
+                        elif pwidth > totwidth:
+                            pwidth = totwidth
+                        br = pwidth / totwidth**nb_mult
+                        br *= math.factorial(nb_mult)
+                    else:
+                        evt_decayfile[pdg],pwidth = self.generate_events(pdg, nb_needed, mg5, cumul=True, output_width=True)
+                        if pwidth > 1.01*totwidth:
+                            logger.warning('partial width (%s) larger than total width (%s) --from param_card--')
+                        elif pwidth > totwidth:
+                            pwidth = totwidth
+                        br *= (pwidth / totwidth)**nb_mult
+                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 onshell 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
+        self.efficiency = 1
+        self.cross, self.error = self.banner.get_cross(witherror=True)
+        self.cross *= self.branching_ratio
+        self.error *= self.branching_ratio
+        
+        # 3. generate the various matrix-element
+        self.update_status('generating Madspin matrix element')
+        self.generate_all = madspin.decay_all_events_onshell(self, self.banner, self.events_file, 
+                                                    self.options)
+        self.generate_all.compile()
+        self.all_me = self.generate_all.all_me
+        self.all_f2py = {}
+        
+        #4. determine the maxwgt
+        maxwgt = self.get_maxwgt_for_onshell(orig_lhe, evt_decayfile)
+        
+        
+        #5. generate the decay 
+        orig_lhe.seek(0)
+        output_lhe = lhe_parser.EventFile(orig_lhe.name.replace('.lhe', '_decayed.lhe'), 'w')
+        self.banner.scale_init_cross(self.branching_ratio)
+        self.banner.write(output_lhe, close_tag=False)
+        
+        
+        self.efficiency =1.
+        nb_try, nb_event = 0, len(orig_lhe)
+        
+        start = time.time()
+        for curr_event,production in enumerate(orig_lhe):
+            if curr_event and curr_event % 1000 == 0 and float(str(curr_event)[1:]) ==0:
+                print "decaying event number %s. Efficiency: %s [%s s]" % (curr_event, 1/self.efficiency, time.time()-start)
+            while 1:
+                nb_try +=1
+                decays = self.get_decay_from_file(production, evt_decayfile, nb_event-curr_event)
+                full_evt, wgt = self.get_onshell_evt_and_wgt(production, decays)
+                if random.random()*maxwgt < wgt:
+                    break
+            self.efficiency = curr_event/nb_try
+            # change the weight associate to the event
+            full_evt.wgt *= self.branching_ratio
+            wgts = full_evt.parse_reweight()
+            for key in wgts:
+                wgts[key] *= self.branching_ratio            
+            output_lhe.write(str(full_evt))
+            
+        output_lhe.write('</LesHouchesEvents>\n')    
+        self.efficiency = 1 # to let me5 to write the correct number of events
+#        misc.sprint('Done so far. output written in %s' % output_lhe.name)
+        
+    
+    def get_decay_from_file(self,production, evt_decayfile, nb_remain):
+        """return a dictionary PDG -> list of associated decay"""
+        
+        out = collections.defaultdict(list)
+        particles = [p for p in production if int(p.status) == 1.0]
+        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 evt_decayfile:
+                continue # nothing to do for this particle
+            # check how the decay need to be done
+            nb_decay = len(evt_decayfile[particle.pdg])
+            if nb_decay == 0:
+                continue #nothing to do for this particle
+            # Determine the file to read in order to get the decay [decay_file]
+            if nb_decay == 1:
+                decay_file = evt_decayfile[particle.pdg][0]
+                decay_file_nb = 0
+            elif ids.count(particle.pdg) == nb_decay:
+                decay_file = evt_decayfile[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(evt_decayfile[particle.pdg][key].cross for key in evt_decayfile[particle.pdg])
+                r = r * tot
+                cumul = 0
+                for j,events in evt_decayfile[particle.pdg].items():
+                    
+                    cumul += events.cross
+                    if r < cumul:
+                        decay_file = events
+                        decay_file_nb = j
+                        break
+                    else:
+                        continue
+                else:
+                    raise Exception
+            # So now we know which file to read. Do it and re-generate events for that 
+            # file if needed.
+            while 1:
+                try:
+                    decay = decay_file.next()
+                    break
+                except StopIteration:
+                    eff = self.efficiency
+                    needed = 1.05 * nb_remain / eff
+                    needed = min(50000, max(needed, 1000))
+                    with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
+                        new_file = self.generate_events(particle.pdg, needed, self.mg5cmd, [decay_file_nb])
+                    evt_decayfile[particle.pdg].update(new_file)
+                    decay_file = evt_decayfile[particle.pdg][decay_file_nb]
+                    continue
+            out[particle.pdg].append(decay)
+                        
+        return out
+        
+    
+    def get_maxwgt_for_onshell(self, orig_lhe, evt_decayfile):
+        """determine the maximum weight for the onshell (or similar) strategy"""
+
+        # event_decay is a dict pdg -> list of event file (contain the decay)
+        
+        
+        nevents = self.options['Nevents_for_max_weigth']
+        if nevents == 0 :
+            nevents = 75
+        
+        all_maxwgt = []
+        logger.info("Estimating the maximum weight")
+        logger.info("*****************************")
+        logger.info("Probing the first %s events with %s phase space points" % (nevents, self.options['max_weight_ps_point']))
+
+
+        self.efficiency = 1. / self.options['max_weight_ps_point']
+        start = time.time()
+        for i in range(nevents):
+            if i % 5 ==1:
+                logger.info( "Event %s/%s :  %2fs" % (i, nevents, time.time()-start))
+            maxwgt = 0
+            orig_lhe.seek(0)
+            base_event = orig_lhe.next()
+            for j in range(self.options['max_weight_ps_point']):
+                decays = self.get_decay_from_file(base_event, evt_decayfile, nevents-i)   
+                #carefull base_event is modified by the following function                  
+                _, wgt = self.get_onshell_evt_and_wgt(base_event, decays)
+                maxwgt = max(wgt, maxwgt)
+            all_maxwgt.append(maxwgt)
+            
+        all_maxwgt.sort(reverse=True)
+        assert all_maxwgt[0] >= all_maxwgt[1]
+        decay_tools=madspin.decay_misc()
+        ave_weight, std_weight = decay_tools.get_mean_sd(all_maxwgt)
+        std_weight=math.sqrt(std_weight)
+        base_max_weight = 1.05 * (ave_weight+self.options['nb_sigma']*std_weight)
+
+#        misc.sprint(all_maxwgt)
+        for i in [20, 30, 40, 50]:
+            if len(all_maxwgt) < i:
+                break
+            ave_weight, std_weight = decay_tools.get_mean_sd(all_maxwgt[:i])
+            std_weight=math.sqrt(std_weight)
+            base_max_weight = max(base_max_weight, 1.05 * (ave_weight+self.options['nb_sigma']*std_weight))
+                
+            if all_maxwgt[1] > base_max_weight:
+                base_max_weight = 1.05 * all_maxwgt[1]
+        return base_max_weight
+            
+            
+            
+
+            
+    def get_onshell_evt_and_wgt(self, production, decays):
+        """ return the onshell wgt for the production event associated to the decays
+            return also the full event with decay. 
+            Carefull this modifies production event (pass to the full one)"""
+        
+        tag, order = production.get_tag_and_order()
+        try:
+            info = self.generate_all.all_me[tag]
+        except:
+            misc.sprint(self.generate_all.all_me)
+            raise
+        import copy
+        
+        
+        if hasattr(production, 'me_wgt'):
+            production_me = production.me_wgt
+        else:
+            production_me = self.calculate_matrix_element(production)
+            production.me_wgt = production_me
+            
+        decay_me = 1.0
+        for pdg in decays:
+            for dec in decays[pdg]:
+                decay_me *= self.calculate_matrix_element(dec)
+            random.shuffle(decays[pdg])
+
+        full_event = lhe_parser.Event(str(production))
+        full_event = full_event.add_decays(decays)
+        full_me = self.calculate_matrix_element(full_event)
+        
+#        misc.sprint([p.pdg for p in production])
+#        misc.sprint([p.pdg for p in full_event])
+#        misc.sprint(full_me, production_me, decay_me)
+
+        return full_event, full_me/(production_me*decay_me)
+        
+        
+    def calculate_matrix_element(self, event):
+        """routine to return the matrix element"""        
+        
+        tag, order = event.get_tag_and_order()
+        orig_order = self.all_me[tag]['order']
+        pdir = self.all_me[tag]['pdir']
+        if pdir in self.all_f2py:
+            p = event.get_momenta(orig_order)
+            p = rwgt_interface.ReweightInterface.invert_momenta(p)
+            return self.all_f2py[pdir](p, event.aqcd, 0)
+        else:
+            if sys.path[0] != pjoin(self.path_me, 'madspin_me', 'SubProcesses'):
+                sys.path.insert(0, pjoin(self.path_me, 'madspin_me', 'SubProcesses'))
+            
+            mymod = __import__("%s.matrix2py" % (pdir))
+            reload(mymod)
+            mymod = getattr(mymod, 'matrix2py')  
+            with misc.chdir(pjoin(self.path_me, 'madspin_me', 'SubProcesses', pdir)):
+                with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                    mymod.initialise(pjoin(self.path_me, 'Cards','param_card.dat'))
+            self.all_f2py[pdir] = mymod.get_me  
+            return self.calculate_matrix_element(event)
+        
+        
+    def generate_all_matrix_element(self):
+        
+        # 1. compute the production matrix element -----------------------------
+        processes = [line[9:].strip() for line in self.banner.proc_card
+                     if line.startswith('generate')]
+        processes += [' '.join(line.split()[2:]) for line in self.banner.proc_card
+                      if re.search('^\s*add\s+process', line)]
+        # 2. compute the decay matrix-element
+        decay_text = []
+        processes_decay = []
+        for decays in self.list_branches.values():
+            for decay in  decays:
+                if '=' not in decay:
+                    decay += ' QCD=99'
+                if ',' in decay:
+                    decay_text.append('(%s)' % decay)
+                else:
+                    decay_text.append(decay)
+                processes_decay.append(decay)
+        decay_text = ', '.join(decay_text)
+        processes += []
+        
+        #handle NLO
+        new_processes = []
+        for proc in processes:
+            # deal with @ syntax need to move it after the decay specification
+            if '@' in proc:
+                proc, proc_nb = proc.split('@')
+                try:
+                    proc_nb = int(proc_nb)
+                except ValueError:
+                    raise MadSpinError, 'MadSpin didn\'t allow order restriction after the @ comment: \"%s\" not valid' % proc_nb
+                proc_nb = '@ %i' % proc_nb 
+            else:
+                proc_nb = ''     
+                
+            rwgt_interface.ReweightInterface.get_LO_definition_from_NLO()        
+        
+        
+
+
 if __name__ == '__main__':
     
     a = MadSpinInterface()

=== modified file 'Template/Common/Cards/madspin_card_default.dat'
--- Template/Common/Cards/madspin_card_default.dat	2013-11-13 03:38:29 +0000
+++ Template/Common/Cards/madspin_card_default.dat	2017-03-08 10:27:27 +0000
@@ -7,14 +7,18 @@
 #*    The MadGraph5_aMC@NLO Development Team - Find us at   *
 #*    https://server06.fynu.ucl.ac.be/projects/madgraph     *
 #*                                                          *
+#*    Manual:                                               *
+#*    cp3.irmp.ucl.ac.be/projects/madgraph/wiki/MadSpin     *
+#*                                                          *
 #************************************************************
 #Some options (uncomment to apply)
 #
 # set seed 1
 # set Nevents_for_max_weigth 75 # number of events for the estimate of the max. weight
-# set BW_cut 15                # cut on how far the particle can be off-shell
+# set BW_cut 15                 # cut on how far the particle can be off-shell
+# set spinmode onshell          # Use one of the madspin special mode
  set max_weight_ps_point 400  # number of PS to estimate the maximum for each event
-# 
+
 # specify the decay for the final state particles
 decay t > w+ b, w+ > all all
 decay t~ > w- b~, w- > all all

=== modified file 'Template/LO/Cards/pythia8_card_default.dat'
--- Template/LO/Cards/pythia8_card_default.dat	2016-09-07 21:34:46 +0000
+++ Template/LO/Cards/pythia8_card_default.dat	2017-03-08 10:27:27 +0000
@@ -57,12 +57,18 @@
 Merging:TMS              = -1.0
 ! This must be set manually, according to Pythia8 directives.
 ! An example of possible value is 'pp>LEPTONS,NEUTRINOS'
+! Alternatively, from Pythia v8.223 onwards, the value 'guess' can be 
+! used to instruct Pythia to guess the hard process. The guess would mean 
+! that all particles apart from light partons will be considered as a part 
+! of the hard process. This guess is prone to errors if the desired hard 
+! process is complicated (i.e. contains light partons). The user should
+! then be wary of suspicious error messages in the Pythia log file. 
 Merging:Process          = <set_by_user>
 ! A value of -1 means that it is automatically guessed by MadGraph.
 ! It is however always safer to explicitly set it.
 Merging:nJetMax  		 = -1
 !
-! For all merging schemes, decide wehter you want the merging scale
+! For all merging schemes, decide whehter you want the merging scale
 ! variation computed for only the central weights or all other 
 ! PDF and scale variation weights as well
 SysCalc:fullCutVariation = off

=== modified file 'Template/LO/SubProcesses/makefile'
--- Template/LO/SubProcesses/makefile	2016-11-16 23:59:06 +0000
+++ Template/LO/SubProcesses/makefile	2017-03-08 10:27:27 +0000
@@ -43,8 +43,8 @@
 $(PROG): $(PROCESS) auto_dsig.o $(LIBS)
 	$(FC) -o $(PROG) $(PROCESS) $(LINKLIBS) $(LDFLAGS) $(BIASDEPENDENCIES)
 
-gensym: $(SYMMETRY) configs.inc $(LIBDIR)libmodel.$(libext) $(LIBDIR)libgeneric.$(libext)
-	$(FC) -o gensym $(SYMMETRY) -L../../lib/ -lmodel -lgeneric -lpdf $(LDFLAGS)
+gensym: $(SYMMETRY) configs.inc $(LIBS)
+	$(FC) -o gensym $(SYMMETRY) $(LINKLIBS) $(LDFLAGS) $(BIASDEPENDENCIES)
 
 $(LIBDIR)libmodel.$(libext): ../../Cards/param_card.dat
 	cd ../../Source/MODEL; make

=== modified file 'Template/LO/SubProcesses/myamp.f'
--- Template/LO/SubProcesses/myamp.f	2016-10-03 09:56:24 +0000
+++ Template/LO/SubProcesses/myamp.f	2017-03-08 10:27:27 +0000
@@ -273,9 +273,6 @@
             if (onshell .and. (lbw(nbw).eq. 2) .or.
      $          .not. onshell .and. (lbw(nbw).eq. 1)) then
                cut_bw=.true.
-               if (gForceBW(i, iconfig).eq.1) then
-                  return
-               endif
 c               write(*,*) 'cut_bw: ',nbw,xmass,onshell,lbw(nbw),cut_bw
             endif
          endif

=== modified file 'Template/LO/bin/internal/gen_jpeg-pl'
--- Template/LO/bin/internal/gen_jpeg-pl	2014-03-07 12:07:33 +0000
+++ Template/LO/bin/internal/gen_jpeg-pl	2017-03-08 10:27:27 +0000
@@ -26,7 +26,7 @@
   }
   close(OUT);
   close(IN);
-  system "/bin/bash -c \"nice gs  \-sDEVICE\=jpeg \-sOutputFile\=matrix$imatrix\%00d.jpg \-q \-dNOPAUSE \-dBATCH matrix-1.ps\"";
+  system "/bin/bash -c \"nice gs  \-sDEVICE\=jpeg \-sOutputFile\=matrix$imatrix\%00d.jpg \-q \-dNOPAUSE \-dBATCH matrix-1.ps > /dev/null\"";
   system "rm -f matrix-1.ps";
 
 # Determine how many jpg files we have
@@ -72,7 +72,7 @@
 
     system ("/bin/bash -c \"cat matrix$imatrix.ps | sed 1,352d >> junk.ps\" ");
 
-    system "/bin/bash -c \"nice gs  \-sDEVICE\=jpeg \-sOutputFile\=card.jpg \-q \-dNOPAUSE \-dBATCH \-g180x150 ./junk.ps; rm -f junk.ps; cp -p card.jpg ../../HTML/card.jpg\" ";
+    system "/bin/bash -c \"nice gs  \-sDEVICE\=jpeg \-sOutputFile\=card.jpg \-q \-dNOPAUSE \-dBATCH \-g180x150 ./junk.ps; rm -f junk.ps; cp -p card.jpg ../../HTML/card.jpg > /dev/null\" ";
   }
   if ($imatrix eq "") {$imatrix = 0;}
   $imatrix = $imatrix + 1;

=== modified file 'Template/NLO/Cards/FO_analyse_card.dat'
--- Template/NLO/Cards/FO_analyse_card.dat	2015-03-06 18:09:36 +0000
+++ Template/NLO/Cards/FO_analyse_card.dat	2017-03-08 10:27:27 +0000
@@ -9,7 +9,8 @@
 #                                                                               
 #######################################################################
 #
-# Analysis format. Can either be 'topdrawer', 'root', 'HwU' or 'none'.
+# Analysis format.
+# Can either be 'topdrawer', 'root', 'HwU', 'LHE' 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'
@@ -17,6 +18,7 @@
 # to be set empty.
 FO_ANALYSIS_FORMAT = HwU
 #
+#
 # Needed extra-libraries (FastJet is already linked):
 FO_EXTRALIBS = 
 #

=== modified file 'Template/NLO/FixedOrderAnalysis/HwU.f'
--- Template/NLO/FixedOrderAnalysis/HwU.f	2016-02-18 14:05:45 +0000
+++ Template/NLO/FixedOrderAnalysis/HwU.f	2017-03-08 10:27:27 +0000
@@ -47,11 +47,12 @@
 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
+c     input 3: Same as input 2, but weighted average is same as from MINT
       implicit none
       integer input
       integer error_estimation
       common /HwU_common2/ error_estimation
-      if (input.ge.0 .and. input.le.2) then
+      if (input.ge.0 .and. input.le.3) then
          error_estimation=input
       else
          write (*,*) 'unknown error estimation',input
@@ -68,7 +69,7 @@
 c (approx. 110 MB).
       integer error_estimation
       common /HwU_common2/ error_estimation
-      data error_estimation /1/
+      data error_estimation /3/
       end
       
 c Book the histograms at the start of the run. Give a 'label' (an
@@ -190,19 +191,20 @@
 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)
+      subroutine HwU_accum_iter(inclde,nPSpoints,values)
       implicit none
       include "HwU.inc"
       logical inclde
       integer nPSpoints,label,i,j
       double precision nPSinv,etot,vtot(max_wgts),niter,y_squared
+     $     ,values(2)
       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)
+            call accumulate_results(label,nPSinv,niter,values)
          endif
 c     Reset the histo of the current iteration to zero so that they are
 c     ready for the next iteration.
@@ -227,7 +229,7 @@
       implicit none
       include "HwU.inc"
       integer label,nPSpoints,i,j
-      double precision nPSinv,niter
+      double precision nPSinv,niter,dummy(2)
       nPSinv=1d0/dble(nPSpoints)
       niter=1d0
       do label=1,max_plots
@@ -239,7 +241,7 @@
             enddo
             histy_err(label,i)=0d0
          enddo
-         call accumulate_results(label,nPSinv,niter)
+         call accumulate_results(label,nPSinv,niter,dummy)
       enddo
       return
       end
@@ -254,11 +256,12 @@
 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)
+      subroutine accumulate_results(label,nPSinv,niter,values)
       implicit none
       include "HwU.inc"
       integer label,i,j
       double precision nPSinv,etot,vtot(max_wgts),niter,y_squared
+     $     ,values(2),a1,a2
       integer error_estimation
       common /HwU_common2/ error_estimation
       if (error_estimation.eq.2) then
@@ -281,7 +284,7 @@
                etot=etot* sqrt(dble(histi(label,i))
      &                    /(dble(histi(label,i))-1.5d0))
             else
-               etot=0.999d49
+               etot=abs(vtot(1))*10d0 ! multiply by 10 to make it large
             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
@@ -303,6 +306,47 @@
      $              1d0/sqrt(1d0/histy_err(label,i)**2+1d0/etot**2)
             endif
          enddo
+      elseif (error_estimation.eq.3) then
+         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=abs(vtot(1))*10d0 ! multiply by 10 to make it large
+            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)
+     &                 /values(2)+vtot(j)/values(1))/(1d0
+     &                 /values(2) + 1d0/values(1))
+               enddo
+               a1=((1d0/values(1))/((1d0/values(1))+1d0/values(2)))**2
+               a2=((1d0/values(2))/((1d0/values(1))+1d0/values(2)))**2
+               histy_err(label,i)=sqrt(a2*histy_err(label,i)**2 +
+     &              a1*etot**2)
+            endif
+         enddo 
       elseif(error_estimation.eq.1) then
 c     simply sum the weights in the bins
          do i=1,nbin(label)

=== modified file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_lplm.f'
--- Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_lplm.f	2016-02-17 14:32:27 +0000
+++ Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_lplm.f	2017-03-08 10:27:27 +0000
@@ -136,67 +136,67 @@
       detallb=etal-etalb
 c
       l=0
-      call HwU_fill(l+1,(ptv),(WGTS))
-      call HwU_fill(l+2,(ptv),(WGTS))
-      if(ptv.gt.0.d0)call HwU_fill(l+3,(log10(ptv)),(WGTS))
-      call HwU_fill(l+4,(yv),(WGTS))
-      call HwU_fill(l+5,(etav),(WGTS))
-      call HwU_fill(l+6,(xmv),(WGTS))
-c
-      call HwU_fill(l+7,(ptl),(WGTS))
-      call HwU_fill(l+8,(ptl),(WGTS))
-      if(ptl.gt.0.d0)call HwU_fill(l+9,(log10(ptl)),(WGTS))
-      call HwU_fill(l+10,(etal),(WGTS))
-      call HwU_fill(l+11,(ptlb),(WGTS))
-      call HwU_fill(l+12,(ptlb),(WGTS))
-      if(ptlb.gt.0.d0)call HwU_fill(l+13,(log10(ptlb)),(WGTS))
-      call HwU_fill(l+14,(etalb),(WGTS))
-c
-      call HwU_fill(l+15,(detallb),(WGTS))
-      call HwU_fill(l+16,(azi),(WGTS))
-      if(azinorm.gt.0.d0) call HwU_fill(l+17,(log10(azinorm)),(WGTS))
-      call HwU_fill(l+18,(xmll),(WGTS))
-      call HwU_fill(l+19,(ptpair),(WGTS))
-      if(ptpair.gt.0)call HwU_fill(l+20,(log10(ptpair)),(WGTS))
-      call HwU_fill(l+21,(0d0),(WGTS))
+      call HwU_fill(l+1,(ptv),WGTS)
+      call HwU_fill(l+2,(ptv),WGTS)
+      if(ptv.gt.0.d0)call HwU_fill(l+3,(log10(ptv)),WGTS)
+      call HwU_fill(l+4,(yv),WGTS)
+      call HwU_fill(l+5,(etav),WGTS)
+      call HwU_fill(l+6,(xmv),WGTS)
+c
+      call HwU_fill(l+7,(ptl),WGTS)
+      call HwU_fill(l+8,(ptl),WGTS)
+      if(ptl.gt.0.d0)call HwU_fill(l+9,(log10(ptl)),WGTS)
+      call HwU_fill(l+10,(etal),WGTS)
+      call HwU_fill(l+11,(ptlb),WGTS)
+      call HwU_fill(l+12,(ptlb),WGTS)
+      if(ptlb.gt.0.d0)call HwU_fill(l+13,(log10(ptlb)),WGTS)
+      call HwU_fill(l+14,(etalb),WGTS)
+c
+      call HwU_fill(l+15,(detallb),WGTS)
+      call HwU_fill(l+16,(azi),WGTS)
+      if(azinorm.gt.0.d0) call HwU_fill(l+17,(log10(azinorm)),WGTS)
+      call HwU_fill(l+18,(xmll),WGTS)
+      call HwU_fill(l+19,(ptpair),WGTS)
+      if(ptpair.gt.0)call HwU_fill(l+20,(log10(ptpair)),WGTS)
+      call HwU_fill(l+21,(0d0),WGTS)
 c
       l=l+21
 
       if(abs(etav).lt.ycut)then
-        call HwU_fill(l+1,(ptv),(WGTS))
-        call HwU_fill(l+2,(ptv),(WGTS))
-        if(ptv.gt.0.d0)call HwU_fill(l+3,(log10(ptv)),(WGTS))
+        call HwU_fill(l+1,(ptv),WGTS)
+        call HwU_fill(l+2,(ptv),WGTS)
+        if(ptv.gt.0.d0)call HwU_fill(l+3,(log10(ptv)),WGTS)
       endif
       if(ptv.gt.20.d0)then
-        call HwU_fill(l+4,(yv),(WGTS))
-        call HwU_fill(l+5,(etav),(WGTS))
+        call HwU_fill(l+4,(yv),WGTS)
+        call HwU_fill(l+5,(etav),WGTS)
       endif
       if(abs(etav).lt.ycut.and.ptv.gt.20.d0)then
-         call HwU_fill(l+6,(xmv),(WGTS))
-         call HwU_fill(l+21,(0d0),(WGTS))
+         call HwU_fill(l+6,(xmv),WGTS)
+         call HwU_fill(l+21,(0d0),WGTS)
       endif
 c
       if(abs(etal).lt.ycut)then
-        call HwU_fill(l+7,(ptl),(WGTS))
-        call HwU_fill(l+8,(ptl),(WGTS))
-        if(ptl.gt.0.d0)call HwU_fill(l+9,(log10(ptl)),(WGTS))
+        call HwU_fill(l+7,(ptl),WGTS)
+        call HwU_fill(l+8,(ptl),WGTS)
+        if(ptl.gt.0.d0)call HwU_fill(l+9,(log10(ptl)),WGTS)
       endif
-      if(ptl.gt.20.d0)call HwU_fill(l+10,(etal),(WGTS))
+      if(ptl.gt.20.d0)call HwU_fill(l+10,(etal),WGTS)
       if(abs(etalb).lt.ycut)then
-        call HwU_fill(l+11,(ptlb),(WGTS))
-        call HwU_fill(l+12,(ptlb),(WGTS))
-        if(ptlb.gt.0.d0)call HwU_fill(l+13,(log10(ptlb)),(WGTS))
+        call HwU_fill(l+11,(ptlb),WGTS)
+        call HwU_fill(l+12,(ptlb),WGTS)
+        if(ptlb.gt.0.d0)call HwU_fill(l+13,(log10(ptlb)),WGTS)
       endif
-      if(ptlb.gt.20.d0)call HwU_fill(l+14,(etalb),(WGTS))
+      if(ptlb.gt.20.d0)call HwU_fill(l+14,(etalb),WGTS)
 c
       if( abs(etal).lt.ycut.and.abs(etalb).lt.ycut .and.
      &     ptl.gt.20.d0.and.ptlb.gt.20.d0)then
-        call HwU_fill(l+15,(detallb),(WGTS))
-        call HwU_fill(l+16,(azi),(WGTS))
-        if(azinorm.gt.0.d0) call HwU_fill(l+17,(log10(azinorm)),(WGTS))
-        call HwU_fill(l+18,(xmll),(WGTS))
-        call HwU_fill(l+19,(ptpair),(WGTS))
-        if(ptpair.gt.0) call HwU_fill(l+20,(log10(ptpair)),(WGTS))
+        call HwU_fill(l+15,(detallb),WGTS)
+        call HwU_fill(l+16,(azi),WGTS)
+        if(azinorm.gt.0.d0) call HwU_fill(l+17,(log10(azinorm)),WGTS)
+        call HwU_fill(l+18,(xmll),WGTS)
+        call HwU_fill(l+19,(ptpair),WGTS)
+        if(ptpair.gt.0) call HwU_fill(l+20,(log10(ptpair)),WGTS)
       endif
 
  999  return      

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_V.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_V.f	2016-03-21 08:45:42 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_V.f	2017-03-08 10:27:27 +0000
@@ -125,7 +125,7 @@
       ENDIF
 c
 c CHOOSE IDENT=24 FOR W+, IDENT=-24 FOR W-, IDENT=23 FOR Z0
-      IDENT=23
+      IDENT=24
 C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT''S A POWER-SUPPRESSED
 C EFFECT, SO THROW THE EVENT AWAY
       IF(SIGN(1.D0,PHEP(3,4)).EQ.SIGN(1.D0,PHEP(3,5)))THEN

=== modified file 'Template/NLO/MCatNLO/shower_template.sh'
--- Template/NLO/MCatNLO/shower_template.sh	2016-10-27 14:53:25 +0000
+++ Template/NLO/MCatNLO/shower_template.sh	2017-03-08 10:27:27 +0000
@@ -7,7 +7,7 @@
 RUN_NAME=$3
 NFILE=$4
 
-export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:%(extralibs)s
+export %(ld_library_path)s=$%(ld_library_path)s:%(extralibs)s
 
 # this is for py8 runs
 export PYTHIA8DATA=`pwd`/xmldoc

=== modified file 'Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f'
--- Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f	2016-11-13 13:28:20 +0000
+++ Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f	2017-03-08 10:27:27 +0000
@@ -85,7 +85,6 @@
                 call read_rwgt_line_wgt(iunit,ww(iww))
              enddo
              read(iunit,'(a)')string ! </rwgt>
-             write(*,*) string
           endif
            do while(index(string,'</event>').eq.0)
              read(iunit,'(a)')string ! <rwgt>

=== modified file 'Template/NLO/Source/ranmar.f'
--- Template/NLO/Source/ranmar.f	2014-06-12 08:52:03 +0000
+++ Template/NLO/Source/ranmar.f	2017-03-08 10:27:27 +0000
@@ -1,10 +1,9 @@
       function ran2()
 c Wrapper for the random numbers; needed for the NLO stuff
       implicit none
+      include '../SubProcesses/mint.inc' ! includes iconfig common
       double precision ran2,x,a,b
       integer ii,jconfig
-      integer           iconfig
-      common/to_configs/iconfig
       a=0d0           ! min allowed value for x
       b=1d0           ! max allowed value for x
       ii=0            ! dummy argument of ntuple
@@ -194,6 +193,8 @@
 c     Arguments
 c
       integer iseed
+      integer random_offset_split
+      common /c_random_offset_split/ random_offset_split
 c
 c     Local
 c
@@ -207,7 +208,7 @@
          close(lun)
          return
  14   close(lun)
- 25   iseed = 0
+ 25   iseed = random_offset_split
       end
 
       subroutine ranmar(rvec)

=== modified file 'Template/NLO/Source/run.inc'
--- Template/NLO/Source/run.inc	2016-09-08 23:15:34 +0000
+++ Template/NLO/Source/run.inc	2017-03-08 10:27:27 +0000
@@ -79,4 +79,8 @@
       common/to_rwgt/ do_rwgt_scale, rw_Fscale_down, rw_Fscale_up, rw_Rscale_down, rw_Rscale_up,
      # do_rwgt_pdf, pdf_set_min, pdf_set_max,
      # store_rwgt_info 
-
+c
+c    For FO run (with lhe type of analysis
+c
+      double precision FO_LHE_weight_ratio
+      common /FO_ANALYSIS_LHW/FO_LHE_weight_ratio
\ No newline at end of file

=== modified file 'Template/NLO/SubProcesses/HwU_dummy.f'
--- Template/NLO/SubProcesses/HwU_dummy.f	2015-03-05 16:12:37 +0000
+++ Template/NLO/SubProcesses/HwU_dummy.f	2017-03-08 10:27:27 +0000
@@ -8,7 +8,8 @@
       subroutine HwU_add_points
       end
 
-      subroutine HwU_accum_iter(ldummy,idummy)
+      subroutine HwU_accum_iter(ldummy,idummy,dummy)
       logical ldummy
       integer idummy
+      double precision dummy(2)
       end

=== modified file 'Template/NLO/SubProcesses/ajob_template'
--- Template/NLO/SubProcesses/ajob_template	2015-09-25 14:53:36 +0000
+++ Template/NLO/SubProcesses/ajob_template	2017-03-08 10:27:27 +0000
@@ -34,7 +34,11 @@
 TAGTAGTAGTAGTAGTAGTAG for i in 1 ; do
     
     if [[ $run_mode == 'all' || $run_mode == 'born' ]] ; then
-	j=$run_mode\_G$i
+	if [[ $runnumber == '0' ]] ; then
+	    j=$run_mode\_G$i
+	else
+	    j=$run_mode\_G$i\_$runnumber
+	fi
     else
 	if [[ $runnumber == '0' ]] ; then
 	    j=G$run_mode$i

=== added file 'Template/NLO/SubProcesses/analysis_lhe.f'
--- Template/NLO/SubProcesses/analysis_lhe.f	1970-01-01 00:00:00 +0000
+++ Template/NLO/SubProcesses/analysis_lhe.f	2017-03-08 10:27:27 +0000
@@ -0,0 +1,222 @@
+c LHE analysis routines used when having fo_analysis_format=lhe in the
+c FO_analyse_card.
+c that card also defines FO_LHE_WEIGHT which is a minimal weight to write the 
+c event (event with lower weights will be un-weighted)
+
+      subroutine analysis_begin(nwgt,weights_info)
+      implicit none
+      integer nwgt
+      character*(*) weights_info(*)
+      integer nwgts,nevents
+      double precision sum_of_wgts
+      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
+      logical lopen
+      nwgts=nwgt
+      sum_of_wgts=0d0
+      nevents=0
+      inquire(41,OPENED=lopen) ! in case run is reset due to 2 bad iterations...
+      if (lopen) close(41)
+      open(41,file= 'events.lhe',status='UNKNOWN')
+      write (41,'(a)') '<eventgroup>'
+      end
+
+      subroutine analysis_end(xnorm)
+      implicit none
+      double precision xnorm
+      integer nwgts,nevents
+      double precision sum_of_wgts
+      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
+      logical lopen
+      integer                                   npoints
+      double precision            cross_section
+      common /for_FixedOrder_lhe/ cross_section,npoints
+      inquire(41,OPENED=lopen)  ! safety
+      if (lopen) then
+         backspace(41)          ! overwrite the final </eventgroup> tag
+         write (41,*) nevents,sum_of_wgts,cross_section
+         close(41)
+         open(41, file='header.txt')
+         call write_lhef_header(41, 0, 'FO')
+         close(41)
+      endif
+      end
+
+      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+      implicit none
+      include 'nexternal.inc'
+      include 'reweight.inc'
+      include 'run.inc'
+      integer nwgt,max_weight
+      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
+c      include 'genps.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision p(0:4,nexternal)
+      double precision wgts(max_weight)
+      integer ibody
+      integer nwgts,nevents
+      double precision sum_of_wgts
+      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
+c
+      integer i,j,npart
+      double precision ran2
+      external ran2
+      double precision R,twgt
+      double precision zero
+      integer izero 
+      parameter (zero=0.d0)
+      parameter (izero=0)
+      include 'nFKSconfigs.inc'
+      INTEGER NFKSPROCESS
+      COMMON/C_NFKSPROCESS/NFKSPROCESS
+      integer iSorH_lhe,ifks_lhe(fks_configs) ,jfks_lhe(fks_configs)
+     &     ,fksfather_lhe(fks_configs) ,ipartner_lhe(fks_configs)
+      double precision scale1_lhe(fks_configs),scale2_lhe(fks_configs)
+      common/cto_LHE1/iSorH_lhe,ifks_lhe,jfks_lhe, fksfather_lhe,ipartner_lhe
+      common/cto_LHE2/scale1_lhe,scale2_lhe
+c
+c     forbid unweighting close to the fks pole (should not happen but better safe than sorry)
+      double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
+      double precision xi_i_fks_ev,y_ij_fks_ev
+      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
+c
+c Auxiliary quantities used when writing events
+      integer kwgtinfo
+      integer i_wgt, kk, ii, jj, n, nn
+      character*140 buff
+      INTEGER MAXNUP
+      PARAMETER (MAXNUP=500)
+      INTEGER NUP,IDPRUP,IDUP(MAXNUP),ISTUP(MAXNUP),
+     & MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP)
+      DOUBLE PRECISION XWGTUP,AQEDUP,AQCDUP,
+     & PUP(5,MAXNUP),VTIMUP(MAXNUP),SPINUP(MAXNUP)
+c maximum weight ratio for the partial unweighting
+      integer                                   npoints
+      double precision            cross_section
+      common /for_FixedOrder_lhe/ cross_section,npoints
+c --- do the partial unweighting (don't do it when to close to singular
+c --- region)
+      twgt=abs(cross_section)*FO_LHE_weight_ratio !/npoints
+      if(abs(wgts(1)).lt.abs(twgt) .and. xi_i_fks_ev.gt.1d-3 .and.
+     $     1d0-y_ij_fks_ev.gt.1d-2)then
+         R = ran2()*abs(twgt)
+         if (R.gt.abs(wgts(1)))then
+            return
+         else
+            do i=2,nwgts
+               wgts(i) = wgts(i)*abs(twgt/wgts(1))
+            enddo
+            wgts(1) = sign(twgt,wgts(1))
+         endif
+      endif
+c --- accumulate the total weights of all the events in the event file
+      sum_of_wgts=sum_of_wgts+wgts(1)
+
+c --- fill the multi-weight common blocks with the scale and PDF
+c --- variation weights
+      i_wgt=1
+      if (do_rwgt_scale) then
+         do kk=1,dyn_scale(0)
+            if (lscalevar(kk)) then
+               do ii=1,nint(scalevarF(0))
+                  do jj=1,nint(scalevarR(0))
+                     i_wgt=i_wgt+1
+                     wgtxsecmu(jj,ii,kk)= wgts(i_wgt)
+                  enddo
+               enddo
+            else
+               i_wgt=i_wgt+1
+               wgtxsecmu(1,1,kk)= wgts(i_wgt)
+            endif
+         enddo
+      endif
+      if (do_rwgt_pdf) then
+         do nn=1,lhaPDFid(0)
+            if (lpdfvar(nn)) then
+               do n=0,nmemPDF(nn)
+                  i_wgt=i_wgt+1
+                  wgtxsecPDF(n,nn) = wgts(i_wgt)
+               enddo
+            else
+               i_wgt=i_wgt+1
+               wgtxsecPDF(0,nn) = wgts(i_wgt)
+            endif
+         enddo
+      endif
+
+c --- prepare the event file info according to the LesHouches standard
+      npart = nexternal
+      do i=1,nexternal
+        IDUP(i)= ipdg(i)
+        ISTUP(i)= istatus(i)
+        MOTHUP(1,i)=0
+        MOTHUP(2,i)=0
+        ICOLUP(1,i)=599
+        ICOLUP(2,i)=599
+        PUP(1,i)=p(1,i)
+        PUP(2,i)=p(2,i)
+        PUP(3,i)=p(3,i)
+        PUP(4,i)=p(0,i)
+        PUP(5,i)=p(4,i)
+        VTIMUP(i)=0.d0
+        SPINUP(i)=9
+      enddo
+
+c --- prepare the buffer information
+      if(.not.doreweight)then
+         write(buff,201)'#aMCatNLO',iSorH_lhe,ifks_lhe(nFKSprocess)
+     &        ,jfks_lhe(nFKSprocess),fksfather_lhe(nFKSprocess)
+     &        ,ipartner_lhe(nFKSprocess),scale1_lhe(nFKSprocess)
+     &        ,scale2_lhe(nFKSprocess),izero,izero,izero,zero,zero
+     &        ,zero,zero,ibody*1d0
+      else
+         if(iwgtinfo.ne.-5)then
+            write(*,*)'Error in write_events_lhe'
+            write(*,*)'  Inconsistency in reweight parameters'
+            write(*,*)doreweight,iwgtinfo
+            stop
+         endif
+         kwgtinfo= 9
+         write(buff,201)'#aMCatNLO',iSorH_lhe,ifks_lhe(nFKSprocess)
+     &        ,jfks_lhe(nFKSprocess),fksfather_lhe(nFKSprocess)
+     &        ,ipartner_lhe(nFKSprocess),scale1_lhe(nFKSprocess)
+     &        ,scale2_lhe(nFKSprocess),kwgtinfo,nexternal,iwgtnumpartn
+     &        ,zero,zero,zero,zero,ibody*1d0
+      endif
+
+c --- write the event
+      call write_lhef_event(41,
+     &     npart,IDPRUP,wgts(1),0d0,0d0,0d0,
+     &     IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
+
+ 201  format(a9,1x,i1,4(1x,i2),2(1x,d14.8),2x,i2,2(1x,i2),5(1x,d14.8))
+      end
+
+c This we can use for the event grouping!      
+      subroutine HwU_add_points
+      implicit none
+      integer nwgts,nevents
+      double precision sum_of_wgts
+      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
+      nevents=nevents+1
+      write (41,'(a)') '</eventgroup>'
+      write (41,'(a)') '<eventgroup>'
+      end
+
+c Dummy routines
+      subroutine HwU_accum_iter(ldummy,idummy)
+      logical ldummy
+      integer idummy
+      end
+      subroutine HwU_output(idummy,dummy)
+      integer idummy
+      double precision dummy
+      write (*,*) 'HwU_output should not be called',idummy
+      stop 1
+      end
+      subroutine addfil(string)
+      character*(*) string
+      end
+      subroutine accum(ldummy)
+      logical ldummy
+      end

=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2016-01-22 15:02:25 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2017-03-08 10:27:27 +0000
@@ -32,10 +32,6 @@
       include 'run.inc'
       include 'coupl.inc'
       
-      integer           iconfig
-      common/to_configs/iconfig
-
-
       double precision twgt, maxwgt,swgt(maxevents)
       integer                             lun, nw
       common/to_unwgt/twgt, maxwgt, swgt, lun, nw
@@ -59,9 +55,11 @@
       integer n_mp, n_disc
 c For MINT:
       include "mint.inc"
-      real* 8 xgrid(0:nintervals,ndimmax),ymax(nintervals,ndimmax)
-     $     ,ymax_virt,ans(nintegrals),unc(nintegrals),chi2(nintegrals)
-     $     ,x(ndimmax),itmax_fl
+      integer nhits_in_grids(maxchannels)
+      real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
+     $     ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals
+     $     ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
+     $     ,0:maxchannels),x(ndimmax),itmax_fl
       integer ixi_i,iphi_i,iy_ij,vn
       integer ifold(ndimmax) 
       common /cifold/ifold
@@ -72,13 +70,6 @@
       logical unwgt
       double precision evtsgn
       common /c_unwgt/evtsgn,unwgt
-      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
-     $     ,ndimmax)
-      double precision ave_virt(nintervals_virt,ndimmax)
-     $     ,ave_virt_acc(nintervals_virt,ndimmax)
-     $     ,ave_born_acc(nintervals_virt ,ndimmax)
-      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
-     $     ,nvirt_acc
 
       logical SHsep
       logical Hevents
@@ -90,7 +81,8 @@
 
       double precision virtual_over_born
       common/c_vob/virtual_over_born
-      double precision average_virtual,virtual_fraction
+      double precision average_virtual(maxchannels)
+     $     ,virtual_fraction(maxchannels)
       common/c_avg_virt/average_virtual,virtual_fraction
 
 c timing statistics
@@ -121,8 +113,10 @@
 c     Read general MadFKS parameters
 c
       call FKSParamReader(paramFileName,.TRUE.,.FALSE.)
-      average_virtual=0d0
-      virtual_fraction=max(virt_fraction,min_virt_fraction)
+      do kchan=1,maxchannels
+         average_virtual(kchan)=0d0
+         virtual_fraction(kchan)=max(virt_fraction,min_virt_fraction)
+      enddo
       
       
 c
@@ -164,7 +158,7 @@
 c     Get user input
 c
       write(*,*) "getting user params"
-      call get_user_params(ncall,itmax,iconfig,imode)
+      call get_user_params(ncall,itmax,imode)
       if(imode.eq.0)then
         flat_grid=.true.
       else
@@ -183,7 +177,7 @@
         write(*,*)'NLO computations require muF1=muF2'
         stop
       endif
-      write(*,*) "about to integrate ", ndim,ncall,itmax,iconfig
+      write(*,*) "about to integrate ", ndim,ncall,itmax
 c APPLgrid
       if (imode.eq.0) iappl=0 ! overwrite when starting completely fresh
       if(iappl.ne.0) then
@@ -203,23 +197,27 @@
          if(imode.eq.0)then
 c Don't safe the reweight information when just setting up the grids.
             doreweight=.false.
-            do j=0,nintervals
+            do kchan=1,nchans
                do i=1,ndimmax
-                  xgrid(j,i)=0.d0
+                  do j=0,nintervals
+                     xgrid(j,i,kchan)=0.d0
+                  enddo
                enddo
             enddo
          else
             doreweight=do_rwgt_scale.or.do_rwgt_pdf
 c to restore grids:
             open (unit=12, file='mint_grids',status='old')
-            do j=0,nintervals
-               read (12,*) (xgrid(j,i),i=1,ndim)
-            enddo
-            do j=1,nintervals_virt
-               read (12,*) (ave_virt(j,i),i=1,ndim)
-            enddo
-            read (12,*) ans(1),unc(1),dummy,dummy
-            read (12,*) virtual_fraction,average_virtual
+            do kchan=1,nchans
+               do j=0,nintervals
+                  read (12,*) (xgrid(j,i,kchan),i=1,ndim)
+               enddo
+               do j=1,nintervals_virt
+                  read (12,*) (ave_virt(j,i,kchan),i=1,ndim)
+               enddo
+               read(12,*) ans(1,kchan),unc(1,kchan),dummy,dummy,nhits_in_grids(kchan)
+               read(12,*) virtual_fraction(kchan),average_virtual(kchan)
+            enddo
             close (12)
             write (*,*) "Update iterations and points to",itmax,ncall
          endif
@@ -228,31 +226,38 @@
 
          if (ickkw.eq.-1) then
             min_virt_fraction=1d0
-            virtual_fraction=1d0
+            do kchan=1,nchans
+               virtual_fraction(kchan)=1d0
+            enddo
          endif
 c
 c Setup for parton-level NLO reweighting
          if(do_rwgt_scale.or.do_rwgt_pdf) call setup_fill_rwgt_NLOplot()
-         call mint(sigint,ndim,ncall,itmax,imode,xgrid,ymax,ymax_virt
-     $        ,ans,unc,chi2)
+         call mint(sigint,ndim,ncall,itmax,imode,xgrid,ymax
+     $        ,ymax_virt,ans,unc,chi2,nhits_in_grids)
          call topout
-         write(*,*)'Final result [ABS]:',ans(1),' +/-',unc(1)
-         write(*,*)'Final result:',ans(2),' +/-',unc(2)
-         write(*,*)'chi**2 per D.o.F.:',chi2(1)
+         write(*,*)'Final result [ABS]:',ans(1,0),' +/-',unc(1,0)
+         write(*,*)'Final result:',ans(2,0),' +/-',unc(2,0)
+         write(*,*)'chi**2 per D.o.F.:',chi2(1,0)
          open(unit=58,file='results.dat',status='unknown')
-         write(58,*) ans(1),unc(2),0d0,0,0,0,0,0d0,0d0,ans(2)
+         do kchan=0,nchans
+            write(58,*) ans(1,kchan),unc(2,kchan),0d0,0,0,0,0,0d0,0d0
+     $           ,ans(2,kchan)
+         enddo
          close(58)
 c
 c to save grids:
          open (unit=12, file='mint_grids',status='unknown')
-         do j=0,nintervals
-            write (12,*) (xgrid(j,i),i=1,ndim)
-         enddo
-         do j=1,nintervals_virt
-            write (12,*) (ave_virt(j,i),i=1,ndim)
-         enddo
-         write (12,*) ans(1),unc(1),ncall,itmax
-         write (12,*) virtual_fraction,average_virtual
+         do kchan=1,nchans
+            do j=0,nintervals
+               write (12,*) (xgrid(j,i,kchan),i=1,ndim)
+            enddo
+            do j=1,nintervals_virt
+               write (12,*) (ave_virt(j,i,kchan),i=1,ndim)
+            enddo
+            write (12,*) ans(1,kchan),unc(1,kchan),ncall,itmax,nhits_in_grids(kchan)
+            write (12,*) virtual_fraction(kchan),average_virtual(kchan)
+         enddo
          close (12)
       else
          write (*,*) 'Unknown imode',imode
@@ -317,7 +322,10 @@
       write(*,*) 'Time spent in Total : ',tTot
 
       open (unit=12, file='res.dat',status='unknown')
-      write (12,*)ans(1),unc(1),ans(2),unc(2),itmax,ncall,tTot
+      do kchan=0,nchans
+         write (12,*)ans(1,kchan),unc(1,kchan),ans(2,kchan),unc(2,kchan)
+     $        ,itmax,ncall,tTot
+      enddo
       close(12)
 
 
@@ -378,8 +386,6 @@
       common/tosigint/ndim
       logical       nbody
       common/cnbody/nbody
-      integer           iconfig
-      common/to_configs/iconfig
       double precision p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
      $     ,pswgt_cnt(-2:2),jac_cnt(-2:2)
       common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
@@ -612,7 +618,7 @@
       end
 
 c
-      subroutine get_user_params(ncall,itmax,iconfig,irestart)
+      subroutine get_user_params(ncall,itmax,irestart)
 c**********************************************************************
 c     Routine to get user specified parameters for run
 c**********************************************************************
@@ -625,23 +631,22 @@
       include 'nFKSconfigs.inc'
       include 'fks_info.inc'
       include 'run.inc'
+      include 'mint.inc'
 c
 c     Arguments
 c
-      integer ncall,itmax,iconfig, jconfig
+      integer ncall,itmax
 c
 c     Local
 c
       integer i, j
-      double precision dconfig
+      double precision dconfig(maxchannels)
 c
 c     Global
 c
       integer           isum_hel
       logical                   multi_channel
       common/to_matrix/isum_hel, multi_channel
-      double precision    accuracy
-      common /to_accuracy/accuracy
       integer           use_cut
       common /to_weight/use_cut
 
@@ -679,9 +684,10 @@
       integer mc_hel,ihel
       double precision volh
       common/mc_int2/volh,mc_hel,ihel,fillh
-
+      integer random_offset_split
+      common /c_random_offset_split/ random_offset_split
       logical done
-      character*100 buffer
+      character*140 buffer
 c-----
 c  Begin Code
 c-----
@@ -689,6 +695,7 @@
       unwgt=.false.
       open (unit=83,file='input_app.txt',status='old')
       done=.false.
+      nchans=0
       do while (.not. done)
          read(83,'(a)',err=222,end=222) buffer
          if (buffer(1:7).eq.'NPOINTS') then
@@ -728,19 +735,40 @@
                write(*,*) 'Do MC over helicities for the virtuals'
             endif
             isum_hel=0
+         elseif(buffer(1:6).eq.'NCHANS') then
+            read(buffer(9:),*) nchans
+            write (*,*) 'Number of channels to integrate together:'
+     $           ,nchans
+            if (nchans.gt.maxchannels) then
+               write (*,*) 'Too many integration channels to be '/
+     $              /'integrated together. Increase maxchannels',nchans
+     $              ,maxchannels
+               stop 1
+            endif
          elseif(buffer(1:7).eq.'CHANNEL') then
-            read(buffer(10:),*) dconfig
-            iconfig = int(dconfig)
-            do i=1,mapconfig(0)
-               if (iconfig.eq.mapconfig(i)) then
-                  iconfig=i
-                  exit
-               endif
+            if (nchans.le.0) then
+               write (*,*) '"NCHANS" missing in input files'/
+     $              /' (still zero)',nchans
+               stop
+            endif
+            read(buffer(10:),*) (dconfig(kchan),kchan=1,nchans)
+            do kchan=1,nchans
+               iconfigs(kchan) = int(dconfig(kchan))
+               do i=1,mapconfig(0)
+                  if (iconfigs(kchan).eq.mapconfig(i)) then
+                     iconfigs(kchan)=i
+                     exit
+                  endif
+               enddo
             enddo
-            write(*,12) 'Running Configuration Number: ',iconfig
+            write(*,*) 'Running Configuration Number(s): '
+     $           ,(iconfigs(kchan),kchan=1,nchans)
          elseif(buffer(1:5).eq.'SPLIT') then
-            read(buffer(8:),*) i
-            write (*,*) 'Splitting channel:',i
+            read(buffer(8:),*) random_offset_split
+            write (*,*) 'Splitting channel:',random_offset_split
+         elseif(buffer(1:8).eq.'WGT_MULT') then
+            read(buffer(11:),*) wgt_mult
+            write (*,*) 'Weight multiplier:',wgt_mult
          elseif(buffer(1:8).eq.'RUN_MODE') then
             read(buffer(11:),*) abrvinput
             if(abrvinput(5:5).eq.'0')then
@@ -785,20 +813,9 @@
          endif
       endif
 c
-c
 c     Here I want to set up with B.W. we map and which we don't
 c
-      dconfig = dconfig-iconfig
-      if (dconfig .eq. 0) then
-         write(*,*) 'Not subdividing B.W.'
-         lbw(0)=0
-      else
-         lbw(0)=1
-         jconfig=dconfig*1000.1
-         write(*,*) 'Using dconfig=',jconfig
-         call DeCode(jconfig,lbw(1),3,nexternal)
-         write(*,*) 'BW Setting ', (lbw(j),j=1,nexternal-2)
-      endif
+      lbw(0)=0
  10   format( a)
  12   format( a,i4)
       end

=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2016-01-22 15:02:25 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2017-03-08 10:27:27 +0000
@@ -30,10 +30,8 @@
 cc
       include 'run.inc'
       include 'coupl.inc'
+      include "mint.inc"
       
-      integer           iconfig
-      common/to_configs/iconfig
-
 c Vegas stuff
       common/tosigint/ndim
 
@@ -49,16 +47,16 @@
 
       double precision virtual_over_born
       common/c_vob/virtual_over_born
-      double precision average_virtual,virtual_fraction
+      double precision average_virtual(maxchannels),virtual_fraction(maxchannels)
       common/c_avg_virt/average_virtual,virtual_fraction
 
       double precision weight
 c For MINT:
-      include "mint.inc"
-      real* 8 xgrid(0:nintervals,ndimmax),ymax(nintervals,ndimmax)
-     $     ,ymax_virt,ans(nintegrals),unc(nintegrals),chi2(nintegrals)
-     $     ,x(ndimmax)
-      integer ixi_i,iphi_i,iy_ij,vn
+      real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
+     $     ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals
+     $     ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
+     $     ,0:maxchannels),x(ndimmax)
+      integer ixi_i,iphi_i,iy_ij,vn,nhits_in_grids(maxchannels)
       integer ifold(ndimmax) 
       common /cifold/ifold
       integer ifold_energy,ifold_phi,ifold_yij
@@ -70,13 +68,6 @@
       logical unwgt
       double precision evtsgn
       common /c_unwgt/evtsgn,unwgt
-      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
-     $     ,ndimmax)
-      double precision ave_virt(nintervals_virt,ndimmax)
-     $     ,ave_virt_acc(nintervals_virt,ndimmax)
-     $     ,ave_born_acc(nintervals_virt ,ndimmax)
-      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
-     $     ,nvirt_acc
       double precision ran2
       external ran2
       
@@ -118,7 +109,7 @@
 c
       call FKSParamReader(paramFileName,.TRUE.,.FALSE.)
       average_virtual=0d0
-      virtual_fraction=virt_fraction
+      virtual_fraction(1)=virt_fraction
 
       ntot=0
       nsun=0
@@ -144,7 +135,7 @@
 c     Get user input
 c
       write(*,*) "getting user params"
-      call get_user_params(ncall,itmax,iconfig,imode,
+      call get_user_params(ncall,itmax,imode,
      &     ixi_i,iphi_i,iy_ij,SHsep)
 c Only do the reweighting when actually generating the events
       if (imode.eq.2) then
@@ -189,49 +180,49 @@
 c initialize grids
             do j=0,nintervals
                do i=1,ndimmax
-                  xgrid(j,i)=0.d0
+                  xgrid(j,i,1)=0.d0
                enddo
             enddo
          else
 c to restore grids:
             open (unit=12, file='mint_grids',status='old')
             do j=0,nintervals
-               read (12,*) (xgrid(j,i),i=1,ndim)
+               read (12,*) (xgrid(j,i,1),i=1,ndim)
             enddo
             do j=1,nintervals_virt
-               read (12,*) (ave_virt(j,i),i=1,ndim)
+               read (12,*) (ave_virt(j,i,1),i=1,ndim)
             enddo
-            read (12,*) (ans(i),i=1,nintegrals)
+            read (12,*) (ans(i,1),i=1,nintegrals)
             read (12,*) ifold_energy,ifold_phi,ifold_yij
-            read (12,*) virtual_fraction,average_virtual
+            read (12,*) virtual_fraction(1),average_virtual(1)
             close (12)
          endif
 c
          write (*,*) 'imode is ',imode
          call mint(sigintF,ndim,ncall,itmax,imode,xgrid,ymax,ymax_virt
-     $        ,ans,unc,chi2)
+     $        ,ans,unc,chi2,nhits_in_grids)
          open(unit=58,file='res_0',status='unknown')
-         write(58,*)'Final result [ABS]:',ans(1),' +/-',unc(1)
-         write(58,*)'Final result:',ans(2),' +/-',unc(2)
+         write(58,*)'Final result [ABS]:',ans(1,1),' +/-',unc(1,1)
+         write(58,*)'Final result:',ans(2,1),' +/-',unc(2,1)
          close(58)
-         write(*,*)'Final result [ABS]:',ans(1),' +/-',unc(1)
-         write(*,*)'Final result:',ans(2),' +/-',unc(2)
-         write(*,*)'chi**2 per D.o.F.:',chi2(1)
+         write(*,*)'Final result [ABS]:',ans(1,1),' +/-',unc(1,1)
+         write(*,*)'Final result:',ans(2,1),' +/-',unc(2,1)
+         write(*,*)'chi**2 per D.o.F.:',chi2(1,1)
          open(unit=58,file='results.dat',status='unknown')
-         write(58,*) ans(1),unc(2),0d0,0,0,0,0,0d0,0d0,ans(2)
+         write(58,*) ans(1,1),unc(2,1),0d0,0,0,0,0,0d0,0d0,ans(2,1)
          close(58)
 c
 c to save grids:
          open (unit=12, file='mint_grids',status='unknown')
          do j=0,nintervals
-            write (12,*) (xgrid(j,i),i=1,ndim)
+            write (12,*) (xgrid(j,i,1),i=1,ndim)
          enddo
          do j=1,nintervals_virt
-            write (12,*) (ave_virt(j,i),i=1,ndim)
+            write (12,*) (ave_virt(j,i,1),i=1,ndim)
          enddo
-         write (12,*) (ans(i),i=1,nintegrals)
+         write (12,*) (ans(i,1),i=1,nintegrals)
          write (12,*) ifold_energy,ifold_phi,ifold_yij
-         write (12,*) virtual_fraction,average_virtual
+         write (12,*) virtual_fraction(1),average_virtual(1)
          close (12)
 
 c*************************************************************
@@ -241,14 +232,14 @@
 c to restore grids:
          open (unit=12, file='mint_grids',status='old')
          do j=0,nintervals
-            read (12,*) (xgrid(j,i),i=1,ndim)
+            read (12,*) (xgrid(j,i,1),i=1,ndim)
          enddo
          do j=1,nintervals_virt
-            read (12,*) (ave_virt(j,i),i=1,ndim)
+            read (12,*) (ave_virt(j,i,1),i=1,ndim)
          enddo
-         read (12,*) (ans(i),i=1,nintegrals)
+         read (12,*) (ans(i,1),i=1,nintegrals)
          read (12,*) ifold_energy,ifold_phi,ifold_yij
-         read (12,*) virtual_fraction,average_virtual
+         read (12,*) virtual_fraction(1),average_virtual(1)
          close (12)
 
 c Prepare the MINT folding
@@ -265,45 +256,45 @@
          
          write (*,*) 'imode is ',imode
          call mint(sigintF,ndim,ncall,itmax,imode,xgrid,ymax,ymax_virt
-     $        ,ans,unc,chi2)
+     $        ,ans,unc,chi2,nhits_in_grids)
          
 c If integrating the virtuals alone, we include the virtuals in
 c ans(1). Therefore, no need to have them in ans(5) and we have to set
 c them to zero.
          if (only_virt) then
-            ans(3)=0d0 ! virtual Xsec
-            ans(5)=0d0 ! ABS virtual Xsec
+            ans(3,1)=0d0 ! virtual Xsec
+            ans(5,1)=0d0 ! ABS virtual Xsec
          endif
 
          open(unit=58,file='res_1',status='unknown')
-         write(58,*)'Final result [ABS]:',ans(1)+ans(5),' +/-'
-     $        ,sqrt(unc(1)**2+unc(5)**2)
-         write(58,*)'Final result:',ans(2),' +/-',unc(2)
+         write(58,*)'Final result [ABS]:',ans(1,1)+ans(5,1),' +/-'
+     $        ,sqrt(unc(1,1)**2+unc(5,1)**2)
+         write(58,*)'Final result:',ans(2,1),' +/-',unc(2,1)
          close(58)
-         write(*,*)'Final result [ABS]:',ans(1)+ans(5),' +/-'
-     $        ,sqrt(unc(1)**2+unc(5)**2)
-         write(*,*)'Final result:',ans(2),' +/-',unc(2)
-         write(*,*)'chi**2 per D.o.F.:',chi2(1)
+         write(*,*)'Final result [ABS]:',ans(1,1)+ans(5,1),' +/-'
+     $        ,sqrt(unc(1,1)**2+unc(5,1)**2)
+         write(*,*)'Final result:',ans(2,1),' +/-',unc(2,1)
+         write(*,*)'chi**2 per D.o.F.:',chi2(1,1)
 c write the results.dat file 
          open(unit=58,file='results.dat',status='unknown')
-         write(58,*)ans(1)+ans(5), unc(2), 0d0, 0, 0, 0, 0, 0d0 ,0d0, ans(2) 
+         write(58,*)ans(1,1)+ans(5,1), unc(2,1), 0d0, 0, 0, 0, 0, 0d0 ,0d0, ans(2,1) 
          close(58)
 
 c to save grids:
          open (unit=12, file='mint_grids',status='unknown')
-         write (12,*) (xgrid(0,i),i=1,ndim)
+         write (12,*) (xgrid(0,i,1),i=1,ndim)
          do j=1,nintervals
-            write (12,*) (xgrid(j,i),i=1,ndim)
-            write (12,*) (ymax(j,i),i=1,ndim)
+            write (12,*) (xgrid(j,i,1),i=1,ndim)
+            write (12,*) (ymax(j,i,1),i=1,ndim)
          enddo
          do j=1,nintervals_virt
-            write (12,*) (ave_virt(j,i),i=1,ndim)
+            write (12,*) (ave_virt(j,i,1),i=1,ndim)
          enddo
-         write (12,*) ymax_virt
+         write (12,*) ymax_virt(1)
          write (12,*) (ifold(i),i=1,ndim)
-         write (12,*) (ans(i),i=1,nintegrals)
-         write (12,*) (unc(i),i=1,nintegrals)
-         write (12,*) virtual_fraction,average_virtual
+         write (12,*) (ans(i,1),i=1,nintegrals)
+         write (12,*) (unc(i,1),i=1,nintegrals)
+         write (12,*) virtual_fraction(1),average_virtual(1)
          close (12)
 
 c*************************************************************
@@ -328,23 +319,23 @@
 
 c to restore grids:
          open (unit=12, file='mint_grids',status='unknown')
-         read (12,*) (xgrid(0,i),i=1,ndim)
+         read (12,*) (xgrid(0,i,1),i=1,ndim)
          do j=1,nintervals
-            read (12,*) (xgrid(j,i),i=1,ndim)
-            read (12,*) (ymax(j,i),i=1,ndim)
+            read (12,*) (xgrid(j,i,1),i=1,ndim)
+            read (12,*) (ymax(j,i,1),i=1,ndim)
          enddo
          do j=1,nintervals_virt
-            read (12,*) (ave_virt(j,i),i=1,ndim)
+            read (12,*) (ave_virt(j,i,1),i=1,ndim)
          enddo
-         read (12,*) ymax_virt
+         read (12,*) ymax_virt(1)
          read (12,*) (ifold(i),i=1,ndim)
-         read (12,*) (ans(i),i=1,nintegrals)
-         read (12,*) (unc(i),i=1,nintegrals)
-         read (12,*) virtual_fraction,average_virtual
+         read (12,*) (ans(i,1),i=1,nintegrals)
+         read (12,*) (unc(i,1),i=1,nintegrals)
+         read (12,*) virtual_fraction(1),average_virtual(1)
          close (12)
 
 c determine how many events for the virtual and how many for the no-virt
-         ncall_virt=int(ans(5)/(ans(1)+ans(5)) * ncall)
+         ncall_virt=int(ans(5,1)/(ans(1,1)+ans(5,1)) * ncall)
          ncall_novi=ncall-ncall_virt
 
          write (*,*) "Generating virt :: novi approx.",ncall_virt
@@ -355,11 +346,11 @@
 c fill the information for the write_header_init common block
          ifile=lunlhe
          ievents=ncall
-         inter=ans(2)
-         absint=ans(1)+ans(5)
-         uncer=unc(2)
+         inter=ans(2,1)
+         absint=ans(1,1)+ans(5,1)
+         uncer=unc(2,1)
 
-         weight=(ans(1)+ans(5))/ncall
+         weight=(ans(1,1)+ans(5,1))/ncall
 
          if (abrv(1:3).ne.'all' .and. abrv(1:4).ne.'born' .and.
      $        abrv(1:4).ne.'virt') then
@@ -375,7 +366,7 @@
                vn=3
                call gen(sigintF,ndim,xgrid,ymax,ymax_virt,1,x,vn)
             else
-               if (ran2().lt.ans(5)/(ans(1)+ans(5)) .or. only_virt) then
+               if (ran2().lt.ans(5,1)/(ans(1,1)+ans(5,1)) .or. only_virt) then
                   abrv='virt'
                   if (only_virt) then
                      vn=2
@@ -474,10 +465,10 @@
 
       open (unit=12, file='res.dat',status='unknown')
       if (imode.eq.0) then
-         write (12,*)ans(1),unc(1),ans(2),unc(2),itmax,ncall,tTot
+         write (12,*)ans(1,1),unc(1,1),ans(2,1),unc(2,1),itmax,ncall,tTot
       else
-         write (12,*)ans(1)+ans(5),sqrt(unc(1)**2+unc(5)**2),ans(2)
-     $        ,unc(2),itmax,ncall,tTot
+         write (12,*)ans(1,1)+ans(5,1),sqrt(unc(1,1)**2+unc(5,1)**2),ans(2,1)
+     $        ,unc(2,1),itmax,ncall,tTot
       endif
       close(12)
 
@@ -515,7 +506,7 @@
       end
 
 
-      subroutine get_user_params(ncall,itmax,iconfig,
+      subroutine get_user_params(ncall,itmax,
      &     imode,ixi_i,iphi_i,iy_ij,SHsep)
 c**********************************************************************
 c     Routine to get user specified parameters for run
@@ -533,7 +524,7 @@
 c
 c     Arguments
 c
-      integer ncall,itmax,iconfig, jconfig
+      integer ncall,itmax,jconfig
 c
 c     Local
 c
@@ -667,6 +658,9 @@
          endif
       enddo
       write(*,12) 'Running Configuration Number: ',iconfig
+      nchans=1
+      iconfigs(1)=iconfig
+      wgt_mult=1d0
 
       write (*,'(a)') 'Enter running mode for MINT:'
       write (*,'(a)') '0 to set-up grids, 1 to integrate,'//
@@ -766,8 +760,6 @@
       external passcuts
       parameter (izero=0,ione=1,itwo=2,mohdr=-100)
       data firsttime/.true./
-      integer           iconfig
-      common/to_configs/iconfig
       double precision p_born(0:3,nexternal-1)
       common /pborn/   p_born
       integer     fold

=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2016-10-12 13:52:46 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2017-03-08 10:27:27 +0000
@@ -4404,6 +4404,7 @@
       include "run.inc"
       include "fks_powers.inc"
       include 'reweight.inc'
+      include "mint.inc"
       double precision p(0:3,nexternal),bsv_wgt,born_wgt,avv_wgt
       double precision pp(0:3,nexternal)
       
@@ -4460,7 +4461,7 @@
       integer iminmax
       common/cExceptPSpoint/iminmax,ExceptPSpoint
 
-      double precision average_virtual,virtual_fraction
+      double precision average_virtual(maxchannels),virtual_fraction(maxchannels)
       common/c_avg_virt/average_virtual,virtual_fraction
       double precision virtual_over_born
       common/c_vob/virtual_over_born
@@ -4635,7 +4636,7 @@
 c convert to Binoth Les Houches Accord standards
          virt_wgt=0d0
          if (fold.eq.0) then
-            if ((ran2().le.virtual_fraction .and.
+            if ((ran2().le.virtual_fraction(ichan) .and.
      $           abrv(1:3).ne.'nov').or.abrv(1:4).eq.'virt') then
                call cpu_time(tBefore)
                Call BinothLHA(p_born,born_wgt,virt_wgt)
@@ -4644,9 +4645,9 @@
                tOLP=tOLP+(tAfter-tBefore)
                virtual_over_born=virt_wgt/(born_wgt*ao2pi)
                if (ickkw.ne.-1)
-     &              virt_wgt=virt_wgt-average_virtual*born_wgt*ao2pi
+     &              virt_wgt=virt_wgt-average_virtual(ichan)*born_wgt*ao2pi
                if (abrv.ne.'virt') then
-                  virt_wgt=virt_wgt/virtual_fraction
+                  virt_wgt=virt_wgt/virtual_fraction(ichan)
                endif
                virt_wgt_save=virt_wgt
 c$$$               bsv_wgt=bsv_wgt+virt_wgt_save
@@ -4656,7 +4657,7 @@
 c$$$            bsv_wgt=bsv_wgt+virt_wgt_save
          endif
          if (abrv(1:4).ne.'virt' .and. ickkw.ne.-1)
-     &        avv_wgt=average_virtual*born_wgt*ao2pi
+     &        avv_wgt=average_virtual(ichan)*born_wgt*ao2pi
 
 c eq.(MadFKS.C.13)
          if(abrv.eq.'viSA'.or.abrv.eq.'viSB')then
@@ -5415,9 +5416,11 @@
       end
 
 
-      subroutine setfksfactor(iconfig,match_to_shower)
+      subroutine setfksfactor(match_to_shower)
       implicit none
 
+      include 'mint.inc'
+      
       double precision CA,CF, PI
       parameter (CA=3d0,CF=4d0/3d0)
       parameter (pi=3.1415926535897932385d0)
@@ -5431,7 +5434,7 @@
       logical softtest,colltest
       common/sctests/softtest,colltest
 
-      integer config_fks,i,j,iconfig,fac1,fac2
+      integer config_fks,i,j,fac1,fac2
       logical match_to_shower
 
       double precision fkssymmetryfactor,fkssymmetryfactorBorn,
@@ -5474,6 +5477,8 @@
       common /cxiScut_used/xiScut_used,xiBSVcut_used
       logical rotategranny
       common/crotategranny/rotategranny
+      double precision diagramsymmetryfactor_save(maxchannels)
+      save diagramsymmetryfactor_save
       double precision diagramsymmetryfactor
       common /dsymfactor/diagramsymmetryfactor
 
@@ -5715,6 +5720,7 @@
 c 'FKS_EW' STUFF
       iden_comp=dble(iden_born_FKS(nFKSprocess))/
      &          dble(iden_real_FKS(nFKSprocess))
+
       
       
 c Set matrices used by MC counterterms
@@ -5749,27 +5755,35 @@
 
       if (firsttime) then
 c Check to see if this channel needs to be included in the multi-channeling
-         diagramsymmetryfactor=0d0
+         do kchan=1,nchans
+            diagramsymmetryfactor_save(kchan)=0d0
+         enddo
          if (multi_channel) then
             open (unit=19,file="symfact.dat",status="old",err=14)
             do i=1,mapconfig(0)
                read (19,*,err=23) fac1,fac2
-               if (i.eq.iconfig) then
-                  if (mapconfig(iconfig).ne.fac1) then
-                     write (*,*) 'inconsistency in symfact.dat',i
-     $                    ,iconfig,mapconfig(iconfig),fac1
-                     stop
+               do kchan=1,nchans
+                  if (i.eq.iconfigs(kchan)) then
+                     if (mapconfig(iconfigs(kchan)).ne.fac1) then
+                        write (*,*) 'inconsistency in symfact.dat',i
+     $                       ,kchan,iconfigs(kchan)
+     $                       ,mapconfig(iconfigs(kchan)),fac1
+                        stop
+                     endif
+                     diagramsymmetryfactor_save(kchan)=dble(fac2)
                   endif
-                  diagramsymmetryfactor=dble(fac2)
-               endif
+               enddo
             enddo
             close(19)
          else                   ! no multi_channel
-            diagramsymmetryfactor=1d0
+            do kchan=1,nchans
+               diagramsymmetryfactor_save(kchan)=1d0
+            enddo
          endif
  12      continue
          firsttime=.false.
       endif
+      diagramsymmetryfactor=diagramsymmetryfactor_save(ichan)
 
       return
 
@@ -5781,7 +5795,9 @@
       write (*,*) '"symfact.dat" is not of the correct format'
       stop
  14   continue
-      diagramsymmetryfactor=1d0
+      do kchan=1,nchans
+         diagramsymmetryfactor_save(kchan)=1d0
+      enddo
       goto 12
       end
 

=== modified file 'Template/NLO/SubProcesses/madfks_plot.f'
--- Template/NLO/SubProcesses/madfks_plot.f	2016-11-23 21:35:25 +0000
+++ Template/NLO/SubProcesses/madfks_plot.f	2017-03-08 10:27:27 +0000
@@ -236,14 +236,9 @@
       integer istatus(nexternal),iPDG(nexternal)
       double precision pmass(nexternal)
       common/to_mass/pmass
-      integer maxflow
-      parameter (maxflow=999)
-      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
-     &     icolup(2,nexternal,maxflow),niprocs
-      common /c_leshouche_inc/idup,mothup,icolup,niprocs
-      integer nwgt,max_weight
+      integer max_weight
       parameter (max_weight=maxscales*maxscales+maxpdfs+1)
-      double precision www(max_weight),wgtden,ratio
+      double precision www(max_weight)
       double precision xsecScale_acc(maxscales,maxscales,maxdynscales)
      $     ,xsecPDFr_acc(0:maxPDFs,maxPDFsets)
       common /scale_pdf_print/xsecScale_acc,xsecPDFr_acc
@@ -265,7 +260,7 @@
       chybst=cosh(ybst_til_tolab)
       shybst=sinh(ybst_til_tolab)
       chybstmo=chybst-1.d0
-      do i=nincoming+1,nexternal
+      do i=1,nexternal
          call boostwdir2(chybst,shybst,chybstmo,xd,pp(0,i),pplab(0,i))
       enddo
 c Fill the arrays (momenta, status and PDG):

=== modified file 'Template/NLO/SubProcesses/mint-integrator2.f'
--- Template/NLO/SubProcesses/mint-integrator2.f	2016-09-27 15:07:53 +0000
+++ Template/NLO/SubProcesses/mint-integrator2.f	2017-03-08 10:27:27 +0000
@@ -63,7 +63,7 @@
 c nintegrals=6 : born*alpha_S/2Pi
 c
       subroutine mint(fun,ndim,ncalls0,nitmax,imode,xgrid,ymax,ymax_virt
-     $     ,ans,unc,chi2)
+     $     ,ans,unc,chi2,nhits_in_grids)
 c imode= 0: integrate and adapt the grid
 c imode= 1: frozen grid, compute the integral and the upper bounds
 c imode=-1: same as imode=0, but use previously generated grids
@@ -72,28 +72,33 @@
       include "mint.inc"
       include "FKSParams.inc"
       integer i,j,ncalls0,ndim,nitmax,imode
-      real * 8 fun,xgrid(0:nintervals,ndimmax),xint,ymax(nintervals
-     $     ,ndimmax),ans(nintegrals),unc(nintegrals),ans3(nintegrals,3)
-     $     ,unc3(nintegrals,3),ans_l3(nintegrals),unc_l3(nintegrals)
-     $     ,chi2_l3(nintegrals)
-      real * 8 xint_virt,ymax_virt
+      real * 8 fun,xgrid(0:nintervals,ndimmax,maxchannels)
+     $     ,xint(maxchannels),ymax(nintervals,ndimmax,maxchannels)
+     $     ,ans(nintegrals,0:maxchannels),unc(nintegrals,0:maxchannels)
+     $     ,ans3(nintegrals,3),unc3(nintegrals,3),ans_l3(nintegrals)
+     $     ,unc_l3(nintegrals),chi2_l3(nintegrals),vol_chan
+      real * 8 xint_virt(maxchannels),ymax_virt(maxchannels)
+     $     ,ans_chan(0:maxchannels)
       real * 8 x(ndimmax),vol
-      real * 8 xacc(0:nintervals,ndimmax)
+      real * 8 xacc(0:nintervals,ndimmax,maxchannels)
       integer icell(ndimmax),ncell(ndimmax),ncell_virt
       integer ifold(ndimmax),kfold(ndimmax)
       common/cifold/ifold
-      integer nhits(nintervals,ndimmax)
-      real * 8 rand(ndimmax)
-      real * 8 dx(ndimmax),f(nintegrals),vtot(nintegrals)
-     $     ,etot(nintegrals),prod,f1(nintegrals),chi2(nintegrals)
-     $     ,efrac(nintegrals),dummy
+      integer nhits(nintervals,ndimmax,maxchannels)
+     $     ,nhits_in_grids(maxchannels)
+      real * 8 rand(ndimmax),HwU_values(2)
+      real * 8 dx(ndimmax),f(nintegrals),vtot(nintegrals,0:maxchannels)
+     $     ,etot(nintegrals,0:maxchannels),prod,f1(nintegrals)
+     $     ,chi2(nintegrals,0:maxchannels),efrac(nintegrals),dummy
       integer kdim,kint,kpoint,nit,ncalls,iret,nintcurr,nintcurr_virt
      $     ,ifirst,nit_included,kpoint_iter,non_zero_point(nintegrals)
-     $     ,ntotcalls(nintegrals),nint_used,nint_used_virt,min_it
+     $     ,ntotcalls(nintegrals),nint_used,nint_used_virt,min_it,np
       real * 8 ran3
       external ran3,fun
-      logical even,double_events,bad_iteration
-      double precision average_virtual,virtual_fraction
+      logical even,double_events,bad_iteration,regridded(maxchannels)
+     $     ,reset
+      double precision average_virtual(maxchannels)
+     $     ,virtual_fraction(maxchannels)
       common/c_avg_virt/average_virtual,virtual_fraction
       character*13 title(nintegrals)
       logical new_point
@@ -109,6 +114,11 @@
       common /for_applgrid/ iappl
       logical              fixed_order,nlo_ps
       common /c_fnlo_nlops/fixed_order,nlo_ps
+      integer                                   npoints
+      double precision            cross_section
+      common /for_FixedOrder_lhe/ cross_section,npoints
+      reset=.false.
+
 c if ncalls0 is greater than 0, use the default running, i.e. do not
 c double the events after each iteration as well as use a fixed number
 c of intervals in the grids.
@@ -118,7 +128,7 @@
          nint_used_virt=nintervals_virt
       else
 c if ncalls0.le.0, reset it and double the events per iteration
-         ncalls0=80*ndim
+         ncalls0=80*ndim*(nchans/3+1)
          double_events=.true.
          if (imode.eq.1 .or. imode.eq.-1) then
             nint_used=nintervals
@@ -139,21 +149,46 @@
          do kdim=1,ndim
             ifold(kdim)=1
          enddo
+         ans_chan(0)=0d0
+         do kchan=1,nchans
+            ans_chan(kchan)=ans(1,kchan)
+            ans_chan(0)=ans_chan(0)+ans(1,kchan)
+         enddo
       elseif(imode.eq.0) then
 c Initialize grids
          even=.true.
          min_it=min_it0
          do kdim=1,ndim
             ifold(kdim)=1
-            do kint=0,nint_used
-               xgrid(kint,kdim)=dble(kint)/nint_used
+            do kchan=1,nchans
+               do kint=0,nint_used
+                  xgrid(kint,kdim,kchan)=dble(kint)/nint_used
+                  nhits(kint,kdim,kchan)=0
+               enddo
+               regridded(kchan)=.true.
             enddo
          enddo
+         do kchan=1,maxchannels
+            nhits_in_grids(kchan)=0
+         enddo
          call init_ave_virt(nint_used_virt,ndim)
+         do kchan=0,nchans
+            ans_chan(kchan)=0d0
+         enddo
+         if (double_events) then
+c           when double events, start with the very first channel
+c           only. For the first iteration, we compute each channel
+c           separately.
+            ans_chan(0)=1d0
+            ans_chan(1)=1d0
+            ncalls0=ncalls0/nchans
+         endif
       elseif(imode.eq.1) then
 c Initialize upper bounding envelope
-         xint=ans(1)
-         xint_virt=ans(5)
+         do kchan=1,nchans
+            xint(kchan)=ans(1,kchan)
+            xint_virt(kchan)=ans(5,kchan)
+         enddo
          even=.false.
          min_it=min_it1
          do kdim=1,ndim
@@ -166,36 +201,59 @@
      $              ,nint_used_virt
                stop
             endif
-            do kint=1,nintcurr
-               ymax(kint,kdim)=xint**(1d0/ndim)
+            do kchan=1,nchans
+               do kint=1,nintcurr
+                  ymax(kint,kdim,kchan)=xint(kchan)**(1d0/ndim)
+               enddo
             enddo
          enddo
-         ymax_virt=xint_virt
+         do kchan=1,nchans
+            ymax_virt(kchan)=xint_virt(kchan)
+         enddo
+         ans_chan(0)=0d0
+         do kchan=1,nchans
+            ans_chan(kchan)=ans(1,kchan)
+            ans_chan(0)=ans_chan(0)+ans(1,kchan)
+         enddo
       endif
+      cross_section=ans_chan(0)
       nit=0
       nit_included=0
       do i=1,nintegrals
-         ans(i)=0d0
-         unc(i)=0d0
+         do kchan=0,nchans
+            ans(i,kchan)=0d0
+            unc(i,kchan)=0d0
+         enddo
          do j=1,3
             ans3(i,j)=0d0
             unc3(i,j)=0d0
          enddo
       enddo
+      do i=1,2
+         HwU_values(i)=0d0
+      enddo
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c Main loop over the iterations
  10   continue
-c This makes sure that current stdout is written to log.txt and cache is
-c emptied.
+      do kchan=1,nchans
+         np=0
+         do kint=1,nint_used
+            np=np+nhits(kint,1,kchan)
+         enddo
+         write (*,250) 'channel',kchan,':',iconfigs(kchan)
+     $        ,regridded(kchan),np,nhits_in_grids(kchan),ans_chan(kchan)
+     $        ,ans(2,kchan),virtual_fraction(kchan)
+      enddo
       call flush(6)
       if(nit.ge.nitmax) then
 c We did enough iterations, update arguments and return
-         if(imode.eq.0) xint=ans(1)
-         if(imode.eq.0) xint_virt=ans(5)
-         if (nit_included.ge.2) then
-            chi2(1)=chi2(1)/dble(nit_included-1)
-         else
-            chi2(1)=0d0
-         endif
+         do kchan=0,nchans
+            if (nit_included.ge.2) then
+               chi2(1,kchan)=chi2(1,kchan)/dble(nit_included-1)
+            else
+               chi2(1,kchan)=0d0
+            endif
+         enddo
          write (*,*) '-------'
          ncalls0=ncalls*kpoint_iter ! return number of points used
          if (double_events) then
@@ -203,6 +261,17 @@
          else
             nitmax=nit_included
          endif
+         cross_section=ans(2,0)
+         do kchan=1,nchans
+            if (regridded(kchan)) then
+               np=0
+               do kint=1,nint_used
+                  np=np+nhits(kint,1,kchan)
+               enddo
+               nhits_in_grids(kchan)=np ! set equal to number of points
+                                        ! used for last update
+            endif
+         enddo
          return
       endif
       nit=nit+1
@@ -217,29 +286,51 @@
          ncalls=ncalls0
          write (*,*) 'Update # PS points: ',ncalls0,' --> ',ncalls
       endif
+      npoints=ncalls
 c Reset the accumulated results for grid updating
       if(imode.eq.0) then
-         do kdim=1,ndim
-            do kint=0,nint_used
-               xacc(kint,kdim)=0
-               if(kint.gt.0) then
-                  nhits(kint,kdim)=0
+         do kchan=1,nchans
+            if (regridded(kchan).or.reset) then ! only reset if grids were 
+                                                ! updated (or forced reset)
+               if (regridded(kchan) .and. .not. reset) then
+                  np=0
+                  do kint=1,nint_used
+                     np=np+nhits(kint,1,kchan)
+                  enddo
+                  nhits_in_grids(kchan)=np ! set equal to number of
+                                           ! points used for last update
+               elseif (regridded(kchan) .and. reset) then
+                  nhits_in_grids(kchan)=0
                endif
-            enddo
+               do kdim=1,ndim
+                  do kint=0,nint_used
+                     xacc(kint,kdim,kchan)=0
+                     if(kint.gt.0) then
+                        nhits(kint,kdim,kchan)=0
+                     endif
+                  enddo
+               enddo
+            endif
          enddo
+         reset=.false.
       endif
-      do i=1,nintegrals
-         vtot(i)=0
-         etot(i)=0
+      do kchan=0,nchans
+         do i=1,nintegrals
+            vtot(i,kchan)=0
+            etot(i,kchan)=0
+         enddo
       enddo
       kpoint_iter=0
       do i=1,nintegrals
          non_zero_point(i)=0
       enddo
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c Loop over PS points
  2    kpoint_iter=kpoint_iter+1
       do kpoint=1,ncalls
          new_point=.true.
+c first determine which integration channel we want to pick this point from         
+         call get_channel(ans_chan,vol_chan)
 c find random x, and its random cell
          do kdim=1,ndim
             kfold(kdim)=1
@@ -260,16 +351,17 @@
          enddo
          ifirst=0
  1       continue
-         vol=1
+         vol=1d0/vol_chan * wgt_mult
 c c convert 'flat x' ('rand') to 'vegas x' ('x') and include jacobian ('vol')
          do kdim=1,ndim
             nintcurr=nint_used/ifold(kdim)
             icell(kdim)=ncell(kdim)+(kfold(kdim)-1)*nintcurr
-            dx(kdim)=xgrid(icell(kdim),kdim)-xgrid(icell(kdim)-1,kdim)
+            dx(kdim)=xgrid(icell(kdim),kdim,ichan)-xgrid(icell(kdim)-1
+     $           ,kdim,ichan)
             vol=vol*dx(kdim)*nintcurr
-            x(kdim)=xgrid(icell(kdim)-1,kdim)+rand(kdim)*dx(kdim)
-            if(imode.eq.0)
-     &           nhits(icell(kdim),kdim)=nhits(icell(kdim),kdim)+1
+            x(kdim)=xgrid(icell(kdim)-1,kdim,ichan)+rand(kdim)*dx(kdim)
+            if(imode.eq.0) nhits(icell(kdim),kdim,ichan)=
+     $           nhits(icell(kdim),kdim,ichan)+1
          enddo
          call get_ave_virt(x,nint_used_virt,ndim,average_virtual)
 c contribution to integral
@@ -297,13 +389,14 @@
          if(imode.eq.0) then
 c accumulate the function in xacc(icell(kdim),kdim) to adjust the grid later
             do kdim=1,ndim
-               xacc(icell(kdim),kdim)=xacc(icell(kdim),kdim)+f(1)
+               xacc(icell(kdim),kdim,ichan)= 
+     $              xacc(icell(kdim),kdim,ichan) + f(1)
             enddo
 c Set the Born contribution (to compute the average_virtual) to zero if
 c the virtual was not computed for this phase-space point. Compensate by
 c including the virtual_fraction.
             if (f(3).ne.0d0) then
-               f(6)=f(6)/virtual_fraction
+               f(6)=f(6)/virtual_fraction(ichan)
                call fill_ave_virt(x,nint_used_virt,ndim,f(3),f(6))
             else
                f(6)=0d0
@@ -312,7 +405,7 @@
 c update the upper bounding envelope total rate
             prod=1d0
             do kdim=1,ndim
-               prod=prod*ymax(ncell(kdim),kdim)
+               prod=prod*ymax(ncell(kdim),kdim,ichan)
             enddo
             prod=(f(1)/prod)
             if (prod.gt.1d0) then
@@ -324,14 +417,16 @@
                prod=min(2d0,prod)
                prod=prod**(1d0/dble(ndim))
                do kdim=1,ndim
-                  ymax(ncell(kdim),kdim)=ymax(ncell(kdim),kdim)*prod
+                  ymax(ncell(kdim),kdim,ichan)= 
+     $                 ymax(ncell(kdim),kdim,ichan)*prod
                enddo
             endif
 c Update the upper bounding envelope virtual. Do not include the
 c enhancement due to the virtual_fraction. (And again limit by factor 2
 c at most).
-            if (f(5)*virtual_fraction.gt.ymax_virt) ymax_virt=min(f(5)
-     $           *virtual_fraction,ymax_virt*2d0)
+            if (f(5)*virtual_fraction(ichan).gt.ymax_virt(ichan))
+     $           ymax_virt(ichan)=min(f(5)*virtual_fraction(ichan)
+     $           ,ymax_virt(ichan)*2d0)
 c for consistent printing in the log files (in particular when doing LO
 c runs), set also f(6) to zero when imode.eq.1 and the virtuals are not
 c included.
@@ -342,8 +437,8 @@
          enddo
 c Add the PS point to the result of this iteration
          do i=1,nintegrals
-            vtot(i)=vtot(i)+f(i)
-            etot(i)=etot(i)+f(i)**2
+            vtot(i,ichan)=vtot(i,ichan)+f(i)
+            etot(i,ichan)=etot(i,ichan)+f(i)**2
          enddo
          if (f(1).ne.0d0) call HwU_add_points
       enddo
@@ -365,32 +460,74 @@
 c that pass cuts.
       if (non_zero_point(1).lt.int(0.99*ncalls)
      &                        .and. double_events) goto 2
-
+c This is the loop over all the channels for the very first iteration.
+      if (imode.eq.0 .and. nit.eq.1 .and. double_events) then
+         do kchan=nchans,1,-1
+            if (ans_chan(kchan).eq.1d0) then
+               do i=1,nintegrals
+                  vtot(i,kchan)=vtot(i,kchan)/dble(ntotcalls(i))
+                  etot(i,kchan)=etot(i,kchan)/dble(ntotcalls(i))
+                  etot(i,kchan)=sqrt(abs(etot(i,kchan)-vtot(i,kchan)**2)
+     $                 /dble(ntotcalls(i)))
+               enddo
+               if (kchan.eq.nchans) exit
+               ans_chan(kchan)=0d0
+               ans_chan(kchan+1)=1d0
+               do i=1,nintegrals
+                  ntotcalls(i)=0
+                  non_zero_point(i)=0
+               enddo
+               kpoint_iter=0
+               goto 2
+            endif
+         enddo
+c set the total result for the first iteration as the sum over all the channels
+         do i=1,nintegrals
+            do kchan=1,nchans
+               vtot(i,0)=vtot(i,0)+vtot(i,kchan)
+               etot(i,0)=etot(i,0)+etot(i,kchan)**2
+            enddo
+            etot(i,0)=sqrt(etot(i,0))
+         enddo
+         ncalls0=ncalls0*nchans
+      endif
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c Iteration done. Update the accumulated results and print them to the
 c screen
       do i=1,nintegrals
-         vtot(i)=vtot(i)/dble(ntotcalls(i))
-         etot(i)=etot(i)/dble(ntotcalls(i))
+         if (imode.eq.0 .and. nit.eq.1 .and. double_events) then
+            continue
+         else
+            do kchan=nchans,0,-1
+               if (kchan.ne.0) then
+                  vtot(i,0)=vtot(i,0)+vtot(i,kchan)
+                  etot(i,0)=etot(i,0)+etot(i,kchan)
+               endif
+               vtot(i,kchan)=vtot(i,kchan)/dble(ntotcalls(i))
+               etot(i,kchan)=etot(i,kchan)/dble(ntotcalls(i))
 c the abs is to avoid tiny negative values
-         etot(i)=sqrt(abs(etot(i)-vtot(i)**2)
-     $        /dble(ntotcalls(i)))
-         if (vtot(i).ne.0d0) then
-            efrac(i)=abs(etot(i)/vtot(i))
+               etot(i,kchan)=sqrt(abs(etot(i,kchan)-vtot(i,kchan)**2)
+     $              /dble(ntotcalls(i)))
+            enddo
+         endif
+         if (vtot(i,0).ne.0d0) then
+            efrac(i)=abs(etot(i,0)/vtot(i,0))
          else
             efrac(i)=0d0
          endif
       enddo
       do i=1,nintegrals
          write(*,'(a,1x,e10.4,1x,a,1x,e10.4,1x,a,1x,f7.3,1x,a)')
-     $        title(i)//' =',vtot(i),' +/- ',etot(i),' (',efrac(i)*100d0
-     $        ,'%)'
+     $        title(i)//' =',vtot(i,0),' +/- ',etot(i,0),' (',efrac(i)
+     $        *100d0 ,'%)'
       enddo
 C If there was a large fluctation in this iteration, be careful with
 C including it in the accumalated results and plots.
-      if (efrac(1).gt.0.3d0 .and. iappl.eq.0) then
+      if (efrac(1).gt.0.3d0 .and. iappl.eq.0 .and. nit.gt.3) then
 c Do not include the results in the plots
          if (fixed_order) call accum(.false.)
-         if (fixed_order) call HwU_accum_iter(.false.,ntotcalls(1))
+         if (fixed_order) call HwU_accum_iter(.false.,ntotcalls(1)
+     $        ,HwU_values)
       endif
       if (efrac(1).gt.0.3d0 .and. nit.gt.3 .and. iappl.eq.0) then
 c Do not include the results in the updating of the grids.
@@ -401,7 +538,7 @@
 c empty the accumalated results for the MINT grids
          if (imode.eq.0) then
 c emptying accum. results is done above when the iteration starts
-            continue
+            reset=.true.
          elseif (imode.eq.1) then
 c Cannot really skip the increase of the upper bounding envelope. So,
 c simply continue here. Note that no matter how large the integrand for
@@ -423,33 +560,56 @@
             nit_included=0
 c Reset the MINT grids
             if (imode.eq.0) then
-               do kdim=1,ndim
-                  do kint=0,nint_used
-                     xgrid(kint,kdim)=dble(kint)/nint_used
+               do kchan=1,nchans
+                  do kdim=1,ndim
+                     do kint=0,nint_used
+                        xgrid(kint,kdim,kchan)=dble(kint)/nint_used
+                     enddo
+                     regridded(kchan)=.true.
                   enddo
                enddo
+               do kchan=1,maxchannels
+                  nhits_in_grids(kchan)=0
+               enddo
                call init_ave_virt(nint_used_virt,ndim)
+               do kchan=0,nchans
+                  ans_chan(kchan)=0d0
+               enddo
+               if (double_events) then
+                  ans_chan(0)=1d0
+                  ans_chan(1)=1d0
+                  ncalls0=ncalls0/nchans
+               endif
             elseif (imode.eq.1) then
                do kdim=1,ndim
                   nintcurr=nint_used/ifold(kdim)
-                  do kint=1,nintcurr
-                     ymax(kint,kdim)=xint**(1d0/ndim)
+                  do kchan=1,nchans
+                     do kint=1,nintcurr
+                        ymax(kint,kdim,kchan)=xint(kchan)**(1d0/ndim)
+                     enddo
                   enddo
                   nintcurr_virt=nint_used_virt/ifold(kdim)
                enddo
-               ymax_virt=xint_virt
+               do kchan=1,nchans
+                  ymax_virt(kchan)=xint_virt(kchan)
+               enddo
             endif
             call reset_MC_grid  ! reset the grid for the integers
             if (fixed_order) call initplot  ! Also reset all the plots
             do i=1,nintegrals
-               ans(i)=0d0
-               unc(i)=0d0
-               chi2(i)=0d0
+               do kchan=0,nchans
+                  ans(i,kchan)=0d0
+                  unc(i,kchan)=0d0
+                  chi2(i,kchan)=0d0
+               enddo
                do j=1,3
                   ans3(i,j)=0d0
                   unc3(i,j)=0d0
                enddo
             enddo
+            do i=1,2
+               HwU_values(i)=0d0
+            enddo
             bad_iteration=.false.
          else
             bad_iteration=.true.
@@ -458,81 +618,82 @@
       else
          bad_iteration=.false.
       endif
+      HwU_values(1)=etot(1,0)
+      HwU_values(2)=unc(1,0)
       if(nit.eq.1) then
-         do i=1,nintegrals
-            ans(i)=vtot(i)
-            unc(i)=etot(i)
+         do kchan=0,nchans
+            do i=1,nintegrals
+               ans(i,kchan)=vtot(i,kchan)
+               unc(i,kchan)=etot(i,kchan)
+            enddo
+            ans_chan(kchan)=ans(1,kchan)
          enddo
          write (*,'(a,1x,e10.4)') 'Chi^2 per d.o.f.',0d0
       else
 c prevent annoying division by zero for nearly zero
 c integrands
+         do kchan=nchans,0,-1 ! go backwards so that kchan=0 goes last
+                              ! (this makes sure central value is
+                              ! correctly updated).
          do i=1,nintegrals
-            if(etot(i).eq.0.and.unc(i).eq.0) then
-               if(ans(i).eq.vtot(i) .and. i.eq.1) then
+            if(etot(i,0).eq.0.and.unc(i,0).eq.0) then
+               if(ans(i,0).eq.vtot(i,0) .and. i.eq.1) then
 c double the number of points for the next iteration
                   if (double_events) ncalls0=ncalls0*2
                   goto 10
                else
-                  unc(i)=abs(vtot(i)-ans(i))
-                  etot(i)=abs(vtot(i)-ans(i))
+                  unc(i,kchan)=abs(vtot(i,kchan)-ans(i,kchan))
+                  etot(i,kchan)=abs(vtot(i,kchan)-ans(i,kchan))
                endif
-            elseif(etot(i).eq.0) then
-               etot(i)=unc(i)
-            elseif(unc(i).eq.0) then ! 1st iteration; set to a large value
-               unc(i)=etot(i)*1d99
+            elseif(etot(i,0).eq.0) then
+               etot(i,kchan)=unc(i,kchan)
+            elseif(unc(i,0).eq.0) then ! 1st iteration; set to a large value
+               unc(i,kchan)=etot(i,kchan)*1d99
             endif
-            if (i.ne.1 .and. (etot(i).eq.0 .or. unc(i).eq.0)) then
-               ans(i)=0d0
-               unc(i)=0d0
-               chi2(i)=0d0
+            if (i.ne.1 .and. (etot(i,0).eq.0 .or. unc(i,0).eq.0)) then
+               ans(i,kchan)=0d0
+               unc(i,kchan)=0d0
+               chi2(i,kchan)=0d0
             else
-               ans(i)=(ans(i)/unc(i)+vtot(i)/etot(i))/
-     &              (1/unc(i)+1/etot(i))
-               unc(i)=1/sqrt(1/unc(i)**2+1/etot(i)**2)
-               chi2(i)=chi2(i)+(vtot(i)-ans(i))**2/etot(i)**2
+               ans(i,kchan)=(ans(i,kchan)/unc(i,0)+vtot(i,kchan)/etot(i
+     $              ,0))/(1/unc(i,0)+1/etot(i,0))
+               unc(i,kchan)=1/sqrt(1/unc(i,kchan)**2+1/etot(i,kchan)**2)
+               chi2(i,kchan)=chi2(i,kchan)+(vtot(i,kchan)-ans(i,kchan))
+     $              **2/etot(i,kchan)**2
             endif
          enddo
-         write (*,'(a,1x,e10.4)') 'Chi^2=',(vtot(1)-ans(1))**2
-     $        /etot(1)**2
+         ans_chan(kchan)=ans(1,kchan)
+         enddo
+         write (*,'(a,1x,e10.4)') 'Chi^2=',(vtot(1,0)-ans(1,0))**2
+     $        /etot(1,0)**2
       endif
-
       nit_included=nit_included+1
       do i=1,nintegrals
-         if (ans(i).ne.0d0) then
-            efrac(i)=abs(unc(i)/ans(i))
+         if (ans(i,0).ne.0d0) then
+            efrac(i)=abs(unc(i,0)/ans(i,0))
          else
             efrac(i)=0d0
          endif
          write(*,'(a,1x,e10.4,1x,a,1x,e10.4,1x,a,1x,f7.3,1x,a)')
-     $        'accumulated results '//title(i)//' =',ans(i),' +/- '
-     $        ,unc(i) ,' (',efrac(i)*100d0,'%)'
+     $        'accumulated results '//title(i)//' =',ans(i,0),' +/- '
+     $        ,unc(i,0) ,' (',efrac(i)*100d0,'%)'
       enddo
+      cross_section=ans(1,0)
       if (nit_included.le.1) then
          write (*,'(a,1x,e10.4)') 'accumulated result Chi^2 per DoF ='
      $        ,0d0
       else
          write (*,'(a,1x,e10.4)') 'accumulated result Chi^2 per DoF =',
-     &        chi2(1)/dble(nit_included-1)
+     &        chi2(1,0)/dble(nit_included-1)
       endif
       if (imode.eq.0) then
-c     Update the average_virtual: a_new=(virt+a_old*born)/born
-         if (vtot(6).ne.0d0) then
-            if (average_virtual.eq.0d0) then ! i.e. first iteration
-               average_virtual=vtot(3)/vtot(6)+average_virtual
-            else  ! give some importance to the iterations already done
-               average_virtual=(vtot(3)/vtot(6)+average_virtual*2d0)/2d0
-            endif
-         endif
 c Update the fraction of the events for which we include the virtual corrections
 c in the calculation
-         virtual_fraction=max(min(virtual_fraction*max(min(2d0*etot(3)
-     $        /etot(1),2d0),0.25d0),1d0),Min_virt_fraction)
-         write (*,'(a,1x,f7.3,1x,f7.3)') 'update virtual fraction to:'
-     $        ,virtual_fraction,average_virtual
-      elseif (imode.eq.1) then
-         write (*,'(a,1x,f7.3,1x,f7.3)') 'virtual fraction is:'
-     $        ,virtual_fraction,average_virtual
+         do kchan=1,nchans
+            virtual_fraction(kchan)=max(min(virtual_fraction(kchan)
+     $           *max(min(2d0*etot(3,kchan)/etot(1,kchan),2d0),0.25d0)
+     $           ,1d0),Min_virt_fraction)
+         enddo
       endif
 c Update the results of the last tree iterations
       do i=1,nintegrals
@@ -540,8 +701,8 @@
             ans3(i,j)=ans3(i,j+1)
             unc3(i,j)=unc3(i,j+1)
          enddo
-         ans3(i,3)=vtot(i)
-         unc3(i,3)=etot(i)
+         ans3(i,3)=vtot(i,0)
+         unc3(i,3)=etot(i,0)
       enddo
 c Compute the results of the last three iterations
       if (nit_included.ge.4) then
@@ -582,9 +743,12 @@
       endif
       if(imode.eq.0) then
 c Iteration is finished; now rearrange the grid
-         do kdim=1,ndim
-            call regrid(xacc(0,kdim),xgrid(0,kdim),nhits(1,kdim)
-     $           ,nint_used)
+         do kchan=1,nchans
+            do kdim=1,ndim
+               call regrid(xacc(0,kdim,kchan),xgrid(0,kdim,kchan)
+     $              ,nhits(1,kdim,kchan),nint_used,nhits_in_grids(kchan)
+     $              ,regridded(kchan))
+            enddo
          enddo
          call regrid_ave_virt(nint_used_virt,ndim)
 c Regrid the MC over integers (used for the MC over FKS dirs)
@@ -592,13 +756,14 @@
       endif
 c Quit if the desired accuracy has been reached
       if (nit_included.ge.min_it .and. accuracy.gt.0d0) then
-         if (unc(1)/ans(1)*max(1d0,chi2(1)/dble(nit_included-1))
+         if (unc(1,0)/ans(1,0)*max(1d0,chi2(1,0)/dble(nit_included-1))
      $        .lt.accuracy) then
             write (*,*) 'Found desired accuracy'
             nit=nitmax
 c Improve the stats in the plots
             if (fixed_order) call accum(.true.)
-            if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1))
+            if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1)
+     $           ,HwU_values)
             goto 10
          elseif(unc_l3(1)/ans_l3(1)*max(1d0,chi2_l3(1)).lt.accuracy)
      $           then
@@ -606,20 +771,24 @@
      &           'Found desired accuracy in last 3 iterations'
             nit=nitmax
             do i=1,nintegrals
-               ans(i)=ans_l3(i)
-               unc(i)=unc_l3(i)
-               chi2(i)=chi2_l3(i)*dble(nit_included-1)
+               ans(i,0)=ans_l3(i)
+               unc(i,0)=unc_l3(i)
+               chi2(i,0)=chi2_l3(i)*dble(nit_included-1)
             enddo
 c Improve the stats in the plots
             if (fixed_order) call accum(.true.)
-            if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1))
+            if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1)
+     $           ,HwU_values)
             goto 10
          endif
       endif
 c Double the number of intervals in the grids if not yet reach the maximum
       if (2*nint_used.le.nintervals .and. double_events) then
-         do kdim=1,ndim
-            call double_grid(xgrid(0,kdim),nint_used)
+         do kchan=1,nchans
+            do kdim=1,ndim
+               call double_grid(xgrid(0,kdim,kchan),nhits(1,kdim,kchan)
+     $              ,nint_used)
+            enddo
          enddo
          nint_used=2*nint_used
       endif
@@ -632,33 +801,49 @@
       if (double_events) ncalls0=ncalls0*2
 c Also improve stats in plots
       if (fixed_order) call accum(.true.)
-      if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1))
+      if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1)
+     $     ,HwU_values)
 c Do next iteration
       goto 10
+      return
+ 250  format(a7,i5,1x,a1,1x,i5,1x,l,1x,i8,1x,i8,2x,e10.4,2x,e10.4,2x
+     $     ,e10.4)
       end
 
-      subroutine double_grid(xgrid,ninter)
+      subroutine double_grid(xgrid,nhits,ninter)
       implicit none
       include "mint.inc"
-      integer  ninter
+      integer  ninter,nhits(nintervals)
       real * 8 xgrid(0:nintervals)
       integer i
       do i=ninter,1,-1
          xgrid(i*2)=xgrid(i)
          xgrid(i*2-1)=(xgrid(i)+xgrid(i-1))/2d0
+         nhits(i*2)=nhits(i)/2
+         nhits(i*2-1)=nhits(i)-nhits(i*2)
       enddo
       return
       end
 
 
-      subroutine regrid(xacc,xgrid,nhits,ninter)
+      subroutine regrid(xacc,xgrid,nhits,ninter,nhits_in_grid,regridded)
       implicit none
       include "mint.inc"
-      integer  ninter,nhits(nintervals)
+      logical regridded
+      integer  ninter,nhits(nintervals),np,nhits_in_grid
       real * 8 xacc(0:nintervals),xgrid(0:nintervals)
       real * 8 xn(nintervals),r,tiny,xl,xu,nl,nu,sum
       parameter ( tiny=1d-8 )
       integer kint,jint
+c compute total number of points and update grids if large
+      np=0
+      do kint=1,ninter
+         np=np+nhits(kint)
+      enddo
+      regridded=.false.
+      if (np.lt.nhits_in_grid) return
+      regridded=.true.
+      
 c Use the same smoothing as in VEGAS uses for the grids, i.e. use the
 c average of the central and the two neighbouring grid points: (Only do
 c this if we are already at the maximum intervals, because the doubling
@@ -797,33 +982,39 @@
       implicit none
       integer ndim,imode
       include "mint.inc"
-      real * 8 fun,xgrid(0:nintervals,ndimmax),ymax(nintervals,ndimmax)
-     $     ,ymax_virt,x(ndimmax)
-      real * 8 dx(ndimmax),xx(ndimmax)
+      real * 8 fun,xgrid(0:nintervals,ndimmax,maxchannels)
+     $     ,ymax(nintervals,ndimmax,maxchannels),ymax_virt(maxchannels)
+     $     ,x(ndimmax)
+      real * 8 dx(ndimmax),xx(ndimmax),vol_chan,dummy
       integer icell(ndimmax),ncell(ndimmax),ncell_virt
       integer ifold(ndimmax),kfold(ndimmax)
       common/cifold/ifold
       real * 8 r,f(nintegrals),f1(nintegrals),ubound,vol,ran3
-     $     ,xmmm(nintervals,ndimmax),dummy
+     $     ,xmmm(nintervals,ndimmax,maxchannels)
       real * 8 rand(ndimmax)
       external fun,ran3
       integer icalls,mcalls,kdim,kint,nintcurr,nintcurr_virt,iret,ifirst
      $     ,i,vn,icalls_virt,mcalls_virt,icalls_nz,icalls_virt_nz
-      double precision average_virtual,virtual_fraction
+      double precision average_virtual(maxchannels)
+     $     ,virtual_fraction(maxchannels)
       common/c_avg_virt/average_virtual,virtual_fraction
       save icalls,mcalls,icalls_virt,mcalls_virt,xmmm,icalls_nz
      $     ,icalls_virt_nz
       logical new_point
       common /c_new_point/ new_point
       if(imode.eq.0) then
-         do kdim=1,ndim
-            nintcurr=nintervals/ifold(kdim)
-            xmmm(1,kdim)=ymax(1,kdim)
-            do kint=2,nintcurr
-               xmmm(kint,kdim)=xmmm(kint-1,kdim)+ymax(kint,kdim)
-            enddo
-            do kint=1,nintcurr
-               xmmm(kint,kdim)=xmmm(kint,kdim)/xmmm(nintcurr,kdim)
+         do kchan=1,nchans
+            do kdim=1,ndim
+               nintcurr=nintervals/ifold(kdim)
+               xmmm(1,kdim,kchan)=ymax(1,kdim,kchan)
+               do kint=2,nintcurr
+                  xmmm(kint,kdim,kchan)=xmmm(kint-1,kdim,kchan)
+     $                 +ymax(kint,kdim,kchan)
+               enddo
+               do kint=1,nintcurr
+                  xmmm(kint,kdim,kchan)=xmmm(kint,kdim,kchan)
+     $                 /xmmm(nintcurr,kdim,kchan)
+               enddo
             enddo
          enddo
          icalls=0
@@ -868,6 +1059,7 @@
          stop
       endif
  10   continue
+      call get_channel(ymax_virt,vol_chan)
       new_point=.true.
       if (vn.eq.1) then
          icalls_virt=icalls_virt+1
@@ -885,7 +1077,7 @@
             nintcurr=nintervals/ifold(kdim)
             r=ran3(.false.)
             do kint=1,nintcurr
-               if(r.lt.xmmm(kint,kdim)) then
+               if(r.lt.xmmm(kint,kdim,ichan)) then
                   ncell(kdim)=kint
                   exit
                endif
@@ -896,7 +1088,7 @@
       if (vn.eq.2 .or. vn.eq.3) then
          ubound=1
          do kdim=1,ndim
-            ubound=ubound*ymax(ncell(kdim),kdim)
+            ubound=ubound*ymax(ncell(kdim),kdim,ichan)
          enddo
       endif
       do kdim=1,ndim
@@ -907,18 +1099,19 @@
       enddo
       ifirst=0
  5    continue
-      vol=1
+      vol=1d0/vol_chan
       do kdim=1,ndim
          nintcurr=nintervals/ifold(kdim)
          nintcurr_virt=nintervals_virt/ifold(kdim)
          icell(kdim)=ncell(kdim)+(kfold(kdim)-1)*nintcurr
-         dx(kdim)=xgrid(icell(kdim),kdim)-xgrid(icell(kdim)-1,kdim)
+         dx(kdim)=xgrid(icell(kdim),kdim,ichan)-xgrid(icell(kdim)-1,kdim
+     $        ,ichan)
          vol=vol*dx(kdim)*nintervals/ifold(kdim)
-         x(kdim)=xgrid(icell(kdim)-1,kdim)+rand(kdim)*dx(kdim)
+         x(kdim)=xgrid(icell(kdim)-1,kdim,ichan)+rand(kdim)*dx(kdim)
       enddo
       call get_ave_virt(x,nintcurr_virt,ndim,average_virtual)
       if (vn.eq.1) then
-         ubound=ymax_virt
+         ubound=ymax_virt(ichan)
       endif
       dummy=fun(x,vol,ifirst,f1)
       do i=1,nintegrals
@@ -1108,20 +1301,15 @@
       implicit none
       include "mint.inc"
       integer kdim,ndim,ninter,i
-      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
-     $     ,ndimmax)
-      double precision ave_virt(nintervals_virt,ndimmax)
-     $     ,ave_virt_acc(nintervals_virt,ndimmax)
-     $     ,ave_born_acc(nintervals_virt ,ndimmax)
-      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
-     $     ,nvirt_acc
-      do kdim=1,ndim
-         do i=1,ninter
-            nvirt(i,kdim)=0
-            ave_virt(i,kdim)=0d0
-            nvirt_acc(i,kdim)=0
-            ave_virt_acc(i,kdim)=0d0
-            ave_born_acc(i,kdim)=0d0
+      do kchan=1,nchans
+         do kdim=1,ndim
+            do i=1,ninter
+               nvirt(i,kdim,kchan)=0
+               ave_virt(i,kdim,kchan)=0d0
+               nvirt_acc(i,kdim,kchan)=0
+               ave_virt_acc(i,kdim,kchan)=0d0
+               ave_born_acc(i,kdim,kchan)=0d0
+            enddo
          enddo
       enddo
       return
@@ -1131,20 +1319,14 @@
       implicit none
       include "mint.inc"
       integer kdim,ndim,ninter,ncell
-      double precision x(ndimmax),average_virtual
-      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
-     $     ,ndimmax)
-      double precision ave_virt(nintervals_virt,ndimmax)
-     $     ,ave_virt_acc(nintervals_virt,ndimmax)
-     $     ,ave_born_acc(nintervals_virt ,ndimmax)
-      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
-     $     ,nvirt_acc
-      average_virtual=0d0
+      double precision x(ndimmax),average_virtual(maxchannels)
+      average_virtual(ichan)=0d0
       do kdim=1,ndim
          ncell=min(int(x(kdim)*ninter)+1,ninter)
-         average_virtual=average_virtual+ave_virt(ncell,kdim)
+         average_virtual(ichan)=average_virtual(ichan)+
+     &        ave_virt(ncell,kdim,ichan)
       enddo
-      average_virtual=average_virtual/ndim
+      average_virtual(ichan)=average_virtual(ichan)/ndim
       return
       end
 
@@ -1153,18 +1335,13 @@
       include "mint.inc"
       integer kdim,ndim,ninter,ncell
       double precision x(ndimmax),virtual,born
-      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
-     $     ,ndimmax)
-      double precision ave_virt(nintervals_virt,ndimmax)
-     $     ,ave_virt_acc(nintervals_virt,ndimmax)
-     $     ,ave_born_acc(nintervals_virt ,ndimmax)
-      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
-     $     ,nvirt_acc
       do kdim=1,ndim
          ncell=min(int(x(kdim)*ninter)+1,ninter)
-         nvirt_acc(ncell,kdim)=nvirt_acc(ncell,kdim)+1
-         ave_virt_acc(ncell,kdim)=ave_virt_acc(ncell,kdim)+virtual
-         ave_born_acc(ncell,kdim)=ave_born_acc(ncell,kdim)+born
+         nvirt_acc(ncell,kdim,ichan)=nvirt_acc(ncell,kdim,ichan)+1
+         ave_virt_acc(ncell,kdim,ichan)=ave_virt_acc(ncell,kdim,ichan)
+     $        +virtual
+         ave_born_acc(ncell,kdim,ichan)=ave_born_acc(ncell,kdim,ichan)
+     $        +born
       enddo
       return
       end
@@ -1173,35 +1350,38 @@
       implicit none
       include "mint.inc"
       integer ninter,ndim,kdim,i
-      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
-     $     ,ndimmax)
-      double precision ave_virt(nintervals_virt,ndimmax)
-     $     ,ave_virt_acc(nintervals_virt,ndimmax)
-     $     ,ave_born_acc(nintervals_virt,ndimmax)
-      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
-     $     ,nvirt_acc
-c need to solve for k_new = (virt+k_old*born)/born
+c need to solve for k_new = (virt+k_old*born)/born = virt/born + k_old
+      do kchan=1,nchans
       do kdim=1,ndim
          do i=1,ninter
-            if (ave_born_acc(i,kdim).eq.0d0) cycle
-            if (ave_virt(i,kdim).eq.0d0) then ! i.e. first iteration
-               ave_virt(i,kdim)= ave_virt_acc(i,kdim)/ave_born_acc(i
-     $              ,kdim)+ave_virt(i,kdim)
-            else  ! give some importance to the iterations already done
-               ave_virt(i,kdim)=(ave_virt_acc(i,kdim)/ave_born_acc(i
-     $              ,kdim)+ave_virt(i,kdim)*2d0)/2d0
+            if (ave_born_acc(i,kdim,kchan).eq.0d0) cycle
+            if (ave_virt(i,kdim,kchan).eq.0d0) then ! i.e. first iteration
+               ave_virt(i,kdim,kchan)= ave_virt_acc(i,kdim,kchan)
+     $              /ave_born_acc(i,kdim,kchan)+ave_virt(i,kdim,kchan)
+            else  ! give some importance to the iterations already
+                  ! done: simply base it on the number of points already
+                  ! computed in that bin. Give half the weight to old iterations
+                  ! k_new = [N_new*(V/B+k_old) + (N_old/2)*k_old] / [N_new+(N_old/2)]
+               ave_virt(i,kdim,kchan)=(nvirt_acc(i,kdim,kchan)
+     $              *ave_virt_acc(i,kdim,kchan)/ave_born_acc(i,kdim
+     $              ,kchan)+(nvirt_acc(i,kdim,kchan)+nint(nvirt(i,kdim
+     $              ,kchan)/2d0))*ave_virt(i,kdim,kchan))
+     $              /dble(nvirt_acc(i,kdim,kchan)+nint(nvirt(i,kdim
+     $              ,kchan)/2d0))
             endif
          enddo
       enddo
 c reset the acc values
       do kdim=1,ndim
          do i=1,ninter
-            nvirt(i,kdim)=nvirt(i,kdim)+nvirt_acc(i,kdim)
-            nvirt_acc(i,kdim)=0
-            ave_born_acc(i,kdim)=0d0
-            ave_virt_acc(i,kdim)=0d0
+            nvirt(i,kdim,kchan)=nvirt(i,kdim,kchan)+nvirt_acc(i,kdim
+     $           ,kchan)
+            nvirt_acc(i,kdim,kchan)=0
+            ave_born_acc(i,kdim,kchan)=0d0
+            ave_virt_acc(i,kdim,kchan)=0d0
          enddo
       enddo
+      enddo
       return
       end
 
@@ -1210,33 +1390,83 @@
       implicit none
       include "mint.inc"
       integer kdim,ndim,i,ninter
-      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
-     $     ,ndimmax)
-      double precision ave_virt(nintervals_virt,ndimmax)
-     $     ,ave_virt_acc(nintervals_virt,ndimmax)
-     $     ,ave_born_acc(nintervals_virt ,ndimmax)
-      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
-     $     ,nvirt_acc
+      do kchan=1,nchans
       do kdim=1,ndim
          do i=ninter,1,-1
-            ave_virt(i*2,kdim)=ave_virt(i,kdim)
-            if (nvirt(i,kdim).ne.0) then
-               nvirt(i*2,kdim)=max(nvirt(i,kdim)/2,1)
+            ave_virt(i*2,kdim,kchan)=ave_virt(i,kdim,kchan)
+            if (nvirt(i,kdim,kchan).ne.0) then
+               nvirt(i*2,kdim,kchan)=max(nvirt(i,kdim,kchan)/2,1)
             else
-               nvirt(i*2,kdim)=0
+               nvirt(i*2,kdim,kchan)=0
             endif
             if (i.ne.1) then
-               ave_virt(i*2-1,kdim)=(ave_virt(i,kdim)
-     $              +ave_virt(i-1,kdim))/2d0
-               if (nvirt(i,kdim)+nvirt(i-1,kdim).ne.0) then
-                  nvirt(i*2-1,kdim)=
-     &                 max((nvirt(i,kdim)+nvirt(i-1,kdim))/4,1)
-               else
-                  nvirt(i*2-1,kdim)=0
+               ave_virt(i*2-1,kdim,kchan)=(ave_virt(i,kdim,kchan)
+     $              +ave_virt(i-1,kdim,kchan))/2d0
+               if (nvirt(i,kdim,kchan)+nvirt(i-1,kdim,kchan).ne.0) then
+                  nvirt(i*2-1,kdim,kchan)=
+     &              max((nvirt(i,kdim,kchan)+nvirt(i-1,kdim,kchan))/4,1)
+               else
+                  nvirt(i*2-1,kdim,kchan)=0
+               endif
+            else
+               if (nvirt(1,kdim,kchan).ne.0) then
+                  nvirt(1,kdim,kchan)=max(nvirt(1,kdim,kchan)/2,1)
+               else
+                  nvirt(1,kdim,kchan)=0
                endif
             endif
          enddo
       enddo
-      return
-      end
-
+      enddo
+      return
+      end
+
+
+
+      subroutine get_channel(ans,vol_chan)
+c Picks one random 'ichan' among the 'nchans' integration channels and
+c fills the channels common block in mint.inc.
+      implicit none
+      include 'mint.inc'
+      double precision ans(0:maxchannels),vol_chan,ran3,target,sum
+      external ran3
+      if (nchans.eq.1) then
+         ichan=1
+         iconfig=iconfigs(ichan)
+         vol_chan=1d0
+      elseif (nchans.gt.1) then
+         if (ans(0).le.0d0) then
+c     pick one at random (flat)
+            ichan=int(ran3(.false.)*nchans)+1
+            iconfig=iconfigs(ichan)
+            vol_chan=1d0/dble(nchans)
+         else
+c     pick one at random (weighted by cross section)
+            sum=0d0
+            do kchan=1,nchans
+               sum=sum+ans(kchan)
+            enddo
+            if (abs(sum-ans(0))/(sum+ans(0)).gt.1d-8) then
+               write (*,*) 'ERROR: sum should be equal to ans',
+     &              sum,ans(0)
+               stop 1
+            endif
+            target=ans(0)*ran3(.false.)
+            sum=0d0
+            ichan=0
+            do while (sum.lt.target)
+               ichan=ichan+1
+               sum=sum+ans(ichan)
+            enddo
+            if (ichan.eq.0 .or. ichan.gt.nchans) then
+               write (*,*) 'ERROR: ichan cannot be zero or'/
+     $              /' larger than nchans',ichan,nchans
+               stop
+            endif
+            iconfig=iconfigs(ichan)
+            vol_chan=ans(ichan)/ans(0)
+         endif
+      endif
+      return
+      end
+      

=== modified file 'Template/NLO/SubProcesses/mint.inc'
--- Template/NLO/SubProcesses/mint.inc	2016-09-28 15:25:10 +0000
+++ Template/NLO/SubProcesses/mint.inc	2017-03-08 10:27:27 +0000
@@ -1,3 +1,11 @@
+*     -*-fortran-*-
+
+      integer maxchannels
+      parameter (maxchannels=20) ! same as in amcatnlo_run_interface
+      integer nchans,iconfigs(maxchannels),kchan,iconfig,ichan
+      double precision wgt_mult
+      common /channels/ wgt_mult,nchans,iconfigs,iconfig,ichan
+
 c Note that the number of intervals in the integration
 c grids, 'nintervals', should be equal to
 c "nintervals = min_inter * 2^n",
@@ -16,3 +24,11 @@
 
       double precision accuracy
       common /to_accuracy/ accuracy
+
+      integer nvirt(nintervals_virt,ndimmax,maxchannels)
+     $     ,nvirt_acc(nintervals_virt,ndimmax,maxchannels)
+      double precision ave_virt(nintervals_virt,ndimmax,maxchannels)
+     $     ,ave_virt_acc(nintervals_virt,ndimmax,maxchannels)
+     $     ,ave_born_acc(nintervals_virt ,ndimmax,maxchannels)
+      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
+     $     ,nvirt_acc

=== modified file 'Template/NLO/SubProcesses/montecarlocounter.f'
--- Template/NLO/SubProcesses/montecarlocounter.f	2016-10-25 07:08:18 +0000
+++ Template/NLO/SubProcesses/montecarlocounter.f	2017-03-08 10:27:27 +0000
@@ -203,7 +203,7 @@
      #      ( abs(particle_type(ipart)).ne.3 .and.
      #        particle_type(ipart).ne.8 ) )then
           write(*,*)'Error #2 in check_mc_matrices',i,ipart,
-     #              particle_type(ipart)
+     #particle_type(ipart)
           stop
         endif
       enddo
@@ -2501,6 +2501,8 @@
       double precision shattmp,dot,emsca_bare,ref_scale,scalemin,
      &scalemax,rrnd,ran2,emscainv,dum(5),xm12,qMC,ptresc
       integer ileg
+      double precision p_born(0:3,nexternal-1)
+      common/pborn/p_born
 
       logical emscasharp
       double precision emsca
@@ -2521,7 +2523,7 @@
      &                       xm12,dum(1),dum(2),dum(3),dum(4),dum(5),qMC,.true.)
 
       emsca=2d0*sqrt(ebeam(1)*ebeam(2))
-      scalemax=sqrt((1-xi_i_fks)*shat)
+      call assign_ref_scale(p_born,xi_i_fks,shat,scalemax)
       if(dampMCsubt)then
          call assign_scaleminmax(shat,xi_i_fks,scalemin,scalemax,ileg,xm12)
          emscasharp=(scalemax-scalemin).lt.(1d-3*scalemax)
@@ -2545,14 +2547,17 @@
 
       subroutine assign_scaleminmax(shat,xi,xscalemin,xscalemax,ileg,xm12)
       implicit none
+      include "nexternal.inc"
       include "run.inc"
       include "madfks_mcatnlo.inc"
       integer i,ileg
       double precision shat,xi,ref_scale,xscalemax,xscalemin,xm12
       character*4 abrv
       common/to_abrv/abrv
+      double precision p_born(0:3,nexternal-1)
+      common/pborn/p_born
 
-      ref_scale=sqrt((1-xi)*shat)
+      call assign_ref_scale(p_born,xi,shat,ref_scale)
       xscalemin=max(shower_scale_factor*frac_low*ref_scale,scaleMClow)
       xscalemax=max(shower_scale_factor*frac_upp*ref_scale,
      &              xscalemin+scaleMCdelta)
@@ -2568,6 +2573,33 @@
       end
 
 
+      subroutine assign_ref_scale(p,xii,sh,ref_sc)
+      implicit none
+      include "nexternal.inc"
+      double precision p(0:3,nexternal-1),xii,sh,ref_sc
+      integer i_scale,i
+      parameter(i_scale=1)
+
+      ref_sc=0d0
+      if(i_scale.eq.0)then
+c Born-level CM energy squared
+         ref_sc=dsqrt(max(0d0,(1-xii)*sh))
+      elseif(i_scale.eq.1)then
+c Sum of final-state transverse masses
+         do i=3,nexternal-1
+            ref_sc=ref_sc+dsqrt(max(0d0,(p(0,i)+p(3,i))*(p(0,i)-p(3,i))))
+         enddo
+         ref_sc=ref_sc/2d0
+      else
+         write(*,*)'Wrong i_scale in assign_ref_scale',i_scale
+         stop
+      endif
+c Safety threshold for the reference scale
+      ref_sc=max(ref_sc,30d0)
+
+      return
+      end
+
 
       subroutine dinvariants_dFKS(ileg,s,x,yi,yj,xm12,xm22,dw1dx,dw1dy,dw2dx,dw2dy)
 c Returns derivatives of Mandelstam invariants with respect to FKS variables

=== modified file 'Template/NLO/SubProcesses/setcuts.f'
--- Template/NLO/SubProcesses/setcuts.f	2016-11-14 10:30:46 +0000
+++ Template/NLO/SubProcesses/setcuts.f	2017-03-08 10:27:27 +0000
@@ -144,6 +144,7 @@
       include 'coupl.inc'
       include 'nFKSconfigs.inc'
       include "fks_info.inc"
+      include "mint.inc"
       LOGICAL  IS_A_J(NEXTERNAL),IS_A_LP(NEXTERNAL),IS_A_LM(NEXTERNAL)
       LOGICAL  IS_A_PH(NEXTERNAL)
       COMMON /TO_SPECISA/IS_A_J,IS_A_LP,IS_A_LM,IS_A_PH
@@ -151,13 +152,14 @@
       double precision pmass(-nexternal:0,lmaxconfigs)
       double precision pwidth(-nexternal:0,lmaxconfigs)
       integer pow(-nexternal:0,lmaxconfigs)
-      integer itree(2,-max_branch:-1),iconfig
-      common /to_itree/itree,iconfig
+      integer itree(2,-max_branch:-1),iconf
+      common /to_itree/itree,iconf
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
-      double precision taumin(fks_configs),taumin_s(fks_configs)
-     &     ,taumin_j(fks_configs),stot,xk(nexternal)
-      save  taumin,taumin_s,taumin_j,stot
+      double precision taumin(fks_configs,maxchannels)
+     $     ,taumin_s(fks_configs,maxchannels),taumin_j(fks_configs
+     $     ,maxchannels),stot,xk(nexternal)
+      save  taumin,taumin_s,taumin_j
       integer i,j,k,d1,d2,iFKS,nt
       double precision xm(-nexternal:nexternal),xm1,xm2,xmi
       double precision xw(-nexternal:nexternal),xw1,xw2,xwi
@@ -170,11 +172,12 @@
       double precision mass_min(-nexternal:nexternal),masslow(
      $     -nexternal:-1),widthlow(-nexternal:-1),sum_all_s
       integer t_channel
-      integer cBW_FKS_level_max(fks_configs),
-     &     cBW_FKS(fks_configs,-nexternal:-1),
-     &     cBW_FKS_level(fks_configs,-nexternal:-1)
-      double precision cBW_FKS_mass(fks_configs,-1:1,-nexternal:-1),
-     &     cBW_FKS_width(fks_configs,-1:1,-nexternal:-1)
+      integer cBW_FKS_level_max(fks_configs,maxchannels),
+     &     cBW_FKS(fks_configs,-nexternal:-1,maxchannels),
+     &     cBW_FKS_level(fks_configs,-nexternal:-1,maxchannels)
+      double precision cBW_FKS_mass(fks_configs,-1:1,-nexternal:-1
+     $     ,maxchannels),cBW_FKS_width(fks_configs,-1:1,-nexternal:-1
+     $     ,maxchannels)
       save cBW_FKS_level_max,cBW_FKS,cBW_FKS_level,cBW_FKS_mass
      $     ,cBW_FKS_width
       integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
@@ -183,7 +186,7 @@
       common/c_conflictingBW/cBW_mass,cBW_width,cBW_level_max,cBW
      $     ,cBW_level
       double precision s_mass(-nexternal:nexternal)
-     $     ,s_mass_FKS(fks_configs,-nexternal:nexternal)
+     $     ,s_mass_FKS(fks_configs,-nexternal:nexternal,maxchannels)
       save s_mass_FKS
       common/to_phase_space_s_channel/s_mass
 c Les Houches common block
@@ -192,18 +195,19 @@
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
       common /c_leshouche_inc/idup,mothup,icolup,niprocs
-c
       real*8         emass(nexternal)
       common/to_mass/emass
-      logical firsttime
+      logical firsttime,firsttime_chans(maxchannels)
       data firsttime /.true./
+      data firsttime_chans/maxchannels*.true./
       if (firsttime) then
          do i = 1,lmaxconfigs
             do j = -nexternal,0
                pmass(j,i) = 0d0
                pwidth(j,i) = 0d0
-            end do
-         end do
+            enddo
+         enddo
+         firsttime=.false.
       endif
       include "born_props.inc"
 
@@ -221,29 +225,29 @@
 c not necessarily the softest.  Therefore, it could be that even though
 c the Born does not have enough energy to pass the cuts set by ptj, the
 c event could.
-      if (firsttime) then
+      if (firsttime_chans(ichan)) then
          do i=-nexternal,nexternal
             xm(i)=0d0
             xw(i)=0d0
             mass_min(i)=0d0
          end do
-         firsttime=.false.
+         firsttime_chans(ichan)=.false.
          do iFKS=1,fks_configs
             j_fks=FKS_J_D(iFKS)
-            taumin(iFKS)=0.d0
-            taumin_s(iFKS)=0.d0
-            taumin_j(iFKS)=0.d0
+            taumin(iFKS,ichan)=0.d0
+            taumin_s(iFKS,ichan)=0.d0
+            taumin_j(iFKS,ichan)=0.d0
             do i=nincoming+1,nexternal
 c Add the minimal jet pTs to tau
                if(IS_A_J(i) .and. i.ne.nexternal)then
                   if  (j_fks.gt.nincoming .and. j_fks.lt.nexternal) then
-                     taumin(iFKS)=taumin(iFKS)+max(ptj,emass(i))
-                     taumin_s(iFKS)=taumin_s(iFKS)+max(ptj,emass(i))
-                     taumin_j(iFKS)=taumin_j(iFKS)+max(ptj,emass(i))
+                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+max(ptj,emass(i))
+                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+max(ptj,emass(i))
+                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+max(ptj,emass(i))
                   elseif (j_fks.ge.1 .and. j_fks.le.nincoming) then
-                     taumin(iFKS)=taumin(iFKS)+emass(i)
-                     taumin_s(iFKS)=taumin_s(iFKS)+max(ptj,emass(i))
-                     taumin_j(iFKS)=taumin_j(iFKS)+max(ptj,emass(i))
+                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
+                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+max(ptj,emass(i))
+                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+max(ptj,emass(i))
                   elseif (j_fks.eq.nexternal) then
                      write (*,*)
      &                    'ERROR, j_fks cannot be the final parton'
@@ -264,17 +268,17 @@
                      stop
                   endif
                   if  (j_fks.gt.nincoming)
-     &                 taumin(iFKS)=taumin(iFKS)+ptgmin
-                  taumin_s(iFKS)=taumin_s(iFKS)+ptgmin
-                  taumin_j(iFKS)=taumin_j(iFKS)+ptgmin
+     &                 taumin(iFKS,ichan)=taumin(iFKS,ichan)+ptgmin
+                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+ptgmin
+                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+ptgmin
                   xm(i)=emass(i)+ptgmin
                elseif (is_a_lp(i)) then
 c Add the postively charged lepton pTs to tau
-                  taumin(iFKS)=taumin(iFKS)+emass(i)
+                  taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
                   if (j_fks.gt.nincoming)
-     &                 taumin(iFKS)=taumin(iFKS)+ptl
-                  taumin_s(iFKS)=taumin_s(iFKS)+emass(i)+ptl
-                  taumin_j(iFKS)=taumin_j(iFKS)+emass(i)+ptl
+     &                 taumin(iFKS,ichan)=taumin(iFKS,ichan)+ptl
+                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+emass(i)+ptl
+                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+emass(i)+ptl
                   xm(i)=emass(i)+ptl
 c Add the lepton invariant mass to tau if there is at least another
 c lepton of opposite charge. (Only add half of it, i.e. 'the part
@@ -284,22 +288,22 @@
                      if (is_a_lm(j) .and. idup(i,1).eq.-idup(j,1) .and.
      $                    (mll_sf.ne.0d0 .or. mll.ne.0d0) ) then
                         if (j_fks.gt.nincoming)
-     &                       taumin(iFKS) = taumin(iFKS)-ptl-emass(i) +
+     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-ptl-emass(i) +
      &                              max(mll/2d0,mll_sf/2d0,ptl+emass(i))
-                        taumin_s(iFKS) = taumin_s(iFKS)-ptl-emass(i)
+                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
      $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
-                        taumin_j(iFKS) = taumin_j(iFKS)-ptl-emass(i)
+                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
      $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
                         xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,mll_sf/2d0
      $                       ,ptl+emass(i))
                         exit
                      elseif (is_a_lm(j) .and. mll.ne.0d0) then
                         if (j_fks.gt.nincoming)
-     &                       taumin(iFKS)= taumin(iFKS)-ptl-emass(i) +
+     &                       taumin(iFKS,ichan)= taumin(iFKS,ichan)-ptl-emass(i) +
      &                                     max(mll/2d0,ptl+emass(i))
-                        taumin_s(iFKS) = taumin_s(iFKS)-ptl-emass(i)
+                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
      $                       + max(mll/2d0,ptl+emass(i))
-                        taumin_j(iFKS) = taumin_j(iFKS)-ptl-emass(i)
+                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
      $                       + max(mll/2d0,ptl+emass(i))
                         xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,ptl
      $                       +emass(i))
@@ -308,11 +312,11 @@
                   enddo
                elseif (is_a_lm(i)) then
 c Add the negatively charged lepton pTs to tau
-                  taumin(iFKS)=taumin(iFKS)+emass(i)
+                  taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
                   if (j_fks.gt.nincoming)
-     &                 taumin(iFKS)=taumin(iFKS)+ptl
-                  taumin_s(iFKS)=taumin_s(iFKS)+emass(i)+ptl
-                  taumin_j(iFKS)=taumin_j(iFKS)+emass(i)+ptl
+     &                 taumin(iFKS,ichan)=taumin(iFKS,ichan)+ptl
+                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+emass(i)+ptl
+                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+emass(i)+ptl
                   xm(i)=emass(i)+ptl
 c Add the lepton invariant mass to tau if there is at least another
 c lepton of opposite charge. (Only add half of it, i.e. 'the part
@@ -322,22 +326,22 @@
                      if (is_a_lp(j) .and. idup(i,1).eq.-idup(j,1) .and.
      $                    (mll_sf.ne.0d0 .or. mll.ne.0d0) ) then
                         if (j_fks.gt.nincoming)
-     &                       taumin(iFKS) = taumin(iFKS)-ptl-emass(i) +
+     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-ptl-emass(i) +
      &                              max(mll/2d0,mll_sf/2d0,ptl+emass(i))
-                        taumin_s(iFKS) = taumin_s(iFKS)-ptl-emass(i)
+                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
      $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
-                        taumin_j(iFKS) = taumin_j(iFKS)-ptl-emass(i)
+                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
      $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
                         xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,mll_sf/2d0
      $                       ,ptl+emass(i))
                         exit
                      elseif (is_a_lp(j) .and. mll.ne.0d0) then
                         if (j_fks.gt.nincoming)
-     &                       taumin(iFKS) = taumin(iFKS)-ptl-emass(i) +
+     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-ptl-emass(i) +
      &                                      max(mll/2d0,ptl+emass(i))
-                        taumin_s(iFKS) = taumin_s(iFKS)-ptl-emass(i)
+                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
      $                       + max(mll/2d0,ptl+emass(i))
-                        taumin_j(iFKS) = taumin_j(iFKS)-ptl-emass(i)
+                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
      $                       + max(mll/2d0,ptl+emass(i))
                         xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,ptl
      $                       +emass(i))
@@ -345,16 +349,16 @@
                      endif
                   enddo
                else
-                  taumin(iFKS)=taumin(iFKS)+emass(i)
-                  taumin_s(iFKS)=taumin_s(iFKS)+emass(i)
-                  taumin_j(iFKS)=taumin_j(iFKS)+emass(i)
+                  taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
+                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+emass(i)
+                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+emass(i)
                   xm(i)=emass(i)
                endif
                xw(i)=0d0
             enddo
             stot = 4d0*ebeam(1)*ebeam(2)
-            tau_Born_lower_bound=taumin(iFKS)**2/stot
-            tau_lower_bound=taumin_j(iFKS)**2/stot
+            tau_Born_lower_bound=taumin(iFKS,ichan)**2/stot
+            tau_lower_bound=taumin_j(iFKS,ichan)**2/stot
 c         
 c Also find the minimum lower bound if all internal s-channel particles
 c were on-shell
@@ -373,9 +377,9 @@
                   xw1=xw(d1)
                   xw2=xw(d2)
 c On-shell mass of the intermediate resonance
-                  xmi=pmass(i,iconfig)
+                  xmi=pmass(i,iconf)
 c Width of the intermediate resonance
-                  xwi=pwidth(i,iconfig)
+                  xwi=pwidth(i,iconf)
 c Set the intermediate mass equal to the max of its actual mass and
 c the sum of the masses of the two daugters.
                   if (xmi.gt.xm1+xm2) then
@@ -388,7 +392,7 @@
 c Add the new mass to the bound. To avoid double counting, we should
 c subtract the daughters, because they are already included above or in
 c the previous iteration of the loop
-                  taumin_s(iFKS)=taumin_s(iFKS)+xm(i)-xm1-xm2
+                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+xm(i)-xm1-xm2
                else             ! t-channels
                   if (i.eq.-(nexternal-3)) then
                      xm(i)=0d0
@@ -409,10 +413,10 @@
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c Determine the "minimal" s-channel invariant masses
             do i=nincoming+1,nexternal-1
-               s_mass_FKS(iFKS,i)=xm(i)**2
+               s_mass_FKS(iFKS,i,ichan)=xm(i)**2
             enddo
             do i=-1,-(nexternal-3),-1 ! All propagators 
-               s_mass_FKS(iFKS,i)=xm(i)**2
+               s_mass_FKS(iFKS,i,ichan)=xm(i)**2
             enddo
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c Determine the conflicting Breit-Wigner's. Note that xm(i) contains the
@@ -420,13 +424,13 @@
             do i=nincoming+1,nexternal-1
                mass_min(i)=xm(i) ! minimal allowed resonance mass (including masses set by cuts)
             enddo
-            cBW_FKS_level_max(iFKS)=0
+            cBW_FKS_level_max(iFKS,ichan)=0
             t_channel=0
             do i=-1,-(nexternal-3),-1 ! All propagators
-               cBW_FKS_mass(iFKS,1,i)=0d0
-               cBW_FKS_width(iFKS,1,i)=0d0
-               cBW_FKS_mass(iFKS,-1,i)=0d0
-               cBW_FKS_width(iFKS,-1,i)=0d0
+               cBW_FKS_mass(iFKS,1,i,ichan)=0d0
+               cBW_FKS_width(iFKS,1,i,ichan)=0d0
+               cBW_FKS_mass(iFKS,-1,i,ichan)=0d0
+               cBW_FKS_width(iFKS,-1,i,ichan)=0d0
                masslow(i)=9d99
                widthlow(i)=0d0
                if ( itree(1,i).eq.1 .or. itree(1,i).eq.2 ) t_channel=i
@@ -438,39 +442,39 @@
      $                 ,xm(i),mass_min(i)
                   stop
                endif
-               if (pmass(i,iconfig).lt.xm(i) .and.
-     $              pwidth(i,iconfig).gt.0d0) then
+               if (pmass(i,iconf).lt.xm(i) .and.
+     $              pwidth(i,iconf).gt.0d0) then
 c     Possible conflict in BW
-                  if (pmass(i,iconfig).lt.mass_min(i)) then
+                  if (pmass(i,iconf).lt.mass_min(i)) then
 c     Resonance can never go on-shell due to the kinematics of the event
-                     cBW_FKS(iFKS,i)=2
-                     cBW_FKS_level(iFKS,i)=0
-                  elseif(pmass(i,iconfig).lt.xm(i)) then
+                     cBW_FKS(iFKS,i,ichan)=2
+                     cBW_FKS_level(iFKS,i,ichan)=0
+                  elseif(pmass(i,iconf).lt.xm(i)) then
 c     Conflicting Breit-Wigner
-                     cBW_FKS(iFKS,i)=1
-                     cBW_FKS_level(iFKS,i)=1
-                     cBW_FKS_level_max(iFKS)=max(cBW_FKS_level_max(iFKS)
-     $                    ,cBW_FKS_level(iFKS,i))
+                     cBW_FKS(iFKS,i,ichan)=1
+                     cBW_FKS_level(iFKS,i,ichan)=1
+                     cBW_FKS_level_max(iFKS,ichan)=max(cBW_FKS_level_max(iFKS,ichan)
+     $                    ,cBW_FKS_level(iFKS,i,ichan))
 c     Set here the mass (and width) of the alternative mass; it's the
 c     sum of daughter masses. (2nd argument is '1', because this
 c     alternative mass is LARGER than the resonance mass).
-                     cBW_FKS_mass(iFKS,1,i)=xm(i)
-                     cBW_FKS_width(iFKS,1,i)=xw(i)
+                     cBW_FKS_mass(iFKS,1,i,ichan)=xm(i)
+                     cBW_FKS_width(iFKS,1,i,ichan)=xw(i)
                   endif
 c     set the daughters also as conflicting (recursively)
-                  masslow(i)=pmass(i,iconfig)
-                  widthlow(i)=pwidth(i,iconfig)
+                  masslow(i)=pmass(i,iconf)
+                  widthlow(i)=pwidth(i,iconf)
                   do j=i,-1
-                     if (cBW_FKS(iFKS,j).eq.0) cycle
+                     if (cBW_FKS(iFKS,j,ichan).eq.0) cycle
                      do k=1,2   ! loop over the 2 daughters
                         if (itree(k,j).ge.0) cycle
-                        if (cBW_FKS(iFKS,itree(k,j)).eq.2) cycle
-                        cBW_FKS(iFKS,itree(k,j))=1
-                        cBW_FKS_level(iFKS,itree(k,j))=
-     $                       cBW_FKS_level(iFKS,j)+1
-                        cBW_FKS_level_max(iFKS)=
-     $                       max(cBW_FKS_level_max(iFKS)
-     $                       ,cBW_FKS_level(iFKS,itree(k,j)))
+                        if (cBW_FKS(iFKS,itree(k,j),ichan).eq.2) cycle
+                        cBW_FKS(iFKS,itree(k,j),ichan)=1
+                        cBW_FKS_level(iFKS,itree(k,j),ichan)=
+     $                       cBW_FKS_level(iFKS,j,ichan)+1
+                        cBW_FKS_level_max(iFKS,ichan)=
+     $                       max(cBW_FKS_level_max(iFKS,ichan)
+     $                       ,cBW_FKS_level(iFKS,itree(k,j),ichan))
 c     Set here the mass (and width) of the alternative mass; it's the
 c     difference between the mother and the sister masses. (3rd argument
 c     is '-1', because this alternative mass is SMALLER than the
@@ -479,18 +483,18 @@
      &                       max(masslow(j)-xm(itree(3-k,j)),0d0)) ! mass difference
                         widthlow(itree(k,j))=max(widthlow(itree(k,j)),
      &                       widthlow(j)+xw(itree(3-k,j))) ! sum of widths
-                        if (pwidth(itree(k,j),iconfig).eq.0d0 .or.
+                        if (pwidth(itree(k,j),iconf).eq.0d0 .or.
      $                       masslow(itree(k,j)).ge.pmass(itree(k,j)
-     $                       ,iconfig)) cycle
-                        cBW_FKS_mass(iFKS,-1,itree(k,j))=
+     $                       ,iconf)) cycle
+                        cBW_FKS_mass(iFKS,-1,itree(k,j),ichan)=
      $                       masslow(itree(k,j))
-                        cBW_FKS_width(iFKS,-1,itree(k,j))=
+                        cBW_FKS_width(iFKS,-1,itree(k,j),ichan)=
      $                       widthlow(itree(k,j))
                      enddo
                   enddo
                else
 c     Normal Breit-Wigner
-                  cBW_FKS(iFKS,i)=0
+                  cBW_FKS(iFKS,i,ichan)=0
                endif
             enddo
 c loop over t-channel to make sure that s-hat is consistent with sum of
@@ -500,9 +504,9 @@
                do i=t_channel,-(nexternal-3),-1
 c Breit-wigner can never go on-shell:
                   if (itree(2,i).gt.0) cycle
-                  if ( pmass(itree(2,i),iconfig).gt.sqrt(stot) .and.
-     $                 pwidth(itree(2,i),iconfig).gt.0d0) then
-                     cBW_FKS(iFKS,itree(2,i))=2
+                  if ( pmass(itree(2,i),iconf).gt.sqrt(stot) .and.
+     $                 pwidth(itree(2,i),iconf).gt.0d0) then
+                     cBW_FKS(iFKS,itree(2,i),ichan)=2
                   endif
 c     s-channel is always 2nd argument of itree, sum it to sum_all_s
                   sum_all_s=sum_all_s+xm(itree(2,i))
@@ -511,10 +515,10 @@
 c     conflicting BWs: set all s-channels as conflicting
                   do i=t_channel,-(nexternal-3),-1
                      if (itree(2,i).gt.0) cycle
-                     if (cBW_FKS(iFKS,itree(2,i)).ne.2) then
-                        cBW_FKS(iFKS,itree(2,i))=1
-                        cBW_FKS_mass(iFKS,-1,itree(2,i))=sqrt(stot)/2d0
-                        cBW_FKS_width(iFKS,-1,itree(2,i))=xw(itree(2,i))
+                     if (cBW_FKS(iFKS,itree(2,i),ichan).ne.2) then
+                        cBW_FKS(iFKS,itree(2,i),ichan)=1
+                        cBW_FKS_mass(iFKS,-1,itree(2,i),ichan)=sqrt(stot)/2d0
+                        cBW_FKS_width(iFKS,-1,itree(2,i),ichan)=xw(itree(2,i))
                      endif
                   enddo
                endif
@@ -526,41 +530,38 @@
 c
 c If the lower bound found here is smaller than the hard bound,
 c simply set the soft bound equal to the hard bound.
-            taumin_s(iFKS)=
-     &           max(taumin_j(iFKS),taumin_s(iFKS))
+            taumin_s(iFKS,ichan)=
+     &           max(taumin_j(iFKS,ichan),taumin_s(iFKS,ichan))
 c
 c For the bound, we have to square and divide by stot.
-            tau_lower_bound_resonance=taumin_s(iFKS)**2/stot
+            tau_lower_bound_resonance=taumin_s(iFKS,ichan)**2/stot
 c
-            write (*,'(a,i3,a,3(e12.5,x))') 'nFKSprocess:',iFKS
-     &           ,'. Absolute lower bound for tau at the Born is'
-     &           ,tau_Born_lower_bound,taumin(iFKS),dsqrt(stot) 
-            if (j_fks.le.nincoming) then
-               write (*,'(a,i3,a,3(e12.5,x))') 'nFKSprocess:',iFKS
-     &              ,'. Lower bound for tau is',tau_lower_bound
-     &              ,taumin_j(iFKS),dsqrt(stot)
+            if (j_fks.gt.nincoming) then
+               write (*,'(a7,x,i3,x,i5,x,a1,3(e12.5,x)))') 'tau_min'
+     $              ,iFKS,ichan,':',taumin(iFKS,ichan),taumin_j(iFKS
+     $              ,ichan),taumin_s(iFKS,ichan)
+            else
+               write (*,'(a7,x,i3,x,i5,x,a1,e12.5,x,a13,e12.5,x))')
+     $              'tau_min',iFKS,ichan,':',taumin(iFKS,ichan)
+     $              ,'     --      ',taumin_s(iFKS,ichan)
             endif
-            write (*,'(a,i3,a,3(e12.5,x))') 'nFKSprocess:',iFKS
-     &           ,'. Lower bound for tau is (taking resonances'/
-     &           /' into account)' ,tau_lower_bound_resonance
-     &           ,taumin_s(iFKS) ,dsqrt(stot)
          enddo
       endif
-      tau_Born_lower_bound=taumin(nFKSprocess)**2/stot
-      tau_lower_bound=taumin_j(nFKSprocess)**2/stot
-      tau_lower_bound_resonance=taumin_s(nFKSprocess)**2/stot
+      tau_Born_lower_bound=taumin(nFKSprocess,ichan)**2/stot
+      tau_lower_bound=taumin_j(nFKSprocess,ichan)**2/stot
+      tau_lower_bound_resonance=taumin_s(nFKSprocess,ichan)**2/stot
       do i=-nexternal,-1
-         cBW(i)=cBW_FKS(nFKSprocess,i)
-         cBW_level(i)=cBW_FKS_level(nFKSprocess,i)
+         cBW(i)=cBW_FKS(nFKSprocess,i,ichan)
+         cBW_level(i)=cBW_FKS_level(nFKSprocess,i,ichan)
          do j=-1,1,2
-            cBW_mass(j,i)=cBW_FKS_mass(nFKSprocess,j,i)
-            cBW_width(j,i)=cBW_FKS_width(nFKSprocess,j,i)
+            cBW_mass(j,i)=cBW_FKS_mass(nFKSprocess,j,i,ichan)
+            cBW_width(j,i)=cBW_FKS_width(nFKSprocess,j,i,ichan)
          enddo
       enddo
       do i=-nexternal,nexternal
-         s_mass(i)=s_mass_FKS(nFKSprocess,i)
+         s_mass(i)=s_mass_FKS(nFKSprocess,i,ichan)
       enddo
-      cBW_level_max=cBW_FKS_level_max(nFKSprocess)
+      cBW_level_max=cBW_FKS_level_max(nFKSprocess,ichan)
       return
       end
 
@@ -574,8 +575,8 @@
       double precision pmass(-nexternal:0,lmaxconfigs)
       double precision pwidth(-nexternal:0,lmaxconfigs)
       integer pow(-nexternal:0,lmaxconfigs)
-      integer itree(2,-max_branch:-1),iconfig
-      common /to_itree/itree,iconfig
+      integer itree(2,-max_branch:-1),iconf
+      common /to_itree/itree,iconf
       logical new_point
       common /c_new_point/new_point
       double precision ran2,rnd

=== modified file 'Template/NLO/SubProcesses/symmetry_fks_test_MC.f'
--- Template/NLO/SubProcesses/symmetry_fks_test_MC.f	2016-09-27 15:07:53 +0000
+++ Template/NLO/SubProcesses/symmetry_fks_test_MC.f	2017-03-08 10:27:27 +0000
@@ -14,6 +14,7 @@
       include 'fks_info.inc'
       include 'run.inc'
       include 'cuts.inc'
+      include 'mint.inc'
       
       double precision ZERO,one
       parameter       (ZERO = 0d0)
@@ -68,7 +69,7 @@
       double precision p(0:3,99), wgt, x(99), fx
       double complex wgt1(2)
       double precision p1(0:3,99),xx(maxinvar)
-      integer ninvar, ndim, iconfig, minconfig, maxconfig
+      integer ninvar, ndim,  minconfig, maxconfig
       common/tosigint/ndim
       integer ncall,itmax,nconfigs,ntry, ngraphs
       integer ic(nexternal,maxswitch), jc(12),nswitch
@@ -333,7 +334,9 @@
       endif
 
       do iconfig=bs_min,bs_max       ! Born configurations
-
+         ichan=1
+         iconfigs(1)=iconfig
+         
       call set_mc_matrices
 
       wgt=1d0

=== modified file 'Template/NLO/SubProcesses/symmetry_fks_test_ME.f'
--- Template/NLO/SubProcesses/symmetry_fks_test_ME.f	2016-09-27 15:07:53 +0000
+++ Template/NLO/SubProcesses/symmetry_fks_test_ME.f	2017-03-08 10:27:27 +0000
@@ -14,6 +14,7 @@
       include 'fks_info.inc'
       include 'run.inc'
       include 'cuts.inc'
+      include 'mint.inc'
       
       double precision ZERO,one
       parameter       (ZERO = 0d0)
@@ -59,7 +60,7 @@
       double precision p(0:3,99), wgt, x(99), fx
       double complex wgt1(2)
       double precision p1(0:3,99),xx(maxinvar)
-      integer ninvar, ndim, iconfig, minconfig, maxconfig
+      integer ninvar, ndim,  minconfig, maxconfig
       common/tosigint/ndim
       integer ncall,itmax,nconfigs,ntry, ngraphs
       integer ic(nexternal,maxswitch), jc(12),nswitch
@@ -279,7 +280,9 @@
          bs_max=iconfig_in
       endif
 
-      do iconfig=bs_min,bs_max       ! Born configurations
+      do iconfig=bs_min,bs_max  ! Born configurations
+         ichan=1
+         iconfigs(1)=iconfig
       call setcuts
       call setfksfactor(iconfig,.false.)
       wgt=1d0

=== modified file 'Template/NLO/SubProcesses/symmetry_fks_v3.f'
--- Template/NLO/SubProcesses/symmetry_fks_v3.f	2016-09-27 15:07:53 +0000
+++ Template/NLO/SubProcesses/symmetry_fks_v3.f	2017-03-08 10:27:27 +0000
@@ -10,6 +10,7 @@
       include 'genps.inc'      
       include 'nexternal.inc'
       include '../../Source/run_config.inc'
+      include 'mint.inc'
       
       double precision ZERO
       parameter       (ZERO = 0d0)
@@ -49,7 +50,7 @@
      &     -1),p_ev_red_save(0:3,nexternal-1)
       double precision p1_cnt_save(0:3,nexternal,-2:2),p_born_save(0:3
      &     ,nexternal-1),p_ev_red1(0:3,nexternal-1)
-      integer ninvar, ndim, iconfig, minconfig, maxconfig
+      integer ninvar, ndim, minconfig, maxconfig
       common/tosigint/ndim
       integer ncall,itmax,nconfigs,ntry, ngraphs
       integer icb(nexternal-1,maxswitch),jc(12),nswitch
@@ -180,6 +181,8 @@
       call printout
       call run_printout
       iconfig=1
+      ichan=1
+      iconfigs(1)=iconfig
       call setfksfactor(iconfig,.false.)
 c
       ndim = 55

=== modified file 'Template/NLO/SubProcesses/write_event.f'
--- Template/NLO/SubProcesses/write_event.f	2016-05-13 12:40:54 +0000
+++ Template/NLO/SubProcesses/write_event.f	2017-03-08 10:27:27 +0000
@@ -5,16 +5,14 @@
       include "unlops.inc"
       include "run.inc"
       include 'timing_variables.inc'
+      include 'mint.inc'
       integer ndim
       common/tosigint/ndim
-      integer           iconfig
-      common/to_configs/iconfig
       integer itmax,ncall
       common/citmax/itmax,ncall
       logical Hevents
       common/SHevents/Hevents
       integer i,j,lunlhe
-      include 'mint.inc'
       real*8 xx(ndimmax),weight,evnt_wgt
       logical putonshell
       double precision wgt,unwgtfun
@@ -263,7 +261,7 @@
             write(*,*)doreweight,iwgtinfo
             stop
           endif
-          kwgtinfo=iwgtinfo
+          kwgtinfo= iwgtinfo
           write(buff,201)'#aMCatNLO',iSorH_lhe,ifks_lhe(nFKSprocess)
      &         ,jfks_lhe(nFKSprocess),fksfather_lhe(nFKSprocess)
      &         ,ipartner_lhe(nFKSprocess),scale1_lhe(nFKSprocess)

=== modified file 'Template/loop_material/StandAlone/SubProcesses/makefile'
--- Template/loop_material/StandAlone/SubProcesses/makefile	2016-06-24 18:12:15 +0000
+++ Template/loop_material/StandAlone/SubProcesses/makefile	2017-03-08 10:27:27 +0000
@@ -81,9 +81,11 @@
 # This is the core of madloop computationally wise, so make sure to turn optimizations on and bound checks off.
 # We use %olynomial.o and not directly polynomial.o because we want it to match when both doing make check here
 # or make OLP one directory above
+%oloop_matrix.o : %olynomial.o %oloop_matrix.f
 %olynomial.o : %olynomial.f
 	$(FC) $(patsubst -O%,, $(subst -fbounds-check,,$(FFLAGS))) $(POLYNOMIAL_OPTIMIZATION) $(POLYNOMIAL_BOUNDS_CHECK) -c $< -o $@ $(LOOP_INCLUDE)
 
+
 $(DOTO) : $(DOTF)
 	$(FC) $(FFLAGS) -c $< -o $@ $(LOOP_INCLUDE)
 

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2016-12-10 02:08:09 +0000
+++ UpdateNotes.txt	2017-03-08 10:27:27 +0000
@@ -1,5 +1,25 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
+2.5.3(XX/XX/17)
+      PT: Modified the default shower starting scale in montecarlocounter.f.
+          The new reference scale from which the dampening profile is computed is sum_i mt_i/2, i
+          being Born level final-state momenta.
+      OM: New "special mode" for madspin accessible via "set spinmode onshell".
+          This mode allow for full spin-correlation with 3 (or more) body decay but the decaying particle
+          remains exactly onshell in this mode. Loop-induced production/decay are not allowed in this mode.
+      OM+RF: Allowing for creation of LHE like output of the fixed order run at NLO.
+          This LHEF file is unvalid for parton-shower (all PS should crash on such file). It will be 
+          unphysical to shower such sample anyway.
+          Two hidden parameters of the  FO_analyse_card.dat allow some control on the LHEF creation
+       	  "fo_lhe_weight_ratio" allows to control the strength of a partial unweighting [default:1e-3]
+             increasing this number reduce the LHEF size.
+          "fo_lhe_postprocessing" can take value like nogrouping, norandom, noidentification.
+             nogrouping forbids the appearance of the LHEF(version2) tag <eventgroup>
+             norandom   does not apply the randomization of the events.
+             noidentification does not merge born event with other born like counter-event
+      RF: Better job handling for fNLO runs.
+      VH: Fixing various problem with the pythia8 interface (especially for MLM merging)
+
 2.5.2(10/12/16)
       OM: improve systematics (thanks to Philipp Pigard)
       OM: new syntax to modify the run_card: set no_parton_cut

=== modified file 'aloha/aloha_lib.py'
--- aloha/aloha_lib.py	2015-10-22 22:04:45 +0000
+++ aloha/aloha_lib.py	2017-03-08 10:27:27 +0000
@@ -50,6 +50,11 @@
 import re
 import aloha # define mode of writting
 
+try:
+    import madgraph.various.misc as misc
+except Exception:
+    import aloha.misc as misc
+
 class defaultdict(collections.defaultdict):
 
     def __call__(self, *args):
@@ -158,7 +163,6 @@
         # check if the function is a pure numerical function.
         if (fct_tag.startswith('cmath.') or fct_tag in self.known_fct) and \
                     all(isinstance(x, numbers.Number) for x in argument):
-            print "start the hack"
             import cmath
             if fct_tag.startswith('cmath.'):
                 module = ''
@@ -1056,7 +1060,8 @@
         self.spin_ind = spin_indices #spin indices
         self.nb_spin = len(spin_indices) #their number
         self.nb_ind = self.nb_lor + self.nb_spin #total number of indices
-                
+        
+        
         #store the representation
         if self.lorentz_ind or self.spin_ind:
             dict.__init__(self, representation) 

=== modified file 'aloha/aloha_parsers.py'
--- aloha/aloha_parsers.py	2015-10-01 16:00:08 +0000
+++ aloha/aloha_parsers.py	2017-03-08 10:27:27 +0000
@@ -36,6 +36,11 @@
 from aloha.aloha_lib import KERNEL
 logger = logging.getLogger('aloha.parsers')
 
+try:
+    import madgraph.various.misc as misc
+except Exception:
+    import aloha.misc as misc
+
 
 # PLY lexer class
 class UFOExpressionParser(object):
@@ -268,7 +273,6 @@
         if re_groups:
             p1 = re_groups.group("name")
         new = aloha_lib.KERNEL.add_function_expression(p1, eval(p[3]))
-        
         p[0] = str(new)
     
     def p_expression_function2(self, p):

=== modified file 'aloha/create_aloha.py'
--- aloha/create_aloha.py	2016-04-15 11:14:25 +0000
+++ aloha/create_aloha.py	2017-03-08 10:27:27 +0000
@@ -211,9 +211,10 @@
         new_id = _conjugate_gap + old_id
         
         self.kernel_tag = set()
-        if not self.routine_kernel or isinstance(self.routine_kernel, str):            
+        aloha_lib.KERNEL.use_tag = set()
+        if not self.routine_kernel or isinstance(self.routine_kernel, str):    
             self.routine_kernel = eval(self.parse_expression(self.lorentz_expr))
-        
+            self.kernel_tag = aloha_lib.KERNEL.use_tag
         # We need to compute C Gamma^T C^-1 = C_ab G_cb (-1) C_cd 
         #                  = C_ac G_bc (-1) C_bd = C_ac G_bc C_db
         self.routine_kernel = \
@@ -399,7 +400,7 @@
         
         if factorize:
             lorentz = lorentz.factorize()
-            
+        
         lorentz.tag = set(aloha_lib.KERNEL.use_tag)
         return lorentz     
 

=== modified file 'bin/compile.py'
--- bin/compile.py	2015-06-24 12:55:44 +0000
+++ bin/compile.py	2017-03-08 10:27:27 +0000
@@ -17,6 +17,8 @@
 import sys
 import logging
 import time
+import shutil
+import subprocess
 # Get the parent directory (mg root) of the script real path (bin)
 # and add it to the current PYTHONPATH
 root_path = os.path.split(os.path.dirname(os.path.realpath( __file__ )))[0]
@@ -37,21 +39,44 @@
 pjoin = os.path.join
 
 class Compile_MG5:
-    
-    def __init__(self):
+
+    autorun = True
+    overwrite_configuration=True
+
+    def __init__(self, ext_programs):
         """ launch all the compilation """
-        
-        self.make_UFO_pkl()
-        self.make_v4_pkl()
-        self.make_stdHep()
-        self.make_CutTools()
-        self.make_IREGI()
-
-        #important for UCL cluster
-        files.cp(pjoin(MG5DIR,'input','.mg5_configuration_default.txt'),
-                 pjoin(MG5DIR,'input','mg5_configuration.txt'))
-        self.cmd = interface.MasterCmd()        
-        self.install_package()
+
+        # important to uclclus
+        if self.overwrite_configuration:
+            files.cp(pjoin(MG5DIR,'input','.mg5_configuration_default.txt'),
+                     pjoin(MG5DIR,'input','mg5_configuration.txt'))
+            
+        self.cmd = interface.MasterCmd()                
+
+        if self.autorun:
+            self.make_UFO_pkl()
+            self.make_v4_pkl()
+            self.make_stdHep()
+            self.make_CutTools()
+            self.make_IREGI()
+            self.install_package(ext_programs)
+            self.test_output_LO()
+            self.test_output_NLO()
+            self.precompilation(debug=True)
+            self.precompilation(debug=False)
+
+
+    def test_output_LO(self):
+        """do the output of a simple LO process to ensure that LO is correctly configure."""
+        self.cmd.exec_cmd('generate p p > t t~')
+        self.cmd.exec_cmd('output %s/TESTLO' %root_path)
+        shutil.rmtree('%s/TESTLO' % root_path)
+
+    def test_output_NLO(self):
+        """do the output of a simple LO process to ensure that LO is correctly configure."""
+        self.cmd.exec_cmd('generate p p > e+ ve [QCD]')
+        self.cmd.exec_cmd('output %s/TESTNLO' %root_path)
+        shutil.rmtree('%s/TESTNLO' % root_path)
 
     @staticmethod
     def make_v4_pkl():
@@ -250,13 +275,22 @@
 
         misc.compile(cwd = os.path.join(iregi_path,'src'))        
 
-    def install_package(self):
+    def install_package(self, programs=[]):
         print "installing external package"
-        self.cmd.exec_cmd('install pythia-pgs')
-        self.cmd.exec_cmd('install Delphes')
-        self.cmd.exec_cmd('install ExRootAnalysis')
-        self.cmd.exec_cmd('install MadAnalysis')
-        self.cmd.exec_cmd('install SysCalc')
+        if not programs:
+            programs = ['pythia-pgs','Delphes','ExRootAnalysis','MadAnalysis4','SysCalc']
+            
+        for prog in programs:
+            self.cmd.exec_cmd('install %s' % prog)
+
+    def precompilation(self, debug=False):
+        if debug:
+            subprocess.call('python -m compileall .', shell=True, cwd=root_path)            
+        else:
+            subprocess.call('python -O -m compileall .', shell=True, cwd=root_path)
 
 if __name__ == '__main__':
-    Compile_MG5()
+    Compile_MG5(sys.argv[1:])
+
+
+

=== modified file 'bin/create_release.py'
--- bin/create_release.py	2016-09-07 17:11:30 +0000
+++ bin/create_release.py	2017-03-08 10:27:27 +0000
@@ -176,7 +176,11 @@
 shutil.rmtree(path.join(filepath, '.bzr'))
 for data in glob.glob(path.join(filepath, 'bin', '*')):
     if not data.endswith('mg5') and not data.endswith('mg5_aMC'):
-        os.remove(data)
+        if 'compile.py' not in data:
+            os.remove(data)
+        else:
+            os.rename(data, data.replace('compile.py','.compile.py'))
+
 os.remove(path.join(filepath, 'README.developer'))
 shutil.move(path.join(filepath, 'README.release'), path.join(filepath, 'README'))
 

=== modified file 'bin/mg5_aMC'
--- bin/mg5_aMC	2016-05-31 16:43:42 +0000
+++ bin/mg5_aMC	2017-03-08 10:27:27 +0000
@@ -1,4 +1,6 @@
 #! /usr/bin/env python
+import time
+start = time.time()
 
 ################################################################################
 #
@@ -103,7 +105,6 @@
 if __debug__:
         print 'Running MG5 in debug mode'
 
-
 # Set logging level according to the logging level given by options
 #logging.basicConfig(level=vars(logging)[options.logging])
 try:
@@ -120,8 +121,12 @@
     logging.getLogger('madevent').setLevel(level)
 except:
     pass
+
 import madgraph.interface.master_interface as interface
 
+if __debug__ and time.time() - start > 0.5:
+    print 'WARNING: loading of madgraph too slow!!!', time.time() - start
+
 if options.plugin:
     if not os.path.exists(os.path.join(root_path, 'PLUGIN', options.plugin)):
         print "ERROR: %s is not present in the PLUGIN directory. Please install it first"
@@ -140,6 +145,8 @@
 else: 
     cmd_line = interface.MasterCmd(mgme_dir = options.mgme_dir)
 
+
+
 # Call the cmd interface main loop
 try:
     if options.file or args:

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2016-12-06 11:25:57 +0000
+++ madgraph/core/base_objects.py	2017-03-08 10:27:27 +0000
@@ -1223,6 +1223,10 @@
             self['order_hierarchy'] = {}
             self['expansion_order'] = None
 
+        if name == 'name2pdg':
+            self['name2pgg'] = value
+            return
+
         result = Model.__bases__[0].set(self, name, value, force) # call the mother routine
 
         if name == 'particles':

=== modified file 'madgraph/core/diagram_generation.py'
--- madgraph/core/diagram_generation.py	2016-08-07 07:16:07 +0000
+++ madgraph/core/diagram_generation.py	2017-03-08 10:27:27 +0000
@@ -1333,7 +1333,7 @@
         self['decay_chains'] = DecayChainAmplitudeList()
 
     def __init__(self, argument = None, collect_mirror_procs = False,
-                 ignore_six_quark_processes = False, loop_filter=None):
+                 ignore_six_quark_processes = False, loop_filter=None, diagram_filter=False):
         """Allow initialization with Process and with ProcessDefinition"""
 
         if isinstance(argument, base_objects.Process):
@@ -1348,11 +1348,13 @@
                   MultiProcessClass.generate_multi_amplitudes(argument,
                                                     collect_mirror_procs,
                                                     ignore_six_quark_processes,
-                                                    loop_filter=loop_filter))
+                                                    loop_filter=loop_filter,
+                                                    diagram_filter=diagram_filter))
             else:
                 self['amplitudes'].append(\
                   MultiProcessClass.get_amplitude_from_proc(argument,
-                                                       loop_filter=loop_filter))
+                                                       loop_filter=loop_filter,
+                                                       diagram_filter=diagram_filter))
                 # Clean decay chains from process, since we haven't
                 # combined processes with decay chains yet
                 process = copy.copy(self.get('amplitudes')[0].get('process'))
@@ -1630,7 +1632,8 @@
                     self['amplitudes'].append(\
                         DecayChainAmplitude(process_def,
                                        self.get('collect_mirror_procs'),
-                                       self.get('ignore_six_quark_processes')))
+                                       self.get('ignore_six_quark_processes'),
+                                       diagram_filter=self['diagram_filter']))
                 else:
                     self['amplitudes'].extend(\
                        self.generate_multi_amplitudes(process_def,
@@ -1666,7 +1669,8 @@
 
         # Set automatic coupling orders
         process_definition.set('orders', MultiProcess.\
-                               find_optimal_process_orders(process_definition))
+                               find_optimal_process_orders(process_definition,
+                               diagram_filter))
         # Check for maximum orders from the model
         process_definition.check_expansion_orders()
 
@@ -1848,7 +1852,7 @@
         
 
     @staticmethod
-    def find_optimal_process_orders(process_definition):
+    def find_optimal_process_orders(process_definition, diagram_filter=False):
         """Find the minimal WEIGHTED order for this set of processes.
 
         The algorithm:
@@ -2022,7 +2026,7 @@
 
                     amplitude = Amplitude({'process': process})
                     try:
-                        amplitude.generate_diagrams()
+                        amplitude.generate_diagrams(diagram_filter=diagram_filter)
                     except InvalidCmd:
                         failed_procs.append(tuple(sorted_legs))
                     else:

=== modified file 'madgraph/core/helas_objects.py'
--- madgraph/core/helas_objects.py	2016-03-22 20:44:25 +0000
+++ madgraph/core/helas_objects.py	2017-03-08 10:27:27 +0000
@@ -1379,6 +1379,18 @@
         # wavefunction number
         return new_wf, wf_number
 
+    def has_multifermion(self):
+        """check the presence of 4 fermion vertex"""
+        
+        mothers = self.get('mothers')
+        if len(mothers) >2:
+            nb_fermion = len([1 for wf in mothers if wf.is_fermion()])
+            if nb_fermion>2:
+                return True
+        
+        return any(wf.has_multifermion() for wf in self.get('mothers'))
+        
+
     def get_fermion_order(self):
         """Recursive function to get a list of fermion numbers
         corresponding to the order of fermions along fermion lines
@@ -1399,7 +1411,6 @@
 
         other_fermions = [wf for wf in self.get('mothers') if \
                           wf.is_fermion() and wf != fermion_mother]
-
         # Pick out bosons
         bosons = filter(lambda wf: wf.is_boson(), self.get('mothers'))
 
@@ -1426,7 +1437,7 @@
 
         if self.is_fermion():
             return [mother_list[0], fermion_number_list]
-
+        
         return fermion_number_list
 
     def needs_hermitian_conjugate(self):
@@ -2618,7 +2629,43 @@
         mystr = mystr + '\n}'
 
         return mystr
-
+    
+    def has_multifermion(self):
+
+        return any(wf.has_multifermion() for wf in self.get('mothers'))
+        
+
+    def nice_string(self):
+        """ simple way to check which FD is related to this amplitude"""
+        def get_structure(wf):
+            """ funtion to allow to loop over helaswavefunction"""
+            mothers = []
+            try:
+                mothers = wf.get('mothers')
+            except:
+                if wf['is_loop']:
+                    return '%s*' % wf['particle'].get('pdg_code')
+                else:
+                    return  wf['particle'].get('pdg_code')
+                
+            struct = [get_structure(w) for w in mothers]
+            if struct:
+                if 'is_loop' in wf:
+                    if  wf['is_loop']:
+                        return (struct,'>%s*'%wf.get('pdg_code') )
+                    else: