Hadrons calo entry merge

CMakeLists.txt

@@ -58,6 +58,7 @@ endif()
 
 option(INSTRUMENT_MODEL "Instrument the steps of the showerModel call" OFF)
 
+
 option(DOWNLOAD_MODELS "Download and install the models that are stored externally" ON)
 
 #---------------------------------------------------

include/DDML/DDMLEventAction.h

@@ -0,0 +1,62 @@
+#include "G4UserEventAction.hh" // for G4UserEventAction
+#include <G4Types.hh>           // for G4int, G4double
+#include <vector>               // for vector
+// #include "G4Timer.hh"            // for G4Timer
+class G4Event;
+#include "CLHEP/Units/PhysicalConstants.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "DDG4/Geant4EventAction.h"
+#include "DDG4/Geant4Handle.h"
+#include "DDG4/Geant4Kernel.h"
+
+namespace ddml {
+
+/** Event action for ddml.
+ *
+ *  @author P. McKeown, CERN
+ *  @date Feb 2025
+ *
+ */
+
+class FastSimTriggerInfoRecordingEventAction : public dd4hep::sim::Geant4EventAction {
+public:
+  /// Standard constructor with initializing arguments
+  FastSimTriggerInfoRecordingEventAction(dd4hep::sim::Geant4Context* c, const std::string& n);
+  /// Default destructor
+  virtual ~FastSimTriggerInfoRecordingEventAction();
+  /// 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 ddml

include/DDML/DDMLRunAction.h

@@ -0,0 +1,32 @@
+#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 FastSimTriggerInfoRecordingRunAction : public dd4hep::sim::Geant4RunAction {
+public:
+  FastSimTriggerInfoRecordingRunAction() = delete;
+  /// Standard constructor with initializing arguments
+  FastSimTriggerInfoRecordingRunAction(dd4hep::sim::Geant4Context* c, const std::string& n);
+  /// Default destructor
+  virtual ~FastSimTriggerInfoRecordingRunAction();
+  /// begin-of-run callback
+  void begin(const G4Run*) override;
+  /// End-of-run callback
+  void end(const G4Run*) override;
+};
+
+} // namespace ddml

include/DDML/FastMLShower.h

@@ -41,6 +41,9 @@ inline auto run_void_member_timed(Obj& obj, MemberFunc func, Args&&... args) {
 }
 #endif
 
+// Required to record calo entry position
+#include "G4AnalysisManager.hh"
+
 namespace ddml {
 /** The templated base class for running fast shower simulation with ML.
  *  The actual implementation is provided by the templated class
@@ -138,6 +141,28 @@ class FastMLShower : public dd4hep::sim::Geant4FastSimShowerModel {
     podio::UserDataCollection<uint64_t> nHits;
 #endif
 
+    if (m_fastsimML.m_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);
+    }
+
     for (auto& invec : m_input) {
       invec.clear();
     }
@@ -223,10 +248,19 @@ struct FastMLModel {
 
   FastMLModel() : hitMaker(std::make_unique<Geant4FastHitMakerGlobal>()) {}
 
+  const bool has_constructGeo = false;
+  const bool has_constructField = false;
+  const bool has_constructSensitives = false;
+  const bool has_check_applicability = false;
+  const bool has_check_trigger = false;
+
+  bool m_record_calo_impact = false;
+
   void declareProperties(dd4hep::sim::Geant4Action* plugin) {
     model.declareProperties(plugin);
     inference.declareProperties(plugin);
     geometry.declareProperties(plugin);
+    plugin->declareProperty("Record_Calo_Impact", this->m_record_calo_impact);
   }
 };
 

include/DDML/LoadHdf5.h

@@ -49,7 +49,7 @@ class LoadHdf5 {
   std::vector<float> m_library{};
 
   // shower library dimensions
-  std::vector<unsigned long> m_dimsOut{};
+  std::vector<hsize_t> m_dimsOut{};
 
   // dimensions of shower library
   int m_rank;

scripts/ddsim_steer.py

@@ -396,6 +396,23 @@ def aiDanceTorch(kernel):
     old_DD4hep = False  ## use for DD4hep versions/commits before ~ Apr 21st 2023
     hadrons = False
 
+    ## Record calo entry position?
+    Record_Calo_Entry = True
+
+    if Record_Calo_Entry:
+        from DDG4 import RunAction
+        from DDG4 import EventAction
+
+        ### Configure Run actions
+        run1 = RunAction(kernel, "FastSimTriggerInfoRecordingRunAction/runaction")
+        kernel.registerGlobalAction(run1)
+        kernel.runAction().add(run1)
+        event1 = EventAction(
+            kernel, "FastSimTriggerInfoRecordingEventAction/eventaction"
+        )
+        kernel.registerGlobalAction(event1)
+        kernel.eventAction().add(event1)
+
     if ild == True:
         ml_barrel_name = "EcalBarrel"
         ml_barrel_symmetry = 8
@@ -487,6 +504,10 @@ def aiDanceTorch(kernel):
     model.OptimizeFlag = 1
     model.IntraOpNumThreads = 1
 
+    ## Optional: record entry point of particles to calorimeter face
+    if Record_Calo_Entry:
+        model.Record_Calo_Impact = True
+
     model.enableUI()
     seq.adopt(model)
     # -------------------

scripts/test_onnx.mac

@@ -10,13 +10,13 @@
 
 #/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
+/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
+#/gps/particle  pi+
+#/gps/energy    50 GeV
 
 /gps/pos/centre -5 0 -15
 /gps/direction -0.037864591009775746 0.9992828792427412 0 # 90 deg imact for 50 GeV p+
@@ -75,5 +75,5 @@
 #/gps/direction  0.1 -0.8 .1
 
 
-/run/beamOn 2000 #100
+/run/beamOn 100 #2000 #100
 

src/CMakeLists.txt

@@ -12,6 +12,8 @@ set(sources
   EndcapTriggerTwoAngleBIBAE.cc
   CaloCloudsTwoAngleModel.cc
   L2LFlowsModel.cc
+  DDMLRunAction.cc
+  DDMLEventAction.cc
   )
 
 add_dd4hep_plugin(${PROJECT_NAME} SHARED ${sources})
@@ -97,6 +99,8 @@ set(headers
   ${PROJECT_SOURCE_DIR}/include/DDML/DDML.h
   ${PROJECT_SOURCE_DIR}/include/DDML/CaloCloudsTwoAngleModel.h
   ${PROJECT_SOURCE_DIR}/include/DDML/L2LFlowsModel.h
+  ${PROJECT_SOURCE_DIR}/include/DDML/DDMLRunAction.h
+  ${PROJECT_SOURCE_DIR}/include/DDML/DDMLEventAction.h
   )
 
 if(OnnxRuntime_FOUND)

src/DDMLEventAction.cc

@@ -0,0 +1,54 @@
+#include "DDML/DDMLEventAction.h"
+#include "DD4hep/InstanceCount.h"
+#include "DDG4/Geant4Data.h"
+#include "G4AnalysisManager.hh"
+#include "G4Event.hh"
+#include "G4EventManager.hh"
+#include "G4PrimaryParticle.hh"
+#include "G4PrimaryVertex.hh"
+#include "G4ThreeVector.hh"
+
+#define DEBUGPRINT 0
+
+namespace ddml {
+
+FastSimTriggerInfoRecordingEventAction::FastSimTriggerInfoRecordingEventAction(dd4hep::sim::Geant4Context* c,
+                                                                               const std::string& n)
+    : Geant4EventAction(c, n) {
+  dd4hep::InstanceCount::increment(this);
+}
+
+/// Default destructor
+FastSimTriggerInfoRecordingEventAction::~FastSimTriggerInfoRecordingEventAction() {
+  dd4hep::InstanceCount::decrement(this);
+}
+
+void FastSimTriggerInfoRecordingEventAction::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 ddml
+
+#include "DDG4/Factories.h"
+DECLARE_GEANT4ACTION_NS(ddml, FastSimTriggerInfoRecordingEventAction)

src/DDMLRunAction.cc

@@ -0,0 +1,80 @@
+#include "DDML/DDMLRunAction.h"
+#include "DD4hep/InstanceCount.h"
+#include "DDML/DDMLEventAction.h"
+#include "G4AnalysisManager.hh"
+#include "G4EventManager.hh"
+#include "G4Run.hh"
+
+#define DEBUGPRINT 0
+
+namespace ddml {
+
+/// Standard constructor with initializing arguments
+FastSimTriggerInfoRecordingRunAction::FastSimTriggerInfoRecordingRunAction(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
+FastSimTriggerInfoRecordingRunAction::~FastSimTriggerInfoRecordingRunAction() {
+  dd4hep::InstanceCount::decrement(this);
+}
+
+/// begin-of-run callback
+void FastSimTriggerInfoRecordingRunAction::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();
+
+  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 FastSimTriggerInfoRecordingRunAction::end(const G4Run* aRun) {
+  auto analysisManager = G4AnalysisManager::Instance();
+  analysisManager->FillNtupleDColumn(2, 0, aRun->GetNumberOfEvent());
+  analysisManager->AddNtupleRow(2);
+  analysisManager->Write();
+  analysisManager->CloseFile();
+}
+
+} // namespace ddml
+
+#include "DDG4/Factories.h"
+DECLARE_GEANT4ACTION_NS(ddml, FastSimTriggerInfoRecordingRunAction)

src/Geant4FastHitMakerGlobal.cc

@@ -48,7 +48,7 @@ Geant4FastHitMakerGlobal::~Geant4FastHitMakerGlobal() { delete m_navigator; }
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 void Geant4FastHitMakerGlobal::make(const G4FastHit& aHit, const G4FastTrack& aTrack) {
-  // do not make empty deposit
+  //  do not make empty deposit
   if (aHit.GetEnergy() <= 0) {
     return;
   }