Merge lp:~maddevelopers/mg5amcnlo/2.3.3_test into lp:mg5amcnlo/lts

[Olivier Mattelaer]

This is just a test to see if the procedure that I applied fix the conflict.
I will not perform this merge (if it works, I will push this insice 2.3.3)

[Olivier Mattelaer (olivier-mattelaer)]

The test seems succesfull.
So closing it and I'm going to apply it to 2.3.3

review: Disapprove

Unmerged revisions

326. By Olivier Mattelaer on 2015-10-20

merge with 2.3.2.2

=== modified file 'MadSpin/__init__.py' (properties changed: -x to +x)
=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2015-08-20 10:00:00 +0000
+++ MadSpin/decay.py	2015-10-20 14:02:27 +0000
@@ -906,7 +906,7 @@
             if init[0] in self.banner.param_card['decay'].decay_table:
                 br *= self.banner.param_card['decay'].decay_table[init[0]].get(lhaid).value
                 br *= self.get_br(decay)
-            else:
+            elif -init[0] in self.banner.param_card['decay'].decay_table:
                 init = -init[0]
                 lhaid=[x if self.model.get_particle(x)['self_antipart'] else -x
                        for x in final]
@@ -914,6 +914,13 @@
                 lhaid = tuple([len(final)] + lhaid)
                 br *= self.banner.param_card['decay'].decay_table[init].get(lhaid).value
                 br *= self.get_br(decay)
+            elif init[0] not in self.decay_ids and -init[0] not in self.decay_ids:
+                logger.warning("No Branching ratio applied for %s. Please check if this is expected" % init[0])
+                br *= self.get_br(decay)
+            else:
+                raise MadGraph5Error,"No valid decay for %s. No 2 body decay for that particle. (three body are not supported by MadSpin)" % init[0]
+
+                
 
         for decays in ids.values():
             if len(decays) == 1:
@@ -929,7 +936,7 @@
                     except ValueError:
                         break
                 br /= math.factorial(nb)
-
+                
         return br
 
             
@@ -1013,7 +1020,7 @@
         self[tag]['decaying'] = tuple(decaying)
                 
         # sanity check
-        assert self[tag]['total_br'] <= 1.01, self[tag]['total_br']
+        assert self[tag]['total_br'] <= 1.01, "wrong BR for %s: %s " % (tag,self[tag]['total_br'])
 
 
         

=== modified file 'MadSpin/interface_madspin.py' (properties changed: -x to +x)
--- MadSpin/interface_madspin.py	2015-08-04 14:28:22 +0000
+++ MadSpin/interface_madspin.py	2015-10-20 14:02:27 +0000
@@ -34,6 +34,7 @@
 import madgraph.interface.master_interface as master_interface
 import madgraph.various.misc as misc
 import madgraph.iolibs.files as files
+import madgraph.iolibs.export_v4 as export_v4
 import madgraph.various.banner as banner
 import madgraph.various.lhe_parser as lhe_parser
 
@@ -229,15 +230,24 @@
             args.remove('--bypass_check')
             bypass_check = True
             
-        if len(args) != 2:     
+        if len(args) == 1:  
+            logger.warning("""No param_card defined for the new model. We will use the default one but this might completely wrong.""")
+        elif len(args) != 2:
             return self.InvalidCmd, 'import model requires two arguments'
         
         model_name = args[0]
-        card = args[1]
-        if not os.path.exists(card):
-            raise self.InvalidCmd('%s: no such file' % card)
-
         self.load_model(model_name, False, False)
+        
+        if len(args) == 2:
+            card = args[1]
+            if not os.path.exists(card):
+                raise self.InvalidCmd('%s: no such file' % card)
+        else:
+            card = "madspin_param_card.dat"
+            export_v4.UFO_model_to_mg4.create_param_card_static(self.model,
+                                card, rule_card_path=None)
+
+        
 
         #Check the param_card
         if not bypass_check:
@@ -255,11 +265,15 @@
             if diff:
                 raise self.InvalidCmd('''Original param_card differs on some parameters:
     %s
-    Due to those preferences, we prefer not to proceed. If you are sure about what you are doing, 
+    Due to those differences, we prefer not to proceed. If you are sure about what you are doing, 
     you can use the command \'import model MODELNAME PARAM_CARD_PATH --bypass_check\''''
     % diff.replace('\n','\n    '))
+   
+   
                 
-        #OK load the new param_card
+        #OK load the new param_card (but back up the old one)
+        if 'slha' in self.banner:
+            self.banner['slha_original'] = self.banner['slha']
         self.banner['slha'] = open(card).read()
         if hasattr(self.banner, 'param_card'):
             del self.banner.param_card
@@ -749,7 +763,7 @@
         # 2. Generate the events requested
         with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
             mg5 = self.mg5cmd
-            modelpath = self.model.get('modelpath')
+            modelpath = self.model.get('modelpath+restriction')
             mg5.exec_cmd("import model %s" % modelpath)      
             to_event = {}
             for pdg, nb_needed in to_decay.items():
@@ -923,12 +937,11 @@
         #base_model = import_ufo.import_model(model_path)
 
         # Import model
-        base_model = import_ufo.import_model(name, decay=True)
+        base_model = import_ufo.import_model(name, decay=True,
+                                               complex_mass_scheme=complex_mass)
 
         if use_mg_default:
             base_model.pass_particles_name_in_mg_default()
-        if complex_mass:
-            base_model.change_mass_to_complex_scheme()
         
         self.model = base_model
         self.mg5cmd._curr_model = self.model

=== modified file 'Template/LO/Source/dsample.f'
--- Template/LO/Source/dsample.f	2015-06-08 16:01:11 +0000
+++ Template/LO/Source/dsample.f	2015-10-20 14:02:27 +0000
@@ -1608,7 +1608,7 @@
       if (iteration .eq. cur_it) then
 c        Add the current point to the DiscreteSamplerGrid
          call add_entry_to_discrete_dimensions(wgt)
-         if (kn.eq.0) then
+         if (kn.eq.0.and.(iteration.eq.1.or.use_cut.eq.-2)) then
             ! ensure that all cumulative variable are at zero (usefull for reset)
             twgt1 = 0d0         !
             iavg = 0            !Vars for averging to increase err estimate
@@ -2156,7 +2156,7 @@
 c     $                 access='append',err=129)
 c                  write(22, 80) real(tmean), real(tsigma), real(chi2)
 c 129              close(22)
-                  tsigma = tsigma*sqrt(chi2) !This gives the 68% confidence cross section
+                  tsigma = tsigma*sqrt(max(0d0,chi2)) !This gives the 68% confidence cross section
                   cur_it = itm+20
                   return
                endif

=== modified file 'Template/LO/Source/kin_functions.f'
--- Template/LO/Source/kin_functions.f	2014-09-11 15:22:26 +0000
+++ Template/LO/Source/kin_functions.f	2015-10-20 14:02:27 +0000
@@ -127,7 +127,11 @@
 c      pm=dsqrt(p(1)**2+p(2)**2+p(3)**2)
       
       pm = p(0)
-      rap = .5d0*dlog((pm+p(3))/(pm-p(3)))+cm_rap
+      if (pm.gt.p(3)) then
+        rap = .5d0*dlog((pm+p(3))/(pm-p(3)))+cm_rap
+      else
+        rap = -1d99
+      endif
       end
       DOUBLE PRECISION  FUNCTION rap2(p)
 c************************************************************************

=== modified file 'Template/LO/bin/internal/Gridpack/compile'
--- Template/LO/bin/internal/Gridpack/compile	2013-11-13 02:21:10 +0000
+++ Template/LO/bin/internal/Gridpack/compile	2015-10-20 14:02:27 +0000
@@ -37,8 +37,8 @@
 fi
 
 # Check for LHAPDF
-c=`awk '/^[^#].*=.*pdlabel/{print $1}' Cards/run_card.dat`
-if [[ $c == "'lhapdf'" ]]; then
+c=`awk '/^[^#].*=[^#!]*pdlabel/{print $1}' Cards/run_card.dat`
+if [[ $c == "'lhapdf'" || $c == "lhapdf" ]]; then
     echo Using LHAPDF interface!
     export lhapdf=true
 else

=== renamed file 'Template/LO/bin/internal/addmasses.py' => 'Template/LO/bin/internal/addmasses_optional.py'
=== modified file 'Template/LO/bin/internal/create_matching_plots.C'
--- Template/LO/bin/internal/create_matching_plots.C	2011-06-29 14:15:21 +0000
+++ Template/LO/bin/internal/create_matching_plots.C	2015-10-20 14:02:27 +0000
@@ -1,4 +1,14 @@
 {
+  const char* file="pythia.root";
+  TFile *f = new TFile(file,"RECREATE");
+  TTree *events = new TTree("events","Events");
+  Long64_t nlines1 = events->ReadFile("events.tree","Npart:DJR1:DJR2:DJR3:DJR4:PTISR:PTFSR:PT2MX:NCJET:IFILE");
+  cout << "Found "<< nlines1 << " events"<<endl;
+  TTree *xsecs = new TTree("xsecs","Xsecs");
+  Long64_t nlines2 = xsecs->ReadFile("xsecs.tree","Xsecfact");
+  cout << "Found "<< nlines2 << " files"<<endl;
+  f->Write();
+
   gROOT->ProcessLine(".x ../bin/internal/plot_tree.C(\"DJR1\")");
   gROOT->ProcessLine(".x ../bin/internal/plot_tree.C(\"DJR2\")");
   gROOT->ProcessLine(".x ../bin/internal/plot_tree.C(\"DJR3\")");

=== modified file 'Template/LO/bin/internal/create_matching_plots.sh'
--- Template/LO/bin/internal/create_matching_plots.sh	2012-07-20 06:20:21 +0000
+++ Template/LO/bin/internal/create_matching_plots.sh	2015-10-20 14:02:27 +0000
@@ -24,9 +24,7 @@
     exit
 fi
 
-echo Running root
-root -q -b -l ../bin/internal/read_tree_files.C &> /dev/null
-echo Creating plots
+echo Creating matching plots
 root -q -b -l ../bin/internal/create_matching_plots.C &> /dev/null
 mv pythia.root $1/$2_pythia.root
 

=== modified file 'Template/LO/bin/internal/plot_tree.C'
--- Template/LO/bin/internal/plot_tree.C	2014-09-24 19:41:18 +0000
+++ Template/LO/bin/internal/plot_tree.C	2015-10-20 14:02:27 +0000
@@ -1,7 +1,11 @@
 #include <string>
 
-bool plot_tree(char* quantity,char* plotdim="(100,0.,3.)",bool log=true)
+bool plot_tree(const char* quantity, const char* plotdim = 0, bool log = true)
 {
+  const char* plotdim_default = "(100,0.,3.)";
+  if(plotdim == 0){
+    plotdim = plotdim_default;
+  }
   char tmp1[250];
   char tmp2[300];
   char tmp3[100];

=== removed file 'Template/LO/bin/internal/read_tree_files.C'
--- Template/LO/bin/internal/read_tree_files.C	2011-06-29 14:15:21 +0000
+++ Template/LO/bin/internal/read_tree_files.C	1970-01-01 00:00:00 +0000
@@ -1,12 +0,0 @@
-{
-  char* file="pythia.root";
-  TFile *f = new TFile(file,"RECREATE");    
-  TTree *events = new TTree("events","Events");
-  Long64_t nlines1 = events->ReadFile("events.tree","Npart:DJR1:DJR2:DJR3:DJR4:PTISR:PTFSR:PT2MX:NCJET:IFILE");
-  cout << "Found "<< nlines1 << " events"<<endl;
-  TTree *xsecs = new TTree("xsecs","Xsecs");                   
-  Long64_t nlines2 = xsecs->ReadFile("xsecs.tree","Xsecfact");
-  cout << "Found "<< nlines2 << " files"<<endl;
-
-  f->Write();
-}

=== modified file 'Template/MadWeight/Python/Info.py'
--- Template/MadWeight/Python/Info.py	2011-10-25 14:42:02 +0000
+++ Template/MadWeight/Python/Info.py	2015-10-20 14:02:27 +0000
@@ -1,11 +1,11 @@
 #!/usr/bin/env python
 
 def giveInfo(class_):
-	if type(class_)!=str:
-		class_=class_.__class__.__name__
-	for info in dir(eval(class_)):
-		
-		print class_+'.'+info+' : ',eval(class_+'.'+info+'.__doc__')
+        if type(class_)!=str:
+                class_=class_.__class__.__name__
+        for info in dir(eval(class_)):
+                
+                print class_+'.'+info+' : ',eval(class_+'.'+info+'.__doc__')
 
 
 

=== modified file 'Template/MadWeight/Python/clean.py'
--- Template/MadWeight/Python/clean.py	2011-10-25 14:42:02 +0000
+++ Template/MadWeight/Python/clean.py	2015-10-20 14:02:27 +0000
@@ -70,7 +70,7 @@
     ls=os.listdir('.')
     for element in ls:
         if os.path.isdir(element) and element[0]=='P':
-	    status,mess=clean.suppress_dir(os.path.join(element,run_name))
+            status,mess=clean.suppress_dir(os.path.join(element,run_name))
             #if not status:
             #    print 'supress ',element,' failed:' 
             #    print mess
@@ -93,7 +93,7 @@
     ls=os.listdir('.')
     for element in ls:
         if os.path.isdir(element) and len(element)>4 and element[:4]=='MW_P':
-	    status,mess=clean.suppress_dir(os.path.join(element,run_name))
+            status,mess=clean.suppress_dir(os.path.join(element,run_name))
             #if not status:
             #    print 'supress ',element,' failed:' 
             #    print mess

=== modified file 'Template/NLO/Cards/run_card.dat'
--- Template/NLO/Cards/run_card.dat	2015-07-30 15:57:44 +0000
+++ Template/NLO/Cards/run_card.dat	2015-10-20 14:02:27 +0000
@@ -66,6 +66,8 @@
 # WARNING: PYTHIA6PT works only for processes without FSR!!!!          *
 #***********************************************************************
   %(parton_shower)s   = parton_shower
+  %(shower_scale_factor)s       = shower_scale_factor ! multiply default shower starting
+                                  ! scale by this factor
 #***********************************************************************
 # Renormalization and factorization scales                             *
 # (Default functional form for the non-fixed scales is the sum of      *

=== modified file 'Template/NLO/Cards/shower_card.dat'
--- Template/NLO/Cards/shower_card.dat	2015-03-09 17:19:46 +0000
+++ Template/NLO/Cards/shower_card.dat	2015-10-20 14:02:27 +0000
@@ -85,15 +85,15 @@
 #                                                                      *
 # Examples of syntax:                                                  *
 # Z -> e+ e- or mu+ mu- with BR = 0.5 each                             *
-# DM_1 = 23 > -11 11 @ 0.5d0 @ 100                                     *
-# DM_2 = 23 > -13 13 @ 0.5d0 @ 100                                     *
+# DM_1 = 23 > -11 11 @ 0.5d0 @ 100
+# DM_2 = 23 > -13 13 @ 0.5d0 @ 100
 # H -> tau+ tau- with BR = 1                                           *
-# DM_3 = 25 > -15 15 @ 1.0d0 @ 0                                       *
+# DM_3 = 25 > -15 15 @ 1.0d0 @ 0
 # t -> nu_e e+ b with BR = 1 (HERWIG)                                  *
-# DM_4 = 6 > 12 -11 5 @ 1d0 @ 100                                      *
+# DM_4 = 6 > 12 -11 5 @ 1d0 @ 100
 # t -> nu_e e+ b with BR = 1 (PYTHIA)                                  *
-# DM_5 = 6 > 24 5 @ 1d0 @ 100                                          *
-# DM_6 = 24 > 12 -11 @ 1d0 @ 100                                       *
+# DM_5 = 6 > 24 5 @ 1d0 @ 100
+# DM_6 = 24 > 12 -11 @ 1d0 @ 100
 #***********************************************************************
 # Extra Libraries/analyses                                             *
 # The following lines need to be changed if the user does not want to  *

=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script	2015-04-24 09:57:53 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script	2015-10-20 14:02:27 +0000
@@ -517,6 +517,7 @@
 TimeShower:MEafterFirst = off      ! No Matrix-element corrections after first emission
 TimeShower:phiPolAsym = on         ! Azimuthal asymmetry induced by gluon polarization
 TimeShower:alphaSuseCMW = false    ! Use the CMW prescription in FSR
+TimeShower:weightGluonToQuark = 1  ! Use normal Altarelli-Parisi kernels for g -> q qbar
 
 ! 6) Initial-state shower.
 SpaceShower:pTmaxMatch = 1         ! Use scalup (re-check)

=== modified file 'Template/NLO/Source/run.inc'
--- Template/NLO/Source/run.inc	2015-04-02 22:56:24 +0000
+++ Template/NLO/Source/run.inc	2015-10-20 14:02:27 +0000
@@ -39,6 +39,9 @@
       common/cscales_current_values/muR2_current,muF12_current,
      #                              muF22_current,QES2_current
 
+c Re-scale of the shower starting scale
+      double precision shower_scale_factor
+      common/cshower_scale_factor/shower_scale_factor
 c
 c     Collider
 c

=== modified file 'Template/NLO/Source/run_config.inc'
--- Template/NLO/Source/run_config.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/run_config.inc	2015-10-20 14:02:27 +0000
@@ -34,7 +34,7 @@
       parameter     (PBS_QUE = 'madgraph')
 
       integer    ChanPerJob
-      parameter (ChanPerJob=1)  !Number of channels / job for survey
+      parameter (ChanPerJob=100000000)  !Number of channels / job for survey
 
 c      integer    max_np
 c      parameter (max_np=1)  !Number of channels / job for refine

=== modified file 'Template/NLO/SubProcesses/ajob_template'
--- Template/NLO/SubProcesses/ajob_template	2015-03-18 14:50:37 +0000
+++ Template/NLO/SubProcesses/ajob_template	2015-10-20 14:02:27 +0000
@@ -26,57 +26,38 @@
     exit
 fi
 
+channel=$1
+run_mode=$2
+runnumber=$3
+integration_step=$4
+
 TAGTAGTAGTAGTAGTAGTAG for i in 1 ; do
     
-    runnumber=0
-    if [[ $1 == '0' ]]; then
-        j=$2\_G$i
-        if [[ ! -e $j ]]; then
-            mkdir $j
-        fi
-        cd $j
-        if [[ "$4" != "" ]]; then
-            if [[ -e ../$4\_G$i ]]; then
-                if [[ $1 == '0' ]]; then
-                    cp -f ../$4\_G$i/mint_grids . >/dev/null 2>&1
-                    cp -f ../$4\_G$i/grid.MC_integer . >/dev/null 2>&1
-                elif [[ $1 == '1' ]]; then
-                    cp -f ../$4\_G$i/mint_grids . >/dev/null 2>&1
-                    cp -f ../$4\_G$i/grid.MC_integer . >/dev/null 2>&1
-                fi
-            else
-                echo "Cannot find directory ../$4\_G$i/" > log.txt
-                exit
-            fi
-    	fi
-    elif [[ $1 == '2' ]]; then
-        j=G$2$i
-        if [[ ! -e $j ]]; then
-            mkdir $j
-        fi
-        
-        cd $j
-        if [[ "$4" != "" ]]; then
-            if [[ "$4" == "H" ||"$4" == "S" || "$4" == "V" || "$4" == "B" || "$4" == "F" ]]; then
-                if [[ -e ../G$4$i ]]; then
-                    cp -f ../G$4$i/mint_grids ./preset_mint_grids >/dev/null 2>&1
-                    cp -f ../G$4$i/grid.MC_integer . >/dev/null 2>&1
-                else
-                    echo "Cannot find direcotry ../G$4$i/" > log.txt
-                    exit
-                fi
-            else
-                runnumber=$4
-                if [[ ! -e ../${j}_$4 ]]; then
-                    mkdir ../${j}_$4
-                fi
-                cd ../${j}_$4
-                ln -sf ../${j}/mint_grids 
-                ln -sf ../${j}/mint_grids_NLO 
-                ln -sf ../${j}/grid.MC_integer
-                ln -sf ../${j}/res_1
-            fi
-        fi
+    if [[ $run_mode == 'all' || $run_mode == 'born' ]] ; then
+	j=$run_mode\_G$i
+    else
+	if [[ $runnumber == '0' ]] ; then
+	    j=G$run_mode$i
+	else
+	    j=G$run_mode$i\_$runnumber
+	fi
+    fi
+    cd $j
+    
+    if  [[ -e res.dat ]] ; then
+	rm -f res.dat
+    fi
+    if  [[ -e log.txt ]] ; then
+	rm -f log.txt
+    fi
+    if  [[ -e MADatNLO.top ]] ; then
+	rm -f MADatNLO.top
+    fi
+    if  [[ -e MADatNLO.HwU ]] ; then
+	rm -f MADatNLO.HwU
+    fi
+    if  [[ -e MADatNLO.root ]] ; then
+	rm -f MADatNLO.root
     fi
     if  [[ -e randinit ]] ; then
 	rm -f randinit
@@ -97,67 +78,25 @@
     link1up FKS_params.dat
     link1up configs_and_props_info.dat
     link1up leshouche_info.dat
-#   Not necessary anymore
-#    link1up MadLoop5_resources
     link1up OLE_order.olc
     link1up param_card.dat
     link1up initial_states_map.dat
     
-#   check where is the basic file for the creation of input_app.txt
-#       
-    if [[ $1 == '0' ]]; then
-	if [[ $3 == '-1' && -e ./madinM1 ]] ; then
-	    input_template=./madinM1
-	else
- 	    if  [[  -e ../madin.$2 ]] ; then
- 		input_template=../madin.$2
- 	    else
- 		input_template=../../madin.$2
- 	    fi 
-	fi
-    elif [[ $1 == '2' ]]; then
-    	if [[ $3 == '0' || $3 == '2' ]]; then
-    	    if [[  -e ../madinMMC_$2.2 ]] ; then
- 		input_template=../madinMMC_$2.2
- 	    else
- 		input_template=../../madinMMC_$2.2
- 	    fi
- 	else
-	    input_template=./madinM1
-    	fi
-    fi
     
-    if [[ $1 == '0' ]]; then
-        head -n 5 $input_template >& input_app.txt
-        echo $i >> input_app.txt
-        tail -n 4 $input_template >> input_app.txt
-        T="$(date +%s)"
+    T="$(date +%s)"
+    if [[ $run_mode == 'all' || $run_mode == 'born' ]]; then
         ../madevent_mintFO > log.txt <input_app.txt 2>&1
-        status=$?
-        T="$(($(date +%s)-T))"
-        echo "Time in seconds: ${T}" >>log.txt
-    elif [[ $1 == '2' ]]; then
+    else
         if [[ $runnumber != 0 ]]; then
-            tar --extract --file=../nevents.tar nevts_${j}_$runnumber
-            mv  nevts_${j}_$runnumber nevts
             echo "$runnumber" >& moffset.dat
         fi
-        if [[ $3 == '0' || $3 == '2' ]]; then
-            head -n 6 $input_template > input_app.txt
-            echo $i >> input_app.txt
-            tail -n 3 $input_template >> input_app.txt
-        elif [[ $3 == '1' ]]; then
-            head -n 6 $input_template > input_app.txt
-            echo $i >> input_app.txt
-            tail -n 3 $input_template >> input_app.txt
-        fi
-        T="$(date +%s)"
         ../madevent_mintMC > log.txt <input_app.txt 2>&1
-        status=$?
-        T="$(($(date +%s)-T))"
-        echo "Time in seconds: ${T}" >>log.txt
-        cp -f log.txt log_MINT$3.txt  >/dev/null 2>&1
     fi
+    status=$?
+    T="$(($(date +%s)-T))"
+    echo "Time in seconds: ${T}" >>log.txt
+    cp -f log.txt log_MINT$integration_step.txt  >/dev/null 2>&1
+    cp -f res.dat res_$integration_step.dat  >/dev/null 2>&1
     exit $status
 done
     

=== removed file 'Template/NLO/SubProcesses/combine_results.sh'
--- Template/NLO/SubProcesses/combine_results.sh	2013-06-28 13:57:17 +0000
+++ Template/NLO/SubProcesses/combine_results.sh	1970-01-01 00:00:00 +0000
@@ -1,55 +0,0 @@
-#!/bin/bash
-
-# find the correct directory
-if [[ ! -d ./SubProcesses ]]; then
-    cd ../
-fi
-if [[ -d ./SubProcesses ]]; then
-    cd SubProcesses
-fi
-
-if [[ -e res.txt ]]; then
-    rm -f res.txt
-fi
-if [[ -e dirs.txt ]]; then
-    rm -f dirs.txt
-fi
-if [[ -e nevents_unweighted ]]; then
-    rm -f nevents_unweighted
-fi
-
-arg1=$1
-arg2=$2
-arg3=$3
-# shift the list of arguments by 3
-shift
-shift
-shift
-if [[ "$@" == "" ]]; then
-    echo "Please give the G directories that should be combined,"
-    echo "e.g. 'GF* GV*', as final arguments of this script" 
-    exit
-fi
-
-touch res.txt
-touch dirs.txt
-NTOT=0
-for dir in "$@" ; do
-    N=`ls -d P*/$dir | wc -l`
-    NTOT=`expr $NTOT + $N`
-    ls -d P*/$dir >> dirs.txt
-    grep -H 'Final result' P*/$dir/res_$arg1 >> res.txt
-done
-echo N of directories: $NTOT
-if [[ $arg1 == '0' ]] ; then
-    echo 'Determining the number of unweighted events per channel'
-elif [[ $arg1 == '1' ]] ; then
-    echo 'Updating the number of unweighted events per channel'
-fi
-./sumres.py $NTOT $arg2 $arg3
-
-echo 'Integrated abs(cross-section)' 
-tail -n2 res.txt | head -n1
-echo 'Integrated cross-section' 
-tail -n1 res.txt
-mv res.txt res_$arg1.txt

=== removed file 'Template/NLO/SubProcesses/combine_results_FO.sh'
--- Template/NLO/SubProcesses/combine_results_FO.sh	2013-06-28 13:57:17 +0000
+++ Template/NLO/SubProcesses/combine_results_FO.sh	1970-01-01 00:00:00 +0000
@@ -1,52 +0,0 @@
-#!/bin/bash
-
-# find the correct directory
-if [[ ! -d ./SubProcesses ]]; then
-    cd ../
-fi
-if [[ -d ./SubProcesses ]]; then
-    cd SubProcesses
-fi
-
-if [[ $1 == "0" ]] ; then
-    mint_mode=0
-    shift
-elif [[ $1 == "1" ]] ; then
-    mint_mode=1
-    shift
-elif [[ $1 == "2" ]] ; then
-    echo "Cannot combine results for mint_mode 2"
-    exit
-else
-    mint_mode=0
-fi
-
-if [[ -e res.txt ]]; then
-    rm -f res.txt
-fi
-if [[ -e dirs.txt ]]; then
-    rm -f dirs.txt
-fi
-
-req_acc=$1
-shift
-
-touch res.txt
-touch dirs.txt
-NTOT=0
-for dir in "$@" ; do
-    N=`ls -d P*/$dir | wc -l`
-    NTOT=`expr $NTOT + $N`
-    ls -d P*/$dir >> dirs.txt
-    grep -H 'Final result' P*/$dir/res_$mint_mode >> res.txt
-done
-
-sed -i.bak s/"\+\/\-"/" \+\/\-"/ res.txt
-
-echo N of directories: $NTOT
-
-./sumres.py $NTOT -1 $req_acc
-
-rm -r res.txt.bak
-
-tail -n1 res.txt

=== modified file 'Template/NLO/SubProcesses/combine_root.sh'
--- Template/NLO/SubProcesses/combine_root.sh	2013-12-09 10:15:28 +0000
+++ Template/NLO/SubProcesses/combine_root.sh	2015-10-20 14:02:27 +0000
@@ -65,7 +65,7 @@
 echo ".x combine_root.C" >> rootinput.txt
 echo ".q" >> rootinput.txt
 
-root < rootinput.txt
+root -b < rootinput.txt
 
 rm -f rootinput.txt
 }

=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2015-03-09 18:27:17 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2015-10-20 14:02:27 +0000
@@ -18,7 +18,6 @@
 C     LOCAL
 C
       integer i,j,l,l1,l2,ndim
-      integer npoints
       character*130 buf
 c
 c     Global
@@ -171,7 +170,7 @@
       else
         flat_grid=.false.
       endif
-      ndim = 3*(nexternal-2)-4
+      ndim = 3*(nexternal-nincoming)-4
       if (abs(lpp(1)) .ge. 1) ndim=ndim+1
       if (abs(lpp(2)) .ge. 1) ndim=ndim+1
 c Don't proceed if muF1#muF2 (we need to work out the relevant formulae
@@ -219,23 +218,7 @@
             do j=1,nintervals_virt
                read (12,*) (ave_virt(j,i),i=1,ndim)
             enddo
-            if (ncall.gt.0 .and. accuracy.ne.0d0) then
-               read (12,*) ans(1),unc(1),ncall,itmax
-c Update the number of PS points based on unc(1), ncall and accuracy
-               itmax_fl=itmax*(unc(1)/accuracy)**2
-               if (itmax_fl.le.4d0) then
-                  itmax=max(nint(itmax_fl),2)
-               elseif (itmax_fl.gt.4d0 .and. itmax_fl.le.16d0) then
-                  ncall=nint(ncall*itmax_fl/4d0)
-                  itmax=4
-               else
-                  itmax=nint(sqrt(itmax_fl))
-                  ncall=nint(ncall*itmax_fl/nint(sqrt(itmax_fl)))
-               endif
-               accuracy=accuracy/ans(1) ! relative accuracy on the ABS X-section
-            else
-               read (12,*) ans(1),unc(1),dummy,dummy
-            endif
+            read (12,*) ans(1),unc(1),dummy,dummy
             read (12,*) virtual_fraction,average_virtual
             close (12)
             write (*,*) "Update iterations and points to",itmax,ncall
@@ -253,10 +236,6 @@
          call mint(sigint,ndim,ncall,itmax,imode,xgrid,ymax,ymax_virt
      $        ,ans,unc,chi2)
          call topout
-         open(unit=58,file='res_0',status='unknown')
-         write(58,*)'Final result [ABS]:',ans(1),' +/-',unc(1)
-         write(58,*)'Final result:',ans(2),' +/-',unc(2)
-         close(58)
          write(*,*)'Final result [ABS]:',ans(1),' +/-',unc(1)
          write(*,*)'Final result:',ans(2),' +/-',unc(2)
          write(*,*)'chi**2 per D.o.F.:',chi2(1)
@@ -337,6 +316,11 @@
       write(*,*) 'Time spent in Other_tasks : ',tOther
       write(*,*) 'Time spent in Total : ',tTot
 
+      open (unit=12, file='res.dat',status='unknown')
+      write (12,*)ans(1),unc(1),ans(2),unc(2),itmax,ncall,tTot
+      close(12)
+
+
       if(i_momcmp_count.ne.0)then
         write(*,*)'     '
         write(*,*)'WARNING: genps_fks code 555555'
@@ -680,8 +664,6 @@
       character * 70 idstring
       logical savegrid
 
-      character * 80 runstr
-      common/runstr/runstr
       logical usexinteg,mint
       common/cusexinteg/usexinteg,mint
       logical unwgt
@@ -692,84 +674,98 @@
       double precision volh
       common/mc_int2/volh,mc_hel,ihel,fillh
 
-
+      logical done
+      character*100 buffer
 c-----
 c  Begin Code
 c-----
       mint=.true.
       unwgt=.false.
-      write(*,'(a)') 'Enter number of events and iterations: '
-      read(*,*) ncall,itmax
-      write(*,*) 'Number of events and iterations ',ncall,itmax
-      write(*,'(a)') 'Enter desired accuracy: '
-      read(*,*) accuracy
-      write(*,*) 'Desired absolute accuracy: ',accuracy
-
-      write(*,'(a)') 'Enter 0 for fixed, 2 for adjustable grid: '
-      read(*,*) use_cut
-      if (use_cut .lt. 0 .or. use_cut .gt. 2) then
-         write(*,*) 'Bad choice, using 2',use_cut
-         use_cut = 2
-      endif
-
-      write(*,10) 'Suppress amplitude (0 no, 1 yes)? '
-      read(*,*) i
-      if (i .eq. 1) then
-         multi_channel = .true.
-         write(*,*) 'Using suppressed amplitude.'
-      else
-         multi_channel = .false.
-         write(*,*) 'Using full amplitude.'
-      endif
-
-      write(*,10) 'Exact helicity sum (0 yes, n = number/event)? '
-      read(*,*) i
-      if (i .eq. 0) then
-         mc_hel = 0
-         write(*,*) 'Explicitly summing over helicities for virt'
-      else
-         mc_hel= i
-         write(*,*) 'Summing over',i,' helicities/event for virt'
-      endif
-      isum_hel=0
-
-      write(*,10) 'Enter Configuration Number: '
-      read(*,*) dconfig
-      iconfig = int(dconfig)
-      do i=1,mapconfig(0)
-         if (iconfig.eq.mapconfig(i)) then
-            iconfig=i
-            exit
+      open (unit=83,file='input_app.txt',status='old')
+      done=.false.
+      do while (.not. done)
+         read(83,'(a)',err=222,end=222) buffer
+         if (buffer(1:7).eq.'NPOINTS') then
+            buffer=buffer(10:100)
+            read(buffer,*) ncall
+            write (*,*) 'Number of phase-space points per iteration:',ncall
+         elseif(buffer(1:11).eq.'NITERATIONS') then
+            read(buffer(14:),*) itmax
+            write (*,*) 'Maximum number of iterations is:',itmax
+         elseif(buffer(1:8).eq.'ACCURACY') then
+            read(buffer(11:),*) accuracy
+            write (*,*) 'Desired accuracy is:',accuracy
+         elseif(buffer(1:10).eq.'ADAPT_GRID') then
+            read(buffer(13:),*) use_cut
+            write (*,*) 'Using adaptive grids:',use_cut
+         elseif(buffer(1:12).eq.'MULTICHANNEL') then
+            read(buffer(15:),*) i
+            if (i.eq.1) then
+               multi_channel=.true.
+               write (*,*) 'Using Multi-channel integration'
+            else
+               multi_channel=.false.
+               write (*,*) 'Not using Multi-channel integration'
+            endif
+         elseif(buffer(1:12).eq.'SUM_HELICITY') then
+            read(buffer(15:),*) i
+            if (nincoming.eq.1) then
+               write (*,*) 'Sum over helicities in the virtuals'/
+     $              /' for decay process'
+               mc_hel=0
+            elseif (i.eq.0) then
+               mc_hel=0
+               write (*,*) 'Explicitly summing over helicities'/
+     $              /' for the virtuals'
+            else
+               mc_hel=1
+               write(*,*) 'Do MC over helicities for the virtuals'
+            endif
+            isum_hel=0
+         elseif(buffer(1:7).eq.'CHANNEL') then
+            read(buffer(10:),*) dconfig
+            iconfig = int(dconfig)
+            do i=1,mapconfig(0)
+               if (iconfig.eq.mapconfig(i)) then
+                  iconfig=i
+                  exit
+               endif
+            enddo
+            write(*,12) 'Running Configuration Number: ',iconfig
+         elseif(buffer(1:5).eq.'SPLIT') then
+            read(buffer(8:),*) i
+            write (*,*) 'Splitting channel:',i
+         elseif(buffer(1:8).eq.'RUN_MODE') then
+            read(buffer(11:),*) abrvinput
+            if(abrvinput(5:5).eq.'0')then
+               nbody=.true.
+            else
+               nbody=.false.
+            endif
+            abrv=abrvinput(1:4)
+            write (*,*) "doing the ",abrv," of this channel"
+            if(nbody)then
+               write (*,*) "integration Born/virtual with Sfunction=1"
+            else
+               write (*,*) "Normal integration (Sfunction != 1)"
+            endif
+         elseif(buffer(1:7).eq.'RESTART') then
+            read(buffer(10:),*) irestart
+            if (irestart.eq.0) then
+               write (*,*) 'RESTART: Fresh run'
+            elseif(irestart.eq.-1) then
+               write (*,*) 'RESTART: Use old grids, but refil plots'
+            elseif(irestart.eq.1) then
+               write (*,*) 'RESTART: continue with existing run'
+            else
+               write (*,*) 'RESTART:',irestart
+            endif
          endif
+         cycle
+ 222     done=.true.
       enddo
-      write(*,12) 'Running Configuration Number: ',iconfig
-c
-c Enter parameters that control Vegas grids
-c
-      write(*,*)'enter id string for this run'
-      read(*,*) idstring
-      runstr=idstring
-      write(*,*)'enter 1 if you want restart files'
-      read (*,*) itmp
-      if(itmp.eq.1) then
-         savegrid = .true.
-      else
-         savegrid = .false.
-      endif
-      write(*,*)'enter 0 to exclude, 1 for new run, 2 to restart'
-      read(5,*)irestart
+      close(83)
 
-      abrvinput='     '
-      write (*,*) "'all ', 'born', 'real', 'virt', 'novi' or 'grid'?"
-      write (*,*) "Enter 'born0' or 'virt0' to perform"
-      write (*,*) " a pure n-body integration (no S functions)"
-      read(5,*) abrvinput
-      if(abrvinput(5:5).eq.'0')then
-        nbody=.true.
-      else
-        nbody=.false.
-      endif
-      abrv=abrvinput(1:4)
       if (fks_configs.eq.1) then
          if (pdg_type_d(1,fks_i_d(1)).eq.-21) then
             write (*,*) 'Process generated with [LOonly=QCD]. '/
@@ -782,26 +778,6 @@
             endif
          endif
       endif
-c Options are way too many: make sure we understand all of them
-      if ( abrv.ne.'all '.and.abrv.ne.'born'.and.abrv.ne.'real'.and.
-     &     abrv.ne.'virt'.and.
-     &     abrv.ne.'viSC'.and.abrv.ne.'viLC'.and.abrv.ne.'novA'.and.
-     &     abrv.ne.'novB'.and.abrv.ne.'viSA'.and.abrv.ne.'viSB') then
-        write(*,*)'Error in input: abrv is:',abrv
-        stop
-      endif
-      if(nbody.and.abrv.ne.'born'.and.abrv(1:2).ne.'vi'
-     &     .and. abrv.ne.'grid')then
-        write(*,*)'Error in driver: inconsistent input',abrvinput
-        stop
-      endif
-
-      write (*,*) "doing the ",abrv," of this channel"
-      if(nbody)then
-        write (*,*) "integration Born/virtual with Sfunction=1"
-      else
-        write (*,*) "Normal integration (Sfunction != 1)"
-      endif
 c
 c
 c     Here I want to set up with B.W. we map and which we don't

=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2015-08-13 12:43:02 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2015-10-20 14:02:27 +0000
@@ -108,6 +108,11 @@
       call cpu_time(tBefore)
       fixed_order=.false.
       nlo_ps=.true.
+      if (nincoming.ne.2) then
+         write (*,*) 'Decay processes not supported for'/
+     &        /' event generation'
+         stop 1
+      endif
 
 c     Read general MadFKS parameters
 c
@@ -158,7 +163,7 @@
       else
         flat_grid=.false.
       endif
-      ndim = 3*(nexternal-2)-4
+      ndim = 3*(nexternal-nincoming)-4
       if (abs(lpp(1)) .ge. 1) ndim=ndim+1
       if (abs(lpp(2)) .ge. 1) ndim=ndim+1
 c Don''t proceed if muF1#muF2 (we need to work out the relevant formulae
@@ -189,7 +194,7 @@
             enddo
          else
 c to restore grids:
-            open (unit=12, file='preset_mint_grids',status='old')
+            open (unit=12, file='mint_grids',status='old')
             do j=0,nintervals
                read (12,*) (xgrid(j,i),i=1,ndim)
             enddo
@@ -285,7 +290,7 @@
          close(58)
 
 c to save grids:
-         open (unit=12, file='mint_grids_NLO',status='unknown')
+         open (unit=12, file='mint_grids',status='unknown')
          write (12,*) (xgrid(0,i),i=1,ndim)
          do j=1,nintervals
             write (12,*) (xgrid(j,i),i=1,ndim)
@@ -301,7 +306,6 @@
          write (12,*) virtual_fraction,average_virtual
          close (12)
 
-
 c*************************************************************
 c     event generation
 c*************************************************************
@@ -323,7 +327,7 @@
          ncall=nevts ! Update ncall with the number found in 'nevts'
 
 c to restore grids:
-         open (unit=12, file='mint_grids_NLO',status='unknown')
+         open (unit=12, file='mint_grids',status='unknown')
          read (12,*) (xgrid(0,i),i=1,ndim)
          do j=1,nintervals
             read (12,*) (xgrid(j,i),i=1,ndim)
@@ -468,6 +472,15 @@
       write(*,*) 'Time spent in Other_tasks : ',tOther
       write(*,*) 'Time spent in Total : ',tTot
 
+      open (unit=12, file='res.dat',status='unknown')
+      if (imode.eq.0) then
+         write (12,*)ans(1),unc(1),ans(2),unc(2),itmax,ncall,tTot
+      else
+         write (12,*)ans(1)+ans(5),sqrt(unc(1)**2+unc(5)**2),ans(2)
+     $        ,unc(2),itmax,ncall,tTot
+      endif
+      close(12)
+
       return
  999  write (*,*) 'nevts file not found'
       stop
@@ -630,12 +643,17 @@
 
       write(*,10) 'Exact helicity sum (0 yes, n = number/event)? '
       read(*,*) i
-      if (i .eq. 0) then
-         mc_hel= 0
-         write(*,*) 'Explicitly summing over helicities for virt'
+      if (nincoming.eq.1) then
+         write (*,*) 'Sum over helicities in the virtuals'/
+     $        /' for decay process'
+         mc_hel=0
+      elseif (i.eq.0) then
+         mc_hel=0
+         write (*,*) 'Explicitly summing over helicities'/
+     $        /' for the virtuals'
       else
-         mc_hel= i
-         write(*,*) 'Summing over',i,' helicities/event for virt'
+         mc_hel=1
+         write(*,*) 'Do MC over helicities for the virtuals'
       endif
       isum_hel = 0
 

=== modified file 'Template/NLO/SubProcesses/fks_Sij.f'
--- Template/NLO/SubProcesses/fks_Sij.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/SubProcesses/fks_Sij.f	2015-10-20 14:02:27 +0000
@@ -90,9 +90,13 @@
 
 c Consistency check -- call to set_cms_stuff() must be done prior to
 c entering this function
-      shattmp=2d0*dot(p(0,1),p(0,2))
+      if (nincoming.eq.2) then
+         shattmp=2d0*dot(p(0,1),p(0,2))
+      else
+         shattmp=p(0,1)**2
+      endif
       if(abs(shattmp/shat-1.d0).gt.1.d-5)then
-        write(*,*)'Error in fks_Sij: inconsistent shat'
+        write(*,*)'Error in fks_Sij: inconsistent shat #1'
         write(*,*)shattmp,shat
         stop
       endif
@@ -148,6 +152,7 @@
            do j=1,fks_j_from_i(i,0)
             kk = i
             ll = fks_j_from_i(i,j)
+            if (nincoming.ne.2 .and. ll.le.nincoming) cycle
             if     ( ijskip(kk,ll).eq.0 .and. ijskip(ll,kk).eq.0 ) then
                ijskip(kk,ll) = 1
             elseif ( ijskip(kk,ll).eq.0 .and. ijskip(ll,kk).eq.1 ) then
@@ -178,6 +183,7 @@
         do j=1,fks_j_from_i(i,0)
          kk = i
          ll = fks_j_from_i(i,j)
+         if (nincoming.ne.2 .and. ll.le.nincoming) cycle
          if(ijskip(kk,ll).ne.1)goto 222
          if(particle_type(ll).eq.8.and.particle_type(kk).ne.8.and.
      #      ll.gt.nincoming)then
@@ -567,7 +573,11 @@
 
 c Consistency check -- call to set_cms_stuff() must be done prior to
 c entering this function
-      shattmp=2d0*dot(p(0,1),p(0,2))
+      if (nincoming.eq.2) then
+         shattmp=2d0*dot(p(0,1),p(0,2))
+      else
+         shattmp=p(0,1)**2
+      endif
       if(abs(shattmp/shat-1.d0).gt.1.d-5)then
         write(*,*)'Error in fks_Hij: inconsistent shat'
         write(*,*)shattmp,shat

=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2015-08-19 12:50:52 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2015-10-20 14:02:27 +0000
@@ -1034,6 +1034,10 @@
      $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
       common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
       if (wgt1.eq.0d0 .and. wgt2.eq.0d0 .and. wgt3.eq.0d0) return
+c Check for NaN's and INF's. Simply skip the contribution
+      if (wgt1.ne.wgt1) return
+      if (wgt2.ne.wgt2) return
+      if (wgt3.ne.wgt3) return
       icontr=icontr+1
       if (icontr.gt.max_contr) then
          write (*,*) 'ERROR in add_wgt: too many contributions'
@@ -1226,7 +1230,12 @@
 c PDG codes
       niproc(ict)=iproc
       do j=1,iproc
-         parton_iproc(j,ict)=pd(j)*conv
+         if (nincoming.eq.2) then
+            parton_iproc(j,ict)=pd(j)*conv
+         else
+c           Keep GeV's for decay processes (no conv. factor needed)
+            parton_iproc(j,ict)=pd(j)
+         endif
          do k=1,nexternal
             parton_pdg(k,j,ict)=idup_d(iFKS,k,j)
             if (k.lt.fks_j_d(iFKS)) then
@@ -2195,9 +2204,15 @@
             do ii=1,iproc_save(nFKS(ict))
                if (eto(ii,nFKS(ict)).ne.ipr) cycle
                n_ctr_found=n_ctr_found+1
-               write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
-     &              (wgt(j,ict)*conv,j=1,3),
-     &              nexternal
+               if (nincoming.eq.2) then
+                  write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
+     &                 (wgt(j,ict)*conv,j=1,3),
+     &                 nexternal
+               else
+                  write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
+     &                 (wgt(j,ict),j=1,3),
+     &                 nexternal
+               endif
                procid=''
                do j=1,nexternal
                   write (str_temp,*) parton_pdg(j,ii,ict)
@@ -2223,9 +2238,15 @@
 c H-event
             ipr=iproc_picked
             n_ctr_found=n_ctr_found+1
-            write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
-     &           (wgt(j,ict)*conv,j=1,3),
-     &           nexternal
+            if (nincoming.eq.2) then
+               write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
+     &              (wgt(j,ict)*conv,j=1,3),
+     &              nexternal
+            else
+               write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
+     &              (wgt(j,ict),j=1,3),
+     &              nexternal
+            endif
             procid=''
             do j=1,nexternal
                write (str_temp,*) parton_pdg(j,ipr,ict)
@@ -2338,7 +2359,7 @@
       do i=0,3
         xsum(i)=0.d0
         xsuma(i)=0.d0
-        do j=3,npart
+        do j=nincoming+1,npart
           xsum(i)=xsum(i)+xmom(i,j)
           xsuma(i)=xsuma(i)+abs(xmom(i,j))
         enddo
@@ -2405,9 +2426,14 @@
       pass=.true.
       jflag=0
       do i=0,3
-        xsum(i)=-xmom(i,1)-xmom(i,2)
-        xsuma(i)=abs(xmom(i,1))+abs(xmom(i,2))
-        do j=3,npart
+        if (nincoming.eq.2) then
+          xsum(i)=-xmom(i,1)-xmom(i,2)
+          xsuma(i)=abs(xmom(i,1))+abs(xmom(i,2))
+        elseif(nincoming.eq.1) then
+          xsum(i)=-xmom(i,1)
+          xsuma(i)=abs(xmom(i,1))
+        endif
+        do j=nincoming+1,npart
           xsum(i)=xsum(i)+xmom(i,j)
           xsuma(i)=xsuma(i)+abs(xmom(i,j))
         enddo
@@ -2420,14 +2446,14 @@
           write(*,*)'Momentum is not conserved [nocms]'
           write(*,*)'i=',i
           do j=1,npart
-            write(*,'(4(d14.8,1x))') (xmom(jj,j),jj=0,3)
+            write(*,'(i2,1x,4(d14.8,1x))') j,(xmom(jj,j),jj=0,3)
           enddo
           jflag=1
         endif
       enddo
       if(jflag.eq.1)then
-        write(*,'(4(d14.8,1x))') (xsum(jj),jj=0,3)
-        write(*,'(4(d14.8,1x))') (xrat(jj),jj=0,3)
+        write(*,'(a3,1x,4(d14.8,1x))') 'sum',(xsum(jj),jj=0,3)
+        write(*,'(a3,1x,4(d14.8,1x))') 'rat',(xrat(jj),jj=0,3)
         pass=.false.
         return
       endif
@@ -2454,7 +2480,11 @@
         endif
       enddo
 c
-      ecmtmp=sqrt(2d0*dot(xmom(0,1),xmom(0,2)))
+      if (nincoming.eq.2) then
+         ecmtmp=sqrt(2d0*dot(xmom(0,1),xmom(0,2)))
+      elseif (nincoming.eq.1) then
+         ecmtmp=xmom(0,1)
+      endif
       if(abs(ecm-ecmtmp).gt.vtiny)then
         write(*,*)'Inconsistent shat [nocms]'
         write(*,*)'ecm given=   ',ecm
@@ -2687,6 +2717,9 @@
 
       double precision zero,tiny
       parameter (zero=0d0)
+      
+      integer icount
+      data icount /0/
 
 c Particle types (=color) of i_fks, j_fks and fks_mother
       integer i_type,j_type,m_type
@@ -2709,7 +2742,11 @@
 
 c Consistency check -- call to set_cms_stuff() must be done prior to
 c entering this function
-      shattmp=2d0*dot(pp(0,1),pp(0,2))
+      if (nincoming.eq.2) then
+         shattmp=2d0*dot(pp(0,1),pp(0,2))
+      else
+         shattmp=pp(0,1)**2
+      endif
       if(abs(shattmp/shat-1.d0).gt.1.d-5)then
         write(*,*)'Error in sreal: inconsistent shat'
         write(*,*)shattmp,shat
@@ -2717,9 +2754,9 @@
       endif
 
       if (1d0-y_ij_fks.lt.tiny)then
-         if (pmass(j_fks).eq.zero.and.j_fks.le.2)then
+         if (pmass(j_fks).eq.zero.and.j_fks.le.nincoming)then
             call sborncol_isr(pp,xi_i_fks,y_ij_fks,wgt)
-         elseif (pmass(j_fks).eq.zero.and.j_fks.ge.3)then
+         elseif (pmass(j_fks).eq.zero.and.j_fks.ge.nincoming+1)then
             call sborncol_fsr(pp,xi_i_fks,y_ij_fks,wgt)
          else
             wgt=0d0
@@ -2741,11 +2778,20 @@
       endif
 
       if(wgt.lt.0.d0)then
-         write(*,*) 'Fatal error #2 in sreal',wgt,xi_i_fks,y_ij_fks
-         do i=1,nexternal
-            write(*,*) 'particle ',i,', ',(pp(j,i),j=0,3)
-         enddo
-         stop
+         icount=icount+1
+         if (icount.le.10) then
+            write(*,*) 'Warning, numerical problem found in sreal. '/
+     $           /'Setting weight to zero',wgt,xi_i_fks,y_ij_fks
+            do i=1,nexternal
+               write(*,*) 'particle ',i,', ',(pp(j,i),j=0,3)
+            enddo
+            if (icount.eq.25) then
+               write (*,*) 'ERROR 25 problems found... '/
+     $              /'stopping the code'
+               stop
+            endif
+         endif
+         wgt=0d0
       endif
 
       return
@@ -2814,7 +2860,7 @@
       E_i_fks = p(0,i_fks)
       z = 1d0 - E_i_fks/(E_i_fks+E_j_fks)
       t = z * shat/4d0
-      if(rotategranny .and. nexternal-1.ne.3)then
+      if(rotategranny .and. nexternal-1.ne.3 .and. nincoming.eq.2)then
 c Exclude 2->1 (at the Born level) processes: matrix elements are
 c independent of the PS point, but non-zero helicity configurations
 c might flip when rotating the momenta.
@@ -2934,7 +2980,7 @@
 c Thus, an extra factor z (implicit in the flux of the reduced Born
 c in FKS) has to be inserted here
       t = z*shat/4d0
-      if(j_fks.eq.2 .and. nexternal-1.ne.3)then
+      if(j_fks.eq.2 .and. nexternal-1.ne.3 .and. nincoming.eq.2)then
 c Rotation according to innerpin.m. Use rotate_invar() if a more 
 c general rotation is needed.
 c Exclude 2->1 (at the Born level) processes: matrix elements are
@@ -2965,7 +3011,7 @@
                pi(i)=p_i_fks_ev(i)
                pj(i)=p(i,j_fks)
             enddo
-            if(j_fks.eq.2)then
+            if(j_fks.eq.2 .and. nincoming.eq.2)then
 c Rotation according to innerpin.m. Use rotate_invar() if a more 
 c general rotation is needed
                pi(1)=-pi(1)
@@ -2986,7 +3032,7 @@
             azifact=Wij_angle/Wij_recta
          endif
 c Insert the extra factor due to Madgraph convention for polarization vectors
-         if(j_fks.eq.2)then
+         if(j_fks.eq.2 .and. nincoming.eq.2)then
            cphi_mother=-1.d0
            sphi_mother=0.d0
          else
@@ -3461,7 +3507,11 @@
 
 c Consistency check -- call to set_cms_stuff() must be done prior to
 c entering this function
-      shattmp=2d0*dot(p(0,1),p(0,2))
+      if (nincoming.eq.2) then
+         shattmp=2d0*dot(p(0,1),p(0,2))
+      else
+         shattmp=p(0,1)**2
+      endif
       if(abs(shattmp/shat-1.d0).gt.1.d-5)then
         write(*,*)'Error in sreal: inconsistent shat'
         write(*,*)shattmp,shat
@@ -3770,7 +3820,7 @@
         enddo
       enddo
       do i=0,3
-        if(j_fks.gt.2)then
+        if(j_fks.gt.nincoming)then
           xnum=p1_cnt(i,i_fks,inum)+p1_cnt(i,j_fks,inum)
           xden=p1_cnt(i,i_fks,iden)+p1_cnt(i,j_fks,iden)
         else
@@ -4370,7 +4420,11 @@
 
 c Consistency check -- call to set_cms_stuff() must be done prior to
 c entering this function
-         shattmp=2d0*dot(p(0,1),p(0,2))
+         if (nincoming.eq.2) then
+            shattmp=2d0*dot(p(0,1),p(0,2))
+         else
+            shattmp=p(0,1)**2
+         endif
          if(abs(shattmp/shat-1.d0).gt.1.d-5)then
            write(*,*)'Error in sreal: inconsistent shat'
            write(*,*)shattmp,shat
@@ -4447,7 +4501,7 @@
          enddo
 c
          do i=1,nincoming
-            if (particle_type(i).ne.1)then
+            if (particle_type(i).ne.1 .and. pmass(i).eq.ZERO) then
                if (particle_type(i).eq.8) then
                   aj=0
                elseif(abs(particle_type(i)).eq.3) then
@@ -4972,8 +5026,13 @@
           elseif(abrv.ne.'virt' .and. abrv.ne.'viSC' .and.
      #           abrv.ne.'viLC')then
 c 1+2+3+4
-            tmp=dlog(xicut_used**2*shat/QES2)-
-     #          1/betai*dlog((1+betai)/(1-betai))
+            if (betai.gt.1d-6) then
+               tmp=dlog(xicut_used**2*shat/QES2)-
+     &              1/betai*dlog((1+betai)/(1-betai))
+            else
+               tmp=dlog(xicut_used**2*shat/QES2)-
+     &              2d0*(1d0+betai**2/3d0+betai**4/5d0)
+            endif
           else
              write(*,*)'Error #14 in eikonal_Ireg',abrv
              stop
@@ -5200,7 +5259,11 @@
               elseif(pmass(m).ne.zero.and.pmass(n).ne.zero)then
                 kikj=dot(p(0,n),p(0,m))
                 vij=sqrt(1-(pmass(n)*pmass(m)/kikj)**2)
-                single=single+0.5d0*1/vij*log((1+vij)/(1-vij))*wgt
+                if (vij .gt. 1d-6) then
+                   single=single+0.5d0*1/vij*log((1+vij)/(1-vij))*wgt
+                else
+                   single=single+(1d0+vij**2/3d0+vij**4/5d0)*wgt
+                endif
               else
                 write(*,*)'Error in getpoles',i,j,n,m,pmass(n),pmass(m)
                 stop

=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
--- Template/NLO/SubProcesses/genps_fks.f	2015-02-04 13:26:19 +0000
+++ Template/NLO/SubProcesses/genps_fks.f	2015-10-20 14:02:27 +0000
@@ -182,15 +182,19 @@
          endif
       enddo
       if( firsttime .or. iconfig0.ne.iconfigsave ) then
-         stot = 4d0*ebeam(1)*ebeam(2)
+         if (nincoming.eq.2) then
+            stot = 4d0*ebeam(1)*ebeam(2)
+         else
+            stot=pmass(1)**2
+         endif
 c Make sure have enough mass for external particles
          totmassin=0d0
-         do i=3-nincoming,2
+         do i=1,nincoming
             totmassin=totmassin+m(i)
          enddo
          totmass=0d0
          nbranch = nexternal-3 ! nexternal is for n+1-body, while itree uses n-body
-         do i=3,nbranch+2
+         do i=nincoming+1,nexternal-1
             totmass=totmass+m(i)
          enddo
          fksmass=totmass
@@ -294,7 +298,7 @@
       else
          pb(0,1)=sqrtshat_born
          do i=1,2
-            pb(0,1)=0d0
+            pb(i,1)=0d0
          enddo
          p(3,1)=1e-14           ! For HELAS routine ixxxxx for neg. mass
       endif
@@ -453,22 +457,22 @@
             isolsign=1
             call generate_momenta_massless_final(icountevts,i_fks,j_fks
      &           ,p_born(0,imother),shat,sqrtshat ,x(ixEi),xmrec2,xp
-     &           ,phi_i_fks,   xiimax,xinorm,xi_i_fks,y_ij_fks,p_i_fks
+     &           ,phi_i_fks,xiimax,xinorm,xi_i_fks,y_ij_fks,p_i_fks
      &           ,xjac,xpswgt,pass)
             if (.not.pass) goto 112
          elseif(m_j_fks.gt.0d0) then
             call generate_momenta_massive_final(icountevts,isolsign
-     &           ,i_fks,j_fks,p_born(0,imother),shat,sqrtshat ,m_j_fks
-     &           ,x(ixEi),xmrec2,xp,phi_i_fks,    xiimax,xinorm
+     &           ,i_fks,j_fks,p_born(0,imother),shat,sqrtshat,m_j_fks
+     &           ,x(ixEi),xmrec2,xp,phi_i_fks,xiimax,xinorm
      &           ,xi_i_fks,y_ij_fks,p_i_fks,xjac,xpswgt,pass)
             if (.not.pass) goto 112
          endif
       elseif(j_fks.le.nincoming) then
          isolsign=1
          call generate_momenta_initial(icountevts,i_fks,j_fks,xbjrk_born
-     &        ,tau_born,ycm_born,ycmhat,shat_born,phi_i_fks ,xp ,x(ixEi)
-     &        ,shat,stot ,sqrtshat,tau,ycm,xbjrk ,p_i_fks,xiimax,xinorm
-     &        ,xi_i_fks,y_ij_fks ,xpswgt,xjac ,pass)
+     &        ,tau_born,ycm_born,ycmhat,shat_born,phi_i_fks,xp ,x(ixEi)
+     &        ,shat,stot,sqrtshat,tau,ycm,xbjrk,p_i_fks,xiimax,xinorm
+     &        ,xi_i_fks,y_ij_fks ,xpswgt,xjac,pass)
          if (.not.pass) goto 112
       else
          write (*,*) 'Error #2 in genps_fks.f',j_fks
@@ -750,6 +754,7 @@
 c remove the following if no importance sampling towards soft
 c singularity is performed when integrating over xi_i_hat
       xjac=xjac*2d0*x(1)
+
 c Check that xii is in the allowed range
       if( icountevts.eq.-100 .or. abs(icountevts).eq.1 )then
          if(xi_i_fks.gt.(1-xmrec2/shat))then
@@ -829,7 +834,11 @@
 c mother four momenta
       do i=0,3
          xp_mother(i)=xp(i,i_fks)+xp(i,j_fks)
-         recoil(i)=xp(i,1)+xp(i,2)-xp_mother(i)
+         if (nincoming.eq.2) then
+            recoil(i)=xp(i,1)+xp(i,2)-xp_mother(i)
+         else
+            recoil(i)=xp(i,1)-xp_mother(i)
+         endif
       enddo
       sumrec=recoil(0)+rho(recoil)
       sumrec2=sumrec**2
@@ -843,7 +852,7 @@
          xdir(j)=xp_mother(j)/x3len_fks_mother
       enddo
 c Perform the boost here
-      do i=3,nexternal
+      do i=nincoming+1,nexternal
          if(i.ne.i_fks.and.i.ne.j_fks.and.shybst.ne.0.d0)
      &      call boostwdir2(chybst,shybst,chybstmo,xdir,xp(0,i),xp(0,i))
       enddo
@@ -1145,7 +1154,11 @@
 c mother four momenta
       do i=0,3
          xp_mother(i)=xp(i,i_fks)+xp(i,j_fks)
-         recoil(i)=xp(i,1)+xp(i,2)-xp_mother(i)
+         if (nincoming.eq.2) then
+            recoil(i)=xp(i,1)+xp(i,2)-xp_mother(i)
+         else
+            recoil(i)=xp(i,1)-xp_mother(i)
+         endif
       enddo
 c
       sumrec=recoil(0)+rho(recoil)
@@ -1168,7 +1181,7 @@
          xdir(j)=xp_mother(j)/x3len_fks_mother
       enddo
 c Boost the momenta
-      do i=3,nexternal
+      do i=nincoming+1,nexternal
          if(i.ne.i_fks.and.i.ne.j_fks.and.shybst.ne.0.d0)
      &      call boostwdir2(chybst,shybst,chybstmo,xdir,xp(0,i),xp(0,i))
       enddo
@@ -2357,7 +2370,11 @@
       external dot
       pass=.true.
       do i=0,3
-         recoilbar(i)=p_born(i,1)+p_born(i,2)-p_born(i,imother)
+         if (nincoming.eq.2) then
+            recoilbar(i)=p_born(i,1)+p_born(i,2)-p_born(i,imother)
+         else
+            recoilbar(i)=p_born(i,1)-p_born(i,imother)
+         endif
       enddo
       xmrec2=dot(recoilbar,recoilbar)
       if(xmrec2.lt.0.d0)then

=== modified file 'Template/NLO/SubProcesses/iproc_map.f'
--- Template/NLO/SubProcesses/iproc_map.f	2015-03-09 18:27:17 +0000
+++ Template/NLO/SubProcesses/iproc_map.f	2015-10-20 14:02:27 +0000
@@ -211,7 +211,8 @@
       include "appl_common.inc"
 *
       character*200 buffer
-      integer procnum,found,i,l,j,ll
+      integer procnum,i,l,j,ll,found_a,found_m
+      logical found_appl(mxpdflumi),found_mg(maxproc)
       integer  nmatch_total
 
 *     npdflumi is the number of pdf luminosities in this particular process
@@ -296,7 +297,7 @@
      1           ((pdgs(i,j,kpdflumi),i=1,2),j=1,nproc(kpdflumi))
 *
             if(kpdflumi.le.0.or.kpdflumi.gt.mxpdflumi) then
-               write(6,*) "In fks_singular.f"
+               write(6,*) "In iproc_map.f"
                write(6,*) "Invalid value of kpdflumi = ",kpdflumi
                stop
             endif
@@ -341,33 +342,29 @@
 *     Be careful with all possible permutations in initial_parton_map!
 *     Check all possible npdflumi conbinations until a match is found
          do kpdflumi=1,npdflumi
-*     First, check that the same number of 
-*     parton subprocesses is used
-            if(maxproc_used.ne.nproc(kpdflumi)) continue
-            if(flav_map_debug) then
-               write(6,*) " "
-               write(6,*) " kpdflumi = ",kpdflumi
-               write(6,*) " nproc(kpdflumi) = ",nproc(kpdflumi)
-               do l=1,nproc(kpdflumi)
-                  write(6,*) l,pdgs(1,l,kpdflumi),pdgs(2,l,kpdflumi)
-               enddo
-               write(6,*) " "
-            endif
-            
+
 *     Initialization
-            found = 0
+            do l=1,nproc(kpdflumi)
+               found_appl(l) = .false.
+            enddo
+            do ll=1,maxproc_used
+               found_mg(ll) = .false.
+            enddo
+            found_a=0
+            found_m=0
 
 *     Look for all possible pairs of
 *     a) pdgs(1,l,kpdflumi),pdgs(2,l,kpdflumi) and
 *     b) IDUP(1,ll),  IDUP(2,ll)
 *     with kpdflumi fixed
             do l=1,nproc(kpdflumi)
-               do ll=1,nproc(kpdflumi)
+               do ll=1,maxproc_used
                   if ( ( pdgs(1,l,kpdflumi).eq.
      1                 IDUP(1,ll).and.pdgs(2,l,kpdflumi).
      2                 eq.IDUP(2,ll) ) ) 
      4                 then
-                     found = found + 1
+                     found_appl(l)=.true.
+                     found_mg(ll)=.true.
                      if(flav_map_debug) then
                         write(6,*) "match found!"
                         write(6,*) "pdgs = ",pdgs(1,l,kpdflumi),
@@ -377,8 +374,15 @@
                   endif
                enddo
             enddo
+            do l=1,nproc(kpdflumi)
+               if (found_appl(l)) found_a=found_a+1
+            enddo
+            do ll=1,maxproc_used
+               if (found_mg(ll)) found_m=found_m+1
+            enddo
 *
-            if (found.eq.nproc(kpdflumi)) then
+            if ( found_a.eq.nproc(kpdflumi) .and.
+     &           found_m.eq.maxproc_used ) then
                if(flav_map_debug) then
                   write(6,*) " ------------------------------- "
                   write(6,*) "         Match found!"
@@ -389,9 +393,9 @@
                nmatch_total = nmatch_total+1
             endif
          enddo
-         if(flavour_map(nFKSprocess).eq.0) then
+         if (nmatch_total.ne.nFKSprocess) then
             write(6,*) 
-     1           "Problem with setup_flavourmap in fks_singular.f"
+     1           "Problem with setup_flavourmap in iproc_map.f"
             write(6,*) "nFKSprocess = ",nFKSprocess
             write(6,*)" flavour_map(nFKSprocess)= ",
      1           flavour_map(nFKSprocess)
@@ -399,17 +403,6 @@
          endif
       enddo
       
-*     write(6,*) " "
-*     write(6,*) "  nmatch_total = ", nmatch_total
-*     write(6,*) " "
-      if( nmatch_total.ne.fks_configs) then
-         write(6,*) 
-     1        "Problem with setup_flavourmap  in fks_singular.f"
-         write(6,*) "nmatch_total, fks_configs"
-         write(6,*) nmatch_total, fks_configs
-         stop
-      endif
-      
 *     Check flavor map properly initialized
 *     All the entries of flavour_map(1:nFKSprocess) must be from
 *     1 to npdflumi

=== modified file 'Template/NLO/SubProcesses/leshouche_inc_chooser.f'
--- Template/NLO/SubProcesses/leshouche_inc_chooser.f	2015-03-09 18:27:17 +0000
+++ Template/NLO/SubProcesses/leshouche_inc_chooser.f	2015-10-20 14:02:27 +0000
@@ -91,8 +91,8 @@
                icolup_d(1,1,nexternal,j)=-99999 ! should not be used
                icolup_d(1,2,nexternal,j)=-99999
             enddo
+            return
          endif
-         return
       endif
 
       open(unit=78, file='leshouche_info.dat', status='old')

=== modified file 'Template/NLO/SubProcesses/montecarlocounter.f'
--- Template/NLO/SubProcesses/montecarlocounter.f	2014-10-06 16:45:08 +0000
+++ Template/NLO/SubProcesses/montecarlocounter.f	2015-10-20 14:02:27 +0000
@@ -2549,8 +2549,9 @@
       common/to_abrv/abrv
 
       ref_scale=sqrt((1-xi)*shat)
-      xscalemin=max(frac_low*ref_scale,scaleMClow)
-      xscalemax=max(frac_upp*ref_scale,xscalemin+scaleMCdelta)
+      xscalemin=max(shower_scale_factor*frac_low*ref_scale,scaleMClow)
+      xscalemax=max(shower_scale_factor*frac_upp*ref_scale,
+     &              xscalemin+scaleMCdelta)
       xscalemax=min(xscalemax,2d0*sqrt(ebeam(1)*ebeam(2)))
       xscalemin=min(xscalemin,xscalemax)
 c

=== modified file 'Template/NLO/SubProcesses/setcuts.f'
--- Template/NLO/SubProcesses/setcuts.f	2015-07-24 18:59:18 +0000
+++ Template/NLO/SubProcesses/setcuts.f	2015-10-20 14:02:27 +0000
@@ -91,8 +91,6 @@
             lpp(2)=0
             ebeam(1)=pmass(1)/2d0
             ebeam(2)=pmass(1)/2d0
-            scale=pmass(1)
-            fixed_ren_scale=.true.
          endif
 c-check consistency of maxjetflavor in the run_card and with Nf
 c specified in coupl.inc

=== modified file 'Template/NLO/SubProcesses/setscales.f'
--- Template/NLO/SubProcesses/setscales.f	2015-07-30 15:57:44 +0000
+++ Template/NLO/SubProcesses/setscales.f	2015-10-20 14:02:27 +0000
@@ -237,7 +237,10 @@
      $     ,FxFx_fac_scale
 c
       tmp=0
-      if(ickkw.eq.3)then
+      if (nincoming.eq.1) then
+         tmp=pp(0,1) ! mass of the decaying particle
+         temp_scale_id='Mass of decaying particle'
+      elseif(ickkw.eq.3)then
 c FxFx merging scale:
 c     Note that nFxFx_ren_scales includes the one scale that corresponds
 c     to the real-emission one (and is zero for the n-body conf.). Skip
@@ -411,7 +414,7 @@
       tmp=0
       if(ickkw.eq.3)then
 c FxFx merging scale:
-        tmp=min(FxFx_fac_scale(1),FxFx_fac_scale(2))
+        tmp=(FxFx_fac_scale(1)+FxFx_fac_scale(2))/2d0
         temp_scale_id='FxFx merging scale'
       elseif(imuftype.eq.1)then
         tmp=scale_global_reference(pp)
@@ -490,7 +493,10 @@
       parameter (iQEStype=1)
 c
       tmp=0
-      if(iQEStype.eq.1)then
+      if (nincoming.eq.1) then
+         tmp=pp(0,1) ! mass of the decaying particle
+         temp_scale_id='Mass of decaying particle'
+      elseif(iQEStype.eq.1)then
         tmp=scale_global_reference(pp)
       elseif(iQEStype.eq.2)then
         do i=nincoming+1,nexternal

=== removed file 'Template/NLO/SubProcesses/sumres.py'
--- Template/NLO/SubProcesses/sumres.py	2014-07-23 10:33:38 +0000
+++ Template/NLO/SubProcesses/sumres.py	1970-01-01 00:00:00 +0000
@@ -1,241 +0,0 @@
-#!/usr/bin/env python
-
-#script to combine reults and tell the number of events that need
-#  to be generated in each channel. 
-#  Replaces the sumres.f and sumres2.f files
-#  MZ, 2011-10-22
-
-from __future__ import division
-import math
-import sys
-import random
-import os
-
-nexpected=int(sys.argv[1])
-nevents=int(sys.argv[2])
-req_acc=float(sys.argv[3])
-# if nevents is >=0 the script will also determine the 
-# number of events required for each process
-
-
-def Mirrorprocs(p1, p2):
-    """determine if the folder names p1, p2 (with the _N already taken out)
-    correspond to the same process with
-    mirrror initial state. Returns true/false"""
-    return False
-
-file=open("res.txt")
-content = file.read()
-file.close()
-lines = content.split("\n")
-processes=[]
-tot=0
-err=0
-totABS=0
-errABS=0
-
-# open the file containing the list of directories
-file=open("dirs.txt")
-dirs = file.read().split("\n")
-file.close()
-dirs.remove('')
-
-# The syntax of lines should be first the ABS cross section for the
-# channel and the line after that the cross section for the same
-# channel.
-for line in range(0,len(lines),2):
-    list = lines[line].split()
-    if list:
-        proc={}
-        proc['folder'] = list[0].split('/')[0]
-        proc['subproc'] = proc['folder'][0:proc['folder'].rfind('_')]
-        proc['channel'] = list[0].split('/')[1]
-        dirs.remove(os.path.join(proc['folder'], proc['channel']))
-        proc['resultABS'] = float(list[4])
-        proc['errorABS'] = float(list[6])
-        proc['err_percABS'] = proc['errorABS']/proc['resultABS']*100.
-        processes.append(proc)
-        totABS+= proc['resultABS']
-        errABS+= math.pow(proc['errorABS'],2)
-        list = lines[line+1].split()
-        if list:
-            proc['result'] = float(list[3])
-            proc['error'] = float(list[5])
-            proc['err_perc'] = proc['error']/proc['result']*100.
-            tot+= proc['result']
-            err+= math.pow(proc['error'],2)
-if dirs:
-    print "%d jobs did not terminate correctly: " % len(dirs)
-    print '\n'.join(dirs)
-    print "The results are probably not correct. Please check the relevant log files corresponding to the above jobs for more information."
-
-processes.sort(key = lambda proc: -proc['errorABS'])
-
-correct = len(processes) == nexpected
-print "Found %d correctly terminated jobs " %len(processes) 
-if not len(processes)==nexpected:
-    print len(processes), nexpected
-
-subprocs_string=[]
-for proc in processes:
-    subprocs_string.append(proc['subproc'])
-subprocs_string=set(subprocs_string)
-
-content+='\n\nCross-section per integration channel:\n'
-for proc in processes:
-    content+='%(folder)20s  %(channel)15s   %(result)10.8e    %(error)6.4e       %(err_perc)6.4f%%  \n' %  proc
-
-content+='\n\nABS cross-section per integration channel:\n'
-for proc in processes:
-    content+='%(folder)20s  %(channel)15s   %(resultABS)10.8e    %(errorABS)6.4e       %(err_percABS)6.4f%%  \n' %  proc
-
-content+='\n\nCross-section per subprocess:\n'
-#for subpr in sorted(set(subprocs)):
-subprocesses=[]
-for sub in subprocs_string:
-    subpr={}
-    subpr['subproc']=sub
-    subpr['xsect']=0.
-    subpr['err']=0.
-    for proc in processes:
-        if proc['subproc'] == sub:
-            subpr['xsect'] += proc['result']
-            subpr['err'] += math.pow(proc['error'],2)
-    subpr['err']=math.sqrt(subpr['err'])
-    subprocesses.append(subpr)
-
-
-#find and combine mirror configurations (if in v4)
-for i1, s1 in enumerate(subprocesses):
-    for i2, s2 in enumerate(subprocesses):
-        if Mirrorprocs(s1['subproc'], s2['subproc']) and i1 >= i2:
-            s1['xsect'] += s2['xsect']
-            s1['err'] = math.sqrt(math.pow(s1['err'],2)+ math.pow(s2['err'],2))
-            s2['toremove'] = True
-
-new = []
-for s in subprocesses:
-    try:
-        a= s['toremove']
-    except KeyError:
-        new.append(s)
-subprocesses= new
-
-
-subprocesses.sort(key = lambda proc: -proc['xsect'])
-for subpr in subprocesses:
-    content+=  '%(subproc)20s    %(xsect)10.8e   %(err)6.4e\n' % subpr
- 
-
-content+='\nTotal ABS and \nTotal: \n                      %10.8e +- %6.4e  (%6.4e%%)\n                      %10.8e +- %6.4e  (%6.4e%%)\n' %\
-        (totABS, math.sqrt(errABS), math.sqrt(errABS)/totABS *100.,tot, math.sqrt(err), math.sqrt(err)/tot *100.)
-
-if not correct:
-    sys.exit('ERROR: not all jobs terminated correctly\n')
-
-file=open("res.txt", 'w')
-
-file.write(content)
-file.close()
-
-#determine the events for each process:
-if nevents>=0:
-    if req_acc<0: 
-        req_acc2_inv=nevents
-    else:
-        req_acc2_inv=1/(req_acc*req_acc)
-    #get the random number seed from the randinit file
-    file=open("randinit")
-    exec file
-    file.close
-    print "random seed found in 'randinit' is", r
-    random.seed(r)
-    totevts=nevents
-    for proc in processes:
-        proc['lhefile'] = os.path.join(proc['folder'], proc['channel'], 'events.lhe')
-        proc['nevents'] = 0
-    while totevts :
-        target = random.random() * totABS
-        crosssum = 0.
-        i = 0
-        while i<len(processes) and crosssum < target:
-            proc = processes[i]
-            crosssum += proc['resultABS']
-            i += 1            
-        totevts -= 1
-        i -= 1
-        processes[i]['nevents'] += 1
-        
-#check that we now have all the events in the channels
-    totevents = sum(proc['nevents'] for proc in processes)
-    if totevents != nevents:
-        sys.exit('failed to obtain the correct number of events. Required: %d,   Obtained: %d' \
-                 % (nevents, totevents))
-
-    content_evts = ''
-    for proc in processes:
-        content_evts+= ' '+proc['lhefile']+'              %(nevents)10d            %(resultABS)10.8e        1.0 \n' %  proc
-        nevts_file = open(os.path.join(proc['folder'], proc['channel'], 'nevts'),'w')
-        nevts_file.write('%10d\n' % proc['nevents'])
-        nevts_file.close()
-        if proc['channel'][1] == 'B':
-            fileinputs = open("madinMMC_B.2")
-        elif proc['channel'][1] == 'F':
-            fileinputs = open("madinMMC_F.2")
-        elif proc['channel'][1] == 'V':
-            fileinputs = open("madinMMC_V.2")
-        else:
-            sys.exit("ERROR, DONT KNOW WHICH INPUTS TO USE")
-        fileinputschannel = open(os.path.join(proc['folder'], proc['channel'], 'madinM1'),'w')
-        i=0
-        for line in fileinputs:
-            i += 1
-            if i == 2:
-                accuracy=min(math.sqrt(totABS/(req_acc2_inv*proc['resultABS'])),0.2)
-                fileinputschannel.write('%10.8e\n' % accuracy)
-            elif i == 8:
-                fileinputschannel.write('1        ! MINT mode\n')
-            else:
-                fileinputschannel.write(line)
-        fileinputschannel.close()
-        fileinputs.close()
-
-    evts_file = open('nevents_unweighted', 'w')
-    evts_file.write(content_evts)
-    evts_file.close()
-
-# if nevents = -1 and req_acc >= 0, we need to determine the required
-# accuracy in each of the channels: this is for fixed order running!
-elif req_acc>=0 and nevents==-1:
-    req_accABS=req_acc*abs(tot)/totABS
-    content_evts = ''
-    for proc in processes:
-        if proc['channel'][0:3] == 'all':
-            fileinputs = open("madin.all")
-        elif proc['channel'][0:4] == 'novB':
-            fileinputs = open("madin.novB")
-        elif proc['channel'][0:4] == 'born':
-            fileinputs = open("madin.born")
-        elif proc['channel'][0:4] == 'grid':
-            fileinputs = open("madin.grid")
-        elif proc['channel'][0:4] == 'viSB':
-            fileinputs = open("madin.viSB")
-        elif proc['channel'][0:4] == 'virt':
-            fileinputs = open("madin.virt")
-        elif proc['channel'][0:4] == 'novi':
-            fileinputs = open("madin.novi")
-        else:
-            sys.exit("ERROR, DONT KNOW WHICH INPUTS TO USE")
-        fileinputschannel = open(os.path.join(proc['folder'], proc['channel'], 'madinM1'),'w')
-        i=0
-        for line in fileinputs:
-            i += 1
-            if i == 2:
-                accuracy=req_accABS*math.sqrt(totABS*proc['resultABS'])
-                fileinputschannel.write('%10.8e\n' % accuracy)
-            elif i == 9:
-                fileinputschannel.write('-1        ! restart from existing grids\n')
-            else:
-                fileinputschannel.write(line)
-        fileinputschannel.close()
-        fileinputs.close()

=== modified file 'Template/NLO/SubProcesses/symmetry_fks_v3.f'
--- Template/NLO/SubProcesses/symmetry_fks_v3.f	2014-06-26 08:45:41 +0000
+++ Template/NLO/SubProcesses/symmetry_fks_v3.f	2015-10-20 14:02:27 +0000
@@ -541,6 +541,8 @@
       lname=4
       mname='mg'
       call open_bash_file(26,fname,lname,mname)
+      call close_bash_file(26)
+      open(unit=26,file='channels.txt',status='unknown')
       ic = 0      
       do i=1,mapconfig(0)
          if (use_config(i) .gt. 0) then
@@ -567,15 +569,6 @@
             done = .false.
             do while (.not. done)
                call enCode(icode,iarray,ibase,imax)
-               ic=ic+1
-               if (ic .gt. ChanPerJob) then
-                  call close_bash_file(26)
-                  fname='ajob'
-                  lname=4
-                  mname='mg'
-                  call open_bash_file(26,fname,lname,mname)
-                  ic = 1
-               endif
 c               write(*,*) 'mapping',ic,mapconfig(i)
 c$$$               if (r_from_b(mapconfig(i)) .lt. 10) then
 c$$$                  write(26,'(i1$)') r_from_b(mapconfig(i))
@@ -611,7 +604,7 @@
             enddo
          endif
       enddo
-      call close_bash_file(26)
+      close(26)
       if (mapconfig(0) .gt. 9999) then
          write(*,*) 'Only writing first 9999 jobs',mapconfig(0)
       endif

=== modified file 'Template/NLO/SubProcesses/write_ajob.f'
--- Template/NLO/SubProcesses/write_ajob.f	2012-10-18 06:17:30 +0000
+++ Template/NLO/SubProcesses/write_ajob.f	2015-10-20 14:02:27 +0000
@@ -53,7 +53,7 @@
          if (index(buff,'TAGTAGTAGTAGTAG').ne.0) exit
          write(lun,15) buff
       enddo
-      write(lun,'(a$)') 'for i in '
+      write(lun,'(a$)') 'for i in $channel '
       return
  99   write (*,*) 'ajob_template or ajob_template_cluster '/
      &     /'does not have the correct format'

=== modified file 'Template/loop_material/StandAlone/Cards/MadLoopParams.dat'
--- Template/loop_material/StandAlone/Cards/MadLoopParams.dat	2014-10-27 08:40:30 +0000
+++ Template/loop_material/StandAlone/Cards/MadLoopParams.dat	2015-10-20 14:02:27 +0000
@@ -4,8 +4,8 @@
 
 !     
 #MLReductionLib
-!1|4|3|2
-! Default :: 1|4|3|2
+!1|3|2
+! Default :: 1|3|2
 !     The tensor integral reduction library.The current choices are:
 !     1 | CutTools
 !     2 | PJFry++

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2015-09-07 18:05:30 +0000
+++ UpdateNotes.txt	2015-10-20 14:02:27 +0000
@@ -1,6 +1,29 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
-2.3.2.2 (06/09/15)
+<<<<<<< TREE
+2.3.2.2 (06/09/15)
+=======
+
+2.3.3(XX/XX/XX)
+        OM: Allow new syntax for the param_card: instead of an entry you can enter scan:[val1, val2,...]
+            To perform a scan on this parameter.
+        OM: Having two mode for "output pythia8" one (default) for pythia8.2 and one for pythia8.1 (with --version=8.1)
+	RF: Rewriting of job-control for NLO processes. Better accuracy estimates for FO processes
+	RF: Fix for factorisation scale setting in FxFx merging when very large difference in scale in the
+ 	    non-QCD part of a process. 
+	RF: Better discarding of numerical instabilities in the real-emission matrix elements. Only of interested for
+	    processes which have jets at Born level, but do not require generation cut (like t-channel single-top). 
+	RF: Added an option to the run_card to allow for easier variation of the shower starting scale (NLO only).
+	RF: Fixed a problem in the setting of the flavour map used for runs with iAPPL >= 1. 
+	RF: Allow for decay processes to compute (partial) decay widths at NLO accuracy (fixed order only).
+        OM: (SysCalc interface) Allow to bypass the pdf reweighting/alpsfact reweighting
+        MZ: fixed bug related to slurm clusters
+	OM: remove the addmasses.py script of running by default on gridpack mode. 
+            if you want to have it running, you just have to rename the file madevent/bin/internal/addmasses_optional.py to
+ 	    madevent/bin/internal/addmasses_optional.py and it will work as before. (Do not work with SysCalc tag)
+
+2.3.2.2 (06/09/15)
+>>>>>>> MERGE-SOURCE
         VH: Finalized the MG5aMC-GoSam interface
 
 2.3.2(20/08/15)   

=== modified file 'VERSION'
--- VERSION	2015-09-07 18:04:08 +0000
+++ VERSION	2015-10-20 14:02:27 +0000
@@ -1,5 +1,11 @@
+<<<<<<< TREE
 version = 2.3.2.2
 date = 2015-09-06
+=======
+version = 2.3.3
+date = 2015-10-20
+
+>>>>>>> MERGE-SOURCE
 
 
 

=== modified file 'aloha/__init__.py' (properties changed: -x to +x)
=== modified file 'aloha/aloha_fct.py' (properties changed: -x to +x)
=== modified file 'aloha/aloha_lib.py' (properties changed: -x to +x)
--- aloha/aloha_lib.py	2015-07-30 15:43:54 +0000
+++ aloha/aloha_lib.py	2015-10-20 14:02:27 +0000
@@ -1073,9 +1073,9 @@
                     self[(0,)] = representation[(0,)]
                 except Exception:
                     if representation:
-		   	            raise LorentzObjectRepresentation.LorentzObjectRepresentationError("There is no key of (0,) in representation.")
+                        raise LorentzObjectRepresentation.LorentzObjectRepresentationError("There is no key of (0,) in representation.")
                     else:
-			            self[(0,)] = 0
+                        self[(0,)] = 0
             else:
                 self[(0,)] = representation
 

=== modified file 'aloha/aloha_object.py' (properties changed: -x to +x)
--- aloha/aloha_object.py	2014-02-18 21:45:48 +0000
+++ aloha/aloha_object.py	2015-10-20 14:02:27 +0000
@@ -434,12 +434,12 @@
         self.sub11 = aloha_lib.Variable('T%s_6' % self.particle)
         self.sub12 = aloha_lib.Variable('T%s_7' % self.particle)
         self.sub13 = aloha_lib.Variable('T%s_8' % self.particle)
-	
+
         self.sub20 = aloha_lib.Variable('T%s_9' % self.particle)
         self.sub21 = aloha_lib.Variable('T%s_10' % self.particle)
         self.sub22 = aloha_lib.Variable('T%s_11' % self.particle)
         self.sub23 = aloha_lib.Variable('T%s_12' % self.particle)
-	
+
         self.sub30 = aloha_lib.Variable('T%s_13' % self.particle)
         self.sub31 = aloha_lib.Variable('T%s_14' % self.particle)
         self.sub32 = aloha_lib.Variable('T%s_15' % self.particle)

=== modified file 'aloha/aloha_parsers.py' (properties changed: -x to +x)
=== modified file 'aloha/aloha_writers.py' (properties changed: -x to +x)
--- aloha/aloha_writers.py	2014-12-11 15:23:43 +0000
+++ aloha/aloha_writers.py	2015-10-20 14:02:27 +0000
@@ -689,7 +689,7 @@
                     if number.imag == 1:
                         out = 'CI'
                     elif number.imag == -1:
-			            out = '-CI'
+                        out = '-CI'
                     else: 
                         out = '%s * CI' % self.change_number_format(number.imag)
             else:

=== modified file 'aloha/create_aloha.py' (properties changed: -x to +x)
--- aloha/create_aloha.py	2015-06-28 14:44:22 +0000
+++ aloha/create_aloha.py	2015-10-20 14:02:27 +0000
@@ -854,7 +854,7 @@
         # reorganize the data (in order to use optimization for a given lorentz
         #structure
         aloha.loop_mode = False
-	    # self.explicit_combine = False
+        # self.explicit_combine = False
         request = {}
 
         for list_l_name, tag, outgoing in data:

=== modified file 'aloha/template_files/__init__.py' (properties changed: -x to +x)
=== modified file 'aloha/template_files/aloha_functions_loop.f'
--- aloha/template_files/aloha_functions_loop.f	2013-11-29 07:28:53 +0000
+++ aloha/template_files/aloha_functions_loop.f	2015-10-20 14:02:27 +0000
@@ -86,6 +86,7 @@
 
          pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
 
+
          if ( pp.eq.rZero ) then
 
             sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
@@ -477,13 +478,13 @@
 
             sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
             sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
-            ip = -((1+nh)/2)
-            im =  (1-nh)/2
 
-            fo(5) = im     * sqm(im)
-            fo(6) = ip*nsf * sqm(im)
-            fo(7) = im*nsf * sqm(-ip)
-            fo(8) = ip     * sqm(-ip)
+            im = nhel * (1+nh)/2
+            ip = nhel * -1 * ((1-nh)/2)
+            fo(5) = im     * sqm(abs(ip))
+            fo(6) = ip*nsf * sqm(abs(ip))
+            fo(7) = im*nsf * sqm(abs(im))
+            fo(8) = ip     * sqm(abs(im))
 
          else
 

=== modified file 'aloha/template_files/wavefunctions.py' (properties changed: -x to +x)
--- aloha/template_files/wavefunctions.py	2013-07-23 16:48:52 +0000
+++ aloha/template_files/wavefunctions.py	2015-10-20 14:02:27 +0000
@@ -316,7 +316,7 @@
     elif nhel == 1:
         for j in range(4):
             for i in range(4): 
-                 ft[(i,j)] = sqh*( ep[i]*e0[j] + e0[i]*ep[j] )
+                ft[(i,j)] = sqh*( ep[i]*e0[j] + e0[i]*ep[j] )
     elif nhel == 0:
         for j in range(4):
             for i in range(4):       
@@ -324,10 +324,10 @@
     elif nhel == -1:
         for j in range(4):
             for i in range(4): 
-            	ft[(i,j)] = sqh*( em[i]*e0[j] + e0[i]*em[j] )
+                ft[(i,j)] = sqh*( em[i]*e0[j] + e0[i]*em[j] )
 
     else:
-    	raise Exception, 'invalid helicity TXXXXXX' 
+        raise Exception, 'invalid helicity TXXXXXX' 
 
 
 

=== modified file 'bin/mg5'
--- bin/mg5	2014-07-17 13:34:15 +0000
+++ bin/mg5	2015-10-20 14:02:27 +0000
@@ -28,6 +28,6 @@
 
 sys.argv.pop(0)
 if __debug__:
-    os.system('%s %s %s' %(sys.executable, str(exe_path) , ' '.join(sys.argv) ))
+    os.system('%s -tt %s %s' %(sys.executable, str(exe_path) , ' '.join(sys.argv) ))
 else:
     os.system('%s -O %s %s' %(sys.executable, str(exe_path) , ' '.join(sys.argv) ))

=== modified file 'madgraph/__init__.py' (properties changed: -x to +x)
=== modified file 'madgraph/core/__init__.py' (properties changed: -x to +x)
=== modified file 'madgraph/core/base_objects.py' (properties changed: -x to +x)
--- madgraph/core/base_objects.py	2015-08-16 22:19:30 +0000
+++ madgraph/core/base_objects.py	2015-10-20 14:02:27 +0000
@@ -202,8 +202,8 @@
 
     sorted_keys = ['name', 'antiname', 'spin', 'color',
                    'charge', 'mass', 'width', 'pdg_code',
-                   'texname', 'antitexname', 'line', 'propagating', 'propagator',
-                   'is_part', 'self_antipart', 'ghost', 'counterterm']
+                  'line', 'propagator',
+                   'is_part', 'self_antipart', 'type', 'counterterm']
 
     def default_setup(self):
         """Default values for all properties"""
@@ -216,19 +216,55 @@
         self['mass'] = 'ZERO'
         self['width'] = 'ZERO'
         self['pdg_code'] = 0
-        self['texname'] = 'none'
-        self['antitexname'] = 'none'
+        #self['texname'] = 'none'
+        #self['antitexname'] = 'none'
         self['line'] = 'dashed'
-        self['propagating'] = True
+        #self['propagating'] = True -> removed in favor or 'line' = None
         self['propagator'] = ''
         self['is_part'] = True
         self['self_antipart'] = False
         # True if ghost, False otherwise
-        self['ghost'] = False
+        #self['ghost'] = False
+        self['type'] = '' # empty means normal can also be ghost or goldstone 
         # Counterterm defined as a dictionary with format:
         # ('ORDER_OF_COUNTERTERM',((Particle_list_PDG))):{laurent_order:CTCouplingName}
         self['counterterm'] = {}
 
+    def get(self, name):
+        
+        if name == 'ghost':
+            return self['type'] == 'ghost'
+        elif name == 'goldstone':
+            return self['type'] == 'goldstone'
+        elif name == 'propagating':
+            return self['line'] is not None
+        else:
+            return super(Particle, self).get(name)
+
+    def set(self, name, value, force=False):
+        
+        if name in ['texname', 'antitexname']:
+            return True
+        elif name == 'propagating':
+            if not value:
+                return self.set('line', None, force=force)
+            elif not self.get('line'):
+                return self.set('line', 'dashed',force=force)
+            return True
+        elif name  in ['ghost', 'goldstone']:
+            if self.get('type') == name:
+                if value:
+                    return True
+                else:
+                    return self.set('type', '', force=force)
+            else:
+                if value:
+                    return self.set('type', name, force=force)
+                else:
+                    return True
+        return super(Particle, self).set(name, value,force=force)
+        
+
     def filter(self, name, value):
         """Filter for valid particle property values."""
 
@@ -1213,7 +1249,7 @@
             if isinstance(id, int):
                 try:
                     return self.get("particle_dict")[id]
-                except Exception:
+                except Exception,error:
                     return None
             else:
                 if not hasattr(self, 'name2part'):
@@ -1547,9 +1583,9 @@
     def write_param_card(self):
         """Write out the param_card, and return as string."""
         
-        import models.write_param_card as writter
+        import models.write_param_card as writer
         out = StringIO.StringIO() # it's suppose to be written in a file
-        param = writter.ParamCardWriter(self)
+        param = writer.ParamCardWriter(self)
         param.define_output_file(out)
         param.write_card()
         return out.getvalue()
@@ -1577,55 +1613,118 @@
     def change_electroweak_mode(self, mode):
         """Change the electroweak mode. The only valid mode now is external.
         Where in top of the default MW and sw2 are external parameters."""
-        
-        assert mode == "external"
+
+        assert mode in ["external",set(['mz','mw','alpha'])]
         
         try:
             W = self.get('particle_dict')[24]
         except KeyError:
-            raise InvalidCmd('No W particle in the model impossible to change the EW scheme!')
+            raise InvalidCmd('No W particle in the model impossible to '+
+                                                        'change the EW scheme!')
         
-        MW = self.get_parameter(W.get('mass'))
-        if not isinstance(MW, ParamCardVariable):
-            newMW = ParamCardVariable(MW.name, MW.value, 'MASS', [24])
-            if not newMW.value:
-                newMW.value = 80.385
-            #remove the old definition
-            self.get('parameters')[MW.depend].remove(MW)
-            # add the new one
-            self.add_param(newMW, ['external'])
-            
-        # Now check for sw2. if not define bypass this
-        try:
-            sw2 = self.get_parameter('sw2')
-        except KeyError:
+        if mode=='external':
+            MW = self.get_parameter(W.get('mass'))
+            if not isinstance(MW, ParamCardVariable):
+                newMW = ParamCardVariable(MW.name, MW.value, 'MASS', [24])
+                if not newMW.value:
+                    newMW.value = 80.385
+                #remove the old definition
+                self.get('parameters')[MW.depend].remove(MW)
+                # add the new one
+                self.add_param(newMW, ['external'])
+                
+            # Now check for sw2. if not define bypass this
             try:
-                sw2 = self.get_parameter('mdl_sw2')
+                sw2 = self.get_parameter('sw2')
             except KeyError:
-                sw2=None
-        
-        if sw2:
-            newsw2 = ParamCardVariable(sw2.name,sw2.value, 'SMINPUTS', [4])
-            if not newsw2.value:
-                newsw2.value = 0.222246485786
-            #remove the old definition
-            self.get('parameters')[sw2.depend].remove(sw2)
-            # add the new one
-            self.add_param(newsw2, ['external'])            
-
-    def change_mass_to_complex_scheme(self):
+                try:
+                    sw2 = self.get_parameter('mdl_sw2')
+                except KeyError:
+                    sw2=None
+    
+            if sw2:
+                newsw2 = ParamCardVariable(sw2.name,sw2.value, 'SMINPUTS', [4])
+                if not newsw2.value:
+                    newsw2.value = 0.222246485786
+                #remove the old definition
+                self.get('parameters')[sw2.depend].remove(sw2)
+                # add the new one
+                self.add_param(newsw2, ['external'])
+            # Force a refresh of the parameter dictionary
+            self.parameters_dict = None
+            return true
+
+        elif mode==set(['mz','mw','alpha']):
+            # For now, all we support is to go from mz, Gf, alpha to mz, mw, alpha
+            W = self.get('particle_dict')[24]
+            mass = self.get_parameter(W.get('mass'))
+            mass_expr = 'cmath.sqrt(%(prefix)sMZ__exp__2/2. + cmath.sqrt('+\
+              '%(prefix)sMZ__exp__4/4. - (%(prefix)saEW*cmath.pi*%(prefix)s'+\
+              'MZ__exp__2)/(%(prefix)sGf*%(prefix)ssqrt__2)))'
+            if 'external' in mass.depend:
+                # Nothing to be done
+                return True
+            match = False
+            if mass.expr == mass_expr%{'prefix':''}:
+                prefix = ''
+                match = True
+            elif mass.expr == mass_expr%{'prefix':'mdl_'}:
+                prefix = 'mdl_'
+                match = True
+            if match:
+                MW = ParamCardVariable(mass.name, mass.value, 'MASS', [24])
+                if not MW.value:
+                    MW.value = 80.385                
+                self.get('parameters')[('external',)].append(MW)
+                self.get('parameters')[mass.depend].remove(mass)
+                # Make Gf an internal parameter
+                new_param = ModelVariable('Gf',
+                '-%(prefix)saEW*%(prefix)sMZ**2*cmath.pi/(cmath.sqrt(2)*%(MW)s**2*(%(MW)s**2 - %(prefix)sMZ**2))' %\
+                {'MW': mass.name,'prefix':prefix}, 'complex', mass.depend)
+                Gf = self.get_parameter('%sGf'%prefix)
+                self.get('parameters')[('external',)].remove(Gf)
+                self.add_param(new_param, ['%saEW'%prefix])
+                # Force a refresh of the parameter dictionary
+                self.parameters_dict = None
+                return True
+            else:
+                return False
+
+    def change_mass_to_complex_scheme(self, toCMS=True):
         """modify the expression changing the mass to complex mass scheme"""
         
-        # 1) Find All input parameter mass and width associated
+        # 1) Change the 'CMSParam' of loop_qcd_qed model to 1.0 so as to remove
+        #    the 'real' prefix fromall UVCT wf renormalization expressions.
+        #    If toCMS is False, then it makes sure CMSParam is 0.0 and returns
+        #    immediatly.
+        # 2) Find All input parameter mass and width associated
         #   Add a internal parameter and replace mass with that param
-        # 2) Find All mass fixed by the model and width associated
+        # 3) Find All mass fixed by the model and width associated
         #   -> Both need to be fixed with a real() /Imag()
-        # 3) Find All width fixed by the model
-        #   -> Need to be fixed with a real()
-        # 4) Fix the Yukawa mass to the value of the complex mass/ real mass
-        # 5) Loop through all expression and modify those accordingly
+        # 4) Find All width set by the model
+        #   -> Need to be set with a real()
+        # 5) Fix the Yukawa mass to the value of the complex mass/ real mass
+        # 6) Loop through all expression and modify those accordingly
         #    Including all parameter expression as complex
-
+        
+        try:
+            CMSParam = self.get_parameter('CMSParam')
+        except KeyError:
+            try:
+                CMSParam = self.get_parameter('mdl_CMSParam')
+            except KeyError:
+                CMSParam = None
+                
+        # Handle the case where we want to make sure the CMS is turned off
+        if not toCMS:
+            if CMSParam:
+                CMSParam.expr = '0.0'
+            return            
+        
+        # Now handle the case where we want to turn to CMS.    
+        if CMSParam:
+            CMSParam.expr = '1.0'
+        
         to_change = {}
         mass_widths = [] # parameter which should stay real
         for particle in self.get('particles'):
@@ -1634,52 +1733,29 @@
                 continue
             mass_widths.append(particle.get('width'))
             mass_widths.append(particle.get('mass'))
-            if particle.get('width') == 'ZERO':
+            width = self.get_parameter(particle.get('width'))
+            if (isinstance(width.value, (complex,float)) and abs(width.value)==0.0) or \
+                                                    width.name.lower() =='zero':
                 #everything is fine since the width is zero
                 continue
-            width = self.get_parameter(particle.get('width'))
             if not isinstance(width, ParamCardVariable):
                 width.expr = 're(%s)' % width.expr
-            if particle.get('mass') != 'ZERO':
-                mass = self.get_parameter(particle.get('mass'))
+            mass = self.get_parameter(particle.get('mass'))
+            if (isinstance(width.value, (complex,float)) and abs(width.value)!=0.0) or \
+                                                    mass.name.lower() != 'zero':
                 # special SM treatment to change the gauge scheme automatically.
-                if particle.get('pdg_code') == 24:
-                    if hasattr(mass, 'expr') and mass.expr == 'cmath.sqrt(MZ__exp__2/2. + cmath.sqrt(MZ__exp__4/4. - (aEW*cmath.pi*MZ__exp__2)/(Gf*sqrt__2)))':
-                        # Make MW an external parameter
-                        MW = ParamCardVariable(mass.name, mass.value, 'MASS', [24])
-                        if not MW.value:
-                            MW.value = 80.385
-                        self.get('parameters')[('external',)].append(MW)
-                        self.get('parameters')[mass.depend].remove(mass)
-                        # Make Gf an internal parameter
-                        new_param = ModelVariable('Gf',
-                        '-aEW*MZ**2*cmath.pi/(cmath.sqrt(2)*%(MW)s**2*(%(MW)s**2 - MZ**2))' %\
-                        {'MW': mass.name}, 'complex', mass.depend)
-                        Gf = self.get_parameter('Gf')
-                        self.get('parameters')[('external',)].remove(Gf)
-                        self.add_param(new_param, ['aEW'])
-                        # Use the new mass for the future modification
-                        mass = MW
-                    #option with prefixing
-                    elif hasattr(mass, 'expr') and mass.expr == 'cmath.sqrt(mdl_MZ__exp__2/2. + cmath.sqrt(mdl_MZ__exp__4/4. - (mdl_aEW*cmath.pi*mdl_MZ__exp__2)/(mdl_Gf*mdl_sqrt__2)))':
-                        # Make MW an external parameter
-                        MW = ParamCardVariable(mass.name, mass.value, 'MASS', [24])
-                        if not MW.value:
-                            MW.value = 80.385
-                        self.get('parameters')[('external',)].append(MW)
-                        self.get('parameters')[mass.depend].remove(mass)
-                        # Make Gf an internal parameter
-                        new_param = ModelVariable('mdl_Gf',
-                        '-mdl_aEW*mdl_MZ**2*cmath.pi/(cmath.sqrt(2)*%(MW)s**2*(%(MW)s**2 - mdl_MZ**2))' %\
-                        {'MW': mass.name}, 'complex', mass.depend)
-                        Gf = self.get_parameter('mdl_Gf')
-                        self.get('parameters')[('external',)].remove(Gf)
-                        self.add_param(new_param, ['mdl_aEW'])
-                        # Use the new mass for the future modification
-                        mass = MW
-                    elif isinstance(mass, ModelVariable):
-                        logger.warning('W mass is not an external parameter. This is not adviced for the complex mass scheme.')
-                
+                if particle.get('pdg_code') == 24 and isinstance(mass, 
+                                                                 ModelVariable):
+                    status = self.change_electroweak_mode(
+                                                   set(['mz','mw','alpha']))
+                    # Use the newly defined parameter for the W mass
+                    mass = self.get_parameter(particle.get('mass'))
+                    if not status:
+                        logger.warning('The W mass is not an external '+
+                        'parameter in this model and the automatic change of'+
+                        ' electroweak scheme changed. This is not advised for '+
+                                            'applying the complex mass scheme.')
+
                 # Add A new parameter CMASS
                 #first compute the dependencies (as,...)
                 depend = list(set(mass.depend + width.depend))
@@ -3638,7 +3714,7 @@
         my_isids = [leg.get('ids') for leg in self.get('legs') \
               if not leg.get('state')]
         my_fsids = [leg.get('ids') for leg in self.get('legs') \
-             if leg.get('state')]
+             if leg.get('state')]            
         for i, is_id in enumerate(initial_state_ids):
             assert is_id in my_isids[i]
         for i, fs_id in enumerate(final_state_ids):

=== modified file 'madgraph/core/color_algebra.py' (properties changed: -x to +x)
=== modified file 'madgraph/core/color_amp.py' (properties changed: -x to +x)
=== modified file 'madgraph/core/diagram_generation.py' (properties changed: -x to +x)
--- madgraph/core/diagram_generation.py	2015-06-28 14:44:22 +0000
+++ madgraph/core/diagram_generation.py	2015-10-20 14:02:27 +0000
@@ -26,7 +26,8 @@
 
 import madgraph.core.base_objects as base_objects
 import madgraph.various.misc as misc
-from madgraph import InvalidCmd
+from madgraph import InvalidCmd, MadGraph5Error
+
 logger = logging.getLogger('madgraph.diagram_generation')
 
 
@@ -1552,7 +1553,7 @@
                         "%s is not a valid ProcessDefinitionList object" % str(value)
 
         if name == 'amplitudes':
-            if not isinstance(value, diagram_generation.AmplitudeList):
+            if not isinstance(value, AmplitudeList):
                 raise self.PhysicsObjectError, \
                         "%s is not a valid AmplitudeList object" % str(value)
 

=== modified file 'madgraph/core/drawing.py' (properties changed: -x to +x)
--- madgraph/core/drawing.py	2014-10-07 05:52:01 +0000
+++ madgraph/core/drawing.py	2015-10-20 14:02:27 +0000
@@ -1131,8 +1131,8 @@
             return []
 
         assert self.min_level <= level <= self.max_level , \
-        	'Incorrect value of min/max level: %s <= %s <= %s' % \
-                                        (self.min_level, level, self.max_level)  
+                    'Incorrect value of min/max level: %s <= %s <= %s' % \
+                    (self.min_level, level, self.max_level)  
                 
 
         # At final level we should authorize min=0 and max=1 position    
@@ -1188,7 +1188,7 @@
         # Assign position to each vertex
         for i, vertex in enumerate(vertex_at_level):
             vertex.def_position((level - self.min_level) / self.nb_level, 
-            							     min + dist * (begin_gap + i))
+                                min + dist * (begin_gap + i))
 
         return vertex_at_level
 

=== modified file 'madgraph/core/helas_objects.py' (properties changed: -x to +x)
--- madgraph/core/helas_objects.py	2015-07-02 02:29:22 +0000
+++ madgraph/core/helas_objects.py	2015-10-20 14:02:27 +0000
@@ -3164,7 +3164,7 @@
         """Calculate the actual coupling orders of this diagram"""
 
         wavefunctions = HelasWavefunctionList.extract_wavefunctions(\
-            self.get('amplitudes')[0].get('mothers'))
+                                       self.get('amplitudes')[0].get('mothers'))
 
         coupling_orders = {}
         for wf in wavefunctions + [self.get('amplitudes')[0]]:

=== modified file 'madgraph/fks/__init__.py' (properties changed: -x to +x)
=== modified file 'madgraph/fks/fks_base.py' (properties changed: -x to +x)
--- madgraph/fks/fks_base.py	2015-05-25 20:43:18 +0000
+++ madgraph/fks/fks_base.py	2015-10-20 14:02:27 +0000
@@ -132,6 +132,9 @@
         perturbation = []
         for procdef in self['process_definitions']:
             soft_particles = []
+            # do not warn for decay processes
+            if [ i['state'] for i in procdef['legs']].count(False) == 1:
+                continue
             for pert in procdef['perturbation_couplings']:
                 if pert not in perturbation:
                     perturbation.append(pert)
@@ -550,10 +553,19 @@
         """finds the FKS real configurations for a given process"""
         if range(len(self.leglist)) != [l['number']-1 for l in self.leglist]:
             raise fks_common.FKSProcessError('Disordered numbers of leglist')
+
+        if [ i['state'] for i in self.leglist].count(False) == 1:
+            decay_process=True
+        else:
+            decay_process=False
+
         for i in self.leglist:
             i_i = i['number'] - 1
             self.reals.append([])
-            self.splittings[i_i] = fks_common.find_splittings(i, self.born_proc['model'], {}, pert_order)
+            if decay_process and not i['state']:
+                self.splittings[i_i]=[]
+            else:
+                self.splittings[i_i] = fks_common.find_splittings(i, self.born_proc['model'], {}, pert_order)
             for split in self.splittings[i_i]:
                 self.reals[i_i].append(
                             fks_common.insert_legs(self.leglist, i, split,pert=pert_order))

=== modified file 'madgraph/fks/fks_common.py' (properties changed: -x to +x)
=== modified file 'madgraph/fks/fks_helas_objects.py' (properties changed: -x to +x)
=== modified file 'madgraph/interface/__init__.py' (properties changed: -x to +x)
=== modified file 'madgraph/interface/amcatnlo_interface.py' (properties changed: -x to +x)
--- madgraph/interface/amcatnlo_interface.py	2015-08-12 20:50:38 +0000
+++ madgraph/interface/amcatnlo_interface.py	2015-10-20 14:02:27 +0000
@@ -30,6 +30,7 @@
 import madgraph.interface.extended_cmd as cmd
 import madgraph.interface.madgraph_interface as mg_interface
 import madgraph.interface.madevent_interface as me_interface
+import madgraph.interface.extended_cmd as extended_cmd
 import madgraph.interface.amcatnlo_run_interface as run_interface
 import madgraph.interface.launch_ext_program as launch_ext
 import madgraph.interface.loop_interface as Loop_interface
@@ -542,8 +543,13 @@
                         me.get('processes')[0].set('uid', uid)
                         for fksreal in me.real_processes:
                         # Pick out all initial state particles for the two beams
-                            initial_states.append(sorted(list(set((p.get_initial_pdg(1),p.get_initial_pdg(2)) for \
-                                                             p in fksreal.matrix_element.get('processes')))))
+                            try:
+                                initial_states.append(sorted(list(set((p.get_initial_pdg(1),p.get_initial_pdg(2)) for \
+                                                             p in fksreal.matrix_element.get('processes')))))
+                            except IndexError:
+                                initial_states.append(sorted(list(set((p.get_initial_pdg(1)) for \
+                                                             p in fksreal.matrix_element.get('processes')))))
+                                
                         
                     # remove doubles from the list
                     checked = []
@@ -569,10 +575,10 @@
 
             #_curr_matrix_element is a FKSHelasMultiProcess Object 
             self._fks_directories = []
-            proc_charac = banner_mod.ProcCharacteristic()
+            proc_charac = self._curr_exporter.proc_characteristic 
             for charac in ['has_isr', 'has_fsr', 'has_loops']:
                 proc_charac[charac] = self._curr_matrix_elements[charac]
-            proc_charac.write(pjoin(path, 'proc_characteristics'))
+
 
             for ime, me in \
                 enumerate(self._curr_matrix_elements.get('matrix_elements')):
@@ -624,7 +630,7 @@
 #            self.options['automatic_html_opening'] = False
 
         if options['interactive']:
-            if hasattr(self, 'do_shell'):
+            if isinstance(self, extended_cmd.CmdShell):
                 ME = run_interface.aMCatNLOCmdShell(me_dir=argss[0], options=self.options)
             else:
                 ME = run_interface.aMCatNLOCmd(me_dir=argss[0],options=self.options)
@@ -636,7 +642,9 @@
             stop = self.define_child_cmd_interface(ME)                
             return stop
 
-        ext_program = launch_ext.aMCatNLOLauncher(argss[0], self, run_mode=argss[1], **options)
+        ext_program = launch_ext.aMCatNLOLauncher(argss[0], self, run_mode=argss[1],
+                                                  shell = isinstance(self, extended_cmd.CmdShell),
+                                                  **options)
         ext_program.run()
         
                     

=== modified file 'madgraph/interface/amcatnlo_run_interface.py' (properties changed: -x to +x)
--- madgraph/interface/amcatnlo_run_interface.py	2015-09-05 17:49:40 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2015-10-20 14:02:27 +0000
@@ -119,6 +119,11 @@
             #this can be improved/better written to handle the output
             misc.call(['./%s' % (test)], cwd=this_dir, 
                     stdin = open(input), stdout=open(pjoin(this_dir, '%s.log' % test), 'w'))
+            if test == 'check_poles' and os.path.exists(pjoin(this_dir,'MadLoop5_resources')) :
+                tf=tarfile.open(pjoin(this_dir,'MadLoop5_resources.tar.gz'),'w:gz',
+                                                                 dereference=True)
+                tf.add(pjoin(this_dir,'MadLoop5_resources'),arcname='MadLoop5_resources')
+                tf.close()
             
         if not options['reweightonly']:
             misc.compile(['gensym'], cwd=this_dir, job_specs = False)
@@ -1155,7 +1160,7 @@
         functions, such as generate_events or calculate_xsect
         mode gives the list of switch needed for the computation (usefull for banner_run)
         """
-        
+
         if not argss and not options:
             self.start_time = time.time()
             argss = self.split_arg(line)
@@ -1164,6 +1169,7 @@
             options = options.__dict__
             self.check_launch(argss, options)
 
+        
         if 'run_name' in options.keys() and options['run_name']:
             self.run_name = options['run_name']
             # if a dir with the given run_name already exists
@@ -1183,6 +1189,7 @@
         
         if not switch:
             mode = argss[0]
+
             if mode in ['LO', 'NLO']:
                 options['parton'] = True
             mode = self.ask_run_configuration(mode, options)
@@ -1202,7 +1209,7 @@
         self.compile(mode, options) 
         evt_file = self.run(mode, options)
         
-        if int(self.run_card['nevents']) == 0 and not mode in ['LO', 'NLO']:
+        if self.run_card['nevents'] == 0 and not mode in ['LO', 'NLO']:
             logger.info('No event file generated: grids have been set-up with a '\
                             'relative precision of %s' % self.run_card['req_acc'])
             return
@@ -1222,14 +1229,40 @@
 
 
         self.update_status('', level='all', update_results=True)
-        if int(self.run_card['ickkw']) == 3 and mode in ['noshower', 'aMC@NLO']:
+        if self.run_card['ickkw'] == 3 and mode in ['noshower', 'aMC@NLO']:
             logger.warning("""You are running with FxFx merging enabled.
 To be able to merge samples of various multiplicities without double counting,
 you have to remove some events after showering 'by hand'.
 Please read http://amcatnlo.cern.ch/FxFx_merging.htm for more details.""")
 
 
-
+        #check if the param_card defines a scan.
+        if self.param_card_iterator:
+            param_card_iterator = self.param_card_iterator
+            self.param_card_iterator = [] #avoid to next generate go trough here
+            param_card_iterator.store_entry(self.run_name, self.results.current['cross'])
+            orig_name = self.run_name
+            #go trough the scal
+            for i,card in enumerate(param_card_iterator):
+                card.write(pjoin(self.me_dir,'Cards','param_card.dat'))
+                if not options['force']:
+                    options['force'] = True
+                if options['run_name']:
+                    options['run_name'] = '%s_%s' % (orig_name, i+1)
+                if not argss:
+                    argss = [mode, "-f"]
+                elif argss[0] == "auto":
+                    argss[0] = mode
+                self.do_launch("", options=options, argss=argss, switch=switch)
+                #self.exec_cmd("launch -f ",precmd=True, postcmd=True,errorhandling=False)
+                param_card_iterator.store_entry(self.run_name, self.results.current['cross'])
+            #restore original param_card
+            param_card_iterator.write(pjoin(self.me_dir,'Cards','param_card.dat'))
+            name = misc.get_scan_name(orig_name, self.run_name)
+            path = pjoin(self.me_dir, 'Events','scan_%s.txt' % name)
+            logger.info("write all cross-section results in %s" % path, '$MG:color:BLACK')
+            param_card_iterator.write_summary(path)
+            
     ############################################################################      
     def do_compile(self, line):
         """Advanced commands: just compile the executables """
@@ -1248,59 +1281,11 @@
 
         self.update_status('', level='all', update_results=True)
 
-    def print_results_in_shell(self, data):
-        """Have a nice results prints in the shell,
-        data should be of type: gen_crossxhtml.OneTagResults"""
-        if not data:
-            return
-        logger.info("  === Results Summary for run: %s tag: %s ===\n" % (data['run_name'],data['tag']))
-        if self.ninitial == 1:
-            logger.info("     Width :   %.4g +- %.4g GeV" % (data['cross'], data['error']))
-        else:
-            logger.info("     Cross-section :   %.4g +- %.4g pb" % (data['cross'], data['error']))
-        logger.info("     Nb of events :  %s" % data['nb_event'] )
-        #if data['cross_pythia'] and data['nb_event_pythia']:
-        #    if self.ninitial == 1:
-        #        logger.info("     Matched Width :   %.4g +- %.4g GeV" % (data['cross_pythia'], data['error_pythia']))
-        #    else:
-        #        logger.info("     Matched Cross-section :   %.4g +- %.4g pb" % (data['cross_pythia'], data['error_pythia']))            
-        #    logger.info("     Nb of events after Matching :  %s" % data['nb_event_pythia'])
-        #    if self.run_card['use_syst'] in self.true:
-        #        logger.info("     Be carefull that matched information are here NOT for the central value. Refer to SysCalc output for it")    
-        logger.info(" " )
-
-    def print_results_in_file(self, data, path, mode='w'):
-        """Have a nice results prints in the shell,
-        data should be of type: gen_crossxhtml.OneTagResults"""
-        if not data:
-            return
-        
-        fsock = open(path, mode)
-        
-        fsock.write("  === Results Summary for run: %s tag: %s  process: %s ===\n" % \
-                    (data['run_name'],data['tag'], os.path.basename(self.me_dir)))
-        
-        if self.ninitial == 1:
-            fsock.write("     Width :   %.4g +- %.4g GeV\n" % (data['cross'], data['error']))
-        else:
-            fsock.write("     Cross-section :   %.4g +- %.4g pb\n" % (data['cross'], data['error']))
-        fsock.write("     Nb of events :  %s\n" % data['nb_event'] )
-        #if data['cross_pythia'] and data['nb_event_pythia']:
-        #    if self.ninitial == 1:
-        #        fsock.write("     Matched Width :   %.4g +- %.4g GeV\n" % (data['cross_pythia'], data['error_pythia']))
-        #    else:
-        #        fsock.write("     Matched Cross-section :   %.4g +- %.4g pb\n" % (data['cross_pythia'], data['error_pythia']))            
-        #    fsock.write("     Nb of events after Matching :  %s\n" % data['nb_event_pythia'])
-        fsock.write(" \n" )
-
-
-
-
 
     def update_random_seed(self):
         """Update random number seed with the value from the run_card. 
         If this is 0, update the number according to a fresh one"""
-        iseed = int(self.run_card['iseed'])
+        iseed = self.run_card['iseed']
         if iseed == 0:
             randinit = open(pjoin(self.me_dir, 'SubProcesses', 'randinit'))
             iseed = int(randinit.read()[2:]) + 1
@@ -1317,206 +1302,95 @@
         if not 'only_generation' in options.keys():
             options['only_generation'] = False
 
+        # for second step in applgrid mode, do only the event generation step
         if mode in ['LO', 'NLO'] and self.run_card['iappl'] == 2 and not options['only_generation']:
             options['only_generation'] = True
         self.get_characteristics(pjoin(self.me_dir, 'SubProcesses', 'proc_characteristics'))
-
-        if self.cluster_mode == 1:
-            cluster_name = self.options['cluster_type']
-            self.cluster = cluster.from_name[cluster_name](**self.options)
-        if self.cluster_mode == 2:
-            try:
-                import multiprocessing
-                if not self.nb_core:
-                    try:
-                        self.nb_core = int(self.options['nb_core'])
-                    except TypeError:
-                        self.nb_core = multiprocessing.cpu_count()
-                logger.info('Using %d cores' % self.nb_core)
-            except ImportError:
-                self.nb_core = 1
-                logger.warning('Impossible to detect the number of cores => Using One.\n'+
-                        'Use set nb_core X in order to set this number and be able to'+
-                        'run in multicore.')
-
-            self.cluster = cluster.MultiCore(**self.options)
+        self.setup_cluster_or_multicore()
         self.update_random_seed()
         #find and keep track of all the jobs
         folder_names = {'LO': ['born_G*'], 'NLO': ['all_G*'],
                     'aMC@LO': ['GB*'], 'aMC@NLO': ['GF*']}
         folder_names['noshower'] = folder_names['aMC@NLO']
         folder_names['noshowerLO'] = folder_names['aMC@LO']
-        job_dict = {}
         p_dirs = [d for d in \
                 open(pjoin(self.me_dir, 'SubProcesses', 'subproc.mg')).read().split('\n') if d]
-        #find jobs and clean previous results
-        if not options['only_generation'] and not options['reweightonly']:
-            self.update_status('Cleaning previous results', level=None)
-        for dir in p_dirs:
-            job_dict[dir] = [file for file in \
-                                 os.listdir(pjoin(self.me_dir, 'SubProcesses', dir)) \
-                                 if file.startswith('ajob')] 
-            #find old folders to be removed
-            for obj in folder_names[mode]:
-                to_rm = [file for file in \
-                             os.listdir(pjoin(self.me_dir, 'SubProcesses', dir)) \
-                             if file.startswith(obj[:-1]) and \
-                            (os.path.isdir(pjoin(self.me_dir, 'SubProcesses', dir, file)) or \
-                             os.path.exists(pjoin(self.me_dir, 'SubProcesses', dir, file)))] 
-                #always clean dirs for the splitted event generation
-                # do not include the born_G/ grid_G which should be kept when
-                # doing a f.o. run keeping old grids
-                to_always_rm = [file for file in \
-                             os.listdir(pjoin(self.me_dir, 'SubProcesses', dir)) \
-                             if file.startswith(obj[:-1]) and
-                             '_' in file and not '_G' in file and \
-                            (os.path.isdir(pjoin(self.me_dir, 'SubProcesses', dir, file)) or \
-                             os.path.exists(pjoin(self.me_dir, 'SubProcesses', dir, file)))]
-
-                if not options['only_generation'] and not options['reweightonly']:
-                    to_always_rm.extend(to_rm)
-                    if os.path.exists(pjoin(self.me_dir, 'SubProcesses', dir,'MadLoop5_resources.tar.gz')):
-                        to_always_rm.append(pjoin(self.me_dir, 'SubProcesses', dir,'MadLoop5_resources.tar.gz'))
-                files.rm([pjoin(self.me_dir, 'SubProcesses', dir, d) for d in to_always_rm])
-
-        mcatnlo_status = ['Setting up grid', 'Computing upper envelope', 'Generating events']
-
-        if self.run_card['iappl'] == 2:
-            self.applgrid_distribute(options,mode,p_dirs)
+        #Clean previous results
+        self.clean_previous_results(options,p_dirs,folder_names[mode])
+
+        mcatnlo_status = ['Setting up grids', 'Computing upper envelope', 'Generating events']
+
 
         if options['reweightonly']:
             event_norm=self.run_card['event_norm']
-            nevents=int(self.run_card['nevents'])
+            nevents=self.run_card['nevents']
             return self.reweight_and_collect_events(options, mode, nevents, event_norm)
 
         devnull = os.open(os.devnull, os.O_RDWR) 
+
         if mode in ['LO', 'NLO']:
             # this is for fixed order runs
             mode_dict = {'NLO': 'all', 'LO': 'born'}
             logger.info('Doing fixed order %s' % mode)
             req_acc = self.run_card['req_acc_FO']
-            if not options['only_generation'] and req_acc != -1:
-                self.write_madin_file(pjoin(self.me_dir, 'SubProcesses'), mode_dict[mode], 0, '-1', '6','0.10') 
-                self.update_status('Setting up grids', level=None)
-                self.run_all(job_dict, [['0', mode_dict[mode], '0']], 'Setting up grids')
-            elif not options['only_generation']:
-                npoints = self.run_card['npoints_FO_grid']
-                niters = self.run_card['niters_FO_grid']
-                self.write_madin_file(pjoin(self.me_dir, 'SubProcesses'), mode_dict[mode], 0, npoints, niters) 
-                self.update_status('Setting up grids', level=None)
-                self.run_all(job_dict, [['0', mode_dict[mode], '0']], 'Setting up grids')
-
-            npoints = self.run_card['npoints_FO']
-            niters = self.run_card['niters_FO']
-            self.write_madin_file(pjoin(self.me_dir, 'SubProcesses'), mode_dict[mode], -1, npoints, niters) 
-            # collect the results and logs
-            self.collect_log_files(folder_names[mode], 0)
-            p = misc.Popen(['./combine_results_FO.sh', str(req_acc), '%s_G*' % mode_dict[mode]], \
-                               stdout=subprocess.PIPE, \
-                               cwd=pjoin(self.me_dir, 'SubProcesses'))
-            output = p.communicate()
-
-            self.cross_sect_dict = self.read_results(output, mode)
-            self.print_summary(options, 0, mode)
-            cross, error = sum_html.make_all_html_results(self, ['%s*' % mode_dict[mode]])
-            self.results.add_detail('cross', cross)
-            self.results.add_detail('error', error) 
-
-            self.update_status('Computing cross-section', level=None)
-            self.run_all(job_dict, [['0', mode_dict[mode], '0', mode_dict[mode]]], 'Computing cross-section')
-
-            # collect the results and logs
-            self.collect_log_files(folder_names[mode], 1)
-            p = misc.Popen(['./combine_results_FO.sh', '-1'] + folder_names[mode], \
-                                stdout=subprocess.PIPE, 
-                                cwd=pjoin(self.me_dir, 'SubProcesses'))
-            output = p.communicate()
-            self.cross_sect_dict = self.read_results(output, mode)
-
-            # collect the scale and PDF uncertainties
-            scale_pdf_info={}
-            if self.run_card['reweight_scale'] or self.run_card['reweight_PDF']:
-                data_files=[]
-                for dir in p_dirs:
-                    for obj in folder_names[mode]:
-                        for file in os.listdir(pjoin(self.me_dir, 'SubProcesses', dir)):
-                            if file.startswith(obj[:-1]) and \
-                                    (os.path.exists(pjoin(self.me_dir, 'SubProcesses', dir, file,'scale_pdf_dependence.dat'))):
-                                data_files.append(pjoin(dir,file,'scale_pdf_dependence.dat'))
-                scale_pdf_info = self.pdf_scale_from_reweighting(data_files)
-            # print the results:
-            self.print_summary(options, 1, mode, scale_pdf_info)
-
-            files.cp(pjoin(self.me_dir, 'SubProcesses', 'res.txt'),
-                     pjoin(self.me_dir, 'Events', self.run_name))
-            
-            if self.analyse_card['fo_analysis_format'].lower() == 'topdrawer':
-                misc.call(['./combine_plots_FO.sh'] + folder_names[mode], \
-                                stdout=devnull, 
-                                cwd=pjoin(self.me_dir, 'SubProcesses'))
-                files.cp(pjoin(self.me_dir, 'SubProcesses', 'MADatNLO.top'),
-                                pjoin(self.me_dir, 'Events', self.run_name))
-                logger.info('The results of this run and the TopDrawer file with the plots' + \
-                        ' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
-            elif self.analyse_card['fo_analysis_format'].lower() == 'hwu':
-                self.combine_plots_HwU(folder_names[mode])
-                files.cp(pjoin(self.me_dir, 'SubProcesses', 'MADatNLO.HwU'),
-                                pjoin(self.me_dir, 'Events', self.run_name))
-                files.cp(pjoin(self.me_dir, 'SubProcesses', 'MADatNLO.gnuplot'),
-                                pjoin(self.me_dir, 'Events', self.run_name))
-                try:
-                    misc.call(['gnuplot','MADatNLO.gnuplot'],\
-                              stdout=os.open(os.devnull, os.O_RDWR),\
-                              stderr=os.open(os.devnull, os.O_RDWR),\
-                              cwd=pjoin(self.me_dir, 'Events', self.run_name))
-                except Exception:
-                    pass
-
-
-                logger.info('The results of this run and the HwU and GnuPlot files with the plots' + \
-                        ' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
-            elif self.analyse_card['fo_analysis_format'].lower() == 'root':
-                misc.call(['./combine_root.sh'] + folder_names[mode], \
-                                stdout=devnull, 
-                                cwd=pjoin(self.me_dir, 'SubProcesses'))
-                files.cp(pjoin(self.me_dir, 'SubProcesses', 'MADatNLO.root'),
-                                pjoin(self.me_dir, 'Events', self.run_name))
-                logger.info('The results of this run and the ROOT file with the plots' + \
-                        ' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
-            else:
-                logger.info('The results of this run' + \
-                            ' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
-                
-            cross, error = sum_html.make_all_html_results(self, folder_names[mode])
-            self.results.add_detail('cross', cross)
-            self.results.add_detail('error', error)
-            if self.run_card['iappl'] != 0:
-                self.applgrid_combine(cross,error)
+
+            # Re-distribute the grids for the 2nd step of the applgrid
+            # running
+            if self.run_card['iappl'] == 2:
+                self.applgrid_distribute(options,mode_dict[mode],p_dirs)
+
+            # create a list of dictionaries "jobs_to_run" with all the
+            # jobs that need to be run
+            integration_step=-1
+            jobs_to_run,jobs_to_collect,integration_step = self.create_jobs_to_run(options,p_dirs, \
+                                req_acc,mode_dict[mode],integration_step,mode,fixed_order=True)
+            self.prepare_directories(jobs_to_run,mode)
+
+            # loop over the integration steps. After every step, check
+            # if we have the required accuracy. If this is the case,
+            # stop running, else do another step.
+            while True:
+                integration_step=integration_step+1
+                self.run_all_jobs(jobs_to_run,integration_step)
+                self.collect_log_files(jobs_to_run,integration_step)
+                jobs_to_run,jobs_to_collect=self.collect_the_results(options,req_acc,jobs_to_run, \
+                                  jobs_to_collect,integration_step,mode,mode_dict[mode])
+                if not jobs_to_run:
+                    # there are no more jobs to run (jobs_to_run is empty)
+                    break
+            # We are done.
+            self.finalise_run_FO(folder_names[mode],jobs_to_collect)
             self.update_status('Run complete', level='parton', update_results=True)
-
             return
 
         elif mode in ['aMC@NLO','aMC@LO','noshower','noshowerLO']:
+            if self.ninitial == 1:
+                raise aMCatNLOError('Decay processes can only be run at fixed order.')
+            mode_dict = {'aMC@NLO': 'all', 'aMC@LO': 'born',\
+                         'noshower': 'all', 'noshowerLO': 'born'}
             shower = self.run_card['parton_shower'].upper()
-            nevents = int(self.run_card['nevents'])
+            nevents = self.run_card['nevents']
             req_acc = self.run_card['req_acc']
-            if nevents == 0 and float(req_acc) < 0 :
+            if nevents == 0 and req_acc < 0 :
                 raise aMCatNLOError('Cannot determine the required accuracy from the number '\
                                         'of events, because 0 events requested. Please set '\
-                                        'the "req_acc" parameter in the run_card to a value between 0 and 1')
-            elif float(req_acc) >1 or float(req_acc) == 0 :
+                                        'the "req_acc" parameter in the run_card to a value '\
+                                        'between 0 and 1')
+            elif req_acc >1 or req_acc == 0 :
                 raise aMCatNLOError('Required accuracy ("req_acc" in the run_card) should '\
                                         'be between larger than 0 and smaller than 1, '\
-                                        'or set to -1 for automatic determination. Current value is %s' % req_acc)
+                                        'or set to -1 for automatic determination. Current '\
+                                        'value is %f' % req_acc)
 # For more than 1M events, set req_acc to 0.001 (except when it was explicitly set in the run_card)
-            elif float(req_acc) < 0 and nevents > 1000000 :
-                req_acc='0.001'
+            elif req_acc < 0 and nevents > 1000000 :
+                req_acc=0.001
 
             shower_list = ['HERWIG6', 'HERWIGPP', 'PYTHIA6Q', 'PYTHIA6PT', 'PYTHIA8']
 
             if not shower in shower_list:
-                raise aMCatNLOError('%s is not a valid parton shower. Please use one of the following: %s' \
-                    % (shower, ', '.join(shower_list)))
+                raise aMCatNLOError('%s is not a valid parton shower. '\
+                                    'Please use one of the following: %s' \
+                                    % (shower, ', '.join(shower_list)))
 
 # check that PYTHIA6PT is not used for processes with FSR
             if shower == 'PYTHIA6PT' and self.proc_characteristics['has_fsr']:
@@ -1529,98 +1403,556 @@
             elif options['only_generation']:
                 logger.info('Generating events starting from existing results')
             
-
-            for i, status in enumerate(mcatnlo_status):
-                #check if need to split jobs 
-                # at least one channel must have enough events
-                try:
-                    nevents_unweighted = open(pjoin(self.me_dir, 
-                                                'SubProcesses', 
-                                                'nevents_unweighted')).read().split('\n')
-                except IOError:
-                    nevents_unweighted = []
-
-                split = i == 2 and \
-                        int(self.run_card['nevt_job']) > 0 
-
-                if i == 2 or not options['only_generation']:
-                    # if the number of events requested is zero,
-                    # skip mint step 2
-                    if i==2 and nevents==0:
-                        self.print_summary(options, 2,mode)
-                        return
-
-                    if split:
-                        # split the event generation
-                        misc.call([pjoin(self.me_dir, 'bin', 'internal', 'split_jobs.py')] + \
-                                   [str(self.run_card['nevt_job'])],
-                                   stdout = devnull,
-                                   cwd = pjoin(self.me_dir, 'SubProcesses'))
-                        assert os.path.exists(pjoin(self.me_dir, 'SubProcesses', 
-                            'nevents_unweighted_splitted'))
-
-                    self.update_status(status, level='parton')
-                    if mode in ['aMC@NLO', 'noshower']:
-                        self.write_madinMMC_file(pjoin(self.me_dir, 'SubProcesses'), 'all', i) 
-                        self.run_all(job_dict, [['2', 'F', '%d' % i]], status, split_jobs = split)
-                        
-                    elif mode in ['aMC@LO', 'noshowerLO']:
-                        self.write_madinMMC_file(
-                            pjoin(self.me_dir, 'SubProcesses'), 'born', i) 
-                        self.run_all(job_dict, 
-                                     [['2', 'B', '%d' % i]], 
-                                     '%s at LO' % status, split_jobs = split)
-
-                if (i < 2 and not options['only_generation']) or i == 1 :
-                    # collect the results and logs
-                    self.collect_log_files(folder_names[mode], i)
-                    p = misc.Popen(['./combine_results.sh'] + \
-                                   ['%d' % i,'%d' % nevents, '%s' % req_acc ] + \
-                                   folder_names[mode],
-                                   stdout=subprocess.PIPE, 
-                                   cwd = pjoin(self.me_dir, 'SubProcesses'))
-                    output = p.communicate()
-                    files.cp(pjoin(self.me_dir, 'SubProcesses', 'res_%d.txt' % i), \
-                             pjoin(self.me_dir, 'Events', self.run_name))
-
-                    self.cross_sect_dict = self.read_results(output, mode)
-                    self.print_summary(options, i, mode)
-
-                    cross, error = sum_html.make_all_html_results(self, folder_names[mode])
-                    self.results.add_detail('cross', cross)
-                    self.results.add_detail('error', error) 
-
-                #check that split jobs are all correctly terminated
-                if split:
-                    self.check_event_files()
-
-        if self.cluster_mode == 1:
-            #if cluster run, wait 15 sec so that event files are transferred back
-            self.update_status(
+            jobs_to_run,jobs_to_collect,integration_step = self.create_jobs_to_run(options,p_dirs, \
+                                            req_acc,mode_dict[mode],1,mode,fixed_order=False)
+
+            # Make sure to update all the jobs to be ready for the event generation step
+            if options['only_generation']:
+                jobs_to_run,jobs_to_collect=self.collect_the_results(options,req_acc,jobs_to_run, \
+                                jobs_to_collect,1,mode,mode_dict[mode],fixed_order=False)
+            else:
+                self.prepare_directories(jobs_to_run,mode,fixed_order=False)
+
+
+            # Main loop over the three MINT generation steps:
+            for mint_step, status in enumerate(mcatnlo_status):
+                if options['only_generation'] and mint_step < 2:
+                    continue
+                self.update_status(status, level='parton')
+                self.run_all_jobs(jobs_to_run,mint_step,fixed_order=False)
+                self.collect_log_files(jobs_to_run,mint_step)
+                if mint_step+1==2 and nevents==0:
+                    self.print_summary(options,2,mode)
+                    return
+                jobs_to_run,jobs_to_collect=self.collect_the_results(options,req_acc,jobs_to_run, \
+                                jobs_to_collect,mint_step,mode,mode_dict[mode],fixed_order=False)
+            # Sanity check on the event files. If error the jobs are resubmitted
+            self.check_event_files(jobs_to_collect)
+
+            if self.cluster_mode == 1:
+            #if cluster run, wait 10 sec so that event files are transferred back
+                self.update_status(
                     'Waiting while files are transferred back from the cluster nodes',
                     level='parton')
-            time.sleep(10)
-        if split:
-            files.cp(pjoin(self.me_dir, 'SubProcesses', 'nevents_unweighted_splitted'), \
-                     pjoin(self.me_dir, 'SubProcesses', 'nevents_unweighted'))
-
-
-        event_norm=self.run_card['event_norm']
-        self.collect_log_files(folder_names[mode], 2)
-        return self.reweight_and_collect_events(options, mode, nevents, event_norm)
-
-    def combine_plots_HwU(self,folder_names):
+                time.sleep(10)
+
+            event_norm=self.run_card['event_norm']
+            return self.reweight_and_collect_events(options, mode, nevents, event_norm)
+
+    def create_jobs_to_run(self,options,p_dirs,req_acc,run_mode,\
+                           integration_step,mode,fixed_order=True):
+        """Creates a list of dictionaries with all the jobs to be run"""
+        jobs_to_run=[]
+        if not options['only_generation']:
+            # Fresh, new run. Check all the P*/channels.txt files
+            # (created by the 'gensym' executable) to set-up all the
+            # jobs using the default inputs.
+            npoints = self.run_card['npoints_FO_grid']
+            niters = self.run_card['niters_FO_grid']
+            for p_dir in p_dirs:
+                try:
+                    with open(pjoin(self.me_dir,'SubProcesses',p_dir,'channels.txt')) as chan_file:
+                        channels=chan_file.readline().split()
+                except IOError:
+                    logger.warning('No integration channels found for contribution %s' % p_dir)
+                    continue
+                for channel in channels:
+                    job={}
+                    job['p_dir']=p_dir
+                    job['channel']=channel
+                    job['split']=0
+                    if fixed_order and req_acc == -1:
+                        job['accuracy']=0
+                        job['niters']=niters
+                        job['npoints']=npoints
+                    elif fixed_order and req_acc > 0:
+                        job['accuracy']=0.10
+                        job['niters']=6
+                        job['npoints']=-1
+                    elif not fixed_order:
+                        job['accuracy']=0.03
+                        job['niters']=12
+                        job['npoints']=-1
+                    else:
+                        raise aMCatNLOError('No consistent "req_acc_FO" set. Use a value '+
+                                            'between 0 and 1 or set it equal to -1.')
+                    job['mint_mode']=0
+                    job['run_mode']=run_mode
+                    job['wgt_frac']=1.0
+                    jobs_to_run.append(job)
+            jobs_to_collect=copy.copy(jobs_to_run) # These are all jobs
+        else:
+            # if options['only_generation'] is true, we need to loop
+            # over all the existing G* directories and create the jobs
+            # from there.
+            name_suffix={'born' :'B', 'all':'F'}
+            for p_dir in p_dirs:
+                for chan_dir in os.listdir(pjoin(self.me_dir,'SubProcesses',p_dir)):
+                    if ((chan_dir.startswith(run_mode+'_G') and fixed_order) or\
+                        (chan_dir.startswith('G'+name_suffix[run_mode]) and (not fixed_order))) and \
+                       (os.path.isdir(pjoin(self.me_dir, 'SubProcesses', p_dir, chan_dir)) or \
+                        os.path.exists(pjoin(self.me_dir, 'SubProcesses', p_dir, chan_dir))):
+                        job={}
+                        job['p_dir']=p_dir
+                        if fixed_order:
+                            channel=chan_dir.split('_')[1]
+                            job['channel']=channel[1:] # remove the 'G'
+                            if len(chan_dir.split('_')) == 3:
+                                split=int(chan_dir.split('_')[2])
+                            else:
+                                split=0
+                        else:
+                            if len(chan_dir.split('_')) == 2:
+                                split=int(chan_dir.split('_')[1])
+                                channel=chan_dir.split('_')[0]
+                                job['channel']=channel[2:] # remove the 'G'
+                            else:
+                                job['channel']=chan_dir[2:] # remove the 'G'
+                                split=0
+                        job['split']=split
+                        job['run_mode']=run_mode
+                        job['dirname']=pjoin(self.me_dir, 'SubProcesses', p_dir, chan_dir)
+                        job['wgt_frac']=1.0
+                        if not fixed_order: job['mint_mode']=1
+                        jobs_to_run.append(job)
+            jobs_to_collect=copy.copy(jobs_to_run) # These are all jobs
+            if fixed_order:
+                jobs_to_run,jobs_to_collect=self.collect_the_results(options,req_acc,jobs_to_run,
+                                                 jobs_to_collect,integration_step,mode,run_mode)
+                # Update the integration_step to make sure that nothing will be overwritten
+                integration_step=1
+                for job in jobs_to_run:
+                    while os.path.exists(pjoin(job['dirname'],'res_%s.dat' % integration_step)):
+                        integration_step=integration_step+1
+                integration_step=integration_step-1
+            else:
+                self.append_the_results(jobs_to_collect,integration_step)
+        return jobs_to_run,jobs_to_collect,integration_step
+
+    def prepare_directories(self,jobs_to_run,mode,fixed_order=True):
+        """Set-up the G* directories for running"""
+        name_suffix={'born' :'B' , 'all':'F'}
+        for job in jobs_to_run:
+            if job['split'] == 0:
+                if fixed_order :
+                    dirname=pjoin(self.me_dir,'SubProcesses',job['p_dir'],
+                                  job['run_mode']+'_G'+job['channel'])
+                else:
+                    dirname=pjoin(self.me_dir,'SubProcesses',job['p_dir'],
+                                  'G'+name_suffix[job['run_mode']]+job['channel'])
+            else:
+                if fixed_order :
+                    dirname=pjoin(self.me_dir,'SubProcesses',job['p_dir'],
+                            job['run_mode']+'_G'+job['channel']+'_'+str(job['split']))
+                else:
+                    dirname=pjoin(self.me_dir,'SubProcesses',job['p_dir'],
+                            'G'+name_suffix[job['run_mode']]+job['channel']+'_'+str(job['split']))
+            job['dirname']=dirname
+            if not os.path.isdir(dirname):
+                os.makedirs(dirname)
+            self.write_input_file(job,fixed_order)
+            if not fixed_order:
+                # copy the grids from the base directory to the split directory:
+                if job['split'] != 0:
+                    for f in ['grid.MC_integer','mint_grids','res_1']:
+                        if not os.path.isfile(pjoin(job['dirname'],f)):
+                            files.ln(pjoin(job['dirname'].rsplit("_",1)[0],f),job['dirname'])
+
+
+    def write_input_file(self,job,fixed_order):
+        """write the input file for the madevent_mint* executable in the appropriate directory"""
+        if fixed_order:
+            content= \
+"""NPOINTS = %(npoints)s
+NITERATIONS = %(niters)s
+ACCURACY = %(accuracy)s
+ADAPT_GRID = 2
+MULTICHANNEL = 1
+SUM_HELICITY = 1
+CHANNEL = %(channel)s
+SPLIT = %(split)s
+RUN_MODE = %(run_mode)s
+RESTART = %(mint_mode)s
+""" \
+              % job
+        else:
+            content = \
+"""-1 12      ! points, iterations
+%(accuracy)s       ! desired fractional accuracy
+1 -0.1     ! alpha, beta for Gsoft
+-1 -0.1    ! alpha, beta for Gazi
+1          ! Suppress amplitude (0 no, 1 yes)?
+1          ! Exact helicity sum (0 yes, n = number/event)?
+%(channel)s          ! Enter Configuration Number:
+%(mint_mode)s          ! MINT imode: 0 to set-up grids, 1 to perform integral, 2 generate events
+1 1 1      ! if imode is 1: Folding parameters for xi_i, phi_i and y_ij
+%(run_mode)s        ! all, born, real, virt
+""" \
+                    % job
+        with open(pjoin(job['dirname'], 'input_app.txt'), 'w') as input_file:
+            input_file.write(content)
+
+
+    def run_all_jobs(self,jobs_to_run,integration_step,fixed_order=True):
+        """Loops over the jobs_to_run and executes them using the function 'run_exe'"""
+        if fixed_order:
+            if integration_step == 0:
+                self.update_status('Setting up grids', level=None)
+            else:
+                self.update_status('Refining results, step %i' % integration_step, level=None)
+        self.ijob = 0
+        name_suffix={'born' :'B', 'all':'F'}
+        if fixed_order:
+            run_type="Fixed order integration step %s" % integration_step
+        else:
+            run_type="MINT step %s" % integration_step
+        for job in jobs_to_run:
+            executable='ajob1'
+            if fixed_order:
+                arguments=[job['channel'],job['run_mode'], \
+                                    str(job['split']),str(integration_step)]
+            else:
+                arguments=[job['channel'],name_suffix[job['run_mode']], \
+                                    str(job['split']),str(integration_step)]
+            self.run_exe(executable,arguments,run_type,
+                         cwd=pjoin(self.me_dir,'SubProcesses',job['p_dir']))
+
+        if self.cluster_mode == 2:
+            time.sleep(1) # security to allow all jobs to be launched
+        self.njobs=len(jobs_to_run)
+        self.wait_for_complete(run_type)
+
+
+    def collect_the_results(self,options,req_acc,jobs_to_run,jobs_to_collect,\
+                            integration_step,mode,run_mode,fixed_order=True):
+        """Collect the results, make HTML pages, print the summary and
+           determine if there are more jobs to run. Returns the list
+           of the jobs that still need to be run, as well as the
+           complete list of jobs that need to be collected to get the
+           final answer.
+        """
+# Get the results of the current integration/MINT step
+        self.append_the_results(jobs_to_run,integration_step)
+        self.cross_sect_dict = self.write_res_txt_file(jobs_to_collect,integration_step)
+# Update HTML pages
+        if fixed_order:
+            cross, error = sum_html.make_all_html_results(self, ['%s*' % run_mode])
+        else:
+            name_suffix={'born' :'B' , 'all':'F'}
+            cross, error = sum_html.make_all_html_results(self, ['G%s*' % name_suffix[run_mode]])
+        self.results.add_detail('cross', cross)
+        self.results.add_detail('error', error) 
+# Set-up jobs for the next iteration/MINT step
+        jobs_to_run_new=self.update_jobs_to_run(req_acc,integration_step,jobs_to_run,fixed_order)
+        # if there are no more jobs, we are done!
+# Print summary
+        if (not jobs_to_run_new) and fixed_order:
+            # print final summary of results (for fixed order)
+            scale_pdf_info=self.collect_scale_pdf_info(options,jobs_to_collect)
+            self.print_summary(options,integration_step,mode,scale_pdf_info,done=True)
+            return jobs_to_run_new,jobs_to_collect
+        elif jobs_to_run_new:
+            # print intermediate summary of results
+            scale_pdf_info={}
+            self.print_summary(options,integration_step,mode,scale_pdf_info,done=False)
+        else:
+            # When we are done for (N)LO+PS runs, do not print
+            # anything yet. This will be done after the reweighting
+            # and collection of the events
+            scale_pdf_info={}
+# Prepare for the next integration/MINT step
+        if (not fixed_order) and integration_step+1 == 2 :
+            # next step is event generation (mint_step 2)
+            jobs_to_run_new,jobs_to_collect_new= \
+                    self.check_the_need_to_split(jobs_to_run_new,jobs_to_collect)
+            self.prepare_directories(jobs_to_run_new,mode,fixed_order)
+            self.write_nevents_unweighted_file(jobs_to_collect_new)
+            self.write_nevts_files(jobs_to_run_new)
+        else:
+            self.prepare_directories(jobs_to_run_new,mode,fixed_order)
+            jobs_to_collect_new=jobs_to_collect
+        return jobs_to_run_new,jobs_to_collect_new
+
+
+    def write_nevents_unweighted_file(self,jobs):
+        """writes the nevents_unweighted file in the SubProcesses directory"""
+        content=[]
+        for job in jobs:
+            path=pjoin(job['dirname'].split('/')[-2],job['dirname'].split('/')[-1])
+            lhefile=pjoin(path,'events.lhe')
+            content.append(' %s     %d     %9e     %9e' % \
+                (lhefile.ljust(40),job['nevents'],job['resultABS']*job['wgt_frac'],job['wgt_frac']))
+        with open(pjoin(self.me_dir,'SubProcesses',"nevents_unweighted"),'w') as f:
+            f.write('\n'.join(content)+'\n')
+
+    def write_nevts_files(self,jobs):
+        """write the nevts files in the SubProcesses/P*/G*/ directories"""
+        for job in jobs:
+            with open(pjoin(job['dirname'],'nevts'),'w') as f:
+                f.write('%i\n' % job['nevents'])
+
+    def check_the_need_to_split(self,jobs_to_run,jobs_to_collect):
+        """Looks in the jobs_to_run to see if there is the need to split the
+           event generation step. Updates jobs_to_run and
+           jobs_to_collect to replace the split-job by its
+           splits. Also removes jobs that do not need any events.
+        """
+        nevt_job=self.run_card['nevt_job']
+        if nevt_job > 0:
+            jobs_to_collect_new=copy.copy(jobs_to_collect)
+            for job in jobs_to_run:
+                nevents=job['nevents']
+                if nevents == 0:
+                    jobs_to_collect_new.remove(job)
+                elif nevents > nevt_job:
+                    jobs_to_collect_new.remove(job)
+                    if nevents % nevt_job != 0 :
+                        nsplit=int(nevents/nevt_job)+1
+                    else:
+                        nsplit=int(nevents/nevt_job)
+                    for i in range(1,nsplit+1):
+                        job_new=copy.copy(job)
+                        left_over=nevents % nsplit
+                        if i <= left_over:
+                            job_new['nevents']=int(nevents/nsplit)+1
+                            job_new['wgt_frac']=float(job_new['nevents'])/float(nevents)
+                        else:
+                            job_new['nevents']=int(nevents/nsplit)
+                            job_new['wgt_frac']=float(job_new['nevents'])/float(nevents)
+                        job_new['split']=i
+                        job_new['dirname']=job['dirname']+'_%i' % job_new['split']
+                        jobs_to_collect_new.append(job_new)
+            jobs_to_run_new=copy.copy(jobs_to_collect_new)
+        else:
+            jobs_to_run_new=copy.copy(jobs_to_collect)
+            for job in jobs_to_collect:
+                if job['nevents'] == 0:
+                    jobs_to_run_new.remove(job)
+            jobs_to_collect_new=copy.copy(jobs_to_run_new)
+
+        return jobs_to_run_new,jobs_to_collect_new
+    
+
+    def update_jobs_to_run(self,req_acc,step,jobs,fixed_order=True):
+        """
+        For (N)LO+PS:    determines the number of events and/or the required
+                         accuracy per job.
+        For fixed order: determines which jobs need higher precision and
+                         returns those with the newly requested precision.
+        """
+        err=self.cross_sect_dict['errt']
+        tot=self.cross_sect_dict['xsect']
+        errABS=self.cross_sect_dict['erra']
+        totABS=self.cross_sect_dict['xseca']
+        jobs_new=[]
+        if fixed_order:
+            if req_acc == -1:
+                if step+1 == 1:
+                    npoints = self.run_card['npoints_FO']
+                    niters = self.run_card['niters_FO']
+                    for job in jobs:
+                        job['mint_mode']=-1
+                        job['niters']=niters
+                        job['npoints']=npoints
+                        jobs_new.append(job)
+                elif step+1 == 2:
+                    pass
+                elif step+1 > 2:
+                    raise aMCatNLOError('Cannot determine number of iterations and PS points '+
+                                        'for integration step %i' % step )
+            elif ( req_acc > 0 and err/tot > req_acc*1.2 ) or step == 0:
+                req_accABS=req_acc*abs(tot)/totABS # overal relative required accuracy on ABS Xsec.
+                for job in jobs:
+                    job['mint_mode']=-1
+                    # Determine relative required accuracy on the ABS for this job
+                    job['accuracy']=req_accABS*math.sqrt(totABS/job['resultABS'])
+                    # If already accurate enough, skip running
+                    if job['accuracy'] > job['errorABS']/job['resultABS'] and step != 0:
+                        continue
+                    # Update the number of PS points based on errorABS, ncall and accuracy
+                    itmax_fl=job['niters_done']*math.pow(job['errorABS']/
+                                                         (job['accuracy']*job['resultABS']),2)
+                    if itmax_fl <= 4.0 :
+                        job['niters']=max(int(round(itmax_fl)),2)
+                        job['npoints']=job['npoints_done']*2
+                    elif itmax_fl > 4.0 and itmax_fl <= 16.0 :
+                        job['niters']=4
+                        job['npoints']=int(round(job['npoints_done']*itmax_fl/4.0))*2
+                    else:
+                        if itmax_fl > 100.0 : itmax_fl=50.0
+                        job['niters']=int(round(math.sqrt(itmax_fl)))
+                        job['npoints']=int(round(job['npoints_done']*itmax_fl/
+                                                 round(math.sqrt(itmax_fl))))*2
+                    # Add the job to the list of jobs that need to be run
+                    jobs_new.append(job)
+            return jobs_new
+        elif step+1 <= 2:
+            nevents=self.run_card['nevents']
+            # Total required accuracy for the upper bounding envelope
+            if req_acc<0: 
+                req_acc2_inv=nevents
+            else:
+                req_acc2_inv=1/(req_acc*req_acc)
+            if step+1 == 1 or step+1 == 2 :
+                # determine the req. accuracy for each of the jobs for Mint-step = 1
+                for job in jobs:
+                    accuracy=min(math.sqrt(totABS/(req_acc2_inv*job['resultABS'])),0.2)
+                    job['accuracy']=accuracy
+            if step+1 == 2:
+                # Randomly (based on the relative ABS Xsec of the job) determine the 
+                # number of events each job needs to generate for MINT-step = 2.
+                r=self.get_randinit_seed()
+                random.seed(r)
+                totevts=nevents
+                for job in jobs:
+                    job['nevents'] = 0
+                while totevts :
+                    target = random.random() * totABS
+                    crosssum = 0.
+                    i = 0
+                    while i<len(jobs) and crosssum < target:
+                        job = jobs[i]
+                        crosssum += job['resultABS']
+                        i += 1            
+                    totevts -= 1
+                    i -= 1
+                    jobs[i]['nevents'] += 1
+            for job in jobs:
+                job['mint_mode']=step+1 # next step
+            return jobs
+        else:
+            return []
+
+
+    def get_randinit_seed(self):
+        """ Get the random number seed from the randinit file """
+        with open(pjoin(self.me_dir,"SubProcesses","randinit")) as randinit:
+            # format of the file is "r=%d".
+            iseed = int(randinit.read()[2:])
+        return iseed
+
+
+    def append_the_results(self,jobs,integration_step):
+        """Appends the results for each of the jobs in the job list"""
+        error_found=False
+        for job in jobs:
+            try:
+                if integration_step >= 0 :
+                    with open(pjoin(job['dirname'],'res_%s.dat' % integration_step)) as res_file:
+                        results=res_file.readline().split()
+                else:
+                # should only be here when doing fixed order with the 'only_generation'
+                # option equal to True. Take the results from the final run done.
+                    with open(pjoin(job['dirname'],'res.dat')) as res_file:
+                        results=res_file.readline().split()
+            except IOError:
+                if not error_found:
+                    error_found=True
+                    error_log=[]
+                error_log.append(pjoin(job['dirname'],'log.txt'))
+                continue
+            job['resultABS']=float(results[0])
+            job['errorABS']=float(results[1])
+            job['result']=float(results[2])
+            job['error']=float(results[3])
+            job['niters_done']=int(results[4])
+            job['npoints_done']=int(results[5])
+            job['time_spend']=float(results[6])
+            job['err_percABS'] = job['errorABS']/job['resultABS']*100.
+            job['err_perc'] = job['error']/job['result']*100.
+        if error_found:
+            raise aMCatNLOError('An error occurred during the collection of results.\n' + 
+                   'Please check the .log files inside the directories which failed:\n' +
+                                '\n'.join(error_log)+'\n')
+
+
+
+    def write_res_txt_file(self,jobs,integration_step):
+        """writes the res.txt files in the SubProcess dir"""
+        jobs.sort(key = lambda job: -job['errorABS'])
+        content=[]
+        content.append('\n\nCross-section per integration channel:')
+        for job in jobs:
+            content.append('%(p_dir)20s  %(channel)15s   %(result)10.8e    %(error)6.4e       %(err_perc)6.4f%%  ' %  job)
+        content.append('\n\nABS cross-section per integration channel:')
+        for job in jobs:
+            content.append('%(p_dir)20s  %(channel)15s   %(resultABS)10.8e    %(errorABS)6.4e       %(err_percABS)6.4f%%  ' %  job)
+        totABS=0
+        errABS=0
+        tot=0
+        err=0
+        for job in jobs:
+            totABS+= job['resultABS']
+            errABS+= math.pow(job['errorABS'],2)
+            tot+= job['result']
+            err+= math.pow(job['error'],2)
+        if jobs:
+            content.append('\nTotal ABS and \nTotal: \n                      %10.8e +- %6.4e  (%6.4e%%)\n                      %10.8e +- %6.4e  (%6.4e%%) \n' %\
+                           (totABS, math.sqrt(errABS), math.sqrt(errABS)/totABS *100.,\
+                            tot, math.sqrt(err), math.sqrt(err)/tot *100.))
+        with open(pjoin(self.me_dir,'SubProcesses','res_%s.txt' % integration_step),'w') as res_file:
+            res_file.write('\n'.join(content))
+        randinit=self.get_randinit_seed()
+        return {'xsect':tot,'xseca':totABS,'errt':math.sqrt(err),\
+                'erra':math.sqrt(errABS),'randinit':randinit}
+        
+
+    def collect_scale_pdf_info(self,options,jobs):
+        """read the scale_pdf_dependence.dat files and collects there results"""
+        scale_pdf_info={}
+        if self.run_card['reweight_scale'] or self.run_card['reweight_PDF']:
+            data_files=[]
+            for job in jobs:
+                data_files.append(pjoin(job['dirname'],'scale_pdf_dependence.dat'))
+            scale_pdf_info = self.pdf_scale_from_reweighting(data_files)
+        return scale_pdf_info
+
+
+    def combine_plots_FO(self,folder_name,jobs):
+        """combines the plots and puts then in the Events/run* directory"""
+        devnull = os.open(os.devnull, os.O_RDWR) 
+        if self.analyse_card['fo_analysis_format'].lower() == 'topdrawer':
+            misc.call(['./combine_plots_FO.sh'] + folder_name, \
+                      stdout=devnull, 
+                      cwd=pjoin(self.me_dir, 'SubProcesses'))
+            files.cp(pjoin(self.me_dir, 'SubProcesses', 'MADatNLO.top'),
+                     pjoin(self.me_dir, 'Events', self.run_name))
+            logger.info('The results of this run and the TopDrawer file with the plots' + \
+                        ' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
+        elif self.analyse_card['fo_analysis_format'].lower() == 'hwu':
+            self.combine_plots_HwU(jobs)
+            files.cp(pjoin(self.me_dir, 'SubProcesses', 'MADatNLO.HwU'),
+                     pjoin(self.me_dir, 'Events', self.run_name))
+            files.cp(pjoin(self.me_dir, 'SubProcesses', 'MADatNLO.gnuplot'),
+                     pjoin(self.me_dir, 'Events', self.run_name))
+            try:
+                misc.call(['gnuplot','MADatNLO.gnuplot'],\
+                          stdout=devnull,stderr=devnull,\
+                          cwd=pjoin(self.me_dir, 'Events', self.run_name))
+            except Exception:
+                pass
+            logger.info('The results of this run and the HwU and GnuPlot files with the plots' + \
+                        ' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
+        elif self.analyse_card['fo_analysis_format'].lower() == 'root':
+            misc.call(['./combine_root.sh'] + folder_name, \
+                      stdout=devnull, 
+                      cwd=pjoin(self.me_dir, 'SubProcesses'))
+            files.cp(pjoin(self.me_dir, 'SubProcesses', 'MADatNLO.root'),
+                     pjoin(self.me_dir, 'Events', self.run_name))
+            logger.info('The results of this run and the ROOT file with the plots' + \
+                        ' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
+        else:
+            logger.info('The results of this run' + \
+                        ' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
+
+
+    def combine_plots_HwU(self,jobs):
         """Sums all the plots in the HwU format."""
-        
         logger.debug('Combining HwU plots.')
-        
-        with open(pjoin(self.me_dir,'SubProcesses','dirs.txt')) as dirf:
-            all_histo_paths = dirf.readlines()
-        all_histo_paths = [pjoin(self.me_dir,'SubProcesses',
-                      path.rstrip(),"MADatNLO.HwU") for path in all_histo_paths]
-        
+        all_histo_paths=[]
+        for job in jobs:
+            all_histo_paths.append(pjoin(job['dirname'],"MADatNLO.HwU"))
         histogram_list = histograms.HwUList(all_histo_paths[0])
-
         for histo_path in all_histo_paths[1:]:
             for i, histo in enumerate(histograms.HwUList(histo_path)):
                 # First make sure the plots have the same weight labels and such
@@ -1632,19 +1964,18 @@
         histogram_list.output(pjoin(self.me_dir,'SubProcesses',"MADatNLO"),
                                                              format = 'gnuplot')
 
-    def applgrid_combine(self,cross,error):
+    def applgrid_combine(self,cross,error,jobs):
         """Combines the APPLgrids in all the SubProcess/P*/all_G*/ directories"""
         logger.debug('Combining APPLgrids \n')
         applcomb=pjoin(self.options['applgrid'].rstrip('applgrid-config'),
                                                             'applgrid-combine')
-        with open(pjoin(self.me_dir,'SubProcesses','dirs.txt')) as dirf:
-            all_jobs=dirf.readlines()
+        all_jobs=[]
+        for job in jobs:
+            all_jobs.append(job['dirname'])
         ngrids=len(all_jobs)
-        nobs  =len([name for name in os.listdir(pjoin(self.me_dir,'SubProcesses',
-                          all_jobs[0].rstrip())) if name.endswith("_out.root")])
+        nobs  =len([name for name in os.listdir(all_jobs[0]) if name.endswith("_out.root")])
         for obs in range(0,nobs):
-            gdir = [pjoin(self.me_dir,'SubProcesses',job.rstrip(),"grid_obs_"+
-                                     str(obs)+"_out.root") for job in all_jobs]
+            gdir = [pjoin(job,"grid_obs_"+str(obs)+"_out.root") for job in all_jobs]
             # combine APPLgrids from different channels for observable 'obs'
             if self.run_card["iappl"] == 1:
                 misc.call([applcomb,'-o', pjoin(self.me_dir,"Events",self.run_name,
@@ -1656,8 +1987,7 @@
                         self.run_name,"aMCfast_obs_"+str(obs)+".root"),'-s',
                                   str(unc2_inv),'--weight',str(unc2_inv)]+ gdir)
                 for job in all_jobs:
-                    os.remove(pjoin(self.me_dir,'SubProcesses',job.rstrip(),
-                                               "grid_obs_"+str(obs)+"_in.root"))
+                    os.remove(pjoin(job,"grid_obs_"+str(obs)+"_in.root"))
             else:
                 raise aMCatNLOError('iappl parameter can only be 0, 1 or 2')
             # after combining, delete the original grids
@@ -1698,14 +2028,10 @@
         if not hasattr(self, 'appl_start_grid') or not self.appl_start_grid:
             raise self.InvalidCmd('No APPLgrid name currently defined.'+
                                              'Please provide this information.')             
-        if mode == 'NLO':
-            gdir='all_G'
-        elif mode == 'LO':
-            gdir='born_G'
         #copy the grid to all relevant directories
         for pdir in p_dirs:
             g_dirs = [file for file in os.listdir(pjoin(self.me_dir,
-              "SubProcesses",pdir)) if file.startswith(gdir) and 
+                        "SubProcesses",pdir)) if file.startswith(mode+'_G') and 
                    os.path.isdir(pjoin(self.me_dir,"SubProcesses",pdir, file))]
             for g_dir in g_dirs:
                 for grid in all_grids:
@@ -1714,28 +2040,20 @@
                                                     'grid_obs_'+obs+'_in.root'))
 
 
-    def collect_log_files(self, folders, istep):
+
+
+    def collect_log_files(self, jobs, integration_step):
         """collect the log files and put them in a single, html-friendly file
-        inside the run_... directory"""
-        step_list = ['Grid setting', 'Cross-section computation', 
-                                                             'Event generation']
+        inside the Events/run_.../ directory"""
         log_file = pjoin(self.me_dir, 'Events', self.run_name, 
-                'alllogs_%d.html' % istep)
-        # this keeps track of which step has been computed for which channel
-        channel_dict = {}
-        log_files = []
-        for folder in folders:
-            log_files += glob.glob(pjoin(self.me_dir, 'SubProcesses', 'P*', 
-                                                             folder, 'log.txt'))
+                         'alllogs_%d.html' % integration_step)
+        outfile = open(log_file, 'w')
 
         content = ''
-
-        outfile = open(log_file, 'w')
-
         content += '<HTML><BODY>\n<font face="courier" size=2>'
-        for log in log_files:
-            channel_dict[os.path.dirname(log)] = [istep]
+        for job in jobs:
             # put an anchor
+            log=pjoin(job['dirname'],'log_MINT%s.txt' % integration_step)
             content += '<a name=%s></a>\n' % (os.path.dirname(log).replace(
                                           pjoin(self.me_dir,'SubProcesses'),''))
             # and put some nice header
@@ -1743,7 +2061,7 @@
             content += '<br>LOG file for integration channel %s, %s <br>' % \
                     (os.path.dirname(log).replace(pjoin(self.me_dir,
                                                            'SubProcesses'), ''), 
-                     step_list[istep])
+                     integration_step)
             content += '</font>\n'
             #then just flush the content of the small log inside the big log
             #the PRE tag prints everything verbatim
@@ -1756,53 +2074,80 @@
         outfile.close()
 
 
-    def read_results(self, output, mode):
-        """extract results (cross-section, absolute cross-section and errors)
-        from output, which should be formatted as
-            Found 4 correctly terminated jobs 
-            random seed found in 'randinit' is 33
-            Integrated abs(cross-section)
-            7.94473937e+03 +- 2.9953e+01  (3.7702e-01%)
-             Integrated cross-section
-            6.63392298e+03 +- 3.7669e+01  (5.6782e-01%)
-        for aMC@NLO/aMC@LO, and as
-
-        for NLO/LO
-        The cross_sect_dict is returned"""
-        res = {}
-        if mode in ['aMC@LO', 'aMC@NLO', 'noshower', 'noshowerLO']:
-            pat = re.compile(\
-'''Found (\d+) correctly terminated jobs 
-random seed found in 'randinit' is (\d+)
-Integrated abs\(cross-section\)
-\s*(\d+\.\d+e[+-]\d+) \+\- (\d+\.\d+e[+-]\d+)  \((\d+\.\d+e[+-]\d+)\%\)
-Integrated cross-section
-\s*(\-?\d+\.\d+e[+-]\d+) \+\- (\d+\.\d+e[+-]\d+)  \((\-?\d+\.\d+e[+-]\d+)\%\)''')
-        else:
-            pat = re.compile(\
-'''Found (\d+) correctly terminated jobs 
-\s*(\-?\d+\.\d+e[+-]\d+) \+\- (\d+\.\d+e[+-]\d+)  \((\-?\d+\.\d+e[+-]\d+)\%\)''')
-            pass
-
-        match = re.search(pat, output[0])
-        if not match or output[1]:
-            logger.info('Return code of the event collection: '+str(output[1]))
-            logger.info('Output of the event collection:\n'+output[0])
-            raise aMCatNLOError('An error occurred during the collection of results.\n' + 
-            'Please check the .log files inside the directories which failed.')
-#        if int(match.groups()[0]) != self.njobs:
-#            raise aMCatNLOError('Not all jobs terminated successfully')
-        if mode in ['aMC@LO', 'aMC@NLO', 'noshower', 'noshowerLO']:
-            return {'randinit' : int(match.groups()[1]),
-                    'xseca' : float(match.groups()[2]),
-                    'erra' : float(match.groups()[3]),
-                    'xsect' : float(match.groups()[5]),
-                    'errt' : float(match.groups()[6])}
-        else:
-            return {'xsect' : float(match.groups()[1]),
-                    'errt' : float(match.groups()[2])}
-
-    def print_summary(self, options, step, mode, scale_pdf_info={}):
+    def finalise_run_FO(self,folder_name,jobs):
+        """Combine the plots and put the res*.txt files in the Events/run.../ folder."""
+        # Copy the res_*.txt files to the Events/run* folder
+        res_files=glob.glob(pjoin(self.me_dir, 'SubProcesses', 'res_*.txt'))
+        for res_file in res_files:
+            files.mv(res_file,pjoin(self.me_dir, 'Events', self.run_name))
+        # Collect the plots and put them in the Events/run* folder
+        self.combine_plots_FO(folder_name,jobs)
+        # If doing the applgrid-stuff, also combine those grids
+        # and put those in the Events/run* folder
+        if self.run_card['iappl'] != 0:
+            cross=self.cross_sect_dict['xsect']
+            error=self.cross_sect_dict['errt']
+            self.applgrid_combine(cross,error)
+
+
+    def setup_cluster_or_multicore(self):
+        """setup the number of cores for multicore, and the cluster-type for cluster runs"""
+        if self.cluster_mode == 1:
+            cluster_name = self.options['cluster_type']
+            self.cluster = cluster.from_name[cluster_name](**self.options)
+        if self.cluster_mode == 2:
+            try:
+                import multiprocessing
+                if not self.nb_core:
+                    try:
+                        self.nb_core = int(self.options['nb_core'])
+                    except TypeError:
+                        self.nb_core = multiprocessing.cpu_count()
+                logger.info('Using %d cores' % self.nb_core)
+            except ImportError:
+                self.nb_core = 1
+                logger.warning('Impossible to detect the number of cores => Using One.\n'+
+                        'Use set nb_core X in order to set this number and be able to'+
+                        'run in multicore.')
+
+            self.cluster = cluster.MultiCore(**self.options)
+
+
+    def clean_previous_results(self,options,p_dirs,folder_name):
+        """Clean previous results.
+             o.  If doing only the reweighting step, do not delete anything and return directlty.
+             o.  Always remove all the G*_* files (from split event generation).
+             o.  Remove the G* (or born_G* or all_G*) only when NOT doing only_generation or reweight_only."""
+        if options['reweightonly']:
+            return
+        if not options['only_generation']:
+            self.update_status('Cleaning previous results', level=None)
+        for dir in p_dirs:
+            #find old folders to be removed
+            for obj in folder_name:
+                # list all the G* (or all_G* or born_G*) directories
+                to_rm = [file for file in \
+                             os.listdir(pjoin(self.me_dir, 'SubProcesses', dir)) \
+                             if file.startswith(obj[:-1]) and \
+                            (os.path.isdir(pjoin(self.me_dir, 'SubProcesses', dir, file)) or \
+                             os.path.exists(pjoin(self.me_dir, 'SubProcesses', dir, file)))] 
+                # list all the G*_* directories (from split event generation)
+                to_always_rm = [file for file in \
+                             os.listdir(pjoin(self.me_dir, 'SubProcesses', dir)) \
+                             if file.startswith(obj[:-1]) and
+                             '_' in file and not '_G' in file and \
+                            (os.path.isdir(pjoin(self.me_dir, 'SubProcesses', dir, file)) or \
+                             os.path.exists(pjoin(self.me_dir, 'SubProcesses', dir, file)))]
+
+                if not options['only_generation']:
+                    to_always_rm.extend(to_rm)
+                    if os.path.exists(pjoin(self.me_dir, 'SubProcesses', dir,'MadLoop5_resources.tar.gz')):
+                        to_always_rm.append(pjoin(self.me_dir, 'SubProcesses', dir,'MadLoop5_resources.tar.gz'))
+                files.rm([pjoin(self.me_dir, 'SubProcesses', dir, d) for d in to_always_rm])
+        return
+
+
+    def print_summary(self, options, step, mode, scale_pdf_info={}, done=True):
         """print a summary of the results contained in self.cross_sect_dict.
         step corresponds to the mintMC step, if =2 (i.e. after event generation)
         some additional infos are printed"""
@@ -1813,25 +2158,35 @@
             if line.startswith('generate') or line.startswith('add process'):
                 process = process+(line.replace('generate ', '')).replace('add process ','')+' ; '
         lpp = {0:'l', 1:'p', -1:'pbar'}
-        proc_info = '\n      Process %s\n      Run at %s-%s collider (%s + %s GeV)' % \
-        (process[:-3], lpp[self.run_card['lpp1']], lpp[self.run_card['lpp2']], 
-                self.run_card['ebeam1'], self.run_card['ebeam2'])
-        
+        if self.ninitial == 1:
+            proc_info = '\n      Process %s' % process[:-3]
+        else:
+            proc_info = '\n      Process %s\n      Run at %s-%s collider (%s + %s GeV)' % \
+                (process[:-3], lpp[self.run_card['lpp1']], lpp[self.run_card['lpp2']], 
+                 self.run_card['ebeam1'], self.run_card['ebeam2'])
+
+        if self.ninitial == 1:
+            self.cross_sect_dict['unit']='GeV'
+            self.cross_sect_dict['xsec_string']='(Partial) decay width'
+            self.cross_sect_dict['axsec_string']='(Partial) abs(decay width)'
+        else:
+            self.cross_sect_dict['unit']='pb'
+            self.cross_sect_dict['xsec_string']='Total cross-section'
+            self.cross_sect_dict['axsec_string']='Total abs(cross-section)'
         # Gather some basic statistics for the run and extracted from the log files.
         if mode in ['aMC@NLO', 'aMC@LO', 'noshower', 'noshowerLO']: 
             log_GV_files =  glob.glob(pjoin(self.me_dir, \
                                     'SubProcesses', 'P*','G*','log_MINT*.txt'))
-            all_log_files = glob.glob(pjoin(self.me_dir, \
-                                          'SubProcesses', 'P*','G*','log*.txt'))
+            all_log_files = log_GV_files
         elif mode == 'NLO':
             log_GV_files =  glob.glob(pjoin(self.me_dir, \
-                                    'SubProcesses', 'P*','all_G*','log*.txt'))
-            all_log_files = sum([glob.glob(pjoin(self.me_dir,'SubProcesses', 'P*',
-              '%sG*'%foldName,'log*.txt')) for foldName in ['all_']],[])
+                                    'SubProcesses', 'P*','all_G*','log_MINT*.txt'))
+            all_log_files = log_GV_files
+
         elif mode == 'LO':
             log_GV_files = ''
-            all_log_files = sum([glob.glob(pjoin(self.me_dir,'SubProcesses', 'P*',
-              '%sG*'%foldName,'log*.txt')) for foldName in ['born_']],[])
+            all_log_files = glob.glob(pjoin(self.me_dir, \
+                                    'SubProcesses', 'P*','born_G*','log_MINT*.txt'))
         else:
             raise aMCatNLOError, 'Running mode %s not supported.'%mode
             
@@ -1843,13 +2198,13 @@
             if step != 2:
                 message = status[step] + '\n\n      Intermediate results:' + \
                     ('\n      Random seed: %(randinit)d' + \
-                     '\n      Total cross-section:      %(xsect)8.3e +- %(errt)6.1e pb' + \
-                     '\n      Total abs(cross-section): %(xseca)8.3e +- %(erra)6.1e pb \n') \
+                     '\n      %(xsec_string)s:      %(xsect)8.3e +- %(errt)6.1e %(unit)s' + \
+                     '\n      %(axsec_string)s: %(xseca)8.3e +- %(erra)6.1e %(unit)s \n') \
                      % self.cross_sect_dict
             else:
         
                 message = '\n      ' + status[step] + proc_info + \
-                          '\n      Total cross-section: %(xsect)8.3e +- %(errt)6.1e pb' % \
+                          '\n      %(xsec_string)s: %(xsect)8.3e +- %(errt)6.1e %(unit)s' % \
                         self.cross_sect_dict
 
                 if self.run_card['nevents']>=10000 and self.run_card['reweight_scale']:
@@ -1874,18 +2229,22 @@
                          misc.format_timer(time.time()-self.start_time))
 
         elif mode in ['NLO', 'LO']:
-            status = ['Results after grid setup (cross-section is non-physical):',
+            status = ['Results after grid setup:','Current results:',
                       'Final results and run summary:']
-            if step == 0:
-                message = '\n      ' + status[step] + \
-                     '\n      Total cross-section:      %(xsect)8.3e +- %(errt)6.1e pb' % \
-                             self.cross_sect_dict
-            elif step == 1:
-                message = '\n      ' + status[step] + proc_info + \
-                     '\n      Total cross-section:      %(xsect)8.3e +- %(errt)6.1e pb' % \
+            if (not done) and (step == 0):
+                message = '\n      ' + status[0] + \
+                     '\n      %(xsec_string)s:      %(xsect)8.3e +- %(errt)6.1e %(unit)s' % \
+                             self.cross_sect_dict
+            elif not done:
+                message = '\n      ' + status[1] + \
+                     '\n      %(xsec_string)s:      %(xsect)8.3e +- %(errt)6.1e %(unit)s' % \
+                             self.cross_sect_dict
+            elif done:
+                message = '\n      ' + status[2] + proc_info + \
+                     '\n      %(xsec_string)s:      %(xsect)8.3e +- %(errt)6.1e %(unit)s' % \
                              self.cross_sect_dict
                 if self.run_card['reweight_scale']:
-                    if int(self.run_card['ickkw'])!=-1:
+                    if self.run_card['ickkw'] != -1:
                         message = message + \
                             ('\n      Ren. and fac. scale uncertainty: +%0.1f%% -%0.1f%%') % \
                             (scale_pdf_info['scale_upp'], scale_pdf_info['scale_low'])
@@ -1898,7 +2257,7 @@
                         ('\n      PDF uncertainty: +%0.1f%% -%0.1f%%') % \
                         (scale_pdf_info['pdf_upp'], scale_pdf_info['pdf_low'])
         
-        if (mode in ['NLO', 'LO'] and step!=1) or \
+        if (mode in ['NLO', 'LO'] and not done) or \
            (mode in ['aMC@NLO', 'aMC@LO', 'noshower', 'noshowerLO'] and step!=2):
             logger.info(message+'\n')
             return
@@ -2359,7 +2718,6 @@
         scale_pdf_info={}
         if self.run_card['reweight_scale'] or self.run_card['reweight_PDF'] :
             scale_pdf_info = self.run_reweight(options['reweightonly'])
-
         self.update_status('Collecting events', level='parton', update_results=True)
         misc.compile(['collect_events'], 
                     cwd=pjoin(self.me_dir, 'SubProcesses'))
@@ -2383,6 +2741,10 @@
         misc.gzip(pjoin(self.me_dir, 'SubProcesses', filename), stdout=evt_file)
         if not options['reweightonly']:
             self.print_summary(options, 2, mode, scale_pdf_info)
+            res_files=glob.glob(pjoin(self.me_dir, 'SubProcesses', 'res*.txt'))
+            for res_file in res_files:
+                files.mv(res_file,pjoin(self.me_dir, 'Events', self.run_name))
+
         logger.info('The %s file has been generated.\n' % (evt_file))
         self.results.add_detail('nb_event', nevents)
         self.update_status('Events generated', level='parton', update_results=True)
@@ -2403,9 +2765,9 @@
 
         #check that the number of split event files divides the number of
         # events, otherwise set it to 1
-        if int(int(self.banner.get_detail('run_card', 'nevents')) / \
+        if int(self.banner.get_detail('run_card', 'nevents') / \
                 self.shower_card['nsplit_jobs']) * self.shower_card['nsplit_jobs'] \
-                != int(self.banner.get_detail('run_card', 'nevents')):
+                != self.banner.get_detail('run_card', 'nevents'):
             logger.warning(\
                 'nsplit_jobs in the shower card is not a divisor of the number of events.\n' + \
                 'Setting it to 1.')
@@ -2413,7 +2775,7 @@
 
         # don't split jobs if the user asks to shower only a part of the events
         if self.shower_card['nevents'] > 0 and \
-           self.shower_card['nevents'] < int(self.banner.get_detail('run_card', 'nevents')) and \
+           self.shower_card['nevents'] < self.banner.get_detail('run_card', 'nevents') and \
            self.shower_card['nsplit_jobs'] != 1:
             logger.warning(\
                 'Only a part of the events will be showered.\n' + \
@@ -2503,6 +2865,18 @@
         
         mcatnlo_log = pjoin(self.me_dir, 'mcatnlo.log')
         self.update_status('Compiling MCatNLO for %s...' % shower, level='shower') 
+
+        
+        # libdl may be needded for pythia 82xx
+        if shower == 'PYTHIA8' and not \
+           os.path.exists(pjoin(self.options['pythia8_path'], 'xmldoc')) and \
+           'dl' not in self.shower_card['extralibs'].split():
+            # 'dl' has to be linked with the extralibs
+            self.shower_card['extralibs'] += ' dl'
+            logger.warning("'dl' was added to extralibs from the shower_card.dat.\n" + \
+                          "It is needed for the correct running of PY8.2xx.\n" + \
+                          "If this library cannot be found on your system, a crash will occur.")
+
         misc.call(['./MCatNLO_MadFKS.inputs'], stdout=open(mcatnlo_log, 'w'),
                     stderr=open(mcatnlo_log, 'w'), 
                     cwd=pjoin(self.me_dir, 'MCatNLO'))
@@ -2554,10 +2928,10 @@
         # special treatment for pythia8
             files.mv(pjoin(self.me_dir, 'MCatNLO', 'Pythia8.cmd'), rundir)
             files.mv(pjoin(self.me_dir, 'MCatNLO', 'Pythia8.exe'), rundir)
-            if os.path.exists(pjoin(self.options['pythia8_path'], 'xmldoc')):
+            if os.path.exists(pjoin(self.options['pythia8_path'], 'xmldoc')): # this is PY8.1xxx
                 files.ln(pjoin(self.options['pythia8_path'], 'examples', 'config.sh'), rundir)
                 files.ln(pjoin(self.options['pythia8_path'], 'xmldoc'), rundir)
-            else:
+            else: # this is PY8.2xxx
                 files.ln(pjoin(self.options['pythia8_path'], 'share/Pythia8/xmldoc'), rundir)
         #link the hwpp exe in the rundir
         if shower == 'HERWIGPP':
@@ -2987,8 +3361,8 @@
         init_dict = self.get_init_dict(evt_file)
 
         if nevents < 0 or \
-           nevents > int(self.banner.get_detail('run_card', 'nevents')):
-            nevents = int(self.banner.get_detail('run_card', 'nevents'))
+           nevents > self.banner.get_detail('run_card', 'nevents'):
+            nevents = self.banner.get_detail('run_card', 'nevents')
 
         nevents = nevents / self.shower_card['nsplit_jobs']
 
@@ -3000,7 +3374,7 @@
 
         content = 'EVPREFIX=%s\n' % pjoin(os.path.split(evt_file)[1])
         content += 'NEVENTS=%d\n' % nevents
-        content += 'NEVENTS_TOT=%d\n' % (int(self.banner.get_detail('run_card', 'nevents')) /\
+        content += 'NEVENTS_TOT=%d\n' % (self.banner.get_detail('run_card', 'nevents') /\
                                              self.shower_card['nsplit_jobs'])
         content += 'MCMODE=%s\n' % shower
         content += 'PDLABEL=%s\n' % pdlabel
@@ -3113,7 +3487,7 @@
 
 
     def run_reweight(self, only):
-        """runs the reweight_xsec_events eecutables on each sub-event file generated
+        """runs the reweight_xsec_events executables on each sub-event file generated
         to compute on the fly scale and/or PDF uncertainities"""
         logger.info('   Doing reweight')
 
@@ -3231,7 +3605,7 @@
                 scale_pdf_info['scale_low'] = 0.0
 
         # get the pdf uncertainty in percent (according to the Hessian method)
-        lhaid=int(self.run_card['lhaid'])
+        lhaid=self.run_card['lhaid']
         pdf_upp=0.0
         pdf_low=0.0
         if lhaid <= 90000:
@@ -3246,7 +3620,6 @@
                 else:
                     scale_pdf_info['pdf_upp'] = 0.0
                     scale_pdf_info['pdf_low'] = 0.0
-
         else:
             # use Gaussian method (NNPDF)
             pdf_stdev=0.0
@@ -3263,7 +3636,6 @@
 
     def wait_for_complete(self, run_type):
         """this function waits for jobs on cluster to complete their run."""
-
         starttime = time.time()
         #logger.info('     Waiting for submitted jobs to complete')
         update_status = lambda i, r, f: self.update_status((i, r, f, run_type), 
@@ -3276,29 +3648,15 @@
 
     def run_all(self, job_dict, arg_list, run_type='monitor', split_jobs = False):
         """runs the jobs in job_dict (organized as folder: [job_list]), with arguments args"""
-        njob_split = 0
         self.ijob = 0
-
-        #  this is to keep track, if splitting evt generation, of the various 
-        # folders/args in order to resubmit the jobs if some of them fail
-        self.split_folders = {}
-
         if run_type != 'shower':
             self.njobs = sum(len(jobs) for jobs in job_dict.values()) * len(arg_list)
             for args in arg_list:
                 for Pdir, jobs in job_dict.items():
                     for job in jobs:
-                        if not split_jobs:
-                            self.run_exe(job, args, run_type, cwd=pjoin(self.me_dir, 'SubProcesses', Pdir) )
-                        else:
-                            for n in self.find_jobs_to_split(Pdir, job, args[1]):
-                                self.run_exe(job, args + [n], run_type, cwd=pjoin(self.me_dir, 'SubProcesses', Pdir) )
-                                njob_split += 1
-                        # print some statistics if running serially
+                        self.run_exe(job, args, run_type, cwd=pjoin(self.me_dir, 'SubProcesses', Pdir) )
             if self.cluster_mode == 2:
                 time.sleep(1) # security to allow all jobs to be launched
-            if njob_split > 0:
-                self.njobs = njob_split
         else:
             self.njobs = len(arg_list)
             for args in arg_list:
@@ -3309,37 +3667,27 @@
 
 
 
-    def check_event_files(self):
+    def check_event_files(self,jobs):
         """check the integrity of the event files after splitting, and resubmit 
         those which are not nicely terminated"""
-        to_resubmit = []
-        for dir in self.split_folders.keys():
+        jobs_to_resubmit = []
+        for job in jobs:
             last_line = ''
             try:
                 last_line = subprocess.Popen(
-                        ['tail', '-n1', pjoin(dir, 'events.lhe')], \
+                        ['tail', '-n1', pjoin(job['dirname'], 'events.lhe')], \
                     stdout = subprocess.PIPE).stdout.read().strip()
             except IOError:
                 pass
-
             if last_line != "</LesHouchesEvents>":
-                to_resubmit.append(dir)
-
+                jobs_to_resubmit.append(job)
         self.njobs = 0
-        if to_resubmit:
+        if jobs_to_resubmit:
             run_type = 'Resubmitting broken jobs'
             logger.info('Some event files are broken, corresponding jobs will be resubmitted.')
-            logger.debug('Resubmitting\n' + '\n'.join(to_resubmit) + '\n')
-            for dir in to_resubmit:
-                files.rm([dir])
-                job = self.split_folders[dir][0]
-                args = self.split_folders[dir][1:]
-                run_type = 'monitor'
-                cwd = os.path.split(dir)[0]
-                self.run_exe(job, args, run_type, cwd=cwd )
-                self.njobs +=1
-
-            self.wait_for_complete(run_type)
+            for job in jobs_to_resubmit:
+                logger.debug('Resubmitting ' + job['dirname'] + '\n')
+            self.run_all_jobs(jobs_to_resubmit,2,fixed_order=False)
 
 
     def find_jobs_to_split(self, pdir, job, arg):
@@ -3412,16 +3760,16 @@
             # the 'standard' amcatnlo job
             # check if args is a list of string 
             if type(args[0]) == str:
-                input_files, output_files, required_output, args = self.getIO_ajob(exe,cwd, args)
+                input_files, output_files, required_output, args = self.getIO_ajob(exe,cwd,args)
                 #submitting
                 self.cluster.submit2(exe, args, cwd=cwd, 
                              input_files=input_files, output_files=output_files,
                              required_output=required_output)
 
-                # keep track of folders and arguments for splitted evt gen
-                subfolder=output_files[-1].split('/')[0]
-                if len(args) == 4 and '_' in subfolder:
-                    self.split_folders[pjoin(cwd,subfolder)] = [exe] + args
+#                # keep track of folders and arguments for splitted evt gen
+#                subfolder=output_files[-1].split('/')[0]
+#                if len(args) == 4 and '_' in subfolder:
+#                    self.split_folders[pjoin(cwd,subfolder)] = [exe] + args
 
         elif 'shower' in exe:
             # a shower job
@@ -3487,7 +3835,6 @@
         # use local disk if possible => need to stands what are the 
         # input/output files
         
-        keep_fourth_arg = False
         output_files = []
         required_output = []
         input_files = [pjoin(self.me_dir, 'SubProcesses', 'randinit'),
@@ -3509,91 +3856,58 @@
             input_files.append(pjoin(cwd, 'OLE_order.olc'))
 
         # File for the loop (might not be present if MadLoop is not used)
-        if os.path.exists(pjoin(cwd,'MadLoop5_resources')) and \
+        if os.path.exists(pjoin(cwd,'MadLoop5_resources.tar.gz')) and \
                                             cluster.need_transfer(self.options):
             input_files.append(pjoin(cwd, 'MadLoop5_resources.tar.gz'))
+        elif os.path.exists(pjoin(cwd,'MadLoop5_resources')) and \
+                                            cluster.need_transfer(self.options):
             tf=tarfile.open(pjoin(cwd,'MadLoop5_resources.tar.gz'),'w:gz',
                                                            dereference=True)
             tf.add(pjoin(cwd,'MadLoop5_resources'),arcname='MadLoop5_resources')
             tf.close()
-
-        Ire = re.compile("for i in ([\d\s]*) ; do")
-        try : 
-            fsock = open(exe)
-        except IOError:
-            fsock = open(pjoin(cwd,exe))
-        text = fsock.read()
-        data = Ire.findall(text)
-        subdir = ' '.join(data).split()
+            input_files.append(pjoin(cwd, 'MadLoop5_resources.tar.gz'))
                
-        if args[0] == '0':
+        if args[1] == 'born' or args[1] == 'all':
             # MADEVENT MINT FO MODE
             input_files.append(pjoin(cwd, 'madevent_mintFO'))
-            input_files.append(pjoin(self.me_dir, 'SubProcesses','madin.%s' % args[1]))
-            #j=$2\_G$i
-            for i in subdir:
-                current = '%s_G%s' % (args[1],i)
-                if os.path.exists(pjoin(cwd,current)):
-                    input_files.append(pjoin(cwd, current))
-                output_files.append(current)
+            if args[2] == '0':
+                current = '%s_G%s' % (args[1],args[0])
+            else:
+                current = '%s_G%s_%s' % (args[1],args[0],args[2])
+            if os.path.exists(pjoin(cwd,current)):
+                input_files.append(pjoin(cwd, current))
+            output_files.append(current)
 
-                required_output.append('%s/results.dat' % current)
-                required_output.append('%s/log.txt' % current)
-                required_output.append('%s/mint_grids' % current)
-                required_output.append('%s/grid.MC_integer' % current)
-                if len(args) == 4:
-                    required_output.append('%s/scale_pdf_dependence.dat' % current)
-                    args[2] = '-1'
-                    # use a grid train on another part
-                    base = '%s_G%s' % (args[3],i)
-                    if args[0] == '0':
-                        to_move = ['grid.MC_integer','mint_grids']
-                    elif args[0] == '1':
-                        to_move = ['mint_grids', 'grid.MC_integer']
-                    else: 
-                        to_move  = []
-                    if self.run_card['iappl'] == 2:
-                        for grid in glob.glob(pjoin(cwd,base,'grid_obs_*_in.root')):
-                            to_move.append(grid)
-                    if not os.path.exists(pjoin(cwd,current)):
-                        os.mkdir(pjoin(cwd,current))
-                        input_files.append(pjoin(cwd, current))
-                    for name in to_move:
-                        files.cp(pjoin(cwd,base, name), 
-                                        pjoin(cwd,current))
-                    files.cp(pjoin(cwd,base, 'grid.MC_integer'), 
-                                        pjoin(cwd,current))
+            required_output.append('%s/results.dat' % current)
+            required_output.append('%s/res_%s.dat' % (current,args[3]))
+            required_output.append('%s/log_MINT%s.txt' % (current,args[3]))
+            required_output.append('%s/mint_grids' % current)
+            required_output.append('%s/grid.MC_integer' % current)
+            if args[3] != '0':
+                required_output.append('%s/scale_pdf_dependence.dat' % current)
                             
-        elif args[0] == '2':
+        elif args[1] == 'F' or args[1] == 'B':
             # MINTMC MODE
             input_files.append(pjoin(cwd, 'madevent_mintMC'))
-            if args[2] in ['0','2']:
-                input_files.append(pjoin(self.me_dir, 'SubProcesses','madinMMC_%s.2' % args[1]))
-
-            for i in subdir:
-                current = 'G%s%s' % (args[1], i)
-                if os.path.exists(pjoin(cwd,current)):
-                    input_files.append(pjoin(cwd, current))
-                output_files.append(current)
-                if len(args) == 4 and args[3] in ['H','S','V','B','F']:
-                    # use a grid train on another part
-                    base = '%s_%s' % (args[3],i)
-                    files.ln(pjoin(cwd,base,'mint_grids'), name = 'preset_mint_grids', 
-                                            starting_dir=pjoin(cwd,current))
-                    files.ln(pjoin(cwd,base,'grid.MC_integer'), 
-                                          starting_dir=pjoin(cwd,current))
-                elif len(args) ==4:
-                    keep_fourth_arg = True
-                    # this is for the split event generation
-                    output_files.append('G%s%s_%s' % (args[1], i, args[3]))
-                    required_output.append('G%s%s_%s/log_MINT%s.txt' % (args[1], i, args[3],args[2]))
-
-                else:
-                    required_output.append('%s/log_MINT%s.txt' % (current,args[2]))
-                if args[2] in ['0','1']:
-                    required_output.append('%s/results.dat' % current)
-                if args[2] == '1':
-                    output_files.append('%s/results.dat' % current)
+
+            if args[2] == '0':
+                current = 'G%s%s' % (args[1],args[0])
+            else:
+                current = 'G%s%s_%s' % (args[1],args[0],args[2])
+            if os.path.exists(pjoin(cwd,current)):
+                input_files.append(pjoin(cwd, current))
+            output_files.append(current)
+            if args[2] > '0':
+                # this is for the split event generation
+                output_files.append('G%s%s_%s' % (args[1], args[0], args[2]))
+                required_output.append('G%s%s_%s/log_MINT%s.txt' % (args[1],args[0],args[2],args[3]))
+
+            else:
+                required_output.append('%s/log_MINT%s.txt' % (current,args[3]))
+            if args[3] in ['0','1']:
+                required_output.append('%s/results.dat' % current)
+            if args[3] == '1':
+                output_files.append('%s/results.dat' % current)
 
         else:
             raise aMCatNLOError, 'not valid arguments: %s' %(', '.join(args))
@@ -3601,73 +3915,9 @@
         #Find the correct PDF input file
         pdfinput = self.get_pdf_input_filename()
         if os.path.exists(pdfinput):
-            input_files.append(pdfinput)
-
-        if len(args) == 4 and not keep_fourth_arg:
-            args = args[:3]
-            
+            input_files.append(pdfinput)            
         return input_files, output_files, required_output,  args
-            
-    def write_madinMMC_file(self, path, run_mode, mint_mode):
-        """writes the madinMMC_?.2 file"""
-        #check the validity of the arguments
-        run_modes = ['born', 'virt', 'novi', 'all', 'viSB', 'novB']
-        if run_mode not in run_modes:
-            raise aMCatNLOError('%s is not a valid mode for run. Please use one of the following: %s' \
-                    % (run_mode, ', '.join(run_modes)))
-        mint_modes = [0, 1, 2]
-        if mint_mode not in mint_modes:
-            raise aMCatNLOError('%s is not a valid mode for mintMC. Please use one of the following: %s' \
-                    % (mint_mode, ', '.join(mint_modes)))
-        if run_mode in ['born']:
-            name_suffix = 'B'
-        elif run_mode in ['virt', 'viSB']:
-            name_suffix = 'V'
-        else:
-            name_suffix = 'F'
-
-        content = \
-"""-1 12      ! points, iterations
-0.03       ! desired fractional accuracy
-1 -0.1     ! alpha, beta for Gsoft
--1 -0.1    ! alpha, beta for Gazi
-1          ! Suppress amplitude (0 no, 1 yes)?
-1          ! Exact helicity sum (0 yes, n = number/event)?
-1          ! Enter Configuration Number:
-%1d          ! MINT imode: 0 to set-up grids, 1 to perform integral, 2 generate events
-1 1 1      ! if imode is 1: Folding parameters for xi_i, phi_i and y_ij
-%s        ! all, born, real, virt
-""" \
-                    % (mint_mode, run_mode)
-        file = open(pjoin(path, 'madinMMC_%s.2' % name_suffix), 'w')
-        file.write(content)
-        file.close()
-
-    def write_madin_file(self, path, run_mode, vegas_mode, npoints, niters, accuracy='0'):
-        """writes the madin.run_mode file"""
-        #check the validity of the arguments
-        run_modes = ['born', 'virt', 'novi', 'all', 'viSB', 'novB', 'grid']
-        if run_mode not in run_modes:
-            raise aMCatNLOError('%s is not a valid mode for run. Please use one of the following: %s' \
-                    % (run_mode, ', '.join(run_modes)))
-        name_suffix = run_mode
-
-        content = \
-"""%s %s  ! points, iterations
-%s ! accuracy
-2 ! 0 fixed grid 2 adjust
-1 ! 1 suppress amp, 0 doesnt
-1 ! 0 for exact hel sum
-1 ! hel configuration numb
-'test'
-1 ! 1 to save grids
-%s ! 0 to exclude, 1 for new run, 2 to restart, 3 to reset w/ keeping grid
-%s        ! all, born, real, virt
-""" \
-                    % (npoints,niters,accuracy,vegas_mode,run_mode)
-        file = open(pjoin(path, 'madin.%s' % name_suffix), 'w')
-        file.write(content)
-        file.close()
+
 
     def compile(self, mode, options):
         """compiles aMC@NLO to compute either NLO or NLO matched to shower, as
@@ -3730,10 +3980,10 @@
 
             self.link_lhapdf(libdir, [pjoin('SubProcesses', p) for p in p_dirs])
             pdfsetsdir = self.get_lhapdf_pdfsetsdir()
-            lhaid_list = [int(self.run_card['lhaid'])]
+            lhaid_list = [self.run_card['lhaid']]
             if self.run_card['reweight_PDF']:
-                lhaid_list.append(int(self.run_card['PDF_set_min']))
-                lhaid_list.append(int(self.run_card['PDF_set_max']))
+                lhaid_list.append(self.run_card['PDF_set_min'])
+                lhaid_list.append(self.run_card['PDF_set_max'])
             self.copy_lhapdf_set(lhaid_list, pdfsetsdir)
 
         else:
@@ -4045,18 +4295,23 @@
         void = 'NOT INSTALLED'
         switch_order = ['order', 'fixed_order', 'shower','madspin', 'reweight']
         switch_default = {'order': 'NLO', 'fixed_order': 'OFF', 'shower': void,
-                  'madspin': void}
+                  'madspin': void,'reweight':'OFF'}
         if not switch:
             switch = switch_default
         else:
             switch.update(dict((k,value) for k,v in switch_default.items() if k not in switch))
-
         default_switch = ['ON', 'OFF']
+        
+
         allowed_switch_value = {'order': ['LO', 'NLO'],
                                 'fixed_order': default_switch,
                                 'shower': default_switch,
                                 'madspin': default_switch,
                                 'reweight': default_switch}
+
+            
+            
+        
         
         description = {'order':  'Perturbative order of the calculation:',
                        'fixed_order': 'Fixed order (no event generation and no MC@[N]LO matching):',
@@ -4072,22 +4327,35 @@
         special_values = ['LO', 'NLO', 'aMC@NLO', 'aMC@LO', 'noshower', 'noshowerLO']
 
         assign_switch = lambda key, value: switch.__setitem__(key, value if switch[key] != void else void )
-        
+
+        if self.proc_characteristics['ninitial'] == 1:
+            switch['fixed_order'] = 'ON'
+            switch['shower'] = 'Not available for decay'
+            switch['madspin'] = 'Not available for decay'
+            switch['reweight'] = 'Not available for decay'
+            allowed_switch_value['fixed_order'] = ['ON']
+            allowed_switch_value['shower'] = ['OFF']
+            allowed_switch_value['madspin'] = ['OFF']
+            allowed_switch_value['reweight'] = ['OFF']
+            available_mode = ['0','1']
+            special_values = ['LO', 'NLO']
+        else: 
+            # Init the switch value according to the current status
+            available_mode = ['0', '1', '2','3']
 
         if mode == 'auto': 
             mode = None
         if not mode and (options['parton'] or options['reweightonly']):
             mode = 'noshower'         
         
-        # Init the switch value according to the current status
-        available_mode = ['0', '1', '2']
-        available_mode.append('3')
-        if os.path.exists(pjoin(self.me_dir, 'Cards', 'shower_card.dat')):
-            switch['shower'] = 'ON'
-        else:
-            switch['shower'] = 'OFF'
+
+        if '3' in available_mode:
+            if os.path.exists(pjoin(self.me_dir, 'Cards', 'shower_card.dat')):
+                switch['shower'] = 'ON'
+            else:
+                switch['shower'] = 'OFF' 
                 
-        if not aMCatNLO or self.options['mg5_path']:
+        if (not aMCatNLO or self.options['mg5_path']) and '3' in available_mode:
             available_mode.append('4')
             if os.path.exists(pjoin(self.me_dir,'Cards','madspin_card.dat')):
                 switch['madspin'] = 'ON'
@@ -4102,8 +4370,7 @@
             else:
                 switch['reweight'] = 'Not available (requires NumPy)'
 
-        
-        if 'do_reweight' in options and options['do_reweight']:
+        if 'do_reweight' in options and options['do_reweight'] and '3' in available_mode:
             if switch['reweight'] == "OFF":
                 switch['reweight'] = "ON"
             elif switch['reweight'] != "ON":
@@ -4113,12 +4380,12 @@
                 switch['madspin'] = 'ON'
             elif switch['madspin'] != "ON":
                 logger.critical("Cannot run MadSpin module: %s" % switch['reweight'])
-                    
-                    
+                        
         answers = list(available_mode) + ['auto', 'done']
         alias = {}
         for id, key in enumerate(switch_order):
-            if switch[key] != void and switch[key] in allowed_switch_value[key]:
+            if switch[key] != void and switch[key] in allowed_switch_value[key] and \
+                len(allowed_switch_value[key]) >1:
                 answers += ['%s=%s' % (key, s) for s in allowed_switch_value[key]]
                 #allow lower case for on/off
                 alias.update(dict(('%s=%s' % (key, s.lower()), '%s=%s' % (key, s))
@@ -4149,7 +4416,7 @@
             elif answer in ['0', 'auto', 'done']:
                 return 
             elif answer in special_values:
-                logger.info('Enter mode value: Go to the related mode', '$MG:color:BLACK')
+                logger.info('Enter mode value: %s. Go to the related mode' % answer, '$MG:color:BLACK')
                 #assign_switch('reweight', 'OFF')
                 #assign_switch('madspin', 'OFF')
                 if answer == 'LO':
@@ -4180,7 +4447,6 @@
                     return
             return switch
 
-
         modify_switch(mode, self.last_mode, switch)
         if switch['madspin'] == 'OFF' and  os.path.exists(pjoin(self.me_dir,'Cards','madspin_card.dat')):
             assign_switch('madspin', 'ON')
@@ -4272,9 +4538,9 @@
                     if mode in ['LO','aMC@LO','noshowerLO']:
                         self.run_name += '_LO' 
             self.set_run_name(self.run_name, self.run_tag, 'parton')
-            if int(self.run_card['ickkw']) == 3 and mode in ['LO', 'aMC@LO', 'noshowerLO']:
+            if self.run_card['ickkw'] == 3 and mode in ['LO', 'aMC@LO', 'noshowerLO']:
                 raise self.InvalidCmd("""FxFx merging (ickkw=3) not allowed at LO""")
-            elif int(self.run_card['ickkw']) == 3 and mode in ['aMC@NLO', 'noshower']:
+            elif self.run_card['ickkw'] == 3 and mode in ['aMC@NLO', 'noshower']:
                 logger.warning("""You are running with FxFx merging enabled.  To be able to merge
     samples of various multiplicities without double counting, you
     have to remove some events after showering 'by hand'.  Please
@@ -4290,7 +4556,7 @@
                         error = '''Stop opertation'''
                         self.ask_run_configuration(mode, options)
     #                    raise aMCatNLOError(error)
-            elif int(self.run_card['ickkw']) == -1 and mode in ['aMC@NLO', 'noshower']:
+            elif self.run_card['ickkw'] == -1 and mode in ['aMC@NLO', 'noshower']:
                     # NNLL+NLO jet-veto only possible for LO event generation or fNLO runs.
                 raise self.InvalidCmd("""NNLL+NLO jet veto runs (ickkw=-1) only possible for fNLO or LO.""")
         if 'aMC@' in mode or mode == 'onlyshower':

=== modified file 'madgraph/interface/coloring_logging.py' (properties changed: -x to +x)
=== modified file 'madgraph/interface/common_run_interface.py' (properties changed: -x to +x)
--- madgraph/interface/common_run_interface.py	2015-08-27 19:32:34 +0000
+++ madgraph/interface/common_run_interface.py	2015-10-20 14:02:27 +0000
@@ -221,19 +221,17 @@
             #restrict_file = None
             #if os.path.exists(pjoin(ufo_path, 'restrict_default.dat')):
             #    restrict_file = pjoin(ufo_path, 'restrict_default.dat')
+            
+            force_CMS = self.mother and self.mother.options['complex_mass_scheme']
             model = import_ufo.import_model(modelname, decay=True, 
-                        restrict=True)
-            if self.mother and self.mother.options['complex_mass_scheme']:
-                model.change_mass_to_complex_scheme()
+                                   restrict=True, complex_mass_scheme=force_CMS)
         else:
-            model = import_ufo.import_model(pjoin(
-                           self.me_dir,'bin','internal', 'ufomodel'),decay=True)
             #pattern for checking complex mass scheme.
             has_cms = re.compile(r'''set\s+complex_mass_scheme\s*(True|T|1|true|$|;)''')
-            if has_cms.search(open(pjoin(self.me_dir,'Cards','proc_card_mg5.dat')\
-                                   ).read()):
-                model.change_mass_to_complex_scheme()
-   
+            force_CMS =  has_cms.search(open(pjoin(self.me_dir,'Cards',
+                                                   'proc_card_mg5.dat')).read())
+            model = import_ufo.import_model(pjoin(self.me_dir,'bin','internal',
+                         'ufomodel'), decay=True, complex_mass_scheme=force_CMS)
             
 #        if not hasattr(model.get('particles')[0], 'partial_widths'):
 #            raise self.InvalidCmd, 'The UFO model does not include partial widths information. Impossible to compute widths automatically'
@@ -527,6 +525,8 @@
 
         self.me_dir = me_dir
         self.options = options
+        
+        self.param_card_iterator = [] #an placeholder containing a generator of paramcard for scanning
 
         # usefull shortcut
         self.status = pjoin(self.me_dir, 'status')
@@ -658,7 +658,7 @@
 
             if amcatnlo and not keepwidth:
                 # force particle in final states to have zero width
-                pids = self.get_pid_final_states()
+                pids = self.get_pid_final_initial_states()
                 # check those which are charged under qcd
                 if not MADEVENT and pjoin(self.me_dir,'bin','internal') not in sys.path:
                         sys.path.insert(0,pjoin(self.me_dir,'bin','internal'))
@@ -1351,8 +1351,8 @@
         self.update_status('delphes done', level='delphes', makehtml=False)
 
     ############################################################################
-    def get_pid_final_states(self):
-        """Find the pid of all particles in the final states"""
+    def get_pid_final_initial_states(self):
+        """Find the pid of all particles in the final and initial states"""
         pids = set()
         subproc = [l.strip() for l in open(pjoin(self.me_dir,'SubProcesses',
                                                                  'subproc.mg'))]
@@ -1363,7 +1363,7 @@
             group = pat.findall(text)
             for particles in group:
                 particles = particles.split(',')
-                pids.update(set(particles[nb_init:]))
+                pids.update(set(particles))
 
         return pids
 
@@ -1594,12 +1594,29 @@
 
 
     def check_param_card(self, path, run=True):
-        """Check that all the width are define in the param_card.
-        If some width are set on 'Auto', call the computation tools."""
+        """
+        1) Check that no scan parameter are present
+        2) Check that all the width are define in the param_card.
+        - If a scan parameter is define. create the iterator and recall this fonction 
+          on the first element.
+        - If some width are set on 'Auto', call the computation tools."""
         
-        pattern = re.compile(r'''decay\s+(\+?\-?\d+)\s+auto(@NLO|)''',re.I)
+        pattern_scan = re.compile(r'''^(decay)?[\s\d]*scan''', re.I+re.M)  
+        pattern_width = re.compile(r'''decay\s+(\+?\-?\d+)\s+auto(@NLO|)''',re.I)
         text = open(path).read()
-        pdg_info = pattern.findall(text)
+               
+        if pattern_scan.search(text):
+            if not isinstance(self, cmd.CmdShell):
+                # we are in web mode => forbid scan due to security risk
+                raise Exception, "Scan are not allowed in web mode"
+            # at least one scan parameter found. create an iterator to go trough the cards
+            main_card = check_param_card.ParamCardIterator(text)
+            self.param_card_iterator = main_card
+            first_card = main_card.next(autostart=True)
+            first_card.write(path)
+            return self.check_param_card(path, run)
+        
+        pdg_info = pattern_width.findall(text)
         if pdg_info:
             if run:
                 logger.info('Computing the width set on auto in the param_card.dat')
@@ -2320,40 +2337,8 @@
         # Read the comment of the param_card_default to find name variable for
         # the param_card also check which value seems to be constrained in the
         # model.
-        for bname, block in default_param.items():
-            for lha_id, param in block.param_dict.items():
-                all_var = []
-                comment = param.comment
-                # treat merge parameter
-                if comment.strip().startswith('set of param :'):
-                    all_var = list(re.findall(r'''[^-]1\*(\w*)\b''', comment))
-                # just the variable name as comment
-                elif len(comment.split()) == 1:
-                    all_var = [comment.strip().lower()]
-                # either contraction or not formatted
-                else:
-                    split = comment.split()
-                    if len(split) >2 and split[1] == ':':
-                        # NO VAR associated
-                        self.restricted_value[(bname, lha_id)] = ' '.join(split[1:])
-                    elif len(split) == 2:
-                        if re.search(r'''\[[A-Z]\]eV\^''', split[1]):
-                            all_var = [comment.strip().lower()]
-                    elif len(split) >=2 and split[1].startswith('('):
-                        all_var = [split[0].strip().lower()]
-                    else:
-                        if not bname.startswith('qnumbers'):
-                            logger.debug("not recognize information for %s %s : %s",
-                                      bname, lha_id, comment)
-                        # not recognized format
-                        continue
-
-                for var in all_var:
-                    var = var.lower()
-                    if var in self.pname2block:
-                        self.pname2block[var].append((bname, lha_id))
-                    else:
-                        self.pname2block[var] = [(bname, lha_id)]
+        self.pname2block, self.restricted_value = \
+                                              default_param.analyze_param_card()
 
         if run_card_def:
             self.run_set = run_card_def.keys() + self.run_card.hidden_param
@@ -2622,12 +2607,15 @@
         """ edit the value of one parameter in the card"""
         
         args = self.split_arg(line)
+        # fix some formatting problem
         if '=' in args[-1]:
             arg1, arg2 = args.pop(-1).split('=')
             args += [arg1, arg2]
         if '=' in args:
             args.remove('=')
-        args[:-1] = [ a.lower() for a in args[:-1]]
+        # do not set lowercase the case-sensitive parameters from the shower_card
+        if args[0].lower() not in ['analyse', 'extralibs', 'extrapaths', 'includepaths']:
+            args[:-1] = [ a.lower() for a in args[:-1]]
         # special shortcut:
         if args[0] in self.special_shortcut:
             if len(args) == 1:
@@ -2793,6 +2781,11 @@
         ### PARAM_CARD WITH BLOCK NAME -----------------------------------------
         elif (args[start] in self.param_card or args[start] == 'width') \
                                                   and card in ['','param_card']:
+            #special treatment for scan
+            if any(t.startswith('scan') for t in args):
+                index = [i for i,t in enumerate(args) if t.startswith('scan')][0]
+                args = args[:index] + [' '.join(args[index:])]
+                
             if args[start] in self.conflict and card == '':
                 text  = 'ambiguous name (present in more than one card). Please specify which card to edit'
                 text += ' in the format < set card parameter value>'
@@ -2836,7 +2829,7 @@
                     text += "You need to match this expression for external program (such pythia)."
                     logger.warning(text)
 
-                if args[-1].lower() in ['default', 'auto', 'auto@nlo']:
+                if args[-1].lower() in ['default', 'auto', 'auto@nlo'] or args[-1].startswith('scan'):
                     self.setP(args[start], key, args[-1])
                 else:
                     try:
@@ -3027,6 +3020,18 @@
                 if block != 'decay':
                     logger.warning('Invalid input: \'Auto\' value only valid for DECAY')
                     return
+            elif value.startswith('scan'):
+                if ':' not in value:
+                    logger.warning('Invalid input: \'scan\' mode requires a \':\' before the definition.')
+                    return
+                tag = value.split(':')[0]
+                tag = tag[4:].strip()
+                if tag and not tag.isdigit():
+                    logger.warning('Invalid input: scan tag need to be integer and not "%s"' % tag)
+                    return
+                
+                
+                pass
             else:
                 try:
                     value = float(value)
@@ -3111,7 +3116,6 @@
             logger.warning("invalid command for decay. Line ignored")
             return
         
-        misc.sprint( line, "-add" in line)
         if "-add" in line:
             # just to have to add the line to the end of the file
             particle = line.split('>')[0].strip()

=== modified file 'madgraph/interface/extended_cmd.py' (properties changed: -x to +x)
--- madgraph/interface/extended_cmd.py	2015-08-13 17:08:17 +0000
+++ madgraph/interface/extended_cmd.py	2015-10-20 14:02:27 +0000
@@ -80,7 +80,7 @@
         if readline and not 'libedit' in readline.__doc__:
             readline.set_completion_display_matches_hook(self.print_suggestions)
 
-    def deal_multiple_categories(self, dico):
+    def deal_multiple_categories(self, dico, forceCategory=False):
         """convert the multiple category in a formatted list understand by our
         specific readline parser"""
 
@@ -92,7 +92,7 @@
             return out
 
         # check if more than one categories but only one value:
-        if all(len(s) <= 1 for s in dico.values() ):
+        if not forceCategory and all(len(s) <= 1 for s in dico.values() ):
             values = set((s[0] for s in dico.values() if len(s)==1))
             if len(values) == 1:
                 return values
@@ -115,9 +115,9 @@
             opt.sort()
             out += opt
 
-            
-        if valid == 1:
+        if not forceCategory and valid == 1:
             out = out[1:]
+            
         return out
     
     @debug()

=== modified file 'madgraph/interface/launch_ext_program.py' (properties changed: -x to +x)
--- madgraph/interface/launch_ext_program.py	2015-08-18 09:16:14 +0000
+++ madgraph/interface/launch_ext_program.py	2015-10-20 14:02:27 +0000
@@ -481,7 +481,7 @@
         assert hasattr(self, 'cluster')
         assert hasattr(self, 'multicore')
         assert hasattr(self, 'name')
-#        assert hasattr(self, 'shell')
+        assert hasattr(self, 'shell')
 
         self.unit = unit
         self.run_mode = run_mode
@@ -519,7 +519,8 @@
                 nb_node=max_node
                 
         import madgraph.interface.amcatnlo_run_interface as run_int
-        if hasattr(self, 'shell'):
+        
+        if hasattr(self, 'shell') and self.shell:
             usecmd = run_int.aMCatNLOCmdShell(me_dir=self.running_dir, options = self.cmd_int.options)
         else:
             usecmd = run_int.aMCatNLOCmd(me_dir=self.running_dir, options = self.cmd_int.options)
@@ -542,7 +543,7 @@
                      usecmd, interface=False)
         #launch.me_dir = self.running_dir
         option_line = ' '.join([' --%s' % opt for opt in self.options.keys() \
-                if self.options[opt] and not opt in ['cluster', 'multicore', 'name', 'appl_start_grid']])
+                if self.options[opt] and not opt in ['cluster', 'multicore', 'name', 'appl_start_grid','shell']])
         if self.options['name']:
             option_line += ' --name %s' %  self.options['name']
         if 'appl_start_grid' in self.options and  self.options['appl_start_grid']:

=== modified file 'madgraph/interface/loop_interface.py' (properties changed: -x to +x)
--- madgraph/interface/loop_interface.py	2015-08-14 16:02:08 +0000
+++ madgraph/interface/loop_interface.py	2015-10-20 14:02:27 +0000
@@ -223,17 +223,20 @@
         if self._curr_amps and self._curr_amps[0].get_ninitial() != \
             proc.get_ninitial():
             raise self.InvalidCmd("Can not mix processes with different number of initial states.")               
-            
+
+#        It is partially supported for now if the initial state is not charged
+#        under the gauge group perturbed.
 #        if proc.get_ninitial()==1 and tool=='aMC@NLO':            
 #            raise self.InvalidCmd("At this stage %s cannot handle decay process."%tool+\
 #                                  "\nIt is however a straight-forward extension which "+\
 #                                  "will come out with the next release.")                           
 
-        if isinstance(proc, base_objects.ProcessDefinition) and mode.startswith('ML5'):
-            if proc.has_multiparticle_label():        
-                raise self.InvalidCmd, \
+#        Now all checks should support multi-particle label for loops as well.
+        if isinstance(proc, base_objects.ProcessDefinition) and mode=='ML5':
+            if proc.has_multiparticle_label():
+                raise self.InvalidCmd(
                   "When running ML5 standalone, multiparticle labels cannot be"+\
-                  " employed."
+                  " employed.")
         
         if proc['decay_chains']:
             raise self.InvalidCmd(
@@ -303,7 +306,7 @@
                 model_path = self._curr_model.get('modelpath')
                 model_name = self._curr_model.get('name')
                 if model_name.split('-')[0]=='loop_sm':
-		           model_name = model_name[5:]
+                    model_name = model_name[5:]