Merge lp:~maddevelopers/mg5amcnlo/3.0.3 into lp:mg5amcnlo

[Olivier Mattelaer]

Hi,

I think that we should release this quite fast.
One important point is that this include the compatibility with gcc10.
So I would not wait for the 4FS branch to be merged (but if it is release before the end of next week). We can still release the 4FS independently if needed (in 3.0.4)

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/3.0.3]

966. By marco zaro on 2020-06-30

fix for some unit tests (more needed)

967. By marco zaro on 2020-06-30

trivial updates for iotests

968. By marco zaro on 2020-06-30

fixed acceptance tests

969. By marco zaro on 2020-06-30

reverted to rev 363 (other commits were for another branch)

970. By olivier-mattelaer on 2020-07-04

fixing HWU issue

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2020-03-11 09:28:14 +0000
+++ MadSpin/decay.py	2020-07-04 19:03:17 +0000
@@ -544,7 +544,7 @@
                     # record the mass for the reshuffling phase, 
                     # in case the point passes the reweighting creteria
                     tree[tag]["mass"] = m
-                    #update the weigth of the phase-space point
+                    #update the weight of the phase-space point
                     weight=weight*jac
                 # for checking conservation of energy
                 mass_sum -= m
@@ -2749,7 +2749,7 @@
             for proc in processes:
                 if not proc.strip().startswith(('add','generate')):
                     proc = 'add process %s' % proc
-                commandline += self.get_proc_with_decay(proc, decay_text, mgcmd._curr_model)
+                commandline += self.get_proc_with_decay(proc, decay_text, mgcmd._curr_model, self.options)
                 
             commandline = commandline.replace('add process', 'generate',1)
             logger.info(commandline)
@@ -2835,7 +2835,7 @@
 #                    assert decay.shell_string() in self.all_decay
             
     @staticmethod
-    def get_proc_with_decay(proc, decay_text, model):
+    def get_proc_with_decay(proc, decay_text, model, msoptions=None):
         
         commands = []
         if '[' in proc:
@@ -2880,6 +2880,13 @@
             else:
                 baseproc = new_proc
                 proc_nb = ''      
+
+            if msoptions and msoptions['global_order_coupling']:
+                if '@' in proc_nb:
+                    proc_nb += " %s" % msoptions['global_order_coupling']
+                else:
+                    proc_nb += " @0 %s" %  msoptions['global_order_coupling']                
+
             nb_comma = baseproc.count(',')
             if nb_comma == 0:
                 commands.append("%s, %s %s %s" % (baseproc, decay_text, proc_nb, options))
@@ -3087,7 +3094,7 @@
         for production in self.all_ME.values():
             decay_set.add(production['decaying'])
 
-        numberev = self.options['Nevents_for_max_weigth'] # number of events
+        numberev = self.options['Nevents_for_max_weight'] # number of events
         numberps = self.options['max_weight_ps_point'] # number of phase pace points per event
         
         logger.info('  ')

=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2020-02-25 14:05:36 +0000
+++ MadSpin/interface_madspin.py	2020-07-04 19:03:17 +0000
@@ -60,14 +60,14 @@
 
         self.add_param("max_weight", -1)
         self.add_param('curr_dir', os.path.realpath(os.getcwd()))
-        self.add_param('Nevents_for_max_weigth', 0)
+        self.add_param('Nevents_for_max_weight', 0)
         self.add_param("max_weight_ps_point", 400)
         self.add_param('BW_cut', -1)
         self.add_param('nb_sigma', 0.)
         self.add_param('ms_dir', '')
         self.add_param('max_running_process', 100)
         self.add_param('onlyhelicity', False)
-        self.add_param('spinmode', "madspin", allowed=['madspin','none','onshell'])
+        self.add_param('spinmode', "madspin", allowed=['full','madspin','none','onshell'])
         self.add_param('use_old_dir', False, comment='should be use only for faster debugging')
         self.add_param('run_card', '' , comment='define cut for spinmode==none. Path to run_card to use')
         self.add_param('fixed_order', False, comment='to activate fixed order handling of counter-event')
@@ -76,6 +76,7 @@
         self.add_param('new_wgt', 'cross-section' ,allowed=['cross-section', 'BR'], comment="if not consistent number of particles, choose what to do for the weight. (BR: means local according to number of part, cross use the force cross-section")
         self.add_param('input_format', 'auto', allowed=['auto','lhe', 'hepmc', 'lhe_no_banner'])
         self.add_param('frame_id', 6)
+        self.add_param('global_order_coupling', '')
         
     ############################################################################
     ##  Special post-processing of the options                                ## 
@@ -137,8 +138,8 @@
         
         self.decay = madspin.decay_misc()
         self.model = None
-        self.mode = "madspin" # can be flat/bridge change the way the decay is done.
-                              # note amc@nlo does not support bridge.
+        #self.mode = "madspin" # can be flat/bridge change the way the decay is done.
+        #                      # note amc@nlo does not support bridge.
         
         self.options = MadSpinOptions()
         
@@ -227,10 +228,10 @@
         
         if 'mgruncard' in self.banner:
             run_card = self.banner.charge_card('run_card')
-            if not self.options['Nevents_for_max_weigth']:
+            if not self.options['Nevents_for_max_weight']:
                 nevents = run_card['nevents']
                 N_weight = max([75, int(3*nevents**(1/3))])
-                self.options['Nevents_for_max_weigth'] = N_weight
+                self.options['Nevents_for_max_weight'] = N_weight
                 N_sigma = max(4.5, math.log(nevents,7.7))
                 self.options['nb_sigma'] = N_sigma
             if self.options['BW_cut'] == -1:
@@ -242,8 +243,8 @@
             else:
                 self.options['frame_id'] = 6
         else:
-            if not self.options['Nevents_for_max_weigth']:
-                self.options['Nevents_for_max_weigth'] = 75
+            if not self.options['Nevents_for_max_weight']:
+                self.options['Nevents_for_max_weight'] = 75
                 self.options['nb_sigma'] = 4.5
             if self.options['BW_cut'] == -1:
                 self.options['BW_cut'] = 15.0
@@ -314,7 +315,7 @@
                         key = line.split()[1]
                         if key in self.multiparticles_ms:
                             del self.multiparticles_ms[key]            
-            elif line.startswith('set'):
+            elif line.startswith('set') and not line.startswith('set gauge'):
                 self.mg5cmd.exec_cmd(line, printcmd=False, precmd=False, postcmd=False)
             elif line.startswith('import model'):
                 if model_name in line:
@@ -458,7 +459,7 @@
                 raise self.InvalidCmd('second argument should be a path to a existing directory')
         
         elif args[0] == "spinmode":
-            if args[1].lower() not in ["full", "onshell", "none"]:
+            if args[1].lower() not in ["full", "onshell", "none", "madspin"]:
                 raise self.InvalidCmd("spinmode can only take one of those 3 value: full/onshell/none")
              
         elif args[0] == "run_card":
@@ -471,6 +472,8 @@
                 arg, value = data.split("=")
                 args.append(arg)
                 args.append(value)
+        elif args[0] == 'Nevents_for_max_weigth':
+            args[0] = 'Nevents_for_max_weight'
         
     def do_set(self, line):
         """ add one of the options """
@@ -1204,7 +1207,6 @@
                         run_card = self.run_card
                     else:
                         run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))                        
-                    
                     run_card["iseed"] = self.options['seed']
                     run_card['gridpack'] = True
                     run_card['systematics_program'] = 'False'
@@ -1213,6 +1215,7 @@
                     param_card = self.banner['slha']
                     open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
                     self.options['seed'] += 1
+                    self.seed = self.options['seed'] 
                     # actually creation
                     me5_cmd.exec_cmd("generate_events run_01 -f")
                     if output_width:
@@ -1222,9 +1225,11 @@
                             width *= me5_cmd.results.current['cross']
                     me5_cmd.exec_cmd("exit")                        
                     #remove pointless informat
-                    misc.call(["rm", "Cards", "bin", 'Source', 'SubProcesses'], cwd=decay_dir)
-                    misc.call(['tar', '-xzpvf', 'run_01_gridpack.tar.gz'], cwd=decay_dir)
-            
+                    if not os.path.exists(pjoin(decay_dir, 'run.sh')):
+                        devnull = open('/dev/null','w')
+                        misc.call(["rm", "Cards", "bin", 'Source', 'SubProcesses'], cwd=decay_dir,stdout=devnull, stderr=-2)
+                        misc.call(['tar', '-xzpvf', 'run_01_gridpack.tar.gz'], cwd=decay_dir,stdout=devnull, stderr=-2)
+                        devnull.close()
             # Now generate the events
             if not self.options['ms_dir']:
                 if decay_dir in self.me_int:
@@ -1364,7 +1369,7 @@
 
 
         # 2. Generate the events requested
-        nevents_for_max = self.options['Nevents_for_max_weigth']
+        nevents_for_max = self.options['Nevents_for_max_weight']
         if nevents_for_max == 0 :
             nevents_for_max = 75
         nevents_for_max *= self.options['max_weight_ps_point']
@@ -1460,8 +1465,8 @@
                 production, counterevt= production[0], production[1:]
             if curr_event and self.efficiency and curr_event % 10 == 0 and float(str(curr_event)[1:]) ==0:
                 logger.info("decaying event number %s. Efficiency: %s [%s s]" % (curr_event, 1/self.efficiency, time.time()-start))
-            else:
-                logger.info("next event [%s]", time.time()-start)
+            #else:
+            #    logger.info("next event [%s]", time.time()-start)
             while 1:
                 nb_try +=1
                 decays = self.get_decay_from_file(production, evt_decayfile, nb_event-curr_event)
@@ -1559,7 +1564,7 @@
         if self.options['ms_dir'] and os.path.exists(pjoin(self.options['ms_dir'], 'max_wgt')):
             return float(open(pjoin(self.options['ms_dir'], 'max_wgt'),'r').read())
         
-        nevents = self.options['Nevents_for_max_weigth']
+        nevents = self.options['Nevents_for_max_weight']
         if nevents == 0 :
             nevents = 75
         
@@ -1700,7 +1705,7 @@
                     decay_text.append('(%s)' % decay)
                 else:
                     decay_text.append(decay)
-                processes_decay.append(decay)
+                processes_decay.append(decay)            
         decay_text = ', '.join(decay_text)
         processes += []
         
@@ -1715,12 +1720,15 @@
                 except ValueError:
                     raise MadSpinError, 'MadSpin didn\'t allow order restriction after the @ comment: \"%s\" not valid' % proc_nb
                 proc_nb = '@ %i' % proc_nb 
+                if self.options['global_order_coupling']:
+                    proc_nb = '%s %s' % (proc_nb, self.options['global_order_coupling'])
             else:
-                proc_nb = ''     
+                if self.options['global_order_coupling']:
+                    proc_nb = '@0 %s ' % self.options['global_order_coupling']    
                 
             rwgt_interface.ReweightInterface.get_LO_definition_from_NLO()        
         
-
+        raise Exception
 
 if __name__ == '__main__':
     

=== modified file 'Template/Common/Cards/madspin_card_default.dat'
--- Template/Common/Cards/madspin_card_default.dat	2017-01-13 09:56:03 +0000
+++ Template/Common/Cards/madspin_card_default.dat	2020-07-04 19:03:17 +0000
@@ -14,7 +14,7 @@
 #Some options (uncomment to apply)
 #
 # set seed 1
-# set Nevents_for_max_weigth 75 # number of events for the estimate of the max. weight
+# set Nevents_for_max_weight 75 # number of events for the estimate of the max. weight
 # set BW_cut 15                 # cut on how far the particle can be off-shell
 # set spinmode onshell          # Use one of the madspin special mode
  set max_weight_ps_point 400  # number of PS to estimate the maximum for each event

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2020-01-10 20:02:44 +0000
+++ Template/LO/Cards/run_card.dat	2020-07-04 19:03:17 +0000
@@ -165,7 +165,8 @@
  %(mmnl)s   = mmnl    ! min invariant mass for all letpons (l+- and vl) 
  %(mmnlmax)s  = mmnlmax ! max invariant mass for all letpons (l+- and vl) 
 #IF(LL)# #*********************************************************************
-#IF(LL)# # Minimum and maximum pt for 4-momenta sum of leptons                *
+#IF(LL)# # Minimum and maximum pt for 4-momenta sum of leptons / neutrino     *
+#IF(LL)# #  for pair of lepton includes only same flavor, opposite charge
 #IF(LL)# #*********************************************************************
  %(ptllmin)s   = ptllmin  ! Minimum pt for 4-momenta sum of leptons(l and vl)
  %(ptllmax)s  = ptllmax  ! Maximum pt for 4-momenta sum of leptons(l and vl)

=== modified file 'Template/LO/Source/PDF/PhotonFlux.f'
--- Template/LO/Source/PDF/PhotonFlux.f	2019-12-04 21:34:25 +0000
+++ Template/LO/Source/PDF/PhotonFlux.f	2020-07-04 19:03:17 +0000
@@ -1,7 +1,8 @@
 c/* ********************************************************* */
 c/*  Equivalent photon approximation structure function.   * */
+c/*   V.M.Budnev et al., Phys.Rep. 15C (1975) 181          * */
 c/*     Improved Weizsaecker-Williams formula              * */
-c/*   V.M.Budnev et al., Phys.Rep. 15C (1975) 181          * */
+c/*     http://inspirehep.net/record/359425                * */  
 c/* ********************************************************* */
 c   provided by Tomasz Pierzchala - UCL
 

=== modified file 'Template/LO/Source/PDF/makefile'
--- Template/LO/Source/PDF/makefile	2018-01-26 23:06:46 +0000
+++ Template/LO/Source/PDF/makefile	2020-07-04 19:03:17 +0000
@@ -12,9 +12,14 @@
 
 ifdef lhapdf
   ifeq ($(lhapdfversion),5)
-    PDF         = pdfwrap_lhapdf.o pdf_lhapdf.o pdg2pdf_lhapdf.o opendata.o  PhotonFlux.o
+    $(error Bad lhadpfversion version 6 is now required)
   else
-    PDF         = pdfwrap_lhapdf.o pdf_lhapdf6.o pdg2pdf_lhapdf6.o opendata.o PhotonFlux.o
+    ifeq ($(lhapdfsubversion),1) # 6.1.X
+	PDF         = pdfwrap_lhapdf.o pdf_lhapdf6.o pdg2pdf_lhapdf6.o opendata.o PhotonFlux.o
+    else # 6.2.X
+        CXXFLAGS+=-std=c++11
+	PDF         = pdfwrap_lhapdf.o pdf_lhapdf62.o pdg2pdf_lhapdf6.o opendata.o PhotonFlux.o
+    endif    
   endif
 else
   PDF         = Ctq6Pdf.o pdfwrap.o opendata.o pdf.o PhotonFlux.o pdg2pdf.o NNPDFDriver.o

=== added file 'Template/LO/Source/PDF/pdf_lhapdf62.cc'
--- Template/LO/Source/PDF/pdf_lhapdf62.cc	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/PDF/pdf_lhapdf62.cc	2020-07-04 19:03:17 +0000
@@ -0,0 +1,1569 @@
+// -*- C++ -*-
+//
+// This file is part of LHAPDF
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
+//
+#include "LHAPDF/PDF.h"
+#include "LHAPDF/PDFSet.h"
+#include "LHAPDF/PDFIndex.h"
+#include "LHAPDF/Factories.h"
+#include "LHAPDF/Utils.h"
+#include "LHAPDF/Paths.h"
+#include "LHAPDF/Version.h"
+#include "LHAPDF/LHAGlue.h"
+#include <cstring>
+
+using namespace std;
+
+
+// We have to create and initialise some common blocks here for backwards compatibility
+struct w50512 {
+  double qcdl4, qcdl5;
+};
+w50512 w50512_;
+
+struct w50513 {
+  double xmin, xmax, q2min, q2max;
+};
+w50513 w50513_;
+
+struct lhapdfr {
+  double qcdlha4, qcdlha5;
+  int nfllha;
+};
+lhapdfr lhapdfr_;
+
+
+
+namespace   { //< Unnamed namespace to restrict visibility to this file
+
+
+  /// @brief PDF object storage here is a smart pointer to ensure deletion of created PDFs
+  typedef std::shared_ptr<LHAPDF::PDF> PDFPtr;
+
+
+  /// @brief A struct for handling the active PDFs for the Fortran interface.
+  ///
+  /// We operate in a string-based way, since maybe there will be sets with names, but no
+  /// index entry in pdfsets.index.
+  ///
+  /// @todo Can we avoid the strings and just work via the LHAPDF ID and factory construction?
+  ///
+  /// Smart pointers are used in the native map used for PDF member storage so
+  /// that they auto-delete if the PDFSetHandler that holds them goes out of
+  /// scope (i.e. is overwritten).
+  struct PDFSetHandler {
+
+    /// Default constructor
+    ///
+    /// It'll be stored in a map so we need one of these...
+    PDFSetHandler() : currentmem(0)
+    { }
+
+    /// Constructor from a PDF set name
+    ///
+    /// @note If the set name contains a member specification, i.e. myname/2,
+    /// that member rather than the central one will be initialised and made
+    /// current.
+    PDFSetHandler(const string& name) {
+      pair<string, int> set_mem = LHAPDF::lookupPDF(name);
+      // First check that the lookup was successful, i.e. it was a valid ID for the LHAPDF6 set collection
+      if (set_mem.first.empty() || set_mem.second < 0)
+        throw LHAPDF::UserError("Could not find a valid PDF with string = " + name);
+      // Try to load this PDF
+      setname = set_mem.first;
+      loadMember(set_mem.second);
+    }
+
+    /// Constructor from a PDF set's LHAPDF ID code
+    ///
+    /// @note The set member given by the ID (rather than the central one) will
+    /// be initialised and made current.
+    PDFSetHandler(int lhaid) {
+      pair<string,int> set_mem = LHAPDF::lookupPDF(lhaid);
+      // First check that the lookup was successful, i.e. it was a valid ID for the LHAPDF6 set collection
+      if (set_mem.first.empty() || set_mem.second < 0)
+        throw LHAPDF::UserError("Could not find a valid PDF with LHAPDF ID = " + LHAPDF::to_str(lhaid));
+      // Try to load this PDF
+      setname = set_mem.first;
+      loadMember(set_mem.second);
+    }
+
+    /// @brief Load a new PDF member, set it to be active
+    ///
+    /// If it's already loaded, the existing object will not be reloaded.
+    void loadMember(int mem) {
+      if (mem < 0)
+        throw LHAPDF::UserError("Tried to load a negative PDF member ID: " + LHAPDF::to_str(mem) + " in set " + setname);
+      if (members.find(mem) == members.end())
+        members[mem] = PDFPtr(LHAPDF::mkPDF(setname, mem));
+      currentmem = mem;
+      //return members[mem];
+    }
+
+    /// Actively delete a PDF member to save memory, set the active member to be the next available, or 0
+    void unloadMember(int mem) {
+      members.erase(mem);
+      const int nextmem = (!members.empty()) ? members.begin()->first : 0;
+      loadMember(nextmem);
+    }
+
+    /// @brief Get a PDF member, making it active
+    ///
+    /// Non-const because it can secretly load the member. Not that constness
+    /// matters in a Fortran interface utility function!
+    const PDFPtr member(int mem) {
+      loadMember(mem);
+      return members.find(mem)->second;
+    }
+
+    /// Get the currently active PDF member
+    ///
+    /// Non-const because it can secretly load the member. Not that constness
+    /// matters in a Fortran interface utility function!
+    const PDFPtr activeMember() {
+      return member(currentmem);
+    }
+
+    /// Get the currently active PDF member
+    ///
+    /// Non-const because it can secretly load the member. Not that constness
+    /// matters in a Fortran interface utility function!
+    void setActiveMember(int mem) {
+      loadMember(mem);
+    }
+
+    /// The currently active member in this set
+    int currentmem;
+
+    /// Name of this set
+    string setname;
+
+    /// Map of pointers to selected member PDFs
+    ///
+    // /// It's mutable so that a "const" member-getting operation can implicitly
+    // /// load a new PDF object. Good idea / bad idea? Disabled for now.
+    // mutable map<int, PDFPtr> members;
+    map<int, PDFPtr> members;
+  };
+
+
+  /// Collection of active sets
+  static map<int, PDFSetHandler> ACTIVESETS;
+
+  /// The currently active set
+  int CURRENTSET = 0;
+
+}
+
+
+
+string lhaglue_get_current_pdf(int nset) {
+  if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    return "NONE";
+  CURRENTSET = nset;
+  return ACTIVESETS[nset].activeMember()->set().name() + " (" +
+    LHAPDF::to_str(ACTIVESETS[nset].activeMember()->lhapdfID()) + ")";
+}
+
+
+
+namespace {
+
+
+  /// C-string -> Fortran-string converter
+  ///
+  /// Credit: https://stackoverflow.com/questions/10163485/passing-char-arrays-from-c-to-fortran
+  void cstr_to_fstr(const char* cstring, char* fstring, std::size_t fstring_len) {
+    std::size_t inlen = std::strlen(cstring);
+    std::size_t cpylen = std::min(inlen, fstring_len);
+    // TODO: truncation error or warning
+    //if (inlen > fstring_len) FOOOOO();
+    std::copy(cstring, cstring+cpylen, fstring);
+    std::fill(fstring+cpylen, fstring+fstring_len, ' ');
+  }
+
+
+  /// C++-string -> Fortran-string converter
+  void ccstr_to_fstr(const string& ccstring, char* fstring, std::size_t fstring_len) {
+    const char* cstring = ccstring.c_str();
+    cstr_to_fstr(cstring, fstring, fstring_len);
+  }
+
+
+  /// Fortran-string -> C++-string converter
+  string fstr_to_ccstr(const char* fstring, const std::size_t fstring_len, bool spcpad=false) {
+    // Allocate space for an equivalent C-string (with an extra terminating null byte)
+    char* s = new char[fstring_len+1];
+    // Copy all characters and add the terminating null byte
+    strncpy(s, fstring, fstring_len);
+    s[fstring_len] = '\0';
+    // Replace all trailing spaces with null bytes unless explicitly stopped
+    if (!spcpad) {
+      for (int i = fstring_len-1; i >= 0; --i) {
+        if (s[i] != ' ') break;
+        s[i] = '\0';
+      }
+    }
+    string rtn(s); //< copy the result to a C++ string
+    delete[] s; //< clean up the dynamic array
+    return rtn;
+  }
+
+
+}
+
+
+extern "C" {
+
+
+  // NEW FORTRAN INTERFACE FUNCTIONS
+
+  /// Get the LHAPDF library version as a string
+  void lhapdf_getversion_(char* s, size_t len) {
+    cstr_to_fstr(LHAPDF_VERSION, s, len);
+  }
+
+
+  /// List of available PDF sets, returned as a space-separated string
+  void lhapdf_getpdfsetlist_(char* s, size_t len) {
+    string liststr;
+    for (const string& setname : LHAPDF::availablePDFSets()) {
+      if (!liststr.empty()) liststr += " ";
+      liststr += setname;
+    }
+    ccstr_to_fstr(liststr, s, len);
+  }
+
+
+  /// Get PDF data path (colon-separated if there is more than one element)
+  void lhapdf_getdatapath_(char* s, size_t len) {
+    string pathstr;
+    for (const string& path : LHAPDF::paths()) {
+      if (!pathstr.empty()) pathstr += ":";
+      pathstr += path;
+    }
+    ccstr_to_fstr(pathstr, s, len);
+  }
+
+  /// Set PDF data path(s)
+  void lhapdf_setdatapath_(const char* s, size_t len) {
+    LHAPDF::setPaths(fstr_to_ccstr(s, len));
+  }
+
+  /// Prepend to PDF data path
+  void lhapdf_prependdatapath_(const char* s, size_t len) {
+    LHAPDF::pathsPrepend(fstr_to_ccstr(s, len));
+  }
+
+  /// Append to PDF data path
+  void lhapdf_appenddatapath_(const char* s, size_t len) {
+    LHAPDF::pathsAppend(fstr_to_ccstr(s, len));
+  }
+
+
+  //------------------
+
+
+  void lhapdf_initpdfset_byname_(const int& nset, const char* name, int namelength) {
+    const string cname = fstr_to_ccstr(name, namelength);
+    ACTIVESETS[nset] = PDFSetHandler(cname);
+    CURRENTSET = nset;
+  }
+
+  void lhapdf_initpdfset_byid_(const int& nset, const int& lhaid) {
+    ACTIVESETS[nset] = PDFSetHandler(lhaid);
+    CURRENTSET = nset;
+  }
+
+  void lhapdf_delpdfset_(const int& nset) {
+    ACTIVESETS.erase(nset);
+    CURRENTSET = 0;
+  }
+
+  void lhapdf_delpdf_(const int& nset, const int& nmem) {
+    CURRENTSET = nset;
+    ACTIVESETS[CURRENTSET].unloadMember(nmem);
+  }
+
+
+  //------------------
+
+
+  void lhapdf_hasflavor(const int& nset, const int& nmem, const int& pid, int& rtn) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use set slot " + LHAPDF::to_str(nset) + " but it is not initialised");
+    rtn = ACTIVESETS[nset].member(nmem)->hasFlavor(pid) ? 1 : 0;
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+
+
+  void lhapdf_xfxq2_(const int& nset, const int& nmem, const int& pid, const double& x, const double& q2, double& xf) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use set slot " + LHAPDF::to_str(nset) + " but it is not initialised");
+    try {
+      xf = ACTIVESETS[nset].member(nmem)->xfxQ2(pid, x, q2);
+    } catch (const exception& e) {
+      xf = 0;
+    }
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+
+  void lhapdf_xfxq_(const int& nset, const int& nmem, const int& pid, const double& x, const double& q, double& xf) {
+    const double q2 = q*q;
+    lhapdf_xfxq2_(nset, nmem, pid, x, q2, xf);
+  }
+
+
+  void lhapdf_xfxq2_stdpartons_(const int& nset, const int& nmem, const int& pid, const double& x, const double& q2, double* xfs) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    // Evaluate for the 13 LHAPDF5 standard partons (-6..6)
+    for (size_t i = 0; i < 13; ++i) {
+      try {
+        xfs[i] = ACTIVESETS[nset].member(nmem)->xfxQ2(i-6, x, q2);
+      } catch (const exception& e) {
+        xfs[i] = 0;
+      }
+    }
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+
+  void lhapdf_xfxq_stdpartons_(const int& nset, const int& nmem, const int& pid, const double& x, const double& q, double* xfs) {
+    const double q2 = q*q;
+    lhapdf_xfxq2_stdpartons_(nset, nmem, pid, x, q2, xfs);
+  }
+
+
+  //-----------------
+
+
+  /// Get the alpha_s order for the set
+  void lhapdf_getorderas_(const int& nset, const int& nmem, int& oas) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    oas = ACTIVESETS[nset].member(nmem)->info().get_entry_as<int>("AlphaS_OrderQCD");
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+
+  /// Get the alpha_s(Q2) value for set nset
+  void lhapdf_alphasq2_(const int& nset, const int& nmem, const double& q2, double& alphas) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    alphas = ACTIVESETS[nset].member(nmem)->alphasQ2(q2);
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+
+  /// Get the alpha_s(Q) value for set nset
+  /// @todo Return value rather than return arg? Can we do that elsewhere, too, e.g. single-value PDF xf functions?
+  void lhapdf_alphasq_(const int& nset, const int& nmem, const double& q, double& alphas) {
+    const double q2 = q*q;
+    lhapdf_alphasq2_(nset, nmem, q2, alphas);
+  }
+
+
+  // Metadata functions
+
+  // /// Get the number of error members in the set (with special treatment for single member sets)
+  // void numberpdfm_(const int& nset, int& numpdf) {
+  //   if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //     throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  //   // Set equal to the number of members  for the requested set
+  //   numpdf=  ACTIVESETS[nset].activeMember()->info().get_entry_as<int>("NumMembers");
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+
+  // /// Get the max number of active flavours
+  // void getnfm_(const int& nset, int& nf) {
+  //   //nf = ACTIVESETS[nset].activeMember()->info().get_entry_as<int>("AlphaS_NumFlavors");
+  //   nf = ACTIVESETS[nset].activeMember()->info().get_entry_as<int>("NumFlavors");
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+
+  // /// Get nf'th quark mass
+  // void getqmassm_(const int& nset, const int& nf, double& mass) {
+  //   if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //     throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  //   if      (nf*nf ==  1) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MDown");
+  //   else if (nf*nf ==  4) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MUp");
+  //   else if (nf*nf ==  9) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MStrange");
+  //   else if (nf*nf == 16) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MCharm");
+  //   else if (nf*nf == 25) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MBottom");
+  //   else if (nf*nf == 36) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MTop");
+  //   else throw LHAPDF::UserError("Trying to get quark mass for invalid quark ID #" + LHAPDF::to_str(nf));
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+
+  // /// Get the nf'th quark threshold
+  // void getthresholdm_(const int& nset, const int& nf, double& Q) {
+  //   try {
+  //     if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //       throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  //     if      (nf*nf ==  1) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdDown");
+  //     else if (nf*nf ==  4) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdUp");
+  //     else if (nf*nf ==  9) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdStrange");
+  //     else if (nf*nf == 16) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdCharm");
+  //     else if (nf*nf == 25) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdBottom");
+  //     else if (nf*nf == 36) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdTop");
+  //     //else throw LHAPDF::UserError("Trying to get quark threshold for invalid quark ID #" + LHAPDF::to_str(nf));
+  //   } catch (...) {
+  //     getqmassm_(nset, nf, Q);
+  //   }
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+
+  // void getxminm_(const int& nset, const int& nmem, double& xmin) {
+  //   if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //     throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  //   const int activemem = ACTIVESETS[nset].currentmem;
+  //   ACTIVESETS[nset].loadMember(nmem);
+  //   xmin = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("XMin");
+  //   ACTIVESETS[nset].loadMember(activemem);
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+
+  // void getxmaxm_(const int& nset, const int& nmem, double& xmax) {
+  //   if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //     throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  //   const int activemem = ACTIVESETS[nset].currentmem;
+  //   ACTIVESETS[nset].loadMember(nmem);
+  //   xmax = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("XMax");
+  //   ACTIVESETS[nset].loadMember(activemem);
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+
+  // void getq2minm_(const int& nset, const int& nmem, double& q2min) {
+  //   if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //     throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  //   const int activemem = ACTIVESETS[nset].currentmem;
+  //   ACTIVESETS[nset].loadMember(nmem);
+  //   q2min = LHAPDF::sqr(ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("QMin"));
+  //   ACTIVESETS[nset].loadMember(activemem);
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+
+  // void getq2maxm_(const int& nset, const int& nmem, double& q2max) {
+  //   if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //     throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  //   const int activemem = ACTIVESETS[nset].currentmem;
+  //   ACTIVESETS[nset].loadMember(nmem);
+  //   q2max = LHAPDF::sqr(ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("QMax"));
+  //   ACTIVESETS[nset].loadMember(activemem);
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+
+  // void getminmaxm_(const int& nset, const int& nmem, double& xmin, double& xmax, double& q2min, double& q2max) {
+  //   if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //     throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  //   const int activemem = ACTIVESETS[nset].currentmem;
+  //   ACTIVESETS[nset].loadMember(nmem);
+  //   xmin = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("XMin");
+  //   xmax = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("XMax");
+  //   q2min = LHAPDF::sqr(ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("QMin"));
+  //   q2max = LHAPDF::sqr(ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("QMax"));
+  //   ACTIVESETS[nset].loadMember(activemem);
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+
+
+  // /// Backwards compatibility functions for LHAPDF5 calculations of
+  // /// PDF uncertainties and PDF correlations (G. Watt, March 2014).
+
+  // // subroutine GetPDFUncTypeM(nset,lMonteCarlo,lSymmetric)
+  // void getpdfunctypem_(const int& nset, int& lmontecarlo, int& lsymmetric) {
+  //   if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //     throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  //   const string errorType = ACTIVESETS[nset].activeMember()->set().errorType();
+  //   if (errorType == "replicas") { // Monte Carlo PDF sets
+  //     lmontecarlo = 1;
+  //     lsymmetric = 1;
+  //   } else if (errorType == "symmhessian") { // symmetric eigenvector PDF sets
+  //     lmontecarlo = 0;
+  //     lsymmetric = 1;
+  //   } else { // default: assume asymmetric Hessian eigenvector PDF sets
+  //     lmontecarlo = 0;
+  //     lsymmetric = 0;
+  //   }
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+  // // subroutine GetPDFUncType(lMonteCarlo,lSymmetric)
+  // void getpdfunctype_(int& lmontecarlo, int& lsymmetric) {
+  //   int nset1 = 1;
+  //   getpdfunctypem_(nset1, lmontecarlo, lsymmetric);
+  // }
+
+
+  // // subroutine GetPDFuncertaintyM(nset,values,central,errplus,errminus,errsym)
+  // void getpdfuncertaintym_(const int& nset, const double* values, double& central, double& errplus, double& errminus, double& errsymm) {
+  //   if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //     throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  //   const size_t nmem = ACTIVESETS[nset].activeMember()->set().size()-1;
+  //   const vector<double> vecvalues(values, values + nmem + 1);
+  //   LHAPDF::PDFUncertainty err = ACTIVESETS[nset].activeMember()->set().uncertainty(vecvalues, -1);
+  //   central = err.central;
+  //   errplus = err.errplus;
+  //   errminus = err.errminus;
+  //   errsymm = err.errsymm;
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+  // // subroutine GetPDFuncertainty(values,central,errplus,errminus,errsym)
+  // void getpdfuncertainty_(const double* values, double& central, double& errplus, double& errminus, double& errsymm) {
+  //   int nset1 = 1;
+  //   getpdfuncertaintym_(nset1, values, central, errplus, errminus, errsymm);
+  // }
+
+
+  // // subroutine GetPDFcorrelationM(nset,valuesA,valuesB,correlation)
+  // void getpdfcorrelationm_(const int& nset, const double* valuesA, const double* valuesB, double& correlation) {
+  //   if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //     throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  //   const size_t nmem = ACTIVESETS[nset].activeMember()->set().size()-1;
+  //   const vector<double> vecvaluesA(valuesA, valuesA + nmem + 1);
+  //   const vector<double> vecvaluesB(valuesB, valuesB + nmem + 1);
+  //   correlation = ACTIVESETS[nset].activeMember()->set().correlation(vecvaluesA,vecvaluesB);
+  //   // Update current set focus
+  //   CURRENTSET = nset;
+  // }
+  // // subroutine GetPDFcorrelation(valuesA,valuesB,correlation)
+  // void getpdfcorrelation_(const double* valuesA, const double* valuesB, double& correlation) {
+  //   int nset1 = 1;
+  //   getpdfcorrelationm_(nset1, valuesA, valuesB, correlation);
+  // }
+
+
+
+
+
+  //////////////////
+
+  // LHAPDF5 / PDFLIB COMPATIBILITY INTERFACE FUNCTIONS
+
+
+  // System-level info
+
+  /// LHAPDF library version
+  void getlhapdfversion_(char* s, size_t len) {
+    // strncpy(s, LHAPDF_VERSION, len);
+    cstr_to_fstr(LHAPDF_VERSION, s, len);
+  }
+
+
+  /// Does nothing, only provided for backward compatibility
+  void lhaprint_(int& a) {  }
+
+
+  /// Set LHAPDF parameters
+  ///
+  /// @note Only the verbosity parameters have any effect: PDF behaviour is not
+  /// controlled globally in LHAPDF6.
+  void setlhaparm_(const char* par, int parlength) {
+    const string cpar = LHAPDF::to_upper(fstr_to_ccstr(par, parlength));
+    if (cpar == "NOSTAT" || cpar == "16") {
+      cerr << "WARNING: Fortran call to control LHAPDF statistics collection has no effect" << endl;
+    } else if (cpar == "LHAPDF" || cpar == "17") {
+      cerr << "WARNING: Fortran call to globally control alpha_s calculation has no effect" << endl;
+    } else if (cpar == "EXTRAPOLATE" || cpar == "18") {
+      cerr << "WARNING: Fortran call to globally control PDF extrapolation has no effect" << endl;
+    } else if (cpar == "SILENT" || cpar == "LOWKEY") {
+      LHAPDF::setVerbosity(0);
+    } else if (cpar == "19") {
+      LHAPDF::setVerbosity(1);
+    }
+  }
+  /// Get LHAPDF parameters -- does nothing in LHAPDF6!
+  void getlhaparm_(int dummy, char* par, int parlength) {
+    cstr_to_fstr("", par, parlength);
+  }
+
+
+  /// Return a dummy max number of sets (there is no limitation now)
+  void getmaxnumsets_(int& nmax) {
+    nmax = 1000;
+  }
+
+
+  /// Set PDF data path
+  void setpdfpath_(const char* s, size_t len) {
+    /// @todo Works? Need to check C-string copying, null termination
+    char s2[1024];
+    s2[len] = '\0';
+    strncpy(s2, s, len);
+    LHAPDF::pathsPrepend(s2);
+  }
+
+  /// Get PDF data path (colon-separated if there is more than one element)
+  void getdatapath_(char* s, size_t len) {
+    /// @todo Works? Need to check Fortran string return, string macro treatment, etc.
+    string pathstr;
+    for (const string& path : LHAPDF::paths()) {
+      if (!pathstr.empty()) pathstr += ":";
+      pathstr += path;
+    }
+    // strncpy(s, pathstr.c_str(), len);
+    cstr_to_fstr(pathstr.c_str(), s, len);
+  }
+
+
+  // PDF initialisation and focus-switching
+
+  /// Load a PDF set
+  ///
+  /// @todo Does this version actually take a *path*? What to do?
+  void initpdfsetm_(const int& nset, const char* setpath, int setpathlength) {
+    // Strip file extension for backward compatibility
+    string fullp = string(setpath, setpathlength);
+    // Remove trailing whitespace
+    fullp.erase( std::remove_if( fullp.begin(), fullp.end(), ::isspace ), fullp.end() );
+    // Use only items after the last /
+    const string pap = LHAPDF::dirname(fullp);
+    const string p = LHAPDF::basename(fullp);
+    // Prepend path to search area
+    LHAPDF::pathsPrepend(pap);
+    // Handle extensions
+    string path = LHAPDF::file_extn(p).empty() ? p : LHAPDF::file_stem(p);
+    /// @note We correct the misnamed CTEQ6L1/CTEQ6ll set name as a backward compatibility special case.
+    if (LHAPDF::to_lower(path) == "cteq6ll") path = "cteq6l1";
+    // Create the PDF set with index nset
+    // if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    if (path != ACTIVESETS[nset].setname)
+      ACTIVESETS[nset] = PDFSetHandler(path); ///< @todo Will be wrong if a structured path is given
+    CURRENTSET = nset;
+  }
+  /// Load a PDF set (non-multiset version)
+  void initpdfset_(const char* setpath, int setpathlength) {
+    int nset1 = 1;
+    initpdfsetm_(nset1, setpath, setpathlength);
+  }
+
+
+  /// Load a PDF set by name
+  void initpdfsetbynamem_(const int& nset, const char* setname, int setnamelength) {
+    // Truncate input to size setnamelength
+    string p = setname;
+    p.erase(setnamelength, std::string::npos);
+    // Strip file extension for backward compatibility
+    string name = LHAPDF::file_extn(p).empty() ? p : LHAPDF::file_stem(p);
+    // Remove trailing whitespace
+    name.erase( std::remove_if( name.begin(), name.end(), ::isspace ), name.end() );
+    /// @note We correct the misnamed CTEQ6L1/CTEQ6ll set name as a backward compatibility special case.
+    if (LHAPDF::to_lower(name) == "cteq6ll") name = "cteq6l1";
+    // Create the PDF set with index nset
+    // if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    if (name != ACTIVESETS[nset].setname)
+      ACTIVESETS[nset] = PDFSetHandler(name);
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  /// Load a PDF set by name (non-multiset version)
+  void initpdfsetbyname_(const char* setname, int setnamelength) {
+    int nset1 = 1;
+    initpdfsetbynamem_(nset1, setname, setnamelength);
+  }
+
+
+  /// Load a PDF in current set
+  void initpdfm_(const int& nset, const int& nmember) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    ACTIVESETS[nset].loadMember(nmember);
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  /// Load a PDF in current set (non-multiset version)
+  void initpdf_(const int& nmember) {
+    int nset1 = 1;
+    initpdfm_(nset1, nmember);
+  }
+
+
+  /// Get the current set number (i.e. allocation slot index)
+  void getnset_(int& nset) {
+    nset = CURRENTSET;
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  }
+
+  /// Explicitly set the current set number (i.e. allocation slot index)
+  void setnset_(const int& nset) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    CURRENTSET = nset;
+  }
+
+
+  /// Get the current member number in slot nset
+  void getnmem_(int& nset, int& nmem) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    nmem = ACTIVESETS[nset].currentmem;
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+
+  /// Set the current member number in slot nset
+  void setnmem_(const int& nset, const int& nmem) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" +
+                              LHAPDF::to_str(nset) + " but it is not initialised");
+    ACTIVESETS[nset].loadMember(nmem);
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+
+
+
+  // PDF evolution functions
+
+  // NEW BY MZ to evolve one single parton
+
+  /// Get xf(x) values for common partons from current PDF
+  void evolvepartm_(const int& nset, const int& ipart, const double& x, const double& q, double& fxq) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    int ipart_copy; // this is to deal with photons, which are labeled 7 in MG5aMC
+    ipart_copy = ipart;
+    if (ipart==7) ipart_copy = 22;
+    try {
+      fxq = ACTIVESETS[nset].activeMember()->xfxQ(ipart_copy, x, q);
+    } catch (const exception& e) {
+      fxq = 0;
+    }
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  /// Get xf(x) values for common partons from current PDF (non-multiset version)
+  void evolvepart_( const int& ipart, const double& x, const double& q, double& fxq) {
+    int nset1 = 1;
+    evolvepartm_(nset1, ipart, x, q, fxq);
+  }
+
+  /// Get xf(x) values for common partons from current PDF
+  void evolvepdfm_(const int& nset, const double& x, const double& q, double* fxq) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    // Evaluate for the 13 LHAPDF5 standard partons (-6..6)
+    for (size_t i = 0; i < 13; ++i) {
+      try {
+        fxq[i] = ACTIVESETS[nset].activeMember()->xfxQ(i-6, x, q);
+      } catch (const exception& e) {
+        fxq[i] = 0;
+      }
+    }
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  /// Get xf(x) values for common partons from current PDF (non-multiset version)
+  void evolvepdf_(const double& x, const double& q, double* fxq) {
+    int nset1 = 1;
+    evolvepdfm_(nset1, x, q, fxq);
+  }
+
+
+  /// Determine if the current PDF has a photon flavour (historically only MRST2004QED)
+  /// @todo Function rather than subroutine?
+  /// @note There is no multiset version. has_photon will respect the current set slot.
+  bool has_photon_() {
+    return ACTIVESETS[CURRENTSET].activeMember()->hasFlavor(22);
+  }
+
+
+  /// Get xfx values from current PDF, including an extra photon flavour
+  void evolvepdfphotonm_(const int& nset, const double& x, const double& q, double* fxq, double& photonfxq) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    // First evaluate the "normal" partons
+    evolvepdfm_(nset, x, q, fxq);
+    // Then evaluate the photon flavor (historically only for MRST2004QED)
+    try {
+      photonfxq = ACTIVESETS[nset].activeMember()->xfxQ(22, x, q);
+    } catch (const exception& e) {
+      photonfxq = 0;
+    }
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  /// Get xfx values from current PDF, including an extra photon flavour (non-multiset version)
+  void evolvepdfphoton_(const double& x, const double& q, double* fxq, double& photonfxq) {
+    int nset1 = 1;
+    evolvepdfphotonm_(nset1, x, q, fxq, photonfxq);
+  }
+
+
+  /// Get xf(x) values for common partons from a photon PDF
+  void evolvepdfpm_(const int& nset, const double& x, const double& q, const double& p2, const int& ip2, double& fxq) {
+    // Update current set focus
+    CURRENTSET = nset;
+    throw LHAPDF::NotImplementedError("Photon structure functions are not yet supported in LHAPDF6");
+  }
+  /// Get xf(x) values for common partons from a photon PDF (non-multiset version)
+  void evolvepdfp_(const double& x, const double& q, const double& p2, const int& ip2, double& fxq) {
+    int nset1 = 1;
+    evolvepdfpm_(nset1, x, q, p2, ip2, fxq);
+  }
+
+
+  // alpha_s evolution
+
+  /// Get the alpha_s order for the set
+  void getorderasm_(const int& nset, int& oas) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    // Set equal to the number of members for the requested set
+    oas = ACTIVESETS[nset].activeMember()->info().get_entry_as<int>("AlphaS_OrderQCD");
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  /// Get the alpha_s order for the set (non-multiset version)
+  void getorderas_(int& oas) {
+    int nset1 = 1;
+    getorderasm_(nset1, oas);
+  }
+
+
+  /// Get the alpha_s(Q) value for set nset
+  double alphaspdfm_(const int& nset, const double& Q){
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    return ACTIVESETS[nset].activeMember()->alphasQ(Q);
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  /// Get the alpha_s(Q) value for the set (non-multiset version)
+  double alphaspdf_(const double& Q){
+    int nset1 = 1;
+    return alphaspdfm_(nset1, Q);
+  }
+
+
+  // Metadata functions
+
+  /// Get the number of error members in the set
+  void numberpdfm_(const int& nset, int& numpdf) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    // Set equal to the number of members  for the requested set
+    numpdf=  ACTIVESETS[nset].activeMember()->info().get_entry_as<int>("NumMembers");
+    // Reproduce old LHAPDF v5 behaviour, i.e. subtract 1
+    numpdf -= 1;
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  /// Get the number of error members in the set (non-multiset version)
+  void numberpdf_(int& numpdf) {
+    int nset1 = 1;
+    numberpdfm_(nset1, numpdf);
+  }
+
+
+  /// Get the max number of active flavours
+  void getnfm_(const int& nset, int& nf) {
+    //nf = ACTIVESETS[nset].activeMember()->info().get_entry_as<int>("AlphaS_NumFlavors");
+    nf = ACTIVESETS[nset].activeMember()->info().get_entry_as<int>("NumFlavors");
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  /// Get the max number of active flavours (non-multiset version)
+  void getnf_(int& nf) {
+    int nset1 = 1;
+    getnfm_(nset1, nf);
+  }
+
+
+  /// Get nf'th quark mass
+  void getqmassm_(const int& nset, const int& nf, double& mass) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    if      (nf*nf ==  1) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MDown");
+    else if (nf*nf ==  4) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MUp");
+    else if (nf*nf ==  9) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MStrange");
+    else if (nf*nf == 16) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MCharm");
+    else if (nf*nf == 25) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MBottom");
+    else if (nf*nf == 36) mass = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("MTop");
+    else throw LHAPDF::UserError("Trying to get quark mass for invalid quark ID #" + LHAPDF::to_str(nf));
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  /// Get nf'th quark mass (non-multiset version)
+  void getqmass_(const int& nf, double& mass) {
+    int nset1 = 1;
+    getqmassm_(nset1, nf, mass);
+  }
+
+
+  /// Get the nf'th quark threshold
+  void getthresholdm_(const int& nset, const int& nf, double& Q) {
+    try {
+      if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+        throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+      if      (nf*nf ==  1) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdDown");
+      else if (nf*nf ==  4) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdUp");
+      else if (nf*nf ==  9) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdStrange");
+      else if (nf*nf == 16) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdCharm");
+      else if (nf*nf == 25) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdBottom");
+      else if (nf*nf == 36) Q = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("ThresholdTop");
+      //else throw LHAPDF::UserError("Trying to get quark threshold for invalid quark ID #" + LHAPDF::to_str(nf));
+    } catch (...) {
+      getqmassm_(nset, nf, Q);
+    }
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  /// Get the nf'th quark threshold
+  void getthreshold_(const int& nf, double& Q) {
+    int nset1 = 1;
+    getthresholdm_(nset1, nf, Q);
+  }
+
+
+  /// Print PDF set's description to stdout
+  void getdescm_(const int& nset) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    cout << ACTIVESETS[nset].activeMember()->description() << endl;
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  void getdesc_() {
+    int nset1 = 1;
+    getdescm_(nset1);
+  }
+
+
+  void getxminm_(const int& nset, const int& nmem, double& xmin) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    const int activemem = ACTIVESETS[nset].currentmem;
+    ACTIVESETS[nset].loadMember(nmem);
+    xmin = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("XMin");
+    ACTIVESETS[nset].loadMember(activemem);
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  void getxmin_(const int& nmem, double& xmin) {
+    int nset1 = 1;
+    getxminm_(nset1, nmem, xmin);
+  }
+
+
+  void getxmaxm_(const int& nset, const int& nmem, double& xmax) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    const int activemem = ACTIVESETS[nset].currentmem;
+    ACTIVESETS[nset].loadMember(nmem);
+    xmax = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("XMax");
+    ACTIVESETS[nset].loadMember(activemem);
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  void getxmax_(const int& nmem, double& xmax) {
+    int nset1 = 1;
+    getxmaxm_(nset1, nmem, xmax);
+  }
+
+
+  void getq2minm_(const int& nset, const int& nmem, double& q2min) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    const int activemem = ACTIVESETS[nset].currentmem;
+    ACTIVESETS[nset].loadMember(nmem);
+    q2min = LHAPDF::sqr(ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("QMin"));
+    ACTIVESETS[nset].loadMember(activemem);
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  void getq2min_(const int& nmem, double& q2min) {
+    int nset1 = 1;
+    getq2minm_(nset1, nmem, q2min);
+  }
+
+
+  void getq2maxm_(const int& nset, const int& nmem, double& q2max) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    const int activemem = ACTIVESETS[nset].currentmem;
+    ACTIVESETS[nset].loadMember(nmem);
+    q2max = LHAPDF::sqr(ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("QMax"));
+    ACTIVESETS[nset].loadMember(activemem);
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  void getq2max_(const int& nmem, double& q2max) {
+    int nset1 = 1;
+    getq2maxm_(nset1, nmem, q2max);
+  }
+
+
+  void getminmaxm_(const int& nset, const int& nmem, double& xmin, double& xmax, double& q2min, double& q2max) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    const int activemem = ACTIVESETS[nset].currentmem;
+    ACTIVESETS[nset].loadMember(nmem);
+    xmin = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("XMin");
+    xmax = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("XMax");
+    q2min = LHAPDF::sqr(ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("QMin"));
+    q2max = LHAPDF::sqr(ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("QMax"));
+    ACTIVESETS[nset].loadMember(activemem);
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  void getminmax_(const int& nmem, double& xmin, double& xmax, double& q2min, double& q2max) {
+    int nset1 = 1;
+     getminmaxm_(nset1, nmem, xmin, xmax, q2min, q2max);
+   }
+
+
+
+  void getlam4m_(const int& nset, const int& nmem, double& qcdl4) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    CURRENTSET = nset;
+    ACTIVESETS[nset].loadMember(nmem);
+    qcdl4 = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("AlphaS_Lambda4", -1.0);
+  }
+  void getlam4_(const int& nmem, double& qcdl4) {
+    int nset1 = 1;
+    getlam4m_(nset1, nmem, qcdl4);
+  }
+
+
+  void getlam5m_(const int& nset, const int& nmem, double& qcdl5) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    CURRENTSET = nset;
+    ACTIVESETS[nset].loadMember(nmem);
+    qcdl5 = ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("AlphaS_Lambda5", -1.0);
+  }
+  void getlam5_(const int& nmem, double& qcdl5) {
+    int nset1 = 1;
+    getlam5m_(nset1, nmem, qcdl5);
+  }
+
+
+
+
+
+  /// Backwards compatibility functions for LHAPDF5 calculations of
+  /// PDF uncertainties and PDF correlations (G. Watt, March 2014).
+
+  // subroutine GetPDFUncTypeM(nset,lMonteCarlo,lSymmetric)
+  void getpdfunctypem_(const int& nset, int& lmontecarlo, int& lsymmetric) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    const string errorType = ACTIVESETS[nset].activeMember()->set().errorType();
+    if (LHAPDF::startswith(errorType, "replicas")) { // Monte Carlo PDF sets
+      lmontecarlo = 1;
+      lsymmetric = 1;
+    } else if (LHAPDF::startswith(errorType, "symmhessian")) { // symmetric eigenvector PDF sets
+      lmontecarlo = 0;
+      lsymmetric = 1;
+    } else { // default: assume asymmetric Hessian eigenvector PDF sets
+      lmontecarlo = 0;
+      lsymmetric = 0;
+    }
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  // subroutine GetPDFUncType(lMonteCarlo,lSymmetric)
+  void getpdfunctype_(int& lmontecarlo, int& lsymmetric) {
+    int nset1 = 1;
+    getpdfunctypem_(nset1, lmontecarlo, lsymmetric);
+  }
+
+
+  // subroutine GetPDFuncertaintyM(nset,values,central,errplus,errminus,errsym)
+  void getpdfuncertaintym_(const int& nset, const double* values, double& central, double& errplus, double& errminus, double& errsymm) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    const size_t nmem = ACTIVESETS[nset].activeMember()->set().size()-1;
+    const vector<double> vecvalues(values, values + nmem + 1);
+    LHAPDF::PDFUncertainty err = ACTIVESETS[nset].activeMember()->set().uncertainty(vecvalues, -1);
+    central = err.central;
+    // For a combined set, the PDF and parameter variation uncertainties will be added in quadrature.
+    errplus = err.errplus;
+    errminus = err.errminus;
+    errsymm = err.errsymm;
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  // subroutine GetPDFuncertainty(values,central,errplus,errminus,errsym)
+  void getpdfuncertainty_(const double* values, double& central, double& errplus, double& errminus, double& errsymm) {
+    int nset1 = 1;
+    getpdfuncertaintym_(nset1, values, central, errplus, errminus, errsymm);
+  }
+
+
+  // subroutine GetPDFcorrelationM(nset,valuesA,valuesB,correlation)
+  void getpdfcorrelationm_(const int& nset, const double* valuesA, const double* valuesB, double& correlation) {
+    if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+      throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+    const size_t nmem = ACTIVESETS[nset].activeMember()->set().size()-1;
+    const vector<double> vecvaluesA(valuesA, valuesA + nmem + 1);
+    const vector<double> vecvaluesB(valuesB, valuesB + nmem + 1);
+    correlation = ACTIVESETS[nset].activeMember()->set().correlation(vecvaluesA,vecvaluesB);
+    // Update current set focus
+    CURRENTSET = nset;
+  }
+  // subroutine GetPDFcorrelation(valuesA,valuesB,correlation)
+  void getpdfcorrelation_(const double* valuesA, const double* valuesB, double& correlation) {
+    int nset1 = 1;
+    getpdfcorrelationm_(nset1, valuesA, valuesB, correlation);
+  }
+
+
+  ///////////////////////////////////////
+
+
+  /// REALLY OLD PDFLIB COMPATIBILITY FUNCTIONS
+
+  /// PDFLIB initialisation function
+  void pdfset_(const char* par, const double* value, int parlength) {
+
+    string my_par(par), message;
+    int id;
+    // Identify the calling program (yuck!)
+    if (my_par.find("NPTYPE") != string::npos) {
+      message = "==== LHAPDF6 USING PYTHIA-TYPE LHAGLUE INTERFACE ====";
+      // Take PDF ID from value[2]
+      id = value[2]+1000*value[1];
+    } else if (my_par.find("HWLHAPDF") != string::npos) {
+      message = "==== LHAPDF6 USING HERWIG-TYPE LHAGLUE INTERFACE ====";
+      // Take PDF ID from value[0]
+      id = value[0];
+    } else if (my_par.find("DEFAULT") != string::npos) {
+      message = "==== LHAPDF6 USING DEFAULT-TYPE LHAGLUE INTERFACE ====";
+      // Take PDF ID from value[0]
+      id = value[0];
+    } else {
+      message = "==== LHAPDF6 USING PDFLIB-TYPE LHAGLUE INTERFACE ====";
+      // Take PDF ID from value[2]
+      id = value[2]+1000*value[1];
+    }
+    pair<string, int> set_id = LHAPDF::lookupPDF(id);
+    if (set_id.first != ACTIVESETS[1].setname || set_id.second != ACTIVESETS[1].currentmem) {
+      if (LHAPDF::verbosity() > 0) cout << message << endl;
+      ACTIVESETS[1] = PDFSetHandler(id);
+    }
+
+    CURRENTSET = 1;
+
+    // Extract parameters for common blocks (with sensible fallback values)
+    PDFPtr pdf = ACTIVESETS[1].activeMember();
+    w50513_.xmin = pdf->info().get_entry_as<double>("XMin", 0.0);
+    w50513_.xmax = pdf->info().get_entry_as<double>("XMax", 1.0);
+    w50513_.q2min = LHAPDF::sqr(pdf->info().get_entry_as<double>("QMin", 1.0));
+    w50513_.q2max = LHAPDF::sqr(pdf->info().get_entry_as<double>("QMax", 1.0e5));
+    w50512_.qcdl4 = pdf->info().get_entry_as<double>("AlphaS_Lambda4", 0.0);
+    w50512_.qcdl5 = pdf->info().get_entry_as<double>("AlphaS_Lambda5", 0.0);
+    lhapdfr_.qcdlha4 = pdf->info().get_entry_as<double>("AlphaS_Lambda4", 0.0);
+    lhapdfr_.qcdlha5 = pdf->info().get_entry_as<double>("AlphaS_Lambda5", 0.0);
+    lhapdfr_.nfllha = 4;
+    // Activate legacy/compatibility LHAPDF5-type behaviour re. broken Lambda values
+    if (pdf->info().get_entry_as<bool>("Pythia6LambdaV5Compat", true)) {
+      w50512_.qcdl4 = 0.192;
+      w50512_.qcdl5 = 0.192;
+      lhapdfr_.qcdlha4 = 0.192;
+      lhapdfr_.qcdlha5 = 0.192;
+    }
+  }
+
+  /// PDFLIB nucleon structure function querying
+  void structm_(const double& x, const double& q,
+                double& upv, double& dnv, double& usea, double& dsea,
+                double& str, double& chm, double& bot, double& top, double& glu) {
+    CURRENTSET = 1;
+    /// Fill (partial) parton return variables
+    PDFPtr pdf = ACTIVESETS[1].activeMember();
+    dsea = pdf->xfxQ(-1, x, q);
+    usea = pdf->xfxQ(-2, x, q);
+    dnv = pdf->xfxQ(1, x, q) - dsea;
+    upv = pdf->xfxQ(2, x, q) - usea;
+    str = pdf->xfxQ(3, x, q);
+    chm = (pdf->hasFlavor(4)) ? pdf->xfxQ(4, x, q) : 0;
+    bot = (pdf->hasFlavor(5)) ? pdf->xfxQ(5, x, q) : 0;
+    top = (pdf->hasFlavor(6)) ? pdf->xfxQ(6, x, q) : 0;
+    glu = pdf->xfxQ(21, x, q);
+  }
+
+  /// PDFLIB photon structure function querying
+  void structp_(const double& x, const double& q2, const double& p2, const double& ip2,
+                double& upv, double& dnv, double& usea, double& dsea,
+                double& str, double& chm, double& bot, double& top, double& glu) {
+    throw LHAPDF::NotImplementedError("Photon structure functions are not yet supported");
+  }
+
+  /// PDFLIB statistics on PDF under/overflows
+  void pdfsta_() {
+    /// @note Can't do anything...
+  }
+
+
+}
+
+
+// LHAPDF namespace C++ compatibility code
+#ifdef ENABLE_LHAGLUE_CXX
+
+
+void LHAPDF::setVerbosity(LHAPDF::Verbosity noiselevel) {
+  LHAPDF::setVerbosity((int) noiselevel);
+}
+
+void LHAPDF::setPDFPath(const string& path) {
+  pathsPrepend(path);
+}
+
+string LHAPDF::pdfsetsPath() {
+  return paths()[0];
+}
+
+int LHAPDF::numberPDF() {
+  int nmem;
+  numberpdf_(nmem);
+  return nmem;
+}
+int LHAPDF::numberPDF(int nset) {
+  int nmem;
+  numberpdfm_(nset,nmem);
+  return nmem;
+}
+
+void LHAPDF::initPDF( int memset) {
+  int nset1 = 1;
+  initpdfm_(nset1, memset);
+}
+void LHAPDF::initPDF(int nset, int memset) {
+  initpdfm_(nset, memset);
+}
+
+
+double LHAPDF::xfx(double x, double Q, int fl) {
+  vector<double> r(13);
+  evolvepdf_(x, Q, &r[0]);
+  return r[fl+6];
+}
+double LHAPDF::xfx(int nset, double x, double Q, int fl) {
+  vector<double> r(13);
+  evolvepdfm_(nset, x, Q, &r[0]);
+  return r[fl+6];
+}
+
+vector<double> LHAPDF::xfx(double x, double Q) {
+  vector<double> r(13);
+  evolvepdf_(x, Q, &r[0]);
+  return r;
+}
+vector<double> LHAPDF::xfx(int nset, double x, double Q) {
+  vector<double> r(13);
+  evolvepdfm_(nset, x, Q, &r[0]);
+  return r;
+}
+
+void LHAPDF::xfx(double x, double Q, double* results) {
+  evolvepdf_(x, Q, results);
+}
+void LHAPDF::xfx(int nset, double x, double Q, double* results) {
+  evolvepdfm_(nset, x, Q, results);
+}
+
+
+vector<double> LHAPDF::xfxphoton(double x, double Q) {
+  vector<double> r(13);
+  double mphoton;
+  evolvepdfphoton_(x, Q, &r[0], mphoton);
+  r.push_back(mphoton);
+  return r;
+}
+vector<double> LHAPDF::xfxphoton(int nset, double x, double Q) {
+  vector<double> r(13);
+  double mphoton;
+  evolvepdfphotonm_(nset, x, Q, &r[0], mphoton);
+  r.push_back(mphoton);
+  return r;
+}
+
+void LHAPDF::xfxphoton(double x, double Q, double* results) {
+  evolvepdfphoton_(x, Q, results, results[13]);
+}
+void LHAPDF::xfxphoton(int nset, double x, double Q, double* results) {
+  evolvepdfphotonm_(nset, x, Q, results, results[13]);
+}
+
+double LHAPDF::xfxphoton(double x, double Q, int fl) {
+  vector<double> r(13);
+  double mphoton;
+  evolvepdfphoton_(x, Q, &r[0], mphoton);
+  if (fl == 7) return mphoton;
+  return r[fl+6];
+}
+double LHAPDF::xfxphoton(int nset, double x, double Q, int fl) {
+  vector<double> r(13);
+  double mphoton;
+  evolvepdfphotonm_(nset, x, Q, &r[0], mphoton);
+  if ( fl == 7 ) return mphoton;
+  return r[fl+6];
+}
+
+
+void LHAPDF::initPDFSet(const string& filename, int nmem) {
+  initPDFSet(1,filename, nmem);
+}
+
+void LHAPDF::initPDFSet(int nset, const string& filename, int nmem) {
+  initPDFSetByName(nset,filename);
+  ACTIVESETS[nset].loadMember(nmem);
+  CURRENTSET = nset;
+}
+
+
+void LHAPDF::initPDFSet(const string& filename, SetType type, int nmem) {
+  // silently ignore type
+  initPDFSet(1,filename, nmem);
+}
+
+void LHAPDF::initPDFSet(int nset, const string& filename, SetType type, int nmem) {
+  // silently ignore type
+  initPDFSetByName(nset,filename);
+  ACTIVESETS[nset].loadMember(nmem);
+  CURRENTSET = nset;
+}
+
+void LHAPDF::initPDFSet(int nset, int setid, int nmem) {
+  pair<string, int> set_id = LHAPDF::lookupPDF(setid+nmem);
+  if (set_id.second != nmem)
+    throw LHAPDF::UserError("Inconsistent member numbers: " + LHAPDF::to_str(set_id.second) + " != " + LHAPDF::to_str(nmem));
+  if (set_id.first != ACTIVESETS[nset].setname || nmem != ACTIVESETS[nset].currentmem)
+      ACTIVESETS[nset] = PDFSetHandler(setid+nmem);
+  CURRENTSET = nset;
+}
+
+void LHAPDF::initPDFSet(int setid, int nmem) {
+  initPDFSet(1,setid,nmem);
+}
+
+#define SIZE 999
+void LHAPDF::initPDFSetByName(const string& filename) {
+  std::cout << "initPDFSetByName: " << filename << std::endl;
+  char cfilename[SIZE+1];
+  strncpy(cfilename, filename.c_str(), SIZE);
+  initpdfsetbyname_(cfilename, filename.length());
+}
+
+void LHAPDF::initPDFSetByName(int nset, const string& filename) {
+  char cfilename[SIZE+1];
+  strncpy(cfilename, filename.c_str(), SIZE);
+  initpdfsetbynamem_(nset, cfilename, filename.length());
+}
+
+void LHAPDF::initPDFSetByName(const string& filename, SetType type) {
+  //silently ignore type
+  std::cout << "initPDFSetByName: " << filename << std::endl;
+  char cfilename[SIZE+1];
+  strncpy(cfilename, filename.c_str(), SIZE);
+  initpdfsetbyname_(cfilename, filename.length());
+}
+
+void LHAPDF::initPDFSetByName(int nset, const string& filename, SetType type) {
+  //silently ignore type
+  char cfilename[SIZE+1];
+  strncpy(cfilename, filename.c_str(), SIZE);
+  initpdfsetbynamem_(nset, cfilename, filename.length());
+}
+
+
+void LHAPDF::getDescription() {
+  getDescription(1);
+}
+
+void LHAPDF::getDescription(int nset) {
+  if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  cout << ACTIVESETS[nset].activeMember()->set().description() << endl;
+}
+
+
+double LHAPDF::alphasPDF(double Q) {
+  return LHAPDF::alphasPDF(1, Q) ;
+}
+
+double LHAPDF::alphasPDF(int nset, double Q) {
+  if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  CURRENTSET = nset;
+  // return alphaS for the requested set
+  return ACTIVESETS[nset].activeMember()->alphasQ(Q);
+}
+
+
+bool LHAPDF::hasPhoton(){
+  return has_photon_();
+}
+
+
+int LHAPDF::getOrderAlphaS() {
+  return LHAPDF::getOrderAlphaS(1) ;
+}
+
+int LHAPDF::getOrderAlphaS(int nset) {
+  if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  CURRENTSET = nset;
+  // return alphaS Order for the requested set
+  return ACTIVESETS[nset].activeMember()->info().get_entry_as<int>("AlphaS_OrderQCD", -1);
+}
+
+
+int LHAPDF::getOrderPDF() {
+  return LHAPDF::getOrderPDF(1) ;
+}
+
+int LHAPDF::getOrderPDF(int nset) {
+  if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  CURRENTSET = nset;
+  // return PDF order for the requested set
+  return ACTIVESETS[nset].activeMember()->info().get_entry_as<int>("OrderQCD", -1);
+}
+
+
+double LHAPDF::getLam4(int nmem) {
+  return LHAPDF::getLam4(1, nmem) ;
+}
+
+double LHAPDF::getLam4(int nset, int nmem) {
+  // if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //   throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  // CURRENTSET = nset;
+  // ACTIVESETS[nset].loadMember(nmem);
+  // return ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("AlphaS_Lambda4", -1.0);
+  double qcdl4;
+  getlam4m_(nset, nmem, qcdl4);
+  return qcdl4;
+}
+
+
+double LHAPDF::getLam5(int nmem) {
+  return LHAPDF::getLam5(1, nmem) ;
+}
+
+double LHAPDF::getLam5(int nset, int nmem) {
+  // if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+  //   throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  // CURRENTSET = nset;
+  // ACTIVESETS[nset].loadMember(nmem);
+  // return ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("AlphaS_Lambda5", -1.0);
+  double qcdl5;
+  getlam5m_(nset, nmem, qcdl5);
+  return qcdl5;
+}
+
+
+int LHAPDF::getNf() {
+  return LHAPDF::getNf(1) ;
+}
+
+int LHAPDF::getNf(int nset) {
+  if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  CURRENTSET = nset;
+  // return alphaS Order for the requested set
+  return ACTIVESETS[nset].activeMember()->info().get_entry_as<int>("NumFlavors");
+}
+
+
+double LHAPDF::getXmin(int nmem) {
+  return LHAPDF::getXmin(1, nmem) ;
+}
+
+double LHAPDF::getXmin(int nset, int nmem) {
+  if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  CURRENTSET = nset;
+  // return alphaS Order for the requested set
+  ACTIVESETS[nset].loadMember(nmem);
+  return ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("XMin");
+}
+
+double LHAPDF::getXmax(int nmem) {
+  return LHAPDF::getXmax(1, nmem) ;
+}
+
+double LHAPDF::getXmax(int nset, int nmem) {
+  if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  CURRENTSET = nset;
+  // return alphaS Order for the requested set
+  ACTIVESETS[nset].loadMember(nmem);
+  return ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("XMax");
+}
+
+double LHAPDF::getQ2min(int nmem) {
+  return LHAPDF::getQ2min(1, nmem) ;
+}
+
+double LHAPDF::getQ2min(int nset, int nmem) {
+  if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  CURRENTSET = nset;
+  // return alphaS Order for the requested set
+  ACTIVESETS[nset].loadMember(nmem);
+  return pow(ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("QMin"),2);
+}
+
+double LHAPDF::getQ2max(int nmem) {
+  return LHAPDF::getQ2max(1,nmem) ;
+}
+
+double LHAPDF::getQ2max(int nset, int nmem) {
+  if (ACTIVESETS.find(nset) == ACTIVESETS.end())
+    throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
+  CURRENTSET = nset;
+  // return alphaS Order for the requested set
+  ACTIVESETS[nset].loadMember(nmem);
+  return pow(ACTIVESETS[nset].activeMember()->info().get_entry_as<double>("QMax"),2);
+}
+
+double LHAPDF::getQMass(int nf) {
+  return LHAPDF::getQMass(1, nf) ;
+}
+
+double LHAPDF::getQMass(int nset, int nf) {
+  double mass;
+  getqmassm_(nset, nf, mass);
+  return mass;
+}
+
+double LHAPDF::getThreshold(int nf) {
+  return LHAPDF::getThreshold(1, nf) ;
+}
+
+double LHAPDF::getThreshold(int nset, int nf) {
+  double thres;
+  getthresholdm_(nset, nf, thres);
+  return thres;
+}
+
+void LHAPDF::usePDFMember(int member) {
+  initpdf_(member);
+}
+
+void LHAPDF::usePDFMember(int nset, int member) {
+  initpdfm_(nset, member);
+}
+
+#endif // ENABLE_LHAGLUE_CXX

=== modified file 'Template/LO/Source/dsample.f'
--- Template/LO/Source/dsample.f	2020-01-20 14:53:46 +0000
+++ Template/LO/Source/dsample.f	2020-07-04 19:03:17 +0000
@@ -141,7 +141,7 @@
                fx =0d0
                wgt=0d0
             endif
-            call sample_put_point(wgt,x(1),iter,ipole,itmin) !Store result
+            call sample_put_point(wgt,x(1),iter,ipole) !Store result
          endif
          if (wgt .ne. 0d0) kevent=kevent+1    
 c
@@ -318,7 +318,7 @@
             endif
             
             if (nzoom .le. 0) then
-               call sample_put_point(wgt,x(1),iter,ipole,itmin) !Store result
+               call sample_put_point(wgt,x(1),iter,ipole) !Store result
             else
                nzoom = nzoom -1
                ievent=ievent-1

=== modified file 'Template/LO/SubProcesses/cuts.f'
--- Template/LO/SubProcesses/cuts.f	2019-07-30 13:51:02 +0000
+++ Template/LO/SubProcesses/cuts.f	2020-07-04 19:03:17 +0000
@@ -1394,10 +1394,10 @@
 C
     2 IF (N.EQ.1)            RETURN
       IF (MODE)    10,20,30
-   10 CALL SORTTI (A,INDEX,N)
+   10 STOP 5 ! CALL SORTTI (A,INDEX,N)
       GO TO 40
 C
-   20 CALL SORTTC(A,INDEX,N)
+   20 STOP 5 ! CALL SSORTTC(A,INDEX,N)
       GO TO 40
 C
    30 CALL SORTTF (A,INDEX,N)

=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2019-09-26 11:09:58 +0000
+++ Template/LO/SubProcesses/genps.f	2020-07-04 19:03:17 +0000
@@ -666,7 +666,8 @@
 c      data nerr/0/
       double precision smin,smax,totmass,totmassin,xa2,xb2,wgt
       double precision costh,phi,tmin,tmax,t
-      double precision ma2,mb2,m12,mn2,s1
+      double precision ma2,mb2,m12,mn2,s1, mi2
+      double precision tmass(-max_branch:-1)
 c
 c     External
 c
@@ -864,7 +865,7 @@
 c         and not to be block by numerical inacuracy
 c         tmax = max(tmax,0d0) !This line if want really t freedom
          call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
-     $        0, -tmin/stot)
+     $        0d0, -tmin/stot)
          t = stot*(-x(-ibranch))
 
       else
@@ -893,7 +894,14 @@
 c     pa+pb -> p1+ p2; t=(pa-p1)**2;   pr = pa-p1
 c     gentcms(pa,pb,t,phi,m1,m2,p1,pr) 
 c
-         call gentcms(p(0,itree(1,ibranch)),p(0,2),t,phi,
+
+         if (itree(1,ibranch).gt.-ns_channel-1)then
+            mi2 = m(itree(1,ibranch))**2
+         else
+            mi2 = tmass(itree(1,ibranch))
+         endif
+         tmass(ibranch) = t
+         call gentcms(p(0,itree(1,ibranch)),p(0,2),t,phi, mi2,
      &        m(itree(2,ibranch)),m(ibranch-1),p(0,itree(2,ibranch)),
      &        p(0,ibranch),jac)
 
@@ -1002,7 +1010,7 @@
       endif
       end
 
-      subroutine gentcms(pa,pb,t,phi,m1,m2,p1,pr,jac)
+      subroutine gentcms(pa,pb,t,phi,ma2,m1,m2,p1,pr,jac)
 c*************************************************************************
 c     Generates 4 momentum for particle 1, and remainder pr
 c     given the values t, and phi
@@ -1038,7 +1046,7 @@
             ptotm(i) = ptot(i)
          endif
       enddo
-      ma2 = dot(pa,pa)
+c      ma2 = dot(pa,pa)
 c
 c     determine magnitude of p1 in cms frame (from dhelas routine mom2cx)
 c

=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f	2020-03-08 20:49:57 +0000
+++ Template/LO/SubProcesses/reweight.f	2020-07-04 19:03:17 +0000
@@ -417,6 +417,12 @@
           ipart(1,imo)=ipart(1,ida2)
           ipart(2,imo)=ipart(1,ida1)
         endif
+      else if (abs(get_color(idmo)).eq.8.and.abs(get_color(idda1)).eq.1.and.abs(get_color(idda2)).eq.8)then
+         ipart(1,imo)=ipart(1,ida2)
+         ipart(2,imo)=ipart(2,ida2)
+      else if (abs(get_color(idmo)).eq.8.and.abs(get_color(idda1)).eq.8.and.abs(get_color(idda2)).eq.1)then
+         ipart(1,imo)=ipart(1,ida1)
+         ipart(2,imo)=ipart(2,ida1)
       else
          write(*,*) idmo,'>', idda1, idda2, 'color', get_color(idmo),'>', get_color(idda1), get_color(idda2)
          write(*,*) "failed for ipartupdate. Please retry without MLM/default dynamical scale"

=== modified file 'Template/LO/SubProcesses/unwgt.f'
--- Template/LO/SubProcesses/unwgt.f	2018-06-20 11:54:42 +0000
+++ Template/LO/SubProcesses/unwgt.f	2020-07-04 19:03:17 +0000
@@ -478,7 +478,8 @@
       double precision beta, get_betaz
       double precision ebi(0:3), ebo(0:3)
       double precision ptcltmp(nexternal), pdum(0:3)
-
+      integer beam_number
+      
       integer idup(nexternal,maxproc,maxsproc)
       integer mothup(2,nexternal)
       integer icolup(2,nexternal,maxflow,maxsproc)
@@ -705,31 +706,42 @@
          else
             write(s_buff(3), '(a)') '<asrwt>0</asrwt>'
          endif
-         if(n_pdfrw(1).gt.0)then
+         beam_number = 1
+         if (flip) then
+            beam_number =2
+         endif
+         
+         if(n_pdfrw(1).gt.0.and.abs(lpp(1)).eq.1)then
             if(2*n_pdfrw(1).lt.10) then
-               write(cfmt,'(a,I1,a,I1,a)') '(a,I3,',
+               write(cfmt,'(a,I1,a,I1,a)') '(a,I1,a,I3,',
      $              n_pdfrw(1),'I9,',2*n_pdfrw(1),'E15.8,a)'
             else
-               write(cfmt,'(a,I1,a,I2,a)') '(a,I3,',
+               write(cfmt,'(a,I1,a,I2,a)') '(a,I1,a,I3,',
      $              n_pdfrw(1),'I9,',2*n_pdfrw(1),'E15.8,a)'
             endif
-            write(s_buff(4), cfmt) '<pdfrwt beam="1">',
+            
+            write(s_buff(4), cfmt) '<pdfrwt beam="', beam_number, '">',
      $           n_pdfrw(1),(i_pdgpdf(i,1),i=1,n_pdfrw(1)),
      $           (s_xpdf(i,1),i=1,n_pdfrw(1)),
      $           (s_qpdf(i,1),i=1,n_pdfrw(1)),
      $           '</pdfrwt>'
          else
-            write(s_buff(4), '(a)') '<pdfrwt beam="1">0</pdfrwt>'
-         endif
-         if(n_pdfrw(2).gt.0)then
+            write(s_buff(4), '(a,I1,a)') '<pdfrwt beam="',
+     $       beam_number,'">0</pdfrwt>'
+         endif
+         beam_number = 2
+         if (flip) then
+            beam_number	= 1
+         endif
+         if(n_pdfrw(2).gt.0.and.abs(lpp(2)).eq.1)then
             if(2*n_pdfrw(2).lt.10) then
-               write(cfmt,'(a,I1,a,I1,a)') '(a,I3,',
+               write(cfmt,'(a,I1,a,I1,a)') '(a,I1,a,I3,',
      $              n_pdfrw(2),'I9,',2*n_pdfrw(2),'E15.8,a)'
             else
-               write(cfmt,'(a,I1,a,I2,a)') '(a,I3,',
+               write(cfmt,'(a,I1,a,I2,a)') '(a,I1,a,I3,',
      $              n_pdfrw(2),'I9,',2*n_pdfrw(2),'E15.8,a)'
             endif
-            write(s_buff(5), cfmt) '<pdfrwt beam="2">',
+            write(s_buff(5), cfmt) '<pdfrwt beam="',beam_number,'">',
      $           n_pdfrw(2),(i_pdgpdf(i,2),i=1,n_pdfrw(2)),
      $           (s_xpdf(i,2),i=1,n_pdfrw(2)),
      $           (s_qpdf(i,2),i=1,n_pdfrw(2)),

=== modified file 'Template/LO/bin/internal/merge.pl'
--- Template/LO/bin/internal/merge.pl	2014-08-07 12:08:24 +0000
+++ Template/LO/bin/internal/merge.pl	2020-07-04 19:03:17 +0000
@@ -52,8 +52,8 @@
   # Keep track if we are in the init block or not
   $initblock = 0;
 
-  # LHE version extracted from current file ; 1 by default
-  $lhe_version = 1.0;
+  # LHE version extracted from current file ; 3 by default
+  $lhe_version = 3.0;
 
   while (1) {
     $gzbytes = $gzin->gzreadline($gzline);

=== modified file 'Template/LO/bin/internal/restore_data'
--- Template/LO/bin/internal/restore_data	2011-06-29 14:15:21 +0000
+++ Template/LO/bin/internal/restore_data	2020-07-04 19:03:17 +0000
@@ -30,13 +30,21 @@
     exit
 fi
 
-
-cp  $1_results.dat results.dat  >& /dev/null
+if [[ -e $1_results.dat ]]; then
+    cp  $1_results.dat results.dat  >& /dev/null
+else
+    cp  results.dat $1_results.dat  >& /dev/null
+fi
+    
 for i in `cat subproc.mg` ; do
     cd $i
     echo $i
     rm -f ftn25 ftn26 >& /dev/null
-    cp $1_results.dat results.dat  >& /dev/null
+    if [[ -e $1_results.dat ]]; then
+	cp  $1_results.dat results.dat  >& /dev/null
+    else
+	cp  results.dat $1_results.dat  >& /dev/null
+    fi
     for k in G* ; do
 	if [[ ! -d $k ]]; then
 	    continue
@@ -45,6 +53,8 @@
 	for j in $1_results.dat ; do
 	    if [[ -e $j ]] ; then
 		cp  $j results.dat
+	    else
+		cp  results.dat $j
 	    fi
 	done
 	for j in $1_ftn26.gz ; do

=== modified file 'Template/NLO/MCatNLO/include/LHEF.h'
--- Template/NLO/MCatNLO/include/LHEF.h	2014-10-31 07:52:06 +0000
+++ Template/NLO/MCatNLO/include/LHEF.h	2020-07-04 19:03:17 +0000
@@ -1971,7 +1971,7 @@
    * Used internally to read a single line from the stream.
    */
   bool getline() {
-    return ( std::getline(file, currentLine) );
+    return ( (bool) std::getline(file, currentLine) );
   }
 
 protected:

=== modified file 'Template/NLO/SubProcesses/cuts.f'
--- Template/NLO/SubProcesses/cuts.f	2020-05-23 19:29:30 +0000
+++ Template/NLO/SubProcesses/cuts.f	2020-07-04 19:03:17 +0000
@@ -705,10 +705,10 @@
 C
     2 IF (N.EQ.1)            RETURN
       IF (MODE)    10,20,30
-   10 CALL SORTTI (A,INDEX,N)
+   10 STOP 5 ! CALL SORTTI (A,INDEX,N)
       GO TO 40
 C
-   20 CALL SORTTC(A,INDEX,N)
+   20 STOP 5 ! CALL SORTTC(A,INDEX,N)
       GO TO 40
 C
    30 CALL SORTTF (A,INDEX,N)

=== modified file 'Template/NLO/SubProcesses/splitorders_stuff.f'
--- Template/NLO/SubProcesses/splitorders_stuff.f	2016-02-10 13:36:19 +0000
+++ Template/NLO/SubProcesses/splitorders_stuff.f	2020-07-04 19:03:17 +0000
@@ -33,13 +33,35 @@
       integer i,j
       integer base, step
       parameter(base=100)
+      ! this is for the printout of the informations
+      logical firsttime, firsttime_contr(amp_split_size)
+      data firsttime/.true./
+      data firsttime_contr/amp_split_size * .true./
+      integer orders_to_amp_split_pos
+
+      ! print out some extra informations
+      if (firsttime) write(*,fmt='(a)',advance="NO") 
+     $    "INFO: orders_tag_plot is computed as:"
 
       get_orders_tag=0
       step=1
       do i =1, nsplitorders
+        if (firsttime) write(*,fmt='(3a,i8)',advance="NO") 
+     $      "         + ", ordernames(i), " * ", step
         get_orders_tag=get_orders_tag+step*ord(i)
         step=step*100
       enddo
+      if (firsttime) then
+        write(*,*)
+        firsttime=.false.
+      endif
+
+      if (firsttime_contr(orders_to_amp_split_pos(ord))) then
+        write(*,*) 'orders_tag_plot= ', get_orders_tag, ' for ',
+     #     (ordernames(i),",",i=1,nsplitorders), ' = ',
+     #     (ord(i),",",i=1,nsplitorders)
+        firsttime_contr(orders_to_amp_split_pos(ord)) = .false.
+      endif
 
       return 
       end

=== modified file 'Template/NLO/Utilities/collect_plots.py'
--- Template/NLO/Utilities/collect_plots.py	2017-02-23 15:15:43 +0000
+++ Template/NLO/Utilities/collect_plots.py	2020-07-04 19:03:17 +0000
@@ -1,9 +1,10 @@
 import pickle
 import subprocess
+import sys
 
 def combine_plots_HwU(jobs,out,normalisation=None):
     """Sums all the plots in the HwU format."""
-    command =  []
+    command =  [sys.executable]
     command.append('bin/internal/histograms.py')
     for job in jobs:
         if job['dirname'].endswith('.HwU'):

=== removed directory 'Template/RWGTNLO'
=== removed file 'Template/RWGTNLO/__init__.py'
=== removed file 'Template/RWGTNLO/alfas.inc'
=== removed file 'Template/RWGTNLO/alfas_functions_lhapdf.f'
=== removed file 'Template/RWGTNLO/makefile'
=== removed file 'Template/RWGTNLO/rwgt.f'
=== removed file 'Template/RWGTNLO/setrun.f'
=== removed file 'Template/RWGTNLO/timing_variables.inc'
=== modified file 'Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.inc'
--- Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.inc	2016-11-21 21:05:57 +0000
+++ Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.inc	2020-07-04 19:03:17 +0000
@@ -141,7 +141,8 @@
 
       character(512) path
       character(512) dummy      
-
+      character(512) epath ! path of the executable
+      integer pos
       character(512) prefix,fpath
       character(17) nameToCheck
       parameter (nameToCheck='MadLoopParams.dat')
@@ -185,11 +186,42 @@
           close(1)
           prefix='../MadLoop5_resources/'
           call joinPath(prefix,nameToCheck,fpath)
-          OPEN(1, FILE=fpath, ERR=66, STATUS='OLD',ACTION='READ')
-          MLPath=prefix
-          goto 10
-66        continue
-          close(1)
+          OPEN(1, FILE=fpath, ERR=3, STATUS='OLD',ACTION='READ')
+          MLPath=prefix
+          goto 10
+ 3        continue
+          close(1)
+c
+c     Try to automatically find the path from the executable location
+c     particularly usefull in gridpack readonly mode
+c
+	  call getarg(0,path) !path is the PATH to the madevent executable (either global or from launching directory)
+	  pos = index(path,'/',.true.)
+	  prefix = path(:pos)
+          call joinPath(prefix,nameToCheck,fpath)
+	  write(*,*) 'test', fpath
+	  OPEN(1, FILE=fpath, ERR=4, STATUS='OLD',ACTION='READ')
+          MLPath=prefix
+          goto 10
+4         continue
+          close(1)
+          prefix= prefix // '/MadLoop5_resources/'
+          call joinPath(prefix,nameToCheck,fpath)
+	  write(*,*) 'test', fpath
+          OPEN(1, FILE=fpath, ERR=5, STATUS='OLD',ACTION='READ')
+          MLPath=prefix
+          goto 10
+5         continue
+          close(1)
+          prefix= path(:pos) // '/../MadLoop5_resources/'
+          call joinPath(prefix,nameToCheck,fpath)
+	  write(*,*) 'test', fpath
+ 	  OPEN(1, FILE=fpath, ERR=6, STATUS='OLD',ACTION='READ')
+          MLPath=prefix
+          goto 10
+ 6        continue
+          close(1)
+
 c     We could not automatically find the auxiliary files
           write(*,*) '==='
           write(*,*) 'ERROR: MadLoop5 could not automatically find the file MadLoopParams.dat.'
@@ -215,9 +247,9 @@
 
 C     Check that the FilePath set is correct
       call joinPath(MLPath,nameToCheck,fpath)
-      OPEN(1, FILE=fpath, ERR=3, STATUS='OLD',ACTION='READ')
+      OPEN(1, FILE=fpath, ERR=33, STATUS='OLD',ACTION='READ')
       goto 11
-3     continue
+33     continue
       close(1)
       write(*,*) '==='
       write(*,*) 'ERROR: The MadLoop5 auxiliary files could not be found in ',MLPath
@@ -606,7 +638,7 @@
       IMPLICIT  NONE
       INTEGER MAXNREF_EVALS
       PARAMETER (MAXNREF_EVALS=100)
-      INTEGER, DIMENSION(MAXNREF_EVALS), INTENT(IN) :: x
+      Double Precision, DIMENSION(MAXNREF_EVALS), INTENT(IN) :: x
       INTEGER, INTENT(IN)							  :: mStart, mEnd
       INTEGER										  :: Minimum
       INTEGER										  :: Location

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2020-05-25 18:16:07 +0000
+++ UpdateNotes.txt	2020-07-04 19:03:17 +0000
@@ -1,5 +1,7 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
+3.0.3 (21/06/20)
+     include up to 2.7.3
 
 3.0.2 (25/05/20)
      include up to 2.7.2
@@ -34,6 +36,26 @@
      RF: Reduction of the number of jobs -- results in earlier detection of small integration channels for fNLO
      OM: remove ./bin/mg5 executable
 
+2.7.3(21/06/20)
+      OM: Fixing some bug for read-only LO gridpacks (wrong cross-section and shape when generating events).
+          Thanks to Congqiao Li for this
+      OM: Allowing loop-induced process to run on LO gridpack with read-only mode.
+          Thanks to Congqiao Li for this      
+      OM: Fix a bug in the longitudinal polarization for off-shell effect, leading to deviation at large invariant mass.
+      OM: Adding more option to the run_card for fine tuning phase-space integration steps
+          All are hidden by default:
+	    - hard_survey [default=1]: request for more points in survey (and subsequent refine)
+	    - second_refine_treshold [default=1.5]: forbid second refine if cross section after first refine is
+	       is smaller than cross-section of the survey times such treshold
+      OM: new command for the editions of the cards:
+           - set nodecay: remove all decay line from the madspin_card
+	   - set BLOCKNAME all VALUE: set all entry of the param_card "BLOCKNAME" to VALUE
+	   - edit CARDNAME --comment_line='<regular_expression>' : new syntax to comment all lines of a card
+	       that are matching a given regular expression
+      OM: For Mac only, when running in script mode, MG5aMC will now prevent the computer to go to idle sleep.
+          You can prevent this by running with the '-s' option. like ./bin/mg5_aMC -s PATH_TO_CMD
+>>>>>>> MERGE-SOURCE
+
 2.7.2(17/03/20)
       OM: Fix a Bug in pythia8 running on Ubuntu 18.04.4 machine
       OM: Speed up standalone_cpp code by changing compilation flag

=== modified file 'VERSION'
--- VERSION	2020-05-25 18:16:07 +0000
+++ VERSION	2020-07-04 19:03:17 +0000
@@ -1,6 +1,5 @@
-version = 3.0.2
-date = 2020-05-25
-
+version = 3.0.3
+date = 2020-06-21
 
 
 

=== modified file 'aloha/create_aloha.py'
--- aloha/create_aloha.py	2020-03-11 09:28:14 +0000
+++ aloha/create_aloha.py	2020-07-04 19:03:17 +0000
@@ -442,7 +442,7 @@
             denominator = propagator.denominator      
         elif propa == "1L":
             numerator = "EPSL(1,id) * EPSL(2,id)"
-            denominator = "-1*PVec(-2,id)*PVec(-2,id)*Mass(id)**2 * (P(-1,id)**2 - Mass(id) * Mass(id) + complex(0,1) * Mass(id) * Width(id))"
+            denominator = "-1*PVec(-2,id)*PVec(-2,id)*P(-3,id)*P(-3,id) * (P(-1,id)**2 - Mass(id) * Mass(id) + complex(0,1) * Mass(id) * Width(id))"
         elif propa == "1T":
             numerator = "-1*PVec(-2,id)*PVec(-2,id) * EPST2(1,id)*EPST2(2,id) + EPST1(1,id)*EPST1(2,id)"
             denominator = "PVec(-2,id)*PVec(-2,id) * PT(-3,id)*PT(-3,id) * (P(-1,id)**2 - Mass(id) * Mass(id) + complex(0,1) * Mass(id) * Width(id))"

=== modified file 'aloha/template_files/aloha_functions_loop.f'
--- aloha/template_files/aloha_functions_loop.f	2016-01-27 01:59:42 +0000
+++ aloha/template_files/aloha_functions_loop.f	2020-07-04 19:03:17 +0000
@@ -2298,7 +2298,7 @@
       COMPLEX*16 Q(0:3)
       INTEGER CFIG,J
       LOGICAL SCD
-      COMPLEX*16 M
+      double precision M
       COMPLEX*16 W(20)
 
       IF (CFIG.EQ.1) THEN
@@ -2558,7 +2558,7 @@
       COMPLEX*16 Q(0:3)
       INTEGER CFIG
       LOGICAL SCD
-      COMPLEX*16 M
+      double precision M
       COMPLEX*16 W(20)
 
       IF (CFIG.EQ.1) THEN

=== modified file 'bin/mg5_aMC'
--- bin/mg5_aMC	2019-05-05 18:53:50 +0000
+++ bin/mg5_aMC	2020-07-04 19:03:17 +0000
@@ -49,7 +49,8 @@
                  help='force to launch debug mode')
 parser.add_option("-m", "--mode", dest="plugin",
                   help="Define some additional command provide by a PLUGIN")
-                               
+parser.add_option("-s", "--nocaffeinate", action="store_false", default=True, dest='nosleep',
+                  help='For mac user, forbids to use caffeinate when running with a script')                               
 (options, args) = parser.parse_args()
 if len(args) == 0:
     args = ''
@@ -66,6 +67,9 @@
 import logging.config
 import madgraph.interface.coloring_logging
 
+if ' ' in os.getcwd():
+    logging.warning("Path does contains spaces. We advise that you change your current path to avoid to have space in the path.")
+
 try: 
     import readline
 except ImportError:
@@ -157,6 +161,10 @@
 # Call the cmd interface main loop
 try:
     if options.file or args:
+        if sys.platform == "darwin" and options.nosleep:
+            logging.getLogger('madgraph').warning("launching caffeinate to prevent idle sleep when MG5aMC is running. Run './bin/mg5_aMC -s' to prevent this.")
+            pid = os.getpid()
+            subprocess.Popen(['caffeinate', '-i', '-w', str(pid)])
         # They are an input file 
         if args:
             input_file = os.path.realpath(args[0])

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2020-03-11 09:28:14 +0000
+++ madgraph/core/base_objects.py	2020-07-04 19:03:17 +0000
@@ -3546,18 +3546,31 @@
         final.sort()
         return (tuple(initial), tuple(final))
     
-    def get_final_ids_after_decay(self):
+    def get_initial_final_ids_after_decay(self, max_depth=-1):
+        """return a tuple of two tuple containing the id of the initial/final
+           state particles. Each list is ordered"""
+           
+        initial = [l.get('id') for l in self.get('legs')\
+              if not l.get('state')]
+        final = self.get_final_ids_after_decay(max_depth=max_depth)
+        initial.sort()
+        final.sort()
+        return (tuple(initial), tuple(final))
+        
+    
+    def get_final_ids_after_decay(self, max_depth=-1):
         """Give the pdg code of the process including decay"""
         
         finals = self.get_final_ids()
-        for proc in self.get('decay_chains'):
-            init = proc.get_initial_ids()[0]
-            #while 1:
-            try:
-                pos = finals.index(init)
-            except:
-                break
-            finals[pos] = proc.get_final_ids_after_decay()
+        if max_depth !=0 :
+            for proc in self.get('decay_chains'):
+                init = proc.get_initial_ids()[0]
+                #while 1:
+                try:
+                    pos = finals.index(init)
+                except:
+                    break
+                finals[pos] = proc.get_final_ids_after_decay(max_depth-1)
         output = []
         for d in finals:
             if isinstance(d, list):

=== modified file 'madgraph/core/diagram_generation.py'
--- madgraph/core/diagram_generation.py	2020-03-11 09:28:14 +0000
+++ madgraph/core/diagram_generation.py	2020-07-04 19:03:17 +0000
@@ -896,7 +896,7 @@
                                              fcts=['remove_diag'])
         else:
             #example and simple tests
-            def remove_diag(diag):
+            def remove_diag(diag, model=None):
                 for vertex in diag['vertices']: #last 
                     if vertex['id'] == 0: #special final vertex
                         continue 
@@ -907,8 +907,9 @@
 
         res = diag_list.__class__()                
         nb_removed = 0 
+        model = self['process']['model'] 
         for diag in diag_list:
-            if remove_diag(diag):
+            if remove_diag(diag, model):
                 nb_removed +=1
             else:
                 res.append(diag)

=== modified file 'madgraph/core/helas_objects.py'
--- madgraph/core/helas_objects.py	2020-03-09 09:27:32 +0000
+++ madgraph/core/helas_objects.py	2020-07-04 19:03:17 +0000
@@ -5871,6 +5871,7 @@
                     # Identical matrix element found
                     other_processes = identified_matrix_elements[me_index].\
                                       get('processes')
+                    
                     other_processes.append(cls.reorder_process(\
                         amplitude.get('process'),
                         permutations[me_index],
@@ -5911,14 +5912,24 @@
         """Reorder the legs in the process according to the difference
         between org_perm and proc_perm"""
 
+        
+        
         leglist = base_objects.LegList(\
                   [copy.copy(process.get('legs_with_decays')[i]) for i in \
                    diagram_generation.DiagramTag.reorder_permutation(\
                        proc_perm, org_perm)])
         new_proc = copy.copy(process)
+        if org_perm == proc_perm:
+            return new_proc
+
+        if  len(org_perm) != len(process.get('legs_with_decays')):
+            raise Exception, "issue on symmetry between process"
+
         new_proc.set('legs_with_decays', leglist)
-
+        
         if not new_proc.get('decay_chains'):
             new_proc.set('legs', leglist)
+            assert len(process.get('legs')) == len(leglist)
 
+            
         return new_proc

=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
--- madgraph/interface/amcatnlo_run_interface.py	2020-03-11 09:28:14 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2020-07-04 19:03:17 +0000
@@ -174,7 +174,13 @@
             logger.warning(msg % compiler)
     else:
         curr_version = misc.get_gfortran_version(compiler)
-        if not ''.join(curr_version.split('.')) >= '46':
+        curr_version = curr_version.split('.')
+        if len(curr_version) == 1:
+            curr_version.append(0)
+
+        if int(curr_version[0]) < 5:
+            if int(curr_version[0]) == 4 and int(curr_version[1]) > 5:
+                return
             if block:
                 raise aMCatNLOError(msg % (compiler + ' ' + curr_version))
             else:
@@ -2933,7 +2939,7 @@
         """Sums all the plots in the HwU format."""
         logger.debug('Combining HwU plots.')
 
-        command =  []
+        command =  [sys.executable]
         command.append(pjoin(self.me_dir, 'bin', 'internal','histograms.py'))
         for job in jobs:
             if job['dirname'].endswith('.HwU'):

=== modified file 'madgraph/interface/common_run_interface.py'
--- madgraph/interface/common_run_interface.py	2020-03-11 09:28:14 +0000
+++ madgraph/interface/common_run_interface.py	2020-07-04 19:03:17 +0000
@@ -1039,7 +1039,11 @@
             out = ask(question, '0', possible_answer, timeout=int(1.5*timeout),
                               path_msg='enter path', ask_class = AskforEditCard,
                               cards=cards, mode=mode, **opt)
-
+            if 'return_instance' in opt and opt['return_instance']:
+                out, cmd = out
+        if 'return_instance' in opt and opt['return_instance']:
+            return (out, cmd)
+        return out
 
     @staticmethod
     def detect_card_type(path):
@@ -2839,8 +2843,8 @@
                     reco_output = pjoin(self.me_dir,
                            'MA5_%s_ANALYSIS%s_%d'%(mode.upper(),MA5_runtag,i+1))
                     # Look for either a root or .lhe.gz output
-                    reco_event_file = misc.glob('*.lhe.gz',pjoin(reco_output,'Output','_reco_events','lheEvents0_%d'%MA5_run_number))+\
-                                       misc.glob('*.root',pjoin(reco_output,'Output','_reco_events', 'RecoEvents0_%d'%MA5_run_number))
+                    reco_event_file = misc.glob('*.lhe.gz',pjoin(reco_output,'Output','SAF','_reco_events','lheEvents0_%d'%MA5_run_number))+\
+                                       misc.glob('*.root',pjoin(reco_output,'Output','SAF','_reco_events', 'RecoEvents0_%d'%MA5_run_number))
                     if len(reco_event_file)==0:
                         raise MadGraph5Error, "MadAnalysis5 failed to produce the "+\
                   "reconstructed event file for reconstruction '%s'."%MA5_runtag[6:]
@@ -2855,7 +2859,7 @@
                     parent_dir_name = os.path.basename(os.path.dirname(reco_event_file))
                     files.ln(pjoin(self.me_dir,'HTML',self.run_name,
                       '%s_MA5_%s_ANALYSIS%s_%d'%(self.run_tag,mode.upper(),
-                      MA5_runtag,i+1),'Output','_reco_events',parent_dir_name,links_created[-1]),
+                      MA5_runtag,i+1),'Output','SAF','_reco_events',parent_dir_name,links_created[-1]),
                                       pjoin(self.me_dir,'Events',self.run_name))
 
                 logger.info("MadAnalysis5 successfully completed the reconstruction "+
@@ -3331,7 +3335,7 @@
             if mass and abs(width/mass) < 1e-12:
                 if hasattr(interface, 'run_card') and isinstance(interface.run_card, banner_mod.RunCardLO):
                     if interface.run_card['small_width_treatment'] < 1e-12:
-                        logger.error('The width of particle %s is too small for an s-channel resonance (%s) and the small_width_paramer is too small to prevent numerical issues. If you have this particle in an s-channel, this is likely to create numerical instabilities .', param.lhacode[0], width)
+                        logger.error('The width of particle %s is too small for an s-channel resonance (%s) and the small_width_treatment parameter is too small to prevent numerical issues. If you have this particle in an s-channel, this is likely to create numerical instabilities .', param.lhacode[0], width)
                 else:
                     logger.error('The width of particle %s is too small for an s-channel resonance (%s). If you have this particle in an s-channel, this is likely to create numerical instabilities .', param.lhacode[0], width)
                 if CommonRunCmd.sleep_for_error:
@@ -4102,6 +4106,18 @@
             if not require_local and (os.path.exists(pjoin(pdfsets_dir, pdfset)) or \
                                     os.path.isdir(pjoin(pdfsets_dir, pdfset))):
                 continue
+            if not require_local:
+                if 'LHAPDF_DATA_PATH' in os.environ:
+                    found = False
+                    for path in os.environ['LHAPDF_DATA_PATH'].split(":"):
+                        if (os.path.exists(pjoin(path, pdfset)) or \
+                                    os.path.isdir(pjoin(path, pdfset))):
+                            found =True
+                            break
+                    if found:
+                        continue
+                    
+                    
             #check that the pdfset is not already there
             elif not os.path.exists(pjoin(self.me_dir, 'lib', 'PDFsets', pdfset)) and \
                not os.path.isdir(pjoin(self.me_dir, 'lib', 'PDFsets', pdfset)):
@@ -4313,7 +4329,7 @@
                         stdout = subprocess.PIPE).stdout.read().strip()
         except OSError, error:
             if error.errno == 2:
-                raise Exception, 'lhapdf executable (%s) is not found on your system. Please install it and/or indicate the path to the correct executable in input/mg5_configuration.txt' % self.options['lhapdf']
+                raise Exception, 'lhapdf executable (%s) is not found on your system. Please install it and/or indicate the path to the correct executable in input/mg5_configuration.txt' % lhapdf_config
             else:
                 raise
                 
@@ -4499,7 +4515,12 @@
         self.modified_card = set() #set of cards not in sync with filesystem
                               # need to sync them before editing/leaving
         self.init_from_banner(from_banner, banner)
-        
+        self.writting_card = True
+        if 'write_file' in opt:
+            if not opt['write_file']:
+                self.writting_card = False
+                self.param_consistency = False
+                        
         #update default path by custom one if specify in cards
         for card in cards:
             if os.path.exists(card) and os.path.sep in cards:
@@ -4777,10 +4798,12 @@
             return []
         
         self.special_shortcut.update({
-            'spinmode':([str], ['add madspin_card --before_line="launch" set spinmode %(0)s'])
+            'spinmode':([str], ['add madspin_card --before_line="launch" set spinmode %(0)s']),
+            'nodecay':([], ['edit madspin_card --comment_line="decay"'])
             })
         self.special_shortcut_help.update({
-            'spinmode' : 'full|none|onshell. Choose the mode of madspin.\n   - full: spin-correlation and off-shell effect\n  - onshell: only spin-correlation,]\n  - none: no spin-correlation and not offshell effects.'
+            'spinmode' : 'full|none|onshell. Choose the mode of madspin.\n   - full: spin-correlation and off-shell effect\n  - onshell: only spin-correlation,]\n  - none: no spin-correlation and not offshell effects.',
+            'nodecay': 'remove all decay previously defined in madspin',
              })
         return []
     
@@ -5245,7 +5268,7 @@
     def do_set(self, line):
         """ edit the value of one parameter in the card"""
         
-        
+
         args = self.split_arg(line)
         
         
@@ -5508,7 +5531,7 @@
                 try:
                     key = tuple([int(i) for i in args[start+1:-1]])
                 except ValueError:
-                    if args[start] == 'decay' and args[start+1:-1] == ['all']:
+                    if args[start+1:-1] == ['all']:
                         for key in self.param_card[args[start]].param_dict:
                             if (args[start], key) in self.restricted_value:
                                 continue
@@ -5969,8 +5992,9 @@
                     self.do_set('shower_card extrapaths None ') 
                     
         # ensure that all cards are in sync
-        for key in list(self.modified_card):
-            self.write_card(key)
+        if self.writting_card:
+            for key in list(self.modified_card):
+                self.write_card(key)
 
 
     def reask(self, *args, **opt):
@@ -6096,7 +6120,8 @@
         self.run_card.write(self.paths['run'], self.paths['run_default'])
         
     def write_card_param(self):
-        """ write the param_card """        
+        """ write the param_card """    
+    
         self.param_card.write(self.paths['param'])
         
     @staticmethod
@@ -6481,6 +6506,7 @@
         logger.info( '     --before_line="<regular-expression>" write the line before the first line matching the regular expression')
         logger.info( '     --replace_line="<regular-expression>" replace the line matching the regular expression')
         logger.info( '     --clean remove all previously existing line in  the file')
+        logger.info( '     --comment_line="<regular-expression>"  comment all lines matching the regular expression')
         logger.info('')
         logger.info('    Note: all regular-expression will be prefixed by ^\s*')
         logger.info('')
@@ -6626,7 +6652,29 @@
                 logger.info("Replacing the line \"%s\" [line %d of %s] by \"%s\"" %
                          (old_line, posline, card, new_line ),'$MG:BOLD') 
                 self.last_editline_pos = posline               
-                                            
+
+            elif args[1].startswith('--comment_line='):
+                # catch the line/regular expression and replace the associate line
+                # if no line match go to check if args[2] has other instruction starting with --
+                text = open(path).read()
+                split = text.split('\n')
+                search_pattern=r'''comment_line=(?P<quote>["'])(?:(?=(\\?))\2.)*?\1'''
+                pattern = '^\s*' + re.search(search_pattern, line).group()[14:-1]
+                nb_mod = 0
+                for posline,l in enumerate(split):
+                    if re.search(pattern, l):
+                        split[posline] = '#%s' % l
+                        nb_mod +=1
+                        logger.info("Commenting line \"%s\" [line %d of %s]" %
+                         (l, posline, card ),'$MG:BOLD') 
+                        # overwrite the previous line
+                if not nb_mod:
+                    logger.warning('no line commented (no line matching)')
+                ff = open(path,'w')
+                ff.write('\n'.join(split))
+
+                self.last_editline_pos = posline               
+
             
             elif args[1].startswith('--before_line='):
                 # catch the line/regular expression and write before that line

=== modified file 'madgraph/interface/extended_cmd.py'
--- madgraph/interface/extended_cmd.py	2020-03-11 09:28:14 +0000
+++ madgraph/interface/extended_cmd.py	2020-07-04 19:03:17 +0000
@@ -2594,6 +2594,8 @@
         elif line in 'auto':
             self.switch['dynamical'] = True
             return super(ControlSwitch, self).default(line)
+        elif line.startswith('set ') and not hasattr(self.__class__, 'do_set'):
+            raise NotValidInput('unknow command: %s. Did you mean \"%s\"' % (line, line[4:]))
         elif raise_error:
             raise NotValidInput('unknow command: %s' % line)
         else:
@@ -2609,6 +2611,8 @@
             getattr(self, 'ans_%s' % base)(value)
         elif base in self.switch:
             self.set_switch(base, value)
+        elif line.startswith('set ') and not hasattr(self.__class__, 'do_set'):
+            raise NotValidInput('Not valid command: %s. Did you mean \"%s\"' % (line, line[4:]))
         elif raise_error:
             raise NotValidInput('Not valid command: %s' % line)                
         else:

=== modified file 'madgraph/interface/madevent_interface.py'
--- madgraph/interface/madevent_interface.py	2020-03-16 08:37:10 +0000
+++ madgraph/interface/madevent_interface.py	2020-07-04 19:03:17 +0000
@@ -797,7 +797,7 @@
         
         self.allowed_madspin = []
         if 'MadSpin'  in self.available_module:
-            self.allowed_madspin = ['OFF',"ON",'onshell']
+            self.allowed_madspin = ['OFF',"ON",'onshell',"full"]
         return self.allowed_madspin
     
     def check_value_madspin(self, value):
@@ -834,7 +834,7 @@
         if value == 'onshell':
             return ["edit madspin_card --replace_line='set spinmode' --before_line='decay' set spinmode onshell"]
         elif value in ['full', 'madspin']:
-            return ["edit madspin_card --replace_line='set spinmode' --before_line='decay' set spinmode madspin"]
+            return ["edit madspin_card --replace_line='set spinmode' --before_line='decay' set spinmode full"]
         elif value == 'none':
             return ["edit madspin_card --replace_line='set spinmode' --before_line='decay' set spinmode none"]
         else:
@@ -1291,7 +1291,6 @@
                 raise self.InvalidCmd('No run_name currently define. Unable to run refine')
 
         if len(args) > 2:
-            self.help_refine()
             raise self.InvalidCmd('Too many argument for refine command')
         else:
             try:
@@ -2497,7 +2496,8 @@
                           postcmd=False)
             self.exec_cmd('combine_events', postcmd=False)
             self.exec_cmd('store_events', postcmd=False)
-            self.exec_cmd('decay_events -from_cards', postcmd=False)
+            with misc.TMP_variable(self, 'run_name', self.run_name):
+                self.exec_cmd('decay_events -from_cards', postcmd=False)
             self.exec_cmd('create_gridpack', postcmd=False)
         else:
             # Regular run mode
@@ -2526,7 +2526,8 @@
             
             #we can bypass the following if scan and first result is zero
             if not bypass_run:
-                self.exec_cmd('refine %s' % nb_event, postcmd=False)
+                self.exec_cmd('refine %s --treshold=%s' % (nb_event,self.run_card['second_refine_treshold'])
+                              , postcmd=False)
             
                 self.exec_cmd('combine_events', postcmd=False,printcmd=False)
                 self.print_results_in_shell(self.results.current)
@@ -2882,7 +2883,7 @@
                         particle = 0
                 # Read BRs for this decay
                 line = param_card[line_number]
-                while line.startswith('#') or line.startswith(' '):
+                while re.search('^(#|\s|\d)', line):
                     line = param_card.pop(line_number)
                     if not particle or line.startswith('#'):
                         line=param_card[line_number]
@@ -3354,7 +3355,7 @@
         if float(self.run_card['mmjj']) > 0.01 * (float(self.run_card['ebeam1'])+float(self.run_card['ebeam2'])):
             self.pass_in_difficult_integration_mode()
         elif self.run_card['hard_survey']:
-            self.pass_in_difficult_integration_mode()
+            self.pass_in_difficult_integration_mode(self.run_card['hard_survey'])
             
         jobs, P_zero_result = ajobcreator.launch()
         # Check if all or only some fails
@@ -3382,24 +3383,26 @@
         self.update_status('End survey', 'parton', makehtml=False)
 
     ############################################################################
-    def pass_in_difficult_integration_mode(self):
+    def pass_in_difficult_integration_mode(self, rate=1):
         """be more secure for the integration to not miss it due to strong cut"""
         
         # improve survey options if default
         if self.opts['points'] == self._survey_options['points'][1]:
-            self.opts['points'] = 3 * self._survey_options['points'][1]
+            self.opts['points'] = (rate+2) * self._survey_options['points'][1]
         if self.opts['iterations'] == self._survey_options['iterations'][1]:
-            self.opts['iterations'] = 2 + self._survey_options['iterations'][1]
+            self.opts['iterations'] = 1 + rate + self._survey_options['iterations'][1]
         if self.opts['accuracy'] == self._survey_options['accuracy'][1]:
-            self.opts['accuracy'] = self._survey_options['accuracy'][1]/3  
+            self.opts['accuracy'] = self._survey_options['accuracy'][1]/(rate+2)  
         
         # Modify run_config.inc in order to improve the refine
         conf_path = pjoin(self.me_dir, 'Source','run_config.inc')
         files.cp(conf_path, conf_path + '.bk')
         #
         text = open(conf_path).read()
-        text = re.sub('''\(min_events = \d+\)''', '''(min_events = 7500 )''', text)
-        text = re.sub('''\(max_events = \d+\)''', '''(max_events = 40000 )''', text)
+        min_evt, max_evt = 2500 *(2+rate), 10000*(rate+1) 
+        
+        text = re.sub('''\(min_events = \d+\)''', '(min_events = %i )' % min_evt, text)
+        text = re.sub('''\(max_events = \d+\)''', '(max_events = %i )' % max_evt, text)
         fsock = open(conf_path, 'w')
         fsock.write(text)
         fsock.close()
@@ -3415,6 +3418,18 @@
         devnull = open(os.devnull, 'w')  
         self.nb_refine += 1
         args = self.split_arg(line)
+        treshold=None
+        for a in args:
+            if a.startswith('--treshold='):
+                treshold = float(a.split('=',1)[1])
+                old_xsec = self.results.current['prev_cross']
+                new_xsec = self.results.current['cross']
+                if old_xsec > new_xsec * treshold:
+                    logger.info('No need for second refine due to stability of cross-section')
+                    return
+                else:
+                    args.remove(a)
+                    break
         # Check argument's validity
         self.check_refine(args)
         
@@ -3471,8 +3486,7 @@
             cross, error = x_improve.update_html() #update html results for survey
             if  cross == 0:
                 return
-            logger.info("Current estimate of cross-section: %s +- %s" % (cross, error))
-        
+            logger.info("- Current estimate of cross-section: %s +- %s" % (cross, error))
         if isinstance(x_improve, gen_ximprove.gen_ximprove_v4):
             # Non splitted mode is based on writting ajob so need to track them
             # Splitted mode handle the cluster submition internally.
@@ -3811,6 +3825,17 @@
         self.update_status('Creating gridpack', level='parton')
         # compile gen_ximprove
         misc.compile(['../bin/internal/gen_ximprove'], cwd=pjoin(self.me_dir, "Source"))
+        
+        Gdir = self.get_Gdir()
+        Pdir = set([os.path.dirname(G) for G in Gdir])
+        for P in Pdir: 
+            allG = misc.glob('G*', path=P)
+            for G in allG:
+                if pjoin(P, G) not in Gdir:
+                    logger.debug('removing %s', pjoin(P,G))
+                    shutil.rmtree(pjoin(P,G))
+                    
+        
         args = self.split_arg(line)
         self.check_combine_events(args)
         if not self.run_tag: self.run_tag = 'tag_1'
@@ -4668,7 +4693,7 @@
                             devnull.close()
                             if pid == 0:
                                 misc.call('head -n -1 %s | tail -n +%d > %s/tmpfile' %
-                                          (hepmc_file, n_head, os.path.dirname(hepmc_file)), shell=True)
+                                          (hepmc_file, n_head+1, os.path.dirname(hepmc_file)), shell=True)
                                 misc.call(['mv', 'tmpfile', os.path.basename(hepmc_file)], cwd=os.path.dirname(hepmc_file))
                             elif sys.platform == 'darwin':
                                 # sed on MAC has slightly different synthax than on
@@ -5967,10 +5992,16 @@
 
         eradir = self.options['exrootanalysis_path']
         totar = False
+        torm = False
         if input.endswith('.gz'):
-            misc.gunzip(input, keep=True)
-            totar = True
-            input = input[:-3]
+            if not os.path.exists(input) and os.path.exists(input[:-3]):
+                totar = True
+                input = input[:-3]
+            else:
+                misc.gunzip(input, keep=True)
+                totar = False
+                torm = True
+                input = input[:-3]
             
         try:
             misc.call(['%s/ExRootLHEFConverter' % eradir, 
@@ -5982,12 +6013,13 @@
         if totar:
             if os.path.exists('%s.gz' % input):
                 try:
-                    os.remove(input)
+                    os.remove('%s.gz' % input)
                 except:
                     pass
             else:
                 misc.gzip(input)
-            
+        if torm:
+            os.remove(input)
     
     def run_syscalc(self, mode='parton', event_path=None, output=None):
         """create the syscalc output""" 
@@ -6458,7 +6490,10 @@
             os.chdir(self.me_dir)
         else:
             for line in open(pjoin(self.me_dir,'SubProcesses','subproc.mg')):
-                os.mkdir(line.strip())
+                p = line.strip()
+                os.mkdir(p)
+                files.cp(pjoin(self.me_dir,'SubProcesses',p,'symfact.dat'),
+                         pjoin(p, 'symfact.dat'))
             
 
     def launch(self, nb_event, seed):
@@ -6472,6 +6507,12 @@
             misc.call([pjoin(self.me_dir,'bin','internal','restore_data'),
                          'default'], cwd=self.me_dir)
 
+        if self.run_card['python_seed'] == -2:
+            import random
+            random.seed(seed)
+        elif self.run_card['python_seed'] > 0:
+            import random
+            random.seed(self.run_card['python_seed'])            
         # 2) Run the refine for the grid
         self.update_status('Generating Events', level=None)
         #misc.call([pjoin(self.me_dir,'bin','refine4grid'),
@@ -6649,11 +6690,11 @@
                 sum_axsec += result.get('axsec')*gscalefact[Gdir]
                 
                 if len(AllEvent) >= 80: #perform a partial unweighting
-                    AllEvent.unweight(pjoin(self.me_dir, "Events", self.run_name, "partials%s.lhe.gz" % partials),
+                    AllEvent.unweight(pjoin(outdir, self.run_name, "partials%s.lhe.gz" % partials),
                           get_wgt, log_level=5,  trunc_error=1e-2, event_target=self.nb_event)
                     AllEvent = lhe_parser.MultiEventFile()
                     AllEvent.banner = self.banner
-                    AllEvent.add(pjoin(self.me_dir, "Events", self.run_name, "partials%s.lhe.gz" % partials),
+                    AllEvent.add(pjoin(outdir, self.run_name, "partials%s.lhe.gz" % partials),
                                  sum_xsec,
                                  math.sqrt(sum(x**2 for x in sum_xerru)),
                                  sum_axsec) 

=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2020-03-11 09:28:14 +0000
+++ madgraph/interface/madgraph_interface.py	2020-07-04 19:03:17 +0000
@@ -459,7 +459,7 @@
         logger.info("      -d: specify other MG/ME directory")
         logger.info("      -noclean: no cleaning performed in \"path\".")
         logger.info("      -nojpeg: no jpeg diagrams will be generated.")
-        logger.info("      -noeps: no jpeg and eps diagrams will be generated.")
+        logger.info("      --noeps=True: no jpeg and eps diagrams will be generated.")
         logger.info("      -name: the postfix of the main file in pythia8 mode.")
         logger.info("   Examples:",'$MG:color:GREEN')
         logger.info("       output",'$MG:color:GREEN')
@@ -530,6 +530,11 @@
         logger.info(" > Except for the 'gauge' test, all checks above are also")
         logger.info("   available for loop processes with ML5 ('virt=' mode)")
         logger.info("Example: check full p p > j j",'$MG:color:GREEN')
+        logger.info("Using leshouches file as input",'$MG:color:GREEN')
+        logger.info("    use the option --events=PATH")
+        logger.info("      zipped file are not supported")
+        logger.info("      to loop over the file use the option --skip_evt=X")
+        logger.info("")
         logger.info("Options for loop processes only:",'$MG:BOLD')
         logger.info("o timing:",'$MG:color:GREEN')
         logger.info("   Generate and output a process and returns detailed")
@@ -928,7 +933,9 @@
                    '--helicity':'-1','--seed':'-1','--collier_cache':'-1',
                    '--collier_req_acc':'auto',
                    '--collier_internal_stability_test':'False',
-                   '--collier_mode':'1'}  
+                   '--collier_mode':'1',
+                   '--events': None,
+                   '--skip_evt':0}  
 
         if args[0] in ['cms'] or args[0].lower()=='cmsoptions':
             # increase the default energy to 5000
@@ -1077,7 +1084,7 @@
             
         if '[' in process and '{' in process:
             valid = False
-            if 'noborn' in process:
+            if 'noborn' in process or 'sqrvirt' in process:
                 valid = True
             else:
                 raise self.InvalidCmd('Polarization restriction can not be used for NLO processes')
@@ -1673,7 +1680,7 @@
                     continue
                 elif not '=' in arg:
                     raise self.InvalidCmd('Options required an equal (and then the value)')
-                arg, value = arg.split('=')
+                arg, value = arg.split('=',1)
                 if arg[2:] not in options:
                     raise self.InvalidCmd('%s not valid options' % arg)
                 options[arg[2:]] = value
@@ -2836,7 +2843,7 @@
     _import_formats = ['model_v4', 'model', 'proc_v4', 'command', 'banner']
     _install_opts = ['Delphes', 'MadAnalysis4', 'ExRootAnalysis',
                      'update', 'Golem95', 'QCDLoop', 'maddm', 'maddump',
-                     'looptools']
+                     'looptools', 'MadSTR']
     
     # The targets below are installed using the HEPToolsInstaller.py script
     _advanced_install_opts = ['pythia8','zlib','boost','lhapdf6','lhapdf5','collier',
@@ -3876,6 +3883,15 @@
             option = args[i].split('=')
             if option[0] =='--energy':
                 options['energy']=float(option[1])
+            elif option[0] == '--events' and option[1]:
+                if option[1] == 'None':
+                    options['events'] = None
+                elif not os.path.exists(option[1]):
+                    raise Exception, 'path %s does not exists' % option[1]
+                else:
+                    options['events'] = option[1]
+            elif option[0] == '--skip_evt':
+                options['skip_evt']=int(option[1])
             elif option[0]=='--split_orders':
                 options['split_orders']=int(option[1])
             elif option[0]=='--helicity':
@@ -5622,7 +5638,7 @@
                     #self.do_define(line)
                     self.exec_cmd('define %s' % line, printcmd=False, precmd=True)
             except self.InvalidCmd, why:
-                logger_stderr.warning('impossible to set default multiparticles %s because %s' %
+                logger.warning('impossible to set default multiparticles %s because %s' %
                                         (line.split()[0],why))
                 if self.history[-1] == 'define %s' % line.strip():
                     self.history.pop(-1)
@@ -6013,7 +6029,7 @@
          # Return true for successful installation
         return True
 
-    install_plugin = ['maddm', 'maddump']
+    install_plugin = ['maddm', 'maddump', 'MadSTR']
     install_ad = {'pythia-pgs':['arXiv:0603175'],
                           'Delphes':['arXiv:1307.6346'],
                           'Delphes2':['arXiv:0903.2225'],
@@ -6031,7 +6047,8 @@
                           'collier':['arXiv:1604.06792'],
                           'oneloop':['arXiv:1007.4716'],
                           'maddm':['arXiv:1804.00444'],
-                          'maddump':['arXiv:1812.06771']}
+                          'maddump':['arXiv:1812.06771'],
+                          'MadSTR':['arXiv:1612.00440']}
     
     install_server = ['http://madgraph.phys.ucl.ac.be/package_info.dat',
                          'http://madgraph.physics.illinois.edu/package_info.dat']
@@ -6157,7 +6174,10 @@
             name = args[0]
         if args[0] == 'MadAnalysis4':
             args[0] = 'MadAnalysis'
-        
+        elif args[0] in ['madstr', 'madSTR']:
+            args[0] = 'MadSTR'
+            name = 'MadSTR'
+            
         if args[0] in self._advanced_install_opts:
             # Now launch the advanced installation of the tool args[0]
             # path['HEPToolsInstaller'] is the online adress where to downlaod
@@ -6380,6 +6400,7 @@
                                                cwd = os.path.join(MG5DIR, name))
             else:
                 status = misc.call(['make']+make_flags, cwd = os.path.join(MG5DIR, name))
+            devnull.close()
         else:
             try:
                 misc.compile(['clean'], mode='', cwd = os.path.join(MG5DIR, name))
@@ -6428,18 +6449,18 @@
 
             if sys.platform == "darwin":
                 logger.info('Downloading TD for Mac')
-                target = 'http://madgraph.phys.ucl.ac.be/Downloads/td_mac_intel.tar.gz'
+                target = 'https://home.fnal.gov/~parke/TD/td_mac_intel64.tar.gz'
                 misc.wget(target, 'td.tgz', cwd=pjoin(MG5DIR,'td'))
                 misc.call(['tar', '-xzpvf', 'td.tgz'],
                                                   cwd=pjoin(MG5DIR,'td'))
-                files.mv(MG5DIR + '/td/td_mac_intel',MG5DIR+'/td/td')
+                files.mv(MG5DIR + '/td/td_intel_mac64',MG5DIR+'/td/td')
             else:
                 if sys.maxsize > 2**32:
                     logger.info('Downloading TD for Linux 64 bit')
-                    target = 'http://madgraph.phys.ucl.ac.be/Downloads/td64/td'
-                    logger.warning('''td program (needed by MadAnalysis) is not compile for 64 bit computer.
-                In 99% of the case, this is perfectly fine. If you do not have plot, please follow 
-                instruction in https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/TopDrawer .''')
+                    target = 'https://home.fnal.gov/~parke/TD/td_linux_64bit.tar.gz'
+                    #logger.warning('''td program (needed by MadAnalysis) is not compile for 64 bit computer.
+                #In 99% of the case, this is perfectly fine. If you do not have plot, please follow 
+                #instruction in https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/TopDrawer .''')
                 else:                    
                     logger.info('Downloading TD for Linux 32 bit')
                     target = 'http://madgraph.phys.ucl.ac.be/Downloads/td'

=== modified file 'madgraph/interface/master_interface.py'
--- madgraph/interface/master_interface.py	2019-09-26 14:06:19 +0000
+++ madgraph/interface/master_interface.py	2020-07-04 19:03:17 +0000
@@ -215,6 +215,8 @@
                 elif nlo_mode in ['virt', 'sqrvirt']:
                     self.change_principal_cmd('MadLoop', allow_switch)
                 elif nlo_mode == 'noborn': 
+                    if self.current_interface == "MadGraph":
+                        allow_switch = True
                     self.change_principal_cmd('MadLoop', allow_switch)
                     self.cmd.validate_model(self, loop_type=nlo_mode,
                                                             coupling_type=orders)

=== modified file 'madgraph/interface/reweight_interface.py'
--- madgraph/interface/reweight_interface.py	2020-03-11 09:28:14 +0000
+++ madgraph/interface/reweight_interface.py	2020-07-04 19:03:17 +0000
@@ -58,7 +58,7 @@
 nb_f2py_module = 0 # each time the process/model is changed this number is modified to 
                    # forced the python module to re-create an executable
 
-lhapdf = None
+#lhapdf = None
 
 
 class ReweightInterface(extended_cmd.Cmd):
@@ -688,7 +688,8 @@
             for i,card in enumerate(param_card_iterator):
                 if self.options['rwgt_name']:
                     self.options['rwgt_name'] = '%s_%s' % (self.options['rwgt_name'].rsplit('_',1)[0], i+1)
-                card.write(pjoin(rw_dir, 'Cards', 'param_card.dat'))
+                self.new_param_card = card
+                #card.write(pjoin(rw_dir, 'Cards', 'param_card.dat'))
                 self.exec_cmd("launch --keep_card", printcmd=False, precmd=True)
         
         self.options['rwgt_name'] = None
@@ -708,24 +709,23 @@
         
         
         if not '--keep_card' in args:
-            ff = open(pjoin(rw_dir,'Cards', 'param_card.dat'), 'w')
-            ff.write(self.banner['slha'])
-            ff.close()
             if self.has_nlo and self.rwgt_mode != "LO":
                 rwdir_virt = rw_dir.replace('rw_me', 'rw_mevirt')
-                files.ln(ff.name, starting_dir=pjoin(rwdir_virt, 'Cards')) 
-            ff = open(pjoin(path_me, 'rw_me','Cards', 'param_card_orig.dat'), 'w')
-            ff.write(self.banner['slha'])
-            ff.close()      
-            if self.has_nlo and self.rwgt_mode != "LO":
-                files.ln(ff.name, starting_dir=pjoin(path_me, 'rw_mevirt', 'Cards'))
-            cmd = common_run_interface.CommonRunCmd.ask_edit_card_static(cards=['param_card.dat'],
-                                   ask=self.ask, pwd=rw_dir, first_cmd=self.stored_line)
+                
+            out, cmd = common_run_interface.CommonRunCmd.ask_edit_card_static(cards=['param_card.dat'],
+                                   ask=self.ask, pwd=rw_dir, first_cmd=self.stored_line,
+                                   write_file=False, return_instance=True
+                                   )
             self.stored_line = None
-        
+            card = cmd.param_card
+            new_card = card.write()
+        elif self.new_param_card:
+            new_card = self.new_param_card.write()
+        else:
+            new_card = open(pjoin(rw_dir, 'Cards', 'param_card.dat')).read()
+            
         # check for potential scan in the new card 
-        new_card = open(pjoin(rw_dir, 'Cards', 'param_card.dat')).read()
-        pattern_scan = re.compile(r'''^[\s\d]*scan''', re.I+re.M) 
+        pattern_scan = re.compile(r'''^(decay)?[\s\d]*scan''', re.I+re.M) 
         param_card_iterator = []
         if pattern_scan.search(new_card):
             import madgraph.interface.extended_cmd as extended_cmd
@@ -744,11 +744,19 @@
             param_card_iterator = main_card
             first_card = param_card_iterator.next(autostart=True)
             new_card = first_card.write()
-            first_card.write(pjoin(rw_dir, 'Cards', 'param_card.dat'))  
-                          
-        # check if "Auto" is present for a width parameter
+            self.new_param_card = first_card
+            #first_card.write(pjoin(rw_dir, 'Cards', 'param_card.dat'))  
+        
+        # check if "Auto" is present for a width parameter)
         tmp_card = new_card.lower().split('block',1)[1]
         if "auto" in tmp_card: 
+            if param_card_iterator:
+                first_card.write(pjoin(rw_dir, 'Cards', 'param_card.dat'))
+            else:
+                ff = open(pjoin(rw_dir, 'Cards', 'param_card.dat'),'w')
+                ff.write(new_card)
+                ff.close()
+                
             self.mother.check_param_card(pjoin(rw_dir, 'Cards', 'param_card.dat'))
             new_card = open(pjoin(rw_dir, 'Cards', 'param_card.dat')).read()
 
@@ -795,6 +803,7 @@
         # add the reweighting in the banner information:
         #starts by computing the difference in the cards.
         s_orig = self.banner['slha']
+        self.orig_param_card_text = s_orig
         s_new = new_card
         self.new_param_card = check_param_card.ParamCard(s_new.splitlines())
         
@@ -866,15 +875,25 @@
         else:
             tag_name = 'rwgt_%s' % rewgtid
 
+                
         #initialise module.
         for (path,tag), module in self.f2pylib.items():
             with misc.chdir(pjoin(os.path.dirname(rw_dir), path)):
                 with misc.stdchannel_redirected(sys.stdout, os.devnull):
                     if 'second' in path or tag == 3:
-                        module.initialise(pjoin(rw_dir, 'Cards', 'param_card.dat'))
+                        param_card = self.new_param_card
                     else:
-                        module.initialise(pjoin(path_me, 'rw_me', 'Cards', 'param_card_orig.dat'))
+                        param_card = check_param_card.ParamCard(self.orig_param_card_text)
+                    
+                    for block in param_card:
 
+                        for param   in param_card[block]:
+                            lhacode = param.lhacode
+                            value = param.value
+                            name = '%s_%s' % (block.upper(), '_'.join([str(i) for i in lhacode]))
+                            module.change_para(name, value)
+                    module.update_all_coup()
+                        
         return param_card_iterator, tag_name
 
         
@@ -961,8 +980,10 @@
     def calculate_weight(self, event):
         """space defines where to find the calculator (in multicore)"""
         
-        global lhapdf
-        
+        if not hasattr(self,'pdf'):
+            lhapdf = misc.import_python_lhapdf(self.mg5cmd.options['lhapdf'])
+            self.pdf = lhapdf.mkPDF(self.banner.run_card.get_lhapdf_id())
+            
         if self.has_nlo and self.rwgt_mode != "LO":
             return self.calculate_nlo_weight(event)
         
@@ -1075,19 +1096,20 @@
                     
                 base_wgt.append(c_wgt.pwgt[:3])
         
-        #change the ordering to the fortran one:
-        scales2 = self.invert_momenta(scales2)
-        pdg = self.invert_momenta(pdg)
-        bjx = self.invert_momenta(bjx)
+        
+        orig_wgt_check, partial_check = self.combine_wgt_local(scales2, pdg, bjx, base_wgt, gs, qcdpower, self.pdf)
+        #change the ordering to the fortran one: 
+        #scales2_i = self.invert_momenta(scales2)
+        #pdg_i = self.invert_momenta(pdg)
+        #bjx_i = self.invert_momenta(bjx)
         # re-compute original weight to reduce numerical inacurracy
-        base_wgt = self.invert_momenta(base_wgt)
-        
-        orig_wgt_check, partial_check = self.combine_wgt(scales2, pdg, bjx, base_wgt, gs, qcdpower, 1., 1.)
+        #base_wgt_i = self.invert_momenta(base_wgt)
+        #orig_wgt_check, partial_check = self.combine_wgt(scales2_i, pdg_i, bjx_i, base_wgt_i, gs, qcdpower, 1., 1.)
         
         if '_nlo' in type_nlo:
-            wgt = self.invert_momenta(wgt_virt)
+            #wgt = self.invert_momenta(wgt_virt)
             with misc.stdchannel_redirected(sys.stdout, os.devnull):
-                new_out, partial = self.combine_wgt(scales2, pdg, bjx, wgt, gs, qcdpower, 1., 1.)
+                new_out, partial = self.combine_wgt_local(scales2, pdg, bjx, wgt_virt, gs, qcdpower, self.pdf)
             # try to correct for precision issue
             avg = [partial_check[i]/ref_wgts[i] for i in range(len(ref_wgts))]
             out = sum(partial[i]/avg[i] if 0.85<avg[i]<1.15 else 0 \
@@ -1096,9 +1118,9 @@
 
             
         if '_tree' in type_nlo:
-            wgt = self.invert_momenta(wgt_tree)
+            #wgt = self.invert_momenta(wgt_tree)
             with misc.stdchannel_redirected(sys.stdout, os.devnull):
-                out, partial = self.combine_wgt(scales2, pdg, bjx, wgt, gs, qcdpower, 1., 1.)
+                out, partial = self.combine_wgt_local(scales2, pdg, bjx, wgt_tree, gs, qcdpower, self.pdf)
             # try to correct for precision issue
             avg = [partial_check[i]/ref_wgts[i] for i in range(len(ref_wgts))]
             new_out = sum(partial[i]/avg[i] if 0.85<avg[i]<1.15 else partial[i] \
@@ -1126,7 +1148,25 @@
             assert not to_write
             assert not wgt_tree
         return final_weight 
+
+
+    def combine_wgt_local(self, scale2s, pdgs, bjxs, base_wgts, gss, qcdpowers, pdf):
+
+        wgt = 0.
+        wgts = []
+        for (scale2, pdg, bjx, base_wgt, gs, qcdpower) in   zip(scale2s, pdgs, bjxs, base_wgts, gss, qcdpowers):
+            Q2, mur2, muf2 = scale2 #Q2 is Ellis-Sexton scale
+            #misc.sprint(Q2, mur2, muf2, base_wgt, gs, qcdpower)
+            pdf1 = pdf.xfxQ2(pdg[0], bjx[0], muf2)/bjx[0]
+            pdf2 = pdf.xfxQ2(pdg[1], bjx[1], muf2)/bjx[1]
+            alphas = pdf.alphasQ2(mur2)
+            tmp = base_wgt[0] + base_wgt[1] * math.log(mur2/Q2) + base_wgt[2] * math.log(muf2/Q2)
+            tmp *= gs**qcdpower*pdf1*pdf2
+            wgt += tmp
+            wgts.append(tmp)
+        return wgt, wgts
         
+    
      
     @staticmethod   
     def invert_momenta(p):
@@ -1450,10 +1490,9 @@
         commandline = re.sub('@\s*\d+', '', commandline)
         # deactivate golem since it creates troubles
         old_options = dict(mgcmd.options)
-        if mgcmd.options['golem'] or mgcmd.options['pjfry']:
-            logger.info(" When doing NLO reweighting, MG5aMC cannot use the loop reduction algorithms Golem and/or PJFry++")
+        if mgcmd.options['golem']:
+            logger.info(" When doing NLO reweighting, MG5aMC cannot use the loop reduction algorithms Golem")
         mgcmd.options['golem'] = None            
-        mgcmd.options['pjfry'] = None 
         commandline = commandline.replace('add process', 'generate',1)
         logger.info(commandline)
         mgcmd.exec_cmd(commandline, precmd=True)
@@ -1462,22 +1501,10 @@
         
         #put back golem to original value
         mgcmd.options['golem'] = old_options['golem']
-        mgcmd.options['pjfry'] = old_options['pjfry']
         # update make_opts
-        m_opts = {}
-        if mgcmd.options['lhapdf']:
-            #lhapdfversion = subprocess.Popen([mgcmd.options['lhapdf'], '--version'], 
-            #        stdout = subprocess.PIPE).stdout.read().strip()[0]
-            m_opts['lhapdf'] = True
-            m_opts['f2pymode'] = True
-            m_opts['lhapdfversion'] = 5 # 6 always fail on my computer since 5 is compatible but slower always use 5
-            m_opts['llhapdf'] = self.mother.get_lhapdf_libdir()                       
-        else:
+        if not mgcmd.options['lhapdf']:
             raise Exception, "NLO reweighting requires LHAPDF to work correctly"
 
-        path = pjoin(path_me,data['paths'][1], 'Source', 'make_opts')             
-        common_run_interface.CommonRunCmd.update_make_opts_full(path, m_opts)      
-        logger.info('Done %.4g' % (time.time()-start))
 
 
         # Download LHAPDF SET
@@ -1635,7 +1662,7 @@
         elif has_nlo and 'NLO' in self.rwgt_mode:
             self.create_standalone_virt_directory(data, second)
             
-            if not second:
+            if False:#not second:
                 #compile the module to combine the weight
                 misc.compile(cwd=pjoin(path_me, data['paths'][1], 'Source'))
                 #link it 
@@ -1666,8 +1693,7 @@
             commandline='import model loop_sm;generate g g > e+ ve [virt=QCD]'
             # deactivate golem since it creates troubles
             old_options = dict(mgcmd.options)
-            mgcmd.options['golem'] = None            
-            mgcmd.options['pjfry'] = None 
+            mgcmd.options['golem'] = None             
             commandline = commandline.replace('add process', 'generate',1)
             logger.info(commandline)
             mgcmd.exec_cmd(commandline, precmd=True)
@@ -1675,41 +1701,30 @@
             mgcmd.exec_cmd(commandline, precmd=True)    
             #put back golem to original value
             mgcmd.options['golem'] = old_options['golem']
-            mgcmd.options['pjfry'] = old_options['pjfry']
             # update make_opts
-            m_opts = {}
-            if mgcmd.options['lhapdf']:
-                #lhapdfversion = subprocess.Popen([mgcmd.options['lhapdf'], '--version'], 
-                #        stdout = subprocess.PIPE).stdout.read().strip()[0]
-                m_opts['lhapdf'] = True
-                m_opts['f2pymode'] = True
-                m_opts['lhapdfversion'] = 5 # 6 always fail on my computer since 5 is compatible but slower always use 5
-                m_opts['llhapdf'] = self.mother.get_lhapdf_libdir()                        
-            else:
+            if not mgcmd.options['lhapdf']:
                 raise Exception, "NLO_tree reweighting requires LHAPDF to work correctly"
  
-            path = pjoin(path_me,data['paths'][1], 'Source', 'make_opts')             
-            common_run_interface.CommonRunCmd.update_make_opts_full(path, m_opts)      
-            logger.info('Done %.4g' % (time.time()-start))
-
             # Download LHAPDF SET
             common_run_interface.CommonRunCmd.install_lhapdf_pdfset_static(\
                 mgcmd.options['lhapdf'], None, self.banner.run_card.get_lhapdf_id())
             
             #compile the module to combine the weight
-            misc.compile(cwd=pjoin(path_me, data['paths'][1], 'Source'))
-            #link it 
-            with misc.chdir(pjoin(path_me)):
-                if path_me not in sys.path:
-                    sys.path.insert(0, path_me)
-                mymod = __import__('%s.Source.rwgt2py' % data['paths'][1], globals(), locals(), [],-1)
-                mymod =  mymod.Source.rwgt2py
-                with misc.stdchannel_redirected(sys.stdout, os.devnull):
-                    mymod.initialise([self.banner.run_card['lpp1'], 
-                                  self.banner.run_card['lpp2']],
-                                 self.banner.run_card.get_lhapdf_id())
-                self.combine_wgt = mymod.get_wgt
-            
+            if False:
+                #use python module instead
+                misc.compile(cwd=pjoin(path_me, data['paths'][1], 'Source'))
+                #link it 
+                with misc.chdir(pjoin(path_me)):
+                    if path_me not in sys.path:
+                        sys.path.insert(0, path_me)
+                    mymod = __import__('%s.Source.rwgt2py' % data['paths'][1], globals(), locals(), [],-1)
+                    mymod =  mymod.Source.rwgt2py
+                    with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                        mymod.initialise([self.banner.run_card['lpp1'], 
+                                      self.banner.run_card['lpp2']],
+                                     self.banner.run_card.get_lhapdf_id())
+                    self.combine_wgt = mymod.get_wgt
+                
              
         # 6. If we need a new model/process-------------------------------------
         if (self.second_model or self.second_process or self.dedicated_path) and not second :
@@ -1895,6 +1910,7 @@
         to_save['has_nlo'] = self.has_nlo
         to_save['rwgt_mode'] = self.rwgt_mode
         to_save['rwgt_name'] = self.options['rwgt_name']
+        to_save['allow_missing_finalstate'] = self.options['allow_missing_finalstate']
 
         name = pjoin(self.rwgt_dir, 'rw_me', 'rwgt.pkl')
         save_load_object.save_to_file(name, to_save)
@@ -1910,7 +1926,7 @@
                         'rwgt_name': None}
         if keep_name:
             self.options['rwgt_name'] = obj['rwgt_name']
-        
+        self.options['allow_missing_finalstate'] = obj['allow_missing_finalstate']
         old_rwgt = obj['rwgt_dir']
            
         # path to fortran executable
@@ -1933,7 +1949,8 @@
         if not self.rwgt_mode:
             self.rwgt_mode = obj['rwgt_mode']
             logger.info("mode set to %s" % self.rwgt_mode)
-        if self.has_nlo and 'NLO' in self.rwgt_mode:
+        if False:#self.has_nlo and 'NLO' in self.rwgt_mode:
+            #use python version
             path = pjoin(obj['rwgt_dir'], 'rw_mevirt','Source')
             sys.path.insert(0, path)
             try:

=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2020-05-11 16:36:37 +0000
+++ madgraph/iolibs/export_v4.py	2020-07-04 19:03:17 +0000
@@ -12,7 +12,6 @@
 # For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
 #
 ################################################################################
-from madgraph.iolibs.helas_call_writers import HelasCallWriter
 """Methods and classes to export matrix elements to v4 format."""
 
 import copy
@@ -881,7 +880,7 @@
         # Make sure aloha is in quadruple precision if needed
         old_aloha_mp=aloha.mp_precision
         aloha.mp_precision=self.opt['mp']
-
+        self.model = model
         # create the MODEL
         write_dir=pjoin(self.dir_path, 'Source', 'MODEL')
         model_builder = UFO_model_to_mg4(model, write_dir, self.opt + self.proc_characteristic)
@@ -2040,7 +2039,7 @@
         open(pjoin(self.dir_path,'__init__.py'),'w')
         open(pjoin(self.dir_path,'SubProcesses','__init__.py'),'w')
 
-        if 'mode' in self.opt and self.opt['mode'] == "reweight":
+        if False:#'mode' in self.opt and self.opt['mode'] == "reweight":
             #add the module to hande the NLO weight
             files.copytree(pjoin(MG5DIR, 'Template', 'RWGTNLO'),
                           pjoin(self.dir_path, 'Source'))
@@ -2048,7 +2047,7 @@
                            pjoin(self.dir_path, 'Source', 'PDF'))
             self.write_pdf_opendata()
             
-        if self.prefix_info:
+        if self.prefix_info: 
             self.write_f2py_splitter()
             self.write_f2py_makefile()
             self.write_f2py_check_sa(matrix_elements,
@@ -2105,6 +2104,34 @@
       CALL SETPARA(PATH)  !first call to setup the paramaters
       RETURN
       END
+      
+      
+      subroutine CHANGE_PARA(name, value)
+      implicit none
+CF2PY intent(in) :: name
+CF2PY intent(in) :: value
+
+      character*512 name
+      double precision value
+
+      include '../Source/MODEL/input.inc'
+      include '../Source/MODEL/coupl.inc'
+
+      SELECT CASE (name)
+         %(parameter_setup)s
+         CASE DEFAULT
+            write(*,*) 'no parameter matching', name, value
+      END SELECT
+
+      return
+      end
+      
+    subroutine update_all_coup()
+    implicit none
+     call coup()
+    return 
+    end
+      
 
     subroutine get_pdg_order(PDG)
   IMPLICIT NONE
@@ -2159,20 +2186,46 @@
         if min_nexternal != max_nexternal:
             text.append('endif')
 
+        params = self.get_model_parameter(self.model)
+        parameter_setup =[]
+        for key, var in params.items():
+            parameter_setup.append('        CASE ("%s")\n          %s = value' 
+                                   % (key, var))
+
         formatting = {'python_information':'\n'.join(info), 
                           'smatrixhel': '\n'.join(text),
                           'maxpart': max_nexternal,
                           'nb_me': len(allids),
                           'pdgs': ','.join(str(pdg[i]) if i<len(pdg) else '0' 
                                              for i in range(max_nexternal) for pdg in allids),
-                          'prefix':'\',\''.join(allprefix)
+                          'prefix':'\',\''.join(allprefix),
+                          'parameter_setup': '\n'.join(parameter_setup),
                           }
         formatting['lenprefix'] = len(formatting['prefix'])
         text = template % formatting
         fsock = writers.FortranWriter(pjoin(self.dir_path, 'SubProcesses', 'all_matrix.f'),'w')
         fsock.writelines(text)
         fsock.close()
+    
+    def get_model_parameter(self, model):
+        """ returns all the model parameter
+        """
+        params = {}
+        for p in model.get('parameters')[('external',)]:
+            name = p.name
+            nopref = name[4:] if name.startswith('mdl_') else name
+            params[nopref] = name
             
+            block = p.lhablock
+            lha = '_'.join([str(i) for i in p.lhacode])
+            params['%s_%s' % (block.upper(), lha)] = name
+
+        return params                      
+                                        
+        
+        
+        
+         
     def write_f2py_check_sa(self, matrix_element, writer):
         """ Write the general check_sa.py in SubProcesses that calls all processes successively."""
         # To be implemented. It is just an example file, i.e. not crucial.

=== modified file 'madgraph/iolibs/file_writers.py'
--- madgraph/iolibs/file_writers.py	2018-05-03 09:00:52 +0000
+++ madgraph/iolibs/file_writers.py	2020-07-04 19:03:17 +0000
@@ -342,19 +342,28 @@
         columns. Split in preferential order according to
         split_characters, and start each new line with line_start."""
 
-        res_lines = [line]
-
-        while len(res_lines[-1]) > self.line_length:
+        def get_split_index(line, max_length, max_split, split_characters):
             split_at = 0
             for character in split_characters:
-                index = res_lines[-1][(self.line_length - self.max_split): \
-                                      self.line_length].rfind(character)
+                index = line[(max_length - max_split): \
+                                      max_length].rfind(character)
                 if index >= 0:
-                    split_at_tmp = self.line_length - self.max_split + index
+                    split_at_tmp = max_length - max_split + index
                     if split_at_tmp > split_at:
                         split_at = split_at_tmp
+            return split_at
+                
+
+        res_lines = [line]
+
+        while len(res_lines[-1]) > self.line_length:
+            split_at = get_split_index(res_lines[-1], self.line_length, 
+                                       self.max_split, split_characters)
             if split_at == 0:
-                split_at = self.line_length
+                split_at = get_split_index(res_lines[-1], self.line_length, 
+                                       self.max_split + 30, split_characters)
+                if split_at == 0:
+                    split_at = self.line_length
                 
             newline = res_lines[-1][split_at:]
             nquotes = self.count_number_of_quotes(newline)

=== modified file 'madgraph/iolibs/files.py'
--- madgraph/iolibs/files.py	2018-05-18 11:55:37 +0000
+++ madgraph/iolibs/files.py	2020-07-04 19:03:17 +0000
@@ -36,7 +36,8 @@
             ret_value = myfunct(sock, *args)
         finally:
             sock.close()
-    except IOError, (errno, strerror):
+    except IOError as error:
+        errno, strerror = error.errno, str(error)
         if opt.has_key('print_error'):
             if not opt['print_error']:
                 return None

=== modified file 'madgraph/iolibs/template_files/madevent_symmetry.f'
--- madgraph/iolibs/template_files/madevent_symmetry.f	2018-03-16 12:13:34 +0000
+++ madgraph/iolibs/template_files/madevent_symmetry.f	2020-07-04 19:03:17 +0000
@@ -103,7 +103,7 @@
          if (abs(lpp(1)) .eq. 3) m1 = 0.000511d0
          if (abs(lpp(2)) .eq. 3) m2 = 0.000511d0
          if (mass_ion(1).ge.0d0) m1 = mass_ion(1)
-         if (mass_ion(2).ge.0d0) m1 = mass_ion(2)
+         if (mass_ion(2).ge.0d0) m2 = mass_ion(2)
          if(ebeam(1).lt.m1.and.lpp(1).ne.9) ebeam(1)=m1
          if(ebeam(2).lt.m2.and.lpp(2).ne.9) ebeam(2)=m2
          pi1(0)=ebeam(1)

=== modified file 'madgraph/loop/loop_exporters.py'
--- madgraph/loop/loop_exporters.py	2019-08-28 08:02:11 +0000
+++ madgraph/loop/loop_exporters.py	2020-07-04 19:03:17 +0000
@@ -297,6 +297,36 @@
       RETURN
       END
 
+      subroutine CHANGE_PARA(name, value)
+      implicit none
+CF2PY intent(in) :: name
+CF2PY intent(in) :: value
+
+      character*512 name
+      double precision value
+
+      include '../Source/MODEL/input.inc'
+      include '../Source/MODEL/coupl.inc'
+      include '../Source/MODEL/mp_coupl.inc'
+      include '../Source/MODEL/mp_input.inc'
+      
+      SELECT CASE (name)   
+         %(parameter_setup)s
+         CASE DEFAULT
+            write(*,*) 'no parameter matching', name
+      END SELECT
+
+      return
+      end
+      
+    subroutine update_all_coup()
+    implicit none
+     call coup()
+     call printout()
+    return 
+    end
+
+
       SUBROUTINE SET_MADLOOP_PATH(PATH)
 C     Routine to set the path of the folder 'MadLoop5_resources' to MadLoop
         CHARACTER(512) PATH
@@ -377,6 +407,14 @@
         #close the function
         if min_nexternal != max_nexternal:
             text.append('endif')
+            
+        params = self.get_model_parameter(self.model)
+        parameter_setup =[]
+        for key, var in params.items():
+            parameter_setup.append('        CASE ("%s")\n          %s = value\n        MP__%s = value' 
+                                   % (key, var, var))
+            
+            
 
         formatting = {'python_information':'\n'.join(info), 
                           'smatrixhel': '\n'.join(text),
@@ -385,7 +423,8 @@
                           'pdgs': ','.join([str(pdg[i]) if i<len(pdg) else '0' 
                                              for i in range(max_nexternal) \
                                              for pdg in allids]),
-                      'prefix':'\',\''.join(allprefix)
+                      'prefix':'\',\''.join(allprefix),
+                      'parameter_setup': '\n'.join(parameter_setup),
                       }
     
     

=== modified file 'madgraph/madevent/gen_crossxhtml.py'
--- madgraph/madevent/gen_crossxhtml.py	2018-05-25 09:13:24 +0000
+++ madgraph/madevent/gen_crossxhtml.py	2020-07-04 19:03:17 +0000
@@ -375,6 +375,9 @@
             run[name] = int(value)
         elif name in ['run_mode','run_statistics']:
             run[name] = value
+        elif name == 'cross' and run[name] != 0:
+            run['prev_' + name] = run[name]
+            run[name] = float(value) 
         else:    
             run[name] = float(value)    
     

=== modified file 'madgraph/madevent/gen_ximprove.py'
--- madgraph/madevent/gen_ximprove.py	2020-02-25 14:05:36 +0000
+++ madgraph/madevent/gen_ximprove.py	2020-07-04 19:03:17 +0000
@@ -108,6 +108,10 @@
             self.splitted_grid = self.run_card['survey_splitting']
         if self.run_card['survey_nchannel_per_job'] != -1:
             self.combining_job = self.run_card['survey_nchannel_per_job']        
+        elif self.run_card['hard_survey'] > 1:
+            self.combining_job = 1
+            
+            
         
         self.splitted_Pdir = {}
         self.splitted_for_dir = lambda x,y: self.splitted_grid
@@ -971,7 +975,7 @@
         super(gen_ximprove_v4, self).__init__(cmd, opt)
         
         if cmd.opts['accuracy'] < cmd._survey_options['accuracy'][1]:
-            self.increase_precision()
+            self.increase_precision(cmd._survey_options['accuracy'][1]/cmd.opts['accuracy'])
 
     def reset_multijob(self):
 
@@ -986,15 +990,23 @@
         f.write('%i\n' % nb_split)
         f.close()
     
-    def increase_precision(self):
-        
-        self.max_event_in_iter = 20000
-        self.min_events = 7500
+    def increase_precision(self, rate=3):
+        misc.sprint(rate)
+        if rate < 3:
+            self.max_event_in_iter = 20000
+            self.min_events = 7500
+            self.gen_events_security = 1.3
+        else:
+            rate = rate -2
+            self.max_event_in_iter = int((rate+1) * 10000)
+            self.min_events = int(rate+2) * 2500
+            self.gen_events_security = 1 + 0.1 * (rate+2) 
+                        
         if int(self.nhel) == 1:
             self.min_event_in_iter *= 2**(self.cmd.proc_characteristics['nexternal']//3)
             self.max_event_in_iter *= 2**(self.cmd.proc_characteristics['nexternal']//2)
 
-        self.gen_events_security = 1.3
+        
             
     alphabet = "abcdefghijklmnopqrstuvwxyz"
     def get_job_for_event(self):

=== modified file 'madgraph/various/banner.py'
--- madgraph/various/banner.py	2020-05-23 19:22:37 +0000
+++ madgraph/various/banner.py	2020-07-04 19:03:17 +0000
@@ -3093,7 +3093,8 @@
         self.add_param('issgridfile', '', hidden=True)
         #job handling of the survey/ refine
         self.add_param('job_strategy', 0, hidden=True, include=False, allowed=[0,1,2], comment='see appendix of 1507.00020 (page 26)')
-        self.add_param('hard_survey', False, hidden=True, include=False, comment='force to have better estimate of the integral at survey for difficult mode like VBF')
+        self.add_param('hard_survey', 0, hidden=True, include=False, comment='force to have better estimate of the integral at survey for difficult mode like VBF')
+        self.add_param("second_refine_treshold", 0.9, hidden=True, include=False, comment="set a treshold to bypass the use of a second refine. if the ratio of cross-section after survey by the one of the first refine is above the treshold, the  second refine will not be done.")
         self.add_param('survey_splitting', -1, hidden=True, include=False, comment="for loop-induced control how many core are used at survey for the computation of a single iteration.")
         self.add_param('survey_nchannel_per_job', 2, hidden=True, include=False, comment="control how many Channel are integrated inside a single job on cluster/multicore")
         self.add_param('refine_evt_by_job', -1, hidden=True, include=False, comment="control the maximal number of events for the first iteration of the refine (larger means less jobs)")
@@ -4072,7 +4073,7 @@
         super(RunCardNLO, self).check_validity()
 
         # for lepton-lepton collisions, ignore 'pdlabel' and 'lhaid'
-        if self['lpp1']!=1 or self['lpp2']!=1:
+        if abs(self['lpp1'])!=1 or abs(self['lpp2'])!=1:
             if self['lpp1'] == 1 or self['lpp2']==1:
                 raise InvalidRunCard('Process like Deep Inelastic scattering not supported at NLO accuracy.')
             

=== modified file 'madgraph/various/histograms.py'
--- madgraph/various/histograms.py	2019-08-28 08:02:11 +0000
+++ madgraph/various/histograms.py	2020-07-04 19:03:17 +0000
@@ -2397,7 +2397,11 @@
             gnuplot_output_list=gnuplot_output_list_v5
         else:
             output, _ = p.communicate()
-            if float(output.split()[1]) < 5. :
+            try:
+                version = float(output.split()[1])
+            except:
+                version = 5
+            if version < 5. :
                 gnuplot_output_list=gnuplot_output_list_v4
             else:
                 gnuplot_output_list=gnuplot_output_list_v5

=== modified file 'madgraph/various/lhe_parser.py'
--- madgraph/various/lhe_parser.py	2020-01-09 23:14:56 +0000
+++ madgraph/various/lhe_parser.py	2020-07-04 19:03:17 +0000
@@ -1449,8 +1449,8 @@
         if not hasattr(Event, 'loweight_pattern'):
             Event.loweight_pattern = re.compile('''<rscale>\s*(?P<nqcd>\d+)\s+(?P<ren_scale>[\d.e+-]+)\s*</rscale>\s*\n\s*
                                     <asrwt>\s*(?P<asrwt>[\s\d.+-e]+)\s*</asrwt>\s*\n\s*
-                                    <pdfrwt\s+beam=["']?1["']?\>\s*(?P<beam1>[\s\d.e+-]*)\s*</pdfrwt>\s*\n\s*
-                                    <pdfrwt\s+beam=["']?2["']?\>\s*(?P<beam2>[\s\d.e+-]*)\s*</pdfrwt>\s*\n\s*
+                                    <pdfrwt\s+beam=["']?(?P<idb1>1|2)["']?\>\s*(?P<beam1>[\s\d.e+-]*)\s*</pdfrwt>\s*\n\s*
+                                    <pdfrwt\s+beam=["']?(?P<idb2>1|2)["']?\>\s*(?P<beam2>[\s\d.e+-]*)\s*</pdfrwt>\s*\n\s*
                                     <totfact>\s*(?P<totfact>[\d.e+-]*)\s*</totfact>
             ''',re.X+re.I+re.M)
         
@@ -1468,13 +1468,22 @@
             self.loweight['asrwt'] =[float(x) for x in info.group('asrwt').split()[1:]]
             self.loweight['tot_fact'] = float(info.group('totfact'))
             
-            args = info.group('beam1').split()
+            if info.group('idb1') == info.group('idb2'):
+                raise Exception, '%s not parsed'% text
+            
+            if info.group('idb1') =="1":
+                args = info.group('beam1').split()
+            else:
+                args = info.group('beam2').split()
             npdf = int(args[0])
             self.loweight['n_pdfrw1'] = npdf
             self.loweight['pdf_pdg_code1'] = [int(i) for i in args[1:1+npdf]]
             self.loweight['pdf_x1'] = [float(i) for i in args[1+npdf:1+2*npdf]]
             self.loweight['pdf_q1'] = [float(i) for i in args[1+2*npdf:1+3*npdf]]
-            args = info.group('beam2').split()
+            if info.group('idb2') =="2":
+                args = info.group('beam2').split()
+            else:
+                args = info.group('beam1').split()
             npdf = int(args[0])
             self.loweight['n_pdfrw2'] = npdf
             self.loweight['pdf_pdg_code2'] = [int(i) for i in args[1:1+npdf]]

=== modified file 'madgraph/various/misc.py'
--- madgraph/various/misc.py	2020-01-08 14:13:32 +0000
+++ madgraph/various/misc.py	2020-07-04 19:03:17 +0000
@@ -488,7 +488,7 @@
             error_text += "In general this means that your computer is not able to compile."
             if sys.platform == "darwin":
                 error_text += "Note that MacOSX doesn\'t have gmake/gfortan install by default.\n"
-                error_text += "Xcode3 contains those required programs"
+                error_text += "Xcode contains gmake. For gfortran we advise: http://hpc.sourceforge.net/"
             raise MadGraph5Error, error_text
 
     if p.returncode:
@@ -539,7 +539,7 @@
         p = Popen([compiler, '-dumpversion'], stdout=subprocess.PIPE, 
                     stderr=subprocess.PIPE)
         output, error = p.communicate()
-        version_finder=re.compile(r"(?P<version>(\d.)*\d)")
+        version_finder=re.compile(r"(?P<version>\d[\d.]*)")
         version = version_finder.search(output).group('version')
         return version
     except Exception:
@@ -1772,7 +1772,8 @@
                         import random
                         msg = choices[random.randint(0,len(choices)-2)]
                     EasterEgg.message_aprilfirst[msgtype].remove(msg)
-                    
+                else:
+                    return
             elif msgtype=='loading' and date in self.special_banner:
                 self.change_banner(date)
                 return
@@ -2065,7 +2066,9 @@
                 logger.warning("Failed to access python version of LHAPDF: "\
                                    "If the python interface to LHAPDF is available on your system, try "\
                                    "adding its location to the PYTHONPATH environment variable and the"\
-                                   "LHAPDF library location to LD_LIBRARY_PATH (linux) or DYLD_LIBRARY_PATH (mac os x).")
+                                   "LHAPDF library location to LD_LIBRARY_PATH (linux) or DYLD_LIBRARY_PATH (mac os x)."\
+                                   "The required LD_LIBRARY_PATH is "+ lhapdf_libdir 
+                                   )
         
     if use_lhapdf:
         python_lhapdf = lhapdf

=== modified file 'madgraph/various/process_checks.py'
--- madgraph/various/process_checks.py	2019-03-01 08:53:21 +0000
+++ madgraph/various/process_checks.py	2020-07-04 19:03:17 +0000
@@ -375,7 +375,7 @@
         else:
             energy = options['energy']
             events = options['events']
-            to_skip = 0
+            to_skip = options['skip_evt']
             
         if not (isinstance(process, base_objects.Process) and \
                 isinstance(energy, (float,int))):

=== modified file 'models/model_reader.py'
--- models/model_reader.py	2018-09-21 13:46:51 +0000
+++ models/model_reader.py	2020-07-04 19:03:17 +0000
@@ -180,6 +180,11 @@
                                     value = 0.0
                                 else:
                                     raise
+                            except OverflowError, err:
+                                if scale < 1:
+                                    value = 0.0
+                                else:
+                                    raise
                         exec("locals()[\'%s\'] = %s" % (parameter_dict[block][pid].name,
                                           value))
                         parameter_dict[block][pid].value = float(value)

=== modified file 'tests/acceptance_tests/test_cmd_reweight.py'
--- tests/acceptance_tests/test_cmd_reweight.py	2018-06-01 12:14:52 +0000
+++ tests/acceptance_tests/test_cmd_reweight.py	2020-07-04 19:03:17 +0000
@@ -136,16 +136,22 @@
         ff.close()
         
 
-        with misc.stdchannel_redirected(sys.stdout, os.devnull):
-            me_cmd.run_cmd('reweight run_01 --from_cards')
+        if logger.level <= 10:
+                me_cmd.run_cmd('reweight run_01 --from_cards')
+        else:
+            with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                me_cmd.run_cmd('reweight run_01 --from_cards')
         
         lhe = lhe_parser.EventFile(pjoin(self.run_dir,'Events','run_01', 'unweighted_events.lhe.gz'))
         
-        solutions = [12.905371, 12.271068, 12.271753, 12.436969, 12.330121, 12.441459, 12.290359, 12.406486, 12.352107, 11.862291, 12.346005, 12.603302, 12.483697, 12.597262, 12.352434, 11.841136, 12.559252, 12.693822, 12.073114, 12.080773, 12.088188, 12.169921, 12.36124, 12.620677, 12.641769]
+        #solutions = [12.905371, 12.271068, 12.271753, 12.436969, 12.330121, 12.441459, 12.290359, 12.406486, 12.352107, 11.862291, 12.346005, 12.603302, 12.483697, 12.597262, 12.352434, 11.841136, 12.559252, 12.693822, 12.073114, 12.080773, 12.088188, 12.169921, 12.36124, 12.620677, 12.641769]
+        solutions = [14.288834, 13.551237, 13.564131, 13.745311, 13.620339, 13.750088, 13.611093, 13.709763, 13.670902, 13.071603, 13.658932, 13.939627, 13.827959, 13.932638, 13.669482, 13.045913, 13.888161, 14.045136, 13.319103, 13.327884, 13.336799, 13.430042, 13.658435, 13.959817, 13.984111]
+        #solutions = []
         for i,event in enumerate(lhe):
             
             rwgt_data = event.parse_reweight()
             #solutions.append(rwgt_data['rwgt_1'])
+            #continue
             self.assertTrue('rwgt_1' in rwgt_data)
             self.assertTrue(misc.equal(rwgt_data['rwgt_1'], solutions[i]))
         #misc.sprint(solutions)
@@ -239,7 +245,8 @@
         if 1: #with misc.stdchannel_redirected(sys.stdout, os.devnull):
             me_cmd.run_cmd('reweight run_01 --from_cards')
         
-        solutions = [41.511565, 41.930505, 41.511565, 41.511565, 41.586169, 41.511565, 41.511565, 41.511565, 41.511565, 42.046806, 41.511565, 44.164503, 41.511565, -41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.778368, 41.511565, 42.05448, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 42.115494, 41.511854, 41.511567, 41.511565, 42.00028, 42.120605, 41.514867, -41.511565, 41.511565, 45.125706, 41.511565, 42.180208, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.997509, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.513926, 41.511565, 41.882499, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 47.532736, 41.511565, 41.809063, 41.511565, 41.511565, 41.511565, 41.927695, 41.511565, 41.511565, 41.555692, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, -41.511565, 41.511565, 42.029675, 41.725129, 41.511565, 41.511565, 41.511778, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 42.050439, 41.511565, 41.511565, 41.511565, -41.511565, 42.071105, 41.511565, 41.511565]
+        #solutions = [41.511565, 41.930505, 41.511565, 41.511565, 41.586169, 41.511565, 41.511565, 41.511565, 41.511565, 42.046806, 41.511565, 44.164503, 41.511565, -41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.778368, 41.511565, 42.05448, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 42.115494, 41.511854, 41.511567, 41.511565, 42.00028, 42.120605, 41.514867, -41.511565, 41.511565, 45.125706, 41.511565, 42.180208, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.997509, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.513926, 41.511565, 41.882499, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 47.532736, 41.511565, 41.809063, 41.511565, 41.511565, 41.511565, 41.927695, 41.511565, 41.511565, 41.555692, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, -41.511565, 41.511565, 42.029675, 41.725129, 41.511565, 41.511565, 41.511778, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 42.050439, 41.511565, 41.511565, 41.511565, -41.511565, 42.071105, 41.511565, 41.511565]
+        solutions =  [41.511565, 41.75164, 41.511565, 41.511565, 41.557009, 41.511565, 41.511565, 41.511565, 41.511565, 41.759881, 41.511565, 44.717244, 41.511565, -41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 42.041827, 41.511565, 41.663092, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.920917, 41.511736, 41.511566, 41.511565, 205.29791, 42.538389, 41.513573, -41.511565, 41.511565, 44.650974, 41.511565, 42.549187, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.431553, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.513265, 41.511565, 43.902643, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 8853.0397, 41.511565, 333.38491, 41.511565, 41.511565, 41.511565, 75.701362, 41.511565, 41.511565, 41.54695, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, -41.511565, 41.511565, 41.454923, 41.53984, 41.511565, 41.511565, 41.511743, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.732627, 41.511565, 41.511565, 41.511565, -41.511565, 41.395865, 41.511565, 41.511565]
         lhe = lhe_parser.EventFile(pjoin(self.run_dir,'Events','run_01', 'rwgt_events_tree_rwgt_1.lhe.gz'))
         lhe_orig = lhe_parser.EventFile(pjoin(self.run_dir,'Events','run_01', 'events.lhe.gz'))
         from itertools import izip
@@ -249,9 +256,9 @@
             i+=1
             rwgt_data = event_orig.parse_reweight()
             #solutions.append(rwgt_data['rwgt_1_tree'])
+            #continue
             self.assertTrue('rwgt_1_tree' in rwgt_data)
             self.assertEqual(event.wgt, rwgt_data['rwgt_1_tree'])
-            #misc.sprint(abs(event.wgt))
             self.assertTrue(misc.equal(event.wgt, solutions[i]))
             nlo1 = event.parse_nlo_weight()
             nlo2 = event_orig.parse_nlo_weight()
@@ -259,6 +266,7 @@
             self.assertNotEqual(nlo1.total_wgt, nlo2.total_wgt)
 
         #misc.sprint(solutions)
+        #raise Exception
 
     def test_loop_improved_reweighting(self):
         """ check identical re-weighting in ttbar 
@@ -281,8 +289,11 @@
         ff.write(cmd_lines)
         ff.close()
         
-        with misc.stdchannel_redirected(sys.stdout, os.devnull):
+        if logger.level <=10:
             me_cmd.run_cmd('reweight run_01 --from_cards')
+        else:
+            with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                me_cmd.run_cmd('reweight run_01 --from_cards')
         
         lhe = lhe_parser.EventFile(pjoin(self.run_dir,'Events','run_01', 'events.lhe.gz'))
         lhe2=  lhe_parser.EventFile(pjoin(self.run_dir,'Events','run_01', 'rwgt_events_tree_rwgt_1.lhe.gz'))
@@ -335,11 +346,12 @@
             self.assertTrue(misc.equal(rwgt_data['NAME_1'], solutions2[i]))
             
 
-    def test_nlo_reweighting_comb(self):
+    def old_test_nlo_reweighting_comb(self):
         """check that nlo reweighting is working.
            The main point is to check the recombination of the weights
            Since the rest should be either checked by the lhe_parser class
-           or by the various check of the standalone checks
+           or by the various check of the standalone checks.
+           This way of combining weights is now outdated
         """
         
         # create a reweight directory   

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_001.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_001.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_001.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_002.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_002.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_002.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_003.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_003.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_003.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_004.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_004.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_004.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_005.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_005.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_005.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_006.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_006.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_006.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_007.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_007.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_007.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_008.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_008.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%b_sf_008.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%born.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%born.f	2020-07-04 19:03:17 +0000
@@ -131,8 +131,8 @@
           DO J = 1, NSPLITORDERS
             AMP_ORDERS(J) = GETORDPOWFROMINDEX_B(J, I)
           ENDDO
-          IF(ABS(ANS(1,I)).GT.MAX_VAL*TINY) AMP_SPLIT(ORDERS_TO_AMP_SPL
-     $IT_POS(AMP_ORDERS)) = ANS(1,I)
+          IF(ABS(ANS(1,I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS)) = ANS(1,I)
         ENDIF
       ENDDO
 C     this is to avoid fake non-zero contributions 
@@ -169,10 +169,12 @@
               IF (K.EQ.J) AMP_ORDERS(K) = AMP_ORDERS(K) + 2
             ENDDO
 C           this is to avoid fake non-zero contributions 
-            IF (ABS(ANS(1,I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_
-     $AMP_SPLIT_POS(AMP_ORDERS),1,J) = ANS(1,I)
-            IF (ABS(ANS(2,I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_
-     $AMP_SPLIT_POS(AMP_ORDERS),2,J) = ANS(2,I)
+            IF (ABS(ANS(1,I)).GT.MAX_VAL*TINY)
+     $        AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1,J) =
+     $        ANS(1,I)
+            IF (ABS(ANS(2,I)).GT.MAX_VAL*TINY)
+     $        AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),2,J) =
+     $        ANS(2,I)
           ENDIF
         ENDDO
 C       this is to avoid fake non-zero contributions 
@@ -599,8 +601,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSPLITORDERS
-        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)+AMPSPLITORDERS(AMPORDER
-     $B,I)
+        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)
+     $   +AMPSPLITORDERS(AMPORDERB,I)
       ENDDO
       SQSOINDEXB=SQSOINDEXB_FROM_ORDERS(SQORDERS)
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%sborn_sf.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%sborn_sf.f	2017-07-31 08:48:20 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_gg_ttx%sborn_sf.f	2020-07-04 19:03:17 +0000
@@ -64,8 +64,8 @@
 
         WGT = WGT_COL * G**2
 C       update the amp_split_soft, which is summed in sbornsoft
-        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) = DBLE(AMP_SPLIT_CNT(1:AMP_SPL
-     $IT_SIZE,1,QCD_POS)) * G**2
+        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) =
+     $    DBLE(AMP_SPLIT_CNT(1:AMP_SPLIT_SIZE,1,QCD_POS)) * G**2
 
       ELSE IF (NEED_CHARGE_LINKS) THEN
         CHARGEPROD = CHARGES_BORN(M) * CHARGES_BORN(N)
@@ -76,8 +76,9 @@
         IF (M.EQ.N) CHARGEPROD = CHARGEPROD / 2D0
         WGT = DBLE(ANS_CNT(1, QED_POS)) * CHARGEPROD * DBLE(GAL(1))**2
 C       update the amp_split_soft, which is summed in sbornsoft
-        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) = DBLE(AMP_SPLIT_CNT(1:AMP_SPL
-     $IT_SIZE,1,QED_POS)) * CHARGEPROD * DBLE(GAL(1))**2
+        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) =
+     $    DBLE(AMP_SPLIT_CNT(1:AMP_SPLIT_SIZE,1,QED_POS)) * CHARGEPROD
+     $    * DBLE(GAL(1))**2
       ENDIF
 
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_001.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_001.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_001.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_002.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_002.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_002.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_003.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_003.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_003.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_004.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_004.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_004.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_005.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_005.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_005.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_006.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_006.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_006.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_007.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_007.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_007.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_008.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_008.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%b_sf_008.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%born.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%born.f	2020-07-04 19:03:17 +0000
@@ -134,8 +134,8 @@
           DO J = 1, NSPLITORDERS
             AMP_ORDERS(J) = GETORDPOWFROMINDEX_B(J, I)
           ENDDO
-          IF(ABS(ANS(1,I)).GT.MAX_VAL*TINY) AMP_SPLIT(ORDERS_TO_AMP_SPL
-     $IT_POS(AMP_ORDERS)) = ANS(1,I)
+          IF(ABS(ANS(1,I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS)) = ANS(1,I)
         ENDIF
       ENDDO
 C     this is to avoid fake non-zero contributions 
@@ -172,10 +172,12 @@
               IF (K.EQ.J) AMP_ORDERS(K) = AMP_ORDERS(K) + 2
             ENDDO
 C           this is to avoid fake non-zero contributions 
-            IF (ABS(ANS(1,I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_
-     $AMP_SPLIT_POS(AMP_ORDERS),1,J) = ANS(1,I)
-            IF (ABS(ANS(2,I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_
-     $AMP_SPLIT_POS(AMP_ORDERS),2,J) = ANS(2,I)
+            IF (ABS(ANS(1,I)).GT.MAX_VAL*TINY)
+     $        AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1,J) =
+     $        ANS(1,I)
+            IF (ABS(ANS(2,I)).GT.MAX_VAL*TINY)
+     $        AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),2,J) =
+     $        ANS(2,I)
           ENDIF
         ENDDO
 C       this is to avoid fake non-zero contributions 
@@ -602,8 +604,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSPLITORDERS
-        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)+AMPSPLITORDERS(AMPORDER
-     $B,I)
+        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)
+     $   +AMPSPLITORDERS(AMPORDERB,I)
       ENDDO
       SQSOINDEXB=SQSOINDEXB_FROM_ORDERS(SQORDERS)
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%sborn_sf.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%sborn_sf.f	2017-07-31 08:48:20 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uux_ttx%sborn_sf.f	2020-07-04 19:03:17 +0000
@@ -64,8 +64,8 @@
 
         WGT = WGT_COL * G**2
 C       update the amp_split_soft, which is summed in sbornsoft
-        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) = DBLE(AMP_SPLIT_CNT(1:AMP_SPL
-     $IT_SIZE,1,QCD_POS)) * G**2
+        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) =
+     $    DBLE(AMP_SPLIT_CNT(1:AMP_SPLIT_SIZE,1,QCD_POS)) * G**2
 
       ELSE IF (NEED_CHARGE_LINKS) THEN
         CHARGEPROD = CHARGES_BORN(M) * CHARGES_BORN(N)
@@ -76,8 +76,9 @@
         IF (M.EQ.N) CHARGEPROD = CHARGEPROD / 2D0
         WGT = DBLE(ANS_CNT(1, QED_POS)) * CHARGEPROD * DBLE(GAL(1))**2
 C       update the amp_split_soft, which is summed in sbornsoft
-        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) = DBLE(AMP_SPLIT_CNT(1:AMP_SPL
-     $IT_SIZE,1,QED_POS)) * CHARGEPROD * DBLE(GAL(1))**2
+        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) =
+     $    DBLE(AMP_SPLIT_CNT(1:AMP_SPLIT_SIZE,1,QED_POS)) * CHARGEPROD
+     $    * DBLE(GAL(1))**2
       ENDIF
 
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_001.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_001.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_001.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_002.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_002.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_002.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_003.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_003.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_003.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_004.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_004.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_004.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_005.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_005.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_005.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_006.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_006.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_006.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_007.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_007.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_007.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_008.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_008.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%b_sf_008.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%born.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%born.f	2020-07-04 19:03:17 +0000
@@ -134,8 +134,8 @@
           DO J = 1, NSPLITORDERS
             AMP_ORDERS(J) = GETORDPOWFROMINDEX_B(J, I)
           ENDDO
-          IF(ABS(ANS(1,I)).GT.MAX_VAL*TINY) AMP_SPLIT(ORDERS_TO_AMP_SPL
-     $IT_POS(AMP_ORDERS)) = ANS(1,I)
+          IF(ABS(ANS(1,I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS)) = ANS(1,I)
         ENDIF
       ENDDO
 C     this is to avoid fake non-zero contributions 
@@ -172,10 +172,12 @@
               IF (K.EQ.J) AMP_ORDERS(K) = AMP_ORDERS(K) + 2
             ENDDO
 C           this is to avoid fake non-zero contributions 
-            IF (ABS(ANS(1,I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_
-     $AMP_SPLIT_POS(AMP_ORDERS),1,J) = ANS(1,I)
-            IF (ABS(ANS(2,I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_
-     $AMP_SPLIT_POS(AMP_ORDERS),2,J) = ANS(2,I)
+            IF (ABS(ANS(1,I)).GT.MAX_VAL*TINY)
+     $        AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1,J) =
+     $        ANS(1,I)
+            IF (ABS(ANS(2,I)).GT.MAX_VAL*TINY)
+     $        AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),2,J) =
+     $        ANS(2,I)
           ENDIF
         ENDDO
 C       this is to avoid fake non-zero contributions 
@@ -602,8 +604,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSPLITORDERS
-        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)+AMPSPLITORDERS(AMPORDER
-     $B,I)
+        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)
+     $   +AMPSPLITORDERS(AMPORDERB,I)
       ENDDO
       SQSOINDEXB=SQSOINDEXB_FROM_ORDERS(SQORDERS)
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%sborn_sf.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%sborn_sf.f	2017-07-31 08:48:20 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fks_loonly/%SubProcesses%P0_uxu_ttx%sborn_sf.f	2020-07-04 19:03:17 +0000
@@ -64,8 +64,8 @@
 
         WGT = WGT_COL * G**2
 C       update the amp_split_soft, which is summed in sbornsoft
-        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) = DBLE(AMP_SPLIT_CNT(1:AMP_SPL
-     $IT_SIZE,1,QCD_POS)) * G**2
+        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) =
+     $    DBLE(AMP_SPLIT_CNT(1:AMP_SPLIT_SIZE,1,QCD_POS)) * G**2
 
       ELSE IF (NEED_CHARGE_LINKS) THEN
         CHARGEPROD = CHARGES_BORN(M) * CHARGES_BORN(N)
@@ -76,8 +76,9 @@
         IF (M.EQ.N) CHARGEPROD = CHARGEPROD / 2D0
         WGT = DBLE(ANS_CNT(1, QED_POS)) * CHARGEPROD * DBLE(GAL(1))**2
 C       update the amp_split_soft, which is summed in sbornsoft
-        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) = DBLE(AMP_SPLIT_CNT(1:AMP_SPL
-     $IT_SIZE,1,QED_POS)) * CHARGEPROD * DBLE(GAL(1))**2
+        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) =
+     $    DBLE(AMP_SPLIT_CNT(1:AMP_SPLIT_SIZE,1,QED_POS)) * CHARGEPROD
+     $    * DBLE(GAL(1))**2
       ENDIF
 
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_001.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_001.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_001.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_002.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_002.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_002.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_003.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_003.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_003.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_004.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_004.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_004.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_005.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_005.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_005.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_006.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_006.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_006.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_007.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_007.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_007.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_008.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_008.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_008.f	2020-07-04 19:03:17 +0000
@@ -69,8 +69,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%born.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%born.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%born.f	2020-07-04 19:03:17 +0000
@@ -131,8 +131,8 @@
           DO J = 1, NSPLITORDERS
             AMP_ORDERS(J) = GETORDPOWFROMINDEX_B(J, I)
           ENDDO
-          IF(ABS(ANS(1,I)).GT.MAX_VAL*TINY) AMP_SPLIT(ORDERS_TO_AMP_SPL
-     $IT_POS(AMP_ORDERS)) = ANS(1,I)
+          IF(ABS(ANS(1,I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS)) = ANS(1,I)
         ENDIF
       ENDDO
 C     this is to avoid fake non-zero contributions 
@@ -169,10 +169,12 @@
               IF (K.EQ.J) AMP_ORDERS(K) = AMP_ORDERS(K) + 2
             ENDDO
 C           this is to avoid fake non-zero contributions 
-            IF (ABS(ANS(1,I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_
-     $AMP_SPLIT_POS(AMP_ORDERS),1,J) = ANS(1,I)
-            IF (ABS(ANS(2,I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_
-     $AMP_SPLIT_POS(AMP_ORDERS),2,J) = ANS(2,I)
+            IF (ABS(ANS(1,I)).GT.MAX_VAL*TINY)
+     $        AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1,J) =
+     $        ANS(1,I)
+            IF (ABS(ANS(2,I)).GT.MAX_VAL*TINY)
+     $        AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),2,J) =
+     $        ANS(2,I)
           ENDIF
         ENDDO
 C       this is to avoid fake non-zero contributions 
@@ -599,8 +601,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSPLITORDERS
-        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)+AMPSPLITORDERS(AMPORDER
-     $B,I)
+        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)
+     $   +AMPSPLITORDERS(AMPORDERB,I)
       ENDDO
       SQSOINDEXB=SQSOINDEXB_FROM_ORDERS(SQORDERS)
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_1.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_1.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_1.f	2020-07-04 19:03:17 +0000
@@ -84,8 +84,8 @@
           DO J = 1, NSPLITORDERS
             AMP_ORDERS(J) = GETORDPOWFROMINDEX1(J, I)
           ENDDO
-          IF (ABS(ANS(I)).GT.ANS_MAX*TINY) AMP_SPLIT(ORDERS_TO_AMP_SPLI
-     $T_POS(AMP_ORDERS)) = ANS(I)
+          IF (ABS(ANS(I)).GT.ANS_MAX*TINY)
+     $      AMP_SPLIT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS)) = ANS(I)
         ENDIF
       ENDDO
 
@@ -420,8 +420,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSPLITORDERS
-        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)+AMPSPLITORDERS(AMPORDER
-     $B,I)
+        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)
+     $   +AMPSPLITORDERS(AMPORDERB,I)
       ENDDO
       SQSOINDEX1=SQSOINDEX_FROM_ORDERS1(SQORDERS)
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_2.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_2.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_2.f	2020-07-04 19:03:17 +0000
@@ -87,8 +87,8 @@
           DO J = 1, NSPLITORDERS
             AMP_ORDERS(J) = GETORDPOWFROMINDEX2(J, I)
           ENDDO
-          IF (ABS(ANS(I)).GT.ANS_MAX*TINY) AMP_SPLIT(ORDERS_TO_AMP_SPLI
-     $T_POS(AMP_ORDERS)) = ANS(I)
+          IF (ABS(ANS(I)).GT.ANS_MAX*TINY)
+     $      AMP_SPLIT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS)) = ANS(I)
         ENDIF
       ENDDO
 
@@ -384,8 +384,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSPLITORDERS
-        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)+AMPSPLITORDERS(AMPORDER
-     $B,I)
+        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)
+     $   +AMPSPLITORDERS(AMPORDERB,I)
       ENDDO
       SQSOINDEX2=SQSOINDEX_FROM_ORDERS2(SQORDERS)
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_3.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_3.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_3.f	2020-07-04 19:03:17 +0000
@@ -87,8 +87,8 @@
           DO J = 1, NSPLITORDERS
             AMP_ORDERS(J) = GETORDPOWFROMINDEX3(J, I)
           ENDDO
-          IF (ABS(ANS(I)).GT.ANS_MAX*TINY) AMP_SPLIT(ORDERS_TO_AMP_SPLI
-     $T_POS(AMP_ORDERS)) = ANS(I)
+          IF (ABS(ANS(I)).GT.ANS_MAX*TINY)
+     $      AMP_SPLIT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS)) = ANS(I)
         ENDIF
       ENDDO
 
@@ -384,8 +384,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSPLITORDERS
-        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)+AMPSPLITORDERS(AMPORDER
-     $B,I)
+        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)
+     $   +AMPSPLITORDERS(AMPORDERB,I)
       ENDDO
       SQSOINDEX3=SQSOINDEX_FROM_ORDERS3(SQORDERS)
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_4.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_4.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_4.f	2020-07-04 19:03:17 +0000
@@ -87,8 +87,8 @@
           DO J = 1, NSPLITORDERS
             AMP_ORDERS(J) = GETORDPOWFROMINDEX4(J, I)
           ENDDO
-          IF (ABS(ANS(I)).GT.ANS_MAX*TINY) AMP_SPLIT(ORDERS_TO_AMP_SPLI
-     $T_POS(AMP_ORDERS)) = ANS(I)
+          IF (ABS(ANS(I)).GT.ANS_MAX*TINY)
+     $      AMP_SPLIT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS)) = ANS(I)
         ENDIF
       ENDDO
 
@@ -384,8 +384,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSPLITORDERS
-        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)+AMPSPLITORDERS(AMPORDER
-     $B,I)
+        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)
+     $   +AMPSPLITORDERS(AMPORDERB,I)
       ENDDO
       SQSOINDEX4=SQSOINDEX_FROM_ORDERS4(SQORDERS)
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_5.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_5.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%matrix_5.f	2020-07-04 19:03:17 +0000
@@ -87,8 +87,8 @@
           DO J = 1, NSPLITORDERS
             AMP_ORDERS(J) = GETORDPOWFROMINDEX5(J, I)
           ENDDO
-          IF (ABS(ANS(I)).GT.ANS_MAX*TINY) AMP_SPLIT(ORDERS_TO_AMP_SPLI
-     $T_POS(AMP_ORDERS)) = ANS(I)
+          IF (ABS(ANS(I)).GT.ANS_MAX*TINY)
+     $      AMP_SPLIT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS)) = ANS(I)
         ENDIF
       ENDDO
 
@@ -384,8 +384,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSPLITORDERS
-        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)+AMPSPLITORDERS(AMPORDER
-     $B,I)
+        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)
+     $   +AMPSPLITORDERS(AMPORDERB,I)
       ENDDO
       SQSOINDEX5=SQSOINDEX_FROM_ORDERS5(SQORDERS)
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%sborn_sf.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%sborn_sf.f	2017-07-31 08:48:20 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%sborn_sf.f	2020-07-04 19:03:17 +0000
@@ -64,8 +64,8 @@
 
         WGT = WGT_COL * G**2
 C       update the amp_split_soft, which is summed in sbornsoft
-        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) = DBLE(AMP_SPLIT_CNT(1:AMP_SPL
-     $IT_SIZE,1,QCD_POS)) * G**2
+        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) =
+     $    DBLE(AMP_SPLIT_CNT(1:AMP_SPLIT_SIZE,1,QCD_POS)) * G**2
 
       ELSE IF (NEED_CHARGE_LINKS) THEN
         CHARGEPROD = CHARGES_BORN(M) * CHARGES_BORN(N)
@@ -76,8 +76,9 @@
         IF (M.EQ.N) CHARGEPROD = CHARGEPROD / 2D0
         WGT = DBLE(ANS_CNT(1, QED_POS)) * CHARGEPROD * DBLE(GAL(1))**2
 C       update the amp_split_soft, which is summed in sbornsoft
-        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) = DBLE(AMP_SPLIT_CNT(1:AMP_SPL
-     $IT_SIZE,1,QED_POS)) * CHARGEPROD * DBLE(GAL(1))**2
+        AMP_SPLIT_SOFT(1:AMP_SPLIT_SIZE) =
+     $    DBLE(AMP_SPLIT_CNT(1:AMP_SPLIT_SIZE,1,QED_POS)) * CHARGEPROD
+     $    * DBLE(GAL(1))**2
       ENDIF
 
       RETURN

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%b_sf_001.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%b_sf_001.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%b_sf_001.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%b_sf_002.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%b_sf_002.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%b_sf_002.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
-          IF(ABS(ANS(I)).GT.MAX_VAL*TINY) AMP_SPLIT_CNT(ORDERS_TO_AMP_S
-     $PLIT_POS(AMP_ORDERS),1,QCD_POS) = ANS(I)
+          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
+     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
+     $     ,QCD_POS) = ANS(I)
         ENDIF
       ENDDO
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%b_sf_003.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%b_sf_003.f	2020-03-11 09:28:14 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%b_sf_003.f	2020-07-04 19:03:17 +0000
@@ -72,8 +72,9 @@
             IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
           ENDDO
             !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)