source: trunk/source/processes/biasing/src/G4WeightWindowProcess.cc@ 1101

Last change on this file since 1101 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 14.6 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: G4WeightWindowProcess.cc,v 1.4 2008/04/21 09:10:29 ahoward Exp $
28// GEANT4 tag $Name: geant4-09-02 $
29//
30// ----------------------------------------------------------------------
31// GEANT 4 class source file
32//
33// G4WeightWindowProcess.cc
34//
35// ----------------------------------------------------------------------
36
37#include "G4WeightWindowProcess.hh"
38#include "G4VWeightWindowAlgorithm.hh"
39#include "G4GeometryCell.hh"
40#include "G4SamplingPostStepAction.hh"
41#include "G4VTrackTerminator.hh"
42#include "G4PlaceOfAction.hh"
43#include "G4VWeightWindowStore.hh"
44
45#include "G4Step.hh"
46#include "G4Navigator.hh"
47#include "G4VTouchable.hh"
48#include "G4VPhysicalVolume.hh"
49#include "G4ParticleChange.hh"
50#include "G4PathFinder.hh"
51#include "G4TransportationManager.hh"
52#include "G4StepPoint.hh"
53#include "G4FieldTrackUpdator.hh"
54
55
56G4WeightWindowProcess::G4WeightWindowProcess(
57 const G4VWeightWindowAlgorithm &aWeightWindowAlgorithm,
58 const G4VWeightWindowStore &aWWStore,
59 const G4VTrackTerminator *TrackTerminator,
60 G4PlaceOfAction placeOfAction,
61 const G4String &aName, G4bool para)
62 : G4VProcess(aName),
63 fParticleChange(new G4ParticleChange),
64 fWeightWindowAlgorithm(aWeightWindowAlgorithm),
65 fWeightWindowStore(aWWStore),
66 fPostStepAction(0),
67 fPlaceOfAction(placeOfAction),
68 fGhostNavigator(0), fNavigatorID(-1), fFieldTrack('0')
69{
70 if (TrackTerminator)
71 {
72 fPostStepAction = new G4SamplingPostStepAction(*TrackTerminator);
73 }
74 else
75 {
76 fPostStepAction = new G4SamplingPostStepAction(*this);
77 }
78 if (!fParticleChange)
79 {
80 G4Exception("G4WeightWindowProcess::G4WeightWindowProcess()",
81 "FatalError", FatalException,
82 "Failed allocation of G4ParticleChange !");
83 }
84 G4VProcess::pParticleChange = fParticleChange;
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 paraflag = para;
99
100}
101
102G4WeightWindowProcess::~G4WeightWindowProcess()
103{
104
105 delete fPostStepAction;
106 delete fParticleChange;
107 delete fGhostStep;
108
109}
110
111
112//------------------------------------------------------
113//
114// SetParallelWorld
115//
116//------------------------------------------------------
117void G4WeightWindowProcess::
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
129void G4WeightWindowProcess::
130SetParallelWorld(G4VPhysicalVolume* parallelWorld)
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//------------------------------------------------------
146void G4WeightWindowProcess::StartTracking(G4Track* trk)
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("G4WeightWindowProcess::StartTracking",
159 "ProcParaWorld000",FatalException,
160 "G4WeightWindowProcess 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
186G4double G4WeightWindowProcess::
187PostStepGetPhysicalInteractionLength(const G4Track& ,
188 G4double ,
189 G4ForceCondition* condition)
190{
191// *condition = Forced;
192// return kInfinity;
193
194// *condition = StronglyForced;
195 *condition = Forced;
196 return DBL_MAX;
197}
198
199G4VParticleChange *
200G4WeightWindowProcess::PostStepDoIt(const G4Track &aTrack,
201 const G4Step &aStep)
202{
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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
220 }
221 else
222 {
223 // Do I need this ??????????????????????????????????????????????????????????
224 // fGhostNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
225 // ?????????????????????????????????????????????????????????????????????????
226
227 // fPathFinder->ReLocate(track.GetPosition());
228
229 // reuse the touchable
230 fNewGhostTouchable = fOldGhostTouchable;
231 }
232
233 fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
234 fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
235
236 }
237
238 if (aStep.GetStepLength() > kCarTolerance)
239 {
240// if ( ( fPlaceOfAction == onBoundaryAndCollision)
241// || ( (fPlaceOfAction == onBoundary) &&
242// (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary) )
243// || ( (fPlaceOfAction == onCollision) &&
244// (aStep.GetPostStepPoint()->GetStepStatus() != fGeomBoundary) ) )
245 if(paraflag) {
246 if ( ( fPlaceOfAction == onBoundaryAndCollision)
247 || ( (fPlaceOfAction == onBoundary) &&
248 (fGhostPostStepPoint->GetStepStatus() == fGeomBoundary) )
249 || ( (fPlaceOfAction == onCollision) &&
250 (fGhostPostStepPoint->GetStepStatus() != fGeomBoundary) ) )
251 {
252
253// G4StepPoint *postpoint = aStep.GetPostStepPoint();
254
255// G4GeometryCell postCell(*(postpoint->GetPhysicalVolume()),
256// postpoint->GetTouchable()->GetReplicaNumber());
257
258 G4GeometryCell postCell(*(fGhostPostStepPoint->GetPhysicalVolume()),
259 fGhostPostStepPoint->GetTouchable()->GetReplicaNumber());
260 G4Nsplit_Weight nw =
261 fWeightWindowAlgorithm.Calculate(aTrack.GetWeight(),
262 fWeightWindowStore.GetLowerWeitgh(postCell,
263 aTrack.GetKineticEnergy()));
264 fPostStepAction->DoIt(aTrack, fParticleChange, nw);
265 }
266 } else {
267 if ( ( fPlaceOfAction == onBoundaryAndCollision)
268 || ( (fPlaceOfAction == onBoundary) &&
269 (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary) )
270 || ( (fPlaceOfAction == onCollision) &&
271 (aStep.GetPostStepPoint()->GetStepStatus() != fGeomBoundary) ) )
272 {
273
274 G4StepPoint *postpoint = aStep.GetPostStepPoint();
275
276 G4GeometryCell postCell(*(postpoint->GetPhysicalVolume()),
277 postpoint->GetTouchable()->GetReplicaNumber());
278
279 G4Nsplit_Weight nw =
280 fWeightWindowAlgorithm.Calculate(aTrack.GetWeight(),
281 fWeightWindowStore.GetLowerWeitgh(postCell,
282 aTrack.GetKineticEnergy()));
283 fPostStepAction->DoIt(aTrack, fParticleChange, nw);
284 }
285 }
286 }
287 return fParticleChange;
288}
289
290void G4WeightWindowProcess::KillTrack() const
291{
292 fParticleChange->ProposeTrackStatus(fStopAndKill);
293}
294
295const G4String &G4WeightWindowProcess::GetName() const
296{
297 return G4VProcess::GetProcessName();
298}
299
300G4double G4WeightWindowProcess::
301AlongStepGetPhysicalInteractionLength(
302 const G4Track& track, G4double previousStepSize, G4double currentMinimumStep,
303 G4double& proposedSafety, G4GPILSelection* selection)
304{
305 if(paraflag) {
306 static G4FieldTrack endTrack('0');
307 static ELimited eLimited;
308
309 *selection = NotCandidateForSelection;
310 G4double returnedStep = DBL_MAX;
311
312 if (previousStepSize > 0.)
313 { fGhostSafety -= previousStepSize; }
314 // else
315 // { fGhostSafety = -1.; }
316 if (fGhostSafety < 0.) fGhostSafety = 0.0;
317
318 // ------------------------------------------
319 // Determination of the proposed STEP LENGTH:
320 // ------------------------------------------
321 if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.)
322 {
323 // I have no chance to limit
324 returnedStep = currentMinimumStep;
325 fOnBoundary = false;
326 proposedSafety = fGhostSafety - currentMinimumStep;
327 }
328 else // (currentMinimumStep > fGhostSafety: I may limit the Step)
329 {
330 G4FieldTrackUpdator::Update(&fFieldTrack,&track);
331 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
332 // ComputeStep
333 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
334 returnedStep
335 = fPathFinder->ComputeStep(fFieldTrack,currentMinimumStep,fNavigatorID,
336 track.GetCurrentStepNumber(),fGhostSafety,eLimited,
337 endTrack,track.GetVolume());
338 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
339 if(eLimited == kDoNot)
340 {
341 // Track is not on the boundary
342 fOnBoundary = false;
343 fGhostSafety = fGhostNavigator->ComputeSafety(endTrack.GetPosition());
344 }
345 else
346 {
347 // Track is on the boundary
348 fOnBoundary = true;
349 proposedSafety = fGhostSafety;
350 }
351 //xbug? proposedSafety = fGhostSafety;
352 if(eLimited == kUnique || eLimited == kSharedOther) {
353 *selection = CandidateForSelection;
354 }else if (eLimited == kSharedTransport) {
355 returnedStep *= (1.0 + 1.0e-9);
356 // Expand to disable its selection in Step Manager comparison
357 }
358
359 }
360
361 // ----------------------------------------------
362 // Returns the fGhostSafety as the proposedSafety
363 // The SteppingManager will take care of keeping
364 // the smallest one.
365 // ----------------------------------------------
366 return returnedStep;
367
368 } else {
369 return DBL_MAX;
370 // not sensible - time goes backwards! return -1.0;
371 }
372
373}
374
375G4double G4WeightWindowProcess::
376AtRestGetPhysicalInteractionLength(const G4Track& ,
377 G4ForceCondition*)
378{
379 return -1.0;
380}
381
382G4VParticleChange* G4WeightWindowProcess::
383AtRestDoIt(const G4Track&, const G4Step&)
384{
385 return 0;
386}
387
388G4VParticleChange* G4WeightWindowProcess::
389AlongStepDoIt(const G4Track& track, const G4Step&)
390{
391 // Dummy ParticleChange ie: does nothing
392 // Expecting G4Transportation to move the track
393 pParticleChange->Initialize(track);
394 return pParticleChange;
395
396 // return 0;
397}
398
399void G4WeightWindowProcess::CopyStep(const G4Step & step)
400{
401 fGhostStep->SetTrack(step.GetTrack());
402 fGhostStep->SetStepLength(step.GetStepLength());
403 fGhostStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
404 fGhostStep->SetControlFlag(step.GetControlFlag());
405
406 *fGhostPreStepPoint = *(step.GetPreStepPoint());
407 *fGhostPostStepPoint = *(step.GetPostStepPoint());
408
409//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
410// Set StepStatus for ghost world
411//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
412 if(fOnBoundary)
413 { fGhostPostStepPoint->SetStepStatus(fGeomBoundary); }
414 else if(fGhostPostStepPoint->GetStepStatus()==fGeomBoundary)
415 { fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc); }
416//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
417}
Note: See TracBrowser for help on using the repository browser.