Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 88804 2015-03-10 17:12:21Z 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  fParaflag(para), fEndTrack('0'), feLimited(kDoNot)
71 {
72  if (!fParticleChange)
73  {
74  G4Exception("G4WeightCutOffProcess::G4WeightCutOffProcess()",
75  "FatalError", FatalException,
76  "Failed to allocate G4ParticleChange !");
77  }
78 
79  G4VProcess::pParticleChange = fParticleChange;
80 
81  fGhostStep = new G4Step();
82  fGhostPreStepPoint = fGhostStep->GetPreStepPoint();
83  fGhostPostStepPoint = fGhostStep->GetPostStepPoint();
84 
85  fTransportationManager = G4TransportationManager::GetTransportationManager();
86  fPathFinder = G4PathFinder::GetInstance();
87 
88  if (verboseLevel>0)
89  {
90  G4cout << GetProcessName() << " is created " << G4endl;
91  }
92 }
93 
95 {
96  delete fParticleChange;
97  // delete fGhostStep;
98 }
99 
100 
101 //------------------------------------------------------
102 //
103 // SetParallelWorld
104 //
105 //------------------------------------------------------
107 SetParallelWorld(const G4String &parallelWorldName)
108 {
109 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
110 // Get pointers of the parallel world and its navigator
111 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
112  fGhostWorldName = parallelWorldName;
113  fGhostWorld = fTransportationManager->GetParallelWorld(fGhostWorldName);
114  fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
115 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
116 }
117 
120 {
121 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
122 // Get pointer of navigator
123 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
124  fGhostWorldName = parallelWorld->GetName();
125  fGhostWorld = parallelWorld;
126  fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
127 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
128 }
129 
130 //------------------------------------------------------
131 //
132 // StartTracking
133 //
134 //------------------------------------------------------
136 {
137 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
138 // Activate navigator and get the navigator ID
139 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
140 // G4cout << " G4ParallelWorldScoringProcess::StartTracking" << G4endl;
141 
142  if(fParaflag) {
143  if(fGhostNavigator)
144  { fNavigatorID = fTransportationManager->ActivateNavigator(fGhostNavigator); }
145  else
146  {
147  G4Exception("G4WeightCutOffProcess::StartTracking",
148  "ProcParaWorld000",FatalException,
149  "G4WeightCutOffProcess is used for tracking without having a parallel world assigned");
150  }
151 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
152 
153 // G4cout << "G4ParallelWorldScoringProcess::StartTracking <<<<<<<<<<<<<<<<<< " << G4endl;
154 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
155 // Let PathFinder initialize
156 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
157  fPathFinder->PrepareNewTrack(trk->GetPosition(),trk->GetMomentumDirection());
158 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
159 
160 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
161 // Setup initial touchables for the first step
162 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
163  fOldGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
164  fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
165  fNewGhostTouchable = fOldGhostTouchable;
166  fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
167 
168  // Initialize
169  fGhostSafety = -1.;
170  fOnBoundary = false;
171  }
172 }
173 
174 
178 {
179 // *condition = Forced;
180 // return kInfinity;
181 
182 // *condition = StronglyForced;
183  *condition = Forced;
184  return DBL_MAX;
185 }
186 
189  const G4Step &aStep)
190 {
191  fParticleChange->Initialize(aTrack);
192 
193  if(fParaflag) {
194  fOldGhostTouchable = fGhostPostStepPoint->GetTouchableHandle();
195  //xbug? fOnBoundary = false;
196  CopyStep(aStep);
197 
198  if(fOnBoundary)
199  {
200 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
201 // Locate the point and get new touchable
202 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
203  //?? fPathFinder->Locate(step.GetPostStepPoint()->GetPosition(),
204  //?? step.GetPostStepPoint()->GetMomentumDirection());
205  fNewGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
206 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
207  }
208  else
209  {
210  // Do I need this ??????????????????????????????????????????????????????????
211  // fGhostNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
212  // ?????????????????????????????????????????????????????????????????????????
213 
214  // fPathFinder->ReLocate(track.GetPosition());
215 
216  // reuse the touchable
217  fNewGhostTouchable = fOldGhostTouchable;
218  }
219 
220  fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
221  fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
222 
223  }
224 
225  if(fParaflag) {
226  G4GeometryCell postCell(*(fGhostPostStepPoint->GetPhysicalVolume()),
227  fGhostPostStepPoint->GetTouchable()->GetReplicaNumber());
228 
229 
230  // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(aStep);
231  // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(fGhostStep);
232  G4double R = fSourceImportance;
233  if (fIStore)
234  {
235  G4double i = fIStore->GetImportance(postCell);
236  if (i>0)
237  {
238  R/=i;
239  }
240  }
241  G4double w = aTrack.GetWeight();
242  if (w<R*fWeightLimit)
243  {
244  G4double ws = fWeightSurvival*R;
245  G4double p = w/(ws);
246  if (G4UniformRand()<p)
247  {
248  fParticleChange->ProposeTrackStatus(fStopAndKill);
249  }
250  else
251  {
252  fParticleChange->ProposeWeight(ws);
253  }
254  }
255  } else {
256 
257  G4GeometryCell postCell(*(aStep.GetPostStepPoint()->GetPhysicalVolume()),
259 
260  // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(aStep);
261  // G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(fGhostStep);
262  G4double R = fSourceImportance;
263  if (fIStore)
264  {
265  G4double i = fIStore->GetImportance(postCell);
266  if (i>0)
267  {
268  R/=i;
269  }
270  }
271  G4double w = aTrack.GetWeight();
272  if (w<R*fWeightLimit)
273  {
274  G4double ws = fWeightSurvival*R;
275  G4double p = w/(ws);
276  if (G4UniformRand()<p)
277  {
278  fParticleChange->ProposeTrackStatus(fStopAndKill);
279  }
280  else
281  {
282  fParticleChange->ProposeWeight(ws);
283  }
284  }
285  }
286 
287  return fParticleChange;
288 
289 }
290 
292 {
293  return theProcessName;
294 }
295 
298  const G4Track& track, G4double previousStepSize, G4double currentMinimumStep,
299  G4double& proposedSafety, G4GPILSelection* selection)
300 {
301  if(fParaflag) {
302 
303  *selection = NotCandidateForSelection;
304  G4double returnedStep = DBL_MAX;
305 
306  if (previousStepSize > 0.)
307  { fGhostSafety -= previousStepSize; }
308  // else
309  // { fGhostSafety = -1.; }
310  if (fGhostSafety < 0.) fGhostSafety = 0.0;
311 
312  // ------------------------------------------
313  // Determination of the proposed STEP LENGTH:
314  // ------------------------------------------
315  if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.)
316  {
317  // I have no chance to limit
318  returnedStep = currentMinimumStep;
319  fOnBoundary = false;
320  proposedSafety = fGhostSafety - currentMinimumStep;
321  }
322  else // (currentMinimumStep > fGhostSafety: I may limit the Step)
323  {
324  G4FieldTrackUpdator::Update(&fFieldTrack,&track);
325  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
326  // ComputeStep
327  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
328  returnedStep
329  = fPathFinder->ComputeStep(fFieldTrack,currentMinimumStep,fNavigatorID,
330  track.GetCurrentStepNumber(),fGhostSafety,feLimited,
331  fEndTrack,track.GetVolume());
332  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
333  if(feLimited == kDoNot)
334  {
335  // Track is not on the boundary
336  fOnBoundary = false;
337  fGhostSafety = fGhostNavigator->ComputeSafety(fEndTrack.GetPosition());
338  }
339  else
340  {
341  // Track is on the boundary
342  fOnBoundary = true;
343  proposedSafety = fGhostSafety;
344  }
345  //xbug? proposedSafety = fGhostSafety;
346  if(feLimited == kUnique || feLimited == kSharedOther) {
347  *selection = CandidateForSelection;
348  }else if (feLimited == kSharedTransport) {
349  returnedStep *= (1.0 + 1.0e-9);
350  // Expand to disable its selection in Step Manager comparison
351  }
352 
353  }
354 
355  // ----------------------------------------------
356  // Returns the fGhostSafety as the proposedSafety
357  // The SteppingManager will take care of keeping
358  // the smallest one.
359  // ----------------------------------------------
360  return returnedStep;
361 
362  } else {
363  return DBL_MAX;
364  //not sensible! return -1.0;
365  }
366 
367 }
368 
369 
373 {
374  return -1.0;
375 }
376 
379 {
380  return 0;
381 }
382 
385 {
386  // Dummy ParticleChange ie: does nothing
387  // Expecting G4Transportation to move the track
388  pParticleChange->Initialize(track);
389  return pParticleChange;
390 
391  // return 0;
392 }
393 
394 void G4WeightCutOffProcess::CopyStep(const G4Step & step)
395 {
396  fGhostStep->SetTrack(step.GetTrack());
397  fGhostStep->SetStepLength(step.GetStepLength());
398  fGhostStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
399  fGhostStep->SetControlFlag(step.GetControlFlag());
400 
401  *fGhostPreStepPoint = *(step.GetPreStepPoint());
402  *fGhostPostStepPoint = *(step.GetPostStepPoint());
403 
404 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
405 // Set StepStatus for ghost world
406 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
407  if(fOnBoundary)
408  { fGhostPostStepPoint->SetStepStatus(fGeomBoundary); }
409  else if(fGhostPostStepPoint->GetStepStatus()==fGeomBoundary)
410  { fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc); }
411 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
412 }
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 &)
void SetTrack(G4Track *value)
G4int verboseLevel
Definition: G4VProcess.hh:368
G4double GetStepLength() const
const char * p
Definition: xmltok.h:285
G4StepStatus GetStepStatus() const
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
void SetParallelWorld(const G4String &parallelWorldName)
const G4VTouchable * GetTouchable() const
void SetStepStatus(const G4StepStatus aValue)
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:97
G4GLOB_DLL std::ostream G4cout
G4int GetCurrentStepNumber() const
G4VPhysicalVolume * GetPhysicalVolume() const
const G4String & GetName() 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 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)
#define DBL_MAX
Definition: templates.hh:83
const G4TouchableHandle & GetTouchableHandle() const
static void Update(G4FieldTrack *, const G4Track *)
G4GPILSelection