Added additional MadSpin configuration options. Fixed XML reweighting and...

Added additional MadSpin configuration options. Fixed XML reweighting and added fall-back to old-style reweighting for anything other than scale/PDF weights. 2017-10-26 13:09:37. (PowhegControl-00-03-04)

Generators/PowhegControl/python/algorithms/postprocessors/madspin.py

@@ -35,18 +35,14 @@ def __construct_inputs(input_LHE_events, process):
     """
     # Find insertion point for MadSpin header
     logger.info("Constructing MadSpin runcard header")
-    opening_string = LHE.get_opening_string(input_LHE_events)
-    closing_string = LHE.get_closing_string(input_LHE_events)
-    if opening_string.find("<header>") != -1:
-        pre_header = opening_string[: opening_string.find("<header>")]
-        post_header = opening_string[opening_string.find("<header>") + 8: opening_string.find("</header>") + 10]
-    else:
-        pre_header = opening_string[: opening_string.find("<init>")]
-        post_header = "</header>\n" + opening_string[opening_string.find("<init>"): opening_string.find("</init>") + 8]
+    pre_header  = "\n".join([elem for elem in [LHE.opening_tag(input_LHE_events), LHE.comment_block(input_LHE_events)] if elem])
+    header_contents = LHE.header_block(input_LHE_events).replace("<header>", "").replace("</header>", "").strip("\n")
+    post_header  = LHE.init_block(input_LHE_events)
+    postamble  = LHE.postamble(input_LHE_events)
 
     # Write events to LHE file with MadSpin information
     with open("madspin_LHE_input.lhe", "wb") as f_madspin_LHE:
-        f_madspin_LHE.write(pre_header)
+        f_madspin_LHE.write("{}\n".format(pre_header))
         f_madspin_LHE.write("<header>\n")
         f_madspin_LHE.write("<mgversion>\n")
         f_madspin_LHE.write("{}\n".format(os.environ["MADPATH"].split("/")[-1].split("v")[-1].split("_p")[0].replace("_", ".")))
@@ -58,12 +54,24 @@ def __construct_inputs(input_LHE_events, process):
         f_madspin_LHE.write("set gauge unitary\n")
         f_madspin_LHE.write("set complex_mass_scheme False\n")
         f_madspin_LHE.write("import model {}\n".format(process.MadSpin_model))
-        f_madspin_LHE.write("define p = g u c d s u~ c~ d~ s~\n")
-        f_madspin_LHE.write("define j = g u c d s u~ c~ d~ s~\n")
-        f_madspin_LHE.write("define l+ = e+ mu+ ta+\n")
-        f_madspin_LHE.write("define l- = e- mu- ta-\n")
-        f_madspin_LHE.write("define vl = ve vm vt\n")
-        f_madspin_LHE.write("define vl~ = ve~ vm~ vt~\n")
+        if int(process.MadSpin_nFlavours) == 4:
+            f_madspin_LHE.write("define p = g u c d s u~ c~ d~ s~\n")
+            f_madspin_LHE.write("define j = g u c d s u~ c~ d~ s~\n")
+        elif int(process.MadSpin_nFlavours) == 5:
+            f_madspin_LHE.write("define p = g u c d s b u~ c~ d~ s~ b~\n")
+            f_madspin_LHE.write("define j = g u c d s b u~ c~ d~ s~ b~\n")
+        else:
+            raise ValueError("'MadSpin_nFlavours' must be set to 4 or 5")
+        if process.MadSpin_taus_are_leptons:
+            f_madspin_LHE.write("define l+ = e+ mu+ ta+\n")
+            f_madspin_LHE.write("define l- = e- mu- ta-\n")
+            f_madspin_LHE.write("define vl = ve vm vt\n")
+            f_madspin_LHE.write("define vl~ = ve~ vm~ vt~\n")
+        else:
+            f_madspin_LHE.write("define l+ = e+ mu+\n")
+            f_madspin_LHE.write("define l- = e- mu-\n")
+            f_madspin_LHE.write("define vl = ve vm\n")
+            f_madspin_LHE.write("define vl~ = ve~ vm~\n")
         f_madspin_LHE.write("{}\n".format(process.MadSpin_process))
         f_madspin_LHE.write("output tchan\n")
         f_madspin_LHE.write("</mg5proccard>\n")
@@ -216,10 +224,12 @@ def __construct_inputs(input_LHE_events, process):
         f_madspin_LHE.write("#      PDG        Width\n")
         f_madspin_LHE.write("DECAY  82   0.000000e+00\n")
         f_madspin_LHE.write("</slha>\n")
-        f_madspin_LHE.write(post_header)
+        f_madspin_LHE.write("{}\n".format(header_contents))
+        f_madspin_LHE.write("</header>\n")
+        f_madspin_LHE.write("{}\n".format(post_header))
         for event in LHE.event_iterator(input_LHE_events):
             f_madspin_LHE.write(event)
-        f_madspin_LHE.write(closing_string)
+        f_madspin_LHE.write(postamble)
 
     # Rename LHE files
     shutil.move(input_LHE_events, "{}.undecayed".format(input_LHE_events))
@@ -254,32 +264,25 @@ def __prepare_outputs(input_LHE_events):
     subprocess.call("gunzip pwgevents_decayed.lhe.gz 2> /dev/null", shell=True)
     shutil.move("pwgevents_decayed.lhe", "{}.decayed".format(input_LHE_events))
 
-    # Combine headers
-    pwg_opening_string = LHE.get_opening_string("{}.undecayed".format(input_LHE_events))
-    pwg_closing_string = LHE.get_closing_string("{}.undecayed".format(input_LHE_events))
-    madspin_opening_string = LHE.get_opening_string("{}.decayed".format(input_LHE_events))
-
-    # Split Powheg headers to insert MadSpin information
-    if pwg_opening_string.find("<header>") != -1:
-        pwg_pre_header = pwg_opening_string[: pwg_opening_string.find("<header>") + 8]
-        pwg_post_header = pwg_opening_string[pwg_opening_string.find("</header>") + 10:]
-    else:
-        pwg_pre_header = pwg_opening_string[: pwg_opening_string.find("<init>")] + "<header>"
-        pwg_post_header = pwg_opening_string[pwg_opening_string.find("<init>"):]
-    madspin_header = madspin_opening_string[madspin_opening_string.find("<header>") + 8: madspin_opening_string.find("</header>") + 10]
+    # Get appropriate header sections from Powheg and MadSpin LHE
+    powheg_LHE, madspin_LHE = "{}.undecayed".format(input_LHE_events), "{}.decayed".format(input_LHE_events)
+    pre_header = "\n".join([elem for elem in [LHE.opening_tag(powheg_LHE), LHE.comment_block(powheg_LHE)] if elem])
+    header = LHE.header_block(madspin_LHE)
+    post_header = LHE.init_block(powheg_LHE)
+    postamble = LHE.postamble(powheg_LHE)
 
     # Write events to LHE file with MadSpin information
     with open("pwgevents_with_MadSpin.lhe", "wb") as f_combined_LHE:
-        f_combined_LHE.write(pwg_pre_header)
-        f_combined_LHE.write(madspin_header)
-        f_combined_LHE.write(pwg_post_header)
+        f_combined_LHE.write("{}\n".format(pre_header))
+        f_combined_LHE.write("{}\n".format(header))
+        f_combined_LHE.write("{}\n".format(post_header))
         for event in LHE.event_iterator("{}.decayed".format(input_LHE_events)):
             f_combined_LHE.write(event)
-        f_combined_LHE.write(pwg_closing_string)
+        f_combined_LHE.write(postamble)
 
     # Rename output file
     if os.path.isfile(input_LHE_events):
-        shutil.move(input_LHE_events, "madspin_LHE.input")
+        shutil.move(input_LHE_events, "{}.undecayed_ready_for_madspin".format(input_LHE_events))
     shutil.move("pwgevents_with_MadSpin.lhe", input_LHE_events)
 
     for LHE_tarball in list(set(glob.glob("*lhe*.gz") + glob.glob("*events*.gz"))):

Generators/PowhegControl/python/algorithms/postprocessors/nnlo_reweighter.py

@@ -60,7 +60,6 @@ def acquire_NNLO_inputs(aux_directory, NNLO_reweighting_inputs):
 def run_NNLO_executable(process, powheg_LHE_output, output_LHE_name):
     """! Run the NNLO executable."""
     # Run NNLOPS
-    # process.executable, process.NNLO_reweighting_inputs,
     NNLO_executable = process.executable
     if "nnlopsreweighter" in NNLO_executable:
         process.NNLO_weight_list = construct_NNLOPS_weight_list(process.NNLO_output_weights, powheg_LHE_output)
@@ -158,20 +157,15 @@ def reformat_NNLOPS_events(NNLO_weight_list, powheg_LHE_output, output_LHE_name)
                     f_output.write(output_line)
 
 
-def reformat_DYNNLO_events(powheg_LHE_output, output_LHE_name):
+def reformat_DYNNLO_events(powheg_LHE_output, input_LHE_name):
     """! Reformat DYNNLO events to fit ATLAS conventions."""
     logger.info("Converting output to LHEv3 format")
 
     # Extract intro, header and footer
-    opening_string = LHE.get_opening_string(output_LHE_name)
-    closing_string = LHE.get_closing_string(output_LHE_name)
-    if opening_string.find("<header>") != -1:
-        intro = opening_string[: opening_string.find("<header>")]
-        header = opening_string[opening_string.find("<header>"): opening_string.find("</header>") + 9]
-    else:
-        intro = opening_string[: opening_string.find("<init>")]
-        header = "<header></header>"
-    init = opening_string[opening_string.find("<init>"): opening_string.find("</init>") + 8]
+    intro = "\n".join([elem for elem in [LHE.opening_tag(input_LHE_name), LHE.comment_block(input_LHE_name)] if elem])
+    header = LHE.header_block(input_LHE_events)
+    init = LHE.init_block(input_LHE_events)
+    closing_string = LHE.postamble(input_LHE_name)
     if closing_string.find("<LesHouchesEvents>") != -1:
         footer = closing_string[closing_string.find("</LesHouchesEvents>"):]
     else:
@@ -179,8 +173,6 @@ def reformat_DYNNLO_events(powheg_LHE_output, output_LHE_name):
 
     # Add nominal weight to header
     logger.info("Adding LHEv3 weight: nominal")
-    # self._nominal_group_name = [x[1] for x in configurator.fixed_parameters if x[0] == "lhrwgt_group_name"][0]
-    # header_elem = LHE.add_weight_to_header(header, self._nominal_group_name, "nominal", 0)
     header_elem = LHE.add_weight_to_header(header, "nominal", "nominal", 0)
 
     # Get list of any weights already used in the range 4000-5000
@@ -192,7 +184,7 @@ def reformat_DYNNLO_events(powheg_LHE_output, output_LHE_name):
     weight_number_offset = max(filter(lambda w: 4000 <= w < 5000, [int(w_elem.attrib["id"]) for w_elem in w_elems]) + [4000])
 
     # Add DYNNLO weights to header
-    dynnlo_weight_names = [x[0] for x in LHE.string_to_weight(LHE.get_first_event(output_LHE_name))]
+    dynnlo_weight_names = [x[0] for x in LHE.string_to_weight(LHE.get_first_event(input_LHE_name))]
     for idx, weight_name in enumerate(dynnlo_weight_names, start=1):
         logger.info("Adding LHEv3 weight: {}".format(weight_name))
         weight_name_to_id[weight_name] = weight_number_offset + idx
@@ -201,10 +193,10 @@ def reformat_DYNNLO_events(powheg_LHE_output, output_LHE_name):
     # Convert Powheg input events to LHEv3 output ones
     logger.info("Converting Powheg output into LHEv3 format")
     with open(powheg_LHE_output, "wb") as f_output:
-        f_output.write(intro)
+        f_output.write("{}\n".format(intro))
         ElementTree.ElementTree(header_elem).write(f_output)
-        f_output.write(init)
-        for event in LHE.event_iterator(output_LHE_name):
+        f_output.write("{}\n".format(init))
+        for event in LHE.event_iterator(input_LHE_name):
             f_output.write(LHE.Powheg2LHEv3(event, weight_name_to_id))
         f_output.write(footer)
 

Generators/PowhegControl/python/algorithms/postprocessors/quark_colour_fixer.py

@@ -23,20 +23,20 @@ def quark_colour_fixer(process, powheg_LHE_output):
         logger.info("Checking for presence of colourless quarks!")
 
         # Get opening and closing strings
-        opening_string = LHE.get_opening_string(powheg_LHE_output)
-        closing_string = LHE.get_closing_string(powheg_LHE_output)
+        preamble = LHE.preamble(powheg_LHE_output)
+        postamble = LHE.postamble(powheg_LHE_output)
 
         # Check events for colourless final-state quarks and colour them
         n_events = 0
         logger.info("Checking events for colourless final-state quarks")
         powheg_LHE_recoloured = "{}.recoloured".format(powheg_LHE_output)
         with open(powheg_LHE_recoloured, "wb") as f_output:
-            f_output.write(opening_string)
+            f_output.write("{}\n".format(preamble))
             for input_event in LHE.event_iterator(powheg_LHE_output):
                 is_event_changed, output_event = LHE.ensure_coloured_quarks(input_event)
                 f_output.write(output_event)
                 n_events += [0, 1][is_event_changed]
-            f_output.write(closing_string)
+            f_output.write(postamble)
         logger.info("Re-coloured final-state quarks in {} events!".format(n_events))
 
         # Make a backup of the original events

Generators/PowhegControl/python/algorithms/postprocessors/reweighter.py

@@ -12,7 +12,7 @@ import shutil
 logger = Logging.logging.getLogger("PowhegControl")
 
 ## Tuple to contain arbitrary grouped weights
-WeightTuple = collections.namedtuple("WeightTuple", ["parameter_settings", "keywords", "ID", "name", "combine", "group"])
+WeightTuple = collections.namedtuple("WeightTuple", ["parameter_settings", "keywords", "ID", "name", "combine", "group", "xml_compatible"])
 
 
 @timed("reweighting")
@@ -23,80 +23,105 @@ def reweighter(process, weight_groups, powheg_LHE_output):
     @param weight_groups      OrderedDict containing groups of weights to add.
     @param powheg_LHE_output  Name of LHE file produced by PowhegBox.
     """
-    # List of tuples holding reweighting information
+    ## List of tuples holding reweighting information
     weight_list = []
 
+    ## List of attributes of a weight_group that are not themselves a weight
+    non_weight_attributes = ["parameter_names", "combination_method", "keywords"]
+
+    ## Dictionary of available keywords for faster XML-reweighting: from rwl_setup_param_weights.f
+    xml_kwds = {"lhans1": "lhapdf", "lhans2": "lhapdf", "renscfact": "renscfact", "facscfact": "facscfact"}
+
     # Initial values for scale, PDF and other weight groups
     _idx_scale_start, _idx_PDF_start, _idx_other_start = 1001, 2001, 3001
 
     # Construct ordered list of weights
     for group_name, weight_group in weight_groups.items():
-        keyword_dict = dict((n, k) for n, k in zip(weight_group["parameter_names"], weight_group["keywords"]))
-        method = weight_group["combination_method"]
-        logger.info("Preparing weight group: {:<19} with {} weights".format(group_name, len(weight_group)))
+        logger.info("Preparing weight group: {:<19} with {} weights".format(group_name, len(weight_group) - len(non_weight_attributes)))
+
+        # Name -> keyword dictionary is different if this is an XML reweighting
+        is_xml_compatible = all([k in xml_kwds.keys() for _kw_set in weight_group["keywords"] for k in _kw_set])
+        if is_xml_compatible:
+            _keywords = [list(set([xml_kwds[k] for k in _kw_set])) for _kw_set in weight_group["keywords"]]
+        else:
+            logger.warning("... this weight group is incompatible with XML-style reweighting. Will use (slow) old method.")
+            _keywords = weight_group["keywords"]
+        keyword_dict = dict((n, k) for n, k in zip(weight_group["parameter_names"], _keywords))
+
+        # Common arguments for all weights
+        tuple_kwargs = {"keywords": keyword_dict, "combine": weight_group["combination_method"], "group": group_name, "xml_compatible": is_xml_compatible}
 
         # Nominal variation: ID=0
         if group_name == "nominal":
-            weight_list.append(WeightTuple(parameter_settings=weight_group["nominal"], keywords=keyword_dict, ID=0, name="nominal", combine=method, group=group_name))
+            weight_list.append(WeightTuple(parameter_settings=weight_group["nominal"], ID=0, name="nominal", **tuple_kwargs))
 
         # Scale variations: ID=1001+
         elif group_name == "scale_variation":
-            for idx, name in enumerate([k for k in weight_group.keys() if k not in ["parameter_names", "combination_method", "keywords"]], start=_idx_scale_start):
-                weight_list.append(WeightTuple(parameter_settings=weight_group[name], keywords=keyword_dict, ID=idx, name=name, combine=method, group=group_name))
+            for idx, name in enumerate([k for k in weight_group.keys() if k not in non_weight_attributes], start=_idx_scale_start):
+                weight_list.append(WeightTuple(parameter_settings=weight_group[name], ID=idx, name=name, **tuple_kwargs))
             _idx_scale_start = idx + 1
 
         # PDF variations: ID=2001+
         elif group_name == "PDF_variation":
-            for idx, name in enumerate([k for k in weight_group.keys() if k not in ["parameter_names", "combination_method", "keywords"]], start=_idx_PDF_start):
-                weight_list.append(WeightTuple(parameter_settings=weight_group[name], keywords=keyword_dict, ID=idx, name=name, combine=method, group=group_name))
+            for idx, name in enumerate([k for k in weight_group.keys() if k not in non_weight_attributes], start=_idx_PDF_start):
+                weight_list.append(WeightTuple(parameter_settings=weight_group[name], ID=idx, name=name, **tuple_kwargs))
             _idx_PDF_start = idx + 1
 
         # Other variations: ID=3001+
         else:
-            for idx, name in enumerate([k for k in weight_group.keys() if k not in ["parameter_names", "combination_method", "keywords"]], start=_idx_other_start):
-                weight_list.append(WeightTuple(parameter_settings=weight_group[name], keywords=keyword_dict, ID=idx, name=name, combine=method, group=group_name))
+            for idx, name in enumerate([k for k in weight_group.keys() if k not in non_weight_attributes], start=_idx_other_start):
+                weight_list.append(WeightTuple(parameter_settings=weight_group[name], ID=idx, name=name, **tuple_kwargs))
             _idx_other_start = idx + 1
 
     # Construct xml output
-    xml_lines, weightgroup = [], ""
+    xml_lines, non_xml_weight_list, current_weightgroup = [], [], ""
     for weight in weight_list:
         # Group elements
-        if weight.group != weightgroup:
+        if weight.group != current_weightgroup:
             xml_lines.append("</weightgroup>")
-            weightgroup = weight.group
-            xml_lines.append("<weightgroup name='{}'>".format(weightgroup))
+            current_weightgroup = weight.group
+            xml_lines.append("<weightgroup name='{}'>".format(current_weightgroup))
         # Weight elements
-        keyword_pairs = ["{}={}".format(k, v) for p, v in weight.parameter_settings for k in weight.keywords[p]]
-        xml_lines.append("<weight id='{}'> {} </weight>".format(weight.ID, " ".join(keyword_pairs)))
-    xml_lines.append(xml_lines.pop(0))
+        if weight.xml_compatible:
+            keyword_pairs = ["{}={}".format(k, v) for p, v in weight.parameter_settings for k in weight.keywords[p]]
+            xml_lines.append("<weight id='{}'> {} </weight>".format(weight.ID, " ".join(keyword_pairs)))
+            print "group {}, <weight id='{}'> {} </weight>".format(current_weightgroup, weight.ID, " ".join(keyword_pairs))
+        else:
+            non_xml_weight_list.append(weight)
+    xml_lines.append(xml_lines.pop(0)) # move the closing tag to the end
 
     # Write xml output
-    with open("reweighting_input.xml", "wb") as f_rwgt:
-        f_rwgt.write("<initrwgt>\n")
-        for xml_line in xml_lines:
-            f_rwgt.write("{}\n".format(xml_line))
-        f_rwgt.write("</initrwgt>")
-
-    # Make backup of generation statistics
-    if os.path.isfile("pwgcounters.dat"):
-        shutil.copy("pwgcounters.dat", "pwgcounters.dat.bak")
-
-    # # Iterate over variations
-    # for idx_weight, weight in enumerate(weight_list, start=1):
-    #     add_single_weight(process, weight, idx_weight, len(weight_list))
-
-    # Add reweighting lines to runcard
-    FileParser("powheg.input").text_replace("rwl_file .*", "rwl_file 'reweighting_input.xml'")
-    FileParser("powheg.input").text_replace("rwl_add .*", "rwl_add 1")
-    FileParser("powheg.input").text_replace("clobberlhe .*", "clobberlhe 1")
-
-    # Copy the old events and then run the reweighter
-    shutil.copy(powheg_LHE_output, "{}.unweighted".format(powheg_LHE_output))
-    logger.info("Preparing to add {} weights to generated events.".format(len(weight_list)))
-    singlecore_untimed(process)
-
-    # Move the reweighted file back
-    shutil.move(powheg_LHE_output.replace(".lhe", "-rwgt.lhe"), powheg_LHE_output)
+    n_xml_weights = len(weight_list) - len(non_xml_weight_list)
+    if n_xml_weights > 0:
+        with open("reweighting_input.xml", "wb") as f_rwgt:
+            f_rwgt.write("<initrwgt>\n")
+            for xml_line in xml_lines:
+                f_rwgt.write("{}\n".format(xml_line))
+            f_rwgt.write("</initrwgt>")
+
+        # Make backup of generation statistics
+        if os.path.isfile("pwgcounters.dat"):
+            shutil.copy("pwgcounters.dat", "pwgcounters.dat.bak")
+
+        # Add reweighting lines to runcard
+        FileParser("powheg.input").text_replace("rwl_file .*", "rwl_file 'reweighting_input.xml'")
+        FileParser("powheg.input").text_replace("rwl_add .*", "rwl_add 1")
+        FileParser("powheg.input").text_replace("clobberlhe .*", "clobberlhe 1")
+
+        # Copy the old events and then run the reweighter
+        shutil.copy(powheg_LHE_output, "{}.unweighted".format(powheg_LHE_output))
+        logger.info("Preparing simultaneous calculation of {} additional weights for generated events.".format(n_xml_weights))
+        singlecore_untimed(process)
+
+        # Move the reweighted file back
+        shutil.move(powheg_LHE_output.replace(".lhe", "-rwgt.lhe"), powheg_LHE_output)
+
+    # Iterate over any variations which require old-style reweighting
+    if len(non_xml_weight_list) > 0:
+        logger.info("Preparing individual calculation of {} additional weights for generated events.".format(len(non_xml_weight_list)))
+        for idx_weight, weight in enumerate(non_xml_weight_list, start=1):
+            add_single_weight(process, weight, idx_weight, len(non_xml_weight_list))
+            shutil.move(powheg_LHE_output.replace(".lhe", "-rwgt.lhe"), powheg_LHE_output)
 
     # Remove rwgt and pdf lines, which crash Pythia
     FileParser(powheg_LHE_output).text_remove("^#pdf")
@@ -116,45 +141,44 @@ def reweighter(process, weight_groups, powheg_LHE_output):
         shutil.move("pwgcounters.dat.bak", "pwgcounters.dat")
 
 
-# def add_single_weight(process, weight, idx_weight, n_total):
-#     """! Add a single additional weight to an LHE file.
-#
-#     @param process       PowhegBox process.
-#     @param weight     Tuple with details of the weight to be added.
-#     @param idx_weight Which number this weight is.
-#     @param n_total    Total number of weights to add.
-#     """
-#     @timed("weight variation {}/{}".format(idx_weight, n_total))
-#     def __timed_inner_fn():
-#         logger.info("... weight name is:     {}".format(weight.name))
-#         logger.info("... weight index ID is: {}".format(weight.ID))
-#         shutil.copy("powheg_nominal.input", "powheg.input")
-#
-#         # As the nominal process has already been run, turn on compute_rwgt
-#         FileParser("powheg.input").text_replace("compute_rwgt 0", "compute_rwgt 1")
-#
-#         # Ensure that manyseeds is turned off,  as this would cause the reweighting to crash
-#         FileParser("powheg.input").text_replace("manyseeds .*", "manyseeds 0")
-#         FileParser("powheg.input").text_replace("parallelstage .*", "parallelstage -1")
-#
-#         # Apply weight settings
-#         for (parameter, value) in weight.parameter_settings:
-#             try:
-#                 for keyword in weight.keywords[parameter]:
-#                     FileParser("powheg.input").text_replace("^{}.*".format(keyword), "{} {}".format(keyword, value))
-#                 logger.info("... setting {} to {}".format(parameter, value))
-#             except KeyError:
-#                 logger.warning("Parameter '{}' not recognised. Cannot reweight!".format(parameter))
-#
-#         # Set reweighting metadata
-#         FileParser("powheg.input").text_replace("lhrwgt_descr .*", "lhrwgt_descr '{}'".format(weight.name))
-#         FileParser("powheg.input").text_replace("lhrwgt_id .*", "lhrwgt_id '{}'".format(weight.ID))
-#         FileParser("powheg.input").text_replace("lhrwgt_group_combine .*", "lhrwgt_group_combine '{}'".format(weight.combine))
-#         FileParser("powheg.input").text_replace("lhrwgt_group_name .*", "lhrwgt_group_name '{}'".format(weight.group))
-#
-#         # Run the process until termination
-#         untimed_generation_singlecore(process)
-#         shutil.move("pwgevents-rwgt.lhe", "pwgevents.lhe")
-#
-#     # Run single weight variation
-#     __timed_inner_fn()
+def add_single_weight(process, weight, idx_weight, n_total):
+    """! Add a single additional weight to an LHE file.
+
+    @param process    PowhegBox process.
+    @param weight     Tuple with details of the weight to be added.
+    @param idx_weight Which number this weight is.
+    @param n_total    Total number of weights to add.
+    """
+    @timed("weight variation {}/{}".format(idx_weight, n_total))
+    def __timed_inner_fn():
+        logger.info("... weight name is:     {}".format(weight.name))
+        logger.info("... weight index ID is: {}".format(weight.ID))
+        shutil.copy("powheg_nominal.input", "powheg.input")
+
+        # As the nominal process has already been run, turn on compute_rwgt
+        FileParser("powheg.input").text_replace("compute_rwgt 0", "compute_rwgt 1")
+
+        # Ensure that manyseeds is turned off,  as this would cause the reweighting to crash
+        FileParser("powheg.input").text_replace("manyseeds .*", "manyseeds 0")
+        FileParser("powheg.input").text_replace("parallelstage .*", "parallelstage -1")
+
+        # Apply weight settings
+        for (parameter, value) in weight.parameter_settings:
+            try:
+                for keyword in weight.keywords[parameter]:
+                    FileParser("powheg.input").text_replace("^{} .*".format(keyword), "{} {}".format(keyword, value))
+                logger.info("... setting {} to {}".format(parameter, value))
+            except KeyError:
+                logger.warning("Parameter '{}' not recognised. Cannot reweight!".format(parameter))
+
+        # Set reweighting metadata
+        FileParser("powheg.input").text_replace("lhrwgt_descr .*", "lhrwgt_descr '{}'".format(weight.name))
+        FileParser("powheg.input").text_replace("lhrwgt_id .*", "lhrwgt_id '{}'".format(weight.ID))
+        FileParser("powheg.input").text_replace("lhrwgt_group_combine .*", "lhrwgt_group_combine '{}'".format(weight.combine))
+        FileParser("powheg.input").text_replace("lhrwgt_group_name .*", "lhrwgt_group_name '{}'".format(weight.group))
+
+        # Run the process until termination
+        singlecore_untimed(process)
+
+    # Run single weight variation
+    __timed_inner_fn()

Generators/PowhegControl/python/parameters/registry.py

@@ -207,11 +207,13 @@ class Registry(object):
         self.add_default("lhrwgt_id", 0, frozen=True, description="weight ID.")
         self.add_default("LOevents", 0, description="(0:disabled; 1:enabled) produce LOPS events (scalup=ptj); in this case bornonly should also be enabled.")
         self.add_default("m2bb", -1, description="(-1:use Powheg default)")
-        self.add_default("MadSpin_decays", [], description="Decays allowed by MadSpin")
-        self.add_default("MadSpin_enabled", True, description="(False:disabled; True:enabled) use MadSpin for top decays (only if decays are disabled in Powheg)")
+        self.add_default("MadSpin_decays", [], description="decays allowed by MadSpin")
+        self.add_default("MadSpin_enabled", True, description="use MadSpin for top decays (only if decays are disabled in Powheg). [False:disabled; True:enabled]")
         self.add_default("MadSpin_model", "loop_sm-ckm", description="which model to import in MadSpin.")
-        self.add_default("MadSpin_mode", "full", description="('full'; 'bridge'; 'none') which spin mode to use in MadSpin")
-        self.add_default("MadSpin_process", "", description="Process that MadSpin is operating on")
+        self.add_default("MadSpin_mode", "full", description="which spin mode to use in MadSpin. ['full'; 'bridge'; 'none']")
+        self.add_default("MadSpin_nFlavours", "4", description="which flavour scheme to use. [4; 5]")
+        self.add_default("MadSpin_process", "", description="process that MadSpin is operating on")
+        self.add_default("MadSpin_taus_are_leptons", True, description="whether lepton definitions should include taus. [False: do not include taus, True: include taus]")
         self.add_default("manyseeds", 0, description="read multiple seeds for the random number generator from pwgseeds.dat. [1:enabled]")
         self.add_default("mass_b", atlas_common.mass.b, description="b-quark mass in GeV")
         self.add_default("mass_c", atlas_common.mass.c, name="mass_c", description="c-quark mass in GeV")

Generators/PowhegControl/python/processes/external/external_madspin.py

@@ -36,7 +36,9 @@ class ExternalMadSpin(ExternalBase):
         self.add_keyword("MadSpin_enabled")
         self.add_keyword("MadSpin_model")
         self.add_keyword("MadSpin_mode")
+        self.add_keyword("MadSpin_nFlavours")
         self.add_keyword("MadSpin_process", process)
+        self.add_keyword("MadSpin_taus_are_leptons")
         self.add_keyword("mass_b")
         self.add_keyword("mass_H")
         self.add_keyword("mass_t")
@@ -70,7 +72,7 @@ class ExternalMadSpin(ExternalBase):
             return False
         # Check that decay list has at least one entry otherwise add defaults
         if len(self.MadSpin_decays) == 0:
-            logger.info("No MadSpin decays specified, so defaults will be used (only relevant if MadSpin is enabled)")
+            logger.warning("No MadSpin decays specified, so defaults will be used.")
             for decay in ("t > w+ b, w+ > l+ vl", "t~ > w- b~, w- > l- vl~", "t > w+ b, w+ > j j", "t~ > w- b~, w- > j j"):
                 logger.info("... adding MadSpin decay: 'decay {0}'".format(decay))
                 self.parameters_by_name("MadSpin_decays")[0].value.append("decay {0}".format(decay))

Generators/PowhegControl/python/processes/powheg/Wj.py

@@ -36,6 +36,7 @@ class Wj(PowhegV2):
 
         # Add all keywords for this process, overriding defaults if required
         self.add_keyword("alphaem")
+        self.add_keyword("alphas_from_lhapdf")
         self.add_keyword("bornktmin", 5.0)
         self.add_keyword("bornonly")
         self.add_keyword("bornsuppfact")

Generators/PowhegControl/python/utility/LHE.py

@@ -27,9 +27,9 @@ def merge(input_file_pattern, output_file):
     with open(output_file, "ab") as f_output:
         logger.info("... reading metadata from {}".format(input_file_list[0]))
         # Start with the first file and extract opening/closing string
-        opening_string = get_opening_string(input_file_list[0])
-        closing_string = get_closing_string(input_file_list[0])
-        f_output.write(opening_string)
+        s_preamble = preamble(input_file_list[0])
+        s_postamble = postamble(input_file_list[0])
+        f_output.write(s_preamble)
 
         # Now write events from all files to output
         # Use sanitised list in case output_file matches the pattern
@@ -37,7 +37,7 @@ def merge(input_file_pattern, output_file):
             f_output.write(event)
 
         # Finally add the footer
-        f_output.write(closing_string)
+        f_output.write(s_postamble)
 
     # Write the event count to the logger
     logger.info("Wrote {} events to {}".format(nEvents + 1, output_file))
@@ -57,9 +57,10 @@ def event_iterator(input_files, verbose=True):
         # Group all lines inside an XML event element
         with open(file_name, "rb") as f_input:
             for line in f_input:
-                if "<event>" in line:
+                # Both <event ...> and <event> are permitted
+                if "<event" in line:
                     in_event = True
-                    line = line[line.index("<event>"):]  # catch cases like "</init><event>"
+                    line = line[line.index("<event"):]  # catch cases like "</init><event>"
                 if in_event:
                     event_lines += line
                 if "</event>" in line:
@@ -105,7 +106,7 @@ def add_weight_to_header(header, weightgroup_name, weight_name, weight_id):
     return header_elem
 
 
-def get_opening_string(input_LHE_file):
+def preamble(input_LHE_file):
     """! Get opening lines from file as a string."""
     with open(input_LHE_file, "rb") as f_input:
         s_input = mmap.mmap(f_input.fileno(), 0, access=mmap.ACCESS_READ)
@@ -113,19 +114,55 @@ def get_opening_string(input_LHE_file):
     return s_output
 
 
-def get_first_event(input_LHE_file):
-    """! Get first event from file as a string."""
+def postamble(input_LHE_file):
+    """! Get closing lines from file as a string."""
     with open(input_LHE_file, "rb") as f_input:
         s_input = mmap.mmap(f_input.fileno(), 0, access=mmap.ACCESS_READ)
-        s_output = s_input[s_input.find("<event>"): s_input.find("</event>") + 8]
+        s_output = s_input[s_input.rfind("</event>") + 9:]
     return s_output
 
 
-def get_closing_string(input_LHE_file):
-    """! Get closing lines from file as a string."""
+def opening_tag(input_LHE_file):
+    """! Get <LesHouchesEvents> opening tag from file as a string."""
+    s_opening = preamble(input_LHE_file)
+    if s_opening.find("<LesHouchesEvents") != -1:
+        return s_opening[s_opening.find("<LesHouchesEvents"): s_opening.find(">") + 1].strip("\n")
+    return ""
+
+
+def comment_block(input_LHE_file):
+    """! Get comment block from file as a string."""
+    s_opening = preamble(input_LHE_file)
+    if s_opening.find("<!--") != -1:
+        return s_opening[s_opening.find("<!--"): s_opening.find("-->") + 3].strip("\n")
+    return ""
+
+
+def header_block(input_LHE_file):
+    """! Get <header> block from file as a string."""
+    s_opening = preamble(input_LHE_file)
+    if s_opening.find("<header>") != -1:
+        return s_opening[s_opening.find("<header>"): s_opening.find("</header>") + 9].strip("\n")
+    return "<header>\n</header>"
+
+
+def init_block(input_LHE_file):
+    """! Get <init> block from file as a string."""
+    s_opening = preamble(input_LHE_file)
+    # Both <init ...> and <init> are permitted
+    if s_opening.find("<init>") != -1:
+        return s_opening[s_opening.find("<init>"): s_opening.find("</init>") + 7].strip("\n")
+    if s_opening.find("<init ") != -1:
+        return s_opening[s_opening.find("<init "): s_opening.find("</init>") + 7].strip("\n")
+    return "<init>\n</init>"
+
+
+def get_first_event(input_LHE_file):
+    """! Get first event from file as a string."""
     with open(input_LHE_file, "rb") as f_input:
         s_input = mmap.mmap(f_input.fileno(), 0, access=mmap.ACCESS_READ)
-        s_output = s_input[s_input.rfind("</event>") + 9:]
+        # Both <event ...> and <event> are permitted
+        s_output = s_input[s_input.find("<event"): s_input.find("</event>") + 8]
     return s_output