Merge lp:~mg5core1/mg5amcnlo/2.7.3 into lp:mg5amcnlo/lts

[Olivier Mattelaer]

I think that this is time to merge this. (and then update one by one our four release with those fixes).

[Rikkert Frederix (frederix)]

Wait... we have four different releases of our code? Shouldn't we try to get those in line? B.t.w., what's the 4th release: there is the 2.7.2, 3.0.1 and python3, right? Note that PY8meetsMG5aMC_release has not been released: that is still very much a development branch that is being made ready to be merged with 2.7.2 or 2.7.3.

Best,
Rikkert

[Olivier Mattelaer (olivier-mattelaer)]

The four versions are

2.7.2
2.7.2.py3
3.0.2
3.0.2.py3

I guess that we should soon be able to drop the non ".py3" branches (Remember that they do support python2.7, but that's a discussion for our meeting).

Cheers,

Olivier

[marco zaro (marco-zaro)]

Hi Olivier,
the py3 versions can be still run with python2?
If so, I agree we can drop the others (we should remember to upgrade all the various plugins to py3)…

Cheers,

Marco

> On 8 Jun 2020, at 18:03, Olivier Mattelaer <email address hidden> wrote:
>
> The four versions are
>
> 2.7.2
> 2.7.2.py3
> 3.0.2
> 3.0.2.py3
>
> I guess that we should soon be able to drop the non ".py3" branches (Remember that they do support python2.7, but that's a discussion for our meeting).
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~mg5core1/mg5amcnlo/2.7.3/+merge/385288
> Your team mg5core1 is subscribed to branch lp:~mg5core1/mg5amcnlo/2.7.3.

[Olivier Mattelaer (olivier-mattelaer)]

Yes,

the releases:
>> 2.7.2 and 3.02
can run with python 2.6 and 2.7

while
>> 2.7.2.py3 and 3.0.2.py3

can run with python 2.7, 3.7 and 3.8

Cheers,

Olivier

> On 8 Jun 2020, at 18:09, marco zaro <email address hidden> wrote:
>
> Hi Olivier,
> the py3 versions can be still run with python2?
> If so, I agree we can drop the others (we should remember to upgrade all the various plugins to py3)…
>
> Cheers,
>
> Marco
>
>> On 8 Jun 2020, at 18:03, Olivier Mattelaer <email address hidden> wrote:
>>
>> The four versions are
>>
>> 2.7.2
>> 2.7.2.py3
>> 3.0.2
>> 3.0.2.py3
>>
>> I guess that we should soon be able to drop the non ".py3" branches (Remember that they do support python2.7, but that's a discussion for our meeting).
>>
>> Cheers,
>>
>> Olivier
>> --
>> https://code.launchpad.net/~mg5core1/mg5amcnlo/2.7.3/+merge/385288
>> Your team mg5core1 is subscribed to branch lp:~mg5core1/mg5amcnlo/2.7.3.
>
>
> --
> https://code.launchpad.net/~mg5core1/mg5amcnlo/2.7.3/+merge/385288
> Your team mg5core1 is subscribed to branch lp:~mg5core1/mg5amcnlo/2.7.3.

[lp:~mg5core1/mg5amcnlo/2.7.3]

331. By olivier-mattelaer on 2020-06-08

fix numerical stability issue for T-channel

332. By olivier-mattelaer on 2020-06-09

add the install option for the MadSTR plugin

333. By olivier-mattelaer on 2020-06-11

allow new function for standalone f2py mode

334. By olivier-mattelaer on 2020-06-12

fix an issue with alphas runninng

335. By olivier-mattelaer on 2020-06-12

fix gcc version issue

336. By olivier-mattelaer on 2020-06-15

fix a bug introduce in the wrong merging of feature for increase numerical stability

337. By olivier-mattelaer on 2020-06-15

change default lhef version for merge.pl (from 1 to 3)

338. By olivier-mattelaer on 2020-06-16

change a comment referencing a wrong parameter

339. By olivier-mattelaer on 2020-06-16

add global_order_coupling option for MadSPin

340. By olivier-mattelaer on 2020-06-17

fix an issue for auto width when the path contains a equal sign

341. By olivier-mattelaer on 2020-06-17

merge with rwgt_fix to support multi-thread reweighting

342. By olivier-mattelaer on 2020-06-17

remove a debug raise Exception that should not have been committed

[Rikkert Frederix (frederix)]

Hi Olivier,

I just learned that python 2.7 was released in July 2010 -- that's 10 years ago! I think it's okay to assume that users have upgraded from 2.6 to 2.7 (if they do not want to use 3.X).

Cheers,
Rikkert

[lp:~mg5core1/mg5amcnlo/2.7.3]

343. By olivier-mattelaer on 2020-06-18

fix an issue for gridpack in readonly

344. By olivier-mattelaer on 2020-06-19

fixing a test for one-loop reweighting (change rwgt model without card now use default model param)

345. By olivier-mattelaer on 2020-06-20

1. Fixing some iotest (changing of the fortranwriter)
2. changing diagram filter function at LO
3. fixing some issue with re-weighting

346. By olivier-mattelaer on 2020-06-21

update tests

347. By olivier-mattelaer on 2020-06-21

update version/update note for release

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2020-01-08 14:13:32 +0000
+++ MadSpin/decay.py	2020-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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	2017-11-17 16:07:39 +0000
+++ Template/NLO/SubProcesses/cuts.f	2020-06-21 12:17:46 +0000
@@ -553,10 +553,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)

=== 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-06-21 12:17:46 +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-03-17 07:12:13 +0000
+++ UpdateNotes.txt	2020-06-21 12:17:46 +0000
@@ -1,5 +1,24 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
+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
+
 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-03-17 07:12:13 +0000
+++ VERSION	2020-06-21 12:17:46 +0000
@@ -1,5 +1,5 @@
-version = 2.7.2
-date = 2020-03-17
+version = 2.7.3
+date = 2020-06-21
 
 
 

=== modified file 'aloha/create_aloha.py'
--- aloha/create_aloha.py	2020-01-09 16:39:38 +0000
+++ aloha/create_aloha.py	2020-06-21 12:17:46 +0000
@@ -443,7 +443,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-06-21 12:17:46 +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-06-21 12:17:46 +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-04 20:39:24 +0000
+++ madgraph/core/base_objects.py	2020-06-21 12:17:46 +0000
@@ -3484,18 +3484,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	2019-05-13 20:41:22 +0000
+++ madgraph/core/diagram_generation.py	2020-06-21 12:17:46 +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-06-21 12:17:46 +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-07 20:56:41 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2020-06-21 12:17:46 +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:
@@ -2785,7 +2791,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-04 09:57:24 +0000
+++ madgraph/interface/common_run_interface.py	2020-06-21 12:17:46 +0000
@@ -1036,7 +1036,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):
@@ -2836,8 +2840,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:]
@@ -2852,7 +2856,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 "+
@@ -3328,7 +3332,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:
@@ -4098,6 +4102,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)):
@@ -4309,7 +4325,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
                 
@@ -4495,7 +4511,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:
@@ -4773,10 +4794,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 []
     
@@ -5241,7 +5264,7 @@
     def do_set(self, line):
         """ edit the value of one parameter in the card"""
         
-        
+
         args = self.split_arg(line)
         
         
@@ -5504,7 +5527,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
@@ -5965,8 +5988,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):
@@ -6092,7 +6116,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
@@ -6477,6 +6502,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('')
@@ -6622,7 +6648,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-02-25 14:05:36 +0000
+++ madgraph/interface/extended_cmd.py	2020-06-21 12:17:46 +0000
@@ -2586,6 +2586,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:
@@ -2601,6 +2603,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-06-21 12:17:46 +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-06 20:47:59 +0000
+++ madgraph/interface/madgraph_interface.py	2020-06-21 12:17:46 +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")
@@ -926,7 +931,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
@@ -1075,7 +1082,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')
@@ -1670,7 +1677,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
@@ -2832,7 +2839,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',
@@ -3869,6 +3876,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':
@@ -5615,7 +5631,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)
@@ -5985,7 +6001,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'],
@@ -6003,7 +6019,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']
@@ -6129,7 +6146,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
@@ -6352,6 +6372,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))
@@ -6400,18 +6421,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-06-21 12:17:46 +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-01-15 18:23:15 +0000
+++ madgraph/interface/reweight_interface.py	2020-06-21 12:17:46 +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):
             try:
@@ -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-03-10 15:07:48 +0000
+++ madgraph/iolibs/export_v4.py	2020-06-21 12:17:46 +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)
@@ -2035,7 +2034,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'))
@@ -2043,7 +2042,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,
@@ -2100,6 +2099,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
@@ -2154,20 +2181,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-06-21 12:17:46 +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/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-06-21 12:17:46 +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-06-05 11:03:48 +0000
+++ madgraph/loop/loop_exporters.py	2020-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-03-08 20:49:57 +0000
+++ madgraph/various/banner.py	2020-06-21 12:17:46 +0000
@@ -3091,7 +3091,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)")
@@ -3195,8 +3196,6 @@
                 if self['mmjj'] > self['xqcut']:
                     logger.warning('mmjj > xqcut (and auto_ptj_mjj = F). MMJJ set to 0')
                     self['mmjj'] = 0.0 
-
-
     
         # check validity of the pdf set
         if self['pdlabel'] == 'lhapdf':
@@ -4067,7 +4066,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/lhe_parser.py'
--- madgraph/various/lhe_parser.py	2020-01-09 23:14:56 +0000
+++ madgraph/various/lhe_parser.py	2020-06-21 12:17:46 +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-06-21 12:17:46 +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-06-21 12:17:46 +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-06-20 08:18:36 +0000
+++ models/model_reader.py	2020-06-21 12:17:46 +0000
@@ -181,6 +181,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-06-21 12:17:46 +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_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%COLLIER_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%COLLIER_interface.f	2016-08-07 07:16:07 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%COLLIER_interface.f	2020-06-21 12:17:46 +0000
@@ -164,8 +164,9 @@
             IF (NCALLS_IN_CACHE(I).NE.-1.AND..NOT. ( NCOLLIERCALLS(I)
      $       .EQ.NCALLS_IN_CACHE(I).OR. ( CTMODERUN.EQ.
      $       -1.AND..NOT.COLLIERUSEINTERNALSTABILITYTEST.AND.NROTATIONS
-     $_DP.GT.0.AND.MOD(NCALLS_IN_CACHE(I),2).EQ.0.AND.(NCALLS_IN_CACHE(
-     $I)/2).EQ.NCOLLIERCALLS(I) ) ) ) THEN
+     $_DP.GT.0.AND.MOD(NCALLS_IN_CACHE(I),2)
+     $       .EQ.0.AND.(NCALLS_IN_CACHE(I)/2).EQ.NCOLLIERCALLS(I) ) ) )
+     $        THEN
               WRITE(*,*) 'WARNING: A consistency check in MadLoop'
      $         //' failed and, for safety, forced MadLoop to'
      $         //' reinitialize the global cache of COLLIER. Report'
@@ -256,8 +257,8 @@
 C               The expression below is like taking the absolute value
 C                when summing errors linearly
 C               TIRCOEFSERRORS(I,K)=(TIR_COEFS_ERROR(CURR_RANK,ID)/MaxC
-C               oefForRank(CURR_RANK))*DCMPLX( ABS(DBLE(TIRCOEFS(I,K)))
-C               ,ABS(DIMAG(TIRCOEFS(I,K))) )
+C               oefForRank(CURR_RANK))*DCMPLX(
+C                ABS(DBLE(TIRCOEFS(I,K))),ABS(DIMAG(TIRCOEFS(I,K))) )
 C               But empirically, I observed that retaining the
 C                original complex phase leads to slightly more
 C                accurate estimates
@@ -479,8 +480,8 @@
 C             The expression below is like taking the absolute value
 C              when summing errors linearly
 C             TIRCOEFSERRORS(I,K)=(TNtenerr(CURR_RANK)/MaxCoefForRank(C
-C             URR_RANK))*DCMPLX( ABS(DBLE(TIRCOEFS(I,K))),ABS(DIMAG(TIR
-C             COEFS(I,K))) )			
+C             URR_RANK))*DCMPLX(
+C              ABS(DBLE(TIRCOEFS(I,K))),ABS(DIMAG(TIRCOEFS(I,K))) )			
 C             But empirically, I observed that retaining the original
 C              complex phase leads to slightly more accurate estimates
               TIRCOEFSERRORS(I,K)=(TNTENERR(CURR_RANK)
@@ -592,8 +593,8 @@
 
 C     The common blocks below are to retrieve the necessary
 C      information about
-C     MadLoop running mode and store it in the sCOLLIER_CACHE_RELEVANT_
-C     PARAMS common block.
+C     MadLoop running mode and store it in the
+C      sCOLLIER_CACHE_RELEVANT_PARAMS common block.
 
       INCLUDE 'MadLoopParams.inc'
       INCLUDE 'unique_id.inc'
@@ -662,8 +663,8 @@
         USERHEL_BU         = USERHEL
         SQSO_TARGET_BU     = SQSO_TARGET
         COLLIERMODE_BU     = COLLIERMODE
-        COLLIERUSEINTERNALSTABILITYTEST_BU = COLLIERUSEINTERNALSTABILIT
-     $YTEST
+        COLLIERUSEINTERNALSTABILITYTEST_BU =
+     $    COLLIERUSEINTERNALSTABILITYTEST
         CTMODERUN_BU       = CTMODERUN
         CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*4+1)
         COLLIER_CACHE_ACTIVE = 1
@@ -710,8 +711,8 @@
         USERHEL_BU         = USERHEL
         SQSO_TARGET_BU     = SQSO_TARGET
         COLLIERMODE_BU     = COLLIERMODE
-        COLLIERUSEINTERNALSTABILITYTEST_BU = COLLIERUSEINTERNALSTABILIT
-     $YTEST
+        COLLIERUSEINTERNALSTABILITYTEST_BU =
+     $    COLLIERUSEINTERNALSTABILITYTEST
         CTMODERUN_BU       = CTMODERUN
         IF (COLLIERGLOBALCACHE.EQ.-1) THEN
           CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS,MAXNEXTERNAL)

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f	2016-08-07 07:16:07 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f	2020-06-21 12:17:46 +0000
@@ -663,8 +663,8 @@
 C     Determine it uses qp or not
       DOING_QP = (CTMODE.GE.4)
 
-      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(SQUAREDSOIND
-     $EX,LOOPNUM)) THEN
+      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)
+     $ .OR.GOODAMP(SQUAREDSOINDEX,LOOPNUM)) THEN
         WE(1)=W1
         WE(2)=W2
         WE(3)=W3

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f	2018-11-08 22:24:49 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f	2020-06-21 12:17:46 +0000
@@ -54,8 +54,8 @@
       INTEGER LAST_LIB_USED
       DATA LAST_LIB_USED/-1/
 
-      COMPLEX*16 TIRCOEFS(0:LOOPMAXCOEFS-1,3),TIRCOEFSERRORS(0:LOOPMAXC
-     $OEFS-1,3)
+      COMPLEX*16 TIRCOEFS(0:LOOPMAXCOEFS-1,3)
+     $ ,TIRCOEFSERRORS(0:LOOPMAXCOEFS-1,3)
       COMPLEX*16 PJCOEFS(0:LOOPMAXCOEFS-1,3)
 C     
 C     EXTERNAL FUNCTIONS
@@ -398,8 +398,9 @@
 C     
       INCLUDE 'MadLoopParams.inc'
       INCLUDE 'process_info.inc'
-C     Change the list 'LOOPLIBS_QPAVAILABLE' in loop_matrix_standalone.
-C     inc to change the list of QPTools availables
+C     Change the list 'LOOPLIBS_QPAVAILABLE' in
+C      loop_matrix_standalone.inc to change the list of QPTools
+C      availables
       LOGICAL QP_TOOLS_AVAILABLE
       INTEGER INDEX_QP_TOOLS(QP_NLOOPLIB+1)
       COMMON/LOOP_TOOLS/QP_TOOLS_AVAILABLE,INDEX_QP_TOOLS

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f	2018-05-18 12:22:47 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f	2020-06-21 12:17:46 +0000
@@ -461,8 +461,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSO
-        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)+AMPSPLITORDERS(ORDERI
-     $NDEXB,I)
+        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)
+     $   +AMPSPLITORDERS(ORDERINDEXB,I)
       ENDDO
       SQSOINDEX=SOINDEX_FOR_SQUARED_ORDERS(SQORDERS)
       END
@@ -556,8 +556,8 @@
         RETURN
       ENDIF
 
-      WRITE(*,*) 'ERROR:: Stopping function GET_SQUARED_ORDERS_FOR_SOIN'
-     $ //'DEX'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' GET_SQUARED_ORDERS_FOR_SOINDEX'
       WRITE(*,*) 'Could not find squared orders index ',SOINDEX
       STOP
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa_born_splitOrders.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa_born_splitOrders.f	2016-03-31 10:53:12 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa_born_splitOrders.f	2020-06-21 12:17:46 +0000
@@ -170,8 +170,8 @@
 C     CALL SMATRIX(P,MATELEM)
 C     
 C     write (*,*) "-------------------------------------------------"
-C     write (*,*) "Matrix element = ", MATELEM, " GeV^",-(2*nexternal-8)
-C     	
+C     write (*,*) "Matrix element = ", MATELEM, "
+C      GeV^",-(2*nexternal-8)	
 C     write (*,*) "-------------------------------------------------"
 
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f	2018-11-08 22:24:49 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f	2020-06-21 12:17:46 +0000
@@ -78,8 +78,8 @@
       INTEGER NSQUAREDSOP1
       PARAMETER (NSQUAREDSOP1=NSQUAREDSO+1)
 C     The total number of loop reduction libraries
-C     At present, there are only CutTools,PJFry++,IREGI,Golem95,Samurai
-C     , Ninja and COLLIER
+C     At present, there are only
+C      CutTools,PJFry++,IREGI,Golem95,Samurai, Ninja and COLLIER
       INTEGER NLOOPLIB
       PARAMETER (NLOOPLIB=7)
 C     Only CutTools or possibly Ninja (if installed with qp support)
@@ -404,14 +404,15 @@
 C      reading the parameters
       LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
      $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPU
-     $TATION_CHOICE
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE,
+     $  COLLIER_IR_POLE_COMPUTATION_CHOICE
       DATA  FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION
      $ ,FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION/.FALSE.,.FALSE./
-      COMMON/COLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_
-     $COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_
-     $CHOICE
+      COMMON/COLLIERPOLESFORCEDCHOICE
+     $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
+     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE
+     $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE
 
 C     This variable controls the general initialization which is
 C      *common* between all MadLoop SubProcesses.
@@ -880,8 +881,8 @@
      $ +P(2,2))**2-(P(3,1)+P(3,2))**2))
 
       CTCALL_REQ_SO_DONE=.FALSE.
-      FILTER_SO = (.NOT.CHECKPHASE).AND.HELDOUBLECHECKED.AND.(SQSO_TARG
-     $ET.NE.-1)
+      FILTER_SO = (.NOT.CHECKPHASE)
+     $ .AND.HELDOUBLECHECKED.AND.(SQSO_TARGET.NE.-1)
 
       DO I=1,NLOOPGROUPS
         DO J=0,LOOPMAXCOEFS-1
@@ -920,13 +921,14 @@
       ENDIF
 
       DO H=1,NCOMB
-        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1).AND.(CHECKPHASE.OR.(
-     $.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H).GT.-HELOFFSET.AND.GOODHEL(H)
-     $   .NE.0)))) THEN
+        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1)
+     $   .AND.(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H)
+     $   .GT.-HELOFFSET.AND.GOODHEL(H).NE.0)))) THEN
 
 C         Handle the possible requirement of specific polarizations
-          IF ((.NOT.CHECKPHASE).AND.HELDOUBLECHECKED.AND.POLARIZATIONS(
-     $0,0).EQ.0.AND.(.NOT.IS_HEL_SELECTED(H))) THEN
+          IF ((.NOT.CHECKPHASE)
+     $     .AND.HELDOUBLECHECKED.AND.POLARIZATIONS(0,0)
+     $     .EQ.0.AND.(.NOT.IS_HEL_SELECTED(H))) THEN
             CYCLE
           ENDIF
 
@@ -1164,8 +1166,8 @@
 C             Of course if it is one, then we do not need to do
 C              anything (because with HELINITSTARTOVER=.FALSE. we only
 C              support exactly identical Hels.)
-              IF(GOODHEL(HELPICKED).GT.-HELOFFSET.AND.GOODHEL(HELPICKED)
-     $         .NE.1) THEN
+              IF(GOODHEL(HELPICKED).GT.
+     $         -HELOFFSET.AND.GOODHEL(HELPICKED).NE.1) THEN
                 NEWHELREF=-1
                 DO H=1,NCOMB
                   IF (GOODHEL(H).EQ.(-HELOFFSET-HELPICKED)) THEN
@@ -1359,16 +1361,18 @@
 
         CTMODE=BASIC_CT_MODE
 
-        IF(.NOT.EVAL_DONE(3).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
-     $.1).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.1)) ) THEN
+        IF(.NOT.EVAL_DONE(3).AND.
+     $    ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE.1)
+     $   .OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.1)) ) THEN
           EVAL_DONE(3)=.TRUE.
           CALL ROTATE_PS(PS,P,1)
           IF (DOING_QP_EVALS) CALL MP_ROTATE_PS(MP_PS,MP_P,1)
           GOTO 200
         ENDIF
 
-        IF(.NOT.EVAL_DONE(4).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
-     $.2).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.2)) ) THEN
+        IF(.NOT.EVAL_DONE(4).AND.
+     $    ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE.2)
+     $   .OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.2)) ) THEN
           EVAL_DONE(4)=.TRUE.
           CALL ROTATE_PS(PS,P,2)
           IF (DOING_QP_EVALS) CALL MP_ROTATE_PS(MP_PS,MP_P,2)
@@ -1924,8 +1928,8 @@
 C       When using COLLIER with the internal stability test, the first
 C        evaluation is typically more reliable so we do not want to
 C        use the average but rather the first evaluation.
-        IF (MLREDUCTIONLIB(I_LIB).EQ.7.AND.COLLIERUSEINTERNALSTABILITYT
-     $EST) THEN
+        IF (MLREDUCTIONLIB(I_LIB)
+     $   .EQ.7.AND.COLLIERUSEINTERNALSTABILITYTEST) THEN
           DO I=1,3
             ESTIMATE(I,K) = FULLLIST(I,K,1)
           ENDDO
@@ -2070,8 +2074,8 @@
  1009   CONTINUE
       ENDDO
 
-      WRITE(*,*) 'ERROR:: Stopping function ML5SOINDEX_FOR_SQUARED_ORDE'
-     $ //'RS'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' ML5SOINDEX_FOR_SQUARED_ORDERS'
       WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO)
       STOP
 
@@ -2170,8 +2174,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSO
-        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)+AMPSPLITORDERS(ORDERI
-     $NDEXB,I)
+        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)
+     $   +AMPSPLITORDERS(ORDERINDEXB,I)
       ENDDO
       ML5SQSOINDEX=ML5SOINDEX_FOR_SQUARED_ORDERS(SQORDERS)
       END
@@ -2209,8 +2213,8 @@
         RETURN
       ENDIF
 
-      WRITE(*,*) 'ERROR:: Stopping function ML5GET_SQUARED_ORDERS_FOR_S'
-     $ //'OINDEX'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' ML5GET_SQUARED_ORDERS_FOR_SOINDEX'
       WRITE(*,*) 'Could not find squared orders index ',SOINDEX
       STOP
 
@@ -2249,8 +2253,8 @@
         RETURN
       ENDIF
 
-      WRITE(*,*) 'ERROR:: Stopping function ML5GET_ORDERS_FOR_AMPSOINDE'
-     $ //'X'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' ML5GET_ORDERS_FOR_AMPSOINDEX'
       WRITE(*,*) 'Could not find amplitude split orders index ',SOINDEX
       STOP
 
@@ -2315,12 +2319,13 @@
 
       LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
      $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPU
-     $TATION_CHOICE
-      COMMON/COLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_
-     $COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_
-     $CHOICE
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE,
+     $  COLLIER_IR_POLE_COMPUTATION_CHOICE
+      COMMON/COLLIERPOLESFORCEDCHOICE
+     $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
+     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE
+     $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE
 
       COLLIERCOMPUTEUVPOLES                        = ONOFF
 C     This is just so that if we read the param again, we don't
@@ -2342,12 +2347,13 @@
 
       LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
      $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPU
-     $TATION_CHOICE
-      COMMON/COLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_
-     $COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_
-     $CHOICE
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE,
+     $  COLLIER_IR_POLE_COMPUTATION_CHOICE
+      COMMON/COLLIERPOLESFORCEDCHOICE
+     $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
+     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE
+     $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE
 
       COLLIERCOMPUTEIRPOLES         = ONOFF
 C     This is just so that if we read the param again, we don't

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%mp_compute_loop_coefs.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%mp_compute_loop_coefs.f	2017-07-07 10:47:20 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%mp_compute_loop_coefs.f	2020-06-21 12:17:46 +0000
@@ -268,13 +268,14 @@
 
       DPW_COPIED = .FALSE.
       DO H=1,NCOMB
-        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1).AND.(CHECKPHASE.OR.(
-     $.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H).GT.-HELOFFSET.AND.GOODHEL(H)
-     $   .NE.0)))) THEN
+        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1)
+     $   .AND.(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H)
+     $   .GT.-HELOFFSET.AND.GOODHEL(H).NE.0)))) THEN
 
 C         Handle the possible requirement of specific polarizations
-          IF ((.NOT.CHECKPHASE).AND.HELDOUBLECHECKED.AND.POLARIZATIONS(
-     $0,0).EQ.0.AND.(.NOT.IS_HEL_SELECTED(H))) THEN
+          IF ((.NOT.CHECKPHASE)
+     $     .AND.HELDOUBLECHECKED.AND.POLARIZATIONS(0,0)
+     $     .EQ.0.AND.(.NOT.IS_HEL_SELECTED(H))) THEN
             CYCLE
           ENDIF
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%polynomial.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%polynomial.f	2016-04-06 18:45:52 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%polynomial.f	2020-06-21 12:17:46 +0000
@@ -113,8 +113,8 @@
         CFTOT=CMPLX(CF_N(COLOR_ID,I)/(ONE*ABS(CF_D(COLOR_ID,I))),ZERO
      $   ,KIND=8)
         IF(CF_D(COLOR_ID,I).LT.0) CFTOT=CFTOT*IMAG1
-        CONST(ML5SOINDEX_FOR_BORN_AMP(I))=CONST(ML5SOINDEX_FOR_BORN_AMP
-     $(I))+CFTOT*CONJG(AMP(I))
+        CONST(ML5SOINDEX_FOR_BORN_AMP(I))
+     $   =CONST(ML5SOINDEX_FOR_BORN_AMP(I))+CFTOT*CONJG(AMP(I))
       ENDDO
 
       DO I=1,NAMPSO
@@ -232,8 +232,8 @@
         CFTOT=CMPLX(CF_N(COLOR_ID,I)/(ONE*ABS(CF_D(COLOR_ID,I))),ZERO
      $   ,KIND=16)
         IF(CF_D(COLOR_ID,I).LT.0) CFTOT=CFTOT*IMAG1
-        CONST(ML5SOINDEX_FOR_BORN_AMP(I))=CONST(ML5SOINDEX_FOR_BORN_AMP
-     $(I))+CFTOT*CONJG(AMP(I))
+        CONST(ML5SOINDEX_FOR_BORN_AMP(I))
+     $   =CONST(ML5SOINDEX_FOR_BORN_AMP(I))+CFTOT*CONJG(AMP(I))
       ENDDO
 
       DO I=1,NAMPSO

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%COLLIER_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%COLLIER_interface.f	2016-08-07 07:16:07 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%COLLIER_interface.f	2020-06-21 12:17:46 +0000
@@ -164,8 +164,9 @@
             IF (NCALLS_IN_CACHE(I).NE.-1.AND..NOT. ( NCOLLIERCALLS(I)
      $       .EQ.NCALLS_IN_CACHE(I).OR. ( CTMODERUN.EQ.
      $       -1.AND..NOT.COLLIERUSEINTERNALSTABILITYTEST.AND.NROTATIONS
-     $_DP.GT.0.AND.MOD(NCALLS_IN_CACHE(I),2).EQ.0.AND.(NCALLS_IN_CACHE(
-     $I)/2).EQ.NCOLLIERCALLS(I) ) ) ) THEN
+     $_DP.GT.0.AND.MOD(NCALLS_IN_CACHE(I),2)
+     $       .EQ.0.AND.(NCALLS_IN_CACHE(I)/2).EQ.NCOLLIERCALLS(I) ) ) )
+     $        THEN
               WRITE(*,*) 'WARNING: A consistency check in MadLoop'
      $         //' failed and, for safety, forced MadLoop to'
      $         //' reinitialize the global cache of COLLIER. Report'
@@ -256,8 +257,8 @@
 C               The expression below is like taking the absolute value
 C                when summing errors linearly
 C               TIRCOEFSERRORS(I,K)=(TIR_COEFS_ERROR(CURR_RANK,ID)/MaxC
-C               oefForRank(CURR_RANK))*DCMPLX( ABS(DBLE(TIRCOEFS(I,K)))
-C               ,ABS(DIMAG(TIRCOEFS(I,K))) )
+C               oefForRank(CURR_RANK))*DCMPLX(
+C                ABS(DBLE(TIRCOEFS(I,K))),ABS(DIMAG(TIRCOEFS(I,K))) )
 C               But empirically, I observed that retaining the
 C                original complex phase leads to slightly more
 C                accurate estimates
@@ -479,8 +480,8 @@
 C             The expression below is like taking the absolute value
 C              when summing errors linearly
 C             TIRCOEFSERRORS(I,K)=(TNtenerr(CURR_RANK)/MaxCoefForRank(C
-C             URR_RANK))*DCMPLX( ABS(DBLE(TIRCOEFS(I,K))),ABS(DIMAG(TIR
-C             COEFS(I,K))) )			
+C             URR_RANK))*DCMPLX(
+C              ABS(DBLE(TIRCOEFS(I,K))),ABS(DIMAG(TIRCOEFS(I,K))) )			
 C             But empirically, I observed that retaining the original
 C              complex phase leads to slightly more accurate estimates
               TIRCOEFSERRORS(I,K)=(TNTENERR(CURR_RANK)
@@ -592,8 +593,8 @@
 
 C     The common blocks below are to retrieve the necessary
 C      information about
-C     MadLoop running mode and store it in the sCOLLIER_CACHE_RELEVANT_
-C     PARAMS common block.
+C     MadLoop running mode and store it in the
+C      sCOLLIER_CACHE_RELEVANT_PARAMS common block.
 
       INCLUDE 'MadLoopParams.inc'
       INCLUDE 'unique_id.inc'
@@ -662,8 +663,8 @@
         USERHEL_BU         = USERHEL
         SQSO_TARGET_BU     = SQSO_TARGET
         COLLIERMODE_BU     = COLLIERMODE
-        COLLIERUSEINTERNALSTABILITYTEST_BU = COLLIERUSEINTERNALSTABILIT
-     $YTEST
+        COLLIERUSEINTERNALSTABILITYTEST_BU =
+     $    COLLIERUSEINTERNALSTABILITYTEST
         CTMODERUN_BU       = CTMODERUN
         CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*4+1)
         COLLIER_CACHE_ACTIVE = 1
@@ -710,8 +711,8 @@
         USERHEL_BU         = USERHEL
         SQSO_TARGET_BU     = SQSO_TARGET
         COLLIERMODE_BU     = COLLIERMODE
-        COLLIERUSEINTERNALSTABILITYTEST_BU = COLLIERUSEINTERNALSTABILIT
-     $YTEST
+        COLLIERUSEINTERNALSTABILITYTEST_BU =
+     $    COLLIERUSEINTERNALSTABILITYTEST
         CTMODERUN_BU       = CTMODERUN
         IF (COLLIERGLOBALCACHE.EQ.-1) THEN
           CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS,MAXNEXTERNAL)

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%CT_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%CT_interface.f	2016-08-07 07:16:07 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%CT_interface.f	2020-06-21 12:17:46 +0000
@@ -663,8 +663,8 @@
 C     Determine it uses qp or not
       DOING_QP = (CTMODE.GE.4)
 
-      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(SQUAREDSOIND
-     $EX,LOOPNUM)) THEN
+      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)
+     $ .OR.GOODAMP(SQUAREDSOINDEX,LOOPNUM)) THEN
         WE(1)=W1
         WE(2)=W2
         WE(3)=W3

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%TIR_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%TIR_interface.f	2018-11-08 22:24:49 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%TIR_interface.f	2020-06-21 12:17:46 +0000
@@ -54,8 +54,8 @@
       INTEGER LAST_LIB_USED
       DATA LAST_LIB_USED/-1/
 
-      COMPLEX*16 TIRCOEFS(0:LOOPMAXCOEFS-1,3),TIRCOEFSERRORS(0:LOOPMAXC
-     $OEFS-1,3)
+      COMPLEX*16 TIRCOEFS(0:LOOPMAXCOEFS-1,3)
+     $ ,TIRCOEFSERRORS(0:LOOPMAXCOEFS-1,3)
       COMPLEX*16 PJCOEFS(0:LOOPMAXCOEFS-1,3)
 C     
 C     EXTERNAL FUNCTIONS
@@ -398,8 +398,9 @@
 C     
       INCLUDE 'MadLoopParams.inc'
       INCLUDE 'process_info.inc'
-C     Change the list 'LOOPLIBS_QPAVAILABLE' in loop_matrix_standalone.
-C     inc to change the list of QPTools availables
+C     Change the list 'LOOPLIBS_QPAVAILABLE' in
+C      loop_matrix_standalone.inc to change the list of QPTools
+C      availables
       LOGICAL QP_TOOLS_AVAILABLE
       INTEGER INDEX_QP_TOOLS(QP_NLOOPLIB+1)
       COMMON/LOOP_TOOLS/QP_TOOLS_AVAILABLE,INDEX_QP_TOOLS

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%born_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%born_matrix.f	2018-05-18 12:22:47 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%born_matrix.f	2020-06-21 12:17:46 +0000
@@ -461,8 +461,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSO
-        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)+AMPSPLITORDERS(ORDERI
-     $NDEXB,I)
+        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)
+     $   +AMPSPLITORDERS(ORDERINDEXB,I)
       ENDDO
       SQSOINDEX=SOINDEX_FOR_SQUARED_ORDERS(SQORDERS)
       END
@@ -556,8 +556,8 @@
         RETURN
       ENDIF
 
-      WRITE(*,*) 'ERROR:: Stopping function GET_SQUARED_ORDERS_FOR_SOIN'
-     $ //'DEX'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' GET_SQUARED_ORDERS_FOR_SOINDEX'
       WRITE(*,*) 'Could not find squared orders index ',SOINDEX
       STOP
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%check_sa_born_splitOrders.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%check_sa_born_splitOrders.f	2016-03-31 10:53:12 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%check_sa_born_splitOrders.f	2020-06-21 12:17:46 +0000
@@ -170,8 +170,8 @@
 C     CALL SMATRIX(P,MATELEM)
 C     
 C     write (*,*) "-------------------------------------------------"
-C     write (*,*) "Matrix element = ", MATELEM, " GeV^",-(2*nexternal-8)
-C     	
+C     write (*,*) "Matrix element = ", MATELEM, "
+C      GeV^",-(2*nexternal-8)	
 C     write (*,*) "-------------------------------------------------"
 
       END

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_matrix.f	2018-11-08 22:24:49 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_matrix.f	2020-06-21 12:17:46 +0000
@@ -78,8 +78,8 @@
       INTEGER NSQUAREDSOP1
       PARAMETER (NSQUAREDSOP1=NSQUAREDSO+1)
 C     The total number of loop reduction libraries
-C     At present, there are only CutTools,PJFry++,IREGI,Golem95,Samurai
-C     , Ninja and COLLIER
+C     At present, there are only
+C      CutTools,PJFry++,IREGI,Golem95,Samurai, Ninja and COLLIER
       INTEGER NLOOPLIB
       PARAMETER (NLOOPLIB=7)
 C     Only CutTools or possibly Ninja (if installed with qp support)
@@ -404,14 +404,15 @@
 C      reading the parameters
       LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
      $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPU
-     $TATION_CHOICE
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE,
+     $  COLLIER_IR_POLE_COMPUTATION_CHOICE
       DATA  FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION
      $ ,FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION/.FALSE.,.FALSE./
-      COMMON/COLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_
-     $COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_
-     $CHOICE
+      COMMON/COLLIERPOLESFORCEDCHOICE
+     $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
+     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE
+     $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE
 
 C     This variable controls the general initialization which is
 C      *common* between all MadLoop SubProcesses.
@@ -880,8 +881,8 @@
      $ +P(2,2))**2-(P(3,1)+P(3,2))**2))
 
       CTCALL_REQ_SO_DONE=.FALSE.
-      FILTER_SO = (.NOT.CHECKPHASE).AND.HELDOUBLECHECKED.AND.(SQSO_TARG
-     $ET.NE.-1)
+      FILTER_SO = (.NOT.CHECKPHASE)
+     $ .AND.HELDOUBLECHECKED.AND.(SQSO_TARGET.NE.-1)
 
       DO I=1,NLOOPGROUPS
         DO J=0,LOOPMAXCOEFS-1
@@ -920,13 +921,14 @@
       ENDIF
 
       DO H=1,NCOMB
-        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1).AND.(CHECKPHASE.OR.(
-     $.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H).GT.-HELOFFSET.AND.GOODHEL(H)
-     $   .NE.0)))) THEN
+        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1)
+     $   .AND.(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H)
+     $   .GT.-HELOFFSET.AND.GOODHEL(H).NE.0)))) THEN
 
 C         Handle the possible requirement of specific polarizations
-          IF ((.NOT.CHECKPHASE).AND.HELDOUBLECHECKED.AND.POLARIZATIONS(
-     $0,0).EQ.0.AND.(.NOT.IS_HEL_SELECTED(H))) THEN
+          IF ((.NOT.CHECKPHASE)
+     $     .AND.HELDOUBLECHECKED.AND.POLARIZATIONS(0,0)
+     $     .EQ.0.AND.(.NOT.IS_HEL_SELECTED(H))) THEN
             CYCLE
           ENDIF
 
@@ -1164,8 +1166,8 @@
 C             Of course if it is one, then we do not need to do
 C              anything (because with HELINITSTARTOVER=.FALSE. we only
 C              support exactly identical Hels.)
-              IF(GOODHEL(HELPICKED).GT.-HELOFFSET.AND.GOODHEL(HELPICKED)
-     $         .NE.1) THEN
+              IF(GOODHEL(HELPICKED).GT.
+     $         -HELOFFSET.AND.GOODHEL(HELPICKED).NE.1) THEN
                 NEWHELREF=-1
                 DO H=1,NCOMB
                   IF (GOODHEL(H).EQ.(-HELOFFSET-HELPICKED)) THEN
@@ -1359,16 +1361,18 @@
 
         CTMODE=BASIC_CT_MODE
 
-        IF(.NOT.EVAL_DONE(3).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
-     $.1).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.1)) ) THEN
+        IF(.NOT.EVAL_DONE(3).AND.
+     $    ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE.1)
+     $   .OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.1)) ) THEN
           EVAL_DONE(3)=.TRUE.
           CALL ROTATE_PS(PS,P,1)
           IF (DOING_QP_EVALS) CALL MP_ROTATE_PS(MP_PS,MP_P,1)
           GOTO 200
         ENDIF
 
-        IF(.NOT.EVAL_DONE(4).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
-     $.2).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.2)) ) THEN
+        IF(.NOT.EVAL_DONE(4).AND.
+     $    ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE.2)
+     $   .OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.2)) ) THEN
           EVAL_DONE(4)=.TRUE.
           CALL ROTATE_PS(PS,P,2)
           IF (DOING_QP_EVALS) CALL MP_ROTATE_PS(MP_PS,MP_P,2)
@@ -1924,8 +1928,8 @@
 C       When using COLLIER with the internal stability test, the first
 C        evaluation is typically more reliable so we do not want to
 C        use the average but rather the first evaluation.
-        IF (MLREDUCTIONLIB(I_LIB).EQ.7.AND.COLLIERUSEINTERNALSTABILITYT
-     $EST) THEN
+        IF (MLREDUCTIONLIB(I_LIB)
+     $   .EQ.7.AND.COLLIERUSEINTERNALSTABILITYTEST) THEN
           DO I=1,3
             ESTIMATE(I,K) = FULLLIST(I,K,1)
           ENDDO
@@ -2070,8 +2074,8 @@
  1009   CONTINUE
       ENDDO
 
-      WRITE(*,*) 'ERROR:: Stopping function ML5SOINDEX_FOR_SQUARED_ORDE'
-     $ //'RS'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' ML5SOINDEX_FOR_SQUARED_ORDERS'
       WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO)
       STOP
 
@@ -2170,8 +2174,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSO
-        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)+AMPSPLITORDERS(ORDERI
-     $NDEXB,I)
+        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)
+     $   +AMPSPLITORDERS(ORDERINDEXB,I)
       ENDDO
       ML5SQSOINDEX=ML5SOINDEX_FOR_SQUARED_ORDERS(SQORDERS)
       END
@@ -2209,8 +2213,8 @@
         RETURN
       ENDIF
 
-      WRITE(*,*) 'ERROR:: Stopping function ML5GET_SQUARED_ORDERS_FOR_S'
-     $ //'OINDEX'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' ML5GET_SQUARED_ORDERS_FOR_SOINDEX'
       WRITE(*,*) 'Could not find squared orders index ',SOINDEX
       STOP
 
@@ -2249,8 +2253,8 @@
         RETURN
       ENDIF
 
-      WRITE(*,*) 'ERROR:: Stopping function ML5GET_ORDERS_FOR_AMPSOINDE'
-     $ //'X'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' ML5GET_ORDERS_FOR_AMPSOINDEX'
       WRITE(*,*) 'Could not find amplitude split orders index ',SOINDEX
       STOP
 
@@ -2315,12 +2319,13 @@
 
       LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
      $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPU
-     $TATION_CHOICE
-      COMMON/COLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_
-     $COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_
-     $CHOICE
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE,
+     $  COLLIER_IR_POLE_COMPUTATION_CHOICE
+      COMMON/COLLIERPOLESFORCEDCHOICE
+     $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
+     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE
+     $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE
 
       COLLIERCOMPUTEUVPOLES                        = ONOFF
 C     This is just so that if we read the param again, we don't
@@ -2342,12 +2347,13 @@
 
       LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
      $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPU
-     $TATION_CHOICE
-      COMMON/COLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_
-     $COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_
-     $CHOICE
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE,
+     $  COLLIER_IR_POLE_COMPUTATION_CHOICE
+      COMMON/COLLIERPOLESFORCEDCHOICE
+     $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
+     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE
+     $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE
 
       COLLIERCOMPUTEIRPOLES         = ONOFF
 C     This is just so that if we read the param again, we don't

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%mp_compute_loop_coefs.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%mp_compute_loop_coefs.f	2017-07-07 10:47:20 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%mp_compute_loop_coefs.f	2020-06-21 12:17:46 +0000
@@ -268,13 +268,14 @@
 
       DPW_COPIED = .FALSE.
       DO H=1,NCOMB
-        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1).AND.(CHECKPHASE.OR.(
-     $.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H).GT.-HELOFFSET.AND.GOODHEL(H)
-     $   .NE.0)))) THEN
+        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1)
+     $   .AND.(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H)
+     $   .GT.-HELOFFSET.AND.GOODHEL(H).NE.0)))) THEN
 
 C         Handle the possible requirement of specific polarizations
-          IF ((.NOT.CHECKPHASE).AND.HELDOUBLECHECKED.AND.POLARIZATIONS(
-     $0,0).EQ.0.AND.(.NOT.IS_HEL_SELECTED(H))) THEN
+          IF ((.NOT.CHECKPHASE)
+     $     .AND.HELDOUBLECHECKED.AND.POLARIZATIONS(0,0)
+     $     .EQ.0.AND.(.NOT.IS_HEL_SELECTED(H))) THEN
             CYCLE
           ENDIF
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%polynomial.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%polynomial.f	2016-04-06 18:45:52 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%polynomial.f	2020-06-21 12:17:46 +0000
@@ -113,8 +113,8 @@
         CFTOT=CMPLX(CF_N(COLOR_ID,I)/(ONE*ABS(CF_D(COLOR_ID,I))),ZERO
      $   ,KIND=8)
         IF(CF_D(COLOR_ID,I).LT.0) CFTOT=CFTOT*IMAG1
-        CONST(ML5SOINDEX_FOR_BORN_AMP(I))=CONST(ML5SOINDEX_FOR_BORN_AMP
-     $(I))+CFTOT*CONJG(AMP(I))
+        CONST(ML5SOINDEX_FOR_BORN_AMP(I))
+     $   =CONST(ML5SOINDEX_FOR_BORN_AMP(I))+CFTOT*CONJG(AMP(I))
       ENDDO
 
       DO I=1,NAMPSO
@@ -232,8 +232,8 @@
         CFTOT=CMPLX(CF_N(COLOR_ID,I)/(ONE*ABS(CF_D(COLOR_ID,I))),ZERO
      $   ,KIND=16)
         IF(CF_D(COLOR_ID,I).LT.0) CFTOT=CFTOT*IMAG1
-        CONST(ML5SOINDEX_FOR_BORN_AMP(I))=CONST(ML5SOINDEX_FOR_BORN_AMP
-     $(I))+CFTOT*CONJG(AMP(I))
+        CONST(ML5SOINDEX_FOR_BORN_AMP(I))
+     $   =CONST(ML5SOINDEX_FOR_BORN_AMP(I))+CFTOT*CONJG(AMP(I))
       ENDDO
 
       DO I=1,NAMPSO

=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/matrix1.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/matrix1.f	2019-09-26 11:09:58 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/matrix1.f	2020-06-21 12:17:46 +0000
@@ -142,11 +142,13 @@
 
         !   If the helicity grid status is 0, this means that it is not yet initialized.
         !   If HEL_PICKED==-1, this means that calls to other matrix<i> where in initialization mode as well for the helicity.
-      IF ((ISHEL(IMIRROR).EQ.0.AND.ISUM_HEL.EQ.0).OR.(DS_GET_DIM_STATUS
-     $('Helicity').EQ.0).OR.(HEL_PICKED.EQ.-1)) THEN
+      IF ((ISHEL(IMIRROR).EQ.0.AND.ISUM_HEL.EQ.0)
+     $ .OR.(DS_GET_DIM_STATUS('Helicity').EQ.0).OR.(HEL_PICKED.EQ.-1))
+     $  THEN
         DO I=1,NCOMB
-          IF (GOODHEL(I,IMIRROR) .OR. NTRY(IMIRROR).LE.MAXTRIES.OR.(ISU
-     $M_HEL.NE.0).OR.THIS_NTRY(IMIRROR).LE.2) THEN
+          IF (GOODHEL(I,IMIRROR) .OR. NTRY(IMIRROR)
+     $     .LE.MAXTRIES.OR.(ISUM_HEL.NE.0).OR.THIS_NTRY(IMIRROR).LE.2)
+     $      THEN
             T=MATRIX1(P ,NHEL(1,I),JC(1))
             DO JJ=1,NINCOMING
               IF(POL(JJ).NE.1D0.AND.NHEL(JJ,I).EQ.INT(SIGN(1D0,POL(JJ))
@@ -174,8 +176,8 @@
           HEL_JACOBIAN   = 1.0D0
             !         We don't want to re-update the helicity grid if it was already updated by another matrix<i>, so we make sure that the reference grid is empty.
           REF_HELICITY_GRID = DS_GET_DIMENSION(REF_GRID,'Helicity')
-          IF((DS_GET_DIM_STATUS('Helicity').EQ.1).AND.(REF_HELICITY_GRI
-     $D%%N_TOT_ENTRIES.EQ.0)) THEN
+          IF((DS_GET_DIM_STATUS('Helicity').EQ.1)
+     $     .AND.(REF_HELICITY_GRID%%N_TOT_ENTRIES.EQ.0)) THEN
               !           If we finished the initialization we can update the grid so as to start sampling over it.
               !           However the grid will now be filled by dsample with different kind of weights (including pdf, flux, etc...) so by setting the grid_mode of the reference grid to 'initialization' we make sure it will be overwritten (as opposed to 'combined') by the running grid at the next update.
             CALL DS_UPDATE_GRID('Helicity')
@@ -436,8 +438,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSO
-        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)+AMPSPLITORDERS(ORDERI
-     $NDEXB,I)
+        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)
+     $   +AMPSPLITORDERS(ORDERINDEXB,I)
       ENDDO
       SQSOINDEX1=SOINDEX_FOR_SQUARED_ORDERS1(SQORDERS)
       END
@@ -531,8 +533,8 @@
         RETURN
       ENDIF
 
-      WRITE(*,*) 'ERROR:: Stopping function GET_SQUARED_ORDERS_FOR_SOIN'
-     $ //'DEX1'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' GET_SQUARED_ORDERS_FOR_SOINDEX1'
       WRITE(*,*) 'Could not find squared orders index ',SOINDEX
       STOP
 

=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f	2019-09-26 11:09:58 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f	2020-06-21 12:17:46 +0000
@@ -387,8 +387,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSO
-        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)+AMPSPLITORDERS(ORDERI
-     $NDEXB,I)
+        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)
+     $   +AMPSPLITORDERS(ORDERINDEXB,I)
       ENDDO
       SQSOINDEX=SOINDEX_FOR_SQUARED_ORDERS(SQORDERS)
       END
@@ -482,8 +482,8 @@
         RETURN
       ENDIF
 
-      WRITE(*,*) 'ERROR:: Stopping function GET_SQUARED_ORDERS_FOR_SOIN'
-     $ //'DEX'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' GET_SQUARED_ORDERS_FOR_SOINDEX'
       WRITE(*,*) 'Could not find squared orders index ',SOINDEX
       STOP
 

=== modified file 'tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_eq_4.f'
--- tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_eq_4.f	2018-05-25 09:13:24 +0000
+++ tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_eq_4.f	2020-06-21 12:17:46 +0000
@@ -83,8 +83,8 @@
       INTEGER NSQUAREDSOP1
       PARAMETER (NSQUAREDSOP1=NSQUAREDSO+1)
 C     The total number of loop reduction libraries
-C     At present, there are only CutTools,PJFry++,IREGI,Golem95,Samurai
-C     , Ninja and COLLIER
+C     At present, there are only
+C      CutTools,PJFry++,IREGI,Golem95,Samurai, Ninja and COLLIER
       INTEGER NLOOPLIB
       PARAMETER (NLOOPLIB=7)
 C     Only CutTools or possibly Ninja (if installed with qp support)
@@ -411,14 +411,15 @@
 C      reading the parameters
       LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
      $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPU
-     $TATION_CHOICE
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE,
+     $  COLLIER_IR_POLE_COMPUTATION_CHOICE
       DATA  FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION
      $ ,FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION/.FALSE.,.FALSE./
-      COMMON/ML5_0_COLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV
-     $_POLE_COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
-     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_
-     $CHOICE
+      COMMON/ML5_0_COLLIERPOLESFORCEDCHOICE
+     $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
+     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE
+     $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE
 
 C     This variable controls the general initialization which is
 C      *common* between all MadLoop SubProcesses.
@@ -888,8 +889,8 @@
      $ +P(2,2))**2-(P(3,1)+P(3,2))**2))
 
       CTCALL_REQ_SO_DONE=.FALSE.
-      FILTER_SO = (.NOT.CHECKPHASE).AND.HELDOUBLECHECKED.AND.(SQSO_TARG
-     $ET.NE.-1)
+      FILTER_SO = (.NOT.CHECKPHASE)
+     $ .AND.HELDOUBLECHECKED.AND.(SQSO_TARGET.NE.-1)
 
       DO I=1,NLOOPGROUPS
         DO J=0,LOOPMAXCOEFS-1
@@ -928,13 +929,14 @@
       ENDIF
 
       DO H=1,NCOMB
-        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1).AND.(CHECKPHASE.OR.(
-     $.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H).GT.-HELOFFSET.AND.GOODHEL(H)
-     $   .NE.0)))) THEN
+        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1)
+     $   .AND.(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H)
+     $   .GT.-HELOFFSET.AND.GOODHEL(H).NE.0)))) THEN
 
 C         Handle the possible requirement of specific polarizations
-          IF ((.NOT.CHECKPHASE).AND.HELDOUBLECHECKED.AND.POLARIZATIONS(
-     $0,0).EQ.0.AND.(.NOT.ML5_0_IS_HEL_SELECTED(H))) THEN
+          IF ((.NOT.CHECKPHASE)
+     $     .AND.HELDOUBLECHECKED.AND.POLARIZATIONS(0,0)
+     $     .EQ.0.AND.(.NOT.ML5_0_IS_HEL_SELECTED(H))) THEN
             CYCLE
           ENDIF
 
@@ -975,8 +977,9 @@
             DO I=1,NCTAMPS
               CFTOT=DCMPLX(CF_N(I,J)/DBLE(ABS(CF_D(I,J))),0.0D0)
               IF(CF_D(I,J).LT.0) CFTOT=CFTOT*IMAG1
-              ITEMP = ML5_0_ML5SQSOINDEX(ML5_0_ML5SOINDEX_FOR_LOOP_AMP(
-     $I),ML5_0_ML5SOINDEX_FOR_BORN_AMP(J))
+              ITEMP =
+     $          ML5_0_ML5SQSOINDEX(ML5_0_ML5SOINDEX_FOR_LOOP_AMP(I)
+     $         ,ML5_0_ML5SOINDEX_FOR_BORN_AMP(J))
               IF (.NOT.FILTER_SO.OR.SQSO_TARGET.EQ.ITEMP) THEN
                 DO K=1,3
                   TEMP2 = DBLE(CFTOT*AMPL(K,I)*CTEMP)
@@ -1102,8 +1105,8 @@
         IF((USERHEL.EQ.-1).OR.(USERHEL.EQ.HELPICKED)) THEN
 C         Make sure that that no polarization constraint filters out
 C          this helicity
-          IF (POLARIZATIONS(0,0).EQ.-1.OR.ML5_0_IS_HEL_SELECTED(HELPICK
-     $ED)) THEN
+          IF (POLARIZATIONS(0,0).EQ.
+     $     -1.OR.ML5_0_IS_HEL_SELECTED(HELPICKED)) THEN
 C           TO KEEP TRACK OF THE FINAL ANSWER TO BE RETURNED DURING
 C            CHECK PHASE
             DO I=0,NSQUAREDSO
@@ -1194,8 +1197,8 @@
 C             Of course if it is one, then we do not need to do
 C              anything (because with HELINITSTARTOVER=.FALSE. we only
 C              support exactly identical Hels.)
-              IF(GOODHEL(HELPICKED).GT.-HELOFFSET.AND.GOODHEL(HELPICKED)
-     $         .NE.1) THEN
+              IF(GOODHEL(HELPICKED).GT.
+     $         -HELOFFSET.AND.GOODHEL(HELPICKED).NE.1) THEN
                 NEWHELREF=-1
                 DO H=1,NCOMB
                   IF (GOODHEL(H).EQ.(-HELOFFSET-HELPICKED)) THEN
@@ -1389,16 +1392,18 @@
 
         CTMODE=BASIC_CT_MODE
 
-        IF(.NOT.EVAL_DONE(3).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
-     $.1).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.1)) ) THEN
+        IF(.NOT.EVAL_DONE(3).AND.
+     $    ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE.1)
+     $   .OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.1)) ) THEN
           EVAL_DONE(3)=.TRUE.
           CALL ML5_0_ROTATE_PS(PS,P,1)
           IF (DOING_QP_EVALS) CALL ML5_0_MP_ROTATE_PS(MP_PS,MP_P,1)
           GOTO 200
         ENDIF
 
-        IF(.NOT.EVAL_DONE(4).AND. ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE
-     $.2).OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.2)) ) THEN
+        IF(.NOT.EVAL_DONE(4).AND.
+     $    ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE.2)
+     $   .OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.2)) ) THEN
           EVAL_DONE(4)=.TRUE.
           CALL ML5_0_ROTATE_PS(PS,P,2)
           IF (DOING_QP_EVALS) CALL ML5_0_MP_ROTATE_PS(MP_PS,MP_P,2)
@@ -1957,8 +1962,8 @@
 C       When using COLLIER with the internal stability test, the first
 C        evaluation is typically more reliable so we do not want to
 C        use the average but rather the first evaluation.
-        IF (MLREDUCTIONLIB(I_LIB).EQ.7.AND.COLLIERUSEINTERNALSTABILITYT
-     $EST) THEN
+        IF (MLREDUCTIONLIB(I_LIB)
+     $   .EQ.7.AND.COLLIERUSEINTERNALSTABILITYTEST) THEN
           DO I=1,3
             ESTIMATE(I,K) = FULLLIST(I,K,1)
           ENDDO
@@ -2105,8 +2110,8 @@
  1009   CONTINUE
       ENDDO
 
-      WRITE(*,*) 'ERROR:: Stopping function ML5_0_ML5SOINDEX_FOR_SQUARE'
-     $ //'D_ORDERS'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS'
       WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO)
       STOP
 
@@ -2233,8 +2238,8 @@
 C     BEGIN CODE
 C     
       DO I=1,NSO
-        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)+AMPSPLITORDERS(ORDERI
-     $NDEXB,I)
+        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)
+     $   +AMPSPLITORDERS(ORDERINDEXB,I)
       ENDDO
       ML5_0_ML5SQSOINDEX=ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS(SQORDERS)
       END
@@ -2273,8 +2278,8 @@
         RETURN
       ENDIF
 
-      WRITE(*,*) 'ERROR:: Stopping function ML5_0_ML5GET_SQUARED_ORDERS'
-     $ //'_FOR_SOINDEX'
+      WRITE(*,*) 'ERROR:: Stopping function'
+     $ //' ML5_0_ML5GET_SQUARED_ORDERS_FOR_SOINDEX'
       WRITE(*,*) 'Could not find squared orders index ',SOINDEX
       STOP