From b66a3004b8f9462aeb0263dae52ee1e5f04b08f6 Mon Sep 17 00:00:00 2001 From: Peter McKeown Date: Fri, 8 Nov 2024 17:06:07 +0100 Subject: [PATCH 01/12] initial progress --- include/DDML/PolyhedraBarrelGeometry.h | 4 ++- scripts/ddsim_steer.py | 41 +++++++++++++++++++++++++- src/PolyhedraBarrelGeometry.cc | 14 +++++++++ 3 files changed, 57 insertions(+), 2 deletions(-) diff --git a/include/DDML/PolyhedraBarrelGeometry.h b/include/DDML/PolyhedraBarrelGeometry.h index ee08de2..8dbce2b 100644 --- a/include/DDML/PolyhedraBarrelGeometry.h +++ b/include/DDML/PolyhedraBarrelGeometry.h @@ -27,8 +27,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", m_isHadShower) plugin->declareProperty("Detector", this->m_detector); + plugin->declareProperty("HadDetector", this->m_HadDetector); plugin->declareProperty("Symmetry", this->m_nSymmetry); + plugin->declareProperty("HadSymmetry", this->m_nHadSymmetry); plugin->declareProperty("CorrectForAngles", this->m_correctForAngles); } diff --git a/scripts/ddsim_steer.py b/scripts/ddsim_steer.py index 94dec80..5a7ec2b 100644 --- a/scripts/ddsim_steer.py +++ b/scripts/ddsim_steer.py @@ -515,16 +515,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 +545,10 @@ def LoadHdf5(kernel): ml_model_1 = "LoadHDF5RegularGridTwoAngleBIBAEModelEndcap/EndcapModelTorch" ml_correct_angles = False + if hadrons == True: + ml_model_had = "LoadHDF5HadronPointCloudModelPolyhedraBarrel/BarrelModelTorch" + ml_correct_angles = False + from g4units import GeV, MeV # DO NOT REMOVE OR MOVE!!!!! (EXCLAMATION MARK) from DDG4 import DetectorConstruction, Geant4, PhysicsList @@ -557,6 +573,7 @@ def LoadHdf5(kernel): seq.adopt(sensitives) # ----------------- + ## EM in Barrel model = DetectorConstruction(kernel, str(ml_model)) ## # Mandatory model parameters @@ -579,6 +596,7 @@ def LoadHdf5(kernel): model.enableUI() seq.adopt(model) # ------------------- + ## EM in Endcap model1 = DetectorConstruction(kernel, str(ml_model_1)) ## # Mandatory model parameters @@ -599,12 +617,33 @@ def LoadHdf5(kernel): model1.enableUI() seq.adopt(model1) + + # ------------------- + ## Hadrons in Barrel + modelHad1 = DetectorConstruction(kernel, str()) + + ## # Mandatory model parameters + model.isHadShower = True + model.RegionName = "EcalBarrelRegion" ## hadron model triggers in ecal + model.Detector = ml_barrel_name + model.HadDetector = ml_had_barrel_name + model.Symmetry = ml_barrel_symmetry + model.HadSymmetry = ml_had_barrel_symmetry + model.Enable = True + model.CorrectForAngles = ml_correct_angles + # Energy boundaries are optional: Units are GeV + model.ApplicableParticles = {"pi+"} + model.Etrigger = {"pi+": 10.0 * GeV} # trigger on lower training threshold + model.FilePath = ml_file + # model.OptimizeFlag = 1 + # model.IntraOpNumThreads = 1 + # ------------------- # 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/src/PolyhedraBarrelGeometry.cc b/src/PolyhedraBarrelGeometry.cc index 9b2ef65..f7f4f6a 100644 --- a/src/PolyhedraBarrelGeometry.cc +++ b/src/PolyhedraBarrelGeometry.cc @@ -16,6 +16,12 @@ void PolyhedraBarrelGeometry::initialize() { auto det = theDetector.detector(m_detector); auto* cal = det.extension(); + // For hadronic shower simulation + if (m_isHadShower == true) { + auto det_had = theDetector.detector(m_HadDetector); + auto* cal_had = det_had.extension(); + } + if (cal) { for (auto l : cal->layers) { m_caloLayerDistances.push_back((l.distance + l.inner_thickness) / dd4hep::mm); @@ -23,6 +29,14 @@ void PolyhedraBarrelGeometry::initialize() { } else { std::cout << " ###### error: detector " << m_detector << " not found !" << std::endl; } + + if (m_isHadShower == true) { + if (cal_had) { + for (auto l_had : cal_had->layers) { + m_caloLayerDistances.push_back((l_had.distance + l_had.inner_thickness) / dd4hep::mm); + } + } + } } int PolyhedraBarrelGeometry::phiSector(G4ThreeVector const& position) const { From a58fd8adeae4032d2faffb1ab2ef4c8bed2c1bfe Mon Sep 17 00:00:00 2001 From: Peter McKeown Date: Tue, 19 Nov 2024 22:59:33 +0100 Subject: [PATCH 02/12] Loading and placement in ECAL successful, pending correct placement in HCAL --- include/DDML/PolyhedraBarrelGeometry.h | 6 +++++- scripts/ddsim_steer.py | 30 +++++++++++++++----------- scripts/test_onnx.mac | 6 +++--- src/CMakeLists.txt | 2 ++ src/LoadHdf5.cc | 21 +++++++++++++----- src/MLModelActions.cc | 16 ++++++++++++++ src/PolyhedraBarrelGeometry.cc | 12 +++++------ 7 files changed, 65 insertions(+), 28 deletions(-) diff --git a/include/DDML/PolyhedraBarrelGeometry.h b/include/DDML/PolyhedraBarrelGeometry.h index 8dbce2b..6a23c1b 100644 --- a/include/DDML/PolyhedraBarrelGeometry.h +++ b/include/DDML/PolyhedraBarrelGeometry.h @@ -27,7 +27,8 @@ class PolyhedraBarrelGeometry : public GeometryInterface { /// declare the proerties needed for the plugin void declareProperties(dd4hep::sim::Geant4Action* plugin) { - plugin->declareProperty("isHadShower", m_isHadShower) plugin->declareProperty("Detector", this->m_detector); + 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); @@ -54,6 +55,9 @@ class PolyhedraBarrelGeometry : public GeometryInterface { std::string m_detector = {"EcalBarrel"}; int m_nSymmetry = 8; bool m_correctForAngles = false; + bool m_isHadShower = false; + std::string m_HadDetector = {"HcalBarrel"}; + int m_nHadSymmetry = m_nSymmetry; }; } // namespace ddml diff --git a/scripts/ddsim_steer.py b/scripts/ddsim_steer.py index 5a7ec2b..71b192e 100644 --- a/scripts/ddsim_steer.py +++ b/scripts/ddsim_steer.py @@ -546,7 +546,8 @@ def LoadHdf5(kernel): ml_correct_angles = False if hadrons == True: - ml_model_had = "LoadHDF5HadronPointCloudModelPolyhedraBarrel/BarrelModelTorch" + ml_model_had = "LoadHDF5PionCloudsPCHadronModelPolyhedraBarrel/BarrelModelTorch" + ml_had_file = "../models/PionClouds_50GeV_sp.h5" ml_correct_angles = False from g4units import GeV, MeV # DO NOT REMOVE OR MOVE!!!!! (EXCLAMATION MARK) @@ -620,24 +621,27 @@ def LoadHdf5(kernel): # ------------------- ## Hadrons in Barrel - modelHad1 = DetectorConstruction(kernel, str()) + modelHad1 = DetectorConstruction(kernel, str(ml_model_had)) ## # Mandatory model parameters - model.isHadShower = True - model.RegionName = "EcalBarrelRegion" ## hadron model triggers in ecal - model.Detector = ml_barrel_name - model.HadDetector = ml_had_barrel_name - model.Symmetry = ml_barrel_symmetry - model.HadSymmetry = ml_had_barrel_symmetry - model.Enable = True - model.CorrectForAngles = ml_correct_angles + modelHad1.isHadShower = True + modelHad1.RegionName = "EcalBarrelRegion" ## 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 - model.ApplicableParticles = {"pi+"} - model.Etrigger = {"pi+": 10.0 * GeV} # trigger on lower training threshold - model.FilePath = ml_file + 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: diff --git a/scripts/test_onnx.mac b/scripts/test_onnx.mac index b7a22b4..10f6964 100644 --- a/scripts/test_onnx.mac +++ b/scripts/test_onnx.mac @@ -4,12 +4,12 @@ # --- the particle gun /gps/pos/centre 0 0 0 -/gps/particle gamma -#/gps/particle pi- +#/gps/particle gamma +/gps/particle pi+ /gps/energy 50 GeV # 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 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 ) # phi barrel #1 0.414213562 0 #(pi/8) 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..2cb8360 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"; //"layers"; H5::DataSet dataset = m_file.openDataSet(datasetName); dd4hep::printout(dd4hep::DEBUG, "LoadHdf5::initialize", "Accessed HDF5 dataset"); H5::DataSpace dataspace = dataset.getSpace(); @@ -39,10 +39,15 @@ void LoadHdf5::initialize() { m_dimsOut = dims_out; - assert(rank == 4); // assuming 4 dimensional input +// assert(rank == 4); // assuming 4 dimensional input // index 0: shower number // index 1, 2, 3: x, y, z cell number + assert(rank == 3); // assuming 3 dimensional PionCloud input + // 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", rank); for (int i = 0, N = rank; i < N; ++i) { @@ -80,11 +85,17 @@ void LoadHdf5::runInference(const InputVecs&, const TensorDimVecs&, std::vector< } // 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]); + //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]); + + // select shower from library + 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]); // enforce length of shower - assert(shower.size() == m_dimsOut[1] * m_dimsOut[2] * m_dimsOut[3]); + //assert(shower.size() == m_dimsOut[1] * m_dimsOut[2] * m_dimsOut[3]); + + assert(shower.size() == m_dimsOut[1] * m_dimsOut[2]); // write output output = std::move(shower); diff --git a/src/MLModelActions.cc b/src/MLModelActions.cc index ab56dd1..5e88b6b 100644 --- a/src/MLModelActions.cc +++ b/src/MLModelActions.cc @@ -25,6 +25,7 @@ #include "DDML/RegularGridBIBAEModel.h" #include "DDML/RegularGridGANModel.h" #include "DDML/RegularGridTwoAngleBIBAEModel.h" +#include "DDML/PionCloudsModel.h" #include "DDML/TriggerInterface.h" namespace ddml { @@ -127,8 +128,21 @@ typedef FastMLShower< FastMLModel> // add ML trigger LoadHDF5RegularGridTwoAngleBIBAEModelEndcap; + +/// Load from HDF5 file- as an example for Hadron showers from PionClouds +// Barrel +typedef FastMLShower> // add ML trigger + LoadHDF5PionCloudsPCHadronModelPolyhedraBarrel; +// Endcap // ENDCAP IS CURRENTLY NOT IMPLEMENTED!!!! +typedef FastMLShower< + FastMLModel> // add ML trigger + LoadHDF5PionCloudsPCHadronModelEndcap; #endif + } // namespace ddml #include @@ -156,4 +170,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/PolyhedraBarrelGeometry.cc b/src/PolyhedraBarrelGeometry.cc index f7f4f6a..829ddb0 100644 --- a/src/PolyhedraBarrelGeometry.cc +++ b/src/PolyhedraBarrelGeometry.cc @@ -16,12 +16,6 @@ void PolyhedraBarrelGeometry::initialize() { auto det = theDetector.detector(m_detector); auto* cal = det.extension(); - // For hadronic shower simulation - if (m_isHadShower == true) { - auto det_had = theDetector.detector(m_HadDetector); - auto* cal_had = det_had.extension(); - } - if (cal) { for (auto l : cal->layers) { m_caloLayerDistances.push_back((l.distance + l.inner_thickness) / dd4hep::mm); @@ -30,12 +24,18 @@ void PolyhedraBarrelGeometry::initialize() { std::cout << " ###### error: detector " << m_detector << " not found !" << std::endl; } + // 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); } } + else { + std::cout << " ###### error: detector " << m_HadDetector << " not found !" << std::endl; + } } } From 424f0a4e75cc02df7831b8247075e568db6cf38a Mon Sep 17 00:00:00 2001 From: Peter McKeown Date: Mon, 17 Feb 2025 17:25:38 +0100 Subject: [PATCH 03/12] Working implmentation of hadron shower loading --- include/DDML/PionCloudsModel.h | 74 +++++++++++++ include/DDML/PolyhedraBarrelGeometry.h | 7 +- scripts/ddsim_steer.py | 4 +- scripts/test_onnx.mac | 4 +- src/Geant4FastHitMakerGlobal.cc | 7 ++ src/PionCloudsModel.cc | 142 +++++++++++++++++++++++++ src/PolyhedraBarrelGeometry.cc | 8 ++ 7 files changed, 242 insertions(+), 4 deletions(-) create mode 100644 include/DDML/PionCloudsModel.h create mode 100644 src/PionCloudsModel.cc diff --git a/include/DDML/PionCloudsModel.h b/include/DDML/PionCloudsModel.h new file mode 100644 index 0000000..0bb9bfc --- /dev/null +++ b/include/DDML/PionCloudsModel.h @@ -0,0 +1,74 @@ +#ifndef PionClouds_H +#define PionClouds_H + +#include "DDML/ModelInterface.h" +#include "DDML/FastMLShower.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() {} ; + + virtual ~PionCloudsModel(){}; + + /// 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 + // These grid sizes were used for the two angle BIBAE + int m_numPoints = 2600; + 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) + int m_nLayer = 78; + + struct Vector3d{double x; double y; double z;}; + + Vector3d crossProduct(const Vector3d& v1, const Vector3d& v2); + + Vector3d normalize(const Vector3d& v); + + 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 6a23c1b..aa8cedb 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 { @@ -55,7 +60,7 @@ class PolyhedraBarrelGeometry : public GeometryInterface { std::string m_detector = {"EcalBarrel"}; int m_nSymmetry = 8; bool m_correctForAngles = false; - bool m_isHadShower = false; + bool m_isHadShower = true; std::string m_HadDetector = {"HcalBarrel"}; int m_nHadSymmetry = m_nSymmetry; }; diff --git a/scripts/ddsim_steer.py b/scripts/ddsim_steer.py index 71b192e..23450fb 100644 --- a/scripts/ddsim_steer.py +++ b/scripts/ddsim_steer.py @@ -574,6 +574,7 @@ def LoadHdf5(kernel): seq.adopt(sensitives) # ----------------- + ''' ## EM in Barrel model = DetectorConstruction(kernel, str(ml_model)) @@ -596,6 +597,7 @@ def LoadHdf5(kernel): model.enableUI() seq.adopt(model) + ''' # ------------------- ## EM in Endcap model1 = DetectorConstruction(kernel, str(ml_model_1)) @@ -625,7 +627,7 @@ def LoadHdf5(kernel): ## # Mandatory model parameters modelHad1.isHadShower = True - modelHad1.RegionName = "EcalBarrelRegion" ## hadron model triggers in ecal + 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 diff --git a/scripts/test_onnx.mac b/scripts/test_onnx.mac index 10f6964..3d5a716 100644 --- a/scripts/test_onnx.mac +++ b/scripts/test_onnx.mac @@ -3,7 +3,7 @@ # example script to run inference on a GAN w/ ONNX runtime # --- the particle gun -/gps/pos/centre 0 0 0 +/gps/pos/centre 30 180.48 0 #120 180.48 120 #30 180.48 0 #30 0 0 #/gps/particle gamma /gps/particle pi+ /gps/energy 50 GeV @@ -64,5 +64,5 @@ #/gps/direction 0.1 -0.8 .1 -/run/beamOn 10 #100 +/run/beamOn 2000 #100 diff --git a/src/Geant4FastHitMakerGlobal.cc b/src/Geant4FastHitMakerGlobal.cc index 9740761..bd7108f 100644 --- a/src/Geant4FastHitMakerGlobal.cc +++ b/src/Geant4FastHitMakerGlobal.cc @@ -48,6 +48,9 @@ Geant4FastHitMakerGlobal::~Geant4FastHitMakerGlobal() { //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void Geant4FastHitMakerGlobal::make(const G4FastHit& aHit, const G4FastTrack& aTrack) { + + std::cout << "Geant4FastHitMakerGlobal::make - aHit.GetPosition().x() " << aHit.GetPosition().x() << " aHit.GetPosition().y() " << + aHit.GetPosition().y() << " aHit.GetPosition().z() " << aHit.GetPosition().z() << std::endl; // do not make empty deposit if (aHit.GetEnergy() <= 0) { return; @@ -80,9 +83,13 @@ void Geant4FastHitMakerGlobal::make(const G4FastHit& aHit, const G4FastTrack& aT G4VSensitiveDetector* sensitive; if (currentVolume != 0) { + //std::cout << "IN FULL Sim Sensitive Geant4FastHitMakerGlobal::make - aHit.GetPosition().x() " << aHit.GetPosition().x() << " aHit.GetPosition().y() " << + // aHit.GetPosition().y() << " aHit.GetPosition().z() " << aHit.GetPosition().z() << " Energy: " << aHit.GetEnergy() << std::endl; sensitive = currentVolume->GetLogicalVolume()->GetSensitiveDetector(); G4VFastSimSensitiveDetector* fastSimSensitive = dynamic_cast(sensitive); if (fastSimSensitive) { + std::cout << "IN FAST Sim Sensitive Geant4FastHitMakerGlobal::make - aHit.GetPosition().x() " << aHit.GetPosition().x() << " aHit.GetPosition().y() " << + aHit.GetPosition().y() << " aHit.GetPosition().z() " << aHit.GetPosition().z() << " Energy: " << aHit.GetEnergy() << std::endl; fastSimSensitive->Hit(&aHit, &aTrack, &m_touchableHandle); } else if (sensitive && currentVolume->GetLogicalVolume()->GetFastSimulationManager()) { G4cerr << "ERROR - Geant4FastHitMakerGlobal::make()" << G4endl << " It is required to derive from the " diff --git a/src/PionCloudsModel.cc b/src/PionCloudsModel.cc new file mode 100644 index 0000000..6b04e9e --- /dev/null +++ b/src/PionCloudsModel.cc @@ -0,0 +1,142 @@ +#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_ = rotZ * (rotX * direction) ; + + 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()); + // std::cout << "PionClouds::localDir:" << "(" << localDir_.x() << "," << localDir_.y() << "," << localDir_.z() << ")" + // << std::endl; + + // 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); + // std::cout << "PionClouds::localDir:" << " theta = " << theta << " phi = " << phi << std::endl; + + // 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; + + // if (DEBUGPRINT) { + // std::cout << " Input_energy_tensor : " << inputs[0][0] << std::endl; + // std::cout << " Input_theta_tensor : " << inputs[1][0] << std::endl; + // std::cout << " Input_phi_tensor : " << inputs[2][0] << std::endl; + // } + + 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); +} + + + + void PionCloudsModel::convertOutput(G4FastTrack const& aFastTrack, + G4ThreeVector const& localDir, + const std::vector& output, + std::vector& spacepoints ){ + + int nPoints = m_numPoints ; // number of points in shower + + dd4hep::printout(dd4hep::DEBUG, "PionCloudsModel::convertOutput", "m_numPoints : %i", m_numPoints); + + spacepoints.resize( m_nLayer ) ; + + int numSP = 0; + int layerNum = 0; + int numElements = 0; + for (int i = 0; i < m_nLayer; i++) { + numSP = output[i] + 1; + spacepoints[i].reserve(numSP); + numElements += output[i] * 4; + } + + std::cout << " PionCloudsModel::convertOutput DONE numElements, numElements = " << numElements << std::endl; + + /* + for (int i = m_nLayer; i < m_nLayer + numElements; i += 4) { + std::cout << " PionCloudsModel::convertOutput Layer Loop, i = " << i << std::endl; + std::cout << " PionCloudsModel::convertOutput Layer Loop, output[i] " << output[i] << std::endl; + std::cout << " PionCloudsModel::convertOutput Layer Loop, output[i+1] " << output[i] << std::endl; + std::cout << " PionCloudsModel::convertOutput Layer Loop, output[i+3] " << output[i] << std::endl; + ddml::SpacePoint sp( + output[i], // x // *(-1) to align local to global convention in ddml + output[i+1], // y // *(-1) to align local to global convention in ddml + 0., // z + output[i+3], // energy + 0. // time + ); + layerNum = output[i+2]; + std::cout << " PionCloudsModel::convertOutput Layer Loop, layerNum " << layerNum << std::endl; + spacepoints[layerNum].emplace_back( sp ) ; + */ + + + float reshaped[2600][4]; + + // Fill the 3D array using the flattened vector + size_t index = 0; + for (size_t j = 0; j < 2600; ++j) { + for (size_t k = 0; k < 4; ++k) { + reshaped[j][k] = output[index++]; + } + } + + for (int i = 0; i < nPoints; i++) { + ddml::SpacePoint sp( + reshaped[i][0]*1e3, // x // *(-1) to align local to global convention in ddml + reshaped[i][2]*1e3, // y // *(-1) to align local to global convention in ddml + 0., // z + reshaped[i][3], // energy + 0. // time + ); + layerNum = reshaped[i][1]; + std::cout << "PionCloudsModel::convertOutput - layerNum" << layerNum <layers) { m_caloLayerDistances.push_back((l.distance + l.inner_thickness) / dd4hep::mm); + std::cout << " ECAL Layer distances " << l.distance + l.inner_thickness << std::endl; } } else { std::cout << " ###### error: detector " << m_detector << " not found !" << std::endl; @@ -31,6 +32,7 @@ void PolyhedraBarrelGeometry::initialize() { if (cal_had) { for (auto l_had : cal_had->layers) { m_caloLayerDistances.push_back((l_had.distance + l_had.inner_thickness) / dd4hep::mm); + std::cout << " HCAL Layer distances " << l_had.distance + l_had.inner_thickness << std::endl; } } else { @@ -138,6 +140,8 @@ void PolyhedraBarrelGeometry::localToGlobal(G4FastTrack const& aFastTrack, int firstLayer = 0; int nLayer = m_caloLayerDistances.size(); + std::cout << " PolyhedraBarrelGeometry::localToGlobal nLayer " << nLayer << std::endl; + for (int l = 0; l < nLayer; ++l) { double r = m_caloLayerDistances[l]; firstLayer = l; @@ -181,6 +185,10 @@ void PolyhedraBarrelGeometry::localToGlobal(G4FastTrack const& aFastTrack, sp.X = global.x(); sp.Y = global.y(); sp.Z = global.z(); + + std::cout << " PolyhedraBarrelGeometry::localToGlobal - global.x() = " << global.x() << " global.y() = " << global.y() + << " global.z() " << global.z() << " Energy: " << sp.E << std::endl; + } } } From dc199587fed34864759333cf47b9c0acbd4ceef2 Mon Sep 17 00:00:00 2001 From: Peter McKeown Date: Wed, 19 Feb 2025 09:30:31 +0100 Subject: [PATCH 04/12] Add option to record information about MC particle entry into calorimter --- CMakeLists.txt | 2 + include/DDML/DDMLEventAction.h | 70 +++ include/DDML/DDMLRunAction.h | 39 ++ include/DDML/FastMLShower.h | 45 ++ scripts/ddsim_steer.py | 2 +- scripts/ddsim_steer_Record_Calo_entry.py | 656 +++++++++++++++++++++++ scripts/test_onnx.mac | 4 +- src/CMakeLists.txt | 16 + src/DDMLEventAction.cc | 58 ++ src/DDMLRunAction.cc | 94 ++++ src/PionCloudsModel.cc | 4 +- 11 files changed, 985 insertions(+), 5 deletions(-) create mode 100644 include/DDML/DDMLEventAction.h create mode 100644 include/DDML/DDMLRunAction.h create mode 100644 scripts/ddsim_steer_Record_Calo_entry.py create mode 100644 src/DDMLEventAction.cc create mode 100644 src/DDMLRunAction.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index d6dceb3..2b58d4b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -58,6 +58,8 @@ endif() option(INSTRUMENT_MODEL "Instrument the steps of the showerModel call" OFF) +option(RECORD_CALO_IMPACT "Record properties of particle that triggers modelShower call" OFF) + option(DOWNLOAD_MODELS "Download and install the models that are stored externally" ON) #--------------------------------------------------- diff --git a/include/DDML/DDMLEventAction.h b/include/DDML/DDMLEventAction.h new file mode 100644 index 0000000..1aa857c --- /dev/null +++ b/include/DDML/DDMLEventAction.h @@ -0,0 +1,70 @@ +#include // for G4int, G4double +#include // for vector +#include "G4UserEventAction.hh" // for G4UserEventAction +//#include "G4Timer.hh" // for G4Timer +class G4Event; +#include "DDG4/Geant4Handle.h" +#include "DDG4/Geant4Kernel.h" +#include "DDG4/Geant4EventAction.h" +#include "CLHEP/Units/SystemOfUnits.h" +#include "CLHEP/Units/PhysicalConstants.h" + + +namespace ddml { + +/** Event action for ddml. +* +* @author P. McKeown, CERN +* @date Feb 2025 +* +*/ + +class DDMLEventAction: public dd4hep::sim::Geant4EventAction{ +public: + /// Standard constructor with initializing arguments + DDMLEventAction(dd4hep::sim::Geant4Context* c, const std::string& n); + /// Default destructor + virtual ~DDMLEventAction(); + /// begin-of-event callback + inline virtual void begin(const G4Event*) override; + /// End-of-event callback + virtual void end(const G4Event*) override; + /// begin-of-run callback + void beginRun(const G4Run*); + /// End-of-run callback + void endRun(const G4Run*); + + //// Get and Set methods for Calo Face info + inline std::vector& GetCaloMC_PDG() {return m_CaloMCPDG;} + inline std::vector& GetCaloMC_E() {return m_CaloMCE;} + inline std::vector& GetCaloMC_PosX() {return m_CaloMCPosX;} + inline std::vector& GetCaloMC_PosY() {return m_CaloMCPosY;} + inline std::vector& GetCaloMC_PosZ() {return m_CaloMCPosZ;} + inline std::vector& GetCaloMC_DirX() {return m_CaloMCDirX;} + inline std::vector& GetCaloMC_DirY() {return m_CaloMCDirY;} + inline std::vector& GetCaloMC_DirZ() {return m_CaloMCDirZ;} + + // Set methods to push back vector + inline void SetElCaloMC_PDG(G4int aValue) {m_CaloMCPDG.push_back(aValue);} + inline void SetElCaloMC_E(G4double aValue) {m_CaloMCE.push_back(aValue);} + inline void SetElCaloMC_PosX(G4double aValue) {m_CaloMCPosX.push_back(aValue);} + inline void SetElCaloMC_PosY(G4double aValue) {m_CaloMCPosY.push_back(aValue);} + inline void SetElCaloMC_PosZ(G4double aValue) {m_CaloMCPosZ.push_back(aValue);} + inline void SetElCaloMC_DirX(G4double aValue) {m_CaloMCDirX.push_back(aValue);} + inline void SetElCaloMC_DirY(G4double aValue) {m_CaloMCDirY.push_back(aValue);} + inline void SetElCaloMC_DirZ(G4double aValue) {m_CaloMCDirZ.push_back(aValue);} + +private: + // Fast Sim Calo entrace particle properties to store in ntuple + std::vector m_CaloMCPDG; + std::vector m_CaloMCE; + std::vector m_CaloMCPosX; + std::vector m_CaloMCPosY; + std::vector m_CaloMCPosZ; + std::vector m_CaloMCDirX; + std::vector m_CaloMCDirY; + std::vector m_CaloMCDirZ; +}; + +} //namespace + diff --git a/include/DDML/DDMLRunAction.h b/include/DDML/DDMLRunAction.h new file mode 100644 index 0000000..49c5dc6 --- /dev/null +++ b/include/DDML/DDMLRunAction.h @@ -0,0 +1,39 @@ +#include // for GeV +#include // for G4String +#include // for G4ThreeVector +#include // for G4int +//#include "G4Timer.hh" // for G4Timer +class G4ParticleDefinition; +#include "DDG4/Geant4Handle.h" +#include "DDG4/Geant4Kernel.h" +#include "DDG4/Geant4RunAction.h" + + + +namespace ddml{ +/* Create analysis files in standard G4 manner +* +* @author P. McKeown +* @date Feb 2025 +* +*/ + +class DDMLRunAction: public dd4hep::sim::Geant4RunAction { + public: + DDMLRunAction() = delete; + /// Standard constructor with initializing arguments + DDMLRunAction(dd4hep::sim::Geant4Context* c, const std::string& n) ; + /// Default destructor + virtual ~DDMLRunAction() ; + /// begin-of-run callback + void begin(const G4Run*) override; + /// End-of-run callback + void end(const G4Run*) override; + /// begin-of-event callback + void beginEvent(const G4Event*) ; + /// End-of-event callback + void endEvent(const G4Event*) ; + +}; + +} // namespace \ No newline at end of file diff --git a/include/DDML/FastMLShower.h b/include/DDML/FastMLShower.h index e48a35a..0615eb2 100644 --- a/include/DDML/FastMLShower.h +++ b/include/DDML/FastMLShower.h @@ -14,6 +14,10 @@ #define DDML_INSTRUMENT_MODEL_SHOWER 0 #endif +#ifndef DDML_RECORD_CALO_IMPACT + #define DDML_RECORD_CALO_IMPACT 0 +#endif + #if DDML_INSTRUMENT_MODEL_SHOWER #include "podio/Frame.h" #include "podio/ROOTFrameWriter.h" @@ -35,6 +39,12 @@ inline auto run_void_member_timed(Obj& obj, MemberFunc func, Args&&... args) { } #endif +#if DDML_RECORD_CALO_IMPACT + //#include "G4EventManager.hh" + //#include "DDML/DDMLEventAction.h" + #include "G4AnalysisManager.hh" +#endif + namespace ddml { /** The templated base class for running fast shower simulation with ML. * The actual implementation is provided by the templated class @@ -144,6 +154,41 @@ class FastMLShower : public dd4hep::sim::Geant4FastSimShowerModel { podio::UserDataCollection nHits; #endif +#if DDML_RECORD_CALO_IMPACT + // Create analysis manager + G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); + + // Store information about particle impact on calorimter face from track + G4int PDG = track.GetPrimaryTrack()->GetParticleDefinition()->GetPDGEncoding(); + G4ThreeVector position = track.GetPrimaryTrack()->GetPosition(); + G4ThreeVector direction = track.GetPrimaryTrack()->GetMomentumDirection(); + + // Fill analysis manager + analysisManager->FillNtupleIColumn(1, 0, PDG); + analysisManager->FillNtupleDColumn(1, 1, energy); + analysisManager->FillNtupleDColumn(1, 2, position.x()); + analysisManager->FillNtupleDColumn(1, 3, position.y()); + analysisManager->FillNtupleDColumn(1, 4, position.z()); + analysisManager->FillNtupleDColumn(1, 5, direction.x()); + analysisManager->FillNtupleDColumn(1, 6, direction.y()); + analysisManager->FillNtupleDColumn(1, 7, direction.z()); + + analysisManager->AddNtupleRow(1); + + /* + DDMLEventAction* mEventAction = dynamic_cast(G4EventManager::GetEventManager()->GetUserEventAction()); + mEventAction->SetElCaloMC_PDG(PDG); + mEventAction->SetElCaloMC_E(energy); + mEventAction->SetElCaloMC_PosX(position.x()); + mEventAction->SetElCaloMC_PosY(position.y()); + mEventAction->SetElCaloMC_PosZ(position.z()); + mEventAction->SetElCaloMC_DirX(direction.x()); + mEventAction->SetElCaloMC_DirY(direction.y()); + mEventAction->SetElCaloMC_DirZ(direction.z()); + */ + +#endif + for (auto& invec : m_input) { invec.clear(); } diff --git a/scripts/ddsim_steer.py b/scripts/ddsim_steer.py index 23450fb..9c7f8a0 100644 --- a/scripts/ddsim_steer.py +++ b/scripts/ddsim_steer.py @@ -547,7 +547,7 @@ def LoadHdf5(kernel): if hadrons == True: ml_model_had = "LoadHDF5PionCloudsPCHadronModelPolyhedraBarrel/BarrelModelTorch" - ml_had_file = "../models/PionClouds_50GeV_sp.h5" + ml_had_file = "../models/PionClouds_50GeV_sp.h5" #"../models/PionClouds_50GeV_sp_scaled.h5" ml_correct_angles = False from g4units import GeV, MeV # DO NOT REMOVE OR MOVE!!!!! (EXCLAMATION MARK) diff --git a/scripts/ddsim_steer_Record_Calo_entry.py b/scripts/ddsim_steer_Record_Calo_entry.py new file mode 100644 index 0000000..6107367 --- /dev/null +++ b/scripts/ddsim_steer_Record_Calo_entry.py @@ -0,0 +1,656 @@ +###################################################################### +# +# standard steering file for ILD simulation +# +# +# +###################################################################### +from DDSim.DD4hepSimulation import DD4hepSimulation +import DDG4 +from g4units import m, mm, GeV, MeV, rad +import os + +SIM = DD4hepSimulation() + +## The compact XML file +SIM.compactFile = "" +## Lorentz boost for the crossing angle, in radian! +SIM.crossingAngleBoost = 7.0e-3 * rad +SIM.enableDetailedShowerMode = True +SIM.enableG4GPS = True +SIM.enableG4Gun = False +SIM.enableGun = False +## InputFiles for simulation .stdhep, .slcio, .HEPEvt, .hepevt, .hepmc, .pairs files are supported +SIM.inputFiles = [] +## Macro file to execute for runType 'run' or 'vis' +SIM.macroFile = "./test_onnx.mac" +## number of events to simulate, used in batch mode +SIM.numberOfEvents = 100 +## Outputfile from the simulation,only lcio output is supported +# SIM.outputFile = "dummyOutput_edm4hep.root" ##"dummyOutput.slcio" +SIM.outputFile = "dummyOutput.slcio" + +## Physics list to use in simulation +SIM.physicsList = None +## Verbosity use integers from 1(most) to 7(least) verbose +## or strings: VERBOSE, DEBUG, INFO, WARNING, ERROR, FATAL, ALWAYS +SIM.printLevel = "INFO" +## The type of action to do in this invocation +## batch: just simulate some events, needs numberOfEvents, and input file or gun +## vis: enable visualisation, run the macroFile if it is set +## run: run the macroFile and exit +## shell: enable interactive session +SIM.runType = "run" # "batch" +## Skip first N events when reading a file +SIM.skipNEvents = 0 +## Steering file to change default behaviour +SIM.steeringFile = None +## FourVector of translation for the Smearing of the Vertex position: x y z t +SIM.vertexOffset = [0.0, 0.0, 0.0, 0.0] +## FourVector of the Sigma for the Smearing of the Vertex position: x y z t +SIM.vertexSigma = [0.0, 0.0, 0.0, 0.0] + + +################################################################################ +## Action holding sensitive detector actions +## The default tracker and calorimeter actions can be set with +## +## >>> SIM = DD4hepSimulation() +## >>> SIM.action.tracker = ('Geant4TrackerWeightedAction', {'HitPositionCombination': 2, 'CollectSingleDeposits': False}) +## >>> SIM.action.calo = "Geant4CalorimeterAction" +## +## for specific subdetectors specific sensitive detectors can be set based on pattern matching +## +## >>> SIM = DD4hepSimulation() +## >>> SIM.action.mapActions['tpc'] = "TPCSDAction" +## +## and additional parameters for the sensitive detectors can be set when the map is given a tuple +## +## >>> SIM = DD4hepSimulation() +## >>> SIM.action.mapActions['ecal'] =( "CaloPreShowerSDAction", {"FirstLayerNumber": 1} ) +## +## +################################################################################ + +## set the default calorimeter action +SIM.action.calo = "Geant4ScintillatorCalorimeterAction" + +## create a map of patterns and actions to be applied to sensitive detectors +## example: SIM.action.mapActions['tpc'] = "TPCSDAction" +SIM.action.mapActions = {} + +SIM.action.mapActions["tpc"] = "TPCSDAction" + +## set the default tracker action +SIM.action.tracker = ( + "Geant4TrackerWeightedAction", + {"HitPositionCombination": 2, "CollectSingleDeposits": False}, +) + +### Configure Run actions +kernel = DDG4.Kernel() +run1 = DDG4.RunAction(kernel, 'DDMLRunAction/runaction') +kernel.registerGlobalAction(run1) +kernel.runAction().add(run1) +event1 = DDG4.EventAction(kernel, 'DDMLEventAction/eventaction') +kernel.registerGlobalAction(event1) +kernel.eventAction().add(event1) + +################################################################################ +## Configuration for the magnetic field (stepper) +################################################################################ +## --- used in v01-19-05 : +SIM.field.delta_chord = 1e-05 +SIM.field.delta_intersection = 1e-05 +SIM.field.delta_one_step = 0.5e-03 * mm +SIM.field.eps_max = 1e-04 +SIM.field.eps_min = 1e-05 +SIM.field.equation = "Mag_UsualEqRhs" +SIM.field.largest_step = 10.0 * m +SIM.field.min_chord_step = 1.0e-2 * mm +SIM.field.stepper = "HelixSimpleRunge" + +## --- default values in ddsim +##SIM.field.delta_chord = 0.25 +##SIM.field.delta_intersection = 0.001 +##SIM.field.delta_one_step = 0.01 +##SIM.field.eps_max = 0.001 +##SIM.field.eps_min = 5e-05 +##SIM.field.equation = "Mag_UsualEqRhs" +##SIM.field.largest_step = 10000.0 +##SIM.field.min_chord_step = 0.01 +##SIM.field.stepper = "G4ClassicalRK4" + +################################################################################ +## Configuration for sensitive detector filters +## +## Set the default filter for tracker or caliromter +## >>> SIM.filter.tracker = "edep1kev" +## >>> SIM.filter.calo = "" +## +## Assign a filter to a sensitive detector via pattern matching +## >>> SIM.filter.mapDetFilter['FTD'] = "edep1kev" +## +## Or more than one filter: +## >>> SIM.filter.mapDetFilter['FTD'] = ["edep1kev", "geantino"] +## +## Don't use the default filter or anything else: +## >>> SIM.filter.mapDetFilter['TPC'] = None ## or "" or [] +## +## Create a custom filter. The dictionary is used to instantiate the filter later on +## >>> SIM.filter.filters['edep3kev'] = dict(name="EnergyDepositMinimumCut/3keV", parameter={"Cut": 3.0*keV} ) +## +## +################################################################################ + +## default filter for calorimeter sensitive detectors; this is applied if no other filter is used for a calorimeter +SIM.filter.calo = "edep0" + +## list of filter objects: map between name and parameter dictionary +SIM.filter.filters = { + "edep0": {"parameter": {"Cut": 0.0}, "name": "EnergyDepositMinimumCut/Cut0"}, + "geantino": {"parameter": {}, "name": "GeantinoRejectFilter/GeantinoRejector"}, + "edep1kev": {"parameter": {"Cut": 0.001}, "name": "EnergyDepositMinimumCut"}, +} + +## a map between patterns and filter objects, using patterns to attach filters to sensitive detector +SIM.filter.mapDetFilter = {} + +SIM.filter.mapDetFilter["TPC"] = None + +## default filter for tracking sensitive detectors; this is applied if no other filter is used for a tracker +SIM.filter.tracker = "edep1kev" + + +################################################################################ +## Configuration for the GuineaPig InputFiles +################################################################################ + +## Set the number of pair particles to simulate per event. +## Only used if inputFile ends with ".pairs" +## If "-1" all particles will be simulated in a single event +## +SIM.guineapig.particlesPerEvent = "-1" + + +################################################################################ +## Configuration for the DDG4 ParticleGun +################################################################################ + +## direction of the particle gun, 3 vector +SIM.gun.direction = (0, 0, 1) + +## choose the distribution of the random direction for theta +## +## Options for random distributions: +## +## 'uniform' is the default distribution, flat in theta +## 'cos(theta)' is flat in cos(theta) +## 'eta', or 'pseudorapidity' is flat in pseudorapity +## 'ffbar' is distributed according to 1+cos^2(theta) +## +## Setting a distribution will set isotrop = True +## +SIM.gun.distribution = None +SIM.gun.energy = 10000.0 + +## isotropic distribution for the particle gun +## +## use the options phiMin, phiMax, thetaMin, and thetaMax to limit the range of randomly distributed directions +## if one of these options is not None the random distribution will be set to True and cannot be turned off! +## +SIM.gun.isotrop = False +SIM.gun.multiplicity = 1 +SIM.gun.particle = "mu-" +SIM.gun.phiMax = None + +## Minimal azimuthal angle for random distribution +SIM.gun.phiMin = None + +## position of the particle gun, 3 vector +SIM.gun.position = (0.0, 0.0, 0.0) +SIM.gun.thetaMax = None +SIM.gun.thetaMin = None + + +################################################################################ +## Configuration for the output levels of DDG4 components +################################################################################ + +## Output level for input sources +SIM.output.inputStage = 3 + +## Output level for Geant4 kernel +SIM.output.kernel = 3 + +## Output level for ParticleHandler +SIM.output.part = 3 + +## Output level for Random Number Generator setup +SIM.output.random = 6 + + +################################################################################ +## Configuration for the Particle Handler/ MCTruth treatment +################################################################################ + +## Enable lots of printout on simulated hits and MC-truth information +SIM.part.enableDetailedHitsAndParticleInfo = False + +## Keep all created particles +SIM.part.keepAllParticles = False + +## Minimal distance between particle vertex and endpoint of parent after +## which the vertexIsNotEndpointOfParent flag is set +## +SIM.part.minDistToParentVertex = 2.2e-14 + +## MinimalKineticEnergy to store particles created in the tracking region +SIM.part.minimalKineticEnergy = 1 * MeV + +## Printout at End of Tracking +SIM.part.printEndTracking = False + +## Printout at Start of Tracking +SIM.part.printStartTracking = False + +## List of processes to save, on command line give as whitespace separated string in quotation marks +SIM.part.saveProcesses = ["Decay"] + + +################################################################################ +## Configuration for the PhysicsList +################################################################################ +# this needs to be set to False if any standard physics list is used: +SIM.physics.decays = False +SIM.physics.list = "QGSP_BERT" # "FTFP_BERT" + +## location of particle.tbl file containing extra particles and their lifetime information +## +SIM.physics.pdgfile = os.path.join( + os.environ.get("DD4HEP"), "DDG4/examples/particle.tbl" +) + +## The global geant4 rangecut for secondary production +## +## Default is 0.7 mm as is the case in geant4 10 +## +## To disable this plugin and be absolutely sure to use the Geant4 default range cut use "None" +## +## Set printlevel to DEBUG to see a printout of all range cuts, +## but this only works if range cut is not "None" +## +SIM.physics.rangecut = 0.1 * mm + + +################################################################################ +## Properties for the random number generator +################################################################################ + +## If True, calculate random seed for each event based on eventID and runID +## allows reproducibility even when SkippingEvents +SIM.random.enableEventSeed = True #False +SIM.random.file = None +SIM.random.luxury = 1 +SIM.random.replace_gRandom = True +SIM.random.seed = 42 #None +SIM.random.type = None + +# --------------------------------------------- +# +# Configure ML inference +# +# --------------------------------------------- + + +def aiDance(kernel): + ild = True + par04 = True # False + old_DD4hep = False ## use for DD4hep versions/commits before ~ Apr 21st 2023 + + if ild == True: + ml_barrel_name = "EcalBarrel" + ml_barrel_symmetry = 8 + ml_endcap_name = "EcalEndcap" + else: + ml_barrel_name = "ECalBarrel" + ml_barrel_symmetry = 12 + ml_endcap_name = "ECalEndcap" + + if par04 == True: + ml_file = "../models/Generator.onnx" + ml_model = "Par04ExampleVAEPolyhedraBarrelONNXModel/ShowerModel" + ml_model1 = "Par04ExampleVAEEndcapONNXModel/ShowerModel" + else: + ml_file = "../models/francisca_gan.onnx" + ml_model = "RegularGridGANPolyhedraBarrelONNXModel/ShowerModel" + ml_model1 = "RegularGridGANEndcapONNXModel/ShowerModel" + + ml_correct_angles = True + + from g4units import GeV, MeV # DO NOT REMOVE OR MOVE!!!!! (EXCLAMATION MARK) + from DDG4 import DetectorConstruction, Geant4, PhysicsList + + geant4 = Geant4(kernel) + seq = geant4.detectorConstruction() + + if old_DD4hep: # this is now done in DD4hepSimulations.py, i.e. in ddsim + seq, act = geant4.addDetectorConstruction( + "Geant4DetectorGeometryConstruction/ConstructGeo" + ) + act.DebugMaterials = True + act.DebugElements = False + act.DebugVolumes = True + act.DebugShapes = True + # Apply sensitive detectors + sensitives = DetectorConstruction( + kernel, str("Geant4DetectorSensitivesConstruction/ConstructSD") + ) + sensitives.enableUI() + seq.adopt(sensitives) + + # ----------------- + model = DetectorConstruction(kernel, str(ml_model)) + + ## # Mandatory model parameters + model.RegionName = "EcalBarrelRegion" + model.Detector = ml_barrel_name + model.Symmetry = ml_barrel_symmetry + model.Enable = True + model.CorrectForAngles = ml_correct_angles + # Energy boundaries are optional: Units are GeV + model.ApplicableParticles = {"e+", "e-", "gamma"} + model.Etrigger = {"e+": 5.0 * GeV, "e-": 5.0 * GeV, "gamma": 5.0 * GeV} + model.ModelPath = ml_file + model.OptimizeFlag = 1 + + model.enableUI() + seq.adopt(model) + # ------------------- + model1 = DetectorConstruction(kernel, str(ml_model1)) + + ## # Mandatory model parameters + model1.RegionName = "EcalEndcapRegion" + model1.Detector = ml_endcap_name + model1.Enable = True + model1.CorrectForAngles = ml_correct_angles + # Energy boundaries are optional: Units are GeV + model1.ApplicableParticles = {"e+", "e-", "gamma"} + model1.Etrigger = {"e+": 5.0 * GeV, "e-": 5.0 * GeV, "gamma": 5.0 * GeV} + model1.ModelPath = ml_file + model.OptimizeFlag = 1 + + model1.enableUI() + seq.adopt(model1) + # ------------------- + + # Now build the physics list: + phys = kernel.physicsList() + ph = PhysicsList(kernel, str("Geant4FastPhysics/FastPhysicsList")) + ph.EnabledParticles = ["e+", "e-", "gamma"] + ph.BeVerbose = True + ph.enableUI() + phys.adopt(ph) + phys.dump() + + +def aiDanceTorch(kernel): + ild = True + BIBAE = True + Two_Angle = True #True + old_DD4hep = False ## use for DD4hep versions/commits before ~ Apr 21st 2023 + + if ild == True: + ml_barrel_name = "EcalBarrel" + ml_barrel_symmetry = 8 + ml_endcap_name = "EcalEndcap" + else: + ml_barrel_name = "ECalBarrel" + ml_barrel_symmetry = 12 + ml_endcap_name = "ECalEndcap" + + if BIBAE == True and Two_Angle == False: + ml_file = "../models/BIBAE_Full_PP_cut.pt" + ml_model = "RegularGridBIBAEPolyhedraBarrelTorchModel/BarrelModelTorch" + ml_model_1 = "RegularGridBIBAEEndcapTorchModel/EndcapModelTorch" + ml_correct_angles = False + elif BIBAE == True and Two_Angle == True: + ml_file = "../models/BIBAE_Two_Angle_Full_PP_cut.pt" + ml_model = ( + "RegularGridTwoAngleBIBAEModelPolyhedraBarrelTorchModel/BarrelModelTorch" + ) + ml_model_1 = "RegularGridTwoAngleBIBAEModelEndcapTorchModel/EndcapModelTorch" + ml_correct_angles = False + else: + ml_file = "../models/francisca_gan_jit.pt" + ml_model = "RegularGridGANPolyhedraBarrelTorchModel/BarrelModelTorch" + ml_model_1 = "RegularGridGANEndcapTorchModel/EndcapModelTorch" + ml_correct_angles = True + + from g4units import GeV, MeV # DO NOT REMOVE OR MOVE!!!!! (EXCLAMATION MARK) + from DDG4 import DetectorConstruction, Geant4, PhysicsList + + geant4 = Geant4(kernel) + + seq = geant4.detectorConstruction() + + if old_DD4hep: # this is now done in DD4hepSimulations.py, i.e. in ddsim + seq, act = geant4.addDetectorConstruction( + "Geant4DetectorGeometryConstruction/ConstructGeo" + ) + act.DebugMaterials = True + act.DebugElements = False + act.DebugVolumes = True + act.DebugShapes = True + + # Apply sensitive detectors + sensitives = DetectorConstruction( + kernel, str("Geant4DetectorSensitivesConstruction/ConstructSD") + ) + sensitives.enableUI() + seq.adopt(sensitives) + + # ----------------- + model = DetectorConstruction(kernel, str(ml_model)) + + ## # Mandatory model parameters + model.RegionName = "EcalBarrelRegion" + model.Detector = ml_barrel_name + model.Symmetry = ml_barrel_symmetry + model.Enable = True + model.CorrectForAngles = ml_correct_angles + # Energy boundaries are optional: Units are GeV + model.ApplicableParticles = {"e+", "e-", "gamma"} + model.Etrigger = { + "e+": 10.0 * GeV, + "e-": 10.0 * GeV, + "gamma": 10.0 * GeV, + } # trigger on lower training threshold + model.ModelPath = ml_file + model.OptimizeFlag = 1 + model.IntraOpNumThreads = 1 + + model.enableUI() + seq.adopt(model) + # ------------------- + model1 = DetectorConstruction(kernel, str(ml_model_1)) + + ## # Mandatory model parameters + model1.RegionName = "EcalEndcapRegion" + model1.Detector = ml_endcap_name + model1.Enable = True + model1.CorrectForAngles = ml_correct_angles + # Energy boundaries are optional: Units are GeV + model1.ApplicableParticles = {"e+", "e-", "gamma"} + model1.Etrigger = { + "e+": 10.0 * GeV, + "e-": 10.0 * GeV, + "gamma": 10.0 * GeV, + } # trigger on lower training threshold + model1.ModelPath = ml_file + model1.OptimizeFlag = 1 + model1.IntraOpNumThreads = 1 + + model1.enableUI() + seq.adopt(model1) + # ------------------- + + # Now build the physics list: + phys = kernel.physicsList() + ph = PhysicsList(kernel, str("Geant4FastPhysics/FastPhysicsList")) + ph.EnabledParticles = ["e+", "e-", "gamma"] + ph.BeVerbose = True + ph.enableUI() + phys.adopt(ph) + phys.dump() + + +def LoadHdf5(kernel): + ild = True + 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 = ( + "LoadHDF5RegularGridTwoAngleBIBAEModelPolyhedraBarrel/BarrelModelTorch" + ) + ml_model_1 = "LoadHDF5RegularGridTwoAngleBIBAEModelEndcap/EndcapModelTorch" + ml_correct_angles = False + + if hadrons == True: + ml_model_had = "LoadHDF5PionCloudsPCHadronModelPolyhedraBarrel/BarrelModelTorch" + ml_had_file = "../models/PionClouds_50GeV_sp.h5" #"../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 + + geant4 = Geant4(kernel) + + seq = geant4.detectorConstruction() + + if old_DD4hep: # this is now done in DD4hepSimulations.py, i.e. in ddsim + seq, act = geant4.addDetectorConstruction( + "Geant4DetectorGeometryConstruction/ConstructGeo" + ) + act.DebugMaterials = True + act.DebugElements = False + act.DebugVolumes = True + act.DebugShapes = True + + # Apply sensitive detectors + sensitives = DetectorConstruction( + kernel, str("Geant4DetectorSensitivesConstruction/ConstructSD") + ) + sensitives.enableUI() + seq.adopt(sensitives) + + # ----------------- + ''' + ## EM in Barrel + model = DetectorConstruction(kernel, str(ml_model)) + + ## # Mandatory model parameters + model.RegionName = "EcalBarrelRegion" + model.Detector = ml_barrel_name + model.Symmetry = ml_barrel_symmetry + model.Enable = True + model.CorrectForAngles = ml_correct_angles + # Energy boundaries are optional: Units are GeV + model.ApplicableParticles = {"e+", "e-", "gamma"} + model.Etrigger = { + "e+": 10.0 * GeV, + "e-": 10.0 * GeV, + "gamma": 10.0 * GeV, + } # trigger on lower training threshold + model.FilePath = ml_file + # model.OptimizeFlag = 1 + # model.IntraOpNumThreads = 1 + + model.enableUI() + seq.adopt(model) + ''' + # ------------------- + ## EM in Endcap + model1 = DetectorConstruction(kernel, str(ml_model_1)) + + ## # Mandatory model parameters + model1.RegionName = "EcalEndcapRegion" + model1.Detector = ml_endcap_name + model1.Enable = True + model1.CorrectForAngles = ml_correct_angles + # Energy boundaries are optional: Units are GeV + model1.ApplicableParticles = {"e+", "e-", "gamma"} + model1.Etrigger = { + "e+": 10.0 * GeV, + "e-": 10.0 * GeV, + "gamma": 10.0 * GeV, + } # trigger on lower training threshold + model1.FilePath = ml_file + # model1.OptimizeFlag = 1 + # model1.IntraOpNumThreads = 1 + + 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", "pi+"] + ph.BeVerbose = True + ph.enableUI() + phys.adopt(ph) + phys.dump() + + +#SIM.physics.setupUserPhysics( aiDance) +#SIM.physics.setupUserPhysics(aiDanceTorch) +SIM.physics.setupUserPhysics(LoadHdf5) diff --git a/scripts/test_onnx.mac b/scripts/test_onnx.mac index 3d5a716..4226261 100644 --- a/scripts/test_onnx.mac +++ b/scripts/test_onnx.mac @@ -3,7 +3,7 @@ # example script to run inference on a GAN w/ ONNX runtime # --- the particle gun -/gps/pos/centre 30 180.48 0 #120 180.48 120 #30 180.48 0 #30 0 0 +/gps/pos/centre 30 0 0 #30 180.48 0 #### this is ILD calo Face #120 180.48 120 #30 180.48 0 #30 0 0 #/gps/particle gamma /gps/particle pi+ /gps/energy 50 GeV @@ -64,5 +64,5 @@ #/gps/direction 0.1 -0.8 .1 -/run/beamOn 2000 #100 +/run/beamOn 100 #2000 #100 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2882c79..0a0114c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -35,6 +35,18 @@ if(INSTRUMENT_MODEL) ) endif() +if(RECORD_CALO_IMPACT) + target_sources(${PROJECT_NAME} PRIVATE + DDMLRunAction.cc + DDMLEventAction.cc + ) + target_compile_definitions(${PROJECT_NAME} PUBLIC + DDML_RECORD_CALO_IMPACT=1 + ) +endif() + + + if(OnnxRuntime_FOUND) target_sources(${PROJECT_NAME} PRIVATE ONNXInference.cc @@ -111,6 +123,10 @@ if(HDF5_FOUND) list(APPEND headers ${PROJECT_SOURCE_DIR}/include/DDML/LoadHdf5.h) endif() +if(RECORD_CALO_IMPACT) + list(APPEND headers ${PROJECT_SOURCE_DIR}/include/DDML/DDMLRunAction.h ${PROJECT_SOURCE_DIR}/include/DDML/DDMLEventAction.h) +endif() + #install(FILES # ${headers} # DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) diff --git a/src/DDMLEventAction.cc b/src/DDMLEventAction.cc new file mode 100644 index 0000000..7decd72 --- /dev/null +++ b/src/DDMLEventAction.cc @@ -0,0 +1,58 @@ +#include "DDML/DDMLEventAction.h" +#include "G4AnalysisManager.hh" +#include "DD4hep/InstanceCount.h" +#include "G4Event.hh" +#include "G4EventManager.hh" +#include "G4ThreeVector.hh" +#include "G4PrimaryVertex.hh" +#include "G4PrimaryParticle.hh" +#include "DDG4/Geant4Data.h" + +#define DEBUGPRINT 0 + +namespace ddml { + +DDMLEventAction::DDMLEventAction(dd4hep::sim::Geant4Context* c, const std::string& n): +Geant4EventAction(c,n) + { + dd4hep::InstanceCount::increment(this); + } + +/// Default destructor +DDMLEventAction::~DDMLEventAction() { + dd4hep::InstanceCount::decrement(this); + } + +void DDMLEventAction::begin(const G4Event*){} + +void DDMLEventAction::end(const G4Event*) { + // Get analysis manager + auto analysisManager = G4AnalysisManager::Instance(); + + // Retrieve information from primary vertex and primary particle + auto primaryVertex = + G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetPrimaryVertex(); + G4ThreeVector primaryVertexPos = primaryVertex->GetPosition(); + auto primaryParticle = primaryVertex->GetPrimary(0); + G4double primaryEnergy = primaryParticle->GetTotalEnergy(); + G4ThreeVector primaryDirection = primaryParticle->GetMomentumDirection(); + G4int primaryPDG = primaryParticle->GetPDGcode(); + + // Fill analysis manager + analysisManager->FillNtupleIColumn(0, 0, primaryPDG); + analysisManager->FillNtupleDColumn(0, 1, primaryEnergy); + analysisManager->FillNtupleDColumn(0, 2, primaryVertexPos.x()); + analysisManager->FillNtupleDColumn(0, 3, primaryVertexPos.y()); + analysisManager->FillNtupleDColumn(0, 4, primaryVertexPos.z()); + analysisManager->FillNtupleDColumn(0, 5, primaryDirection.x()); + analysisManager->FillNtupleDColumn(0, 6, primaryDirection.y()); + analysisManager->FillNtupleDColumn(0, 7, primaryDirection.z()); + + analysisManager->AddNtupleRow(0); + +} + +} //namespace + +#include "DDG4/Factories.h" +DECLARE_GEANT4ACTION_NS(ddml, DDMLEventAction) diff --git a/src/DDMLRunAction.cc b/src/DDMLRunAction.cc new file mode 100644 index 0000000..c7cd807 --- /dev/null +++ b/src/DDMLRunAction.cc @@ -0,0 +1,94 @@ +#include "DDML/DDMLRunAction.h" +#include "DDML/DDMLEventAction.h" +#include "G4AnalysisManager.hh" +#include "G4EventManager.hh" +#include "DD4hep/InstanceCount.h" +#include "G4Run.hh" + +#define DEBUGPRINT 0 + +namespace ddml { + +/// Standard constructor with initializing arguments +DDMLRunAction::DDMLRunAction(dd4hep::sim::Geant4Context* c, const std::string& n): + Geant4RunAction(c, n) { + // Create analysis manager + G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); + analysisManager->SetDefaultFileType("root"); + + // Default filename, can be overriden with /analysis/setFileName + analysisManager->SetFileName("Output"); + dd4hep::InstanceCount::increment(this); +} + +/// Default destructor +DDMLRunAction::~DDMLRunAction() { + dd4hep::InstanceCount::decrement(this); +} + +/// begin-of-run callback +void DDMLRunAction::begin(const G4Run*) { + // Get analysis manager +G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); + +// Create directories +analysisManager->SetVerboseLevel(0); + +// Create ntuples +analysisManager->CreateNtuple("global", "Event data"); +analysisManager->CreateNtupleIColumn("MC_PDG"); +analysisManager->CreateNtupleDColumn("MC_Energy"); +analysisManager->CreateNtupleDColumn("MC_PosX"); +analysisManager->CreateNtupleDColumn("MC_PosY"); +analysisManager->CreateNtupleDColumn("MC_PosZ"); +analysisManager->CreateNtupleDColumn("MC_DirX"); +analysisManager->CreateNtupleDColumn("MC_DirY"); +analysisManager->CreateNtupleDColumn("MC_DirZ"); +analysisManager->FinishNtuple(); + +/* +DDMLEventAction* mEventAction = dynamic_cast(G4EventManager::GetEventManager()->GetUserEventAction()); +analysisManager->CreateNtuple("Fast Sim Info", "MC info at calo face"); +analysisManager->CreateNtupleIColumn("Calo_MC_PDG", mEventAction->GetCaloMC_PDG()); +analysisManager->CreateNtupleDColumn("Calo_MC_Energy", mEventAction->GetCaloMC_E()); +analysisManager->CreateNtupleDColumn("Calo_MC_PosX", mEventAction->GetCaloMC_PosX()); +analysisManager->CreateNtupleDColumn("Calo_MC_PosY", mEventAction->GetCaloMC_PosY()); +analysisManager->CreateNtupleDColumn("Calo_MC_PosZ", mEventAction->GetCaloMC_PosZ()); +analysisManager->CreateNtupleDColumn("Calo_MC_DirX", mEventAction->GetCaloMC_DirX()); +analysisManager->CreateNtupleDColumn("Calo_MC_DirY", mEventAction->GetCaloMC_DirY()); +analysisManager->CreateNtupleDColumn("Calo_MC_DirZ", mEventAction->GetCaloMC_DirZ()); +analysisManager->FinishNtuple(); +*/ + +analysisManager->CreateNtuple("Fast_Sim_Info", "MC info at calo face"); +analysisManager->CreateNtupleIColumn("Calo_MC_PDG"); +analysisManager->CreateNtupleDColumn("Calo_MC_Energy"); +analysisManager->CreateNtupleDColumn("Calo_MC_PosX"); +analysisManager->CreateNtupleDColumn("Calo_MC_PosY"); +analysisManager->CreateNtupleDColumn("Calo_MC_PosZ"); +analysisManager->CreateNtupleDColumn("Calo_MC_DirX"); +analysisManager->CreateNtupleDColumn("Calo_MC_DirY"); +analysisManager->CreateNtupleDColumn("Calo_MC_DirZ"); +analysisManager->FinishNtuple(); + + +analysisManager->CreateNtuple("run", "Run data"); +analysisManager->CreateNtupleDColumn("N"); +analysisManager->FinishNtuple(); + +analysisManager->OpenFile(); +} + +/// End-of-run callback +void DDMLRunAction::end(const G4Run* aRun){ + auto analysisManager = G4AnalysisManager::Instance(); + analysisManager->FillNtupleDColumn(2, 0, aRun->GetNumberOfEvent()); + analysisManager->AddNtupleRow(2); + analysisManager->Write(); + analysisManager->CloseFile(); +} + +} // namespace + +#include "DDG4/Factories.h" +DECLARE_GEANT4ACTION_NS(ddml, DDMLRunAction) \ No newline at end of file diff --git a/src/PionCloudsModel.cc b/src/PionCloudsModel.cc index 6b04e9e..0dbad53 100644 --- a/src/PionCloudsModel.cc +++ b/src/PionCloudsModel.cc @@ -124,8 +124,8 @@ namespace ddml { for (int i = 0; i < nPoints; i++) { ddml::SpacePoint sp( - reshaped[i][0]*1e3, // x // *(-1) to align local to global convention in ddml - reshaped[i][2]*1e3, // y // *(-1) to align local to global convention in ddml + reshaped[i][0], // x // *(-1) to align local to global convention in ddml + reshaped[i][2], // y // *(-1) to align local to global convention in ddml 0., // z reshaped[i][3], // energy 0. // time From 40d2356d70163f9bda69433a30af360246fadca9 Mon Sep 17 00:00:00 2001 From: Peter McKeown Date: Tue, 25 Feb 2025 11:24:47 +0100 Subject: [PATCH 05/12] Updates ready for obtaining calo entries --- include/DDML/PionCloudsModel.h | 2 +- scripts/ddsim_steer_Record_Calo_entry.py | 5 +- scripts/test_onnx.mac | 17 ++- src/Geant4FastHitMakerGlobal.cc | 8 +- src/LoadHdf5.cc | 2 +- src/PionCloudsModel.cc | 129 ++++++++++++----------- 6 files changed, 93 insertions(+), 70 deletions(-) diff --git a/include/DDML/PionCloudsModel.h b/include/DDML/PionCloudsModel.h index 0bb9bfc..111a602 100644 --- a/include/DDML/PionCloudsModel.h +++ b/include/DDML/PionCloudsModel.h @@ -55,7 +55,7 @@ namespace ddml { /// model properties for plugin // These grid sizes were used for the two angle BIBAE - int m_numPoints = 2600; + int m_numPoints = 2600; //4148; //2600; //number of points in the shower 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) int m_nLayer = 78; diff --git a/scripts/ddsim_steer_Record_Calo_entry.py b/scripts/ddsim_steer_Record_Calo_entry.py index 6107367..1f21c95 100644 --- a/scripts/ddsim_steer_Record_Calo_entry.py +++ b/scripts/ddsim_steer_Record_Calo_entry.py @@ -542,7 +542,8 @@ def LoadHdf5(kernel): if hadrons == True: ml_model_had = "LoadHDF5PionCloudsPCHadronModelPolyhedraBarrel/BarrelModelTorch" - ml_had_file = "../models/PionClouds_50GeV_sp.h5" #"../models/PionClouds_50GeV_sp_scaled.h5" + ml_had_file = "../models/PionClouds_50GeV_sp_scaled.h5" #"../models/gen_showers_for_reco_50GeV_2000.hdf5" + #"../models/gen_showers_50GeV_2000_Martina.hdf5" #"../models/PionClouds_50GeV_sp_scaled.h5" #"../models/PionClouds_50GeV_sp.h5" # ml_correct_angles = False from g4units import GeV, MeV # DO NOT REMOVE OR MOVE!!!!! (EXCLAMATION MARK) @@ -631,7 +632,7 @@ def LoadHdf5(kernel): 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.Etrigger = {"pi+": 30.0 * GeV} #{"pi+": 10.0 * GeV} # trigger on lower training threshold modelHad1.FilePath = ml_had_file # model.OptimizeFlag = 1 # model.IntraOpNumThreads = 1 diff --git a/scripts/test_onnx.mac b/scripts/test_onnx.mac index 4226261..a748499 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 30 0 0 #30 180.48 0 #### this is ILD calo Face #120 180.48 120 #30 180.48 0 #30 0 0 + +# 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 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 -# for SimpleCalo -/gps/direction 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/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) diff --git a/src/Geant4FastHitMakerGlobal.cc b/src/Geant4FastHitMakerGlobal.cc index bd7108f..2395447 100644 --- a/src/Geant4FastHitMakerGlobal.cc +++ b/src/Geant4FastHitMakerGlobal.cc @@ -49,8 +49,8 @@ Geant4FastHitMakerGlobal::~Geant4FastHitMakerGlobal() { void Geant4FastHitMakerGlobal::make(const G4FastHit& aHit, const G4FastTrack& aTrack) { - std::cout << "Geant4FastHitMakerGlobal::make - aHit.GetPosition().x() " << aHit.GetPosition().x() << " aHit.GetPosition().y() " << - aHit.GetPosition().y() << " aHit.GetPosition().z() " << aHit.GetPosition().z() << std::endl; + //std::cout << "Geant4FastHitMakerGlobal::make - aHit.GetPosition().x() " << aHit.GetPosition().x() << " aHit.GetPosition().y() " << + // aHit.GetPosition().y() << " aHit.GetPosition().z() " << aHit.GetPosition().z() << std::endl; // do not make empty deposit if (aHit.GetEnergy() <= 0) { return; @@ -88,8 +88,8 @@ void Geant4FastHitMakerGlobal::make(const G4FastHit& aHit, const G4FastTrack& aT sensitive = currentVolume->GetLogicalVolume()->GetSensitiveDetector(); G4VFastSimSensitiveDetector* fastSimSensitive = dynamic_cast(sensitive); if (fastSimSensitive) { - std::cout << "IN FAST Sim Sensitive Geant4FastHitMakerGlobal::make - aHit.GetPosition().x() " << aHit.GetPosition().x() << " aHit.GetPosition().y() " << - aHit.GetPosition().y() << " aHit.GetPosition().z() " << aHit.GetPosition().z() << " Energy: " << aHit.GetEnergy() << std::endl; + // std::cout << "IN FAST Sim Sensitive Geant4FastHitMakerGlobal::make - aHit.GetPosition().x() " << aHit.GetPosition().x() << " aHit.GetPosition().y() " << + // aHit.GetPosition().y() << " aHit.GetPosition().z() " << aHit.GetPosition().z() << " Energy: " << aHit.GetEnergy() << std::endl; fastSimSensitive->Hit(&aHit, &aTrack, &m_touchableHandle); } else if (sensitive && currentVolume->GetLogicalVolume()->GetFastSimulationManager()) { G4cerr << "ERROR - Geant4FastHitMakerGlobal::make()" << G4endl << " It is required to derive from the " diff --git a/src/LoadHdf5.cc b/src/LoadHdf5.cc index 2cb8360..5007782 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 = "spase_points"; //"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(); diff --git a/src/PionCloudsModel.cc b/src/PionCloudsModel.cc index 0dbad53..888afc2 100644 --- a/src/PionCloudsModel.cc +++ b/src/PionCloudsModel.cc @@ -70,73 +70,84 @@ namespace ddml { } - - void PionCloudsModel::convertOutput(G4FastTrack const& aFastTrack, - G4ThreeVector const& localDir, - const std::vector& output, - std::vector& spacepoints ){ - - int nPoints = m_numPoints ; // number of points in shower - - dd4hep::printout(dd4hep::DEBUG, "PionCloudsModel::convertOutput", "m_numPoints : %i", m_numPoints); +/* For array structure: (No. showers, dimensions(4), No. points) +void PionCloudsModel::convertOutput(G4FastTrack const& aFastTrack, + G4ThreeVector const& localDir, + const std::vector& output, + std::vector& spacepoints ){ - spacepoints.resize( m_nLayer ) ; + //int nPoints = m_numPoints ; // number of points in shower - int numSP = 0; int layerNum = 0; - int numElements = 0; - for (int i = 0; i < m_nLayer; i++) { - numSP = output[i] + 1; - spacepoints[i].reserve(numSP); - numElements += output[i] * 4; - } - std::cout << " PionCloudsModel::convertOutput DONE numElements, numElements = " << numElements << std::endl; - - /* - for (int i = m_nLayer; i < m_nLayer + numElements; i += 4) { - std::cout << " PionCloudsModel::convertOutput Layer Loop, i = " << i << std::endl; - std::cout << " PionCloudsModel::convertOutput Layer Loop, output[i] " << output[i] << std::endl; - std::cout << " PionCloudsModel::convertOutput Layer Loop, output[i+1] " << output[i] << std::endl; - std::cout << " PionCloudsModel::convertOutput Layer Loop, output[i+3] " << output[i] << std::endl; - ddml::SpacePoint sp( - output[i], // x // *(-1) to align local to global convention in ddml - output[i+1], // y // *(-1) to align local to global convention in ddml - 0., // z - output[i+3], // energy - 0. // time - ); - layerNum = output[i+2]; - std::cout << " PionCloudsModel::convertOutput Layer Loop, layerNum " << layerNum << std::endl; - spacepoints[layerNum].emplace_back( sp ) ; - */ - - - float reshaped[2600][4]; - - // Fill the 3D array using the flattened vector - size_t index = 0; - for (size_t j = 0; j < 2600; ++j) { - for (size_t k = 0; k < 4; ++k) { + /// This is too C-like - once model is concrete use std::vector + //float reshaped[m_numPoints][4]; + std::vector> reshaped(4, std::vector(m_numPoints)); + + // Fill the 3D array-like vector using the flattened vector + int index = 0; + for (int j = 0; j < 4; ++j) { + for (int k = 0; k < m_numPoints; ++k) { reshaped[j][k] = output[index++]; } } - for (int i = 0; i < nPoints; i++) { - ddml::SpacePoint sp( - reshaped[i][0], // x // *(-1) to align local to global convention in ddml - reshaped[i][2], // y // *(-1) to align local to global convention in ddml - 0., // z - reshaped[i][3], // energy - 0. // time - ); - layerNum = reshaped[i][1]; - std::cout << "PionCloudsModel::convertOutput - layerNum" << layerNum <& output, + std::vector& spacepoints ){ + + //int nPoints = m_numPoints ; // number of points in shower + + int layerNum = 0; + + /// This is too C-like - use std::vector + //float reshaped[m_numPoints][4]; + std::vector> reshaped(m_numPoints, std::vector(4)); + + // Fill the 3D array-like vector using the flattened vector + int index = 0; + for (int j = 0; j < m_numPoints; ++j) { + for (int k = 0; k < 4; ++k) { + reshaped[j][k] = output[index]; + index++; + } + } + + spacepoints.resize( m_nLayer ) ; + + for (int i = 0; i < m_numPoints; i++) { + ddml::SpacePoint sp( + reshaped[i][0], // x // *(-1) to align local to global convention in ddml + reshaped[i][2], // y // *(-1) to align local to global convention in ddml + 0., // z + reshaped[i][3], // energy + 0. // time + ); + layerNum = reshaped[i][1]; + //std::cout << "PionCloudsModel::convertOutput - layerNum" << layerNum < Date: Tue, 25 Feb 2025 14:31:08 +0100 Subject: [PATCH 06/12] Update incorrect naming in ONNX inference --- src/ONNXInference.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ONNXInference.cc b/src/ONNXInference.cc index 09e9c4a..261d7ae 100644 --- a/src/ONNXInference.cc +++ b/src/ONNXInference.cc @@ -56,7 +56,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 @@ -87,7 +87,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 From a8372214592a046c67fd20632e36314e863c525f Mon Sep 17 00:00:00 2001 From: Peter McKeown Date: Tue, 25 Feb 2025 14:35:43 +0100 Subject: [PATCH 07/12] Correct capitalisation --- src/ONNXInference.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ONNXInference.cc b/src/ONNXInference.cc index 261d7ae..1f7b30e 100644 --- a/src/ONNXInference.cc +++ b/src/ONNXInference.cc @@ -56,7 +56,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(m_Session->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 @@ -87,7 +87,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(m_Session->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 From a2546234d31432e6e25598ead89b8708f3d84e62 Mon Sep 17 00:00:00 2001 From: Peter McKeown Date: Tue, 25 Feb 2025 14:52:24 +0100 Subject: [PATCH 08/12] Hot fix type of m_dimsOut for older HDF5 versions --- include/DDML/LoadHdf5.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/DDML/LoadHdf5.h b/include/DDML/LoadHdf5.h index 74bdf02..7fdfd44 100644 --- a/include/DDML/LoadHdf5.h +++ b/include/DDML/LoadHdf5.h @@ -49,7 +49,8 @@ class LoadHdf5 : public InferenceInterface { std::vector m_library{}; // shower library dimensions - std::vector m_dimsOut{}; + std::vector m_dimsOut{}; // This is a hot fix for older HDF5 versions! + //std::vector m_dimsOut{}; // properties for plugin std::string m_filePath = {}; From d4d497d7533ba6fc3a6be643140fb5193d58a8da Mon Sep 17 00:00:00 2001 From: Peter McKeown Date: Thu, 8 May 2025 11:19:58 +0200 Subject: [PATCH 09/12] code clean up --- scripts/ddsim_steer.py | 2 +- src/DDMLRunAction.cc | 14 ----------- src/PionCloudsModel.cc | 55 ------------------------------------------ 3 files changed, 1 insertion(+), 70 deletions(-) diff --git a/scripts/ddsim_steer.py b/scripts/ddsim_steer.py index 9c7f8a0..6d24a50 100644 --- a/scripts/ddsim_steer.py +++ b/scripts/ddsim_steer.py @@ -547,7 +547,7 @@ def LoadHdf5(kernel): if hadrons == True: ml_model_had = "LoadHDF5PionCloudsPCHadronModelPolyhedraBarrel/BarrelModelTorch" - ml_had_file = "../models/PionClouds_50GeV_sp.h5" #"../models/PionClouds_50GeV_sp_scaled.h5" + ml_had_file = "../models/PionClouds_50GeV_sp.h5" ml_correct_angles = False from g4units import GeV, MeV # DO NOT REMOVE OR MOVE!!!!! (EXCLAMATION MARK) diff --git a/src/DDMLRunAction.cc b/src/DDMLRunAction.cc index c7cd807..60c92eb 100644 --- a/src/DDMLRunAction.cc +++ b/src/DDMLRunAction.cc @@ -46,20 +46,6 @@ analysisManager->CreateNtupleDColumn("MC_DirY"); analysisManager->CreateNtupleDColumn("MC_DirZ"); analysisManager->FinishNtuple(); -/* -DDMLEventAction* mEventAction = dynamic_cast(G4EventManager::GetEventManager()->GetUserEventAction()); -analysisManager->CreateNtuple("Fast Sim Info", "MC info at calo face"); -analysisManager->CreateNtupleIColumn("Calo_MC_PDG", mEventAction->GetCaloMC_PDG()); -analysisManager->CreateNtupleDColumn("Calo_MC_Energy", mEventAction->GetCaloMC_E()); -analysisManager->CreateNtupleDColumn("Calo_MC_PosX", mEventAction->GetCaloMC_PosX()); -analysisManager->CreateNtupleDColumn("Calo_MC_PosY", mEventAction->GetCaloMC_PosY()); -analysisManager->CreateNtupleDColumn("Calo_MC_PosZ", mEventAction->GetCaloMC_PosZ()); -analysisManager->CreateNtupleDColumn("Calo_MC_DirX", mEventAction->GetCaloMC_DirX()); -analysisManager->CreateNtupleDColumn("Calo_MC_DirY", mEventAction->GetCaloMC_DirY()); -analysisManager->CreateNtupleDColumn("Calo_MC_DirZ", mEventAction->GetCaloMC_DirZ()); -analysisManager->FinishNtuple(); -*/ - analysisManager->CreateNtuple("Fast_Sim_Info", "MC info at calo face"); analysisManager->CreateNtupleIColumn("Calo_MC_PDG"); analysisManager->CreateNtupleDColumn("Calo_MC_Energy"); diff --git a/src/PionCloudsModel.cc b/src/PionCloudsModel.cc index 888afc2..7604f3a 100644 --- a/src/PionCloudsModel.cc +++ b/src/PionCloudsModel.cc @@ -2,8 +2,6 @@ #include // for G4FastTrack -//#include - #define DEBUGPRINT 0 namespace ddml { @@ -31,8 +29,6 @@ namespace ddml { dd4hep::printout(dd4hep::DEBUG, "PionCloudsModel::prepareInput", "DDML::localDir: (%f, %f, %f)", localDir_.x(), localDir_.y(), localDir_.z()); - // std::cout << "PionClouds::localDir:" << "(" << localDir_.x() << "," << localDir_.y() << "," << localDir_.z() << ")" - // << std::endl; // compute local incident angles double r = sqrt(localDir_.x() * localDir_.x() + localDir_.y() * localDir_.y() + localDir_.z() * localDir_.z()); @@ -54,12 +50,6 @@ namespace ddml { inputs[1][0] = theta; // 89.*(M_PI/180.) ; //Theta_vec[0]/(90.*(M_PI/180.)); inputs[2][0] = phi; - // if (DEBUGPRINT) { - // std::cout << " Input_energy_tensor : " << inputs[0][0] << std::endl; - // std::cout << " Input_theta_tensor : " << inputs[1][0] << std::endl; - // std::cout << " Input_phi_tensor : " << inputs[2][0] << std::endl; - // } - 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]); @@ -70,46 +60,6 @@ namespace ddml { } -/* For array structure: (No. showers, dimensions(4), No. points) -void PionCloudsModel::convertOutput(G4FastTrack const& aFastTrack, - G4ThreeVector const& localDir, - const std::vector& output, - std::vector& spacepoints ){ - - //int nPoints = m_numPoints ; // number of points in shower - - int layerNum = 0; - - /// This is too C-like - once model is concrete use std::vector - //float reshaped[m_numPoints][4]; - std::vector> reshaped(4, std::vector(m_numPoints)); - - // Fill the 3D array-like vector using the flattened vector - int index = 0; - for (int j = 0; j < 4; ++j) { - for (int k = 0; k < m_numPoints; ++k) { - reshaped[j][k] = output[index++]; - } - } - - spacepoints.resize( m_nLayer ) ; - - for (int i = 0; i < m_numPoints; i++) { - ddml::SpacePoint sp( - reshaped[0][i], // x // *(-1) to align local to global convention in ddml - reshaped[2][i], // y // *(-1) to align local to global convention in ddml - 0., // z - reshaped[3][i], // energy - 0. // time - ); - layerNum = reshaped[1][i]; - spacepoints[layerNum].emplace_back( sp ) ; - } - } -} -*/ - - // For array structure: (No. showers, No. points, dimensions(4)) void PionCloudsModel::convertOutput(G4FastTrack const& aFastTrack, G4ThreeVector const& localDir, @@ -120,8 +70,6 @@ void PionCloudsModel::convertOutput(G4FastTrack const& aFastTrack, int layerNum = 0; - /// This is too C-like - use std::vector - //float reshaped[m_numPoints][4]; std::vector> reshaped(m_numPoints, std::vector(4)); // Fill the 3D array-like vector using the flattened vector @@ -144,9 +92,6 @@ void PionCloudsModel::convertOutput(G4FastTrack const& aFastTrack, 0. // time ); layerNum = reshaped[i][1]; - //std::cout << "PionCloudsModel::convertOutput - layerNum" << layerNum < Date: Wed, 4 Jun 2025 11:00:46 +0200 Subject: [PATCH 10/12] Update include/DDML/DDMLEventAction.h Co-authored-by: Thomas Madlener --- include/DDML/DDMLEventAction.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/include/DDML/DDMLEventAction.h b/include/DDML/DDMLEventAction.h index 1aa857c..7f81c78 100644 --- a/include/DDML/DDMLEventAction.h +++ b/include/DDML/DDMLEventAction.h @@ -25,8 +25,6 @@ class DDMLEventAction: public dd4hep::sim::Geant4EventAction{ DDMLEventAction(dd4hep::sim::Geant4Context* c, const std::string& n); /// Default destructor virtual ~DDMLEventAction(); - /// begin-of-event callback - inline virtual void begin(const G4Event*) override; /// End-of-event callback virtual void end(const G4Event*) override; /// begin-of-run callback From 5a394563a11b8c168b9c56cf6188538532dee4ef Mon Sep 17 00:00:00 2001 From: peter-mckeown <74307317+peter-mckeown@users.noreply.github.com> Date: Wed, 4 Jun 2025 11:01:03 +0200 Subject: [PATCH 11/12] Update include/DDML/DDMLEventAction.h Co-authored-by: Thomas Madlener --- include/DDML/DDMLEventAction.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/DDML/DDMLEventAction.h b/include/DDML/DDMLEventAction.h index 7f81c78..ef1c6c5 100644 --- a/include/DDML/DDMLEventAction.h +++ b/include/DDML/DDMLEventAction.h @@ -26,7 +26,7 @@ class DDMLEventAction: public dd4hep::sim::Geant4EventAction{ /// Default destructor virtual ~DDMLEventAction(); /// End-of-event callback - virtual void end(const G4Event*) override; + void end(const G4Event*) override; /// begin-of-run callback void beginRun(const G4Run*); /// End-of-run callback From 455bd8d1d497e59c484ba70a891bbf8833540bac Mon Sep 17 00:00:00 2001 From: peter-mckeown <74307317+peter-mckeown@users.noreply.github.com> Date: Wed, 4 Jun 2025 11:01:24 +0200 Subject: [PATCH 12/12] Update include/DDML/DDMLEventAction.h Co-authored-by: Thomas Madlener --- include/DDML/DDMLEventAction.h | 4 ---- 1 file changed, 4 deletions(-) diff --git a/include/DDML/DDMLEventAction.h b/include/DDML/DDMLEventAction.h index ef1c6c5..290dc09 100644 --- a/include/DDML/DDMLEventAction.h +++ b/include/DDML/DDMLEventAction.h @@ -27,10 +27,6 @@ class DDMLEventAction: public dd4hep::sim::Geant4EventAction{ virtual ~DDMLEventAction(); /// End-of-event callback void end(const G4Event*) override; - /// begin-of-run callback - void beginRun(const G4Run*); - /// End-of-run callback - void endRun(const G4Run*); //// Get and Set methods for Calo Face info inline std::vector& GetCaloMC_PDG() {return m_CaloMCPDG;}