vtkGeantEvent

app/gcompose/src/main.cpp

@@ -42,7 +42,7 @@ int main(int argc, char** argv) {
     
 
     Geant::Scene scene;
-    scene.ConstructWorldBox(&world_box, "G4_AIR");
+    scene.ConstructWorldBox(world_box.GetSize(), "G4_AIR");
     scene.Initialize();
 
     // 2. Initialize MainWindow (contains embedded VTK QViewport)

src/HEP/Geant/ActionInitialization.cpp

@@ -1,30 +1,38 @@
 #include "ActionInitialization.hh"
 #include "EmitterPrimary.hh"
+#include "SteppingAction.hh"
 
 namespace uLib {
 namespace Geant {
 
-ActionInitialization::ActionInitialization() : G4VUserActionInitialization() {}
+ActionInitialization::ActionInitialization(EmitterPrimary *emitter,
+                                           Vector<GeantEvent> *output)
+    : G4VUserActionInitialization(),
+      m_Emitter(emitter),
+      m_Output(output)
+{}
 
 ActionInitialization::~ActionInitialization() {}
 
 void ActionInitialization::BuildForMaster() const {
-  // Questo metodo viene usato in modalità Multi-Threading.
-  // Serve per le azioni che devono esistere solo nel thread Master
-  // (tipicamente solo per inizializzare file di output o il RunAction globale).
-
-  // Esempio: SetUserAction(new RunAction());
+  // Master thread: no per-event actions needed
 }
 
 void ActionInitialization::Build() const {
-  // Questo è il cuore dell'inizializzazione per i thread di lavoro.
-  // Qui passiamo il nostro generatore di muoni a Geant4.
-  SetUserAction(new EmitterPrimary());
+  // Register the primary generator
+  if (m_Emitter) {
+    SetUserAction(m_Emitter);
+  } else {
+    // Fallback: default EmitterPrimary
+    SetUserAction(new EmitterPrimary());
+  }
 
-  // In una simulazione completa, qui registreresti anche le altre classi:
-  // SetUserAction(new RunAction());
-  // SetUserAction(new EventAction());
-  // SetUserAction(new SteppingAction());
+  // Register actions
+  if (m_Output) {
+    SteppingAction *sa = new SteppingAction(m_Output);
+    SetUserAction(static_cast<G4UserSteppingAction*>(sa));
+    SetUserAction(static_cast<G4UserEventAction*>(sa));
+  }
 }
 
 } // namespace Geant

src/HEP/Geant/ActionInitialization.hh

@@ -2,13 +2,20 @@
 #define ActionInitialization_h
 
 #include "G4VUserActionInitialization.hh"
+#include "Core/Vector.h"
 
 namespace uLib {
 namespace Geant {
 
+class EmitterPrimary;
+class GeantEvent;
+
 class ActionInitialization : public G4VUserActionInitialization {
 public:
-  ActionInitialization();
+  /// @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);
   ~ActionInitialization();
 
   // Metodo chiamato solo dal thread principale (Master)
@@ -16,6 +23,10 @@ public:
 
   // Metodo chiamato dai thread di lavoro (Worker) o in modalità sequenziale
   virtual void Build() const;
+
+private:
+  EmitterPrimary *m_Emitter;
+  Vector<GeantEvent> *m_Output;
 };
 
 } // namespace Geant

src/HEP/Geant/CMakeLists.txt

@@ -19,6 +19,7 @@ set(HEADERS
     DetectorConstruction.hh
     PhysicsList.hh
     ActionInitialization.hh
+    SteppingAction.hh
 )
 
 set(SOURCES
@@ -28,6 +29,7 @@ set(SOURCES
     DetectorConstruction.cpp
     PhysicsList.cpp
     ActionInitialization.cpp
+    SteppingAction.cpp
 )
 
 set(libname ${PACKAGE_LIBPREFIX}Geant)

src/HEP/Geant/GeantEvent.cpp

@@ -0,0 +1,25 @@
+
+#include "HEP/Geant/GeantEvent.h"
+
+#include <cstddef>
+#include <iostream>
+
+using namespace uLib;
+
+void GeantEvent::Print(const size_t size) const {
+  std::cout << "Event " << m_EventID << ":" << std::endl;
+  std::cout << "  Momentum: " << m_Momentum << std::endl;
+  std::cout << "  GenVector: " << m_GenVector << std::endl;
+  std::cout << "  Path: " << std::endl;
+  
+  size_t limit = m_Path.size();
+  if(size > 0 && size < m_Path.size()) {
+    limit = size;
+  }
+  for (size_t i = 0; i < limit; ++i) {
+    std::cout << "    " << i << ": " << m_Path[i].m_Length << " " << m_Path[i].m_Momentum << " " << m_Path[i].m_Direction << " " << m_Path[i].m_SolidName << std::endl;
+  }
+  if (limit < m_Path.size()) {
+    std::cout << "    ... (" << m_Path.size() - limit << " more deltas)" << std::endl;
+  }
+}

src/HEP/Geant/GeantEvent.h

@@ -30,34 +30,63 @@
 #include "Core/Vector.h"
 #include "Math/Dense.h"
 
-#include "Detectors/ChamberHitEvent.h"
+#include <string>
 
 namespace uLib {
 namespace Geant {
 
-class GeantEventData {
+// Forward declaration for friend access
+class SteppingAction;
+
+///////////////////////////////////////////////////////////////////////////////
+/// GeantEvent — output of a Geant4 simulation run.
+///
+/// m_Momentum and m_GenVector are set from the EmitterPrimary at generation.
+/// During simulation, each scattering interaction deposits a Delta in m_Path,
+/// recording the change of momentum and direction at each step boundary.
+///////////////////////////////////////////////////////////////////////////////
+
+class GeantEvent {
+  
 public:
-  uLibGetMacro(EventID, Id_t) uLibGetMacro(Momentum, Scalarf)
-  uLibConstRefMacro(GenPos, Vector3f) uLibConstRefMacro(GenDir, Vector3f)
-  uLibConstRefMacro(ChEvents, Vector<ChamberHitEventData>)
+
+  /// A single interaction step along the muon path.
+  struct Delta {
+    Scalarf     m_Length;     ///< step length through the solid
+    Scalarf     m_Momentum;  ///< momentum magnitude at this step
+    HVector3f   m_Direction; ///< direction after scattering
+    std::string m_SolidName; ///< name of the solid where interaction occurred
+
+    uLibGetMacro(Length, Scalarf)
+    uLibGetMacro(Momentum, Scalarf)
+    uLibConstRefMacro(Direction, HVector3f)
+    uLibConstRefMacro(SolidName, std::string)
+
+    Delta() : m_Length(0), m_Momentum(0), m_Direction(HVector3f::Zero()) {}
+  };
+
+  // --- Read-only accessors (public) --- //
+  uLibGetMacro(EventID, Id_t) 
+  uLibGetMacro(Momentum, Scalarf)
+  uLibConstRefMacro(GenVector, HLine3f)
+  uLibConstRefMacro(Path, Vector<Delta>)
+
+  void Print(const size_t size = 10) const;
 
 private:
-  friend class GeantEvent;
   Id_t m_EventID;
   Scalarf m_Momentum;
-  Vector3f m_GenPos;
-  Vector3f m_GenDir;
-  Vector<ChamberHitEventData> m_ChEvents;
-};
+  HLine3f m_GenVector;
+  Vector<Delta> m_Path;
 
-class GeantEvent : public GeantEventData {
 public:
-  uLibSetMacro(EventID, Id_t) uLibSetMacro(Momentum, Scalarf)
-  uLibRefMacro(GenPos, Vector3f) uLibRefMacro(GenDir, Vector3f)
-  uLibRefMacro(ChEvents, Vector<ChamberHitEventData>)
+  GeantEvent() : m_EventID(0), m_Momentum(0), m_GenVector(HLine3f()), m_Path() {}
+
+  // SteppingAction can populate the private fields during simulation
+  friend class SteppingAction;
 };
 
 } // namespace Geant
 } // namespace uLib
 
-#endif // GEANTEVENT_H
+#endif // U_GEANTEVENT_H

src/HEP/Geant/Scene.cpp

@@ -54,10 +54,27 @@ private:
   class SceneImpl *m_Owner;
 };
 
+static void CheckGeant4Environment() {
+  static bool checked = false;
+  if (checked) return;
+  checked = true;
+  if (!std::getenv("G4ENSDFSTATEDATA")) {
+    std::cerr << "********************************************************" << std::endl;
+    std::cerr << " WARNING: Geant4 environment variables are not set!" << std::endl;
+    std::cerr << " Please activate the environment before running:" << std::endl;
+    std::cerr << "   micromamba activate mutom" << std::endl;
+    std::cerr << "********************************************************" << std::endl;
+  }
+}
+
 class SceneImpl {
 public:
   // constructor //
-  SceneImpl() : m_RunManager(G4RunManagerFactory::CreateRunManager()) {}
+  SceneImpl() : m_RunManager(G4RunManagerFactory::CreateRunManager(G4RunManagerType::Serial)),
+                m_Emitter(nullptr),
+                m_Output(nullptr) {
+    m_RunManager->SetUserInitialization(new PhysicsList);
+  }
 
   // destructor //
   ~SceneImpl() { 
@@ -68,8 +85,8 @@ public:
   void Initialize() {
     // Set mandatory initialization classes for Geant4
     m_RunManager->SetUserInitialization(new SceneDetectorConstruction(this));
-    m_RunManager->SetUserInitialization(new PhysicsList);
-    m_RunManager->SetUserInitialization(new ActionInitialization);
+    m_RunManager->SetUserInitialization(
+        new ActionInitialization(m_Emitter, m_Output));
 
     // Initialize Geant4
     m_RunManager->Initialize();
@@ -78,7 +95,10 @@ public:
   // members //
   Vector<Solid *> m_Solids;
   Solid *m_World = nullptr;
+  ContainerBox m_WorldBox;
   G4RunManager *m_RunManager;
+  EmitterPrimary *m_Emitter;
+  Vector<GeantEvent> *m_Output;
 };
 
 SceneDetectorConstruction::SceneDetectorConstruction(SceneImpl *owner) 
@@ -91,10 +111,10 @@ G4VPhysicalVolume *SceneDetectorConstruction::Construct() {
 
 
 
-
-
-
-Scene::Scene() : d(new SceneImpl()) {}
+Scene::Scene() {
+  CheckGeant4Environment();
+  d = new SceneImpl();
+}
 Scene::~Scene() { delete d; }
 
 void Scene::AddSolid(Solid *solid, Solid *parent) { 
@@ -106,16 +126,22 @@ void Scene::AddSolid(Solid *solid, Solid *parent) {
   }
 }
 
-void Scene::ConstructWorldBox(const ContainerBox *box, const char *material) {
+const Solid* Scene::GetWorld() const { return d->m_World; }
+
+ContainerBox* Scene::GetWorldBox() const { return &d->m_WorldBox; }
+
+void Scene::ConstructWorldBox(const Vector3f &size, const char *material) {
   // Get nist material manager
+
+  d->m_WorldBox.Scale(size);
+  d->m_WorldBox.SetPosition(-size/2.0f);
+
   if (!d->m_World) {
     d->m_World = new Solid("World");
     d->m_World->SetNistMaterial(material);
     AddSolid(d->m_World);
   }
 
-  Vector3f size = box->GetSize();
-
   G4Box *solidWorld = new G4Box("World", 
                                 0.5 * size(0), 
                                 0.5 * size(1),
@@ -138,14 +164,32 @@ void Scene::ConstructWorldBox(const ContainerBox *box, const char *material) {
     true);
 
   d->m_World->SetPhysical(physWorld);
-                                                    
-  Matrix4f transform = box->GetMatrix(); 
-  d->m_World->SetTransform(transform);
+
+  // no transforms are allowed for the world box
+  // Matrix4f transform = box->GetMatrix(); 
+  // d->m_World->SetTransform(transform);
+}
+
+void Scene::SetEmitter(EmitterPrimary *emitter) {
+  d->m_Emitter = emitter;
 }
 
 void Scene::Initialize() {
   d->Initialize();
 }
 
+void Scene::RunSimulation(int nEvents, Vector<GeantEvent> &results) {
+  d->m_Output = &results;
+
+  // Re-initialize ActionInitialization with the output buffer
+  // (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));
+
+  d->m_RunManager->BeamOn(nEvents);
+}
+
 } // namespace Geant
 } // namespace uLib
+

src/HEP/Geant/Scene.h

@@ -31,12 +31,15 @@
 #include "Core/Object.h"
 #include "Core/Vector.h"
 #include "Solid.h"
+#include "GeantEvent.h"
 
 class G4VPhysicalVolume;
 
 namespace uLib {
 namespace Geant {
 
+class EmitterPrimary;
+
 class Scene : public Object {
 public:
   Scene();
@@ -44,10 +47,24 @@ public:
 
   void AddSolid(Solid *solid, Solid *parent = nullptr);
   
-  void ConstructWorldBox(const ContainerBox *box, const char *material);
+  void ConstructWorldBox(const Vector3f &size, const char *material);
 
+  /// Get the world box
+  const Solid* GetWorld() const;
+  
+  ContainerBox* GetWorldBox() const;
+
+  /// Set the primary generator (emitter) for the simulation.
+  /// The Scene does NOT take ownership of the emitter.
+  void SetEmitter(EmitterPrimary *emitter);
+
+  /// Initialize the Geant4 run manager with detector, physics, and action.
   void Initialize();
 
+  /// Run the simulation for nEvents muons.
+  /// Results are appended to the provided vector.
+  void RunSimulation(int nEvents, Vector<GeantEvent> &results);
+
 private:
   class SceneImpl *d;
 };

src/HEP/Geant/Solid.cpp

@@ -55,23 +55,32 @@ public:
 };
 
 Solid::Solid()
-    : m_Name("unnamed_solid"), m_Material(NULL), m_Logical(NULL), m_Physical(NULL) {}
+    : m_Name("unnamed_solid"), m_Material(NULL), m_Logical(NULL), m_Physical(NULL), 
+      m_Position(new G4ThreeVector(0,0,0)), m_Rotation(NULL) {}
 
 Solid::Solid(const char *name)
-    : m_Name(name), m_Material(NULL), m_Logical(NULL), m_Physical(NULL) {}
+    : m_Name(name), m_Material(NULL), m_Logical(NULL), m_Physical(NULL), 
+      m_Position(new G4ThreeVector(0,0,0)), m_Rotation(NULL) {}
+
+Solid::~Solid() {
+  if (m_Position) delete m_Position;
+  if (m_Rotation) delete m_Rotation;
+}
 
 void Solid::SetNistMaterial(const char *name) {
   G4NistManager *nist = G4NistManager::Instance();
-  m_Material = nist->FindOrBuildMaterial(name);
-  if (m_Logical)
-    m_Logical->SetMaterial(m_Material);
+  G4Material *mat = nist->FindOrBuildMaterial(name);
+  if (mat) SetMaterial(mat);
 }
 
 void Solid::SetMaterial(G4Material *material) {
   if (material) {
     m_Material = material;
-    if (m_Logical)
+    if (m_Logical) {
       m_Logical->SetMaterial(material);
+    } else if (GetG4Solid()) {
+      m_Logical = new G4LogicalVolume(GetG4Solid(), m_Material, GetName());
+    }
   }
 }
 
@@ -79,21 +88,22 @@ void Solid::SetTransform(Matrix4f transform) {
   uLib::AffineTransform t;
   t.SetMatrix(transform);
 
-    // 2. Extracto position and rotation for Geant4
+    // 2. Extract position and rotation for Geant4
     Vector3f pos = t.GetPosition();
-    G4ThreeVector g4pos(pos(0), pos(1), pos(2));
+    if (!m_Position) m_Position = new G4ThreeVector();
+    *m_Position = G4ThreeVector(pos(0), pos(1), pos(2));
 
     // Create a G4 rotation matrix from the 4x4 matrix
     Matrix3f m = t.GetRotation();
-    G4RotationMatrix* rot = new G4RotationMatrix();
-    rot->set(G4ThreeVector(m(0,0), m(1,0), m(2,0)), 
-             G4ThreeVector(m(0,1), m(1,1), m(2,1)), 
-             G4ThreeVector(m(0,2), m(1,2), m(2,2))); 
+    if (!m_Rotation) m_Rotation = new G4RotationMatrix();
+    m_Rotation->set(G4ThreeVector(m(0,0), m(1,0), m(2,0)), 
+                    G4ThreeVector(m(0,1), m(1,1), m(2,1)), 
+                    G4ThreeVector(m(0,2), m(1,2), m(2,2))); 
 
-    // 3. Se l'oggetto è già stato piazzato, aggiorniamo la sua trasformazione
+    // 3. If object is already placed, update its transformation
     if (m_Physical) {
-        m_Physical->SetTranslation(g4pos);
-        m_Physical->SetRotation(rot);
+        m_Physical->SetTranslation(*m_Position);
+        m_Physical->SetRotation(m_Rotation);
     }
 }
 
@@ -119,8 +129,8 @@ void Solid::SetParent(Solid *parent) {
 
   // G4PVPlacement
   m_Physical = new G4PVPlacement(
-      nullptr,               // Rotation
-      G4ThreeVector(0,0,0),  // Position (translation) inside the parent
+      m_Rotation,            // Rotation
+      *m_Position,           // Position (translation) inside the parent
       m_Logical,             // The logical volume of this solid (the child)
       GetName(),             // Name of the physical volume
       parentLogical,         // The logical volume of the parent (nullptr if it's the World volume)
@@ -137,7 +147,8 @@ void Solid::SetParent(Solid *parent) {
 
 
 TessellatedSolid::TessellatedSolid(const char *name)
-    : BaseClass(name), m_Solid(new G4TessellatedSolid(name)) {}
+    : BaseClass(name), m_Solid(new G4TessellatedSolid(name)) {
+    }
 
 void TessellatedSolid::SetMesh(TriangleMesh &mesh) {
   G4TessellatedSolid *ts = this->m_Solid;
@@ -149,7 +160,9 @@ void TessellatedSolid::SetMesh(TriangleMesh &mesh) {
         DetectorsSolidImpl::getG4Vector3f(mesh.Points().at(trg(2))), ABSOLUTE);
     ts->AddFacet((G4VFacet *)facet);
   }
-  this->m_Logical->SetSolid(ts);
+  if (this->m_Logical) {
+    this->m_Logical->SetSolid(ts);
+  }
 }
 
 
@@ -157,10 +170,13 @@ void TessellatedSolid::SetMesh(TriangleMesh &mesh) {
 
 
 BoxSolid::BoxSolid(const char *name, ContainerBox *box) : BaseClass(name) {
-  m_Solid = new G4Box(name, 0.5, 0.5, 0.5);
+  m_Solid = new G4Box(name, 1,1,1);
   m_Object = box;
   Object::connect(box, &ContainerBox::Updated, this, &BoxSolid::Update);
-  this->m_Logical->SetSolid(m_Solid);
+  if (m_Logical) {
+    m_Logical->SetSolid(m_Solid);
+  }
+  Update();
 }
 
 void BoxSolid::Update() {
@@ -170,8 +186,20 @@ void BoxSolid::Update() {
     m_Solid->SetYHalfLength(size(1) * 0.5);
     m_Solid->SetZHalfLength(size(2) * 0.5);
     
-    this->SetTransform(m_Object->GetMatrix());
-    // this->m_Logical->SetSolid(m_Solid);
+    // Geant4 placement is relative to center. uLib Box is anchored at corner.
+    // 1. Get position and rotation (clean, without scale)
+    Vector3f pos = m_Object->GetPosition();
+    Matrix3f rot = m_Object->GetRotation();
+
+    // 2. Center = Corner + Rotation * (Half-Size)
+    // We must rotate the offset vector because uLib box can be rotated.
+    Vector3f center = pos + rot * (size * 0.5);
+
+    uLib::AffineTransform t;
+    t.SetPosition(center);
+    t.SetRotation(rot);
+    
+    this->SetTransform(t.GetMatrix());
   }
 }
 

src/HEP/Geant/Solid.h

@@ -45,6 +45,7 @@ class Solid : public Object {
 public:
   Solid();
   Solid(const char *name);
+  virtual ~Solid();
 
   void SetNistMaterial(const char *name);
   void SetMaterial(G4Material *material);
@@ -59,6 +60,8 @@ public:
   uLibGetSetMacro(Logical, G4LogicalVolume *)
   uLibGetSetMacro(Physical, G4VPhysicalVolume *)
 
+  virtual G4VSolid* GetG4Solid() const { return nullptr; }
+
   inline const char *GetName() const { 
     return m_Logical ? m_Logical->GetName().c_str() : m_Name.c_str(); 
   }
@@ -83,6 +86,7 @@ public:
   TessellatedSolid(const char *name);
   void SetMesh(TriangleMesh &mesh);
   uLibGetMacro(Solid, G4TessellatedSolid *) 
+  virtual G4VSolid* GetG4Solid() const override { return (G4VSolid*)m_Solid; }
 
 public slots:
   void Update();
@@ -101,6 +105,7 @@ class BoxSolid : public Solid {
 
 public:
   BoxSolid(const char *name, ContainerBox *box);
+  virtual G4VSolid* GetG4Solid() const override { return (G4VSolid*)m_Solid; }
 
 public slots:
   void Update();

src/HEP/Geant/SteppingAction.cpp

@@ -0,0 +1,104 @@
+#include "SteppingAction.hh"
+
+#include "G4Step.hh"
+#include "G4Track.hh"
+#include "G4Event.hh"
+#include "G4RunManager.hh"
+#include "G4LogicalVolume.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4ParticleDefinition.hh"
+#include <set>
+#include <iostream>
+
+namespace uLib {
+namespace Geant {
+
+SteppingAction::SteppingAction(Vector<GeantEvent> *output)
+    : G4UserSteppingAction(),
+      m_Output(output),
+      m_Current(),
+      m_LastEventID(-1)
+{}
+
+SteppingAction::~SteppingAction() {}
+
+void SteppingAction::BeginOfEventAction(const G4Event *event) {
+  if (!event || !m_Output) return;
+
+  // Start a new GeantEvent
+  m_Current = GeantEvent();
+  m_Current.m_EventID = static_cast<Id_t>(event->GetEventID());
+
+  // Set initial momentum and generation vector from primary vertex
+  if (event->GetNumberOfPrimaryVertex() > 0) {
+    G4PrimaryVertex *vtx = event->GetPrimaryVertex(0);
+    G4ThreeVector pos = vtx->GetPosition();
+    m_Current.m_GenVector.origin = HPoint3f(pos.x(), pos.y(), pos.z());
+
+    if (vtx->GetNumberOfParticle() > 0) {
+      G4PrimaryParticle *prim = vtx->GetPrimary(0);
+      G4ThreeVector mom = prim->GetMomentumDirection();
+      m_Current.m_GenVector.direction = HVector3f(mom.x(), mom.y(), mom.z());
+      m_Current.m_Momentum = static_cast<Scalarf>(prim->GetTotalMomentum() / MeV);
+    }
+  }
+}
+
+void SteppingAction::EndOfEventAction(const G4Event *event) {
+  if (m_Output && !m_Current.m_Path.empty()) {
+    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
+    std::set<std::string> volumes;
+    for (const auto& delta : m_Current.m_Path) {
+      if (!delta.m_SolidName.empty()) volumes.insert(delta.m_SolidName);
+    }
+    if (volumes.size() > 1) {
+      std::cout << "        - Hit volumes: ";
+      for (const auto& v : volumes) std::cout << v << " ";
+      std::cout << std::endl;
+    }
+
+    m_Output->push_back(m_Current);
+  }
+}
+
+void SteppingAction::UserSteppingAction(const G4Step *step) {
+  if (!step || !m_Output) return;
+
+  const G4Track *track = step->GetTrack();
+  if (!track) return;
+
+  // Only record primary particle (muon)
+  if (track->GetParentID() != 0) return;
+
+  // Record a Delta for this step
+  GeantEvent::Delta delta;
+
+  // Step length
+  delta.m_Length = static_cast<Scalarf>(step->GetStepLength() / mm);
+
+  // Post-step momentum
+  G4ThreeVector postMom = track->GetMomentum();
+  delta.m_Momentum = static_cast<Scalarf>(postMom.mag() / MeV);
+
+  // Post-step direction
+  G4ThreeVector dir = track->GetMomentumDirection();
+  delta.m_Direction = HVector3f(static_cast<float>(dir.x()),
+                                 static_cast<float>(dir.y()),
+                                 static_cast<float>(dir.z()));
+
+  // Solid name where the step occurred
+  const G4LogicalVolume *vol = track->GetVolume() 
+      ? track->GetVolume()->GetLogicalVolume() 
+      : nullptr;
+  if (vol) {
+    delta.m_SolidName = vol->GetName();
+  }
+
+  m_Current.m_Path.push_back(delta);
+}
+
+} // namespace Geant
+} // namespace uLib

src/HEP/Geant/SteppingAction.hh

@@ -0,0 +1,33 @@
+#ifndef U_GEANT_STEPPINGACTION_HH
+#define U_GEANT_STEPPINGACTION_HH
+
+#include "G4UserSteppingAction.hh"
+#include "G4UserEventAction.hh"
+#include "Core/Vector.h"
+#include "GeantEvent.h"
+
+namespace uLib {
+namespace Geant {
+
+/// SteppingAction collects scattering data at each Geant4 step and
+/// builds GeantEvent objects in the output buffer.
+class SteppingAction : public G4UserSteppingAction, public G4UserEventAction {
+public:
+  /// @param output  pointer to the results vector owned by the Scene
+  SteppingAction(Vector<GeantEvent> *output);
+  virtual ~SteppingAction();
+
+  virtual void UserSteppingAction(const G4Step *step) override;
+  virtual void BeginOfEventAction(const G4Event *event) override;
+  virtual void EndOfEventAction(const G4Event *event) override;
+
+private:
+  Vector<GeantEvent> *m_Output;   ///< destination for finished events
+  GeantEvent          m_Current;  ///< event being built
+  int                 m_LastEventID; ///< track event transitions
+};
+
+} // namespace Geant
+} // namespace uLib
+
+#endif

src/HEP/Geant/testing/EventTest.cpp

@@ -1,6 +1,11 @@
 #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/TriangleMesh.h"
+#include "Math/Dense.h"
+#include "Math/Units.h"
 #include "testing-prototype.h"
 #include <Geant4/G4Material.hh>
 #include <Geant4/G4NistManager.hh>
@@ -13,59 +18,72 @@ using namespace uLib;
 int main() {
     BEGIN_TESTING(Geant Event);
 
-    // Test Solid initialization and NIST material //
+    // Test: Scene with iron cube in air, launch muons, collect events //
     {
-        Geant::Solid solid("test_solid");
-        TEST1(solid.GetLogical() != nullptr);
-        
-        solid.SetNistMaterial("G4_AIR");
-        TEST1(solid.GetMaterial() != nullptr);
-        TEST1(solid.GetMaterial()->GetName() == "G4_AIR");
-    }
+        // 1. Create world box (air, 30m x 30m x 30m)
+        Geant::Scene scene;
+        scene.ConstructWorldBox(Vector3f(30_m, 30_m, 30_m), "G4_AIR");
 
-    // Test TessellatedSolid with a simple mesh //
-    {
-        Geant::TessellatedSolid tsolid("test_tessellated");
-        TEST1(tsolid.GetLogical() != nullptr);
-        TEST1(tsolid.GetSolid() != nullptr);
+        // 2. Create iron cube (1m x 1m x 1m) at center
+        ContainerBox iron_box(Vector3f(1000, 1000, 1000)); // mm
+        Geant::BoxSolid *iron_cube = new Geant::BoxSolid("IronCube", &iron_box);
+        iron_cube->SetNistMaterial("G4_Fe");
+        iron_cube->Update(); // apply dimensions
+        scene.AddSolid(iron_cube);
 
-        // cube mesh //
-        TriangleMesh mesh;
-        mesh.AddPoint(Vector3f(0,0,0));
-        mesh.AddPoint(Vector3f(1,0,0));
-        mesh.AddPoint(Vector3f(0,1,0));
-        mesh.AddPoint(Vector3f(1,1,0));
-        mesh.AddPoint(Vector3f(0,0,1));
-        mesh.AddPoint(Vector3f(1,0,1));
-        mesh.AddPoint(Vector3f(0,1,1));
-        mesh.AddPoint(Vector3f(1,1,1));
-        
-        // create triangles (consistent outward winding) //
-        // bottom (z=0)
-        mesh.AddTriangle(Vector3i(0,2,3));
-        mesh.AddTriangle(Vector3i(0,3,1));
-        // top (z=1)
-        mesh.AddTriangle(Vector3i(4,5,7));
-        mesh.AddTriangle(Vector3i(4,7,6));
-        // left (x=0)
-        mesh.AddTriangle(Vector3i(0,4,6));
-        mesh.AddTriangle(Vector3i(0,6,2));
-        // right (x=1)
-        mesh.AddTriangle(Vector3i(1,3,7));
-        mesh.AddTriangle(Vector3i(1,7,5));
-        // front (y=0)
-        mesh.AddTriangle(Vector3i(0,1,5));
-        mesh.AddTriangle(Vector3i(0,5,4));
-        // back (y=1)
-        mesh.AddTriangle(Vector3i(2,6,7));
-        mesh.AddTriangle(Vector3i(2,7,3));
+        // 3. Set up emitter (default: mu- at 1 GeV, from z=+10m downward)
+        Geant::EmitterPrimary *emitter = new Geant::EmitterPrimary();
+        scene.SetEmitter(emitter);
 
+        // 4. Initialize Geant4
+        scene.Initialize();
 
+        // 5. Run simulation: 10 muons
+        int nEvents = 10;
+        Vector<Geant::GeantEvent> results;
+        scene.RunSimulation(nEvents, results);
 
-        tsolid.SetMesh(mesh);
-        // GeantEvent geant_event;
-        
-        
+        // 6. Check results
+        printf("  Collected %zu events\n", results.size());
+        TEST1(results.size() > 0);
+
+        for (size_t i = 0; i < results.size(); ++i) {
+            const Geant::GeantEvent &ev = results[i];
+            bool hitIron = false;
+            for (const auto &d : ev.Path()) {
+                if (d.SolidName() == "IronCube") {
+                    hitIron = true;
+                    break;
+                }
+            }
+
+            printf("  Event %d: momentum=%.1f MeV, path steps=%zu, hitIron=%s\n",
+                   ev.GetEventID(),
+                   ev.GetMomentum(),
+                   ev.Path().size(),
+                   hitIron ? "YES" : "NO");
+            
+            // Each event should have at least one step
+            TEST1(ev.Path().size() > 0);
+
+            // Print first few deltas
+            const auto &path = ev.Path();
+            for (size_t j = 0; j < path.size() && j < 5; ++j) {
+                const auto &d = path[j];
+                printf("    Delta[%zu]: solid=%s len=%.2f mm, p=%.1f MeV, "
+                       "dir=(%.3f, %.3f, %.3f)\n",
+                       j,
+                       d.SolidName().c_str(),
+                       d.GetLength(),
+                       d.GetMomentum(),
+                       d.Direction()(0),
+                       d.Direction()(1),
+                       d.Direction()(2));
+            }
+            if (path.size() > 5) {
+                printf("    ... (%zu more deltas)\n", path.size() - 5);
+            }
+        }
     }
 
     END_TESTING

src/HEP/Geant/testing/GeantApp.cpp

@@ -8,10 +8,8 @@ using namespace uLib;
 
 int main() {
 
-  uLib::ContainerBox world_box(Vector3f(100, 100, 100));
   uLib::Geant::Scene scene;
-
-  scene.ConstructWorldBox(&world_box, "G4_AIR");
+  scene.ConstructWorldBox(Vector3f(100, 100, 100), "G4_AIR");
   scene.Initialize();
 
   return 0;

src/HEP/Geant/testing/SolidTest.cpp

@@ -15,9 +15,12 @@ int main() {
     // Test Solid initialization and NIST material //
     {
         Geant::Solid solid("test_solid");
-        TEST1(solid.GetLogical() != nullptr);
+        // Logical volume is not created until material and solid are set
+        TEST1(solid.GetLogical() == nullptr);
         
         solid.SetNistMaterial("G4_AIR");
+        // Still null because base Solid has no GetG4Solid()
+        TEST1(solid.GetLogical() == nullptr);
         TEST1(solid.GetMaterial() != nullptr);
         TEST1(solid.GetMaterial()->GetName() == "G4_AIR");
     }
@@ -25,6 +28,7 @@ int main() {
     // Test TessellatedSolid with a simple mesh //
     {
         Geant::TessellatedSolid tsolid("test_tessellated");
+        tsolid.SetNistMaterial("G4_AIR");
         TEST1(tsolid.GetLogical() != nullptr);
         TEST1(tsolid.GetSolid() != nullptr);
 

src/Math/ContainerBox.h

@@ -101,7 +101,11 @@ public:
    * @brief Gets the current size (scale) of the box.
    * @return The size vector.
    */
-  inline Vector3f GetSize() const { return m_LocalT.GetScale(); }
+  inline Vector3f GetSize() const { 
+    Vector3f s = this->GetScale();
+    Vector3f ls = m_LocalT.GetScale();
+    return Vector3f(s(0) * ls(0), s(1) * ls(1), s(2) * ls(2));
+  }
 
   /**
    * @brief Swaps two local axes of the box.

src/Math/Transform.h

@@ -105,7 +105,11 @@ public:
 
     inline void Scale(const Vector3f v) { this->m_T.scale(v); }
 
-    inline Vector3f GetScale() const { return this->m_T.linear() * Vector3f(1,1,1); } // FIXXXXXXX
+    inline Vector3f GetScale() const { 
+        return Vector3f(m_T.linear().col(0).norm(), 
+                        m_T.linear().col(1).norm(), 
+                        m_T.linear().col(2).norm()); 
+    }
 
 
     inline void Rotate(const Matrix3f m) { this->m_T.rotate(m); }

src/Math/testing/ContainerBoxTest.cpp

@@ -122,7 +122,7 @@ int main()
 
   {
     ContainerBox Box;
-    Box.SetPosition(Vector3f(1,1,1));
+    Box.SetOrigin(Vector3f(1,1,1));
     Box.SetSize(Vector3f(2,2,2));
     Box.EulerYZYRotate(Vector3f(0,0,0));
     HPoint3f pt = Box.GetLocalPoint(HPoint3f(2,3,2));
@@ -130,6 +130,39 @@ int main()
     TEST0( Vector4f0(wp - HPoint3f(2,3,2)) );
   }
 
+  {
+    ContainerBox Box;
+    Box.SetPosition(Vector3f(-0.5,-0.5,-0.5));
+    Box.SetSize(Vector3f(1,1,1));
+    HPoint3f pt = Box.GetLocalPoint(HPoint3f(0,0,0));
+    HPoint3f wp = Box.GetWorldPoint(pt);
+    TEST0( Vector4f0(wp - HPoint3f(0,0,0)) );
+  }
+
+  {
+    ContainerBox Box;
+    Box.SetOrigin(Vector3f(-1.5,-1.5,-1.5));
+    Box.SetSize(Vector3f(3,3,3));
+    HPoint3f pt = Box.GetLocalPoint(HPoint3f(0,0,0));
+    HPoint3f wp = Box.GetWorldPoint(pt);
+    TEST0( Vector4f0(wp - HPoint3f(0,0,0)) );
+
+    pt = HPoint3f(1,1,1);
+    wp = Box.GetWorldPoint(pt);
+    TEST0( Vector4f0(wp - HPoint3f(1.5, 1.5, 1.5)) ); 
+
+    Box.EulerYZYRotate(Vector3f(M_PI,0,0));
+    pt = HPoint3f(1,1,1);
+    wp = Box.GetWorldPoint(pt);
+    TEST0( Vector4f0(wp - HPoint3f(-1.5, 1.5, -1.5)) ); 
+
+  }
+
+
+
+
+
+
   END_TESTING;
 }
 

src/Math/testing/testing-prototype.h

@@ -31,8 +31,8 @@
 static int _fail = 0;                      \
 printf("..:: Testing " #name " ::..\n");
 
-#define TEST1(val) _fail += (val)==0
-#define TEST0(val) _fail += (val)!=0
+#define TEST1(val) if ((val)==0) { printf("Assertion failed: %s != 0\n", #val); _fail++; }
+#define TEST0(val) if ((val)!=0) { printf("Assertion failed: %s != 0\n", #val); _fail++; }
 #define END_TESTING return _fail;
 
 #define ASSERT_EQ(a, b) if ((a) != (b)) { printf("Assertion failed: %s != %s\n", #a, #b); _fail++; }

src/Vtk/CMakeLists.txt

@@ -3,6 +3,7 @@ set(HEADERS uLibVtkInterface.h
             vtkContainerBox.h
             vtkHandlerWidget.h
             vtkQViewport.h
+            vtkPolydata.h
             )
 
 set(SOURCES uLibVtkInterface.cxx
@@ -10,6 +11,7 @@ set(SOURCES uLibVtkInterface.cxx
             vtkContainerBox.cpp
             vtkHandlerWidget.cpp
             vtkQViewport.cpp
+            vtkPolydata.cpp
             )
 
 ## Pull in Math VTK wrappers (sets MATH_SOURCES / MATH_HEADERS)
@@ -27,11 +29,17 @@ add_subdirectory(HEP/MuonTomography)
 list(APPEND SOURCES ${HEP_MUONTOMOGRAPHY_SOURCES})
 list(APPEND HEADERS ${HEP_MUONTOMOGRAPHY_HEADERS})
 
+## Pull in HEP/Geant VTK wrappers (sets HEP_GEANT_SOURCES / HEADERS)
+add_subdirectory(HEP/Geant)
+list(APPEND SOURCES ${HEP_GEANT_SOURCES})
+list(APPEND HEADERS ${HEP_GEANT_HEADERS})
+
 set(LIBRARIES Eigen3::Eigen
               ${ROOT_LIBRARIES}
               ${VTK_LIBRARIES}
               ${PACKAGE_LIBPREFIX}Math
-              ${PACKAGE_LIBPREFIX}Detectors)
+              ${PACKAGE_LIBPREFIX}Detectors
+              ${PACKAGE_LIBPREFIX}Geant)
 
 if(USE_CUDA)
   find_package(CUDAToolkit REQUIRED)
@@ -59,6 +67,7 @@ install(FILES ${HEADERS} DESTINATION ${INSTALL_INC_DIR}/Vtk)
 install(FILES ${MATH_HEADERS} DESTINATION ${INSTALL_INC_DIR}/Vtk/Math)
 install(FILES ${HEP_DETECTORS_HEADERS} DESTINATION ${INSTALL_INC_DIR}/Vtk/HEP/Detectors)
 install(FILES ${HEP_MUONTOMOGRAPHY_HEADERS} DESTINATION ${INSTALL_INC_DIR}/Vtk/HEP/MuonTomography)
+install(FILES ${HEP_GEANT_HEADERS} DESTINATION ${INSTALL_INC_DIR}/Vtk/HEP/Geant)
 
 if(BUILD_TESTING)
   include(uLibTargetMacros)

src/Vtk/HEP/Detectors/vtkMuonEvent.h

@@ -45,6 +45,7 @@
 
 #include "HEP/Detectors/MuonEvent.h"
 #include "Vtk/uLibVtkInterface.h"
+#include "Vtk/vtkPolydata.h"
 
 namespace uLib {
 namespace Vtk {

src/Vtk/HEP/Detectors/vtkMuonScatter.h

@@ -46,6 +46,7 @@
 
 #include "HEP/Detectors/MuonScatter.h"
 #include "Vtk/uLibVtkInterface.h"
+#include "Vtk/vtkPolydata.h"
 
 class vtkRenderWindowInteractor;
 

src/Vtk/HEP/Geant/CMakeLists.txt

@@ -0,0 +1,17 @@
+
+################################################################################
+##### Vtk/HEP/Geant - VTK wrappers for HEP Geant objects #######################
+################################################################################
+
+set(HEP_GEANT_SOURCES
+    ${CMAKE_CURRENT_SOURCE_DIR}/vtkGeantEvent.cpp
+    PARENT_SCOPE)
+
+set(HEP_GEANT_HEADERS
+    ${CMAKE_CURRENT_SOURCE_DIR}/vtkGeantEvent.h
+    PARENT_SCOPE)
+
+if(BUILD_TESTING)
+  include(uLibTargetMacros)
+  add_subdirectory(testing)
+endif()

src/Vtk/HEP/Geant/testing/CMakeLists.txt

@@ -0,0 +1,14 @@
+# TESTS
+set(TESTS
+       vtkGeantEventTest
+)
+
+set(LIBRARIES
+       ${PACKAGE_LIBPREFIX}Core
+       ${PACKAGE_LIBPREFIX}Math
+       ${PACKAGE_LIBPREFIX}Geant
+       ${PACKAGE_LIBPREFIX}Vtk
+       Eigen3::Eigen
+)
+
+uLib_add_tests(VtkGeant)

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

@@ -0,0 +1,150 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// CMT Cosmic Muon Tomography project //////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+  Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova
+  All rights reserved
+
+  Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it >
+
+//////////////////////////////////////////////////////////////////////////////*/
+
+#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 "Vtk/uLibVtkViewer.h"
+#include "Vtk/HEP/Geant/vtkGeantEvent.h"
+#include "Vtk/vtkContainerBox.h"
+
+#include <vtkSmartPointer.h>
+#include <vtkCallbackCommand.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkProperty.h>
+#include <vtkRenderer.h>
+#include <vtkRenderWindow.h>
+
+#include <Geant4/Randomize.hh>
+#include <Geant4/G4ParticleGun.hh>
+#include <Geant4/G4SystemOfUnits.hh>
+
+#include <iostream>
+
+using namespace uLib;
+
+// Custom emitter to fire random muons towards the cube
+class RandomEmitter : public Geant::EmitterPrimary {
+public:
+    virtual void GeneratePrimaries(G4Event* anEvent) override {
+        // Start from a random point on the top face of the world box (z = 15m)
+        double x = 0_m; 
+        double y = 0_m; 
+        double z = 14.9_m; // Top face
+        
+        fParticleGun->SetParticlePosition(G4ThreeVector(x, y, z));
+        
+        // Aim at a random point on the bottom face (z = -15m)
+        double tx = 0_m;
+        double ty = 0_m;
+        double tz = -14.9_m; // Bottom face
+        
+        G4ThreeVector dir(tx - x, ty - y, tz - z);
+        fParticleGun->SetParticleMomentumDirection(dir.unit());
+        fParticleGun->SetParticleEnergy(15_GeV);
+        fParticleGun->GeneratePrimaryVertex(anEvent);
+    }
+};
+
+struct AppState {
+    Geant::Scene* scene;
+    Vtk::Viewer* viewer;
+    Vector<Geant::GeantEvent> results;
+};
+
+void KeyPressCallbackFunction(vtkObject* caller, long unsigned int eventId, void* clientData, void* callData) {
+    auto* interactor = static_cast<vtkRenderWindowInteractor*>(caller);
+    auto* state = static_cast<AppState*>(clientData);
+    
+    std::string key = interactor->GetKeySym();
+    if (key == "Return") {
+        std::cout << "--> Firing random muon..." << std::endl;
+        
+        // Run one event
+        state->scene->RunSimulation(1, state->results);
+        
+        if (!state->results.empty()) {
+            // Get the last event
+            Geant::GeantEvent* lastEvent = &state->results.back();
+            std::cout << "    Collected event " << lastEvent->GetEventID() 
+                      << " with " << lastEvent->Path().size() << " steps." << std::endl;
+                      
+            // Wrap it for VTK
+            Vtk::vtkGeantEvent* vtkEvent = new Vtk::vtkGeantEvent(lastEvent);
+            state->viewer->AddPuppet(*vtkEvent);
+            
+            // Re-render
+            state->viewer->GetRenderer()->Render();
+            state->viewer->GetRenderWindow()->Render();
+        }
+        else {
+            std::cout << "    No event collected." << std::endl;
+        }
+    }
+}
+
+int main(int argc, char** argv) {
+    // 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(10_m, 10_m, 10_m));
+    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);
+
+    RandomEmitter* emitter = new RandomEmitter();
+    scene.SetEmitter(emitter);
+    scene.Initialize();
+
+    // 2. Setup VTK Viewer
+    Vtk::Viewer viewer;
+    viewer.GetRenderer()->SetBackground(0.05, 0.05, 0.1);
+    
+    // Visualize world box
+    Vtk::vtkContainerBox* vtkWorld = new Vtk::vtkContainerBox(scene.GetWorldBox());
+    // vtkWorld->ShowBoundingBox(true);
+    vtkWorld->ShowScaleMeasures(true);
+    viewer.AddPuppet(*vtkWorld);
+    
+    // Visualize iron cube
+    Vtk::vtkContainerBox* vtkIron = new Vtk::vtkContainerBox(&iron_box);
+    vtkIron->SetOpacity(0.2);
+    vtkIron->SetRepresentation(Vtk::Puppet::Surface);
+    viewer.AddPuppet(*vtkIron);
+    
+    // 3. Event Handling
+    AppState state = { &scene, &viewer, {} };
+    
+    vtkSmartPointer<vtkCallbackCommand> keyCallback = vtkSmartPointer<vtkCallbackCommand>::New();
+    keyCallback->SetCallback(KeyPressCallbackFunction);
+    keyCallback->SetClientData(&state);
+    
+    viewer.GetInteractor()->AddObserver(vtkCommand::KeyPressEvent, keyCallback);
+
+    std::cout << "=================================================" << std::endl;
+    std::cout << " Geant Muon Simulation Viewer" << std::endl;
+    std::cout << " Press [ENTER] to fire a new random muon" << std::endl;
+    std::cout << " Press [q] to exit" << std::endl;
+    std::cout << "=================================================" << std::endl;
+
+    viewer.ZoomAuto();
+    viewer.Start();
+
+    return 0;
+}

src/Vtk/HEP/Geant/vtkGeantEvent.cpp

@@ -0,0 +1,108 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// CMT Cosmic Muon Tomography project //////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+  Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova
+  All rights reserved
+
+  Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it >
+
+  ------------------------------------------------------------------
+  This library is free software;  you  can  redistribute  it  and/or
+  modify it  under the  terms  of  the  GNU  Lesser  General  Public
+  License as published  by  the  Free  Software  Foundation;  either
+  version 3.0 of the License, or (at your option) any later version.
+
+  This library is  distributed in  the hope that it will  be useful,
+  but  WITHOUT ANY WARRANTY;  without  even  the implied warranty of
+  MERCHANTABILITY  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+  Lesser General Public License for more details.
+
+  You should have received a copy of  the GNU Lesser General  Public
+  License along with this library.
+
+//////////////////////////////////////////////////////////////////////////////*/
+
+#include "vtkGeantEvent.h"
+
+#include <vtkPolyData.h>
+#include <vtkPoints.h>
+#include <vtkCellArray.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkProperty.h>
+#include <vtkSmartPointer.h>
+
+namespace uLib {
+namespace Vtk {
+
+vtkGeantEvent::vtkGeantEvent(Content *content)
+    : m_MuonPath(vtkActor::New()), m_Content(content) {
+  this->InstallPipe();
+  this->Update();
+}
+
+vtkGeantEvent::~vtkGeantEvent() {
+  m_MuonPath->Delete();
+}
+
+vtkPolyData *vtkGeantEvent::GetPolyData() const {
+  if (!m_MuonPath || !m_MuonPath->GetMapper())
+    return NULL;
+  return vtkPolyData::SafeDownCast(m_MuonPath->GetMapper()->GetInput());
+}
+
+void vtkGeantEvent::Update() {
+  if (!m_Content)
+    return;
+
+  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
+  vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();
+
+  const HLine3f &gen = m_Content->GenVector();
+  const Vector<Geant::GeantEvent::Delta> &path = m_Content->Path();
+
+  // Start point from generation vector
+  Vector3f currentPos(gen.origin(0), gen.origin(1), gen.origin(2));
+  points->InsertNextPoint(currentPos(0), currentPos(1), currentPos(2));
+
+  for (size_t i = 0; i < path.size(); ++i) {
+    const Geant::GeantEvent::Delta &delta = path[i];
+    
+    // P_{i+1} = P_i + Length * Direction
+    // Note: HVector3f is stored as (x,y,z,0) in HPoint3f template
+    Vector3f dir(delta.Direction()(0), delta.Direction()(1), delta.Direction()(2));
+    currentPos += delta.GetLength() * dir;
+    
+    points->InsertNextPoint(currentPos(0), currentPos(1), currentPos(2));
+
+    vtkIdType line[2] = {static_cast<vtkIdType>(i),
+                          static_cast<vtkIdType>(i + 1)};
+    lines->InsertNextCell(2, line);
+  }
+
+  vtkPolyData *polyData = GetPolyData();
+  if (polyData) {
+    polyData->SetPoints(points);
+    polyData->SetLines(lines);
+    polyData->Modified();
+  }
+}
+
+void vtkGeantEvent::InstallPipe() {
+  vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
+  vtkSmartPointer<vtkPolyDataMapper> mapper =
+      vtkSmartPointer<vtkPolyDataMapper>::New();
+  mapper->SetInputData(polyData);
+
+  m_MuonPath->SetMapper(mapper);
+  
+  // Set default look: Red line
+  m_MuonPath->GetProperty()->SetColor(1.0, 1.0, 1.0);
+  m_MuonPath->GetProperty()->SetLineWidth(1.0);
+  m_MuonPath->GetProperty()->SetAmbient(1.0);
+
+  this->SetProp(m_MuonPath);
+}
+
+} // namespace Vtk
+} // namespace uLib

src/Vtk/HEP/Geant/vtkGeantEvent.h

@@ -0,0 +1,58 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// CMT Cosmic Muon Tomography project //////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+  Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova
+  All rights reserved
+
+  Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it >
+
+  ------------------------------------------------------------------
+  This library is free software;  you  can  redistribute  it  and/or
+  modify it  under the  terms  of  the  GNU  Lesser  General  Public
+  License as published  by  the  Free  Software  Foundation;  either
+  version 3.0 of the License, or (at your option) any later version.
+
+  This library is  distributed in  the hope that it will  be useful,
+  but  WITHOUT ANY WARRANTY;  without  even  the implied warranty of
+  MERCHANTABILITY  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+  Lesser General Public License for more details.
+
+  You should have received a copy of  the GNU Lesser General  Public
+  License along with this library.
+
+//////////////////////////////////////////////////////////////////////////////*/
+
+#ifndef U_VTKGEANTEVENT_H
+#define U_VTKGEANTEVENT_H
+
+#include "HEP/Geant/GeantEvent.h"
+#include "uLibVtkInterface.h"
+#include "vtkPolydata.h"
+#include <vtkActor.h>
+
+namespace uLib {
+namespace Vtk {
+
+class vtkGeantEvent : public Puppet, public Polydata {
+  typedef Geant::GeantEvent Content;
+
+public:
+  vtkGeantEvent(Content *content);
+  ~vtkGeantEvent();
+
+  virtual class vtkPolyData *GetPolyData() const override;
+
+  virtual void Update();
+
+protected:
+  virtual void InstallPipe();
+
+  vtkActor *m_MuonPath;
+  Content *m_Content;
+};
+
+} // namespace Vtk
+} // namespace uLib
+
+#endif // U_VTKGEANTEVENT_H

src/Vtk/Math/vtkQuadMesh.h

@@ -28,6 +28,7 @@
 
 #include "Math/QuadMesh.h"
 #include "Vtk/uLibVtkInterface.h"
+#include "Vtk/vtkPolydata.h"
 
 class vtkPolyData;
 class vtkActor;

src/Vtk/Math/vtkTriangleMesh.h

@@ -28,6 +28,7 @@
 
 #include "Math/TriangleMesh.h"
 #include "Vtk/uLibVtkInterface.h"
+#include "Vtk/vtkPolydata.h"
 
 class vtkPolyData;
 class vtkActor;

src/Vtk/uLibVtkInterface.cxx

@@ -37,6 +37,7 @@
 #endif
 
 
+#include <string>
 #include <vtkVersion.h>
 #include <vtkProp.h>
 #include <vtkActor.h>
@@ -46,10 +47,13 @@
 #include <vtkPropCollection.h>
 #include <vtkRendererCollection.h>
 #include <vtkPropAssembly.h>
+#include <vtkOutlineSource.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkCubeAxesActor.h>
+#include <vtkRenderer.h>
+#include <vtkProperty.h>
+#include <vtkCamera.h>
 
-#include <vtkPolyData.h>
-#include <vtkPolyDataWriter.h>
-#include <vtkXMLPolyDataWriter.h>
 
 #include "uLibVtkInterface.h"
 
@@ -64,46 +68,6 @@ namespace Vtk {
 
 
 
-////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////
-// POLYDATA //
-
-
-void Polydata::SaveToFile(const char *vtk_file)
-{
-    vtkSmartPointer<vtkPolyDataWriter> writer =
-      vtkSmartPointer<vtkPolyDataWriter>::New();
-    writer->SetFileName(vtk_file);
-    vtkPolyData * data = GetPolyData();
-    if(data) {
-#       if VTK_MAJOR_VERSION <= 5
-        writer->SetInputConnection(data->GetProducerPort());
-#       else
-        writer->SetInputData(data);
-#       endif
-        writer->Update();
-    }
-}
-
-void Polydata::SaveToXMLFile(const char *vtp_file)
-{
-    vtkSmartPointer<vtkXMLPolyDataWriter> writer =
-      vtkSmartPointer<vtkXMLPolyDataWriter>::New();
-    writer->SetFileName(vtp_file);
-    vtkPolyData * data = GetPolyData();
-    if(data) {
-#       if VTK_MAJOR_VERSION <= 5
-        writer->SetInputConnection(data->GetProducerPort());
-#       else
-        writer->SetInputData(data);
-#       endif
-        writer->Update();
-    }
-}
-
-
-
 ////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////
@@ -115,19 +79,59 @@ class PuppetData {
 public:
     PuppetData() :
         m_Renderers(vtkRendererCollection::New()),
-        m_Assembly(vtkPropAssembly::New())
-    {}
+        m_Assembly(vtkPropAssembly::New()),
+        m_OutlineSource(NULL),
+        m_OutlineActor(NULL),
+        m_CubeAxesActor(NULL),
+        m_ShowBoundingBox(false),
+        m_ShowScaleMeasures(false),
+        m_Representation(-1),
+        m_Opacity(-1.0)
+    {
+        m_Color[0] = m_Color[1] = m_Color[2] = -1.0;
+    }
 
     ~PuppetData() {
         m_Renderers->RemoveAllItems();
         m_Assembly->GetParts()->RemoveAllItems();
         m_Renderers->Delete();
         m_Assembly->Delete();
+        if (m_OutlineSource) m_OutlineSource->Delete();
+        if (m_OutlineActor) m_OutlineActor->Delete();
+        if (m_CubeAxesActor) m_CubeAxesActor->Delete();
     }
 
     // members //
     vtkRendererCollection *m_Renderers;
     vtkPropAssembly       *m_Assembly;
+
+    vtkOutlineSource      *m_OutlineSource;
+    vtkActor              *m_OutlineActor;
+    vtkCubeAxesActor      *m_CubeAxesActor;
+
+    bool m_ShowBoundingBox;
+    bool m_ShowScaleMeasures;
+
+    int m_Representation;
+    double m_Color[3];
+    double m_Opacity;
+
+    void ApplyAppearance(vtkProp *p) {
+        vtkActor *actor = vtkActor::SafeDownCast(p);
+        if (!actor) return;
+
+        if (m_Representation != -1) {
+            actor->GetProperty()->SetRepresentation(m_Representation);
+        }
+
+        if (m_Color[0] != -1.0) {
+            actor->GetProperty()->SetColor(m_Color);
+        }
+
+        if (m_Opacity != -1.0) {
+            actor->GetProperty()->SetOpacity(m_Opacity);
+        }
+    }
 };
 
 // -------------------------------------------------------------------------- //
@@ -155,7 +159,10 @@ vtkProp *Puppet::GetProp()
 
 void Puppet::SetProp(vtkProp *prop)
 {
-    if(prop) d->m_Assembly->AddPart(prop);
+    if(prop) {
+        d->m_Assembly->AddPart(prop);
+        d->ApplyAppearance(prop);
+    }
 }
 
 void Puppet::RemoveProp(vtkProp *prop)
@@ -182,6 +189,12 @@ void Puppet::ConnectRenderer(vtkRenderer *renderer)
             for (int i=0; i<props->GetNumberOfItems(); ++i)
                 renderer->AddActor(props->GetNextProp());
         }
+
+        if (d->m_ShowBoundingBox && d->m_OutlineActor) renderer->AddActor(d->m_OutlineActor);
+        if (d->m_ShowScaleMeasures && d->m_CubeAxesActor) {
+            d->m_CubeAxesActor->SetCamera(renderer->GetActiveCamera());
+            renderer->AddActor(d->m_CubeAxesActor);
+        }
     }
 }
 
@@ -197,6 +210,10 @@ void Puppet::DisconnectRenderer(vtkRenderer *renderer)
             for (int i=0; i<props->GetNumberOfItems(); ++i)
                 renderer->RemoveViewProp(props->GetNextProp());
         }
+
+        if (d->m_ShowBoundingBox && d->m_OutlineActor) renderer->RemoveActor(d->m_OutlineActor);
+        if (d->m_ShowScaleMeasures && d->m_CubeAxesActor) renderer->RemoveActor(d->m_CubeAxesActor);
+
         this->GetRenderers()->RemoveItem(renderer);
     }
 }
@@ -225,6 +242,121 @@ void Puppet::PrintSelf(std::ostream &o) const
     d->m_Renderers->PrintSelf(o,vtkIndent(1));
 }
 
+void Puppet::ShowBoundingBox(bool show)
+{
+    if (d->m_ShowBoundingBox == show) return;
+    d->m_ShowBoundingBox = show;
+    if (show) {
+        if (!d->m_OutlineActor) {
+            d->m_OutlineSource = vtkOutlineSource::New();
+            d->m_OutlineActor = vtkActor::New();
+            vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
+            mapper->SetInputConnection(d->m_OutlineSource->GetOutputPort());
+            d->m_OutlineActor->SetMapper(mapper);
+            d->m_OutlineActor->GetProperty()->SetColor(1.0, 1.0, 1.0);
+        }
+
+        double* bounds = d->m_Assembly->GetBounds();
+        d->m_OutlineSource->SetBounds(bounds);
+        d->m_OutlineSource->Update();
+
+        d->m_Renderers->InitTraversal();
+        for (int i = 0; i < d->m_Renderers->GetNumberOfItems(); ++i) {
+            vtkRenderer *renderer = d->m_Renderers->GetNextItem();
+            renderer->AddActor(d->m_OutlineActor);
+        }
+    } else {
+        if (d->m_OutlineActor) {
+            d->m_Renderers->InitTraversal();
+            for (int i = 0; i < d->m_Renderers->GetNumberOfItems(); ++i) {
+                vtkRenderer *renderer = d->m_Renderers->GetNextItem();
+                renderer->RemoveActor(d->m_OutlineActor);
+            }
+        }
+    }
+}
+
+void Puppet::ShowScaleMeasures(bool show)
+{
+    if (d->m_ShowScaleMeasures == show) return;
+    d->m_ShowScaleMeasures = show;
+    if (show) {
+        if (!d->m_CubeAxesActor) {
+            d->m_CubeAxesActor = vtkCubeAxesActor::New();
+            d->m_CubeAxesActor->SetFlyModeToOuterEdges();
+            d->m_CubeAxesActor->SetUseTextActor3D(1);
+            d->m_CubeAxesActor->GetProperty()->SetColor(1.0, 1.0, 1.0);
+        }
+
+        double* bounds = d->m_Assembly->GetBounds();
+        d->m_CubeAxesActor->SetBounds(bounds);
+
+        d->m_Renderers->InitTraversal();
+        for (int i = 0; i < d->m_Renderers->GetNumberOfItems(); ++i) {
+            vtkRenderer *renderer = d->m_Renderers->GetNextItem();
+            d->m_CubeAxesActor->SetCamera(renderer->GetActiveCamera());
+            renderer->AddActor(d->m_CubeAxesActor);
+        }
+    } else {
+        if (d->m_CubeAxesActor) {
+            d->m_Renderers->InitTraversal();
+            for (int i = 0; i < d->m_Renderers->GetNumberOfItems(); ++i) {
+                vtkRenderer *renderer = d->m_Renderers->GetNextItem();
+                renderer->RemoveActor(d->m_CubeAxesActor);
+            }
+        }
+    }
+}
+
+void Puppet::SetRepresentation(Representation mode)
+{
+    int rep = VTK_SURFACE;
+    switch (mode) {
+        case Points:    rep = VTK_POINTS; break;
+        case Wireframe: rep = VTK_WIREFRAME; break;
+        case Surface:   rep = VTK_SURFACE; break;
+    }
+    d->m_Representation = rep;
+
+    vtkPropCollection *props = d->m_Assembly->GetParts();
+    props->InitTraversal();
+    for (int i = 0; i < props->GetNumberOfItems(); ++i) {
+        d->ApplyAppearance(props->GetNextProp());
+    }
+}
+
+void Puppet::SetRepresentation(const char *mode)
+{
+    std::string s(mode);
+    if (s == "points") SetRepresentation(Points);
+    else if (s == "wireframe") SetRepresentation(Wireframe);
+    else if (s == "shaded" || s == "surface") SetRepresentation(Surface);
+}
+
+void Puppet::SetColor(double r, double g, double b)
+{
+    d->m_Color[0] = r;
+    d->m_Color[1] = g;
+    d->m_Color[2] = b;
+
+    vtkPropCollection *props = d->m_Assembly->GetParts();
+    props->InitTraversal();
+    for (int i = 0; i < props->GetNumberOfItems(); ++i) {
+        d->ApplyAppearance(props->GetNextProp());
+    }
+}
+
+void Puppet::SetOpacity(double alpha)
+{
+    d->m_Opacity = alpha;
+
+    vtkPropCollection *props = d->m_Assembly->GetParts();
+    props->InitTraversal();
+    for (int i = 0; i < props->GetNumberOfItems(); ++i) {
+        d->ApplyAppearance(props->GetNextProp());
+    }
+}
+
 
 
 

src/Vtk/uLibVtkInterface.h

@@ -35,6 +35,8 @@ class vtkPolyData;
 class vtkPropCollection;
 class vtkRenderer;
 class vtkRendererCollection;
+class vtkRenderWindowInteractor;
+
 
 namespace uLib {
 namespace Vtk {
@@ -58,8 +60,28 @@ public:
 
   vtkRendererCollection *GetRenderers() const;
 
+  enum Representation {
+      Points,
+      Wireframe,
+      Surface
+  };
+
   virtual void PrintSelf(std::ostream &o) const;
 
+  void ShowBoundingBox(bool show = true);
+
+  void ShowScaleMeasures(bool show = true); 
+
+  void SetRepresentation(Representation mode);
+
+  void SetRepresentation(const char *mode);
+
+  void SetColor(double r, double g, double b);
+
+  void SetOpacity(double alpha);
+
+
+  
   virtual void ConnectInteractor(class vtkRenderWindowInteractor *interactor) {
     (void)interactor;
   }
@@ -74,18 +96,6 @@ private:
   class PuppetData *d;
 };
 
-class Polydata {
-public:
-  virtual vtkPolyData *GetPolyData() const { return NULL; }
-
-  virtual void SaveToFile(const char *vtk_file);
-
-  virtual void SaveToXMLFile(const char *vtp_file);
-
-protected:
-  virtual ~Polydata() {}
-};
-
 } // namespace Vtk
 } // namespace uLib
 

src/Vtk/uLibVtkViewer.cpp

@@ -56,6 +56,8 @@ public:
 };
 vtkStandardNewMacro(vtkInteractorStyleNoSpin);
 
+
+
 namespace uLib {
 namespace Vtk {
 

src/Vtk/vtkContainerBox.h

@@ -28,6 +28,7 @@
 
 #include "Math/ContainerBox.h"
 #include "uLibVtkInterface.h"
+#include "vtkPolydata.h"
 #include <vtkActor.h>
 
 namespace uLib {

src/Vtk/vtkPolydata.cpp

@@ -0,0 +1,70 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// CMT Cosmic Muon Tomography project //////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+  Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova
+  All rights reserved
+
+  Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it >
+
+  ------------------------------------------------------------------
+  This library is free software;  you  can  redistribute  it  and/or
+  modify it  under the  terms  of  the  GNU  Lesser  General  Public
+  License as published  by  the  Free  Software  Foundation;  either
+  version 3.0 of the License, or (at your option) any later version.
+
+  This library is  distributed in  the hope that it will  be useful,
+  but  WITHOUT ANY WARRANTY;  without  even  the implied warranty of
+  MERCHANTABILITY  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+  Lesser General Public License for more details.
+
+  You should have received a copy of  the GNU Lesser General  Public
+  License along with this library.
+
+//////////////////////////////////////////////////////////////////////////////*/
+
+#include <vtkVersion.h>
+#include <vtkPolyData.h>
+#include <vtkPolyDataWriter.h>
+#include <vtkXMLPolyDataWriter.h>
+#include <vtkSmartPointer.h>
+
+#include "vtkPolydata.h"
+
+namespace uLib {
+namespace Vtk {
+
+void Polydata::SaveToFile(const char *vtk_file)
+{
+    vtkSmartPointer<vtkPolyDataWriter> writer =
+      vtkSmartPointer<vtkPolyDataWriter>::New();
+    writer->SetFileName(vtk_file);
+    vtkPolyData * data = GetPolyData();
+    if(data) {
+#       if VTK_MAJOR_VERSION <= 5
+        writer->SetInputConnection(data->GetProducerPort());
+#       else
+        writer->SetInputData(data);
+#       endif
+        writer->Update();
+    }
+}
+
+void Polydata::SaveToXMLFile(const char *vtp_file)
+{
+    vtkSmartPointer<vtkXMLPolyDataWriter> writer =
+      vtkSmartPointer<vtkXMLPolyDataWriter>::New();
+    writer->SetFileName(vtp_file);
+    vtkPolyData * data = GetPolyData();
+    if(data) {
+#       if VTK_MAJOR_VERSION <= 5
+        writer->SetInputConnection(data->GetProducerPort());
+#       else
+        writer->SetInputData(data);
+#       endif
+        writer->Update();
+    }
+}
+
+} // Vtk
+} // uLib

src/Vtk/vtkPolydata.h

@@ -28,22 +28,25 @@
 #ifndef VTKPOLYDATA_H
 #define VTKPOLYDATA_H
 
-#include "uLibVtkInterface.h"
-
+class vtkPolyData;
 
 namespace uLib {
+namespace Vtk {
 
-class vtkPolyData : public Abstract::uLibVtkPolydata {
-    typedef ::vtkPolyData Content;
+class Polydata {
 public:
+  virtual vtkPolyData *GetPolyData() const { return NULL; }
 
+  virtual void SaveToFile(const char *vtk_file);
 
+  virtual void SaveToXMLFile(const char *vtp_file);
 
+protected:
+  virtual ~Polydata() {}
 };
 
-
-
-}
+} // namespace Vtk
+} // namespace uLib
 
 
 #endif // VTKPOLYDATA_H