diff --git a/src/HEP/Geant/ActionInitialization.cpp b/src/HEP/Geant/ActionInitialization.cpp index cdbc4b2..2447b39 100644 --- a/src/HEP/Geant/ActionInitialization.cpp +++ b/src/HEP/Geant/ActionInitialization.cpp @@ -1,30 +1,36 @@ #include "ActionInitialization.hh" #include "EmitterPrimary.hh" +#include "SteppingAction.hh" namespace uLib { namespace Geant { -ActionInitialization::ActionInitialization() : G4VUserActionInitialization() {} +ActionInitialization::ActionInitialization(EmitterPrimary *emitter, + Vector *output) + : G4VUserActionInitialization(), + m_Emitter(emitter), + m_Output(output) +{} ActionInitialization::~ActionInitialization() {} void ActionInitialization::BuildForMaster() const { - // Questo metodo viene usato in modalità Multi-Threading. - // Serve per le azioni che devono esistere solo nel thread Master - // (tipicamente solo per inizializzare file di output o il RunAction globale). - - // Esempio: SetUserAction(new RunAction()); + // Master thread: no per-event actions needed } void ActionInitialization::Build() const { - // Questo è il cuore dell'inizializzazione per i thread di lavoro. - // Qui passiamo il nostro generatore di muoni a Geant4. - SetUserAction(new EmitterPrimary()); + // Register the primary generator + if (m_Emitter) { + SetUserAction(m_Emitter); + } else { + // Fallback: default EmitterPrimary + SetUserAction(new EmitterPrimary()); + } - // In una simulazione completa, qui registreresti anche le altre classi: - // SetUserAction(new RunAction()); - // SetUserAction(new EventAction()); - // SetUserAction(new SteppingAction()); + // Register stepping action to collect scattering data + if (m_Output) { + SetUserAction(new SteppingAction(m_Output)); + } } } // namespace Geant diff --git a/src/HEP/Geant/ActionInitialization.hh b/src/HEP/Geant/ActionInitialization.hh index e5e570f..1a872f6 100644 --- a/src/HEP/Geant/ActionInitialization.hh +++ b/src/HEP/Geant/ActionInitialization.hh @@ -2,13 +2,20 @@ #define ActionInitialization_h #include "G4VUserActionInitialization.hh" +#include "Core/Vector.h" namespace uLib { namespace Geant { +class EmitterPrimary; +class GeantEvent; + class ActionInitialization : public G4VUserActionInitialization { public: - ActionInitialization(); + /// @param emitter the primary generator to use (owned by caller) + /// @param output pointer to the results vector (owned by caller) + ActionInitialization(EmitterPrimary *emitter = nullptr, + Vector *output = nullptr); ~ActionInitialization(); // Metodo chiamato solo dal thread principale (Master) @@ -16,6 +23,10 @@ public: // Metodo chiamato dai thread di lavoro (Worker) o in modalità sequenziale virtual void Build() const; + +private: + EmitterPrimary *m_Emitter; + Vector *m_Output; }; } // namespace Geant diff --git a/src/HEP/Geant/CMakeLists.txt b/src/HEP/Geant/CMakeLists.txt index 8c8d852..08d33e1 100644 --- a/src/HEP/Geant/CMakeLists.txt +++ b/src/HEP/Geant/CMakeLists.txt @@ -19,6 +19,7 @@ set(HEADERS DetectorConstruction.hh PhysicsList.hh ActionInitialization.hh + SteppingAction.hh ) set(SOURCES @@ -28,6 +29,7 @@ set(SOURCES DetectorConstruction.cpp PhysicsList.cpp ActionInitialization.cpp + SteppingAction.cpp ) set(libname ${PACKAGE_LIBPREFIX}Geant) diff --git a/src/HEP/Geant/GeantEvent.h b/src/HEP/Geant/GeantEvent.h index 2c28928..9042e7a 100644 --- a/src/HEP/Geant/GeantEvent.h +++ b/src/HEP/Geant/GeantEvent.h @@ -30,34 +30,61 @@ #include "Core/Vector.h" #include "Math/Dense.h" -#include "Detectors/ChamberHitEvent.h" +#include namespace uLib { namespace Geant { -class GeantEventData { +// Forward declaration for friend access +class SteppingAction; + +/////////////////////////////////////////////////////////////////////////////// +/// GeantEvent — output of a Geant4 simulation run. +/// +/// m_Momentum and m_GenVector are set from the EmitterPrimary at generation. +/// During simulation, each scattering interaction deposits a Delta in m_Path, +/// recording the change of momentum and direction at each step boundary. +/////////////////////////////////////////////////////////////////////////////// + +class GeantEvent { + public: - uLibGetMacro(EventID, Id_t) uLibGetMacro(Momentum, Scalarf) - uLibConstRefMacro(GenPos, Vector3f) uLibConstRefMacro(GenDir, Vector3f) - uLibConstRefMacro(ChEvents, Vector) + + /// A single interaction step along the muon path. + struct Delta { + Scalarf m_Length; ///< step length through the solid + Scalarf m_Momentum; ///< momentum magnitude at this step + HVector3f m_Direction; ///< direction after scattering + std::string m_SolidName; ///< name of the solid where interaction occurred + + uLibGetMacro(Length, Scalarf) + uLibGetMacro(Momentum, Scalarf) + uLibConstRefMacro(Direction, HVector3f) + uLibConstRefMacro(SolidName, std::string) + + Delta() : m_Length(0), m_Momentum(0), m_Direction(HVector3f::Zero()) {} + }; + + // --- Read-only accessors (public) --- // + uLibGetMacro(EventID, Id_t) + uLibGetMacro(Momentum, Scalarf) + uLibConstRefMacro(GenVector, HLine3f) + uLibConstRefMacro(Path, Vector) private: - friend class GeantEvent; Id_t m_EventID; Scalarf m_Momentum; - Vector3f m_GenPos; - Vector3f m_GenDir; - Vector m_ChEvents; -}; + HLine3f m_GenVector; + Vector m_Path; -class GeantEvent : public GeantEventData { public: - uLibSetMacro(EventID, Id_t) uLibSetMacro(Momentum, Scalarf) - uLibRefMacro(GenPos, Vector3f) uLibRefMacro(GenDir, Vector3f) - uLibRefMacro(ChEvents, Vector) + GeantEvent() : m_EventID(0), m_Momentum(0), m_GenVector(HLine3f()), m_Path() {} + + // SteppingAction can populate the private fields during simulation + friend class SteppingAction; }; } // namespace Geant } // namespace uLib -#endif // GEANTEVENT_H +#endif // U_GEANTEVENT_H diff --git a/src/HEP/Geant/Scene.cpp b/src/HEP/Geant/Scene.cpp index 45071c0..db5bf2b 100644 --- a/src/HEP/Geant/Scene.cpp +++ b/src/HEP/Geant/Scene.cpp @@ -54,10 +54,27 @@ private: class SceneImpl *m_Owner; }; +static void CheckGeant4Environment() { + static bool checked = false; + if (checked) return; + checked = true; + if (!std::getenv("G4ENSDFSTATEDATA")) { + std::cerr << "********************************************************" << std::endl; + std::cerr << " WARNING: Geant4 environment variables are not set!" << std::endl; + std::cerr << " Please activate the environment before running:" << std::endl; + std::cerr << " micromamba activate mutom" << std::endl; + std::cerr << "********************************************************" << std::endl; + } +} + class SceneImpl { public: // constructor // - SceneImpl() : m_RunManager(G4RunManagerFactory::CreateRunManager()) {} + SceneImpl() : m_RunManager(G4RunManagerFactory::CreateRunManager(G4RunManagerType::Serial)), + m_Emitter(nullptr), + m_Output(nullptr) { + m_RunManager->SetUserInitialization(new PhysicsList); + } // destructor // ~SceneImpl() { @@ -68,8 +85,8 @@ public: void Initialize() { // Set mandatory initialization classes for Geant4 m_RunManager->SetUserInitialization(new SceneDetectorConstruction(this)); - m_RunManager->SetUserInitialization(new PhysicsList); - m_RunManager->SetUserInitialization(new ActionInitialization); + m_RunManager->SetUserInitialization( + new ActionInitialization(m_Emitter, m_Output)); // Initialize Geant4 m_RunManager->Initialize(); @@ -79,6 +96,8 @@ public: Vector m_Solids; Solid *m_World = nullptr; G4RunManager *m_RunManager; + EmitterPrimary *m_Emitter; + Vector *m_Output; }; SceneDetectorConstruction::SceneDetectorConstruction(SceneImpl *owner) @@ -91,10 +110,10 @@ G4VPhysicalVolume *SceneDetectorConstruction::Construct() { - - - -Scene::Scene() : d(new SceneImpl()) {} +Scene::Scene() { + CheckGeant4Environment(); + d = new SceneImpl(); +} Scene::~Scene() { delete d; } void Scene::AddSolid(Solid *solid, Solid *parent) { @@ -143,9 +162,26 @@ void Scene::ConstructWorldBox(const ContainerBox *box, const char *material) { d->m_World->SetTransform(transform); } +void Scene::SetEmitter(EmitterPrimary *emitter) { + d->m_Emitter = emitter; +} + void Scene::Initialize() { d->Initialize(); } +void Scene::RunSimulation(int nEvents, Vector &results) { + d->m_Output = &results; + + // Re-initialize ActionInitialization with the output buffer + // (ActionInitialization was already set during Initialize, but we need + // to ensure the output pointer is current) + d->m_RunManager->SetUserInitialization( + new ActionInitialization(d->m_Emitter, &results)); + + d->m_RunManager->BeamOn(nEvents); +} + } // namespace Geant } // namespace uLib + diff --git a/src/HEP/Geant/Scene.h b/src/HEP/Geant/Scene.h index 5d031b1..3467dee 100644 --- a/src/HEP/Geant/Scene.h +++ b/src/HEP/Geant/Scene.h @@ -31,12 +31,15 @@ #include "Core/Object.h" #include "Core/Vector.h" #include "Solid.h" +#include "GeantEvent.h" class G4VPhysicalVolume; namespace uLib { namespace Geant { +class EmitterPrimary; + class Scene : public Object { public: Scene(); @@ -46,8 +49,17 @@ public: void ConstructWorldBox(const ContainerBox *box, const char *material); + /// Set the primary generator (emitter) for the simulation. + /// The Scene does NOT take ownership of the emitter. + void SetEmitter(EmitterPrimary *emitter); + + /// Initialize the Geant4 run manager with detector, physics, and action. void Initialize(); + /// Run the simulation for nEvents muons. + /// Results are appended to the provided vector. + void RunSimulation(int nEvents, Vector &results); + private: class SceneImpl *d; }; diff --git a/src/HEP/Geant/Solid.cpp b/src/HEP/Geant/Solid.cpp index c8e8cc9..1e66c90 100644 --- a/src/HEP/Geant/Solid.cpp +++ b/src/HEP/Geant/Solid.cpp @@ -62,16 +62,18 @@ Solid::Solid(const char *name) void Solid::SetNistMaterial(const char *name) { G4NistManager *nist = G4NistManager::Instance(); - m_Material = nist->FindOrBuildMaterial(name); - if (m_Logical) - m_Logical->SetMaterial(m_Material); + G4Material *mat = nist->FindOrBuildMaterial(name); + if (mat) SetMaterial(mat); } void Solid::SetMaterial(G4Material *material) { if (material) { m_Material = material; - if (m_Logical) + if (m_Logical) { m_Logical->SetMaterial(material); + } else if (GetG4Solid()) { + m_Logical = new G4LogicalVolume(GetG4Solid(), m_Material, GetName()); + } } } @@ -149,7 +151,9 @@ void TessellatedSolid::SetMesh(TriangleMesh &mesh) { DetectorsSolidImpl::getG4Vector3f(mesh.Points().at(trg(2))), ABSOLUTE); ts->AddFacet((G4VFacet *)facet); } - this->m_Logical->SetSolid(ts); + if (this->m_Logical) { + this->m_Logical->SetSolid(ts); + } } @@ -160,7 +164,9 @@ BoxSolid::BoxSolid(const char *name, ContainerBox *box) : BaseClass(name) { m_Solid = new G4Box(name, 0.5, 0.5, 0.5); m_Object = box; Object::connect(box, &ContainerBox::Updated, this, &BoxSolid::Update); - this->m_Logical->SetSolid(m_Solid); + if (m_Logical) { + m_Logical->SetSolid(m_Solid); + } } void BoxSolid::Update() { diff --git a/src/HEP/Geant/Solid.h b/src/HEP/Geant/Solid.h index 65d5bfc..3f7a966 100644 --- a/src/HEP/Geant/Solid.h +++ b/src/HEP/Geant/Solid.h @@ -59,6 +59,8 @@ public: uLibGetSetMacro(Logical, G4LogicalVolume *) uLibGetSetMacro(Physical, G4VPhysicalVolume *) + virtual G4VSolid* GetG4Solid() const { return nullptr; } + inline const char *GetName() const { return m_Logical ? m_Logical->GetName().c_str() : m_Name.c_str(); } @@ -83,6 +85,7 @@ public: TessellatedSolid(const char *name); void SetMesh(TriangleMesh &mesh); uLibGetMacro(Solid, G4TessellatedSolid *) + virtual G4VSolid* GetG4Solid() const override { return (G4VSolid*)m_Solid; } public slots: void Update(); @@ -101,6 +104,7 @@ class BoxSolid : public Solid { public: BoxSolid(const char *name, ContainerBox *box); + virtual G4VSolid* GetG4Solid() const override { return (G4VSolid*)m_Solid; } public slots: void Update(); diff --git a/src/HEP/Geant/SteppingAction.cpp b/src/HEP/Geant/SteppingAction.cpp new file mode 100644 index 0000000..a12f125 --- /dev/null +++ b/src/HEP/Geant/SteppingAction.cpp @@ -0,0 +1,95 @@ +#include "SteppingAction.hh" + +#include "G4Step.hh" +#include "G4Track.hh" +#include "G4Event.hh" +#include "G4RunManager.hh" +#include "G4LogicalVolume.hh" +#include "G4SystemOfUnits.hh" +#include "G4ParticleDefinition.hh" + +namespace uLib { +namespace Geant { + +SteppingAction::SteppingAction(Vector *output) + : G4UserSteppingAction(), + m_Output(output), + m_Current(), + m_LastEventID(-1) +{} + +SteppingAction::~SteppingAction() { + // Push the last event if any steps were collected + if (m_LastEventID >= 0 && !m_Current.m_Path.empty()) { + m_Output->push_back(m_Current); + } +} + +void SteppingAction::UserSteppingAction(const G4Step *step) { + if (!step || !m_Output) return; + + const G4Track *track = step->GetTrack(); + if (!track) return; + + // Only record primary particle (muon) + if (track->GetParentID() != 0) return; + + // 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; + + // Step length + delta.m_Length = static_cast(step->GetStepLength() / mm); + + // Post-step momentum + G4ThreeVector postMom = track->GetMomentum(); + delta.m_Momentum = static_cast(postMom.mag() / MeV); + + // Post-step direction + G4ThreeVector dir = track->GetMomentumDirection(); + delta.m_Direction = HVector3f(static_cast(dir.x()), + static_cast(dir.y()), + static_cast(dir.z())); + + // Solid name where the step occurred + const G4LogicalVolume *vol = track->GetVolume() + ? track->GetVolume()->GetLogicalVolume() + : nullptr; + if (vol) { + delta.m_SolidName = vol->GetName(); + } + + m_Current.m_Path.push_back(delta); +} + +} // namespace Geant +} // namespace uLib diff --git a/src/HEP/Geant/SteppingAction.hh b/src/HEP/Geant/SteppingAction.hh new file mode 100644 index 0000000..920113c --- /dev/null +++ b/src/HEP/Geant/SteppingAction.hh @@ -0,0 +1,30 @@ +#ifndef U_GEANT_STEPPINGACTION_HH +#define U_GEANT_STEPPINGACTION_HH + +#include "G4UserSteppingAction.hh" +#include "Core/Vector.h" +#include "GeantEvent.h" + +namespace uLib { +namespace Geant { + +/// SteppingAction collects scattering data at each Geant4 step and +/// builds GeantEvent objects in the output buffer. +class SteppingAction : public G4UserSteppingAction { +public: + /// @param output pointer to the results vector owned by the Scene + SteppingAction(Vector *output); + virtual ~SteppingAction(); + + virtual void UserSteppingAction(const G4Step *step) override; + +private: + Vector *m_Output; ///< destination for finished events + GeantEvent m_Current; ///< event being built + int m_LastEventID; ///< track event transitions +}; + +} // namespace Geant +} // namespace uLib + +#endif diff --git a/src/HEP/Geant/testing/EventTest.cpp b/src/HEP/Geant/testing/EventTest.cpp index 58c268e..e653228 100644 --- a/src/HEP/Geant/testing/EventTest.cpp +++ b/src/HEP/Geant/testing/EventTest.cpp @@ -1,6 +1,10 @@ #include "Geant/Solid.h" #include "HEP/Geant/GeantEvent.h" +#include "HEP/Geant/Scene.h" +#include "HEP/Geant/EmitterPrimary.hh" +#include "Math/ContainerBox.h" #include "Math/TriangleMesh.h" +#include "Math/Dense.h" #include "testing-prototype.h" #include #include @@ -13,59 +17,73 @@ using namespace uLib; int main() { BEGIN_TESTING(Geant Event); - // Test Solid initialization and NIST material // + // Test: Scene with iron cube in air, launch muons, collect events // { - Geant::Solid solid("test_solid"); - TEST1(solid.GetLogical() != nullptr); - - solid.SetNistMaterial("G4_AIR"); - TEST1(solid.GetMaterial() != nullptr); - TEST1(solid.GetMaterial()->GetName() == "G4_AIR"); - } + // 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"); - // Test TessellatedSolid with a simple mesh // - { - Geant::TessellatedSolid tsolid("test_tessellated"); - TEST1(tsolid.GetLogical() != nullptr); - TEST1(tsolid.GetSolid() != nullptr); + // 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); - // 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)); + // 3. Set up emitter (default: mu- at 1 GeV, from z=+10m downward) + Geant::EmitterPrimary *emitter = new Geant::EmitterPrimary(); + scene.SetEmitter(emitter); + // 4. Initialize Geant4 + scene.Initialize(); + // 5. Run simulation: 10 muons + int nEvents = 10; + Vector results; + scene.RunSimulation(nEvents, results); - tsolid.SetMesh(mesh); - // GeantEvent geant_event; - - + // 6. Check results + printf(" Collected %zu events\n", results.size()); + TEST1(results.size() > 0); + + for (size_t i = 0; i < results.size(); ++i) { + const Geant::GeantEvent &ev = results[i]; + bool hitIron = false; + for (const auto &d : ev.Path()) { + if (d.SolidName() == "IronCube") { + hitIron = true; + break; + } + } + + printf(" Event %d: momentum=%.1f MeV, path steps=%zu, hitIron=%s\n", + ev.GetEventID(), + ev.GetMomentum(), + ev.Path().size(), + hitIron ? "YES" : "NO"); + + // Each event should have at least one step + TEST1(ev.Path().size() > 0); + + // Print first few deltas + const auto &path = ev.Path(); + for (size_t j = 0; j < path.size() && j < 5; ++j) { + const auto &d = path[j]; + printf(" Delta[%zu]: solid=%s len=%.2f mm, p=%.1f MeV, " + "dir=(%.3f, %.3f, %.3f)\n", + j, + d.SolidName().c_str(), + d.GetLength(), + d.GetMomentum(), + d.Direction()(0), + d.Direction()(1), + d.Direction()(2)); + } + if (path.size() > 5) { + printf(" ... (%zu more deltas)\n", path.size() - 5); + } + } } END_TESTING diff --git a/src/Vtk/CMakeLists.txt b/src/Vtk/CMakeLists.txt index a8e43f4..cbc9bba 100644 --- a/src/Vtk/CMakeLists.txt +++ b/src/Vtk/CMakeLists.txt @@ -27,11 +27,17 @@ add_subdirectory(HEP/MuonTomography) list(APPEND SOURCES ${HEP_MUONTOMOGRAPHY_SOURCES}) list(APPEND HEADERS ${HEP_MUONTOMOGRAPHY_HEADERS}) +## Pull in HEP/Geant VTK wrappers (sets HEP_GEANT_SOURCES / HEADERS) +add_subdirectory(HEP/Geant) +list(APPEND SOURCES ${HEP_GEANT_SOURCES}) +list(APPEND HEADERS ${HEP_GEANT_HEADERS}) + set(LIBRARIES Eigen3::Eigen ${ROOT_LIBRARIES} ${VTK_LIBRARIES} ${PACKAGE_LIBPREFIX}Math - ${PACKAGE_LIBPREFIX}Detectors) + ${PACKAGE_LIBPREFIX}Detectors + ${PACKAGE_LIBPREFIX}Geant) if(USE_CUDA) find_package(CUDAToolkit REQUIRED) @@ -59,6 +65,7 @@ install(FILES ${HEADERS} DESTINATION ${INSTALL_INC_DIR}/Vtk) install(FILES ${MATH_HEADERS} DESTINATION ${INSTALL_INC_DIR}/Vtk/Math) install(FILES ${HEP_DETECTORS_HEADERS} DESTINATION ${INSTALL_INC_DIR}/Vtk/HEP/Detectors) install(FILES ${HEP_MUONTOMOGRAPHY_HEADERS} DESTINATION ${INSTALL_INC_DIR}/Vtk/HEP/MuonTomography) +install(FILES ${HEP_GEANT_HEADERS} DESTINATION ${INSTALL_INC_DIR}/Vtk/HEP/Geant) if(BUILD_TESTING) include(uLibTargetMacros) diff --git a/src/Vtk/HEP/Geant/CMakeLists.txt b/src/Vtk/HEP/Geant/CMakeLists.txt new file mode 100644 index 0000000..48228fa --- /dev/null +++ b/src/Vtk/HEP/Geant/CMakeLists.txt @@ -0,0 +1,17 @@ + +################################################################################ +##### Vtk/HEP/Geant - VTK wrappers for HEP Geant objects ####################### +################################################################################ + +set(HEP_GEANT_SOURCES + ${CMAKE_CURRENT_SOURCE_DIR}/vtkGeantEvent.cpp + PARENT_SCOPE) + +set(HEP_GEANT_HEADERS + ${CMAKE_CURRENT_SOURCE_DIR}/vtkGeantEvent.h + PARENT_SCOPE) + +if(BUILD_TESTING) + include(uLibTargetMacros) + add_subdirectory(testing) +endif() diff --git a/src/Vtk/HEP/Geant/testing/CMakeLists.txt b/src/Vtk/HEP/Geant/testing/CMakeLists.txt new file mode 100644 index 0000000..8086ed4 --- /dev/null +++ b/src/Vtk/HEP/Geant/testing/CMakeLists.txt @@ -0,0 +1,14 @@ +# TESTS +set(TESTS + vtkGeantEventTest +) + +set(LIBRARIES + ${PACKAGE_LIBPREFIX}Core + ${PACKAGE_LIBPREFIX}Math + ${PACKAGE_LIBPREFIX}Geant + ${PACKAGE_LIBPREFIX}Vtk + Eigen3::Eigen +) + +uLib_add_tests(VtkGeant) diff --git a/src/Vtk/HEP/Geant/testing/vtkGeantEventTest.cpp b/src/Vtk/HEP/Geant/testing/vtkGeantEventTest.cpp new file mode 100644 index 0000000..337429e --- /dev/null +++ b/src/Vtk/HEP/Geant/testing/vtkGeantEventTest.cpp @@ -0,0 +1,140 @@ +/*////////////////////////////////////////////////////////////////////////////// +// CMT Cosmic Muon Tomography project ////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova + All rights reserved + + Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it > + +//////////////////////////////////////////////////////////////////////////////*/ + +#include "Geant/Solid.h" +#include "HEP/Geant/GeantEvent.h" +#include "HEP/Geant/Scene.h" +#include "HEP/Geant/EmitterPrimary.hh" +#include "Math/ContainerBox.h" +#include "Math/Dense.h" +#include "Vtk/uLibVtkViewer.h" +#include "Vtk/HEP/Geant/vtkGeantEvent.h" +#include "Vtk/vtkContainerBox.h" + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +using namespace uLib; + +// Custom emitter to fire random muons towards the cube +class RandomEmitter : public Geant::EmitterPrimary { +public: + virtual void GeneratePrimaries(G4Event* anEvent) override { + // Start from a random point 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 + + 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; + + G4ThreeVector dir(tx - x, ty - y, tz - z); + dir = dir.unit(); + fParticleGun->SetParticleMomentumDirection(dir); + + fParticleGun->GeneratePrimaryVertex(anEvent); + } +}; + +struct AppState { + Geant::Scene* scene; + Vtk::Viewer* viewer; + Vector results; +}; + +void KeyPressCallbackFunction(vtkObject* caller, long unsigned int eventId, void* clientData, void* callData) { + auto* interactor = static_cast(caller); + auto* state = static_cast(clientData); + + std::string key = interactor->GetKeySym(); + if (key == "Return") { + std::cout << "--> Firing random muon..." << std::endl; + + // Run one event + state->scene->RunSimulation(1, state->results); + + if (!state->results.empty()) { + // Get the last event + Geant::GeantEvent* lastEvent = &state->results.back(); + std::cout << " Collected event " << lastEvent->GetEventID() + << " with " << lastEvent->Path().size() << " steps." << std::endl; + + // Wrap it for VTK + Vtk::vtkGeantEvent* vtkEvent = new Vtk::vtkGeantEvent(lastEvent); + state->viewer->AddPuppet(*vtkEvent); + + // Re-render + state->viewer->GetRenderer()->Render(); + state->viewer->GetRenderWindow()->Render(); + } + } +} + +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"); + + ContainerBox iron_box(Vector3f(1000, 1000, 1000)); + Geant::BoxSolid* iron_cube = new Geant::BoxSolid("IronCube", &iron_box); + iron_cube->SetNistMaterial("G4_Fe"); + iron_cube->Update(); + scene.AddSolid(iron_cube); + + RandomEmitter* emitter = new RandomEmitter(); + scene.SetEmitter(emitter); + scene.Initialize(); + + // 2. Setup VTK Viewer + Vtk::Viewer viewer; + viewer.GetRenderer()->SetBackground(0.05, 0.05, 0.1); + + // Visualize iron cube + Vtk::vtkContainerBox* vtkIron = new Vtk::vtkContainerBox(&iron_box); + viewer.AddPuppet(*vtkIron); + + // 3. Event Handling + AppState state = { &scene, &viewer, {} }; + + vtkSmartPointer keyCallback = vtkSmartPointer::New(); + keyCallback->SetCallback(KeyPressCallbackFunction); + keyCallback->SetClientData(&state); + + viewer.GetInteractor()->AddObserver(vtkCommand::KeyPressEvent, keyCallback); + + std::cout << "=================================================" << std::endl; + std::cout << " Geant Muon Simulation Viewer" << std::endl; + std::cout << " Press [ENTER] to fire a new random muon" << std::endl; + std::cout << " Press [q] to exit" << std::endl; + std::cout << "=================================================" << std::endl; + + viewer.ZoomAuto(); + viewer.Start(); + + return 0; +} diff --git a/src/Vtk/HEP/Geant/vtkGeantEvent.cpp b/src/Vtk/HEP/Geant/vtkGeantEvent.cpp new file mode 100644 index 0000000..070d704 --- /dev/null +++ b/src/Vtk/HEP/Geant/vtkGeantEvent.cpp @@ -0,0 +1,108 @@ +/*////////////////////////////////////////////////////////////////////////////// +// CMT Cosmic Muon Tomography project ////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + + Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova + All rights reserved + + Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it > + + ------------------------------------------------------------------ + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 3.0 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library. + +//////////////////////////////////////////////////////////////////////////////*/ + +#include "vtkGeantEvent.h" + +#include +#include +#include +#include +#include +#include + +namespace uLib { +namespace Vtk { + +vtkGeantEvent::vtkGeantEvent(Content *content) + : m_MuonPath(vtkActor::New()), m_Content(content) { + this->InstallPipe(); + this->Update(); +} + +vtkGeantEvent::~vtkGeantEvent() { + m_MuonPath->Delete(); +} + +vtkPolyData *vtkGeantEvent::GetPolyData() const { + if (!m_MuonPath || !m_MuonPath->GetMapper()) + return NULL; + return vtkPolyData::SafeDownCast(m_MuonPath->GetMapper()->GetInput()); +} + +void vtkGeantEvent::Update() { + if (!m_Content) + return; + + vtkSmartPointer points = vtkSmartPointer::New(); + vtkSmartPointer lines = vtkSmartPointer::New(); + + const HLine3f &gen = m_Content->GenVector(); + const Vector &path = m_Content->Path(); + + // Start point from generation vector + Vector3f currentPos(gen.origin(0), gen.origin(1), gen.origin(2)); + points->InsertNextPoint(currentPos(0), currentPos(1), currentPos(2)); + + for (size_t i = 0; i < path.size(); ++i) { + const Geant::GeantEvent::Delta &delta = path[i]; + + // P_{i+1} = P_i + Length * Direction + // Note: HVector3f is stored as (x,y,z,0) in HPoint3f template + Vector3f dir(delta.Direction()(0), delta.Direction()(1), delta.Direction()(2)); + currentPos += delta.GetLength() * dir; + + points->InsertNextPoint(currentPos(0), currentPos(1), currentPos(2)); + + vtkIdType line[2] = {static_cast(i), + static_cast(i + 1)}; + lines->InsertNextCell(2, line); + } + + vtkPolyData *polyData = GetPolyData(); + if (polyData) { + polyData->SetPoints(points); + polyData->SetLines(lines); + polyData->Modified(); + } +} + +void vtkGeantEvent::InstallPipe() { + vtkSmartPointer polyData = vtkSmartPointer::New(); + vtkSmartPointer mapper = + vtkSmartPointer::New(); + mapper->SetInputData(polyData); + + 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()->SetAmbient(1.0); + + this->SetProp(m_MuonPath); +} + +} // namespace Vtk +} // namespace uLib diff --git a/src/Vtk/HEP/Geant/vtkGeantEvent.h b/src/Vtk/HEP/Geant/vtkGeantEvent.h new file mode 100644 index 0000000..c9feebd --- /dev/null +++ b/src/Vtk/HEP/Geant/vtkGeantEvent.h @@ -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 U_VTKGEANTEVENT_H +#define U_VTKGEANTEVENT_H + +#include "HEP/Geant/GeantEvent.h" +#include "uLibVtkInterface.h" +#include + +namespace uLib { +namespace Vtk { + +class vtkGeantEvent : public Puppet, public Polydata { + typedef Geant::GeantEvent Content; + +public: + vtkGeantEvent(Content *content); + ~vtkGeantEvent(); + + virtual class vtkPolyData *GetPolyData() const override; + + virtual void Update(); + +protected: + virtual void InstallPipe(); + + vtkActor *m_MuonPath; + Content *m_Content; +}; + +} // namespace Vtk +} // namespace uLib + +#endif // U_VTKGEANTEVENT_H diff --git a/src/Vtk/uLibVtkInterface.h b/src/Vtk/uLibVtkInterface.h index ae6ff26..dd473e1 100644 --- a/src/Vtk/uLibVtkInterface.h +++ b/src/Vtk/uLibVtkInterface.h @@ -35,6 +35,7 @@ class vtkPolyData; class vtkPropCollection; class vtkRenderer; class vtkRendererCollection; +class vtkRenderWindowInteractor; namespace uLib { namespace Vtk {