diff --git a/include/DDML/EndcapTriggerTwoAngleBIBAE.h b/include/DDML/EndcapTriggerTwoAngleBIBAE.h index d38bd14..8dfc3ec 100644 --- a/include/DDML/EndcapTriggerTwoAngleBIBAE.h +++ b/include/DDML/EndcapTriggerTwoAngleBIBAE.h @@ -28,4 +28,4 @@ class EndcapTriggerTwoAngleBIBAE : public TriggerInterface { }; } // namespace ddml -#endif \ No newline at end of file +#endif diff --git a/include/DDML/LoadHdf5.h b/include/DDML/LoadHdf5.h index 74bdf02..e104b7a 100644 --- a/include/DDML/LoadHdf5.h +++ b/include/DDML/LoadHdf5.h @@ -51,6 +51,9 @@ class LoadHdf5 : public InferenceInterface { // shower library dimensions std::vector m_dimsOut{}; + // dimensions of shower library + int m_rank; + // properties for plugin std::string m_filePath = {}; diff --git a/include/DDML/OctogonalBarrelTrigger.h b/include/DDML/OctogonalBarrelTrigger.h index 4851d0f..5e66849 100644 --- a/include/DDML/OctogonalBarrelTrigger.h +++ b/include/DDML/OctogonalBarrelTrigger.h @@ -28,4 +28,4 @@ class OctogonalBarrelTrigger : public TriggerInterface { }; } // namespace ddml -#endif \ No newline at end of file +#endif diff --git a/include/DDML/PionCloudsModel.h b/include/DDML/PionCloudsModel.h new file mode 100644 index 0000000..3830e59 --- /dev/null +++ b/include/DDML/PionCloudsModel.h @@ -0,0 +1,57 @@ +#ifndef PionClouds_H +#define PionClouds_H + +#include "DDML/FastMLShower.h" +#include "DDML/ModelInterface.h" + +namespace ddml { + +/** Class for running a point cloud based ML model for fast shower simulation. + * Assumes a cartesian (x,y) coordinates defining the calorimeter planes (layers) and z the depth + * of the calorimeter. + * + * Based on BiBAETwoAngleModel. + * + * Implemented here for the PionClouds model intended for hadron shower simulation (ECAL+HCAL). + * + * @author A.Korol, DESY + * @author P. McKeown, CERN + * @date Feb. 2025 + */ + +class PionCloudsModel : public ModelInterface { +public: + PionCloudsModel() = default; + + virtual ~PionCloudsModel() = default; + + /// declare the proerties needed for the plugin + void declareProperties(dd4hep::sim::Geant4Action* plugin) { + plugin->declareProperty("LatentVectorSize", this->m_latentSize); + } + + /** prepare the input vector and resize the output vector for this model + * based on the current FastTrack (e.g. extract kinetic energy and incident + * angles.) + */ + virtual void prepareInput(G4FastTrack const& aFastTrack, G4ThreeVector const& localDir, InputVecs& inputs, + TensorDimVecs& tensDims, std::vector& output); + + /** create a vector of spacepoints per layer interpreting the model output + */ + virtual void convertOutput(G4FastTrack const& aFastTrack, G4ThreeVector const& localDir, + const std::vector& output, std::vector& spacepoints); + +private: + /// model properties for plugin + int m_numPoints = 2600; // 4148; //2600; //number of points in the shower + size_t m_nDims = 4; // Number of dimensions + int m_latentSize = 3; // number of input features (energy, theta, phi) + int m_maxNumElements = m_numPoints * 4; // number of space points in the output multiplied by 4 (x,y,z,energy) + // number of layers for ILD is 78: int m_nLayer = 78; + + TensorDimVecs m_tensDims = {{1, 1}, {1, 1}, {1, 1}, {1, 3}}; +}; + +} // namespace ddml +#endif diff --git a/include/DDML/PolyhedraBarrelGeometry.h b/include/DDML/PolyhedraBarrelGeometry.h index 5879a96..520d2ec 100644 --- a/include/DDML/PolyhedraBarrelGeometry.h +++ b/include/DDML/PolyhedraBarrelGeometry.h @@ -12,6 +12,11 @@ namespace ddml { * * @author F.Gaede, DESY * @date Mar 2023 + * + * Addiional option included to support Hadronic shower simulation in ECAL + HCAL + * @author P. McKeown, CERN + * @date Feb 2025 + * */ class PolyhedraBarrelGeometry : public GeometryInterface { @@ -26,7 +31,10 @@ class PolyhedraBarrelGeometry : public GeometryInterface { /// declare the proerties needed for the plugin void declareProperties(dd4hep::sim::Geant4Action* plugin) { plugin->declareProperty("Detector", this->m_detector); + plugin->declareProperty("isHadShower", this->m_isHadShower); + plugin->declareProperty("HadDetector", this->m_hadDetector); plugin->declareProperty("Symmetry", this->m_nSymmetry); + plugin->declareProperty("HadSymmetry", this->m_nHadSymmetry); plugin->declareProperty("CorrectForAngles", this->m_correctForAngles); } @@ -50,6 +58,9 @@ class PolyhedraBarrelGeometry : public GeometryInterface { std::string m_detector = {"EcalBarrel"}; int m_nSymmetry = 8; bool m_correctForAngles = false; + bool m_isHadShower; //= true; + std::string m_hadDetector = {"HcalBarrel"}; + int m_nHadSymmetry = m_nSymmetry; }; } // namespace ddml diff --git a/include/DDML/TriggerInterface.h b/include/DDML/TriggerInterface.h index 1141165..801e780 100644 --- a/include/DDML/TriggerInterface.h +++ b/include/DDML/TriggerInterface.h @@ -22,4 +22,4 @@ class TriggerInterface { } // namespace ddml -#endif \ No newline at end of file +#endif diff --git a/scripts/ddsim_steer.py b/scripts/ddsim_steer.py index 94dec80..785d160 100644 --- a/scripts/ddsim_steer.py +++ b/scripts/ddsim_steer.py @@ -394,16 +394,27 @@ def aiDanceTorch(kernel): CaloClouds = True L2LFlows = False old_DD4hep = False ## use for DD4hep versions/commits before ~ Apr 21st 2023 + hadrons = False if ild == True: ml_barrel_name = "EcalBarrel" ml_barrel_symmetry = 8 ml_endcap_name = "EcalEndcap" + + ## For hadron shower fast simulation + ml_had_barrel_name = "HcalBarrel" + ml_had_barrel_symmetry = 8 + ml_had_endcap_name = "HcalEndcap" else: ml_barrel_name = "ECalBarrel" ml_barrel_symmetry = 12 ml_endcap_name = "ECalEndcap" + ## For hadron shower fast simulation is needed + ml_had_barrel_name = "HCalBarrel" + ml_had_barrel_symmetry = 12 + ml_had_endcap_name = "HCalEndcap" + if BIBAE == True and Two_Angle == False: ml_file = "../models/BIBAE_Full_PP_cut.pt" ml_model = "RegularGridBIBAEPolyhedraBarrelTorchModel/BarrelModelTorch" @@ -459,6 +470,7 @@ def aiDanceTorch(kernel): model = DetectorConstruction(kernel, str(ml_model)) ## # Mandatory model parameters + model.isHadShower = False model.RegionName = "EcalBarrelRegion" model.Detector = ml_barrel_name model.Symmetry = ml_barrel_symmetry @@ -515,16 +527,28 @@ def LoadHdf5(kernel): BIBAE = True Two_Angle = True old_DD4hep = False ## use for DD4hep versions/commits before ~ Apr 21st 2023 + hadrons = True if ild == True: ml_barrel_name = "EcalBarrel" ml_barrel_symmetry = 8 ml_endcap_name = "EcalEndcap" + + ## For hadron shower fast simulation + ml_had_barrel_name = "HcalBarrel" + ml_had_barrel_symmetry = 8 + ml_had_endcap_name = "HcalEndcap" + else: ml_barrel_name = "ECalBarrel" ml_barrel_symmetry = 12 ml_endcap_name = "ECalEndcap" + ## For hadron shower fast simulation is needed + ml_had_barrel_name = "HCalBarrel" + ml_had_barrel_symmetry = 12 + ml_had_endcap_name = "HCalEndcap" + if BIBAE == True and Two_Angle == True: ml_file = "../models/photons-E5050A-theta9090A-phi9090-p1.hdf5" ml_model = ( @@ -533,6 +557,11 @@ def LoadHdf5(kernel): ml_model_1 = "LoadHDF5RegularGridTwoAngleBIBAEModelEndcap/EndcapModelTorch" ml_correct_angles = False + if hadrons == True: + ml_model_had = "LoadHDF5PionCloudsPCHadronModelPolyhedraBarrel/BarrelModelTorch" + ml_had_file = "../models/PionClouds_50GeV_sp_scaled.h5" + ml_correct_angles = False + from g4units import GeV, MeV # DO NOT REMOVE OR MOVE!!!!! (EXCLAMATION MARK) from DDG4 import DetectorConstruction, Geant4, PhysicsList @@ -557,6 +586,8 @@ def LoadHdf5(kernel): seq.adopt(sensitives) # ----------------- + """ + ## EM in Barrel model = DetectorConstruction(kernel, str(ml_model)) ## # Mandatory model parameters @@ -578,7 +609,9 @@ def LoadHdf5(kernel): model.enableUI() seq.adopt(model) + """ # ------------------- + ## EM in Endcap model1 = DetectorConstruction(kernel, str(ml_model_1)) ## # Mandatory model parameters @@ -599,12 +632,38 @@ def LoadHdf5(kernel): model1.enableUI() seq.adopt(model1) + + # ------------------- + ## Hadrons in Barrel + modelHad1 = DetectorConstruction(kernel, str(ml_model_had)) + + ## # Mandatory model parameters + modelHad1.isHadShower = True + modelHad1.RegionName = ( + "EcalBarrelRegion" # or "HcalBarrelRegion" ## hadron model triggers in ecal + ) + modelHad1.Detector = ml_barrel_name + modelHad1.HadDetector = ml_had_barrel_name + modelHad1.Symmetry = ml_barrel_symmetry + modelHad1.HadSymmetry = ml_had_barrel_symmetry + modelHad1.Enable = True + modelHad1.CorrectForAngles = ml_correct_angles + # Energy boundaries are optional: Units are GeV + modelHad1.ApplicableParticles = {"pi+"} + modelHad1.Etrigger = {"pi+": 10.0 * GeV} # trigger on lower training threshold + modelHad1.FilePath = ml_had_file + # model.OptimizeFlag = 1 + # model.IntraOpNumThreads = 1 + + modelHad1.enableUI() + seq.adopt(modelHad1) + # ------------------- # Now build the physics list: phys = kernel.physicsList() ph = PhysicsList(kernel, str("Geant4FastPhysics/FastPhysicsList")) - ph.EnabledParticles = ["e+", "e-", "gamma"] + ph.EnabledParticles = ["e+", "e-", "gamma", "pi+"] ph.BeVerbose = True ph.enableUI() phys.adopt(ph) diff --git a/scripts/test_onnx.mac b/scripts/test_onnx.mac index b7a22b4..c1f4695 100644 --- a/scripts/test_onnx.mac +++ b/scripts/test_onnx.mac @@ -3,13 +3,24 @@ # example script to run inference on a GAN w/ ONNX runtime # --- the particle gun -/gps/pos/centre 0 0 0 -/gps/particle gamma -#/gps/particle pi- -/gps/energy 50 GeV + +# At right side face (perpendicular to x=0) +#/gps/pos/centre 18.47 0 0 #0 180.47 150 #### this is ILD calo Face #120 180.48 120 #30 180.48 0 #30 0 0 +#30 0 0 + +#/gps/pos/centre 0 0 0 #0 180.47 150 #### this is ILD calo Face #120 180.48 120 #30 180.48 0 #30 0 0 +#30 0 0 +#/gps/particle gamma # for SimpleCalo -/gps/direction -0.4152 1 0 #(pi/8 + pi/4) #1 0.414213562 0 #(pi/8) #-1 2.414213562 0 #(pi/8 + 1 deg) # 0.25 0.25 0.5 #0 0.5 0.655 # ( 37 deg in theta- most you can get in endcap from IP ) +#/gps/direction 1 0 0 #0 1 0 #-0.4152 1 0 #(pi/8 + pi/4) #1 0.414213562 0 #(pi/8) #-1 2.414213562 0 #(pi/8 + 1 deg) # 0.25 0.25 0.5 #0 0.5 0.655 # ( 37 deg in theta- most you can get in endcap from IP ) + +/gps/particle pi+ +/gps/energy 50 GeV + +/gps/pos/centre -5 0 -15 +/gps/direction -0.037864591009775746 0.9992828792427412 0 # 90 deg imact for 50 GeV p+ + # phi barrel #1 0.414213562 0 #(pi/8) @@ -64,5 +75,5 @@ #/gps/direction 0.1 -0.8 .1 -/run/beamOn 10 #100 +/run/beamOn 2000 #100 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b56fb29..2882c79 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,6 +7,7 @@ set(sources Par04ExampleVAE.cc RegularGridBIBAEModel.cc RegularGridTwoAngleBIBAEModel.cc + PionCloudsModel.cc OctogonalBarrelTrigger.cc EndcapTriggerTwoAngleBIBAE.cc CaloCloudsTwoAngleModel.cc @@ -82,6 +83,7 @@ set(headers ${PROJECT_SOURCE_DIR}/include/DDML/Par04ExampleVAE.h ${PROJECT_SOURCE_DIR}/include/DDML/RegularGridBIBAEModel.h ${PROJECT_SOURCE_DIR}/include/DDML/RegularGridTwoAngleBIBAEModel.h + ${PROJECT_SOURCE_DIR}/include/DDML/PionCloudsModel.h ${PROJECT_SOURCE_DIR}/include/DDML/PolyhedraBarrelGeometry.h ${PROJECT_SOURCE_DIR}/include/DDML/ModelInterface.h ${PROJECT_SOURCE_DIR}/include/DDML/InferenceInterface.h diff --git a/src/LoadHdf5.cc b/src/LoadHdf5.cc index 7bfebe1..1e27139 100644 --- a/src/LoadHdf5.cc +++ b/src/LoadHdf5.cc @@ -19,7 +19,7 @@ void LoadHdf5::initialize() { // inputs and TensorDimVecs unused // Open dataset + dataspace - std::string datasetName = "layers"; + std::string datasetName = "spase_points"; //"events"; //"spase_points"; //"layers"; H5::DataSet dataset = m_file.openDataSet(datasetName); dd4hep::printout(dd4hep::DEBUG, "LoadHdf5::initialize", "Accessed HDF5 dataset"); H5::DataSpace dataspace = dataset.getSpace(); @@ -28,8 +28,8 @@ void LoadHdf5::initialize() { // bool H5::DataSpace::isSimple() const; // Get dimensions - int rank = dataspace.getSimpleExtentNdims(); - std::vector dims_out(rank); + m_rank = dataspace.getSimpleExtentNdims(); + std::vector dims_out(m_rank); dataspace.getSimpleExtentDims(dims_out.data(), nullptr); // Calculate the total number of elements @@ -39,13 +39,19 @@ void LoadHdf5::initialize() { m_dimsOut = dims_out; - assert(rank == 4); // assuming 4 dimensional input + assert(m_rank == 3 || m_rank == 4); // ensure either 3 dimensional PionCloud input or 4 dimensional input + // if m_rank == 4 // index 0: shower number // index 1, 2, 3: x, y, z cell number - dd4hep::printout(dd4hep::DEBUG, "LoadHdf5::initialize", "Rank %i", rank); + // else if m_rank == 3 + // index 0: shower number + // index 1: number of points + // index 2: 4 point dimensions (x, y, z, E) - currently in ILD coordinates + + dd4hep::printout(dd4hep::DEBUG, "LoadHdf5::initialize", "Rank %i", m_rank); - for (int i = 0, N = rank; i < N; ++i) { + for (int i = 0, N = m_rank; i < N; ++i) { dd4hep::printout(dd4hep::DEBUG, "LoadHdf5::initialize", "dimension %i: %lu", i, (unsigned long)(m_dimsOut[i])); } @@ -71,23 +77,29 @@ void LoadHdf5::runInference(const InputVecs&, const TensorDimVecs&, std::vector< m_isInitialized = true; } - // keep track of shower index in file - // LoadHdf5::incrementCounter(); - // If counter exceeds number of showers in file, reset if (m_count > m_totalSize) { m_count = 0; } // select shower from library - std::vector shower(m_library.begin() + m_count * m_dimsOut[1] * m_dimsOut[2] * m_dimsOut[3], - m_library.begin() + (m_count + 1) * m_dimsOut[1] * m_dimsOut[2] * m_dimsOut[3]); - - // enforce length of shower - assert(shower.size() == m_dimsOut[1] * m_dimsOut[2] * m_dimsOut[3]); - - // write output - output = std::move(shower); + if (m_rank == 4) { + std::vector shower(m_library.begin() + m_count * m_dimsOut[1] * m_dimsOut[2] * m_dimsOut[3], + m_library.begin() + (m_count + 1) * m_dimsOut[1] * m_dimsOut[2] * m_dimsOut[3]); + assert(shower.size() == m_dimsOut[1] * m_dimsOut[2] * m_dimsOut[3]); + + // write output + output = std::move(shower); + } else if (m_rank == 3) { + std::vector shower(m_library.begin() + m_count * m_dimsOut[1] * m_dimsOut[2], + m_library.begin() + (m_count + 1) * m_dimsOut[1] * m_dimsOut[2]); + assert(shower.size() == m_dimsOut[1] * m_dimsOut[2]); + + // write output + output = std::move(shower); + } else { + throw std::runtime_error("Shower library does not have the correct dimensions!"); + } ++m_count; } diff --git a/src/MLModelActions.cc b/src/MLModelActions.cc index 9ce1aae..64d601d 100644 --- a/src/MLModelActions.cc +++ b/src/MLModelActions.cc @@ -21,6 +21,7 @@ #include "DDML/L2LFlowsModel.h" #include "DDML/OctogonalBarrelTrigger.h" #include "DDML/Par04ExampleVAE.h" +#include "DDML/PionCloudsModel.h" #include "DDML/PolyhedraBarrelGeometry.h" #include "DDML/RegularGridBIBAEModel.h" #include "DDML/RegularGridGANModel.h" @@ -127,6 +128,17 @@ typedef FastMLShower< FastMLModel> // add ML trigger LoadHDF5RegularGridTwoAngleBIBAEModelEndcap; + +/// Load from HDF5 file- as an example for Hadron showers from PionClouds +// Barrel +typedef FastMLShower< + FastMLModel> // add ML trigger + LoadHDF5PionCloudsPCHadronModelPolyhedraBarrel; +// Endcap // ENDCAP IS CURRENTLY NOT IMPLEMENTED!!!! +typedef FastMLShower> // add ML trigger + LoadHDF5PionCloudsPCHadronModelEndcap; #endif } // namespace ddml @@ -156,4 +168,6 @@ DECLARE_GEANT4ACTION_NS(ddml, L2LFlowsModelEndcapTorchModel) #ifdef DDML_USE_LOAD_HDF5 DECLARE_GEANT4ACTION_NS(ddml, LoadHDF5RegularGridTwoAngleBIBAEModelPolyhedraBarrel) DECLARE_GEANT4ACTION_NS(ddml, LoadHDF5RegularGridTwoAngleBIBAEModelEndcap) +DECLARE_GEANT4ACTION_NS(ddml, LoadHDF5PionCloudsPCHadronModelPolyhedraBarrel) +DECLARE_GEANT4ACTION_NS(ddml, LoadHDF5PionCloudsPCHadronModelEndcap) #endif diff --git a/src/ONNXInference.cc b/src/ONNXInference.cc index d489bc7..9de0469 100644 --- a/src/ONNXInference.cc +++ b/src/ONNXInference.cc @@ -55,7 +55,7 @@ void ONNXInference::initialize() { std::vector input_node_names(num_input_nodes); for (std::size_t i = 0; i < num_input_nodes; i++) { #if ORT_API_VERSION < 13 - const auto input_name = AllocatedStringPtr(fSession->GetInputName(i, allocator), allocDeleter).release(); + const auto input_name = AllocatedStringPtr(m_session->GetInputName(i, allocator), allocDeleter).release(); #else const auto input_name = m_session->GetInputNameAllocated(i, allocator).release(); #endif @@ -81,7 +81,7 @@ void ONNXInference::initialize() { std::vector output_node_names(num_output_nodes); for (std::size_t i = 0; i < num_output_nodes; i++) { #if ORT_API_VERSION < 12 - const auto output_name = AllocatedStringPtr(fSession->GetOutputName(i, allocator), allocDeleter).release(); + const auto output_name = AllocatedStringPtr(m_session->GetOutputName(i, allocator), allocDeleter).release(); #else const auto output_name = m_session->GetOutputNameAllocated(i, allocator).release(); #endif diff --git a/src/PionCloudsModel.cc b/src/PionCloudsModel.cc new file mode 100644 index 0000000..ba005c0 --- /dev/null +++ b/src/PionCloudsModel.cc @@ -0,0 +1,80 @@ +#include "DDML/PionCloudsModel.h" + +#include // for G4FastTrack +#include + +#define DEBUGPRINT 0 + +namespace ddml { + +void PionCloudsModel::prepareInput(G4FastTrack const& aFastTrack, G4ThreeVector const& localDir, InputVecs& inputs, + TensorDimVecs& tensDims, std::vector& output) { + tensDims = m_tensDims; + + G4double energy = aFastTrack.GetPrimaryTrack()->GetKineticEnergy(); + + G4ThreeVector direction = aFastTrack.GetPrimaryTrack()->GetMomentumDirection(); + G4RotationMatrix rotZ; + rotZ.rotateZ(M_PI / 2.); + G4RotationMatrix rotX; + rotX.rotateX(M_PI / 2.); + // this convention is used for the local coordinates in the dataset (model was trained in this convention) + + G4ThreeVector localDir_ = localDir; + localDir_.setX(-1. * localDir_.x()); // *(-1) to align local to global convention in ddml + localDir_.setY(-1. * localDir_.y()); // *(-1) to align local to global convention in ddml + + dd4hep::printout(dd4hep::DEBUG, "PionCloudsModel::prepareInput", "DDML::localDir: (%f, %f, %f)", localDir_.x(), + localDir_.y(), localDir_.z()); + + // compute local incident angles + double r = sqrt(localDir_.x() * localDir_.x() + localDir_.y() * localDir_.y() + localDir_.z() * localDir_.z()); + double theta = acos(localDir_.z() / r) / M_PI * 180.; + double phi = atan2(localDir_.y(), localDir_.x()) / M_PI * 180.; + + dd4hep::printout(dd4hep::DEBUG, "PionCloudsModel::prepareInput", "DDML::localDir: (%f, %f)", theta, phi); + + // the input for the PionClouds is one energy and two angles (local Theta and Phi) + inputs.resize(m_latentSize); + + inputs[0].resize(1); // Energy + inputs[1].resize(1); // Theta + inputs[2].resize(1); // Phi + + // For now, assume batch size one, and just assign values + inputs[0][0] = energy / CLHEP::GeV; // E_vec[0]/100.; + inputs[1][0] = theta; // 89.*(M_PI/180.) ; //Theta_vec[0]/(90.*(M_PI/180.)); + inputs[2][0] = phi; + + dd4hep::printout(dd4hep::DEBUG, "PionCloudsModel::prepareInput", "Input_energy_tensor : %f", inputs[0][0]); + dd4hep::printout(dd4hep::DEBUG, "PionCloudsModel::prepareInput", "Input_theta_tensor : %f", inputs[1][0]); + dd4hep::printout(dd4hep::DEBUG, "PionCloudsModel::prepareInput", "Input_phi_tensor : %f", inputs[2][0]); + + // ---- resize output vector + + output.assign(m_maxNumElements, 0); +} + +// For array structure: (No. showers, No. points, dimensions(4)) +void PionCloudsModel::convertOutput(G4FastTrack const&, G4ThreeVector const&, const std::vector& output, + std::vector& spacepoints) { + const int nPoints = output.size() / m_nDims; + int layerNum = 0; + // Reshape into intermediate representation + auto reshaped = std::views::iota(0, nPoints) | std::views::transform([&output](int i) { + return std::span{output.data() + i * 4, 4}; + }); + + spacepoints.resize(nPoints); + for (const auto& values : reshaped) { + ddml::SpacePoint sp(values[0], // x // *(-1) to align local to global convention in ddml + values[2], // y // *(-1) to align local to global convention in ddml + 0., // z + values[3], // energy + 0. // time + ); + layerNum = int(values[1]); + spacepoints[layerNum].emplace_back(sp); + } +} +} // namespace ddml diff --git a/src/PolyhedraBarrelGeometry.cc b/src/PolyhedraBarrelGeometry.cc index 96e064d..fbb811d 100644 --- a/src/PolyhedraBarrelGeometry.cc +++ b/src/PolyhedraBarrelGeometry.cc @@ -26,6 +26,22 @@ void PolyhedraBarrelGeometry::initialize() { dd4hep::printout(dd4hep::ERROR, "PolyhedraBarrelGeometry::initialize", "Detector %s not found!", m_detector.c_str()); } + + // For hadronic shower simulation + if (m_isHadShower == true) { + auto det_had = theDetector.detector(m_hadDetector); + auto* cal_had = det_had.extension(); + if (cal_had) { + for (auto l_had : cal_had->layers) { + m_caloLayerDistances.push_back((l_had.distance + l_had.inner_thickness) / dd4hep::mm); + dd4hep::printout(dd4hep::INFO, "PolyhedraBarrelGeometry::initialize", "HCAL Layer distances %f", + l_had.distance + l_had.inner_thickness); + } + } else { + dd4hep::printout(dd4hep::ERROR, "PolyhedraBarrelGeometry::initialize", "Detector %s not found!", + m_hadDetector.c_str()); + } + } } int PolyhedraBarrelGeometry::phiSector(G4ThreeVector const& position) const {