Geant4  10.00.p03
G4WeightCutOffProcess.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: G4WeightCutOffProcess.cc 77999 2013-12-02 08:18:01Z gcosmo $
28 //
29 // ----------------------------------------------------------------------
30 // GEANT 4 class source file
31 //
32 // G4WeightCutOffProcess.cc
33 //
34 // ----------------------------------------------------------------------
35 
36 #include "G4WeightCutOffProcess.hh"
37 //#include "G4VScorer.hh"
38 #include "G4GeometryCellStep.hh"
39 //#include "G4GCellFinder.hh"
40 #include "G4TouchableHandle.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 
56  G4double wlimit,
57  G4double isource,
58  G4VIStore *istore,
59  // const G4VGCellFinder &aGCellFinder,
60  const G4String &aName, G4bool para)
61  : G4VProcess(aName),
62  fParticleChange(new G4ParticleChange),
63  fWeightSurvival(wsurvival),
64  fWeightLimit(wlimit),
65  fSourceImportance(isource),
66  fIStore(istore),
67  // fGCellFinder(aGCellFinder),
68  fGhostWorldName("NoParallelWorld"), fGhostWorld(0),
69  fGhostNavigator(0), fNavigatorID(-1), fFieldTrack('0')
70 {
71  if (!fParticleChange)
72  {
73  G4Exception("G4WeightCutOffProcess::G4WeightCutOffProcess()",
74  "FatalError", FatalException,
75  "Failed to allocate G4ParticleChange !");
76  }
77 
79 
80  fGhostStep = new G4Step();
83 
86 
87  if (verboseLevel>0)
88  {
89  G4cout << GetProcessName() << " is created " << G4endl;
90  }
91 
92  paraflag = para;
93 
94 }
95 
97 {
98  delete fParticleChange;
99  // delete fGhostStep;
100 }
101 
102 
103 //------------------------------------------------------
104 //
105 // SetParallelWorld
106 //
107 //------------------------------------------------------
109 SetParallelWorld(G4String parallelWorldName)
110 {
111 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
112 // Get pointers of the parallel world and its navigator
113 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
114  fGhostWorldName = parallelWorldName;
117 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
118 }
119 
122 {
123 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
124 // Get pointer of navigator
125 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
126  fGhostWorldName = parallelWorld->GetName();
127  fGhostWorld = parallelWorld;
129 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
130 }
131 
132 //------------------------------------------------------
133 //
134 // StartTracking
135 //
136 //------------------------------------------------------
138 {
139 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
140 // Activate navigator and get the navigator ID
141 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
142 // G4cout << " G4ParallelWorldScoringProcess::StartTracking" << G4endl;
143 
144  if(paraflag) {
145  if(fGhostNavigator)
147  else
148  {
149  G4Exception("G4WeightCutOffProcess::StartTracking",
150  "ProcParaWorld000",FatalException,
151  "G4WeightCutOffProcess is used for tracking without having a parallel world assigned");
152  }
153 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
154 
155 // G4cout << "G4ParallelWorldScoringProcess::StartTracking <<<<<<<<<<<<<<<<<< " << G4endl;
156 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
157 // Let PathFinder initialize
158 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
160 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
161 
162 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
163 // Setup initial touchables for the first step
164 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
169 
170  // Initialize
171  fGhostSafety = -1.;
172  fOnBoundary = false;
173  }
174 }
175 
176 
180 {
181 // *condition = Forced;
182 // return kInfinity;
183 
184 // *condition = StronglyForced;
185  *condition = Forced;
186  return DBL_MAX;
187 }
188 
191  const G4Step &aStep)
192 {
193  fParticleChange->Initialize(aTrack);
194 
195  if(paraflag) {
197  //xbug? fOnBoundary = false;
198  CopyStep(aStep);
199 
200  if(fOnBoundary)
201  {
202 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
203 // Locate the point and get new touchable
204 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
205  //?? fPathFinder->Locate(step.GetPostStepPoint()->GetPosition(),
206  //?? step.GetPostStepPoint()->GetMomentumDirection());
208 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
209  }
210  else
211  {
212  // Do I need this ??????????????????????????????????????????????????????????
213  // fGhostNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
214  // ?????????????????????????????????????????????????????????????????????????
215 
216  // fPathFinder->ReLocate(track.GetPosition());
217 
218  // reuse the touchable
220  }
221 
224 
225  }
226 
227  if(paraflag) {
230 
231 
232  // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(aStep);
233  // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(fGhostStep);
235  if (fIStore)
236  {
237  G4double i = fIStore->GetImportance(postCell);
238  if (i>0)
239  {
240  R/=i;
241  }
242  }
243  G4double w = aTrack.GetWeight();
244  if (w<R*fWeightLimit)
245  {
246  G4double ws = fWeightSurvival*R;
247  G4double p = w/(ws);
248  if (G4UniformRand()<p)
249  {
251  }
252  else
253  {
255  }
256  }
257  } else {
258 
259  G4GeometryCell postCell(*(aStep.GetPostStepPoint()->GetPhysicalVolume()),
261 
262  // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(aStep);
263  // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(fGhostStep);
265  if (fIStore)
266  {
267  G4double i = fIStore->GetImportance(postCell);
268  if (i>0)
269  {
270  R/=i;
271  }
272  }
273  G4double w = aTrack.GetWeight();
274  if (w<R*fWeightLimit)
275  {
276  G4double ws = fWeightSurvival*R;
277  G4double p = w/(ws);
278  if (G4UniformRand()<p)
279  {
281  }
282  else
283  {
285  }
286  }
287  }
288 
289  return fParticleChange;
290 
291 }
292 
294 {
295  return theProcessName;
296 }
297 
300  const G4Track& track, G4double previousStepSize, G4double currentMinimumStep,
301  G4double& proposedSafety, G4GPILSelection* selection)
302 {
303  if(paraflag) {
304  static G4ThreadLocal G4FieldTrack *endTrack_G4MT_TLS_ = 0 ; if (!endTrack_G4MT_TLS_) endTrack_G4MT_TLS_ = new G4FieldTrack ('0') ; G4FieldTrack &endTrack = *endTrack_G4MT_TLS_;
305  static G4ThreadLocal ELimited *eLimited_G4MT_TLS_ = 0 ; if (!eLimited_G4MT_TLS_) eLimited_G4MT_TLS_ = new ELimited ; ELimited &eLimited = *eLimited_G4MT_TLS_;
306 
307  *selection = NotCandidateForSelection;
308  G4double returnedStep = DBL_MAX;
309 
310  if (previousStepSize > 0.)
311  { fGhostSafety -= previousStepSize; }
312  // else
313  // { fGhostSafety = -1.; }
314  if (fGhostSafety < 0.) fGhostSafety = 0.0;
315 
316  // ------------------------------------------
317  // Determination of the proposed STEP LENGTH:
318  // ------------------------------------------
319  if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.)
320  {
321  // I have no chance to limit
322  returnedStep = currentMinimumStep;
323  fOnBoundary = false;
324  proposedSafety = fGhostSafety - currentMinimumStep;
325  }
326  else // (currentMinimumStep > fGhostSafety: I may limit the Step)
327  {
329  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
330  // ComputeStep
331  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
332  returnedStep
333  = fPathFinder->ComputeStep(fFieldTrack,currentMinimumStep,fNavigatorID,
334  track.GetCurrentStepNumber(),fGhostSafety,eLimited,
335  endTrack,track.GetVolume());
336  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
337  if(eLimited == kDoNot)
338  {
339  // Track is not on the boundary
340  fOnBoundary = false;
341  fGhostSafety = fGhostNavigator->ComputeSafety(endTrack.GetPosition());
342  }
343  else
344  {
345  // Track is on the boundary
346  fOnBoundary = true;
347  proposedSafety = fGhostSafety;
348  }
349  //xbug? proposedSafety = fGhostSafety;
350  if(eLimited == kUnique || eLimited == kSharedOther) {
351  *selection = CandidateForSelection;
352  }else if (eLimited == kSharedTransport) {
353  returnedStep *= (1.0 + 1.0e-9);
354  // Expand to disable its selection in Step Manager comparison
355  }
356 
357  }
358 
359  // ----------------------------------------------
360  // Returns the fGhostSafety as the proposedSafety
361  // The SteppingManager will take care of keeping
362  // the smallest one.
363  // ----------------------------------------------
364  return returnedStep;
365 
366  } else {
367  return DBL_MAX;
368  //not sensible! return -1.0;
369  }
370 
371 }
372 
373 
377 {
378  return -1.0;
379 }
380 
383 {
384  return 0;
385 }
386 
389 {
390  // Dummy ParticleChange ie: does nothing
391  // Expecting G4Transportation to move the track
392  pParticleChange->Initialize(track);
393  return pParticleChange;
394 
395  // return 0;
396 }
397 
399 {
400  fGhostStep->SetTrack(step.GetTrack());
404 
405  *fGhostPreStepPoint = *(step.GetPreStepPoint());
407 
408 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
409 // Set StepStatus for ghost world
410 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
411  if(fOnBoundary)
415 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
416 }
void PrepareNewTrack(const G4ThreeVector &position, const G4ThreeVector &direction, G4VPhysicalVolume *massStartVol=0)
static G4PathFinder * GetInstance()
Definition: G4PathFinder.cc:57
G4double condition(const G4ErrorSymMatrix &m)
G4VPhysicalVolume * fGhostWorld
void SetStepLength(G4double value)
virtual void Initialize(const G4Track &)
void SetTrack(G4Track *value)
G4int verboseLevel
Definition: G4VProcess.hh:368
G4double GetStepLength() const
G4StepStatus GetStepStatus() const
ELimited
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
const G4ThreeVector & GetPosition() const
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
G4TouchableHandle CreateTouchableHandle(G4int navId) const
G4SteppingControl GetControlFlag() const
const G4VTouchable * GetTouchable() const
#define G4ThreadLocal
Definition: tls.hh:52
void SetStepStatus(const G4StepStatus aValue)
G4TransportationManager * fTransportationManager
virtual G4VParticleChange * AtRestDoIt(const G4Track &, const G4Step &)
const G4String & GetName() const
G4StepPoint * GetPreStepPoint() const
void ProposeWeight(G4double finalWeight)
G4WeightCutOffProcess(G4double wsurvival, G4double wlimit, G4double isource, G4VIStore *istore, const G4String &aName="WeightCutOffProcess", G4bool para=false)
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4int GetCurrentStepNumber() const
G4VPhysicalVolume * GetPhysicalVolume() const
G4TouchableHandle fOldGhostTouchable
const G4String & GetName() const
void SetControlFlag(G4SteppingControl StepControlFlag)
bool G4bool
Definition: G4Types.hh:79
G4ParticleChange * fParticleChange
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 Initialize(const G4Track &)
G4VPhysicalVolume * GetParallelWorld(const G4String &worldName)
G4int ActivateNavigator(G4Navigator *aNavigator)
const G4ThreeVector & GetMomentumDirection() const
virtual G4VParticleChange * AlongStepDoIt(const G4Track &, const G4Step &)
G4Navigator * GetNavigator(const G4String &worldName)
virtual G4double AtRestGetPhysicalInteractionLength(const G4Track &, G4ForceCondition *)
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
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
#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
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4TouchableHandle fNewGhostTouchable
#define DBL_MAX
Definition: templates.hh:83
void SetParallelWorld(G4String parallelWorldName)
const G4TouchableHandle & GetTouchableHandle() const
static void Update(G4FieldTrack *, const G4Track *)
G4GPILSelection
void CopyStep(const G4Step &step)