source: trunk/source/processes/biasing/src/G4ImportanceProcess.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: 16.4 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: G4ImportanceProcess.cc,v 1.4 2008/04/21 09:10:28 ahoward Exp $
28// GEANT4 tag $Name: geant4-09-03 $
29//
30// ----------------------------------------------------------------------
31// GEANT 4 class source file
32//
33// G4ImportanceProcess.cc
34//
35// ----------------------------------------------------------------------
36
37#include "G4ImportanceProcess.hh"
38#include "G4VImportanceAlgorithm.hh"
39#include "G4GeometryCell.hh"
40#include "G4SamplingPostStepAction.hh"
41#include "G4VTrackTerminator.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
55G4ImportanceProcess::
56G4ImportanceProcess(const G4VImportanceAlgorithm &aImportanceAlgorithm,
57                        const G4VIStore &aIstore,
58                        const G4VTrackTerminator *TrackTerminator,
59                        const G4String &aName, G4bool para)
60 : G4VProcess(aName),
61   fParticleChange(new G4ParticleChange),
62   fImportanceAlgorithm(aImportanceAlgorithm),
63   fIStore(aIstore),
64   fPostStepAction(0),
65   fGhostNavigator(0), fNavigatorID(-1), fFieldTrack('0'),
66   paraflag(para)
67 
68{
69  if (TrackTerminator)
70  {
71    fPostStepAction = new G4SamplingPostStepAction(*TrackTerminator);
72  }
73  else
74  {
75    fPostStepAction = new G4SamplingPostStepAction(*this);
76  }
77  if (!fParticleChange)
78  {
79    G4Exception("G4ImportanceProcess::G4ImportanceProcess()",
80                "FatalError", FatalException,
81                "Failed allocation of G4ParticleChange !");
82  }
83  G4VProcess::pParticleChange = fParticleChange;
84
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  G4cout << " importance process paraflag is: " << paraflag << G4endl;
99
100}
101
102G4ImportanceProcess::~G4ImportanceProcess()
103{
104
105  delete fPostStepAction;
106  delete fParticleChange;
107  delete fGhostStep;
108
109}
110
111
112
113//------------------------------------------------------
114//
115// SetParallelWorld
116//
117//------------------------------------------------------
118void G4ImportanceProcess::
119SetParallelWorld(G4String parallelWorldName)
120{
121//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
122// Get pointers of the parallel world and its navigator
123//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
124  fGhostWorldName = parallelWorldName;
125  fGhostWorld = fTransportationManager->GetParallelWorld(fGhostWorldName);
126  fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
127//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
128}
129
130void G4ImportanceProcess::
131SetParallelWorld(G4VPhysicalVolume* parallelWorld)
132{
133//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
134// Get pointer of navigator
135//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
136  fGhostWorldName = parallelWorld->GetName();
137  fGhostWorld = parallelWorld;
138  fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
139//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
140}
141
142//------------------------------------------------------
143//
144// StartTracking
145//
146//------------------------------------------------------
147void G4ImportanceProcess::StartTracking(G4Track* trk)
148{
149//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
150// Activate navigator and get the navigator ID
151//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
152// G4cout << " G4ParallelWorldScoringProcess::StartTracking" << G4endl;
153
154  if(paraflag) {
155    if(fGhostNavigator)
156      { fNavigatorID = fTransportationManager->ActivateNavigator(fGhostNavigator); }
157    else
158      {
159        G4Exception("G4ImportanceProcess::StartTracking",
160                    "ProcParaWorld000",FatalException,
161                    "G4ImportanceProcess is used for tracking without having a parallel world assigned");
162      }
163//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
164
165// G4cout << "G4ParallelWorldScoringProcess::StartTracking <<<<<<<<<<<<<<<<<< " << G4endl;
166//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
167// Let PathFinder initialize
168//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
169    fPathFinder->PrepareNewTrack(trk->GetPosition(),trk->GetMomentumDirection());
170//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
171
172//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
173// Setup initial touchables for the first step
174//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
175    fOldGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
176    fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
177    fNewGhostTouchable = fOldGhostTouchable;
178    fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
179
180  // Initialize
181    fGhostSafety = -1.;
182    fOnBoundary = false;
183  }
184
185}
186
187
188G4double G4ImportanceProcess::
189PostStepGetPhysicalInteractionLength(const G4Track& ,
190                                     G4double   ,
191                                     G4ForceCondition* condition)
192{
193//   *condition = Forced;
194//   return kInfinity;
195
196//  *condition = StronglyForced;
197  *condition = Forced;
198  return DBL_MAX;
199}
200 
201G4VParticleChange *
202G4ImportanceProcess::PostStepDoIt(const G4Track &aTrack,
203                                      const G4Step &aStep)
204{
205  fParticleChange->Initialize(aTrack);
206
207  if(paraflag) {
208    fOldGhostTouchable = fGhostPostStepPoint->GetTouchableHandle();
209    //xbug?    fOnBoundary = false;
210    CopyStep(aStep);
211   
212    if(fOnBoundary)
213      {
214//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
215// Locate the point and get new touchable
216//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
217  //??  fPathFinder->Locate(step.GetPostStepPoint()->GetPosition(),
218  //??                      step.GetPostStepPoint()->GetMomentumDirection());
219        fNewGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
220        //AH    G4cout << " on boundary " << G4endl;
221//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
222      }
223    else
224      {
225        // Do I need this ??????????????????????????????????????????????????????????
226        // fGhostNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
227        // ?????????????????????????????????????????????????????????????????????????
228       
229        // fPathFinder->ReLocate(track.GetPosition());
230       
231        // reuse the touchable
232        fNewGhostTouchable = fOldGhostTouchable;
233        //AH    G4cout << " NOT on boundary " << G4endl;
234      }
235   
236    fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
237    fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
238   
239//   if ( (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
240//     && (aStep.GetStepLength() > kCarTolerance) )
241//   {
242//     if (aTrack.GetTrackStatus()==fStopAndKill)
243//     {
244//       G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
245//              << "          StopAndKill track." << G4endl;
246//     }
247
248//     G4StepPoint *prepoint  = aStep.GetPreStepPoint();
249//     G4StepPoint *postpoint = aStep.GetPostStepPoint();
250//     G4GeometryCell prekey(*(prepoint->GetPhysicalVolume()),
251//                          prepoint->GetTouchable()->GetReplicaNumber());
252//     G4GeometryCell postkey(*(postpoint->GetPhysicalVolume()),
253//                           postpoint->GetTouchable()->GetReplicaNumber());
254
255
256    if ( (fGhostPostStepPoint->GetStepStatus() == fGeomBoundary)
257         && (aStep.GetStepLength() > kCarTolerance) )
258      {
259        if (aTrack.GetTrackStatus()==fStopAndKill)
260          {
261            G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
262                   << "          StopAndKill track. on boundary" << G4endl;
263          }
264       
265        G4GeometryCell prekey(*(fGhostPreStepPoint->GetPhysicalVolume()), 
266                              fGhostPreStepPoint->GetTouchable()->GetReplicaNumber());
267        G4GeometryCell postkey(*(fGhostPostStepPoint->GetPhysicalVolume()), 
268                               fGhostPostStepPoint->GetTouchable()->GetReplicaNumber());
269       
270        //AH
271        /*
272        G4cout << G4endl;
273        G4cout << G4endl;
274        G4cout << " inside parallel importance process " << aTrack.GetCurrentStepNumber() << G4endl;
275        G4cout << G4endl;
276        G4cout << G4endl;
277        G4cout << " prekey: " << fGhostPreStepPoint->GetPhysicalVolume()->GetName() << " replica:"
278               << fGhostPreStepPoint->GetTouchable()->GetReplicaNumber() << G4endl;
279        G4cout << " prekey ISTORE: " << fIStore.GetImportance(prekey) << G4endl;
280        G4cout << " postkey: " << G4endl;
281        G4cout << " postkey ISTORE: " << fIStore.GetImportance(postkey) << G4endl;
282        */
283        //AH
284        G4Nsplit_Weight nw = fImportanceAlgorithm.
285          Calculate(fIStore.GetImportance(prekey),
286                    fIStore.GetImportance(postkey), 
287                    aTrack.GetWeight());
288        //AH
289        /*
290        G4cout << " prekey weight: " << fIStore.GetImportance(prekey)
291               << " postkey weight: " << fIStore.GetImportance(postkey)
292               << " split weight: " << nw << G4endl;
293        G4cout << " before poststepaction " << G4endl;
294        */
295        //AH
296        fPostStepAction->DoIt(aTrack, fParticleChange, nw);
297        //AH    G4cout << " after post step do it " << G4endl;
298      }
299  } else {
300    if ( (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
301         && (aStep.GetStepLength() > kCarTolerance) )
302      {
303        //AH    G4cout << " inside non-parallel importance process " << G4endl;
304        if (aTrack.GetTrackStatus()==fStopAndKill)
305          {
306            G4cout << "WARNING - G4ImportanceProcess::PostStepDoIt()"
307                   << "          StopAndKill track. on boundary non-parallel" << G4endl;
308          }
309       
310        G4StepPoint *prepoint  = aStep.GetPreStepPoint();
311        G4StepPoint *postpoint = aStep.GetPostStepPoint();
312       
313        G4GeometryCell prekey(*(prepoint->GetPhysicalVolume()), 
314                              prepoint->GetTouchable()->GetReplicaNumber());
315        G4GeometryCell postkey(*(postpoint->GetPhysicalVolume()), 
316                               postpoint->GetTouchable()->GetReplicaNumber());
317       
318        G4Nsplit_Weight nw = fImportanceAlgorithm.
319          Calculate(fIStore.GetImportance(prekey),
320                    fIStore.GetImportance(postkey), 
321                    aTrack.GetWeight());
322        //AH
323        /*
324        G4cout << " prekey weight: " << fIStore.GetImportance(prekey)
325               << " postkey weight: " << fIStore.GetImportance(postkey)
326               << " split weight: " << nw << G4endl;
327        G4cout << " before poststepaction 2 " << G4endl;
328        */
329        //AH
330        fPostStepAction->DoIt(aTrack, fParticleChange, nw);
331        //AH    G4cout << " after poststepaction 2 " << G4endl;
332      }
333  }
334  return fParticleChange;
335}
336
337void G4ImportanceProcess::KillTrack() const
338{
339  fParticleChange->ProposeTrackStatus(fStopAndKill);
340}
341
342const G4String &G4ImportanceProcess::GetName() const
343{
344  return theProcessName;
345}
346
347G4double G4ImportanceProcess::
348AlongStepGetPhysicalInteractionLength(
349            const G4Track& track, G4double  previousStepSize, G4double  currentMinimumStep,
350            G4double& proposedSafety, G4GPILSelection* selection)
351{
352
353  //AH  G4cout << " along step physical interaction length " << G4endl;
354
355  if(paraflag) {
356    static G4FieldTrack endTrack('0');
357    static ELimited eLimited;
358   
359    *selection = NotCandidateForSelection;
360    G4double returnedStep = DBL_MAX;
361   
362    if (previousStepSize > 0.)
363      { fGhostSafety -= previousStepSize; }
364    //  else
365    //  { fGhostSafety = -1.; }
366    if (fGhostSafety < 0.) fGhostSafety = 0.0;
367   
368    // ------------------------------------------
369    // Determination of the proposed STEP LENGTH:
370    // ------------------------------------------
371    if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.)
372      {
373        // I have no chance to limit
374        returnedStep = currentMinimumStep;
375        fOnBoundary = false;
376        proposedSafety = fGhostSafety - currentMinimumStep;
377        //AH    G4cout << " step not limited, why? " << G4endl;
378      }
379    else // (currentMinimumStep > fGhostSafety: I may limit the Step)
380      {
381        G4FieldTrackUpdator::Update(&fFieldTrack,&track);
382        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
383        // ComputeStep
384        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
385        returnedStep
386          = fPathFinder->ComputeStep(fFieldTrack,currentMinimumStep,fNavigatorID,
387                                     track.GetCurrentStepNumber(),fGhostSafety,eLimited,
388                                     endTrack,track.GetVolume());
389        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
390        if(eLimited == kDoNot)
391          {
392            //AH            G4cout << " computestep came back with not-boundary " << G4endl;
393            // Track is not on the boundary
394            fOnBoundary = false;
395            fGhostSafety = fGhostNavigator->ComputeSafety(endTrack.GetPosition());
396          }
397        else
398          {
399            // Track is on the boundary
400            //AH            G4cout << " FOUND A BOUNDARY ! " << track.GetCurrentStepNumber() << G4endl;
401            fOnBoundary = true;
402            // proposedSafety = fGhostSafety;
403          }
404        proposedSafety = fGhostSafety;
405        if(eLimited == kUnique || eLimited == kSharedOther) {
406          *selection = CandidateForSelection;
407        }else if (eLimited == kSharedTransport) { 
408          returnedStep *= (1.0 + 1.0e-9); 
409          // Expand to disable its selection in Step Manager comparison
410        }
411       
412      }
413
414  // ----------------------------------------------
415  // Returns the fGhostSafety as the proposedSafety
416  // The SteppingManager will take care of keeping
417  // the smallest one.
418  // ----------------------------------------------
419    return returnedStep;
420
421  } else {
422
423    return DBL_MAX;
424
425  }
426
427}
428
429G4double G4ImportanceProcess::
430AtRestGetPhysicalInteractionLength(const G4Track& ,
431                                   G4ForceCondition*) 
432{
433  return -1.0;
434}
435 
436G4VParticleChange* G4ImportanceProcess::
437AtRestDoIt(const G4Track&, const G4Step&) 
438{
439  return 0;
440}
441
442G4VParticleChange* G4ImportanceProcess::
443AlongStepDoIt(const G4Track& aTrack, const G4Step& )
444{
445  // Dummy ParticleChange ie: does nothing
446  // Expecting G4Transportation to move the track
447  //AH  G4cout << " along step do it " << G4endl;
448  pParticleChange->Initialize(aTrack);
449  return pParticleChange; 
450  //  return 0;
451}
452
453void G4ImportanceProcess::CopyStep(const G4Step & step)
454{
455  fGhostStep->SetTrack(step.GetTrack());
456  fGhostStep->SetStepLength(step.GetStepLength());
457  fGhostStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
458  fGhostStep->SetControlFlag(step.GetControlFlag());
459
460  *fGhostPreStepPoint = *(step.GetPreStepPoint());
461  *fGhostPostStepPoint = *(step.GetPostStepPoint());
462
463//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
464// Set StepStatus for ghost world
465//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
466  if(fOnBoundary)
467  { fGhostPostStepPoint->SetStepStatus(fGeomBoundary); }
468  else if(fGhostPostStepPoint->GetStepStatus()==fGeomBoundary)
469  { fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc); }
470//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
471}
Note: See TracBrowser for help on using the repository browser.