source: trunk/source/processes/biasing/src/G4WeightCutOffProcess.cc @ 1253

Last change on this file since 1253 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
File size: 13.7 KB
Line 
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,v 1.3 2008/04/21 09:10:29 ahoward Exp $
28// GEANT4 tag $Name: geant4-09-03 $
29//
30// ----------------------------------------------------------------------
31// GEANT 4 class source file
32//
33// G4WeightCutOffProcess.cc
34//
35// ----------------------------------------------------------------------
36
37#include "G4WeightCutOffProcess.hh"
38//#include "G4VScorer.hh"
39#include "G4GeometryCellStep.hh"
40//#include "G4GCellFinder.hh"
41#include "G4TouchableHandle.hh"
42#include "G4VIStore.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"
50#include "G4TransportationManager.hh"
51#include "G4StepPoint.hh"
52#include "G4FieldTrackUpdator.hh"
53
54
55G4WeightCutOffProcess::
56G4WeightCutOffProcess(G4double wsurvival,
57                      G4double wlimit,
58                      G4double isource,
59                      G4VIStore *istore,
60                      //                      const G4VGCellFinder &aGCellFinder,
61                      const G4String &aName, G4bool para)
62  : G4VProcess(aName), 
63    fParticleChange(new G4ParticleChange),
64    fWeightSurvival(wsurvival),
65    fWeightLimit(wlimit),
66    fSourceImportance(isource),
67    fIStore(istore),
68    //    fGCellFinder(aGCellFinder),
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
78  G4VProcess::pParticleChange = fParticleChange;
79
80  fGhostStep = new G4Step();
81  fGhostPreStepPoint = fGhostStep->GetPreStepPoint();
82  fGhostPostStepPoint = fGhostStep->GetPostStepPoint();
83
84  fTransportationManager = G4TransportationManager::GetTransportationManager();
85  fPathFinder = G4PathFinder::GetInstance();
86
87  if (verboseLevel>0)
88  {
89    G4cout << GetProcessName() << " is created " << G4endl;
90  }
91
92  paraflag = para;
93
94}
95
96G4WeightCutOffProcess::~G4WeightCutOffProcess()
97{
98  delete fParticleChange;
99  delete fGhostStep;
100}
101
102
103//------------------------------------------------------
104//
105// SetParallelWorld
106//
107//------------------------------------------------------
108void G4WeightCutOffProcess::
109SetParallelWorld(G4String parallelWorldName)
110{
111//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
112// Get pointers of the parallel world and its navigator
113//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
114  fGhostWorldName = parallelWorldName;
115  fGhostWorld = fTransportationManager->GetParallelWorld(fGhostWorldName);
116  fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
117//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
118}
119
120void G4WeightCutOffProcess::
121SetParallelWorld(G4VPhysicalVolume* parallelWorld)
122{
123//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
124// Get pointer of navigator
125//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
126  fGhostWorldName = parallelWorld->GetName();
127  fGhostWorld = parallelWorld;
128  fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
129//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
130}
131
132//------------------------------------------------------
133//
134// StartTracking
135//
136//------------------------------------------------------
137void G4WeightCutOffProcess::StartTracking(G4Track* trk)
138{
139//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
140// Activate navigator and get the navigator ID
141//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
142// G4cout << " G4ParallelWorldScoringProcess::StartTracking" << G4endl;
143
144  if(paraflag) {
145    if(fGhostNavigator)
146      { fNavigatorID = fTransportationManager->ActivateNavigator(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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
159    fPathFinder->PrepareNewTrack(trk->GetPosition(),trk->GetMomentumDirection());
160//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
161
162//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
163// Setup initial touchables for the first step
164//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
165    fOldGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
166    fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
167    fNewGhostTouchable = fOldGhostTouchable;
168    fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
169   
170    // Initialize
171    fGhostSafety = -1.;
172    fOnBoundary = false;
173  }
174}
175
176
177G4double G4WeightCutOffProcess::
178PostStepGetPhysicalInteractionLength(const G4Track &,
179                                     G4double, G4ForceCondition* condition)
180{
181//   *condition = Forced;
182//   return kInfinity;
183
184//  *condition = StronglyForced;
185  *condition = Forced;
186  return DBL_MAX;
187}
188 
189G4VParticleChange * 
190G4WeightCutOffProcess::PostStepDoIt(const G4Track& aTrack, 
191                                    const G4Step &aStep)
192{
193  fParticleChange->Initialize(aTrack);
194
195  if(paraflag) {
196    fOldGhostTouchable = fGhostPostStepPoint->GetTouchableHandle();
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());
207        fNewGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
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
219        fNewGhostTouchable = fOldGhostTouchable;
220      }
221   
222    fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
223    fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
224
225  }
226
227  if(paraflag) {
228    G4GeometryCell postCell(*(fGhostPostStepPoint->GetPhysicalVolume()), 
229                            fGhostPostStepPoint->GetTouchable()->GetReplicaNumber());
230
231
232    //  G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(aStep);
233    //  G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(fGhostStep);
234    G4double R = fSourceImportance;
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          {
250            fParticleChange->ProposeTrackStatus(fStopAndKill);
251          }
252        else
253          {
254            fParticleChange->ProposeWeight(ws);
255          }                 
256      }
257  } else {
258
259    G4GeometryCell postCell(*(aStep.GetPostStepPoint()->GetPhysicalVolume()), 
260                            aStep.GetPostStepPoint()->GetTouchable()->GetReplicaNumber());
261
262    //  G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(aStep);
263    //  G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(fGhostStep);
264    G4double R = fSourceImportance;
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          {
280            fParticleChange->ProposeTrackStatus(fStopAndKill);
281          }
282        else
283          {
284            fParticleChange->ProposeWeight(ws);
285          }                 
286      }
287  }
288
289  return fParticleChange;
290
291}
292
293const G4String &G4WeightCutOffProcess::GetName() const
294{
295  return theProcessName;
296}
297
298G4double G4WeightCutOffProcess::
299AlongStepGetPhysicalInteractionLength(
300            const G4Track& track, G4double  previousStepSize, G4double  currentMinimumStep,
301            G4double& proposedSafety, G4GPILSelection* selection)
302{
303  if(paraflag) {
304    static G4FieldTrack endTrack('0');
305    static ELimited eLimited;
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      {
328        G4FieldTrackUpdator::Update(&fFieldTrack,&track);
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 
374G4double G4WeightCutOffProcess::
375AtRestGetPhysicalInteractionLength(const G4Track& ,
376                                   G4ForceCondition*)
377{
378  return -1.0;
379}
380 
381G4VParticleChange*
382G4WeightCutOffProcess::AtRestDoIt(const G4Track&, const G4Step&)
383{
384  return 0;
385}
386
387G4VParticleChange*
388G4WeightCutOffProcess::AlongStepDoIt(const G4Track& track, const G4Step&)
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
398void G4WeightCutOffProcess::CopyStep(const G4Step & step)
399{
400  fGhostStep->SetTrack(step.GetTrack());
401  fGhostStep->SetStepLength(step.GetStepLength());
402  fGhostStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
403  fGhostStep->SetControlFlag(step.GetControlFlag());
404
405  *fGhostPreStepPoint = *(step.GetPreStepPoint());
406  *fGhostPostStepPoint = *(step.GetPostStepPoint());
407
408//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
409// Set StepStatus for ghost world
410//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
411  if(fOnBoundary)
412  { fGhostPostStepPoint->SetStepStatus(fGeomBoundary); }
413  else if(fGhostPostStepPoint->GetStepStatus()==fGeomBoundary)
414  { fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc); }
415//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
416}
Note: See TracBrowser for help on using the repository browser.