From d8ef413216b13c2a1b8ff36ddcc8710c106fa033 Mon Sep 17 00:00:00 2001 From: AndreaRigoni Date: Mon, 16 Mar 2026 17:51:53 +0000 Subject: [PATCH] vtkGeantEvent --- app/gcompose/src/main.cpp | 2 +- src/HEP/Geant/ActionInitialization.cpp | 6 +- src/HEP/Geant/GeantEvent.cpp | 25 ++ src/HEP/Geant/GeantEvent.h | 2 + src/HEP/Geant/Scene.cpp | 20 +- src/HEP/Geant/Scene.h | 7 +- src/HEP/Geant/Solid.cpp | 56 +++-- src/HEP/Geant/Solid.h | 1 + src/HEP/Geant/SteppingAction.cpp | 75 +++--- src/HEP/Geant/SteppingAction.hh | 5 +- src/HEP/Geant/testing/EventTest.cpp | 4 +- src/HEP/Geant/testing/GeantApp.cpp | 4 +- src/HEP/Geant/testing/SolidTest.cpp | 6 +- src/Math/ContainerBox.h | 6 +- src/Math/Transform.h | 6 +- src/Math/testing/ContainerBoxTest.cpp | 35 ++- src/Math/testing/testing-prototype.h | 4 +- src/Vtk/CMakeLists.txt | 2 + src/Vtk/HEP/Detectors/vtkMuonEvent.h | 1 + src/Vtk/HEP/Detectors/vtkMuonScatter.h | 1 + .../HEP/Geant/testing/vtkGeantEventTest.cpp | 42 ++-- src/Vtk/HEP/Geant/vtkGeantEvent.cpp | 4 +- src/Vtk/HEP/Geant/vtkGeantEvent.h | 1 + src/Vtk/Math/vtkQuadMesh.h | 1 + src/Vtk/Math/vtkTriangleMesh.h | 1 + src/Vtk/uLibVtkInterface.cxx | 224 ++++++++++++++---- src/Vtk/uLibVtkInterface.h | 33 ++- src/Vtk/uLibVtkViewer.cpp | 2 + src/Vtk/vtkContainerBox.h | 1 + src/Vtk/vtkPolydata.cpp | 70 ++++++ src/Vtk/vtkPolydata.h | 17 +- 31 files changed, 509 insertions(+), 155 deletions(-) create mode 100644 src/HEP/Geant/GeantEvent.cpp create mode 100644 src/Vtk/vtkPolydata.cpp diff --git a/app/gcompose/src/main.cpp b/app/gcompose/src/main.cpp index 3ca8645..afa7f68 100644 --- a/app/gcompose/src/main.cpp +++ b/app/gcompose/src/main.cpp @@ -42,7 +42,7 @@ int main(int argc, char** argv) { Geant::Scene scene; - scene.ConstructWorldBox(&world_box, "G4_AIR"); + scene.ConstructWorldBox(world_box.GetSize(), "G4_AIR"); scene.Initialize(); // 2. Initialize MainWindow (contains embedded VTK QViewport) diff --git a/src/HEP/Geant/ActionInitialization.cpp b/src/HEP/Geant/ActionInitialization.cpp index 2447b39..10901db 100644 --- a/src/HEP/Geant/ActionInitialization.cpp +++ b/src/HEP/Geant/ActionInitialization.cpp @@ -27,9 +27,11 @@ void ActionInitialization::Build() const { SetUserAction(new EmitterPrimary()); } - // Register stepping action to collect scattering data + // Register actions if (m_Output) { - SetUserAction(new SteppingAction(m_Output)); + SteppingAction *sa = new SteppingAction(m_Output); + SetUserAction(static_cast(sa)); + SetUserAction(static_cast(sa)); } } diff --git a/src/HEP/Geant/GeantEvent.cpp b/src/HEP/Geant/GeantEvent.cpp new file mode 100644 index 0000000..571423a --- /dev/null +++ b/src/HEP/Geant/GeantEvent.cpp @@ -0,0 +1,25 @@ + +#include "HEP/Geant/GeantEvent.h" + +#include +#include + +using namespace uLib; + +void GeantEvent::Print(const size_t size) const { + std::cout << "Event " << m_EventID << ":" << std::endl; + std::cout << " Momentum: " << m_Momentum << std::endl; + std::cout << " GenVector: " << m_GenVector << std::endl; + std::cout << " Path: " << std::endl; + + size_t limit = m_Path.size(); + if(size > 0 && size < m_Path.size()) { + limit = size; + } + for (size_t i = 0; i < limit; ++i) { + std::cout << " " << i << ": " << m_Path[i].m_Length << " " << m_Path[i].m_Momentum << " " << m_Path[i].m_Direction << " " << m_Path[i].m_SolidName << std::endl; + } + if (limit < m_Path.size()) { + std::cout << " ... (" << m_Path.size() - limit << " more deltas)" << std::endl; + } +} diff --git a/src/HEP/Geant/GeantEvent.h b/src/HEP/Geant/GeantEvent.h index 9042e7a..35ce0cb 100644 --- a/src/HEP/Geant/GeantEvent.h +++ b/src/HEP/Geant/GeantEvent.h @@ -71,6 +71,8 @@ public: uLibConstRefMacro(GenVector, HLine3f) uLibConstRefMacro(Path, Vector) + void Print(const size_t size = 10) const; + private: Id_t m_EventID; Scalarf m_Momentum; diff --git a/src/HEP/Geant/Scene.cpp b/src/HEP/Geant/Scene.cpp index db5bf2b..ac2900e 100644 --- a/src/HEP/Geant/Scene.cpp +++ b/src/HEP/Geant/Scene.cpp @@ -95,6 +95,7 @@ public: // members // Vector m_Solids; Solid *m_World = nullptr; + ContainerBox m_WorldBox; G4RunManager *m_RunManager; EmitterPrimary *m_Emitter; Vector *m_Output; @@ -125,16 +126,22 @@ void Scene::AddSolid(Solid *solid, Solid *parent) { } } -void Scene::ConstructWorldBox(const ContainerBox *box, const char *material) { +const Solid* Scene::GetWorld() const { return d->m_World; } + +ContainerBox* Scene::GetWorldBox() const { return &d->m_WorldBox; } + +void Scene::ConstructWorldBox(const Vector3f &size, const char *material) { // Get nist material manager + + d->m_WorldBox.Scale(size); + d->m_WorldBox.SetPosition(-size/2.0f); + if (!d->m_World) { d->m_World = new Solid("World"); d->m_World->SetNistMaterial(material); AddSolid(d->m_World); } - Vector3f size = box->GetSize(); - G4Box *solidWorld = new G4Box("World", 0.5 * size(0), 0.5 * size(1), @@ -157,9 +164,10 @@ void Scene::ConstructWorldBox(const ContainerBox *box, const char *material) { true); d->m_World->SetPhysical(physWorld); - - Matrix4f transform = box->GetMatrix(); - d->m_World->SetTransform(transform); + + // no transforms are allowed for the world box + // Matrix4f transform = box->GetMatrix(); + // d->m_World->SetTransform(transform); } void Scene::SetEmitter(EmitterPrimary *emitter) { diff --git a/src/HEP/Geant/Scene.h b/src/HEP/Geant/Scene.h index 3467dee..7f4f689 100644 --- a/src/HEP/Geant/Scene.h +++ b/src/HEP/Geant/Scene.h @@ -47,7 +47,12 @@ public: void AddSolid(Solid *solid, Solid *parent = nullptr); - void ConstructWorldBox(const ContainerBox *box, const char *material); + void ConstructWorldBox(const Vector3f &size, const char *material); + + /// Get the world box + const Solid* GetWorld() const; + + ContainerBox* GetWorldBox() const; /// Set the primary generator (emitter) for the simulation. /// The Scene does NOT take ownership of the emitter. diff --git a/src/HEP/Geant/Solid.cpp b/src/HEP/Geant/Solid.cpp index 1e66c90..cf78862 100644 --- a/src/HEP/Geant/Solid.cpp +++ b/src/HEP/Geant/Solid.cpp @@ -55,10 +55,17 @@ public: }; Solid::Solid() - : m_Name("unnamed_solid"), m_Material(NULL), m_Logical(NULL), m_Physical(NULL) {} + : m_Name("unnamed_solid"), m_Material(NULL), m_Logical(NULL), m_Physical(NULL), + m_Position(new G4ThreeVector(0,0,0)), m_Rotation(NULL) {} Solid::Solid(const char *name) - : m_Name(name), m_Material(NULL), m_Logical(NULL), m_Physical(NULL) {} + : m_Name(name), m_Material(NULL), m_Logical(NULL), m_Physical(NULL), + m_Position(new G4ThreeVector(0,0,0)), m_Rotation(NULL) {} + +Solid::~Solid() { + if (m_Position) delete m_Position; + if (m_Rotation) delete m_Rotation; +} void Solid::SetNistMaterial(const char *name) { G4NistManager *nist = G4NistManager::Instance(); @@ -81,21 +88,22 @@ void Solid::SetTransform(Matrix4f transform) { uLib::AffineTransform t; t.SetMatrix(transform); - // 2. Extracto position and rotation for Geant4 + // 2. Extract position and rotation for Geant4 Vector3f pos = t.GetPosition(); - G4ThreeVector g4pos(pos(0), pos(1), pos(2)); + if (!m_Position) m_Position = new G4ThreeVector(); + *m_Position = G4ThreeVector(pos(0), pos(1), pos(2)); // Create a G4 rotation matrix from the 4x4 matrix Matrix3f m = t.GetRotation(); - G4RotationMatrix* rot = new G4RotationMatrix(); - rot->set(G4ThreeVector(m(0,0), m(1,0), m(2,0)), - G4ThreeVector(m(0,1), m(1,1), m(2,1)), - G4ThreeVector(m(0,2), m(1,2), m(2,2))); + if (!m_Rotation) m_Rotation = new G4RotationMatrix(); + m_Rotation->set(G4ThreeVector(m(0,0), m(1,0), m(2,0)), + G4ThreeVector(m(0,1), m(1,1), m(2,1)), + G4ThreeVector(m(0,2), m(1,2), m(2,2))); - // 3. Se l'oggetto è già stato piazzato, aggiorniamo la sua trasformazione + // 3. If object is already placed, update its transformation if (m_Physical) { - m_Physical->SetTranslation(g4pos); - m_Physical->SetRotation(rot); + m_Physical->SetTranslation(*m_Position); + m_Physical->SetRotation(m_Rotation); } } @@ -121,8 +129,8 @@ void Solid::SetParent(Solid *parent) { // G4PVPlacement m_Physical = new G4PVPlacement( - nullptr, // Rotation - G4ThreeVector(0,0,0), // Position (translation) inside the parent + m_Rotation, // Rotation + *m_Position, // Position (translation) inside the parent m_Logical, // The logical volume of this solid (the child) GetName(), // Name of the physical volume parentLogical, // The logical volume of the parent (nullptr if it's the World volume) @@ -139,7 +147,8 @@ void Solid::SetParent(Solid *parent) { TessellatedSolid::TessellatedSolid(const char *name) - : BaseClass(name), m_Solid(new G4TessellatedSolid(name)) {} + : BaseClass(name), m_Solid(new G4TessellatedSolid(name)) { + } void TessellatedSolid::SetMesh(TriangleMesh &mesh) { G4TessellatedSolid *ts = this->m_Solid; @@ -161,12 +170,13 @@ void TessellatedSolid::SetMesh(TriangleMesh &mesh) { BoxSolid::BoxSolid(const char *name, ContainerBox *box) : BaseClass(name) { - m_Solid = new G4Box(name, 0.5, 0.5, 0.5); + m_Solid = new G4Box(name, 1,1,1); m_Object = box; Object::connect(box, &ContainerBox::Updated, this, &BoxSolid::Update); if (m_Logical) { m_Logical->SetSolid(m_Solid); } + Update(); } void BoxSolid::Update() { @@ -176,8 +186,20 @@ void BoxSolid::Update() { m_Solid->SetYHalfLength(size(1) * 0.5); m_Solid->SetZHalfLength(size(2) * 0.5); - this->SetTransform(m_Object->GetMatrix()); - // this->m_Logical->SetSolid(m_Solid); + // Geant4 placement is relative to center. uLib Box is anchored at corner. + // 1. Get position and rotation (clean, without scale) + Vector3f pos = m_Object->GetPosition(); + Matrix3f rot = m_Object->GetRotation(); + + // 2. Center = Corner + Rotation * (Half-Size) + // We must rotate the offset vector because uLib box can be rotated. + Vector3f center = pos + rot * (size * 0.5); + + uLib::AffineTransform t; + t.SetPosition(center); + t.SetRotation(rot); + + this->SetTransform(t.GetMatrix()); } } diff --git a/src/HEP/Geant/Solid.h b/src/HEP/Geant/Solid.h index 3f7a966..4f76adf 100644 --- a/src/HEP/Geant/Solid.h +++ b/src/HEP/Geant/Solid.h @@ -45,6 +45,7 @@ class Solid : public Object { public: Solid(); Solid(const char *name); + virtual ~Solid(); void SetNistMaterial(const char *name); void SetMaterial(G4Material *material); diff --git a/src/HEP/Geant/SteppingAction.cpp b/src/HEP/Geant/SteppingAction.cpp index a12f125..3d1157c 100644 --- a/src/HEP/Geant/SteppingAction.cpp +++ b/src/HEP/Geant/SteppingAction.cpp @@ -7,6 +7,8 @@ #include "G4LogicalVolume.hh" #include "G4SystemOfUnits.hh" #include "G4ParticleDefinition.hh" +#include +#include namespace uLib { namespace Geant { @@ -18,9 +20,46 @@ SteppingAction::SteppingAction(Vector *output) m_LastEventID(-1) {} -SteppingAction::~SteppingAction() { - // Push the last event if any steps were collected - if (m_LastEventID >= 0 && !m_Current.m_Path.empty()) { +SteppingAction::~SteppingAction() {} + +void SteppingAction::BeginOfEventAction(const G4Event *event) { + if (!event || !m_Output) return; + + // Start a new GeantEvent + m_Current = GeantEvent(); + m_Current.m_EventID = static_cast(event->GetEventID()); + + // Set initial momentum and generation vector from primary vertex + if (event->GetNumberOfPrimaryVertex() > 0) { + G4PrimaryVertex *vtx = event->GetPrimaryVertex(0); + G4ThreeVector pos = vtx->GetPosition(); + m_Current.m_GenVector.origin = HPoint3f(pos.x(), pos.y(), pos.z()); + + if (vtx->GetNumberOfParticle() > 0) { + G4PrimaryParticle *prim = vtx->GetPrimary(0); + G4ThreeVector mom = prim->GetMomentumDirection(); + m_Current.m_GenVector.direction = HVector3f(mom.x(), mom.y(), mom.z()); + m_Current.m_Momentum = static_cast(prim->GetTotalMomentum() / MeV); + } + } +} + +void SteppingAction::EndOfEventAction(const G4Event *event) { + if (m_Output && !m_Current.m_Path.empty()) { + std::cout << "[Geant] Finished Event " << m_Current.m_EventID + << " with " << m_Current.m_Path.size() << " steps." << std::endl; + + // Check if we hit anything other than World + std::set volumes; + for (const auto& delta : m_Current.m_Path) { + if (!delta.m_SolidName.empty()) volumes.insert(delta.m_SolidName); + } + if (volumes.size() > 1) { + std::cout << " - Hit volumes: "; + for (const auto& v : volumes) std::cout << v << " "; + std::cout << std::endl; + } + m_Output->push_back(m_Current); } } @@ -34,36 +73,6 @@ void SteppingAction::UserSteppingAction(const G4Step *step) { // Only record primary particle (muon) if (track->GetParentID() != 0) return; - // Get current event ID - const G4Event *event = G4RunManager::GetRunManager()->GetCurrentEvent(); - if (!event) return; - int eventID = event->GetEventID(); - - // Detect new event — push completed event and start fresh - if (eventID != m_LastEventID) { - if (m_LastEventID >= 0 && !m_Current.m_Path.empty()) { - m_Output->push_back(m_Current); - } - // Start a new GeantEvent - m_Current = GeantEvent(); - m_Current.m_EventID = static_cast(eventID); - m_LastEventID = eventID; - - // Set initial momentum and generation vector from primary vertex - if (event->GetNumberOfPrimaryVertex() > 0) { - G4PrimaryVertex *vtx = event->GetPrimaryVertex(0); - G4ThreeVector pos = vtx->GetPosition(); - m_Current.m_GenVector.origin = HPoint3f(pos.x(), pos.y(), pos.z()); - - if (vtx->GetNumberOfParticle() > 0) { - G4PrimaryParticle *prim = vtx->GetPrimary(0); - G4ThreeVector mom = prim->GetMomentumDirection(); - m_Current.m_GenVector.direction = HVector3f(mom.x(), mom.y(), mom.z()); - m_Current.m_Momentum = static_cast(prim->GetTotalMomentum() / MeV); - } - } - } - // Record a Delta for this step GeantEvent::Delta delta; diff --git a/src/HEP/Geant/SteppingAction.hh b/src/HEP/Geant/SteppingAction.hh index 920113c..84fbdc4 100644 --- a/src/HEP/Geant/SteppingAction.hh +++ b/src/HEP/Geant/SteppingAction.hh @@ -2,6 +2,7 @@ #define U_GEANT_STEPPINGACTION_HH #include "G4UserSteppingAction.hh" +#include "G4UserEventAction.hh" #include "Core/Vector.h" #include "GeantEvent.h" @@ -10,13 +11,15 @@ namespace Geant { /// SteppingAction collects scattering data at each Geant4 step and /// builds GeantEvent objects in the output buffer. -class SteppingAction : public G4UserSteppingAction { +class SteppingAction : public G4UserSteppingAction, public G4UserEventAction { public: /// @param output pointer to the results vector owned by the Scene SteppingAction(Vector *output); virtual ~SteppingAction(); virtual void UserSteppingAction(const G4Step *step) override; + virtual void BeginOfEventAction(const G4Event *event) override; + virtual void EndOfEventAction(const G4Event *event) override; private: Vector *m_Output; ///< destination for finished events diff --git a/src/HEP/Geant/testing/EventTest.cpp b/src/HEP/Geant/testing/EventTest.cpp index e653228..ab4f38d 100644 --- a/src/HEP/Geant/testing/EventTest.cpp +++ b/src/HEP/Geant/testing/EventTest.cpp @@ -5,6 +5,7 @@ #include "Math/ContainerBox.h" #include "Math/TriangleMesh.h" #include "Math/Dense.h" +#include "Math/Units.h" #include "testing-prototype.h" #include #include @@ -20,9 +21,8 @@ int main() { // Test: Scene with iron cube in air, launch muons, collect events // { // 1. Create world box (air, 30m x 30m x 30m) - ContainerBox world_box(Vector3f(30000, 30000, 30000)); // mm Geant::Scene scene; - scene.ConstructWorldBox(&world_box, "G4_AIR"); + 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 diff --git a/src/HEP/Geant/testing/GeantApp.cpp b/src/HEP/Geant/testing/GeantApp.cpp index 2ef35e2..39980fa 100644 --- a/src/HEP/Geant/testing/GeantApp.cpp +++ b/src/HEP/Geant/testing/GeantApp.cpp @@ -8,10 +8,8 @@ using namespace uLib; int main() { - uLib::ContainerBox world_box(Vector3f(100, 100, 100)); uLib::Geant::Scene scene; - - scene.ConstructWorldBox(&world_box, "G4_AIR"); + scene.ConstructWorldBox(Vector3f(100, 100, 100), "G4_AIR"); scene.Initialize(); return 0; diff --git a/src/HEP/Geant/testing/SolidTest.cpp b/src/HEP/Geant/testing/SolidTest.cpp index 181edc7..bc7859c 100644 --- a/src/HEP/Geant/testing/SolidTest.cpp +++ b/src/HEP/Geant/testing/SolidTest.cpp @@ -15,9 +15,12 @@ int main() { // Test Solid initialization and NIST material // { Geant::Solid solid("test_solid"); - TEST1(solid.GetLogical() != nullptr); + // Logical volume is not created until material and solid are set + TEST1(solid.GetLogical() == nullptr); solid.SetNistMaterial("G4_AIR"); + // Still null because base Solid has no GetG4Solid() + TEST1(solid.GetLogical() == nullptr); TEST1(solid.GetMaterial() != nullptr); TEST1(solid.GetMaterial()->GetName() == "G4_AIR"); } @@ -25,6 +28,7 @@ int main() { // Test TessellatedSolid with a simple mesh // { Geant::TessellatedSolid tsolid("test_tessellated"); + tsolid.SetNistMaterial("G4_AIR"); TEST1(tsolid.GetLogical() != nullptr); TEST1(tsolid.GetSolid() != nullptr); diff --git a/src/Math/ContainerBox.h b/src/Math/ContainerBox.h index fc20538..af3a8f3 100644 --- a/src/Math/ContainerBox.h +++ b/src/Math/ContainerBox.h @@ -101,7 +101,11 @@ public: * @brief Gets the current size (scale) of the box. * @return The size vector. */ - inline Vector3f GetSize() const { return m_LocalT.GetScale(); } + inline Vector3f GetSize() const { + Vector3f s = this->GetScale(); + Vector3f ls = m_LocalT.GetScale(); + return Vector3f(s(0) * ls(0), s(1) * ls(1), s(2) * ls(2)); + } /** * @brief Swaps two local axes of the box. diff --git a/src/Math/Transform.h b/src/Math/Transform.h index 9cca0df..3746a23 100644 --- a/src/Math/Transform.h +++ b/src/Math/Transform.h @@ -105,7 +105,11 @@ public: inline void Scale(const Vector3f v) { this->m_T.scale(v); } - inline Vector3f GetScale() const { return this->m_T.linear() * Vector3f(1,1,1); } // FIXXXXXXX + inline Vector3f GetScale() const { + return Vector3f(m_T.linear().col(0).norm(), + m_T.linear().col(1).norm(), + m_T.linear().col(2).norm()); + } inline void Rotate(const Matrix3f m) { this->m_T.rotate(m); } diff --git a/src/Math/testing/ContainerBoxTest.cpp b/src/Math/testing/ContainerBoxTest.cpp index 2c0b391..dec6b09 100644 --- a/src/Math/testing/ContainerBoxTest.cpp +++ b/src/Math/testing/ContainerBoxTest.cpp @@ -122,7 +122,7 @@ int main() { ContainerBox Box; - Box.SetPosition(Vector3f(1,1,1)); + Box.SetOrigin(Vector3f(1,1,1)); Box.SetSize(Vector3f(2,2,2)); Box.EulerYZYRotate(Vector3f(0,0,0)); HPoint3f pt = Box.GetLocalPoint(HPoint3f(2,3,2)); @@ -130,6 +130,39 @@ int main() TEST0( Vector4f0(wp - HPoint3f(2,3,2)) ); } + { + ContainerBox Box; + Box.SetPosition(Vector3f(-0.5,-0.5,-0.5)); + Box.SetSize(Vector3f(1,1,1)); + HPoint3f pt = Box.GetLocalPoint(HPoint3f(0,0,0)); + HPoint3f wp = Box.GetWorldPoint(pt); + TEST0( Vector4f0(wp - HPoint3f(0,0,0)) ); + } + + { + ContainerBox Box; + Box.SetOrigin(Vector3f(-1.5,-1.5,-1.5)); + Box.SetSize(Vector3f(3,3,3)); + HPoint3f pt = Box.GetLocalPoint(HPoint3f(0,0,0)); + HPoint3f wp = Box.GetWorldPoint(pt); + TEST0( Vector4f0(wp - HPoint3f(0,0,0)) ); + + pt = HPoint3f(1,1,1); + wp = Box.GetWorldPoint(pt); + TEST0( Vector4f0(wp - HPoint3f(1.5, 1.5, 1.5)) ); + + Box.EulerYZYRotate(Vector3f(M_PI,0,0)); + pt = HPoint3f(1,1,1); + wp = Box.GetWorldPoint(pt); + TEST0( Vector4f0(wp - HPoint3f(-1.5, 1.5, -1.5)) ); + + } + + + + + + END_TESTING; } diff --git a/src/Math/testing/testing-prototype.h b/src/Math/testing/testing-prototype.h index 5bfc9dd..2e46ed8 100644 --- a/src/Math/testing/testing-prototype.h +++ b/src/Math/testing/testing-prototype.h @@ -31,8 +31,8 @@ static int _fail = 0; \ printf("..:: Testing " #name " ::..\n"); -#define TEST1(val) _fail += (val)==0 -#define TEST0(val) _fail += (val)!=0 +#define TEST1(val) if ((val)==0) { printf("Assertion failed: %s != 0\n", #val); _fail++; } +#define TEST0(val) if ((val)!=0) { printf("Assertion failed: %s != 0\n", #val); _fail++; } #define END_TESTING return _fail; #define ASSERT_EQ(a, b) if ((a) != (b)) { printf("Assertion failed: %s != %s\n", #a, #b); _fail++; } \ No newline at end of file diff --git a/src/Vtk/CMakeLists.txt b/src/Vtk/CMakeLists.txt index cbc9bba..2254b96 100644 --- a/src/Vtk/CMakeLists.txt +++ b/src/Vtk/CMakeLists.txt @@ -3,6 +3,7 @@ set(HEADERS uLibVtkInterface.h vtkContainerBox.h vtkHandlerWidget.h vtkQViewport.h + vtkPolydata.h ) set(SOURCES uLibVtkInterface.cxx @@ -10,6 +11,7 @@ set(SOURCES uLibVtkInterface.cxx vtkContainerBox.cpp vtkHandlerWidget.cpp vtkQViewport.cpp + vtkPolydata.cpp ) ## Pull in Math VTK wrappers (sets MATH_SOURCES / MATH_HEADERS) diff --git a/src/Vtk/HEP/Detectors/vtkMuonEvent.h b/src/Vtk/HEP/Detectors/vtkMuonEvent.h index 9f6d329..0c3abf6 100644 --- a/src/Vtk/HEP/Detectors/vtkMuonEvent.h +++ b/src/Vtk/HEP/Detectors/vtkMuonEvent.h @@ -45,6 +45,7 @@ #include "HEP/Detectors/MuonEvent.h" #include "Vtk/uLibVtkInterface.h" +#include "Vtk/vtkPolydata.h" namespace uLib { namespace Vtk { diff --git a/src/Vtk/HEP/Detectors/vtkMuonScatter.h b/src/Vtk/HEP/Detectors/vtkMuonScatter.h index 10393a1..ac7c85f 100644 --- a/src/Vtk/HEP/Detectors/vtkMuonScatter.h +++ b/src/Vtk/HEP/Detectors/vtkMuonScatter.h @@ -46,6 +46,7 @@ #include "HEP/Detectors/MuonScatter.h" #include "Vtk/uLibVtkInterface.h" +#include "Vtk/vtkPolydata.h" class vtkRenderWindowInteractor; diff --git a/src/Vtk/HEP/Geant/testing/vtkGeantEventTest.cpp b/src/Vtk/HEP/Geant/testing/vtkGeantEventTest.cpp index 337429e..6ef69af 100644 --- a/src/Vtk/HEP/Geant/testing/vtkGeantEventTest.cpp +++ b/src/Vtk/HEP/Geant/testing/vtkGeantEventTest.cpp @@ -15,6 +15,7 @@ #include "HEP/Geant/EmitterPrimary.hh" #include "Math/ContainerBox.h" #include "Math/Dense.h" +#include "Math/Units.h" #include "Vtk/uLibVtkViewer.h" #include "Vtk/HEP/Geant/vtkGeantEvent.h" #include "Vtk/vtkContainerBox.h" @@ -38,24 +39,21 @@ using namespace uLib; class RandomEmitter : public Geant::EmitterPrimary { public: virtual void GeneratePrimaries(G4Event* anEvent) override { - // Start from a random point in a square at z = 5m - // Note: unit in Geant4 is mm by default if not specified, - // but here we use the constants from G4SystemOfUnits.h - double x = (G4UniformRand() - 0.5) * 2000; // -1m to 1m - double y = (G4UniformRand() - 0.5) * 2000; // -1m to 1m - double z = 5000.0; // 5m above origin + // Start from a random point on the top face of the world box (z = 15m) + double x = 0_m; + double y = 0_m; + double z = 14.9_m; // Top face fParticleGun->SetParticlePosition(G4ThreeVector(x, y, z)); - // Aim at a random point within the cube (which is at origin, size 1m) - double tx = (G4UniformRand() - 0.5) * 1000; - double ty = (G4UniformRand() - 0.5) * 1000; - double tz = (G4UniformRand() - 0.5) * 1000; + // Aim at a random point on the bottom face (z = -15m) + double tx = 0_m; + double ty = 0_m; + double tz = -14.9_m; // Bottom face G4ThreeVector dir(tx - x, ty - y, tz - z); - dir = dir.unit(); - fParticleGun->SetParticleMomentumDirection(dir); - + fParticleGun->SetParticleMomentumDirection(dir.unit()); + fParticleGun->SetParticleEnergy(15_GeV); fParticleGun->GeneratePrimaryVertex(anEvent); } }; @@ -91,16 +89,20 @@ void KeyPressCallbackFunction(vtkObject* caller, long unsigned int eventId, void state->viewer->GetRenderer()->Render(); state->viewer->GetRenderWindow()->Render(); } + else { + std::cout << " No event collected." << std::endl; + } } } int main(int argc, char** argv) { // 1. Setup Geant4 Scene - ContainerBox world_box(Vector3f(30000, 30000, 30000)); Geant::Scene scene; - scene.ConstructWorldBox(&world_box, "G4_AIR"); + scene.ConstructWorldBox(Vector3f(30_m, 30_m, 30_m), "G4_AIR"); - ContainerBox iron_box(Vector3f(1000, 1000, 1000)); + ContainerBox iron_box; + iron_box.Scale(Vector3f(10_m, 10_m, 10_m)); + iron_box.SetPosition(Vector3f(-5_m, -5_m, -5_m)); Geant::BoxSolid* iron_cube = new Geant::BoxSolid("IronCube", &iron_box); iron_cube->SetNistMaterial("G4_Fe"); iron_cube->Update(); @@ -114,8 +116,16 @@ int main(int argc, char** argv) { Vtk::Viewer viewer; viewer.GetRenderer()->SetBackground(0.05, 0.05, 0.1); + // Visualize world box + Vtk::vtkContainerBox* vtkWorld = new Vtk::vtkContainerBox(scene.GetWorldBox()); + // vtkWorld->ShowBoundingBox(true); + vtkWorld->ShowScaleMeasures(true); + viewer.AddPuppet(*vtkWorld); + // Visualize iron cube Vtk::vtkContainerBox* vtkIron = new Vtk::vtkContainerBox(&iron_box); + vtkIron->SetOpacity(0.2); + vtkIron->SetRepresentation(Vtk::Puppet::Surface); viewer.AddPuppet(*vtkIron); // 3. Event Handling diff --git a/src/Vtk/HEP/Geant/vtkGeantEvent.cpp b/src/Vtk/HEP/Geant/vtkGeantEvent.cpp index 070d704..d70c7a7 100644 --- a/src/Vtk/HEP/Geant/vtkGeantEvent.cpp +++ b/src/Vtk/HEP/Geant/vtkGeantEvent.cpp @@ -97,8 +97,8 @@ void vtkGeantEvent::InstallPipe() { m_MuonPath->SetMapper(mapper); // Set default look: Red line - m_MuonPath->GetProperty()->SetColor(1.0, 0.0, 0.0); - m_MuonPath->GetProperty()->SetLineWidth(2.0); + m_MuonPath->GetProperty()->SetColor(1.0, 1.0, 1.0); + m_MuonPath->GetProperty()->SetLineWidth(1.0); m_MuonPath->GetProperty()->SetAmbient(1.0); this->SetProp(m_MuonPath); diff --git a/src/Vtk/HEP/Geant/vtkGeantEvent.h b/src/Vtk/HEP/Geant/vtkGeantEvent.h index c9feebd..c904baf 100644 --- a/src/Vtk/HEP/Geant/vtkGeantEvent.h +++ b/src/Vtk/HEP/Geant/vtkGeantEvent.h @@ -28,6 +28,7 @@ #include "HEP/Geant/GeantEvent.h" #include "uLibVtkInterface.h" +#include "vtkPolydata.h" #include namespace uLib { diff --git a/src/Vtk/Math/vtkQuadMesh.h b/src/Vtk/Math/vtkQuadMesh.h index b238ed7..a49e2ab 100644 --- a/src/Vtk/Math/vtkQuadMesh.h +++ b/src/Vtk/Math/vtkQuadMesh.h @@ -28,6 +28,7 @@ #include "Math/QuadMesh.h" #include "Vtk/uLibVtkInterface.h" +#include "Vtk/vtkPolydata.h" class vtkPolyData; class vtkActor; diff --git a/src/Vtk/Math/vtkTriangleMesh.h b/src/Vtk/Math/vtkTriangleMesh.h index 8899838..8920c1e 100644 --- a/src/Vtk/Math/vtkTriangleMesh.h +++ b/src/Vtk/Math/vtkTriangleMesh.h @@ -28,6 +28,7 @@ #include "Math/TriangleMesh.h" #include "Vtk/uLibVtkInterface.h" +#include "Vtk/vtkPolydata.h" class vtkPolyData; class vtkActor; diff --git a/src/Vtk/uLibVtkInterface.cxx b/src/Vtk/uLibVtkInterface.cxx index 47ca752..186a024 100644 --- a/src/Vtk/uLibVtkInterface.cxx +++ b/src/Vtk/uLibVtkInterface.cxx @@ -37,6 +37,7 @@ #endif +#include #include #include #include @@ -46,10 +47,13 @@ #include #include #include +#include +#include +#include +#include +#include +#include -#include -#include -#include #include "uLibVtkInterface.h" @@ -64,46 +68,6 @@ namespace Vtk { -//////////////////////////////////////////////////////////////////////////////// -//////////////////////////////////////////////////////////////////////////////// -//////////////////////////////////////////////////////////////////////////////// -// POLYDATA // - - -void Polydata::SaveToFile(const char *vtk_file) -{ - vtkSmartPointer writer = - vtkSmartPointer::New(); - writer->SetFileName(vtk_file); - vtkPolyData * data = GetPolyData(); - if(data) { -# if VTK_MAJOR_VERSION <= 5 - writer->SetInputConnection(data->GetProducerPort()); -# else - writer->SetInputData(data); -# endif - writer->Update(); - } -} - -void Polydata::SaveToXMLFile(const char *vtp_file) -{ - vtkSmartPointer writer = - vtkSmartPointer::New(); - writer->SetFileName(vtp_file); - vtkPolyData * data = GetPolyData(); - if(data) { -# if VTK_MAJOR_VERSION <= 5 - writer->SetInputConnection(data->GetProducerPort()); -# else - writer->SetInputData(data); -# endif - writer->Update(); - } -} - - - //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// @@ -115,19 +79,59 @@ class PuppetData { public: PuppetData() : m_Renderers(vtkRendererCollection::New()), - m_Assembly(vtkPropAssembly::New()) - {} + m_Assembly(vtkPropAssembly::New()), + m_OutlineSource(NULL), + m_OutlineActor(NULL), + m_CubeAxesActor(NULL), + m_ShowBoundingBox(false), + m_ShowScaleMeasures(false), + m_Representation(-1), + m_Opacity(-1.0) + { + m_Color[0] = m_Color[1] = m_Color[2] = -1.0; + } ~PuppetData() { m_Renderers->RemoveAllItems(); m_Assembly->GetParts()->RemoveAllItems(); m_Renderers->Delete(); m_Assembly->Delete(); + if (m_OutlineSource) m_OutlineSource->Delete(); + if (m_OutlineActor) m_OutlineActor->Delete(); + if (m_CubeAxesActor) m_CubeAxesActor->Delete(); } // members // vtkRendererCollection *m_Renderers; vtkPropAssembly *m_Assembly; + + vtkOutlineSource *m_OutlineSource; + vtkActor *m_OutlineActor; + vtkCubeAxesActor *m_CubeAxesActor; + + bool m_ShowBoundingBox; + bool m_ShowScaleMeasures; + + int m_Representation; + double m_Color[3]; + double m_Opacity; + + void ApplyAppearance(vtkProp *p) { + vtkActor *actor = vtkActor::SafeDownCast(p); + if (!actor) return; + + if (m_Representation != -1) { + actor->GetProperty()->SetRepresentation(m_Representation); + } + + if (m_Color[0] != -1.0) { + actor->GetProperty()->SetColor(m_Color); + } + + if (m_Opacity != -1.0) { + actor->GetProperty()->SetOpacity(m_Opacity); + } + } }; // -------------------------------------------------------------------------- // @@ -155,7 +159,10 @@ vtkProp *Puppet::GetProp() void Puppet::SetProp(vtkProp *prop) { - if(prop) d->m_Assembly->AddPart(prop); + if(prop) { + d->m_Assembly->AddPart(prop); + d->ApplyAppearance(prop); + } } void Puppet::RemoveProp(vtkProp *prop) @@ -182,6 +189,12 @@ void Puppet::ConnectRenderer(vtkRenderer *renderer) for (int i=0; iGetNumberOfItems(); ++i) renderer->AddActor(props->GetNextProp()); } + + if (d->m_ShowBoundingBox && d->m_OutlineActor) renderer->AddActor(d->m_OutlineActor); + if (d->m_ShowScaleMeasures && d->m_CubeAxesActor) { + d->m_CubeAxesActor->SetCamera(renderer->GetActiveCamera()); + renderer->AddActor(d->m_CubeAxesActor); + } } } @@ -197,6 +210,10 @@ void Puppet::DisconnectRenderer(vtkRenderer *renderer) for (int i=0; iGetNumberOfItems(); ++i) renderer->RemoveViewProp(props->GetNextProp()); } + + if (d->m_ShowBoundingBox && d->m_OutlineActor) renderer->RemoveActor(d->m_OutlineActor); + if (d->m_ShowScaleMeasures && d->m_CubeAxesActor) renderer->RemoveActor(d->m_CubeAxesActor); + this->GetRenderers()->RemoveItem(renderer); } } @@ -225,6 +242,121 @@ void Puppet::PrintSelf(std::ostream &o) const d->m_Renderers->PrintSelf(o,vtkIndent(1)); } +void Puppet::ShowBoundingBox(bool show) +{ + if (d->m_ShowBoundingBox == show) return; + d->m_ShowBoundingBox = show; + if (show) { + if (!d->m_OutlineActor) { + d->m_OutlineSource = vtkOutlineSource::New(); + d->m_OutlineActor = vtkActor::New(); + vtkSmartPointer mapper = vtkSmartPointer::New(); + mapper->SetInputConnection(d->m_OutlineSource->GetOutputPort()); + d->m_OutlineActor->SetMapper(mapper); + d->m_OutlineActor->GetProperty()->SetColor(1.0, 1.0, 1.0); + } + + double* bounds = d->m_Assembly->GetBounds(); + d->m_OutlineSource->SetBounds(bounds); + d->m_OutlineSource->Update(); + + d->m_Renderers->InitTraversal(); + for (int i = 0; i < d->m_Renderers->GetNumberOfItems(); ++i) { + vtkRenderer *renderer = d->m_Renderers->GetNextItem(); + renderer->AddActor(d->m_OutlineActor); + } + } else { + if (d->m_OutlineActor) { + d->m_Renderers->InitTraversal(); + for (int i = 0; i < d->m_Renderers->GetNumberOfItems(); ++i) { + vtkRenderer *renderer = d->m_Renderers->GetNextItem(); + renderer->RemoveActor(d->m_OutlineActor); + } + } + } +} + +void Puppet::ShowScaleMeasures(bool show) +{ + if (d->m_ShowScaleMeasures == show) return; + d->m_ShowScaleMeasures = show; + if (show) { + if (!d->m_CubeAxesActor) { + d->m_CubeAxesActor = vtkCubeAxesActor::New(); + d->m_CubeAxesActor->SetFlyModeToOuterEdges(); + d->m_CubeAxesActor->SetUseTextActor3D(1); + d->m_CubeAxesActor->GetProperty()->SetColor(1.0, 1.0, 1.0); + } + + double* bounds = d->m_Assembly->GetBounds(); + d->m_CubeAxesActor->SetBounds(bounds); + + d->m_Renderers->InitTraversal(); + for (int i = 0; i < d->m_Renderers->GetNumberOfItems(); ++i) { + vtkRenderer *renderer = d->m_Renderers->GetNextItem(); + d->m_CubeAxesActor->SetCamera(renderer->GetActiveCamera()); + renderer->AddActor(d->m_CubeAxesActor); + } + } else { + if (d->m_CubeAxesActor) { + d->m_Renderers->InitTraversal(); + for (int i = 0; i < d->m_Renderers->GetNumberOfItems(); ++i) { + vtkRenderer *renderer = d->m_Renderers->GetNextItem(); + renderer->RemoveActor(d->m_CubeAxesActor); + } + } + } +} + +void Puppet::SetRepresentation(Representation mode) +{ + int rep = VTK_SURFACE; + switch (mode) { + case Points: rep = VTK_POINTS; break; + case Wireframe: rep = VTK_WIREFRAME; break; + case Surface: rep = VTK_SURFACE; break; + } + d->m_Representation = rep; + + vtkPropCollection *props = d->m_Assembly->GetParts(); + props->InitTraversal(); + for (int i = 0; i < props->GetNumberOfItems(); ++i) { + d->ApplyAppearance(props->GetNextProp()); + } +} + +void Puppet::SetRepresentation(const char *mode) +{ + std::string s(mode); + if (s == "points") SetRepresentation(Points); + else if (s == "wireframe") SetRepresentation(Wireframe); + else if (s == "shaded" || s == "surface") SetRepresentation(Surface); +} + +void Puppet::SetColor(double r, double g, double b) +{ + d->m_Color[0] = r; + d->m_Color[1] = g; + d->m_Color[2] = b; + + vtkPropCollection *props = d->m_Assembly->GetParts(); + props->InitTraversal(); + for (int i = 0; i < props->GetNumberOfItems(); ++i) { + d->ApplyAppearance(props->GetNextProp()); + } +} + +void Puppet::SetOpacity(double alpha) +{ + d->m_Opacity = alpha; + + vtkPropCollection *props = d->m_Assembly->GetParts(); + props->InitTraversal(); + for (int i = 0; i < props->GetNumberOfItems(); ++i) { + d->ApplyAppearance(props->GetNextProp()); + } +} + diff --git a/src/Vtk/uLibVtkInterface.h b/src/Vtk/uLibVtkInterface.h index dd473e1..c3127fc 100644 --- a/src/Vtk/uLibVtkInterface.h +++ b/src/Vtk/uLibVtkInterface.h @@ -37,6 +37,7 @@ class vtkRenderer; class vtkRendererCollection; class vtkRenderWindowInteractor; + namespace uLib { namespace Vtk { @@ -59,8 +60,28 @@ public: vtkRendererCollection *GetRenderers() const; + enum Representation { + Points, + Wireframe, + Surface + }; + virtual void PrintSelf(std::ostream &o) const; + void ShowBoundingBox(bool show = true); + + void ShowScaleMeasures(bool show = true); + + void SetRepresentation(Representation mode); + + void SetRepresentation(const char *mode); + + void SetColor(double r, double g, double b); + + void SetOpacity(double alpha); + + + virtual void ConnectInteractor(class vtkRenderWindowInteractor *interactor) { (void)interactor; } @@ -75,18 +96,6 @@ private: class PuppetData *d; }; -class Polydata { -public: - virtual vtkPolyData *GetPolyData() const { return NULL; } - - virtual void SaveToFile(const char *vtk_file); - - virtual void SaveToXMLFile(const char *vtp_file); - -protected: - virtual ~Polydata() {} -}; - } // namespace Vtk } // namespace uLib diff --git a/src/Vtk/uLibVtkViewer.cpp b/src/Vtk/uLibVtkViewer.cpp index b956c64..e25318d 100644 --- a/src/Vtk/uLibVtkViewer.cpp +++ b/src/Vtk/uLibVtkViewer.cpp @@ -56,6 +56,8 @@ public: }; vtkStandardNewMacro(vtkInteractorStyleNoSpin); + + namespace uLib { namespace Vtk { diff --git a/src/Vtk/vtkContainerBox.h b/src/Vtk/vtkContainerBox.h index 64dc701..103175a 100644 --- a/src/Vtk/vtkContainerBox.h +++ b/src/Vtk/vtkContainerBox.h @@ -28,6 +28,7 @@ #include "Math/ContainerBox.h" #include "uLibVtkInterface.h" +#include "vtkPolydata.h" #include namespace uLib { diff --git a/src/Vtk/vtkPolydata.cpp b/src/Vtk/vtkPolydata.cpp new file mode 100644 index 0000000..c5a9860 --- /dev/null +++ b/src/Vtk/vtkPolydata.cpp @@ -0,0 +1,70 @@ +/*////////////////////////////////////////////////////////////////////////////// +// CMT Cosmic Muon Tomography project ////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova + All rights reserved + + Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it > + + ------------------------------------------------------------------ + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 3.0 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library. + +//////////////////////////////////////////////////////////////////////////////*/ + +#include +#include +#include +#include +#include + +#include "vtkPolydata.h" + +namespace uLib { +namespace Vtk { + +void Polydata::SaveToFile(const char *vtk_file) +{ + vtkSmartPointer writer = + vtkSmartPointer::New(); + writer->SetFileName(vtk_file); + vtkPolyData * data = GetPolyData(); + if(data) { +# if VTK_MAJOR_VERSION <= 5 + writer->SetInputConnection(data->GetProducerPort()); +# else + writer->SetInputData(data); +# endif + writer->Update(); + } +} + +void Polydata::SaveToXMLFile(const char *vtp_file) +{ + vtkSmartPointer writer = + vtkSmartPointer::New(); + writer->SetFileName(vtp_file); + vtkPolyData * data = GetPolyData(); + if(data) { +# if VTK_MAJOR_VERSION <= 5 + writer->SetInputConnection(data->GetProducerPort()); +# else + writer->SetInputData(data); +# endif + writer->Update(); + } +} + +} // Vtk +} // uLib diff --git a/src/Vtk/vtkPolydata.h b/src/Vtk/vtkPolydata.h index 5aa5b27..d5f63a6 100644 --- a/src/Vtk/vtkPolydata.h +++ b/src/Vtk/vtkPolydata.h @@ -28,22 +28,25 @@ #ifndef VTKPOLYDATA_H #define VTKPOLYDATA_H -#include "uLibVtkInterface.h" - +class vtkPolyData; namespace uLib { +namespace Vtk { -class vtkPolyData : public Abstract::uLibVtkPolydata { - typedef ::vtkPolyData Content; +class Polydata { public: + virtual vtkPolyData *GetPolyData() const { return NULL; } + virtual void SaveToFile(const char *vtk_file); + virtual void SaveToXMLFile(const char *vtp_file); +protected: + virtual ~Polydata() {} }; - - -} +} // namespace Vtk +} // namespace uLib #endif // VTKPOLYDATA_H