Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ImportanceProcess.cc
Go to the documentation of this file.
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$
28 //
29 // ----------------------------------------------------------------------
30 // GEANT 4 class source file
31 //
32 // G4ImportanceProcess.cc
33 //
34 // ----------------------------------------------------------------------
35 
36 #include "G4ImportanceProcess.hh"
37 #include "G4VImportanceAlgorithm.hh"
38 #include "G4GeometryCell.hh"
39 #include "G4SamplingPostStepAction.hh"
40 #include "G4VTrackTerminator.hh"
41 #include "G4VIStore.hh"
42 
43 #include "G4Step.hh"
44 #include "G4Navigator.hh"
45 #include "G4VTouchable.hh"
46 #include "G4VPhysicalVolume.hh"
47 #include "G4ParticleChange.hh"
48 #include "G4PathFinder.hh"
49 #include "G4TransportationManager.hh"
50 #include "G4StepPoint.hh"
51 #include "G4FieldTrackUpdator.hh"
52 
53 
54 G4ImportanceProcess::
55 G4ImportanceProcess(const G4VImportanceAlgorithm &aImportanceAlgorithm,
56  const G4VIStore &aIstore,
57  const G4VTrackTerminator *TrackTerminator,
58  const G4String &aName, G4bool para)
59  : G4VProcess(aName),
60  fParticleChange(new G4ParticleChange),
61  fImportanceAlgorithm(aImportanceAlgorithm),
62  fIStore(aIstore),
63  fPostStepAction(0),
64  fGhostNavigator(0), fNavigatorID(-1), fFieldTrack('0'),
65  paraflag(para)
66 
67 {
68  if (TrackTerminator)
69  {
70  fPostStepAction = new G4SamplingPostStepAction(*TrackTerminator);
71  }
72  else
73  {
74  fPostStepAction = new G4SamplingPostStepAction(*this);
75  }
76  if (!fParticleChange)
77  {
78  G4Exception("G4ImportanceProcess::G4ImportanceProcess()",
79  "FatalError", FatalException,
80  "Failed allocation of G4ParticleChange !");
81  }
82  G4VProcess::pParticleChange = fParticleChange;
83 
84 
85  fGhostStep = new G4Step();
86  fGhostPreStepPoint = fGhostStep->GetPreStepPoint();
87  fGhostPostStepPoint = fGhostStep->GetPostStepPoint();
88 
89  fTransportationManager = G4TransportationManager::GetTransportationManager();
90  fPathFinder = G4PathFinder::GetInstance();
91 
92  if (verboseLevel>0)
93  {
94  G4cout << GetProcessName() << " is created " << G4endl;
95  }
96 
97  G4cout << " importance process paraflag is: " << paraflag << G4endl;
98 
99 }
100 
101 G4ImportanceProcess::~G4ImportanceProcess()
102 {
103 
104  delete fPostStepAction;
105  delete fParticleChange;
106  delete fGhostStep;
107 
108 }
109 
110 
111 
112 //------------------------------------------------------
113 //
114 // SetParallelWorld
115 //
116 //------------------------------------------------------
117 void G4ImportanceProcess::
118 SetParallelWorld(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 
129 void G4ImportanceProcess::
130 SetParallelWorld(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 //------------------------------------------------------
146 void G4ImportanceProcess::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("G4ImportanceProcess::StartTracking",
159  "ProcParaWorld000",FatalException,
160  "G4ImportanceProcess 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 
186 
187 G4double G4ImportanceProcess::
188 PostStepGetPhysicalInteractionLength(const G4Track& ,
189  G4double ,
190  G4ForceCondition* condition)
191 {
192 // *condition = Forced;
193 // return kInfinity;
194 
195 // *condition = StronglyForced;
196  *condition = Forced;
197  return DBL_MAX;
198 }
199 
200 G4VParticleChange *
201 G4ImportanceProcess::PostStepDoIt(const G4Track &aTrack,
202  const G4Step &aStep)
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  //AH G4cout << " on boundary " << G4endl;
220 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
221  }
222  else
223  {
224  // Do I need this ??????????????????????????????????????????????????????????
225  // fGhostNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
226  // ?????????????????????????????????????????????????????????????????????????
227 
228  // fPathFinder->ReLocate(track.GetPosition());
229 
230  // reuse the touchable
231  fNewGhostTouchable = fOldGhostTouchable;
232  //AH G4cout << " NOT on boundary " << G4endl;
233  }
234 
235  fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
236  fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
237 
238 // if ( (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
239 // && (aStep.GetStepLength() > kCarTolerance) )
240 // {
241 // if (aTrack.GetTrackStatus()==fStopAndKill)
242 // {
243 // G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
244 // << " StopAndKill track." << G4endl;
245 // }
246 
247 // G4StepPoint *prepoint = aStep.GetPreStepPoint();
248 // G4StepPoint *postpoint = aStep.GetPostStepPoint();
249 // G4GeometryCell prekey(*(prepoint->GetPhysicalVolume()),
250 // prepoint->GetTouchable()->GetReplicaNumber());
251 // G4GeometryCell postkey(*(postpoint->GetPhysicalVolume()),
252 // postpoint->GetTouchable()->GetReplicaNumber());
253 
254 
255  if ( (fGhostPostStepPoint->GetStepStatus() == fGeomBoundary)
256  && (aStep.GetStepLength() > kCarTolerance) )
257  {
258  if (aTrack.GetTrackStatus()==fStopAndKill)
259  {
260  G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
261  << " StopAndKill track. on boundary" << G4endl;
262  }
263 
264  G4GeometryCell prekey(*(fGhostPreStepPoint->GetPhysicalVolume()),
265  fGhostPreStepPoint->GetTouchable()->GetReplicaNumber());
266  G4GeometryCell postkey(*(fGhostPostStepPoint->GetPhysicalVolume()),
267  fGhostPostStepPoint->GetTouchable()->GetReplicaNumber());
268 
269  //AH
270  /*
271  G4cout << G4endl;
272  G4cout << G4endl;
273  G4cout << " inside parallel importance process " << aTrack.GetCurrentStepNumber() << G4endl;
274  G4cout << G4endl;
275  G4cout << G4endl;
276  G4cout << " prekey: " << fGhostPreStepPoint->GetPhysicalVolume()->GetName() << " replica:"
277  << fGhostPreStepPoint->GetTouchable()->GetReplicaNumber() << G4endl;
278  G4cout << " prekey ISTORE: " << fIStore.GetImportance(prekey) << G4endl;
279  G4cout << " postkey: " << G4endl;
280  G4cout << " postkey ISTORE: " << fIStore.GetImportance(postkey) << G4endl;
281  */
282  //AH
283  G4Nsplit_Weight nw = fImportanceAlgorithm.
284  Calculate(fIStore.GetImportance(prekey),
285  fIStore.GetImportance(postkey),
286  aTrack.GetWeight());
287  //AH
288  /*
289  G4cout << " prekey weight: " << fIStore.GetImportance(prekey)
290  << " postkey weight: " << fIStore.GetImportance(postkey)
291  << " split weight: " << nw << G4endl;
292  G4cout << " before poststepaction " << G4endl;
293  */
294  //AH
295  fPostStepAction->DoIt(aTrack, fParticleChange, nw);
296  //AH G4cout << " after post step do it " << G4endl;
297  }
298  } else {
299  if ( (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
300  && (aStep.GetStepLength() > kCarTolerance) )
301  {
302  //AH G4cout << " inside non-parallel importance process " << G4endl;
303  if (aTrack.GetTrackStatus()==fStopAndKill)
304  {
305  G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
306  << " StopAndKill track. on boundary non-parallel" << G4endl;
307  }
308 
309  G4StepPoint *prepoint = aStep.GetPreStepPoint();
310  G4StepPoint *postpoint = aStep.GetPostStepPoint();
311 
312  G4GeometryCell prekey(*(prepoint->GetPhysicalVolume()),
313  prepoint->GetTouchable()->GetReplicaNumber());
314  G4GeometryCell postkey(*(postpoint->GetPhysicalVolume()),
315  postpoint->GetTouchable()->GetReplicaNumber());
316 
317  G4Nsplit_Weight nw = fImportanceAlgorithm.
318  Calculate(fIStore.GetImportance(prekey),
319  fIStore.GetImportance(postkey),
320  aTrack.GetWeight());
321  //AH
322  /*
323  G4cout << " prekey weight: " << fIStore.GetImportance(prekey)
324  << " postkey weight: " << fIStore.GetImportance(postkey)
325  << " split weight: " << nw << G4endl;
326  G4cout << " before poststepaction 2 " << G4endl;
327  */
328  //AH
329  fPostStepAction->DoIt(aTrack, fParticleChange, nw);
330  //AH G4cout << " after poststepaction 2 " << G4endl;
331  }
332  }
333  return fParticleChange;
334 }
335 
336 void G4ImportanceProcess::KillTrack() const
337 {
338  fParticleChange->ProposeTrackStatus(fStopAndKill);
339 }
340 
341 const G4String &G4ImportanceProcess::GetName() const
342 {
343  return theProcessName;
344 }
345 
346 G4double G4ImportanceProcess::
347 AlongStepGetPhysicalInteractionLength(
348  const G4Track& track, G4double previousStepSize, G4double currentMinimumStep,
349  G4double& proposedSafety, G4GPILSelection* selection)
350 {
351 
352  //AH G4cout << " along step physical interaction length " << G4endl;
353 
354  if(paraflag) {
355  static G4FieldTrack endTrack('0');
356  static ELimited eLimited;
357 
358  *selection = NotCandidateForSelection;
359  G4double returnedStep = DBL_MAX;
360 
361  if (previousStepSize > 0.)
362  { fGhostSafety -= previousStepSize; }
363  // else
364  // { fGhostSafety = -1.; }
365  if (fGhostSafety < 0.) fGhostSafety = 0.0;
366 
367  // ------------------------------------------
368  // Determination of the proposed STEP LENGTH:
369  // ------------------------------------------
370  if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.)
371  {
372  // I have no chance to limit
373  returnedStep = currentMinimumStep;
374  fOnBoundary = false;
375  proposedSafety = fGhostSafety - currentMinimumStep;
376  //AH G4cout << " step not limited, why? " << G4endl;
377  }
378  else // (currentMinimumStep > fGhostSafety: I may limit the Step)
379  {
380  G4FieldTrackUpdator::Update(&fFieldTrack,&track);
381  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
382  // ComputeStep
383  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
384  returnedStep
385  = fPathFinder->ComputeStep(fFieldTrack,currentMinimumStep,fNavigatorID,
386  track.GetCurrentStepNumber(),fGhostSafety,eLimited,
387  endTrack,track.GetVolume());
388  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
389  if(eLimited == kDoNot)
390  {
391  //AH G4cout << " computestep came back with not-boundary " << G4endl;
392  // Track is not on the boundary
393  fOnBoundary = false;
394  fGhostSafety = fGhostNavigator->ComputeSafety(endTrack.GetPosition());
395  }
396  else
397  {
398  // Track is on the boundary
399  //AH G4cout << " FOUND A BOUNDARY ! " << track.GetCurrentStepNumber() << G4endl;
400  fOnBoundary = true;
401  // proposedSafety = fGhostSafety;
402  }
403  proposedSafety = fGhostSafety;
404  if(eLimited == kUnique || eLimited == kSharedOther) {
405  *selection = CandidateForSelection;
406  }else if (eLimited == kSharedTransport) {
407  returnedStep *= (1.0 + 1.0e-9);
408  // Expand to disable its selection in Step Manager comparison
409  }
410 
411  }
412 
413  // ----------------------------------------------
414  // Returns the fGhostSafety as the proposedSafety
415  // The SteppingManager will take care of keeping
416  // the smallest one.
417  // ----------------------------------------------
418  return returnedStep;
419 
420  } else {
421 
422  return DBL_MAX;
423 
424  }
425 
426 }
427 
428 G4double G4ImportanceProcess::
429 AtRestGetPhysicalInteractionLength(const G4Track& ,
430  G4ForceCondition*)
431 {
432  return -1.0;
433 }
434 
435 G4VParticleChange* G4ImportanceProcess::
436 AtRestDoIt(const G4Track&, const G4Step&)
437 {
438  return 0;
439 }
440 
441 G4VParticleChange* G4ImportanceProcess::
442 AlongStepDoIt(const G4Track& aTrack, const G4Step& )
443 {
444  // Dummy ParticleChange ie: does nothing
445  // Expecting G4Transportation to move the track
446  //AH G4cout << " along step do it " << G4endl;
447  pParticleChange->Initialize(aTrack);
448  return pParticleChange;
449  // return 0;
450 }
451 
452 void G4ImportanceProcess::CopyStep(const G4Step & step)
453 {
454  fGhostStep->SetTrack(step.GetTrack());
455  fGhostStep->SetStepLength(step.GetStepLength());
456  fGhostStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
457  fGhostStep->SetControlFlag(step.GetControlFlag());
458 
459  *fGhostPreStepPoint = *(step.GetPreStepPoint());
460  *fGhostPostStepPoint = *(step.GetPostStepPoint());
461 
462 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
463 // Set StepStatus for ghost world
464 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
465  if(fOnBoundary)
466  { fGhostPostStepPoint->SetStepStatus(fGeomBoundary); }
467  else if(fGhostPostStepPoint->GetStepStatus()==fGeomBoundary)
468  { fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc); }
469 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
470 }
Generated on Sat May 25 2013 14:33:29 for Geant4 by   doxygen 1.8.4