From 80952cc706bc05b23b5317a8b381e9a9e57b14e9 Mon Sep 17 00:00:00 2001 From: AndreaRigoni Date: Thu, 19 Mar 2026 16:19:36 +0000 Subject: [PATCH] add cylinder emitter --- src/HEP/Geant/EmitterPrimary.cpp | 61 ++++++++ src/HEP/Geant/EmitterPrimary.hh | 22 +++ src/Vtk/HEP/Geant/testing/CMakeLists.txt | 1 + .../testing/vtkCylinderEmitterPrimaryTest.cpp | 136 ++++++++++++++++++ src/Vtk/HEP/Geant/vtkEmitterPrimary.cpp | 51 +++++++ src/Vtk/HEP/Geant/vtkEmitterPrimary.h | 15 ++ 6 files changed, 286 insertions(+) create mode 100644 src/Vtk/HEP/Geant/testing/vtkCylinderEmitterPrimaryTest.cpp diff --git a/src/HEP/Geant/EmitterPrimary.cpp b/src/HEP/Geant/EmitterPrimary.cpp index 22c949a..8831eba 100644 --- a/src/HEP/Geant/EmitterPrimary.cpp +++ b/src/HEP/Geant/EmitterPrimary.cpp @@ -11,6 +11,7 @@ #include "Randomize.hh" #include "EcoMug.hh" +#include "Math/Cylinder.h" namespace uLib { @@ -126,6 +127,66 @@ void SkyPlaneEmitterPrimary::GeneratePrimaries(G4Event* anEvent) { fParticleGun->GeneratePrimaryVertex(anEvent); } + +// -------------------------------------------------------------------------- // +// CylinderEmitterPrimary using EcoMug + +CylinderEmitterPrimary::CylinderEmitterPrimary() + : EmitterPrimary(), m_EcoMug(new EcoMug()), m_Radius(1000.0), m_Height(1000.0) { + m_EcoMug->SetUseCylinder(); + m_EcoMug->SetCylinderRadius(m_Radius/1000.0); + m_EcoMug->SetCylinderHeight(m_Height/1000.0); + m_EcoMug->SetCylinderCenterPosition({0.0, 0.0, m_Height/2000.0}); +} + +CylinderEmitterPrimary::~CylinderEmitterPrimary() { + delete m_EcoMug; +} + +void CylinderEmitterPrimary::SetRadius(float r) { + m_Radius = r; + m_EcoMug->SetCylinderRadius(m_Radius/1000.0); + this->Updated(); +} + +void CylinderEmitterPrimary::SetHeight(float h) { + m_Height = h; + m_EcoMug->SetCylinderHeight(m_Height/1000.0); + m_EcoMug->SetCylinderCenterPosition({0.0, 0.0, m_Height/2000.0}); + this->Updated(); +} + +void CylinderEmitterPrimary::GeneratePrimaries(G4Event* anEvent) { + if (!m_EcoMug) return; + m_EcoMug->Generate(); + + std::array pos_m = m_EcoMug->GetGenerationPosition(); + G4ThreeVector local_pos(pos_m[0]*1000.0, pos_m[1]*1000.0, pos_m[2]*1000.0); + + double p_mag = m_EcoMug->GetGenerationMomentum(); + double theta = m_EcoMug->GetGenerationTheta(); + double phi = m_EcoMug->GetGenerationPhi(); + + G4ThreeVector local_dir( + sin(theta) * cos(phi), + sin(theta) * sin(phi), + cos(theta) + ); + + Matrix4f world_mat = this->GetWorldMatrix(); + Vector3f world_pos = (world_mat * Vector4f(local_pos.x(), local_pos.y(), local_pos.z(), 1.0)).head<3>(); + Vector3f world_dir = (world_mat * Vector4f(local_dir.x(), local_dir.y(), local_dir.z(), 0.0)).head<3>().normalized(); + + G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable(); + G4String particleName = (m_EcoMug->GetCharge() > 0) ? "mu+" : "mu-"; + fParticleGun->SetParticleDefinition(particleTable->FindParticle(particleName)); + + fParticleGun->SetParticlePosition(G4ThreeVector(world_pos(0), world_pos(1), world_pos(2))); + fParticleGun->SetParticleMomentumDirection(G4ThreeVector(world_dir(0), world_dir(1), world_dir(2))); + fParticleGun->SetParticleEnergy(p_mag * GeV); + + fParticleGun->GeneratePrimaryVertex(anEvent); +} // -------------------------------------------------------------------------- // diff --git a/src/HEP/Geant/EmitterPrimary.hh b/src/HEP/Geant/EmitterPrimary.hh index 02fe333..a70f6f0 100644 --- a/src/HEP/Geant/EmitterPrimary.hh +++ b/src/HEP/Geant/EmitterPrimary.hh @@ -55,6 +55,28 @@ class SkyPlaneEmitterPrimary : public EmitterPrimary uLib::Vector2f m_Size; }; +class CylinderEmitterPrimary : public EmitterPrimary +{ + public: + CylinderEmitterPrimary(); + virtual ~CylinderEmitterPrimary(); + + virtual void GeneratePrimaries(G4Event*); + + void SetRadius(float r); + float GetRadius() const { return m_Radius; } + void SetHeight(float h); + float GetHeight() const { return m_Height; } + + private: + EcoMug* m_EcoMug; + float m_Radius; + float m_Height; +}; + + + + class QuadMeshEmitterPrimary : public EmitterPrimary diff --git a/src/Vtk/HEP/Geant/testing/CMakeLists.txt b/src/Vtk/HEP/Geant/testing/CMakeLists.txt index 62595e2..ee0fd69 100644 --- a/src/Vtk/HEP/Geant/testing/CMakeLists.txt +++ b/src/Vtk/HEP/Geant/testing/CMakeLists.txt @@ -3,6 +3,7 @@ set(TESTS vtkGeantEventTest vtkEmitterPrimaryTest vtkSkyPlaneEmitterPrimaryTest + vtkCylinderEmitterPrimaryTest ) set(LIBRARIES diff --git a/src/Vtk/HEP/Geant/testing/vtkCylinderEmitterPrimaryTest.cpp b/src/Vtk/HEP/Geant/testing/vtkCylinderEmitterPrimaryTest.cpp new file mode 100644 index 0000000..f794b76 --- /dev/null +++ b/src/Vtk/HEP/Geant/testing/vtkCylinderEmitterPrimaryTest.cpp @@ -0,0 +1,136 @@ +/*////////////////////////////////////////////////////////////////////////////// +// 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/Cylinder.h" +#include "Math/Dense.h" +#include "Math/Units.h" +#include "Vtk/uLibVtkViewer.h" +#include "Vtk/HEP/Geant/vtkGeantEvent.h" +#include "Vtk/HEP/Geant/vtkEmitterPrimary.h" +#include "Vtk/vtkContainerBox.h" + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +using namespace uLib; + +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 muon from cylinder emitter..." << std::endl; + + // Run one event + state->scene->RunSimulation(1, state->results); + + if (!state->results.empty()) { + Geant::GeantEvent* lastEvent = &state->results.back(); + std::cout << " Collected event " << lastEvent->GetEventID() + << " with " << lastEvent->Path().size() << " steps." << std::endl; + + Vtk::vtkGeantEvent* vtkEvent = new Vtk::vtkGeantEvent(lastEvent); + state->viewer->AddPuppet(*vtkEvent); + + 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 + Geant::Scene scene; + scene.ConstructWorldBox(Vector3f(30_m, 30_m, 30_m), "G4_AIR"); + + ContainerBox iron_box; + iron_box.Scale(Vector3f(10_m, 10_m, 10_m)); + iron_box.SetPosition(Vector3f(0, 0, 0)); + Geant::BoxSolid* iron_cube = new Geant::BoxSolid("IronCube", &iron_box); + iron_cube->SetNistMaterial("G4_Fe"); + iron_cube->Update(); + scene.AddSolid(iron_cube); + + // Use CylinderEmitterPrimary + Geant::CylinderEmitterPrimary* emitter = new Geant::CylinderEmitterPrimary(); + emitter->SetPosition(Vector3f(0, -15_m, 0)); // Vertically offset so base is at -10m + emitter->Rotate(-90_deg, Vector3f(1, 0, 0)); + emitter->SetRadius(15_m); + emitter->SetHeight(30_m); + + scene.SetEmitter(emitter); + scene.Initialize(); + + // 2. Setup VTK Viewer + Vtk::Viewer viewer; + viewer.GetRenderer()->SetBackground(0.05, 0.05, 0.1); + + Vtk::vtkContainerBox* vtkWorld = new Vtk::vtkContainerBox(scene.GetWorldBox()); + vtkWorld->ShowScaleMeasures(true); + vtkWorld->SetRepresentation(Vtk::Puppet::Wireframe); + vtkWorld->SetSelectable(false); + viewer.AddPuppet(*vtkWorld); + + Vtk::vtkContainerBox* vtkIron = new Vtk::vtkContainerBox(&iron_box); + vtkIron->SetOpacity(0.2); + vtkIron->SetRepresentation(Vtk::Puppet::Surface); + viewer.AddPuppet(*vtkIron); + + // Use vtkCylinderEmitterPrimary + Vtk::vtkCylinderEmitterPrimary* vtkEmitter = new Vtk::vtkCylinderEmitterPrimary(*emitter); + vtkEmitter->SetSelectable(false); + viewer.AddPuppet(*vtkEmitter); + + // 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 << " Cylinder Emitter Interactive Test" << std::endl; + std::cout << " Visualize the generation volume (blue cylinder)" << std::endl; + std::cout << " Press [ENTER] to fire a muon using EcoMug" << 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/vtkEmitterPrimary.cpp b/src/Vtk/HEP/Geant/vtkEmitterPrimary.cpp index 66e2b89..c074012 100644 --- a/src/Vtk/HEP/Geant/vtkEmitterPrimary.cpp +++ b/src/Vtk/HEP/Geant/vtkEmitterPrimary.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include "Math/vtkDense.h" namespace uLib { @@ -121,6 +122,56 @@ void vtkQuadMeshEmitterPrimary::contentUpdate() { // For now stick with the arrow. In the future visualize the mesh if useful. vtkEmitterPrimary::contentUpdate(); } + +// -------------------------------------------------------------------------- // +// vtkCylinderEmitterPrimary + +vtkCylinderEmitterPrimary::vtkCylinderEmitterPrimary(Geant::CylinderEmitterPrimary &emitter) + : vtkEmitterPrimary(emitter), m_cylinderEmitter(emitter), m_CylinderSource(vtkCylinderSource::New()) { + + // vtkCylinderSource is along Y by default. + // We will update its actual dimensions in contentUpdate(). + m_CylinderSource->SetResolution(64); + + // Rotate it to be along local Z + vtkNew tcyl; + tcyl->RotateX(90.0); + + vtkNew cylFilter; + cylFilter->SetTransform(tcyl); + cylFilter->SetInputConnection(m_CylinderSource->GetOutputPort()); + + vtkNew append; + // Keep the directional arrow from base class + append->AddInputData(vtkPolyData::SafeDownCast(m_Actor->GetMapper()->GetInput())); + append->AddInputConnection(cylFilter->GetOutputPort()); + + vtkNew mapper; + mapper->SetInputConnection(append->GetOutputPort()); + m_Actor->SetMapper(mapper); + + m_Actor->GetProperty()->SetOpacity(0.5); + m_Actor->GetProperty()->SetColor(0.2, 0.6, 0.9); + + this->contentUpdate(); +} + +vtkCylinderEmitterPrimary::~vtkCylinderEmitterPrimary() { + m_CylinderSource->Delete(); +} + +void vtkCylinderEmitterPrimary::contentUpdate() { + float r = m_cylinderEmitter.GetRadius(); + float h = m_cylinderEmitter.GetHeight(); + + m_CylinderSource->SetRadius(r); + m_CylinderSource->SetHeight(h); + // vtkCylinderSource is along Y. To have the base at Y=0, we center it at h/2. + m_CylinderSource->SetCenter(0, h/2.0, 0); + m_CylinderSource->Update(); + + vtkEmitterPrimary::contentUpdate(); +} } // namespace Vtk } // namespace uLib diff --git a/src/Vtk/HEP/Geant/vtkEmitterPrimary.h b/src/Vtk/HEP/Geant/vtkEmitterPrimary.h index 583df5f..b350eb1 100644 --- a/src/Vtk/HEP/Geant/vtkEmitterPrimary.h +++ b/src/Vtk/HEP/Geant/vtkEmitterPrimary.h @@ -9,6 +9,7 @@ class vtkLineSource; class vtkPolyData; class vtkActor; class vtkPlaneSource; +class vtkCylinderSource; namespace uLib { namespace Vtk { @@ -38,6 +39,8 @@ private: vtkPlaneSource *m_PlaneSource; }; + + class vtkQuadMeshEmitterPrimary : public vtkEmitterPrimary { public: vtkQuadMeshEmitterPrimary(Geant::QuadMeshEmitterPrimary &emitter); @@ -47,6 +50,18 @@ private: Geant::QuadMeshEmitterPrimary &m_meshEmitter; }; +class vtkCylinderEmitterPrimary : public vtkEmitterPrimary { +public: + vtkCylinderEmitterPrimary(Geant::CylinderEmitterPrimary &emitter); + virtual ~vtkCylinderEmitterPrimary(); + + virtual void contentUpdate(); + +private: + Geant::CylinderEmitterPrimary &m_cylinderEmitter; + vtkCylinderSource *m_CylinderSource; +}; +