From 8c83cd4df59acdd3d11bfb13890bd1247e3e448c Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Sat, 11 Jun 2022 13:31:41 -0700
Subject: [PATCH] Handle heavy flavor and tau decays in G4 with Pythia8
---
.../Pythia8Decayer/CMakeLists.txt | 31 ++++
.../python/Pythia8DecayerConfigNew.py | 11 ++
.../Pythia8Decayer/src/Pythia8Decayer.cxx | 136 ++++++++++++++++++
.../Pythia8Decayer/src/Pythia8Decayer.h | 63 ++++++++
.../src/Pythia8DecayerPhysicsTool.cxx | 108 ++++++++++++++
.../src/Pythia8DecayerPhysicsTool.h | 43 ++++++
.../src/components/Pythia8Decayer_entries.cxx | 4 +
.../python/G4FaserServicesConfigNew.py | 6 +-
8 files changed, 400 insertions(+), 2 deletions(-)
create mode 100644 Simulation/G4Extensions/Pythia8Decayer/CMakeLists.txt
create mode 100644 Simulation/G4Extensions/Pythia8Decayer/python/Pythia8DecayerConfigNew.py
create mode 100644 Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.cxx
create mode 100644 Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.h
create mode 100644 Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.cxx
create mode 100644 Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.h
create mode 100644 Simulation/G4Extensions/Pythia8Decayer/src/components/Pythia8Decayer_entries.cxx
diff --git a/Simulation/G4Extensions/Pythia8Decayer/CMakeLists.txt b/Simulation/G4Extensions/Pythia8Decayer/CMakeLists.txt
new file mode 100644
index 00000000..901e0e3e
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/CMakeLists.txt
@@ -0,0 +1,31 @@
+################################################################################
+# Package: Pythia8Decayer
+################################################################################
+
+# Declare the package name:
+atlas_subdir( Pythia8Decayer )
+
+#if( NOT GENERATIONBASE )
+
+ # External dependencies:
+ find_package( CLHEP )
+ find_package( Geant4 )
+ find_package( XercesC )
+ find_package( Pythia8 )
+
+ # Component(s) in the package:
+ atlas_add_library( Pythia8DecayerLib
+ src/*.cxx
+ NO_PUBLIC_HEADERS
+ PRIVATE_INCLUDE_DIRS ${GEANT4_INCLUDE_DIRS} ${XERCESC_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${PYTHIA8_INCLUDE_DIRS}
+ PRIVATE_LINK_LIBRARIES ${GEANT4_LIBRARIES} ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} ${PYTHIA8_LIBRARIES} GaudiKernel AthenaBaseComps G4AtlasInterfaces G4AtlasToolsLib Pythia8_iLib )
+
+ atlas_add_component( Pythia8Decayer
+ src/components/*.cxx
+ INCLUDE_DIRS ${GEANT4_INCLUDE_DIRS} ${XERCESC_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${PYTHIA8_INCLUDE_DIRS}
+ LINK_LIBRARIES ${GEANT4_LIBRARIES} ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} ${PYTHIA8_LIBRARIES} Pythia8DecayerLib AthenaBaseComps G4AtlasInterfaces )
+#endif()
+
+# Install files from the package:
+atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
+#atlas_install_joboptions( share/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} --extend-ignore=F401,F821 )
diff --git a/Simulation/G4Extensions/Pythia8Decayer/python/Pythia8DecayerConfigNew.py b/Simulation/G4Extensions/Pythia8Decayer/python/Pythia8DecayerConfigNew.py
new file mode 100644
index 00000000..6085f5cc
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/python/Pythia8DecayerConfigNew.py
@@ -0,0 +1,11 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+
+from AthenaConfiguration.ComponentFactory import CompFactory
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaCommon import Logging
+import os
+
+def Pythia8DecayerPhysicsToolCfg(flags, name='Pythia8DecayerPhysicsTool', **kwargs):
+ result = ComponentAccumulator()
+ result.setPrivateTools( CompFactory.Pythia8DecayerPhysicsTool(name,**kwargs) )
+ return result
diff --git a/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.cxx b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.cxx
new file mode 100644
index 00000000..7047c554
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.cxx
@@ -0,0 +1,136 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// Header for my class
+#include "Pythia8Decayer.h"
+
+// The actual decayers. Singleton classes, not toolhandles
+// #include "Pythia8Instance.h"
+
+// For passing things around
+#include "CLHEP/Vector/LorentzVector.h"
+#include "G4Track.hh"
+#include "G4DynamicParticle.hh"
+#include "G4ParticleTable.hh"
+#include "G4DecayProducts.hh"
+
+Pythia8Decayer::Pythia8Decayer( const std::string s )
+ : G4VExtDecayer(s)
+{
+ // In the constructor, make a decayer instance, so that it's initialized here and not in the event loop
+ std::string docstring = Pythia8_i::xmlpath();
+ m_decayer = std::make_unique<Pythia8::Pythia>(docstring);
+
+ m_decayer->readString("ProcessLevel:all = off");
+
+ m_decayer->readString("ProcessLevel:resonanceDecays=on");
+
+ // shut off Pythia8 (default) verbosty
+ //
+ m_decayer->readString("Init:showAllSettings=false");
+ m_decayer->readString("Init:showChangedSettings=false");
+ m_decayer->readString("Init:showAllParticleData=false");
+ m_decayer->readString("Init:showChangedParticleData=false");
+ //
+ // specify how many Py8 events to print out, at either level
+ // in this particular case print out a maximum of 10 events
+ //
+ m_decayer->readString("Next:numberShowProcess = 0" );
+ m_decayer->readString("Next:numberShowEvent = 10");
+
+ m_decayer->init();
+
+ // shut off decays of pi0's as we want Geant4 to handle them
+ // if other immediate decay products should be handled by Geant4,
+ // their respective decay modes should be shut off as well
+ //
+ m_decayer->readString("111:onMode = off");
+}
+
+G4DecayProducts* Pythia8Decayer::ImportDecayProducts(const G4Track& aTrack){
+
+ m_decayer->event.reset();
+
+ G4DecayProducts* dproducts = nullptr;
+
+ G4ParticleDefinition* pd = aTrack.GetDefinition();
+ int pdgid = pd->GetPDGEncoding();
+
+ // check if pdgid is consistent with Pythia8 particle table
+ //
+ if ( !m_decayer->particleData.findParticle( pdgid ) )
+ {
+ G4cout << " can NOT find pdgid = " << pdgid
+ << " in Pythia8::ParticleData" << G4endl;
+ return dproducts;
+ }
+
+ if ( !m_decayer->particleData.canDecay(pdgid) )
+ {
+ G4cout << " Particle of pdgid = " << pdgid
+ << " can NOT be decayed by Pythia8" << G4endl;
+ return dproducts;
+ }
+
+ // NOTE: Energy should be in GeV
+
+ m_decayer->event.append( pdgid, 1, 0, 0,
+ aTrack.GetMomentum().x() / CLHEP::GeV,
+ aTrack.GetMomentum().y() / CLHEP::GeV,
+ aTrack.GetMomentum().z() / CLHEP::GeV,
+ aTrack.GetDynamicParticle()->GetTotalEnergy() / CLHEP::GeV,
+ pd->GetPDGMass() / CLHEP::GeV );
+
+ // specify polarization, if any
+
+ // NOTE: while in Py8 polarization is a double variable ,
+ // in reality it's expected to be -1, 0., or 1 in case of "external" tau's,
+ // similar to LHA SPINUP; see Particle Decays, Hadron and Tau Decays in docs at
+ // https://pythia.org/manuals/pythia8305/Welcome.html
+ // so it's not able to handle anything like 0.99, thus we're rounding off
+ m_decayer->event.back().pol(
+ round( std::cos( aTrack.GetPolarization().angle( aTrack.GetMomentumDirection() ) )
+ )
+ );
+
+ int npart_before_decay = m_decayer->event.size();
+
+ m_decayer->next();
+
+ int npart_after_decay = m_decayer->event.size();
+
+ // create & fill up decay products
+ //
+ dproducts = new G4DecayProducts(*(aTrack.GetDynamicParticle()));
+
+ // create G4DynamicParticle out of each m_decayer->event entry (except the 1st one)
+ // and push into dproducts
+
+ for ( int ip=npart_before_decay; ip<npart_after_decay; ++ip )
+ {
+
+ // only select final state decay products (direct or via subsequent decays);
+ // skip all others
+ //
+ // NOTE: in general, final state decays products will have
+ // positive status code between 91 and 99
+ // (in case such information could be of interest in the future)
+ //
+ if ( m_decayer->event[ip].status() < 0 ) continue;
+
+ // should we also skip neutrinos ???
+ // if so, skip products with abs(m_decayer->event[ip].id()) of 12, 14, or 16
+
+ G4ParticleDefinition* pddec =
+ G4ParticleTable::GetParticleTable()->FindParticle( m_decayer->event[ip].id() );
+ if ( !pddec ) continue; // maybe we should print out a warning !
+ G4ThreeVector momentum = G4ThreeVector( m_decayer->event[ip].px() * CLHEP::GeV,
+ m_decayer->event[ip].py() * CLHEP::GeV,
+ m_decayer->event[ip].pz() * CLHEP::GeV );
+ dproducts->PushProducts( new G4DynamicParticle( pddec, momentum) );
+ }
+ dproducts->DumpInfo();
+ return dproducts;
+
+}
diff --git a/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.h b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.h
new file mode 100644
index 00000000..61b52afc
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.h
@@ -0,0 +1,63 @@
+//
+// ********************************************************************
+// * License and Disclaimer *
+// * *
+// * The Geant4 software is copyright of the Copyright Holders of *
+// * the Geant4 Collaboration. It is provided under the terms and *
+// * conditions of the Geant4 Software License, included in the file *
+// * LICENSE and available at http://cern.ch/geant4/license . These *
+// * include a list of copyright holders. *
+// * *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work make any representation or warranty, express or implied, *
+// * regarding this software system or assume any liability for its *
+// * use. Please see the license in the file LICENSE and URL above *
+// * for the full disclaimer and the limitation of liability. *
+// * *
+// * This code implementation is the result of the scientific and *
+// * technical work of the GEANT4 collaboration. *
+// * By using, copying, modifying or distributing the software (or *
+// * any work based on the software) you agree to acknowledge its *
+// * use in resulting scientific publications, and indicate your *
+// * acceptance of all terms of the Geant4 Software license. *
+// ********************************************************************
+//
+//
+/// \file eventgenerator/pythia/pythia8decayer/include/Py8Decayer.hh
+/// \brief Definition of the Py8Decayer class
+///
+/// \author J. Yarba; FNAL
+///
+
+#ifndef Pythia8Decayer_H
+#define Pythia8Decayer_H
+
+#include "G4VExtDecayer.hh"
+#include "globals.hh"
+
+#include "Pythia8_i/Pythia8_i.h"
+
+class G4Track;
+class G4DecayProducts;
+
+class Pythia8Decayer : public G4VExtDecayer
+{
+
+ public:
+
+ //ctor & dtor
+ Pythia8Decayer( const std::string s );
+
+ virtual G4DecayProducts* ImportDecayProducts(const G4Track&);
+
+ private:
+
+ // data members
+ std::unique_ptr<Pythia8::Pythia> m_decayer;
+
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
diff --git a/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.cxx b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.cxx
new file mode 100644
index 00000000..750272b2
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.cxx
@@ -0,0 +1,108 @@
+/*
+ Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+*/
+
+// class header
+#include "Pythia8DecayerPhysicsTool.h"
+// package headers
+#include "Pythia8Decayer.h"
+// Geant4 headers
+#include "globals.hh"
+#include "G4ParticleTable.hh"
+#include "G4VProcess.hh"
+#include "G4ProcessManager.hh"
+#include "G4Decay.hh"
+#include "G4DecayTable.hh"
+
+// STL headers
+#include <string>
+
+//=============================================================================
+// Standard constructor, initializes variables
+//=============================================================================
+Pythia8DecayerPhysicsTool::Pythia8DecayerPhysicsTool( const std::string& type,
+ const std::string& nam,
+ const IInterface* parent )
+ : base_class ( type, nam , parent )
+{
+}
+
+//=============================================================================
+// Destructor
+//=============================================================================
+
+Pythia8DecayerPhysicsTool::~Pythia8DecayerPhysicsTool()
+{
+}
+
+//=============================================================================
+// Initialize
+//=============================================================================
+StatusCode Pythia8DecayerPhysicsTool::initialize( )
+{
+ ATH_MSG_DEBUG("Pythia8DecayerPhysicsTool initialize( )");
+ this->SetPhysicsName(name());
+ return StatusCode::SUCCESS;
+}
+
+Pythia8DecayerPhysicsTool* Pythia8DecayerPhysicsTool::GetPhysicsOption()
+{
+ return this;
+}
+
+
+void Pythia8DecayerPhysicsTool::ConstructParticle()
+{
+ // nothing to do here
+}
+void Pythia8DecayerPhysicsTool::ConstructProcess()
+{
+ ATH_MSG_DEBUG("Pythia8DecayerPhysicsTool::ConstructProcess() called ");
+
+ Pythia8Decayer* extDecayer = new Pythia8Decayer(name());
+
+ auto particleIterator=GetParticleIterator();
+ particleIterator->reset();
+ while ((*particleIterator)())
+ {
+ G4ParticleDefinition* particle = particleIterator->value();
+ G4ProcessManager* pmanager = particle->GetProcessManager();
+ G4ProcessVector* processVector = pmanager->GetProcessList();
+ for ( size_t i=0; i<processVector->length(); ++i )
+ {
+ G4Decay* decay = dynamic_cast<G4Decay*>((*processVector)[i]);
+ if ( decay )
+ {
+ // remove native/existing decay table for
+ // a)tau's
+ // b) B+/-
+ // and replace with external decayer
+ if ( std::abs(particle->GetPDGEncoding()) == 15 ||
+ std::abs(particle->GetPDGEncoding()) == 521 )
+ {
+ if ( particle->GetDecayTable() )
+ {
+ delete particle->GetDecayTable();
+ particle->SetDecayTable(nullptr);
+ }
+ decay->SetExtDecayer(extDecayer);
+ G4cout << "Setting ext decayer for: "
+ << particleIterator->value()->GetParticleName()
+ << G4endl;
+ }
+ // now set external decayer to all particles
+ // that don't yet have a decay table
+ if ( !particle->GetDecayTable() )
+ {
+ decay->SetExtDecayer(extDecayer);
+ G4cout << "Setting ext decayer for: "
+ << particleIterator->value()->GetParticleName()
+ << G4endl;
+ }
+ }
+ }
+ }
+
+ return;
+
+}
diff --git a/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.h b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.h
new file mode 100644
index 00000000..9acd13bd
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.h
@@ -0,0 +1,43 @@
+/*
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS and FASER collaborations
+*/
+
+#ifndef PYTHIA8DECAYER_PYTHIA8DECAYERPHYSICSTOOL_H
+#define PYTHIA8DECAYER_PYTHIA8DECAYERPHYSICSTOOL_H
+
+// Include files
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "G4AtlasInterfaces/IPhysicsOptionTool.h"
+#include "G4VPhysicsConstructor.hh"
+
+/** @class Pythia8DecayerPhysicsTool Pythia8DecayerPhysicsTool.h "Pythia8Decayer/Pythia8DecayerPhysicsTool.h"
+ *
+ *
+ *
+ * @author Edoardo Farina (modified for FASER by D. Casper)
+ * @date 2015-05-14
+ */
+class Pythia8DecayerPhysicsTool : public G4VPhysicsConstructor, public extends<AthAlgTool, IPhysicsOptionTool>
+{
+public:
+ /// Standard constructor
+ Pythia8DecayerPhysicsTool( const std::string& type , const std::string& name,
+ const IInterface* parent ) ;
+
+ virtual ~Pythia8DecayerPhysicsTool( ); ///< Destructor
+
+ /// Initialize method
+ virtual StatusCode initialize( ) ;
+ virtual void ConstructParticle();
+ virtual void ConstructProcess();
+
+ /** Implements
+ */
+
+ virtual Pythia8DecayerPhysicsTool* GetPhysicsOption();
+
+};
+
+
+
+#endif
diff --git a/Simulation/G4Extensions/Pythia8Decayer/src/components/Pythia8Decayer_entries.cxx b/Simulation/G4Extensions/Pythia8Decayer/src/components/Pythia8Decayer_entries.cxx
new file mode 100644
index 00000000..8bbd4660
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/src/components/Pythia8Decayer_entries.cxx
@@ -0,0 +1,4 @@
+#include "../Pythia8DecayerPhysicsTool.h"
+
+DECLARE_COMPONENT( Pythia8DecayerPhysicsTool )
+
diff --git a/Simulation/G4Faser/G4FaserServices/python/G4FaserServicesConfigNew.py b/Simulation/G4Faser/G4FaserServices/python/G4FaserServicesConfigNew.py
index 839a4fda..1b17627a 100644
--- a/Simulation/G4Faser/G4FaserServices/python/G4FaserServicesConfigNew.py
+++ b/Simulation/G4Faser/G4FaserServices/python/G4FaserServicesConfigNew.py
@@ -5,6 +5,8 @@ from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
DetectorGeometrySvc, G4AtlasSvc, G4GeometryNotifierSvc, PhysicsListSvc=CompFactory.getComps("DetectorGeometrySvc","G4AtlasSvc","G4GeometryNotifierSvc","PhysicsListSvc",)
from G4FaserTools.G4GeometryToolConfig import G4AtlasDetectorConstructionToolCfg
+from G4StepLimitation.G4StepLimitationConfigNew import G4StepLimitationToolCfg
+from Pythia8Decayer.Pythia8DecayerConfigNew import Pythia8DecayerPhysicsToolCfg
def DetectorGeometrySvcCfg(ConfigFlags, name="DetectorGeometrySvc", **kwargs):
result = ComponentAccumulator()
@@ -34,9 +36,9 @@ def G4GeometryNotifierSvcCfg(ConfigFlags, name="G4GeometryNotifierSvc", **kwargs
def PhysicsListSvcCfg(ConfigFlags, name="PhysicsListSvc", **kwargs):
result = ComponentAccumulator()
- G4StepLimitationTool = CompFactory.G4StepLimitationTool
- PhysOptionList = [G4StepLimitationTool("G4StepLimitationTool")]
+ PhysOptionList = [ result.popToolsAndMerge(G4StepLimitationToolCfg(ConfigFlags)) ]
#PhysOptionList += ConfigFlags.Sim.PhysicsOptions # FIXME Missing functionality
+ PhysOptionList += [ result.popToolsAndMerge(Pythia8DecayerPhysicsToolCfg(ConfigFlags)) ]
PhysDecaysList = []
kwargs.setdefault("PhysOption", PhysOptionList)
kwargs.setdefault("PhysicsDecay", PhysDecaysList)
--
GitLab