Merge lp:~maddevelopers/mg5amcnlo/2.2.3 into lp:mg5amcnlo/lts

[Olivier Mattelaer]

[lp:~maddevelopers/mg5amcnlo/2.2.3]

341. By Olivier Mattelaer on 2015-02-03

fix a problem in add model with particle/anti-particle identification

342. By Olivier Mattelaer on 2015-02-03

Improve MadSpin in helicity only mode

343. By Olivier Mattelaer on 2015-02-03

fix misc.gunzip

344. By marco zaro on 2015-02-03

 small fix/change (suggested by J. Bendavid) about when jobs are split

345. By Rikkert Frederix on 2015-02-04

Fix in the definition of the total cross section when generating a HepMC file with pythia8.

346. By Rikkert Frederix on 2015-02-04

The NLO code now stops if trying to compute cross section at fixed
s-hat (i.e. without PDFs) but with initial state j_fks.

[Olivier Mattelaer (olivier-mattelaer)]

Hi guys, I'm starting the merge.

All the following revision in this branch will not be part of the official release.

Cheers,

Olivier

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rik,

Could you be more specific for the following two input in the UpdateNotes:
 RF: If a NAN is found, the code now skips that PS point and continues.
 RF: For fNLO runs the virtuals were included twice in the setting of the integration grids.

Especially how critical is the second. (i.e. can previous computation can be trusted)

Cheers,

Olivier

[Rikkert Frederix (frederix)]

Hi Olivier,

> RF: If a NAN is found, the code now skips that PS point and continues.

Due to numerical inaccuracies, very rarely you might get a NAN from the phase-space generation or somewhere else for a single PS point. Before the code added this NAN phase-space points to the result, which meant that all became NANs. Now it simply skips this very rare PS point and continues correctly.

> RF: For fNLO runs the virtuals were included twice in the setting of the
> integration grids.

This is nothing serious: it was only in the setting of the grids, not in the actual results. Moreover, due to the virt-tricks the actual contribution from the virtuals is nearly always negligible anyway.

Cheers,
Rik

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2014-10-22 13:41:57 +0000
+++ MadSpin/decay.py	2015-02-04 13:27:21 +0000
@@ -2003,9 +2003,10 @@
         # generate BR and all the square matrix element based on the banner.
         self.generate_all_matrix_element()
 
-        resonances = self.width_estimator.resonances
-        logger.debug('List of resonances: %s' % resonances)
-        self.extract_resonances_mass_width(resonances) 
+        if not self.options["onlyhelicity"]:
+            resonances = self.width_estimator.resonances
+            logger.debug('List of resonances: %s' % resonances)
+            self.extract_resonances_mass_width(resonances) 
 
         self.compile()
         
@@ -2056,6 +2057,8 @@
             for key in self.all_ME:
                 for mode in self.all_ME['decays']:
                     mode['max_weight'] = max_weight_arg
+        elif self.options["onlyhelicity"]:
+            logger.info("not needed in helicity mode")
         else:
             #self.get_max_weight_from_1toN()
             self.get_max_weight_from_event(decay_mapping)  
@@ -2093,7 +2096,8 @@
         if not self.options['ms_dir']:
             shutil.rmtree(pjoin(self.path_me,'production_me'))
             shutil.rmtree(pjoin(self.path_me,'full_me'))
-            shutil.rmtree(pjoin(self.path_me,'decay_me'))
+            if not self.options["onlyhelicity"]:
+                shutil.rmtree(pjoin(self.path_me,'decay_me'))
         # set the environment variable GFORTRAN_UNBUFFERED_ALL 
         # to its original value
         #os.environ['GFORTRAN_UNBUFFERED_ALL']='n'
@@ -2137,7 +2141,7 @@
         self.write_banner_information()
         
         
-        event_nb = 0
+        event_nb, fail_nb = 0, 0
         nb_skip = 0 
         trial_nb_all_events=0
         starttime = time.time()
@@ -2146,34 +2150,23 @@
             production_tag, event_map = self.load_event()
             if production_tag == 0 == event_map: #end of file
                 break
- 
-            numdecays=0
-            for decay in self.all_ME[production_tag]['decays']:
-              if decay['decay_tag']:
-                  numdecays = numdecays + 1
-
-            if not numdecays:
-                #no decays for this production mode, run in passthrough mode, only adding the helicities to the events
-                nb_mc_masses=0
-                p, p_str=self.curr_event.give_momenta(event_map)
-                stdin_text=' %s %s %s %s \n' % ('2', self.options['BW_cut'], self.Ecollider, 1.0)
-                stdin_text+=p_str
-                # here I also need to specify the Monte Carlo Masses
-                stdin_text+=" %s \n" % nb_mc_masses
-                
-                mepath = self.all_ME[production_tag]['path']
-                mepath = mepath.replace('production_me','full_me');
-                                
-                trial_nb, BWvalue, weight, momenta, failed, use_mc_masses, helicities = self.loadfortran( 'unweighting', mepath, stdin_text)
-                                
-                self.reset_helicityonly_in_prod_event(event_map, helicities)
-                event_nb+=1
-                decayed_event = self.curr_event
-                self.outputfile.write(decayed_event.string_event())
-                #print "number of trials: "+str(trial_nb)
+
+            if  event_nb and \
+                (event_nb % max(int(10**int(math.log10(float(event_nb)))),1000)==0): 
+                running_time = misc.format_timer(time.time()-starttime)
+                logger.info('Event nb %s %s' % (event_nb, running_time))
+                self.mscmd.update_status(('$events',1,event_nb, 'decaying events'), 
+                                         force=False, print_log=False)
+            if (event_nb==10001): logger.info('reducing number of print status. Next status update in 10000 events')
+
+
+            if self.options["onlyhelicity"]:
+                trial_nb, fail = self.adding_only_helicity(event_map, production_tag)
                 trial_nb_all_events+=trial_nb
+                fail_nb += fail
+                event_nb += 1
                 continue
-            
+
             # Here we need to select a decay configuration on a random basis:
             decay = self.all_ME.get_random_decay(production_tag)
             if not decay['decay_tag']:
@@ -2191,12 +2184,6 @@
                     
                 event_nb+=1
                 report[decay['decay_tag']] += 1 
-                if (event_nb % max(int(10**int(math.log10(float(event_nb)))),1000)==0): 
-                    running_time = misc.format_timer(time.time()-starttime)
-                    logger.info('Event nb %s %s' % (event_nb, running_time))
-                    self.mscmd.update_status(('$events',1,event_nb, 'decaying events'), 
-                                             force=False, print_log=False)
-                if (event_nb==10001): logger.info('reducing number of print status. Next status update in 10000 events')
 
             indices_for_mc_masses, values_for_mc_masses=self.get_montecarlo_masses_from_event(decay['decay_struct'], event_map, decay['prod2full'])
             nb_mc_masses=len(indices_for_mc_masses)
@@ -2211,7 +2198,12 @@
                 stdin_text+='%s  \n' % str(values_for_mc_masses).strip('[]').replace(',', ' ')
             
 #            here apply the reweighting procedure in fortran
-            trial_nb, BWvalue, weight, momenta, failed, use_mc_masses, helicities = self.loadfortran( 'unweighting', decay_me['path'], stdin_text)
+            output = self.loadfortran( 'unweighting', decay_me['path'], stdin_text)
+            if not output:
+                fail_nb +=1
+                continue
+            trial_nb, BWvalue, weight, momenta, failed, use_mc_masses, helicities = output 
+            
             # next: need to fill all intermediate momenta
             if nb_mc_masses>0 and use_mc_masses==0:nb_fail_mc_mass+=1
             
@@ -2312,9 +2304,43 @@
         logger.info('Number of events with weights larger than max_weight: %s' % report['over_weight'])
         logger.info('Number of subprocesses '+str(len(self.calculator)))
         logger.info('Number of failures when restoring the Monte Carlo masses: %s ' % nb_fail_mc_mass)
+        if fail_nb:
+            logger.info('Number of failures in reshuffling (event skipped): %s ' % fail_nb)
         
         return  event_nb/(event_nb+nb_skip)       
 
+
+    def adding_only_helicity(self, event_map, production_tag):
+        """no decays for this production mode, run in passthrough mode,
+           only adding the helicities to the events """
+
+        #no decays for this production mode, run in passthrough mode, only adding the helicities to the events
+        nb_mc_masses=0
+        p, p_str=self.curr_event.give_momenta(event_map)
+        stdin_text=' %s %s %s %s \n' % ('2', self.options['BW_cut'], self.Ecollider, 1.0)
+        stdin_text+=p_str
+        # here I also need to specify the Monte Carlo Masses
+        stdin_text+=" %s \n" % nb_mc_masses
+        
+        mepath = self.all_ME[production_tag]['path']
+        decay = self.all_ME[production_tag]['decays'][0]
+        decay_me=self.all_ME.get_decay_from_tag(production_tag, decay['decay_tag'])
+        mepath = decay_me['path']
+                        
+        output = self.loadfortran( 'unweighting', mepath, stdin_text)
+        if not output:
+            # Event fail
+            return 0, 1
+        trial_nb, BWvalue, weight, momenta, failed, use_mc_masses, helicities = output                
+        self.reset_helicityonly_in_prod_event(event_map, helicities)
+
+        decayed_event = self.curr_event
+        self.outputfile.write(decayed_event.string_event())
+        #print "number of trials: "+str(trial_nb)
+        
+        return trial_nb, 0
+
+
     def get_int_mom_in_decay(self,decay_struct,ext_mom):
         """  fill  """
         momenta_in_decay={}
@@ -2569,7 +2595,8 @@
         path_me = self.path_me        
         
         # 1. compute the partial width------------------------------------------
-        self.get_branching_ratio()
+        if not self.options["onlyhelicity"]:
+            self.get_branching_ratio()
         
         # 2. compute the production matrix element -----------------------------
         processes = [line[9:].strip() for line in self.banner.proc_card
@@ -2592,11 +2619,22 @@
         start = time.time()
         commandline=''
         for proc in processes:
+            # deal with @ syntax need to move it after the decay specification
+            if '@' in proc:
+                proc, proc_nb = proc.split('@')
+                try:
+                    proc_nb = int(proc_nb)
+                except ValueError:
+                    raise MadSpinError, 'MadSpin didn\'t allow order restriction after the @ comment: \"%s\" not valid' % proc_nb
+                proc_nb = '@ %i' % proc_nb 
+            else:
+                proc_nb = ''
+            
             if '[' not in proc:
                 commandline+="add process %s  --no_warning=duplicate;" % proc
             else:
                 process, order, final = re.split('\[\s*(.*)\s*\]', proc)
-                commandline+="add process %s  --no_warning=duplicate;" % (process)
+                commandline+="add process %s %s --no_warning=duplicate;" % (process, proc_nb)
                 if not order:
                     continue
                 elif not order.startswith('virt='):
@@ -2620,14 +2658,14 @@
                     result = re.split('([/$@]|\w+=\w+)', process, 1)                    
                     if len(result) ==3:
                         process, split, rest = result
-                        commandline+="add process %s pert_%s %s%s  --no_warning=duplicate;" % (process, order ,split, rest)
+                        commandline+="add process %s pert_%s %s%s %s --no_warning=duplicate;" % (process, order ,split, rest, proc_nb)
                     else:
-                        commandline +='add process %s pert_%s --no_warning=duplicate;' % (process,order)                                       
+                        commandline +='add process %s pert_%s  %s --no_warning=duplicate;' % (process,order, proc_nb)                                       
         commandline = commandline.replace('add process', 'generate',1)
         logger.info(commandline)
         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(self.list_branches.keys()))        
         mgcmd.exec_cmd(commandline, precmd=True)        
         logger.info('Done %.4g' % (time.time()-start))
 
@@ -2651,90 +2689,97 @@
         for key in 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]
-                else:
-                  logger.info('keeping dummy decay for passthrough mode')
+                    del self.list_branches[key]
+                elif not self.options["onlyhelicity"]:
+                    raise Exception, " No decay define for process."
+                    logger.info('keeping dummy decay for passthrough mode')
 
         # 4. compute the full matrix element -----------------------------------
-        logger.info('generating the full square matrix element (with decay)')
-        start = time.time()
-        to_decay = self.mscmd.list_branches.keys()
-        decay_text = []
-        for decays in self.mscmd.list_branches.values():
-            for decay in  decays:
-                if '=' not in decay:
-                    decay += ' QCD=99'
-                if ',' in decay:
-                    decay_text.append('(%s)' % decay)
-                else:
-                    decay_text.append(decay)
-        decay_text = ', '.join(decay_text)
-        commandline = ''
-        
-        
-        for proc in processes:
-            # deal with @ syntax need to move it after the decay specification
-            if '@' in proc:
-                proc, proc_nb = proc.split('@')
-                try:
-                    proc_nb = int(proc_nb)
-                except ValueError:
-                    raise MadSpinError, 'MadSpin didn\'t allow order restriction after the @ comment: \"%s\" not valid' % proc_nb
-                proc_nb = '@ %i' % proc_nb 
-            else:
-                proc_nb = '' 
-            
-            if '[' not in proc:
-                nb_comma = proc.count(',')
-                if nb_comma == 0:
-                    commandline+="add process %s, %s %s  --no_warning=duplicate;" % (proc, decay_text, proc_nb)
-                elif nb_comma == 1:
-                    before, after = proc.split(',')
-                    commandline+="add process %s, %s, (%s, %s) %s  --no_warning=duplicate;" % (before, decay_text, after, decay_text, proc_nb)
-                else:
-                    raise Exception, 'too much decay at MG level. this can not be done for the moment)'
-            else:
-                process, order, final = re.split('\[\s*(.*)\s*\]', proc)
-                commandline+="add process %s, %s %s  --no_warning=duplicate;" % (process, decay_text, proc_nb)
-                if not order:
-                    continue
-                elif not order.startswith('virt='):
-                    if '=' in order:
-                        order = order.split('=',1)[1]
-                    # define the list of particles that are needed for the radiateion
-                    pert = fks_common.find_pert_particles_interactions(
-                         mgcmd._curr_model,pert_order = order)['soft_particles']
-                    commandline += "define pert_%s = %s;" % (order, ' '.join(map(str,pert)) )
-                    
-                    # check if we have to increase by one the born order
-                    if '%s=' % order in process:
-                        result=re.split(' ',process)
-                        process=''
-                        for r in result:
-                            if '%s=' % order in r:
-                                ior=re.split('=',r)
-                                r='QCD=%i' % (int(ior[1])+1)
-                            process=process+r+' '
-                    #handle special tag $ | / @
-                    result = re.split('([/$@]|\w+=\w+)', process, 1)                    
-                    if len(result) ==3:
-                        process, split, rest = result
-                        commandline+="add process %s pert_%s %s%s , %s %s --no_warning=duplicate;" % \
-                              (process, order, split, rest, decay_text, proc_nb)
-                    else:
-                        commandline +='add process %s pert_%s, %s %s  --no_warning=duplicate;' % \
-                                           (process, order, decay_text, proc_nb)
-                    
+        if not self.options["onlyhelicity"]:
+            logger.info('generating the full square matrix element (with decay)')
+            start = time.time()
+            to_decay = self.mscmd.list_branches.keys()
+            decay_text = []
+            for decays in self.mscmd.list_branches.values():
+                for decay in  decays:
+                    if '=' not in decay:
+                        decay += ' QCD=99'
+                    if ',' in decay:
+                        decay_text.append('(%s)' % decay)
+                    else:
+                        decay_text.append(decay)
+            decay_text = ', '.join(decay_text)
+            commandline = ''
+            
+            
+            for proc in processes:
+                # deal with @ syntax need to move it after the decay specification
+                if '@' in proc:
+                    proc, proc_nb = proc.split('@')
+                    try:
+                        proc_nb = int(proc_nb)
+                    except ValueError:
+                        raise MadSpinError, 'MadSpin didn\'t allow order restriction after the @ comment: \"%s\" not valid' % proc_nb
+                    proc_nb = '@ %i' % proc_nb 
+                else:
+                    proc_nb = '' 
+                
+                if '[' not in proc:
+                    nb_comma = proc.count(',')
+                    if nb_comma == 0:
+                        commandline+="add process %s, %s %s  --no_warning=duplicate;" % (proc, decay_text, proc_nb)
+                    elif nb_comma == 1:
+                        before, after = proc.split(',')
+                        commandline+="add process %s, %s, (%s, %s) %s  --no_warning=duplicate;" % (before, decay_text, after, decay_text, proc_nb)
+                    else:
+                        raise Exception, 'too much decay at MG level. this can not be done for the moment)'
+                else:
+                    process, order, final = re.split('\[\s*(.*)\s*\]', proc)
+                    commandline+="add process %s, %s %s  --no_warning=duplicate;" % (process, decay_text, proc_nb)
+                    if not order:
+                        continue
+                    elif not order.startswith('virt='):
+                        if '=' in order:
+                            order = order.split('=',1)[1]
+                        # define the list of particles that are needed for the radiateion
+                        pert = fks_common.find_pert_particles_interactions(
+                             mgcmd._curr_model,pert_order = order)['soft_particles']
+                        commandline += "define pert_%s = %s;" % (order, ' '.join(map(str,pert)) )
                         
-        commandline = commandline.replace('add process', 'generate',1)
-        logger.info(commandline)
-        mgcmd.exec_cmd(commandline, precmd=True)
-        # 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()))
-        mgcmd.exec_cmd(commandline, precmd=True)
-        logger.info('Done %.4g' % (time.time()-start))
+                        # check if we have to increase by one the born order
+                        if '%s=' % order in process:
+                            result=re.split(' ',process)
+                            process=''
+                            for r in result:
+                                if '%s=' % order in r:
+                                    ior=re.split('=',r)
+                                    r='QCD=%i' % (int(ior[1])+1)
+                                process=process+r+' '
+                        #handle special tag $ | / @
+                        result = re.split('([/$@]|\w+=\w+)', process, 1)                    
+                        if len(result) ==3:
+                            process, split, rest = result
+                            commandline+="add process %s pert_%s %s%s , %s %s --no_warning=duplicate;" % \
+                                  (process, order, split, rest, decay_text, proc_nb)
+                        else:
+                            commandline +='add process %s pert_%s, %s %s  --no_warning=duplicate;' % \
+                                               (process, order, decay_text, proc_nb)
+            commandline = commandline.replace('add process', 'generate',1)
+            logger.info(commandline)
+            mgcmd.exec_cmd(commandline, precmd=True)
+            # 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()))
+            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()))        
+            mgcmd.exec_cmd(commandline, precmd=True)        
+            logger.info('Done %.4g' % (time.time()-start))                    
+
 
 
         # 5. add the decay information to the all_topology object --------------                        
@@ -2762,6 +2807,9 @@
                 self.generate_configs_file(nfinal,dico,full_path)
                 
 
+        if self.options["onlyhelicity"]:
+            return
+
         # 6. generate decay only part ------------------------------------------
         logger.info('generate matrix element for decay only (1 - > N).')
         start = time.time()
@@ -2822,9 +2870,10 @@
 
     def compile(self):
         logger.info('Compiling code')
-        self.compile_fortran(self.path_me, mode="production_me")
         self.compile_fortran(self.path_me, mode="full_me")
-        self.compile_fortran(self.path_me, mode="decay_me")
+        if not self.options["onlyhelicity"]:
+            self.compile_fortran(self.path_me, mode="production_me")
+            self.compile_fortran(self.path_me, mode="decay_me")
 
     def compile_fortran(self, path_me, mode='production_me'):
         """ Compile the fortran executables associated with the evalutation of the 
@@ -3212,25 +3261,38 @@
             output = me_value
         elif mode == 'unweighting':
             firstline=external.stdout.readline().split()
-            nexternal=int(firstline[0])
-            trials= int(firstline[1])
-            BWvalue= float(firstline[2])
-            weight= float(firstline[3])
-            failed= float(firstline[4])
-            use_mc_masses=int(firstline[5])
+            try:
+                nexternal=int(firstline[0])
+                trials= int(firstline[1])
+                BWvalue= float(firstline[2])
+                weight= float(firstline[3])
+                failed= float(firstline[4])
+                use_mc_masses=int(firstline[5])
+            except ValueError:
+                logger.debug(firstline)
+                return
             momenta=[external.stdout.readline() for i in range(nexternal)]
             lastline=external.stdout.readline().split()
             helicities=[lastline[i] for i in range(len(lastline))]
             output = trials, BWvalue, weight, momenta, failed, use_mc_masses, helicities
 
-        if len(self.calculator) > 100:
-            logger.debug('more than 100 calculator. Perform cleaning')
+        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.sort()
             cut = max([nb_calls[len(nb_calls)//2], 0.001 * nb_calls[-1]])
             for key, external in list(self.calculator.items()):
                 nb = self.calculator_nbcall[key]
                 if nb < cut:
+                    if key[0]=='full':
+                        path=key[1]
+                        end_signal="5 0 0 0 \n"  # before closing, write down the seed 
+                        external.stdin.write(end_signal)
+                        ranmar_state=external.stdout.readline()
+                        ranmar_file=pjoin(path,'ranmar_state.dat')
+                        ranmar=open(ranmar_file, 'w')
+                        ranmar.write(ranmar_state)
+                        ranmar.close()
                     external.stdin.close()
                     external.stdout.close()
                     external.terminate()
@@ -3244,6 +3306,9 @@
     def calculate_matrix_element(self, mode, production, stdin_text):
         """routine to return the matrix element"""
 
+        if mode != "decay":
+            raise Exception, "This function is only secure in mode decay."
+
         tmpdir = ''
         if (mode, production) in self.calculator:
             external = self.calculator[(mode, production)]

=== modified file 'MadSpin/src/driver.f'
--- MadSpin/src/driver.f	2014-08-20 09:50:23 +0000
+++ MadSpin/src/driver.f	2015-02-04 13:27:21 +0000
@@ -229,10 +229,10 @@
            if (jac.lt.0d0) then
              counter2=counter2+1 
              counter3=counter3+1
-c             if (counter3.gt.500) then 
-c               write(*,*) "500_pts_failed_stop_executation"
-c               stop
-c             endif  
+             if (counter3.gt.500) then
+               write(*,*) "500_pts_failed_stop_executation"
+               goto 1
+             endif
              if (counter2.ge.8) then ! use another topology to generate PS points
                do k=1,n_max_cg
                  amp2(k)=0d0

=== modified file 'Template/LO/Source/cuts.inc'
--- Template/LO/Source/cuts.inc	2013-11-23 21:12:58 +0000
+++ Template/LO/Source/cuts.inc	2015-02-04 13:27:21 +0000
@@ -87,4 +87,7 @@
 C COMMON BLOCK FOR JET MEASURE CUTS: DURHAM KT
       REAL*8 KT_DURHAM, D_PARAMETER
       LOGICAL DO_KT_DURHAM
-      COMMON /JET_MEASURE_CUTS/ KT_DURHAM, D_PARAMETER
\ No newline at end of file
+      COMMON /JET_MEASURE_CUTS/ KT_DURHAM, D_PARAMETER
+      
+      
+      
\ No newline at end of file

=== modified file 'Template/LO/Source/setrun.f'
--- Template/LO/Source/setrun.f	2014-02-23 03:15:09 +0000
+++ Template/LO/Source/setrun.f	2015-02-04 13:27:21 +0000
@@ -194,8 +194,8 @@
      $   10000,
      $   10041,
      $   10042,
-     $   200200,
-     $   200400,
+     $   246800,
+     $   247000,
      $   244600/
 
 

=== modified file 'Template/LO/SubProcesses/myamp.f'
--- Template/LO/SubProcesses/myamp.f	2014-02-20 11:27:50 +0000
+++ Template/LO/SubProcesses/myamp.f	2015-02-04 13:27:21 +0000
@@ -259,9 +259,14 @@
             endif
 c
 c     Here we set onshell for phase space integration (JA 4/8/11)
-c
-            onshell = (abs(xmass - prmass(i,iconfig)) .lt.
+c     For decay-chain syntax use BWcutoff here too (22/12/14)
+            if (gForceBW(i, iconfig).eq.1) then
+               onshell = (abs(xmass - prmass(i,iconfig)) .lt.
+     $           bwcutoff*prwidth(i,iconfig))
+            else
+               onshell = (abs(xmass - prmass(i,iconfig)) .lt.
      $           5d0*prwidth(i,iconfig))
+            endif
 
             if (onshell .and. (lbw(nbw).eq. 2) .or.
      $          .not. onshell .and. (lbw(nbw).eq. 1)) then
@@ -297,6 +302,7 @@
 c
       double precision  xm(-nexternal:nexternal)
       double precision  xe(-nexternal:nexternal)
+      double precision bwcut_for_PS(-nexternal:0)
       double precision tsgn, xo, a
       double precision x1,x2,xk(nexternal)
       double precision dr,mtot,etot,xqfact
@@ -440,25 +446,28 @@
             if (prwidth(i,iconfig) .gt. 0 ) then
                nbw=nbw+1
 c              JA 6/8/2011 Set xe(i) for resonances
-               if (lbw(nbw).eq.1) then
+               if (gforcebw(i,iconfig).eq.1) then
+                  xm(i) = max(xm(i), prmass(i,iconfig)-bwcutoff*prwidth(i,iconfig))
+                  bwcut_for_PS(i) = bwcutoff
+               else if (lbw(nbw).eq.1) then
                   xm(i) = max(xm(i), prmass(i,iconfig)-5d0*prwidth(i,iconfig))
-               else if (gforcebw(i,iconfig).eq.1) then
-                  xm(i) = max(xm(i), prmass(i,iconfig)-bwcutoff*prwidth(i,iconfig))
+                  bwcut_for_PS(i) = 5d0
+               else
+                  bwcut_for_PS(i) = 5d0
                endif
             endif
             xe(i)=max(xe(i),xm(i))
 c     Check for impossible onshell configurations
 c     Either: required onshell and daughter masses too large
 c     Or: forced and daughter masses too large
-c     Or: required offshell and forced, with bwcutoff.le.5
+c     Or: required offshell and forced
             if(prwidth(i,iconfig) .gt. 0.and.
      $         (lbw(nbw).eq.1.and.
-     $          (prmass(i,iconfig)+5d0*prwidth(i,iconfig).lt.xm(i)
-     $           .or.prmass(i,iconfig)-5d0*prwidth(i,iconfig).gt.dsqrt(stot))
+     $          (prmass(i,iconfig)+bwcut_for_PS(i)*prwidth(i,iconfig).lt.xm(i)
+     $           .or.prmass(i,iconfig)-bwcut_for_PS(i)*prwidth(i,iconfig).gt.dsqrt(stot))
      $          .or.gforcebw(i,iconfig).eq.1.and.
      $              prmass(i,iconfig)+bwcutoff*prwidth(i,iconfig).lt.xm(i)
-     $          .or.lbw(nbw).eq.2.and.gforcebw(i,iconfig).eq.1 .and.
-     $              bwcutoff.le.5d0))
+     $          .or.lbw(nbw).eq.2.and.gforcebw(i,iconfig).eq.1))
      $        then
 c     Write results.dat and quit
                call write_null_results()
@@ -478,7 +487,7 @@
                      spole(j)=prmass(i,iconfig)*prmass(i,iconfig)/stot
                      swidth(j) = prwidth(i,iconfig)*prmass(i,iconfig)/stot
                   endif
-               else if((prmass(i,iconfig)+5d0*prwidth(i,iconfig)).ge.xm(i)
+               else if((prmass(i,iconfig)+bwcut_for_PS(i)*prwidth(i,iconfig)).ge.xm(i)
      $                  .and. iden_part(i).eq.0 .or. lbw(nbw).eq.1) then
 c              JA 02/13 Only allow BW if xm below M+5*Gamma
                   write(*,*) 'Setting BW',i,nbw,prmass(i,iconfig)
@@ -496,7 +505,7 @@
 c     Set spmass for BWs
                if (swidth(-i) .ne. 0d0)
      $              spmass=spmass-xm(i) +
-     $              max(xm(i),prmass(i,iconfig)-5d0*prwidth(i,iconfig))
+     $              max(xm(i),prmass(i,iconfig)-bwcut_for_PS(i)*prwidth(i,iconfig))
             else                                  !1/x^pow
               a=prmass(i,iconfig)**2/stot
 c     JA 4/1/2011 always set grid

=== modified file 'Template/MadWeight/Python/MW_param_default.inc'
--- Template/MadWeight/Python/MW_param_default.inc	2014-01-26 14:08:26 +0000
+++ Template/MadWeight/Python/MW_param_default.inc	2015-02-04 13:27:21 +0000
@@ -14,8 +14,8 @@
 mw_run  use_sobol        logical   F                    # don t use sobol generator by default
 mw_run  inputfile        string   './Events/input.lhco' # path to the inputfile
 
-mw_parameter mode logical 1                         
-mw_parameter 1    logical mw_parameter@@mode # choose the mode for the generation of the param_card.dat
+mw_parameter mode integer 1                         
+mw_parameter 1    integer mw_parameter@@mode # choose the mode for the generation of the param_card.dat
 mw_parameter 2    logical 0                  # put on 1 to add new param_card.dat (for step 1) creates the corresponding new events dir if step 5 already performs
 
 mw_perm  permutation    logical T                       # make the permutation between identical particles/jets

=== modified file 'Template/MadWeight/src/driver.f'
--- Template/MadWeight/src/driver.f	2014-01-26 14:08:26 +0000
+++ Template/MadWeight/src/driver.f	2015-02-04 13:27:21 +0000
@@ -36,6 +36,7 @@
       double precision order_value(nb_channel)
       double precision order_error(nb_channel)
       double precision xi_by_channel(100, 20, nb_channel)
+      integer ndo_by_channel(nb_channel)
 
       double precision normalize_perm
 c
@@ -226,6 +227,7 @@
                 enddo
 
              else
+                ndo = ndo_by_channel(integral_index)
                 do i = 1, ndo
                     do j = 1, NDIM
                         xi(i,j) = xi_by_channel(i, j, integral_index)
@@ -330,6 +332,7 @@
                         xi_by_channel(i, j, integral_index) = xi(i,j)
                     enddo
                 enddo
+                ndo_by_channel(integral_index) = ndo
              else
                 order_value(integral_index)=0d0
              endif

=== modified file 'Template/NLO/Cards/shower_card.dat'
--- Template/NLO/Cards/shower_card.dat	2014-08-08 09:58:08 +0000
+++ Template/NLO/Cards/shower_card.dat	2015-02-04 13:27:21 +0000
@@ -103,6 +103,7 @@
 #***********************************************************************
 EXTRALIBS    = stdhep Fmcfio     # Extra-libraries (not LHAPDF) 
                                  # Default: "stdhep Fmcfio"
+                                 # PYTHIA > 8.200 may require library dl
 EXTRAPATHS   = ../lib            # Path to the extra-libraries
                                  # Default: "../lib"
 INCLUDEPATHS =                   # Path to header files needed by c++

=== added symlink 'Template/NLO/MCatNLO/HWPPAnalyzer/fjcore.cc'
=== target is u'../../SubProcesses/fjcore.cc'
=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIGPP.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIGPP.Script	2014-10-29 14:44:37 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIGPP.Script	2015-02-04 13:27:21 +0000
@@ -611,6 +611,8 @@
 #
 #553
 set /Herwig/Particles/Upsilon:Stable Stable
+#551
+set /Herwig/Particles/eta_b:Stable Stable
 
 EOF
 fi

=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script	2014-10-30 15:11:46 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script	2015-02-04 13:27:21 +0000
@@ -581,6 +581,7 @@
 521:maydecay = false               ! stable B hadrons
 531:maydecay = false               ! stable B hadrons
 541:maydecay = false               ! stable B hadrons
+551:maydecay = false               ! stable B hadrons
 553:maydecay = false               ! stable B hadrons
 5112:maydecay = false              ! stable B hadrons
 5122:maydecay = false              ! stable B hadrons

=== added symlink 'Template/NLO/MCatNLO/include/fjcore.hh'
=== target is u'../../SubProcesses/fjcore.hh'
=== modified file 'Template/NLO/MCatNLO/shower_template.sh'
--- Template/NLO/MCatNLO/shower_template.sh	2014-10-30 15:48:03 +0000
+++ Template/NLO/MCatNLO/shower_template.sh	2015-02-04 13:27:21 +0000
@@ -18,7 +18,8 @@
     cd run_$NFILE
     cp -H ../events_$NFILE.lhe events.lhe
     if [ $SHOWER == "PYTHIA8" ] ; then
-        cp ../Pythia8.exe ../Pythia8.cmd ../config.sh .
+        cp ../Pythia8.exe ../Pythia8.cmd .
+        if [ -f ../config.sh ] ; then cp ../config.sh . ; fi
     else
         if [ $SHOWER == "HERWIGPP" ] ; then
             cp ../Herwig++ ../HepMCFortran.so .

=== modified file 'Template/NLO/MCatNLO/srcCommon/mcatnlo_hbook.f'
--- Template/NLO/MCatNLO/srcCommon/mcatnlo_hbook.f	2013-11-20 10:02:06 +0000
+++ Template/NLO/MCatNLO/srcCommon/mcatnlo_hbook.f	2015-02-04 13:27:21 +0000
@@ -1396,6 +1396,8 @@
       YTIT=YU+YTIT0
       WRITE(99,101) XL,YTIT,TITLE(NH)(1:40)
 101   FORMAT('  TITLE ',2(F8.4,1X),'"',A,'"')
+      WRITE(99,151) TITLE(NH)
+151   FORMAT(' (  TITLE    TOP "',A,'"')
       YTIT=YTIT-2.*YTIT0
       WRITE(99,102) XTIT,YTIT,HINT(NH)
 102   FORMAT('  TITLE ',2(F8.4,1X),'" INT=',1PE10.3,'"')

=== modified file 'Template/NLO/MCatNLO/srcCommon/mcatnlo_hbook_gfortran4.f'
--- Template/NLO/MCatNLO/srcCommon/mcatnlo_hbook_gfortran4.f	2013-11-20 10:02:06 +0000
+++ Template/NLO/MCatNLO/srcCommon/mcatnlo_hbook_gfortran4.f	2015-02-04 13:27:21 +0000
@@ -1390,6 +1390,8 @@
       YTIT=YU+YTIT0
       WRITE(99,101) XL,YTIT,TITLE(NH)(1:40)
 101   FORMAT('  TITLE ',2(F8.4,1X),'"',A,'"')
+      WRITE(99,151) TITLE(NH)
+151   FORMAT(' (  TITLE    TOP "',A,'"')
       YTIT=YTIT-2.*YTIT0
       WRITE(99,102) XTIT,YTIT,HINT(NH)
 102   FORMAT('  TITLE ',2(F8.4,1X),'" INT=',1PE10.3,'"')

=== modified file 'Template/NLO/MCatNLO/srcCommon/mcatnlo_hbook_gfortran8.f'
--- Template/NLO/MCatNLO/srcCommon/mcatnlo_hbook_gfortran8.f	2013-11-20 10:02:06 +0000
+++ Template/NLO/MCatNLO/srcCommon/mcatnlo_hbook_gfortran8.f	2015-02-04 13:27:21 +0000
@@ -1327,6 +1327,8 @@
       YTIT=YU+YTIT0
       WRITE(99,101) XL,YTIT,TITLE(NH)(1:40)
 101   FORMAT('  TITLE ',2(F8.4,1X),'"',A,'"')
+      WRITE(99,151) TITLE(NH)
+151   FORMAT(' (  TITLE    TOP "',A,'"')
       YTIT=YTIT-2.*YTIT0
       WRITE(99,102) XTIT,YTIT,HINT(NH)
 102   FORMAT('  TITLE ',2(F8.4,1X),'" INT=',1PE10.3,'"')

=== modified file 'Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f'
--- Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f	2014-08-25 15:05:30 +0000
+++ Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f	2015-02-04 13:27:21 +0000
@@ -236,6 +236,7 @@
         CALL HWUSTA('XI_B-   ')
         CALL HWUSTA('OMEGA_B-')
         CALL HWUSTA('B_C-    ')
+        CALL HWUSTA('ETA_B   ')
         CALL HWUSTA('UPSLON1S')
         CALL HWUSTA('SGM_BBR-')
         CALL HWUSTA('LMD_BBR0')
@@ -350,7 +351,6 @@
       CALL HWABEG
 C---LOOP OVER EVENTS
       DO 100 N=1,MAXEV
- 50      CONTINUE
 C---INITIALISE EVENT
          CALL HWUINE
 C---GENERATE HARD SUBPROCESS
@@ -383,6 +383,7 @@
             WRITE(*,*)'# of events processed=',NEVHEP
 c            CALL HWAEND
          ENDIF
+ 50      CONTINUE
   100 CONTINUE
 C---TERMINATE ELEMENTARY PROCESS
       CALL HWEFIN

=== modified file 'Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f'
--- Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f	2014-06-03 08:44:17 +0000
+++ Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f	2015-02-04 13:27:21 +0000
@@ -155,6 +155,7 @@
             MDCY(PYCOMP(ISIGN*521),1)=0
             MDCY(PYCOMP(ISIGN*531),1)=0
             MDCY(PYCOMP(ISIGN*541),1)=0
+            MDCY(PYCOMP(551),1)=0
             MDCY(PYCOMP(553),1)=0
             MDCY(PYCOMP(ISIGN*5112),1)=0
             MDCY(PYCOMP(ISIGN*5122),1)=0

=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia82.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia82.cc	2014-10-28 13:16:07 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia82.cc	2015-02-04 13:27:21 +0000
@@ -39,10 +39,8 @@
   pythia.readFile(inputname.c_str());
 
   //Create UserHooks pointer for the FxFX matching. Stop if it failed. Pass pointer to Pythia.
-  CombineMatchingInput combined;
-  UserHooks* matching = combined.getHook(pythia);
-  if (!matching) return 1;
-  pythia.setUserHooksPtr(matching);
+  CombineMatchingInput* combined = NULL;
+  UserHooks* matching            = NULL;
 
   pythia.init();
   string filename = pythia.word("Beams:LHEF");
@@ -67,6 +65,13 @@
   //FxFx merging
   bool isFxFx=pythia.flag("JetMatching:doFxFx");
   if (isFxFx) {
+    matching = combined->getHook(pythia);
+    if (!matching) {
+      std::cout << " Failed to initialise jet matching structures.\n"
+                << " Program stopped.";
+      return 1;
+    }
+    pythia.setUserHooksPtr(matching);
     int nJmax=pythia.mode("JetMatching:nJetMax");
     double Qcut=pythia.parm("JetMatching:qCut");
     double PTcut=pythia.parm("JetMatching:qCutME");
@@ -147,7 +152,6 @@
     std::cout << "*********************************************************************** \n";
     std::cout << "*********************************************************************** \n";
   }
-  delete matching;
 
   return 0;
 }

=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia82_hep.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia82_hep.cc	2014-10-28 13:16:07 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia82_hep.cc	2015-02-04 13:27:21 +0000
@@ -17,10 +17,8 @@
   pythia.readFile(inputname.c_str());
 
   //Create UserHooks pointer for the FxFX matching. Stop if it failed. Pass pointer to Pythia.
-  CombineMatchingInput combined;
-  UserHooks* matching = combined.getHook(pythia);
-  if (!matching) return 1;
-  pythia.setUserHooksPtr(matching);
+  CombineMatchingInput* combined = NULL;
+  UserHooks* matching            = NULL;
 
   pythia.init();
 
@@ -35,6 +33,13 @@
   //FxFx merging
   bool isFxFx=pythia.flag("JetMatching:doFxFx");
   if (isFxFx) {
+    matching = combined->getHook(pythia);
+    if (!matching) {
+      std::cout << " Failed to initialise jet matching structures.\n"
+                << " Program stopped.";
+      return 1;
+    }
+    pythia.setUserHooksPtr(matching);
     int nJmax=pythia.mode("JetMatching:nJetMax");
     double Qcut=pythia.parm("JetMatching:qCut");
     double PTcut=pythia.parm("JetMatching:qCutME");
@@ -83,6 +88,7 @@
     } else {
       normhepmc = double(iEventtot) / double(iEventshower);
     }
+    sigmaTotal += evtweight*normhepmc;
     hepmcevt->weights().push_back(evtweight*normhepmc);
     ToHepMC.fill_next_event( pythia, hepmcevt );
     // Add the weight of the current event to the cross section.
@@ -106,7 +112,5 @@
     std::cout << "*********************************************************************** \n";
   }
 
-  delete matching;
-
   return 0;
 }

=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia8_hep.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia8_hep.cc	2014-10-28 14:04:42 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia8_hep.cc	2015-02-04 13:27:21 +0000
@@ -82,6 +82,7 @@
     } else {
       normhepmc = double(iEventtot) / double(iEventshower);
     }
+    sigmaTotal += evtweight*normhepmc;
     hepmcevt->weights().push_back(evtweight*normhepmc);
     ToHepMC.fill_next_event( pythia, hepmcevt );
     // Add the weight of the current event to the cross section.

=== modified file 'Template/NLO/Source/PDF/pdg2pdf.f'
--- Template/NLO/Source/PDF/pdg2pdf.f	2013-12-20 05:39:24 +0000
+++ Template/NLO/Source/PDF/pdg2pdf.f	2015-02-04 13:27:21 +0000
@@ -35,17 +35,18 @@
       include 'pdf.inc'
 C      
       double precision Ctq3df,Ctq4Fn,Ctq5Pdf,Ctq6Pdf,Ctq5L
-      integer mode,Irt,i,j
-      double precision xlast(2),xmulast(2),pdflast(-7:7,2),q2max
-      character*7 pdlabellast(2)
+      integer mode,Irt,i,j,i_replace,ii
+      double precision xlast(20),xmulast(20),pdflast(-7:7,20),q2max
+      character*7 pdlabellast(20)
       double precision epa_electron,epa_proton
-      integer ipart,ireuse,iporg,ihlast(2)
+      integer ipart,ireuse,iporg,ihlast(20)
       save xlast,xmulast,pdflast,pdlabellast,ihlast
-      data xlast/2*-99d9/
-      data xmulast/2*-99d9/
-      data pdflast/30*-99d9/
-      data pdlabellast/2*'abcdefg'/
-      data ihlast/2*-99/
+      data xlast/20*-99d9/
+      data xmulast/20*-99d9/
+      data pdflast/300*-99d9/
+      data pdlabellast/20*'abcdefg'/
+      data ihlast/20*-99/
+      data i_replace/20/
 
 c     Make sure we have a reasonable Bjorken x. Note that even though
 c     x=0 is not reasonable, we prefer to simply return pdg2pdf=0
@@ -76,12 +77,22 @@
       endif
 
       ireuse = 0
-      do i=1,2
-c     Check if result can be reused since any of last two calls
-         if (x.eq.xlast(i) .and. xmu.eq.xmulast(i) .and.
-     $        pdlabel.eq.pdlabellast(i) .and. ih.eq.ihlast(i)) then
-            ireuse = i
+      ii=i_replace
+      do i=1,20
+c     Check if result can be reused since any of last twenty
+c     calls. Start checking with the last call and move back in time
+         if (ih.eq.ihlast(ii)) then
+            if (x.eq.xlast(ii)) then
+               if (xmu.eq.xmulast(ii)) then
+                  if (pdlabel.eq.pdlabellast(ii)) then
+                     ireuse = ii
+                     exit
+                  endif
+               endif
+            endif
          endif
+         ii=ii-1
+         if (ii.eq.0) ii=ii+20
       enddo
 
 c     Reuse previous result, if possible
@@ -92,41 +103,15 @@
          endif
       endif
 
-c     Bjorken x and/or facrorization scale and/or PDF set are not
-c     identical to the saved values: this means a new event and we
-c     should reset everything to compute new PDF values. Also, determine
-c     if we should fill ireuse=1 or ireuse=2.
-      if (ireuse.eq.0.and.xlast(1).ne.-99d9.and.xlast(2).ne.-99d9)then
-         do i=1,2
-            xlast(i)=-99d9
-            xmulast(i)=-99d9
-            do j=-7,7
-               pdflast(j,i)=-99d9
-            enddo
-            pdlabellast(i)='abcdefg'
-            ihlast(i)=-99
-         enddo
-c     everything has been reset. Now set ireuse=1 to fill the first
-c     arrays of saved values below
-         ireuse=1
-      else if(ireuse.eq.0.and.xlast(1).ne.-99d9)then
-c     This is first call after everything has been reset, so the first
-c     arrays are already filled with the saved values (hence
-c     xlast(1).ne.-99d9). Fill the second arrays of saved values (done
-c     below) by setting ireuse=2
-         ireuse=2
-      else if(ireuse.eq.0)then
-c     Special: only used for the very first call to this function:
-c     xlast(i) are initialized as data statements to be equal to -99d9
-         ireuse=1
-      endif
-
+c Calculated a new value: replace the value computed longest ago.
+      i_replace=mod(i_replace,20)+1
+      
 c     Give the current values to the arrays that should be
 c     saved. 'pdflast' is filled below.
-      xlast(ireuse)=x
-      xmulast(ireuse)=xmu
-      pdlabellast(ireuse)=pdlabel
-      ihlast(ireuse)=ih
+      xlast(i_replace)=x
+      xmulast(i_replace)=xmu
+      pdlabellast(i_replace)=pdlabel
+      ihlast(i_replace)=ih
 
       if(iabs(ipart).eq.7.and.ih.gt.1) then
          q2max=xmu*xmu
@@ -135,7 +120,7 @@
          elseif(ih .eq. 2) then !from a proton without breaking
             pdg2pdf=epa_proton(x,q2max)
          endif 
-         pdflast(iporg,ireuse)=pdg2pdf
+         pdflast(iporg,i_replace)=pdg2pdf
          return
       endif
       
@@ -210,11 +195,11 @@
 
          pdg2pdf=Ctq6Pdf(ipart,x,xmu)
       else
-         call pftopdg(ih,x,xmu,pdflast(-7,ireuse))
-         pdg2pdf=pdflast(iporg,ireuse)
+         call pftopdg(ih,x,xmu,pdflast(-7,i_replace))
+         pdg2pdf=pdflast(iporg,i_replace)
       endif      
 
-      pdflast(iporg,ireuse)=pdg2pdf
+      pdflast(iporg,i_replace)=pdg2pdf
       return
       end
 

=== modified file 'Template/NLO/Source/PDF/pdg2pdf_lhapdf.f'
--- Template/NLO/Source/PDF/pdg2pdf_lhapdf.f	2013-12-20 05:39:24 +0000
+++ Template/NLO/Source/PDF/pdg2pdf_lhapdf.f	2015-02-04 13:27:21 +0000
@@ -34,14 +34,16 @@
 C
       include 'pdf.inc'
 C      
-      integer i,j,ihlast(2),ipart,iporg,ireuse,imemlast(2),iset,imem
-      double precision xlast(2),xmulast(2),pdflast(-7:7,2)
+      integer i,j,ihlast(20),ipart,iporg,ireuse,imemlast(20),iset,imem
+     &     ,i_replace,ii
+      double precision xlast(20),xmulast(20),pdflast(-7:7,20)
       save ihlast,xlast,xmulast,pdflast,imemlast
-      data ihlast/2*-99/
-      data xlast/2*-99d9/
-      data xmulast/2*-99d9/
-      data pdflast/30*-99d9/
-      data imemlast/2*-99/
+      data ihlast/20*-99/
+      data xlast/20*-99d9/
+      data xmulast/20*-99d9/
+      data pdflast/300*-99d9/
+      data imemlast/20*-99/
+      data i_replace/20/
 
 c     Make sure we have a reasonable Bjorken x. Note that even though
 c     x=0 is not reasonable, we prefer to simply return pdg2pdf=0
@@ -84,12 +86,22 @@
       call getnmem(iset,imem)
 
       ireuse = 0
-      do i=1,2
-c     Check if result can be reused since any of last two calls
-         if (x.eq.xlast(i) .and. xmu.eq.xmulast(i) .and.
-     $        imem.eq.imemlast(i) .and. ih.eq.ihlast(i)) then
-            ireuse = i
+      ii=i_replace
+      do i=1,20
+c     Check if result can be reused since any of last twenty
+c     calls. Start checking with the last call and move back in time
+         if (ih.eq.ihlast(ii)) then
+            if (x.eq.xlast(ii)) then
+               if (xmu.eq.xmulast(ii)) then
+                  if (imem.eq.imemlast(ii)) then
+                     ireuse = ii
+                     exit
+                  endif
+               endif
+            endif
          endif
+         ii=ii-1
+         if (ii.eq.0) ii=ii+20
       enddo
 
 c     Reuse previous result, if possible
@@ -100,44 +112,18 @@
          endif
       endif
 
-c     Bjorken x and/or facrorization scale and/or PDF set are not
-c     identical to the saved values: this means a new event and we
-c     should reset everything to compute new PDF values. Also, determine
-c     if we should fill ireuse=1 or ireuse=2.
-      if (ireuse.eq.0.and.xlast(1).ne.-99d9.and.xlast(2).ne.-99d9)then
-         do i=1,2
-            xlast(i)=-99d9
-            xmulast(i)=-99d9
-            do j=-7,7
-               pdflast(j,i)=-99d9
-            enddo
-            imemlast(i)=-99
-            ihlast(i)=-99
-         enddo
-c     everything has been reset. Now set ireuse=1 to fill the first
-c     arrays of saved values below
-         ireuse=1
-      else if(ireuse.eq.0.and.xlast(1).ne.-99d9)then
-c     This is first call after everything has been reset, so the first
-c     arrays are already filled with the saved values (hence
-c     xlast(1).ne.-99d9). Fill the second arrays of saved values (done
-c     below) by setting ireuse=2
-         ireuse=2
-      else if(ireuse.eq.0)then
-c     Special: only used for the very first call to this function:
-c     xlast(i) are initialized as data statements to be equal to -99d9
-         ireuse=1
-      endif
+c Calculated a new value: replace the value computed longest ago
+      i_replace=mod(i_replace,20)+1
 
 c     Call lhapdf and give the current values to the arrays that should
 c     be saved
-      call pftopdglha(ih,x,xmu,pdflast(-7,ireuse))
-      xlast(ireuse)=x
-      xmulast(ireuse)=xmu
-      ihlast(ireuse)=ih
-      imemlast(ireuse)=imem
+      call pftopdglha(ih,x,xmu,pdflast(-7,i_replace))
+      xlast(i_replace)=x
+      xmulast(i_replace)=xmu
+      ihlast(i_replace)=ih
+      imemlast(i_replace)=imem
 c
-      pdg2pdf=pdflast(ipart,ireuse);
+      pdg2pdf=pdflast(ipart,i_replace)
       return
       end
 

=== modified file 'Template/NLO/Source/rw_routines.f'
--- Template/NLO/Source/rw_routines.f	2013-08-13 11:57:15 +0000
+++ Template/NLO/Source/rw_routines.f	2015-02-04 13:27:21 +0000
@@ -400,7 +400,7 @@
 c
       integer i,k
 
-      do i=1,20
+      do i=1,len(name)
          k=ichar(name(i:i))
          if(k.ge.65.and.k.le.90) then  !upper case A-Z
             k=ichar(name(i:i))+32   

=== modified file 'Template/NLO/Source/setrun.f'
--- Template/NLO/Source/setrun.f	2014-10-30 11:17:06 +0000
+++ Template/NLO/Source/setrun.f	2015-02-04 13:27:21 +0000
@@ -56,7 +56,6 @@
      &     xmaxup(maxpup),lprup(maxpup)
 c
       include 'nexternal.inc'
-      include 'leshouche_decl.inc'
       logical gridrun,gridpack
       integer          iseed
       common /to_seed/ iseed
@@ -65,6 +64,13 @@
       common /event_normalisation/event_norm
       integer iappl
       common /for_applgrid/ iappl
+C      
+      integer    maxflow
+      parameter (maxflow=999)
+      integer idup(nexternal,maxproc)
+      integer mothup(2,nexternal,maxproc)
+      integer icolup(2,nexternal,maxflow)
+      include 'born_leshouche.inc'
 c
 c----------
 c     start
@@ -74,7 +80,6 @@
 c MZ add the possibility to have shower_MC input lowercase
       call to_upper(shower_MC)
 C
-      call read_leshouche_info(idup_d,mothup_d,icolup_d)
 
 c merging cuts
       xqcut=0d0
@@ -160,7 +165,7 @@
         elseif(lpp(i).eq.-3) then
           idbmup(i)=-11
         elseif(lpp(i).eq.0) then
-          idbmup(i)=idup_d(1,i,1)
+          idbmup(i)=idup(i,1)
         else
           idbmup(i)=lpp(i)
         endif
@@ -219,8 +224,8 @@
      $   10000,
      $   10041,
      $   10042,
-     $   200200,
-     $   200400,
+     $   246800,
+     $   247000,
      $   244600/
 
 

=== modified file 'Template/NLO/SubProcesses/add_write_info.f'
--- Template/NLO/SubProcesses/add_write_info.f	2014-10-30 11:17:06 +0000
+++ Template/NLO/SubProcesses/add_write_info.f	2015-02-04 13:27:21 +0000
@@ -31,8 +31,8 @@
       data firsttime2/.true./
 
 c The process chosen to write
-      integer i_process
-      common/c_addwrite/i_process
+      integer i_process_addwrite
+      common/c_addwrite/i_process_addwrite
       
 c Random numbers
       double precision ran2
@@ -205,7 +205,7 @@
 c This is an (n+1)-body process (see update_unwgt_table in
 c driver_mintMC.f). For S events it corresponds to the underlying Born
 c process chosen
-      ip=i_process
+      ip=i_process_addwrite
       if (ip.lt.1 .or. ip.gt.maxproc_used) then
          write (*,*)'ERROR #12 in add_write_info,'/
      &        /' not a well-defined process',ip,Hevents

=== added file 'Template/NLO/SubProcesses/c_weight.inc'
--- Template/NLO/SubProcesses/c_weight.inc	1970-01-01 00:00:00 +0000
+++ Template/NLO/SubProcesses/c_weight.inc	2015-02-04 13:27:21 +0000
@@ -0,0 +1,12 @@
+*     -*-fortran-*-
+
+      integer max_contr,max_wgt
+      parameter (max_contr=128,max_wgt=256)
+      integer nFKS(max_contr),itype(max_contr),QCDpower(max_contr)
+     &     ,pdg(nexternal,max_contr),icontr,iwgt,plot_id(max_contr)
+      double precision momenta(0:3,nexternal,max_contr),wgt(3,max_contr)
+     &     ,bjx(2,max_contr),scales2(3,max_contr),g_strong(max_contr)
+     &     ,wgts(max_wgt,max_contr),y_bst(max_contr),plot_wgts(max_wgt
+     &     ,max_contr)
+      common /c_weights/momenta,wgt,wgts,plot_wgts,bjx,scales2,g_strong
+     &     ,y_bst,nFKS,itype,QCDpower,pdg,plot_id,icontr,iwgt

=== modified file 'Template/NLO/SubProcesses/cuts.f'
--- Template/NLO/SubProcesses/cuts.f	2014-10-17 19:39:20 +0000
+++ Template/NLO/SubProcesses/cuts.f	2015-02-04 13:27:21 +0000
@@ -399,6 +399,7 @@
       include "nexternal.inc"
       include 'run.inc'
       include 'genps.inc'
+      include 'timing_variables.inc'
       REAL*8 P(0:3,nexternal),rwgt
       integer i,j,istatus(nexternal),iPDG(nexternal)
 c For boosts
@@ -421,6 +422,7 @@
       common /c_leshouche_inc/idup,mothup,icolup
       logical passcuts_user
       external passcuts_user
+      call cpu_time(tBefore)
 c Make sure have reasonable 4-momenta
       if (p(0,1) .le. 0d0) then
          passcuts=.false.
@@ -459,6 +461,8 @@
       enddo
 c Call the actual cuts function  
       passcuts = passcuts_user(pp,istatus,ipdg)
+      call cpu_time(tAfter)
+      t_cuts=t_cuts+(tAfter-tBefore)
       RETURN
       END
 

=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2014-09-02 12:20:56 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2015-02-04 13:27:21 +0000
@@ -347,14 +347,23 @@
 
       call cpu_time(tAfter)
       tTot = tAfter-tBefore
-      tOther = tTot - tOLP - tPDF - tFastJet - tGenPS - tDSigI - tDSigR
-      write(*,*) 'Time spent in clustering : ',tFastJet      
-      write(*,*) 'Time spent in PDF_Engine : ',tPDF
+      tOther = tTot - (tBorn+tGenPS+tReal+tCount+tIS+tFxFx+tf_nb+tf_all
+     &     +t_as+tr_s+tr_pdf+t_plot+t_cuts)
+      write(*,*) 'Time spent in Born : ',tBorn
       write(*,*) 'Time spent in PS_Generation : ',tGenPS
-      write(*,*) 'Time spent in Reals_evaluation: ',tDSigR
-      write(*,*) 'Time spent in IS_evaluation : ',tDSigI
-      write(*,*) 'Time spent in OneLoop_Engine : ',tOLP      
-      write(*,*) 'Time spent in other_tasks : ',tOther
+      write(*,*) 'Time spent in Reals_evaluation: ',tReal
+      write(*,*) 'Time spent in Counter_terms : ',tCount
+      write(*,*) 'Time spent in Integrated_CT : ',tIS-tOLP
+      write(*,*) 'Time spent in Virtuals : ',tOLP      
+      write(*,*) 'Time spent in FxFx_cluster : ',tFxFx
+      write(*,*) 'Time spent in Nbody_prefactor : ',tf_nb
+      write(*,*) 'Time spent in N1body_prefactor : ',tf_all
+      write(*,*) 'Time spent in Adding_alphas_pdf : ',t_as
+      write(*,*) 'Time spent in Reweight_scale : ',tr_s
+      write(*,*) 'Time spent in Reweight_pdf : ',tr_pdf
+      write(*,*) 'Time spent in Filling_plots : ',t_plot
+      write(*,*) 'Time spent in Applying_cuts : ',t_cuts
+      write(*,*) 'Time spent in Other_tasks : ',tOther
       write(*,*) 'Time spent in Total : ',tTot
 
       if(i_momcmp_count.ne.0)then
@@ -375,195 +384,255 @@
       data tDSigI/0.0/
       data tDSigR/0.0/
       data tGenPS/0.0/
+      data tBorn/0.0/
+      data tIS/0.0/
+      data tReal/0.0/
+      data tCount/0.0/
+      data tFxFx/0.0/
+      data tf_nb/0.0/
+      data tf_all/0.0/
+      data t_as/0.0/
+      data tr_s/0.0/
+      data tr_pdf/0.0/
+      data t_plot/0.0/
       end
 
 
-      function sigint(xx,peso,ifl,f)
-c From dsample_fks
+      double precision function sigint(xx,vegas_wgt,ifl,f)
       implicit none
       include 'nexternal.inc'
       include 'mint.inc'
-      real*8 sigint,peso,xx(ndimmax),f(nintegrals)
-      integer ione
-      parameter (ione=1)
-      integer ifl
-      integer ndim
-      common/tosigint/ndim
-      integer           iconfig
-      common/to_configs/iconfig
-      integer i
-      double precision wgt,dsig,ran2,rnd
-      external ran2
-      double precision x(99),p(0:3,nexternal)
       include 'nFKSconfigs.inc'
-      include 'reweight_all.inc'
-      integer nfksprocess_all
-      INTEGER NFKSPROCESS
-      COMMON/C_NFKSPROCESS/NFKSPROCESS
-      character*4 abrv
-      common /to_abrv/ abrv
-      logical nbody
-      common/cnbody/nbody
-      integer fks_j_from_i(nexternal,0:nexternal)
-     &     ,particle_type(nexternal),pdg_type(nexternal)
-      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
-      integer i_fks,j_fks
-      common/fks_indices/i_fks,j_fks
-      logical sum,firsttime
+      include 'c_weight.inc'
+      include 'reweight.inc'
+      include 'run.inc'
+      double precision xx(ndimmax),vegas_wgt,f(nintegrals),jac,p(0:3
+     $     ,nexternal),rwgt,vol,sig,x(99),MC_int_wgt
+      integer ifl,nFKS_born,nFKS_picked,iFKS,nFKS_min
+     $     ,nFKS_max,izero,ione,itwo,mohdr
+      parameter (izero=0,ione=1,itwo=2,mohdr=-100)
+      logical passcuts,passcuts_nbody,passcuts_n1body,sum
+      external passcuts
       parameter (sum=.false.)
+      integer         ndim,ipole
+      common/tosigint/ndim,ipole
+      logical       nbody
+      common/cnbody/nbody
+      integer           iconfig
+      common/to_configs/iconfig
+      double precision p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
+     $     ,pswgt_cnt(-2:2),jac_cnt(-2:2)
+      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
+      double precision p_born(0:3,nexternal-1)
+      common /pborn/   p_born
+      double precision           virt_wgt_mint,born_wgt_mint
+      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
+      double precision virtual_over_born
+      common/c_vob/virtual_over_born
+      logical                calculatedBorn
+      common/ccalculatedBorn/calculatedBorn
+      character*4      abrv
+      common /to_abrv/ abrv
+      integer iappl
+      common /for_applgrid/ iappl
+      if (ifl.ne.0) then
+         write (*,*) 'ERROR ifl not equal to zero in sigint',ifl
+         stop 1
+      endif
+      sigint=0d0
+      icontr=0
+      virt_wgt_mint=0d0
+      born_wgt_mint=0d0
+      virtual_over_born=0d0
+      if (ickkw.eq.3) call set_FxFx_scale(-1,p)
+      call update_vegas_x(xx,x)
+      call get_MC_integer(1,fks_configs,nFKS_picked,vol)
+
+c The nbody contributions
+      if (abrv.eq.'real') goto 11
+      nbody=.true.
+      calculatedBorn=.false.
+      call get_born_nFKSprocess(nFKS_picked,nFKS_born)
+      call update_fks_dir(nFKS_born,iconfig)
+      jac=1d0
+      call generate_momenta(ndim,iconfig,jac,x,p)
+      if (p_born(0,1).lt.0d0) goto 12
+      call compute_prefactors_nbody(vegas_wgt)
+      call set_cms_stuff(izero)
+      passcuts_nbody=passcuts(p1_cnt(0,1,0),rwgt)
+      if (passcuts_nbody) then
+         if (ickkw.eq.3) call set_FxFx_scale(izero,p1_cnt(0,1,0))
+         call set_alphaS(p1_cnt(0,1,0))
+         if (abrv(1:2).ne.'vi') then
+            call compute_born
+         endif
+         if (abrv.ne.'born') then
+            call compute_nbody_noborn
+         endif
+      endif
+
+ 11   continue
+c The n+1-body contributions (including counter terms)
+      if (abrv.eq.'born'.or.abrv(1:2).eq.'vi') goto 12
+      nbody=.false.
+      if (sum) then
+         nFKS_min=1
+         nFKS_max=fks_configs
+         MC_int_wgt=1d0
+      else
+         nFKS_min=nFKS_picked
+         nFKS_max=nFKS_picked
+         MC_int_wgt=1d0/vol
+      endif
+      do iFKS=nFKS_min,nFKS_max
+         jac=MC_int_wgt
+         call update_fks_dir(iFKS,iconfig)
+         call generate_momenta(ndim,iconfig,jac,x,p)
+         if (p_born(0,1).lt.0d0) cycle
+         call compute_prefactors_n1body(vegas_wgt,jac)
+         call set_cms_stuff(izero)
+         passcuts_nbody =passcuts(p1_cnt(0,1,0),rwgt)
+         call set_cms_stuff(mohdr)
+         passcuts_n1body=passcuts(p,rwgt)
+         if (passcuts_nbody .and. abrv.ne.'real') then
+            call set_cms_stuff(izero)
+            if (ickkw.eq.3) call set_FxFx_scale(izero,p1_cnt(0,1,0))
+            call set_alphaS(p1_cnt(0,1,0))
+            call compute_soft_counter_term
+            call set_cms_stuff(ione)
+            call compute_collinear_counter_term
+            call set_cms_stuff(itwo)
+            call compute_soft_collinear_counter_term
+         endif
+         if (passcuts_n1body) then
+            call set_cms_stuff(mohdr)
+            if (ickkw.eq.3) call set_FxFx_scale(mohdr,p)
+            call set_alphaS(p)
+            call compute_real_emission(p)
+         endif
+      enddo
+      
+ 12   continue
+c Include PDFs and alpha_S and reweight to include the uncertainties
+      call include_PDF_and_alphas
+      if (doreweight) then
+         if (do_rwgt_scale) call reweight_scale
+         if (do_rwgt_pdf) call reweight_pdf
+      endif
+      
+      if (iappl.ne.0) then
+         if (sum) then
+            write (*,*) 'ERROR: applgrid only possible '/
+     &           /'with MC over FKS directories',iappl,sum
+            stop 1
+         endif
+         call fill_applgrid_weights(vegas_wgt)
+      endif
+
+c Importance sampling for FKS configurations
+      if (sum) then
+         call get_wgt_nbody(sig)
+         call fill_MC_integer(1,nFKS_picked,abs(sig))
+      else
+         call get_wgt_no_nbody(sig)
+         call fill_MC_integer(1,nFKS_picked,abs(sig)*vol)
+      endif
+
+c Finalize PS point
+      call fill_plots
+      call fill_mint_function(f)
+      return
+      end
+
+      subroutine update_fks_dir(nFKS,iconfig)
+      implicit none
+      integer nFKS,iconfig
+      integer              nFKSprocess
+      common/c_nFKSprocess/nFKSprocess
+      nFKSprocess=nFKS
+      call fks_inc_chooser()
+      call leshouche_inc_chooser()
+      call setcuts
+      call setfksfactor(iconfig)
+      return
+      end
+      
+      subroutine get_born_nFKSprocess(nFKS_in,nFKS_out)
+      implicit none
+      include 'nexternal.inc'
+      include 'nFKSconfigs.inc'
+      include 'fks_info.inc'
+      integer nFKS_in,nFKS_out,iFKS,nFKSprocessBorn(2)
+      logical firsttime,foundB(2)
       data firsttime /.true./
-      logical foundB(2)
-      integer nFKSprocessBorn(2)
       save nFKSprocessBorn,foundB
-      double precision vol,sigintR
-      integer itotalpoints
-      common/ctotalpoints/itotalpoints
-      double precision virtual_over_born
-      common/c_vob/virtual_over_born
-      double precision virt_wgt_current,born_wgt_ao2pi_current
-      double precision virt_wgt,born_wgt_ao2pi
-      common/c_fks_singular/virt_wgt,born_wgt_ao2pi
-      logical fillh
-      integer mc_hel,ihel
-      double precision volh
-      common/mc_int2/volh,mc_hel,ihel,fillh
-c
-      if (ifl.ne.0) then
-         write (*,*) 'ifl not equal to zero in sigint()',ifl
-         stop
-      endif
-      virtual_over_born=0d0
-
-      do i=1,99
-         if (abrv.eq.'grid'.or.abrv.eq.'born'.or.abrv(1:2).eq.'vi')
-     &        then
-            if(i.le.ndim-3)then
-               x(i)=xx(i)
-            elseif(i.le.ndim) then
-               x(i)=ran2()      ! Choose them flat when not including real-emision
-            else
-               x(i)=0.d0
-            endif
-         else
-            if(i.le.ndim)then
-               x(i)=xx(i)
-            else
-               x(i)=0.d0
-            endif
-         endif
-      enddo
-
       if (firsttime) then
          firsttime=.false.
          foundB(1)=.false.
          foundB(2)=.false.
-         do nFKSprocess=1,fks_configs
-            call fks_inc_chooser()
-            if (particle_type(i_fks).eq.8) then
-               if (j_fks.le.nincoming) then
+         do iFKS=1,fks_configs
+            if (particle_type_D(iFKS,fks_i_D(iFKS)).eq.8) then
+               if (fks_j_D(iFKS).le.nincoming) then
                   foundB(1)=.true.
-                  nFKSprocessBorn(1)=nFKSprocess
+                  nFKSprocessBorn(1)=iFKS
                else
                   foundB(2)=.true.
-                  nFKSprocessBorn(2)=nFKSprocess
+                  nFKSprocessBorn(2)=iFKS
                endif
             endif
          enddo
          write (*,*) 'Total number of FKS directories is', fks_configs
          write (*,*) 'For the Born we use nFKSprocesses  #',
-     &        nFKSprocessBorn
+     $        nFKSprocessBorn
       endif
-
-      sigint=0d0
-c
-c Compute the Born-like contributions with nbody=.true.
-c THIS CAN BE OPTIMIZED
-c
-      call get_MC_integer(1,fks_configs,nFKSprocess,vol)
-      nFKSprocess_all=nFKSprocess
-      call fks_inc_chooser()
-      if (j_fks.le.nincoming) then
+      if (fks_j_D(nFKS_in).le.nincoming) then
          if (.not.foundB(1)) then
             write(*,*) 'Trying to generate Born momenta with '/
      &           /'initial state j_fks, but there is no '/
      &           /'configuration with i_fks a gluon and j_fks '/
      &           /'initial state'
-            stop
+            stop 1
          endif
-         nFKSprocess=nFKSprocessBorn(1)
+         nFKS_out=nFKSprocessBorn(1)
       else
          if (.not.foundB(2)) then
             write(*,*) 'Trying to generate Born momenta with '/
      &           /'final state j_fks, but there is no configuration'/
      &           /' with i_fks a gluon and j_fks final state'
-            stop
-         endif
-         nFKSprocess=nFKSprocessBorn(2)
-      endif
-      nbody=.true.
-      fillh=.false.  ! this is set to true in BinothLHA if doing MC over helicities
-      nFKSprocess_used=nFKSprocess
-      nFKSprocess_used_Born=nFKSprocess
-      call fks_inc_chooser()
-      call leshouche_inc_chooser()
-      call setcuts
-      call setfksfactor(iconfig)
-      wgt=1d0
-      call generate_momenta(ndim,iconfig,wgt,x,p)
-      sigint = sigint+dsig(p,wgt,peso)*peso
-      if (mc_hel.ne.0 .and. fillh) then
-c Fill the importance sampling array
-         call fill_MC_integer(2,ihel,abs(sigint*volh))
-      endif
-
-      virt_wgt_current=virt_wgt
-      born_wgt_ao2pi_current=born_wgt_ao2pi
-c
-c Compute the subtracted real-emission corrections either as an explicit
-c sum or a Monte Carlo sum.
-c      
-      if (abrv.ne.'born' .and. abrv.ne.'grid' .and.
-     &     abrv(1:2).ne.'vi') then
-         nbody=.false.
-         if (sum) then
-c THIS CAN BE OPTIMIZED
-            do nFKSprocess=1,fks_configs
-               nFKSprocess_used=nFKSprocess
-               call fks_inc_chooser()
-               call leshouche_inc_chooser()
-               call setcuts
-               call setfksfactor(iconfig)
-               wgt=1d0
-               call generate_momenta(ndim,iconfig,wgt,x,p)
-               sigint = sigint+dsig(p,wgt,peso)*peso
-            enddo
-         else                   ! Monte Carlo over nFKSprocess
-            nFKSprocess=nFKSprocess_all
-            nFKSprocess_used=nFKSprocess
-c THIS CAN BE OPTIMIZED
-            call fks_inc_chooser()
-            call leshouche_inc_chooser()
-            call setcuts
-            call setfksfactor(iconfig)
-c     The variable 'vol' is the size of the cell for the MC over
-c     nFKSprocess. Need to divide by it here to correctly take into
-c     account this Jacobian
-            wgt=1d0/vol
-            call generate_momenta(ndim,iconfig,wgt,x,p)
-            sigintR = dsig(p,wgt,peso)*peso
-            call fill_MC_integer(1,nFKSprocess,abs(sigintR)*vol)
-            sigint = sigint+ sigintR
-         endif
-      endif
-
-      f(1)=abs(sigint+virt_wgt_current)
-      f(2)=sigint+virt_wgt_current
-      f(3)=virt_wgt_current
-      f(4)=virtual_over_born
-      f(5)=abs(virt_wgt_current)
-      f(6)=born_wgt_ao2pi_current
-
-      if (sigint.ne.0d0)itotalpoints=itotalpoints+1
+            stop 1
+         endif
+         nFKS_out=nFKSprocessBorn(2)
+      endif
+      return
+      end
+
+      subroutine update_vegas_x(xx,x)
+      implicit none
+      include 'mint.inc'
+      integer i
+      double precision xx(ndimmax),x(99),ran2
+      external ran2
+      integer ndim,ipole
+      common/tosigint/ndim,ipole
+      character*4 abrv
+      common /to_abrv/ abrv
+      do i=1,99
+         if (abrv.eq.'born'.or.abrv(1:2).eq.'vi') then
+            if(i.le.ndim-3)then
+               x(i)=xx(i)
+            elseif(i.le.ndim) then
+               x(i)=ran2()      ! Choose them flat when not including real-emision
+            else
+               x(i)=0.d0
+            endif
+         else
+            if(i.le.ndim)then
+               x(i)=xx(i)
+            else
+               x(i)=0.d0
+            endif
+         endif
+      enddo
       return
       end
 
@@ -719,7 +788,7 @@
       abrv=abrvinput(1:4)
 c Options are way too many: make sure we understand all of them
       if ( abrv.ne.'all '.and.abrv.ne.'born'.and.abrv.ne.'real'.and.
-     &     abrv.ne.'virt'.and.abrv.ne.'novi'.and.abrv.ne.'grid'.and.
+     &     abrv.ne.'virt'.and.
      &     abrv.ne.'viSC'.and.abrv.ne.'viLC'.and.abrv.ne.'novA'.and.
      &     abrv.ne.'novB'.and.abrv.ne.'viSA'.and.abrv.ne.'viSB') then
         write(*,*)'Error in input: abrv is:',abrv

=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2014-08-05 13:04:11 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2015-02-04 13:27:21 +0000
@@ -1304,7 +1304,9 @@
       logical Hevents
       common/SHevents/Hevents
       integer i_process
-      common/c_addwrite/i_process
+      common/c_i_process/i_process
+      integer i_process_addwrite
+      common/c_addwrite/i_process_addwrite
       logical unwgt
       double precision evtsgn_save,evtsgn_target
       double precision evtsgn
@@ -1327,6 +1329,7 @@
       logical only_virt
       integer imode
       common /c_imode/imode,only_virt
+
 c Trivial check on the Born contribution
       do i=1,iproc_save(nFKSprocess_used_born)
          if (unwgt_table(0,2,i).ne.0d0) then
@@ -1568,6 +1571,7 @@
                      evtsgn=sign(1d0,unwgt_table(nFKSprocess,2
      $                    ,i_process))
                      evtsgn_save=evtsgn_save+evtsgn
+                     i_process_addwrite=i_process
                      if (doreweight) then
                         nFKSprocess_reweight(1)=nFKSprocess
                      endif
@@ -1593,9 +1597,9 @@
                      endif
                      evtsgn=sign(1d0,f_unwgt(nFKSprocess,i_process))
                      evtsgn_save=evtsgn_save+evtsgn
-c Set the i_process to one of the (n+1)-body configurations that leads
+c Set the i_process_addwrite to one of the (n+1)-body configurations that leads
 c to this Born configuration. Needed for add_write_info to work properly
-                     i_process=etoi(i_process,nFKSprocess)
+                     i_process_addwrite=etoi(i_process,nFKSprocess)
                      if (doreweight) then
 c for the reweight info, do not write the ones that gave a zero
 c contribution
@@ -1715,9 +1719,9 @@
                   evtsgn=sign(1d0,f_unwgt(nFKSprocess_used_born
      $                 ,i_process))
                   evtsgn_save=evtsgn_save + evtsgn
-c Set the i_process to one of the (n+1)-body configurations that leads
+c Set the i_process_addwrite to one of the (n+1)-body configurations that leads
 c to this Born configuration. Needed for add_write_info to work properly
-                  i_process=etoi(i_process,nFKSprocess_used_born)
+                  i_process_addwrite=etoi(i_process,nFKSprocess_used_born)
                   if (evtsgn.eq.evtsgn_target) return
                enddo
 c Pick the sign randomly according to all the signs accumulated. 

=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2014-10-17 19:39:20 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2015-02-04 13:27:21 +0000
@@ -1,3 +1,1230 @@
+      subroutine compute_born
+c This subroutine computes the Born matrix elements and adds its value
+c to the list of weights using the add_wgt subroutine
+      implicit none
+      include 'nexternal.inc'
+      include 'reweight0.inc'
+      include 'coupl.inc'
+      include 'timing_variables.inc'
+      double complex wgt_c(2)
+      double precision wgt1
+      double precision p_born(0:3,nexternal-1)
+      common /pborn/   p_born
+      double precision      f_b,f_nb
+      common /factor_nbody/ f_b,f_nb
+      double precision    xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev(0:3)
+     $                    ,p_i_fks_cnt(0:3,-2:2)
+      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
+      double precision     xiScut_used,xiBSVcut_used
+      common /cxiScut_used/xiScut_used,xiBSVcut_used
+      call cpu_time(tBefore)
+      if (f_b.eq.0d0) return
+      if (xi_i_fks_ev .gt. xiBSVcut_used) return
+      call sborn(p_born,wgt_c)
+      wgt1=dble(wgt_c(1))*f_b/g**(nint(2*wgtbpower))
+      call add_wgt(2,wgt1,0d0,0d0)
+      call cpu_time(tAfter)
+      tBorn=tBorn+(tAfter-tBefore)
+      return
+      end
+
+      subroutine compute_nbody_noborn
+c This subroutine computes the soft-virtual matrix elements and adds its
+c value to the list of weights using the add_wgt subroutine
+      implicit none
+      include 'nexternal.inc'
+      include 'reweight.inc'
+      include 'coupl.inc'
+      include 'run.inc'
+      include 'timing_variables.inc'
+      double precision wgt1,wgt2,wgt3,bsv_wgt,virt_wgt,born_wgt,pi,g2,g22
+      parameter (pi=3.1415926535897932385d0)
+      double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
+     $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
+      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
+      double precision           virt_wgt_mint,born_wgt_mint
+      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
+      double precision      f_b,f_nb
+      common /factor_nbody/ f_b,f_nb
+      double precision    xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev(0:3)
+     $                    ,p_i_fks_cnt(0:3,-2:2)
+      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
+      double precision     xiScut_used,xiBSVcut_used
+      common /cxiScut_used/xiScut_used,xiBSVcut_used
+      double precision fxfx_exp_rewgt
+      common /c_fxfx_exp_regt/ fxfx_exp_rewgt
+      call cpu_time(tBefore)
+      if (f_nb.eq.0d0) return
+      if (xi_i_fks_ev .gt. xiBSVcut_used) return
+      call bornsoftvirtual(p1_cnt(0,1,0),bsv_wgt,virt_wgt,born_wgt)
+      g2=g**(nint(2*wgtbpower))
+      g22=g**(nint(2*wgtbpower+2))
+      wgt1=wgtnstmp*f_nb/g22
+      if (ickkw.eq.3 .and. fxfx_exp_rewgt.ne.0d0) then
+         wgt1=wgt1 - fxfx_exp_rewgt*born_wgt*f_nb/g2/(4d0*pi)
+      endif
+      wgt2=wgtwnstmpmur*f_nb/g22
+      wgt3=wgtwnstmpmuf*f_nb/g22
+      call add_wgt(3,wgt1,wgt2,wgt3)
+c Special for the soft-virtual needed for the virt-tricks. The
+c *_wgt_mint variable should be directly passed to the mint-integrator
+c and not be part of the plots nor computation of the cross section.
+      virt_wgt_mint=virt_wgt*f_nb/g22
+      born_wgt_mint=born_wgt*f_b/g2
+      call cpu_time(tAfter)
+      tIS=tIS+(tAfter-tBefore)
+      return
+      end
+
+      subroutine compute_real_emission(p)
+c This subroutine computes the real-emission matrix elements and adds
+c its value to the list of weights using the add_wgt subroutine
+      implicit none
+      include 'nexternal.inc'
+      include 'coupl.inc'
+      include 'reweight0.inc'
+      include 'timing_variables.inc'
+      double precision x,dot,f_damp,ffact,s_ev,fks_Sij,p(0:3,nexternal)
+     $     ,wgt1,fx_ev
+      external dot,f_damp,fks_Sij
+      double precision        ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
+      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
+      integer            i_fks,j_fks
+      common/fks_indices/i_fks,j_fks
+      double precision    xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev(0:3)
+     $                    ,p_i_fks_cnt(0:3,-2:2)
+      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
+      double precision     f_r,f_s,f_c,f_dc,f_sc,f_dsc(4)
+      common/factor_n1body/f_r,f_s,f_c,f_dc,f_sc,f_dsc
+      call cpu_time(tBefore)
+      if (f_r.eq.0d0) return
+      x = abs(2d0*dot(p(0,i_fks),p(0,j_fks))/shat)
+      ffact = f_damp(x)
+      if (ffact.le.0d0) return
+      s_ev = fks_Sij(p,i_fks,j_fks,xi_i_fks_ev,y_ij_fks_ev)
+      if (s_ev.le.0.d0) return
+      call sreal(p,xi_i_fks_ev,y_ij_fks_ev,fx_ev)
+      wgt1=fx_ev*s_ev*f_r/g**(nint(2*wgtbpower+2))
+      call add_wgt(1,wgt1,0d0,0d0)
+      call cpu_time(tAfter)
+      tReal=tReal+(tAfter-tBefore)
+      return
+      end
+
+      subroutine compute_soft_counter_term
+c This subroutine computes the soft counter term and adds its value to
+c the list of weights using the add_wgt subroutine
+      implicit none
+      include 'nexternal.inc'
+      include 'coupl.inc'
+      include 'reweight0.inc'
+      include 'timing_variables.inc'
+      double precision wgt1,s_s,fks_Sij,fx_s,zero
+      parameter (zero=0d0)
+      external fks_Sij
+      double precision     p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
+     $                     ,pswgt_cnt(-2:2),jac_cnt(-2:2)
+      common/counterevnts/ p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
+      double precision     xiScut_used,xiBSVcut_used
+      common /cxiScut_used/xiScut_used,xiBSVcut_used
+      integer            i_fks,j_fks
+      common/fks_indices/i_fks,j_fks
+      double precision    xi_i_fks_ev,y_ij_fks_ev
+      double precision    p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
+      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
+      double precision     f_r,f_s,f_c,f_dc,f_sc,f_dsc(4)
+      common/factor_n1body/f_r,f_s,f_c,f_dc,f_sc,f_dsc
+      call cpu_time(tBefore)
+      if (f_s.eq.0d0) return
+      if (xi_i_fks_ev .gt. xiScut_used) return
+      s_s = fks_Sij(p1_cnt(0,1,0),i_fks,j_fks,zero,y_ij_fks_ev)
+      if (s_s.le.0d0) return
+      call sreal(p1_cnt(0,1,0),0d0,y_ij_fks_ev,fx_s)
+      wgt1=-fx_s*s_s*f_s/g**(nint(2*wgtbpower+2))
+      call add_wgt(4,wgt1,0d0,0d0)
+      call cpu_time(tAfter)
+      tCount=tCount+(tAfter-tBefore)
+      return
+      end
+
+      subroutine compute_collinear_counter_term
+c This subroutine computes the collinear counter term and adds its value
+c to the list of weights using the add_wgt subroutine
+      implicit none
+      include 'nexternal.inc'
+      include 'coupl.inc'
+      include 'fks_powers.inc'
+      include 'reweight.inc'
+      include 'timing_variables.inc'
+      double precision zero,one,s_c,fks_Sij,fx_c,deg_xi_c,deg_lxi_c,wgt1
+     &     ,wgt3,g22
+      external fks_Sij
+      parameter (zero=0d0,one=1d0)
+      double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
+     $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
+      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
+      integer            i_fks,j_fks
+      common/fks_indices/i_fks,j_fks
+      double precision    xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev(0:3)
+     $                    ,p_i_fks_cnt(0:3,-2:2)
+      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
+      double precision   xi_i_fks_cnt(-2:2)
+      common /cxiifkscnt/xi_i_fks_cnt
+      double precision     f_r,f_s,f_c,f_dc,f_sc,f_dsc(4)
+      common/factor_n1body/f_r,f_s,f_c,f_dc,f_sc,f_dsc
+      double precision pmass(nexternal)
+      call cpu_time(tBefore)
+      include 'pmass.inc'
+      if (f_c.eq.0d0 .and. f_dc.eq.0d0)return
+      if (y_ij_fks_ev.le.1d0-deltaS .or. pmass(j_fks).ne.0.d0) return
+      s_c = fks_Sij(p1_cnt(0,1,1),i_fks,j_fks,xi_i_fks_cnt(1),one)
+      if (s_c.le.0d0) return
+      g22=g**(nint(2*wgtbpower+2))
+      call sreal(p1_cnt(0,1,1),xi_i_fks_cnt(1),one,fx_c)
+      wgt1=-fx_c*s_c*f_c/g22
+      call sreal_deg(p1_cnt(0,1,1),xi_i_fks_cnt(1),one,deg_xi_c
+     $     ,deg_lxi_c)
+      wgt1=wgt1+ ( wgtdegrem_xi+wgtdegrem_lxi*log(xi_i_fks_cnt(1)) )*
+     $     f_dc/g22
+      wgt3=wgtdegrem_muF*f_dc/g22
+      call add_wgt(5,wgt1,0d0,wgt3)
+      call cpu_time(tAfter)
+      tCount=tCount+(tAfter-tBefore)
+      return
+      end
+
+      subroutine compute_soft_collinear_counter_term
+c This subroutine computes the soft-collinear counter term and adds its
+c value to the list of weights using the add_wgt subroutine
+      implicit none
+      include 'nexternal.inc'
+      include 'coupl.inc'
+      include 'reweight.inc'
+      include 'fks_powers.inc'
+      include 'timing_variables.inc'
+      double precision zero,one,s_sc,fks_Sij,fx_sc,wgt1,wgt3,deg_xi_sc
+     $     ,deg_lxi_sc,g22
+      external fks_Sij
+      parameter (zero=0d0,one=1d0)
+      double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
+     $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
+      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
+      double precision    xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev(0:3)
+     $                    ,p_i_fks_cnt(0:3,-2:2)
+      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
+      integer            i_fks,j_fks
+      common/fks_indices/i_fks,j_fks
+      double precision     xiScut_used,xiBSVcut_used
+      common /cxiScut_used/xiScut_used,xiBSVcut_used
+      double precision   xi_i_fks_cnt(-2:2)
+      common /cxiifkscnt/xi_i_fks_cnt
+      double precision     f_r,f_s,f_c,f_dc,f_sc,f_dsc(4)
+      common/factor_n1body/f_r,f_s,f_c,f_dc,f_sc,f_dsc
+      double precision pmass(nexternal)
+      include 'pmass.inc'
+      call cpu_time(tBefore)
+      if (f_sc.eq.0d0 .and. f_dsc(1).eq.0d0 .and. f_dsc(2).eq.0d0 .and.
+     $     f_dsc(3).eq.0d0 .and. f_dsc(4).eq.0d0) return
+      if (xi_i_fks_cnt(1).ge.xiScut_used .or. y_ij_fks_ev.le.1d0-deltaS
+     $     .or. pmass(j_fks).ne.0.d0 ) return
+      s_sc = fks_Sij(p1_cnt(0,1,2),i_fks,j_fks,zero,one)
+      if (s_sc.le.0d0) return
+      g22=g**(nint(2*wgtbpower+2))
+      call sreal(p1_cnt(0,1,2),zero,one,fx_sc)
+      wgt1=fx_sc*s_sc*f_sc/g22
+      call sreal_deg(p1_cnt(0,1,2),zero,one, deg_xi_sc,deg_lxi_sc)
+      wgt1=wgt1+(-(wgtdegrem_xi+wgtdegrem_lxi*log(xi_i_fks_cnt(1)))
+     &     *f_dsc(1)-(wgtdegrem_xi*f_dsc(2)+wgtdegrem_lxi*f_dsc(3)))/g22
+      wgt3=-wgtdegrem_muF*f_dsc(4)/g22
+      call add_wgt(6,wgt1,0d0,wgt3)
+      call cpu_time(tAfter)
+      tCount=tCount+(tAfter-tBefore)
+      return
+      end
+
+      logical function pdg_equal(pdg1,pdg2)
+c Returns .true. if the lists of PDG codes --'pdg1' and 'pdg2'-- are
+c equal.
+      implicit none
+      include 'nexternal.inc'
+      integer i,pdg1(nexternal),pdg2(nexternal)
+      pdg_equal=.true.
+      do i=1,nexternal
+         if (pdg1(i).ne.pdg2(i)) then
+            pdg_equal=.false.
+            return
+         endif
+      enddo
+      end
+      
+      logical function momenta_equal(p1,p2)
+c Returns .true. only if the momenta p1 and p2 are equal. To save time,
+c it only checks the 0th and 3rd components (energy and z-direction).
+      implicit none
+      include 'nexternal.inc'
+      integer i,j
+      double precision p1(0:3,nexternal),p2(0:3,nexternal),vtiny
+      parameter (vtiny=1d-8)
+      momenta_equal=.true.
+      do i=1,nexternal
+         do j=0,3,3
+            if (p1(j,i).eq.0d0 .or. p2(j,i).eq.0d0) then
+               if (abs(p1(j,i)-p2(j,i)).gt.vtiny) then
+                  momenta_equal=.false.
+                  return
+               endif
+            else
+               if (abs((p1(j,i)-p2(j,i))/max(p1(j,i),p2(j,i))).gt.vtiny)
+     &              then
+                  momenta_equal=.false.
+                  return
+               endif
+            endif
+         enddo
+      enddo
+      end
+      
+      subroutine set_FxFx_scale(iterm,p)
+c Sets the FxFx cluster scale and multiplies the f_* factors (computed
+c by 'compute_prefactors_nbody' and 'compute_prefactors_n1body') by the
+c Sudakov suppression. If called more than once with the same momenta
+c and iterm, skip setting of the scales, and multiply the f_* factors by
+c the cached Sudakovs.
+c     iterm= -1: reset the computation of the Sudakovs
+c     iterm=  0: Sudakov for n-body kinematics
+c     iterm=100: Sudakov for n+1-body kinematics
+      implicit none
+      include 'nexternal.inc'
+      include 'run.inc'
+      include 'timing_variables.inc'
+      integer iterm,iterm_last,i,j
+      double precision p(0:3,nexternal),p_last(0:3,nexternal),rewgt
+     &     ,rewgt_izero,rewgt_mohdr,rewgt_exp_izero,rewgt_exp_mohdr
+      logical setclscales,rewgt_izero_calculated,rewgt_mohdr_calculated
+     &     ,momenta_equal,already_set
+      external setclscales,rewgt,momenta_equal
+      double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
+     $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
+      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
+      double precision     f_r,f_s,f_c,f_dc,f_sc,f_dsc(4)
+      common/factor_n1body/f_r,f_s,f_c,f_dc,f_sc,f_dsc
+      double precision      f_b,f_nb
+      common /factor_nbody/ f_b,f_nb
+      double precision fxfx_exp_rewgt
+      common /c_fxfx_exp_regt/ fxfx_exp_rewgt
+      call cpu_time(tBefore)
+      ktscheme=1
+      if (iterm.eq.-1) then
+         rewgt_mohdr_calculated=.false.
+         rewgt_izero_calculated=.false.
+         fxfx_exp_rewgt=0d0
+         return
+      endif
+      already_set=.false.
+      if (iterm.eq.0) then
+         if (rewgt_izero_calculated) then
+            if (iterm.eq.iterm_last) then
+               if (momenta_equal(p1_cnt(0,1,0),p_last)) then
+                  already_set=.true.
+               endif
+            endif
+         endif
+         if (.not.already_set) then
+            if (.not. setclscales(p1_cnt(0,1,0))) then
+               write (*,*) 'ERROR in setclscales izero'
+               stop 1
+            endif
+            rewgt_izero=rewgt(p1_cnt(0,1,0),rewgt_exp_izero)
+            fxfx_exp_rewgt=rewgt_exp_izero
+         endif
+         rewgt_izero_calculated=.true.
+         iterm_last=iterm
+         do i=1,nexternal
+            do j=0,3
+               p_last(j,i)=p1_cnt(j,i,0)
+            enddo
+         enddo
+         f_b =f_b *rewgt_izero
+         f_nb=f_nb*rewgt_izero
+         f_s =f_s *rewgt_izero
+         f_c =f_c *rewgt_izero
+         f_dc=f_dc*rewgt_izero
+         f_sc=f_sc*rewgt_izero
+         do i=1,4
+            f_dsc(i)=f_dsc(i)*rewgt_izero
+         enddo
+         call cpu_time(tAfter)
+         tFxFx=tFxFx+(tAfter-tBefore)
+         return
+      endif
+      if (iterm.eq.-100) then
+         if (rewgt_mohdr_calculated) then
+            if (iterm.eq.iterm_last) then
+               if (momenta_equal(p,p_last)) then
+                  already_set=.true.
+               endif
+            endif
+         endif
+         if (.not. already_set) then
+            if (.not. setclscales(p)) then
+               write (*,*) 'ERROR in setclscales mohdr'
+               stop 1
+            endif
+            rewgt_mohdr=rewgt(p,rewgt_exp_mohdr)
+         endif
+         rewgt_mohdr_calculated=.true.
+         iterm_last=iterm
+         do i=1,nexternal
+            do j=0,3
+               p_last(j,i)=p(j,i)
+            enddo
+         enddo
+         f_r=f_r*rewgt_mohdr
+         call cpu_time(tAfter)
+         tFxFx=tFxFx+(tAfter-tBefore)
+         return
+      endif
+      end
+      
+      
+      subroutine compute_prefactors_nbody(vegas_wgt)
+c Compute all the relevant prefactors for the Born and the soft-virtual,
+c i.e. all the nbody contributions. Also initialises the plots and
+c bpower.
+      implicit none
+      include 'nexternal.inc'
+      include 'run.inc'
+      include 'genps.inc'
+      include 'reweight0.inc'
+      include 'timing_variables.inc'
+      double precision pi,unwgtfun,vegas_wgt,enhance,xnoborn_cnt,xtot
+     $     ,bpower,cpower,tiny
+      data xnoborn_cnt /0d0/
+      integer inoborn_cnt,i
+      data inoborn_cnt /0/
+      double complex wgt_c(2)
+      logical firsttime
+      data firsttime /.true./
+      parameter (pi=3.1415926535897932385d0)
+      parameter (tiny=1d-6)
+      double precision p_born(0:3,nexternal-1)
+      common/pborn/    p_born
+      double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
+     $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
+      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
+      double precision  xinorm_ev
+      common /cxinormev/xinorm_ev
+      double precision  xiimax_ev
+      common /cxiimaxev/xiimax_ev
+      double precision     xiScut_used,xiBSVcut_used
+      common /cxiScut_used/xiScut_used,xiBSVcut_used
+      double precision        ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
+      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
+      double precision         fkssymmetryfactor,fkssymmetryfactorBorn,
+     $                         fkssymmetryfactorDeg
+      integer                  ngluons,nquarks(-6:6)
+      common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
+     &                         fkssymmetryfactorDeg,ngluons,nquarks
+      integer            mapconfig(0:lmaxconfigs), iconfig
+      common/to_mconfigs/mapconfig,                iconfig
+      Double Precision amp2(maxamps), jamp2(0:maxamps)
+      common/to_amps/  amp2,          jamp2
+      double precision   diagramsymmetryfactor
+      common /dsymfactor/diagramsymmetryfactor
+      double precision      f_b,f_nb
+      common /factor_nbody/ f_b,f_nb
+      integer iappl
+      common /for_applgrid/ iappl
+      logical needrndec
+      parameter (needrndec=.true.)
+      real*8 ran2
+      external ran2
+      real*8 rndec(10)
+      common/crndec/rndec
+      include "appl_common.inc" 
+      call cpu_time(tBefore)
+c Random numbers to be used in the plotting routine
+      if(needrndec)then
+         do i=1,10
+            rndec(i)=ran2()
+         enddo
+      endif
+      if (firsttime) then
+c Put here call to compute bpower
+         call compute_bpower(p_born,bpower)
+         wgtbpower=bpower
+c Store the power of alphas of the Born events in the appl common block.
+         if(iappl.ne.0) appl_bpower = wgtbpower
+c Initialize hiostograms
+         call initplot
+c Compute cpower done for bottom Yukawa, routine needs to be adopted
+c for other muR-dependendent factors
+         call compute_cpower(p_born,cpower)
+         if(dabs(cpower+1d0).lt.tiny) then
+            wgtcpower=0d0
+         else
+            wgtcpower=cpower
+         endif
+c Check that things are done consistently
+         if(wgtcpower.ne.cpowerinput.and.dabs(cpower+1d0).gt.tiny)then
+           write(*,*)'Inconsistency in the computation of cpower',
+     #               wgtcpower,cpowerinput
+           write(*,*)'Check value in reweight0.inc'
+           stop
+         endif
+         firsttime=.false.
+      endif
+c Compute the multi-channel enhancement factor 'enhance'.
+      enhance=1.d0
+      if (p_born(0,1).gt.0d0) then
+         call sborn(p_born,wgt_c)
+      elseif(p_born(0,1).lt.0d0)then
+         enhance=0d0
+      endif
+      if (enhance.eq.0d0)then
+         xnoborn_cnt=xnoborn_cnt+1.d0
+         if(log10(xnoborn_cnt).gt.inoborn_cnt)then
+            write (*,*) 'WARNING: no Born momenta more than 10**',
+     $           inoborn_cnt,'times'
+            inoborn_cnt=inoborn_cnt+1
+         endif
+      else
+         xtot=0d0
+         if (mapconfig(0).eq.0) then
+            write (*,*) 'Fatal error in compute_prefactor_nbody:'/
+     &           /' no Born diagrams ',mapconfig,
+     &           '. Check bornfromreal.inc'
+            write (*,*) 'Is fks_singular compiled correctly?'
+            stop 1
+         endif
+         do i=1, mapconfig(0)
+            xtot=xtot+amp2(mapconfig(i))
+         enddo
+         if (xtot.ne.0d0) then
+            enhance=amp2(mapconfig(iconfig))/xtot
+            enhance=enhance*diagramsymmetryfactor
+         else
+            enhance=0d0
+         endif
+      endif
+      call unweight_function(p_born,unwgtfun)
+      call set_cms_stuff(0)
+c f_* multiplication factors for Born and nbody
+      f_b=jac_cnt(0)*xinorm_ev/(min(xiimax_ev,xiBSVcut_used)*shat/(16
+     $     *pi**2))*enhance*unwgtfun *fkssymmetryfactorBorn*vegas_wgt
+      f_nb=f_b
+      call cpu_time(tAfter)
+      tf_nb=tf_nb+(tAfter-tBefore)
+      return
+      end
+
+      subroutine compute_prefactors_n1body(vegas_wgt,jac_ev)
+c Compute all relevant prefactors for the real emission and counter
+c terms.
+      implicit none
+      include 'nexternal.inc'
+      include 'run.inc'
+      include 'genps.inc'
+      include 'fks_powers.inc'
+      include 'coupl.inc'
+      include 'timing_variables.inc'
+      double precision unwgtfun,vegas_wgt,enhance,xnoborn_cnt,xtot
+     &     ,prefact,prefact_cnt_ssc,prefact_deg,prefact_c,prefact_coll
+     &     ,jac_ev,pi,prefact_cnt_ssc_c,prefact_coll_c,prefact_deg_slxi
+     &     ,prefact_deg_sxi,zero
+      parameter (pi=3.1415926535897932385d0, zero=0d0)
+      data xnoborn_cnt /0d0/
+      integer inoborn_cnt,i
+      data inoborn_cnt /0/
+      double complex wgt_c(2)
+      double precision p_born(0:3,nexternal-1)
+      common/pborn/    p_born
+      double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
+     $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
+      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
+      double precision    xi_i_fks_ev,y_ij_fks_ev
+      double precision    p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
+      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
+      integer            i_fks,j_fks
+      common/fks_indices/i_fks,j_fks
+      double precision   xi_i_fks_cnt(-2:2)
+      common /cxiifkscnt/xi_i_fks_cnt
+      double precision  xinorm_ev
+      common /cxinormev/xinorm_ev
+      double precision  xiimax_ev
+      common /cxiimaxev/xiimax_ev
+      double precision   xiimax_cnt(-2:2)
+      common /cxiimaxcnt/xiimax_cnt
+      double precision   xinorm_cnt(-2:2)
+      common /cxinormcnt/xinorm_cnt
+      double precision    delta_used
+      common /cdelta_used/delta_used
+      double precision    xicut_used
+      common /cxicut_used/xicut_used
+      double precision     xiScut_used,xiBSVcut_used
+      common /cxiScut_used/xiScut_used,xiBSVcut_used
+      double precision        ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
+      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
+      double precision         fkssymmetryfactor,fkssymmetryfactorBorn,
+     &                         fkssymmetryfactorDeg
+      integer                  ngluons,nquarks(-6:6)
+      common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
+     &                         fkssymmetryfactorDeg,ngluons,nquarks
+      integer            mapconfig(0:lmaxconfigs), iconfig
+      common/to_mconfigs/mapconfig,                iconfig
+      Double Precision amp2(maxamps), jamp2(0:maxamps)
+      common/to_amps/  amp2,          jamp2
+      double precision   diagramsymmetryfactor
+      common /dsymfactor/diagramsymmetryfactor
+      logical nocntevents
+      common/cnocntevents/nocntevents
+      double precision     f_r,f_s,f_c,f_dc,f_sc,f_dsc(4)
+      common/factor_n1body/f_r,f_s,f_c,f_dc,f_sc,f_dsc
+      double precision pmass(nexternal)
+      include 'pmass.inc'
+      call cpu_time(tBefore)
+      enhance=1.d0
+      if (p_born(0,1).gt.0d0) then
+         call sborn(p_born,wgt_c)
+      elseif(p_born(0,1).lt.0d0)then
+         enhance=0d0
+      endif
+c Compute the multi-channel enhancement factor 'enhance'.
+      if (enhance.eq.0d0)then
+         xnoborn_cnt=xnoborn_cnt+1.d0
+         if(log10(xnoborn_cnt).gt.inoborn_cnt)then
+            write (*,*) 'WARNING: no Born momenta more than 10**',
+     $           inoborn_cnt,'times'
+            inoborn_cnt=inoborn_cnt+1
+         endif
+      else
+         xtot=0d0
+         if (mapconfig(0).eq.0) then
+            write (*,*) 'Fatal error in compute_prefactor_n1body,'/
+     &           /' no Born diagrams ',mapconfig
+     &           ,'. Check bornfromreal.inc'
+            write (*,*) 'Is fks_singular compiled correctly?'
+            stop 1
+         endif
+         do i=1, mapconfig(0)
+            xtot=xtot+amp2(mapconfig(i))
+         enddo
+         if (xtot.ne.0d0) then
+            enhance=amp2(mapconfig(iconfig))/xtot
+            enhance=enhance*diagramsymmetryfactor
+         else
+            enhance=0d0
+         endif
+      endif
+      call unweight_function(p_born,unwgtfun)
+      prefact=xinorm_ev/xi_i_fks_ev/(1-y_ij_fks_ev)
+
+c f_* multiplication factors for real-emission, soft counter, ... etc.       
+      f_r=prefact*jac_ev*enhance*unwgtfun*fkssymmetryfactor*vegas_wgt
+      if (.not.nocntevents) then
+         prefact_cnt_ssc=xinorm_ev/min(xiimax_ev,xiScut_used)*
+     &        log(xicut_used/min(xiimax_ev,xiScut_used))/(1
+     &        -y_ij_fks_ev)
+         f_s=(prefact+prefact_cnt_ssc)*jac_cnt(0)*enhance
+     $        *unwgtfun*fkssymmetryfactor*vegas_wgt
+         
+         if (pmass(j_fks).eq.0d0) then
+            prefact_c=xinorm_cnt(1)/xi_i_fks_cnt(1)/(1-y_ij_fks_ev)
+            prefact_coll=xinorm_cnt(1)/xi_i_fks_cnt(1)*log(delta_used
+     $           /deltaS)/deltaS
+            f_c=(prefact_c+prefact_coll)*jac_cnt(1)
+     $           *enhance*unwgtfun*fkssymmetryfactor*vegas_wgt
+
+            call set_cms_stuff(1)
+            prefact_deg=xinorm_cnt(1)/xi_i_fks_cnt(1)/deltaS
+            prefact_cnt_ssc_c=xinorm_cnt(1)/min(xiimax_cnt(1)
+     &           ,xiScut_used)*log(xicut_used/min(xiimax_cnt(1)
+     &           ,xiScut_used))*1/(1-y_ij_fks_ev)
+            prefact_coll_c=xinorm_cnt(1)/min(xiimax_cnt(1),xiScut_used)
+     $           *log(xicut_used/min(xiimax_cnt(1),xiScut_used))
+     $           *log(delta_used/deltaS)/deltaS
+            f_dc=jac_cnt(1)*prefact_deg/(shat/(32*pi**2))*enhance
+     $           *unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+            f_sc=(prefact_c+prefact_coll+prefact_cnt_ssc_c
+     &           +prefact_coll_c)*jac_cnt(2)*enhance*unwgtfun
+     &           *fkssymmetryfactorDeg*vegas_wgt
+
+            call set_cms_stuff(2)
+            prefact_deg_sxi=xinorm_cnt(1)/min(xiimax_cnt(1),xiScut_used)
+     &           *log(xicut_used/min(xiimax_cnt(1),xiScut_used))*1
+     &           /deltaS
+            prefact_deg_slxi=xinorm_cnt(1)/min(xiimax_cnt(1)
+     &           ,xiScut_used)*( log(xicut_used)**2
+     &           -log(min(xiimax_cnt(1),xiScut_used))**2 )*1/(2.d0
+     &           *deltaS)
+            f_dsc(1)=prefact_deg*jac_cnt(2)/(shat/(32*pi**2))*enhance
+     &           *unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+            f_dsc(2)=prefact_deg_sxi*jac_cnt(2)/(shat/(32*pi**2))
+     &           *enhance*unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+            f_dsc(3)=prefact_deg_slxi*jac_cnt(2)/(shat/(32*pi**2))
+     &           *enhance*unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+            f_dsc(4)=( prefact_deg+prefact_deg_sxi )*jac_cnt(2)/(shat
+     &           /(32*pi**2))*enhance*unwgtfun*fkssymmetryfactorDeg
+     &           *vegas_wgt
+         else
+            f_c=0d0
+            f_dc=0d0
+            f_sc=0d0
+            do i=1,4
+               f_dsc(i)=0d0
+            enddo
+         endif
+      else
+         f_s=0d0
+         f_c=0d0
+         f_dc=0d0
+         f_sc=0d0
+         do i=1,4
+            f_dsc(i)=0d0
+         enddo
+      endif
+      call cpu_time(tAfter)
+      tf_all=tf_all+(tAfter-tBefore)
+      return
+      end
+
+      
+      subroutine add_wgt(type,wgt1,wgt2,wgt3)
+c Adds a contribution to the list in c_weight.inc. 'type' sets the type
+c of the contribution and wgt1..wgt3 are the coefficients multiplying
+c the logs. The arguments are:
+c     type=1 : real-emission
+c     type=2 : Born
+c     type=3 : soft-virtual
+c     type=4 : soft counter-term
+c     type=5 : collinear counter-term
+c     type=6 : soft-collinear counter-term
+c     wgt1 : weight of the contribution not multiplying a scale log
+c     wgt2 : coefficient of the weight multiplying the log[mu_R^2/Q^2]
+c     wgt3 : coefficient of the weight multiplying the log[mu_F^2/Q^2]
+c
+c This subroutine increments the 'icontr' counter: each new call to this
+c function makes sure that it's considered a new contribution. For each
+c contribution, we save the
+c     The type: itype(icontr)
+c     The weights: wgt(1,icontr),wgt(2,icontr) and wgt(3,icontr) for
+c         wgt1, wgt2 and wgt3, respectively.
+c     The Bjorken x's: bjx(1,icontr), bjx(2,icontr)
+c     The Ellis-Sexton scale squared used to compute the weight:
+c        scales2(1,icontr)
+c     The renormalisation scale squared used to compute the weight:
+c        scales2(2,icontr)
+c     The factorisation scale squared used to compute the weight:
+c       scales2(3,icontr)
+c     The value of the strong coupling: g_strong(icontr)
+c     The FKS configuration: nFKS(icontr)
+c     The boost to go from the momenta in the C.o.M. frame to the
+c         laboratory frame: y_bst(icontr)      
+c     The power of the strong coupling (g_strong) for the current
+c       weight: QCDpower(icontr)
+c     The momenta: momenta(j,i,icontr). For the Born contribution, the
+c        counter-term momenta are used. This is okay for any IR-safe
+c        observables.
+c     The PDG codes: pdg(i,icontr). Always the ones with length
+c        'nexternal' are used, because the momenta are also the 
+c        'nexternal' ones. This is okay for IR-safe observables.
+c
+c Not set in this subroutine, but included in the c_weights common block
+c are the
+c     wgts(iwgt,icontr) : weights including scale/PDFs/logs. These are
+c        normalised so that they can be used directly to compute cross
+c        sections and fill plots. 'iwgt' goes from 1 to the maximum
+c        number of weights obtained from scale and PDF reweighting, with
+c        the iwgt=1 element being the central value.
+c     plot_id(icontr) : =20 for Born, 11 for real-emission and 12 for
+c        anything else.
+c     plot_wgts(iwgt,icontr) : same as wgts(), but only non-zero for
+c        unique contributions and non-unique are added to the unique
+c        ones. 'Unique' here is defined that they would be identical in
+c        an analysis routine (i.e. same momenta and PDG codes)
+      implicit none
+      include 'nexternal.inc'
+      include 'run.inc'
+      include 'genps.inc'
+      include 'coupl.inc'
+      include 'fks_info.inc'
+      include 'c_weight.inc'
+      include 'q_es.inc'
+      include 'reweight0.inc'
+      integer type,i,j
+      double precision wgt1,wgt2,wgt3
+      integer              nFKSprocess
+      common/c_nFKSprocess/nFKSprocess
+      double precision p_born(0:3,nexternal-1)
+      common/pborn/    p_born
+      double precision p_ev(0:3,nexternal)
+      common/pev/      p_ev
+      integer    maxflow
+      parameter (maxflow=999)
+      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
+     $     icolup(2,nexternal,maxflow)
+      common /c_leshouche_inc/idup,mothup,icolup
+      double precision        ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
+      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
+      double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
+     $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
+      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
+      if (wgt1.eq.0d0 .and. wgt2.eq.0d0 .and. wgt3.eq.0d0) return
+      icontr=icontr+1
+      if (icontr.gt.max_contr) then
+         write (*,*) 'ERROR in add_wgt: too many contributions'
+     &        ,max_contr
+         stop 1
+      endif
+      itype(icontr)=type
+      wgt(1,icontr)=wgt1
+      wgt(2,icontr)=wgt2
+      wgt(3,icontr)=wgt3
+      bjx(1,icontr)=xbk(1)
+      bjx(2,icontr)=xbk(2)
+      scales2(1,icontr)=QES2
+      scales2(2,icontr)=scale**2
+      scales2(3,icontr)=q2fact(1)
+      g_strong(icontr)=g
+      nFKS(icontr)=nFKSprocess
+      y_bst(icontr)=ybst_til_tolab
+      if(type.eq.1) then
+c real emission
+         QCDpower(icontr)=nint(2*wgtbpower+2)
+         do i=1,nexternal
+            do j=0,3
+               momenta(j,i,icontr)=p_ev(j,i)
+            enddo
+            pdg(i,icontr)=idup(i,1)
+         enddo
+      elseif(type.ge.2 .and. type.le.6) then
+c Born, counter term, or soft-virtual
+         if (type.eq.2) then
+            QCDpower(icontr)=nint(2*wgtbpower)
+         else
+            QCDpower(icontr)=nint(2*wgtbpower+2)
+         endif
+         do i=1,nexternal
+            do j=0,3
+               if (p1_cnt(0,1,0).gt.0d0) then
+                  momenta(j,i,icontr)=p1_cnt(j,i,0)
+               elseif (p1_cnt(0,1,1).gt.0d0) then
+                  momenta(j,i,icontr)=p1_cnt(j,i,1)
+               elseif (p1_cnt(0,1,2).gt.0d0) then
+                  momenta(j,i,icontr)=p1_cnt(j,i,2)
+               else
+                  write (*,*) 'ERROR in add_wgt: no valid momenta'
+                  stop 1
+               endif
+            enddo
+            pdg(i,icontr)=idup(i,1)
+         enddo
+      else
+         write (*,*) 'ERROR: unknown type in add_wgt',type
+         stop 1
+      endif
+      return
+      end
+
+      subroutine include_PDF_and_alphas
+c Multiply the saved wgt() info by the PDFs, alpha_S and the scale
+c dependence and saves the weights in the wgts() array. The weights in
+c this array are now correctly normalised to compute the cross section
+c or to fill histograms.
+      implicit none
+      include 'nexternal.inc'
+      include 'run.inc'
+      include 'c_weight.inc'
+      include 'coupl.inc'
+      include 'timing_variables.inc'
+      integer i
+      double precision xlum,dlum,pi,mu2_r,mu2_f,mu2_q,rwgt_muR_dep_fac
+      external rwgt_muR_dep_fac
+      parameter (pi=3.1415926535897932385d0)
+      external dlum
+      integer              nFKSprocess
+      common/c_nFKSprocess/nFKSprocess
+      double precision           virt_wgt_mint,born_wgt_mint
+      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
+      call cpu_time(tBefore)
+      if (icontr.eq.0) return
+      do i=1,icontr
+         nFKSprocess=nFKS(i)
+         xbk(1) = bjx(1,i)
+         xbk(2) = bjx(2,i)
+         mu2_q=scales2(1,i)
+         mu2_r=scales2(2,i)
+         mu2_f=scales2(3,i)
+         q2fact(1)=mu2_f
+         q2fact(2)=mu2_f
+c call the PDFs
+         xlum = dlum()
+c iwgt=1 is the central value (i.e. no scale/PDF reweighting).
+         iwgt=1
+         wgts(iwgt,i)=xlum * (wgt(1,i) + wgt(2,i)*log(mu2_r/mu2_q) +
+     &        wgt(3,i)*log(mu2_f/mu2_q))*g_strong(i)**QCDpower(i)
+         wgts(iwgt,i)=wgts(iwgt,i)*rwgt_muR_dep_fac(sqrt(mu2_r))
+         if (itype(i).eq.3) then
+c Special for the soft-virtual needed for the virt-tricks. The
+c *_wgt_mint variable should be directly passed to the mint-integrator
+c and not be part of the plots nor computation of the cross section.
+            virt_wgt_mint=virt_wgt_mint*xlum*g_strong(i)**QCDpower(i)
+     &           *rwgt_muR_dep_fac(sqrt(scales2(2,i)))
+            born_wgt_mint=born_wgt_mint*xlum*g_strong(i)**QCDpower(i)
+     &           /(8d0*Pi**2)*rwgt_muR_dep_fac(sqrt(mu2_r))
+         endif
+      enddo
+      call cpu_time(tAfter)
+      t_as=t_as+(tAfter-tBefore)
+      return
+      end
+
+      subroutine reweight_scale
+c Use the saved c_weight info to perform scale reweighting. Extends the
+c wgts() array to include the weights.
+      implicit none
+      include 'nexternal.inc'
+      include 'run.inc'
+      include 'c_weight.inc'
+      include 'reweight.inc'
+      include 'reweightNLO.inc'
+      include 'timing_variables.inc'
+      integer i,kr,kf,iwgt_save
+      double precision xlum(maxscales),dlum,pi,mu2_r(maxscales)
+     &     ,mu2_f(maxscales),mu2_q,alphas,g(maxscales),rwgt_muR_dep_fac
+      external rwgt_muR_dep_fac
+      parameter (pi=3.1415926535897932385d0)
+      external dlum,alphas
+      integer              nFKSprocess
+      common/c_nFKSprocess/nFKSprocess
+      double precision           virt_wgt_mint,born_wgt_mint
+      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
+      call cpu_time(tBefore)
+      if (icontr.eq.0) return
+c currently we have 'iwgt' weights in the wgts() array.
+      iwgt_save=iwgt
+c loop over all the contributions in the c_weights common block
+      do i=1,icontr
+         iwgt=iwgt_save
+         nFKSprocess=nFKS(i)
+         xbk(1) = bjx(1,i)
+         xbk(2) = bjx(2,i)
+         mu2_q=scales2(1,i)
+c renormalisation scale variation (requires recomputation of the strong
+c coupling)
+         do kr=1,numscales
+            mu2_r(kr)=scales2(2,i)*yfactR(kr)**2
+            g(kr)=sqrt(4d0*pi*alphas(sqrt(mu2_r(kr))))
+         enddo
+c factorisation scale variation (require recomputation of the PDFs)
+         do kf=1,numscales
+            mu2_f(kf)=scales2(3,i)*yfactF(kf)**2
+            q2fact(1)=mu2_f(kf)
+            q2fact(2)=mu2_f(kf)
+            xlum(kf) = dlum()
+         enddo
+         do kr=1,numscales
+            do kf=1,numscales
+               iwgt=iwgt+1 ! increment the iwgt for the wgts() array
+               if (iwgt.gt.max_wgt) then
+                  write (*,*) 'ERROR too many weights in reweight_scale'
+     &                 ,iwgt,max_wgt
+                  stop 1
+               endif
+c add the weights to the array
+               wgts(iwgt,i)=xlum(kf) * (wgt(1,i)+wgt(2,i)*log(mu2_r(kr)
+     &              /mu2_q)+wgt(3,i)*log(mu2_f(kf)/mu2_q))*g(kr)
+     &              **QCDpower(i)
+               wgts(iwgt,i)=wgts(iwgt,i)
+     &              *rwgt_muR_dep_fac(sqrt(mu2_r(kr)))
+            enddo
+         enddo
+      enddo
+      call cpu_time(tAfter)
+      tr_s=tr_s+(tAfter-tBefore)
+      return
+      end
+
+      subroutine reweight_pdf
+c Use the saved c_weight info to perform PDF reweighting. Extends the
+c wgts() array to include the weights.
+      implicit none
+      include 'nexternal.inc'
+      include 'run.inc'
+      include 'c_weight.inc'
+      include 'reweight.inc'
+      include 'reweightNLO.inc'
+      include 'timing_variables.inc'
+      integer n,izero,i
+      parameter (izero=0)
+      double precision xlum,dlum,pi,mu2_r,mu2_f,mu2_q,rwgt_muR_dep_fac
+      external rwgt_muR_dep_fac
+      parameter (pi=3.1415926535897932385d0)
+      external dlum,alphas
+      integer              nFKSprocess
+      common/c_nFKSprocess/nFKSprocess
+      call cpu_time(tBefore)
+      if (icontr.eq.0) return
+c Use as external loop the one over the PDF sets and as internal the one
+c over the icontr. This reduces the number of calls to InitPDF and
+c allows for better caching of the PDFs
+      do n=1,numPDFs-1
+         iwgt=iwgt+1
+         if (iwgt.gt.max_wgt) then
+            write (*,*) 'ERROR too many weights in reweight_pdf',iwgt
+     &           ,max_wgt
+            stop 1
+         endif
+         call InitPDF(n)
+         do i=1,icontr
+            nFKSprocess=nFKS(i)
+            xbk(1) = bjx(1,i)
+            xbk(2) = bjx(2,i)
+            mu2_q=scales2(1,i)
+            mu2_r=scales2(2,i)
+            mu2_f=scales2(3,i)
+            q2fact(1)=mu2_f
+            q2fact(2)=mu2_f
+            xlum = dlum()
+c add the weights to the array
+            wgts(iwgt,i)=xlum * (wgt(1,i) + wgt(2,i)*log(mu2_r/mu2_q) +
+     &           wgt(3,i)*log(mu2_f/mu2_q))*g_strong(i)**QCDpower(i)
+            wgts(iwgt,i)=wgts(iwgt,i)*rwgt_muR_dep_fac(sqrt(mu2_r))
+         enddo
+         call InitPDF(izero)
+      enddo
+      call cpu_time(tAfter)
+      tr_pdf=tr_pdf+(tAfter-tBefore)
+      return
+      end
+
+      subroutine fill_applgrid_weights(vegas_wgt)
+c Fills the ApplGrid weights of appl_common.inc. This subroutine assumes
+c that there is an unique PS configuration: at most one Born, one real
+c and one set of counter terms. Among other things, this means that one
+c must do MC over FKS directories.
+      implicit none
+      include 'nexternal.inc'
+      include 'c_weight.inc'
+      include 'appl_common.inc'
+      include 'nFKSconfigs.inc'
+      include 'genps.inc'
+      integer i
+      double precision final_state_rescaling,vegas_wgt
+      integer              flavour_map(fks_configs)
+      common/c_flavour_map/flavour_map
+      integer iproc_save(fks_configs),eto(maxproc,fks_configs),
+     &     etoi(maxproc,fks_configs),maxproc_found
+      common/cproc_combination/iproc_save,eto,etoi,maxproc_found
+      if (icontr.gt.6) then
+         write (*,*) 'ERROR: too many applgrid weights. '/
+     &        /'Should have at most one of each itype.',icontr
+         stop 1
+      endif
+      do i=1,4
+         appl_w0(i)=0d0
+         appl_wR(i)=0d0
+         appl_wF(i)=0d0
+         appl_wB(i)=0d0
+         appl_x1(i)=0d0
+         appl_x2(i)=0d0
+         appl_QES2(i)=0d0
+         appl_muR2(i)=0d0
+         appl_muF2(i)=0d0
+      enddo
+      appl_event_weight = 0d0
+      appl_vegaswgt = vegas_wgt
+      if (icontr.eq.0) return
+      do i=1,icontr
+         appl_event_weight=appl_event_weight+wgts(1,i)
+         final_state_rescaling = dble(iproc_save(nFKS(i))) /
+     &        dble(appl_nproc(flavour_map(nFKS(i))))
+         if (itype(i).eq.1) then
+c     real
+            appl_w0(1)=appl_w0(1)+wgt(1,i)*final_state_rescaling
+            appl_x1(1)=bjx(1,i)
+            appl_x2(1)=bjx(2,i)
+            appl_flavmap(1) = flavour_map(nFKS(i))
+            appl_QES2(1)=scales2(1,i)
+            appl_muR2(1)=scales2(2,i)
+            appl_muF2(1)=scales2(3,i)
+         elseif (itype(i).eq.2) then
+c     born
+            appl_wB(2)=appl_wB(2)+wgt(1,i)*final_state_rescaling
+            appl_x1(2)=bjx(1,i)
+            appl_x2(2)=bjx(2,i)
+            appl_flavmap(2) = flavour_map(nFKS(i))
+            appl_QES2(2)=scales2(1,i)
+            appl_muR2(2)=scales2(2,i)
+            appl_muF2(2)=scales2(3,i)
+         elseif (itype(i).eq.3 .or. itype(i).eq.4) then
+c     soft-virtual or soft-counter
+            appl_w0(2)=appl_w0(2)+wgt(1,i)*final_state_rescaling
+            appl_wR(2)=appl_wR(2)+wgt(2,i)*final_state_rescaling
+            appl_wF(2)=appl_wF(2)+wgt(3,i)*final_state_rescaling
+            appl_x1(2)=bjx(1,i)
+            appl_x2(2)=bjx(2,i)
+            appl_flavmap(2) = flavour_map(nFKS(i))
+            appl_QES2(2)=scales2(1,i)
+            appl_muR2(2)=scales2(2,i)
+            appl_muF2(2)=scales2(3,i)
+         elseif (itype(i).eq.5) then
+c     collinear counter            
+            appl_w0(3)=appl_w0(3)+wgt(1,i)*final_state_rescaling
+            appl_wF(3)=appl_wF(3)+wgt(3,i)*final_state_rescaling
+            appl_x1(3)=bjx(1,i)
+            appl_x2(3)=bjx(2,i)
+            appl_flavmap(3) = flavour_map(nFKS(i))
+            appl_QES2(3)=scales2(1,i)
+            appl_muR2(3)=scales2(2,i)
+            appl_muF2(3)=scales2(3,i)
+         elseif (itype(i).eq.6) then
+c     soft-collinear counter            
+            appl_w0(4)=appl_w0(4)+wgt(1,i)*final_state_rescaling
+            appl_wF(4)=appl_wF(4)+wgt(3,i)*final_state_rescaling
+            appl_x1(4)=bjx(1,i)
+            appl_x2(4)=bjx(2,i)
+            appl_flavmap(4) = flavour_map(nFKS(i))
+            appl_QES2(4)=scales2(1,i)
+            appl_muR2(4)=scales2(2,i)
+            appl_muF2(4)=scales2(3,i)
+         endif
+      enddo
+      return
+      end
+      
+      
+      subroutine get_wgt_nbody(sig)
+c Sums all the central weights that contribution to the nbody cross
+c section
+      implicit none
+      include 'nexternal.inc'
+      include 'c_weight.inc'
+      double precision sig
+      integer i
+      sig=0d0
+      if (icontr.eq.0) return
+      do i=1,icontr
+         if (itype(i).eq.2 .or. itype(i).eq.3) then
+            sig=sig+wgts(1,i)
+         endif
+      enddo
+      return
+      end
+
+      subroutine get_wgt_no_nbody(sig)
+c Sums all the central weights that contribution to the cross section
+c excluding the nbody contributions.
+      implicit none
+      include 'nexternal.inc'
+      include 'c_weight.inc'
+      double precision sig
+      integer i
+      sig=0d0
+      if (icontr.eq.0) return
+      do i=1,icontr
+         if (itype(i).eq.1 .or. itype(i).eq.4 .or. itype(i).eq.5 .or.
+     &        itype(i).eq.6) then
+            sig=sig+wgts(1,i)
+         endif
+      enddo
+      return
+      end
+
+      subroutine fill_plots
+c Calls the analysis routine (which fill plots) for all the
+c contributions in the c_weight common block. Instead of really calling
+c it for all, it first checks if weights can be summed (i.e. they have
+c the same PDG codes and the same momenta) before calling the analysis
+c to greatly reduce the calls to the analysis routines.
+      implicit none
+      include 'nexternal.inc'
+      include 'c_weight.inc'
+      include 'reweight0.inc'
+      include 'timing_variables.inc'
+      integer i,ii,j,max_weight
+      logical momenta_equal,pdg_equal
+      external momenta_equal,pdg_equal
+      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
+      double precision www(max_weight)
+      call cpu_time(tBefore)
+      if (icontr.eq.0) return
+c fill the plots_wgts. Check if we can sum weights together before
+c calling the analysis routines. This is the case if the PDG codes and
+c the momenta are identical.
+      do i=1,icontr
+         do j=1,iwgt
+            plot_wgts(j,i)=0d0
+         enddo
+         if (itype(i).eq.2) then
+            plot_id(i)=20 ! Born
+         elseif(itype(1).eq.1) then
+            plot_id(i)=11 ! real-emission
+         else
+            plot_id(i)=12 ! soft-virtual and counter terms
+         endif
+c Loop over all previous icontr. If the plot_id, PDGs and momenta are
+c equal to a previous icountr, add the current weight to the plot_wgts
+c of that contribution and exit the do-loop. This loop extends to 'i',
+c so if the current weight cannot be summed to a previous one, the ii=i
+c contribution makes sure that it is added as a new element.
+         do ii=1,i
+            if (plot_id(ii).eq.plot_id(i)) then
+               if (pdg_equal(pdg(1,ii),pdg(1,i))) then
+                  if (momenta_equal(momenta(0,1,ii),momenta(0,1,i)))then
+                     do j=1,iwgt
+                        plot_wgts(j,ii)=plot_wgts(j,ii)+wgts(j,i)
+                     enddo
+                     exit
+                  endif
+               endif
+            endif
+         enddo
+      enddo
+      do i=1,icontr
+         if (plot_wgts(1,i).ne.0d0) then
+            if (iwgt.gt.max_weight) then
+               write (*,*) 'ERROR too many weights in fill_plots',iwgt
+     &              ,max_weight
+               stop 1
+            endif
+            do j=1,iwgt
+               www(j)=plot_wgts(j,i)
+            enddo
+c call the analysis/histogramming routines
+            call outfun(momenta(0,1,i),y_bst(i),www,pdg(1,i),plot_id(i))
+         endif
+      enddo
+      call cpu_time(tAfter)
+      t_plot=t_plot+(tAfter-tBefore)
+      return
+      end
+
+      subroutine fill_mint_function(f)
+      implicit none
+      include 'nexternal.inc'
+      include 'c_weight.inc'
+      include 'mint.inc'
+      integer i
+      double precision f(nintegrals),sigint
+      double precision           virt_wgt_mint,born_wgt_mint
+      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
+      double precision virtual_over_born
+      common /c_vob/   virtual_over_born
+      sigint=0d0
+      do i=1,icontr
+         sigint=sigint+wgts(1,i)
+      enddo
+      f(1)=abs(sigint)
+      f(2)=sigint
+      f(3)=virt_wgt_mint
+      f(4)=virtual_over_born
+      f(5)=abs(virt_wgt_mint)
+      f(6)=born_wgt_mint
+      return
+      end
+      
+
       subroutine rotate_invar(pin,pout,cth,sth,cphi,sphi)
 c Given the four momentum pin, returns the four momentum pout (in the same
 c Lorentz frame) by performing a three-rotation of an angle theta 
@@ -222,796 +1449,6 @@
       end
 
 
-      double precision function dsig(pp,wgt,vegaswgt)
-c Here are the subtraction terms, the Sij function, 
-c the f-damping function, and the single diagram
-c enhanced multi-channel factor included
-      implicit none
-      include "genps.inc"
-      include 'nexternal.inc'
-c      include "fks.inc"
-      integer fks_j_from_i(nexternal,0:nexternal)
-     &     ,particle_type(nexternal),pdg_type(nexternal)
-      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
-      include "fks_powers.inc"
-      include 'coupl.inc'
-      include 'q_es.inc'
-      include 'run.inc'
-      include 'reweight.inc'
-      include 'reweightNLO.inc'
-
-c     timing statistics
-      include 'timing_variables.inc'
-      real deltaTOLP,deltaTPDF,deltaTFJ
-
-      double precision pp(0:3,nexternal),wgt,vegaswgt
-
-      double precision fks_Sij,f_damp,dot,dlum
-      external fks_Sij,f_damp,dot,dlum
-
-      double precision x,xtot,s_ev,s_c,s_s,s_sc,ffact,fx_ev,fx_c,
-     #                 fx_s,fx_sc,ev_wgt,cnt_wgt_c,cnt_wgt_s,cnt_wgt_sc,
-     #                 bsv_wgt,plot_wgt,cnt_swgt_s,cnt_swgt_sc,cnt_sc,cnt_s,
-     #                 prefact_cnt_ssc,prefact_cnt_ssc_c,prefact_coll,
-     #                 prefact_coll_c,born_wgt,prefact_deg,prefact,prefact_c,
-     #                 prefact_deg_sxi,prefact_deg_slxi,deg_wgt,deg_swgt,
-     #                 deg_xi_c,deg_lxi_c,deg_xi_sc,deg_lxi_sc,
-     #                 cnt_swgt,cnt_wgt,xlum_ev,xlum_c,xlum_s,xlum_sc,xsec,
-     #                 bpower,cpower
-      integer i,j,iplot
-
-      integer izero,ione,itwo,mohdr,iplot_ev,iplot_cnt,iplot_born
-      integer ithree,ifour,ifill2,ifill3,ifill4
-      parameter (izero=0)
-      parameter (ione=1)
-      parameter (itwo=2)
-      parameter (ithree=3)
-      parameter (ifour=4)
-      parameter (mohdr=-100)
-      parameter (iplot_ev=11)
-      parameter (iplot_cnt=12)
-      parameter (iplot_born=20)
-
-      double precision ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
-      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,
-     #                        sqrtshat,shat
-
-      double precision shattmp
-      double precision pi
-      parameter (pi=3.1415926535897932385d0)
-
-      logical nocntevents
-      common/cnocntevents/nocntevents
-
-      double precision xnormsv
-      common/cxnormsv/xnormsv
-
-c Multi channel stuff:
-      Double Precision amp2(maxamps), jamp2(0:maxamps)
-      common/to_amps/  amp2,       jamp2
-
-      integer          isum_hel
-      logical                    multi_channel
-      common/to_matrix/isum_hel, multi_channel
-      INTEGER MAPCONFIG(0:LMAXCONFIGS), ICONFIG
-      common/to_mconfigs/mapconfig, iconfig
-
-      double complex wgt1(2)
-      double precision p_born(0:3,nexternal-1)
-      common/pborn/p_born
-
-      logical calculatedBorn
-      common/ccalculatedBorn/calculatedBorn
-
-      double precision ev_enh,enhance,rwgt,unwgtfun
-      logical firsttime,passcuts
-      data firsttime /.true./
-      integer inoborn_ev,inoborn_cnt
-      double precision xnoborn_ev,xnoborn_cnt
-
-c FKS stuff:
-      integer i_fks,j_fks
-      common/fks_indices/i_fks,j_fks
-
-      double precision xi_i_fks_ev,y_ij_fks_ev
-      double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
-      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
-
-      double precision xi_i_fks_cnt(-2:2)
-      common /cxiifkscnt/xi_i_fks_cnt
-
-      double precision xiimax_ev
-      common /cxiimaxev/xiimax_ev
-      double precision xiimax_cnt(-2:2)
-      common /cxiimaxcnt/xiimax_cnt
-
-      double precision xinorm_ev
-      common /cxinormev/xinorm_ev
-      double precision xinorm_cnt(-2:2)
-      common /cxinormcnt/xinorm_cnt
-
-      double precision p1_cnt(0:3,nexternal,-2:2)
-      double precision wgt_cnt(-2:2)
-      double precision pswgt_cnt(-2:2)
-      double precision jac_cnt(-2:2)
-      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
-
-      double precision zero,one
-      parameter (zero=0d0,one=1d0)
-      double precision tiny
-      parameter (tiny=1d-6)
-
-      double precision xicut_used
-      common /cxicut_used/xicut_used
-      double precision delta_used
-      common /cdelta_used/delta_used
-      double precision xiScut_used,xiBSVcut_used
-      common /cxiScut_used/xiScut_used,xiBSVcut_used
-      double precision fkssymmetryfactor,fkssymmetryfactorBorn,
-     &     fkssymmetryfactorDeg
-      integer ngluons,nquarks(-6:6)
-      common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
-     &                         fkssymmetryfactorDeg,ngluons,nquarks
-      double precision diagramsymmetryfactor
-      common /dsymfactor/diagramsymmetryfactor
-
-      character*4 abrv
-      common /to_abrv/ abrv
-
-      logical nbody
-      common/cnbody/nbody
-
-      double precision vegas_weight
-      common/cvegas_weight/vegas_weight
-
-c For the MINT folding
-      integer fold
-      common /cfl/fold
-
-c Random numbers to be used in the plotting routine
-      logical needrndec
-      parameter (needrndec=.true.)
-
-      real*8 ran2
-      external ran2
-
-      real*8 rndec(10)
-      common/crndec/rndec
-
-c For tests
-      real*8 fksmaxwgt,xisave,ysave
-      common/cfksmaxwgt/fksmaxwgt,xisave,ysave
-
-      logical ExceptPSpoint
-      integer iminmax
-      common/cExceptPSpoint/iminmax,ExceptPSpoint
-
-      double precision central_wgt_saved
-      save central_wgt_saved
-
-      double precision dsig_max,dsig_min
-      double precision total_wgt_sum,total_wgt_sum_max,
-     &                 total_wgt_sum_min
-      common/csum_of_wgts/total_wgt_sum,total_wgt_sum_max,
-     &                 total_wgt_sum_min
-      double precision virt_wgt,born_wgt_ao2pi
-      common/c_fks_singular/virt_wgt,born_wgt_ao2pi
-
-c APPLgrid stuff
-      integer iappl
-      common /for_applgrid/ iappl
-      include "appl_common.inc" 
-
-c FxFx merging
-      logical rewgt_mohdr_calculated,rewgt_izero_calculated
-      double precision rewgt_mohdr,rewgt_izero,rewgt_exp_mohdr
-     $     ,rewgt_exp_izero,enhanceS,enhanceH
-      logical setclscales
-      double precision rewgt
-      external setclscales,rewgt
-
-      double precision pmass(nexternal)
-      double precision rwgt_muR_dep_fac
-      include "pmass.inc"
-
-      vegas_weight=vegaswgt
-
-c FxFx merging
-      ktscheme=1
-      rewgt_mohdr_calculated=.false.
-      rewgt_izero_calculated=.false.
-
-
-      if (fold.eq.0) then
-         calculatedBorn=.false.
-         call get_helicity(i_fks,j_fks)
-      endif
-
-c Make sure that the result can be non-zero. If the jacobian from the
-c PS-setup or vegas are zero, we can skip this PS point and 'return'.
-c Note that all the wgts and jacs should be positive.
-      dsig=0d0
-      if ( (wgt.le.0d0 .and. jac_cnt(0).le.0d0 .and. jac_cnt(1).le.0d0
-     &     .and. jac_cnt(2).le.0d0) .or. vegaswgt.le.0d0) return
-
-      if (firsttime)then
-         inoborn_ev=0
-         xnoborn_ev=0.d0
-         inoborn_cnt=0
-         xnoborn_cnt=0.d0
-         fksmaxwgt=0.d0
-c Put here call to compute bpower
-         call compute_bpower(p_born,bpower)
-         wgtbpower=bpower
-c Store the power of alphas of the Born events in the appl common block.
-         if(iappl.ne.0) appl_bpower = wgtbpower
-c Initialize histograms
-         call initplot
-c Compute cpower done for bottom Yukawa, routine needs to be adopted
-c for other muR-dependendent factors
-         call compute_cpower(p_born,cpower)
-         if(dabs(cpower+1d0).lt.tiny) then
-            wgtcpower=0d0
-         else
-            wgtcpower=cpower
-         endif
-c Check that things are done consistently
-         if(wgtcpower.ne.cpowerinput.and.dabs(cpower+1d0).gt.tiny)then
-           write(*,*)'Inconsistency in the computation of cpower',
-     #               wgtcpower,cpowerinput
-           write(*,*)'Check value in reweight0.inc'
-           stop
-         endif
-
-         firsttime=.false.
-      endif
-
-      prefact=xinorm_ev/xi_i_fks_ev*
-     #        1/(1-y_ij_fks_ev)
-      if( (.not.nocntevents) .and. (.not.(abrv.eq.'born' .or. abrv.eq
-     &     .'grid' .or. abrv(1:2).eq.'vi' .or. nbody))
-     &     )then
-        prefact_cnt_ssc=xinorm_ev/min(xiimax_ev,xiScut_used)*
-     #                  log(xicut_used/min(xiimax_ev,xiScut_used))*
-     #                  1/(1-y_ij_fks_ev)
-        if(pmass(j_fks).eq.0.d0)then
-          prefact_c=xinorm_cnt(ione)/xi_i_fks_cnt(ione)*
-     #              1/(1-y_ij_fks_ev)
-          prefact_cnt_ssc_c=xinorm_cnt(ione)/min(xiimax_cnt(ione),xiScut_used)*
-     #                      log(xicut_used/min(xiimax_cnt(ione),xiScut_used))*
-     #                      1/(1-y_ij_fks_ev)
-          prefact_coll=xinorm_cnt(ione)/xi_i_fks_cnt(ione)*
-     #                 log(delta_used/deltaS)/deltaS
-          prefact_coll_c=xinorm_cnt(ione)/min(xiimax_cnt(ione),xiScut_used)*
-     #                   log(xicut_used/min(xiimax_cnt(ione),xiScut_used))*
-     #                   log(delta_used/deltaS)/deltaS
-          prefact_deg=xinorm_cnt(ione)/xi_i_fks_cnt(ione)*
-     #                1/deltaS
-          prefact_deg_sxi=xinorm_cnt(ione)/min(xiimax_cnt(ione),xiScut_used)*
-     #                    log(xicut_used/min(xiimax_cnt(ione),xiScut_used))*
-     #                    1/deltaS
-          prefact_deg_slxi=xinorm_cnt(ione)/min(xiimax_cnt(ione),xiScut_used)*
-     #                     ( log(xicut_used)**2 -
-     #                       log(min(xiimax_cnt(ione),xiScut_used))**2 )*
-     #                     1/(2.d0*deltaS)
-        endif
-      endif
-
-      if(needrndec)then
-        do i=1,10
-          rndec(i)=ran2()
-        enddo
-      endif
-
-c If there was an exceptional phase-space point found for the 
-c virtual corrections, at the end of this subroutine, goto 44
-c and compute also the "possible" minimal and maximal weight
-c these points could have gotton (based upon previous PS
-c points)
-      ExceptPSpoint=.false.
-      iminmax=-1
- 44   continue
-      iminmax=iminmax+1
-
-      ev_wgt=0.d0
-      cnt_wgt=0.d0
-      cnt_wgt_s=0.d0
-      cnt_wgt_c=0.d0
-      cnt_wgt_sc=0.d0
-      bsv_wgt=0.d0
-      virt_wgt=0d0
-      born_wgt=0.d0
-      born_wgt_ao2pi=0.d0
-      cnt_swgt=0.d0
-      cnt_swgt_s=0.d0
-      cnt_swgt_sc=0.d0
-      deg_wgt=0.d0
-      deg_swgt=0.d0
-      plot_wgt=0.d0
-      iplot=-3
-
-      if(doreweight)then
-        call reweight_settozero()
-        ifill2=0
-        ifill3=0
-        ifill4=0
-      endif
-      if (abrv.eq.'born' .or. abrv.eq.'grid' .or. abrv(1:2).eq.'vi' .or.
-     &     nbody)goto 540
-
-      call cpu_time(tBefore)
-      deltaTPDF = tPDF
-      deltaTFJ  = tFastJet
-c Real contribution:
-c Set the ybst_til_tolab before applying the cuts. 
-      call set_cms_stuff(mohdr)
-      if(doreweight)then
-        wgtxbj(1,1)=xbk(1)
-        wgtxbj(2,1)=xbk(2)
-      endif
-      if (passcuts(pp,rwgt)) then
-c       Compute the scales and sudakov-reweighting for the FxFx merging
-        if (ickkw.eq.3) then
-           if (.not. setclscales(pp)) then
-               write (*,*) 'ERROR in setclscales mohdr'
-              stop
-           endif
-           rewgt_mohdr=rewgt(pp,rewgt_exp_mohdr)
-           rewgt_mohdr_calculated=.true.
-        endif
-        call set_alphaS(pp)
-        x = abs(2d0*dot(pp(0,i_fks),pp(0,j_fks))/shat)
-        ffact = f_damp(x)
-        s_ev = fks_Sij(pp,i_fks,j_fks,xi_i_fks_ev,y_ij_fks_ev)
-        if(s_ev.gt.0.d0)then
-           call sreal(pp,xi_i_fks_ev,y_ij_fks_ev,fx_ev)
-           xlum_ev = dlum()
-           xsec = fx_ev*s_ev*ffact*wgt*prefact*rwgt
-           ev_wgt = xlum_ev*xsec * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-           if(doreweight)then
-             wgtmuR2(1)=muR2_current/muR_over_ref**2
-             wgtmuF12(1)=muF12_current/muF1_over_ref**2
-             wgtmuF22(1)=muF22_current/muF2_over_ref**2
-             call reweight_fillkin(pp,ione)
-             wgtwreal(1)=xsec/g**(nint(2*wgtbpower+2.d0))
-           endif
-        endif
-      endif
-      call cpu_time(tAfter)
-      tDSigR=tDSigR + (tAfter-tBefore) - (tPDF-deltaTPDF) -
-     &     (tFastJet-deltaTFJ)
-c
-c All counterevent have the same final-state kinematics. Check that
-c one of them passes the hard cuts, and they exist at all
- 540  continue
-      call cpu_time(tBefore)
-      deltaTPDF = tPDF
-      deltaTFJ  = tFastJet
-      deltaTOLP = tOLP
-c Set the ybst_til_tolab before applying the cuts. Update below
-c for the collinear, soft and/or soft-collinear subtraction terms
-      call set_cms_stuff(izero)
-      if ( (.not.passcuts(p1_cnt(0,1,0),rwgt)) .or.
-     #     nocntevents ) goto 547
-
-c Compute the scales and sudakov-reweighting for the FxFx merging
-      if (ickkw.eq.3) then
-         if (.not. setclscales(p1_cnt(0,1,0))) then
-            write (*,*) 'ERROR in setclscales izero'
-            stop
-         endif
-         rewgt_izero=rewgt(p1_cnt(0,1,0),rewgt_exp_izero)
-         rewgt_izero_calculated=.true.
-      endif
-      call set_alphaS(p1_cnt(0,1,0))
-      if(doreweight)then
-        wgtqes2(2)=QES2
-        wgtqes2(3)=QES2
-        wgtqes2(4)=QES2
-        wgtmuR2(2)=muR2_current/muR_over_ref**2
-        wgtmuF12(2)=muF12_current/muF1_over_ref**2
-        wgtmuF22(2)=muF22_current/muF2_over_ref**2
-        call reweight_fillkin(pp,itwo)
-        ifill2=1
-      endif
-
-      if (abrv.eq.'born' .or. abrv.eq.'grid' .or. abrv(1:2).eq.'vi' .or.
-     &     nbody)goto 545
-c
-c Collinear subtraction term:
-      if( y_ij_fks_ev.gt.1d0-deltaS .and.
-     #    pmass(j_fks).eq.0.d0 )then
-         call set_cms_stuff(ione)
-         if(doreweight)then
-           wgtxbj(1,3)=xbk(1)
-           wgtxbj(2,3)=xbk(2)
-         endif
-         s_c = fks_Sij(p1_cnt(0,1,1),i_fks,j_fks,xi_i_fks_cnt(ione),one)
-         if(s_c.gt.0.d0)then
-            if(abs(s_c-1.d0).gt.1.d-6.and.j_fks.le.nincoming)then
-              write(*,*)'Wrong S function in dsig[c]',s_c
-              stop
-            endif
-            call sreal(p1_cnt(0,1,1),xi_i_fks_cnt(ione),one,fx_c)
-            xlum_c = dlum()
-            xsec = fx_c*s_c*jac_cnt(1)*(prefact_c+prefact_coll)*rwgt
-            cnt_wgt_c=cnt_wgt_c-xlum_c*xsec * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-            call sreal_deg(p1_cnt(0,1,1),xi_i_fks_cnt(ione),one,
-     #                     deg_xi_c,deg_lxi_c)
-            deg_wgt=deg_wgt+( deg_xi_c+deg_lxi_c*log(xi_i_fks_cnt(ione)) )*
-     #                      jac_cnt(1)*prefact_deg*rwgt/(shat/(32*pi**2))*
-     #                      xlum_c * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-            iplot=1
-            if(doreweight)then
-              call reweight_fillkin(pp,ithree)
-              ifill3=1
-              wgtwreal(3)=-xsec/g**(nint(2*wgtbpower+2.d0))
-              wgtwdeg(3)=
-     #              ( wgtdegrem_xi+wgtdegrem_lxi*log(xi_i_fks_cnt(ione)) )*
-     #                      jac_cnt(1)*prefact_deg*rwgt/(shat/(32*pi**2))/
-     #              g**(nint(2*wgtbpower+2.d0))
-              wgtwdegmuf(3)=wgtdegrem_muF *
-     #                      jac_cnt(1)*prefact_deg*rwgt/(shat/(32*pi**2))/
-     #              g**(nint(2*wgtbpower+2.d0))
-            endif
-         endif
-      endif
-c Soft subtraction term:
- 545  continue
-      if (xi_i_fks_ev .lt. max(xiScut_used,xiBSVcut_used)) then
-         call set_cms_stuff(izero)
-         if(doreweight)then
-           wgtxbj(1,2)=xbk(1)
-           wgtxbj(2,2)=xbk(2)
-           if(ifill2.ne.1)then
-             write(*,*)'Error #1[wg] in dsig',ifill2
-             stop
-           endif
-         endif
-         s_s = fks_Sij(p1_cnt(0,1,0),i_fks,j_fks,zero,y_ij_fks_ev)
-         if(nbody)s_s=1.d0
-         if(s_s.gt.0.d0)then
-            xlum_s = dlum()
-            if (abrv.eq.'born' .or. abrv.eq.'grid' .or.
-     &           abrv(1:2).eq.'vi' .or. nbody) goto 546
-            if (xi_i_fks_ev .lt. xiScut_used) then
-              call sreal(p1_cnt(0,1,0),zero,y_ij_fks_ev,fx_s)
-              xsec=fx_s*s_s*jac_cnt(0)
-              cnt_s=xlum_s*xsec
-              cnt_wgt_s=cnt_wgt_s-cnt_s*prefact*rwgt
-     f             * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-              cnt_swgt_s=cnt_swgt_s-cnt_s*prefact_cnt_ssc*rwgt
-     f             * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-              if(doreweight)
-     #          wgtwreal(2)=-xsec*(prefact+prefact_cnt_ssc)*rwgt/
-     #                      g**(nint(2*wgtbpower+2.d0))
-            endif
- 546        continue
-            if (abrv.eq.'real' .or. .not.nbody) goto 548
-            if (xi_i_fks_ev .lt. xiBSVcut_used) then
-              xsec=s_s*jac_cnt(0)*xinorm_ev/
-     #             (min(xiimax_ev,xiBSVcut_used)*shat/(16*pi**2))*
-     #             rwgt
-              xnormsv=xlum_s*xsec
-              call bornsoftvirtual(p1_cnt(0,1,0),bsv_wgt,virt_wgt
-     $             ,born_wgt)
-              born_wgt_ao2pi=born_wgt*g**2/(8d0*Pi**2)
-c For FxFx merging, include the compensation term
-              if (rewgt_izero_calculated.and.rewgt_izero.lt.1d0) then
-                 bsv_wgt=bsv_wgt-g**2/(4d0*Pi)*rewgt_exp_izero*born_wgt
-              endif
-              if(doreweight)then
-                if(wgtbpower.gt.0)then
-                  wgtwborn(2)=born_wgt*xsec/g**(nint(2*wgtbpower))
-                else
-                  wgtwborn(2)=born_wgt*xsec
-                endif
-                wgtwns(2)=wgtnstmp*xsec/g**(nint(2*wgtbpower+2.d0))
-                 if (rewgt_izero_calculated.and.rewgt_izero.lt.1d0) then
-                    wgtwns(2)=wgtwns(2)-rewgt_exp_izero*wgtwborn(2)
-     $                   /(4d0*pi)
-                 endif
-                wgtwnsmuf(2)=wgtwnstmpmuf*xsec/g**(nint(2*wgtbpower
-     &               +2.d0))
-                wgtwnsmur(2)=wgtwnstmpmur*xsec/g**(nint(2*wgtbpower
-     &               +2.d0))
-              endif
-              bsv_wgt=bsv_wgt*xnormsv * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-              virt_wgt=virt_wgt*xnormsv * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-              born_wgt=born_wgt*xnormsv * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-              born_wgt_ao2pi=born_wgt_ao2pi*xnormsv
-     f             * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-            endif
- 548        continue
-            iplot=0
-         endif
-      endif
-c Soft-Collinear subtraction term:
-      if (abrv.eq.'born' .or. abrv.eq.'grid' .or. abrv(1:2).eq.'vi' .or.
-     &     nbody)goto 547
-      if (xi_i_fks_cnt(ione) .lt. xiScut_used .and.
-     #    y_ij_fks_ev .gt. 1d0-deltaS .and.
-     #    pmass(j_fks).eq.0.d0 )then
-         call set_cms_stuff(itwo)
-         if(doreweight)then
-           wgtxbj(1,4)=xbk(1)
-           wgtxbj(2,4)=xbk(2)
-         endif
-         s_sc = fks_Sij(p1_cnt(0,1,2),i_fks,j_fks,zero,one)
-         if(s_sc.gt.0.d0)then
-            if(abs(s_sc-1.d0).gt.1.d-6.and.j_fks.le.nincoming)then
-              write(*,*)'Wrong S function in dsig[sc]',s_sc
-              stop
-            endif
-            xlum_sc = dlum()
-            call sreal(p1_cnt(0,1,2),zero,one,fx_sc)
-            xsec=fx_sc*s_sc*jac_cnt(2)
-            cnt_sc=xlum_sc*xsec
-            cnt_wgt_sc=cnt_wgt_sc+cnt_sc*(prefact_c+prefact_coll)*rwgt
-     f           * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-            cnt_swgt_sc=cnt_swgt_sc+
-     &           cnt_sc*(prefact_cnt_ssc_c+prefact_coll_c)*rwgt
-     f           * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-            call sreal_deg(p1_cnt(0,1,2),zero,one,
-     #                     deg_xi_sc,deg_lxi_sc)
-            deg_wgt=deg_wgt-( deg_xi_sc+deg_lxi_sc*log(xi_i_fks_cnt(ione)) )*
-     #                     jac_cnt(2)*prefact_deg*rwgt/(shat/(32*pi**2))*
-     #                     xlum_sc * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-            deg_swgt=deg_swgt-( deg_xi_sc*prefact_deg_sxi +
-     #                     deg_lxi_sc*prefact_deg_slxi )*
-     #                     jac_cnt(2)*rwgt/(shat/(32*pi**2))*
-     #                     xlum_sc * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-            if(iplot.ne.0)iplot=2
-            if(doreweight)then
-              call reweight_fillkin(pp,ifour)
-              ifill4=1
-              wgtwreal(4)=xsec*(prefact_c+prefact_coll+
-     #                          prefact_cnt_ssc_c+prefact_coll_c)*rwgt/
-     #                    g**(nint(2*wgtbpower+2.d0))
-              wgtwdeg(4)=(
-     # -( wgtdegrem_xi+wgtdegrem_lxi*log(xi_i_fks_cnt(ione)) )*prefact_deg
-     # -( wgtdegrem_xi*prefact_deg_sxi+wgtdegrem_lxi*prefact_deg_slxi ) )*
-     #                     jac_cnt(2)*rwgt/(shat/(32*pi**2))/
-     #              g**(nint(2*wgtbpower+2.d0))
-              wgtwdegmuf(4)=
-     #          -wgtdegrem_muF*( prefact_deg+prefact_deg_sxi )*
-     #                     jac_cnt(2)*rwgt/(shat/(32*pi**2))/
-     #              g**(nint(2*wgtbpower+2.d0))
-            endif
-         endif
-      endif
- 547  continue
-c
-      call cpu_time(tAfter)
-      tDSigI=tDSigI + (tAfter-tBefore) - (tPDF-deltaTPDF) -
-     &     (tFastJet-deltaTFJ) - (tOLP-deltaTOLP)
-c
-c Enhance the one channel for multi-channel integration
-c
-      enhance=1.d0
-      if ((ev_wgt.ne.0d0 .or. cnt_wgt_c.ne.0d0 .or. cnt_wgt_s.ne.0d0
-     $     .or.cnt_wgt_sc.ne.0d0 .or. bsv_wgt.ne.0d0 .or.
-     $     virt_wgt.ne.0d0 .or. deg_wgt.ne.0d0 .or.deg_swgt.ne.0d0 .or.
-     $     cnt_swgt_s.ne.0d0 .or. cnt_swgt_sc.ne.0d0).and.
-     $     multi_channel) then
-         if (bsv_wgt.eq.0d0 .and. virt_wgt.eq.0d0 .and. deg_wgt.eq.0d0
-     $        .and. deg_swgt.eq.0d0 .and. cnt_wgt_c.eq.0d0 )
-     $        CalculatedBorn=.false.
-         if (.not.calculatedBorn .and. p_born(0,1).gt.0d0)then
-            call sborn(p_born,wgt1)
-         elseif(p_born(0,1).lt.0d0)then
-            enhance=0d0
-         endif
-
-         if (enhance.eq.0d0)then
-            xnoborn_cnt=xnoborn_cnt+1.d0
-            if(log10(xnoborn_cnt).gt.inoborn_cnt)then
-               write (*,*) 
-     #           'Function dsig: no Born momenta more than 10**',
-     #           inoborn_cnt,'times'
-               inoborn_cnt=inoborn_cnt+1
-            endif
-         else
-            xtot=0d0
-            if (mapconfig(0).eq.0) then
-               write (*,*) 'Fatal error in dsig, no Born diagrams '
-     &           ,mapconfig,'. Check bornfromreal.inc'
-               write (*,*) 'Is fks_singular compiled correctly?'
-               stop
-            endif
-            do i=1, mapconfig(0)
-               xtot=xtot+amp2(mapconfig(i))
-            enddo
-            if (xtot.ne.0d0) then
-               enhance=amp2(mapconfig(iconfig))/xtot
-               enhance=enhance*diagramsymmetryfactor
-            else
-               enhance=0d0
-            endif
-         endif
-      endif
-
-      cnt_wgt = cnt_wgt_c + cnt_wgt_s + cnt_wgt_sc
-      cnt_swgt = cnt_swgt_s + cnt_swgt_sc
-
-c Apply the FxFx Sudakov damping on the H events
-      if (ev_wgt.ne.0d0 .and. ickkw.eq.3..and.
-     $     .not.rewgt_mohdr_calculated) then
-         write (*,*) 'Error rewgt_mohdr_calculated'
-         stop
-      elseif(rewgt_mohdr_calculated) then
-         if (rewgt_mohdr.gt.1d0) rewgt_mohdr=1d0
-         enhanceH=enhance*rewgt_mohdr
-      else
-         enhanceH=enhance
-      endif
-
-      ev_wgt = ev_wgt * enhanceH
-
-c Apply the FxFx Sudakov damping on the S events
-      if(.not.(cnt_wgt.eq.0d0 .and. cnt_swgt.eq.0d0 .and. bsv_wgt.eq.0d0
-     $     .and. born_wgt.eq.0d0 .and. deg_wgt.eq.0d0 .and.
-     $     deg_swgt.eq.0d0 )
-     $     .and. ickkw.eq.3 .and. .not.rewgt_izero_calculated) then
-         write (*,*) 'Error rewgt_izero_calculated'
-         stop
-      elseif(rewgt_izero_calculated) then
-         if (rewgt_izero.gt.1d0) rewgt_izero=1d0
-         enhanceS=enhance*rewgt_izero
-      else
-         enhanceS=enhance
-      endif
-
-      cnt_wgt = cnt_wgt * enhanceS
-      cnt_swgt = cnt_swgt * enhanceS
-      bsv_wgt = bsv_wgt * enhanceS
-      virt_wgt = virt_wgt * enhanceS
-      born_wgt = born_wgt * enhanceS
-      born_wgt_ao2pi = born_wgt_ao2pi * enhanceS
-      deg_wgt = deg_wgt * enhanceS
-      deg_swgt = deg_swgt * enhanceS
-
-      if(iminmax.eq.0) then
-         dsig = (ev_wgt+cnt_wgt)*fkssymmetryfactor +
-     &        cnt_swgt*fkssymmetryfactor +
-     &        bsv_wgt*fkssymmetryfactorBorn +
-     &        virt_wgt*fkssymmetryfactorBorn +
-     &        deg_wgt*fkssymmetryfactorDeg +
-     &        deg_swgt*fkssymmetryfactorDeg
-
-         if (dsig.ne.dsig) then
-            write (*,*) 'ERROR #51 in dsig:',dsig,'skipping event'
-            dsig=0d0
-            return
-         endif
-
-         total_wgt_sum=total_wgt_sum+dsig*vegaswgt
-         central_wgt_saved=dsig
-
-         call unweight_function(p_born,unwgtfun)
-         dsig=dsig*unwgtfun
-         virt_wgt=virt_wgt*unwgtfun*fkssymmetryfactorBorn*vegaswgt
-         born_wgt_ao2pi=born_wgt_ao2pi*unwgtfun*fkssymmetryfactorBorn
-     $        *vegaswgt
-         if(doreweight)then
-           if(ifill2.eq.0.and.(ifill3.ne.0.or.ifill4.ne.0))then
-             write(*,*)'Error #2[wg] in dsig',ifill2,ifill3,ifill4
-             stop
-           endif
-           wgtref = dsig
-           xsec = enhanceS*unwgtfun
-           do i=1,4
-              if (i.eq.1) then
-                 wgtwreal(i)=wgtwreal(i) * xsec*fkssymmetryfactor
-     $                *enhanceH/enhanceS
-              else
-                 wgtwreal(i)=wgtwreal(i) * xsec*fkssymmetryfactor
-              endif
-             wgtwdeg(i)=wgtwdeg(i) * xsec*fkssymmetryfactorDeg
-             wgtwdegmuf(i)=wgtwdegmuf(i) * xsec*fkssymmetryfactorDeg
-             if(i.eq.2)then
-               wgtwborn(i)=wgtwborn(i) * xsec*fkssymmetryfactorBorn
-               wgtwns(i)=wgtwns(i) * xsec*fkssymmetryfactorBorn
-               wgtwnsmuf(i)=wgtwnsmuf(i) * xsec*fkssymmetryfactorBorn
-               wgtwnsmur(i)=wgtwnsmur(i) * xsec*fkssymmetryfactorBorn
-             endif
-           enddo
-           call reweight_fill_extra()
-           if(check_reweight.and.doreweight)
-     #       call check_rwgt_wgt("NLO")
-           if(do_rwgt_scale.or.do_rwgt_pdf)call fill_rwgt_NLOplot()
-c Example of reweighted cross section (scale changed)
-c           dsig_new=compute_rwgt_wgt_NLO(new_muR_fact,new_muF1_fact,
-c     #                                   new_muF2_fact,new_QES_fact,
-c     #                                   iwgtinfo)
-         endif      
-
-c For tests
-        if(abs(dsig).gt.fksmaxwgt)then
-            fksmaxwgt=abs(dsig)
-            xisave=xi_i_fks_ev
-            ysave=y_ij_fks_ev
-         endif
-
-c Plot observables for event
-         plot_wgt=ev_wgt*fkssymmetryfactor*vegaswgt
-         if(abs(plot_wgt).gt.1.d-20.and.pp(0,1).ne.-99d0)
-     &        call outfun(pp,ybst_til_tolab,plot_wgt,iplot_ev)
-c Plot observables for counterevents and Born
-         plot_wgt=( cnt_wgt*fkssymmetryfactor +
-     &              cnt_swgt*fkssymmetryfactor +
-     &              (bsv_wgt-born_wgt)*fkssymmetryfactorBorn +
-     &              deg_wgt*fkssymmetryfactorDeg +
-     &              deg_swgt*fkssymmetryfactorDeg )*vegaswgt +
-     &              virt_wgt/unwgtfun
-         if(abs(plot_wgt).gt.1.d-20) then
-            if(iplot.eq.-3)then
-               write(*,*)'Error #1 in dsig'
-               stop
-            elseif (p1_cnt(0,1,iplot).ne.-99d0)then
-               call outfun(p1_cnt(0,1,iplot),ybst_til_tolab,plot_wgt,
-     &              iplot_cnt)
-            endif
-         endif
-c Plot observables for Born; pass cnt momenta assuming they are
-c identical to Born ones
-         plot_wgt=born_wgt*fkssymmetryfactorBorn*vegaswgt
-         if(abs(plot_wgt).gt.1.d-20) then
-            if(iplot.eq.-3)then
-               write(*,*)'Error #1 in dsig'
-               stop
-            elseif (p1_cnt(0,1,iplot).ne.-99d0)then
-               call outfun(p1_cnt(0,1,iplot),ybst_til_tolab,plot_wgt,
-     &              iplot_born)
-            endif
-         endif
-      elseif (iminmax.eq.1 .and. ExceptPSpoint) then
-c for except PS points, this is the maximal approx for the virtual         
-         call unweight_function(p_born,unwgtfun)
-         dsig_max = ((ev_wgt+cnt_wgt)*fkssymmetryfactor +
-     &        cnt_swgt*fkssymmetryfactor +
-     &        bsv_wgt*fkssymmetryfactorBorn +
-     &        deg_wgt*fkssymmetryfactorDeg +
-     &        deg_swgt*fkssymmetryfactorDeg)*unwgtfun+virt_wgt
-         total_wgt_sum_max=total_wgt_sum_max+
-     &        ((dsig_max - central_wgt_saved)*vegaswgt)**2
-
-      elseif (iminmax.eq.2 .and. ExceptPSpoint) then
-c for except PS points, this is the minimal approx for the virtual         
-         call unweight_function(p_born,unwgtfun)
-         dsig_min = ((ev_wgt+cnt_wgt)*fkssymmetryfactor +
-     &        cnt_swgt*fkssymmetryfactor +
-     &        bsv_wgt*fkssymmetryfactorBorn +
-     &        deg_wgt*fkssymmetryfactorDeg +
-     &        deg_swgt*fkssymmetryfactorDeg)*unwgtfun+virt_wgt
-         total_wgt_sum_min=total_wgt_sum_min+
-     &        ((central_wgt_saved - dsig_min)*vegaswgt)**2
-      else
-         write (*,*) 'Error #12 in dsig',iminmax
-         stop
-      endif
-
-c If exceptional PS point found, go back to beginning recompute
-c the weight for this PS point using an approximation
-c based on previous PS points (done in BinothLHA.f)
-      if (ExceptPSpoint .and. iminmax.le.1) goto 44
-
-      return
-      end
-
 
       subroutine dsigF(pp,wgt,vegaswgt,dsigS,dsigH)
 c Here are the subtraction terms, the Sij function, 
@@ -1220,6 +1657,8 @@
       double precision ximin
       parameter(ximin=0.05d0)
 
+      double precision virtual_over_born
+      common/c_vob/virtual_over_born
 c
 c This is the table that will be used to unweight. (It contains for
 c arguments, 1st argument: nFKSproces; 2nd argument: S or H events; 3rd
@@ -1233,7 +1672,7 @@
       DOUBLE PRECISION       CONV
       PARAMETER (CONV=389379660d0)  !CONV TO PICOBARNS
       integer i_process
-      common/c_addwrite/i_process
+      common/c_i_process/i_process
       integer iproc_save(fks_configs),eto(maxproc,fks_configs)
      $     ,etoi(maxproc,fks_configs),maxproc_found
       common/cproc_combination/iproc_save,eto,etoi,maxproc_found
@@ -1401,8 +1840,16 @@
      #                     1/(2.d0*deltaS)
         endif
       endif
-c
-      probne=1.d0
+
+c For UNLOPS all real-emission contributions need to be added to the
+c S-events. Do this by setting probne to 0. For UNLOPS, no MC counter
+c events are called, so this will remain 0.
+      if (ickkw.eq.4) then
+         probne=0d0
+      else
+         probne=1.d0
+      endif
+
 c All counterevent have the same final-state kinematics. Check that
 c one of them passes the hard cuts, and they exist at all
 c
@@ -1995,15 +2442,8 @@
              unwgt_table(nFKSprocess,1,j)=unwgt_table(nFKSprocess,1,j)
      &            +PD(j)*xsec*(1-probne)*CONV * rwgt_muR_dep_fac(scale)
 com-- muR-dependent fac is reweighted here
-             if (ickkw.ne.4) then
-                unwgt_table(nFKSprocess,2,j)=unwgt_table(nFKSprocess,2,j)
-     &               +PD(j)*xsec*probne*CONV * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-             else         ! for UNLOPS, add the H-events to the S-events
-                unwgt_table(nFKSprocess,1,j)=unwgt_table(nFKSprocess,1,j)
-     &               +PD(j)*xsec*probne*CONV * rwgt_muR_dep_fac(scale)
-com-- muR-dependent fac is reweighted here
-             endif
+             unwgt_table(nFKSprocess,2,j)=unwgt_table(nFKSprocess,2,j)
+     &            +PD(j)*xsec*probne*CONV * rwgt_muR_dep_fac(scale)
           enddo
           if(doreweight)then
              if(ifill1H.eq.0)then
@@ -2164,11 +2604,6 @@
      &        deg_wgt*fkssymmetryfactorDeg +
      &        deg_swgt*fkssymmetryfactorDeg
 
-c Add the H-events to the S-events for UNLOPS
-         if (ickkw.eq.4) then
-            dsigS=dsigS+totH_wgt*fkssymmetryfactor
-         endif
-
          call unweight_function(p_born,unwgtfun)
          dsigS=dsigS*unwgtfun
 
@@ -2183,6 +2618,7 @@
                   unwgt_table(0,1,j)=0d0
                   unwgt_table(0,3,j)=0d0
                   unwgt_table(1,3,j)=0d0
+                  virtual_over_born=0d0
                endif
             enddo
          endif
@@ -2283,7 +2719,7 @@
      &                 ,nFKSprocess*2-1) * xsec*fkssymmetryfactor
                enddo
             endif
-            if(check_reweight.and.doreweight .and. ickkw.ne.4) then
+            if(check_reweight.and.doreweight) then
                do i_process=1,iproc_save(nFKSprocess)
                   if (nbody) then
                      call fill_reweight0inc_nbody(i_process)
@@ -2309,25 +2745,6 @@
             xisave=xi_i_fks_ev
             ysave=y_ij_fks_ev
          endif
-c Plot observables for counterevents and Born
-         if (.not.unwgt) then
-            plot_wgt=( (Sev_wgt+Sxmc_wgt+cnt_wgt)*fkssymmetryfactor +
-     &           cnt_swgt*fkssymmetryfactor +
-     &           (bsv_wgt-born_wgt)*fkssymmetryfactorBorn +
-     &           virt_wgt*fkssymmetryfactorBorn +
-     &           deg_wgt*fkssymmetryfactorDeg +
-     &           deg_swgt*fkssymmetryfactorDeg )*vegaswgt
-            if( abs(plot_wgt).gt.1.d-20.and.p1_cnt(0,1,0).ne.-99d0 .and.
-     &           (plotEv.or.plotKin) )
-     &           call outfun(p1_cnt(0,1,0),ybst_til_tolab,plot_wgt,iplot_cnt)
-            plot_wgt=(born_wgt*fkssymmetryfactorBorn)*vegaswgt
-            if( abs(plot_wgt).gt.1.d-20.and.p1_cnt(0,1,0).ne.-99d0 .and.
-     &           (plotEv.or.plotKin) )
-     &           call outfun(p1_cnt(0,1,0),ybst_til_tolab,plot_wgt,iplot_born)
-         endif
-
-c For UNLOPS, the H-events are already added to the S-events
-         if (ickkw.eq.4) return
 
          dsigH = totH_wgt*fkssymmetryfactor
          call unweight_function(p_born,unwgtfun)
@@ -2341,7 +2758,7 @@
                do j=1,iproc_save(nFKSprocess)
                   unwgt_table(nFKSprocess,2,j)=0d0
                enddo
-         endif
+            endif
          endif
 
          if (nbody.and.dsigH.ne.0d0) then
@@ -2419,15 +2836,6 @@
             ysave=y_ij_fks_ev
          endif
 
-c Plot observables for event
-         if (.not.unwgt) then
-            plot_wgt=totH_wgt*fkssymmetryfactor*vegaswgt 
-            if( abs(plot_wgt).gt.1.d-20.and.pp(0,1).ne.-99d0. and.
-     &           (plotEv.or.plotKin) )
-     &           call outfun(pp,ybst_til_tolab,plot_wgt,iplot_ev)
-         endif
-
-
       elseif (iminmax.eq.1 .and. ExceptPSpoint) then
 c for except PS points, this is the maximal approx for the virtual         
          call unweight_function(p_born,unwgtfun)
@@ -3426,11 +3834,9 @@
 c Do not include this contribution for final-state branchings
          collrem_xi=0.d0
          collrem_lxi=0.d0
-         if(doreweight)then
-           wgtdegrem_xi=0.d0
-           wgtdegrem_lxi=0.d0
-           wgtdegrem_muF=0.d0
-         endif
+         wgtdegrem_xi=0.d0
+         wgtdegrem_lxi=0.d0
+         wgtdegrem_muF=0.d0
          return
       endif
 
@@ -3439,11 +3845,9 @@
          write (*,*) "No born momenta in sreal_deg"
          collrem_xi=0.d0
          collrem_lxi=0.d0
-         if(doreweight)then
-           wgtdegrem_xi=0.d0
-           wgtdegrem_lxi=0.d0
-           wgtdegrem_muF=0.d0
-         endif
+         wgtdegrem_xi=0.d0
+         wgtdegrem_lxi=0.d0
+         wgtdegrem_muF=0.d0
          return
       endif
 
@@ -3483,13 +3887,11 @@
       collrem_xi=oo2pi * born_wgt * collrem_xi * xnorm
       collrem_lxi=oo2pi * born_wgt * collrem_lxi * xnorm
 
-      if(doreweight)then
-        wgtdegrem_xi=ap*log(shat*delta_used/(2*QES2)) -
+      wgtdegrem_xi=ap*log(shat*delta_used/(2*QES2)) -
      #               apprime - xkkern 
-        wgtdegrem_xi=oo2pi * born_wgt * wgtdegrem_xi * xnorm
-        wgtdegrem_lxi=collrem_lxi
-        wgtdegrem_muF= - oo2pi * born_wgt * ap * xnorm
-      endif
+      wgtdegrem_xi=oo2pi * born_wgt * wgtdegrem_xi * xnorm
+      wgtdegrem_lxi=collrem_lxi
+      wgtdegrem_muF= - oo2pi * born_wgt * ap * xnorm
 
       return
       end
@@ -4572,39 +4974,36 @@
 
  549     continue
 
-         if(doreweight)then
-           wgtwnstmpmuf=0.d0
-           if(abrv.ne.'born' .and. abrv.ne.'grid')then
-             if(abrv(1:2).eq.'vi')then
+         wgtwnstmpmuf=0.d0
+         if(abrv.ne.'born' .and. abrv.ne.'grid')then
+            if(abrv(1:2).eq.'vi')then
                wgtwnstmpmur=0.d0
-             else
+            else
                do i=1,nincoming
-                 if (particle_type(i).ne.1)then
-                   if (particle_type(i).eq.8) then
-                     aj=0
-                   elseif(abs(particle_type(i)).eq.3) then
-                     aj=1
-                   endif
-                   wgtwnstmpmuf=wgtwnstmpmuf-
+                  if (particle_type(i).ne.1)then
+                     if (particle_type(i).eq.8) then
+                        aj=0
+                     elseif(abs(particle_type(i)).eq.3) then
+                        aj=1
+                     endif
+                     wgtwnstmpmuf=wgtwnstmpmuf-
      #                   ( gamma(aj)+2d0*c(aj)*dlog(xicut_used) )
-                 endif
+                  endif
                enddo
                wgtwnstmpmuf=ao2pi*wgtwnstmpmuf*dble(wgt1(1))
                wgtwnstmpmur=2*pi*(beta0*wgtbpower
-     #      +ren_group_coeff*wgtcpower)*ao2pi*dble(wgt1(1))
-             endif
+     #         +ren_group_coeff*wgtcpower)*ao2pi*dble(wgt1(1))
+            endif
 c bsv_wgt here always contains the Born; must subtract it, since 
 c we need the pure NLO terms only
-             wgtnstmp=bsv_wgt+virt_wgt-born_wgt-
+            wgtnstmp=bsv_wgt+virt_wgt-born_wgt-
      #                wgtwnstmpmuf*log(q2fact(1)/QES2)-
      #                wgtwnstmpmur*log(scale**2/QES2)
-           else
-             wgtnstmp=0d0
-             wgtwnstmpmur=0.d0
-           endif
+         else
+            wgtnstmp=0d0
+            wgtwnstmpmur=0.d0
          endif
 
-
          if (abrv(1:2).eq.'vi') then
             bsv_wgt=bsv_wgt-born_wgt
 

=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
--- Template/NLO/SubProcesses/genps_fks.f	2014-03-17 09:43:37 +0000
+++ Template/NLO/SubProcesses/genps_fks.f	2015-02-04 13:27:21 +0000
@@ -93,6 +93,8 @@
       common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
       double precision p_born(0:3,nexternal-1)
       common/pborn/p_born
+      double precision p_ev(0:3,nexternal)
+      common/pev/p_ev
       double precision xi_i_fks_ev,y_ij_fks_ev
       double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
       common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
@@ -267,7 +269,12 @@
       else
 c No PDFs (also use fixed energy when performing tests)
          call compute_tau_y_epem(j_fks,one_body,fksmass,stot,
-     &                          tau_born,ycm_born,ycmhat)
+     &        tau_born,ycm_born,ycmhat)
+         if (j_fks.le.nincoming .and. .not.(softtest.or.colltest)) then
+            write (*,*) 'Process has incoming j_fks, but fixed shat: '/
+     &           /'not allowed for processes generated at NLO.'
+            stop 1
+         endif
       endif
 c Compute Bjorken x's from tau and y
       xbjrk_born(1)=sqrt(tau_born)*exp(ycm_born)
@@ -523,6 +530,7 @@
          do i=1,nexternal
             do j=0,3
                p(j,i)=xp(j,i)
+               p_ev(j,i)=xp(j,i)
             enddo
          enddo
          jac=xjac

=== modified file 'Template/NLO/SubProcesses/handling_lhe_events.f'
--- Template/NLO/SubProcesses/handling_lhe_events.f	2014-10-31 07:52:06 +0000
+++ Template/NLO/SubProcesses/handling_lhe_events.f	2015-02-04 13:27:21 +0000
@@ -483,7 +483,7 @@
       integer event_id
       common /c_event_id/ event_id
       integer i_process
-      common/c_addwrite/i_process
+      common/c_i_process/i_process
       integer nattr,npNLO,npLO
       common/event_attributes/nattr,npNLO,npLO
       include 'reweight_all.inc'
@@ -616,9 +616,6 @@
                  do i=1,mexternal
                     write(ifile,405)(wgtkin_all(j,i,1,iFKS),j=0,3)
                  enddo
-c$$$                 do i=1,mexternal
-c$$$                    write(ifile,405)(wgtkin_all(j,i,2,iFKS),j=0,3)
-c$$$                 enddo
                  write(ifile,402)
      &                wgtxbj_all(1,1,iFKS),wgtxbj_all(2,1,iFKS),
      &                wgtxbj_all(1,2,iFKS),wgtxbj_all(2,2,iFKS),
@@ -746,7 +743,7 @@
       integer ii,j,nps,nng,iFKS,idwgt
       double precision wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
       integer i_process
-      common/c_addwrite/i_process
+      common/c_i_process/i_process
       integer nattr,npNLO,npLO
       common/event_attributes/nattr,npNLO,npLO
       include 'reweight_all.inc'
@@ -976,7 +973,7 @@
       integer ii,j,nps,nng,iFKS,idwgt
       double precision wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
       integer i_process
-      common/c_addwrite/i_process
+      common/c_i_process/i_process
       integer nattr,npNLO,npLO
       common/event_attributes/nattr,npNLO,npLO
       include 'reweight_all.inc'

=== modified file 'Template/NLO/SubProcesses/madfks_plot.f'
--- Template/NLO/SubProcesses/madfks_plot.f	2014-06-27 12:41:14 +0000
+++ Template/NLO/SubProcesses/madfks_plot.f	2015-02-04 13:27:21 +0000
@@ -131,7 +131,7 @@
       end
 
 
-      subroutine outfun(pp,ybst_til_tolab,www,itype)
+      subroutine outfun(pp,ybst_til_tolab,www,iPDG,itype)
 C
 C *WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING*
 C
@@ -156,11 +156,11 @@
       include 'genps.inc'
       include 'reweight.inc'
       include 'reweightNLO.inc'
-      double precision pp(0:3,nexternal),ybst_til_tolab,www
+      double precision pp(0:3,nexternal),ybst_til_tolab
       integer itype
       double precision p(0:4,nexternal),pplab(0:3,nexternal),chybst
      $     ,shybst,chybstmo
-      integer i,j,ibody
+      integer i,j,ibody,i_wgt
       double precision xd(3)
       data (xd(i),i=1,3) /0d0,0d0,1d0/
       integer istatus(nexternal),iPDG(nexternal)
@@ -173,7 +173,7 @@
       common /c_leshouche_inc/idup,mothup,icolup
       integer nwgt,max_weight
       parameter (max_weight=maxscales*maxscales+maxpdfs+1)
-      double precision wgts(max_weight),wgtden,ratio
+      double precision www(max_weight),wgtden,ratio
       double precision xsecScale_acc(maxscales,maxscales)
      $     ,xsecPDFr_acc(0:maxPDFs)
       common /scale_pdf_print/xsecScale_acc,xsecPDFr_acc
@@ -183,13 +183,10 @@
 c Born, n-body or (n+1)-body contribution:
       if(itype.eq.11) then
          ibody=1 ! (n+1)-body
-         wgtden=wgtrefNLO11
       elseif(itype.eq.12)then
          ibody=2 ! n-body
-         wgtden=wgtrefNLO12
       elseif(itype.eq.20)then
          ibody=3 ! Born
-         wgtden=wgtrefNLO20
       else
          write(*,*)'Error in outfun: unknown itype',itype
          stop
@@ -199,8 +196,7 @@
       shybst=sinh(ybst_til_tolab)
       chybstmo=chybst-1.d0
       do i=3,nexternal
-        call boostwdir2(chybst,shybst,chybstmo,xd,
-     #                  pp(0,i),pplab(0,i))
+         call boostwdir2(chybst,shybst,chybstmo,xd,pp(0,i),pplab(0,i))
       enddo
 c Fill the arrays (momenta, status and PDG):
       do i=1,nexternal
@@ -213,57 +209,27 @@
             p(j,i)=pplab(j,i)
          enddo
          p(4,i)=pmass(i)
-         ipdg(i)=idup(i,1)
       enddo
-c The weights comming from reweighting:
-      nwgt=1
-      wgts(1)=www
-      if (wgtden.eq.0d0) then
-         if (www.eq.0d0) then
-            ratio=0d0
-         else
-            if (doreweight) then
-               write (*,*) 'ERROR in madfks_plot.f', wgtden,www
-            else
-               ratio=1d0
-            endif
-         endif
-      else
-c this ratio should essentially be the weight from vegas
-         ratio=www/wgtden
-      endif
-      if (do_rwgt_scale) then
-         do i=1,numscales
-            do j=1,numscales
-               nwgt=nwgt+1
-               wgts(nwgt)=wgtNLOxsecmu(ibody,i,j)*ratio
-            enddo
-         enddo
-      endif
-      if (do_rwgt_pdf) then
-         do i=1,numPDFs-1 ! exclude the central set, so 'numPDFs-1'
-            nwgt=nwgt+1
-            wgts(nwgt)=wgtNLOxsecPDF(ibody,i)*ratio
-         enddo
-      endif
       if(iappl.ne.0)then
          appl_itype     = ibody
-         appl_www_histo = wgts(1)
+         appl_www_histo = www(1)
       endif
-      call analysis_fill(p,istatus,ipdg,wgts,ibody)
+      call analysis_fill(p,istatus,ipdg,www,ibody)
 c Fill the accumulated results
+      i_wgt=1
       if (do_rwgt_scale) then
          do i=1,numscales
             do j=1,numscales
-               xsecScale_acc(i,j)=xsecScale_acc(i,j)+
-     &              wgtNLOxsecmu(ibody,i,j)*ratio
+               i_wgt=i_wgt+1
+               xsecScale_acc(i,j)=xsecScale_acc(i,j)+www(i_wgt)
             enddo
          enddo
       endif
       if (do_rwgt_pdf) then
-         xsecPDFr_acc(0)=xsecPDFr_acc(0)+www
+         xsecPDFr_acc(0)=xsecPDFr_acc(0)+www(1)
          do i=1,numPDFs-1
-            xsecPDFr_acc(i)=xsecPDFr_acc(i)+wgtNLOxsecPDF(ibody,i)*ratio
+            i_wgt=i_wgt+1
+            xsecPDFr_acc(i)=xsecPDFr_acc(i)+www(i_wgt)
          enddo
       endif
  999  return      

=== modified file 'Template/NLO/SubProcesses/reweight_xsec.f'
--- Template/NLO/SubProcesses/reweight_xsec.f	2014-09-24 17:48:39 +0000
+++ Template/NLO/SubProcesses/reweight_xsec.f	2015-02-04 13:27:21 +0000
@@ -975,7 +975,7 @@
       DOUBLE PRECISION       CONV
       PARAMETER (CONV=389379660D0)  !CONV TO PICOBARNS             
       integer i_process
-      common/c_addwrite/i_process
+      common/c_i_process/i_process
 
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
@@ -1206,7 +1206,7 @@
       DOUBLE PRECISION       CONV
       PARAMETER (CONV=389379660D0)  !CONV TO PICOBARNS             
       integer i_process
-      common/c_addwrite/i_process
+      common/c_i_process/i_process
 
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
@@ -1511,7 +1511,7 @@
       DOUBLE PRECISION       CONV
       PARAMETER (CONV=389379660D0)  !CONV TO PICOBARNS             
       integer i_process
-      common/c_addwrite/i_process
+      common/c_i_process/i_process
 
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS

=== modified file 'Template/NLO/SubProcesses/reweight_xsec_events.f'
--- Template/NLO/SubProcesses/reweight_xsec_events.f	2014-05-13 07:36:22 +0000
+++ Template/NLO/SubProcesses/reweight_xsec_events.f	2015-02-04 13:27:21 +0000
@@ -49,7 +49,7 @@
       external compute_rwgt_wgt_Sev,compute_rwgt_wgt_Sev_nbody
      &     ,compute_rwgt_wgt_Hev
       integer i_process
-      common/c_addwrite/i_process
+      common/c_i_process/i_process
 c
       call setrun                !Sets up run parameters
 

=== modified file 'Template/NLO/SubProcesses/setcuts.f'
--- Template/NLO/SubProcesses/setcuts.f	2014-09-02 12:20:56 +0000
+++ Template/NLO/SubProcesses/setcuts.f	2015-02-04 13:27:21 +0000
@@ -489,6 +489,7 @@
                sum_all_s=0d0
                do i=t_channel,-(nexternal-3),-1
 c Breit-wigner can never go on-shell:
+                  if (itree(2,i).gt.0) cycle
                   if ( pmass(itree(2,i),iconfig).gt.sqrt(stot) .and.
      $                 pwidth(itree(2,i),iconfig).gt.0d0) then
                      cBW_FKS(iFKS,itree(2,i))=2
@@ -499,6 +500,7 @@
                if (sum_all_s.gt.sqrt(stot)) then
 c     conflicting BWs: set all s-channels as conflicting
                   do i=t_channel,-(nexternal-3),-1
+                     if (itree(2,i).gt.0) cycle
                      if (cBW_FKS(iFKS,itree(2,i)).ne.2) then
                         cBW_FKS(iFKS,itree(2,i))=1
                         cBW_FKS_mass(iFKS,itree(2,i),-1)=sqrt(stot)/2d0

=== modified file 'Template/NLO/SubProcesses/timing_variables.inc'
--- Template/NLO/SubProcesses/timing_variables.inc	2013-12-20 07:28:29 +0000
+++ Template/NLO/SubProcesses/timing_variables.inc	2015-02-04 13:27:21 +0000
@@ -1,5 +1,10 @@
+*     -*-fortran-*-
+
       real tBefore, tAfter
-      real tOLP, tFastJet, tPDF,
-     & tGenPS, tDSigI, tDSigR
-      common/timings/tOLP, tFastJet, tPDF,
-     & tGenPS, tDSigI, tDSigR
+      real tOLP, tFastJet, tPDF, tGenPS, tDSigI, tDSigR
+      common/timings/tOLP, tFastJet, tPDF, tGenPS, tDSigI, tDSigR
+
+      real tBorn,tIS,tReal,tCount,tFxFx,tf_nb,tf_all,t_as,tr_s,tr_pdf
+     &     ,t_plot,t_cuts
+      common/timings_fNLO/tBorn,tIS,tReal,tCount,tFxFx,tf_nb,tf_all,t_as
+     &     ,tr_s,tr_pdf,t_plot,t_cuts

=== modified file 'Template/NLO/SubProcesses/write_event.f'
--- Template/NLO/SubProcesses/write_event.f	2014-10-31 07:52:06 +0000
+++ Template/NLO/SubProcesses/write_event.f	2015-02-04 13:27:21 +0000
@@ -268,16 +268,7 @@
             write(*,*)doreweight,iwgtinfo
             stop
           endif
-          if (ickkw.eq.4) then
-             if (iwgtinfo.ne.5) then
-                write (*,*)'if using ickkw=4, iwgtinfo should be 5'
-                stop
-             else
-                kwgtinfo=15
-             endif
-          else
-             kwgtinfo=iwgtinfo
-          endif
+          kwgtinfo=iwgtinfo
           write(buff,201)'#aMCatNLO',iSorH_lhe,ifks_lhe(nFKSprocess)
      &         ,jfks_lhe(nFKSprocess),fksfather_lhe(nFKSprocess)
      &         ,ipartner_lhe(nFKSprocess),scale1_lhe(nFKSprocess)

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2014-11-06 16:26:15 +0000
+++ UpdateNotes.txt	2015-02-04 13:27:21 +0000
@@ -1,6 +1,29 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
-2.2.2(06/11/14) OM: Correct a bug in the integration grid (introduces in 2.1.2). This was biasing the cross-section of 
+<<<<<<< TREE
+2.2.2(06/11/14) OM: Correct a bug in the integration grid (introduces in 2.1.2). This was biasing the cross-section of 
+=======
+2.2.3(XX/XX/XX) RF: If a NAN is found, the code now skips that PS point and continues.
+                OM: Fix a bug in MadWeight (correlated param_card was not creating the correct input file)
+                RF: For fNLO runs the virtuals were included twice in the setting of the integration grids. 
+		RF: When requiring more than 1M events for (N)LO+PS runs, do not go to higher precision than 0.001
+                    for the grids and cross section (can be overwritten with the req_acc run_card parameter).
+		RF: Make sure that reweight info (for PDF and scale uncertainties) also works for UNLOPS events.
+		RF: When setting the B's stable in the shower_card, also set the eta_b (PDG=551) stable.
+		OM: Change the Breit-Wigner splitting for the multi-channel integration, use the bwcutoff instead of
+                    the hardcoded value 5.
+		RF: Re-factoring of the structure of the code for fNLO computations.
+                MZ: Fix to bug 1406000 (segfault appearing when doing FxFx merging). Thanks to Josh Bendavid for
+                    having reported it
+                MZ: Fix to a bug occurring when generating event in the "split" mode: the required output was 
+                    not correctly specified
+                MZ: Fixes to different small bugs / improvement in the error and warning messages
+                OM: The built-in pdf "nn23lo" and "nn23lo1" where associate to the wrong lhapdfid in the lhe file
+                    This was creating bias in using SysCalc. (Thanks Alexis)
+                OM: Fix a bug in the re-weighing which was removing the SysCalc weight from the lhe file (thanks Shin-Shan)
+
+2.2.2(06/11/14) OM: Correct a bug in the integration grid (introduces in 2.1.2). This was biasing the cross-section of 
+>>>>>>> MERGE-SOURCE
                     processes like a a > mu+ mu- in the Effective Photon Approximation by three order of magnitude.
                     For LHC processes no sizeable effect have been observe so far.
                 MZ: some informations for aMC@NLO runs which were before passed via include files are

=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py	2014-09-08 11:03:40 +0000
+++ aloha/aloha_writers.py	2015-02-04 13:27:21 +0000
@@ -654,8 +654,10 @@
         """Formatting the variable name to Fortran format"""
         
         if isinstance(name, aloha_lib.ExtVariable):
-            # external parameter nothing to do
+            # external parameter nothing to do but handling model prefix
             self.has_model_parameter = True
+            if name.lower() in ['pi', 'as', 'mu_r', 'aewm1','g']:
+                return name
             return '%s%s' % (aloha.aloha_prefix, name)
         
         if '_' in name:
@@ -727,8 +729,8 @@
         for name in keys:
             fct, objs = self.routine.fct[name]
 
-	    format = ' %s = %s\n' % (name, self.get_fct_format(fct))
-	    try:
+            format = ' %s = %s\n' % (name, self.get_fct_format(fct))
+            try:
                 text = format % ','.join([self.write_obj(obj) for obj in objs])
             except TypeError:
                 text = format % tuple([self.write_obj(obj) for obj in objs])
@@ -1827,6 +1829,32 @@
         name = re.sub('(?P<var>\w*)_(?P<num>\d+)$', self.shift_indices , name)
         
         return name
+
+    def get_fct_format(self, fct):
+        """Put the function in the correct format"""
+        if not hasattr(self, 'fct_format'):
+            one = self.change_number_format(1)
+            self.fct_format = {'csc' : '{0}/cmath.cos(%s)'.format(one),
+                   'sec': '{0}/cmath.sin(%s)'.format(one),
+                   'acsc': 'cmath.asin({0}/(%s))'.format(one),
+                   'asec': 'cmath.acos({0}/(%s))'.format(one),
+                   're': ' complex(%s).real',
+                   'im': 'complex(%s).imag',
+                   'cmath.sqrt': 'cmath.sqrt(%s)',
+                   'sqrt': 'cmath.sqrt(%s)',
+                   'pow': 'pow(%s, %s)',
+                   'complexconjugate': 'complex(%s).conjugate()',
+                   '/' : '{0}/%s'.format(one),
+                   'abs': 'cmath.fabs(%s)'
+                   }
+            
+        if fct in self.fct_format:
+            return self.fct_format[fct]
+        elif hasattr(cmath, fct):
+            self.declaration.add(('fct', fct))
+            return 'cmath.{0}(%s)'.format(fct)
+        else:
+            raise Exception, "Unable to handle function name %s (no special rule defined and not in cmath)" % fct
     
     def define_expression(self):
         """Define the functions in a 100% way """
@@ -1837,6 +1865,30 @@
             for name,obj in self.routine.contracted.items():
                 out.write('    %s = %s\n' % (name, self.write_obj(obj)))
 
+        def sort_fct(a, b):
+            if len(a) < len(b):
+                return -1
+            elif len(a) > len(b):
+                return 1
+            elif a < b:
+                return -1
+            else:
+                return +1
+            
+        keys = self.routine.fct.keys()        
+        keys.sort(sort_fct)
+        for name in keys:
+            fct, objs = self.routine.fct[name]
+            format = '    %s = %s\n' % (name, self.get_fct_format(fct))
+            try:
+                text = format % ','.join([self.write_obj(obj) for obj in objs])
+            except TypeError:
+                text = format % tuple([self.write_obj(obj) for obj in objs])
+            finally:
+                out.write(text)
+
+
+
         numerator = self.routine.expr
         if not 'Coup(1)' in self.routine.infostr:
             coup_name = 'COUP'
@@ -1892,7 +1944,7 @@
             name = self.name
            
         out = StringIO()
-        
+        out.write("import cmath\n")
         if self.mode == 'mg5':
             out.write('import aloha.template_files.wavefunctions as wavefunctions\n')
         else:

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2014-08-07 21:10:49 +0000
+++ madgraph/core/base_objects.py	2015-02-04 13:27:21 +0000
@@ -1513,7 +1513,8 @@
                 if str(part.get('width')) in one_change:
                     part.set('width', rep_pattern.sub(replace, str(part.get('width'))))  
                 if  hasattr(part, 'partial_widths'):
-                    for key, value in part.partial_widths.items():                        part.partial_widths[key] = rep_pattern.sub(replace, value)
+                    for key, value in part.partial_widths.items():    
+                        part.partial_widths[key] = rep_pattern.sub(replace, value)
                 
         #ensure that the particle_dict is up-to-date
         self['particle_dict'] =''

=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
--- madgraph/interface/amcatnlo_run_interface.py	2014-11-05 17:34:19 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2015-02-04 13:27:21 +0000
@@ -140,6 +140,8 @@
         compiler = options['fortran_compiler']
     elif misc.which('gfortran'):
         compiler = 'gfortran'
+    else: 
+        compiler = ''
         
     if 'gfortran' not in compiler:
         if block:
@@ -1487,6 +1489,9 @@
                 raise aMCatNLOError('Required accuracy ("req_acc" in the run_card) should '\
                                         'be between larger than 0 and smaller than 1, '\
                                         'or set to -1 for automatic determination. Current value is %s' % req_acc)
+# For more than 1M events, set req_acc to 0.001 (except when it was explicitly set in the run_card)
+            elif float(req_acc) < 0 and nevents > 1000000 :
+                req_acc='0.001'
 
             shower_list = ['HERWIG6', 'HERWIGPP', 'PYTHIA6Q', 'PYTHIA6PT', 'PYTHIA8']
 
@@ -1518,9 +1523,7 @@
                     nevents_unweighted = []
 
                 split = i == 2 and \
-                        int(self.run_card['nevt_job']) > 0 and \
-                        any([int(l.split()[1]) > int(self.run_card['nevt_job']) \
-                            for l in nevents_unweighted if l])
+                        int(self.run_card['nevt_job']) > 0 
 
                 if i == 2 or not options['only_generation']:
                     # if the number of events requested is zero,
@@ -1787,11 +1790,11 @@
                           '\n      Total cross-section: %(xsect)8.3e +- %(errt)6.1e pb' % \
                         self.cross_sect_dict
 
-                if int(self.run_card['nevents'])>=10000 and self.run_card['reweight_scale']=='.true.' and int(self.run_card['ickkw']) != 4:
+                if int(self.run_card['nevents'])>=10000 and self.run_card['reweight_scale']=='.true.':
                    message = message + \
                        ('\n      Ren. and fac. scale uncertainty: +%0.1f%% -%0.1f%%') % \
                        (scale_pdf_info['scale_upp'], scale_pdf_info['scale_low'])
-                if int(self.run_card['nevents'])>=10000 and self.run_card['reweight_PDF']=='.true.' and int(self.run_card['ickkw']) != 4:
+                if int(self.run_card['nevents'])>=10000 and self.run_card['reweight_PDF']=='.true.':
                    message = message + \
                        ('\n      PDF uncertainty: +%0.1f%% -%0.1f%%') % \
                        (scale_pdf_info['pdf_upp'], scale_pdf_info['pdf_low'])
@@ -2287,7 +2290,7 @@
         Event dir. Return the name of the event file created
         """
         scale_pdf_info={}
-        if (self.run_card['reweight_scale'] == '.true.' or self.run_card['reweight_PDF'] == '.true.') and int(self.run_card['ickkw']) != 4 :
+        if (self.run_card['reweight_scale'] == '.true.' or self.run_card['reweight_PDF'] == '.true.') :
             scale_pdf_info = self.run_reweight(options['reweightonly'])
 
         self.update_status('Collecting events', level='parton', update_results=True)
@@ -2398,6 +2401,18 @@
                     '\n'.join(fjwrapper_lines) + '\n')
 
         extrapaths = self.shower_card['extrapaths'].split()
+
+        # check that the path needed by HW++ and PY8 are set if one uses these shower
+        if shower in ['HERWIGPP', 'PYTHIA8']:
+            path_dict = {'HERWIGPP': ['hepmc_path',
+                                      'thepeg_path',
+                                      'hwpp_path'],
+                         'PYTHIA8': ['pythia8_path']}
+
+            if not all([self.options[ppath] for ppath in path_dict[shower]]):
+                raise aMCatNLOError('Some paths are missing in the configuration file.\n' + \
+                        ('Please make sure you have set these variables: %s' % ', '.join(path_dict[shower])))
+
         if shower == 'HERWIGPP':
             extrapaths.append(pjoin(self.options['hepmc_path'], 'lib'))
 
@@ -3088,8 +3103,13 @@
         scale_upp=0.0
         scale_low=0.0
         if numofscales>0:
-            scale_pdf_info['scale_upp'] = (max(scales)/cntrl_val-1)*100
-            scale_pdf_info['scale_low'] = (1-min(scales)/cntrl_val)*100
+            if cntrl_val != 0.0:
+                scale_pdf_info['scale_upp'] = (max(scales)/cntrl_val-1)*100
+                scale_pdf_info['scale_low'] = (1-min(scales)/cntrl_val)*100
+            else:
+                scale_pdf_info['scale_upp'] = 0.0
+                scale_pdf_info['scale_low'] = 0.0
+
 
         # get the pdf uncertainty in percent (according to the Hessian method)
         lhaid=int(self.run_card['lhaid'])
@@ -3101,15 +3121,23 @@
                 for i in range(int(numofpdf/2)):
                     pdf_upp=pdf_upp+math.pow(max(0.0,pdfs[2*i+1]-cntrl_val,pdfs[2*i+2]-cntrl_val),2)
                     pdf_low=pdf_low+math.pow(max(0.0,cntrl_val-pdfs[2*i+1],cntrl_val-pdfs[2*i+2]),2)
-                scale_pdf_info['pdf_upp'] = math.sqrt(pdf_upp)/cntrl_val*100
-                scale_pdf_info['pdf_low'] = math.sqrt(pdf_low)/cntrl_val*100
+                if cntrl_val != 0.0:
+                    scale_pdf_info['pdf_upp'] = math.sqrt(pdf_upp)/cntrl_val*100
+                    scale_pdf_info['pdf_low'] = math.sqrt(pdf_low)/cntrl_val*100
+                else:
+                    scale_pdf_info['pdf_upp'] = 0.0
+                    scale_pdf_info['pdf_low'] = 0.0
+
         else:
             # use Gaussian method (NNPDF)
             pdf_stdev=0.0
             for i in range(int(numofpdf-1)):
                 pdf_stdev = pdf_stdev + pow(pdfs[i+1] - cntrl_val,2)
             pdf_stdev = math.sqrt(pdf_stdev/int(numofpdf-2))
-            scale_pdf_info['pdf_upp'] = pdf_stdev/cntrl_val*100
+            if cntrl_val != 0.0:
+                scale_pdf_info['pdf_upp'] = pdf_stdev/cntrl_val*100
+            else:
+                scale_pdf_info['pdf_upp'] = 0.0
             scale_pdf_info['pdf_low'] = scale_pdf_info['pdf_upp']
         return scale_pdf_info
 
@@ -3286,8 +3314,8 @@
             if shower == 'PYTHIA8':
                 input_files.append(pjoin(cwd, 'Pythia8.exe'))
                 input_files.append(pjoin(cwd, 'Pythia8.cmd'))
-                input_files.append(pjoin(cwd, 'config.sh'))
                 if os.path.exists(pjoin(self.options['pythia8_path'], 'xmldoc')):
+                    input_files.append(pjoin(cwd, 'config.sh'))
                     input_files.append(pjoin(self.options['pythia8_path'], 'xmldoc'))
                 else:
                     input_files.append(pjoin(self.options['pythia8_path'], 'share/Pythia8/xmldoc'))
@@ -3437,7 +3465,10 @@
                     keep_fourth_arg = True
                     # this is for the split event generation
                     output_files.append('G%s%s_%s' % (args[1], i, args[3]))
-                required_output.append('%s/log_MINT%s.txt' % (current,args[2]))
+                    required_output.append('G%s%s_%s/log_MINT%s.txt' % (args[1], i, args[3],args[2]))
+
+                else:
+                    required_output.append('%s/log_MINT%s.txt' % (current,args[2]))
                 if args[2] in ['0','1']:
                     required_output.append('%s/results.dat' % current)
                 if args[2] == '1':
@@ -3601,21 +3632,23 @@
             for code in ['applgrid','amcfast']:
                 try:
                     p = subprocess.Popen([self.options[code], '--version'], \
-                                             stdout=subprocess.PIPE, stderr=subprocess.PIPE)
-                    output, error = p.communicate()
+                                          stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+                except OSError:
+                    raise aMCatNLOError(('No valid %s installation found. \n' + \
+                       'Please set the path to %s-config by using \n' + \
+                       'MG5_aMC> set <absolute-path-to-%s>/bin/%s-config \n') % (code,code,code,code))
+                else:
+                    output, _ = p.communicate()
                     if code is 'applgrid' and output < '1.4.63':
                         raise aMCatNLOError('Version of APPLgrid is too old. Use 1.4.69 or later.'\
-                                                +' You are using %s',output)
+                                             +' You are using %s',output)
                     if code is 'amcfast' and output < '1.1.1':
                         raise aMCatNLOError('Version of aMCfast is too old. Use 1.1.1 or later.'\
-                                                +' You are using %s',output)
-                except Exception:
-                    raise aMCatNLOError(('No valid %s installation found. \n' + \
-                          'Please set the path to %s-config by using \n' + \
-                          'MG5_aMC> set <absolute-path-to-%s>/bin/%s-config \n') % (code,code,code,code))
+                                             +' You are using %s',output)
+                
             # set-up the Source/make_opts with the correct applgrid-config file
-            appllibs="  APPLLIBS=$(shell %s --ldcflags) $(shell %s --ldflags) \n" \
-                             % (self.options['applgrid'],self.options['amcfast'])
+            appllibs="  APPLLIBS=$(shell %s --ldflags) $(shell %s --ldcflags) \n" \
+                             % (self.options['amcfast'],self.options['applgrid'])
             text=open(pjoin(self.me_dir,'Source','make_opts'),'r').readlines()
             text_out=[]
             for line in text:

=== modified file 'madgraph/interface/common_run_interface.py'
--- madgraph/interface/common_run_interface.py	2014-11-05 17:34:19 +0000
+++ madgraph/interface/common_run_interface.py	2015-02-04 13:27:21 +0000
@@ -731,6 +731,7 @@
            delphes_card.dat
            delphes_trigger.dat
            shower_card.dat [aMCatNLO]
+           FO_analyse_card.dat [aMCatNLO]
            madspin_card.dat [MS]
            transfer_card.dat [MW]
            madweight_card.dat [MW]
@@ -740,7 +741,7 @@
         if text == '':
             logger.warning('File %s is empty' % path)
             return 'unknown'
-        text = re.findall('(<MGVersion>|ParticlePropagator|<mg5proccard>|CEN_max_tracker|#TRIGGER CARD|parameter set name|muon eta coverage|QES_over_ref|MSTP|b_stable|MSTU|Begin Minpts|gridpack|ebeam1|block\s+mw_run|BLOCK|DECAY|launch|madspin|transfer_card\.dat|set)', text, re.I)
+        text = re.findall('(<MGVersion>|ParticlePropagator|<mg5proccard>|CEN_max_tracker|#TRIGGER CARD|parameter set name|muon eta coverage|QES_over_ref|MSTP|b_stable|FO_ANALYSIS_FORMAT|MSTU|Begin Minpts|gridpack|ebeam1|block\s+mw_run|BLOCK|DECAY|launch|madspin|transfer_card\.dat|set)', text, re.I)
         text = [t.lower() for t in text]
         if '<mgversion>' in text or '<mg5proccard>' in text:
             return 'banner'
@@ -769,6 +770,8 @@
             return 'param_card.dat'
         elif 'b_stable' in text:
             return 'shower_card.dat'
+        elif 'fo_analysis_format' in text:
+            return 'FO_analyse_card.dat'
         elif 'decay' in text and 'launch' in text and 'madspin' in text:
             return 'madspin_card.dat'
         elif 'launch' in text and 'set' in text:
@@ -1701,8 +1704,9 @@
         """ find a valid run_name for the current job """
 
         name = 'run_%02d'
-        data = [int(s[4:]) for s in os.listdir(pjoin(me_dir,'Events')) if
-                        s.startswith('run_') and len(s)>5 and s[4:].isdigit()]
+        data = [int(s[4:j]) for s in os.listdir(pjoin(me_dir,'Events')) for 
+                j in range(4,len(s)+1) if \
+                s.startswith('run_') and s[4:j].isdigit()]
         return name % (max(data+[0])+1)
 
 

=== modified file 'madgraph/interface/madevent_interface.py'
--- madgraph/interface/madevent_interface.py	2014-10-20 14:26:15 +0000
+++ madgraph/interface/madevent_interface.py	2015-02-04 13:27:21 +0000
@@ -1610,6 +1610,7 @@
             model = self.find_model_name()
             process = self.process # define in find_model_name
             self.results = gen_crossxhtml.AllResults(model, process, self.me_dir)
+            self.results.resetall(self.me_dir)
         self.results.def_web_mode(self.web)
         
         self.prompt = "%s>"%os.path.basename(pjoin(self.me_dir))
@@ -2508,7 +2509,15 @@
             os.remove(pjoin(self.me_dir,'SubProcesses', 'combine.log'))
         except Exception:
             pass
-        self.cluster.launch_and_wait('../bin/internal/run_combine', 
+        
+        if self.options['run_mode'] ==1 and self.options['cluster_tmp_path']:
+            tmpcluster = cluster.MultiCore(nb_core=1)
+            tmpcluster.launch_and_wait('../bin/internal/run_combine', 
+                                        cwd=pjoin(self.me_dir,'SubProcesses'),
+                                        stdout=pjoin(self.me_dir,'SubProcesses', 'combine.log'),
+                                        required_output=[pjoin(self.me_dir,'SubProcesses', 'combine.log')])
+        else:
+            self.cluster.launch_and_wait('../bin/internal/run_combine', 
                                         cwd=pjoin(self.me_dir,'SubProcesses'),
                                         stdout=pjoin(self.me_dir,'SubProcesses', 'combine.log'),
                                         required_output=[pjoin(self.me_dir,'SubProcesses', 'combine.log')])
@@ -2573,6 +2582,8 @@
         self.check_combine_events(args)
         self.update_status('Storing parton level results', level='parton')
 
+
+
         run = self.run_name
         tag = self.run_card['run_tag']
         devnull = open(os.devnull, 'w')
@@ -2588,36 +2599,54 @@
         files.cp(input, output) 
 
         # 2) Treat the files present in the P directory
-        for P_path in SubProcesses.get_subP(self.me_dir):
-            G_dir = [G for G in os.listdir(P_path) if G.startswith('G') and 
-                                                os.path.isdir(pjoin(P_path,G))]
-            for G in G_dir:
-                G_path = pjoin(P_path,G)
-                # Remove events file (if present)
-                if os.path.exists(pjoin(G_path, 'events.lhe')):
-                    os.remove(pjoin(G_path, 'events.lhe'))
-                # Store results.dat
-                if os.path.exists(pjoin(G_path, 'results.dat')):
-                    input = pjoin(G_path, 'results.dat')
-                    output = pjoin(G_path, '%s_results.dat' % run)
-                    files.cp(input, output) 
-                # Store log
-                if os.path.exists(pjoin(G_path, 'log.txt')):
-                    input = pjoin(G_path, 'log.txt')
-                    output = pjoin(G_path, '%s_log.txt' % run)
-                    files.mv(input, output) 
-                # Grid
-                for name in ['ftn26']:
-                    if os.path.exists(pjoin(G_path, name)):
-                        if os.path.exists(pjoin(G_path, '%s_%s.gz'%(run,name))):
-                            os.remove(pjoin(G_path, '%s_%s.gz'%(run,name)))
-                        input = pjoin(G_path, name)
-                        output = pjoin(G_path, '%s_%s' % (run,name))
-                        files.mv(input, output)
-                        misc.gzip(pjoin(G_path, output), error=None) 
-                # Delete ftn25 to ensure reproducible runs
-                if os.path.exists(pjoin(G_path, 'ftn25')):
-                    os.remove(pjoin(G_path, 'ftn25'))
+        # Ensure that the number of events is different of 0 
+        if self.results.current['nb_event'] == 0:
+            logger.warning("No event detected. No cleaning performed! This should allow to run:\n" +
+                           "    cd Subprocesses; ../bin/internal/combine_events\n"+
+                           "  to have your events if those one are missing.")
+        else:
+            for P_path in SubProcesses.get_subP(self.me_dir):
+                G_dir = [G for G in os.listdir(P_path) if G.startswith('G') and 
+                                                    os.path.isdir(pjoin(P_path,G))]
+                for G in G_dir:
+                    G_path = pjoin(P_path,G)
+                    try:
+                        # Remove events file (if present)
+                        if os.path.exists(pjoin(G_path, 'events.lhe')):
+                            os.remove(pjoin(G_path, 'events.lhe'))
+                    except Exception:
+                        continue
+                    try:
+                        # Store results.dat
+                        if os.path.exists(pjoin(G_path, 'results.dat')):
+                            input = pjoin(G_path, 'results.dat')
+                            output = pjoin(G_path, '%s_results.dat' % run)
+                            files.cp(input, output) 
+                    except Exception:
+                        continue                    
+                    # Store log
+                    try:
+                        if os.path.exists(pjoin(G_path, 'log.txt')):
+                            input = pjoin(G_path, 'log.txt')
+                            output = pjoin(G_path, '%s_log.txt' % run)
+                            files.mv(input, output) 
+                    except Exception:
+                        continue
+                    try:   
+                        # Grid
+                        for name in ['ftn26']:
+                            if os.path.exists(pjoin(G_path, name)):
+                                if os.path.exists(pjoin(G_path, '%s_%s.gz'%(run,name))):
+                                    os.remove(pjoin(G_path, '%s_%s.gz'%(run,name)))
+                                input = pjoin(G_path, name)
+                                output = pjoin(G_path, '%s_%s' % (run,name))
+                                files.mv(input, output)
+                                misc.gzip(pjoin(G_path, output), error=None) 
+                    except Exception:
+                        continue
+                    # Delete ftn25 to ensure reproducible runs
+                    if os.path.exists(pjoin(G_path, 'ftn25')):
+                        os.remove(pjoin(G_path, 'ftn25'))
 
         # 3) Update the index.html
         misc.call(['%s/gen_cardhtml-pl' % self.dirbin],

=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2014-11-05 17:34:19 +0000
+++ madgraph/interface/madgraph_interface.py	2015-02-04 13:27:21 +0000
@@ -899,13 +899,14 @@
             return
 
         # request that we have one or two > in the process
-        if process.count('>') not in [1,2]:
+        nbsep = len(re.findall('>\D', process)) # not use process.count because of QCD^2>2
+        if nbsep not in [1,2]:
             raise self.InvalidCmd(
                'wrong format for \"%s\" this part requires one or two symbols \'>\', %s found'
-               % (process, process.count('>')))
+               % (process, nbsep))
 
         # we need at least one particles in each pieces
-        particles_parts = process.split('>')
+        particles_parts = re.split('>\D', process)
         for particles in particles_parts:
             if re.match(r'^\s*$', particles):
                 raise self.InvalidCmd(
@@ -3494,7 +3495,7 @@
         a ProcessDefinition."""
 
         # Check basic validity of the line
-        if not line.count('>') in [1,2]:
+        if not len(re.findall('>\D', line)) in [1,2]:
             self.do_help('generate')
             raise self.InvalidCmd('Wrong use of \">\" special character.')
 
@@ -5018,8 +5019,8 @@
                     continue
                 path = self.options[key]
                 #this is for pythia8
-                if key == 'pythia8_path' and not os.path.isfile(pjoin(MG5DIR, path, 'include', 'Pythia.h')):
-                    if not os.path.isfile(pjoin(path, 'include', 'Pythia.h')):
+                if key == 'pythia8_path' and not os.path.isfile(pjoin(MG5DIR, path, 'include', 'Pythia8', 'Pythia.h')):
+                    if not os.path.isfile(pjoin(path, 'include', 'Pythia8', 'Pythia.h')):
                         self.options['pythia8_path'] = None
                     else:
                         continue
@@ -6062,8 +6063,8 @@
                 logger.info('The value for lhapdf in the current configuration does not ' + \
                         'correspond to a valid executable.\nPlease set it correctly either in ' + \
                         'input/mg5_configuration or with "set lhapdf /path/to/lhapdf-config" ' + \
-                        'and regenrate the process. \nTo avoid regeneration, manually edit the ' + \
-                        ('%s/Source/fj_lhapdf_opts file.\n' % self._export_dir ) + \
+                        'and regenrate the process. \nTo avoid regeneration, edit the ' + \
+                        ('%s/Cards/amcatnlo_configuration.txt file.\n' % self._export_dir ) + \
                         'Note that you can still compile and run aMC@NLO with the built-in PDFs\n')
 
             compiler_dict = {'fortran': self.options['fortran_compiler'],
@@ -6213,9 +6214,12 @@
                     tmp_mass = mass
                     for p in mode:
                         try:
-                            tmp_mass -= abs(eval(str(p.mass), data))
+                            value_mass = eval(str(p.mass), data)
                         except Exception:
-                            tmp_mass -= abs(eval("mdl_"+str(p.mass), data))
+                            # the p object can still be UFO reference. since the 
+                            # mass name might hve change load back the MG5 one.
+                            value_mass = eval(str(model.get_particle(p.pdg_code).get('mass')), data)
+                        tmp_mass -= abs(value_mass)             
                     if tmp_mass <=0:
                         continue
     

=== modified file 'madgraph/iolibs/export_fks.py'
--- madgraph/iolibs/export_fks.py	2014-11-04 11:30:53 +0000
+++ madgraph/iolibs/export_fks.py	2015-02-04 13:27:21 +0000
@@ -562,6 +562,7 @@
                      'fks_Sij.f',
                      'fks_powers.inc',
                      'fks_singular.f',
+                     'c_weight.inc',
                      'fks_inc_chooser.f',
                      'leshouche_inc_chooser.f',
                      'configs_and_props_inc_chooser.f',
@@ -630,7 +631,7 @@
 
         #import nexternal/leshouches in Source
         ln('nexternal.inc', '../../Source', log=False)
-        ln('leshouche_decl.inc', '../../Source', log=False)
+        ln('born_leshouche.inc', '../../Source', log=False)
 
 
         # Return to SubProcesses dir
@@ -2599,19 +2600,34 @@
                 for initial_state in init_states:
                     if initial_state in pdf_codes.keys():
                         if subproc_group:
-                            pdf_lines = pdf_lines + \
-                                        ("%s%d=PDG2PDF(ABS(LPP(IB(%d))),%d*LP," + \
+                            if abs(pdgtopdf[initial_state]) <= 7:  
+                                pdf_lines = pdf_lines + \
+                                    ("%s%d=PDG2PDF(ABS(LPP(IB(%d))),%d*LP," + \
                                          "XBK(IB(%d)),DSQRT(Q2FACT(%d)))\n") % \
                                          (pdf_codes[initial_state],
                                           i + 1, ibeam, pdgtopdf[initial_state],
                                           ibeam, ibeam)
+                            else:
+                                # setting other partons flavours outside quark, gluon, photon to be 0d0
+                                pdf_lines = pdf_lines + \
+                                    ("c settings other partons flavours outside quark, gluon, photon to 0d0\n" + \
+                                     "%s%d=0d0\n") % \
+                                         (pdf_codes[initial_state],i + 1)                                
                         else:
-                            pdf_lines = pdf_lines + \
-                                        ("%s%d=PDG2PDF(ABS(LPP(%d)),%d*LP," + \
+                            if abs(pdgtopdf[initial_state]) <= 7:  
+                                pdf_lines = pdf_lines + \
+                                    ("%s%d=PDG2PDF(ABS(LPP(%d)),%d*LP," + \
                                          "XBK(%d),DSQRT(Q2FACT(%d)))\n") % \
                                          (pdf_codes[initial_state],
                                           i + 1, ibeam, pdgtopdf[initial_state],
                                           ibeam, ibeam)
+                            else:
+                                # setting other partons flavours outside quark, gluon, photon to be 0d0
+                                pdf_lines = pdf_lines + \
+                                    ("c settings other partons flavours outside quark, gluon, photon to 0d0\n" + \
+                                     "%s%d=0d0\n") % \
+                                         (pdf_codes[initial_state],i + 1)                                
+
                 pdf_lines = pdf_lines + "ENDIF\n"
 
             # Add up PDFs for the different initial state particles

=== modified file 'madgraph/iolibs/template_files/madevent_symmetry.f'
--- madgraph/iolibs/template_files/madevent_symmetry.f	2013-12-09 02:09:49 +0000
+++ madgraph/iolibs/template_files/madevent_symmetry.f	2015-02-04 13:27:21 +0000
@@ -588,14 +588,12 @@
          xwidth(-i)=prwidth(-i,iconfig)
          if (xwidth(-i) .gt. 0d0) then
             nbw=nbw+1
-            if (iarray(nbw) .eq. 1) then
-               if(xmass(-i).gt.prmass(-i,iconfig)+5d0*xwidth(-i)) then
-                  failConfig=.true.
-                  return
-               else
-                  xmass(-i)=max(xmass(-i),prmass(-i,iconfig)-5d0*xwidth(-i))
-               endif
-            else if(forcebw(-i) .eq. 1) then
+            if(forcebw(-i) .eq. 1) then
+c               if (iarray(nbw) .ne. 1) then
+c                  write(*,*) "fail due to iarray", iarray(nbw)
+c                  failConfig=.true.
+c                  return
+c               endif
                if(xmass(-i).gt.prmass(-i,iconfig)+bwcutoff*xwidth(-i)) then
                   failConfig=.true.
                   return
@@ -603,6 +601,14 @@
                   xmass(-i)=max(xmass(-i),prmass(-i,iconfig)-
      $                 bwcutoff*xwidth(-i))
                endif