OpenCMT/uLib
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

geant events for multiple scattering

Browse Source
This commit is contained in:
AndreaRigoni
2026-03-14 23:33:31 +00:00
parent 692cdf7ae3
commit c63a1ae047
18 changed files with 681 additions and 90 deletions
Download Patch File Download Diff File Expand all files Collapse all files

32
src/HEP/Geant/ActionInitialization.cpp
View File

@@ -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<GeantEvent> *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

13
src/HEP/Geant/ActionInitialization.hh
View File

@@ -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<GeantEvent> *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<GeantEvent> *m_Output;
};
} // namespace Geant

2
src/HEP/Geant/CMakeLists.txt
View File

@@ -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)

57
src/HEP/Geant/GeantEvent.h
View File

@@ -30,34 +30,61 @@
#include "Core/Vector.h"
#include "Math/Dense.h"
#include "Detectors/ChamberHitEvent.h"
#include <string>
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<ChamberHitEventData>)
/// 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<Delta>)
private:
friend class GeantEvent;
Id_t m_EventID;
Scalarf m_Momentum;
Vector3f m_GenPos;
Vector3f m_GenDir;
Vector<ChamberHitEventData> m_ChEvents;
};
HLine3f m_GenVector;
Vector<Delta> m_Path;
class GeantEvent : public GeantEventData {
public:
uLibSetMacro(EventID, Id_t) uLibSetMacro(Momentum, Scalarf)
uLibRefMacro(GenPos, Vector3f) uLibRefMacro(GenDir, Vector3f)
uLibRefMacro(ChEvents, Vector<ChamberHitEventData>)
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

50
src/HEP/Geant/Scene.cpp
View File

@@ -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<Solid *> m_Solids;
Solid *m_World = nullptr;
G4RunManager *m_RunManager;
EmitterPrimary *m_Emitter;
Vector<GeantEvent> *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<GeantEvent> &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

12
src/HEP/Geant/Scene.h
View File

@@ -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<GeantEvent> &results);
private:
class SceneImpl *d;
};

18
src/HEP/Geant/Solid.cpp
View File

@@ -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() {

4
src/HEP/Geant/Solid.h
View File

@@ -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();

95
src/HEP/Geant/SteppingAction.cpp Normal file
View File

@@ -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<GeantEvent> *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<Id_t>(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<Scalarf>(prim->GetTotalMomentum() / MeV);
}
}
}
// Record a Delta for this step
GeantEvent::Delta delta;
// Step length
delta.m_Length = static_cast<Scalarf>(step->GetStepLength() / mm);
// Post-step momentum
G4ThreeVector postMom = track->GetMomentum();
delta.m_Momentum = static_cast<Scalarf>(postMom.mag() / MeV);
// Post-step direction
G4ThreeVector dir = track->GetMomentumDirection();
delta.m_Direction = HVector3f(static_cast<float>(dir.x()),
static_cast<float>(dir.y()),
static_cast<float>(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

30
src/HEP/Geant/SteppingAction.hh Normal file
View File

@@ -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<GeantEvent> *output);
virtual ~SteppingAction();
virtual void UserSteppingAction(const G4Step *step) override;
private:
Vector<GeantEvent> *m_Output; ///< destination for finished events
GeantEvent m_Current; ///< event being built
int m_LastEventID; ///< track event transitions
};
} // namespace Geant
} // namespace uLib
#endif

112
src/HEP/Geant/testing/EventTest.cpp
View File

@@ -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 <Geant4/G4Material.hh>
#include <Geant4/G4NistManager.hh>
@@ -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<Geant::GeantEvent> 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