| 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 // G4ParallelGeometriesLimiterProcess
27 // --------------------------------------------------------------------
28
29 #include "G4ios.hh"
30 #include "G4ParallelGeometriesLimiterProcess.hh"
31 #include "G4BiasingProcessSharedData.hh"
32 #include "G4ProcessManager.hh"
33 #include "G4TransportationManager.hh"
34 #include "G4PathFinder.hh"
35 #include "G4FieldTrackUpdator.hh"
36
37 #include "G4SystemOfUnits.hh"
38
39 G4ParallelGeometriesLimiterProcess::
40 G4ParallelGeometriesLimiterProcess(const G4String& processName)
41 : G4VProcess(processName, fParallel)
42 {
43 // -- Process Sub Type ?
44
45 fPathFinder = G4PathFinder::GetInstance();
46 fTransportationManager = G4TransportationManager::GetTransportationManager();
47 }
48
49 // ----------------------------
50 // -- Add/Remove world volumes:
51 // ----------------------------
52 void G4ParallelGeometriesLimiterProcess::
53 AddParallelWorld(const G4String& parallelWorldName)
54 {
55 // -- Refuse adding parallel geometry during tracking time:
56 if (fIsTrackingTime)
57 {
58 G4ExceptionDescription ed;
59 ed << "G4ParallelGeometriesLimiterProcess `" << GetProcessName()
60 << "': adding a parallel world volume at tracking time is not allowed."
61 << G4endl;
62 G4Exception("G4ParallelGeometriesLimiterProcess::AddParallelWorld(..)",
63 "BIAS.GEN.21", JustWarning, ed, "Call ignored.");
64 return;
65 }
66 else
67 {
68 G4VPhysicalVolume* newWorld = fTransportationManager->IsWorldExisting( parallelWorldName );
69
70 // -- Fatal exception if requested world does not exist:
71 if (newWorld == nullptr)
72 {
73 G4ExceptionDescription tellWhatIsWrong;
74 tellWhatIsWrong << "Volume `" << parallelWorldName
75 << "' is not a parallel world nor the mass world volume."
76 << G4endl;
77 G4Exception("G4ParallelGeometriesLimiterProcess::SetWorldVolume(..)",
78 "BIAS.GEN.22", FatalException, tellWhatIsWrong);
79 }
80
81 // -- Protection against adding the mass geometry world as parallel world:
82 if ( newWorld == fTransportationManager->GetNavigatorForTracking()->GetWorldVolume() )
83 {
84 G4ExceptionDescription ed;
85 ed << "G4ParallelGeometriesLimiterProcess `" << GetProcessName()
86 << "': trying to add the world volume for tracking as a parallel world."
87 << G4endl;
88 G4Exception("G4ParallelGeometriesLimiterProcess::AddParallelWorld(..)",
89 "BIAS.GEN.23", JustWarning, ed, "Call ignored.");
90 return;
91 }
92
93 // -- Add parallel world, taking care it is not in the list yet:
94 G4bool isNew = true;
95 for ( auto knownWorld : fParallelWorlds )
96 {
97 if ( knownWorld == newWorld ) { isNew = false; }
98 }
99 if ( isNew )
100 {
101 fParallelWorlds.push_back( newWorld );
102 }
103 else
104 {
105 G4ExceptionDescription ed;
106 ed << "G4ParallelGeometriesLimiterProcess `" << GetProcessName()
107 << "': trying to re-add the parallel world volume `"
108 << parallelWorldName << "'." << G4endl;
109 G4Exception("G4ParallelGeometriesLimiterProcess::AddParallelWorld(..)",
110 "BIAS.GEN.24", JustWarning, ed, "Call ignored.");
111 return;
112 }
113 }
114 }
115
116 void G4ParallelGeometriesLimiterProcess::
117 RemoveParallelWorld(const G4String& parallelWorldName)
118 {
119 // -- Refuse refuse removing parallel geometry during tracking time:
120 if (fIsTrackingTime)
121 {
122 G4ExceptionDescription ed;
123 ed << "G4ParallelGeometriesLimiterProcess `" << GetProcessName()
124 << "': removing a parallel world volume at tracking time is not allowed."
125 << G4endl;
126 G4Exception("G4ParallelGeometriesLimiterProcess::RemoveParallelWorld(..)",
127 "BIAS.GEN.25", JustWarning, ed, "Call ignored.");
128 return;
129 }
130 else
131 {
132 G4VPhysicalVolume* newWorld = fTransportationManager->IsWorldExisting( parallelWorldName );
133 if (newWorld == nullptr)
134 {
135 G4ExceptionDescription ed;
136 ed << "G4ParallelGeometriesLimiterProcess `" << GetProcessName()
137 << "': trying to remove an inexisting parallel world '"
138 << parallelWorldName << "'." << G4endl;
139 G4Exception("G4ParallelGeometriesLimiterProcess::RemoveParallelWorld(..)",
140 "BIAS.GEN.26", JustWarning, ed, "Call ignored.");
141 return;
142 }
143
144 // -- get position of world volume in list:
145 std::size_t iWorld = 0;
146 for ( auto knownWorld : fParallelWorlds )
147 {
148 if ( knownWorld == newWorld ) break;
149 ++iWorld;
150 }
151
152 if ( iWorld == fParallelWorlds.size() )
153 {
154 G4ExceptionDescription ed;
155 ed << "G4ParallelGeometriesLimiterProcess `" << GetProcessName()
156 << "': trying to remove an non-registerered parallel world '"
157 << parallelWorldName << "'." << G4endl;
158 G4Exception("G4ParallelGeometriesLimiterProcess::RemoveParallelWorld(..)",
159 "BIAS.GEN.27", JustWarning, ed, "Call ignored.");
160 return;
161 }
162 // -- remove from vector:
163 fParallelWorlds.erase( fParallelWorlds.begin() + iWorld );
164 }
165 }
166
167 // --------------------
168 // Start/End tracking:
169 // --------------------
170 void G4ParallelGeometriesLimiterProcess::StartTracking(G4Track* track)
171 {
172 fIsTrackingTime = true;
173
174 // -- fetch the navigators, their indeces, and activate:
175 fParallelWorldNavigators.clear();
176 fParallelWorldNavigatorIndeces.clear();
177 fParallelWorldSafeties.clear();
178 fParallelWorldIsLimiting.clear();
179 fParallelWorldWasLimiting.clear();
180 fCurrentVolumes.clear();
181 fPreviousVolumes.clear();
182 for ( auto parallelWorld : fParallelWorlds )
183 {
184 fParallelWorldNavigators.push_back( fTransportationManager->GetNavigator( parallelWorld ) );
185 fParallelWorldNavigatorIndeces.push_back( fTransportationManager->ActivateNavigator( fParallelWorldNavigators.back() ) );
186 fParallelWorldSafeties.push_back( 0.0 );
187 fParallelWorldIsLimiting.push_back( false );
188 fParallelWorldWasLimiting.push_back( false );
189 }
190
191 fPathFinder->PrepareNewTrack( track->GetPosition(), track->GetMomentumDirection() );
192 // -- Does it work at this level, after "PrepareNewTrack" above ?
193 for ( auto navigatorIndex : fParallelWorldNavigatorIndeces )
194 {
195 fPreviousVolumes.push_back( nullptr );
196 fCurrentVolumes .push_back( fPathFinder->GetLocatedVolume( navigatorIndex ) );
197 }
198
199 // -- will force updating safety:
200 fParallelWorldSafety = 0.0;
201 for ( std::size_t i = 0 ; i < fParallelWorldNavigatorIndeces.size() ; ++i )
202 {
203 fParallelWorldSafeties[i] = 0.0;
204 }
205 }
206
207 void G4ParallelGeometriesLimiterProcess::EndTracking()
208 {
209 fIsTrackingTime = false;
210 for ( auto parallelWorldNavigator : fParallelWorldNavigators )
211 {
212 fTransportationManager->DeActivateNavigator( parallelWorldNavigator );
213 }
214 }
215
216 G4double G4ParallelGeometriesLimiterProcess::
217 PostStepGetPhysicalInteractionLength(const G4Track&, G4double,
218 G4ForceCondition* condition)
219 {
220 // -- push previous step limitation flags and volumes:
221 // -- consider switching pointers insteads of making copies of std::vector's:
222 fParallelWorldWasLimiting = fParallelWorldIsLimiting;
223 fPreviousVolumes = fCurrentVolumes;
224
225 // -- update volumes:
226 std::size_t i = 0;
227 for ( auto navigatorIndex : fParallelWorldNavigatorIndeces )
228 {
229 fCurrentVolumes[i++] = fPathFinder->GetLocatedVolume( navigatorIndex );
230 }
231
232 *condition = NotForced;
233 return DBL_MAX;
234 }
235
236 G4double G4ParallelGeometriesLimiterProcess::
237 AlongStepGetPhysicalInteractionLength(const G4Track& track,
238 G4double previousStepSize,
239 G4double currentMinimumStep,
240 G4double& proposedSafety,
241 G4GPILSelection* selection)
242 {
243 // -- Init:
244 // -- Note that the returnedStep must be physically meaningful,
245 // -- even if we return NotCandidateForSelection as condition;
246 // -- the reason is that the stepping manager always takes the smallest
247 // -- alongstep among the returned ones (point related
248 // -- to geometry step length wrt to true path length).
249 *selection = NotCandidateForSelection;
250 G4double returnedStep = DBL_MAX;
251
252 // -- G4FieldTrack and ELimited:
253 static G4ThreadLocal G4FieldTrack* endTrack_MT = nullptr;
254 if (!endTrack_MT) endTrack_MT = new G4FieldTrack ('0');
255 G4FieldTrack& endTrack = *endTrack_MT;
256
257 static G4ThreadLocal ELimited* eLimited_MT = nullptr;
258 if (!eLimited_MT) eLimited_MT = new ELimited;
259 ELimited &eLimited = *eLimited_MT;
260
261 // -------------------
262 // -- Update safeties:
263 // -------------------
264 if ( previousStepSize > 0.0 )
265 {
266 for ( auto& parallelWorldSafety : fParallelWorldSafeties )
267 {
268 parallelWorldSafety -= previousStepSize;
269 if ( parallelWorldSafety < 0. ) { parallelWorldSafety = 0.0; }
270 fParallelWorldSafety = parallelWorldSafety < fParallelWorldSafety
271 ? parallelWorldSafety : fParallelWorldSafety;
272 }
273 }
274
275 // ------------------------------------------
276 // Determination of the proposed step length:
277 // ------------------------------------------
278 if ( ( currentMinimumStep <= fParallelWorldSafety )
279 && ( currentMinimumStep > 0. ) )
280 {
281 // -- No chance to limit the step, as proposed move inside safety
282
283 returnedStep = currentMinimumStep;
284 proposedSafety = fParallelWorldSafety - currentMinimumStep;
285 }
286 else
287 {
288 // -- Proposed move exceeds common safety, need to state
289 G4double smallestReturnedStep = -1.0;
290 ELimited eLimitedForSmallestStep = kDoNot;
291 for ( std::size_t i = 0 ; i < fParallelWorldNavigatorIndeces.size() ; ++i )
292 {
293 // -- Update safety of geometries having safety smaller than current minimum step
294 if ( currentMinimumStep >= fParallelWorldSafeties[i] )
295 {
296 G4FieldTrackUpdator::Update(&fFieldTrack, &track);
297 G4double tmpReturnedStep = fPathFinder->ComputeStep(fFieldTrack,
298 currentMinimumStep,
299 fParallelWorldNavigatorIndeces[i],
300 track.GetCurrentStepNumber(),
301 fParallelWorldSafeties[i],
302 eLimited, endTrack, track.GetVolume());
303
304 if ( ( smallestReturnedStep < 0.0 ) || ( tmpReturnedStep <= smallestReturnedStep ) )
305 {
306 smallestReturnedStep = tmpReturnedStep;
307 eLimitedForSmallestStep = eLimited;
308 }
309
310 if (eLimited == kDoNot)
311 {
312 // -- Step not limited by this geometry
313 fParallelWorldSafeties[i] = fParallelWorldNavigators[i]->ComputeSafety(endTrack.GetPosition());
314 fParallelWorldIsLimiting[i] = false;
315 }
316 else
317 {
318 fParallelWorldIsLimiting[i] = true;
319 }
320 }
321
322 // -- update with smallest safety:
323 fParallelWorldSafety = fParallelWorldSafeties[i] < fParallelWorldSafety
324 ? fParallelWorldSafeties[i] : fParallelWorldSafety;
325 }
326
327 // -- no geometry limitation among all geometries, can return currentMinimumStep (or DBL_MAX):
328 // -- Beware : the returnedStep must be physically meaningful, even if we say "NotCandidateForSelection" !
329 if ( eLimitedForSmallestStep == kDoNot )
330 {
331 returnedStep = currentMinimumStep;
332 }
333 // -- proposed step length of limiting geometry:
334 if ( eLimitedForSmallestStep == kUnique ||
335 eLimitedForSmallestStep == kSharedOther )
336 {
337 *selection = CandidateForSelection;
338 returnedStep = smallestReturnedStep;
339 }
340 else if ( eLimitedForSmallestStep == kSharedTransport )
341 {
342 // -- Expand to disable its selection in Step Manager comparison
343 returnedStep = smallestReturnedStep* (1.0 + 1.0e-9);
344 }
345
346 // -- and smallest safety among geometries:
347 proposedSafety = fParallelWorldSafety ;
348 }
349
350 // -- returns step length, and proposedSafety
351 return returnedStep;
352 }
353
354 G4VParticleChange* G4ParallelGeometriesLimiterProcess::
355 AlongStepDoIt( const G4Track& track, const G4Step& )
356 {
357 fDummyParticleChange.Initialize(track);
358 return &fDummyParticleChange;
359 }
360
361 void G4ParallelGeometriesLimiterProcess::
362 SetProcessManager(const G4ProcessManager* mgr)
363 {
364 G4BiasingProcessSharedData *sharedData(nullptr);
365
366 // -- initialize sharedData pointer:
367 if ( G4BiasingProcessSharedData::fSharedDataMap.Find(mgr) == G4BiasingProcessSharedData::fSharedDataMap.End() )
368 {
369 sharedData = new G4BiasingProcessSharedData( mgr );
370 G4BiasingProcessSharedData::fSharedDataMap[mgr] = sharedData;
371 }
372 else
373 {
374 sharedData = G4BiasingProcessSharedData::fSharedDataMap[mgr] ;
375 }
376
377 // -- add itself to the shared data:
378 if ( sharedData->fParallelGeometriesLimiterProcess == nullptr )
379 {
380 sharedData->fParallelGeometriesLimiterProcess = this;
381 }
382 else
383 {
384 G4ExceptionDescription ed;
385 ed << " Trying to add more than one G4ParallelGeometriesLimiterProcess process to the process manager "
386 << mgr
387 << " (process manager for `" << mgr->GetParticleType()->GetParticleName()
388 << "'). Only one is needed. Call ignored." << G4endl;
389 G4Exception(" G4ParallelGeometriesLimiterProcess::SetProcessManager(..)",
390 "BIAS.GEN.29", JustWarning, ed);
391 }
392 }
393
394 G4int G4ParallelGeometriesLimiterProcess::
395 GetParallelWorldIndex( const G4VPhysicalVolume* parallelWorld ) const
396 {
397 G4int toReturn = -1;
398 G4int iWorld = 0;
399 for ( auto world : fParallelWorlds )
400 {
401 if ( world == parallelWorld )
402 {
403 toReturn = iWorld;
404 break;
405 }
406 ++iWorld;
407 }
408 return toReturn;
409 }
410
411 G4int G4ParallelGeometriesLimiterProcess::
412 GetParallelWorldIndex( const G4String& parallelWorldName ) const
413 {
414 G4VPhysicalVolume* aWorld = fTransportationManager->IsWorldExisting( parallelWorldName );
415 // note aWorld might be nullptr
416 return GetParallelWorldIndex( aWorld );
417 }
418