Skip to content
Snippets Groups Projects
Commit abeeab77 authored by Graeme Stewart's avatar Graeme Stewart
Browse files

Simulation/ISF/ISF_HepMC/ISF_HepMC_Services deleted from 20.20

parent 93b98091
No related merge requests found
################################################################################
# Package: ISF_HepMC_Services
################################################################################
# Declare the package name:
atlas_subdir( ISF_HepMC_Services )
# Declare the package's dependencies:
atlas_depends_on_subdirs( PUBLIC
GaudiKernel
PRIVATE
Control/AthenaBaseComps
Control/StoreGate
DetectorDescription/AtlasDetDescr
Generators/GeneratorObjects
Simulation/Barcode/BarcodeInterfaces
Simulation/ISF/ISF_Core/ISF_Event
Simulation/ISF/ISF_Core/ISF_Interfaces
Simulation/ISF/ISF_HepMC/ISF_HepMC_Interfaces )
# External dependencies:
find_package( CLHEP )
find_package( HepMC )
# Component(s) in the package:
atlas_add_component( ISF_HepMC_Services
src/*.cxx
src/components/*.cxx
INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS}
LINK_LIBRARIES ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} GaudiKernel AthenaBaseComps StoreGateLib SGtests AtlasDetDescr GeneratorObjects ISF_Event ISF_Interfaces )
# Install files from the package:
atlas_install_headers( ISF_HepMC_Services )
atlas_install_python_modules( python/*.py )
package ISF_HepMC_Services
author <Andreas.Salzburger@cern.ch>
manager Andreas Salzburger <Andreas.Salzburger@cern.ch>
manager Elmar Ritsch <Elmar.Ritsch@cern.ch>
manager Wolfgang Lukas <Wolfgang.Lukas@cern.ch>
#################################################################
# public use statements
use AtlasPolicy AtlasPolicy-*
use GaudiInterface GaudiInterface-* External
#################################################################
# private use statements
private
use AthenaBaseComps AthenaBaseComps-* Control
use AtlasDetDescr AtlasDetDescr-* DetectorDescription
use AtlasHepMC AtlasHepMC-* External
use AtlasCLHEP AtlasCLHEP-* External
use BarcodeInterfaces BarcodeInterfaces-* Simulation/Barcode
use GeneratorObjects GeneratorObjects-* Generators
use ISF_Event ISF_Event-* Simulation/ISF/ISF_Core
use ISF_HepMC_Interfaces ISF_HepMC_Interfaces-* Simulation/ISF/ISF_HepMC
use ISF_Interfaces ISF_Interfaces-* Simulation/ISF/ISF_Core
use StoreGate StoreGate-* Control
public
library ISF_HepMC_Services *.cxx components/*.cxx
apply_pattern component_library
apply_pattern declare_python_modules files="*.py"
private
# use this to enable debugging for this package
#macro cppdebugflags '$(cppdebugflags_s)'
#macro_remove componentshr_linkopts "-Wl,-s"
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
"""
Tools configurations for ISF
KG Tan, 17/06/2012
"""
from AthenaCommon.CfgGetter import getPrivateTool,getPrivateToolClone,getPublicTool,getPublicToolClone,\
getService,getServiceClone,getAlgorithm,getAlgorithmClone
from AthenaCommon.Constants import * # FATAL,ERROR etc.
from AthenaCommon.SystemOfUnits import *
from AthenaCommon.DetFlags import DetFlags
from ISF_Config.ISF_jobProperties import ISF_Flags # IMPORTANT: Flags must be set before tools are retrieved
def getGenericTruthService(name="ISF_TruthService", **kwargs):
kwargs.setdefault('BarcodeSvc' , ISF_Flags.BarcodeService() )
kwargs.setdefault('McEventCollection' , 'TruthEvent')
kwargs.setdefault('SkipIfNoChildren' , True)
kwargs.setdefault('SkipIfNoParentBarcode' , True)
kwargs.setdefault('StoreExtraBarcodes' , False)
kwargs.setdefault('ForceEndVtxInRegions' , [ ] )
if 'longLived' in ISF_Flags.Simulator(): #FIXME this should be configured in a nicer way. ATLASSIM-526
kwargs.setdefault('QuasiStableParticlesIncluded', True)
from ISF_HepMC_Services.ISF_HepMC_ServicesConf import ISF__HepMC_TruthSvc
return ISF__HepMC_TruthSvc(name, **kwargs);
def getMC15TruthService(name="ISF_MC15TruthService", **kwargs):
# importing Reflex dictionary to access AtlasDetDescr::AtlasRegion enum
import ROOT, cppyy
cppyy.loadDictionary('AtlasDetDescrDict')
AtlasRegion = ROOT.AtlasDetDescr
kwargs.setdefault('BeamPipeTruthStrategies' , [ 'ISF_MCTruthStrategyGroupID_MC15' ] ) # this is used for beam pipe but not BeamPipeCentral which uses same as ID
kwargs.setdefault('IDTruthStrategies' , [ 'ISF_MCTruthStrategyGroupID_MC15', 'ISF_MCTruthStrategyGroupIDHadInt_MC15' ] )
kwargs.setdefault('CaloTruthStrategies' , [ 'ISF_MCTruthStrategyGroupCaloMuBrem', 'ISF_MCTruthStrategyGroupCaloDecay_MC15' ])
kwargs.setdefault('MSTruthStrategies' , [])
kwargs.setdefault('IgnoreUndefinedBarcodes' , False)
kwargs.setdefault('PassWholeVertices' , False) # new for MC15 - can write out partial vertices.
kwargs.setdefault('ForceEndVtxInRegions' , [ AtlasRegion.fAtlasID ] )
return getGenericTruthService(name, **kwargs);
def getMC15aTruthService(name="ISF_MC15aTruthService", **kwargs):
kwargs.setdefault('ForceEndVtxInRegions' , [ ] )
return getMC15TruthService(name, **kwargs);
def getMC15aPlusTruthService(name="ISF_MC15aPlusTruthService", **kwargs):
# importing Reflex dictionary to access AtlasDetDescr::AtlasRegion enum
import ROOT, cppyy
cppyy.loadDictionary('AtlasDetDescrDict')
AtlasRegion = ROOT.AtlasDetDescr
kwargs.setdefault('ForceEndVtxInRegions' , [ AtlasRegion.fAtlasID ] )
return getMC15TruthService(name, **kwargs);
def getMC12TruthService(name="ISF_MC12TruthService", **kwargs):
kwargs.setdefault('BeamPipeTruthStrategies' , [ 'ISF_MCTruthStrategyGroupID' ] ) # this is used for beam pipe but not BeamPipeCentral which uses same as ID
kwargs.setdefault('IDTruthStrategies' , [ 'ISF_MCTruthStrategyGroupID', 'ISF_MCTruthStrategyGroupIDHadInt' ] )
kwargs.setdefault('CaloTruthStrategies' , [ 'ISF_MCTruthStrategyGroupCaloMuBrem' ])
kwargs.setdefault('MSTruthStrategies' , [])
kwargs.setdefault('IgnoreUndefinedBarcodes' , False)
kwargs.setdefault('PassWholeVertices' , True)
return getGenericTruthService(name, **kwargs);
def getMC12LLPTruthService(name="ISF_MC12TruthLLPService", **kwargs):
kwargs.setdefault('BeamPipeTruthStrategies' , [ 'ISF_MCTruthStrategyGroupID', 'ISF_LLPTruthStrategy' ] ) # this is used for beam pipe but not BeamPipeCentral which uses same as ID
kwargs.setdefault('IDTruthStrategies' , [ 'ISF_MCTruthStrategyGroupID', 'ISF_MCTruthStrategyGroupIDHadInt', 'ISF_LLPTruthStrategy' ] )
kwargs.setdefault('CaloTruthStrategies' , [ 'ISF_MCTruthStrategyGroupCaloMuBrem', 'ISF_LLPTruthStrategy' ])
kwargs.setdefault('MSTruthStrategies' , ['ISF_LLPTruthStrategy'])
return getMC12TruthService(name, **kwargs);
def getMC12PlusTruthService(name="ISF_MC12PlusTruthService", **kwargs):
# importing Reflex dictionary to access AtlasDetDescr::AtlasRegion enum
import ROOT, cppyy
cppyy.loadDictionary('AtlasDetDescrDict')
AtlasRegion = ROOT.AtlasDetDescr
kwargs.setdefault('ForceEndVtxInRegions' , [ AtlasRegion.fAtlasID ] )
return getMC12TruthService(name, **kwargs);
def getValidationTruthService(name="ISF_ValidationTruthService", **kwargs):
kwargs.setdefault('BeamPipeTruthStrategies' , [])
kwargs.setdefault('IDTruthStrategies' , [ 'ISF_ValidationTruthStrategy' ] )
kwargs.setdefault('CaloTruthStrategies' , [ 'ISF_ValidationTruthStrategy' ] )
kwargs.setdefault('MSTruthStrategies' , [])
kwargs.setdefault('IgnoreUndefinedBarcodes' , True)
kwargs.setdefault('PassWholeVertices' , True)
return getGenericTruthService(name, **kwargs);
def getTruthService(name="ISF_TruthService", **kwargs):
if ISF_Flags.ValidationMode() :
return getValidationTruthService( name, **kwargs )
else :
return getMC12TruthService( name, **kwargs );
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
from AthenaCommon.CfgGetter import addTool, addService, addAlgorithm
addService("ISF_HepMC_Services.ISF_HepMC_ServicesConfig.getTruthService" , "ISF_TruthService")
addService("ISF_HepMC_Services.ISF_HepMC_ServicesConfig.getMC12TruthService" , "ISF_MC12TruthService")
addService("ISF_HepMC_Services.ISF_HepMC_ServicesConfig.getMC12PlusTruthService" , "ISF_MC12PlusTruthService")
addService("ISF_HepMC_Services.ISF_HepMC_ServicesConfig.getMC12LLPTruthService" , "ISF_MC12LLPTruthService")
addService("ISF_HepMC_Services.ISF_HepMC_ServicesConfig.getMC15aTruthService" , "ISF_MC15aTruthService")
addService("ISF_HepMC_Services.ISF_HepMC_ServicesConfig.getMC15aPlusTruthService" , "ISF_MC15aPlusTruthService")
addService("ISF_HepMC_Services.ISF_HepMC_ServicesConfig.getMC15TruthService" , "ISF_MC15TruthService")
addService("ISF_HepMC_Services.ISF_HepMC_ServicesConfig.getValidationTruthService" , "ISF_ValidationTruthService")
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// HepMC_TruthSvc.cxx, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
// class header
#include "HepMC_TruthSvc.h"
// other ISF_HepMC includes
#include "ISF_HepMC_Interfaces/ITruthStrategy.h"
// ISF includes
#include "ISF_Event/ITruthIncident.h"
// Framework
#include "GaudiKernel/ISvcLocator.h"
#include "StoreGate/StoreGateSvc.h"
#include "GaudiKernel/SystemOfUnits.h"
// Barcode
#include "BarcodeInterfaces/IBarcodeSvc.h"
// HepMC includes
#include "HepMC/SimpleVector.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"
// CLHEP includes
#include "CLHEP/Geometry/Point3D.h"
// DetectorDescription
#include "AtlasDetDescr/AtlasRegionHelper.h"
/** Constructor **/
ISF::HepMC_TruthSvc::HepMC_TruthSvc(const std::string& name,ISvcLocator* svc) :
AthService(name,svc),
m_storeGate(0),
m_barcodeSvc("BarcodeSvc",name),
m_barcodeSvcQuick(0),
m_collectionName("TruthEvent"),
m_mcEventCollection(0),
m_mcEvent(0),
m_geoStrategies(),
m_numStrategies(),
m_skipIfNoChildren(true),
m_skipIfNoParentBarcode(true),
m_ignoreUndefinedBarcodes(false),
m_screenOutputPrefix("isf >> "),
m_screenEmptyPrefix(),
m_storeExtraBCs(true),
m_passWholeVertex(true),
m_forceEndVtxRegionsVec(),
m_forceEndVtx(),
m_quasiStableParticlesIncluded(false),
m_secondaryParticleBcOffset(Barcode::fUndefinedBarcode),
m_myLowestVertexBC(Barcode::fUndefinedBarcode)
{
// the particle stack filler tool
declareProperty("McEventCollection", m_collectionName );
// refine the screen output for debugging
declareProperty("ScreenOutputPrefix", m_screenOutputPrefix );
// the barcode service (used to compute Vertex Barco des)
declareProperty("BarcodeSvc", m_barcodeSvc );
// MCTruth writing strategies
declareProperty("SkipIfNoChildren", m_skipIfNoChildren );
declareProperty("SkipIfNoParentBarcode", m_skipIfNoParentBarcode );
declareProperty("IgnoreUndefinedBarcodes", m_ignoreUndefinedBarcodes );
declareProperty("PassWholeVertices", m_passWholeVertex );
// the truth strategies for the different AtlasDetDescr regions
declareProperty("BeamPipeTruthStrategies", m_geoStrategyHandles[AtlasDetDescr::fAtlasForward] );
declareProperty("IDTruthStrategies", m_geoStrategyHandles[AtlasDetDescr::fAtlasID] );
declareProperty("CaloTruthStrategies", m_geoStrategyHandles[AtlasDetDescr::fAtlasCalo] );
declareProperty("MSTruthStrategies", m_geoStrategyHandles[AtlasDetDescr::fAtlasMS] );
declareProperty("CavernTruthStrategies", m_geoStrategyHandles[AtlasDetDescr::fAtlasCavern] );
// attach end-vertex if parent particle dies for the different AtlasDetDescr regions
declareProperty("ForceEndVtxInRegions", m_forceEndVtxRegionsVec );
declareProperty("StoreExtraBarcodes", m_storeExtraBCs);
declareProperty("QuasiStableParticlesIncluded", m_quasiStableParticlesIncluded);
}
ISF::HepMC_TruthSvc::~HepMC_TruthSvc()
{}
/** Query the interfaces. */
StatusCode ISF::HepMC_TruthSvc::queryInterface(const InterfaceID& riid, void** ppvInterface)
{
if ( IID_ITruthSvc == riid )
*ppvInterface = (ITruthSvc*)this;
else {
// Interface is not directly available: try out a base class
return Service::queryInterface(riid, ppvInterface);
}
addRef();
return StatusCode::SUCCESS;
}
/** framework methods */
StatusCode ISF::HepMC_TruthSvc::initialize()
{
ATH_MSG_VERBOSE( "initialize()" );
// Screen output
for (size_t prl = 0; prl < m_screenOutputPrefix.size(); ++prl) m_screenEmptyPrefix += " ";
ATH_MSG_DEBUG ( m_screenOutputPrefix << "--------------------------------------------------------" );
// retrieve StoreGateSvc
ISvcLocator* svcLocator = Gaudi::svcLocator(); // from Bootstrap
if ( svcLocator->service("StoreGateSvc", m_storeGate).isFailure()) {
ATH_MSG_WARNING("AthenaHitsCollectionHelper: could not accessStoreGateSvc!");
return StatusCode::FAILURE;
}
// retrieve BarcodeSvc
if ( m_barcodeSvc.retrieve().isFailure() ) {
ATH_MSG_FATAL( m_screenOutputPrefix << "Could not retrieve BarcodeService. Abort.");
return StatusCode::FAILURE;
}
// store a point directly to the BarcodeSvc class
// (removes the gaudi overhead in each call)
m_barcodeSvcQuick = &(*m_barcodeSvc);
// retrieve the strategies for each atlas region (Athena Alg Tools)
// and setup whether we want to write end-vertices in this region whenever a truth particle dies
for ( unsigned short geoID=AtlasDetDescr::fFirstAtlasRegion; geoID<AtlasDetDescr::fNumAtlasRegions; ++geoID) {
if ( m_geoStrategyHandles[geoID].retrieve().isFailure() ) {
ATH_MSG_FATAL( m_screenOutputPrefix << "Could not retrieve TruthStrategy Tool Array for SimGeoID="
<< AtlasDetDescr::AtlasRegionHelper::getName(geoID) << ". Abort.");
return StatusCode::FAILURE;
}
// copy a pointer to the strategy instance to the local
// array of pointers (for faster access)
unsigned short curNumStrategies = m_geoStrategyHandles[geoID].size();
m_numStrategies[geoID] = curNumStrategies;
m_geoStrategies[geoID] = new ISF::ITruthStrategy*[curNumStrategies];
for ( unsigned short i = 0; i < curNumStrategies; ++i) {
m_geoStrategies[geoID][i] = &(*m_geoStrategyHandles[geoID][i]);
}
// create an end-vertex for all truth particles ending in the current AtlasRegion?
bool forceEndVtx = std::find( m_forceEndVtxRegionsVec.begin(),
m_forceEndVtxRegionsVec.end(),
geoID ) != m_forceEndVtxRegionsVec.end();
m_forceEndVtx[geoID] = forceEndVtx;
}
ATH_MSG_VERBOSE( "initialize() successful" );
return StatusCode::SUCCESS;
}
/** framework methods */
StatusCode ISF::HepMC_TruthSvc::finalize()
{
ATH_MSG_VERBOSE ( m_screenOutputPrefix << "Finalizing ..." );
return StatusCode::SUCCESS;
}
/** Initialize the HepMC_TruthSvc and the truthSvc */
StatusCode ISF::HepMC_TruthSvc::initializeTruthCollection()
{
m_mcEventCollection = 0;
// retrieve McEventCollection from storegate
if ( m_collectionName.empty()) {
StatusCode status = m_storeGate->retrieve( m_mcEventCollection);
if (status.isFailure())
ATH_MSG_WARNING( "Unable to retrieve McEventCollection with name=" << m_collectionName);
else
ATH_MSG_DEBUG( "Sucessfully retrieved McEventCollection with name=" << m_collectionName);
} else {
if (m_storeGate->contains<McEventCollection>( m_collectionName)) {
StatusCode status = m_storeGate->retrieve( m_mcEventCollection, m_collectionName);
if (status.isFailure())
ATH_MSG_WARNING( "Unable to retrieve McEventCollection with name=" << m_collectionName
<< ". Will create new collection.");
else
ATH_MSG_DEBUG( "Sucessfully retrieved McEventCollection with name=" << m_collectionName);
}
}
// in case no McEventCollection is present in storegate -> create new one and record it on SG
if (! m_mcEventCollection) {
McEventCollection* newcoll = new McEventCollection();
HepMC::GenEvent *evt = new HepMC::GenEvent();
newcoll->push_back(evt);
StatusCode status = m_storeGate->record(newcoll, m_collectionName, true);
if (status.isFailure()) {
ATH_MSG_ERROR( "Unable to record newly created McEventCollection with name=" << m_collectionName);
return StatusCode::FAILURE;
}
status = m_storeGate->retrieve( m_mcEventCollection, m_collectionName);
if (status.isFailure()) {
ATH_MSG_ERROR( "Unable to retrieve newly created McEventCollection with name=" << m_collectionName);
return StatusCode::FAILURE;
}
}
// collect last GenEvent from McEventCollection
m_mcEvent = m_mcEventCollection->back();
if (!m_mcEvent) {
// no GenEvent present
// -> create new (emtpy) GenEvent and add it to McEventCollection
m_mcEvent = new HepMC::GenEvent();
m_mcEventCollection->push_back( m_mcEvent);
}
// retrieve secondary particle barcode and vertex barcode offsets from the BarcodeService
m_myLowestVertexBC = m_barcodeSvcQuick->secondaryVertexBcOffset();
m_secondaryParticleBcOffset = m_barcodeSvcQuick->secondaryParticleBcOffset();
return StatusCode::SUCCESS;
}
StatusCode ISF::HepMC_TruthSvc::releaseEvent() {
// set the output collection to const
return m_storeGate->setConst(m_mcEventCollection);
}
/** Register a truth incident */
void ISF::HepMC_TruthSvc::registerTruthIncident( ISF::ITruthIncident& ti) {
// pass whole vertex or individual child particles
ti.setPassWholeVertices(m_passWholeVertex);
// the GeoID
AtlasDetDescr::AtlasRegion geoID = ti.geoID();
// check geoID assigned to the TruthIncident
if ( !validAtlasRegion(geoID) ) {
ATH_MSG_ERROR("Unable to register truth incident with unknown SimGeoID="<< geoID);
return;
}
ATH_MSG_VERBOSE( "Registering TruthIncident for SimGeoID="
<< AtlasDetDescr::AtlasRegionHelper::getName(geoID) );
// number of child particles
unsigned short numSec = ti.numberOfChildren();
if ( m_skipIfNoChildren && (numSec==0) ) {
ATH_MSG_VERBOSE( "No child particles present in the TruthIncident,"
<< " will not record this TruthIncident.");
return;
}
// the parent particle -> get its barcode
Barcode::ParticleBarcode parentBC = ti.parentBarcode();
if ( m_skipIfNoParentBarcode && (parentBC==Barcode::fUndefinedBarcode) ) {
ATH_MSG_VERBOSE( "Parent particle in TruthIncident does not have a barcode,"
<< " will not record this TruthIncident.");
return;
}
// loop over registered truth strategies for given geoID
bool pass = false;
for ( unsigned short stratID=0; (!pass) && (stratID<m_numStrategies[geoID]); stratID++) {
// (*) test if given TruthIncident passes current strategy
pass = m_geoStrategies[geoID][stratID]->pass(ti);
}
if (pass) {
// at least one truth stategy returend true
// -> record incident
recordIncidentToMCTruth( ti);
} else {
// none of the truth strategies returned true
// -> child particles will NOT be added to the TruthEvent collection
// attach parent particle end vertex if it gets killed by this interaction
if ( m_forceEndVtx[geoID] && !ti.parentSurvivesIncident() ) {
HepMC::GenVertex *vtx = createGenVertexFromTruthIncident( ti);
m_mcEvent->add_vertex( vtx);
}
// -> assign shared barcode to all child particles (if barcode service supports it)
setSharedChildParticleBarcode( ti);
}
return;
}
/** Record the given truth incident to the MC Truth */
void ISF::HepMC_TruthSvc::recordIncidentToMCTruth( ISF::ITruthIncident& ti) {
Barcode::PhysicsProcessCode processCode = ti.physicsProcessCode();
Barcode::ParticleBarcode parentBC = ti.parentBarcode();
// record the GenVertex
HepMC::GenVertex *vtx = createGenVertexFromTruthIncident(ti);
ATH_MSG_VERBOSE ( "Outgoing particles:" );
// update parent barcode and add it to the vertex as outgoing particle
Barcode::ParticleBarcode newPrimBC = m_barcodeSvcQuick->incrementBarcode( parentBC, processCode);
if ( newPrimBC == Barcode::fUndefinedBarcode) {
if (m_ignoreUndefinedBarcodes) {
ATH_MSG_WARNING("Unable to generate new Particle Barcode. Continuing due to 'IgnoreUndefinedBarcodes'==True");
} else {
ATH_MSG_ERROR("Unable to generate new Particle Barcode. Aborting");
abort();
}
}
HepMC::GenParticle *parentAfterIncident = ti.parentParticleAfterIncident( newPrimBC );
if(parentAfterIncident) {
ATH_MSG_VERBOSE ( "Parent After Incident: " << *parentAfterIncident);
vtx->add_particle_out( parentAfterIncident );
}
// update all _extra_ barcodes of child particles with parent info
// MB: sofar extra barcode only contains parent info, so can be the same for each child
if (m_storeExtraBCs) {
// the parent particle -> get its extra barcode
Barcode::ParticleBarcode parentExtraBC = ti.parentExtraBarcode();
ti.setAllChildrenExtraBarcodes( parentExtraBC );
}
// add child particles to the vertex
unsigned short numSec = ti.numberOfChildren();
for ( unsigned short i=0; i<numSec; ++i) {
bool writeOutChild = m_passWholeVertex || ti.childPassedFilters(i);
if (writeOutChild) {
// generate a new barcode for the child particle
Barcode::ParticleBarcode secBC = m_barcodeSvcQuick->newSecondary( parentBC, processCode);
if ( secBC == Barcode::fUndefinedBarcode) {
if (m_ignoreUndefinedBarcodes)
ATH_MSG_WARNING("Unable to generate new Secondary Particle Barcode. Continuing due to 'IgnoreUndefinedBarcodes'==True");
else {
ATH_MSG_ERROR("Unable to generate new Secondary Particle Barcode. Aborting");
abort();
}
}
HepMC::GenParticle *p = ti.childParticle(i, secBC );
ATH_MSG_VERBOSE ( "Writing out " << i << "th child particle: " << *p);
// add particle to vertex
vtx->add_particle_out( p);
// Check to see if this is meant to be a "parent" vertex
if ( p && p->barcode() < m_secondaryParticleBcOffset ){
vtx->suggest_barcode( m_myLowestVertexBC );
++m_myLowestVertexBC;
if (m_quasiStableParticlesIncluded){
ATH_MSG_VERBOSE( "Found a case of low barcode (" << p->barcode() << " < " <<
m_secondaryParticleBcOffset << " changing vtx barcode to " << m_myLowestVertexBC+1 );
} else {
ATH_MSG_WARNING( "Modifying vertex barcode, but no apparent quasi-stable particle simulation enabled." );
ATH_MSG_WARNING( "This means that you have encountered a very strange configuration. Watch out!" );
}
}
} // <-- if write out child particle
else {
ATH_MSG_VERBOSE ( "Not writing out " << i << "th child particle." );
}
} // <-- loop over all child particles
// finally add the vertex to the current GenEvent
m_mcEvent->add_vertex( vtx);
}
/** Record the given truth incident to the MC Truth */
HepMC::GenVertex *ISF::HepMC_TruthSvc::createGenVertexFromTruthIncident( ISF::ITruthIncident& ti) {
Barcode::PhysicsProcessCode processCode = ti.physicsProcessCode();
Barcode::ParticleBarcode parentBC = ti.parentBarcode();
std::vector<double> weights(1);
Barcode::ParticleBarcode primaryBC = parentBC % m_barcodeSvcQuick->particleGenerationIncrement();
weights[0] = static_cast<double>( primaryBC );
// Check for a previous end vertex on this particle. If one existed, then we should put down next to this
// a new copy of the particle. This is the agreed upon version of the quasi-stable particle truth, where
// the vertex at which we start Q-S simulation no longer conserves energy, but we keep both copies of the
// truth particles
HepMC::GenParticle *parent = ti.parentParticle();
if (!parent) {
ATH_MSG_ERROR("Unable to write particle interaction to MC truth due to missing parent HepMC::GenParticle instance");
abort();
}
// generate vertex
Barcode::VertexBarcode vtxbcode = m_barcodeSvcQuick->newVertex( parentBC, processCode );
if ( vtxbcode == Barcode::fUndefinedBarcode) {
if (m_ignoreUndefinedBarcodes) {
ATH_MSG_WARNING("Unable to generate new Truth Vertex Barcode. Continuing due to 'IgnoreUndefinedBarcodes'==True");
} else {
ATH_MSG_ERROR("Unable to generate new Truth Vertex Barcode. Aborting");
abort();
}
}
int vtxID = 1000 + static_cast<int>(processCode);
HepMC::GenVertex *vtx = new HepMC::GenVertex( ti.position(), vtxID, weights );
vtx->suggest_barcode( vtxbcode );
if (parent->end_vertex()){
if(!m_quasiStableParticlesIncluded) {
ATH_MSG_WARNING("Parent particle found with an end vertex attached. This should only happen");
ATH_MSG_WARNING("in the case of simulating quasi-stable particles. That functionality");
ATH_MSG_WARNING("is not yet validated in ISF, so you'd better know what you're doing.");
ATH_MSG_WARNING("Will delete the old vertex and swap in the new one.");
}
// Remove the old vertex from the event
parent->parent_event()->remove_vertex( parent->end_vertex() );
// Now add the new vertex to the new parent
vtx->add_particle_in( parent );
ATH_MSG_VERBOSE ( "QS End Vertex representing process: " << processCode << ", for parent with barcode "<<parentBC<<". Creating." );
ATH_MSG_VERBOSE ( "Parent: " << *parent);
} else { // Normal simulation
// add parent particle to vtx
vtx->add_particle_in( parent );
ATH_MSG_VERBOSE ( "End Vertex representing process: " << processCode << ", for parent with barcode "<<parentBC<<". Creating." );
ATH_MSG_VERBOSE ( "Parent: " << *parent);
}
return vtx;
}
/** Set shared barcode for child particles particles */
void ISF::HepMC_TruthSvc::setSharedChildParticleBarcode( ISF::ITruthIncident& ti) {
Barcode::PhysicsProcessCode processCode = ti.physicsProcessCode();
Barcode::ParticleBarcode parentBC = ti.parentBarcode();
ATH_MSG_VERBOSE ( "End Vertex representing process: " << processCode << ". TruthIncident failed cuts. Skipping.");
// generate one new barcode for all child particles
Barcode::ParticleBarcode childBC = m_barcodeSvcQuick->sharedChildBarcode( parentBC, processCode);
// propagate this barcode into the TruthIncident only if
// it is a proper barcode, ie !=fUndefinedBarcode
if (childBC != Barcode::fUndefinedBarcode) {
ti.setAllChildrenBarcodes( childBC );
}
// update all _extra_ barcodes of child particles with parent info
// MB: sofar extra barcode only contains parent info, so can be the same for each child particle
if (m_storeExtraBCs) {
// the parent particle -> get its extra barcode
Barcode::ParticleBarcode parentExtraBC = ti.parentExtraBarcode();
ti.setAllChildrenExtraBarcodes( parentExtraBC );
}
}
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// HepMC_TruthSvc.h, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
#ifndef ISF_SERVICES_HEPMC_TRUTHSVC_H
#define ISF_SERVICES_HEPMC_TRUTHSVC_H 1
// STL includes
#include <string>
// FrameWork includes
#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/ServiceHandle.h"
#include "AthenaBaseComps/AthService.h"
// ISF include
#include "ISF_Interfaces/ITruthSvc.h"
// DetectorDescription
#include "AtlasDetDescr/AtlasRegion.h"
// Barcode
#include "BarcodeInterfaces/Barcode.h"
// McEventCollection
#include "GeneratorObjects/McEventCollection.h"
// forward declarations
class StoreGateSvc;
namespace Barcode {
class IBarcodeSvc;
}
namespace HepMC {
class GenEvent;
}
namespace ISF {
class ITruthStrategy;
typedef ToolHandleArray<ITruthStrategy> TruthStrategyArray;
/** @class HepMC_TruthSvc
HepMC based version of the ISF::ITruthSvc,
currently it takes an ITruthIncident base class
@author Andreas.Salzburger -at- cern.ch , Elmar.Ritsch -at- cern.ch
*/
class HepMC_TruthSvc : public AthService, public ITruthSvc {
public:
//** Constructor with parameters */
HepMC_TruthSvc( const std::string& name, ISvcLocator* pSvcLocator );
/** Destructor */
virtual ~HepMC_TruthSvc();
/** Athena algorithm's interface method initialize() */
StatusCode initialize() override final;
/** Athena algorithm's interface method finalize() */
StatusCode finalize() override final;
/** Register a truth incident */
void registerTruthIncident( ITruthIncident& truthincident) override final;
/** Initialize the Truth Svc at the beginning of each event */
StatusCode initializeTruthCollection() override final;
/** Finalize the Truth Svc at the end of each event*/
StatusCode releaseEvent() override final;
/** Query the interfaces. **/
StatusCode queryInterface( const InterfaceID& riid, void** ppvInterface );
private:
/** Record the given truth incident to the MC Truth */
void recordIncidentToMCTruth( ITruthIncident& truthincident);
/** Record and end vertex to the MC Truth for the parent particle */
HepMC::GenVertex *createGenVertexFromTruthIncident( ITruthIncident& truthincident);
/** Set shared barcode for child particles */
void setSharedChildParticleBarcode( ITruthIncident& truthincident);
StoreGateSvc *m_storeGate; //!< The storegate svc
ServiceHandle<Barcode::IBarcodeSvc> m_barcodeSvc; //!< The Barcode service
Barcode::IBarcodeSvc *m_barcodeSvcQuick; //!< The Barcode service for quick access
std::string m_collectionName; //!< name of the output McEventCollection
McEventCollection *m_mcEventCollection; //!< pointer to the McEventCollection
HepMC::GenEvent *m_mcEvent;
/** the truth strategie applied (as AthenaToolHandle Array) */
TruthStrategyArray m_geoStrategyHandles[AtlasDetDescr::fNumAtlasRegions];
/** for faster access: using an internal pointer to the actual ITruthStrategy instances */
ITruthStrategy** m_geoStrategies[AtlasDetDescr::fNumAtlasRegions];
unsigned short m_numStrategies[AtlasDetDescr::fNumAtlasRegions];
/** MCTruth steering */
bool m_skipIfNoChildren; //!< do not record incident if numChildren==0
bool m_skipIfNoParentBarcode; //!< do not record if parentBarcode==fUndefinedBarcode
bool m_ignoreUndefinedBarcodes;//!< do/don't abort if retrieve an undefined barcode
std::string m_screenOutputPrefix;
std::string m_screenEmptyPrefix;
bool m_storeExtraBCs;
bool m_passWholeVertex;
std::vector<bool> m_forceEndVtxRegionsVec; //!< property containing AtlasRegions for which
// to write end-vtx
bool m_forceEndVtx[AtlasDetDescr::fNumAtlasRegions]; //!< attach end vertex to
// all parent particles if they die
bool m_quasiStableParticlesIncluded; //!< does this job simulate quasi-stable particles.
Barcode::ParticleBarcode m_secondaryParticleBcOffset;
Barcode::VertexBarcode m_myLowestVertexBC;
};
}
#endif //> !ISF_SERVICES_HEPMC_TRUTHSVC_H
#include "GaudiKernel/DeclareFactoryEntries.h"
#include "../HepMC_TruthSvc.h"
DECLARE_NAMESPACE_SERVICE_FACTORY( ISF , HepMC_TruthSvc )
DECLARE_FACTORY_ENTRIES( ISF_HepMC_Services ) {
DECLARE_NAMESPACE_SERVICE( ISF , HepMC_TruthSvc )
}
#include "GaudiKernel/LoadFactoryEntries.h"
LOAD_FACTORY_ENTRIES( ISF_HepMC_Services )
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