Merge lp:~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig_merged2.3.3 into lp:~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig

[marco zaro]

[Olivier Mattelaer (olivier-mattelaer)]

Hi Marco,

Can you comment on this change?
For example, you made a couple of change to MadSpin, why?

Thanks,

Olivier

[marco zaro (marco-zaro)]

Hi Olivier,
the changes in MS should be due to the merge with the main 2.3.3 or 2.3.4 trunk…
I am trying to reorganise a bit the EW branches and in particular to remove the ones which are not used any more…
Cheers,

Marco
On 18 Jan 2016, at 12:07, Olivier Mattelaer <email address hidden> wrote:

> Hi Marco,
>
> Can you comment on this change?
> For example, you made a couple of change to MadSpin, why?
>
>
> Thanks,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig_merged2.3.3/+merge/282921
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig_merged2.3.3 into lp:~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig.

[lp:~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig_merged2.3.3]

766. By Marco Zaro <mzaro@ingrid-ui1> on 2016-01-20

ct14PDF w QED efefcts should be included now (Thanks to Davide Pagani!)

767. By Marco Zaro <mzaro@ingrid-ui1> on 2016-01-20

added other PDF sets

[Olivier Mattelaer (olivier-mattelaer)]

Ok then,

I guess this is pointless to review this right.
We will do this review when the merge with the stable version occurs.

Except if you want a review on a specific point.

Olivier

[lp:~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig_merged2.3.3]

768. By Marco Zaro <email address hidden> on 2016-01-22

fixes for test_ME

769. By Marco Zaro <email address hidden> on 2016-01-22

fixed a bug affecting symmetry factors for extra cnt's with final state splittings

770. By Marco Zaro <mzaro@ingrid-ui1> on 2016-01-25

another bug fixed (related to chaching and extra_cnt's)

[Olivier Mattelaer (olivier-mattelaer)]

Then I guess go ahead and this will be review in a future merge no?
Otherwise this sounds like a lot of work to review stuff that need to be review at different stage anyway.

Don't you agree?

Olivier

[marco zaro (marco-zaro)]

Ciao Olivier,
anyway, this is not for 2.3.4…
Cheers,

Marco
On 27 Jan 2016, at 12:08, Olivier Mattelaer <email address hidden> wrote:

> Then I guess go ahead and this will be review in a future merge no?
> Otherwise this sounds like a lot of work to review stuff that need to be review at different stage anyway.
>
> Don't you agree?
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig_merged2.3.3/+merge/282921
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig_merged2.3.3 into lp:~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig.

[lp:~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig_merged2.3.3]

771. By Rikkert Frederix on 2016-02-02

Make sure the code doesn't crash when there is a channel with NO
points passing cuts.

772. By Marco Zaro <email address hidden> on 2016-02-04

fix number of points in test_ME (soft)

773. By Rikkert Frederix on 2016-02-05

small fix in the printing of the limits in the fort.77 file in the case
some of the Collinear tests are failing.

774. By Rikkert Frederix on 2016-02-05

fixed a problem with the split-orders collinear counter terms (final
state radiation).

775. By marco zaro on 2016-02-09

fix in the run_interface to have the stats in the full summary

776. By marco zaro on 2016-02-10

implementing the splitorder stuff in applgrid

777. By marco zaro on 2016-02-10

fix to previous push

778. By marco zaro on 2016-02-10

another fix for applgrid

779. By Rikkert Frederix on 2016-02-11

Also set all the split orders to zero if there is an unstable,
non-rescued PS point in the virtual corrections.

780. By Rikkert Frederix on 2016-02-11

Made the checking of the accuracy in the BinothLHA order dependent.

781. By marco zaro on 2016-02-11

fix in fks_singular

782. By marco zaro on 2016-02-11

fix to rev 767 (some files have not been updated)

783. By Rikkert Frederix on 2016-02-12

fixed a problem in the BinothLHA which was related to the skipping of
the unstable points introduced in revision 779.

784. By Marco Zaro <mzaro@ingrid-ui1> on 2016-02-22

BUG the behaviour when extra CNTs exist should be fixed.
Also a bug in find_reals_to_integrate has been found and fixed

785. By marco zaro on 2016-02-23

fix in the generation of the extra counterterms (the amplitude was not
generated)

786. By Marco Zaro <mzaro@ingrid-ui1> on 2016-02-25

fix stupid bug in common_run_interfact leading to KeyError: mg5_path on ingrid

787. By Rikkert Frederix on 2016-02-25

Better updating of the fraction of PS points to compute the virtual
corrections. Removes some run-away cases where the virtual_fraction
became one even though there was no sign of an largely fluctuating
virtual.

788. By marco zaro on 2016-03-09

res_*.txt now includes the dir by dir info

789. By marco zaro on 2016-03-09

the 6fto5f implementation in fks_singular seems to work (at least for
the total rate)

790. By marco zaro on 2016-03-11

fix error for plots in folders with 0 xsec

791. By marco zaro on 2016-03-18

fix in 6to5fcounterterm

792. By marco-zaro <email address hidden> on 2016-04-05

fix in one test inside fks_singular.f (no effect on results)

793. By marco zaro on 2016-04-13

      changes in BinothLHA.f to be able (in principle) to use MC/hel

794. By Valentin Hirschi on 2016-04-14

1. Modified MadLoop so that UV counterterms are added in an order which
   is as canonical as possible, so that Identical loop ME which differ
   only in their massless flavors can be merged together.

795. By Valentin Hirschi on 2016-04-15

1. Improved the loop diagram generation algorithm so that it picks an
   order of l-cut particles making it likely that loop MatrixElement which
   are equivalent but with different flavors are indeed flagged as identical
   when compared.

796. By Valentin Hirschi on 2016-05-03

1. Fixed some virtual 'equal' test issue when doing mixed QCD_QED corrections.

797. By marco zaro on 2016-05-12

changes related to event generation, still not working

798. By marco zaro on 2016-05-24

merged with 2.4.1 rev 410

799. By marco zaro on 2016-05-24

merged with 2.4.1 rev415

800. By marco zaro on 2016-05-24

merged with FKS_EW_flattened_dsig_merged2.3.3 r799

801. By marco zaro on 2016-05-24

fix in fks_singular 6to5cnt;
fix in loop_interface ( perturbations other than QCD were forbidden)
fix in test_extra_ew

802. By marco zaro on 2016-05-24

fixing unit tests

803. By marco zaro on 2016-05-25

fixes towards having event generation working

804. By Valentin Hirschi on 2016-05-26

1. Fixed the missing color specification of the interactions of the hardcoded loop model in the test 'test_loop_diagram_generation.py'

805. By Valentin Hirschi on 2016-05-26

1. Fixed the other issue in the test_loop_diagram_generation.py which was that the EW model was not loaded in Feynman gauge for some reason.

806. By marco zaro on 2016-05-26

fix in histograms.py to correctly do the 68% for NNPDFQED

807. By marco zaro on 2016-05-30

fixing unit tests

808. By marco zaro on 2016-06-23

fix in
  Template/NLO/SubProcesses/reweight0.inc

809. By marco zaro on 2016-11-16

merged with 2.5.1 r342

810. By marco zaro on 2016-11-22

fix in common_run_interface

[Olivier Mattelaer (olivier-mattelaer)]

just closing this review since they are no point to keep it open (right?)

[marco zaro (marco-zaro)]

yes, close it
cheers,
M
On 07 Aug 2017, at 11:04, Olivier Mattelaer <email address hidden> wrote:

> just closing this review since they are no point to keep it open (right?)
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig_merged2.3.3/+merge/282921
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig_merged2.3.3 into lp:~maddevelopers/mg5amcnlo/FKS_EW_flattened_dsig.

Unmerged revisions

810. By marco zaro on 2016-11-22

fix in common_run_interface

809. By marco zaro on 2016-11-16

merged with 2.5.1 r342

808. By marco zaro on 2016-06-23

fix in
  Template/NLO/SubProcesses/reweight0.inc

807. By marco zaro on 2016-05-30

fixing unit tests

806. By marco zaro on 2016-05-26

fix in histograms.py to correctly do the 68% for NNPDFQED

805. By Valentin Hirschi on 2016-05-26

1. Fixed the other issue in the test_loop_diagram_generation.py which was that the EW model was not loaded in Feynman gauge for some reason.

804. By Valentin Hirschi on 2016-05-26

1. Fixed the missing color specification of the interactions of the hardcoded loop model in the test 'test_loop_diagram_generation.py'

803. By marco zaro on 2016-05-25

fixes towards having event generation working

802. By marco zaro on 2016-05-24

fixing unit tests

801. By marco zaro on 2016-05-24

fix in fks_singular 6to5cnt;
fix in loop_interface ( perturbations other than QCD were forbidden)
fix in test_extra_ew

=== modified file 'MadSpin/__init__.py' (properties changed: -x to +x)
=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2015-05-22 16:18:38 +0000
+++ MadSpin/decay.py	2016-11-23 08:16:42 +0000
@@ -64,7 +64,7 @@
 from madgraph import MG5DIR, MadGraph5Error
 import madgraph.various.misc as misc
 #import time
-import tests.parallel_tests.test_aloha as test_aloha
+
 
 class MadSpinError(MadGraph5Error):
     pass
@@ -111,16 +111,16 @@
             # change the wgt associate to the additional weight
             start, stop = self.rwgt.find('<rwgt>'), self.rwgt.find('</rwgt>')
             if start != -1 != stop :
-                pattern = re.compile(r'''<\s*wgt id=\'(?P<id>[^\']+)\'\s*>\s*(?P<val>[\ded+-.]*)\s*</wgt>''')
+                pattern = re.compile(r'''<\s*wgt id=[\'\"](?P<id>[^\'\"]+)[\'\"]\s*>\s*(?P<val>[\ded+-.]*)\s*</wgt>''')
                 data = pattern.findall(self.rwgt)
+                if len(data)==0:
+                    print self.rwgt
                 try:
                     text = ''.join('   <wgt id=\'%s\'> %+15.7e </wgt>\n' % (pid, float(value) * factor)
                                      for (pid,value) in data) 
                 except ValueError, error:
                     raise Exception, 'Event File has unvalid weight. %s' % error
-                self.rwgt = self.rwgt[:start] + '<rwgt>\n'+ text + self.rwgt[stop:]
-            
-            
+                self.rwgt = self.rwgt[:start] + '<rwgt>\n'+ text + self.rwgt[stop:]          
 
     def string_event_compact(self):
         """ return a string with the momenta of the event written 
@@ -906,7 +906,7 @@
             if init[0] in self.banner.param_card['decay'].decay_table:
                 br *= self.banner.param_card['decay'].decay_table[init[0]].get(lhaid).value
                 br *= self.get_br(decay)
-            else:
+            elif -init[0] in self.banner.param_card['decay'].decay_table:
                 init = -init[0]
                 lhaid=[x if self.model.get_particle(x)['self_antipart'] else -x
                        for x in final]
@@ -914,6 +914,13 @@
                 lhaid = tuple([len(final)] + lhaid)
                 br *= self.banner.param_card['decay'].decay_table[init].get(lhaid).value
                 br *= self.get_br(decay)
+            elif init[0] not in self.decay_ids and -init[0] not in self.decay_ids:
+                logger.warning("No Branching ratio applied for %s. Please check if this is expected" % init[0])
+                br *= self.get_br(decay)
+            else:
+                raise MadGraph5Error,"No valid decay for %s. No 2 body decay for that particle. (three body are not supported by MadSpin)" % init[0]
+
+                
 
         for decays in ids.values():
             if len(decays) == 1:
@@ -929,7 +936,7 @@
                     except ValueError:
                         break
                 br /= math.factorial(nb)
-
+                
         return br
 
             
@@ -1013,7 +1020,7 @@
         self[tag]['decaying'] = tuple(decaying)
                 
         # sanity check
-        assert self[tag]['total_br'] <= 1.01, self[tag]['total_br']
+        assert self[tag]['total_br'] <= 1.01, "wrong BR for %s: %s " % (tag,self[tag]['total_br'])
 
 
         
@@ -1217,6 +1224,7 @@
                                   'P%s' % me.get('processes')[0].shell_string())
             topo['decays'] = []
             topo['total_br'] = 0 
+            topo['Pid'] = proc.get('id')
         
             if skip:
                 self.add_me_symmetric(skip, topo)
@@ -1629,7 +1637,7 @@
                 (' '.join([str(i) for i in opts['particles']]),
                  ' '.join('--%s=%s' % (key,value) for (key,value) in opts.items()
                         if key not in ['model', 'force', 'particles'] and value))
-        cmd.exec_cmd('import model %s' % model.get('name'))
+        cmd.exec_cmd('import model %s' % model.get('modelpath+restriction'))
 #        cmd._curr_model = model
 #        cmd._curr_fortran_model = helas_call_writers.FortranUFOHelasCallWriter(model)
         cmd.exec_cmd(line)
@@ -1742,7 +1750,7 @@
 
         return finalfound
    
-         
+        
     def reorder_branch(self,branch):
         """ branch is a string with the definition of a decay chain
                 If branch contains " A > B C , B > ... " 
@@ -1757,19 +1765,31 @@
             if list_branch[index]==' ' or list_branch[index]=='': del list_branch[index]
         #print list_branch
         for index, item in enumerate(list_branch):
-            if item =="," and list_branch[index+1]!="(": 
-                if list_branch[index-2]==list_branch[index+1]:
-                    # swap the two particles before the comma:
-                    temp=list_branch[index-2]
-                    list_branch[index-2]=list_branch[index-1]
-                    list_branch[index-1]=temp
-            if item =="," and list_branch[index+1]=="(":
-                if list_branch[index-2]==list_branch[index+2]:
-                    # swap the two particles before the comma:
-                    temp=list_branch[index-2]
-                    list_branch[index-2]=list_branch[index-1]
-                    list_branch[index-1]=temp
 
+            if item[-1] =="," and list_branch[index+1]!="(":
+                # search pos of B and C 
+                counter=1
+                while 1:
+                  if list_branch[index-counter].find("=")<0:
+                     break
+                  counter+=1
+                if list_branch[index-counter-1]==list_branch[index+1]:
+                    # swap the two particles before the comma:
+                    temp=list_branch[index-counter-1]
+                    list_branch[index-counter-1]=list_branch[index-counter]
+                    list_branch[index-counter]=temp
+            if item[-1] =="," and list_branch[index+1]=="(":
+                # search pos of B and C 
+                counter=1
+                while 1:
+                  if list_branch[index-counter].find("=")<0:
+                     break
+                  counter+=1
+                if list_branch[index-counter -1]==list_branch[index+2]:
+                    # swap the two particles before the comma:
+                    temp=list_branch[index-counter-1]
+                    list_branch[index-counter-1]=list_branch[index-counter]
+                    list_branch[index-counter]=temp
 
         new_branch=""
         for item in list_branch:
@@ -1777,7 +1797,6 @@
 
         return new_branch, list_branch[0]
 
-
     def set_light_parton_massless(self,topo):
         """ masses of light partons are set to zero for 
             the evaluation of the matrix elements
@@ -2097,7 +2116,8 @@
         # launch the decay and reweighting
         self.mscmd.update_status('MadSpin: Decaying Events')
         efficiency = self.decaying_events(self.inverted_decay_mapping)
-        if  efficiency != 1:
+        self.efficiency = efficiency
+        if  efficiency != 1 and any(v==-1 for v in self.br_per_id.values()):
             # need to change the banner information [nb_event/cross section]
             files.cp(self.outputfile.name, '%s_tmp' % self.outputfile.name)
             self.outputfile = open(self.outputfile.name, 'w')
@@ -2398,7 +2418,13 @@
                     self.curr_event.particle[part_for_curr_evt]['momentum']=ext_mom[prod2full[part-1]-1]
                     self.curr_event.particle[part_for_curr_evt]['helicity']=helicities[prod2full[part-1]-1]
                     if not use_mc_masses or abs(pid) not in self.MC_masses:
-                        self.curr_event.particle[part_for_curr_evt]['mass']=self.banner.get('param_card','mass', abs(pid)).value
+                        try:
+                            self.curr_event.particle[part_for_curr_evt]['mass']=self.banner.get('param_card','mass', abs(pid)).value
+                        except KeyError:
+                            if self.model.get_particle(abs(pid)).get('mass').lower() == 'zero':
+                                self.curr_event.particle[part_for_curr_evt]['mass'] = 0
+                            else:
+                                raise
                     else:
                         self.curr_event.particle[part_for_curr_evt]['mass']=self.MC_masses[abs(pid)]
 
@@ -2589,7 +2615,7 @@
     
 
     @misc.mute_logger()
-    @test_aloha.set_global()
+    @misc.set_global()
     def generate_all_matrix_element(self):
         """generate the full series of matrix element needed by Madspin.
         i.e. the undecayed and the decay one. And associate those to the 
@@ -2627,14 +2653,23 @@
                       if re.search('^\s*add\s+process', line)]
         
         mgcmd = self.mgcmd
-        modelpath = self.model.get('modelpath')
-        if os.path.basename(modelpath) != mgcmd._curr_model['name']:
-            name, restrict = mgcmd._curr_model['name'].rsplit('-',1)
-            if os.path.exists(pjoin(os.path.dirname(modelpath),name, 'restrict_%s.dat' % restrict)):
-                modelpath = pjoin(os.path.dirname(modelpath), mgcmd._curr_model['name'])
+        modelpath = self.model.get('modelpath+restriction')
+
+        commandline="import model %s" % modelpath
+        if not self.model.mg5_name:
+            commandline += ' --modelname'
             
-        commandline="import model %s " % modelpath
         mgcmd.exec_cmd(commandline)
+        # Handle the multiparticle of the banner        
+        #for name, definition in self.mscmd.multiparticles:
+        if hasattr(self.mscmd, 'multiparticles_ms'):
+            for name, pdgs in  self.mscmd.multiparticles_ms.items():
+                if name == 'all':
+                    continue
+                #self.banner.get('proc_card').get('multiparticles'):
+                mgcmd.do_define("%s = %s" % (name, ' '.join(`i` for i in pdgs)))
+            
+
         mgcmd.exec_cmd("set group_subprocesses False")
 
         logger.info('generating the production square matrix element')
@@ -2651,9 +2686,13 @@
                 proc_nb = '@ %i' % proc_nb 
             else:
                 proc_nb = ''
+                        
+            if ',' in proc:
+                raise MadSpinError, 'MadSpin can not decay event which comes from a decay chain.'+\
+                        '\n  The full decay chain should either be handle by MadGraph or by Masdspin.'
             
             if '[' not in proc:
-                commandline+="add process %s  --no_warning=duplicate;" % proc
+                commandline+="add process %s %s  --no_warning=duplicate;" % (proc, proc_nb)
             else:
                 process, order, final = re.split('\[\s*(.*)\s*\]', proc)
                 commandline+="add process %s %s --no_warning=duplicate;" % (process, proc_nb)
@@ -2685,6 +2724,7 @@
                         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()))        
@@ -3072,7 +3112,7 @@
 
         
         probe_weight = []
-        
+       
 
         starttime = time.time()
         ev = -1
@@ -3110,6 +3150,7 @@
             atleastonedecay=False
             for decay in self.all_ME[production_tag]['decays']:
                 tag = decay['decay_tag']
+
                 if decay_mapping and not tag in decay_mapping:
                     continue
                 if not tag:
@@ -3123,10 +3164,11 @@
                 else:
                     max_decay[tag] = weight
                     #print weight, max_decay[name]
-                    #raise Exception   
+                    #raise Exception 
+                      
             if not atleastonedecay:
                 # NO decay [one possibility is all decay are identical to their particle]
-                logger.info('No independent decay for one type of final states -> skip those events')
+                logger.info('No independent decay for one type of final states -> skip those events for the maximum weight computation')
                 nb_decay[decaying] = numberev
                 ev += numberev -1
                 continue
@@ -3198,14 +3240,16 @@
                         logger.debug('Decay channel %s :Using maximum weight %s (BR: %s)' % \
                                     (','.join(associated_decay), max_weight, br/nb_finals)) 
 
-        if __debug__: 
+#        if __debug__: 
         # check that all decay have a max_weight and fix it if not the case.
-            for prod in self.all_ME.values():
-                for dec in prod['decays']:
-                    if dec['decay_tag']:                                                
-                        assert 'max_weight' in dec and dec['max_weight'] ,\
-                                  'fail for %s (%s)' % (str(dec['decay_tag']), \
-                                                  os.path.basename(prod['path']))
+        for prod in self.all_ME.values():
+            for dec in prod['decays']:
+                if dec['decay_tag'] and not 'max_weight' in dec:
+                    dec['max_weight'] = 0. 
+
+#                        assert 'max_weight' in dec and dec['max_weight'] ,\
+#                                  'fail for %s (%s)' % (str(dec['decay_tag']), \
+#                                                  os.path.basename(prod['path']))
         self.evtfile.seek(0)
         return
 
@@ -3255,7 +3299,7 @@
 
         return max_weight
         
-    def loadfortran(self, mode, path, stdin_text):
+    def loadfortran(self, mode, path, stdin_text, first=True):
         """ call the fortran executable """
 
         tmpdir = ''
@@ -3272,8 +3316,18 @@
             self.calculator[('full',path,)] = external 
             self.calculator_nbcall[('full',path)] = 1 
 
-        external.stdin.write(stdin_text)
-        
+        try:
+            external.stdin.write(stdin_text)
+        except IOError:
+            if not first:
+                raise
+            try:
+                external.terminate()
+            except:
+                pass
+            del self.calculator[('full',path,)]
+            return self.loadfortran(mode, path, stdin_text, first=False)
+
         if mode == 'maxweight':
             maxweight=float(external.stdout.readline())
             output = maxweight
@@ -3649,6 +3703,45 @@
                             d1colup1=maxcol
                             d2colup2=maxcol
                             d2colup1=0
+                        elif colord1==-3 and colord2==-3 and colormother == 3:
+                            maxcol+=2
+                            d1colup1=0
+                            d1colup2=maxcol
+                            d2colup1=0
+                            d2colup2=maxcol-1
+                        elif (colord1==-3 and colord2==3 and colormother == 3) or\
+                             (colord1==-3 and colord2==3 and colormother == -3):
+                            maxcol+=2
+                            d1colup1 = 0
+                            d1colup2 = maxcol
+                            d2colup1 = maxcol-1
+                            d2colup2 = 0
+                        elif (colord1==3 and colord2==-3 and colormother == 3) or\
+                            (colord1==3 and colord2==-3 and colormother == -3):
+                            maxcol+=2
+                            d1colup1=maxcol
+                            d1colup2=0
+                            d2colup1=0
+                            d2colup2=maxcol-1
+                        elif colord1==3 and colord2==3 and colormother == -3:
+                            maxcol+=2
+                            d1colup1=maxcol
+                            d1colup2=0
+                            d2colup1=maxcol-1
+                            d2colup2=0    
+                        elif colord2==8 and colord1==8 and colormother ==8:
+                            maxcol+=1
+                            ran = random.random()
+                            if ran> 0.5:
+                                d1colup2=colup2
+                                d1colup1=maxcol
+                                d2colup2=maxcol
+                                d2colup1=colup1
+                            else:                            
+                                d1colup2=maxcol
+                                d1colup1=colup1
+                                d2colup2=colup2
+                                d2colup1=maxcol                        
                         else:
                             raise Exception, 'color combination not treated by MadSpin (yet). (%s,%s,%s)' \
                                 % (colord1,colord2,colormother)
@@ -3785,16 +3878,24 @@
         ms_banner = ""
         cross_section = True # tell if possible to write the cross-section in advance
         total_br = []
+        self.br_per_id = {}
         for production in self.all_ME.values():
             one_br = 0
+            partial_br = 0
             for decay in production['decays']:
                 if not decay['decay_tag']:
                     cross_section = False
                     one_br += decay['br']
                     continue
+                partial_br += decay['br']
                 ms_banner += "# %s\n" % ','.join(decay['decay_tag']).replace('\n',' ')
                 ms_banner += "# BR: %s\n# max_weight: %s\n" % (decay['br'], decay['max_weight'])
                 one_br += decay['br']
+            
+            if production['Pid'] not in self.br_per_id:
+                self.br_per_id[production['Pid']] = partial_br
+            elif self.br_per_id[production['Pid']] != partial_br:
+                self.br_per_id[production['Pid']] = -1
             total_br.append(one_br)
         
         if __debug__:
@@ -3802,7 +3903,6 @@
                 assert production['total_br'] - min(total_br) < 1e-4
         
         self.branching_ratio = max(total_br) * eff
-
         #self.banner['madspin'] += ms_banner
         # Update cross-section in the banner
         if 'mggenerationinfo' in self.banner:
@@ -3833,18 +3933,36 @@
                     mg_info[i] = '%s : %s' % (info, value * self.branching_ratio)
                 self.banner['mggenerationinfo'] = '\n'.join(mg_info)
                 
-            
-        
-        if 'init' in self.banner:
+        self.cross = 0
+        self.error = 0
+        if 'init' in self.banner and (eff!=1 or not any(v==-1 for v in self.br_per_id.values())):
             new_init =''
+            curr_proc = 0
+            has_missing=False
             for line in self.banner['init'].split('\n'):
                 if len(line.split()) != 4:
                     new_init += '%s\n' % line
                 else:
+                    curr_proc += 1 
                     data = [float(nb) for nb in line.split()]
-                    data[:3] = [ data[i] * self.branching_ratio for i  in range(3)]
+                    id = int(data[-1])
+                    if id in self.br_per_id and not any(v==-1 for v in self.br_per_id.values()):
+                        data[:3] = [data[i] * self.br_per_id[id] for i  in range(3)]
+                    else:
+                        data[:3] = [ data[i] * self.branching_ratio for i  in range(3)]
+                        has_missing=True
                     new_init += ' %.12E %.12E %.12E %i\n' % tuple(data)
+                    cross, error = [float(d) for d in data[:2]]
+                    self.cross += cross
+                    self.error += error**2
+                    
+                    
             self.banner['init'] = new_init
+            self.error = math.sqrt(self.error)
+            if has_missing and curr_proc not in [0,1]:
+                logger.warning('''The partial cross section for each subprocess can not be determine. due
+    Reason: multiple final state in the same subprocess (and the presence of multiple BR)
+    Consequence: the <init> information of the lhe will therefore be incorrect. Please correct it if needed.''')
         self.banner.write(self.outputfile, close_tag=False)        
         
     def terminate_fortran_executables(self, path_to_decay=0 ):

=== modified file 'MadSpin/interface_madspin.py' (properties changed: -x to +x)
--- MadSpin/interface_madspin.py	2015-08-04 14:28:22 +0000
+++ MadSpin/interface_madspin.py	2016-11-23 08:16:42 +0000
@@ -32,8 +32,10 @@
 import madgraph.interface.extended_cmd as extended_cmd
 import madgraph.interface.madgraph_interface as mg_interface
 import madgraph.interface.master_interface as master_interface
+import madgraph.interface.madevent_interface as madevent_interface
 import madgraph.various.misc as misc
 import madgraph.iolibs.files as files
+import madgraph.iolibs.export_v4 as export_v4
 import madgraph.various.banner as banner
 import madgraph.various.lhe_parser as lhe_parser
 
@@ -176,7 +178,13 @@
         logger.info("process: %s" % process)
         logger.info("options: %s" % option)
 
-
+        if not hasattr(self,'multiparticles_ms'):
+            for key, value in self.banner.get_detail('proc_card','multiparticles'):
+                try:
+                    self.do_define('%s = %s' % (key, value))
+                except self.mg5cmd.InvalidCmd:  
+                    pass
+                
         # Read the final state of the production process:
         #     "_full" means with the complete decay chain syntax 
         #     "_compact" means without the decay chain syntax 
@@ -201,12 +209,13 @@
             mg_names = True
         model_name = self.banner.get('proc_card', 'model')
         if model_name:
+            model_name = os.path.expanduser(model_name)
             self.load_model(model_name, mg_names, complex_mass)
         else:
             raise self.InvalidCmd('Only UFO model can be loaded in MadSpin.')
-        
         # check particle which can be decayed:
         self.final_state = set()
+        final_model = False
         for line in self.banner.proc_card:
             line = ' '.join(line.strip().split())
             if line.startswith('generate'):
@@ -214,30 +223,50 @@
             elif line.startswith('add process'):
                 self.final_state.update(self.mg5cmd.get_final_part(line[11:]))
             elif line.startswith('define'):
-                self.mg5cmd.exec_cmd(line, printcmd=False, precmd=False, postcmd=False)            
+                try:
+                    self.mg5cmd.exec_cmd(line, printcmd=False, precmd=False, postcmd=False)
+                except self.mg5cmd.InvalidCmd:
+                    if final_model:
+                        raise
+                    else:
+                        key = line.split()[1]
+                        if key in self.multiparticles_ms:
+                            del self.multiparticles_ms[key]            
             elif line.startswith('set'):
                 self.mg5cmd.exec_cmd(line, printcmd=False, precmd=False, postcmd=False)
+            elif line.startswith('import model'):
+                if model_name in line:
+                    final_model = True
                     
             
                 
     def import_model(self, args):
         """syntax: import model NAME CARD_PATH
             args didn't include import model"""
-            
+        
         bypass_check = False
         if '--bypass_check' in args:
             args.remove('--bypass_check')
             bypass_check = True
             
-        if len(args) != 2:     
+        if len(args) == 1:  
+            logger.warning("""No param_card defined for the new model. We will use the default one but this might completely wrong.""")
+        elif len(args) != 2:
             return self.InvalidCmd, 'import model requires two arguments'
         
         model_name = args[0]
-        card = args[1]
-        if not os.path.exists(card):
-            raise self.InvalidCmd('%s: no such file' % card)
-
         self.load_model(model_name, False, False)
+        
+        if len(args) == 2:
+            card = args[1]
+            if not os.path.exists(card):
+                raise self.InvalidCmd('%s: no such file' % card)
+        else:
+            card = "madspin_param_card.dat"
+            export_v4.UFO_model_to_mg4.create_param_card_static(self.model,
+                                card, rule_card_path=None)
+
+        
 
         #Check the param_card
         if not bypass_check:
@@ -255,11 +284,15 @@
             if diff:
                 raise self.InvalidCmd('''Original param_card differs on some parameters:
     %s
-    Due to those preferences, we prefer not to proceed. If you are sure about what you are doing, 
+    Due to those differences, we prefer not to proceed. If you are sure about what you are doing, 
     you can use the command \'import model MODELNAME PARAM_CARD_PATH --bypass_check\''''
     % diff.replace('\n','\n    '))
+   
+   
                 
-        #OK load the new param_card
+        #OK load the new param_card (but back up the old one)
+        if 'slha' in self.banner:
+            self.banner['slha_original'] = self.banner['slha']
         self.banner['slha'] = open(card).read()
         if hasattr(self.banner, 'param_card'):
             del self.banner.param_card
@@ -359,12 +392,14 @@
         args = self.split_arg(line)
         self.check_set(args)
         
+        if args[0] not in ['ms_dir', 'run_card']:
+            args[1] = args[1].lower()
+        
         if args[0] in  ['max_weight', 'BW_effect','ms_dir', 'spinmode']:
             self.options[args[0]] = args[1]
             if args[0] == 'ms_dir':
                 self.options['curr_dir'] = self.options['ms_dir']
         elif args[0] == 'seed':
-            import random
             random.seed(int(args[1]))
             self.seed = int(args[1])
         elif args[0] == 'BW_cut':
@@ -426,7 +461,19 @@
     
     def do_define(self, line):
         """ """
-        return self.mg5cmd.do_define(line)
+
+        try:
+            self.mg5cmd.exec_cmd('define %s' % line)
+        except:
+            #cleaning if the error is recover later
+            key = line.split()[0]
+            if hasattr(self, 'multiparticles_ms') and key in self.multiparticles_ms:
+                del self.multiparticles_ms[key]
+            raise
+           
+        self.multiparticles_ms = dict([(k,list(pdgs)) for k, pdgs in \
+                                        self.mg5cmd._multiparticles.items()])
+    
     
     def update_status(self, *args, **opts):
         """ """
@@ -500,7 +547,6 @@
         model_line = self.banner.get('proc_card', 'full_model_line')
 
         if not self.seed:
-            import random
             self.seed = random.randint(0, int(30081*30081))
             self.do_set('seed %s' % self.seed)
             logger.info('Will use seed %s' % self.seed)
@@ -523,6 +569,9 @@
         generate_all.run()
                         
         self.branching_ratio = generate_all.branching_ratio
+        self.cross = generate_all.cross
+        self.error = generate_all.error
+        self.efficiency = generate_all.efficiency
         try:
             self.err_branching_ratio = generate_all.err_branching_ratio
         except Exception:
@@ -619,11 +668,15 @@
             #seed is specified need to use that one:
             open(pjoin(self.options['ms_dir'],'seeds.dat'),'w').write('%s\n'%self.seed)
             #remove all ranmar_state
-            for name in glob.glob(pjoin(self.options['ms_dir'], '*', 'SubProcesses','*','ranmar_state.dat')):
+            for name in misc.glob(pjoin('*', 'SubProcesses','*','ranmar_state.dat'), 
+                                                        self.options['ms_dir']):
                 os.remove(name)    
         
         generate_all.ending_run()
         self.branching_ratio = generate_all.branching_ratio
+        self.cross = generate_all.cross
+        self.error = generate_all.error
+        self.efficiency = generate_all.efficiency
         try:
             self.err_branching_ratio = generate_all.err_branching_ratio
         except Exception:
@@ -648,8 +701,7 @@
         #   of event to generate for each type of particle.
         # 2. Generate the events requested
         # 3. perform the merge of the events.
-        #    if not enough events. re-generate the missing one.
-                
+        #    if not enough events. re-generate the missing one.  
         # First define an utility function for generating events when needed
         def generate_events(pdg, nb_event, mg5, restrict_file=None, cumul=False):
             """generate new events for this particle
@@ -679,24 +731,54 @@
                     else:
                         mg5.exec_cmd("generate %s" % proc)
                         mg5.exec_cmd("output %s -f" % decay_dir)
-                options = dict(mg5.options)
-                import madgraph.interface.madevent_interface as madevent_interface
-                me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
+                    
+                    options = dict(mg5.options)
+                    if self.options['ms_dir']:
+                        misc.sprint("start gridpack!")
+                        # we are in gridpack mode -> create it
+                        me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
                                                 decay_dir), options=options)
-                me5_cmd.options["automatic_html_opening"] = False
-                if self.options["run_card"]:
-                    run_card = self.options["run_card"]
+                        me5_cmd.options["automatic_html_opening"] = False
+                        if self.options["run_card"]:
+                            run_card = self.options["run_card"]
+                        else:
+                            run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))                        
+                        
+                        run_card["iseed"] = self.seed
+                        run_card['gridpack'] = True
+                        run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
+                        param_card = self.banner['slha']
+                        open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
+                        self.seed += 1
+                        # actually creation
+                        me5_cmd.exec_cmd("generate_events run_01 -f")
+                        me5_cmd.exec_cmd("exit")                        
+                        #remove pointless informat
+                        misc.call(["rm", "Cards", "bin", 'Source', 'SubProcesses'], cwd=decay_dir)
+                        misc.call(['tar', '-xzpvf', 'run_01_gridpack.tar.gz'], cwd=decay_dir)
+                
+                # Now generate the events
+
+                if not self.options['ms_dir']:
+                    me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
+                                                    decay_dir), options=mg5.options)
+                    me5_cmd.options["automatic_html_opening"] = False
+                    if self.options["run_card"]:
+                        run_card = self.options["run_card"]
+                    else:
+                        run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))
+                    run_card["nevents"] = int(1.2*nb_event)
+                    run_card["iseed"] = self.seed
+                    run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
+                    param_card = self.banner['slha']
+                    open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
+                    self.seed += 1
+                    me5_cmd.exec_cmd("generate_events run_01 -f")
+                    me5_cmd.exec_cmd("exit")
+                    out[i] = lhe_parser.EventFile(pjoin(decay_dir, "Events", 'run_01', 'unweighted_events.lhe.gz'))            
                 else:
-                    run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))
-                run_card["nevents"] = int(1.2*nb_event)
-                run_card["iseed"] = self.seed
-                run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
-                param_card = self.banner['slha']
-                open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
-                self.seed += 1
-                me5_cmd.exec_cmd("generate_events run_01 -f")
-                me5_cmd.exec_cmd("exit")     
-                out[i] = lhe_parser.EventFile(pjoin(decay_dir, "Events", 'run_01', 'unweighted_events.lhe.gz'))            
+                    misc.call(['run.sh', str(int(1.2*nb_event)), str(self.seed)], cwd=decay_dir)     
+                    out[i] = lhe_parser.EventFile(pjoin(decay_dir, 'events.lhe.gz'))            
                 if cumul:
                     break
             
@@ -713,7 +795,7 @@
                 os.mkdir(self.path_me) 
         else:
             # cleaning
-            for name in glob.glob(pjoin(self.path_me, "decay_*_*")):
+            for name in misc.glob("decay_*_*", self.path_me):
                 shutil.rmtree(name)
 
         self.events_file.close()
@@ -730,7 +812,6 @@
 
         # Handle the banner of the output file
         if not self.seed:
-            import random
             self.seed = random.randint(0, int(30081*30081))
             self.do_set('seed %s' % self.seed)
             logger.info('Will use seed %s' % self.seed)
@@ -749,8 +830,9 @@
         # 2. Generate the events requested
         with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
             mg5 = self.mg5cmd
-            modelpath = self.model.get('modelpath')
-            mg5.exec_cmd("import model %s" % modelpath)      
+            if not self.model:
+                modelpath = self.model.get('modelpath+restriction')
+                mg5.exec_cmd("import model %s" % modelpath)      
             to_event = {}
             for pdg, nb_needed in to_decay.items():
                 #check if a splitting is needed
@@ -808,6 +890,12 @@
             else:
                 raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")
         self.branching_ratio = br
+        self.efficiency = 1
+        self.cross, self.error = self.banner.get_cross(witherror=True)
+        self.cross *= br
+        self.error *= br
+        
+        
         # modify the cross-section in the init block of the banner
         self.banner.scale_init_cross(self.branching_ratio)
                     
@@ -852,10 +940,10 @@
                 else:
                     #need to select the file according to the associate cross-section
                     r = random.random()
-                    tot = sum(events.cross for events in to_event[particle.pdg])
+                    tot = sum(to_event[particle.pdg][key].cross for key in to_event[particle.pdg])
                     r = r * tot
                     cumul = 0
-                    for j,events in enumerate(to_event[particle.pdg]):
+                    for j,events in to_event[particle.pdg].items():
                         cumul += events.cross
                         if r < cumul:
                             decay_file = events
@@ -923,17 +1011,15 @@
         #base_model = import_ufo.import_model(model_path)
 
         # Import model
-        base_model = import_ufo.import_model(name, decay=True)
+        base_model = import_ufo.import_model(name, decay=True,
+                                               complex_mass_scheme=complex_mass)
 
         if use_mg_default:
             base_model.pass_particles_name_in_mg_default()
-        if complex_mass:
-            base_model.change_mass_to_complex_scheme()
         
         self.model = base_model
         self.mg5cmd._curr_model = self.model
         self.mg5cmd.process_model()
-        
 
 if __name__ == '__main__':
     

=== modified file 'MadSpin/src/driver.f'
--- MadSpin/src/driver.f	2015-01-29 14:35:06 +0000
+++ MadSpin/src/driver.f	2016-11-23 08:16:42 +0000
@@ -907,7 +907,7 @@
                 if (index_p2.gt.0) then
                    m2_tchan(i)=m(index_p2)
                 else
-        write(*,*) 'Warning: m_2^2 is negative in t-channel branching '
+        write(*,*) 'Warning: m_2^2 is negative in t-channel branching ',m2_tchan(i) 
                 endif
             endif
          ! extract phi

=== added directory 'PLUGIN'
=== added file 'PLUGIN/__init__.py'
=== modified file 'Template/Common/Cards/delphes_card_ATLAS.dat'
--- Template/Common/Cards/delphes_card_ATLAS.dat	2014-05-13 14:23:18 +0000
+++ Template/Common/Cards/delphes_card_ATLAS.dat	2016-11-23 08:16:42 +0000
@@ -10,7 +10,7 @@
   MuonTrackingEfficiency
 
   ChargedHadronMomentumSmearing
-  ElectronEnergySmearing
+  ElectronMomentumSmearing
   MuonMomentumSmearing
 
   TrackMerger
@@ -20,6 +20,7 @@
   PhotonEfficiency
   PhotonIsolation
 
+  ElectronFilter
   ElectronEfficiency
   ElectronIsolation
 
@@ -28,10 +29,15 @@
 
   MissingET
 
+  NeutrinoFilter
   GenJetFinder
+  GenMissingET
+  
   FastJetFinder
 
   JetEnergyScale
+  
+  JetFlavorAssociation
 
   BTagging
   TauTagging
@@ -75,11 +81,11 @@
   # add EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # tracking efficiency formula for charged hadrons
-  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.60) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.85) + \
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.60) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.85) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
 
@@ -94,13 +100,13 @@
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # tracking efficiency formula for electrons
-  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.73) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e2) * (0.95) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0e2)                * (0.99) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.50) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e2) * (0.83) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e2)                * (0.90) + \
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.73) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e2) * (0.95) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0e2)                * (0.99) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.50) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e2) * (0.83) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e2)                * (0.90) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
 
@@ -115,11 +121,11 @@
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # tracking efficiency formula for muons
-  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) + \
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
 
@@ -134,31 +140,25 @@
   # set ResolutionFormula {resolution formula as a function of eta and pt}
 
   # resolution formula for charged hadrons
-  set ResolutionFormula {                  (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.02) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e1) * (0.01) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0e1 && pt <= 2.0e2) * (0.03) + \
-                                           (abs(eta) <= 1.5) * (pt > 2.0e2)                * (0.05) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.03) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e1) * (0.02) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e1 && pt <= 2.0e2) * (0.04) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 2.0e2)                * (0.05)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
 }
 
-#################################
-# Energy resolution for electrons
-#################################
+###################################
+# Momentum resolution for electrons
+###################################
 
-module EnergySmearing ElectronEnergySmearing {
+module MomentumSmearing ElectronMomentumSmearing {
   set InputArray ElectronTrackingEfficiency/electrons
   set OutputArray electrons
 
   # set ResolutionFormula {resolution formula as a function of eta and energy}
 
-  set ResolutionFormula {                  (abs(eta) <= 2.5) * (energy > 0.1   && energy <= 2.5e1) * (energy*0.015) + \
-                                           (abs(eta) <= 2.5) * (energy > 2.5e1)                    * sqrt(energy^2*0.005^2 + energy*0.05^2 + 0.25^2) + \
-                         (abs(eta) > 2.5 && abs(eta) <= 3.0)                                       * sqrt(energy^2*0.005^2 + energy*0.05^2 + 0.25^2) + \
-                         (abs(eta) > 3.0 && abs(eta) <= 5.0)                                       * sqrt(energy^2*0.107^2 + energy*2.08^2)}
-
+  # resolution formula for electrons
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
 }
 
 ###############################
@@ -170,16 +170,10 @@
   set OutputArray muons
 
   # set ResolutionFormula {resolution formula as a function of eta and pt}
-
   # resolution formula for muons
-  set ResolutionFormula {                  (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.03) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 5.0e1) * (0.03) + \
-                                           (abs(eta) <= 1.5) * (pt > 5.0e1 && pt <= 1.0e2) * (0.04) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0e2)                * (0.07) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.04) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 5.0e1) * (0.04) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 5.0e1 && pt <= 1.0e2) * (0.05) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e2)                * (0.10)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
 }
 
 ##############
@@ -189,7 +183,7 @@
 module Merger TrackMerger {
 # add InputArray InputArray
   add InputArray ChargedHadronMomentumSmearing/chargedHadrons
-  add InputArray ElectronEnergySmearing/electrons
+  add InputArray ElectronMomentumSmearing/electrons
   add InputArray MuonMomentumSmearing/muons
   set OutputArray tracks
 }
@@ -209,6 +203,14 @@
   set EFlowPhotonOutputArray eflowPhotons
   set EFlowNeutralHadronOutputArray eflowNeutralHadrons
 
+  set ECalEnergyMin 0.5
+  set HCalEnergyMin 1.0
+
+  set ECalEnergySignificanceMin 1.0
+  set HCalEnergySignificanceMin 1.0
+
+  set SmearTowerCenter true
+
   set pi [expr {acos(-1)}]
 
   # lists of the edges of each tower in eta and phi
@@ -257,15 +259,15 @@
   # http://arxiv.org/pdf/physics/0608012v1 jinst8_08_s08003
   # http://villaolmo.mib.infn.it/ICATPP9th_2005/Calorimetry/Schram.p.pdf
   # http://www.physics.utoronto.ca/~krieger/procs/ComoProceedings.pdf
-  set ECalResolutionFormula {                  (abs(eta) <= 3.2) * sqrt(energy^2*0.0017^2 + energy*0.101^2) + \
+  set ECalResolutionFormula {                  (abs(eta) <= 3.2) * sqrt(energy^2*0.0017^2 + energy*0.101^2) +
                              (abs(eta) > 3.2 && abs(eta) <= 4.9) * sqrt(energy^2*0.0350^2 + energy*0.285^2)}
 
   # set HCalResolutionFormula {resolution formula as a function of eta and energy}
   # http://arxiv.org/pdf/hep-ex/0004009v1
   # http://villaolmo.mib.infn.it/ICATPP9th_2005/Calorimetry/Schram.p.pdf
-  set HCalResolutionFormula {                  (abs(eta) <= 1.7) * sqrt(energy^2*0.0302^2 + energy*0.5205^2 + 1.59^2) + \
-                             (abs(eta) > 1.7 && abs(eta) <= 3.2) * sqrt(energy^2*0.0500^2 + energy*0.706^2) + \
-                             (abs(eta) > 3.2 && abs(eta) <= 4.9) * sqrt(energy^2*0.9420^2 + energy*0.075^2)}
+  set HCalResolutionFormula {                  (abs(eta) <= 1.7) * sqrt(energy^2*0.0302^2 + energy*0.5205^2 + 1.59^2) +
+                             (abs(eta) > 1.7 && abs(eta) <= 3.2) * sqrt(energy^2*0.0500^2 + energy*0.706^2) +
+                             (abs(eta) > 3.2 && abs(eta) <= 4.9) * sqrt(energy^2*0.09420^2 + energy*1.00^2)}
 }
 
 ####################
@@ -285,15 +287,15 @@
 ###################
 
 module Efficiency PhotonEfficiency {
-  set InputArray Calorimeter/photons
+  set InputArray Calorimeter/eflowPhotons
   set OutputArray photons
 
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # efficiency formula for photons
-  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) + \
+  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) +
                          (abs(eta) > 2.5)                                   * (0.00)}
 }
 
@@ -311,7 +313,19 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
+}
+
+#################
+# Electron filter
+#################
+
+module PdgCodeFilter ElectronFilter {
+  set InputArray Calorimeter/eflowTracks
+  set OutputArray electrons
+  set Invert true
+  add PdgCode {11}
+  add PdgCode {-11}
 }
 
 #####################
@@ -319,15 +333,15 @@
 #####################
 
 module Efficiency ElectronEfficiency {
-  set InputArray ElectronEnergySmearing/electrons
+  set InputArray ElectronFilter/electrons
   set OutputArray electrons
 
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # efficiency formula for electrons
-  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) + \
+  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) +
                          (abs(eta) > 2.5)                                   * (0.00)}
 }
 
@@ -345,7 +359,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
 }
 
 #################
@@ -359,9 +373,9 @@
   # set EfficiencyFormula {efficiency as a function of eta and pt}
 
   # efficiency formula for muons
-  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.7) * (pt > 10.0)  * (0.85) + \
+  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.7) * (pt > 10.0)  * (0.85) +
                          (abs(eta) > 2.7)                                   * (0.00)}
 }
 
@@ -379,7 +393,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.25
 }
 
 ###################
@@ -405,12 +419,33 @@
   set EnergyOutputArray energy
 }
 
+
+#####################
+# Neutrino Filter
+#####################
+
+module PdgCodeFilter NeutrinoFilter {
+
+  set InputArray Delphes/stableParticles
+  set OutputArray filteredParticles
+
+  set PTMin 0.0
+
+  add PdgCode {12}
+  add PdgCode {14}
+  add PdgCode {16}
+  add PdgCode {-12}
+  add PdgCode {-14}
+  add PdgCode {-16}
+
+}
+
 #####################
 # MC truth jet finder
 #####################
 
 module FastJetFinder GenJetFinder {
-  set InputArray Delphes/stableParticles
+  set InputArray NeutrinoFilter/filteredParticles
 
   set OutputArray jets
 
@@ -421,6 +456,19 @@
   set JetPTMin 20.0
 }
 
+
+#########################
+# Gen Missing ET merger
+########################
+
+module Merger GenMissingET {
+# add InputArray InputArray
+  add InputArray NeutrinoFilter/filteredParticles
+  set MomentumOutputArray momentum
+}
+
+
+
 ############
 # Jet finder
 ############
@@ -445,8 +493,25 @@
   set InputArray FastJetFinder/jets
   set OutputArray jets
 
- # scale formula for jets
-  set ScaleFormula {1.00}
+  # scale formula for jets
+  set ScaleFormula {  sqrt( (3.0 - 0.2*(abs(eta)))^2 / pt + 1.0 )  }
+}
+
+########################
+# Jet Flavor Association
+########################
+
+module JetFlavorAssociation JetFlavorAssociation {
+  
+  set PartonInputArray Delphes/partons
+  set ParticleInputArray Delphes/allParticles
+  set ParticleLHEFInputArray Delphes/allParticlesLHEF
+  set JetInputArray JetEnergyScale/jets
+  
+  set DeltaR 0.5
+  set PartonPTMin 1.0
+  set PartonEtaMax 2.5
+
 }
 
 ###########
@@ -454,54 +519,60 @@
 ###########
 
 module BTagging BTagging {
-  set PartonInputArray Delphes/partons
   set JetInputArray JetEnergyScale/jets
 
   set BitNumber 0
 
-  set DeltaR 0.5
-
-  set PartonPTMin 1.0
-
-  set PartonEtaMax 2.5
-
   # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt}
   # PDG code = the highest PDG code of a quark or gluon inside DeltaR cone around jet axis
   # gluon's PDG code has the lowest priority
 
+  # based on ATL-PHYS-PUB-2015-022
+
   # default efficiency formula (misidentification rate)
-  add EfficiencyFormula {0} {0.001}
+  add EfficiencyFormula {0} {0.002+7.3e-06*pt}
 
   # efficiency formula for c-jets (misidentification rate)
-  add EfficiencyFormula {4} {                                      (pt <= 15.0) * (0.000) + \
-                                                (abs(eta) <= 1.2) * (pt > 15.0) * (0.2*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 1.2 && abs(eta) <= 2.5) * (pt > 15.0) * (0.1*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 2.5)                                  * (0.000)}
+  add EfficiencyFormula {4} {0.20*tanh(0.02*pt)*(1/(1+0.0034*pt))}
 
   # efficiency formula for b-jets
-  add EfficiencyFormula {5} {                                      (pt <= 15.0) * (0.000) + \
-                                                (abs(eta) <= 1.2) * (pt > 15.0) * (0.5*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 1.2 && abs(eta) <= 2.5) * (pt > 15.0) * (0.4*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 2.5)                                  * (0.000)}
+  add EfficiencyFormula {5} {0.80*tanh(0.003*pt)*(30/(1+0.086*pt))}
 }
 
-module TauTagging TauTagging {
+#############
+# tau-tagging
+#############
+
+module TrackCountingTauTagging TauTagging {
+ 
   set ParticleInputArray Delphes/allParticles
   set PartonInputArray Delphes/partons
+  set TrackInputArray TrackMerger/tracks
   set JetInputArray JetEnergyScale/jets
 
-  set DeltaR 0.5
+  set DeltaR 0.2
+  set DeltaRTrack 0.2
 
+  set TrackPTMin 1.0
+ 
   set TauPTMin 1.0
-
   set TauEtaMax 2.5
 
-  # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt}
+  # instructions: {n-prongs} {eff} 
+  
+  # 1 - one prong efficiency
+  # 2 - two or more efficiency
+  # -1 - one prong mistag rate
+  # -2 - two or more mistag rate
+ 
+  set BitNumber 0
+ 
+  # taken from ATL-PHYS-PUB-2015-045 (medium working point)
+  add EfficiencyFormula {1} {0.70}
+  add EfficiencyFormula {2} {0.60}
+  add EfficiencyFormula {-1} {0.02}
+  add EfficiencyFormula {-2} {0.01}
 
-  # default efficiency formula (misidentification rate)
-  add EfficiencyFormula {0} {0.001}
-  # efficiency formula for tau-jets
-  add EfficiencyFormula {15} {0.4}
 }
 
 #####################################################
@@ -537,6 +608,8 @@
   add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
 
   add Branch GenJetFinder/jets GenJet Jet
+  add Branch GenMissingET/momentum GenMissingET MissingET
+
   add Branch UniqueObjectFinder/jets Jet Jet
   add Branch UniqueObjectFinder/electrons Electron Electron
   add Branch UniqueObjectFinder/photons Photon Photon

=== modified file 'Template/Common/Cards/delphes_card_CMS.dat'
--- Template/Common/Cards/delphes_card_CMS.dat	2014-05-13 14:23:18 +0000
+++ Template/Common/Cards/delphes_card_CMS.dat	2016-11-23 08:16:42 +0000
@@ -10,7 +10,7 @@
   MuonTrackingEfficiency
 
   ChargedHadronMomentumSmearing
-  ElectronEnergySmearing
+  ElectronMomentumSmearing
   MuonMomentumSmearing
 
   TrackMerger
@@ -20,6 +20,7 @@
   PhotonEfficiency
   PhotonIsolation
 
+  ElectronFilter
   ElectronEfficiency
   ElectronIsolation
 
@@ -28,11 +29,16 @@
 
   MissingET
 
+  NeutrinoFilter
   GenJetFinder
+  GenMissingET
+  
   FastJetFinder
 
   JetEnergyScale
 
+  JetFlavorAssociation
+
   BTagging
   TauTagging
 
@@ -75,11 +81,11 @@
   # add EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # tracking efficiency formula for charged hadrons
-  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.60) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.85) + \
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.60) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.85) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
 
@@ -94,13 +100,13 @@
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # tracking efficiency formula for electrons
-  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.73) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e2) * (0.95) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0e2)                * (0.99) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.50) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e2) * (0.83) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e2)                * (0.90) + \
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.73) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e2) * (0.95) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0e2)                * (0.99) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.50) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e2) * (0.83) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e2)                * (0.90) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
 
@@ -115,11 +121,11 @@
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # tracking efficiency formula for muons
-  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) + \
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
 
@@ -134,32 +140,27 @@
   # set ResolutionFormula {resolution formula as a function of eta and pt}
 
   # resolution formula for charged hadrons
-  set ResolutionFormula {                  (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.02) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e1) * (0.01) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0e1 && pt <= 2.0e2) * (0.03) + \
-                                           (abs(eta) <= 1.5) * (pt > 2.0e2)                * (0.05) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.03) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e1) * (0.02) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e1 && pt <= 2.0e2) * (0.04) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 2.0e2)                * (0.05)}
+  # based on arXiv:1405.6569
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
 }
 
-#################################
-# Energy resolution for electrons
-#################################
+###################################
+# Momentum resolution for electrons
+###################################
 
-module EnergySmearing ElectronEnergySmearing {
+module MomentumSmearing ElectronMomentumSmearing {
   set InputArray ElectronTrackingEfficiency/electrons
   set OutputArray electrons
 
   # set ResolutionFormula {resolution formula as a function of eta and energy}
 
   # resolution formula for electrons
-  set ResolutionFormula {                  (abs(eta) <= 2.5) * (energy > 0.1   && energy <= 2.0e1) * (energy*0.0225) + \
-                                           (abs(eta) <= 2.5) * (energy > 2.0e1)                    * sqrt(energy^2*0.007^2 + energy*0.07^2 + 0.35^2) + \
-                         (abs(eta) > 2.5 && abs(eta) <= 3.0)                                       * sqrt(energy^2*0.007^2 + energy*0.07^2 + 0.35^2) + \
-                         (abs(eta) > 3.0 && abs(eta) <= 5.0)                                       * sqrt(energy^2*0.107^2 + energy*2.08^2)}
-
+  # based on arXiv:1405.6569
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
 }
 
 ###############################
@@ -173,18 +174,9 @@
   # set ResolutionFormula {resolution formula as a function of eta and pt}
 
   # resolution formula for muons
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1   && pt <= 5.0)   * (0.02) + \
-                                           (abs(eta) <= 0.5) * (pt > 5.0   && pt <= 1.0e2) * (0.015) + \
-                                           (abs(eta) <= 0.5) * (pt > 1.0e2 && pt <= 2.0e2) * (0.03) + \
-                                           (abs(eta) <= 0.5) * (pt > 2.0e2)                * (0.05 + pt*1.e-4) + \
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1   && pt <= 5.0)   * (0.03) + \
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 5.0   && pt <= 1.0e2) * (0.02) + \
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 1.0e2 && pt <= 2.0e2) * (0.04) + \
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 2.0e2)                * (0.05 + pt*1.e-4) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 5.0)   * (0.04) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 5.0   && pt <= 1.0e2) * (0.035) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e2 && pt <= 2.0e2) * (0.05) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 2.0e2)                * (0.05 + pt*1.e-4)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
 }
 
 ##############
@@ -194,7 +186,7 @@
 module Merger TrackMerger {
 # add InputArray InputArray
   add InputArray ChargedHadronMomentumSmearing/chargedHadrons
-  add InputArray ElectronEnergySmearing/electrons
+  add InputArray ElectronMomentumSmearing/electrons
   add InputArray MuonMomentumSmearing/muons
   set OutputArray tracks
 }
@@ -214,6 +206,14 @@
   set EFlowPhotonOutputArray eflowPhotons
   set EFlowNeutralHadronOutputArray eflowNeutralHadrons
 
+  set ECalEnergyMin 0.5
+  set HCalEnergyMin 1.0
+
+  set ECalEnergySignificanceMin 1.0
+  set HCalEnergySignificanceMin 1.0
+
+  set SmearTowerCenter true
+
   set pi [expr {acos(-1)}]
 
   # lists of the edges of each tower in eta and phi
@@ -268,11 +268,13 @@
   add EnergyFraction {3122} {0.3 0.7}
 
   # set ECalResolutionFormula {resolution formula as a function of eta and energy}
-  set ECalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.007^2 + energy*0.07^2 + 0.35^2)  + \
-                             (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
+  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
+  set ECalResolutionFormula {                  (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
 
   # set HCalResolutionFormula {resolution formula as a function of eta and energy}
-  set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) + \
+  set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
                              (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.130^2 + energy*2.70^2)}
 }
 
@@ -299,9 +301,9 @@
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # efficiency formula for photons
-  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) + \
+  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) +
                          (abs(eta) > 2.5)                                   * (0.00)}
 }
 
@@ -319,7 +321,19 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
+}
+
+#################
+# Electron filter
+#################
+
+module PdgCodeFilter ElectronFilter {
+  set InputArray Calorimeter/eflowTracks
+  set OutputArray electrons
+  set Invert true
+  add PdgCode {11}
+  add PdgCode {-11}
 }
 
 #####################
@@ -327,15 +341,15 @@
 #####################
 
 module Efficiency ElectronEfficiency {
-  set InputArray ElectronEnergySmearing/electrons
+  set InputArray ElectronFilter/electrons
   set OutputArray electrons
 
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # efficiency formula for electrons
-  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) + \
+  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) +
                          (abs(eta) > 2.5)                                   * (0.00)}
 }
 
@@ -353,7 +367,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
 }
 
 #################
@@ -367,11 +381,11 @@
   # set EfficiencyFormula {efficiency as a function of eta and pt}
 
   # efficiency formula for muons
-  set EfficiencyFormula {                                      (pt <= 10.0)               * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 10.0 && pt <= 1.0e3) * (0.95) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0 && pt <= 1.0e3) * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) + \
+  set EfficiencyFormula {                                      (pt <= 10.0)               * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
                          (abs(eta) > 2.4)                                                 * (0.00)}
 }
 
@@ -389,7 +403,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.25
 }
 
 ###################
@@ -415,12 +429,34 @@
   set EnergyOutputArray energy
 }
 
+
+#####################
+# Neutrino Filter
+#####################
+
+module PdgCodeFilter NeutrinoFilter {
+
+  set InputArray Delphes/stableParticles
+  set OutputArray filteredParticles
+
+  set PTMin 0.0
+
+  add PdgCode {12}
+  add PdgCode {14}
+  add PdgCode {16}
+  add PdgCode {-12}
+  add PdgCode {-14}
+  add PdgCode {-16}
+
+}
+
+
 #####################
 # MC truth jet finder
 #####################
 
 module FastJetFinder GenJetFinder {
-  set InputArray Delphes/stableParticles
+  set InputArray NeutrinoFilter/filteredParticles
 
   set OutputArray jets
 
@@ -431,6 +467,18 @@
   set JetPTMin 20.0
 }
 
+#########################
+# Gen Missing ET merger
+########################
+
+module Merger GenMissingET {
+# add InputArray InputArray
+  add InputArray NeutrinoFilter/filteredParticles
+  set MomentumOutputArray momentum
+}
+
+
+
 ############
 # Jet finder
 ############
@@ -456,8 +504,25 @@
   set InputArray FastJetFinder/jets
   set OutputArray jets
 
- # scale formula for jets
-  set ScaleFormula {1.00}
+  # scale formula for jets
+  set ScaleFormula {sqrt( (2.5 - 0.15*(abs(eta)))^2 / pt + 1.0 )}
+}
+
+########################
+# Jet Flavor Association
+########################
+
+module JetFlavorAssociation JetFlavorAssociation {
+
+  set PartonInputArray Delphes/partons
+  set ParticleInputArray Delphes/allParticles
+  set ParticleLHEFInputArray Delphes/allParticlesLHEF
+  set JetInputArray JetEnergyScale/jets
+
+  set DeltaR 0.5
+  set PartonPTMin 1.0
+  set PartonEtaMax 2.5
+
 }
 
 ###########
@@ -465,38 +530,30 @@
 ###########
 
 module BTagging BTagging {
-  set PartonInputArray Delphes/partons
   set JetInputArray JetEnergyScale/jets
 
   set BitNumber 0
 
-  set DeltaR 0.5
-
-  set PartonPTMin 1.0
-
-  set PartonEtaMax 2.5
-
   # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt}
   # PDG code = the highest PDG code of a quark or gluon inside DeltaR cone around jet axis
   # gluon's PDG code has the lowest priority
 
-  # https://twiki.cern.ch/twiki/bin/view/CMSPublic/PhysicsResultsBTV
+  # based on arXiv:1211.4462
+  
   # default efficiency formula (misidentification rate)
-  add EfficiencyFormula {0} {0.001}
+  add EfficiencyFormula {0} {0.01+0.000038*pt}
 
   # efficiency formula for c-jets (misidentification rate)
-  add EfficiencyFormula {4} {                                      (pt <= 15.0) * (0.000) + \
-                                                (abs(eta) <= 1.2) * (pt > 15.0) * (0.2*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 1.2 && abs(eta) <= 2.5) * (pt > 15.0) * (0.1*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 2.5)                                  * (0.000)}
+  add EfficiencyFormula {4} {0.25*tanh(0.018*pt)*(1/(1+ 0.0013*pt))}
 
   # efficiency formula for b-jets
-  add EfficiencyFormula {5} {                                      (pt <= 15.0) * (0.000) + \
-                                                (abs(eta) <= 1.2) * (pt > 15.0) * (0.5*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 1.2 && abs(eta) <= 2.5) * (pt > 15.0) * (0.4*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 2.5)                                  * (0.000)}
+  add EfficiencyFormula {5} {0.85*tanh(0.0025*pt)*(25.0/(1+0.063*pt))}
 }
 
+#############
+# tau-tagging
+#############
+
 module TauTagging TauTagging {
   set ParticleInputArray Delphes/allParticles
   set PartonInputArray Delphes/partons
@@ -511,9 +568,9 @@
   # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt}
 
   # default efficiency formula (misidentification rate)
-  add EfficiencyFormula {0} {0.001}
+  add EfficiencyFormula {0} {0.01}
   # efficiency formula for tau-jets
-  add EfficiencyFormula {15} {0.4}
+  add EfficiencyFormula {15} {0.6}
 }
 
 #####################################################
@@ -549,6 +606,8 @@
   add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
 
   add Branch GenJetFinder/jets GenJet Jet
+  add Branch GenMissingET/momentum GenMissingET MissingET
+
   add Branch UniqueObjectFinder/jets Jet Jet
   add Branch UniqueObjectFinder/electrons Electron Electron
   add Branch UniqueObjectFinder/photons Photon Photon
@@ -556,4 +615,3 @@
   add Branch MissingET/momentum MissingET MissingET
   add Branch ScalarHT/energy ScalarHT ScalarHT
 }
-

=== modified file 'Template/Common/Cards/delphes_card_default.dat'
--- Template/Common/Cards/delphes_card_default.dat	2015-06-17 13:28:26 +0000
+++ Template/Common/Cards/delphes_card_default.dat	2016-11-23 08:16:42 +0000
@@ -10,7 +10,7 @@
   MuonTrackingEfficiency
 
   ChargedHadronMomentumSmearing
-  ElectronEnergySmearing
+  ElectronMomentumSmearing
   MuonMomentumSmearing
 
   TrackMerger
@@ -20,6 +20,7 @@
   PhotonEfficiency
   PhotonIsolation
 
+  ElectronFilter
   ElectronEfficiency
   ElectronIsolation
 
@@ -28,11 +29,16 @@
 
   MissingET
 
+  NeutrinoFilter
   GenJetFinder
+  GenMissingET
+  
   FastJetFinder
 
   JetEnergyScale
 
+  JetFlavorAssociation
+
   BTagging
   TauTagging
 
@@ -75,11 +81,11 @@
   # add EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # tracking efficiency formula for charged hadrons
-  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.60) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.85) + \
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.60) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.85) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
 
@@ -94,13 +100,13 @@
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # tracking efficiency formula for electrons
-  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.73) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e2) * (0.95) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0e2)                * (0.99) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.50) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e2) * (0.83) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e2)                * (0.90) + \
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.73) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e2) * (0.95) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0e2)                * (0.99) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.50) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e2) * (0.83) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e2)                * (0.90) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
 
@@ -115,11 +121,11 @@
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # tracking efficiency formula for muons
-  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) + \
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
 
@@ -134,32 +140,27 @@
   # set ResolutionFormula {resolution formula as a function of eta and pt}
 
   # resolution formula for charged hadrons
-  set ResolutionFormula {                  (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.02) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e1) * (0.01) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0e1 && pt <= 2.0e2) * (0.03) + \
-                                           (abs(eta) <= 1.5) * (pt > 2.0e2)                * (0.05) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.03) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e1) * (0.02) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e1 && pt <= 2.0e2) * (0.04) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 2.0e2)                * (0.05)}
+  # based on arXiv:1405.6569
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
 }
 
-#################################
-# Energy resolution for electrons
-#################################
+###################################
+# Momentum resolution for electrons
+###################################
 
-module EnergySmearing ElectronEnergySmearing {
+module MomentumSmearing ElectronMomentumSmearing {
   set InputArray ElectronTrackingEfficiency/electrons
   set OutputArray electrons
 
   # set ResolutionFormula {resolution formula as a function of eta and energy}
 
   # resolution formula for electrons
-  set ResolutionFormula {                  (abs(eta) <= 2.5) * (energy > 0.1   && energy <= 2.0e1) * (energy*0.0225) + \
-                                           (abs(eta) <= 2.5) * (energy > 2.0e1)                    * sqrt(energy^2*0.007^2 + energy*0.07^2 + 0.35^2) + \
-                         (abs(eta) > 2.5 && abs(eta) <= 3.0)                                       * sqrt(energy^2*0.007^2 + energy*0.07^2 + 0.35^2) + \
-                         (abs(eta) > 3.0 && abs(eta) <= 5.0)                                       * sqrt(energy^2*0.107^2 + energy*2.08^2)}
-
+  # based on arXiv:1405.6569
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
 }
 
 ###############################
@@ -173,18 +174,9 @@
   # set ResolutionFormula {resolution formula as a function of eta and pt}
 
   # resolution formula for muons
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1   && pt <= 5.0)   * (0.02) + \
-                                           (abs(eta) <= 0.5) * (pt > 5.0   && pt <= 1.0e2) * (0.015) + \
-                                           (abs(eta) <= 0.5) * (pt > 1.0e2 && pt <= 2.0e2) * (0.03) + \
-                                           (abs(eta) <= 0.5) * (pt > 2.0e2)                * (0.05 + pt*1.e-4) + \
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1   && pt <= 5.0)   * (0.03) + \
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 5.0   && pt <= 1.0e2) * (0.02) + \
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 1.0e2 && pt <= 2.0e2) * (0.04) + \
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 2.0e2)                * (0.05 + pt*1.e-4) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 5.0)   * (0.04) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 5.0   && pt <= 1.0e2) * (0.035) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e2 && pt <= 2.0e2) * (0.05) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 2.0e2)                * (0.05 + pt*1.e-4)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
 }
 
 ##############
@@ -194,7 +186,7 @@
 module Merger TrackMerger {
 # add InputArray InputArray
   add InputArray ChargedHadronMomentumSmearing/chargedHadrons
-  add InputArray ElectronEnergySmearing/electrons
+  add InputArray ElectronMomentumSmearing/electrons
   add InputArray MuonMomentumSmearing/muons
   set OutputArray tracks
 }
@@ -214,6 +206,14 @@
   set EFlowPhotonOutputArray eflowPhotons
   set EFlowNeutralHadronOutputArray eflowNeutralHadrons
 
+  set ECalEnergyMin 0.5
+  set HCalEnergyMin 1.0
+
+  set ECalEnergySignificanceMin 1.0
+  set HCalEnergySignificanceMin 1.0
+
+  set SmearTowerCenter true
+
   set pi [expr {acos(-1)}]
 
   # lists of the edges of each tower in eta and phi
@@ -268,11 +268,13 @@
   add EnergyFraction {3122} {0.3 0.7}
 
   # set ECalResolutionFormula {resolution formula as a function of eta and energy}
-  set ECalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.007^2 + energy*0.07^2 + 0.35^2)  + \
-                             (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
+  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
+  set ECalResolutionFormula {                  (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
 
   # set HCalResolutionFormula {resolution formula as a function of eta and energy}
-  set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) + \
+  set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
                              (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.130^2 + energy*2.70^2)}
 }
 
@@ -299,9 +301,9 @@
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # efficiency formula for photons
-  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) + \
+  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) +
                          (abs(eta) > 2.5)                                   * (0.00)}
 }
 
@@ -319,7 +321,19 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
+}
+
+#################
+# Electron filter
+#################
+
+module PdgCodeFilter ElectronFilter {
+  set InputArray Calorimeter/eflowTracks
+  set OutputArray electrons
+  set Invert true
+  add PdgCode {11}
+  add PdgCode {-11}
 }
 
 #####################
@@ -327,15 +341,15 @@
 #####################
 
 module Efficiency ElectronEfficiency {
-  set InputArray ElectronEnergySmearing/electrons
+  set InputArray ElectronFilter/electrons
   set OutputArray electrons
 
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
 
   # efficiency formula for electrons
-  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) + \
+  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) +
                          (abs(eta) > 2.5)                                   * (0.00)}
 }
 
@@ -353,7 +367,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.12
 }
 
 #################
@@ -367,11 +381,11 @@
   # set EfficiencyFormula {efficiency as a function of eta and pt}
 
   # efficiency formula for muons
-  set EfficiencyFormula {                                      (pt <= 10.0)               * (0.00) + \
-                                           (abs(eta) <= 1.5) * (pt > 10.0 && pt <= 1.0e3) * (0.95) + \
-                                           (abs(eta) <= 1.5) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0 && pt <= 1.0e3) * (0.95) + \
-                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) + \
+  set EfficiencyFormula {                                      (pt <= 10.0)               * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
                          (abs(eta) > 2.4)                                                 * (0.00)}
 }
 
@@ -389,7 +403,7 @@
 
   set PTMin 0.5
 
-  set PTRatioMax 0.1
+  set PTRatioMax 0.25
 }
 
 ###################
@@ -415,12 +429,34 @@
   set EnergyOutputArray energy
 }
 
+
+#####################
+# Neutrino Filter
+#####################
+
+module PdgCodeFilter NeutrinoFilter {
+
+  set InputArray Delphes/stableParticles
+  set OutputArray filteredParticles
+
+  set PTMin 0.0
+
+  add PdgCode {12}
+  add PdgCode {14}
+  add PdgCode {16}
+  add PdgCode {-12}
+  add PdgCode {-14}
+  add PdgCode {-16}
+
+}
+
+
 #####################
 # MC truth jet finder
 #####################
 
 module FastJetFinder GenJetFinder {
-  set InputArray Delphes/stableParticles
+  set InputArray NeutrinoFilter/filteredParticles
 
   set OutputArray jets
 
@@ -431,6 +467,18 @@
   set JetPTMin 20.0
 }
 
+#########################
+# Gen Missing ET merger
+########################
+
+module Merger GenMissingET {
+# add InputArray InputArray
+  add InputArray NeutrinoFilter/filteredParticles
+  set MomentumOutputArray momentum
+}
+
+
+
 ############
 # Jet finder
 ############
@@ -456,8 +504,25 @@
   set InputArray FastJetFinder/jets
   set OutputArray jets
 
- # scale formula for jets
-  set ScaleFormula {1.00}
+  # scale formula for jets
+  set ScaleFormula {sqrt( (2.5 - 0.15*(abs(eta)))^2 / pt + 1.0 )}
+}
+
+########################
+# Jet Flavor Association
+########################
+
+module JetFlavorAssociation JetFlavorAssociation {
+
+  set PartonInputArray Delphes/partons
+  set ParticleInputArray Delphes/allParticles
+  set ParticleLHEFInputArray Delphes/allParticlesLHEF
+  set JetInputArray JetEnergyScale/jets
+
+  set DeltaR 0.5
+  set PartonPTMin 1.0
+  set PartonEtaMax 2.5
+
 }
 
 ###########
@@ -465,38 +530,30 @@
 ###########
 
 module BTagging BTagging {
-  set PartonInputArray Delphes/partons
   set JetInputArray JetEnergyScale/jets
 
   set BitNumber 0
 
-  set DeltaR 0.5
-
-  set PartonPTMin 1.0
-
-  set PartonEtaMax 2.5
-
   # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt}
   # PDG code = the highest PDG code of a quark or gluon inside DeltaR cone around jet axis
   # gluon's PDG code has the lowest priority
 
-  # https://twiki.cern.ch/twiki/bin/view/CMSPublic/PhysicsResultsBTV
+  # based on arXiv:1211.4462
+  
   # default efficiency formula (misidentification rate)
-  add EfficiencyFormula {0} {0.001}
+  add EfficiencyFormula {0} {0.01+0.000038*pt}
 
   # efficiency formula for c-jets (misidentification rate)
-  add EfficiencyFormula {4} {                                      (pt <= 15.0) * (0.000) + \
-                                                (abs(eta) <= 1.2) * (pt > 15.0) * (0.2*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 1.2 && abs(eta) <= 2.5) * (pt > 15.0) * (0.1*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 2.5)                                  * (0.000)}
+  add EfficiencyFormula {4} {0.25*tanh(0.018*pt)*(1/(1+ 0.0013*pt))}
 
   # efficiency formula for b-jets
-  add EfficiencyFormula {5} {                                      (pt <= 15.0) * (0.000) + \
-                                                (abs(eta) <= 1.2) * (pt > 15.0) * (0.5*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 1.2 && abs(eta) <= 2.5) * (pt > 15.0) * (0.4*tanh(pt*0.03 - 0.4)) + \
-                              (abs(eta) > 2.5)                                  * (0.000)}
+  add EfficiencyFormula {5} {0.85*tanh(0.0025*pt)*(25.0/(1+0.063*pt))}
 }
 
+#############
+# tau-tagging
+#############
+
 module TauTagging TauTagging {
   set ParticleInputArray Delphes/allParticles
   set PartonInputArray Delphes/partons
@@ -511,9 +568,9 @@
   # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt}
 
   # default efficiency formula (misidentification rate)
-  add EfficiencyFormula {0} {0.001}
+  add EfficiencyFormula {0} {0.01}
   # efficiency formula for tau-jets
-  add EfficiencyFormula {15} {0.4}
+  add EfficiencyFormula {15} {0.6}
 }
 
 #####################################################
@@ -549,6 +606,8 @@
   add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
 
   add Branch GenJetFinder/jets GenJet Jet
+  add Branch GenMissingET/momentum GenMissingET MissingET
+
   add Branch UniqueObjectFinder/jets Jet Jet
   add Branch UniqueObjectFinder/electrons Electron Electron
   add Branch UniqueObjectFinder/photons Photon Photon
@@ -556,4 +615,3 @@
   add Branch MissingET/momentum MissingET MissingET
   add Branch ScalarHT/energy ScalarHT ScalarHT
 }
-

=== renamed file 'Template/LO/Cards/reweight_card_default.dat' => 'Template/Common/Cards/reweight_card_default.dat'
--- Template/LO/Cards/reweight_card_default.dat	2013-11-15 12:29:09 +0000
+++ Template/Common/Cards/reweight_card_default.dat	2016-11-23 08:16:42 +0000
@@ -1,25 +1,69 @@
-#******************************************************************
-#                       Reweight Module                           *
-#******************************************************************
+#*************************************************************************
+#                          Reweight Module                               *
+#                Matrix-Element reweighting at LO/NLO                    *
+#  Mattelaer Olivier                                    arxiv:xxxx.xxxx  *  
+#*************************************************************************
+#          
+# Note: 
+#   1) the value of alpha_s will be used from the event so the value in
+#      the param_card is not taken into account.
+#   2) It is (in general) dangerous/wrong to change parameters by a large
+#      amount, if this changes the shape of the matrix elements a lot.
+#      (For example, changing a particle's mass by much more than its
+#      width leads to very inaccurate result). In such a case, separate 
+#      event generation runs are needed.
+#
+#************************************************************************
+# ENTER YOUR COMMANDS BELOW. 
+#************************************************************************
+
+change mode NLO    # Define type of Reweighting. For LO sample this command
+                   # has no effect since only "LO" mode is allowed.
+                   
 launch
-#* Use the set command to specify the new set of parameter 
-#* Or specify a path to a valid param_card/banner
-#* Example of valid command:
-#*     set aewm1 137
-#*     ~/param_card.dat
-#*
-#* Note: 
-#*   1) the value of alphas will be used from the event
-#*      so the value of the param_card is not taken into account.
-#*   2) It is dangerous to change a mass of any particle. 
-
-
-#* If you want to compute the weight for more than one hyppothesis
-#* you need first to uncomment the following line:
-# launch
-# and then use the set command to specify your parameter.
-# All modification will start from the ORIGINAL card not from the 
-# last define one.
-#* You can have as many weight as you want.
-
+#  SPECIFY A PATH OR USE THE SET COMMAND LIKE THIS:
+# set sminputs 1 130 # modify 1/alpha_EW
+
+
+
+#************************************************************************
+# Manual:      https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Reweight
+#************************************************************************
+#
+# Example of (standard) code for the computation of two weights:
+#
+# launch                    ! tag to start the computation of the first weight
+#  set BLOCKNAME ID VALUE   ! rule to modify the current param_card
+#  set BLOCKNAME ID VALUE   ! rule to modify the current param_card
+# launch                    ! start to compute a second weight
+#  /home/Cards/param_card_2.dat ! you can also enter a path to a valid card
+#
+# Note: The command to specify the parameter are AFTER the associated "launch" 
+#
+# Possible options: 
+#   You can enter one of the following lines to customize the reweighting 
+#   procedure. These need to be given before the 'launch' command.
+#
+# change model NAME : use another model for the matrix-elements to reweight 
+#        with. In this case you need to provide the path to a correct 
+#        param_card for the new model; you cannot modify the original one 
+#        with the 'set' command.
+# change process DEF [--add]: change the process by which you reweight. 
+#        The initial and final state particles of the new process should 
+#        be exactly identical to the ones in the original process.
+# change helicity False: perform the reweighting by helicity summed 
+#        matrix-elements even if the events have been written with a 
+#        single helicity state.       
+# change mode XXX: change the type of reweighting performed.
+#        allowed values: LO, NLO, LO+NLO
+#        - This command has no effect for reweighting an .lhe event file with LO accuracy. 
+#          In that case LO mode is always used (whatever entry is set).
+#        - When the .lhe file reweighted is at NLO accuracy, then all modes are allowed.
+#          * "LO" is an approximate leading order method
+#          * "NLO" is the NLO accurate method
+#          * "LO+NLO" runs both
+#       - "NLO" and "LO+NLO" modes requires 'store_rwgt_info' equals True (run_card.dat)
+#          If the reweighting is done at generation level this parameter will
+#          automatically be set on True.
+#************************************************************************
 

=== removed file 'Template/Common/bin/internal/run_pythia'
--- Template/Common/bin/internal/run_pythia	2013-11-05 10:16:20 +0000
+++ Template/Common/bin/internal/run_pythia	1970-01-01 00:00:00 +0000
@@ -1,23 +0,0 @@
-#!/bin/bash
-#
-#  This runs pythia on the unweighted_events.dat
-#
-#  Usage: run_pythia [pydir]
-#  where pydir is the path to the pythia executable
-
-pydir=$1
-main=`pwd`
-
-if [ ! -e ../Cards/pythia_card.dat ]; then
-  echo "No pythia_card.dat found. Quitting..."
-  exit
-fi
-
-echo $$ >> ../myprocid
-# shower and hadronize event through Pythia
-echo " " >> ../Cards/pythia_card.dat
-echo "      LHAPATH=$pydir/PDFsets" >> ../Cards/pythia_card.dat
-export PDG_MASS_TBL=$pydir/mass_width_2004.mc
-rm -rf Events/hepmv_event.hepmc &> /dev/null
-$pydir/pythia && touch pythia.done
-

=== added file 'Template/LO/Cards/madanalysis5_hadron_card_default.dat'
--- Template/LO/Cards/madanalysis5_hadron_card_default.dat	1970-01-01 00:00:00 +0000
+++ Template/LO/Cards/madanalysis5_hadron_card_default.dat	2016-11-23 08:16:42 +0000
@@ -0,0 +1,3 @@
+# This card is used only if MA5 failed to create a default for this run
+# We therefore use as default: do nothing
+@MG5aMC skip_analysis
\ No newline at end of file

=== added file 'Template/LO/Cards/madanalysis5_parton_card_default.dat'
--- Template/LO/Cards/madanalysis5_parton_card_default.dat	1970-01-01 00:00:00 +0000
+++ Template/LO/Cards/madanalysis5_parton_card_default.dat	2016-11-23 08:16:42 +0000
@@ -0,0 +1,3 @@
+# This card is used only if MA5 failed to create a default for this run
+# We therefore use as default: do nothing
+@MG5aMC skip_analysis
\ No newline at end of file

=== added file 'Template/LO/Cards/pythia8_card_default.dat'
--- Template/LO/Cards/pythia8_card_default.dat	1970-01-01 00:00:00 +0000
+++ Template/LO/Cards/pythia8_card_default.dat	2016-11-23 08:16:42 +0000
@@ -0,0 +1,77 @@
+!
+! Pythia8 cmd card automatically generated by MadGraph5_aMC@NLO
+! For more information on the use of the MG5aMC / Pythia8 interface, visit
+!    https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/LOPY8Merging
+!
+! ==================
+! General parameters 
+! ==================
+!
+Main:numberOfEvents      = -1
+!
+! -------------------------------------------------------------------
+! Specify the HEPMC output of the Pythia8 shower. You can set it to:
+!   auto      : MG5aMC will automatically place it the run_<i> directory
+!   /dev/null : to turn off the HEPMC output.
+!   <path>    : to select where the HEPMC file must written. It will 
+!               therefore not be placed in the run_<i> directory. The
+!               specified path, if not absolute, will be relative to 
+!               the Event/run_<i> directory of the process output.
+!   fifo      : to have MG5aMC setup the piping of the PY8 output to 
+!               analysis tools such as MadAnalysis5.  
+!   fifo@<fifo_path> :
+!               Same as 'fifo', but selecting a custom path to create the
+!               fifo pipe. (useful to select a mounted drive that supports 
+!               fifo). Note that the fifo file extension *must* be '.hepmc.fifo'.
+! -------------------------------------------------------------------
+!
+HEPMCoutput:file         = auto
+!
+! --------------------------------------------------------------------
+! Parameters relevant only when performing MLM merging, which can be
+! turned on by setting ickkw to '1' in the run_card and chosing a 
+! positive value for the parameter xqcut.
+! For details, see section 'Jet Matching' on the left-hand menu of 
+!    http://home.thep.lu.se/~torbjorn/pythia81html/Welcome.html
+! --------------------------------------------------------------------
+! If equal to -1.0, MadGraph5_aMC@NLO will set it automatically based 
+! on the parameter 'xqcut' of the run_card.dat 
+JetMatching:qCut         = -1.0
+! Use default kt-MLM to match parton level jets to those produced by the
+! shower. But the other Shower-kt scheme is available too with this option.
+JetMatching:doShowerKt   = off
+! A value of -1 means that it is automatically guessed by MadGraph.
+! It is however always safer to explicitly set it.
+JetMatching:nJetMax      = -1
+!
+! --------------------------------------------------------------------
+! Parameters relevant only when performing CKKW-L merging, which can
+! be turned on by setting the parameter 'ptlund' *or* 'ktdurham' to
+! a positive value. 
+! For details, see section 'CKKW-L Merging' on the left-hand menu of 
+!    http://home.thep.lu.se/~torbjorn/pythia81html/Welcome.html
+! --------------------------------------------------------------------
+! Central merging scale values you want to be used.
+! If equal to -1.0, then MadGraph5_aMC@NLO will set this automatically
+! based on the parameter 'ktdurham' of the run_card.dat
+Merging:TMS              = -1.0
+! This must be set manually, according to Pythia8 directives.
+! An example of possible value is 'pp>LEPTONS,NEUTRINOS'
+Merging:Process          = <set_by_user>
+! A value of -1 means that it is automatically guessed by MadGraph.
+! It is however always safer to explicitly set it.
+Merging:nJetMax  		 = -1
+!
+! For all merging schemes, decide wehter you want the merging scale
+! variation computed for only the central weights or all other 
+! PDF and scale variation weights as well
+SysCalc:fullCutVariation = off
+!
+! ==========================
+! User customized parameters 
+! ==========================
+!
+! By default, Pythia8 generates multi-parton interaction events. This is
+! often irrelevant for phenomenology and very slow. You can turn this 
+! feature off by uncommenting the line below if so desired.
+!partonlevel:mpi = off

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2015-07-30 15:57:44 +0000
+++ Template/LO/Cards/run_card.dat	2016-11-23 08:16:42 +0000
@@ -21,10 +21,6 @@
 #*********************************************************************
   %(run_tag)s     = run_tag ! name of the run 
 #*********************************************************************
-# Run to generate the grid pack                                      *
-#*********************************************************************
-  %(gridpack)s     = gridpack  !True = setting up the grid pack
-#*********************************************************************
 # Number of events and rnd seed                                      *
 # Warning: Do not generate more than 1M events in a single run       *
 # If you want to run Pythia, avoid more than 50k events in a run.    *
@@ -61,51 +57,55 @@
  %(dynamical_scale_choice)s = dynamical_scale_choice ! Choose one of the preselected dynamical choices
  %(scalefact)s  = scalefact        ! scale factor for event-by-event scales
 #*********************************************************************
-# Time of flight information. (-1 means not run)
-#*********************************************************************
- %(time_of_flight)s = time_of_flight ! threshold below which info is not written
-#*********************************************************************
-# Matching - Warning! ickkw > 1 is still beta
-#*********************************************************************
- %(ickkw)s = ickkw            ! 0 no matching, 1 MLM, 2 CKKW matching
- %(highestmult)s = highestmult      ! for ickkw=2, highest mult group
- %(ktscheme)s = ktscheme         ! for ickkw=1, 1 Durham kT, 2 Pythia pTE
+# Type and output format
+#*********************************************************************
+  %(gridpack)s     = gridpack  !True = setting up the grid pack
+  %(time_of_flight)s = time_of_flight ! threshold (in mm) below which the invariant livetime is not written (-1 means not written)
+  %(lhe_version)s = lhe_version       ! Change the way clustering information pass to shower.        
+  %(clusinfo)s = clusinfo         ! include clustering tag in output
+  %(event_norm)s =  event_norm       ! average/sum. Normalization of the weight in the LHEF
+
+#*********************************************************************
+# Matching parameter (MLM only)
+#*********************************************************************
+ %(ickkw)s = ickkw            ! 0 no matching, 1 MLM
  %(alpsfact)s = alpsfact         ! scale factor for QCD emission vx
  %(chcluster)s = chcluster        ! cluster only according to channel diag
- %(pdfwgt)s = pdfwgt           ! for ickkw=1, perform pdf reweighting
  %(asrwgtflavor)s = asrwgtflavor     ! highest quark flavor for a_s reweight
- %(clusinfo)s = clusinfo         ! include clustering tag in output
- %(lhe_version)s = lhe_version       ! Change the way clustering information pass to shower.        
-#*********************************************************************
-#**********************************************************
-#
-#**********************************************************
-# Automatic ptj and mjj cuts if xqcut > 0
-# (turn off for VBF and single top processes)
-#**********************************************************
-   %(auto_ptj_mjj)s  = auto_ptj_mjj  ! Automatic setting of ptj and mjj
-#**********************************************************
+ %(auto_ptj_mjj)s  = auto_ptj_mjj  ! Automatic setting of ptj and mjj if xqcut >0
+                                   ! (turn off for VBF and single top processes) 
+ %(xqcut)s   = xqcut   ! minimum kt jet measure between partons
+#*********************************************************************
+#
+#*********************************************************************
+# handling of the helicities:
+#  0: sum over all helicities
+#  1: importance sampling over helicities
+#*********************************************************************
+   %(nhel)s  = nhel          ! using helicities importance sampling or not.
+#*********************************************************************
+# Generation bias, check the wiki page below for more information:   *
+#  'cp3.irmp.ucl.ac.be/projects/madgraph/wiki/LOEventGenerationBias' *
+#*********************************************************************
+ %(bias_module)s = bias_module  ! Bias type of bias, [None, ptj_bias, -custom_folder-]
+ %(bias_parameters)s = bias_parameters ! Specifies the parameters of the module.
+#
+#*******************************                                                 
+# Parton level cuts definition *
+#*******************************                                     
 #                                                                    
-#**********************************
-# BW cutoff (M+/-bwcutoff*Gamma)
-#**********************************
+#
+#*********************************************************************
+# BW cutoff (M+/-bwcutoff*Gamma) ! Define on/off-shell for "$" and decay  
+#*********************************************************************
   %(bwcutoff)s  = bwcutoff      ! (M+/-bwcutoff*Gamma)
-#**********************************************************
-# Apply pt/E/eta/dr/mij cuts on decay products or not
+#*********************************************************************
+# Apply pt/E/eta/dr/mij/kt_durham cuts on decay products or not
 # (note that etmiss/ptll/ptheavy/ht/sorted cuts always apply)
-#**********************************************************
+#*********************************************************************
    %(cut_decays)s  = cut_decays    ! Cut decay products 
-#*************************************************************
-# Number of helicities to sum per event (0 = all helicities)
-# 0 gives more stable result, but longer run time (needed for
-# long decay chains e.g.).
-# Use >=2 if most helicities contribute, e.g. pure QCD.
-#*************************************************************
-   %(nhel)s  = nhel          ! Number of helicities used per event
-#*******************                                                 
-# Standard Cuts
-#*******************                                                 
-#                                                                    
+#*********************************************************************
+# Standard Cuts                                                      *
 #*********************************************************************
 # Minimum and maximum pt's (for max, -1 means no cut)                *
 #*********************************************************************
@@ -246,27 +246,24 @@
 #*********************************************************************
  %(xetamin)s   = xetamin ! minimum rapidity for two jets in the WBF case  
  %(deltaeta)s   = deltaeta ! minimum rapidity for two jets in the WBF case 
-#*********************************************************************
-# KT DURHAM CUT                                                      *
-#*********************************************************************
- %(ktdurham)s    =  ktdurham        
- %(dparameter)s  =  dparameter 
+#***********************************************************************
+# Turn on either the ktdurham or ptlund cut to activate                *
+# CKKW(L) merging with Pythia8 [arXiv:1410.3012, arXiv:1109.4829]      *
+#***********************************************************************
+ %(ktdurham)s  =  ktdurham        
+ %(dparameter)s   =  dparameter
+ %(ptlund)s  =  ptlund
+ %(pdgs_for_merging_cut)s  =  pdgs_for_merging_cut ! PDGs for two cuts above   
 #*********************************************************************
 # maximal pdg code for quark to be considered as a light jet         *
 # (otherwise b cuts are applied)                                     *
 #*********************************************************************
  %(maxjetflavor)s = maxjetflavor    ! Maximum jet pdg code
 #*********************************************************************
-# Jet measure cuts                                                   *
-#*********************************************************************
- %(xqcut)s   = xqcut   ! minimum kt jet measure between partons
-#*********************************************************************
 #
 #*********************************************************************
 # Store info for systematics studies                                 *
-# WARNING: If use_syst is T, matched Pythia output is                *
-#          meaningful ONLY if plotted taking matchscale              *
-#          reweighting into account!                                 *
+# WARNING: Do not use for interference type of computation           *
 #*********************************************************************
    %(use_syst)s  = use_syst      ! Enable systematics studies
 #
@@ -279,5 +276,6 @@
 %(sys_alpsfact)s = sys_alpsfact  # \alpha_s emission scale factors
 %(sys_matchscale)s = sys_matchscale # variation of merging scale
 # PDF sets and number of members (0 or none for all members).
-%(sys_pdf)s = sys_pdf # matching scales
+%(sys_pdf)s = sys_pdf # separate by && if more than one set.
 # MSTW2008nlo68cl.LHgrid 1  = sys_pdf
+#

=== added directory 'Template/LO/Source/BIAS'
=== added directory 'Template/LO/Source/BIAS/dummy'
=== added file 'Template/LO/Source/BIAS/dummy/dummy.f'
--- Template/LO/Source/BIAS/dummy/dummy.f	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/dummy/dummy.f	2016-11-23 08:16:42 +0000
@@ -0,0 +1,44 @@
+C ************************************************************
+C Source for the library implementing a dummt bias function 
+C always returns one
+C ************************************************************
+      
+      subroutine bias_wgt(p, original_weight, bias_weight)
+          implicit none
+C
+C Parameters
+C
+          include '../../nexternal.inc'
+C
+C Arguments
+C
+          double precision p(0:3,nexternal)
+          double precision original_weight, bias_weight
+C
+C local variables
+C
+C
+C Global variables
+C
+C Mandatory common block to be defined in bias modules
+C
+          double precision stored_bias_weight
+          data stored_bias_weight/1.0d0/          
+          logical impact_xsec, requires_full_event_info
+C         Not impacting the xsec since the bias is 1.0. Therefore
+C         bias_wgt will not be written in the lhe event file.
+C         Setting it to .True. makes sure that it will not be written.
+          data impact_xsec/.True./
+C         Of course this module does not require the full event
+C         information (color, resonances, helicities, etc..)
+          data requires_full_event_info/.False./ 
+          common/bias/stored_bias_weight,impact_xsec,
+     &                requires_full_event_info
+
+C --------------------
+C BEGIN IMPLEMENTATION
+C --------------------
+
+          bias_weight = 1.0d0
+
+      end subroutine bias_wgt

=== added file 'Template/LO/Source/BIAS/dummy/makefile'
--- Template/LO/Source/BIAS/dummy/makefile	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/dummy/makefile	2016-11-23 08:16:42 +0000
@@ -0,0 +1,21 @@
+
+include ../../make_opts
+
+all: dummy 
+
+clean:
+	$(RM) *.o $(BIASLIBDIR)$(BIASLIBRARY)
+
+#
+# Compilation of the module dummy
+#
+
+dummy: dummy.o
+	$(call CREATELIB, $(BIASLIBDIR)$(BIASLIBRARY), $^)
+
+#
+# List of the requirements for this module.
+# 'VALID' is the keyword that *must* be returned if everything is in order.
+#
+requirements:	
+	@echo "VALID"

=== added directory 'Template/LO/Source/BIAS/ptj_bias'
=== added file 'Template/LO/Source/BIAS/ptj_bias/makefile'
--- Template/LO/Source/BIAS/ptj_bias/makefile	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/ptj_bias/makefile	2016-11-23 08:16:42 +0000
@@ -0,0 +1,27 @@
+include ../../make_opts
+
+all: ptj_bias
+
+clean:
+	$(RM) *.o $(BIASLIBDIR)$(BIASLIBRARY)
+
+#
+# Compilation of the module ptj_bias
+#
+ptj_bias.o: ptj_bias.f ../bias.inc
+	$(FC) $(FFLAGS) $(LDFLAGS) -c -o ptj_bias.o ptj_bias.f
+
+ptj_bias: ptj_bias.o
+	$(call CREATELIB, $(BIASLIBDIR)$(BIASLIBRARY), $^)
+
+#
+# List of the requirements for this module.
+# 'VALID' is the keyword that *must* be returned if everything is in order.
+#
+requirements:	
+ifeq ($(shell $(call CHECK_MG5AMC_VERSION,2.4.2)),True)
+	@echo "VALID"
+else
+	@echo "Error:: MG5aMC is not recent enough (found "$(MG5AMC_VERSION)")"
+	@echo "FAIL"
+endif

=== added file 'Template/LO/Source/BIAS/ptj_bias/ptj_bias.f'
--- Template/LO/Source/BIAS/ptj_bias/ptj_bias.f	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/BIAS/ptj_bias/ptj_bias.f	2016-11-23 08:16:42 +0000
@@ -0,0 +1,101 @@
+C ************************************************************
+C Source for the library implementing a bias function that 
+C populates the large pt tale of the leading jet. 
+C
+C The two options of this subroutine, that can be set in
+C the run card are:
+C    > (double precision) ptj_bias_target_ptj : target ptj value
+C    > (double precision) ptj_bias_enhancement_power : exponent
+C
+C Schematically, the functional form of the enhancement is
+C    bias_wgt = [ptj(evt)/mean_ptj]^enhancement_power
+C ************************************************************
+C
+C The following lines are read by MG5aMC to set what are the 
+C relevant parameters for this bias module.
+C
+C  parameters = {'ptj_bias_target_ptj': 1000.0,
+C               'ptj_bias_enhancement_power': 4.0}
+C
+
+      subroutine bias_wgt(p, original_weight, bias_weight)
+          implicit none
+C
+C Parameters
+C
+          include '../../maxparticles.inc'          
+          include '../../nexternal.inc'
+
+C
+C Arguments
+C
+          double precision p(0:3,nexternal)
+          double precision original_weight, bias_weight
+C
+C local variables
+C
+          integer i
+          double precision ptj(nexternal)
+          double precision max_ptj
+c
+c local variables defined in the run_card
+c
+          double precision ptj_bias_target_ptj
+          double precision ptj_bias_enhancement_power
+C
+C Global variables
+C
+C
+C Mandatory common block to be defined in bias modules
+C
+          double precision stored_bias_weight
+          data stored_bias_weight/1.0d0/          
+          logical impact_xsec, requires_full_event_info
+C         We only want to bias distributions, but not impact the xsec. 
+          data impact_xsec/.False./
+C         Of course this module does not require the full event
+C         information (color, resonances, helicities, etc..)
+          data requires_full_event_info/.False./ 
+          common/bias/stored_bias_weight,impact_xsec,
+     &                requires_full_event_info
+C
+C Accessingt the details of the event
+C
+          logical is_a_j(nexternal),is_a_l(nexternal),
+     &            is_a_b(nexternal),is_a_a(nexternal),
+     &            is_a_onium(nexternal),is_a_nu(nexternal),
+     &            is_heavy(nexternal),do_cuts(nexternal)
+          common/to_specisa/is_a_j,is_a_a,is_a_l,is_a_b,is_a_nu,
+     &                      is_heavy,is_a_onium,do_cuts
+
+C
+C    Setup the value of the parameters from the run_card    
+C
+      include '../bias.inc'
+
+C --------------------
+C BEGIN IMPLEMENTATION
+C --------------------
+          
+          do i=1,nexternal
+            ptj(i)=-1.0d0
+            if (is_a_j(i)) then
+              ptj(i)=sqrt(p(1,i)**2+p(2,i)**2)
+            endif
+          enddo
+
+          max_ptj=-1.0d0
+          do i=1,nexternal
+            max_ptj = max(max_ptj,ptj(i))
+          enddo
+          if (max_ptj.lt.0.0d0) then
+            bias_weight = 1.0d0
+            return
+          endif
+
+          bias_weight = (max_ptj/ptj_bias_target_ptj)
+     &                                      **ptj_bias_enhancement_power
+
+          return
+
+      end subroutine bias_wgt

=== modified file 'Template/LO/Source/PDF/makefile'
--- Template/LO/Source/PDF/makefile	2014-09-05 11:30:25 +0000
+++ Template/LO/Source/PDF/makefile	2016-11-23 08:16:42 +0000
@@ -11,7 +11,11 @@
 include ../make_opts
 
 ifdef lhapdf
-  PDF         = pdfwrap_lhapdf.o pdf_lhapdf.o pdg2pdf_lhapdf.o opendata.o
+  ifeq ($(lhapdfversion),5)
+    PDF         = pdfwrap_lhapdf.o pdf_lhapdf.o pdg2pdf_lhapdf.o opendata.o  PhotonFlux.o
+  else
+    PDF         = pdfwrap_lhapdf.o pdf_lhapdf6.o pdg2pdf_lhapdf6.o opendata.o PhotonFlux.o
+  endif
 else
   PDF         = Ctq4Fn.o Ctq5Par.o Ctq5Pdf.o Partonx5.o Ctq6Pdf.o cteq3.o \
                 mrs98.o mrs98lo.o mrs98ht.o mrs99.o mrst2001.o mrst2002.o jeppe02.o \

=== modified file 'Template/LO/Source/PDF/pdf_lhapdf.f'
--- Template/LO/Source/PDF/pdf_lhapdf.f	2014-06-19 11:20:33 +0000
+++ Template/LO/Source/PDF/pdf_lhapdf.f	2016-11-23 08:16:42 +0000
@@ -11,6 +11,7 @@
       INTEGER IH,I
       double precision  photon
       LOGICAL has_photon
+      double precision epa_electron,epa_proton
 C
 C     Include
 C
@@ -28,11 +29,11 @@
          do i=-6,6
             pdf(i)=f(i)/x
          enddo
+      elseif(ih .eq. 2) then  !from a proton without breaking
+          pdf(7) = epa_proton(x, q * q)
       else
          write (*,*) 'beam type not supported in lhadpf'
-         do i=-6,6
-            pdf(i)=0d0
-         enddo
+         stop 1
       endif
 
       return	

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

=== added file 'Template/LO/Source/PDF/pdg2pdf_lhapdf6.f'
--- Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2016-11-23 08:16:42 +0000
@@ -0,0 +1,120 @@
+      double precision function pdg2pdf(ih,ipdg,x,xmu)
+c***************************************************************************
+c     Based on pdf.f, wrapper for calling the pdf of MCFM
+c***************************************************************************
+      implicit none
+c
+c     Arguments
+c
+      DOUBLE  PRECISION x,xmu
+      INTEGER IH,ipdg
+C
+C     Include
+C
+      include 'pdf.inc'
+C      
+      integer i,j,ihlast(20),ipart,iporg,ireuse,imemlast(20),iset,imem
+     &     ,i_replace,ii,ipartlast(20)
+      double precision xlast(20),xmulast(20),pdflast(20)
+      double precision epa_proton
+      save ihlast,xlast,xmulast,pdflast,imemlast,ipartlast
+      data ihlast/20*-99/
+      data ipartlast/20*-99/
+      data xlast/20*-99d9/
+      data xmulast/20*-99d9/
+      data pdflast/20*-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
+c     instead of stopping the code, as this might accidentally happen.
+      if (x.eq.0d0) then
+         pdg2pdf=0d0
+         return
+      elseif (x.lt.0d0 .or. x.gt.1d0) then
+         write (*,*) 'PDF not supported for Bjorken x ', x
+         open(unit=26,file='../../../error',status='unknown')
+         write(26,*) 'Error: PDF not supported for Bjorken x ',x
+         stop 1
+      endif
+
+      ipart=ipdg
+      if(iabs(ipart).eq.21) ipart=0
+      if(iabs(ipart).eq.22) ipart=7
+      iporg=ipart
+
+c     This will be called for any PDG code, but we only support up to 7
+      if(iabs(ipart).gt.7)then
+         write(*,*) 'PDF not supported for pdg ',ipdg
+         write(*,*) 'For lepton colliders, please set the lpp* '//
+     $    'variables to 0 in the run_card'  
+         open(unit=26,file='../../../error',status='unknown')
+         write(26,*) 'Error: PDF not supported for pdg ',ipdg
+         stop 1
+      endif
+
+c     Determine the iset used in lhapdf
+      call getnset(iset)
+      if (iset.ne.1) then
+         write (*,*) 'PDF not supported for Bjorken x ', x
+         open(unit=26,file='../../../error',status='unknown')
+         write(26,*) 'Error: PDF not supported for Bjorken x ',x
+         stop 1
+      endif
+
+c     Determine the member of the set (function of lhapdf)
+      call getnmem(iset,imem)
+
+      ireuse = 0
+      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 (ipart.eq.ipartlast(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
+         endif
+         ii=ii-1
+         if (ii.eq.0) ii=ii+20
+      enddo
+
+c     Reuse previous result, if possible
+      if (ireuse.gt.0) then
+         if (pdflast(ireuse).ne.-99d9) then
+            pdg2pdf=pdflast(ireuse)
+            return 
+         endif
+      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
+      if(ih.eq.1) then
+         call evolvepart(ipart,x,xmu,pdg2pdf)
+         pdg2pdf=pdg2pdf/x
+      else if(ih.eq.2) then ! photon from a proton without breaking
+          pdg2pdf = epa_proton(x,xmu*xmu)
+      else
+         write (*,*) 'beam type not supported in lhadpf'
+         stop 1
+      endif
+      pdflast(i_replace)=pdg2pdf
+      xlast(i_replace)=x
+      xmulast(i_replace)=xmu
+      ihlast(i_replace)=ih
+      imemlast(i_replace)=imem
+c
+      return
+      end
+

=== modified file 'Template/LO/Source/cuts.inc'
--- Template/LO/Source/cuts.inc	2015-04-09 01:31:09 +0000
+++ Template/LO/Source/cuts.inc	2016-11-23 08:16:42 +0000
@@ -88,6 +88,10 @@
       REAL*8 KT_DURHAM, D_PARAMETER
       LOGICAL DO_KT_DURHAM
       COMMON /JET_MEASURE_CUTS/ KT_DURHAM, D_PARAMETER
-      
-      
-      
\ No newline at end of file
+
+C COMMON BLOCK FOR PYTHIA JET MEASURE CUT
+      REAL*8 PT_LUND
+      COMMON /PYTHIA_MEASURE_CUTS/ PT_LUND
+      
+      
+      

=== modified file 'Template/LO/Source/dsample.f'
--- Template/LO/Source/dsample.f	2015-06-08 16:01:11 +0000
+++ Template/LO/Source/dsample.f	2016-11-23 08:16:42 +0000
@@ -1608,7 +1608,7 @@
       if (iteration .eq. cur_it) then
 c        Add the current point to the DiscreteSamplerGrid
          call add_entry_to_discrete_dimensions(wgt)
-         if (kn.eq.0) then
+         if (kn.eq.0.and.(iteration.eq.1.or.use_cut.eq.-2)) then
             ! ensure that all cumulative variable are at zero (usefull for reset)
             twgt1 = 0d0         !
             iavg = 0            !Vars for averging to increase err estimate
@@ -2156,7 +2156,7 @@
 c     $                 access='append',err=129)
 c                  write(22, 80) real(tmean), real(tsigma), real(chi2)
 c 129              close(22)
-                  tsigma = tsigma*sqrt(chi2) !This gives the 68% confidence cross section
+                  tsigma = tsigma*sqrt(max(0d0,chi2)) !This gives the 68% confidence cross section
                   cur_it = itm+20
                   return
                endif

=== modified file 'Template/LO/Source/genps.inc'
--- Template/LO/Source/genps.inc	2014-11-24 00:33:21 +0000
+++ Template/LO/Source/genps.inc	2016-11-23 08:16:42 +0000
@@ -30,9 +30,9 @@
 c	Parameters for helicity sums in matrixN.f
 c*************************************************************************
       REAL*8 LIMHEL
-      PARAMETER(LIMHEL=1e-6)
+      PARAMETER(LIMHEL=1e-8)
       INTEGER MAXTRIES
-      PARAMETER(MAXTRIES=10)
+      PARAMETER(MAXTRIES=25)
 C     To pass the helicity configuration chosen by the DiscreteSampler to 
 C     matrix<i>.f
       double precision hel_jacobian

=== modified file 'Template/LO/Source/kin_functions.f'
--- Template/LO/Source/kin_functions.f	2014-09-11 15:22:26 +0000
+++ Template/LO/Source/kin_functions.f	2016-11-23 08:16:42 +0000
@@ -127,7 +127,11 @@
 c      pm=dsqrt(p(1)**2+p(2)**2+p(3)**2)
       
       pm = p(0)
-      rap = .5d0*dlog((pm+p(3))/(pm-p(3)))+cm_rap
+      if (pm.gt.abs(p(3))) then
+        rap = .5d0*dlog((pm+p(3))/(pm-p(3)))+cm_rap
+      else
+        rap = -1d99
+      endif
       end
       DOUBLE PRECISION  FUNCTION rap2(p)
 c************************************************************************

=== added file 'Template/LO/Source/lhe_event_infos.inc'
--- Template/LO/Source/lhe_event_infos.inc	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/lhe_event_infos.inc	2016-11-23 08:16:42 +0000
@@ -0,0 +1,16 @@
+      integer jpart(7,-nexternal+3:2*nexternal-3)
+      double precision pb(0:4,-nexternal+3:2*nexternal-3)
+      integer isym(nexternal,99),jsym, npart
+      double precision sscale,aaqcd,aaqed
+      character*1000 buff
+      character*(s_bufflen) s_buff(7)
+      integer nclus
+      character*(clus_bufflen) buffclus(nexternal)
+	  character*(maxEventLength) event_record
+      logical AlreadySetInBiasModule
+
+	  common/to_lhe_event_info/jpart,pb,s_buff,buff,nclus,buffclus,event_record,
+     &  sscale,aaqcd,aaqed,isym,jsym,npart,AlreadySetInBiasModule
+
+      integer ngroup
+      common/to_group/ngroup

=== modified file 'Template/LO/Source/make_opts'
--- Template/LO/Source/make_opts	2014-06-03 07:42:21 +0000
+++ Template/LO/Source/make_opts	2016-11-23 08:16:42 +0000
@@ -1,59 +1,103 @@
-# Environment variables
-
+DEFAULT_F2PY_COMPILER=f2py
+DEFAULT_F_COMPILER=gfortran
+MACFLAG=-mmacosx-version-min=10.7
+DEFAULT_CPP_COMPILER=clang
+MG5AMC_VERSION=SpecifiedByMG5aMCAtRunTime
+STDLIB=-lc++
+PYTHIA8_PATH=NotInstalled
+STDLIB_FLAG=-stdlib=libc++
+#end_of_make_opts_variables
+
+BIASLIBDIR=../../../lib/
+BIASLIBRARY=libbias.$(libext)
+
+# Rest of the makefile
 ifeq ($(origin FFLAGS),undefined)
-  FFLAGS= -O -w -fbounds-check -fPIC
-  #FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
-endif
+FFLAGS= -O -w -fbounds-check -fPIC
+#FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
+endif
+
+
+# REMOVE MACFLAG IF NOT ON MAC OR FOR F2PY
+UNAME := $(shell uname -s)
+ifdef f2pymode
+MACFLAG=
+else
+ifneq ($(UNAME), Darwin)
+MACFLAG=
+endif
+endif
+
+
+ifeq ($(origin CXXFLAGS),undefined)
+CXXFLAGS= -O $(STDLIB_FLAG) $(MACFLAG)
+endif
+
+ifeq ($(origin CFLAGS),undefined)
+CFLAGS= -O $(STDLIB_FLAG) $(MACFLAG)
+endif
+
+# Increase the number of allowed charcters in a Fortran line
 FFLAGS+= -ffixed-line-length-132
 
 # Set FC unless it's defined by an environment variable
 ifeq ($(origin FC),default)
-#  FC=gfortran
-   FC=gfortran
+FC=$(DEFAULT_F_COMPILER)
+endif
+ifeq ($(origin F2PY), undefined)
+F2PY=$(DEFAULT_F2PY_COMPILER)
+endif
+
+UNAME := $(shell uname -s)
+ifeq ($(origin LDFLAGS), undefined)
+LDFLAGS=$(STDLIB) $(MACFLAG)
 endif
 
 # Options: dynamic, lhapdf
-
-UNAME := $(shell uname)
-
 # Option dynamic
 
 ifeq ($(UNAME), Darwin)
-  dylibext=dylib
+dylibext=dylib
 else
-  dylibext=so  
+dylibext=so
 endif
 
 ifdef dynamic
-  ifeq ($(UNAME), Darwin)
-    libext=dylib
-    FFLAGS+= -fno-common
-    LDFLAGS += -bundle
-    define CREATELIB
-      $(FC) -dynamiclib -undefined dynamic_lookup -o $(1) $(2)
-    endef
-  else
-    libext=so
-    FFLAGS+= -fPIC
-    LDFLAGS += -shared
-    define CREATELIB
-      $(FC) $(FFLAGS) $(LDFLAGS) -o $(1) $(2)
-    endef
-  endif
-else
-  libext=a
-  define CREATELIB
-    $(AR) cru $(1) $(2)
-    ranlib $(1)
-  endef
+ifeq ($(UNAME), Darwin)
+libext=dylib
+FFLAGS+= -fno-common
+LDFLAGS += -bundle
+define CREATELIB
+$(FC) -dynamiclib -undefined dynamic_lookup -o $(1) $(2)
+endef
+else
+libext=so
+FFLAGS+= -fPIC
+LDFLAGS += -shared
+define CREATELIB
+$(FC) $(FFLAGS) $(LDFLAGS) -o $(1) $(2)
+endef
+endif
+else
+libext=a
+define CREATELIB
+$(AR) cru $(1) $(2)
+ranlib $(1)
+endef
 endif
 
 # Option lhapdf
 
-ifdef lhapdf
-  alfas_functions=alfas_functions_lhapdf
-  lhapdf=-lLHAPDF -lstdc++
+ifneq ($(lhapdf),)
+CXXFLAGS += $(shell $(lhapdf) --cppflags)
+alfas_functions=alfas_functions_lhapdf
+llhapdf+= -lLHAPDF
 else
-  alfas_functions=alfas_functions
-  lhapdf=
+alfas_functions=alfas_functions
+llhapdf=
 endif
+
+# Helper function to check MG5 version
+define CHECK_MG5AMC_VERSION
+python -c 'import re; from distutils.version import StrictVersion; print StrictVersion("$(MG5AMC_VERSION)") >= StrictVersion("$(1)") if re.match("^[\d\.]+$$","$(MG5AMC_VERSION)") else True;'
+endef
\ No newline at end of file

=== modified file 'Template/LO/Source/run.inc'
--- Template/LO/Source/run.inc	2015-04-02 22:56:24 +0000
+++ Template/LO/Source/run.inc	2016-11-23 08:16:42 +0000
@@ -63,3 +63,11 @@
 C
       logical MC_grouped_subproc
       common/to_MC_grouped_subproc/MC_grouped_subproc
+
+C
+C Controls what are the PDGs included in the CKKWl merging procedure, i.e. what
+C are the PDGs subject to the ktdurham cut
+C
+      integer pdgs_for_merging_cut(0:1000)
+      common/TO_MERGE/pdgs_for_merging_cut
+

=== modified file 'Template/LO/Source/rw_events.f'
--- Template/LO/Source/rw_events.f	2014-10-07 13:31:57 +0000
+++ Template/LO/Source/rw_events.f	2016-11-23 08:16:42 +0000
@@ -30,7 +30,6 @@
       character*(s_bufflen) s_buff(*)
       integer nclus
       character*(clus_bufflen) buffclus(*)
-
 c
 c     Local
 c
@@ -46,6 +45,10 @@
 
       data lun_ban/37/
       data banner_open/.false./
+
+      double precision bias_weight
+      logical impact_xsec
+      common/bias/bias_weight,impact_xsec
 c-----
 c  Begin Code
 c-----     
@@ -80,6 +83,21 @@
          buff=''
         endif
       endif
+
+c     Reading the bias weight (if present)
+      read(lun,'(a300)',end=99,err=99) buftmp
+      if(buftmp(1:6).ne.'<rwgt>') then
+         backspace(lun)
+         bias_weight = 1.0d0
+      else
+         do while(buftmp(1:7).ne.'</rwgt>')
+           read(lun,'(a300)',end=99,err=99) buftmp
+           if (buftmp(1:16).eq."<wgt id='bias'> ") then
+              read(buftmp(17:31),'(1e15.7)') bias_weight
+           endif
+         enddo
+      endif
+
 c     Systematics info
       read(lun,'(a)',end=99,err=99) s_buff(1)
       if(s_buff(1).ne.'<mgrwt>') then
@@ -114,9 +132,103 @@
  55   format(i3,5e19.11)         
       end
 
+      subroutine write_event_to_stream(evt_record,P,wgt,nexternal,ic,
+     &     ievent,scale,aqcd, aqed,buff,u_syst,s_buff,nclus,buffclus)
+c********************************************************************
+C     This an *exact* copy of write_event, except that it writes it 
+C     to a character array argument as opposed to an I/O stream.
+c********************************************************************
+      implicit none
+
+      include 'maxparticles.inc'
+      include 'run_config.inc'
+c
+c     parameters
+c
+      double precision pi
+      parameter (pi = 3.1415926d0)
+c
+c     Arguments
+c      
+      character*(maxEventLength) evt_record
+      integer ievent
+      integer nexternal, ic(7,*)
+      double precision P(0:4,*),wgt
+      double precision aqcd, aqed, scale
+      character*1000 buff
+      logical u_syst
+      character*(s_bufflen) s_buff(*)
+      integer nclus
+      character*(clus_bufflen) buffclus(*)
+c
+c     Local
+c
+      integer i,j,k
+      character*(maxEventLength) largeBuff
+c
+c     Global
+c
+      double precision bias_weight
+      logical impact_xsec
+      common/bias/bias_weight,impact_xsec
+
+c-----
+c  Begin Code
+c-----     
+c      aqed= gal(1)*gal(1)/4d0/pi
+c      aqcd = g*g/4d0/pi
+      write(largeBuff,'(a)') '<event>'
+      evt_record=trim(evt_record)//trim(largeBuff)
+      write(largeBuff,'(i2,i4,e16.7e3,3e15.7)') nexternal,ievent,wgt,scale,
+     $                                   aqed,aqcd
+      evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      do i=1,nexternal
+         write(largeBuff,51) ic(1,i),ic(6,i),(ic(j,i),j=2,5),
+     $     (p(j,i),j=1,3),p(0,i),p(4,i),0.,real(ic(7,i))
+         evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      enddo
+      if(buff(1:7).eq.'<scales') then
+        write(largeBuff,'(a)') buff(1:len_trim(buff))
+        evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      endif
+      if(buff(1:1).eq.'#') then
+        write(largeBuff,'(a)') buff(1:len_trim(buff))
+        evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      endif
+      if(.not.impact_xsec) then
+          write(largeBuff,'(a)') '<rwgt>'
+          evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+          write(largeBuff,'(a16,1e15.7,a6)') "<wgt id='bias'> ",
+     $                                              bias_weight,"</wgt>"
+          evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+          write(largeBuff,'(a)') '</rwgt>'
+          evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      endif 
+      if(u_syst)then
+         do i=1,7
+            write(largeBuff,'(a)') s_buff(i)(1:len_trim(s_buff(i)))
+            evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+         enddo
+      endif
+      do i=1,nclus
+         write(largeBuff,'(a)') buffclus(i)(1:len_trim(buffclus(i)))
+         evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      enddo
+      write(largeBuff,'(a)') '</event>'
+      evt_record=trim(evt_record)//CHAR(13)//CHAR(10)//trim(largeBuff)
+      return
+ 51   format(i9,5i5,5e19.11,f3.0,f4.0)
+      end
+
+
       subroutine write_event(lun,P,wgt,nexternal,ic,ievent,scale,aqcd,
      $     aqed,buff,u_syst,s_buff,nclus,buffclus)
 c********************************************************************
+c
+c /!\ When making changes to this subroutine, make sure to accordingly
+c     update write_event_to_stream
+c
+c********************************************************************
 c     Writes one event from data file #lun according to LesHouches
 c     ic(1,*) = Particle ID
 c     ic(2.*) = Mothup(1)
@@ -154,6 +266,9 @@
 c
 c     Global
 c
+      double precision bias_weight
+      logical impact_xsec
+      common/bias/bias_weight,impact_xsec
 
 c-----
 c  Begin Code
@@ -162,13 +277,20 @@
 c      aqcd = g*g/4d0/pi
 
       write(lun,'(a)') '<event>'
-      write(lun,'(i2,i4,4e15.7)') nexternal,ievent,wgt,scale,aqed,aqcd
+      write(lun,'(i2,i4,e16.7e3,3e15.7)') nexternal,ievent,wgt,scale,aqed,aqcd
+      write(*,'(i2,i4,e15.7e3,3e15.7)') nexternal,ievent,wgt,scale,aqed,aqcd
       do i=1,nexternal
          write(lun,51) ic(1,i),ic(6,i),(ic(j,i),j=2,5),
      $     (p(j,i),j=1,3),p(0,i),p(4,i),0.,real(ic(7,i))
       enddo
       if(buff(1:7).eq.'<scales') write(lun,'(a)') buff(1:len_trim(buff))
-      if(buff(1:1).eq.'#') write(lun,'(a)') buff(1:len_trim(buff))
+      if(buff(1:1).eq.'#') write(lun,'(a)') buff(1:len_trim(buff))      
+      if(.not.impact_xsec) then
+          write(lun,'(a)') '<rwgt>'
+          write(lun,'(a16,1e15.7,a6)') "<wgt id='bias'> ",bias_weight,
+     $                                                          "</wgt>"
+          write(lun,'(a)') '</rwgt>'
+      endif 
       if(u_syst)then
          do i=1,7
             write(lun,'(a)') s_buff(i)(1:len_trim(s_buff(i)))

=== modified file 'Template/LO/Source/setrun.f'
--- Template/LO/Source/setrun.f	2015-04-09 01:31:09 +0000
+++ Template/LO/Source/setrun.f	2016-11-23 08:16:42 +0000
@@ -199,7 +199,7 @@
      $   10042,
      $   246800,
      $   247000,
-     $   244600/
+     $   244800/
 
 
       if(pdfin.eq."lhapdf") then

=== modified file 'Template/LO/SubProcesses/cuts.f'
--- Template/LO/SubProcesses/cuts.f	2015-04-07 22:51:03 +0000
+++ Template/LO/SubProcesses/cuts.f	2016-11-23 08:16:42 +0000
@@ -1,1316 +1,1750 @@
-      logical function pass_point(p)
-c************************************************************************
-c     This function is called from sample to see if it needs to 
-c     bother calculating the weight from all the different conficurations
-c     You can either just return true, or have it call passcuts
-c************************************************************************
-      implicit none
-c
-c     Arguments
-c
-      double precision p
-c
-c     External
-c
-      logical passcuts
-      external passcuts
-c-----
-c  Begin Code
-c-----
-      pass_point = .true.
-c      pass_point = passcuts(p)
-      end
-C 
-      LOGICAL FUNCTION PASSCUTS(P)
-C**************************************************************************
-C     INPUT:
-C            P(0:3,1)           MOMENTUM OF INCOMING PARTON
-C            P(0:3,2)           MOMENTUM OF INCOMING PARTON
-C            P(0:3,3)           MOMENTUM OF ...
-C            ALL MOMENTA ARE IN THE REST FRAME!!
-C            COMMON/JETCUTS/   CUTS ON JETS
-C     OUTPUT:
-C            TRUE IF EVENTS PASSES ALL CUTS LISTED
-C**************************************************************************
-      IMPLICIT NONE
-c
-c     Constants
-c
-      include 'genps.inc'
-      include 'nexternal.inc'
-C
-C     ARGUMENTS
-C
-      REAL*8 P(0:3,nexternal)
-
-C
-C     LOCAL
-C
-      LOGICAL FIRSTTIME,FIRSTTIME2,pass_bw,notgood,good,foundheavy
-      LOGICAL DEBUG
-      integer i,j,njets,nheavyjets,nleptons,hardj1,hardj2
-      REAL*8 XVAR,ptmax1,ptmax2,htj,tmp,inclht
-      real*8 ptemp(0:3), ptemp2(0:3)
-      character*20 formstr
-C
-C     PARAMETERS
-C
-      real*8 PI
-      parameter( PI = 3.14159265358979323846d0 )
-C
-C     EXTERNAL
-C
-      REAL*8 R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,PtDot
-      logical cut_bw,setclscales
-      external R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,cut_bw,setclscales,PtDot
-C
-C     GLOBAL
-C
-      include 'run.inc'
-      include 'cuts.inc'
-
-      double precision ptjet(nexternal)
-      double precision ptheavyjet(nexternal)
-      double precision ptlepton(nexternal)
-      double precision temp
-
-C VARIABLES TO SPECIFY JETS
-      DOUBLE PRECISION PJET(NEXTERNAL,0:3)
-      DOUBLE PRECISION PTMIN
-      DOUBLE PRECISION PT1,PT2
-      INTEGER K,J1,J2
-
-C VARIABLES FOR KT CUT
-      DOUBLE PRECISION PTNOW,COSTH,PABS1,PABS2
-      DOUBLE PRECISION ETA1,ETA2,COSH_DETA,COS_DPHI,KT1SQ,KT2SQ, DPHI
-
-      double precision etmin(nincoming+1:nexternal),etamax(nincoming+1:nexternal)
-      double precision emin(nincoming+1:nexternal)
-      double precision                    r2min(nincoming+1:nexternal,nincoming+1:nexternal)
-      double precision s_min(nexternal,nexternal)
-      double precision etmax(nincoming+1:nexternal),etamin(nincoming+1:nexternal)
-      double precision emax(nincoming+1:nexternal)
-      double precision r2max(nincoming+1:nexternal,nincoming+1:nexternal)
-      double precision s_max(nexternal,nexternal)
-      double precision ptll_min(nexternal,nexternal),ptll_max(nexternal,nexternal)
-      double precision inclHtmin,inclHtmax
-      common/to_cuts/  etmin, emin, etamax, r2min, s_min,
-     $     etmax, emax, etamin, r2max, s_max, ptll_min, ptll_max, inclHtmin,inclHtmax
-
-      double precision ptjmin4(4),ptjmax4(4),htjmin4(2:4),htjmax4(2:4)
-      logical jetor
-      common/to_jet_cuts/ ptjmin4,ptjmax4,htjmin4,htjmax4,jetor
-
-      double precision ptlmin4(4),ptlmax4(4)
-      common/to_lepton_cuts/ ptlmin4,ptlmax4
-
-c
-c     Special cuts
-c
-
-      integer        lbw(0:nexternal)  !Use of B.W.
-      common /to_BW/ lbw
-C
-C     SPECIAL CUTS
-C
-      LOGICAL  IS_A_J(NEXTERNAL),IS_A_L(NEXTERNAL)
-      LOGICAL  IS_A_B(NEXTERNAL),IS_A_A(NEXTERNAL),IS_A_ONIUM(NEXTERNAL)
-      LOGICAL  IS_A_NU(NEXTERNAL),IS_HEAVY(NEXTERNAL)
-      COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,
-     . IS_A_ONIUM
-C
-C     Keep track of whether cuts already calculated for this event
-C
-      LOGICAL CUTSDONE,CUTSPASSED
-      COMMON/TO_CUTSDONE/CUTSDONE,CUTSPASSED
-      DATA CUTSDONE,CUTSPASSED/.FALSE.,.FALSE./
-
-C $B$ MW_NEW_DEF $E$ !this is a tag for MadWeight
-
-      double precision xqcutij(nexternal,nexternal),xqcuti(nexternal)
-      common/to_xqcuts/xqcutij,xqcuti
-
-c jet cluster algorithm
-      integer nQCD !,NJET,JET(nexternal)
-c      double precision plab(0:3, nexternal)
-      double precision pQCD(0:3,nexternal)!,PJET(0:3,nexternal)
-c      double precision rfj,sycut,palg,fastjetdmerge
-c      integer njet_eta
-c     Photon isolation
-      integer nph,nem,nin
-      double precision ptg,chi_gamma_iso,iso_getdrv40
-      double precision Etsum(0:nexternal)
-      real drlist(nexternal)
-      double precision pgamma(0:3,nexternal),pem(0:3,nexternal)
-      logical alliso
-C     Sort array of results: ismode>0 for real, isway=0 for ascending order
-      integer ismode,isway,izero,isorted(nexternal)
-      parameter (ismode=1)
-      parameter (isway=0)
-      parameter (izero=0)
-
-      include 'coupl.inc'
-C
-C
-c
-      DATA FIRSTTIME,FIRSTTIME2/.TRUE.,.TRUE./
-
-c put momenta in common block for couplings.f
-      double precision pp(0:3,max_particles)
-      common /momenta_pp/pp
-
-      DATA DEBUG/.FALSE./
-
-C-----
-C  BEGIN CODE
-C-----
-
-
-
-      PASSCUTS=.TRUE.             !EVENT IS OK UNLESS OTHERWISE CHANGED
-      IF (FIRSTTIME) THEN
-         FIRSTTIME=.FALSE.
-c      Preparation for reweighting by setting up clustering by diagrams
-         call initcluster()
-c
-c
-         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'i8)'
-         write(*,formstr) 'Particle',(i,i=nincoming+1,nexternal)
-         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'f8.1)'
-         write(*,formstr) 'Et >',(etmin(i),i=nincoming+1,nexternal)
-         write(*,formstr) 'E >',(emin(i),i=nincoming+1,nexternal)
-         write(*,formstr) 'Eta <',(etamax(i),i=nincoming+1,nexternal)
-         write(*,formstr) 'xqcut: ',(xqcuti(i),i=nincoming+1,nexternal)
-         write(formstr,'(a,i2.2,a)')'(a,i2,a,',nexternal,'f8.1)'
-         do j=nincoming+1,nexternal-1
-            write(*,formstr) 'd R #',j,'  >',(-0.0,i=nincoming+1,j),
-     &           (r2min(i,j),i=j+1,nexternal)
-            do i=j+1,nexternal
-               r2min(i,j)=r2min(i,j)*dabs(r2min(i,j))    !Since r2 returns distance squared
-               r2max(i,j)=r2max(i,j)*dabs(r2max(i,j))
-            enddo
-         enddo
-         do j=nincoming+1,nexternal-1
-            write(*,formstr) 's min #',j,'>',
-     &           (s_min(i,j),i=nincoming+1,nexternal)
-         enddo
-         do j=nincoming+1,nexternal-1
-            write(*,formstr) 'xqcutij #',j,'>',
-     &           (xqcutij(i,j),i=nincoming+1,nexternal)
-         enddo
-
-cc
-cc     Set the strong coupling
-cc
-c         call set_ren_scale(P,scale)
-c
-cc     Check that the user funtions for setting the scales
-cc     have been edited if the choice of an event-by-event
-cc     scale choice has been made 
-c
-c         if(.not.fixed_ren_scale) then
-c            if(scale.eq.0d0) then
-c               write(6,*) 
-c               write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
-c               write(6,*) ' Dynamical renormalization scale choice '
-c               write(6,*) ' selected but user subroutine' 
-c               write(6,*) ' set_ren_scale not edited in file:setpara.f'
-c               write(6,*) ' Switching to a fixed_ren_scale choice'
-c               write(6,*) ' with scale=zmass'
-c               scale=91.2d0
-c               write(6,*) 'scale=',scale
-c               fixed_ren_scale=.true.
-c               call set_ren_scale(P,scale)
-c            endif
-c         endif
-         
-c         if(.not.fixed_fac_scale) then
-c            call set_fac_scale(P,q2fact)
-c            if(q2fact(1).eq.0d0.or.q2fact(2).eq.0d0) then
-c               write(6,*) 
-c               write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
-c               write(6,*) ' Dynamical renormalization scale choice '
-c               write(6,*) ' selected but user subroutine' 
-c               write(6,*) ' set_fac_scale not edited in file:setpara.f'
-c               write(6,*) ' Switching to a fixed_fac_scale choice'
-c               write(6,*) ' with q2fact(i)=zmass**2'
-c               fixed_fac_scale=.true.
-c               q2fact(1)=91.2d0**2
-c               q2fact(2)=91.2d0**2
-c               write(6,*) 'scales=',q2fact(1),q2fact(2)
-c            endif
-c         endif
-
-         if(fixed_ren_scale) then
-            G = SQRT(4d0*PI*ALPHAS(scale))
-            call update_as_param()
-         endif
-
-c     Put momenta in the common block to zero to start
-         do i=0,3
-            do j=1,max_particles
-               pp(i,j) = 0d0
-            enddo
-         enddo
-         
-      ENDIF ! IF FIRSTTIME
-
-      IF (CUTSDONE) THEN
-         PASSCUTS=CUTSPASSED
-         RETURN
-      ENDIF
-      CUTSDONE=.TRUE.
-c      CUTSPASSED=.FALSE.
-
-c
-c     Make sure have reasonable 4-momenta
-c
-      if (p(0,1) .le. 0d0) then
-         passcuts=.false.
-         return
-      endif
-
-c     Also make sure there's no INF or NAN
-      do i=1,nexternal
-         do j=0,3
-            if(p(j,i).gt.1d32.or.p(j,i).ne.p(j,i))then
-               passcuts=.false.
-               return
-            endif
-         enddo
-      enddo
-      
-c
-c     Limit S_hat
-c
-c      if (x1*x2*stot .gt. 500**2) then
-c         passcuts=.false.
-c         return
-c      endif
-
-C $B$ DESACTIVATE_CUT $E$ !This is a tag for MadWeight
-
-      if(debug) write (*,*) '============================='
-      if(debug) write (*,*) ' EVENT STARTS TO BE CHECKED  '
-      if(debug) write (*,*) '============================='
-c     
-c     p_t min & max cuts
-c     
-      do i=nincoming+1,nexternal
-         if(debug) write (*,*) 'pt(',i,')=',pt(p(0,i)),'   ',etmin(i),
-     $        ':',etmax(i)
-         notgood=(pt(p(0,i)) .lt. etmin(i)).or.
-     &        (etmax(i).ge.0d0.and.pt(p(0,i)) .gt. etmax(i))
-         if (notgood) then
-            if(debug) write (*,*) i,' -> fails'
-            passcuts=.false.
-            return
-         endif
-      enddo
-c
-c    missing ET min & max cut + Invariant mass of leptons and neutrino 
-c    nb: missing Et defined as the vector sum over the neutrino's pt
-c
-c-- reset ptemp(0:3)
-      do j=0,3
-         ptemp(j)=0 ! for the neutrino
-         ptemp2(j)=0 ! for the leptons
-      enddo
-c-  sum over the momenta
-      do i=nincoming+1,nexternal
-         if(is_a_nu(i)) then            
-         if(debug) write (*,*) i,' -> neutrino '
-            do j=0,3
-               ptemp(j)=ptemp(j)+p(j,i)
-            enddo
-         elseif(is_a_l(i)) then            
-         if(debug) write (*,*) i,' -> lepton '
-            do j=0,3
-               ptemp2(j)=ptemp2(j)+p(j,i)
-            enddo
-         endif
-
-      enddo
-c-  check the et
-      if(debug.and.ptemp(0).eq.0d0) write (*,*) 'No et miss in event'
-      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Et miss =',pt(ptemp(0)),'   ',misset,':',missetmax
-      if(debug.and.ptemp2(0).eq.0d0) write (*,*) 'No leptons in event'
-      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Energy of leptons =',pt(ptemp2(0))
-      if(ptemp(0).gt.0d0) then
-         notgood=(pt(ptemp(0)) .lt. misset).or.
-     &        (missetmax.ge.0d0.and.pt(ptemp(0)) .gt. missetmax)
-         if (notgood) then
-            if(debug) write (*,*) ' missing et cut -> fails'
-            passcuts=.false.
-            return
-         endif
-      endif
-      if (mmnl.gt.0d0.or.mmnlmax.ge.0d0)then
-         if(dsqrt(SumDot(ptemp,ptemp2,1d0)).lt.mmnl.or.mmnlmax.ge.0d0.and.dsqrt(SumDot(ptemp, ptemp2,1d0)).gt.mmnlmax) then
-            if(debug) write (*,*) 'lepton invariant mass -> fails'
-            passcuts=.false.
-            return
-         endif
-      endif
-c
-c     pt cut on heavy particles
-c     gives min(pt) for (at least) one heavy particle
-c
-      if(ptheavy.gt.0d0)then
-         passcuts=.false.
-         foundheavy=.false.
-         do i=nincoming+1,nexternal
-            if(is_heavy(i)) then            
-               if(debug) write (*,*) i,' -> heavy '
-               foundheavy=.true.
-               if(pt(p(0,i)).gt.ptheavy) passcuts=.true.
-            endif
-         enddo
-         
-         if(.not.passcuts.and.foundheavy)then
-            if(debug) write (*,*) ' heavy particle cut -> fails'
-            return
-         else
-            passcuts=.true.
-         endif
-      endif
-c     
-c     E min & max cuts
-c     
-      do i=nincoming+1,nexternal
-         if(debug) write (*,*) 'p(0,',i,')=',p(0,i),'   ',emin(i),':',emax(i)
-         notgood=(p(0,i) .le. emin(i)).or.
-     &        (emax(i).ge.0d0 .and. p(0,i) .gt. emax(i))
-         if (notgood) then
-            if(debug) write (*,*) i,' -> fails'
-            passcuts=.false.
-            return
-         endif
-      enddo
-c     
-c     Rapidity  min & max cuts
-c     
-      do i=nincoming+1,nexternal
-         if(debug) write (*,*) 'abs(rap(',i,'))=',abs(rap(p(0,i))),'   ',etamin(i),':',etamax(i)
-         notgood=(etamax(i).ge.0.and.abs(rap(p(0,i))) .gt. etamax(i)).or.
-     &        (abs(rap(p(0,i))) .lt. etamin(i))
-         if (notgood) then
-            if(debug) write (*,*) i,' -> fails'
-            passcuts=.false.
-            return
-         endif
-      enddo
-c     
-c     DeltaR min & max cuts
-c     
-      do i=nincoming+1,nexternal-1
-         do j=i+1,nexternal
-            if(debug) write (*,*) 'r2(',i, ',' ,j,')=',dsqrt(r2(p(0,i),p(0,j)))
-            if(debug) write (*,*) dsqrt(r2min(j,i)),dsqrt(r2max(j,i))
-            if(r2min(j,i).gt.0.or.r2max(j,i).ge.0d0) then
-               tmp=r2(p(0,i),p(0,j))
-               notgood=(tmp .lt. r2min(j,i)).or.
-     $              (r2max(j,i).ge.0d0 .and. tmp .gt. r2max(j,i))
-               if (notgood) then
-                  if(debug) write (*,*) i,j,' -> fails'
-                  passcuts=.false.
-                  return
-               endif
-            endif
-         enddo
-      enddo
-
-
-c     s-channel min & max pt of sum of 4-momenta
-c     
-      do i=nincoming+1,nexternal-1
-         do j=i+1,nexternal
-            if(debug)write (*,*) 'ptll(',i,',',j,')=',PtDot(p(0,i),p(0,j))
-            if(debug)write (*,*) ptll_min(j,i),ptll_max(j,i)
-            if(ptll_min(j,i).gt.0.or.ptll_max(j,i).ge.0d0) then
-               tmp=PtDot(p(0,i),p(0,j))
-               notgood=(tmp .lt. ptll_min(j,i).or.
-     $              ptll_max(j,i).ge.0d0 .and. tmp.gt.ptll_max(j,i))
-               if (notgood) then
-                  if(debug) write (*,*) i,j,' -> fails'
-                  passcuts=.false.
-                  return
-               endif
-            endif
-         enddo
-      enddo
-
-
-
-
-c     
-c     s-channel min & max invariant mass cuts
-c     
-      do i=nincoming+1,nexternal-1
-         do j=i+1,nexternal
-            if(debug) write (*,*) 's(',i,',',j,')=',Sumdot(p(0,i),p(0,j),+1d0)
-            if(debug) write (*,*) s_min(j,i),s_max(j,i)
-            if(s_min(j,i).gt.0.or.s_max(j,i).ge.0d0) then
-               tmp=SumDot(p(0,i),p(0,j),+1d0)
-               if(s_min(j,i).le.s_max(j,i) .or. s_max(j,i).lt.0d0)then
-                  notgood=(tmp .lt. s_min(j,i).or.
-     $                 s_max(j,i).ge.0d0 .and. tmp .gt. s_max(j,i)) 
-                  if (notgood) then
-                     if(debug) write (*,*) i,j,' -> fails'
-                     passcuts=.false.
-                     return
-                  endif
-               else
-                  notgood=(tmp .lt. s_min(j,i).and.tmp .gt. s_max(j,i)) 
-                  if (notgood) then
-                     if(debug) write (*,*) i,j,' -> fails'
-                     passcuts=.false.
-                     return
-                  endif
-               endif
-            endif
-         enddo
-      enddo
-C     $B$DESACTIVATE_BW_CUT$B$ This is a Tag for MadWeight
-c     
-c     B.W. phase space cuts
-c     
-      pass_bw=cut_bw(p)
-c     JA 4/8/11 always check pass_bw
-      if ( pass_bw.and..not.CUTSPASSED) then
-         passcuts=.false.
-         if(debug) write (*,*) ' pass_bw -> fails'