Force collision

core/opengate_core/g4_bindings/pyPhysicsLists.cpp

@@ -90,6 +90,14 @@ namespace py = pybind11;
       .def("PhysicsBias",                                                      \
            py::overload_cast<const G4String &, const std::vector<G4String> &>( \
                &G4GenericBiasingPhysics::PhysicsBias),                         \
+           py::return_value_policy::reference_internal)                        \
+      .def("Bias",                                                             \
+           py::overload_cast<const G4String &, const std::vector<G4String> &>( \
+               &G4GenericBiasingPhysics::Bias),                                \
+           py::return_value_policy::reference_internal)                        \
+      .def("NonPhysicsBias",                                                   \
+           py::overload_cast<const G4String &>(                                \
+               &G4GenericBiasingPhysics::NonPhysicsBias),                      \
            py::return_value_policy::reference_internal);
 
 namespace pyPhysicsLists {

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

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;
+}

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

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 G4VProcess *, G4double>::const_iterator it =
+           fCrossSections.begin();
+       it != fCrossSections.end(); it++) {
+    sigmaSelect += (*it).second;
+    if (sigmaRand <= sigmaSelect) {
+      fProcessToApply = (*it).first;
+      break;
+    }
+  }
+}