refactor: Simplify vtkDetectorChamber by removing redundant transform management and improve vtkHandlerWidget rotation calculation using ray-plane intersection.

CMakeConfig.in.h

@@ -34,7 +34,9 @@
 #cmakedefine HAVE_FLOOR
 
 /* Having Geant4 installed */
+#ifndef HAVE_GEANT4
 #cmakedefine HAVE_GEANT4
+#endif
 
 /* Define to 1 if you have the <inttypes.h> header file. */
 #cmakedefine HAVE_INTTYPES_H

CMakeLists.txt

@@ -3,7 +3,14 @@
 ##### CMAKE LISTS ##############################################################
 ################################################################################
 
+if(EXISTS "${CMAKE_BINARY_DIR}/conan_toolchain.cmake")
+    include("${CMAKE_BINARY_DIR}/conan_toolchain.cmake")
+endif()
+
 cmake_minimum_required (VERSION 3.26)
+
+set(QT_NO_VERSION_CHECK TRUE)
+
 if(POLICY CMP0167)
     cmake_policy(SET CMP0167 NEW)
 endif()
@@ -107,6 +114,8 @@ set(Boost_USE_MULTITHREADED ON)
 set(Boost_USE_STATIC_RUNTIME OFF)
 message(STATUS "CMAKE_PREFIX_PATH is ${CMAKE_PREFIX_PATH}")
 
+find_package(HDF5 REQUIRED CONFIG)
+
 find_package(Boost 1.45.0 COMPONENTS program_options serialization unit_test_framework REQUIRED)
 include_directories(${Boost_INCLUDE_DIRS})
 
@@ -118,15 +127,12 @@ find_package(ROOT CONFIG REQUIRED)
 include(${ROOT_USE_FILE})
 
 find_package(VTK REQUIRED)
-# include(${VTK_USE_FILE})
-
 find_package(pybind11 REQUIRED)
 
-
 option(CENTOS_SUPPORT "VTK definitions for CentOS" OFF)
 if(CENTOS_SUPPORT)
     find_package(VTK CONFIG REQUIRED)
-    include(${VTK_USE_FILE})
+    # include(${VTK_USE_FILE})
 else()
     find_package(VTK REQUIRED
                      COMPONENTS CommonColor
@@ -146,7 +152,33 @@ else()
                                 RenderingFreeType
                                 RenderingGL2PSOpenGL2
                                 RenderingOpenGL2
-                                RenderingVolumeOpenGL2)
+                                RenderingVolumeOpenGL2
+                                IOGeometry
+                                GUISupportQt)
+endif()
+
+find_package(Qt6 COMPONENTS Widgets REQUIRED)
+
+find_package(Geant4)
+if(Geant4_FOUND)
+    message(STATUS "Geant4 libs: ${Geant4_LIBRARIES}")
+    add_compile_definitions(HAVE_GEANT4)
+    set(HAVE_GEANT4 1)
+
+    # Sanitize Geant4 targets to remove Qt5 dependencies that conflict with VTK/Qt6
+    if(TARGET Geant4::G4interfaces)
+        set_target_properties(Geant4::G4interfaces PROPERTIES 
+            INTERFACE_LINK_LIBRARIES "Geant4::G4global;Geant4::G4graphics_reps;Geant4::G4intercoms"
+        )
+    endif()
+    if(TARGET Geant4::G4OpenGL)
+        set_target_properties(Geant4::G4OpenGL PROPERTIES 
+            INTERFACE_LINK_LIBRARIES "Geant4::G4vis_management;Geant4::G4graphics_reps;Geant4::G4geometry;Geant4::G4materials;Geant4::G4intercoms;Geant4::G4global;OpenGL::GL;OpenGL::GLU"
+        )
+    endif()
+else()
+    message(STATUS "Geant4 NOT found - optional features will be disabled")
+    set(HAVE_GEANT4 0)
 endif()
 
 set(CMAKE_REQUIRED_INCLUDES CMAKE_REQUIRED_INCLUDES math.h)
@@ -204,8 +236,8 @@ add_subdirectory(${SRC_DIR}/Core)
 include_directories(${SRC_DIR}/Math)
 add_subdirectory(${SRC_DIR}/Math)
 
-include_directories(${SRC_DIR}/Detectors)
+include_directories(${SRC_DIR}/HEP)
-add_subdirectory(${SRC_DIR}/Detectors)
+add_subdirectory(${SRC_DIR}/HEP)
 
 include_directories(${SRC_DIR}/Root)
 add_subdirectory(${SRC_DIR}/Root)
@@ -215,7 +247,7 @@ add_subdirectory(${SRC_DIR}/Vtk)
 
 add_subdirectory(${SRC_DIR}/Python)
 
-#add_subdirectory("${SRC_DIR}/utils/make_recipe")
+add_subdirectory(app)
 
 ## Documentation and packages
 

app/CMakeLists.txt

@@ -0,0 +1,3 @@
+if(HAVE_GEANT4)
+    add_subdirectory(gcompose)
+endif()

app/gcompose/CMakeLists.txt

@@ -0,0 +1,38 @@
+
+add_executable(gcompose 
+    src/main.cpp
+    src/MainWindow.h
+    src/MainWindow.cpp
+)
+
+set_target_properties(gcompose PROPERTIES
+    AUTOMOC ON
+    AUTOUIC ON
+    AUTORCC ON
+)
+
+target_include_directories(gcompose PRIVATE 
+    ${SRC_DIR}
+    ${PROJECT_BINARY_DIR}
+    ${Geant4_INCLUDE_DIRS}
+    ${VTK_INCLUDE_DIRS}
+)
+
+# Filter Geant4 libraries to remove Qt-dependent ones
+set(Geant4_LIBS_FILTERED ${Geant4_LIBRARIES})
+if(Geant4_LIBS_FILTERED)
+    list(REMOVE_ITEM Geant4_LIBS_FILTERED Geant4::G4interfaces Geant4::G4OpenGL Geant4::G4visQt3D)
+endif()
+
+target_link_libraries(gcompose
+    mutomCore
+    mutomMath
+    mutomGeant
+    mutomVtk
+    ${Geant4_LIBS_FILTERED}
+    ${VTK_LIBRARIES}
+    Qt6::Widgets
+    VTK::GUISupportQt
+)
+
+install(TARGETS gcompose RUNTIME DESTINATION bin)

app/gcompose/src/MainWindow.cpp

@@ -0,0 +1,15 @@
+#include "MainWindow.h"
+#include <Vtk/vtkQViewport.h>
+
+using namespace uLib;
+
+MainWindow::MainWindow(QWidget* parent) : QMainWindow(parent) {
+    m_viewport = new Vtk::QViewport(this);
+    setCentralWidget(m_viewport);
+
+    setWindowTitle("gcompose - Qt VTK Interface");
+    resize(1200, 800);
+}
+
+MainWindow::~MainWindow() {
+}

app/gcompose/src/MainWindow.h

@@ -0,0 +1,25 @@
+#ifndef MAINWINDOW_H
+#define MAINWINDOW_H
+
+#include <QMainWindow>
+#include <QVTKOpenGLNativeWidget.h>
+
+namespace uLib {
+namespace Vtk {
+class QViewport;
+}
+}
+
+class MainWindow : public QMainWindow {
+    Q_OBJECT
+public:
+    MainWindow(QWidget* parent = nullptr);
+    virtual ~MainWindow();
+
+    uLib::Vtk::QViewport* getViewport() { return m_viewport; }
+
+private:
+    uLib::Vtk::QViewport* m_viewport;
+};
+
+#endif

app/gcompose/src/main.cpp

@@ -0,0 +1,67 @@
+#include <QApplication>
+#include "MainWindow.h"
+
+#include "Math/ContainerBox.h"
+#include <HEP/Geant/Scene.h>
+#include "HEP/Detectors/DetectorChamber.h"
+#include "Vtk/HEP/Detectors/vtkDetectorChamber.h"
+
+#include <Vtk/vtkContainerBox.h>
+#include <Vtk/vtkQViewport.h>
+
+#include <vtkSmartPointer.h>
+#include <vtkCubeSource.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkActor.h>
+#include <vtkRenderer.h>
+
+
+#include "Math/Units.h"
+#include <iostream>
+
+using namespace uLib;
+using namespace uLib::literals;
+
+int main(int argc, char** argv) {
+    QApplication app(argc, argv);
+    std::cout << "Starting gcompose Qt application..." << std::endl;
+
+    ContainerBox world_box(Vector3f(1, 1, 1));
+    world_box.Scale(Vector3f(20_mm, 20_mm, 20_mm));
+
+    DetectorChamber d1, d2;
+    d1.SetSize(Vector3f(1, 1, 1));
+    d1.SetPosition(Vector3f(0, 0, 0));
+    d1.Scale(Vector3f(5, 10, 2));
+    d1.Translate(Vector3f(0, 0, 0));
+
+    d2.SetSize(Vector3f(1, 1, 1));
+    d2.SetPosition(Vector3f(0, 0, 0));
+    d2.Scale(Vector3f(5, 10, 2));
+    d2.Translate(Vector3f(0, 0, 10));
+    
+
+    Geant::Scene scene;
+    scene.ConstructWorldBox(world_box.GetSize(), "G4_AIR");
+    scene.Initialize();
+
+    // 2. Initialize MainWindow (contains embedded VTK QViewport)
+    MainWindow window;
+    Vtk::QViewport* viewport = window.getViewport();
+
+    Vtk::vtkDetectorChamber vtk_d1(&d1);
+    viewport->AddPuppet(vtk_d1);
+
+    Vtk::vtkDetectorChamber vtk_d2(&d2);
+    viewport->AddPuppet(vtk_d2);
+
+    Vtk::vtkContainerBox vtk_box(&world_box);
+    viewport->AddPuppet(vtk_box);
+
+    viewport->ZoomAuto();
+
+    std::cout << "Geant4 and VTK scenes are ready." << std::endl;
+    
+    window.show();
+    return app.exec();
+}

conanfile.txt

@@ -1,8 +1,13 @@
 [requires]
 eigen/3.4.0
 boost/1.83.0
-pybind11/3.0.2
+# pybind11/3.0.2
+hdf5/1.14.3
 
 [generators]
 CMakeDeps
 CMakeToolchain
+
+[options]
+hdf5/*:threadsafe=True
+hdf5/*:enable_unsupported=True

src/Core/CMakeLists.txt

@@ -49,7 +49,7 @@ endif()
 target_link_libraries(${libname} ${LIBRARIES})
 
 install(TARGETS ${libname}
-        EXPORT "${PROJECT_NAME}Targets"
+        EXPORT "uLibTargets"
         RUNTIME DESTINATION ${INSTALL_BIN_DIR} COMPONENT bin
         LIBRARY DESTINATION ${INSTALL_LIB_DIR} COMPONENT lib)
 

src/Core/Debug.h

@@ -119,8 +119,8 @@ public:
     void AddAdapter(AdapterInterface &ad) { m_a.push_back(Adapter(ad)); }
 
     void Update() {
-        foreach(Adapter &ad, m_a) {
+        for(Adapter &ad : m_a) {
-            foreach(DItem &item, m_v) {
+            for(DItem &item : m_v) {
                 item.m_adapter->operator()(ad, item.m_value);
             }
         }

src/Detectors/CMakeLists.txt

@@ -1,12 +0,0 @@
-set(HEADERS MuonScatter.h MuonError.h MuonEvent.h)
-
-set(ULIB_SELECTED_MODULES ${ULIB_SELECTED_MODULES} Detectors PARENT_SCOPE)
-
-
-install(FILES ${HEADERS}
-        DESTINATION ${INSTALL_INC_DIR}/Detectors)
-
-if(BUILD_TESTING)
-  include(uLibTargetMacros)
-  add_subdirectory(testing)
-endif()

src/Detectors/Solid.cpp

@@ -1,114 +0,0 @@
-/*//////////////////////////////////////////////////////////////////////////////
-// 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.
-
-//////////////////////////////////////////////////////////////////////////////*/
-
-
-
-// G4 Solid //
-#include <Geant4/G4Material.hh>
-#include <Geant4/G4NistManager.hh>
-#include <Geant4/G4LogicalVolume.hh>
-
-// Tessellated solid //
-#include <Geant4/G4TessellatedSolid.hh>
-#include <Geant4/G4TriangularFacet.hh>
-#include <Geant4/G4ThreeVector.hh>
-
-
-#include "Math/Dense.h"
-
-#include "Solid.h"
-
-namespace uLib {
-
-class DetectorsSolidPimpl {
-public:
-    static G4ThreeVector getG4Vector3f(const Vector3f &vector) {
-        return G4ThreeVector( vector(0), vector(1), vector(2) );
-    }
-};
-
-
-
-Solid::Solid() :
-    m_Logical (new G4LogicalVolume(NULL,NULL,"unnamed_solid")),
-    m_Material(NULL)
-{}
-
-Solid::Solid(const char *name) :
-    m_Logical(new G4LogicalVolume(NULL,NULL,name)),
-    m_Material(NULL)
-{}
-
-
-
-void Solid::SetNistMaterial(const char *name)
-{
-    G4NistManager *nist = G4NistManager::Instance();
-    if (m_Material) delete m_Material;
-    m_Material = nist->FindOrBuildMaterial(name);
-    m_Logical->SetMaterial(m_Material);
-}
-
-void Solid::SetMaterial(G4Material *material)
-{
-    if(material)
-    {
-        m_Material = material;
-        m_Logical->SetMaterial(material);
-    }
-}
-
-
-
-
-
-
-
-
-
-
-TessellatedSolid::TessellatedSolid(const char *name) :
-    BaseClass(name),
-    m_Solid(new G4TessellatedSolid(name))
-{}
-
-
-void TessellatedSolid::SetMesh(TriangleMesh &mesh)
-{
-    G4TessellatedSolid *ts = this->m_Solid;
-    for (int i=0; i<mesh.Triangles().size(); ++i) {
-        const Vector3i &trg = mesh.Triangles().at(i);
-        G4TriangularFacet *facet = new G4TriangularFacet(
-                    DetectorsSolidPimpl::getG4Vector3f(mesh.Points().at(trg(0))),
-                    DetectorsSolidPimpl::getG4Vector3f(mesh.Points().at(trg(1))),
-                    DetectorsSolidPimpl::getG4Vector3f(mesh.Points().at(trg(2))),
-                    ABSOLUTE);
-        ts->AddFacet((G4VFacet *)facet);
-    }
-    this->m_Logical->SetSolid(ts);
-}
-
-
-}

src/HEP/CMakeLists.txt

@@ -0,0 +1,13 @@
+
+################################################################################
+##### HEP - High Energy Physics modules ########################################
+################################################################################
+
+include_directories(${SRC_DIR}/HEP)
+
+add_subdirectory(Detectors)
+add_subdirectory(Geant)
+# add_subdirectory(MuonTomography)
+
+set(ULIB_SHARED_LIBRARIES ${ULIB_SHARED_LIBRARIES} PARENT_SCOPE)
+set(ULIB_SELECTED_MODULES ${ULIB_SELECTED_MODULES} PARENT_SCOPE)

src/HEP/Detectors/CMakeLists.txt

@@ -0,0 +1,34 @@
+set(HEADERS
+    ChamberHitEvent.h
+    DetectorChamber.h
+    ExperimentFitEvent.h
+    HierarchicalEncoding.h
+    Hit.h
+    HitMC.h
+    LinearFit.h
+    MuonError.h
+    MuonEvent.h
+    MuonScatter.h
+)
+
+set(libname ${PACKAGE_LIBPREFIX}Detectors)
+set(ULIB_SHARED_LIBRARIES ${ULIB_SHARED_LIBRARIES} ${libname} PARENT_SCOPE)
+set(ULIB_SELECTED_MODULES ${ULIB_SELECTED_MODULES} Detectors PARENT_SCOPE)
+
+## Headers-only INTERFACE library
+add_library(${libname} INTERFACE)
+target_include_directories(${libname} INTERFACE
+    $<BUILD_INTERFACE:${SRC_DIR}>
+    $<INSTALL_INTERFACE:${INSTALL_INC_DIR}>
+)
+
+install(TARGETS ${libname}
+        EXPORT "uLibTargets")
+
+install(FILES ${HEADERS}
+        DESTINATION ${INSTALL_INC_DIR}/HEP/Detectors)
+
+if(BUILD_TESTING)
+  include(uLibTargetMacros)
+  add_subdirectory(testing)
+endif()

src/HEP/Detectors/ChamberHitEvent.h

@@ -29,8 +29,8 @@
 #define U_CHAMBERHITEVENT_H
 
 #include "Core/Vector.h"
-#include "Hit.h"
+#include "Detectors/HitMC.h"
-#include "ChamberDetector.h"
+#include "Detectors/DetectorChamber.h"
 
 namespace uLib {
 
@@ -38,17 +38,17 @@ class ChamberHitEventData
 {
 public:
     uLibConstRefMacro (Hits, Vector<HitData>    )
-    uLibGetMacro      (Idv,  ChamberDetector::ID)
+    uLibGetMacro      (Idv,  uint)
 private:
     friend class ChamberHitEvent;
     Vector<HitData>     m_Hits;
-    DetectorChamber::ID m_Idv; // -> chamber/view
+    uint                m_Idv; // -> chamber/view
 };
 
 class ChamberHitEvent : public ChamberHitEventData {
 public:
     uLibRefMacro  (Hits, Vector<HitData>    )
-    uLibSetMacro  (Idv,  ChamberDetector::ID)
+    uLibSetMacro  (Idv,  uint)
 };
 
 }

src/HEP/Detectors/Hit.h

@@ -33,9 +33,6 @@
 namespace uLib {
 
 
-
-
-
 class HitRawCode_CMSDrift :
         public BitCode4<unsigned short,6,3,2,5>
 {
@@ -59,8 +56,14 @@ public:
 };
 
 
+class HitData {
+public:
+    HitData() {}
+    ~HitData() {}
     
     
+};
+
 
 
 

src/HEP/Geant/ActionInitialization.cpp

@@ -0,0 +1,39 @@
+#include "ActionInitialization.hh"
+#include "EmitterPrimary.hh"
+#include "SteppingAction.hh"
+
+namespace uLib {
+namespace Geant {
+
+ActionInitialization::ActionInitialization(EmitterPrimary *emitter,
+                                           Vector<GeantEvent> *output)
+    : G4VUserActionInitialization(),
+      m_Emitter(emitter),
+      m_Output(output)
+{}
+
+ActionInitialization::~ActionInitialization() {}
+
+void ActionInitialization::BuildForMaster() const {
+  // Master thread: no per-event actions needed
+}
+
+void ActionInitialization::Build() const {
+  // Register the primary generator
+  if (m_Emitter) {
+    SetUserAction(m_Emitter);
+  } else {
+    // Fallback: default EmitterPrimary
+    SetUserAction(new EmitterPrimary());
+  }
+
+  // Register actions
+  if (m_Output) {
+    SteppingAction *sa = new SteppingAction(m_Output);
+    SetUserAction(static_cast<G4UserSteppingAction*>(sa));
+    SetUserAction(static_cast<G4UserEventAction*>(sa));
+  }
+}
+
+} // namespace Geant
+} // namespace uLib

src/HEP/Geant/ActionInitialization.hh

@@ -0,0 +1,35 @@
+#ifndef ActionInitialization_h
+#define ActionInitialization_h
+
+#include "G4VUserActionInitialization.hh"
+#include "Core/Vector.h"
+
+namespace uLib {
+namespace Geant {
+
+class EmitterPrimary;
+class GeantEvent;
+
+class ActionInitialization : public G4VUserActionInitialization {
+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);
+  ~ActionInitialization();
+
+  // Metodo chiamato solo dal thread principale (Master)
+  virtual void BuildForMaster() const;
+
+  // 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
+} // namespace uLib
+
+#endif

src/HEP/Geant/CMakeLists.txt

@@ -0,0 +1,73 @@
+
+################################################################################
+##### HEP/Geant - Geant4 integration library ###################################
+################################################################################
+
+if(NOT Geant4_FOUND)
+    message(STATUS "Geant4 not found - skipping mutomGeant library")
+    return()
+endif()
+
+message(STATUS "Geant4 found: ${Geant4_VERSION}")
+include(${Geant4_USE_FILE})
+
+set(HEADERS
+    GeantEvent.h
+    Matter.h
+    Scene.h
+    Solid.h
+    DetectorConstruction.hh
+    PhysicsList.hh
+    ActionInitialization.hh
+    SteppingAction.hh
+)
+
+set(SOURCES
+    Scene.cpp
+    Solid.cpp
+    EmitterPrimary.cpp
+    DetectorConstruction.cpp
+    PhysicsList.cpp
+    ActionInitialization.cpp
+    SteppingAction.cpp
+)
+
+set(libname ${PACKAGE_LIBPREFIX}Geant)
+set(ULIB_SHARED_LIBRARIES ${ULIB_SHARED_LIBRARIES} ${libname} PARENT_SCOPE)
+set(ULIB_SELECTED_MODULES ${ULIB_SELECTED_MODULES} Geant PARENT_SCOPE)
+
+add_library(${libname} SHARED ${SOURCES})
+set_target_properties(${libname} PROPERTIES
+                      VERSION ${PROJECT_VERSION}
+                      SOVERSION ${PROJECT_SOVERSION})
+
+target_include_directories(${libname} PRIVATE ${Geant4_INCLUDE_DIRS})
+
+target_link_libraries(${libname}
+    ${PACKAGE_LIBPREFIX}Core
+    ${PACKAGE_LIBPREFIX}Math
+    ${PACKAGE_LIBPREFIX}Detectors
+)
+
+# Filter Geant4 libraries to remove Qt-dependent ones
+set(Geant4_LIBS_FILTERED ${Geant4_LIBRARIES})
+if(Geant4_LIBS_FILTERED)
+    list(REMOVE_ITEM Geant4_LIBS_FILTERED Geant4::G4interfaces Geant4::G4OpenGL Geant4::G4visQt3D)
+endif()
+
+target_link_libraries(${libname}
+    ${Geant4_LIBS_FILTERED}
+)
+
+install(TARGETS ${libname}
+        EXPORT "uLibTargets"
+        RUNTIME DESTINATION ${INSTALL_BIN_DIR} COMPONENT bin
+        LIBRARY DESTINATION ${INSTALL_LIB_DIR} COMPONENT lib)
+
+install(FILES ${HEADERS}
+        DESTINATION ${INSTALL_INC_DIR}/HEP/Geant)
+
+if(BUILD_TESTING)
+  include(uLibTargetMacros)
+  add_subdirectory(testing)
+endif()

src/HEP/Geant/DetectorConstruction.cpp

@@ -0,0 +1,26 @@
+
+#include "DetectorConstruction.hh"
+
+#include "Core/Object.h"
+#include "Math/ContainerBox.h"
+
+#include "G4Box.hh"
+#include "G4LogicalVolume.hh"
+#include "G4NistManager.hh"
+#include "G4PVPlacement.hh"
+#include "G4RunManager.hh"
+#include "G4SystemOfUnits.hh"
+
+namespace uLib {
+namespace Geant {
+
+DetectorConstruction::DetectorConstruction(const char *name) : G4VUserDetectorConstruction() {}
+
+DetectorConstruction::~DetectorConstruction() {}
+
+G4VPhysicalVolume *DetectorConstruction::Construct() { return nullptr; }
+
+void DetectorConstruction::ConstructSDandField() {}
+
+} // namespace Geant
+} // namespace uLib

src/HEP/Geant/DetectorConstruction.hh

@@ -0,0 +1,30 @@
+#ifndef DetectorConstruction_h
+#define DetectorConstruction_h
+
+#include "Core/Object.h"
+#include "G4VUserDetectorConstruction.hh"
+#include "globals.hh"
+
+class G4VPhysicalVolume;
+class G4LogicalVolume;
+
+namespace uLib {
+namespace Geant {
+
+class DetectorConstruction : public G4VUserDetectorConstruction {
+public:
+  DetectorConstruction(const char *name);
+  virtual ~DetectorConstruction();
+
+  virtual G4VPhysicalVolume *Construct();
+
+  virtual void ConstructSDandField();
+};
+
+} // namespace Geant
+} // namespace uLib
+
+#endif
+
+
+

src/HEP/Geant/EmitterPrimary.cpp

@@ -0,0 +1,155 @@
+#include "EmitterPrimary.hh"
+
+#include "G4Box.hh"
+#include "G4LogicalVolume.hh"
+#include "G4LogicalVolumeStore.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4ParticleGun.hh"
+#include "G4ParticleTable.hh"
+#include "G4RunManager.hh"
+#include "G4SystemOfUnits.hh"
+#include "Randomize.hh"
+
+namespace uLib {
+namespace Geant {
+
+EmitterPrimary::EmitterPrimary()
+    : G4VUserPrimaryGeneratorAction(), fParticleGun(nullptr) {
+  // Creiamo il ParticleGun impostandolo per sparare 1 particella alla volta
+  G4int n_particle = 1;
+  fParticleGun = new G4ParticleGun(n_particle);
+
+  // Otteniamo la tabella delle particelle di Geant4
+  G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
+
+  // Cerchiamo il muone negativo (usa "mu+" per l'antimuone)
+  G4String particleName = "mu-";
+  G4ParticleDefinition *particle = particleTable->FindParticle(particleName);
+
+  // Configuriamo le proprietà iniziali della particella
+  fParticleGun->SetParticleDefinition(particle);
+
+  // Impostiamo la direzione della quantità di moto (es. lungo l'asse Z)
+  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., -1.));
+
+  // Impostiamo l'energia cinetica a 1 GeV
+  fParticleGun->SetParticleEnergy(1.0 * GeV);
+
+  // Impostiamo la posizione di partenza (origine)
+  fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., 10. * m));
+}
+
+EmitterPrimary::~EmitterPrimary() {
+  // Importante: liberare la memoria
+  delete fParticleGun;
+}
+
+void EmitterPrimary::GeneratePrimaries(G4Event *anEvent) {
+  // Questo metodo viene invocato all'inizio di ogni evento.
+  // Qui potresti anche aggiungere una randomizzazione della posizione o
+  // dell'energia.
+
+  fParticleGun->GeneratePrimaryVertex(anEvent);
+}
+
+// -------------------------------------------------------------------------- //
+
+QuadMeshEmitterPrimary::QuadMeshEmitterPrimary()
+    : EmitterPrimary(), m_Mesh(nullptr), m_TotalArea(0.0) {
+}
+
+QuadMeshEmitterPrimary::~QuadMeshEmitterPrimary() {
+}
+
+void QuadMeshEmitterPrimary::SetMesh(uLib::QuadMesh *mesh) {
+  m_Mesh = mesh;
+  CalculateAreas();
+}
+
+void QuadMeshEmitterPrimary::CalculateAreas() {
+  if (!m_Mesh) return;
+  m_CumulativeAreas.clear();
+  m_TotalArea = 0.0;
+
+  const auto &points = m_Mesh->Points();
+  const auto &quads = m_Mesh->Quads();
+
+  for (const auto &q : quads) {
+    uLib::Vector3f v0 = points[q(0)];
+    uLib::Vector3f v1 = points[q(1)];
+    uLib::Vector3f v2 = points[q(2)];
+    uLib::Vector3f v3 = points[q(3)];
+
+    double a1 = 0.5 * (v1 - v0).cross(v2 - v0).norm();
+    double a2 = 0.5 * (v2 - v0).cross(v3 - v0).norm();
+    m_TotalArea += (a1 + a2);
+    m_CumulativeAreas.push_back(m_TotalArea);
+  }
+}
+
+void QuadMeshEmitterPrimary::GeneratePrimaries(G4Event *anEvent) {
+  if (!m_Mesh || m_TotalArea <= 0.0) return;
+
+  // 1. Choose a quad
+  double r = G4UniformRand() * m_TotalArea;
+  auto it = std::lower_bound(m_CumulativeAreas.begin(), m_CumulativeAreas.end(), r);
+  int quadIdx = std::distance(m_CumulativeAreas.begin(), it);
+
+  const auto &q = m_Mesh->Quads()[quadIdx];
+  const auto &points = m_Mesh->Points();
+  uLib::Vector3f v0 = points[q(0)];
+  uLib::Vector3f v1 = points[q(1)];
+  uLib::Vector3f v2 = points[q(2)];
+  uLib::Vector3f v3 = points[q(3)];
+
+  // 2. Choose a point on the quad
+  double a1 = 0.5 * (v1 - v0).cross(v2 - v0).norm();
+  double a2 = 0.5 * (v2 - v0).cross(v3 - v0).norm();
+  
+  G4ThreeVector pos;
+  uLib::Vector3f normal = m_Mesh->GetNormal(quadIdx);
+  
+  if (G4UniformRand() < a1 / (a1 + a2)) {
+    double u = G4UniformRand();
+    double v = G4UniformRand();
+    if (u + v > 1.0) { u = 1.0 - u; v = 1.0 - v; }
+    uLib::Vector3f p = v0 + u * (v1 - v0) + v * (v2 - v0);
+    pos.set(p(0), p(1), p(2));
+  } else {
+    double u = G4UniformRand();
+    double v = G4UniformRand();
+    if (u + v > 1.0) { u = 1.0 - u; v = 1.0 - v; }
+    uLib::Vector3f p = v0 + u * (v2 - v0) + v * (v3 - v0);
+    pos.set(p(0), p(1), p(2));
+  }
+
+  // 3. Choose a direction (Cosmic Muon: cos^2(theta))
+  G4ThreeVector dir;
+  bool accepted = false;
+  int tries = 0;
+  
+  while (!accepted && tries < 1000) {
+    tries++;
+    double cosTheta = std::pow(G4UniformRand(), 1.0/3.0);
+    double sinTheta = std::sqrt(1.0 - cosTheta * cosTheta);
+    double phi = 2.0 * M_PI * G4UniformRand();
+    
+    // Incoming from above (+Z towards -Z)
+    dir.set(sinTheta * std::cos(phi), sinTheta * std::sin(phi), -cosTheta);
+    
+    // Filtering: pointing on the same side of the face normal
+    if (dir.x() * normal(0) + dir.y() * normal(1) + dir.z() * normal(2) > 0) {
+      accepted = true;
+    }
+  }
+
+  if (accepted) {
+    fParticleGun->SetParticlePosition(pos);
+    fParticleGun->SetParticleMomentumDirection(dir);
+    // Keep energy from base class or set here if needed
+    fParticleGun->GeneratePrimaryVertex(anEvent);
+  }
+}
+
+} // namespace Geant
+} // namespace uLib

src/HEP/Geant/EmitterPrimary.hh

@@ -0,0 +1,54 @@
+#ifndef U_GEANT_EMITTERPRIMARY_HH
+#define U_GEANT_EMITTERPRIMARY_HH 1
+
+#include "G4VUserPrimaryGeneratorAction.hh"
+#include "globals.hh"
+
+#include "Math/QuadMesh.h"
+#include <vector> // Added for std::vector
+
+class G4ParticleGun;
+class G4Event;
+
+namespace uLib {
+namespace Geant {
+
+class EmitterPrimary : public G4VUserPrimaryGeneratorAction
+{
+  public:
+    EmitterPrimary();    
+    virtual ~EmitterPrimary();
+
+    // Metodo principale chiamato all'inizio di ogni evento
+    virtual void GeneratePrimaries(G4Event*);         
+
+  protected:
+    G4ParticleGun* fParticleGun; // Puntatore al cannone di particelle
+};
+
+
+class QuadMeshEmitterPrimary : public EmitterPrimary
+{
+  public:
+    QuadMeshEmitterPrimary();    
+    virtual ~QuadMeshEmitterPrimary();
+
+    // Metodo principale chiamato all'inizio di ogni evento
+    virtual void GeneratePrimaries(G4Event*);
+
+    void SetMesh(uLib::QuadMesh* mesh);
+
+  private:
+    uLib::QuadMesh* m_Mesh;
+    std::vector<double> m_CumulativeAreas;
+    double m_TotalArea;
+
+    void CalculateAreas();
+};
+
+} // namespace Geant
+} // namespace uLib
+
+
+
+#endif

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

@@ -0,0 +1,92 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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_GEANTEVENT_H
+#define U_GEANTEVENT_H
+
+#include "Core/Types.h"
+#include "Core/Vector.h"
+#include "Math/Dense.h"
+
+#include <string>
+
+namespace uLib {
+namespace Geant {
+
+// 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:
+
+  /// 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:
+  Id_t m_EventID;
+  Scalarf m_Momentum;
+  HLine3f m_GenVector;
+  Vector<Delta> m_Path;
+
+public:
+  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 // U_GEANTEVENT_H

src/HEP/Geant/Matter.h

@@ -34,6 +34,7 @@ class G4Element;
 class G4Material;
 
 namespace uLib {
+namespace Geant {
 
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -65,6 +66,7 @@ private:
 
 
 
+}
 }
 
 

src/HEP/Geant/PhysicsList.cpp

@@ -0,0 +1,28 @@
+#include "PhysicsList.hh"
+
+#include "G4DecayPhysics.hh"
+#include "G4EmStandardPhysics.hh"
+#include "G4RadioactiveDecayPhysics.hh"
+
+namespace uLib {
+namespace Geant {
+
+PhysicsList::PhysicsList() : G4VModularPhysicsList() {
+  SetVerboseLevel(1);
+
+  // Default physics
+  RegisterPhysics(new G4DecayPhysics());
+
+  // EM physics
+  RegisterPhysics(new G4EmStandardPhysics());
+
+  // Radioactive decay
+  RegisterPhysics(new G4RadioactiveDecayPhysics());
+}
+
+PhysicsList::~PhysicsList() {}
+
+void PhysicsList::SetCuts() { G4VModularPhysicsList::SetCuts(); }
+
+} // namespace Geant
+} // namespace uLib

src/HEP/Geant/PhysicsList.hh

@@ -0,0 +1,20 @@
+#ifndef PhysicsList_h
+#define PhysicsList_h
+
+#include "G4VModularPhysicsList.hh"
+
+namespace uLib {
+namespace Geant {
+
+class PhysicsList : public G4VModularPhysicsList {
+public:
+  PhysicsList();
+  virtual ~PhysicsList();
+
+  virtual void SetCuts();
+};
+
+} // namespace Geant
+} // namespace uLib
+
+#endif

src/HEP/Geant/Scene.cpp

@@ -0,0 +1,195 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 <Geant4/G4Box.hh>
+#include <Geant4/G4LogicalVolume.hh>
+#include <Geant4/G4Material.hh>
+#include <Geant4/G4NistManager.hh>
+#include <Geant4/G4PVPlacement.hh>
+#include <Geant4/G4RunManager.hh>
+#include <Geant4/G4RunManagerFactory.hh>
+#include <Geant4/G4SystemOfUnits.hh>
+#include <Geant4/G4VPhysicalVolume.hh>
+
+#include "Core/Vector.h"
+#include "HEP/Geant/DetectorConstruction.hh"
+#include "Math/ContainerBox.h"
+#include "Math/Dense.h"
+#include "Solid.h"
+
+#include "Scene.h"
+#include "PhysicsList.hh"
+#include "ActionInitialization.hh"
+
+namespace uLib {
+namespace Geant {
+
+class SceneDetectorConstruction : public DetectorConstruction {
+public:
+  SceneDetectorConstruction(class SceneImpl *owner);
+  G4VPhysicalVolume *Construct() override;
+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(G4RunManagerType::Serial)),
+                m_Emitter(nullptr),
+                m_Output(nullptr) {
+    m_RunManager->SetUserInitialization(new PhysicsList);
+  }
+
+  // destructor //
+  ~SceneImpl() { 
+    if (m_RunManager) delete m_RunManager; 
+    if (m_World) delete m_World; 
+  }
+
+  void Initialize() {
+    // Set mandatory initialization classes for Geant4
+    m_RunManager->SetUserInitialization(new SceneDetectorConstruction(this));
+    m_RunManager->SetUserInitialization(
+        new ActionInitialization(m_Emitter, m_Output));
+
+    // Initialize Geant4
+    m_RunManager->Initialize();
+  }
+
+  // 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) 
+  : DetectorConstruction("Scene"), m_Owner(owner) {}
+
+G4VPhysicalVolume *SceneDetectorConstruction::Construct() {
+  return m_Owner->m_World->GetPhysical();
+}
+
+
+
+
+Scene::Scene() {
+  CheckGeant4Environment();
+  d = new SceneImpl();
+}
+Scene::~Scene() { delete d; }
+
+void Scene::AddSolid(Solid *solid, Solid *parent) { 
+  d->m_Solids.push_back(solid);
+  if (!d->m_World) {
+    d->m_World = solid;
+  } else {
+    solid->SetParent(parent ? parent : d->m_World);
+  }
+}
+
+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);
+  }
+
+  G4Box *solidWorld = new G4Box("World", 
+                                0.5 * size(0), 
+                                0.5 * size(1),
+                                0.5 * size(2));
+
+  G4LogicalVolume *logicWorld = new G4LogicalVolume(solidWorld,
+                                                    d->m_World->GetMaterial(),
+                                                    d->m_World->GetName());
+
+  d->m_World->SetLogical(logicWorld);
+  
+  G4PVPlacement *physWorld = new G4PVPlacement(
+    nullptr,
+    G4ThreeVector(0, 0, 0),
+    logicWorld,
+    d->m_World->GetName(),
+    0,
+    false,
+    0,
+    true);
+
+  d->m_World->SetPhysical(physWorld);
+
+  // 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

@@ -23,64 +23,52 @@
 
 //////////////////////////////////////////////////////////////////////////////*/
 
+#ifndef SCENE_H
+#define SCENE_H
 
 
-#ifndef SOLID_H
-#define SOLID_H
 
 #include "Core/Object.h"
-#include "Math/Dense.h"
+#include "Core/Vector.h"
-#include "Math/TriangleMesh.h"
+#include "Solid.h"
-#include "Detectors/Matter.h"
+#include "GeantEvent.h"
-
-class G4Material;
-class G4LogicalVolume;
-class G4TessellatedSolid;
 
+class G4VPhysicalVolume;
 
 namespace uLib {
+namespace Geant {
 
-class Solid : public Object {
+class EmitterPrimary;
+
+class Scene : public Object {
 public:
+  Scene();
+  ~Scene();
 
-    Solid();
+  void AddSolid(Solid *solid, Solid *parent = nullptr);
-    Solid(const char *name);
   
-    void SetNistMaterial(const char *name);
+  void ConstructWorldBox(const Vector3f &size, const char *material);
 
-    void SetMaterial(G4Material *material);
+  /// Get the world box
+  const Solid* GetWorld() const;
   
-    uLibGetMacro(Material,G4Material *)
+  ContainerBox* GetWorldBox() const;
-    uLibGetMacro(Logical,G4LogicalVolume *)
 
-protected:                
+  /// Set the primary generator (emitter) for the simulation.
-    G4Material         *m_Material;
+  /// The Scene does NOT take ownership of the emitter.
-    G4LogicalVolume    *m_Logical;
+  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);
 
-
-
-class TessellatedSolid : public Solid {
-    typedef Solid BaseClass;
-
-public:
-    TessellatedSolid(const char *name);
-
-    void SetMesh(TriangleMesh &mesh);
-
-    uLibGetMacro(Solid,G4TessellatedSolid *)
 private:
-    G4TessellatedSolid *m_Solid;
+  class SceneImpl *d;
 };
 
-
+} // namespace Geant
-
+} // namespace uLib
-
+#endif // SCENE_H
-}
-
-
-
-
-#endif // SOLID_H

src/HEP/Geant/Solid.cpp

@@ -0,0 +1,208 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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.
+
+//////////////////////////////////////////////////////////////////////////////*/
+
+// G4 Solid //
+#include <CLHEP/Units/SystemOfUnits.h>
+#include <Geant4/G4LogicalVolume.hh>
+#include <Geant4/G4Material.hh>
+#include <Geant4/G4NistManager.hh>
+
+// Tessellated solid //
+#include <Geant4/G4TessellatedSolid.hh>
+#include <Geant4/G4ThreeVector.hh>
+#include <Geant4/G4TriangularFacet.hh>
+#include <Geant4/G4Box.hh>
+#include <Geant4/G4PVPlacement.hh>
+
+
+#include "Math/Dense.h"
+#include "Math/Transform.h"
+
+#include "Solid.h"
+
+namespace uLib {
+namespace Geant {
+
+class DetectorsSolidImpl {
+public:
+  static G4ThreeVector getG4Vector3f(const Vector3f &vector) {
+    return G4ThreeVector(vector(0), vector(1), vector(2));
+  }
+
+  
+};
+
+Solid::Solid()
+    : 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_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();
+  G4Material *mat = nist->FindOrBuildMaterial(name);
+  if (mat) SetMaterial(mat);
+}
+
+void Solid::SetMaterial(G4Material *material) {
+  if (material) {
+    m_Material = material;
+    if (m_Logical) {
+      m_Logical->SetMaterial(material);
+    } else if (GetG4Solid()) {
+      m_Logical = new G4LogicalVolume(GetG4Solid(), m_Material, GetName());
+    }
+  }
+}
+
+void Solid::SetTransform(Matrix4f transform) {
+  uLib::AffineTransform t;
+  t.SetMatrix(transform);
+
+    // 2. Extract position and rotation for Geant4
+    Vector3f pos = t.GetPosition();
+    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();
+    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. If object is already placed, update its transformation
+    if (m_Physical) {
+        m_Physical->SetTranslation(*m_Position);
+        m_Physical->SetRotation(m_Rotation);
+    }
+}
+
+void Solid::SetParent(Solid *parent) {
+  if (!m_Logical) {
+    std::cerr << "logical volume not created for solid " << GetName() << std::endl;
+    return;
+  }
+
+  if(m_Physical) {
+    std::cerr << "physical volume already created for solid " << GetName() << std::endl;
+    return;
+  }
+
+  G4LogicalVolume* parentLogical = nullptr;
+  if (parent) {
+    parentLogical = parent->GetLogical();
+    if (!parentLogical) {
+      std::cerr << "parent logical volume not created for solid " << parent->GetName() << std::endl;
+      return;
+    }
+  }
+
+  // G4PVPlacement
+  m_Physical = new G4PVPlacement(
+      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)
+      false,                 // Boolean operations (usually false)
+      0,                     // Copy number
+      true                   // Check overlaps (useful to enable in debug phase)
+  );
+}
+
+
+
+
+
+
+
+TessellatedSolid::TessellatedSolid(const char *name)
+    : BaseClass(name), m_Solid(new G4TessellatedSolid(name)) {
+    }
+
+void TessellatedSolid::SetMesh(TriangleMesh &mesh) {
+  G4TessellatedSolid *ts = this->m_Solid;
+  for (int i = 0; i < mesh.Triangles().size(); ++i) {
+    const Vector3i &trg = mesh.Triangles().at(i);
+    G4TriangularFacet *facet = new G4TriangularFacet(
+        DetectorsSolidImpl::getG4Vector3f(mesh.Points().at(trg(0))),
+        DetectorsSolidImpl::getG4Vector3f(mesh.Points().at(trg(1))),
+        DetectorsSolidImpl::getG4Vector3f(mesh.Points().at(trg(2))), ABSOLUTE);
+    ts->AddFacet((G4VFacet *)facet);
+  }
+  if (this->m_Logical) {
+    this->m_Logical->SetSolid(ts);
+  }
+}
+
+
+
+
+
+BoxSolid::BoxSolid(const char *name, ContainerBox *box) : BaseClass(name) {
+  m_Solid = new G4Box(name, 1,1,1);
+  m_Object = box;
+  Object::connect(box, &ContainerBox::Updated, this, &BoxSolid::Update);
+  if (m_Logical) {
+    m_Logical->SetSolid(m_Solid);
+  }
+  Update();
+}
+
+void BoxSolid::Update() {
+  if (m_Object) {
+    Vector3f size = m_Object->GetSize();
+    m_Solid->SetXHalfLength(size(0) * 0.5);
+    m_Solid->SetYHalfLength(size(1) * 0.5);
+    m_Solid->SetZHalfLength(size(2) * 0.5);
+    
+    // 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());
+  }
+}
+
+
+} // namespace Geant
+} // namespace uLib

src/HEP/Geant/Solid.h

@@ -0,0 +1,121 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 SOLID_H
+#define SOLID_H
+
+#include "Core/Object.h"
+#include "Geant/Matter.h"
+#include <Geant4/G4LogicalVolume.hh>
+#include "Math/ContainerBox.h"
+#include "Math/Dense.h"
+#include "Math/TriangleMesh.h"
+
+class G4Material;
+class G4LogicalVolume;
+class G4TessellatedSolid;
+class G4Box;
+
+namespace uLib {
+namespace Geant {
+
+class Solid : public Object {
+public:
+  Solid();
+  Solid(const char *name);
+  virtual ~Solid();
+
+  void SetNistMaterial(const char *name);
+  void SetMaterial(G4Material *material);
+  void SetSizeUnit(const char *unit);
+  // Implementiamo SetParent qui, per tutti.
+  virtual void SetParent(Solid *parent); 
+
+  // Setters per la posizione (necessari per il piazzamento)
+  void SetTransform(Matrix4f transform); 
+
+  uLibGetMacro(Material, G4Material *) 
+  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(); 
+  }
+
+protected:
+
+  std::string m_Name;
+  G4Material *m_Material;
+  G4LogicalVolume *m_Logical;
+  G4VPhysicalVolume *m_Physical; // <-- Memorizza l'oggetto posizionato
+  G4ThreeVector *m_Position;      // <-- Offset rispetto al centro del padre
+  G4RotationMatrix* m_Rotation;  // <-- Rotazione rispetto al padre
+};
+
+
+
+
+
+class TessellatedSolid : public Solid {
+  typedef Solid BaseClass;
+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();
+
+private : 
+  G4TessellatedSolid *m_Solid;
+};
+
+
+
+
+
+
+class BoxSolid : public Solid {
+  typedef Solid BaseClass;
+
+public:
+  BoxSolid(const char *name, ContainerBox *box);
+  virtual G4VSolid* GetG4Solid() const override { return (G4VSolid*)m_Solid; }
+
+public slots:
+  void Update();
+
+private:
+  ContainerBox *m_Object;
+  G4Box *m_Solid;
+};
+
+} // namespace Geant
+} // namespace uLib
+
+#endif // SOLID_H

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/ActionInitialization.cpp

@@ -0,0 +1,23 @@
+#include "HEP/Geant/ActionInitialization.hh" // Il file appena creato
+#include "G4RunManagerFactory.hh"            // Per il RunManager moderno
+// ... altri include (DetectorConstruction, PhysicsList, ecc.)
+
+int main(int argc, char **argv) {
+  // Creazione del Run Manager
+  auto *runManager = G4RunManagerFactory::CreateRunManager();
+
+  // 1. Inizializzazione della Geometria
+  // runManager->SetUserInitialization(new DetectorConstruction());
+
+  // 2. Inizializzazione della Fisica
+  // runManager->SetUserInitialization(new PhysicsList());
+
+  // 3. INIZIALIZZAZIONE DELLE AZIONI (Il nostro generatore!)
+  runManager->SetUserInitialization(new uLib::Geant::ActionInitialization());
+
+  // ... Inizializzazione del kernel ( runManager->Initialize(); ), UI manager,
+  // vis manager, ecc.
+
+  delete runManager;
+  return 0;
+}

src/HEP/Geant/testing/CMakeLists.txt

@@ -0,0 +1,15 @@
+# TESTS
+set(TESTS
+       SolidTest
+       EventTest
+       GeantApp
+       ActionInitialization
+)
+
+set(LIBRARIES
+       ${PACKAGE_LIBPREFIX}Core
+       ${PACKAGE_LIBPREFIX}Math
+       ${PACKAGE_LIBPREFIX}Geant
+       Eigen3::Eigen
+)
+uLib_add_tests(Geant)

src/HEP/Geant/testing/EventTest.cpp

@@ -0,0 +1,90 @@
+#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>
+#include <Geant4/G4LogicalVolume.hh>
+#include <Geant4/G4TessellatedSolid.hh>
+#include <string.h>
+
+using namespace uLib;
+
+int main() {
+    BEGIN_TESTING(Geant Event);
+
+    // Test: Scene with iron cube in air, launch muons, collect events //
+    {
+        // 1. Create world box (air, 30m x 30m x 30m)
+        Geant::Scene scene;
+        scene.ConstructWorldBox(Vector3f(30_m, 30_m, 30_m), "G4_AIR");
+
+        // 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);
+
+        // 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);
+
+        // 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

@@ -0,0 +1,16 @@
+
+
+#include "Math/ContainerBox.h"
+#include "Math/Dense.h"
+#include "HEP/Geant/Scene.h"
+
+using namespace uLib;
+
+int main() {
+
+  uLib::Geant::Scene scene;
+  scene.ConstructWorldBox(Vector3f(100, 100, 100), "G4_AIR");
+  scene.Initialize();
+
+  return 0;
+}

src/HEP/Geant/testing/SolidTest.cpp

@@ -0,0 +1,73 @@
+#include "Geant/Solid.h"
+#include "Math/TriangleMesh.h"
+#include "testing-prototype.h"
+#include <Geant4/G4Material.hh>
+#include <Geant4/G4NistManager.hh>
+#include <Geant4/G4LogicalVolume.hh>
+#include <Geant4/G4TessellatedSolid.hh>
+#include <string.h>
+
+using namespace uLib;
+
+int main() {
+    BEGIN_TESTING(Geant Solid);
+
+    // Test Solid initialization and NIST material //
+    {
+        Geant::Solid solid("test_solid");
+        // 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");
+    }
+
+    // Test TessellatedSolid with a simple mesh //
+    {
+        Geant::TessellatedSolid tsolid("test_tessellated");
+        tsolid.SetNistMaterial("G4_AIR");
+        TEST1(tsolid.GetLogical() != nullptr);
+        TEST1(tsolid.GetSolid() != nullptr);
+
+        // 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));
+
+
+
+        tsolid.SetMesh(mesh);
+        TEST1(tsolid.GetSolid()->GetNumberOfFacets() == 12);
+    }
+
+    END_TESTING
+}

src/HEP/Geant/testing/testing-prototype.h

@@ -25,28 +25,13 @@
 
 
 
-#ifndef SCENE_H
+#include <stdio.h>
-#define SCENE_H
 
-#include "Core/Object.h"
+#define BEGIN_TESTING(name)                \
-#include "Core/Vector.h"
+static int _fail = 0;                      \
-#include "Solid.h"
+printf("..:: Testing " #name " ::..\n");
 
-namespace uLib {
+#define TEST1(val) _fail += (val)==0
+#define TEST0(val) _fail += (val)!=0
+#define END_TESTING return _fail;
 
-
-class DetectorsScene : public Object {
-public:
-    DetectorsScene();
-    ~DetectorsScene();
-
-
-    void AddSolid(const Solid &solid);
-
-private:
-    class DetectorsScenePimpl *d;
-};
-
-
-}
-#endif // SCENE_H

src/HEP/MuonTomography/Hit.h

@@ -39,12 +39,6 @@ class Hit {
     Type m_DriftTime;
 };
 
-
-
-
-
-
-
 class HitMC {
 public:
     virtual Id_t     GetChamber() = 0;

src/Math/CLHEP_defs.hh

@@ -0,0 +1,373 @@
+#if !defined(HAVE_GEANT4) && !defined(HEP_SYSTEM_OF_UNITS_H) && !defined(HEP_PHYSICAL_CONSTANTS_H)
+#ifndef CLHEP_defs_h
+#define CLHEP_defs_h
+
+
+
+namespace CLHEP {
+
+
+// UNITS //
+//
+//
+//
+static constexpr double     pi  = 3.14159265358979323846;
+static constexpr double  twopi  = 2*pi;
+static constexpr double halfpi  = pi/2;
+static constexpr double     pi2 = pi*pi;
+
+//
+//
+//
+static constexpr double millimeter  = 1.;                        
+static constexpr double millimeter2 = millimeter*millimeter;
+static constexpr double millimeter3 = millimeter*millimeter*millimeter;
+
+static constexpr double centimeter  = 10.*millimeter;   
+static constexpr double centimeter2 = centimeter*centimeter;
+static constexpr double centimeter3 = centimeter*centimeter*centimeter;
+
+static constexpr double meter  = 1000.*millimeter;                  
+static constexpr double meter2 = meter*meter;
+static constexpr double meter3 = meter*meter*meter;
+
+static constexpr double kilometer = 1000.*meter;                   
+static constexpr double kilometer2 = kilometer*kilometer;
+static constexpr double kilometer3 = kilometer*kilometer*kilometer;
+
+static constexpr double parsec = 3.0856775807e+16*meter;
+
+static constexpr double micrometer = 1.e-6 *meter;             
+static constexpr double  nanometer = 1.e-9 *meter;
+static constexpr double  angstrom  = 1.e-10*meter;
+static constexpr double  fermi     = 1.e-15*meter;
+
+static constexpr double      barn = 1.e-28*meter2;
+static constexpr double millibarn = 1.e-3 *barn;
+static constexpr double microbarn = 1.e-6 *barn;
+static constexpr double  nanobarn = 1.e-9 *barn;
+static constexpr double  picobarn = 1.e-12*barn;
+
+// symbols
+static constexpr double nm  = nanometer;                        
+static constexpr double um  = micrometer;                        
+
+static constexpr double mm  = millimeter;                        
+static constexpr double mm2 = millimeter2;
+static constexpr double mm3 = millimeter3;
+
+static constexpr double cm  = centimeter;   
+static constexpr double cm2 = centimeter2;
+static constexpr double cm3 = centimeter3;
+
+static constexpr double liter = 1.e+3*cm3;
+static constexpr double  L = liter;
+static constexpr double dL = 1.e-1*liter;
+static constexpr double cL = 1.e-2*liter;
+static constexpr double mL = 1.e-3*liter;       
+
+static constexpr double m  = meter;                  
+static constexpr double m2 = meter2;
+static constexpr double m3 = meter3;
+
+static constexpr double km  = kilometer;                   
+static constexpr double km2 = kilometer2;
+static constexpr double km3 = kilometer3;
+
+static constexpr double pc = parsec;
+
+//
+// Angle
+//
+static constexpr double radian      = 1.;                  
+static constexpr double milliradian = 1.e-3*radian;
+static constexpr double degree = (pi/180.0)*radian;
+
+static constexpr double   steradian = 1.;
+  
+// symbols
+static constexpr double rad  = radian;
+static constexpr double mrad = milliradian;
+static constexpr double sr   = steradian;
+static constexpr double deg  = degree;
+
+//
+// Time [T]
+//
+static constexpr double nanosecond  = 1.;
+static constexpr double second      = 1.e+9 *nanosecond;
+static constexpr double millisecond = 1.e-3 *second;
+static constexpr double microsecond = 1.e-6 *second;
+static constexpr double  picosecond = 1.e-12*second;
+
+static constexpr double minute = 60*second;
+static constexpr double hour   = 60*minute;
+static constexpr double day    = 24*hour;
+static constexpr double year   = 365*day;  
+
+static constexpr double hertz = 1./second;
+static constexpr double kilohertz = 1.e+3*hertz;
+static constexpr double megahertz = 1.e+6*hertz;
+
+// symbols
+static constexpr double ns = nanosecond;
+static constexpr double  s = second;
+static constexpr double ms = millisecond;
+static constexpr double us = microsecond;
+static constexpr double ps = picosecond;
+
+//
+// Electric charge [Q]
+//
+static constexpr double eplus = 1. ;// positron charge
+static constexpr double e_SI  = 1.602176634e-19;// positron charge in coulomb
+static constexpr double coulomb = eplus/e_SI;// coulomb = 6.24150 e+18 * eplus
+
+//
+// Energy [E]
+//
+static constexpr double megaelectronvolt = 1. ;
+static constexpr double     electronvolt = 1.e-6*megaelectronvolt;
+static constexpr double kiloelectronvolt = 1.e-3*megaelectronvolt;
+static constexpr double gigaelectronvolt = 1.e+3*megaelectronvolt;
+static constexpr double teraelectronvolt = 1.e+6*megaelectronvolt;
+static constexpr double petaelectronvolt = 1.e+9*megaelectronvolt;
+static constexpr double millielectronvolt = 1.e-9*megaelectronvolt;  
+
+static constexpr double joule = electronvolt/e_SI;// joule = 6.24150 e+12 * MeV
+
+// symbols
+static constexpr double MeV = megaelectronvolt;
+static constexpr double  eV = electronvolt;
+static constexpr double keV = kiloelectronvolt;
+static constexpr double GeV = gigaelectronvolt;
+static constexpr double TeV = teraelectronvolt;
+static constexpr double PeV = petaelectronvolt;
+
+//
+// Mass [E][T^2][L^-2]
+//
+static constexpr double  kilogram = joule*second*second/(meter*meter);   
+static constexpr double      gram = 1.e-3*kilogram;
+static constexpr double milligram = 1.e-3*gram;
+
+// symbols
+static constexpr double  kg = kilogram;
+static constexpr double   g = gram;
+static constexpr double  mg = milligram;
+
+//
+// Power [E][T^-1]
+//
+static constexpr double watt = joule/second;// watt = 6.24150 e+3 * MeV/ns
+
+//
+// Force [E][L^-1]
+//
+static constexpr double newton = joule/meter;// newton = 6.24150 e+9 * MeV/mm
+
+//
+// Pressure [E][L^-3]
+//
+#define pascal hep_pascal                          // a trick to avoid warnings 
+static constexpr double hep_pascal = newton/m2;   // pascal = 6.24150 e+3 * MeV/mm3
+static constexpr double bar        = 100000*pascal; // bar    = 6.24150 e+8 * MeV/mm3
+static constexpr double atmosphere = 101325*pascal; // atm    = 6.32420 e+8 * MeV/mm3
+
+//
+// Electric current [Q][T^-1]
+//
+static constexpr double      ampere = coulomb/second; // ampere = 6.24150 e+9 * eplus/ns
+static constexpr double milliampere = 1.e-3*ampere;
+static constexpr double microampere = 1.e-6*ampere;
+static constexpr double  nanoampere = 1.e-9*ampere;
+
+//
+// Electric potential [E][Q^-1]
+//
+static constexpr double megavolt = megaelectronvolt/eplus;
+static constexpr double kilovolt = 1.e-3*megavolt;
+static constexpr double     volt = 1.e-6*megavolt;
+
+//
+// Electric resistance [E][T][Q^-2]
+//
+static constexpr double ohm = volt/ampere;// ohm = 1.60217e-16*(MeV/eplus)/(eplus/ns)
+
+//
+// Electric capacitance [Q^2][E^-1]
+//
+static constexpr double farad = coulomb/volt;// farad = 6.24150e+24 * eplus/Megavolt
+static constexpr double millifarad = 1.e-3*farad;
+static constexpr double microfarad = 1.e-6*farad;
+static constexpr double  nanofarad = 1.e-9*farad;
+static constexpr double  picofarad = 1.e-12*farad;
+
+//
+// Magnetic Flux [T][E][Q^-1]
+//
+static constexpr double weber = volt*second;// weber = 1000*megavolt*ns
+
+//
+// Magnetic Field [T][E][Q^-1][L^-2]
+//
+static constexpr double tesla     = volt*second/meter2;// tesla =0.001*megavolt*ns/mm2
+
+static constexpr double gauss     = 1.e-4*tesla;
+static constexpr double kilogauss = 1.e-1*tesla;
+
+//
+// Inductance [T^2][E][Q^-2]
+//
+static constexpr double henry = weber/ampere;// henry = 1.60217e-7*MeV*(ns/eplus)**2
+
+//
+// Temperature
+//
+static constexpr double kelvin = 1.;
+
+//
+// Amount of substance
+//
+static constexpr double mole = 1.;
+
+//
+// Activity [T^-1]
+//
+static constexpr double becquerel = 1./second ;
+static constexpr double curie = 3.7e+10 * becquerel;
+static constexpr double kilobecquerel = 1.e+3*becquerel;
+static constexpr double megabecquerel = 1.e+6*becquerel;
+static constexpr double gigabecquerel = 1.e+9*becquerel;
+static constexpr double millicurie = 1.e-3*curie;
+static constexpr double microcurie = 1.e-6*curie;
+static constexpr double Bq = becquerel;
+static constexpr double kBq = kilobecquerel;
+static constexpr double MBq = megabecquerel;
+static constexpr double GBq = gigabecquerel;
+static constexpr double Ci = curie;
+static constexpr double mCi = millicurie;
+static constexpr double uCi = microcurie;
+
+//
+// Absorbed dose [L^2][T^-2]
+//
+static constexpr double      gray = joule/kilogram ;
+static constexpr double  kilogray = 1.e+3*gray;
+static constexpr double milligray = 1.e-3*gray;
+static constexpr double microgray = 1.e-6*gray;
+
+//
+// Luminous intensity [I]
+//
+static constexpr double candela = 1.;
+
+//
+// Luminous flux [I]
+//
+static constexpr double lumen = candela*steradian;
+
+//
+// Illuminance [I][L^-2]
+//
+static constexpr double lux = lumen/meter2;
+
+//
+// Miscellaneous
+//
+static constexpr double perCent     = 0.01 ;
+static constexpr double perThousand = 0.001;
+static constexpr double perMillion  = 0.000001;
+
+
+
+// CONSTANTS //
+
+//
+// 
+//
+static constexpr double Avogadro = 6.02214076e+23/mole;
+
+//
+// c   = 299.792458 mm/ns
+// c^2 = 898.7404 (mm/ns)^2 
+//
+static constexpr double c_light   = 2.99792458e+8 * m/s;
+static constexpr double c_squared = c_light * c_light;
+
+//
+// h     = 4.13566e-12 MeV*ns
+// hbar  = 6.58212e-13 MeV*ns
+// hbarc = 197.32705e-12 MeV*mm
+//
+static constexpr double h_Planck      = 6.62607015e-34 * joule*s;
+static constexpr double hbar_Planck   = h_Planck/twopi;
+static constexpr double hbarc         = hbar_Planck * c_light;
+static constexpr double hbarc_squared = hbarc * hbarc;
+
+//
+//
+//
+static constexpr double electron_charge = - eplus; // see SystemOfUnits.h
+static constexpr double e_squared = eplus * eplus;
+
+//
+// amu_c2 - atomic equivalent mass unit
+//        - AKA, unified atomic mass unit (u)
+// amu    - atomic mass unit
+//
+static constexpr double electron_mass_c2 = 0.510998910 * MeV;
+static constexpr double   proton_mass_c2 = 938.272013 * MeV;
+static constexpr double  neutron_mass_c2 = 939.56536 * MeV;
+static constexpr double           amu_c2 = 931.494028 * MeV;
+static constexpr double              amu = amu_c2/c_squared;
+
+//
+// permeability of free space mu0    = 2.01334e-16 Mev*(ns*eplus)^2/mm
+// permittivity of free space epsil0 = 5.52636e+10 eplus^2/(MeV*mm)
+//
+static constexpr double mu0      = 4*pi*1.e-7 * henry/m;
+static constexpr double epsilon0 = 1./(c_squared*mu0);
+
+//
+// electromagnetic coupling = 1.43996e-12 MeV*mm/(eplus^2)
+//
+static constexpr double elm_coupling           = e_squared/(4*pi*epsilon0);
+static constexpr double fine_structure_const   = elm_coupling/hbarc;
+static constexpr double classic_electr_radius  = elm_coupling/electron_mass_c2;
+static constexpr double electron_Compton_length = hbarc/electron_mass_c2;
+static constexpr double Bohr_radius = electron_Compton_length/fine_structure_const;
+
+static constexpr double alpha_rcl2 = fine_structure_const
+                                   *classic_electr_radius
+                                   *classic_electr_radius;
+
+static constexpr double twopi_mc2_rcl2 = twopi*electron_mass_c2
+                                             *classic_electr_radius
+                                             *classic_electr_radius;
+
+static constexpr double Bohr_magneton = (eplus*hbarc*c_light)/(2*electron_mass_c2);
+static constexpr double nuclear_magneton = (eplus*hbarc*c_light)/(2*proton_mass_c2);
+
+//
+//
+//
+static constexpr double k_Boltzmann = 8.617333e-11 * MeV/kelvin;
+
+//
+//
+//
+static constexpr double STP_Temperature = 273.15*kelvin;
+static constexpr double STP_Pressure    = 1.*atmosphere;
+static constexpr double kGasThreshold   = 10.*mg/cm3;
+
+//
+//
+//
+static constexpr double universe_mean_density = 1.e-25*g/cm3;
+
+}  // namespace CLHEP
+
+#endif // CLHEP_defs_h
+#endif // HAVE_GEANT4 / CLHEP checks
+

src/Math/CMakeLists.txt

@@ -19,20 +19,25 @@ set(HEADERS ContainerBox.h
             VoxImageFilterCustom.hpp
             Accumulator.h
             TriangleMesh.h
+            QuadMesh.h
             BitCode.h
             Structured2DGrid.h
-            Structured4DGrid.h)
+            Structured4DGrid.h
+            Units.h
+            CLHEP_defs.hh)
 
 set(SOURCES VoxRaytracer.cpp
             StructuredData.cpp
             StructuredGrid.cpp
             VoxImage.cpp
             TriangleMesh.cpp
+            QuadMesh.cpp
             Dense.cpp
             Structured2DGrid.cpp
             Structured4DGrid.cpp)
 
-set(LIBRARIES Eigen3::Eigen
+set(LIBRARIES ${PACKAGE_LIBPREFIX}Core
+              Eigen3::Eigen
               ${ROOT_LIBRARIES}
               ${VTK_LIBRARIES})
 
@@ -54,7 +59,7 @@ endif()
 
 
 install(TARGETS ${libname}
-        EXPORT "${PROJECT_NAME}Targets"
+        EXPORT "uLibTargets"
         RUNTIME DESTINATION ${INSTALL_BIN_DIR} COMPONENT bin
         LIBRARY DESTINATION ${INSTALL_LIB_DIR} COMPONENT lib)
 

src/Math/ContainerBox.h

@@ -23,76 +23,169 @@
 
 //////////////////////////////////////////////////////////////////////////////*/
 
-
-
 #ifndef U_CONTAINERBOX_H
 #define U_CONTAINERBOX_H
 
 #include "Geometry.h"
-
+#include "Core/Object.h"
+#include "Math/Dense.h"
+#include "Math/Transform.h"
+#include <utility>
 
 namespace uLib {
 
+/**
+ * @brief Represents an oriented bounding box (OBB) within a hierarchical
+ * transformation system.
+ *
+ * ContainerBox inherits from AffineTransform, which defines its parent
+ * coordinate system. It contains an internal local transformation (m_LocalT)
+ * that defines the box's specific origin and size relative to its own
+ * coordinate system.
+ */
+class ContainerBox : public AffineTransform, public Object {
 
+  typedef AffineTransform BaseClass;
 
-class ContainerBox : public AffineTransform {
 public:
-    ContainerBox() : m_LocalT(this) {}
+  /**
+   * @brief Default constructor.
+   * Initializes the local transformation with this instance as its parent.
+   */
+  ContainerBox()
+      : m_LocalT(this) // BaseClass is Parent of m_LocalTransform
+  {}
 
-    ContainerBox(const ContainerBox &copy) :
+  /**
-        m_LocalT(this),
+   * @brief Constructor with size.
-        AffineTransform(copy)
+   * @param size The size vector.
-    {
+   */
-        // FIX for performance //
+  ContainerBox(const Vector3f &size) : m_LocalT(this) { this->SetSize(size); }
-        this->SetOrigin(copy.GetOrigin());
-        this->SetSize(copy.GetSize());
-    }
 
-    inline void SetOrigin(const Vector3f &v) { m_LocalT.SetPosition(v); }
+  /**
+   * @brief Copy constructor.
+   * @param copy The ContainerBox instance to copy from.
+   */
+  ContainerBox(const ContainerBox &copy)
+      : m_LocalT(this), // BaseClass is Parent of m_LocalTransform
+        AffineTransform(copy) {
+    this->SetOrigin(copy.GetOrigin());
+    this->SetSize(copy.GetSize());
+  }
 
-    inline Vector3f GetOrigin() const { return m_LocalT.GetPosition(); }
+  /**
+   * @brief Sets the box origin relative to its coordinate system.
+   * @param v The origin position vector.
+   */
+  inline void SetOrigin(const Vector3f &v) { m_LocalT.SetPosition(v); }
 
-    void SetSize(const Vector3f &v) {
+  /**
-        Vector3f pos = this->GetOrigin();
+   * @brief Gets the box origin relative to its coordinate system.
-        m_LocalT = AffineTransform(this);
+   * @return The origin position vector.
-        m_LocalT.Scale(v);
+   */
-        m_LocalT.SetPosition(pos);
+  inline Vector3f GetOrigin() const { return m_LocalT.GetPosition(); }
-    }
 
-    inline Vector3f GetSize() const { return m_LocalT.GetScale(); }
+  /**
+   * @brief Sets the size of the box.
+   * Re-initializes the local transformation and applies the new scale.
+   * @param v The size vector (width, height, depth).
+   */
+  void SetSize(const Vector3f &v) {
+    Vector3f pos = this->GetOrigin();
+    m_LocalT = AffineTransform(this); // regenerate local transform
+    m_LocalT.Scale(v);
+    m_LocalT.SetPosition(pos);
+  }
 
-    // FIX... //
+  /**
-    inline void FlipLocalAxes(int first, int second)
+   * @brief Gets the current size (scale) of the box.
-    { m_LocalT.FlipAxes(first,second); }
+   * @return The size vector.
+   */
+  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));
+  }
 
-    Matrix4f GetWorldMatrix() const { return m_LocalT.GetWorldMatrix(); }
+  /**
+   * @brief Swaps two local axes of the box.
+   * @param first Index of the first axis (0=X, 1=Y, 2=Z).
+   * @param second Index of the second axis (0=X, 1=Y, 2=Z).
+   */
+  inline void FlipLocalAxes(int first, int second) {
+    m_LocalT.FlipAxes(first, second);
+  }
 
-    inline Vector4f GetWorldPoint(const Vector4f &v) const {
+  /**
-        return m_LocalT.GetWorldMatrix() * v;
+   * @brief Returns the world transformation matrix of the box's volume.
-    }
+   * @return A 4x4 transformation matrix.
+   */
+  Matrix4f GetWorldMatrix() const { return m_LocalT.GetWorldMatrix(); }
 
-    inline Vector4f GetWorldPoint(const float x, const float y, const float z) {
+  /**
-        return this->GetWorldPoint(Vector4f(x,y,z,1));
+   * @brief Returns the local transformation matrix of the box's volume.
-    }
+   * @return A 4x4 transformation matrix.
+   */
+  Matrix4f GetLocalMatrix() const { return m_LocalT.GetMatrix(); }
 
-    inline Vector4f GetLocalPoint(const Vector4f &v) const {
+  /**
-        return m_LocalT.GetWorldMatrix().inverse() * v;
+   * @brief Transforms a point from box-local space to world space.
-    }
+   * @param v The local point (4D homogeneous vector).
+   * @return The transformed point in world space.
+   */
+  inline Vector4f GetWorldPoint(const Vector4f &v) const {
+    return m_LocalT.GetWorldMatrix() * v;
+  }
 
-    inline Vector4f GetLocalPoint(const float x, const float y, const float z) {
+  /**
-         return this->GetLocalPoint(Vector4f(x,y,z,1));
+   * @brief Transforms a point from box-local space coordinates to world space.
-    }
+   * @param x X coordinate in local space.
+   * @param y Y coordinate in local space.
+   * @param z Z coordinate in local space.
+   * @return The transformed point in world space.
+   */
+  inline Vector4f GetWorldPoint(const float x, const float y, const float z) {
+    return this->GetWorldPoint(Vector4f(x, y, z, 1));
+  }
 
+  /**
+   * @brief Transforms a point from world space to box-local space.
+   * @param v The world point (4D homogeneous vector).
+   * @return The transformed point in box-local space.
+   */
+  inline Vector4f GetLocalPoint(const Vector4f &v) const {
+    return m_LocalT.GetWorldMatrix().inverse() * v;
+  }
 
-protected:
+  /**
+   * @brief Transforms a point from world space coordinates to box-local space.
+   * @param x X coordinate in world space.
+   * @param y Y coordinate in world space.
+   * @param z Z coordinate in world space.
+   * @return The transformed point in box-local space.
+   */
+  inline Vector4f GetLocalPoint(const float x, const float y, const float z) {
+    return this->GetLocalPoint(Vector4f(x, y, z, 1));
+  }
+
+  /** Translate using transformation chain */
+  using BaseClass::Translate;
+
+  /** Rotate using transformation chain */
+  using BaseClass::Rotate;
+
+  /** Scale using transformation chain */
+  using BaseClass::Scale;
+
+signals:
+
+  // signal to emit when the box is updated //
+  void Updated() { ULIB_SIGNAL_EMIT(ContainerBox::Updated); }
 
 private:
-    AffineTransform  m_LocalT;
+  AffineTransform m_LocalT;
 };
 
-
+} // namespace uLib
-}
-
-
 
 #endif // CONTAINERBOX_H

src/Math/Geometry.h

@@ -36,7 +36,7 @@ namespace uLib {
 class Geometry : public AffineTransform {
 public:
 
-    inline Vector4f GetWorldPoint(const Vector4f &v) const {
+    inline Vector4f GetWorldPoint(const Vector4f v) const {
         return this->GetWorldMatrix() * v;
     }
 
@@ -44,7 +44,7 @@ public:
         return this->GetWorldPoint(Vector4f(x,y,z,1));
     }
 
-    inline Vector4f GetLocalPoint(const Vector4f &v) const {
+    inline Vector4f GetLocalPoint(const Vector4f v) const {
         return this->GetWorldMatrix().inverse() * v;
     }
 

src/Math/QuadMesh.cpp

@@ -0,0 +1,77 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 "QuadMesh.h"
+
+
+namespace uLib {
+
+void QuadMesh::PrintSelf(std::ostream &o)
+{
+    o << " // ------- QUAD MESH ------- // \n" ;
+    o << "    #Points : " << m_Points.size() << "\n";
+    o << "    #Quads  : " << m_Quads.size() << "\n";
+    for(int i=0; i < m_Quads.size(); ++i ) {
+        o << "     - quad[" << i << "]" <<
+             " " << m_Quads[i](0) << "->(" << m_Points[m_Quads[i](0)].transpose() << ") " <<
+             " " << m_Quads[i](1) << "->(" << m_Points[m_Quads[i](1)].transpose() << ") " <<
+             " " << m_Quads[i](2) << "->(" << m_Points[m_Quads[i](2)].transpose() << ") " <<
+             " " << m_Quads[i](3) << "->(" << m_Points[m_Quads[i](3)].transpose() << ") " <<
+             " \n";
+    }
+    o << " // ------------------------- // \n";
+}
+
+void QuadMesh::AddPoint(const Vector3f &pt)
+{
+    this->m_Points.push_back(pt);
+}
+
+void QuadMesh::AddQuad(const Id_t *id)
+{
+    Vector4i quad(id[0],id[1],id[2],id[3]);
+    this->m_Quads.push_back(quad);
+}
+
+void QuadMesh::AddQuad(const Vector4i &id)
+{
+    this->m_Quads.push_back(id);
+}
+
+Vector3f QuadMesh::GetNormal(const Id_t id) const
+{
+    const Vector4i &quad = m_Quads.at(id);
+    const Vector3f &v0 = m_Points.at(quad(0));
+    const Vector3f &v1 = m_Points.at(quad(1));
+    const Vector3f &v3 = m_Points.at(quad(3));
+    
+    Vector3f edge1 = v1 - v0;
+    Vector3f edge2 = v3 - v0;
+    
+    return edge1.cross(edge2).normalized();
+}
+
+
+}

src/Math/QuadMesh.h

@@ -0,0 +1,60 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 QUADMESH_H
+#define QUADMESH_H
+
+#include <vector>
+
+#include "Math/Dense.h"
+
+namespace uLib {
+
+class QuadMesh
+{
+public:
+    void PrintSelf(std::ostream &o);
+
+    void AddPoint(const Vector3f &pt);
+
+    void AddQuad(const Id_t *id);
+    void AddQuad(const Vector4i &id);
+
+    inline std::vector<Vector3f> & Points() { return this->m_Points; }
+    inline std::vector<Vector4i> & Quads() { return this->m_Quads; }
+
+    const Vector4i & GetQuad(const Id_t id) const { return m_Quads.at(id); }
+    Vector3f GetNormal(const Id_t id) const;
+
+private:
+    std::vector<Vector3f> m_Points;
+    std::vector<Vector4i> m_Quads;
+
+};
+
+}
+
+
+#endif // QUADMESH_H

src/Math/Transform.h

@@ -84,8 +84,8 @@ public:
 
     inline void SetParent(AffineTransform *name) { this->m_Parent = name; }
 
-    inline void SetMatrix (Matrix4f &mat) { m_T.matrix() = mat; }
+    inline void SetMatrix (Matrix4f mat) { m_T.matrix() = mat; }
-    inline Matrix4f& GetMatrix () { return m_T.matrix(); }
+    inline Matrix4f GetMatrix() const { return m_T.matrix(); }
 
     Matrix4f GetWorldMatrix() const
     {
@@ -93,22 +93,26 @@ public:
         else return  m_Parent->GetWorldMatrix() * m_T.matrix(); // T = B * A //
     }
 
-    inline void SetPosition(const Vector3f &v) { this->m_T.translation() = v; }
+    inline void SetPosition(const Vector3f v) { this->m_T.translation() = v; }
 
     inline Vector3f GetPosition() const { return this->m_T.translation(); }
 
-    inline void SetRotation(const Matrix3f &m) { this->m_T.linear() = m; }
+    inline void SetRotation(const Matrix3f m) { this->m_T.linear() = m; }
 
     inline Matrix3f GetRotation() const { return this->m_T.rotation(); }
 
-    inline void Translate(const Vector3f &v) { this->m_T.translate(v); }
+    inline void Translate(const Vector3f v) { this->m_T.translate(v); }
 
-    inline void Scale(const Vector3f &v) { this->m_T.scale(v); }
+    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); }
+    inline void Rotate(const Matrix3f m) { this->m_T.rotate(m); }
 
     inline void Rotate(const float angle, Vector3f axis)
     {
@@ -122,12 +126,12 @@ public:
         Rotate(angle,euler_axis);
     }
 
-    inline void PreRotate(const Matrix3f &m) { this->m_T.prerotate(m); }
+    inline void PreRotate(const Matrix3f m) { this->m_T.prerotate(m); }
 
-    inline void QuaternionRotate(const Vector4f &q)
+    inline void QuaternionRotate(const Vector4f q)
     { this->m_T.rotate(Eigen::Quaternion<float>(q)); }
 
-    inline void EulerYZYRotate(const Vector3f &e) {
+    inline void EulerYZYRotate(const Vector3f e) {
         Matrix3f mat;
         mat =     Eigen::AngleAxisf(e.x(), Vector3f::UnitY())
                 * Eigen::AngleAxisf(e.y(), Vector3f::UnitZ())

src/Math/TriangleMesh.cpp

@@ -24,8 +24,6 @@
 //////////////////////////////////////////////////////////////////////////////*/
 
 
-
-
 #include "TriangleMesh.h"
 
 
@@ -65,5 +63,18 @@ void TriangleMesh::AddTriangle(const Vector3i &id)
     this->m_Triangles.push_back(id);
 }
 
+Vector3f TriangleMesh::GetNormal(const Id_t id) const
+{
+    const Vector3i &trg = m_Triangles.at(id);
+    const Vector3f &v0 = m_Points.at(trg(0));
+    const Vector3f &v1 = m_Points.at(trg(1));
+    const Vector3f &v2 = m_Points.at(trg(2));
+    
+    Vector3f edge1 = v1 - v0;
+    Vector3f edge2 = v2 - v0;
+    
+    return edge1.cross(edge2).normalized();
+}
+
 
 }

src/Math/TriangleMesh.h

@@ -47,6 +47,9 @@ public:
     inline std::vector<Vector3f> & Points() { return this->m_Points; }
     inline std::vector<Vector3i> & Triangles() { return this->m_Triangles; }
 
+    const Vector3i & GetTriangle(const Id_t id) const { return m_Triangles.at(id); }
+    Vector3f GetNormal(const Id_t id) const;
+
 private:
     std::vector<Vector3f> m_Points;
     std::vector<Vector3i> m_Triangles;

src/Math/Units.h

@@ -0,0 +1,56 @@
+#ifndef ULIB_MATH_UNITS_H
+#define ULIB_MATH_UNITS_H
+
+
+#ifdef HAVE_GEANT4
+#include <CLHEP/Units/PhysicalConstants.h>
+#include <CLHEP/Units/SystemOfUnits.h>
+#else
+#include "CLHEP_defs.hh"
+#endif
+
+namespace uLib {
+  using namespace CLHEP;
+
+  inline namespace literals {
+    constexpr double operator"" _m(long double v) { return static_cast<double>(v) * CLHEP::meter; }
+    constexpr double operator"" _cm(long double v) { return static_cast<double>(v) * CLHEP::centimeter; }
+    constexpr double operator"" _mm(long double v) { return static_cast<double>(v) * CLHEP::millimeter; }
+    constexpr double operator"" _um(long double v) { return static_cast<double>(v) * CLHEP::micrometer; }
+    constexpr double operator"" _nm(long double v) { return static_cast<double>(v) * CLHEP::nanometer; }
+    constexpr double operator"" _km(long double v) { return static_cast<double>(v) * CLHEP::kilometer; }
+
+    constexpr double operator"" _m(unsigned long long v) { return static_cast<double>(v) * CLHEP::meter; }
+    constexpr double operator"" _cm(unsigned long long v) { return static_cast<double>(v) * CLHEP::centimeter; }
+    constexpr double operator"" _mm(unsigned long long v) { return static_cast<double>(v) * CLHEP::millimeter; }
+    constexpr double operator"" _um(unsigned long long v) { return static_cast<double>(v) * CLHEP::micrometer; }
+    constexpr double operator"" _nm(unsigned long long v) { return static_cast<double>(v) * CLHEP::nanometer; }
+    constexpr double operator"" _km(unsigned long long v) { return static_cast<double>(v) * CLHEP::kilometer; }
+
+    constexpr double operator"" _deg(long double v) { return static_cast<double>(v) * CLHEP::degree; }
+    constexpr double operator"" _rad(long double v) { return static_cast<double>(v) * CLHEP::radian; }
+    constexpr double operator"" _deg(unsigned long long v) { return static_cast<double>(v) * CLHEP::degree; }
+    constexpr double operator"" _rad(unsigned long long v) { return static_cast<double>(v) * CLHEP::radian; }
+
+    constexpr double operator"" _ns(long double v) { return static_cast<double>(v) * CLHEP::nanosecond; }
+    constexpr double operator"" _s(long double v) { return static_cast<double>(v) * CLHEP::second; }
+    constexpr double operator"" _ms(long double v) { return static_cast<double>(v) * CLHEP::millisecond; }
+    constexpr double operator"" _ns(unsigned long long v) { return static_cast<double>(v) * CLHEP::nanosecond; }
+    constexpr double operator"" _s(unsigned long long v) { return static_cast<double>(v) * CLHEP::second; }
+    constexpr double operator"" _ms(unsigned long long v) { return static_cast<double>(v) * CLHEP::millisecond; }
+
+    constexpr double operator"" _MeV(long double v) { return static_cast<double>(v) * CLHEP::megaelectronvolt; }
+    constexpr double operator"" _eV(long double v) { return static_cast<double>(v) * CLHEP::electronvolt; }
+    constexpr double operator"" _keV(long double v) { return static_cast<double>(v) * CLHEP::kiloelectronvolt; }
+    constexpr double operator"" _GeV(long double v) { return static_cast<double>(v) * CLHEP::gigaelectronvolt; }
+    constexpr double operator"" _TeV(long double v) { return static_cast<double>(v) * CLHEP::teraelectronvolt; }
+    constexpr double operator"" _MeV(unsigned long long v) { return static_cast<double>(v) * CLHEP::megaelectronvolt; }
+    constexpr double operator"" _eV(unsigned long long v) { return static_cast<double>(v) * CLHEP::electronvolt; }
+    constexpr double operator"" _keV(unsigned long long v) { return static_cast<double>(v) * CLHEP::kiloelectronvolt; }
+    constexpr double operator"" _GeV(unsigned long long v) { return static_cast<double>(v) * CLHEP::gigaelectronvolt; }
+  }
+}
+
+#endif
+
+

src/Math/VoxRaytracer.h

@@ -63,9 +63,16 @@ public:
 
     inline size_t Count() const { return this->m_Count; }
 
+    inline size_t size() const { return this->m_Count; }
+
     inline const Scalarf &TotalLength() const { return this->m_TotalLength; }
 
-    inline void SetCount(size_t c) { this->m_Count = c; }
+    inline void SetCount(size_t c) {
+      this->m_Count = c;
+      if (this->m_Data.size() != c) {
+        this->m_Data.resize(c);
+      }
+    }
 
     inline void SetTotalLength(Scalarf tl) { this->m_TotalLength = tl; }
 

src/Math/testing/CMakeLists.txt

@@ -12,7 +12,9 @@ set(TESTS
         AccumulatorTest
         VoxImageCopyTest
         TriangleMeshTest
+        QuadMeshTest
         BitCodeTest
+        UnitsTest
 )
 
 set(LIBRARIES

src/Math/testing/ContainerBoxTest.cpp

@@ -31,6 +31,7 @@
 #include "Math/Dense.h"
 #include "Math/ContainerBox.h"
 
+#include <cmath>
 #include <iostream>
 #include <math.h>
 
@@ -52,38 +53,112 @@ int main()
 
   BEGIN_TESTING(Math ContainerBox);
 
-  ContainerBox Cnt;
+  {
-
+    ContainerBox Cnt;
-//  // Local transform:
+    Cnt.SetOrigin(Vector3f(0,0,0));
-  Cnt.SetOrigin(Vector3f(-1,-1,-1));
+    Cnt.SetSize(Vector3f(2,2,2));
-  Cnt.SetSize(Vector3f(2,2,2)); // scaling //
+    TEST0( Vector4f0(Cnt.GetOrigin().homogeneous() - HVector3f(0,0,0)) );
-  std::cout << "Container scale is: " << Cnt.GetSize().transpose() << "\n";
+    TEST0( Vector4f0(Cnt.GetSize().homogeneous() - HVector3f(2,2,2)) );
-  std::cout << "Container scale is: " << Cnt.GetSize().transpose() << "\n";
-  TEST0( Vector4f0(Cnt.GetSize().homogeneous() - HVector3f(2,2,2)) );
 
 
-  ContainerBox Box;
+    HPoint3f pt = Cnt.GetLocalPoint(HPoint3f(0,0,0));
+    HPoint3f wp = Cnt.GetWorldPoint(pt);  
+    TEST0( Vector4f0(wp - HPoint3f(0,0,0)) );
 
-  Box.SetPosition(Vector3f(1,1,1));
+    HPoint3f pt2 = Cnt.GetLocalPoint(HPoint3f(2,2,2));
-  Box.SetSize(Vector3f(2,2,2));
+    HPoint3f wp2 = Cnt.GetWorldPoint(pt2);  
-  Box.EulerYZYRotate(Vector3f(0,0,0));
+    TEST0( Vector4f0(wp2 - HPoint3f(2,2,2)) );
-  HPoint3f pt = Box.GetLocalPoint(HPoint3f(2,3,2));
+
-  HPoint3f wp = Box.GetWorldPoint(pt);
+    HPoint3f pt3 = Cnt.GetLocalPoint(HPoint3f(1,1,1));
-  TEST0( Vector4f0(wp - HPoint3f(2,3,2)) );
+    HPoint3f wp3 = Cnt.GetWorldPoint(pt3);  
+    TEST0( Vector4f0(wp3 - HPoint3f(1,1,1)) );
+
+    HPoint3f pt4 = Cnt.GetLocalPoint(HPoint3f(1,2,3));
+    HPoint3f wp4 = Cnt.GetWorldPoint(pt4);  
+    TEST0( Vector4f0(wp4 - HPoint3f(1,2,3)) );
+  }
+
+  {
+    ContainerBox Cnt;
+    Cnt.SetOrigin(Vector3f(0,0,0));
+    Cnt.SetSize(Vector3f(2,2,2));
+    Cnt.EulerYZYRotate(Vector3f(M_PI,0,0));
+    
+    HPoint3f pt = Cnt.GetLocalPoint(HPoint3f(0,0,0));
+    HPoint3f wp = Cnt.GetWorldPoint(pt);  
+    TEST0( Vector4f0(wp - HPoint3f(0,0,0)) );
+
+    HPoint3f pt2 = Cnt.GetLocalPoint(HPoint3f(2,2,2));
+    HPoint3f wp2 = Cnt.GetWorldPoint(pt2);  
+    TEST0( Vector4f0(wp2 - HPoint3f(2,2,2)) );
+
+    pt2 = HPoint3f(1,1,1);
+    wp2 = Cnt.GetWorldPoint(pt2);  
+    TEST0( Vector4f0(wp2 - HPoint3f(-2,2,-2)) );
+
+    pt2 = HPoint3f(1,2,3);
+    wp2 = Cnt.GetWorldPoint(pt2);  
+    TEST0( Vector4f0(wp2 - HPoint3f(-2,4,-6)) );
+  }
+
+  {
+    ContainerBox Cnt;
+    Cnt.SetOrigin(Vector3f(-1,-1,-1));
+    Cnt.SetSize(Vector3f(2,2,2)); // scaling //
+
+    HPoint3f pt2 = HPoint3f(.5,.5,.5);
+    HPoint3f wp2 = Cnt.GetWorldPoint(pt2);  
+    TEST0( Vector4f0(wp2 - HPoint3f(0,0,0)) );
+
+    pt2 = HPoint3f(0,0,0);
+    wp2 = Cnt.GetWorldPoint(pt2);  
+    TEST0( Vector4f0(wp2 - HPoint3f(-1,-1,-1)) );
+
+    Cnt.EulerYZYRotate(Vector3f(M_PI,0,0));
+    pt2 = HPoint3f(0,0,0);
+    wp2 = Cnt.GetWorldPoint(pt2);  
+    TEST0( Vector4f0(wp2 - HPoint3f(1,-1,1)) );
+  }
+
+  {
+    ContainerBox Box;
+    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));
+    HPoint3f wp = Box.GetWorldPoint(pt);
+    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)) ); 
+
+  }
 
 
-////  // Global
-//  Cnt.SetPosition(Vector3f(1,1,1));
-//  Cnt.EulerYZYRotate(Vector3f(M_PI_2,M_PI_2,0));
-//  HPoint3f p = Cnt.GetWorldPoint(1,1,1);
-//  //std::cout << p.transpose() << "\n";
-//  TEST0( Vector4f0(p - HVector3f(2,1,2)) );
-//  p = Cnt.GetWorldPoint(1,2,3);
-//  //std::cout << p.transpose() << "\n";
-//  TEST0( Vector4f0(p - HVector3f(4,1,3)) );
-
-
-//  // scaling //
 
 
 

src/Math/testing/Makefile.am

@@ -14,7 +14,9 @@ TESTS = MathVectorTest  	\
 	PolicyTest		\
 	AccumulatorTest		\
 	VoxImageCopyTest	\
-	TriangleMeshTest
+	TriangleMeshTest	\
+	BitCodeTest		\
+	UnitsTest
 
 # else
 # TEST =

src/Math/testing/QuadMeshTest.cpp

@@ -23,52 +23,29 @@
 
 //////////////////////////////////////////////////////////////////////////////*/
 
+#include "testing-prototype.h"
 
+#include "Math/QuadMesh.h"
+#include <iostream>
 
+using namespace uLib;
 
+int main() {
+  BEGIN_TESTING(Quad Mesh);
 
-#include <Geant4/G4Material.hh>
+  QuadMesh mesh;
-#include <Geant4/G4NistManager.hh>
-#include <Geant4/G4LogicalVolume.hh>
-
-#include "Core/Vector.h"
-#include "Matter.h"
-#include "Solid.h"
-
-#include "Scene.h"
-
-
-namespace uLib {
-
-
-class DetectorsScenePimpl {
-
-public:
-
-
-    // members //
-    //Vector<Solid> m_Solids;
-};
-
-
-
-
-DetectorsScene::DetectorsScene() :
-    d(new DetectorsScenePimpl())
-{}
-
-DetectorsScene::~DetectorsScene()
-{
-    delete d;
-}
-
-void DetectorsScene::AddSolid(const Solid &solid)
-{
-//    d->m_Solids.push_back(solid);
-}
 
+  mesh.AddPoint(Vector3f(0, 0, 0));
+  mesh.AddPoint(Vector3f(1, 0, 0));
+  mesh.AddPoint(Vector3f(1, 1, 0));
+  mesh.AddPoint(Vector3f(0, 1, 0));
 
+  mesh.AddQuad(Vector4i(0, 1, 2, 3));
 
+  mesh.PrintSelf(std::cout);
 
+  ASSERT_EQ(mesh.Points().size(), 4);
+  ASSERT_EQ(mesh.Quads().size(), 1);
 
+  END_TESTING;
 }

src/Math/testing/UnitsTest.cpp

@@ -0,0 +1,71 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 "testing-prototype.h"
+#include "Math/Units.h"
+#include <cmath>
+
+using namespace uLib;
+
+int main()
+{
+  BEGIN_TESTING(Math Units);
+
+  // Length
+  TEST0(1.0_m  != 1000.0);
+  TEST0(1_m    != 1000.0);
+  TEST0(1.0_cm != 10.0);
+  TEST0(1_cm   != 10.0);
+  TEST0(1.0_mm != 1.0);
+  TEST0(1_mm   != 1.0);
+  TEST0(1.0_km != 1000000.0);
+  TEST0(1_km   != 1000000.0);
+  TEST0(std::abs(1.0_um - 1e-3) > 1e-12);
+  TEST0(std::abs(1.0_nm - 1e-6) > 1e-15);
+
+  // Angle
+  TEST0(std::abs(180.0_deg - M_PI) > 1e-9);
+  TEST0(std::abs(180_deg - M_PI) > 1e-9);
+  TEST0(std::abs(1.0_rad - 1.0) > 1e-9);
+
+  // Time
+  TEST0(1.0_s  != 1e9);
+  TEST0(1_s    != 1e9);
+  TEST0(1.0_ms != 1e6);
+  TEST0(1_ms   != 1e6);
+  TEST0(1.0_ns != 1.0);
+  TEST0(1_ns   != 1.0);
+
+  // Energy
+  TEST0(1.0_MeV != 1.0);
+  TEST0(1_MeV   != 1.0);
+  TEST0(1.0_GeV != 1000.0);
+  TEST0(1_GeV   != 1000.0);
+  TEST0(1.0_TeV != 1000000.0);
+  TEST0(1.0_keV != 1e-3);
+  TEST0(1.0_eV  != 1e-6);
+
+  END_TESTING;
+}

src/Math/testing/VoxRaytracerTestExtended.cpp

@@ -48,7 +48,7 @@ int main() {
 
   Raytracer rt(img);
 
-  const size_t NUM_RAYS = 1000000;
+  const size_t NUM_RAYS = 100000;
   std::cout << "Generating " << NUM_RAYS
             << " random ray pairs across a 100x100x100 grid...\n";
 

src/Math/testing/testing-prototype.h

@@ -31,7 +31,8 @@
 static int _fail = 0;                      \
 printf("..:: Testing " #name " ::..\n");
 
-#define TEST1(val) _fail += (val)==0
+#define TEST1(val) if ((val)==0) { printf("Assertion failed: %s != 0\n", #val); _fail++; }
-#define TEST0(val) _fail += (val)!=0
+#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/ParticlePhysics/Geant/Makefile.am

@@ -1,14 +0,0 @@
-include $(top_srcdir)/Common.am
-
-library_includedir = $(includedir)/libmutom-${PACKAGE_VERSION}/ParticlePhysics/Geant
-
-library_include_HEADERS =
-
-_PPGEANT_SOURCES =
-
-
-
-noinst_LTLIBRARIES = libmutomppgeant.la
-libmutomppgeant_la_SOURCES = ${_PPGEANT_SOURCES}
-
-

src/ParticlePhysics/MuonTomography/Makefile.am

@@ -1,11 +0,0 @@
-include $(top_srcdir)/Common.am
-
-library_includedir = $(includedir)/libmutom-${PACKAGE_VERSION}/ParticlePhysics/MuonTomography
-library_include_HEADERS =
-
-_PPMUTOM_SOURCES =
-
-noinst_LTLIBRARIES = libmutomppmutom.la
-libmutomppmutom_la_SOURCES = ${_PPMUTOM_SOURCES}
-
-

src/Python/CMakeLists.txt

@@ -25,7 +25,7 @@ target_include_directories(uLib_python PRIVATE
 
 # Install uLib_python within the uLib install target
 install(TARGETS uLib_python
-        EXPORT "${PROJECT_NAME}Targets"
+        EXPORT "uLibTargets"
         RUNTIME DESTINATION ${INSTALL_BIN_DIR} COMPONENT bin
         LIBRARY DESTINATION ${INSTALL_LIB_DIR} COMPONENT lib
         ARCHIVE DESTINATION ${INSTALL_LIB_DIR} COMPONENT lib

src/Python/math_bindings.cpp

@@ -1,22 +1,23 @@
-#include <pybind11/pybind11.h>
 #include <pybind11/eigen.h>
+#include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 #include <pybind11/stl_bind.h>
 
 #include <pybind11/numpy.h>
 
-#include "Math/Dense.h"
+#include "Math/Accumulator.h"
-#include "Math/Transform.h"
-#include "Math/Geometry.h"
 #include "Math/ContainerBox.h"
-#include "Math/StructuredData.h"
+#include "Math/Dense.h"
-#include "Math/StructuredGrid.h"
+#include "Math/Geometry.h"
 #include "Math/Structured2DGrid.h"
 #include "Math/Structured4DGrid.h"
+#include "Math/StructuredData.h"
+#include "Math/StructuredGrid.h"
+#include "Math/Transform.h"
 #include "Math/TriangleMesh.h"
-#include "Math/VoxRaytracer.h"
+#include "Math/QuadMesh.h"
-#include "Math/Accumulator.h"
 #include "Math/VoxImage.h"
+#include "Math/VoxRaytracer.h"
 
 namespace py = pybind11;
 using namespace uLib;
@@ -40,346 +41,435 @@ PYBIND11_MAKE_OPAQUE(uLib::Vector<VoxRaytracer::RayData::Element>);
 
 template <typename MatrixType>
 void bind_eigen_type(py::module_ &m, const char *name) {
-    using Scalar = typename MatrixType::Scalar;
+  using Scalar = typename MatrixType::Scalar;
-    constexpr bool is_vector = MatrixType::IsVectorAtCompileTime;
+  constexpr bool is_vector = MatrixType::IsVectorAtCompileTime;
 
-    // Default constructor (zeros)
+  // Default constructor (zeros)
-    m.def(name, []() -> MatrixType {
+  m.def(name, []() -> MatrixType {
-        if constexpr (MatrixType::RowsAtCompileTime == Eigen::Dynamic || MatrixType::ColsAtCompileTime == Eigen::Dynamic) {
+    if constexpr (MatrixType::RowsAtCompileTime == Eigen::Dynamic ||
-            return MatrixType(); // Empty dynamic matrix
+                  MatrixType::ColsAtCompileTime == Eigen::Dynamic) {
-        } else {
+      return MatrixType(); // Empty dynamic matrix
-            return MatrixType::Zero(); // Zero static matrix
+    } else {
-        }
+      return MatrixType::Zero(); // Zero static matrix
-    });
-
-    // Specialized constructor for dynamic matrices
-    if constexpr (MatrixType::RowsAtCompileTime == Eigen::Dynamic || MatrixType::ColsAtCompileTime == Eigen::Dynamic) {
-        m.def(name, [](int rows, int cols) -> MatrixType {
-            MatrixType mat;
-            mat.setZero(rows, cols);
-            return mat;
-        });
     }
+  });
 
-    // Initialize from list
+  // Specialized constructor for dynamic matrices
-    m.def(name, [](py::list l) -> MatrixType {
+  if constexpr (MatrixType::RowsAtCompileTime == Eigen::Dynamic ||
-        MatrixType mat;
+                MatrixType::ColsAtCompileTime == Eigen::Dynamic) {
-        if constexpr (is_vector) {
+    m.def(name, [](int rows, int cols) -> MatrixType {
-            mat.setZero(l.size());
+      MatrixType mat;
-            for (size_t i = 0; i < l.size(); ++i) {
+      mat.setZero(rows, cols);
-                mat(i) = l[i].cast<Scalar>();
+      return mat;
-            }
-        } else {
-            int rows = MatrixType::RowsAtCompileTime == Eigen::Dynamic ? (int)std::sqrt(l.size()) : MatrixType::RowsAtCompileTime;
-            int cols = MatrixType::ColsAtCompileTime == Eigen::Dynamic ? (int)std::sqrt(l.size()) : MatrixType::ColsAtCompileTime;
-            mat.setZero(rows, cols);
-            for (size_t i = 0; i < (size_t)l.size(); ++i) {
-                mat(i / cols, i % cols) = l[i].cast<Scalar>();
-            }
-        }
-        return mat;
     });
+  }
 
-    // Initialize from py::array
+  // Initialize from list
-    m.def(name, [](py::array_t<Scalar, py::array::c_style | py::array::forcecast> arr) -> MatrixType {
+  m.def(name, [](py::list l) -> MatrixType {
-        auto buf = arr.request();
+    MatrixType mat;
-        MatrixType mat;
+    if constexpr (is_vector) {
-        if constexpr (is_vector) {
+      mat.setZero(l.size());
+      for (size_t i = 0; i < l.size(); ++i) {
+        mat(i) = l[i].cast<Scalar>();
+      }
+    } else {
+      int rows = MatrixType::RowsAtCompileTime == Eigen::Dynamic
+                     ? (int)std::sqrt(l.size())
+                     : MatrixType::RowsAtCompileTime;
+      int cols = MatrixType::ColsAtCompileTime == Eigen::Dynamic
+                     ? (int)std::sqrt(l.size())
+                     : MatrixType::ColsAtCompileTime;
+      mat.setZero(rows, cols);
+      for (size_t i = 0; i < (size_t)l.size(); ++i) {
+        mat(i / cols, i % cols) = l[i].cast<Scalar>();
+      }
+    }
+    return mat;
+  });
+
+  // Initialize from py::array
+  m.def(name,
+        [](py::array_t<Scalar, py::array::c_style | py::array::forcecast> arr)
+            -> MatrixType {
+          auto buf = arr.request();
+          MatrixType mat;
+          if constexpr (is_vector) {
             mat.setZero(buf.size);
-            Scalar* ptr = static_cast<Scalar*>(buf.ptr);
+            Scalar *ptr = static_cast<Scalar *>(buf.ptr);
-            for (ssize_t i = 0; i < buf.size; ++i) mat(i) = ptr[i];
+            for (ssize_t i = 0; i < buf.size; ++i)
-        } else {
+              mat(i) = ptr[i];
+          } else {
             int rows = buf.shape.size() > 0 ? (int)buf.shape[0] : 1;
             int cols = buf.shape.size() > 1 ? (int)buf.shape[1] : 1;
             mat.setZero(rows, cols);
-            Scalar* ptr = static_cast<Scalar*>(buf.ptr);
+            Scalar *ptr = static_cast<Scalar *>(buf.ptr);
             for (int i = 0; i < rows; ++i) {
-                for (int j = 0; j < cols; ++j) {
+              for (int j = 0; j < cols; ++j) {
-                    mat(i, j) = ptr[i * cols + j];
+                mat(i, j) = ptr[i * cols + j];
-                }
+              }
             }
-        }
+          }
-        return mat;
+          return mat;
-    });
+        });
 }
 
 void init_math(py::module_ &m) {
 
-    // 1. Basic Eigen Types (Vectors and Matrices)
+  // 1. Basic Eigen Types (Vectors and Matrices)
-    bind_eigen_type<Vector1f>(m, "Vector1f");
+  bind_eigen_type<Vector1f>(m, "Vector1f");
-    bind_eigen_type<Vector2f>(m, "Vector2f");
+  bind_eigen_type<Vector2f>(m, "Vector2f");
-    bind_eigen_type<Vector3f>(m, "Vector3f");
+  bind_eigen_type<Vector3f>(m, "Vector3f");
-    bind_eigen_type<Vector4f>(m, "Vector4f");
+  bind_eigen_type<Vector4f>(m, "Vector4f");
-    bind_eigen_type<Vector1i>(m, "Vector1i");
+  bind_eigen_type<Vector1i>(m, "Vector1i");
-    bind_eigen_type<Vector2i>(m, "Vector2i");
+  bind_eigen_type<Vector2i>(m, "Vector2i");
-    bind_eigen_type<Vector3i>(m, "Vector3i");
+  bind_eigen_type<Vector3i>(m, "Vector3i");
-    bind_eigen_type<Vector4i>(m, "Vector4i");
+  bind_eigen_type<Vector4i>(m, "Vector4i");
-    bind_eigen_type<Vector1d>(m, "Vector1d");
+  bind_eigen_type<Vector1d>(m, "Vector1d");
-    bind_eigen_type<Vector2d>(m, "Vector2d");
+  bind_eigen_type<Vector2d>(m, "Vector2d");
-    bind_eigen_type<Vector3d>(m, "Vector3d");
+  bind_eigen_type<Vector3d>(m, "Vector3d");
-    bind_eigen_type<Vector4d>(m, "Vector4d");
+  bind_eigen_type<Vector4d>(m, "Vector4d");
 
-    bind_eigen_type<Matrix2f>(m, "Matrix2f");
+  bind_eigen_type<Matrix2f>(m, "Matrix2f");
-    bind_eigen_type<Matrix3f>(m, "Matrix3f");
+  bind_eigen_type<Matrix3f>(m, "Matrix3f");
-    bind_eigen_type<Matrix4f>(m, "Matrix4f");
+  bind_eigen_type<Matrix4f>(m, "Matrix4f");
-    bind_eigen_type<Matrix2i>(m, "Matrix2i");
+  bind_eigen_type<Matrix2i>(m, "Matrix2i");
-    bind_eigen_type<Matrix3i>(m, "Matrix3i");
+  bind_eigen_type<Matrix3i>(m, "Matrix3i");
-    bind_eigen_type<Matrix4i>(m, "Matrix4i");
+  bind_eigen_type<Matrix4i>(m, "Matrix4i");
-    bind_eigen_type<Matrix2d>(m, "Matrix2d");
+  bind_eigen_type<Matrix2d>(m, "Matrix2d");
-    bind_eigen_type<Matrix3d>(m, "Matrix3d");
+  bind_eigen_type<Matrix3d>(m, "Matrix3d");
-    bind_eigen_type<Matrix4d>(m, "Matrix4d");
+  bind_eigen_type<Matrix4d>(m, "Matrix4d");
 
-    bind_eigen_type<MatrixXi>(m, "MatrixXi");
+  bind_eigen_type<MatrixXi>(m, "MatrixXi");
-    bind_eigen_type<MatrixXf>(m, "MatrixXf");
+  bind_eigen_type<MatrixXf>(m, "MatrixXf");
-    bind_eigen_type<MatrixXd>(m, "MatrixXd");
+  bind_eigen_type<MatrixXd>(m, "MatrixXd");
 
-    // 2. Homogeneous types
+  // 2. Homogeneous types
-    py::class_<HPoint3f>(m, "HPoint3f")
+  py::class_<HPoint3f>(m, "HPoint3f")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def(py::init<float, float, float>())
+      .def(py::init<float, float, float>())
-        .def(py::init<Vector3f &>());
+      .def(py::init<Vector3f &>());
-    py::class_<HVector3f>(m, "HVector3f")
+  py::class_<HVector3f>(m, "HVector3f")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def(py::init<float, float, float>())
+      .def(py::init<float, float, float>())
-        .def(py::init<Vector3f &>());
+      .def(py::init<Vector3f &>());
-    py::class_<HLine3f>(m, "HLine3f")
+  py::class_<HLine3f>(m, "HLine3f")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def_readwrite("origin", &HLine3f::origin)
+      .def_readwrite("origin", &HLine3f::origin)
-        .def_readwrite("direction", &HLine3f::direction);
+      .def_readwrite("direction", &HLine3f::direction);
-    py::class_<HError3f>(m, "HError3f")
+  py::class_<HError3f>(m, "HError3f")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def_readwrite("position_error", &HError3f::position_error)
+      .def_readwrite("position_error", &HError3f::position_error)
-        .def_readwrite("direction_error", &HError3f::direction_error);
+      .def_readwrite("direction_error", &HError3f::direction_error);
 
-    // 3. Dynamic Vectors (uLib::Vector)
+  // 3. Dynamic Vectors (uLib::Vector)
-    py::bind_vector<uLib::Vector<Scalari>>(m, "Vector_i")
+  py::bind_vector<uLib::Vector<Scalari>>(m, "Vector_i")
-        .def("MoveToVRAM", &uLib::Vector<Scalari>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Scalari>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Scalari>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Scalari>::MoveToRAM);
-    py::bind_vector<uLib::Vector<Scalarui>>(m, "Vector_ui")
+  py::bind_vector<uLib::Vector<Scalarui>>(m, "Vector_ui")
-        .def("MoveToVRAM", &uLib::Vector<Scalarui>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Scalarui>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Scalarui>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Scalarui>::MoveToRAM);
-    py::bind_vector<uLib::Vector<Scalarl>>(m, "Vector_l")
+  py::bind_vector<uLib::Vector<Scalarl>>(m, "Vector_l")
-        .def("MoveToVRAM", &uLib::Vector<Scalarl>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Scalarl>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Scalarl>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Scalarl>::MoveToRAM);
-    py::bind_vector<uLib::Vector<Scalarul>>(m, "Vector_ul")
+  py::bind_vector<uLib::Vector<Scalarul>>(m, "Vector_ul")
-        .def("MoveToVRAM", &uLib::Vector<Scalarul>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Scalarul>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Scalarul>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Scalarul>::MoveToRAM);
-    py::bind_vector<uLib::Vector<Scalarf>>(m, "Vector_f")
+  py::bind_vector<uLib::Vector<Scalarf>>(m, "Vector_f")
-        .def("MoveToVRAM", &uLib::Vector<Scalarf>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Scalarf>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Scalarf>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Scalarf>::MoveToRAM);
-    py::bind_vector<uLib::Vector<Scalard>>(m, "Vector_d")
+  py::bind_vector<uLib::Vector<Scalard>>(m, "Vector_d")
-        .def("MoveToVRAM", &uLib::Vector<Scalard>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Scalard>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Scalard>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Scalard>::MoveToRAM);
 
-    py::bind_vector<uLib::Vector<Vector3f>>(m, "Vector_Vector3f")
+  py::bind_vector<uLib::Vector<Vector3f>>(m, "Vector_Vector3f")
-        .def("MoveToVRAM", &uLib::Vector<Vector3f>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Vector3f>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Vector3f>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Vector3f>::MoveToRAM);
-    py::bind_vector<uLib::Vector<Vector3i>>(m, "Vector_Vector3i")
+  py::bind_vector<uLib::Vector<Vector3i>>(m, "Vector_Vector3i")
-        .def("MoveToVRAM", &uLib::Vector<Vector3i>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Vector3i>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Vector3i>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Vector3i>::MoveToRAM);
-    py::bind_vector<uLib::Vector<Vector4f>>(m, "Vector_Vector4f")
+  py::bind_vector<uLib::Vector<Vector4f>>(m, "Vector_Vector4f")
-        .def("MoveToVRAM", &uLib::Vector<Vector4f>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Vector4f>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Vector4f>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Vector4f>::MoveToRAM);
-    py::bind_vector<uLib::Vector<Vector4i>>(m, "Vector_Vector4i")
+  py::bind_vector<uLib::Vector<Vector4i>>(m, "Vector_Vector4i")
-        .def("MoveToVRAM", &uLib::Vector<Vector4i>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Vector4i>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Vector4i>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Vector4i>::MoveToRAM);
-    py::bind_vector<uLib::Vector<Vector3d>>(m, "Vector_Vector3d")
+  py::bind_vector<uLib::Vector<Vector3d>>(m, "Vector_Vector3d")
-        .def("MoveToVRAM", &uLib::Vector<Vector3d>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Vector3d>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Vector3d>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Vector3d>::MoveToRAM);
-    py::bind_vector<uLib::Vector<Vector4d>>(m, "Vector_Vector4d")
+  py::bind_vector<uLib::Vector<Vector4d>>(m, "Vector_Vector4d")
-        .def("MoveToVRAM", &uLib::Vector<Vector4d>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Vector4d>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Vector4d>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Vector4d>::MoveToRAM);
-    py::bind_vector<uLib::Vector<Voxel>>(m, "Vector_Voxel")
+  py::bind_vector<uLib::Vector<Voxel>>(m, "Vector_Voxel")
-        .def("MoveToVRAM", &uLib::Vector<Voxel>::MoveToVRAM)
+      .def("MoveToVRAM", &uLib::Vector<Voxel>::MoveToVRAM)
-        .def("MoveToRAM", &uLib::Vector<Voxel>::MoveToRAM);
+      .def("MoveToRAM", &uLib::Vector<Voxel>::MoveToRAM);
-    py::bind_vector<uLib::Vector<VoxRaytracer::RayData::Element>>(m, "Vector_VoxRaytracerRayDataElement")
+  py::bind_vector<uLib::Vector<VoxRaytracer::RayData::Element>>(
-        .def("MoveToVRAM", &uLib::Vector<VoxRaytracer::RayData::Element>::MoveToVRAM)
+      m, "Vector_VoxRaytracerRayDataElement")
-        .def("MoveToRAM", &uLib::Vector<VoxRaytracer::RayData::Element>::MoveToRAM);
+      .def("MoveToVRAM",
+           &uLib::Vector<VoxRaytracer::RayData::Element>::MoveToVRAM)
+      .def("MoveToRAM",
+           &uLib::Vector<VoxRaytracer::RayData::Element>::MoveToRAM);
 
-    // 4. Accumulators
+  // 4. Accumulators
-    py::class_<Accumulator_Mean<float>>(m, "Accumulator_Mean_f")
+  py::class_<Accumulator_Mean<float>>(m, "Accumulator_Mean_f")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def("AddPass", &Accumulator_Mean<float>::AddPass)
+      .def("AddPass", &Accumulator_Mean<float>::AddPass)
-        .def("__call__", py::overload_cast<const float>(&Accumulator_Mean<float>::operator()))
+      .def("__call__",
-        .def("__call__", py::overload_cast<>(&Accumulator_Mean<float>::operator(), py::const_));
+           py::overload_cast<const float>(&Accumulator_Mean<float>::operator()))
+      .def("__call__", py::overload_cast<>(&Accumulator_Mean<float>::operator(),
+                                           py::const_));
 
-    py::class_<Accumulator_Mean<double>>(m, "Accumulator_Mean_d")
+  py::class_<Accumulator_Mean<double>>(m, "Accumulator_Mean_d")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def("AddPass", &Accumulator_Mean<double>::AddPass)
+      .def("AddPass", &Accumulator_Mean<double>::AddPass)
-        .def("__call__", py::overload_cast<const double>(&Accumulator_Mean<double>::operator()))
+      .def("__call__", py::overload_cast<const double>(
-        .def("__call__", py::overload_cast<>(&Accumulator_Mean<double>::operator(), py::const_));
+                           &Accumulator_Mean<double>::operator()))
+      .def("__call__", py::overload_cast<>(
+                           &Accumulator_Mean<double>::operator(), py::const_));
 
-    py::class_<Accumulator_ABTrim<float>>(m, "Accumulator_ABTrim_f")
+  py::class_<Accumulator_ABTrim<float>>(m, "Accumulator_ABTrim_f")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def("SetABTrim", &Accumulator_ABTrim<float>::SetABTrim)
+      .def("SetABTrim", &Accumulator_ABTrim<float>::SetABTrim)
-        .def("__iadd__", [](Accumulator_ABTrim<float> &self, float val) { self += val; return &self; })
+      .def("__iadd__",
-        .def("__call__", &Accumulator_ABTrim<float>::operator());
+           [](Accumulator_ABTrim<float> &self, float val) {
+             self += val;
+             return &self;
+           })
+      .def("__call__", &Accumulator_ABTrim<float>::operator());
 
-    py::class_<Accumulator_ABTrim<double>>(m, "Accumulator_ABTrim_d")
+  py::class_<Accumulator_ABTrim<double>>(m, "Accumulator_ABTrim_d")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def("SetABTrim", &Accumulator_ABTrim<double>::SetABTrim)
+      .def("SetABTrim", &Accumulator_ABTrim<double>::SetABTrim)
-        .def("__iadd__", [](Accumulator_ABTrim<double> &self, double val) { self += val; return &self; })
+      .def("__iadd__",
-        .def("__call__", &Accumulator_ABTrim<double>::operator());
+           [](Accumulator_ABTrim<double> &self, double val) {
+             self += val;
+             return &self;
+           })
+      .def("__call__", &Accumulator_ABTrim<double>::operator());
 
-    py::class_<Accumulator_ABClip<float>>(m, "Accumulator_ABClip_f")
+  py::class_<Accumulator_ABClip<float>>(m, "Accumulator_ABClip_f")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def("SetABTrim", &Accumulator_ABClip<float>::SetABTrim)
+      .def("SetABTrim", &Accumulator_ABClip<float>::SetABTrim)
-        .def("__iadd__", [](Accumulator_ABClip<float> &self, float val) { self += val; return &self; })
+      .def("__iadd__",
-        .def("__call__", &Accumulator_ABClip<float>::operator());
+           [](Accumulator_ABClip<float> &self, float val) {
+             self += val;
+             return &self;
+           })
+      .def("__call__", &Accumulator_ABClip<float>::operator());
 
-    py::class_<Accumulator_ABClip<double>>(m, "Accumulator_ABClip_d")
+  py::class_<Accumulator_ABClip<double>>(m, "Accumulator_ABClip_d")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def("SetABTrim", &Accumulator_ABClip<double>::SetABTrim)
+      .def("SetABTrim", &Accumulator_ABClip<double>::SetABTrim)
-        .def("__iadd__", [](Accumulator_ABClip<double> &self, double val) { self += val; return &self; })
+      .def("__iadd__",
-        .def("__call__", &Accumulator_ABClip<double>::operator());
+           [](Accumulator_ABClip<double> &self, double val) {
+             self += val;
+             return &self;
+           })
+      .def("__call__", &Accumulator_ABClip<double>::operator());
 
-    // 5. Core Math Structures
+  // 5. Core Math Structures
-    py::class_<AffineTransform>(m, "AffineTransform")
+  py::class_<AffineTransform, std::shared_ptr<AffineTransform>>(m,
-        .def(py::init<>())
+                                                              "AffineTransform")
-        .def("GetWorldMatrix", &AffineTransform::GetWorldMatrix)
+      .def(py::init<>())
-        .def("SetPosition", &AffineTransform::SetPosition)
+      .def("GetWorldMatrix", &AffineTransform::GetWorldMatrix)
-        .def("GetPosition", &AffineTransform::GetPosition)
+      .def("SetPosition", &AffineTransform::SetPosition)
-        .def("Translate", &AffineTransform::Translate)
+      .def("GetPosition", &AffineTransform::GetPosition)
-        .def("Scale", &AffineTransform::Scale)
+      .def("Translate", &AffineTransform::Translate)
-        .def("SetRotation", &AffineTransform::SetRotation)
+      .def("Scale", &AffineTransform::Scale)
-        .def("GetRotation", &AffineTransform::GetRotation)
+      .def("SetRotation", &AffineTransform::SetRotation)
-        .def("Rotate", py::overload_cast<const Matrix3f &>(&AffineTransform::Rotate))
+      .def("GetRotation", &AffineTransform::GetRotation)
-        .def("Rotate", py::overload_cast<const Vector3f>(&AffineTransform::Rotate))
+      .def("Rotate",
-        .def("EulerYZYRotate", &AffineTransform::EulerYZYRotate)
+           py::overload_cast<const Matrix3f>(&AffineTransform::Rotate))
-        .def("FlipAxes", &AffineTransform::FlipAxes);
+      .def("Rotate",
+           py::overload_cast<float, Vector3f>(&AffineTransform::Rotate))
+      .def("Rotate", py::overload_cast<Vector3f>(&AffineTransform::Rotate))
+      .def("EulerYZYRotate", &AffineTransform::EulerYZYRotate)
+      .def("FlipAxes", &AffineTransform::FlipAxes);
 
-    py::class_<Geometry, AffineTransform>(m, "Geometry")
+  py::class_<Geometry, AffineTransform, std::shared_ptr<Geometry>>(m, "Geometry")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def("GetWorldPoint", py::overload_cast<const Vector4f &>(&Geometry::GetWorldPoint, py::const_))
+      .def("GetWorldPoint", py::overload_cast<const Vector4f>(
-        .def("GetLocalPoint", py::overload_cast<const Vector4f &>(&Geometry::GetLocalPoint, py::const_));
+                                &Geometry::GetWorldPoint, py::const_))
+      .def("GetWorldPoint",
+           py::overload_cast<float, float, float>(&Geometry::GetWorldPoint))
+      .def("GetLocalPoint", py::overload_cast<const Vector4f>(
+                                &Geometry::GetLocalPoint, py::const_))
+      .def("GetLocalPoint",
+           py::overload_cast<float, float, float>(&Geometry::GetLocalPoint));
 
-    py::class_<ContainerBox, AffineTransform>(m, "ContainerBox")
+  py::class_<ContainerBox, AffineTransform, Object, std::shared_ptr<ContainerBox>>(
-        .def(py::init<>())
+      m, "ContainerBox")
-        .def("SetOrigin", &ContainerBox::SetOrigin)
+      .def(py::init<>())
-        .def("GetOrigin", &ContainerBox::GetOrigin)
+      .def("SetOrigin", &ContainerBox::SetOrigin)
-        .def("SetSize", &ContainerBox::SetSize)
+      .def("GetOrigin", &ContainerBox::GetOrigin)
-        .def("GetSize", &ContainerBox::GetSize)
+      .def("SetSize", &ContainerBox::SetSize)
-        .def("GetWorldMatrix", &ContainerBox::GetWorldMatrix)
+      .def("GetSize", &ContainerBox::GetSize)
-        .def("GetWorldPoint", py::overload_cast<const Vector4f &>(&ContainerBox::GetWorldPoint, py::const_))
+      .def("GetLocalMatrix", &ContainerBox::GetLocalMatrix)
-        .def("GetLocalPoint", py::overload_cast<const Vector4f &>(&ContainerBox::GetLocalPoint, py::const_));
+      .def("GetWorldMatrix", &ContainerBox::GetWorldMatrix)
+      .def("GetWorldPoint", py::overload_cast<const Vector4f &>(
+                                &ContainerBox::GetWorldPoint, py::const_))
+      .def("GetWorldPoint",
+           py::overload_cast<float, float, float>(&ContainerBox::GetWorldPoint))
+      .def("GetLocalPoint", py::overload_cast<const Vector4f &>(
+                                &ContainerBox::GetLocalPoint, py::const_))
+      .def("GetLocalPoint",
+           py::overload_cast<float, float, float>(&ContainerBox::GetLocalPoint))
+      .def("FlipLocalAxes", &ContainerBox::FlipLocalAxes);
 
-    py::enum_<StructuredData::_Order>(m, "StructuredDataOrder")
+  py::enum_<StructuredData::_Order>(m, "StructuredDataOrder")
-        .value("CustomOrder", StructuredData::CustomOrder)
+      .value("CustomOrder", StructuredData::CustomOrder)
-        .value("XYZ", StructuredData::XYZ)
+      .value("XYZ", StructuredData::XYZ)
-        .value("XZY", StructuredData::XZY)
+      .value("XZY", StructuredData::XZY)
-        .value("YXZ", StructuredData::YXZ)
+      .value("YXZ", StructuredData::YXZ)
-        .value("YZX", StructuredData::YZX)
+      .value("YZX", StructuredData::YZX)
-        .value("ZXY", StructuredData::ZXY)
+      .value("ZXY", StructuredData::ZXY)
-        .value("ZYX", StructuredData::ZYX)
+      .value("ZYX", StructuredData::ZYX)
-        .export_values();
+      .export_values();
 
-    py::class_<StructuredData>(m, "StructuredData")
+  py::class_<StructuredData, std::shared_ptr<StructuredData>>(m, "StructuredData")
-        .def(py::init<const Vector3i &>())
+      .def(py::init<const Vector3i &>())
-        .def("GetDims", &StructuredData::GetDims)
+      .def("GetDims", &StructuredData::GetDims)
-        .def("SetDims", &StructuredData::SetDims)
+      .def("SetDims", &StructuredData::SetDims)
-        .def("GetIncrements", &StructuredData::GetIncrements)
+      .def("GetIncrements", &StructuredData::GetIncrements)
-        .def("SetIncrements", &StructuredData::SetIncrements)
+      .def("SetIncrements", &StructuredData::SetIncrements)
-        .def("SetDataOrder", &StructuredData::SetDataOrder)
+      .def("SetDataOrder", &StructuredData::SetDataOrder)
-        .def("GetDataOrder", &StructuredData::GetDataOrder)
+      .def("GetDataOrder", &StructuredData::GetDataOrder)
-        .def("IsInsideGrid", &StructuredData::IsInsideGrid)
+      .def("IsInsideGrid", &StructuredData::IsInsideGrid)
-        .def("Map", &StructuredData::Map)
+      .def("Map", &StructuredData::Map)
-        .def("UnMap", &StructuredData::UnMap);
+      .def("UnMap", &StructuredData::UnMap);
 
-    py::class_<StructuredGrid, ContainerBox, StructuredData>(m, "StructuredGrid")
+  py::class_<StructuredGrid, ContainerBox, StructuredData,
-        .def(py::init<const Vector3i &>())
+             std::shared_ptr<StructuredGrid>>(m, "StructuredGrid")
-        .def("SetSpacing", &StructuredGrid::SetSpacing)
+      .def(py::init<const Vector3i &>())
-        .def("GetSpacing", &StructuredGrid::GetSpacing)
+      .def("SetSpacing", &StructuredGrid::SetSpacing)
-        .def("IsInsideBounds", &StructuredGrid::IsInsideBounds)
+      .def("GetSpacing", &StructuredGrid::GetSpacing)
-        .def("Find", [](StructuredGrid &self, Vector3f pt) {
+      .def("IsInsideBounds", &StructuredGrid::IsInsideBounds)
-            return self.Find(HPoint3f(pt));
+      .def("Find", [](StructuredGrid &self, Vector3f pt) {
-        });
+        return self.Find(HPoint3f(pt));
+      });
 
-    py::class_<Structured2DGrid>(m, "Structured2DGrid")
+  py::class_<Structured2DGrid>(m, "Structured2DGrid")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def("SetDims", &Structured2DGrid::SetDims)
+      .def("SetDims", &Structured2DGrid::SetDims)
-        .def("GetDims", &Structured2DGrid::GetDims)
+      .def("GetDims", &Structured2DGrid::GetDims)
-        .def("IsInsideGrid", &Structured2DGrid::IsInsideGrid)
+      .def("IsInsideGrid", &Structured2DGrid::IsInsideGrid)
-        .def("Map", &Structured2DGrid::Map)
+      .def("Map", &Structured2DGrid::Map)
-        .def("UnMap", &Structured2DGrid::UnMap)
+      .def("UnMap", &Structured2DGrid::UnMap)
-        .def("SetPhysicalSpace", &Structured2DGrid::SetPhysicalSpace)
+      .def("SetPhysicalSpace", &Structured2DGrid::SetPhysicalSpace)
-        .def("GetSpacing", &Structured2DGrid::GetSpacing)
+      .def("GetSpacing", &Structured2DGrid::GetSpacing)
-        .def("GetOrigin", &Structured2DGrid::GetOrigin)
+      .def("GetOrigin", &Structured2DGrid::GetOrigin)
-        .def("IsInsideBounds", &Structured2DGrid::IsInsideBounds)
+      .def("IsInsideBounds", &Structured2DGrid::IsInsideBounds)
-        .def("PhysicsToUnitSpace", &Structured2DGrid::PhysicsToUnitSpace)
+      .def("PhysicsToUnitSpace", &Structured2DGrid::PhysicsToUnitSpace)
-        .def("UnitToPhysicsSpace", &Structured2DGrid::UnitToPhysicsSpace)
+      .def("UnitToPhysicsSpace", &Structured2DGrid::UnitToPhysicsSpace)
-        .def("SetDebug", &Structured2DGrid::SetDebug);
+      .def("SetDebug", &Structured2DGrid::SetDebug);
 
-    py::class_<Structured4DGrid>(m, "Structured4DGrid")
+  py::class_<Structured4DGrid>(m, "Structured4DGrid")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def("SetDims", &Structured4DGrid::SetDims)
+      .def("SetDims", &Structured4DGrid::SetDims)
-        .def("GetDims", &Structured4DGrid::GetDims)
+      .def("GetDims", &Structured4DGrid::GetDims)
-        .def("IsInsideGrid", &Structured4DGrid::IsInsideGrid)
+      .def("IsInsideGrid", &Structured4DGrid::IsInsideGrid)
-        .def("Map", &Structured4DGrid::Map)
+      .def("Map", &Structured4DGrid::Map)
-        .def("UnMap", &Structured4DGrid::UnMap)
+      .def("UnMap", &Structured4DGrid::UnMap)
-        .def("SetPhysicalSpace", &Structured4DGrid::SetPhysicalSpace)
+      .def("SetPhysicalSpace", &Structured4DGrid::SetPhysicalSpace)
-        .def("GetSpacing", &Structured4DGrid::GetSpacing)
+      .def("GetSpacing", &Structured4DGrid::GetSpacing)
-        .def("GetOrigin", &Structured4DGrid::GetOrigin)
+      .def("GetOrigin", &Structured4DGrid::GetOrigin)
-        .def("IsInsideBounds", &Structured4DGrid::IsInsideBounds)
+      .def("IsInsideBounds", &Structured4DGrid::IsInsideBounds)
-        .def("PhysicsToUnitSpace", &Structured4DGrid::PhysicsToUnitSpace)
+      .def("PhysicsToUnitSpace", &Structured4DGrid::PhysicsToUnitSpace)
-        .def("UnitToPhysicsSpace", &Structured4DGrid::UnitToPhysicsSpace)
+      .def("UnitToPhysicsSpace", &Structured4DGrid::UnitToPhysicsSpace)
-        .def("SetDebug", &Structured4DGrid::SetDebug);
+      .def("SetDebug", &Structured4DGrid::SetDebug);
 
-    // 6. High-level Structures
+  // 6. High-level Structures
-    py::class_<Voxel>(m, "Voxel")
+  py::class_<Voxel>(m, "Voxel")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def_readwrite("Value", &Voxel::Value)
+      .def_readwrite("Value", &Voxel::Value)
-        .def_readwrite("Count", &Voxel::Count);
+      .def_readwrite("Count", &Voxel::Count);
 
-    py::class_<Abstract::VoxImage, StructuredGrid>(m, "AbstractVoxImage")
+  py::class_<Abstract::VoxImage, StructuredGrid,
-        .def("GetValue", py::overload_cast<const Vector3i &>(&Abstract::VoxImage::GetValue, py::const_))
+             std::shared_ptr<Abstract::VoxImage>>(m, "AbstractVoxImage")
-        .def("GetValue", py::overload_cast<const int>(&Abstract::VoxImage::GetValue, py::const_))
+      .def("GetValue", py::overload_cast<const Vector3i &>(
-        .def("SetValue", py::overload_cast<const Vector3i &, float>(&Abstract::VoxImage::SetValue))
+                           &Abstract::VoxImage::GetValue, py::const_))
-        .def("SetValue", py::overload_cast<const int, float>(&Abstract::VoxImage::SetValue))
+      .def("GetValue", py::overload_cast<const int>(
-        .def("ExportToVtk", &Abstract::VoxImage::ExportToVtk)
+                           &Abstract::VoxImage::GetValue, py::const_))
-        .def("ExportToVti", &Abstract::VoxImage::ExportToVti)
+      .def("SetValue", py::overload_cast<const Vector3i &, float>(
-        .def("ImportFromVtk", &Abstract::VoxImage::ImportFromVtk)
+                           &Abstract::VoxImage::SetValue))
-        .def("ImportFromVti", &Abstract::VoxImage::ImportFromVti);
+      .def("SetValue",
+           py::overload_cast<const int, float>(&Abstract::VoxImage::SetValue))
+      .def("ExportToVtk", &Abstract::VoxImage::ExportToVtk)
+      .def("ExportToVti", &Abstract::VoxImage::ExportToVti)
+      .def("ImportFromVtk", &Abstract::VoxImage::ImportFromVtk)
+      .def("ImportFromVti", &Abstract::VoxImage::ImportFromVti);
 
-    py::class_<VoxImage<Voxel>, Abstract::VoxImage>(m, "VoxImage")
+  py::class_<VoxImage<Voxel>, Abstract::VoxImage,
-        .def(py::init<>())
+             std::shared_ptr<VoxImage<Voxel>>>(m, "VoxImage")
-        .def(py::init<const Vector3i &>())
+      .def(py::init<>())
-        .def("Data", &VoxImage<Voxel>::Data, py::return_value_policy::reference_internal)
+      .def(py::init<const Vector3i &>())
-        .def("InitVoxels", &VoxImage<Voxel>::InitVoxels)
+      .def("Data", &VoxImage<Voxel>::Data,
-        .def("Abs", &VoxImage<Voxel>::Abs)
+           py::return_value_policy::reference_internal)
-        .def("clipImage", py::overload_cast<const Vector3i, const Vector3i>(&VoxImage<Voxel>::clipImage, py::const_))
+      .def("InitVoxels", &VoxImage<Voxel>::InitVoxels)
-        .def("clipImage", py::overload_cast<const HPoint3f, const HPoint3f>(&VoxImage<Voxel>::clipImage, py::const_))
+      .def("Abs", &VoxImage<Voxel>::Abs)
-        .def("clipImage", py::overload_cast<const float>(&VoxImage<Voxel>::clipImage, py::const_))
+      .def("clipImage", py::overload_cast<const Vector3i, const Vector3i>(
-        .def("maskImage", py::overload_cast<const HPoint3f, const HPoint3f, float>(&VoxImage<Voxel>::maskImage, py::const_))
+                            &VoxImage<Voxel>::clipImage, py::const_))
-        .def("maskImage", py::overload_cast<const float, float, float>(&VoxImage<Voxel>::maskImage, py::const_), py::arg("threshold"), py::arg("belowValue") = 0, py::arg("aboveValue") = 0)
+      .def("clipImage", py::overload_cast<const HPoint3f, const HPoint3f>(
-        .def("fixVoxels", py::overload_cast<const float, float>(&VoxImage<Voxel>::fixVoxels, py::const_))
+                            &VoxImage<Voxel>::clipImage, py::const_))
-        .def("__getitem__", py::overload_cast<unsigned int>(&VoxImage<Voxel>::operator[]))
+      .def("clipImage", py::overload_cast<const float>(
-        .def("__getitem__", py::overload_cast<const Vector3i &>(&VoxImage<Voxel>::operator[]));
+                            &VoxImage<Voxel>::clipImage, py::const_))
+      .def("maskImage",
+           py::overload_cast<const HPoint3f, const HPoint3f, float>(
+               &VoxImage<Voxel>::maskImage, py::const_))
+      .def("maskImage",
+           py::overload_cast<const float, float, float>(
+               &VoxImage<Voxel>::maskImage, py::const_),
+           py::arg("threshold"), py::arg("belowValue") = 0,
+           py::arg("aboveValue") = 0)
+      .def("fixVoxels", py::overload_cast<const float, float>(
+                            &VoxImage<Voxel>::fixVoxels, py::const_))
+      .def("__getitem__",
+           py::overload_cast<unsigned int>(&VoxImage<Voxel>::operator[]))
+      .def("__getitem__",
+           py::overload_cast<const Vector3i &>(&VoxImage<Voxel>::operator[]));
 
-    py::class_<TriangleMesh>(m, "TriangleMesh")
+  py::class_<TriangleMesh>(m, "TriangleMesh")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def("AddPoint", &TriangleMesh::AddPoint)
+      .def("AddPoint", &TriangleMesh::AddPoint)
-        .def("AddTriangle", py::overload_cast<const Vector3i &>(&TriangleMesh::AddTriangle))
+      .def("AddTriangle",
-        .def("Points", &TriangleMesh::Points, py::return_value_policy::reference_internal)
+           py::overload_cast<const Vector3i &>(&TriangleMesh::AddTriangle))
-        .def("Triangles", &TriangleMesh::Triangles, py::return_value_policy::reference_internal);
+      .def("Points", &TriangleMesh::Points,
+           py::return_value_policy::reference_internal)
+      .def("Triangles", &TriangleMesh::Triangles,
+           py::return_value_policy::reference_internal)
+      .def("GetTriangle", &TriangleMesh::GetTriangle)
+      .def("GetNormal", &TriangleMesh::GetNormal);
 
-    py::class_<VoxRaytracer::RayData::Element>(m, "VoxRaytracerRayDataElement")
+  py::class_<QuadMesh>(m, "QuadMesh")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def_readwrite("vox_id", &VoxRaytracer::RayData::Element::vox_id)
+      .def("AddPoint", &QuadMesh::AddPoint)
-        .def_readwrite("L", &VoxRaytracer::RayData::Element::L);
+      .def("AddQuad",
+           py::overload_cast<const Vector4i &>(&QuadMesh::AddQuad))
+      .def("Points", &QuadMesh::Points,
+           py::return_value_policy::reference_internal)
+      .def("Quads", &QuadMesh::Quads,
+           py::return_value_policy::reference_internal)
+      .def("GetQuad", &QuadMesh::GetQuad)
+      .def("GetNormal", &QuadMesh::GetNormal);
 
-    py::class_<VoxRaytracer::RayData>(m, "VoxRaytracerRayData")
+  py::class_<VoxRaytracer::RayData::Element>(m, "VoxRaytracerRayDataElement")
-        .def(py::init<>())
+      .def(py::init<>())
-        .def("AppendRay", &VoxRaytracer::RayData::AppendRay)
+      .def_readwrite("vox_id", &VoxRaytracer::RayData::Element::vox_id)
-        .def("Data", py::overload_cast<>(&VoxRaytracer::RayData::Data), py::return_value_policy::reference_internal)
+      .def_readwrite("L", &VoxRaytracer::RayData::Element::L);
-        .def("Count", &VoxRaytracer::RayData::Count)
-        .def("TotalLength", &VoxRaytracer::RayData::TotalLength)
-        .def("SetCount", &VoxRaytracer::RayData::SetCount)
-        .def("SetTotalLength", &VoxRaytracer::RayData::SetTotalLength);
 
-    py::class_<VoxRaytracer>(m, "VoxRaytracer")
+  py::class_<VoxRaytracer::RayData>(m, "VoxRaytracerRayData")
-        .def(py::init<StructuredGrid &>(), py::keep_alive<1, 2>())
+      .def(py::init<>())
-        .def("GetImage", &VoxRaytracer::GetImage, py::return_value_policy::reference_internal)
+      .def("AppendRay", &VoxRaytracer::RayData::AppendRay)
-        .def("TraceLine", &VoxRaytracer::TraceLine)
+      .def("Data", py::overload_cast<>(&VoxRaytracer::RayData::Data),
-        .def("TraceBetweenPoints", &VoxRaytracer::TraceBetweenPoints);
+           py::return_value_policy::reference_internal)
+      .def("Count", &VoxRaytracer::RayData::Count)
+      .def("TotalLength", &VoxRaytracer::RayData::TotalLength)
+      .def("SetCount", &VoxRaytracer::RayData::SetCount)
+      .def("SetTotalLength", &VoxRaytracer::RayData::SetTotalLength);
 
+  py::class_<VoxRaytracer>(m, "VoxRaytracer")
+      .def(py::init<StructuredGrid &>(), py::keep_alive<1, 2>())
+      .def("GetImage", &VoxRaytracer::GetImage,
+           py::return_value_policy::reference_internal)
+      .def("TraceLine", &VoxRaytracer::TraceLine)
+      .def("TraceBetweenPoints", &VoxRaytracer::TraceBetweenPoints);
 }

src/Python/testing/math_pybind_test.py

@@ -155,6 +155,18 @@ class TestMathNewTypes(unittest.TestCase):
         # SetValue(id, value) sets At(id).Value = value
         self.assertAlmostEqual(img.GetValue([0, 0, 0]), 10.5)
 
+    def test_quad_mesh(self):
+        mesh = uLib.Math.QuadMesh()
+        mesh.AddPoint(uLib.Math.Vector3f([0, 0, 0]))
+        mesh.AddPoint(uLib.Math.Vector3f([1, 0, 0]))
+        mesh.AddPoint(uLib.Math.Vector3f([1, 1, 0]))
+        mesh.AddPoint(uLib.Math.Vector3f([0, 1, 0]))
+        
+        mesh.AddQuad(uLib.Math.Vector4i([0, 1, 2, 3]))
+        
+        self.assertEqual(len(mesh.Points()), 4)
+        self.assertEqual(len(mesh.Quads()), 1)
+
 class TestMathVoxRaytracer(unittest.TestCase):
     def test_raytracer(self):
         grid = uLib.Math.StructuredGrid([10, 10, 10])

src/Python/uLib/__init__.py

@@ -1,7 +1,10 @@
 try:
     from .uLib_python import Core, Math
 except ImportError:
-    # Handle cases where the binary extension is not yet built
+    try:
-    pass
+        from uLib_python import Core, Math
+    except ImportError:
+        # Handle cases where the binary extension is not yet built
+        pass
 
 __all__ = ["Core", "Math"]

src/Root/CMakeLists.txt

@@ -56,7 +56,7 @@ set_target_properties(${libname} PROPERTIES
 target_link_libraries(${libname} ${LIBRARIES})
 
 install(TARGETS ${libname}
-        EXPORT "${PROJECT_NAME}Targets"
+        EXPORT "uLibTargets"
         RUNTIME DESTINATION ${INSTALL_BIN_DIR} COMPONENT bin
         LIBRARY DESTINATION ${INSTALL_LIB_DIR} COMPONENT lib)
 

src/Vtk/CMakeLists.txt

@@ -1,23 +1,47 @@
 set(HEADERS uLibVtkInterface.h
             uLibVtkViewer.h
             vtkContainerBox.h
-            vtkMuonScatter.h
+            vtkHandlerWidget.h
-            vtkStructuredGrid.h
+            vtkQViewport.h
-            vtkVoxRaytracerRepresentation.h
+            vtkViewport.h
-            vtkVoxImage.h)
+            vtkPolydata.h
+            )
 
 set(SOURCES uLibVtkInterface.cxx
             uLibVtkViewer.cpp
             vtkContainerBox.cpp
-            vtkMuonScatter.cxx
+            vtkHandlerWidget.cpp
-            vtkStructuredGrid.cpp
+            vtkQViewport.cpp
-            vtkVoxRaytracerRepresentation.cpp
+            vtkViewport.cpp
-            vtkVoxImage.cpp)
+            vtkPolydata.cpp
+            )
+
+## Pull in Math VTK wrappers (sets MATH_SOURCES / MATH_HEADERS)
+add_subdirectory(Math)
+list(APPEND SOURCES ${MATH_SOURCES})
+list(APPEND HEADERS ${MATH_HEADERS})
+
+## Pull in HEP/Detectors VTK wrappers (sets HEP_DETECTORS_SOURCES / HEADERS)
+add_subdirectory(HEP/Detectors)
+list(APPEND SOURCES ${HEP_DETECTORS_SOURCES})
+list(APPEND HEADERS ${HEP_DETECTORS_HEADERS})
+
+## Pull in HEP/MuonTomography VTK wrappers (sets HEP_MUONTOMOGRAPHY_SOURCES / HEADERS)
+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}Math
+              ${PACKAGE_LIBPREFIX}Detectors
+              ${PACKAGE_LIBPREFIX}Geant)
 
 if(USE_CUDA)
   find_package(CUDAToolkit REQUIRED)
@@ -30,16 +54,22 @@ set(ULIB_SELECTED_MODULES ${ULIB_SELECTED_MODULES} Vtk PARENT_SCOPE)
 
 add_library(${libname} SHARED ${SOURCES})
 set_target_properties(${libname} PROPERTIES
-                      VERSION ${PROJECT_VERSION}
+                                             SOVERSION ${PROJECT_SOVERSION}
-                      SOVERSION ${PROJECT_SOVERSION})
+                       AUTOMOC ON
-target_link_libraries(${libname} ${LIBRARIES})
+                       AUTOUIC ON
+                       AUTORCC ON)
+target_link_libraries(${libname} ${LIBRARIES} Qt6::Widgets)
 
 install(TARGETS ${libname}
-        EXPORT "${PROJECT_NAME}Targets"
+        EXPORT "uLibTargets"
         RUNTIME DESTINATION ${INSTALL_BIN_DIR} COMPONENT bin
         LIBRARY DESTINATION ${INSTALL_LIB_DIR} COMPONENT lib)
 
 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/CMakeLists.txt

@@ -0,0 +1,24 @@
+
+################################################################################
+##### Vtk/HEP/Detectors - VTK wrappers for HEP Detectors objects ##############
+################################################################################
+
+## Sources and headers are exported to parent scope so they get compiled
+## into the single mutomVtk shared library.
+
+set(HEP_DETECTORS_SOURCES
+    ${CMAKE_CURRENT_SOURCE_DIR}/vtkMuonScatter.cxx
+    ${CMAKE_CURRENT_SOURCE_DIR}/vtkMuonEvent.cxx
+    ${CMAKE_CURRENT_SOURCE_DIR}/vtkDetectorChamber.cxx
+    PARENT_SCOPE)
+
+set(HEP_DETECTORS_HEADERS
+    ${CMAKE_CURRENT_SOURCE_DIR}/vtkMuonScatter.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/vtkMuonEvent.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/vtkDetectorChamber.h
+    PARENT_SCOPE)
+
+if(BUILD_TESTING)
+  include(uLibTargetMacros)
+  add_subdirectory(testing)
+endif()

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

@@ -0,0 +1,16 @@
+# TESTS
+set(TESTS
+    vtkMuonScatterTest
+    vtkDetectorChamberTest
+)
+
+set(LIBRARIES
+    ${PACKAGE_LIBPREFIX}Core
+    ${PACKAGE_LIBPREFIX}Detectors
+    ${PACKAGE_LIBPREFIX}Vtk
+    ${VTK_LIBRARIES}
+    Boost::unit_test_framework
+)
+
+uLib_add_tests(VtkDetectors)
+

src/Vtk/HEP/Detectors/testing/vtkDetectorChamberTest.cpp

@@ -0,0 +1,68 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 "Vtk/HEP/Detectors/vtkDetectorChamber.h"
+#include "HEP/Detectors/DetectorChamber.h"
+#include "Math/Units.h"
+
+#include "Vtk/uLibVtkViewer.h"
+
+#define BOOST_TEST_MODULE vtkDetectorChamberTest
+#include <boost/test/unit_test.hpp>
+
+using namespace uLib;
+
+BOOST_AUTO_TEST_CASE(vtkDetectorChamberTest) {
+  DetectorChamber d1, d2;
+  // d1.SetSize(Vector3f(1, 1, 1));
+  // d1.SetPosition(Vector3f(0, 0, 0));
+  d1.Scale(Vector3f(5_m, 10_m, 2_m));
+  d1.Translate(Vector3f(0, 0, 0));
+
+  // d2.SetSize(Vector3f(1, 1, 1));
+  // d2.SetPosition(Vector3f(0, 0, 0));
+  d2.Scale(Vector3f(5_m, 10_m, 2_m));
+  d2.Translate(Vector3f(0, 0, 10_m));
+
+
+  Vtk::vtkDetectorChamber v_d1(&d1);
+  Vtk::vtkDetectorChamber v_d2(&d2);
+  v_d1.SetRepresentation(Vtk::Puppet::Surface);
+  v_d2.SetRepresentation(Vtk::Puppet::Surface);
+
+  if (!v_d1.GetProp()) {
+    BOOST_FAIL("vtkDetectorChamber::GetProp() returned NULL");
+  }
+
+  if (std::getenv("CTEST_PROJECT_NAME") == nullptr) {
+    Vtk::Viewer viewer;
+    viewer.SetGridAxis(Vtk::Viewport::Y);
+    viewer.AddPuppet(v_d1);
+    viewer.AddPuppet(v_d2);
+    viewer.Start();
+  }
+
+  BOOST_CHECK(true); // reached here without crash
+}

src/Vtk/HEP/Detectors/testing/vtkMuonScatterTest.cpp

@@ -0,0 +1,62 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 "Vtk/HEP/Detectors/vtkMuonScatter.h"
+#include "HEP/Detectors/MuonScatter.h"
+
+#include "Vtk/uLibVtkViewer.h"
+
+#define BOOST_TEST_MODULE VtkMuonScatterTest
+#include <boost/test/unit_test.hpp>
+
+using namespace uLib;
+
+BOOST_AUTO_TEST_CASE(vtkMuonScatterTest) {
+  MuonScatter event;
+  event.LineIn().direction << 0, -1, 1, 0;
+  event.LineIn().origin << 0, 1, -1, 1;
+
+  event.LineOut().direction << 0, -1, 0, 0;
+  event.LineOut().origin << 0, -1, 0, 1;
+
+  Vtk::vtkMuonScatter v_event(event);
+  v_event.AddPocaPoint(HPoint3f(0, 0, 0));
+
+  v_event.SaveToXMLFile("vtk_testing_muonevent.vtp");
+
+  if (std::getenv("CTEST_PROJECT_NAME") == nullptr) {
+    Vtk::Viewer viewer;
+
+    //    Vtk::Tie<Vtk::vtkMuonScatter> tms;
+    //    tms.DoAction();
+    //    Vtk::Tie<Vtk::Viewer> vms;
+    //    vms.DoAction();
+
+    viewer.AddPuppet(v_event);
+    viewer.Start();
+  }
+
+  BOOST_CHECK(true); // reached here without crash
+}

src/Vtk/HEP/Detectors/vtkDetectorChamber.cxx

@@ -0,0 +1,60 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 <vtkAbstractTransform.h>
+#include <vtkAxes.h>
+#include <vtkCubeSource.h>
+#include <vtkLineSource.h>
+#include <vtkMatrix4x4.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkPropPicker.h>
+#include <vtkProperty.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkRenderer.h>
+#include <vtkRendererCollection.h>
+#include <vtkSmartPointer.h>
+#include <vtkTransform.h>
+
+#include "Vtk/HEP/Detectors/vtkDetectorChamber.h"
+#include <vtkBoxWidget.h>
+#include <vtkTransformPolyDataFilter.h>
+
+namespace uLib {
+namespace Vtk {
+
+vtkDetectorChamber::vtkDetectorChamber(DetectorChamber *content)
+    : vtkContainerBox(content) {
+}
+
+vtkDetectorChamber::~vtkDetectorChamber() {
+}
+
+DetectorChamber *vtkDetectorChamber::GetContent() {
+  return static_cast<DetectorChamber *>(m_Content);
+}
+
+} // namespace Vtk
+} // namespace uLib

src/Vtk/HEP/Detectors/vtkDetectorChamber.h

@@ -0,0 +1,61 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 VTK_DETECTOR_CHAMBER_H
+#define VTK_DETECTOR_CHAMBER_H
+
+#include <vtkActor.h>
+#include <vtkCommand.h>
+#include <vtkSmartPointer.h>
+#include <vtkTransform.h>
+
+#include "HEP/Detectors/DetectorChamber.h"
+#include "Math/Dense.h"
+#include "Vtk/uLibVtkInterface.h"
+#include "Vtk/vtkContainerBox.h"
+#include <vtkActor.h>
+#include <vtkBoxWidget.h>
+#include <vtkTransformPolyDataFilter.h>
+
+namespace uLib {
+namespace Vtk {
+
+// class vtkHandlerWidget; // Removed as we use vtkBoxWidget now
+
+class vtkDetectorChamber : public vtkContainerBox {
+  typedef DetectorChamber Content;
+  typedef vtkContainerBox BaseClass;
+
+public:
+  vtkDetectorChamber(DetectorChamber *content);
+  virtual ~vtkDetectorChamber();
+
+  Content *GetContent();
+};
+
+} // namespace Vtk
+} // namespace uLib
+
+#endif // VTK_DETECTOR_CHAMBER_H

src/Vtk/HEP/Detectors/vtkMuonEvent.cxx

@@ -0,0 +1,145 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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.
+
+//////////////////////////////////////////////////////////////////////////////*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "Math/Dense.h"
+#include "Vtk/HEP/Detectors/vtkMuonEvent.h"
+
+namespace uLib {
+namespace Vtk {
+
+///// CALLBACK /////////////////////////////////////////////////////////////////
+
+class vtkWidgetCallback : public vtkCommand {
+public:
+  static vtkWidgetCallback *New() { return new vtkWidgetCallback; }
+  void SetParent(uLib::Vtk::vtkMuonEvent *parent) { this->parent = parent; }
+  virtual void Execute(vtkObject *caller, unsigned long, void *) {
+    vtkSmartPointer<vtkTransform> t = vtkSmartPointer<vtkTransform>::New();
+    vtkBoxWidget *widget = reinterpret_cast<vtkBoxWidget *>(caller);
+    widget->GetTransform(t);
+    // parent->SetTransform(t);
+    // std::cout << "event\n";
+  }
+
+private:
+  uLib::Vtk::vtkMuonEvent *parent;
+};
+
+////////////////////////////////////////////////////////////////////////////////
+/////  VTK MUON EVENT  /////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+vtkMuonEvent::vtkMuonEvent(MuonEventData &content)
+    : m_PolyData(NULL), m_Appender(vtkAppendPolyData::New()),
+      content(&content) {
+  InstallPipe();
+}
+
+vtkMuonEvent::vtkMuonEvent(const MuonEventData &content)
+    : m_PolyData(NULL), m_Appender(vtkAppendPolyData::New()),
+      content(const_cast<MuonEventData *>(&content)) {
+  InstallPipe();
+}
+
+vtkMuonEvent::~vtkMuonEvent() {}
+
+vtkMuonEvent::Content &vtkMuonEvent::GetContent() { return *content; }
+
+void vtkMuonEvent::PrintSelf(std::ostream &o) const {
+  o << "..:: MuonEvent ::..\n"
+       "\t[in]  Origin    > "
+    << content->LineIn().origin.transpose() << "\n"
+    << "\t[in]  Direction > " << content->LineIn().direction.transpose() << "\n"
+    << "\t[out] Origin    > " << content->LineOut().origin.transpose() << "\n"
+    << "\t[out] Direction > " << content->LineOut().direction.transpose()
+    << "\n"
+    << "\tMomentum        > " << content->GetMomentum() << "\n"
+    << "...................\n";
+}
+
+void vtkMuonEvent::InstallPipe() {
+
+  vtkAppendPolyData *appender = m_Appender;
+  vtkSmartPointer<vtkPolyDataMapper> mapper =
+      vtkSmartPointer<vtkPolyDataMapper>::New();
+
+  if (content) {
+    vtkSmartPointer<vtkLineSource> line_in =
+        vtkSmartPointer<vtkLineSource>::New();
+    vtkSmartPointer<vtkLineSource> line_out =
+        vtkSmartPointer<vtkLineSource>::New();
+
+    float distance =
+        (content->LineIn().origin - content->LineOut().origin).norm() / 10;
+
+    HPoint3f pt;
+    pt = content->LineIn().origin;
+    line_in->SetPoint1(pt(0), pt(1), pt(2));
+    pt = content->LineIn().origin + content->LineIn().direction * distance;
+    line_in->SetPoint2(pt(0), pt(1), pt(2));
+    pt = content->LineOut().origin;
+    line_out->SetPoint1(pt(0), pt(1), pt(2));
+    pt = content->LineOut().origin + content->LineOut().direction * distance;
+    line_out->SetPoint2(pt(0), pt(1), pt(2));
+
+    appender->AddInputConnection(line_in->GetOutputPort());
+    appender->AddInputConnection(line_out->GetOutputPort());
+  }
+
+  appender->Update();
+  mapper->SetInputConnection(appender->GetOutputPort());
+  mapper->Update();
+
+  vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
+  actor->SetMapper(mapper);
+  this->SetProp(actor);
+}
+
+vtkPolyData *vtkMuonEvent::GetPolyData() const {
+  return m_Appender->GetOutput();
+}
+
+void vtkMuonEvent::AddPocaPoint(HPoint3f poca) {
+  m_Poca = poca;
+  vtkSmartPointer<vtkSphereSource> sphere =
+      vtkSmartPointer<vtkSphereSource>::New();
+  float size =
+      (content->LineIn().origin - content->LineOut().origin).head(3).norm();
+  size /= 100;
+  sphere->SetRadius(size);
+  sphere->SetCenter(poca(0), poca(1), poca(2));
+  sphere->Update();
+  m_Appender->AddInputConnection(sphere->GetOutputPort());
+  m_Appender->Update();
+}
+
+HPoint3f vtkMuonEvent::GetPocaPoint() { return m_Poca; }
+
+} // namespace Vtk
+} // namespace uLib

src/Vtk/HEP/Detectors/vtkMuonEvent.h

@@ -23,17 +23,15 @@
 
 //////////////////////////////////////////////////////////////////////////////*/
 
-
-
 #ifndef VTKMUONEVENT_H
 #define VTKMUONEVENT_H
 
+#include <vtkActor.h>
+#include <vtkAppendPolyData.h>
+#include <vtkLineSource.h>
+#include <vtkPolyDataMapper.h>
 #include <vtkSmartPointer.h>
 #include <vtkSphereSource.h>
-#include <vtkAppendPolyData.h>
-#include <vtkPolyDataMapper.h>
-#include <vtkLineSource.h>
-#include <vtkActor.h>
 
 #include <vtk3DWidget.h>
 #include <vtkBoxWidget.h>
@@ -45,50 +43,46 @@
 
 #include "Math/Dense.h"
 
-#include "uLibVtkInterface.h"
+#include "HEP/Detectors/MuonEvent.h"
-#include "Detectors/MuonEvent.h"
+#include "Vtk/uLibVtkInterface.h"
-
+#include "Vtk/vtkPolydata.h"
-
-
 
 namespace uLib {
+namespace Vtk {
+
+class vtkMuonEvent : public Puppet, public Polydata {
+  typedef MuonEventData Content;
 
-class vtkMuonEvent : public Abstract::uLibVtkPolydata {
-    typedef MuonEventData Content;
 public:
-    vtkMuonEvent(const MuonEventData &content);
+  vtkMuonEvent(const MuonEventData &content);
-    vtkMuonEvent(MuonEventData &content);
+  vtkMuonEvent(MuonEventData &content);
-    ~vtkMuonEvent();
+  ~vtkMuonEvent();
 
-    Content& GetContent();
+  Content &GetContent();
 
-    void PrintSelf(std::ostream &o) const;
+  void PrintSelf(std::ostream &o) const;
 
-    virtual vtkProp *GetProp();
+  virtual vtkPolyData *GetPolyData() const;
 
-    virtual vtkPolyData* GetPolyData() const;
+  void AddPocaPoint(HPoint3f poca);
 
-    void AddPocaPoint(HPoint3f poca);
+  HPoint3f GetPocaPoint();
-
-    HPoint3f GetPocaPoint();
-
-    void vtkStartInteractive();
 
+  void vtkStartInteractive();
 
 private:
-    void InstallPipe();
+  void InstallPipe();
 
-    vtkMuonEvent::Content *content;
+  vtkMuonEvent::Content *content;
-    vtkPolyData           *m_PolyData;
+  vtkPolyData *m_PolyData;
-    vtkActor              *m_Prop;
 
-    vtkAppendPolyData     *m_Appender;
+  vtkAppendPolyData *m_Appender;
-    vtkBoxWidget          *m_WidgetIN;
+  vtkBoxWidget *m_WidgetIN;
-    vtkBoxWidget          *m_WidgetOUT;
+  vtkBoxWidget *m_WidgetOUT;
-    HPoint3f               m_Poca;
+  HPoint3f m_Poca;
 };
 
+} // namespace Vtk
+} // namespace uLib
 
-}
+#endif // VTKMUONEVENT_H
-
-#endif // VTKMUONSCATTER_H

src/Vtk/HEP/Detectors/vtkMuonScatter.cxx

@@ -0,0 +1,144 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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.
+
+//////////////////////////////////////////////////////////////////////////////*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "Math/Dense.h"
+#include "Vtk/HEP/Detectors/vtkMuonScatter.h"
+
+namespace uLib {
+namespace Vtk {
+
+////////////////////////////////////////////////////////////////////////////////
+/////  VTK MUON Scatter  ///////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+vtkMuonScatter::vtkMuonScatter(MuonScatter &content)
+    : m_Content(&content), m_LineIn(vtkLineSource::New()),
+      m_LineOut(vtkLineSource::New()), m_PolyData(vtkPolyData::New()),
+      m_SpherePoca(NULL) {
+  InstallPipe();
+}
+
+vtkMuonScatter::vtkMuonScatter(const MuonScatter &content)
+    : m_Content(const_cast<MuonScatter *>(&content)),
+      m_LineIn(vtkLineSource::New()), m_LineOut(vtkLineSource::New()),
+      m_PolyData(vtkPolyData::New()), m_SpherePoca(NULL) {
+  InstallPipe();
+}
+
+vtkMuonScatter::~vtkMuonScatter() {
+  m_LineIn->Delete();
+  m_LineOut->Delete();
+  if (m_SpherePoca)
+    m_SpherePoca->Delete();
+}
+
+vtkMuonScatter::Content &vtkMuonScatter::GetContent() { return *m_Content; }
+
+void vtkMuonScatter::PrintSelf(std::ostream &o) const {}
+
+void vtkMuonScatter::InstallPipe() {
+  if (m_Content) {
+    vtkLineSource *line_in = m_LineIn;
+    vtkLineSource *line_out = m_LineOut;
+
+    float distance =
+        (m_Content->LineIn().origin - m_Content->LineOut().origin).norm() / 10;
+
+    HPoint3f pt;
+    pt = m_Content->LineIn().origin;
+    line_in->SetPoint1(pt(0), pt(1), pt(2));
+    pt = m_Content->LineIn().origin + m_Content->LineIn().direction * distance;
+    line_in->SetPoint2(pt(0), pt(1), pt(2));
+    pt = m_Content->LineOut().origin;
+    line_out->SetPoint1(pt(0), pt(1), pt(2));
+    pt =
+        m_Content->LineOut().origin + m_Content->LineOut().direction * distance;
+    line_out->SetPoint2(pt(0), pt(1), pt(2));
+  }
+
+  vtkSmartPointer<vtkPolyDataMapper> mapper =
+      vtkSmartPointer<vtkPolyDataMapper>::New();
+  mapper->SetInputConnection(m_LineIn->GetOutputPort());
+  vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
+  actor->SetMapper(mapper);
+  this->SetProp(actor);
+
+  mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
+  mapper->SetInputConnection(m_LineOut->GetOutputPort());
+  actor = vtkSmartPointer<vtkActor>::New();
+  actor->SetMapper(mapper);
+  this->SetProp(actor);
+}
+
+vtkPolyData *vtkMuonScatter::GetPolyData() const {
+  vtkSmartPointer<vtkAppendPolyData> append =
+      vtkSmartPointer<vtkAppendPolyData>::New();
+  append->AddInputConnection(m_LineIn->GetOutputPort());
+  append->AddInputConnection(m_LineOut->GetOutputPort());
+  if (m_SpherePoca)
+    append->AddInputConnection(m_SpherePoca->GetOutputPort());
+  append->Update();
+  m_PolyData->DeepCopy(append->GetOutput());
+  return m_PolyData;
+}
+
+void vtkMuonScatter::AddPocaPoint(HPoint3f poca) {
+  vtkSphereSource *sphere = vtkSphereSource::New();
+  float size =
+      (m_Content->LineIn().origin - m_Content->LineOut().origin).head(3).norm();
+  size /= 100;
+  sphere->SetRadius(size);
+  sphere->SetCenter(poca(0), poca(1), poca(2));
+  sphere->Update();
+  m_SpherePoca = sphere;
+
+  vtkSmartPointer<vtkPolyDataMapper> mapper =
+      vtkSmartPointer<vtkPolyDataMapper>::New();
+  mapper->SetInputConnection(m_SpherePoca->GetOutputPort());
+  vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
+  actor->SetMapper(mapper);
+  this->SetProp(actor);
+}
+
+HPoint3f vtkMuonScatter::GetPocaPoint() {
+  double center[3];
+  if (m_SpherePoca) {
+    m_SpherePoca->GetCenter(center);
+    return HPoint3f(center[0], center[1], center[2]);
+  } else {
+    return HPoint3f(0, 0, 0);
+  }
+}
+
+void vtkMuonScatter::ConnectInteractor(vtkRenderWindowInteractor *interactor) {
+  // TODO
+}
+
+} // namespace Vtk
+} // namespace uLib

src/Vtk/HEP/Detectors/vtkMuonScatter.h

@@ -23,18 +23,16 @@
 
 //////////////////////////////////////////////////////////////////////////////*/
 
-
-
 #ifndef VTKMUONSCATTER_H
 #define VTKMUONSCATTER_H
 
-#include <vtkVersion.h>
+#include <vtkActor.h>
+#include <vtkAppendPolyData.h>
+#include <vtkLineSource.h>
+#include <vtkPolyDataMapper.h>
 #include <vtkSmartPointer.h>
 #include <vtkSphereSource.h>
-#include <vtkAppendPolyData.h>
+#include <vtkVersion.h>
-#include <vtkPolyDataMapper.h>
-#include <vtkLineSource.h>
-#include <vtkActor.h>
 
 #include <vtk3DWidget.h>
 #include <vtkBoxWidget.h>
@@ -46,8 +44,9 @@
 
 #include "Math/Dense.h"
 
-#include "uLibVtkInterface.h"
+#include "HEP/Detectors/MuonScatter.h"
-#include "Detectors/MuonScatter.h"
+#include "Vtk/uLibVtkInterface.h"
+#include "Vtk/vtkPolydata.h"
 
 class vtkRenderWindowInteractor;
 
@@ -55,40 +54,40 @@ namespace uLib {
 namespace Vtk {
 
 class vtkMuonScatter : public Puppet, public Polydata {
-    typedef MuonScatter Content;
+  typedef MuonScatter Content;
+
 public:
+  vtkMuonScatter(const MuonScatter &content);
+  vtkMuonScatter(MuonScatter &content);
 
-    vtkMuonScatter(const MuonScatter &content);
+  ~vtkMuonScatter();
-    vtkMuonScatter(MuonScatter &content);
 
-    ~vtkMuonScatter();
+  Content &GetContent();
 
-    Content &GetContent();
+  void PrintSelf(std::ostream &o) const;
 
-    void PrintSelf(std::ostream &o) const;
+  virtual vtkPolyData *GetPolyData() const;
 
-    virtual vtkPolyData* GetPolyData() const;
+  void AddPocaPoint(HPoint3f poca);
 
-    void AddPocaPoint(HPoint3f poca);
+  HPoint3f GetPocaPoint();
 
-    HPoint3f GetPocaPoint();
+  void vtkStartInteractive();
-
-    void vtkStartInteractive();
 
 protected:
-    void ConnectInteractor(vtkRenderWindowInteractor *interactor);
+  void ConnectInteractor(vtkRenderWindowInteractor *interactor);
 
 private:
-    void InstallPipe();
+  void InstallPipe();
 
-    vtkMuonScatter::Content *m_Content;
+  vtkMuonScatter::Content *m_Content;
-    vtkLineSource           *m_LineIn;
+  vtkLineSource *m_LineIn;
-    vtkLineSource           *m_LineOut;
+  vtkLineSource *m_LineOut;
-    vtkSphereSource         *m_SpherePoca;
+  vtkSphereSource *m_SpherePoca;
-    vtkPolyData             *m_PolyData;
+  vtkPolyData *m_PolyData;
 };
 
-} // vtk
+} // namespace Vtk
-} // uLib
+} // namespace uLib
 
 #endif // VTKMUONSCATTER_H

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/HEP/MuonTomography/CMakeLists.txt

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