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