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

[Olivier Mattelaer]

First non beta release

Target time for the release: Friday.
Please perform all the checks that you want on this version.

[marco zaro (marco-zaro)]

Hi Olivier,
wait one second, we still need to merge the 4FS branch into 3.0.4…
I will update everybody once this is done
Cheers,

marco

> On 24 Mar 2021, at 11:13, Olivier Mattelaer <email address hidden> wrote:
>
> Olivier Mattelaer has proposed merging lp:~maddevelopers/mg5amcnlo/3.0.4 into lp:mg5amcnlo/3.0.
>
> Commit message:
> First non beta release
>
> Requested reviews:
> MadTeam (madteam)
>
> For more details, see:
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.4/+merge/400095
>
> Target time for the release: Friday.
> Please perform all the checks that you want on this version.
> --
> The attached diff has been truncated due to its size.
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.0.4.
> === modified file 'MadSpin/decay.py'
> --- MadSpin/decay.py 2020-06-21 21:31:33 +0000
> +++ MadSpin/decay.py 2021-03-24 10:12:56 +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"]...

[lp:~maddevelopers/mg5amcnlo/3.0.4]

991. By olivier-mattelaer on 2021-03-24

fixing one IOtest

992. By marco zaro on 2021-03-24

real-born links are written for processees with only QCD splittings

993. By marco zaro on 2021-03-24

merged 4FS branch

994. By olivier-mattelaer on 2021-03-25

merge with 2.9.3 -> forcing DP for CF computationn

995. By olivier-mattelaer on 2021-03-25

fixed nlo_multicore_lowmem

996. By olivier-mattelaer on 2021-03-25

merge 3.0.4-qcdonly

997. By olivier-mattelaer on 2021-03-25

update VERSION/UPDATE

998. By olivier-mattelaer on 2021-03-25

remove print statement

999. By olivier-mattelaer on 2021-03-25

includ fxfx patch of Rik

[Olivier Mattelaer (olivier-mattelaer)]

approve will be release as 3.1.0

review: Approve

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2020-06-21 21:31:33 +0000
+++ MadSpin/decay.py	2021-03-25 20:41:08 +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()):
@@ -3355,8 +3371,8 @@
                     if key[0]=='full':
                         path=key[1]
                         end_signal="5 0 0 0 0\n"  # before closing, write down the seed 
-                        external.stdin.write(end_signal)
-                        ranmar_state=external.stdout.readline()
+                        external.stdin.write(end_signal.encode())
+                        ranmar_state=external.stdout.readline().decode()
                         ranmar_file=pjoin(path,'ranmar_state.dat')
                         ranmar=open(ranmar_file, 'w')
                         ranmar.write(ranmar_state)
@@ -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)]
@@ -3997,12 +4037,21 @@
                 pass
             else:
                 stdin_text="5 0 0 0 0"
-                external.stdin.write(stdin_text)
-                external.stdin.close()
-                external.stdout.close()
+                try:
+                    external.stdin.write(stdin_text)
+                except Exception:
+                    pass
+                try:
+                    external.stdin.close()
+                except Exception:
+                    pass
+                try:
+                    external.stdout.close()
+                except Exception:
+                    pass
                 external.terminate()       
                 del external
-
+            
         self.calculator = {}
 
 
@@ -4058,7 +4107,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 +4136,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 +4221,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-03-25 20:41:08 +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    
@@ -446,14 +451,6 @@
             except ValueError:
                 raise self.InvalidCmd('second argument should be a real number.')
         
-        elif args[0] == 'BW_effect':
-            if args[1] in [0, False,'.false.', 'F', 'f', 'False', 'no']:
-                args[1] = 0
-            elif args[1] in [1, True,'.true.', 'T', 't', 'True', 'yes']:
-                args[1] = 1
-            else:
-                raise self.InvalidCmd('second argument should be either T or F.')
-        
         elif args[0] == 'curr_dir':
             if not os.path.isdir(args[1]):
                 raise self.InvalidCmd('second argument should be a path to a existing directory')
@@ -491,11 +488,9 @@
 
         # 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':
-                return self.list_completion(text, ['True', 'False']) 
             if args[1] == 'ms_dir':
                 return self.path_completion(text, '.', only_dirs = True)
         elif args[1] == 'ms_dir':
@@ -508,20 +503,18 @@
     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 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 +540,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 +561,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 +639,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 +733,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 +825,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 +866,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 +1006,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 +1075,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 +1086,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 +1216,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 +1263,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 +1300,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 +1388,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 +1566,7 @@
             # file if needed.
             while 1:
                 try:
-                    decay = decay_file.next()
+                    decay = next(decay_file)
                     break
                 except StopIteration:
                     eff = self.efficiency
@@ -1580,7 +1608,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 +1702,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 +1750,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 'MadSpin/src/lha_read_ms.f'
--- MadSpin/src/lha_read_ms.f	2012-12-12 17:38:07 +0000
+++ MadSpin/src/lha_read_ms.f	2021-03-25 20:41:08 +0000
@@ -261,7 +261,7 @@
 
 
 
-      subroutine LHA_get_real(npara,param,value,name,var,def_value_num)
+      subroutine LHA_get_real_silent(npara,param,value,name,var,def_value_num)
 c----------------------------------------------------------------------------------
 c     finds the parameter named "name" in param and associate to "value" in value
 c----------------------------------------------------------------------------------
@@ -284,12 +284,17 @@
 c
 c     local
 c
-      logical found
+      logical found, log
       integer i
+
+      log = .false.
+      goto 10
+      entry LHA_get_real(npara,param,value,name,var,def_value_num)
 c
 c     start
 c
-      i=1
+      
+ 10   i=1
       found=.false.
       do while(.not.found.and.i.le.npara)
          ctemp=param(i)
@@ -305,8 +310,10 @@
          i=i+1
       enddo
       if (.not.found) then
-         write (*,*) "Warning: parameter ",name," not found"
-         write (*,*) "         setting it to default value ",def_value_num
+         if (log) then
+            write (*,*) "Warning: parameter ",name," not found"
+            write (*,*) "         setting it to default value ",def_value_num
+         endif
          var=def_value_num
       endif
 

=== 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-03-25 20:41:08 +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-03-25 20:41:08 +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
@@ -58,16 +58,19 @@
   %(gridpack)s     = gridpack  !True = setting up the grid pack
   %(time_of_flight)s = time_of_flight ! threshold (in mm) below which the invariant livetime is not written (-1 means not written)
   %(event_norm)s =  event_norm       ! average/sum. Normalization of the weight in the LHEF
-$MLM
-$CKKW
+$mlm
+$ckkw
 #*********************************************************************
 #
 #*********************************************************************
-# handling of the helicities:
-#  0: sum over all helicities
-#  1: importance sampling over helicities
+# Phase-Space Optimization strategy (basic options)
 #*********************************************************************
    %(nhel)s  = nhel          ! using helicities importance sampling or not.
+                             ! 0: sum over helicity, 1: importance sampling
+   %(sde_strategy)s  = sde_strategy  ! default integration strategy (hep-ph/2021.xxxxx)
+                             ! 1 is old strategy (using amp square)
+			     ! 2 is new strategy (using only the denominator)
+$psoptim			     
 #*********************************************************************
 # Generation bias, check the wiki page below for more information:   *
 #  'cp3.irmp.ucl.ac.be/projects/madgraph/wiki/LOEventGenerationBias' *

=== modified file 'Template/LO/Source/.make_opts'
--- Template/LO/Source/.make_opts	2019-01-23 16:03:38 +0000
+++ Template/LO/Source/.make_opts	2021-03-25 20:41:08 +0000
@@ -13,10 +13,11 @@
 
 # Rest of the makefile
 ifeq ($(origin FFLAGS),undefined)
-FFLAGS= -O -w -fbounds-check -fPIC
+FFLAGS= -w -fPIC
 #FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
 endif
 
+FFLAGS += $(GLOBAL_FLAG)
 
 # REMOVE MACFLAG IF NOT ON MAC OR FOR F2PY
 UNAME := $(shell uname -s)

=== 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-03-25 20:41:08 +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-03-25 20:41:08 +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,9 +297,14 @@
 c  a "diffractive" photon
 c      
       q2max=xmu*xmu
-      if(ih .eq. 3) then  !from the electron
-          fx(7)=epa_electron(x,q2max)
-      elseif(ih .eq. 2) then  !from a proton without breaking
+      if(abs(ih) .eq. 3.or.abs(ih) .eq. 4) then !from the electron
+         write(*,*) 'impossible call (or was it) to pdf-> please reporrt'
+         stop 23
+         fx(7)=epa_lepton(x,q2max, ih)
+      elseif(ih .eq. 2) then    !from a proton without breaking
+         write(*,*) 'impossible call (or was it) to pdf-> please reporrt'
+         stop 23
+         ! isssue with next call since this function takes now a third argument (beamid not define here)
           fx(7)=epa_proton(x,q2max)
       endif      
       

=== removed file 'Template/LO/Source/PDF/pdf_lhapdf.f'
=== 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-03-25 20:41:08 +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-03-25 20:41:08 +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,10 +139,11 @@
             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)
+          pdg2pdf = epa_proton(x,xmu*xmu, beamid)
+
       else
          write (*,*) 'beam type not supported in lhadpf'
          stop 1

=== modified file 'Template/LO/Source/dsample.f'
--- Template/LO/Source/dsample.f	2020-05-13 09:45:34 +0000
+++ Template/LO/Source/dsample.f	2021-03-25 20:41:08 +0000
@@ -257,6 +257,7 @@
 c     Need to start from scratch. This is clunky but I'll just
 c     remove the grid, so we are clean
 c
+      goto 200
       write(*,*) "Trying w/ fresh grid"
       open(unit=25,file='ftn25',status='unknown',err=102)
       write(25,*) ' '
@@ -329,7 +330,7 @@
 c
 c     All done
 c
-      open(unit=66,file='results.dat',status='unknown')
+200   open(unit=66,file='results.dat',status='unknown')
       i=1
       do while(xmean(i) .ne. 0 .and. i .lt. cur_it)
          i=i+1
@@ -453,6 +454,10 @@
       CUMULATED_TIMING = t_after - CUMULATED_TIMING
 
       if (N_EVALS.eq.0) then
+         write(outUnit,*) '<lo_statistics> '
+         write(outUnit,*) '<cumulated_time>'//trim(toStr_real(CUMULATED_TIMING))
+     &        //'</cumulated_time>'
+         write(outUnit,*) '</lo_statistics>'
         return
       endif
       
@@ -2270,7 +2275,11 @@
       integer                                      nsteps
       character*40          result_file,where_file
       common /sample_status/result_file,where_file,nsteps
-
+c
+c     
+c
+       logical init_mode
+       common/to_determine_zero_hel/init_mode
 c----
 c  Begin Code
 c----
@@ -2278,6 +2287,9 @@
          write(*,*) nb_pass_cuts, 
      &    ' points passed the cut but all returned zero'
          write(*,*) 'therefore considering this contribution as zero'
+         if (init_mode) then
+            call print_zero_amp()
+         endif
       else if (nb_pass_cuts.gt.0.and.nb_pass_cuts.lt.1000)then
          write(*,*) 'only', nb_pass_cuts, 
      &    ' points passed the cut and they all returned zero'

=== modified file 'Template/LO/Source/invarients.f'
--- Template/LO/Source/invarients.f	2011-04-16 05:25:14 +0000
+++ Template/LO/Source/invarients.f	2021-03-25 20:41:08 +0000
@@ -185,10 +185,10 @@
       end
 
 
-      subroutine map_invarients(Minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
+      subroutine map_invarients(Minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming,nb_tchannel)
 c****************************************************************************
 c     Determines mappings for each structure of invarients onto integration
-c     variables.  Input: Ninvar, iforest.  Output: Minvar, ninvar
+c     variables.  Input: Ninvar, iforest.  Output: Minvar, ninvar, nb_tchannel
 c****************************************************************************
       implicit none
 c
@@ -201,11 +201,13 @@
 c
       integer Minvar(maxdim,lmaxconfigs),nconfigs,ninvar,nexternal,nincoming
       integer mincfig,maxcfig
+      integer nb_tchannel
 c
 c     Local
 c
       integer iconfig, jgrid,j, nbranch
       logical found,tchannel
+      integer ns_channel
 c
 c     Global
 c
@@ -224,20 +226,30 @@
 c     
 c     Try simple mapping if nconfigs = 1
 c
-
+      nb_tchannel = 0
       if (nconfigs .eq. 1) then
 c         do j=1,3*nbranch-4+2
          do j=1,maxdim
             minvar(j,mincfig)=j
          enddo
          jgrid=j-1
+
+         write(*,*) 'Determine nb_t'
+         ns_channel=1
+         do while((iforest(1,-ns_channel,mincfig) .ne. 1.and.iforest(1,-ns_channel,mincfig) .ne. 2).and.ns_channel.lt.nbranch)
+            ns_channel=ns_channel+1
+         enddo
+         ns_channel=ns_channel - 1
+         nb_tchannel=nbranch-ns_channel-1
+         write(*,*) 'T-channel found: ',nb_tchannel 
+         
       else
 c      if (ep) jgrid=1
 c      if (pp) jgrid=2
       do iconfig=mincfig,maxcfig
          tchannel = .false.         
          do j=1,nbranch-1
-            if (iforest(1,-j,iconfig) .eq. 1) then
+            if (iforest(1,-j,iconfig) .eq. 1.or.(nincoming.eq.2.and.iforest(1,-j,iconfig) .eq. 2)) then
                tchannel=.true.
             endif
             jgrid=jgrid+1

=== 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-03-25 20:41:08 +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/run.inc'
--- Template/LO/Source/run.inc	2019-03-23 22:04:23 +0000
+++ Template/LO/Source/run.inc	2021-03-25 20:41:08 +0000
@@ -97,4 +97,9 @@
      &etamax4pdg, mxxmin4pdg,mxxpart_antipart
 
        double precision small_width_treatment
-       common/narrow_width/small_width_treatment
\ No newline at end of file
+       common/narrow_width/small_width_treatment
+
+       double precision tmin_for_channel
+       integer sde_strat ! 1 means standard single diagram enhancement strategy,
+c                          2 means approximation by the denominator of the propagator
+       common/TO_CHANNEL_STRAT/tmin_for_channel, sde_strat
\ No newline at end of file

=== modified file 'Template/LO/Source/run_printout.f'
--- Template/LO/Source/run_printout.f	2012-03-22 06:16:35 +0000
+++ Template/LO/Source/run_printout.f	2021-03-25 20:41:08 +0000
@@ -22,12 +22,17 @@
       write(6,*)  'Collider parameters:'
       write(6,*)  '--------------------'
 
+      ab(1) = '?'
+      ab(2) = '?'
       do i=1,2
          IF(LPP(i).EQ. 0) ab(i)='e'
          IF(LPP(i).EQ. 1) ab(i)='P'
          IF(LPP(i).EQ.-1) ab(i)='Pb'
-         IF(LPP(i).EQ.2) ab(i)='a2'
-         IF(LPP(i).EQ.3) ab(i)='a3'
+         IF(LPP(i).EQ.2) ab(i)='a'
+         IF(LPP(i).EQ.3) ab(i)='e-'
+         IF(LPP(i).EQ.-3) ab(i)='e+'
+         IF(LPP(i).EQ.4) ab(i)='m-'
+         IF(LPP(i).EQ.-4) ab(i)='m+'
       enddo
 
       ene=2d0*dsqrt(ebeam(1)*ebeam(2))

=== 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-03-25 20:41:08 +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/cluster.f'
--- Template/LO/SubProcesses/cluster.f	2013-03-21 12:45:09 +0000
+++ Template/LO/SubProcesses/cluster.f	2021-03-25 20:41:08 +0000
@@ -210,7 +210,8 @@
       integer ignum, ipnum, ipids(nexternal,4,2:nexternal)
 C $B$ IFOREST $B$ !this is a tag for MadWeight
       integer i, iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest, tstrategy
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
       common/to_sprop/sprop,tprid
@@ -404,7 +405,8 @@
       common/to_mconfigs/mapconfig, this_config
 C $B$ IFOREST $B$ !this is a tag for MadWeight
       integer i, iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest,	tstrategy
 C $E$ IFOREST $E$ !this is a tag for MadWeight
 C $B$ DECAYBW $E$ !this is a tag for MadWeight
 

=== removed file 'Template/LO/SubProcesses/dipole.inc'
=== removed file 'Template/LO/SubProcesses/dipolesub.f'
=== removed file 'Template/LO/SubProcesses/epsterms.f'
=== removed file 'Template/LO/SubProcesses/finiteterms.f'
=== 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-03-25 20:41:08 +0000
@@ -163,7 +163,8 @@
 
 
       integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/ iforest, tstrategy
 
       integer            mapconfig(0:lmaxconfigs), this_config
       common/to_mconfigs/mapconfig, this_config
@@ -321,7 +322,7 @@
       endif
       pswgt = 1d0
       jac   = 1d0
-      call one_tree(iforest(1,-max_branch,iconfig),mincfig,
+      call one_tree(iforest(1,-max_branch,iconfig), tstrategy(iconfig),mincfig,
      &     nbranch,P,M,S,X,jac,pswgt)
 c
 c     Add what I think are the essentials
@@ -520,7 +521,7 @@
       integer nbranch,ndim,nconfigs
       integer ninvar
       integer nparticles,nfinal
-
+      integer nb_tchannel
 
 c
 c     Arguments
@@ -572,6 +573,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
@@ -593,7 +596,7 @@
          nconfigs = 1
          mincfig=iconfig
          maxcfig=iconfig
-         call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
+         call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming,nb_tchannel)
          maxwgt=0d0
          nparticles   = nexternal
          nfinal       = nparticles-nincoming
@@ -630,7 +633,7 @@
       return
       end
 
-      subroutine one_tree(itree,iconfig,nbranch,P,M,S,X,jac,pswgt)
+      subroutine one_tree(itree,tstrategy,iconfig,nbranch,P,M,S,X,jac,pswgt)
 c************************************************************************
 c     Calculates the momentum for everything below in the tree until
 c     it reaches the end.
@@ -650,8 +653,10 @@
 c     Arguments
 c
       integer itree(2,-max_branch:-1) !Structure of configuration
+      integer tstrategy ! current strategy for t-channel
       integer iconfig                 !Which configuration working on
       double precision P(0:3,-max_branch:max_particles)
+      double precision pother(0:3), ptemp(0:3), pboost(0:3), ptemp2(0:3)
       double precision M(-max_branch:max_particles)
       double precision S(-max_branch:0)
 c      double precision spole(-max_branch:0),swidth(-max_branch:0)
@@ -662,12 +667,13 @@
 c     Local
 c
       logical pass
+      double precision tmass(-max_branch:-1)
       integer ibranch,i,ns_channel,nt_channel,ix  !,nerr
-c      data nerr/0/
+      integer iopposite ! index for t-channel mapping for the part not handle by itree
+c     data nerr/0/
       double precision smin,smax,totmass,totmassin,xa2,xb2,wgt
       double precision costh,phi,tmin,tmax,t
       double precision ma2,mb2,m12,mn2,s1, mi2
-      double precision tmass(-max_branch:-1)
 c
 c     External
 c
@@ -701,7 +707,12 @@
 c     the s channel p1+p2
 c
       ns_channel=1
-      do while(itree(1,-ns_channel) .ne. 1 .and.ns_channel.lt.nbranch)
+      iopposite = 1
+      if (abs(tstrategy).eq.1) then
+         iopposite = 2
+      endif
+      
+      do while(itree(1,-ns_channel) .ne. iopposite .and.ns_channel.lt.nbranch)
          m(-ns_channel)=0d0                 
          ns_channel=ns_channel+1         
       enddo
@@ -711,7 +722,6 @@
       if (nt_channel .eq. 0 .and. nincoming .eq. 2) then
          ns_channel=ns_channel-1
       endif
-
 c
 c     Determine masses for all intermediate states.  Starting
 c     from outer most (real particle) states
@@ -809,7 +819,7 @@
          totmass=totmass-m(itree(2,ibranch))      !for remaining particles
          smin = totmass**2                        !This affects t_min/max
          smax = (m(ibranch) - m(itree(2,ibranch)))**2
-
+         
          if (smin .gt. smax) then
             jac=-3d0
             return
@@ -835,18 +845,239 @@
 c
 c     From here, on we can just pretend this is a 2->2 scattering with
 c     Pa                    + Pb     -> P1          + P2
+
+      if (tstrategy.eq.-2.or.tstrategy.eq.-1) then
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
+cc       T-channel ping-pong strategy starting with 2
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+         
+c     No -flipping case:      
+c      p(0,itree(ibranch,1)) + p(0,2) -> p(0,ibranch)+ p(0,itree(ibranch,2))
+c       -  M(ibranch) is the total mass available (Pa+Pb)^2
+c       - M(ibranch-1) is the mass of P2  (all the remaining particles)
+c
+c     With flipping case
+c      p(0,itree(ibranch,1)) + pother -> p(0,ibranch)+ p(0,itree(ibranch,2))
+c       - pother = p(0,itree(ibranch,1)) -p(0,itree(ibranch,2))         
+c       - M(ibranch) is the total mass available (Pa+Pb)^2
+c       - M(ibranch-1) is the mass of P2  (all the remaining particles)      
+c
+c     This assumes that P(0, ibranch) is set to the T-channel propa (likely)
+c      do ibranch = -ns_channel-1,-nbranch,-1
+c         totmass=totmass+m(itree(2,ibranch))
+c      enddo
+      do ibranch=-ns_channel-1,-nbranch+1,-1
+c         totmass=totmass-m(itree(2,ibranch))      !for remaining particles
+c         smin = totmass**2                        !This affects t_min/max
+c         smax = (m(ibranch) - m(itree(2,ibranch)))**2
+
+         if (ibranch.ne.-ns_channel-1)then
+            pother(:) = P(:,ibranch+1)
+            iopposite = ibranch +1
+         else
+            pother(:) = p(:,2)
+            iopposite = abs(tstrategy)
+         endif
+         s1  = m(ibranch)**2                        !Total mass available
+         ma2 = dot(pother, pother)
+         mb2 = dot(P(0,itree(1,ibranch)),P(0,itree(1,ibranch)))
+         m12 = m(itree(2,ibranch))**2
+         mn2 = m(ibranch-1)**2
+c
+c$$$         write(*,*) itree(1, ibranch), '-----------T----------', itree(2, ibranch), 'm=', m(itree(2, ibranch))
+c$$$         write(*,*)        '                       |           '
+c$$$         write(*,*)   	   '                       |           '
+c$$$         write(*,*)   	   '                       | ',ibranch
+c$$$         write(*,*)   	   '                       |           '
+c$$$         write(*,*)        '                       |           '
+c$$$         if (ibranch.ne.-ns_channel-1)then
+c$$$            write(*,*) iopposite, '-----------T---------- m=', m(ibranch-1)
+c$$$         else
+c$$$            write(*,*) 2, '-----------T---------- m=', m(ibranch-1)
+c$$$            write(*,*) m(1), m(2), m(3), m(4), m(5)
+c$$$         endif
+c$$$         write(*,*) 'Pa', P(0,itree(1, ibranch)),P(1,itree(1, ibranch)),P(2,itree(1, ibranch)),P(3,itree(1, ibranch))
+c$$$         if (ibranch.ne.-ns_channel-1) then
+c$$$            write(*,*) 'Pb', P(0,iopposite),P(1,iopposite),P(2,iopposite), P(3,iopposite)
+c$$$            do i=0,3
+c$$$               pother(i) = P(i,itree(1, ibranch)) + P(i,ibranch+1) 
+c$$$            enddo
+c$$$         else
+c$$$            write(*,*) 'Pb', P(0,2),P(1,2),P(2,2),P(3,2)
+c$$$            do i=0,3
+c$$$               pother(i) = P(i,1) + P(i,2)
+c$$$            enddo
+c$$$         endif
+c$$$         do i=0,3
+c$$$            pother(i) = P(i,itree(1, ibranch)) + P(i,iopposite)
+c$$$         enddo
+c$$$         write(*,*) 'DSQRT(s1) = ', m(ibranch), DSQRT(dot(pother, pother))
+c$$$c         if (m(ibranch)**2.ne.dot(pother, pother)) stop 1
+c$$$         write(*,*) 'm12= Pd**2 = ', m12 ,DSQRT(m12)
+c$$$         write(*,*) 'mn2 = Pc**2 =', mn2, DSQRT(mn2)
+         
+C     WRITE(*,*) 'Enertering yminmax',sqrt(s1),sqrt(m12),sqrt(mn2)
+         
+         call yminmax(s1,0d0,m12,ma2,mb2,mn2,tmin,tmax)
+c         call yminmax(s1,0d0,m12,ma2,mb2,smax,tmin_temp,tmax_temp)
+c         if (tmin_temp.lt.tmin) tmin = tmin_temp
+c         if (tmax_temp.gt.tmax) tmax = tmax_temp
+         
+c
+c     Call for 0<x<1
+c
+c         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+c     &        .5d0*(tmin/stot+1d0),
+c     &        .5d0*(tmax/stot+1d0))
+c         t   = Stot*(x(-ibranch)*2d0-1d0)
+c
+c     call for -1<x<1
+c
+
+c         write(*,*) 'tmin, tmax/ temp',tmin,tmax, tmin_temp, tmax_temp
+
+c         if (nt_channel.ge.2)then
+c            tmin = max(tmin, -stot)
+c         endif
+c      if ((tmax-tmin)/stot.gt.0.1)then
+c            call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+c     $           0d0, 1d0)
+c         t = stot*(-x(-ibranch))
+      
+c      else if (tmax/stot.gt.-0.01.and.tmin/stot.lt.-0.02)then
+c         set tmax to 0. The idea is to be sure to be able to hit zero
+c         and not to be block by numerical inacuracy
+c         tmax = max(tmax,0d0) !This line if want really t freedom
+c         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+c     $        0d0, -tmin/stot)
+c         t = stot*(-x(-ibranch))
+
+c      else
+         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+     $        -tmax/stot, -tmin/stot)
+         t = stot*(-x(-ibranch))
+c      endif
+
+c      call yminmax(s1,0d0,m12,ma2,mb2,mn2,tmin_temp,tmax_temp)
+      if (t .lt. tmin .or. t .gt. tmax) then
+         jac=-3d0
+         return
+      endif
+c
+c     tmin and tmax set to -s,+s for jacobian because part of jacobian
+c     was determined from choosing the point x from the grid based on
+c     tmin and tmax.  (ie wgt contains some of the jacobian)
+c
+         tmin=-stot
+         tmax= stot
+         call sample_get_x(wgt,x(nbranch+(-ibranch-1)*2),
+     &        nbranch+(-ibranch-1)*2,iconfig,0d0,1d0)
+         phi = 2d0*pi*x(nbranch+(-ibranch-1)*2)
+         jac = jac*(tmax-tmin)*2d0*pi /2d0 ! I need /2d0 if -1<x<1
+
+         
+c
+c     Finally generate the momentum. The call is of the form
+c     pa+pb -> p1+ p2; t=(pa-p1)**2;   pr = pa-p1
+c     gentcms(pa,pb,t,phi,m1,m2,p1,pr) 
+c
+         if (itree(1,ibranch).gt.-ns_channel-1)then
+            mi2 = m(itree(1,ibranch))**2
+         else
+            mi2 = tmass(itree(1,ibranch))
+         endif
+         tmass(ibranch) = t
+         call gentcms(p(0,itree(1,ibranch)),p(0,iopposite),t,phi,mi2,
+     &        m(itree(2,ibranch)),m(ibranch-1),p(0,itree(2,ibranch)),
+     &        p(0,ibranch),jac)
+c$$$         write(*,*) 'RESULT'
+c$$$         write(*,*) 'pa', p(0,itree(1,ibranch)),p(1,itree(1,ibranch)),p(2,itree(1,ibranch)),p(3,itree(1,ibranch))
+c$$$         write(*,*) 'pb', p(0,iopposite),p(1,iopposite),p(2,iopposite),p(3,iopposite)
+c$$$         write(*,*) '->'
+c$$$         write(*,*) 'pc', p(0,itree(2,ibranch)),p(1,itree(2,ibranch)),p(2,itree(2,ibranch)),p(3,itree(2,ibranch))
+c$$$         do i =0,3
+c$$$            pother(i) = p(i,itree(1,ibranch)) + p(i,iopposite) - p(i,itree(2,ibranch))
+c$$$         enddo
+c$$$         write(*,*) 'pc', p(0,itree(2,ibranch)),p(1,itree(2,ibranch)),p(2,itree(2,ibranch)),p(3,itree(2,ibranch))
+c$$$         write(*,*) 'pd', pother(0), pother(1), pother(2), pother(3) , DSQRT(dot(pother,pother))
+c$$$         write(*,*) 'T channel'
+c$$$         write(*,*) 'pr', p(0,ibranch),p(1,ibranch),p(2,ibranch),p(3,ibranch)
+c$$$         write(*,*) 'pa-pc', p(0,itree(1,ibranch))-p(0,itree(2,ibranch)),p(1,itree(1,ibranch))-p(1,itree(2,ibranch))
+
+
+         
+         if (jac .lt. 0d0) then
+c            nerr=nerr+1
+c            if(nerr.le.5)
+c     $           write(*,*) 'Failed gentcms',iconfig,ibranch
+            return              !Failed, probably due to negative x
+         endif
+
+         pswgt = pswgt/(4d0*dsqrt(lambda(s1,ma2,mb2)))
+      enddo
+
+c     
+c     We need to get the momentum of the last external particle.
+c     This should just be the sum of p(0,2) and the remaining
+c     momentum from our last t channel 2->2
+c
+      if (nt_channel.eq.1) then
+c$$$         write(*,*) 'need to assign last', itree(2,-nbranch)
+c$$$         write(*,*) 'nbranch is at', nbranch
+c$$$         do i=-nbranch,nexternal
+c$$$            write(*,*) 'p',i, p(0,i),p(1,i),p(2,i),p(3,i)
+c$$$         enddo
+         do i=0,3
+            p(i,itree(2,-nbranch)) = p(i,-nbranch+1)+p(i,2)
+         enddo
+      else
+c$$$                  write(*,*) 'need to assign last', itree(2,-nbranch)
+c$$$         write(*,*) 'nbranch is at', nbranch
+c$$$         do i=-nbranch,nexternal
+c$$$            write(*,*) 'p',i, p(0,i),p(1,i),p(2,i),p(3,i)
+c$$$         enddo
+         do i=0,3
+            p(i,itree(2,-nbranch)) = p(i,-nbranch+1)+p(i,-nbranch+2)
+         enddo
+      endif
+
+
+      else if (tstrategy.eq.2.or.tstrategy.eq.1) then
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+cc       T-channel One side eat all strategy ending with 2
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c     Now perform the t-channel decay sequence. Most of this comes from: 
+c     Particle Kinematics Chapter 6 section 3 page 166
+c
+c     From here, on we can just pretend this is a 2->2 scattering with
+c     Pa                    + Pb     -> P1          + P2
 c     p(0,itree(ibranch,1)) + p(0,2) -> p(0,ibranch)+ p(0,itree(ibranch,2))
 c     M(ibranch) is the total mass available (Pa+Pb)^2
 c     M(ibranch-1) is the mass of P2  (all the remaining particles)
 c
       do ibranch=-ns_channel-1,-nbranch+1,-1
          s1  = m(ibranch)**2                        !Total mass available
-         ma2 = m(2)**2
+         ma2 = m(tstrategy)**2
          mb2 = dot(P(0,itree(1,ibranch)),P(0,itree(1,ibranch)))
          m12 = m(itree(2,ibranch))**2
          mn2 = m(ibranch-1)**2
 c         write(*,*) 'Enertering yminmax',sqrt(s1),sqrt(m12),sqrt(mn2)
          call yminmax(s1,0d0,m12,ma2,mb2,mn2,tmin,tmax)
+
+         if(.false.) then
+             write(*,*) itree(1, ibranch), 'a----------T----------', itree(2, ibranch), 'm=', m(itree(2, ibranch))
+             write(*,*)        '                       |           '
+             write(*,*)   	   '                       |           '
+             write(*,*)   	   '                       | ',ibranch
+             write(*,*)   	   '                       |           '
+             write(*,*)        '                       |           '
+             write(*,*) tstrategy, '-----------T---------- m=', m(ibranch-1), ibranch-1
+             write(*,*) 'S', dsqrt(s1), 'm_top=', dsqrt(mb2), 'M_bottom', dsqrt(ma2)
+             
+c     write(*,*) m(1), m(2), m(3), m(4), m(5)
+            write(*,*) 'Pa', P(0,itree(1, ibranch)),P(1,itree(1, ibranch)),P(2,itree(1, ibranch)),P(3,itree(1, ibranch))
+         endif
 c
 c     Call for 0<x<1
 c
@@ -859,20 +1090,92 @@
 c
 
 c         write(*,*) 'tmin, tmax',tmin,tmax
-
-      if (tmax.gt.-0.01.and.tmin.lt.-0.02)then
+         if(.false.) then
+            ! NOT VALIDATED METHOD, momenta are ok but not jacobian
+            
+            call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+     $        0d0, 1d0)
+            costh= 2d0*x(-ibranch)-1d0
+            call sample_get_x(wgt,x(nbranch+(-ibranch-1)*2),
+     &        nbranch+(-ibranch-1)*2,iconfig,0d0,1d0)
+            phi = 2d0*pi*x(nbranch+(-ibranch-1)*2)
+            jac = jac * 4d0*pi
+                     m12 = m(itree(2,ibranch))**2
+         mn2 = m(ibranch-1)**2
+         call mom2cx(dsqrt(s1),m(itree(2,ibranch)),m(ibranch-1),costh,phi,
+     &        p(0,itree(2,ibranch)),pother)
+
+         I= itree(2,ibranch)
+         DO I=0,3
+            pboost(I) = P(I,tstrategy) + P(I,itree(1, ibranch))
+         ENDDO
+         
+         call boostm(p(0,itree(2,ibranch)),pboost,m(itree(2,ibranch)),p(0,itree(2,ibranch)))
+         call boostm(pother,pboost,m(ibranch),pother)
+
+         do I=0,3
+            p(I,ibranch) = pother(i) - p(i, tstrategy)
+         enddo
+         if(.false.)then
+         write(*,*) 'input'
+         write(*,*) 'p(tstrategy=',tstrategy,')', p(0,tstrategy), p(1,tstrategy), p(2,tstrategy), p(3,tstrategy)
+         I=itree(1, ibranch)
+         write(*,*) 'p(',I,',)', p(0,i), p(1,i), p(2,i), p(3,i)
+         write(*,*) 'output'
+         I= itree(2,ibranch)
+         write(*,*) 'p(',I,',)', p(0,i), p(1,i), p(2,i), p(3,i)
+         write(*,*) 'pother', pother(0),pother(1),pother(2),pother(3)
+         write(*,*) 'check'
+         do i=0,3
+            ptemp(i) = p(i,tstrategy) + p(i, itree(1, ibranch))
+         enddo
+         write(*,*) 'pa+pb', ptemp(0), ptemp(1), ptemp(2), ptemp(3), dsqrt(dot(ptemp, ptemp))
+         do i=0,3
+            ptemp(i) = pother(i) + p(i,itree(2,ibranch))
+         enddo
+         write(*,*) 'p1+p2', ptemp(0), ptemp(1), ptemp(2), ptemp(3), dsqrt(dot(ptemp, ptemp))
+         
+         write(*,*) 'Tchannel'
+         I = ibranch
+         write(*,*) 'p(',I,',)', p(0,i), p(1,i), p(2,i), p(3,i), dot(p(0,I), p(0,I))
+         endif
+         
+         pswgt = pswgt/(4d0*dsqrt(lambda(s1,ma2,mb2)))
+
+      else
+
+c     test of impact of low t part
+c         if (nt_channel.ge.2)then
+c            tmin = max(tmin,  -stot)
+c         endif
+
+
+         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+     $        -tmax/stot, -tmin/stot)
+         t = stot*(-x(-ibranch))
+         
+
+c         if ((tmax-tmin)/stot.gt.0.1)then
+c            call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+c     $           0d0, 1d0)
+c         t = stot*(-x(-ibranch))
+c     if (dabs(tmax - tmin)/stot.gt.0.05d0) then
+c         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+c     $        0d0,  1d0)
+c     set tmax to 0 and tmin to -1 The idea is to avoid dimension correlation
+c     the condition ensure a minimum efficiency in the generation of events
+c         t = stot*(-x(-ibranch))
+c      else if (tmax.gt.-0.01.and.tmin.lt.-0.02)then
 c         set tmax to 0. The idea is to be sure to be able to hit zero
 c         and not to be block by numerical inacuracy
-c         tmax = max(tmax,0d0) !This line if want really t freedom
-         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
-     $        0d0, -tmin/stot)
-         t = stot*(-x(-ibranch))
-
-      else
-         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
-     $        -tmax/stot, -tmin/stot)
-         t = stot*(-x(-ibranch))
-      endif
+c         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+c     $        0d0, -tmin/stot)
+c         t = stot*(-x(-ibranch))
+c      else
+c         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+c     $        -tmax/stot, -tmin/stot)
+c         t = stot*(-x(-ibranch))
+c      endif
 
          if (t .lt. tmin .or. t .gt. tmax) then
             jac=-3d0
@@ -901,7 +1204,7 @@
             mi2 = tmass(itree(1,ibranch))
          endif
          tmass(ibranch) = t
-         call gentcms(p(0,itree(1,ibranch)),p(0,2),t,phi, mi2,
+         call gentcms(p(0,itree(1,ibranch)),p(0,tstrategy),t,phi, mi2,
      &        m(itree(2,ibranch)),m(ibranch-1),p(0,itree(2,ibranch)),
      &        p(0,ibranch),jac)
 
@@ -913,6 +1216,7 @@
          endif
 
          pswgt = pswgt/(4d0*dsqrt(lambda(s1,ma2,mb2)))
+      endif
       enddo
 c
 c     We need to get the momentum of the last external particle.
@@ -920,10 +1224,16 @@
 c     momentum from our last t channel 2->2
 c
       do i=0,3
-         p(i,itree(2,-nbranch)) = p(i,-nbranch+1)+p(i,2)
+         p(i,itree(2,-nbranch)) = p(i,-nbranch+1)+p(i,tstrategy)
       enddo
 
 
+      
+      else
+         write(*,*) 'not supported tstrategy'
+         stop 2
+      endif
+
       endif                     !t-channel stuff
 
 
@@ -958,6 +1268,18 @@
          call boostm(p(0,itree(1,i)),p(0,i),m(i),p(0,itree(1,i)))
          call boostm(p(0,itree(2,i)),p(0,i),m(i),p(0,itree(2,i)))
       enddo
+c$$$      write(*,*) '****'
+c$$$      do i=-nbranch,nexternal
+c$$$         write(*,*) 'mass', i, m(i)
+c$$$      enddo
+c$$$      do i=-nbranch,nexternal
+c$$$         write(*,*) 'p',i, p(0,i),p(1,i),p(2,i),p(3,i)
+c$$$      enddo
+c$$$      do i =0,3
+c$$$         pother(i) = p(i,1) + p(i,2) - p(i,3) -p(i,4)
+c$$$      enddo
+c$$$      write(*,*) 'p5 expected', pother(0), pother(1),pother(2),pother(3)
+
       jac = jac*wgt
       if (.not. pass) jac = -99
       end
@@ -1022,7 +1344,7 @@
 c
 c     Arguments
 c
-      double precision t,phi,m1,m2               !inputs
+      double precision t,phi,m1,m2,ma               !inputs
       double precision pa(0:3),pb(0:3),jac
       double precision p1(0:3),pr(0:3)           !outputs
 c
@@ -1046,7 +1368,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
@@ -1313,6 +1634,146 @@
       return
       end
 
-
+      double precision function get_channel_cut(p, config)
+      implicit none
+
+      include 'maxconfigs.inc'
+      include 'nexternal.inc'
+      include 'genps.inc'
+      include 'maxamps.inc'
+      include 'coupl.inc'
+c     include 'run.inc'
+
+      double precision p(0:3, nexternal)
+      integer config
+      
+
+      integer iforest(2,-max_branch:-1,lmaxconfigs)
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/ iforest, tstrategy
+
+      integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
+      integer tprid(-max_branch:-1,lmaxconfigs)
+      common/to_sprop/sprop,tprid
+
+      double precision stot,m1,m2
+      common/to_stot/stot,m1,m2
+
+      double precision tmin_for_channel
+       integer sde_strat ! 1 means standard single diagram enhancement strategy,
+c      	      	      	   2 means approximation by the	denominator of the propagator
+       common/TO_CHANNEL_STRAT/tmin_for_channel,	sde_strat
+      
+      integer            mapconfig(0:lmaxconfigs), this_config
+      common/to_mconfigs/mapconfig, this_config
+
+      double precision      spole(maxinvar),swidth(maxinvar),bwjac
+      common/to_brietwigner/spole          ,swidth          ,bwjac
+
+      double precision ptemp(0:3, -nexternal:nexternal)
+      integer i,j
+      integer d1, d2
+      double precision t
+      double precision dot
+      external dot
+      integer ns_channel
+      integer nb_tchannel
+      integer nbranch
+      double precision tmp, tmp2
+      
+      double precision ZERO
+      parameter (ZERO=0d0)
+      double precision prmass(-nexternal:0,lmaxconfigs)
+      double precision prwidth(-nexternal:0,lmaxconfigs)
+      integer pow(-nexternal:0,lmaxconfigs)
+      logical first_time
+      save prmass,prwidth,pow
+      data first_time /.true./
+
+      double precision Mass, Width
+      
+      include 'configs.inc'
+
+      if(sde_strat.eq.1.and.tmin_for_channel.eq.-1)then
+         get_channel_cut = 1d0
+         return
+      endif
+      
+      if (first_time) then
+         include 'props.inc'
+         first_time=.false.
+      endif
+      
+      do i = 1, nexternal
+         do j =0,3
+            ptemp(j,i) = p(j,i)
+            ptemp(j,-i) = 0d0
+         enddo
+      enddo
+
+      nbranch = nexternal -2
+      ns_channel=1
+      do while((iforest(1,-ns_channel,config) .ne. 1.and.iforest(1,-ns_channel,config) .ne. 2).and.ns_channel.lt.nbranch)
+         ns_channel=ns_channel+1
+      enddo
+      ns_channel=ns_channel - 1
+      nb_tchannel=nbranch-ns_channel-1
+c      write(*,*) 'T-channel found: ',nb_tchannel
+
+
+
+
+
+      
+      get_channel_cut = 1.
+      if (nb_tchannel.lt.2.and.sde_strat.eq.1)then
+         get_channel_cut = 1.
+         return
+      endif
+      
+      do i = 1, nexternal-3
+         d1 = iforest(1, -i, config)
+         d2 = iforest(2, -i, config)
+         do j=0,3
+            if (d1.gt.0.and.d1.le.2) then
+               ptemp(j,-i) = ptemp(j,-i) - ptemp(j, d1)
+            else
+               ptemp(j,-i) = ptemp(j,-i)+ptemp(j, d1)
+            endif
+            if (d2.gt.0.and.d2.le.2) then
+               ptemp(j,-i) = ptemp(j,-i) - ptemp(j, d2)
+            else
+               ptemp(j,-i) = ptemp(j,-i)+ptemp(j, d2)
+            endif
+         enddo
+         if (tprid(-i,config).ne.0)then
+            if(sde_strat.eq.2)then
+               t = dot(ptemp(0,-i), ptemp(0,-i))
+               Mass  = prmass(-i, config)
+               get_channel_cut = get_channel_cut / ((t-Mass)*(t+Mass))**2
+            endif
+c            write(*,*) i, "t, Mass, fact", t, Mass, ((t-Mass)*(t+Mass))**2,get_channel_cut
+            t = t/stot 
+            if (t.lt.tmin_for_channel)then
+                get_channel_cut = get_channel_cut * exp((t-tmin_for_channel)/(t+1))
+c               get_channel_cut = get_channel_cut * (t+1)/(1+tmin_for_channel)
+c            else if(t.gt.2*tmin_for_channel)then
+c               get_channel_cut = get_channel_cut * (2*tmin_for_channel-t)/tmin_for_channel
+            endif
+         else
+            if(sde_strat.eq.2)then
+               t = dot(ptemp(0,-i), ptemp(0,-i))
+               Mass  = prmass(-i, config)
+               Width = prwidth(-i, config)
+               tmp = (t-Mass)*(t+Mass)
+               tmp2 = Mass*Width
+               get_channel_cut = get_channel_cut* (tmp**2 - tmp2**2)/(tmp**2 + tmp2**2)**2 
+            endif
+c            write(*,*) i, "s, Mass, Width, fact", t, Mass, Width, (((t-Mass)*(t+Mass) )**2 + Width**2*Mass**2), get_channel_cut
+         endif
+      enddo
+c      write(*,*) 'final for config', config, get_channel_cut
+      return
+      end
 
 

=== modified file 'Template/LO/SubProcesses/makefile'
--- Template/LO/SubProcesses/makefile	2017-07-13 10:04:07 +0000
+++ Template/LO/SubProcesses/makefile	2021-03-25 20:41:08 +0000
@@ -30,18 +30,27 @@
 
 # Source files
 
+MATRIX_HEL = $(patsubst %.f,%.o,$(wildcard matrix*_orig.f))
+MATRIX = $(patsubst %.f,%.o,$(wildcard matrix*_optim.f))
+ifeq ($(strip $(MATRIX_HEL)),)
+        MATRIX = $(patsubst %.f,%.o,$(wildcard matrix*.f))
+endif
+
+
 PROCESS= driver.o myamp.o genps.o unwgt.o setcuts.o get_color.o \
          cuts.o cluster.o reweight.o initcluster.o addmothers.o setscales.o \
 	 idenparts.o dummy_fct.o \
          $(patsubst %.f,%.o,$(wildcard auto_dsig*.f)) \
-         $(patsubst %.f,%.o,$(wildcard matrix*.f))
 
 SYMMETRY = symmetry.o idenparts.o 
 
 # Binaries
 
-$(PROG): $(PROCESS) auto_dsig.o $(LIBS)
-	$(FC) -o $(PROG) $(PROCESS) $(LINKLIBS) $(LDFLAGS) $(BIASDEPENDENCIES)
+$(PROG): $(PROCESS) auto_dsig.o $(LIBS) $(MATRIX)
+	$(FC) -o $(PROG) $(PROCESS) $(MATRIX) $(LINKLIBS) $(LDFLAGS) $(BIASDEPENDENCIES)
+
+$(PROG)_forhel: $(PROCESS) auto_dsig.o $(LIBS) $(MATRIX_HEL)
+	$(FC) -o $(PROG)_forhel $(PROCESS) $(MATRIX_HEL) $(LINKLIBS) $(LDFLAGS) $(BIASDEPENDENCIES)
 
 gensym: $(SYMMETRY) configs.inc $(LIBDIR)libmodel.$(libext) $(LIBDIR)libgeneric.$(libext)
 	$(FC) -o gensym $(SYMMETRY) -L../../lib/ -lmodel -lgeneric -lpdf $(llhapdf) $(LDFLAGS)
@@ -56,6 +65,8 @@
 	cd ../../Source/PDF; make
 
 # Add source so that the compiler finds the DiscreteSampler module.
+$(MATRIX): %.o: %.f
+	$(FC) $(FFLAGS) $(MATRIX_FLAG) -c $< -I../../Source/
 %.o: %.f
 	$(FC) $(FFLAGS) -c $< -I../../Source/
 
@@ -78,4 +89,4 @@
 initcluster.o: message.inc
 
 clean:
-	$(RM) *.o gensym madevent
+	$(RM) *.o gensym madevent madevent_forhel

=== modified file 'Template/LO/SubProcesses/myamp.f'
--- Template/LO/SubProcesses/myamp.f	2019-12-04 21:34:25 +0000
+++ Template/LO/SubProcesses/myamp.f	2021-03-25 20:41:08 +0000
@@ -38,7 +38,8 @@
 c
       include 'maxamps.inc'
       integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest,	tstrategy
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
       common/to_sprop/sprop,tprid
@@ -77,7 +78,8 @@
          include 'props.inc'
          nbw = 0
          do i=-1,-(nexternal-3),-1
-            if (iforest(1,i,iconfig) .eq. 1 .or. prwidth(i,iconfig).le.0) then
+            if (iforest(1,i,iconfig) .eq. 1 .or. prwidth(i,iconfig).le.0.or.
+     &                       (nincoming.eq.2.and.iforest(1,i,iconfig) .eq. 2)) then
               cycle
             endif
             nbw=nbw+1
@@ -110,7 +112,7 @@
 c     Start loop over propagators
       do i=-1,-(nexternal-3),-1
          onbw(i) = .false.
-         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
+         if (iforest(1,i,iconfig) .eq. 1.or.(nincoming.eq.2.and.iforest(1,i,iconfig).eq.2)) tsgn=-1d0
          do j=0,3
             xp(j,i) = xp(j,iforest(1,i,iconfig))
      $           +tsgn*xp(j,iforest(2,i,iconfig))
@@ -123,9 +125,12 @@
 c            write(*,*) 'xmass',xmass,prmass(i,iconfig)
 c
 c           Here we set if the BW is "on-shell" for LesHouches
-c            
-            prwidth_tmp(i,iconfig) = max(prwidth(i,iconfig), prmass(i,iconfig)*small_width_treatment)
- 
+c
+            if (prwidth(i,iconfig).gt.0) then
+               prwidth_tmp(i,iconfig) = max(prwidth(i,iconfig), prmass(i,iconfig)*small_width_treatment)
+            else
+               prwidth_tmp(i,iconfig) = 0d0
+            endif
             onshell = (abs(xmass - prmass(i,iconfig)) .lt.
      $           bwcutoff*prwidth_tmp(i,iconfig).and.
      $           (prwidth_tmp(i,iconfig)/prmass(i,iconfig).lt.0.1d0.or.
@@ -247,7 +252,8 @@
       common/to_smin/ Smin
 
       integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest,	tstrategy
 
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
@@ -310,7 +316,11 @@
 c      etmin = 10
       nt = 0
       do i = -nexternal,-1
-         prwidth_tmp(i,iconfig) = max(prwidth(i,iconfig), prmass(i,iconfig)*small_width_treatment)
+         if (prwidth(i,iconfig) .gt.0d0)then
+            prwidth_tmp(i,iconfig) = max(prwidth(i,iconfig), prmass(i,iconfig)*small_width_treatment)
+         else
+            prwidth_tmp(i,iconfig) = 0d0
+         endif
       enddo
 
 
@@ -354,8 +364,10 @@
 
 c     Start loop over propagators
       do i=-1,-(nexternal-3),-1
-         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
-         if (tsgn .eq. 1d0) then                         !s channel
+         if (iforest(1,i,iconfig) .eq. 1.or.iforest(1,i,iconfig) .eq. 2)then
+              tsgn=-1d0
+         endif
+         if (tsgn .eq. 1d0) then !s channel
             xm(i) = xm(iforest(1,i,iconfig))+xm(iforest(2,i,iconfig))
             xe(i) = xe(iforest(1,i,iconfig))+xe(iforest(2,i,iconfig))
             mtot = mtot - xm(i)

=== removed file 'Template/LO/SubProcesses/projection.f'
=== 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-03-25 20:41:08 +0000
@@ -590,6 +590,8 @@
       common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
       double precision stot,m1,m2
       common/to_stot/stot,m1,m2
+      logical init_mode
+      common/to_determine_zero_hel/init_mode
 
 C   local variables
       integer i, j, idi, idj, k,m
@@ -597,6 +599,7 @@
       real*8 PI
       parameter( PI = 3.14159265358979323846d0 )
       integer iforest(2,-max_branch:-1,lmaxconfigs)
+      integer tstrategy
       double precision asref, pt2prev(n_max_cl),pt2min
       integer n, ibeam(2), iqcd(0:2)
       integer idfl, idmap(-nexternal:nexternal)
@@ -650,11 +653,12 @@
       endif
  100  clustered = cluster(p(0,1))
       if(.not.clustered) then
+         if(init_mode) goto 999
          open(unit=26,file='../../../error',status='unknown',err=999)
          write(26,*) 'Error: Clustering failed in cluster.f.'
          write(*,*) 'Error: Clustering failed in cluster.f.'
          stop
- 999     write(*,*) 'error'
+ 999     write(*,*) 'error for clustering'
          setclscales=.false.
          clustered = .false.
          return
@@ -1349,6 +1353,7 @@
       logical setclscales
       integer mapconfig(0:lmaxconfigs), this_config
       integer iforest(2,-max_branch:-1,lmaxconfigs)
+      integer tstrategy(lmaxconfigs)
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
       include 'configs.inc'
@@ -1738,9 +1743,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/setcuts.f'
--- Template/LO/SubProcesses/setcuts.f	2019-06-17 14:27:09 +0000
+++ Template/LO/SubProcesses/setcuts.f	2021-03-25 20:41:08 +0000
@@ -879,7 +879,8 @@
       double precision pmass(nexternal)
       common/to_mass/  pmass
       integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest,	tstrategy
       integer mapconfig(0:lmaxconfigs), this_config
       common/to_mconfigs/mapconfig, this_config
 
@@ -946,7 +947,8 @@
 c     Global
 c
       integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest,	tstrategy
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
       common/to_sprop/sprop,tprid

=== modified file 'Template/LO/SubProcesses/setscales.f'
--- Template/LO/SubProcesses/setscales.f	2020-03-08 20:49:57 +0000
+++ Template/LO/SubProcesses/setscales.f	2021-03-25 20:41:08 +0000
@@ -24,6 +24,7 @@
       double precision pmass(nexternal)
       common/to_mass/  pmass
 
+      
       real*8 xptj,xptb,xpta,xptl,xmtc
       real*8 xetamin,xqcut,deltaeta
       common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
@@ -69,8 +70,8 @@
           enddo
           rscale=rscale/2d0
       elseif(dynamical_scale_choice.eq.4) then
-c         \sqrt(s), partonic energy
-          rscale=dsqrt(max(0d0,2d0*dot(P(0,1),P(0,2))))
+c     \sqrt(s), partonic energy
+         rscale = dsqrt(max(0d0,sumdot(p(0,1), p(0,2), 1d0)))
       elseif(dynamical_scale_choice.eq.5) then
 c         \decaying particle mass, for decays
           rscale=dsqrt(max(0d0,dot(P(0,1),P(0,1))))

=== removed file 'Template/LO/SubProcesses/shrinktops.f'
=== removed file 'Template/LO/SubProcesses/transform.f'
=== removed file 'Template/LO/SubProcesses/transformint.f'
=== 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-03-25 20:41:08 +0000
@@ -145,11 +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
-         get_betaz = (pin(3) * pout(0) - pout(3) * pin(0)) / denom
+         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
       else if (pin(0).eq.pin(3)) then
          get_betaz = (pin(0)**2 - pout(0)**2)/(pin(0)**2 + pout(0)**2)
       else if (pin(0).eq.abs(pin(3))) then
@@ -460,6 +468,8 @@
       include 'run.inc'
       include 'run_config.inc'
 
+
+      
 c
 c     Arguments
 c
@@ -496,6 +506,9 @@
 C     
 C     GLOBAL
 C
+      logical init_mode
+      common/to_determine_zero_hel/init_mode
+      
       double precision twgt, maxwgt,swgt(maxevents)
       integer                             lun, nw, itmin
       common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
@@ -543,7 +556,7 @@
 C-----
       
       if ((nw .ge. maxevents).and.do_write_events) return
-
+      if (init_mode) return
 C     if all the necessary inputs to write the events have already been
 C     computed in the bias module, then directly jump to write_events
       if (AlreadySetInBiasModule) then

=== 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-03-25 20:41:08 +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-03-25 20:41:08 +0000
@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 ################################################################################
 #
@@ -18,11 +18,20 @@
 """This is the main executable, a simple frontend to set up the PYTHONPATH
 and call immediately the command line interface scripts"""
 
+from __future__ import absolute_import
+from __future__ import print_function
 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
 
@@ -62,7 +71,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
@@ -81,13 +90,12 @@
 # Call the cmd interface main loop
 
 
-
 try:
     cmd_line = cmd_interface.GridPackCmd(me_dir=root_path, nb_event=args[0], seed=args[1], gran=args[2])            
 except KeyboardInterrupt:
-    print 'Quit on KeyboardInterrupt' 
+    print('Quit on KeyboardInterrupt') 
 
-print 'DONE'
+print('DONE')
 
 try:
     # Remove lock file

=== 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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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/FKS_params.dat'
--- Template/NLO/Cards/FKS_params.dat	2020-05-13 10:13:31 +0000
+++ Template/NLO/Cards/FKS_params.dat	2021-03-25 20:41:08 +0000
@@ -136,7 +136,7 @@
 ! this to '-1' if you want to forbid the use of this method
 ! altogether.
 #NHelForMCoverHels
--1
+4
 ! Default :: 4
 !
 ! This parameter sets for which fraction of the events the virtual

=== 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-03-25 20:41:08 +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-interpolation grid     *
+# (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/Cards/shower_card.dat'
--- Template/NLO/Cards/shower_card.dat	2017-08-15 19:34:40 +0000
+++ Template/NLO/Cards/shower_card.dat	2021-03-25 20:41:08 +0000
@@ -57,7 +57,8 @@
 # FxFx merging parameters                             !ONLY FOR PYTHIA8!
 #***********************************************************************
 Qcut       = -1.0     # Merging scale
-njmax      = 0        # Maximal multiplicity in the merging
+njmax      = -1       # Maximal multiplicity in the merging.
+#                     # -1 means guessed  from the process definition
 #***********************************************************************
 # Decay channels                                                       *
 #***********************************************************************
@@ -127,6 +128,8 @@
                                  # Dir names separated by white spaces
 ANALYSE      =                   # User's analysis and histogramming
                                  # routines (please use .o as extension
-                                 # and use spaces to separate files)
+                                 # and use spaces to separate files).
+                                 # If the HwU.o files has to be linked,
+                                 # it should be put *first*.
 
 

=== 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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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/Scripts/MCatNLO_MadFKS_HERWIG6.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIG6.Script	2014-11-10 09:33:16 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIG6.Script	2021-03-25 20:41:08 +0000
@@ -141,7 +141,7 @@
 EXEC_NAME=
 EXEC_NAMEF=
 EXEC_NAMEFF=
-HWUTI+=' jimmy.o'
+HWUTI+=' pineappl_interface_dummy.o jimmy.o HwU.o'
 read -a HWUTIA <<< "$HWUTI"
 HWUTI=
 for i in "${HWUTIA[@]}"; do HWUTI+=`pwd`"/objects/$i "; done

=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6PT.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6PT.Script	2014-11-10 09:33:16 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6PT.Script	2021-03-25 20:41:08 +0000
@@ -141,6 +141,7 @@
 EXEC_NAME=
 EXEC_NAMEF=
 EXEC_NAMEFF=
+PYUTI+=' pineappl_interface_dummy.o HwU.o'
 read -a PYUTIA <<< "$PYUTI"
 PYUTI=
 for i in "${PYUTIA[@]}"; do PYUTI+=`pwd`"/objects/$i "; done

=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6Q.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6Q.Script	2014-11-10 09:33:16 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6Q.Script	2021-03-25 20:41:08 +0000
@@ -141,6 +141,7 @@
 EXEC_NAME=
 EXEC_NAMEF=
 EXEC_NAMEFF=
+PYUTI+=' pineappl_interface_dummy.o HwU.o'
 read -a PYUTIA <<< "$PYUTI"
 PYUTI=
 for i in "${PYUTIA[@]}"; do PYUTI+=`pwd`"/objects/$i "; done

=== 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-03-25 20:41:08 +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;
 	  }
 	}

=== added file 'Template/NLO/MCatNLO/include/orders.inc'
--- Template/NLO/MCatNLO/include/orders.inc	1970-01-01 00:00:00 +0000
+++ Template/NLO/MCatNLO/include/orders.inc	2021-03-25 20:41:08 +0000
@@ -0,0 +1,5 @@
+C WARNING: this is a dummy file to make the code for showering compile
+C DO NOT USE ELSEWHERE
+      integer amp_split_size
+      parameter (amp_split_size = 2)
+C =================================================================

=== added symlink 'Template/NLO/MCatNLO/include/pineappl_common.inc'
=== target is '../../SubProcesses/pineappl_common.inc'
=== added symlink 'Template/NLO/MCatNLO/include/pineappl_maxproc.inc'
=== target is '../../SubProcesses/pineappl_maxproc.inc'
=== added symlink 'Template/NLO/MCatNLO/include/reweight_pineappl.inc'
=== target is '../../SubProcesses/reweight_pineappl.inc'
=== added symlink 'Template/NLO/MCatNLO/srcCommon/pineappl_interface_dummy.f'
=== target is '../../SubProcesses/pineappl_interface_dummy.f'
=== 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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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-03-25 20:41:08 +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	2017-07-31 08:48:20 +0000
+++ Template/NLO/Source/extra_weights.f	2021-03-25 20:41:08 +0000
@@ -1,18 +1,18 @@
       module extra_weights
       
          integer,parameter :: iwgtinfo=-5,maxscales=9,maxPDFs=200
-     $     ,maxPDFsets=25,maxdynscales=10
+     $     ,maxPDFsets=25,maxdynscales=10,maxorders=10
          integer :: max_mom_str=1,max_mext=1,max_n_ctr=1
          logical :: doreweight,lscalevar(maxdynscales)
      $        ,lpdfvar(maxPDFsets)
          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)
-     $        ,wgtxsecmu(maxscales,maxscales,maxdynscales)
+     $        ,wgtxsecmu(maxorders,maxscales,maxscales,maxdynscales)
      $        ,wgtxsecPDF(0:maxPDFs,maxPDFsets),wgtbpower,wgtcpower
      $        ,veto_multiplier,H1_factor_virt,veto_compensating_factor
      $        ,born_wgt_veto

=== 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-03-25 20:41:08 +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/run_printout.f'
--- Template/NLO/Source/run_printout.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/run_printout.f	2021-03-25 20:41:08 +0000
@@ -21,12 +21,17 @@
       write(6,*)  'Collider parameters:'
       write(6,*)  '--------------------'
 
+      ab(1) = '?'
+      ab(2) = '?'
       do i=1,2
          IF(LPP(i).EQ. 0) ab(i)='e'
          IF(LPP(i).EQ. 1) ab(i)='P'
          IF(LPP(i).EQ.-1) ab(i)='Pb'
          IF(LPP(i).EQ.2) ab(i)='a2'
-         IF(LPP(i).EQ.3) ab(i)='a3'
+         IF(LPP(i).EQ.3) ab(i)='e-'
+         IF(LPP(i).EQ.-3) ab(i)='e+'
+         IF(LPP(i).EQ.4) ab(i)='m-'
+         IF(LPP(i).EQ.-4) ab(i)='m+'         
       enddo
 
       ene=2d0*dsqrt(ebeam(1)*ebeam(2))

=== 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-03-25 20:41:08 +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
@@ -114,7 +114,7 @@
           write(*,*) 'The default order of alpha_s running is fixed to '
      &         ,nloop
       endif
-      if (nlo_ps) then