add cylinder emitter

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<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);
+
+    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);
+}
 
 // -------------------------------------------------------------------------- //
 

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

src/Vtk/HEP/Geant/testing/CMakeLists.txt

@@ -3,6 +3,7 @@ set(TESTS
        vtkGeantEventTest
        vtkEmitterPrimaryTest
        vtkSkyPlaneEmitterPrimaryTest
+       vtkCylinderEmitterPrimaryTest
 )
 
 set(LIBRARIES

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 <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 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<vtkCallbackCommand> keyCallback = vtkSmartPointer<vtkCallbackCommand>::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;
+}

src/Vtk/HEP/Geant/vtkEmitterPrimary.cpp

@@ -10,6 +10,7 @@
 #include <vtkAppendPolyData.h>
 #include <vtkPolyData.h>
 #include <vtkProperty.h>
+#include <vtkCylinderSource.h>
 #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<vtkTransform> tcyl;
+    tcyl->RotateX(90.0);
+    
+    vtkNew<vtkTransformPolyDataFilter> cylFilter;
+    cylFilter->SetTransform(tcyl);
+    cylFilter->SetInputConnection(m_CylinderSource->GetOutputPort());
+
+    vtkNew<vtkAppendPolyData> append;
+    // Keep the directional arrow from base class
+    append->AddInputData(vtkPolyData::SafeDownCast(m_Actor->GetMapper()->GetInput()));
+    append->AddInputConnection(cylFilter->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();
+}
+
+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

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;
+};
+