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 << 135 G4EmBiasingManager(G4EmBiasingManager &) = d << 136 G4EmBiasingManager & operator=(const G4EmBia << 137 << 138 private: 134 private: 139 135 140 void ApplyRangeCut(std::vector<G4DynamicPart 136 void ApplyRangeCut(std::vector<G4DynamicParticle*>& vd, 141 const G4Track& track, 137 const G4Track& track, 142 G4double& eloss, 138 G4double& eloss, 143 G4double safety); 139 G4double safety); 144 140 145 G4double ApplySplitting(std::vector<G4Dynami 141 G4double ApplySplitting(std::vector<G4DynamicParticle*>& vd, 146 const G4Track& track, 142 const G4Track& track, 147 G4VEmModel* currentModel, 143 G4VEmModel* currentModel, 148 G4int index, 144 G4int index, 149 G4double tcut); 145 G4double tcut); 150 146 151 G4double ApplyDirectionalSplitting(std::vect << 147 void ApplyDirectionalSplitting(std::vector<G4DynamicParticle*>& vd, 152 const G4Track& track, 148 const G4Track& track, 153 G4VEmModel* currentModel, 149 G4VEmModel* currentModel, 154 G4int index, 150 G4int index, 155 G4double tcut, 151 G4double tcut, 156 G4ParticleChangeForGamma* partChange); 152 G4ParticleChangeForGamma* partChange); 157 153 158 G4double ApplyDirectionalSplitting(std::vect << 154 void ApplyDirectionalSplitting(std::vector<G4DynamicParticle*>& vd, 159 const G4Track& track, 155 const G4Track& track, 160 G4VEmModel* currentModel, 156 G4VEmModel* currentModel, 161 G4int index, 157 G4int index, 162 G4double tcut); 158 G4double tcut); 163 159 164 inline G4double ApplyRussianRoulette(std::ve 160 inline G4double ApplyRussianRoulette(std::vector<G4DynamicParticle*>& vd, 165 G4int index); 161 G4int index); 166 162 167 G4VEnergyLossProcess* eIonisation = nullptr; << 163 // hide copy constructor and assignment operator 168 << 164 G4EmBiasingManager(G4EmBiasingManager &) = delete; 169 const G4ParticleDefinition* theElectron; << 165 G4EmBiasingManager & operator=(const G4EmBiasingManager &right) = delete; 170 const G4ParticleDefinition* theGamma; << 171 << 172 G4double fSafetyMin; << 173 G4double currentStepLimit = 0.0; << 174 G4double fDirectionalSplittingRadius = 0.0; << 175 << 176 G4int nForcedRegions = 0; << 177 G4int nSecBiasedRegions = 0; << 178 << 179 G4bool startTracking = true; << 180 G4bool fDirectionalSplitting = false; << 181 << 182 G4ThreeVector fDirectionalSplittingTarget; << 183 std::vector<G4double> fDirectionalSplittingW << 184 166 >> 167 G4int nForcedRegions; >> 168 G4int nSecBiasedRegions; >> 169 std::vector<const G4Region*> forcedRegions; 185 std::vector<G4double> lengthForRegion 170 std::vector<G4double> lengthForRegion; >> 171 std::vector<const G4Region*> secBiasedRegions; 186 std::vector<G4double> secBiasedWeight 172 std::vector<G4double> secBiasedWeight; 187 std::vector<G4double> secBiasedEnegry 173 std::vector<G4double> secBiasedEnegryLimit; 188 << 189 std::vector<const G4Region*> forcedRegions; << 190 std::vector<const G4Region*> secBiasedRegion << 191 << 192 std::vector<G4int> nBremSplitting; 174 std::vector<G4int> nBremSplitting; >> 175 193 std::vector<G4int> idxForcedCouple 176 std::vector<G4int> idxForcedCouple; 194 std::vector<G4int> idxSecBiasedCou 177 std::vector<G4int> idxSecBiasedCouple; 195 178 196 std::vector<G4DynamicParticle*> tmpSecondari 179 std::vector<G4DynamicParticle*> tmpSecondaries; >> 180 >> 181 G4VEnergyLossProcess* eIonisation; >> 182 >> 183 const G4ParticleDefinition* theElectron; >> 184 const G4ParticleDefinition* theGamma; >> 185 >> 186 G4double fSafetyMin; >> 187 G4double currentStepLimit; >> 188 G4bool startTracking; >> 189 >> 190 G4double fWeight; >> 191 >> 192 G4bool fDirectionalSplitting; >> 193 G4ThreeVector fDirectionalSplittingTarget; >> 194 G4double fDirectionalSplittingRadius; >> 195 std::vector<G4double> fDirectionalSplittingWeights; 197 }; 196 }; 198 197 199 inline G4bool 198 inline G4bool 200 G4EmBiasingManager::SecondaryBiasingRegion(G4i 199 G4EmBiasingManager::SecondaryBiasingRegion(G4int coupleIdx) 201 { 200 { 202 G4bool res = false; 201 G4bool res = false; 203 if(nSecBiasedRegions > 0) { 202 if(nSecBiasedRegions > 0) { 204 if(idxSecBiasedCouple[coupleIdx] >= 0) { r 203 if(idxSecBiasedCouple[coupleIdx] >= 0) { res = true; } 205 } 204 } 206 return res; 205 return res; 207 } 206 } 208 207 209 inline G4bool G4EmBiasingManager::ForcedIntera 208 inline G4bool G4EmBiasingManager::ForcedInteractionRegion(G4int coupleIdx) 210 { 209 { 211 G4bool res = false; 210 G4bool res = false; 212 if(nForcedRegions > 0) { 211 if(nForcedRegions > 0) { 213 if(idxForcedCouple[coupleIdx] >= 0) { res 212 if(idxForcedCouple[coupleIdx] >= 0) { res = true; } 214 } 213 } 215 return res; 214 return res; 216 } 215 } 217 216 218 inline void G4EmBiasingManager::ResetForcedInt 217 inline void G4EmBiasingManager::ResetForcedInteraction() 219 { 218 { 220 startTracking = true; 219 startTracking = true; 221 } 220 } 222 221 223 //....oooOO0OOooo........oooOO0OOooo........oo 222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 224 223 225 inline G4double 224 inline G4double 226 G4EmBiasingManager::ApplyRussianRoulette(std:: 225 G4EmBiasingManager::ApplyRussianRoulette(std::vector<G4DynamicParticle*>& vd, 227 G4int index) 226 G4int index) 228 { 227 { 229 size_t n = vd.size(); 228 size_t n = vd.size(); 230 G4double weight = secBiasedWeight[index]; 229 G4double weight = secBiasedWeight[index]; 231 for(size_t k=0; k<n; ++k) { 230 for(size_t k=0; k<n; ++k) { 232 if(G4UniformRand()*weight > 1.0) { 231 if(G4UniformRand()*weight > 1.0) { 233 const G4DynamicParticle* dp = vd[k]; 232 const G4DynamicParticle* dp = vd[k]; 234 delete dp; 233 delete dp; 235 vd[k] = nullptr; 234 vd[k] = nullptr; 236 } 235 } 237 } 236 } 238 return weight; 237 return weight; 239 } 238 } 240 239 241 //....oooOO0OOooo........oooOO0OOooo........oo 240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 242 241 243 #endif 242 #endif 244 243