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 6a4fee2ba..38e57c843 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -347,7 +347,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); @@ -573,7 +573,7 @@ PYBIND11_MODULE(opengate_core, m) { init_GateProductionAndStoppingActor(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..9242484b5 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnCloning.cpp @@ -0,0 +1,45 @@ +// +// ******************************************************************** +// * 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..8a27ad2fa --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnCloning.h @@ -0,0 +1,92 @@ +// +// ******************************************************************** +// * 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 "G4ParticleChange.hh" +#include "G4VBiasingOperation.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..99da7b38b --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.cpp @@ -0,0 +1,165 @@ +// +// ******************************************************************** +// * 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 "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); + + 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_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..491a6f337 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnForceCommonTruncatedExp.h @@ -0,0 +1,124 @@ +// +// ******************************************************************** +// * 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 "G4ParticleChangeForNothing.hh" +#include "G4ThreeVector.hh" +#include "G4VBiasingOperation.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 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..f9c07d2d1 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.cpp @@ -0,0 +1,98 @@ +// +// ******************************************************************** +// * 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 "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); +} + +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..ab8e0641d --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptnForceFreeFlight.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. * +// ******************************************************************** +// +// +// +//--------------------------------------------------------------- +// +// 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 "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 &); + + // -- 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..c0b94b60b --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptrForceCollision.cpp @@ -0,0 +1,431 @@ +// +// ******************************************************************** +// * 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 "G4BiasingProcessInterface.hh" +#include "G4PhysicsModelCatalog.hh" +#include "GateOptrForceCollisionTrackData.h" + +#include "G4ILawCommonTruncatedExp.hh" +#include "GateOptnCloning.h" +#include "GateOptnForceCommonTruncatedExp.h" +#include "GateOptnForceFreeFlight.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::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..e70bc3b3d --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptrForceCollision.h @@ -0,0 +1,114 @@ +// +// ******************************************************************** +// * 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 "G4ThreeVector.hh" +#include "GateVActor.h" +#include +#include +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 + 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..162c6e4f4 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateOptrForceCollisionTrackData.cpp @@ -0,0 +1,77 @@ +// +// ******************************************************************** +// * 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..a8d335023 --- /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 d8a0c3873..39a980d04 100644 --- a/opengate/actors/miscactors.py +++ b/opengate/actors/miscactors.py @@ -403,7 +403,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 @@ -416,7 +416,6 @@ class SplittingActorBase(ActorBase): splitting_factor: int bias_primary_only: bool bias_only_once: bool - particles: list user_info_defaults = { "splitting_factor": ( @@ -428,28 +427,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 @@ -464,6 +454,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": ( @@ -499,21 +491,23 @@ class ComptSplittingActor(SplittingActorBase, g4.GateOptrComptSplittingActor): } 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 @@ -522,27 +516,47 @@ 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.", - }, - ), - } - - processes = ("eBrem",) + mode = "physics" + particles = ( + "e+", + "e-", + ) + processes = (("eBrem",), ("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() @@ -552,6 +566,7 @@ def initialize(self): 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 2a79cbaca..7e65b5011 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 613540cb0..8aaa05c33 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -97,9 +97,10 @@ SimulationStatisticsActor, KillActor, KillAccordingProcessesActor, - SplittingActorBase, + GenericBiasingActorBase, ComptSplittingActor, BremSplittingActor, + ForceCollisionActor, ) from .actors.digitizers import ( DigitizerAdderActor, @@ -135,6 +136,7 @@ "KillAccordingProcessesActor": KillAccordingProcessesActor, "BremSplittingActor": BremSplittingActor, "ComptSplittingActor": ComptSplittingActor, + "ForceCollisionActor": ForceCollisionActor, "DigitizerAdderActor": DigitizerAdderActor, "DigitizerBlurringActor": DigitizerBlurringActor, "DigitizerSpatialBlurringActor": DigitizerSpatialBlurringActor, @@ -875,36 +877,45 @@ 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]) + activate_biasing = False for actor in self.simulation.actor_manager.actors.values(): - if isinstance(actor, SplittingActorBase): + if isinstance(actor, GenericBiasingActorBase): + activate_biasing = True 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 - return dict( - [ - (particle, list(processes)) - 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): if volume_name == self.simulation.world.name: