OpenCMT/uLib
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add quadmesh

Browse Source
This commit is contained in:
AndreaRigoni
2026-03-14 10:28:16 +00:00
parent 6bf9eaf309
commit 20d4967356
13 changed files with 600 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

31
CMakeLists.txt
View File

@@ -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_GEANT)
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)

4
app/CMakeLists.txt
View File

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

3
src/HEP/Geant/EmitterPrimary.hh
View File

@@ -20,4 +20,7 @@ class EmitterPrimary : public G4VUserPrimaryGeneratorAction
G4ParticleGun* fParticleGun; // Puntatore al cannone di particelle
};
#endif

368
src/Math/CLHEP_defs.hh Normal file
View File

@@ -0,0 +1,368 @@
#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

6
src/Math/CMakeLists.txt
View File

@@ -19,15 +19,19 @@ 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)

64
src/Math/QuadMesh.cpp Normal file
View File

@@ -0,0 +1,64 @@
/*//////////////////////////////////////////////////////////////////////////////
// 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);
}
}

57
src/Math/QuadMesh.h Normal file
View File

@@ -0,0 +1,57 @@
/*//////////////////////////////////////////////////////////////////////////////
// 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; }
private:
std::vector<Vector3f> m_Points;
std::vector<Vector4i> m_Quads;
};
}
#endif // QUADMESH_H

5
src/Math/Units.h
View File

@@ -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;

1
src/Math/testing/CMakeLists.txt
View File

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

51
src/Math/testing/QuadMeshTest.cpp Normal file
View File

@@ -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;
}

1
src/Math/testing/testing-prototype.h
View File

@@ -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++; }

11
src/Python/math_bindings.cpp
View File

@@ -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"
@@ -431,6 +432,16 @@ void init_math(py::module_ &m) {
.def("Triangles", &TriangleMesh::Triangles,
py::return_value_policy::reference_internal);
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);
py::class_<VoxRaytracer::RayData::Element>(m, "VoxRaytracerRayDataElement")
.def(py::init<>())
.def_readwrite("vox_id", &VoxRaytracer::RayData::Element::vox_id)

12
src/Python/testing/math_pybind_test.py
View File

@@ -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])