add Geant namespace

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

@@ -159,19 +159,26 @@ endif()
 
 find_package(Qt6 COMPONENTS Widgets REQUIRED)
 
-find_package(Geant4 REQUIRED)
-message(STATUS "Geant4 libs: ${Geant4_LIBRARIES}")
+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"
-    )
+    # 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)

app/CMakeLists.txt

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

app/gcompose/src/main.cpp

@@ -41,7 +41,7 @@ int main(int argc, char** argv) {
     d2.Translate(Vector3f(0, 0, 10));
     
 
-    Scene scene;
+    Geant::Scene scene;
     scene.ConstructWorldBox(&world_box, "G4_AIR");
     scene.Initialize();
 

src/HEP/Geant/ActionInitialization.cpp

@@ -1,6 +1,9 @@
 #include "ActionInitialization.hh"
 #include "EmitterPrimary.hh"
 
+namespace uLib {
+namespace Geant {
+
 ActionInitialization::ActionInitialization() : G4VUserActionInitialization() {}
 
 ActionInitialization::~ActionInitialization() {}
@@ -22,4 +25,7 @@ void ActionInitialization::Build() const {
   // SetUserAction(new RunAction());
   // SetUserAction(new EventAction());
   // SetUserAction(new SteppingAction());
-}
+}
+
+} // namespace Geant
+} // namespace uLib

src/HEP/Geant/ActionInitialization.hh

@@ -3,6 +3,9 @@
 
 #include "G4VUserActionInitialization.hh"
 
+namespace uLib {
+namespace Geant {
+
 class ActionInitialization : public G4VUserActionInitialization {
 public:
   ActionInitialization();
@@ -15,4 +18,7 @@ public:
   virtual void Build() const;
 };
 
+} // namespace Geant
+} // namespace uLib
+
 #endif

src/HEP/Geant/DetectorConstruction.cpp

@@ -12,6 +12,7 @@
 #include "G4SystemOfUnits.hh"
 
 namespace uLib {
+namespace Geant {
 
 DetectorConstruction::DetectorConstruction(const char *name) : G4VUserDetectorConstruction() {}
 
@@ -21,9 +22,5 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { return nullptr; }
 
 void DetectorConstruction::ConstructSDandField() {}
 
-
-
-
-
-
-}
+} // namespace Geant
+} // namespace uLib

src/HEP/Geant/DetectorConstruction.hh

@@ -9,6 +9,7 @@ class G4VPhysicalVolume;
 class G4LogicalVolume;
 
 namespace uLib {
+namespace Geant {
 
 class DetectorConstruction : public G4VUserDetectorConstruction {
 public:
@@ -20,6 +21,7 @@ public:
   virtual void ConstructSDandField();
 };
 
+} // namespace Geant
 } // namespace uLib
 
 #endif

src/HEP/Geant/EmitterPrimary.cpp

@@ -10,6 +10,9 @@
 #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
@@ -47,4 +50,106 @@ void EmitterPrimary::GeneratePrimaries(G4Event *anEvent) {
   // 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

@@ -4,9 +4,15 @@
 #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:
@@ -16,8 +22,33 @@ class EmitterPrimary : public G4VUserPrimaryGeneratorAction
     // Metodo principale chiamato all'inizio di ogni evento
     virtual void GeneratePrimaries(G4Event*);         
 
-  private:
+  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.h

@@ -33,6 +33,7 @@
 #include "Detectors/ChamberHitEvent.h"
 
 namespace uLib {
+namespace Geant {
 
 class GeantEventData {
 public:
@@ -56,6 +57,7 @@ public:
   uLibRefMacro(ChEvents, Vector<ChamberHitEventData>)
 };
 
+} // namespace Geant
 } // namespace uLib
 
 #endif // 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

@@ -4,6 +4,9 @@
 #include "G4EmStandardPhysics.hh"
 #include "G4RadioactiveDecayPhysics.hh"
 
+namespace uLib {
+namespace Geant {
+
 PhysicsList::PhysicsList() : G4VModularPhysicsList() {
   SetVerboseLevel(1);
 
@@ -20,3 +23,6 @@ PhysicsList::PhysicsList() : G4VModularPhysicsList() {
 PhysicsList::~PhysicsList() {}
 
 void PhysicsList::SetCuts() { G4VModularPhysicsList::SetCuts(); }
+
+} // namespace Geant
+} // namespace uLib

src/HEP/Geant/PhysicsList.hh

@@ -3,6 +3,9 @@
 
 #include "G4VModularPhysicsList.hh"
 
+namespace uLib {
+namespace Geant {
+
 class PhysicsList : public G4VModularPhysicsList {
 public:
   PhysicsList();
@@ -11,4 +14,7 @@ public:
   virtual void SetCuts();
 };
 
+} // namespace Geant
+} // namespace uLib
+
 #endif

src/HEP/Geant/Scene.cpp

@@ -44,6 +44,7 @@
 #include "ActionInitialization.hh"
 
 namespace uLib {
+namespace Geant {
 
 class SceneDetectorConstruction : public DetectorConstruction {
 public:
@@ -146,4 +147,5 @@ void Scene::Initialize() {
   d->Initialize();
 }
 
+} // namespace Geant
 } // namespace uLib

src/HEP/Geant/Scene.h

@@ -35,6 +35,7 @@
 class G4VPhysicalVolume;
 
 namespace uLib {
+namespace Geant {
 
 class Scene : public Object {
 public:
@@ -51,5 +52,6 @@ private:
   class SceneImpl *d;
 };
 
+} // namespace Geant
 } // namespace uLib
 #endif // SCENE_H

src/HEP/Geant/Solid.cpp

@@ -43,6 +43,7 @@
 #include "Solid.h"
 
 namespace uLib {
+namespace Geant {
 
 class DetectorsSolidImpl {
 public:
@@ -54,22 +55,23 @@ public:
 };
 
 Solid::Solid()
-    : m_Logical(new G4LogicalVolume(NULL, NULL, "unnamed_solid")),
-      m_Material(NULL) {}
+    : m_Name("unnamed_solid"), m_Material(NULL), m_Logical(NULL), m_Physical(NULL) {}
 
 Solid::Solid(const char *name)
-    : m_Logical(new G4LogicalVolume(NULL, NULL, name)), m_Material(NULL) {}
+    : m_Name(name), m_Material(NULL), m_Logical(NULL), m_Physical(NULL) {}
 
 void Solid::SetNistMaterial(const char *name) {
   G4NistManager *nist = G4NistManager::Instance();
   m_Material = nist->FindOrBuildMaterial(name);
-  m_Logical->SetMaterial(m_Material);
+  if (m_Logical)
+    m_Logical->SetMaterial(m_Material);
 }
 
 void Solid::SetMaterial(G4Material *material) {
   if (material) {
     m_Material = material;
-    m_Logical->SetMaterial(material);
+    if (m_Logical)
+      m_Logical->SetMaterial(material);
   }
 }
 
@@ -174,4 +176,5 @@ void BoxSolid::Update() {
 }
 
 
+} // namespace Geant
 } // namespace uLib

src/HEP/Geant/Solid.h

@@ -9,7 +9,7 @@
 
   ------------------------------------------------------------------
   This library is free software;  you  can  redistribute  it  and/or
-  modify it  under the  terms  of  the  GNU  Lesser  General  Public
+  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.
 
@@ -39,6 +39,7 @@ class G4TessellatedSolid;
 class G4Box;
 
 namespace uLib {
+namespace Geant {
 
 class Solid : public Object {
 public:
@@ -58,10 +59,13 @@ public:
   uLibGetSetMacro(Logical, G4LogicalVolume *)
   uLibGetSetMacro(Physical, G4VPhysicalVolume *)
 
-  inline const char *GetName() const { return m_Logical->GetName().c_str(); }
+  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
@@ -106,6 +110,7 @@ private:
   G4Box *m_Solid;
 };
 
+} // namespace Geant
 } // namespace uLib
 
 #endif // SOLID_H

src/HEP/Geant/testing/ActionInitialization.cpp

@@ -13,7 +13,7 @@ int main(int argc, char **argv) {
   // runManager->SetUserInitialization(new PhysicsList());
 
   // 3. INIZIALIZZAZIONE DELLE AZIONI (Il nostro generatore!)
-  runManager->SetUserInitialization(new ActionInitialization());
+  runManager->SetUserInitialization(new uLib::Geant::ActionInitialization());
 
   // ... Inizializzazione del kernel ( runManager->Initialize(); ), UI manager,
   // vis manager, ecc.

src/HEP/Geant/testing/CMakeLists.txt

@@ -3,6 +3,7 @@ set(TESTS
        SolidTest
        EventTest
        GeantApp
+       ActionInitialization
 )
 
 set(LIBRARIES

src/HEP/Geant/testing/EventTest.cpp

@@ -15,7 +15,7 @@ int main() {
 
     // Test Solid initialization and NIST material //
     {
-        Solid solid("test_solid");
+        Geant::Solid solid("test_solid");
         TEST1(solid.GetLogical() != nullptr);
         
         solid.SetNistMaterial("G4_AIR");
@@ -25,7 +25,7 @@ int main() {
 
     // Test TessellatedSolid with a simple mesh //
     {
-        TessellatedSolid tsolid("test_tessellated");
+        Geant::TessellatedSolid tsolid("test_tessellated");
         TEST1(tsolid.GetLogical() != nullptr);
         TEST1(tsolid.GetSolid() != nullptr);
 

src/HEP/Geant/testing/GeantApp.cpp

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

src/HEP/Geant/testing/SolidTest.cpp

@@ -14,7 +14,7 @@ int main() {
 
     // Test Solid initialization and NIST material //
     {
-        Solid solid("test_solid");
+        Geant::Solid solid("test_solid");
         TEST1(solid.GetLogical() != nullptr);
         
         solid.SetNistMaterial("G4_AIR");
@@ -24,7 +24,7 @@ int main() {
 
     // Test TessellatedSolid with a simple mesh //
     {
-        TessellatedSolid tsolid("test_tessellated");
+        Geant::TessellatedSolid tsolid("test_tessellated");
         TEST1(tsolid.GetLogical() != nullptr);
         TEST1(tsolid.GetSolid() != nullptr);
 

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

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/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

@@ -2,7 +2,12 @@
 #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;

src/Math/testing/CMakeLists.txt

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

src/Math/testing/QuadMeshTest.cpp

@@ -0,0 +1,51 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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/QuadMesh.h"
+#include <iostream>
+
+using namespace uLib;
+
+int main() {
+  BEGIN_TESTING(Quad Mesh);
+
+  QuadMesh mesh;
+
+  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/testing-prototype.h

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

src/Python/math_bindings.cpp

@@ -15,6 +15,7 @@
 #include "Math/StructuredGrid.h"
 #include "Math/Transform.h"
 #include "Math/TriangleMesh.h"
+#include "Math/QuadMesh.h"
 #include "Math/VoxImage.h"
 #include "Math/VoxRaytracer.h"
 
@@ -266,7 +267,8 @@ void init_math(py::module_ &m) {
       .def("__call__", &Accumulator_ABClip<double>::operator());
 
   // 5. Core Math Structures
-  py::class_<AffineTransform>(m, "AffineTransform")
+  py::class_<AffineTransform, std::shared_ptr<AffineTransform>>(m,
+                                                              "AffineTransform")
       .def(py::init<>())
       .def("GetWorldMatrix", &AffineTransform::GetWorldMatrix)
       .def("SetPosition", &AffineTransform::SetPosition)
@@ -283,7 +285,7 @@ void init_math(py::module_ &m) {
       .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("GetWorldPoint", py::overload_cast<const Vector4f>(
                                 &Geometry::GetWorldPoint, py::const_))
@@ -294,7 +296,8 @@ void init_math(py::module_ &m) {
       .def("GetLocalPoint",
            py::overload_cast<float, float, float>(&Geometry::GetLocalPoint));
 
-  py::class_<ContainerBox, AffineTransform, Object>(m, "ContainerBox")
+  py::class_<ContainerBox, AffineTransform, Object, std::shared_ptr<ContainerBox>>(
+      m, "ContainerBox")
       .def(py::init<>())
       .def("SetOrigin", &ContainerBox::SetOrigin)
       .def("GetOrigin", &ContainerBox::GetOrigin)
@@ -322,7 +325,7 @@ void init_math(py::module_ &m) {
       .value("ZYX", StructuredData::ZYX)
       .export_values();
 
-  py::class_<StructuredData>(m, "StructuredData")
+  py::class_<StructuredData, std::shared_ptr<StructuredData>>(m, "StructuredData")
       .def(py::init<const Vector3i &>())
       .def("GetDims", &StructuredData::GetDims)
       .def("SetDims", &StructuredData::SetDims)
@@ -334,7 +337,8 @@ void init_math(py::module_ &m) {
       .def("Map", &StructuredData::Map)
       .def("UnMap", &StructuredData::UnMap);
 
-  py::class_<StructuredGrid, ContainerBox, StructuredData>(m, "StructuredGrid")
+  py::class_<StructuredGrid, ContainerBox, StructuredData,
+             std::shared_ptr<StructuredGrid>>(m, "StructuredGrid")
       .def(py::init<const Vector3i &>())
       .def("SetSpacing", &StructuredGrid::SetSpacing)
       .def("GetSpacing", &StructuredGrid::GetSpacing)
@@ -379,7 +383,8 @@ void init_math(py::module_ &m) {
       .def_readwrite("Value", &Voxel::Value)
       .def_readwrite("Count", &Voxel::Count);
 
-  py::class_<Abstract::VoxImage, StructuredGrid>(m, "AbstractVoxImage")
+  py::class_<Abstract::VoxImage, StructuredGrid,
+             std::shared_ptr<Abstract::VoxImage>>(m, "AbstractVoxImage")
       .def("GetValue", py::overload_cast<const Vector3i &>(
                            &Abstract::VoxImage::GetValue, py::const_))
       .def("GetValue", py::overload_cast<const int>(
@@ -393,7 +398,8 @@ void init_math(py::module_ &m) {
       .def("ImportFromVtk", &Abstract::VoxImage::ImportFromVtk)
       .def("ImportFromVti", &Abstract::VoxImage::ImportFromVti);
 
-  py::class_<VoxImage<Voxel>, Abstract::VoxImage>(m, "VoxImage")
+  py::class_<VoxImage<Voxel>, Abstract::VoxImage,
+             std::shared_ptr<VoxImage<Voxel>>>(m, "VoxImage")
       .def(py::init<>())
       .def(py::init<const Vector3i &>())
       .def("Data", &VoxImage<Voxel>::Data,
@@ -429,7 +435,21 @@ void init_math(py::module_ &m) {
       .def("Points", &TriangleMesh::Points,
            py::return_value_policy::reference_internal)
       .def("Triangles", &TriangleMesh::Triangles,
-           py::return_value_policy::reference_internal);
+           py::return_value_policy::reference_internal)
+      .def("GetTriangle", &TriangleMesh::GetTriangle)
+      .def("GetNormal", &TriangleMesh::GetNormal);
+
+  py::class_<QuadMesh>(m, "QuadMesh")
+      .def(py::init<>())
+      .def("AddPoint", &QuadMesh::AddPoint)
+      .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::Element>(m, "VoxRaytracerRayDataElement")
       .def(py::init<>())

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
-    pass
+    try:
+        from uLib_python import Core, Math
+    except ImportError:
+        # Handle cases where the binary extension is not yet built
+        pass
 
 __all__ = ["Core", "Math"]

src/Vtk/Math/CMakeLists.txt

@@ -6,6 +6,7 @@
 set(MATH_SOURCES
     ${CMAKE_CURRENT_SOURCE_DIR}/vtkStructuredGrid.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/vtkTriangleMesh.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/vtkQuadMesh.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/vtkVoxImage.cpp
     PARENT_SCOPE)
 
@@ -13,6 +14,7 @@ set(MATH_HEADERS
     ${CMAKE_CURRENT_SOURCE_DIR}/vtkHLineRepresentation.h
     ${CMAKE_CURRENT_SOURCE_DIR}/vtkStructuredGrid.h
     ${CMAKE_CURRENT_SOURCE_DIR}/vtkTriangleMesh.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/vtkQuadMesh.h
     ${CMAKE_CURRENT_SOURCE_DIR}/vtkVoxImage.h
     PARENT_SCOPE)
 

src/Vtk/Math/testing/CMakeLists.txt

@@ -2,6 +2,7 @@
 set(TESTS
     vtkStructuredGridTest
     vtkTriangleMeshTest
+    vtkQuadMeshTest
     vtkVoxImageTest
 )
 

src/Vtk/Math/testing/vtkQuadMeshTest.cpp

@@ -0,0 +1,56 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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/Math/vtkQuadMesh.h"
+#include "Math/QuadMesh.h"
+#include "Vtk/uLibVtkViewer.h"
+
+#define BOOST_TEST_MODULE VtkQuadMeshTest
+#include <boost/test/unit_test.hpp>
+
+using namespace uLib;
+
+BOOST_AUTO_TEST_CASE(vtkQuadMeshConstruction) {
+  QuadMesh mesh;
+
+  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));
+
+  Vtk::vtkQuadMesh v_mesh(mesh);
+  v_mesh.Update();
+
+  if (std::getenv("CTEST_PROJECT_NAME") == nullptr) {
+    Vtk::Viewer viewer;
+    viewer.AddPuppet(v_mesh);
+    viewer.Start();
+  }
+
+  BOOST_CHECK_EQUAL(mesh.Points().size(), 4u);
+  BOOST_CHECK_EQUAL(mesh.Quads().size(), 1u);
+}

src/Vtk/Math/testing/vtkTriangleMeshTest.cpp

@@ -66,6 +66,6 @@ BOOST_AUTO_TEST_CASE(vtkTriangleMeshConstruction2) {
     viewer.Start();
   }
 
-  BOOST_CHECK_EQUAL(mesh.Points().size(), 3u);
-  BOOST_CHECK_EQUAL(mesh.Triangles().size(), 1u);
+  BOOST_CHECK_EQUAL(mesh.Points().size(), 9430u);
+  BOOST_CHECK_EQUAL(mesh.Triangles().size(), 18753u);
 }

src/Vtk/Math/vtkQuadMesh.cpp

@@ -0,0 +1,171 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 <vtkSmartPointer.h>
+
+#include <vtkOBJReader.h>
+#include <vtkPolyDataReader.h>
+#include <vtkSTLReader.h>
+#include <vtkXMLPolyDataReader.h>
+
+#include <vtkActor.h>
+#include <vtkCellArray.h>
+#include <vtkPoints.h>
+#include <vtkPolyData.h>
+#include <vtkPolyDataMapper.h>
+
+#include "Vtk/Math/vtkQuadMesh.h"
+#include <iostream>
+
+namespace uLib {
+namespace Vtk {
+
+void vtkQuadMesh::vtk2uLib_update() {
+  vtkIdType number_of_points = m_Poly->GetNumberOfPoints();
+  vtkIdType number_of_quads = m_Poly->GetNumberOfPolys();
+
+  std::cout << "//////\n"
+            << "number of points = " << number_of_points << "\n"
+            << "number of quads  = " << number_of_quads << "\n"
+            << "//////\n";
+
+  m_content.Points().resize(number_of_points);
+  for (int i = 0; i < number_of_points; ++i) {
+    double *point = m_Poly->GetPoint(i);
+    m_content.Points()[i](0) = point[0];
+    m_content.Points()[i](1) = point[1];
+    m_content.Points()[i](2) = point[2];
+  }
+
+  m_content.Quads().resize(number_of_quads);
+  m_Poly->GetPolys()->InitTraversal();
+  vtkSmartPointer<vtkIdList> idList = vtkSmartPointer<vtkIdList>::New();
+  for (int i = 0; i < number_of_quads; ++i) {
+    m_Poly->GetPolys()->GetNextCell(idList);
+    if (idList->GetNumberOfIds() == 4) {
+      m_content.Quads()[i](0) = idList->GetId(0);
+      m_content.Quads()[i](1) = idList->GetId(1);
+      m_content.Quads()[i](2) = idList->GetId(2);
+      m_content.Quads()[i](3) = idList->GetId(3);
+    }
+  }
+  m_Poly->Modified();
+  m_Actor->GetMapper()->Update();
+}
+
+void vtkQuadMesh::uLib2vtk_update() {
+  vtkIdType number_of_points = m_content.Points().size();
+  vtkIdType number_of_quads = m_content.Quads().size();
+
+  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
+  points->SetNumberOfPoints(number_of_points);
+  for (vtkIdType i = 0; i < number_of_points; i++) {
+    double x, y, z;
+    x = m_content.Points().at(i)(0);
+    y = m_content.Points().at(i)(1);
+    z = m_content.Points().at(i)(2);
+    points->SetPoint(i, x, y, z);
+  }
+
+  vtkSmartPointer<vtkCellArray> polys = vtkSmartPointer<vtkCellArray>::New();
+  for (vtkIdType i = 0; i < number_of_quads; i++) {
+    vtkIdType a, b, c, d;
+    a = m_content.Quads().at(i)(0);
+    b = m_content.Quads().at(i)(1);
+    c = m_content.Quads().at(i)(2);
+    d = m_content.Quads().at(i)(3);
+    polys->InsertNextCell(4);
+    polys->InsertCellPoint(a);
+    polys->InsertCellPoint(b);
+    polys->InsertCellPoint(c);
+    polys->InsertCellPoint(d);
+  }
+
+  m_Poly->SetPoints(points);
+  m_Poly->SetPolys(polys);
+  m_Poly->Modified();
+  m_Actor->GetMapper()->Update();
+}
+
+// -------------------------------------------------------------------------- //
+
+vtkQuadMesh::vtkQuadMesh(vtkQuadMesh::Content &content)
+    : m_content(content), m_Poly(vtkPolyData::New()), m_Actor(vtkActor::New()) {
+  vtkSmartPointer<vtkPolyDataMapper> mapper =
+      vtkSmartPointer<vtkPolyDataMapper>::New();
+  mapper->SetInputData(m_Poly);
+  m_Actor->SetMapper(mapper);
+  this->SetProp(m_Actor);
+}
+
+vtkQuadMesh::~vtkQuadMesh() {
+  m_Poly->Delete();
+  m_Actor->Delete();
+}
+
+void vtkQuadMesh::ReadFromFile(const char *filename) {
+  vtkSmartPointer<vtkPolyDataReader> reader =
+      vtkSmartPointer<vtkPolyDataReader>::New();
+  reader->SetFileName(filename);
+  reader->Update();
+  m_Poly->DeepCopy(reader->GetOutput());
+  vtk2uLib_update();
+}
+
+void vtkQuadMesh::ReadFromXMLFile(const char *filename) {
+  vtkSmartPointer<vtkXMLPolyDataReader> reader =
+      vtkSmartPointer<vtkXMLPolyDataReader>::New();
+  reader->SetFileName(filename);
+  reader->Update();
+  m_Poly->DeepCopy(reader->GetOutput());
+  vtk2uLib_update();
+}
+
+void vtkQuadMesh::ReadFromObjFile(const char *filename) {
+  vtkSmartPointer<vtkOBJReader> reader = vtkSmartPointer<vtkOBJReader>::New();
+  reader->SetFileName(filename);
+  reader->Update();
+  m_Poly->DeepCopy(reader->GetOutput());
+  vtk2uLib_update();
+}
+
+void vtkQuadMesh::ReadFromStlFile(const char *filename) {
+  vtkSmartPointer<vtkSTLReader> reader = vtkSmartPointer<vtkSTLReader>::New();
+  reader->SetFileName(filename);
+  reader->Update();
+  m_Poly->DeepCopy(reader->GetOutput());
+  vtk2uLib_update();
+}
+
+vtkPolyData *vtkQuadMesh::GetPolyData() const { return m_Poly; }
+
+void vtkQuadMesh::Update() { uLib2vtk_update(); }
+
+} // namespace Vtk
+} // namespace uLib

src/Vtk/Math/vtkQuadMesh.h

@@ -0,0 +1,69 @@
+/*//////////////////////////////////////////////////////////////////////////////
+// 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 VTKQUADMESH_H
+#define VTKQUADMESH_H
+
+#include "Math/QuadMesh.h"
+#include "Vtk/uLibVtkInterface.h"
+
+class vtkPolyData;
+class vtkActor;
+
+namespace uLib {
+namespace Vtk {
+
+class vtkQuadMesh : public Puppet, public Polydata {
+  typedef QuadMesh Content;
+
+public:
+  vtkQuadMesh(Content &content);
+  ~vtkQuadMesh();
+
+  void ReadFromFile(const char *filename);
+
+  void ReadFromXMLFile(const char *filename);
+
+  void ReadFromObjFile(const char *filename);
+
+  void ReadFromStlFile(const char *filename);
+
+  virtual class vtkPolyData *GetPolyData() const;
+
+  void Update();
+
+private:
+  void vtk2uLib_update();
+  void uLib2vtk_update();
+
+  QuadMesh &m_content;
+  vtkPolyData *m_Poly;
+  vtkActor *m_Actor;
+};
+
+} // namespace Vtk
+} // namespace uLib
+
+#endif // VTKQUADMESH_H

src/Vtk/Math/vtkVoxImage.cpp

@@ -50,33 +50,43 @@
 
 #include "Vtk/Math/vtkVoxImage.h"
 
+#include <vtkAutoInit.h>
+VTK_MODULE_INIT(vtkRenderingVolumeOpenGL2);
+VTK_MODULE_INIT(vtkRenderingOpenGL2);
+VTK_MODULE_INIT(vtkInteractionStyle);
+
 namespace uLib {
 namespace Vtk {
 
 void vtkVoxImage::GetContent() {
-  const int *dims = static_cast<const int *>(m_Content.GetDims().data());
-  m_Image->SetDimensions(dims);
+  Vector3i ev_dims = m_Content.GetDims();
+  m_Image->SetDimensions(ev_dims.data());
 
-  float *spacing = m_Content.GetSpacing().data();
-  m_Image->SetSpacing(spacing[0], spacing[1], spacing[2]);
+  Vector3f ev_spacing = m_Content.GetSpacing();
+  m_Image->SetSpacing(ev_spacing(0), ev_spacing(1), ev_spacing(2));
 
-  float *pos = m_Content.GetPosition().data();
-  m_Image->SetOrigin(pos[0], pos[1], pos[2]);
+  Vector3f ev_pos = m_Content.GetPosition();
+  m_Image->SetOrigin(ev_pos(0), ev_pos(1), ev_pos(2));
 
   vtkFloatArray *array =
       vtkFloatArray::SafeDownCast(m_Image->GetPointData()->GetScalars());
+  if (!array) {
+    array = vtkFloatArray::New();
+    m_Image->GetPointData()->SetScalars(array);
+    array->Delete();
+  }
+
   array->SetNumberOfTuples(m_Content.GetDims().prod());
   Vector3i index(0, 0, 0);
   int i = 0;
-  for (int zv = 0; zv < dims[2]; ++zv) {
-    for (int yv = 0; yv < dims[1]; ++yv) {
-      for (int xv = 0; xv < dims[0]; ++xv) {
+  for (int zv = 0; zv < ev_dims(2); ++zv) {
+    for (int yv = 0; yv < ev_dims(1); ++yv) {
+      for (int xv = 0; xv < ev_dims(0); ++xv) {
         index << xv, yv, zv;
         array->SetValue(i++, m_Content.GetValue(index));
       }
     }
   }
-  m_Image->GetPointData()->SetScalars(array);
 }
 
 void vtkVoxImage::SetContent() {