source: trunk/source/processes/biasing/src/G4WeightCutOffProcess.cc@ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 15 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 13.7 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: G4WeightCutOffProcess.cc,v 1.3 2008/04/21 09:10:29 ahoward Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30// ----------------------------------------------------------------------
31// GEANT 4 class source file
32//
33// G4WeightCutOffProcess.cc
34//
35// ----------------------------------------------------------------------
36
37#include "G4WeightCutOffProcess.hh"
38//#include "G4VScorer.hh"
39#include "G4GeometryCellStep.hh"
40//#include "G4GCellFinder.hh"
41#include "G4TouchableHandle.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
55G4WeightCutOffProcess::
56G4WeightCutOffProcess(G4double wsurvival,
57 G4double wlimit,
58 G4double isource,
59 G4VIStore *istore,
60 // const G4VGCellFinder &aGCellFinder,
61 const G4String &aName, G4bool para)
62 : G4VProcess(aName),
63 fParticleChange(new G4ParticleChange),
64 fWeightSurvival(wsurvival),
65 fWeightLimit(wlimit),
66 fSourceImportance(isource),
67 fIStore(istore),
68 // fGCellFinder(aGCellFinder),
69 fGhostNavigator(0), fNavigatorID(-1), fFieldTrack('0')
70{
71 if (!fParticleChange)
72 {
73 G4Exception("G4WeightCutOffProcess::G4WeightCutOffProcess()",
74 "FatalError", FatalException,
75 "Failed to allocate G4ParticleChange !");
76 }
77
78 G4VProcess::pParticleChange = fParticleChange;
79
80 fGhostStep = new G4Step();
81 fGhostPreStepPoint = fGhostStep->GetPreStepPoint();
82 fGhostPostStepPoint = fGhostStep->GetPostStepPoint();
83
84 fTransportationManager = G4TransportationManager::GetTransportationManager();
85 fPathFinder = G4PathFinder::GetInstance();
86
87 if (verboseLevel>0)
88 {
89 G4cout << GetProcessName() << " is created " << G4endl;
90 }
91
92 paraflag = para;
93
94}
95
96G4WeightCutOffProcess::~G4WeightCutOffProcess()
97{
98 delete fParticleChange;
99 delete fGhostStep;
100}
101
102
103//------------------------------------------------------
104//
105// SetParallelWorld
106//
107//------------------------------------------------------
108void G4WeightCutOffProcess::
109SetParallelWorld(G4String parallelWorldName)
110{
111//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
112// Get pointers of the parallel world and its navigator
113//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
114 fGhostWorldName = parallelWorldName;
115 fGhostWorld = fTransportationManager->GetParallelWorld(fGhostWorldName);
116 fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
117//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
118}
119
120void G4WeightCutOffProcess::
121SetParallelWorld(G4VPhysicalVolume* parallelWorld)
122{
123//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
124// Get pointer of navigator
125//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
126 fGhostWorldName = parallelWorld->GetName();
127 fGhostWorld = parallelWorld;
128 fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
129//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
130}
131
132//------------------------------------------------------
133//
134// StartTracking
135//
136//------------------------------------------------------
137void G4WeightCutOffProcess::StartTracking(G4Track* trk)
138{
139//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
140// Activate navigator and get the navigator ID
141//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
142// G4cout << " G4ParallelWorldScoringProcess::StartTracking" << G4endl;
143
144 if(paraflag) {
145 if(fGhostNavigator)
146 { fNavigatorID = fTransportationManager->ActivateNavigator(fGhostNavigator); }
147 else
148 {
149 G4Exception("G4WeightCutOffProcess::StartTracking",
150 "ProcParaWorld000",FatalException,
151 "G4WeightCutOffProcess is used for tracking without having a parallel world assigned");
152 }
153//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
154
155// G4cout << "G4ParallelWorldScoringProcess::StartTracking <<<<<<<<<<<<<<<<<< " << G4endl;
156//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
157// Let PathFinder initialize
158//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
159 fPathFinder->PrepareNewTrack(trk->GetPosition(),trk->GetMomentumDirection());
160//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
161
162//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
163// Setup initial touchables for the first step
164//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
165 fOldGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
166 fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
167 fNewGhostTouchable = fOldGhostTouchable;
168 fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
169
170 // Initialize
171 fGhostSafety = -1.;
172 fOnBoundary = false;
173 }
174}
175
176
177G4double G4WeightCutOffProcess::
178PostStepGetPhysicalInteractionLength(const G4Track &,
179 G4double, G4ForceCondition* condition)
180{
181// *condition = Forced;
182// return kInfinity;
183
184// *condition = StronglyForced;
185 *condition = Forced;
186 return DBL_MAX;
187}
188
189G4VParticleChange *
190G4WeightCutOffProcess::PostStepDoIt(const G4Track& aTrack,
191 const G4Step &aStep)
192{
193 fParticleChange->Initialize(aTrack);
194
195 if(paraflag) {
196 fOldGhostTouchable = fGhostPostStepPoint->GetTouchableHandle();
197 //xbug? fOnBoundary = false;
198 CopyStep(aStep);
199
200 if(fOnBoundary)
201 {
202//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
203// Locate the point and get new touchable
204//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
205 //?? fPathFinder->Locate(step.GetPostStepPoint()->GetPosition(),
206 //?? step.GetPostStepPoint()->GetMomentumDirection());
207 fNewGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
208//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
209 }
210 else
211 {
212 // Do I need this ??????????????????????????????????????????????????????????
213 // fGhostNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
214 // ?????????????????????????????????????????????????????????????????????????
215
216 // fPathFinder->ReLocate(track.GetPosition());
217
218 // reuse the touchable
219 fNewGhostTouchable = fOldGhostTouchable;
220 }
221
222 fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
223 fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
224
225 }
226
227 if(paraflag) {
228 G4GeometryCell postCell(*(fGhostPostStepPoint->GetPhysicalVolume()),
229 fGhostPostStepPoint->GetTouchable()->GetReplicaNumber());
230
231
232 // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(aStep);
233 // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(fGhostStep);
234 G4double R = fSourceImportance;
235 if (fIStore)
236 {
237 G4double i = fIStore->GetImportance(postCell);
238 if (i>0)
239 {
240 R/=i;
241 }
242 }
243 G4double w = aTrack.GetWeight();
244 if (w<R*fWeightLimit)
245 {
246 G4double ws = fWeightSurvival*R;
247 G4double p = w/(ws);
248 if (G4UniformRand()<p)
249 {
250 fParticleChange->ProposeTrackStatus(fStopAndKill);
251 }
252 else
253 {
254 fParticleChange->ProposeWeight(ws);
255 }
256 }
257 } else {
258
259 G4GeometryCell postCell(*(aStep.GetPostStepPoint()->GetPhysicalVolume()),
260 aStep.GetPostStepPoint()->GetTouchable()->GetReplicaNumber());
261
262 // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(aStep);
263 // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(fGhostStep);
264 G4double R = fSourceImportance;
265 if (fIStore)
266 {
267 G4double i = fIStore->GetImportance(postCell);
268 if (i>0)
269 {
270 R/=i;
271 }
272 }
273 G4double w = aTrack.GetWeight();
274 if (w<R*fWeightLimit)
275 {
276 G4double ws = fWeightSurvival*R;
277 G4double p = w/(ws);
278 if (G4UniformRand()<p)
279 {
280 fParticleChange->ProposeTrackStatus(fStopAndKill);
281 }
282 else
283 {
284 fParticleChange->ProposeWeight(ws);
285 }
286 }
287 }
288
289 return fParticleChange;
290
291}
292
293const G4String &G4WeightCutOffProcess::GetName() const
294{
295 return theProcessName;
296}
297
298G4double G4WeightCutOffProcess::
299AlongStepGetPhysicalInteractionLength(
300 const G4Track& track, G4double previousStepSize, G4double currentMinimumStep,
301 G4double& proposedSafety, G4GPILSelection* selection)
302{
303 if(paraflag) {
304 static G4FieldTrack endTrack('0');
305 static ELimited eLimited;
306
307 *selection = NotCandidateForSelection;
308 G4double returnedStep = DBL_MAX;
309
310 if (previousStepSize > 0.)
311 { fGhostSafety -= previousStepSize; }
312 // else
313 // { fGhostSafety = -1.; }
314 if (fGhostSafety < 0.) fGhostSafety = 0.0;
315
316 // ------------------------------------------
317 // Determination of the proposed STEP LENGTH:
318 // ------------------------------------------
319 if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.)
320 {
321 // I have no chance to limit
322 returnedStep = currentMinimumStep;
323 fOnBoundary = false;
324 proposedSafety = fGhostSafety - currentMinimumStep;
325 }
326 else // (currentMinimumStep > fGhostSafety: I may limit the Step)
327 {
328 G4FieldTrackUpdator::Update(&fFieldTrack,&track);
329 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
330 // ComputeStep
331 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
332 returnedStep
333 = fPathFinder->ComputeStep(fFieldTrack,currentMinimumStep,fNavigatorID,
334 track.GetCurrentStepNumber(),fGhostSafety,eLimited,
335 endTrack,track.GetVolume());
336 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
337 if(eLimited == kDoNot)
338 {
339 // Track is not on the boundary
340 fOnBoundary = false;
341 fGhostSafety = fGhostNavigator->ComputeSafety(endTrack.GetPosition());
342 }
343 else
344 {
345 // Track is on the boundary
346 fOnBoundary = true;
347 proposedSafety = fGhostSafety;
348 }
349 //xbug? proposedSafety = fGhostSafety;
350 if(eLimited == kUnique || eLimited == kSharedOther) {
351 *selection = CandidateForSelection;
352 }else if (eLimited == kSharedTransport) {
353 returnedStep *= (1.0 + 1.0e-9);
354 // Expand to disable its selection in Step Manager comparison
355 }
356
357 }
358
359 // ----------------------------------------------
360 // Returns the fGhostSafety as the proposedSafety
361 // The SteppingManager will take care of keeping
362 // the smallest one.
363 // ----------------------------------------------
364 return returnedStep;
365
366 } else {
367 return DBL_MAX;
368 //not sensible! return -1.0;
369 }
370
371}
372
373
374G4double G4WeightCutOffProcess::
375AtRestGetPhysicalInteractionLength(const G4Track& ,
376 G4ForceCondition*)
377{
378 return -1.0;
379}
380
381G4VParticleChange*
382G4WeightCutOffProcess::AtRestDoIt(const G4Track&, const G4Step&)
383{
384 return 0;
385}
386
387G4VParticleChange*
388G4WeightCutOffProcess::AlongStepDoIt(const G4Track& track, const G4Step&)
389{
390 // Dummy ParticleChange ie: does nothing
391 // Expecting G4Transportation to move the track
392 pParticleChange->Initialize(track);
393 return pParticleChange;
394
395 // return 0;
396}
397
398void G4WeightCutOffProcess::CopyStep(const G4Step & step)
399{
400 fGhostStep->SetTrack(step.GetTrack());
401 fGhostStep->SetStepLength(step.GetStepLength());
402 fGhostStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
403 fGhostStep->SetControlFlag(step.GetControlFlag());
404
405 *fGhostPreStepPoint = *(step.GetPreStepPoint());
406 *fGhostPostStepPoint = *(step.GetPostStepPoint());
407
408//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
409// Set StepStatus for ghost world
410//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
411 if(fOnBoundary)
412 { fGhostPostStepPoint->SetStepStatus(fGeomBoundary); }
413 else if(fGhostPostStepPoint->GetStepStatus()==fGeomBoundary)
414 { fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc); }
415//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
416}
Note: See TracBrowser for help on using the repository browser.