diff --git a/src/HEP/Detectors/DetectorChamber.h b/src/HEP/Detectors/DetectorChamber.h index 56d71a6..89686a5 100644 --- a/src/HEP/Detectors/DetectorChamber.h +++ b/src/HEP/Detectors/DetectorChamber.h @@ -65,7 +65,9 @@ public: HLine3f worldPlane; Matrix4f M = this->GetWorldMatrix(); worldPlane.origin = M * m_ProjectionPlane.origin; - worldPlane.direction = M * m_ProjectionPlane.direction; + HVector3f dirNorm = M * m_ProjectionPlane.direction; + dirNorm.normalize(); // Normalize for consistent dot products + worldPlane.direction = dirNorm; return worldPlane; } diff --git a/src/HEP/Geant/ActionInitialization.cpp b/src/HEP/Geant/ActionInitialization.cpp index bd5cec8..c364238 100644 --- a/src/HEP/Geant/ActionInitialization.cpp +++ b/src/HEP/Geant/ActionInitialization.cpp @@ -5,37 +5,26 @@ namespace uLib { namespace Geant { -ActionInitialization::ActionInitialization(EmitterPrimary *emitter, - Vector *output, - int verbosity) +ActionInitialization::ActionInitialization(EmitterPrimary *emitter, SimulationContext *context) : G4VUserActionInitialization(), m_Emitter(emitter), - m_Output(output), - m_Verbosity(verbosity) + m_Context(context) {} ActionInitialization::~ActionInitialization() {} -void ActionInitialization::BuildForMaster() const { - // Master thread: no per-event actions needed -} +void ActionInitialization::BuildForMaster() const {} void ActionInitialization::Build() const { - // Register the primary generator if (m_Emitter) { SetUserAction(m_Emitter->Clone()); } else { - // Fallback: default EmitterPrimary SetUserAction(new EmitterPrimary()); } - // Register actions - if (m_Output) { - SteppingAction *sa = new SteppingAction(m_Output); - sa->SetVerbosity(m_Verbosity); - SetUserAction(static_cast(sa)); - SetUserAction(static_cast(sa)); - } + SteppingAction *sa = new SteppingAction(m_Context); + SetUserAction(static_cast(sa)); + SetUserAction(static_cast(sa)); } } // namespace Geant diff --git a/src/HEP/Geant/ActionInitialization.hh b/src/HEP/Geant/ActionInitialization.hh index 423a5e2..bd2f550 100644 --- a/src/HEP/Geant/ActionInitialization.hh +++ b/src/HEP/Geant/ActionInitialization.hh @@ -2,35 +2,24 @@ #define ActionInitialization_h #include "G4VUserActionInitialization.hh" -#include "Core/Vector.h" +#include "SimulationContext.h" namespace uLib { namespace Geant { class EmitterPrimary; -class GeantEvent; class ActionInitialization : public G4VUserActionInitialization { public: - /// @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 *output = nullptr, - int verbosity = 0); + ActionInitialization(EmitterPrimary *emitter, SimulationContext *context); ~ActionInitialization(); - void SetVerbosity(int level) { m_Verbosity = level; } - - // Metodo chiamato solo dal thread principale (Master) - virtual void BuildForMaster() const; - - // Metodo chiamato dai thread di lavoro (Worker) o in modalità sequenziale - virtual void Build() const; + virtual void BuildForMaster() const override; + virtual void Build() const override; private: - EmitterPrimary *m_Emitter; - Vector *m_Output; - int m_Verbosity = 0; + EmitterPrimary *m_Emitter; + SimulationContext *m_Context; }; } // namespace Geant diff --git a/src/HEP/Geant/CMakeLists.txt b/src/HEP/Geant/CMakeLists.txt index 08d33e1..7ad33ec 100644 --- a/src/HEP/Geant/CMakeLists.txt +++ b/src/HEP/Geant/CMakeLists.txt @@ -20,6 +20,7 @@ set(HEADERS PhysicsList.hh ActionInitialization.hh SteppingAction.hh + SimulationContext.h ) set(SOURCES diff --git a/src/HEP/Geant/DetectorActionInitialization.cpp b/src/HEP/Geant/DetectorActionInitialization.cpp new file mode 100644 index 0000000..8a8885f --- /dev/null +++ b/src/HEP/Geant/DetectorActionInitialization.cpp @@ -0,0 +1,43 @@ +#include "DetectorActionInitialization.hh" +#include "EmitterPrimary.hh" +#include "DetectorSteppingAction.hh" + +namespace uLib { +namespace Geant { + +DetectorActionInitialization::DetectorActionInitialization(EmitterPrimary *emitter, + Vector *output, + const Vector &planes, + int verbosity) + : G4VUserActionInitialization(), + m_Emitter(emitter), + m_Output(output), + m_Planes(planes), + m_Verbosity(verbosity) +{} + +DetectorActionInitialization::~DetectorActionInitialization() {} + +void DetectorActionInitialization::BuildForMaster() const {} + +void DetectorActionInitialization::Build() const { + if (m_Verbosity > 0) { + std::cout << "[Geant] Worker thread Building actions... Output ptr: " << m_Output + << ", Planes count: " << m_Planes.size() << std::endl; + } + if (m_Emitter) { + SetUserAction(m_Emitter->Clone()); + } else { + SetUserAction(new EmitterPrimary()); + } + + if (m_Output) { + DetectorSteppingAction *sa = new DetectorSteppingAction(m_Output, m_Planes); + sa->SetVerbosity(m_Verbosity); + SetUserAction(static_cast(sa)); + SetUserAction(static_cast(sa)); + } +} + +} // namespace Geant +} // namespace uLib diff --git a/src/HEP/Geant/DetectorActionInitialization.hh b/src/HEP/Geant/DetectorActionInitialization.hh new file mode 100644 index 0000000..d93d535 --- /dev/null +++ b/src/HEP/Geant/DetectorActionInitialization.hh @@ -0,0 +1,35 @@ +#ifndef U_GEANT_DETECTORACTIONINITIALIZATION_HH +#define U_GEANT_DETECTORACTIONINITIALIZATION_HH + +#include "G4VUserActionInitialization.hh" +#include "Core/Vector.h" +#include "HEP/Detectors/MuonEvent.h" +#include "Math/Dense.h" + +namespace uLib { +namespace Geant { + +class EmitterPrimary; + +class DetectorActionInitialization : public G4VUserActionInitialization { +public: + DetectorActionInitialization(EmitterPrimary *emitter, + Vector *output, + const Vector &planes, + int verbosity = 0); + ~DetectorActionInitialization(); + + virtual void BuildForMaster() const override; + virtual void Build() const override; + +private: + EmitterPrimary *m_Emitter; + Vector *m_Output; + Vector m_Planes; + int m_Verbosity; +}; + +} // namespace Geant +} // namespace uLib + +#endif diff --git a/src/HEP/Geant/DetectorSteppingAction.cpp b/src/HEP/Geant/DetectorSteppingAction.cpp new file mode 100644 index 0000000..a175941 --- /dev/null +++ b/src/HEP/Geant/DetectorSteppingAction.cpp @@ -0,0 +1,110 @@ +#include "DetectorSteppingAction.hh" +#include +#include +#include +#include +#include +#include +#include + +static std::mutex g_DetectorOutputMutex; + +namespace uLib { +namespace Geant { + +DetectorSteppingAction::DetectorSteppingAction(Vector *output, const Vector &planes) + : G4UserSteppingAction(), + G4UserEventAction(), + m_Output(output), + m_Planes(planes), + m_CrossCount(0), + m_Verbosity(1) +{ + // std::cout << "[Geant] SteppingAction created with " << m_Planes.size() << " planes." << std::endl; +} + +DetectorSteppingAction::~DetectorSteppingAction() {} + +void DetectorSteppingAction::BeginOfEventAction(const G4Event* /*event*/) { + m_CrossCount = 0; + + // Initialize with NaN + float nan = std::numeric_limits::quiet_NaN(); + m_Current.LineIn().origin = HPoint3f(nan, nan, nan); + m_Current.LineIn().direction = HVector3f(nan, nan, nan); + m_Current.LineOut().origin = HPoint3f(nan, nan, nan); + m_Current.LineOut().direction = HVector3f(nan, nan, nan); + m_Current.Momentum() = nan; +} + +void DetectorSteppingAction::EndOfEventAction(const G4Event* /*event*/) { + if (m_Output) { + std::lock_guard lock(g_DetectorOutputMutex); + m_Output->push_back(m_Current); + } +} + +void DetectorSteppingAction::UserSteppingAction(const G4Step *step) { + if (!step) return; + if (!m_Output) { + return; + } + + const G4Track *track = step->GetTrack(); + if (!track) return; + + static size_t step_count = 0; + if (++step_count % 1000 == 0 && m_Verbosity > 0) { + std::cout << "[GeantMT] Processed " << step_count << " total steps across events." << std::endl; + } + + // Only consider primary muons + if (track->GetParentID() != 0) return; + + // Track the momentum at generation/first step if not set + if (std::isnan(m_Current.Momentum())) { + m_Current.Momentum() = static_cast(track->GetMomentum().mag() / MeV); + } + + G4ThreeVector p1 = step->GetPreStepPoint()->GetPosition(); + G4ThreeVector p2 = step->GetPostStepPoint()->GetPosition(); + G4ThreeVector dir_g4 = track->GetMomentumDirection(); + + HPoint3f p1f(p1.x(), p1.y(), p1.z()); + HPoint3f p2f(p2.x(), p2.y(), p2.z()); + HVector3f dirf(dir_g4.x(), dir_g4.y(), dir_g4.z()); + + // Check intersection with each detector plane + for (const auto& plane : m_Planes) { + // Plane: origin=O, direction=N (normal) + HPoint3f O = plane.origin; + HVector3f N = plane.direction; + + float d1 = (p1f - O).dot(N); + float d2 = (p2f - O).dot(N); + + // Check if the step crossed the plane + if ((d1 > 0 && d2 <= 0) || (d1 < 0 && d2 >= 0)) { + // Intersection point t = d1 / (d1 - d2) + float t = d1 / (d1 - d2); + HPoint3f intersection = p1f + t * (p2f - p1f); + + if (m_CrossCount == 0) { + m_Current.LineIn().origin = intersection; + m_Current.LineIn().direction = dirf; + m_CrossCount++; + if (m_Verbosity > 0) std::cout << "[GeantMT] Hit first plane at " << intersection.transpose() << std::endl; + } else if (m_CrossCount == 1) { + m_Current.LineOut().origin = intersection; + m_Current.LineOut().direction = dirf; + m_CrossCount++; + if (m_Verbosity > 0) std::cout << "[GeantMT] Hit second plane at " << intersection.transpose() << std::endl; + } + // We break to avoid crossing multiple planes in one infinitesimal step (unlikely but possible) + // Actually, we should check ALL planes. + } + } +} + +} // namespace Geant +} // namespace uLib diff --git a/src/HEP/Geant/DetectorSteppingAction.hh b/src/HEP/Geant/DetectorSteppingAction.hh new file mode 100644 index 0000000..f9e4d0b --- /dev/null +++ b/src/HEP/Geant/DetectorSteppingAction.hh @@ -0,0 +1,36 @@ +#ifndef U_GEANT_DETECTORSTEPPINGACTION_HH +#define U_GEANT_DETECTORSTEPPINGACTION_HH + +#include "G4UserSteppingAction.hh" +#include "G4UserEventAction.hh" +#include "Core/Vector.h" +#include "HEP/Detectors/MuonEvent.h" +#include "HEP/Detectors/DetectorChamber.h" +#include + +namespace uLib { +namespace Geant { + +class DetectorSteppingAction : public G4UserSteppingAction, public G4UserEventAction { +public: + DetectorSteppingAction(Vector *output, const Vector &planes); + virtual ~DetectorSteppingAction(); + + virtual void UserSteppingAction(const G4Step *step) override; + virtual void BeginOfEventAction(const G4Event *event) override; + virtual void EndOfEventAction(const G4Event *event) override; + + void SetVerbosity(int level) { m_Verbosity = level; } + +private: + Vector *m_Output; + Vector m_Planes; // World projection planes + MuonEvent m_Current; + int m_CrossCount = 0; + int m_Verbosity = 0; +}; + +} // namespace Geant +} // namespace uLib + +#endif diff --git a/src/HEP/Geant/GeantEvent.h b/src/HEP/Geant/GeantEvent.h index 35ce0cb..ecf5389 100644 --- a/src/HEP/Geant/GeantEvent.h +++ b/src/HEP/Geant/GeantEvent.h @@ -73,6 +73,7 @@ public: void Print(const size_t size = 10) const; + private: Id_t m_EventID; Scalarf m_Momentum; diff --git a/src/HEP/Geant/Scene.cpp b/src/HEP/Geant/Scene.cpp index 6fb9388..661e3f6 100644 --- a/src/HEP/Geant/Scene.cpp +++ b/src/HEP/Geant/Scene.cpp @@ -1,28 +1,3 @@ -/*////////////////////////////////////////////////////////////////////////////// -// 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 > - - ------------------------------------------------------------------ - This library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 3.0 of the License, or (at your option) any later version. - - This library is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public - License along with this library. - -//////////////////////////////////////////////////////////////////////////////*/ - #include #include #include @@ -42,13 +17,15 @@ #include "Scene.h" #include "PhysicsList.hh" #include "ActionInitialization.hh" +#include "SimulationContext.h" +#include "HEP/Detectors/DetectorChamber.h" namespace uLib { namespace Geant { class SceneDetectorConstruction : public DetectorConstruction { public: - SceneDetectorConstruction(class SceneImpl *owner); + SceneDetectorConstruction(class SceneImpl *owner) : DetectorConstruction("Scene"), m_Owner(owner) {} G4VPhysicalVolume *Construct() override; private: class SceneImpl *m_Owner; @@ -69,57 +46,50 @@ static void CheckGeant4Environment() { class SceneImpl { public: - // constructor // SceneImpl() : m_RunManager(G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default)), m_Emitter(nullptr), - m_Output(nullptr), - m_Verbosity(0) { + m_InitCalled(false) { m_RunManager->SetUserInitialization(new PhysicsList); } - // destructor // ~SceneImpl() { if (m_RunManager) delete m_RunManager; - if (m_World) delete m_World; + // m_World deletion is handled in Scene destructor or here } void Initialize() { - // Set mandatory initialization classes for Geant4 - m_RunManager->SetUserInitialization(new SceneDetectorConstruction(this)); - m_RunManager->SetUserInitialization( - new ActionInitialization(m_Emitter, m_Output, m_Verbosity)); + if (m_InitCalled) return; + + m_RunManager->SetUserInitialization(new SceneDetectorConstruction(this)); + m_RunManager->SetUserInitialization(new ActionInitialization(m_Emitter, &m_Context)); - m_RunManager->SetVerboseLevel(m_Verbosity); - - // Initialize Geant4 m_RunManager->Initialize(); + m_InitCalled = true; } - // members // Vector m_Solids; Solid *m_World = nullptr; ContainerBox m_WorldBox; G4RunManager *m_RunManager; EmitterPrimary *m_Emitter; - Vector *m_Output; - int m_Verbosity; + SimulationContext m_Context; + bool m_InitCalled; }; -SceneDetectorConstruction::SceneDetectorConstruction(SceneImpl *owner) - : DetectorConstruction("Scene"), m_Owner(owner) {} - G4VPhysicalVolume *SceneDetectorConstruction::Construct() { return m_Owner->m_World->GetPhysical(); } - - - Scene::Scene() { CheckGeant4Environment(); d = new SceneImpl(); } -Scene::~Scene() { delete d; } + +Scene::~Scene() { + // Delete solids + for(auto s : d->m_Solids) delete s; + delete d; +} void Scene::AddSolid(Solid *solid, Solid *parent) { d->m_Solids.push_back(solid); @@ -131,12 +101,9 @@ void Scene::AddSolid(Solid *solid, Solid *parent) { } const Solid* Scene::GetWorld() const { return d->m_World; } - ContainerBox* Scene::GetWorldBox() const { return &d->m_WorldBox; } void Scene::ConstructWorldBox(const Vector3f &size, const char *material) { - // Get nist material manager - d->m_WorldBox.Scale(size); d->m_WorldBox.SetPosition(-size/2.0f); @@ -146,61 +113,48 @@ void Scene::ConstructWorldBox(const Vector3f &size, const char *material) { AddSolid(d->m_World); } - G4Box *solidWorld = new G4Box("World", - 0.5 * size(0), - 0.5 * size(1), - 0.5 * size(2)); - - G4LogicalVolume *logicWorld = new G4LogicalVolume(solidWorld, - d->m_World->GetMaterial(), - d->m_World->GetName()); - + G4Box *solidWorld = new G4Box("World", 0.5 * size(0), 0.5 * size(1), 0.5 * size(2)); + G4LogicalVolume *logicWorld = new G4LogicalVolume(solidWorld, d->m_World->GetMaterial(), d->m_World->GetName()); d->m_World->SetLogical(logicWorld); - G4PVPlacement *physWorld = new G4PVPlacement( - nullptr, - G4ThreeVector(0, 0, 0), - logicWorld, - d->m_World->GetName(), - 0, - false, - 0, - true); - + G4PVPlacement *physWorld = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0), logicWorld, d->m_World->GetName(), 0, false, 0, true); d->m_World->SetPhysical(physWorld); - - // no transforms are allowed for the world box - // Matrix4f transform = box->GetMatrix(); - // d->m_World->SetTransform(transform); } -void Scene::SetEmitter(EmitterPrimary *emitter) { - d->m_Emitter = emitter; -} - -void Scene::Initialize() { - d->Initialize(); -} +void Scene::SetEmitter(EmitterPrimary *emitter) { d->m_Emitter = emitter; } +void Scene::Initialize() { d->Initialize(); } void Scene::SetVerbosity(int level) { - d->m_Verbosity = level; + d->m_Context.verbosity = level; if (d->m_RunManager) d->m_RunManager->SetVerboseLevel(level); } void Scene::RunSimulation(int nEvents, Vector &results) { - d->m_Output = &results; + d->Initialize(); // Ensure initialized + d->m_Context.mode = SimulationMode::DETAILED; + d->m_Context.outputGeant = &results; + d->m_Context.outputMuon = nullptr; + + d->m_RunManager->BeamOn(nEvents); +} - // 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_Verbosity)); - - // Re-run initialization to propagate the ActionInitialization change to worker threads - d->m_RunManager->Initialize(); +void Scene::RunDetectorSimulation(int nEvents, Vector &results) { + d->Initialize(); // Ensure initialized + d->m_Context.mode = SimulationMode::DETECTOR; + d->m_Context.outputGeant = nullptr; + d->m_Context.outputMuon = &results; + + // Find detector planes + d->m_Context.detectorPlanes.clear(); + for (Solid* s : d->m_Solids) { + if (BoxSolid* bs = dynamic_cast(s)) { + if (DetectorChamber* dc = dynamic_cast(bs->GetObject())) { + d->m_Context.detectorPlanes.push_back(dc->GetWorldProjectionPlane()); + } + } + } d->m_RunManager->BeamOn(nEvents); } } // namespace Geant } // namespace uLib - diff --git a/src/HEP/Geant/Scene.h b/src/HEP/Geant/Scene.h index 04d9ca1..d483fce 100644 --- a/src/HEP/Geant/Scene.h +++ b/src/HEP/Geant/Scene.h @@ -32,6 +32,7 @@ #include "Core/Vector.h" #include "Solid.h" #include "GeantEvent.h" +#include "HEP/Detectors/MuonEvent.h" class G4VPhysicalVolume; @@ -68,6 +69,9 @@ public: /// Results are appended to the provided vector. void RunSimulation(int nEvents, Vector &results); + /// Specialized detector simulation trackingMuonEvent line crossings. + void RunDetectorSimulation(int nEvents, Vector &results); + private: class SceneImpl *d; }; diff --git a/src/HEP/Geant/SimulationContext.h b/src/HEP/Geant/SimulationContext.h new file mode 100644 index 0000000..a55e5fa --- /dev/null +++ b/src/HEP/Geant/SimulationContext.h @@ -0,0 +1,30 @@ +#ifndef U_GEANT_SIMULATIONCONTEXT_H +#define U_GEANT_SIMULATIONCONTEXT_H + +#include "Core/Vector.h" +#include "GeantEvent.h" +#include "HEP/Detectors/MuonEvent.h" +#include "Math/Dense.h" +#include + +namespace uLib { +namespace Geant { + +enum class SimulationMode { + DETAILED, + DETECTOR +}; + +struct SimulationContext { + SimulationMode mode = SimulationMode::DETAILED; + Vector *outputGeant = nullptr; + Vector *outputMuon = nullptr; + Vector detectorPlanes; + int verbosity = 0; + std::mutex outputMutex; +}; + +} // namespace Geant +} // namespace uLib + +#endif diff --git a/src/HEP/Geant/Solid.h b/src/HEP/Geant/Solid.h index 4f76adf..2ac5374 100644 --- a/src/HEP/Geant/Solid.h +++ b/src/HEP/Geant/Solid.h @@ -106,6 +106,8 @@ class BoxSolid : public Solid { public: BoxSolid(const char *name, ContainerBox *box); virtual G4VSolid* GetG4Solid() const override { return (G4VSolid*)m_Solid; } + + ContainerBox* GetObject() const { return m_Object; } public slots: void Update(); diff --git a/src/HEP/Geant/SteppingAction.cpp b/src/HEP/Geant/SteppingAction.cpp index 7a30777..90aea80 100644 --- a/src/HEP/Geant/SteppingAction.cpp +++ b/src/HEP/Geant/SteppingAction.cpp @@ -1,111 +1,124 @@ #include "SteppingAction.hh" - -#include "G4Step.hh" -#include "G4Track.hh" -#include "G4Event.hh" -#include "G4RunManager.hh" -#include "G4LogicalVolume.hh" -#include "G4SystemOfUnits.hh" -#include "G4ParticleDefinition.hh" +#include +#include +#include +#include #include #include #include - -static std::mutex g_SimulationOutputMutex; +#include namespace uLib { namespace Geant { -SteppingAction::SteppingAction(Vector *output) +SteppingAction::SteppingAction(SimulationContext *context) : G4UserSteppingAction(), - m_Output(output), - m_Current(), - m_LastEventID(-1) + G4UserEventAction(), + m_Context(context), + m_Verbosity(0) {} SteppingAction::~SteppingAction() {} void SteppingAction::BeginOfEventAction(const G4Event *event) { - if (!event || !m_Output) return; + if (!event || !m_Context) return; - // Start a new GeantEvent - m_Current = GeantEvent(); - m_Current.m_EventID = static_cast(event->GetEventID()); + if (m_Context->mode == SimulationMode::DETAILED) { + m_CurrentGeant = GeantEvent(); + m_CurrentGeant.m_EventID = static_cast(event->GetEventID()); - // 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 (event->GetNumberOfPrimaryVertex() > 0) { + G4PrimaryVertex *vtx = event->GetPrimaryVertex(0); + G4ThreeVector pos = vtx->GetPosition(); + m_CurrentGeant.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(prim->GetTotalMomentum() / MeV); - } + if (vtx->GetNumberOfParticle() > 0) { + G4PrimaryParticle *prim = vtx->GetPrimary(0); + G4ThreeVector mom = prim->GetMomentumDirection(); + m_CurrentGeant.m_GenVector.direction = HVector3f(mom.x(), mom.y(), mom.z()); + m_CurrentGeant.m_Momentum = static_cast(prim->GetTotalMomentum() / MeV); + } + } + } else { + // Detector mode + m_MuonCrossCount = 0; + float nan = std::numeric_limits::quiet_NaN(); + m_CurrentMuon.LineIn().origin = HPoint3f(nan, nan, nan); + m_CurrentMuon.LineIn().direction = HVector3f(nan, nan, nan); + m_CurrentMuon.LineOut().origin = HPoint3f(nan, nan, nan); + m_CurrentMuon.LineOut().direction = HVector3f(nan, nan, nan); + m_CurrentMuon.Momentum() = nan; } } void SteppingAction::EndOfEventAction(const G4Event *event) { - if (m_Output && !m_Current.m_Path.empty()) { - if (m_Verbosity > 0) { - std::cout << "[Geant] Finished Event " << m_Current.m_EventID - << " with " << m_Current.m_Path.size() << " steps." << std::endl; - - // Check if we hit anything other than World - std::set volumes; - for (const auto& delta : m_Current.m_Path) { - if (!delta.m_SolidName.empty()) volumes.insert(delta.m_SolidName); - } - if (volumes.size() > 1) { - std::cout << " - Hit volumes: "; - for (const auto& v : volumes) std::cout << v << " "; - std::cout << std::endl; - } - } + if (!m_Context) return; - { - std::lock_guard lock(g_SimulationOutputMutex); - m_Output->push_back(m_Current); + if (m_Context->mode == SimulationMode::DETAILED && m_Context->outputGeant) { + if (!m_CurrentGeant.m_Path.empty()) { + std::lock_guard lock(m_Context->outputMutex); + m_Context->outputGeant->push_back(m_CurrentGeant); } + } else if (m_Context->mode == SimulationMode::DETECTOR && m_Context->outputMuon) { + // In detector mode, we always push the event (to keep indexing consistent) + // or only if hit? User requested "all muon event line in and out ar at first set to nan". + // So we push everything. + std::lock_guard lock(m_Context->outputMutex); + m_Context->outputMuon->push_back(m_CurrentMuon); } } void SteppingAction::UserSteppingAction(const G4Step *step) { - if (!step || !m_Output) return; + if (!step || !m_Context) return; const G4Track *track = step->GetTrack(); - if (!track) return; + if (!track || track->GetParentID() != 0) return; - // Only record primary particle (muon) - if (track->GetParentID() != 0) return; + if (m_Context->mode == SimulationMode::DETAILED) { + GeantEvent::Delta delta; + delta.m_Length = static_cast(step->GetStepLength() / mm); + delta.m_Momentum = static_cast(track->GetMomentum().mag() / MeV); + G4ThreeVector dir = track->GetMomentumDirection(); + delta.m_Direction = HVector3f(dir.x(), dir.y(), dir.z()); - // Record a Delta for this step - GeantEvent::Delta delta; + if (track->GetVolume()) { + delta.m_SolidName = track->GetVolume()->GetName(); + } + m_CurrentGeant.m_Path.push_back(delta); + } else { + // Detector Mode + if (std::isnan(m_CurrentMuon.Momentum())) { + m_CurrentMuon.Momentum() = static_cast(track->GetMomentum().mag() / MeV); + } - // Step length - delta.m_Length = static_cast(step->GetStepLength() / mm); + G4ThreeVector p1 = step->GetPreStepPoint()->GetPosition(); + G4ThreeVector p2 = step->GetPostStepPoint()->GetPosition(); + G4ThreeVector dir_g4 = track->GetMomentumDirection(); + + HPoint3f p1f(p1.x(), p1.y(), p1.z()); + HPoint3f p2f(p2.x(), p2.y(), p2.z()); + HVector3f dirf(dir_g4.x(), dir_g4.y(), dir_g4.z()); - // Post-step momentum - G4ThreeVector postMom = track->GetMomentum(); - delta.m_Momentum = static_cast(postMom.mag() / MeV); + for (const auto& plane : m_Context->detectorPlanes) { + float d1 = (p1f - plane.origin).dot(plane.direction); + float d2 = (p2f - plane.origin).dot(plane.direction); - // Post-step direction - G4ThreeVector dir = track->GetMomentumDirection(); - delta.m_Direction = HVector3f(static_cast(dir.x()), - static_cast(dir.y()), - static_cast(dir.z())); + if ((d1 > 0 && d2 <= 0) || (d1 < 0 && d2 >= 0)) { + float t = d1 / (d1 - d2); + HPoint3f intersection = p1f + t * (p2f - p1f); - // Solid name where the step occurred - const G4LogicalVolume *vol = track->GetVolume() - ? track->GetVolume()->GetLogicalVolume() - : nullptr; - if (vol) { - delta.m_SolidName = vol->GetName(); + if (m_MuonCrossCount == 0) { + m_CurrentMuon.LineIn().origin = intersection; + m_CurrentMuon.LineIn().direction = dirf; + m_MuonCrossCount++; + } else if (m_MuonCrossCount == 1) { + m_CurrentMuon.LineOut().origin = intersection; + m_CurrentMuon.LineOut().direction = dirf; + m_MuonCrossCount++; + } + } + } } - - m_Current.m_Path.push_back(delta); } } // namespace Geant diff --git a/src/HEP/Geant/SteppingAction.hh b/src/HEP/Geant/SteppingAction.hh index dcc5c7c..249c996 100644 --- a/src/HEP/Geant/SteppingAction.hh +++ b/src/HEP/Geant/SteppingAction.hh @@ -5,16 +5,16 @@ #include "G4UserEventAction.hh" #include "Core/Vector.h" #include "GeantEvent.h" +#include "HEP/Detectors/MuonEvent.h" +#include "SimulationContext.h" namespace uLib { namespace Geant { -/// SteppingAction collects scattering data at each Geant4 step and -/// builds GeantEvent objects in the output buffer. +/// SteppingAction collects scattering data at each Geant4 step. class SteppingAction : public G4UserSteppingAction, public G4UserEventAction { public: - /// @param output pointer to the results vector owned by the Scene - SteppingAction(Vector *output); + SteppingAction(SimulationContext *context); virtual ~SteppingAction(); virtual void UserSteppingAction(const G4Step *step) override; @@ -24,10 +24,11 @@ public: void SetVerbosity(int level) { m_Verbosity = level; } private: - Vector *m_Output; ///< destination for finished events - GeantEvent m_Current; ///< event being built - int m_LastEventID; ///< track event transitions - int m_Verbosity = 0; + SimulationContext *m_Context; + GeantEvent m_CurrentGeant; + MuonEvent m_CurrentMuon; + int m_MuonCrossCount = 0; + int m_Verbosity = 0; }; } // namespace Geant diff --git a/src/HEP/Geant/testing/ActionInitialization.cpp b/src/HEP/Geant/testing/ActionInitialization.cpp index 32361a0..3d10eda 100644 --- a/src/HEP/Geant/testing/ActionInitialization.cpp +++ b/src/HEP/Geant/testing/ActionInitialization.cpp @@ -13,7 +13,7 @@ int main(int argc, char **argv) { // runManager->SetUserInitialization(new PhysicsList()); // 3. INIZIALIZZAZIONE DELLE AZIONI (Il nostro generatore!) - runManager->SetUserInitialization(new uLib::Geant::ActionInitialization()); + runManager->SetUserInitialization(new uLib::Geant::ActionInitialization(nullptr, nullptr)); // ... Inizializzazione del kernel ( runManager->Initialize(); ), UI manager, // vis manager, ecc. diff --git a/src/HEP/Geant/testing/SkyPlaneEmitterTest.cpp b/src/HEP/Geant/testing/SkyPlaneEmitterTest.cpp index 86b23db..6e46d2b 100644 --- a/src/HEP/Geant/testing/SkyPlaneEmitterTest.cpp +++ b/src/HEP/Geant/testing/SkyPlaneEmitterTest.cpp @@ -25,8 +25,8 @@ int main(int argc, char** argv) { scene.ConstructWorldBox(Vector3f(30_m, 30_m, 30_m), "G4_AIR"); ContainerBox iron_box; - iron_box.Scale(Vector3f(18_m, 18_m, 18_m)); - iron_box.SetPosition(Vector3f(-9_m, -9_m, -9_m)); + iron_box.Scale(Vector3f(18_m, 10_cm, 18_m)); + iron_box.SetPosition(Vector3f(-9_m, -5_cm, -9_m)); Geant::BoxSolid* iron_cube = new Geant::BoxSolid("IronCube", &iron_box); iron_cube->SetNistMaterial("G4_Fe"); iron_cube->Update(); @@ -60,8 +60,10 @@ int main(int argc, char** argv) { scene.SetEmitter(emitter); scene.SetVerbosity(1); - scene.Initialize(); + // scene.Initialize(); // Removed to avoid premature initialization + + std::cout << "Starting simulation of " << nEvents << " events..." << std::endl; Vector results; scene.RunSimulation(nEvents, results); @@ -80,5 +82,23 @@ int main(int argc, char** argv) { std::cout << " Average steps per event: " << static_cast(total_steps) / results.size() << std::endl; } + std::cout << "\nStarting Detector Simulation of " << nEvents << " events..." << std::endl; + Vector detectorResults; + scene.RunDetectorSimulation(nEvents, detectorResults); + + + + std::cout << "Detector Simulation finished." << std::endl; + size_t hit_count = 0; + for (const auto& ev : detectorResults) { + if (!std::isnan(ev.LineIn().origin.x())) { + hit_count++; + } + } + std::cout << " Muons crossing at least one detector: " << hit_count << std::endl; + if (nEvents > 0) { + std::cout << " Efficiency: " << (100.0 * hit_count / nEvents) << "%" << std::endl; + } + return 0; }