From 1f9dc893cf7eaa39c35fa1f9583a13b755019e94 Mon Sep 17 00:00:00 2001 From: majacquet Date: Fri, 29 Nov 2024 14:52:49 +0100 Subject: [PATCH 1/3] a first version of the ForceCollision actor --- .../g4_bindings/pyPhysicsLists.cpp | 8 + core/opengate_core/opengate_core.cpp | 4 +- .../opengate_lib/GateOptnCloning.cpp | 51 ++ .../opengate_lib/GateOptnCloning.h | 87 ++++ .../GateOptnForceCommonTruncatedExp.cpp | 182 ++++++++ .../GateOptnForceCommonTruncatedExp.h | 122 +++++ .../opengate_lib/GateOptnForceFreeFlight.cpp | 104 +++++ .../opengate_lib/GateOptnForceFreeFlight.h | 97 ++++ .../opengate_lib/GateOptrForceCollision.cpp | 442 ++++++++++++++++++ .../opengate_lib/GateOptrForceCollision.h | 105 +++++ .../GateOptrForceCollisionTrackData.cpp | 74 +++ .../GateOptrForceCollisionTrackData.h | 80 ++++ .../opengate_lib/pyGateOptrForceCollision.cpp | 19 + opengate/actors/miscactors.py | 75 +-- opengate/engines.py | 14 +- opengate/managers.py | 43 +- 16 files changed, 1452 insertions(+), 55 deletions(-) create mode 100644 core/opengate_core/opengate_lib/GateOptnCloning.cpp create mode 100644 core/opengate_core/opengate_lib/GateOptnCloning.h create mode 100644 core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.cpp create mode 100644 core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.h create mode 100644 core/opengate_core/opengate_lib/GateOptnForceFreeFlight.cpp create mode 100644 core/opengate_core/opengate_lib/GateOptnForceFreeFlight.h create mode 100644 core/opengate_core/opengate_lib/GateOptrForceCollision.cpp create mode 100644 core/opengate_core/opengate_lib/GateOptrForceCollision.h create mode 100644 core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.cpp create mode 100644 core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.h create mode 100644 core/opengate_core/opengate_lib/pyGateOptrForceCollision.cpp diff --git a/core/opengate_core/g4_bindings/pyPhysicsLists.cpp b/core/opengate_core/g4_bindings/pyPhysicsLists.cpp index 354807e88..e365be7c9 100644 --- a/core/opengate_core/g4_bindings/pyPhysicsLists.cpp +++ b/core/opengate_core/g4_bindings/pyPhysicsLists.cpp @@ -90,6 +90,14 @@ namespace py = pybind11; .def("PhysicsBias", \ py::overload_cast &>( \ &G4GenericBiasingPhysics::PhysicsBias), \ + py::return_value_policy::reference_internal) \ + .def("Bias", \ + py::overload_cast &>( \ + &G4GenericBiasingPhysics::Bias), \ + py::return_value_policy::reference_internal) \ + .def("NonPhysicsBias", \ + py::overload_cast( \ + &G4GenericBiasingPhysics::NonPhysicsBias), \ py::return_value_policy::reference_internal); namespace pyPhysicsLists { diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index d34ba5f80..9020a1859 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -345,7 +345,7 @@ void init_GateSimulationStatisticsActor(py::module &); void init_GatePhaseSpaceActor(py::module &); -// void init_GateComptonSplittingActor(py::module &); +void init_GateOptrForceCollision(py::module &m); void init_GateOptrComptSplittingActor(py::module &m); @@ -570,7 +570,7 @@ PYBIND11_MODULE(opengate_core, m) { init_GateLETActor(m); init_GateSimulationStatisticsActor(m); init_GatePhaseSpaceActor(m); - // init_GateComptonSplittingActor(m); + init_GateOptrForceCollision(m); init_GateBOptrBremSplittingActor(m); init_GateOptrComptSplittingActor(m); init_GateHitsCollectionActor(m); diff --git a/core/opengate_core/opengate_lib/GateOptnCloning.cpp b/core/opengate_core/opengate_lib/GateOptnCloning.cpp new file mode 100644 index 000000000..91e5e5803 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnCloning.cpp @@ -0,0 +1,51 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +#include "GateOptnCloning.h" + + +GateOptnCloning::GateOptnCloning(G4String name) + : G4VBiasingOperation( name ), + fClone1W ( -1.0 ), + fClone2W ( -1.0 ), + fParticleChange(), + fCloneTrack ( nullptr ) +{} + +GateOptnCloning::~GateOptnCloning() +{} + +G4VParticleChange* GateOptnCloning::GenerateBiasingFinalState( const G4Track* track, + const G4Step* ) +{ + fParticleChange.Initialize(*track); + fParticleChange.ProposeParentWeight( fClone1W ); + fParticleChange.SetSecondaryWeightByProcess(true); + fParticleChange.SetNumberOfSecondaries(1); + fCloneTrack = new G4Track( *track ); + fCloneTrack->SetWeight( fClone2W ); + fParticleChange.AddSecondary( fCloneTrack ); + return &fParticleChange; +} diff --git a/core/opengate_core/opengate_lib/GateOptnCloning.h b/core/opengate_core/opengate_lib/GateOptnCloning.h new file mode 100644 index 000000000..fd6033fbd --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnCloning.h @@ -0,0 +1,87 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// +//--------------------------------------------------------------- +// +// GateOptnCloning +// +// Class Description: +// A G4VBiasingOperation to clone a track, allowing to set +// weight arbitrary weights. +// +// +//--------------------------------------------------------------- +// Initial version Nov. 2013 M. Verderi + + +#ifndef GateOptnCloning_h +#define GateOptnCloning_h 1 + +#include "G4VBiasingOperation.hh" +#include "G4ParticleChange.hh" + +class GateOptnCloning : public G4VBiasingOperation { +public: + // -- Constructor : + GateOptnCloning(G4String name); + // -- destructor: + virtual ~GateOptnCloning(); + +public: + // -- Methods from G4VBiasingOperation interface: + // ------------------------------------------- + // -- Unsed: + virtual const G4VBiasingInteractionLaw* ProvideOccurenceBiasingInteractionLaw( const G4BiasingProcessInterface*, G4ForceCondition& ) {return 0;} + virtual G4VParticleChange* ApplyFinalStateBiasing( const G4BiasingProcessInterface*, + const G4Track*, + const G4Step*, + G4bool& ) {return 0;} + // -- Used: + virtual G4double DistanceToApplyOperation( const G4Track*, + G4double, + G4ForceCondition* condition) + { + *condition = NotForced; return 0; // -- acts immediately + } + virtual G4VParticleChange* GenerateBiasingFinalState( const G4Track*, + const G4Step* ); + +public: + // -- Additional methods, specific to this class: + // ---------------------------------------------- + void SetCloneWeights(G4double clone1Weight, G4double clone2Weight) {fClone1W = clone1Weight ; fClone2W = clone2Weight;} + + G4Track* GetCloneTrack() const { return fCloneTrack; } + +private: + G4double fClone1W, + fClone2W; + G4ParticleChange fParticleChange; + G4Track* fCloneTrack; +}; + +#endif diff --git a/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.cpp b/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.cpp new file mode 100644 index 000000000..3da3410b3 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.cpp @@ -0,0 +1,182 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +#include "GateOptnForceCommonTruncatedExp.h" +#include "G4ILawCommonTruncatedExp.hh" +#include "G4ILawForceFreeFlight.hh" +#include "G4TransportationManager.hh" + +#include "Randomize.hh" +#include "G4BiasingProcessInterface.hh" + +GateOptnForceCommonTruncatedExp::GateOptnForceCommonTruncatedExp(G4String name) + : G4VBiasingOperation(name), + fNumberOfSharing(0), + fProcessToApply(nullptr), + fInteractionOccured(false), + fMaximumDistance(-1.0) +{ + fCommonTruncatedExpLaw = new G4ILawCommonTruncatedExp("ExpLawForOperation"+name); + fForceFreeFlightLaw = new G4ILawForceFreeFlight ("FFFLawForOperation"+name); + + fTotalCrossSection = 0.0; +} + +GateOptnForceCommonTruncatedExp::~GateOptnForceCommonTruncatedExp() +{ + if ( fCommonTruncatedExpLaw ) delete fCommonTruncatedExpLaw; + if ( fForceFreeFlightLaw ) delete fForceFreeFlightLaw; +} + +const G4VBiasingInteractionLaw* GateOptnForceCommonTruncatedExp:: +ProvideOccurenceBiasingInteractionLaw( const G4BiasingProcessInterface* callingProcess, + G4ForceCondition& proposeForceCondition ) +{ + if ( callingProcess->GetWrappedProcess() == fProcessToApply ) + { + proposeForceCondition = Forced; + return fCommonTruncatedExpLaw; + } + else + { + proposeForceCondition = Forced; + return fForceFreeFlightLaw; + } +} + + +G4GPILSelection GateOptnForceCommonTruncatedExp::ProposeGPILSelection( const G4GPILSelection ) +{ + return NotCandidateForSelection; +} + + +G4VParticleChange* GateOptnForceCommonTruncatedExp::ApplyFinalStateBiasing( const G4BiasingProcessInterface* callingProcess, + const G4Track* track, + const G4Step* step, + G4bool& forceFinalState ) +{ + if ( callingProcess->GetWrappedProcess() != fProcessToApply ) + { + forceFinalState = true; + fDummyParticleChange.Initialize( *track ); + return &fDummyParticleChange; + } + if ( fInteractionOccured ) + { + forceFinalState = true; + fDummyParticleChange.Initialize( *track ); + return &fDummyParticleChange; + } + + // -- checks if process won the GPIL race: + G4double processGPIL = callingProcess->GetPostStepGPIL() < callingProcess->GetAlongStepGPIL() ? + callingProcess->GetPostStepGPIL() : callingProcess->GetAlongStepGPIL() ; + if ( processGPIL <= step->GetStepLength() ) + { + // -- if process won, wrapped process produces the final state. + // -- In this case, the weight for occurrence biasing is applied + // -- by the callingProcess, at exit of present method. This is + // -- selected by "forceFinalState = false": + forceFinalState = false; + fInteractionOccured = true; + return callingProcess->GetWrappedProcess()->PostStepDoIt( *track, *step ); + } + else + { + forceFinalState = true; + fDummyParticleChange.Initialize( *track ); + return &fDummyParticleChange; + } +} + + +void GateOptnForceCommonTruncatedExp::AddCrossSection( const G4VProcess* process, G4double crossSection ) +{ + fTotalCrossSection += crossSection; + fCrossSections[process] = crossSection; + fNumberOfSharing = fCrossSections.size(); +} + + +void GateOptnForceCommonTruncatedExp::Initialize( const G4Track* track ) +{ + fCrossSections.clear(); + fTotalCrossSection = 0.0; + fNumberOfSharing = 0; + fProcessToApply = 0; + fInteractionOccured = false; + fInitialMomentum = track->GetMomentum(); + + G4VSolid* currentSolid = track->GetVolume()->GetLogicalVolume()->GetSolid(); + G4ThreeVector localPosition = (G4TransportationManager::GetTransportationManager()-> + GetNavigatorForTracking()-> + GetGlobalToLocalTransform()).TransformPoint(track->GetPosition()); + G4ThreeVector localDirection = (G4TransportationManager::GetTransportationManager()-> + GetNavigatorForTracking()-> + GetGlobalToLocalTransform()).TransformAxis(track->GetMomentumDirection()); + fMaximumDistance = currentSolid->DistanceToOut(localPosition, localDirection); + if ( fMaximumDistance <= DBL_MIN ) fMaximumDistance = 0.0; + fCommonTruncatedExpLaw->SetMaximumDistance( fMaximumDistance ); +} + + +void GateOptnForceCommonTruncatedExp::UpdateForStep( const G4Step* step ) +{ + fCrossSections.clear(); + fTotalCrossSection = 0.0; + fNumberOfSharing = 0; + fProcessToApply = 0; + + fCommonTruncatedExpLaw->UpdateForStep( step->GetStepLength() ); + fMaximumDistance = fCommonTruncatedExpLaw->GetMaximumDistance(); +} + + +void GateOptnForceCommonTruncatedExp::Sample() +{ + fCommonTruncatedExpLaw->SetForceCrossSection( fTotalCrossSection ); + fCommonTruncatedExpLaw->Sample(); + ChooseProcessToApply(); + fCommonTruncatedExpLaw->SetSelectedProcessXSfraction(fCrossSections[fProcessToApply] / fTotalCrossSection); +} + + +void GateOptnForceCommonTruncatedExp::ChooseProcessToApply() +{ + G4double sigmaRand = G4UniformRand() * fTotalCrossSection; + G4double sigmaSelect = 0.0; + for ( std::map< const G4VProcess*, G4double>::const_iterator it = fCrossSections.begin(); + it != fCrossSections.end(); + it++) + { + sigmaSelect += (*it).second; + if ( sigmaRand <= sigmaSelect ) + { + fProcessToApply = (*it).first; + break; + } + } +} diff --git a/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.h b/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.h new file mode 100644 index 000000000..ec15f1745 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.h @@ -0,0 +1,122 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// +//--------------------------------------------------------------- +// +// GateOptnForceCommonTruncatedExp +// +// Class Description: +// A G4VBiasingOperation physics-based biasing operation. It +// handles several processes together, biasing on the total +// cross-section of these processes (instead of biasing them +// individually). +// The biasing interaction law is a truncated exponential +// one, driven by the total cross-section and which extends +// in the range [0,L]. +// Process are registered with the AddCrossSection method. +// As cross-sections are all known at the end of the +// PostStepGPIL loop, the step limitation is made at the +// AlongStepGPIL level. +// +//--------------------------------------------------------------- +// Initial version Nov. 2013 M. Verderi + + +#ifndef GateOptnForceCommonTruncatedExp_h +#define GateOptnForceCommonTruncatedExp_h 1 + +#include "G4VBiasingOperation.hh" +#include "G4ThreeVector.hh" +#include "G4ParticleChangeForNothing.hh" +class G4ILawCommonTruncatedExp; +class G4ILawForceFreeFlight; +#include + +class GateOptnForceCommonTruncatedExp : public G4VBiasingOperation { +public: + // -- Constructor : + GateOptnForceCommonTruncatedExp(G4String name); + // -- destructor: + virtual ~GateOptnForceCommonTruncatedExp(); + +public: + // -- Methods from G4VBiasingOperation interface: + // ------------------------------------------- + // -- Used: + virtual const G4VBiasingInteractionLaw* ProvideOccurenceBiasingInteractionLaw( const G4BiasingProcessInterface*, G4ForceCondition& ); + virtual G4double ProposeAlongStepLimit( const G4BiasingProcessInterface* ) { return DBL_MAX; } + virtual G4GPILSelection ProposeGPILSelection( const G4GPILSelection processSelection ); + virtual G4VParticleChange* ApplyFinalStateBiasing( const G4BiasingProcessInterface*, + const G4Track*, + const G4Step*, + G4bool& ); + // -- Unused: + virtual G4double DistanceToApplyOperation( const G4Track*, + G4double, + G4ForceCondition*) {return DBL_MAX;} + virtual G4VParticleChange* GenerateBiasingFinalState( const G4Track*, + const G4Step* ) {return 0;} + +public: + // -- Additional methods, specific to this class: + // ---------------------------------------------- + // -- return concrete type of interaction laws: + G4ILawCommonTruncatedExp* GetCommonTruncatedExpLaw() + { + return fCommonTruncatedExpLaw; + } + G4ILawForceFreeFlight* GetForceFreeFlightLaw() + { + return fForceFreeFlightLaw; + } + + void Initialize( const G4Track* ); + void UpdateForStep( const G4Step* ); + void Sample(); + const G4ThreeVector& GetInitialMomentum() const { return fInitialMomentum; } + G4double GetMaximumDistance() const { return fMaximumDistance; } + void ChooseProcessToApply(); + const G4VProcess* GetProcessToApply() const { return fProcessToApply; } + void AddCrossSection( const G4VProcess*, G4double ); + size_t GetNumberOfSharing() const { return fNumberOfSharing;} + void SetInteractionOccured( G4bool b ) { fInteractionOccured = b; } + G4bool GetInteractionOccured() const { return fInteractionOccured; } + +private: + G4ILawCommonTruncatedExp* fCommonTruncatedExpLaw; + G4ILawForceFreeFlight* fForceFreeFlightLaw; + G4double fTotalCrossSection; + std::map < const G4VProcess*, G4double > fCrossSections; + size_t fNumberOfSharing; + const G4VProcess* fProcessToApply; + G4bool fInteractionOccured; + G4ThreeVector fInitialMomentum; + G4double fMaximumDistance; + G4ParticleChangeForNothing fDummyParticleChange; +}; + +#endif diff --git a/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.cpp b/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.cpp new file mode 100644 index 000000000..686501dd6 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.cpp @@ -0,0 +1,104 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +#include "GateOptnForceFreeFlight.h" +#include "G4ILawForceFreeFlight.hh" +#include "G4BiasingProcessInterface.hh" +#include "G4Step.hh" + + + +GateOptnForceFreeFlight::GateOptnForceFreeFlight(G4String name) + : G4VBiasingOperation ( name ), + fCumulatedWeightChange ( -1.0 ), + fInitialTrackWeight ( -1.0 ), + fOperationComplete ( true ) +{ + fForceFreeFlightInteractionLaw = new G4ILawForceFreeFlight("LawForOperation"+name); +} + +GateOptnForceFreeFlight::~GateOptnForceFreeFlight() +{ + if ( fForceFreeFlightInteractionLaw ) delete fForceFreeFlightInteractionLaw; +} + +const G4VBiasingInteractionLaw* GateOptnForceFreeFlight::ProvideOccurenceBiasingInteractionLaw( const G4BiasingProcessInterface*, G4ForceCondition& proposeForceCondition ) +{ + fOperationComplete = false; + proposeForceCondition = Forced; + return fForceFreeFlightInteractionLaw; +} + + +G4VParticleChange* GateOptnForceFreeFlight::ApplyFinalStateBiasing( const G4BiasingProcessInterface* callingProcess, + const G4Track* track, + const G4Step* step, + G4bool& forceFinalState) +{ + // -- If the track is reaching the volume boundary, its free flight ends. In this case, its zero + // -- weight is brought back to non-zero value: its initial weight is restored by the first + // -- ApplyFinalStateBiasing operation called, and the weight for force free flight is applied + // -- is applied by each operation. + // -- If the track is not reaching the volume boundary, it zero weight flight continues. + + fParticleChange.Initialize( *track ); + forceFinalState = true; + if ( step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary ) + { + // -- Sanity checks: + if ( fInitialTrackWeight <= DBL_MIN ) + { + G4ExceptionDescription ed; + ed << " Initial track weight is null ! " << G4endl; + G4Exception(" GateOptnForceFreeFlight::ApplyFinalStateBiasing(...)", + "BIAS.GEN.05", + JustWarning, + ed); + } + if ( fCumulatedWeightChange <= DBL_MIN ) + { + G4ExceptionDescription ed; + ed << " Cumulated weight is null ! " << G4endl; + G4Exception(" GateOptnForceFreeFlight::ApplyFinalStateBiasing(...)", + "BIAS.GEN.06", + JustWarning, + ed); + } + + G4double proposedWeight = track->GetWeight(); + if ( callingProcess->GetIsFirstPostStepDoItInterface() ) proposedWeight = fCumulatedWeightChange * fInitialTrackWeight; + else proposedWeight *= fCumulatedWeightChange; + fParticleChange.ProposeWeight(proposedWeight); + fOperationComplete = true; + } + + return &fParticleChange; +} + + +void GateOptnForceFreeFlight::AlongMoveBy( const G4BiasingProcessInterface*, const G4Step*, G4double weightChange ) +{ + fCumulatedWeightChange *= weightChange; +} diff --git a/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.h b/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.h new file mode 100644 index 000000000..e7dffb49a --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.h @@ -0,0 +1,97 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// +//--------------------------------------------------------------- +// +// GateOptnForceFreeFlight +// +// Class Description: +// A G4VBiasingOperation physics-based biasing operation. +// If forces the physics process to not act on the track. +// In this implementation (meant for the ForceCollision +// operator) the free flight is done under zero weight for +// the track, and the action is meant to accumulate the weight +// change for making this uninteracting flight, +// cumulatedWeightChange. +// When the track reaches the current volume boundary, its +// weight is restored with value : +// initialWeight * cumulatedWeightChange +// +//--------------------------------------------------------------- +// Initial version Nov. 2013 M. Verderi +#ifndef GateOptnForceFreeFlight_h +#define GateOptnForceFreeFlight_h 1 + +#include "G4VBiasingOperation.hh" +#include "G4ForceCondition.hh" +#include "G4ParticleChange.hh" // -- §§ should add a dedicated "weight change only" particle change +class G4ILawForceFreeFlight; + + +class GateOptnForceFreeFlight : public G4VBiasingOperation { +public: + // -- Constructor : + GateOptnForceFreeFlight(G4String name); + // -- destructor: + virtual ~GateOptnForceFreeFlight(); + +public: + // -- Methods from G4VBiasingOperation interface: + // ------------------------------------------- + // -- Used: + virtual const G4VBiasingInteractionLaw* ProvideOccurenceBiasingInteractionLaw( const G4BiasingProcessInterface*, G4ForceCondition& ); + virtual void AlongMoveBy( const G4BiasingProcessInterface*, const G4Step*, G4double ); + virtual G4VParticleChange* ApplyFinalStateBiasing( const G4BiasingProcessInterface*, const G4Track*, const G4Step*, G4bool&); + + // -- Unused: + virtual G4double DistanceToApplyOperation( const G4Track*, + G4double, + G4ForceCondition*) {return DBL_MAX;} + virtual G4VParticleChange* GenerateBiasingFinalState( const G4Track*, + const G4Step* ) {return 0;} + + +public: + // -- Additional methods, specific to this class: + // ---------------------------------------------- + // -- return concrete type of interaction law: + G4ILawForceFreeFlight* GetForceFreeFlightLaw() { + return fForceFreeFlightInteractionLaw; + } + // -- initialization for weight: + void ResetInitialTrackWeight(G4double w) {fInitialTrackWeight = w; fCumulatedWeightChange = 1.0;} + G4bool OperationComplete() const { return fOperationComplete; } + +private: + G4ILawForceFreeFlight* fForceFreeFlightInteractionLaw; + G4double fCumulatedWeightChange, + fInitialTrackWeight; + G4ParticleChange fParticleChange; + G4bool fOperationComplete; +}; + +#endif diff --git a/core/opengate_core/opengate_lib/GateOptrForceCollision.cpp b/core/opengate_core/opengate_lib/GateOptrForceCollision.cpp new file mode 100644 index 000000000..cc5294f8e --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptrForceCollision.cpp @@ -0,0 +1,442 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +#include "GateOptrForceCollision.h" +#include "GateOptrForceCollisionTrackData.h" +#include "G4BiasingProcessInterface.hh" +#include "G4PhysicsModelCatalog.hh" + +#include "GateOptnForceCommonTruncatedExp.h" +#include "G4ILawCommonTruncatedExp.hh" +#include "GateOptnForceFreeFlight.h" +#include "GateOptnCloning.h" + +#include "G4Step.hh" +#include "G4StepPoint.hh" +#include "G4VProcess.hh" + +#include "G4ParticleDefinition.hh" +#include "G4ParticleTable.hh" + +#include "G4SystemOfUnits.hh" +#include "GateHelpersDict.h" +#include "GateHelpersImage.h" + +// -- §§ consider calling other constructor, thanks to C++11 + +GateOptrForceCollision::GateOptrForceCollision(py::dict &user_info) + : G4VBiasingOperator("forceCollisionActor"),GateVActor(user_info, false), + fForceCollisionModelID(G4PhysicsModelCatalog::GetModelID("model_GenBiasForceCollision")), + fCurrentTrack(nullptr), + fCurrentTrackData(nullptr), + fInitialTrackWeight(-1.0), + fSetup(true) +{ + fSharedForceInteractionOperation = new GateOptnForceCommonTruncatedExp("SharedForceInteraction"); + fCloningOperation = new GateOptnCloning("Cloning"); +} + + +void GateOptrForceCollision::AttachAllLogicalDaughtersVolumes( + G4LogicalVolume *volume) { + AttachTo(volume); + G4int nbOfDaughters = volume->GetNoDaughters(); + if (nbOfDaughters > 0) { + for (int i = 0; i < nbOfDaughters; i++) { + G4LogicalVolume *logicalDaughtersVolume = + volume->GetDaughter(i)->GetLogicalVolume(); + AttachAllLogicalDaughtersVolumes(logicalDaughtersVolume); + } + } +} + +GateOptrForceCollision::~GateOptrForceCollision() +{ + for ( std::map< const G4BiasingProcessInterface*, GateOptnForceFreeFlight* >::iterator it = fFreeFlightOperations.begin() ; + it != fFreeFlightOperations.end() ; + it++ ) delete (*it).second; + delete fSharedForceInteractionOperation; + delete fCloningOperation; +} + +void GateOptrForceCollision::InitializeUserInfo(py::dict &user_info) { + // IMPORTANT: call the base class method + GateVActor::InitializeUserInfo(user_info); + G4String particleToBias = "gamma"; + G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable(); + for (G4int i = 0; i < particleTable->size(); ++i) { + G4ParticleDefinition *particle = particleTable->GetParticle(i); + G4String particleName = particle->GetParticleName(); + if (particleName == particleToBias) + fParticleToBias = particle; + } +} + + + +void GateOptrForceCollision::Configure() +{ + // -- build free flight operations: + ConfigureForWorker(); +} + + +void GateOptrForceCollision::ConfigureForWorker() +{ + // -- start by remembering processes under biasing, and create needed biasing operations: + if ( fSetup ) + { + const G4ProcessManager* processManager = fParticleToBias->GetProcessManager(); + const G4BiasingProcessSharedData* interfaceProcessSharedData = G4BiasingProcessInterface::GetSharedData( processManager ); + if ( interfaceProcessSharedData ) // -- sharedData tested, as is can happen a user attaches an operator + // -- to a volume but without defining BiasingProcessInterface processes. + { + for ( size_t i = 0 ; i < (interfaceProcessSharedData->GetPhysicsBiasingProcessInterfaces()).size(); i++ ) + { + const G4BiasingProcessInterface* wrapperProcess = + (interfaceProcessSharedData->GetPhysicsBiasingProcessInterfaces())[i]; + G4String operationName = "FreeFlight-"+wrapperProcess->GetWrappedProcess()->GetProcessName(); + fFreeFlightOperations[wrapperProcess] = new GateOptnForceFreeFlight(operationName); + } + } + fSetup = false; + } +} + + +void GateOptrForceCollision::StartRun() +{ + G4LogicalVolume *biasingVolume = + G4LogicalVolumeStore::GetInstance()->GetVolume(fAttachedToVolumeName); + + // Here we need to attach all the daughters and daughters of daughters (...) + // to the biasing operator. To do that, I use the function + // AttachAllLogicalDaughtersVolumes. + AttachAllLogicalDaughtersVolumes(biasingVolume); +} + + +G4VBiasingOperation* GateOptrForceCollision::ProposeOccurenceBiasingOperation(const G4Track* track, const G4BiasingProcessInterface* callingProcess) +{ + // -- does nothing if particle is not of requested type: + if ( track->GetDefinition() != fParticleToBias ) return 0; + + // -- trying to get auxiliary track data... + if ( fCurrentTrackData == nullptr ) + { + // ... and if the track has no aux. track data, it means its biasing is not started yet (note that cloning is to happen first): + fCurrentTrackData = (GateOptrForceCollisionTrackData*)(track->GetAuxiliaryTrackInformation(fForceCollisionModelID)); + if ( fCurrentTrackData == nullptr ) return nullptr; + } + + + // -- Send force free flight to the callingProcess: + // ------------------------------------------------ + // -- The track has been cloned in the previous step, it has now to be + // -- forced for a free flight. + // -- This track will fly with 0.0 weight during its forced flight: + // -- it is to forbid the double counting with the force interaction track. + // -- Its weight is restored at the end of its free flight, this weight + // -- being its initial weight * the weight for the free flight travel, + // -- this last one being per process. The initial weight is common, and is + // -- arbitrary asked to the first operation to take care of it. + if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeFreeFlight ) + { + GateOptnForceFreeFlight* operation = fFreeFlightOperations[callingProcess]; + if ( callingProcess->GetWrappedProcess()->GetCurrentInteractionLength() < DBL_MAX/10. ) + { + // -- the initial track weight will be restored only by the first DoIt free flight: + operation->ResetInitialTrackWeight(fInitialTrackWeight); + return operation; + } + else + { + return nullptr; + } + } + + + // -- Send force interaction operation to the callingProcess: + // ---------------------------------------------------------- + // -- at this level, a copy of the track entering the volume was + // -- generated (borned) earlier. This copy will make the forced + // -- interaction in the volume. + if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced ) + { + // -- Remember if this calling process is the first of the physics wrapper in + // -- the PostStepGPIL loop (using default argument of method below): + G4bool isFirstPhysGPIL = callingProcess-> GetIsFirstPostStepGPILInterface(); + + // -- [*first process*] Initialize or update force interaction operation: + if ( isFirstPhysGPIL ) + { + // -- first step of cloned track, initialize the forced interaction operation: + if ( track->GetCurrentStepNumber() == 1 ) fSharedForceInteractionOperation->Initialize( track ); + else + { + if ( fSharedForceInteractionOperation->GetInitialMomentum() != track->GetMomentum() ) + { + // -- means that some other physics process, not under control of the forced interaction operation, + // -- has occured, need to re-initialize the operation as distance to boundary has changed. + // -- [ Note the re-initialization is only possible for a Markovian law. ] + fSharedForceInteractionOperation->Initialize( track ); + } + else + { + // -- means that some other non-physics process (biasing or not, like step limit), has occured, + // -- but track conserves its momentum direction, only need to reduced the maximum distance for + // -- forced interaction. + // -- [ Note the update is only possible for a Markovian law. ] + fSharedForceInteractionOperation->UpdateForStep( track->GetStep() ); + } + } + } + + // -- [*all processes*] Sanity check : it may happen in limit cases that distance to + // -- out is zero, weight would be infinite in this case: abort forced interaction + // -- and abandon biasing. + if ( fSharedForceInteractionOperation->GetMaximumDistance() < DBL_MIN ) + { + fCurrentTrackData->Reset(); + return 0; + } + + // -- [* first process*] collect cross-sections and sample force law to determine interaction length + // -- and winning process: + if ( isFirstPhysGPIL ) + { + // -- collect cross-sections: + // -- ( Remember that the first of the G4BiasingProcessInterface triggers the update + // -- of these cross-sections ) + const G4BiasingProcessSharedData* sharedData = callingProcess->GetSharedData(); + for ( size_t i = 0 ; i < (sharedData->GetPhysicsBiasingProcessInterfaces()).size() ; i++ ) + { + const G4BiasingProcessInterface* wrapperProcess = ( sharedData->GetPhysicsBiasingProcessInterfaces() )[i]; + G4double interactionLength = wrapperProcess->GetWrappedProcess()->GetCurrentInteractionLength(); + // -- keep only well defined cross-sections, other processes are ignored. These are not pathological cases + // -- but cases where a threhold effect par example (eg pair creation) may be at play. (**!**) + if ( interactionLength < DBL_MAX/10. ) + fSharedForceInteractionOperation->AddCrossSection( wrapperProcess->GetWrappedProcess(), 1.0/interactionLength ); + } + // -- sample the shared law (interaction length, and winning process): + if ( fSharedForceInteractionOperation->GetNumberOfSharing() > 0 ) fSharedForceInteractionOperation->Sample(); + } + + // -- [*all processes*] Send operation for processes with well defined XS (see "**!**" above): + G4VBiasingOperation* operationToReturn = nullptr; + if ( callingProcess->GetWrappedProcess()->GetCurrentInteractionLength() < DBL_MAX/10. ) operationToReturn = fSharedForceInteractionOperation; + return operationToReturn; + + + } // -- end of "if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced )" + + + // -- other cases here: particle appearing in the volume by some + // -- previous interaction : we decide to not bias these. + return 0; + +} + + +G4VBiasingOperation* GateOptrForceCollision::ProposeNonPhysicsBiasingOperation(const G4Track* track, + const G4BiasingProcessInterface* /* callingProcess */) +{ + if ( track->GetDefinition() != fParticleToBias ) return nullptr; + + // -- Apply biasing scheme only to tracks entering the volume. + // -- A "cloning" is done: + // -- - the primary will be forced flight under a zero weight up to volume exit, + // -- where the weight will be restored with proper weight for free flight + // -- - the clone will be forced to interact in the volume. + if ( track->GetStep()->GetPreStepPoint()->GetStepStatus() == fGeomBoundary ) // -- §§§ extent to case of a track shoot on the boundary ? + { + // -- check that track is free of undergoing biasing scheme ( no biasing data, or no active active ) + // -- Get possible track data: + fCurrentTrackData = (GateOptrForceCollisionTrackData*)(track->GetAuxiliaryTrackInformation(fForceCollisionModelID)); + if ( fCurrentTrackData != nullptr ) + { + if ( fCurrentTrackData->IsFreeFromBiasing() ) + { + // -- takes "ownership" (some track data created before, left free, reuse it): + fCurrentTrackData->fForceCollisionOperator = this ; + } + else + { + // §§§ Would something be really wrong in this case ? Could this be that a process made a zero step ? + } + } + else + { + fCurrentTrackData = new GateOptrForceCollisionTrackData( this ); + track->SetAuxiliaryTrackInformation(fForceCollisionModelID, fCurrentTrackData); + } + fCurrentTrackData->fForceCollisionState = ForceCollisionState::toBeCloned; + fInitialTrackWeight = track->GetWeight(); + fCloningOperation->SetCloneWeights(0.0, fInitialTrackWeight); + return fCloningOperation; + } + + // -- + return nullptr; +} + + +G4VBiasingOperation* GateOptrForceCollision::ProposeFinalStateBiasingOperation(const G4Track*, const G4BiasingProcessInterface* callingProcess) +{ + // -- Propose at final state generation the same operation which was proposed at GPIL level, + // -- (which is either the force free flight or the force interaction operation). + // -- count on the process interface to collect this operation: + return callingProcess->GetCurrentOccurenceBiasingOperation(); +} + + +void GateOptrForceCollision::StartTracking( const G4Track* track ) +{ + fCurrentTrack = track; + fCurrentTrackData = nullptr; +} + + +void GateOptrForceCollision::EndTracking() +{ + // -- check for consistency, operator should have cleaned the track: + if ( fCurrentTrackData != nullptr ) + { + if ( !fCurrentTrackData->IsFreeFromBiasing() ) + { + if ( (fCurrentTrack->GetTrackStatus() == fStopAndKill) || (fCurrentTrack->GetTrackStatus() == fKillTrackAndSecondaries) ) + { + G4ExceptionDescription ed; + ed << "Current track deleted while under biasing by " << GetName() << ". Will result in inconsistencies."; + G4Exception(" GateOptrForceCollision::EndTracking()", + "BIAS.GEN.18", + JustWarning, + ed); + } + } + } +} + + +void GateOptrForceCollision::OperationApplied( const G4BiasingProcessInterface* callingProcess, + G4BiasingAppliedCase BAC, + G4VBiasingOperation* operationApplied, + const G4VParticleChange* ) +{ + + if ( fCurrentTrackData == nullptr ) + { + if ( BAC != BAC_None ) + { + G4ExceptionDescription ed; + ed << " Internal inconsistency : please submit bug report. " << G4endl; + G4Exception(" GateOptrForceCollision::OperationApplied(...)", + "BIAS.GEN.20.1", + JustWarning, + ed); + } + return; + } + + if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeCloned ) + { + fCurrentTrackData->fForceCollisionState = ForceCollisionState::toBeFreeFlight; + auto cloneData = new GateOptrForceCollisionTrackData( this ); + cloneData->fForceCollisionState = ForceCollisionState::toBeForced; + fCloningOperation->GetCloneTrack()->SetAuxiliaryTrackInformation(fForceCollisionModelID, cloneData); + } + else if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeFreeFlight ) + { + if ( fFreeFlightOperations[callingProcess]->OperationComplete() ) fCurrentTrackData->Reset(); // -- off biasing for this track + } + else if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced ) + { + if ( operationApplied != fSharedForceInteractionOperation ) + { + G4ExceptionDescription ed; + ed << " Internal inconsistency : please submit bug report. " << G4endl; + G4Exception(" GateOptrForceCollision::OperationApplied(...)", + "BIAS.GEN.20.2", + JustWarning, + ed); + } + if ( fSharedForceInteractionOperation->GetInteractionOccured() ) + { + if ( operationApplied != fSharedForceInteractionOperation ) + { + G4ExceptionDescription ed; + ed << " Internal inconsistency : please submit bug report. " << G4endl; + G4Exception(" GateOptrForceCollision::OperationApplied(...)", + "BIAS.GEN.20.3", + JustWarning, + ed); + } + } + } + else + { + if ( fCurrentTrackData->fForceCollisionState != ForceCollisionState::free ) + { + G4ExceptionDescription ed; + ed << " Internal inconsistency : please submit bug report. " << G4endl; + G4Exception(" GateOptrForceCollision::OperationApplied(...)", + "BIAS.GEN.20.4", + JustWarning, + ed); + } + } +} + + +void GateOptrForceCollision::OperationApplied( const G4BiasingProcessInterface* /*callingProcess*/, G4BiasingAppliedCase /*biasingCase*/, + G4VBiasingOperation* /*occurenceOperationApplied*/, G4double /*weightForOccurenceInteraction*/, + G4VBiasingOperation* finalStateOperationApplied, const G4VParticleChange* /*particleChangeProduced*/ ) +{ + + if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced ) + { + if ( finalStateOperationApplied != fSharedForceInteractionOperation ) + { + G4ExceptionDescription ed; + ed << " Internal inconsistency : please submit bug report. " << G4endl; + G4Exception(" GateOptrForceCollision::OperationApplied(...)", + "BIAS.GEN.20.5", + JustWarning, + ed); + } + if ( fSharedForceInteractionOperation->GetInteractionOccured() ) fCurrentTrackData->Reset(); // -- off biasing for this track + } + else + { + G4ExceptionDescription ed; + ed << " Internal inconsistency : please submit bug report. " << G4endl; + G4Exception(" GateOptrForceCollision::OperationApplied(...)", + "BIAS.GEN.20.6", + JustWarning, + ed); + } + +} + diff --git a/core/opengate_core/opengate_lib/GateOptrForceCollision.h b/core/opengate_core/opengate_lib/GateOptrForceCollision.h new file mode 100644 index 000000000..d720955bc --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptrForceCollision.h @@ -0,0 +1,105 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// +// --------------------------------------------------------------- +// +// GateOptrForceCollision +// +// Class Description: +// A G4VBiasingOperator that implements a "force collision" a la +// MCNP. This is meant for neutral particles. +// When the track enters the volume, it is cloned. One copy makes +// a forced free flight up to the volume exit. The other copy makes +// a forced collision inside the volume. +// +// --------------------------------------------------------------- +// Initial version Nov. 2013 M. Verderi + + +#ifndef GateOptrForceCollision_h +#define GateOptrForceCollision_h 1 + +#include "G4VBiasingOperator.hh" +class GateOptnForceFreeFlight; +class GateOptnForceCommonTruncatedExp; +class GateOptnCloning; +class G4VProcess; +class G4BiasingProcessInterface; +class G4ParticleDefinition; +#include +#include +#include "G4ThreeVector.hh" +#include "GateVActor.h" +class GateOptrForceCollisionTrackData; + +class GateOptrForceCollision : public G4VBiasingOperator,public GateVActor { +public: + GateOptrForceCollision(py::dict &user_info); + ~GateOptrForceCollision(); + +public: + // -- Mandatory from base class : + virtual G4VBiasingOperation* ProposeNonPhysicsBiasingOperation(const G4Track* track, const G4BiasingProcessInterface* callingProcess) final; + virtual G4VBiasingOperation* ProposeOccurenceBiasingOperation(const G4Track* track, const G4BiasingProcessInterface* callingProcess) final; + virtual G4VBiasingOperation* ProposeFinalStateBiasingOperation(const G4Track* track, const G4BiasingProcessInterface* callingProcess) final; + // -- optional methods from base class: +public: + virtual void Configure() final; + virtual void ConfigureForWorker() final; + virtual void StartRun() final; + virtual void StartTracking( const G4Track* track ) final; + virtual void ExitBiasing( const G4Track*, const G4BiasingProcessInterface* ) final {}; + virtual void EndTracking() final; + void InitializeUserInfo(py::dict &user_info) override; + +void AttachAllLogicalDaughtersVolumes(G4LogicalVolume *volume); + + // -- operation applied: + void OperationApplied( const G4BiasingProcessInterface* callingProcess, G4BiasingAppliedCase biasingCase, + G4VBiasingOperation* operationApplied, const G4VParticleChange* particleChangeProduced ) final; + void OperationApplied( const G4BiasingProcessInterface* callingProcess, G4BiasingAppliedCase biasingCase, + G4VBiasingOperation* occurenceOperationApplied, G4double weightForOccurenceInteraction, + G4VBiasingOperation* finalStateOperationApplied, const G4VParticleChange* particleChangeProduced ) final; + + +const G4String GetName(){ + return GateVActor::GetName(); +} + +public: + G4int fForceCollisionModelID; + const G4Track* fCurrentTrack; + GateOptrForceCollisionTrackData* fCurrentTrackData; + std::map< const G4BiasingProcessInterface*, GateOptnForceFreeFlight* > fFreeFlightOperations; + GateOptnForceCommonTruncatedExp* fSharedForceInteractionOperation; + GateOptnCloning* fCloningOperation; + G4double fInitialTrackWeight; + G4bool fSetup; + const G4ParticleDefinition* fParticleToBias; +}; + +#endif diff --git a/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.cpp b/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.cpp new file mode 100644 index 000000000..26d8374e4 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.cpp @@ -0,0 +1,74 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +#include "GateOptrForceCollisionTrackData.h" +#include "GateOptrForceCollision.h" + +GateOptrForceCollisionTrackData::GateOptrForceCollisionTrackData( const GateOptrForceCollision* optr ) +: G4VAuxiliaryTrackInformation(), + fForceCollisionOperator( optr ) +{ + fForceCollisionState = ForceCollisionState::free; +} + +GateOptrForceCollisionTrackData::~GateOptrForceCollisionTrackData() +{ + if ( fForceCollisionState != ForceCollisionState::free ) + { + G4ExceptionDescription ed; + ed << "Track deleted while under G4BOptrForceCollision biasing scheme of operator `"; + if ( fForceCollisionOperator == nullptr ) ed << "(none)"; else ed << "forceCollisionActor"; + ed <<"'. Will result in inconsistencies."; + G4Exception(" GateOptrForceCollisionTrackData::~GateOptrForceCollisionTrackData()", + "BIAS.GEN.19", + JustWarning, + ed); + } +} + +void GateOptrForceCollisionTrackData::Print() const +{ + G4cout << " GateOptrForceCollisionTrackData object : " << this << G4endl; + G4cout << " Force collision operator : "; if ( fForceCollisionOperator == nullptr ) G4cout << "(none)"; else G4cout << "forceCollisionActor"; G4cout << G4endl; + G4cout << " Force collision state : "; + switch ( fForceCollisionState ) + { + case ForceCollisionState::free : + G4cout << "free from biasing "; + break; + case ForceCollisionState::toBeCloned : + G4cout << "to be cloned "; + break; + case ForceCollisionState::toBeForced : + G4cout << "to be interaction forced "; + break; + case ForceCollisionState::toBeFreeFlight : + G4cout << "to be free flight forced (under weight = 0) "; + break; + default: + break; + } + G4cout << G4endl; +} diff --git a/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.h b/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.h new file mode 100644 index 000000000..8e0ecaf4f --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.h @@ -0,0 +1,80 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// +// -------------------------------------------------------------------- +// GEANT 4 class header file +// +// Class Description: +// Extends G4Track properties with information needed for the +// force collision biasing operator. +// The G4BOptrForceCollision class is made friend of this one in +// order to keep unexposed to public most of the data members, as +// they are used to control the logic. +// +// ------------------ GateOptrForceCollisionTrackData ------------------ +// +// Author: M.Verderi (LLR), October 2015 +// +// -------------------------------------------------------------------- + +#ifndef GateOptrForceCollisionTrackData_h +#define GateOptrForceCollisionTrackData_h + +class GateOptrForceCollision; +#include "G4VAuxiliaryTrackInformation.hh" + +enum class ForceCollisionState { free, toBeCloned, toBeForced, toBeFreeFlight }; + +class GateOptrForceCollisionTrackData : public G4VAuxiliaryTrackInformation { + +friend class GateOptrForceCollision; + +public: + GateOptrForceCollisionTrackData( const GateOptrForceCollision* ); + ~GateOptrForceCollisionTrackData(); + + // -- from base class: + void Print() const; + + // -- Get methods: + G4bool IsFreeFromBiasing() const + { return ( fForceCollisionState == ForceCollisionState::free);} + // -- no set methods are provided : sets are made under exclusive control of G4BOptrForceCollision objects through friendness. + +private: + const GateOptrForceCollision* fForceCollisionOperator; + ForceCollisionState fForceCollisionState; + + void Reset() + { + fForceCollisionOperator = nullptr; + fForceCollisionState = ForceCollisionState::free; + } + +}; + +#endif diff --git a/core/opengate_core/opengate_lib/pyGateOptrForceCollision.cpp b/core/opengate_core/opengate_lib/pyGateOptrForceCollision.cpp new file mode 100644 index 000000000..245ddadd8 --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateOptrForceCollision.cpp @@ -0,0 +1,19 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ +#include + +namespace py = pybind11; +#include "G4VBiasingOperator.hh" +#include "GateOptrForceCollision.h" + +void init_GateOptrForceCollision(py::module &m) { + + py::class_>( + m, "GateOptrForceCollision") + .def(py::init()); +} diff --git a/opengate/actors/miscactors.py b/opengate/actors/miscactors.py index b263a16ea..dab6e7ca2 100644 --- a/opengate/actors/miscactors.py +++ b/opengate/actors/miscactors.py @@ -396,7 +396,7 @@ def _setter_hook_particles(self, value): return list(value) -class SplittingActorBase(ActorBase): +class GenericBiasingActorBase(ActorBase): """ Actors based on the G4GenericBiasing class of GEANT4. This class provides tools to interact with GEANT4 processes during a simulation, allowing direct modification of process properties. Additionally, it enables non-physics-based @@ -409,7 +409,6 @@ class SplittingActorBase(ActorBase): splitting_factor: int bias_primary_only: bool bias_only_once: bool - particles: list user_info_defaults = { "splitting_factor": ( @@ -421,28 +420,19 @@ class SplittingActorBase(ActorBase): "bias_primary_only": ( True, { - "doc": "If true, the splitting mechanism is applied only to particles with a ParentID of 1", + "doc": "If true, the biasing mechanism is applied only to particles with a ParentID of 1", }, ), "bias_only_once": ( True, { - "doc": "If true, the splitting mechanism is applied only once per particle history", - }, - ), - "particles": ( - [ - "all", - ], - { - "doc": "Specifies the particles to split. The default value, all, includes all particles", - "setter_hook": _setter_hook_particles, + "doc": "If true, the biasing mechanism is applied only once per particle history", }, ), } -class ComptSplittingActor(SplittingActorBase, g4.GateOptrComptSplittingActor): +class ComptSplittingActor(GenericBiasingActorBase, g4.GateOptrComptSplittingActor): """ This splitting actor enables process-based splitting specifically for Compton interactions. Each time a Compton process occurs, its behavior is modified by generating multiple Compton scattering tracks @@ -457,6 +447,8 @@ class ComptSplittingActor(SplittingActorBase, g4.GateOptrComptSplittingActor): rotation_vector_director: bool vector_director: list max_theta: float + mode : str + particles: list user_info_defaults = { "min_weight_of_particle": ( @@ -491,22 +483,24 @@ class ComptSplittingActor(SplittingActorBase, g4.GateOptrComptSplittingActor): ), } - processes = ("compt",) + processes = (("compt",)) + particles = ("gamma") + mode = "physics" def __init__(self, *args, **kwargs): - SplittingActorBase.__init__(self, *args, **kwargs) + GenericBiasingActorBase.__init__(self, *args, **kwargs) self.__initcpp__() def __initcpp__(self): g4.GateOptrComptSplittingActor.__init__(self, {"name": self.name}) def initialize(self): - SplittingActorBase.initialize(self) + GenericBiasingActorBase.initialize(self) self.InitializeUserInfo(self.user_info) self.InitializeCpp() -class BremSplittingActor(SplittingActorBase, g4.GateBOptrBremSplittingActor): +class BremSplittingActor(GenericBiasingActorBase, g4.GateBOptrBremSplittingActor): """ This splitting actor enables process-based splitting specifically for bremsstrahlung process. Each time a Brem process occurs, its behavior is modified by generating multiple secondary Brem scattering tracks @@ -515,36 +509,57 @@ class BremSplittingActor(SplittingActorBase, g4.GateBOptrBremSplittingActor): # hints for IDE processes: list + particles: list + mode :list - user_info_defaults = { - "processes": ( - ["eBrem"], - { - "doc": "Specifies the process split by this actor. This parameter is set to eBrem, as the actor is specifically developed for this process. It is recommended not to modify this setting.", - }, - ), - } + mode = "physics" + particles = ("e+","e-",) + processes = (("eBrem",), ("eBrem",)) - processes = ("eBrem",) def __init__(self, *args, **kwargs): - SplittingActorBase.__init__(self, *args, **kwargs) + GenericBiasingActorBase.__init__(self, *args, **kwargs) self.__initcpp__() def __initcpp__(self): g4.GateBOptrBremSplittingActor.__init__(self, {"name": self.name}) def initialize(self): - SplittingActorBase.initialize(self) + GenericBiasingActorBase.initialize(self) self.InitializeUserInfo(self.user_info) self.InitializeCpp() +class ForceCollisionActor(GenericBiasingActorBase, g4.GateOptrForceCollision): + particles: str + processes = list + mode : str + + + mode = "both" + particles = ("gamma",) + processes = [("compt","phot","conv","Rayl")] + + def __init__(self, *args, **kwargs): + GenericBiasingActorBase.__init__(self, *args, **kwargs) + self.__initcpp__() + + def __initcpp__(self): + g4.GateOptrForceCollision.__init__(self, {"name": self.name}) + + def initialize(self): + GenericBiasingActorBase.initialize(self) + self.InitializeUserInfo(self.user_info) + self.InitializeCpp() + + + process_cls(ActorOutputStatisticsActor) process_cls(SimulationStatisticsActor) process_cls(KillActor) process_cls(ActorOutputKillAccordingProcessesActor) process_cls(KillAccordingProcessesActor) -process_cls(SplittingActorBase) +process_cls(GenericBiasingActorBase) process_cls(ComptSplittingActor) process_cls(BremSplittingActor) +process_cls(ForceCollisionActor) diff --git a/opengate/engines.py b/opengate/engines.py index eb8fe6954..e3274483c 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -390,13 +390,21 @@ def initialize_g4_em_parameters(self): def initialize_physics_biasing(self): # get a dictionary {particle:[processes]} particles_processes = self.physics_manager.get_biasing_particles_and_processes() - # check if there are any processes requested for any particle + print(particles_processes) if any([len(v) > 0 for v in particles_processes.values()]): g4_biasing_physics = g4.G4GenericBiasingPhysics() for particle, processes in particles_processes.items(): - if len(processes) > 0: - g4_biasing_physics.PhysicsBias(particle, processes) + process_list = processes[0] + mode = processes[1] + if mode != "non physics": + if len(process_list) > 0: + if mode == "physics" : + g4_biasing_physics.PhysicsBias(particle, process_list) + elif mode == "both": + g4_biasing_physics.Bias(particle,process_list) + else : + g4_biasing_physics.NonPhysicsBias(process_list) self.g4_physics_list.RegisterPhysics(g4_biasing_physics) # This function deals with calling the parse function diff --git a/opengate/managers.py b/opengate/managers.py index e2590f9ca..dff98faf2 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -91,9 +91,10 @@ SimulationStatisticsActor, KillActor, KillAccordingProcessesActor, - SplittingActorBase, + GenericBiasingActorBase, ComptSplittingActor, BremSplittingActor, + ForceCollisionActor, ) from .actors.digitizers import ( DigitizerAdderActor, @@ -128,6 +129,7 @@ "KillAccordingProcessesActor": KillAccordingProcessesActor, "BremSplittingActor": BremSplittingActor, "ComptSplittingActor": ComptSplittingActor, + "ForceCollisionActor": ForceCollisionActor, "DigitizerAdderActor": DigitizerAdderActor, "DigitizerBlurringActor": DigitizerBlurringActor, "DigitizerSpatialBlurringActor": DigitizerSpatialBlurringActor, @@ -863,32 +865,33 @@ def get_biasing_particles_and_processes(self): all_particles = charged_particles.union({"gamma"}) # create a dictionary with sets as entries (to ensure uniqueness) - particles_processes = dict([(p, set()) for p in all_particles]) + particles_processes= dict([(p, set()) for p in all_particles]) for actor in self.simulation.actor_manager.actors.values(): - if isinstance(actor, SplittingActorBase): + if isinstance(actor, GenericBiasingActorBase): particles = set() - if "all" in actor.particles: - particles.update(all_particles) - elif "all_charged" in actor.particles: - particles.update(charged_particles) - else: - for particle in actor.particles: - p_ = translate_particle_name_gate_to_geant4(particle) - if p_ in all_particles: - particles.add(p_) - else: - fatal( - f"Biasing actor {actor.name} wants to apply a bias to particle '{p_}'. " - f"This is not possible. Compatible particles are: {list(all_particles)}. " - ) - for p in particles: - particles_processes[p].update(actor.processes) + mode = actor.mode + for particle in actor.particles: + p_ = translate_particle_name_gate_to_geant4(particle) + if p_ in all_particles: + particles.add(p_) + else: + fatal( + f"Biasing actor {actor.name} wants to apply a bias to particle '{p_}'. " + f"This is not possible. Compatible particles are: {list(all_particles)}. " + ) + for i,p in enumerate(particles): + if mode != "non physics": + print(actor.processes) + particles_processes[p].update(actor.processes[i]) + else : + particles_processes[p].update([None]) # convert the dictionary entries back from set to list + print(particles_processes.items()) return dict( [ - (particle, list(processes)) + (particle, [list(processes),mode]) for particle, processes in particles_processes.items() ] ) From 6981544215aa021ef0a9025379775d75d5f6641c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 5 Dec 2024 15:35:29 +0000 Subject: [PATCH 2/3] [pre-commit.ci] Automatic python and c++ formatting --- .../opengate_lib/GateOptnCloning.cpp | 26 +- .../opengate_lib/GateOptnCloning.h | 49 +- .../GateOptnForceCommonTruncatedExp.cpp | 217 +++---- .../GateOptnForceCommonTruncatedExp.h | 94 +-- .../opengate_lib/GateOptnForceFreeFlight.cpp | 108 ++-- .../opengate_lib/GateOptnForceFreeFlight.h | 48 +- .../opengate_lib/GateOptrForceCollision.cpp | 573 +++++++++--------- .../opengate_lib/GateOptrForceCollision.h | 75 ++- .../GateOptrForceCollisionTrackData.cpp | 77 +-- .../GateOptrForceCollisionTrackData.h | 30 +- opengate/actors/miscactors.py | 20 +- opengate/engines.py | 6 +- opengate/managers.py | 8 +- 13 files changed, 659 insertions(+), 672 deletions(-) diff --git a/core/opengate_core/opengate_lib/GateOptnCloning.cpp b/core/opengate_core/opengate_lib/GateOptnCloning.cpp index 91e5e5803..9242484b5 100644 --- a/core/opengate_core/opengate_lib/GateOptnCloning.cpp +++ b/core/opengate_core/opengate_lib/GateOptnCloning.cpp @@ -25,27 +25,21 @@ // #include "GateOptnCloning.h" - GateOptnCloning::GateOptnCloning(G4String name) - : G4VBiasingOperation( name ), - fClone1W ( -1.0 ), - fClone2W ( -1.0 ), - fParticleChange(), - fCloneTrack ( nullptr ) -{} + : G4VBiasingOperation(name), fClone1W(-1.0), fClone2W(-1.0), + fParticleChange(), fCloneTrack(nullptr) {} -GateOptnCloning::~GateOptnCloning() -{} +GateOptnCloning::~GateOptnCloning() {} -G4VParticleChange* GateOptnCloning::GenerateBiasingFinalState( const G4Track* track, - const G4Step* ) -{ +G4VParticleChange * +GateOptnCloning::GenerateBiasingFinalState(const G4Track *track, + const G4Step *) { fParticleChange.Initialize(*track); - fParticleChange.ProposeParentWeight( fClone1W ); + fParticleChange.ProposeParentWeight(fClone1W); fParticleChange.SetSecondaryWeightByProcess(true); fParticleChange.SetNumberOfSecondaries(1); - fCloneTrack = new G4Track( *track ); - fCloneTrack->SetWeight( fClone2W ); - fParticleChange.AddSecondary( fCloneTrack ); + fCloneTrack = new G4Track(*track); + fCloneTrack->SetWeight(fClone2W); + fParticleChange.AddSecondary(fCloneTrack); return &fParticleChange; } diff --git a/core/opengate_core/opengate_lib/GateOptnCloning.h b/core/opengate_core/opengate_lib/GateOptnCloning.h index fd6033fbd..8a27ad2fa 100644 --- a/core/opengate_core/opengate_lib/GateOptnCloning.h +++ b/core/opengate_core/opengate_lib/GateOptnCloning.h @@ -32,17 +32,16 @@ // Class Description: // A G4VBiasingOperation to clone a track, allowing to set // weight arbitrary weights. -// +// // //--------------------------------------------------------------- // Initial version Nov. 2013 M. Verderi - #ifndef GateOptnCloning_h #define GateOptnCloning_h 1 -#include "G4VBiasingOperation.hh" #include "G4ParticleChange.hh" +#include "G4VBiasingOperation.hh" class GateOptnCloning : public G4VBiasingOperation { public: @@ -50,38 +49,44 @@ class GateOptnCloning : public G4VBiasingOperation { GateOptnCloning(G4String name); // -- destructor: virtual ~GateOptnCloning(); - + public: // -- Methods from G4VBiasingOperation interface: // ------------------------------------------- // -- Unsed: - virtual const G4VBiasingInteractionLaw* ProvideOccurenceBiasingInteractionLaw( const G4BiasingProcessInterface*, G4ForceCondition& ) {return 0;} - virtual G4VParticleChange* ApplyFinalStateBiasing( const G4BiasingProcessInterface*, - const G4Track*, - const G4Step*, - G4bool& ) {return 0;} + virtual const G4VBiasingInteractionLaw * + ProvideOccurenceBiasingInteractionLaw(const G4BiasingProcessInterface *, + G4ForceCondition &) { + return 0; + } + virtual G4VParticleChange * + ApplyFinalStateBiasing(const G4BiasingProcessInterface *, const G4Track *, + const G4Step *, G4bool &) { + return 0; + } // -- Used: - virtual G4double DistanceToApplyOperation( const G4Track*, - G4double, - G4ForceCondition* condition) - { - *condition = NotForced; return 0; // -- acts immediately + virtual G4double DistanceToApplyOperation(const G4Track *, G4double, + G4ForceCondition *condition) { + *condition = NotForced; + return 0; // -- acts immediately } - virtual G4VParticleChange* GenerateBiasingFinalState( const G4Track*, - const G4Step* ); - + virtual G4VParticleChange *GenerateBiasingFinalState(const G4Track *, + const G4Step *); + public: // -- Additional methods, specific to this class: // ---------------------------------------------- - void SetCloneWeights(G4double clone1Weight, G4double clone2Weight) {fClone1W = clone1Weight ; fClone2W = clone2Weight;} + void SetCloneWeights(G4double clone1Weight, G4double clone2Weight) { + fClone1W = clone1Weight; + fClone2W = clone2Weight; + } - G4Track* GetCloneTrack() const { return fCloneTrack; } + G4Track *GetCloneTrack() const { return fCloneTrack; } private: - G4double fClone1W, - fClone2W; + G4double fClone1W, fClone2W; G4ParticleChange fParticleChange; - G4Track* fCloneTrack; + G4Track *fCloneTrack; }; #endif diff --git a/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.cpp b/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.cpp index 3da3410b3..99da7b38b 100644 --- a/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.cpp +++ b/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.cpp @@ -28,155 +28,138 @@ #include "G4ILawForceFreeFlight.hh" #include "G4TransportationManager.hh" -#include "Randomize.hh" #include "G4BiasingProcessInterface.hh" +#include "Randomize.hh" GateOptnForceCommonTruncatedExp::GateOptnForceCommonTruncatedExp(G4String name) - : G4VBiasingOperation(name), - fNumberOfSharing(0), - fProcessToApply(nullptr), - fInteractionOccured(false), - fMaximumDistance(-1.0) -{ - fCommonTruncatedExpLaw = new G4ILawCommonTruncatedExp("ExpLawForOperation"+name); - fForceFreeFlightLaw = new G4ILawForceFreeFlight ("FFFLawForOperation"+name); - + : G4VBiasingOperation(name), fNumberOfSharing(0), fProcessToApply(nullptr), + fInteractionOccured(false), fMaximumDistance(-1.0) { + fCommonTruncatedExpLaw = + new G4ILawCommonTruncatedExp("ExpLawForOperation" + name); + fForceFreeFlightLaw = new G4ILawForceFreeFlight("FFFLawForOperation" + name); + fTotalCrossSection = 0.0; } -GateOptnForceCommonTruncatedExp::~GateOptnForceCommonTruncatedExp() -{ - if ( fCommonTruncatedExpLaw ) delete fCommonTruncatedExpLaw; - if ( fForceFreeFlightLaw ) delete fForceFreeFlightLaw; +GateOptnForceCommonTruncatedExp::~GateOptnForceCommonTruncatedExp() { + if (fCommonTruncatedExpLaw) + delete fCommonTruncatedExpLaw; + if (fForceFreeFlightLaw) + delete fForceFreeFlightLaw; } -const G4VBiasingInteractionLaw* GateOptnForceCommonTruncatedExp:: -ProvideOccurenceBiasingInteractionLaw( const G4BiasingProcessInterface* callingProcess, - G4ForceCondition& proposeForceCondition ) -{ - if ( callingProcess->GetWrappedProcess() == fProcessToApply ) - { - proposeForceCondition = Forced; - return fCommonTruncatedExpLaw; - } - else - { - proposeForceCondition = Forced; - return fForceFreeFlightLaw; - } +const G4VBiasingInteractionLaw * +GateOptnForceCommonTruncatedExp::ProvideOccurenceBiasingInteractionLaw( + const G4BiasingProcessInterface *callingProcess, + G4ForceCondition &proposeForceCondition) { + if (callingProcess->GetWrappedProcess() == fProcessToApply) { + proposeForceCondition = Forced; + return fCommonTruncatedExpLaw; + } else { + proposeForceCondition = Forced; + return fForceFreeFlightLaw; + } } - -G4GPILSelection GateOptnForceCommonTruncatedExp::ProposeGPILSelection( const G4GPILSelection ) -{ +G4GPILSelection +GateOptnForceCommonTruncatedExp::ProposeGPILSelection(const G4GPILSelection) { return NotCandidateForSelection; } +G4VParticleChange *GateOptnForceCommonTruncatedExp::ApplyFinalStateBiasing( + const G4BiasingProcessInterface *callingProcess, const G4Track *track, + const G4Step *step, G4bool &forceFinalState) { + if (callingProcess->GetWrappedProcess() != fProcessToApply) { + forceFinalState = true; + fDummyParticleChange.Initialize(*track); + return &fDummyParticleChange; + } + if (fInteractionOccured) { + forceFinalState = true; + fDummyParticleChange.Initialize(*track); + return &fDummyParticleChange; + } -G4VParticleChange* GateOptnForceCommonTruncatedExp::ApplyFinalStateBiasing( const G4BiasingProcessInterface* callingProcess, - const G4Track* track, - const G4Step* step, - G4bool& forceFinalState ) -{ - if ( callingProcess->GetWrappedProcess() != fProcessToApply ) - { - forceFinalState = true; - fDummyParticleChange.Initialize( *track ); - return &fDummyParticleChange; - } - if ( fInteractionOccured ) - { - forceFinalState = true; - fDummyParticleChange.Initialize( *track ); - return &fDummyParticleChange; - } - // -- checks if process won the GPIL race: - G4double processGPIL = callingProcess->GetPostStepGPIL() < callingProcess->GetAlongStepGPIL() ? - callingProcess->GetPostStepGPIL() : callingProcess->GetAlongStepGPIL() ; - if ( processGPIL <= step->GetStepLength() ) - { - // -- if process won, wrapped process produces the final state. - // -- In this case, the weight for occurrence biasing is applied - // -- by the callingProcess, at exit of present method. This is - // -- selected by "forceFinalState = false": - forceFinalState = false; - fInteractionOccured = true; - return callingProcess->GetWrappedProcess()->PostStepDoIt( *track, *step ); - } - else - { - forceFinalState = true; - fDummyParticleChange.Initialize( *track ); - return &fDummyParticleChange; - } + G4double processGPIL = + callingProcess->GetPostStepGPIL() < callingProcess->GetAlongStepGPIL() + ? callingProcess->GetPostStepGPIL() + : callingProcess->GetAlongStepGPIL(); + if (processGPIL <= step->GetStepLength()) { + // -- if process won, wrapped process produces the final state. + // -- In this case, the weight for occurrence biasing is applied + // -- by the callingProcess, at exit of present method. This is + // -- selected by "forceFinalState = false": + forceFinalState = false; + fInteractionOccured = true; + return callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); + } else { + forceFinalState = true; + fDummyParticleChange.Initialize(*track); + return &fDummyParticleChange; + } } - -void GateOptnForceCommonTruncatedExp::AddCrossSection( const G4VProcess* process, G4double crossSection ) -{ - fTotalCrossSection += crossSection; - fCrossSections[process] = crossSection; - fNumberOfSharing = fCrossSections.size(); +void GateOptnForceCommonTruncatedExp::AddCrossSection(const G4VProcess *process, + G4double crossSection) { + fTotalCrossSection += crossSection; + fCrossSections[process] = crossSection; + fNumberOfSharing = fCrossSections.size(); } - -void GateOptnForceCommonTruncatedExp::Initialize( const G4Track* track ) -{ +void GateOptnForceCommonTruncatedExp::Initialize(const G4Track *track) { fCrossSections.clear(); - fTotalCrossSection = 0.0; - fNumberOfSharing = 0; - fProcessToApply = 0; + fTotalCrossSection = 0.0; + fNumberOfSharing = 0; + fProcessToApply = 0; fInteractionOccured = false; - fInitialMomentum = track->GetMomentum(); - - G4VSolid* currentSolid = track->GetVolume()->GetLogicalVolume()->GetSolid(); - G4ThreeVector localPosition = (G4TransportationManager::GetTransportationManager()-> - GetNavigatorForTracking()-> - GetGlobalToLocalTransform()).TransformPoint(track->GetPosition()); - G4ThreeVector localDirection = (G4TransportationManager::GetTransportationManager()-> - GetNavigatorForTracking()-> - GetGlobalToLocalTransform()).TransformAxis(track->GetMomentumDirection()); + fInitialMomentum = track->GetMomentum(); + + G4VSolid *currentSolid = track->GetVolume()->GetLogicalVolume()->GetSolid(); + G4ThreeVector localPosition = + (G4TransportationManager::GetTransportationManager() + ->GetNavigatorForTracking() + ->GetGlobalToLocalTransform()) + .TransformPoint(track->GetPosition()); + G4ThreeVector localDirection = + (G4TransportationManager::GetTransportationManager() + ->GetNavigatorForTracking() + ->GetGlobalToLocalTransform()) + .TransformAxis(track->GetMomentumDirection()); fMaximumDistance = currentSolid->DistanceToOut(localPosition, localDirection); - if ( fMaximumDistance <= DBL_MIN ) fMaximumDistance = 0.0; - fCommonTruncatedExpLaw->SetMaximumDistance( fMaximumDistance ); + if (fMaximumDistance <= DBL_MIN) + fMaximumDistance = 0.0; + fCommonTruncatedExpLaw->SetMaximumDistance(fMaximumDistance); } - -void GateOptnForceCommonTruncatedExp::UpdateForStep( const G4Step* step ) -{ +void GateOptnForceCommonTruncatedExp::UpdateForStep(const G4Step *step) { fCrossSections.clear(); - fTotalCrossSection = 0.0; - fNumberOfSharing = 0; - fProcessToApply = 0; - - fCommonTruncatedExpLaw->UpdateForStep( step->GetStepLength() ); + fTotalCrossSection = 0.0; + fNumberOfSharing = 0; + fProcessToApply = 0; + + fCommonTruncatedExpLaw->UpdateForStep(step->GetStepLength()); fMaximumDistance = fCommonTruncatedExpLaw->GetMaximumDistance(); } - -void GateOptnForceCommonTruncatedExp::Sample() -{ - fCommonTruncatedExpLaw->SetForceCrossSection( fTotalCrossSection ); +void GateOptnForceCommonTruncatedExp::Sample() { + fCommonTruncatedExpLaw->SetForceCrossSection(fTotalCrossSection); fCommonTruncatedExpLaw->Sample(); ChooseProcessToApply(); - fCommonTruncatedExpLaw->SetSelectedProcessXSfraction(fCrossSections[fProcessToApply] / fTotalCrossSection); + fCommonTruncatedExpLaw->SetSelectedProcessXSfraction( + fCrossSections[fProcessToApply] / fTotalCrossSection); } - -void GateOptnForceCommonTruncatedExp::ChooseProcessToApply() -{ - G4double sigmaRand = G4UniformRand() * fTotalCrossSection; +void GateOptnForceCommonTruncatedExp::ChooseProcessToApply() { + G4double sigmaRand = G4UniformRand() * fTotalCrossSection; G4double sigmaSelect = 0.0; - for ( std::map< const G4VProcess*, G4double>::const_iterator it = fCrossSections.begin(); - it != fCrossSections.end(); - it++) - { - sigmaSelect += (*it).second; - if ( sigmaRand <= sigmaSelect ) - { - fProcessToApply = (*it).first; - break; - } + for (std::map::const_iterator it = + fCrossSections.begin(); + it != fCrossSections.end(); it++) { + sigmaSelect += (*it).second; + if (sigmaRand <= sigmaSelect) { + fProcessToApply = (*it).first; + break; } + } } diff --git a/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.h b/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.h index ec15f1745..491a6f337 100644 --- a/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.h +++ b/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.h @@ -45,13 +45,12 @@ //--------------------------------------------------------------- // Initial version Nov. 2013 M. Verderi - #ifndef GateOptnForceCommonTruncatedExp_h #define GateOptnForceCommonTruncatedExp_h 1 -#include "G4VBiasingOperation.hh" -#include "G4ThreeVector.hh" #include "G4ParticleChangeForNothing.hh" +#include "G4ThreeVector.hh" +#include "G4VBiasingOperation.hh" class G4ILawCommonTruncatedExp; class G4ILawForceFreeFlight; #include @@ -62,61 +61,64 @@ class GateOptnForceCommonTruncatedExp : public G4VBiasingOperation { GateOptnForceCommonTruncatedExp(G4String name); // -- destructor: virtual ~GateOptnForceCommonTruncatedExp(); - + public: // -- Methods from G4VBiasingOperation interface: // ------------------------------------------- // -- Used: - virtual const G4VBiasingInteractionLaw* ProvideOccurenceBiasingInteractionLaw( const G4BiasingProcessInterface*, G4ForceCondition& ); - virtual G4double ProposeAlongStepLimit( const G4BiasingProcessInterface* ) { return DBL_MAX; } - virtual G4GPILSelection ProposeGPILSelection( const G4GPILSelection processSelection ); - virtual G4VParticleChange* ApplyFinalStateBiasing( const G4BiasingProcessInterface*, - const G4Track*, - const G4Step*, - G4bool& ); + virtual const G4VBiasingInteractionLaw * + ProvideOccurenceBiasingInteractionLaw(const G4BiasingProcessInterface *, + G4ForceCondition &); + virtual G4double ProposeAlongStepLimit(const G4BiasingProcessInterface *) { + return DBL_MAX; + } + virtual G4GPILSelection + ProposeGPILSelection(const G4GPILSelection processSelection); + virtual G4VParticleChange * + ApplyFinalStateBiasing(const G4BiasingProcessInterface *, const G4Track *, + const G4Step *, G4bool &); // -- Unused: - virtual G4double DistanceToApplyOperation( const G4Track*, - G4double, - G4ForceCondition*) {return DBL_MAX;} - virtual G4VParticleChange* GenerateBiasingFinalState( const G4Track*, - const G4Step* ) {return 0;} - + virtual G4double DistanceToApplyOperation(const G4Track *, G4double, + G4ForceCondition *) { + return DBL_MAX; + } + virtual G4VParticleChange *GenerateBiasingFinalState(const G4Track *, + const G4Step *) { + return 0; + } + public: // -- Additional methods, specific to this class: // ---------------------------------------------- // -- return concrete type of interaction laws: - G4ILawCommonTruncatedExp* GetCommonTruncatedExpLaw() - { + G4ILawCommonTruncatedExp *GetCommonTruncatedExpLaw() { return fCommonTruncatedExpLaw; } - G4ILawForceFreeFlight* GetForceFreeFlightLaw() - { - return fForceFreeFlightLaw; - } - - void Initialize( const G4Track* ); - void UpdateForStep( const G4Step* ); - void Sample(); - const G4ThreeVector& GetInitialMomentum() const { return fInitialMomentum; } - G4double GetMaximumDistance() const { return fMaximumDistance; } - void ChooseProcessToApply(); - const G4VProcess* GetProcessToApply() const { return fProcessToApply; } - void AddCrossSection( const G4VProcess*, G4double ); - size_t GetNumberOfSharing() const { return fNumberOfSharing;} - void SetInteractionOccured( G4bool b ) { fInteractionOccured = b; } - G4bool GetInteractionOccured() const { return fInteractionOccured; } - + G4ILawForceFreeFlight *GetForceFreeFlightLaw() { return fForceFreeFlightLaw; } + + void Initialize(const G4Track *); + void UpdateForStep(const G4Step *); + void Sample(); + const G4ThreeVector &GetInitialMomentum() const { return fInitialMomentum; } + G4double GetMaximumDistance() const { return fMaximumDistance; } + void ChooseProcessToApply(); + const G4VProcess *GetProcessToApply() const { return fProcessToApply; } + void AddCrossSection(const G4VProcess *, G4double); + size_t GetNumberOfSharing() const { return fNumberOfSharing; } + void SetInteractionOccured(G4bool b) { fInteractionOccured = b; } + G4bool GetInteractionOccured() const { return fInteractionOccured; } + private: - G4ILawCommonTruncatedExp* fCommonTruncatedExpLaw; - G4ILawForceFreeFlight* fForceFreeFlightLaw; - G4double fTotalCrossSection; - std::map < const G4VProcess*, G4double > fCrossSections; - size_t fNumberOfSharing; - const G4VProcess* fProcessToApply; - G4bool fInteractionOccured; - G4ThreeVector fInitialMomentum; - G4double fMaximumDistance; - G4ParticleChangeForNothing fDummyParticleChange; + G4ILawCommonTruncatedExp *fCommonTruncatedExpLaw; + G4ILawForceFreeFlight *fForceFreeFlightLaw; + G4double fTotalCrossSection; + std::map fCrossSections; + size_t fNumberOfSharing; + const G4VProcess *fProcessToApply; + G4bool fInteractionOccured; + G4ThreeVector fInitialMomentum; + G4double fMaximumDistance; + G4ParticleChangeForNothing fDummyParticleChange; }; #endif diff --git a/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.cpp b/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.cpp index 686501dd6..f9c07d2d1 100644 --- a/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.cpp +++ b/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.cpp @@ -24,81 +24,75 @@ // ******************************************************************** // #include "GateOptnForceFreeFlight.h" -#include "G4ILawForceFreeFlight.hh" #include "G4BiasingProcessInterface.hh" +#include "G4ILawForceFreeFlight.hh" #include "G4Step.hh" - - GateOptnForceFreeFlight::GateOptnForceFreeFlight(G4String name) - : G4VBiasingOperation ( name ), - fCumulatedWeightChange ( -1.0 ), - fInitialTrackWeight ( -1.0 ), - fOperationComplete ( true ) -{ - fForceFreeFlightInteractionLaw = new G4ILawForceFreeFlight("LawForOperation"+name); + : G4VBiasingOperation(name), fCumulatedWeightChange(-1.0), + fInitialTrackWeight(-1.0), fOperationComplete(true) { + fForceFreeFlightInteractionLaw = + new G4ILawForceFreeFlight("LawForOperation" + name); } -GateOptnForceFreeFlight::~GateOptnForceFreeFlight() -{ - if ( fForceFreeFlightInteractionLaw ) delete fForceFreeFlightInteractionLaw; +GateOptnForceFreeFlight::~GateOptnForceFreeFlight() { + if (fForceFreeFlightInteractionLaw) + delete fForceFreeFlightInteractionLaw; } -const G4VBiasingInteractionLaw* GateOptnForceFreeFlight::ProvideOccurenceBiasingInteractionLaw( const G4BiasingProcessInterface*, G4ForceCondition& proposeForceCondition ) -{ +const G4VBiasingInteractionLaw * +GateOptnForceFreeFlight::ProvideOccurenceBiasingInteractionLaw( + const G4BiasingProcessInterface *, + G4ForceCondition &proposeForceCondition) { fOperationComplete = false; proposeForceCondition = Forced; return fForceFreeFlightInteractionLaw; } - -G4VParticleChange* GateOptnForceFreeFlight::ApplyFinalStateBiasing( const G4BiasingProcessInterface* callingProcess, - const G4Track* track, - const G4Step* step, - G4bool& forceFinalState) -{ - // -- If the track is reaching the volume boundary, its free flight ends. In this case, its zero - // -- weight is brought back to non-zero value: its initial weight is restored by the first - // -- ApplyFinalStateBiasing operation called, and the weight for force free flight is applied +G4VParticleChange *GateOptnForceFreeFlight::ApplyFinalStateBiasing( + const G4BiasingProcessInterface *callingProcess, const G4Track *track, + const G4Step *step, G4bool &forceFinalState) { + // -- If the track is reaching the volume boundary, its free flight ends. In + // this case, its zero + // -- weight is brought back to non-zero value: its initial weight is restored + // by the first + // -- ApplyFinalStateBiasing operation called, and the weight for force free + // flight is applied // -- is applied by each operation. - // -- If the track is not reaching the volume boundary, it zero weight flight continues. + // -- If the track is not reaching the volume boundary, it zero weight flight + // continues. - fParticleChange.Initialize( *track ); - forceFinalState = true; - if ( step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary ) - { - // -- Sanity checks: - if ( fInitialTrackWeight <= DBL_MIN ) - { - G4ExceptionDescription ed; - ed << " Initial track weight is null ! " << G4endl; - G4Exception(" GateOptnForceFreeFlight::ApplyFinalStateBiasing(...)", - "BIAS.GEN.05", - JustWarning, - ed); - } - if ( fCumulatedWeightChange <= DBL_MIN ) - { - G4ExceptionDescription ed; - ed << " Cumulated weight is null ! " << G4endl; - G4Exception(" GateOptnForceFreeFlight::ApplyFinalStateBiasing(...)", - "BIAS.GEN.06", - JustWarning, - ed); - } - - G4double proposedWeight = track->GetWeight(); - if ( callingProcess->GetIsFirstPostStepDoItInterface() ) proposedWeight = fCumulatedWeightChange * fInitialTrackWeight; - else proposedWeight *= fCumulatedWeightChange; - fParticleChange.ProposeWeight(proposedWeight); - fOperationComplete = true; + fParticleChange.Initialize(*track); + forceFinalState = true; + if (step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary) { + // -- Sanity checks: + if (fInitialTrackWeight <= DBL_MIN) { + G4ExceptionDescription ed; + ed << " Initial track weight is null ! " << G4endl; + G4Exception(" GateOptnForceFreeFlight::ApplyFinalStateBiasing(...)", + "BIAS.GEN.05", JustWarning, ed); + } + if (fCumulatedWeightChange <= DBL_MIN) { + G4ExceptionDescription ed; + ed << " Cumulated weight is null ! " << G4endl; + G4Exception(" GateOptnForceFreeFlight::ApplyFinalStateBiasing(...)", + "BIAS.GEN.06", JustWarning, ed); } - + + G4double proposedWeight = track->GetWeight(); + if (callingProcess->GetIsFirstPostStepDoItInterface()) + proposedWeight = fCumulatedWeightChange * fInitialTrackWeight; + else + proposedWeight *= fCumulatedWeightChange; + fParticleChange.ProposeWeight(proposedWeight); + fOperationComplete = true; + } + return &fParticleChange; } - -void GateOptnForceFreeFlight::AlongMoveBy( const G4BiasingProcessInterface*, const G4Step*, G4double weightChange ) -{ +void GateOptnForceFreeFlight::AlongMoveBy(const G4BiasingProcessInterface *, + const G4Step *, + G4double weightChange) { fCumulatedWeightChange *= weightChange; } diff --git a/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.h b/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.h index e7dffb49a..ab8e0641d 100644 --- a/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.h +++ b/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.h @@ -46,52 +46,60 @@ #ifndef GateOptnForceFreeFlight_h #define GateOptnForceFreeFlight_h 1 -#include "G4VBiasingOperation.hh" #include "G4ForceCondition.hh" #include "G4ParticleChange.hh" // -- §§ should add a dedicated "weight change only" particle change +#include "G4VBiasingOperation.hh" class G4ILawForceFreeFlight; - class GateOptnForceFreeFlight : public G4VBiasingOperation { public: // -- Constructor : GateOptnForceFreeFlight(G4String name); // -- destructor: virtual ~GateOptnForceFreeFlight(); - + public: // -- Methods from G4VBiasingOperation interface: // ------------------------------------------- // -- Used: - virtual const G4VBiasingInteractionLaw* ProvideOccurenceBiasingInteractionLaw( const G4BiasingProcessInterface*, G4ForceCondition& ); - virtual void AlongMoveBy( const G4BiasingProcessInterface*, const G4Step*, G4double ); - virtual G4VParticleChange* ApplyFinalStateBiasing( const G4BiasingProcessInterface*, const G4Track*, const G4Step*, G4bool&); + virtual const G4VBiasingInteractionLaw * + ProvideOccurenceBiasingInteractionLaw(const G4BiasingProcessInterface *, + G4ForceCondition &); + virtual void AlongMoveBy(const G4BiasingProcessInterface *, const G4Step *, + G4double); + virtual G4VParticleChange * + ApplyFinalStateBiasing(const G4BiasingProcessInterface *, const G4Track *, + const G4Step *, G4bool &); // -- Unused: - virtual G4double DistanceToApplyOperation( const G4Track*, - G4double, - G4ForceCondition*) {return DBL_MAX;} - virtual G4VParticleChange* GenerateBiasingFinalState( const G4Track*, - const G4Step* ) {return 0;} - + virtual G4double DistanceToApplyOperation(const G4Track *, G4double, + G4ForceCondition *) { + return DBL_MAX; + } + virtual G4VParticleChange *GenerateBiasingFinalState(const G4Track *, + const G4Step *) { + return 0; + } public: // -- Additional methods, specific to this class: // ---------------------------------------------- // -- return concrete type of interaction law: - G4ILawForceFreeFlight* GetForceFreeFlightLaw() { + G4ILawForceFreeFlight *GetForceFreeFlightLaw() { return fForceFreeFlightInteractionLaw; } // -- initialization for weight: - void ResetInitialTrackWeight(G4double w) {fInitialTrackWeight = w; fCumulatedWeightChange = 1.0;} + void ResetInitialTrackWeight(G4double w) { + fInitialTrackWeight = w; + fCumulatedWeightChange = 1.0; + } G4bool OperationComplete() const { return fOperationComplete; } - + private: - G4ILawForceFreeFlight* fForceFreeFlightInteractionLaw; - G4double fCumulatedWeightChange, - fInitialTrackWeight; - G4ParticleChange fParticleChange; - G4bool fOperationComplete; + G4ILawForceFreeFlight *fForceFreeFlightInteractionLaw; + G4double fCumulatedWeightChange, fInitialTrackWeight; + G4ParticleChange fParticleChange; + G4bool fOperationComplete; }; #endif diff --git a/core/opengate_core/opengate_lib/GateOptrForceCollision.cpp b/core/opengate_core/opengate_lib/GateOptrForceCollision.cpp index cc5294f8e..c0b94b60b 100644 --- a/core/opengate_core/opengate_lib/GateOptrForceCollision.cpp +++ b/core/opengate_core/opengate_lib/GateOptrForceCollision.cpp @@ -24,14 +24,14 @@ // ******************************************************************** // #include "GateOptrForceCollision.h" -#include "GateOptrForceCollisionTrackData.h" #include "G4BiasingProcessInterface.hh" #include "G4PhysicsModelCatalog.hh" +#include "GateOptrForceCollisionTrackData.h" -#include "GateOptnForceCommonTruncatedExp.h" #include "G4ILawCommonTruncatedExp.hh" -#include "GateOptnForceFreeFlight.h" #include "GateOptnCloning.h" +#include "GateOptnForceCommonTruncatedExp.h" +#include "GateOptnForceFreeFlight.h" #include "G4Step.hh" #include "G4StepPoint.hh" @@ -47,18 +47,16 @@ // -- §§ consider calling other constructor, thanks to C++11 GateOptrForceCollision::GateOptrForceCollision(py::dict &user_info) - : G4VBiasingOperator("forceCollisionActor"),GateVActor(user_info, false), - fForceCollisionModelID(G4PhysicsModelCatalog::GetModelID("model_GenBiasForceCollision")), - fCurrentTrack(nullptr), - fCurrentTrackData(nullptr), - fInitialTrackWeight(-1.0), - fSetup(true) -{ - fSharedForceInteractionOperation = new GateOptnForceCommonTruncatedExp("SharedForceInteraction"); - fCloningOperation = new GateOptnCloning("Cloning"); + : G4VBiasingOperator("forceCollisionActor"), GateVActor(user_info, false), + fForceCollisionModelID( + G4PhysicsModelCatalog::GetModelID("model_GenBiasForceCollision")), + fCurrentTrack(nullptr), fCurrentTrackData(nullptr), + fInitialTrackWeight(-1.0), fSetup(true) { + fSharedForceInteractionOperation = + new GateOptnForceCommonTruncatedExp("SharedForceInteraction"); + fCloningOperation = new GateOptnCloning("Cloning"); } - void GateOptrForceCollision::AttachAllLogicalDaughtersVolumes( G4LogicalVolume *volume) { AttachTo(volume); @@ -72,11 +70,12 @@ void GateOptrForceCollision::AttachAllLogicalDaughtersVolumes( } } -GateOptrForceCollision::~GateOptrForceCollision() -{ - for ( std::map< const G4BiasingProcessInterface*, GateOptnForceFreeFlight* >::iterator it = fFreeFlightOperations.begin() ; - it != fFreeFlightOperations.end() ; - it++ ) delete (*it).second; +GateOptrForceCollision::~GateOptrForceCollision() { + for (std::map::iterator it = + fFreeFlightOperations.begin(); + it != fFreeFlightOperations.end(); it++) + delete (*it).second; delete fSharedForceInteractionOperation; delete fCloningOperation; } @@ -94,41 +93,45 @@ void GateOptrForceCollision::InitializeUserInfo(py::dict &user_info) { } } - - -void GateOptrForceCollision::Configure() -{ +void GateOptrForceCollision::Configure() { // -- build free flight operations: ConfigureForWorker(); } - -void GateOptrForceCollision::ConfigureForWorker() -{ - // -- start by remembering processes under biasing, and create needed biasing operations: - if ( fSetup ) +void GateOptrForceCollision::ConfigureForWorker() { + // -- start by remembering processes under biasing, and create needed biasing + // operations: + if (fSetup) { + const G4ProcessManager *processManager = + fParticleToBias->GetProcessManager(); + const G4BiasingProcessSharedData *interfaceProcessSharedData = + G4BiasingProcessInterface::GetSharedData(processManager); + if (interfaceProcessSharedData) // -- sharedData tested, as is can happen a + // user attaches an operator + // -- to a volume but without defining + // BiasingProcessInterface processes. { - const G4ProcessManager* processManager = fParticleToBias->GetProcessManager(); - const G4BiasingProcessSharedData* interfaceProcessSharedData = G4BiasingProcessInterface::GetSharedData( processManager ); - if ( interfaceProcessSharedData ) // -- sharedData tested, as is can happen a user attaches an operator - // -- to a volume but without defining BiasingProcessInterface processes. - { - for ( size_t i = 0 ; i < (interfaceProcessSharedData->GetPhysicsBiasingProcessInterfaces()).size(); i++ ) - { - const G4BiasingProcessInterface* wrapperProcess = - (interfaceProcessSharedData->GetPhysicsBiasingProcessInterfaces())[i]; - G4String operationName = "FreeFlight-"+wrapperProcess->GetWrappedProcess()->GetProcessName(); - fFreeFlightOperations[wrapperProcess] = new GateOptnForceFreeFlight(operationName); - } - } - fSetup = false; + for (size_t i = 0; + i < + (interfaceProcessSharedData->GetPhysicsBiasingProcessInterfaces()) + .size(); + i++) { + const G4BiasingProcessInterface *wrapperProcess = + (interfaceProcessSharedData + ->GetPhysicsBiasingProcessInterfaces())[i]; + G4String operationName = + "FreeFlight-" + + wrapperProcess->GetWrappedProcess()->GetProcessName(); + fFreeFlightOperations[wrapperProcess] = + new GateOptnForceFreeFlight(operationName); + } } + fSetup = false; + } } - -void GateOptrForceCollision::StartRun() -{ - G4LogicalVolume *biasingVolume = +void GateOptrForceCollision::StartRun() { + G4LogicalVolume *biasingVolume = G4LogicalVolumeStore::GetInstance()->GetVolume(fAttachedToVolumeName); // Here we need to attach all the daughters and daughters of daughters (...) @@ -137,21 +140,23 @@ void GateOptrForceCollision::StartRun() AttachAllLogicalDaughtersVolumes(biasingVolume); } - -G4VBiasingOperation* GateOptrForceCollision::ProposeOccurenceBiasingOperation(const G4Track* track, const G4BiasingProcessInterface* callingProcess) -{ +G4VBiasingOperation *GateOptrForceCollision::ProposeOccurenceBiasingOperation( + const G4Track *track, const G4BiasingProcessInterface *callingProcess) { // -- does nothing if particle is not of requested type: - if ( track->GetDefinition() != fParticleToBias ) return 0; + if (track->GetDefinition() != fParticleToBias) + return 0; // -- trying to get auxiliary track data... - if ( fCurrentTrackData == nullptr ) - { - // ... and if the track has no aux. track data, it means its biasing is not started yet (note that cloning is to happen first): - fCurrentTrackData = (GateOptrForceCollisionTrackData*)(track->GetAuxiliaryTrackInformation(fForceCollisionModelID)); - if ( fCurrentTrackData == nullptr ) return nullptr; - } + if (fCurrentTrackData == nullptr) { + // ... and if the track has no aux. track data, it means its biasing is not + // started yet (note that cloning is to happen first): + fCurrentTrackData = + (GateOptrForceCollisionTrackData *)(track->GetAuxiliaryTrackInformation( + fForceCollisionModelID)); + if (fCurrentTrackData == nullptr) + return nullptr; + } - // -- Send force free flight to the callingProcess: // ------------------------------------------------ // -- The track has been cloned in the previous step, it has now to be @@ -162,281 +167,265 @@ G4VBiasingOperation* GateOptrForceCollision::ProposeOccurenceBiasingOperation(co // -- being its initial weight * the weight for the free flight travel, // -- this last one being per process. The initial weight is common, and is // -- arbitrary asked to the first operation to take care of it. - if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeFreeFlight ) - { - GateOptnForceFreeFlight* operation = fFreeFlightOperations[callingProcess]; - if ( callingProcess->GetWrappedProcess()->GetCurrentInteractionLength() < DBL_MAX/10. ) - { - // -- the initial track weight will be restored only by the first DoIt free flight: - operation->ResetInitialTrackWeight(fInitialTrackWeight); - return operation; - } - else - { - return nullptr; - } + if (fCurrentTrackData->fForceCollisionState == + ForceCollisionState::toBeFreeFlight) { + GateOptnForceFreeFlight *operation = fFreeFlightOperations[callingProcess]; + if (callingProcess->GetWrappedProcess()->GetCurrentInteractionLength() < + DBL_MAX / 10.) { + // -- the initial track weight will be restored only by the first DoIt + // free flight: + operation->ResetInitialTrackWeight(fInitialTrackWeight); + return operation; + } else { + return nullptr; } - + } // -- Send force interaction operation to the callingProcess: // ---------------------------------------------------------- // -- at this level, a copy of the track entering the volume was // -- generated (borned) earlier. This copy will make the forced // -- interaction in the volume. - if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced ) - { - // -- Remember if this calling process is the first of the physics wrapper in - // -- the PostStepGPIL loop (using default argument of method below): - G4bool isFirstPhysGPIL = callingProcess-> GetIsFirstPostStepGPILInterface(); - - // -- [*first process*] Initialize or update force interaction operation: - if ( isFirstPhysGPIL ) - { - // -- first step of cloned track, initialize the forced interaction operation: - if ( track->GetCurrentStepNumber() == 1 ) fSharedForceInteractionOperation->Initialize( track ); - else - { - if ( fSharedForceInteractionOperation->GetInitialMomentum() != track->GetMomentum() ) - { - // -- means that some other physics process, not under control of the forced interaction operation, - // -- has occured, need to re-initialize the operation as distance to boundary has changed. - // -- [ Note the re-initialization is only possible for a Markovian law. ] - fSharedForceInteractionOperation->Initialize( track ); - } - else - { - // -- means that some other non-physics process (biasing or not, like step limit), has occured, - // -- but track conserves its momentum direction, only need to reduced the maximum distance for - // -- forced interaction. - // -- [ Note the update is only possible for a Markovian law. ] - fSharedForceInteractionOperation->UpdateForStep( track->GetStep() ); - } - } - } - - // -- [*all processes*] Sanity check : it may happen in limit cases that distance to - // -- out is zero, weight would be infinite in this case: abort forced interaction - // -- and abandon biasing. - if ( fSharedForceInteractionOperation->GetMaximumDistance() < DBL_MIN ) - { - fCurrentTrackData->Reset(); - return 0; - } - - // -- [* first process*] collect cross-sections and sample force law to determine interaction length - // -- and winning process: - if ( isFirstPhysGPIL ) - { - // -- collect cross-sections: - // -- ( Remember that the first of the G4BiasingProcessInterface triggers the update - // -- of these cross-sections ) - const G4BiasingProcessSharedData* sharedData = callingProcess->GetSharedData(); - for ( size_t i = 0 ; i < (sharedData->GetPhysicsBiasingProcessInterfaces()).size() ; i++ ) - { - const G4BiasingProcessInterface* wrapperProcess = ( sharedData->GetPhysicsBiasingProcessInterfaces() )[i]; - G4double interactionLength = wrapperProcess->GetWrappedProcess()->GetCurrentInteractionLength(); - // -- keep only well defined cross-sections, other processes are ignored. These are not pathological cases - // -- but cases where a threhold effect par example (eg pair creation) may be at play. (**!**) - if ( interactionLength < DBL_MAX/10. ) - fSharedForceInteractionOperation->AddCrossSection( wrapperProcess->GetWrappedProcess(), 1.0/interactionLength ); - } - // -- sample the shared law (interaction length, and winning process): - if ( fSharedForceInteractionOperation->GetNumberOfSharing() > 0 ) fSharedForceInteractionOperation->Sample(); - } - - // -- [*all processes*] Send operation for processes with well defined XS (see "**!**" above): - G4VBiasingOperation* operationToReturn = nullptr; - if ( callingProcess->GetWrappedProcess()->GetCurrentInteractionLength() < DBL_MAX/10. ) operationToReturn = fSharedForceInteractionOperation; - return operationToReturn; - - - } // -- end of "if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced )" - - + if (fCurrentTrackData->fForceCollisionState == + ForceCollisionState::toBeForced) { + // -- Remember if this calling process is the first of the physics wrapper + // in + // -- the PostStepGPIL loop (using default argument of method below): + G4bool isFirstPhysGPIL = callingProcess->GetIsFirstPostStepGPILInterface(); + + // -- [*first process*] Initialize or update force interaction operation: + if (isFirstPhysGPIL) { + // -- first step of cloned track, initialize the forced interaction + // operation: + if (track->GetCurrentStepNumber() == 1) + fSharedForceInteractionOperation->Initialize(track); + else { + if (fSharedForceInteractionOperation->GetInitialMomentum() != + track->GetMomentum()) { + // -- means that some other physics process, not under control of the + // forced interaction operation, + // -- has occured, need to re-initialize the operation as distance to + // boundary has changed. + // -- [ Note the re-initialization is only possible for a Markovian + // law. ] + fSharedForceInteractionOperation->Initialize(track); + } else { + // -- means that some other non-physics process (biasing or not, like + // step limit), has occured, + // -- but track conserves its momentum direction, only need to reduced + // the maximum distance for + // -- forced interaction. + // -- [ Note the update is only possible for a Markovian law. ] + fSharedForceInteractionOperation->UpdateForStep(track->GetStep()); + } + } + } + + // -- [*all processes*] Sanity check : it may happen in limit cases that + // distance to + // -- out is zero, weight would be infinite in this case: abort forced + // interaction + // -- and abandon biasing. + if (fSharedForceInteractionOperation->GetMaximumDistance() < DBL_MIN) { + fCurrentTrackData->Reset(); + return 0; + } + + // -- [* first process*] collect cross-sections and sample force law to + // determine interaction length + // -- and winning process: + if (isFirstPhysGPIL) { + // -- collect cross-sections: + // -- ( Remember that the first of the G4BiasingProcessInterface triggers + // the update + // -- of these cross-sections ) + const G4BiasingProcessSharedData *sharedData = + callingProcess->GetSharedData(); + for (size_t i = 0; + i < (sharedData->GetPhysicsBiasingProcessInterfaces()).size(); i++) { + const G4BiasingProcessInterface *wrapperProcess = + (sharedData->GetPhysicsBiasingProcessInterfaces())[i]; + G4double interactionLength = + wrapperProcess->GetWrappedProcess()->GetCurrentInteractionLength(); + // -- keep only well defined cross-sections, other processes are + // ignored. These are not pathological cases + // -- but cases where a threhold effect par example (eg pair creation) + // may be at play. (**!**) + if (interactionLength < DBL_MAX / 10.) + fSharedForceInteractionOperation->AddCrossSection( + wrapperProcess->GetWrappedProcess(), 1.0 / interactionLength); + } + // -- sample the shared law (interaction length, and winning process): + if (fSharedForceInteractionOperation->GetNumberOfSharing() > 0) + fSharedForceInteractionOperation->Sample(); + } + + // -- [*all processes*] Send operation for processes with well defined XS + // (see "**!**" above): + G4VBiasingOperation *operationToReturn = nullptr; + if (callingProcess->GetWrappedProcess()->GetCurrentInteractionLength() < + DBL_MAX / 10.) + operationToReturn = fSharedForceInteractionOperation; + return operationToReturn; + + } // -- end of "if ( fCurrentTrackData->fForceCollisionState == + // ForceCollisionState::toBeForced )" + // -- other cases here: particle appearing in the volume by some // -- previous interaction : we decide to not bias these. return 0; - } +G4VBiasingOperation *GateOptrForceCollision::ProposeNonPhysicsBiasingOperation( + const G4Track *track, + const G4BiasingProcessInterface * /* callingProcess */) { + if (track->GetDefinition() != fParticleToBias) + return nullptr; -G4VBiasingOperation* GateOptrForceCollision::ProposeNonPhysicsBiasingOperation(const G4Track* track, - const G4BiasingProcessInterface* /* callingProcess */) -{ - if ( track->GetDefinition() != fParticleToBias ) return nullptr; - // -- Apply biasing scheme only to tracks entering the volume. // -- A "cloning" is done: - // -- - the primary will be forced flight under a zero weight up to volume exit, + // -- - the primary will be forced flight under a zero weight up to volume + // exit, // -- where the weight will be restored with proper weight for free flight // -- - the clone will be forced to interact in the volume. - if ( track->GetStep()->GetPreStepPoint()->GetStepStatus() == fGeomBoundary ) // -- §§§ extent to case of a track shoot on the boundary ? - { - // -- check that track is free of undergoing biasing scheme ( no biasing data, or no active active ) - // -- Get possible track data: - fCurrentTrackData = (GateOptrForceCollisionTrackData*)(track->GetAuxiliaryTrackInformation(fForceCollisionModelID)); - if ( fCurrentTrackData != nullptr ) - { - if ( fCurrentTrackData->IsFreeFromBiasing() ) - { - // -- takes "ownership" (some track data created before, left free, reuse it): - fCurrentTrackData->fForceCollisionOperator = this ; - } - else - { - // §§§ Would something be really wrong in this case ? Could this be that a process made a zero step ? - } - } - else - { - fCurrentTrackData = new GateOptrForceCollisionTrackData( this ); - track->SetAuxiliaryTrackInformation(fForceCollisionModelID, fCurrentTrackData); - } - fCurrentTrackData->fForceCollisionState = ForceCollisionState::toBeCloned; - fInitialTrackWeight = track->GetWeight(); - fCloningOperation->SetCloneWeights(0.0, fInitialTrackWeight); - return fCloningOperation; + if (track->GetStep()->GetPreStepPoint()->GetStepStatus() == + fGeomBoundary) // -- §§§ extent to case of a track shoot on the boundary ? + { + // -- check that track is free of undergoing biasing scheme ( no biasing + // data, or no active active ) + // -- Get possible track data: + fCurrentTrackData = + (GateOptrForceCollisionTrackData *)(track->GetAuxiliaryTrackInformation( + fForceCollisionModelID)); + if (fCurrentTrackData != nullptr) { + if (fCurrentTrackData->IsFreeFromBiasing()) { + // -- takes "ownership" (some track data created before, left free, + // reuse it): + fCurrentTrackData->fForceCollisionOperator = this; + } else { + // §§§ Would something be really wrong in this case ? Could this be that + // a process made a zero step ? + } + } else { + fCurrentTrackData = new GateOptrForceCollisionTrackData(this); + track->SetAuxiliaryTrackInformation(fForceCollisionModelID, + fCurrentTrackData); } - - // -- + fCurrentTrackData->fForceCollisionState = ForceCollisionState::toBeCloned; + fInitialTrackWeight = track->GetWeight(); + fCloningOperation->SetCloneWeights(0.0, fInitialTrackWeight); + return fCloningOperation; + } + + // -- return nullptr; } - -G4VBiasingOperation* GateOptrForceCollision::ProposeFinalStateBiasingOperation(const G4Track*, const G4BiasingProcessInterface* callingProcess) -{ - // -- Propose at final state generation the same operation which was proposed at GPIL level, - // -- (which is either the force free flight or the force interaction operation). +G4VBiasingOperation *GateOptrForceCollision::ProposeFinalStateBiasingOperation( + const G4Track *, const G4BiasingProcessInterface *callingProcess) { + // -- Propose at final state generation the same operation which was proposed + // at GPIL level, + // -- (which is either the force free flight or the force interaction + // operation). // -- count on the process interface to collect this operation: return callingProcess->GetCurrentOccurenceBiasingOperation(); } - -void GateOptrForceCollision::StartTracking( const G4Track* track ) -{ - fCurrentTrack = track; - fCurrentTrackData = nullptr; +void GateOptrForceCollision::StartTracking(const G4Track *track) { + fCurrentTrack = track; + fCurrentTrackData = nullptr; } - -void GateOptrForceCollision::EndTracking() -{ +void GateOptrForceCollision::EndTracking() { // -- check for consistency, operator should have cleaned the track: - if ( fCurrentTrackData != nullptr ) - { - if ( !fCurrentTrackData->IsFreeFromBiasing() ) - { - if ( (fCurrentTrack->GetTrackStatus() == fStopAndKill) || (fCurrentTrack->GetTrackStatus() == fKillTrackAndSecondaries) ) - { - G4ExceptionDescription ed; - ed << "Current track deleted while under biasing by " << GetName() << ". Will result in inconsistencies."; - G4Exception(" GateOptrForceCollision::EndTracking()", - "BIAS.GEN.18", - JustWarning, - ed); - } - } - } + if (fCurrentTrackData != nullptr) { + if (!fCurrentTrackData->IsFreeFromBiasing()) { + if ((fCurrentTrack->GetTrackStatus() == fStopAndKill) || + (fCurrentTrack->GetTrackStatus() == fKillTrackAndSecondaries)) { + G4ExceptionDescription ed; + ed << "Current track deleted while under biasing by " << GetName() + << ". Will result in inconsistencies."; + G4Exception(" GateOptrForceCollision::EndTracking()", "BIAS.GEN.18", + JustWarning, ed); + } + } + } } +void GateOptrForceCollision::OperationApplied( + const G4BiasingProcessInterface *callingProcess, G4BiasingAppliedCase BAC, + G4VBiasingOperation *operationApplied, const G4VParticleChange *) { -void GateOptrForceCollision::OperationApplied( const G4BiasingProcessInterface* callingProcess, - G4BiasingAppliedCase BAC, - G4VBiasingOperation* operationApplied, - const G4VParticleChange* ) -{ - - if ( fCurrentTrackData == nullptr ) - { - if ( BAC != BAC_None ) - { - G4ExceptionDescription ed; - ed << " Internal inconsistency : please submit bug report. " << G4endl; - G4Exception(" GateOptrForceCollision::OperationApplied(...)", - "BIAS.GEN.20.1", - JustWarning, - ed); - } - return; - } - - if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeCloned ) - { - fCurrentTrackData->fForceCollisionState = ForceCollisionState::toBeFreeFlight; - auto cloneData = new GateOptrForceCollisionTrackData( this ); - cloneData->fForceCollisionState = ForceCollisionState::toBeForced; - fCloningOperation->GetCloneTrack()->SetAuxiliaryTrackInformation(fForceCollisionModelID, cloneData); + if (fCurrentTrackData == nullptr) { + if (BAC != BAC_None) { + G4ExceptionDescription ed; + ed << " Internal inconsistency : please submit bug report. " << G4endl; + G4Exception(" GateOptrForceCollision::OperationApplied(...)", + "BIAS.GEN.20.1", JustWarning, ed); } - else if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeFreeFlight ) - { - if ( fFreeFlightOperations[callingProcess]->OperationComplete() ) fCurrentTrackData->Reset(); // -- off biasing for this track + return; + } + + if (fCurrentTrackData->fForceCollisionState == + ForceCollisionState::toBeCloned) { + fCurrentTrackData->fForceCollisionState = + ForceCollisionState::toBeFreeFlight; + auto cloneData = new GateOptrForceCollisionTrackData(this); + cloneData->fForceCollisionState = ForceCollisionState::toBeForced; + fCloningOperation->GetCloneTrack()->SetAuxiliaryTrackInformation( + fForceCollisionModelID, cloneData); + } else if (fCurrentTrackData->fForceCollisionState == + ForceCollisionState::toBeFreeFlight) { + if (fFreeFlightOperations[callingProcess]->OperationComplete()) + fCurrentTrackData->Reset(); // -- off biasing for this track + } else if (fCurrentTrackData->fForceCollisionState == + ForceCollisionState::toBeForced) { + if (operationApplied != fSharedForceInteractionOperation) { + G4ExceptionDescription ed; + ed << " Internal inconsistency : please submit bug report. " << G4endl; + G4Exception(" GateOptrForceCollision::OperationApplied(...)", + "BIAS.GEN.20.2", JustWarning, ed); } - else if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced ) - { - if ( operationApplied != fSharedForceInteractionOperation ) - { - G4ExceptionDescription ed; - ed << " Internal inconsistency : please submit bug report. " << G4endl; - G4Exception(" GateOptrForceCollision::OperationApplied(...)", - "BIAS.GEN.20.2", - JustWarning, - ed); - } - if ( fSharedForceInteractionOperation->GetInteractionOccured() ) - { - if ( operationApplied != fSharedForceInteractionOperation ) - { - G4ExceptionDescription ed; - ed << " Internal inconsistency : please submit bug report. " << G4endl; - G4Exception(" GateOptrForceCollision::OperationApplied(...)", - "BIAS.GEN.20.3", - JustWarning, - ed); - } - } + if (fSharedForceInteractionOperation->GetInteractionOccured()) { + if (operationApplied != fSharedForceInteractionOperation) { + G4ExceptionDescription ed; + ed << " Internal inconsistency : please submit bug report. " << G4endl; + G4Exception(" GateOptrForceCollision::OperationApplied(...)", + "BIAS.GEN.20.3", JustWarning, ed); + } } - else - { - if ( fCurrentTrackData->fForceCollisionState != ForceCollisionState::free ) - { - G4ExceptionDescription ed; - ed << " Internal inconsistency : please submit bug report. " << G4endl; - G4Exception(" GateOptrForceCollision::OperationApplied(...)", - "BIAS.GEN.20.4", - JustWarning, - ed); - } + } else { + if (fCurrentTrackData->fForceCollisionState != ForceCollisionState::free) { + G4ExceptionDescription ed; + ed << " Internal inconsistency : please submit bug report. " << G4endl; + G4Exception(" GateOptrForceCollision::OperationApplied(...)", + "BIAS.GEN.20.4", JustWarning, ed); } + } } - -void GateOptrForceCollision::OperationApplied( const G4BiasingProcessInterface* /*callingProcess*/, G4BiasingAppliedCase /*biasingCase*/, - G4VBiasingOperation* /*occurenceOperationApplied*/, G4double /*weightForOccurenceInteraction*/, - G4VBiasingOperation* finalStateOperationApplied, const G4VParticleChange* /*particleChangeProduced*/ ) -{ - - if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced ) - { - if ( finalStateOperationApplied != fSharedForceInteractionOperation ) - { - G4ExceptionDescription ed; - ed << " Internal inconsistency : please submit bug report. " << G4endl; - G4Exception(" GateOptrForceCollision::OperationApplied(...)", - "BIAS.GEN.20.5", - JustWarning, - ed); - } - if ( fSharedForceInteractionOperation->GetInteractionOccured() ) fCurrentTrackData->Reset(); // -- off biasing for this track - } - else - { +void GateOptrForceCollision::OperationApplied( + const G4BiasingProcessInterface * /*callingProcess*/, + G4BiasingAppliedCase /*biasingCase*/, + G4VBiasingOperation * /*occurenceOperationApplied*/, + G4double /*weightForOccurenceInteraction*/, + G4VBiasingOperation *finalStateOperationApplied, + const G4VParticleChange * /*particleChangeProduced*/) { + + if (fCurrentTrackData->fForceCollisionState == + ForceCollisionState::toBeForced) { + if (finalStateOperationApplied != fSharedForceInteractionOperation) { G4ExceptionDescription ed; ed << " Internal inconsistency : please submit bug report. " << G4endl; G4Exception(" GateOptrForceCollision::OperationApplied(...)", - "BIAS.GEN.20.6", - JustWarning, - ed); + "BIAS.GEN.20.5", JustWarning, ed); } - + if (fSharedForceInteractionOperation->GetInteractionOccured()) + fCurrentTrackData->Reset(); // -- off biasing for this track + } else { + G4ExceptionDescription ed; + ed << " Internal inconsistency : please submit bug report. " << G4endl; + G4Exception(" GateOptrForceCollision::OperationApplied(...)", + "BIAS.GEN.20.6", JustWarning, ed); + } } - diff --git a/core/opengate_core/opengate_lib/GateOptrForceCollision.h b/core/opengate_core/opengate_lib/GateOptrForceCollision.h index d720955bc..e70bc3b3d 100644 --- a/core/opengate_core/opengate_lib/GateOptrForceCollision.h +++ b/core/opengate_core/opengate_lib/GateOptrForceCollision.h @@ -39,7 +39,6 @@ // --------------------------------------------------------------- // Initial version Nov. 2013 M. Verderi - #ifndef GateOptrForceCollision_h #define GateOptrForceCollision_h 1 @@ -50,56 +49,66 @@ class GateOptnCloning; class G4VProcess; class G4BiasingProcessInterface; class G4ParticleDefinition; -#include -#include #include "G4ThreeVector.hh" #include "GateVActor.h" +#include +#include class GateOptrForceCollisionTrackData; -class GateOptrForceCollision : public G4VBiasingOperator,public GateVActor { +class GateOptrForceCollision : public G4VBiasingOperator, public GateVActor { public: GateOptrForceCollision(py::dict &user_info); ~GateOptrForceCollision(); - + public: // -- Mandatory from base class : - virtual G4VBiasingOperation* ProposeNonPhysicsBiasingOperation(const G4Track* track, const G4BiasingProcessInterface* callingProcess) final; - virtual G4VBiasingOperation* ProposeOccurenceBiasingOperation(const G4Track* track, const G4BiasingProcessInterface* callingProcess) final; - virtual G4VBiasingOperation* ProposeFinalStateBiasingOperation(const G4Track* track, const G4BiasingProcessInterface* callingProcess) final; + virtual G4VBiasingOperation *ProposeNonPhysicsBiasingOperation( + const G4Track *track, + const G4BiasingProcessInterface *callingProcess) final; + virtual G4VBiasingOperation *ProposeOccurenceBiasingOperation( + const G4Track *track, + const G4BiasingProcessInterface *callingProcess) final; + virtual G4VBiasingOperation *ProposeFinalStateBiasingOperation( + const G4Track *track, + const G4BiasingProcessInterface *callingProcess) final; // -- optional methods from base class: public: - virtual void Configure() final; - virtual void ConfigureForWorker() final; - virtual void StartRun() final; - virtual void StartTracking( const G4Track* track ) final; - virtual void ExitBiasing( const G4Track*, const G4BiasingProcessInterface* ) final {}; - virtual void EndTracking() final; + virtual void Configure() final; + virtual void ConfigureForWorker() final; + virtual void StartRun() final; + virtual void StartTracking(const G4Track *track) final; + virtual void ExitBiasing(const G4Track *, + const G4BiasingProcessInterface *) final {}; + virtual void EndTracking() final; void InitializeUserInfo(py::dict &user_info) override; -void AttachAllLogicalDaughtersVolumes(G4LogicalVolume *volume); + void AttachAllLogicalDaughtersVolumes(G4LogicalVolume *volume); // -- operation applied: - void OperationApplied( const G4BiasingProcessInterface* callingProcess, G4BiasingAppliedCase biasingCase, - G4VBiasingOperation* operationApplied, const G4VParticleChange* particleChangeProduced ) final; - void OperationApplied( const G4BiasingProcessInterface* callingProcess, G4BiasingAppliedCase biasingCase, - G4VBiasingOperation* occurenceOperationApplied, G4double weightForOccurenceInteraction, - G4VBiasingOperation* finalStateOperationApplied, const G4VParticleChange* particleChangeProduced ) final; - + void OperationApplied(const G4BiasingProcessInterface *callingProcess, + G4BiasingAppliedCase biasingCase, + G4VBiasingOperation *operationApplied, + const G4VParticleChange *particleChangeProduced) final; + void OperationApplied(const G4BiasingProcessInterface *callingProcess, + G4BiasingAppliedCase biasingCase, + G4VBiasingOperation *occurenceOperationApplied, + G4double weightForOccurenceInteraction, + G4VBiasingOperation *finalStateOperationApplied, + const G4VParticleChange *particleChangeProduced) final; -const G4String GetName(){ - return GateVActor::GetName(); -} + const G4String GetName() { return GateVActor::GetName(); } public: - G4int fForceCollisionModelID; - const G4Track* fCurrentTrack; - GateOptrForceCollisionTrackData* fCurrentTrackData; - std::map< const G4BiasingProcessInterface*, GateOptnForceFreeFlight* > fFreeFlightOperations; - GateOptnForceCommonTruncatedExp* fSharedForceInteractionOperation; - GateOptnCloning* fCloningOperation; - G4double fInitialTrackWeight; - G4bool fSetup; - const G4ParticleDefinition* fParticleToBias; + G4int fForceCollisionModelID; + const G4Track *fCurrentTrack; + GateOptrForceCollisionTrackData *fCurrentTrackData; + std::map + fFreeFlightOperations; + GateOptnForceCommonTruncatedExp *fSharedForceInteractionOperation; + GateOptnCloning *fCloningOperation; + G4double fInitialTrackWeight; + G4bool fSetup; + const G4ParticleDefinition *fParticleToBias; }; #endif diff --git a/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.cpp b/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.cpp index 26d8374e4..162c6e4f4 100644 --- a/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.cpp +++ b/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.cpp @@ -26,49 +26,52 @@ #include "GateOptrForceCollisionTrackData.h" #include "GateOptrForceCollision.h" -GateOptrForceCollisionTrackData::GateOptrForceCollisionTrackData( const GateOptrForceCollision* optr ) -: G4VAuxiliaryTrackInformation(), - fForceCollisionOperator( optr ) -{ +GateOptrForceCollisionTrackData::GateOptrForceCollisionTrackData( + const GateOptrForceCollision *optr) + : G4VAuxiliaryTrackInformation(), fForceCollisionOperator(optr) { fForceCollisionState = ForceCollisionState::free; } -GateOptrForceCollisionTrackData::~GateOptrForceCollisionTrackData() -{ - if ( fForceCollisionState != ForceCollisionState::free ) - { - G4ExceptionDescription ed; - ed << "Track deleted while under G4BOptrForceCollision biasing scheme of operator `"; - if ( fForceCollisionOperator == nullptr ) ed << "(none)"; else ed << "forceCollisionActor"; - ed <<"'. Will result in inconsistencies."; - G4Exception(" GateOptrForceCollisionTrackData::~GateOptrForceCollisionTrackData()", - "BIAS.GEN.19", - JustWarning, - ed); - } +GateOptrForceCollisionTrackData::~GateOptrForceCollisionTrackData() { + if (fForceCollisionState != ForceCollisionState::free) { + G4ExceptionDescription ed; + ed << "Track deleted while under G4BOptrForceCollision biasing scheme of " + "operator `"; + if (fForceCollisionOperator == nullptr) + ed << "(none)"; + else + ed << "forceCollisionActor"; + ed << "'. Will result in inconsistencies."; + G4Exception( + " GateOptrForceCollisionTrackData::~GateOptrForceCollisionTrackData()", + "BIAS.GEN.19", JustWarning, ed); + } } -void GateOptrForceCollisionTrackData::Print() const -{ +void GateOptrForceCollisionTrackData::Print() const { G4cout << " GateOptrForceCollisionTrackData object : " << this << G4endl; - G4cout << " Force collision operator : "; if ( fForceCollisionOperator == nullptr ) G4cout << "(none)"; else G4cout << "forceCollisionActor"; G4cout << G4endl; + G4cout << " Force collision operator : "; + if (fForceCollisionOperator == nullptr) + G4cout << "(none)"; + else + G4cout << "forceCollisionActor"; + G4cout << G4endl; G4cout << " Force collision state : "; - switch ( fForceCollisionState ) - { - case ForceCollisionState::free : - G4cout << "free from biasing "; - break; - case ForceCollisionState::toBeCloned : - G4cout << "to be cloned "; - break; - case ForceCollisionState::toBeForced : - G4cout << "to be interaction forced "; - break; - case ForceCollisionState::toBeFreeFlight : - G4cout << "to be free flight forced (under weight = 0) "; - break; - default: - break; - } + switch (fForceCollisionState) { + case ForceCollisionState::free: + G4cout << "free from biasing "; + break; + case ForceCollisionState::toBeCloned: + G4cout << "to be cloned "; + break; + case ForceCollisionState::toBeForced: + G4cout << "to be interaction forced "; + break; + case ForceCollisionState::toBeFreeFlight: + G4cout << "to be free flight forced (under weight = 0) "; + break; + default: + break; + } G4cout << G4endl; } diff --git a/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.h b/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.h index 8e0ecaf4f..a8d335023 100644 --- a/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.h +++ b/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.h @@ -26,7 +26,7 @@ // // // -------------------------------------------------------------------- -// GEANT 4 class header file +// GEANT 4 class header file // // Class Description: // Extends G4Track properties with information needed for the @@ -51,30 +51,30 @@ enum class ForceCollisionState { free, toBeCloned, toBeForced, toBeFreeFlight }; class GateOptrForceCollisionTrackData : public G4VAuxiliaryTrackInformation { -friend class GateOptrForceCollision; - + friend class GateOptrForceCollision; + public: - GateOptrForceCollisionTrackData( const GateOptrForceCollision* ); + GateOptrForceCollisionTrackData(const GateOptrForceCollision *); ~GateOptrForceCollisionTrackData(); - + // -- from base class: void Print() const; // -- Get methods: - G4bool IsFreeFromBiasing() const - { return ( fForceCollisionState == ForceCollisionState::free);} - // -- no set methods are provided : sets are made under exclusive control of G4BOptrForceCollision objects through friendness. - + G4bool IsFreeFromBiasing() const { + return (fForceCollisionState == ForceCollisionState::free); + } + // -- no set methods are provided : sets are made under exclusive control of + // G4BOptrForceCollision objects through friendness. + private: - const GateOptrForceCollision* fForceCollisionOperator; - ForceCollisionState fForceCollisionState; + const GateOptrForceCollision *fForceCollisionOperator; + ForceCollisionState fForceCollisionState; - void Reset() - { + void Reset() { fForceCollisionOperator = nullptr; - fForceCollisionState = ForceCollisionState::free; + fForceCollisionState = ForceCollisionState::free; } - }; #endif diff --git a/opengate/actors/miscactors.py b/opengate/actors/miscactors.py index 9b8cb0ac6..39a980d04 100644 --- a/opengate/actors/miscactors.py +++ b/opengate/actors/miscactors.py @@ -454,7 +454,7 @@ class ComptSplittingActor(GenericBiasingActorBase, g4.GateOptrComptSplittingActo rotation_vector_director: bool vector_director: list max_theta: float - mode : str + mode: str particles: list user_info_defaults = { @@ -490,8 +490,8 @@ class ComptSplittingActor(GenericBiasingActorBase, g4.GateOptrComptSplittingActo ), } - processes = (("compt",)) - particles = ("gamma") + processes = ("compt",) + particles = "gamma" mode = "physics" def __init__(self, *args, **kwargs): @@ -517,13 +517,15 @@ class BremSplittingActor(GenericBiasingActorBase, g4.GateBOptrBremSplittingActor # hints for IDE processes: list particles: list - mode :list + mode: list mode = "physics" - particles = ("e+","e-",) + particles = ( + "e+", + "e-", + ) processes = (("eBrem",), ("eBrem",)) - def __init__(self, *args, **kwargs): GenericBiasingActorBase.__init__(self, *args, **kwargs) self.__initcpp__() @@ -540,12 +542,11 @@ def initialize(self): class ForceCollisionActor(GenericBiasingActorBase, g4.GateOptrForceCollision): particles: str processes = list - mode : str - + mode: str mode = "both" particles = ("gamma",) - processes = [("compt","phot","conv","Rayl")] + processes = [("compt", "phot", "conv", "Rayl")] def __init__(self, *args, **kwargs): GenericBiasingActorBase.__init__(self, *args, **kwargs) @@ -560,7 +561,6 @@ def initialize(self): self.InitializeCpp() - process_cls(ActorOutputStatisticsActor) process_cls(SimulationStatisticsActor) process_cls(KillActor) diff --git a/opengate/engines.py b/opengate/engines.py index d260aee21..7e65b5011 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -399,11 +399,11 @@ def initialize_physics_biasing(self): mode = processes[1] if mode != "non physics": if len(process_list) > 0: - if mode == "physics" : + if mode == "physics": g4_biasing_physics.PhysicsBias(particle, process_list) elif mode == "both": - g4_biasing_physics.Bias(particle,process_list) - else : + g4_biasing_physics.Bias(particle, process_list) + else: g4_biasing_physics.NonPhysicsBias(process_list) self.g4_physics_list.RegisterPhysics(g4_biasing_physics) diff --git a/opengate/managers.py b/opengate/managers.py index 829b7e9a0..a9c4868c4 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -877,7 +877,7 @@ def get_biasing_particles_and_processes(self): all_particles = charged_particles.union({"gamma"}) # create a dictionary with sets as entries (to ensure uniqueness) - particles_processes= dict([(p, set()) for p in all_particles]) + particles_processes = dict([(p, set()) for p in all_particles]) for actor in self.simulation.actor_manager.actors.values(): if isinstance(actor, GenericBiasingActorBase): @@ -892,18 +892,18 @@ def get_biasing_particles_and_processes(self): f"Biasing actor {actor.name} wants to apply a bias to particle '{p_}'. " f"This is not possible. Compatible particles are: {list(all_particles)}. " ) - for i,p in enumerate(particles): + for i, p in enumerate(particles): if mode != "non physics": print(actor.processes) particles_processes[p].update(actor.processes[i]) - else : + else: particles_processes[p].update([None]) # convert the dictionary entries back from set to list print(particles_processes.items()) return dict( [ - (particle, [list(processes),mode]) + (particle, [list(processes), mode]) for particle, processes in particles_processes.items() ] ) From 4ffdfc29657d2d17ae59c2652d0d88f6013abba9 Mon Sep 17 00:00:00 2001 From: majacquet Date: Tue, 17 Dec 2024 17:54:33 +0100 Subject: [PATCH 3/3] bias mechanism bug correction --- opengate/managers.py | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/opengate/managers.py b/opengate/managers.py index 829b7e9a0..c1f91cf20 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -878,9 +878,10 @@ def get_biasing_particles_and_processes(self): # create a dictionary with sets as entries (to ensure uniqueness) particles_processes= dict([(p, set()) for p in all_particles]) - + activate_biasing = False for actor in self.simulation.actor_manager.actors.values(): if isinstance(actor, GenericBiasingActorBase): + activate_biasing = True particles = set() mode = actor.mode for particle in actor.particles: @@ -900,13 +901,21 @@ def get_biasing_particles_and_processes(self): particles_processes[p].update([None]) # convert the dictionary entries back from set to list - print(particles_processes.items()) - return dict( - [ - (particle, [list(processes),mode]) - for particle, processes in particles_processes.items() - ] - ) + if activate_biasing : + print(particles_processes.items()) + return dict( + [ + (particle, [list(processes),mode]) + for particle, processes in particles_processes.items() + ] + ) + else : + return dict( + [ + (particle, list(processes)) + for particle, processes in particles_processes.items() + ] + ) # New name, more specific def set_production_cut(self, volume_name, particle_name, value):