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