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