manual sweep into master of 16464 from athenaprivate1

Simulation/G4Utilities/G4DebuggingTools/python/G4DebuggingToolsConfig.py

@@ -60,3 +60,16 @@ def getGeant4SetupCheckerTool(name="G4UA::Geant4SetupCheckerTool", **kwargs):
             kwargs.setdefault(prop,value)
     # Set up the user action
     return CfgMgr.G4UA__Geant4SetupCheckerTool(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)
+

Simulation/G4Utilities/G4DebuggingTools/python/G4DebuggingToolsConfigDb.py

@@ -10,3 +10,4 @@ addTool("G4DebuggingTools.G4DebuggingToolsConfig.getEnergyConservationTestTool",
 addTool("G4DebuggingTools.G4DebuggingToolsConfig.getHyperspaceCatcherTool", "G4UA::HyperspaceCatcherTool")
 addTool("G4DebuggingTools.G4DebuggingToolsConf.G4UA__CheckActivationTool", "G4UA::CheckActivationTool")
 addTool("G4DebuggingTools.G4DebuggingToolsConfig.getGeant4SetupCheckerTool","G4UA::Geant4SetupCheckerTool")
+addTool("G4DebuggingTools.G4DebuggingToolsConfig.getStepHistogramTool", "G4UA::StepHistogramTool")

Simulation/G4Utilities/G4DebuggingTools/share/preInclude.StepHistogram.py

+## 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')
+simFlags.UserActionConfig.addConfig('G4UA::StepHistogramTool', 'doGeneralHistograms', True)
+simFlags.UserActionConfig.addConfig('G4UA::StepHistogramTool', 'do2DHistograms', False)

Simulation/G4Utilities/G4DebuggingTools/src/G4DebuggingHelper.cxx

+#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

Simulation/G4Utilities/G4DebuggingTools/src/G4DebuggingHelper.h

+#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

Simulation/G4Utilities/G4DebuggingTools/src/StepHistogram.cxx

+/*
+  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 

Simulation/G4Utilities/G4DebuggingTools/src/StepHistogram.h

+/*
+  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