diff --git a/Simulation/G4Utilities/G4DebuggingTools/python/G4DebuggingToolsConfig.py b/Simulation/G4Utilities/G4DebuggingTools/python/G4DebuggingToolsConfig.py index f52b57fe2059fda006a2cf456c9959bd184320f5..a57217a0ede98ecc27ebb47913d4e6311ccff76e 100644 --- a/Simulation/G4Utilities/G4DebuggingTools/python/G4DebuggingToolsConfig.py +++ b/Simulation/G4Utilities/G4DebuggingTools/python/G4DebuggingToolsConfig.py @@ -55,6 +55,18 @@ def getStepNtupleTool(name="G4UA::StepNtupleTool", **kwargs): return False return CfgMgr.G4UA__StepNtupleTool(name, **kwargs) +def getStepHistogramTool(name="G4UA::StepHistogramTool", **kwargs): + from AthenaCommon.ConcurrencyFlags import jobproperties as concurrencyProps + if concurrencyProps.ConcurrencyFlags.NumThreads() >1: + log=Logging.logging.getLogger(name) + log.fatal('Attempt to run '+name+' with more than one thread, which is not supported') + return False + from G4AtlasApps.SimFlags import simFlags + 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__StepHistogramTool(name, **kwargs) + def getVolumeDebuggerTool(name="G4UA::VolumeDebuggerTool", **kwargs): from AthenaCommon.ConcurrencyFlags import jobproperties as concurrencyProps diff --git a/Simulation/G4Utilities/G4DebuggingTools/python/G4DebuggingToolsConfigDb.py b/Simulation/G4Utilities/G4DebuggingTools/python/G4DebuggingToolsConfigDb.py index 106f99057e4740e1e72832517f5649db553133fe..6e3bd096942d6207fd85d92d7de8963d1d424bd7 100644 --- a/Simulation/G4Utilities/G4DebuggingTools/python/G4DebuggingToolsConfigDb.py +++ b/Simulation/G4Utilities/G4DebuggingTools/python/G4DebuggingToolsConfigDb.py @@ -6,6 +6,7 @@ addTool("G4DebuggingTools.G4DebuggingToolsConfig.getVolumeDebuggerTool", "G4UA:: addTool("G4DebuggingTools.G4DebuggingToolsConfig.getG4AtlantisDumperTool", "G4UA::G4AtlantisDumperTool") addTool("G4DebuggingTools.G4DebuggingToolsConfig.getVerboseSelectorTool", "G4UA::VerboseSelectorTool") addTool("G4DebuggingTools.G4DebuggingToolsConfig.getStepNtupleTool", "G4UA::StepNtupleTool") +addTool("G4DebuggingTools.G4DebuggingToolsConfig.getStepHistogramTool", "G4UA::StepHistogramTool") addTool("G4DebuggingTools.G4DebuggingToolsConfig.getEnergyConservationTestTool", "G4UA::EnergyConservationTestTool") addTool("G4DebuggingTools.G4DebuggingToolsConfig.getHyperspaceCatcherTool", "G4UA::HyperspaceCatcherTool") addTool("G4DebuggingTools.G4DebuggingToolsConf.G4UA__CheckActivationTool", "G4UA::CheckActivationTool") diff --git a/Simulation/G4Utilities/G4DebuggingTools/share/preInclude.StepHistogram.py b/Simulation/G4Utilities/G4DebuggingTools/share/preInclude.StepHistogram.py new file mode 100644 index 0000000000000000000000000000000000000000..9bebc8bb5a5356531a43832f9f411bbc28ab294d --- /dev/null +++ b/Simulation/G4Utilities/G4DebuggingTools/share/preInclude.StepHistogram.py @@ -0,0 +1,13 @@ +## Histogram Service +from AthenaCommon.AppMgr import ServiceMgr +if not hasattr(ServiceMgr, 'THistSvc'): + from GaudiSvc.GaudiSvcConf import THistSvc + ServiceMgr += THistSvc() +ServiceMgr.THistSvc.Output = ["stepHisto DATAFILE='StepHistograms.root' OPT='RECREATE'"]; +print ServiceMgr.THistSvc + +## User Actions +from G4AtlasApps.SimFlags import simFlags +simFlags.OptionalUserActionList.addAction('G4UA::StepHistogramTool', ['Step']) +simFlags.UserActionConfig.addConfig('G4UA::StepHistogramTool', 'doGeneralHistograms', True) +simFlags.UserActionConfig.addConfig('G4UA::StepHistogramTool', 'do2DHistograms', False) diff --git a/Simulation/G4Utilities/G4DebuggingTools/src/G4DebuggingHelper.cxx b/Simulation/G4Utilities/G4DebuggingTools/src/G4DebuggingHelper.cxx new file mode 100644 index 0000000000000000000000000000000000000000..f6568e139876e83b7e6fc7d231861f592476d442 --- /dev/null +++ b/Simulation/G4Utilities/G4DebuggingTools/src/G4DebuggingHelper.cxx @@ -0,0 +1,121 @@ +#include "StepHistogram.h" + +#include "G4Electron.hh" +#include "G4Positron.hh" +#include "G4Gamma.hh" +#include "G4Neutron.hh" +#include "G4Proton.hh" + +namespace G4DebuggingHelpers { + + const G4String ClassifyParticle( const G4ParticleDefinition* def ) { + if (def == G4Electron::Electron()) + return "e-"; + else if (def == G4Positron::Positron()) + return "e+"; + else if (def == G4Gamma::Gamma()) + return "gamma"; + else if (def == G4Neutron::Neutron()) + return "neutron"; + else if (def == G4Proton::Proton()) + return "proton"; + return "other"; + } + + const G4String ClassifyMaterial( const G4String &nom ) { + if (nom == "FCal1Absorber" + || nom == "LiquidArgon" + || nom == "Copper" + || nom == "Lead" + || nom == "Aluminum" + || nom == "FCal23Absorber" + || nom == "Iron" + || nom == "Air" + || nom == "myLead" + || nom == "shieldIron" + || nom == "FCal23Slugs" + || nom == "Glue" + || nom == "KaptonC" + || nom == "Kapton" + || nom == "ShieldSteel" + || nom == "myIron" + || nom == "ShieldBrass" + || nom == "Straw" + || nom == "XeCO2O2" + || nom == "CO2" + || nom == "Valmat" + || nom == "BoratedPolyethelyne" + || nom == "FoilRadiatorB" + || nom == "G10" + || nom == "FoilRadiatorAC" + || nom == "PyrogelXT" + || nom == "Vacuum") + return nom; + else if (nom.substr(0,12)=="pix::IBL_Fwd") + return "IBL_Fwd"; + return "other"; + } + + const G4String ClassifyVolume( const G4String &nom ) { + if ( nom.length() >= 17 && nom.substr(13, 4) == "EMEC" ) { + return "EMEC"; + } + else if ( nom.length() >= 16 && nom.substr(13, 3) == "EMB" ) { + return "EMB"; + } + else if ( nom.length() >= 25 && nom.substr(21, 4) == "Cryo" ) { + return "Cryo"; + } + else if ( nom.length() >= 26 && nom.substr(13, 13) == "FCAL::Module1" ) { + return "FC1"; + } + else if ( nom.length() >= 25 && nom.substr(13, 12) == "FCAL::Module" ) { + return "FC23"; + } + else if ( nom.length() >= 17 && nom.substr(13, 4) == "FCAL" ) { + return "FCOther"; + } + else if ( nom.length() >= 16 && nom.substr(13, 3) == "HEC" ) { + return "HEC"; + } + else if ( nom.length() >= 31 && nom.substr(21, 10) == "Presampler" ) { + return "Presampler"; + } + else if ( nom.length() >= 3 && nom.substr(0, 3) == "LAr" ) { + return "LAr"; + } + else if ( ( nom.substr(0, 4) == "MUON" ) + || ( nom.length() >= 4 && nom.substr(0, 4) == "Muon" ) + || ( nom.length() >= 9 && nom.substr(0, 9) == "DriftTube" ) + || ( nom.length() >= 12 && nom.substr(0, 12) == "SensitiveGas" ) + || nom.contains("MDT") + || nom.contains("station") ) { + return "Muon"; + } + else if ( ( nom.length() >= 5 && nom.substr(0, 5) == "Pixel" ) + || nom == "Outside Barrel Service") { + return "Pixel"; + } + else if ( nom.length() >= 3 && nom.substr(0, 3) == "SCT" ) { + return "SCT"; + } + else if ( ( nom.length() >= 3 && nom.substr(0, 3) == "TRT" ) + || nom == "GasMANeg" ) { + return "TRT"; + } + else if ( nom.length() >= 4 && nom.substr(0, 4) == "Tile" ) { + return "Tile"; + } + else if (nom.length() >= 7 && nom.substr(0, 7) == "Section" ) + return "Section"; + else if ( ( nom.length() >= 12 && nom.substr(0, 12) == "InDetServMat" ) + || ( nom.length() >= 4 && nom.substr(0, 4) == "IDET" ) + || ( nom.length() >= 8 && nom.substr(0, 8) == "BeamPipe" ) + || ( nom.length() >= 3 && + ( nom.substr(0, 3) == "BLM" || nom.substr(0, 3) == "BCM" ) ) ) { + return "Service"; + } + return "other"; + } + +} // end namespace G4DebuggingHelpers diff --git a/Simulation/G4Utilities/G4DebuggingTools/src/G4DebuggingHelper.h b/Simulation/G4Utilities/G4DebuggingTools/src/G4DebuggingHelper.h new file mode 100644 index 0000000000000000000000000000000000000000..272cddfa61f406131b6bd194a00e3be5235443a0 --- /dev/null +++ b/Simulation/G4Utilities/G4DebuggingTools/src/G4DebuggingHelper.h @@ -0,0 +1,19 @@ +#ifndef G4DEBUGGINGTOOLS_G4DEBUGGINGHELPER_H +#define G4DEBUGGINGTOOLS_G4DEBUGGINGHELPER_H + +#include "G4String.hh" +#include "G4ParticleDefinition.hh" + +#include <map> + +namespace G4DebuggingHelpers { + + const G4String ClassifyParticle( const G4ParticleDefinition* def ); + + const G4String ClassifyMaterial( const G4String &nom ); + + const G4String ClassifyVolume( const G4String &nom ); + +} + +#endif \ No newline at end of file diff --git a/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogram.cxx b/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogram.cxx new file mode 100644 index 0000000000000000000000000000000000000000..8c603afaccf2a4d0f95123a1cb27215120fdace9 --- /dev/null +++ b/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogram.cxx @@ -0,0 +1,261 @@ +/* + Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +*/ + +#include "StepHistogram.h" + +#include "GaudiKernel/IDataProviderSvc.h" +#include "GaudiKernel/INTupleSvc.h" +#include "GaudiKernel/NTuple.h" +#include "GaudiKernel/SmartDataPtr.h" +#include "GaudiKernel/Bootstrap.h" +#include "GaudiKernel/ISvcLocator.h" +#include "GaudiKernel/IMessageSvc.h" + +#include "SimHelpers/ServiceAccessor.h" + +#include "G4Step.hh" +#include "G4VProcess.hh" + +#include <iostream> + +namespace G4UA{ + + StepHistogram::StepHistogram(const Config& config): + AthMessaging(Gaudi::svcLocator()->service<IMessageSvc>("MessageSvc"),"StepHistogram"), + m_config(config), + m_initialKineticEnergyOfStep() + {} + + void StepHistogram::UserSteppingAction(const G4Step * aStep){ + + G4Track *tr = aStep->GetTrack(); + G4ThreeVector myPos = aStep->GetPostStepPoint()->GetPosition(); + + G4String particleName = "nucleus"; + if (!(tr->GetDefinition()->GetParticleType() == "nucleus")) + particleName = G4DebuggingHelpers::ClassifyParticle(tr->GetParticleDefinition()); + + G4String volumeName = tr->GetVolume()->GetName(); + volumeName = G4DebuggingHelpers::ClassifyVolume(volumeName); + + G4String materialName = tr->GetMaterial()->GetName(); + materialName = G4DebuggingHelpers::ClassifyMaterial(materialName); + + G4String processName = aStep->GetPostStepPoint()->GetProcessDefinedStep() ? + aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() : "Unknown"; + + const std::vector<const G4Track*>* secondaries = aStep->GetSecondaryInCurrentStep(); + + // 2D map + if (m_config.do2DHistograms) { + InitializeFillHistogram2D(m_report.histoMapMap2D_vol_RZ, "vol_RZ", particleName, volumeName, + 1000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), 1.); + InitializeFillHistogram2D(m_report.histoMapMap2D_mat_RZ, "mat_RZ", particleName, materialName, + 1000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), 1.); + InitializeFillHistogram2D(m_report.histoMapMap2D_prc_RZ, "prc_RZ", particleName, processName, + 1000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), 1.); + } + + // step length + InitializeFillHistogram(m_report.histoMapMap_vol_stepSize, "vol_stepLength", particleName, volumeName, + 1000, -12, 4, log10(aStep->GetStepLength()), 1.); + InitializeFillHistogram(m_report.histoMapMap_mat_stepSize, "mat_stepLength", particleName, materialName, + 1000, -12, 4, log10(aStep->GetStepLength()), 1.); + InitializeFillHistogram(m_report.histoMapMap_prc_stepSize, "prc_stepLength", particleName, processName, + 1000, -12, 4, log10(aStep->GetStepLength()), 1.); + + // step pseudorapidity + InitializeFillHistogram(m_report.histoMapMap_vol_stepPseudorapidity, "vol_stepPseudorapidity", particleName, volumeName, + 200, -10, 10, myPos.eta(), 1.); + InitializeFillHistogram(m_report.histoMapMap_mat_stepPseudorapidity, "mat_stepPseudorapidity", particleName, materialName, + 200, -10, 10, myPos.eta(), 1.); + InitializeFillHistogram(m_report.histoMapMap_prc_stepPseudorapidity, "prc_stepPseudorapidity", particleName, processName, + 200, -10, 10, myPos.eta(), 1.); + + // step kinetic energy + InitializeFillHistogram(m_report.histoMapMap_vol_stepKineticEnergy, "vol_stepKineticEnergy", particleName, volumeName, + 1000, -9, 7, log10(tr->GetKineticEnergy()), 1.); + InitializeFillHistogram(m_report.histoMapMap_mat_stepKineticEnergy, "mat_stepKineticEnergy", particleName, materialName, + 1000, -9, 7, log10(tr->GetKineticEnergy()), 1.); + InitializeFillHistogram(m_report.histoMapMap_prc_stepKineticEnergy, "prc_stepKineticEnergy", particleName, processName, + 1000, -9, 7, log10(tr->GetKineticEnergy()), 1.); + InitializeFillHistogram(m_report.histoMapMap_stepSecondaryKinetic, "stepKineticEnergy", particleName, "AllATLAS", + 1000, -9, 7, log10(tr->GetKineticEnergy()), 1.); + + // step energy deposit + InitializeFillHistogram(m_report.histoMapMap_vol_stepEnergyDeposit, "vol_stepEnergyDeposit", particleName, volumeName, + 1000, -11, 3, log10(aStep->GetTotalEnergyDeposit()), 1.); + InitializeFillHistogram(m_report.histoMapMap_mat_stepEnergyDeposit, "mat_stepEnergyDeposit", particleName, materialName, + 1000, -11, 3, log10(aStep->GetTotalEnergyDeposit()), 1.); + InitializeFillHistogram(m_report.histoMapMap_prc_stepEnergyDeposit, "prc_stepEnergyDeposit", particleName, processName, + 1000, -11, 3, log10(aStep->GetTotalEnergyDeposit()), 1.); + + // step non-ionizing energy deposit + InitializeFillHistogram(m_report.histoMapMap_vol_stepEnergyNonIonDeposit, "vol_stepEnergyNonIonDeposit", particleName, volumeName, + 1000, -11, 1, log10(aStep->GetNonIonizingEnergyDeposit()), 1.); + InitializeFillHistogram(m_report.histoMapMap_mat_stepEnergyNonIonDeposit, "mat_stepEnergyNonIonDeposit", particleName, materialName, + 1000, -11, 1, log10(aStep->GetNonIonizingEnergyDeposit()), 1.); + InitializeFillHistogram(m_report.histoMapMap_prc_stepEnergyNonIonDeposit, "prc_stepEnergyNonIonDeposit", particleName, processName, + 1000, -11, 1, log10(aStep->GetNonIonizingEnergyDeposit()), 1.); + + // secondary kinetic energy + for (const auto &track : *secondaries) { + G4String secondary_particleName = G4DebuggingHelpers::ClassifyParticle(track->GetParticleDefinition()); + InitializeFillHistogram(m_report.histoMapMap_vol_stepSecondaryKinetic, "vol_stepSecondaryKinetic", secondary_particleName, volumeName, + 1000, -7, 5, log10(track->GetKineticEnergy()), 1.); + InitializeFillHistogram(m_report.histoMapMap_mat_stepSecondaryKinetic, "mat_stepSecondaryKinetic", secondary_particleName, materialName, + 1000, -7, 5, log10(track->GetKineticEnergy()), 1.); + InitializeFillHistogram(m_report.histoMapMap_prc_stepSecondaryKinetic, "prc_stepSecondaryKinetic", secondary_particleName, processName, + 1000, -7, 5, log10(track->GetKineticEnergy()), 1.); + } + + // stop here if 'general' histograms not activated + // _______________________________________________ + if (!m_config.doGeneralHistograms) + return; + + // first step (after initial step) + if (tr->GetCurrentStepNumber()==1) { + m_initialKineticEnergyOfStep = tr->GetKineticEnergy(); + m_initialKineticEnergyOfStep += aStep->GetTotalEnergyDeposit(); + for (const auto &track : *secondaries) { + m_initialKineticEnergyOfStep += track->GetKineticEnergy(); + } + // save track ID for checking if we later have the same track + m_trackID = tr->GetTrackID(); + // initial energy + InitializeFillHistogram(m_report.histoMapMap_InitialE, "InitialE", particleName, "AllATLAS", + 1000, -9, 7, log10(m_initialKineticEnergyOfStep), 1.0); + } + + // last step + if ( tr->GetTrackStatus() == 2 ) { + // assert to check if we have the correct track + if (not (tr->GetTrackID() == m_trackID)) { + ATH_MSG_ERROR("Track ID changed between the assumed first step and the last."); + throw std::exception(); + } + // number of steps + int nSteps = tr->GetCurrentStepNumber() + 1; + InitializeFillHistogram(m_report.histoMapMap_numberOfSteps, "numberOfSteps", particleName, "AllATLAS", + 10000, 0.5, 10000.5, nSteps, 1.); + // number of steps vs initial energy + InitializeFillHistogram(m_report.histoMapMap_numberOfStepsPerInitialE, "numberOfStepsPerInitialE", particleName, "AllATLAS", + 1000, -9, 7, log10(m_initialKineticEnergyOfStep), nSteps); + } + } + + void StepHistogram::InitializeFillHistogram2D(HistoMapMap_t &hMapMap, const char* suffix, + G4String particleName, G4String vol, + int nbinsx, double xmin, double xmax, + int nbinsy, double ymin, double ymax, + double valuex, double valuey, double weight) + { + if ( hMapMap.find(vol) == hMapMap.end() ) { + // initialize HistoMap_t if not yet exist + hMapMap.emplace(vol,HistoMap_t()); + } + HistoMap_t &hMap = hMapMap[vol]; + if ( hMap.find(particleName) == hMap.end() ) { + // initialize histogram if not yet exist + std::ostringstream stringStream; + stringStream << vol << "_" << particleName << "_" << suffix; + hMap[particleName] = new TH2F(stringStream.str().c_str(), stringStream.str().c_str(), nbinsx, xmin, xmax, nbinsy, ymin, ymax); + } + ((TH2*)hMap[particleName])->Fill(valuex, valuey, weight); + } + + void StepHistogram::InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix, + G4String particleName, G4String vol, + int nbins, double xmin, double xmax, double value, double weight) + { + if ( hMapMap.find(vol) == hMapMap.end() ) { + // initialize HistoMap_t if not yet exist + hMapMap.emplace(vol,HistoMap_t()); + } + HistoMap_t &hMap = hMapMap[vol]; + if ( hMap.find(particleName) == hMap.end() ) { + // initialize histogram if not yet exist + std::ostringstream stringStream; + stringStream << vol << "_" << particleName << "_" << suffix; + hMap[particleName] = new TH1F(stringStream.str().c_str(), stringStream.str().c_str(), nbins, xmin, xmax); + } + hMap[particleName]->Fill(value, weight); + } + + void StepHistogram::InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix, + G4String particleName, G4String vol, + int nbins, double *edges, double value, double weight) + { + if ( hMapMap.find(vol) == hMapMap.end() ) { + // initialize HistoMap_t if not yet exist + hMapMap.emplace(vol,HistoMap_t()); + } + HistoMap_t &hMap = hMapMap[vol]; + if ( hMap.find(particleName) == hMap.end() ) { + // initialize histogram if not yet exist + std::ostringstream stringStream; + stringStream << vol << "_" << particleName << "_" << suffix; + hMap[particleName] = new TH1F(stringStream.str().c_str(), stringStream.str().c_str(), nbins, edges); + } + hMap[particleName]->Fill(value, weight); + } + + void StepHistogram::Report::mergeMaps(HistoMapMap_t &selfMap, HistoMapMap_t refMap) { + for (auto const& ref : refMap) + { + if ( selfMap.find(ref.first) == selfMap.end() ) { + // HistoMap_t does not yet exist + selfMap.emplace(ref.first, ref.second); + } + else { + HistoMap_t &target = selfMap[ref.first]; + for (auto const& hm : ref.second) + { + if ( target.find(hm.first) == target.end() ) { + // histogram does not yet exist + target.emplace(hm.first, hm.second); + } + else { + // add histograms + target[hm.first]->Add(hm.second); + } + } + } + } + } + + void StepHistogram::Report::merge(const Report & rep) { + mergeMaps(histoMapMap_vol_stepSize, rep.histoMapMap_vol_stepSize); + mergeMaps(histoMapMap_vol_stepKineticEnergy, rep.histoMapMap_vol_stepKineticEnergy); + mergeMaps(histoMapMap_vol_stepPseudorapidity, rep.histoMapMap_vol_stepPseudorapidity); + mergeMaps(histoMapMap_vol_stepEnergyDeposit, rep.histoMapMap_vol_stepEnergyDeposit); + mergeMaps(histoMapMap_vol_stepEnergyNonIonDeposit, rep.histoMapMap_vol_stepEnergyNonIonDeposit); + mergeMaps(histoMapMap_vol_stepSecondaryKinetic, rep.histoMapMap_vol_stepSecondaryKinetic); + + mergeMaps(histoMapMap_mat_stepSize, rep.histoMapMap_mat_stepSize); + mergeMaps(histoMapMap_mat_stepKineticEnergy, rep.histoMapMap_mat_stepKineticEnergy); + mergeMaps(histoMapMap_mat_stepPseudorapidity, rep.histoMapMap_mat_stepPseudorapidity); + mergeMaps(histoMapMap_mat_stepEnergyDeposit, rep.histoMapMap_mat_stepEnergyDeposit); + mergeMaps(histoMapMap_mat_stepEnergyNonIonDeposit, rep.histoMapMap_mat_stepEnergyNonIonDeposit); + mergeMaps(histoMapMap_mat_stepSecondaryKinetic, rep.histoMapMap_mat_stepSecondaryKinetic); + + mergeMaps(histoMapMap_prc_stepSize, rep.histoMapMap_prc_stepSize); + mergeMaps(histoMapMap_prc_stepKineticEnergy, rep.histoMapMap_prc_stepKineticEnergy); + mergeMaps(histoMapMap_prc_stepPseudorapidity, rep.histoMapMap_prc_stepPseudorapidity); + mergeMaps(histoMapMap_prc_stepEnergyDeposit, rep.histoMapMap_prc_stepEnergyDeposit); + mergeMaps(histoMapMap_prc_stepEnergyNonIonDeposit, rep.histoMapMap_prc_stepEnergyNonIonDeposit); + mergeMaps(histoMapMap_prc_stepSecondaryKinetic, rep.histoMapMap_prc_stepSecondaryKinetic); + + mergeMaps(histoMapMap_numberOfSteps, rep.histoMapMap_numberOfSteps); + mergeMaps(histoMapMap_numberOfStepsPerInitialE, rep.histoMapMap_numberOfStepsPerInitialE); + mergeMaps(histoMapMap_InitialE, rep.histoMapMap_InitialE); + mergeMaps(histoMapMap_stepSecondaryKinetic, rep.histoMapMap_stepSecondaryKinetic); + + mergeMaps(histoMapMap2D_vol_RZ, rep.histoMapMap2D_vol_RZ); + mergeMaps(histoMapMap2D_mat_RZ, rep.histoMapMap2D_mat_RZ); + mergeMaps(histoMapMap2D_prc_RZ, rep.histoMapMap2D_prc_RZ); + } + +} // namespace G4UA diff --git a/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogram.h b/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogram.h new file mode 100644 index 0000000000000000000000000000000000000000..9414b971e623e33b641dfa0739f0e9465db37786 --- /dev/null +++ b/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogram.h @@ -0,0 +1,127 @@ +/* + Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef G4DEBUGGINGTOOLS_StepHistogram_H +#define G4DEBUGGINGTOOLS_StepHistogram_H + +//C++ +#include <map> +#include <string> + +//ROOT +#include "TH1.h" +#include "TH2.h" + +//G4 +#include "G4UserEventAction.hh" +#include "G4UserSteppingAction.hh" +#include "G4UserRunAction.hh" +#include "G4String.hh" + +//Athena +#include "AthenaBaseComps/AthMessaging.h" +#include "G4DebuggingHelper.h" + +namespace G4UA{ + + class StepHistogram : public AthMessaging, + public G4UserEventAction, + public G4UserRunAction, + public G4UserSteppingAction + { + public: + /// the hooks for G4 UA handling + virtual void UserSteppingAction(const G4Step*) override; + + // maps to hold info per volume/process/material per particle type + typedef std::map<G4String, TH1*> HistoMap_t; + typedef std::map<G4String, HistoMap_t> HistoMapMap_t; + + /// this holds all the data from individual threads that needs to be merged at EoR + struct Report + { + // distributions per volume per particle type + HistoMapMap_t histoMapMap_vol_stepSize; + HistoMapMap_t histoMapMap_vol_stepKineticEnergy; + HistoMapMap_t histoMapMap_vol_stepPseudorapidity; + HistoMapMap_t histoMapMap_vol_stepEnergyDeposit; + HistoMapMap_t histoMapMap_vol_stepEnergyNonIonDeposit; + HistoMapMap_t histoMapMap_vol_stepSecondaryKinetic; + + // distributions per material per particle type + HistoMapMap_t histoMapMap_mat_stepSize; + HistoMapMap_t histoMapMap_mat_stepKineticEnergy; + HistoMapMap_t histoMapMap_mat_stepPseudorapidity; + HistoMapMap_t histoMapMap_mat_stepEnergyDeposit; + HistoMapMap_t histoMapMap_mat_stepEnergyNonIonDeposit; + HistoMapMap_t histoMapMap_mat_stepSecondaryKinetic; + + // distributions per process per particle type + HistoMapMap_t histoMapMap_prc_stepSize; + HistoMapMap_t histoMapMap_prc_stepKineticEnergy; + HistoMapMap_t histoMapMap_prc_stepPseudorapidity; + HistoMapMap_t histoMapMap_prc_stepEnergyDeposit; + HistoMapMap_t histoMapMap_prc_stepEnergyNonIonDeposit; + HistoMapMap_t histoMapMap_prc_stepSecondaryKinetic; + + // n steps + HistoMapMap_t histoMapMap_numberOfSteps; + HistoMapMap_t histoMapMap_numberOfStepsPerInitialE; + HistoMapMap_t histoMapMap_InitialE; + HistoMapMap_t histoMapMap_stepSecondaryKinetic; + + // 2D maps + HistoMapMap_t histoMapMap2D_vol_RZ; + HistoMapMap_t histoMapMap2D_mat_RZ; + HistoMapMap_t histoMapMap2D_prc_RZ; + + // rather complicated function that merges two maps + void mergeMaps(HistoMapMap_t &selfMap, HistoMapMap_t refMap); + + // function needed by ActionToolBaseReport base class + void merge(const Report & rep); + }; + + // configurable properties + struct Config + { + bool do2DHistograms = false; + bool doGeneralHistograms = false; + }; + + /// ctor + StepHistogram(const Config&); + + const Report& getReport() const { return m_report; } + + private: + // report + Report m_report; + + /// configuration data + Config m_config; + + // initialize and fill histogram in a map + void InitializeFillHistogram2D(HistoMapMap_t &hMapMap, const char* suffix, + G4String pdgId, G4String vol, + int nbinsx, double xmin, double xmax, + int nbinsy, double ymin, double ymax, + double valuex, double valuey, double weight); + + void InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix, + G4String pdgId, G4String vol, + int nbins, double xmin, double xmax, double value, double weight); + + void InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix, + G4String pdgId, G4String vol, + int nbins, double *edges, double value, double weight); + + float m_initialKineticEnergyOfStep; + int m_trackID; + + }; // class StepHistogram + +} // namespace G4UA + +#endif // G4DEBUGGINGTOOLS_StepHistogram_H diff --git a/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogramTool.cxx b/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogramTool.cxx new file mode 100644 index 0000000000000000000000000000000000000000..33f889b5788fddec2964eda00440a07fc0944623 --- /dev/null +++ b/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogramTool.cxx @@ -0,0 +1,120 @@ +/* + Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +*/ + +#include "CxxUtils/make_unique.h" + +#include "GaudiKernel/ITHistSvc.h" + +#include "StepHistogramTool.h" + +namespace G4UA{ + + + StepHistogramTool::StepHistogramTool(const std::string& type, + const std::string& name, + const IInterface* parent): + ActionToolBaseReport<StepHistogram>(type, name, parent), + m_config(), + m_histSvc("THistSvc", name){ + + declareProperty("do2DHistograms", m_config.do2DHistograms); + declareProperty("doGeneralHistograms", m_config.doGeneralHistograms); + } + + std::unique_ptr<StepHistogram> StepHistogramTool::makeAction(){ + ATH_MSG_DEBUG("makeAction"); + auto action = CxxUtils::make_unique<StepHistogram>(m_config); + return std::move(action); + } + + StatusCode StepHistogramTool::initialize(){ + ATH_CHECK(m_histSvc.retrieve()); + return StatusCode::SUCCESS; + } + + StatusCode StepHistogramTool::queryInterface(const InterfaceID& riid, void** ppvIf){ + if(riid == IG4EventActionTool::interfaceID()) { + *ppvIf = (IG4EventActionTool*) this; + addRef(); + return StatusCode::SUCCESS; + } + if(riid == IG4SteppingActionTool::interfaceID()) { + *ppvIf = (IG4SteppingActionTool*) this; + addRef(); + return StatusCode::SUCCESS; + } + if(riid == IG4RunActionTool::interfaceID()) { + *ppvIf = (IG4RunActionTool*) this; + addRef(); + return StatusCode::SUCCESS; + } + return ActionToolBase<StepHistogram>::queryInterface(riid, ppvIf); + } + + void StepHistogramTool::BookHistograms(StepHistogram::HistoMapMap_t &hMap, const char* suffix, const char* subfolder) + { + for (auto const& x : hMap) + { + ATH_MSG_INFO("Currently in volume:\t" << x.first << " got histoMap " << x.second); + for (auto const& hm : x.second) + { + ATH_MSG_INFO("Registering histogram:\t" << hm.first); + std::ostringstream stringStream; + stringStream << "/stepHisto/" << subfolder << x.first << "/" << suffix << hm.first; + if ( m_histSvc->regHist(stringStream.str().c_str(), hm.second).isFailure() ) { + ATH_MSG_ERROR("Could not register histogram!"); + } + } + } + } + + StatusCode StepHistogramTool::finalize() { + ATH_MSG_INFO("Preparing histograms..."); + + mergeReports(); + + if (m_histSvc) { + BookHistograms(m_report.histoMapMap_vol_stepSize, "stepLength/", "volumes/"); + BookHistograms(m_report.histoMapMap_vol_stepKineticEnergy, "stepKineticEnergy/", "volumes/"); + BookHistograms(m_report.histoMapMap_vol_stepPseudorapidity, "stepPseudorapidity/", "volumes/"); + BookHistograms(m_report.histoMapMap_vol_stepEnergyDeposit, "stepEnergyDeposit/", "volumes/"); + BookHistograms(m_report.histoMapMap_vol_stepEnergyNonIonDeposit, "stepEnergyNonIonDeposit/", "volumes/"); + BookHistograms(m_report.histoMapMap_vol_stepSecondaryKinetic, "stepSecondaryKinetic/", "volumes/"); + + BookHistograms(m_report.histoMapMap_mat_stepSize, "stepLength/", "materials/"); + BookHistograms(m_report.histoMapMap_mat_stepKineticEnergy, "stepKineticEnergy/", "materials/"); + BookHistograms(m_report.histoMapMap_mat_stepPseudorapidity, "stepPseudorapidity/", "materials/"); + BookHistograms(m_report.histoMapMap_mat_stepEnergyDeposit, "stepEnergyDeposit/", "materials/"); + BookHistograms(m_report.histoMapMap_mat_stepEnergyNonIonDeposit, "stepEnergyNonIonDeposit/", "materials/"); + BookHistograms(m_report.histoMapMap_mat_stepSecondaryKinetic, "stepSecondaryKinetic/", "materials/"); + + BookHistograms(m_report.histoMapMap_prc_stepSize, "stepLength/", "processes/"); + BookHistograms(m_report.histoMapMap_prc_stepKineticEnergy, "stepKineticEnergy/", "processes/"); + BookHistograms(m_report.histoMapMap_prc_stepPseudorapidity, "stepPseudorapidity/", "processes/"); + BookHistograms(m_report.histoMapMap_prc_stepEnergyDeposit, "stepEnergyDeposit/", "processes/"); + BookHistograms(m_report.histoMapMap_prc_stepEnergyNonIonDeposit, "stepEnergyNonIonDeposit/", "processes/"); + BookHistograms(m_report.histoMapMap_prc_stepSecondaryKinetic, "stepSecondaryKinetic/", "processes/"); + + if (m_config.doGeneralHistograms) { + BookHistograms(m_report.histoMapMap_numberOfSteps, "numberOfSteps/", "nSteps/"); + BookHistograms(m_report.histoMapMap_numberOfStepsPerInitialE, "numberOfStepsPerInitialE/", "nSteps/"); + BookHistograms(m_report.histoMapMap_InitialE, "InitialE/", "nSteps/"); + BookHistograms(m_report.histoMapMap_stepSecondaryKinetic, "stepKineticEnergy/", "nSteps/"); + } + + if (m_config.do2DHistograms) { + BookHistograms(m_report.histoMapMap2D_vol_RZ, "2DMaps/", "volumes/"); + BookHistograms(m_report.histoMapMap2D_mat_RZ, "2DMaps/", "materials/"); + BookHistograms(m_report.histoMapMap2D_prc_RZ, "2DMaps/", "processes/"); + } + } + else { + ATH_MSG_WARNING("HistSvc not initialized..."); + } + + return StatusCode::SUCCESS; + } + + +} // namespace G4UA diff --git a/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogramTool.h b/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogramTool.h new file mode 100644 index 0000000000000000000000000000000000000000..e17d60b06f8271838325bb81b3a6d4c504fcc191 --- /dev/null +++ b/Simulation/G4Utilities/G4DebuggingTools/src/StepHistogramTool.h @@ -0,0 +1,67 @@ +/* + Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef G4DEBUGGINGTOOLS_G4UA__STEPHISTOGRAMTOOL_H +#define G4DEBUGGINGTOOLS_G4UA__STEPHISTOGRAMTOOL_H +#include "G4AtlasInterfaces/IG4EventActionTool.h" +#include "G4AtlasInterfaces/IG4SteppingActionTool.h" +#include "G4AtlasInterfaces/IG4RunActionTool.h" +#include "G4AtlasTools/ActionToolBase.h" +#include "StepHistogram.h" + +class ITHistSvc; + +namespace G4UA{ + + + /// @class StepHistogramTool + /// @brief Tool which manages the StepHistogram action. + /// + /// Create the StepHistogram for each worker thread + /// + /// @author Miha Muskinja + /// + + class StepHistogramTool: public ActionToolBaseReport<StepHistogram>, + public IG4EventActionTool, + public IG4RunActionTool, + public IG4SteppingActionTool + { + + public: + /// standard tool ctor + StepHistogramTool(const std::string& type, const std::string& name,const IInterface* parent); + /// return the event action + virtual G4UserEventAction* getEventAction() override final + { return static_cast<G4UserEventAction*>( getAction() ); } + /// return the stepping action + virtual G4UserSteppingAction* getSteppingAction() override final + { return static_cast<G4UserSteppingAction*>( getAction() ); } + /// return the run action + virtual G4UserRunAction* getRunAction() override final + { return static_cast<G4UserRunAction*>( getAction() ); } + /// gaudi's interface handling + virtual StatusCode queryInterface(const InterfaceID& riid, void** ppvInterface) override; + + // initialize + virtual StatusCode initialize() override; + // finalize + virtual StatusCode finalize() override; + + protected: + /// Create action for this thread + virtual std::unique_ptr<StepHistogram> makeAction() override final; + + private: + StepHistogram::Config m_config; + + ServiceHandle<ITHistSvc> m_histSvc; + + void BookHistograms(StepHistogram::HistoMapMap_t &hMap, const char* suffix, const char* subfolder = ""); + + }; // class StepHistogramTool + + +} // namespace G4UA +#endif diff --git a/Simulation/G4Utilities/G4DebuggingTools/src/components/G4DebuggingTools_entries.cxx b/Simulation/G4Utilities/G4DebuggingTools/src/components/G4DebuggingTools_entries.cxx index 1da7b4e69a98cff5e698f206696934f4a249a9ce..6472f06cda90f247d65c50e093518efd5ac80ca6 100644 --- a/Simulation/G4Utilities/G4DebuggingTools/src/components/G4DebuggingTools_entries.cxx +++ b/Simulation/G4Utilities/G4DebuggingTools/src/components/G4DebuggingTools_entries.cxx @@ -6,6 +6,7 @@ #include "../VerboseSelectorTool.h" #include "../CheckActivationTool.h" #include "../StepNtupleTool.h" +#include "../StepHistogramTool.h" #include "../VolumeDebuggerTool.h" #include "../Geant4SetupCheckerTool.h" @@ -15,6 +16,7 @@ DECLARE_TOOL_FACTORY( G4UA::G4AtlantisDumperTool ) DECLARE_TOOL_FACTORY( G4UA::VerboseSelectorTool ) DECLARE_TOOL_FACTORY( G4UA::CheckActivationTool ) DECLARE_TOOL_FACTORY( G4UA::StepNtupleTool ) +DECLARE_TOOL_FACTORY( G4UA::StepHistogramTool ) DECLARE_TOOL_FACTORY( G4UA::VolumeDebuggerTool ) DECLARE_TOOL_FACTORY( G4UA::Geant4SetupCheckerTool ) @@ -25,6 +27,7 @@ DECLARE_FACTORY_ENTRIES( G4DebuggingTools ) { DECLARE_TOOL( G4UA::VerboseSelectorTool ) DECLARE_TOOL( G4UA::CheckActivationTool ) DECLARE_TOOL( G4UA::StepNtupleTool ) + DECLARE_TOOL( G4UA::StepHistogramTool ) DECLARE_TOOL( G4UA::VolumeDebuggerTool ) DECLARE_TOOL( G4UA::Geant4SetupCheckerTool ) }