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