Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/biasing/GB06/src/GB06BOptnSplitAndKillByImportance.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 /// \file GB06BOptnSplitAndKillByImportance.cc
 28 /// \brief Implementation of the GB06BOptnSplitAndKillByImportance class
 29 
 30 #include "GB06BOptnSplitAndKillByImportance.hh"
 31 
 32 #include "G4BiasingProcessInterface.hh"
 33 #include "G4BiasingProcessSharedData.hh"
 34 #include "G4ParallelGeometriesLimiterProcess.hh"
 35 #include "G4ProcessManager.hh"
 36 #include "Randomize.hh"
 37 
 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 39 
 40 GB06BOptnSplitAndKillByImportance::GB06BOptnSplitAndKillByImportance(G4String name)
 41   : G4VBiasingOperation(name),
 42     fParallelWorldIndex(-1),
 43     fBiasingSharedData(nullptr),
 44     fParticleChange(),
 45     fDummyParticleChange()
 46 {}
 47 
 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 49 
 50 GB06BOptnSplitAndKillByImportance::~GB06BOptnSplitAndKillByImportance() {}
 51 
 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 53 
 54 G4double GB06BOptnSplitAndKillByImportance::DistanceToApplyOperation(const G4Track*, G4double,
 55                                                                      G4ForceCondition* condition)
 56 {
 57   // -- Remember the touchable history (ie geometry state) at the beginning of the step:
 58   // -- Start by getting the process handling the step limitation in parallel
 59   // -- geometries for the generic biasing:
 60   auto biasingLimiterProcess = fBiasingSharedData->GetParallelGeometriesLimiterProcess();
 61   fPreStepTouchableHistory =
 62     biasingLimiterProcess
 63       ->GetNavigator(fParallelWorldIndex)  // -- get the navigator of the geometry...
 64       ->CreateTouchableHistoryHandle();  // -- ...and collect the geometry state.
 65 
 66   // -- We request to be "forced" : meaning GenerateBiasingFinalState(...) will be called
 67   // -- anyway at the PostStepDoIt(...) stage (ie, when the track will have been moved to
 68   // -- its end point position) and, there, we will have to handle the decision to
 69   // -- split/kill if we are on a volume boundary, or do nothing, if we're not:
 70   *condition = Forced;
 71 
 72   // -- As we're forced, we make this returned length void:
 73   return DBL_MAX;
 74 }
 75 
 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 77 
 78 G4VParticleChange*
 79 GB06BOptnSplitAndKillByImportance::GenerateBiasingFinalState(const G4Track* track, const G4Step*)
 80 {
 81   // -- Given we used the "Forced" condition, this method is always called.
 82   // -- We check the status of the step in the parallel geometry, and apply
 83   // -- splitting/killing if the step has been limited by the geometry.
 84 
 85   // -- We first check if we limit the step in the parallel geometry:
 86   G4bool isLimiting =
 87     fBiasingSharedData->GetParallelGeometriesLimiterProcess()->GetIsLimiting(fParallelWorldIndex);
 88 
 89   // -- if not limiting, we leave the track unchanged:
 90   if (!isLimiting) {
 91     fDummyParticleChange.Initialize(*track);
 92     return &fDummyParticleChange;
 93   }
 94 
 95   // -- We collect the geometry state at the end point step:
 96   // -- Note this is the same call than in the DistanceToApplyOperation(...) for the
 97   // -- fPreStepTouchableHistory, but the navigator has pushed the track in the next
 98   // -- volume since then (even if the track is still on the boundary), and hence the
 99   // -- geometry state has changed.
100   auto biasingLimiterProcess = fBiasingSharedData->GetParallelGeometriesLimiterProcess();
101   fPostStepTouchableHistory =
102     biasingLimiterProcess->GetNavigator(fParallelWorldIndex)->CreateTouchableHistoryHandle();
103 
104   // -- We verify we are still in the "same" physical volume:
105   // -- remember : using a replica, we have "one" physical volume
106   // -- but representing many different placements, distinguished
107   // -- by replica number. By checking we are in the same physical
108   // -- volume, we verify the end step point has not reached the
109   // -- world volume of the parallel geometry:
110   if (fPreStepTouchableHistory->GetVolume() != fPostStepTouchableHistory->GetVolume()) {
111     // -- the track is leaving the volumes in which biasing is applied,
112     // -- we leave this track unchanged:
113     fDummyParticleChange.Initialize(*track);
114     return &fDummyParticleChange;
115   }
116 
117   // -------------------------------------------------------------------------------------
118   // -- At this stage, we know we have a particle crossing a boundary between two slices,
119   // -- each of this slice has a well defined importance : we apply the biasing.
120   // -- We will split if the track is entering a volume with higher importance, and
121   // -- kill (applying Russian roulette) in the other case.
122   // -------------------------------------------------------------------------------------
123 
124   // -- We start by getting the replica numbers:
125   G4int preReplicaNumber = fPreStepTouchableHistory->GetReplicaNumber();
126   G4int postReplicaNumber = fPostStepTouchableHistory->GetReplicaNumber();
127 
128   // -- and get the related importance we defined in the importance map:
129   G4int preImportance = (*fImportanceMap)[preReplicaNumber];
130   G4int postImportance = (*fImportanceMap)[postReplicaNumber];
131 
132   // -- Get track weight:
133   G4double initialWeight = track->GetWeight();
134 
135   // -- Initialize the "particle change" (it will communicate the new track state to
136   // -- the tracking):
137   fParticleChange.Initialize(*track);
138 
139   if (postImportance > preImportance) {
140     // -- We split :
141     G4int splittingFactor = postImportance / preImportance;
142 
143     // Define the tracks weight:
144     G4double weightOfTrack = initialWeight / splittingFactor;
145 
146     // Ask currect track weight to be changed to the new value:
147     fParticleChange.ProposeParentWeight(weightOfTrack);
148     // Now we clone this track (this is the actual splitting):
149     // we will then have the primary and clone of it, hence the
150     // splitting by a factor 2:
151     G4Track* clone = new G4Track(*track);
152     clone->SetWeight(weightOfTrack);
153     fParticleChange.AddSecondary(clone);
154     // -- Below's call added for safety & illustration : inform particle change to not
155     // -- modify the clone (ie : daughter) weight to make it that of the primary.
156     // -- Here this call is not mandatory as both tracks have same weights.
157     fParticleChange.SetSecondaryWeightByProcess(true);
158   }
159   else {
160     // -- We apply Russian roulette:
161     G4double survivingProbability = G4double(postImportance) / G4double(preImportance);
162 
163     // Shoot a random number (in ]0,1[ segment):
164     G4double random = G4UniformRand();
165     if (random > survivingProbability) {
166       // We ask for the the track to be killed:
167       fParticleChange.ProposeTrackStatus(fStopAndKill);
168     }
169     else {
170       // In this case, the track survives. We change its weight
171       // to conserve weight among killed and survival tracks:
172       fParticleChange.ProposeParentWeight(initialWeight / survivingProbability);
173     }
174   }
175 
176   return &fParticleChange;
177 }
178 
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180