Fixed missing bracket and added workaround for colourless quarks from hadronic...

Fixed missing bracket and added workaround for colourless quarks from hadronic Z decays. 2016-12-01 13:57:05. (PowhegControl-00-02-21)


Former-commit-id: 9d71f300f160b6a8f75b51b818cb03451275d9d4

Generators/PowhegControl/python/PowhegConfig_base.py

@@ -183,7 +183,7 @@ class PowhegConfig_base(object):
                 f.write("{:<30}! [ATLAS default: {}] {}\n".format("{} {}".format(name, value), default, desc))
                 # Print warnings for specific parameters
                 if name == "bornsuppfact" and value > 0:
-                    self.logger.warning("Born-level suppression is enabled:"
+                    self.logger.warning("Born-level suppression is enabled:")
                     self.logger.warning("-> the cross-section passed to the parton shower will be inaccurate.")
                     self.logger.warning("-> please use the cross-section printed in the log file before showering begins.")
                     self.logger.warning("-> in addition, using this in conjunction with J-slicing may give problems.")
@@ -291,6 +291,11 @@ class PowhegConfig_base(object):
         if self.bornsuppfact > 0.0:
             strategies.afterburner_mean_event_weight_calculator(self.powheg_LHE_output)
 
+        # Run quark colour fix for gg_HZ
+        if hasattr(self, "_require_quark_colour_check") and hasattr(self, "vdecaymode"):
+            if (self.vdecaymode == 0 or self.vdecaymode == 10): # hadronic decay modes
+                strategies.afterburner_ensure_coloured_quarks(self.powheg_LHE_output)
+
         # Run MadSpin afterburner if requested
         if hasattr(self, "_MadSpin_executable"):
             strategies.afterburner_MadSpin(self)

Generators/PowhegControl/python/processes/PowhegConfig_ggF_HZ.py

@@ -21,6 +21,7 @@ class PowhegConfig_ggF_HZ(PowhegConfig_base):
     def __init__(self, runArgs=None, opts=None, LO_process=True):
         super(PowhegConfig_ggF_HZ, self).__init__(runArgs, opts)
         self._powheg_executable += "/ggHZ/pwhg_main"
+        self._require_quark_colour_check = True
 
         # Decorate with generic option sets
         self.add_parameter_set("Higgs + vector boson", boson="Z")
@@ -39,5 +40,6 @@ class PowhegConfig_ggF_HZ(PowhegConfig_base):
         self.nubound  = 60000
         self.xupbound = 6
         self.itmx1    = 1
+        # vdecaymode => 1: e, 2: mu, 3: tau, 4:nu_e, 5: nu_mu, 6: nu_tau, 0: to hadrons, 10: inclusive (all hadrons and leptons), 11: inclusive leptons, 12: inclusive neutrinos
 
         self.populate_default_strings()

Generators/PowhegControl/python/strategies/AfterburnerEnsureColouredQuarks.py

+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+# @PowhegControl AfterburnerEnsureColouredQuarks
+#  Ensure that all final-state quarks have a defined colour
+#
+#  Authors: James Robinson  <james.robinson@cern.ch>
+
+#! /usr/bin/env python
+from AthenaCommon.Logging import logging
+from ..utility import LHEUtils
+import shutil
+
+# Initialise logging handler
+logger = logging.getLogger("PowhegControl")
+
+# Post-process existing events to ensure that quarks are coloured
+def afterburner_ensure_coloured_quarks(LHE_file_name):
+    logger.info("Checking for presence of colourless quarks!")
+
+    # Move file to backup location
+    LHE_file_name_backup = "{}.colourless_quarks_backup".format(LHE_file_name)
+    shutil.move(LHE_file_name, LHE_file_name_backup)
+
+    # Get opening and closing strings
+    opening_string = LHEUtils.get_opening_string(LHE_file_name_backup)
+    closing_string = LHEUtils.get_closing_string(LHE_file_name_backup)
+
+    # Checking events for colourless final-state quarks
+    n_events = 0
+    logger.info("Checking events for colourless final-state quarks")
+    with open(LHE_file_name, "wb") as f_output:
+        f_output.write(opening_string)
+        for input_event in LHEUtils.event_iterator(LHE_file_name_backup):
+            is_event_changed, output_event = LHEUtils.ensure_coloured_quarks(input_event)
+            f_output.write(output_event)
+            if is_event_changed : n_events += 1
+        f_output.write(closing_string)
+    logger.info("Re-coloured final-state quarks in {0} events!".format(n_events))

Generators/PowhegControl/python/strategies/AfterburnerMeanEventWeightCalculator.py

@@ -2,7 +2,6 @@
 
 # @PowhegControl AfterburnerMeanEventWeightCalculator
 #  Calculate the mean event weight of pre-generated Powheg events
-#  Apply PHOTOS afterburner to pre-generated Powheg events
 #
 #  Authors: James Robinson  <james.robinson@cern.ch>
 
@@ -13,7 +12,7 @@ from ..utility import LHEUtils
 # Initialise logging handler
 logger = logging.getLogger("PowhegControl")
 
-# Run PHOTOS as an afterburner to existing events
+# Read existing events to calculate mean event weights
 def afterburner_mean_event_weight_calculator(LHE_file_name):
     logger.info("Born-level suppression is enabled so the cross-section MUST be recalculated!")
     # LHEUtils returns a generator, since the input file may be large

Generators/PowhegControl/python/strategies/__init__.py

 # Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 
+from AfterburnerEnsureColouredQuarks import *
 from AfterburnerMadSpin import *
 from AfterburnerMeanEventWeightCalculator import *
 from AfterburnerNNLOReweighting import *

Generators/PowhegControl/python/utility/LHEUtils.py

@@ -6,6 +6,7 @@
 #  Authors: James Robinson  <james.robinson@cern.ch>
 
 #! /usr/bin/env python
+import re
 import glob
 import mmap
 from AthenaCommon.Logging import logging
@@ -136,6 +137,30 @@ def string_to_weight(input_event):
     return [(line.split(" ")[-1], line.split(" ")[2]) for line in weight_lines]
 
 
+# Ensure that all final-state quarks in the event are coloured
+def ensure_coloured_quarks(input_event):
+    initial_colour_flow, is_event_changed = -1, False
+    event_lines = ""
+    for input_line in input_event.splitlines(True):
+        output_line = None
+        try: # interpret line as a particle
+            IDUP, ISTUP, MOTHUP0, MOTHUP1, ICOLUP0, ICOLUP1, PUP0, PUP1, PUP2, PUPU3, PUP4, VTIMUP, SPINUP = input_line.split()
+            if int(IDUP) == 21 and int(ISTUP) == -1: # this is an initial state gluon
+                initial_colour_flow = max(initial_colour_flow, int(ICOLUP0), int(ICOLUP1))
+            if int(ICOLUP0) == 0 and int(ICOLUP0) == 0: # this is a colourless particle
+                tokens = re.split(r"(\s+)", input_line)
+                if int(IDUP) in range(1, 7): # this is a quark
+                    output_line = "".join(tokens[:9])+" {0:>5d} {1:>5d}".format(initial_colour_flow+1, 0)+"".join(tokens[13:])
+                    is_event_changed = True
+                if int(IDUP) in range(-6, 0): # this is an anti-quark
+                    output_line = "".join(tokens[:9])+" {0:>5d} {1:>5d}".format(0, initial_colour_flow+1)+"".join(tokens[13:])
+                    is_event_changed = True
+        except ValueError: # this is not a particle line
+            pass
+        event_lines += output_line if output_line is not None else input_line
+    return (is_event_changed, event_lines)
+
+
 # Get new-style event weights from an input event string
 def Powheg2LHEv3(input_event, name_to_ID):
     # Return event as-is if there are no Powheg-style weights