| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /*
28 * Based on 'g4pbc'.
29 * Copyright (c) 2020 Amentum Pty Ltd
30 * team@amentum.space
31 * The original open-source version of this code
32 * may be found at https://github.com/amentumspace/g4pbc
33 * Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
34 * associated documentation files (the "Software"), to deal in the Software without restriction,
35 * including without limitation the rights to use, copy, modify, merge, publish, distribute,
36 * sublicense, and/or sell copies of the Software, and to permit persons to whom the Software
37 * is furnished to do so, subject to the following conditions:
38 * The above copyright notice and this permission notice shall be included in all copies
39 * or substantial portions of the Software.
40 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
41 * NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
42 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
43 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
44 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
45 *
46 *
47 */
48 #include "PeriodicBoundaryProcess.hh"
49
50 #include "LogicalVolumePeriodic.hh"
51 #include "ParticleChangeForPeriodic.hh"
52
53 #include "G4EventManager.hh"
54 #include "G4GeometryTolerance.hh"
55 #include "G4Navigator.hh"
56 #include "G4PathFinder.hh"
57 #include "G4TransportationManager.hh"
58 #include "G4UnitsTable.hh"
59 #include "G4VTrajectory.hh"
60 #include "G4ios.hh"
61
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63
64 PeriodicBoundaryProcess::PeriodicBoundaryProcess(const G4String& processName, G4ProcessType type,
65 G4bool per_x, G4bool per_y, G4bool per_z)
66 : G4VDiscreteProcess(processName, type), fPeriodicX(per_x), fPeriodicY(per_y), fPeriodicZ(per_z)
67 {
68 fkCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
69 pParticleChange = &fParticleChange;
70 G4PathFinder::GetInstance()->SetVerboseLevel(0);
71 }
72
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74
75 G4VParticleChange* PeriodicBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
76 {
77 fTheStatus = Undefined;
78
79 fParticleChange.InitializeForPostStep(aTrack);
80
81 const G4Step* pStep = &aStep;
82
83 G4bool isOnBoundary = (pStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
84
85 if (!isOnBoundary) {
86 fTheStatus = NotAtBoundary;
87 if (verboseLevel > 0) {
88 BoundaryProcessVerbose();
89 }
90 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
91 }
92
93 if (aTrack.GetParentID() == 0) // not primary particle (electron ?)
94 {
95 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
96 }
97
98 auto thePrePV = pStep->GetPreStepPoint()->GetPhysicalVolume();
99 auto thePostPV = pStep->GetPostStepPoint()->GetPhysicalVolume();
100
101 if (verboseLevel > 0) {
102 G4cout << " Particle at Boundary! " << G4endl;
103 if (thePrePV) {
104 G4cout << " thePrePV: " << thePrePV->GetName() << G4endl;
105 }
106 if (thePostPV) {
107 G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
108 }
109 G4cout << "step length " << aTrack.GetStepLength() << G4endl;
110 G4cout << "ParentID : " << aTrack.GetParentID() << " TrackID : " << aTrack.GetTrackID()
111 << " Position : " << aTrack.GetPosition()
112 << " aTrack Energy : " << G4BestUnit(aTrack.GetKineticEnergy(), "Energy") << G4endl;
113 }
114
115 // avoid trapped particles at boundaries by testing for minimum step length
116 if (aTrack.GetStepLength() <= fkCarTolerance / 2) {
117 fTheStatus = StepTooSmall;
118 if (verboseLevel > 0) {
119 BoundaryProcessVerbose();
120 }
121 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
122 }
123
124 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
125
126 // store the current values
127 fOldMomentum = aParticle->GetMomentumDirection();
128 fOldPolarization = aParticle->GetPolarization();
129 fOldPosition = pStep->GetPostStepPoint()->GetPosition();
130 fNewPosition = fOldPosition;
131
132 if (verboseLevel > 0) {
133 G4cout << " Old Momentum Direction: " << fOldMomentum << G4endl;
134 G4cout << " Old Position: " << fNewPosition << G4endl;
135 G4cout << "Get aTrack.GetParentID() : " << aTrack.GetParentID() << G4endl;
136 }
137
138 auto theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
139
140 G4bool valid = false;
141 // Use the new method for Exit Normal in global coordinates,
142 // which provides the normal more reliably.
143 // get from the real-world navigator (process not applied to parallel worlds)
144 fTheGlobalNormal = G4TransportationManager::GetTransportationManager()
145 ->GetNavigatorForTracking()
146 ->GetGlobalExitNormal(theGlobalPoint, &valid);
147
148 if (valid) {
149 fTheGlobalNormal = -fTheGlobalNormal;
150 }
151 else {
152 G4cout << "global normal " << fTheGlobalNormal << G4endl;
153 G4ExceptionDescription ed;
154 ed << " PeriodicBoundaryProcess::PostStepDoIt(): "
155 << " The Navigator reports that it returned an invalid normal" << G4endl;
156 G4Exception("G4PeriodicBoundaryProcess::PostStepDoIt", "PerBoun01", FatalException, ed,
157 "Invalid Surface Normal - Geometry must return valid surface normal");
158 }
159
160 G4bool isWrongDirection = fOldMomentum * fTheGlobalNormal > 0.0;
161 if (isWrongDirection) {
162 if (verboseLevel > 0) {
163 G4cout << "ftheGlobalNormal points in a wrong direction." << G4endl;
164 G4cout << "Invalid Surface Normal - Geometry must return valid surface \
165 normal pointing in the right direction"
166 << G4endl;
167 }
168
169 fTheGlobalNormal = -fTheGlobalNormal;
170 }
171
172 if (thePostPV == nullptr) {
173 G4ExceptionDescription ed;
174 ed << " PeriodicBoundaryProcess::PostStepDoIt(): "
175 << " thePostPV == nullptr" << G4endl;
176 G4Exception("G4PeriodicBoundaryProcess::PostStepDoIt", "PB14", FatalException, ed,
177 "Invalid thePostPV");
178 }
179 auto lvol = thePostPV->GetLogicalVolume();
180
181 if (lvol == nullptr) {
182 G4ExceptionDescription ed;
183 ed << " PeriodicBoundaryProcess::PostStepDoIt(): "
184 << " lvol == nullptr" << G4endl;
185 G4Exception("G4PeriodicBoundaryProcess::PostStepDoIt", "PB12", FatalException, ed,
186 "Invalid lvol");
187 }
188
189 if (verboseLevel > 0) {
190 G4cout << "Post step logical " << lvol->GetName() << G4endl;
191 }
192
193 G4LogicalVolume* dlvol = nullptr;
194
195 if (lvol->GetNoDaughters() > 0) {
196 if (verboseLevel > 0) {
197 G4cout << "eldest daughter " << lvol->GetDaughter(0)->GetName() << G4endl;
198 }
199
200 dlvol = lvol->GetDaughter(0)->GetLogicalVolume();
201 }
202
203 if (dlvol && dlvol->IsExtended()) {
204 if (verboseLevel > 0) {
205 G4cout << " Logical surface, periodic " << G4endl;
206 }
207
208 G4bool on_x = fTheGlobalNormal.isParallel(G4ThreeVector(1, 0, 0));
209 G4bool on_y = fTheGlobalNormal.isParallel(G4ThreeVector(0, 1, 0));
210 G4bool on_z = fTheGlobalNormal.isParallel(G4ThreeVector(0, 0, 1));
211
212 // make sure that we are at a plane
213 G4bool on_plane = (on_x || on_y || on_z);
214
215 if (!on_plane) {
216 G4ExceptionDescription ed;
217 ed << " G4PeriodicBoundaryProcess::ostStepDoIt(): "
218 << " The particle is not on a surface of the cyclic world" << G4endl;
219 G4Exception(
220 "G4PeriodicBoundaryProcess::PostStepDoIt", "Periodic01", FatalException, ed,
221 "Periodic boundary process must only occur for particle on periodic world surface");
222 }
223 else {
224 G4bool on_x_and_periodic = (on_x && fPeriodicX);
225 G4bool on_y_and_periodic = (on_y && fPeriodicY);
226 G4bool on_z_and_periodic = (on_z && fPeriodicZ);
227
228 G4bool on_a_periodic_plane = (on_x_and_periodic || on_y_and_periodic || on_z_and_periodic);
229
230 if (on_a_periodic_plane) {
231 if (verboseLevel > 0) {
232 G4cout << " on periodic plane " << G4endl;
233 }
234
235 fTheStatus = Cycling;
236
237 if (verboseLevel > 0) {
238 G4cout << " periodic " << G4endl;
239 }
240
241 if (verboseLevel > 0) {
242 G4cout << "Global normal " << fTheGlobalNormal << G4endl;
243 }
244
245 // translate a component of the position vector according to which plane we are on
246 if (on_x_and_periodic) {
247 fNewPosition.setX(-fNewPosition.x());
248 }
249 else if (on_y_and_periodic) {
250 fNewPosition.setY(-fNewPosition.y());
251 }
252 else if (on_z_and_periodic) {
253 fNewPosition.setZ(-fNewPosition.z());
254 }
255 else {
256 G4cout << "global normal does not belong to periodic plane!!" << G4endl;
257 }
258
259 fNewMomentum = fOldMomentum.unit();
260 fNewPolarization = fOldPolarization.unit();
261
262 if (verboseLevel > 0) {
263 G4cout << " New Position: " << fNewPosition << G4endl;
264 G4cout << " New Momentum Direction: " << fNewMomentum << G4endl;
265 G4cout << " New Polarization: " << fNewPolarization << G4endl;
266 BoundaryProcessVerbose();
267 }
268
269 fParticleChange.ProposeMomentumDirection(fNewMomentum);
270 fParticleChange.ProposePolarization(fNewPolarization);
271 fParticleChange.ProposePosition(fNewPosition);
272
273 G4PathFinder::GetInstance()->ReLocate(fNewPosition);
274 G4PathFinder::GetInstance()->ComputeSafety(fNewPosition);
275
276 // we must notify the navigator that we have moved the particle artificially
277 G4Navigator* gNavigator =
278 G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
279
280 // Inform the navigator that the previous Step calculated
281 // by the geometry was taken in its entirety.
282 gNavigator->SetGeometricallyLimitedStep();
283
284 // Search the geometrical hierarchy for the volumes deepest in the hierarchy
285 // containing the point in the global coordinate space.
286 gNavigator->LocateGlobalPointAndSetup(fNewPosition, &fNewMomentum, false,
287 false); // do not ignore direction
288
289 // Calculate the isotropic distance to the nearest boundary from the
290 // specified point in the global coordinate system.
291 gNavigator->ComputeSafety(fNewPosition);
292
293 // Locates the volume containing the specified global point.
294 // force drawing of the step prior to periodic the particle
295 auto evtm = G4EventManager::GetEventManager();
296 auto tckm = evtm->GetTrackingManager();
297 auto pTrajectory = tckm->GimmeTrajectory();
298 if (pTrajectory) {
299 pTrajectory->AppendStep(pStep);
300 }
301 }
302 }
303 }
304 return &fParticleChange;
305 }
306
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308
309 G4double PeriodicBoundaryProcess::GetMeanFreePath(const G4Track&, G4double,
310 G4ForceCondition* condition)
311 {
312 *condition = Forced;
313 return DBL_MAX;
314 }
315
316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
317
318 void PeriodicBoundaryProcess::BoundaryProcessVerbose()
319 {
320 if (fStatusMessages.count(fTheStatus) > 0) {
321 G4cout << PeriodicBoundaryProcess::fStatusMessages[fTheStatus] << G4endl;
322 }
323 }
324
325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
326