Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4ImportanceProcess.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$
28//
29// ----------------------------------------------------------------------
30// GEANT 4 class source file
31//
32// G4ImportanceProcess.cc
33//
34// ----------------------------------------------------------------------
35
38#include "G4GeometryCell.hh"
40#include "G4VTrackTerminator.hh"
41#include "G4VIStore.hh"
42
43#include "G4Step.hh"
44#include "G4Navigator.hh"
45#include "G4VTouchable.hh"
46#include "G4VPhysicalVolume.hh"
47#include "G4ParticleChange.hh"
48#include "G4PathFinder.hh"
50#include "G4StepPoint.hh"
52
53
55G4ImportanceProcess(const G4VImportanceAlgorithm &aImportanceAlgorithm,
56 const G4VIStore &aIstore,
57 const G4VTrackTerminator *TrackTerminator,
58 const G4String &aName, G4bool para)
59 : G4VProcess(aName),
60 fParticleChange(new G4ParticleChange),
61 fImportanceAlgorithm(aImportanceAlgorithm),
62 fIStore(aIstore),
63 fPostStepAction(0),
64 fGhostNavigator(0), fNavigatorID(-1), fFieldTrack('0'),
65 paraflag(para)
66
67{
68 if (TrackTerminator)
69 {
70 fPostStepAction = new G4SamplingPostStepAction(*TrackTerminator);
71 }
72 else
73 {
74 fPostStepAction = new G4SamplingPostStepAction(*this);
75 }
76 if (!fParticleChange)
77 {
78 G4Exception("G4ImportanceProcess::G4ImportanceProcess()",
79 "FatalError", FatalException,
80 "Failed allocation of G4ParticleChange !");
81 }
82 G4VProcess::pParticleChange = fParticleChange;
83
84
85 fGhostStep = new G4Step();
86 fGhostPreStepPoint = fGhostStep->GetPreStepPoint();
87 fGhostPostStepPoint = fGhostStep->GetPostStepPoint();
88
90 fPathFinder = G4PathFinder::GetInstance();
91
92 if (verboseLevel>0)
93 {
94 G4cout << GetProcessName() << " is created " << G4endl;
95 }
96
97 G4cout << " importance process paraflag is: " << paraflag << G4endl;
98
99}
100
102{
103
104 delete fPostStepAction;
105 delete fParticleChange;
106 delete fGhostStep;
107
108}
109
110
111
112//------------------------------------------------------
113//
114// SetParallelWorld
115//
116//------------------------------------------------------
118SetParallelWorld(G4String parallelWorldName)
119{
120//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
121// Get pointers of the parallel world and its navigator
122//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
123 fGhostWorldName = parallelWorldName;
124 fGhostWorld = fTransportationManager->GetParallelWorld(fGhostWorldName);
125 fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
126//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
127}
128
131{
132//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
133// Get pointer of navigator
134//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
135 fGhostWorldName = parallelWorld->GetName();
136 fGhostWorld = parallelWorld;
137 fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
138//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
139}
140
141//------------------------------------------------------
142//
143// StartTracking
144//
145//------------------------------------------------------
147{
148//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
149// Activate navigator and get the navigator ID
150//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
151// G4cout << " G4ParallelWorldScoringProcess::StartTracking" << G4endl;
152
153 if(paraflag) {
154 if(fGhostNavigator)
155 { fNavigatorID = fTransportationManager->ActivateNavigator(fGhostNavigator); }
156 else
157 {
158 G4Exception("G4ImportanceProcess::StartTracking",
159 "ProcParaWorld000",FatalException,
160 "G4ImportanceProcess is used for tracking without having a parallel world assigned");
161 }
162//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
163
164// G4cout << "G4ParallelWorldScoringProcess::StartTracking <<<<<<<<<<<<<<<<<< " << G4endl;
165//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
166// Let PathFinder initialize
167//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
168 fPathFinder->PrepareNewTrack(trk->GetPosition(),trk->GetMomentumDirection());
169//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
170
171//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
172// Setup initial touchables for the first step
173//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
174 fOldGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
175 fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
176 fNewGhostTouchable = fOldGhostTouchable;
177 fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
178
179 // Initialize
180 fGhostSafety = -1.;
181 fOnBoundary = false;
182 }
183
184}
185
186
189 G4double ,
191{
192// *condition = Forced;
193// return kInfinity;
194
195// *condition = StronglyForced;
196 *condition = Forced;
197 return DBL_MAX;
198}
199
202 const G4Step &aStep)
203{
204 fParticleChange->Initialize(aTrack);
205
206 if(paraflag) {
207 fOldGhostTouchable = fGhostPostStepPoint->GetTouchableHandle();
208 //xbug? fOnBoundary = false;
209 CopyStep(aStep);
210
211 if(fOnBoundary)
212 {
213//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
214// Locate the point and get new touchable
215//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
216 //?? fPathFinder->Locate(step.GetPostStepPoint()->GetPosition(),
217 //?? step.GetPostStepPoint()->GetMomentumDirection());
218 fNewGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
219 //AH G4cout << " on boundary " << G4endl;
220//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
221 }
222 else
223 {
224 // Do I need this ??????????????????????????????????????????????????????????
225 // fGhostNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
226 // ?????????????????????????????????????????????????????????????????????????
227
228 // fPathFinder->ReLocate(track.GetPosition());
229
230 // reuse the touchable
231 fNewGhostTouchable = fOldGhostTouchable;
232 //AH G4cout << " NOT on boundary " << G4endl;
233 }
234
235 fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
236 fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
237
238// if ( (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
239// && (aStep.GetStepLength() > kCarTolerance) )
240// {
241// if (aTrack.GetTrackStatus()==fStopAndKill)
242// {
243// G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
244// << " StopAndKill track." << G4endl;
245// }
246
247// G4StepPoint *prepoint = aStep.GetPreStepPoint();
248// G4StepPoint *postpoint = aStep.GetPostStepPoint();
249// G4GeometryCell prekey(*(prepoint->GetPhysicalVolume()),
250// prepoint->GetTouchable()->GetReplicaNumber());
251// G4GeometryCell postkey(*(postpoint->GetPhysicalVolume()),
252// postpoint->GetTouchable()->GetReplicaNumber());
253
254
255 if ( (fGhostPostStepPoint->GetStepStatus() == fGeomBoundary)
256 && (aStep.GetStepLength() > kCarTolerance) )
257 {
258 if (aTrack.GetTrackStatus()==fStopAndKill)
259 {
260 G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
261 << " StopAndKill track. on boundary" << G4endl;
262 }
263
264 G4GeometryCell prekey(*(fGhostPreStepPoint->GetPhysicalVolume()),
265 fGhostPreStepPoint->GetTouchable()->GetReplicaNumber());
266 G4GeometryCell postkey(*(fGhostPostStepPoint->GetPhysicalVolume()),
267 fGhostPostStepPoint->GetTouchable()->GetReplicaNumber());
268
269 //AH
270 /*
271 G4cout << G4endl;
272 G4cout << G4endl;
273 G4cout << " inside parallel importance process " << aTrack.GetCurrentStepNumber() << G4endl;
274 G4cout << G4endl;
275 G4cout << G4endl;
276 G4cout << " prekey: " << fGhostPreStepPoint->GetPhysicalVolume()->GetName() << " replica:"
277 << fGhostPreStepPoint->GetTouchable()->GetReplicaNumber() << G4endl;
278 G4cout << " prekey ISTORE: " << fIStore.GetImportance(prekey) << G4endl;
279 G4cout << " postkey: " << G4endl;
280 G4cout << " postkey ISTORE: " << fIStore.GetImportance(postkey) << G4endl;
281 */
282 //AH
283 G4Nsplit_Weight nw = fImportanceAlgorithm.
284 Calculate(fIStore.GetImportance(prekey),
285 fIStore.GetImportance(postkey),
286 aTrack.GetWeight());
287 //AH
288 /*
289 G4cout << " prekey weight: " << fIStore.GetImportance(prekey)
290 << " postkey weight: " << fIStore.GetImportance(postkey)
291 << " split weight: " << nw << G4endl;
292 G4cout << " before poststepaction " << G4endl;
293 */
294 //AH
295 fPostStepAction->DoIt(aTrack, fParticleChange, nw);
296 //AH G4cout << " after post step do it " << G4endl;
297 }
298 } else {
299 if ( (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
300 && (aStep.GetStepLength() > kCarTolerance) )
301 {
302 //AH G4cout << " inside non-parallel importance process " << G4endl;
303 if (aTrack.GetTrackStatus()==fStopAndKill)
304 {
305 G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
306 << " StopAndKill track. on boundary non-parallel" << G4endl;
307 }
308
309 G4StepPoint *prepoint = aStep.GetPreStepPoint();
310 G4StepPoint *postpoint = aStep.GetPostStepPoint();
311
312 G4GeometryCell prekey(*(prepoint->GetPhysicalVolume()),
313 prepoint->GetTouchable()->GetReplicaNumber());
314 G4GeometryCell postkey(*(postpoint->GetPhysicalVolume()),
315 postpoint->GetTouchable()->GetReplicaNumber());
316
317 G4Nsplit_Weight nw = fImportanceAlgorithm.
318 Calculate(fIStore.GetImportance(prekey),
319 fIStore.GetImportance(postkey),
320 aTrack.GetWeight());
321 //AH
322 /*
323 G4cout << " prekey weight: " << fIStore.GetImportance(prekey)
324 << " postkey weight: " << fIStore.GetImportance(postkey)
325 << " split weight: " << nw << G4endl;
326 G4cout << " before poststepaction 2 " << G4endl;
327 */
328 //AH
329 fPostStepAction->DoIt(aTrack, fParticleChange, nw);
330 //AH G4cout << " after poststepaction 2 " << G4endl;
331 }
332 }
333 return fParticleChange;
334}
335
337{
338 fParticleChange->ProposeTrackStatus(fStopAndKill);
339}
340
342{
343 return theProcessName;
344}
345
348 const G4Track& track, G4double previousStepSize, G4double currentMinimumStep,
349 G4double& proposedSafety, G4GPILSelection* selection)
350{
351
352 //AH G4cout << " along step physical interaction length " << G4endl;
353
354 if(paraflag) {
355 static G4FieldTrack endTrack('0');
356 static ELimited eLimited;
357
358 *selection = NotCandidateForSelection;
359 G4double returnedStep = DBL_MAX;
360
361 if (previousStepSize > 0.)
362 { fGhostSafety -= previousStepSize; }
363 // else
364 // { fGhostSafety = -1.; }
365 if (fGhostSafety < 0.) fGhostSafety = 0.0;
366
367 // ------------------------------------------
368 // Determination of the proposed STEP LENGTH:
369 // ------------------------------------------
370 if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.)
371 {
372 // I have no chance to limit
373 returnedStep = currentMinimumStep;
374 fOnBoundary = false;
375 proposedSafety = fGhostSafety - currentMinimumStep;
376 //AH G4cout << " step not limited, why? " << G4endl;
377 }
378 else // (currentMinimumStep > fGhostSafety: I may limit the Step)
379 {
380 G4FieldTrackUpdator::Update(&fFieldTrack,&track);
381 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
382 // ComputeStep
383 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
384 returnedStep
385 = fPathFinder->ComputeStep(fFieldTrack,currentMinimumStep,fNavigatorID,
386 track.GetCurrentStepNumber(),fGhostSafety,eLimited,
387 endTrack,track.GetVolume());
388 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
389 if(eLimited == kDoNot)
390 {
391 //AH G4cout << " computestep came back with not-boundary " << G4endl;
392 // Track is not on the boundary
393 fOnBoundary = false;
394 fGhostSafety = fGhostNavigator->ComputeSafety(endTrack.GetPosition());
395 }
396 else
397 {
398 // Track is on the boundary
399 //AH G4cout << " FOUND A BOUNDARY ! " << track.GetCurrentStepNumber() << G4endl;
400 fOnBoundary = true;
401 // proposedSafety = fGhostSafety;
402 }
403 proposedSafety = fGhostSafety;
404 if(eLimited == kUnique || eLimited == kSharedOther) {
405 *selection = CandidateForSelection;
406 }else if (eLimited == kSharedTransport) {
407 returnedStep *= (1.0 + 1.0e-9);
408 // Expand to disable its selection in Step Manager comparison
409 }
410
411 }
412
413 // ----------------------------------------------
414 // Returns the fGhostSafety as the proposedSafety
415 // The SteppingManager will take care of keeping
416 // the smallest one.
417 // ----------------------------------------------
418 return returnedStep;
419
420 } else {
421
422 return DBL_MAX;
423
424 }
425
426}
427
431{
432 return -1.0;
433}
434
436AtRestDoIt(const G4Track&, const G4Step&)
437{
438 return 0;
439}
440
442AlongStepDoIt(const G4Track& aTrack, const G4Step& )
443{
444 // Dummy ParticleChange ie: does nothing
445 // Expecting G4Transportation to move the track
446 //AH G4cout << " along step do it " << G4endl;
448 return pParticleChange;
449 // return 0;
450}
451
452void G4ImportanceProcess::CopyStep(const G4Step & step)
453{
454 fGhostStep->SetTrack(step.GetTrack());
455 fGhostStep->SetStepLength(step.GetStepLength());
457 fGhostStep->SetControlFlag(step.GetControlFlag());
458
459 *fGhostPreStepPoint = *(step.GetPreStepPoint());
460 *fGhostPostStepPoint = *(step.GetPostStepPoint());
461
462//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
463// Set StepStatus for ghost world
464//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
465 if(fOnBoundary)
466 { fGhostPostStepPoint->SetStepStatus(fGeomBoundary); }
467 else if(fGhostPostStepPoint->GetStepStatus()==fGeomBoundary)
468 { fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc); }
469//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
470}
G4double condition(const G4ErrorSymMatrix &m)
@ FatalException
G4ForceCondition
@ Forced
G4GPILSelection
@ CandidateForSelection
@ NotCandidateForSelection
ELimited
@ kDoNot
@ kUnique
@ kSharedOther
@ kSharedTransport
@ fGeomBoundary
Definition: G4StepStatus.hh:54
@ fPostStepDoItProc
Definition: G4StepStatus.hh:60
@ fStopAndKill
double G4double
Definition: G4Types.hh:64
bool G4bool
Definition: G4Types.hh:67
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
static void Update(G4FieldTrack *, const G4Track *)
G4ThreeVector GetPosition() const
virtual G4VParticleChange * AlongStepDoIt(const G4Track &, const G4Step &)
virtual G4VParticleChange * AtRestDoIt(const G4Track &, const G4Step &)
virtual void KillTrack() const
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
void SetParallelWorld(G4String parallelWorldName)
G4ImportanceProcess(const G4VImportanceAlgorithm &aImportanceAlgorithm, const G4VIStore &aIstore, const G4VTrackTerminator *TrackTerminator, const G4String &aName="ImportanceProcess", G4bool para=false)
void StartTracking(G4Track *)
virtual G4double AtRestGetPhysicalInteractionLength(const G4Track &, G4ForceCondition *)
virtual const G4String & GetName() const
virtual G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4double pProposedMaxLength=DBL_MAX, const G4bool keepState=true)
virtual void Initialize(const G4Track &)
G4double ComputeStep(const G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4int navigatorId, G4int stepNo, G4double &pNewSafety, ELimited &limitedStep, G4FieldTrack &EndState, G4VPhysicalVolume *currentVolume)
void PrepareNewTrack(const G4ThreeVector &position, const G4ThreeVector &direction, G4VPhysicalVolume *massStartVol=0)
static G4PathFinder * GetInstance()
Definition: G4PathFinder.cc:57
G4TouchableHandle CreateTouchableHandle(G4int navId) const
void DoIt(const G4Track &aTrack, G4ParticleChange *aParticleChange, const G4Nsplit_Weight &nw)
G4StepStatus GetStepStatus() const
const G4VTouchable * GetTouchable() const
void SetStepStatus(const G4StepStatus aValue)
void SetTouchableHandle(const G4TouchableHandle &apValue)
const G4TouchableHandle & GetTouchableHandle() const
G4VPhysicalVolume * GetPhysicalVolume() const
Definition: G4Step.hh:78
G4SteppingControl GetControlFlag() const
G4Track * GetTrack() const
void SetStepLength(G4double value)
G4StepPoint * GetPreStepPoint() const
G4double GetStepLength() const
G4double GetTotalEnergyDeposit() const
void SetControlFlag(G4SteppingControl StepControlFlag)
void SetTotalEnergyDeposit(G4double value)
G4StepPoint * GetPostStepPoint() const
void SetTrack(G4Track *value)
G4TrackStatus GetTrackStatus() const
G4VPhysicalVolume * GetVolume() const
G4double GetWeight() const
const G4ThreeVector & GetPosition() const
G4int GetCurrentStepNumber() const
const G4ThreeVector & GetMomentumDirection() const
G4VPhysicalVolume * GetParallelWorld(const G4String &worldName)
static G4TransportationManager * GetTransportationManager()
G4int ActivateNavigator(G4Navigator *aNavigator)
G4Navigator * GetNavigator(const G4String &worldName)
virtual G4double GetImportance(const G4GeometryCell &gCell) const =0
void ProposeTrackStatus(G4TrackStatus status)
virtual void Initialize(const G4Track &)
const G4String & GetName() const
G4int verboseLevel
Definition: G4VProcess.hh:368
G4String theProcessName
Definition: G4VProcess.hh:335
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
const G4String & GetProcessName() const
Definition: G4VProcess.hh:379
virtual G4int GetReplicaNumber(G4int depth=0) const
Definition: G4VTouchable.cc:58
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
#define DBL_MAX
Definition: templates.hh:83