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 // ------------------------------------------- 27 // ------------------------------------------------------------------- 28 // 28 // 29 // GEANT4 Class header file 29 // GEANT4 Class header file 30 // 30 // 31 // 31 // 32 // File name: G4EmBiasingManager 32 // File name: G4EmBiasingManager 33 // 33 // 34 // Author: Vladimir Ivanchenko 34 // Author: Vladimir Ivanchenko 35 // 35 // 36 // Creation date: 28.07.2011 36 // Creation date: 28.07.2011 37 // 37 // 38 // Modifications: 38 // Modifications: 39 // 39 // 40 // Class Description: 40 // Class Description: 41 // 41 // 42 // It is a class providing step limit for forc 42 // It is a class providing step limit for forced process biasing 43 43 44 // ------------------------------------------- 44 // ------------------------------------------------------------------- 45 // 45 // 46 46 47 #ifndef G4EmBiasingManager_h 47 #ifndef G4EmBiasingManager_h 48 #define G4EmBiasingManager_h 1 48 #define G4EmBiasingManager_h 1 49 49 50 #include "globals.hh" 50 #include "globals.hh" 51 #include "G4ParticleDefinition.hh" 51 #include "G4ParticleDefinition.hh" 52 #include "G4DynamicParticle.hh" 52 #include "G4DynamicParticle.hh" 53 #include "Randomize.hh" 53 #include "Randomize.hh" 54 #include <vector> 54 #include <vector> 55 55 56 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 57 57 58 class G4Region; 58 class G4Region; 59 class G4Track; 59 class G4Track; 60 class G4VEnergyLossProcess; 60 class G4VEnergyLossProcess; 61 class G4VEmModel; 61 class G4VEmModel; 62 class G4MaterialCutsCouple; 62 class G4MaterialCutsCouple; 63 class G4ParticleChangeForLoss; 63 class G4ParticleChangeForLoss; 64 class G4ParticleChangeForGamma; 64 class G4ParticleChangeForGamma; 65 65 66 class G4EmBiasingManager 66 class G4EmBiasingManager 67 { 67 { 68 public: 68 public: 69 69 70 G4EmBiasingManager(); 70 G4EmBiasingManager(); 71 71 72 ~G4EmBiasingManager(); 72 ~G4EmBiasingManager(); 73 73 74 void Initialise(const G4ParticleDefinition& 74 void Initialise(const G4ParticleDefinition& part, 75 const G4String& procName, G4int verbose) 75 const G4String& procName, G4int verbose); 76 76 77 // default parameters are possible 77 // default parameters are possible 78 void ActivateForcedInteraction(G4double leng 78 void ActivateForcedInteraction(G4double length = 0.0, 79 const G4String& r = ""); 79 const G4String& r = ""); 80 80 81 // no default parameters 81 // no default parameters 82 void ActivateSecondaryBiasing(const G4String 82 void ActivateSecondaryBiasing(const G4String& region, G4double factor, 83 G4double energyLimit); 83 G4double energyLimit); 84 84 85 // return forced step limit 85 // return forced step limit 86 G4double GetStepLimit(G4int coupleIdx, G4dou 86 G4double GetStepLimit(G4int coupleIdx, G4double previousStep); 87 87 88 // return weight of splitting or Russian rou 88 // return weight of splitting or Russian roulette 89 // G4DynamicParticle may be deleted 89 // G4DynamicParticle may be deleted 90 // two functions are required because of the 90 // two functions are required because of the different ParticleChange 91 // ApplySecondaryBiasing() are wrappers 91 // ApplySecondaryBiasing() are wrappers 92 92 93 // for G4VEmProcess 93 // for G4VEmProcess 94 G4double ApplySecondaryBiasing(std::vector<G 94 G4double ApplySecondaryBiasing(std::vector<G4DynamicParticle*>&, 95 const G4Track& track, 95 const G4Track& track, 96 G4VEmModel* currentModel, 96 G4VEmModel* currentModel, 97 G4ParticleChangeForGamma* pParticleCh 97 G4ParticleChangeForGamma* pParticleChange, 98 G4double& eloss, 98 G4double& eloss, 99 G4int coupleIdx, 99 G4int coupleIdx, 100 G4double tcut, 100 G4double tcut, 101 G4double safety = 0.0); 101 G4double safety = 0.0); 102 102 103 // for G4VEnergyLossProcess 103 // for G4VEnergyLossProcess 104 G4double ApplySecondaryBiasing(std::vector<G 104 G4double ApplySecondaryBiasing(std::vector<G4DynamicParticle*>&, 105 const G4Track& track, 105 const G4Track& track, 106 G4VEmModel* currentModel, 106 G4VEmModel* currentModel, 107 G4ParticleChangeForLoss* pParticleCha 107 G4ParticleChangeForLoss* pParticleChange, 108 G4double& eloss, 108 G4double& eloss, 109 G4int coupleIdx, << 109 G4int coupleIdx, 110 G4double tcut, 110 G4double tcut, 111 G4double safety = 0.0); 111 G4double safety = 0.0); 112 112 113 // for G4VEnergyLossProcess 113 // for G4VEnergyLossProcess 114 G4double ApplySecondaryBiasing(std::vector<G 114 G4double ApplySecondaryBiasing(std::vector<G4Track*>&, 115 G4int coupleIdx); 115 G4int coupleIdx); 116 116 117 inline G4bool SecondaryBiasingRegion(G4int c 117 inline G4bool SecondaryBiasingRegion(G4int coupleIdx); 118 118 119 inline G4bool ForcedInteractionRegion(G4int 119 inline G4bool ForcedInteractionRegion(G4int coupleIdx); 120 120 121 inline void ResetForcedInteraction(); 121 inline void ResetForcedInteraction(); 122 122 123 G4bool CheckDirection(G4ThreeVector pos, G4T 123 G4bool CheckDirection(G4ThreeVector pos, G4ThreeVector momdir) const; 124 124 125 G4bool GetDirectionalSplitting() { return f 125 G4bool GetDirectionalSplitting() { return fDirectionalSplitting; } 126 void SetDirectionalSplitting(G4bool v) { 126 void SetDirectionalSplitting(G4bool v) { fDirectionalSplitting = v; } 127 127 128 void SetDirectionalSplittingTarget(G4ThreeVe 128 void SetDirectionalSplittingTarget(G4ThreeVector v) 129 { fDirectionalSplittingTarget = v; } 129 { fDirectionalSplittingTarget = v; } 130 void SetDirectionalSplittingRadius(G4double 130 void SetDirectionalSplittingRadius(G4double r) 131 { fDirectionalSplittingRadius = r; } 131 { fDirectionalSplittingRadius = r; } 132 G4double GetWeight(G4int i); 132 G4double GetWeight(G4int i); 133 133 134 // hide copy constructor and assignment oper 134 // hide copy constructor and assignment operator 135 G4EmBiasingManager(G4EmBiasingManager &) = d 135 G4EmBiasingManager(G4EmBiasingManager &) = delete; 136 G4EmBiasingManager & operator=(const G4EmBia 136 G4EmBiasingManager & operator=(const G4EmBiasingManager &right) = delete; 137 137 138 private: 138 private: 139 139 140 void ApplyRangeCut(std::vector<G4DynamicPart 140 void ApplyRangeCut(std::vector<G4DynamicParticle*>& vd, 141 const G4Track& track, 141 const G4Track& track, 142 G4double& eloss, 142 G4double& eloss, 143 G4double safety); 143 G4double safety); 144 144 145 G4double ApplySplitting(std::vector<G4Dynami 145 G4double ApplySplitting(std::vector<G4DynamicParticle*>& vd, 146 const G4Track& track, 146 const G4Track& track, 147 G4VEmModel* currentModel, 147 G4VEmModel* currentModel, 148 G4int index, 148 G4int index, 149 G4double tcut); 149 G4double tcut); 150 150 151 G4double ApplyDirectionalSplitting(std::vect 151 G4double ApplyDirectionalSplitting(std::vector<G4DynamicParticle*>& vd, 152 const G4Track& track, 152 const G4Track& track, 153 G4VEmModel* currentModel, 153 G4VEmModel* currentModel, 154 G4int index, 154 G4int index, 155 G4double tcut, 155 G4double tcut, 156 G4ParticleChangeForGamma* partChange); 156 G4ParticleChangeForGamma* partChange); 157 157 158 G4double ApplyDirectionalSplitting(std::vect 158 G4double ApplyDirectionalSplitting(std::vector<G4DynamicParticle*>& vd, 159 const G4Track& track, 159 const G4Track& track, 160 G4VEmModel* currentModel, 160 G4VEmModel* currentModel, 161 G4int index, 161 G4int index, 162 G4double tcut); 162 G4double tcut); 163 163 164 inline G4double ApplyRussianRoulette(std::ve 164 inline G4double ApplyRussianRoulette(std::vector<G4DynamicParticle*>& vd, 165 G4int index); 165 G4int index); 166 166 167 G4VEnergyLossProcess* eIonisation = nullptr; 167 G4VEnergyLossProcess* eIonisation = nullptr; 168 168 169 const G4ParticleDefinition* theElectron; 169 const G4ParticleDefinition* theElectron; 170 const G4ParticleDefinition* theGamma; 170 const G4ParticleDefinition* theGamma; 171 171 172 G4double fSafetyMin; 172 G4double fSafetyMin; 173 G4double currentStepLimit = 0.0; 173 G4double currentStepLimit = 0.0; 174 G4double fDirectionalSplittingRadius = 0.0; 174 G4double fDirectionalSplittingRadius = 0.0; 175 175 176 G4int nForcedRegions = 0; 176 G4int nForcedRegions = 0; 177 G4int nSecBiasedRegions = 0; 177 G4int nSecBiasedRegions = 0; 178 178 179 G4bool startTracking = true; 179 G4bool startTracking = true; 180 G4bool fDirectionalSplitting = false; 180 G4bool fDirectionalSplitting = false; 181 181 182 G4ThreeVector fDirectionalSplittingTarget; 182 G4ThreeVector fDirectionalSplittingTarget; 183 std::vector<G4double> fDirectionalSplittingW 183 std::vector<G4double> fDirectionalSplittingWeights; 184 184 185 std::vector<G4double> lengthForRegion 185 std::vector<G4double> lengthForRegion; 186 std::vector<G4double> secBiasedWeight 186 std::vector<G4double> secBiasedWeight; 187 std::vector<G4double> secBiasedEnegry 187 std::vector<G4double> secBiasedEnegryLimit; 188 188 189 std::vector<const G4Region*> forcedRegions; 189 std::vector<const G4Region*> forcedRegions; 190 std::vector<const G4Region*> secBiasedRegion 190 std::vector<const G4Region*> secBiasedRegions; 191 191 192 std::vector<G4int> nBremSplitting; 192 std::vector<G4int> nBremSplitting; 193 std::vector<G4int> idxForcedCouple 193 std::vector<G4int> idxForcedCouple; 194 std::vector<G4int> idxSecBiasedCou 194 std::vector<G4int> idxSecBiasedCouple; 195 195 196 std::vector<G4DynamicParticle*> tmpSecondari 196 std::vector<G4DynamicParticle*> tmpSecondaries; 197 }; 197 }; 198 198 199 inline G4bool 199 inline G4bool 200 G4EmBiasingManager::SecondaryBiasingRegion(G4i 200 G4EmBiasingManager::SecondaryBiasingRegion(G4int coupleIdx) 201 { 201 { 202 G4bool res = false; 202 G4bool res = false; 203 if(nSecBiasedRegions > 0) { 203 if(nSecBiasedRegions > 0) { 204 if(idxSecBiasedCouple[coupleIdx] >= 0) { r 204 if(idxSecBiasedCouple[coupleIdx] >= 0) { res = true; } 205 } 205 } 206 return res; 206 return res; 207 } 207 } 208 208 209 inline G4bool G4EmBiasingManager::ForcedIntera 209 inline G4bool G4EmBiasingManager::ForcedInteractionRegion(G4int coupleIdx) 210 { 210 { 211 G4bool res = false; 211 G4bool res = false; 212 if(nForcedRegions > 0) { 212 if(nForcedRegions > 0) { 213 if(idxForcedCouple[coupleIdx] >= 0) { res 213 if(idxForcedCouple[coupleIdx] >= 0) { res = true; } 214 } 214 } 215 return res; 215 return res; 216 } 216 } 217 217 218 inline void G4EmBiasingManager::ResetForcedInt 218 inline void G4EmBiasingManager::ResetForcedInteraction() 219 { 219 { 220 startTracking = true; 220 startTracking = true; 221 } 221 } 222 222 223 //....oooOO0OOooo........oooOO0OOooo........oo 223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 224 224 225 inline G4double 225 inline G4double 226 G4EmBiasingManager::ApplyRussianRoulette(std:: 226 G4EmBiasingManager::ApplyRussianRoulette(std::vector<G4DynamicParticle*>& vd, 227 G4int index) 227 G4int index) 228 { 228 { 229 size_t n = vd.size(); 229 size_t n = vd.size(); 230 G4double weight = secBiasedWeight[index]; 230 G4double weight = secBiasedWeight[index]; 231 for(size_t k=0; k<n; ++k) { 231 for(size_t k=0; k<n; ++k) { 232 if(G4UniformRand()*weight > 1.0) { 232 if(G4UniformRand()*weight > 1.0) { 233 const G4DynamicParticle* dp = vd[k]; 233 const G4DynamicParticle* dp = vd[k]; 234 delete dp; 234 delete dp; 235 vd[k] = nullptr; 235 vd[k] = nullptr; 236 } 236 } 237 } 237 } 238 return weight; 238 return weight; 239 } 239 } 240 240 241 //....oooOO0OOooo........oooOO0OOooo........oo 241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 242 242 243 #endif 243 #endif 244 244