Merge lp:~maddevelopers/mg5amcnlo/LI_NLO into lp:~olivier-mattelaer/mg5amcnlo/2.5.3

[Olivier Mattelaer]

Those are deep change in order to improve re-weighting (memory consumption and speed are highly improved due to better structure of the code (i.e. only use a single library))

[lp:~maddevelopers/mg5amcnlo/LI_NLO]

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

small fix of bug in reweighting

345. By olivier-mattelaer on 2017-03-29

merge with some old version of 2.5.4

346. By olivier-mattelaer on 2017-03-30

allow symmetric matrix element to go trough for the re-weighting

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

small fix

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

new fix (not very save but ok)

349. By Valentin Hirschi on 2017-04-21

1. Fixed an issue in the algorithm for numerically solving the rescaling factor for the PS point accuracy improvement using the parton-shower technique. This issue was affecting only the case where all external particles are massless.
   I solved it by rescaling the first momentum by x**2 instead of just x, so that the derivative of the function to optimize is not
   also zero when the optimzed function is zero.

350. By Valentin Hirschi on 2017-04-25

1. Added NaN's detection when computing the accuracy of the evaluation.
   If a NaN is detected the resulting accuracy is hardset to 1.0d99.

351. By Valentin Hirschi on 2017-05-25

1. Fixed the scale setting in the NLO all_matrix.f wrapper.

352. By Manoj Kumar Mandal on 2017-05-31

improve onshell particle when merging particle

353. By olivier-mattelaer on 2017-07-27

use the born configuration related to the weight rather than merging fks particles in presence of too soft/colinear sector

354. By olivier-mattelaer on 2017-07-27

remove one file

[Valentin Hirschi (valentin-hirschi)]

a) For:

- raise MadGraph5Error, 'user filter required to be defined in PLUGIN/user_filter.py with the function remove_diag(ONEDIAG) which returns True if the daigram has to be removed'
189 + try:
190 + from MG5aMC_PLUGIN.user_filter import remove_diag
191 + except ImportError:
192 + raise MadGraph5Error, 'user filter required to be defined in PLUGIN/user_filter.py with the function remove_diag(ONEDIAG) which returns True if the daigram has to be removed'

It would be nice to wrap the look-up into both PLUGIN and MG5aMC_PLUGIN in a misc.function so that we won't have to code up this dual plugin look-up everywhere.

b) What is the difference between

elif args[0] in ['virtual_path', 'tree_path']:

and change rwgt_me_dir

c) Remove
os.environ['GFORTRAN_UNBUFFERED_ALL'] =

d) Use a similar wrapper like in handle_param for the lines (460 rwgt interface):

 if self.rwgt_dir:
442 + path_me =self.rwgt_dir
443 + else:
444 + path_me = self.me_dir
445 +
446 + if self.second_model or self.second_process or self.dedicated_path:
447 + rw_dir = pjoin(path_me, 'rw_me_second')
448 + else:
449 + rw_dir = pjoin(path_me, 'rw_me')

e) Make sure to use nhel=-1 (and not 0!) for the convention 'sum over helicity' for MadLoop in the f2py splitter.

f) IMPORTANT: Make sure WIKI is updated with *ALL* the new fancy options of the reweight interface!! :)

review: Approve

[Olivier Mattelaer (olivier-mattelaer)]

> a)
> It would be nice to wrap the look-up into both PLUGIN and MG5aMC_PLUGIN in a
> misc.function so that we won't have to code up this dual plugin look-up
> everywhere.

replaced by:
            remove_diag = misc.plugin_import('user_filter',
                                             'user filter required to be defined in PLUGIN/user_filter.py with the function remove_diag(ONEDIAG) which returns True if the diagram has to be removed',
                                             fcts=['remove_diag'])

> b) What is the difference between
>
> elif args[0] in ['virtual_path', 'tree_path']:
>
> and change rwgt_me_dir

If the path to "rwgt_dir" does not exists,then this path will be created and all the code required for the code will be created at that position (and kept at the end of the run)
if the path exists, then those file will be re-used. (This include the standalone library for the original hyppothesis)

While 'virtual_path', 'tree_path' defines only path for the new hyppothesis if someone needs to specify them in advance (in principle only if someone needs to modify them like in the case of the two loop re-weighting).

> c) Remove
> os.environ['GFORTRAN_UNBUFFERED_ALL'] =

done

>
> d) Use a similar wrapper like in handle_param for the lines (460 rwgt
> interface):

maybe later.

>
> e) Make sure to use nhel=-1 (and not 0!) for the convention 'sum over
> helicity' for MadLoop in the f2py splitter.

For this I'm a bit confused. In the file f2py_wrapper, I see this:

      IF (NHEL.EQ.0) THEN
        CALL ML5_0_SLOOPMATRIX_THRES(P,MATELEM,-1.0D0, PREC_FOUND,
     $ RETURNCODE)
      ELSE
        CALL ML5_0_SLOOPMATRIXHEL_THRES(P,NHEL, MATELEM,-1.0D0,
     $ PREC_FOUND, RETURNCODE)
      ENDIF

So this sounds wrong according to your comment (at the same time, I did not modify such lines).
Now I'm a bit confused by the handling of nhel in madloop.
USERHEL is never reset to -1 inside the code at the end of the ML5_0_SLOOPMATRIXHEL
   [Like it is at tree level]
Do you want me to change
1) the above line to (NHEL.EQ.-1)
2) systematically reset USERHEL to -1 at the end of ML5_0_SLOOPMATRIXHEL

>
> f) IMPORTANT: Make sure WIKI is updated with *ALL* the new fancy options of
> the reweight interface!! :)

DONE

[lp:~maddevelopers/mg5amcnlo/LI_NLO]

355. By olivier-mattelaer on 2017-07-29

implement Valentin comment but the one on nhel

[Valentin Hirschi (valentin-hirschi)]

> > a)
> > It would be nice to wrap the look-up into both PLUGIN and MG5aMC_PLUGIN in a
> > misc.function so that we won't have to code up this dual plugin look-up
> > everywhere.
>
> replaced by:
> remove_diag = misc.plugin_import('user_filter',
> 'user filter required to be
> defined in PLUGIN/user_filter.py with the function remove_diag(ONEDIAG) which
> returns True if the diagram has to be removed',
> fcts=['remove_diag'])
>
> > b) What is the difference between
> >
> > elif args[0] in ['virtual_path', 'tree_path']:
> >
> > and change rwgt_me_dir
>
> If the path to "rwgt_dir" does not exists,then this path will be created and
> all the code required for the code will be created at that position (and kept
> at the end of the run)
> if the path exists, then those file will be re-used. (This include the
> standalone library for the original hyppothesis)
>
> While 'virtual_path', 'tree_path' defines only path for the new hyppothesis if
> someone needs to specify them in advance (in principle only if someone needs
> to modify them like in the case of the two loop re-weighting).

Got it. Can you confirm that the directories set with 'change virtual_path' and 'change tree_path' have indeed precedence over 'change rwgt_dir' if it is set as well?

>
> > c) Remove
> > os.environ['GFORTRAN_UNBUFFERED_ALL'] =
>
> done
>
> >
> > d) Use a similar wrapper like in handle_param for the lines (460 rwgt
> > interface):
>
> maybe later.
>
> >
> > e) Make sure to use nhel=-1 (and not 0!) for the convention 'sum over
> > helicity' for MadLoop in the f2py splitter.
>
> For this I'm a bit confused. In the file f2py_wrapper, I see this:
>
> IF (NHEL.EQ.0) THEN
> CALL ML5_0_SLOOPMATRIX_THRES(P,MATELEM,-1.0D0, PREC_FOUND,
> $ RETURNCODE)
> ELSE
> CALL ML5_0_SLOOPMATRIXHEL_THRES(P,NHEL, MATELEM,-1.0D0,
> $ PREC_FOUND, RETURNCODE)
> ENDIF
>
> So this sounds wrong according to your comment (at the same time, I did not
> modify such lines).
> Now I'm a bit confused by the handling of nhel in madloop.
> USERHEL is never reset to -1 inside the code at the end of the
> ML5_0_SLOOPMATRIXHEL
> [Like it is at tree level]
> Do you want me to change
> 1) the above line to (NHEL.EQ.-1)
> 2) systematically reset USERHEL to -1 at the end of ML5_0_SLOOPMATRIXHEL

Exactly, I had already noticed this "bug" of me not resetting USERHEL to -1 at the end of ML5_0_SLOOPMATRIXHEL and I have changed it in the MadEvent7 branch already.
So sorry for this; What you suggested is exactly what I would see appropriate indeed. Otherwise it is a bit weird to have two different conventions for the summation (i.e. NHEL = 0 or NHEL = -1) at two different levels of the code.

>
>
> >
> > f) IMPORTANT: Make sure WIKI is updated with *ALL* the new fancy options of
> > the reweight interface!! :)
>
> DONE

Seems good, I didn't re-read it carefully but I trust you that all what new functionalities that need to be documented now indeed are :)

=== modified file 'Template/loop_material/StandAlone/SubProcesses/makefile'
--- Template/loop_material/StandAlone/SubProcesses/makefile	2017-02-27 15:22:19 +0000
+++ Template/loop_material/StandAlone/SubProcesses/makefile	2017-07-29 19:33:46 +0000
@@ -68,6 +68,11 @@
  $(patsubst $(DOTF),$(DOTO),$(wildcard $(LOOP_PREFIX)*/COLLIER_interface.f)) \
  $(patsubst $(DOTF),$(DOTO),$(wildcard $(LOOP_PREFIX)*/compute_color_flows.f))
 
+POLYNOMIAL = $(patsubst $(DOTF),$(DOTO),$(wildcard polynomial.f))
+OLP_POLYNOMIAL = $(patsubst $(DOTF),$(DOTO),$(wildcard $(LOOP_PREFIX)*/polynomial.f))
+
+
+
 $(PROG): check_sa.o $(PROCESS) makefile $(LIBS)
 	$(FC) $(FFLAGS) -o $(PROG) check_sa.o $(PROCESS) $(LINKLIBS)
 
@@ -86,7 +91,7 @@
 	$(FC) $(patsubst -O%,, $(subst -fbounds-check,,$(FFLAGS))) $(POLYNOMIAL_OPTIMIZATION) $(POLYNOMIAL_BOUNDS_CHECK) -c $< -o $@ $(LOOP_INCLUDE)
 
 
-$(DOTO) : $(DOTF)
+$(DOTO) : $(DOTF) $(POLYNOMIAL) $(OLP_POLYNOMIAL)
 	$(FC) $(FFLAGS) -c $< -o $@ $(LOOP_INCLUDE)
 
 $(OLP): $(OLP_PROCESS) $(LIBS)
@@ -113,5 +118,9 @@
 	touch __init__.py
 	$(F2PY) $(MADLOOP_LIB) -m matrix$(MENUM)py -c f2py_wrapper.f --fcompiler=$(FC)   -L../../lib/ -ldhelas -lmodel $(LINK_LOOP_LIBS) $(STDLIB)
 
+allmatrix$(MENUM)py.so: $(OLP)_static all_matrix.f $(LIBS) $(WRAPPER)
+	touch __init__.py
+	$(F2PY) $(MADLOOP_LIB)  -m allmatrix$(MENUM)py -c all_matrix.f $(wildcard $(LOOP_PREFIX)*/f2py_wrapper.f) --fcompiler=$(FC)   -L../lib/ -ldhelas -lmodel $(LINK_LOOP_LIBS) $(STDLIB)
+
 clean:
 	@rm -f *.o *.so *.$(libext) *.$(dylibext)

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2017-01-11 21:56:50 +0000
+++ madgraph/core/base_objects.py	2017-07-29 19:33:46 +0000
@@ -2950,9 +2950,11 @@
                 if order not in self['sqorders_types']:
                     # Then assign its type to the default '='
                     self['sqorders_types'][order]='='
-
+                    
         return super(Process, self).get(name) # call the mother routine
 
+    
+
     def get_sorted_keys(self):
         """Return process property names as a nicely sorted list."""
 
@@ -3436,6 +3438,18 @@
         self['legs_with_decays'] = LegList(legs)
 
         return self['legs_with_decays']
+    
+    def get_tag(self):
+        """return the tag for standalone call"""
+        
+        initial = []    #filled in the next line
+        final = [l.get('id') for l in self.get('legs')\
+              if l.get('state') or initial.append(l.get('id'))]
+        decay_finals = self.get_final_ids_after_decay()
+        decay_finals.sort()
+        tag = (tuple(initial), tuple(decay_finals))
+        return tag
+        
 
     def list_for_sort(self):
         """Output a list that can be compared to other processes as:

=== modified file 'madgraph/core/diagram_generation.py'
--- madgraph/core/diagram_generation.py	2017-02-02 14:05:47 +0000
+++ madgraph/core/diagram_generation.py	2017-07-29 19:33:46 +0000
@@ -892,10 +892,9 @@
         list in argument."""
 
         if True:
-            try:
-                from PLUGIN.user_filter import remove_diag
-            except ImportError:
-                raise MadGraph5Error, 'user filter required to be defined in PLUGIN/user_filter.py with the function remove_diag(ONEDIAG) which returns True if the daigram has to be removed'
+            remove_diag = misc.plugin_import('user_filter', 
+                                             'user filter required to be defined in PLUGIN/user_filter.py with the function remove_diag(ONEDIAG) which returns True if the diagram has to be removed',
+                                             fcts=['remove_diag'])
         else:
             #example and simple tests
             def remove_diag(diag):

=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
--- madgraph/interface/amcatnlo_run_interface.py	2017-03-07 15:43:44 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2017-07-29 19:33:46 +0000
@@ -4784,8 +4784,6 @@
             file.write(content)
         file.close()
 
-
-
     ############################################################################
     def find_model_name(self):
         """ return the model name """

=== modified file 'madgraph/interface/common_run_interface.py'
--- madgraph/interface/common_run_interface.py	2017-03-03 23:41:19 +0000
+++ madgraph/interface/common_run_interface.py	2017-07-29 19:33:46 +0000
@@ -1878,7 +1878,7 @@
                 if 'nevt_job' in self.run_card and self.run_card['nevt_job'] !=-1:
                     nevt_job = self.run_card['nevt_job']
                 else:
-                    nevt_job = max(5000, self.run_card['nevents']/50)
+                    nevt_job = max(2500, self.run_card['nevents']/self.options['nb_core'])
                 logger.info("split the event file in bunch of %s events" % nevt_job)
                 nb_file = lhe_parser.EventFile(new_args[0]).split(nevt_job)
                 starttime = time.time()
@@ -5094,7 +5094,7 @@
             isinstance(self.run_card,banner_mod.RunCardNLO) and \
             not self.run_card['store_rwgt_info']:
             #check if a NLO reweighting is required
-                re_pattern = re.compile(r'''^\s*change\s*mode\s* (LO\+NLO|LO|NLO)\s*(?:#|$)''', re.M+re.I)
+                re_pattern = re.compile(r'''^\s*change\s*mode\s* (LO\+NLO|LO|NLO|NLO_tree)\s*(?:#|$)''', re.M+re.I)
                 text = open(self.paths['reweight']).read()
                 options = re_pattern.findall(text)
                 if any(o in ['NLO', 'LO+NLO'] for o in options):

=== modified file 'madgraph/interface/extended_cmd.py'
--- madgraph/interface/extended_cmd.py	2017-02-02 14:05:47 +0000
+++ madgraph/interface/extended_cmd.py	2017-07-29 19:33:46 +0000
@@ -36,7 +36,7 @@
 
 try:
     import madgraph.various.misc as misc
-    from madgraph import MG5DIR
+    from madgraph import MG5DIR, MadGraph5Error
     MADEVENT = False
 except ImportError, error:
     try:
@@ -862,7 +862,6 @@
     keyboard_stop_msg = """stopping all current operation
             in order to quit the program please enter exit"""
 
-
     if MADEVENT:
         plugin_path = []
     else:
@@ -1135,7 +1134,8 @@
         
         # Execute the card
         self.import_command_file(args[1])
-         
+    
+        
     def check_import(self, args):
         """check import command"""
         

=== modified file 'madgraph/interface/loop_interface.py'
--- madgraph/interface/loop_interface.py	2016-09-01 10:26:58 +0000
+++ madgraph/interface/loop_interface.py	2017-07-29 19:33:46 +0000
@@ -428,6 +428,7 @@
             main_file_name = args[args.index('-name') + 1]
         except Exception:
             pass
+        line_options = dict(arg[2:].split('=') for arg in args if arg.startswith('--') and '=' in arg)
 
         # Whatever the format we always output the quadruple precision routines
         # to allow for curing possible unstable points.
@@ -463,7 +464,8 @@
             output_type = 'madloop_matchbox'
 
         self._curr_exporter = export_v4.ExportV4Factory(self, \
-                     noclean, output_type=output_type, group_subprocesses=False)
+                     noclean, output_type=output_type, group_subprocesses=False,
+                     cmd_options=line_options)
 
         if self._export_format in ['standalone', 'matchbox']:
             self._curr_exporter.copy_template(self._curr_model)

=== modified file 'madgraph/interface/madevent_interface.py'
--- madgraph/interface/madevent_interface.py	2017-03-06 20:39:01 +0000
+++ madgraph/interface/madevent_interface.py	2017-07-29 19:33:46 +0000
@@ -42,7 +42,8 @@
 
 root_path = os.path.split(os.path.dirname(os.path.realpath( __file__ )))[0]
 root_path = os.path.split(root_path)[0]
-sys.path.insert(0, os.path.join(root_path,'bin'))
+if __name__ == '__main__':
+    sys.path.insert(0, os.path.join(root_path,'bin'))
 
 # usefull shortcut
 pjoin = os.path.join
@@ -52,7 +53,7 @@
  
 try:
     import madgraph
-except ImportError: 
+except ImportError,error: 
     # import from madevent directory
     MADEVENT = True
     import internal.extended_cmd as cmd
@@ -5938,7 +5939,7 @@
                                                                      
         # Now that we know in which mode we are check that all the card
         #exists (copy default if needed)
-
+    
         cards = ['param_card.dat', 'run_card.dat']
         if switch['shower'] in ['PY6', 'PYTHIA6']:
             cards.append('pythia_card.dat')
@@ -6732,9 +6733,9 @@
             if '--web' in args:
                 i = args.index('--web') 
                 args.pop(i)                                                                                                                                                                     
-                cmd_line = MadEventCmd(force_run=True)
+                cmd_line = MadEventCmd(os.path.dirname(root_path),force_run=True)
             else:
-                cmd_line = MadEventCmdShell(force_run=True)
+                cmd_line = MadEventCmdShell(os.path.dirname(root_path),force_run=True)
             if not hasattr(cmd_line, 'do_%s' % args[0]):
                 if parser_error:
                     print parser_error

=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2017-03-07 15:43:44 +0000
+++ madgraph/interface/madgraph_interface.py	2017-07-29 19:33:46 +0000
@@ -7288,6 +7288,7 @@
         if '--postpone_model' in args:
             flaglist.append('store_model')
         
+        line_options = dict(arg[2:].split('=') for arg in args if arg.startswith('--') and '=' in arg)
         main_file_name = ""
         try:
             main_file_name = args[args.index('-name') + 1]
@@ -7416,9 +7417,11 @@
         #Exporter + Template
         if options['exporter'] == 'v4':
             self._curr_exporter = export_v4.ExportV4Factory(self, noclean, 
-                                             group_subprocesses=group_processes)
+                                             group_subprocesses=group_processes,
+                                             cmd_options=line_options)
         elif options['exporter'] == 'cpp':
-            self._curr_exporter = export_cpp.ExportCPPFactory(self, group_subprocesses=group_processes)
+            self._curr_exporter = export_cpp.ExportCPPFactory(self, group_subprocesses=group_processes,
+                                                              cmd_options=line_options)
         
         self._curr_exporter.pass_information_from_cmd(self)
         

=== modified file 'madgraph/interface/reweight_interface.py'
--- madgraph/interface/reweight_interface.py	2017-03-03 10:07:55 +0000
+++ madgraph/interface/reweight_interface.py	2017-07-29 19:33:46 +0000
@@ -66,6 +66,7 @@
     
     prompt = 'Reweight>'
     debug_output = 'Reweight_debug'
+    rwgt_dir_possibility = ['rw_me','rw_me_second','rw_mevirt','rw_mevirt_second']
     
     @misc.mute_logger()
     def __init__(self, event_path=None, allow_madspin=False, mother=None, *completekey, **stdin):
@@ -89,8 +90,10 @@
 
         self.events_file = None
         self.processes = {}
+        self.f2pylib = {}
         self.second_model = None
         self.second_process = None
+        self.dedicated_path = {}
         self.mg5cmd = master_interface.MasterCmd()
         if mother:
             self.mg5cmd.options.update(mother.options)
@@ -198,6 +201,7 @@
             raise Exception, msg
         process, option = mg_interface.MadGraphCmd.split_process_line(process)
         self.proc_option = option
+        self.is_decay = len(process.split('>',1)[0].split()) == 1 
         
         logger.info("process: %s" % process)
         logger.info("options: %s" % option)
@@ -242,7 +246,7 @@
                         r='QCD<=%i' % (int(ior[1])+1)
                     process=process+r+' '
             #handle special tag $ | / @
-            result = re.split('([/$@]|\w+(?:^2)?(?:=|<=|>)?\w+)', process, 1)                    
+            result = re.split('([/$@]|\w+(?:^2)?(?:=|<=|>)+\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.replace(' ','') ,split, rest, final)
@@ -344,6 +348,8 @@
                 self.second_process.append(" ".join(args[1:-1]))
             else:
                 self.second_process = [" ".join(args[1:])]
+        elif args[0] in ['virtual_path', 'tree_path']:
+            self.dedicated_path[args[0]] = args[1]
         elif args[0] == "output":
             if args[1] in ['default', '2.0', 'unweight']:
                 self.output_type = args[1]
@@ -429,36 +435,244 @@
 
         model_line = self.banner.get('proc_card', 'full_model_line')
 
-        if not self.has_standalone_dir:                
+        if not self.has_standalone_dir:                           
             if self.rwgt_dir and os.path.exists(pjoin(self.rwgt_dir,'rw_me','rwgt.pkl')):
                 self.load_from_pickle()
                 if not self.rwgt_dir:
                     self.me_dir = self.rwgt_dir
+                self.load_module()
             elif self.multicore == 'wait':
+                i=0
                 while not os.path.exists(pjoin(self.me_dir,'rw_me','rwgt.pkl')):
-                    time.sleep(60)
+                    time.sleep(10+i)
+                    i+=5
                     print 'wait for pickle'                  
                 print "loading from pickle"
                 if not self.rwgt_dir:
                     self.rwgt_dir = self.me_dir
                 self.load_from_pickle(keep_name=True)
+                self.load_module()
             else:
-                self.create_standalone_directory()  
+                self.create_standalone_directory()
+                self.compile()
+                self.load_module()  
                 if self.multicore == 'create':
+                    self.load_module()
                     if not self.rwgt_dir:
                         self.rwgt_dir = self.me_dir
                     self.save_to_pickle()      
         
-        if self.rwgt_dir:
-            path_me =self.rwgt_dir
-        else:
-            path_me = self.me_dir 
-            
-        if self.second_model or self.second_process:
-            rw_dir = pjoin(path_me, 'rw_me_second')
-        else:
-            rw_dir = pjoin(path_me, 'rw_me')
-
+                    
+        type_rwgt = self.get_weight_names()
+        param_card_iterator, tag_name, output2 = self.handle_param_card(model_line, args, type_rwgt)
+        
+        if self.rwgt_dir:
+            path_me =self.rwgt_dir
+        else:
+            path_me = self.me_dir 
+            
+        if self.second_model or self.second_process or self.dedicated_path:
+            rw_dir = pjoin(path_me, 'rw_me_second')
+        else:
+            rw_dir = pjoin(path_me, 'rw_me')
+        
+        
+        start = time.time()
+        cross, ratio, ratio_square,error = {},{},{}, {}
+        for name in type_rwgt + ['orig']:
+            cross[name], error[name] = 0.,0.
+            ratio[name],ratio_square[name] = 0., 0.# to compute the variance and associate error
+
+        if self.output_type == "default":
+            output = open( self.lhe_input.name +'rw', 'w')
+            #write the banner to the output file
+            self.banner.write(output, close_tag=False)
+        else:
+            output = {}
+            for name in type_rwgt:
+                output[name] = open( self.lhe_input.name +'rw'+name, 'w')
+                #write the banner to the output file
+                self.banner.write(output[name], close_tag=False)
+        
+                
+        if self.lhe_input.closed:
+            self.lhe_input = lhe_parser.EventFile(self.lhe_input.name)
+
+        self.lhe_input.seek(0)
+        for event_nb,event in enumerate(self.lhe_input):
+            #control logger
+            if (event_nb % max(int(10**int(math.log10(float(event_nb)+1))),10)==0): 
+                    running_time = misc.format_timer(time.time()-start)
+                    logger.info('Event nb %s %s' % (event_nb, running_time))
+            if (event_nb==10001): logger.info('reducing number of print status. Next status update in 10000 events')
+
+
+            weight = self.calculate_weight(event)
+            if not isinstance(weight, dict):
+                weight = {'':weight}
+            
+            for name in weight:
+                cross[name] += weight[name]
+                ratio[name] += weight[name]/event.wgt
+                ratio_square[name] += (weight[name]/event.wgt)**2
+                
+            # ensure to have a consistent order of the weights. new one are put 
+            # at the back, remove old position if already defines
+            for tag in type_rwgt:
+                try:
+                    event.reweight_order.remove('%s%s'  % (tag_name,tag))
+                except ValueError:
+                    continue
+            
+            event.reweight_order += ['%s%s' % (tag_name,name) for name in type_rwgt]  
+            if self.output_type == "default":
+                for name in weight:
+                    if 'orig' in name:
+                        continue             
+                    event.reweight_data['%s%s' % (tag_name,name)] = weight[name]
+                    #write this event with weight
+                output.write(str(event))
+                if self.mother:
+                    event.wgt = weight[type_rwgt[0]]
+                    event.reweight_data = {}
+                    output2.write(str(event))
+            else:
+                for i,name in enumerate(weight):
+                    event.wgt = weight[name]
+                    event.reweight_data = {}
+                    if self.mother and len(weight)==1:
+                        output2.write(str(event))
+                    elif self.mother and i == 0:
+                        output[name].write(str(event))
+                        output2.write(str(event))
+                    else:
+                        output[name].write(str(event))
+                
+        # check normalisation of the events:
+        if 'event_norm' in self.run_card:
+            if self.run_card['event_norm'] == 'average':
+                for key, value in cross.items():
+                    cross[key] = value / (event_nb+1)
+                
+        running_time = misc.format_timer(time.time()-start)
+        logger.info('All event done  (nb_event: %s) %s' % (event_nb+1, running_time))        
+        
+        
+        if self.output_type == "default":
+            output.write('</LesHouchesEvents>\n')
+            output.close()
+        else:
+            for key in output:
+                output[key].write('</LesHouchesEvents>\n')
+                output.close()
+        
+        if self.mother:
+            output2.write('</LesHouchesEvents>\n')
+            output2.close()        
+            # add output information
+            if hasattr(self.mother, 'results'):
+                run_name = self.mother.run_name
+                results = self.mother.results
+                results.add_run(run_name, self.run_card, current=True)
+                results.add_detail('nb_event', event_nb+1)
+                name = type_rwgt[0]
+                results.add_detail('cross', cross[name])
+                event_nb +=1
+                for name in type_rwgt:
+                    variance = ratio_square[name]/event_nb - (ratio[name]/event_nb)**2
+                    orig_cross, orig_error = self.orig_cross
+                    error[name] = variance/math.sqrt(event_nb) * orig_cross + ratio[name]/event_nb * orig_error
+                results.add_detail('error', error[type_rwgt[0]])
+                import madgraph.interface.madevent_interface as ME_interface
+                if isinstance(self.mother, ME_interface.MadEventCmd):
+                    self.mother.create_plot(mode='reweight', event_path=output2.name,
+                                        tag=self.run_card['run_tag'])
+                #modify the html output to add the original run
+                if 'plot' in results.current.reweight:
+                    html_dir = pjoin(self.mother.me_dir, 'HTML', run_name)
+                    td = pjoin(self.mother.options['td_path'], 'td') 
+                    MA = pjoin(self.mother.options['madanalysis_path'])
+                    path1 = pjoin(html_dir, 'plots_parton')
+                    path2 = pjoin(html_dir, 'plots_%s' % self.run_card['run_tag'])
+                    outputplot = path2
+                    combine_plots.merge_all_plots(path2, path1, outputplot, td, MA)
+                #results.update_status(level='reweight')
+                #results.update(status, level, makehtml=True, error=False)
+                
+                #old_name = self.mother.results.current['run_name']
+                #new_run = '%s_rw_%s' % (old_name, rewgtid)
+                #self.mother.results.add_run( new_run, self.run_card)
+                #self.mother.results.add_detail('nb_event', event_nb+1)
+                #self.mother.results.add_detail('cross', cross)
+                #self.mother.results.add_detail('error', 'nan')
+                #self.mother.do_plot('%s -f' % new_run)
+                #self.mother.update_status('Reweight %s done' % rewgtid, 'madspin')
+                #self.mother.results.def_current(old_name)
+                #self.run_card['run_tag'] = self.run_card['run_tag'][9:]
+                #self.mother.run_name = old_name
+        self.lhe_input.close()
+        if not self.mother or self.output_type != "default" :
+            target = pjoin(self.mother.me_dir, 'Events', run_name, 'events.lhe')
+        else:
+            target = self.lhe_input.name
+        
+        if self.output_type == "default":
+            files.mv(output.name, target)
+            logger.info('Event %s have now the additional weight' % self.lhe_input.name)
+        elif self.output_type == "unweight":
+            output2.close()
+            lhe = lhe_parser.EventFile(output2.name)
+            nb_event = lhe.unweight(target)
+            if self.mother and  hasattr(self.mother, 'results'):
+                results = self.mother.results
+                results.add_detail('nb_event', nb_event)
+                results.current.parton.append('lhe')
+            logger.info('Event %s is now unweighted under the new theory' % output2.name)                
+        else:
+            files.mv(output2.name, self.lhe_input.name)     
+            if self.mother and  hasattr(self.mother, 'results'):
+                results = self.mother.results
+                results.current.parton.append('lhe')       
+            logger.info('Event %s is now created with new central weight' % output2.name)
+        
+        if self.multicore != 'create':
+            for name in cross:
+                if name == 'orig':
+                    continue
+                logger.info('new cross-section is %s: %g pb (indicative error: %g pb)' %\
+                        ('(%s)' %name if name else '',cross[name], error[name]))
+            
+        self.terminate_fortran_executables(new_card_only=True)
+        #store result
+        for name in cross:
+            if name == 'orig':
+                self.all_cross_section[name] = (cross[name], error[name])
+            else:
+                self.all_cross_section[(tag_name,name)] = (cross[name], error[name])
+        
+        # perform the scanning
+        if param_card_iterator:
+            for i,card in enumerate(param_card_iterator):
+                if self.options['rwgt_name']:
+                    self.options['rwgt_name'] = '%s_%s' % (self.options['rwgt_name'].rsplit('_',1)[0], i+1)
+                card.write(pjoin(rw_dir, 'Cards', 'param_card.dat'))
+                self.exec_cmd("launch --keep_card", printcmd=False, precmd=True)
+        
+        self.options['rwgt_name'] = None
+    
+    def handle_param_card(self, model_line, args, type_rwgt):
+        
+        if self.rwgt_dir:
+            path_me =self.rwgt_dir
+        else:
+            path_me = self.me_dir 
+            
+        if self.second_model or self.second_process or self.dedicated_path:
+            rw_dir = pjoin(path_me, 'rw_me_second')
+        else:
+            rw_dir = pjoin(path_me, 'rw_me')
+        
+        
         if not '--keep_card' in args:
             ff = open(pjoin(rw_dir,'Cards', 'param_card.dat'), 'w')
             ff.write(self.banner['slha'])
@@ -476,7 +690,7 @@
             self.stored_line = None
         
         # get the names of type of reweighting requested
-        type_rwgt = self.get_weight_names()
+        #type_rwgt = self.get_weight_names()
 
         # check for potential scan in the new card 
         new_card = open(pjoin(rw_dir, 'Cards', 'param_card.dat')).read()
@@ -499,7 +713,8 @@
             param_card_iterator = main_card
             first_card = param_card_iterator.next(autostart=True)
             new_card = first_card.write()
-            first_card.write(pjoin(rw_dir, 'Cards', 'param_card.dat'))                
+            first_card.write(pjoin(rw_dir, 'Cards', 'param_card.dat'))  
+                          
         # check if "Auto" is present for a width parameter
         if "auto" in new_card.lower():            
             self.mother.check_param_card(pjoin(rw_dir, 'Cards', 'param_card.dat'))
@@ -556,16 +771,15 @@
         else:
             tag = str(rewgtid)
         
-        if not self.second_model:
+        if not self.second_model and not self.dedicated_path:
             old_param = check_param_card.ParamCard(s_orig.splitlines())
             new_param =  check_param_card.ParamCard(s_new.splitlines())
             card_diff = old_param.create_diff(new_param)
             if card_diff == '' and not self.second_process:
-                if not __debug__:
-                    logger.warning(' REWEIGHTING: original card and new card are identical. Bypass this run')
-                    return
-                else:
-                    logger.warning(' REWEIGHTING: original card and new card are identical. Run it due to debug mode')
+                    logger.warning(' REWEIGHTING: original card and new card are identical.')
+#                    return
+#                else:
+#                    logger.warning(' REWEIGHTING: original card and new card are identical. Run it due to debug mode')
                 #raise self.InvalidCmd, 'original card and new card are identical'
             try:
                 if old_param['sminputs'].get(3)- new_param['sminputs'].get(3) > 1e-3 * new_param['sminputs'].get(3):
@@ -581,9 +795,15 @@
                 for name in type_rwgt:
                     mg_rwgt_info.append((tag, name, str_proc + '\n'+ card_diff))
         else:
-            str_info = "change model %s" % self.second_model
+            if self.second_model:
+                str_info = "change model %s" % self.second_model
+            else:
+                str_info =''
             if self.second_process:
                 str_info += "\n change process  ".join([""]+self.second_process)
+            if self.dedicated_path:
+                for k,v in self.dedicated_path.items():
+                    str_info += "\n change %s %s" % (k,v)
             card_diff = str_info
             str_info += '\n' + s_new
             for name in type_rwgt:
@@ -602,23 +822,6 @@
         self.banner['initrwgt'] = self.banner['initrwgt'].replace('\n\n', '\n')
 
 
-        start = time.time()
-        cross, ratio, ratio_square,error = {},{},{}, {}
-        for name in type_rwgt + ['orig']:
-            cross[name], error[name] = 0.,0.
-            ratio[name],ratio_square[name] = 0., 0.# to compute the variance and associate error
-
-        if self.output_type == "default":
-            output = open( self.lhe_input.name +'rw', 'w')
-            #write the banner to the output file
-            self.banner.write(output, close_tag=False)
-        else:
-            output = {}
-            for name in type_rwgt:
-                output[name] = open( self.lhe_input.name +'rw'+name, 'w')
-                #write the banner to the output file
-                self.banner.write(output[name], close_tag=False)
-        
         logger.info('starts to compute weight for events with the following modification to the param_card:')
         logger.info(card_diff.replace('\n','\nKEEP:'))
         # prepare the output file for the weight plot
@@ -640,195 +843,26 @@
                           (self.mother.run_name, self.run_card['run_tag'])),'w') 
             new_banner.write(ff)
             ff.close()
-        
+
         # Loop over all events
         if self.options['rwgt_name']:
             tag_name = self.options['rwgt_name']
         else:
             tag_name = 'rwgt_%s' % rewgtid
-                
-        os.environ['GFORTRAN_UNBUFFERED_ALL'] = 'y'
-        if self.lhe_input.closed:
-            self.lhe_input = lhe_parser.EventFile(self.lhe_input.name)
-
-#        Multicore option not really stable -> not use it
-        nb_core = 1
-#        if nb_core >1:
-#            multicore = cluster.MultiCore(nb_core)
-
-        self.lhe_input.seek(0)
-        for event_nb,event in enumerate(self.lhe_input):
-            #control logger
-            if (event_nb % max(int(10**int(math.log10(float(event_nb)+1))),10)==0): 
-                    running_time = misc.format_timer(time.time()-start)
-                    logger.info('Event nb %s %s' % (event_nb, running_time))
-            if (event_nb==10001): logger.info('reducing number of print status. Next status update in 10000 events')
-
-            if nb_core > 1:
-                #        Multicore option not really stable -> not use it
-                while 1:
-                    if multicore.queue.qsize() < 100 * nb_core:
-                        multicore.submit(self.write_reweighted_event, argument=[event, tag_name])
-                        break
-                    #else:
-                    #    time.sleep(0.001)
-                continue
-            else:
-                weight = self.calculate_weight(event)
-                if not isinstance(weight, dict):
-                    weight = {'':weight}
-                
-                for name in weight:
-                    cross[name] += weight[name]
-                    ratio[name] += weight[name]/event.wgt
-                    ratio_square[name] += (weight[name]/event.wgt)**2
-                    
-                # ensure to have a consistent order of the weights. new one are put 
-                # at the back, remove old position if already defines
-                for tag in type_rwgt:
-                    try:
-                        event.reweight_order.remove('%s%s'  % (tag_name,tag))
-                    except ValueError:
-                        continue
-                
-                event.reweight_order += ['%s%s' % (tag_name,name) for name in type_rwgt]  
-                if self.output_type == "default":
-                    for name in weight:
-                        if 'orig' in name:
-                            continue             
-                        event.reweight_data['%s%s' % (tag_name,name)] = weight[name]
-                        #write this event with weight
-                    output.write(str(event))
-                    if self.mother:
-                        event.wgt = weight[type_rwgt[0]]
-                        event.reweight_data = {}
-                        output2.write(str(event))
+
+
+        #initialise module.
+        
+        for (path,tag), module in self.f2pylib.items():
+            with misc.chdir(pjoin(os.path.dirname(rw_dir), path)):
+#            with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                if 'second' in path or tag == 3:
+                    module.initialise(pjoin(rw_dir, 'Cards', 'param_card.dat'))
                 else:
-                    for i,name in enumerate(weight):
-                        event.wgt = weight[name]
-                        event.reweight_data = {}
-                        if self.mother and len(weight)==1:
-                            output2.write(str(event))
-                        elif self.mother and i == 0:
-                            output[name].write(str(event))
-                            output2.write(str(event))
-                        else:
-                            output[name].write(str(event))
-                
-        # check normalisation of the events:
-        if 'event_norm' in self.run_card:
-            if self.run_card['event_norm'] == 'average':
-                for key, value in cross.items():
-                    cross[key] = value / (event_nb+1)
-                
-        running_time = misc.format_timer(time.time()-start)
-        logger.info('All event done  (nb_event: %s) %s' % (event_nb+1, running_time))        
-        
-        
-        if self.output_type == "default":
-            output.write('</LesHouchesEvents>\n')
-            output.close()
-        else:
-            for key in output:
-                output[key].write('</LesHouchesEvents>\n')
-                output.close()
-            
-        os.environ['GFORTRAN_UNBUFFERED_ALL'] = 'n'
-        
-        if self.mother:
-            output2.write('</LesHouchesEvents>\n')
-            output2.close()        
-            # add output information
-            if hasattr(self.mother, 'results'):
-                run_name = self.mother.run_name
-                results = self.mother.results
-                results.add_run(run_name, self.run_card, current=True)
-                results.add_detail('nb_event', event_nb+1)
-                name = type_rwgt[0]
-                results.add_detail('cross', cross[name])
-                event_nb +=1
-                for name in type_rwgt:
-                    variance = ratio_square[name]/event_nb - (ratio[name]/event_nb)**2
-                    orig_cross, orig_error = self.orig_cross
-                    error[name] = variance/math.sqrt(event_nb) * orig_cross + ratio[name]/event_nb * orig_error
-                results.add_detail('error', error[type_rwgt[0]])
-                import madgraph.interface.madevent_interface as ME_interface
-                if isinstance(self.mother, ME_interface.MadEventCmd):
-                    self.mother.create_plot(mode='reweight', event_path=output2.name,
-                                        tag=self.run_card['run_tag'])
-                #modify the html output to add the original run
-                if 'plot' in results.current.reweight:
-                    html_dir = pjoin(self.mother.me_dir, 'HTML', run_name)
-                    td = pjoin(self.mother.options['td_path'], 'td') 
-                    MA = pjoin(self.mother.options['madanalysis_path'])
-                    path1 = pjoin(html_dir, 'plots_parton')
-                    path2 = pjoin(html_dir, 'plots_%s' % self.run_card['run_tag'])
-                    outputplot = path2
-                    combine_plots.merge_all_plots(path2, path1, outputplot, td, MA)
-                #results.update_status(level='reweight')
-                #results.update(status, level, makehtml=True, error=False)
-                
-                #old_name = self.mother.results.current['run_name']
-                #new_run = '%s_rw_%s' % (old_name, rewgtid)
-                #self.mother.results.add_run( new_run, self.run_card)
-                #self.mother.results.add_detail('nb_event', event_nb+1)
-                #self.mother.results.add_detail('cross', cross)
-                #self.mother.results.add_detail('error', 'nan')
-                #self.mother.do_plot('%s -f' % new_run)
-                #self.mother.update_status('Reweight %s done' % rewgtid, 'madspin')
-                #self.mother.results.def_current(old_name)
-                #self.run_card['run_tag'] = self.run_card['run_tag'][9:]
-                #self.mother.run_name = old_name
-        self.lhe_input.close()
-        if not self.mother or self.output_type != "default" :
-            target = pjoin(self.mother.me_dir, 'Events', run_name, 'events.lhe')
-        else:
-            target = self.lhe_input.name
-        
-        if self.output_type == "default":
-            files.mv(output.name, target)
-            logger.info('Event %s have now the additional weight' % self.lhe_input.name)
-        elif self.output_type == "unweight":
-            output2.close()
-            lhe = lhe_parser.EventFile(output2.name)
-            nb_event = lhe.unweight(target)
-            if self.mother and  hasattr(self.mother, 'results'):
-                results = self.mother.results
-                results.add_detail('nb_event', nb_event)
-                results.current.parton.append('lhe')
-            logger.info('Event %s is now unweighted under the new theory' % output2.name)                
-        else:
-            files.mv(output2.name, self.lhe_input.name)     
-            if self.mother and  hasattr(self.mother, 'results'):
-                results = self.mother.results
-                results.current.parton.append('lhe')       
-            logger.info('Event %s is now created with new central weight' % output2.name)
-        
-        if self.multicore != 'create':
-            for name in cross:
-                if name == 'orig':
-                    continue
-                logger.info('new cross-section is %s: %g pb (indicative error: %g pb)' %\
-                        ('(%s)' %name if name else '',cross[name], error[name]))
-            
-        self.terminate_fortran_executables(new_card_only=True)
-        #store result
-        for name in cross:
-            if name == 'orig':
-                self.all_cross_section[name] = (cross[name], error[name])
-            else:
-                self.all_cross_section[(tag_name,name)] = (cross[name], error[name])
-        
-        # perform the scanning
-        if param_card_iterator:
-            for i,card in enumerate(param_card_iterator):
-                if self.options['rwgt_name']:
-                    self.options['rwgt_name'] = '%s_%s' % (self.options['rwgt_name'].rsplit('_',1)[0], i+1)
-                card.write(pjoin(rw_dir, 'Cards', 'param_card.dat'))
-                self.exec_cmd("launch --keep_card", printcmd=False, precmd=True)
-        
-        self.options['rwgt_name'] = None
+                    module.initialise(pjoin(path_me, 'rw_me', 'Cards', 'param_card_orig.dat'))
+        return param_card_iterator, tag_name, output2
     
+        
     def do_set(self, line):
         "Not in help"
         
@@ -882,19 +916,18 @@
     def do_compute_widths(self, line):
         return self.mother.do_compute_widths(line)
     
-    def calculate_weight(self, event, space=None):
+    def calculate_weight(self, event):
         """space defines where to find the calculator (in multicore)"""
         
         if self.has_nlo and self.rwgt_mode != "LO":
-            return self.calculate_nlo_weight(event, space)
+            return self.calculate_nlo_weight(event)
         
-        if not space: 
-            space = self
         event.parse_reweight()                    
            
         # LO reweighting    
-        w_orig = self.calculate_matrix_element(event, 0, space)
-        w_new =  self.calculate_matrix_element(event, 1, space)
+        w_orig = self.calculate_matrix_element(event, 0)
+        w_new =  self.calculate_matrix_element(event, 1)
+        
         if w_orig == 0:
             tag, order = event.get_tag_and_order()
             orig_order, Pdir, hel_dict = self.id_to_path[tag]
@@ -913,14 +946,11 @@
 
         return {'orig': event.wgt, '': w_new/w_orig*event.wgt}
      
-    def calculate_nlo_weight(self, event, space=None):
+    def calculate_nlo_weight(self, event):
 
 
         type_nlo = self.get_weight_names()
         final_weight = {'orig': event.wgt}
-        
-        if not space: 
-            space = self #for multicore: not use so far
             
         event.parse_reweight()
         event.parse_nlo_weight() 
@@ -944,20 +974,20 @@
             if '_nlo' in type_nlo and any(c in all_ctype for c in [2,14,15]):
                 need_V =True
                 
-            w_orig = self.calculate_matrix_element(cevent, 0, space)
-            w_new =  self.calculate_matrix_element(cevent, 1, space)
+            w_orig = self.calculate_matrix_element(cevent, 0)
+            w_new =  self.calculate_matrix_element(cevent, 1)
             ratio_T = w_new/w_orig
             if need_V:
                 scale2 = cevent.wgts[0].scales2[0]
                 #for scale2 in set(c.scales2[1] for c in cevent.wgts): 
-                w_origV = self.calculate_matrix_element(cevent, 'V0', space, scale2=scale2)
-                w_newV =  self.calculate_matrix_element(cevent, 'V1', space, scale2=scale2)                    
+                w_origV = self.calculate_matrix_element(cevent, 'V0', scale2=scale2)
+                w_newV =  self.calculate_matrix_element(cevent, 'V1', scale2=scale2)                    
                 ratio_BV = (w_newV + w_new) / (w_origV + w_orig)
                 ratio_V = w_newV/w_origV
             else:
                 ratio_V = "should not be used"
                 ratio_BV = "should not be used"
-                
+
             for c_wgt in cevent.wgts:
                 orig_wgt += c_wgt.ref_wgt
                 #add the information to the input
@@ -1018,8 +1048,8 @@
             final_weight['_tree'] = new_out/orig_wgt*event.wgt            
              
         if '_lo' in type_nlo:
-            w_orig = self.calculate_matrix_element(event, 0, space)
-            w_new =  self.calculate_matrix_element(event, 1, space)            
+            w_orig = self.calculate_matrix_element(event, 0)
+            w_new =  self.calculate_matrix_element(event, 1)            
             final_weight['_lo'] = w_new/w_orig*event.wgt
             
         return final_weight 
@@ -1045,7 +1075,7 @@
             open(pjoin(Pdir, 'matrix%spy.so' % tag),'w').write(open(pjoin(Pdir, 'matrix2py.so')
                                         ).read().replace('matrix2py', 'matrix%spy' % tag))
     
-    def calculate_matrix_element(self, event, hypp_id, space, scale2=0):
+    def calculate_matrix_element(self, event, hypp_id, scale2=0):
         """routine to return the matrix element"""
         
         if self.has_nlo:
@@ -1057,130 +1087,23 @@
         if isinstance(hypp_id, str) and hypp_id.startswith('V'):
             tag = (tag,'V')
             hypp_id = int(hypp_id[1:])
-            base = "rw_mevirt"
-        else:
-            base = "rw_me"
+        #    base = "rw_mevirt"
+        #else:
+        #    base = "rw_me"
         
-        if (not self.second_model and not self.second_process) or hypp_id==0:
+        if (not self.second_model and not self.second_process and not self.dedicated_path) or hypp_id==0:
             orig_order, Pdir, hel_dict = self.id_to_path[tag]
         else:
             orig_order, Pdir, hel_dict = self.id_to_path_second[tag] 
             
-
-        run_id = (tag, hypp_id)
-
-        assert space == self
-        start = False
-        if run_id in space.calculator:
-            external = space.calculator[run_id]
-            #    mod = space.calculator[(run_id,'module')]
-            #    #with misc.chdir(Pdir):
-            #    #    if hypp_id==0:
-            #    #        mod.initialise('param_card_orig.dat')
-            #    #    else:
-            #    #        mod.initialise('param_card.dat')
-        elif (not self.second_model and not self.second_process) or hypp_id==0:
-            # create the executable for this param_card  
-            subdir = pjoin(self.me_dir, base, 'SubProcesses')
-            if self.me_dir not in sys.path:
-                sys.path.insert(0,self.me_dir)
-            if self.rwgt_dir and self.rwgt_dir not in sys.path:
-                sys.path.insert(0,self.rwgt_dir)
-            Pname = os.path.basename(Pdir)
-            if hypp_id == 0:
-                if (Pdir, 0) not in dir_to_f2py_free_mod:
-                    metag = 1
-                    dir_to_f2py_free_mod[(Pdir,0)] = (metag, nb_f2py_module)
-                else:
-                    metag, old_module = dir_to_f2py_free_mod[(Pdir,0)]
-                    if old_module != nb_f2py_module:
-                        metag += 1
-                        dir_to_f2py_free_mod[(Pdir,0)] = (metag, nb_f2py_module)
-                os.environ['MENUM'] = '2'
-                if not self.rwgt_dir or not os.path.exists(pjoin(Pdir, 'matrix2py.so')):
-                    misc.multiple_try(nb_retry,sleep)(misc.compile)(['matrix2py.so'], cwd=Pdir)                
-                self.rename_f2py_lib(Pdir, 2*metag)
-                try:
-                    mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, 2*metag), globals(), locals(), [],-1)
-                except:
-                    import platform
-                    if platform.system() == 'Darwin':
-                        os.system('install_name_tool -change libMadLoop.dylib %s/libMadLoop.dylib matrix%spy.so' % (Pdir,2*metag))
-                        mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, 2*metag), globals(), locals(), [],-1)
-                    else:
-                        misc.sprint("fail compilation")
-                        raise
-                S = mymod.SubProcesses
-                P = getattr(S, Pname)
-                mymod = getattr(P, 'matrix%spy' % (2*metag))
-                with misc.chdir(Pdir):
-                    with misc.stdchannel_redirected(sys.stdout, os.devnull):
-                        mymod.initialise('param_card_orig.dat')
-                 
-                     
-            if hypp_id == 1:
-                #incorrect line
-                metag = dir_to_f2py_free_mod[(Pdir,0)][0]
-                newtag = 2*metag+1
-                self.rename_f2py_lib(Pdir, newtag)
-                try:
-                    mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, newtag), globals(), locals(), [],-1)
-                except Exception, error:
-                    newtag  = "L%s" % newtag
-                    os.environ['MENUM'] = newtag
-                    misc.multiple_try(nb_retry,sleep)(misc.compile)(['matrix%spy.so' % newtag], cwd=Pdir)
-                    mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, newtag), globals(), locals(), [],-1)
-                 
-                S = mymod.SubProcesses
-                P = getattr(S, Pname)
-                mymod = getattr(P, 'matrix%spy' % newtag)
-                with misc.chdir(Pdir):
-                    with misc.stdchannel_redirected(sys.stdout, os.devnull):
-                        mymod.initialise('param_card.dat') 
-                     
-            space.calculator[run_id] = mymod.get_me
-            space.calculator[(run_id,'module')] = mymod
-            external = space.calculator[run_id]                      
+        base = os.path.basename(os.path.dirname(Pdir))
+        if '_second' in base:
+            moduletag = (base, 2)
         else:
-            subdir = pjoin(self.me_dir,'%s_second' % base, 'SubProcesses')
-            if self.me_dir not in sys.path:
-                sys.path.append(self.me_dir)
- 
-            assert hypp_id == 1
-            Pname = os.path.basename(Pdir)
-            os.environ['MENUM'] = '2'
-            if not self.rwgt_dir or not os.path.exists(pjoin(Pdir, 'matrix2py.so')):
-                misc.multiple_try(nb_retry,sleep)(misc.compile)(['matrix2py.so'], cwd=pjoin(subdir, Pdir))
-            if (Pdir, 1) not in dir_to_f2py_free_mod:
-                metag = 1
-                dir_to_f2py_free_mod[(Pdir,1)] = (metag, nb_f2py_module)
-            else:
-                metag, old_module = dir_to_f2py_free_mod[(Pdir,1)]
-                if old_module != nb_f2py_module:
-                    metag += 1
-                    dir_to_f2py_free_mod[(Pdir,1)] = (metag, nb_f2py_module)
-            self.rename_f2py_lib(Pdir, metag)
-            try:
-                mymod = __import__("%s_second.SubProcesses.%s.matrix%spy" % (base, Pname, metag))
-            except ImportError:
-                metag = "L%s" % metag
-                os.environ['MENUM'] = str(metag)
-                misc.multiple_try(nb_retry,sleep)(misc.compile)(['matrix%spy.so' % metag], cwd=pjoin(subdir, Pdir))
-                mymod = __import__("%s_second.SubProcesses.%s.matrix%spy" % (base, Pname, metag))
-                     
-            reload(mymod)
-            S = mymod.SubProcesses
-            P = getattr(S, Pname)
-            mymod = getattr(P, 'matrix%spy' % metag)
-            with misc.chdir(Pdir):   
-                with misc.stdchannel_redirected(sys.stdout, os.devnull):   
-                    mymod.initialise('param_card.dat')              
-            space.calculator[run_id] = mymod.get_me
-            space.calculator[(run_id,'module')] = mymod
-            external = space.calculator[run_id]                
+            moduletag = (base, 2+hypp_id)
         
+        module = self.f2pylib[moduletag]
 
-        
         p = event.get_momenta(orig_order)
         # add helicity information
         
@@ -1192,26 +1115,18 @@
                 pboost = lhe_parser.FourMomentum(p[0]) + lhe_parser.FourMomentum(p[1])
                 for i,thisp in enumerate(p):
                     p[i] = lhe_parser.FourMomentum(thisp).zboost(pboost).get_tuple()
-                assert p[0][1] == p[0][2] == 0 == p[1][2] == p[1][2] == 0 
-
-                
+                assert p[0][1] == p[0][2] == 0 == p[1][2] == p[1][2] == 0                 
         else:
-            nhel = 0
+            nhel = -1
         pold = list(p)
         p = self.invert_momenta(p)
+        pdg = list(orig_order[0])+list(orig_order[1])
 
         with misc.chdir(Pdir):
             with misc.stdchannel_redirected(sys.stdout, os.devnull):
-                if 'V' in tag or \
-                   (hypp_id ==1  and self.second_process and any('sqrvirt' in l for l in self.second_process)):
-                    me_value = external(p,event.aqcd, math.sqrt(scale2), nhel)
-                else:
-                    try:
-                        me_value = external(p,event.aqcd, nhel)
-                    except TypeError:
-                        me_value = external(p,event.aqcd, math.sqrt(scale2), nhel)                    
-
-        # for NLO we have also the stability status code
+                me_value = module.smatrixhel(pdg,p, event.aqcd, scale2, nhel)
+                                
+        # for loop we have also the stability status code
         if isinstance(me_value, tuple):
             me_value, code = me_value
             #if code points unstability -> returns 0
@@ -1277,6 +1192,226 @@
         
 
     @misc.mute_logger()
+    def create_standalone_tree_directory(self, data ,second=False):
+        """generate the various directory for the weight evaluation"""
+        
+        mgcmd = self.mg5cmd         
+        path_me = data['path'] 
+        # 2. compute the production matrix element -----------------------------
+        has_nlo = False  
+        mgcmd.exec_cmd("set group_subprocesses False")
+
+        if not second:
+            logger.info('generating the square matrix element for reweighting')
+        else:
+            logger.info('generating the square matrix element for reweighting (second model and/or processes)')
+        start = time.time()
+        commandline=''
+        for i,proc in enumerate(data['processes']):
+            if '[' not in proc:
+                commandline += "add process %s ;" % proc
+            else:
+                has_nlo = True
+                if self.banner.get('run_card','ickkw') == 3:
+                    if len(proc) == min([len(p.strip()) for p in data['processes']]):
+                        commandline += self.get_LO_definition_from_NLO(proc, self.model)
+                    else:
+                        commandline += self.get_LO_definition_from_NLO(proc,
+                                                     self.model, real_only=True)
+                else:
+                    commandline += self.get_LO_definition_from_NLO(proc, self.model)
+        
+        commandline = commandline.replace('add process', 'generate',1)
+        logger.info(commandline)
+        try:
+            mgcmd.exec_cmd(commandline, precmd=True, errorhandling=False)
+        except diagram_generation.NoDiagramException:
+            commandline=''
+            for proc in data['processes']:
+                if '[' not in proc:
+                    raise
+                # pass to virtsq=
+                base, post = proc.split('[',1)
+                nlo_order, post = post.split(']',1)
+                if '=' not in nlo_order:
+                    nlo_order = 'virt=%s' % nlo_order
+                elif 'noborn' in nlo_order:
+                    nlo_order = nlo_order.replace('noborn', 'virt')
+                commandline += "add process %s [%s] %s;" % (base,nlo_order,post)
+            commandline = commandline.replace('add process', 'generate',1)
+            logger.info("RETRY with %s", commandline)
+            mgcmd.exec_cmd(commandline, precmd=True)
+            has_nlo = False
+        except Exception, error:
+            raise
+        
+        commandline = 'output standalone_rw %s --prefix=int' % pjoin(path_me,data['paths'][0])
+        mgcmd.exec_cmd(commandline, precmd=True)        
+        logger.info('Done %.4g' % (time.time()-start))
+        self.has_standalone_dir = True
+        
+
+        # 3. Store id to directory information ---------------------------------
+        if False:
+            # keep this for debugging
+            matrix_elements = mgcmd._curr_matrix_elements.get_matrix_elements()
+            
+            to_check = [] # list of tag that do not have a Pdir at creation time.
+            for me in matrix_elements:
+                for proc in me.get('processes'):
+                    initial = []    #filled in the next line
+                    final = [l.get('id') for l in proc.get('legs')\
+                          if l.get('state') or initial.append(l.get('id'))]
+                    order = (initial, final)
+                    tag = proc.get_initial_final_ids()
+                    decay_finals = proc.get_final_ids_after_decay()
+    
+                    if tag[1] != decay_finals:
+                        order = (initial, list(decay_finals))
+                        decay_finals.sort()
+                        tag = (tag[0], tuple(decay_finals))
+                    Pdir = pjoin(path_me, data['paths'][0], 'SubProcesses', 
+                                      'P%s' % me.get('processes')[0].shell_string())
+    
+                    if not os.path.exists(Pdir):
+                        to_check.append(tag)
+                        continue                        
+                    if tag in data['id2path']:
+                        if not Pdir == data['id2path'][tag][1]:
+                            misc.sprint(tag, Pdir, data['id2path'][tag][1])
+                            raise self.InvalidCmd, '2 different process have the same final states. This module can not handle such situation'
+                        else:
+                            continue
+                    # build the helicity dictionary
+                    hel_nb = 0
+                    hel_dict = {9:0} # unknown helicity -> use full ME
+                    for helicities in me.get_helicity_matrix():
+                        hel_nb +=1 #fortran starts at 1
+                        hel_dict[tuple(helicities)] = hel_nb
+    
+                    data['id2path'][tag] = [order, Pdir, hel_dict]        
+     
+            for tag in to_check:
+                if tag not in self.id_to_path:
+                    logger.warning("no valid path for %s" % (tag,))
+                    #raise self.InvalidCmd, "no valid path for %s" % (tag,)
+        
+        # 4. Check MadLoopParam for Loop induced
+        if os.path.exists(pjoin(path_me, data['paths'][0], 'Cards', 'MadLoopParams.dat')):
+            MLCard = banner.MadLoopParam(pjoin(path_me, data['paths'][0], 'Cards', 'MadLoopParams.dat'))
+            MLCard.set('WriteOutFilters', False)
+            MLCard.set('UseLoopFilter', False)
+            MLCard.set("DoubleCheckHelicityFilter", False)
+            MLCard.set("HelicityFilterLevel", 0)
+            MLCard.write(pjoin(path_me, data['paths'][0], 'SubProcesses', 'MadLoopParams.dat'),
+                         pjoin(path_me, data['paths'][0], 'Cards', 'MadLoopParams.dat'), 
+                         commentdefault=False)
+            
+            #if self.multicore == 'create':
+            #    print "compile OLP", data['paths'][0]
+            #    misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][0],'SubProcesses'),
+            #                 nb_core=self.mother.options['nb_core'])
+        
+        if os.path.exists(pjoin(path_me, data['paths'][1], 'Cards', 'MadLoopParams.dat')):
+            if self.multicore == 'create':
+                print "compile OLP", data['paths'][1]
+                # It is potentially unsafe to use several cores, We limit ourself to one for now
+                # n_cores = self.mother.options['nb_core']
+                n_cores = 1
+                misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][1],'SubProcesses'),
+                             nb_core=self.mother.options['nb_core'])
+                
+        return has_nlo
+
+                
+    @misc.mute_logger()
+    def create_standalone_virt_directory(self, data ,second=False):
+        """generate the various directory for the weight evaluation"""
+                
+        mgcmd = self.mg5cmd
+        path_me = data['path'] 
+        # Do not pass here for LO/NLO_tree
+        start = time.time()
+        commandline=''
+        for proc in data['processes']:
+            if '[' not in proc:
+                pass
+            else:
+                proc = proc.replace('[', '[ virt=')
+                commandline += "add process %s ;" % proc
+        # deactivate golem since it creates troubles
+        old_options = dict(mgcmd.options)
+        if mgcmd.options['golem'] or mgcmd.options['pjfry']:
+            logger.info(" When doing NLO reweighting, MG5aMC cannot use the loop reduction algorithms Golem and/or PJFry++")
+        mgcmd.options['golem'] = None            
+        mgcmd.options['pjfry'] = None 
+        commandline = commandline.replace('add process', 'generate',1)
+        logger.info(commandline)
+        mgcmd.exec_cmd(commandline, precmd=True)
+        commandline = 'output standalone_rw %s --prefix=int -f' % pjoin(path_me, data['paths'][1])
+        mgcmd.exec_cmd(commandline, precmd=True) 
+        
+        #put back golem to original value
+        mgcmd.options['golem'] = old_options['golem']
+        mgcmd.options['pjfry'] = old_options['pjfry']
+        # update make_opts
+        m_opts = {}
+        if mgcmd.options['lhapdf']:
+            #lhapdfversion = subprocess.Popen([mgcmd.options['lhapdf'], '--version'], 
+            #        stdout = subprocess.PIPE).stdout.read().strip()[0]
+            m_opts['lhapdf'] = True
+            m_opts['f2pymode'] = True
+            m_opts['lhapdfversion'] = 5 # 6 always fail on my computer since 5 is compatible but slower always use 5
+            m_opts['llhapdf'] = self.mother.get_lhapdf_libdir()                       
+        else:
+            raise Exception, "NLO reweighting requires LHAPDF to work correctly"
+
+        path = pjoin(path_me,data['paths'][1], 'Source', 'make_opts')             
+        common_run_interface.CommonRunCmd.update_make_opts_full(path, m_opts)      
+        logger.info('Done %.4g' % (time.time()-start))
+
+
+        # Download LHAPDF SET
+        common_run_interface.CommonRunCmd.install_lhapdf_pdfset_static(\
+            mgcmd.options['lhapdf'], None, self.banner.run_card.get_lhapdf_id())
+        
+        # now store the id information             
+        if False:
+            # keep it for debugging purposes
+            matrix_elements = mgcmd._curr_matrix_elements.get_matrix_elements()            
+            for me in matrix_elements:
+                for proc in me.get('processes'):
+                    initial = []    #filled in the next line
+                    final = [l.get('id') for l in proc.get('legs')\
+                          if l.get('state') or initial.append(l.get('id'))]
+                    order = (initial, final)
+                    tag = proc.get_initial_final_ids()
+                    decay_finals = proc.get_final_ids_after_decay()
+    
+                    if tag[1] != decay_finals:
+                        order = (initial, list(decay_finals))
+                        decay_finals.sort()
+                        tag = (tag[0], tuple(decay_finals))
+                    Pdir = pjoin(path_me, data['paths'][1], 'SubProcesses', 
+                                      'P%s' % me.get('processes')[0].shell_string())
+                    assert os.path.exists(Pdir), "Pdir %s do not exists" % Pdir                        
+                    if (tag,'V') in data['id2path']:
+                        if not Pdir == data['id2path'][(tag,'V')][1]:
+                            misc.sprint(tag, Pdir, self.id_to_path[(tag,'V')][1])
+                            raise self.InvalidCmd, '2 different process have the same final states. This module can not handle such situation'
+                        else:
+                            continue
+                    # build the helicity dictionary
+                    hel_nb = 0
+                    hel_dict = {9:0} # unknown helicity -> use full ME
+                    for helicities in me.get_helicity_matrix():
+                        hel_nb +=1 #fortran starts at 1
+                        hel_dict[tuple(helicities)] = hel_nb
+    
+                    data['id2path'][(tag,'V')] = [order, Pdir, hel_dict]
+
+
+    @misc.mute_logger()
     def create_standalone_directory(self, second=False):
         """generate the various directory for the weight evaluation"""
                 
@@ -1296,8 +1431,8 @@
             data['processes'] += [' '.join(line.split()[2:]) for line in self.banner.proc_card
                       if re.search('^\s*add\s+process', line)]  
             #object_collector
-            self.id_to_path = {}
-            data['id2path'] = self.id_to_path
+            #self.id_to_path = {}
+            #data['id2path'] = self.id_to_path
         else:
             data['paths'] = ['rw_me_second', 'rw_mevirt_second']
             # model
@@ -1322,14 +1457,15 @@
                                       for line in self.banner.proc_card
                                       if re.search('^\s*add\s+process', line)]
             #object_collector
-            self.id_to_path_second = {}   
-            data['id2path'] = self.id_to_path_second 
+            #self.id_to_path_second = {}   
+            #data['id2path'] = self.id_to_path_second 
         
         # 0. clean previous run ------------------------------------------------
         if not self.rwgt_dir:
             path_me = self.me_dir
         else:
             path_me = self.rwgt_dir
+        data['path'] = path_me
         try:
             shutil.rmtree(pjoin(path_me,data['paths'][0]))
         except Exception: 
@@ -1374,229 +1510,45 @@
             for name, content in self.banner.get('proc_card', 'multiparticles'):
                 mgcmd.exec_cmd("define %s = %s" % (name, content))
         
-        # 2. compute the production matrix element -----------------------------
-        has_nlo = False  
-        mgcmd.exec_cmd("set group_subprocesses False")
-
-        if not second:
-            logger.info('generating the square matrix element for reweighting')
+        if  second and 'tree_path' in self.dedicated_path:
+            files.ln(self.dedicated_path['tree_path'], path_me,name=data['paths'][0])
+            if 'virtual_path' in self.dedicated_path:
+                has_nlo=True
+            else:
+                has_nlo=False
         else:
-            logger.info('generating the square matrix element for reweighting (second model and/or processes)')
-        start = time.time()
-        commandline=''
-        for i,proc in enumerate(data['processes']):
-            if '[' not in proc:
-                commandline += "add process %s ;" % proc
-            else:
-                has_nlo = True
-                if self.banner.get('run_card','ickkw') == 3:
-                    if len(proc) == min([len(p.strip()) for p in data['processes']]):
-                        commandline += self.get_LO_definition_from_NLO(proc,self.model)
-                    else:
-                        commandline += self.get_LO_definition_from_NLO(proc,self.model, real_only=True)
-                else:
-                    commandline += self.get_LO_definition_from_NLO(proc,self.model)
-        
-        commandline = commandline.replace('add process', 'generate',1)
-        logger.info(commandline)
-        try:
-            mgcmd.exec_cmd(commandline, precmd=True, errorhandling=False)
-        except diagram_generation.NoDiagramException:
-            commandline=''
-            for proc in data['processes']:
-                if '[' not in proc:
-                    raise
-                # pass to virtsq=
-                base, post = proc.split('[',1)
-                nlo_order, post = post.split(']',1)
-                if '=' not in nlo_order:
-                    nlo_order = 'virt=%s' % nlo_order
-                elif 'noborn' in nlo_order:
-                    nlo_order = nlo_order.replace('noborn', 'virt')
-                commandline += "add process %s [%s] %s;" % (base,nlo_order,post)
-            commandline = commandline.replace('add process', 'generate',1)
-            logger.info("RETRY with %s", commandline)
-            mgcmd.exec_cmd(commandline, precmd=True)
-            has_nlo = False
-        except Exception, error:
-            raise
-        
-        commandline = 'output standalone_rw %s' % pjoin(path_me,data['paths'][0])
-        mgcmd.exec_cmd(commandline, precmd=True)        
-        logger.info('Done %.4g' % (time.time()-start))
-        self.has_standalone_dir = True
-        
-
-        # 3. Store id to directory information ---------------------------------
-        matrix_elements = mgcmd._curr_matrix_elements.get_matrix_elements()
-        
-        to_check = [] # list of tag that do not have a Pdir at creation time.
-        for me in matrix_elements:
-            for proc in me.get('processes'):
-                initial = []    #filled in the next line
-                final = [l.get('id') for l in proc.get('legs')\
-                      if l.get('state') or initial.append(l.get('id'))]
-                order = (initial, final)
-                tag = proc.get_initial_final_ids()
-                decay_finals = proc.get_final_ids_after_decay()
-
-                if tag[1] != decay_finals:
-                    order = (initial, list(decay_finals))
-                    decay_finals.sort()
-                    tag = (tag[0], tuple(decay_finals))
-                Pdir = pjoin(path_me, data['paths'][0], 'SubProcesses', 
-                                  'P%s' % me.get('processes')[0].shell_string())
-
-                if not os.path.exists(Pdir):
-                    to_check.append(tag)
-                    continue                        
-                if tag in data['id2path']:
-                    if not Pdir == data['id2path'][tag][1]:
-                        misc.sprint(tag, Pdir, data['id2path'][tag][1])
-                        raise self.InvalidCmd, '2 different process have the same final states. This module can not handle such situation'
-                    else:
-                        continue
-                # build the helicity dictionary
-                hel_nb = 0
-                hel_dict = {9:0} # unknown helicity -> use full ME
-                for helicities in me.get_helicity_matrix():
-                    hel_nb +=1 #fortran starts at 1
-                    hel_dict[tuple(helicities)] = hel_nb
-
-                data['id2path'][tag] = [order, Pdir, hel_dict]        
- 
-        for tag in to_check:
-            if tag not in self.id_to_path:
-                logger.warning("no valid path for %s" % (tag,))
-                #raise self.InvalidCmd, "no valid path for %s" % (tag,)
-        
-        # 4. Check MadLoopParam for Loop induced
-        if os.path.exists(pjoin(path_me, data['paths'][0], 'Cards', 'MadLoopParams.dat')):
-            MLCard = banner.MadLoopParam(pjoin(path_me, data['paths'][0], 'Cards', 'MadLoopParams.dat'))
-            MLCard.set('WriteOutFilters', False)
-            MLCard.set('UseLoopFilter', False)
-            MLCard.set("DoubleCheckHelicityFilter", False)
-            MLCard.set("HelicityFilterLevel", 0)
-            MLCard.write(pjoin(path_me, data['paths'][0], 'SubProcesses', 'MadLoopParams.dat'),
-                         pjoin(path_me, data['paths'][0], 'Cards', 'MadLoopParams.dat'), 
-                         commentdefault=False)
-            
-            if self.multicore == 'create':
-                try:
-                    misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][0],'SubProcesses'),
-                             nb_core=self.mother.options['nb_core'])
-                except:
-                    misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][0],'SubProcesses'))                    
-        
-        if os.path.exists(pjoin(path_me, data['paths'][1], 'Cards', 'MadLoopParams.dat')):
-            if self.multicore == 'create':
-                print "compile OLP", data['paths'][1]
-                try:
-                    misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][1],'SubProcesses'),
-                             nb_core=self.mother.options['nb_core'])
-                except:
-                    misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][1],'SubProcesses'))                    
-                
+            has_nlo = self.create_standalone_tree_directory(data, second)
+
         
         # 5. create the virtual for NLO reweighting  ---------------------------
-        if has_nlo and 'NLO' in self.rwgt_mode:
-            # Do not pass here for LO/NLO_tree
-            start = time.time()
-            commandline=''
-            for proc in data['processes']:
-                if '[' not in proc:
-                    pass
-                else:
-                    proc = proc.replace('[', '[ virt=')
-                    commandline += "add process %s ;" % proc
-            # deactivate golem since it creates troubles
-            old_options = dict(mgcmd.options)
-            if mgcmd.options['golem'] or mgcmd.options['pjfry']:
-                logger.info(" When doing NLO reweighting, MG5aMC cannot use the loop reduction algorithms Golem and/or PJFry++")
-            mgcmd.options['golem'] = None            
-            mgcmd.options['pjfry'] = None 
-            commandline = commandline.replace('add process', 'generate',1)
-            logger.info(commandline)
-            mgcmd.exec_cmd(commandline, precmd=True)
-            commandline = 'output standalone_rw %s -f' % pjoin(path_me, data['paths'][1])
-            mgcmd.exec_cmd(commandline, precmd=True) 
-            
-            #put back golem to original value
-            mgcmd.options['golem'] = old_options['golem']
-            mgcmd.options['pjfry'] = old_options['pjfry']
-            # update make_opts
-            m_opts = {}
-            if mgcmd.options['lhapdf']:
-                #lhapdfversion = subprocess.Popen([mgcmd.options['lhapdf'], '--version'], 
-                #        stdout = subprocess.PIPE).stdout.read().strip()[0]
-                m_opts['lhapdf'] = True
-                m_opts['f2pymode'] = True
-                m_opts['lhapdfversion'] = 5 # 6 always fail on my computer since 5 is compatible but slower always use 5
-                m_opts['llhapdf'] = self.mother.get_lhapdf_libdir()                       
-            else:
-                raise Exception, "NLO reweighting requires LHAPDF to work correctly"
- 
-            path = pjoin(path_me,data['paths'][1], 'Source', 'make_opts')             
-            common_run_interface.CommonRunCmd.update_make_opts_full(path, m_opts)      
-            logger.info('Done %.4g' % (time.time()-start))
-
-
-            # Download LHAPDF SET
-            common_run_interface.CommonRunCmd.install_lhapdf_pdfset_static(\
-                mgcmd.options['lhapdf'], None, self.banner.run_card.get_lhapdf_id())
-            
-            # now store the id information             
-            matrix_elements = mgcmd._curr_matrix_elements.get_matrix_elements()            
-            for me in matrix_elements:
-                for proc in me.get('processes'):
-                    initial = []    #filled in the next line
-                    final = [l.get('id') for l in proc.get('legs')\
-                          if l.get('state') or initial.append(l.get('id'))]
-                    order = (initial, final)
-                    tag = proc.get_initial_final_ids()
-                    decay_finals = proc.get_final_ids_after_decay()
-    
-                    if tag[1] != decay_finals:
-                        order = (initial, list(decay_finals))
-                        decay_finals.sort()
-                        tag = (tag[0], tuple(decay_finals))
-                    Pdir = pjoin(path_me, data['paths'][1], 'SubProcesses', 
-                                      'P%s' % me.get('processes')[0].shell_string())
-                    assert os.path.exists(Pdir), "Pdir %s do not exists" % Pdir                        
-                    if (tag,'V') in data['id2path']:
-                        if not Pdir == data['id2path'][(tag,'V')][1]:
-                            misc.sprint(tag, Pdir, self.id_to_path[(tag,'V')][1])
-                            raise self.InvalidCmd, '2 different process have the same final states. This module can not handle such situation'
-                        else:
-                            continue
-                    # build the helicity dictionary
-                    hel_nb = 0
-                    hel_dict = {9:0} # unknown helicity -> use full ME
-                    for helicities in me.get_helicity_matrix():
-                        hel_nb +=1 #fortran starts at 1
-                        hel_dict[tuple(helicities)] = hel_nb
-    
-                    data['id2path'][(tag,'V')] = [order, Pdir, hel_dict] 
-                
-            #compile the module to combine the weight
-            misc.compile(cwd=pjoin(path_me, data['paths'][1], 'Source'))
-            #link it 
-            if path_me not in sys.path:
-                sys.path.insert(0, os.path.realpath(path_me))
-            with misc.chdir(pjoin(path_me)):
-                mymod = __import__('%s.Source.rwgt2py' % data['paths'][1], globals(), locals(), [],-1)
-                mymod =  mymod.Source.rwgt2py
-                with misc.stdchannel_redirected(sys.stdout, os.devnull):
-                    mymod.initialise([self.banner.run_card['lpp1'], 
-                                  self.banner.run_card['lpp2']],
-                                 self.banner.run_card.get_lhapdf_id())
-                self.combine_wgt = mymod.get_wgt
+        if second and 'virtual_path' in self.dedicated_path:
+            files.ln(self.dedicated_path['virtual_path'], path_me, name=data['paths'][1])
+        elif has_nlo and 'NLO' in self.rwgt_mode:
+            self.create_standalone_virt_directory(data, second)
+            
+            if not second:
+                #compile the module to combine the weight
+                misc.compile(cwd=pjoin(path_me, data['paths'][1], 'Source'))
+                #link it 
+                if path_me not in sys.path:
+                    sys.path.insert(0, os.path.realpath(path_me))
+                with misc.chdir(pjoin(path_me)):
+                    mymod = __import__('%s.Source.rwgt2py' % data['paths'][1], globals(), locals(), [],-1)
+                    mymod =  mymod.Source.rwgt2py
+                    with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                        mymod.initialise([self.banner.run_card['lpp1'], 
+                                      self.banner.run_card['lpp2']],
+                                     self.banner.run_card.get_lhapdf_id())
+                    self.combine_wgt = mymod.get_wgt
 
             if self.multicore == 'create':
                 print "compile OLP", data['paths'][1]
-                misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][1],'SubProcesses'),
+                try:
+                    misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][1],'SubProcesses'),
                              nb_core=self.mother.options['nb_core'])
-                
+                except:
+                    misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][1],'SubProcesses'),
+                             nb_core=1)
         elif has_nlo and not second and self.rwgt_mode == ['NLO_tree']:
             # We do not have any virtual reweighting to do but we still have to
             #combine the weights.
@@ -1610,7 +1562,7 @@
             commandline = commandline.replace('add process', 'generate',1)
             logger.info(commandline)
             mgcmd.exec_cmd(commandline, precmd=True)
-            commandline = 'output standalone_rw %s -f' % pjoin(path_me, data['paths'][1])
+            commandline = 'output standalone_rw %s --prefix=int -f' % pjoin(path_me, data['paths'][1])
             mgcmd.exec_cmd(commandline, precmd=True)    
             #put back golem to original value
             mgcmd.options['golem'] = old_options['golem']
@@ -1651,7 +1603,7 @@
             
              
         # 6. If we need a new model/process-------------------------------------
-        if (self.second_model or self.second_process) and not second:
+        if (self.second_model or self.second_process or self.dedicated_path) and not second :
             self.create_standalone_directory(second=True)    
 
         if not second:
@@ -1659,7 +1611,123 @@
             
         
         
-
+    def compile(self):
+        """compile the code"""
+        
+        if self.multicore=='wait':
+            return
+        
+        if not self.rwgt_dir:
+            path_me = self.me_dir
+        else:
+            path_me = self.rwgt_dir
+        for onedir in self.rwgt_dir_possibility:
+            if not os.path.isdir(pjoin(path_me,onedir)):
+                continue
+            pdir = pjoin(path_me, onedir, 'SubProcesses')
+            nb_core = self.mother.options['nb_core'] if self.mother.options['run_mode'] !=0 else 1
+            os.environ['MENUM'] = '2'
+            misc.compile(['allmatrix2py.so'], cwd=pdir, nb_core=nb_core)
+            if not (self.second_model or self.second_process or self.dedicated_path):
+                os.environ['MENUM'] = '3'
+                misc.compile(['allmatrix3py.so'], cwd=pdir, nb_core=nb_core)
+
+    def load_module(self, metag=1):
+        """load the various module and load the associate information"""
+        
+        
+        if not self.rwgt_dir:
+            path_me = self.me_dir
+        else:
+            path_me = self.rwgt_dir        
+        
+        self.id_to_path = {}
+        self.id_to_path_second = {}
+        for onedir in self.rwgt_dir_possibility:
+            if not os.path.exists(pjoin(path_me,onedir)):
+                continue 
+            pdir = pjoin(path_me, onedir, 'SubProcesses')
+            for tag in [2*metag,2*metag+1]:
+                with misc.TMP_variable(sys, 'path', [pjoin(path_me)]+sys.path):
+                    mymod = __import__('%s.SubProcesses.allmatrix%spy' % (onedir, tag), globals(), locals(), [],-1)
+                    reload(mymod)
+                    S = mymod.SubProcesses
+                    mymod = getattr(S, 'allmatrix%spy' % tag)
+                
+                # Param card not available -> no initialisation
+                self.f2pylib[(onedir,tag)] = mymod
+                if hasattr(mymod, 'set_madloop_path'):
+                    mymod.set_madloop_path(pjoin(path_me,onedir,'SubProcesses','MadLoop5_resources'))
+                if (self.second_model or self.second_process or self.dedicated_path):
+                    break
+
+            data = self.id_to_path
+            if '_second' in onedir:
+                data = self.id_to_path_second
+                
+            # get all the information
+            all_pdgs = mymod.get_pdg_order()
+            all_pdgs = [[pdg for pdg in pdgs if pdg!=0] for pdgs in  mymod.get_pdg_order()]
+            all_prefix = [''.join(j).strip().lower() for j in mymod.get_prefix()]
+            prefix_set = set(all_prefix)
+
+
+            hel_dict={}
+            for prefix in prefix_set:
+                if hasattr(mymod,'%sprocess_nhel' % prefix):
+                    nhel = getattr(mymod, '%sprocess_nhel' % prefix).nhel    
+                    hel_dict[prefix] = {}
+                    for i, onehel in enumerate(zip(*nhel)):
+                        hel_dict[prefix][tuple(onehel)] = i+1
+                elif hasattr(mymod, 'set_madloop_path') and \
+                     os.path.exists(pjoin(path_me,onedir,'SubProcesses','MadLoop5_resources', '%sHelConfigs.dat' % prefix.upper())):
+                    hel_dict[prefix] = {}
+                    for i,line in enumerate(open(pjoin(path_me,onedir,'SubProcesses','MadLoop5_resources', '%sHelConfigs.dat' % prefix.upper()))):
+                        onehel = [int(h) for h in line.split()]
+                        hel_dict[prefix][tuple(onehel)] = i+1
+                else:
+                    misc.sprint(pjoin(path_me,onedir,'SubProcesses','MadLoop5_resources', '%sHelConfigs.dat' % prefix.upper() ))
+                    misc.sprint(os.path.exists(pjoin(path_me,onedir,'SubProcesses','MadLoop5_resources', '%sHelConfigs.dat' % prefix.upper())))
+                    continue
+
+            for i,pdg in enumerate(all_pdgs):
+                if self.is_decay:
+                    incoming = pdg[0]
+                    outgoing = pdg[1:]
+                else:
+                    incoming = pdg[0:2]
+                    outgoing = pdg[2:]
+                order = (list(incoming), list(outgoing))
+                incoming.sort()
+                outgoing.sort()
+                tag = (tuple(incoming), tuple(outgoing))
+                if 'virt' in onedir:
+                    tag = (tag, 'V')
+                prefix = all_prefix[i]
+                hel = hel_dict[prefix]
+                if tag in data:
+                    oldpdg = data[tag][0][0]+data[tag][0][1]
+                    if all_prefix[all_pdgs.index(pdg)] == all_prefix[all_pdgs.index(oldpdg)]:
+                        for i in range(len(pdg)):
+                            if pdg[i] == oldpdg[i]:
+                                continue
+                            if not self.model:
+                                continue
+                            if self.model.get_mass(int(pdg[i])) == self.model.get_mass(int(oldpdg[i])):
+                                continue
+                            misc.sprint(tag, onedir)
+                            misc.sprint(data[tag][:-1])
+                            misc.sprint(order, pdir,)
+                            raise Exception                                
+                    else: 
+                        misc.sprint(tag, onedir)
+                        misc.sprint(data[tag][:-1])
+                        misc.sprint(order, pdir,)
+                        raise Exception
+                
+                data[tag] = order, pdir, hel
+             
+             
     def load_model(self, name, use_mg_default, complex_mass=False):
         """load the model"""
         

=== modified file 'madgraph/iolibs/export_cpp.py'
--- madgraph/iolibs/export_cpp.py	2016-08-26 23:45:51 +0000
+++ madgraph/iolibs/export_cpp.py	2017-07-29 19:33:46 +0000
@@ -2512,12 +2512,13 @@
          
         return OneProcessExporterPythia8.read_template_file(*args, **opts)
 
-def ExportCPPFactory(cmd, group_subprocesses=False):
+def ExportCPPFactory(cmd, group_subprocesses=False, cmd_options={}):
     """ Determine which Export class is required. cmd is the command 
         interface containing all potential usefull information.
     """
 
-    opt = cmd.options
+    opt = dict(cmd.options)
+    opt['output_options'] = cmd_options
     cformat = cmd._export_format
     
     if cformat == 'pythia8':

=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2017-03-03 10:07:55 +0000
+++ madgraph/iolibs/export_v4.py	2017-07-29 19:33:46 +0000
@@ -96,6 +96,7 @@
     
     
     def __init__(self, dir_path = "", opt=None):
+        # cmd_options is a dictionary with all the optional argurment passed at output time 
         return
 
     def copy_template(self, model):
@@ -138,7 +139,8 @@
 
     default_opt = {'clean': False, 'complex_mass':False,
                         'export_format':'madevent', 'mp': False,
-                        'v5_model': True
+                        'v5_model': True,
+                        'output_options':{}
                         }
     grouped_mode = False
 
@@ -152,6 +154,8 @@
         if opt:
             self.opt.update(opt)
         
+        self.cmd_options = self.opt['output_options']
+        
         #place holder to pass information to the run_interface
         self.proc_characteristic = banner_mod.ProcCharacteristic()
         
@@ -1848,7 +1852,7 @@
 
     matrix_template = "matrix_standalone_v4.inc"
 
-    def __init__(self, *args, **opts):
+    def __init__(self, *args,**opts):
         """add the format information compare to standard init"""
         
         if 'format' in opts:
@@ -1856,6 +1860,8 @@
             del opts['format']
         else:
             self.format = 'standalone'
+        
+        self.prefix_info = {}
         ProcessExporterFortran.__init__(self, *args, **opts)
 
     def copy_template(self, model):
@@ -1894,7 +1900,7 @@
         
         # Add file in SubProcesses
         shutil.copy(pjoin(self.mgme_dir, 'madgraph', 'iolibs', 'template_files', 'makefile_sa_f_sp'), 
-                    pjoin(self.dir_path, 'SubProcesses', 'makefile'))
+                    pjoin(self.dir_path, 'SubProcesses', 'makefileP'))
         
         if self.format == 'standalone':
             shutil.copy(pjoin(self.mgme_dir, 'madgraph', 'iolibs', 'template_files', 'check_sa.f'), 
@@ -1982,24 +1988,141 @@
                            pjoin(self.dir_path, 'Source', 'PDF'))
             self.write_pdf_opendata()
             
-        # create a single makefile to compile all the subprocesses
-        text = '''\n# For python linking (require f2py part of numpy)\nifeq ($(origin MENUM),undefined)\n  MENUM=2\nendif\n''' 
-        deppython = ''
-        for Pdir in os.listdir(pjoin(self.dir_path,'SubProcesses')):
-            if os.path.isdir(pjoin(self.dir_path, 'SubProcesses', Pdir)):
-                text += '%(0)s/matrix$(MENUM)py.so:\n\tcd %(0)s;make matrix$(MENUM)py.so\n'% {'0': Pdir}
-                deppython += ' %(0)s/matrix$(MENUM)py.so ' % {'0': Pdir}
-        
-        text+='all: %s\n\techo \'done\'' % deppython
-            
-        ff = open(pjoin(self.dir_path, 'SubProcesses', 'makefile'),'a')
-        ff.write(text)
-        ff.close()
-        
-        
-        
-        
-            
+        if self.prefix_info:
+            self.write_f2py_splitter()
+            self.write_f2py_makefile()
+            self.write_f2py_check_sa(matrix_elements,
+                            pjoin(self.dir_path,'SubProcesses','check_sa.py'))
+        else:
+            # create a single makefile to compile all the subprocesses
+            text = '''\n# For python linking (require f2py part of numpy)\nifeq ($(origin MENUM),undefined)\n  MENUM=2\nendif\n''' 
+            deppython = ''
+            for Pdir in os.listdir(pjoin(self.dir_path,'SubProcesses')):
+                if os.path.isdir(pjoin(self.dir_path, 'SubProcesses', Pdir)):
+                    text += '%(0)s/matrix$(MENUM)py.so:\n\tcd %(0)s;make matrix$(MENUM)py.so\n'% {'0': Pdir}
+                    deppython += ' %(0)s/matrix$(MENUM)py.so ' % {'0': Pdir}
+            text+='all: %s\n\techo \'done\'' % deppython
+            
+            ff = open(pjoin(self.dir_path, 'SubProcesses', 'makefile'),'a')
+            ff.write(text)
+            ff.close()
+                    
+    def write_f2py_splitter(self):
+        """write a function to call the correct matrix element"""
+        
+        template = """
+%(python_information)s
+  subroutine smatrixhel(pdgs, npdg, p, ALPHAS, SCALE2, nhel, ANS)
+  IMPLICIT NONE
+
+CF2PY double precision, intent(in), dimension(0:3,npdg) :: p
+CF2PY integer, intent(in), dimension(npdg) :: pdgs
+CF2PY integer, intent(in) :: npdg
+CF2PY double precision, intent(out) :: ANS
+CF2PY double precision, intent(in) :: ALPHAS
+CF2PY double precision, intent(in) :: SCALE2
+  integer pdgs(*)
+  integer npdg, nhel
+  double precision p(*)
+  double precision ANS, ALPHAS, PI,SCALE2
+  include 'coupl.inc'
+  
+  PI = 3.141592653589793D0
+  G = 2* DSQRT(ALPHAS*PI)
+  CALL UPDATE_AS_PARAM()
+  if (scale2.ne.0d0) stop 1
+
+%(smatrixhel)s
+
+      return
+      end
+  
+      SUBROUTINE INITIALISE(PATH)
+C     ROUTINE FOR F2PY to read the benchmark point.
+      IMPLICIT NONE
+      CHARACTER*512 PATH
+CF2PY INTENT(IN) :: PATH
+      CALL SETPARA(PATH)  !first call to setup the paramaters
+      RETURN
+      END
+
+    subroutine get_pdg_order(PDG)
+  IMPLICIT NONE
+CF2PY INTEGER, intent(out) :: PDG(%(nb_me)i,%(maxpart)i)  
+  INTEGER PDG(%(nb_me)i,%(maxpart)i), PDGS(%(nb_me)i,%(maxpart)i)
+  DATA PDGS/ %(pdgs)s /
+  PDG = PDGS
+  RETURN
+  END 
+
+    subroutine get_prefix(PREFIX)
+  IMPLICIT NONE
+CF2PY CHARACTER*20, intent(out) :: PREFIX(%(nb_me)i)
+  character*20 PREFIX(%(nb_me)i),PREF(%(nb_me)i)
+  DATA PREF / '%(prefix)s'/
+  PREFIX = PREF
+  RETURN
+  END 
+ 
+  
+        """
+         
+        allids = self.prefix_info.keys()
+        allprefix = [self.prefix_info[key][0] for key in allids]
+        min_nexternal = min([len(ids) for ids in allids])
+        max_nexternal = max([len(ids) for ids in allids])
+
+        info = []
+        for key, (prefix, tag) in self.prefix_info.items():
+            info.append('#PY %s : %s # %s' % (tag, key, prefix))
+            
+
+        text = []
+        for n_ext in range(min_nexternal, max_nexternal+1):
+            current = [ids for ids in allids if len(ids)==n_ext]
+            if not current:
+                continue
+            if min_nexternal != max_nexternal:
+                if n_ext == min_nexternal:
+                    text.append('       if (npdg.eq.%i)then' % n_ext)
+                else:
+                    text.append('       else if (npdg.eq.%i)then' % n_ext)
+            for ii,pdgs in enumerate(current):
+                condition = '.and.'.join(['%i.eq.pdgs(%i)' %(pdg, i+1) for i, pdg in enumerate(pdgs)])
+                if ii==0:
+                    text.append( ' if(%s) then ! %i' % (condition, i))
+                else:
+                    text.append( ' else if(%s) then ! %i' % (condition,i))
+                text.append(' call %ssmatrixhel(p, nhel, ans)' % self.prefix_info[pdgs][0])
+            text.append(' endif')
+        #close the function
+        if min_nexternal != max_nexternal:
+            text.append('endif')
+
+        formatting = {'python_information':'\n'.join(info), 
+                          'smatrixhel': '\n'.join(text),
+                          'maxpart': max_nexternal,
+                          'nb_me': len(allids),
+                          'pdgs': ','.join(str(pdg[i]) if i<len(pdg) else '0' 
+                                             for i in range(max_nexternal) for pdg in allids),
+                          'prefix':'\',\''.join(allprefix)
+                          }
+        formatting['lenprefix'] = len(formatting['prefix'])
+        text = template % formatting
+        fsock = writers.FortranWriter(pjoin(self.dir_path, 'SubProcesses', 'all_matrix.f'),'w')
+        fsock.writelines(text)
+        fsock.close()
+            
+    def write_f2py_check_sa(self, matrix_element, writer):
+        """ Write the general check_sa.py in SubProcesses that calls all processes successively."""
+        # To be implemented. It is just an example file, i.e. not crucial.
+        return
+    
+    def write_f2py_makefile(self):
+        """ """
+        # Add file in SubProcesses
+        shutil.copy(pjoin(self.mgme_dir, 'madgraph', 'iolibs', 'template_files', 'makefile_sa_f2py'), 
+                    pjoin(self.dir_path, 'SubProcesses', 'makefile'))
 
     def create_MA5_cards(self,*args,**opts):
         """ Overload the function of the mother so as to bypass this in StandAlone."""
@@ -2027,13 +2150,8 @@
         if self.opt['sa_symmetry']:
             # avoid symmetric output
             for i,proc in enumerate(matrix_element.get('processes')):
-                        
-                initial = []    #filled in the next line
-                final = [l.get('id') for l in proc.get('legs')\
-                      if l.get('state') or initial.append(l.get('id'))]
-                decay_finals = proc.get_final_ids_after_decay()
-                decay_finals.sort()
-                tag = (tuple(initial), tuple(decay_finals))
+                   
+                tag = proc.get_tag()     
                 legs = proc.get('legs')[:]
                 leg0 = proc.get('legs')[0]
                 leg1 = proc.get('legs')[1]
@@ -2074,10 +2192,24 @@
             filename = pjoin(dirpath, 'matrix_prod.f')
         else:
             filename = pjoin(dirpath, 'matrix.f')
+            
+        proc_prefix = ''
+        if 'prefix' in self.cmd_options:
+            if self.cmd_options['prefix'] == 'int':
+                proc_prefix = 'M%s_' % number
+            elif self.cmd_options['prefix'] == 'proc':
+                proc_prefix = matrix_element.get('processes')[0].shell_string().split('_',1)[1]
+            else:
+                raise Exception, '--prefix options supports only \'int\' and \'proc\''
+            for proc in matrix_element.get('processes'):
+                ids = [l.get('id') for l in proc.get('legs_with_decays')]
+                self.prefix_info[tuple(ids)] = [proc_prefix, proc.get_tag()] 
+                
         calls = self.write_matrix_element_v4(
             writers.FortranWriter(filename),
             matrix_element,
-            fortran_model)
+            fortran_model,
+            proc_prefix=proc_prefix)
 
         if self.opt['export_format'] == 'standalone_msP':
             filename =  pjoin(dirpath,'configs_production.inc')
@@ -2124,11 +2256,18 @@
                      matrix_element.get('processes')[0].nice_string())
         plot.draw()
 
-        linkfiles = ['check_sa.f', 'coupl.inc', 'makefile']
+        linkfiles = ['check_sa.f', 'coupl.inc']
+
+        if proc_prefix and os.path.exists(pjoin(dirpath, '..', 'check_sa.f')):
+            text = open(pjoin(dirpath, '..', 'check_sa.f')).read()
+            new_text, n  = re.subn('smatrix', '%ssmatrix' % proc_prefix, text, flags=re.I)
+            with open(pjoin(dirpath, 'check_sa.f'),'w') as f:
+                f.write(new_text)
+            linkfiles.pop(0)
 
         for file in linkfiles:
             ln('../%s' % file, cwd=dirpath)
-
+        ln('../makefileP', name='makefile', cwd=dirpath)
         # Return to original PWD
         #os.chdir(cwd)
 
@@ -2176,7 +2315,6 @@
             self.opt['sa_symmetry']=False
 
 
-
         # The proc_id is for MadEvent grouping which is never used in SA.
         replace_dict = {'global_variable':'', 'amp2_lines':'',
                                        'proc_prefix':proc_prefix, 'proc_id':''}
@@ -6445,14 +6583,15 @@
                                       rule_card_path=rule_card, 
                                       mssm_convert=True)
         
-def ExportV4Factory(cmd, noclean, output_type='default', group_subprocesses=True):
+def ExportV4Factory(cmd, noclean, output_type='default', group_subprocesses=True, cmd_options={}):
     """ Determine which Export_v4 class is required. cmd is the command 
         interface containing all potential usefull information.
         The output_type argument specifies from which context the output
         is called. It is 'madloop' for MadLoop5, 'amcatnlo' for FKS5 output
         and 'default' for tree-level outputs."""
 
-    opt = cmd.options
+    opt = dict(cmd.options)
+    opt['output_options'] = cmd_options
 
     # ==========================================================================
     # First check whether Ninja must be installed.
@@ -6501,10 +6640,10 @@
       'SubProc_prefix':'P',
       'compute_color_flows':cmd.options['loop_color_flows'],
       'mode': 'reweight' if cmd._export_format == "standalone_rw" else '',
-      'cluster_local_path': cmd.options['cluster_local_path']
+      'cluster_local_path': cmd.options['cluster_local_path'],
+      'output_options': cmd_options
       }
 
-
     if output_type.startswith('madloop'):        
         import madgraph.loop.loop_exporters as loop_exporters
         if os.path.isdir(os.path.join(cmd._mgme_dir, 'Template/loop_material')):

=== added symlink 'madgraph/iolibs/template_files/loop/check_sa_all.py.inc'
=== target is u'../loop_optimized/check_sa_all.py.inc'
=== modified file 'madgraph/iolibs/template_files/loop/improve_ps.inc'
--- madgraph/iolibs/template_files/loop/improve_ps.inc	2014-11-17 15:32:46 +0000
+++ madgraph/iolibs/template_files/loop/improve_ps.inc	2017-07-29 19:33:46 +0000
@@ -733,7 +733,14 @@
 C     Apply the rescaling
       DO I=1,NEXTERNAL
 	    DO J=1,3
-		  NEWP(J,I)=NEWP(J,I)*XSCALE
+C         Consider scaling by x**2 for the first particle so that
+C         the algorithm for numerically solving for XSCALE has a non-vanishing
+C         derivative in the case that all particle are massless.
+		  IF (I.eq.1) THEN
+		    NEWP(J,I)=NEWP(J,I)*XSCALE**2
+	      ELSE
+		    NEWP(J,I)=NEWP(J,I)*XSCALE		  
+		  ENDIF
 		ENDDO
 	  ENDDO
 
@@ -912,20 +919,32 @@
 	  RES=ZERO
 	  BUFF=ZERO
 
+C     Consider scaling by x**2 for the first particle so that
+C     the algorithm for numerically solving for XSCALE has a non-vanishing
+C     derivative in the case that all particle are massless.
+
 	  DO I=1,NEXTERNAL
 	    IF (I.LE.NINITIAL) THEN
 		  FACTOR=-ONE
 		ELSE
 		  FACTOR=ONE
 		ENDIF
-		BUFF=MASSES(I)**2+PVECSQ(I)*X**2
+		IF (I.eq.1) THEN
+		  BUFF=MASSES(I)**2+PVECSQ(I)*X**4
+		ELSE
+		  BUFF=MASSES(I)**2+PVECSQ(I)*X**2
+		ENDIF
         IF (BUFF.LT.ZERO) THEN
 		  RES=ZERO
 		  ERROR = 1
 		  GOTO 800
 		ENDIF
 		IF (DERIVATIVE) THEN
-		  RES=RES + FACTOR*((X*PVECSQ(I))/SQRT(BUFF))
+		  IF (I.eq.1) THEN
+		    RES=RES + FACTOR*((2*X*PVECSQ(I))/SQRT(BUFF))
+		  ELSE
+		    RES=RES + FACTOR*((X*PVECSQ(I))/SQRT(BUFF))
+		  ENDIF
 		ELSE
 		  RES=RES + FACTOR*SQRT(BUFF)		  
 		ENDIF

=== modified file 'madgraph/iolibs/template_files/loop/loop_matrix_standalone.inc'
--- madgraph/iolibs/template_files/loop/loop_matrix_standalone.inc	2016-03-22 20:44:25 +0000
+++ madgraph/iolibs/template_files/loop/loop_matrix_standalone.inc	2017-07-29 19:33:46 +0000
@@ -1001,6 +1001,15 @@
       ENDDO
       ACC  = ACC / ( ABS(AVG) / 3.0d0)
 
+c     If NaN are present in the evaluation, automatically set the accuracy to 1.0d99.
+      DO I=1,3
+        DO J=1,MAXSTABILITYLENGTH
+          IF (isnan(FULLLIST(I,J))) then
+            ACC = 1.0d99
+          ENDIF
+        ENDDO
+      ENDDO
+
       end
 
       SUBROUTINE %(proc_prefix)sSET_N_EVALS(N_DP_EVALS,N_QP_EVALS)

=== modified file 'madgraph/iolibs/template_files/loop_optimized/check_py.f.inc'
--- madgraph/iolibs/template_files/loop_optimized/check_py.f.inc	2016-03-03 22:14:52 +0000
+++ madgraph/iolibs/template_files/loop_optimized/check_py.f.inc	2017-07-29 19:33:46 +0000
@@ -1,4 +1,4 @@
-      Subroutine Initialise(path)
+      Subroutine %(prefix_routine)sInitialise(path)
 
       CHARACTER(128) path
 CF2PY intent(in)::path
@@ -12,7 +12,7 @@
       return 
       end
 
-      SUBROUTINE GET_ME(P, ALPHAS, SCALE2, NHEL , ANS,RETURNCODE)
+      SUBROUTINE %(prefix_routine)sGET_ME(P, ALPHAS, SCALE2, NHEL , ANS,RETURNCODE)
       IMPLICIT NONE
 C     
 C     CONSTANTS  

=== modified file 'madgraph/iolibs/template_files/loop_optimized/check_sa.py.inc'
--- madgraph/iolibs/template_files/loop_optimized/check_sa.py.inc	2016-03-04 17:56:39 +0000
+++ madgraph/iolibs/template_files/loop_optimized/check_sa.py.inc	2017-07-29 19:33:46 +0000
@@ -44,7 +44,7 @@
 # print help(matrix2py)
 
 # Read the model parameters
-matrix2py.initialise(os.path.abspath(pjoin(root_path,os.pardir,os.pardir,'Cards','param_card.dat')))
+matrix2py.{2}initialise(os.path.abspath(pjoin(root_path,os.pardir,os.pardir,'Cards','param_card.dat')))
 
 
 def invert_momenta(p):
@@ -86,7 +86,7 @@
 # Choice of renormalization scale
 renormalization_scale = 91.188
 
-finite_loop_me, return_code = matrix2py.get_me(P, alphas, renormalization_scale, nhel)
+finite_loop_me, return_code = matrix2py.{2}get_me(P, alphas, renormalization_scale, nhel)
 
 
 print '='*112

=== added file 'madgraph/iolibs/template_files/loop_optimized/check_sa_all.py.inc'
--- madgraph/iolibs/template_files/loop_optimized/check_sa_all.py.inc	1970-01-01 00:00:00 +0000
+++ madgraph/iolibs/template_files/loop_optimized/check_sa_all.py.inc	2017-07-29 19:33:46 +0000
@@ -0,0 +1,139 @@
+#! /usr/bin/env python
+
+# This is an example of how to run MadLoop from Python using the f2py compilation of the wrapper file 'f2py_wrapper.f'.
+
+import os
+import sys
+import subprocess
+
+pjoin = os.path.join
+root_path = os.path.dirname(os.path.realpath( __file__ ))
+sys.path.insert(0, root_path)
+
+try:
+    import allmatrix2py
+except:
+    if os.path.isfile(pjoin(root_path,'makefile')) and \
+       os.path.isfile(pjoin(root_path,'all_matrix.f')) and \
+       not os.path.isfile(pjoin(root_path,'allmatrix2py.so')):
+        print "Trying to automatically generate the python module 'allmatrix2py.so' with f2py..."
+        p = subprocess.Popen(['make','allmatrix2py.so'], stdout=subprocess.PIPE, 
+                             stderr=subprocess.PIPE, cwd=root_path)
+        (out, err) = p.communicate()
+        if p.returncode or not os.path.isfile(pjoin(root_path,'allmatrix2py.so')):
+            print "ERROR: Failed to produce 'allmatrix2py.so' with 'make allmatrix2py.so' in '%s'. The error was:\n%s"%(root_path,err)
+            sys.exit(0)
+        try:
+            import allmatrix2py
+        except:
+            print "ERROR: Could not load the f2py module 'allmatrix2py.so'. The following error occurred:\n",sys.exc_info()[0]
+            sys.exit(0)
+    else:
+        if os.path.exists(pjoin(root_path,'allmatrix2py.so')):
+            print "ERROR: Could not load the f2py module 'matrix2py.so'. The following error occurred:\n",sys.exc_info()[0]
+            sys.exit(0)
+        else:
+            print "ERROR: Could not find the 'matrix2py.so' f2py module. Please generate it by running:\n"+\
+                  "  > make matrix2py.so\n"+\
+                  "in the <PROC_OUTPUT>/SubProcesses/P<chosen_proc> directory."
+            sys.exit(0)
+
+# Now we can use this MadLoop python module.
+
+# This is a handy way of looking at what is available in the imported f2py module
+# print help(matrix2py)
+
+# Read the model parameters
+allmatrix2py.initialise(os.path.abspath(pjoin(root_path,os.pardir,'Cards','param_card.dat')))
+# Specify where MadLoop5_resources is
+allmatrix2py.set_madloop_path(os.path.abspath(pjoin(root_path,'MadLoop5_resources')))
+# List all processes available in this module
+
+all_procs_available = [[pdg for pdg in proc_pdgs if pdg!=0] for proc_pdgs in allmatrix2py.get_pdg_order()]
+print "Available processes are:"
+for process_pdgs in all_procs_available:
+    print "   --> %s"%str(process_pdgs)
+
+def invert_momenta(p):
+    """ fortran/C-python do not order table in the same order"""
+    return [[p[j][i] for j in range(len(p))] for i in range(len(p[0]))]
+
+# Now choose MadLoop's inputs for this evaluation
+
+# The kinematic configuration in the convention (E, px, py, pz) and with particles ordered as in the process definition.
+# This is the structure of this dictionary.
+p= dict( (tuple(process_pdgs),[[None,]*4]*len(process_pdgs)) for process_pdgs in all_procs_available)
+
+if (p.values()[0])[0][0] is None:
+    if not os.path.isfile(pjoin(root_path,'PS.input')):
+        print "\n\n==================================================================================================================="
+        print "*                                           No kinematics defined!                                                *"
+        print "*                                   -------------------------------------                                         *"
+        print "* Please either define your kinematic configuration directly in check_sa.py or in a file 'PS.input'. Exiting now. *"
+        print " Here is an example of the content of such a file for a process output containing g g > h and g g > h g."
+        print \
+"""
+
+\"\"\"
+Process=[21,21,25]
+6.250000000000000E+001        0.000000000000000E+000        0.000000000000000E+000        6.250000000000000E+001
+6.250000000000000E+001        0.000000000000000E+000        0.000000000000000E+000       -6.250000000000000E+001
+1.250000000000000E+002        0.000000000000000E+000        0.000000000000000E+000        0.000000000000000E+000
+
+Process=[21,21,25,21]
+5.000000000000000E+002        0.000000000000000E+000        0.000000000000000E+000        5.000000000000000E+002
+5.000000000000000E+002        0.000000000000000E+000        0.000000000000000E+000       -5.000000000000000E+002
+5.078125000000001E+002       -1.567603625592635E+002        2.039552324827553E+002        4.196152845517468E+002
+4.921874999999999E+002        1.567603625592635E+002       -2.039552324827552E+002       -4.196152845517467E+002
+\"\"\"
+"""
+        print "===================================================================================================================\n\n"
+        sys.exit(0)
+    try:
+        p               = {}
+        current_process = None
+        current_p       = []
+        for line in open(pjoin(root_path,'PS.input'),'r').readlines():
+            if line.startswith('Process='):
+                if not current_process is None:
+                    p[current_process] = current_p
+                    current_p = []
+                current_process = tuple(eval(line[8:]))
+                continue
+            if len(line.split())>0:
+                current_p.append([float(line.split()[j]) for j in range(4)])
+        p[current_process] = current_p
+    except:
+        print "ERROR: File PS.input is malformed. Error was:\n",sys.exc_info()[0]
+        sys.exit(0)
+
+
+for iProc, (process_pdgs, kinematic_config) in enumerate(p.items()):
+    print '='*112
+    P =invert_momenta(kinematic_config)
+    # Alpha_s value
+    alphas = 0.118
+    # -1 means that the MadLoop averages/sum over helicity configuration. 
+    # If a positive integer is picked it corresponds to the position of the helicity configuration listed in the file
+    # 'MadLoop5_resouces/ML5_<id>_HelConfigs.dat'
+    nhel = -1 # means sum over all helicity
+    # Choice of renormalization scale
+    renormalization_scale = 91.188
+
+    finite_loop_me, return_code = allmatrix2py.smatrixhel(pdgs=list(process_pdgs),
+		                          p=P, alphas=alphas,scales2=renormalization_scale,nhel=nhel)
+
+    print '* %-108s *'%' MadLoop evaluation for the process '
+    print '* %-108s *'%('  %s'%str(process_pdgs))
+    print '* %-108s *'%' and the kinematic configuration:'
+    print '* %-108s *'%((' %-3s'+' %-25s'*4)%('#','E',' p_x',' p_y',' p_z'))
+    for i,k in enumerate(kinematic_config):
+        # The complicated printout below is just so as to align the digits negative numbers with positive ones
+        print '* %-108s *'%((' %-3d%s')%(i,
+          ''.join([' %-25.15e'%e if j==0 or e<0.0 else '  %-24.15e'%e for j,e in enumerate(k)])))
+    print '* %-108s *'%('-'*108)
+    print '* %-108s *'%(' Finite part obtained for the loop matrix element (Madloop return_code=%d )'%return_code)
+    print '* %-108s *'%''
+    print '* %-108s *'%('      %.18e'%finite_loop_me)
+    print '* %-108s *'%''
+    print '='*112

=== modified file 'madgraph/iolibs/template_files/loop_optimized/loop_matrix_standalone.inc'
--- madgraph/iolibs/template_files/loop_optimized/loop_matrix_standalone.inc	2016-09-08 23:15:34 +0000
+++ madgraph/iolibs/template_files/loop_optimized/loop_matrix_standalone.inc	2017-07-29 19:33:46 +0000
@@ -2149,6 +2149,15 @@
 		  ENDDO
 		ENDIF
 
+c       If NaN are present in the evaluation, automatically set the accuracy to 1.0d99.
+        DO I=1,3
+          DO J=1,MAXSTABILITYLENGTH
+            IF (isnan(FULLLIST(I,K,J))) then
+              ACC(K) = 1.0d99
+            ENDIF
+          ENDDO
+        ENDDO
+
 	  enddo
 
       end

=== added file 'madgraph/iolibs/template_files/makefile_sa_f2py'
--- madgraph/iolibs/template_files/makefile_sa_f2py	1970-01-01 00:00:00 +0000
+++ madgraph/iolibs/template_files/makefile_sa_f2py	2017-07-29 19:33:46 +0000
@@ -0,0 +1,23 @@
+include ../Source/make_opts
+SHELL = /bin/bash
+LIBDIR = ../lib/
+PROG   = check
+LINKLIBS =  -L../lib/ -ldhelas -lmodel
+LIBS = $(LIBDIR)libdhelas.$(libext) $(LIBDIR)libmodel.$(libext)
+PROCESS=  allmatrix.o
+LIBRARY = libmatrix.a
+
+# For python linking (require f2py part of numpy)
+ifeq ($(origin MENUM),undefined)
+  MENUM=2
+endif
+
+
+allmatrix$(MENUM)py.so: $(LIBDIR)/$(LIBRARY)  all_matrix.o
+	$(F2PY) $(LINKLIBS) -lmatrix  -c all_matrix.f */matrix.f -m allmatrix$(MENUM)py
+
+
+
+$(LIBDIR)/$(LIBRARY): $(patsubst %.f,%.o,$(wildcard */matrix.f))  all_matrix.o
+	$(call CREATELIB, $@, $^)  
+

=== modified file 'madgraph/iolibs/template_files/matrix_standalone_splitOrders_v4.inc'
--- madgraph/iolibs/template_files/matrix_standalone_splitOrders_v4.inc	2016-03-22 20:44:25 +0000
+++ madgraph/iolibs/template_files/matrix_standalone_splitOrders_v4.inc	2017-07-29 19:33:46 +0000
@@ -292,7 +292,7 @@
 
       END
 
-      SUBROUTINE %(proc_prefix)sGET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE %(proc_prefix)sGET_value(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE   
 C
 C CONSTANT
@@ -325,10 +325,10 @@
       return 
       end
 
-      SUBROUTINE %(proc_prefix)sINITIALISE(PATH)
+      SUBROUTINE %(proc_prefix)sINITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE   
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH 
       call setpara(PATH)  !first call to setup the paramaters    
       return 

=== modified file 'madgraph/iolibs/template_files/matrix_standalone_v4.inc'
--- madgraph/iolibs/template_files/matrix_standalone_v4.inc	2016-03-22 20:44:25 +0000
+++ madgraph/iolibs/template_files/matrix_standalone_v4.inc	2017-07-29 19:33:46 +0000
@@ -66,6 +66,8 @@
 C LOCAL VARIABLES 
 C  
       INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+C  put in common block to expose this variable to python interface
+      COMMON/%(proc_prefix)sPROCESS_NHEL/NHEL 
       REAL*8 T
       REAL*8 %(proc_prefix)sMATRIX
       INTEGER IHEL,IDEN, I, J
@@ -208,7 +210,7 @@
 %(amp2_lines)s
       END
 
-      SUBROUTINE %(proc_prefix)sGET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE %(proc_prefix)sGET_value(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE   
 C
 C CONSTANT
@@ -241,10 +243,10 @@
       return 
       end
 
-      SUBROUTINE %(proc_prefix)sINITIALISE(PATH)
+      SUBROUTINE %(proc_prefix)sINITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE   
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH 
       call setpara(PATH)  !first call to setup the paramaters    
       return 

=== modified file 'madgraph/loop/loop_diagram_generation.py'
--- madgraph/loop/loop_diagram_generation.py	2016-08-25 12:35:03 +0000
+++ madgraph/loop/loop_diagram_generation.py	2017-07-29 19:33:46 +0000
@@ -384,7 +384,7 @@
         # By default the user filter does nothing if filter is not set, 
         # if you want to turn it on and edit it by hand, then set the 
         # variable edit_filter_manually to True
-        edit_filter_manually = False
+        edit_filter_manually = False 
         if not edit_filter_manually and filter in [None,'None']:
             return
         if isinstance(filter,str) and  filter.lower() == 'true':
@@ -414,11 +414,19 @@
                 except Exception as e:
                     raise InvalidCmd("The user-defined filter '%s' did not"%filter+
                                  " returned the following error:\n       > %s"%str(e))
+
+#            if any([abs(pdg) not in range(1,7) for pdg in diag.get_loop_lines_pdgs()]):
+#                valid_diag = False
+
 #            if any([abs(i)!=1000021 for i in diag.get_loop_lines_pdgs()]):
 #                valid_diag=False
 #            if len(diag.get_loop_lines_pdgs())<4:
 #                    valid_diag = False
 
+#            connected_id = diag.get_pdgs_attached_to_loop(structs)
+#            if connected_id.count(22)!=2 or not all(abs(pdg) in range(7) for pdg in diag.get_loop_lines_pdgs()):
+#                valid_diag=False
+
              # Ex. 0: Chose a specific diagram number, here the 8th one for ex.     
 #            if i not in [31]:
 #                valid_diag = False                

=== modified file 'madgraph/loop/loop_exporters.py'
--- madgraph/loop/loop_exporters.py	2017-03-07 15:43:44 +0000
+++ madgraph/loop/loop_exporters.py	2017-07-29 19:33:46 +0000
@@ -80,14 +80,15 @@
         ProcessExporterFortran but give access to arguments like dir_path and
         clean using options. This avoids method resolution object ambiguity"""
 
-    default_opt = {'clean': False, 'complex_mass':False,
+    default_opt = dict(export_v4.ProcessExporterFortran.default_opt)
+    default_opt.update({'clean': False, 'complex_mass':False,
                         'export_format':'madloop', 'mp':True,
                         'loop_dir':'', 'cuttools_dir':'', 
                         'fortran_compiler':'gfortran',
                         'SubProc_prefix': 'P',
                         'output_dependencies': 'external',
                         'compute_color_flows': False,
-                        'mode':''}
+                        'mode':''})
 
     include_names    = {'ninja' : 'mninja.mod',
                         'golem' : 'generic_function_1p.mod',
@@ -242,6 +243,133 @@
 
         self.loop_additional_template_setup()
     
+    def write_f2py_makefile(self):
+        return
+    
+    def write_f2py_check_sa(self, matrix_element, output_path):
+        """ Write the general check_sa.py in SubProcesses that calls all processes successively."""
+ 
+        # No need to further edit this file for now.
+        file = open(os.path.join(self.template_dir,\
+                                       'check_sa_all.py.inc')).read()
+        open(output_path,'w').writelines(file)
+        # Make it executable
+        os.chmod(output_path, os.stat(output_path).st_mode | stat.S_IEXEC)
+    
+    def write_f2py_splitter(self):
+        """write a function to call the correct matrix element"""
+        
+        template = """
+%(python_information)s
+
+      SUBROUTINE INITIALISE(PATH)
+C     ROUTINE FOR F2PY to read the benchmark point.
+      IMPLICIT NONE
+      CHARACTER*512 PATH
+CF2PY INTENT(IN) :: PATH
+      CALL SETPARA(PATH)  !first call to setup the paramaters
+      RETURN
+      END
+
+      SUBROUTINE SET_MADLOOP_PATH(PATH)
+C     Routine to set the path of the folder 'MadLoop5_resources' to MadLoop
+        CHARACTER(512) PATH
+CF2PY intent(in)::path
+        CALL SETMADLOOPPATH(PATH)
+      END
+
+  subroutine smatrixhel(pdgs, npdg, p, ALPHAS, SCALES2, nhel, ANS, RETURNCODE)
+  IMPLICIT NONE
+
+CF2PY double precision, intent(in), dimension(0:3,npdg) :: p
+CF2PY integer, intent(in), dimension(npdg) :: pdgs
+CF2PY integer, intent(in) :: npdg
+CF2PY double precision, intent(out) :: ANS
+CF2PY integer, intent(out) :: RETURNCODE
+CF2PY double precision, intent(in) :: ALPHAS
+CF2PY double precision, intent(in) :: SCALES2
+
+  integer pdgs(*)
+  integer npdg, nhel, RETURNCODE
+  double precision p(*)
+  double precision ANS, ALPHAS, PI,SCALES2
+
+%(smatrixhel)s
+
+      return
+      end
+  
+  subroutine get_pdg_order(OUT)
+  IMPLICIT NONE
+CF2PY INTEGER, intent(out) :: OUT(%(nb_me)i,%(maxpart)i)  
+  
+  INTEGER OUT(%(nb_me)i,%(maxpart)i), PDGS(%(nb_me)i,%(maxpart)i)
+  DATA PDGS/ %(pdgs)s /
+  OUT=PDGS
+  RETURN
+  END
+  
+  subroutine get_prefix(PREFIX)
+  IMPLICIT NONE
+CF2PY CHARACTER*20, intent(out) :: PREFIX(%(nb_me)i)
+  character*20 PREFIX(%(nb_me)i),PREF(%(nb_me)i)
+  DATA PREF / '%(prefix)s'/
+  PREFIX = PREF
+  RETURN
+  END 
+  
+        """
+         
+        allids = self.prefix_info.keys()
+        allprefix = [self.prefix_info[key][0] for key in allids]
+        min_nexternal = min([len(ids) for ids in allids])
+        max_nexternal = max([len(ids) for ids in allids])
+
+        info = []
+        for key, (prefix, tag) in self.prefix_info.items():
+            info.append('#PY %s : %s # %s' % (tag, key, prefix))
+            
+
+        text = []
+        for n_ext in range(min_nexternal, max_nexternal+1):
+            current = [ids for ids in allids if len(ids)==n_ext]
+            if not current:
+                continue
+            if min_nexternal != max_nexternal:
+                if n_ext == min_nexternal:
+                    text.append('       if (npdg.eq.%i)then' % n_ext)
+                else:
+                    text.append('       else if (npdg.eq.%i)then' % n_ext)
+            for ii,pdgs in enumerate(current):
+                condition = '.and.'.join(['%i.eq.pdgs(%i)' %(pdg, i+1) for i, pdg in enumerate(pdgs)])
+                if ii==0:
+                    text.append( ' if(%s) then ! %i' % (condition, i))
+                else:
+                    text.append( ' else if(%s) then ! %i' % (condition,i))
+                text.append(' call %sget_me(p, ALPHAS, DSQRT(SCALES2), NHEL, ANS, RETURNCODE)' % self.prefix_info[pdgs][0])
+            text.append(' endif')
+        #close the function
+        if min_nexternal != max_nexternal:
+            text.append('endif')
+
+        formatting = {'python_information':'\n'.join(info), 
+                          'smatrixhel': '\n'.join(text),
+                          'maxpart': max_nexternal,
+                          'nb_me': len(allids),
+                          'pdgs': ','.join([str(pdg[i]) if i<len(pdg) else '0' 
+                                             for i in range(max_nexternal) \
+                                             for pdg in allids]),
+                      'prefix':'\',\''.join(allprefix)
+                      }
+    
+    
+        text = template % formatting
+        fsock = writers.FortranWriter(pjoin(self.dir_path, 'SubProcesses', 'all_matrix.f'),'w')
+        fsock.writelines(text)
+        fsock.close()
+        
+    
+    
     def loop_additional_template_setup(self, copy_Source_makefile = True):
         """ Perform additional actions specific for this class when setting
         up the template with the copy_template function."""
@@ -832,6 +960,11 @@
         dict['proc_prefix'] = self.get_ME_identifier(matrix_element,
                        group_number = group_number, group_elem_number = proc_id)
 
+        if 'prefix' in self.cmd_options and self.cmd_options['prefix'] in ['int','proc']:
+            for proc in matrix_element.get('processes'):
+                ids = [l.get('id') for l in proc.get('legs_with_decays')]
+                self.prefix_info[tuple(ids)] = [dict['proc_prefix'], proc.get_tag()]
+
         # The proc_id is used for MadEvent grouping, so none of our concern here
         # and it is simply set to an empty string.        
         dict['proc_id'] = ''
@@ -1071,9 +1204,15 @@
         # We can always write the f2py wrapper if present (in loop optimized mode, it is)
         if not os.path.isfile(pjoin(self.template_dir,'check_py.f.inc')):
             return
+        
         file = open(os.path.join(self.template_dir,\
                                       'check_py.f.inc')).read()
-        file=file%replace_dict
+
+        if 'prefix' in self.cmd_options and self.cmd_options['prefix'] in ['int','proc']:
+            replace_dict['prefix_routine'] = replace_dict['proc_prefix']
+        else:
+            replace_dict['prefix_routine'] = ''
+        file=file%replace_dict        
         new_path = writer.name.replace('check_sa.f', 'f2py_wrapper.f')
         new_writer = writer.__class__(new_path, 'w')
         new_writer.writelines(file)
@@ -1088,7 +1227,8 @@
 p= [[None,]*4]*%d"""%len(curr_proc.get('legs'))
 
         process_definition_string = curr_proc.nice_string().replace('Process:','')
-        file=file.format(random_PSpoint_python_formatted,process_definition_string)
+        file=file.format(random_PSpoint_python_formatted,process_definition_string,
+                         replace_dict['proc_prefix'].lower())
         new_path = writer.name.replace('check_sa.f', 'check_sa.py')
         new_writer = open(new_path, 'w')
         new_writer.writelines(file)

=== modified file 'madgraph/various/banner.py'
--- madgraph/various/banner.py	2017-03-03 10:01:29 +0000
+++ madgraph/various/banner.py	2017-07-29 19:33:46 +0000
@@ -333,7 +333,6 @@
 
     def modify_init_cross(self, cross):
         """modify the init information with the associate cross-section"""
-
         assert isinstance(cross, dict)
 #        assert "all" in cross
         assert "init" in self

=== modified file 'madgraph/various/lhe_parser.py'
--- madgraph/various/lhe_parser.py	2017-03-04 20:26:57 +0000
+++ madgraph/various/lhe_parser.py	2017-07-29 19:33:46 +0000
@@ -2019,7 +2019,7 @@
 
     @property
     def mass(self):
-        """return the mass"""    
+        """return the mass"""
         return math.sqrt(self.E**2 - self.px**2 - self.py**2 - self.pz**2)
 
     @property
@@ -2029,7 +2029,7 @@
 
     @property
     def pt(self):
-        return math.sqrt(max(0, self.pt2()))
+        return math.sqrt(max(0, self.pt2))
     
     @property
     def pseudorapidity(self):
@@ -2041,7 +2041,7 @@
         return  0.5* math.log((self.E +self.pz) / (self.E - self.pz))
     
     
-    
+    @property
     def pt2(self):
         """ return the pt square """
         
@@ -2314,11 +2314,13 @@
 class NLO_PARTIALWEIGHT(object):
 
     class BasicEvent(list):
+
         
         def __init__(self, momenta, wgts, event, real_type=(1,11)):
+            
             list.__init__(self, momenta)
-            
             assert self
+            self.soft = False
             self.wgts = wgts
             self.pdgs = list(wgts[0].pdgs)
             self.event = event
@@ -2338,58 +2340,130 @@
                 self.pop(ind2)
                 self.pdgs.pop(ind1) 
                 self.pdgs.insert(ind1, wgts[0].merge_new_pdg )
-                self.pdgs.pop(ind2)                 
+                self.pdgs.pop(ind2)
+                #ensure the onshell of the particle
+                self.put_onshell(ind1, 0)
                 # DO NOT update the pdgs of the partial weight!
+
             elif any(w.type in self.real_type for w in wgts):
                 if any(w.type not in self.real_type for w in wgts):
                     raise Exception
-                # check if this is too soft/colinear if so use the born
-                ind1, ind2 = [ind-1 for ind in wgts[0].to_merge_pdg] 
-                if ind1> ind2: 
-                    ind1, ind2 = ind2, ind1                
-                if ind1 >= sum(1 for p in event if p.status==-1):
-                    new_p = self[ind1] + self[ind2]
-                else:
-                    new_p = self[ind1] - self[ind2]
-
-                if __debug__:
-                    ptot = FourMomentum()
-                    for i in xrange(len(self)):
-                        if i <2:
-                            ptot += self[i]
-                        else:
-                            ptot -= self[i]
-                    if ptot.mass_sqr > 1e-16:
-                        misc.sprint(ptot, ptot.mass_sqr)
-                
-                inv_mass = new_p.mass_sqr
-                shat = (self[0]+self[1]).mass_sqr
-                if (abs(inv_mass)/shat < 1e-6):
-#                    misc.sprint(abs(inv_mass)/shat)
-                    self.pop(ind1) 
-                    self.insert(ind1, new_p)
-                    self.pop(ind2)
-                    self.pdgs.pop(ind1) 
-                    self.pdgs.insert(ind1, wgts[0].merge_new_pdg )
-                    self.pdgs.pop(ind2)                 
-                    # DO NOT update the pdgs of the partial weight!                    
-                
-                if __debug__:
-                    ptot = FourMomentum()
-                    for i in xrange(len(self)):
-                        if i <2:
-                            ptot += self[i]
-                        else:
-                            ptot -= self[i]
-                    if ptot.mass_sqr > 1e-16:
-                        misc.sprint(ptot, ptot.mass_sqr)
-#                            raise Exception
+                # Do nothing !!!
+                # previously (commented we were checking here if the particle 
+                # were too soft this is done  later now
+#                    The comment line below allow to convert this event 
+#                    to a born one (old method)    
+#                    self.pop(ind1) 
+#                    self.insert(ind1, new_p)
+#                    self.pop(ind2)
+#                    self.pdgs.pop(ind1) 
+#                    self.pdgs.insert(ind1, wgts[0].merge_new_pdg )
+#                    self.pdgs.pop(ind2)                 
+#                    # DO NOT update the pdgs of the partial weight!                    
             else:
                 raise Exception
+
+        def check_fks_singularity(self, ind1, ind2, nb_init=2, threshold=1e-6):
+            """check that the propagator associated to ij is not too light 
+               [related to soft-collinear singularity]"""
+
+            if ind1> ind2: 
+                ind1, ind2 = ind2, ind1                
+            if ind1 >= nb_init:
+                new_p = self[ind1] + self[ind2]
+            else:
+                new_p = self[ind1] - self[ind2]
+                
+            inv_mass = new_p.mass_sqr
+            shat = (self[0]+self[1]).mass_sqr
+            if (abs(inv_mass)/shat < threshold):
+                return True
+            else:
+                return False
+ 
  
         def get_pdg_code(self):
             return self.pdgs
-        
+
+        def put_onshell(self, ind, mass):
+            """ check that particle mass of particle number 'ind' 
+                is onshell and if not use algorithm from VH thesis
+                to restore them."""
+            
+            if misc.equal(self[ind].mass_sqr, mass**2, 20):
+                return # nothing to do
+
+            if sum(1 for p in self.event if p.status==-1) == 1:
+                logger.debug('1 to N can not be put back onshell. Continue')
+                return # does not work for 1 to N: keep it offshell
+
+            assert self.pdgs[ind] in (range(-6,7)+[21,22]) ,"onshell routine works only for zero mass particle. PDG is not assume massless."
+            
+            if ind<2:
+                self.put_onshell_initial(ind, mass)
+            else:
+                raise Exception, 'code not implemented for final state'
+#                self.put_onshell_final(ind, mass)
+            
+        def put_onshell_initial(self, ind, mass):
+            """ """
+
+            assert mass == 0
+
+            ptot = self[0] + self[1]
+                
+            # particle 2 should not have any pt. The one reshuffle has a pt
+            if ind == 0:
+                p1 = self[0] # the one we modify with pt
+                p2 = self[1] # original one no pt
+            else:
+                p1 = self[1] # the one we modify with pt
+                p2 = self[0] # original one no pt
+                
+            assert misc.equal(p2.px, 0 , 8)
+            assert misc.equal(p2.py, 0 , 8)
+
+            discr = -ptot.mass_sqr
+            
+            shifte1 = ( ptot.E*(ptot.pt2 -2*p1.E*ptot.E +ptot.E**2) +
+                        (2*p1.E-ptot.E)*ptot.pz**2 + ptot.pz*discr
+                      )/ (2*(ptot.E**2-ptot.pz**2))
+             
+            shifte2 = -( ptot.E*(ptot.pt2 + 2*p2.E*ptot.E - ptot.E**2) +
+                        (-2*p2.E + ptot.E)*ptot.pz**2 + ptot.pz*discr
+                      )/ (2*(ptot.E**2-ptot.pz**2))
+            
+            shiftz1 = ( -2*p1.pz*(ptot.E**2 - ptot.pz**2) 
+                        + ptot.pz*(ptot.pt2 + ptot.E**2 - ptot.pz**2)
+                        + ptot.E*discr
+                      )/ (2*(ptot.E**2 - ptot.pz**2))
+
+
+            shiftz2 = -( 2*p2.pz*(ptot.E**2 - ptot.pz**2) 
+                         + ptot.pz*(ptot.pt2 - ptot.E**2 + ptot.pz**2)
+                         + ptot.E*discr
+                       )/ (2*(ptot.E**2 - ptot.pz**2))
+
+     
+            p1 = FourMomentum(p1.E + shifte1,
+                                   p1.px,
+                                   p1.py,
+                                   p1.pz + shiftz1)
+                                
+            p2 = FourMomentum(p2.E + shifte2,
+                                   0.,
+                                   0.,
+                                   p2.pz + shiftz2)
+
+            if ind == 0:
+                self[0], self[1] = p1, p2
+            else:
+                self[0], self[1] = p2, p1
+                
+            assert misc.equal(self[0].mass_sqr, 0, 8)
+            assert misc.equal(self[1].mass_sqr, 0, 8)
+                 
+            
         def get_tag_and_order(self):
             """ return the tag and order for this basic event""" 
             (initial, _), _ = self.event.get_tag_and_order()
@@ -2403,7 +2477,7 @@
         
         def get_momenta(self, get_order, allow_reversed=True):
             """return the momenta vector in the order asked for"""
-        
+             
             #avoid to modify the input
             order = [list(get_order[0]), list(get_order[1])] 
             out = [''] *(len(order[0])+len(order[1]))
@@ -2518,6 +2592,7 @@
 #0.359257891118d-05 0.000000000000d+00 0.000000000000d+00  5 21 3 -11 11 3 2 0.12292838d-02 0.43683926d-01 0.58606724d+03 0.58606724d+03 0.58606724d+03  1 12  3 0.334254554762d+00
 #0.929944817736d-03 0.000000000000d+00 0.000000000000d+00  5 21 3 -11 11 3 2 0.12112732d-02 0.45047393d-01 0.58606724d+03 0.58606724d+03 0.58606724d+03  2 11  3 0.835109616010d+02
         
+        
         text = text.lower().replace('d','e')
         all_line = text.split('\n')
         #get global information
@@ -2541,27 +2616,47 @@
         
         assert len(wgts) == int(nb_wgt)
         
-        get_weights_for_momenta = {}
+        get_weights_for_momenta = dict( (i,[]) for i in range(1,nb_event+1)  )
         size_momenta = 0
         for wgt in wgts:
             if wgt.momenta_config in get_weights_for_momenta:
                 get_weights_for_momenta[wgt.momenta_config].append(wgt)
             else: 
+                misc.sprint(wgt.type, wgt.real_related,wgt.momenta_config)
                 if size_momenta == 0: size_momenta = wgt.nexternal
                 assert size_momenta == wgt.nexternal
                 get_weights_for_momenta[wgt.momenta_config] = [wgt]
     
         assert sum(len(c) for c in get_weights_for_momenta.values()) == int(nb_wgt), "%s != %s" % (sum(len(c) for c in get_weights_for_momenta.values()), nb_wgt)
     
+        # check singular behavior
+        for key in range(1, nb_event+1):
+            wgts = get_weights_for_momenta[key]
+            if not wgts:
+                continue
+            if size_momenta == 0: size_momenta = wgts[0].nexternal
+            p = momenta[size_momenta*(key-1):key*size_momenta]
+            evt = self.BasicEvent(p, wgts, self.event, self.real_type) 
+            if len(evt) == size_momenta: #real type 
+                for wgt in wgts:
+                    if evt.check_fks_singularity(wgt.to_merge_pdg[0]-1,
+                                                 wgt.to_merge_pdg[1]-1,
+                                                 nb_init=sum(1 for p in self.event if p.status==-1)):
+                        get_weights_for_momenta[wgt.momenta_config].remove(wgt)
+                        get_weights_for_momenta[wgt.born_related].append(wgt)
+     
          
-    
+        assert sum(len(c) for c in get_weights_for_momenta.values()) == int(nb_wgt), "%s != %s" % (sum(len(c) for c in get_weights_for_momenta.values()), nb_wgt)
+           
         self.cevents = []   
         for key in range(1, nb_event+1): 
             if key in get_weights_for_momenta:
                 wgt = get_weights_for_momenta[key]
-                evt = self.BasicEvent(momenta[:size_momenta], get_weights_for_momenta[key], self.event, self.real_type) 
+                if not wgt:
+                    continue
+                p = momenta[size_momenta*(key-1):key*size_momenta]
+                evt = self.BasicEvent(p, get_weights_for_momenta[key], self.event, self.real_type)                         
                 self.cevents.append(evt)
-            momenta = momenta[size_momenta:]
            
         nb_wgt_check = 0 
         for cevt in self.cevents:

=== modified file 'madgraph/various/misc.py'
--- madgraph/various/misc.py	2017-03-03 23:41:19 +0000
+++ madgraph/various/misc.py	2017-07-29 19:33:46 +0000
@@ -1689,7 +1689,23 @@
     
     
 
-
+def plugin_import(module, error_msg, fcts=[]):
+    """convenient way to import a plugin file/function"""
+    
+    try:
+        _temp = __import__('PLUGIN.%s' % module, globals(), locals(), fcts, -1)
+    except ImportError:
+        try:
+            _temp = __import__('MG5aMC_PLUGIN.%s' % module, globals(), locals(), fcts, -1)
+        except ImportError:
+            raise MadGraph5Error, error_msg
+    
+    if not fcts:
+        return _temp
+    elif len(fcts) == 1:
+        return getattr(_temp,fcts[0])
+    else:
+        return [getattr(_temp,name) for name in fcts]
 
 
 python_lhapdf=None

=== modified file 'models/template_files/fortran/lha_read.f'
--- models/template_files/fortran/lha_read.f	2014-10-23 00:33:10 +0000
+++ models/template_files/fortran/lha_read.f	2017-07-29 19:33:46 +0000
@@ -259,12 +259,12 @@
       enddo
       
       npara=npara-1
-c      close(iunit)
+99      close(iunit)
       if (WriteParamLog) then
         close(logfile)
       endif
 
- 99   return
+      return
  
       end
 

=== modified file 'models/template_files/fortran/lha_read_mp.f'
--- models/template_files/fortran/lha_read_mp.f	2014-10-23 00:33:10 +0000
+++ models/template_files/fortran/lha_read_mp.f	2017-07-29 19:33:46 +0000
@@ -279,12 +279,12 @@
       enddo
       
       npara=npara-1
-c      close(iunit)
+ 99   close(iunit)
       if (WriteParamLog) then
         close(logfile)
       endif
       
- 99   return
+      return
  
       end
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f	2017-07-29 19:33:46 +0000
@@ -325,7 +325,7 @@
 
       END
 
-      SUBROUTINE GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -358,10 +358,10 @@
       RETURN
       END
 
-      SUBROUTINE INITIALISE(PATH)
+      SUBROUTINE INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%born_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%born_matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%born_matrix.f	2017-07-29 19:33:46 +0000
@@ -325,7 +325,7 @@
 
       END
 
-      SUBROUTINE GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -358,10 +358,10 @@
       RETURN
       END
 
-      SUBROUTINE INITIALISE(PATH)
+      SUBROUTINE INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_standalone/matrix.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_standalone/matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_standalone/matrix.f	2017-07-29 19:33:46 +0000
@@ -68,6 +68,8 @@
 C     LOCAL VARIABLES 
 C     
       INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+C     put in common block to expose this variable to python interface
+      COMMON/PROCESS_NHEL/NHEL
       REAL*8 T
       REAL*8 MATRIX
       INTEGER IHEL,IDEN, I, J
@@ -275,7 +277,7 @@
 
       END
 
-      SUBROUTINE GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -308,10 +310,10 @@
       RETURN
       END
 
-      SUBROUTINE INITIALISE(PATH)
+      SUBROUTINE INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%born_matrix.f'
--- tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%born_matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%born_matrix.f	2017-07-29 19:33:46 +0000
@@ -339,7 +339,7 @@
 
       END
 
-      SUBROUTINE ML5_0_GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE ML5_0_GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -372,10 +372,10 @@
       RETURN
       END
 
-      SUBROUTINE ML5_0_INITIALISE(PATH)
+      SUBROUTINE ML5_0_INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_NoSQSO.f'
--- tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_NoSQSO.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_NoSQSO.f	2017-07-29 19:33:46 +0000
@@ -68,6 +68,8 @@
 C     LOCAL VARIABLES 
 C     
       INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+C     put in common block to expose this variable to python interface
+      COMMON/PROCESS_NHEL/NHEL
       REAL*8 T
       REAL*8 MATRIX
       INTEGER IHEL,IDEN, I, J
@@ -451,7 +453,7 @@
 
       END
 
-      SUBROUTINE GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -484,10 +486,10 @@
       RETURN
       END
 
-      SUBROUTINE INITIALISE(PATH)
+      SUBROUTINE INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_QCDsq_le_6.f'
--- tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_QCDsq_le_6.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_QCDsq_le_6.f	2017-07-29 19:33:46 +0000
@@ -537,7 +537,7 @@
 
       END
 
-      SUBROUTINE GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -570,10 +570,10 @@
       RETURN
       END
 
-      SUBROUTINE INITIALISE(PATH)
+      SUBROUTINE INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_ampOrderQED2_eq_2_WGTsq_le_14_QCDsq_gt_4.f'
--- tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_ampOrderQED2_eq_2_WGTsq_le_14_QCDsq_gt_4.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_ampOrderQED2_eq_2_WGTsq_le_14_QCDsq_gt_4.f	2017-07-29 19:33:46 +0000
@@ -537,7 +537,7 @@
 
       END
 
-      SUBROUTINE GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -570,10 +570,10 @@
       RETURN
       END
 
-      SUBROUTINE INITIALISE(PATH)
+      SUBROUTINE INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/long_ML_SMQCD_default/dux_mumvmxg/born_matrix.f'
--- tests/input_files/IOTestsComparison/long_ML_SMQCD_default/dux_mumvmxg/born_matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/long_ML_SMQCD_default/dux_mumvmxg/born_matrix.f	2017-07-29 19:33:46 +0000
@@ -68,6 +68,8 @@
 C     LOCAL VARIABLES 
 C     
       INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+C     put in common block to expose this variable to python interface
+      COMMON/ML5_0_PROCESS_NHEL/NHEL
       REAL*8 T
       REAL*8 ML5_0_MATRIX
       INTEGER IHEL,IDEN, I, J
@@ -264,7 +266,7 @@
 
       END
 
-      SUBROUTINE ML5_0_GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE ML5_0_GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -297,10 +299,10 @@
       RETURN
       END
 
-      SUBROUTINE ML5_0_INITIALISE(PATH)
+      SUBROUTINE ML5_0_INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/long_ML_SMQCD_default/gg_wmtbx/born_matrix.f'
--- tests/input_files/IOTestsComparison/long_ML_SMQCD_default/gg_wmtbx/born_matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/long_ML_SMQCD_default/gg_wmtbx/born_matrix.f	2017-07-29 19:33:46 +0000
@@ -68,6 +68,8 @@
 C     LOCAL VARIABLES 
 C     
       INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+C     put in common block to expose this variable to python interface
+      COMMON/ML5_0_PROCESS_NHEL/NHEL
       REAL*8 T
       REAL*8 ML5_0_MATRIX
       INTEGER IHEL,IDEN, I, J
@@ -300,7 +302,7 @@
 
       END
 
-      SUBROUTINE ML5_0_GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE ML5_0_GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -333,10 +335,10 @@
       RETURN
       END
 
-      SUBROUTINE ML5_0_INITIALISE(PATH)
+      SUBROUTINE ML5_0_INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/dux_mumvmxg/born_matrix.f'
--- tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/dux_mumvmxg/born_matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/dux_mumvmxg/born_matrix.f	2017-07-29 19:33:46 +0000
@@ -68,6 +68,8 @@
 C     LOCAL VARIABLES 
 C     
       INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+C     put in common block to expose this variable to python interface
+      COMMON/ML5_0_PROCESS_NHEL/NHEL
       REAL*8 T
       REAL*8 ML5_0_MATRIX
       INTEGER IHEL,IDEN, I, J
@@ -264,7 +266,7 @@
 
       END
 
-      SUBROUTINE ML5_0_GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE ML5_0_GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -297,10 +299,10 @@
       RETURN
       END
 
-      SUBROUTINE ML5_0_INITIALISE(PATH)
+      SUBROUTINE ML5_0_INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/gg_wmtbx/born_matrix.f'
--- tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/gg_wmtbx/born_matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/gg_wmtbx/born_matrix.f	2017-07-29 19:33:46 +0000
@@ -68,6 +68,8 @@
 C     LOCAL VARIABLES 
 C     
       INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+C     put in common block to expose this variable to python interface
+      COMMON/ML5_0_PROCESS_NHEL/NHEL
       REAL*8 T
       REAL*8 ML5_0_MATRIX
       INTEGER IHEL,IDEN, I, J
@@ -300,7 +302,7 @@
 
       END
 
-      SUBROUTINE ML5_0_GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE ML5_0_GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -333,10 +335,10 @@
       RETURN
       END
 
-      SUBROUTINE ML5_0_INITIALISE(PATH)
+      SUBROUTINE ML5_0_INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/short_ML_SMQCD_default/ddx_ttx/born_matrix.f'
--- tests/input_files/IOTestsComparison/short_ML_SMQCD_default/ddx_ttx/born_matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/short_ML_SMQCD_default/ddx_ttx/born_matrix.f	2017-07-29 19:33:46 +0000
@@ -68,6 +68,8 @@
 C     LOCAL VARIABLES 
 C     
       INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+C     put in common block to expose this variable to python interface
+      COMMON/ML5_0_PROCESS_NHEL/NHEL
       REAL*8 T
       REAL*8 ML5_0_MATRIX
       INTEGER IHEL,IDEN, I, J
@@ -247,7 +249,7 @@
 
       END
 
-      SUBROUTINE ML5_0_GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE ML5_0_GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -280,10 +282,10 @@
       RETURN
       END
 
-      SUBROUTINE ML5_0_INITIALISE(PATH)
+      SUBROUTINE ML5_0_INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/short_ML_SMQCD_default/gg_ttx/born_matrix.f'
--- tests/input_files/IOTestsComparison/short_ML_SMQCD_default/gg_ttx/born_matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/short_ML_SMQCD_default/gg_ttx/born_matrix.f	2017-07-29 19:33:46 +0000
@@ -68,6 +68,8 @@
 C     LOCAL VARIABLES 
 C     
       INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+C     put in common block to expose this variable to python interface
+      COMMON/ML5_0_PROCESS_NHEL/NHEL
       REAL*8 T
       REAL*8 ML5_0_MATRIX
       INTEGER IHEL,IDEN, I, J
@@ -253,7 +255,7 @@
 
       END
 
-      SUBROUTINE ML5_0_GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE ML5_0_GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -286,10 +288,10 @@
       RETURN
       END
 
-      SUBROUTINE ML5_0_INITIALISE(PATH)
+      SUBROUTINE ML5_0_INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/ddx_ttx/born_matrix.f'
--- tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/ddx_ttx/born_matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/ddx_ttx/born_matrix.f	2017-07-29 19:33:46 +0000
@@ -68,6 +68,8 @@
 C     LOCAL VARIABLES 
 C     
       INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+C     put in common block to expose this variable to python interface
+      COMMON/ML5_0_PROCESS_NHEL/NHEL
       REAL*8 T
       REAL*8 ML5_0_MATRIX
       INTEGER IHEL,IDEN, I, J
@@ -247,7 +249,7 @@
 
       END
 
-      SUBROUTINE ML5_0_GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE ML5_0_GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -280,10 +282,10 @@
       RETURN
       END
 
-      SUBROUTINE ML5_0_INITIALISE(PATH)
+      SUBROUTINE ML5_0_INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/gg_ttx/born_matrix.f'
--- tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/gg_ttx/born_matrix.f	2016-09-15 20:20:30 +0000
+++ tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/gg_ttx/born_matrix.f	2017-07-29 19:33:46 +0000
@@ -68,6 +68,8 @@
 C     LOCAL VARIABLES 
 C     
       INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+C     put in common block to expose this variable to python interface
+      COMMON/ML5_0_PROCESS_NHEL/NHEL
       REAL*8 T
       REAL*8 ML5_0_MATRIX
       INTEGER IHEL,IDEN, I, J
@@ -253,7 +255,7 @@
 
       END
 
-      SUBROUTINE ML5_0_GET_ME(P, ALPHAS, NHEL ,ANS)
+      SUBROUTINE ML5_0_GET_VALUE(P, ALPHAS, NHEL ,ANS)
       IMPLICIT NONE
 C     
 C     CONSTANT
@@ -286,10 +288,10 @@
       RETURN
       END
 
-      SUBROUTINE ML5_0_INITIALISE(PATH)
+      SUBROUTINE ML5_0_INITIALISEMODEL(PATH)
 C     ROUTINE FOR F2PY to read the benchmark point.    
       IMPLICIT NONE
-      CHARACTER*180 PATH
+      CHARACTER*512 PATH
 CF2PY INTENT(IN) :: PATH
       CALL SETPARA(PATH)  !first call to setup the paramaters    
       RETURN