| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 /// \file GB06BOptnSplitAndKillByImportance.cc 27 /// \file GB06BOptnSplitAndKillByImportance.cc
28 /// \brief Implementation of the GB06BOptnSpli 28 /// \brief Implementation of the GB06BOptnSplitAndKillByImportance class
29 29
30 #include "GB06BOptnSplitAndKillByImportance.hh 30 #include "GB06BOptnSplitAndKillByImportance.hh"
>> 31 #include "Randomize.hh"
>> 32
31 33
32 #include "G4BiasingProcessInterface.hh" 34 #include "G4BiasingProcessInterface.hh"
33 #include "G4BiasingProcessSharedData.hh" <<
34 #include "G4ParallelGeometriesLimiterProcess.h 35 #include "G4ParallelGeometriesLimiterProcess.hh"
>> 36 #include "G4BiasingProcessSharedData.hh"
35 #include "G4ProcessManager.hh" 37 #include "G4ProcessManager.hh"
36 #include "Randomize.hh" << 38
37 39
38 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
39 41
40 GB06BOptnSplitAndKillByImportance::GB06BOptnSp 42 GB06BOptnSplitAndKillByImportance::GB06BOptnSplitAndKillByImportance(G4String name)
41 : G4VBiasingOperation(name), << 43 : G4VBiasingOperation (name),
42 fParallelWorldIndex(-1), << 44 fParallelWorldIndex ( -1 ),
43 fBiasingSharedData(nullptr), << 45 fBiasingSharedData ( nullptr ),
44 fParticleChange(), << 46 fParticleChange(),
45 fDummyParticleChange() << 47 fDummyParticleChange()
46 {} 48 {}
47 49
48 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 51
50 GB06BOptnSplitAndKillByImportance::~GB06BOptnS << 52 GB06BOptnSplitAndKillByImportance::~GB06BOptnSplitAndKillByImportance()
>> 53 {}
51 54
52 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 56
54 G4double GB06BOptnSplitAndKillByImportance::Di << 57 G4double GB06BOptnSplitAndKillByImportance::
55 << 58 DistanceToApplyOperation( const G4Track*,
>> 59 G4double,
>> 60 G4ForceCondition* condition)
56 { 61 {
>> 62
57 // -- Remember the touchable history (ie geo 63 // -- Remember the touchable history (ie geometry state) at the beginning of the step:
58 // -- Start by getting the process handling 64 // -- Start by getting the process handling the step limitation in parallel
59 // -- geometries for the generic biasing: 65 // -- geometries for the generic biasing:
60 auto biasingLimiterProcess = fBiasingSharedD << 66 auto biasingLimiterProcess = fBiasingSharedData->GetParallelGeometriesLimiterProcess();
61 fPreStepTouchableHistory = << 67 fPreStepTouchableHistory =
62 biasingLimiterProcess 68 biasingLimiterProcess
63 ->GetNavigator(fParallelWorldIndex) // << 69 ->GetNavigator( fParallelWorldIndex ) // -- get the navigator of the geometry...
64 ->CreateTouchableHistoryHandle(); // -- << 70 ->CreateTouchableHistoryHandle(); // -- ...and collect the geometry state.
65 << 71
66 // -- We request to be "forced" : meaning Ge 72 // -- We request to be "forced" : meaning GenerateBiasingFinalState(...) will be called
67 // -- anyway at the PostStepDoIt(...) stage 73 // -- anyway at the PostStepDoIt(...) stage (ie, when the track will have been moved to
68 // -- its end point position) and, there, we 74 // -- its end point position) and, there, we will have to handle the decision to
69 // -- split/kill if we are on a volume bound 75 // -- split/kill if we are on a volume boundary, or do nothing, if we're not:
70 *condition = Forced; 76 *condition = Forced;
71 77
72 // -- As we're forced, we make this returned 78 // -- As we're forced, we make this returned length void:
73 return DBL_MAX; 79 return DBL_MAX;
74 } 80 }
75 81
76 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
77 83
78 G4VParticleChange* << 84 G4VParticleChange*
79 GB06BOptnSplitAndKillByImportance::GenerateBia << 85 GB06BOptnSplitAndKillByImportance::GenerateBiasingFinalState( const G4Track* track,
>> 86 const G4Step* )
80 { 87 {
81 // -- Given we used the "Forced" condition, 88 // -- Given we used the "Forced" condition, this method is always called.
82 // -- We check the status of the step in the 89 // -- We check the status of the step in the parallel geometry, and apply
83 // -- splitting/killing if the step has been 90 // -- splitting/killing if the step has been limited by the geometry.
84 91
85 // -- We first check if we limit the step in 92 // -- We first check if we limit the step in the parallel geometry:
86 G4bool isLimiting = << 93 G4bool isLimiting = fBiasingSharedData
87 fBiasingSharedData->GetParallelGeometriesL << 94 ->GetParallelGeometriesLimiterProcess()
>> 95 ->GetIsLimiting( fParallelWorldIndex );
88 96
89 // -- if not limiting, we leave the track un 97 // -- if not limiting, we leave the track unchanged:
90 if (!isLimiting) { << 98 if ( !isLimiting )
91 fDummyParticleChange.Initialize(*track); << 99 {
92 return &fDummyParticleChange; << 100 fDummyParticleChange.Initialize( *track );
93 } << 101 return &fDummyParticleChange;
94 << 102 }
>> 103
95 // -- We collect the geometry state at the e 104 // -- We collect the geometry state at the end point step:
96 // -- Note this is the same call than in the 105 // -- Note this is the same call than in the DistanceToApplyOperation(...) for the
97 // -- fPreStepTouchableHistory, but the navi 106 // -- fPreStepTouchableHistory, but the navigator has pushed the track in the next
98 // -- volume since then (even if the track i 107 // -- volume since then (even if the track is still on the boundary), and hence the
99 // -- geometry state has changed. 108 // -- geometry state has changed.
100 auto biasingLimiterProcess = fBiasingSharedD 109 auto biasingLimiterProcess = fBiasingSharedData->GetParallelGeometriesLimiterProcess();
101 fPostStepTouchableHistory = << 110 fPostStepTouchableHistory =
102 biasingLimiterProcess->GetNavigator(fParal << 111 biasingLimiterProcess
>> 112 ->GetNavigator( fParallelWorldIndex )
>> 113 ->CreateTouchableHistoryHandle();
103 114
>> 115
104 // -- We verify we are still in the "same" p 116 // -- We verify we are still in the "same" physical volume:
105 // -- remember : using a replica, we have "o 117 // -- remember : using a replica, we have "one" physical volume
106 // -- but representing many different placem 118 // -- but representing many different placements, distinguished
107 // -- by replica number. By checking we are 119 // -- by replica number. By checking we are in the same physical
108 // -- volume, we verify the end step point h 120 // -- volume, we verify the end step point has not reached the
109 // -- world volume of the parallel geometry: 121 // -- world volume of the parallel geometry:
110 if (fPreStepTouchableHistory->GetVolume() != << 122 if ( fPreStepTouchableHistory ->GetVolume() !=
111 // -- the track is leaving the volumes in << 123 fPostStepTouchableHistory->GetVolume() )
112 // -- we leave this track unchanged: << 124 {
113 fDummyParticleChange.Initialize(*track); << 125 // -- the track is leaving the volumes in which biasing is applied,
114 return &fDummyParticleChange; << 126 // -- we leave this track unchanged:
115 } << 127 fDummyParticleChange.Initialize( *track );
>> 128 return &fDummyParticleChange;
>> 129 }
116 130
117 // ----------------------------------------- 131 // -------------------------------------------------------------------------------------
118 // -- At this stage, we know we have a parti 132 // -- At this stage, we know we have a particle crossing a boundary between two slices,
119 // -- each of this slice has a well defined 133 // -- each of this slice has a well defined importance : we apply the biasing.
120 // -- We will split if the track is entering 134 // -- We will split if the track is entering a volume with higher importance, and
121 // -- kill (applying Russian roulette) in th 135 // -- kill (applying Russian roulette) in the other case.
122 // ----------------------------------------- 136 // -------------------------------------------------------------------------------------
123 << 137
124 // -- We start by getting the replica number << 138 // -- We start by getting the replica numbers:
125 G4int preReplicaNumber = fPreStepTouchableHi << 139 G4int preReplicaNumber = fPreStepTouchableHistory->GetReplicaNumber();
126 G4int postReplicaNumber = fPostStepTouchable 140 G4int postReplicaNumber = fPostStepTouchableHistory->GetReplicaNumber();
127 141
128 // -- and get the related importance we defi 142 // -- and get the related importance we defined in the importance map:
129 G4int preImportance = (*fImportanceMap)[preR << 143 G4int preImportance = (*fImportanceMap)[ preReplicaNumber];
130 G4int postImportance = (*fImportanceMap)[pos 144 G4int postImportance = (*fImportanceMap)[postReplicaNumber];
131 145
>> 146
132 // -- Get track weight: 147 // -- Get track weight:
133 G4double initialWeight = track->GetWeight(); 148 G4double initialWeight = track->GetWeight();
134 << 149
135 // -- Initialize the "particle change" (it w 150 // -- Initialize the "particle change" (it will communicate the new track state to
136 // -- the tracking): 151 // -- the tracking):
137 fParticleChange.Initialize(*track); 152 fParticleChange.Initialize(*track);
138 153
139 if (postImportance > preImportance) { << 154 if ( postImportance > preImportance )
140 // -- We split : << 155 {
141 G4int splittingFactor = postImportance / p << 156 // -- We split :
142 << 157 G4int splittingFactor = postImportance/preImportance;
143 // Define the tracks weight: << 158
144 G4double weightOfTrack = initialWeight / s << 159 // Define the tracks weight:
145 << 160 G4double weightOfTrack = initialWeight/splittingFactor;
146 // Ask currect track weight to be changed << 161
147 fParticleChange.ProposeParentWeight(weight << 162 // Ask currect track weight to be changed to the new value:
148 // Now we clone this track (this is the ac << 163 fParticleChange.ProposeParentWeight( weightOfTrack );
149 // we will then have the primary and clone << 164 // Now we clone this track (this is the actual splitting):
150 // splitting by a factor 2: << 165 // we will then have the primary and clone of it, hence the
151 G4Track* clone = new G4Track(*track); << 166 // splitting by a factor 2:
152 clone->SetWeight(weightOfTrack); << 167 G4Track* clone = new G4Track( *track );
153 fParticleChange.AddSecondary(clone); << 168 clone->SetWeight( weightOfTrack );
154 // -- Below's call added for safety & illu << 169 fParticleChange.AddSecondary( clone );
155 // -- modify the clone (ie : daughter) wei << 170 // -- Below's call added for safety & illustration : inform particle change to not
156 // -- Here this call is not mandatory as b << 171 // -- modify the clone (ie : daughter) weight to make it that of the primary.
157 fParticleChange.SetSecondaryWeightByProces << 172 // -- Here this call is not mandatory as both tracks have same weights.
158 } << 173 fParticleChange.SetSecondaryWeightByProcess(true);
159 else { <<
160 // -- We apply Russian roulette: <<
161 G4double survivingProbability = G4double(p <<
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(fStop <<
168 } 174 }
169 else { << 175 else
170 // In this case, the track survives. We << 176 {
171 // to conserve weight among killed and s << 177 // -- We apply Russian roulette:
172 fParticleChange.ProposeParentWeight(init << 178 G4double survivingProbability = G4double(postImportance) / G4double(preImportance);
>> 179
>> 180 // Shoot a random number (in ]0,1[ segment):
>> 181 G4double random = G4UniformRand();
>> 182 if ( random > survivingProbability )
>> 183 {
>> 184 // We ask for the the track to be killed:
>> 185 fParticleChange.ProposeTrackStatus(fStopAndKill);
>> 186 }
>> 187 else
>> 188 {
>> 189 // In this case, the track survives. We change its weight
>> 190 // to conserve weight among killed and survival tracks:
>> 191 fParticleChange.ProposeParentWeight( initialWeight/survivingProbability );
>> 192 }
173 } 193 }
174 } <<
175 194
176 return &fParticleChange; 195 return &fParticleChange;
>> 196
177 } 197 }
178 198
179 //....oooOO0OOooo........oooOO0OOooo........oo 199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180 200