Skip to content
Snippets Groups Projects
Commit 1d2ac3bd authored by John Chapman's avatar John Chapman
Browse files

Add new ZeroLifetimePositioner Service

The `ZeroLifetimePositioner` works around the case where a neutral
particle oscillates into its anti-particle in one `GenVertex` then
immediately decays in a second `GenVertex` at the same position and
time. This is currently a problem seen when B0/B0bar (511/-511) particles
with pre-defined oscillations and decays are fed into Geant4, which cannot
currently handle the concept of a zero-lifetime particle.
As Geant4 does not implement any processes for these particles other than
the pre-defined decays currently and the particles are neutral, it is safe
to shift the point at which the particle oscillates into its anti-particle
back along the (straight-line) trajectory a bit. For simplicity halfway
between the production point and the oscillation/decay point has been used.
The decay position is left unchanged.

With this change Geant4 then successfully simulates the decay in the
required position.

After simulation is complete the `ZeroLifetimePositioner` can then be used
to move the oscillation vertex back to the decay vertex position.


Former-commit-id: 3abe0a79af1d016dc1776952305ff54786a9d756
parent 8f375840
No related branches found
No related tags found
No related merge requests found
......@@ -61,6 +61,11 @@ def getCrabKissingVertexPositioner(name="CrabKissingVertexPositioner", **kwargs)
def getGenEventValidityChecker(name="GenEventValidityChecker", **kwargs):
return CfgMgr.Simulation__GenEventValidityChecker(name, **kwargs)
def getZeroLifetimePositioner(name="ZeroLifetimePositioner", **kwargs):
kwargs.setdefault('ApplyPatch', True)
kwargs.setdefault('RemovePatch', True)
return CfgMgr.Simulation__ZeroLifetimePositioner(name, **kwargs)
def getGenEventVertexPositioner(name="GenEventVertexPositioner", **kwargs):
from G4AtlasApps.SimFlags import simFlags
readVtxPosFromFile = simFlags.VertexOverrideFile.statusOn or simFlags.VertexOverrideEventFile.statusOn
......
......@@ -8,6 +8,7 @@ addTool("BeamEffects.BeamEffectsConfig.getLongBeamspotVertexPositioner", "L
addTool("BeamEffects.BeamEffectsConfig.getCrabKissingVertexPositioner", "CrabKissingVertexPositioner")
## GenEvent Manipulators
addTool("BeamEffects.BeamEffectsConfig.getGenEventValidityChecker", "GenEventValidityChecker")
addTool("BeamEffects.BeamEffectsConfig.getZeroLifetimePositioner", "ZeroLifetimePositioner")
addTool("BeamEffects.BeamEffectsConfig.getGenEventVertexPositioner", "GenEventVertexPositioner")
addTool("BeamEffects.BeamEffectsConfig.getGenEventBeamEffectBooster", "GenEventBeamEffectBooster")
addTool("BeamEffects.BeamEffectsConfig.getGenEventRotator", "GenEventRotator")
......
/*
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
// class header include
#include "ZeroLifetimePositioner.h"
// HepMC includes
#include "HepMC/GenEvent.h"
// CLHEP includes
#include "CLHEP/Vector/LorentzVector.h"
#include <limits>
#include <algorithm>
/** Constructor **/
Simulation::ZeroLifetimePositioner::ZeroLifetimePositioner( const std::string& name,
ISvcLocator* pSvcLocator )
: base_class(name, pSvcLocator)
{
declareProperty("ApplyPatch", m_applyPatch, "");
declareProperty("RemovePatch", m_removePatch, "");
declareProperty("PDGCodesToCheck", m_pdgCodesToCheck, "");
}
/** Athena algtool's Hooks */
StatusCode Simulation::ZeroLifetimePositioner::initialize()
{
ATH_MSG_VERBOSE("Initializing ...");
for(auto& pdgcode : m_pdgCodesToCheck) {
pdgcode = std::abs(pdgcode);
ATH_MSG_DEBUG("Will look for verices where " << pdgcode <<
" oscillates into -" << pdgcode << " (and vice versa).");
}
return StatusCode::SUCCESS;
}
/** Athena algtool's Hooks */
StatusCode Simulation::ZeroLifetimePositioner::finalize()
{
ATH_MSG_VERBOSE("Finalizing ...");
return StatusCode::SUCCESS;
}
StatusCode Simulation::ZeroLifetimePositioner::applyWorkaround(HepMC::GenEvent& ge) const
{
ATH_MSG_DEBUG("applyWorkaround");
return this->manipulate(ge, m_applyPatch, false);
}
StatusCode Simulation::ZeroLifetimePositioner::removeWorkaround(HepMC::GenEvent& ge) const
{
ATH_MSG_DEBUG("removeWorkaround");
return this->manipulate(ge, false, m_removePatch);
}
/** modifies (displaces) the given GenEvent */
StatusCode Simulation::ZeroLifetimePositioner::manipulate(HepMC::GenEvent& ge, bool applyPatch, bool removePatch) const
{
// loop over the vertices in the event
HepMC::GenEvent::vertex_iterator vtxIt = ge.vertices_begin();
const HepMC::GenEvent::vertex_iterator vtxItEnd = ge.vertices_end();
const auto pdgCodesBegin = m_pdgCodesToCheck.begin();
const auto pdgCodesEnd = m_pdgCodesToCheck.end();
for (; vtxIt != vtxItEnd; ++vtxIt) {
// quick access:
HepMC::GenVertex *curVtx = (*vtxIt);
if (curVtx->particles_in_size()!=1 || curVtx->particles_out_size()!=1) { continue; }
const int pdgIn=(*(curVtx->particles_in_const_begin()))->pdg_id();
const int pdgOut=(*(curVtx->particles_out_const_begin()))->pdg_id();
if (pdgIn!=-pdgOut ||
std::find(pdgCodesBegin, pdgCodesEnd, std::abs(pdgIn))== pdgCodesEnd) {
continue;
}
HepMC::GenVertex* nextVtx = (*(curVtx->particles_out_const_begin()))->end_vertex();
if(!nextVtx) { continue; }
ATH_MSG_DEBUG("Found a vertex to correct with incoming PDG code = " << pdgIn);
ATH_MSG_VERBOSE("Next Vertex:");
if (ATH_UNLIKELY(this->msgLvl (MSG::VERBOSE))) {
nextVtx->print();
}
const HepMC::FourVector &nextVec = nextVtx->position();
const CLHEP::HepLorentzVector nextPos( nextVec.x(), nextVec.y(), nextVec.z(), nextVec.t() );
ATH_MSG_VERBOSE("Current Vertex:");
if (ATH_UNLIKELY(this->msgLvl (MSG::VERBOSE))) {
curVtx->print();
}
if (applyPatch) {
HepMC::GenVertex* prevVtx = (*(curVtx->particles_in_const_begin()))->production_vertex();
ATH_MSG_VERBOSE("Previous Vertex:");
if (ATH_UNLIKELY(this->msgLvl (MSG::VERBOSE))) {
prevVtx->print();
}
const HepMC::FourVector &prevVec = prevVtx->position();
const CLHEP::HepLorentzVector prevPos( prevVec.x(), prevVec.y(), prevVec.z(), prevVec.t() );
CLHEP::HepLorentzVector newPos = 0.5*(prevPos+nextPos);
curVtx->set_position(newPos);
ATH_MSG_DEBUG("Revised current Vertex");
if (ATH_UNLIKELY(this->msgLvl (MSG::VERBOSE))) {
curVtx->print();
}
}
if (removePatch) {
CLHEP::HepLorentzVector newPos = nextPos;
curVtx->set_position(newPos);
ATH_MSG_DEBUG("Revised current Vertex");
if (ATH_UNLIKELY(this->msgLvl (MSG::VERBOSE))) {
curVtx->print();
}
}
}
return StatusCode::SUCCESS;
}
/*
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// ZeroLifetimePositioner.h, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
#ifndef BEAMEFFECTS_ZEROLIFETIMEPOSITIONER_H
#define BEAMEFFECTS_ZEROLIFETIMEPOSITIONER_H 1
// FrameWork includes
#include "GaudiKernel/ToolHandle.h"
#include "AthenaBaseComps/AthService.h"
#include "HepMC_Interfaces/IZeroLifetimePatcher.h"
#include <vector>
namespace Simulation {
/** @class ZeroLifetimePositioner
This tool works around the issue of zero-lifetime particles in Geant4.
*/
class ZeroLifetimePositioner : public extends<AthService, Simulation::IZeroLifetimePatcher> {
public:
/** Constructor with parameters */
ZeroLifetimePositioner( const std::string& name, ISvcLocator* pSvcLocator );
/** Athena algtool's Hooks */
StatusCode initialize() override final;
StatusCode finalize() override final;
/** Applies the workaround for zero-lifetime particles to the GenEvent */
virtual StatusCode applyWorkaround(HepMC::GenEvent& ge) const override final;
/** Removes the workaround for zero-lifetime particles from the GenEvent */
virtual StatusCode removeWorkaround(HepMC::GenEvent& ge) const override final;
private:
StatusCode manipulate(HepMC::GenEvent& ge, bool applyPatch, bool removePatch) const;
bool m_applyPatch{false};
bool m_removePatch{false};
std::vector<int> m_pdgCodesToCheck{511};
};
}
#endif //> !BEAMEFFECTS_ZEROLIFETIMEPOSITIONER_H
#include "GaudiKernel/DeclareFactoryEntries.h"
#include "../GenEventValidityChecker.h"
#include "../ZeroLifetimePositioner.h"
#include "../GenEventVertexPositioner.h"
#include "../VertexBeamCondPositioner.h"
#include "../LongBeamspotVertexPositioner.h"
......@@ -8,6 +9,7 @@
#include "../GenEventBeamEffectBooster.h"
#include "../GenEventRotator.h"
#include "../BeamEffectsAlg.h"
DECLARE_NAMESPACE_SERVICE_FACTORY( Simulation , ZeroLifetimePositioner )
DECLARE_NAMESPACE_TOOL_FACTORY( Simulation , GenEventValidityChecker )
DECLARE_NAMESPACE_TOOL_FACTORY( Simulation , GenEventVertexPositioner )
DECLARE_NAMESPACE_TOOL_FACTORY( Simulation , VertexBeamCondPositioner )
......
/*
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
#ifndef HEPMC_INTERFACES_IZEROLIFETIMEPATCHER_H
#define HEPMC_INTERFACES_IZEROLIFETIMEPATCHER_H 1
// Gaudi
// framework includes
#include "GaudiKernel/IInterface.h"
#include "GaudiKernel/StatusCode.h"
namespace HepMC {
class GenEvent;
}
namespace Simulation {
/**
@class IZeroLifetimePatcher
Interface for a service which works around the issue of zero-lifetime particles in Geant4.
*/
class IZeroLifetimePatcher : virtual public IInterface {
public:
/** Virtual destructor */
virtual ~IZeroLifetimePatcher(){}
/// Creates the InterfaceID and interfaceID() method
DeclareInterfaceID(IZeroLifetimePatcher, 1, 0);
/** Applies the workaround for zero-lifetime particles to the GenEvent */
virtual StatusCode applyWorkaround(HepMC::GenEvent& ge) const = 0;
/** Removes the workaround for zero-lifetime particles from the GenEvent */
virtual StatusCode removeWorkaround(HepMC::GenEvent& ge) const = 0;
};
} // end of namespace
#endif // HEPMC_INTERFACES_IZEROLIFETIMEPATCHER_H
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment