Merge lp:~maddevelopers/mg5amcnlo/3.0.1_4FS_merged3.0.2_latest into lp:~maddevelopers/mg5amcnlo/3.0.3

[marco zaro]

Enables the computation of EW corrections in the 4FS

the possibility of computing EW corrections in the 4FS is added, including in the complex-mass scheme.

[lp:~maddevelopers/mg5amcnlo/3.0.1_4FS_merged3.0.2_latest]

953. By marco zaro on 2020-06-26

merged with 3.0.3

[Rikkert Frederix (frederix)]

Hi Marco,

It looks like you have turned on the compute_6to5flav_cnt() stuff (which is now relevant for the 5to4flav). This should be fixed. (Also, these terms should not be added to the Born -> rather, a new 'type' should be created (current list is in add_wgt() subroutine of fks_singular.f). This will probably require a couple of changes throughout the code, but will be the cleanest and most easy to maintain in the future.

Cheers,
Rikkert

review: Needs Fixing

[lp:~maddevelopers/mg5amcnlo/3.0.1_4FS_merged3.0.2_latest]

954. By marco zaro on 2020-06-26

restored 6to5 cnt function

[Olivier Mattelaer (olivier-mattelaer)]

If youu addd at new 'type'
Do not forget that it impacts reweighting/systematics (off-line version).
So please ping me such that I modify those accordingly

Olivier

[lp:~maddevelopers/mg5amcnlo/3.0.1_4FS_merged3.0.2_latest]

955. By marco zaro on 2020-08-25

merged with 3.0.4

956. By olivier-mattelaer on 2020-12-02

remove loop_qcd_qed_sm_Gmu_4FS and loop_qcd_qed_sm_4FS from model directory -> put online

957. By marco zaro on 2021-01-08

merged with 3.0.4 rev 973

Unmerged revisions

957. By marco zaro on 2021-01-08

merged with 3.0.4 rev 973

956. By olivier-mattelaer on 2020-12-02

remove loop_qcd_qed_sm_Gmu_4FS and loop_qcd_qed_sm_4FS from model directory -> put online

955. By marco zaro on 2020-08-25

merged with 3.0.4

954. By marco zaro on 2020-06-26

restored 6to5 cnt function

953. By marco zaro on 2020-06-26

merged with 3.0.3

952. By marco zaro on 2020-06-26

merged with 2.7.3

951. By marco zaro on 2020-05-21

merged with 3.0.2 rev 984

950. By marco zaro on 2020-05-20

update notes

949. By marco zaro on 2019-03-01

fix in the function that changes the number-of-flavour scheme.
      Was raising an error if a given coupling order could not be found.
            Now it has been fixed

948. By Hua-Sheng Shao on 2018-10-05

A crucial fix for B0F, where the tiny in sqrt should be smaller than that in log. The fix allows a correct Riemann sheet for B0F(MW**2,MB**2,MT**2).

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2020-06-21 21:31:33 +0000
+++ MadSpin/decay.py	2021-01-08 10:23:50 +0000
@@ -1,7 +1,12 @@
 #!/usr/bin/env python
 
 from __future__ import division
+from __future__ import absolute_import
+from __future__ import print_function
 from madgraph.interface import reweight_interface
+from six.moves import map
+from six.moves import range
+from six.moves import zip
 
 ################################################################################
 #
@@ -94,7 +99,7 @@
                 
         p=[]
         string=""
-        for id in xrange(len(self.particle)):
+        for id in range(len(self.particle)):
             particle = self.particle[map_event[id] + 1]
             if particle["istup"] < 2:
                 mom = particle["momentum"]
@@ -115,28 +120,26 @@
                 pattern = re.compile(r'''<\s*wgt id=[\'\"](?P<id>[^\'\"]+)[\'\"]\s*>\s*(?P<val>[\ded+-.]*)\s*</wgt>''')
                 data = pattern.findall(self.rwgt)
                 if len(data)==0:
-                    print self.rwgt
+                    print(self.rwgt)
                 try:
                     text = ''.join('   <wgt id=\'%s\'> %+15.7e </wgt>\n' % (pid, float(value) * factor)
                                      for (pid,value) in data) 
-                except ValueError, error:
-                    raise Exception, 'Event File has unvalid weight. %s' % error
+                except ValueError as error:
+                    raise Exception('Event File has unvalid weight. %s' % error)
                 self.rwgt = self.rwgt[:start] + '<rwgt>\n'+ text + self.rwgt[stop:]          
 
     def string_event_compact(self):
         """ return a string with the momenta of the event written 
                 in an easy readable way
         """
-        line=""
+        line=[]
         for part in range(1,len(self.particle)+1):
-            line+=str(self.particle[part]["pid"])+" "
-            line+=str(self.particle[part]["momentum"].px)+" "
-            line+=str(self.particle[part]["momentum"].py)+" "
-            line+=str(self.particle[part]["momentum"].pz)+" "
-            line+=str(self.particle[part]["momentum"].E)+"    " 
-            line+=str(self.particle[part]["momentum"].m)+"    " 
-            line+="\n"
-        return line
+            pid = self.particle[part]["pid"]
+            m = self.particle[part]["momentum"]
+            line.append("%i %s %.7g %.7g %.7g %.7g" % 
+                        (pid, m.px, m.py, m.pz, m.E, m.m))
+        line.append('')
+        return "\n".join(line)
     
     def get_tag(self):
         
@@ -170,7 +173,7 @@
         (self.nexternal,self.ievent,self.wgt,self.scale,self.aqed,self.aqcd)
         line+=line1+"\n"
         scales= []
-        for item in range(1,len(self.event2mg.keys())+1):
+        for item in range(1,len(list(self.event2mg.keys()))+1):
             part=self.event2mg[item]
             if part>0:
                 particle_line=self.get_particle_line(self.particle[part])
@@ -556,7 +559,7 @@
                     # I don't expect this to be inefficient, since there is a BW cut                                    
 
             if tree["nbody"] > 2:
-                raise Exception, 'Phase Space generator not yet ready for 3 body decay'
+                raise Exception('Phase Space generator not yet ready for 3 body decay')
 
             if ran==1:
                 decay_mom=generate_2body_decay(index2mom[res]["momentum"],mA, all_mass[0],all_mass[1])
@@ -874,7 +877,7 @@
             line+=str(d1)+" + "
             line+=str(d2)+" ,    type="
             line+=branch["type"]
-            print line
+            print(line)
 
 class AllMatrixElement(dict):
     """Object containing all the production topologies required for event to decay.
@@ -921,7 +924,7 @@
                 logger.warning("No Branching ratio applied for %s. Please check if this is expected" % init[0])
                 br *= self.get_br(decay)
             else:
-                raise MadGraph5Error,"No valid decay for %s. No 2 body decay for that particle. (three body are not supported by MadSpin)" % init[0]
+                raise MadGraph5Error("No valid decay for %s. No 2 body decay for that particle. (three body are not supported by MadSpin)" % init[0])
 
                 
 
@@ -1127,7 +1130,7 @@
 
                     # Reorganize s-channel vertices to get a list of all
                     # subprocesses for each vertex
-                    schannels = zip(*[s for s,t in stchannels])
+                    schannels = list(zip(*[s for s,t in stchannels]))
             else:
                     schannels = []
 
@@ -1168,7 +1171,7 @@
             if decay['decay_tag']==decay_tag: return decay
 
         msg = 'Unable to retrieve decay from decay_tag\n%s\n%s' %(production_tag, decay_tag)
-        raise Exception, msg
+        raise Exception(msg)
  
     def get_random_decay(self, production_tag,first=[]):
         """select randomly a decay channel"""
@@ -1384,32 +1387,32 @@
                     initial=item
                     if item not in [ particle['name'] for particle in base_model['particles'] ] \
                         and item not in [ particle['antiname'] for particle in base_model['particles'] ]:
-                        raise Exception, "No particle "+item+ " in the model "+model
+                        raise Exception("No particle "+item+ " in the model "+model)
                     continue
                 elif item=='>': continue       # case 2: we have the > symbole
                 elif item not in ponctuation : # case 3: we have a particle originating from a branching
                     final.append(item)
                     if next_symbol=='' or next_symbol in ponctuation:
                         #end of a splitting, verify that it exists
-                        if initial not in self.br.keys():
+                        if initial not in list(self.br.keys()):
                             logger.debug('Branching fractions of particle '+initial+' are unknown')
                             return 0
                         if len(final)>2:
-                            raise Exception, 'splittings different from A > B +C are currently not implemented '
+                            raise Exception('splittings different from A > B +C are currently not implemented ')
 
-                        if final[0] in multiparticles.keys():
+                        if final[0] in list(multiparticles.keys()):
                             set_B=[pid2label[pid] for pid in multiparticles[final[0]]]
                         else:
                             if final[0] not in [ particle['name'] for particle in base_model['particles'] ] \
                                and final[0] not in [ particle['antiname'] for particle in base_model['particles'] ]:
-                                raise Exception, "No particle "+item+ " in the model "
+                                raise Exception("No particle "+item+ " in the model ")
                             set_B=[final[0]]
-                        if final[1] in multiparticles.keys():
+                        if final[1] in list(multiparticles.keys()):
                             set_C=[pid2label[pid] for pid in multiparticles[final[1]]]
                         else:
                             if final[1] not in [ particle['name'] for particle in base_model['particles'] ] \
                                and final[1] not in [ particle['antiname'] for particle in base_model['particles'] ]:
-                                raise Exception, "No particle "+item+ " in the model "+model
+                                raise Exception("No particle "+item+ " in the model "+model)
                             set_C=[final[1]]
 
                         splittings={}
@@ -1538,7 +1541,7 @@
         ff.close()
         
         lhapdf = False
-        if os.environ.has_key('lhapdf'):
+        if 'lhapdf' in os.environ:
             lhapdf = os.environ['lhapdf']
             del os.environ['lhapdf']
         
@@ -1566,7 +1569,7 @@
         
         label2pid = self.label2pid
         pid2label = self.label2pid
-        for res in self.br.keys():
+        for res in list(self.br.keys()):
             particle=self.model.get_particle(label2pid[res])
             if particle['self_antipart']: 
                 continue
@@ -1588,7 +1591,7 @@
                 self.br[anti_res][chan]['daughters']=[]
                 self.br[anti_res][chan]['daughters'].append(d1bar)
                 self.br[anti_res][chan]['daughters'].append(d2bar)
-                if decay.has_key('width'):
+                if 'width' in decay:
                     self.br[anti_res][chan]['width']=decay['width']                  
 
     def launch_width_evaluation(self,resonances, mgcmd):
@@ -1847,7 +1850,7 @@
             if line =="": break
             list_line=line.split()
             if len(list_line)>2:
-                if list_line[0]=="DECAY" and int(list_line[1]) in pid2widths.keys():
+                if list_line[0]=="DECAY" and int(list_line[1]) in list(pid2widths.keys()):
                     list_line[2]=str(pid2widths[int(list_line[1])]) 
                     line=""
                     for item in list_line:
@@ -2049,7 +2052,7 @@
         else:
             try:
                 self.all_ME, self.all_decay,self.width_estimator = save_load_object.load_from_file(pjoin(self.path_me,"production_me", "all_ME.pkl"))
-            except Exception,error:
+            except Exception as error:
                 logger.debug(str(error))
                 self.generate_all_matrix_element()
                 self.save_to_file(pickle_info,
@@ -2310,7 +2313,7 @@
     Please relaunch MS with more events/PS point by event in the
     computation of the maximum_weight.
                     """ % (report['over_weight'], event_nb, 100 * report['over_weight']/event_nb )  
-                    raise MadSpinError, error
+                    raise MadSpinError(error)
                         
                     error = True
                 elif report['%s_f' % (decay['decay_tag'],)] > max(0.01*report[decay['decay_tag']],3):
@@ -2321,7 +2324,7 @@
                             report['%s' % (decay['decay_tag'],)],\
                             100 * report['%s_f' % (decay['decay_tag'],)] / report[ decay['decay_tag']] ,\
                             decay['decay_tag'])  
-                    raise MadSpinError, error
+                    raise MadSpinError(error)
                     
              
             decayed_event.change_wgt(factor= self.branching_ratio) 
@@ -2491,7 +2494,6 @@
             return decay_mapping
         
         BW_cut = self.options['BW_cut']       
-        
         #class the decay by class (nbody/pid)
         nbody_to_decay = collections.defaultdict(list)
         for decay in self.all_decay.values():
@@ -2525,7 +2527,6 @@
                     '\n'.join(str(tree[i]['momentum']) for i in range(1, len(tree))
                                                                   if i in tree))
                 
-                
                 values = {}                
                 for i in range(len(decays)):
                     if any([valid[(i,j)] for j in range(len(decays)) if i !=j]):
@@ -2559,7 +2560,6 @@
                         comment+=  "%4e " % valid[(i,j)]
                     comment+= "|"+ os.path.basename(decays[i]['path'])                     
                     logger.debug(comment)
-            
             # store the result in the relation object. (using tag as key)
             for i in range(len(decays)):
                 tag_i = decays[i]['tag'][2:]
@@ -2654,7 +2654,7 @@
             pass
         try:
             shutil.rmtree(pjoin(path_me,'production_me'))
-        except Exception, error:
+        except Exception as error:
             pass
         path_me = self.path_me        
         
@@ -2683,7 +2683,7 @@
                 if name == 'all':
                     continue
                 #self.banner.get('proc_card').get('multiparticles'):
-                mgcmd.do_define("%s = %s" % (name, ' '.join(`i` for i in pdgs)))
+                mgcmd.do_define("%s = %s" % (name, ' '.join(repr(i) for i in pdgs)))
             
         
         mgcmd.exec_cmd("set group_subprocesses False")
@@ -2701,7 +2701,7 @@
         
         mgcmd.exec_cmd(commandline, precmd=True)
         commandline = 'output standalone_msP %s %s' % \
-        (pjoin(path_me,'production_me'), ' '.join(self.list_branches.keys()))        
+        (pjoin(path_me,'production_me'), ' '.join(list(self.list_branches.keys())))        
         mgcmd.exec_cmd(commandline, precmd=True)        
         logger.info('Done %.4g' % (time.time()-start))
 
@@ -2722,19 +2722,19 @@
                 if self.all_ME.has_particles_ambiguity:
                     final_states.add(self.pid2label[-1*self.pid2label[label]])
                 final_states.add(label)
-        for key in self.list_branches.keys():
+        for key in list(self.list_branches.keys()):
             if key not in final_states and key not in self.mgcmd._multiparticles:
                 if (len(self.list_branches)>1):
                     del self.list_branches[key]
                 elif not self.options["onlyhelicity"]:
-                    raise Exception, " No decay define for process."
+                    raise Exception(" No decay define for process.")
                     logger.info('keeping dummy decay for passthrough mode')
 
         # 4. compute the full matrix element -----------------------------------
         if not self.options["onlyhelicity"]:
             logger.info('generating the full matrix element squared (with decay)')
             start = time.time()
-            to_decay = self.mscmd.list_branches.keys()
+            to_decay = list(self.mscmd.list_branches.keys())
             decay_text = []
             for decays in self.mscmd.list_branches.values():
                 for decay in  decays:
@@ -2757,13 +2757,13 @@
             # remove decay with 0 branching ratio.
             mgcmd.remove_pointless_decay(self.banner.param_card)
             commandline = 'output standalone_msF %s %s' % (pjoin(path_me,'full_me'),
-                                                          ' '.join(self.list_branches.keys()))
+                                                          ' '.join(list(self.list_branches.keys())))
             mgcmd.exec_cmd(commandline, precmd=True)
             logger.info('Done %.4g' % (time.time()-start))
         elif self.options["onlyhelicity"]:
             logger.info("Helicity Matrix-Element")      
             commandline = 'output standalone_msF %s %s' % \
-            (pjoin(path_me,'full_me'), ' '.join(self.list_branches.keys()))        
+            (pjoin(path_me,'full_me'), ' '.join(list(self.list_branches.keys())))        
             mgcmd.exec_cmd(commandline, precmd=True)        
             logger.info('Done %.4g' % (time.time()-start))                    
 
@@ -2875,7 +2875,7 @@
                 try:
                     proc_nb = int(proc_nb)
                 except ValueError:
-                    raise MadSpinError, 'MadSpin didn\'t allow order restriction after the @ comment: \"%s\" not valid' % proc_nb
+                    raise MadSpinError('MadSpin didn\'t allow order restriction after the @ comment: \"%s\" not valid' % proc_nb)
                 proc_nb = '@%i' % proc_nb 
             else:
                 baseproc = new_proc
@@ -2896,7 +2896,7 @@
             else:
                 part = baseproc.split(',')
                 if any('(' in p for p in part):
-                    raise Exception, 'too much decay at MG level. this can not be done for the moment)'            
+                    raise Exception('too much decay at MG level. this can not be done for the moment)')            
                 else:
                     decay_part = []
                     for p in part[1:]:
@@ -2912,7 +2912,7 @@
         # Compute the width branching ratio. Doing this at this point allows
         #to remove potential pointless decay in the diagram generation.
         resonances = decay_misc.get_all_resonances(self.banner, 
-                         self.mgcmd, self.mscmd.list_branches.keys())
+                         self.mgcmd, list(self.mscmd.list_branches.keys()))
 
         logger.debug('List of resonances:%s' % resonances)
         path_me = os.path.realpath(self.path_me) 
@@ -3030,9 +3030,9 @@
 
                 if __debug__:
                     if(os.path.getsize(pjoin(path_me, mode,'SubProcesses', 'parameters.inc'))<10):
-                        print pjoin(path_me, mode,'SubProcesses', 'parameters.inc')
-                        raise Exception, "Parameters of the model were not written correctly ! %s " %\
-                            os.path.getsize(pjoin(path_me, mode,'SubProcesses', 'parameters.inc'))
+                        print(pjoin(path_me, mode,'SubProcesses', 'parameters.inc'))
+                        raise Exception("Parameters of the model were not written correctly ! %s " %\
+                            os.path.getsize(pjoin(path_me, mode,'SubProcesses', 'parameters.inc')))
 
 
     def extract_resonances_mass_width(self, resonances):
@@ -3050,7 +3050,7 @@
                 part=abs(self.pid2label[particle_label])
                 #mass = self.banner.get('param_card','mass', abs(part))
                 width = self.banner.get('param_card','decay', abs(part))
-            except ValueError, error:
+            except ValueError as error:
                 continue
             else:
                 if (width.value > 0.001):  
@@ -3153,7 +3153,6 @@
                     continue # No decay for this process
                 atleastonedecay = True
                 weight = self.get_max_weight_from_fortran(decay['path'], event_map,numberps,self.options['BW_cut'])
-                
                     #weight=mg5_me_full*BW_weight_prod*BW_weight_decay/mg5_me_prod
                 if tag in max_decay:
                     max_decay[tag] = max([max_decay[tag], weight])
@@ -3291,10 +3290,12 @@
                                path, std_in)
 
         return max_weight
-        
+    
+    nb_load = 0
     def loadfortran(self, mode, path, stdin_text, first=True):
         """ call the fortran executable """
 
+        self.nb_load +=1
         tmpdir = ''
         if ('full',path) in self.calculator:
             external = self.calculator[('full',path)]
@@ -3304,17 +3305,32 @@
             tmpdir = path
 
             executable_prod="./check"
+            my_env = os.environ.copy()
+            my_env["GFORTRAN_UNBUFFERED_ALL"] = "y"
             external = Popen(executable_prod, stdout=PIPE, stdin=PIPE, 
-                                                      stderr=STDOUT, cwd=tmpdir)
+                                          stderr=STDOUT, cwd=tmpdir, env=my_env)
             self.calculator[('full',path,)] = external 
             self.calculator_nbcall[('full',path)] = 1 
-
+            
         try:
-            external.stdin.write(stdin_text)
-        except IOError:
+            external.stdin.write(stdin_text.encode())
+            external.stdin.flush()
+        except IOError as error:
             if not first:
                 raise
             try:
+                external.stdin.close()
+            except Exception as  error:
+                misc.sprint(error)
+            try:
+                external.stdout.close()
+            except Exception as error:
+                misc.sprint(error)
+            try:
+                external.stderr.close()
+            except Exception as error:
+                misc.sprint(error)
+            try:
                 external.terminate()
             except:
                 pass
@@ -3346,7 +3362,7 @@
 
         if len(self.calculator) > self.options['max_running_process']:
             logger.debug('more than %s calculators. Perform cleaning' % self.options['max_running_process'])
-            nb_calls = self.calculator_nbcall.values()
+            nb_calls = list(self.calculator_nbcall.values())
             nb_calls.sort()
             cut = max([nb_calls[len(nb_calls)//2], 0.001 * nb_calls[-1]])
             for key, external in list(self.calculator.items()):
@@ -3375,7 +3391,7 @@
         """routine to return the matrix element"""
 
         if mode != "decay":
-            raise Exception, "This function is only secure in mode decay."
+            raise Exception("This function is only secure in mode decay.")
 
         tmpdir = ''
         if (mode, production) in self.calculator:
@@ -3390,30 +3406,40 @@
                 tmpdir = pjoin(self.path_me,'%s_me' % mode, 'SubProcesses',
                            production)
             executable_prod="./check"
-            external = Popen(executable_prod, stdout=PIPE, stdin=PIPE, 
-                                                      stderr=STDOUT, cwd=tmpdir)
+            my_env = os.environ.copy()
+            my_env["GFORTRAN_UNBUFFERED_ALL"] = "y"
+            external = Popen(executable_prod, stdout=PIPE, stdin=PIPE,
+                                                      stderr=STDOUT, cwd=tmpdir,
+                                                      env=my_env,
+                                                      bufsize=0)
+            assert (mode, production) not in self.calculator
             self.calculator[(mode, production)] = external 
             self.calculator_nbcall[(mode, production)] = 1       
 
-        external.stdin.write(stdin_text)
+
+        external.stdin.write(stdin_text.encode())
         if mode == 'prod':
-            info = int(external.stdout.readline())
+            info = int(external.stdout.readline().decode())
             nb_output = abs(info)+1
         else:
             info = 1
             nb_output = 1
-         
-        prod_values = ' '.join([external.stdout.readline() for i in range(nb_output)])
+        std = []
+        for i in range(nb_output):
+            external.stdout.flush()
+            line = external.stdout.readline().decode()
+            std.append(line)
+        prod_values = ' '.join(std)
+        #prod_values = ' '.join([external.stdout.readline().decode() for i in range(nb_output)])
         if info < 0:
-            print 'ZERO DETECTED'
-            print prod_values
-            print stdin_text
+            print('ZERO DETECTED')
+            print(prod_values)
+            print(stdin_text)
             os.system('lsof -p %s' % external.pid)
             return ' '.join(prod_values.split()[-1*(nb_output-1):])
-        
         if len(self.calculator) > self.options['max_running_process']:
             logger.debug('more than 100 calculator. Perform cleaning')
-            nb_calls = self.calculator_nbcall.values()
+            nb_calls = list(self.calculator_nbcall.values())
             nb_calls.sort()
             cut = max([nb_calls[len(nb_calls)//2], 0.001 * nb_calls[-1]])
             for key, external in list(self.calculator.items()):
@@ -3426,7 +3452,6 @@
                     del self.calculator_nbcall[key]
                 else:
                     self.calculator_nbcall[key] = self.calculator_nbcall[key] //10
-        
         if mode == 'prod':
             return prod_values
         else:
@@ -3457,7 +3482,7 @@
                 #print part
                 decay_struct[part]['mg_tree']=[]
 
-                nb_res=len(decay_struct[part]["tree"].keys())
+                nb_res=len(list(decay_struct[part]["tree"].keys()))
                 for res in range(-1,-nb_res-1,-1):
                     label=abs(decay_struct[part]["tree"][res]['label'])
                     mass=self.pid2massvar[label]
@@ -3533,7 +3558,7 @@
                 else:
                     # now we need to write the decay products in the event
                     # follow the decay chain order, so that we can easily keep track of the mother index                                           
-                    for res in range(-1,-len(decay_struct[part]["tree"].keys())-1,-1):
+                    for res in range(-1,-len(list(decay_struct[part]["tree"].keys()))-1,-1):
                         index_d1=decay_struct[part]['mg_tree'][-res-1][1]
                         index_d2=decay_struct[part]['mg_tree'][-res-1][2]
                         
@@ -3596,7 +3621,7 @@
                     # follow the decay chain order, so that we can easily keep track of the mother index
                        
                     map_to_part_number={}
-                    for res in range(-1,-len(decay_struct[part]["tree"].keys())-1,-1):
+                    for res in range(-1,-len(list(decay_struct[part]["tree"].keys()))-1,-1):
                         index_res_for_mom=decay_struct[part]['mg_tree'][-res-1][0]
                         if (res==-1):
                             part_number+=1
@@ -3736,8 +3761,8 @@
                                 d2colup2=colup2
                                 d2colup1=maxcol                        
                         else:
-                            raise Exception, 'color combination not treated by MadSpin (yet). (%s,%s,%s)' \
-                                % (colord1,colord2,colormother)
+                            raise Exception('color combination not treated by MadSpin (yet). (%s,%s,%s)' \
+                                % (colord1,colord2,colormother))
                         part_number+=1
                         index_d1_for_mom=decay_struct[part]['mg_tree'][-res-1][1]
                         mom=momenta_in_decay[index_d1_for_mom].copy()
@@ -3819,14 +3844,14 @@
                     mothup1=curr_event.resonance[part]["mothup1"]         
                     mothup2=curr_event.resonance[part]["mothup2"] 
                     if mothup1==index:
-                        if mothup2!=index: print "Warning: mothup1!=mothup2"
+                        if mothup2!=index: print("Warning: mothup1!=mothup2")
                         curr_event.resonance[part]["mothup1"]=part_number
                         curr_event.resonance[part]["mothup2"]=part_number
                 for part in curr_event.particle.keys():
                     mothup1=curr_event.particle[part]["mothup1"]         
                     mothup2=curr_event.particle[part]["mothup2"] 
                     if mothup1==index:
-                        if mothup2!=index: print "Warning: mothup1!=mothup2"
+                        if mothup2!=index: print("Warning: mothup1!=mothup2")
                         curr_event.particle[part]["mothup1"]=part_number
                         curr_event.particle[part]["mothup2"]=part_number
 
@@ -3844,7 +3869,7 @@
         max_br = max([m['total_br'] for m in self.all_ME.values()])
         if max_br >= 1:
             if max_br > 1.0001:
-                raise MadSpinError, 'BR is larger than one.'
+                raise MadSpinError('BR is larger than one.')
             max_br = 1
         for production in self.all_ME.values():
             if production['total_br'] < max_br:
@@ -3966,30 +3991,45 @@
             for (mode, path) in self.calculator:
                 if mode=='decay':
                     external = self.calculator[(mode, path)]
+                    try:
+                        external.stdin.close()
+                    except Exception as  error:
+                        misc.sprint(error)
+                        continue
+                    try:
+                        external.stdout.close()
+                    except Exception as error:
+                        misc.sprint(error)
+                        continue
                     external.terminate()
                     del external
                 elif mode=='full':
-                    stdin_text="5 0 0 0 0\n"  # before closing, write down the seed 
+                    stdin_text="5 0 0 0 0\n".encode()  # before closing, write down the seed 
                     external = self.calculator[('full',path)]
                     try:
                         external.stdin.write(stdin_text)
-                    except Exception:
+                        external.stdin.flush() 
+                    except Exception as error:
+                        misc.sprint(error)
+                        raise
                         continue
-                    ranmar_state=external.stdout.readline()
+                    ranmar_state=external.stdout.readline().decode()
                     ranmar_file=pjoin(path,'ranmar_state.dat')
                     ranmar=open(ranmar_file, 'w')
                     ranmar.write(ranmar_state)
                     ranmar.close()
                     try:
                         external.stdin.close()
-                    except Exception:
-                        continue
+                    except Exception as  error:
+                        misc.sprint(error)
                     try:
                         external.stdout.close()
-                    except Exception:
-                        continue
+                    except Exception as error:
+                        misc.sprint(error)                   
                     external.terminate()
                     del external
+                else:
+                    misc.sprint('not closed', mode, type(mode))
         else:
             try:
                 external = self.calculator[('full', path_to_decay)]
@@ -4002,7 +4042,7 @@
                 external.stdout.close()
                 external.terminate()       
                 del external
-
+            
         self.calculator = {}
 
 
@@ -4058,7 +4098,7 @@
                 if name == 'all':
                     continue
                 #self.banner.get('proc_card').get('multiparticles'):
-                mgcmd.do_define("%s = %s" % (name, ' '.join(`i` for i in pdgs)))
+                mgcmd.do_define("%s = %s" % (name, ' '.join(repr(i) for i in pdgs)))
             
         
         mgcmd.exec_cmd("set group_subprocesses False")
@@ -4087,7 +4127,7 @@
         # 4. compute the full matrix element -----------------------------------
         logger.info('generating the full matrix element squared (with decay)')
 #        start = time.time()
-        to_decay = self.mscmd.list_branches.keys()
+        to_decay = list(self.mscmd.list_branches.keys())
         decay_text = []
         for decays in self.mscmd.list_branches.values():
             for decay in  decays:
@@ -4172,6 +4212,8 @@
 
     def compile(self):
         logger.info('Compiling code')
+        #my_env = os.environ.copy()
+        #os.environ["GFORTRAN_UNBUFFERED_ALL"] = "y"
         misc.compile(cwd=pjoin(self.path_me,'madspin_me', 'Source'),
                      nb_core=self.mgcmd.options['nb_core'])        
         misc.compile(['all'],cwd=pjoin(self.path_me,'madspin_me', 'SubProcesses'),
 
=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2020-06-16 19:24:01 +0000
+++ MadSpin/interface_madspin.py	2021-01-08 10:23:50 +0000
@@ -14,6 +14,8 @@
 ################################################################################
 """ Command interface for MadSpin """
 from __future__ import division
+from __future__ import absolute_import
+from __future__ import print_function
 import collections
 import logging
 import math
@@ -24,6 +26,8 @@
 import sys
 import time
 import glob
+import six
+from six.moves import range
 
 
 pjoin = os.path.join
@@ -150,6 +154,7 @@
         self.mg5cmd = master_interface.MasterCmd()
         self.seed = None
         self.err_branching_ratio = 0
+        self.me_run_name = "" # Events diretory name where to stotre the events (used by madevent) not use internally
         
         
         if event_path:
@@ -254,7 +259,7 @@
         process = self.banner.get_detail('proc_card', 'generate')
         if not process:
             msg = 'Invalid proc_card information in the file (no generate line):\n %s' % self.banner['mg5proccard']
-            raise Exception, msg
+            raise Exception(msg)
         process, option = mg_interface.MadGraphCmd.split_process_line(process)
         self.proc_option = option
         
@@ -414,7 +419,7 @@
                 logger.warning("polarization option used with spinmode=onshell. This combination is not validated and is by construction using sub-optimal method which can likely lead to bias in some situation. Use at your own risk.")
 
         decay_process, init_part = self.decay.reorder_branch(decaybranch)
-        if not self.list_branches.has_key(init_part):
+        if init_part not in self.list_branches:
             self.list_branches[init_part] = []
         self.list_branches[init_part].append(decay_process)
         del decay_process, init_part    
@@ -491,7 +496,7 @@
 
         # Format
         if len(args) == 1:
-            opts = self.options.keys() + ['seed', "spinmode"]
+            opts = list(self.options.keys()) + ['seed', "spinmode"]
             return self.list_completion(text, opts) 
         elif len(args) == 2:
             if args[1] == 'BW_effect':
@@ -508,20 +513,20 @@
     def help_set(self):
         """help the set command"""
         
-        print 'syntax: set OPTION VALUE'
-        print ''
-        print '-- assign to a given option a given value'
-        print '   - set max_weight VALUE: pre-define the maximum_weight for the reweighting'
-        print '   - set BW_effect True|False: [default:True] reshuffle the momenta to describe'
-        print '       corrrectly the Breit-Wigner of the decayed particle'
-        print '   - set seed VALUE: fix the value of the seed to a given value.'
-        print '       by default use the current time to set the seed. random number are'
-        print '       generated by the python module random using the Mersenne Twister generator.'
-        print '       It has a period of 2**19937-1.'
-        print '   - set max_running_process VALUE: allow to limit the number of open file used by the code'
-        print '       The number of running is raising like 2*VALUE'
-        print '   - set spinmode=none: mode with simple file merging. No spin correlation attempt.'
-        print '       This mode allows 3 (and more) body decay.'
+        print('syntax: set OPTION VALUE')
+        print('')
+        print('-- assign to a given option a given value')
+        print('   - set max_weight VALUE: pre-define the maximum_weight for the reweighting')
+        print('   - set BW_effect True|False: [default:True] reshuffle the momenta to describe')
+        print('       corrrectly the Breit-Wigner of the decayed particle')
+        print('   - set seed VALUE: fix the value of the seed to a given value.')
+        print('       by default use the current time to set the seed. random number are')
+        print('       generated by the python module random using the Mersenne Twister generator.')
+        print('       It has a period of 2**19937-1.')
+        print('   - set max_running_process VALUE: allow to limit the number of open file used by the code')
+        print('       The number of running is raising like 2*VALUE')
+        print('   - set spinmode=none: mode with simple file merging. No spin correlation attempt.')
+        print('       This mode allows 3 (and more) body decay.')
     
     def do_define(self, line):
         """ """
@@ -547,14 +552,14 @@
         """ """
         try:
             return self.mg5cmd.complete_define(*args)
-        except Exception,error:
+        except Exception as error:
             misc.sprint(error)
             
     def complete_decay(self, *args):
         """ """
         try:
             return self.mg5cmd.complete_generate(*args)
-        except Exception,error:
+        except Exception as error:
             misc.sprint(error)
             
     def check_launch(self, args):
@@ -568,29 +573,53 @@
         # Validity check. Need lhe version 3 if matching is on
         if self.banner.get("run_card", "lhe_version") < 3 and \
             self.banner.get("run_card", "ickkw") > 0:
-            raise Exception, "MadSpin requires LHEF version 3 when running with matching/merging"
+            raise Exception("MadSpin requires LHEF version 3 when running with matching/merging")
 
     def help_launch(self):
         """help for the launch command"""
         
-        print '''Running Madspin on the loaded events, following the decays enter
+        print('''Running Madspin on the loaded events, following the decays enter
         An example of a full run is the following:
         import ../mssm_events.lhe.gz
         define sq = ur ur~
         decay go > sq j
         launch
-        '''
+        ''')
+        
+        self.parser_launch.print_help()
+
+    def parser_launch(self):
+        usage = """launch [-n RUN_NAME]   
+        """
+        parser = misc.OptionParser(usage=usage)
+        parser.add_option("-n", "--name",
+                  default="",
+                  help="When NOT run in standalone instruct MG5aMC where to store the events file")
+        return parser
+    
+    def parse_launch(self, line):
+        
+        args = self.split_arg(line)
+        return self.parser_launch().parse_args(args)
+        
 
     @misc.mute_logger()
     def do_launch(self, line):
         """end of the configuration launched the code"""
         
+        (options, args) = self.parse_launch(line)
+        
+        if options.name:
+            self.me_run_name = options.name # Only use by MG5aMC
+        else:
+            self.me_run_name = ''
+        
         if self.options["spinmode"] in ["none"]:
             return self.run_bridge(line)
         elif self.options["spinmode"] == "onshell":
             return self.run_onshell(line)
         elif self.options["spinmode"] == "bridge":
-            raise Exception, "Bridge mode not available."
+            raise Exception("Bridge mode not available.")
         
         if self.options['ms_dir'] and os.path.exists(pjoin(self.options['ms_dir'], 'madspin.pkl')):
             return self.run_from_pickle()
@@ -622,7 +651,7 @@
 
         if self.options['seed'] > 30081*30081: # can't use too big random number
             msg = 'Random seed too large ' + str(self.options['seed']) + ' > 30081*30081'
-            raise Exception, msg
+            raise Exception(msg)
 
         #self.options['seed'] = self.seed
         text = '%s\n' % '\n'.join([ line for line in self.history if line])
@@ -716,13 +745,13 @@
                         
             orig_block = generate_all.banner.param_card[name]
             if block != orig_block:                
-                raise Exception, """The directory %s is specific to a mass spectrum. 
+                raise Exception("""The directory %s is specific to a mass spectrum. 
                 Your event file is not compatible with this one. (Different param_card: %s different)
                 orig block:
                 %s
                 new block:
                 %s""" \
-                % (self.options['ms_dir'], name, orig_block, block)
+                % (self.options['ms_dir'], name, orig_block, block))
 
         #replace init information
         generate_all.banner['init'] = self.banner['init']
@@ -808,7 +837,7 @@
             elif '.hepmc' in filename:
                 self.options['input_format']  = 'hepmc'
             else:
-                raise Exception, "fail to recognized input format automatically"
+                raise Exception("fail to recognized input format automatically")
                 
         if self.options['input_format'] in ['lhe', 'lhe_no_banner']:
             orig_lhe = lhe_parser.EventFile(filename)
@@ -849,7 +878,7 @@
 
         if self.options['seed'] > 30081*30081: # can't use too big random number
             msg = 'Random seed too large ' + str(self.options['seed']) + ' > 30081*30081'
-            raise Exception, msg
+            raise Exception(msg)
 
         #self.options['seed'] = self.options['seed']
         
@@ -989,7 +1018,7 @@
 
         for event in orig_lhe:
             if counter and counter % 100 == 0 and float(str(counter)[1:]) ==0:
-                print "decaying event number %s [%s s]" % (counter, time.time()-start)
+                print("decaying event number %s [%s s]" % (counter, time.time()-start))
             counter +=1
             
             # use random order for particles to avoid systematics when more than 
@@ -1058,7 +1087,7 @@
                     # read the event file up to completion
                     while 1:
                         try:
-                            decay = decay_file.next()
+                            decay = next(decay_file)
                         except StopIteration:
                             # check how far we are
                             ratio = counter / nb_event 
@@ -1069,6 +1098,7 @@
                             evt_decayfile[particle.pdg].update(new_file)
                             decay_file = evt_decayfile[particle.pdg][decay_file_nb]
                             continue
+
                         if helicity == decay[0].helicity or helicity==9 or \
                                             self.options["spinmode"] == "none":
                             break # use that event
@@ -1198,7 +1228,7 @@
                         try:
                             os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
                             os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
-                        except Exception,error:
+                        except Exception as error:
                             logger.debug(error)
                             pass 
                         self.me_int[decay_dir] = me5_cmd
@@ -1245,7 +1275,7 @@
                     try:
                         os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
                         os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
-                    except Exception,error:
+                    except Exception as error:
                         logger.debug(error)
                         pass                 
                     self.me_int[decay_dir] = me5_cmd
@@ -1282,6 +1312,16 @@
                 me5_cmd.exec_cmd("exit")
                 out[i] = lhe_parser.EventFile(pjoin(decay_dir, "Events", 'run_01', 'unweighted_events.lhe.gz'))            
             else:
+                if not self.seed:
+                    if hasattr(self, 'mother'):
+                        try:
+                            self.seed = 100 + self.mother.run_card['iseed']
+                        except:
+                            self.seed = random.randint(0, int(30081*30081))
+                self.seed += 1
+                if self.seed > 30081*30081:
+                    self.seed -= 30081*30081        
+                logger.info('Will use seed %s' % (self.seed))
                 misc.call(['run.sh', str(int(1.2*nb_event)), str(self.seed)], cwd=decay_dir)     
                 out[i] = lhe_parser.EventFile(pjoin(decay_dir, 'events.lhe.gz'))            
             if cumul:
@@ -1360,7 +1400,7 @@
 
         if self.seed > 30081*30081: # can't use too big random number
             msg = 'Random seed too large ' + str(self.seed) + ' > 30081*30081'
-            raise Exception, msg
+            raise Exception(msg)
 
         self.options['seed'] = self.seed
         
@@ -1538,7 +1578,7 @@
             # file if needed.
             while 1:
                 try:
-                    decay = decay_file.next()
+                    decay = next(decay_file)
                     break
                 except StopIteration:
                     eff = self.efficiency
@@ -1580,7 +1620,7 @@
                 logger.info( "Event %s/%s :  %2fs" % (i, nevents, time.time()-start))
             maxwgt = 0
             orig_lhe.seek(0)
-            base_event = orig_lhe.next()
+            base_event = next(orig_lhe)
             if self.options['fixed_order']:
                 base_event = base_event[0]
             for j in range(self.options['max_weight_ps_point']):
@@ -1674,8 +1714,12 @@
                 sys.path.insert(0, pjoin(self.path_me, 'madspin_me', 'SubProcesses'))
             
             mymod = __import__("%s.matrix2py" % (pdir))
+            if six.PY3:
+                from importlib import reload
+            else:
+                from imp import reload
             reload(mymod)
-            mymod = getattr(mymod, 'matrix2py')
+            mymod = getattr(mymod, 'matrix2py')  
             with misc.chdir(pjoin(self.path_me, 'madspin_me', 'SubProcesses', pdir)):
                 with misc.stdchannel_redirected(sys.stdout, os.devnull):
                     if not os.path.exists(pjoin(self.path_me, 'Cards','param_card.dat')) and \
@@ -1718,7 +1762,7 @@
                 try:
                     proc_nb = int(proc_nb)
                 except ValueError:
-                    raise MadSpinError, 'MadSpin didn\'t allow order restriction after the @ comment: \"%s\" not valid' % proc_nb
+                    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'])
 
=== modified file 'Template/LO/Cards/pythia8_card_default.dat'
--- Template/LO/Cards/pythia8_card_default.dat	2017-02-06 18:35:42 +0000
+++ Template/LO/Cards/pythia8_card_default.dat	2021-01-08 10:23:50 +0000
@@ -12,6 +12,8 @@
 ! -------------------------------------------------------------------
 ! Specify the HEPMC output of the Pythia8 shower. You can set it to:
 !   auto      : MG5aMC will automatically place it the run_<i> directory
+!   autoremove: MG5aMC will automatically remove the file at the end of the run.
+!                (usefull when running with Delphes)
 !   /dev/null : to turn off the HEPMC output.
 !   <path>    : to select where the HEPMC file must written. It will 
 !               therefore not be placed in the run_<i> directory. The
 
=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2020-05-27 07:52:34 +0000
+++ Template/LO/Cards/run_card.dat	2021-01-08 10:23:50 +0000
@@ -28,7 +28,7 @@
 #*********************************************************************
 # Collider type and energy                                           *
 # lpp: 0=No PDF, 1=proton, -1=antiproton, 2=photon from proton,      *
-#                                         3=photon from electron     *
+#                3=photon from electron, 4=photon from muon          *
 #*********************************************************************
      %(lpp1)s        = lpp1    ! beam 1 type 
      %(lpp2)s        = lpp2    ! beam 2 type
 
=== modified file 'Template/LO/Source/PDF/PhotonFlux.f'
--- Template/LO/Source/PDF/PhotonFlux.f	2020-03-27 13:59:33 +0000
+++ Template/LO/Source/PDF/PhotonFlux.f	2021-01-08 10:23:50 +0000
@@ -6,25 +6,28 @@
 c/* ********************************************************* */
 c   provided by Tomasz Pierzchala - UCL
 
-      real*8 function epa_electron(x,q2max)
-      integer i
+      real*8 function epa_lepton(x,q2max, mode)
+      implicit none
+      integer i, mode, imode
+c     mode is +3/-3  for electron and +4/-4 for muon      
       real*8 x,phi_f
-      real*8 xin
+      real*8 xin(3:4)
       real*8 alpha
       real*8 f, q2min,q2max
       real*8 PI
       data PI/3.14159265358979323846/
 
-      data xin/0.511d-3/ !electron mass in GeV
+      data xin/0.511d-3, 0.105658d0/ !electron mass in GeV
 
       alpha = .0072992701
-
+      imode = abs(mode)
+      
 C     // x = omega/E = (E-E')/E
       if (x.lt.1) then
-         q2min= xin*xin*x*x/(1-x)
+         q2min= xin(imode)*xin(imode)*x*x/(1-x)
          if(q2min.lt.q2max) then 
-            f = alpha/2d0/PI*
-     &           (2d0*xin*xin*x*(-1/q2min+1/q2max)+
+             f = alpha/2d0/PI*
+     &           (2d0*xin(imode)*xin(imode)*x*(-1/q2min+1/q2max)+
      &           (2-2d0*x+x*x)/x*dlog(q2max/q2min))
             
          else
@@ -35,7 +38,7 @@
       endif
 c      write (*,*) x,dsqrt(q2min),dsqrt(q2max),f
       if (f .lt. 0) f = 0
-      epa_electron= f
+      epa_lepton = f
 
       end
 
 
=== modified file 'Template/LO/Source/PDF/pdf.f'
--- Template/LO/Source/PDF/pdf.f	2018-02-01 23:17:18 +0000
+++ Template/LO/Source/PDF/pdf.f	2021-01-08 10:23:50 +0000
@@ -29,7 +29,7 @@
       double precision u_val,d_val,u_sea,d_sea,s_sea,c_sea,b_sea,gluon
       double precision Ctq3df,Ctq4Fn,Ctq5Pdf,Ctq6Pdf,Ctq5L
       double precision q2max
-      double precision epa_electron,epa_proton
+      double precision epa_lepton,epa_proton
       include 'pdf.inc'
 
       integer mode,Iprtn,Irt
@@ -297,8 +297,8 @@
 c  a "diffractive" photon
 c      
       q2max=xmu*xmu
-      if(ih .eq. 3) then  !from the electron
-          fx(7)=epa_electron(x,q2max)
+      if(abs(ih) .eq. 3.or.abs(ih) .eq. 4) then  !from the electron
+          fx(7)=epa_lepton(x,q2max, ih)
       elseif(ih .eq. 2) then  !from a proton without breaking
           fx(7)=epa_proton(x,q2max)
       endif      
 
=== modified file 'Template/LO/Source/PDF/pdg2pdf.f'
--- Template/LO/Source/PDF/pdg2pdf.f	2019-12-04 21:34:25 +0000
+++ Template/LO/Source/PDF/pdg2pdf.f	2021-01-08 10:23:50 +0000
@@ -24,7 +24,7 @@
       integer mode,Irt,i,j
       double precision xlast(2),xmulast(2),pdflast(-7:7,2),q2max
       character*7 pdlabellast(2)
-      double precision epa_electron,epa_proton
+      double precision epa_lepton,epa_proton
       integer ipart,ireuse,iporg,ihlast(2)
       save xlast,xmulast,pdflast,pdlabellast,ihlast
       data xlast/2*-99d9/
@@ -131,8 +131,8 @@
 
       if(iabs(ipart).eq.7.and.ih.gt.1) then
          q2max=xmu*xmu
-         if(ih.eq.3) then       !from the electron
-            pdg2pdf=epa_electron(x,q2max)
+         if(abs(ih).eq.3.or.abs(ih).eq.4) then       !from the electron or muonn
+            pdg2pdf=epa_lepton(x,q2max, ih)
          elseif(ih .eq. 2) then !from a proton without breaking
             pdg2pdf=epa_proton(x,q2max,beamid)
          endif 
 
=== removed file 'Template/LO/Source/PDF/pdg2pdf_lhapdf.f'
=== modified file 'Template/LO/Source/PDF/pdg2pdf_lhapdf6.f'
--- Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2020-03-06 20:47:59 +0000
+++ Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2021-01-08 10:23:50 +0000
@@ -24,7 +24,7 @@
       integer i,j,ihlast(20),ipart,iporg,ireuse,imemlast(20),iset,imem
      &     ,i_replace,ii,ipartlast(20)
       double precision xlast(20),xmulast(20),pdflast(-7:7,20)
-      double precision epa_proton, epa_electron
+      double precision epa_proton, epa_lepton
       save ihlast,xlast,xmulast,pdflast,imemlast,ipartlast
       data ihlast/20*-99/
       data ipartlast/20*-99/
@@ -139,8 +139,8 @@
             pdg2pdf = get_ion_pdf(pdflast(-7, i_replace), ipart, nb_proton(beamid), nb_neutron(beamid))
          endif
          pdg2pdf=pdg2pdf/x
-      else if(ih.eq.3) then       !from the electron
-            pdg2pdf=epa_electron(x,xmu*xmu)
+      else if(abs(ih).eq.3.or.abs(ih).eq.4) then       !from the electron
+            pdg2pdf=epa_lepton(x,xmu*xmu, ih)
       else if(ih.eq.2) then ! photon from a proton without breaking
           pdg2pdf = epa_proton(x,xmu*xmu)
       else
 
=== modified file 'Template/LO/Source/kin_functions.f'
--- Template/LO/Source/kin_functions.f	2016-02-10 14:56:17 +0000
+++ Template/LO/Source/kin_functions.f	2021-01-08 10:23:50 +0000
@@ -590,11 +590,13 @@
 C     4-Vector Dot product
 C****************************************************************************
       implicit none
-      double precision p1(0:3),p2(0:3)
+      double precision p1(0:3),p2(0:3), dot2
       dot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
 
-      if(dabs(dot).lt.1d-6)then ! solve numerical problem 
-         dot=0d0
+      if(dabs(dot).lt.1d-6)then ! solve numerical problem
+         dot2 = p1(0)*p2(0)+p1(1)*p2(1)+p1(2)*p2(2)+p1(3)*p2(3)
+         dot2 = max(1e-99,dot2)
+         if (dot/dot2.lt.1e-6) dot=0d0
       endif
 
       end
 
=== modified file 'Template/LO/Source/setrun.f'
--- Template/LO/Source/setrun.f	2016-03-04 16:16:23 +0000
+++ Template/LO/Source/setrun.f	2021-01-08 10:23:50 +0000
@@ -83,6 +83,22 @@
       q2fact(1) = sf1**2      ! fact scale**2 for pdf1
       q2fact(2) = sf2**2      ! fact scale**2 for pdf2     
 
+      if(pb1.ne.0d0)then
+         if (abs(lpp(1)).eq.1.or.abs(lpp(1)).eq.2)then
+            write(*,*) 'proton/anti-proton beam polarization are not allowed'
+            stop 1
+         endif
+         pol(1)=sign(1+abs(pb1)/100d0,pb1)
+      endif
+      if(pb2.ne.0d0)then
+         if (abs(lpp(2)).eq.1.or.abs(lpp(2)).eq.2)then
+            write(*,*) 'proton/anti-proton beam polarization are not allowed'
+            stop 1
+         endif
+         pol(2)=sign(1+abs(pb2)/100d0,pb2)
+      endif
+
+      
       if(pb1.ne.0d0.and.lpp(1).eq.0) pol(1)=sign(1+abs(pb1)/100d0,pb1)
       if(pb2.ne.0d0.and.lpp(2).eq.0) pol(2)=sign(1+abs(pb2)/100d0,pb2)
 
 
=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2020-06-15 07:16:51 +0000
+++ Template/LO/SubProcesses/genps.f	2021-01-08 10:23:50 +0000
@@ -572,6 +572,8 @@
             if (abs(lpp(2)) .eq. 1 .or. abs(lpp(2)) .eq. 2) m2 = 0.938d0
             if (abs(lpp(1)) .eq. 3) m1 = 0.000511d0
             if (abs(lpp(2)) .eq. 3) m2 = 0.000511d0
+            if (abs(lpp(1)) .eq. 4) m1 = 0.105658d0
+            if (abs(lpp(2)) .eq. 4) m2 = 0.105658d0
             if (mass_ion(1).ge.0d0) m1 = mass_ion(1)
             if (mass_ion(2).ge.0d0) m2 = mass_ion(2)
             if(ebeam(1).lt.m1.and.lpp(1).ne.9) ebeam(1)=m1
@@ -1046,7 +1048,6 @@
             ptotm(i) = ptot(i)
          endif
       enddo
-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-05-24 20:16:56 +0000
+++ Template/LO/SubProcesses/reweight.f	2021-01-08 10:23:50 +0000
@@ -1738,9 +1738,9 @@
          if (btest(mlevel,3))
      $        write(*,*)' set fact scales for PS to ',
      $        sqrt(q2fact(1)),sqrt(q2fact(2))
-      else if (abs(lpp(1)).eq.2.or.abs(lpp(1)).eq.3) then
+      else if (abs(lpp(1)).ge.2.and.abs(lpp(1)).le.4) then
          q2fact(1)=q2bck(1)
-      else if (abs(lpp(2)).eq.2.or.abs(lpp(2)).eq.3) then
+      else if (abs(lpp(2)).ge.2.or.abs(lpp(2)).le.4) then
          q2fact(2)=q2bck(2)
       endif
 
 
=== modified file 'Template/LO/SubProcesses/unwgt.f'
--- Template/LO/SubProcesses/unwgt.f	2020-05-20 11:52:45 +0000
+++ Template/LO/SubProcesses/unwgt.f	2021-01-08 10:23:50 +0000
@@ -145,10 +145,19 @@
 C     compute the boost for the requested transformation
       implicit none
       double precision pin(0:3), pout(0:3)
-      double precision denom
+      double precision denom, control
 
       denom = pin(0)*pout(0) - pin(3)*pout(3)
       if (denom.ne.0d0) then
+         control = pin(0)*pout(0) + pin(3)*pout(3)
+         if (denom/control.gt.1d-12)then
+            get_betaz = (pin(3) * pout(0) - pout(3) * pin(0)) / denom
+c        pass in massless case
+         else if (abs(pin(0)+pin(3)).gt.abs(pin(0))) then
+            get_betaz = (pin(0)**2 - pout(0)**2)/(pin(0)**2 + pout(0)**2)
+         else
+            get_betaz = (pout(0)**2 - pin(0)**2)/(pin(0)**2 + pout(0)**2)
+         endif
          get_betaz = (pin(3) * pout(0) - pout(3) * pin(0)) / denom
       else if (pin(0).eq.pin(3)) then
          get_betaz = (pin(0)**2 - pout(0)**2)/(pin(0)**2 + pout(0)**2)
 
=== modified file 'Template/LO/bin/generate_events'
--- Template/LO/bin/generate_events	2018-03-19 21:56:39 +0000
+++ Template/LO/bin/generate_events	2021-01-08 10:23:50 +0000
@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 ################################################################################
 #
 # Copyright (c) 2011 The MadGraph5_aMC@NLO Development team and Contributors
@@ -27,10 +27,17 @@
 root_path = os.path.split(os.path.dirname(os.path.realpath( __file__ )))[0]
 pjoin = os.path.join
 
-if not sys.version_info[0] == 2 or sys.version_info[1] < 6:
-    sys.exit('MadEvent works with python 2.6 or higher (but not python 3.X).\n\
+if sys.version_info[1] < 7:
+    sys.exit('MadEvent works with python 2.7/3.7 or higher.\n\
                Please upgrade your version of python.')
 
+try:
+    import six
+except ImportError:
+    message = 'madgraph requires the six module. The easiest way to install it is to run "pip install six --user"\n'
+    message += 'in case of problem with pip, you can download the file at https://pypi.org/project/six/ . It has a single python file that you just need to put inside a directory of your $PYTHONPATH environment variable.'
+    sys.exit(message)
+    
 # Check if optimize mode is (and should be) activated
 if __debug__ and (not os.path.exists(pjoin(root_path,'../..', 'bin','create_release.py'))):
     subprocess.call([sys.executable] + ['-O'] + sys.argv)
@@ -50,7 +57,7 @@
     try:
         import pyreadline as readline
     except:
-        print "For tab completion and history, install module readline."
+        print("For tab completion and history, install module readline.")
 else:
     import rlcompleter
 
@@ -81,7 +88,7 @@
    pass
 
 if __debug__:
-        print 'Running MG5 in debug mode'
+        print('Running MG5 in debug mode')
 
 
 
@@ -99,25 +106,25 @@
     try:
         mode = int(argument[1])
     except:
-        mode = int(raw_input('Enter 2 for multi-core, 1 for parallel, 0 for serial run\n'))
+        mode = int(six.moves.input('Enter 2 for multi-core, 1 for parallel, 0 for serial run\n'))
     if mode == 0:
         try:
             name = argument[2]
         except:
-            name = raw_input('Enter run name\n')
+            name = six.moves.input('Enter run name\n')
     else:
         try:
             opt = argument[2]
         except:
             if mode == 1: 
-                opt = raw_input('Enter name for jobs on pbs queue\n')
+                opt = six.moves.input('Enter name for jobs on pbs queue\n')
             else:
-                opt = int(raw_input('Enter number of cores\n'))
+                opt = int(six.moves.input('Enter number of cores\n'))
         
         try:
             name = argument[3]
         except:
-            name = raw_input('enter run name\n')
+            name = six.moves.input('enter run name\n')
 
 #    launch = ME.MadEventCmd(me_dir=root_path)
         
@@ -158,7 +165,7 @@
             launch = ME.MadEventCmdShell(me_dir=root_path, force_run=True)
             launch.run_cmd('generate_events %s' % ' '.join(argument[1:]))
             launch.run_cmd('quit')
-    except ME.MadEventAlreadyRunning, message:
+    except ME.MadEventAlreadyRunning as message:
         logging.getLogger('madgraph').log(40, str(message))
         sys.exit(1)
     except KeyboardInterrupt:
@@ -166,7 +173,7 @@
             launch.run_cmd('quit')
         except:
             pass
-    except Exception, error:
+    except Exception as error:
         logging.error(str(error))
         sys.exit()
             
 
=== modified file 'Template/LO/bin/internal/Gridpack/gridrun'
--- Template/LO/bin/internal/Gridpack/gridrun	2017-06-08 18:43:40 +0000
+++ Template/LO/bin/internal/Gridpack/gridrun	2021-01-08 10:23:50 +0000
@@ -19,10 +19,17 @@
 and call immediately the command line interface scripts"""
 
 import sys
-if not sys.version_info[0] == 2 or sys.version_info[1] < 6:
-    sys.exit('MadGraph/MadEvent 5 works only with python 2.6 or later (but not python 3.X).\n\
+if sys.version_info[1] < 7:
+    sys.exit('MadGraph/MadEvent 5 works only with python 2.7/3.7 or later.\n\
                Please upgrate your version of python.')
 
+try:
+    import six
+except ImportError:
+    message = 'madgraph requires the six module. The easiest way to install it is to run "pip install six --user"\n'
+    message += 'in case of problem with pip, you can download the file at https://pypi.org/project/six/ . It has a single python file that you just need to put inside a directory of your $PYTHONPATH environment variable.'
+    sys.exit(message)
+    
 import os
 import optparse
 
 
=== modified file 'Template/LO/bin/internal/addmasses_optional.py'
--- Template/LO/bin/internal/addmasses_optional.py	2015-10-15 20:22:39 +0000
+++ Template/LO/bin/internal/addmasses_optional.py	2021-01-08 10:23:50 +0000
@@ -5,8 +5,11 @@
 # Add masses to charged leptons, fix kinematics
 # Insert W/Z when missing
 
+from __future__ import absolute_import
+from __future__ import print_function
 import math,sys,re
 from xml.dom import minidom
+from six.moves import range
 
 #tolerance for energy momentum conservation
 toler = 1e-4
@@ -70,7 +73,7 @@
         self.pz = self.pz/pi*po
     def printMe(self):
         li = [self.px,self.py, self.pz, self.E, self.m]
-        print "| %18.10E %18.10E %18.10E %18.10E %18.10E |" % tuple(li)    
+        print("| %18.10E %18.10E %18.10E %18.10E %18.10E |" % tuple(li))    
 
 #useful class to describe a particle
 class Particle:
@@ -87,7 +90,7 @@
         self.polar = l[12]
     def printMe(self):
         li = [self.no, self.id, self.status,self.mo1, self.mo2, self.co1, self.co2, self.mom.px,self.mom.py, self.mom.pz, self.mom.E, self.mom.m, self.life, self.polar]
-        print "%2i | %9i | %4i | %4i %4i | %4i %4i | %18.10E %18.10E %18.10E %18.10E %18.10E | %1.0f. %2.0f" % tuple(li)
+        print("%2i | %9i | %4i | %4i %4i | %4i %4i | %18.10E %18.10E %18.10E %18.10E %18.10E | %1.0f. %2.0f" % tuple(li))
     def writeMe(self):
         li = [self.id, self.status,self.mo1, self.mo2, self.co1, self.co2, self.mom.px,self.mom.py, self.mom.pz, self.mom.E, self.mom.m, self.life, self.polar]
         return "%9i %4i %4i %4i %4i %4i %18.10E %18.10E %18.10E %18.10E %18.10E  %1.0f. %2.0f\n" % tuple(li)        
@@ -108,7 +111,7 @@
         try:
             line=f.readline()
         except IOError:
-            print "Problem reading from file ",sys.argv[1]
+            print("Problem reading from file ",sys.argv[1])
             sys.exit(0)
         if line.find("<event>")==-1:
             g.write(line)
@@ -121,7 +124,7 @@
     try:
         xmldoc = minidom.parse(sys.argv[1])
     except IOError:
-        print " could not open file for xml parsing ",sys.argv[1]
+        print(" could not open file for xml parsing ",sys.argv[1])
         sys.exit(0)
         
         
@@ -178,7 +181,7 @@
             if len(noMotherList)==0:
                 pass
             elif len(noMotherList)%2 != 0:
-                print "single orphan; do not know how to process"
+                print("single orphan; do not know how to process")
             else:
                 ki=0
                 while ki<len(noMotherList)-1:
@@ -265,7 +268,7 @@
                 pSum = pSum + p.mom
     
         if abs(pSum.px)>toler or abs(pSum.py)>toler or abs(pSum.pz)>toler or abs(pSum.E)>toler:
-            print "Event does not pass tolerance ",toler
+            print("Event does not pass tolerance ",toler)
             pSum.printMe()
     
         if 1:
@@ -286,28 +289,28 @@
 
     #main part of analysis
     if len(sys.argv)!=3:
-        print "Usage: addmasses.py <infile> <outfile>    "
-        print " Last modified: Fri Nov 21 10:49:14 CST 2008 "
+        print("Usage: addmasses.py <infile> <outfile>    ")
+        print(" Last modified: Fri Nov 21 10:49:14 CST 2008 ")
         sys.exit(1)
     else:
-        print "Running addmasses.py to add masses and correct kinematics of fixed particles"
+        print("Running addmasses.py to add masses and correct kinematics of fixed particles")
     
     #first print out leading information
     try:
         f=open(sys.argv[1],'r')
     except IOError:
-        print "need a file for reading"
+        print("need a file for reading")
         sys.exit(1)
 
     try:
         g=open(sys.argv[2],'w')
     except IOError:
-        print "need a file for writing"
+        print("need a file for writing")
         sys.exit(1)
     
     try:
         add_masses(f,g)
-    except Exception, error:
-        print "addmasses failed with error, %s" % error
+    except Exception as error:
+        print("addmasses failed with error, %s" % error)
         sys.exit(1)
     
 
=== modified file 'Template/LO/bin/madevent'
--- Template/LO/bin/madevent	2018-03-30 07:54:37 +0000
+++ Template/LO/bin/madevent	2021-01-08 10:23:50 +0000
@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 ################################################################################
 #
@@ -19,10 +19,18 @@
 and call immediately the command line interface scripts"""
 
 import sys
-if not sys.version_info[0] == 2 or sys.version_info[1] < 6:
-    sys.exit('MadGraph/MadEvent 5 works only with python 2.6 or later (but not python 3.X).\n\
+if sys.version_info[1] < 7:
+    sys.exit('MadGraph/MadEvent 5 works only with python 2.7/3.7 or later.\n\
                Please upgrate your version of python.')
 
+try:
+    import six
+except ImportError:
+    message = 'madgraph requires the six module. The easiest way to install it is to run "pip install six --user"\n'
+    message += 'in case of problem with pip, you can download the file at https://pypi.org/project/six/ . It has a single python file that you just need to put inside a directory of your $PYTHONPATH environment variable.'
+    sys.exit(message)
+
+    
 import os
 import optparse
 
@@ -59,7 +67,7 @@
     try:
         (options, args) = parser.parse_args(sys.argv[1:len(sys.argv)-i])
         done = True
-    except MyOptParser.InvalidOption, error:
+    except MyOptParser.InvalidOption as error:
         pass 
     else:
         args += sys.argv[len(sys.argv)-i:]
@@ -67,8 +75,8 @@
     # raise correct error:
     try:
         (options, args) = parser.parse_args()
-    except MyOptParser.InvalidOption, error:
-        print error
+    except MyOptParser.InvalidOption as error:
+        print(error)
         sys.exit(2)  
 
 if len(args) == 0:
@@ -79,7 +87,7 @@
 # Check if optimize mode is (and should be) activated
 if __debug__ and not options.debug and \
     (not os.path.exists(os.path.join(root_path,'../..', 'bin','create_release.py')) or options.web):
-        print 'launch in debug mode'
+        print( 'launch in debug mode')
         subprocess.call([sys.executable] + ['-O'] + sys.argv)
         sys.exit()
 
@@ -92,7 +100,7 @@
     try:
         import pyreadline as readline
     except:
-        print "For tab completion and history, install module readline."
+        print( "For tab completion and history, install module readline.")
 else:
     import rlcompleter
 
@@ -123,7 +131,7 @@
    pass
 
 if __debug__:
-        print 'Running MG5 in debug mode'
+        print( 'Running MG5 in debug mode')
 
 
 # Set logging level according to the logging level given by options
@@ -175,10 +183,10 @@
                 cmd_line = cmd_interface.MadEventCmdShell(force_run=True)
             if not hasattr(cmd_line, 'do_%s' % args[0]):
                 if parser_error:
-                    print parser_error
-                    print 'and %s  can not be interpreted as a valid command.' % args[0]
+                    print( parser_error)
+                    print( 'and %s  can not be interpreted as a valid command.' % args[0])
                 else:
-                    print 'ERROR: %s  not a valid command. Please retry' % args[0]
+                    print( 'ERROR: %s  not a valid command. Please retry' % args[0])
             else:
                 cmd_line.use_rawinput = False    
                 cmd_line.run_cmd(' '.join(args))
@@ -192,9 +200,9 @@
                 cmd_line = cmd_interface.MadEventCmdShell(force_run=True)
                 cmd_line.cmdloop()
 except KeyboardInterrupt:
-    print 'writting history and directory and quit on KeyboardInterrupt' 
+    print( 'writting history and directory and quit on KeyboardInterrupt' )
     pass
-except cmd_interface.MadEventAlreadyRunning, error:
+except cmd_interface.MadEventAlreadyRunning  as  error:
     logging.error(str(error))
     sys.exit(1)
 
 
=== modified file 'Template/MadWeight/Python/Info.py'
--- Template/MadWeight/Python/Info.py	2015-10-01 16:00:08 +0000
+++ Template/MadWeight/Python/Info.py	2021-01-08 10:23:50 +0000
@@ -1,18 +1,21 @@
-#!/usr/bin/env python
-
-def giveInfo(class_):
-        if type(class_)!=str:
-                class_=class_.__class__.__name__
-        for info in dir(eval(class_)):
-                
-                print class_+'.'+info+' : ',eval(class_+'.'+info+'.__doc__')
-
-
-
-if __name__=='__main__':
-    class_=raw_input('enter the name of the class')
-    try:
-        import class_
-    except:
-        pass
-    giveInfo(class_)
+#!/usr/bin/env python
+
+from __future__ import absolute_import
+from __future__ import print_function
+from six.moves import input
+def giveInfo(class_):
+        if type(class_)!=str:
+                class_=class_.__class__.__name__
+        for info in dir(eval(class_)):
+                
+                print(class_+'.'+info+' : ',eval(class_+'.'+info+'.__doc__'))
+
+
+
+if __name__=='__main__':
+    class_=input('enter the name of the class')
+    try:
+        import class_
+    except:
+        pass
+    giveInfo(class_)
 
=== modified file 'Template/MadWeight/Python/clean.py'
--- Template/MadWeight/Python/clean.py	2015-10-01 16:00:08 +0000
+++ Template/MadWeight/Python/clean.py	2021-01-08 10:23:50 +0000
@@ -33,6 +33,8 @@
 ##
 ## BEGIN INCLUDE
 ##
+from __future__ import absolute_import
+from __future__ import print_function
 import os
 
 #1 #########################################################################
@@ -48,8 +50,8 @@
         if os.path.isdir(element) and (element[0]=='P' or element[:4]=='MW_P'):
             status,mess=clean.suppress_dir(os.path.join(element,run_name))
             if not status:
-                print 'supress ',element,' failed:' 
-                print mess
+                print('supress ',element,' failed:') 
+                print(mess)
 
     os.chdir(os.pardir)
     
@@ -161,14 +163,14 @@
 
         if suppres and os.path.isfile(file):
             if 'n' in opt:
-                print 'schedulle deleting file',file
+                print('schedulle deleting file',file)
                 return 1,mess
             try:
                 os.remove(file)
             except:
                 pass
             if 's' not in opt:
-                print 'delete file',file
+                print('delete file',file)
             return 1,''
         elif suppres:
             return 0, 'not a file'+file
@@ -221,9 +223,9 @@
             if 'n' not in opt:
                 os.rmdir(pos)
                 if 's' not in opt:
-                    print 'delete dir',pos
+                    print('delete dir',pos)
             else:
-                print 'schedulle removing directory',pos
+                print('schedulle removing directory',pos)
                 
             return 1,''
         else:
@@ -237,4 +239,4 @@
     go_to_main_dir()
     #Clean_event('fermi')
     #Clean_run('fermi')
-    print 'no cleaning by default'
+    print('no cleaning by default')
 
=== modified file 'Template/MadWeight/Python/expand_MadWeight.py'
--- Template/MadWeight/Python/expand_MadWeight.py	2011-10-25 14:42:02 +0000
+++ Template/MadWeight/Python/expand_MadWeight.py	2021-01-08 10:23:50 +0000
@@ -1,6 +1,8 @@
 #!/usr/bin/env python
 
 # Module
+from __future__ import absolute_import
+from __future__ import print_function
 import string
 import os
 import sys
@@ -11,10 +13,10 @@
 
 def expand_all(echap=[]):
     go_to_main_dir()
-    print "copying files"
+    print("copying files")
     copy_file(echap)
     del_file()
-    print "modifying files"
+    print("modifying files")
     opt={}
     opt['nowarning']="""['DESACTIVATE_CUT','DESACTIVATE_BW_CUT','get_user_params','main_make','obj_for_MW']"""
     mod_file.mod_file('./Source/MadWeight/mod_file/MW_pos',opt=opt)
 
=== modified file 'Template/MadWeight/Python/madweight.py'
--- Template/MadWeight/Python/madweight.py	2013-03-29 20:38:55 +0000
+++ Template/MadWeight/Python/madweight.py	2021-01-08 10:23:50 +0000
@@ -6,12 +6,15 @@
 
 
 #Extension
+from __future__ import absolute_import
+from __future__ import print_function
 import string
 import sys
 import os
 import re
 import time
 import stat
+from six.moves import input
 # patch if symbolic directory replace by real file
 sys.path.append('./Source/MadWeight/Python')
 sys.path.append('../Source/MadWeight/Python')
@@ -37,7 +40,7 @@
 def Launch_all_SubProcess(MWparam):
 
     name=MWparam.name
-    print 'name :',name
+    print('name :',name)
     P_proclist,MW_proclist=MWparam.P_listdir,MWparam.MW_listdir
     #create banner
     if MWparam.run_opt['launch']:
@@ -46,14 +49,14 @@
         banner.write()
 
     if MWparam.run_opt['compilation']:
-        print 'starting program compilation'
+        print('starting program compilation')
         compile_SubProcesses(MW_proclist)
             
     if MWparam.run_opt['event']:
         verif_event(MWparam)
 
     if MWparam.run_opt['refine']:
-        print "collecting data to find data with a precision less than",MWparam.run_opt['refine']
+        print("collecting data to find data with a precision less than",MWparam.run_opt['refine'])
         collect_schedular(MWparam)          
 
 
@@ -62,7 +65,7 @@
     cluster.driver()
 
     if MWparam.run_opt['collect']:
-        print "collecting data"
+        print("collecting data")
         collect_schedular(MWparam)          
 
 
@@ -90,7 +93,7 @@
         if  os.path.isfile("./comp_madweight") and exit_status==0 :
             os.chdir("..")
         else:
-            print "fortran compilation error"
+            print("fortran compilation error")
             sys.exit()
     os.chdir("..")   
     return
@@ -109,7 +112,7 @@
         if  os.path.isfile("./madevent") and exit_status==0:
             os.chdir("..")
         else:
-            print "fortran compilation error"
+            print("fortran compilation error")
             sys.exit()
     os.chdir("..")   
     return  
@@ -142,15 +145,15 @@
 
     for dir in MWparam.MW_listdir:
       if not os.path.exists("SubProcesses/"+dir+"/call_TF.f"): 
-        print "Currently no transfer function loaded ..."
+        print("Currently no transfer function loaded ...")
         listdir=os.listdir('./Source/MadWeight/transfer_function/data')
-        print 'Available transfer functions:\n   ',
-        print '\n    '.join([content[3:-4] for content in listdir if (content.startswith('TF') and content.endswith('dat'))])
-        name=raw_input('Choose your transfer Function\n')
+        print('Available transfer functions:\n   ', end=' ')
+        print('\n    '.join([content[3:-4] for content in listdir if (content.startswith('TF') and content.endswith('dat'))]))
+        name=input('Choose your transfer Function\n')
         P_dir,MW_dir=detect_SubProcess(P_mode=1)
         os.chdir('./Source/MadWeight/transfer_function')
         change_tf.create_TF_main(name,0,MW_dir)
-        print os.getcwd()
+        print(os.getcwd())
         os.chdir('../../..')
         break
 
@@ -173,7 +176,7 @@
             plot.Differential_Graph(MWparam,auto=1)        
 
     if MWparam.run_opt['clean']:
-        print 'cleaning in progress ....'
+        print('cleaning in progress ....')
         from clean import Clean_run
         if MWparam.run_opt['clean']==1:
             Clean_run(MWparam.name)
 
=== modified file 'Template/MadWeight/Python/put_banner.py'
--- Template/MadWeight/Python/put_banner.py	2011-10-25 14:42:02 +0000
+++ Template/MadWeight/Python/put_banner.py	2021-01-08 10:23:50 +0000
@@ -99,7 +99,11 @@
 
 
 
+from __future__ import absolute_import
+from __future__ import print_function
 import os,re,sys
+from six.moves import range
+from six.moves import input
 sys.path+=['../'*i+'./Source/MadWeight/Python' for i in range(1,6)]
 import MW_param
 
@@ -140,7 +144,7 @@
         """    |   put in variable full_banner_txt the content of the file      ##
         ##     |   defined in position self.header_file                         ##
         """
-        self.full_banner_txt+=file(self.header_file,'rU').read()
+        self.full_banner_txt+=open(self.header_file,'rU').read()
 
     #2 #######################################################################
     def put_version_info(self):
@@ -149,7 +153,7 @@
         self.full_banner_txt+="<MGVersion>\n"
         for key,pos in self.version_info.items():
             try:
-                self.full_banner_txt+="# "+key+' '*(25-len(key))+':'+file(pos,'rU').read()
+                self.full_banner_txt+="# "+key+' '*(25-len(key))+':'+open(pos,'rU').read()
             except IOError:
                 self.full_banner_txt+="# "+key+' '*(25-len(key))+':'
             if self.full_banner_txt[-1]!='\n':
@@ -170,7 +174,7 @@
             try:
                 self.full_banner_txt+=eval('self.mod_'+key+'(\"'+pos+'\")')
             except:
-                self.full_banner_txt+=file(pos,'rU').read()
+                self.full_banner_txt+=open(pos,'rU').read()
             self.full_banner_txt+="</"+key+">\n"
 
     #2 #######################################################################
@@ -179,7 +183,7 @@
         ##     |   self.input_file                                              ##
         """
         
-        self.full_banner_txt+=file(self.input_file,'rU').read()
+        self.full_banner_txt+=open(self.input_file,'rU').read()
         if self.full_banner_txt[-1]!='\n':
             self.full_banner_txt+='\n'
         self.full_banner_txt+="</LesHouchesEvents>\n"
@@ -192,7 +196,7 @@
         """
         suppress_empty_line=re.compile(r'''\n\s*\n''')
         self.full_banner_txt=suppress_empty_line.sub('\n',self.full_banner_txt)
-        output=file(self.output_file,'w')
+        output=open(self.output_file,'w')
         output.writelines(self.full_banner_txt)
         output.close()
 
@@ -255,11 +259,11 @@
         ##     |   modify (or original) run_card.dat                            ##
         """
 
-        text=file(inputfile).read()
+        text=open(inputfile).read()
         if os.path.isfile('./SubProcesses/randinit'):
             pattern_random=re.compile(r'''r=\s*(?P<val>\d*)''')
-            print file('./SubProcesses/randinit').read()
-            seed=pattern_random.search(file('./SubProcesses/randinit').read()).group('val')
+            print(open('./SubProcesses/randinit').read())
+            seed=pattern_random.search(open('./SubProcesses/randinit').read()).group('val')
             pattern_seed=re.compile(r'''0\s*=\s*iseed''')
             text=pattern_seed.sub(' '+seed+' = iseed',text)#
 
@@ -306,7 +310,7 @@
     def_pos.to_main()
     opt=sys.argv
     if len(opt)==1:
-        filename=raw_input('Enter file with events (in directory Events)')
+        filename=input('Enter file with events (in directory Events)')
     else:
         filename=opt[1]
     
 
=== modified file 'Template/MadWeight/Python/splitbanner.py'
--- Template/MadWeight/Python/splitbanner.py	2011-10-25 14:42:02 +0000
+++ Template/MadWeight/Python/splitbanner.py	2021-01-08 10:23:50 +0000
@@ -14,7 +14,10 @@
 ##   last-modif:24/01/09                                                ##
 ##                                                                      ##
 ##########################################################################
+from __future__ import absolute_import
+from __future__ import print_function
 import sys
+from six.moves import input
 # patch if symbolic directory replace by real file
 sys.path.append('./Source/MadWeight/Python')
 sys.path.append('../Source/MadWeight/Python')
@@ -52,7 +55,7 @@
 
     def write_card(self,pos):
 
-        print 'writing', pos
+        print('writing', pos)
         ff=open(pos,'w')
         for line in self.file:
             if self.pat_end.search(line):
@@ -62,7 +65,7 @@
                 ff.writelines(line)
 
     def FindCardName(self,tag):
-        if tag in dico.keys():
+        if tag in list(dico.keys()):
             if '.' not in dico[tag]:
                 return dico[tag]+'.dat'
             else:
@@ -77,9 +80,9 @@
     import MW_param
     MW_param.go_to_main_dir()
     opt=sys.argv
-    print opt
+    print(opt)
     if len(opt)<2:
-       name=raw_input('enter the run name to restore or the position of a banner file\n')
+       name=input('enter the run name to restore or the position of a banner file\n')
     else:
         name=opt[1]
 
@@ -91,15 +94,15 @@
         pos='/'.join([part for part in name.split('/')[:-1]])
         
     ban.split()
-    print 'done'
+    print('done')
     
     if len(opt)==3:
         import_lhco=opt[2]
     else:
-        import_lhco=raw_input('Do you want import the input.lhco file? (1/0)\n')
+        import_lhco=input('Do you want import the input.lhco file? (1/0)\n')
 
     
     if int(import_lhco):
         os.system('cp '+pos+'/input.lhco ./Events')
-        print 'replace Events/input.lhco by '+pos+'/input.lhco'
+        print('replace Events/input.lhco by '+pos+'/input.lhco')
 
 
=== modified file 'Template/MadWeight/Python/tests.py'
--- Template/MadWeight/Python/tests.py	2011-10-25 14:42:02 +0000
+++ Template/MadWeight/Python/tests.py	2021-01-08 10:23:50 +0000
@@ -20,7 +20,10 @@
 #########################################################################################################
 #  TEST #################################################################################################
 #########################################################################################################
+from __future__ import absolute_import
+from __future__ import print_function
 import sys
+from six.moves import input
 sys.path.append('./Source/MadWeight/Python')
 import filecmp
 import unittest
@@ -279,10 +282,10 @@
 
 
 
-print 'WARNING: This test modify this version of the Template. The correct way to launch it is to first made a copy of the Template'
-y =raw_input('Do you want to continue? [y/n]')
+print('WARNING: This test modify this version of the Template. The correct way to launch it is to first made a copy of the Template')
+y =input('Do you want to continue? [y/n]')
 if y == 'y':
     unittest.main()
 else:
-    print 'stopped'
+    print('stopped')
 
 
=== modified file 'Template/MadWeight/bin/madweight.py'
--- Template/MadWeight/bin/madweight.py	2013-03-30 07:45:35 +0000
+++ Template/MadWeight/bin/madweight.py	2021-01-08 10:23:50 +0000
@@ -15,6 +15,8 @@
 ################################################################################
 """ This is the main script in order to generate events in MadEvent """
 
+from __future__ import absolute_import
+from __future__ import print_function
 import sys
 import os
 import subprocess 
@@ -33,7 +35,7 @@
 
 # Check if optimize mode is (and should be) activated
 if __debug__ and (not os.path.exists(os.path.join(root_path,'../..', 'bin','create_release.py'))):
-    print 'launch in debug mode'
+    print('launch in debug mode')
     subprocess.call([sys.executable] + ['-O'] + sys.argv)
     sys.exit()
 
@@ -50,7 +52,7 @@
     try:
         import pyreadline as readline
     except:
-        print "For tab completion and history, install module readline."
+        print("For tab completion and history, install module readline.")
 else:
     import rlcompleter
 
@@ -81,7 +83,7 @@
    pass
 
 if __debug__:
-        print 'Running MG5 in debug mode'
+        print('Running MG5 in debug mode')
 
 
 def set_configuration():
@@ -163,10 +165,10 @@
             launch.run_cmd('quit')
         except:
             pass
-    except MW.AlreadyRunning, error:
+    except MW.AlreadyRunning as error:
         logging.error(str(error))
         sys.exit()
-    except Exception, error:
+    except Exception as error:
         if os.path.exists(pjoin(root_path, 'RunWeb')): 
             os.remove(pjoin(root_path, 'RunWeb'))
         raise 
 
=== modified file 'Template/MadWeight/mod_file/check_model.py'
--- Template/MadWeight/mod_file/check_model.py	2011-10-25 14:42:02 +0000
+++ Template/MadWeight/mod_file/check_model.py	2021-01-08 10:23:50 +0000
@@ -3,6 +3,8 @@
 #test
 
 #Extension
+from __future__ import absolute_import
+from __future__ import print_function
 import string
 import os
 import sys
@@ -29,7 +31,7 @@
     if 'bin' in list_dir:
         return
     else:
-        print 'script must be launched from main or bin/Source directory'
+        print('script must be launched from main or bin/Source directory')
         sys.exit()
             
 
 
=== modified file 'Template/MadWeight/mod_file/mod_file.py'
--- Template/MadWeight/mod_file/mod_file.py	2011-10-25 14:42:02 +0000
+++ Template/MadWeight/mod_file/mod_file.py	2021-01-08 10:23:50 +0000
@@ -1,11 +1,14 @@
 #!/usr/bin/env python
 
 #Extension
+from __future__ import absolute_import
+from __future__ import print_function
 import string
 import os
 import sys
 import time
 import re
+from six.moves import range
 
 ############################################################################
 ##                          PARAMETER
@@ -21,11 +24,11 @@
     
 #SECURITY POINT
     if type(list_input)!=list:
-        print 'subroutine need at least two files for fuse'
+        print('subroutine need at least two files for fuse')
         return 0
     
     if output in list_input:
-        print 'output file cann\'t have the same name than an input file'
+        print('output file cann\'t have the same name than an input file')
         return 0
 
 #initialisation
@@ -179,7 +182,7 @@
 
     file_in.close()
     file_out.close()
-    print 'no need of mod_photon_flux'
+    print('no need of mod_photon_flux')
     sys.exit()
     return
 
 
=== modified file 'Template/NLO/Cards/run_card.dat'
--- Template/NLO/Cards/run_card.dat	2018-04-25 06:47:06 +0000
+++ Template/NLO/Cards/run_card.dat	2021-01-08 10:23:50 +0000
@@ -179,7 +179,8 @@
   %(pt_max_pdg)s = pt_max_pdg ! Max pT for a massive particle
   %(mxx_min_pdg)s = mxx_min_pdg ! inv. mass for any pair of (anti)particles
 #***********************************************************************
-# For aMCfast+APPLGRID use in PDF fitting (http://amcfast.hepforge.org)*
+# Use PineAPPL to generate PDF-independent fast-interpolationuse       *
+# grids (https://zenodo.org/record/3992765#.X2EWy5MzbVo)               *
 #***********************************************************************
- %(iappl)s = iappl ! aMCfast switch (0=OFF, 1=prepare grids, 2=fill grids)
+ %(pineappl)s = pineappl ! PineAPPL switch 
 #***********************************************************************
 
=== modified file 'Template/NLO/FixedOrderAnalysis/HwU.f'
--- Template/NLO/FixedOrderAnalysis/HwU.f	2017-06-21 11:53:16 +0000
+++ Template/NLO/FixedOrderAnalysis/HwU.f	2021-01-08 10:23:50 +0000
@@ -10,16 +10,29 @@
 C                                                                          C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
+      module HwU_wgts_info_len
+         use iso_c_binding
+         integer, parameter :: wgts_info_len=80 
+         contains 
+             integer function get_wgts_info_len() bind(c,name="get_wgts_info_len")
+                 get_wgts_info_len = wgts_info_len
+                 return
+             end function get_wgts_info_len
+
+      end module HwU_wgts_info_len
+
+
 c The module contains effectively the common block with allocatable
 c variables (something not possible in old fortran version)
       module HwU_variables
+         use HwU_wgts_info_len
          implicit none
          integer :: max_plots,max_points,max_bins,nwgts,np
          integer :: error_estimation=3
          logical, allocatable :: booked(:)
          integer, allocatable :: nbin(:),histi(:,:),p_bin(:),p_label(:)
          character(len=50), allocatable :: title(:)
-         character(len=50), allocatable :: wgts_info(:)
+         character(len=wgts_info_len), allocatable :: wgts_info(:)
          double precision, allocatable :: histy(:,:,:),histy_acc(:,:,:)
      $        ,histy2(:,:),histy_err(:,:),histxl(:,:),histxm(:,:)
      $        ,step(:),p_wgts(:,:)
@@ -34,6 +47,16 @@
       implicit none
       integer i,nweights
       character*(*) wgt_info(*)
+
+c     PineAPPL commons
+      include "reweight_pineappl.inc"
+      include "pineappl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
+
+C     Initialize the number of bins of the aMCfast grids
+      if(pineappl) appl_obs_nbins = 0
+
       call HwU_deallocate_all
       max_plots=0
       max_points=0
@@ -85,6 +108,43 @@
       integer label,nbin_l,i,j
       character*(*) title_l
       double precision xmin,xmax
+c     PineAPPL commons
+      include "reweight_pineappl.inc"
+      include "pineappl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
+      double precision del
+
+c     Initialize the grids only if the switch "pineappl" is set to True
+c     and if the title  does not contain the word "Born". 
+      if(pineappl.and.index(title_l,"Born").eq.0)then
+c     Observable parameters
+c     Compute number of bins and edges only if they have not been given by the user.
+         if(appl_obs_nbins.eq.0)then
+            appl_obs_nbins = nbin_l
+            ! bin width
+            del = (xmax - xmin) / nbin_l
+c     compute bin edges
+            do i=0,appl_obs_nbins
+               appl_obs_bins(i) = xmin + i * del
+            enddo
+         endif
+         appl_obs_min = appl_obs_bins(0)
+         appl_obs_max = appl_obs_bins(appl_obs_nbins)
+         if(abs(appl_obs_max-xmax).gt.0.00000001d0)then
+            write(*,*) 'PineAPPL Histogram: ', 
+     1                 'Change of the upper limit:',xmax,'-->',
+     2                  appl_obs_max
+         endif
+c     Initialize PineAPPL routines
+         call APPL_init
+c     Keep track of the position of this histogram
+         nh_obs = nh_obs + 1
+         ih_obs(nh_obs) = label
+c     Reset number of bins to zero
+         appl_obs_nbins = 0
+      endif
+
 c     Allocate space for new histograms if needed      
       call HwU_allocate_histo(label,nbin_l)
 c     Setup the histogram
@@ -119,6 +179,24 @@
       implicit none
       integer label,i,j,bin
       double precision x, wgts(*)
+
+c     PineAPPL commons
+      include "reweight_pineappl.inc"
+      include "pineappl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
+      if(pineappl)then
+         do j=1,nh_obs
+            if(label.eq.ih_obs(j))then
+               appl_obs_num   = j
+               appl_obs_histo = x
+c     Fill the reference PineAPPL histograms
+c     Fill the PineAPPL files
+               call APPL_fill
+            endif
+         enddo
+      endif
+
 c     If central weight is zero do not add this point.
       if (wgts(1).eq.0d0) return
 c     Check if point is within plotting range
 
=== added file 'Template/NLO/FixedOrderAnalysis/analysis_td_hw_orders.f'
--- Template/NLO/FixedOrderAnalysis/analysis_td_hw_orders.f	1970-01-01 00:00:00 +0000
+++ Template/NLO/FixedOrderAnalysis/analysis_td_hw_orders.f	2021-01-08 10:23:50 +0000
@@ -0,0 +1,190 @@
+c
+c This file contains the default histograms for fixed order runs: it
+c only plots the total rate as an example. It can be used as a template
+c to make distributions for other observables.
+c
+c This uses the hbook package and generates histograms in the top-drawer
+c format. This format is human-readable. After running, the histograms
+c can be found in the Events/run_XX/ directory.
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_begin(nwgt,weights_info)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c This subroutine is called once at the start of each run. Here the
+c histograms should be declared. 
+c
+c Declare the histograms using 'bookup'.
+c     o) The first argument is an integer that labels the histogram. In
+c     the analysis_end and analysis_fill subroutines this label is used
+c     to keep track of the histogram. The label should be a number
+c     between 1 and NPLOTS/4=5000 (can be increased in dbook.inc).
+c     o) The second argument is a string that will apear above the
+c     histogram. Do not use brackets "(" or ")" inside this string.
+c     o) The third, forth and fifth arguments are the bin size, the
+c     lower edge of the first bin and the upper edge of the last
+c     bin. There is a maximum of 100 bins per histogram.
+c     o) When including scale and/or PDF uncertainties, declare a
+c     histogram for each weight, and compute the uncertainties from the
+c     final set of histograms
+c
+      implicit none
+c When including scale and/or PDF uncertainties the total number of
+c weights considered is nwgt
+      integer nwgt
+c In the weights_info, there is an text string that explains what each
+c weight will mean. The size of this array of strings is equal to nwgt.
+      character*(*) weights_info(*)
+c Local variables
+      integer kk,l,nwgt_analysis
+      common/c_analysis/nwgt_analysis
+c Initialize the histogramming package (hbook):
+      call inihist
+c Fill the c_analysis common block with the number of weights that will
+c be computed
+      nwgt_analysis=nwgt
+c
+c loop over all the weights that are computed (depends on run_card
+c parameters do_rwgt_scale and do_rwgt_pdf):
+      do kk=1,nwgt_analysis
+c make sure that there is a separate histogram initialized for each
+c weight
+         l=(kk-1)*4
+c declare (i.e. book) the histograms
+         call bookup(l+1,'total rate      '//weights_info(kk),
+     &        1.0d0,0.5d0,5.5d0)
+         call bookup(l+2,'total rate Born '//weights_info(kk),
+     &        1.0d0,0.5d0,5.5d0)
+         call bookup(l+3,'total rate QCD'//weights_info(kk),
+     &        1.0d0,0.5d0,5.5d0)
+         call bookup(l+4,'total rate EW '//weights_info(kk),
+     &        1.0d0,0.5d0,5.5d0)
+      enddo
+      return
+      end
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_end(xnorm)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c This subroutine is called once at the end of the run. Here the
+c histograms are written to disk. Note that this is done for each
+c integration channel separately. There is an external script that will
+c read the top drawer files in each of the integration channels and
+c combines them by summing all the bins in a final single top-drawer
+c file to be put in the Events/run_XX directory.
+c      
+c The histograms are put in a format to be written to file. Use the
+c multitop() subroutine.
+c     o) The first argument is the histogram label
+c     o) The second and third arguments are not used (keep them to the
+c     default 3,2)
+c     o) Fourth argument is the label for the x-axis. Do not use
+c     brackets "(" or ")" inside this string.
+c     o) Fifth argument is the y-axis
+c     o) Final argument declares if the y-axis should be a linear 'LIN'
+c     or logarithmic 'LOG' scale.
+      implicit none
+      character*14 ytit
+      double precision xnorm
+      integer i
+c Local variables
+      integer kk,l,nwgt_analysis
+      common/c_analysis/nwgt_analysis
+c This defines NPLOTS:
+      include 'dbook.inc'
+c Open the topdrawer file to write the histograms (do not remove the
+c next to calls)
+      call open_topdrawer_file
+      call mclear
+c Do not touch the folloing 4 lines. These lines make sure that the
+c histograms will have the correct overall normalisation: cross section
+c (in pb) per bin.
+      do i=1,NPLOTS
+         call mopera(i,'+',i,i,xnorm,0.d0)
+         call mfinal(i)
+      enddo
+      ytit='sigma per bin '
+c Loop over the plots that are declared with bookup
+      do kk=1,nwgt_analysis
+         l=(kk-1)*4
+         do i=1,8
+c convert them to a format suitable for writing
+            call multitop(l+i,3,2,'total rate',ytit,'LIN')
+         enddo
+      enddo
+c Close the topdrawer file (Do not remove this call)
+      call close_topdrawer_file
+      return
+      end
+
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c This subroutine is called for each n-body and (n+1)-body configuration
+c that passes the generation cuts. Here the histrograms are filled.
+      implicit none
+c This includes the 'nexternal' parameter that labels the number of
+c particles in the (n+1)-body process
+      include 'nexternal.inc'
+c This is an array which is '-1' for initial state and '1' for final
+c state particles
+      integer istatus(nexternal)
+c This is an array with (simplified) PDG codes for the particles. Note
+c that channels that are combined (i.e. they have the same matrix
+c elements) are given only 1 set of PDG codes. This means, e.g., that
+c when using a 5-flavour scheme calculation (massless b quark), no
+c b-tagging can be applied.
+      integer iPDG(nexternal)
+c The array of the momenta and masses of the initial and final state
+c particles in the lab frame. The format is "E, px, py, pz, mass", while
+c the second dimension loops over the particles in the process. Note
+c that these are the (n+1)-body particles; for the n-body there is one
+c momenta equal to all zero's (this is not necessarily the last particle
+c in the list). If one uses IR-safe obserables only, there should be no
+c difficulty in using this.
+      double precision p(0:4,nexternal)
+c The weight of the current phase-space point is wgts(1). If scale
+c and/or PDF uncertainties are included through reweighting, the rest of
+c the array contains the list of weights in the same order as described
+c by the weigths_info strings in analysis_begin
+      double precision wgts(*)
+c The ibody variable is:
+c     ibody=1 : (n+1)-body contribution
+c     ibody=2 : n-body contribution (excluding the Born)
+c     ibody=3 : Born contribution
+c The histograms need to be filled for all these contribution to get the
+c physics NLO results. (Note that the adaptive phase-space integration
+c is optimized using the sum of the contributions, therefore plotting
+c them separately might lead to larger than expected statistical
+c fluctuations).
+      integer ibody
+c local variables
+      double precision wgt,var
+      integer kk,l,nwgt_analysis
+      common/c_analysis/nwgt_analysis
+      ! stuff for plotting the different splitorders
+      integer orders_tag_plot
+      common /corderstagplot/ orders_tag_plot
+c
+c Fill the histograms here using a call to the mfill() subroutine. The
+c first argument is the histogram label, the second is the numerical
+c value of the variable to plot for the current phase-space point and
+c the final argument is the weight of the current phase-space point.
+      var=1d0
+c loop over all the weights that are computed (depends on run_card
+c parameters do_rwgt_scale and do_rwgt_pdf):
+      do kk=1,nwgt_analysis
+         wgt=wgts(kk)
+         l=(kk-1)*4
+c always fill the total rate
+         call mfill(l+1,var,wgt)
+c only fill the total rate for the Born when ibody=3
+         if (ibody.eq.3) call mfill(l+2,var,wgt)
+C fill the different partonic subchannels         
+         if (orders_tag_plot.eq.402) call mfill(l+3,var,wgt)
+         if (orders_tag_plot.eq.600) call mfill(l+4,var,wgt)
+      enddo
+      return
+      end
 
=== modified file 'Template/NLO/FixedOrderAnalysis/dbook.f'
--- Template/NLO/FixedOrderAnalysis/dbook.f	2015-03-04 10:22:46 +0000
+++ Template/NLO/FixedOrderAnalysis/dbook.f	2021-01-08 10:23:50 +0000
@@ -124,11 +124,11 @@
 c     initialization
       IMPLICIT NONE
       include 'dbook.inc'
-c     APPLgrid commons
-      include "reweight_appl.inc"
-      include "appl_common.inc"
-      integer iappl
-      common /for_applgrid/ iappl
+c     PineAPPL commons
+      include "reweight_pineappl.inc"
+      include "pineappl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
 c
 c     LOCAL
 c
@@ -139,8 +139,8 @@
       NHIST=0
       DO 1, I=1,NPLOTS             
    1  BOOK(I)=' NO'
-C     Initialize the number of bins of the aMCfast grids
-      if(iappl.ne.0) appl_obs_nbins = 0
+C     Initialize the number of bins of the amcblast grids
+      if(pineappl) appl_obs_nbins = 0
       END  
  
 
@@ -1001,20 +1001,19 @@
       implicit none
       integer n,n_by4
       double precision var,www
-c     APPLgrid commons
-      include "reweight_appl.inc"
-      include "appl_common.inc"
-      integer iappl
-      common /for_applgrid/ iappl
+c     PineAPPL commons
+      include "reweight_pineappl.inc"
+      include "pineappl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
       integer j
-      if(iappl.ne.0)then
+      if(pineappl)then
          do j=1,nh_obs
             if(n.eq.ih_obs(j))then
                appl_obs_num   = j
                appl_obs_histo = var 
-c     Fill the reference APPLgrid histograms
-               call APPL_fill_ref
-c     Fill the APPLgrid files
+c     Fill the reference PineAPPL histograms
+c     Fill the PineAPPL files
                call APPL_fill
             endif
          enddo
@@ -1029,15 +1028,15 @@
       integer i,n,n_by4
       character*(*) string
       double precision del,xl,xu
-c     APPLgrid commons
-      include "reweight_appl.inc"
-      include "appl_common.inc"
-      integer iappl
-      common /for_applgrid/ iappl
-c     Initialize the grids only if the switch "iappl" is different from zero
+c     PineAPPL commons
+      include "reweight_pineappl.inc"
+      include "pineappl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
+c     Initialize the grids only if the switch "pineappl" is different from zero
 c     and if the title string containes the word "central" and does not contain
 c     the word "Born". 
-      if(iappl.ne.0.and.index(string,"central").ne.0.and.
+      if(pineappl.and.index(string,"central").ne.0.and.
      1                  index(string,"Born").eq.0)then
 c     Observable parameters
 c     Compute number of bins and edges only if they have not been given by the user.
@@ -1051,11 +1050,11 @@
          appl_obs_min = appl_obs_bins(0)
          appl_obs_max = appl_obs_bins(appl_obs_nbins)
          if(abs(appl_obs_max-xu).gt.0.00000001d0)then
-            write(*,*) 'APPLgrid Histogram: ', 
+            write(*,*) 'PineAPPL Histogram: ', 
      1                 'Change of the upper limit:',xu,'-->',
      2                  appl_obs_max
          endif
-c     Initialize APPLgrid routines
+c     Initialize PineAPPL routines
          call APPL_init
 c     Keep track of the position of this histogram
          nh_obs = nh_obs + 1
@@ -1072,17 +1071,16 @@
       implicit none
       integer n,n_by4,m_by4
       character*(*) string1,string2,string3
-c     APPLgrid commons
-      include "reweight_appl.inc"
-      include "appl_common.inc"
-      integer iappl
-      common /for_applgrid/ iappl
+c     PineAPPL commons
+      include "reweight_pineappl.inc"
+      include "pineappl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
       integer j
-      if(iappl.ne.0)then
+      if(pineappl)then
          do j=1,nh_obs
             if(n.eq.ih_obs(j))then
                appl_obs_num = j
-               call APPL_fill_ref_out
                call APPL_term
             endif
          enddo
@@ -1098,17 +1096,16 @@
       implicit none
       integer n,n_by4,lr,lh,m_by4
       character*(*) string1,string2,string3
-c     APPLgrid commons
-      include "reweight_appl.inc"
-      include "appl_common.inc"
-      integer iappl
-      common /for_applgrid/ iappl
+c     PineAPPL commons
+      include "reweight_pineappl.inc"
+      include "pineappl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
       integer j
-      if(iappl.ne.0)then
+      if(pineappl)then
          do j=1,nh_obs
             if(n.eq.ih_obs(j))then
                appl_obs_num = j
-               call APPL_fill_ref_out
                call APPL_term
             endif
          enddo
 
=== modified file 'Template/NLO/FixedOrderAnalysis/rbook_fe8.f'
--- Template/NLO/FixedOrderAnalysis/rbook_fe8.f	2014-06-30 10:49:23 +0000
+++ Template/NLO/FixedOrderAnalysis/rbook_fe8.f	2021-01-08 10:23:50 +0000
@@ -95,24 +95,24 @@
        integer id,nbin
        double precision xbin,xlo,xhi
 
-c      APPLgrid commons
-       include "reweight_appl.inc"
-       include "appl_common.inc"
-       integer iappl
-       common /for_applgrid/ iappl
+c      PineAPPL commons
+       include "reweight_pineappl.inc"
+       include "pineappl_common.inc"
+       logical pineappl
+       common /for_pineappl/ pineappl
 
        nbin = int((xhi-xlo)/(xbin*0.9999d0))
 
-c      Initialize the grids only if the switch "iappl" is different from zero
+c      Initialize the grids only if the switch "pineappl" is different from zero
 c      and if the title string containes the word "central" and does not contain
 c      the word "Born".
-       if(iappl.ne.0.and.index(name,"central").ne.0.and.
+       if(pineappl.and.index(name,"central").ne.0.and.
      1                   index(name,"Born").eq.0)then
 c      Observable parameters
           appl_obs_nbins = nbin
           appl_obs_min   = xlo
           appl_obs_max   = xhi
-c      Initialize APPLgrid routines
+c      Initialize PineAPPL routines
           call APPL_init
 c      Keep track of the position of this histogram
           nh_obs = nh_obs + 1
@@ -150,21 +150,20 @@
        integer ihisto
        double precision xval,wgt
 
-c      APPLgrid commons
-       include "reweight_appl.inc"
-       include "appl_common.inc"
-       integer iappl
-       common /for_applgrid/ iappl
+c      PineAPPL commons
+       include "reweight_pineappl.inc"
+       include "pineappl_common.inc"
+       logical pineappl
+       common /for_pineappl/ pineappl
        integer j
 c
-       if(iappl.ne.0)then
+       if(pineappl)then
           do j=1,nh_obs
              if(ihisto.eq.ih_obs(j))then
                 appl_obs_num   = j
                 appl_obs_histo = xval
-c      Fill the reference APPLgrid histograms
-                call APPL_fill_ref
-c      Fill the APPLgrid files
+c      Fill the reference PineAPPL histograms
+c      Fill the PineAPPL files
                 call APPL_fill
              endif
           enddo
@@ -198,17 +197,16 @@
       character*(*) oper
       double precision x,y
 
-c     APPLgrid commons
-      include "reweight_appl.inc"
-      include "appl_common.inc"
-      integer iappl
-      common /for_applgrid/ iappl
+c     PineAPPL commons
+      include "reweight_pineappl.inc"
+      include "pineappl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
       integer j
-      if(iappl.ne.0)then
+      if(pineappl)then
          do j=1,nh_obs
             if(ih1.eq.ih_obs(j))then
                appl_obs_num = j
-               call APPL_fill_ref_out
                call APPL_term
             endif
          enddo
 
=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_V.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_V.f	2017-02-02 08:59:42 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_V.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE HWABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
       include 'reweight0.inc'
       real * 8 xm0,gam,xmlow,xmupp,bin
@@ -49,7 +50,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_h.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_h.f	2016-03-21 08:45:42 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_h.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE HWABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
       include 'reweight0.inc'
 c
@@ -46,7 +47,7 @@
       integer nwgt,max_weight,nwgt_analysis,kk,l
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_hjj.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_hjj.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_hjj.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE HWABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
       include 'reweight0.inc'
 c
@@ -51,7 +52,7 @@
       common /to_veto_hist/vetomin,vetomax,nbinveto
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_leptons.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_leptons.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_leptons.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE HWABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
       include 'reweight0.inc'
 c
@@ -60,7 +61,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_lplm.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_lplm.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_lplm.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE HWABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
       include 'reweight0.inc'
       real * 8 bin,xmi,xms,pi
@@ -48,7 +49,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_lvl.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_lvl.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_lvl.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE HWABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
       include 'reweight0.inc'
       real * 8 bin,xmi,xms,pi
@@ -48,7 +49,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_taptam.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_taptam.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_taptam.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE HWABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
       include 'reweight0.inc'
 c
@@ -46,7 +47,7 @@
       integer nwgt,max_weight,nwgt_analysis,kk,l
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_tj.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_tj.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_tj.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE HWABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
       include 'reweight0.inc'
       REAL*8 pi
@@ -46,7 +47,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_ttx.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_ttx.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_ttx.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE HWABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
       include 'reweight0.inc'
       REAL*8 pi
@@ -46,7 +47,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_ttx_v2.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_ttx_v2.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_ttx_v2.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE HWABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
       include 'reweight0.inc'
       REAL*8 pi
@@ -46,7 +47,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_rates.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_rates.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_rates.f	2021-01-08 10:23:50 +0000
@@ -23,6 +23,7 @@
       SUBROUTINE HWABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
       include 'reweight0.inc'
 c
@@ -42,7 +43,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/Makefile_MadFKS'
--- Template/NLO/MCatNLO/Makefile_MadFKS	2018-07-06 13:15:54 +0000
+++ Template/NLO/MCatNLO/Makefile_MadFKS	2021-01-08 10:23:50 +0000
@@ -43,10 +43,10 @@
 PDFFILES=$(OBJDIR)/mcatnlo_pdftomlm.o $(OBJDIR)/mcatnlo_libofpdf.o $(OBJDIR)/dummies.o
 CPDFFILES=$(OBJDIR)/mcatnlo_mlmtopdf.o $(OBJDIR)/dummies.o
 LPDFFILES=$(OBJDIR)/mcatnlo_mlmtolha.o $(OBJDIR)/dummies.o
-HWFILES=$(OBJDIR)/$(HERWIGVER) $(OBJDIR)/madfks_hwdriver.o $(OBJDIR)/madfks_hwlhin.o \
-        $(OBJDIR)/mcatnlo_str.o $(HWUTI)
-PYFILES=$(OBJDIR)/$(PYTHIAVER) $(OBJDIR)/madfks_pydriver.o $(OBJDIR)/madfks_pylhin.o \
-        $(OBJDIR)/mcatnlo_str.o $(PYUTI)
+HWFILES=$(HWUTI) $(OBJDIR)/$(HERWIGVER) $(OBJDIR)/madfks_hwdriver.o $(OBJDIR)/madfks_hwlhin.o \
+        $(OBJDIR)/mcatnlo_str.o 
+PYFILES=$(PYUTI) $(OBJDIR)/$(PYTHIAVER) $(OBJDIR)/madfks_pydriver.o $(OBJDIR)/madfks_pylhin.o \
+        $(OBJDIR)/mcatnlo_str.o
 
 
 HW_EXE_DEFAULT : $(HWFILES) $(PDFFILES)
 
=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_V.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_V.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_V.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG(nnn,wwwi)
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HEPMC.INC'
       include 'reweight0.inc'
       real * 8 xm0,gam,xmlow,xmupp,bin
@@ -48,7 +49,7 @@
       data cc/'     ','Born '/
       integer nwgt_analysis
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
      $     ,wwwi(max_weight_shower)
 c
       do i=1,nnn
 
=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_h.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_h.f	2016-07-21 14:37:54 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_h.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG(nnn,wwwi)
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HEPMC.INC'
       include 'reweight0.inc'
       REAL*8 pi
@@ -45,7 +46,7 @@
       data HwUtype/'|T@NOCUT','|T@CUT  '/
       integer nwgt_analysis
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
      $     ,wwwi(max_weight_shower)
 c
       do i=1,nnn
 
=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_hjj.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_hjj.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_hjj.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG(nnn,wwwi)
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HEPMC.INC'
       include 'reweight0.inc'
       REAL*8 pi
@@ -45,7 +46,7 @@
       data HwUtype/'|T@NOCUT','|T@CUT  '/
       integer nwgt_analysis
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
      $     ,wwwi(max_weight_shower)
       character*8 cc(2)
       data cc/'        ','vbfcuts '/
 
=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_leptons.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_leptons.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_leptons.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG(nnn,wwwi)
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HEPMC.INC'
       include 'reweight0.inc'
 c
@@ -43,7 +44,7 @@
       integer j,jpr
       integer nwgt_analysis
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
      $     ,wwwi(max_weight_shower)
       integer nsingle,ncorr,nlepton,nplots,ncuts
       common/cplots/nsingle,ncorr,nlepton,nplots,ncuts
 
=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_lplm.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_lplm.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_lplm.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG(nnn,wwwi)
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HEPMC.INC'
       include 'reweight0.inc'
       real * 8 bin,xmi,xms,pi
@@ -47,7 +48,7 @@
       data cc/'     ','Born '/
       integer nwgt_analysis
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
      $     ,wwwi(max_weight_shower)
 c
       do i=1,nnn
 
=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_lvl.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_lvl.f	2016-12-07 21:05:19 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_lvl.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG(nnn,wwwi)
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HEPMC.INC'
       include 'reweight0.inc'
       real * 8 bin,xmi,xms,pi
@@ -47,7 +48,7 @@
       data cc/'     ','Born '/
       integer nwgt_analysis
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
      $     ,wwwi(max_weight_shower)
 c
       do i=1,nnn
 
=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_taptam.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_taptam.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_taptam.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG(nnn,wwwi)
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HEPMC.INC'
       include 'reweight0.inc'
       integer j,kk,l,i,nnn
@@ -43,7 +44,7 @@
       data HwUtype/'|T@NOCUT','|T@CUT  '/
       integer nwgt_analysis
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
      $     ,wwwi(max_weight_shower)
 c
       do i=1,nnn
 
=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_tj.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_tj.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_tj.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG(nnn,wwwi)
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HEPMC.INC'
       include 'reweight0.inc'
       REAL*8 pi
@@ -45,7 +46,7 @@
       data HwUtype/'|T@NOCUT','|T@CUT  '/
       integer nwgt_analysis
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
      $     ,wwwi(max_weight_shower)
 c
       do i=1,nnn
 
=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_ttx.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_ttx.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_ttx.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG(nnn,wwwi)
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HEPMC.INC'
       include 'reweight0.inc'
       REAL*8 pi
@@ -45,7 +46,7 @@
       data HwUtype/'|T@NOCUT','|T@CUT  '/
       integer nwgt_analysis
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
      $     ,wwwi(max_weight_shower)
 c
       do i=1,nnn
 
=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_ttx_v2.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_ttx_v2.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_ttx_v2.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG(nnn,wwwi)
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HEPMC.INC'
       include 'reweight0.inc'
       REAL*8 pi
@@ -45,7 +46,7 @@
       data HwUtype/'|T@NOCUT','|T@CUT  '/
       integer nwgt_analysis
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
      $     ,wwwi(max_weight_shower)
 c
       do i=1,nnn
 
=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_rates.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_rates.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_rates.f	2021-01-08 10:23:50 +0000
@@ -23,6 +23,7 @@
       SUBROUTINE PYABEG(nnn,wwwi)
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HEPMC.INC'
       include 'reweight0.inc'
 c
@@ -41,7 +42,7 @@
       data cc/'     ','Born '/
       integer nwgt_analysis
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
      $     ,wwwi(max_weight_shower)
 c
       do i=1,nnn
 
=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_V.f'
--- Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_V.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_V.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       implicit double precision(a-h, o-z)
       implicit integer(i-n)
       common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
@@ -51,7 +52,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_h.f'
--- Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_h.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_h.f	2021-01-08 10:23:50 +0000
@@ -28,6 +28,7 @@
       SUBROUTINE PYABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       implicit none
       include 'reweight0.inc'
       REAL*8 pi
@@ -48,7 +49,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c
 c Initialize histograms
 
=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_hjj.f'
--- Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_hjj.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_hjj.f	2021-01-08 10:23:50 +0000
@@ -28,6 +28,7 @@
       SUBROUTINE PYABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       implicit none
       include 'reweight0.inc'
 c
@@ -52,7 +53,7 @@
       common /to_veto_hist/vetomin,vetomax,nbinveto
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c
 c Initialize histograms
 
=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_leptons.f'
--- Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_leptons.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_leptons.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       implicit none
       include 'reweight0.inc'
 c
@@ -59,7 +60,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_lplm.f'
--- Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_lplm.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_lplm.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       implicit none
       include 'reweight0.inc'
       real * 8 bin,xmi,xms,pi
@@ -48,7 +49,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_lvl.f'
--- Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_lvl.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_lvl.f	2021-01-08 10:23:50 +0000
@@ -27,6 +27,7 @@
       SUBROUTINE PYABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       implicit double precision(a-h, o-z)
       implicit integer(i-n)
       common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
@@ -50,7 +51,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_taptam.f'
--- Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_taptam.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_taptam.f	2021-01-08 10:23:50 +0000
@@ -28,6 +28,7 @@
       SUBROUTINE PYABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       implicit none
       include 'reweight0.inc'
       integer j,kk,l,i
@@ -46,7 +47,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c
 c Initialize histograms
 
=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_tj.f'
--- Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_tj.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_tj.f	2021-01-08 10:23:50 +0000
@@ -28,6 +28,7 @@
       SUBROUTINE PYABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       implicit none
       include 'reweight0.inc'
       REAL*8 pi
@@ -48,7 +49,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c
 c Initialize histograms
 
=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_ttx.f'
--- Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_ttx.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_ttx.f	2021-01-08 10:23:50 +0000
@@ -28,6 +28,7 @@
       SUBROUTINE PYABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       implicit none
       include 'reweight0.inc'
       REAL*8 pi
@@ -48,7 +49,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c
 c Initialize histograms
 
=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_ttx_v2.f'
--- Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_ttx_v2.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_pp_ttx_v2.f	2021-01-08 10:23:50 +0000
@@ -28,6 +28,7 @@
       SUBROUTINE PYABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       implicit none
       include 'reweight0.inc'
       REAL*8 pi
@@ -48,7 +49,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c
 c Initialize histograms
 
=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_rates.f'
--- Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_rates.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/py6an_HwU_rates.f	2021-01-08 10:23:50 +0000
@@ -23,6 +23,7 @@
       SUBROUTINE PYABEG
 C     USER''S ROUTINE FOR INITIALIZATION
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       implicit none
       include 'reweight0.inc'
 c
@@ -42,7 +43,7 @@
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 c Initialize histograms
       call HwU_inithist(nwgt,weights_info)
 
=== modified file 'Template/NLO/MCatNLO/include/LHEFRead.h'
--- Template/NLO/MCatNLO/include/LHEFRead.h	2016-02-23 11:25:37 +0000
+++ Template/NLO/MCatNLO/include/LHEFRead.h	2021-01-08 10:23:50 +0000
@@ -2,6 +2,13 @@
 #include <iomanip>
 #include <string>
 
+
+extern "C" {
+    int get_wgts_info_len(void);
+}
+
+const int wgts_info_len_used = 80;
+
 class MyReader {
 
   private:
@@ -16,12 +23,22 @@
   LHEF::Reader reader;
 
 
-  void lhef_read_wgtsinfo_(int &cwgtinfo_nn, char (cwgtinfo_weights_info[1024][50])) {
+  void lhef_read_wgtsinfo_(int &cwgtinfo_nn, char (cwgtinfo_weights_info[1024][wgts_info_len_used])) {
+
+    // check that wgts_info_len_used > wgts_info_len from the fortran module
+    if (wgts_info_len_used != get_wgts_info_len()) {
+      std::cout << " wgts_info_len_used is different from the corresponding quantity in HwU_wgts_info_len\n"
+                << " Please set wgts_info_len_used in LHEFRead.h to "<< get_wgts_info_len() <<"\n"
+                << " Program stopped.\n";
+      exit(1);
+    }
 
     // Read header of event file
     std::stringstream hss;
     std::string hs;
-    sprintf(cwgtinfo_weights_info[0], "%50s","central value");
+    std::stringstream format_str;
+    format_str <<"%"<<wgts_info_len_used<<"s";
+    sprintf(cwgtinfo_weights_info[0],format_str.str().c_str(),"central value");
     cwgtinfo_nn=1;
     while (true){
       hss << reader.headerBlock;
@@ -37,7 +54,7 @@
 	  if (hs.find("<weight id") != std::string::npos) {
 	    std::string sRWGT = hs.substr(hs.find("'>")+2,hs.find("</w")-3);
 	    sRWGT = sRWGT.substr(0,sRWGT.find("<"));
-	    sprintf(cwgtinfo_weights_info[cwgtinfo_nn],"%50s",sRWGT.c_str());
+	    sprintf(cwgtinfo_weights_info[cwgtinfo_nn],format_str.str().c_str(),sRWGT.c_str());
 	    ++cwgtinfo_nn;
 	  }
 	}
 
=== modified file 'Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f'
--- Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f	2017-06-19 14:46:09 +0000
+++ Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f	2021-01-08 10:23:50 +0000
@@ -125,6 +125,7 @@
 C  Reads MC@NLO input headers and fills Les Houches run common HEPRUP
 C  Event file is written by MadFKS
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
       INCLUDE 'HERWIG65.INC'
 C--Les Houches Common Blocks
       INTEGER MAXPUP
@@ -150,7 +151,7 @@
       include 'reweight0.inc'
       integer nwgt
       common/cnwgt/nwgt
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
       
       nwgt=1
 
=== modified file 'Template/NLO/MCatNLO/srcPythia/madfks_pylhin.f'
--- Template/NLO/MCatNLO/srcPythia/madfks_pylhin.f	2017-06-19 14:46:09 +0000
+++ Template/NLO/MCatNLO/srcPythia/madfks_pylhin.f	2021-01-08 10:23:50 +0000
@@ -124,6 +124,7 @@
 C  Reads MC@NLO input headers and fills Les Houches run common HEPRUP
 C  Event file is written by MadFKS
 C----------------------------------------------------------------------
+      use HwU_wgts_info_len
 C--Les Houches Common Blocks
       INTEGER MAXPUP
       PARAMETER(MAXPUP=100)
@@ -149,7 +150,7 @@
       include 'reweight0.inc'
       integer nwgt
       common/cnwgt/nwgt
-      character*50 weights_info(max_weight_shower)
+      character*(wgts_info_len) weights_info(max_weight_shower)
       common/cwgtsinfo/weights_info
 C
       numscales=0
 
=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia8.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia8.cc	2017-07-25 14:13:14 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia8.cc	2021-01-08 10:23:50 +0000
@@ -21,7 +21,7 @@
 #define cevwgt cevwgt_
 
 extern "C" { 
-  void pyabeg_(int&,char(*)[50]);
+  void pyabeg_(int&,char(*)[wgts_info_len_used]);
   void pyaend_(double&);
   void pyanal_(int&,double(*));
 }
@@ -30,7 +30,7 @@
   Pythia pythia;
 
   int cwgtinfo_nn;
-  char cwgtinfo_weights_info[1024][50];
+  char cwgtinfo_weights_info[1024][wgts_info_len_used];
   double cwgt_ww[1024];
 
   string inputname="Pythia8.cmd",outputname="Pythia8.hep";
 
=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia82.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia82.cc	2018-05-01 13:46:54 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia82.cc	2021-01-08 10:23:50 +0000
@@ -22,7 +22,7 @@
 #define cevwgt cevwgt_
 
 extern "C" { 
-  void pyabeg_(int&,char(*)[50]);
+  void pyabeg_(int&,char(*)[wgts_info_len_used]);
   void pyaend_(double&);
   void pyanal_(int&,double(*));
 }
@@ -31,7 +31,7 @@
   Pythia pythia;
 
   int cwgtinfo_nn;
-  char cwgtinfo_weights_info[1024][50];
+  char cwgtinfo_weights_info[1024][wgts_info_len_used];
   double cwgt_ww[1024];
 
   string inputname="Pythia8.cmd",outputname="Pythia8.hep";
 
=== modified file 'Template/NLO/Source/PDF/makefile'
--- Template/NLO/Source/PDF/makefile	2018-04-16 14:08:47 +0000
+++ Template/NLO/Source/PDF/makefile	2021-01-08 10:23:50 +0000
@@ -13,16 +13,12 @@
 ifdef lhapdf
   ifeq ($(lhapdfversion),5)
     PDF         = pdfwrap_lhapdf.o pdf_lhapdf.o pdg2pdf_lhapdf.o opendata.o
-  else
-    ifdef applgrid
-      PDF         = pdfwrap_lhapdf.o pdf_lhapdf.o pdg2pdf_lhapdf.o opendata.o
-    else # lhapdf6
-      ifeq ($(lhapdfsubversion),1) # 6.1.X
-        PDF         = pdfwrap_lhapdf.o pdf_lhapdf6.o pdg2pdf_lhapdf6.o opendata.o
-      else # 6.2.X
-        CXXFLAGS+=-std=c++11    
-        PDF         = pdfwrap_lhapdf.o pdf_lhapdf62.o pdg2pdf_lhapdf6.o opendata.o
-      endif
+  else # lhapdf v6
+    ifeq ($(lhapdfsubversion),1) # 6.1.X
+      PDF         = pdfwrap_lhapdf.o pdf_lhapdf6.o pdg2pdf_lhapdf6.o opendata.o
+    else # 6.2.X
+      CXXFLAGS+=-std=c++11    
+      PDF         = pdfwrap_lhapdf.o pdf_lhapdf62.o pdg2pdf_lhapdf6.o opendata.o
     endif
   endif
 else
 
=== modified file 'Template/NLO/Source/extra_weights.f'
--- Template/NLO/Source/extra_weights.f	2020-07-28 10:05:22 +0000
+++ Template/NLO/Source/extra_weights.f	2021-01-08 10:23:50 +0000
@@ -8,7 +8,7 @@
          integer :: iwgtnumpartn,jwgtinfo,mexternal
      $        ,lhaPDFid(0:maxPDFsets),nmemPDF(maxPDFsets)
      $        ,dyn_scale(0:maxdynscales),n_ctr_found,n_mom_conf
-     $        ,QCD_power,orders_tag
+     $        ,QCD_power,orders_tag,amp_pos
          double precision :: wgtdegrem_xi,wgtdegrem_lxi,wgtdegrem_muF
      $        ,wgtnstmp,wgtwnstmpmuf,wgtwnstmpmur,wgtnstmp_avgvirt
      $        ,wgtref,scalevarR(0:maxscales),scalevarF(0:maxscales)
 
=== modified file 'Template/NLO/Source/make_opts.inc'
--- Template/NLO/Source/make_opts.inc	2018-04-16 14:08:47 +0000
+++ Template/NLO/Source/make_opts.inc	2021-01-08 10:23:50 +0000
@@ -111,13 +111,14 @@
   endif
 endif
 
-# Option APPLGrid
-ifneq ($(applgrid),)
-  APPLLIBS=$(shell applgrid-config --ldcflags) $(shell amcfast-config --ldflags)
-  applgrid_interface=appl_interface.o
+# Option PineAPPL
+ifneq ($(pineappl),)
+  APPLLIBS=$(shell pkg-config pineappl_capi --libs)
+  CXXFLAGS+=$(shell pkg-config pineappl_capi --cflags)
+  pineappl_interface=pineappl_interface.o
 else
   APPLLIBS=
-  applgrid_interface=appl_interface_dummy.o
+  pineappl_interface=pineappl_interface_dummy.o
 endif
 
 # Madloop
 
=== modified file 'Template/NLO/Source/setrun.f'
--- Template/NLO/Source/setrun.f	2018-04-16 14:08:47 +0000
+++ Template/NLO/Source/setrun.f	2021-01-08 10:23:50 +0000
@@ -29,8 +29,8 @@
       integer nevents
       character*7 event_norm
       common /event_normalisation/event_norm
-      integer iappl
-      common /for_applgrid/ iappl
+      logical pineappl
+      common /for_pineappl/ pineappl
       integer idum
       logical              fixed_order,nlo_ps
       common /c_fnlo_nlops/fixed_order,nlo_ps
 
=== modified file 'Template/NLO/SubProcesses/ajob_template'
--- Template/NLO/SubProcesses/ajob_template	2017-02-13 10:00:30 +0000
+++ Template/NLO/SubProcesses/ajob_template	2021-01-08 10:23:50 +0000
@@ -31,6 +31,10 @@
 runnumber=$3
 integration_step=$4
 
+FO_EXTRAPATHS=
+export LD_LIBRARY_PATH=$FO_EXTRAPATHS:$LD_LIBRARY_PATH
+export DYLD_LIBRARY_PATH=$FO_EXTRAPATHS:$DYLD_LIBRARY_PATH
+
 TAGTAGTAGTAGTAGTAGTAG for i in 1 ; do
     
     if [[ $run_mode == 'all' || $run_mode == 'born' ]] ; then
 
=== modified file 'Template/NLO/SubProcesses/cuts.f'
--- Template/NLO/SubProcesses/cuts.f	2020-06-21 21:31:33 +0000
+++ Template/NLO/SubProcesses/cuts.f	2021-01-08 10:23:50 +0000
@@ -175,7 +175,7 @@
 
 c If we do not require a mimimum jet energy, there's no need to apply
 c jet clustering and all that.
-      if (ptj.ne.0d0.or.ptgmin.ne.0d0) then
+      if (ptj.gt.0d0.or.ptgmin.gt.0d0) then
 c Put all (light) QCD partons in momentum array for jet clustering.
 c From the run_card.dat, maxjetflavor defines if b quark should be
 c considered here (via the logical variable 'is_a_jet').  nQCD becomes
@@ -284,7 +284,7 @@
             is_a_ph(i)=.false.
          endif
       enddo
-      if (ptgmin.ne.0d0) then
+      if (ptgmin.gt.0d0) then
          nph=0
          do j=nincoming+1,nexternal
             if (is_a_ph(j)) then
@@ -442,150 +442,6 @@
 
 
 
-
-      subroutine recombine_momenta(R, etaph, reco_l, reco_q, p_in, pdg_in, p_out, pdg_out)
-      implicit none
-      ! recombine photons with the closest fermion if the distance is
-      ! less than R and if the rapidity of photons is < etaph (etaph < 0
-      ! means no cut). Output a new set of momenta and pdgs corresponding
-      ! to the recombined particles. If recombination occurs the photon
-      ! disappears from the output particles
-      ! arguments
-      include 'nexternal.inc'
-      double precision R, etaph, p_in(0:4,nexternal), p_out(0:4,nexternal)
-      logical reco_l, reco_q
-      integer pdg_in(nexternal), pdg_out(nexternal)
-      ! local variables
-      integer nq, nl
-      integer id_ph
-      parameter (id_ph=22)
-      integer n_ph, i_ph
-      integer i,j
-      integer ifreco
-      double precision dreco, dthis
-      integer skip
-      logical is_light_charged_fermion
-      double precision R2_04, eta_04
-      ! 
-      integer times_reco
-      common/to_times_reco/ times_reco
-      ! reset everything
-      do j=1,nexternal
-        pdg_out(j)=0
-        do i=0,4
-          p_out(i,j)=0d0
-        enddo
-      enddo
-
-      ! check if we want to recombine with leptons
-      if (reco_l) then
-          nl = 3
-      else 
-          nl = 0
-      endif
-
-      ! check if we want to recombine with quarks
-      if (reco_q) then
-          nq = 5
-      else 
-          nq = 0
-      endif
-
-      ! count the photons
-      n_ph=0
-      do i=nincoming+1, nexternal
-        if (pdg_in(i).eq.id_ph.and.
-     $   (abs(eta_04(p_in(0,i))).lt.etaph.or.etaph.lt.0d0)) then
-            n_ph=n_ph+1
-            i_ph=i
-        endif
-      enddo
-      if (n_ph.eq.0 .or. (nl.eq.0 .and. nq.eq.0)) then
-        ! do nothing
-        do j=1,nexternal
-          pdg_out(j)=pdg_in(j)
-          do i=0,4
-            p_out(i,j)=p_in(i,j)
-          enddo
-        enddo
-        return
-      elseif (n_ph.eq.1) then
-        ! do nothing for initial states
-        do j=1,nincoming
-          pdg_out(j)=pdg_in(j)
-          do i=0,4
-            p_out(i,j)=p_in(i,j)
-          enddo
-        enddo
-        ! find the closest fermion to the photon
-        ifreco=0
-        dreco=R
-        if (i_ph.gt.0) then
-          do i = nincoming+1, nexternal
-            if (is_light_charged_fermion(pdg_in(i),nq,nl)) then
-              dthis=dsqrt(R2_04(p_in(0,i_ph),p_in(0,i)))
-              if (dthis.le.dreco) then
-                dreco=dthis
-                ifreco=i
-              endif
-            endif
-          enddo
-        endif
-        if (ifreco.eq.0) then
-        ! do nothing also for final states
-          do j=nincoming+1,nexternal
-            pdg_out(j)=pdg_in(j)
-            do i=0,4
-              p_out(i,j)=p_in(i,j)
-            enddo
-          enddo
-        else
-          times_reco=times_reco+1
-          skip=0
-          do j=nincoming+1,nexternal
-            if (j.ne.i_ph.and.j.ne.ifreco) then
-              pdg_out(j-skip)=pdg_in(j)
-              do i=0,4
-                p_out(i,j-skip)=p_in(i,j)
-              enddo
-            elseif (j.eq.ifreco) then
-              pdg_out(j-skip)=pdg_in(j)
-              do i=0,3
-                p_out(i,j-skip)=p_in(i,j)+p_in(i,i_ph)
-              enddo
-              p_out(4,j-skip)=p_in(4,j)
-            elseif (j.eq.i_ph) then
-              skip=skip+1
-            endif
-          enddo
-        endif
-      else
-        write(*,*) 'ERROR, too many photons', n_ph
-        stop 1
-      endif
-
-      return 
-      end
-
-
-      logical function is_light_charged_fermion(id, nf, nl)
-      implicit none
-      integer id, nf, nl
-      if (abs(id).le.nf) then
-          is_light_charged_fermion = .true.
-      elseif ((abs(id).eq.11.and.nl.ge.1).or.
-     $        (abs(id).eq.13.and.nl.ge.2).or.
-     $        (abs(id).eq.15.and.nl.ge.3)) then
-          is_light_charged_fermion = .true.
-      else
-          is_light_charged_fermion = .false.
-      endif
-      return
-      end
-
-
-
-
 C***************************************************************
 C***************************************************************
 C NO NEED TO CHANGE ANY OF THE FUNCTIONS BELOW
@@ -1128,3 +984,4 @@
       end
 
 
+
 
=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2020-03-26 15:16:08 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2021-01-08 10:23:50 +0000
@@ -67,9 +67,9 @@
       include "timing_variables.inc"
       real*4 tOther, tTot
 
-c applgrid
-      integer iappl
-      common /for_applgrid/ iappl
+c PineAPPL
+      logical pineappl
+      common /for_pineappl/ pineappl
 c stats for granny_is_res
       double precision deravg,derstd,dermax,xi_i_fks_ev_der_max
      &     ,y_ij_fks_ev_der_max
@@ -160,10 +160,10 @@
         stop
       endif
       write(*,*) "about to integrate ", ndim,ncalls0,itmax
-c APPLgrid
-      if (imode.eq.0) iappl=0 ! overwrite when starting completely fresh
-      if(iappl.ne.0) then
-         write(6,*) "Initializing aMCfast ..."
+c PineAPPL
+      if (imode.eq.0) pineappl=.False. ! overwrite when starting completely fresh
+      if(pineappl) then
+         write(6,*) "Initializing PineAPPL ..."
 c     Set flavor map, starting from all possible
 c     parton lumi configurations defined in initial_states_map.dat
          call setup_flavourmap
@@ -365,8 +365,9 @@
       common/ccalculatedBorn/calculatedBorn
       character*4      abrv
       common /to_abrv/ abrv
-      integer iappl
-      common /for_applgrid/ iappl
+c PineAPPL
+      logical pineappl
+      common /for_pineappl/ pineappl
       double precision       wgt_ME_born,wgt_ME_real
       common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
       integer ini_fin_fks_map(0:2,0:fks_configs)
@@ -386,9 +387,9 @@
          write (*,*) 'ERROR ifl not equal to zero in sigint',ifl
          stop 1
       endif
-      if (iappl.ne.0 .and. sum) then
-         write (*,*) 'WARNING: applgrid only possible '/
-     &        /'with MC over FKS directories',iappl,sum
+      if (pineappl .and. sum) then
+         write (*,*) 'WARNING: PineAPPL only possible '/
+     &        /'with MC over FKS directories',pineappl,sum
          write (*,*) 'Switching to MC over FKS directories'
          sum=.false.
       endif
@@ -501,13 +502,13 @@
          if (do_rwgt_pdf) call reweight_pdf
       endif
       
-      if (iappl.ne.0) then
+      if (pineappl) then
          if (sum) then
-            write (*,*) 'ERROR: applgrid only possible '/
-     &           /'with MC over FKS directories',iappl,sum
+            write (*,*) 'ERROR: PineAPPL only possible '/
+     &           /'with MC over FKS directories',pineappl,sum
             stop 1
          endif
-         call fill_applgrid_weights(vegas_wgt)
+         call fill_pineappl_weights(vegas_wgt)
       endif
 
 c Importance sampling for FKS configurations
 
=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2020-07-28 10:05:22 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2021-01-08 10:23:50 +0000
@@ -43,6 +43,7 @@
         wgtcpower=0d0
         if (cpower_pos.gt.0) wgtcpower=dble(orders(cpower_pos))
         orders_tag=get_orders_tag(orders)
+        amp_pos=iamp
         wgt1=amp_split(iamp)*f_b/g**(qcd_power)
         call add_wgt(2,orders,wgt1,0d0,0d0)
       enddo
@@ -62,6 +63,7 @@
         wgtcpower=0d0
         if (cpower_pos.gt.0) wgtcpower=dble(orders(cpower_pos))
         orders_tag=get_orders_tag(orders)
+        amp_pos=iamp
         wgt6f1=amp_split_6to5f(iamp)*f_b/g**(qcd_power)
         wgt6f2=amp_split_6to5f_mur(iamp)*f_b/g**(qcd_power)
         wgt6f3=amp_split_6to5f_muf(iamp)*f_b/g**(qcd_power)
@@ -150,8 +152,10 @@
         alphasbpow = orders(qcd_pos)/2
         if (niglu.ne.0 .or. alphasbpow.ne.0) then
           ! this contribution will end up with one extra power
-          ! of alpha_s
+          ! of alpha_s. Check that we are including the corresponding
+          ! order in the computation, otherwise skip the contribution
           orders(qcd_pos) = orders(qcd_pos) + 2
+          if (orders(qcd_pos).gt.nlo_orders(qcd_pos)) cycle
 
           amp_split_6to5f_muf(orders_to_amp_split_pos(orders)) = 
      &     alphas / 3d0 / pi * TF * dble(niglu) * amp_split(iamp)  
@@ -254,6 +258,7 @@
         wgtcpower=0d0
         if (cpower_pos.gt.0) wgtcpower=dble(orders(cpower_pos))
         orders_tag=get_orders_tag(orders)
+        amp_pos=iamp
         g22=g**(QCD_power)
         wgt1=amp_split_wgtnstmp(iamp)*f_nb/g22
         wgt2=amp_split_wgtwnstmpmur(iamp)*f_nb/g22
@@ -278,6 +283,7 @@
         wgtcpower=0d0
         if (cpower_pos.gt.0) wgtcpower=dble(orders(cpower_pos))
         orders_tag=get_orders_tag(orders)
+        amp_pos=iamp
         virt_wgt_mint(iamp)=amp_split_virt(iamp)*f_nb
         born_wgt_mint(iamp)=amp_split_born_for_virt(iamp)*f_nb
         wgt1=virt_wgt_mint(iamp)/g**(QCD_power)
@@ -322,6 +328,7 @@
         wgtcpower=0d0
         if (cpower_pos.gt.0) wgtcpower=dble(orders(cpower_pos))
         orders_tag=get_orders_tag(orders)
+        amp_pos=iamp
         wgt1=amp_split(iamp)*s_ev*f_r/g**(qcd_power)
         if (sudakov_damp.gt.0d0) then
           call add_wgt(1,orders,wgt1*sudakov_damp,0d0,0d0)
@@ -385,6 +392,7 @@
         wgtcpower=0d0
         if (cpower_pos.gt.0) wgtcpower=dble(orders(cpower_pos))
         orders_tag=get_orders_tag(orders)
+        amp_pos=iamp
         g22=g**(QCD_power)
         if (replace_MC_subt.gt.0d0) then
           wgt1=amp_split(iamp)*s_s/g22*replace_MC_subt
@@ -478,6 +486,7 @@
         wgtcpower=0d0
         if (cpower_pos.gt.0) wgtcpower=dble(orders(cpower_pos))
         orders_tag=get_orders_tag(orders)
+        amp_pos=iamp
         g22=g**(QCD_power)
         if (replace_MC_subt.gt.0d0) then
           wgt1=amp_split(iamp)*s_c/g22*replace_MC_subt
@@ -587,6 +596,7 @@
         wgtcpower=0d0
         if (cpower_pos.gt.0) wgtcpower=dble(orders(cpower_pos))
         orders_tag=get_orders_tag(orders)
+        amp_pos=iamp
         g22=g**(QCD_power)
         if (replace_MC_subt.gt.0d0) then
           wgt1=-amp_split(iamp)*s_sc/g22*replace_MC_subt
@@ -692,6 +702,7 @@
                 wgtcpower=0d0
                 if (cpower_pos.gt.0) wgtcpower=dble(orders(cpower_pos))
                 orders_tag=get_orders_tag(orders)
+                amp_pos=iamp
                 g22=g**(QCD_power)
                 wgt1=sevmc*f_MC_S*xlum_mc_fact*
      &               amp_split_xmcxsec(iamp,i)/g22
@@ -980,7 +991,7 @@
       include 'genps.inc'
       include 'timing_variables.inc'
       double precision pi,vegas_wgt
-      integer i
+      integer i,j
       logical firsttime
       data firsttime /.true./
       parameter (pi=3.1415926535897932385d0)
@@ -1002,16 +1013,17 @@
      &                  fkssymmetryfactorDeg,ngluons,nquarks,nphotons
       double precision      f_b,f_nb
       common /factor_nbody/ f_b,f_nb
-      integer iappl
-      common /for_applgrid/ iappl
+      logical pineappl
+      common /for_pineappl/ pineappl
       logical needrndec
       parameter (needrndec=.false.)
       real*8 ran2
       external ran2
       real*8 rndec(10)
       common/crndec/rndec
-      include "appl_common.inc" 
-      !!!include "orders.inc"
+      include "pineappl_common.inc" 
+      integer orders(nsplitorders)
+C
       call cpu_time(tBefore)
 c Random numbers to be used in the plotting routine: these numbers will
 c not change between events, counter events and n-body contributions.
@@ -1021,7 +1033,15 @@
          enddo
       endif
       if (firsttime) then
-         if (iappl.ne.0) appl_amp_split_size = amp_split_size
+         if (pineappl) then 
+         ! PineAPPL stuff
+         appl_amp_split_size = amp_split_size
+           do j=1,amp_split_size
+             call amp_split_pos_to_orders(j, orders)
+             appl_qcdpower(j) = orders(qcd_pos)
+             appl_qedpower(j) = orders(qed_pos)
+           enddo
+         endif
 c Initialize hiostograms for fixed order runs
          if (fixed_order) call initplot
          firsttime=.false.
@@ -1565,6 +1585,7 @@
       qcdpower(icontr)=QCD_power
       cpower(icontr)=wgtcpower
       orderstag(icontr)=orders_tag
+      amppos(icontr)=amp_pos
       ipr(icontr)=0
       call set_pdg(icontr,nFKSprocess)
 
@@ -2264,17 +2285,18 @@
       return
       end
 
-      subroutine fill_applgrid_weights(vegas_wgt)
-c Fills the ApplGrid weights of appl_common.inc. This subroutine assumes
+      subroutine fill_pineappl_weights(vegas_wgt)
+c Fills the FineAPPL weights of pineappl_common.inc. This subroutine assumes
 c that there is an unique PS configuration: at most one Born, one real
 c and one set of counter terms. Among other things, this means that one
 c must do MC over FKS directories.
       use weight_lines
       implicit none
       include 'nexternal.inc'
-      include 'appl_common.inc'
+      include 'pineappl_common.inc'
       include 'nFKSconfigs.inc'
       include 'genps.inc'
+      integer orders(nsplitorders)
       integer i,j 
       double precision final_state_rescaling,vegas_wgt
       integer              flavour_map(fks_configs)
@@ -2284,24 +2306,18 @@
       common/cproc_combination/iproc_save,eto,etoi,maxproc_found
       integer lo_qcd_to_amp_pos, nlo_qcd_to_amp_pos
       integer pos
-      if (icontr.gt.8) then
-         write (*,*) 'ERROR: too many applgrid weights. '/
-     &        /'Should have at most one of each itype.',icontr
-         stop 1
-      endif
-      do j=1,amp_split_size
-        do i=1,4
+      do i=1,4
+        do j=1,amp_split_size
           appl_w0(i,j)=0d0
           appl_wR(i,j)=0d0
           appl_wF(i,j)=0d0
           appl_wB(i,j)=0d0
-          appl_x1(i,j)=0d0
-          appl_x2(i,j)=0d0
-          appl_QES2(i,j)=0d0
-          appl_muR2(i,j)=0d0
-          appl_muF2(i,j)=0d0
-          appl_qcdpower(i,j) = -1
         enddo
+        appl_x1(i)=0d0
+        appl_x2(i)=0d0
+        appl_QES2(i)=0d0
+        appl_muR2(i)=0d0
+        appl_muF2(i)=0d0
         appl_flavmap(i)=0
       enddo
       appl_event_weight = 0d0
@@ -2316,62 +2332,63 @@
          else
             pos = nlo_qcd_to_amp_pos(qcdpower(i))
          endif
+         ! consistency check
+         if (appl_qcdpower(pos).ne.qcdpower(i)) then
+           write(*,*) 'ERROR in fill_pineappl_weights, QCDpower',
+     %        appl_qcdpower(pos), qcdpower(i)  
+           stop 1
+         endif
 
          if (itype(i).eq.1) then
 c     real
             appl_w0(1,pos)=appl_w0(1,pos)+wgt(1,i)*final_state_rescaling
-            appl_qcdpower(1,pos)=qcdpower(i)
-            appl_x1(1,pos)=bjx(1,i)
-            appl_x2(1,pos)=bjx(2,i)
+            appl_x1(1)=bjx(1,i)
+            appl_x2(1)=bjx(2,i)
             appl_flavmap(1) = flavour_map(nFKS(i))
-            appl_QES2(1,pos)=scales2(1,i)
-            appl_muR2(1,pos)=scales2(2,i)
-            appl_muF2(1,pos)=scales2(3,i)
+            appl_QES2(1)=scales2(1,i)
+            appl_muR2(1)=scales2(2,i)
+            appl_muF2(1)=scales2(3,i)
          elseif (itype(i).eq.2) then
 c     born
             appl_wB(2,pos)=appl_wB(2,pos)+wgt(1,i)*final_state_rescaling
-            appl_qcdpower(2,pos)=qcdpower(i)
-            appl_x1(2,pos)=bjx(1,i)
-            appl_x2(2,pos)=bjx(2,i)
+            appl_x1(2)=bjx(1,i)
+            appl_x2(2)=bjx(2,i)
             appl_flavmap(2) = flavour_map(nFKS(i))
-            appl_QES2(2,pos)=scales2(1,i)
-            appl_muR2(2,pos)=scales2(2,i)
-            appl_muF2(2,pos)=scales2(3,i)
+            appl_QES2(2)=scales2(1,i)
+            appl_muR2(2)=scales2(2,i)
+            appl_muF2(2)=scales2(3,i)
          elseif (itype(i).eq.3 .or. itype(i).eq.4 .or. itype(i).eq.14
      &           .or. itype(i).eq.15)then
 c     virtual, soft-virtual or soft-counter
             appl_w0(2,pos)=appl_w0(2,pos)+wgt(1,i)*final_state_rescaling
             appl_wR(2,pos)=appl_wR(2,pos)+wgt(2,i)*final_state_rescaling
             appl_wF(2,pos)=appl_wF(2,pos)+wgt(3,i)*final_state_rescaling
-            appl_qcdpower(2,pos)=qcdpower(i)
-            appl_x1(2,pos)=bjx(1,i)
-            appl_x2(2,pos)=bjx(2,i)
+            appl_x1(2)=bjx(1,i)
+            appl_x2(2)=bjx(2,i)
             appl_flavmap(2) = flavour_map(nFKS(i))
-            appl_QES2(2,pos)=scales2(1,i)
-            appl_muR2(2,pos)=scales2(2,i)
-            appl_muF2(2,pos)=scales2(3,i)
+            appl_QES2(2)=scales2(1,i)
+            appl_muR2(2)=scales2(2,i)
+            appl_muF2(2)=scales2(3,i)
          elseif (itype(i).eq.5) then
 c     collinear counter            
             appl_w0(3,pos)=appl_w0(3,pos)+wgt(1,i)*final_state_rescaling
             appl_wF(3,pos)=appl_wF(3,pos)+wgt(3,i)*final_state_rescaling
-            appl_qcdpower(3,pos)=qcdpower(i)
-            appl_x1(3,pos)=bjx(1,i)
-            appl_x2(3,pos)=bjx(2,i)
+            appl_x1(3)=bjx(1,i)
+            appl_x2(3)=bjx(2,i)
             appl_flavmap(3) = flavour_map(nFKS(i))
-            appl_QES2(3,pos)=scales2(1,i)
-            appl_muR2(3,pos)=scales2(2,i)
-            appl_muF2(3,pos)=scales2(3,i)
+            appl_QES2(3)=scales2(1,i)
+            appl_muR2(3)=scales2(2,i)
+            appl_muF2(3)=scales2(3,i)
          elseif (itype(i).eq.6) then
 c     soft-collinear counter            
             appl_w0(4,pos)=appl_w0(4,pos)+wgt(1,i)*final_state_rescaling
             appl_wF(4,pos)=appl_wF(4,pos)+wgt(3,i)*final_state_rescaling
-            appl_qcdpower(4,pos)=qcdpower(i)
-            appl_x1(4,pos)=bjx(1,i)
-            appl_x2(4,pos)=bjx(2,i)
+            appl_x1(4)=bjx(1,i)
+            appl_x2(4)=bjx(2,i)
             appl_flavmap(4) = flavour_map(nFKS(i))
-            appl_QES2(4,pos)=scales2(1,i)
-            appl_muR2(4,pos)=scales2(2,i)
-            appl_muF2(4,pos)=scales2(3,i)
+            appl_QES2(4)=scales2(1,i)
+            appl_muR2(4)=scales2(2,i)
+            appl_muF2(4)=scales2(3,i)
          endif
       enddo
       return
@@ -2516,6 +2533,8 @@
       ! stuff for plotting the different splitorders
       integer orders_tag_plot
       common /corderstagplot/ orders_tag_plot
+      integer amp_pos_plot
+      common /campposplot/ amp_pos_plot
       save max_weight
       call cpu_time(tBefore)
       if (icontr.eq.0) return
@@ -2526,12 +2545,27 @@
          do j=1,iwgt
             plot_wgts(j,i)=0d0
          enddo
+         ! The following if lines have been changed with respect to the
+         ! usual (with just 3 plot ids: 20, 11 and 12):
+         !  The kinematics of soft and collinear counterterms may
+         !  be different, for those processes without soft singularities
+         !  from initial(final)-state configurations when the 
+         !  final(initial) confs are integrated (e.g. a a > e+ e-)
+         !  This gives no problem for normal histogramming (and in
+         !  fact plot_id 11 13 and 14 are merged into ibody=2 in
+         !  outfun, but it gives troubles e.g. with applgrid/pineappl. 
+         !  Note that the separation between soft and soft-virtual
+         !  may not be needed in reality
          if (itype(i).eq.2) then
             plot_id(i)=20 ! Born
          elseif(itype(i).eq.1) then
             plot_id(i)=11 ! real-emission
+         elseif(itype(i).eq.5) then
+            plot_id(i)=13 ! collinear counter term
+         elseif(itype(i).eq.6) then
+            plot_id(i)=14 ! soft collinear counter term
          else
-            plot_id(i)=12 ! soft-virtual and counter terms
+            plot_id(i)=12 ! soft-virtual and soft counter term
          endif
 c Loop over all previous icontr. If the plot_id, PDGs and momenta are
 c equal to a previous icountr, add the current weight to the plot_wgts
@@ -2541,12 +2575,12 @@
          do ii=1,i
             if (orderstag(ii).ne.orderstag(i)) cycle
             if (plot_id(ii).ne.plot_id(i)) cycle
-            if (plot_id(i).eq.20 .or. plot_id(i).eq.12) then
+            if (plot_id(i).ne.11) then
                if (.not.pdg_equal(pdg_uborn(1,ii),pdg_uborn(1,i))) cycle
             else
                if (.not.pdg_equal(pdg(1,ii),pdg(1,i))) cycle
             endif
-            if (plot_id(i).eq.20 .or. plot_id(i).eq.12) then
+            if (plot_id(i).ne.11) then
                if (.not.momenta_equal_uborn(momenta(0,1,ii),momenta(0,1
      $              ,i),fks_j_d(nFKS(ii)),fks_i_d(nFKS(ii))
      $                 ,fks_j_d(nFKS(i)) ,fks_i_d(nFKS(i)))) cycle
@@ -2576,6 +2610,7 @@
             enddo
 c call the analysis/histogramming routines
             orders_tag_plot=orderstag(i)
+            amp_pos_plot=amppos(i)
             call outfun(momenta(0,1,i),y_bst(i),www,pdg(1,i),plot_id(i))
          endif
       enddo
 
=== modified file 'Template/NLO/SubProcesses/handling_lhe_events.f'
--- Template/NLO/SubProcesses/handling_lhe_events.f	2020-07-28 10:05:22 +0000
+++ Template/NLO/SubProcesses/handling_lhe_events.f	2021-01-08 10:23:50 +0000
@@ -113,7 +113,7 @@
       include './cuts.inc'
       integer lhaid
       character*20 pdlabel
-      integer iappl
+      logical pineappl
       character*7 event_norm
       character*13 temp
 c     other parameter
 
=== modified file 'Template/NLO/SubProcesses/iproc_map.f'
--- Template/NLO/SubProcesses/iproc_map.f	2017-07-31 08:48:20 +0000
+++ Template/NLO/SubProcesses/iproc_map.f	2021-01-08 10:23:50 +0000
@@ -205,7 +205,7 @@
 
 ************************************************************************
 *     The following routine sets up the flavour map that needs to be
-*     feeded to APPLgrid in the initialization stage.
+*     feeded to PineAPPL in the initialization stage.
 ************************************************************************
       subroutine setup_flavourmap
 *
@@ -215,8 +215,8 @@
       include 'nexternal.inc'
       include 'genps.inc'
       include "leshouche_decl.inc"
-      include "reweight_appl.inc"
-      include "appl_common.inc"
+      include "reweight_pineappl.inc"
+      include "pineappl_common.inc"
 *
       character*200 buffer
       integer procnum,i,l,j,ll,found_a,found_m
@@ -251,8 +251,18 @@
 *     Read the file using a buffer
       do
          read (71,'(a)',err=100,end=100) buffer ! Jump to line 100 when all lines read
+         write(*,*) buffer
          read (buffer,*) kpdflumi,nproc(kpdflumi),
      1        ((pdgs(i,j,kpdflumi),i=1,2),j=1,nproc(kpdflumi))
+         ! check that the allocated arrays are big enough
+         if (kpdflumi.gt.mxpdflumi) then
+            write(*,*) 'ERROR in iproc_map.f, too many processes'
+            write(*,*) 'increase mxpdflumi and max_nproc' //
+     %                 ' inside appl_comon.inc'
+            write(*,*) 'and __max_nproc__ in pineappl_interface.cc'
+            write(*,*) 'Make sure to assign all variables the same value!'
+            stop 1
+         endif
          appl_nproc(kpdflumi) = nproc(kpdflumi)
       enddo
  100  continue
 
=== modified file 'Template/NLO/SubProcesses/madfks_plot.f'
--- Template/NLO/SubProcesses/madfks_plot.f	2018-12-21 15:27:32 +0000
+++ Template/NLO/SubProcesses/madfks_plot.f	2021-01-08 10:23:50 +0000
@@ -11,9 +11,9 @@
       double precision xsecScale_acc(maxscales,maxscales,maxdynscales)
      $     ,xsecPDFr_acc(0:maxPDFs,maxPDFsets)
       common /scale_pdf_print/xsecScale_acc,xsecPDFr_acc
-      integer iappl
-      common /for_applgrid/ iappl
-      include "appl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
+      include "pineappl_common.inc"
       nwgt=1
       if (.not.allocated(weights_info)) allocate(weights_info(1))
       weights_info(nwgt)="central value               "
@@ -109,7 +109,7 @@
 c start with central member of the first set
          call InitPDFm(1,0)
       endif
-      if(iappl.ne.0)then
+      if(pineappl)then
 c Initialize grid parameters to negative values.
          appl_Q2min   = -1d0
          appl_Q2max   = -1d0
@@ -159,14 +159,14 @@
       double precision xsecScale_acc(maxscales,maxscales,maxdynscales)
      $     ,xsecPDFr_acc(0:maxPDFs,maxPDFsets)
       common /scale_pdf_print/xsecScale_acc,xsecPDFr_acc
-      integer iappl
-      common /for_applgrid/ iappl
-      include "appl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
+      include "pineappl_common.inc"
 c
       xnorm=1.d0/float(ncalls0)
       if(useitmax)xnorm=xnorm/float(itmax)
-c Normalization factor for the APPLgrid grids
-      if(iappl.ne.0) appl_norm_histo = 1d0 / dble(ncalls0*itmax)
+c Normalization factor for the PineAPPL grids
+      if(pineappl) appl_norm_histo = 1d0 / dble(ncalls0*itmax)
       call analysis_end(xnorm)
 c Write the accumulated results to a file
       open (unit=34,file='scale_pdf_dependence.dat',status='unknown')
@@ -242,13 +242,17 @@
       double precision xsecScale_acc(maxscales,maxscales,maxdynscales)
      $     ,xsecPDFr_acc(0:maxPDFs,maxPDFsets)
       common /scale_pdf_print/xsecScale_acc,xsecPDFr_acc
-      integer iappl
-      common /for_applgrid/ iappl
-      include "appl_common.inc"
-c Born, n-body or (n+1)-body contribution:
+      integer amp_pos_plot
+      common /campposplot/ amp_pos_plot
+      logical pineappl
+      common /for_pineappl/ pineappl
+      include "pineappl_common.inc"
+c Born, n-body or (n+1)-body contribution. For more information
+C  about itype, see the comments about plot_id inside
+C  the fill_plot subroutine of fks_singular.f
       if(itype.eq.11) then
          ibody=1 ! (n+1)-body
-      elseif(itype.eq.12)then
+      elseif(itype.ge.12 .and. itype.le.14)then
          ibody=2 ! n-body
       elseif(itype.eq.20)then
          ibody=3 ! Born
@@ -275,8 +279,21 @@
          enddo
          p(4,i)=pmass(i)
       enddo
-      if(iappl.ne.0)then
-         appl_itype     = ibody
+      if(pineappl)then
+         appl_itype = ibody
+         if(ibody.eq.2)then
+c      special treatment for collinear
+c      and soft collinear counterterms
+c      (see comments inside the fill_plots subroutine)
+           if(itype.eq.13)then
+            !collinear counterterm
+             appl_itype = 4
+           elseif(itype.eq.14)then
+            !solf-collinear counterterm
+             appl_itype = 5
+           endif
+         endif
+         appl_amp_pos   = amp_pos_plot
          appl_www_histo = www(1)
       endif
       call analysis_fill(p,istatus,ipdg,www,ibody)
 
=== modified file 'Template/NLO/SubProcesses/makefile_fks_dir'
--- Template/NLO/SubProcesses/makefile_fks_dir	2020-07-28 10:05:22 +0000
+++ Template/NLO/SubProcesses/makefile_fks_dir	2021-01-08 10:23:50 +0000
@@ -25,17 +25,19 @@
 $(LIBDIR)libpdf.a $(LIBDIR)libcernlib.a
 
 # Files for all executables
-FILES= $(patsubst %.f,%.o,$(wildcard parton_lum_*.f)) $(patsubst	\
-	%.f,%.o,$(wildcard matrix_*.f)) FKSParams.o real_me_chooser.o	\
-	chooser_functions.o genps_fks.o setcuts.o setscales.o		\
-	veto_xsec.o $(patsubst %.f,%.o,$(wildcard b_sf_???.f))		\
-	born.o sborn_sf.o extra_cnt_wrapper.o $(patsubst		\
-	%.f,%.o,$(wildcard born_cnt_*.f)) fks_Sij.o			\
-	$(fastjetfortran_madfks) fks_singular.o montecarlocounter.o	\
-	reweight_xsec.o boostwdir2.o initcluster.o cluster.o		\
-	splitorders_stuff.o reweight.o get_color.o 			\
-	iproc_map.o MC_integer.o $(reweight_xsec_events_pdf_dummy)	\
-	$(applgrid_interface) weight_lines.o mint_module.o polfit.o
+FILES= $(patsubst %.f,%.o,$(wildcard parton_lum_*.f)) 			\
+	 $(patsubst %.f,%.o,$(wildcard matrix_*.f)) 			\
+	 FKSParams.o real_me_chooser.o chooser_functions.o recmom.o	\
+	 genps_fks.o setcuts.o setscales.o veto_xsec.o		 	\
+	 $(patsubst %.f,%.o,$(wildcard b_sf_???.f)) 			\
+	 born.o sborn_sf.o extra_cnt_wrapper.o				\
+	 $(patsubst %.f,%.o,$(wildcard born_cnt_*.f))			\
+	 fks_Sij.o $(fastjetfortran_madfks) fks_singular.o		\
+	 montecarlocounter.o reweight_xsec.o boostwdir2.o       	\
+	 initcluster.o cluster.o splitorders_stuff.o 			\
+	 reweight.o get_color.o iproc_map.o				\
+	 MC_integer.o $(reweight_xsec_events_pdf_dummy)			\
+	 $(pineappl_interface) weight_lines.o mint_module.o polfit.o
 
 # Files needed for mintFO & mintMC
 RUN= $(FO_ANALYSE) $(FILES) cuts.o pythia_unlops.o recluster.o 		\
@@ -58,7 +60,7 @@
 SYM= $(RUN) BinothLHADummy.o symmetry_fks_v3.o write_ajob.o
 
 # Files for event reweighting to get e.g. scale and PDF uncertainties
-RWGFILES=reweight_xsec_events.o handling_lhe_events.o			\
+RWGFILES=reweight_xsec_events.o handling_lhe_events.o recmom.o 		\
 	$(patsubst %.f,%.o,$(wildcard parton_lum_*.f)) setscales.o	\
 	cuts.o $(fastjetfortran_madfks) boostwdir2.o veto_xsec.o	\
 	$(reweight_xsec_events_pdf_dummy) fill_MC_mshell.o cluster.o	\
@@ -86,17 +88,17 @@
 	$(CXX) $(CXXFLAGS) -c -I. -I$(ANADIR) $(FO_INCLUDEPATHS) $^
 
 gensym: $(SYM)
-	$(FC) $(LDFLAGS) -o gensym $(SYM) $(APPLLIBS) $(LINKLIBS) $(FJLIBS)
+	$(FC) $(LDFLAGS) -o gensym $(SYM) $(APPLLIBS) $(LINKLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS) $(FJLIBS)
 	rm fks_singular.o
 	strip gensym
 
 test_soft_col_limits: $(TEST)
-	$(FC) $(LDFLAGS) -o test_soft_col_limits $(TEST) $(APPLLIBS) $(LINKLIBS) $(FJLIBS)
+	$(FC) $(LDFLAGS) -o test_soft_col_limits $(TEST) $(APPLLIBS) $(LINKLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS) $(FJLIBS)
 	rm test_soft_col_limits.o
 	strip test_soft_col_limits
 
 check_poles: $(POLES) $(libmadloop)
-	$(FC) -o check_poles $(POLES) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(LDFLAGS)
+	$(FC) -o check_poles $(POLES) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS) $(LDFLAGS)
 	strip check_poles
 
 madevent_mintMC: $(MINTMC) $(libmadloop) makefile $(LIBS)
 
=== modified file 'Template/NLO/SubProcesses/open_output_files.f'
--- Template/NLO/SubProcesses/open_output_files.f	2015-03-05 16:12:37 +0000
+++ Template/NLO/SubProcesses/open_output_files.f	2021-01-08 10:23:50 +0000
@@ -32,6 +32,18 @@
       subroutine HwU_write_file
       implicit none
       double precision xnorm
+c     PineAPPL commons (this may not be the best place to put it)
+      include "reweight_pineappl.inc"
+      include "pineappl_common.inc"
+      logical pineappl
+      common /for_pineappl/ pineappl
+      integer j
+      if(pineappl)then
+         do j=1,nh_obs
+           appl_obs_num = j
+           call APPL_term
+         enddo
+      endif
       open (unit=99,file='MADatNLO.HwU',status='unknown')
       xnorm=1d0
       call HwU_output(99,xnorm)
 
=== renamed file 'Template/NLO/SubProcesses/appl_common.inc' => 'Template/NLO/SubProcesses/pineappl_common.inc'
--- Template/NLO/SubProcesses/appl_common.inc	2016-02-10 14:18:55 +0000
+++ Template/NLO/SubProcesses/pineappl_common.inc	2021-01-08 10:23:50 +0000
@@ -1,37 +1,36 @@
 *     -*-fortran-*-
 ************************************************************************
 *     NLO mode of aMC@NLO
+*     Common blocks needed for the PineAPPL interface
 ************************************************************************
       include 'orders.inc'
 *     Amp split size
       integer appl_amp_split_size 
-      common / appl_common_fixed / appl_amp_split_size
+      integer appl_qcdpower(amp_split_size), appl_qedpower(amp_split_size)
+      common / appl_common_fixed / appl_amp_split_size, appl_qcdpower, appl_qedpower
 
 *     Initial state map to be used together with the info in appl_flavmap
       integer appl_nlumi
-      integer mxpdflumi
-      integer max_nproc
-      parameter(mxpdflumi=121,max_nproc=121)
+      include 'pineappl_maxproc.inc'
       integer appl_nproc(mxpdflumi)
       integer appl_lumimap(2,max_nproc,mxpdflumi)
       common / appl_common_lumi / appl_lumimap,appl_nproc,appl_nlumi
 
 *     Event weights, kinematics, observables etc. that are different event by event
-      double precision appl_x1(4, amp_split_size), 
-     1                 appl_x2(4, amp_split_size)
-      double precision appl_muF2(4, amp_split_size),
-     1                 appl_muR2(4, amp_split_size),
-     2                 appl_QES2(4, amp_split_size)
+      double precision appl_x1(4), 
+     1                 appl_x2(4)
+      double precision appl_muF2(4),
+     1                 appl_muR2(4),
+     2                 appl_QES2(4)
       double precision appl_w0(4, amp_split_size),
      1                 appl_wR(4, amp_split_size),
      2                 appl_wF(4, amp_split_size),
      3                 appl_wB(4, amp_split_size)
-      integer appl_qcdpower(4, amp_split_size)
       integer appl_flavmap(4)
       common / appl_common_weights / appl_x1,appl_x2,
      1                               appl_muF2,appl_muR2,appl_QES2,
      2                               appl_w0,appl_wR,appl_wF,appl_wB,
-     3                               appl_flavmap,appl_qcdpower
+     3                               appl_flavmap
 
 *     Parameters of the grids.
 *     These parameters can optionally be singularly specified by the user,
@@ -52,6 +51,7 @@
       double precision appl_obs_bins(0:100)
       integer appl_obs_nbins
       integer appl_itype
+      integer appl_amp_pos
       integer appl_obs_num
       common / appl_common_histokin / appl_www_histo,appl_norm_histo,
      1                                appl_obs_histo,
@@ -59,10 +59,11 @@
      3                                appl_obs_bins,
      4                                appl_obs_nbins,
      5                                appl_itype,
+     5                                appl_amp_pos,
      6                                appl_obs_num
 
 *     This is the actual hadronic event weight, computed with the LHAPDF set
-*     defined above, needed for checks in the C++ applgrid interface
+*     defined above, needed for checks in the C++ pineappl interface
 *     Save also the vegas weight, needed to fill histograms and for checks
 *     appl_www_histo is the same weight used to will the amcatnlo
 *     histograms in the madfks routine
 
=== renamed file 'Template/NLO/SubProcesses/appl_interface.cc' => 'Template/NLO/SubProcesses/pineappl_interface.cc'
--- Template/NLO/SubProcesses/appl_interface.cc	2014-07-01 08:21:03 +0000
+++ Template/NLO/SubProcesses/pineappl_interface.cc	2021-01-08 10:23:50 +0000
@@ -1,45 +1,376 @@
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstddef>
+#include <cstdio>
+#include <cstdlib>
+#include <iomanip>
 #include <iostream>
-#include <stdlib.h>
+#include <iterator>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include <pineappl_capi.h>
+
+#include "orders.h"
+
+#include "pineappl_maxproc.h"
 
 /*
   fNLO mode of aMCatNLO
 */
 
-void (*appl_initptr)() = 0;
-void (*appl_fillptr)() = 0;
-void (*appl_fillrefptr)() = 0;
-void (*appl_fillrefoutptr)() = 0;
-void (*appl_recoptr)() = 0;
-void (*appl_termptr)() = 0;
-
-extern "C" void appl_init_() {
-  // std::cerr<<"I am in appl_init (C++ version)"<<std::endl;
-  if (appl_initptr) appl_initptr(); 
-}
-
-extern "C" void appl_fill_() {
-  //std::cerr<<"I am in appl_fill (C++ version)"<<std::endl;
-  if (appl_fillptr) appl_fillptr(); 
-}
-
-extern "C" void appl_fill_ref_() {
-  //std::cerr<<"I am in appl_fill_ref (C++ version)"<<std::endl;
-  if (appl_fillrefptr) appl_fillrefptr(); 
-}
-
-extern "C" void appl_fill_ref_out_() {
-  // std::cerr<<"I am in appl_fill_ref_out (C++ version)"<<std::endl;
-  if (appl_fillrefoutptr) appl_fillrefoutptr(); 
-}
-
-extern "C" void appl_reco_() {
-  // std::cerr<<"I am in appl_reco (C++ version)"<<std::endl;
-  if (appl_recoptr) appl_recoptr(); 
-}
-
-extern "C" void appl_term_() {
-  //std::cerr<<"I am in appl_term (C++ version)"<<std::endl;
-  if (appl_termptr) appl_termptr();  
-}
-
-
+// Declare grids
+std::vector<pineappl_grid*> grid_obs;
+
+// translates an index from the range [0, __amp_split_size) to the indices need by `fill_grid`
+std::vector<std::vector<int>> translation_tables;
+
+// Information defined at the generation (configuration) step, that does
+// not vary event by event
+extern "C" struct
+{
+    int amp_split_size; // Maximum number of coupling-combinations
+    int qcdpower[__amp_split_size]; // Power of alpha_s for each amp_split
+    int qedpower[__amp_split_size]; // Power of alpha for each amp_split
+} appl_common_fixed_;
+
+// Map of the PDF combinations from aMC@NLO - structure for each
+// "subprocess" i, has some number nproc[i] pairs of parton
+// combinations. To be used together with the info in appl_flavmap.
+extern "C" struct
+{
+    int lumimap[__max_nproc__][__max_nproc__][2]; // (paired) subprocesses per combination
+    int nproc[__max_nproc__]; // number of separate (pairwise) subprocesses for this combination
+    int nlumi; // overall number of combinations ( 0 < nlumi <= __mxpdflumi__ )
+} appl_common_lumi_;
+
+// Event weights, kinematics, etc. that are different event by event
+extern "C" struct
+{
+    double x1[4], x2[4];
+    double muF2[4], muR2[4], muQES2[4];
+    double W0[__amp_split_size][4], WR[__amp_split_size][4];
+    double WF[__amp_split_size][4], WB[__amp_split_size][4];
+    int flavmap[4];
+} appl_common_weights_;
+
+// Parameters of the grids.
+// These parameters can optionally be singularly specified by the user,
+// but if no specification is given, the code will use the default values.
+extern "C" struct
+{
+    double Q2min, Q2max;
+    double xmin, xmax;
+    int nQ2, Q2order;
+    int nx, xorder;
+} appl_common_grid_;
+
+// Parameters of the histograms
+extern "C" struct
+{
+    double www_histo, norm_histo;
+    double obs_histo;
+    double obs_min, obs_max;
+    double obs_bins[101];
+    int obs_nbins;
+    int itype_histo;
+    int amp_pos;
+    int obs_num;
+} appl_common_histokin_;
+
+// Event weight and cross section
+extern "C" struct
+{
+    double event_weight, vegaswgt;
+    double xsec12, xsec11, xsec20;
+} appl_common_reco_;
+
+// Check if a file exists
+bool file_exists(const std::string& s)
+{
+    if (std::FILE* testfile = std::fopen(s.c_str(), "r"))
+    {
+        std::fclose(testfile);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+// Banner
+std::string Banner()
+{
+    return "\n"
+        "    █████╗ ███╗   ███╗ ██████╗██████╗ ██╗      █████╗ ███████╗████████╗\n"
+        "   ██╔══██╗████╗ ████║██╔════╝██╔══██╗██║     ██╔══██╗██╔════╝╚══██╔══╝\n"
+        "   ███████║██╔████╔██║██║     ██████╔╝██║     ███████║███████╗   ██║\n"
+        "   ██╔══██║██║╚██╔╝██║██║     ██╔══██╗██║     ██╔══██║╚════██║   ██║\n"
+        "   ██║  ██║██║ ╚═╝ ██║╚██████╗██████╔╝███████╗██║  ██║███████║   ██║\n"
+        "   ╚═╝  ╚═╝╚═╝     ╚═╝ ╚═════╝╚═════╝ ╚══════╝╚═╝  ╚═╝╚══════╝   ╚═╝\n";
+}
+
+extern "C" void appl_init_()
+{
+    // Grid Initialization and definition of the observables.
+    // Construct the input file name according to its position in the
+    // vector "grid_obs".
+    std::size_t const index = grid_obs.size();
+
+    int lo_power = 9999;
+    int nlo_power = 0;
+
+    for (int i = 0; i != appl_common_fixed_.amp_split_size; ++i)
+    {
+        int sum = appl_common_fixed_.qcdpower[i] + appl_common_fixed_.qedpower[i];
+
+        lo_power = std::min(lo_power, sum);
+        nlo_power = std::max(nlo_power, sum);
+    }
+
+    // TODO: are there any situations is which there are NLOs but no LOs?
+
+    // we assume that there is always at least one LO, and zero or one NLO
+    assert((nlo_power == (lo_power + 2)) || (nlo_power == lo_power));
+
+    std::vector<uint32_t> subgrid_params;
+
+    translation_tables.emplace_back();
+    translation_tables.back().reserve(appl_common_fixed_.amp_split_size);
+
+    for (int i = 0; i != appl_common_fixed_.amp_split_size; ++i)
+    {
+        int const qcd = appl_common_fixed_.qcdpower[i];
+        int const qed = appl_common_fixed_.qedpower[i];
+        int const sum = qcd + qed;
+
+        translation_tables.back().push_back(subgrid_params.size() / 4);
+
+        if (sum == lo_power)
+        {
+            // WB
+            subgrid_params.insert(subgrid_params.end(), { qcd / 2, qed / 2, 0, 0 });
+        }
+        else if (sum == nlo_power)
+        {
+            // W0
+            subgrid_params.insert(subgrid_params.end(), { qcd / 2, qed / 2, 0, 0 });
+            // WR
+            subgrid_params.insert(subgrid_params.end(), { qcd / 2, qed / 2, 1, 0 });
+            // WF
+            subgrid_params.insert(subgrid_params.end(), { qcd / 2, qed / 2, 0, 1 });
+        }
+    }
+
+    // Define the settings for the interpolation in x and Q2.
+    // These are common to all the grids computed.
+    // If values larger than zero (i.e. set by the user) are found the default
+    // settings are replaced with the new ones.
+    auto* key_vals = pineappl_keyval_new();
+
+    if (appl_common_grid_.nQ2 > 0)
+    {
+        pineappl_keyval_set_int(key_vals, "nq2", appl_common_grid_.nQ2);
+    }
+    // Max and min value of Q2
+    if (appl_common_grid_.Q2min > 0.0)
+    {
+        pineappl_keyval_set_double(key_vals, "q2min", appl_common_grid_.Q2min);
+    }
+    if (appl_common_grid_.Q2max > 0.0)
+    {
+        pineappl_keyval_set_double(key_vals, "q2max", appl_common_grid_.Q2max);
+    }
+    // Order of the polynomial interpolation in Q2
+    if (appl_common_grid_.Q2order > 0)
+    {
+        pineappl_keyval_set_int(key_vals, "q2order", appl_common_grid_.Q2order);
+    }
+    // Number of points for the x interpolation
+    if (appl_common_grid_.nx > 0)
+    {
+        pineappl_keyval_set_int(key_vals, "nx", appl_common_grid_.nx);
+    }
+    // Min and max value of x
+    if (appl_common_grid_.xmin > 0.0)
+    {
+        pineappl_keyval_set_double(key_vals, "xmin", appl_common_grid_.xmin);
+    }
+    if (appl_common_grid_.xmax > 0.0)
+    {
+        pineappl_keyval_set_double(key_vals, "xmax", appl_common_grid_.xmax);
+    }
+    // Order of the polynomial interpolation in x
+    if (appl_common_grid_.xorder > 0)
+    {
+        pineappl_keyval_set_int(key_vals, "xorder", appl_common_grid_.xorder);
+    }
+
+    // Set up the PDF luminosities
+    auto* lumi = pineappl_lumi_new();
+
+    // Loop over parton luminosities
+    for (int ilumi = 0; ilumi < appl_common_lumi_.nlumi; ilumi++)
+    {
+        int nproc = appl_common_lumi_.nproc[ilumi];
+
+        std::vector<int32_t> pdg_ids;
+        pdg_ids.reserve(2 * nproc);
+
+        for (int iproc = 0; iproc != nproc; ++iproc)
+        {
+            pdg_ids.push_back(appl_common_lumi_.lumimap[ilumi][iproc][0]);
+            pdg_ids.push_back(appl_common_lumi_.lumimap[ilumi][iproc][1]);
+        }
+
+        pineappl_lumi_add(lumi, nproc, pdg_ids.data(), nullptr);
+    }
+
+    // Use the reweighting function
+    pineappl_keyval_set_bool(key_vals, "reweight", true);
+
+    // valid choices are: "LagrangeSubgrid", "NtupleSubgrid"
+    pineappl_keyval_set_string(key_vals, "subgrid_type", "LagrangeSubgrid");
+
+    // Create a grid with the binning given in the "obsbins[Nbins+1]" array
+    grid_obs.push_back(pineappl_grid_new(
+        lumi,
+        subgrid_params.size() / 4,
+        subgrid_params.data(),
+        appl_common_histokin_.obs_nbins,
+        appl_common_histokin_.obs_bins,
+        key_vals
+    ));
+
+    pineappl_keyval_delete(key_vals);
+    pineappl_lumi_delete(lumi);
+}
+
+extern "C" void appl_fill_()
+{
+    // Check that itype ranges from 1 to 5.
+    int itype = appl_common_histokin_.itype_histo;
+    if ((itype < 1) || (itype > 5))
+    {
+        std::cout << "amcblast ERROR: Invalid value of itype = " << itype << std::endl;
+        std::exit(-10);
+    }
+
+    // this is the second index of the WB/R/F/0 arrays
+    int index = appl_common_histokin_.amp_pos - 1;
+
+    // aMC@NLO weights. Four grids, ordered as {W0,WR,WF,WB}.
+    double(&W0)[4] = appl_common_weights_.W0[index];
+    double(&WR)[4] = appl_common_weights_.WR[index];
+    double(&WF)[4] = appl_common_weights_.WF[index];
+    double(&WB)[4] = appl_common_weights_.WB[index];
+
+    int ilumi;
+    int nlumi = appl_common_lumi_.nlumi;
+    double ttol = 1e-100;
+    double x1, x2;
+    double scale2;
+    double obs = appl_common_histokin_.obs_histo;
+
+    // Histogram number
+    int nh = appl_common_histokin_.obs_num - 1;
+
+    // translate (index,nh) -> index of the PineAPPL grid
+    int const grid_index = translation_tables.at(nh).at(index);
+
+    int k;
+
+    switch (itype)
+    {
+    case 1:
+        k = 0;
+        break;
+
+    case 2:
+    case 3:
+        k = 1;
+        break;
+
+    case 4:
+        k = 2;
+        break;
+
+    case 5:
+        k = 3;
+        break;
+
+    default:
+        assert( false );
+    }
+
+    x1 = appl_common_weights_.x1[k];
+    x2 = appl_common_weights_.x2[k];
+
+    static std::vector<std::vector<std::vector<double>>> x1Saved(5, std::vector<std::vector<double>>(
+        grid_obs.size(), std::vector<double>(pineappl_grid_order_count(grid_obs[nh]), 0.0)));
+    static std::vector<std::vector<std::vector<double>>> x2Saved(5, std::vector<std::vector<double>>(
+        grid_obs.size(), std::vector<double>(pineappl_grid_order_count(grid_obs[nh]), 0.0)));
+
+    if (x1 == x1Saved[itype - 1][nh][grid_index] && x2 == x2Saved[itype - 1][nh][grid_index])
+    {
+        return;
+    }
+    else
+    {
+        x1Saved[itype - 1][nh][grid_index] = x1;
+        x2Saved[itype - 1][nh][grid_index] = x2;
+    }
+
+    scale2 = appl_common_weights_.muF2[k];
+    ilumi = appl_common_weights_.flavmap[k] - 1;
+
+    if (x1 < 0.0 || x1 > 1.0 || x2 < 0.0 || x2 > 1.0)
+    {
+        std::cout << "amcblast ERROR: Invalid value of x1 and/or x2 = " << x1 << " " << x2
+                  << std::endl;
+        std::exit(-10);
+    }
+
+    if (x1 == 0.0 && x2 == 0.0)
+    {
+        return;
+    }
+
+    if (std::fabs(W0[k]) > ttol)
+    {
+        pineappl_grid_fill(grid_obs[nh], x1, x2, scale2, grid_index + 0, obs, ilumi, W0[k]);
+    }
+
+    if (std::fabs(WR[k]) > ttol)
+    {
+        pineappl_grid_fill(grid_obs[nh], x1, x2, scale2, grid_index + 1, obs, ilumi, WR[k]);
+    }
+
+    if (std::fabs(WF[k]) > ttol)
+    {
+        pineappl_grid_fill(grid_obs[nh], x1, x2, scale2, grid_index + 2, obs, ilumi, WF[k]);
+    }
+
+    if (std::fabs(WB[k]) > ttol)
+    {
+        pineappl_grid_fill(grid_obs[nh], x1, x2, scale2, grid_index, obs, ilumi, WB[k]);
+    }
+}
+
+extern "C" void appl_term_()
+{
+    // Histogram number
+    int const nh = appl_common_histokin_.obs_num - 1;
+
+    // convert between gs^2 and alphas, and normalize the grid by (hbarc)^2 and the number of runs
+    pineappl_grid_scale_by_order(grid_obs[nh], 4.0 * std::acos(-1.0), 1.0, 1.0, 1.0,
+        389379660.0 * appl_common_histokin_.norm_histo);
+
+    // Write grid to file
+    pineappl_grid_write(grid_obs[nh], ("grid_obs_" + std::to_string(nh) + "_out.pineappl").c_str());
+
+    pineappl_grid_delete(grid_obs[nh]);
+}
 
=== renamed file 'Template/NLO/SubProcesses/appl_interface_dummy.f' => 'Template/NLO/SubProcesses/pineappl_interface_dummy.f'
--- Template/NLO/SubProcesses/appl_interface_dummy.f	2014-06-26 08:45:41 +0000
+++ Template/NLO/SubProcesses/pineappl_interface_dummy.f	2021-01-08 10:23:50 +0000
@@ -5,20 +5,6 @@
       return
       end
 c
-      subroutine APPL_fill_ref_out
-      implicit none
-      write(6,*) "You should not be here in APPL_fill dummy!"
-      stop
-      return
-      end
-c
-      subroutine APPL_fill_ref
-      implicit none
-      write(6,*) "You should not be here in APPL_fill dummy!"
-      stop
-      return
-      end
-c
       subroutine APPL_init
       implicit none
       write(6,*) "You should not be here in APPL_init dummy!"
 
=== added file 'Template/NLO/SubProcesses/recmom.f'
--- Template/NLO/SubProcesses/recmom.f	1970-01-01 00:00:00 +0000
+++ Template/NLO/SubProcesses/recmom.f	2021-01-08 10:23:50 +0000
@@ -0,0 +1,139 @@
+      subroutine recombine_momenta(R, etaph, reco_l, reco_q, p_in, pdg_in, p_out, pdg_out)
+      implicit none
+      ! recombine photons with the closest fermion if the distance is
+      ! less than R and if the rapidity of photons is < etaph (etaph < 0
+      ! means no cut). Output a new set of momenta and pdgs corresponding
+      ! to the recombined particles. If recombination occurs the photon
+      ! disappears from the output particles
+      ! arguments
+      include 'nexternal.inc'
+      double precision R, etaph, p_in(0:4,nexternal), p_out(0:4,nexternal)
+      logical reco_l, reco_q
+      integer pdg_in(nexternal), pdg_out(nexternal)
+      ! local variables
+      integer nq, nl
+      integer id_ph
+      parameter (id_ph=22)
+      integer n_ph, i_ph
+      integer i,j
+      integer ifreco
+      double precision dreco, dthis
+      integer skip
+      logical is_light_charged_fermion
+      double precision R2, eta
+      ! 
+      integer times_reco
+      common/to_times_reco/ times_reco
+      ! reset everything
+      do j=1,nexternal
+        pdg_out(j)=0
+        do i=0,4
+          p_out(i,j)=0d0
+        enddo
+      enddo
+