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 GB03BOptnSplitOrKillOnBoundary.cc 28 /// \file GB03BOptnSplitOrKillOnBoundary.cc 28 /// \brief Implementation of the GB03BOptnSpli 29 /// \brief Implementation of the GB03BOptnSplitOrKillOnBoundary class 29 30 30 #include "GB03BOptnSplitOrKillOnBoundary.hh" << 31 << 32 #include "Randomize.hh" 31 #include "Randomize.hh" >> 32 #include "GB03BOptnSplitOrKillOnBoundary.hh" 33 33 34 //....oooOO0OOooo........oooOO0OOooo........oo 34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 35 35 36 GB03BOptnSplitOrKillOnBoundary::GB03BOptnSplit 36 GB03BOptnSplitOrKillOnBoundary::GB03BOptnSplitOrKillOnBoundary(G4String name) 37 : G4VBiasingOperation(name), fParticleChange << 37 : G4VBiasingOperation(name), >> 38 fParticleChange(), >> 39 fParticleChangeForNothing() 38 {} 40 {} 39 41 40 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 41 43 42 GB03BOptnSplitOrKillOnBoundary::~GB03BOptnSpli << 44 GB03BOptnSplitOrKillOnBoundary::~GB03BOptnSplitOrKillOnBoundary() >> 45 {} 43 46 44 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 45 48 46 G4double GB03BOptnSplitOrKillOnBoundary::Dista << 49 G4double GB03BOptnSplitOrKillOnBoundary:: 47 << 50 DistanceToApplyOperation( const G4Track*, >> 51 G4double, >> 52 G4ForceCondition* condition) 48 { 53 { 49 // -- return "infinite" distance for interac 54 // -- return "infinite" distance for interaction, but asks for GenerateBiasingFinalState 50 // -- being called anyway at the end of the 55 // -- being called anyway at the end of the step, by returning the "Forced" condition 51 // -- flag. 56 // -- flag. 52 *condition = Forced; 57 *condition = Forced; 53 return DBL_MAX; 58 return DBL_MAX; 54 } 59 } 55 60 56 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 57 62 58 G4VParticleChange* GB03BOptnSplitOrKillOnBound << 63 G4VParticleChange* 59 << 64 GB03BOptnSplitOrKillOnBoundary::GenerateBiasingFinalState( const G4Track* track, >> 65 const G4Step* step ) 60 { 66 { >> 67 61 // Check if step is limited by the geometry: 68 // Check if step is limited by the geometry: as we attached the biasing operator 62 // to the absorber layer, this volume bounda 69 // to the absorber layer, this volume boundary is the one of the absorber. 63 // (check of current step # of track is inel 70 // (check of current step # of track is inelegant, but is to fix a "feature" 64 // that a cloned track can wrongly be seen i 71 // that a cloned track can wrongly be seen in the wrong volume, because of numerical 65 // precision issue. In this case it makes a 72 // precision issue. In this case it makes a tiny step, which should be disregarded). 66 if ((step->GetPostStepPoint()->GetStepStatus << 73 if ( ( step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary ) && 67 && (track->GetCurrentStepNumber() != 1)) << 74 ( track->GetCurrentStepNumber() != 1 ) ) 68 { << 75 { 69 // -- Before deciding for killing or split << 76 70 // -- the technique or not: << 77 // -- Before deciding for killing or splitting, we make decision on applying 71 G4double trial = G4UniformRand(); // -- N << 78 // -- the technique or not: 72 // -- e << 79 G4double trial = G4UniformRand(); // -- Note: G4UniformRand() is thread-safe 73 if (trial <= fApplyProbability) { << 80 // -- engine for random numbers 74 // -- Get z component of track, to see i << 81 if ( trial <= fApplyProbability ) 75 G4double pz = track->GetMomentum().z(); << 82 { 76 << 83 // -- Get z component of track, to see if it moves forward or backward: 77 if (pz > 0.0) { << 84 G4double pz = track->GetMomentum().z(); 78 // ----------------------------------- << 85 79 // Here, we are moving "forward". We d << 86 if ( pz > 0.0 ) 80 // ----------------------------------- << 87 { 81 << 88 // ------------------------------------------------- 82 // Get track weight: << 89 // Here, we are moving "forward". We do "splitting": 83 G4double initialWeight = track->GetWei << 90 // ------------------------------------------------- 84 // Define the tracks weight: << 91 85 G4double weightOfTrack = initialWeight << 92 // Get track weight: 86 << 93 G4double initialWeight = track->GetWeight(); 87 // The "particle change" is the object << 94 // Define the tracks weight: 88 // the tracking the update of the prim << 95 G4double weightOfTrack = initialWeight/fSplittingFactor; 89 // secondary tracks. << 96 90 fParticleChange.Initialize(*track); << 97 // The "particle change" is the object to be used to communicate to 91 << 98 // the tracking the update of the primary state and/or creation 92 // ask currect track weight to be chan << 99 // secondary tracks. 93 fParticleChange.ProposeParentWeight(we << 100 fParticleChange.Initialize(*track); 94 << 101 95 // Now make clones of this track (this << 102 // ask currect track weight to be changed to new value: 96 // we will then have the primary and N << 103 fParticleChange.ProposeParentWeight( weightOfTrack ); 97 // splitting by a factor N: << 104 98 fParticleChange.SetNumberOfSecondaries << 105 // Now make clones of this track (this is the actual splitting): 99 for (G4int iSplit = 1; iSplit < fSplit << 106 // we will then have the primary and N-1 clones of it, hence the 100 G4Track* clone = new G4Track(*track) << 107 // splitting by a factor N: 101 clone->SetWeight(weightOfTrack); << 108 fParticleChange.SetNumberOfSecondaries( fSplittingFactor-1 ); 102 fParticleChange.AddSecondary(clone); << 109 for ( G4int iSplit = 1 ; iSplit < fSplittingFactor ; iSplit++ ) 103 } << 110 { 104 fParticleChange.SetSecondaryWeightByPr << 111 G4Track* clone = new G4Track( *track ); 105 // -- take it as is ;) [though not nec << 112 clone->SetWeight( weightOfTrack ); 106 << 113 fParticleChange.AddSecondary( clone ); 107 // this new final state is returned to << 114 } 108 return &fParticleChange; << 115 fParticleChange.SetSecondaryWeightByProcess(true); // -- tricky 109 } << 116 // -- take it as is ;) [though not necessary here, put for safety] 110 << 117 111 else { << 118 // this new final state is returned to the tracking; 112 // ----------------------------------- << 119 return &fParticleChange; 113 // Here, we are moving backward. We do << 120 114 // roulette, killing 1/fSplittingFacto << 121 } 115 // ----------------------------------- << 122 116 << 123 else 117 // Get track weight: << 124 118 G4double initialWeight = track->GetWei << 125 { 119 << 126 // -------------------------------------------------------------- 120 // The "particle change" is the object << 127 // Here, we are moving backward. We do killing, playing a russian 121 // the tracking the update of the prim << 128 // roulette, killing 1/fSplittingFactor of the tracks in average: 122 // secondary tracks. << 129 // -------------------------------------------------------------- 123 fParticleChange.Initialize(*track); << 130 124 << 131 // Get track weight: 125 // Shoot a random number (in ]0,1[ seg << 132 G4double initialWeight = track->GetWeight(); 126 G4double random = G4UniformRand(); << 133 127 << 134 // The "particle change" is the object to be used to communicate to 128 // Decide to kill, keeping 1/fSplittin << 135 // the tracking the update of the primary state and/or creation 129 G4double survivingProbability = 1.0 / << 136 // secondary tracks. 130 if (random > survivingProbability) { << 137 fParticleChange.Initialize(*track); 131 // We ask for the the track to be ki << 138 132 fParticleChange.ProposeTrackStatus(f << 139 // Shoot a random number (in ]0,1[ segment): 133 } << 140 G4double random = G4UniformRand(); 134 else { << 141 135 // In this case, the track survives. << 142 // Decide to kill with 1/fSplittingFactor probability: 136 // to conserve weight among killed a << 143 G4double killingProbability = 1.0/fSplittingFactor; 137 fParticleChange.ProposeParentWeight( << 144 if ( random < killingProbability ) 138 } << 145 { 139 << 146 // We ask for the the track to be killed: 140 // this new final state is returned to << 147 fParticleChange.ProposeTrackStatus(fStopAndKill); 141 return &fParticleChange; << 148 } 142 } << 149 else 143 } // -- end of : if ( trial > probaForApp << 150 { 144 } // -- end of : if ( ( step->GetPostStepPo << 151 // In this case, the track survives. We change its weight 145 // << 152 // to conserve weight among killed and survival tracks: >> 153 fParticleChange.ProposeParentWeight( initialWeight*fSplittingFactor ); >> 154 } >> 155 >> 156 // this new final state is returned to the tracking; >> 157 return &fParticleChange; >> 158 } >> 159 } // -- end of : if ( trial > probaForApplying ) >> 160 } // -- end of : if ( ( step->GetPostStepPoint()->GetStepStatus() == >> 161 // fGeomBoundary ) ... >> 162 146 163 147 // Here, the step was not limited by the geo 164 // Here, the step was not limited by the geometry (but certainly by a physics 148 // process). We do nothing: ie we make no ch 165 // process). We do nothing: ie we make no change to the current track. 149 fParticleChangeForNothing.Initialize(*track) 166 fParticleChangeForNothing.Initialize(*track); 150 return &fParticleChangeForNothing; 167 return &fParticleChangeForNothing; >> 168 151 } 169 } 152 170 153 //....oooOO0OOooo........oooOO0OOooo........oo 171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 154 172