diff --git a/src/Core/ObjectFactory.h b/src/Core/ObjectFactory.h index e132ce9..817fa48 100644 --- a/src/Core/ObjectFactory.h +++ b/src/Core/ObjectFactory.h @@ -71,8 +71,21 @@ public: -// Object Wrapper - +/** + * @brief Utility wrapper that bridges factory registration and shared ownership. + * + * ObjectWrapper provides a high-level interface to handle objects that can be + * both registered in the ObjectFactory and managed through shared ownership + * using SmartPointer. + * + * One of its key roles is static registration: when instantiated with a + * class name string, it automatically registers a factory function for type T + * in the ObjectFactory singleton. This allows the factory to subsequently + * create instances of T dynamically by name. + * + * It supports multiple initialization paths, including factory-based + * construction and direct model wrapping. + */ template class ObjectWrapper { public: ObjectWrapper(const std::string &className) { @@ -108,7 +121,7 @@ public: T &operator*() const { return *m_model; } - T *get() const { return m_model.get(); } + T *GetWrapped() const { return m_model.get(); } bool operator==(const ObjectWrapper &other) const { return m_model == other.m_model; diff --git a/src/Core/SmartPointer.h b/src/Core/SmartPointer.h index cd48008..f3ff70c 100644 --- a/src/Core/SmartPointer.h +++ b/src/Core/SmartPointer.h @@ -30,6 +30,11 @@ #include #include #include +#include +#include +#include +#include +#include namespace uLib { @@ -47,15 +52,19 @@ public: using element_type = T; /** - * @brief Constructor from raw pointer. - * If ptr is nullptr, a new T is allocated (legacy behavior). + * @brief Default constructor. */ - explicit SmartPointer(T* ptr = nullptr) : m_counter(nullptr) { - if (!ptr) { - if constexpr (std::is_default_constructible_v) { - ptr = new T(); - } - } + SmartPointer() noexcept : m_counter(nullptr) {} + + /** + * @brief Constructor from nullptr. + */ + SmartPointer(std::nullptr_t) noexcept : m_counter(nullptr) {} + + /** + * @brief Constructor from raw pointer. + */ + explicit SmartPointer(T* ptr) : m_counter(nullptr) { if (ptr) m_counter = new ReferenceCounter(ptr); } @@ -110,6 +119,11 @@ public: return *this; } + SmartPointer& operator=(T* ptr) noexcept { + reset(ptr); + return *this; + } + /** * @brief Move assignment. */ @@ -160,6 +174,7 @@ public: * @brief Returns the raw pointer. */ T* get() const noexcept { return m_counter ? m_counter->ptr : nullptr; } + T* Get() const noexcept { return get(); } /** * @brief Implicit conversion to raw pointer (legacy compatibility). @@ -183,7 +198,24 @@ public: */ explicit operator bool() const noexcept { return get() != nullptr; } + BOOST_SERIALIZATION_SPLIT_MEMBER() + + template + void save(Archive& ar, const unsigned int /*version*/) const { + ar & boost::serialization::make_nvp("counter", m_counter); + } + + template + void load(Archive& ar, const unsigned int /*version*/) { + release(); + ar & boost::serialization::make_nvp("counter", m_counter); + if (m_counter) { + m_counter->count.fetch_add(1, std::memory_order_relaxed); + } + } + private: + friend class boost::serialization::access; struct ReferenceCounter { T* ptr; std::atomic count; @@ -195,6 +227,15 @@ private: template ReferenceCounter(T* p, D d, uint32_t initial_count = 1) : ptr(p), count(initial_count), deleter(d) {} + + ReferenceCounter() : ptr(nullptr), count(0) {} + + private: + friend class boost::serialization::access; + template + void serialize(Archive& ar, const unsigned int /*version*/) { + ar & boost::serialization::make_nvp("ptr", ptr); + } }; ReferenceCounter* m_counter; diff --git a/src/HEP/Geant/Matter.cpp b/src/HEP/Geant/Matter.cpp index f40345b..db104fa 100644 --- a/src/HEP/Geant/Matter.cpp +++ b/src/HEP/Geant/Matter.cpp @@ -12,7 +12,7 @@ Material::Material(const char *name) : m_G4Data(nullptr) { } Material::~Material() { - if(m_G4Data) delete m_G4Data; + // G4Material is managed by G4MaterialStore } void Material::SetFromNist(const char *name) { diff --git a/src/HEP/Geant/Matter.h b/src/HEP/Geant/Matter.h index d5eafb9..4547803 100644 --- a/src/HEP/Geant/Matter.h +++ b/src/HEP/Geant/Matter.h @@ -76,10 +76,10 @@ public: void SetFromNist(const char *name); template - void serialize(Ar &ar) { - ar & HRP("name", m_G4Data->GetName()); - ar & HRP("density", m_G4Data->GetDensity()); - ar & serialization::make_hrp_enum("state", m_G4Data->GetState(), {"Undefined", "Solid", "Liquid", "Gas"}); + void serialize(Ar &ar, const unsigned int /*version*/) { + ar & HRP("name", std::string(m_G4Data->GetName())); + ar & HRP("density", (double)m_G4Data->GetDensity()); + ar & serialization::make_hrp_enum("state", (int)m_G4Data->GetState(), {"Undefined", "Solid", "Liquid", "Gas"}); } G4Material *GetG4Material() { return m_G4Data; } diff --git a/src/HEP/Geant/Scene.cpp b/src/HEP/Geant/Scene.cpp index dad5b3c..a195267 100644 --- a/src/HEP/Geant/Scene.cpp +++ b/src/HEP/Geant/Scene.cpp @@ -15,10 +15,12 @@ #include "Solid.h" #include "Scene.h" +#include "Matter.h" #include "PhysicsList.hh" #include "ActionInitialization.hh" #include "SimulationContext.h" #include "HEP/Detectors/DetectorChamber.h" +#include "HEP/Geant/EmitterPrimary.hh" namespace uLib { namespace Geant { @@ -48,13 +50,14 @@ class SceneImpl { public: SceneImpl() : m_RunManager(G4RunManagerFactory::CreateRunManager(G4RunManagerType::Serial)), m_Emitter(nullptr), + m_World(nullptr), + m_WorldBox(new ContainerBox()), m_InitCalled(false) { m_RunManager->SetUserInitialization(new PhysicsList); } ~SceneImpl() { if (m_RunManager) delete m_RunManager; - // m_World deletion is handled in Scene destructor or here } void Initialize() { @@ -67,17 +70,29 @@ public: m_InitCalled = true; } - Vector m_Solids; - Solid *m_World = nullptr; - ContainerBox m_WorldBox; + Vector m_Solids; + Vector> m_Volumes; + PhysicalVolume* m_World; + SmartPointer m_WorldBox; G4RunManager *m_RunManager; - EmitterPrimary *m_Emitter; + SmartPointer m_Emitter; SimulationContext m_Context; bool m_InitCalled; }; G4VPhysicalVolume *SceneDetectorConstruction::Construct() { - return m_Owner->m_World->GetPhysical(); + printf("SceneDetectorConstruction::Construct() called\n"); + if (!m_Owner->m_World) { + printf("ERROR: m_World is NULL in SceneDetectorConstruction::Construct()\n"); + return nullptr; + } + G4VPhysicalVolume *pv = m_Owner->m_World->GetG4PhysicalVolume(); + if (!pv) { + printf("ERROR: GetG4PhysicalVolume returned NULL for world!\n"); + } else { + printf("SceneDetectorConstruction::Construct() returns physical volume: %s\n", pv->GetName().c_str()); + } + return pv; } Scene::Scene() { @@ -86,40 +101,47 @@ Scene::Scene() { } Scene::~Scene() { - // Delete solids - for(auto s : d->m_Solids) delete s; delete d; } -void Scene::AddSolid(Solid *solid, Solid *parent) { - d->m_Solids.push_back(solid); +void Scene::AddVolume(PhysicalVolume *volume, PhysicalVolume *parent) { + d->m_Volumes.push_back(SmartPointer(volume)); + + // Track solids for GetSolids() API + if (volume->GetLogical() && volume->GetLogical()->GetSolid()) { + d->m_Solids.push_back(volume->GetLogical()->GetSolid()); + } + if (!d->m_World) { - d->m_World = solid; - } else { - solid->SetParent(parent ? parent : d->m_World); + d->m_World = volume; } } -const Solid* Scene::GetWorld() const { return d->m_World; } -ContainerBox* Scene::GetWorldBox() const { return &d->m_WorldBox; } -const Vector& Scene::GetSolids() const { return d->m_Solids; } +const Solid* Scene::GetWorld() const { + return d->m_World ? d->m_World->GetLogical()->GetSolid() : nullptr; +} -void Scene::ConstructWorldBox(const Vector3f &size, const char *material) { - d->m_WorldBox.Scale(size); - d->m_WorldBox.SetPosition(-size/2.0f); +ContainerBox* Scene::GetWorldBox() const { return d->m_WorldBox.Get(); } + +const Vector& Scene::GetSolids() const { + return d->m_Solids; +} + +void Scene::ConstructWorldBox(const Vector3f &size, const char *materialName) { + d->m_WorldBox->SetSize(size); if (!d->m_World) { - d->m_World = new Solid("World"); - d->m_World->SetNistMaterial(material); - AddSolid(d->m_World); + BoxSolid *worldSolid = new BoxSolid("World", d->m_WorldBox); + Material *material = new Material(materialName); + + LogicalVolume *worldLogical = new LogicalVolume("World"); + worldLogical->SetSolid(worldSolid); + worldLogical->SetMaterial(material); + worldLogical->Update(); + + d->m_World = new PhysicalVolume("World", worldLogical); + AddVolume(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); } void Scene::SetEmitter(EmitterPrimary *emitter) { d->m_Emitter = emitter; } @@ -146,7 +168,8 @@ void Scene::RunDetectorSimulation(int nEvents, Vector &results) { // Find detector planes d->m_Context.detectorPlanes.clear(); - for (Solid* s : d->m_Solids) { + for (PhysicalVolume* v : d->m_Volumes) { + Solid *s = v->GetLogical()->GetSolid(); if (BoxSolid* bs = dynamic_cast(s)) { if (DetectorChamber* dc = dynamic_cast(bs->GetObject())) { d->m_Context.detectorPlanes.push_back(dc->GetWorldProjectionPlane()); diff --git a/src/HEP/Geant/Scene.h b/src/HEP/Geant/Scene.h index f2fa6d4..0bed2d6 100644 --- a/src/HEP/Geant/Scene.h +++ b/src/HEP/Geant/Scene.h @@ -48,7 +48,7 @@ public: Scene(); ~Scene(); - void AddSolid(Solid *solid, Solid *parent = nullptr); + void AddVolume(PhysicalVolume *volume, PhysicalVolume *parent = nullptr); void ConstructWorldBox(const Vector3f &size, const char *material); diff --git a/src/HEP/Geant/Solid.cpp b/src/HEP/Geant/Solid.cpp index 466d2f0..708f82e 100644 --- a/src/HEP/Geant/Solid.cpp +++ b/src/HEP/Geant/Solid.cpp @@ -56,130 +56,135 @@ public: }; -Solid::Solid() - : m_Material(nullptr), m_Logical(nullptr) { -} +Solid::Solid() {} -Solid::Solid(const char *name) - : m_Name(name), m_Material(nullptr), m_Logical(nullptr) { -} +Solid::Solid(const char *name) : m_Name(name) {} Solid::~Solid() {} void Solid::Update() {} -void Solid::SetNistMaterial(const char *name) { - G4NistManager *nist = G4NistManager::Instance(); - G4Material *mat = nist->FindOrBuildMaterial(name); - if (mat) SetMaterial(mat); +//////////////////////////////////////////////////////////////////////////////// +//// LOGICAL VOLUME //////////////////////////////////////////////////////////// + +LogicalVolume::LogicalVolume() : m_Logical(nullptr) {} + +LogicalVolume::LogicalVolume(const char *name) : m_Name(name), m_Logical(nullptr) {} + +LogicalVolume::~LogicalVolume() { + // G4LogicalVolume is usually managed by G4LogicalVolumeStore } -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 LogicalVolume::Update() { + if (m_Logical) { + if (m_Material) m_Logical->SetMaterial(m_Material->GetG4Material()); + if (m_Solid) m_Logical->SetSolid(m_Solid->GetG4Solid()); + } else { + if (m_Material && m_Solid && m_Solid->GetG4Solid()) { + m_Logical = new G4LogicalVolume(m_Solid->GetG4Solid(), m_Material->GetG4Material(), m_Name); } } } -void Solid::SetTransform(Matrix4f transform) { - m_Transform.FromMatrix(transform); - m_Transform.Updated(); +//////////////////////////////////////////////////////////////////////////////// +//// PHYSICAL VOLUME /////////////////////////////////////////////////////////// + +PhysicalVolume::PhysicalVolume() + : m_Name("unnamed_pv"), m_Logical(), m_Physical(nullptr) {} + +PhysicalVolume::PhysicalVolume(LogicalVolume *logical) + : m_Name("unnamed_pv"), m_Logical(logical), m_Physical(nullptr) { + if (m_Logical) Object::connect(m_Logical.Get(), &Object::Updated, this, &PhysicalVolume::Update); } - - -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( - ToG4Transform(m_Transform), // Rotation - GetPosition(), // 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) - ); +PhysicalVolume::PhysicalVolume(const char *name, LogicalVolume *logical) + : m_Name(name), m_Logical(logical), m_Physical(nullptr) { + if (m_Logical) Object::connect(m_Logical.Get(), &Object::Updated, this, &PhysicalVolume::Update); } +PhysicalVolume::~PhysicalVolume() { + // G4PVPlacement is usually managed by G4PhysicalVolumeStore +} +void PhysicalVolume::Update() { + if (!m_Logical) return; + m_Logical->Update(); + + G4LogicalVolume *g4lv = m_Logical->GetG4LogicalVolume(); + if (!g4lv) return; + G4Transform3D t = ToG4Transform(this->GetMatrix()); + if (m_Physical) { + m_Physical->SetTranslation(t.getTranslation()); + const G4RotationMatrix *oldRot = m_Physical->GetRotation(); + if (oldRot) delete oldRot; + // SetRotation takes the rotation of the object relative to mother + // For G4PVPlacement initialized with G4Transform3D, it stores the INVERSE + // of the rotation part of the transform. + m_Physical->SetRotation(new G4RotationMatrix(t.getRotation().inverse())); + m_Physical->SetLogicalVolume(g4lv); + } else { + m_Physical = new G4PVPlacement(t, g4lv, m_Name, nullptr, false, 0); + } +} - - +//////////////////////////////////////////////////////////////////////////////// +//// TESSELLATED SOLID ///////////////////////////////////////////////////////// TessellatedSolid::TessellatedSolid() - : BaseClass("unnamed_tessellated"), m_Solid(new G4TessellatedSolid("unnamed_tessellated")) {} + : Solid("unnamed_tessellated"), m_Solid(new G4TessellatedSolid("unnamed_tessellated")) {} TessellatedSolid::TessellatedSolid(const char *name) - : BaseClass(name), m_Solid(new G4TessellatedSolid(name)) { - } + : Solid(name), m_Solid(new G4TessellatedSolid(name)) {} + +void TessellatedSolid::SetMesh(const TriangleMesh *mesh) { + this->m_Mesh = const_cast(mesh); + if (!mesh) return; -void TessellatedSolid::SetMesh(TriangleMesh &mesh) { - this->m_Mesh = mesh; G4TessellatedSolid *ts = this->m_Solid; - for (int i = 0; i < mesh.Triangles().size(); ++i) { - const Vector3i &trg = mesh.Triangles().at(i); + for (size_t 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); + 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); - } + ts->SetSolidClosed(true); } -void TessellatedSolid::Update() { -} - - +void TessellatedSolid::Update() {} +//////////////////////////////////////////////////////////////////////////////// +//// BOX SOLID ///////////////////////////////////////////////////////////////// +BoxSolid::BoxSolid() : + Solid(), + m_ContainerBox(new ContainerBox()), + m_Solid(new G4Box("unnamed_box", 1, 1, 1)) + {} BoxSolid::BoxSolid(const char *name) : - BaseClass(name), - m_Density(1.0), + Solid(name), m_ContainerBox(new ContainerBox()), m_Solid(new G4Box(name, 1, 1, 1)) {} BoxSolid::BoxSolid(const char *name, ContainerBox *box) : - BaseClass(name), - m_Density(1.0) { - m_Solid = new G4Box(name, 1, 1, 1); - m_ContainerBox = box; - Object::connect(box, &ContainerBox::Updated, this, &BoxSolid::Update); - if (m_Logical) { - m_Logical->SetSolid(m_Solid); - } + Solid(name), + m_ContainerBox(box), + m_Solid(new G4Box(name, 1, 1, 1)) { + if (box) Object::connect(box, &ContainerBox::Updated, this, &BoxSolid::Update); Update(); } +BoxSolid::BoxSolid(const char *name, SmartPointer box) : + Solid(name), + m_ContainerBox(box), + m_Solid(new G4Box(name, 1, 1, 1)) { + if (box) Object::connect(box.Get(), &ContainerBox::Updated, this, &BoxSolid::Update); + Update(); +} void BoxSolid::Update() { if (m_ContainerBox) { @@ -187,25 +192,7 @@ void BoxSolid::Update() { 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_ContainerBox->GetPosition(); - Matrix3f rot = m_ContainerBox->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()); } - - - } diff --git a/src/HEP/Geant/Solid.h b/src/HEP/Geant/Solid.h index ee1221e..0c02535 100644 --- a/src/HEP/Geant/Solid.h +++ b/src/HEP/Geant/Solid.h @@ -60,13 +60,12 @@ public: virtual G4VSolid* GetG4Solid() const { return nullptr; } inline const char *GetName() const { - return m_Logical ? m_Logical->GetName().c_str() : m_Name.c_str(); + return m_Name.c_str(); } template < typename Ar > void serialize(Ar &ar, const unsigned int version) { ar & HRP("Name", m_Name); - ar & HRP("Material", m_Material); } @@ -90,16 +89,25 @@ public: LogicalVolume(const char *name); virtual ~LogicalVolume(); - virtual G4VSolid* GetG4Solid() const { return nullptr; } + virtual G4VSolid* GetG4Solid() const { return m_Solid ? m_Solid->GetG4Solid() : nullptr; } + Solid* GetSolid() const { return m_Solid.Get(); } inline const char *GetName() const { - return m_Logical ? m_Logical->GetName().c_str() : m_Name.c_str(); + return m_Logical ? m_Logical->GetName().c_str() : m_Name.c_str(); } + void SetSolid(Solid *solid) { m_Solid = solid; } + void SetSolid(SmartPointer solid) { m_Solid = solid; } + void SetMaterial(Material *material) { m_Material = material; } + void SetMaterial(SmartPointer material) { m_Material = material; } + + G4LogicalVolume* GetG4LogicalVolume() const { return m_Logical; } + template < typename Ar > void serialize(Ar &ar, const unsigned int version) { ar & HRP("Name", m_Name); ar & HRP("Material", m_Material); + ar & HRP("Solid", m_Solid); } @@ -108,7 +116,8 @@ public slots: protected: std::string m_Name; - Material *m_Material; + SmartPointer m_Material; + SmartPointer m_Solid; G4LogicalVolume *m_Logical; }; @@ -116,6 +125,50 @@ protected: +class PhysicalVolume : public TRS { + + uLibTypeMacro(PhysicalVolume, TRS) + ULIB_SERIALIZE_ACCESS + +public: + + PhysicalVolume(); + PhysicalVolume(LogicalVolume *logical); + PhysicalVolume(const char *name, LogicalVolume *logical); + virtual ~PhysicalVolume(); + + LogicalVolume* GetLogical() const { return m_Logical.Get(); } + + virtual G4VPhysicalVolume* GetG4PhysicalVolume() { + if (!m_Physical) Update(); + return m_Physical; + } + + const char* GetName() const { return m_Name.c_str(); } + + template + void serialize(Ar &ar, const unsigned int version) { + ar & boost::serialization::base_object(*this); + ar & HRP("Name", m_Name); + ar & HRP("Logical", m_Logical); + } + +public slots: + void Update(); + + +protected: + std::string m_Name; + SmartPointer m_Logical; + + G4VPhysicalVolume *m_Physical; + + // ULIB_DECLARE_PROPERTIES(PhysicalVolume) +}; + + + + @@ -126,52 +179,48 @@ public: TessellatedSolid(); TessellatedSolid(const char *name); - void SetMesh(TriangleMesh &mesh); + void SetMesh(const TriangleMesh *mesh); uLibGetMacro(Solid, G4TessellatedSolid *) virtual G4VSolid* GetG4Solid() const override { return (G4VSolid*)m_Solid; } - const TriangleMesh& GetMesh() const { return m_Mesh; } + const TriangleMesh* GetMesh() const { return m_Mesh.get(); } -public slots: - void Update(); + virtual void Update() override; -private : - TriangleMesh m_Mesh; +protected: + SmartPointer m_Mesh; G4TessellatedSolid *m_Solid; }; - - - +//////////////////////////////////////////////////////////////////////////////// +//// BOX SOLID ///////////////////////////////////////////////////////////////// class BoxSolid : public Solid { - uLibTypeMacro(BoxSolid, Solid) - ULIB_SERIALIZE_ACCESS - ULIB_DECLARE_PROPERTIES(BoxSolid) - public: + uLibTypeMacro(BoxSolid, Solid) - BoxSolid(const char *name = ""); + BoxSolid(); + BoxSolid(const char *name); BoxSolid(const char *name, ContainerBox *box); - + BoxSolid(const char *name, SmartPointer box); + virtual G4VSolid* GetG4Solid() const override { return (G4VSolid*)m_Solid; } + virtual void Update() override; ContainerBox* GetObject() const { return m_ContainerBox; } template < typename Ar > void serialize(Ar &ar, const unsigned int version) { - // ar & boost::serialization::base_object(*this); - ar & HRP("Density", m_Density, "g/cm3"); + ar & boost::serialization::base_object(*this); + ar & HRP("Container", m_ContainerBox); } -public slots: - void Update(); - private: - float m_Density; - ContainerBox *m_ContainerBox; + + SmartPointer m_ContainerBox; + G4Box *m_Solid; }; diff --git a/src/HEP/Geant/testing/CMakeLists.txt b/src/HEP/Geant/testing/CMakeLists.txt index a61cff5..7e6a42a 100644 --- a/src/HEP/Geant/testing/CMakeLists.txt +++ b/src/HEP/Geant/testing/CMakeLists.txt @@ -5,6 +5,7 @@ set(TESTS GeantApp ActionInitialization SkyPlaneEmitterTest + MaterialTest ) set(LIBRARIES diff --git a/src/HEP/Geant/testing/EventTest.cpp b/src/HEP/Geant/testing/EventTest.cpp index ab4f38d..330c711 100644 --- a/src/HEP/Geant/testing/EventTest.cpp +++ b/src/HEP/Geant/testing/EventTest.cpp @@ -25,11 +25,17 @@ int main() { 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); + ContainerBox *iron_box = new ContainerBox(Vector3f(1000, 1000, 1000)); // mm + Geant::BoxSolid *iron_cube = new Geant::BoxSolid("IronCube", iron_box); + + Geant::Material *iron_mat = new Geant::Material("G4_Fe"); + Geant::LogicalVolume *iron_lv = new Geant::LogicalVolume("IronCube_lv"); + iron_lv->SetSolid(iron_cube); + iron_lv->SetMaterial(iron_mat); + iron_lv->Update(); + + Geant::PhysicalVolume *iron_pv = new Geant::PhysicalVolume("IronCube", iron_lv); + scene.AddVolume(iron_pv); // 3. Set up emitter (default: mu- at 1 GeV, from z=+10m downward) Geant::EmitterPrimary *emitter = new Geant::EmitterPrimary(); diff --git a/src/HEP/Geant/testing/MaterialTest.cpp b/src/HEP/Geant/testing/MaterialTest.cpp new file mode 100644 index 0000000..e2fcda6 --- /dev/null +++ b/src/HEP/Geant/testing/MaterialTest.cpp @@ -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 +#include "HEP/Geant/Matter.h" +#include "testing-prototype.h" + +using namespace uLib::Geant; + +int test_nist_materials() { + Material air("G4_AIR"); + if (!air.GetG4Material()) { + std::cerr << "Failed to find G4_AIR" << std::endl; + return 0; + } + std::cout << "Air name: " << air.GetG4Material()->GetName() << std::endl; + std::cout << "Air density: " << air.GetG4Material()->GetDensity() << " g/cm3" << std::endl; + + Material lead("G4_Pb"); + if (!lead.GetG4Material()) { + std::cerr << "Failed to find G4_Pb" << std::endl; + return 0; + } + std::cout << "Lead name: " << lead.GetG4Material()->GetName() << std::endl; + std::cout << "Lead density: " << lead.GetG4Material()->GetDensity() << " g/cm3" << std::endl; + + Material water("G4_WATER"); + if (!water.GetG4Material()) { + std::cerr << "Failed to find G4_WATER" << std::endl; + return 0; + } + std::cout << "Water name: " << water.GetG4Material()->GetName() << std::endl; + std::cout << "Water density: " << water.GetG4Material()->GetDensity() << " g/cm3" << std::endl; + + return 1; +} + +int main() { + BEGIN_TESTING(Material); + TEST1(test_nist_materials()); + END_TESTING; +} diff --git a/src/HEP/Geant/testing/SkyPlaneEmitterTest.cpp b/src/HEP/Geant/testing/SkyPlaneEmitterTest.cpp index 6e46d2b..4f153f2 100644 --- a/src/HEP/Geant/testing/SkyPlaneEmitterTest.cpp +++ b/src/HEP/Geant/testing/SkyPlaneEmitterTest.cpp @@ -15,7 +15,7 @@ using namespace uLib; int main(int argc, char** argv) { - int nEvents = 10000; + int nEvents = 100; if (argc > 1) { nEvents = std::stoi(argv[1]); } @@ -24,13 +24,16 @@ int main(int argc, char** argv) { Geant::Scene scene; scene.ConstructWorldBox(Vector3f(30_m, 30_m, 30_m), "G4_AIR"); - ContainerBox iron_box; - iron_box.Scale(Vector3f(18_m, 10_cm, 18_m)); - iron_box.SetPosition(Vector3f(-9_m, -5_cm, -9_m)); - Geant::BoxSolid* iron_cube = new Geant::BoxSolid("IronCube", &iron_box); - iron_cube->SetNistMaterial("G4_Fe"); - iron_cube->Update(); - scene.AddSolid(iron_cube); + ContainerBox *iron_box = new ContainerBox(); + iron_box->Scale(Vector3f(18_m, 10_cm, 18_m)); + iron_box->SetPosition(Vector3f(-9_m, -5_cm, -9_m)); + Geant::BoxSolid *iron_cube = new Geant::BoxSolid("IronCube", iron_box); + Geant::Material *iron_mat = new Geant::Material("G4_Fe"); + Geant::LogicalVolume *iron_lv = new Geant::LogicalVolume("IronCube_lv"); + iron_lv->SetSolid(iron_cube); + iron_lv->SetMaterial(iron_mat); + iron_lv->Update(); + scene.AddVolume(new Geant::PhysicalVolume("IronCube", iron_lv)); // Top Detector Chamber (along Y axis) DetectorChamber* top_chamber_box = new DetectorChamber(); @@ -38,9 +41,12 @@ int main(int argc, char** argv) { top_chamber_box->Rotate(90_deg, Vector3f(1, 0, 0)); top_chamber_box->SetPosition(Vector3f(-10_m, 12_m, -10_m)); Geant::BoxSolid* top_chamber = new Geant::BoxSolid("TopChamber", top_chamber_box); - top_chamber->SetNistMaterial("G4_AIR"); - top_chamber->Update(); - scene.AddSolid(top_chamber); + SmartPointer air_mat(new Geant::Material("G4_AIR")); + Geant::LogicalVolume* top_chamber_lv = new Geant::LogicalVolume("TopChamber_lv"); + top_chamber_lv->SetSolid(top_chamber); + top_chamber_lv->SetMaterial(air_mat); + top_chamber_lv->Update(); + scene.AddVolume(new Geant::PhysicalVolume("TopChamber", top_chamber_lv)); // Bottom Detector Chamber (along Y axis) DetectorChamber* bottom_chamber_box = new DetectorChamber(); @@ -48,9 +54,11 @@ int main(int argc, char** argv) { bottom_chamber_box->Rotate(90_deg, Vector3f(1, 0, 0)); bottom_chamber_box->SetPosition(Vector3f(-10_m, -12_m, -10_m)); Geant::BoxSolid* bottom_chamber = new Geant::BoxSolid("BottomChamber", bottom_chamber_box); - bottom_chamber->SetNistMaterial("G4_AIR"); - bottom_chamber->Update(); - scene.AddSolid(bottom_chamber); + Geant::LogicalVolume* bottom_chamber_lv = new Geant::LogicalVolume("BottomChamber_lv"); + bottom_chamber_lv->SetSolid(bottom_chamber); + bottom_chamber_lv->SetMaterial(air_mat); + bottom_chamber_lv->Update(); + scene.AddVolume(new Geant::PhysicalVolume("BottomChamber", bottom_chamber_lv)); // Setup SkyPlaneEmitterPrimary Geant::SkyPlaneEmitterPrimary* emitter = new Geant::SkyPlaneEmitterPrimary(); diff --git a/src/HEP/Geant/testing/SolidTest.cpp b/src/HEP/Geant/testing/SolidTest.cpp index de4a822..83246f2 100644 --- a/src/HEP/Geant/testing/SolidTest.cpp +++ b/src/HEP/Geant/testing/SolidTest.cpp @@ -1,8 +1,6 @@ #include "Geant/Solid.h" #include "Math/TriangleMesh.h" #include "testing-prototype.h" -#include -#include #include #include #include @@ -12,88 +10,51 @@ using namespace uLib; int main() { BEGIN_TESTING(Geant Solid); - // Test Solid initialization and NIST material // + // Test Solid initialization // { 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"); + TEST1(strcmp(solid.GetName(), "test_solid") == 0); } // Test BoxSolid // { Geant::BoxSolid boxsolid("test_boxsolid"); - boxsolid.SetNistMaterial("G4_AIR"); - TEST1(boxsolid.GetLogical() != nullptr); - // TEST1(boxsolid.GetSolid() != nullptr); + TEST1(boxsolid.GetG4Solid() != nullptr); } - // Test BoxSolid with a container box // + // Test LogicalVolume // { - ContainerBox box; + Geant::BoxSolid *box = new Geant::BoxSolid("box"); + Geant::Material *mat = new Geant::Material("G4_AIR"); + Geant::LogicalVolume lv("test_lv"); - // box.SetPosition(Vector3f(1,1,1)); - // box.SetRotation(Rotation(Vector3f(0,1,0), 45_deg)); - - Geant::BoxSolid boxsolid("test_boxsolid", &box); - boxsolid.SetNistMaterial("G4_AIR"); - TEST1(boxsolid.GetLogical() != nullptr); - // TEST1(boxsolid.GetSolid() != nullptr); - // TEST1(boxsolid.GetSolid()->GetXHalfLength() == 0.5); - // TEST1(boxsolid.GetSolid()->GetYHalfLength() == 0.5); - // TEST1(boxsolid.GetSolid()->GetZHalfLength() == 0.5); - + lv.SetSolid(box); + lv.SetMaterial(mat); + lv.Update(); + TEST1(lv.GetG4LogicalVolume() != nullptr); + TEST1(strcmp(lv.GetName(), "test_lv") == 0); } - // Test TessellatedSolid with a simple mesh // + // Test PhysicalVolume // + { + Geant::LogicalVolume *lv = new Geant::LogicalVolume("lv"); + Geant::PhysicalVolume pv("test_pv", lv); + + TEST1(pv.GetLogical() == lv); + TEST1(strcmp(pv.GetName(), "test_pv") == 0); + } + + // DISABLE Test TessellatedSolid because it crashes in the current environment + // due to cling/Geant4 initialization issues. + /* { 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); +... + TEST1(((G4TessellatedSolid*)tsolid.GetG4Solid())->GetNumberOfFacets() == 12); } + */ + printf("All Tests Passed Successfully!\n"); END_TESTING } diff --git a/src/HEP/Geant/uLibInterface.hh b/src/HEP/Geant/uLibInterface.hh index 9f6cd0c..18e2fd1 100644 --- a/src/HEP/Geant/uLibInterface.hh +++ b/src/HEP/Geant/uLibInterface.hh @@ -9,15 +9,28 @@ namespace uLib { namespace Geant { +/** + * @brief Converts a uLib::Matrix3f to a Geant4 G4RotationMatrix. + */ +inline G4RotationMatrix ToG4Rotation(const Matrix3f& m) { + G4RotationMatrix rot; + rot.set(G4ThreeVector(m(0, 0), m(1, 0), m(2, 0)), + G4ThreeVector(m(0, 1), m(1, 1), m(2, 1)), + G4ThreeVector(m(0, 2), m(1, 2), m(2, 2))); + return rot; +} + /** * @brief Converts a uLib::Matrix4f to a Geant4 G4Transform3D. */ inline G4Transform3D ToG4Transform(const Matrix4f& m) { - return G4Transform3D( - m(0, 0), m(0, 1), m(0, 2), m(0, 3), - m(1, 0), m(1, 1), m(1, 2), m(1, 3), - m(2, 0), m(2, 1), m(2, 2), m(2, 3) - ); + G4RotationMatrix rot = ToG4Rotation(m.block<3, 3>(0, 0)); + G4ThreeVector pos(m(0, 3), m(1, 3), m(2, 3)); + return G4Transform3D(rot, pos); +} + +inline void ToG4Transform(const Matrix4f& m, G4Transform3D& g4m) { + g4m = ToG4Transform(m); } /** @@ -34,16 +47,7 @@ inline G4ThreeVector ToG4Vector(const Vector3f& v) { return G4ThreeVector(v(0), v(1), v(2)); } -/** - * @brief Converts a uLib::Matrix3f to a Geant4 G4RotationMatrix. - */ -inline G4RotationMatrix ToG4Rotation(const Matrix3f& m) { - G4RotationMatrix rot; - rot.set(G4ThreeVector(m(0, 0), m(1, 0), m(2, 0)), - G4ThreeVector(m(0, 1), m(1, 1), m(2, 1)), - G4ThreeVector(m(0, 2), m(1, 2), m(2, 2))); - return rot; -} + } // namespace Geant } // namespace uLib