make scene able to run parallel simulations

src/HEP/Geant/ActionInitialization.cpp

@@ -6,10 +6,12 @@ namespace uLib {
 namespace Geant {
 
 ActionInitialization::ActionInitialization(EmitterPrimary *emitter,
-                                           Vector<GeantEvent> *output)
+                                           Vector<GeantEvent> *output,
+                                           int verbosity)
     : G4VUserActionInitialization(),
       m_Emitter(emitter),
-      m_Output(output)
+      m_Output(output),
+      m_Verbosity(verbosity)
 {}
 
 ActionInitialization::~ActionInitialization() {}
@@ -21,7 +23,7 @@ void ActionInitialization::BuildForMaster() const {
 void ActionInitialization::Build() const {
   // Register the primary generator
   if (m_Emitter) {
-    SetUserAction(m_Emitter);
+    SetUserAction(m_Emitter->Clone());
   } else {
     // Fallback: default EmitterPrimary
     SetUserAction(new EmitterPrimary());
@@ -30,6 +32,7 @@ void ActionInitialization::Build() const {
   // Register actions
   if (m_Output) {
     SteppingAction *sa = new SteppingAction(m_Output);
+    sa->SetVerbosity(m_Verbosity);
     SetUserAction(static_cast<G4UserSteppingAction*>(sa));
     SetUserAction(static_cast<G4UserEventAction*>(sa));
   }

src/HEP/Geant/ActionInitialization.hh

@@ -15,9 +15,12 @@ public:
   /// @param emitter  the primary generator to use (owned by caller)
   /// @param output   pointer to the results vector (owned by caller)
   ActionInitialization(EmitterPrimary *emitter = nullptr,
-                       Vector<GeantEvent> *output = nullptr);
+                       Vector<GeantEvent> *output = nullptr,
+                       int verbosity = 0);
   ~ActionInitialization();
   
+  void SetVerbosity(int level) { m_Verbosity = level; }
+
   // Metodo chiamato solo dal thread principale (Master)
   virtual void BuildForMaster() const;
 
@@ -27,6 +30,7 @@ public:
 private:
   EmitterPrimary *m_Emitter;
   Vector<GeantEvent> *m_Output;
+  int m_Verbosity = 0;
 };
 
 } // namespace Geant

src/HEP/Geant/EmitterPrimary.cpp

@@ -61,6 +61,12 @@ void EmitterPrimary::GeneratePrimaries(G4Event *anEvent) {
   fParticleGun->GeneratePrimaryVertex(anEvent);
 }
 
+EmitterPrimary* EmitterPrimary::Clone() const {
+  auto* clone = new EmitterPrimary();
+  clone->SetMatrix(this->GetMatrix());
+  return clone;
+}
+
 // -------------------------------------------------------------------------- //
 // SkyPlaneEmitterPrimary using EcoMug
 
@@ -128,6 +134,13 @@ void SkyPlaneEmitterPrimary::GeneratePrimaries(G4Event* anEvent) {
     fParticleGun->GeneratePrimaryVertex(anEvent);
 }
 
+EmitterPrimary* SkyPlaneEmitterPrimary::Clone() const {
+    auto* clone = new SkyPlaneEmitterPrimary();
+    clone->SetSkySize(this->m_Size);
+    clone->SetMatrix(this->GetMatrix());
+    return clone;
+}
+ 
 // -------------------------------------------------------------------------- //
 // CylinderEmitterPrimary using EcoMug
 
@@ -188,6 +201,14 @@ void CylinderEmitterPrimary::GeneratePrimaries(G4Event* anEvent) {
     fParticleGun->GeneratePrimaryVertex(anEvent);
 }
 
+EmitterPrimary* CylinderEmitterPrimary::Clone() const {
+    auto* clone = new CylinderEmitterPrimary();
+    clone->SetRadius(this->m_Radius);
+    clone->SetHeight(this->m_Height);
+    clone->SetMatrix(this->GetMatrix());
+    return clone;
+}
+
 // -------------------------------------------------------------------------- //
 
 
@@ -298,6 +319,12 @@ void QuadMeshEmitterPrimary::GeneratePrimaries(G4Event *anEvent) {
     fParticleGun->GeneratePrimaryVertex(anEvent);
   }
 }
+EmitterPrimary* QuadMeshEmitterPrimary::Clone() const {
+  auto* clone = new QuadMeshEmitterPrimary();
+  if (m_Mesh) clone->SetMesh(m_Mesh);
+  clone->SetMatrix(this->GetMatrix());
+  return clone;
+}
 
 } // namespace Geant
 } // namespace uLib

src/HEP/Geant/EmitterPrimary.hh

@@ -34,6 +34,9 @@ class EmitterPrimary : public G4VUserPrimaryGeneratorAction, public Object, publ
 
     virtual void Updated() override { ULIB_SIGNAL_EMIT(EmitterPrimary::Updated); }
     
+    /// Create a clone of this emitter for multi-threading
+    virtual EmitterPrimary* Clone() const;
+
   protected:
     G4ParticleGun* fParticleGun; // Puntatore al cannone di particelle
 };
@@ -50,6 +53,8 @@ class SkyPlaneEmitterPrimary : public EmitterPrimary
     void SetSkySize(const uLib::Vector2f& size);
     uLib::Vector2f GetSkySize() const { return m_Size; }
     
+    virtual EmitterPrimary* Clone() const override;
+
   private:
     EcoMug* m_EcoMug;
     uLib::Vector2f m_Size;
@@ -68,6 +73,8 @@ class CylinderEmitterPrimary : public EmitterPrimary
     void SetHeight(float h);
     float GetHeight() const { return m_Height; }
     
+    virtual EmitterPrimary* Clone() const override;
+
   private:
     EcoMug* m_EcoMug;
     float m_Radius;
@@ -90,6 +97,8 @@ class QuadMeshEmitterPrimary : public EmitterPrimary
 
     void SetMesh(uLib::QuadMesh* mesh);
     
+    virtual EmitterPrimary* Clone() const override;
+
   private:
     uLib::QuadMesh* m_Mesh;
     std::vector<double> m_CumulativeAreas;

src/HEP/Geant/Scene.cpp

@@ -70,9 +70,10 @@ static void CheckGeant4Environment() {
 class SceneImpl {
 public:
   // constructor //
-  SceneImpl() : m_RunManager(G4RunManagerFactory::CreateRunManager(G4RunManagerType::Serial)),
+  SceneImpl() : m_RunManager(G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default)),
                 m_Emitter(nullptr),
-                m_Output(nullptr) {
+                m_Output(nullptr),
+                m_Verbosity(0) {
     m_RunManager->SetUserInitialization(new PhysicsList);
   }
 
@@ -86,7 +87,9 @@ public:
     // Set mandatory initialization classes for Geant4
     m_RunManager->SetUserInitialization(new SceneDetectorConstruction(this));
     m_RunManager->SetUserInitialization(
-        new ActionInitialization(m_Emitter, m_Output));
+        new ActionInitialization(m_Emitter, m_Output, m_Verbosity));
+    
+    m_RunManager->SetVerboseLevel(m_Verbosity);
 
     // Initialize Geant4
     m_RunManager->Initialize();
@@ -99,6 +102,7 @@ public:
   G4RunManager *m_RunManager;
   EmitterPrimary *m_Emitter;
   Vector<GeantEvent> *m_Output;
+  int m_Verbosity;
 };
 
 SceneDetectorConstruction::SceneDetectorConstruction(SceneImpl *owner) 
@@ -177,6 +181,10 @@ void Scene::SetEmitter(EmitterPrimary *emitter) {
 void Scene::Initialize() {
   d->Initialize();
 }
+void Scene::SetVerbosity(int level) {
+  d->m_Verbosity = level;
+  if (d->m_RunManager) d->m_RunManager->SetVerboseLevel(level);
+}
 
 void Scene::RunSimulation(int nEvents, Vector<GeantEvent> &results) {
   d->m_Output = &results;
@@ -185,7 +193,10 @@ void Scene::RunSimulation(int nEvents, Vector<GeantEvent> &results) {
   // (ActionInitialization was already set during Initialize, but we need
   //  to ensure the output pointer is current)
   d->m_RunManager->SetUserInitialization(
-      new ActionInitialization(d->m_Emitter, &results));
+      new ActionInitialization(d->m_Emitter, &results, d->m_Verbosity));
+
+  // Re-run initialization to propagate the ActionInitialization change to worker threads
+  d->m_RunManager->Initialize();
 
   d->m_RunManager->BeamOn(nEvents);
 }

src/HEP/Geant/Scene.h

@@ -61,6 +61,9 @@ public:
   /// Initialize the Geant4 run manager with detector, physics, and action.
   void Initialize();
   
+  /// Set the verbosity level for console output (default 0)
+  void SetVerbosity(int level);
+
   /// Run the simulation for nEvents muons.
   /// Results are appended to the provided vector.
   void RunSimulation(int nEvents, Vector<GeantEvent> &results);

src/HEP/Geant/SteppingAction.cpp

@@ -9,6 +9,9 @@
 #include "G4ParticleDefinition.hh"
 #include <set>
 #include <iostream>
+#include <mutex>
+
+static std::mutex g_SimulationOutputMutex;
 
 namespace uLib {
 namespace Geant {
@@ -46,21 +49,26 @@ void SteppingAction::BeginOfEventAction(const G4Event *event) {
 
 void SteppingAction::EndOfEventAction(const G4Event *event) {
   if (m_Output && !m_Current.m_Path.empty()) {
-    std::cout << "[Geant] Finished Event " << m_Current.m_EventID 
+    if (m_Verbosity > 0) {
-              << " with " << m_Current.m_Path.size() << " steps." << std::endl;
+      std::cout << "[Geant] Finished Event " << m_Current.m_EventID 
+                << " with " << m_Current.m_Path.size() << " steps." << std::endl;
       
-    // Check if we hit anything other than World
+      // Check if we hit anything other than World
-    std::set<std::string> volumes;
+      std::set<std::string> volumes;
-    for (const auto& delta : m_Current.m_Path) {
+      for (const auto& delta : m_Current.m_Path) {
-      if (!delta.m_SolidName.empty()) volumes.insert(delta.m_SolidName);
+        if (!delta.m_SolidName.empty()) volumes.insert(delta.m_SolidName);
-    }
+      }
-    if (volumes.size() > 1) {
+      if (volumes.size() > 1) {
-      std::cout << "        - Hit volumes: ";
+        std::cout << "        - Hit volumes: ";
-      for (const auto& v : volumes) std::cout << v << " ";
+        for (const auto& v : volumes) std::cout << v << " ";
-      std::cout << std::endl;
+        std::cout << std::endl;
+      }
     }
 
-    m_Output->push_back(m_Current);
+    {
+      std::lock_guard<std::mutex> lock(g_SimulationOutputMutex);
+      m_Output->push_back(m_Current);
+    }
   }
 }
 

src/HEP/Geant/SteppingAction.hh

@@ -21,10 +21,13 @@ public:
   virtual void BeginOfEventAction(const G4Event *event) override;
   virtual void EndOfEventAction(const G4Event *event) override;
   
+  void SetVerbosity(int level) { m_Verbosity = level; }
+
 private:
   Vector<GeantEvent> *m_Output;   ///< destination for finished events
   GeantEvent          m_Current;  ///< event being built
   int                 m_LastEventID; ///< track event transitions
+  int                 m_Verbosity = 0;
 };
 
 } // namespace Geant

src/HEP/Geant/testing/CMakeLists.txt

@@ -4,6 +4,7 @@ set(TESTS
        EventTest
        GeantApp
        ActionInitialization
+       SkyPlaneEmitterTest
 )
 
 set(LIBRARIES

src/HEP/Geant/testing/SkyPlaneEmitterTest.cpp

@@ -0,0 +1,84 @@
+#include "Geant/Solid.h"
+#include "HEP/Geant/GeantEvent.h"
+#include "HEP/Geant/Scene.h"
+#include "HEP/Geant/EmitterPrimary.hh"
+#include "Math/ContainerBox.h"
+#include "Math/Dense.h"
+#include "Math/Units.h"
+#include "HEP/Detectors/DetectorChamber.h"
+
+#include <Geant4/G4SystemOfUnits.hh>
+
+#include <iostream>
+
+using namespace uLib;
+
+int main(int argc, char** argv) {
+
+    int nEvents = 10000;
+    if (argc > 1) {
+        nEvents = std::stoi(argv[1]);
+    }
+
+    // 1. Setup Geant4 Scene
+    Geant::Scene scene;
+    scene.ConstructWorldBox(Vector3f(30_m, 30_m, 30_m), "G4_AIR");
+
+    ContainerBox iron_box;
+    iron_box.Scale(Vector3f(18_m, 18_m, 18_m));
+    iron_box.SetPosition(Vector3f(-9_m, -9_m, -9_m));
+    Geant::BoxSolid* iron_cube = new Geant::BoxSolid("IronCube", &iron_box);
+    iron_cube->SetNistMaterial("G4_Fe");
+    iron_cube->Update();
+    scene.AddSolid(iron_cube);
+
+    // Top Detector Chamber (along Y axis)
+    DetectorChamber* top_chamber_box = new DetectorChamber();
+    top_chamber_box->Scale(Vector3f(20_m, 40_cm, 20_m));
+    top_chamber_box->Rotate(90_deg, Vector3f(1, 0, 0));
+    top_chamber_box->SetPosition(Vector3f(-10_m, 12_m, -10_m));
+    Geant::BoxSolid* top_chamber = new Geant::BoxSolid("TopChamber", top_chamber_box);
+    top_chamber->SetNistMaterial("G4_AIR");
+    top_chamber->Update();
+    scene.AddSolid(top_chamber);
+
+    // Bottom Detector Chamber (along Y axis)
+    DetectorChamber* bottom_chamber_box = new DetectorChamber();
+    bottom_chamber_box->Scale(Vector3f(20_m, 40_cm, 20_m));
+    bottom_chamber_box->Rotate(90_deg, Vector3f(1, 0, 0));
+    bottom_chamber_box->SetPosition(Vector3f(-10_m, -12_m, -10_m));
+    Geant::BoxSolid* bottom_chamber = new Geant::BoxSolid("BottomChamber", bottom_chamber_box);
+    bottom_chamber->SetNistMaterial("G4_AIR");
+    bottom_chamber->Update();
+    scene.AddSolid(bottom_chamber);
+
+    // Setup SkyPlaneEmitterPrimary
+    Geant::SkyPlaneEmitterPrimary* emitter = new Geant::SkyPlaneEmitterPrimary();
+    emitter->SetPosition(Vector3f(0, 14.9_m, 0));
+    emitter->Rotate(-90_deg, Vector3f(1, 0, 0));
+    emitter->SetSkySize(Vector2f(20_m, 20_m)); 
+    
+    scene.SetEmitter(emitter);
+    scene.SetVerbosity(1);
+    scene.Initialize();
+
+    std::cout << "Starting simulation of " << nEvents << " events..." << std::endl;
+    Vector<Geant::GeantEvent> results;
+    scene.RunSimulation(nEvents, results);
+
+    std::cout << "Simulation finished. Collected " << results.size() << " events." << std::endl;
+
+    // Sample output to verify data collection
+    if (!results.empty()) {
+        std::cout << "Summary: " << std::endl;
+        std::cout << "  Total events generated: " << results.size() << std::endl;
+        size_t total_steps = 0;
+        for (const auto& event : results) {
+            total_steps += event.Path().size();
+        }
+        std::cout << "  Total simulation steps: " << total_steps << std::endl;
+        std::cout << "  Average steps per event: " << static_cast<double>(total_steps) / results.size() << std::endl;
+    }
+
+    return 0;
+}

src/Vtk/HEP/Geant/testing/vtkEmitterPrimaryTest.cpp

@@ -9,6 +9,8 @@
 #include "Vtk/HEP/Geant/vtkGeantEvent.h"
 #include "Vtk/HEP/Geant/vtkEmitterPrimary.h"
 #include "Vtk/vtkContainerBox.h"
+#include "HEP/Detectors/DetectorChamber.h"
+#include "Vtk/HEP/Detectors/vtkDetectorChamber.h"
 
 #include <vtkSmartPointer.h>
 #include <vtkCallbackCommand.h>
@@ -69,12 +71,30 @@ int main(int argc, char** argv) {
 
     ContainerBox iron_box;
     iron_box.Scale(Vector3f(10_m, 10_m, 10_m));
-    iron_box.SetPosition(Vector3f(0, 0, 0));
+    iron_box.SetPosition(Vector3f(-5_m, -5_m, -5_m));
     Geant::BoxSolid* iron_cube = new Geant::BoxSolid("IronCube", &iron_box);
     iron_cube->SetNistMaterial("G4_Fe");
     iron_cube->Update();
     scene.AddSolid(iron_cube);
 
+    // Top Detector Chamber
+    DetectorChamber* top_chamber_box = new DetectorChamber();
+    top_chamber_box->Scale(Vector3f(10_m, 10_m, 40_cm));
+    top_chamber_box->SetPosition(Vector3f(-5_m, -5_m, 7_m));
+    Geant::BoxSolid* top_chamber = new Geant::BoxSolid("TopChamber", top_chamber_box);
+    top_chamber->SetNistMaterial("G4_AIR");
+    top_chamber->Update();
+    scene.AddSolid(top_chamber);
+
+    // Bottom Detector Chamber
+    DetectorChamber* bottom_chamber_box = new DetectorChamber();
+    bottom_chamber_box->Scale(Vector3f(10_m, 10_m, 40_cm));
+    bottom_chamber_box->SetPosition(Vector3f(-5_m, -5_m, -7.1_m));
+    Geant::BoxSolid* bottom_chamber = new Geant::BoxSolid("BottomChamber", bottom_chamber_box);
+    bottom_chamber->SetNistMaterial("G4_AIR");
+    bottom_chamber->Update();
+    scene.AddSolid(bottom_chamber);
+
     Geant::EmitterPrimary* emitter = new Geant::EmitterPrimary();
     emitter->SetPosition(Vector3f(0, 0, 14_m));
     scene.SetEmitter(emitter);
@@ -97,6 +117,20 @@ int main(int argc, char** argv) {
     vtkIron->SetRepresentation(Vtk::Puppet::Surface);
     viewer.AddPuppet(*vtkIron);
 
+    // Visualize Top Chamber
+    Vtk::vtkDetectorChamber* vtkTop = new Vtk::vtkDetectorChamber(top_chamber_box);
+    vtkTop->SetOpacity(0.5);
+    vtkTop->SetColor(0.2, 0.8, 0.2);
+    vtkTop->SetRepresentation(Vtk::Puppet::Surface);
+    viewer.AddPuppet(*vtkTop);
+
+    // Visualize Bottom Chamber
+    Vtk::vtkDetectorChamber* vtkBottom = new Vtk::vtkDetectorChamber(bottom_chamber_box);
+    vtkBottom->SetOpacity(0.5);
+    vtkBottom->SetColor(0.2, 0.8, 0.2);
+    vtkBottom->SetRepresentation(Vtk::Puppet::Surface);
+    viewer.AddPuppet(*vtkBottom);
+
     // Visualize Emitter
     Vtk::vtkEmitterPrimary* vtkEmitter = new Vtk::vtkEmitterPrimary(*emitter);
     viewer.AddPuppet(*vtkEmitter);

src/Vtk/HEP/Geant/testing/vtkSkyPlaneEmitterPrimaryTest.cpp

@@ -9,6 +9,8 @@
 #include "Vtk/HEP/Geant/vtkGeantEvent.h"
 #include "Vtk/HEP/Geant/vtkEmitterPrimary.h"
 #include "Vtk/vtkContainerBox.h"
+#include "HEP/Detectors/DetectorChamber.h"
+#include "Vtk/HEP/Detectors/vtkDetectorChamber.h"
 
 #include <vtkSmartPointer.h>
 #include <vtkCallbackCommand.h>
@@ -65,13 +67,33 @@ int main(int argc, char** argv) {
     scene.ConstructWorldBox(Vector3f(30_m, 30_m, 30_m), "G4_AIR");
 
     ContainerBox iron_box;
-    iron_box.Scale(Vector3f(10_m, 10_m, 10_m));
+    iron_box.Scale(Vector3f(18_m, 18_m, 18_m));
-    iron_box.SetPosition(Vector3f(0, 0, 0));
+    iron_box.SetPosition(Vector3f(-9_m, -9_m, -9_m));
     Geant::BoxSolid* iron_cube = new Geant::BoxSolid("IronCube", &iron_box);
     iron_cube->SetNistMaterial("G4_Fe");
     iron_cube->Update();
     scene.AddSolid(iron_cube);
 
+    // Top Detector Chamber (along Y axis)
+    DetectorChamber* top_chamber_box = new DetectorChamber();
+    top_chamber_box->Scale(Vector3f(20_m, 40_cm, 20_m));
+    top_chamber_box->Rotate(90_deg, Vector3f(1, 0, 0));
+    top_chamber_box->SetPosition(Vector3f(-10_m, 12_m, -10_m));
+    Geant::BoxSolid* top_chamber = new Geant::BoxSolid("TopChamber", top_chamber_box);
+    top_chamber->SetNistMaterial("G4_AIR");
+    top_chamber->Update();
+    scene.AddSolid(top_chamber);
+
+    // Bottom Detector Chamber (along Y axis)
+    DetectorChamber* bottom_chamber_box = new DetectorChamber();
+    bottom_chamber_box->Scale(Vector3f(20_m, 40_cm, 20_m));
+    bottom_chamber_box->Rotate(90_deg, Vector3f(1, 0, 0));
+    bottom_chamber_box->SetPosition(Vector3f(-10_m, -12_m, -10_m));
+    Geant::BoxSolid* bottom_chamber = new Geant::BoxSolid("BottomChamber", bottom_chamber_box);
+    bottom_chamber->SetNistMaterial("G4_AIR");
+    bottom_chamber->Update();
+    scene.AddSolid(bottom_chamber);
+
     // Use SkyPlaneEmitterPrimary instead of EmitterPrimary
     Geant::SkyPlaneEmitterPrimary* emitter = new Geant::SkyPlaneEmitterPrimary();
     emitter->SetPosition(Vector3f(0, 14.9_m, 0));
@@ -96,6 +118,20 @@ int main(int argc, char** argv) {
     vtkIron->SetRepresentation(Vtk::Puppet::Surface);
     viewer.AddPuppet(*vtkIron);
 
+    // Visualize Top Chamber
+    Vtk::vtkDetectorChamber* vtkTop = new Vtk::vtkDetectorChamber(top_chamber_box);
+    vtkTop->SetOpacity(0.5);
+    vtkTop->SetColor(0.2, 0.8, 0.2);
+    vtkTop->SetRepresentation(Vtk::Puppet::Surface);
+    viewer.AddPuppet(*vtkTop);
+
+    // Visualize Bottom Chamber
+    Vtk::vtkDetectorChamber* vtkBottom = new Vtk::vtkDetectorChamber(bottom_chamber_box);
+    vtkBottom->SetOpacity(0.5);
+    vtkBottom->SetColor(0.2, 0.8, 0.2);
+    vtkBottom->SetRepresentation(Vtk::Puppet::Surface);
+    viewer.AddPuppet(*vtkBottom);
+
     // Use vtkSkyPlaneEmitterPrimary instead of vtkEmitterPrimary
     Vtk::vtkSkyPlaneEmitterPrimary* vtkEmitter = new Vtk::vtkSkyPlaneEmitterPrimary(*emitter);
     vtkEmitter->SetSelectable(false);