From a8a313e5cf083e5b2adfe3cfe19dc2958e824f93 Mon Sep 17 00:00:00 2001 From: AndreaRigoni Date: Thu, 19 Mar 2026 15:45:48 +0000 Subject: [PATCH] added skyplaneEmitter --- src/HEP/Geant/EmitterPrimary.cpp | 84 ++++++++++++ src/HEP/Geant/EmitterPrimary.hh | 21 +++ src/Vtk/HEP/Geant/testing/CMakeLists.txt | 1 + .../testing/vtkSkyPlaneEmitterPrimaryTest.cpp | 125 ++++++++++++++++++ src/Vtk/HEP/Geant/vtkEmitterPrimary.cpp | 60 ++++++++- src/Vtk/HEP/Geant/vtkEmitterPrimary.h | 29 +++- 6 files changed, 315 insertions(+), 5 deletions(-) create mode 100644 src/Vtk/HEP/Geant/testing/vtkSkyPlaneEmitterPrimaryTest.cpp diff --git a/src/HEP/Geant/EmitterPrimary.cpp b/src/HEP/Geant/EmitterPrimary.cpp index aca33ad..22c949a 100644 --- a/src/HEP/Geant/EmitterPrimary.cpp +++ b/src/HEP/Geant/EmitterPrimary.cpp @@ -10,6 +10,9 @@ #include "G4SystemOfUnits.hh" #include "Randomize.hh" +#include "EcoMug.hh" + + namespace uLib { namespace Geant { @@ -58,6 +61,87 @@ void EmitterPrimary::GeneratePrimaries(G4Event *anEvent) { } // -------------------------------------------------------------------------- // +// SkyPlaneEmitterPrimary using EcoMug + +SkyPlaneEmitterPrimary::SkyPlaneEmitterPrimary() + : EmitterPrimary(), m_EcoMug(new EcoMug()), m_Size(1000.0, 1000.0) { + // Initial configuration for EcoMug in sky mode + m_EcoMug->SetUseSky(); + m_EcoMug->SetSkySize({m_Size(0)/1000.0, m_Size(1)/1000.0}); +} + +SkyPlaneEmitterPrimary::~SkyPlaneEmitterPrimary() { + delete m_EcoMug; +} + +void SkyPlaneEmitterPrimary::SetSkySize(const uLib::Vector2f& size) { + m_Size = size; + // EcoMug units are in meters (m=1), Geant4 units are in mm. + m_EcoMug->SetSkySize({m_Size(0)/1000.0, m_Size(1)/1000.0}); + this->Updated(); +} + +void SkyPlaneEmitterPrimary::SetPlane(const uLib::Vector3f& p0, const uLib::Vector3f& normal) { + this->SetPosition(p0); + // Orient the emitter so that local Z is the normal. + // This is useful for unconventional planes, though EcoMug sky is usually horizontal. + // If we want a truly 'sky' plane, it usually stays horizontal. + this->Updated(); +} + +void SkyPlaneEmitterPrimary::GeneratePrimaries(G4Event* anEvent) { + if (!m_EcoMug) return; + m_EcoMug->Generate(); + + // EcoMug position is relative to its internal sky center in meters. + // Our wrapper uses the AffineTransform for the overall positioning. + 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); + + // EcoMug momentum (direction and magnitude in GeV/c) + double p_mag = m_EcoMug->GetGenerationMomentum(); + double theta = m_EcoMug->GetGenerationTheta(); + double phi = m_EcoMug->GetGenerationPhi(); + + // EcoMug theta is generated in a way that PI means pointing down (-Z) + G4ThreeVector local_dir( + sin(theta) * cos(phi), + sin(theta) * sin(phi), + cos(theta) + ); + + // Transform local coordinates to world + 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(); + + // Set particle charge + 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); +} + +// -------------------------------------------------------------------------- // + + + + + + + + + + + + +// -------------------------------------------------------------------------- // + QuadMeshEmitterPrimary::QuadMeshEmitterPrimary() : EmitterPrimary(), m_Mesh(nullptr), m_TotalArea(0.0) { diff --git a/src/HEP/Geant/EmitterPrimary.hh b/src/HEP/Geant/EmitterPrimary.hh index f8d1fd3..02fe333 100644 --- a/src/HEP/Geant/EmitterPrimary.hh +++ b/src/HEP/Geant/EmitterPrimary.hh @@ -4,10 +4,13 @@ #include "G4VUserPrimaryGeneratorAction.hh" #include "globals.hh" + namespace uLib { class QuadMesh; } +class EcoMug; + #include "Math/QuadMesh.h" #include "Core/Object.h" #include "Math/Transform.h" @@ -16,6 +19,7 @@ class QuadMesh; class G4ParticleGun; class G4Event; + namespace uLib { namespace Geant { @@ -34,6 +38,23 @@ class EmitterPrimary : public G4VUserPrimaryGeneratorAction, public Object, publ G4ParticleGun* fParticleGun; // Puntatore al cannone di particelle }; +class SkyPlaneEmitterPrimary : public EmitterPrimary +{ + public: + SkyPlaneEmitterPrimary(); + virtual ~SkyPlaneEmitterPrimary(); + + virtual void GeneratePrimaries(G4Event*); + + void SetPlane(const uLib::Vector3f& p0, const uLib::Vector3f& normal); + void SetSkySize(const uLib::Vector2f& size); + uLib::Vector2f GetSkySize() const { return m_Size; } + + private: + EcoMug* m_EcoMug; + uLib::Vector2f m_Size; +}; + class QuadMeshEmitterPrimary : public EmitterPrimary diff --git a/src/Vtk/HEP/Geant/testing/CMakeLists.txt b/src/Vtk/HEP/Geant/testing/CMakeLists.txt index 6efcaaa..62595e2 100644 --- a/src/Vtk/HEP/Geant/testing/CMakeLists.txt +++ b/src/Vtk/HEP/Geant/testing/CMakeLists.txt @@ -2,6 +2,7 @@ set(TESTS vtkGeantEventTest vtkEmitterPrimaryTest + vtkSkyPlaneEmitterPrimaryTest ) set(LIBRARIES diff --git a/src/Vtk/HEP/Geant/testing/vtkSkyPlaneEmitterPrimaryTest.cpp b/src/Vtk/HEP/Geant/testing/vtkSkyPlaneEmitterPrimaryTest.cpp new file mode 100644 index 0000000..3a786c9 --- /dev/null +++ b/src/Vtk/HEP/Geant/testing/vtkSkyPlaneEmitterPrimaryTest.cpp @@ -0,0 +1,125 @@ +#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 "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 sky plane 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 SkyPlaneEmitterPrimary instead of EmitterPrimary + Geant::SkyPlaneEmitterPrimary* emitter = new Geant::SkyPlaneEmitterPrimary(); + emitter->SetPosition(Vector3f(0, 14.9_m, 0)); + emitter->Rotate(-90_deg, Vector3f(1, 0, 0)); + emitter->SetSkySize(Vector2f(20_m, 20_m)); // 10m x 10m plane + + 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 vtkSkyPlaneEmitterPrimary instead of vtkEmitterPrimary + Vtk::vtkSkyPlaneEmitterPrimary* vtkEmitter = new Vtk::vtkSkyPlaneEmitterPrimary(*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 << " Sky Plane Emitter Interactive Test" << std::endl; + std::cout << " Use the Widget to move/rotate the Plane (Emitter)" << std::endl; + std::cout << " Visualize the generation area (blue plane)" << 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 c7a3492..66e2b89 100644 --- a/src/Vtk/HEP/Geant/vtkEmitterPrimary.cpp +++ b/src/Vtk/HEP/Geant/vtkEmitterPrimary.cpp @@ -6,18 +6,24 @@ #include #include #include +#include +#include +#include +#include #include "Math/vtkDense.h" namespace uLib { namespace Vtk { vtkEmitterPrimary::vtkEmitterPrimary(Geant::EmitterPrimary &emitter) - : m_emitter(emitter), m_Poly(nullptr), m_Actor(vtkActor::New()) { + : m_emitter(emitter), m_Actor(vtkActor::New()) { vtkNew arrow; - // Default arrow is along X+. Rotate to point towards Z- relative to the origin + // Default arrow is along X+. + // Scale it to 1 meter (1000 mm) and rotate to point towards Z- (local) vtkNew transform; + transform->Scale(1000.0, 1000.0, 1000.0); transform->RotateY(90.0); vtkNew transformFilter; @@ -28,7 +34,7 @@ vtkEmitterPrimary::vtkEmitterPrimary(Geant::EmitterPrimary &emitter) vtkSmartPointer mapper = vtkSmartPointer::New(); mapper->SetInputData(transformFilter->GetOutput()); m_Actor->SetMapper(mapper); - m_Actor->SetScale(1000.0); // 1 meter long + m_Actor->SetScale(1.0); vtkNew vmat; Matrix4fToVtk(m_emitter.GetWorldMatrix(), vmat); @@ -68,5 +74,53 @@ void vtkEmitterPrimary::Update() { m_emitter.Updated(); } +// -------------------------------------------------------------------------- // +// vtkSkyPlaneEmitterPrimary + +vtkSkyPlaneEmitterPrimary::vtkSkyPlaneEmitterPrimary(Geant::SkyPlaneEmitterPrimary &emitter) + : vtkEmitterPrimary(emitter), m_skyEmitter(emitter), m_PlaneSource(vtkPlaneSource::New()) { + + vtkNew append; + // Base class constructor already set an arrow. We keep it as a directional indicator. + append->AddInputData(vtkPolyData::SafeDownCast(m_Actor->GetMapper()->GetInput())); + append->AddInputConnection(m_PlaneSource->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(); +} + +vtkSkyPlaneEmitterPrimary::~vtkSkyPlaneEmitterPrimary() { + m_PlaneSource->Delete(); +} + +void vtkSkyPlaneEmitterPrimary::contentUpdate() { + uLib::Vector2f size = m_skyEmitter.GetSkySize(); + m_PlaneSource->SetOrigin(-size(0)/2.0, -size(1)/2.0, 0.0); + m_PlaneSource->SetPoint1( size(0)/2.0, -size(1)/2.0, 0.0); + m_PlaneSource->SetPoint2(-size(0)/2.0, size(1)/2.0, 0.0); + m_PlaneSource->Update(); + + vtkEmitterPrimary::contentUpdate(); +} + +// -------------------------------------------------------------------------- // +// vtkQuadMeshEmitterPrimary + +vtkQuadMeshEmitterPrimary::vtkQuadMeshEmitterPrimary(Geant::QuadMeshEmitterPrimary &emitter) + : vtkEmitterPrimary(emitter), m_meshEmitter(emitter) { + this->contentUpdate(); +} + +void vtkQuadMeshEmitterPrimary::contentUpdate() { + // For now stick with the arrow. In the future visualize the mesh if useful. + vtkEmitterPrimary::contentUpdate(); +} + } // namespace Vtk } // namespace uLib diff --git a/src/Vtk/HEP/Geant/vtkEmitterPrimary.h b/src/Vtk/HEP/Geant/vtkEmitterPrimary.h index 4780ed3..583df5f 100644 --- a/src/Vtk/HEP/Geant/vtkEmitterPrimary.h +++ b/src/Vtk/HEP/Geant/vtkEmitterPrimary.h @@ -8,6 +8,7 @@ class vtkConeSource; class vtkLineSource; class vtkPolyData; class vtkActor; +class vtkPlaneSource; namespace uLib { namespace Vtk { @@ -20,12 +21,36 @@ public: virtual void contentUpdate(); virtual void Update(); -private: +protected: Geant::EmitterPrimary &m_emitter; - vtkPolyData *m_Poly; vtkActor *m_Actor; }; +class vtkSkyPlaneEmitterPrimary : public vtkEmitterPrimary { +public: + vtkSkyPlaneEmitterPrimary(Geant::SkyPlaneEmitterPrimary &emitter); + virtual ~vtkSkyPlaneEmitterPrimary(); + + virtual void contentUpdate(); + +private: + Geant::SkyPlaneEmitterPrimary &m_skyEmitter; + vtkPlaneSource *m_PlaneSource; +}; + +class vtkQuadMeshEmitterPrimary : public vtkEmitterPrimary { +public: + vtkQuadMeshEmitterPrimary(Geant::QuadMeshEmitterPrimary &emitter); + virtual void contentUpdate(); + +private: + Geant::QuadMeshEmitterPrimary &m_meshEmitter; +}; + + + + + } // namespace Vtk } // namespace uLib