source: trunk/source/processes/biasing/src/G4WeightWindowProcess.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: 14.6 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: G4WeightWindowProcess.cc,v 1.4 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// G4WeightWindowProcess.cc
34//
35// ----------------------------------------------------------------------
36
37#include "G4WeightWindowProcess.hh"
38#include "G4VWeightWindowAlgorithm.hh"
39#include "G4GeometryCell.hh"
40#include "G4SamplingPostStepAction.hh"
41#include "G4VTrackTerminator.hh"
42#include "G4PlaceOfAction.hh"
43#include "G4VWeightWindowStore.hh"
44
45#include "G4Step.hh"
46#include "G4Navigator.hh"
47#include "G4VTouchable.hh"
48#include "G4VPhysicalVolume.hh"
49#include "G4ParticleChange.hh"
50#include "G4PathFinder.hh"
51#include "G4TransportationManager.hh"
52#include "G4StepPoint.hh"
53#include "G4FieldTrackUpdator.hh"
54
55
56G4WeightWindowProcess::G4WeightWindowProcess(
57                    const G4VWeightWindowAlgorithm &aWeightWindowAlgorithm,
58                    const G4VWeightWindowStore &aWWStore,
59                    const G4VTrackTerminator *TrackTerminator,
60                          G4PlaceOfAction placeOfAction,
61                    const G4String &aName, G4bool para)
62 : G4VProcess(aName),
63   fParticleChange(new G4ParticleChange),
64   fWeightWindowAlgorithm(aWeightWindowAlgorithm),
65   fWeightWindowStore(aWWStore),
66   fPostStepAction(0),
67   fPlaceOfAction(placeOfAction),
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
102G4WeightWindowProcess::~G4WeightWindowProcess()
103{
104
105  delete fPostStepAction;
106  delete fParticleChange;
107  delete fGhostStep;
108
109}
110
111
112//------------------------------------------------------
113//
114// SetParallelWorld
115//
116//------------------------------------------------------
117void G4WeightWindowProcess::
118SetParallelWorld(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
129void G4WeightWindowProcess::
130SetParallelWorld(G4VPhysicalVolume* parallelWorld)
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//------------------------------------------------------
146void G4WeightWindowProcess::StartTracking(G4Track* trk)
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
186G4double G4WeightWindowProcess::
187PostStepGetPhysicalInteractionLength(const G4Track& ,
188                                     G4double   ,
189                                     G4ForceCondition* condition)
190{
191//   *condition = Forced;
192//   return kInfinity;
193
194//  *condition = StronglyForced;
195  *condition = Forced;
196  return DBL_MAX;
197}
198 
199G4VParticleChange *
200G4WeightWindowProcess::PostStepDoIt(const G4Track &aTrack,
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.GetLowerWeitgh(postCell,
263                                                    aTrack.GetKineticEnergy()));
264          fPostStepAction->DoIt(aTrack, fParticleChange, nw);
265        }
266    } else {
267      if ( ( fPlaceOfAction == onBoundaryAndCollision)
268           || ( (fPlaceOfAction == onBoundary) && 
269                (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary) )
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.GetLowerWeitgh(postCell,
282                                                                               aTrack.GetKineticEnergy()));
283          fPostStepAction->DoIt(aTrack, fParticleChange, nw);
284        }
285    }
286  }
287  return fParticleChange;
288}
289
290void G4WeightWindowProcess::KillTrack() const
291{
292  fParticleChange->ProposeTrackStatus(fStopAndKill);
293}
294
295const G4String &G4WeightWindowProcess::GetName() const
296{
297  return G4VProcess::GetProcessName();
298}
299
300G4double G4WeightWindowProcess::
301AlongStepGetPhysicalInteractionLength(
302            const G4Track& track, G4double  previousStepSize, G4double  currentMinimumStep,
303            G4double& proposedSafety, G4GPILSelection* selection)
304{
305  if(paraflag) {
306    static G4FieldTrack endTrack('0');
307    static ELimited eLimited;
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
375G4double G4WeightWindowProcess::
376AtRestGetPhysicalInteractionLength(const G4Track& ,
377                                   G4ForceCondition*) 
378{
379  return -1.0;
380}
381 
382G4VParticleChange* G4WeightWindowProcess::
383AtRestDoIt(const G4Track&, const G4Step&) 
384{
385  return 0;
386}
387
388G4VParticleChange* G4WeightWindowProcess::
389AlongStepDoIt(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
399void 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}
Note: See TracBrowser for help on using the repository browser.