diff --git a/include/ChargedParticle.h b/include/ChargedParticle.h index e34ea1c..42f5514 100644 --- a/include/ChargedParticle.h +++ b/include/ChargedParticle.h @@ -67,6 +67,12 @@ class ChargedParticle: public TransientParticle { void PIDReconstruction(CherenkovPID &pid); private: + struct TrajectoryData { + unsigned zdim; + const std::vector>& locations; + }; + + TrajectoryData GetTrajectoryData(CherenkovRadiator* radiator, RadiatorHistory* history) const; // Optical photons produced elsewhere; std::vector m_OrphanPhotons; diff --git a/include/CherenkovRadiator.h b/include/CherenkovRadiator.h index 9580678..e4725a8 100644 --- a/include/CherenkovRadiator.h +++ b/include/CherenkovRadiator.h @@ -21,8 +21,8 @@ class CherenkovRadiator: public TObject { CherenkovRadiator(const G4LogicalVolume *volume = 0, const G4RadiatorMaterial *material = 0): m_LogicalVolume(volume), m_Material(material), m_ReferenceRefractiveIndex(0.0), m_ReferenceAttenuationLength(0.0), - m_ID(0), m_Stat(0), m_AverageTheta(0.0), m_TrajectoryBinCount(1), m_Smearing(0.0), - m_GaussianSmearing(false) {}; + m_ID(0), m_Stat(0), m_AverageTheta(0.0), m_Smearing(0.0), + m_GaussianSmearing(false), m_TrajectoryBinCount_deprecated(1) {}; ~CherenkovRadiator() {}; double n( void ) const { return m_ReferenceRefractiveIndex; }; @@ -66,19 +66,21 @@ class CherenkovRadiator: public TObject { TString m_AlternativeMaterialName; public: - void SetTrajectoryBinCount(unsigned bins) { m_TrajectoryBinCount = bins; }; - double GetTrajectoryBinCount( void) const { return m_TrajectoryBinCount; }; - inline double GetSmearing( void ) const { return m_Smearing; }; void SetGaussianSmearing(double sigma) { m_GaussianSmearing = true; m_Smearing = sigma; } void SetUniformSmearing (double range) { m_GaussianSmearing = false; m_Smearing = range; } bool UseGaussianSmearing( void ) const { return m_GaussianSmearing; }; - // FIXME: memory leak; - void ResetLocations( void ) { m_Locations.clear(); } - void AddLocation(/*unsigned sector,*/ const TVector3 &x, const TVector3 &n) { - m_Locations/*[sector]*/.push_back(std::make_pair(x, n)); - }; + // Backwards compatibility: deprecated methods for old non-thread-safe API + // WARNING: These are NOT thread-safe and should only be used with mutex protection + // For thread-safe code, use RadiatorHistory methods directly + void SetTrajectoryBinCount(unsigned bins) { m_TrajectoryBinCount_deprecated = bins; } + unsigned GetTrajectoryBinCount() const { return m_TrajectoryBinCount_deprecated; } + void ResetLocations() { m_Locations_deprecated.clear(); } + void AddLocation(const TVector3 &x, const TVector3 &n) { + m_Locations_deprecated.push_back(std::make_pair(x, n)); + } + const std::vector>& GetLocations( void ) const { return m_Locations_deprecated; }; // Transient variables for the analysis script convenience; unsigned m_ID; //! @@ -87,17 +89,18 @@ class CherenkovRadiator: public TObject { TVector3 m_AverageVertexPosition; //! double m_AverageRefractiveIndex; //! - unsigned m_TrajectoryBinCount; //! // This is a hack for now; double m_Smearing; //! bool m_GaussianSmearing; //! - //std::map>> _m_Locations; //! - std::vector> m_Locations; //! std::vector> m_ri_lookup_table; //! + // Deprecated members for backwards compatibility (NOT thread-safe) + unsigned m_TrajectoryBinCount_deprecated; //! + std::vector> m_Locations_deprecated; //! + #ifndef DISABLE_ROOT_IO - ClassDef(CherenkovRadiator, 5); + ClassDef(CherenkovRadiator, 6); #endif }; diff --git a/include/RadiatorHistory.h b/include/RadiatorHistory.h index a32fcbd..20560ee 100644 --- a/include/RadiatorHistory.h +++ b/include/RadiatorHistory.h @@ -12,7 +12,7 @@ class RadiatorHistory: public TObject { public: - RadiatorHistory() {}; + RadiatorHistory(): m_TrajectoryBinCount(1) {}; ~RadiatorHistory() { for(auto photon: m_Photons) delete photon; @@ -31,7 +31,16 @@ class RadiatorHistory: public TObject { inline const ChargedParticleStep *GetStep(unsigned id) const { return (id < m_Steps.size() ? m_Steps[id] : 0); }; - inline unsigned StepCount( void ) { return m_Steps.size(); }; + inline unsigned StepCount( void ) { return m_Steps.size(); }; + + // Trajectory location management (moved from CherenkovRadiator for thread safety) + void SetTrajectoryBinCount(unsigned bins) { m_TrajectoryBinCount = bins; }; + unsigned GetTrajectoryBinCount( void ) const { return m_TrajectoryBinCount; }; + void ResetLocations( void ) { m_Locations.clear(); } + void AddLocation(const TVector3 &x, const TVector3 &n) { + m_Locations.push_back(std::make_pair(x, n)); + }; + const std::vector>& GetLocations( void ) const { return m_Locations; }; void AddStepBufferPoint(double time, const TVector3 &x) { // Will be in ascending order of time; @@ -72,8 +81,12 @@ class RadiatorHistory: public TObject { std::map m_StepBuffer; //! + // Per-event trajectory data (moved from CherenkovRadiator for thread safety) + unsigned m_TrajectoryBinCount; //! + std::vector> m_Locations; //! + #ifndef DISABLE_ROOT_IO - ClassDef(RadiatorHistory, 1); + ClassDef(RadiatorHistory, 2); #endif }; diff --git a/source/ChargedParticle.cc b/source/ChargedParticle.cc index fa194b9..75418d1 100644 --- a/source/ChargedParticle.cc +++ b/source/ChargedParticle.cc @@ -11,6 +11,28 @@ // ------------------------------------------------------------------------------------- +ChargedParticle::TrajectoryData ChargedParticle::GetTrajectoryData(CherenkovRadiator* radiator, RadiatorHistory* history) const { + bool history_has_data = !history->GetLocations().empty(); + bool radiator_has_data = !radiator->GetLocations().empty(); + + // Both APIs populated indicates accidental misuse + if (history_has_data && radiator_has_data) { + throw std::runtime_error("Both RadiatorHistory and CherenkovRadiator have trajectory data"); + } + + // Prefer new thread-safe API, fall back to deprecated API for backwards compatibility + if (history_has_data) { + return {history->GetTrajectoryBinCount(), history->GetLocations()}; + } else if (radiator_has_data) { + return {radiator->GetTrajectoryBinCount(), radiator->GetLocations()}; + } else { + static const std::vector> empty_locations; + return {0, empty_locations}; + } +} + +// ------------------------------------------------------------------------------------- + void ChargedParticle::PIDReconstruction(CherenkovPID &pid) { std::vector photons; @@ -42,8 +64,10 @@ void ChargedParticle::PIDReconstruction(CherenkovPID &pid) for(auto rhistory: GetRadiatorHistory()) { auto radiator = GetRadiator(rhistory); - unsigned zdim = radiator->GetTrajectoryBinCount(); - if (radiator->m_Locations.size() != zdim+1) continue; + auto history = GetHistory(rhistory); + + auto [zdim, locations] = GetTrajectoryData(radiator, history); + if (locations.size() != zdim+1) continue; TVector3 phx = photon->GetDetectionPosition(); @@ -63,9 +87,9 @@ void ChargedParticle::PIDReconstruction(CherenkovPID &pid) photon->m_Phi[radiator] = 0.0; for(unsigned iq=0; iqSolve(radiator->m_Locations[iq].first, + auto &solution = solutions[iq] = irt->Solve(locations[iq].first, // FIXME: give beam line as a parameter; - radiator->m_Locations[iq].second.Unit(), phx, TVector3(0,0,1), false); + locations[iq].second.Unit(), phx, TVector3(0,0,1), false); if (!solution.Converged()) { all_converged = false; break; @@ -74,11 +98,11 @@ void ChargedParticle::PIDReconstruction(CherenkovPID &pid) photon->m_Phi[radiator] += solution.GetPhi(); if (attenuation) { - TVector3 from = radiator->m_Locations[iq].first, to; + TVector3 from = locations[iq].first, to; bool ok = rear->GetCrossing(from, solution.m_Direction, &to); if (ok) { double length = (to - from).Mag(); - //printf("%02d -> %7.2f\n", iq, radiator->m_Locations[iq].first.z()); + //printf("%02d -> %7.2f\n", iq, locations[iq].first.z()); //weights[iq] = 1.0; weights[iq] = exp(-length / attenuation); } //if @@ -117,11 +141,13 @@ void ChargedParticle::PIDReconstruction(CherenkovPID &pid) for(auto rhistory: GetRadiatorHistory()) { auto radiator = GetRadiator(rhistory); - unsigned zdim = radiator->GetTrajectoryBinCount(); - if (radiator->m_Locations.size() != zdim+1) continue; + auto history = GetHistory(rhistory); + + auto [zdim, locations] = GetTrajectoryData(radiator, history); + if (locations.size() != zdim+1) continue; // Asume this estimate is good enough; - double pp = radiator->m_Locations[0].second.Mag(); + double pp = locations[0].second.Mag(); double arg = sqrt(pp*pp + m*m)/(radiator->m_AverageRefractiveIndex*pp); // Threshold check; FIXME: do it better?; if (fabs(arg) > 1.0) continue;