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: G4VMscModel.hh 93264 2015-10-14 09:30:04Z gcosmo $ >> 27 // 26 // ------------------------------------------- 28 // ------------------------------------------------------------------- 27 // 29 // 28 // GEANT4 Class header file 30 // GEANT4 Class header file 29 // 31 // 30 // 32 // 31 // File name: G4VMscModel 33 // File name: G4VMscModel 32 // 34 // 33 // Author: Vladimir Ivanchenko 35 // Author: Vladimir Ivanchenko 34 // 36 // 35 // Creation date: 07.03.2008 37 // Creation date: 07.03.2008 36 // 38 // 37 // Modifications: 39 // Modifications: 38 // 07.04.2009 V.Ivanchenko moved msc methods f 40 // 07.04.2009 V.Ivanchenko moved msc methods from G4VEmModel to G4VMscModel 39 // 26.03.2012 V.Ivanchenko added transport x-s 41 // 26.03.2012 V.Ivanchenko added transport x-section pointer and Get?Set methods 40 // 42 // 41 // Class Description: 43 // Class Description: 42 // 44 // 43 // General interface to msc models 45 // General interface to msc models 44 46 45 // ------------------------------------------- 47 // ------------------------------------------------------------------- 46 // 48 // 47 #ifndef G4VMscModel_h 49 #ifndef G4VMscModel_h 48 #define G4VMscModel_h 1 50 #define G4VMscModel_h 1 49 51 50 #include <CLHEP/Units/SystemOfUnits.h> 52 #include <CLHEP/Units/SystemOfUnits.h> 51 53 52 #include "G4VEmModel.hh" 54 #include "G4VEmModel.hh" 53 #include "G4MscStepLimitType.hh" 55 #include "G4MscStepLimitType.hh" 54 #include "globals.hh" 56 #include "globals.hh" 55 #include "G4ThreeVector.hh" 57 #include "G4ThreeVector.hh" 56 #include "G4Track.hh" 58 #include "G4Track.hh" 57 #include "G4SafetyHelper.hh" 59 #include "G4SafetyHelper.hh" >> 60 #include "G4VEnergyLossProcess.hh" >> 61 #include "G4LossTableManager.hh" 58 #include "G4PhysicsTable.hh" 62 #include "G4PhysicsTable.hh" 59 #include "G4ThreeVector.hh" 63 #include "G4ThreeVector.hh" 60 #include <vector> 64 #include <vector> 61 65 62 class G4ParticleChangeForMSC; 66 class G4ParticleChangeForMSC; 63 class G4ParticleDefinition; << 64 class G4VEnergyLossProcess; << 65 67 66 class G4VMscModel : public G4VEmModel 68 class G4VMscModel : public G4VEmModel 67 { 69 { 68 70 69 public: 71 public: 70 72 71 explicit G4VMscModel(const G4String& nam); << 73 G4VMscModel(const G4String& nam); 72 74 73 ~G4VMscModel() override; << 75 virtual ~G4VMscModel(); 74 76 75 virtual G4double ComputeTruePathLengthLimit( 77 virtual G4double ComputeTruePathLengthLimit(const G4Track& track, 76 G4double& stepLimit) = 0; << 78 G4double& stepLimit); 77 79 78 virtual G4double ComputeGeomPathLength(G4dou << 80 virtual G4double ComputeGeomPathLength(G4double truePathLength); 79 81 80 virtual G4double ComputeTrueStepLength(G4dou << 82 virtual G4double ComputeTrueStepLength(G4double geomPathLength); 81 83 82 virtual G4ThreeVector& SampleScattering(cons 84 virtual G4ThreeVector& SampleScattering(const G4ThreeVector&, 83 G4double safety) = 0; << 85 G4double safety); 84 << 85 void InitialiseParameters(const G4ParticleDe << 86 << 87 void DumpParameters(std::ostream& out) const << 88 86 89 // empty method 87 // empty method 90 void SampleSecondaries(std::vector<G4Dynamic << 88 virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*, 91 const G4MaterialCutsCouple*, << 89 const G4MaterialCutsCouple*, 92 const G4DynamicParticle*, << 90 const G4DynamicParticle*, 93 G4double tmin, G4double tmax) override; << 91 G4double tmin, >> 92 G4double tmax); 94 93 95 //========================================== 94 //================================================================ 96 // Set parameters of multiple scattering mo 95 // Set parameters of multiple scattering models 97 //========================================== 96 //================================================================ 98 97 99 inline void SetStepLimitType(G4MscStepLimitT 98 inline void SetStepLimitType(G4MscStepLimitType); 100 99 101 inline void SetLateralDisplasmentFlag(G4bool 100 inline void SetLateralDisplasmentFlag(G4bool val); 102 101 103 inline void SetRangeFactor(G4double); 102 inline void SetRangeFactor(G4double); 104 103 105 inline void SetGeomFactor(G4double); 104 inline void SetGeomFactor(G4double); 106 105 107 inline void SetSkin(G4double); 106 inline void SetSkin(G4double); 108 107 109 inline void SetLambdaLimit(G4double); << 110 << 111 inline void SetSafetyFactor(G4double); << 112 << 113 inline void SetSampleZ(G4bool); 108 inline void SetSampleZ(G4bool); 114 109 115 //========================================== 110 //================================================================ 116 // Get/Set access to Physics Tables 111 // Get/Set access to Physics Tables 117 //========================================== 112 //================================================================ 118 113 119 inline G4VEnergyLossProcess* GetIonisation() 114 inline G4VEnergyLossProcess* GetIonisation() const; 120 115 121 inline void SetIonisation(G4VEnergyLossProce 116 inline void SetIonisation(G4VEnergyLossProcess*, 122 const G4ParticleDefinition* part); 117 const G4ParticleDefinition* part); 123 118 124 //========================================== 119 //================================================================ 125 // Run time methods 120 // Run time methods 126 //========================================== 121 //================================================================ 127 122 128 protected: 123 protected: 129 124 130 // initialisation of the ParticleChange for 125 // initialisation of the ParticleChange for the model 131 // initialisation of interface with geometry 126 // initialisation of interface with geometry and ionisation 132 G4ParticleChangeForMSC* 127 G4ParticleChangeForMSC* 133 GetParticleChangeForMSC(const G4ParticleDefi << 128 GetParticleChangeForMSC(const G4ParticleDefinition* p = 0); 134 129 135 // convert true length to geometry length 130 // convert true length to geometry length 136 inline G4double ConvertTrueToGeom(G4double& 131 inline G4double ConvertTrueToGeom(G4double& tLength, G4double& gLength); 137 132 138 // should be set before initialisation << 139 inline void SetUseSplineForMSC(G4bool val); << 140 << 141 public: 133 public: 142 134 143 // compute safety 135 // compute safety 144 inline G4double ComputeSafety(const G4ThreeV 136 inline G4double ComputeSafety(const G4ThreeVector& position, 145 G4double limit= DBL_MAX); 137 G4double limit= DBL_MAX); 146 138 147 // compute linear distance to a geometry bou 139 // compute linear distance to a geometry boundary 148 inline G4double ComputeGeomLimit(const G4Tra 140 inline G4double ComputeGeomLimit(const G4Track&, G4double& presafety, 149 G4double limit); 141 G4double limit); 150 142 151 G4double GetDEDX(const G4ParticleDefinition* << 143 inline G4double GetDEDX(const G4ParticleDefinition* part, 152 G4double kineticEner << 144 G4double kineticEnergy, 153 const G4MaterialCuts << 145 const G4MaterialCutsCouple* couple); 154 << 155 G4double GetDEDX(const G4ParticleDefinition* << 156 G4double kineticEner << 157 const G4MaterialCuts << 158 G4double logKineticE << 159 << 160 G4double GetRange(const G4ParticleDefinition << 161 G4double kineticEne << 162 const G4MaterialCut << 163 146 164 G4double GetRange(const G4ParticleDefinition << 147 inline G4double GetRange(const G4ParticleDefinition* part, 165 G4double kineticEne 148 G4double kineticEnergy, 166 const G4MaterialCut << 149 const G4MaterialCutsCouple* couple); 167 G4double logKinetic << 168 150 169 G4double GetEnergy(const G4ParticleDefinitio << 151 inline G4double GetEnergy(const G4ParticleDefinition* part, 170 G4double range, 152 G4double range, 171 const G4MaterialCutsCouple* couple); 153 const G4MaterialCutsCouple* couple); 172 154 173 // G4MaterialCutsCouple should be defined be 155 // G4MaterialCutsCouple should be defined before call to this method 174 inline 156 inline 175 G4double GetTransportMeanFreePath(const G4Pa 157 G4double GetTransportMeanFreePath(const G4ParticleDefinition* part, 176 G4double k << 158 G4double kinEnergy); 177 159 178 inline << 160 private: 179 G4double GetTransportMeanFreePath(const G4Pa << 180 G4double k << 181 G4double l << 182 161 183 // hide assignment operator 162 // hide assignment operator 184 G4VMscModel & operator=(const G4VMscModel & << 163 G4VMscModel & operator=(const G4VMscModel &right); 185 G4VMscModel(const G4VMscModel&) = delete; << 164 G4VMscModel(const G4VMscModel&); 186 << 187 private: << 188 165 189 G4SafetyHelper* safetyHelper = nullptr; << 166 G4SafetyHelper* safetyHelper; 190 G4VEnergyLossProcess* ionisation = nullptr; << 167 G4VEnergyLossProcess* ionisation; 191 const G4ParticleDefinition* currentPart = nu << 168 const G4ParticleDefinition* currentPart; 192 << 169 G4LossTableManager* man; 193 G4double dedx = 0.0; << 170 194 G4double localtkin = 0.0; << 171 G4double dedx; 195 G4double localrange = DBL_MAX; << 172 G4double localtkin; >> 173 G4double localrange; 196 174 197 protected: 175 protected: 198 176 199 G4double facrange = 0.04; << 177 G4double facrange; 200 G4double facgeom = 2.5; << 178 G4double facgeom; 201 G4double facsafety = 0.6; << 179 G4double facsafety; 202 G4double skin = 1.0; << 180 G4double skin; 203 G4double dtrl = 0.05; << 181 G4double dtrl; 204 G4double lambdalimit; 182 G4double lambdalimit; 205 G4double geomMin; 183 G4double geomMin; 206 G4double geomMax; 184 G4double geomMax; 207 185 208 G4ThreeVector fDisplacement; << 186 G4ThreeVector fDisplacement; 209 G4MscStepLimitType steppingAlgorithm; 187 G4MscStepLimitType steppingAlgorithm; 210 188 211 G4bool samplez = false; << 189 G4bool samplez; 212 G4bool latDisplasment = true; << 190 G4bool latDisplasment; 213 << 214 private: << 215 191 216 G4bool useSpline = true; << 217 }; 192 }; 218 193 219 //....oooOO0OOooo........oooOO0OOooo........oo 194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 220 //....oooOO0OOooo........oooOO0OOooo........oo 195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 221 196 222 inline void G4VMscModel::SetLateralDisplasment 197 inline void G4VMscModel::SetLateralDisplasmentFlag(G4bool val) 223 { 198 { 224 if(!IsLocked()) { latDisplasment = val; } 199 if(!IsLocked()) { latDisplasment = val; } 225 } 200 } 226 201 227 //....oooOO0OOooo........oooOO0OOooo........oo 202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 228 203 229 inline void G4VMscModel::SetSkin(G4double val) 204 inline void G4VMscModel::SetSkin(G4double val) 230 { 205 { 231 if(!IsLocked()) { skin = val; } 206 if(!IsLocked()) { skin = val; } 232 } 207 } 233 208 234 //....oooOO0OOooo........oooOO0OOooo........oo 209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 235 210 236 inline void G4VMscModel::SetRangeFactor(G4doub 211 inline void G4VMscModel::SetRangeFactor(G4double val) 237 { 212 { 238 if(!IsLocked()) { facrange = val; } 213 if(!IsLocked()) { facrange = val; } 239 } 214 } 240 215 241 //....oooOO0OOooo........oooOO0OOooo........oo 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 242 217 243 inline void G4VMscModel::SetGeomFactor(G4doubl 218 inline void G4VMscModel::SetGeomFactor(G4double val) 244 { 219 { 245 if(!IsLocked()) { facgeom = val; } 220 if(!IsLocked()) { facgeom = val; } 246 } 221 } 247 222 248 //....oooOO0OOooo........oooOO0OOooo........oo 223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 249 224 250 inline void G4VMscModel::SetLambdaLimit(G4doub << 251 { << 252 if(!IsLocked()) { lambdalimit = val; } << 253 } << 254 << 255 //....oooOO0OOooo........oooOO0OOooo........oo << 256 << 257 inline void G4VMscModel::SetSafetyFactor(G4dou << 258 { << 259 if(!IsLocked()) { facsafety = val; } << 260 } << 261 << 262 //....oooOO0OOooo........oooOO0OOooo........oo << 263 << 264 inline void G4VMscModel::SetStepLimitType(G4Ms 225 inline void G4VMscModel::SetStepLimitType(G4MscStepLimitType val) 265 { 226 { 266 if(!IsLocked()) { steppingAlgorithm = val; } 227 if(!IsLocked()) { steppingAlgorithm = val; } 267 } 228 } 268 229 269 //....oooOO0OOooo........oooOO0OOooo........oo 230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 270 231 271 inline void G4VMscModel::SetSampleZ(G4bool val 232 inline void G4VMscModel::SetSampleZ(G4bool val) 272 { 233 { 273 if(!IsLocked()) { samplez = val; } 234 if(!IsLocked()) { samplez = val; } 274 } 235 } 275 236 276 //....oooOO0OOooo........oooOO0OOooo........oo 237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 277 238 278 inline G4double G4VMscModel::ComputeSafety(con 239 inline G4double G4VMscModel::ComputeSafety(const G4ThreeVector& position, 279 G4double limit) 240 G4double limit) 280 { 241 { 281 return safetyHelper->ComputeSafety(position 242 return safetyHelper->ComputeSafety(position, limit); 282 } 243 } 283 244 284 //....oooOO0OOooo........oooOO0OOooo........oo 245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 285 246 286 inline G4double G4VMscModel::ConvertTrueToGeom 247 inline G4double G4VMscModel::ConvertTrueToGeom(G4double& tlength, 287 G4double& glength) 248 G4double& glength) 288 { 249 { 289 glength = ComputeGeomPathLength(tlength); 250 glength = ComputeGeomPathLength(tlength); 290 // should return true length 251 // should return true length 291 return tlength; 252 return tlength; 292 } 253 } 293 254 294 //....oooOO0OOooo........oooOO0OOooo........oo 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 295 256 296 inline G4double G4VMscModel::ComputeGeomLimit( 257 inline G4double G4VMscModel::ComputeGeomLimit(const G4Track& track, 297 G4double& presafety, 258 G4double& presafety, 298 G4double limit) 259 G4double limit) 299 { 260 { >> 261 /* >> 262 G4double res = geomMax; >> 263 if(track.GetVolume() != safetyHelper->GetWorldVolume()) { >> 264 G4double res = safetyHelper->CheckNextStep( >> 265 track.GetStep()->GetPreStepPoint()->GetPosition(), >> 266 track.GetMomentumDirection(), >> 267 limit, presafety); >> 268 } >> 269 return res; >> 270 */ 300 return safetyHelper->CheckNextStep( 271 return safetyHelper->CheckNextStep( 301 track.GetStep()->GetPreStepPoint()-> 272 track.GetStep()->GetPreStepPoint()->GetPosition(), 302 track.GetMomentumDirection(), 273 track.GetMomentumDirection(), 303 limit, presafety); 274 limit, presafety); 304 } 275 } 305 276 306 //....oooOO0OOooo........oooOO0OOooo........oo 277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 307 278 308 inline G4VEnergyLossProcess* G4VMscModel::GetI << 279 inline G4double >> 280 G4VMscModel::GetDEDX(const G4ParticleDefinition* part, >> 281 G4double kinEnergy, const G4MaterialCutsCouple* couple) 309 { 282 { 310 return ionisation; << 283 G4double x; >> 284 if(ionisation) { x = ionisation->GetDEDX(kinEnergy, couple); } >> 285 else { >> 286 G4double q = part->GetPDGCharge()*inveplus; >> 287 x = dedx*q*q; >> 288 } >> 289 return x; 311 } 290 } 312 291 313 //....oooOO0OOooo........oooOO0OOooo........oo 292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 314 293 315 inline void G4VMscModel::SetIonisation(G4VEner << 294 inline G4double 316 const G4ParticleDefinition* par << 295 G4VMscModel::GetRange(const G4ParticleDefinition* part, >> 296 G4double kinEnergy, const G4MaterialCutsCouple* couple) 317 { 297 { 318 ionisation = p; << 298 //G4cout << "G4VMscModel::GetRange E(MeV)= " << kinEnergy << " " 319 currentPart = part; << 299 // << ionisation << " " << part->GetParticleName() >> 300 // << G4endl; >> 301 localtkin = kinEnergy; >> 302 if(ionisation) { >> 303 localrange = ionisation->GetRangeForLoss(kinEnergy, couple); >> 304 } else { >> 305 G4double q = part->GetPDGCharge()*inveplus; >> 306 localrange = kinEnergy/(dedx*q*q*couple->GetMaterial()->GetDensity()); >> 307 } >> 308 //G4cout << "R(mm)= " << localrange << " " << ionisation << G4endl; >> 309 return localrange; 320 } 310 } 321 311 322 //....oooOO0OOooo........oooOO0OOooo........oo 312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 323 313 324 inline G4double 314 inline G4double 325 G4VMscModel::GetTransportMeanFreePath(const G4 << 315 G4VMscModel::GetEnergy(const G4ParticleDefinition* part, 326 G4double << 316 G4double range, const G4MaterialCutsCouple* couple) 327 { 317 { 328 G4double x; << 318 G4double e; 329 if (nullptr != xSectionTable) { << 319 //G4cout << "G4VMscModel::GetEnergy R(mm)= " << range << " " << ionisation 330 x = pFactor*(*xSectionTable)[basedCoupleIn << 320 // << " Rlocal(mm)= " << localrange << " Elocal(MeV)= " << localtkin 331 } else { << 321 // << G4endl; 332 x = pFactor*CrossSectionPerVolume(pBaseMat << 322 if(ionisation) { e = ionisation->GetKineticEnergy(range, couple); } >> 323 else { >> 324 e = localtkin; >> 325 if(localrange > range) { >> 326 G4double q = part->GetPDGCharge()*inveplus; >> 327 e -= (localrange - range)*dedx*q*q*couple->GetMaterial()->GetDensity(); >> 328 } 333 } 329 } 334 return (x > 0.0) ? 1.0/x : DBL_MAX; << 330 return e; 335 } 331 } 336 332 337 //....oooOO0OOooo........oooOO0OOooo........oo 333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 338 334 339 inline G4double << 335 inline G4VEnergyLossProcess* G4VMscModel::GetIonisation() const 340 G4VMscModel::GetTransportMeanFreePath(const G4 << 341 G4double << 342 { 336 { 343 G4double x; << 337 return ionisation; 344 if (nullptr != xSectionTable) { << 338 } 345 x = pFactor*(*xSectionTable)[basedCoupleIn << 339 346 } else { << 340 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 347 x = pFactor*CrossSectionPerVolume(pBaseMat << 341 348 } << 342 inline void G4VMscModel::SetIonisation(G4VEnergyLossProcess* p, 349 return (x > 0.0) ? 1.0/x : DBL_MAX; << 343 const G4ParticleDefinition* part) >> 344 { >> 345 ionisation = p; >> 346 currentPart = part; 350 } 347 } 351 348 352 //....oooOO0OOooo........oooOO0OOooo........oo 349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 353 350 354 inline void G4VMscModel::SetUseSplineForMSC(G4 << 351 inline G4double >> 352 G4VMscModel::GetTransportMeanFreePath(const G4ParticleDefinition* part, >> 353 G4double ekin) 355 { 354 { 356 useSpline = val; << 355 G4double x; >> 356 if(xSectionTable) { >> 357 G4int idx = CurrentCouple()->GetIndex(); >> 358 x = (*xSectionTable)[(*theDensityIdx)[idx]]->Value(ekin, idxTable) >> 359 *(*theDensityFactor)[idx]/(ekin*ekin); >> 360 } else { >> 361 x = CrossSectionPerVolume(CurrentCouple()->GetMaterial(), part, ekin, >> 362 0.0, DBL_MAX); >> 363 } >> 364 if(0.0 >= x) { x = DBL_MAX; } >> 365 else { x = 1.0/x; } >> 366 return x; 357 } 367 } 358 368 359 //....oooOO0OOooo........oooOO0OOooo........oo 369 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 360 370 361 #endif 371 #endif >> 372 362 373