diff --git a/Generators/GeneratorFilters/share/common/xAODDiLeptonMassFilter_Common.py b/Generators/GeneratorFilters/share/common/xAODDiLeptonMassFilter_Common.py new file mode 100644 index 0000000000000000000000000000000000000000..8c5e7801b9a06c200e2d9fc46469441523395b1d --- /dev/null +++ b/Generators/GeneratorFilters/share/common/xAODDiLeptonMassFilter_Common.py @@ -0,0 +1,12 @@ +# common fragment for xAODDiLeptonMass filter +# conversion to XAOD, +# creation of slimmed container containing electrons +# connecting the filter + +include ("GeneratorFilters/CreatexAODSlimContainers.py") +createxAODSlimmedContainer("TruthLightLeptons",prefiltSeq) +prefiltSeq.xAODCnv.AODContainerName = 'GEN_EVENT' + +from GeneratorFilters.GeneratorFiltersConf import xAODDiLeptonMassFilter + + diff --git a/Generators/GeneratorFilters/share/common/xAODLeptonFilter_Common.py b/Generators/GeneratorFilters/share/common/xAODLeptonFilter_Common.py new file mode 100644 index 0000000000000000000000000000000000000000..d32dd5f5e1ba6b7c19cc5b7bed768b403fa8b348 --- /dev/null +++ b/Generators/GeneratorFilters/share/common/xAODLeptonFilter_Common.py @@ -0,0 +1,22 @@ +# common fragment for xAODLepton filter +# conversion to XAOD, +# creation of slimmed container containing electrons and muons +# connecting the filter + +if not hasattr(prefiltSeq, 'xAODCnv'): + from xAODTruthCnv.xAODTruthCnvConf import xAODMaker__xAODTruthCnvAlg + prefiltSeq += xAODMaker__xAODTruthCnvAlg('xAODCnv',WriteTruthMetaData=False) + prefiltSeq.xAODCnv.AODContainerName = 'GEN_EVENT' + +if not hasattr(prefiltSeq, "xAODTruthParticleSlimmerElectron"): + from GeneratorFilters.GeneratorFiltersConf import xAODTruthParticleSlimmerElectron + prefiltSeq += xAODTruthParticleSlimmerElectron('xAODTruthParticleSlimmerElectron') + +if not hasattr(prefiltSeq, "xAODTruthParticleSlimmerMuon"): + from GeneratorFilters.GeneratorFiltersConf import xAODTruthParticleSlimmerMuon + prefiltSeq += xAODTruthParticleSlimmerMuon('xAODTruthParticleSlimmerMuon') + +from GeneratorFilters.GeneratorFiltersConf import xAODLeptonFilter + + + diff --git a/Generators/GeneratorFilters/src/components/GeneratorFilters_entries.cxx b/Generators/GeneratorFilters/src/components/GeneratorFilters_entries.cxx index 4ec367e7a089d84a56d5e93887da7c9c7bc5cbcf..a15bb61149cd93cc76b3aafbb2b99ffe202fc5dc 100644 --- a/Generators/GeneratorFilters/src/components/GeneratorFilters_entries.cxx +++ b/Generators/GeneratorFilters/src/components/GeneratorFilters_entries.cxx @@ -8,6 +8,7 @@ #include "GeneratorFilters/xAODDiLeptonMassFilter.h" #include "GeneratorFilters/xAODDirectPhotonFilter.h" #include "GeneratorFilters/xAODElectronFilter.h" +#include "GeneratorFilters/xAODLeptonFilter.h" #include "GeneratorFilters/xAODLeptonPairFilter.h" #include "GeneratorFilters/xAODM4MuIntervalFilter.h" #include "GeneratorFilters/xAODMETFilter.h" @@ -141,6 +142,7 @@ DECLARE_COMPONENT( xAODDecayTimeFilter ) DECLARE_COMPONENT( xAODDiLeptonMassFilter ) DECLARE_COMPONENT( xAODDirectPhotonFilter ) DECLARE_COMPONENT( xAODElectronFilter ) +DECLARE_COMPONENT( xAODLeptonFilter ) DECLARE_COMPONENT( xAODLeptonPairFilter ) DECLARE_COMPONENT( xAODM4MuIntervalFilter ) DECLARE_COMPONENT( xAODMETFilter) diff --git a/Generators/Hijing_i/share/genPbPbJobOpt.py b/Generators/Hijing_i/share/genPbPbJobOpt.py deleted file mode 100755 index 571ba7c4d2ff29d8676a0b1891bb21be6085ae46..0000000000000000000000000000000000000000 --- a/Generators/Hijing_i/share/genPbPbJobOpt.py +++ /dev/null @@ -1,108 +0,0 @@ -############################################################### -# -# An example of job options file for Hijing generation of -# Pb + Pb collisions at 5520 GeV/(colliding nucleon pair) -# -# Brian Cole -# Mikhail Leltchouk -# Andrzej Olszewski -# -# February 2007 -# -# more examples in https://twiki.cern.ch/twiki/bin/view/Atlas/AtlasHICSC2007 -# -#============================================================== -# General Application Configuration options -#-------------------------------------------------------------- -include( "AthenaPoolCnvSvc/WriteAthenaPool_jobOptions.py" ) -from AthenaCommon.Configurable import Configurable -from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream - -# configuring the Athena application for a 'generator' job -import AthenaCommon.AtlasUnixGeneratorJob - -# make sure we are loading the ParticleProperty service -from PartPropSvc.PartPropSvcConf import PartPropSvc -svcMgr += PartPropSvc() -#-------------------------------------------------------------- -# Private Application Configuration options -#-------------------------------------------------------------- -from AthenaCommon.AlgSequence import AlgSequence -topAlg = AlgSequence() - -from Hijing_i.Hijing_iConf import Hijing -topAlg += Hijing() - -theApp.ExtSvc += ["AtRndmGenSvc"] - -# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL ) -svcMgr.MessageSvc.OutputLevel = 4 -svcMgr.MessageSvc.defaultLimit = 100000 -#-------------------------------------------------------------- -# Event related parameters -#-------------------------------------------------------------- -# Number of events to be processed (default is 10) -theApp.EvtMax = 5 - -# Set run number (default 0 causes problems) -svcMgr.EventSelector.RunNumber = 12345 - -#-------------------------------------------------------------- -# Algorithms Private Options -#-------------------------------------------------------------- - -import random - -# Generate seeds -# -seed1 = random.randint(2<<15-1,1<<31-1) -seed2 = random.randint(2<<15-1,1<<31-1) -seed1 &= 0xfffffffe -seed2 &= 0xfffffffe - -if not hasattr(svcMgr, 'AtRndmGenSvc'): - from AthenaServices.AthenaServicesConf import AtRndmGenSvc - svcMgr += AtRndmGenSvc() - -svcMgr.AtRndmGenSvc.Seeds = \ - ["HIJING "+str(seed1)+" "+str(seed2), "HIJING_INIT 53240261 53827126"] - -Hijing = Algorithm( "Hijing" ) -Hijing.Initialize = ["efrm 5520.", "frame CMS", "proj A", "targ A", - "iap 208", "izp 82", "iat 208", "izt 82", - # Simulation of minimum-bias events - "bmin 0.", "bmax 20.", - # Maximum pt of hard or semihard scatterings in GeV, - # if negative (or commented as below) there is no up-limit. - ### "hipr1 9 10.", - # Turns OFF jet quenching: - "ihpr2 4 0", - # Jan24,06 turns ON decays charm and bottom but not pi0, lambda, ... - "ihpr2 12 2", - # Turns ON retaining of particle history - truth information: - "ihpr2 21 1"] - -#--------------------------------------------------------------- -# Pool Persistency -#--------------------------------------------------------------- -from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream - -#--- Stream1 --- -Stream1 = AthenaPoolOutputStream( "Stream1" ) - -theApp.Dlls += [ "GeneratorObjectsAthenaPoolPoolCnv" ] -PoolSvc = Service( "PoolSvc" ) - -# 2101 == EventInfo -# 133273 == MCTruth (HepMC) -# 54790518 == HijigEventParams -Stream1.ItemList += [ "2101#*" ] -Stream1.ItemList += [ "133273#*" ] -Stream1.ItemList += [ "54790518#*" ] - -Stream1.OutputFile = "hijing.pool.root" -#============================================================== -# -# End of job options file -# -############################################################### diff --git a/Generators/Hijing_i/share/jobOptions.hijing.py b/Generators/Hijing_i/share/jobOptions.hijing.py deleted file mode 100755 index 71a73306892261a718ba7e204f18523d623828f8..0000000000000000000000000000000000000000 --- a/Generators/Hijing_i/share/jobOptions.hijing.py +++ /dev/null @@ -1,72 +0,0 @@ -############################################################### -# -# Job options file -# -#============================================================== -#-------------------------------------------------------------- -# General Application Configuration options -#-------------------------------------------------------------- -include( "AthenaPoolCnvSvc/WriteAthenaPool_jobOptions.py" ) -from AthenaCommon.Configurable import Configurable -from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream - -# configuring the Athena application for a 'generator' job -import AthenaCommon.AtlasUnixGeneratorJob - -# make sure we are loading the ParticleProperty service -from PartPropSvc.PartPropSvcConf import PartPropSvc -svcMgr += PartPropSvc() -#-------------------------------------------------------------- -# Private Application Configuration options -#-------------------------------------------------------------- -from AthenaCommon.AlgSequence import AlgSequence -topAlg = AlgSequence() - -from Hijing_i.Hijing_iConf import Hijing -topAlg += Hijing() -from TruthExamples.TruthExamplesConf import DumpMC -topAlg += DumpMC() - -theApp.ExtSvc += ["AtRndmGenSvc"] -# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL ) -svcMgr.MessageSvc.OutputLevel = 4 -svcMgr.MessageSvc.defaultLimit = 100000 -#-------------------------------------------------------------- -# Event related parameters -#-------------------------------------------------------------- -# Number of events to be processed (default is 10) -theApp.EvtMax = 5 - -# Set run number (default 0 causes problems) -svcMgr.EventSelector.RunNumber = 12345 - -#-------------------------------------------------------------- -# Algorithms Private Options -#-------------------------------------------------------------- -if not hasattr(svcMgr, 'AtRndmGenSvc'): - from AthenaServices.AthenaServicesConf import AtRndmGenSvc - svcMgr += AtRndmGenSvc() - -svcMgr.AtRndmGenSvc.Seeds = ["HIJING 4789899 989240512", "HIJING_INIT 889223465 78782321"] -#svcMgr.AtRndmGenSvc.ReadFromFile = True - -# Hijing = Algorithm( "Hijing" ) -# Hijing.Initialize = ["efrm 200", "frame CMS", "proj A", "targ A", -# "iap 197", "izp 79", "iat 197", "izt 79"] -# Hijing.wide = FALSE # Allow the random events to be anywhere in the beampipe -# Hijing.randomizeVertices = TRUE # Randomize vertices within the beampipe - default x,y = 0,0 -# Hijing.selectVertex = TRUE # Allows manual vertex settings. Please set x, y, and z. -# Hijing.x = 10.0 #mm - X vertex setting - overridden by randomization. -# Hijing.y = 10.0 #mm - Y vertex setting - overridden by randomization. -# Hijing.z = 10000.0 #mm - Z vertex setting - overridden by randomization. -# Hijing.vertexOffsetCut = 1.0E-7 # default -# Hijing.partonStoreMinPt = 5.0 # default -# Hijing.keepSpectators = true # default -#--------------------------------------------------------------- -# Ntuple service output -#--------------------------------------------------------------- -#============================================================== -# -# End of job options file -# -############################################################### diff --git a/Generators/Pythia8_i/CMakeLists.txt b/Generators/Pythia8_i/CMakeLists.txt index bd2d32c7f111ad2d0affad6eb14f43c523988093..2c997cce9d4aae376ff77559c897b7db8c14fdf3 100644 --- a/Generators/Pythia8_i/CMakeLists.txt +++ b/Generators/Pythia8_i/CMakeLists.txt @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration # Declare the package name: atlas_subdir( Pythia8_i ) @@ -47,6 +47,13 @@ atlas_add_library( Pythia8_iLib src/UserHooks/suep_shower.cxx src/UserHooks/main31.cxx src/UserHooks/EnhanceSplittings.cxx + src/UserHooks/EnhanceTest.cxx + src/UserHooks/SettableColourReconnection.cxx + src/UserHooks/PowhegV_EW.cxx + src/UserHooks/PowhegBB4L.cxx + src/UserHooks/PowhegBB4Ltms.cxx + src/UserHooks/TopRecoilHook.cxx + src/UserResonances/ResonanceExcitedCI.cxx src/UserResonances/ResonanceLQ.cxx src/Pythia8Custom/CheckForFinalPartons.cxx diff --git a/Generators/Pythia8_i/Pythia8EnvironmentConfig.cmake b/Generators/Pythia8_i/Pythia8EnvironmentConfig.cmake index 05a3d4715ad42e2356344229fdff7b3d20d6c7ef..98af526880cd7bfba62a3bf3b27ead7bc2ac40c8 100644 --- a/Generators/Pythia8_i/Pythia8EnvironmentConfig.cmake +++ b/Generators/Pythia8_i/Pythia8EnvironmentConfig.cmake @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration # # This module is used to set up the environment for Pythia8 # @@ -11,6 +11,7 @@ find_package( Lhapdf ) if( PYTHIA8_FOUND AND LHAPDF_FOUND ) set( PYTHIA8ENVIRONMENT_ENVIRONMENT FORCESET PYTHIA8VER ${PYTHIA8_LCGVERSION} + FORCESET PY8PATH ${PYTHIA8_LCGROOT} FORCESET LHAPDFVER ${LHAPDF_LCGVERSION} ) endif() diff --git a/Generators/Pythia8_i/README.md b/Generators/Pythia8_i/README.md new file mode 100644 index 0000000000000000000000000000000000000000..9048bf58d09384ddfc6af93e82ae6de9278c37ed --- /dev/null +++ b/Generators/Pythia8_i/README.md @@ -0,0 +1,148 @@ +# Pythia 8 in ATLAS + +For any issues or advice, please contact the support team via atlas-generators-pythia@cern.ch. + +Significant problems requiring code fixes should be submitted via the ATLAS +Generators JIRA: https://its.cern.ch/jira/projects/AGENE/issues/AGENE-1427 + + +## Introduction + +Pythia 8 can be run within Athena by using the `Pythia8_i` interface algorithm. +The resulting HepMC objects are put into the `StoreGate` data store, and may be +used by other packages (e.g full simulation, Atlfast, Rivet, etc.) in the usual +way. The `Pythia_i` interface code is located in [this +package](https://gitlab.cern.ch/atlas/athena/tree/21.6/Generators/Pythia8_i). + + +## Running Pythia 8 in Athena + +`Pythia8_i` can be used from any MC production release, set up with +e.g. `setupATLAS; lsetup asetup; asetup 21.6.20,AthGeneration`. Releases not +used for MC production will not have been tested, but may nevertheless work. + +Pythia8_i should be run via the `Gen_tf.py` transform script in release 21 +(formerly `Generate_tf.py` in MC15/releases 19-20, and `Generate_trf.py` in +MC12/release 17). The transform can download standard "job option" run +configuration scripts automatically, given a job-config/dataset ID number, so if +you just want to run an existing MC process just run like this: + + Gen_tf.py --ecmEnergy=13000 --jobConfig=421113 --maxEvents=10 --outputEVNTFile=test_minbias_inelastic.EVNT.pool.root + +If the process you want does not already exist, you will need to make your own +job option script and run it locally: this is covered in the following section. + + +## Writing job option files + +Athena is steered using job option scripts (JOs), written in Python. To run +"just" Pythia8, i.e. internal hard-process simulation rather than external +partonic events via LHE format (for which special settings are often required) +you should start your JO by `include()`ing a base JO fragment from +https://gitlab.cern.ch/atlas-physics/pmg/infrastructure/mc15joboptions/-/tree/master/share + + include("MC15JobOptions/Pythia8_A14_NNPDF23LO_EvtGen_Common.py") + +This file itself builds upon a non-"user-facing" set of more fundamental setup +fragments, but you should use these versions with an MC tune (A14) and parton +density set (NNPDF23LO) configured, and usually with some dedicated particle +decays handled by the EvtGen afterburner program. + +Any valid Pythia8 command, including those to set the process and generation cuts, may be passed in using + + genSeq.Pythia8.Commands += ['foo=bar'] + +The beam specifications are set on the command line by `Generate_tf.py` and do not need to be specified in the JO. + + +## Using a development version of Pythia8_i + +To use a development version not already in a release, you will need to check +out and build the relevant packages from the ATLAS Git repository. Before doing +this, follow the ATLAS Git workflow instructions at +https://atlassoftwaredocs.web.cern.ch/gittutorial/workflow-quick/ to get a +personal fork of the ATLAS codebase. Then make and a sparse checkout of the +`Pythia8_i` package with your preferred Git branch or tag: + + setupATLAS; lsetup asetup; lsetup git + git atlas init-workdir https://:@gitlab.cern.ch/USERNAME/athena.git + cd athena + git atlas addpkg Pythia8_i + git checkout -b 21.6 release/21.6.20 --no-track #< use your own preferred tag or branch + cd .. + +Now make a build directory separate from the source checkout, set up the release +to build against with your local modifications, and build & set up the run: + + mkdir athena-build && cd athena-build + asetup 21.6.20,AthGeneration + cmake -DATLAS_PACKAGE_FILTER_FILE=../package_filters.txt ../athena/Projects/WorkDir + make + source x86_64-*/setup.sh + + +## Using a custom libPythia8 library + +Using a new copy of the Pythia8 library for testing is a bit more involved, +because you will also need to check out the `Pythia8` externals package and +override some paths to point at the new version. Here is a full set of commands +to run (thanks to Marvin Flores), again starting from a new checkout of your Git +Athena fork: + + setupATLAS; lsetup asetup; lsetup git + git atlas init-workdir https://:@gitlab.cern.ch/USERNAME/athena.git athena-py8custom + cd athena-py8custom + git atlas addpkg Pythia8 Pythia8_i + +You now need to modify the project build instructions to pick up the external +library override, by editing `athena/Projects/WorkDir/CMakeLists.txt`: + +1. Just after `find_package( AtlasCMake QUIET )` add, for example, + + set( PYTHIA8_LCGVERSION 301p3 ) + set( PYTHIA8_LCGROOT /cvmfs/sft.cern.ch/lcg/releases/LCG_88/MCGenerators/pythia8/${PYTHIA8_LCGVERSION}/${LCG_PLATFORM} ) + set( PYTHIA8_VERSION ${PYTHIA8_LCGVERSION}) + set( PYTHIA8_ROOT ${PYTHIA8_LCGROOT}) + +2. Just after `# Find the project that we depend on:` add + + find_package(AthGenerationExternals REQUIRED) + +And in `athena/External/Pythia8/CMakeLists.txt`, add LCG version overrides just +before `find_package( Pythia8 )`: + + set( PYTHIA8_LCGVERSION 301p3 ) + set( PYTHIA8_LCGROOT /cvmfs/sft.cern.ch/lcg/releases/LCG_88/MCGenerators/pythia8/${PYTHIA8_LCGVERSION}/${LCG_PLATFORM}) + +Before building, you should also create a `package_filters.txt` file, perhaps +based on `athena/Projects/WorkDir/package_filters_example.txt`. This should contain the +following lines, specifying which packages to build: + + + External/Pythia8 + + Generators/Pythia8_i + - .* + +Finally, (re)build: + + mkdir athena-py8custom-build && cd athena-py8custom-build + asetup 21.6.20,AthGeneration + # or: asetup 21.6,latest,AthGeneration,slc6 + cmake -DATLAS_PACKAGE_FILTER_FILE=../package_filters.txt ../athena/Projects/WorkDir + make + source x86_64-*/setup.sh + +You may wish to run a `make clean` before `make`, to ensure that everything is +definitely rebuilt as intended. + +## MLM matching within Athena +Where CKKWL cannot be applied (notably, in all loop induced processes) it can be useful to run MLM matching. +The first efforts in implementing MLM matching within Athena are documented in this [JIRA ticket](https://its.cern.ch/jira/browse/AGENE-1879). +The MLM matching base fragment is Pythia8_i/Pythia8_ClassicalMLM_Match.py, which works in close analogy to the corresponding CKKWL base fragment. +This implement the classical, Madgraph inspired, MLM matching technique, checking on an event-by-event basis the event jet multiplicity. + +## Sacrifice standalone steering package + +In addition to the Athena interface, James Monk has written a standalone package +that steers Pythia 8 without needing Athena, and provides Photos++, LHAPDF, LHEF +and HepMC interfaces and can be more easily steered from the command line. The +package is available from the AGILe project on HepForge: https://agile.hepforge.org/ diff --git a/Generators/Pythia8_i/python/Pythia8Util.py b/Generators/Pythia8_i/python/Pythia8Util.py new file mode 100644 index 0000000000000000000000000000000000000000..5e2d15ea55ef47a57f88b6e5f430b4fac86c0735 --- /dev/null +++ b/Generators/Pythia8_i/python/Pythia8Util.py @@ -0,0 +1,137 @@ +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration + +from AthenaCommon.AppMgr import ServiceMgr as svcMgr +from GeneratorModules.EvgenAlg import EvgenAlg +from AthenaPython.PyAthena import StatusCode +import glob,subprocess,os,time,shutil +import xml.dom.minidom +import filecmp + +__author__ = 'Giancarlo Panizzo ' + +class Pythia8Util(EvgenAlg): + """ + A python-only algorithm to monitor Pythia8 behaviour / settings + """ + + def __init__(self, name="Pythia8Util", outputParticleDataFile="ParticleData.local.xml", threshold=50): + super(Pythia8Util,self).__init__(name=name) + self.threshold = threshold + + self.njobs = 1 + self.mode = 0 + self.Pythia8Commands = [] + self.outputParticleDataFile = outputParticleDataFile + + def ensure_toplevel(self,file_name): + """ Insert given string as a new line at the beginning of a file """ + # define name of temporary dummy file + dummy_file = file_name + '.bak' + # open original file in read mode and dummy file in write mode + with open(file_name, 'r') as read_obj, open(dummy_file, 'w') as write_obj: + # Write given line to the dummy file + write_obj.write(""+'\n') + # Read lines from original file one by one and append them to the dummy file + for line in read_obj: + write_obj.write(line) + write_obj.write(""+' \n') + # remove original file + os.remove(file_name) + # Rename dummy file as the original file + os.rename(dummy_file, file_name) + + + def cmpsettings(self, requestedsettings): + if ( filecmp.cmp("Settings_before.log","Settings_after.log") ): + self.msg.info("Settings match before and after initialisation of Pythia8.") + else: + if requestedsettings.has_key("PartonShowers:model"): + if (requestedsettings["PartonShowers:model"] !="1"): + self.msg.info("Settings before and after initialisation of Pythia8 don't match, because you asked PartonShowers:model = %s" % requestedsettings["PartonShowers:model"]) + else: + self.msg.warning("Settings before and after initialisation of Pythia8 don't match.") + with open("Settings_before.log") as f1: + f1_lines = f1.read().splitlines() + with open("Settings_after.log") as f2: + f2_lines = f2.read().splitlines() + # Find and print the diff: + if (len(f1_lines) == len(f2_lines)): + for iline in range(len(f1_lines)): + if (f1_lines[iline]!=f2_lines[iline]): + if ( requestedsettings.has_key(f2_lines[iline].split("=")[0]) ): + self.msg.warning(" >> %s != %s. You requested %s." % (f1_lines[iline],f2_lines[iline].split("=")[1],requestedsettings[f2_lines[iline].split("=")[0]]) ) + else: + self.msg.info(" >> %s != %s" % (f1_lines[iline],f2_lines[iline].split("=")[1]) ) + else: + self.msg.warning("Not even the -number- of settings before and after initialisation of Pythia8 matches, check yourself what's going wrong.") + + def cmpparticledata(self,pids,rpd): + # open particledata before initialisation + f1=self.outputParticleDataFile + f2=self.outputParticleDataFile + if f1.endswith('xml'): + f1 = f1[:-3] + f1 += "orig.xml" + try: + doc1 = xml.dom.minidom.parse(f1); + except: + # not yet a problem, particledata files needs to be decorated to be xml compliant + self.ensure_toplevel(f1) + try: + doc1 = xml.dom.minidom.parse(f1); + except: + self.msg.error("Bad file, exiting") + return; + + # open particledata after initialisation + try: + doc2 = xml.dom.minidom.parse(f2); + except: + # not yet a problem, particledata files needs to be decorated to be xml compliant + self.ensure_toplevel(f2) + try: + doc2 = xml.dom.minidom.parse(f2); + except: + self.msg.error("Bad file, exiting") + return; + + + particles1 = doc1.getElementsByTagName("particle") + particles2 = doc2.getElementsByTagName("particle") + if (particles1.length == particles2.length): + self.msg.info("Number of particles before and after Pythia8 initialisation matches ( %d )" % particles1.length) + else: + self.msg.warning("Mismatch in number of particles before and after Pythia8 initialisation. Before: %d , after: %d" % (particles1.length,particles2.length) ) + for part1 in particles1: + if ( part1.getAttribute("id") in pids): + for part2 in particles2: + if ( part2.getAttribute("id") == part1.getAttribute("id") ): + for attrk2, attrv2 in part2.attributes.items(): + requested="none" + trypartprop = part2.getAttribute("id")+":"+attrk2 + if ( trypartprop in rpd ): + requested = rpd[trypartprop] + if ( not ( attrk2 in part1.attributes.keys() ) ): + self.msg.info("You asked Pythia8 to modify properties for particle %s (%s). Attribute \"%s\" added after initialisation, check it. Requested: %s. After init: %s." % (part2.getAttribute("name"), part2.getAttribute("id"), attrk2, requested, attrv2) ) + for attrk1, attrv1 in part1.attributes.items(): + if (attrk1 == attrk2 and attrv1 != attrv2): + self.msg.warning("You asked Pythia8 to modify properties for particle %s (%s). Attribute \"%s\" modified after initialisation. Requested: %s. Before init: %s. After init: %s." % (part1.getAttribute("name"), part1.getAttribute("id"), attrk1, requested, attrv1, attrv2) ) + + + + def initialize(self): + self.msg.info('Running in serial mode.') + return StatusCode.Success + + def execute(self): + return StatusCode.Success + + def finalize(self): + # get list of PIDs for which user asked some changes + joparticles = [ i.replace(" ", "").split(":")[0] for i in self.Pythia8Commands if (i[0].isdigit()) ] + requestedparticledata = { i.replace(" ", "").split("=")[0]:i.split("=")[1].strip() for i in self.Pythia8Commands if (i[0].isdigit()) } + requestedsettings = { i.replace(" ", "").split("=")[0]:i.split("=")[1].strip() for i in self.Pythia8Commands if (not i[0].isdigit()) } + self.cmpsettings(requestedsettings) + self.cmpparticledata(joparticles,requestedparticledata) + return StatusCode.Success + diff --git a/Generators/Pythia8_i/python/__init__.py b/Generators/Pythia8_i/python/__init__.py new file mode 100755 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/Generators/Pythia8_i/share/common/Pythia8_FxFx_A14mod.py b/Generators/Pythia8_i/share/common/Pythia8_FxFx_A14mod.py new file mode 100644 index 0000000000000000000000000000000000000000..2da27871265fa007d488d1d256b1a2a3822c84c2 --- /dev/null +++ b/Generators/Pythia8_i/share/common/Pythia8_FxFx_A14mod.py @@ -0,0 +1,44 @@ + +# +# This version changes the tuned settings of A14 slightly in favour of +# the recommended FxFx settings +# + +try: + PYTHIA8_nJetMax +except RuntimeError: + raise RuntimeError("Variable \"PYTHIA8_nJetMax\" is not defined, this is needed to configure Pythia8 FxFx matching settings. Please define it in your jobOptions") +else: + print("PYTHIA8_nJetMax = %i"%PYTHIA8_nJetMax) + + +try: + PYTHIA8_qCut +except RuntimeError: + raise RuntimeError("Variable \"PYTHIA8_qCut\" is not defined, this is needed to configure Pythia8 FxFx matching settings. Please define it in your jobOptions") +else: + print("PYTHIA8_qCut = %i"%PYTHIA8_qCut) + + + +genSeq.Pythia8.Commands += ["JetMatching:merge = on", + "JetMatching:scheme = 1", + "JetMatching:setMad = off", + "JetMatching:qCut = %f"%PYTHIA8_qCut, + "JetMatching:coneRadius = 1.0", + "JetMatching:etaJetMax = 10.0", + "JetMatching:doFxFx = on", + "JetMatching:qCutME = 8.0", + "JetMatching:nJetMax = %i"%PYTHIA8_nJetMax, + "JetMatching:jetAlgorithm = 2", #explicit setting of kt-merging for FxFx (also imposed by Py8-FxFx inteface) + "JetMatching:slowJetPower = 1", #explicit setting of kt-merging for FxFx (also imposed by Py8-FxFx inteface) + "JetMatching:nQmatch = 5", #4 corresponds to 4-flavour scheme (no matching of b-quarks), 5 for 5-flavour scheme + "JetMatching:eTjetMin = %f"%PYTHIA8_qCut, #This is 20 in the Pythia default, it should be <= qCut + "SpaceShower:rapidityOrder = off", # FxFx + A14 prescription + "SpaceShower:pTmaxFudge = 1.0", # FxFx + A14 prescription +] + + +genSeq.Pythia8.UserHooks += [ 'JetMatchingMadgraph'] +genSeq.Pythia8.FxFxXS = True + diff --git a/Generators/Pythia8_i/share/common/Pythia8_MLM_ClassicalMatch.py b/Generators/Pythia8_i/share/common/Pythia8_MLM_ClassicalMatch.py new file mode 100644 index 0000000000000000000000000000000000000000..3bb6ae9ade44900480f2be166d523702d2047e1a --- /dev/null +++ b/Generators/Pythia8_i/share/common/Pythia8_MLM_ClassicalMatch.py @@ -0,0 +1,46 @@ +try: + PYTHIA8_nJetMax +except RuntimeError: + raise RuntimeError("Variable \"PYTHIA8_nJetMax\" is not defined, this is needed to configure Pythia8 MLM matching settings. Please define it in your jobOptions") +else: + print("PYTHIA8_nJetMax = %i" % PYTHIA8_nJetMax) + +try: + PYTHIA8_qCut +except RuntimeError: + raise RuntimeError("Variable \"qCut\" is not defined, this is needed to configure Pythia8 MLM matching settings. Please define it in your jobOptions") +else: + print("PYTHIA8_qCut = %f" % PYTHIA8_qCut) + +try: + PYTHIA8_coneRadius +except RuntimeError: + raise RuntimeError("Variable \"coneRadius\" is not defined, this is needed to configure Pythia8 MLM matching settings. Please define it in your jobOptions") +else: + print("PYTHIA8_coneRadius = %f" % PYTHIA8_coneRadius) + +try: + PYTHIA8_nQmatch +except RuntimeError: + raise RuntimeError("Variable \"PYTHIA8_nQmatch\" is not defined, this is needed to configure Pythia8 MLM matching settings. Please define it in your jobOptions") +else: + print("PYTHIA8_nQmatch = %i" % PYTHIA8_nQmatch) + +# decrease relative precision +genSeq.Pythia8.Commands += ["Check:epTolErr = 1e-2"] + +if "UserHooks" in genSeq.Pythia8.__slots__.keys(): + genSeq.Pythia8.UserHooks += ['JetMatchingMadgraph'] +else: + genSeq.Pythia8.UserHook = 'JetMatchingMadgraph' + +genSeq.Pythia8.Commands += ["JetMatching:merge = on", + "JetMatching:scheme = 1", + "JetMatching:setMad = off", # if on, following "qCut" setting can be removed + "JetMatching:nQmatch = %i" % PYTHIA8_nQmatch, + "JetMatching:qCut = %f" % PYTHIA8_qCut, + "JetMatching:eTjetMin = %f" % PYTHIA8_qCut, + "JetMatching:coneRadius = %f" % PYTHIA8_coneRadius, + "JetMatching:etaJetMax = 1000.0", + "JetMatching:nJetMax = %i" % PYTHIA8_nJetMax ] + diff --git a/Generators/Pythia8_i/share/common/Pythia8_Monash_NNPDF23LO_Common.py b/Generators/Pythia8_i/share/common/Pythia8_Monash_NNPDF23LO_Common.py index 656968fca59b9a3367d7fe4abfdf7b6386d1a549..c1fff21689e5e58c1af9c43926d33b86f375177f 100644 --- a/Generators/Pythia8_i/share/common/Pythia8_Monash_NNPDF23LO_Common.py +++ b/Generators/Pythia8_i/share/common/Pythia8_Monash_NNPDF23LO_Common.py @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration ## Config for Py8 Monash tune ## This is not the default recommended tune, but forms the base for A14 @@ -8,21 +8,11 @@ include("Pythia8_i/Pythia8_Base_Fragment.py") genSeq.Pythia8.Commands += [ "Tune:preferLHAPDF = 2", - "Tune:pp = 14"] + "Tune:pp = 14", + "PDF:pSet = LHAPDF6:NNPDF23_lo_as_0130_qed"] # set this AFTER Tune:pp, to avoid crashed due to missing internal Pythia8 PDF data] # set PDF set AFTER Tune:pp, to avoid crashes due to missing internal Pythia8 PDF data rel = os.popen("echo $AtlasVersion").read() -if rel[:2].isdigit() and int(rel[:2])<20: - ver = os.popen("cmt show versions External/Pythia8").read() - print ("Pythia8 version: " + ver) - if 'Pythia8-01' in ver[:50]: - genSeq.Pythia8.Commands += [ - "PDF:useLHAPDF=on", - "PDF:LHAPDFset = NNPDF23_lo_as_0130_qed"] - else: - genSeq.Pythia8.Commands += ["PDF:pSet = LHAPDF6:NNPDF23_lo_as_0130_qed"] -else: - genSeq.Pythia8.Commands += ["PDF:pSet = LHAPDF6:NNPDF23_lo_as_0130_qed"] evgenConfig.tune = "Monash" diff --git a/Generators/Pythia8_i/src/Pythia8_i.cxx b/Generators/Pythia8_i/src/Pythia8_i.cxx index 795a6f88b2950991b187bc9be12694a8a707f715..fa16aa271ca768fa54517a4f8f5e76aff7dfc8e1 100644 --- a/Generators/Pythia8_i/src/Pythia8_i.cxx +++ b/Generators/Pythia8_i/src/Pythia8_i.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration */ #include "Pythia8_i/Pythia8_i.h" #include "Pythia8_i/UserProcessFactory.h" diff --git a/Generators/Pythia8_i/src/UserHooks/GravFlat.cxx b/Generators/Pythia8_i/src/UserHooks/GravFlat.cxx index 5e7c4b32d64bf96a198e0ad56173357a00a2b6e7..5474be970ca77e0008abf23bee5d88cd08a6fa18 100644 --- a/Generators/Pythia8_i/src/UserHooks/GravFlat.cxx +++ b/Generators/Pythia8_i/src/UserHooks/GravFlat.cxx @@ -1,3 +1,7 @@ +/* + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +*/ + #include "UserSetting.h" #include "Pythia8_i/UserHooksFactory.h" #include "Pythia8/PhaseSpace.h" @@ -82,7 +86,7 @@ namespace Pythia8 { const double p4 = 1.45857e-08; const double p5 = -2.13169e-12; const double p6 = 1.26473e-16; - weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); + weightPL = 1/(p0+(p1*rH)+(p2*std::pow(rH,2))+(p3*std::pow(rH,3))+(p4*std::pow(rH,4))+(p5*std::pow(rH,5))+(p6*std::pow(rH,6))); } else if (rH < 5500) { // Completely Flat M=2 TeV, Above 4 TeV, Below 5.5 TeV @@ -94,7 +98,7 @@ namespace Pythia8 { const double p4 = 4.92286e-09; const double p5 = -2.1452e-13; const double p6 = 2.97112e-18; - weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); + weightPL = 1/(p0+(p1*rH)+(p2*std::pow(rH,2))+(p3*std::pow(rH,3))+(p4*std::pow(rH,4))+(p5*std::pow(rH,5))+(p6*std::pow(rH,6))); } else if (rH <= 6500) { // Completely Flat M=2 TeV, Above 5.5 TeV, Below 6.5 TeV diff --git a/Generators/Pythia8_i/src/UserHooks/PowhegBB4Ltms.cxx b/Generators/Pythia8_i/src/UserHooks/PowhegBB4Ltms.cxx new file mode 100644 index 0000000000000000000000000000000000000000..58fb347639dfbb6a7388379d3f91144636fe9cbe --- /dev/null +++ b/Generators/Pythia8_i/src/UserHooks/PowhegBB4Ltms.cxx @@ -0,0 +1,548 @@ +/* + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +*/ + +// PowhegHooksBB4L.h +// Copyright (C) 2017 Silvia Ferrario Ravasio, Tomas Jezo, Paolo Nason, Markus Seidel +// based on PowhegHooks.h by Richard Corke +// Adapted to be Athena and Pythia8_i compliant by Simone Amoroso (amoroso@cern.ch) + +#include "Pythia8/Pythia.h" +#include "UserHooksUtils.h" +#include "Pythia8_i/UserHooksFactory.h" +#include "Pythia8Plugins/PowhegHooks.h" +#include "UserSetting.h" +#include +#include +#include "boost/lexical_cast.hpp" +#include + +struct{ + int radtype; +} radtype_; + + + + +namespace Pythia8 { class PowhegBB4Ltms; } +Pythia8_UserHooks::UserHooksFactory::Creator powhegBB4LtmsCreator("PowhegBB4Ltms"); + + + +namespace Pythia8 { + using namespace std; + + //========================================================================== + // Use userhooks to veto PYTHIA emissions above the POWHEG scale. + class PowhegBB4Ltms : public PowhegHooks { + public: + + // Constructor + PowhegBB4Ltms() + : m_debug("Powheg:bb4l:DEBUG", false), + m_vetoFSREmission("Powheg:bb4l:FSREmission:veto", true), + m_dryRunFSR("Powheg:bb4l:dryRunFSR", false), + m_onlyDistance1("Powheg:bb4l:onlyDistance1", false), + m_vetoAtPL("Powheg:bb4l:vetoAtPL", false), + m_vetoQED("Powheg:bb4l:vetoQED", false), + m_vetoPartonLevel("Powheg:bb4l:PartonLevel:veto", 0), + m_wouldVetoFsr(false), + m_topresscale(-1.), m_vetoDecScale(-1.), m_atopresscale(-1.), + m_nISRveto(0), m_nFSRveto(0), m_nFSRvetoBB4l(0), + m_pTmin("Powheg:bb4l:pTminVeto", 0.5), + m_vetoTopCharge(-1), + m_vetoProduction("Powheg:veto", 0), + m_pTpythiaVeto("Powheg:bb4l:pTpythiaVeto", 0), + m_vetoDipoleFrame("Powheg:bb4l:FSREmission:vetoDipoleFrame", 0), + m_scaleResonanceVeto("Powheg:bb4l:ScaleResonance:veto", 0.) + { + std::cout << "**********************************************************" << std::endl; + std::cout << "* *" << std::endl; + std::cout << "* Applying Powheg BB4L UserHook! *" << std::endl; + std::cout << "* Must run on dedicated Powheg LHE input file *" << std::endl; + std::cout << "* (on your own responsibility) *" << std::endl; + std::cout << "* *" << std::endl; + std::cout << "**********************************************************" << std::endl; + + } + + + ~PowhegBB4Ltms() override { std::cout << "Number of FSR vetoed in BB4l = " << m_nFSRvetoBB4l << std::endl; } + + + //Initialization + bool initAfterBeams() { + // settings of the parent class + PowhegHooks::initAfterBeams(); + // settings of this class + return true; + } + + + //--- PROCESS LEVEL HOOK --------------------------------------------------- + // called at the LHE level + inline bool canVetoProcessLevel() { return true; } + inline bool doVetoProcessLevel(Event &e) { + ///////////// + // count myself final particles .... + //int count = 0; + //int nFinal = settingsPtr->mode("POWHEG:nFinal"); + //double pT1 = 0., pTsum = 0.; + //for (int i = e.size() - 1; i > 0; i--) { + // if (e[i].isFinal()) { + // count++; + // pT1 = e[i].pT(); + // pTsum += e[i].pT(); + // } else break; + //} + //// Extra check that we have the correct final state + //std::cout << "INFO: counted final particles is "<< count<< "; nFinal is "<getEventComments(); + string temp; + ss >> temp >> radtype_.radtype; + assert (temp == "#rwgt"); + + // find last top and the last anti-top in the record + int i_top = -1, i_atop = -1; + for (int i = 0; i < e.size(); i++) { + if (e[i].id() == 6) i_top = i; + if (e[i].id() == -6) i_atop = i; + } + if (i_top != -1) + m_topresscale = findresscale(i_top, e); + else + m_topresscale = 1e30; + if (i_top != -1) + m_atopresscale = findresscale(i_atop, e); + else + m_atopresscale = 1e30; + // initialize stuff + doVetoFSRInit(); + // do not veto, ever + return false; + } + + + //--- PARTON LEVEL HOOK ---------------------------------------------------- + // called after shower + bool retryPartonLevel() { return m_vetoPartonLevel(settingsPtr) || m_vetoAtPL(settingsPtr); } + inline bool canVetoPartonLevel() { return m_vetoPartonLevel(settingsPtr) || m_vetoAtPL(settingsPtr); } + inline bool doVetoPartonLevel(const Event &e) { + if(radtype_.radtype==2) + return false; + if(m_debug(settingsPtr)){ + if (m_dryRunFSR(settingsPtr) && m_wouldVetoFsr) { + double scale = getdechardness(m_vetoTopCharge, e); + std::cout << "FSRdecScale = " << m_vetoDecScale << ", PLdecScale = " << scale << ", ratio = " << m_vetoDecScale/scale << std::endl; + } + } + if (m_vetoPartonLevel(settingsPtr)) { + double topdecscale = getdechardness(1, e); + double atopdecscale = getdechardness(-1, e); + if ((topdecscale > m_topresscale) || (atopdecscale > m_atopresscale)) { + return true; + } + else + return false; + } + if (m_vetoAtPL(settingsPtr)) { + if (m_dryRunFSR(settingsPtr) && m_wouldVetoFsr) return true; + else return false; + } + return false; + } + + + + //--- Internal helper functions -------------------------------------------- + // Calculates the scale of the hardest emission from within the resonance system + // translated by Markus Seidel modified by Tomas Jezo + inline double findresscale( const int iRes, const Event& event) { + double scale = 0.; + + int nDau = event[iRes].daughterList().size(); + + if (nDau == 0) { + // No resonance found, set scale to high value + // Pythia will shower any MC generated resonance unrestricted + scale = 1e30; + } + else if (nDau < 3) { + // No radiating resonance found + scale = m_pTmin(settingsPtr); + } + else if (std::abs(event[iRes].id()) == 6) { + // Find top daughters + int idw = -1, idb = -1, idg = -1; + + for (int i = 0; i < nDau; i++) { + int iDau = event[iRes].daughterList()[i]; + if (std::abs(event[iDau].id()) == 24) idw = iDau; + if (std::abs(event[iDau].id()) == 5) idb = iDau; + if (std::abs(event[iDau].id()) == 21) idg = iDau; + } + + // Get daughter 4-vectors in resonance frame + Vec4 pw(event[idw].p()); + pw.bstback(event[iRes].p()); + + Vec4 pb(event[idb].p()); + pb.bstback(event[iRes].p()); + + Vec4 pg(event[idg].p()); + pg.bstback(event[iRes].p()); + + // Calculate scale + scale = std::sqrt(2*pg*pb*pg.e()/pb.e()); + } + else { + scale = 1e30; + } + + return scale; + } + + // The following routine will match daughters of particle `e[iparticle]` + // to an expected pattern specified via the list of expected particle PDG ID's `ids`, + // id wildcard is specified as 0 if match is obtained, the positions and the momenta + // of these particles are returned in vectors `positions` and `momenta` respectively + // if exitOnExtraLegs==true, it will exit if the decay has more particles than expected, but not less + inline bool match_decay(int iparticle, const Event &e, const vector &ids, + vector &positions, vector &momenta, bool exitOnExtraLegs = true){ + // compare sizes + if (e[iparticle].daughterList().size() != ids.size()) { + if (exitOnExtraLegs && e[iparticle].daughterList().size() > ids.size()) { + std::cout << "extra leg" << std::endl; + exit(-1); + } + return false; + } + + // compare content + for (unsigned int i = 0; i < e[iparticle].daughterList().size(); i++) { + int di = e[iparticle].daughterList()[i]; + if (ids[i] != 0 && e[di].id() != ids[i]) + return false; + } + // reset the positions and momenta vectors (because they may be reused) + positions.clear(); + momenta.clear(); + // construct the array of momenta + for (unsigned int i = 0; i < e[iparticle].daughterList().size(); i++) { + int di = e[iparticle].daughterList()[i]; + positions.push_back(di); + momenta.push_back(e[di].p()); + } + return true; + } + + inline double qSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){ + p1.bstback(pt); + p2.bstback(pt); + return std::sqrt( 2*p1*p2*p2.e()/p1.e() ); + } + + inline double gSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){ + p1.bstback(pt); + p2.bstback(pt); + return std::sqrt( 2*p1*p2*p1.e()*p2.e()/(std::pow(p1.e()+p2.e(),2)) ); + } + // Routines to calculate the pT (according to pTdefMode) in a FS splitting: + // i (radiator before) -> j (emitted after) k (radiator after) + // For the Pythia pT definition, a recoiler (after) must be specified. + // (INSPIRED BY pythia8F77_31.cc double pTpythia) + inline double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, + int RecAfterBranch) + { + + // Convenient shorthands for later + Vec4 radVec = e[RadAfterBranch].p(); + Vec4 emtVec = e[EmtAfterBranch].p(); + Vec4 recVec = e[RecAfterBranch].p(); + int radID = e[RadAfterBranch].id(); + + // Calculate virtuality of splitting + Vec4 Q(radVec + emtVec); + double Qsq = Q.m2Calc(); + + + // Mass term of radiator + double m2Rad = (std::abs(radID) >= 4 && std::abs(radID) < 7) ? + pow2(particleDataPtr->m0(radID)) : 0.; + + // z values for FSR + double z, pTnow; + // Construct 2 -> 3 variables + Vec4 sum = radVec + recVec + emtVec; + double m2Dip = sum.m2Calc(); + double x1 = 2. * (sum * radVec) / m2Dip; + double x3 = 2. * (sum * emtVec) / m2Dip; + z = x1 / (x1 + x3); + pTnow = z * (1. - z); + + // Virtuality + pTnow *= (Qsq - m2Rad); + + if (pTnow < 0.) { + infoPtr->errorMsg("Warning: pTpythia was negative"); + return -1.; + } + else + return(std::sqrt(pTnow)); + } + + + inline double getdechardness(int topcharge, const Event &e){ + int tid = 6*topcharge, wid = 24*topcharge, bid = 5*topcharge, gid = 21, wildcard = 0; + // find last top in the record + int i_top = -1; + Vec4 p_top, p_b, p_g, p_g1, p_g2; + for (int i = 0; i < e.size(); i++) + if (e[i].id() == tid) { + i_top = i; + p_top = e[i].p(); + } + if (i_top == -1) return -1.0; + + // summary of cases + // 1.) t > W b + // a.) b > 3 ... error + // b.) b > b g ... h = std::sqrt(2*p_g*p_b*p_g.e()/p_b.e()) + // c.) b > other ... h = -1 + // return h + // 2.) t > W b g + // a.) b > 3 ... error + // b.) b > b g ... h1 = std::sqrt(2*p_g*p_b*p_g.e()/p_b.e()) + // c.) b > other ... h1 = -1 + // i.) g > 3 ... error + // ii.) g > 2 ... h2 = std::sqrt(2*p_g1*p_g2*p_g1.e()*p_g2.e()/(pow(p_g1.e(),2)+pow(p_g2.e(),2))) ); + // iii.) g > other ... h2 = -1 + // return max(h1,h2) + // 3.) else ... error + + vector momenta; + vector positions; + + // 1.) t > b W + if ( match_decay(i_top, e, vector {wid, bid}, positions, momenta, false) ) { + double h; + int i_b = positions[1]; + // a.+b.) b > 3 or b > b g + if ( match_decay(i_b, e, vector {bid, gid}, positions, momenta) ) + h = qSplittingScale(e[i_top].p(), momenta[0], momenta[1]); + // c.) b > other + else + h = -1; + return h; + } + // 2.) t > b W g + else if ( match_decay(i_top, e, vector {wid, bid, gid}, positions, momenta, false) ) { + double h1, h2; + int i_b = positions[1], i_g = positions[2]; + // a.+b.) b > 3 or b > b g + if ( match_decay(i_b, e, vector {bid, gid}, positions, momenta) ) + h1 = qSplittingScale(e[i_top].p(), momenta[0], momenta[1]); + // c.) b > other + else + h1 = -1; + // i.+ii.) g > 3 or g > 2 + if ( match_decay(i_g, e, vector {wildcard, wildcard}, positions, momenta) ) + h2 = gSplittingScale(e[i_top].p(), momenta[0], momenta[1]); + // c.) b > other + else + h2 = -1; + return max(h1, h2); + } + // 3.) else + else { + std::cout << "getdechardness" << std::endl; + std::cout << "top at position " << i_top << std::endl; + std::cout << "with " << e[i_top].daughterList().size() << " daughters " << std::endl; + for (unsigned int i = 0; i < e[i_top].daughterList().size(); i++) { + int di = e[i_top].daughterList()[i]; + std::cout << "with daughter " << di << ": " << e[di].id() << std::endl; + } + exit(-1); + } + } + + + + //--- SCALE RESONANCE HOOK ------------------------------------------------- + // called before each resonance decay shower + inline bool canSetResonanceScale() { return m_scaleResonanceVeto(settingsPtr); } + // if the resonance is the (anti)top set the scale to: + // ---> (anti)top virtuality if radtype=2 + // ---> (a)topresscale otherwise + // if is not the top, set it to a big number + inline double scaleResonance(int iRes, const Event &e) { + if (e[iRes].id() == 6){ + if(radtype_.radtype == 2) + return std::sqrt(e[iRes].m2Calc()); + else + return m_topresscale; + } + else if (e[iRes].id() == -6){ + if(radtype_.radtype == 2) + return std::sqrt(e[iRes].m2Calc()); + else + return m_atopresscale; + } + else + return std::pow(10.0,30.); + } + + + //--- FSR EMISSION LEVEL HOOK ---------------------------------------------- + + // FSR veto: this should be true if we want PowhegHooksBB4l veto in decay + // OR PowhegHooks veto in production. (The virtual method + // PowhegHooks::canVetoFSREmission has been replaced by + // PowhegHooksBB4L::canVetoFSREmission, so FSR veto in production + // must be handled here. ISR and MPI veto are instead still + // handled by PowhegHooks.) + virtual inline bool canVetoFSREmission() { return m_vetoFSREmission(settingsPtr) || m_vetoProduction(settingsPtr); } + virtual inline bool doVetoFSREmission(int sizeOld, const Event &e, int iSys, bool inResonance) { + ////////////////////////////// + //VETO INSIDE THE RESONANCE // + ////////////////////////////// + if (inResonance && m_vetoFSREmission(settingsPtr)) { + + // get the participants of the splitting: the recoiler, the radiator and the emitted + int iRecAft = e.size() - 1; + int iEmt = e.size() - 2; + int iRadAft = e.size() - 3; + int iRadBef = e[iEmt].mother1(); + + // find the top resonance the radiator originates from + int iTop = e[iRadBef].mother1(); + int distance = 1; + while (std::abs(e[iTop].id()) != 6 && iTop > 0) { + iTop = e[iTop].mother1(); + distance ++; + } + if (iTop == 0) { + infoPtr->errorMsg("Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from top quark, not vetoing"); + return doVetoFSR(false,0,0); + } + int iTopCharge = (e[iTop].id()>0)?1:-1; + + // calculate the scale of the emission + double scale; + //using pythia pT definition ... + if(m_pTpythiaVeto(settingsPtr)) + scale = pTpythia(e, iRadAft, iEmt, iRecAft); + else{ //.. or using POWHEG pT definition + Vec4 pr(e[iRadAft].p()), pe(e[iEmt].p()), pt(e[iTop].p()), prec(e[iRecAft].p()), psystem; + // The computation of the POWHEG pT can be done in the top rest frame or in the dipole one. + // pdipole = pemt +prec +prad (after the emission) + // For the first emission off the top resonance pdipole = pw +pb (before the emission) = ptop + if(m_vetoDipoleFrame(settingsPtr)) + psystem = pr+pe+prec; + else + psystem = pt; + + // gluon splitting into two partons + if (e[iRadBef].id() == 21) + scale = gSplittingScale(psystem, pr, pe); + // quark emitting a gluon (or a photon) + else if (std::abs(e[iRadBef].id()) == 5 && ((e[iEmt].id() == 21) && ! m_vetoQED(settingsPtr)) ) + scale = qSplittingScale(psystem, pr, pe); + // other stuff (which we should not veto) + else { + scale = 0; + } + } + + if (iTopCharge > 0) { + if (m_onlyDistance1(settingsPtr)) { + if ( m_debug(settingsPtr) && (distance == 1) && scale > m_topresscale && ! m_wouldVetoFsr) + std::cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << std::endl; + return doVetoFSR((distance == 1) && scale > m_topresscale,scale,iTopCharge); + } + else { + if ( m_debug(settingsPtr) && scale > m_topresscale && ! m_wouldVetoFsr) + std::cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << std::endl; + return doVetoFSR(scale > m_topresscale,scale,iTopCharge); + } + } + else if (iTopCharge < 0){ + if (m_onlyDistance1(settingsPtr)){ + if ( m_debug(settingsPtr) && (distance == 1) && scale > m_atopresscale && ! m_wouldVetoFsr) + std::cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << std::endl; + return doVetoFSR((distance == 1) && scale > m_atopresscale,scale,iTopCharge); + } + else { + if ( m_debug(settingsPtr) && scale > m_topresscale && ! m_wouldVetoFsr) + std::cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << std::endl; + return doVetoFSR(scale > m_atopresscale,scale,iTopCharge); + } + } + else { + std::cout << "Bug in PohwgeHooksBB4l" << std::endl; + } + } + + ///////////////////////////////// + // VETO THE PRODUCTION PROCESS // + ///////////////////////////////// + else if(!inResonance && m_vetoProduction(settingsPtr)){ + return PowhegHooks::doVetoFSREmission(sizeOld, e, iSys, inResonance); + } + return 0; + } + + + + inline bool doVetoFSR(bool condition, double scale, int iTopCharge) { + if(radtype_.radtype==2) + return false; + if (condition) { + if (!m_wouldVetoFsr) { + m_wouldVetoFsr = true; + m_vetoDecScale = scale; + m_vetoTopCharge = iTopCharge; + } + if (m_dryRunFSR(settingsPtr)) return false; + else { + m_nFSRvetoBB4l++; + return true; + } + } + else return false; + } + + inline void doVetoFSRInit() { + m_wouldVetoFsr = false; + m_vetoDecScale = -1; + m_vetoTopCharge = 0; + } + + + private: + + + Pythia8_UserHooks::UserSetting m_debug; + Pythia8_UserHooks::UserSetting m_vetoFSREmission, m_dryRunFSR, m_onlyDistance1, m_vetoAtPL, m_vetoQED; + Pythia8_UserHooks::UserSetting m_vetoPartonLevel; + bool m_wouldVetoFsr; + double m_topresscale, m_vetoDecScale, m_atopresscale; + unsigned long int m_nISRveto, m_nFSRveto; + unsigned long int m_nFSRvetoBB4l; + + Pythia8_UserHooks::UserSetting m_pTmin; + int m_vetoTopCharge; + Pythia8_UserHooks::UserSetting m_vetoProduction, m_pTpythiaVeto, m_vetoDipoleFrame; + Pythia8_UserHooks::UserSetting m_scaleResonanceVeto; + + }; + + //========================================================================== + +} // end namespace Pythia8 diff --git a/Generators/Pythia8_i/src/UserHooks/PowhegV_EW.cxx b/Generators/Pythia8_i/src/UserHooks/PowhegV_EW.cxx new file mode 100644 index 0000000000000000000000000000000000000000..7842f1dce4c9e5a521bb4bc26654d1f6a81f50c1 --- /dev/null +++ b/Generators/Pythia8_i/src/UserHooks/PowhegV_EW.cxx @@ -0,0 +1,677 @@ +/* + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +*/ + +// This program is based in the example "main31.cc" from the Pythia8 examples, used to interface Pythia shower with Powheg events + +// Copyright (C) 2012 Torbjorn Sjostrand. +// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. +// Please respect the MCnet Guidelines, see GUIDELINES for details. + +// Adapted to be Athena and Pythia8_i compliant by giancarlo.panizzo@cern.ch + +#include "Pythia8/Pythia.h" +#include "UserHooksUtils.h" +#include "UserSetting.h" +#include "Pythia8_i/UserHooksFactory.h" + + +using namespace Pythia8; + +namespace Pythia8{ + struct si_data_type { + int pythiamatching, pytune; + bool use_photos, vetoqed, py8veto, nohad, savelhe, noQEDqopt; + }; + + struct si_event_info_type { + double xphcut; + double vetoscale_isr; + double vetoscale_fsr; + int evtnumber; + }; + + class CustomV_EW; + class PowhegV_EW; +} + +Pythia8_UserHooks::UserHooksFactory::Creator CustomV_EWCreator("CustomV_EW"); +Pythia8_UserHooks::UserHooksFactory::Creator PowhegV_EWCreator("PowhegV_EW"); + + +namespace Pythia8{ + + // Use userhooks to veto PYTHIA emissions above the POWHEG scale. + // CustomV_EW is intended to start from scalup the QED radiation from the + // resonance. It has to be used with SpaceShower:pTmaxMatch = 1 and "TimeShower:pTmaxMatch = 1" + + class CustomV_EW : public UserHooks { + + public: + + CustomV_EW() { + + std::cout<<"**********************************************************"<()["m_si_data_.vetoqed"]=true; + Pythia8_UserHooks::UserHooksFactory::userSettings()["m_si_data_.py8veto"]=true; + Pythia8_UserHooks::UserHooksFactory::userSettings()["m_si_event_info_.vetoscale_fsr"]=10.0; + + } + + ~CustomV_EW() {;} + + // Initialize settings. + + bool initAfterBeams() { + m_si_data_.vetoqed = settingsPtr->mode("m_si_data_.vetoqed"); + m_si_data_.py8veto = settingsPtr->mode("m_si_data_.py8veto"); + m_si_event_info_.vetoscale_fsr = settingsPtr->mode("m_si_event_info_.vetoscale_fsr"); + return true; + } + + + // Allow process cross section to be modified.. + + virtual bool canSetResonanceScale() { + // If we are vetoing QED emissions, and ptmaxmatch = 1, and we are using PYTHIA8 based veto, set scale for QED radiation from leptons + // (the default would be the resonance mass) + std::cout << "**** SI: Allow to set scale to veto QED emissions in PYTHIA" << std::endl; + if ((m_si_data_.vetoqed) && (m_si_data_.py8veto)) return true; + else return false; + } + + virtual double scaleResonance() { + // virtual double scaleResonance( const int iRes, const Event& event) { + // Set scale for the emissions from the resonace (FSR), equal to the scale stored in the LHE file + //std::cout << "SI in scaleResonance. Setting scale: " << si_event_info_.vetoscale_fsr << std::endl; + //std::cout << event[iRes].id(); + //std::cout << "\n"; + + return m_si_event_info_.vetoscale_fsr; + } + + virtual double EventList( const Event& event) + { + event.list(); + return false; + } + + private: + si_data_type m_si_data_; + si_event_info_type m_si_event_info_; + }; + + // PowhegV_EW is intended to veto QCD radiation (ISR and FSR) + // and QED radiation from the resonance. So it has to be used with + // SpaceShower:pTmaxMatch = 2 and "TimeShower:pTmaxMatch = 2" + + class PowhegV_EW : public UserHooks { + + public: + + // Constructor and destructor. + PowhegV_EW() { + + std::cout<<"**********************************************************"<()["m_si_data_.vetoqed"]=true; + Pythia8_UserHooks::UserHooksFactory::userSettings()["m_si_data_.py8veto"]=true; + Pythia8_UserHooks::UserHooksFactory::userSettings()["m_si_event_info_.vetoscale_fsr"]=10.0; + Pythia8_UserHooks::UserHooksFactory::userSettings()["m_si_event_info_.vetoscale_isr"]=10.0; + + }; + + ~PowhegV_EW() {} + + // Initialize settings, detailing merging strategy to use. + bool initAfterBeams() { + m_nFinal = settingsPtr->mode("POWHEG:nFinal"); + m_vetoMode = settingsPtr->mode("POWHEG:veto"); + m_vetoCount = settingsPtr->mode("POWHEG:vetoCount"); + m_pThardMode = settingsPtr->mode("POWHEG:pThard"); + m_pTemtMode = settingsPtr->mode("POWHEG:pTemt"); + m_emittedMode = settingsPtr->mode("POWHEG:emitted"); + m_pTdefMode = settingsPtr->mode("POWHEG:pTdef"); + m_MPIvetoMode = settingsPtr->mode("POWHEG:MPIveto"); + + m_si_data_.vetoqed = settingsPtr->mode("si_data_.vetoqed"); + m_si_data_.py8veto = settingsPtr->mode("si_data_.py8veto"); + m_si_event_info_.vetoscale_fsr = settingsPtr->mode("m_si_event_info_.vetoscale_fsr"); + m_si_event_info_.vetoscale_isr = settingsPtr->mode("m_si_event_info_.vetoscale_isr"); + return true; + } + + + //-------------------------------------------------------------------------- + + // Routines to calculate the pT (according to pTdefMode) in a splitting: + // ISR: i (radiator after) -> j (emitted after) k (radiator before) + // FSR: i (radiator before) -> j (emitted after) k (radiator after) + // For the Pythia pT definition, a recoiler (after) must be specified. + + // Compute the Pythia pT separation. Based on pTLund function in History.cc + double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, + int RecAfterBranch, bool FSR) { + + // Convenient shorthands for later + Vec4 radVec = e[RadAfterBranch].p(); + Vec4 emtVec = e[EmtAfterBranch].p(); + Vec4 recVec = e[RecAfterBranch].p(); + int radID = e[RadAfterBranch].id(); + + // Calculate virtuality of splitting + double sign = (FSR) ? 1. : -1.; + Vec4 Q(radVec + sign * emtVec); + double Qsq = sign * Q.m2Calc(); + + // Mass term of radiator + double m2Rad = (std::abs(radID) >= 4 && std::abs(radID) < 7) ? + pow2(particleDataPtr->m0(radID)) : 0.; + + // z values for FSR and ISR + double z, pTnow; + if (FSR) { + // Construct 2 -> 3 variables + Vec4 sum = radVec + recVec + emtVec; + double m2Dip = sum.m2Calc(); + double x1 = 2. * (sum * radVec) / m2Dip; + double x3 = 2. * (sum * emtVec) / m2Dip; + z = x1 / (x1 + x3); + pTnow = z * (1. - z); + + } else { + // Construct dipoles before/after splitting + Vec4 qBR(radVec - emtVec + recVec); + Vec4 qAR(radVec + recVec); + z = qBR.m2Calc() / qAR.m2Calc(); + pTnow = (1. - z); + } + + // Virtuality with correct sign + pTnow *= (Qsq - sign * m2Rad); + + // Can get negative pT for massive splittings + if (pTnow < 0.) { + std::cout << "Warning: pTpythia was negative" << std::endl; + return -1.; + } + +#ifdef DBGOUTPUT + std::cout << "pTpythia: rad = " << RadAfterBranch << ", emt = " + << EmtAfterBranch << ", rec = " << RecAfterBranch + << ", pTnow = " << std::sqrt(pTnow) << std::endl; +#endif + + // Return pT + return std::sqrt(pTnow); + } + + // Compute the POWHEG pT separation between i and j + // ISR: absolute pT of j + // FSR: pT of j w.r.t. to i + double pTpowheg(const Event &e, int i, int j, bool FSR) { + + // pT value for FSR and ISR + double pTnow = 0.; + if (FSR) { + // POWHEG d_ij (in CM frame). Note that the incoming beams have not + // been updated in the parton systems pointer yet (i.e. prior to any + // potential recoil). + int iInA = partonSystemsPtr->getInA(0); + int iInB = partonSystemsPtr->getInB(0); + double betaZ = - ( e[iInA].pz() + e[iInB].pz() ) / + ( e[iInA].e() + e[iInB].e() ); + Vec4 iVecBst(e[i].p()), jVecBst(e[j].p()); + iVecBst.bst(0., 0., betaZ); + jVecBst.bst(0., 0., betaZ); + + if ( e[i].id() == 21 && e[j].id() == 21) { + pTnow = std::sqrt( (iVecBst + jVecBst).m2Calc() * + iVecBst.e() * jVecBst.e() / + pow2(iVecBst.e() + jVecBst.e()) ); + } else { + pTnow = std::sqrt( (iVecBst + jVecBst).m2Calc() * + jVecBst.e() / iVecBst.e() ); + } + + } else { + // POWHEG pT_ISR is just kinematic pT + pTnow = e[j].pT(); + } + + // Check result + if (pTnow < 0.) { + std::cout << "Warning: pTpowheg was negative" << std::endl; + return -1.; + } + +#ifdef DBGOUTPUT + std::cout << "pTpowheg: i = " << i << ", j = " << j + << ", pTnow = " << pTnow << std::endl; +#endif + + return pTnow; + } + + // Calculate pT for a splitting based on m_pTdefMode. + // If j is -1, all final-state partons are tried. + // If i, k, r and xSR are -1, then all incoming and outgoing + // partons are tried. + // xSR set to 0 means ISR, while xSR set to 1 means FSR + double pTcalc(const Event &e, int i, int j, int k, int r, int xSRin) { + + + // Loop over ISR and FSR if necessary + double pTemt = -1., pTnow; + int xSR1 = (xSRin == -1) ? 0 : xSRin; + int xSR2 = (xSRin == -1) ? 2 : xSRin + 1; + for (int xSR = xSR1; xSR < xSR2; xSR++) { + // FSR flag + bool FSR = (xSR == 0) ? false : true; + + // If all necessary arguments have been given, then directly calculate. + // POWHEG ISR and FSR, need i and j. + if ((m_pTdefMode == 0 || m_pTdefMode == 1) && i > 0 && j > 0) { + pTemt = pTpowheg(e, i, j, (m_pTdefMode == 0) ? false : FSR); + + // Pythia ISR, need i, j and r. + } else if (!FSR && m_pTdefMode == 2 && i > 0 && j > 0 && r > 0) { + pTemt = pTpythia(e, i, j, r, FSR); + + // Pythia FSR, need k, j and r. + } else if (FSR && m_pTdefMode == 2 && j > 0 && k > 0 && r > 0) { + pTemt = pTpythia(e, k, j, r, FSR); + + // Otherwise need to try all possible combintations. + } else { + // Start by finding incoming legs to the hard system after + // branching (radiator after branching, i for ISR). + // Use partonSystemsPtr to find incoming just prior to the + // branching and track mothers. + int iInA = partonSystemsPtr->getInA(0); + int iInB = partonSystemsPtr->getInB(0); + while (e[iInA].mother1() != 1) { iInA = e[iInA].mother1(); } + while (e[iInB].mother1() != 2) { iInB = e[iInB].mother1(); } + + // If we do not have j, then try all final-state partons + int jNow = (j > 0) ? j : 0; + int jMax = (j > 0) ? j + 1 : e.size(); + for (; jNow < jMax; jNow++) { + + // Final-state jNow only + if ( !e[jNow].isFinal() ) continue; + + // POWHEG + if (m_pTdefMode == 0 || m_pTdefMode == 1) { + + // ISR - only done once as just kinematical pT + if (!FSR) { + pTnow = pTpowheg(e, iInA, jNow, (m_pTdefMode == 0) ? false : FSR); + if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow); + + // FSR - try all outgoing partons from system before branching + // as i. Note that for the hard system, there is no + // "before branching" information. + } else { + + int outSize = partonSystemsPtr->sizeOut(0); + for (int iMem = 0; iMem < outSize; iMem++) { + int iNow = partonSystemsPtr->getOut(0, iMem); + + // Coloured only, i != jNow and no carbon copies + if (iNow == jNow) continue; + if (jNow == e[iNow].daughter1() + && jNow == e[iNow].daughter2()) continue; + + pTnow = pTpowheg(e, iNow, jNow, (m_pTdefMode == 0) + ? false : FSR); + if (pTnow > 0.) pTemt = (pTemt < 0) + ? pTnow : min(pTemt, pTnow); + } // for (iMem) + + } // if (!FSR) + + // Pythia + } else if (m_pTdefMode == 2) { + + // ISR - other incoming as recoiler + if (!FSR) { + pTnow = pTpythia(e, iInA, jNow, iInB, FSR); + if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow); + pTnow = pTpythia(e, iInB, jNow, iInA, FSR); + if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow); + + // FSR - try all final-state coloured partons as radiator + // after emission (k). + } else { + for (int kNow = 0; kNow < e.size(); kNow++) { + if (kNow == jNow || !e[kNow].isFinal()) continue; + + // For this kNow, need to have a recoiler. + // Try two incoming. + pTnow = pTpythia(e, kNow, jNow, iInA, FSR); + if (pTnow > 0.) pTemt = (pTemt < 0) + ? pTnow : min(pTemt, pTnow); + pTnow = pTpythia(e, kNow, jNow, iInB, FSR); + if (pTnow > 0.) pTemt = (pTemt < 0) + ? pTnow : min(pTemt, pTnow); + + // Try all other outgoing. + for (int rNow = 0; rNow < e.size(); rNow++) { + if (rNow == kNow || rNow == jNow || + !e[rNow].isFinal()) continue; + pTnow = pTpythia(e, kNow, jNow, rNow, FSR); + if (pTnow > 0.) pTemt = (pTemt < 0) + ? pTnow : min(pTemt, pTnow); + } // for (rNow) + + } // for (kNow) + } // if (!FSR) + } // if (m_pTdefMode) + } // for (j) + } + } // for (xSR) + +#ifdef DBGOUTPUT + std::cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k + << ", r = " << r << ", xSR = " << xSRin + << ", pTemt = " << pTemt << std::endl; +#endif + + return pTemt; + } + + //-------------------------------------------------------------------------- + + // Extraction of m_pThard based on the incoming event. + // Assume that all the final-state particles are in a continuous block + // at the end of the event and the final entry is the POWHEG emission. + // If there is no POWHEG emission, then m_pThard is set to SCALUP. + + bool canVetoMPIStep() { return true; } + int numberVetoMPIStep() { return 1; } + bool doVetoMPIStep(int nMPI, const Event &e) { + + if (nMPI > 1) return false; + + // Find if there is a POWHEG emission. Go backwards through the + // event record until there is a non-final particle. Also sum pT and + // find pT_1 for possible MPI vetoing + int count = 0; + double pT1 = 0., pTsum = 0.; + for (int i = e.size() - 1; i > 0; i--) { + if (e[i].isFinal()) { + count++; + pT1 = e[i].pT(); + pTsum += e[i].pT(); + } else break; + } + // Extra check that we have the correct final state + if (count != m_nFinal && count != m_nFinal + 1) { + std::cout << "Error: wrong number of final state particles in event" << std::endl; + exit(1); + } + // Flag if POWHEG radiation present and index + bool isEmt = (count == m_nFinal) ? false : true; + int iEmt = (isEmt) ? e.size() - 1 : -1; + + // If there is no radiation or if m_pThardMode is 0 then set m_pThard = SCALUP. + m_pThard = -1; + // m_pThardMode is 0 + if (!isEmt || m_pThardMode == 0) { + // This sets the scale to veto emissions in the QCD shower by Pythia + // This scale is used for all emissions, except if they come from the resonance + m_pThard = m_si_event_info_.vetoscale_isr; + // Not using directly scalup, because the special file LHE (two scales) + + // If m_pThardMode is 1 then the pT of the POWHEG emission is checked against + // all other incoming and outgoing partons, with the minimal value taken + } else if (m_pThardMode == 1) { + m_pThard = pTcalc(e, -1, iEmt, -1, -1, -1); + + // If m_pThardMode is 2, then the pT of all final-state partons is checked + // against all other incoming and outgoing partons, with the minimal value + // taken + } else if (m_pThardMode == 2) { + m_pThard = pTcalc(e, -1, -1, -1, -1, -1); + } + + // Find MPI veto pT if necessary + if (m_MPIvetoMode == 1) { + m_pTMPI = (isEmt) ? pTsum / 2. : pT1; + } + +#ifdef DBGOUTPUT + std::cout << "doVetoMPIStep: Qfac = " << infoPtr->scalup() + << ", pThard = " << m_pThard << std::endl; +#endif + + // Initialise other variables + m_accepted = false; + m_nAcceptSeq = m_nISRveto = m_nFSRveto = 0; + + if(m_pThard < 0) + { + std::cout << "something wrong with pThard = " << m_pThard << std::endl; + exit(1); + } + + // Do not veto the event + return false; + } + + //-------------------------------------------------------------------------- + + // ISR veto + + bool canVetoISREmission() { return (m_vetoMode == 0) ? false : true; } + bool doVetoISREmission(int, const Event &e, int iSys) { + + // Must be radiation from the hard system, otherwise return + if (iSys != 0) return false; + + // If m_vetocount != 0 and we already have accepted 'vetoCount' emissions in a row, + // do nothing; if m_vetocount = 0 check all emissions + if (m_vetoCount != 0 && m_nAcceptSeq >= m_vetoCount) return false; + + // Pythia radiator after, emitted and recoiler after. + int iRadAft = -1, iEmt = -1, iRecAft = -1; + for (int i = e.size() - 1; i > 0; i--) { + if (iRadAft == -1 && e[i].status() == -41) iRadAft = i; + else if (iEmt == -1 && e[i].status() == 43) iEmt = i; + else if (iRecAft == -1 && e[i].status() == -42) iRecAft = i; + if (iRadAft != -1 && iEmt != -1 && iRecAft != -1) break; + } + if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) { + e.list(); + std::cout << "Error: couldn't find Pythia ISR emission" << std::endl; + exit(1); + } + + // m_pTemtMode == 0: pT of emitted w.r.t. radiator + // m_pTemtMode == 1: min(pT of emitted w.r.t. all incoming/outgoing) + // m_pTemtMode == 2: min(pT of all outgoing w.r.t. all incoming/outgoing) + int xSR = (m_pTemtMode == 0) ? 0 : -1; + int i = (m_pTemtMode == 0) ? iRadAft : -1; + int j = (m_pTemtMode != 2) ? iEmt : -1; + int k = -1; + int r = (m_pTemtMode == 0) ? iRecAft : -1; + double pTemt = pTcalc(e, i, j, k, r, xSR); + +#ifdef DBGOUTPUT + std::cout << "doVetoISREmission: pTemt = " << pTemt << std::endl; +#endif + + // Veto if pTemt > m_pThard + if (pTemt > m_pThard) { + m_nAcceptSeq = 0; + m_nISRveto++; + return true; + } + + // Else mark that an emission has been accepted and continue + m_nAcceptSeq++; + m_accepted = true; + + return false; + } + + //-------------------------------------------------------------------------- + + // FSR veto + + bool canVetoFSREmission() { return (m_vetoMode == 0) ? false : true; } + bool doVetoFSREmission(int, const Event &e, int iSys, bool inr) { + // radiation from the hard system: isys=0 + // radiation from resonances: isys!=0 and inr=1 + // MPI radiation: isys!=0 and inr=0 + + // we do not veto MPI radiation + // if we veto here gamma from resonance (inr==1), + // we do not have to use canSetResonanceScale + if (iSys != 0 && inr != 1) return false; + + // In case of radiation from resonance we veto + // This is used for ptmaxmatch = 2. + // If py8veto = 1, this method is also used to veto photons, otherwise, use external function + // force the radiation scale, m_pThard, to be equal to the one set in the LHE file + if (inr == 1) { + if ((m_si_data_.vetoqed == false) || (m_si_data_.py8veto == false)) { + return false; + } + else { + // Set scale for FSR from the resonance + m_pThard = m_si_event_info_.vetoscale_fsr; + } + } + + + // If m_vetocount != 0 and we already have accepted 'vetoCount' emissions in a row, + // do nothing; if m_vetocount = 0 check all emissions + if (m_vetoCount != 0 && m_nAcceptSeq >= m_vetoCount) return false; + + // Pythia radiator (before and after), emitted and recoiler (after) + int iRecAft = e.size() - 1; + int iEmt = e.size() - 2; + int iRadAft = e.size() - 3; + int iRadBef = e[iEmt].mother1(); + if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) || + e[iEmt].status() != 51 || e[iRadAft].status() != 51) { + e.list(); + std::cout << "Error: couldn't find Pythia FSR emission" << std::endl; + exit(1); + } + + // Behaviour based on m_pTemtMode: + // 0 - pT of emitted w.r.t. radiator before + // 1 - min(pT of emitted w.r.t. all incoming/outgoing) + // 2 - min(pT of all outgoing w.r.t. all incoming/outgoing) + int xSR = (m_pTemtMode == 0) ? 1 : -1; + + int i = (m_pTemtMode == 0) ? iRadBef : -1; + i = (m_pTdefMode == 1) ? iRadAft : iRadBef; + // using POWHEG pT definition i should be iRadAft (daugther) + int k = (m_pTemtMode == 0) ? iRadAft : -1; + int r = (m_pTemtMode == 0) ? iRecAft : -1; + + // When pTemtMode is 0 or 1, iEmt has been selected + double pTemt = 0.; + if (m_pTemtMode == 0 || m_pTemtMode == 1) { + // Which parton is emitted, based on m_emittedMode: + // 0 - Pythia definition of emitted + // 1 - Pythia definition of radiated after emission + // 2 - Random selection of emitted or radiated after emission + // 3 - Try both emitted and radiated after emission + + // j = radiator after + + int j = iRadAft; + //m_emittedMode = 0 -> j = iRadAft + 1 = iEmt + if (m_emittedMode == 0 || (m_emittedMode == 2 && rndmPtr->flat() < 0.5)) j++; + + for (int jLoop = 0; jLoop < 2; jLoop++) { + if (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR); + else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR)); + + // For m_emittedMode == 3, have tried iRadAft, now try iEmt + if (m_emittedMode != 3) break; + if (k != -1) swap(j, k); else j = iEmt; + } + + // If m_pTemtMode is 2, then try all final-state partons as emitted + } else if (m_pTemtMode == 2) { + pTemt = pTcalc(e, i, -1, k, r, xSR); + } + +#ifdef DBGOUTPUT + std::cout << "doVetoFSREmission: pTemt = " << pTemt << std::endl; +#endif + + // Veto if pTemt > m_pThard + if (pTemt > m_pThard) { + m_nAcceptSeq = 0; + m_nFSRveto++; + return true; + } + + // Else mark that an emission has been accepted and continue + m_nAcceptSeq++; + m_accepted = true; + return false; + } + + //-------------------------------------------------------------------------- + + // MPI veto + + bool canVetoMPIEmission() { return (m_MPIvetoMode == 0) ? false : true; } + bool doVetoMPIEmission(int, const Event &e) { + + if (m_MPIvetoMode == 1) { + if (e[e.size() - 1].pT() > m_pTMPI) { +#ifdef DBGOUTPUT + std::cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT() + << ", pTMPI = " << m_pTMPI << std::endl; +#endif + return true; + } + } + return false; + } + + //-------------------------------------------------------------------------- + + // Functions to return information + + int getNISRveto() { return m_nISRveto; } + int getNFSRveto() { return m_nFSRveto; } + + private: + int m_nFinal, m_vetoMode, m_vetoCount, m_pThardMode, m_pTemtMode, + m_emittedMode, m_pTdefMode, m_MPIvetoMode; + double m_pThard, m_pTMPI; + bool m_accepted; + // The number of accepted emissions (in a row) + int m_nAcceptSeq; + // Statistics on vetos + unsigned long int m_nISRveto, m_nFSRveto; + si_data_type m_si_data_; + si_event_info_type m_si_event_info_; + }; + +} // end namespace Pythia8 + diff --git a/Generators/Pythia8_i/src/UserHooks/SettableColourReconnection.cxx b/Generators/Pythia8_i/src/UserHooks/SettableColourReconnection.cxx index 13d511dcc18ded86360068802dd2f60be69b1d17..ddf9acd24bc897ebf200926fcf94de0ddfaa0050 100644 --- a/Generators/Pythia8_i/src/UserHooks/SettableColourReconnection.cxx +++ b/Generators/Pythia8_i/src/UserHooks/SettableColourReconnection.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration */ #include "Pythia8_i/UserHooksFactory.h" #include "Pythia8Plugins/ColourReconnectionHooks.h" @@ -86,7 +86,7 @@ namespace Pythia8{ private: Pythia8_UserHooks::UserSetting m_modeIn; - Pythia8_UserHooks::UserSetting m_strengthIn; + Pythia8_UserHooks::UserSetting m_strengthIn; }; diff --git a/Generators/Pythia8_i/src/UserHooks/TopRecoilHook.cxx b/Generators/Pythia8_i/src/UserHooks/TopRecoilHook.cxx new file mode 100644 index 0000000000000000000000000000000000000000..abab0e724fae1bbca1c53bd7053f537185845bcb --- /dev/null +++ b/Generators/Pythia8_i/src/UserHooks/TopRecoilHook.cxx @@ -0,0 +1,149 @@ +/* + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +*/ + +// PowhegHooksBB4L.h +// Copyright (C) 2017 Tomas Jezo, Markus Seidel, Ben Nachman + +// Author: Tomas Jezo, Markus Seidel and Ben Nachman based on +// PowhegHooks.h by Richard Corke + +#include "Pythia8_i/UserHooksFactory.h" +#include "Pythia8Plugins/PowhegHooks.h" +#include "UserSetting.h" + +namespace Pythia8 { + + // Create the user hook + class TopRecoilHook; + Pythia8_UserHooks::UserHooksFactory::Creator TopRecoilHookCreator("TopRecoilHook"); + + // Use userhooks to veto PYTHIA emissions above the POWHEG scale + class TopRecoilHook : public UserHooks {//class PowhegHooksBB4L : public PowhegHooks { + + public: + + // Constructor. + // m_doTopRecoil : eikonal correction in GW dipole on/off when no MEC applied. + // m_useOldDipole : in GW dipole, use partons before or after branching. + // m_doList : diagnostic output; set false for production runs. + TopRecoilHook(bool doTopRecoilIn=true, bool useOldDipoleIn=false, bool doListIn = false) { + m_doTopRecoil = doTopRecoilIn; + m_useOldDipole = useOldDipoleIn; + m_doList = doListIn; + // Constructor also creates some histograms for analysis inside User Hook. + m_wtCorr = new Hist("corrective weight", 100, 0., 2.); + std::cout << " enabling TopRecoilHook"<< std::endl; + + } + + // Destructor prints histogram. + ~TopRecoilHook() override { + if (m_doTopRecoil) cout << *m_wtCorr; + delete m_wtCorr; + } + + + // Initialise. Only use hook for simple showers with recoilToColoured = off. + virtual bool initAfterBeams() override { + + // Switch off if not using simple showers or if recoilToColoured = on. + bool recoilToColoured = settingsPtr->flag("TimeShower:recoilToColoured"); + if (recoilToColoured){ + m_doTopRecoil=false; + std::cout << " TopRecoilHook should not be used with RecoilToColoured=on; disabling"<< std::endl; + } + // Flag if W mass term is already accounted for (true) or not (false). + m_recoilDeadCone = settingsPtr->flag("TimeShower:recoilDeadCone"); + // All ok. + return true; + } + + // Allow a veto after an FSR emission + virtual bool canVetoFSREmission() override {return m_doTopRecoil;} + + // Access the event after an FSR emission, specifically inside top decay. + virtual bool doVetoFSREmission( int sizeOld, const Event& event, int iSys, + bool inResonance) override { + + // Check that we are inside a resonance decay. + if (!inResonance) return false; + + // Check that it is a top decay. + int iTop = partonSystemsPtr->getInRes(iSys); + if (iTop == 0 || event[iTop].idAbs() != 6) return false; + + // Skip first emission, where ME corrections are already made. + int sizeOut = partonSystemsPtr->sizeOut(iSys); + if (sizeOut == 2) return false; + + // Location of trial new particles: radiator, emitted, recoiler. + int iRad = sizeOld; + int iEmt = sizeOld + 1; + int iRec = sizeOld + 2; + + // The above partons are after emission; + // alternatively use the ones before. + if (m_useOldDipole) { + iRad = event[iRad].mother1(); + iRec = event[iRec].mother1(); + } + + // Check if newly emitted gluon matches (anti)top colour line. + if (event[iEmt].id() != 21) return false; + if (event[iTop].id() == 6) { + if (event[iEmt].col() != event[iTop].col()) return false; + } else { + if (event[iEmt].acol() != event[iTop].acol()) return false; + } + + // Recoiler should now be a W, else something is wrong. + if (event[iRec].idAbs() != 24) { + cout << " ERROR: recoiler is " << event[iRec].id() << endl; + return false; + } + + // Denominator: eikonal weight with W as recoiler. + double pRadRec = event[iRad].p() * event[iRec].p(); + double pRadEmt = event[iRad].p() * event[iEmt].p(); + double pRecEmt = event[iRec].p() * event[iEmt].p(); + double wtW = 2. * pRadRec / (pRadEmt * pRecEmt) + - pow2(event[iRad].m() / pRadEmt); + // If m_recoilDeadCone = on, include W mass term in denominator. + if (m_recoilDeadCone) wtW -= pow2(event[iRec].m() / pRecEmt); + + // Numerator: eikonal weight with top as recoiler. + double pRadTop = event[iRad].p() * event[iTop].p(); + double pTopEmt = event[iTop].p() * event[iEmt].p(); + double wtT = 2. * pRadTop / (pRadEmt * pTopEmt) + - pow2(event[iRad].m() / pRadEmt) - pow2(event[iTop].m() / pTopEmt); + + // Histogram weight ratio. + m_wtCorr->fill( wtT / wtW ); + + // List relevant properties. + if (m_doList) { + std::cout << "\n now event with sizeOld = " << sizeOld << ", iSys = " + << iSys << ", sizeOut = " << sizeOut << scientific + << setprecision(3) + << ", weight with W = " << wtW << " and with t = " << wtT << std::endl; + partonSystemsPtr->list(); + event.list(); + } + + // Accept/reject emission. Smooth suppression or step function. + return (wtT < wtW * rndmPtr->flat()); + } + + +private: + + // Options and Histograms. + bool m_doTopRecoil, m_useOldDipole, m_doList, m_recoilDeadCone; + Hist *m_wtCorr; + + }; + +//========================================================================== + +} // end namespace Pythia8 diff --git a/Generators/Pythia8_i/src/UserHooks/UserHooksUtils.h b/Generators/Pythia8_i/src/UserHooks/UserHooksUtils.h index 0d73d3125218d9f983d6f77c76ba8c272f9b6db1..ec5417b2e71f37d0ad331107a4b0dbb2ef33488a 100644 --- a/Generators/Pythia8_i/src/UserHooks/UserHooksUtils.h +++ b/Generators/Pythia8_i/src/UserHooks/UserHooksUtils.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration */ #ifndef PYTHIA8_USERHOOKS_USERHOOKSUTILS_H @@ -48,7 +48,7 @@ namespace Pythia8_UserHooks{ * \return the pT of \param leg relative to \param comparison */ inline double pTLeg(const Pythia8::Particle &leg, const Pythia8::Particle &comparison){ - return sqrt(pT2Leg(leg, comparison)); + return std::sqrt(pT2Leg(leg, comparison)); } /** @@ -56,7 +56,7 @@ namespace Pythia8_UserHooks{ * \param comparedIndex in Event \param evt */ inline double pTLeg(size_t legIndex, size_t comparedIndex, const Pythia8::Event &evt){ - return sqrt(pT2Leg(legIndex, comparedIndex, evt)); + return std::sqrt(pT2Leg(legIndex, comparedIndex, evt)); } /** diff --git a/Generators/TruthIO/CMakeLists.txt b/Generators/TruthIO/CMakeLists.txt index 86a910ee1af5590e8c9e1ef59853bc7fc20b048e..4864859a832127ccfa7008c4810becf2e8b10dd7 100644 --- a/Generators/TruthIO/CMakeLists.txt +++ b/Generators/TruthIO/CMakeLists.txt @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration # Declare the package name: atlas_subdir( TruthIO ) @@ -14,4 +14,4 @@ atlas_add_component( TruthIO LINK_LIBRARIES ${HEPPDT_LIBRARIES} AtlasHepMCLib AtlasHepMCfioLib AthenaBaseComps GaudiKernel GeneratorModulesLib StoreGateLib EventInfo GeneratorObjects ) # Install files from the package: -atlas_install_joboptions( share/*.py ) +atlas_install_joboptions( share/common/*.py ) diff --git a/Generators/TruthIO/share/common/HepMCReadFromFile_Common.py b/Generators/TruthIO/share/common/HepMCReadFromFile_Common.py new file mode 100644 index 0000000000000000000000000000000000000000..6be521368094b5732b913303e2ef49aa0c190f07 --- /dev/null +++ b/Generators/TruthIO/share/common/HepMCReadFromFile_Common.py @@ -0,0 +1,6 @@ +from TruthIO.TruthIOConf import HepMCReadFromFile +genSeq += HepMCReadFromFile() +genSeq.HepMCReadFromFile.InputFile="events.hepmc" +evgenConfig.generators += ["HepMCAscii"] + + diff --git a/Generators/TruthIO/share/jobOptions.read_hepmc.py b/Generators/TruthIO/share/jobOptions.read_hepmc.py deleted file mode 100644 index c3319319412b8520ab9af0518cf777dacff48da3..0000000000000000000000000000000000000000 --- a/Generators/TruthIO/share/jobOptions.read_hepmc.py +++ /dev/null @@ -1,50 +0,0 @@ -############################################################### -# -# Job options file -# -#============================================================== -#-------------------------------------------------------------- -# General Application Configuration options -#-------------------------------------------------------------- -import AthenaCommon.AtlasUnixGeneratorJob - -from AthenaCommon.AppMgr import theApp -from AthenaCommon.AppMgr import ServiceMgr - -#-------------------------------------------------------------- -# Private Application Configuration options -#-------------------------------------------------------------- -from AthenaCommon.AlgSequence import AlgSequence -job=AlgSequence() -from TruthIO.TruthIOConf import HepMCReadFromFile -job += HepMCReadFromFile() - -# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL ) -ServiceMgr.MessageSvc.OutputLevel = INFO -#-------------------------------------------------------------- -# Event related parameters -#-------------------------------------------------------------- -# Number of events to be processed (default is 10) -theApp.EvtMax = 5 -#-------------------------------------------------------------- -# Algorithms Private Options -#-------------------------------------------------------------- - -re = job.HepMCReadFromFile -re.AsciiFile= "hepmc.out" - -from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream -Stream1 = AthenaPoolOutputStream( "Stream1" ) -Stream1.ItemList += [ "2101#*", "133273#*" ] -# This line is not really needed. Just to show that PoolSvc is already covered by the ServiceMgr -PoolSvc = ServiceMgr.PoolSvc -Stream1.OutputFile = "McEvent.root" - -#--------------------------------------------------------------- -# Ntuple service output -#--------------------------------------------------------------- -#============================================================== -# -# End of job options file -# -###############################################################