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