source: trunk/source/processes/biasing/src/G4ImportanceProcess.cc@ 1216

Last change on this file since 1216 was 1196, checked in by garnier, 16 years ago

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