Changeset 1196 for trunk/source/processes/electromagnetic/utils
- Timestamp:
- Nov 25, 2009, 5:13:58 PM (15 years ago)
- Location:
- trunk/source/processes/electromagnetic/utils
- Files:
-
- 48 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/electromagnetic/utils/History
r1055 r1196 1 $Id: History,v 1.3 83 2009/05/27 11:55:02vnivanch Exp $1 $Id: History,v 1.398 2009/11/11 23:59:48 vnivanch Exp $ 2 2 ------------------------------------------------------------------- 3 3 … … 17 17 * Reverse chronological order (last date on top), please * 18 18 ---------------------------------------------------------- 19 20 10 November 09: V.Ivant (emutils-V09-02-22) 21 G4EmCalculator - added method FindIon; improve computation for ions 22 allowing to be applied for the model based on ICRU'73 23 24 10 November 09: V.Ivant (emutils-V09-02-21) 25 G4VEmProcess - improved printout 26 27 29 September 09: V.Ivant (emutils-V09-02-20) 28 - G4VMultipleScattering - added SetModel method 29 - G4EmCorrections - fixed bug reported in phys-list forum #532 30 - G4ionEffectiveCharge - return back lost protection for zero energy 31 - G4LossTableManager, G4EmProcessOptions, G4EnergyLossMessenger - added 32 parameter FactorForAngleLimit, Set/Get methods and UI command 33 to be used computation of limit on -t (invariant momentum transfer) 34 inside single scattering and G4WentzelVI models 35 36 29 September 09: V.Ivant (emutils-V09-02-19) 37 - G4EmElementSelector - fixed forgotten migration to updated G4PhysicsVector 38 39 25 September 09: V.Ivant (emutils-V09-02-18) 40 - G4VEmModel - insure definition of pointer to the G4MaterialCutsCouple 41 object both at initialisation and in run time 42 - G4EmSaturation - use PDG encoding instead of pointer to G4ParticleDefinition 43 44 11 August 09: V.Ivant (emutils-V09-02-17) 45 - G4EmModelManager - reduced length of internal arrays, simplified 46 initialisation, in particular, smoothing procedure, the size of 47 executable should be reduced 48 - G4VEmProcess, G4VMultipleScattering, G4VEnergyLossProcess, 49 G4LossTableBuilder - use copy constructors to create Physics Vectors in 50 order to reduce number of calls to log function 51 52 30 July 09: V.Ivant (emutils-V09-02-16) 53 - G4VEmProcess, G4VMultipleScattering, G4VEnergyLossProcess - fixed 54 bug in RetrieveTables - check that each its PhysicsVector was 55 retrieved before set of SplineFlag 56 57 24 July 09: V.Ivant (emutils-V09-02-15) 58 - G4AtomicShell class moved from lowenergy 59 - G4VEmProcess, G4VMultipleScattering, G4VEnergyLossProcess - added 60 initialisation of polarAngleLimit and highEnergyLimit parameters 61 for all models 62 63 22 July 09: V.Ivant (emutils-V09-02-14) 64 - G4VEmProcess - modified method SelectModel required for the 65 G4NuclearStopping process 66 67 20 July 09: V.Ivant (emutils-V09-02-13) 68 - G4VMultipleScattering - added initialisation of generic msc model 69 parameters in this base class (allowing 70 to overwrite default models) 71 - G4EmConfigurator - cleanup 72 - G4VMscModel - set facsafety=0.3 as it is defined in Urban models 73 74 9 July 09: V.Ivant (emutils-V09-02-12) 75 - G4VEnergyLossProcess, G4VEmProcess, G4VMultipleScattering, 76 G4LossTableBuilder - used updated G4PhysicsVector (no hidden bin anymore); 77 used Energy() and Value() methods instead of 78 GetLowEdgeEnergy() and GetValue() 79 - G4VEnergyLossProcess - fixed retrieve from ASCII files 80 - G4VMultipleScattering - remove method obsolete GetMscContinuesStepLimit 81 - G4VAtomDeexcitation - a new header file 19 82 20 83 27 May 09: V.Ivant (emutils-V09-02-11) -
trunk/source/processes/electromagnetic/utils/include/G4DummyModel.hh
r1055 r1196 25 25 // 26 26 // $Id: G4DummyModel.hh,v 1.4 2009/04/07 18:39:47 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/include/G4ElectronIonPair.hh
r1007 r1196 25 25 // 26 26 // $Id: G4ElectronIonPair.hh,v 1.2 2008/10/17 14:46:16 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // -
trunk/source/processes/electromagnetic/utils/include/G4EmCalculator.hh
r1007 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmCalculator.hh,v 1.1 8 2007/03/15 12:34:46vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4EmCalculator.hh,v 1.19 2009/11/11 23:59:48 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // … … 60 60 #include "globals.hh" 61 61 #include "G4DataVector.hh" 62 #include "G4DynamicParticle.hh" 62 63 63 64 class G4LossTableManager; … … 82 83 ~G4EmCalculator(); 83 84 85 //================================================================================== 84 86 // Methods to access precalculated dE/dx and cross sections 85 87 // Materials should exist in the list of the G4MaterialCutsCouple 88 //================================================================================== 86 89 87 90 G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition*, const G4Material*, … … 136 139 void PrintInverseRangeTable(const G4ParticleDefinition*); 137 140 141 //================================================================================== 138 142 // Methods to calculate dE/dx and cross sections "on fly" 139 143 // Existing tables and G4MaterialCutsCouples are not used 144 //================================================================================== 140 145 141 146 G4double ComputeDEDX(G4double kinEnergy, const G4ParticleDefinition*, … … 150 155 const G4String& mat, G4double cut = DBL_MAX); 151 156 152 G4double ComputeNuclearDEDX(G4double kinEnergy, const G4ParticleDefinition*, const G4Material*); 153 G4double ComputeNuclearDEDX(G4double kinEnergy, const G4String& part, const G4String& mat); 157 G4double ComputeNuclearDEDX(G4double kinEnergy, const G4ParticleDefinition*, 158 const G4Material*); 159 G4double ComputeNuclearDEDX(G4double kinEnergy, const G4String& part, 160 const G4String& mat); 154 161 155 162 G4double ComputeTotalDEDX(G4double kinEnergy, const G4ParticleDefinition*, … … 190 197 const G4String&); 191 198 199 //================================================================================== 200 // Methods to access particles, materials, regions 201 //================================================================================== 202 192 203 const G4ParticleDefinition* FindParticle(const G4String&); 193 204 205 const G4ParticleDefinition* FindIon(G4int Z, G4int A); 206 194 207 const G4Material* FindMaterial(const G4String&); 195 208 … … 199 212 200 213 void SetVerbose(G4int val); 214 215 //================================================================================== 216 // Private methods 217 //================================================================================== 201 218 202 219 private: … … 240 257 const G4ParticleDefinition* theGenericIon; 241 258 G4ionEffectiveCharge* ionEffCharge; 259 G4DynamicParticle dynParticle; 242 260 243 261 G4String currentName; -
trunk/source/processes/electromagnetic/utils/include/G4EmConfigurator.hh
r1007 r1196 25 25 // 26 26 // $Id: G4EmConfigurator.hh,v 1.2 2008/11/21 12:30:29 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/include/G4EmCorrections.hh
r1007 r1196 25 25 // 26 26 // $Id: G4EmCorrections.hh,v 1.24 2008/09/12 14:44:48 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/include/G4EmElementSelector.hh
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmElementSelector.hh,v 1. 4 2009/05/10 18:26:43vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4EmElementSelector.hh,v 1.6 2009/09/29 11:31:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 109 109 { 110 110 if (nElmMinusOne > 0) { 111 G4bool b;112 111 G4double x = G4UniformRand(); 113 for(G4int i=0; i<nElmMinusOne; i++) { 114 if (x <= (xSections[i])->GetValue(e,b)) { 112 element = (*theElementVector)[nElmMinusOne]; 113 for(G4int i=0; i<nElmMinusOne; ++i) { 114 if (x <= (xSections[i])->Value(e)) { 115 115 element = (*theElementVector)[i]; 116 116 break; -
trunk/source/processes/electromagnetic/utils/include/G4EmModelManager.hh
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmModelManager.hh,v 1. 28 2009/04/09 15:53:17vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4EmModelManager.hh,v 1.34 2009/08/11 10:29:30 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 54 54 // 15-03-07 Add maxCutInRange (V.Ivanchenko) 55 55 // 08-04-08 Simplify Select method for only one G4RegionModel (VI) 56 // 03-08-09 Removed unused members and simplify model search if only one 57 // model is used (VI) 56 58 // 57 59 // Class Description: … … 96 98 if (nModelsForRegion>1) { 97 99 idx = nModelsForRegion; 98 do { idx--;} while (idx&& e <= lowKineticEnergy[idx]);100 do {--idx;} while (idx > 0 && e <= lowKineticEnergy[idx]); 99 101 } 100 102 return theListOfModelIndexes[idx]; … … 130 132 class G4Region; 131 133 class G4ParticleDefinition; 132 class G4DataVector;133 134 class G4PhysicsVector; 134 135 class G4MaterialCutsCouple; … … 174 175 175 176 private: 177 178 inline G4double ComputeDEDX(G4VEmModel* model, 179 const G4MaterialCutsCouple*, 180 G4double kinEnergy, 181 G4double cutEnergy, 182 G4double minEnergy); 176 183 177 184 // hide assignment operator … … 195 202 G4int nEmModels; 196 203 G4int nRegions; 197 G4int nCouples; 198 199 G4int* idxOfRegionModels; 200 G4RegionModels** setOfRegionModels; 201 202 G4double minSubRange; 203 G4double maxCutInRange; 204 205 std::vector<G4int> idxOfRegionModels; 206 std::vector<G4RegionModels*> setOfRegionModels; 207 204 208 G4double maxSubCutInRange; 205 209 206 210 const G4ParticleDefinition* particle; 207 const G4ParticleDefinition* secondaryParticle;208 const G4ParticleDefinition* theGamma;209 const G4ParticleDefinition* thePositron;210 211 211 212 G4int verboseLevel; 213 G4bool severalModels; 212 214 213 215 // may be changed in run time 214 216 G4RegionModels* currRegionModel; 217 G4VEmModel* currModel; 215 218 }; 216 219 … … 221 224 size_t& index) 222 225 { 223 if(nRegions > 1) currRegionModel = setOfRegionModels[idxOfRegionModels[index]]; 224 return models[currRegionModel->SelectIndex(kinEnergy)]; 226 if(severalModels) { 227 if(nRegions > 1) { 228 currRegionModel = setOfRegionModels[idxOfRegionModels[index]]; 229 } 230 currModel = models[currRegionModel->SelectIndex(kinEnergy)]; 231 } 232 return currModel; 225 233 } 226 234 … … 248 256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 249 257 258 inline G4double 259 G4EmModelManager::ComputeDEDX(G4VEmModel* model, 260 const G4MaterialCutsCouple* couple, 261 G4double e, 262 G4double cut, 263 G4double emin) 264 { 265 G4double dedx = 0.0; 266 if(model && cut > emin) { 267 dedx = model->ComputeDEDX(couple,particle,e,cut); 268 if(emin > 0.0) {dedx -= model->ComputeDEDX(couple,particle,e,emin);} 269 } 270 return dedx; 271 } 272 273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 274 250 275 #endif 251 276 -
trunk/source/processes/electromagnetic/utils/include/G4EmMultiModel.hh
r1007 r1196 25 25 // 26 26 // $Id: G4EmMultiModel.hh,v 1.6 2007/05/22 17:31:57 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/include/G4EmProcessOptions.hh
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmProcessOptions.hh,v 1.1 5 2009/02/18 14:40:10vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4EmProcessOptions.hh,v 1.16 2009/10/29 19:25:28 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // … … 128 128 void SetPolarAngleLimit(G4double val); 129 129 130 void SetFactorForAngleLimit(G4double val); 131 130 132 private: 131 133 -
trunk/source/processes/electromagnetic/utils/include/G4EmProcessSubType.hh
r1007 r1196 26 26 // 27 27 // $Id: G4EmProcessSubType.hh,v 1.9 2008/12/18 13:01:42 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 //--------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/include/G4EmSaturation.hh
r1007 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmSaturation.hh,v 1. 6 2008/03/17 11:27:32vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4EmSaturation.hh,v 1.7 2009/09/25 09:16:40 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // … … 104 104 void Initialise(); 105 105 106 const G4ParticleDefinition* gamma;107 106 const G4ParticleDefinition* electron; 108 const G4ParticleDefinition* neutron;109 107 const G4ParticleDefinition* proton; 110 108 G4LossTableManager* manager; -
trunk/source/processes/electromagnetic/utils/include/G4EmTableType.hh
r1007 r1196 26 26 // 27 27 // $Id: G4EmTableType.hh,v 1.4 2008/10/13 14:56:56 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 //--------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/include/G4EnergyLossMessenger.hh
r1055 r1196 25 25 // 26 26 // 27 // $Id: G4EnergyLossMessenger.hh,v 1.2 3 2009/02/18 14:40:10vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // $Id: G4EnergyLossMessenger.hh,v 1.24 2009/10/29 19:25:28 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 126 126 G4UIcmdWithADouble* frCmd; 127 127 G4UIcmdWithADouble* fgCmd; 128 G4UIcmdWithADouble* mscfCmd; 128 129 G4UIcmdWithADoubleAndUnit* angCmd; 129 130 }; -
trunk/source/processes/electromagnetic/utils/include/G4EnergyLossTables.hh
r1007 r1196 26 26 // 27 27 // $Id: G4EnergyLossTables.hh,v 1.19 2006/06/29 19:54:35 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // $Id: -
trunk/source/processes/electromagnetic/utils/include/G4LossTableBuilder.hh
r1007 r1196 25 25 // 26 26 // $Id: G4LossTableBuilder.hh,v 1.8 2008/07/22 15:55:15 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // -
trunk/source/processes/electromagnetic/utils/include/G4LossTableManager.hh
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4LossTableManager.hh,v 1.5 4 2009/04/09 16:10:57vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4LossTableManager.hh,v 1.55 2009/10/29 19:25:28 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // … … 210 210 void SetBremsstrahlungTh(G4double val); 211 211 212 void SetFactorForAngleLimit(G4double val); 213 212 214 void SetVerbose(G4int val); 213 215 … … 221 223 222 224 G4double BremsstrahlungTh() const; 225 226 G4double FactorForAngleLimit() const; 223 227 224 228 const std::vector<G4VEnergyLossProcess*>& GetEnergyLossProcessVector(); … … 302 306 G4double maxKinEnergyForMuons; 303 307 G4double bremsTh; 308 G4double factorForAngleLimit; 304 309 305 310 G4LossTableBuilder* tableBuilder; -
trunk/source/processes/electromagnetic/utils/include/G4MscStepLimitType.hh
r1007 r1196 26 26 // 27 27 // $Id: G4MscStepLimitType.hh,v 1.2 2007/06/11 14:56:51 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 //--------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/include/G4VEmFluctuationModel.hh
r1055 r1196 25 25 // 26 26 // $Id: G4VEmFluctuationModel.hh,v 1.12 2009/02/19 11:25:50 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/include/G4VEmModel.hh
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VEmModel.hh,v 1. 69 2009/05/26 15:00:49vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4VEmModel.hh,v 1.72 2009/09/23 14:42:47 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 137 137 G4double maxEnergy = DBL_MAX); 138 138 139 // main method to compute cross section depending onatom139 // main method to compute cross section per atom 140 140 virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 141 141 G4double kinEnergy, … … 159 159 G4double kineticEnergy); 160 160 161 // add correction to energy loss and ompute non-ionizing energy loss161 // add correction to energy loss and compute non-ionizing energy loss 162 162 virtual void CorrectionsAlongStep(const G4MaterialCutsCouple*, 163 163 const G4DynamicParticle*, … … 171 171 G4double& eloss); 172 172 173 // add region for the model 174 virtual void DefineForRegion(const G4Region*); 175 176 // initilisation at run time for given material 173 // initilisation at run time for a given material 177 174 virtual void SetupForMaterial(const G4ParticleDefinition*, 178 175 const G4Material*, 179 176 G4double kineticEnergy); 177 178 // add a region for the model 179 virtual void DefineForRegion(const G4Region*); 180 180 181 181 protected: … … 209 209 // cross section per volume 210 210 inline G4double CrossSection(const G4MaterialCutsCouple*, 211 212 213 214 211 const G4ParticleDefinition*, 212 G4double kineticEnergy, 213 G4double cutEnergy = 0.0, 214 G4double maxEnergy = DBL_MAX); 215 215 216 216 // compute mean free path via cross section per volume … … 267 267 inline void SetLowEnergyLimit(G4double); 268 268 269 inline void SetActivationHighEnergyLimit(G4double); 270 271 inline void SetActivationLowEnergyLimit(G4double); 272 273 inline G4bool IsActive(G4double kinEnergy); 274 269 275 inline void SetPolarAngleLimit(G4double); 270 276 … … 308 314 G4double lowLimit; 309 315 G4double highLimit; 316 G4double eMinActive; 317 G4double eMaxActive; 310 318 G4double polarAngleLimit; 311 319 G4double secondaryThreshold; … … 434 442 inline G4int G4VEmModel::SelectIsotopeNumber(const G4Element* elm) 435 443 { 444 currentElement = elm; 436 445 G4int N = G4int(elm->GetN() + 0.5); 437 446 G4int ni = elm->GetNumberOfIsotopes(); … … 439 448 G4int idx = 0; 440 449 if(ni > 1) { 441 G4double* ab = currentElement->GetRelativeAbundanceVector();450 G4double* ab = elm->GetRelativeAbundanceVector(); 442 451 G4double x = G4UniformRand(); 443 452 for(; idx<ni; idx++) { … … 517 526 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 518 527 528 inline void G4VEmModel::SetActivationHighEnergyLimit(G4double val) 529 { 530 eMaxActive = val; 531 } 532 533 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 534 535 inline void G4VEmModel::SetActivationLowEnergyLimit(G4double val) 536 { 537 eMinActive = val; 538 } 539 540 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 541 542 inline G4bool G4VEmModel::IsActive(G4double kinEnergy) 543 { 544 return (kinEnergy >= eMinActive && kinEnergy <= eMaxActive); 545 } 546 547 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 548 519 549 inline void G4VEmModel::SetPolarAngleLimit(G4double val) 520 550 { -
trunk/source/processes/electromagnetic/utils/include/G4VEmProcess.hh
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VEmProcess.hh,v 1.5 1 2009/04/07 18:39:47 vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4VEmProcess.hh,v 1.55 2009/09/23 14:42:47 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 202 202 protected: 203 203 // Select model in run time 204 inline void SelectModel(G4double& kinEnergy);204 inline G4VEmModel* SelectModel(G4double& kinEnergy, size_t index); 205 205 206 206 public: … … 344 344 const G4Material* currentMaterial; 345 345 const G4MaterialCutsCouple* currentCouple; 346 size_t current MaterialIndex;346 size_t currentCoupleIndex; 347 347 348 348 G4double mfpKinEnergy; … … 358 358 G4double kineticEnergy, G4double Z, G4double A, G4double cut) 359 359 { 360 SelectModel(kineticEnergy );360 SelectModel(kineticEnergy, currentCoupleIndex); 361 361 G4double x = 0.0; 362 362 if(currentModel) { … … 468 468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 469 469 470 inline void G4VEmProcess::SelectModel(G4double& kinEnergy) 471 { 472 currentModel = modelManager->SelectModel(kinEnergy, currentMaterialIndex); 470 inline 471 G4VEmModel* G4VEmProcess::SelectModel(G4double& kinEnergy, size_t index) 472 { 473 currentModel = modelManager->SelectModel(kinEnergy, index); 473 474 currentModel->SetCurrentCouple(currentCouple); 474 } 475 476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 477 478 inline G4VEmModel* G4VEmProcess::SelectModelForMaterial( 479 G4double kinEnergy, size_t& idxRegion) const 475 return currentModel; 476 } 477 478 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 479 480 inline 481 G4VEmModel* G4VEmProcess::SelectModelForMaterial(G4double kinEnergy, 482 size_t& idxRegion) const 480 483 { 481 484 return modelManager->SelectModel(kinEnergy, idxRegion); … … 545 548 inline size_t G4VEmProcess::CurrentMaterialCutsCoupleIndex() const 546 549 { 547 return current MaterialIndex;550 return currentCoupleIndex; 548 551 } 549 552 … … 552 555 inline G4double G4VEmProcess::GetGammaEnergyCut() 553 556 { 554 return (*theCutsGamma)[current MaterialIndex];557 return (*theCutsGamma)[currentCoupleIndex]; 555 558 } 556 559 … … 559 562 inline G4double G4VEmProcess::GetElectronEnergyCut() 560 563 { 561 return (*theCutsElectron)[current MaterialIndex];564 return (*theCutsElectron)[currentCoupleIndex]; 562 565 } 563 566 … … 583 586 preStepKinEnergy = track.GetKineticEnergy(); 584 587 DefineMaterial(track.GetMaterialCutsCouple()); 588 SelectModel(preStepKinEnergy, currentCoupleIndex); 585 589 if (theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX; 586 590 } … … 593 597 currentCouple = couple; 594 598 currentMaterial = couple->GetMaterial(); 595 current MaterialIndex = couple->GetIndex();599 currentCoupleIndex = couple->GetIndex(); 596 600 mfpKinEnergy = DBL_MAX; 597 601 } … … 602 606 inline void G4VEmProcess::ComputeIntegralLambda(G4double e) 603 607 { 604 mfpKinEnergy = theEnergyOfCrossSectionMax[current MaterialIndex];608 mfpKinEnergy = theEnergyOfCrossSectionMax[currentCoupleIndex]; 605 609 if (e <= mfpKinEnergy) { 606 610 preStepLambda = GetLambdaFromTable(e); … … 616 620 } 617 621 } else { 618 preStepLambda = theCrossSectionMax[current MaterialIndex];622 preStepLambda = theCrossSectionMax[currentCoupleIndex]; 619 623 } 620 624 } … … 625 629 inline G4double G4VEmProcess::GetLambdaFromTable(G4double e) 626 630 { 627 G4bool b; 628 return (((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b)); 631 return (((*theLambdaTable)[currentCoupleIndex])->Value(e)); 629 632 } 630 633 … … 634 637 { 635 638 G4double x = 0.0; 636 if(theLambdaTable) x = GetLambdaFromTable(e);637 else x = ComputeCurrentLambda(e);639 if(theLambdaTable) { x = GetLambdaFromTable(e); } 640 else { x = ComputeCurrentLambda(e); } 638 641 return x; 639 642 } … … 643 646 inline G4double G4VEmProcess::ComputeCurrentLambda(G4double e) 644 647 { 645 SelectModel(e); 646 G4double x = 0.0; 647 if(currentModel) { 648 x = currentModel->CrossSectionPerVolume(currentMaterial,particle, 649 e,(*theCuts)[currentMaterialIndex]); 650 } 651 return x; 648 SelectModel(e, currentCoupleIndex); 649 return currentModel->CrossSectionPerVolume(currentMaterial,particle, 650 e,(*theCuts)[currentCoupleIndex]); 652 651 } 653 652 -
trunk/source/processes/electromagnetic/utils/include/G4VEnergyLoss.hh
r1007 r1196 26 26 // 27 27 // $Id: G4VEnergyLoss.hh,v 1.18 2006/06/29 19:54:47 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // -
trunk/source/processes/electromagnetic/utils/include/G4VEnergyLossProcess.hh
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VEnergyLossProcess.hh,v 1.8 7 2009/04/07 18:39:47 vnivanch Exp $26 // $Id: G4VEnergyLossProcess.hh,v 1.89 2009/07/03 14:39:17 vnivanch Exp $ 27 27 // GEANT4 tag $Name: 28 28 // … … 355 355 inline void SetLambdaFactor(G4double val); 356 356 inline void SetStepFunction(G4double v1, G4double v2); 357 inline void SetLowestEnergyLimit(G4double); 357 358 358 359 inline G4int NumberOfSubCutoffRegions() const; … … 717 718 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 718 719 720 inline void G4VEnergyLossProcess::SetLowestEnergyLimit(G4double val) 721 { 722 lowestKinEnergy = val; 723 } 724 725 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 726 719 727 inline G4int G4VEnergyLossProcess::NumberOfSubCutoffRegions() const 720 728 { … … 958 966 { 959 967 G4double e = amu_c2*preStepKinEnergy/particle->GetPDGMass(); 960 G4bool b; 961 G4double s = fRange*std::pow(10.,vstrag->GetValue(e,b)); 968 G4double s = fRange*std::pow(10.,vstrag->Value(e)); 962 969 G4double x = fRange + G4RandGauss::shoot(0.0,s); 963 970 if(x > 0.0) fRange = x; … … 992 999 inline G4double G4VEnergyLossProcess::GetDEDXForScaledEnergy(G4double e) 993 1000 { 994 G4bool b; 995 G4double x = 996 ((*theDEDXTable)[currentMaterialIndex]->GetValue(e, b))*chargeSqRatio; 1001 G4double x = ((*theDEDXTable)[currentMaterialIndex]->Value(e))*chargeSqRatio; 997 1002 if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy); 998 1003 return x; … … 1003 1008 inline G4double G4VEnergyLossProcess::GetSubDEDXForScaledEnergy(G4double e) 1004 1009 { 1005 G4bool b; 1006 G4double x = 1007 ((*theDEDXSubTable)[currentMaterialIndex]->GetValue(e, b))*chargeSqRatio; 1010 G4double x = ((*theDEDXSubTable)[currentMaterialIndex]->Value(e))*chargeSqRatio; 1008 1011 if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy); 1009 1012 return x; … … 1014 1017 inline G4double G4VEnergyLossProcess::GetIonisationForScaledEnergy(G4double e) 1015 1018 { 1016 G4bool b; 1017 G4double x = 0.0; 1019 //G4double x = 0.0; 1018 1020 // if(theIonisationTable) { 1019 x = ((*theIonisationTable)[currentMaterialIndex]->GetValue(e, b)) 1020 *chargeSqRatio; 1021 G4double x = ((*theIonisationTable)[currentMaterialIndex]->Value(e))*chargeSqRatio; 1021 1022 if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy); 1022 1023 //} … … 1029 1030 G4double G4VEnergyLossProcess::GetSubIonisationForScaledEnergy(G4double e) 1030 1031 { 1031 G4bool b; 1032 G4double x = 0.0; 1032 // G4double x = 0.0; 1033 1033 //if(theIonisationSubTable) { 1034 x = ((*theIonisationSubTable)[currentMaterialIndex]->GetValue(e, b)) 1035 *chargeSqRatio; 1034 G4double x = ((*theIonisationSubTable)[currentMaterialIndex]->Value(e))*chargeSqRatio; 1036 1035 if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy); 1037 1036 //} … … 1043 1042 inline G4double G4VEnergyLossProcess::GetScaledRangeForScaledEnergy(G4double e) 1044 1043 { 1045 G4bool b; 1046 G4double x = ((*theRangeTableForLoss)[currentMaterialIndex])->GetValue(e, b); 1044 G4double x = ((*theRangeTableForLoss)[currentMaterialIndex])->Value(e); 1047 1045 if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy); 1048 1046 return x; … … 1054 1052 G4double e) 1055 1053 { 1056 G4bool b;1057 1054 G4double x; 1058 1055 1059 1056 if (e < maxKinEnergyCSDA) { 1060 x = ((*theCSDARangeTable)[currentMaterialIndex])-> GetValue(e, b);1057 x = ((*theCSDARangeTable)[currentMaterialIndex])->Value(e); 1061 1058 if(e < minKinEnergy) x *= std::sqrt(e/minKinEnergy); 1062 1059 } else { … … 1072 1069 { 1073 1070 G4PhysicsVector* v = (*theInverseRangeTable)[currentMaterialIndex]; 1074 G4double rmin = v-> GetLowEdgeEnergy(0);1071 G4double rmin = v->Energy(0); 1075 1072 G4double e = 0.0; 1076 if(r >= rmin) { 1077 G4bool b; 1078 e = v->GetValue(r, b); 1079 } else if(r > 0.0) { 1073 if(r >= rmin) { e = v->Value(r); } 1074 else if(r > 0.0) { 1080 1075 G4double x = r/rmin; 1081 1076 e = minKinEnergy*x*x; … … 1088 1083 inline G4double G4VEnergyLossProcess::GetLambdaForScaledEnergy(G4double e) 1089 1084 { 1090 G4bool b; 1091 return 1092 chargeSqRatio*(((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b)); 1085 return chargeSqRatio*(((*theLambdaTable)[currentMaterialIndex])->Value(e)); 1093 1086 } 1094 1087 -
trunk/source/processes/electromagnetic/utils/include/G4VMscModel.hh
r1055 r1196 25 25 // 26 26 // $Id: G4VMscModel.hh,v 1.9 2009/04/07 18:39:47 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/include/G4VMultipleScattering.hh
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VMultipleScattering.hh,v 1. 56 2009/04/07 18:39:47vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4VMultipleScattering.hh,v 1.62 2009/10/29 17:56:04 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 127 127 void PrintInfoDefinition(); 128 128 129 G4VParticleChange* AlongStepDoIt(const G4Track&, const G4Step&);130 131 G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&);132 133 129 // Store PhysicsTable in a file. 134 130 // Return false in case of failure at I/O … … 158 154 // The function overloads the corresponding function of the base 159 155 // class. 160 G4double PostStepGetPhysicalInteractionLength(156 inline G4double PostStepGetPhysicalInteractionLength( 161 157 const G4Track&, 162 158 G4double previousStepSize, 163 159 G4ForceCondition* condition); 160 161 // Along step actions 162 inline G4VParticleChange* AlongStepDoIt(const G4Track&, const G4Step&); 163 164 // Post step actions 165 inline G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&); 164 166 165 167 // This method does not used for tracking, it is intended only for tests … … 209 211 void AddEmModel(G4int order, G4VEmModel*, const G4Region* region = 0); 210 212 213 // Assign a model to a process 214 void SetModel(G4VMscModel*, G4int index = 1); 215 216 // return the assigned model 217 G4VMscModel* Model(G4int index = 1); 218 211 219 // Access to models by index 212 G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false) ;220 G4VEmModel* GetModelByIndex(G4int idx = 0, G4bool ver = false) const; 213 221 214 222 //------------------------------------------------------------------------ … … 255 263 G4double& kineticEnergy); 256 264 257 // This method is used for tracking, it returns step limit258 inline G4double GetMscContinuousStepLimit(const G4Track& track,259 G4double scaledKinEnergy,260 G4double currentMinimalStep,261 G4double& currentSafety);262 263 265 // defines current material in run time 264 266 inline void DefineMaterial(const G4MaterialCutsCouple* couple); … … 278 280 279 281 // ======== Parameters of the class fixed at initialisation ======= 282 283 std::vector<G4VMscModel*> mscModels; 280 284 281 285 G4PhysicsTable* theLambdaTable; … … 294 298 295 299 G4bool latDisplasment; 300 G4bool isIon; 296 301 297 302 // ======== Cashed values - may be state dependent ================ … … 322 327 G4double& currentSafety) 323 328 { 324 return Get MscContinuousStepLimit(track,previousStepSize,currentMinimalStep,325 329 return GetContinuousStepLimit(track,previousStepSize,currentMinimalStep, 330 currentSafety); 326 331 } 327 332 … … 500 505 G4double x; 501 506 if(theLambdaTable) { 502 G4bool b; 503 x = ((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b); 507 x = ((*theLambdaTable)[currentMaterialIndex])->Value(e); 504 508 } else { 505 509 x = currentModel->CrossSection(currentCouple,p,e); … … 512 516 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 513 517 514 inline G4double G4VMultipleScattering::GetMscContinuousStepLimit(515 const G4Track& track,516 G4double scaledKinEnergy,517 G4double currentMinimalStep,518 G4double&)519 {520 G4double x = currentMinimalStep;521 DefineMaterial(track.GetMaterialCutsCouple());522 currentModel = static_cast<G4VMscModel*>(SelectModel(scaledKinEnergy));523 if(x > 0.0 && scaledKinEnergy > 0.0) {524 G4double tPathLength =525 currentModel->ComputeTruePathLengthLimit(track, theLambdaTable, x);526 if (tPathLength < x) valueGPILSelectionMSC = CandidateForSelection;527 x = currentModel->ComputeGeomPathLength(tPathLength);528 // G4cout << "tPathLength= " << tPathLength529 // << " stepLimit= " << x530 // << " currentMinimalStep= " << currentMinimalStep<< G4endl;531 }532 return x;533 }534 535 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....536 537 518 inline 538 519 void G4VMultipleScattering::DefineMaterial(const G4MaterialCutsCouple* couple) … … 553 534 554 535 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 536 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 537 538 // Follwoing methods are virtual, they are inlined because they applied at 539 // each simulation step and some compilers may inline these methods 540 541 inline G4double 542 G4VMultipleScattering::PostStepGetPhysicalInteractionLength( 543 const G4Track&, G4double, G4ForceCondition* condition) 544 { 545 *condition = Forced; 546 return DBL_MAX; 547 } 548 549 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 550 551 inline G4VParticleChange* 552 G4VMultipleScattering::AlongStepDoIt(const G4Track& track, const G4Step& step) 553 { 554 if(currentModel->IsActive(track.GetKineticEnergy())) { 555 fParticleChange.ProposeTrueStepLength(currentModel->ComputeTrueStepLength(step.GetStepLength())); 556 } 557 return &fParticleChange; 558 } 559 560 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 561 562 inline G4VParticleChange* 563 G4VMultipleScattering::PostStepDoIt(const G4Track& track, const G4Step& step) 564 { 565 fParticleChange.Initialize(track); 566 if(currentModel->IsActive(track.GetKineticEnergy())) { 567 currentModel->SampleScattering(track.GetDynamicParticle(), 568 step.GetPostStepPoint()->GetSafety()); 569 } 570 return &fParticleChange; 571 } 572 573 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 555 574 556 575 #endif -
trunk/source/processes/electromagnetic/utils/include/G4ionEffectiveCharge.hh
r1007 r1196 25 25 // 26 26 // $Id: G4ionEffectiveCharge.hh,v 1.12 2008/09/20 19:39:34 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/src/G4DummyModel.cc
r1055 r1196 25 25 // 26 26 // $Id: G4DummyModel.cc,v 1.4 2009/04/07 18:39:47 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/src/G4ElectronIonPair.cc
r1007 r1196 25 25 // 26 26 // $Id: G4ElectronIonPair.cc,v 1.2 2008/10/17 14:46:16 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/src/G4EmCalculator.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmCalculator.cc,v 1.4 6 2009/02/24 09:56:03vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4EmCalculator.cc,v 1.48 2009/11/11 23:59:48 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 127 127 if(couple && UpdateParticle(p, kinEnergy) ) { 128 128 res = manager->GetDEDX(p, kinEnergy, couple); 129 /* 129 130 if(isIon) { 130 G4double eth = 2.0*MeV/massRatio;131 if(kinEnergy > eth) { 132 G4double x1 = corr->ComputeIonCorrections(p,mat,kinEnergy); 133 G4double x2 = corr->ComputeIonCorrections(p,mat,eth);134 res += x1 - x2*eth/kinEnergy;135 /*136 G4cout << "### GetDEDX: E= " << kinEnergy << " res= " << res137 << " x1= " << x1 << " x2= " << x2138 << " del= " << x1 - x2*eth/kinEnergy << G4endl;;139 */140 } 141 } 142 131 if(FindEmModel(p, currentProcessName, kinEnergy)) { 132 G4double length = 1.*um; 133 G4cout << "### GetDEDX: E= " << kinEnergy << " res= " << res; 134 G4double eloss = res*length; 135 G4double niel = 0.0; 136 currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length); 137 res = eloss/length; 138 139 //G4cout << "### GetDEDX: E= " << kinEnergy << " res= " << res; 140 G4cout << " res1= " << res << G4endl;; 141 } 142 } 143 */ 143 144 if(verbose>0) { 144 145 G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV … … 147 148 << " " << p->GetParticleName() 148 149 << " in " << mat->GetName() 150 << " isIon= " << isIon 149 151 << G4endl; 150 152 } … … 307 309 FindLambdaTable(p, processName); 308 310 if(currentLambda) { 309 G4bool b;310 311 G4double e = kinEnergy*massRatio; 311 res = (((*currentLambda)[idx])-> GetValue(e,b))*chargeSquare;312 res = (((*currentLambda)[idx])->Value(e))*chargeSquare; 312 313 if(verbose>0) { 313 314 G4cout << "E(MeV)= " << kinEnergy/MeV … … 436 437 } 437 438 438 // emulate boundary region for different parameterisations439 // emulate smoothing procedure 439 440 G4double eth = currentModel->LowEnergyLimit(); 440 441 // G4cout << "massRatio= " << massRatio << " eth= " << eth << G4endl; 441 if(eth > 0.05*MeV && eth < 10.*MeV && escaled > eth && 442 loweModel != currentModel && loweModel) { 442 if(loweModel) { 443 443 G4double res0 = 0.0; 444 444 G4double res1 = 0.0; … … 462 462 << res1 << " q2= " << chargeSquare << G4endl; 463 463 */ 464 res *= (1.0 + (res0/res1 - 1.0)*eth/escaled); 465 466 if(isIon) { 467 G4double ethscaled = eth/massRatio; 468 if(kinEnergy > ethscaled) { 469 G4double x1 = corr->ComputeIonCorrections(p,mat,kinEnergy); 470 G4double x2 = corr->ComputeIonCorrections(p,mat,ethscaled); 471 res += x1 - x2*ethscaled/kinEnergy; 472 } 464 res += (res0 - res1)*eth/escaled; 465 } 466 467 // low energy correction for ions 468 /* 469 if(isIon) { 470 G4double length = 1.*nm; 471 const G4Region* r = 0; 472 const G4MaterialCutsCouple* couple = FindCouple(mat, r); 473 G4double eloss = res*length; 474 G4double niel = 0.0; 475 currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length); 476 res = eloss/length; 473 477 474 if(verbose > 1) { 475 G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV 476 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity()) 477 << G4endl; 478 } 478 if(verbose > 1) { 479 G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV 480 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity()) 481 << G4endl; 479 482 } 480 483 } 484 */ 485 } 481 486 482 if(verbose > 0) { 483 G4cout << "E(MeV)= " << kinEnergy/MeV 484 << " DEDX(MeV/mm)= " << res*mm/MeV 485 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity()) 486 << " cut(MeV)= " << cut/MeV 487 << " " << p->GetParticleName() 488 << " in " << currentMaterialName 489 << " Zi^2= " << chargeSquare 490 << G4endl; 491 } 487 if(verbose > 0) { 488 G4cout << "E(MeV)= " << kinEnergy/MeV 489 << " DEDX(MeV/mm)= " << res*mm/MeV 490 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity()) 491 << " cut(MeV)= " << cut/MeV 492 << " " << p->GetParticleName() 493 << " in " << currentMaterialName 494 << " Zi^2= " << chargeSquare 495 << G4endl; 492 496 } 493 497 } … … 745 749 { 746 750 if(p != currentParticle) { 751 752 // new particle 747 753 currentParticle = p; 754 dynParticle.SetDefinition(const_cast<G4ParticleDefinition*>(p)); 755 dynParticle.SetKineticEnergy(kinEnergy); 748 756 baseParticle = 0; 749 757 currentParticleName = p->GetParticleName(); … … 753 761 currentProcess = FindEnergyLossProcess(p); 754 762 currentProcessName = ""; 755 if(currentProcess) currentProcessName = currentProcess->GetProcessName(); 756 757 if(p->GetParticleType() == "nucleus" 758 && currentParticleName != "deuteron" 759 && currentParticleName != "triton" 760 && currentParticleName != "alpha+" 761 && currentParticleName != "helium" 762 && currentParticleName != "hydrogen" 763 ) { 764 baseParticle = theGenericIon; 765 massRatio = baseParticle->GetPDGMass()/p->GetPDGMass(); 766 isIon = true; 767 // G4cout << p->GetParticleName() 768 // << " in " << currentMaterial->GetName() 769 // << " e= " << kinEnergy << G4endl; 770 chargeSquare = 771 ionEffCharge->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy); 772 //G4cout << "q2= " << chargeSquare << G4endl; 773 } else { 774 isIon = false; 775 if(currentProcess) { 776 baseParticle = currentProcess->BaseParticle(); 777 778 if(baseParticle) { 779 massRatio = baseParticle->GetPDGMass()/p->GetPDGMass(); 780 G4double q = baseParticle->GetPDGCharge()/eplus; 781 chargeSquare /= (q*q); 782 } 783 } 784 } 785 } 763 isIon = false; 764 765 // ionisation process exist 766 if(currentProcess) { 767 currentProcessName = currentProcess->GetProcessName(); 768 baseParticle = currentProcess->BaseParticle(); 769 770 // base particle is used 771 if(baseParticle) { 772 massRatio = baseParticle->GetPDGMass()/p->GetPDGMass(); 773 G4double q = p->GetPDGCharge()/baseParticle->GetPDGCharge(); 774 chargeSquare = q*q; 775 } 776 777 if(p->GetParticleType() == "nucleus" 778 && currentParticleName != "deuteron" 779 && currentParticleName != "triton" 780 && currentParticleName != "alpha+" 781 && currentParticleName != "helium" 782 && currentParticleName != "hydrogen" 783 ) { 784 isIon = true; 785 massRatio = theGenericIon->GetPDGMass()/p->GetPDGMass(); 786 baseParticle = theGenericIon; 787 // G4cout << p->GetParticleName() 788 // << " in " << currentMaterial->GetName() 789 // << " e= " << kinEnergy << G4endl; 790 } 791 } 792 } 793 794 // Effective charge for ions 786 795 if(isIon) { 787 796 chargeSquare = 788 ionEffCharge->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy)797 corr->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy) 789 798 * corr->EffectiveChargeCorrection(p,currentMaterial,kinEnergy); 790 799 if(currentProcess) { … … 810 819 p = currentParticle; 811 820 } 821 return p; 822 } 823 824 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 825 826 const G4ParticleDefinition* G4EmCalculator::FindIon(G4int Z, G4int A) 827 { 828 const G4ParticleDefinition* p = 829 G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z); 812 830 return p; 813 831 } … … 897 915 G4String partname = p->GetParticleName(); 898 916 const G4ParticleDefinition* part = p; 899 if(isIon) part = theGenericIon; 900 917 if(isIon) { part = theGenericIon; } 918 919 // energy loss process 901 920 G4LossTableManager* lManager = G4LossTableManager::Instance(); 902 921 const std::vector<G4VEnergyLossProcess*> vel = 903 922 lManager->GetEnergyLossProcessVector(); 904 923 G4int n = vel.size(); 905 924 for(G4int i=0; i<n; i++) { 906 925 if((vel[i])->GetProcessName() == currentName && 907 926 (vel[i])->Particle() == part) 908 { 909 currentLambda = (vel[i])->LambdaTable(); 910 isApplicable = true; 911 break; 912 } 913 } 927 { 928 currentLambda = (vel[i])->LambdaTable(); 929 isApplicable = true; 930 break; 931 } 932 } 933 934 // discrete process 914 935 if(!currentLambda) { 915 936 const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector(); … … 925 946 } 926 947 } 948 949 // msc process 927 950 if(!currentLambda) { 928 951 const std::vector<G4VMultipleScattering*> vmsc = … … 957 980 const G4ParticleDefinition* part = p; 958 981 G4double scaledEnergy = kinEnergy*massRatio; 959 if(isIon) part = theGenericIon;982 if(isIon) { part = theGenericIon; } 960 983 961 984 if(verbose > 1) { … … 967 990 } 968 991 969 // Search for theprocess992 // Search for energy loss process 970 993 currentName = processName; 971 994 currentModel = 0; … … 976 999 lManager->GetEnergyLossProcessVector(); 977 1000 G4int n = vel.size(); 1001 G4VEnergyLossProcess* elproc = 0; 978 1002 for(G4int i=0; i<n; i++) { 979 1003 // G4cout << "i= " << i << " part= " 980 1004 // << (vel[i])->Particle()->GetParticleName() 981 1005 // << " proc= " << (vel[i])->GetProcessName() << G4endl; 982 if((vel[i])->GetProcessName() == currentName && 983 (vel[i])->Particle() == part) 984 { 985 const G4ParticleDefinition* bp = (vel[i])->BaseParticle(); 986 // G4cout << "i= " << i << " bp= " << bp << G4endl; 987 if(!bp) { 988 currentModel = (vel[i])->SelectModelForMaterial(scaledEnergy, idx); 989 loweModel = (vel[i])->SelectModelForMaterial(keV, idx); 990 isApplicable = true; 991 break; 1006 if((vel[i])->GetProcessName() == currentName) { 1007 if(baseParticle) { 1008 if((vel[i])->Particle() == baseParticle) { 1009 elproc = vel[i]; 1010 break; 1011 } 992 1012 } else { 993 for(G4int j=0; j<n; j++) { 994 if((vel[j])->Particle() == bp) { 995 currentModel = (vel[j])->SelectModelForMaterial(scaledEnergy, idx); 996 loweModel = (vel[j])->SelectModelForMaterial(keV, idx); 997 isApplicable = true; 998 break; 999 } 1013 if((vel[i])->Particle() == part) { 1014 elproc = vel[i]; 1015 break; 1000 1016 } 1001 1017 } 1002 1018 } 1003 1019 } 1020 if(elproc) { 1021 isApplicable = true; 1022 currentModel = elproc->SelectModelForMaterial(scaledEnergy, idx); 1023 G4double eth = currentModel->LowEnergyLimit(); 1024 if(eth > 0.0) { loweModel = elproc->SelectModelForMaterial(eth-keV, idx); } 1025 } 1026 1027 // Search for discrete process 1004 1028 if(!currentModel) { 1005 1029 const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector(); … … 1010 1034 { 1011 1035 currentModel = (vem[i])->SelectModelForMaterial(kinEnergy, idx); 1012 loweModel = (vem[i])->SelectModelForMaterial(keV, idx); 1036 G4double eth = currentModel->LowEnergyLimit(); 1037 if(eth > 0.0) { loweModel = (vem[i])->SelectModelForMaterial(eth-keV, idx); } 1013 1038 isApplicable = true; 1014 1039 break; … … 1016 1041 } 1017 1042 } 1043 1044 // Search for msc process 1018 1045 if(!currentModel) { 1019 1046 const std::vector<G4VMultipleScattering*> vmsc = … … 1025 1052 { 1026 1053 currentModel = (vmsc[i])->SelectModelForMaterial(kinEnergy, idx); 1027 loweModel = (vmsc[i])->SelectModelForMaterial(keV, idx); 1054 G4double eth = currentModel->LowEnergyLimit(); 1055 if(eth > 0.0) { loweModel = (vmsc[i])->SelectModelForMaterial(eth-keV, idx); } 1028 1056 isApplicable = true; 1029 1057 break; … … 1043 1071 G4String partname = p->GetParticleName(); 1044 1072 const G4ParticleDefinition* part = p; 1073 1045 1074 if(p->GetParticleType() == "nucleus" && 1046 1075 partname != "deuteron" && 1047 partname != "triton") part = G4GenericIon::GenericIon();1076 partname != "triton") { part = theGenericIon; } 1048 1077 1049 1078 G4LossTableManager* lManager = G4LossTableManager::Instance(); … … 1052 1081 G4int n = vel.size(); 1053 1082 for(G4int i=0; i<n; i++) { 1054 if( (vel[i])->Particle() == part) {1083 if( (vel[i])->Particle() == part ) { 1055 1084 elp = vel[i]; 1056 1085 break; -
trunk/source/processes/electromagnetic/utils/src/G4EmConfigurator.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmConfigurator.cc,v 1. 4 2009/02/26 11:33:33vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4EmConfigurator.cc,v 1.5 2009/07/20 17:06:48 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 223 223 // set fluctuation model for ionisation processes 224 224 if(fluc && ptype == eloss) { 225 G4VEnergyLossProcess* p = reinterpret_cast<G4VEnergyLossProcess*>(proc);225 G4VEnergyLossProcess* p = static_cast<G4VEnergyLossProcess*>(proc); 226 226 p->SetFluctModel(fluc); 227 227 … … 262 262 // added model 263 263 if(ptype == eloss) { 264 G4VEnergyLossProcess* p = reinterpret_cast<G4VEnergyLossProcess*>(proc);264 G4VEnergyLossProcess* p = static_cast<G4VEnergyLossProcess*>(proc); 265 265 p->AddEmModel(index,mod,fluc,reg); 266 266 //G4cout << "### Added eloss model order= " << index << " for " 267 267 // << particleName << " and " << processName << " " << mod << G4endl; 268 268 } else if(ptype == discrete) { 269 G4VEmProcess* p = reinterpret_cast<G4VEmProcess*>(proc);269 G4VEmProcess* p = static_cast<G4VEmProcess*>(proc); 270 270 p->AddEmModel(index,mod,reg); 271 271 } else if(ptype == msc) { 272 G4VMultipleScattering* p = reinterpret_cast<G4VMultipleScattering*>(proc);272 G4VMultipleScattering* p = static_cast<G4VMultipleScattering*>(proc); 273 273 p->AddEmModel(index,mod,reg); 274 274 } -
trunk/source/processes/electromagnetic/utils/src/G4EmCorrections.cc
r1007 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmCorrections.cc,v 1.5 1 2008/12/18 13:01:44 gunterExp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4EmCorrections.cc,v 1.54 2009/10/29 17:56:36 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 687 687 massFactor = proton_mass_c2/p->GetPDGMass(); 688 688 idx = -1; 689 G4int dz = 1000;690 689 691 690 for(G4int i=0; i<nIons; i++) { 692 if(materialList[i] == mat) { 693 G4int delz = currentZ - Zion[i]; 694 if(delz < 0) delz = -delz; 695 if(delz < dz) { 696 idx = i; 697 dz = delz; 698 if(0 == delz) break; 699 } 691 if(materialList[i] == mat && currentZ == Zion[i]) { 692 idx = i; 693 break; 700 694 } 701 695 } 702 // G4cout << " idx= " << idx << " dz= " << dz <<G4endl;703 if(idx > 0) {696 // G4cout << " idx= " << idx << " dz= " << G4endl; 697 if(idx >= 0) { 704 698 if(!ionList[idx]) BuildCorrectionVector(); 705 699 if(ionList[idx]) curVector = stopData[idx]; 706 } 700 } else { return factor; } 707 701 } 708 702 if(curVector) { 709 G4bool b; 710 factor = curVector->GetValue(ekin*massFactor,b); 703 factor = curVector->Value(ekin*massFactor); 711 704 if(verbose > 1) { 712 705 G4cout << "E= " << ekin << " factor= " << factor << " massfactor= " … … 771 764 << " massRatio= " << massRatio << G4endl; 772 765 } 773 G4bool b;774 766 775 767 G4PhysicsLogVector* vv = … … 777 769 vv->SetSpline(true); 778 770 G4double e, eion, dedx, dedx1; 779 G4double eth0 = v-> GetLowEdgeEnergy(0);771 G4double eth0 = v->Energy(0); 780 772 G4double escal = eth/massRatio; 781 773 G4double qe = … … 790 782 // << " dedxt1= " << dedx1t << G4endl; 791 783 792 for(G4int i=0; i< nbinCorr; i++) {793 e = vv-> GetLowEdgeEnergy(i);784 for(G4int i=0; i<=nbinCorr; i++) { 785 e = vv->Energy(i); 794 786 escal = e/massRatio; 795 787 eion = escal/A; 796 788 if(eion <= eth0) { 797 dedx = v-> GetValue(eth0, b)*std::sqrt(eion/eth0);789 dedx = v->Value(eth0)*std::sqrt(eion/eth0); 798 790 } else { 799 dedx = v-> GetValue(eion, b);791 dedx = v->Value(eion); 800 792 } 801 793 qe = effCharge.EffectiveChargeSquareRatio(curParticle,curMaterial,escal); -
trunk/source/processes/electromagnetic/utils/src/G4EmElementSelector.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmElementSelector.cc,v 1.1 0 2009/05/26 16:59:35vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4EmElementSelector.cc,v 1.11 2009/09/29 11:31:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 68 68 if(nElmMinusOne > 0) { 69 69 xSections.reserve(n); 70 for(G4int i=0; i<n; i++) {70 for(G4int i=0; i<n; ++i) { 71 71 G4PhysicsLogVector* v = new G4PhysicsLogVector(lowEnergy,highEnergy,nbins); 72 72 //v->SetSpline(spline); … … 82 82 { 83 83 if(nElmMinusOne > 0) { 84 for(G4int i=0; i<=nElmMinusOne; i++) {84 for(G4int i=0; i<=nElmMinusOne; ++i) { 85 85 delete xSections[i]; 86 86 } … … 105 105 106 106 // loop over bins 107 for(G4int j=0; j< nbins; j++) {108 G4double e = (xSections[0])-> GetLowEdgeEnergy(j);107 for(G4int j=0; j<=nbins; ++j) { 108 G4double e = (xSections[0])->Energy(j); 109 109 model->SetupForMaterial(part, material, e); 110 110 cross = 0.0; 111 111 //G4cout << "j= " << j << " e(MeV)= " << e/MeV << G4endl; 112 for (i=0; i<=nElmMinusOne; i++) {112 for (i=0; i<=nElmMinusOne; ++i) { 113 113 cross += theAtomNumDensityVector[i]* 114 114 model->ComputeCrossSectionPerAtom(part, (*theElementVector)[i], e, … … 120 120 // xSections start from null, so use probabilities from the next bin 121 121 if(DBL_MIN >= (*xSections[nElmMinusOne])[0]) { 122 for (i=0; i<=nElmMinusOne; i++) {122 for (i=0; i<=nElmMinusOne; ++i) { 123 123 xSections[i]->PutValue(0, (*xSections[i])[1]); 124 124 } 125 125 } 126 126 // xSections ends with null, so use probabilities from the previous bin 127 if(DBL_MIN >= (*xSections[nElmMinusOne])[nbins -1]) {128 for (i=0; i<=nElmMinusOne; i++) {129 xSections[i]->PutValue(nbins -1, (*xSections[i])[nbins-2]);127 if(DBL_MIN >= (*xSections[nElmMinusOne])[nbins]) { 128 for (i=0; i<=nElmMinusOne; ++i) { 129 xSections[i]->PutValue(nbins, (*xSections[i])[nbins-1]); 130 130 } 131 131 } 132 132 // perform normalization 133 for(G4int j=0; j< nbins; j++) {133 for(G4int j=0; j<=nbins; ++j) { 134 134 cross = (*xSections[nElmMinusOne])[j]; 135 135 // only for positive X-section 136 136 if(cross > DBL_MIN) { 137 for (i=0; i<nElmMinusOne; i++) { 138 xSections[i]->PutValue(j, (*xSections[i])[j]/cross); 137 for (i=0; i<nElmMinusOne; ++i) { 138 G4double x = (*xSections[i])[j]/cross; 139 xSections[i]->PutValue(j, x); 139 140 } 140 141 } … … 158 159 << (*theElementVector)[nElmMinusOne]->GetName() 159 160 << G4endl; 161 G4cout << G4endl; 160 162 } 161 163 -
trunk/source/processes/electromagnetic/utils/src/G4EmModelManager.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmModelManager.cc,v 1. 49 2009/04/17 10:35:32vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4EmModelManager.cc,v 1.58 2009/10/29 18:07:08 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 62 62 // 12-04-07 Add verbosity at destruction (V.Ivanchenko) 63 63 // 08-04-08 Fixed and simplified initialisation of G4RegionModel (VI) 64 // 03-08-09 Create internal vectors only it is needed (VI) 64 65 // 65 66 // Class Description: … … 85 86 theListOfModelIndexes = new G4int [nModelsForRegion]; 86 87 lowKineticEnergy = new G4double [nModelsForRegion+1]; 87 for (G4int i=0; i<nModelsForRegion; i++) {88 for (G4int i=0; i<nModelsForRegion; ++i) { 88 89 theListOfModelIndexes[i] = indx[i]; 89 90 lowKineticEnergy[i] = lowE[i]; … … 114 115 #include "G4Positron.hh" 115 116 #include "G4UnitsTable.hh" 117 #include "G4DataVector.hh" 116 118 117 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... … … 120 122 nEmModels(0), 121 123 nRegions(0), 122 nCouples(0),123 idxOfRegionModels(0),124 setOfRegionModels(0),125 minSubRange(0.1),126 124 particle(0), 127 125 verboseLevel(0) 128 126 { 129 models.clear();130 flucModels.clear();131 regions.clear();132 orderOfModels.clear();133 maxCutInRange = 12.*cm;134 127 maxSubCutInRange = 0.7*mm; 135 theGamma = G4Gamma::Gamma();136 thePositron = G4Positron::Positron();137 128 models.reserve(4); 138 129 flucModels.reserve(4); … … 140 131 orderOfModels.reserve(4); 141 132 isUsed.reserve(4); 133 severalModels = true; 134 currRegionModel = 0; 135 currModel = 0; 142 136 } 143 137 … … 146 140 G4EmModelManager::~G4EmModelManager() 147 141 { 142 verboseLevel = 0; // no verbosity at destruction 148 143 Clear(); 149 144 } … … 156 151 G4cout << "G4EmModelManager::Clear()" << G4endl; 157 152 } 158 159 theCuts.clear(); 160 theSubCuts.clear(); 161 if(idxOfRegionModels) delete [] idxOfRegionModels; 162 if(setOfRegionModels && nRegions) { 163 for(G4int i=0; i<nRegions; i++) { 164 delete (setOfRegionModels[i]); 165 } 166 delete [] setOfRegionModels; 167 } 168 idxOfRegionModels = 0; 169 setOfRegionModels = 0; 153 size_t n = setOfRegionModels.size(); 154 if(n > 0) { 155 for(size_t i=0; i<n; ++i) { 156 delete setOfRegionModels[i]; 157 setOfRegionModels[i] = 0; 158 } 159 } 170 160 } 171 161 … … 195 185 { 196 186 if (nEmModels > 0) { 197 for(G4int i=0; i<nEmModels; i++) {187 for(G4int i=0; i<nEmModels; ++i) { 198 188 if(nam == models[i]->GetName()) { 199 189 models[i]->SetLowEnergyLimit(emin); … … 226 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 227 217 228 const G4DataVector* G4EmModelManager::Initialise(const G4ParticleDefinition* p, 229 const G4ParticleDefinition* sp, 230 G4double theMinSubRange, 231 G4int val) 218 const G4DataVector* 219 G4EmModelManager::Initialise(const G4ParticleDefinition* p, 220 const G4ParticleDefinition* secondaryParticle, 221 G4double minSubRange, 222 G4int val) 232 223 { 233 224 verboseLevel = val; … … 238 229 } 239 230 // Are models defined? 240 if(!nEmModels) { 241 G4Exception("G4EmModelManager::Initialise without any model defined for "+partname); 242 } 231 if(nEmModels < 1) { 232 G4Exception("G4EmModelManager::Initialise: no model defined for " + partname); 233 } 234 243 235 particle = p; 244 secondaryParticle = sp; 245 minSubRange = theMinSubRange; 246 247 Clear(); 236 Clear(); // needed if run is not first 237 248 238 G4RegionStore* regionStore = G4RegionStore::GetInstance(); 249 239 const G4Region* world = … … 256 246 G4bool isWorld = false; 257 247 258 for (G4int ii=0; ii<nEmModels; ii++) {248 for (G4int ii=0; ii<nEmModels; ++ii) { 259 249 const G4Region* r = regions[ii]; 260 250 if ( r == 0 || r == world) { … … 264 254 G4bool newRegion = true; 265 255 if (nRegions>1) { 266 for (G4int j=1; j<nRegions; j++) {256 for (G4int j=1; j<nRegions; ++j) { 267 257 if ( r == setr[j] ) newRegion = false; 268 258 } … … 274 264 } 275 265 } 266 // Are models defined? 267 if(!isWorld) { 268 G4Exception("G4EmModelManager::Initialise: no models defined for " + 269 partname + " in the World volume"); 270 } 276 271 277 272 G4ProductionCutsTable* theCoupleTable= 278 273 G4ProductionCutsTable::GetProductionCutsTable(); 279 G4int numOfCouples = theCoupleTable->GetTableSize(); 280 if(nRegions > 1) idxOfRegionModels = new G4int[numOfCouples]; 281 setOfRegionModels = new G4RegionModels*[nRegions]; 274 size_t numOfCouples = theCoupleTable->GetTableSize(); 275 276 // prepare vectors, shortcut for the case of only 1 model 277 if(nRegions > 1 && nEmModels > 1) { 278 if(numOfCouples > idxOfRegionModels.size()) idxOfRegionModels.resize(numOfCouples); 279 } 280 size_t nr = 1; 281 if(nEmModels > 1) nr = nRegions; 282 if(nr > setOfRegionModels.size()) setOfRegionModels.resize(nr); 282 283 283 284 std::vector<G4int> modelAtRegion(nEmModels); … … 287 288 288 289 // Order models for regions 289 for (G4int reg=0; reg<nRegions; reg++) {290 for (G4int reg=0; reg<nRegions; ++reg) { 290 291 const G4Region* region = setr[reg]; 291 292 G4int n = 0; 292 293 293 if(isWorld || 0 < reg) { 294 295 for (G4int ii=0; ii<nEmModels; ii++) { 296 297 G4VEmModel* model = models[ii]; 298 if ( region == regions[ii] ) { 299 300 G4double tmin = model->LowEnergyLimit(); 301 G4double tmax = model->HighEnergyLimit(); 302 G4int ord = orderOfModels[ii]; 303 G4bool push = true; 304 G4bool insert = false; 305 G4int idx = n; 306 307 if(1 < verboseLevel) { 308 G4cout << "Model #" << ii 309 << " <" << model->GetName() << "> for region <"; 310 if (region) G4cout << region->GetName(); 311 G4cout << "> " 312 << " tmin(MeV)= " << tmin/MeV 313 << "; tmax(MeV)= " << tmax/MeV 314 << "; order= " << ord 315 << G4endl; 316 } 294 for (G4int ii=0; ii<nEmModels; ++ii) { 295 296 G4VEmModel* model = models[ii]; 297 if ( region == regions[ii] ) { 298 299 G4double tmin = model->LowEnergyLimit(); 300 G4double tmax = model->HighEnergyLimit(); 301 G4int ord = orderOfModels[ii]; 302 G4bool push = true; 303 G4bool insert = false; 304 G4int idx = n; 305 306 if(1 < verboseLevel) { 307 G4cout << "Model #" << ii 308 << " <" << model->GetName() << "> for region <"; 309 if (region) G4cout << region->GetName(); 310 G4cout << "> " 311 << " tmin(MeV)= " << tmin/MeV 312 << "; tmax(MeV)= " << tmax/MeV 313 << "; order= " << ord 314 << G4endl; 315 } 317 316 318 if(n > 0) { 319 320 // extend energy range to previous models 321 tmin = std::min(tmin, eHigh[n-1]); 322 tmax = std::max(tmax, eLow[0]); 323 //G4cout << "tmin= " << tmin << " tmax= " 324 // << tmax << " ord= " << ord <<G4endl; 325 // empty energy range 326 if( tmax - tmin <= eV) push = false; 327 // low-energy model 328 else if (tmax == eLow[0]) { 329 push = false; 330 insert = true; 331 idx = 0; 332 // resolve intersections 333 } else if(tmin < eHigh[n-1]) { 334 // compare order 335 for(G4int k=0; k<n; k++) { 336 // new model has lower application 337 if(ord >= modelOrd[k]) { 338 if(tmin < eHigh[k] && tmin >= eLow[k]) tmin = eHigh[k]; 339 if(tmax <= eHigh[k] && tmax > eLow[k]) tmax = eLow[k]; 340 if(tmax > eHigh[k] && tmin < eLow[k]) { 341 if(tmax - eHigh[k] > eLow[k] - tmin) tmin = eHigh[k]; 342 else tmax = eLow[k]; 343 } 344 if( tmax - tmin <= eV) { 345 push = false; 346 break; 347 } 317 if(n > 0) { 318 319 // extend energy range to previous models 320 tmin = std::min(tmin, eHigh[n-1]); 321 tmax = std::max(tmax, eLow[0]); 322 //G4cout << "tmin= " << tmin << " tmax= " 323 // << tmax << " ord= " << ord <<G4endl; 324 // empty energy range 325 if( tmax - tmin <= eV) push = false; 326 // low-energy model 327 else if (tmax == eLow[0]) { 328 push = false; 329 insert = true; 330 idx = 0; 331 // resolve intersections 332 } else if(tmin < eHigh[n-1]) { 333 // compare order 334 for(G4int k=0; k<n; ++k) { 335 // new model has lower application 336 if(ord >= modelOrd[k]) { 337 if(tmin < eHigh[k] && tmin >= eLow[k]) tmin = eHigh[k]; 338 if(tmax <= eHigh[k] && tmax > eLow[k]) tmax = eLow[k]; 339 if(tmax > eHigh[k] && tmin < eLow[k]) { 340 if(tmax - eHigh[k] > eLow[k] - tmin) tmin = eHigh[k]; 341 else tmax = eLow[k]; 342 } 343 if( tmax - tmin <= eV) { 344 push = false; 345 break; 348 346 } 349 347 } 350 //G4cout << "tmin= " << tmin << " tmax= " 351 // << tmax << " push= " << push << " idx= " << idx <<G4endl; 352 if(push) { 353 if (tmax == eLow[0]) { 348 } 349 //G4cout << "tmin= " << tmin << " tmax= " 350 // << tmax << " push= " << push << " idx= " << idx <<G4endl; 351 if(push) { 352 if (tmax == eLow[0]) { 353 push = false; 354 insert = true; 355 idx = 0; 356 // continue resolve intersections 357 } else if(tmin < eHigh[n-1]) { 358 // last energy interval 359 if(tmin > eLow[n-1] && tmax >= eHigh[n-1]) { 360 eHigh[n-1] = tmin; 361 // first energy interval 362 } else if(tmin <= eLow[0] && tmax < eHigh[0]) { 363 eLow[0] = tmax; 354 364 push = false; 355 365 insert = true; 356 366 idx = 0; 357 // continue resolve intersections 358 } else if(tmin < eHigh[n-1]) { 359 // last energy interval 360 if(tmin > eLow[n-1] && tmax >= eHigh[n-1]) { 361 eHigh[n-1] = tmin; 362 // first energy interval 363 } else if(tmin <= eLow[0] && tmax < eHigh[0]) { 364 eLow[0] = tmax; 365 push = false; 366 insert = true; 367 idx = 0; 368 } else { 369 // find energy interval to replace 370 for(G4int k=0; k<n; k++) { 371 if(tmin <= eLow[k] && tmax >= eHigh[k]) { 372 push = false; 373 modelAtRegion[k] = ii; 374 modelOrd[k] = ord; 375 isUsed[ii] = 1; 376 } 377 } 367 } else { 368 // find energy interval to replace 369 for(G4int k=0; k<n; ++k) { 370 if(tmin <= eLow[k] && tmax >= eHigh[k]) { 371 push = false; 372 modelAtRegion[k] = ii; 373 modelOrd[k] = ord; 374 isUsed[ii] = 1; 375 } 378 376 } 379 377 } … … 381 379 } 382 380 } 383 if(insert) { 384 for(G4int k=n-1; k>=idx; k--) { 385 modelAtRegion[k+1] = modelAtRegion[k]; 386 modelOrd[k+1] = modelOrd[k]; 387 eLow[k+1] = eLow[k]; 388 eHigh[k+1] = eHigh[k]; 389 } 390 } 391 //G4cout << "push= " << push << " insert= " << insert 392 //<< " idx= " << idx <<G4endl; 393 if (push || insert) { 394 n++; 395 modelAtRegion[idx] = ii; 396 modelOrd[idx] = ord; 397 eLow[idx] = tmin; 398 eHigh[idx] = tmax; 399 isUsed[ii] = 1; 381 } 382 if(insert) { 383 for(G4int k=n-1; k>=idx; --k) { 384 modelAtRegion[k+1] = modelAtRegion[k]; 385 modelOrd[k+1] = modelOrd[k]; 386 eLow[k+1] = eLow[k]; 387 eHigh[k+1] = eHigh[k]; 400 388 } 401 389 } 402 } 403 } else { 404 n = 1; 405 models.push_back(0); 406 modelAtRegion.push_back(nEmModels); 407 eLow.push_back(0.0); 408 eHigh.push_back(DBL_MAX); 390 //G4cout << "push= " << push << " insert= " << insert 391 //<< " idx= " << idx <<G4endl; 392 if (push || insert) { 393 ++n; 394 modelAtRegion[idx] = ii; 395 modelOrd[idx] = ord; 396 eLow[idx] = tmin; 397 eHigh[idx] = tmax; 398 isUsed[ii] = 1; 399 } 400 } 409 401 } 410 402 eLow[0] = 0.0; … … 415 407 if (region) G4cout << region->GetName(); 416 408 G4cout << "> Elow(MeV)= "; 417 for(G4int ii=0; ii<=n; ii++) {G4cout << eLow[ii]/MeV << " ";}409 for(G4int ii=0; ii<=n; ++ii) {G4cout << eLow[ii]/MeV << " ";} 418 410 G4cout << G4endl; 419 411 } 420 412 G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow, region); 421 413 setOfRegionModels[reg] = rm; 414 if(1 == nEmModels) break; 422 415 } 423 416 … … 425 418 426 419 // Access to materials and build cuts 427 for(G4int i=0; i<numOfCouples; i++) { 420 size_t idx = 1; 421 if(secondaryParticle) { 422 if( secondaryParticle == G4Gamma::Gamma() ) idx = 0; 423 else if( secondaryParticle == G4Positron::Positron()) idx = 2; 424 } 425 426 if(numOfCouples > theCuts.size()) {theCuts.resize(numOfCouples);} 427 if(minSubRange < 1.0 && numOfCouples > theSubCuts.size()) { 428 theSubCuts.resize(numOfCouples); 429 } 430 for(size_t i=0; i<numOfCouples; ++i) { 428 431 429 432 const G4MaterialCutsCouple* couple = … … 433 436 434 437 G4int reg = 0; 435 if(nRegions > 1 ) {438 if(nRegions > 1 && nEmModels > 1) { 436 439 reg = nRegions; 437 do { reg--;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));440 do {--reg;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts())); 438 441 idxOfRegionModels[i] = reg; 439 442 } … … 446 449 } 447 450 448 G4double cut = DBL_MAX;451 G4double cut = (*theCoupleTable->GetEnergyCutsVector(idx))[i]; 449 452 G4double subcut = DBL_MAX; 450 453 if(secondaryParticle) { 451 size_t idx = 1;452 if( secondaryParticle == theGamma ) idx = 0;453 cut = (*theCoupleTable->GetEnergyCutsVector(idx))[i];454 if( secondaryParticle == thePositron && cut < DBL_MAX )455 cut += (*theCoupleTable->GetEnergyCutsVector(2))[i] +456 2.0*electron_mass_c2;457 454 458 455 // compute subcut 459 if( cut < DBL_MAX ) subcut = minSubRange*cut; 460 if(pcuts->GetProductionCut(idx) < maxCutInRange) { 456 if( cut < DBL_MAX && minSubRange < 1.0) { 457 subcut = minSubRange*cut; 458 G4double rcut = std::min(minSubRange*pcuts->GetProductionCut(idx), 459 maxSubCutInRange); 461 460 G4double tcutmax = 462 theCoupleTable->ConvertRangeToEnergy(secondaryParticle, 463 material,maxSubCutInRange); 461 theCoupleTable->ConvertRangeToEnergy(secondaryParticle,material,rcut); 464 462 if(tcutmax < subcut) subcut = tcutmax; 465 463 } … … 467 465 468 466 G4int nm = setOfRegionModels[reg]->NumberOfModels(); 469 for(G4int j=0; j<nm; j++) {467 for(G4int j=0; j<nm; ++j) { 470 468 471 469 G4VEmModel* model = models[setOfRegionModels[reg]->ModelIndex(j)]; … … 484 482 } 485 483 } 486 theCuts.push_back(cut); 487 theSubCuts.push_back(subcut); 488 } 489 490 for(G4int jj=0; jj<nEmModels; jj++) { 484 theCuts[i] = cut; 485 if(minSubRange < 1.0) theSubCuts[i] = subcut; 486 } 487 488 // initialize models 489 G4int nn = 0; 490 severalModels = true; 491 for(G4int jj=0; jj<nEmModels; ++jj) { 491 492 if(1 == isUsed[jj]) { 492 models[jj]->Initialise(particle, theCuts); 493 ++nn; 494 currModel = models[jj]; 495 currModel->Initialise(particle, theCuts); 493 496 if(flucModels[jj]) flucModels[jj]->InitialiseMe(particle); 494 497 } 495 498 } 499 if(1 == nn) severalModels = false; 496 500 497 501 if(1 < verboseLevel) { … … 510 514 G4EmTableType tType) 511 515 { 512 513 G4double e;514 515 516 size_t i = couple->GetIndex(); 516 517 G4double cut = theCuts[i]; 517 G4double subcut= 0.0;518 G4double emin = 0.0; 518 519 519 520 if(fTotal == tType) cut = DBL_MAX; 520 else if(fSubRestricted == tType) subcut = theSubCuts[i]; 521 else if(fSubRestricted == tType) { 522 emin = cut; 523 if(theSubCuts.size() > 0) emin = theSubCuts[i]; 524 } 521 525 522 526 if(1 < verboseLevel) { … … 524 528 << couple->GetMaterial()->GetName() 525 529 << " cut(MeV)= " << cut 526 << " subcut(MeV)= " << subcut530 << " emin(MeV)= " << emin 527 531 << " Type " << tType 528 532 << " for " << particle->GetParticleName() … … 531 535 532 536 G4int reg = 0; 533 if(nRegions > 1 ) reg = idxOfRegionModels[i];537 if(nRegions > 1 && nEmModels > 1) reg = idxOfRegionModels[i]; 534 538 const G4RegionModels* regModels = setOfRegionModels[reg]; 535 539 G4int nmod = regModels->NumberOfModels(); 536 540 537 // vectors to provide continues dE/dx538 G4DataVector factor(nmod);539 G4DataVector eLow(nmod+1);540 G4DataVector dedxLow(nmod);541 G4DataVector dedxHigh(nmod);542 543 if(1 < verboseLevel) {544 G4cout << "There are " << nmod << " models for "545 << couple->GetMaterial()->GetName()546 << " at the region #" << reg547 << G4endl;548 }549 550 // calculate factors to provide continuity of energy loss551 552 553 factor[0] = 1.0;554 G4int j;555 556 G4int totBinsLoss = aVector->GetVectorLength();557 558 dedxLow[0] = 0.0;559 eLow[0] = 0.0;560 561 e = regModels->LowEdgeEnergy(1);562 eLow[1] = e;563 G4VEmModel* model = models[regModels->ModelIndex(0)];564 dedxHigh[0] = 0.0;565 if(model && cut > subcut) {566 dedxHigh[0] = model->ComputeDEDX(couple,particle,e,cut);567 if(subcut > 0.0) {568 dedxHigh[0] -= model->ComputeDEDX(couple,particle,e,subcut);569 }570 }571 if(nmod > 1) {572 for(j=1; j<nmod; j++) {573 574 e = regModels->LowEdgeEnergy(j);575 eLow[j] = e;576 G4int idx = regModels->ModelIndex(j);577 578 dedxLow[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);579 if(subcut > 0.0) {580 dedxLow[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);581 }582 if(subcut == cut) dedxLow[j] = 0.0;583 584 e = regModels->LowEdgeEnergy(j+1);585 eLow[j+1] = e;586 dedxHigh[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);587 if(subcut > 0.0) {588 dedxHigh[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);589 }590 if(subcut == cut) dedxHigh[j] = 0.0;591 }592 if(1 < verboseLevel) {593 G4cout << " model #0"594 << " dedx(" << eLow[0] << ")= " << dedxLow[0]595 << " dedx(" << eLow[1] << ")= " << dedxHigh[0]596 << G4endl;597 }598 599 for(j=1; j<nmod; j++) {600 if(dedxLow[j] > 0.0) {601 factor[j] = (dedxHigh[j-1]/dedxLow[j] - 1.0)*eLow[j];602 } else factor[j] = 0.0;603 if(1 < verboseLevel) {604 G4cout << " model #" << j605 << " dedx(" << eLow[j] << ")= " << dedxLow[j]606 << " dedx(" << eLow[j+1] << ")= " << dedxHigh[j]607 << " factor= " << factor[j]/eLow[j]608 << G4endl;609 }610 }611 612 if(2 < verboseLevel) {613 G4cout << "Loop over " << totBinsLoss << " bins start " << G4endl;614 }615 }616 617 541 // Calculate energy losses vector 618 for(j=0; j<totBinsLoss; j++) { 619 620 G4double e = aVector->GetLowEdgeEnergy(j); 621 G4double fac = 1.0; 622 623 // Choose a model of energy losses 542 543 //G4cout << "nmod= " << nmod << G4endl; 544 size_t totBinsLoss = aVector->GetVectorLength(); 545 for(size_t j=0; j<totBinsLoss; ++j) { 546 547 G4double e = aVector->Energy(j); 548 G4double del = 0.0; 549 550 // Choose a model of energy losses 624 551 G4int k = 0; 625 if (nmod > 1 && e > eLow[1]) { 626 do { 627 k++; 628 fac *= (1.0 + factor[k]/e); 629 } while (k+1 < nmod && e > eLow[k+1]); 630 } 631 632 model = models[regModels->ModelIndex(k)]; 633 G4double dedx = 0.0; 634 G4double dedx0 = 0.0; 635 if(model && cut > subcut) { 636 dedx = model->ComputeDEDX(couple,particle,e,cut); 637 dedx0 = dedx; 638 if(subcut > 0.0) dedx -= model->ComputeDEDX(couple,particle,e,subcut); 639 dedx *= fac; 640 } 552 if (nmod > 1) { 553 k = nmod; 554 do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k)); 555 //G4cout << "k= " << k << G4endl; 556 if(k > 0) { 557 G4double elow = regModels->LowEdgeEnergy(k); 558 G4double dedx1 = ComputeDEDX(models[regModels->ModelIndex(k-1)], 559 couple,elow,cut,emin); 560 G4double dedx2 = ComputeDEDX(models[regModels->ModelIndex(k)], 561 couple,elow,cut,emin); 562 del = (dedx1 - dedx2)*elow/e; 563 //G4cout << "elow= " << elow 564 // << " dedx1= " << dedx1 << " dedx2= " << dedx2 << G4endl; 565 } 566 } 567 G4double dedx = ComputeDEDX(models[regModels->ModelIndex(k)], 568 couple,e,cut,emin) + del; 641 569 642 570 if(dedx < 0.0) dedx = 0.0; 643 571 if(2 < verboseLevel) { 644 645 646 647 << " dEdx0(MeV/mm)= " << dedx0*mm/MeV 648 << " fac= " << fac 649 572 G4cout << "Material= " << couple->GetMaterial()->GetName() 573 << " E(MeV)= " << e/MeV 574 << " dEdx(MeV/mm)= " << dedx*mm/MeV 575 << " del= " << del*mm/MeV<< " k= " << k 576 << " modelIdx= " << regModels->ModelIndex(k) 577 << G4endl; 650 578 } 651 579 aVector->PutValue(j, dedx); … … 660 588 G4EmTableType tType) 661 589 { 662 G4double e;663 664 590 size_t i = couple->GetIndex(); 665 591 G4double cut = theCuts[i]; … … 667 593 if (fSubRestricted == tType) { 668 594 tmax = cut; 669 cut = theSubCuts[i]; 670 } 671 595 if(theSubCuts.size() > 0) cut = theSubCuts[i]; 596 } 597 598 G4int reg = 0; 599 if(nRegions > 1 && nEmModels > 1) reg = idxOfRegionModels[i]; 600 const G4RegionModels* regModels = setOfRegionModels[reg]; 601 G4int nmod = regModels->NumberOfModels(); 672 602 if(1 < verboseLevel) { 673 G4cout << "G4EmModelManager::FillLambdaVector() for particle"603 G4cout << "G4EmModelManager::FillLambdaVector() for " 674 604 << particle->GetParticleName() 675 605 << " in " << couple->GetMaterial()->GetName() 676 << " Ecut(MeV)= " << cut 677 << " Emax(MeV)= " << tmax 678 << " Type " << tType 606 << " Ecut(MeV)= " << cut 607 << " Emax(MeV)= " << tmax 608 << " Type " << tType 609 << " nmod= " << nmod 679 610 << G4endl; 680 611 } 681 612 682 G4int reg = 0;683 if(nRegions > 1) reg = idxOfRegionModels[i];684 const G4RegionModels* regModels = setOfRegionModels[reg];685 G4int nmod = regModels->NumberOfModels();686 687 // vectors to provide continues dE/dx688 G4DataVector factor(nmod);689 G4DataVector eLow(nmod+1);690 G4DataVector sigmaLow(nmod);691 G4DataVector sigmaHigh(nmod);692 693 if(2 < verboseLevel) {694 G4cout << "There are " << nmod << " models for "695 << couple->GetMaterial()->GetName() << G4endl;696 }697 698 // calculate factors to provide continuity of energy loss699 factor[0] = 1.0;700 G4int j;701 G4int totBinsLambda = aVector->GetVectorLength();702 703 sigmaLow[0] = 0.0;704 eLow[0] = 0.0;705 706 e = regModels->LowEdgeEnergy(1);707 eLow[1] = e;708 G4VEmModel* model = models[regModels->ModelIndex(0)];709 sigmaHigh[0] = 0.0;710 if(model) sigmaHigh[0] = model->CrossSection(couple,particle,e,cut,tmax);711 712 if(2 < verboseLevel) {713 G4cout << "### For material " << couple->GetMaterial()->GetName()714 << " " << nmod715 << " models"716 << " Ecut(MeV)= " << cut/MeV717 << " Emax(MeV)= " << e/MeV718 << " nbins= " << totBinsLambda719 << G4endl;720 G4cout << " model #0 eUp= " << e721 << " sigmaUp= " << sigmaHigh[0] << G4endl;722 }723 724 725 if(nmod > 1) {726 727 for(j=1; j<nmod; j++) {728 729 e = regModels->LowEdgeEnergy(j);730 eLow[j] = e;731 G4int idx = regModels->ModelIndex(j);732 733 sigmaLow[j] = models[idx]->CrossSection(couple,particle,e,cut,tmax);734 e = regModels->LowEdgeEnergy(j+1);735 eLow[j+1] = e;736 sigmaHigh[j] = models[idx]->CrossSection(couple,particle,e,cut,tmax);737 }738 if(1 < verboseLevel) {739 G4cout << " model #0"740 << " sigma(" << eLow[0] << ")= " << sigmaLow[0]741 << " sigma(" << eLow[1] << ")= " << sigmaHigh[0]742 << G4endl;743 }744 for(j=1; j<nmod; j++) {745 if(sigmaLow[j] > 0.0) {746 factor[j] = (sigmaHigh[j-1]/sigmaLow[j] - 1.0)*eLow[j];747 } else factor[j] = 0.0;748 if(1 < verboseLevel) {749 G4cout << " model #" << j750 << " sigma(" << eLow[j] << ")= " << sigmaLow[j]751 << " sigma(" << eLow[j+1] << ")= " << sigmaHigh[j]752 << " factor= " << factor[j]/eLow[j]753 << G4endl;754 }755 }756 }757 613 758 614 // Calculate lambda vector 759 for(j=0; j<totBinsLambda; j++) { 760 761 e = aVector->GetLowEdgeEnergy(j); 762 763 // Choose a model of energy losses 615 size_t totBinsLambda = aVector->GetVectorLength(); 616 for(size_t j=0; j<totBinsLambda; ++j) { 617 618 G4double e = aVector->Energy(j); 619 620 G4double del = 0.0; 621 622 // Choose a model 764 623 G4int k = 0; 765 G4double fac = 1.0; 766 if (nmod > 1 && e > eLow[1]) { 767 do { 768 k++; 769 fac *= (1.0 + factor[k]/e); 770 } while ( k+1 < nmod && e > eLow[k+1] ); 771 } 772 773 model = models[regModels->ModelIndex(k)]; 774 G4double cross = 0.0; 775 if(model) cross = model->CrossSection(couple,particle,e,cut,tmax)*fac; 624 G4VEmModel* mod = models[regModels->ModelIndex(0)]; 625 if (nmod > 1) { 626 k = nmod; 627 do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k)); 628 if(k > 0) { 629 G4double elow = regModels->LowEdgeEnergy(k); 630 G4VEmModel* m = models[regModels->ModelIndex(k-1)]; 631 G4double xs1 = m->CrossSection(couple,particle,elow,cut,tmax); 632 mod = models[regModels->ModelIndex(k)]; 633 G4double xs2 = mod->CrossSection(couple,particle,elow,cut,tmax); 634 del = (xs1 - xs2)*elow/e; 635 } 636 } 637 G4double cross = mod->CrossSection(couple,particle,e,cut,tmax) + del; 638 776 639 if(j==0 && startFromNull) cross = 0.0; 777 640 … … 779 642 G4cout << "FillLambdaVector: " << j << ". e(MeV)= " << e/MeV 780 643 << " cross(1/mm)= " << cross*mm 781 << " fac= " << fac<< " k= " << k782 << " model = " << regModels->ModelIndex(k)644 << " del= " << del*mm << " k= " << k 645 << " modelIdx= " << regModels->ModelIndex(k) 783 646 << G4endl; 784 647 } … … 794 657 { 795 658 if(verb == 0) return; 796 for(G4int i=0; i<nRegions; i++) {659 for(G4int i=0; i<nRegions; ++i) { 797 660 G4RegionModels* r = setOfRegionModels[i]; 798 661 const G4Region* reg = r->Region(); 799 // if(verb > -1 || nRegions > 1) {800 // }801 662 G4int n = r->NumberOfModels(); 802 if( verb > -1 ||n > 0) {663 if(n > 0) { 803 664 G4cout << " ===== EM models for the G4Region " << reg->GetName() 804 665 << " ======" << G4endl;; 805 for(G4int j=0; j<n; j++) {666 for(G4int j=0; j<n; ++j) { 806 667 const G4VEmModel* m = models[r->ModelIndex(j)]; 807 668 G4cout << std::setw(20); … … 813 674 } 814 675 } 815 } 816 } 817 818 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 676 if(1 == nEmModels) break; 677 } 678 } 679 680 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... -
trunk/source/processes/electromagnetic/utils/src/G4EmMultiModel.cc
r1007 r1196 25 25 // 26 26 // $Id: G4EmMultiModel.cc,v 1.6 2007/05/22 17:31:58 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/src/G4EmProcessOptions.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmProcessOptions.cc,v 1.2 6 2009/02/18 14:43:27vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4EmProcessOptions.cc,v 1.27 2009/10/29 19:25:28 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 516 516 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 517 517 518 void G4EmProcessOptions::SetFactorForAngleLimit(G4double val) 519 { 520 theManager->SetFactorForAngleLimit(val); 521 } 522 523 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 524 518 525 void G4EmProcessOptions::SetLPMFlag(G4bool val) 519 526 { -
trunk/source/processes/electromagnetic/utils/src/G4EmSaturation.cc
r1007 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4EmSaturation.cc,v 1. 9 2008/11/12 15:37:33vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4EmSaturation.cc,v 1.10 2009/09/25 09:16:40 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 46 46 47 47 #include "G4EmSaturation.hh" 48 #include "G4Gamma.hh"49 #include "G4Electron.hh"50 #include "G4Neutron.hh"51 #include "G4Proton.hh"52 48 #include "G4LossTableManager.hh" 53 49 #include "G4NistManager.hh" 54 50 #include "G4Material.hh" 55 51 #include "G4MaterialCutsCouple.hh" 52 #include "G4Electron.hh" 53 #include "G4Proton.hh" 56 54 57 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... … … 66 64 curChargeSq = 1.0; 67 65 nMaterials = 0; 68 Initialise(); 66 electron = 0; 67 Initialise(); 69 68 } 70 69 … … 90 89 if(bfactor > 0.0) { 91 90 92 // atomic relaxations 93 if(p == gamma) { 91 G4int pdgCode = p->GetPDGEncoding(); 92 // atomic relaxations for gamma incident 93 if(22 == pdgCode) { 94 94 evis /= (1.0 + bfactor*edep/manager->GetRange(electron,edep,couple)); 95 95 … … 101 101 if(nloss < 0.0) nloss = 0.0; 102 102 G4double eloss = edep - nloss; 103 if(p == neutron || eloss < 0.0 || length <= 0.0) { 103 104 // neutrons 105 if(2112 == pdgCode || eloss < 0.0 || length <= 0.0) { 104 106 nloss = edep; 105 107 eloss = 0.0; … … 111 113 // non-ionizing energy loss 112 114 if(nloss > 0.0) { 115 if(!proton) {proton = G4Proton::Proton();} 113 116 G4double escaled = nloss*curRatio; 114 117 G4double s = manager->GetRange(proton,escaled,couple)/curChargeSq; … … 146 149 G4double G4EmSaturation::FindBirksCoefficient(const G4Material* mat) 147 150 { 151 // electron should exist in any case 152 if(!manager) { 153 manager = G4LossTableManager::Instance(); 154 nist = G4NistManager::Instance(); 155 electron= G4Electron::Electron(); 156 proton = 0; 157 } 158 148 159 if(mat == curMaterial) return curBirks; 149 160 … … 161 172 return curBirks; 162 173 } 163 }164 165 if(!manager) {166 manager = G4LossTableManager::Instance();167 nist = G4NistManager::Instance();168 gamma = G4Gamma::Gamma();169 electron= G4Electron::Electron();170 proton = G4Proton::Proton();171 neutron = G4Neutron::Neutron();172 174 } 173 175 -
trunk/source/processes/electromagnetic/utils/src/G4EnergyLossMessenger.cc
r1055 r1196 25 25 // 26 26 // 27 // $Id: G4EnergyLossMessenger.cc,v 1.3 7 2009/02/18 14:43:27vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // $Id: G4EnergyLossMessenger.cc,v 1.38 2009/10/29 19:25:28 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 196 196 dedxCmd->SetGuidance("Set number of bins for DEDX tables"); 197 197 dedxCmd->SetParameterName("binsDEDX",true); 198 dedxCmd->SetDefaultValue( 120);198 dedxCmd->SetDefaultValue(77); 199 199 dedxCmd->AvailableForStates(G4State_PreInit,G4State_Idle); 200 200 … … 202 202 lamCmd->SetGuidance("Set number of bins for Lambda tables"); 203 203 lamCmd->SetParameterName("binsL",true); 204 lamCmd->SetDefaultValue( 120);204 lamCmd->SetDefaultValue(77); 205 205 lamCmd->AvailableForStates(G4State_PreInit,G4State_Idle); 206 206 … … 242 242 frCmd->SetParameterName("Fr",true); 243 243 frCmd->SetRange("Fr>0"); 244 frCmd->SetDefaultValue(0.0 2);244 frCmd->SetDefaultValue(0.04); 245 245 frCmd->AvailableForStates(G4State_PreInit,G4State_Idle); 246 246 … … 252 252 fgCmd->AvailableForStates(G4State_PreInit,G4State_Idle); 253 253 254 mscfCmd = new G4UIcmdWithADouble("/process/msc/FactorForAngleLimit",this); 255 mscfCmd->SetGuidance("Set factor for computation of a limit for -t (invariant trasfer)"); 256 mscfCmd->SetParameterName("Fact",true); 257 mscfCmd->SetRange("Fact>0"); 258 mscfCmd->SetDefaultValue(1.); 259 mscfCmd->AvailableForStates(G4State_PreInit,G4State_Idle); 260 254 261 skinCmd = new G4UIcmdWithADouble("/process/msc/Skin",this); 255 262 skinCmd->SetGuidance("Set skin parameter for msc processes"); … … 258 265 259 266 angCmd = new G4UIcmdWithADoubleAndUnit("/process/msc/ThetaLimit",this); 260 angCmd->SetGuidance("Set the limit on the polar angle ");267 angCmd->SetGuidance("Set the limit on the polar angle for msc and single scattering"); 261 268 angCmd->SetParameterName("theta",true); 262 269 angCmd->SetUnitCategory("Angle"); … … 297 304 delete skinCmd; 298 305 delete angCmd; 306 delete mscfCmd; 299 307 } 300 308 … … 432 440 G4UImanager::GetUIpointer()->ApplyCommand("/run/physicsModified"); 433 441 } 442 if (command == mscfCmd) { 443 opt->SetFactorForAngleLimit(mscfCmd->GetNewDoubleValue(newValue)); 444 G4UImanager::GetUIpointer()->ApplyCommand("/run/physicsModified"); 445 } 434 446 if (command == angCmd) { 435 447 opt->SetPolarAngleLimit(angCmd->GetNewDoubleValue(newValue)); -
trunk/source/processes/electromagnetic/utils/src/G4EnergyLossTables.cc
r1007 r1196 26 26 // 27 27 // $Id: G4EnergyLossTables.cc,v 1.34 2008/07/08 10:57:22 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/src/G4LossTableBuilder.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4LossTableBuilder.cc,v 1. 28 2009/02/18 16:24:47vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4LossTableBuilder.cc,v 1.32 2009/08/11 17:24:53 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 47 47 // 16-01-07 Fill new (not old) DEDX table (V.Ivanchenko) 48 48 // 12-02-07 Use G4LPhysicsFreeVector for the inverse range table (V.Ivanchenko) 49 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko) 49 50 // 50 51 // Class Description: … … 85 86 if(0 >= n_vectors) return; 86 87 87 G4bool b; 88 89 G4PhysicsVector* pv = (*(list[0]))[0]; 90 size_t nbins = pv->GetVectorLength(); 91 G4double elow = pv->GetLowEdgeEnergy(0); 92 G4double ehigh = pv->GetLowEdgeEnergy(nbins); 88 G4PhysicsLogVector* pv0 = static_cast<G4PhysicsLogVector*>((*(list[0]))[0]); 89 size_t npoints = pv0->GetVectorLength(); 93 90 for (size_t i=0; i<n_vectors; i++) { 94 91 95 pv = new G4PhysicsLogVector(elow, ehigh, nbins); 92 G4PhysicsLogVector* pv = new G4PhysicsLogVector(*pv0); 93 // pv = new G4PhysicsLogVector(elow, ehigh, npoints-1); 96 94 pv->SetSpline(splineFlag); 97 for (size_t j=0; j<n bins; j++) {95 for (size_t j=0; j<npoints; j++) { 98 96 G4double dedx = 0.0; 99 G4double energy = pv->GetLowEdgeEnergy(j);100 101 97 for (size_t k=0; k<n_processes; k++) { 102 dedx += ((*(list[k]))[i])->GetValue(energy, b); 98 G4PhysicsVector* pv1 = (*(list[k]))[i]; 99 dedx += (*pv1)[j]; 103 100 } 104 101 pv->PutValue(j, dedx); 105 G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv);106 102 } 103 if(splineFlag) pv->FillSecondDerivatives(); 104 G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv); 107 105 } 108 106 } … … 118 116 if(!n_vectors) return; 119 117 120 G4bool b;121 118 size_t n = 100; 122 119 G4double del = 1.0/(G4double)n; … … 125 122 126 123 if (rangeTable->GetFlag(i) || !isIonisation) { 127 G4PhysicsVector* pv = (*dedxTable)[i]; 128 size_t nbins = pv->GetVectorLength(); 129 size_t bin0 = 0; 130 G4double elow = pv->GetLowEdgeEnergy(0); 131 G4double ehigh = pv->GetLowEdgeEnergy(nbins); 132 G4double dedx1 = pv->GetValue(elow, b); 124 G4PhysicsLogVector* pv = 125 static_cast<G4PhysicsLogVector*>((*dedxTable)[i]); 126 size_t npoints = pv->GetVectorLength(); 127 size_t bin0 = 0; 128 G4double elow = pv->Energy(0); 129 G4double ehigh = pv->Energy(npoints-1); 130 G4double dedx1 = pv->Value(elow); 133 131 134 132 //G4cout << "nbins= " << nbins << " dedx1= " << dedx1 << G4endl; … … 136 134 // protection for specific cases dedx=0 137 135 if(dedx1 == 0.0) { 138 for (size_t k=1; k<n bins-1; k++) {136 for (size_t k=1; k<npoints; k++) { 139 137 bin0++; 140 elow = pv-> GetLowEdgeEnergy(k);141 dedx1 = pv->GetValue(elow, b);138 elow = pv->Energy(k); 139 dedx1 = (*pv)[k]; 142 140 if(dedx1 > 0.0) break; 143 141 } 144 n bins -= bin0;145 } 146 147 // G4cout << "nbins= " << nbins << " elow= " << elow << " ehigh= " << ehigh <<G4endl;142 npoints -= bin0; 143 } 144 // G4cout<<"New Range vector" << G4endl; 145 // G4cout<<"nbins= "<<npoints-1<<" elow= "<<elow<<" ehigh= "<<ehigh<<G4endl; 148 146 // initialisation of a new vector 149 G4PhysicsLogVector* v = new G4PhysicsLogVector(elow, ehigh, nbins); 147 if(npoints < 2) npoints = 2; 148 G4PhysicsLogVector* v; 149 if(0 == bin0) { v = new G4PhysicsLogVector(*pv); } 150 else { v = new G4PhysicsLogVector(elow, ehigh, npoints-1); } 150 151 // dedx is exect zero 151 if( nbins <= 2) {152 if(2 == npoints) { 152 153 v->PutValue(0,1000.); 153 154 v->PutValue(1,2000.); … … 162 163 G4double energy1 = elow; 163 164 164 for (size_t j=1; j<n bins; j++) {165 166 G4double energy2 = pv-> GetLowEdgeEnergy(j+bin0);165 for (size_t j=1; j<npoints; j++) { 166 167 G4double energy2 = pv->Energy(j+bin0); 167 168 G4double de = (energy2 - energy1) * del; 168 169 G4double energy = energy2 + de*0.5; 170 G4double sum = 0.0; 169 171 for (size_t k=0; k<n; k++) { 170 172 energy -= de; 171 dedx1 = pv-> GetValue(energy, b);172 if(dedx1 > 0.0) range+= de/dedx1;173 dedx1 = pv->Value(energy); 174 if(dedx1 > 0.0) sum += de/dedx1; 173 175 } 174 176 range += sum; 175 177 // G4cout << "Range i= " <<i << " j= " << j << G4endl; 176 178 v->PutValue(j,range); 177 179 energy1 = energy2; 178 180 } 181 if(splineFlag) v->FillSecondDerivatives(); 179 182 G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v); 180 183 } … … 191 194 size_t n_vectors = rangeTable->length(); 192 195 if(!n_vectors) return; 193 G4bool b;194 196 195 197 for (size_t i=0; i<n_vectors; i++) { … … 197 199 if (invRangeTable->GetFlag(i) || !isIonisation) { 198 200 G4PhysicsVector* pv = (*rangeTable)[i]; 199 size_t nbins = pv->GetVectorLength(); 200 G4double elow = pv->GetLowEdgeEnergy(0); 201 G4double ehigh = pv->GetLowEdgeEnergy(nbins-1); 202 G4double rlow = pv->GetValue(elow, b); 203 G4double rhigh = pv->GetValue(ehigh, b); 201 size_t npoints = pv->GetVectorLength(); 202 G4double rlow = (*pv)[0]; 203 G4double rhigh = (*pv)[npoints-1]; 204 204 205 G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(n bins,rlow,rhigh);205 G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(npoints,rlow,rhigh); 206 206 v->SetSpline(splineFlag); 207 207 208 for (size_t j=0; j<n bins; j++) {209 G4double e = pv-> GetLowEdgeEnergy(j);210 G4double r = pv->GetValue(e, b);208 for (size_t j=0; j<npoints; j++) { 209 G4double e = pv->Energy(j); 210 G4double r = (*pv)[j]; 211 211 v->PutValues(j,r,e); 212 212 } 213 v->PutValues(nbins,rhigh+rlow,ehigh);213 if(splineFlag) v->FillSecondDerivatives(); 214 214 215 215 G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v); -
trunk/source/processes/electromagnetic/utils/src/G4LossTableManager.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4LossTableManager.cc,v 1.9 6 2009/04/09 16:10:57vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4LossTableManager.cc,v 1.97 2009/10/29 19:25:28 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 176 176 splineFlag = true; 177 177 bremsTh = DBL_MAX; 178 factorForAngleLimit = 1.0; 178 179 verbose = 1; 179 180 tableBuilder->SetSplineFlag(splineFlag); … … 920 921 } 921 922 923 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 924 925 void G4LossTableManager::SetFactorForAngleLimit(G4double val) 926 { 927 if(val > 0.0) factorForAngleLimit = val; 928 } 929 930 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 931 932 G4double G4LossTableManager::FactorForAngleLimit() const 933 { 934 return factorForAngleLimit; 935 } 936 922 937 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 923 938 -
trunk/source/processes/electromagnetic/utils/src/G4VEmFluctuationModel.cc
r1055 r1196 25 25 // 26 26 // $Id: G4VEmFluctuationModel.cc,v 1.4 2009/02/19 11:25:50 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/utils/src/G4VEmModel.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VEmModel.cc,v 1. 27 2009/05/26 15:00:49vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4VEmModel.cc,v 1.30 2009/09/23 14:42:47 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 62 62 G4VEmModel::G4VEmModel(const G4String& nam): 63 63 fluc(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV), 64 eMinActive(0.0),eMaxActive(DBL_MAX), 64 65 polarAngleLimit(0.0),secondaryThreshold(DBL_MAX),theLPMflag(false), 65 66 pParticleChange(0),nuclearStopping(false), … … 94 95 } else { 95 96 p = new G4ParticleChangeForLoss(); 97 pParticleChange = p; 96 98 } 97 99 return p; … … 107 109 } else { 108 110 p = new G4ParticleChangeForGamma(); 111 pParticleChange = p; 109 112 } 110 113 return p; … … 132 135 // initialise vector 133 136 for(G4int i=0; i<numOfCouples; i++) { 134 const G4MaterialCutsCouple* couple = 135 theCoupleTable->GetMaterialCutsCouple(i); 136 const G4Material* material = couple->GetMaterial(); 137 G4int idx = couple->GetIndex(); 137 currentCouple = theCoupleTable->GetMaterialCutsCouple(i); 138 const G4Material* material = currentCouple->GetMaterial(); 139 G4int idx = currentCouple->GetIndex(); 138 140 139 141 // selector already exist check if should be deleted … … 202 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 203 205 206 void G4VEmModel::DefineForRegion(const G4Region*) 207 {} 208 209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 210 204 211 G4double G4VEmModel::MinEnergyCut(const G4ParticleDefinition*, 205 212 const G4MaterialCutsCouple*) … … 249 256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 250 257 251 void G4VEmModel::DefineForRegion(const G4Region*)252 {}253 254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......255 256 258 void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*, 257 259 const G4Material*, G4double) -
trunk/source/processes/electromagnetic/utils/src/G4VEmProcess.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VEmProcess.cc,v 1. 66 2009/04/17 10:35:32vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4VEmProcess.cc,v 1.79 2009/11/10 20:30:55 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 52 52 // 12-04-07 remove double call to Clear model manager (V.Ivanchenko) 53 53 // 27-10-07 Virtual functions moved to source (V.Ivanchenko) 54 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko) 54 55 // 55 56 // Class Description: … … 172 173 { 173 174 G4int n = emModels.size(); 174 if(index >= n) for(G4int i=n; i<index+1; i++) {emModels.push_back(0);}175 if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} } 175 176 emModels[index] = p; 176 177 } … … 181 182 { 182 183 G4VEmModel* p = 0; 183 if(index >= 0 && index < G4int(emModels.size())) p = emModels[index];184 if(index >= 0 && index < G4int(emModels.size())) { p = emModels[index]; } 184 185 return p; 185 186 } … … 218 219 Clear(); 219 220 InitialiseProcess(particle); 221 222 // initialisation of models 223 G4int nmod = modelManager->NumberOfModels(); 224 for(G4int i=0; i<nmod; ++i) { 225 G4VEmModel* mod = modelManager->GetModel(i); 226 mod->SetPolarAngleLimit(polarAngleLimit); 227 if(mod->HighEnergyLimit() > maxKinEnergy) { 228 mod->SetHighEnergyLimit(maxKinEnergy); 229 } 230 } 231 220 232 theCuts = modelManager->Initialise(particle,secondaryParticle,2.,verboseLevel); 221 233 const G4ProductionCutsTable* theCoupleTable= … … 224 236 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut); 225 237 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut); 238 239 // prepare tables 226 240 if(buildLambdaTable){ 227 241 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable); … … 237 251 idxDERegions = new G4bool[numOfCouples]; 238 252 239 for (size_t j=0; j<numOfCouples; j++) {253 for (size_t j=0; j<numOfCouples; ++j) { 240 254 241 255 const G4MaterialCutsCouple* couple = … … 243 257 const G4ProductionCuts* pcuts = couple->GetProductionCuts(); 244 258 G4bool reg = false; 245 for(G4int i=0; i<nDERegions; i++) {259 for(G4int i=0; i<nDERegions; ++i) { 246 260 if(deRegions[i]) { 247 261 if(pcuts == deRegions[i]->GetProductionCuts()) reg = true; … … 253 267 if (1 < verboseLevel && nDERegions>0) { 254 268 G4cout << " Deexcitation is activated for regions: " << G4endl; 255 for (G4int i=0; i<nDERegions; i++) {269 for (G4int i=0; i<nDERegions; ++i) { 256 270 const G4Region* r = deRegions[i]; 257 271 G4cout << " " << r->GetName() << G4endl; … … 264 278 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part) 265 279 { 280 G4String partname = part.GetParticleName(); 266 281 if(1 < verboseLevel) { 267 282 G4cout << "G4VEmProcess::BuildPhysicsTable() for " 268 283 << GetProcessName() 269 << " and particle " << part .GetParticleName()284 << " and particle " << partname 270 285 << " buildLambdaTable= " << buildLambdaTable 271 286 << G4endl; … … 276 291 FindLambdaMax(); 277 292 } 278 if(0 < verboseLevel) PrintInfoDefinition(); 293 294 // reduce printout for nuclear stopping 295 G4bool gproc = true; 296 if(GetProcessName() == "nuclearStopping" && 297 partname != "GenericIon" && partname != "alpha") { gproc = false; } 298 299 if(gproc && 0 < verboseLevel) { PrintInfoDefinition(); } 279 300 280 301 if(1 < verboseLevel) { 281 302 G4cout << "G4VEmProcess::BuildPhysicsTable() done for " 282 303 << GetProcessName() 283 << " and particle " << part .GetParticleName()304 << " and particle " << partname 284 305 << G4endl; 285 306 } … … 301 322 G4ProductionCutsTable::GetProductionCutsTable(); 302 323 size_t numOfCouples = theCoupleTable->GetTableSize(); 303 for(size_t i=0; i<numOfCouples; i++) { 324 325 G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag(); 326 327 G4PhysicsLogVector* aVector = 0; 328 G4PhysicsLogVector* bVector = 0; 329 330 for(size_t i=0; i<numOfCouples; ++i) { 304 331 305 332 if (theLambdaTable->GetFlag(i)) { 306 333 307 334 // create physics vector and fill it 308 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i); 309 G4PhysicsVector* aVector = LambdaPhysicsVector(couple); 335 const G4MaterialCutsCouple* couple = 336 theCoupleTable->GetMaterialCutsCouple(i); 337 if(!bVector) { 338 aVector = 339 static_cast<G4PhysicsLogVector*>(LambdaPhysicsVector(couple)); 340 bVector = aVector; 341 } else { 342 aVector = new G4PhysicsLogVector(*bVector); 343 } 344 // G4PhysicsVector* aVector = LambdaPhysicsVector(couple); 345 aVector->SetSpline(splineFlag); 310 346 modelManager->FillLambdaVector(aVector, couple, startFromNull); 347 if(splineFlag) aVector->FillSecondDerivatives(); 311 348 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector); 312 349 } … … 364 401 if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0; 365 402 InitialiseStep(track); 403 if(!currentModel->IsActive(preStepKinEnergy)) return x; 366 404 367 405 if(preStepKinEnergy < mfpKinEnergy) { … … 444 482 } 445 483 446 SelectModel(finalT); 484 SelectModel(finalT, currentCoupleIndex); 485 if(!currentModel->IsActive(finalT)) return &fParticleChange; 447 486 if(useDeexcitation) { 448 currentModel->SetDeexcitationFlag(idxDERegions[current MaterialIndex]);487 currentModel->SetDeexcitationFlag(idxDERegions[currentCoupleIndex]); 449 488 } 450 489 /* … … 464 503 currentCouple, 465 504 track.GetDynamicParticle(), 466 (*theCuts)[current MaterialIndex]);505 (*theCuts)[currentCoupleIndex]); 467 506 468 507 // save secondaries … … 473 512 G4double edep = fParticleChange.GetLocalEnergyDeposit(); 474 513 475 for (G4int i=0; i<num; i++) {514 for (G4int i=0; i<num; ++i) { 476 515 G4DynamicParticle* dp = secParticles[i]; 477 516 const G4ParticleDefinition* p = dp->GetDefinition(); … … 480 519 if(applyCuts) { 481 520 if (p == theGamma) { 482 if (e < (*theCutsGamma)[current MaterialIndex]) good = false;521 if (e < (*theCutsGamma)[currentCoupleIndex]) good = false; 483 522 484 523 } else if (p == theElectron) { 485 if (e < (*theCutsElectron)[current MaterialIndex]) good = false;524 if (e < (*theCutsElectron)[currentCoupleIndex]) good = false; 486 525 487 526 } else if (p == thePositron) { 488 if (electron_mass_c2 < (*theCutsGamma)[current MaterialIndex] &&489 e < (*theCutsPositron)[current MaterialIndex]) {527 if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] && 528 e < (*theCutsPositron)[currentCoupleIndex]) { 490 529 good = false; 491 530 e += 2.0*electron_mass_c2; … … 539 578 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part, 540 579 const G4String& directory, 541 580 G4bool ascii) 542 581 { 543 582 if(1 < verboseLevel) { … … 558 597 if ( yes ) { 559 598 if (0 < verboseLevel) { 560 G4cout << "Lambda table for " << particleName << " is Retrieved from <" 599 G4cout << "Lambda table for " << particleName 600 << " is Retrieved from <" 561 601 << filename << ">" 562 602 << G4endl; … … 564 604 if((G4LossTableManager::Instance())->SplineFlag()) { 565 605 size_t n = theLambdaTable->length(); 566 for(size_t i=0; i<n; i++) {(* theLambdaTable)[i]->SetSpline(true);} 606 for(size_t i=0; i<n; ++i) { 607 if((* theLambdaTable)[i]) { 608 (* theLambdaTable)[i]->SetSpline(true); 609 } 610 } 567 611 } 568 612 } else { … … 587 631 // the region is in the list 588 632 if (nDERegions) { 589 for (G4int i=0; i<nDERegions; i++) {633 for (G4int i=0; i<nDERegions; ++i) { 590 634 if (reg == deRegions[i]) { 591 635 if(!val) deRegions[i] = 0; … … 613 657 DefineMaterial(couple); 614 658 G4double cross = 0.0; 615 G4bool b;616 659 if(theLambdaTable) { 617 cross = (((*theLambdaTable)[currentMaterialIndex])-> 618 GetValue(kineticEnergy, b)); 660 cross = (((*theLambdaTable)[currentCoupleIndex])->Value(kineticEnergy)); 619 661 } else { 620 SelectModel(kineticEnergy );662 SelectModel(kineticEnergy, currentCoupleIndex); 621 663 cross = currentModel->CrossSectionPerVolume(currentMaterial, 622 664 particle,kineticEnergy); … … 650 692 theEnergyOfCrossSectionMax = new G4double [n]; 651 693 theCrossSectionMax = new G4double [n]; 652 G4bool b; 653 654 for (size_t i=0; i<n; i++) { 694 695 for (size_t i=0; i<n; ++i) { 655 696 pv = (*theLambdaTable)[i]; 656 697 emax = DBL_MAX; … … 658 699 if(pv) { 659 700 size_t nb = pv->GetVectorLength(); 660 emax = pv->GetLowEdgeEnergy(nb);701 emax = DBL_MAX; 661 702 smax = 0.0; 662 for (size_t j=0; j<nb; j++) { 663 e = pv->GetLowEdgeEnergy(j); 664 s = pv->GetValue(e,b); 665 if(s > smax) { 666 smax = s; 667 emax = e; 703 if(nb > 0) { 704 for (size_t j=0; j<nb; ++j) { 705 e = pv->Energy(j); 706 s = (*pv)(j); 707 if(s > smax) { 708 smax = s; 709 emax = e; 710 } 668 711 } 669 712 } -
trunk/source/processes/electromagnetic/utils/src/G4VEnergyLoss.cc
r1007 r1196 26 26 // 27 27 // $Id: G4VEnergyLoss.cc,v 1.46 2006/06/29 19:55:21 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 -
trunk/source/processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VEnergyLossProcess.cc,v 1.1 49 2009/04/17 10:35:32vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4VEnergyLossProcess.cc,v 1.158 2009/10/29 18:07:08 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 107 107 // 25-04-07 move initialisation of safety helper to BuildPhysicsTable (VI) 108 108 // 27-10-07 Virtual functions moved to source (V.Ivanchenko) 109 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko) 109 110 // 110 111 // Class Description: … … 230 231 vstrag->SetSpline(true); 231 232 G4double s[nrbins] = {-0.2, -0.85, -1.3, -1.578, -1.76, -1.85, -1.9}; 232 for(G4int i=0; i<nrbins; i++) {vstrag->PutValue(i, s[i]);}233 for(G4int i=0; i<nrbins; ++i) {vstrag->PutValue(i, s[i]);} 233 234 } 234 235 … … 351 352 { 352 353 G4int n = emModels.size(); 353 if(index >= n) for(G4int i=n; i<index+1; i++) {emModels.push_back(0);}354 if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} } 354 355 emModels[index] = p; 355 356 } … … 360 361 { 361 362 G4VEmModel* p = 0; 362 if(index >= 0 && index < G4int(emModels.size())) p = emModels[index];363 if(index >= 0 && index < G4int(emModels.size())) { p = emModels[index]; } 363 364 return p; 364 365 } … … 466 467 G4double q = initialCharge/baseParticle->GetPDGCharge(); 467 468 chargeSqRatio = q*q; 468 if(chargeSqRatio > 0.0) reduceFactor = 1.0/(chargeSqRatio*massRatio); 469 if(chargeSqRatio > 0.0) { reduceFactor = 1.0/(chargeSqRatio*massRatio); } 470 } 471 472 // initialisation of models 473 G4int nmod = modelManager->NumberOfModels(); 474 for(G4int i=0; i<nmod; ++i) { 475 G4VEmModel* mod = modelManager->GetModel(i); 476 if(mod->HighEnergyLimit() > maxKinEnergy) { 477 mod->SetHighEnergyLimit(maxKinEnergy); 478 } 469 479 } 470 480 … … 483 493 if(nDERegions>0) idxDERegions = new G4bool[numOfCouples]; 484 494 485 for (size_t j=0; j<numOfCouples; j++) {495 for (size_t j=0; j<numOfCouples; ++j) { 486 496 487 497 const G4MaterialCutsCouple* couple = … … 491 501 if(nSCoffRegions>0) { 492 502 G4bool reg = false; 493 for(G4int i=0; i<nSCoffRegions; i++) {503 for(G4int i=0; i<nSCoffRegions; ++i) { 494 504 if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = true; 495 505 } … … 498 508 if(nDERegions>0) { 499 509 G4bool reg = false; 500 for(G4int i=0; i<nDERegions; i++) {510 for(G4int i=0; i<nDERegions; ++i) { 501 511 if( pcuts == deRegions[i]->GetProductionCuts()) reg = true; 502 512 } … … 517 527 if (nSCoffRegions) { 518 528 G4cout << " SubCutoff Regime is ON for regions: " << G4endl; 519 for (G4int i=0; i<nSCoffRegions; i++) {529 for (G4int i=0; i<nSCoffRegions; ++i) { 520 530 const G4Region* r = scoffRegions[i]; 521 531 G4cout << " " << r->GetName() << G4endl; … … 524 534 if (nDERegions) { 525 535 G4cout << " Deexcitation is ON for regions: " << G4endl; 526 for (G4int i=0; i<nDERegions; i++) {536 for (G4int i=0; i<nDERegions; ++i) { 527 537 const G4Region* r = deRegions[i]; 528 538 G4cout << " " << r->GetName() << G4endl; … … 556 566 } 557 567 } 568 569 // Added tracking cut to avoid tracking artifacts 570 if(isIonisation) fParticleChange.SetLowEnergyLimit(lowestKinEnergy); 558 571 559 572 if(1 < verboseLevel) { … … 613 626 if(!table) return table; 614 627 615 for(size_t i=0; i<numOfCouples; i++) { 628 G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag(); 629 G4PhysicsLogVector* aVector = 0; 630 G4PhysicsLogVector* bVector = 0; 631 632 for(size_t i=0; i<numOfCouples; ++i) { 616 633 617 634 if(1 < verboseLevel) { … … 624 641 const G4MaterialCutsCouple* couple = 625 642 theCoupleTable->GetMaterialCutsCouple(i); 626 G4PhysicsVector* aVector = new G4PhysicsLogVector(emin, emax, bin); 627 aVector->SetSpline((G4LossTableManager::Instance())->SplineFlag()); 643 if(!bVector) { 644 aVector = new G4PhysicsLogVector(emin, emax, bin); 645 bVector = aVector; 646 } else { 647 aVector = new G4PhysicsLogVector(*bVector); 648 } 649 // G4PhysicsVector* aVector = new G4PhysicsLogVector(emin, emax, bin); 650 aVector->SetSpline(splineFlag); 628 651 629 652 modelManager->FillDEDXVector(aVector, couple, tType); 653 if(splineFlag) aVector->FillSecondDerivatives(); 630 654 631 655 // Insert vector for this material into the table … … 676 700 size_t numOfCouples = theCoupleTable->GetTableSize(); 677 701 678 for(size_t i=0; i<numOfCouples; i++) { 702 G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag(); 703 704 for(size_t i=0; i<numOfCouples; ++i) { 679 705 680 706 if (table->GetFlag(i)) { … … 686 712 if(fSubRestricted == tType) cut = (*theSubCuts)[i]; 687 713 G4PhysicsVector* aVector = LambdaPhysicsVector(couple, cut); 714 aVector->SetSpline(splineFlag); 715 688 716 modelManager->FillLambdaVector(aVector, couple, true, tType); 717 if(splineFlag) aVector->FillSecondDerivatives(); 689 718 690 719 // Insert vector for this material into the table … … 787 816 // the region is in the list 788 817 if (nSCoffRegions) { 789 for (G4int i=0; i<nSCoffRegions; i++) {818 for (G4int i=0; i<nSCoffRegions; ++i) { 790 819 if (reg == scoffRegions[i]) { 791 820 if(!val) deRegions[i] = 0; … … 799 828 useSubCutoff = true; 800 829 scoffRegions.push_back(reg); 801 nSCoffRegions++;830 ++nSCoffRegions; 802 831 } else { 803 832 useSubCutoff = false; … … 815 844 // the region is in the list 816 845 if (nDERegions) { 817 for (G4int i=0; i<nDERegions; i++) {846 for (G4int i=0; i<nDERegions; ++i) { 818 847 if (reg == deRegions[i]) { 819 848 if(!val) deRegions[i] = 0; … … 827 856 useDeexcitation = true; 828 857 deRegions.push_back(reg); 829 nDERegions++;858 ++nDERegions; 830 859 } else { 831 860 useDeexcitation = false; … … 883 912 preStepScaledEnergy = preStepKinEnergy*massRatio; 884 913 SelectModel(preStepScaledEnergy); 914 if(!currentModel->IsActive(preStepScaledEnergy)) return x; 885 915 886 916 if(isIon) { … … 951 981 fParticleChange.InitializeForAlongStep(track); 952 982 // The process has range table - calculate energy loss 953 if(!isIonisation) return &fParticleChange; 983 if(!isIonisation || !currentModel->IsActive(preStepScaledEnergy)) { 984 return &fParticleChange; 985 } 954 986 955 987 // Get the actual (true) Step length … … 1076 1108 /* 1077 1109 if(nProcesses > 0) { 1078 for(G4int i=0; i<nProcesses; i++) {1110 for(G4int i=0; i<nProcesses; ++i) { 1079 1111 (scProcesses[i])->SampleSubCutSecondaries( 1080 1112 scTracks, step, (scProcesses[i])-> … … 1088 1120 G4ThreeVector mom = dynParticle->GetMomentum(); 1089 1121 fParticleChange.SetNumberOfSecondaries(n); 1090 for(G4int i=0; i<n; i++) {1122 for(G4int i=0; i<n; ++i) { 1091 1123 G4Track* t = scTracks[i]; 1092 1124 G4double e = t->GetKineticEnergy(); … … 1182 1214 const G4DynamicParticle* dp = track->GetDynamicParticle(); 1183 1215 G4double e = dp->GetKineticEnergy()*massRatio; 1184 G4bool b; 1185 G4double cross = chargeSqRatio*(((*theSubLambdaTable)[idx])->GetValue(e,b)); 1216 G4double cross = chargeSqRatio*(((*theSubLambdaTable)[idx])->Value(e)); 1186 1217 G4double length = step.GetStepLength(); 1187 1218 … … 1220 1251 G4ThreeVector r = prepoint + fragment*dr; 1221 1252 std::vector<G4DynamicParticle*>::iterator it; 1222 for(it=secParticles.begin(); it!=secParticles.end(); it++) {1253 for(it=secParticles.begin(); it!=secParticles.end(); ++it) { 1223 1254 1224 1255 G4bool addSec = true; … … 1226 1257 // do not track very low-energy delta-electrons 1227 1258 if(theSecondaryRangeTable && (*it)->GetDefinition() == theElectron) { 1228 G4bool b;1229 1259 G4double ekin = (*it)->GetKineticEnergy(); 1230 G4double rg = ((*theSecondaryRangeTable)[idx]-> GetValue(ekin, b));1260 G4double rg = ((*theSecondaryRangeTable)[idx]->Value(ekin)); 1231 1261 // if(rg < currentMinSafety) { 1232 1262 if(rg < safetyHelper->ComputeSafety(r)) { … … 1265 1295 1266 1296 G4double postStepScaledEnergy = finalT*massRatio; 1297 1298 if(!currentModel->IsActive(postStepScaledEnergy)) return &fParticleChange; 1267 1299 /* 1268 1300 if(-1 < verboseLevel) { … … 1283 1315 << " < " << lx << " (postLambda) " 1284 1316 << G4endl; 1285 nWarnings++;1317 ++nWarnings; 1286 1318 } 1287 1319 */ … … 1308 1340 if(num > 0) { 1309 1341 fParticleChange.SetNumberOfSecondaries(num); 1310 for (G4int i=0; i<num; i++) {1342 for (G4int i=0; i<num; ++i) { 1311 1343 fParticleChange.AddSecondary(secParticles[i]); 1312 1344 } … … 1390 1422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 1391 1423 1392 G4bool G4VEnergyLossProcess::RetrievePhysicsTable( 1393 const G4ParticleDefinition* part, const G4String& directory, 1394 G4bool ascii) 1424 G4bool 1425 G4VEnergyLossProcess::RetrievePhysicsTable(const G4ParticleDefinition* part, 1426 const G4String& directory, 1427 G4bool ascii) 1395 1428 { 1396 1429 G4bool res = true; … … 1405 1438 if(particle == part) { 1406 1439 1407 // G4bool yes = true;1408 1440 if ( !baseParticle ) { 1409 1441 … … 1412 1444 {fpi = false;} 1413 1445 1414 if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false)) 1446 // ionisation table keeps individual dEdx and not sum of sub-processes 1447 if(!RetrieveTable(part,theDEDXTable,ascii,directory,"Ionisation",false)) 1415 1448 {fpi = false;} 1416 1449 … … 1440 1473 1441 1474 if(!fpi) yes = false; 1442 if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes)) 1475 if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes)) 1443 1476 {res = false;} 1444 1477 } … … 1465 1498 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 1466 1499 1467 G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part, 1468 G4PhysicsTable* aTable, G4bool ascii, 1469 const G4String& directory, 1470 const G4String& tname, 1471 G4bool mandatory) 1500 G4bool 1501 G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part, 1502 G4PhysicsTable* aTable, 1503 G4bool ascii, 1504 const G4String& directory, 1505 const G4String& tname, 1506 G4bool mandatory) 1472 1507 { 1473 1508 G4bool res = true; … … 1475 1510 G4bool yes = aTable->ExistPhysicsTable(filename); 1476 1511 if(yes) { 1512 if(!aTable) aTable = G4PhysicsTableHelper::PreparePhysicsTable(0); 1477 1513 yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii); 1514 1478 1515 if((G4LossTableManager::Instance())->SplineFlag()) { 1479 1516 size_t n = aTable->length(); 1480 for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);} 1517 for(size_t i=0; i<n; ++i) { 1518 if((*aTable)[i]) { 1519 (*aTable)[i]->SetSpline(true); 1520 } 1521 } 1481 1522 } 1482 1523 } … … 1525 1566 DefineMaterial(couple); 1526 1567 G4double cross = 0.0; 1527 G4bool b;1528 1568 if(theLambdaTable) { 1529 cross = 1530 ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b); 1569 cross = ((*theLambdaTable)[currentMaterialIndex])->Value(kineticEnergy); 1531 1570 } else { 1532 1571 SelectModel(kineticEnergy); 1533 cross = 1534 currentModel->CrossSectionPerVolume(currentMaterial, 1535 particle, kineticEnergy, 1536 (*theCuts)[currentMaterialIndex]); 1572 cross = currentModel->CrossSectionPerVolume(currentMaterial, 1573 particle, kineticEnergy, 1574 (*theCuts)[currentMaterialIndex]); 1537 1575 } 1538 1576 … … 1587 1625 const G4MaterialCutsCouple* couple, G4double cut) 1588 1626 { 1589 // G4double cut = (*theCuts)[couple->GetIndex()];1590 // G4int nbins = nLambdaBins;1591 1627 G4double tmin = 1592 1628 std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut), … … 1595 1631 G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins); 1596 1632 v->SetSpline((G4LossTableManager::Instance())->SplineFlag()); 1597 1598 1633 return v; 1599 1634 } … … 1607 1642 if(p->GetProcessName() != "eBrem") add = false; 1608 1643 if(add && nProcesses > 0) { 1609 for(G4int i=0; i<nProcesses; i++) {1644 for(G4int i=0; i<nProcesses; ++i) { 1610 1645 if(p == scProcesses[i]) { 1611 1646 add = false; … … 1616 1651 if(add) { 1617 1652 scProcesses.push_back(p); 1618 nProcesses++;1653 ++nProcesses; 1619 1654 if (1 < verboseLevel) { 1620 1655 G4cout << "### The process " << p->GetProcessName() … … 1636 1671 G4PhysicsVector* pv = (*p)[0]; 1637 1672 G4double emax = maxKinEnergyCSDA; 1638 G4bool b;1639 1673 theDEDXAtMaxEnergy = new G4double [n]; 1640 1674 1641 for (size_t i=0; i<n; i++) {1675 for (size_t i=0; i<n; ++i) { 1642 1676 pv = (*p)[i]; 1643 G4double dedx = pv-> GetValue(emax, b);1677 G4double dedx = pv->Value(emax); 1644 1678 theDEDXAtMaxEnergy[i] = dedx; 1645 1679 //G4cout << "i= " << i << " emax(MeV)= " << emax/MeV<< " dedx= " … … 1671 1705 G4PhysicsVector* pv = (*p)[0]; 1672 1706 G4double emax = maxKinEnergyCSDA; 1673 G4bool b;1674 1707 theRangeAtMaxEnergy = new G4double [n]; 1675 1708 1676 for (size_t i=0; i<n; i++) {1709 for (size_t i=0; i<n; ++i) { 1677 1710 pv = (*p)[i]; 1678 G4double r2 = pv-> GetValue(emax, b);1711 G4double r2 = pv->Value(emax); 1679 1712 theRangeAtMaxEnergy[i] = r2; 1680 1713 //G4cout << "i= " << i << " e2(MeV)= " << emax/MeV << " r2= " … … 1744 1777 theEnergyOfCrossSectionMax = new G4double [n]; 1745 1778 theCrossSectionMax = new G4double [n]; 1746 G4bool b; 1747 1748 for (size_t i=0; i<n; i++) { 1779 1780 for (size_t i=0; i<n; ++i) { 1749 1781 pv = (*p)[i]; 1750 1782 emax = DBL_MAX; … … 1752 1784 if(pv) { 1753 1785 size_t nb = pv->GetVectorLength(); 1754 emax = pv->GetLowEdgeEnergy(nb); 1755 for (size_t j=0; j<nb; j++) { 1756 e = pv->GetLowEdgeEnergy(j); 1757 s = pv->GetValue(e,b); 1758 if(s > smax) { 1759 smax = s; 1760 emax = e; 1786 if(nb > 0) { 1787 for (size_t j=0; j<nb; ++j) { 1788 e = pv->Energy(j); 1789 s = (*pv)(j); 1790 if(s > smax) { 1791 smax = s; 1792 emax = e; 1793 } 1761 1794 } 1762 1795 } -
trunk/source/processes/electromagnetic/utils/src/G4VMscModel.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VMscModel.cc,v 1.1 2 2009/05/10 19:36:19vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4VMscModel.cc,v 1.13 2009/07/20 17:32:47 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 60 60 facrange(0.04), 61 61 facgeom(2.5), 62 facsafety(0. 25),62 facsafety(0.3), 63 63 skin(3.0), 64 64 dtrl(0.05), -
trunk/source/processes/electromagnetic/utils/src/G4VMultipleScattering.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4VMultipleScattering.cc,v 1. 66 2009/05/27 11:55:02vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4VMultipleScattering.cc,v 1.77 2009/10/29 18:07:08 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 56 56 // 27-10-07 Virtual functions moved to source (V.Ivanchenko) 57 57 // 11-03-08 Set skin value does not effect step limit type (V.Ivanchenko) 58 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko) 58 59 // 59 60 // Class Description: … … 98 99 facgeom(2.5), 99 100 latDisplasment(true), 101 isIon(false), 100 102 currentParticle(0), 101 103 currentCouple(0) … … 144 146 if(p) p->SetParticleChange(pParticleChange); 145 147 } 146 147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 148 149 G4VEmModel* G4VMultipleScattering::GetModelByIndex(G4int idx, G4bool ver) 148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 149 150 void G4VMultipleScattering::SetModel(G4VMscModel* p, G4int index) 151 { 152 G4int n = mscModels.size(); 153 if(index >= n) { for(G4int i=n; i<=index; ++i) {mscModels.push_back(0);} } 154 mscModels[index] = p; 155 } 156 157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 158 159 G4VMscModel* G4VMultipleScattering::Model(G4int index) 160 { 161 G4VMscModel* p = 0; 162 if(index >= 0 && index < G4int(mscModels.size())) { p = mscModels[index]; } 163 return p; 164 } 165 166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 167 168 G4VEmModel* 169 G4VMultipleScattering::GetModelByIndex(G4int idx, G4bool ver) const 150 170 { 151 171 return modelManager->GetModel(idx, ver); … … 170 190 size_t numOfCouples = theCoupleTable->GetTableSize(); 171 191 172 for (size_t i=0; i<numOfCouples; i++) { 192 G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag(); 193 194 G4PhysicsLogVector* aVector = 0; 195 G4PhysicsLogVector* bVector = 0; 196 197 for (size_t i=0; i<numOfCouples; ++i) { 173 198 174 199 if (theLambdaTable->GetFlag(i)) { … … 176 201 const G4MaterialCutsCouple* couple = 177 202 theCoupleTable->GetMaterialCutsCouple(i); 178 G4PhysicsVector* aVector = PhysicsVector(couple); 203 if(!bVector) { 204 aVector = static_cast<G4PhysicsLogVector*>(PhysicsVector(couple)); 205 bVector = aVector; 206 } else { 207 aVector = new G4PhysicsLogVector(*bVector); 208 } 209 //G4PhysicsVector* aVector = PhysicsVector(couple); 210 aVector->SetSpline(splineFlag); 179 211 modelManager->FillLambdaVector(aVector, couple, false); 212 if(splineFlag) aVector->FillSecondDerivatives(); 180 213 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector); 181 214 } … … 212 245 part.GetParticleSubType() == "generic") { 213 246 firstParticle = G4GenericIon::GenericIon(); 247 isIon = true; 214 248 } else { 215 249 firstParticle = ∂ 250 if(part.GetParticleType() == "nucleus" || 251 part.GetPDGMass() > GeV) {isIon = true;} 216 252 } 217 253 // limitations for ions 254 if(isIon) { 255 SetStepLimitType(fMinimal); 256 SetLateralDisplasmentFlag(false); 257 SetBuildLambdaTable(false); 258 } 218 259 currentParticle = ∂ 219 260 } … … 232 273 233 274 InitialiseProcess(firstParticle); 234 if(buildLambdaTable) 275 276 // initialisation of models 277 G4int nmod = modelManager->NumberOfModels(); 278 for(G4int i=0; i<nmod; ++i) { 279 G4VMscModel* msc = static_cast<G4VMscModel*>(modelManager->GetModel(i)); 280 if(isIon) { 281 msc->SetStepLimitType(fMinimal); 282 msc->SetLateralDisplasmentFlag(false); 283 msc->SetRangeFactor(0.2); 284 } else { 285 msc->SetStepLimitType(StepLimitType()); 286 msc->SetLateralDisplasmentFlag(LateralDisplasmentFlag()); 287 msc->SetSkin(Skin()); 288 msc->SetRangeFactor(RangeFactor()); 289 msc->SetGeomFactor(GeomFactor()); 290 } 291 msc->SetPolarAngleLimit(polarAngleLimit); 292 if(msc->HighEnergyLimit() > maxKinEnergy) { 293 msc->SetHighEnergyLimit(maxKinEnergy); 294 } 295 } 296 297 modelManager->Initialise(firstParticle, G4Electron::Electron(), 298 10.0, verboseLevel); 299 300 // prepare tables 301 if(buildLambdaTable) { 235 302 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable); 236 const G4DataVector* theCuts = 237 modelManager->Initialise(firstParticle, 238 G4Electron::Electron(), 239 10.0, verboseLevel); 240 241 if(2 < verboseLevel) G4cout << theCuts << G4endl; 242 303 } 243 304 } 244 305 } … … 277 338 G4double, 278 339 G4double currentMinimalStep, 279 G4double& currentSafety,340 G4double&, 280 341 G4GPILSelection* selection) 281 342 { 282 343 // get Step limit proposed by the process 283 valueGPILSelectionMSC = NotCandidateForSelection; 284 G4double steplength = GetMscContinuousStepLimit(track, 285 track.GetKineticEnergy(), 286 currentMinimalStep, 287 currentSafety); 288 // G4cout << "StepLimit= " << steplength << G4endl; 289 // set return value for G4GPILSelection 290 *selection = valueGPILSelectionMSC; 291 return steplength; 292 } 293 294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 295 296 G4double G4VMultipleScattering::PostStepGetPhysicalInteractionLength( 297 const G4Track&, G4double, G4ForceCondition* condition) 298 { 299 *condition = Forced; 300 return DBL_MAX; 344 *selection = NotCandidateForSelection; 345 G4double x = currentMinimalStep; 346 DefineMaterial(track.GetMaterialCutsCouple()); 347 G4double ekin = track.GetKineticEnergy(); 348 if(isIon) { ekin *= proton_mass_c2/track.GetDefinition()->GetPDGMass(); } 349 currentModel = static_cast<G4VMscModel*>(SelectModel(ekin)); 350 if(x > 0.0 && ekin > 0.0 && currentModel->IsActive(ekin)) { 351 G4double tPathLength = 352 currentModel->ComputeTruePathLengthLimit(track, theLambdaTable, x); 353 if (tPathLength < x) *selection = CandidateForSelection; 354 x = currentModel->ComputeGeomPathLength(tPathLength); 355 // G4cout << "tPathLength= " << tPathLength 356 // << " stepLimit= " << x 357 // << " currentMinimalStep= " << currentMinimalStep<< G4endl; 358 } 359 return x; 301 360 } 302 361 … … 309 368 G4double& currentSafety) 310 369 { 311 return GetMscContinuousStepLimit(track,previousStepSize,currentMinimalStep, 312 currentSafety); 370 G4GPILSelection* selection = 0; 371 return AlongStepGetPhysicalInteractionLength(track,previousStepSize,currentMinimalStep, 372 currentSafety, selection); 313 373 } 314 374 … … 320 380 *condition = Forced; 321 381 return DBL_MAX; 322 }323 324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....325 326 G4VParticleChange* G4VMultipleScattering::AlongStepDoIt(const G4Track&,327 const G4Step& step)328 {329 fParticleChange.ProposeTrueStepLength(330 currentModel->ComputeTrueStepLength(step.GetStepLength()));331 return &fParticleChange;332 }333 334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....335 336 G4VParticleChange* G4VMultipleScattering::PostStepDoIt(const G4Track& track,337 const G4Step& step)338 {339 fParticleChange.Initialize(track);340 currentModel->SampleScattering(track.GetDynamicParticle(),341 step.GetPostStepPoint()->GetSafety());342 return &fParticleChange;343 382 } 344 383 … … 384 423 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 385 424 386 G4bool G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition* part, 387 const G4String& directory, 388 G4bool ascii) 425 G4bool 426 G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition* part, 427 const G4String& directory, 428 G4bool ascii) 389 429 { 390 430 if(0 < verboseLevel) { 391 // G4cout << "========================================================" << G4endl;392 431 G4cout << "G4VMultipleScattering::RetrievePhysicsTable() for " 393 432 << part->GetParticleName() << " and process " … … 401 440 402 441 G4String filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii); 403 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,filename,ascii); 442 yes = 443 G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,filename,ascii); 404 444 if ( yes ) { 405 445 if (0 < verboseLevel) { 406 G4cout << "Lambda table for " << part->GetParticleName() << " is retrieved from <" 446 G4cout << "Lambda table for " << part->GetParticleName() 447 << " is retrieved from <" 407 448 << filename << ">" 408 449 << G4endl; … … 410 451 if((G4LossTableManager::Instance())->SplineFlag()) { 411 452 size_t n = theLambdaTable->length(); 412 for(size_t i=0; i<n; i++) {(* theLambdaTable)[i]->SetSpline(true);} 453 for(size_t i=0; i<n; ++i) { 454 if((* theLambdaTable)[i]) { 455 (* theLambdaTable)[i]->SetSpline(true); 456 } 457 } 413 458 } 414 459 } else { 415 460 if (1 < verboseLevel) { 416 G4cout << "Lambda table for " << part->GetParticleName() << " in file <" 461 G4cout << "Lambda table for " << part->GetParticleName() 462 << " in file <" 417 463 << filename << "> is not exist" 418 464 << G4endl; -
trunk/source/processes/electromagnetic/utils/src/G4ionEffectiveCharge.cc
r1007 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4ionEffectiveCharge.cc,v 1.2 4 2008/12/18 13:01:46 gunterExp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4ionEffectiveCharge.cc,v 1.25 2009/10/29 16:57:39 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 105 105 // Fast ions or hadrons 106 106 G4double reducedEnergy = kineticEnergy * proton_mass_c2/mass ; 107 108 //G4cout << "e= " << reducedEnergy << " Zi= " << Zi << " " << material->GetName() << G4endl; 109 107 110 if( reducedEnergy > Zi*energyHighLimit || Zi < 1.5 || !material) return charge; 108 111 109 112 G4double z = material->GetIonisation()->GetZeffective(); 110 //reducedEnergy = std::max(reducedEnergy,energyLowLimit);113 reducedEnergy = std::max(reducedEnergy,energyLowLimit); 111 114 112 115 // Helium ion case
Note: See TracChangeset
for help on using the changeset viewer.