diff --git a/Simulation/G4Atlas/G4AtlasApps/python/SimFlags.py b/Simulation/G4Atlas/G4AtlasApps/python/SimFlags.py index f5f64f01e046d222381ddb6b81206ed09724dd29..46d7a160d2e51f5a301c99670ab7825a692e186c 100644 --- a/Simulation/G4Atlas/G4AtlasApps/python/SimFlags.py +++ b/Simulation/G4Atlas/G4AtlasApps/python/SimFlags.py @@ -639,6 +639,14 @@ class ParticleID(JobProperty): allowedTypes = ['bool'] StoredValue = False +class RecordStepInfo(JobProperty): + """ + Should FCS_StepInfoCollections be recorded + """ + statusOn = True + allowedTypes = ['bool'] + StoredValue = False + class RecordFlux(JobProperty): """ Record flux through the entirety of the detector diff --git a/Simulation/G4Atlas/G4AtlasTools/python/G4AtlasToolsConfig.py b/Simulation/G4Atlas/G4AtlasTools/python/G4AtlasToolsConfig.py index 5d5439468b791ec8d18afbff8169994db23c84c4..337d297ff1929c3332063162c132e453d84b4f2f 100644 --- a/Simulation/G4Atlas/G4AtlasTools/python/G4AtlasToolsConfig.py +++ b/Simulation/G4Atlas/G4AtlasTools/python/G4AtlasToolsConfig.py @@ -112,6 +112,9 @@ def generateCaloSensitiveDetectorList(): SensitiveDetectorList += [ 'TileGeoG4CalibSD' ] # mode 1 : With CaloCalibrationHits else: SensitiveDetectorList += [ 'TileGeoG4SD' ] # mode 0 : No CaloCalibrationHits + from G4AtlasApps.SimFlags import simFlags + if simFlags.RecordStepInfo.get_Value(): + SensitiveDetectorList += [ 'FCS_StepInfoSensitiveDetector' ] return SensitiveDetectorList def generateMuonSensitiveDetectorList(): diff --git a/Simulation/ISF/ISF_Example/python/ISF_Output.py b/Simulation/ISF/ISF_Example/python/ISF_Output.py index fcb2bc118614cd3047c515a197939ddf9c07a462..74c397b07b39583c4d4f13a1c48cbfa9c0b9b855 100644 --- a/Simulation/ISF/ISF_Example/python/ISF_Output.py +++ b/Simulation/ISF/ISF_Example/python/ISF_Output.py @@ -70,7 +70,7 @@ def getHITSStreamItemList(): ## TimingAlg hitsItemList +=["RecoTimingObj#EVNTtoHITS_timings"] - if 'G4UA::FastCaloSimParamActionTool' in simFlags.OptionalUserActionList.get_Value()['Step']: + if simFlags.RecordStepInfo.get_Value(): hitsItemList +=["ISF_FCS_Parametrization::FCS_StepInfoCollection#MergedEventSteps"] ## add xAOD::TrackParticles output collection Parametric Simulation diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPostInclude.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPostInclude.py index 67017a344bb85806b0cb4485b6529dbf93939eee..a59ee537a7a600deaaf28d45472527c2b91fb6cd 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPostInclude.py +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPostInclude.py @@ -1,3 +1,23 @@ from ISF_FastCaloSimParametrization.ISF_FastCaloSimParametrizationConf import FastCaloSimParamAlg topSeq += FastCaloSimParamAlg() ISF_HITSStream.stream1.ItemList += ["ISF_FCS_Parametrization::FCS_StepInfoCollection#MergedEventSteps"] +from AthenaCommon.CfgGetter import getPublicTool +stepInfoSDTool = getPublicTool("SensitiveDetectorMasterTool").SensitiveDetectors['FCS_StepInfoSensitiveDetector'] +stepInfoSDTool.shift_lar_subhit=True +stepInfoSDTool.shorten_lar_step=True + +stepInfoSDTool.maxEtaPS=1 +stepInfoSDTool.maxPhiPS=5 +stepInfoSDTool.maxrPS=0 + +stepInfoSDTool.maxEtaEM1=1 +stepInfoSDTool.maxPhiEM1=5 +stepInfoSDTool.maxrEM1=15 + +stepInfoSDTool.maxEtaEM2=1 +stepInfoSDTool.maxPhiEM2=5 +stepInfoSDTool.maxrEM2=60 + +stepInfoSDTool.maxEtaEM3=1 +stepInfoSDTool.maxPhiEM3=5 +stepInfoSDTool.maxrEM3=8 diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPostInclude_1mm.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPostInclude_1mm.py new file mode 100644 index 0000000000000000000000000000000000000000..b42cd82a47bed2666c9ae87fcf77a95df3a2c4e1 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPostInclude_1mm.py @@ -0,0 +1,8 @@ +from ISF_FastCaloSimParametrization.ISF_FastCaloSimParametrizationConf import FastCaloSimParamAlg +topSeq += FastCaloSimParamAlg() +ISF_HITSStream.stream1.ItemList += ["ISF_FCS_Parametrization::FCS_StepInfoCollection#MergedEventSteps"] +from AthenaCommon.CfgGetter import getPublicTool +stepInfoSDTool = getPublicTool("SensitiveDetectorMasterTool").SensitiveDetectors['FCS_StepInfoSensitiveDetector'] +stepInfoSDTool.shift_lar_subhit=True +stepInfoSDTool.shorten_lar_step=True +stepInfoSDTool.maxRadiusLAr=1 diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPreInclude.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPreInclude.py index b10f924dd25a540cad352e7af86af1869323773b..ed3a162c17b790a20a87e26b63ba21df04bb1f51 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPreInclude.py +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPreInclude.py @@ -1,20 +1,2 @@ from G4AtlasApps.SimFlags import simFlags -simFlags.OptionalUserActionList.addAction('G4UA::FastCaloSimParamActionTool',['BeginOfEvent','EndOfEvent','BeginOfRun','EndOfRun','Step']) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"shift_lar_subhit",True) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"shorten_lar_step",True) -# optimised merging scheme parameters -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxEtaPS",1) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxPhiPS",5) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxrPS",0) - -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxEtaEM1",1) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxPhiEM1",5) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxrEM1",15) - -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxEtaEM2",1) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxPhiEM2",5) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxrEM2",60) - -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxEtaEM3",1) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxPhiEM3",5) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxrEM3",8) +simFlags.RecordStepInfo=True diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPreInclude_1mm.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPreInclude_1mm.py deleted file mode 100644 index d4c0a586f46d8c4cba3ff89c2e1e4c717bcc832c..0000000000000000000000000000000000000000 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/share/ISF_FastCaloSimParametrization_SimPreInclude_1mm.py +++ /dev/null @@ -1,6 +0,0 @@ -from G4AtlasApps.SimFlags import simFlags -simFlags.OptionalUserActionList.addAction('G4UA::FastCaloSimParamActionTool',['BeginOfEvent','EndOfEvent','BeginOfRun','EndOfRun','Step']) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"shift_lar_subhit",True) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"shorten_lar_step",True) -simFlags.UserActionConfig.addConfig('G4UA::FastCaloSimParamActionTool',"maxRadiusLAr",1) -# optimised merging scheme parameters are not required for 1mm merging diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/CMakeLists.txt b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..4c03c106111bed7f6a8b90f37a0839758b36b8ca --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/CMakeLists.txt @@ -0,0 +1,41 @@ +################################################################################ +# Package: ISF_FastCaloSimSD +################################################################################ + +# Declare the package name: +atlas_subdir( ISF_FastCaloSimSD ) + +# Declare the package's dependencies: +atlas_depends_on_subdirs( PUBLIC + Control/StoreGate + GaudiKernel + Simulation/G4Atlas/G4AtlasInterfaces + Simulation/G4Atlas/G4AtlasTools + TileCalorimeter/TileG4/TileG4Interfaces + PRIVATE + Calorimeter/CaloDetDescr + Calorimeter/CaloIdentifier + Control/CxxUtils + Control/AthenaBaseComps + Generators/GeneratorObjects + LArCalorimeter/LArG4/LArG4Code + Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent + DetectorDescription/Identifier ) + +# External dependencies: +find_package( CLHEP ) +find_package( Geant4 ) +find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread ) +find_package( XercesC ) + +# Component(s) in the package: +atlas_add_component( ISF_FastCaloSimSD + src/*.cxx + src/components/*.cxx + INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${XERCESC_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} + LINK_LIBRARIES ${ROOT_LIBRARIES} ${XERCESC_LIBRARIES} ${GEANT4_LIBRARIES} ${CLHEP_LIBRARIES} StoreGateLib SGtests GaudiKernel G4AtlasInterfaces G4AtlasToolsLib CaloDetDescrLib CaloIdentifier CxxUtils GeneratorObjects LArG4Code ISF_FastCaloSimEvent Identifier ) + +# Install files from the package: +atlas_install_headers( ISF_FastCaloSimSD ) +atlas_install_python_modules( python/*.py ) + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/python/ISF_FastCaloSimSDConfig.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/python/ISF_FastCaloSimSDConfig.py new file mode 100644 index 0000000000000000000000000000000000000000..a8dc3f98ead00a0b408104c3b30aa00fcf3dc39f --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/python/ISF_FastCaloSimSDConfig.py @@ -0,0 +1,31 @@ +# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + +from AthenaCommon import CfgMgr + +def getFCS_StepInfoSensitiveDetector(name="FCS_StepInfoSensitiveDetector", **kwargs): + ## Main configuration + kwargs.setdefault("StacVolumes",["LArMgr::LAr::EMB::STAC"]) + kwargs.setdefault("PresamplerVolumes",["LArMgr::LAr::Barrel::Presampler::Module"]) + kwargs.setdefault("NegIWVolumes",["LArMgr::LAr::EMEC::Neg::InnerWheel"]) + kwargs.setdefault("NegOWVolumes",["LArMgr::LAr::EMEC::Neg::OuterWheel"]) + kwargs.setdefault("BOBarretteVolumes",["LArMgr::LAr::EMEC::BackOuterBarrette::Module::Phidiv"]) + kwargs.setdefault("MiniVolumes",["LArMgr::MiniFCAL::Wafer"]) + kwargs.setdefault("PosIWVolumes",["LArMgr::LAr::EMEC::Pos::InnerWheel"]) + kwargs.setdefault("PosOWVolumes",["LArMgr::LAr::EMEC::Pos::OuterWheel"]) + kwargs.setdefault("PresVolumes", ["LArMgr::LAr::Endcap::Presampler::LiquidArgon"]) + kwargs.setdefault("SliceVolumes",["LArMgr::LAr::HEC::Module::Depth::Slice"]) + kwargs.setdefault("FCAL1Volumes",["LArMgr::LAr::FCAL::Module1::Gap"]) + kwargs.setdefault("FCAL2Volumes",["LArMgr::LAr::FCAL::Module2::Gap"]) + kwargs.setdefault("FCAL3Volumes",["LArMgr::LAr::FCAL::Module3::Gap"]) + kwargs.setdefault("TileVolumes",["Tile::Scintillator"]) + kwargs.setdefault("OutputCollectionNames", ["EventSteps"]) + # TODO Add extra configuration here!! + return CfgMgr.FCS_Param__FCS_StepInfoSDTool(name, **kwargs) + +##def getFastCaloSimParamActionTool(name='G4UA::FastCaloSimParamActionTool', **kwargs): +## from G4AtlasApps.SimFlags import simFlags +## # example custom configuration +## if name in simFlags.UserActionConfig.get_Value().keys(): +## for prop,value in simFlags.UserActionConfig.get_Value()[name].iteritems(): +## kwargs.setdefault(prop,value) +## return CfgMgr.G4UA__FastCaloSimParamActionTool(name,**kwargs) diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/python/ISF_FastCaloSimSDConfigDb.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/python/ISF_FastCaloSimSDConfigDb.py new file mode 100644 index 0000000000000000000000000000000000000000..4c205253ba04f0718f153a112c49c46d2cce888d --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/python/ISF_FastCaloSimSDConfigDb.py @@ -0,0 +1,5 @@ +# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + +from AthenaCommon.CfgGetter import addTool + +addTool("ISF_FastCaloSimSD.ISF_FastCaloSimSDConfig.getFCS_StepInfoSensitiveDetector", "FCS_StepInfoSensitiveDetector") diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSD.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSD.cxx new file mode 100644 index 0000000000000000000000000000000000000000..9e46bd508c0182f256f3c02e31ce168746bb4068 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSD.cxx @@ -0,0 +1,227 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "FCS_StepInfoSD.h" + +#include "CaloIdentifier/CaloIdManager.h" +#include "CaloIdentifier/LArID_Exception.h" +#include "CaloIdentifier/LArEM_ID.h" +#include "CaloIdentifier/LArFCAL_ID.h" +#include "CaloIdentifier/LArHEC_ID.h" +#include "CaloIdentifier/LArMiniFCAL_ID.h" +#include "CaloDetDescr/CaloDetDescrManager.h" +#include "CaloDetDescr/CaloDetDescrElement.h" + +#include "G4Step.hh" +#include "G4ThreeVector.hh" + +FCS_StepInfoSD::FCS_StepInfoSD(G4String a_name, const FCS_Param::Config& config) + : G4VSensitiveDetector(a_name) + , m_config(config) + , m_larEmID(nullptr) + , m_larFcalID(nullptr) + , m_larHecID(nullptr) + , m_larMiniFcalID(nullptr) + , m_calo_dd_man(nullptr) +{ + m_calo_dd_man = CaloDetDescrManager::instance(); //FIXME Move somewhere else!! +} + +FCS_StepInfoSD::~FCS_StepInfoSD() +{ +} + +G4bool FCS_StepInfoSD::ProcessHits(G4Step*,G4TouchableHistory*) +{ + G4ExceptionDescription description; + description << "ProcessHits: Base class method should not be called!!!"; + G4Exception("FCS_StepInfoSD", "FCSBadCall", FatalException, description); + abort(); + return false; +} + +inline double FCS_StepInfoSD::getMaxTime(const CaloCell_ID::CaloSample& layer) const +{ + /// NB The result of this function should actually be constant for each SD + if (layer >= CaloCell_ID::PreSamplerB && layer <= CaloCell_ID::EME3) { + return m_config.m_maxTimeLAr; + } + if (layer >= CaloCell_ID::HEC0 && layer <= CaloCell_ID::HEC3) { + return m_config.m_maxTimeHEC; + } + if (layer >=CaloCell_ID::FCAL0 && layer <= CaloCell_ID::FCAL2) { + return m_config.m_maxTimeFCAL; + } + return m_config.m_maxTime; +} + +inline double FCS_StepInfoSD::getMaxRadius(const CaloCell_ID::CaloSample& layer) const +{ + /// NB The result of this function should actually be constant for each SD + if (layer >= CaloCell_ID::PreSamplerB && layer <= CaloCell_ID::EME3) { + return m_config.m_maxRadiusLAr; + } + if (layer >= CaloCell_ID::HEC0 && layer <= CaloCell_ID::HEC3) { + return m_config.m_maxRadiusHEC; + } + if (layer >=CaloCell_ID::FCAL0 && layer <= CaloCell_ID::FCAL2) { + return m_config.m_maxRadiusFCAL; + } + return m_config.m_maxRadius; +} + +inline double FCS_StepInfoSD::getMaxDeltaR(const CaloCell_ID::CaloSample& layer) const +{ + /// NB The result of this function should actually be constant for each SD + if(m_config.m_maxRadiusLAr != 25) return 999.; + if (layer >= CaloCell_ID::PreSamplerB && layer <= CaloCell_ID::EME3) { + if (layer==CaloCell_ID::PreSamplerB || layer==CaloCell_ID::PreSamplerE) { + // PS default is 1mm in eta, 5mm in phi, no cut in r + return m_config.m_maxrPS; + } + else if (layer==CaloCell_ID::EMB1 || layer==CaloCell_ID::EME1) { + // EM1 default is 1mm in eta, 5mm in phi, 15mm in r + return m_config.m_maxrEM1; + } + else if (layer==CaloCell_ID::EMB2 || layer==CaloCell_ID::EME2) { + // EM2 default is 1mm in eta, 5mm in phi, 60mm in r + return m_config.m_maxrEM2; + } + else if (layer==CaloCell_ID::EMB3 || layer==CaloCell_ID::EME3) { + // EM3 default is 1mm in eta, 5mm in phi, 8mm in r + return m_config.m_maxrEM3; + } + } + return 999.; +} + +inline double FCS_StepInfoSD::getMaxDeltaEta(const CaloCell_ID::CaloSample& layer) const +{ + /// NB The result of this function should actually be constant for each SD + if(m_config.m_maxRadiusLAr != 25) return 999.; + if (layer >= CaloCell_ID::PreSamplerB && layer <= CaloCell_ID::EME3) { + if (layer==CaloCell_ID::PreSamplerB || layer==CaloCell_ID::PreSamplerE) { + // PS default is 1mm in eta, 5mm in phi, no cut in r + return m_config.m_maxEtaPS; + } + else if (layer==CaloCell_ID::EMB1 || layer==CaloCell_ID::EME1) { + // EM1 default is 1mm in eta, 5mm in phi, 15mm in r + return m_config.m_maxEtaEM1; + } + else if (layer==CaloCell_ID::EMB2 || layer==CaloCell_ID::EME2) { + // EM2 default is 1mm in eta, 5mm in phi, 60mm in r + return m_config.m_maxEtaEM2; + } + else if (layer==CaloCell_ID::EMB3 || layer==CaloCell_ID::EME3) { + // EM3 default is 1mm in eta, 5mm in phi, 8mm in r + return m_config.m_maxEtaEM3; + } + } + return 999.; +} + +inline double FCS_StepInfoSD::getMaxDeltaPhi(const CaloCell_ID::CaloSample& layer) const +{ + /// NB The result of this function should actually be constant for each SD + if(m_config.m_maxRadiusLAr != 25) return 999.; + if (layer >= CaloCell_ID::PreSamplerB && layer <= CaloCell_ID::EME3) { + if (layer==CaloCell_ID::PreSamplerB || layer==CaloCell_ID::PreSamplerE) { + // PS default is 1mm in eta, 5mm in phi, no cut in r + return m_config.m_maxPhiPS; + } + else if (layer==CaloCell_ID::EMB1 || layer==CaloCell_ID::EME1) { + // EM1 default is 1mm in eta, 5mm in phi, 15mm in r + return m_config.m_maxPhiEM1; + } + else if (layer==CaloCell_ID::EMB2 || layer==CaloCell_ID::EME2) { + // EM2 default is 1mm in eta, 5mm in phi, 60mm in r + return m_config.m_maxPhiEM2; + } + else if (layer==CaloCell_ID::EMB3 || layer==CaloCell_ID::EME3) { + // EM3 default is 1mm in eta, 5mm in phi, 8mm in r + return m_config.m_maxPhiEM3; + } + } + return 999.; +} + +void FCS_StepInfoSD::update_map(const CLHEP::Hep3Vector & l_vec, const Identifier & l_cell, double l_energy, double l_time, bool l_valid, int l_detector) +{ + // Drop any hits that don't have a good identifier attached + if (!m_calo_dd_man->get_element(l_cell)) { return; } + + auto map_item = m_hit_map.find( l_cell ); + if (map_item==m_hit_map.end()) { + m_hit_map[l_cell] = new std::vector< ISF_FCS_Parametrization::FCS_StepInfo* >; + m_hit_map[l_cell]->reserve(200); + m_hit_map[l_cell]->push_back( new ISF_FCS_Parametrization::FCS_StepInfo( l_vec , l_cell , l_energy , l_time , l_valid , l_detector ) ); + } + else { + + // Get the appropriate merging limits + const CaloCell_ID::CaloSample& layer = m_calo_dd_man->get_element(l_cell)->getSampling(); + const double tsame(this->getMaxTime(layer)); + const double maxRadius(this->getMaxRadius(layer)); + const double maxDeltaR(this->getMaxDeltaR(layer)); + const double maxDeltaEta(this->getMaxDeltaEta(layer)); + const double maxDeltaPhi(this->getMaxDeltaPhi(layer)); + bool match = false; + for (auto map_it : * map_item->second) { + // Time check is straightforward + const double delta_t = std::fabs(map_it->time()-l_time); + if ( delta_t >= tsame ) { continue; } + + // Distance check is less straightforward + const CLHEP::Hep3Vector & currentPosition = map_it->position(); + const double hit_diff2 = currentPosition.diff2( l_vec ); + + // Overall merging scheme + if (layer >= CaloCell_ID::PreSamplerB && layer <= CaloCell_ID::EME3) { + // Customized merging scheme for LAr barrel and endcap, use only if we're not changing maxRadiusLAr value + if(m_config.m_maxRadiusLAr == 25) { + const CLHEP::Hep3Vector a_diff = l_vec - currentPosition; + const double a_diffmag(a_diff.mag()); + const double a_inv_length = 1./a_diffmag; + const double delta_r = std::fabs(a_diff.dot( l_vec ) * a_inv_length); + if (delta_r >= maxDeltaR) { continue; } + const double delta_eta = std::fabs(std::sin(l_vec.theta()-currentPosition.theta())*a_diffmag); + if (delta_eta >= maxDeltaEta) { continue; } + const double delta_phi = std::fabs(std::sin(l_vec.phi()-currentPosition.phi())*a_diffmag); + if (delta_phi >= maxDeltaPhi) { continue; } + } + else { + if ( hit_diff2 >= maxRadius ) { continue; } + } + } + else if ( hit_diff2 >= maxRadius ) { continue; } + // Found a match. Make a temporary that will be deleted! + const ISF_FCS_Parametrization::FCS_StepInfo my_info( l_vec , l_cell , l_energy , l_time , l_valid , l_detector ); + *map_it += my_info; + match = true; + break; + } // End of search for match in time and space + if (!match) { + map_item->second->push_back( new ISF_FCS_Parametrization::FCS_StepInfo( l_vec , l_cell , l_energy , l_time , l_valid , l_detector ) ); + } // Didn't match + } // ID already in the map + return; +} // That's it for updating the map! + +void FCS_StepInfoSD::EndOfAthenaEvent( ISF_FCS_Parametrization::FCS_StepInfoCollection* hitContainer ) +{ + // Unpack map into vector + for (auto it : m_hit_map) { + for (auto a_s : * it.second) { + // Giving away ownership of the objects! + hitContainer->push_back( a_s ); + } + it.second->clear(); + delete it.second; + } // Vector of IDs in the map + m_hit_map.clear(); + if (m_config.verboseLevel > 4) { + G4cout <<this->GetName()<< " DEBUG EndOfAthenaEvent: After initial cleanup, N=" << hitContainer->size() << G4endl; + } + return; +} diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSD.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSD.h new file mode 100644 index 0000000000000000000000000000000000000000..1a4c3410976fc92b0ea63abb5afeec716db9e01d --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSD.h @@ -0,0 +1,143 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef ISF_FASTCALOSIM_FCS_STEPINFOSD_H +#define ISF_FASTCALOSIM_FCS_STEPINFOSD_H + +#include "G4VSensitiveDetector.hh" + +#include "CaloIdentifier/CaloCell_ID.h" // For CaloCell_ID::CaloSample +#include "LArG4Code/LArG4Identifier.h" + +#include "LArSimEvent/LArHit.h" +#include "CLHEP/Units/SystemOfUnits.h" + +#include "ISF_FastCaloSimEvent/FCS_StepInfoCollection.h" + +#include <map> +#include <vector> + +// Forward declarations +class G4Step; +class G4TouchableHistory; + +class LArEM_ID; +class LArFCAL_ID; +class LArHEC_ID; +class LArMiniFCAL_ID; +class CaloDetDescrManager; + +class ILArCalculatorSvc; +class ITileCalculator; +class LArHitContainer; + +class StoreGateSvc; + +namespace FCS_Param { + + struct Config + { + /** Helper to keep the same verbosity everywhere */ + int verboseLevel=0; + bool shift_lar_subhit=true; + bool shorten_lar_step=false; + + // Merging properties + double m_maxRadius=25.; //!< property, see @link LArG4GenShowerLib::LArG4GenShowerLib @endlink + double m_maxRadiusLAr=25.; //!< property, see @link LArG4GenShowerLib::LArG4GenShowerLib @endlink + double m_maxRadiusHEC=100.; //!< property, see @link LArG4GenShowerLib::LArG4GenShowerLib @endlink + double m_maxRadiusFCAL=100.; //!< property, see @link LArG4GenShowerLib::LArG4GenShowerLib @endlink + double m_maxRadiusTile=100.; //!< property, see @link LArG4GenShowerLib::LArG4GenShowerLib @endlink + + double m_maxTime=25.; + double m_maxTimeLAr=25.; + double m_maxTimeHEC=25.; + double m_maxTimeFCAL=25.; + double m_maxTimeTile=25.; + + // Optimised merging scheme + double m_maxEtaPS=1.; + double m_maxPhiPS=5.; + double m_maxrPS=0.; + + double m_maxEtaEM1=1.; + double m_maxPhiEM1=5.; + double m_maxrEM1=15.; + + double m_maxEtaEM2=1.; + double m_maxPhiEM2=5.; + double m_maxrEM2=60.; + + double m_maxEtaEM3=1.; + double m_maxPhiEM3=5.; + double m_maxrEM3=8.; + + ILArCalculatorSvc *m_LArCalculator=nullptr; + ITileCalculator *m_TileCalculator=nullptr; + }; + +} + +/// @class FCS_StepInfoSD +/// @brief Common sensitive detector class for LAr systems. +/// +/// This SD implementation saves the standard LArHits. +/// See LArG4CalibSD for an SD that handles calibration hits. +/// +class FCS_StepInfoSD : public G4VSensitiveDetector +{ +public: + + /// Constructor + FCS_StepInfoSD(G4String a_name, const FCS_Param::Config& config); + + /// Destructor + virtual ~FCS_StepInfoSD(); + + /// Main processing method + G4bool ProcessHits(G4Step* a_step, G4TouchableHistory*) override; + + /// End of athena event processing + void EndOfAthenaEvent( ISF_FCS_Parametrization::FCS_StepInfoCollection* hitContnainer ); + + /// Sets the ID helper pointers + void setupHelpers( const LArEM_ID* EM , + const LArFCAL_ID* FCAL , + const LArHEC_ID* HEC , + const LArMiniFCAL_ID* mini ) { + m_larEmID = EM; + m_larFcalID = FCAL; + m_larHecID = HEC; + m_larMiniFcalID = mini; + } + +protected: + /// Keep a map instead of trying to keep the full vector. + /// At the end of the event we'll push the map back into the + /// FCS_StepInfoCollection in StoreGate. + virtual void update_map(const CLHEP::Hep3Vector & l_vec, const Identifier & l_cell, double l_energy, double l_time, bool l_valid, int l_detector); + FCS_Param::Config m_config; + /// Pointers to the identifier helpers + const LArEM_ID* m_larEmID; + const LArFCAL_ID* m_larFcalID; + const LArHEC_ID* m_larHecID; + const LArMiniFCAL_ID* m_larMiniFcalID; + const CaloDetDescrManager *m_calo_dd_man; + std::map< Identifier , std::vector< ISF_FCS_Parametrization::FCS_StepInfo* >* > m_hit_map; + +private: + /// + double getMaxTime(const CaloCell_ID::CaloSample& layer) const; + /// + double getMaxRadius(const CaloCell_ID::CaloSample& layer) const; + /// + double getMaxDeltaR(const CaloCell_ID::CaloSample& layer) const; + /// + double getMaxDeltaEta(const CaloCell_ID::CaloSample& layer) const; + /// + double getMaxDeltaPhi(const CaloCell_ID::CaloSample& layer) const; + +}; + +#endif // ISF_FASTCALOSIM_FCS_STEPINFOSD_H diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSDTool.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSDTool.cxx new file mode 100644 index 0000000000000000000000000000000000000000..a4e4660568e51214be86dc03cd1bf0185446c2e9 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSDTool.cxx @@ -0,0 +1,290 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "FCS_StepInfoSDTool.h" + +#include <memory> + +// External includes +#include "CLHEP/Units/SystemOfUnits.h" + +// Athena includes +#include "CaloIdentifier/CaloIdManager.h" +#include "CaloIdentifier/LArEM_ID.h" +#include "CaloIdentifier/LArFCAL_ID.h" +#include "CaloIdentifier/LArHEC_ID.h" +#include "CaloIdentifier/LArMiniFCAL_ID.h" +#include "LArG4Code/ILArCalculatorSvc.h" +#include "LArG4Code/VolumeUtils.h" +#include "TileG4Interfaces/ITileCalculator.h" + +// Local includes +#include "LArFCS_StepInfoSD.h" +#include "TileFCS_StepInfoSD.h" +#include "SDWrapper.h" + +/// Template instantiation for FCS_StepInfoSD +using FCS_StepInfoSDWrapper = FCS_Param::detail::SDWrapper<FCS_StepInfoSD, ISF_FCS_Parametrization::FCS_StepInfoCollection>; + +namespace FCS_Param +{ + + //--------------------------------------------------------------------------- + // Tool constructor + //--------------------------------------------------------------------------- + FCS_StepInfoSDTool::FCS_StepInfoSDTool(const std::string& type, const std::string& name, + const IInterface* parent) + : SensitiveDetectorBase(type, name, parent) + , m_hitCollName("EventSteps") + , m_bpsmodcalc("EMBPresamplerCalculator", name) + , m_embcalc("EMBCalculator", name) + , m_emepiwcalc("EMECPosInnerWheelCalculator", name) + , m_emeniwcalc("EMECNegInnerWheelCalculator", name) + , m_emepowcalc("EMECPosOuterWheelCalculator", name) + , m_emenowcalc("EMECNegOuterWheelCalculator", name) + , m_emepscalc("EMECPresamplerCalculator", name) + , m_emeobarcalc("EMECBackOuterBarretteCalculator", name) + , m_heccalc("HECWheelCalculator", name) + , m_fcal1calc("FCAL1Calculator", name) + , m_fcal2calc("FCAL2Calculator", name) + , m_fcal3calc("FCAL3Calculator", name) + , m_minfcalcalc("MiniFCALCalculator", name) + , m_tileCalculator("TileGeoG4SDCalc", name) + , m_larEmID(nullptr) + , m_larFcalID(nullptr) + , m_larHecID(nullptr) + , m_larMiniFcalID(nullptr) + , m_config() + { + declareProperty("HitCollectionName", m_hitCollName); + declareProperty("StacVolumes", m_stacVolumes); + declareProperty("PresamplerVolumes", m_presBarVolumes); + declareProperty("PosIWVolumes", m_posIWVolumes); + declareProperty("NegIWVolumes", m_negIWVolumes); + declareProperty("PosOWVolumes", m_posOWVolumes); + declareProperty("NegOWVolumes", m_negOWVolumes); + declareProperty("PresVolumes", m_presECVolumes); + declareProperty("BOBarretteVolumes", m_bobVolumes); + declareProperty("FCAL1Volumes", m_fcal1Volumes); + declareProperty("FCAL2Volumes", m_fcal2Volumes); + declareProperty("FCAL3Volumes", m_fcal3Volumes); + declareProperty("SliceVolumes", m_sliceVolumes); + declareProperty("MiniVolumes", m_miniVolumes); + declareProperty("TileVolumes", m_tileVolumes); + + declareProperty("EMBPSCalibrationCalculator",m_bpsmodcalc); + declareProperty("EMBCalibrationCalculator",m_embcalc); + declareProperty("EMECPosIWCalibrationCalculator",m_emepiwcalc); + declareProperty("EMECNegIWCalibrationCalculator",m_emeniwcalc); + declareProperty("EMECPosOWCalibrationCalculator",m_emepowcalc); + declareProperty("EMECNegOWCalibrationCalculator",m_emenowcalc); + declareProperty("EMECPSCalibrationCalculator",m_emepscalc); + declareProperty("EMECBOBCalibrationCalculator",m_emeobarcalc); + declareProperty("HECWActiveCalculator",m_heccalc); + declareProperty("FCAL1CalibCalculator",m_fcal1calc); + declareProperty("FCAL2CalibCalculator",m_fcal2calc); + declareProperty("FCAL3CalibCalculator",m_fcal3calc); + declareProperty("MiniFCALActiveCalibrationCalculator",m_minfcalcalc); + declareProperty("TileCalculator", m_tileCalculator); + + declareProperty("shift_lar_subhit",m_config.shift_lar_subhit, ""); + declareProperty("shorten_lar_step",m_config.shorten_lar_step, ""); + + declareProperty("maxRadius",m_config.m_maxRadius, ""); + declareProperty("maxRadiusLAr",m_config.m_maxRadiusLAr, ""); + declareProperty("maxRadiusHEC",m_config.m_maxRadiusHEC, ""); + declareProperty("maxRadiusFCAL",m_config.m_maxRadiusFCAL, ""); + declareProperty("maxRadiusTile",m_config.m_maxRadiusTile, ""); + + declareProperty("maxTime",m_config.m_maxTime, ""); + declareProperty("maxTimeLAr",m_config.m_maxTimeLAr, ""); + declareProperty("maxTimeHEC",m_config.m_maxTimeHEC, ""); + declareProperty("maxTimeFCAL",m_config.m_maxTimeFCAL, ""); + declareProperty("maxTimeTile",m_config.m_maxTimeTile, ""); + + declareProperty("maxEtaPS", m_config.m_maxEtaPS, ""); + declareProperty("maxPhiPS", m_config.m_maxPhiPS, ""); + declareProperty("maxrPS", m_config.m_maxrPS, ""); + + declareProperty("maxEtaEM1", m_config.m_maxEtaEM1, ""); + declareProperty("maxPhiEM1", m_config.m_maxPhiEM1, ""); + declareProperty("maxrEM1", m_config.m_maxrEM1, ""); + + declareProperty("maxEtaEM2", m_config.m_maxEtaEM2, ""); + declareProperty("maxPhiEM2", m_config.m_maxPhiEM2, ""); + declareProperty("maxrEM2", m_config.m_maxrEM2, ""); + + declareProperty("maxEtaEM3", m_config.m_maxEtaEM3, ""); + declareProperty("maxPhiEM3", m_config.m_maxPhiEM3, ""); + declareProperty("maxrEM3", m_config.m_maxrEM3, ""); + } + + //--------------------------------------------------------------------------- + // Initialize the tool + //--------------------------------------------------------------------------- + StatusCode FCS_StepInfoSDTool::initialize() + { + ATH_MSG_DEBUG( "Initializing " << name() ); + + const CaloIdManager* idMgr = nullptr; + CHECK( detStore()->retrieve(idMgr) ); + if( (m_larEmID = idMgr->getEM_ID()) == nullptr) { + ATH_MSG_ERROR("Invalid LAr EM ID helper"); + return StatusCode::FAILURE; + } + if( (m_larFcalID = idMgr->getFCAL_ID()) == nullptr) { + ATH_MSG_ERROR("Invalid LAr FCAL ID helper"); + return StatusCode::FAILURE; + } + if( (m_larHecID = idMgr->getHEC_ID()) == nullptr) { + ATH_MSG_ERROR("Invalid LAr HEC ID helper"); + return StatusCode::FAILURE; + } + if( (m_larMiniFcalID = idMgr->getMiniFCAL_ID()) == nullptr) { + ATH_MSG_ERROR("Invalid LAr Mini FCAL ID helper"); + return StatusCode::FAILURE; + } + + // No general volume list for SensitiveDetectorBase + m_noVolumes = true; + + ATH_CHECK(this->initializeCalculators()); + + return StatusCode::SUCCESS; + } + + //--------------------------------------------------------------------------- + // Initialization of Athena-components + //--------------------------------------------------------------------------- + StatusCode FCS_StepInfoSDTool::initializeCalculators() + { + // Lots of calculators !!! + ATH_CHECK(m_bpsmodcalc.retrieve()); + ATH_CHECK(m_embcalc.retrieve()); + ATH_CHECK(m_emepiwcalc.retrieve()); + ATH_CHECK(m_emeniwcalc.retrieve()); + ATH_CHECK(m_emepowcalc.retrieve()); + ATH_CHECK(m_emenowcalc.retrieve()); + ATH_CHECK(m_emepscalc.retrieve()); + ATH_CHECK(m_emeobarcalc.retrieve()); + ATH_CHECK(m_heccalc.retrieve()); + ATH_CHECK(m_fcal1calc.retrieve()); + ATH_CHECK(m_fcal2calc.retrieve()); + ATH_CHECK(m_fcal3calc.retrieve()); + ATH_CHECK(m_minfcalcalc.retrieve()); + ATH_CHECK(m_tileCalculator.retrieve()); + + return StatusCode::SUCCESS; + } + + //--------------------------------------------------------------------------- + // Collect hits for this event + //--------------------------------------------------------------------------- + StatusCode FCS_StepInfoSDTool::Gather() + { + ATH_MSG_DEBUG("Gathering hits to write out in " << name()); + auto sdWrapper = dynamic_cast<FCS_StepInfoSDWrapper*>( getSD() ); + if(!sdWrapper) { + ATH_MSG_ERROR("Failed to cast SD to FCS_StepInfoSDWrapper"); + return StatusCode::FAILURE; + } + sdWrapper->EndOfAthenaEvent(); + return StatusCode::SUCCESS; + } + + //--------------------------------------------------------------------------- + // Create SD wrapper for current thread + //--------------------------------------------------------------------------- + G4VSensitiveDetector* FCS_StepInfoSDTool::makeSD() + { + // Create the wrapper + auto sdWrapper = new FCS_StepInfoSDWrapper("FCS_StepInfoSDWrapper", m_hitCollName); + + // Create the SDs. + sdWrapper->addSD( makeOneLArSD( "Barrel::Presampler::Module::StepInfo", &*m_bpsmodcalc, m_presBarVolumes ) ); + sdWrapper->addSD( makeOneLArSD( "EMB::STAC::StepInfo", &*m_embcalc, m_stacVolumes ) ); + sdWrapper->addSD( makeOneLArSD( "EMEC::Pos::InnerWheel::StepInfo", &*m_emepiwcalc, m_posIWVolumes ) ); + sdWrapper->addSD( makeOneLArSD( "EMEC::Neg::InnerWheel::StepInfo", &*m_emeniwcalc, m_negIWVolumes ) ); + sdWrapper->addSD( makeOneLArSD( "EMEC::Pos::OuterWheel::StepInfo", &*m_emepowcalc, m_posOWVolumes ) ); + sdWrapper->addSD( makeOneLArSD( "EMEC::Neg::OuterWheel::StepInfo", &*m_emenowcalc, m_negOWVolumes ) ); + sdWrapper->addSD( makeOneLArSD( "Endcap::Presampler::LiquidArgon::StepInfo", &*m_emepscalc, m_presECVolumes ) ); + sdWrapper->addSD( makeOneLArSD( "EMEC::BackOuterBarrette::StepInfo", &*m_emeobarcalc, m_bobVolumes ) ); + sdWrapper->addSD( makeOneLArSD( "FCAL::Module1::Gap::StepInfo", &*m_fcal1calc, m_fcal1Volumes ) ); + sdWrapper->addSD( makeOneLArSD( "FCAL::Module2::Gap::StepInfo", &*m_fcal2calc, m_fcal2Volumes ) ); + sdWrapper->addSD( makeOneLArSD( "FCAL::Module3::Gap::StepInfo", &*m_fcal3calc, m_fcal3Volumes ) ); + sdWrapper->addSD( makeOneLArSD( "HEC::Module::Depth::Slice::Wheel::StepInfo", &*m_heccalc, m_sliceVolumes ) ); + sdWrapper->addSD( makeOneLArSD( "MiniFCAL::Wafer::StepInfo", &*m_minfcalcalc, m_miniVolumes ) ); + sdWrapper->addSD( makeOneTileSD( "Tile::Scintillator::StepInfo", &*m_tileCalculator, m_tileVolumes ) ); + + return sdWrapper; + } + + //--------------------------------------------------------------------------- + // Create one LAr SD + //--------------------------------------------------------------------------- + std::unique_ptr<FCS_StepInfoSD> + FCS_StepInfoSDTool::makeOneLArSD(const std::string& sdName, ILArCalculatorSvc* calc, + const std::vector<std::string>& volumes) const + { + ATH_MSG_VERBOSE( name() << " makeOneSD" ); + + // Parse the wildcard patterns for existing volume names + auto parsedVolumes = LArG4::findLogicalVolumes(volumes, msg()); // from LArG4Code/VolumeUtils.h + + // Inject the Calculator into m_config + FCS_Param::Config config(m_config); + config.m_LArCalculator = calc; + if(!calc) { + throw GaudiException("nullptr for ILArCalculatorSvc provided to constructor for: " + sdName, + name(), StatusCode::FAILURE); + } + if(msgLvl(MSG::VERBOSE)) { config.verboseLevel = 10; } + else if(msgLvl(MSG::DEBUG)) { config.verboseLevel = 5; } + // Create the simple SD + std::unique_ptr<FCS_StepInfoSD> sd = + std::make_unique<LArFCS_StepInfoSD>(sdName, config); + sd->setupHelpers(m_larEmID, m_larFcalID, m_larHecID, m_larMiniFcalID); + + // Assign the volumes to the SD + if( this->assignSD( sd.get(), parsedVolumes ).isFailure() ) { + // TODO: can I just return NULL here? + throw GaudiException("Failed to assign sd: " + sdName, + name(), StatusCode::FAILURE); + } + return std::move(sd); + } + + //--------------------------------------------------------------------------- + // Create one Tile SD + //--------------------------------------------------------------------------- + std::unique_ptr<FCS_StepInfoSD> + FCS_StepInfoSDTool::makeOneTileSD(const std::string& sdName, ITileCalculator* calc, + const std::vector<std::string>& volumes) const + { + ATH_MSG_VERBOSE( name() << " makeOneSD" ); + + // Inject the Calculator into m_config + FCS_Param::Config config(m_config); + config.m_TileCalculator = calc; + if(!calc) { + throw GaudiException("nullptr for ITileCalculator provided to constructor for: " + sdName, + name(), StatusCode::FAILURE); + } + if(msgLvl(MSG::VERBOSE)) { config.verboseLevel = 10; } + else if(msgLvl(MSG::DEBUG)) { config.verboseLevel = 5; } + // Create the simple SD + std::unique_ptr<FCS_StepInfoSD> sd = + std::make_unique<TileFCS_StepInfoSD>(sdName, config); + sd->setupHelpers(m_larEmID, m_larFcalID, m_larHecID, m_larMiniFcalID); + + // Assign the volumes to the SD + if( this->assignSD( sd.get(), volumes ).isFailure() ) { + // TODO: can I just return NULL here? + throw GaudiException("Failed to assign sd: " + sdName, + name(), StatusCode::FAILURE); + } + return std::move(sd); + } + +} // namespace FCS_Param diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSDTool.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSDTool.h new file mode 100644 index 0000000000000000000000000000000000000000..1bd88ab1e4b149129fa1b54b503e4cb2f74c467f --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/FCS_StepInfoSDTool.h @@ -0,0 +1,123 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef ISF_FASTCALOSIMSD_FCS_STEPINFOSDTOOL_H +#define ISF_FASTCALOSIMSD_FCS_STEPINFOSDTOOL_H + +/// @file FCS_StepInfoSDTool.h +/// @brief Defines the FCS_StepInfoSDTool class +/// @author Steve Farrell <Steven.Farrell@cern.ch> +/// @date 2016-03-26 + +// System includes +#include <string> +#include <vector> + +// G4Atlas includes +#include "G4AtlasTools/SensitiveDetectorBase.h" + +// Local includes +#include "FCS_StepInfoSD.h" + +// Forward declarations +class ILArCalculatorSvc; +class LArEM_ID; +class LArFCAL_ID; +class LArHEC_ID; +class LArMiniFCAL_ID; +class ITileCalculator; + +namespace FCS_Param +{ + + /// @class FCS_StepInfoSDTool + /// @brief A base class for tools that manage FCS_StepInfoSDs. + /// + /// @todo Add more details. + /// + /// @author Steve Farrell <Steven.Farrell@cern.ch> + /// + class FCS_StepInfoSDTool : public SensitiveDetectorBase + { + + public: + + /// Constructor + FCS_StepInfoSDTool(const std::string& type, const std::string& name, + const IInterface* parent); + + /// Initialize the tool + StatusCode initialize() override final; + + /// Calls down to all the SDs to pack their hits into one collection + StatusCode Gather() override final; + + private: + + /// Create the SD wrapper for current worker thread + G4VSensitiveDetector* makeSD() override final; + + /// Initialize Calculator Services + virtual StatusCode initializeCalculators(); + + /// Helper method to create one SD + std::unique_ptr<FCS_StepInfoSD> + makeOneLArSD(const std::string& name, ILArCalculatorSvc* calc, + const std::vector<std::string>& volumes) const; + + /// Helper method to create one SD + std::unique_ptr<FCS_StepInfoSD> + makeOneTileSD(const std::string& name, ITileCalculator* calc, + const std::vector<std::string>& volumes) const; + + /// Hit collection name + std::string m_hitCollName; + + /// @name SD volumes + /// @{ + std::vector<std::string> m_stacVolumes; + std::vector<std::string> m_presBarVolumes; + std::vector<std::string> m_posIWVolumes; + std::vector<std::string> m_negIWVolumes; + std::vector<std::string> m_posOWVolumes; + std::vector<std::string> m_negOWVolumes; + std::vector<std::string> m_presECVolumes; + std::vector<std::string> m_bobVolumes; + std::vector<std::string> m_fcal1Volumes; + std::vector<std::string> m_fcal2Volumes; + std::vector<std::string> m_fcal3Volumes; + std::vector<std::string> m_sliceVolumes; + std::vector<std::string> m_miniVolumes; + std::vector<std::string> m_tileVolumes; + /// @} + + ServiceHandle<ILArCalculatorSvc> m_bpsmodcalc; //LArG4::BarrelPresampler::CalibrationCalculator + ServiceHandle<ILArCalculatorSvc> m_embcalc; //LArG4::Barrel::CalibrationCalculator + ServiceHandle<ILArCalculatorSvc> m_emepiwcalc; //LArG4::EC::CalibrationCalculator(LArWheelCalculator::InnerAbsorberWheel, 1) + ServiceHandle<ILArCalculatorSvc> m_emeniwcalc; //LArG4::EC::CalibrationCalculator(LArWheelCalculator::InnerAbsorberWheel, -1) + ServiceHandle<ILArCalculatorSvc> m_emepowcalc; //LArG4::EC::CalibrationCalculator(LArWheelCalculator::OuterAbsorberWheel, 1) + ServiceHandle<ILArCalculatorSvc> m_emenowcalc; //LArG4::EC::CalibrationCalculator(LArWheelCalculator::OuterAbsorberWheel, -1) + ServiceHandle<ILArCalculatorSvc> m_emepscalc; //LArG4::EC::PresamplerCalibrationCalculator + ServiceHandle<ILArCalculatorSvc> m_emeobarcalc; //LArG4::EC::CalibrationCalculator(LArWheelCalculator::BackOuterBarretteWheelCalib, 1) + ServiceHandle<ILArCalculatorSvc> m_heccalc; //LArG4::HEC::LArHECCalibrationWheelCalculator(LArG4::HEC::kWheelActive) + ServiceHandle<ILArCalculatorSvc> m_fcal1calc; + ServiceHandle<ILArCalculatorSvc> m_fcal2calc; + ServiceHandle<ILArCalculatorSvc> m_fcal3calc; + ServiceHandle<ILArCalculatorSvc> m_minfcalcalc; //LArG4::MiniFCAL::MiniFCALCalibrationCalculator(LArG4::MiniFCAL::kActive) + ServiceHandle<ITileCalculator> m_tileCalculator; + + /// @name Calo identifier helpers + /// @{ + const LArEM_ID* m_larEmID; + const LArFCAL_ID* m_larFcalID; + const LArHEC_ID* m_larHecID; + const LArMiniFCAL_ID* m_larMiniFcalID; + /// @} + FCS_Param::Config m_config; + + }; // class FCS_StepInfoSDTool + +} // namespace FCS_Param + +#endif diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/LArFCS_StepInfoSD.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/LArFCS_StepInfoSD.cxx new file mode 100644 index 0000000000000000000000000000000000000000..db10072268276ec8cfabae3cab0a751ca6abfa93 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/LArFCS_StepInfoSD.cxx @@ -0,0 +1,274 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// class header +#include "LArFCS_StepInfoSD.h" +/// Geant4 headers +#include "G4Step.hh" +#include "G4ThreeVector.hh" +#include "G4TouchableHistory.hh" + +/// Athena headers +#include "CaloDetDescr/CaloDetDescrManager.h" +#include "CaloDetDescr/CaloDetDescrElement.h" +#include "LArG4Code/ILArCalculatorSvc.h" + +LArFCS_StepInfoSD::LArFCS_StepInfoSD(G4String a_name, const FCS_Param::Config& config) + : FCS_StepInfoSD(a_name, config) + , m_calculator(m_config.m_LArCalculator) +{ +} + +LArFCS_StepInfoSD::~LArFCS_StepInfoSD() +{ +} + +G4bool LArFCS_StepInfoSD::ProcessHits(G4Step* a_step,G4TouchableHistory*) +{ + G4bool result(false); + // If there's no energy, there's no hit. (Aside: Isn't this energy + // the same as the energy from the calculator? Not necessarily. + // The calculator may include detector effects such as + // charge-collection which are not modeled by Geant4.) + if(a_step->GetTotalEnergyDeposit() <= 0.) { return result; } + + if (m_calculator) { + const double StepLength = a_step->GetStepLength()/ CLHEP::mm; + const G4ThreeVector preStepPosition = a_step->GetPreStepPoint()->GetPosition(); //pre step is the position we're interested in + const G4ThreeVector postStepPosition = a_step->GetPostStepPoint()->GetPosition(); + double et(0.); // Total collected charge + std::vector<const G4Step*> steps; + bool madeSubSteps(false); + if (m_config.shorten_lar_step && StepLength>0.2) { + //create smaller substeps instead + G4int nsub_step=(int) (StepLength/0.2) + 1; + G4double delta=1./((double) nsub_step); + //G4cout <<"Orig prestep: "<<preStepPosition<<std::endl; + for (G4int i=0;i<nsub_step;i++) { + G4double fraction1 = ((G4double) i)*delta; + G4double fraction2 = (((G4double) i) + 1.)*delta; + G4ThreeVector subpoint1=preStepPosition*(1-fraction1) + postStepPosition*fraction1; + G4ThreeVector subpoint2=preStepPosition*(1-fraction2) + postStepPosition*fraction2; + G4StepPoint *startpoint = new G4StepPoint(*(a_step->GetPreStepPoint())); + G4StepPoint *endpoint = new G4StepPoint(*(a_step->GetPostStepPoint())); + startpoint->SetPosition(subpoint1); + endpoint->SetPosition(subpoint2); + + G4Step* newstep = new G4Step(*a_step); + newstep->SetPreStepPoint(startpoint); + newstep->SetPostStepPoint(endpoint); + newstep->SetStepLength( (subpoint1-subpoint2).mag()); + newstep->SetTotalEnergyDeposit(a_step->GetTotalEnergyDeposit()/nsub_step); + steps.push_back(newstep); + madeSubSteps=true; + } + } + else { + steps.push_back(a_step); + } + + for (const G4Step* substep : steps) { + G4ThreeVector stepPosition = 0.5*(substep->GetPreStepPoint()->GetPosition()+substep->GetPostStepPoint()->GetPosition()); + std::vector<LArHitData> processedHits; + if (m_calculator->Process(substep, processedHits)) { + for (const auto& larhit: processedHits) { + et += larhit.energy; + } + } + else { + if (m_config.verboseLevel > 4) { + //Maybe 0 hits or something like that... + G4cout << this->GetName()<<" WARNING ProcessHits: Call to ILArCalculatorSvc::Process failed! Details:" << G4endl + << " " << "Volume: "<< a_step->GetPreStepPoint()->GetPhysicalVolume()->GetName()<<" "<<m_calculator<< " position: "<<stepPosition<<" SL: "<<StepLength<<G4endl + << " " << "Orig position: "<<substep->GetPreStepPoint()->GetPosition()<<" / "<<substep->GetPostStepPoint()->GetPosition()<<G4endl + << " " << "StepLength: "<<StepLength<<" step: "<<preStepPosition<<" / "<<postStepPosition<<G4endl; + } + return result; + } + + // drop hits with zero deposited energy (could still happen with negative corrections from calculator) + //Or if total energy is <0 + if (et <= 0.) { + if (m_config.verboseLevel > 4) { + G4cout << this->GetName()<<" WARNING ProcessHits: Total negative energy: " << et << " not processing..." << G4endl; + } + return result; + } + + const size_t numberOfProcessedHits = processedHits.size(); + const G4ThreeVector originalStepPosition = stepPosition; + double maxSubHitEnergy = -999.; + int maxSubHitEnergyindex =-1; + if (numberOfProcessedHits>0) { + maxSubHitEnergy = processedHits[0].energy; + maxSubHitEnergyindex = 0; + } + //Figure out the subhit with most energy + for (size_t i=1; i<numberOfProcessedHits; ++i) { + if (maxSubHitEnergy<processedHits[i].energy) { + maxSubHitEnergy = processedHits[i].energy; + maxSubHitEnergyindex = i; + } + } + //Identifier for the subhit with max energy + Identifier maxEnergyIdentifier = this->ConvertID(processedHits[maxSubHitEnergyindex].id); + CaloDetDescrElement *maxEnergyCell = m_calo_dd_man->get_element(maxEnergyIdentifier); + + Identifier invalidIdentifier; + //for (size_t i=0; i<numberOfProcessedHits; ++i) { + for (const auto& larhit: processedHits) { + Identifier id = this->ConvertID(larhit.id); + if (id == invalidIdentifier) { + G4cout << this->GetName()<<" WARNING ProcessHits: Something wrong with LArG4Identifier: "<<(std::string) larhit.id<<G4endl + <<" "<<id<<G4endl + <<" "<<id.getString()<<G4endl + <<" "<< a_step->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName()<<G4endl + <<" numberOfProcessedHits: "<<numberOfProcessedHits<<G4endl; + G4cout <<" "<<invalidIdentifier<<G4endl; + } + //need to get the cell information + if (numberOfProcessedHits>1) { + if (!m_larEmID->is_em_barrel(id)) { + //It didn't seem to happen outside em_barrel, so flag up if it does: + G4cout << this->GetName()<<" WARNING ProcessHits: Outside LAr barrel, but numberOfProcessedHits="<<numberOfProcessedHits + <<", LArG4Identifier: "<<(std::string) larhit.id<<G4endl; + } + else { + if (m_config.shift_lar_subhit) { + //find subhit with largest energy + if (maxSubHitEnergyindex == -1) { + G4cout << this->GetName()<<" WARNING ProcessHits: no subhit index with e>-999??? "<<G4endl; + return result; + } + if(m_config.verboseLevel > 9) { + G4cout << this->GetName()<<" VERBOSE ProcessHits: shifting subhits: largest energy subhit index is "<<maxSubHitEnergyindex<<" E: "<<maxSubHitEnergy<<" identifier: "<<maxEnergyIdentifier.getString()<<G4endl; + } + //from identifier + CaloDetDescrElement *thiscell = m_calo_dd_man->get_element(id); + if (!maxEnergyCell) { + //How often does this happen? Do not shift. + G4cout << this->GetName()<<" WARNING ProcessHits: maxEnergyCell failed: "<<maxEnergyIdentifier.getString()<<G4endl + <<" "<<m_calo_dd_man->get_element(id)->getSampling()<<G4endl + <<" "<<originalStepPosition.eta()<<G4endl + <<" "<< originalStepPosition.phi()<<G4endl; + stepPosition = originalStepPosition; + } + else if (maxEnergyCell == thiscell) { + //The cells match, so do not shift this hit. + if(m_config.verboseLevel > 9) { + G4cout << this->GetName()<<" VERBOSE ProcessHits: Original step position: "<<originalStepPosition.x()<<" "<<originalStepPosition.y()<<" "<<originalStepPosition.z()<<G4endl + <<" "<<"This cell: "<<thiscell->x()<<" "<<thiscell->y()<<" "<<thiscell->z()<<G4endl + <<" "<<"No shift"<<G4endl; + } + stepPosition = originalStepPosition; + } + else { + //the two cells do not match => shift + G4ThreeVector diff(thiscell->x()-maxEnergyCell->x(), thiscell->y()-maxEnergyCell->y(), thiscell->z()-maxEnergyCell->z()); + stepPosition = originalStepPosition+diff; + if(m_config.verboseLevel > 9) { + CaloDetDescrElement *bestcell = m_calo_dd_man->get_element(m_calo_dd_man->get_element(id)->getSampling(),originalStepPosition.eta(), originalStepPosition.phi()); + G4cout << this->GetName()<<" VERBOSE ProcessHits: Original step position: "<<originalStepPosition.x()<<" "<<originalStepPosition.y()<<" "<<originalStepPosition.z()<<G4endl + <<" "<<"This cell: "<<thiscell->x()<<" "<<thiscell->y()<<" "<<thiscell->z()<<G4endl + <<" "<<"Highest E cell: "<<maxEnergyCell->x()<<" "<<maxEnergyCell->y()<<" "<<maxEnergyCell->z()<<G4endl + <<" "<<"(Best cell: "<<bestcell->x()<<" "<<bestcell->y()<<" "<<bestcell->z()<<")"<<G4endl + <<" "<<"Shifted step position: "<<stepPosition.x()<<" "<<stepPosition.y()<<" "<<stepPosition.z()<<G4endl; + } + } + } + } + } + //Finalize time for LAr hits?: NO + //double time = larhit.energy)==0 ? 0. : (double) larhit.time/larhit.energy/CLHEP::ns; + double time = larhit.time; + double energy = larhit.energy/CLHEP::MeV; + this->update_map(stepPosition, id, energy, time, true, numberOfProcessedHits); //store numberOfProcessedHits as info + }//numberOfProcessedHits + } //istep + if (madeSubSteps) { + //only delete steps when doing substeps. Do not delete the original G4Step! + while(!steps.empty()) { delete steps.back(); steps.pop_back(); } + } + } + return result; +} + +Identifier LArFCS_StepInfoSD::ConvertID(const LArG4Identifier& a_ident) const +{ + Identifier id; + if(a_ident[0]==4) { + // is LAr + if(a_ident[1]==1) { + //is LAr EM + try { + id = m_larEmID->channel_id(a_ident[2], // barrel_ec + a_ident[3], // sampling + a_ident[4], // region + a_ident[5], // eta + a_ident[6]); // phi + } + catch (LArID_Exception& e) { + G4ExceptionDescription description; + description << "ConvertID: LArEM_ID error code " << e.code() << " " + << (std::string) e; + G4Exception("LArFCS_StepInfoSD", "ConvertIDEM", FatalException, description); + abort(); + } + } + else if(a_ident[1]==2) { + //is EM HEC + try { + id = m_larHecID->channel_id(a_ident[2], // zSide + a_ident[3], // sampling + a_ident[4], // region + a_ident[5], // eta + a_ident[6]); // phi + } + catch(LArID_Exception& e) { + G4ExceptionDescription description; + description << "ConvertID: LArHEC_ID error code " << e.code() << " " + << (std::string) e; + G4Exception("LArFCS_StepInfoSD", "ConvertIDHEC", FatalException, description); + abort(); + } + } + else if(a_ident[1]==3) { + // FCAL + if(a_ident[3]>0) { + //is EM FCAL + try { + id = m_larFcalID->channel_id(a_ident[2], // zSide + a_ident[3], // sampling + a_ident[4], // eta + a_ident[5]); // phi + } + catch(LArID_Exception& e) { + G4ExceptionDescription description; + description << "ConvertID: LArFCAL_ID error code " << e.code() << " " + << (std::string) e; + G4Exception("LArFCS_StepInfoSD", "ConvertIDFCAL", FatalException, description); + abort(); + } + } + else { + //is Mini FCAL + try { + id = m_larMiniFcalID->channel_id(a_ident[2], // zSide + a_ident[3], // module + a_ident[4], // depth + a_ident[5], // eta + a_ident[6]); // phi + } + catch(LArID_Exception& e) { + G4ExceptionDescription description; + description << "ConvertID: LArMiniFCAL_ID error code " << e.code() << " " + << (std::string) e; + G4Exception("LArFCS_StepInfoSD", "ConvertIDMiniFCAL", FatalException, description); + abort(); + } + } + } + } + return id; +} diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/LArFCS_StepInfoSD.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/LArFCS_StepInfoSD.h new file mode 100644 index 0000000000000000000000000000000000000000..e5756d43a11aecd4e249103047392535d9ba3bcb --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/LArFCS_StepInfoSD.h @@ -0,0 +1,51 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef ISF_FASTCALOSIM_LARFCS_STEPINFOSD_H +#define ISF_FASTCALOSIM_LARFCS_STEPINFOSD_H + +// Base class header +#include "FCS_StepInfoSD.h" + +#include "LArG4Code/LArG4Identifier.h" + +// Forward declarations +class LArEM_ID; +class LArFCAL_ID; +class LArHEC_ID; +class LArMiniFCAL_ID; +class CaloDetDescrManager; + +class G4Step; +class G4TouchableHistory; +class ILArCalculatorSvc; + +/// @class FCS_StepInfoSD +/// @brief Common sensitive detector class for LAr systems. +/// +/// This SD implementation saves the standard LArHits. +/// See LArG4CalibSD for an SD that handles calibration hits. +/// +class LArFCS_StepInfoSD : public FCS_StepInfoSD +{ +public: + + /// Constructor + LArFCS_StepInfoSD(G4String a_name, const FCS_Param::Config& config); + + /// Destructor + virtual ~LArFCS_StepInfoSD(); + + /// Main processing method + G4bool ProcessHits(G4Step* a_step, G4TouchableHistory*) override; + +private: + /// Helper function for making "real" identifiers from LArG4Identifiers + Identifier ConvertID(const LArG4Identifier& a_ident) const; + + /// Member variable - the calculator we'll use + ILArCalculatorSvc * m_calculator; +}; + +#endif // ISF_FASTCALOSIM_LARFCS_STEPINFOSD_H diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/SDWrapper.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/SDWrapper.cxx new file mode 100644 index 0000000000000000000000000000000000000000..7df4b9c974dae8a427ff028af4465e9fbae0085c --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/SDWrapper.cxx @@ -0,0 +1,101 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "SDWrapper.h" + +// Geant4 includes +#include "G4SDManager.hh" + +// Framework utilities +#include "CxxUtils/make_unique.h" + +// Project includes +#include "ISF_FastCaloSimEvent/FCS_StepInfoCollection.h" + +// Local includes +#include "FCS_StepInfoSD.h" + +namespace FCS_Param +{ + + namespace detail + { + + //------------------------------------------------------------------------- + // Construct with a hit collection name + //------------------------------------------------------------------------- + template<class SDType, class HitContainerType> + SDWrapper<SDType, HitContainerType>:: + SDWrapper(const std::string& name, const std::string& hitCollectionName) + : G4VSensitiveDetector(name), + m_hitCollName(hitCollectionName), + m_hitColl(hitCollectionName) + {} + + //------------------------------------------------------------------------- + // Add a unique SD to the list + //------------------------------------------------------------------------- + template<class SDType, class HitContainerType> + void SDWrapper<SDType, HitContainerType>:: + addSD(std::unique_ptr<SDType> sd) + { + m_sdList.push_back( std::move(sd) ); + } + + //------------------------------------------------------------------------- + // Initialize the hit collection at the beginning of the G4 event + //------------------------------------------------------------------------- + template<class SDType, class HitContainerType> + void SDWrapper<SDType, HitContainerType>:: + Initialize(G4HCofThisEvent*) + { + if(!m_hitColl.isValid()) { + if(verboseLevel >= 5) { + G4cout << GetName() << " \tDEBUG\t" << "Initializing hit container: " + << m_hitCollName << G4endl; + } + m_hitColl = CxxUtils::make_unique<HitContainerType>(); + } + } + + //------------------------------------------------------------------------- + // This should not be called + //------------------------------------------------------------------------- + template<class SDType, class HitContainerType> + bool SDWrapper<SDType, HitContainerType>:: + ProcessHits(G4Step*, G4TouchableHistory*) + { + G4ExceptionDescription description; + description << "ProcessHits: this SD shouldn't be assigned to volumes!"; + G4Exception(GetName(), "SDError", FatalException, description); + return false; + } + + //------------------------------------------------------------------------- + // Gather the hits into the WriteHandle from all the SDs + //------------------------------------------------------------------------- + template<class SDType, class HitContainerType> + void SDWrapper<SDType, HitContainerType>:: + EndOfAthenaEvent() + { + if(!m_hitColl.isValid()) { + G4cerr << GetName() << " \tERROR\t" << "Hit collection WriteHandle is " + << "invalid!" << G4endl; + throw std::runtime_error("Invalid hit container WriteHandle: " + + m_hitColl.name()); + } + // Loop over each SD and fill the container + for(auto& sd : m_sdList) { + sd->EndOfAthenaEvent( &*m_hitColl ); + } + } + + //------------------------------------------------------------------------- + // Explit template instantiations + //------------------------------------------------------------------------- + template class SDWrapper<FCS_StepInfoSD, ISF_FCS_Parametrization::FCS_StepInfoCollection>; + + } // namespace detail + +} // namespace FCS_Param diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/SDWrapper.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/SDWrapper.h new file mode 100644 index 0000000000000000000000000000000000000000..ef3a248cf715e95334a6e6befe294f3b84f5c72e --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/SDWrapper.h @@ -0,0 +1,94 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef ISF_FASTCALOSIMSD_SDWRAPPER_H +#define ISF_FASTCALOSIMSD_SDWRAPPER_H + +// System includes +#include <string> +#include <vector> +#include <memory> + +// External includes +#include "G4VSensitiveDetector.hh" + +// Framework includes +#include "StoreGate/WriteHandle.h" + + +// Forward declarations +class FCS_StepInfoSD; +namespace ISF_FCS_Parametrization { + class FCS_StepInfoCollection; +} +namespace FCS_Param +{ + + namespace detail + { + + /// @class SDWrapper + /// @brief A template class which wraps multiple sensitive detectors. + /// + /// Allows for SD tools to manage several SDs which collaborate to fill one + /// hit container in a multi-threading-friendly way. The wrapper owns the + /// WriteHandle for the hit container and gathers hits from each SD at the + /// end of an event. + /// + /// The inheritance from G4VSensitiveDetector is merely a trick so the SD + /// tool can save this object in the SensitiveDetectorBase thread-local + /// container. It also allows to create the hit container at the right time + /// via the SD Initialize method invoked by Geant4. + /// + /// Clients shouldn't use this generic template directly, but should use + /// the explicitly allowed specializations given below. + /// + /// @author Steve Farrell <Steven.Farrell@cern.ch> + /// + template<class SDType, class HitContainerType> + class SDWrapper : public G4VSensitiveDetector + { + + public: + + /// Alias to the SD list type + using SDList_t = std::vector< std::unique_ptr<SDType> >; + + /// Construct the wrapper from the output collection name + SDWrapper(const std::string& name, const std::string& hitCollectionName); + + /// Add an SD to this wrapper + void addSD(std::unique_ptr<SDType> sd); + + /// Beginning of G4 event; initialize the hit collection. + virtual void Initialize(G4HCofThisEvent*) override final; + + /// This method should not be called. It will throw. + virtual bool ProcessHits(G4Step*, G4TouchableHistory*) override final; + + /// Gather the hits into the WriteHandle from all the SDs + void EndOfAthenaEvent(); + + private: + + /// The hit container name + std::string m_hitCollName; + + /// The hit container handle + SG::WriteHandle<HitContainerType> m_hitColl; + + /// The list of sensitive detectors that I own and manage + SDList_t m_sdList; + + }; // class SDWrapper + + } // namespace detail + + + /// Template instantiation for LArG4SimpleSD + using FCS_StepInfoSDWrapper = detail::SDWrapper<FCS_StepInfoSD, ISF_FCS_Parametrization::FCS_StepInfoCollection>; + +} // namespace FCS_Param + +#endif diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/TileFCS_StepInfoSD.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/TileFCS_StepInfoSD.cxx new file mode 100644 index 0000000000000000000000000000000000000000..41ca1bd2cbc1ae42d71150395b899620c0ae4668 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/TileFCS_StepInfoSD.cxx @@ -0,0 +1,100 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +//************************************************************ +// +// Class TileFCS_StepInfoSD +// Sensitive detector for TileCal G4 simulations with TileGeoModel +// +// Author: Vakho Tsulaia <Vakhtang.Tsulaia@cern.ch> +// +// Major updates: July, 2013 (Sergey) +// +//************************************************************ + +// class header +#include "TileFCS_StepInfoSD.h" +/// Geant4 headers +#include "G4Step.hh" +#include "G4ThreeVector.hh" +#include "G4TouchableHistory.hh" + +/// Athena headers +#include "CaloDetDescr/CaloDetDescrManager.h" +#include "TileG4Interfaces/ITileCalculator.h" + +TileFCS_StepInfoSD::TileFCS_StepInfoSD(G4String name, const FCS_Param::Config & config) + : FCS_StepInfoSD(name, config) + , m_calculator(config.m_TileCalculator) +{ +} + +TileFCS_StepInfoSD::~TileFCS_StepInfoSD() { +} + +G4bool TileFCS_StepInfoSD::ProcessHits(G4Step* a_step, G4TouchableHistory* /*ROhist*/) { + G4bool result(false); + // If there's no energy, there's no hit. (Aside: Isn't this energy + // the same as the energy from the calculator? Not necessarily. + // The calculator may include detector effects such as + // charge-collection which are not modeled by Geant4.) + if(a_step->GetTotalEnergyDeposit() <= 0.) { return result; } + + if (m_calculator) { + //calculation of MicroHit with a_step + TileHitData hitData; + TileMicroHit micHit = m_calculator->GetTileMicroHit(a_step, hitData); + Identifier m_invalid_id; + + //Check if MicroHit is not in scintillator + if ((micHit.pmt_up == m_invalid_id) && (micHit.pmt_down == m_invalid_id)) { + G4cout <<this->GetName()<<" WARNING ProcessHits: Invalid hit in Tile??"<<G4endl; + return result; + } + else { + // Store TileHits Information + if ((micHit.pmt_up == m_invalid_id) || (micHit.pmt_down == m_invalid_id)) { + G4cout <<this->GetName()<<" WARNING ProcessHits: Something wrong in identifier: One tile pmt: "<<micHit.pmt_up<<" "<<micHit.pmt_down<<std::endl; + G4cout <<this->GetName()<<" WARNING ProcessHits: E up: "<<micHit.e_up<<" E down: "<<micHit.e_down<<" T up: "<<micHit.time_up<<" T down: "<<micHit.time_down<<std::endl; + } + const G4ThreeVector pos = 0.5*(a_step->GetPreStepPoint()->GetPosition()+a_step->GetPostStepPoint()->GetPosition()); + + this->update_map(pos, micHit.pmt_up, micHit.e_up, micHit.time_up, true,1); + this->update_map(pos, micHit.pmt_down, micHit.e_down,micHit.time_down , true,1); + } + } + return true; +} + +void TileFCS_StepInfoSD::update_map(const CLHEP::Hep3Vector & l_vec, const Identifier & l_cell, double l_energy, double l_time, bool l_valid, int l_detector) +{ + // Drop any hits that don't have a good identifier attached + if (!m_calo_dd_man->get_element(l_cell)) { return; } + + auto map_item = m_hit_map.find( l_cell ); + if (map_item==m_hit_map.end()) { + m_hit_map[l_cell] = new std::vector< ISF_FCS_Parametrization::FCS_StepInfo* >; + m_hit_map[l_cell]->push_back( new ISF_FCS_Parametrization::FCS_StepInfo( l_vec , l_cell , l_energy , l_time , l_valid , l_detector ) ); + } + else { + bool match = false; + for (auto map_it : * map_item->second) { + // Time check + const double delta_t = std::fabs(map_it->time()-l_time); + if ( delta_t >= m_config.m_maxTimeTile ) { continue; } + // Distance check + const double hit_diff2 = map_it->position().diff2( l_vec ); + if ( hit_diff2 >= m_config.m_maxRadiusTile ) { continue; } + // Found a match. Make a temporary that will be deleted! + const ISF_FCS_Parametrization::FCS_StepInfo my_info( l_vec , l_cell , l_energy , l_time , l_valid , l_detector ); + *map_it += my_info; + match = true; + break; + } // End of search for match in time and space + if (!match) { + map_item->second->push_back( new ISF_FCS_Parametrization::FCS_StepInfo( l_vec , l_cell , l_energy , l_time , l_valid , l_detector ) ); + } // Didn't match + } // ID already in the map + +} // That's it for updating the map! diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/TileFCS_StepInfoSD.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/TileFCS_StepInfoSD.h new file mode 100644 index 0000000000000000000000000000000000000000..7664e789390f3e49f8965800974f3357fb5a5fac --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/TileFCS_StepInfoSD.h @@ -0,0 +1,48 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +//************************************************************ +// +// Class TileFCS_StepInfoSD +// Sensitive detector for TileCal G4 simulations using TileGeoModel +// +// Author: Vakho Tsulaia <Vakhtang.Tsulaia@cern.ch> +// +// Major updates: July, 2013 (Sergey) +// +//************************************************************ + +#ifndef ISF_FASTCALOSIMSD_TILEFCS_STEPINFOSD_H +#define ISF_FASTCALOSIMSD_TILEFCS_STEPINFOSD_H + +// Base class header +#include "FCS_StepInfoSD.h" + +class G4Step; +class G4TouchableHistory; +class ITileCalculator; + +class TileFCS_StepInfoSD: public FCS_StepInfoSD { +public: + TileFCS_StepInfoSD(G4String name, const FCS_Param::Config& config); + ~TileFCS_StepInfoSD(); + + G4bool ProcessHits(G4Step*, G4TouchableHistory*) override final; + + // Don't leave copy constructors or assignment operators around + TileFCS_StepInfoSD(const TileFCS_StepInfoSD&) = delete; + TileFCS_StepInfoSD& operator=(const TileFCS_StepInfoSD&) = delete; + +protected: + /// Keep a map instead of trying to keep the full vector. + /// At the end of the event we'll push the map back into the + /// FCS_StepInfoCollection in StoreGate. + virtual void update_map(const CLHEP::Hep3Vector & l_vec, const Identifier & l_cell, double l_energy, double l_time, bool l_valid, int l_detector) override final; + +private: + ITileCalculator* m_calculator; +}; + +#endif // ISF_FASTCALOSIMSD_TILEFCS_STEPINFOSD_H + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/components/ISF_FastCaloSimSD_entries.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/components/ISF_FastCaloSimSD_entries.cxx new file mode 100644 index 0000000000000000000000000000000000000000..0da3aea99724c9970f7a4c6e3a5d20c0e996439a --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimSD/src/components/ISF_FastCaloSimSD_entries.cxx @@ -0,0 +1,5 @@ +#include "GaudiKernel/DeclareFactoryEntries.h" + +#include "../FCS_StepInfoSDTool.h" + +DECLARE_TOOL_FACTORY( FCS_Param::FCS_StepInfoSDTool )