File indexing completed on 2025-02-23 09:22:10
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027 #include "BoundedBrownianAction.hh"
0028
0029 #include "G4DNABoundingBox.hh"
0030 #include "G4Molecule.hh"
0031 #include "G4Track.hh"
0032
0033
0034
0035 BoundedBrownianAction::BoundedBrownianAction() : G4VUserBrownianAction() {}
0036
0037
0038 G4double BoundedBrownianAction::GetDistanceToBoundary(const G4Track& track)
0039 {
0040 if (!fpBoundingBox->contains(track.GetPosition())) {
0041 G4ExceptionDescription errMsg;
0042 errMsg << "Point is out of box : " << *fpBoundingBox
0043 << " of particle : " << GetIT(track)->GetName() << "(" << track.GetTrackID()
0044 << ") : " << track.GetPosition();
0045 G4Exception(
0046 "BoundedBrownianAction::GetDistanceToBoundary"
0047 "BoundedBrownianAction",
0048 "BoundedBrownianAction", FatalErrorInArgument, errMsg);
0049 }
0050
0051 auto dx = std::min(track.GetPosition().getX() - fpBoundingBox->Getxlo(),
0052 fpBoundingBox->Getxhi() - track.GetPosition().getX());
0053 auto dy = std::min(track.GetPosition().getY() - fpBoundingBox->Getylo(),
0054 fpBoundingBox->Getyhi() - track.GetPosition().getY());
0055 auto dz = std::min(track.GetPosition().getZ() - fpBoundingBox->Getzlo(),
0056 fpBoundingBox->Getzhi() - track.GetPosition().getZ());
0057 return std::min({dx, dy, dz});
0058 }
0059
0060
0061
0062 void BoundedBrownianAction::Transport(G4ThreeVector& nextposition, G4Track*)
0063 {
0064 BouncingAction(nextposition);
0065 }
0066
0067
0068
0069 void BoundedBrownianAction::BouncingAction(G4ThreeVector& nextposition) const
0070 {
0071
0072
0073 G4double RxM = fpBoundingBox->Getxhi();
0074 G4double RyM = fpBoundingBox->Getyhi();
0075 G4double RzM = fpBoundingBox->Getzhi();
0076
0077 G4double Rxm = fpBoundingBox->Getxlo();
0078 G4double Rym = fpBoundingBox->Getylo();
0079 G4double Rzm = fpBoundingBox->Getzlo();
0080 G4double Lx = RxM - Rxm;
0081 G4double Ly = RyM - Rym;
0082 G4double Lz = RzM - Rzm;
0083 G4double trunx = trunc(std::abs(nextposition.getX() - Rxm) / Lx);
0084 G4double truny = trunc(std::abs(nextposition.getY() - Rym) / Ly);
0085 G4double trunz = trunc(std::abs(nextposition.getZ() - Rzm) / Lz);
0086 G4double hx = std::fmod(trunx, 2);
0087 G4double hy = std::fmod(truny, 2);
0088 G4double hz = std::fmod(trunz, 2);
0089 G4double x = Rxm + hx * Lx + (1 - 2 * hx) * std::abs(std::fmod(nextposition.getX() - Rxm, Lx));
0090 G4double y = Rym + hy * Ly + (1 - 2 * hy) * std::abs(std::fmod(nextposition.getY() - Rym, Ly));
0091 G4double z = Rzm + hz * Lz + (1 - 2 * hz) * std::abs(std::fmod(nextposition.getZ() - Rzm, Lz));
0092 nextposition.set(x, y, z);
0093 }
0094
0095