added skyplaneEmitter
This commit is contained in:
@@ -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<double, 3> 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) {
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -2,6 +2,7 @@
|
||||
set(TESTS
|
||||
vtkGeantEventTest
|
||||
vtkEmitterPrimaryTest
|
||||
vtkSkyPlaneEmitterPrimaryTest
|
||||
)
|
||||
|
||||
set(LIBRARIES
|
||||
|
||||
125
src/Vtk/HEP/Geant/testing/vtkSkyPlaneEmitterPrimaryTest.cpp
Normal file
125
src/Vtk/HEP/Geant/testing/vtkSkyPlaneEmitterPrimaryTest.cpp
Normal file
@@ -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 <vtkSmartPointer.h>
|
||||
#include <vtkCallbackCommand.h>
|
||||
#include <vtkRenderWindowInteractor.h>
|
||||
#include <vtkProperty.h>
|
||||
#include <vtkRenderer.h>
|
||||
#include <vtkRenderWindow.h>
|
||||
|
||||
#include <Geant4/Randomize.hh>
|
||||
#include <Geant4/G4ParticleGun.hh>
|
||||
#include <Geant4/G4SystemOfUnits.hh>
|
||||
|
||||
#include <iostream>
|
||||
|
||||
using namespace uLib;
|
||||
|
||||
struct AppState {
|
||||
Geant::Scene* scene;
|
||||
Vtk::Viewer* viewer;
|
||||
Vector<Geant::GeantEvent> results;
|
||||
};
|
||||
|
||||
void KeyPressCallbackFunction(vtkObject* caller, long unsigned int eventId, void* clientData, void* callData) {
|
||||
auto* interactor = static_cast<vtkRenderWindowInteractor*>(caller);
|
||||
auto* state = static_cast<AppState*>(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<vtkCallbackCommand> keyCallback = vtkSmartPointer<vtkCallbackCommand>::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;
|
||||
}
|
||||
@@ -6,18 +6,24 @@
|
||||
#include <vtkPolyDataMapper.h>
|
||||
#include <vtkTransform.h>
|
||||
#include <vtkTransformPolyDataFilter.h>
|
||||
#include <vtkPlaneSource.h>
|
||||
#include <vtkAppendPolyData.h>
|
||||
#include <vtkPolyData.h>
|
||||
#include <vtkProperty.h>
|
||||
#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<vtkArrowSource> 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<vtkTransform> transform;
|
||||
transform->Scale(1000.0, 1000.0, 1000.0);
|
||||
transform->RotateY(90.0);
|
||||
|
||||
vtkNew<vtkTransformPolyDataFilter> transformFilter;
|
||||
@@ -28,7 +34,7 @@ vtkEmitterPrimary::vtkEmitterPrimary(Geant::EmitterPrimary &emitter)
|
||||
vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
|
||||
mapper->SetInputData(transformFilter->GetOutput());
|
||||
m_Actor->SetMapper(mapper);
|
||||
m_Actor->SetScale(1000.0); // 1 meter long
|
||||
m_Actor->SetScale(1.0);
|
||||
|
||||
vtkNew<vtkMatrix4x4> 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<vtkAppendPolyData> 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<vtkPolyDataMapper> 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
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
Reference in New Issue
Block a user