Geant4_10
G4WeightWindowProcess.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: G4WeightWindowProcess.cc 77999 2013-12-02 08:18:01Z gcosmo $
28 //
29 // ----------------------------------------------------------------------
30 // GEANT 4 class source file
31 //
32 // G4WeightWindowProcess.cc
33 //
34 // ----------------------------------------------------------------------
35 
36 #include "G4WeightWindowProcess.hh"
38 #include "G4GeometryCell.hh"
40 #include "G4VTrackTerminator.hh"
41 #include "G4PlaceOfAction.hh"
42 #include "G4VWeightWindowStore.hh"
43 
44 #include "G4Step.hh"
45 #include "G4Navigator.hh"
46 #include "G4VTouchable.hh"
47 #include "G4VPhysicalVolume.hh"
48 #include "G4ParticleChange.hh"
49 #include "G4PathFinder.hh"
51 #include "G4StepPoint.hh"
52 #include "G4FieldTrackUpdator.hh"
53 
54 
56  const G4VWeightWindowAlgorithm &aWeightWindowAlgorithm,
57  const G4VWeightWindowStore &aWWStore,
58  const G4VTrackTerminator *TrackTerminator,
59  G4PlaceOfAction placeOfAction,
60  const G4String &aName, G4bool para)
61  : G4VProcess(aName),
62  fParticleChange(new G4ParticleChange),
63  fWeightWindowAlgorithm(aWeightWindowAlgorithm),
64  fWeightWindowStore(aWWStore),
65  fPostStepAction(0),
66  fPlaceOfAction(placeOfAction),
67  fGhostWorldName("NoParallelWorld"),fGhostWorld(0),
68  fGhostNavigator(0), fNavigatorID(-1), fFieldTrack('0')
69 {
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("G4WeightWindowProcess::G4WeightWindowProcess()",
81  "FatalError", FatalException,
82  "Failed allocation of G4ParticleChange !");
83  }
84  G4VProcess::pParticleChange = fParticleChange;
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  paraflag = para;
99 
100 }
101 
103 {
104 
105  delete fPostStepAction;
106  delete fParticleChange;
107  // delete fGhostStep;
108 
109 }
110 
111 
112 //------------------------------------------------------
113 //
114 // SetParallelWorld
115 //
116 //------------------------------------------------------
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 
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 //------------------------------------------------------
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("G4WeightWindowProcess::StartTracking",
159  "ProcParaWorld000",FatalException,
160  "G4WeightWindowProcess 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 
188  G4double ,
190 {
191 // *condition = Forced;
192 // return kInfinity;
193 
194 // *condition = StronglyForced;
195  *condition = Forced;
196  return DBL_MAX;
197 }
198 
201  const G4Step &aStep)
202 {
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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
220  }
221  else
222  {
223  // Do I need this ??????????????????????????????????????????????????????????
224  // fGhostNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
225  // ?????????????????????????????????????????????????????????????????????????
226 
227  // fPathFinder->ReLocate(track.GetPosition());
228 
229  // reuse the touchable
230  fNewGhostTouchable = fOldGhostTouchable;
231  }
232 
233  fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
234  fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
235 
236  }
237 
238  if (aStep.GetStepLength() > kCarTolerance)
239  {
240 // if ( ( fPlaceOfAction == onBoundaryAndCollision)
241 // || ( (fPlaceOfAction == onBoundary) &&
242 // (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary) )
243 // || ( (fPlaceOfAction == onCollision) &&
244 // (aStep.GetPostStepPoint()->GetStepStatus() != fGeomBoundary) ) )
245  if(paraflag) {
246  if ( ( fPlaceOfAction == onBoundaryAndCollision)
247  || ( (fPlaceOfAction == onBoundary) &&
248  (fGhostPostStepPoint->GetStepStatus() == fGeomBoundary) )
249  || ( (fPlaceOfAction == onCollision) &&
250  (fGhostPostStepPoint->GetStepStatus() != fGeomBoundary) ) )
251  {
252 
253 // G4StepPoint *postpoint = aStep.GetPostStepPoint();
254 
255 // G4GeometryCell postCell(*(postpoint->GetPhysicalVolume()),
256 // postpoint->GetTouchable()->GetReplicaNumber());
257 
258  G4GeometryCell postCell(*(fGhostPostStepPoint->GetPhysicalVolume()),
259  fGhostPostStepPoint->GetTouchable()->GetReplicaNumber());
260  G4Nsplit_Weight nw =
261  fWeightWindowAlgorithm.Calculate(aTrack.GetWeight(),
262  fWeightWindowStore.GetLowerWeight(postCell,
263  aTrack.GetKineticEnergy()));
264  fPostStepAction->DoIt(aTrack, fParticleChange, nw);
265  }
266  } else {
267  if ( ( fPlaceOfAction == onBoundaryAndCollision)
268  || ( (fPlaceOfAction == onBoundary) &&
270  || ( (fPlaceOfAction == onCollision) &&
271  (aStep.GetPostStepPoint()->GetStepStatus() != fGeomBoundary) ) )
272  {
273 
274  G4StepPoint *postpoint = aStep.GetPostStepPoint();
275 
276  G4GeometryCell postCell(*(postpoint->GetPhysicalVolume()),
277  postpoint->GetTouchable()->GetReplicaNumber());
278 
279  G4Nsplit_Weight nw =
280  fWeightWindowAlgorithm.Calculate(aTrack.GetWeight(),
281  fWeightWindowStore.GetLowerWeight(postCell,
282  aTrack.GetKineticEnergy()));
283  fPostStepAction->DoIt(aTrack, fParticleChange, nw);
284  }
285  }
286  }
287  return fParticleChange;
288 }
289 
291 {
292  fParticleChange->ProposeTrackStatus(fStopAndKill);
293 }
294 
296 {
298 }
299 
302  const G4Track& track, G4double previousStepSize, G4double currentMinimumStep,
303  G4double& proposedSafety, G4GPILSelection* selection)
304 {
305  if(paraflag) {
306  static G4ThreadLocal G4FieldTrack *endTrack_G4MT_TLS_ = 0 ; if (!endTrack_G4MT_TLS_) endTrack_G4MT_TLS_ = new G4FieldTrack ('0') ; G4FieldTrack &endTrack = *endTrack_G4MT_TLS_;
307  static G4ThreadLocal ELimited *eLimited_G4MT_TLS_ = 0 ; if (!eLimited_G4MT_TLS_) eLimited_G4MT_TLS_ = new ELimited ; ELimited &eLimited = *eLimited_G4MT_TLS_;
308 
309  *selection = NotCandidateForSelection;
310  G4double returnedStep = DBL_MAX;
311 
312  if (previousStepSize > 0.)
313  { fGhostSafety -= previousStepSize; }
314  // else
315  // { fGhostSafety = -1.; }
316  if (fGhostSafety < 0.) fGhostSafety = 0.0;
317 
318  // ------------------------------------------
319  // Determination of the proposed STEP LENGTH:
320  // ------------------------------------------
321  if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.)
322  {
323  // I have no chance to limit
324  returnedStep = currentMinimumStep;
325  fOnBoundary = false;
326  proposedSafety = fGhostSafety - currentMinimumStep;
327  }
328  else // (currentMinimumStep > fGhostSafety: I may limit the Step)
329  {
330  G4FieldTrackUpdator::Update(&fFieldTrack,&track);
331  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
332  // ComputeStep
333  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
334  returnedStep
335  = fPathFinder->ComputeStep(fFieldTrack,currentMinimumStep,fNavigatorID,
336  track.GetCurrentStepNumber(),fGhostSafety,eLimited,
337  endTrack,track.GetVolume());
338  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
339  if(eLimited == kDoNot)
340  {
341  // Track is not on the boundary
342  fOnBoundary = false;
343  fGhostSafety = fGhostNavigator->ComputeSafety(endTrack.GetPosition());
344  }
345  else
346  {
347  // Track is on the boundary
348  fOnBoundary = true;
349  proposedSafety = fGhostSafety;
350  }
351  //xbug? proposedSafety = fGhostSafety;
352  if(eLimited == kUnique || eLimited == kSharedOther) {
353  *selection = CandidateForSelection;
354  }else if (eLimited == kSharedTransport) {
355  returnedStep *= (1.0 + 1.0e-9);
356  // Expand to disable its selection in Step Manager comparison
357  }
358 
359  }
360 
361  // ----------------------------------------------
362  // Returns the fGhostSafety as the proposedSafety
363  // The SteppingManager will take care of keeping
364  // the smallest one.
365  // ----------------------------------------------
366  return returnedStep;
367 
368  } else {
369  return DBL_MAX;
370  // not sensible - time goes backwards! return -1.0;
371  }
372 
373 }
374 
378 {
379  return -1.0;
380 }
381 
383 AtRestDoIt(const G4Track&, const G4Step&)
384 {
385  return 0;
386 }
387 
389 AlongStepDoIt(const G4Track& track, const G4Step&)
390 {
391  // Dummy ParticleChange ie: does nothing
392  // Expecting G4Transportation to move the track
393  pParticleChange->Initialize(track);
394  return pParticleChange;
395 
396  // return 0;
397 }
398 
399 void G4WeightWindowProcess::CopyStep(const G4Step & step)
400 {
401  fGhostStep->SetTrack(step.GetTrack());
402  fGhostStep->SetStepLength(step.GetStepLength());
403  fGhostStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
404  fGhostStep->SetControlFlag(step.GetControlFlag());
405 
406  *fGhostPreStepPoint = *(step.GetPreStepPoint());
407  *fGhostPostStepPoint = *(step.GetPostStepPoint());
408 
409 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
410 // Set StepStatus for ghost world
411 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
412  if(fOnBoundary)
413  { fGhostPostStepPoint->SetStepStatus(fGeomBoundary); }
414  else if(fGhostPostStepPoint->GetStepStatus()==fGeomBoundary)
415  { fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc); }
416 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
417 }
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 G4VParticleChange * AtRestDoIt(const G4Track &, const G4Step &)
virtual void Initialize(const G4Track &)
void SetTrack(G4Track *value)
G4int verboseLevel
Definition: G4VProcess.hh:368
G4double GetStepLength() const
G4StepStatus GetStepStatus() const
ELimited
const G4ThreeVector & GetPosition() const
G4TouchableHandle CreateTouchableHandle(G4int navId) const
G4SteppingControl GetControlFlag() const
const G4VTouchable * GetTouchable() const
#define G4ThreadLocal
Definition: tls.hh:52
void SetStepStatus(const G4StepStatus aValue)
G4StepPoint * GetPreStepPoint() const
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
G4int GetCurrentStepNumber() const
G4VPhysicalVolume * GetPhysicalVolume() const
G4PlaceOfAction
const G4String & GetName() const
virtual G4double AtRestGetPhysicalInteractionLength(const G4Track &, G4ForceCondition *)
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)
virtual const G4String & GetName() const
Definition: G4Step.hh:76
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
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)
void DoIt(const G4Track &aTrack, G4ParticleChange *aParticleChange, const G4Nsplit_Weight &nw)
const G4ThreeVector & GetMomentumDirection() const
virtual G4VParticleChange * AlongStepDoIt(const G4Track &, const G4Step &)
G4Navigator * GetNavigator(const G4String &worldName)
virtual void KillTrack() const
G4StepPoint * GetPostStepPoint() const
virtual G4int GetReplicaNumber(G4int depth=0) const
Definition: G4VTouchable.cc:58
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
G4VPhysicalVolume * GetVolume() const
virtual G4Nsplit_Weight Calculate(G4double init_w, G4double lowerWeightBound) const =0
G4double GetWeight() const
#define G4endl
Definition: G4ios.hh:61
G4WeightWindowProcess(const G4VWeightWindowAlgorithm &aWeightWindowAlgorithm, const G4VWeightWindowStore &aWWStore, const G4VTrackTerminator *TrackTerminator, G4PlaceOfAction placeOfAction, const G4String &aName="WeightWindowProcess", G4bool para=false)
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
void SetTotalEnergyDeposit(G4double value)
virtual G4double GetLowerWeight(const G4GeometryCell &gCell, G4double partEnergy) const =0
double G4double
Definition: G4Types.hh:76
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
void ProposeTrackStatus(G4TrackStatus status)
virtual G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4double pProposedMaxLength=DBL_MAX, const G4bool keepState=true)
G4ForceCondition
G4Track * GetTrack() const
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
void SetTouchableHandle(const G4TouchableHandle &apValue)
#define DBL_MAX
Definition: templates.hh:83
const G4TouchableHandle & GetTouchableHandle() const
static void Update(G4FieldTrack *, const G4Track *)
G4GPILSelection
void SetParallelWorld(G4String parallelWorldName)