Hadrons calo entry merge

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

include/DDML/DDMLEventAction.h

@@ -0,0 +1,64 @@
+#include <G4Types.hh>            // for G4int, G4double
+#include <vector>                // 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();
+    /// End-of-event callback
+    void end(const G4Event*) override;
+
+    //// Get and Set methods for Calo Face info
+    inline std::vector<G4int>& GetCaloMC_PDG() {return m_CaloMCPDG;}
+    inline std::vector<G4double>& GetCaloMC_E() {return m_CaloMCE;}
+    inline std::vector<G4double>& GetCaloMC_PosX() {return m_CaloMCPosX;}
+    inline std::vector<G4double>& GetCaloMC_PosY() {return m_CaloMCPosY;}
+    inline std::vector<G4double>& GetCaloMC_PosZ() {return m_CaloMCPosZ;}
+    inline std::vector<G4double>& GetCaloMC_DirX() {return m_CaloMCDirX;}
+    inline std::vector<G4double>& GetCaloMC_DirY() {return m_CaloMCDirY;}
+    inline std::vector<G4double>& 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<G4int> m_CaloMCPDG;
+    std::vector<G4double> m_CaloMCE;
+    std::vector<G4double> m_CaloMCPosX;
+    std::vector<G4double> m_CaloMCPosY;
+    std::vector<G4double> m_CaloMCPosZ;
+    std::vector<G4double> m_CaloMCDirX;
+    std::vector<G4double> m_CaloMCDirY;
+    std::vector<G4double> m_CaloMCDirZ;
+};
+
+} //namespace
+

include/DDML/DDMLRunAction.h

@@ -0,0 +1,39 @@
+#include <CLHEP/Units/SystemOfUnits.h>       // for GeV
+#include <G4String.hh>                       // for G4String
+#include <G4ThreeVector.hh>                  // for G4ThreeVector
+#include <G4Types.hh>                        // 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

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<uint64_t> 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<DDMLEventAction*>(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();
     }

include/DDML/LoadHdf5.h

@@ -49,7 +49,8 @@ class LoadHdf5 : public InferenceInterface {
   std::vector<float> m_library{};
 
   // shower library dimensions
-  std::vector<unsigned long> m_dimsOut{};
+  std::vector<unsigned long long> m_dimsOut{}; // This is a hot fix for older HDF5 versions!
+  //std::vector<unsigned long> m_dimsOut{};
 
   // properties for plugin
   std::string m_filePath = {};

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<float>& output ) ;
+
+
+    /** create a vector of spacepoints per layer interpreting the model output  
+     */
+    virtual void convertOutput(G4FastTrack const& aFastTrack,
+			       G4ThreeVector const& localDir,
+			       const std::vector<float>& output,
+			       std::vector<SpacePointVec>& spacepoints ) ;
+
+
+
+  private:
+
+    /// model properties for plugin
+    // These grid sizes were used for the two angle BIBAE
+    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;
+
+    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
+

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 {
@@ -28,7 +33,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);
   }
 
@@ -52,6 +60,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

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,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.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 +574,8 @@ def LoadHdf5(kernel):
         seq.adopt(sensitives)
 
     # -----------------
+    '''
+    ## EM in Barrel
     model = DetectorConstruction(kernel, str(ml_model))
 
     ##   # Mandatory model parameters
@@ -578,7 +597,9 @@ def LoadHdf5(kernel):
 
     model.enableUI()
     seq.adopt(model)
+    '''
     # -------------------
+    ## EM in Endcap
     model1 = DetectorConstruction(kernel, str(ml_model_1))
 
     ##   # Mandatory model parameters
@@ -599,12 +620,36 @@ 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)