| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 co
32 * may be found at https://github.com/amentums
33 * Permission is hereby granted, free of charg
34 * associated documentation files (the "Softwa
35 * including without limitation the rights to
36 * sublicense, and/or sell copies of the Softw
37 * is furnished to do so, subject to the follo
38 * The above copyright notice and this permiss
39 * or substantial portions of the Software.
40 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT W
41 * NOT LIMITED TO THE WARRANTIES OF MERCHANTAB
42 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTH
43 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN A
44 * OUT OF OR IN CONNECTION WITH THE SOFTWARE O
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........oo
63
64 PeriodicBoundaryProcess::PeriodicBoundaryProce
65
66 : G4VDiscreteProcess(processName, type), fPe
67 {
68 fkCarTolerance = G4GeometryTolerance::GetIns
69 pParticleChange = &fParticleChange;
70 G4PathFinder::GetInstance()->SetVerboseLevel
71 }
72
73 //....oooOO0OOooo........oooOO0OOooo........oo
74
75 G4VParticleChange* PeriodicBoundaryProcess::Po
76 {
77 fTheStatus = Undefined;
78
79 fParticleChange.InitializeForPostStep(aTrack
80
81 const G4Step* pStep = &aStep;
82
83 G4bool isOnBoundary = (pStep->GetPostStepPoi
84
85 if (!isOnBoundary) {
86 fTheStatus = NotAtBoundary;
87 if (verboseLevel > 0) {
88 BoundaryProcessVerbose();
89 }
90 return G4VDiscreteProcess::PostStepDoIt(aT
91 }
92
93 if (aTrack.GetParentID() == 0) // not prima
94 {
95 return G4VDiscreteProcess::PostStepDoIt(aT
96 }
97
98 auto thePrePV = pStep->GetPreStepPoint()->Ge
99 auto thePostPV = pStep->GetPostStepPoint()->
100
101 if (verboseLevel > 0) {
102 G4cout << " Particle at Boundary! " << G4e
103 if (thePrePV) {
104 G4cout << " thePrePV: " << thePrePV->Ge
105 }
106 if (thePostPV) {
107 G4cout << " thePostPV: " << thePostPV->G
108 }
109 G4cout << "step length " << aTrack.GetStep
110 G4cout << "ParentID : " << aTrack.GetParen
111 << " Position : " << aTrack.GetPosi
112 << " aTrack Energy : " << G4BestUn
113 }
114
115 // avoid trapped particles at boundaries by
116 if (aTrack.GetStepLength() <= fkCarTolerance
117 fTheStatus = StepTooSmall;
118 if (verboseLevel > 0) {
119 BoundaryProcessVerbose();
120 }
121 return G4VDiscreteProcess::PostStepDoIt(aT
122 }
123
124 const G4DynamicParticle* aParticle = aTrack.
125
126 // store the current values
127 fOldMomentum = aParticle->GetMomentumDirecti
128 fOldPolarization = aParticle->GetPolarizatio
129 fOldPosition = pStep->GetPostStepPoint()->Ge
130 fNewPosition = fOldPosition;
131
132 if (verboseLevel > 0) {
133 G4cout << " Old Momentum Direction: " << f
134 G4cout << " Old Position: " << fNewPositio
135 G4cout << "Get aTrack.GetParentID() : " <<
136 }
137
138 auto theGlobalPoint = pStep->GetPostStepPoin
139
140 G4bool valid = false;
141 // Use the new method for Exit Normal in gl
142 // which provides the normal more reliabl
143 // get from the real-world navigator (proces
144 fTheGlobalNormal = G4TransportationManager::
145 ->GetNavigatorForTracki
146 ->GetGlobalExitNormal(t
147
148 if (valid) {
149 fTheGlobalNormal = -fTheGlobalNormal;
150 }
151 else {
152 G4cout << "global normal " << fTheGlobalNo
153 G4ExceptionDescription ed;
154 ed << " PeriodicBoundaryProcess::PostStepD
155 << " The Navigator reports that it retu
156 G4Exception("G4PeriodicBoundaryProcess::Po
157 "Invalid Surface Normal - Geom
158 }
159
160 G4bool isWrongDirection = fOldMomentum * fTh
161 if (isWrongDirection) {
162 if (verboseLevel > 0) {
163 G4cout << "ftheGlobalNormal points in a
164 G4cout << "Invalid Surface Normal - Geom
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::PostStepD
175 << " thePostPV == nullptr" << G4endl;
176 G4Exception("G4PeriodicBoundaryProcess::Po
177 "Invalid thePostPV");
178 }
179 auto lvol = thePostPV->GetLogicalVolume();
180
181 if (lvol == nullptr) {
182 G4ExceptionDescription ed;
183 ed << " PeriodicBoundaryProcess::PostStepD
184 << " lvol == nullptr" << G4endl;
185 G4Exception("G4PeriodicBoundaryProcess::Po
186 "Invalid lvol");
187 }
188
189 if (verboseLevel > 0) {
190 G4cout << "Post step logical " << lvol->Ge
191 }
192
193 G4LogicalVolume* dlvol = nullptr;
194
195 if (lvol->GetNoDaughters() > 0) {
196 if (verboseLevel > 0) {
197 G4cout << "eldest daughter " << lvol->Ge
198 }
199
200 dlvol = lvol->GetDaughter(0)->GetLogicalVo
201 }
202
203 if (dlvol && dlvol->IsExtended()) {
204 if (verboseLevel > 0) {
205 G4cout << " Logical surface, periodic "
206 }
207
208 G4bool on_x = fTheGlobalNormal.isParallel(
209 G4bool on_y = fTheGlobalNormal.isParallel(
210 G4bool on_z = fTheGlobalNormal.isParallel(
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::ostSt
218 << " The particle is not on a surface
219 G4Exception(
220 "G4PeriodicBoundaryProcess::PostStepDo
221 "Periodic boundary process must only o
222 }
223 else {
224 G4bool on_x_and_periodic = (on_x && fPer
225 G4bool on_y_and_periodic = (on_y && fPer
226 G4bool on_z_and_periodic = (on_z && fPer
227
228 G4bool on_a_periodic_plane = (on_x_and_p
229
230 if (on_a_periodic_plane) {
231 if (verboseLevel > 0) {
232 G4cout << " on periodic plane " << G
233 }
234
235 fTheStatus = Cycling;
236
237 if (verboseLevel > 0) {
238 G4cout << " periodic " << G4endl;
239 }
240
241 if (verboseLevel > 0) {
242 G4cout << "Global normal " << fTheGl
243 }
244
245 // translate a component of the positi
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 be
257 }
258
259 fNewMomentum = fOldMomentum.unit();
260 fNewPolarization = fOldPolarization.un
261
262 if (verboseLevel > 0) {
263 G4cout << " New Position: " << fNewP
264 G4cout << " New Momentum Direction:
265 G4cout << " New Polarization:
266 BoundaryProcessVerbose();
267 }
268
269 fParticleChange.ProposeMomentumDirecti
270 fParticleChange.ProposePolarization(fN
271 fParticleChange.ProposePosition(fNewPo
272
273 G4PathFinder::GetInstance()->ReLocate(
274 G4PathFinder::GetInstance()->ComputeSa
275
276 // we must notify the navigator that w
277 G4Navigator* gNavigator =
278 G4TransportationManager::GetTranspor
279
280 // Inform the navigator that the previ
281 // by the geometry was taken in its en
282 gNavigator->SetGeometricallyLimitedSte
283
284 // Search the geometrical hierarchy fo
285 // containing the point in the global
286 gNavigator->LocateGlobalPointAndSetup(
287
288
289 // Calculate the isotropic distance to
290 // specified point in the global coord
291 gNavigator->ComputeSafety(fNewPosition
292
293 // Locates the volume containing the s
294 // force drawing of the step prior to
295 auto evtm = G4EventManager::GetEventMa
296 auto tckm = evtm->GetTrackingManager()
297 auto pTrajectory = tckm->GimmeTrajecto
298 if (pTrajectory) {
299 pTrajectory->AppendStep(pStep);
300 }
301 }
302 }
303 }
304 return &fParticleChange;
305 }
306
307 //....oooOO0OOooo........oooOO0OOooo........oo
308
309 G4double PeriodicBoundaryProcess::GetMeanFreeP
310
311 {
312 *condition = Forced;
313 return DBL_MAX;
314 }
315
316 //....oooOO0OOooo........oooOO0OOooo........oo
317
318 void PeriodicBoundaryProcess::BoundaryProcessV
319 {
320 if (fStatusMessages.count(fTheStatus) > 0) {
321 G4cout << PeriodicBoundaryProcess::fStatus
322 }
323 }
324
325 //....oooOO0OOooo........oooOO0OOooo........oo
326