Merge branch 'gfrattar-powheg-gg4l-swapdecays' into 'master'

Implementation of non-currently-supported decays in gg4l process via change of pdg ID in LHE files AGENE-2156

See merge request !61671

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

@@ -8,6 +8,9 @@ from .integration_gridpack_creator import integration_gridpack_creator
 from .integration_grid_tester import integration_grid_tester
 from .madspin import MadSpin
 from .mu2tau import mu2tau
+from .mu2e import mu2e
+from .e2mu import e2mu
+from .e2tau import e2tau
 from .nnlo_reweighter import NNLO_reweighter
 from .output_file_renamer import output_file_renamer
 from .output_tarball_preparer import output_tarball_preparer

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

+# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
+
+from AthenaCommon import Logging
+from ...decorators import timed
+from ...utility import LHE
+import shutil
+
+## Get handle to Athena logging
+logger = Logging.logging.getLogger("PowhegControl")
+
+
+@timed("e2mu")
+def e2mu(powheg_LHE_output):
+    """! Post-process existing events from electrons to muons
+
+    @param powheg_LHE_output  Name of LHE file produced by PowhegBox.
+
+    @author Jan Kretzschmar <jan.kretzschmar@cern.ch>
+    """
+    logger.warning("Converting LHE events from electron to muon decays.")
+
+    # Get opening and closing strings
+    preamble = LHE.preamble(powheg_LHE_output)
+    postamble = LHE.postamble(powheg_LHE_output)
+    
+    n_events = 0
+    powheg_LHE_mu = "{}.mu".format(powheg_LHE_output)
+    with open(powheg_LHE_mu, "wb") as f_output:
+        f_output.write("{}\n".format(preamble))
+        for input_event in LHE.event_iterator(powheg_LHE_output):
+            is_event_changed, output_event = LHE.e2mu(input_event)
+            f_output.write(output_event)
+            n_events += [0, 1][is_event_changed]
+        f_output.write(postamble)
+    logger.info("Changed e->mu in {} events!".format(n_events))
+
+    # Make a backup of the original events
+    shutil.move(powheg_LHE_output, "{}.e2mu_backup".format(powheg_LHE_output))
+    shutil.move(powheg_LHE_mu, powheg_LHE_output)
\ No newline at end of file

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

+# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
+
+from AthenaCommon import Logging
+from ...decorators import timed
+from ...utility import LHE
+import shutil
+
+## Get handle to Athena logging
+logger = Logging.logging.getLogger("PowhegControl")
+
+
+@timed("e2tau")
+def e2tau(powheg_LHE_output):
+    """! Post-process existing events from electrons to taus
+
+    @param powheg_LHE_output  Name of LHE file produced by PowhegBox.
+
+    @author Jan Kretzschmar <jan.kretzschmar@cern.ch>
+    """
+    logger.warning("Converting LHE events from electron to tau decays. Tau mass must be restored by showering program, ensure to validate physics.")
+
+    # Get opening and closing strings
+    preamble = LHE.preamble(powheg_LHE_output)
+    postamble = LHE.postamble(powheg_LHE_output)
+    
+    n_events = 0
+    powheg_LHE_tau = "{}.tau".format(powheg_LHE_output)
+    with open(powheg_LHE_tau, "wb") as f_output:
+        f_output.write("{}\n".format(preamble))
+        for input_event in LHE.event_iterator(powheg_LHE_output):
+            is_event_changed, output_event = LHE.e2tau(input_event)
+            f_output.write(output_event)
+            n_events += [0, 1][is_event_changed]
+        f_output.write(postamble)
+    logger.info("Changed mu->tau in {} events!".format(n_events))
+
+    # Make a backup of the original events
+    shutil.move(powheg_LHE_output, "{}.e2tau_backup".format(powheg_LHE_output))
+    shutil.move(powheg_LHE_tau, powheg_LHE_output)

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

+# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
+
+from AthenaCommon import Logging
+from ...decorators import timed
+from ...utility import LHE
+import shutil
+
+## Get handle to Athena logging
+logger = Logging.logging.getLogger("PowhegControl")
+
+
+@timed("mu2e")
+def mu2e(powheg_LHE_output):
+    """! Post-process existing events from muons to electrons
+
+    @param powheg_LHE_output  Name of LHE file produced by PowhegBox.
+
+    @author Jan Kretzschmar <jan.kretzschmar@cern.ch>
+    """
+    logger.warning("Converting LHE events from muon to electron decays.")
+
+    # Get opening and closing strings
+    preamble = LHE.preamble(powheg_LHE_output)
+    postamble = LHE.postamble(powheg_LHE_output)
+    
+    n_events = 0
+    powheg_LHE_e = "{}.e".format(powheg_LHE_output)
+    with open(powheg_LHE_e, "wb") as f_output:
+        f_output.write("{}\n".format(preamble))
+        for input_event in LHE.event_iterator(powheg_LHE_output):
+            is_event_changed, output_event = LHE.mu2e(input_event)
+            f_output.write(output_event)
+            n_events += [0, 1][is_event_changed]
+        f_output.write(postamble)
+    logger.info("Changed mu->e in {} events!".format(n_events))
+
+    # Make a backup of the original events
+    shutil.move(powheg_LHE_output, "{}.mu2e_backup".format(powheg_LHE_output))
+    shutil.move(powheg_LHE_e, powheg_LHE_output)
\ No newline at end of file

Generators/PowhegControl/python/algorithms/scheduler.py

@@ -33,6 +33,9 @@ class Scheduler(object):
         "NNLO reweighter",
         "LHE file nominal weight updater",
         "mu2tau",
+        "mu2e",
+        "e2tau",
+        "e2mu",
         "MadSpin",
         "integration grid tester",
         "cross section calculator",
@@ -71,6 +74,9 @@ class Scheduler(object):
         "reweighter": partial(postprocessors.reweighter, powheg_LHE_output=powheg_LHE_output),
         "quark colour fixer": partial(postprocessors.quark_colour_fixer, powheg_LHE_output=powheg_LHE_output),
         "mu2tau": partial(postprocessors.mu2tau, powheg_LHE_output=powheg_LHE_output),
+        "mu2e": partial(postprocessors.mu2e, powheg_LHE_output=powheg_LHE_output),
+        "e2tau": partial(postprocessors.e2tau, powheg_LHE_output=powheg_LHE_output),
+        "e2mu": partial(postprocessors.e2mu, powheg_LHE_output=powheg_LHE_output),
         "LHE file cleaner": partial(postprocessors.lhe_cleaner, powheg_LHE_output=powheg_LHE_output),
         "LHE file nominal weight updater": partial(postprocessors.lhe_nominal_weight_updater, powheg_LHE_output=powheg_LHE_output),
     }

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

@@ -221,6 +221,31 @@ class gg4l(PowhegRES):
         if "\'" not in self.contr:
             self.contr = "\'"+self.contr+"\'"
 
+        #Modify unsupported decay configuration in ZZ production
+        if (self.proc == "'ZZ'" and ((self.vdecaymodeV1 == self.vdecaymodeV2) or self.vdecaymodeV1 == 15 or self.vdecaymodeV2 == 15)): 
+            logger.warning("Powheg/gg4l does support directly 4e, 4mu or tau final states.")
+            if(self.vdecaymodeV1 == 11 and self.vdecaymodeV2 == 11):
+                logger.warning("Ask to generate 2e2mu decays and hack the LHE files to have 4e final states - make sure to validate!")
+                self.add_algorithm("mu2e")
+            elif(self.vdecaymodeV1 == 13 and self.vdecaymodeV2 == 13):
+                logger.warning("Ask to generate 2e2mu decays and hack the LHE files to have 4mu final states - make sure to validate!")
+                self.add_algorithm("e2mu")
+            elif(self.vdecaymodeV1 == 15 and self.vdecaymodeV2 == 15):
+                logger.warning("Ask to generate 2e2mu decays and hack the LHE files to have 4tau final states - make sure to validate!")
+                self.add_algorithm("e2tau")
+                self.add_algorithm("mu2tau")
+            elif(self.vdecaymodeV1 == 11 and self.vdecaymodeV2 == 15) or (self.vdecaymodeV1 == 15 and self.vdecaymodeV2 == 11):
+                logger.warning("Ask to generate 2e2mu decays and hack the LHE files to have 2e2tau final states - make sure to validate!")
+                self.add_algorithm("mu2tau")
+            elif(self.vdecaymodeV1 == 13 and self.vdecaymodeV2 == 15) or (self.vdecaymodeV1 == 15 and self.vdecaymodeV2 == 13):
+                logger.warning("Ask to generate 2e2mu decays and hack the LHE files to have 2mu2tau final states - make sure to validate!")
+                self.add_algorithm("e2tau")
+            
+            self.vdecaymodeV1 = 11
+            self.vdecaymodeV2 = 13
+            self.parameters_by_keyword("vdecaymodeV1")[0].value = self.vdecaymodeV1
+            self.parameters_by_keyword("vdecaymodeV2")[0].value = self.vdecaymodeV2
+
         #check if the setting is allowed
         if self.proc not in self.allowed_process_modes:
             logger.warning("Process mode {} not recognised!".format(self.proc))

Generators/PowhegControl/python/utility/LHE.py

@@ -222,6 +222,81 @@ def mu2tau(input_event):
         event_lines += output_line if output_line is not None else input_line
     return (is_event_changed, event_lines)
 
+def e2tau(input_event):
+    """! 
+    Swap out electrons for taus, and electron neutrinos for tau neutrinos.
+    Note no momentum reshuffling is done, but Pythia appears to restore the correct tau mass.
+    """
+    is_event_changed = False
+    event_lines = ""
+    for input_line in input_event.splitlines(True):
+        output_line = None
+        try:  # interpret line as a particle
+            tokens = re.split(r"(\s+)", input_line)
+            if len(tokens) < 25: raise ValueError
+            IDUP = int(tokens[2])
+            if abs(IDUP) == 11 or abs(IDUP) == 12:  # this is a electron or electron neutrino
+                if IDUP > 0:
+                    IDUP += 4
+                else:
+                    IDUP -= 4
+                is_event_changed = True
+                output_line = "".join("".join(tokens[:2])+str(IDUP)+"".join(tokens[3:]))
+        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)
+
+def mu2e(input_event):
+    """! 
+    Swap out muons for electrons, and muon neutrinos for electron neutrinos.
+    Note no momentum reshuffling is done.
+    """
+    is_event_changed = False
+    event_lines = ""
+    for input_line in input_event.splitlines(True):
+        output_line = None
+        try:  # interpret line as a particle
+            tokens = re.split(r"(\s+)", input_line)
+            if len(tokens) < 25: raise ValueError
+            IDUP = int(tokens[2])
+            if abs(IDUP) == 13 or abs(IDUP) == 14:  # this is a muon or muon neutrino
+                if IDUP > 0:
+                    IDUP -= 2
+                else:
+                    IDUP += 2
+                is_event_changed = True
+                output_line = "".join("".join(tokens[:2])+str(IDUP)+"".join(tokens[3:]))
+        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)
+
+def e2mu(input_event):
+    """! 
+    Swap out electrons for muons, and electron neutrinos for muon neutrinos.
+    Note no momentum reshuffling is done.
+    """
+    is_event_changed = False
+    event_lines = ""
+    for input_line in input_event.splitlines(True):
+        output_line = None
+        try:  # interpret line as a particle
+            tokens = re.split(r"(\s+)", input_line)
+            if len(tokens) < 25: raise ValueError
+            IDUP = int(tokens[2])
+            if abs(IDUP) == 11 or abs(IDUP) == 12:  # this is an electron or electron neutrino
+                if IDUP > 0:
+                    IDUP += 2
+                else:
+                    IDUP -= 2
+                is_event_changed = True
+                output_line = "".join("".join(tokens[:2])+str(IDUP)+"".join(tokens[3:]))
+        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)
+
 
 def update_XWGTUP_with_reweighted_nominal(input_event, wgtid_for_old_XWGTUP_value = None):
     """! Ensure that XWGTUP is equal to the reweighted nominal."""