Commit dbe0ba20 authored by Ewelina Maria Lobodzinska's avatar Ewelina Maria Lobodzinska
Browse files

Merge branch 'jk_updatepowheg' into '21.6'

improve PowhegControl VBF W,Z interface AGENE-1831

See merge request !45492
parents 6dff052b 42cd92a5
......@@ -10,6 +10,7 @@ from output_file_renamer import output_file_renamer
from output_tarball_preparer import output_tarball_preparer
from photos import PHOTOS
from quark_colour_fixer import quark_colour_fixer
from mu2tau import mu2tau
from reweighter import reweighter
from lhe_cleaner import lhe_cleaner
from lhe_nominal_weight_updater import lhe_nominal_weight_updater
# Copyright (C) 2002-2017 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("mu2tau")
def mu2tau(powheg_LHE_output):
"""! Post-process existing events from muons 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 muon 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.mu2tau(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, "{}.mu2tau_backup".format(powheg_LHE_output))
shutil.move(powheg_LHE_tau, powheg_LHE_output)
......@@ -32,6 +32,7 @@ class Scheduler(object):
"PHOTOS", "reweighter",
"NNLO reweighter",
"LHE file nominal weight updater",
"mu2tau",
"MadSpin",
"integration grid tester",
"cross section calculator",
......@@ -69,6 +70,7 @@ class Scheduler(object):
"PHOTOS": partial(postprocessors.PHOTOS, powheg_LHE_output=powheg_LHE_output),
"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),
"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),
}
......
......@@ -49,6 +49,13 @@ class ExternalVBFNLO(ExternalBase):
# Convert allowed decay mode into PROC_ID/DECAYMODE
__vector_boson_type = self.decay_mode.split(">")[0].strip()
__vector_boson_decay = self.decay_mode.split(" ")[2].replace("+", "").replace("-", "")
if __vector_boson_type in ["w+", "w-", "z"] and __vector_boson_decay == "tau":
logger.warning("Powheg/VBFNLO does support directly tau decays for VBF W, Z production")
logger.warning("Ask to generate muon decays and hack the LHE files - make sure to validate!")
__vector_boson_decay = "mu"
process.add_algorithm("mu2tau")
# Add runcard entries
VBF_runcard_entries = [("PROC_ID", {"z": 120, "w+": 130, "w-": 140}[__vector_boson_type], "{} boson".format(__vector_boson_type)),
("DECAYMODE", {"e": 11, "mu": 13}[__vector_boson_decay], self.decay_mode),
......
......@@ -21,7 +21,8 @@ class VBF_W(PowhegV2):
super(VBF_W, self).__init__(base_directory, "VBF_W-Z", **kwargs)
# Add algorithms to the sequence
self.add_algorithm(ExternalVBFNLO("W", ["w+ > e+ ve", "w+ > mu+ vm", "w- > e- ve~", "w- > mu- vm~"]))
self.add_algorithm(ExternalVBFNLO("W", ["w+ > e+ ve", "w+ > mu+ vm", "w- > e- ve~", "w- > mu- vm~",\
"w+ > tau+ vt", "w- > tau- vt~"]))
# Add all keywords for this process, overriding defaults if required
self.add_keyword("bornktmin")
......@@ -60,6 +61,7 @@ class VBF_W(PowhegV2):
self.add_keyword("par_2gsupp")
self.add_keyword("par_diexp")
self.add_keyword("par_dijexp")
self.add_keyword("parallelstage")
self.add_keyword("pdfreweight")
self.add_keyword("Phasespace")
self.add_keyword("ptj_gencut")
......
......@@ -21,7 +21,7 @@ class VBF_Z(PowhegV2):
super(VBF_Z, self).__init__(base_directory, "VBF_W-Z", **kwargs)
# Add algorithms to the sequence
self.add_algorithm(ExternalVBFNLO("Z", ["z > e+ e-", "z > mu+ mu-"]))
self.add_algorithm(ExternalVBFNLO("Z", ["z > e+ e-", "z > mu+ mu-", "z > tau+ tau-"]))
# Add all keywords for this process, overriding defaults if required
self.add_keyword("bornktmin")
......@@ -38,9 +38,9 @@ class VBF_Z(PowhegV2):
self.add_keyword("facscfact", self.default_scales[0])
self.add_keyword("fakevirt")
self.add_keyword("flg_debug")
self.add_keyword("foldcsi")
self.add_keyword("foldphi")
self.add_keyword("foldy")
self.add_keyword("foldcsi", 2)
self.add_keyword("foldphi", 2)
self.add_keyword("foldy", 2)
self.add_keyword("hdamp")
self.add_keyword("hfact")
self.add_keyword("icsimax")
......@@ -54,13 +54,14 @@ class VBF_Z(PowhegV2):
self.add_keyword("lhans1", self.default_PDFs)
self.add_keyword("lhans2", self.default_PDFs)
self.add_keyword("manyseeds")
self.add_keyword("mll_gencut", 20)
self.add_keyword("mll_gencut", 40)
self.add_keyword("ncall1", 1000000)
self.add_keyword("ncall2", 3000000)
self.add_keyword("nubound", 1500000)
self.add_keyword("par_2gsupp")
self.add_keyword("par_diexp")
self.add_keyword("par_dijexp")
self.add_keyword("parallelstage")
self.add_keyword("pdfreweight")
self.add_keyword("Phasespace")
self.add_keyword("ptj_gencut")
......
......@@ -196,10 +196,35 @@ def ensure_coloured_quarks(input_event):
event_lines += output_line if output_line is not None else input_line
return (is_event_changed, event_lines)
def mu2tau(input_event):
"""!
Swap out muons for taus, and muon 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) == 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:]))
# print "LHE mu2tau swap \n\'"+input_line+"\', for \n\'"+output_line+"\'"
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."""
initial_colour_flow, is_event_changed = -1, False
event_lines = ""
rwgt_nominal = None
XWGTUP = None
......
......@@ -28,3 +28,5 @@ PowhegConfig.generate()
#--------------------------------------------------------------
include("Pythia8_i/Pythia8_A14_NNPDF23LO_EvtGen_Common.py")
include("Pythia8_i/Pythia8_Powheg_Main31.py")
genSeq.Pythia8.Commands += ['Powheg:NFinal = 3']
genSeq.Pythia8.Commands += ['SpaceShower:dipoleRecoil = on']
......@@ -28,3 +28,5 @@ PowhegConfig.generate()
#--------------------------------------------------------------
include("Pythia8_i/Pythia8_A14_NNPDF23LO_EvtGen_Common.py")
include("Pythia8_i/Pythia8_Powheg_Main31.py")
genSeq.Pythia8.Commands += ['Powheg:NFinal = 3']
genSeq.Pythia8.Commands += ['SpaceShower:dipoleRecoil = on']
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment