Geant4  10.03.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: G4ImportanceProcess.cc 88804 2015-03-10 17:12:21Z gcosmo $
28 //
29 // ----------------------------------------------------------------------
30 // GEANT 4 class source file
31 //
32 // G4ImportanceProcess.cc
33 //
34 // ----------------------------------------------------------------------
35 
36 #include "G4ImportanceProcess.hh"
38 #include "G4GeometryCell.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"
50 #include "G4StepPoint.hh"
51 #include "G4FieldTrackUpdator.hh"
52 
53 
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  fGhostWorldName("NoParallelWorld"),fGhostWorld(0),
65  fGhostNavigator(0), fNavigatorID(-1), fFieldTrack('0'),
66  fParaflag(para), fEndTrack('0'), feLimited(kDoNot)
67 {
68  G4cout << G4endl << G4endl << G4endl;
69  G4cout << "G4ImportanceProcess:: Creating " << G4endl;
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("G4ImportanceProcess::G4ImportanceProcess()",
81  "FatalError", FatalException,
82  "Failed allocation of G4ParticleChange !");
83  }
84  G4VProcess::pParticleChange = fParticleChange;
85 
86 
87  fGhostStep = new G4Step();
88  fGhostPreStepPoint = fGhostStep->GetPreStepPoint();
89  fGhostPostStepPoint = fGhostStep->GetPostStepPoint();
90 
91  fTransportationManager = G4TransportationManager::GetTransportationManager();
92  fPathFinder = G4PathFinder::GetInstance();
93 
94  if (verboseLevel>0)
95  {
96  G4cout << GetProcessName() << " is created " << G4endl;
97  }
98 
99  G4cout << "G4ImportanceProcess:: importance process paraflag is: " << fParaflag << G4endl;
100 
101 }
102 
104 {
105 
106  delete fPostStepAction;
107  delete fParticleChange;
108  // delete fGhostStep;
109  // delete fGhostWorld;
110  // delete fGhostNavigator;
111 
112 }
113 
114 
115 
116 //------------------------------------------------------
117 //
118 // SetParallelWorld
119 //
120 //------------------------------------------------------
122 SetParallelWorld(const G4String &parallelWorldName)
123 {
124  G4cout << G4endl << G4endl << G4endl;
125  G4cout << "G4ImportanceProcess:: SetParallelWorld name = " << parallelWorldName << G4endl;
126 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
127 // Get pointers of the parallel world and its navigator
128 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
129  fGhostWorldName = parallelWorldName;
130  fGhostWorld = fTransportationManager->GetParallelWorld(fGhostWorldName);
131  fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
132 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
133 }
134 
135 // void G4ImportanceProcess::
136 // SetParallelWorld(const G4VPhysicalVolume* parallelWorld)
137 // {
138 // //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
139 // // Get pointer of navigator
140 // //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
141 // // G4cout << " G4ImportanceProcess:: Got here 1 " << G4endl;
142 // // fGhostWorldName = parallelWorld->GetName();
143 // // G4cout << " G4ImportanceProcess:: Got here 2 ghostName:" << fGhostWorldName << G4endl;
144 // fGhostWorld = parallelWorld;
145 // G4cout << " G4ImportanceProcess:: Got here 3 " << G4endl;
146 // fGhostNavigator = fTransportationManager->GetNavigator(parallelWorld);
147 // G4cout << " G4ImportanceProcess:: Got here 4 " << G4endl;
148 // //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
149 // }
150 
151 //------------------------------------------------------
152 //
153 // StartTracking
154 //
155 //------------------------------------------------------
157 {
158 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
159 // Activate navigator and get the navigator ID
160 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
161 // G4cout << " G4ParallelWorldScoringProcess::StartTracking" << G4endl;
162 
163  if(fParaflag) {
164  if(fGhostNavigator)
165  { fNavigatorID = fTransportationManager->ActivateNavigator(fGhostNavigator); }
166  else
167  {
168  G4Exception("G4ImportanceProcess::StartTracking",
169  "ProcParaWorld000",FatalException,
170  "G4ImportanceProcess is used for tracking without having a parallel world assigned");
171  }
172 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
173 
174 // G4cout << "G4ParallelWorldScoringProcess::StartTracking <<<<<<<<<<<<<<<<<< " << G4endl;
175 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
176 // Let PathFinder initialize
177 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
178  fPathFinder->PrepareNewTrack(trk->GetPosition(),trk->GetMomentumDirection());
179 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
180 
181 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
182 // Setup initial touchables for the first step
183 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
184  fOldGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
185  fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
186  fNewGhostTouchable = fOldGhostTouchable;
187  fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
188 
189  // Initialize
190  fGhostSafety = -1.;
191  fOnBoundary = false;
192  }
193 
194 }
195 
196 
199  G4double ,
201 {
202 // *condition = Forced;
203 // return kInfinity;
204 
205 // *condition = StronglyForced;
206  *condition = Forced;
207  return DBL_MAX;
208 }
209 
212  const G4Step &aStep)
213 {
214  fParticleChange->Initialize(aTrack);
215 
216  if(fParaflag) {
217  fOldGhostTouchable = fGhostPostStepPoint->GetTouchableHandle();
218  //xbug? fOnBoundary = false;
219  CopyStep(aStep);
220 
221  if(fOnBoundary)
222  {
223 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
224 // Locate the point and get new touchable
225 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
226  //?? fPathFinder->Locate(step.GetPostStepPoint()->GetPosition(),
227  //?? step.GetPostStepPoint()->GetMomentumDirection());
228  fNewGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
229  //AH G4cout << " on boundary " << G4endl;
230 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
231  }
232  else
233  {
234  // Do I need this ??????????????????????????????????????????????????????????
235  // fGhostNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
236  // ?????????????????????????????????????????????????????????????????????????
237 
238  // fPathFinder->ReLocate(track.GetPosition());
239 
240  // reuse the touchable
241  fNewGhostTouchable = fOldGhostTouchable;
242  //AH G4cout << " NOT on boundary " << G4endl;
243  }
244 
245  fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
246  fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
247 
248 // if ( (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
249 // && (aStep.GetStepLength() > kCarTolerance) )
250 // {
251 // if (aTrack.GetTrackStatus()==fStopAndKill)
252 // {
253 // G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
254 // << " StopAndKill track." << G4endl;
255 // }
256 
257 // G4StepPoint *prepoint = aStep.GetPreStepPoint();
258 // G4StepPoint *postpoint = aStep.GetPostStepPoint();
259 // G4GeometryCell prekey(*(prepoint->GetPhysicalVolume()),
260 // prepoint->GetTouchable()->GetReplicaNumber());
261 // G4GeometryCell postkey(*(postpoint->GetPhysicalVolume()),
262 // postpoint->GetTouchable()->GetReplicaNumber());
263 
264 
265  if ( (fGhostPostStepPoint->GetStepStatus() == fGeomBoundary)
266  && (aStep.GetStepLength() > kCarTolerance) )
267  {
268  if (aTrack.GetTrackStatus()==fStopAndKill)
269  {
270  G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
271  << " StopAndKill track. on boundary" << G4endl;
272  }
273 
274  G4GeometryCell prekey(*(fGhostPreStepPoint->GetPhysicalVolume()),
275  fGhostPreStepPoint->GetTouchable()->GetReplicaNumber());
276  G4GeometryCell postkey(*(fGhostPostStepPoint->GetPhysicalVolume()),
277  fGhostPostStepPoint->GetTouchable()->GetReplicaNumber());
278 
279  //AH
280  /*
281  G4cout << G4endl;
282  G4cout << G4endl;
283  G4cout << " inside parallel importance process " << aTrack.GetCurrentStepNumber() << G4endl;
284  G4cout << G4endl;
285  G4cout << G4endl;
286  G4cout << " prekey: " << fGhostPreStepPoint->GetPhysicalVolume()->GetName() << " replica:"
287  << fGhostPreStepPoint->GetTouchable()->GetReplicaNumber() << G4endl;
288  G4cout << " prekey ISTORE: " << fIStore.GetImportance(prekey) << G4endl;
289  G4cout << " postkey: " << G4endl;
290  G4cout << " postkey ISTORE: " << fIStore.GetImportance(postkey) << G4endl;
291  */
292  //AH
293  G4Nsplit_Weight nw = fImportanceAlgorithm.
294  Calculate(fIStore.GetImportance(prekey),
295  fIStore.GetImportance(postkey),
296  aTrack.GetWeight());
297  //AH
298  /*
299  G4cout << " prekey weight: " << fIStore.GetImportance(prekey)
300  << " postkey weight: " << fIStore.GetImportance(postkey)
301  << " split weight: " << nw << G4endl;
302  G4cout << " before poststepaction " << G4endl;
303  */
304  //AH
305  fPostStepAction->DoIt(aTrack, fParticleChange, nw);
306  //AH G4cout << " after post step do it " << G4endl;
307  }
308  } else {
309  if ( (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
310  && (aStep.GetStepLength() > kCarTolerance) )
311  {
312  //AH G4cout << " inside non-parallel importance process " << G4endl;
313  if (aTrack.GetTrackStatus()==fStopAndKill)
314  {
315  G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
316  << " StopAndKill track. on boundary non-parallel" << G4endl;
317  }
318 
319  G4StepPoint *prepoint = aStep.GetPreStepPoint();
320  G4StepPoint *postpoint = aStep.GetPostStepPoint();
321 
322  G4GeometryCell prekey(*(prepoint->GetPhysicalVolume()),
323  prepoint->GetTouchable()->GetReplicaNumber());
324  G4GeometryCell postkey(*(postpoint->GetPhysicalVolume()),
325  postpoint->GetTouchable()->GetReplicaNumber());
326 
327  G4Nsplit_Weight nw = fImportanceAlgorithm.
328  Calculate(fIStore.GetImportance(prekey),
329  fIStore.GetImportance(postkey),
330  aTrack.GetWeight());
331  //AH
332  /*
333  G4cout << " prekey weight: " << fIStore.GetImportance(prekey)
334  << " postkey weight: " << fIStore.GetImportance(postkey)
335  << " split weight: " << nw << G4endl;
336  G4cout << " before poststepaction 2 " << G4endl;
337  */
338  //AH
339  fPostStepAction->DoIt(aTrack, fParticleChange, nw);
340  //AH G4cout << " after poststepaction 2 " << G4endl;
341  }
342  }
343  return fParticleChange;
344 }
345 
347 {
348  fParticleChange->ProposeTrackStatus(fStopAndKill);
349 }
350 
352 {
353  return theProcessName;
354 }
355 
358  const G4Track& track, G4double previousStepSize, G4double currentMinimumStep,
359  G4double& proposedSafety, G4GPILSelection* selection)
360 {
361  if(fParaflag) {
362  *selection = NotCandidateForSelection;
363  G4double returnedStep = DBL_MAX;
364 
365  if (previousStepSize > 0.)
366  { fGhostSafety -= previousStepSize; }
367  // else
368  // { fGhostSafety = -1.; }
369  if (fGhostSafety < 0.) fGhostSafety = 0.0;
370 
371  // ------------------------------------------
372  // Determination of the proposed STEP LENGTH:
373  // ------------------------------------------
374  if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.)
375  {
376  // I have no chance to limit
377  returnedStep = currentMinimumStep;
378  fOnBoundary = false;
379  proposedSafety = fGhostSafety - currentMinimumStep;
380  //AH G4cout << " step not limited, why? " << G4endl;
381  }
382  else // (currentMinimumStep > fGhostSafety: I may limit the Step)
383  {
384  G4FieldTrackUpdator::Update(&fFieldTrack,&track);
385  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
386  // ComputeStep
387  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
388  returnedStep
389  = fPathFinder->ComputeStep(fFieldTrack,currentMinimumStep,fNavigatorID,
390  track.GetCurrentStepNumber(),fGhostSafety,feLimited,
391  fEndTrack,track.GetVolume());
392  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
393  if(feLimited == kDoNot)
394  {
395  //AH G4cout << " computestep came back with not-boundary " << G4endl;
396  // Track is not on the boundary
397  fOnBoundary = false;
398  fGhostSafety = fGhostNavigator->ComputeSafety(fEndTrack.GetPosition());
399  }
400  else
401  {
402  // Track is on the boundary
403  //AH G4cout << " FOUND A BOUNDARY ! " << track.GetCurrentStepNumber() << G4endl;
404  fOnBoundary = true;
405  // proposedSafety = fGhostSafety;
406  }
407  proposedSafety = fGhostSafety;
408  if(feLimited == kUnique || feLimited == kSharedOther) {
409  *selection = CandidateForSelection;
410  }else if (feLimited == kSharedTransport) {
411  returnedStep *= (1.0 + 1.0e-9);
412  // Expand to disable its selection in Step Manager comparison
413  }
414 
415  }
416 
417  // ----------------------------------------------
418  // Returns the fGhostSafety as the proposedSafety
419  // The SteppingManager will take care of keeping
420  // the smallest one.
421  // ----------------------------------------------
422  return returnedStep;
423 
424  } else {
425 
426  return DBL_MAX;
427 
428  }
429 
430 }
431 
435 {
436  return -1.0;
437 }
438 
440 AtRestDoIt(const G4Track&, const G4Step&)
441 {
442  return 0;
443 }
444 
446 AlongStepDoIt(const G4Track& aTrack, const G4Step& )
447 {
448  // Dummy ParticleChange ie: does nothing
449  // Expecting G4Transportation to move the track
450  //AH G4cout << " along step do it " << G4endl;
451  pParticleChange->Initialize(aTrack);
452  return pParticleChange;
453  // return 0;
454 }
455 
456 void G4ImportanceProcess::CopyStep(const G4Step & step)
457 {
458  fGhostStep->SetTrack(step.GetTrack());
459  fGhostStep->SetStepLength(step.GetStepLength());
460  fGhostStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
461  fGhostStep->SetControlFlag(step.GetControlFlag());
462 
463  *fGhostPreStepPoint = *(step.GetPreStepPoint());
464  *fGhostPostStepPoint = *(step.GetPostStepPoint());
465 
466 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
467 // Set StepStatus for ghost world
468 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
469  if(fOnBoundary)
470  { fGhostPostStepPoint->SetStepStatus(fGeomBoundary); }
471  else if(fGhostPostStepPoint->GetStepStatus()==fGeomBoundary)
472  { fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc); }
473 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
474 }
void PrepareNewTrack(const G4ThreeVector &position, const G4ThreeVector &direction, G4VPhysicalVolume *massStartVol=0)
static G4PathFinder * GetInstance()
Definition: G4PathFinder.cc:57
G4double condition(const G4ErrorSymMatrix &m)
void SetStepLength(G4double value)
virtual void Initialize(const G4Track &)
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
void SetTrack(G4Track *value)
void SetParallelWorld(const G4String &parallelWorldName)
G4int verboseLevel
Definition: G4VProcess.hh:368
virtual G4VParticleChange * AlongStepDoIt(const G4Track &, const G4Step &)
G4double GetStepLength() const
G4StepStatus GetStepStatus() const
const G4ThreeVector & GetPosition() const
G4TrackStatus GetTrackStatus() const
virtual G4VParticleChange * AtRestDoIt(const G4Track &, const G4Step &)
G4TouchableHandle CreateTouchableHandle(G4int navId) const
G4SteppingControl GetControlFlag() const
const G4VTouchable * GetTouchable() const
void SetStepStatus(const G4StepStatus aValue)
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
G4StepPoint * GetPreStepPoint() const
G4GLOB_DLL std::ostream G4cout
G4int GetCurrentStepNumber() const
virtual const G4String & GetName() const
G4VPhysicalVolume * GetPhysicalVolume() const
void SetControlFlag(G4SteppingControl StepControlFlag)
bool G4bool
Definition: G4Types.hh:79
G4double ComputeStep(const G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4int navigatorId, G4int stepNo, G4double &pNewSafety, ELimited &limitedStep, G4FieldTrack &EndState, G4VPhysicalVolume *currentVolume)
Definition: G4Step.hh:76
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
virtual G4double GetImportance(const G4GeometryCell &gCell) const =0
G4double GetTotalEnergyDeposit() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4TransportationManager * GetTransportationManager()
virtual void KillTrack() const
void StartTracking(G4Track *)
virtual void Initialize(const G4Track &)
G4VPhysicalVolume * GetParallelWorld(const G4String &worldName)
G4int ActivateNavigator(G4Navigator *aNavigator)
void DoIt(const G4Track &aTrack, G4ParticleChange *aParticleChange, const G4Nsplit_Weight &nw)
const G4ThreeVector & GetMomentumDirection() const
G4Navigator * GetNavigator(const G4String &worldName)
G4StepPoint * GetPostStepPoint() const
virtual G4int GetReplicaNumber(G4int depth=0) const
Definition: G4VTouchable.cc:58
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
G4VPhysicalVolume * GetVolume() const
G4double GetWeight() const
G4ImportanceProcess(const G4VImportanceAlgorithm &aImportanceAlgorithm, const G4VIStore &aIstore, const G4VTrackTerminator *TrackTerminator, const G4String &aName="ImportanceProcess", G4bool para=false)
#define G4endl
Definition: G4ios.hh:61
G4String theProcessName
Definition: G4VProcess.hh:335
void SetTotalEnergyDeposit(G4double value)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
virtual G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4double pProposedMaxLength=DBL_MAX, const G4bool keepState=true)
G4ForceCondition
G4Track * GetTrack() const
virtual G4double AtRestGetPhysicalInteractionLength(const G4Track &, G4ForceCondition *)
void SetTouchableHandle(const G4TouchableHandle &apValue)
#define DBL_MAX
Definition: templates.hh:83
const G4TouchableHandle & GetTouchableHandle() const
static void Update(G4FieldTrack *, const G4Track *)
G4GPILSelection
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)