Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // >> 23 // $Id: G4VEnergyLossProcess.hh,v 1.34 2004/12/01 18:01:01 vnivanch Exp $ >> 24 // GEANT4 tag $Name: geant4-07-00-patch-01 $ 26 // 25 // 27 // ------------------------------------------- 26 // ------------------------------------------------------------------- 28 // 27 // 29 // GEANT4 Class header file 28 // GEANT4 Class header file 30 // 29 // 31 // 30 // 32 // File name: G4VEnergyLossProcess 31 // File name: G4VEnergyLossProcess 33 // 32 // 34 // Author: Vladimir Ivanchenko on base 33 // Author: Vladimir Ivanchenko on base of Laszlo Urban code 35 // 34 // 36 // Creation date: 03.01.2002 35 // Creation date: 03.01.2002 37 // 36 // 38 // Modifications: Vladimir Ivanchenko << 37 // Modifications: >> 38 // >> 39 // 26-12-02 Secondary production moved to derived classes (V.Ivanchenko) >> 40 // 20-01-03 Migrade to cut per region (V.Ivanchenko) >> 41 // 24-01-03 Make models region aware (V.Ivanchenko) >> 42 // 05-02-03 Fix compilation warnings (V.Ivanchenko) >> 43 // 13-02-03 SubCutoffProcessors defined for regions (V.Ivanchenko) >> 44 // 17-02-03 Fix problem of store/restore tables (V.Ivanchenko) >> 45 // 26-02-03 Region dependent step limit (V.Ivanchenko) >> 46 // 26-03-03 Add GetDEDXDispersion (V.Ivanchenko) >> 47 // 09-04-03 Fix problem of negative range limit for non integral (V.Ivanchenko) >> 48 // 13-05-03 Add calculation of precise range (V.Ivanchenko) >> 49 // 21-07-03 Add UpdateEmModel method (V.Ivanchenko) >> 50 // 12-11-03 G4EnergyLossSTD -> G4EnergyLossProcess (V.Ivanchenko) >> 51 // 14-01-04 Activate precise range calculation (V.Ivanchenko) >> 52 // 10-03-04 Fix problem of step limit calculation (V.Ivanchenko) >> 53 // 30-06-04 make destructor virtual (V.Ivanchenko) >> 54 // 05-07-04 fix problem of GenericIons seen at small cuts (V.Ivanchenko) >> 55 // 03-08-04 Add DEDX table to all processes for control on integral range(V.Ivanchenko) >> 56 // 06-08-04 Clear up names of member functions (V.Ivanchenko) >> 57 // 27-08-04 Add NeedBuildTables method (V.Ivanchneko) >> 58 // 09-09-04 Bug fix for the integral mode with 2 peaks (V.Ivanchneko) >> 59 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko) 39 // 60 // 40 // Class Description: 61 // Class Description: 41 // 62 // 42 // It is the unified energy loss process it ca 63 // It is the unified energy loss process it calculates the continuous 43 // energy loss for charged particles using a s 64 // energy loss for charged particles using a set of Energy Loss 44 // models valid for different energy regions. 65 // models valid for different energy regions. There are a possibility 45 // to create and access to dE/dx and range tab 66 // to create and access to dE/dx and range tables, or to calculate 46 // that information on fly. 67 // that information on fly. 47 68 48 // ------------------------------------------- 69 // ------------------------------------------------------------------- 49 // 70 // 50 71 51 #ifndef G4VEnergyLossProcess_h 72 #ifndef G4VEnergyLossProcess_h 52 #define G4VEnergyLossProcess_h 1 73 #define G4VEnergyLossProcess_h 1 53 74 54 #include "G4VContinuousDiscreteProcess.hh" 75 #include "G4VContinuousDiscreteProcess.hh" 55 #include "globals.hh" 76 #include "globals.hh" 56 #include "G4Material.hh" 77 #include "G4Material.hh" 57 #include "G4MaterialCutsCouple.hh" 78 #include "G4MaterialCutsCouple.hh" 58 #include "G4Track.hh" 79 #include "G4Track.hh" 59 #include "G4EmModelManager.hh" 80 #include "G4EmModelManager.hh" >> 81 #include "G4UnitsTable.hh" 60 #include "G4ParticleChangeForLoss.hh" 82 #include "G4ParticleChangeForLoss.hh" 61 #include "G4EmTableType.hh" << 62 #include "G4EmSecondaryParticleType.hh" << 63 #include "G4PhysicsTable.hh" << 64 #include "G4PhysicsVector.hh" << 65 83 66 class G4Step; 84 class G4Step; 67 class G4ParticleDefinition; 85 class G4ParticleDefinition; 68 class G4EmParameters; << 69 class G4VEmModel; 86 class G4VEmModel; 70 class G4VEmFluctuationModel; 87 class G4VEmFluctuationModel; 71 class G4DataVector; 88 class G4DataVector; >> 89 class G4PhysicsTable; >> 90 class G4PhysicsVector; >> 91 class G4VSubCutoffProcessor; 72 class G4Region; 92 class G4Region; 73 class G4SafetyHelper; << 74 class G4VAtomDeexcitation; << 75 class G4VSubCutProducer; << 76 class G4EmBiasingManager; << 77 class G4LossTableManager; << 78 class G4EmDataHandler; << 79 93 80 //....oooOO0OOooo........oooOO0OOooo........oo 94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 81 95 82 class G4VEnergyLossProcess : public G4VContinu 96 class G4VEnergyLossProcess : public G4VContinuousDiscreteProcess 83 { 97 { 84 public: 98 public: 85 99 86 G4VEnergyLossProcess(const G4String& name = 100 G4VEnergyLossProcess(const G4String& name = "EnergyLoss", 87 G4ProcessType type = fE << 101 G4ProcessType type = fElectromagnetic); 88 102 89 ~G4VEnergyLossProcess() override; << 103 virtual ~G4VEnergyLossProcess(); 90 104 91 //------------------------------------------ << 105 G4VParticleChange* AlongStepDoIt(const G4Track&, const G4Step&); 92 // Virtual methods to be implemented in conc << 93 //------------------------------------------ << 94 106 95 protected: << 107 G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&); 96 108 97 // description of specific process parameter << 109 virtual std::vector<G4Track*>* SecondariesAlongStep( 98 virtual void StreamProcessInfo(std::ostream& << 110 const G4Step&, >> 111 G4double& tmax, >> 112 G4double& eloss, >> 113 G4double& kinEnergy) = 0; 99 114 100 virtual void InitialiseEnergyLossProcess(con << 115 virtual void SecondariesPostStep( 101 con << 116 G4VEmModel*, >> 117 const G4MaterialCutsCouple*, >> 118 const G4DynamicParticle*, >> 119 G4double& tcut, >> 120 G4double& kinEnergy) = 0; 102 121 103 public: << 122 virtual G4bool IsApplicable(const G4ParticleDefinition& p) = 0; >> 123 // True for all charged particles 104 124 105 // used as low energy limit LambdaTable << 125 virtual void PreparePhysicsTable(const G4ParticleDefinition&); 106 virtual G4double MinPrimaryEnergy(const G4Pa << 126 // Initialise for build of tables 107 const G4Ma << 127 >> 128 virtual void BuildPhysicsTable(const G4ParticleDefinition&); >> 129 // Build physics table during initialisation 108 130 109 // print documentation in html format << 131 virtual void PrintInfoDefinition(); 110 void ProcessDescription(std::ostream& outFil << 132 // Print out of the class parameters 111 133 112 // prepare all tables << 134 G4PhysicsTable* BuildDEDXTable(); 113 void PreparePhysicsTable(const G4ParticleDef << 135 G4PhysicsTable* BuildDEDXTableForPreciseRange(); >> 136 G4PhysicsTable* BuildLambdaTable(); >> 137 G4PhysicsTable* BuildLambdaSubTable(); >> 138 // Build tables 114 139 115 // build all tables << 140 void SetBaseParticle(const G4ParticleDefinition* p); 116 void BuildPhysicsTable(const G4ParticleDefin << 117 141 118 // build a table << 142 const G4ParticleDefinition* Particle() const; 119 G4PhysicsTable* BuildDEDXTable(G4EmTableType << 143 const G4ParticleDefinition* BaseParticle() const; >> 144 const G4ParticleDefinition* SecondaryParticle() const; >> 145 // Particle definition 120 146 121 // build a table << 147 void SetDEDXBinning(G4int nbins); 122 G4PhysicsTable* BuildLambdaTable(G4EmTableTy << 148 // Binning for dEdx, range, and inverse range tables 123 149 124 // Called before tracking of each new G4Trac << 150 void SetDEDXBinningForPreciseRange(G4int nbins); 125 void StartTracking(G4Track*) override; << 151 // Binning for dEdx, range, and inverse range tables 126 152 127 // Step limit from AlongStep << 153 void SetLambdaBinning(G4int nbins); 128 G4double AlongStepGetPhysicalInteractionLeng << 154 // Binning for lambda table 129 const G4Trac << 130 G4double pr << 131 G4double cu << 132 G4double& cu << 133 G4GPILSelect << 134 155 135 // Step limit from cross section << 156 void SetMinKinEnergy(G4double e); 136 G4double PostStepGetPhysicalInteractionLengt << 157 G4double MinKinEnergy() const; 137 const G4Trac << 158 // Min kinetic energy for tables 138 G4double pre << 139 G4ForceCondi << 140 159 141 // AlongStep computations << 160 void SetMaxKinEnergy(G4double e); 142 G4VParticleChange* AlongStepDoIt(const G4Tra << 161 G4double MaxKinEnergy() const; >> 162 // Max kinetic energy for tables 143 163 144 // PostStep sampling of secondaries << 164 void SetMaxKinEnergyForPreciseRange(G4double e); 145 G4VParticleChange* PostStepDoIt(const G4Trac << 165 // Max kinetic energy for tables 146 166 147 // Store all PhysicsTable in files. << 148 // Return false in case of any fatal failure << 149 G4bool StorePhysicsTable(const G4ParticleDef 167 G4bool StorePhysicsTable(const G4ParticleDefinition*, 150 const G4String& dir 168 const G4String& directory, 151 G4bool ascii = fals << 169 G4bool ascii = false); >> 170 // Store PhysicsTable in a file. >> 171 // Return false in case of failure at I/O 152 172 153 // Retrieve all Physics from a files. << 154 // Return true if all the Physics Table are << 155 // Return false if any fatal failure. << 156 G4bool RetrievePhysicsTable(const G4Particle 173 G4bool RetrievePhysicsTable(const G4ParticleDefinition*, 157 const G4String& 174 const G4String& directory, 158 G4bool ascii) ov << 175 G4bool ascii); >> 176 // Retrieve Physics from a file. >> 177 // (return true if the Physics Table can be build by using file) >> 178 // (return false if the process has no functionality or in case of failure) >> 179 // File name should is constructed as processName+particleName and the >> 180 // should be placed under the directory specifed by the argument. 159 181 160 private: << 182 void AddEmModel(G4int, G4VEmModel*, G4VEmFluctuationModel* fluc = 0, >> 183 const G4Region* region = 0); >> 184 // Add EM model coupled with fluctuation model for the region 161 185 162 // summary printout after initialisation << 186 void UpdateEmModel(const G4String&, G4double, G4double); 163 void StreamInfo(std::ostream& out, const G4P << 187 // Define new energy range for thhe model identified by the name 164 G4bool rst=false) const; << 165 << 166 //------------------------------------------ << 167 // Public interface to cross section, mfp an << 168 // These methods are not used in run time << 169 //------------------------------------------ << 170 188 171 public: << 189 void AddSubCutoffProcessor(G4VSubCutoffProcessor*, const G4Region* region = 0); 172 << 190 // Add subcutoff processor for the region 173 // access to dispersion of restricted energy << 174 G4double GetDEDXDispersion(const G4MaterialC << 175 const G4DynamicPa << 176 G4double length); << 177 191 178 // Access to cross section table << 192 virtual void ActivateFluorescence(G4bool, const G4Region* region = 0); 179 G4double CrossSectionPerVolume(G4double kine << 193 virtual void ActivateAugerElectronProduction(G4bool, const G4Region* region = 0); 180 const G4Mater << 194 // Activate deexcitation code 181 G4double CrossSectionPerVolume(G4double kine << 182 const G4Mater << 183 G4double logK << 184 195 185 // access to cross section << 196 virtual void SetSubCutoff(G4bool); 186 G4double MeanFreePath(const G4Track& track); << 187 197 188 // access to step limit << 198 void SetDEDXTable(G4PhysicsTable* p); 189 G4double ContinuousStepLimit(const G4Track& << 199 G4PhysicsTable* DEDXTable() const; 190 G4double previo << 191 G4double curren << 192 G4double& curre << 193 << 194 protected: << 195 << 196 // implementation of the pure virtual method << 197 G4double GetMeanFreePath(const G4Track& trac << 198 G4double previousSt << 199 G4ForceCondition* c << 200 << 201 // implementation of the pure virtual method << 202 G4double GetContinuousStepLimit(const G4Trac << 203 G4double pre << 204 G4double cur << 205 G4double& cu << 206 << 207 // creation of an empty vector for cross sec << 208 G4PhysicsVector* LambdaPhysicsVector(const G << 209 G4doubl << 210 << 211 inline std::size_t CurrentMaterialCutsCouple << 212 << 213 //------------------------------------------ << 214 // Specific methods to set, access, modify m << 215 //------------------------------------------ << 216 << 217 // Select model in run time << 218 inline void SelectModel(G4double kinEnergy); << 219 << 220 public: << 221 // Select model by energy and couple index << 222 // Not for run time processing << 223 inline G4VEmModel* SelectModelForMaterial(G4 << 224 st << 225 << 226 // Add EM model coupled with fluctuation mod << 227 // of order defines which pair of models wil << 228 // energy interval << 229 void AddEmModel(G4int, G4VEmModel*, << 230 G4VEmFluctuationModel* fluc << 231 const G4Region* region = nul << 232 << 233 // Assign a model to a process local list, t << 234 // the derived process should execute AddEmM << 235 void SetEmModel(G4VEmModel*, G4int index=0); << 236 << 237 // Access to models << 238 inline std::size_t NumberOfModels() const; << 239 << 240 // Return a model from the local list << 241 inline G4VEmModel* EmModel(std::size_t index << 242 << 243 // Access to models from G4EmModelManager li << 244 inline G4VEmModel* GetModelByIndex(std::size << 245 200 246 // Assign a fluctuation model to a process << 201 void SetDEDXunRestrictedTable(G4PhysicsTable* p); 247 inline void SetFluctModel(G4VEmFluctuationMo << 202 G4PhysicsTable* DEDXunRestrictedTable() const; 248 << 249 // Return the assigned fluctuation model << 250 inline G4VEmFluctuationModel* FluctModel() c << 251 << 252 //------------------------------------------ << 253 // Define and access particle type << 254 //------------------------------------------ << 255 203 256 protected: << 204 void SetPreciseRangeTable(G4PhysicsTable* pRange); 257 inline void SetParticle(const G4ParticleDefi << 205 G4PhysicsTable* PreciseRangeTable() const; 258 inline void SetSecondaryParticle(const G4Par << 259 206 260 public: << 207 void SetRangeTableForLoss(G4PhysicsTable* p); 261 inline void SetBaseParticle(const G4Particle << 208 G4PhysicsTable* RangeTableForLoss() const; 262 inline const G4ParticleDefinition* Particle( << 263 inline const G4ParticleDefinition* BaseParti << 264 inline const G4ParticleDefinition* Secondary << 265 209 266 // hide assignment operator << 210 void SetInverseRangeTable(G4PhysicsTable* p); 267 G4VEnergyLossProcess(G4VEnergyLossProcess &) << 211 G4PhysicsTable* InverseRangeTable() const; 268 G4VEnergyLossProcess & operator=(const G4VEn << 269 212 270 //------------------------------------------ << 213 void SetSecondaryRangeTable(G4PhysicsTable* p); 271 // Get/set parameters to configure the proce << 272 //------------------------------------------ << 273 214 274 // Add subcut processor for the region << 215 void SetLambdaTable(G4PhysicsTable* p); 275 void ActivateSubCutoff(const G4Region* regio << 216 G4PhysicsTable* LambdaTable(); 276 217 277 // Activate biasing << 218 void SetSubLambdaTable(G4PhysicsTable* p); 278 void SetCrossSectionBiasingFactor(G4double f << 219 G4PhysicsTable* SubLambdaTable(); 279 220 280 void ActivateForcedInteraction(G4double leng << 221 // Return values for particle 281 const G4Strin << 222 G4double GetDEDX(G4double& kineticEnergy, const G4MaterialCutsCouple* couple); 282 G4bool flag = << 223 G4double GetRange(G4double& kineticEnergy, const G4MaterialCutsCouple* couple); >> 224 G4double GetRangeForLoss(G4double& kineticEnergy, const G4MaterialCutsCouple* couple); >> 225 G4double GetKineticEnergy(G4double& range, const G4MaterialCutsCouple* couple); >> 226 G4double GetLambda(G4double& kineticEnergy, const G4MaterialCutsCouple* couple); 283 227 284 void ActivateSecondaryBiasing(const G4String << 228 G4double GetDEDXDispersion(const G4MaterialCutsCouple *couple, 285 G4double energ << 229 const G4DynamicParticle* dp, >> 230 G4double length); 286 231 287 inline void SetLossFluctuations(G4bool val); << 232 G4double MicroscopicCrossSection(G4double kineticEnergy, >> 233 const G4MaterialCutsCouple* couple); 288 234 289 inline void SetSpline(G4bool val); << 235 void SetLossFluctuations(G4bool val); 290 inline void SetCrossSectionType(G4CrossSecti << 291 inline G4CrossSectionType CrossSectionType() << 292 236 293 // Set/Get flag "isIonisation" << 237 void SetIntegral(G4bool val); 294 void SetIonisation(G4bool val); << 238 G4bool IsIntegral() const; 295 inline G4bool IsIonisationProcess() const; << 296 239 297 // Redefine parameteters for stepping contro 240 // Redefine parameteters for stepping control >> 241 // 298 void SetLinearLossLimit(G4double val); 242 void SetLinearLossLimit(G4double val); >> 243 void SetMinSubRange(G4double val); >> 244 void SetStepLimits(G4double v1, G4double v2); 299 void SetStepFunction(G4double v1, G4double v 245 void SetStepFunction(G4double v1, G4double v2); 300 void SetLowestEnergyLimit(G4double); << 246 void SetLambdaFactor(G4double val); 301 247 302 inline G4int NumberOfSubCutoffRegions() cons << 248 G4bool TablesAreBuilt() const; 303 249 304 //------------------------------------------ << 250 G4int NumberOfSubCutoffRegions() const; 305 // Specific methods to path Physics Tables t << 306 //------------------------------------------ << 307 251 308 void SetDEDXTable(G4PhysicsTable* p, G4EmTab << 252 // Helper functions 309 void SetCSDARangeTable(G4PhysicsTable* pRang << 253 G4double MeanFreePath(const G4Track& track, 310 void SetRangeTableForLoss(G4PhysicsTable* p) << 254 G4double previousStepSize, 311 void SetInverseRangeTable(G4PhysicsTable* p) << 255 G4ForceCondition* condition); 312 void SetLambdaTable(G4PhysicsTable* p); << 313 256 314 void SetTwoPeaksXS(std::vector<G4TwoPeaksXS* << 257 G4double ContinuousStepLimit(const G4Track& track, 315 void SetEnergyOfCrossSectionMax(std::vector< << 258 G4double previousStepSize, >> 259 G4double currentMinimumStep, >> 260 G4double& currentSafety); 316 261 317 //------------------------------------------ << 262 void ResetNumberOfInteractionLengthLeft(); 318 // Specific methods to define custom Physics << 263 // reset (determine the value of)NumberOfInteractionLengthLeft 319 //------------------------------------------ << 320 264 321 // Binning for dEdx, range, inverse range an << 265 G4VEmModel* SelectModelForMaterial(G4double kinEnergy, size_t& idx) const; 322 void SetDEDXBinning(G4int nbins); << 323 266 324 // Min kinetic energy for tables << 267 void SetIonisation(G4bool val); 325 void SetMinKinEnergy(G4double e); << 326 inline G4double MinKinEnergy() const; << 327 268 328 // Max kinetic energy for tables << 269 G4bool IsIonisationProcess() const; 329 void SetMaxKinEnergy(G4double e); << 330 inline G4double MaxKinEnergy() const; << 331 270 332 // Biasing parameters << 271 protected: 333 inline G4double CrossSectionBiasingFactor() << 334 272 335 // Return values for given G4MaterialCutsCou << 273 virtual G4double MinPrimaryEnergy(const G4ParticleDefinition*, 336 inline G4double GetDEDX(G4double kineticEner << 274 const G4Material*, G4double cut) = 0; 337 inline G4double GetCSDADEDX(G4double kinetic << 338 const G4Material << 339 inline G4double GetDEDX(G4double kineticEner << 340 G4double logKineticE << 341 inline G4double GetRange(G4double kineticEne << 342 inline G4double GetRange(G4double kineticEne << 343 G4double logKinetic << 344 inline G4double GetCSDARange(G4double kineti << 345 const G4Materia << 346 inline G4double GetKineticEnergy(G4double ra << 347 const G4Mat << 348 inline G4double GetLambda(G4double kineticEn << 349 inline G4double GetLambda(G4double kineticEn << 350 G4double logKineti << 351 << 352 inline G4bool TablesAreBuilt() const; << 353 << 354 // Access to specific tables << 355 inline G4PhysicsTable* DEDXTable() const; << 356 inline G4PhysicsTable* DEDXunRestrictedTable << 357 inline G4PhysicsTable* IonisationTable() con << 358 inline G4PhysicsTable* CSDARangeTable() cons << 359 inline G4PhysicsTable* RangeTableForLoss() c << 360 inline G4PhysicsTable* InverseRangeTable() c << 361 inline G4PhysicsTable* LambdaTable() const; << 362 inline std::vector<G4TwoPeaksXS*>* TwoPeaksX << 363 inline std::vector<G4double>* EnergyOfCrossS << 364 << 365 inline G4bool UseBaseMaterial() const; << 366 << 367 //------------------------------------------ << 368 // Run time method for simulation of ionisat << 369 //------------------------------------------ << 370 275 371 // access atom on which interaction happens << 276 virtual G4double MaxSecondaryEnergy(const G4DynamicParticle* dp) = 0; 372 const G4Element* GetCurrentElement() const; << 373 277 374 // Set scaling parameters for ions is needed << 278 virtual void InitialiseEnergyLossProcess(const G4ParticleDefinition*, 375 void SetDynamicMassCharge(G4double massratio << 279 const G4ParticleDefinition*) = 0; 376 280 377 private: << 281 void SetParticle(const G4ParticleDefinition* p); 378 282 379 void FillSecondariesAlongStep(G4double weigh << 283 void SetSecondaryParticle(const G4ParticleDefinition* p); 380 284 381 void PrintWarning(const G4String&, G4double << 285 virtual G4double GetMeanFreePath(const G4Track& track, >> 286 G4double previousStepSize, >> 287 G4ForceCondition* condition); 382 288 383 // define material and indexes << 289 virtual G4double GetContinuousStepLimit(const G4Track& track, 384 inline void DefineMaterial(const G4MaterialC << 290 G4double previousStepSize, >> 291 G4double currentMinimumStep, >> 292 G4double& currentSafety); 385 293 386 //------------------------------------------ << 294 virtual 387 // Compute values using scaling relation, ma << 295 G4PhysicsVector* DEDXPhysicsVector(const G4MaterialCutsCouple*); 388 //------------------------------------------ << 389 inline G4double GetDEDXForScaledEnergy(G4dou << 390 inline G4double GetDEDXForScaledEnergy(G4dou << 391 G4dou << 392 inline G4double GetIonisationForScaledEnergy << 393 inline G4double GetScaledRangeForScaledEnerg << 394 inline G4double GetScaledRangeForScaledEnerg << 395 << 396 << 397 inline G4double GetLimitScaledRangeForScaled << 398 inline G4double GetLimitScaledRangeForScaled << 399 << 400 << 401 inline G4double ScaledKinEnergyForLoss(G4dou << 402 inline G4double GetLambdaForScaledEnergy(G4d << 403 inline G4double GetLambdaForScaledEnergy(G4d << 404 G4d << 405 296 406 inline G4double LogScaledEkin(const G4Track& << 297 virtual 407 << 298 G4PhysicsVector* DEDXPhysicsVectorForPreciseRange(const G4MaterialCutsCouple*); 408 void ComputeLambdaForScaledEnergy(G4double s << 409 const G4Tr << 410 299 411 G4bool IsRegionForCubcutProcessor(const G4Tr << 300 virtual >> 301 G4PhysicsVector* LambdaPhysicsVector(const G4MaterialCutsCouple*); 412 302 413 protected: << 303 virtual >> 304 G4PhysicsVector* SubLambdaPhysicsVector(const G4MaterialCutsCouple*); 414 305 415 G4ParticleChangeForLoss fParticleChange; << 306 G4VEmModel* SelectModel(G4double kinEnergy); 416 const G4Material* currentMaterial << 417 const G4MaterialCutsCouple* currentCouple = << 418 307 419 private: << 308 G4VSubCutoffProcessor* SubCutoffProcessor(size_t index); 420 309 421 G4LossTableManager* lManager; << 310 size_t CurrentMaterialCutsCoupleIndex() const; 422 G4EmModelManager* modelManager; << 423 G4VEmModel* currentModel = n << 424 G4EmBiasingManager* biasManager = nu << 425 G4SafetyHelper* safetyHelper; << 426 G4EmParameters* theParameters; << 427 G4VEmFluctuationModel* fluctModel = nul << 428 G4VAtomDeexcitation* atomDeexcitation << 429 G4VSubCutProducer* subcutProducer = << 430 << 431 const G4ParticleDefinition* particle = nullp << 432 const G4ParticleDefinition* baseParticle = n << 433 const G4ParticleDefinition* secondaryParticl << 434 G4EmDataHandler* theData = nullptr; << 435 << 436 G4PhysicsTable* theDEDXTable = nullptr; << 437 G4PhysicsTable* theDEDXunRestrictedTable = n << 438 G4PhysicsTable* theIonisationTable = nullptr << 439 G4PhysicsTable* theRangeTableForLoss = nullp << 440 G4PhysicsTable* theCSDARangeTable = nullptr; << 441 G4PhysicsTable* theInverseRangeTable = nullp << 442 G4PhysicsTable* theLambdaTable = nullptr; << 443 << 444 std::vector<const G4Region*>* scoffRegions = << 445 std::vector<G4VEmModel*>* emModels = nul << 446 const std::vector<G4int>* theDensityIdx << 447 const std::vector<G4double>* theDensityFact << 448 const G4DataVector* theCuts = null << 449 311 450 std::vector<G4double>* theEnergyOfCrossSecti << 312 // Set scaling parameters 451 std::vector<G4TwoPeaksXS*>* fXSpeaks = nullp << 313 // >> 314 void SetMassRatio(G4double val); >> 315 void SetReduceFactor(G4double val); >> 316 void SetChargeSquare(G4double val); >> 317 void SetChargeSquareRatio(G4double val); 452 318 453 G4double lowestKinEnergy; << 319 G4double GetCurrentRange() const; 454 G4double minKinEnergy; << 455 G4double maxKinEnergy; << 456 G4double maxKinEnergyCSDA; << 457 320 458 G4double linLossLimit = 0.01; << 321 private: 459 G4double dRoverRange = 0.2; << 460 G4double finalRange; << 461 G4double lambdaFactor = 0.8; << 462 G4double invLambdaFactor; << 463 G4double biasFactor = 1.0; << 464 << 465 G4double massRatio = 1.0; << 466 G4double logMassRatio = 0.0; << 467 G4double fFactor = 1.0; << 468 G4double reduceFactor = 1.0; << 469 G4double chargeSqRatio = 1.0; << 470 G4double fRange = 0.0; << 471 G4double fRangeEnergy = 0.0; << 472 322 473 protected: << 323 // Clear tables >> 324 void Clear(); 474 325 475 G4double preStepLambda = 0.0; << 326 void DefineMaterial(const G4MaterialCutsCouple* couple); 476 G4double preStepKinEnergy = 0.0; << 477 G4double preStepScaledEnergy = 0.0; << 478 G4double mfpKinEnergy = 0.0; << 479 327 480 std::size_t currentCoupleIndex = 0; << 328 // Returnd values for scaled energy and base particles mass >> 329 // >> 330 G4double GetDEDXForScaledEnergy(G4double scaledKinEnergy); >> 331 G4double GetScaledRangeForScaledEnergy(G4double scaledKinEnergy); >> 332 G4double GetLimitScaledRangeForScaledEnergy(G4double scaledKinEnergy); >> 333 G4double GetLambdaForScaledEnergy(G4double scaledKinEnergy); >> 334 G4double ScaledKinEnergyForLoss(G4double range); >> 335 void ComputeLambdaForScaledEnergy(G4double scaledKinEnergy); 481 336 482 private: << 337 // hide assignment operator 483 338 484 G4int nBins; << 339 G4VEnergyLossProcess(G4VEnergyLossProcess &); 485 G4int nBinsCSDA; << 340 G4VEnergyLossProcess & operator=(const G4VEnergyLossProcess &right); 486 G4int numberOfModels = 0; << 487 G4int nSCoffRegions = 0; << 488 G4int secID = _DeltaElectron; << 489 G4int tripletID = _TripletElectron; << 490 G4int biasID = _DeltaEBelowCut; << 491 G4int epixeID = _ePIXE; << 492 G4int gpixeID = _GammaPIXE; << 493 G4int mainSecondaries = 1; << 494 << 495 std::size_t basedCoupleIndex = 0; << 496 std::size_t coupleIdxRange = 0; << 497 std::size_t idxDEDX = 0; << 498 std::size_t idxDEDXunRestricted = 0; << 499 std::size_t idxIonisation = 0; << 500 std::size_t idxRange = 0; << 501 std::size_t idxCSDA = 0; << 502 std::size_t idxSecRange = 0; << 503 std::size_t idxInverseRange = 0; << 504 std::size_t idxLambda = 0; << 505 << 506 G4GPILSelection aGPILSelection; << 507 G4CrossSectionType fXSType = fEmOnePeak; << 508 << 509 G4bool lossFluctuationFlag = true; << 510 G4bool useCutAsFinalRange = false; << 511 G4bool tablesAreBuilt = false; << 512 G4bool spline = true; << 513 G4bool isIon = false; << 514 G4bool isIonisation = false; << 515 G4bool useDeexcitation = false; << 516 G4bool biasFlag = false; << 517 G4bool weightFlag = false; << 518 G4bool isMaster = false; << 519 G4bool baseMat = false; << 520 G4bool actLinLossLimit = false; << 521 G4bool actLossFluc = false; << 522 G4bool actBinning = false; << 523 G4bool actMinKinEnergy = false; << 524 G4bool actMaxKinEnergy = false; << 525 341 526 std::vector<G4DynamicParticle*> secParticles << 342 // ===================================================================== 527 std::vector<G4Track*> scTracks; << 528 }; << 529 343 530 // ======== Run time inline methods ========== << 344 protected: 531 345 532 inline std::size_t G4VEnergyLossProcess::Curre << 346 G4ParticleChangeForLoss fParticleChange; 533 { << 534 return currentCoupleIndex; << 535 } << 536 347 537 //....oooOO0OOooo........oooOO0OOooo........oo << 348 private: 538 349 539 inline void G4VEnergyLossProcess::SelectModel( << 350 G4EmModelManager* modelManager; 540 { << 351 std::vector<G4VSubCutoffProcessor*> scoffProcessors; 541 currentModel = modelManager->SelectModel(kin << 352 std::vector<const G4Region*> scoffRegions; 542 currentModel->SetCurrentCouple(currentCouple << 353 G4int nSCoffRegions; 543 } << 354 std::vector<G4int> idxSCoffRegions; >> 355 >> 356 // tables and vectors >> 357 G4PhysicsTable* theDEDXTable; >> 358 G4PhysicsTable* theRangeTableForLoss; >> 359 G4PhysicsTable* theDEDXunRestrictedTable; >> 360 G4PhysicsTable* thePreciseRangeTable; >> 361 G4PhysicsTable* theSecondaryRangeTable; >> 362 G4PhysicsTable* theInverseRangeTable; >> 363 G4PhysicsTable* theLambdaTable; >> 364 G4PhysicsTable* theSubLambdaTable; >> 365 G4double* theDEDXAtMaxEnergy; >> 366 G4double* theRangeAtMaxEnergy; >> 367 G4double* theEnergyOfCrossSectionMax; >> 368 G4double* theCrossSectionMax; >> 369 >> 370 const G4DataVector* theCuts; >> 371 >> 372 const G4ParticleDefinition* particle; >> 373 const G4ParticleDefinition* baseParticle; >> 374 const G4ParticleDefinition* secondaryParticle; >> 375 >> 376 // cash >> 377 const G4Material* currentMaterial; >> 378 const G4MaterialCutsCouple* currentCouple; >> 379 size_t currentMaterialIndex; >> 380 G4double minStepLimit; >> 381 >> 382 G4int nDEDXBins; >> 383 G4int nDEDXBinsForRange; >> 384 G4int nLambdaBins; 544 385 545 //....oooOO0OOooo........oooOO0OOooo........oo << 386 G4double lowestKinEnergy; >> 387 G4double minKinEnergy; >> 388 G4double maxKinEnergy; >> 389 G4double maxKinEnergyForRange; 546 390 547 inline G4VEmModel* G4VEnergyLossProcess::Selec << 391 G4double massRatio; 548 G4double kinEnergy, std::si << 392 G4double reduceFactor; 549 { << 393 G4double chargeSquare; 550 return modelManager->SelectModel(kinEnergy, << 394 G4double chargeSqRatio; 551 } << 395 >> 396 G4double preStepLambda; >> 397 G4double preStepMFP; >> 398 G4double fRange; >> 399 G4double preStepKinEnergy; >> 400 G4double preStepScaledEnergy; >> 401 G4double linLossLimit; >> 402 G4double minSubRange; >> 403 G4double dRoverRange; >> 404 G4double finalRange; >> 405 G4double defaultRoverRange; >> 406 G4double defaultIntegralRange; >> 407 G4double lambdaFactor; >> 408 G4double mfpKinEnergy; >> 409 >> 410 G4bool lossFluctuationFlag; >> 411 G4bool rndmStepFlag; >> 412 G4bool tablesAreBuilt; >> 413 G4bool integral; >> 414 G4bool meanFreePath; >> 415 G4bool aboveCSmax; >> 416 G4bool isIonisation; >> 417 }; 552 418 553 //....oooOO0OOooo........oooOO0OOooo........oo 419 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 554 421 555 inline void << 422 inline void G4VEnergyLossProcess::DefineMaterial(const G4MaterialCutsCouple* couple) 556 G4VEnergyLossProcess::DefineMaterial(const G4M << 557 { 423 { 558 if(couple != currentCouple) { 424 if(couple != currentCouple) { 559 currentCouple = couple; << 425 currentCouple = couple; 560 currentMaterial = couple->GetMaterial(); 426 currentMaterial = couple->GetMaterial(); 561 basedCoupleIndex = currentCoupleIndex = co << 427 currentMaterialIndex = couple->GetIndex(); 562 fFactor = chargeSqRatio*biasFactor; << 428 minStepLimit = std::min(finalRange, 563 mfpKinEnergy = DBL_MAX; << 429 currentCouple->GetProductionCuts()->GetProductionCut(idxG4ElectronCut)); 564 idxLambda = 0; << 430 if(!meanFreePath) ResetNumberOfInteractionLengthLeft(); 565 if(baseMat) { << 566 basedCoupleIndex = (*theDensityIdx)[curr << 567 fFactor *= (*theDensityFactor)[currentCo << 568 } << 569 reduceFactor = 1.0/(fFactor*massRatio); << 570 } 431 } 571 } 432 } 572 433 573 //....oooOO0OOooo........oooOO0OOooo........oo 434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 574 435 575 inline G4double G4VEnergyLossProcess::GetDEDXF << 436 inline G4double G4VEnergyLossProcess::GetDEDX(G4double& kineticEnergy, 576 { << 437 const G4MaterialCutsCouple* couple) 577 /* << 578 G4cout << "G4VEnergyLossProcess::GetDEDX: Id << 579 << basedCoupleIndex << " E(MeV)= " << 580 << " Emin= " << minKinEnergy << " Fa << 581 << " " << theDEDXTable << G4endl; */ << 582 G4double x = fFactor*(*theDEDXTable)[basedCo << 583 if(e < minKinEnergy) { x *= std::sqrt(e/minK << 584 return x; << 585 } << 586 << 587 //....oooOO0OOooo........oooOO0OOooo........oo << 588 << 589 inline << 590 G4double G4VEnergyLossProcess::GetDEDXForScale << 591 { 438 { 592 /* << 439 DefineMaterial(couple); 593 G4cout << "G4VEnergyLossProcess::GetDEDX: Id << 440 return GetDEDXForScaledEnergy(kineticEnergy*massRatio); 594 << basedCoupleIndex << " E(MeV)= " << 595 << " Emin= " << minKinEnergy << " Fa << 596 << " " << theDEDXTable << G4endl; */ << 597 G4double x = fFactor*(*theDEDXTable)[basedCo << 598 if(e < minKinEnergy) { x *= std::sqrt(e/minK << 599 return x; << 600 } 441 } 601 442 602 //....oooOO0OOooo........oooOO0OOooo........oo 443 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 603 444 604 inline G4double G4VEnergyLossProcess::GetIonis << 445 inline G4double G4VEnergyLossProcess::GetDEDXForScaledEnergy(G4double e) 605 { 446 { 606 G4double x = << 447 G4bool b; 607 fFactor*(*theIonisationTable)[basedCoupleI << 448 G4double x = ((*theDEDXTable)[currentMaterialIndex]->GetValue(e, b))*chargeSqRatio; 608 if(e < minKinEnergy) { x *= std::sqrt(e/minK << 449 if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy); 609 return x; 450 return x; 610 } 451 } 611 452 612 //....oooOO0OOooo........oooOO0OOooo........oo 453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 613 454 614 inline G4double G4VEnergyLossProcess::GetScale << 455 inline G4double G4VEnergyLossProcess::GetRange(G4double& kineticEnergy, >> 456 const G4MaterialCutsCouple* couple) 615 { 457 { 616 //G4cout << "G4VEnergyLossProcess::GetScaled << 458 DefineMaterial(couple); 617 // << basedCoupleIndex << " E(MeV)= << 459 G4double x = DBL_MAX; 618 // << " lastIdx= " << lastIdx << " << 460 if(thePreciseRangeTable) 619 if(currentCoupleIndex != coupleIdxRange || f << 461 x = GetLimitScaledRangeForScaledEnergy(kineticEnergy*massRatio); 620 coupleIdxRange = currentCoupleIndex; << 462 else if(theRangeTableForLoss) 621 fRangeEnergy = e; << 463 x = GetScaledRangeForScaledEnergy(kineticEnergy*massRatio); 622 fRange = reduceFactor*((*theRangeTableForL << 464 return x*reduceFactor; 623 if (fRange < 0.0) { fRange = 0.0; } << 624 else if (e < minKinEnergy) { fRange *= std << 625 } << 626 //G4cout << "G4VEnergyLossProcess::GetScaled << 627 // << basedCoupleIndex << " E(MeV)= << 628 // << " R= " << computedRange << " << 629 return fRange; << 630 } << 631 << 632 inline G4double << 633 G4VEnergyLossProcess::GetScaledRangeForScaledE << 634 { << 635 //G4cout << "G4VEnergyLossProcess::GetScaled << 636 // << basedCoupleIndex << " E(MeV)= << 637 // << " lastIdx= " << lastIdx << " << 638 if(currentCoupleIndex != coupleIdxRange || f << 639 coupleIdxRange = currentCoupleIndex; << 640 fRangeEnergy = e; << 641 fRange = reduceFactor*((*theRangeTableForL << 642 if (fRange < 0.0) { fRange = 0.0; } << 643 else if (e < minKinEnergy) { fRange *= std << 644 } << 645 //G4cout << "G4VEnergyLossProcess::GetScaled << 646 // << basedCoupleIndex << " E(MeV)= << 647 // << " R= " << fRange << " " << t << 648 return fRange; << 649 } 465 } 650 466 651 //....oooOO0OOooo........oooOO0OOooo........oo 467 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 652 468 653 inline G4double << 469 inline G4double G4VEnergyLossProcess::GetLimitScaledRangeForScaledEnergy(G4double e) 654 G4VEnergyLossProcess::GetLimitScaledRangeForSc << 655 { 470 { 656 G4double x = ((*theCSDARangeTable)[basedCoup << 471 G4bool b; 657 if (x < 0.0) { x = 0.0; } << 472 G4double x; 658 else if (e < minKinEnergy) { x *= std::sqrt( << 659 return x; << 660 } << 661 473 662 //....oooOO0OOooo........oooOO0OOooo........oo << 474 if (e < maxKinEnergyForRange) { >> 475 x = ((*thePreciseRangeTable)[currentMaterialIndex])->GetValue(e, b); >> 476 if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy); 663 477 664 inline G4double << 478 } else { 665 G4VEnergyLossProcess::GetLimitScaledRangeForSc << 479 x = theRangeAtMaxEnergy[currentMaterialIndex] + 666 << 480 (e - maxKinEnergyForRange)/theDEDXAtMaxEnergy[currentMaterialIndex]; 667 { << 668 G4double x = ((*theCSDARangeTable)[basedCoup << 669 if (x < 0.0) { x = 0.0; } << 670 else if (e < minKinEnergy) { x *= std::sqrt( << 671 return x; << 672 } << 673 << 674 //....oooOO0OOooo........oooOO0OOooo........oo << 675 << 676 inline G4double G4VEnergyLossProcess::ScaledKi << 677 { << 678 //G4cout << "G4VEnergyLossProcess::GetEnergy << 679 // << basedCoupleIndex << " R(mm)= " << 680 // << theInverseRangeTable << G4endl << 681 G4PhysicsVector* v = (*theInverseRangeTable) << 682 G4double rmin = v->Energy(0); << 683 G4double e = 0.0; << 684 if(r >= rmin) { e = v->Value(r, idxInverseRa << 685 else if(r > 0.0) { << 686 G4double x = r/rmin; << 687 e = minKinEnergy*x*x; << 688 } 481 } 689 return e; << 482 return x; 690 } << 691 << 692 //....oooOO0OOooo........oooOO0OOooo........oo << 693 << 694 inline G4double G4VEnergyLossProcess::GetLambd << 695 { << 696 return fFactor*((*theLambdaTable)[basedCoupl << 697 } 483 } 698 484 699 //....oooOO0OOooo........oooOO0OOooo........oo 485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 700 486 701 inline G4double << 487 inline G4double G4VEnergyLossProcess::GetRangeForLoss(G4double& kineticEnergy, 702 G4VEnergyLossProcess::GetLambdaForScaledEnergy << 488 const G4MaterialCutsCouple* couple) 703 { 489 { 704 return fFactor*((*theLambdaTable)[basedCoupl << 490 DefineMaterial(couple); >> 491 G4double x = DBL_MAX; >> 492 if(theRangeTableForLoss) x = GetScaledRangeForScaledEnergy(kineticEnergy*massRatio); >> 493 return x*reduceFactor; 705 } 494 } 706 495 707 //....oooOO0OOooo........oooOO0OOooo........oo 496 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 708 497 709 inline G4double G4VEnergyLossProcess::LogScale << 498 inline G4double G4VEnergyLossProcess::GetScaledRangeForScaledEnergy(G4double e) 710 { 499 { 711 return track.GetDynamicParticle()->GetLogKin << 500 G4bool b; >> 501 G4double x = ((*theRangeTableForLoss)[currentMaterialIndex])->GetValue(e, b); >> 502 if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy); >> 503 return x; 712 } 504 } 713 505 714 //....oooOO0OOooo........oooOO0OOooo........oo 506 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 715 507 716 inline G4double << 508 inline G4double G4VEnergyLossProcess::GetKineticEnergy(G4double& range, 717 G4VEnergyLossProcess::GetDEDX(G4double kinEner << 509 const G4MaterialCutsCouple* couple) 718 const G4Material << 719 { 510 { 720 DefineMaterial(couple); 511 DefineMaterial(couple); 721 return GetDEDXForScaledEnergy(kinEnergy*mass << 512 G4double r = range/reduceFactor; >> 513 G4double e = ScaledKinEnergyForLoss(r)/massRatio; >> 514 return e; 722 } 515 } 723 516 724 //....oooOO0OOooo........oooOO0OOooo........oo 517 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 725 518 726 inline G4double << 519 inline G4double G4VEnergyLossProcess::ScaledKinEnergyForLoss(G4double r) 727 G4VEnergyLossProcess::GetDEDX(G4double kinEner << 728 const G4Material << 729 G4double logKinE << 730 { 520 { 731 DefineMaterial(couple); << 521 G4PhysicsVector* v = (*theInverseRangeTable)[currentMaterialIndex]; 732 return GetDEDXForScaledEnergy(kinEnergy*mass << 522 G4double rmin = v->GetLowEdgeEnergy(0); >> 523 G4double e = minKinEnergy; >> 524 if(r <= rmin) { >> 525 r /= rmin; >> 526 e *= r*r; >> 527 } else { >> 528 G4bool b; >> 529 e = v->GetValue(r, b); >> 530 } >> 531 return e; 733 } 532 } 734 533 735 //....oooOO0OOooo........oooOO0OOooo........oo 534 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 736 535 737 inline G4double << 536 inline G4double G4VEnergyLossProcess::GetDEDXDispersion( 738 G4VEnergyLossProcess::GetRange(G4double kinEne << 537 const G4MaterialCutsCouple *couple, 739 const G4Materia << 538 const G4DynamicParticle* dp, >> 539 G4double length) 740 { 540 { 741 DefineMaterial(couple); 541 DefineMaterial(couple); 742 return GetScaledRangeForScaledEnergy(kinEner << 542 G4double tmax = MaxSecondaryEnergy(dp); >> 543 tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]); >> 544 return modelManager->GetDEDXDispersion(currentMaterial, dp, tmax, length, >> 545 currentMaterialIndex); 743 } 546 } 744 547 745 //....oooOO0OOooo........oooOO0OOooo........oo 548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 746 549 747 inline G4double << 550 inline G4double G4VEnergyLossProcess::GetLambda(G4double& kineticEnergy, 748 G4VEnergyLossProcess::GetRange(G4double kinEne << 551 const G4MaterialCutsCouple* couple) 749 const G4Materia << 750 G4double logKin << 751 { 552 { 752 DefineMaterial(couple); 553 DefineMaterial(couple); 753 return GetScaledRangeForScaledEnergy(kinEner << 554 G4double x = 0.0; >> 555 if(theLambdaTable) x = GetLambdaForScaledEnergy(kineticEnergy*massRatio); >> 556 return x; 754 } 557 } 755 558 756 //....oooOO0OOooo........oooOO0OOooo........oo 559 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 757 560 758 inline G4double << 561 inline G4double G4VEnergyLossProcess::GetLambdaForScaledEnergy(G4double e) 759 G4VEnergyLossProcess::GetCSDARange(G4double ki << 760 const G4Mat << 761 { 562 { 762 DefineMaterial(couple); << 563 G4bool b; 763 return (nullptr == theCSDARangeTable) ? DBL_ << 564 return chargeSqRatio*(((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b)); 764 GetLimitScaledRangeForScaledEnergy(kinetic << 765 } 565 } 766 566 767 //....oooOO0OOooo........oooOO0OOooo........oo 567 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 768 568 769 inline G4double << 569 inline void G4VEnergyLossProcess::ComputeLambdaForScaledEnergy(G4double e) 770 G4VEnergyLossProcess::GetKineticEnergy(G4doubl << 771 const G << 772 { 570 { 773 DefineMaterial(couple); << 571 meanFreePath = false; 774 return ScaledKinEnergyForLoss(range/reduceFa << 572 aboveCSmax = false; >> 573 mfpKinEnergy = theEnergyOfCrossSectionMax[currentMaterialIndex]; >> 574 if (e <= mfpKinEnergy) { >> 575 preStepLambda = GetLambdaForScaledEnergy(e); >> 576 } else { >> 577 aboveCSmax = true; >> 578 G4double e1 = e*lambdaFactor; >> 579 if(e1 > mfpKinEnergy) { >> 580 preStepLambda = GetLambdaForScaledEnergy(e); >> 581 G4double preStepLambda1 = GetLambdaForScaledEnergy(e1); >> 582 if(preStepLambda1 > preStepLambda) { >> 583 mfpKinEnergy = e1; >> 584 preStepLambda = preStepLambda1; >> 585 } >> 586 } else { >> 587 preStepLambda = chargeSqRatio*theCrossSectionMax[currentMaterialIndex]; >> 588 } >> 589 } 775 } 590 } 776 591 777 //....oooOO0OOooo........oooOO0OOooo........oo 592 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 778 593 779 inline G4double << 594 inline G4double G4VEnergyLossProcess::GetMeanFreePath(const G4Track& track, 780 G4VEnergyLossProcess::GetLambda(G4double kinEn << 595 G4double, 781 const G4Materi << 596 G4ForceCondition* condition) 782 { << 597 { 783 DefineMaterial(couple); << 598 *condition = NotForced; 784 return (nullptr != theLambdaTable) ? << 599 preStepKinEnergy = track.GetKineticEnergy(); 785 GetLambdaForScaledEnergy(kinEnergy*massRat << 600 preStepScaledEnergy = preStepKinEnergy*massRatio; >> 601 if(aboveCSmax && preStepScaledEnergy < mfpKinEnergy) ResetNumberOfInteractionLengthLeft(); >> 602 DefineMaterial(track.GetMaterialCutsCouple()); >> 603 if (meanFreePath) { >> 604 if (integral) ComputeLambdaForScaledEnergy(preStepScaledEnergy); >> 605 else preStepLambda = GetLambdaForScaledEnergy(preStepScaledEnergy); >> 606 if(0.0 < preStepLambda) preStepMFP = 1.0/preStepLambda; >> 607 else preStepMFP = DBL_MAX; >> 608 } >> 609 // G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy<< " eCSmax= " <<mfpKinEnergy<< " mfp= "<<preStepMFP<<G4endl; >> 610 return preStepMFP; 786 } 611 } 787 612 788 //....oooOO0OOooo........oooOO0OOooo........oo 613 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 789 614 790 inline G4double << 615 inline G4double G4VEnergyLossProcess::GetContinuousStepLimit(const G4Track&, 791 G4VEnergyLossProcess::GetLambda(G4double kinEn << 616 G4double, G4double currentMinStep, G4double&) 792 const G4Materi << 793 G4double logKi << 794 { 617 { 795 DefineMaterial(couple); << 618 G4double x = DBL_MAX; 796 return (nullptr != theLambdaTable) ? << 619 if(theRangeTableForLoss) { 797 GetLambdaForScaledEnergy(kinEnergy*massRat << 620 fRange = GetScaledRangeForScaledEnergy(preStepScaledEnergy)*reduceFactor; 798 : 0.0; << 799 } << 800 621 801 // ======== Get/Set inline methods used at ini << 622 x = fRange; >> 623 G4double y = x*dRoverRange; 802 624 803 inline void G4VEnergyLossProcess::SetFluctMode << 625 if(x > minStepLimit && y < currentMinStep ) { 804 { << 626 x = y + minStepLimit*(1.0 - dRoverRange)*(2.0 - minStepLimit/fRange); 805 fluctModel = p; << 627 // G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy<<" range= "<<fRange >> 628 // <<" cMinSt="<<currentMinStep <<" minStepLimit= " << minStepLimit<< G4endl; >> 629 } >> 630 } >> 631 //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy<<" stepLimit= "<<x<<G4endl; >> 632 return x; 806 } 633 } 807 634 808 //....oooOO0OOooo........oooOO0OOooo........oo 635 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 809 636 810 inline G4VEmFluctuationModel* G4VEnergyLossPro << 637 inline void G4VEnergyLossProcess::ResetNumberOfInteractionLengthLeft() 811 { 638 { 812 return fluctModel; << 639 meanFreePath = true; >> 640 aboveCSmax = false; >> 641 G4VProcess::ResetNumberOfInteractionLengthLeft(); 813 } 642 } 814 643 815 //....oooOO0OOooo........oooOO0OOooo........oo << 816 << 817 inline void G4VEnergyLossProcess::SetParticle( << 818 { << 819 particle = p; << 820 } << 821 644 822 //....oooOO0OOooo........oooOO0OOooo........oo 645 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 823 646 824 inline void << 647 inline G4VEmModel* G4VEnergyLossProcess::SelectModel(G4double kinEnergy) 825 G4VEnergyLossProcess::SetSecondaryParticle(con << 826 { 648 { 827 secondaryParticle = p; << 649 return modelManager->SelectModel(kinEnergy, currentMaterialIndex); 828 } 650 } 829 651 830 //....oooOO0OOooo........oooOO0OOooo........oo 652 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 831 653 832 inline void << 654 inline G4VEmModel* G4VEnergyLossProcess::SelectModelForMaterial( 833 G4VEnergyLossProcess::SetBaseParticle(const G4 << 655 G4double kinEnergy, size_t& idx) const 834 { 656 { 835 baseParticle = p; << 657 return modelManager->SelectModel(kinEnergy, idx); 836 } 658 } 837 659 838 //....oooOO0OOooo........oooOO0OOooo........oo 660 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 839 661 840 inline const G4ParticleDefinition* G4VEnergyLo 662 inline const G4ParticleDefinition* G4VEnergyLossProcess::Particle() const 841 { 663 { 842 return particle; 664 return particle; 843 } 665 } 844 666 845 //....oooOO0OOooo........oooOO0OOooo........oo 667 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 846 668 847 inline const G4ParticleDefinition* G4VEnergyLo 669 inline const G4ParticleDefinition* G4VEnergyLossProcess::BaseParticle() const 848 { 670 { 849 return baseParticle; 671 return baseParticle; 850 } 672 } 851 673 852 //....oooOO0OOooo........oooOO0OOooo........oo 674 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 853 675 854 inline const G4ParticleDefinition* << 676 inline const G4ParticleDefinition* G4VEnergyLossProcess::SecondaryParticle() const 855 G4VEnergyLossProcess::SecondaryParticle() cons << 856 { 677 { 857 return secondaryParticle; 678 return secondaryParticle; 858 } 679 } 859 680 860 //....oooOO0OOooo........oooOO0OOooo........oo 681 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 861 682 862 inline void G4VEnergyLossProcess::SetLossFluct << 683 inline G4VSubCutoffProcessor* G4VEnergyLossProcess::SubCutoffProcessor(size_t index) 863 { 684 { 864 lossFluctuationFlag = val; << 685 G4VSubCutoffProcessor* p = 0; 865 actLossFluc = true; << 686 if( nSCoffRegions ) p = scoffProcessors[idxSCoffRegions[index]]; 866 } << 687 return p; 867 << 868 //....oooOO0OOooo........oooOO0OOooo........oo << 869 << 870 inline void G4VEnergyLossProcess::SetSpline(G4 << 871 { << 872 spline = val; << 873 } << 874 << 875 //....oooOO0OOooo........oooOO0OOooo........oo << 876 << 877 inline void G4VEnergyLossProcess::SetCrossSect << 878 { << 879 fXSType = val; << 880 } << 881 << 882 //....oooOO0OOooo........oooOO0OOooo........oo << 883 << 884 inline G4CrossSectionType G4VEnergyLossProcess << 885 { << 886 return fXSType; << 887 } << 888 << 889 //....oooOO0OOooo........oooOO0OOooo........oo << 890 << 891 inline G4bool G4VEnergyLossProcess::IsIonisati << 892 { << 893 return isIonisation; << 894 } << 895 << 896 //....oooOO0OOooo........oooOO0OOooo........oo << 897 << 898 inline G4int G4VEnergyLossProcess::NumberOfSub << 899 { << 900 return nSCoffRegions; << 901 } << 902 << 903 //....oooOO0OOooo........oooOO0OOooo........oo << 904 << 905 inline G4double G4VEnergyLossProcess::MinKinEn << 906 { << 907 return minKinEnergy; << 908 } << 909 << 910 //....oooOO0OOooo........oooOO0OOooo........oo << 911 << 912 inline G4double G4VEnergyLossProcess::MaxKinEn << 913 { << 914 return maxKinEnergy; << 915 } << 916 << 917 //....oooOO0OOooo........oooOO0OOooo........oo << 918 << 919 inline G4double G4VEnergyLossProcess::CrossSec << 920 { << 921 return biasFactor; << 922 } << 923 << 924 //....oooOO0OOooo........oooOO0OOooo........oo << 925 << 926 inline G4bool G4VEnergyLossProcess::TablesAreB << 927 { << 928 return tablesAreBuilt; << 929 } 688 } 930 689 931 //....oooOO0OOooo........oooOO0OOooo........oo 690 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 932 691 933 inline G4PhysicsTable* G4VEnergyLossProcess::D 692 inline G4PhysicsTable* G4VEnergyLossProcess::DEDXTable() const 934 { 693 { 935 return theDEDXTable; 694 return theDEDXTable; 936 } 695 } 937 696 938 //....oooOO0OOooo........oooOO0OOooo........oo 697 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 939 698 940 inline G4PhysicsTable* G4VEnergyLossProcess::D 699 inline G4PhysicsTable* G4VEnergyLossProcess::DEDXunRestrictedTable() const 941 { 700 { 942 return theDEDXunRestrictedTable; 701 return theDEDXunRestrictedTable; 943 } 702 } 944 703 945 //....oooOO0OOooo........oooOO0OOooo........oo 704 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 946 705 947 inline G4PhysicsTable* G4VEnergyLossProcess::I << 706 inline G4PhysicsTable* G4VEnergyLossProcess::PreciseRangeTable() const 948 { 707 { 949 return theIonisationTable; << 708 return thePreciseRangeTable; 950 } << 951 << 952 //....oooOO0OOooo........oooOO0OOooo........oo << 953 << 954 inline G4PhysicsTable* G4VEnergyLossProcess::C << 955 { << 956 return theCSDARangeTable; << 957 } 709 } 958 710 959 //....oooOO0OOooo........oooOO0OOooo........oo 711 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 960 712 961 inline G4PhysicsTable* G4VEnergyLossProcess::R 713 inline G4PhysicsTable* G4VEnergyLossProcess::RangeTableForLoss() const 962 { 714 { 963 return theRangeTableForLoss; 715 return theRangeTableForLoss; 964 } 716 } 965 717 966 //....oooOO0OOooo........oooOO0OOooo........oo 718 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 967 719 968 inline G4PhysicsTable* G4VEnergyLossProcess::I 720 inline G4PhysicsTable* G4VEnergyLossProcess::InverseRangeTable() const 969 { 721 { 970 return theInverseRangeTable; 722 return theInverseRangeTable; 971 } 723 } 972 724 973 //....oooOO0OOooo........oooOO0OOooo........oo 725 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 974 726 975 inline G4PhysicsTable* G4VEnergyLossProcess::L << 727 inline G4PhysicsTable* G4VEnergyLossProcess::LambdaTable() 976 { 728 { 977 return theLambdaTable; 729 return theLambdaTable; 978 } 730 } 979 731 980 //....oooOO0OOooo........oooOO0OOooo........oo 732 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 981 733 982 inline G4bool G4VEnergyLossProcess::UseBaseMat << 734 inline G4PhysicsTable* G4VEnergyLossProcess::SubLambdaTable() 983 { 735 { 984 return baseMat; << 736 return theSubLambdaTable; 985 } 737 } 986 738 987 //....oooOO0OOooo........oooOO0OOooo........oo 739 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 988 << 740 989 inline std::vector<G4double>* << 741 inline G4bool G4VEnergyLossProcess::IsIntegral() const 990 G4VEnergyLossProcess::EnergyOfCrossSectionMax( << 991 { 742 { 992 return theEnergyOfCrossSectionMax; << 743 return integral; 993 } 744 } 994 745 995 //....oooOO0OOooo........oooOO0OOooo........oo 746 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 996 747 997 inline std::vector<G4TwoPeaksXS*>* G4VEnergyLo << 748 inline size_t G4VEnergyLossProcess::CurrentMaterialCutsCoupleIndex() const 998 { 749 { 999 return fXSpeaks; << 750 return currentMaterialIndex; 1000 } 751 } 1001 752 1002 //....oooOO0OOooo........oooOO0OOooo........o 753 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1003 << 754 1004 inline std::size_t G4VEnergyLossProcess::Numb << 755 inline void G4VEnergyLossProcess::SetMassRatio(G4double val) 1005 { 756 { 1006 return numberOfModels; << 757 massRatio = val; 1007 } 758 } 1008 759 1009 //....oooOO0OOooo........oooOO0OOooo........o 760 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 761 >> 762 inline void G4VEnergyLossProcess::SetReduceFactor(G4double val) >> 763 { >> 764 reduceFactor = val; >> 765 } 1010 766 1011 inline G4VEmModel* G4VEnergyLossProcess::EmMo << 767 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 768 >> 769 inline void G4VEnergyLossProcess::SetChargeSquare(G4double val) 1012 { 770 { 1013 return (index < emModels->size()) ? (*emMod << 771 chargeSquare = val; 1014 } 772 } 1015 773 1016 //....oooOO0OOooo........oooOO0OOooo........o 774 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1017 775 1018 inline G4VEmModel* << 776 inline void G4VEnergyLossProcess::SetChargeSquareRatio(G4double val) 1019 G4VEnergyLossProcess::GetModelByIndex(std::si << 1020 { 777 { 1021 return modelManager->GetModel((G4int)idx, v << 778 chargeSqRatio = val; >> 779 } >> 780 >> 781 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 782 >> 783 inline G4double G4VEnergyLossProcess::GetCurrentRange() const >> 784 { >> 785 return fRange; 1022 } 786 } 1023 787 1024 //....oooOO0OOooo........oooOO0OOooo........o 788 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1025 789 1026 #endif 790 #endif 1027 791