Changeset 1055 for trunk/source/processes/electromagnetic/standard
- Timestamp:
- May 28, 2009, 4:26:57 PM (15 years ago)
- Location:
- trunk/source/processes/electromagnetic/standard
- Files:
-
- 59 edited
-
History (modified) (2 diffs)
-
include/G4BetheBlochModel.hh (modified) (5 diffs)
-
include/G4BohrFluctuations.hh (modified) (3 diffs)
-
include/G4BraggIonModel.hh (modified) (4 diffs)
-
include/G4BraggModel.hh (modified) (4 diffs)
-
include/G4ComptonScattering.hh (modified) (3 diffs)
-
include/G4CoulombScattering.hh (modified) (3 diffs)
-
include/G4GammaConversion.hh (modified) (3 diffs)
-
include/G4IonFluctuations.hh (modified) (3 diffs)
-
include/G4MollerBhabhaModel.hh (modified) (4 diffs)
-
include/G4PAIModel.hh (modified) (4 diffs)
-
include/G4PAIPhotonModel.hh (modified) (5 diffs)
-
include/G4PEEffectModel.hh (modified) (4 diffs)
-
include/G4PhotoElectricEffect.hh (modified) (3 diffs)
-
include/G4UniversalFluctuation.hh (modified) (4 diffs)
-
include/G4UrbanMscModel.hh (modified) (12 diffs)
-
include/G4UrbanMscModel2.hh (modified) (6 diffs)
-
include/G4UrbanMscModel90.hh (modified) (2 diffs)
-
include/G4WentzelVIModel.hh (modified) (12 diffs)
-
include/G4eBremsstrahlung.hh (modified) (3 diffs)
-
include/G4eBremsstrahlungModel.hh (modified) (4 diffs)
-
include/G4eBremsstrahlungRelModel.hh (modified) (4 diffs)
-
include/G4eCoulombScatteringModel.hh (modified) (8 diffs)
-
include/G4eIonisation.hh (modified) (4 diffs)
-
include/G4eplusAnnihilation.hh (modified) (4 diffs)
-
include/G4hIonisation.hh (modified) (3 diffs)
-
include/G4ionIonisation.hh (modified) (4 diffs)
-
src/G4BetheBlochModel.cc (modified) (8 diffs)
-
src/G4BetheHeitlerModel.cc (modified) (2 diffs)
-
src/G4BohrFluctuations.cc (modified) (2 diffs)
-
src/G4BraggIonModel.cc (modified) (4 diffs)
-
src/G4BraggModel.cc (modified) (4 diffs)
-
src/G4ComptonScattering.cc (modified) (2 diffs)
-
src/G4CoulombScattering.cc (modified) (2 diffs)
-
src/G4CoulombScatteringModel.cc (modified) (4 diffs)
-
src/G4GammaConversion.cc (modified) (2 diffs)
-
src/G4IonFluctuations.cc (modified) (3 diffs)
-
src/G4KleinNishinaCompton.cc (modified) (2 diffs)
-
src/G4MollerBhabhaModel.cc (modified) (7 diffs)
-
src/G4MscModel71.cc (modified) (2 diffs)
-
src/G4PAIModel.cc (modified) (10 diffs)
-
src/G4PAIPhotonModel.cc (modified) (13 diffs)
-
src/G4PEEffectModel.cc (modified) (5 diffs)
-
src/G4PhotoElectricEffect.cc (modified) (2 diffs)
-
src/G4UniversalFluctuation.cc (modified) (8 diffs)
-
src/G4UrbanMscModel.cc (modified) (30 diffs)
-
src/G4UrbanMscModel2.cc (modified) (13 diffs)
-
src/G4UrbanMscModel90.cc (modified) (7 diffs)
-
src/G4WaterStopping.cc (modified) (4 diffs)
-
src/G4WentzelVIModel.cc (modified) (23 diffs)
-
src/G4eBremsstrahlung.cc (modified) (2 diffs)
-
src/G4eBremsstrahlungModel.cc (modified) (2 diffs)
-
src/G4eBremsstrahlungRelModel.cc (modified) (2 diffs)
-
src/G4eCoulombScatteringModel.cc (modified) (9 diffs)
-
src/G4eIonisation.cc (modified) (2 diffs)
-
src/G4eeToTwoGammaModel.cc (modified) (2 diffs)
-
src/G4eplusAnnihilation.cc (modified) (2 diffs)
-
src/G4hIonisation.cc (modified) (2 diffs)
-
src/G4ionIonisation.cc (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/electromagnetic/standard/History
r1007 r1055 1 $Id: History,v 1.4 30 2008/11/24 18:28:09vnivanch Exp $1 $Id: History,v 1.452 2009/05/18 12:31:46 vnivanch Exp $ 2 2 ------------------------------------------------------------------- 3 3 … … 17 17 * Reverse chronological order (last date on top), please * 18 18 ---------------------------------------------------------- 19 20 18 May 09: V.Ivant (emstand-V09-02-12) 21 - G4UrbanMscModel2 - (L.Urban) correction in tail tuning using MuScat data 22 23 15 May 09: V.Ivant (emstand-V09-02-11) 24 - G4WaterStopping - fixed Ar data 25 - G4PairProductionRelModel - (A.Schaelicke) new relativistic model for 26 gamma conversion 27 - G4UrbanMscModel2 - (L.Urban) new version of theta0 and tail tuning 28 29 10 May 09: V.Ivant (emstand-V09-02-10) 30 - G4WentzelVIModel, G4eCoulombScattering, G4CoulombScattering - added relativistic 31 factor to Reserford cross section; set default limit on 32 kinetic energy of the recoil 100 keV 33 34 28 April 09: V.Ivant (emstand-V09-02-09) 35 - G4UrbanMscModel2 - (L.Urban) new tuning for the central part and for the 36 tail of the angular distribution using the old e- 37 scattering data only (Phys. Rev. 84 (1951) 634; 38 Phys. Rev. 61 (1942) 254) 39 - corrected logic in ComputeTruePathLengthLimit method 40 for type=fUseDistanceToBoundary 41 - G4UrbanMscModel - frozen version of G4UrbanMscModel2 of g4 9.2 42 - G4WentzelVIModel, G4eCoulombScattering - reduce low-limit from 1 keV to 0.1 keV 43 to provide smooth transport cross section table 44 45 23 April 09: V.Ivant (emstand-V09-02-08) 46 - G4BetheBlochModel - do not use pointer to GenericIon introduced in the 47 previous tag due to problem of simple PhysLists without ions 48 49 20 April 09: V.Ivant (emstand-V09-02-07) 50 - G4BetheBlochModel - fixed and simplified initialisation for ions needed for 51 the new G4IonParametrisedLossModel of low-energy package 52 - G4GoudsmitSaundersonMscModel - (O.Kadri) cleanup: discarded no scattering and 53 single scattering theta sampling from SampleCosineTheta() 54 which means the splitting step into two sub-steps occur 55 only for msc regime 56 57 12 April 09: V.Ivant (emstand-V09-02-06) 58 - Simplified initialisation of all models 59 - G4UrbanMscModel, G4UrbanMscModel2, G4UrbanMscModel90 - use methods 60 of G4VMscModel for interface to geometry 61 62 07 April 09: V.Ivant (emstand-V09-02-05) 63 - G4IonFluctuations - T.Toshito removed extra phenomenological factor 64 in fluctuation width 65 - G4HeatedKleinNishinaCompton - V.Grichine added a prototype model for 66 plasma 67 68 21 March 09: V.Ivant (emstand-V09-02-04) 69 - G4UniversalFluctuation - L.Urban introduce modification in width 70 correction, the dependence of the correction on energy deposition 71 at previous steps is removed (addresses T2K report) 72 73 20 March 09: V.Ivant (emstand-V09-02-03) 74 - G4GoudsmitSaundersonMscModel fixed compillation problem 75 76 12 March 09: V.Ivant (emstand-V09-02-02) 77 - G4GoudsmitSaundersonMscModel fixed compillation problem 78 - G4UniversalFluctuation - add temporary fix for T2K report 79 80 05 March 09: V.Ivant (emstand-V09-02-01) 81 - New G4GoudsmitSaundersonMscModel is added 82 - G4WentzelVIModel, G4eCoulombScattringModel: 83 o substitute scaling of low-energy limit by setting 1 keV for 84 all particles; 85 o use EGSnrc form of screening parameter (second order correction) 86 87 20 February 09: V.Ivant (emstand-V09-02-00) 88 - Move all virtual methods from inline to source 89 G4PEEffectModel - substitute ComputeMeanFreePath by CrossSectionPerVolume 90 (minor CPU speadup for compound materials) 91 G4PAIModel, G4PAIPhotonModel - remove usage of random numbers at 92 initialisation (potential non-reproducibility) 93 G4WentzelVIModel - use generic methods of G4VMscModel to access safety 94 and other geometry information 19 95 20 96 24 November 08: V.Ivant (emstand-V09-01-45) -
trunk/source/processes/electromagnetic/standard/include/G4BetheBlochModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BetheBlochModel.hh,v 1. 16 2008/10/22 16:00:57vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4BetheBlochModel.hh,v 1.20 2009/04/23 17:44:43 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 64 64 65 65 #include "G4VEmModel.hh" 66 #include "G4NistManager.hh" 66 67 67 68 class G4EmCorrections; 68 69 class G4ParticleChangeForLoss; 69 class G4NistManager;70 71 70 72 71 class G4BetheBlochModel : public G4VEmModel … … 82 81 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&); 83 82 84 G4double MinEnergyCut(const G4ParticleDefinition*,85 const G4MaterialCutsCouple*);83 virtual G4double MinEnergyCut(const G4ParticleDefinition*, 84 const G4MaterialCutsCouple*); 86 85 87 86 virtual G4double ComputeCrossSectionPerElectron( … … 131 130 protected: 132 131 133 G4double MaxSecondaryEnergy(const G4ParticleDefinition*,134 G4double kinEnergy);132 virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*, 133 G4double kinEnergy); 135 134 136 135 private: 137 136 138 void SetParticle(const G4ParticleDefinition* p); 137 inline void SetupParameters(); 138 139 inline void SetParticle(const G4ParticleDefinition* p); 140 141 inline void SetGenericIon(const G4ParticleDefinition* p); 139 142 140 143 // hide assignment operator … … 166 169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 167 170 168 inline G4double G4BetheBlochModel::MaxSecondaryEnergy( 169 const G4ParticleDefinition* pd, 170 G4double kinEnergy) 171 { 172 if(isIon) SetParticle(pd); 173 G4double tau = kinEnergy/mass; 174 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) / 175 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio); 176 return std::min(tmax,tlimit); 171 inline void G4BetheBlochModel::SetupParameters() 172 { 173 mass = particle->GetPDGMass(); 174 spin = particle->GetPDGSpin(); 175 G4double q = particle->GetPDGCharge()/eplus; 176 chargeSquare = q*q; 177 ratio = electron_mass_c2/mass; 178 G4double magmom = particle->GetPDGMagneticMoment()*mass/(0.5*eplus*hbar_Planck*c_squared); 179 magMoment2 = magmom*magmom - 1.0; 180 formfact = 0.0; 181 if(particle->GetLeptonNumber() == 0) { 182 G4double x = 0.8426*GeV; 183 if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;} 184 else if(mass > GeV) { 185 x /= nist->GetZ13(mass/proton_mass_c2); 186 // tlimit = 51.2*GeV*A13[iz]*A13[iz]; 187 } 188 formfact = 2.0*electron_mass_c2/(x*x); 189 tlimit = 2.0/formfact; 190 } 177 191 } 178 192 179 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 180 194 195 inline void G4BetheBlochModel::SetParticle(const G4ParticleDefinition* p) 196 { 197 if(particle != p) { 198 particle = p; 199 if (p->GetPDGCharge()/eplus > 1.5 && p->GetBaryonNumber() > 2) isIon = true; 200 SetupParameters(); 201 } 202 } 203 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 205 206 inline void G4BetheBlochModel::SetGenericIon(const G4ParticleDefinition* p) 207 { 208 if(p && particle != p) { 209 if(p->GetParticleName() == "GenericIon") isIon = true; 210 } 211 } 212 181 213 #endif -
trunk/source/processes/electromagnetic/standard/include/G4BohrFluctuations.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BohrFluctuations.hh,v 1. 3 2007/09/27 13:53:11vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4BohrFluctuations.hh,v 1.4 2009/02/19 19:17:50 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 79 79 void InitialiseMe(const G4ParticleDefinition*); 80 80 81 protected:82 83 81 private: 84 82 … … 103 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 104 102 105 106 inline G4double G4BohrFluctuations::Dispersion(107 const G4Material* material,108 const G4DynamicParticle* dp,109 G4double& tmax,110 G4double& length)111 {112 if(!particle) InitialiseMe(dp->GetDefinition());113 114 G4double electronDensity = material->GetElectronDensity();115 kineticEnergy = dp->GetKineticEnergy();116 G4double etot = kineticEnergy + particleMass;117 beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);118 G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length119 * electronDensity * chargeSquare;120 121 return siga;122 }123 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....125 126 103 #endif 127 104 -
trunk/source/processes/electromagnetic/standard/include/G4BraggIonModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BraggIonModel.hh,v 1.1 1 2008/10/22 16:00:57 vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4BraggIonModel.hh,v 1.12 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 77 77 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&); 78 78 79 G4double MinEnergyCut(const G4ParticleDefinition*,80 const G4MaterialCutsCouple*);79 virtual G4double MinEnergyCut(const G4ParticleDefinition*, 80 const G4MaterialCutsCouple*); 81 81 82 82 virtual G4double ComputeCrossSectionPerElectron( … … 128 128 protected: 129 129 130 G4double MaxSecondaryEnergy(const G4ParticleDefinition*,131 G4double kinEnergy);130 virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*, 131 G4double kinEnergy); 132 132 133 133 private: … … 180 180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 181 181 182 inline G4double G4BraggIonModel::MaxSecondaryEnergy(183 const G4ParticleDefinition* pd,184 G4double kinEnergy)185 {186 if(pd != particle) SetParticle(pd);187 G4double tau = kinEnergy/mass;188 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /189 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);190 return tmax;191 }192 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......194 195 182 inline void G4BraggIonModel::SetParticle(const G4ParticleDefinition* p) 196 183 { -
trunk/source/processes/electromagnetic/standard/include/G4BraggModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BraggModel.hh,v 1.1 2 2008/09/14 17:11:48vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4BraggModel.hh,v 1.13 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 80 80 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&); 81 81 82 G4double MinEnergyCut(const G4ParticleDefinition*,83 const G4MaterialCutsCouple*);82 virtual G4double MinEnergyCut(const G4ParticleDefinition*, 83 const G4MaterialCutsCouple*); 84 84 85 85 virtual G4double ComputeCrossSectionPerElectron( … … 131 131 protected: 132 132 133 G4double MaxSecondaryEnergy(const G4ParticleDefinition*,134 G4double kinEnergy);133 virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*, 134 G4double kinEnergy); 135 135 136 136 private: 137 137 138 void SetParticle(const G4ParticleDefinition* p);138 inline void SetParticle(const G4ParticleDefinition* p); 139 139 140 140 G4bool HasMaterial(const G4Material* material); … … 182 182 G4bool isInitialised; 183 183 }; 184 185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......186 187 inline G4double G4BraggModel::MaxSecondaryEnergy(188 const G4ParticleDefinition* pd,189 G4double kinEnergy)190 {191 if(pd != particle) SetParticle(pd);192 G4double tau = kinEnergy/mass;193 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /194 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);195 return tmax;196 }197 184 198 185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -
trunk/source/processes/electromagnetic/standard/include/G4ComptonScattering.hh
r1007 r1055 25 25 // 26 26 // 27 // $Id: G4ComptonScattering.hh,v 1.2 0 2007/05/23 08:47:34vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4ComptonScattering.hh,v 1.21 2009/02/20 12:06:37 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 29 29 // 30 30 //------------------ G4ComptonScattering physics process ----------------------- … … 82 82 83 83 // true for Gamma only. 84 G4bool IsApplicable(const G4ParticleDefinition&);84 virtual G4bool IsApplicable(const G4ParticleDefinition&); 85 85 86 86 // Print few lines of informations about the process: validity range, 87 87 virtual void PrintInfo(); 88 89 88 90 89 protected: … … 97 96 }; 98 97 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....101 102 inline103 G4bool G4ComptonScattering::IsApplicable(const G4ParticleDefinition& p)104 {105 return (&p == G4Gamma::Gamma());106 }107 108 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 109 99 -
trunk/source/processes/electromagnetic/standard/include/G4CoulombScattering.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4CoulombScattering.hh,v 1.1 1 2008/06/13 08:19:43vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4CoulombScattering.hh,v 1.13 2009/05/07 18:41:45 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 59 59 public: 60 60 61 G4CoulombScattering(const G4String& name = " eCoulombScat");61 G4CoulombScattering(const G4String& name = "CoulombScat"); 62 62 63 63 virtual ~G4CoulombScattering(); … … 104 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 105 105 106 inline G4bool G4CoulombScattering::IsApplicable(const G4ParticleDefinition& p)107 {108 return (p.GetPDGCharge() != 0.0 && !p.IsShortLived());109 }110 111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....112 113 106 inline void G4CoulombScattering::SetThetaMin(G4double val) 114 107 { -
trunk/source/processes/electromagnetic/standard/include/G4GammaConversion.hh
r1007 r1055 25 25 // 26 26 // 27 // $Id: G4GammaConversion.hh,v 1.2 2 2007/05/23 08:47:34vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4GammaConversion.hh,v 1.23 2009/02/20 12:06:37 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 29 29 // 30 30 // … … 85 85 86 86 // true for Gamma only. 87 G4bool IsApplicable(const G4ParticleDefinition&);87 virtual G4bool IsApplicable(const G4ParticleDefinition&); 88 88 89 89 // Print few lines of informations about the process: validity range, … … 99 99 }; 100 100 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....103 104 inline G4bool G4GammaConversion::IsApplicable(const G4ParticleDefinition& p)105 {106 return (&p == G4Gamma::Gamma());107 }108 109 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 110 102 -
trunk/source/processes/electromagnetic/standard/include/G4IonFluctuations.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4IonFluctuations.hh,v 1. 8 2008/10/22 16:04:33vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4IonFluctuations.hh,v 1.9 2009/02/19 19:17:50 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 84 84 85 85 // Initialisation prestep 86 inlinevoid SetParticleAndCharge(const G4ParticleDefinition*, G4double q2);86 void SetParticleAndCharge(const G4ParticleDefinition*, G4double q2); 87 87 88 88 private: … … 117 117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 118 118 119 inline120 void G4IonFluctuations::SetParticleAndCharge(const G4ParticleDefinition* part,121 G4double q2)122 {123 if(part != particle) {124 particle = part;125 particleMass = part->GetPDGMass();126 charge = part->GetPDGCharge()/eplus;127 chargeSquare = charge*charge;128 }129 effChargeSquare = q2;130 uniFluct.SetParticleAndCharge(part, q2);131 }132 133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....134 135 119 #endif 136 120 -
trunk/source/processes/electromagnetic/standard/include/G4MollerBhabhaModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4MollerBhabhaModel.hh,v 1. 19 2007/05/22 17:34:36vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4MollerBhabhaModel.hh,v 1.20 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 74 74 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&); 75 75 76 G4double MinEnergyCut(const G4ParticleDefinition*,77 const G4MaterialCutsCouple*);76 virtual G4double MinEnergyCut(const G4ParticleDefinition*, 77 const G4MaterialCutsCouple*); 78 78 79 79 virtual G4double ComputeCrossSectionPerElectron( … … 109 109 protected: 110 110 111 G4double MaxSecondaryEnergy(const G4ParticleDefinition*,112 G4double kinEnergy);111 virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*, 112 G4double kinEnergy); 113 113 114 void SetParticle(const G4ParticleDefinition* p);114 inline void SetParticle(const G4ParticleDefinition* p); 115 115 116 116 const G4ParticleDefinition* particle; … … 127 127 G4MollerBhabhaModel & operator=(const G4MollerBhabhaModel &right); 128 128 G4MollerBhabhaModel(const G4MollerBhabhaModel&); 129 130 G4bool isInitialised; 131 129 132 }; 130 133 131 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 132 135 133 inline G4double G4MollerBhabhaModel::MaxSecondaryEnergy( 134 const G4ParticleDefinition*, 135 G4double kinEnergy) 136 inline void G4MollerBhabhaModel::SetParticle(const G4ParticleDefinition* p) 136 137 { 137 G4double tmax = kinEnergy; 138 if(isElectron) tmax *= 0.5; 139 return tmax; 138 particle = p; 139 if(p != theElectron) isElectron = false; 140 140 } 141 141 -
trunk/source/processes/electromagnetic/standard/include/G4PAIModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PAIModel.hh,v 1.22 2009/02/19 19:17:50 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 26 28 // 27 29 // ------------------------------------------------------------------- … … 75 77 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&); 76 78 77 virtual void InitialiseMe(const G4ParticleDefinition*) {};78 79 virtual G4double ComputeDEDX (const G4MaterialCutsCouple*,79 virtual void InitialiseMe(const G4ParticleDefinition*); 80 81 virtual G4double ComputeDEDXPerVolume(const G4Material*, 80 82 const G4ParticleDefinition*, 81 83 G4double kineticEnergy, 82 84 G4double cutEnergy); 83 85 84 virtual G4double CrossSection (const G4MaterialCutsCouple*,86 virtual G4double CrossSectionPerVolume(const G4Material*, 85 87 const G4ParticleDefinition*, 86 88 G4double kineticEnergy, … … 118 120 119 121 void SetVerboseLevel(G4int verbose){fVerbose=verbose;}; 120 121 122 122 123 123 protected: … … 192 192 }; 193 193 194 /////////////////////////////////////////////////////////////////////195 196 inline G4double G4PAIModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,197 G4double kinEnergy)198 {199 G4double tmax = kinEnergy;200 if(p == fElectron) tmax *= 0.5;201 else if(p != fPositron) {202 G4double mass = p->GetPDGMass();203 G4double ratio= electron_mass_c2/mass;204 G4double gamma= kinEnergy/mass + 1.0;205 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /206 (1. + 2.0*gamma*ratio + ratio*ratio);207 }208 return tmax;209 }210 211 ///////////////////////////////////////////////////////////////212 213 inline void G4PAIModel::DefineForRegion(const G4Region* r)214 {215 fPAIRegionVector.push_back(r);216 }217 218 194 #endif 219 195 -
trunk/source/processes/electromagnetic/standard/include/G4PAIPhotonModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PAIPhotonModel.hh,v 1.12 2009/02/19 19:17:50 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 26 28 // 27 29 // ------------------------------------------------------------------- … … 74 76 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&); 75 77 76 virtual void InitialiseMe(const G4ParticleDefinition*) {};77 78 virtual G4double ComputeDEDX (const G4MaterialCutsCouple*,79 80 81 82 83 virtual G4double CrossSection (const G4MaterialCutsCouple*,84 const G4ParticleDefinition*,85 G4double kineticEnergy,86 G4double cutEnergy,87 G4double maxEnergy);78 virtual void InitialiseMe(const G4ParticleDefinition*); 79 80 virtual G4double ComputeDEDXPerVolume(const G4Material*, 81 const G4ParticleDefinition*, 82 G4double kineticEnergy, 83 G4double cutEnergy); 84 85 virtual G4double CrossSectionPerVolume(const G4Material*, 86 const G4ParticleDefinition*, 87 G4double kineticEnergy, 88 G4double cutEnergy, 89 G4double maxEnergy); 88 90 89 91 virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*, … … 122 124 G4double position, G4int iTransfer ); 123 125 124 125 126 126 protected: 127 127 128 128 G4double MaxSecondaryEnergy(const G4ParticleDefinition*, 129 129 G4double kinEnergy); 130 130 131 131 private: … … 145 145 G4int fVerbose; 146 146 G4PhysicsLogVector* fProtonEnergyVector ; 147 148 149 147 150 148 // vectors … … 204 202 }; 205 203 206 /////////////////////////////////////////////////////////////////////207 208 inline G4double G4PAIPhotonModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,209 G4double kinEnergy)210 {211 G4double tmax = kinEnergy;212 if(p == fElectron) tmax *= 0.5;213 else if(p != fPositron) {214 G4double mass = p->GetPDGMass();215 G4double ratio= electron_mass_c2/mass;216 G4double gamma= kinEnergy/mass + 1.0;217 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /218 (1. + 2.0*gamma*ratio + ratio*ratio);219 }220 return tmax;221 }222 223 ///////////////////////////////////////////////////////////////224 225 inline void G4PAIPhotonModel::DefineForRegion(const G4Region* r)226 {227 fPAIRegionVector.push_back(r);228 }229 230 204 #endif 231 205 -
trunk/source/processes/electromagnetic/standard/include/G4PEEffectModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PEEffectModel.hh,v 1. 6 2007/05/22 17:34:36vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4PEEffectModel.hh,v 1.7 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 41 41 // 42 42 // 06.02.2006 : added ComputeMeanFreePath() (mma) 43 // 20.02.2009 : move virtual inline to .cc, substitute 44 // ComputeMeanFreePath() by CrossSectionPerVolume (VI) 43 45 // 44 46 // Class Description: … … 70 72 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&); 71 73 72 G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,73 74 75 76 74 virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 75 G4double kinEnergy, 76 G4double Z, 77 G4double A, 78 G4double, G4double); 77 79 78 G4double ComputeMeanFreePath( const G4ParticleDefinition*, 79 G4double kinEnergy, 80 const G4Material* material, 81 G4double, G4double); 80 virtual G4double CrossSectionPerVolume(const G4Material*, 81 const G4ParticleDefinition*, 82 G4double kineticEnergy, 83 G4double cutEnergy, 84 G4double maxEnergy); 82 85 83 86 virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*, … … 100 103 }; 101 104 102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....103 104 inline G4double G4PEEffectModel::ComputeCrossSectionPerAtom(105 const G4ParticleDefinition*,106 G4double energy,107 G4double Z, G4double,108 G4double, G4double)109 {110 G4double* SandiaCof = G4SandiaTable::GetSandiaCofPerAtom((G4int)Z, energy);111 112 G4double energy2 = energy*energy, energy3 = energy*energy2,113 energy4 = energy2*energy2;114 115 return SandiaCof[0]/energy + SandiaCof[1]/energy2 +116 SandiaCof[2]/energy3 + SandiaCof[3]/energy4;117 }118 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....120 121 inline G4double G4PEEffectModel::ComputeMeanFreePath(122 const G4ParticleDefinition*,123 G4double energy,124 const G4Material* material,125 G4double, G4double)126 {127 G4double* SandiaCof = material->GetSandiaTable()128 ->GetSandiaCofForMaterial(energy);129 130 G4double energy2 = energy*energy, energy3 = energy*energy2,131 energy4 = energy2*energy2;132 133 G4double cross = SandiaCof[0]/energy + SandiaCof[1]/energy2 +134 SandiaCof[2]/energy3 + SandiaCof[3]/energy4;135 136 G4double mfp = DBL_MAX;137 if (cross > 0.) mfp = 1./cross;138 return mfp;139 }140 141 105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 142 106 -
trunk/source/processes/electromagnetic/standard/include/G4PhotoElectricEffect.hh
r1007 r1055 25 25 // 26 26 // 27 // $Id: G4PhotoElectricEffect.hh,v 1.2 4 2007/05/23 08:47:34vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4PhotoElectricEffect.hh,v 1.25 2009/02/20 12:06:37 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 29 29 // 30 30 // … … 91 91 92 92 // true for Gamma only. 93 G4bool IsApplicable(const G4ParticleDefinition&);93 virtual G4bool IsApplicable(const G4ParticleDefinition&); 94 94 95 95 // Print few lines of informations about the process: validity range, 96 v oid PrintInfo();96 virtual void PrintInfo(); 97 97 98 98 protected: 99 99 100 v oid InitialiseProcess(const G4ParticleDefinition*);100 virtual void InitialiseProcess(const G4ParticleDefinition*); 101 101 102 102 private: … … 105 105 }; 106 106 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....109 110 inline111 G4bool G4PhotoElectricEffect::IsApplicable(const G4ParticleDefinition& p)112 {113 return (&p == G4Gamma::Gamma());114 }115 116 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 117 108 -
trunk/source/processes/electromagnetic/standard/include/G4UniversalFluctuation.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4UniversalFluctuation.hh,v 1. 6 2008/10/22 16:04:33vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4UniversalFluctuation.hh,v 1.10 2009/03/19 14:15:17 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 70 70 virtual ~G4UniversalFluctuation(); 71 71 72 G4double SampleFluctuations(const G4Material*,73 74 75 76 72 virtual G4double SampleFluctuations(const G4Material*, 73 const G4DynamicParticle*, 74 G4double&, 75 G4double&, 76 G4double&); 77 77 78 G4double Dispersion( const G4Material*,79 80 81 78 virtual G4double Dispersion( const G4Material*, 79 const G4DynamicParticle*, 80 G4double&, 81 G4double&); 82 82 83 v oid InitialiseMe(const G4ParticleDefinition*);83 virtual void InitialiseMe(const G4ParticleDefinition*); 84 84 85 85 // Initialisation prestep 86 inline void SetParticleAndCharge(const G4ParticleDefinition*, G4double q2); 87 88 protected: 86 virtual void SetParticleAndCharge(const G4ParticleDefinition*, G4double q2); 89 87 90 88 private: … … 113 111 G4double e0; 114 112 115 G4double facwidth;116 G4double oldloss;117 G4double samestep;118 113 G4double e1,e2; 119 114 … … 126 121 }; 127 122 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....129 130 inline void131 G4UniversalFluctuation::SetParticleAndCharge(const G4ParticleDefinition* part,132 G4double q2)133 {134 if(part != particle) {135 particle = part;136 particleMass = part->GetPDGMass();137 }138 chargeSquare = q2;139 }140 141 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 142 124 -
trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4UrbanMscModel.hh,v 1.3 3 2008/03/10 10:39:21vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4UrbanMscModel.hh,v 1.35 2009/04/29 13:30:22 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 37 37 // Author: Laszlo Urban 38 38 // 39 // Creation date: 0 1.03.200139 // Creation date: 06.03.2008 40 40 // 41 41 // Modifications: 42 42 // 43 // 27-03-03 Move model part from G4MultipleScattering (V.Ivanchenko) 44 // 27-03-03 Rename (V.Ivanchenko) 45 // 46 // 05-08-03 angle distribution has been modified (L.Urban) 47 // 26-11-03 new data member currentRange (L.Urban) 48 // 01-03-04 changes in data members + signature changed in SampleCosineTheta 49 // 11-03-04 changes in data members (L.Urban) 50 // 23-04-04 changes in data members and in signature of SampleCosineTheta 51 // (L.Urban) 52 // 17-08-04 name of data member facxsi changed to factail (L.Urban) 53 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko) 54 // 15-04-05 optimize internal interface - add SampleSecondaries method (V.Ivanchenko) 55 // 11-08-05 computation of lateral correlation added (L.Urban) 56 // 02-10-05 nuclear size correction computation removed, the correction 57 // included in the (theoretical) tabulated values (L.Urban) 58 // 16-02-06 data members b and xsi have been removed (L.Urban) 59 // 17-02-06 Save table of transport cross sections not mfp (V.Ivanchenko) 60 // 07-03-06 Create G4UrbanMscModel and move there step limit calculation (V.Ivanchenko) 61 // 10-05-06 SetMscStepLimitation at initialisation (V.Ivantchenko) 62 // 11-05-06 name of data member safety changed to presafety, some new data 63 // members added (frscaling1,frscaling2,tlimitminfix,nstepmax) 64 // (L.Urban) 65 // 13-10-06 data member factail removed, data member tkinlimit changed 66 // to lambdalimit, 67 // new data members tgeom,tnow,skin,skindepth,Zeff,geomlimit 68 // G4double GeomLimit(const G4Track& track) changed to 69 // void GeomLimit(const G4Track& track) (L.Urban) 70 // 20-10-06 parameter theta0 now computed in the (public) 71 // function ComputeTheta0, 72 // single scattering modified allowing not small 73 // angles as well (L.Urban) 74 // 31-01-07 code cleaning (L.Urban) 75 // 06-02-07 Move SetMscStepLimitation method into the source (VI) 76 // 15-02-07 new data member : smallstep (L.Urban) 77 // 10-04-07 remove navigator, smallstep, tnow (V.Ivanchenko) 78 // 43 // 28-04-09 move G4UrbanMscModel2 from the g49.2 to G4UrbanMscModel. 44 // now it is frozen (V.Ivanchenk0) 79 45 // 80 46 // Class Description: … … 111 77 112 78 void Initialise(const G4ParticleDefinition*, const G4DataVector&); 113 79 114 80 G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition* particle, 115 81 G4double KineticEnergy, … … 119 85 G4double emax=DBL_MAX); 120 86 121 void SampleSecondaries(std::vector<G4DynamicParticle*>*,122 const G4MaterialCutsCouple*,123 const G4DynamicParticle*,124 G4double,125 G4double);126 127 87 void SampleScattering(const G4DynamicParticle*, 128 88 G4double safety); … … 141 101 private: 142 102 103 G4double SimpleScattering(G4double xmeanth, G4double x2meanth); 104 143 105 G4double SampleCosineTheta(G4double trueStepLength, G4double KineticEnergy); 144 106 … … 147 109 G4double LatCorrelation(); 148 110 149 void GeomLimit(const G4Track& track);150 151 111 inline G4double GetLambda(G4double kinEnergy); 152 112 153 113 inline void SetParticle(const G4ParticleDefinition*); 114 115 inline void UpdateCache(); 154 116 155 117 // hide assignment operator … … 160 122 G4ParticleChangeForMSC* fParticleChange; 161 123 162 G4SafetyHelper* safetyHelper;163 124 G4PhysicsTable* theLambdaTable; 164 125 const G4MaterialCutsCouple* couple; … … 166 127 167 128 G4double mass; 168 G4double charge; 169 170 G4double masslimite,masslimitmu; 129 G4double charge,ChargeSquare; 130 G4double masslimite,lambdalimit,fr; 171 131 172 132 G4double taubig; … … 174 134 G4double taulim; 175 135 G4double currentTau; 176 G4double frscaling1,frscaling2;177 136 G4double tlimit; 178 137 G4double tlimitmin; 179 138 G4double tlimitminfix; 180 181 G4double nstepmax; 139 G4double tgeom; 140 182 141 G4double geombig; 183 142 G4double geommin; … … 198 157 G4double currentKinEnergy; 199 158 G4double currentRange; 159 G4double rangeinit; 200 160 G4double currentRadLength; 201 161 202 G4double Zeff; 162 G4double theta0max,rellossmax; 163 G4double third; 203 164 204 165 G4int currentMaterialIndex; 166 167 G4double y; 168 G4double Zold; 169 G4double Zeff,Z2,Z23,lnZ; 170 G4double coeffth1,coeffth2; 171 G4double coeffc1,coeffc2; 172 G4double scr1ini,scr2ini,scr1,scr2; 205 173 206 174 G4bool isInitialized; 207 175 G4bool inside; 208 176 G4bool insideskin; 209 210 177 }; 211 178 179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 212 180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 213 181 … … 236 204 mass = p->GetPDGMass(); 237 205 charge = p->GetPDGCharge()/eplus; 206 ChargeSquare = charge*charge; 238 207 } 239 208 } … … 241 210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 242 211 212 inline 213 void G4UrbanMscModel::UpdateCache() 214 { 215 lnZ = std::log(Zeff); 216 coeffth1 = 0.885+lnZ*(0.104-0.0170*lnZ); 217 coeffth2 = 0.028+lnZ*(0.012-0.00125*lnZ); 218 coeffc1 = 2.134-lnZ*(0.1045-0.00602*lnZ); 219 coeffc2 = 0.001126-lnZ*(0.0001089+0.0000247*lnZ); 220 Z2 = Zeff*Zeff; 221 Z23 = std::exp(2.*lnZ/3.); 222 scr1 = scr1ini*Z23; 223 scr2 = scr2ini*Z2*ChargeSquare; 224 // lastMaterial = couple->GetMaterial(); 225 Zold = Zeff; 226 } 227 243 228 #endif 244 229 -
trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel2.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4UrbanMscModel2.hh,v 1.1 1 2008/12/18 13:01:34 gunterExp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4UrbanMscModel2.hh,v 1.17 2009/05/15 09:26:42 urban Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 40 40 // 41 41 // Modifications: 42 // 42 // 23.04.2009 L.Urban updated parameterization in UpdateCache method 43 43 // 44 44 // … … 84 84 G4double emax=DBL_MAX); 85 85 86 void SampleSecondaries(std::vector<G4DynamicParticle*>*,87 const G4MaterialCutsCouple*,88 const G4DynamicParticle*,89 G4double,90 G4double);91 92 86 void SampleScattering(const G4DynamicParticle*, 93 87 G4double safety); … … 114 108 G4double LatCorrelation(); 115 109 116 void GeomLimit(const G4Track& track);117 118 110 inline G4double GetLambda(G4double kinEnergy); 119 111 … … 129 121 G4ParticleChangeForMSC* fParticleChange; 130 122 131 G4SafetyHelper* safetyHelper;132 123 G4PhysicsTable* theLambdaTable; 133 124 const G4MaterialCutsCouple* couple; 134 125 G4LossTableManager* theManager; 135 136 126 137 127 G4double mass; … … 223 213 { 224 214 lnZ = std::log(Zeff); 225 coeffth1 = 0.885+lnZ*(0.104-0.0170*lnZ); 226 coeffth2 = 0.028+lnZ*(0.012-0.00125*lnZ); 227 coeffc1 = 2.134-lnZ*(0.1045-0.00602*lnZ); 228 coeffc2 = 0.001126-lnZ*(0.0001089+0.0000247*lnZ); 215 //new correction in theta0 formula 216 coeffth1 = (1.-8.7780e-2/Zeff)*(0.87+0.03*lnZ); 217 coeffth2 = (4.0780e-2+1.7315e-4*Zeff)*(0.87+0.03*lnZ); 218 // tail parameters 219 G4double lnZ1 = std::log(Zeff+1.); 220 coeffc1 = 2.943-0.197*lnZ1; 221 coeffc2 = 0.0987-0.0143*lnZ1; 222 // for single scattering 229 223 Z2 = Zeff*Zeff; 230 224 Z23 = std::exp(2.*lnZ/3.); 231 225 scr1 = scr1ini*Z23; 232 226 scr2 = scr2ini*Z2*ChargeSquare; 233 // lastMaterial = couple->GetMaterial();227 234 228 Zold = Zeff; 235 229 } -
trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel90.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4UrbanMscModel90.hh,v 1. 4 2008/10/29 14:15:30vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4UrbanMscModel90.hh,v 1.5 2009/04/10 16:34:56 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 110 110 111 111 G4double LatCorrelation(); 112 113 void GeomLimit(const G4Track& track);114 112 115 113 inline G4double GetLambda(G4double kinEnergy); -
trunk/source/processes/electromagnetic/standard/include/G4WentzelVIModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4WentzelVIModel.hh,v 1. 7 2008/08/04 08:49:09 vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4WentzelVIModel.hh,v 1.17 2009/05/10 16:09:29 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 66 66 class G4LossTableManager; 67 67 class G4ParticleChangeForMSC; 68 class G4SafetyHelper;69 68 class G4ParticleDefinition; 70 69 … … 80 79 virtual ~G4WentzelVIModel(); 81 80 82 void Initialise(const G4ParticleDefinition*, const G4DataVector&); 83 84 G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 85 G4double KineticEnergy, 86 G4double AtomicNumber, 87 G4double AtomicWeight=0., 88 G4double cut = DBL_MAX, 89 G4double emax= DBL_MAX); 90 91 void SampleScattering(const G4DynamicParticle*, G4double safety); 92 93 void SampleSecondaries(std::vector<G4DynamicParticle*>*, 94 const G4MaterialCutsCouple*, 95 const G4DynamicParticle*, 96 G4double, 97 G4double); 98 99 G4double ComputeTruePathLengthLimit(const G4Track& track, 100 G4PhysicsTable* theLambdaTable, 101 G4double currentMinimalStep); 102 103 G4double ComputeGeomPathLength(G4double truePathLength); 104 105 G4double ComputeTrueStepLength(G4double geomStepLength); 81 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&); 82 83 virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 84 G4double KineticEnergy, 85 G4double AtomicNumber, 86 G4double AtomicWeight=0., 87 G4double cut = DBL_MAX, 88 G4double emax= DBL_MAX); 89 90 virtual void SampleScattering(const G4DynamicParticle*, G4double safety); 91 92 virtual G4double ComputeTruePathLengthLimit(const G4Track& track, 93 G4PhysicsTable* theLambdaTable, 94 G4double currentMinimalStep); 95 96 virtual G4double ComputeGeomPathLength(G4double truePathLength); 97 98 virtual G4double ComputeTrueStepLength(G4double geomStepLength); 106 99 107 100 private: 108 101 109 G4double ComputeTransportXSectionPer Volume();102 G4double ComputeTransportXSectionPerAtom(); 110 103 111 104 G4double ComputeXSectionPerVolume(); … … 133 126 G4ParticleChangeForMSC* fParticleChange; 134 127 135 G4SafetyHelper* safetyHelper;136 128 G4PhysicsTable* theLambdaTable; 137 129 G4PhysicsTable* theLambda2Table; … … 173 165 // single scattering parameters 174 166 G4double coeff; 175 G4double constn;176 167 G4double cosThetaMin; 177 168 G4double cosThetaMax; … … 182 173 G4double q2Limit; 183 174 G4double alpha2; 184 G4double a0;185 175 186 176 // projectile … … 193 183 G4double mom2; 194 184 G4double invbeta2; 185 G4double kinFactor; 195 186 G4double etag; 196 187 G4double lowEnergyLimit; … … 198 189 // target 199 190 G4double targetZ; 191 G4double targetMass; 200 192 G4double screenZ; 201 193 G4double formfactA; 202 G4double FF[100]; 194 G4int iz; 195 196 static G4double ScreenRSquare[100]; 197 static G4double FormFactor[100]; 203 198 204 199 // flags … … 240 235 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 241 236 242 inline 243 void G4WentzelVIModel::SetupParticle(const G4ParticleDefinition* p) 237 inline void G4WentzelVIModel::SetupParticle(const G4ParticleDefinition* p) 244 238 { 245 239 // Initialise mass and charge … … 251 245 chargeSquare = q*q; 252 246 tkin = 0.0; 253 lowEnergyLimit = keV*mass/electron_mass_c2;254 247 } 255 248 } … … 264 257 invbeta2 = 1.0 + mass*mass/mom2; 265 258 cosTetMaxNuc = cosThetaMax; 266 if( ekin <= 10.*cut && mass < MeV) {259 if(mass < MeV && ekin <= 10.*cut) { 267 260 cosTetMaxNuc = ekin*(cosThetaMax + 1.0)/(10.*cut) - 1.0; 268 261 } … … 278 271 etag = e; 279 272 targetZ = Z; 280 G4int iz= G4int(Z);273 iz = G4int(Z); 281 274 if(iz > 99) iz = 99; 282 G4double x = fNistManager->GetZ13(iz); 283 screenZ = a0*x*x/mom2; 284 if(iz > 1) screenZ *=(1.13 + 3.76*invbeta2*Z*Z*chargeSquare*alpha2); 285 // screenZ = a0*x*x*(1.13 + 3.76*Z*Z*chargeSquare*alpha2)/mom2; 286 // A.V. Butkevich et al., NIM A 488 (2002) 282 287 formfactA = FF[iz]; 288 if(formfactA == 0.0) { 289 x = fNistManager->GetA27(iz); 290 formfactA = constn*x*x; 291 FF[iz] = formfactA; 275 targetMass = fNistManager->GetAtomicMassAmu(iz)*amu_c2; 276 G4double m12 = mass*mass; 277 G4double x = 1.0 + mass/targetMass; 278 kinFactor = coeff*targetZ*chargeSquare*(1.0 + m12/mom2)/mom2; 279 screenZ = ScreenRSquare[iz]/mom2; 280 if(iz > 1) { 281 screenZ *=(1.13 + 3.76*Z*Z*alpha2); 282 kinFactor /= (x*x); 292 283 } 293 formfactA *= mom2; 284 //if(iz > 1) screenZ *=(1.13 + std::min(0.5,3.76*Z*Z*invbeta2*alpha2)); 285 formfactA = FormFactor[iz]*mom2; 294 286 cosTetMaxNuc2 = cosTetMaxNuc; 295 /*296 G4double ee = 10.*eV*Z;297 if(1 == iz) ee *= 2.0;298 G4double z = std::min(cosTetMaxElec, 1.0 - std::max(ecut,ee)*amu_c2299 *fNistManager->GetAtomicMassAmu(iz)/mom2);300 cosTetMaxElec2 = std::max(cosTetMaxNuc2, z);301 */302 287 cosTetMaxElec2 = cosTetMaxElec; 303 288 } -
trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlung.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eBremsstrahlung.hh,v 1.3 6 2007/05/23 08:47:34vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eBremsstrahlung.hh,v 1.37 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 89 89 virtual ~G4eBremsstrahlung(); 90 90 91 G4bool IsApplicable(const G4ParticleDefinition& p);91 virtual G4bool IsApplicable(const G4ParticleDefinition& p); 92 92 93 93 // Print out of the class parameters … … 111 111 112 112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....114 115 inline G4bool G4eBremsstrahlung::IsApplicable(const G4ParticleDefinition& p)116 {117 return (&p == G4Electron::Electron() || &p == G4Positron::Positron());118 }119 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....121 113 122 114 #endif -
trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlungModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eBremsstrahlungModel.hh,v 1.2 5 2008/11/13 19:28:58vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eBremsstrahlungModel.hh,v 1.26 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 75 75 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&); 76 76 77 G4double MinEnergyCut(const G4ParticleDefinition*,78 const G4MaterialCutsCouple*);77 virtual G4double MinEnergyCut(const G4ParticleDefinition*, 78 const G4MaterialCutsCouple*); 79 79 80 80 virtual G4double ComputeDEDXPerVolume(const G4Material*, … … 102 102 103 103 protected: 104 105 inline G4double MaxSecondaryEnergy(const G4ParticleDefinition*,106 G4double kineticEnergy);107 104 108 105 const G4Element* SelectRandomAtom(const G4MaterialCutsCouple* couple); … … 191 188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 192 189 193 inline194 G4double G4eBremsstrahlungModel::MaxSecondaryEnergy(195 const G4ParticleDefinition*,196 G4double kineticEnergy)197 {198 return kineticEnergy;199 }200 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....202 203 190 #endif -
trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlungRelModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eBremsstrahlungRelModel.hh,v 1.1 0 2008/11/14 09:25:19vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eBremsstrahlungRelModel.hh,v 1.11 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 71 71 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&); 72 72 73 G4double MinEnergyCut(const G4ParticleDefinition*,74 const G4MaterialCutsCouple*);73 virtual G4double MinEnergyCut(const G4ParticleDefinition*, 74 const G4MaterialCutsCouple*); 75 75 76 76 virtual G4double ComputeDEDXPerVolume(const G4Material*, … … 97 97 inline G4double LPMconstant() const; 98 98 99 protected:100 101 inline G4double MaxSecondaryEnergy(const G4ParticleDefinition*,102 G4double kineticEnergy);103 104 99 private: 105 100 … … 207 202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 208 203 209 inline G4double210 G4eBremsstrahlungRelModel::MaxSecondaryEnergy(const G4ParticleDefinition*,211 G4double kineticEnergy)212 {213 return kineticEnergy;214 }215 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....217 218 204 219 205 inline G4double G4eBremsstrahlungRelModel::Phi1(G4double gg, G4double) -
trunk/source/processes/electromagnetic/standard/include/G4eCoulombScatteringModel.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eCoulombScatteringModel.hh,v 1. 36 2008/08/04 08:49:09 vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eCoulombScatteringModel.hh,v 1.43 2009/05/10 16:09:29 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 65 65 #include "G4VEmModel.hh" 66 66 #include "G4PhysicsTable.hh" 67 #include "globals.hh" 67 68 #include "G4NistManager.hh" 68 #include "globals.hh"69 69 70 70 class G4ParticleChangeForGamma; … … 137 137 138 138 G4double coeff; 139 G4double constn;140 139 G4double cosThetaMin; 141 140 G4double cosThetaMax; … … 160 159 G4double mom2; 161 160 G4double invbeta2; 161 G4double kinFactor; 162 162 G4double etag; 163 163 G4double lowEnergyLimit; … … 165 165 // target 166 166 G4double targetZ; 167 G4double targetMass; 167 168 G4double screenZ; 168 169 G4double formfactA; 169 170 G4int idxelm; 171 G4int iz; 170 172 171 173 private: 172 174 173 G4double a0;174 175 G4double alpha2; 175 176 G4double faclim; 176 G4double FF[100]; 177 178 static G4double ScreenRSquare[100]; 179 static G4double FormFactor[100]; 177 180 178 181 G4bool isInitialised; … … 204 207 chargeSquare = q*q; 205 208 tkin = 0.0; 206 lowEnergyLimit = keV*mass/electron_mass_c2;207 209 } 208 210 } … … 219 221 cosTetMinNuc = cosThetaMin; 220 222 cosTetMaxNuc = cosThetaMax; 221 if( ekin <= 10.*cut && mass < MeV && cosThetaMin < 1.0) {223 if(mass < MeV && cosThetaMin < 1.0 && ekin <= 10.*cut) { 222 224 cosTetMinNuc = ekin*(cosThetaMin + 1.0)/(10.*cut) - 1.0; 223 225 } … … 233 235 etag = e; 234 236 targetZ = Z; 235 G4intiz= G4int(Z);237 iz= G4int(Z); 236 238 if(iz > 99) iz = 99; 237 G4double x = fNistManager->GetZ13(iz); 238 screenZ = a0*x*x/mom2; 239 if(iz > 1) screenZ *=(1.13 + 3.76*invbeta2*Z*Z*chargeSquare*alpha2); 240 //screenZ = a0*x*x*(1.13 + 3.76*Z*Z*chargeSquare*alpha2)/mom2; 241 // A.V. Butkevich et al., NIM A 488 (2002) 282 242 formfactA = FF[iz]; 243 if(formfactA == 0.0) { 244 x = fNistManager->GetA27(iz); 245 formfactA = constn*x*x; 246 FF[iz] = formfactA; 239 targetMass = fNistManager->GetAtomicMassAmu(iz)*amu_c2; 240 G4double m12 = mass*mass; 241 G4double x = 1.0 + mass/targetMass; 242 kinFactor = (1.0 + m12/mom2)/mom2; 243 244 screenZ = ScreenRSquare[iz]/mom2; 245 if(iz > 1) { 246 screenZ *=(1.13 + 3.76*Z*Z*alpha2); 247 kinFactor /= (x*x); 247 248 } 248 formfactA *= mom2; 249 // if(iz > 1) screenZ *=(1.13 + std::min(0.5,3.76*Z*Z*invbeta2*alpha2)); 250 formfactA = FormFactor[iz]*mom2; 249 251 cosTetMaxNuc2 = cosTetMaxNuc; 250 if( particle == theProton && 1 == iz&& cosTetMaxNuc2 < 0.0) {252 if(1 == iz && particle == theProton && cosTetMaxNuc2 < 0.0) { 251 253 cosTetMaxNuc2 = 0.0; 252 254 } 253 /*254 G4double ee = 10.*eV*Z;255 if(1 == iz) ee *= 2.0;256 G4double z = std::min(cosTetMaxElec, 1.0 - std::max(ecut,ee)*amu_c2257 *fNistManager->GetAtomicMassAmu(iz)/mom2);258 cosTetMaxElec2 = std::max(cosTetMaxNuc2, z);259 */260 255 cosTetMaxElec2 = cosTetMaxElec; 261 256 } -
trunk/source/processes/electromagnetic/standard/include/G4eIonisation.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eIonisation.hh,v 1.3 5 2007/05/23 08:47:34vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eIonisation.hh,v 1.36 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 88 88 virtual ~G4eIonisation(); 89 89 90 G4bool IsApplicable(const G4ParticleDefinition& p);90 virtual G4bool IsApplicable(const G4ParticleDefinition& p); 91 91 92 92 // Print out of the class parameters … … 98 98 const G4ParticleDefinition*); 99 99 100 G4double MinPrimaryEnergy(const G4ParticleDefinition*,101 const G4Material*, G4double cut);100 virtual G4double MinPrimaryEnergy(const G4ParticleDefinition*, 101 const G4Material*, G4double cut); 102 102 103 103 private: … … 114 114 }; 115 115 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....118 119 inline G4double G4eIonisation::MinPrimaryEnergy(const G4ParticleDefinition*,120 const G4Material*,121 G4double cut)122 {123 G4double x = cut;124 if(isElectron) x += cut;125 return x;126 }127 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....129 130 inline G4bool G4eIonisation::IsApplicable(const G4ParticleDefinition& p)131 {132 return (&p == G4Electron::Electron() || &p == G4Positron::Positron());133 }134 135 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 136 117 -
trunk/source/processes/electromagnetic/standard/include/G4eplusAnnihilation.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eplusAnnihilation.hh,v 1.2 3 2007/05/23 08:47:34vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eplusAnnihilation.hh,v 1.24 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 73 73 virtual ~G4eplusAnnihilation(); 74 74 75 G4bool IsApplicable(const G4ParticleDefinition& p);75 virtual G4bool IsApplicable(const G4ParticleDefinition& p); 76 76 77 77 virtual G4VParticleChange* AtRestDoIt( … … 79 79 const G4Step& stepData); 80 80 81 G4double AtRestGetPhysicalInteractionLength(81 virtual G4double AtRestGetPhysicalInteractionLength( 82 82 const G4Track& track, 83 83 G4ForceCondition* condition … … 96 96 }; 97 97 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....100 101 inline G4bool G4eplusAnnihilation::IsApplicable(const G4ParticleDefinition& p)102 {103 return (&p == G4Positron::Positron());104 }105 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....107 108 inline109 G4double G4eplusAnnihilation::AtRestGetPhysicalInteractionLength(110 const G4Track&, G4ForceCondition* condition)111 {112 *condition = NotForced;113 return 0.0;114 }115 116 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 117 99 -
trunk/source/processes/electromagnetic/standard/include/G4hIonisation.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4hIonisation.hh,v 1.4 1 2008/09/14 17:11:48vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4hIonisation.hh,v 1.42 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 93 93 virtual ~G4hIonisation(); 94 94 95 G4bool IsApplicable(const G4ParticleDefinition& p);95 virtual G4bool IsApplicable(const G4ParticleDefinition& p); 96 96 97 G4double MinPrimaryEnergy(const G4ParticleDefinition* p,98 const G4Material*, G4double cut);97 virtual G4double MinPrimaryEnergy(const G4ParticleDefinition* p, 98 const G4Material*, G4double cut); 99 99 100 100 // Print out of the class parameters … … 123 123 124 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....126 127 inline G4bool G4hIonisation::IsApplicable(const G4ParticleDefinition& p)128 {129 return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV &&130 !p.IsShortLived());131 }132 133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....134 135 inline G4double G4hIonisation::MinPrimaryEnergy(const G4ParticleDefinition*,136 const G4Material*,137 G4double cut)138 {139 G4double x = 0.5*cut/electron_mass_c2;140 G4double g = x*ratio + std::sqrt((1. + x)*(1. + x*ratio*ratio));141 return mass*(g - 1.0);142 }143 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....145 125 146 126 inline void G4hIonisation::ActivateNuclearStopping(G4bool val) -
trunk/source/processes/electromagnetic/standard/include/G4ionIonisation.hh
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4ionIonisation.hh,v 1.5 6 2008/09/14 17:11:48vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4ionIonisation.hh,v 1.57 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 85 85 virtual ~G4ionIonisation(); 86 86 87 inlineG4bool IsApplicable(const G4ParticleDefinition& p);87 virtual G4bool IsApplicable(const G4ParticleDefinition& p); 88 88 89 89 // Print out of the class parameters … … 102 102 const G4ParticleDefinition*); 103 103 104 inlineG4double MinPrimaryEnergy(const G4ParticleDefinition* p,104 virtual G4double MinPrimaryEnergy(const G4ParticleDefinition* p, 105 105 const G4Material*, G4double cut); 106 106 … … 127 127 128 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....130 131 inline G4bool G4ionIonisation::IsApplicable(const G4ParticleDefinition& p)132 {133 return (p.GetPDGCharge() != 0.0 && !p.IsShortLived() &&134 p.GetParticleType() == "nucleus");135 }136 137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....138 139 inline G4double G4ionIonisation::MinPrimaryEnergy(140 const G4ParticleDefinition* p, const G4Material*, G4double cut)141 {142 return143 p->GetPDGMass()*(std::sqrt(1. + 0.5*cut/CLHEP::electron_mass_c2) - 1.0);144 }145 146 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 147 130 -
trunk/source/processes/electromagnetic/standard/src/G4BetheBlochModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BetheBlochModel.cc,v 1. 24 2008/10/22 16:00:57vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4BetheBlochModel.cc,v 1.31 2009/04/24 14:24:00 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 65 65 #include "G4EmCorrections.hh" 66 66 #include "G4ParticleChangeForLoss.hh" 67 #include "G4NistManager.hh"68 67 69 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... … … 83 82 { 84 83 fParticleChange = 0; 85 if(p) SetParticle(p); 84 if(p) { 85 SetGenericIon(p); 86 SetParticle(p); 87 } 86 88 theElectron = G4Electron::Electron(); 87 89 corr = G4LossTableManager::Instance()->EmCorrections(); … … 108 110 const G4DataVector&) 109 111 { 110 if (!particle) SetParticle(p); 112 SetGenericIon(p); 113 SetParticle(p); 114 115 //G4cout << "G4BetheBlochModel::Initialise for " << p->GetParticleName() 116 // << " isIon= " << isIon 117 // << G4endl; 111 118 112 119 corrFactor = chargeSquare; 120 // always false before the run 121 SetDeexcitationFlag(false); 113 122 114 123 if(!isInitialised) { 115 124 isInitialised = true; 116 117 if(!fParticleChange) { 118 if (pParticleChange) { 119 fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*> 120 (pParticleChange); 121 } else { 122 fParticleChange = new G4ParticleChangeForLoss(); 123 } 124 } 125 } 126 } 127 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 129 130 void G4BetheBlochModel::SetParticle(const G4ParticleDefinition* p) 131 { 132 if(particle != p) { 133 particle = p; 134 G4String pname = particle->GetParticleName(); 135 if (particle->GetParticleType() == "nucleus" && 136 pname != "deuteron" && pname != "triton") { 137 isIon = true; 138 } 139 140 mass = particle->GetPDGMass(); 141 spin = particle->GetPDGSpin(); 142 G4double q = particle->GetPDGCharge()/eplus; 143 chargeSquare = q*q; 144 ratio = electron_mass_c2/mass; 145 G4double magmom = particle->GetPDGMagneticMoment() 146 *mass/(0.5*eplus*hbar_Planck*c_squared); 147 magMoment2 = magmom*magmom - 1.0; 148 formfact = 0.0; 149 if(particle->GetLeptonNumber() == 0) { 150 G4double x = 0.8426*GeV; 151 if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;} 152 else if(mass > GeV) { 153 x /= nist->GetZ13(mass/proton_mass_c2); 154 // tlimit = 51.2*GeV*A13[iz]*A13[iz]; 155 } 156 formfact = 2.0*electron_mass_c2/(x*x); 157 tlimit = 2.0/formfact; 158 } 125 fParticleChange = GetParticleChangeForLoss(); 159 126 } 160 127 } … … 415 382 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) / 416 383 (deltaMomentum * totMomentum); 384 /* 417 385 if(cost > 1.0) { 418 386 G4cout << "### G4BetheBlochModel WARNING: cost= " … … 429 397 cost = 1.0; 430 398 } 399 */ 431 400 G4double sint = sqrt((1.0 - cost)*(1.0 + cost)); 432 401 … … 440 409 // create G4DynamicParticle object for delta ray 441 410 G4DynamicParticle* delta = new G4DynamicParticle(theElectron, 442 411 deltaDirection,deltaKinEnergy); 443 412 444 413 vdp->push_back(delta); … … 453 422 } 454 423 455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 425 426 G4double G4BetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd, 427 G4double kinEnergy) 428 { 429 // here particle type is checked for any method 430 SetParticle(pd); 431 G4double tau = kinEnergy/mass; 432 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) / 433 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio); 434 return std::min(tmax,tlimit); 435 } 436 437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -
trunk/source/processes/electromagnetic/standard/src/G4BetheHeitlerModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BetheHeitlerModel.cc,v 1.1 2 2008/10/15 15:54:57vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4BetheHeitlerModel.cc,v 1.13 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 94 94 const G4DataVector&) 95 95 { 96 if(!fParticleChange) { 97 if(pParticleChange) { 98 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange); 99 } else { 100 fParticleChange = new G4ParticleChangeForGamma(); 101 } 102 } 96 if(!fParticleChange) fParticleChange = GetParticleChangeForGamma(); 103 97 104 98 if(theCrossSectionTable) { -
trunk/source/processes/electromagnetic/standard/src/G4BohrFluctuations.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BohrFluctuations.cc,v 1. 6 2007/09/27 14:02:41vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4BohrFluctuations.cc,v 1.7 2009/02/19 19:17:50 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 137 137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 138 138 139 G4double G4BohrFluctuations::Dispersion(const G4Material* material, 140 const G4DynamicParticle* dp, 141 G4double& tmax, 142 G4double& length) 143 { 144 if(!particle) InitialiseMe(dp->GetDefinition()); 139 145 146 G4double electronDensity = material->GetElectronDensity(); 147 kineticEnergy = dp->GetKineticEnergy(); 148 G4double etot = kineticEnergy + particleMass; 149 beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot); 150 G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length 151 * electronDensity * chargeSquare; 152 153 return siga; 154 } 155 156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 157 158 -
trunk/source/processes/electromagnetic/standard/src/G4BraggIonModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BraggIonModel.cc,v 1.2 2 2008/10/22 16:00:57vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4BraggIonModel.cc,v 1.24 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 114 114 corrFactor = chargeSquare; 115 115 116 // always false before the run 117 SetDeexcitationFlag(false); 118 116 119 if(!isInitialised) { 117 120 isInitialised = true; … … 123 126 corr = G4LossTableManager::Instance()->EmCorrections(); 124 127 125 if(!fParticleChange) { 126 if(pParticleChange) { 127 fParticleChange = 128 reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange); 129 } else { 130 fParticleChange = new G4ParticleChangeForLoss(); 131 } 132 } 128 fParticleChange = GetParticleChangeForLoss(); 133 129 } 134 130 } … … 352 348 fParticleChange->SetProposedKineticEnergy(kineticEnergy); 353 349 fParticleChange->SetProposedMomentumDirection(finalP); 350 } 351 352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 353 354 G4double G4BraggIonModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd, 355 G4double kinEnergy) 356 { 357 if(pd != particle) SetParticle(pd); 358 G4double tau = kinEnergy/mass; 359 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) / 360 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio); 361 return tmax; 354 362 } 355 363 -
trunk/source/processes/electromagnetic/standard/src/G4BraggModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BraggModel.cc,v 1.2 0 2008/10/22 16:01:46vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4BraggModel.cc,v 1.22 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 112 112 if(p != particle) SetParticle(p); 113 113 114 // always false before the run 115 SetDeexcitationFlag(false); 116 114 117 if(!isInitialised) { 115 118 isInitialised = true; … … 121 124 corr = G4LossTableManager::Instance()->EmCorrections(); 122 125 123 if(pParticleChange) { 124 fParticleChange = 125 reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange); 126 } else { 127 fParticleChange = new G4ParticleChangeForLoss(); 128 } 126 fParticleChange = GetParticleChangeForLoss(); 129 127 } 130 128 } … … 337 335 338 336 vdp->push_back(delta); 337 } 338 339 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 340 341 G4double G4BraggModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd, 342 G4double kinEnergy) 343 { 344 if(pd != particle) SetParticle(pd); 345 G4double tau = kinEnergy/mass; 346 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) / 347 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio); 348 return tmax; 339 349 } 340 350 -
trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4ComptonScattering.cc,v 1.3 0 2008/10/15 17:53:44vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4ComptonScattering.cc,v 1.31 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // … … 84 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 85 86 G4bool G4ComptonScattering::IsApplicable(const G4ParticleDefinition& p) 87 { 88 return (&p == G4Gamma::Gamma()); 89 } 90 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 92 86 93 void G4ComptonScattering::InitialiseProcess(const G4ParticleDefinition*) 87 94 { -
trunk/source/processes/electromagnetic/standard/src/G4CoulombScattering.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4CoulombScattering.cc,v 1. 19 2008/10/15 17:53:44vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4CoulombScattering.cc,v 1.20 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 79 79 G4CoulombScattering::~G4CoulombScattering() 80 80 {} 81 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 83 84 G4bool G4CoulombScattering::IsApplicable(const G4ParticleDefinition& p) 85 { 86 return (p.GetPDGCharge() != 0.0 && !p.IsShortLived()); 87 } 81 88 82 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... -
trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4CoulombScatteringModel.cc,v 1.3 7 2008/07/31 13:11:34vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4CoulombScatteringModel.cc,v 1.39 2009/05/10 16:09:29 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 93 93 94 94 // CM system 95 G4int iz= G4int(Z);95 iz = G4int(Z); 96 96 G4double m2 = fNistManager->GetAtomicMassAmu(iz)*amu_c2; 97 97 G4double etot = tkin + mass; … … 137 137 138 138 G4double Z = currentElement->GetZ(); 139 G4int iz= G4int(Z);139 iz = G4int(Z); 140 140 G4int ia = SelectIsotopeNumber(currentElement); 141 141 G4double m2 = theParticleTable->GetIonTable()->GetNucleusMass(iz, ia); … … 182 182 // recoil 183 183 G4double erec = kinEnergy - ekin; 184 G4double th =185 std::min(recoilThreshold,186 Z*currentElement->GetIonisation()->GetMeanExcitationEnergy());187 184 188 if(erec > th) {185 if(erec > recoilThreshold) { 189 186 G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz); 190 187 G4double plab = sqrt(ekin*(ekin + 2.0*mass)); -
trunk/source/processes/electromagnetic/standard/src/G4GammaConversion.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4GammaConversion.cc,v 1.3 0 2008/10/15 17:53:44vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4GammaConversion.cc,v 1.31 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // … … 90 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 91 91 92 G4bool G4GammaConversion::IsApplicable(const G4ParticleDefinition& p) 93 { 94 return (&p == G4Gamma::Gamma()); 95 } 96 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 98 92 99 void G4GammaConversion::InitialiseProcess(const G4ParticleDefinition*) 93 100 { -
trunk/source/processes/electromagnetic/standard/src/G4IonFluctuations.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4IonFluctuations.cc,v 1.2 4 2008/10/22 16:25:21 vnivanchExp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4IonFluctuations.cc,v 1.26 2009/03/31 13:24:40 toshito Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 190 190 191 191 // heavy ion correction 192 G4double f1 = 1.065e-4*chargeSquare;193 if(beta2 > theBohrBeta2) f1/= beta2;194 else f1/= theBohrBeta2;195 if(f1 > 2.5) f1 = 2.5;196 fac *= (1.0 + f1);192 // G4double f1 = 1.065e-4*chargeSquare; 193 // if(beta2 > theBohrBeta2) f1/= beta2; 194 // else f1/= theBohrBeta2; 195 // if(f1 > 2.5) f1 = 2.5; 196 // fac *= (1.0 + f1); 197 197 198 198 // taking into account the cut 199 if(fac > 1.0) { 200 siga *= (1.0 + (fac - 1.0)*2.0*electron_mass_c2*beta2/(tmax*(1.0 - beta2))); 201 } 199 G4double fac_cut = 1.0 + (fac - 1.0)*2.0*electron_mass_c2*beta2/(tmax*(1.0 - beta2)); 200 if(fac_cut > 0.01 && fac > 0.01) { 201 siga *= fac_cut; 202 } 203 202 204 //G4cout << "siga(keV)= " << sqrt(siga)/keV << " fac= " << fac 203 205 // << " f1= " << f1 << G4endl; … … 420 422 421 423 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 424 425 void G4IonFluctuations::SetParticleAndCharge(const G4ParticleDefinition* part, 426 G4double q2) 427 { 428 if(part != particle) { 429 particle = part; 430 particleMass = part->GetPDGMass(); 431 charge = part->GetPDGCharge()/eplus; 432 chargeSquare = charge*charge; 433 } 434 effChargeSquare = q2; 435 uniFluct.SetParticleAndCharge(part, q2); 436 } 437 438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... -
trunk/source/processes/electromagnetic/standard/src/G4KleinNishinaCompton.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4KleinNishinaCompton.cc,v 1. 9 2007/05/22 17:34:36vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4KleinNishinaCompton.cc,v 1.10 2009/05/15 17:12:33 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 79 79 const G4DataVector&) 80 80 { 81 if(pParticleChange) 82 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange); 83 else 84 fParticleChange = new G4ParticleChangeForGamma(); 81 fParticleChange = GetParticleChangeForGamma(); 85 82 } 86 83 -
trunk/source/processes/electromagnetic/standard/src/G4MollerBhabhaModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4MollerBhabhaModel.cc,v 1.3 0 2007/05/22 17:34:36 vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4MollerBhabhaModel.cc,v 1.34 2009/04/24 17:15:46 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 74 74 const G4String& nam) 75 75 : G4VEmModel(nam), 76 particle(0), 77 isElectron(true), 78 twoln10(2.0*log(10.0)), 79 lowLimit(0.2*keV) 76 particle(0), 77 isElectron(true), 78 twoln10(2.0*log(10.0)), 79 lowLimit(0.2*keV), 80 isInitialised(false) 80 81 { 81 82 theElectron = G4Electron::Electron(); … … 87 88 G4MollerBhabhaModel::~G4MollerBhabhaModel() 88 89 {} 89 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......91 92 void G4MollerBhabhaModel::SetParticle(const G4ParticleDefinition* p)93 {94 particle = p;95 if(p != theElectron) isElectron = false;96 }97 90 98 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... … … 108 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 109 102 103 G4double G4MollerBhabhaModel::MaxSecondaryEnergy(const G4ParticleDefinition*, 104 G4double kinEnergy) 105 { 106 G4double tmax = kinEnergy; 107 if(isElectron) tmax *= 0.5; 108 return tmax; 109 } 110 111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 112 110 113 void G4MollerBhabhaModel::Initialise(const G4ParticleDefinition* p, 111 114 const G4DataVector&) 112 115 { 113 116 if(!particle) SetParticle(p); 114 if(pParticleChange) 115 fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*> 116 (pParticleChange); 117 else 118 fParticleChange = new G4ParticleChangeForLoss(); 117 SetDeexcitationFlag(false); 118 119 if(isInitialised) return; 120 121 isInitialised = true; 122 fParticleChange = GetParticleChangeForLoss(); 119 123 } 120 124 … … 289 293 G4double maxEnergy) 290 294 { 291 G4double tmax = std::min(maxEnergy, MaxSecondaryKinEnergy(dp)); 295 G4double kineticEnergy = dp->GetKineticEnergy(); 296 G4double tmax = kineticEnergy; 297 if(isElectron) tmax *= 0.5; 298 if(maxEnergy < tmax) tmax = maxEnergy; 292 299 if(tmin >= tmax) return; 293 300 294 G4double kineticEnergy = dp->GetKineticEnergy();295 301 G4double energy = kineticEnergy + electron_mass_c2; 296 302 G4double totalMomentum = sqrt(kineticEnergy*(energy + electron_mass_c2)); … … 339 345 G4double b3 = b4 + y122; 340 346 341 y = xmax*xmax; 342 grej = -xmin*b1; 343 grej += y*b2; 344 grej -= xmin*xmin*xmin*b3; 345 grej += y*y*b4; 346 grej *= beta2; 347 grej += 1.0; 347 y = xmax*xmax; 348 grej = 1.0 + (y*y*b4 - xmin*xmin*xmin*b3 + y*b2 - xmin*b1)*beta2; 348 349 do { 349 q = G4UniformRand(); 350 x = xmin*xmax/(xmin*(1.0 - q) + xmax*q); 351 z = -x*b1; 352 y = x*x; 353 z += y*b2; 354 y *= x; 355 z -= y*b3; 356 y *= x; 357 z += y*b4; 358 z *= beta2; 359 z += 1.0; 350 q = G4UniformRand(); 351 x = xmin*xmax/(xmin*(1.0 - q) + xmax*q); 352 y = x*x; 353 z = 1.0 + (y*y*b4 - x*y*b3 + y*b2 - x*b1)*beta2; 360 354 /* 361 355 if(z > grej) { … … 396 390 // create G4DynamicParticle object for delta ray 397 391 G4DynamicParticle* delta = new G4DynamicParticle(theElectron, 398 deltaDirection,deltaKinEnergy);392 deltaDirection,deltaKinEnergy); 399 393 vdp->push_back(delta); 400 394 } -
trunk/source/processes/electromagnetic/standard/src/G4MscModel71.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4MscModel71.cc,v 1. 6 2008/03/13 17:20:07vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4MscModel71.cc,v 1.7 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 140 140 141 141 if(pParticleChange) 142 fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange);142 fParticleChange = static_cast<G4ParticleChangeForMSC*>(pParticleChange); 143 143 else 144 144 fParticleChange = new G4ParticleChangeForMSC(); -
trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PAIModel.cc,v 1.47 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 // 29 // ------------------------------------------------------------------- 30 // 31 // GEANT4 Class 26 32 // File name: G4PAIModel.cc 27 33 // … … 153 159 fTotBin); 154 160 155 if(pParticleChange) 156 fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange); 157 else 158 fParticleChange = new G4ParticleChangeForLoss(); 161 fParticleChange = GetParticleChangeForLoss(); 159 162 160 163 // Prepare initialization … … 177 180 fCutCouple = theCoupleTable->GetMaterialCutsCouple( fMaterial, 178 181 curReg->GetProductionCuts() ); 182 //G4cout << "Reg <" <<curReg->GetName() << "> mat <" 183 // << fMaterial->GetName() << "> fCouple= " 184 // << fCutCouple<<G4endl; 179 185 if( fCutCouple ) { 180 186 fMaterialCutsCoupleVector.push_back(fCutCouple); … … 197 203 } 198 204 } 205 206 ////////////////////////////////////////////////////////////////// 207 208 void G4PAIModel::InitialiseMe(const G4ParticleDefinition*) 209 {} 199 210 200 211 ////////////////////////////////////////////////////////////////// … … 393 404 { 394 405 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 395 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 406 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 407 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ; 396 408 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 397 409 } … … 435 447 { 436 448 // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; 437 if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; 449 // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; 450 if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ; 438 451 else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 439 452 } … … 444 457 ////////////////////////////////////////////////////////////////////////////// 445 458 446 G4double G4PAIModel::ComputeDEDX (const G4MaterialCutsCouple* matCC,447 448 449 459 G4double G4PAIModel::ComputeDEDXPerVolume(const G4Material*, 460 const G4ParticleDefinition* p, 461 G4double kineticEnergy, 462 G4double cutEnergy) 450 463 { 451 464 G4int iTkin,iPlace; 452 465 size_t jMat; 466 467 //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy); 468 G4double cut = cutEnergy; 469 453 470 G4double massRatio = fMass/p->GetPDGMass(); 454 471 G4double scaledTkin = kineticEnergy*massRatio; 455 472 G4double charge = p->GetPDGCharge(); 456 G4double charge2 = charge*charge, dEdx; 473 G4double charge2 = charge*charge; 474 const G4MaterialCutsCouple* matCC = CurrentCouple(); 457 475 458 476 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) … … 470 488 iPlace = iTkin - 1; 471 489 if(iPlace < 0) iPlace = 0; 472 dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cutEnergy) ) ; 473 490 G4double dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) ); 474 491 if( dEdx < 0.) dEdx = 0.; 475 492 return dEdx; … … 478 495 ///////////////////////////////////////////////////////////////////////// 479 496 480 G4double G4PAIModel::CrossSection ( const G4MaterialCutsCouple* matCC,481 482 483 484 497 G4double G4PAIModel::CrossSectionPerVolume( const G4Material*, 498 const G4ParticleDefinition* p, 499 G4double kineticEnergy, 500 G4double cutEnergy, 501 G4double maxEnergy ) 485 502 { 486 503 G4int iTkin,iPlace; 487 504 size_t jMat; 488 G4double tmax = min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy); 505 G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy); 506 if(tmax <= cutEnergy) return 0.0; 489 507 G4double massRatio = fMass/p->GetPDGMass(); 490 508 G4double scaledTkin = kineticEnergy*massRatio; 491 509 G4double charge = p->GetPDGCharge(); 492 510 G4double charge2 = charge*charge, cross, cross1, cross2; 511 const G4MaterialCutsCouple* matCC = CurrentCouple(); 493 512 494 513 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) … … 935 954 } 936 955 956 ///////////////////////////////////////////////////////////////////// 957 958 G4double G4PAIModel::MaxSecondaryEnergy( const G4ParticleDefinition* p, 959 G4double kinEnergy) 960 { 961 G4double tmax = kinEnergy; 962 if(p == fElectron) tmax *= 0.5; 963 else if(p != fPositron) { 964 G4double mass = p->GetPDGMass(); 965 G4double ratio= electron_mass_c2/mass; 966 G4double gamma= kinEnergy/mass + 1.0; 967 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) / 968 (1. + 2.0*gamma*ratio + ratio*ratio); 969 } 970 return tmax; 971 } 972 973 /////////////////////////////////////////////////////////////// 974 975 void G4PAIModel::DefineForRegion(const G4Region* r) 976 { 977 fPAIRegionVector.push_back(r); 978 } 937 979 938 980 // -
trunk/source/processes/electromagnetic/standard/src/G4PAIPhotonModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PAIPhotonModel.cc,v 1.22 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 // 29 // ------------------------------------------------------------------- 30 // 31 // GEANT4 Class 26 32 // File name: G4PAIPhotonModel.cc 27 33 // … … 160 166 if(!fParticle) SetParticle(p); 161 167 162 if(pParticleChange) 163 fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange); 164 else 165 fParticleChange = new G4ParticleChangeForLoss(); 168 fParticleChange = GetParticleChangeForLoss(); 166 169 167 170 const G4ProductionCutsTable* theCoupleTable = … … 214 217 } 215 218 } 219 220 ////////////////////////////////////////////////////////////////// 221 222 void G4PAIPhotonModel::InitialiseMe(const G4ParticleDefinition*) 223 {} 216 224 217 225 ////////////////////////////////////////////////////////////////// … … 487 495 { 488 496 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 489 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 497 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 498 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ; 490 499 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 491 500 } … … 530 539 { 531 540 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 532 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 541 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 542 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ; 533 543 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 534 544 } … … 574 584 { 575 585 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 576 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 586 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 587 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ; 577 588 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 578 589 } … … 617 628 { 618 629 // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; 619 if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; 630 // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; 631 if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ; 620 632 else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 621 633 } … … 626 638 ////////////////////////////////////////////////////////////////////////////// 627 639 628 G4double G4PAIPhotonModel::ComputeDEDX (const G4MaterialCutsCouple* matCC,629 630 631 640 G4double G4PAIPhotonModel::ComputeDEDXPerVolume(const G4Material*, 641 const G4ParticleDefinition* p, 642 G4double kineticEnergy, 643 G4double cutEnergy) 632 644 { 633 645 G4int iTkin,iPlace; 634 646 size_t jMat; 647 648 //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy); 649 G4double cut = cutEnergy; 650 635 651 G4double particleMass = p->GetPDGMass(); 636 652 G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass; … … 638 654 G4double charge2 = charge*charge; 639 655 G4double dEdx = 0.; 656 const G4MaterialCutsCouple* matCC = CurrentCouple(); 640 657 641 658 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) … … 653 670 iPlace = iTkin - 1; 654 671 if(iPlace < 0) iPlace = 0; 655 dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut Energy) ) ;672 dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) ) ; 656 673 657 674 if( dEdx < 0.) dEdx = 0.; … … 661 678 ///////////////////////////////////////////////////////////////////////// 662 679 663 G4double G4PAIPhotonModel::CrossSection ( const G4MaterialCutsCouple* matCC,664 665 666 667 680 G4double G4PAIPhotonModel::CrossSectionPerVolume( const G4Material*, 681 const G4ParticleDefinition* p, 682 G4double kineticEnergy, 683 G4double cutEnergy, 684 G4double maxEnergy ) 668 685 { 669 686 G4int iTkin,iPlace; 670 687 size_t jMat, jMatCC; 671 G4double tmax = min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy); 688 G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy); 689 if(cutEnergy >= tmax) return 0.0; 672 690 G4double particleMass = p->GetPDGMass(); 673 691 G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass; … … 675 693 G4double charge2 = charge*charge, cross, cross1, cross2; 676 694 G4double photon1, photon2, plasmon1, plasmon2; 695 696 const G4MaterialCutsCouple* matCC = CurrentCouple(); 677 697 678 698 const G4ProductionCutsTable* theCoupleTable= … … 1225 1245 } 1226 1246 1247 ///////////////////////////////////////////////////////////////////// 1248 1249 G4double G4PAIPhotonModel::MaxSecondaryEnergy( const G4ParticleDefinition* p, 1250 G4double kinEnergy) 1251 { 1252 G4double tmax = kinEnergy; 1253 if(p == fElectron) tmax *= 0.5; 1254 else if(p != fPositron) { 1255 G4double mass = p->GetPDGMass(); 1256 G4double ratio= electron_mass_c2/mass; 1257 G4double gamma= kinEnergy/mass + 1.0; 1258 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) / 1259 (1. + 2.0*gamma*ratio + ratio*ratio); 1260 } 1261 return tmax; 1262 } 1263 1264 /////////////////////////////////////////////////////////////// 1265 1266 void G4PAIPhotonModel::DefineForRegion(const G4Region* r) 1267 { 1268 fPAIRegionVector.push_back(r); 1269 } 1270 1227 1271 1228 1272 // -
trunk/source/processes/electromagnetic/standard/src/G4PEEffectModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PEEffectModel.cc,v 1. 6 2007/05/22 17:34:36vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4PEEffectModel.cc,v 1.8 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 41 41 // 42 42 // 04.12.05 : SetProposedKineticEnergy(0.) for the killed photon (mma) 43 // 20.02.09 : Added initialisation of deexcitation flag and method 44 // CrossSectionPerVolume instead of mfp (V.Ivanchenko) 43 45 // 44 46 // Class Description: … … 66 68 theGamma = G4Gamma::Gamma(); 67 69 theElectron = G4Electron::Electron(); 70 fminimalEnergy = 1.0*eV; 68 71 } 69 72 … … 71 74 72 75 G4PEEffectModel::~G4PEEffectModel() 73 { 74 } 76 {} 75 77 76 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... … … 79 81 const G4DataVector&) 80 82 { 81 if (isInitialized) return; 82 if (pParticleChange) 83 fParticleChange = 84 reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange); 85 else 86 fParticleChange = new G4ParticleChangeForGamma(); 87 88 fminimalEnergy = 1.0*eV; 83 // always false before the run 84 SetDeexcitationFlag(false); 85 86 if (isInitialized) return; 87 fParticleChange = GetParticleChangeForGamma(); 88 isInitialized = true; 89 } 90 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 92 93 G4double G4PEEffectModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 94 G4double energy, 95 G4double Z, G4double, 96 G4double, G4double) 97 { 98 G4double* SandiaCof = G4SandiaTable::GetSandiaCofPerAtom((G4int)Z, energy); 99 100 G4double energy2 = energy*energy; 101 G4double energy3 = energy*energy2; 102 G4double energy4 = energy2*energy2; 103 104 return SandiaCof[0]/energy + SandiaCof[1]/energy2 + 105 SandiaCof[2]/energy3 + SandiaCof[3]/energy4; 106 } 107 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 109 110 G4double G4PEEffectModel::CrossSectionPerVolume(const G4Material* material, 111 const G4ParticleDefinition*, 112 G4double energy, 113 G4double, G4double) 114 { 115 G4double* SandiaCof = 116 material->GetSandiaTable()->GetSandiaCofForMaterial(energy); 117 118 G4double energy2 = energy*energy; 119 G4double energy3 = energy*energy2; 120 G4double energy4 = energy2*energy2; 121 122 return SandiaCof[0]/energy + SandiaCof[1]/energy2 + 123 SandiaCof[2]/energy3 + SandiaCof[3]/energy4; 89 124 } 90 125 -
trunk/source/processes/electromagnetic/standard/src/G4PhotoElectricEffect.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PhotoElectricEffect.cc,v 1.4 1 2008/10/16 14:12:32vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4PhotoElectricEffect.cc,v 1.42 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // … … 89 89 {} 90 90 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 92 93 G4bool G4PhotoElectricEffect::IsApplicable(const G4ParticleDefinition& p) 94 { 95 return (&p == G4Gamma::Gamma()); 96 } 97 91 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 92 99 -
trunk/source/processes/electromagnetic/standard/src/G4UniversalFluctuation.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4UniversalFluctuation.cc,v 1. 16 2008/10/22 16:04:33 vnivanchExp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4UniversalFluctuation.cc,v 1.22 2009/03/20 18:11:23 urban Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 57 57 // 03-04-07 correction to get better width of eloss distr.(L.Urban) 58 58 // 13-07-07 add protection for very small step or low-density material (VI) 59 // 59 // 19-03-09 new width correction (does not depend on previous steps) (L.Urban) 60 // 20-03-09 modification in the width correction (L.Urban) 61 // 60 62 61 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... … … 84 86 { 85 87 lastMaterial = 0; 86 facwidth = 1.000/keV;87 oldloss = 0.;88 samestep = 0.;89 88 } 90 89 … … 120 119 // 121 120 if (meanLoss < minLoss) 122 {123 oldloss = meanLoss;124 121 return meanLoss; 125 }126 122 127 123 if(!particle) InitialiseMe(dp->GetDefinition()); … … 179 175 e0 = material->GetIonisation()->GetEnergy0fluct(); 180 176 lastMaterial = material; 177 178 // modification of some model parameters 179 // (this part should go to materials later) 180 G4double p = 1.40; 181 f2Fluct *= p; 182 f1Fluct = 1.-f2Fluct; 183 G4double q = 1.00; 184 e2Fluct *= q; 185 e2LogFluct = log(e2Fluct); 186 e1LogFluct = (ipotLogFluct-f2Fluct*e2LogFluct)/f1Fluct; 187 e1Fluct = exp(e1LogFluct); 181 188 } 182 189 … … 185 192 186 193 G4double a1 = 0. , a2 = 0., a3 = 0. ; 187 188 // correction to get better width even using stepmax189 if(abs(meanLoss- oldloss) < 1.*eV)190 samestep += 1;191 else192 samestep = 1.;193 oldloss = meanLoss;194 G4double width = 1.+samestep*facwidth*meanLoss;195 if(width > 4.50) width = 4.50;196 e1 = width*e1Fluct;197 e2 = width*e2Fluct;198 194 199 195 // cut and material dependent rate … … 206 202 rate = 0.03+0.23*log(log(tmax/ipotFluct)); 207 203 G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct); 208 a1 = C*f1Fluct*(w2-e1LogFluct)/e1; 209 a2 = C*f2Fluct*(w2-e2LogFluct)/e2; 204 a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct; 205 a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct; 206 // correction in order to get better FWHM values 207 // ( scale parameters a1 and e1) 208 G4double width = 1.; 209 if(meanLoss > 10.*e1Fluct) 210 { 211 width = 3.1623/sqrt(meanLoss/e1Fluct); 212 if(width < a2/a1) 213 width = a2/a1; 214 } 215 a1 *= width; 216 e1 = e1Fluct/width; 217 e2 = e2Fluct; 210 218 } 211 219 } … … 305 313 } 306 314 307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 316 317 void 318 G4UniversalFluctuation::SetParticleAndCharge(const G4ParticleDefinition* part, 319 G4double q2) 320 { 321 if(part != particle) { 322 particle = part; 323 particleMass = part->GetPDGMass(); 324 } 325 chargeSquare = q2; 326 } 327 328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -
trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4UrbanMscModel.cc,v 1.86 2008/10/29 14:15:30 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02 $ 26 // 27 // $Id: G4UrbanMscModel.cc,v 1.90 2009/04/29 13:30:22 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 29 // 29 30 // ------------------------------------------------------------------- … … 36 37 // Author: Laszlo Urban 37 38 // 38 // Creation date: 0 3.03.200139 // Creation date: 06.03.2008 39 40 // 40 41 // Modifications: 41 42 // 42 // 27-03-03 Move model part from G4MultipleScattering80 (V.Ivanchenko) 43 // 23-05-03 important change in angle distribution for muons/hadrons 44 // the central part now is similar to the Highland parametrization + 45 // minor correction in angle sampling algorithm (for all particles) 46 // (L.Urban) 47 // 30-05-03 misprint in SampleCosineTheta corrected(L.Urban) 48 // 27-03-03 Rename (V.Ivanchenko) 49 // 05-08-03 angle distribution has been modified (L.Urban) 50 // 06-11-03 precision problems solved for high energy (PeV) particles 51 // change in the tail of the angular distribution 52 // highKinEnergy is set to 100 PeV (L.Urban) 53 // 54 // 10-11-03 highKinEnergy is set back to 100 TeV, some tail tuning + 55 // cleaning (L.Urban) 56 // 26-11-03 correction in TrueStepLength : 57 // trueLength <= currentRange (L.Urban) 58 // 01-03-04 signature changed in SampleCosineTheta, 59 // energy dependence calculations has been simplified, 60 // 11-03-04 corrections in GeomPathLength,TrueStepLength, 61 // SampleCosineTheta 62 // 23-04-04 true -> geom and geom -> true transformation has been 63 // rewritten, changes in the angular distribution (L.Urban) 64 // 19-07-04 correction in SampleCosineTheta in order to avoid 65 // num. precision problems at high energy/small step(L.Urban) 66 // 17-08-04 changes in the angle distribution (slightly modified 67 // Highland formula for the width of the central part, 68 // changes in the numerical values of some other parameters) 69 // ---> approximately step independent distribution (L.Urban) 70 // 21-09-04 change in the tail of the angular distribution (L.Urban) 71 // 72 // 03-11-04 precision problem for very high energy ions and small stepsize 73 // solved in SampleCosineTheta (L.Urban). 74 // 15-04-05 optimize internal interface 75 // add SampleSecondaries method (V.Ivanchenko) 76 // 11-08-05 computation of lateral correlation added (L.Urban) 77 // 02-10-05 nuclear size correction computation removed, the correction 78 // included in the (theoretical) tabulated values (L.Urban) 79 // 17-01-06 computation of tail changed in SampleCosineTheta (l.Urban) 80 // 16-02-06 code cleaning + revised 'z' sampling (L.Urban) 81 // 17-02-06 Save table of transport cross sections not mfp (V.Ivanchenko) 82 // 07-03-06 Create G4UrbanMscModel and move there step limit 83 // calculation (V.Ivanchenko) 84 // 23-03-06 Bugfix in SampleCosineTheta method (L.Urban) 85 // 10-05-06 SetMscStepLimitation at initialisation (V.Ivantchenko) 86 // 11-05-06 name of data member safety changed to presafety, some new data 87 // members added (frscaling1,frscaling2,tlimitminfix,nstepmax) 88 // changes in ComputeTruePathLengthLimit,SampleCosineTheta (L.Urban) 89 // 17-05-06 parameters of theta0 in SampleCosineTheta changed 90 // c_highland 13.6*MeV ---> 13.26*MeV, 91 // corr_highland 0.555 ---> 0.54, 92 // value of data member geommin changed (5 nm -> 1 nm) (L.Urban) 93 // 13-10-06 data member factail removed, data member tkinlimit changed 94 // to lambdalimit, 95 // new data members tgeom,tnow,skin,skindepth,Zeff,geomlimit 96 // G4double GeomLimit(const G4Track& track) changed to 97 // void GeomLimit(const G4Track& track) 98 // - important changes in ComputeTruePathLengthLimit: 99 // possibility to have very small step(s) with single scattering 100 // before boundary crossing (with skin > 0) 101 // - changes in SampleCosineTheta : 102 // single scattering if step <= stepmin, parameter theta0 103 // slightly modified, tail modified (L.Urban) 104 // 20-10-06 parameter theta0 now computed in the (public) 105 // function ComputeTheta0, 106 // single scattering modified allowing not small 107 // angles as well (L.Urban) 108 // 23-10-06 correction in SampleSecondaries, now safety update 109 // computed in a simpler/faster way (L.Urban) 110 // 06-11-06 corrections in ComputeTruePathLengthLimit, results are 111 // more stable in calorimeters (L.Urban) 112 // 07-11-06 fix in GeomPathLength and SampleCosineTheta (L.Urban) 113 // 15-11-06 bugfix in SampleCosineTheta (L.Urban) 114 // 20-11-06 bugfix in single scattering part of SampleCosineTheta, 115 // single scattering just before boundary crossing now (L.Urban) 116 // 04-12-06 fix in ComputeTruePathLengthLimit (L.Urban) 117 // 17-01-07 remove LocatePoint from GeomLimit method (V.Ivanchenko) 118 // 19-01-07 fix of true < geom problem (L.Urban) 119 // 25-01-07 add protections from NaN vaues and for zero geometry step (VI) 120 // 31-01-07 correction in SampleCosineTheta: screening parameter 121 // corrected in single/plural scattering + 122 // code cleaning (L.Urban) 123 // 01-02-07 restore logic inside ComputeTrueStepLength (V.Ivanchenko) 124 // 06-02-07 Move SetMscStepLimitation method into the source, add there 125 // reinitialisation of some private members, add protection inside 126 // SampleDisplacement(VI) 127 // 07-02-07 fix single scattering for heavy particles, now skin=1 can be used 128 // for heavy particles as well (L.Urban) 129 // 08-02-07 randomization of tlimit removed (L.Urban) 130 // 11-02-07 modified stepping algorithm for skin=0 131 // 15-02-07 new data member: smallstep, small steps with single scattering 132 // before + after boundary for skin > 1 133 // 23-02-07 use tPathLength inside ComputeStep instead of geombig 134 // 24-02-07 step reduction before boundary for 'small' geomlimit only 135 // 03-03-07 single scattering around boundaries only (L.Urban) 136 // 07-03-07 bugfix in ComputeTruePathLengthLimit (for skin > 0.) (L.Urban) 137 // 10-04-07 optimize logic of ComputeTruePathLengthLimit, remove 138 // unused members, use unique G4SafetyHelper (V.Ivanchenko) 139 // 01-05-07 optimization for skin > 0 (L.Urban) 140 // 05-07-07 modified model functions in SampleCosineTheta (L.Urban) 141 // 06-07-07 theta0 is not the same for e-/e+ as for heavy particles (L.Urban) 142 // 02-08-07 compare safety not with 0. but with tlimitminfix (V.Ivanchenko) 143 // 09-08-07 tail of angular distribution has been modified (L.Urban) 144 // 22-10-07 - corr. in ComputeGeomPathLength in order to get better low 145 // energy behaviour for heavy particles, 146 // - theta0 is slightly modified, 147 // - some old inconsistency/bug is cured in SampleCosineTheta, 148 // now 0 <= prob <= 1 in any case (L.Urban) 149 // 26-10-07 different correction parameters for e/mu/hadrons in ComputeTheta0 150 // 30-11-07 fix in ComputeTheta0 (L.Urban) 151 // 43 // 06-03-2008 starting point : G4UrbanMscModel2 = G4UrbanMscModel 9.1 ref 02 44 // 45 // 13-03-08 Bug in SampleScattering (which caused lateral asymmetry) fixed 46 // (L.Urban) 47 // 48 // 14-03-08 Simplification of step limitation in ComputeTruePathLengthLimit, 49 // + tlimitmin is the same for UseDistancetoBoundary and 50 // UseSafety (L.Urban) 51 // 52 // 16-03-08 Reorganization of SampleCosineTheta + new method SimpleScattering 53 // SimpleScattering is used if the relative energy loss is too big 54 // or theta0 is too big (see data members rellossmax, theta0max) 55 // (L.Urban) 56 // 57 // 17-03-08 tuning of the correction factor in ComputeTheta0 (L.Urban) 58 // 59 // 19-03-08 exponent c of the 'tail' model function is not equal to 2 any more, 60 // value of c has been extracted from some e- scattering data (L.Urban) 61 // 62 // 24-03-08 Step limitation in ComputeTruePathLengthLimit has been 63 // simplified further + some data members have been removed (L.Urban) 64 // 65 // 24-07-08 central part of scattering angle (theta0) has been tuned 66 // tail of the scattering angle distribution has been tuned 67 // using some e- and proton scattering data 68 // 69 // 05-08-08 bugfix in ComputeTruePathLengthLimit (L.Urban) 70 // 71 // 09-10-08 theta0 and tail have been retuned using some e-,mu,proton 72 // scattering data (L.Urban) 73 // + single scattering without path length correction for 74 // small steps (t < tlimitmin, for UseDistanceToBoundary only) 75 // 76 // 15-10-08 Moliere-Bethe screening in the single scattering part(L.Urban) 77 // 78 // 17-10-08 stepping similar to that in model (9.1) for UseSafety case 79 // for e+/e- in order to speed up the code for calorimeters 80 // 81 // 23-10-08 bugfix in the screeningparameter of the single scattering part, 82 // some technical change in order to speed up the code (UpdateCache) 83 // 84 // 27-10-08 bugfix in ComputeTruePathLengthLimit (affects UseDistanceToBoundary 85 // stepping type only) (L.Urban) 86 // 87 // 28-04-09 move G4UrbanMscModel2 from the g49.2 to G4UrbanMscModel. 88 // now it is frozen (V.Ivanchenk0) 152 89 153 90 // Class Description: … … 157 94 158 95 // ------------------------------------------------------------------- 96 // In its present form the model can be used for simulation 97 // of the e-/e+, muon and charged hadron multiple scattering 159 98 // 160 99 … … 166 105 #include "Randomize.hh" 167 106 #include "G4Electron.hh" 168 169 107 #include "G4LossTableManager.hh" 170 108 #include "G4ParticleChangeForMSC.hh" 171 #include "G4TransportationManager.hh"172 #include "G4SafetyHelper.hh"173 109 174 110 #include "G4Poisson.hh" 111 #include "globals.hh" 175 112 176 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... … … 182 119 isInitialized(false) 183 120 { 184 masslimite = 0.6*MeV; 185 masslimitmu = 110.*MeV; 186 121 masslimite = 0.6*MeV; 122 lambdalimit = 1.*mm; 123 fr = 0.02; 124 facsafety = 0.3; 187 125 taubig = 8.0; 188 126 tausmall = 1.e-16; … … 193 131 smallstep = 1.e10; 194 132 currentRange = 0. ; 195 frscaling2 = 0.25; 196 frscaling1 = 1.-frscaling2; 133 rangeinit = 0.; 197 134 tlimit = 1.e10*mm; 198 135 tlimitmin = 10.*tlimitminfix; 199 nstepmax = 25.;136 tgeom = 1.e50*mm; 200 137 geombig = 1.e50*mm; 201 138 geommin = 1.e-3*mm; 202 139 geomlimit = geombig; 203 140 presafety = 0.*mm; 141 142 y = 0.; 143 144 Zold = 0.; 204 145 Zeff = 1.; 146 Z2 = 1.; 147 Z23 = 1.; 148 lnZ = 0.; 149 coeffth1 = 0.; 150 coeffth2 = 0.; 151 coeffc1 = 0.; 152 coeffc2 = 0.; 153 scr1ini = fine_structure_const*fine_structure_const* 154 electron_mass_c2*electron_mass_c2/(0.885*0.885*4.*pi); 155 scr2ini = 3.76*fine_structure_const*fine_structure_const; 156 scr1 = 0.; 157 scr2 = 0.; 158 159 theta0max = pi/6.; 160 rellossmax = 0.50; 161 third = 1./3.; 205 162 particle = 0; 206 163 theManager = G4LossTableManager::Instance(); … … 218 175 219 176 void G4UrbanMscModel::Initialise(const G4ParticleDefinition* p, 220 const G4DataVector&)177 const G4DataVector&) 221 178 { 222 179 skindepth = skin*stepmin; 223 180 if(isInitialized) return; 224 225 181 // set values of some data members 226 182 SetParticle(p); 227 183 228 if (pParticleChange) 229 fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange); 230 else 231 fParticleChange = new G4ParticleChangeForMSC(); 232 233 safetyHelper = G4TransportationManager::GetTransportationManager() 234 ->GetSafetyHelper(); 235 safetyHelper->InitialiseHelper(); 184 fParticleChange = GetParticleChangeForMSC(); 185 InitialiseSafetyHelper(); 236 186 237 187 isInitialized = true; … … 385 335 if(mass > electron_mass_c2) 386 336 { 387 G4double TAU = KineticEnergy/mass ; 388 G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ; 389 G4double w = c-2. ; 390 G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ; 391 eKineticEnergy = electron_mass_c2*tau ; 392 } 393 394 G4double ChargeSquare = charge*charge; 337 G4double TAU = KineticEnergy/mass ; 338 G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ; 339 G4double w = c-2. ; 340 G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ; 341 eKineticEnergy = electron_mass_c2*tau ; 342 } 395 343 396 344 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2 ; … … 490 438 { 491 439 tPathLength = currentMinimalStep; 492 const G4DynamicParticle* dp = track.GetDynamicParticle();493 440 G4StepPoint* sp = track.GetStep()->GetPreStepPoint(); 494 441 G4StepStatus stepStatus = sp->GetStepStatus(); 442 443 const G4DynamicParticle* dp = track.GetDynamicParticle(); 495 444 496 445 if(stepStatus == fUndefined) { … … 516 465 presafety = sp->GetSafety(); 517 466 518 // G4cout << "G4UrbanMscModel::ComputeTruePathLengthLimit tPathLength= "519 // <<tPathLength<<" safety= " << presafety520 // 467 // G4cout << "G4UrbanMscModel::ComputeTruePathLengthLimit tPathLength= " 468 // <<tPathLength<<" safety= " << presafety 469 // << " range= " <<currentRange<<G4endl; 521 470 522 471 // far from geometry boundary … … 532 481 { 533 482 //compute geomlimit and presafety 534 G eomLimit(track);535 536 // is far from boundary537 if(currentRange < =presafety)483 G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange); 484 485 // is it far from boundary ? 486 if(currentRange < presafety) 538 487 { 539 488 inside = true; … … 546 495 if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined)) 547 496 { 497 rangeinit = currentRange; 548 498 if(stepStatus == fUndefined) smallstep = 1.e10; 549 499 else smallstep = 1.; 550 500 551 // facrange scaling in lambda 552 // not so strong step restriction above lambdalimit 553 G4double facr = facrange; 554 if(lambda0 > lambdalimit) 555 facr *= frscaling1+frscaling2*lambda0/lambdalimit; 556 557 // constraint from the physics 558 if (currentRange > lambda0) tlimit = facr*currentRange; 559 else tlimit = facr*lambda0; 560 561 if(tlimit > currentRange) tlimit = currentRange; 562 563 //define stepmin here (it depends on lambda!) 564 //rough estimation of lambda_elastic/lambda_transport 565 G4double rat = currentKinEnergy/MeV ; 566 rat = 1.e-3/(rat*(10.+rat)) ; 567 //stepmin ~ lambda_elastic 568 stepmin = rat*lambda0; 569 skindepth = skin*stepmin; 570 571 //define tlimitmin 572 tlimitmin = 10.*stepmin; 573 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix; 574 575 //lower limit for tlimit 576 if(tlimit < tlimitmin) tlimit = tlimitmin; 577 578 // constraint from the geometry (if tlimit above is too big) 579 G4double tgeom = geombig; 580 501 // constraint from the geometry 581 502 if((geomlimit < geombig) && (geomlimit > geommin)) 582 503 { … … 585 506 else 586 507 tgeom = 2.*geomlimit/facgeom; 587 588 if(tlimit > tgeom) tlimit = tgeom;589 508 } 590 } 591 592 //if track starts far from boundaries increase tlimit! 593 if(tlimit < facsafety*presafety) tlimit = facsafety*presafety ; 509 else 510 tgeom = geombig; 511 512 //define stepmin here (it depends on lambda!) 513 //rough estimation of lambda_elastic/lambda_transport 514 G4double rat = currentKinEnergy/MeV ; 515 rat = 1.e-3/(rat*(10.+rat)) ; 516 //stepmin ~ lambda_elastic 517 stepmin = rat*lambda0; 518 skindepth = skin*stepmin; 519 520 //define tlimitmin 521 tlimitmin = 10.*stepmin; 522 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix; 523 524 } 525 526 //step limit 527 tlimit = facrange*rangeinit; 528 if(tlimit < facsafety*presafety) 529 tlimit = facsafety*presafety; 530 531 //lower limit for tlimit 532 if(tlimit < tlimitmin) tlimit = tlimitmin; 533 534 if(tlimit > tgeom) tlimit = tgeom; 594 535 595 536 // G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit … … 597 538 598 539 // shortcut 599 if((tPathLength < tlimit) && (tPathLength < presafety)) 540 if((tPathLength < tlimit) && (tPathLength < presafety) && 541 (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth)) 600 542 return tPathLength; 601 602 G4double tnow = tlimit;603 // optimization ...604 if(geomlimit < geombig) tnow = max(tlimit,facsafety*geomlimit);605 543 606 544 // step reduction near to boundary 607 545 if(smallstep < skin) 608 546 { 609 t now= stepmin;547 tlimit = stepmin; 610 548 insideskin = true; 611 549 } … … 614 552 if(geomlimit > skindepth) 615 553 { 616 if(t now> geomlimit-0.999*skindepth)617 t now= geomlimit-0.999*skindepth;554 if(tlimit > geomlimit-0.999*skindepth) 555 tlimit = geomlimit-0.999*skindepth; 618 556 } 619 557 else 620 558 { 621 559 insideskin = true; 622 if(t now > stepmin) tnow= stepmin;560 if(tlimit > stepmin) tlimit = stepmin; 623 561 } 624 562 } 625 563 626 if(tnow < stepmin) tnow = stepmin; 627 628 if(tPathLength > tnow) tPathLength = tnow ; 564 if(tlimit < stepmin) tlimit = stepmin; 565 566 if(tPathLength > tlimit) tPathLength = tlimit ; 567 629 568 } 630 569 // for 'normal' simulation with or without magnetic field … … 635 574 // i.e. when it is needed for optimization purposes 636 575 if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 637 presafety = safetyHelper->ComputeSafety(sp->GetPosition());576 presafety = ComputeSafety(sp->GetPosition(),tPathLength); 638 577 639 578 // is far from boundary … … 645 584 646 585 if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined)) 647 { 648 // facrange scaling in lambda 649 // not so strong step restriction above lambdalimit 650 G4double facr = facrange; 651 if(lambda0 > lambdalimit) 652 facr *= frscaling1+frscaling2*lambda0/lambdalimit; 653 654 // constraint from the physics 655 if (currentRange > lambda0) tlimit = facr*currentRange; 656 else tlimit = facr*lambda0; 657 658 //lower limit for tlimit 659 tlimitmin = std::max(tlimitminfix,lambda0/nstepmax); 660 if(tlimit < tlimitmin) tlimit = tlimitmin; 661 } 662 663 //if track starts far from boundaries increase tlimit! 664 if(tlimit < facsafety*presafety) tlimit = facsafety*presafety ; 586 { 587 rangeinit = currentRange; 588 fr = facrange; 589 // 9.1 like stepping for e+/e- only (not for muons,hadrons) 590 if(mass < masslimite) 591 { 592 if(lambda0 > currentRange) 593 rangeinit = lambda0; 594 if(lambda0 > lambdalimit) 595 fr *= 0.75+0.25*lambda0/lambdalimit; 596 } 597 598 //lower limit for tlimit 599 G4double rat = currentKinEnergy/MeV ; 600 rat = 1.e-3/(rat*(10.+rat)) ; 601 tlimitmin = 10.*lambda0*rat; 602 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix; 603 } 604 //step limit 605 tlimit = fr*rangeinit; 606 607 if(tlimit < facsafety*presafety) 608 tlimit = facsafety*presafety; 609 610 //lower limit for tlimit 611 if(tlimit < tlimitmin) tlimit = tlimitmin; 665 612 666 613 if(tPathLength > tlimit) tPathLength = tlimit; … … 679 626 } 680 627 } 681 // 628 // G4cout << "tPathLength= " << tPathLength << " geomlimit= " << geomlimit 682 629 // << " currentMinimalStep= " << currentMinimalStep << G4endl; 683 684 630 return tPathLength ; 685 }686 687 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......688 689 void G4UrbanMscModel::GeomLimit(const G4Track& track)690 {691 geomlimit = geombig;692 693 // no geomlimit for the World volume694 if((track.GetVolume() != 0) &&695 (track.GetVolume() != safetyHelper->GetWorldVolume()))696 {697 G4double cstep = currentRange;698 699 geomlimit = safetyHelper->CheckNextStep(700 track.GetStep()->GetPreStepPoint()->GetPosition(),701 track.GetMomentumDirection(),702 cstep,703 presafety);704 // G4cout << "!!!G4UrbanMscModel::GeomLimit presafety= " << presafety705 // << " limit= " << geomlimit << G4endl;706 }707 631 } 708 632 … … 761 685 if(samplez) 762 686 { 763 const G4double ztmax = 0.99 , onethird = 1./3.;687 const G4double ztmax = 0.99 ; 764 688 G4double zt = zmean/tPathLength ; 765 689 … … 767 691 { 768 692 G4double u,cz1; 769 if(zt >= onethird)693 if(zt >= third) 770 694 { 771 695 G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ; … … 835 759 KineticEnergy*(KineticEnergy+2.*mass)/ 836 760 ((currentKinEnergy+mass)*(KineticEnergy+mass))); 837 G4doubley = trueStepLength/currentRadLength;761 y = trueStepLength/currentRadLength; 838 762 G4double theta0 = c_highland*std::abs(charge)*sqrt(y)/betacp; 839 763 y = log(y); 840 if(mass < masslimite) 841 theta0 *= (1.+0.051*y); 842 else if(mass < masslimitmu) 843 theta0 *= (1.+0.044*y); 844 else 845 theta0 *= (1.+0.038*y); 846 764 // correction factor from e-/proton scattering data 765 G4double corr = coeffth1+coeffth2*y; 766 if(y < -6.5) corr -= 0.011*(6.5+y); 767 theta0 *= corr ; 768 847 769 return theta0; 848 770 } … … 854 776 { 855 777 G4double kineticEnergy = dynParticle->GetKineticEnergy(); 778 856 779 if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix) || 857 (tPathLength/tausmall < lambda0) 780 (tPathLength/tausmall < lambda0)) return; 858 781 859 782 G4double cth = SampleCosineTheta(tPathLength,kineticEnergy); 783 860 784 // protection against 'bad' cth values 861 if( abs(cth) > 1.) return;785 if(std::abs(cth) > 1.) return; 862 786 863 787 G4double sth = sqrt((1.0 - cth)*(1.0 + cth)); … … 874 798 875 799 G4double r = SampleDisplacement(); 876 /*800 /* 877 801 G4cout << "G4UrbanMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy 878 802 << " sinTheta= " << sth << " r(mm)= " << r … … 880 804 << " geomStep(mm)= " << zPathLength 881 805 << G4endl; 882 */806 */ 883 807 if(r > 0.) 884 808 { … … 889 813 // compute it from the lateral correlation 890 814 G4double Phi = 0.; 891 if(std::abs(r*sth) < latcorr) {815 if(std::abs(r*sth) < latcorr) 892 816 Phi = twopi*G4UniformRand(); 893 } else { 817 else 818 { 894 819 G4double psi = std::acos(latcorr/(r*sth)); 895 if(G4UniformRand() < 0.5) Phi = phi+psi; 896 else Phi = phi-psi; 897 } 820 if(G4UniformRand() < 0.5) 821 Phi = phi+psi; 822 else 823 Phi = phi-psi; 824 } 898 825 899 826 dirx = std::cos(Phi); … … 903 830 latDirection.rotateUz(oldDirection); 904 831 905 G4ThreeVector Position = *(fParticleChange->GetProposedPosition()); 906 G4double fac = 1.; 907 if(r > safety) { 908 // ******* so safety is computed at boundary too ************ 909 G4double newsafety = safetyHelper->ComputeSafety(Position); 910 if(r > newsafety) 911 fac = newsafety/r ; 912 } 913 914 if(fac > 0.) 915 { 916 // compute new endpoint of the Step 917 G4ThreeVector newPosition = Position+fac*r*latDirection; 918 919 // definitely not on boundary 920 if(1. == fac) { 921 safetyHelper->ReLocateWithinVolume(newPosition); 922 923 } else { 924 // check safety after displacement 925 G4double postsafety = safetyHelper->ComputeSafety(newPosition); 926 927 // displacement to boundary 928 // if(postsafety < tlimitminfix) { 929 if(postsafety <= 0.0) { 930 safetyHelper->Locate(newPosition, newDirection); 931 932 // not on the boundary 933 } else { 934 safetyHelper->ReLocateWithinVolume(newPosition); 935 } 936 } 937 fParticleChange->ProposePosition(newPosition); 938 } 939 } 832 ComputeDisplacement(fParticleChange, latDirection, r, safety); 833 } 940 834 } 941 835 } … … 951 845 Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/ 952 846 couple->GetMaterial()->GetTotNbOfAtomsPerVolume() ; 847 848 if(Zold != Zeff) 849 UpdateCache(); 953 850 954 851 if(insideskin) … … 960 857 if(n > 0) 961 858 { 962 G4double tm = KineticEnergy/electron_mass_c2; 963 // ascr - screening parameter 964 G4double ascr = exp(log(Zeff)/3.)/(137.*sqrt(tm*(tm+2.))); 965 G4double ascr1 = 1.+0.5*ascr*ascr; 859 //screening (Moliere-Bethe) 860 G4double mom2 = KineticEnergy*(2.*mass+KineticEnergy); 861 G4double beta2 = mom2/((KineticEnergy+mass)*(KineticEnergy+mass)); 862 G4double ascr = scr1/mom2; 863 ascr *= 1.13+scr2/beta2; 864 G4double ascr1 = 1.+2.*ascr; 966 865 G4double bp1=ascr1+1.; 967 866 G4double bm1=ascr1-1.; 867 968 868 // single scattering from screened Rutherford x-section 969 869 G4double ct,st,phi; … … 1001 901 else if (tau >= tausmall) 1002 902 { 1003 G4double b=2.,bp1=3.,bm1=1.;1004 G4double prob = 0.;903 G4double xsi = 3.0; 904 G4double x0 = 1.; 1005 905 G4double a = 1., ea = 0., eaa = 1.; 906 G4double b=2.,b1=3.,bx=1.,eb1=3.,ebx=1.; 907 G4double prob = 1. , qprob = 1. ; 1006 908 G4double xmean1 = 1., xmean2 = 0.; 1007 909 G4double xmeanth = exp(-tau); 910 G4double x2meanth = (1.+2.*exp(-2.5*tau))/3.; 911 912 G4double relloss = 1.-KineticEnergy/currentKinEnergy; 913 if(relloss > rellossmax) 914 return SimpleScattering(xmeanth,x2meanth); 915 1008 916 G4double theta0 = ComputeTheta0(trueStepLength,KineticEnergy); 1009 917 1010 // prote xction for very small angles918 // protection for very small angles 1011 919 if(theta0 < tausmall) return cth; 1012 920 921 if(theta0 > theta0max) 922 return SimpleScattering(xmeanth,x2meanth); 1013 923 G4double sth = sin(0.5*theta0); 1014 924 a = 0.25/(sth*sth); 1015 925 1016 ea = exp(-2.*a); 1017 eaa = 1.-ea ; 1018 xmean1 = (1.+ea)/eaa-1./a; 1019 1020 G4double xmeanth = exp(-tau); 1021 b = 1./xmeanth ; 1022 bp1 = b+1.; 1023 bm1 = b-1.; 1024 // protection 1025 if(bm1 > 0.) 1026 xmean2 = b-0.5*bp1*bm1*log(bp1/bm1); 1027 else 926 ea = exp(-xsi); 927 eaa = 1.-ea ; 928 xmean1 = 1.-(1.-(1.+xsi)*ea)/(a*eaa); 929 x0 = 1.-xsi/a; 930 931 if(xmean1 <= 0.999*xmeanth) 932 return SimpleScattering(xmeanth,x2meanth); 933 934 // from MUSCAT H,Be,Fe data 935 G4double c = coeffc1; 936 if(y > -13.5) 937 c += coeffc2*exp(3.*log(y+13.5)); 938 939 if(abs(c-3.) < 0.001) c = 3.001; 940 if(abs(c-2.) < 0.001) c = 2.001; 941 942 G4double c1 = c-1.; 943 944 //from continuity of derivatives 945 b = 1.+(c-xsi)/a; 946 947 b1 = b+1.; 948 bx = c/a; 949 eb1 = exp(c1*log(b1)); 950 ebx = exp(c1*log(bx)); 951 952 xmean2 = (x0*eb1+ebx-(eb1*bx-b1*ebx)/(c-2.))/(eb1-ebx); 953 954 G4double f1x0 = a*ea/eaa; 955 G4double f2x0 = c1*eb1/(bx*(eb1-ebx)); 956 prob = f2x0/(f1x0+f2x0); 957 958 qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2); 959 960 // sampling of costheta 961 if(G4UniformRand() < qprob) 1028 962 { 1029 b = 1.+tau; 1030 bp1 = 2.+tau; 1031 bm1 = tau; 1032 xmean2 = 1.+tau*(1.-log(2./tau)); 1033 } 1034 1035 if((xmean1 >= xmeanth) && (xmean2 <= xmeanth)) 1036 { 1037 //normal case 1038 prob = (xmeanth-xmean2)/(xmean1-xmean2); 963 if(G4UniformRand() < prob) 964 cth = 1.+log(ea+G4UniformRand()*eaa)/a ; 965 else 966 cth = b-b1*bx/exp(log(ebx+(eb1-ebx)*G4UniformRand())/c1) ; 1039 967 } 1040 968 else 1041 { 1042 // x1 < xth ( x2 < xth automatically if b = 1/xth) 1043 // correct a (xmean1) 1044 if((xmeanth-xmean1)/xmeanth < 1.e-5) 1045 { 1046 // xmean1 is small probably due to precision problems 1047 xmean1 = 0.50*(1.+xmeanth); 1048 prob = (xmeanth-xmean2)/(xmean1-xmean2); 1049 } 1050 else 1051 { 1052 // correct a in order to have x1=xth 1053 G4int i=0, imax=10; 1054 do 1055 { 1056 a = 1./(1.-xmeanth+2.*ea/eaa); 1057 ea = exp(-2.*a); 1058 eaa = 1.-ea; 1059 xmean1 = (1.+ea)/eaa-1./a; 1060 i += 1; 1061 } while ((std::abs((xmeanth-xmean1)/xmeanth) > 0.05) && (i < imax)); 1062 prob = 1.; 1063 } 1064 } 1065 1066 // sampling of costheta 1067 if (G4UniformRand() < prob) 1068 cth = 1.+log(ea+G4UniformRand()*eaa)/a ; 1069 else 1070 cth = b-bp1*bm1/(bm1+2.*G4UniformRand()) ; 969 cth = -1.+2.*G4UniformRand(); 1071 970 } 1072 971 } 1073 972 return cth ; 973 } 974 975 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 976 977 G4double G4UrbanMscModel::SimpleScattering(G4double xmeanth,G4double x2meanth) 978 { 979 // 'large angle scattering' 980 // 2 model functions with correct xmean and x2mean 981 G4double a = (2.*xmeanth+9.*x2meanth-3.)/(2.*xmeanth-3.*x2meanth+1.); 982 G4double prob = (a+2.)*xmeanth/a; 983 984 // sampling 985 G4double cth = 1.; 986 if(G4UniformRand() < prob) 987 cth = -1.+2.*exp(log(G4UniformRand())/(a+1.)); 988 else 989 cth = -1.+2.*G4UniformRand(); 990 return cth; 1074 991 } 1075 992 … … 1137 1054 1138 1055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1139 1140 void G4UrbanMscModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,1141 const G4MaterialCutsCouple*,1142 const G4DynamicParticle*,1143 G4double,1144 G4double)1145 {}1146 1147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -
trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel2.cc
r1007 r1055 25 25 // 26 26 // 27 // $Id: G4UrbanMscModel2.cc,v 1. 18 2008/12/18 13:01:36 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4UrbanMscModel2.cc,v 1.26 2009/05/19 06:26:10 urban Exp $ 28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 104 104 #include "G4LossTableManager.hh" 105 105 #include "G4ParticleChangeForMSC.hh" 106 #include "G4TransportationManager.hh"107 #include "G4SafetyHelper.hh"108 106 109 107 #include "G4Poisson.hh" … … 181 179 SetParticle(p); 182 180 183 if (pParticleChange) 184 fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange); 185 else 186 fParticleChange = new G4ParticleChangeForMSC(); 187 188 safetyHelper = G4TransportationManager::GetTransportationManager() 189 ->GetSafetyHelper(); 190 safetyHelper->InitialiseHelper(); 181 fParticleChange = GetParticleChangeForMSC(); 182 InitialiseSafetyHelper(); 191 183 192 184 isInitialized = true; … … 486 478 { 487 479 //compute geomlimit and presafety 488 G eomLimit(track);480 G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange); 489 481 490 482 // is it far from boundary ? … … 499 491 500 492 if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined)) 501 493 { 502 494 rangeinit = currentRange; 503 if(stepStatus == fUndefined) smallstep = 1.e10; 504 else smallstep = 1.; 505 506 // constraint from the geometry 507 if((geomlimit < geombig) && (geomlimit > geommin)) 508 { 509 if(stepStatus == fGeomBoundary) 510 tgeom = geomlimit/facgeom; 511 else 512 tgeom = 2.*geomlimit/facgeom; 513 } 514 else 515 tgeom = geombig; 495 if(stepStatus == fUndefined) smallstep = 1.e10; 496 else smallstep = 1.; 516 497 517 498 //define stepmin here (it depends on lambda!) … … 522 503 stepmin = rat*lambda0; 523 504 skindepth = skin*stepmin; 524 525 505 //define tlimitmin 526 506 tlimitmin = 10.*stepmin; 527 507 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix; 528 508 509 // constraint from the geometry 510 if((geomlimit < geombig) && (geomlimit > geommin)) 511 { 512 // geomlimit is a geometrical step length 513 // transform it to true path length (estimation) 514 if((1.-geomlimit/lambda0) > 0.) 515 geomlimit = -lambda0*log(1.-geomlimit/lambda0)+tlimitmin ; 516 517 if(stepStatus == fGeomBoundary) 518 tgeom = geomlimit/facgeom; 519 else 520 tgeom = 2.*geomlimit/facgeom; 521 } 522 else 523 tgeom = geombig; 529 524 } 525 530 526 531 527 //step limit … … 579 575 // i.e. when it is needed for optimization purposes 580 576 if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 581 presafety = safetyHelper->ComputeSafety(sp->GetPosition());577 presafety = ComputeSafety(sp->GetPosition(),tPathLength); 582 578 583 579 // is far from boundary … … 634 630 // << " currentMinimalStep= " << currentMinimalStep << G4endl; 635 631 return tPathLength ; 636 }637 638 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......639 640 void G4UrbanMscModel2::GeomLimit(const G4Track& track)641 {642 geomlimit = geombig;643 644 // no geomlimit for the World volume645 if((track.GetVolume() != 0) &&646 (track.GetVolume() != safetyHelper->GetWorldVolume()))647 {648 G4double cstep = currentRange;649 650 geomlimit = safetyHelper->CheckNextStep(651 track.GetStep()->GetPreStepPoint()->GetPosition(),652 track.GetMomentumDirection(),653 cstep,654 presafety);655 /*656 G4cout << "!!!G4UrbanMscModel2::GeomLimit presafety= " << presafety657 << " r= " << currentRange658 << " limit= " << geomlimit << G4endl;659 */660 }661 632 } 662 633 … … 792 763 G4double theta0 = c_highland*std::abs(charge)*sqrt(y)/betacp; 793 764 y = log(y); 794 // correction factor from e- /protonscattering data765 // correction factor from e- scattering data 795 766 G4double corr = coeffth1+coeffth2*y; 796 if(y < -6.5) corr -= 0.011*(6.5+y); 767 797 768 theta0 *= corr ; 798 769 … … 860 831 latDirection.rotateUz(oldDirection); 861 832 862 G4ThreeVector Position = *(fParticleChange->GetProposedPosition()); 863 G4double fac = 1.; 864 if(r > safety) { 865 // ******* so safety is computed at boundary too ************ 866 G4double newsafety = safetyHelper->ComputeSafety(Position); 867 if(r > newsafety) 868 fac = newsafety/r ; 869 } 870 871 if(fac > 0.) 872 { 873 // compute new endpoint of the Step 874 G4ThreeVector newPosition = Position+fac*r*latDirection; 875 876 // definitely not on boundary 877 if(1. == fac) { 878 safetyHelper->ReLocateWithinVolume(newPosition); 879 880 } else { 881 // check safety after displacement 882 G4double postsafety = safetyHelper->ComputeSafety(newPosition); 883 884 // displacement to boundary 885 if(postsafety <= 0.0) { 886 safetyHelper->Locate(newPosition, newDirection); 887 888 // not on the boundary 889 } else { 890 safetyHelper->ReLocateWithinVolume(newPosition); 891 } 892 } 893 fParticleChange->ProposePosition(newPosition); 894 } 895 } 833 ComputeDisplacement(fParticleChange, latDirection, r, safety); 834 } 896 835 } 897 836 } … … 963 902 else if (tau >= tausmall) 964 903 { 965 G4double xsi = 3. 0;904 G4double xsi = 3.; 966 905 G4double x0 = 1.; 967 906 G4double a = 1., ea = 0., eaa = 1.; … … 994 933 return SimpleScattering(xmeanth,x2meanth); 995 934 996 // from MUSCAT H,Be,Fe data 997 G4double c = coeffc1; 998 if(y > -13.5) 999 c += coeffc2*exp(3.*log(y+13.5)); 935 // from e- and muon scattering data 936 G4double c = coeffc1+coeffc2*y; ; 1000 937 1001 938 if(abs(c-3.) < 0.001) c = 3.001; 1002 939 if(abs(c-2.) < 0.001) c = 2.001; 940 if(abs(c-1.) < 0.001) c = 1.001; 1003 941 1004 942 G4double c1 = c-1.; … … 1116 1054 1117 1055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1118 1119 void G4UrbanMscModel2::SampleSecondaries(std::vector<G4DynamicParticle*>*,1120 const G4MaterialCutsCouple*,1121 const G4DynamicParticle*,1122 G4double,1123 G4double)1124 {}1125 1126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -
trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel90.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4UrbanMscModel90.cc,v 1.1 0 2008/10/29 14:15:30vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4UrbanMscModel90.cc,v 1.13 2009/04/10 18:10:58 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 60 60 #include "G4LossTableManager.hh" 61 61 #include "G4ParticleChangeForMSC.hh" 62 #include "G4TransportationManager.hh"63 #include "G4SafetyHelper.hh"64 62 65 63 #include "G4Poisson.hh" … … 113 111 SetParticle(p); 114 112 115 if (pParticleChange) { 116 fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange); 117 } else { 118 fParticleChange = new G4ParticleChangeForMSC(); 119 } 120 safetyHelper = G4TransportationManager::GetTransportationManager() 121 ->GetSafetyHelper(); 122 safetyHelper->InitialiseHelper(); 113 fParticleChange = GetParticleChangeForMSC(); 114 InitialiseSafetyHelper(); 123 115 124 116 isInitialized = true; … … 419 411 { 420 412 //compute geomlimit and presafety 421 G eomLimit(track);413 G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange); 422 414 423 415 // is far from boundary … … 522 514 // i.e. when it is needed for optimization purposes 523 515 if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 524 presafety = safetyHelper->ComputeSafety(sp->GetPosition());516 presafety = ComputeSafety(sp->GetPosition(),tPathLength); 525 517 526 518 // is far from boundary … … 569 561 // << " currentMinimalStep= " << currentMinimalStep << G4endl; 570 562 return tPathLength ; 571 }572 573 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......574 575 void G4UrbanMscModel90::GeomLimit(const G4Track& track)576 {577 geomlimit = geombig;578 579 // no geomlimit for the World volume580 if((track.GetVolume() != 0) &&581 (track.GetVolume() != safetyHelper->GetWorldVolume()))582 {583 G4double cstep = currentRange;584 geomlimit = safetyHelper->CheckNextStep(585 track.GetStep()->GetPreStepPoint()->GetPosition(),586 track.GetMomentumDirection(),587 cstep,588 presafety);589 // G4cout << "!!!G4UrbanMscModel90::GeomLimit presafety= " << presafety590 // << " limit= " << geomlimit << G4endl;591 }592 563 } 593 564 … … 788 759 G4ThreeVector latDirection(dirx,diry,0.0); 789 760 latDirection.rotateUz(oldDirection); 790 791 G4ThreeVector Position = *(fParticleChange->GetProposedPosition()); 792 G4double fac = 1.; 793 if(r > safety) { 794 // ******* so safety is computed at boundary too ************ 795 G4double newsafety = safetyHelper->ComputeSafety(Position); 796 //G4double newsafety = safety; 797 if(r > newsafety) 798 fac = newsafety/r ; 799 } 800 801 if(fac > 0.) 802 { 803 // compute new endpoint of the Step 804 G4ThreeVector newPosition = Position+fac*r*latDirection; 805 806 // definetly not on boundary 807 if(1. == fac) { 808 //if(0. < fac) { 809 safetyHelper->ReLocateWithinVolume(newPosition); 810 811 812 } else { 813 // check safety after displacement 814 G4double postsafety = safetyHelper->ComputeSafety(newPosition); 815 816 // displacement to boundary 817 if(postsafety <= 0.) { 818 safetyHelper->Locate(newPosition, newDirection); 819 820 // not on the boundary 821 } else { 822 safetyHelper->ReLocateWithinVolume(newPosition); 823 } 824 } 825 fParticleChange->ProposePosition(newPosition); 826 } 827 } 761 ComputeDisplacement(fParticleChange, latDirection, r, safety); 762 } 828 763 } 829 764 } -
trunk/source/processes/electromagnetic/standard/src/G4WaterStopping.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4WaterStopping.cc,v 1.1 1 2008/12/18 13:01:38 gunterExp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4WaterStopping.cc,v 1.16 2009/05/15 19:04:21 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 29 29 //--------------------------------------------------------------------------- … … 36 36 // 37 37 // Modifications: 38 // 29.04.2009 A.Ivantchenko added revised data for G4WATER, provided by 39 // Prof.P.Sigmund Univ. Southern Denmark in the framework of 40 // the ESA Technology Research Programme (ESA contract 21435/08/NL/AT) 38 41 // 39 42 //---------------------------------------------------------------------------- … … 98 101 G4int zz[16] = {3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16,17,18}; 99 102 G4int aa[16] = {7, 9, 11, 12, 14, 16, 19, 20, 23, 24, 27, 28,31,32, 35,40}; 100 // G4double A_Ion[1 6] = {6.941,9.0122,10.811,12.011,14.007,15.999,18.998,20.180,22.990,24.305,26.982,28.086,30.974,32.065,35.453,39.948};103 // G4double A_Ion[17] = {6.941,9.0122,10.811,12.011,14.007,15.999,18.998,20.180,22.990,24.305,26.982,28.086,30.974,32.065,35.453,39.948}; 101 104 for(i=0; i<16; i++) { 102 105 Z[i] = zz[i]; … … 109 112 G4double factor = 1000.*MeV/cm; 110 113 111 G4double G4_WATER_Li [53]={2.626, 2.84, 3.191, 3.461, 3.665, 3.817, 3.927, 4.004, 4.056, 4.102, 3.998, 3.853, 3.702, 3.413, 3.158, 2.934, 2.739, 2.567, 2.415, 2.28, 1.782, 1.465, 1.247, 1.087, 0.8706, 0.7299, 0.631, 0.5575, 0.5005, 0.455, 0.4178, 0.3004, 0.2376, 0.1981, 0.1708, 0.1354, 0.1132, 0.09803, 0.08692, 0.07842, 0.07171, 0.06627, 0.0495, 0.04085, 0.03557, 0.03203, 0.0276, 0.02498, 0.02327, 0.0221, 0.02126, 0.02064, 0.02016, }; 112 AddData(E,G4_WATER_Li,factor); 113 114 G4double G4_WATER_Be [53]={3.272, 3.565, 4.061, 4.463, 4.79, 5.052, 5.258, 5.419, 5.542, 5.803, 5.787, 5.675, 5.529, 5.215, 4.912, 4.634, 4.381, 4.152, 3.944, 3.756, 3.026, 2.533, 2.179, 1.913, 1.542, 1.296, 1.122, 0.9911, 0.8898, 0.8087, 0.7423, 0.5335, 0.4219, 0.3518, 0.3034, 0.2406, 0.2013, 0.1743, 0.1545, 0.1394, 0.1275, 0.1178, 0.08805, 0.07266, 0.06328, 0.05698, 0.0491, 0.04444, 0.04141, 0.03933, 0.03783, 0.03672, 0.03588, }; 115 AddData(E,G4_WATER_Be,factor); 116 117 G4double G4_WATER_B [53]={3.773, 4.142, 4.776, 5.304, 5.749, 6.122, 6.431, 6.684, 6.89, 7.432, 7.551, 7.505, 7.391, 7.091, 6.772, 6.463, 6.172, 5.901, 5.65, 5.418, 4.484, 3.817, 3.322, 2.94, 2.392, 2.02, 1.752, 1.549, 1.391, 1.265, 1.161, 0.8332, 0.6587, 0.5492, 0.4737, 0.3757, 0.3144, 0.2723, 0.2415, 0.2179, 0.1993, 0.1842, 0.1376, 0.1136, 0.09894, 0.08909, 0.07678, 0.0695, 0.06477, 0.06151, 0.05916, 0.05743, 0.05611, }; 118 AddData(E,G4_WATER_B,factor); 119 120 G4double G4_WATER_C [53]={4.154, 4.593, 5.358, 6.009, 6.568, 7.049, 7.46, 7.809, 8.103, 8.968, 9.262, 9.311, 9.25, 8.994, 8.68, 8.358, 8.045, 7.747, 7.465, 7.199, 6.093, 5.269, 4.636, 4.137, 3.403, 2.891, 2.516, 2.23, 2.004, 1.823, 1.673, 1.2, 0.9483, 0.7904, 0.6817, 0.5406, 0.4525, 0.392, 0.3477, 0.3138, 0.287, 0.2653, 0.1983, 0.1637, 0.1426, 0.1284, 0.1107, 0.1002, 0.09335, 0.08865, 0.08528, 0.08278, 0.08088, }; 121 AddData(E,G4_WATER_C,factor); 122 123 G4double G4_WATER_N [53]={4.49, 4.984, 5.86, 6.616, 7.276, 7.854, 8.36, 8.799, 9.179, 10.39, 10.89, 11.07, 11.08, 10.9, 10.61, 10.3, 9.974, 9.66, 9.357, 9.068, 7.823, 6.859, 6.097, 5.484, 4.56, 3.9, 3.408, 3.029, 2.727, 2.482, 2.28, 1.636, 1.291, 1.076, 0.9274, 0.7354, 0.6156, 0.5333, 0.4731, 0.427, 0.3906, 0.3611, 0.27, 0.2229, 0.1942, 0.1749, 0.1507, 0.1365, 0.1272, 0.1208, 0.1162, 0.1128, 0.1102, }; 124 AddData(E,G4_WATER_N,factor); 125 126 G4double G4_WATER_O [53]={4.778, 5.321, 6.298, 7.152, 7.907, 8.578, 9.173, 9.7, 10.16, 11.73, 12.46, 12.78, 12.89, 12.81, 12.57, 12.27, 11.95, 11.63, 11.32, 11.01, 9.659, 8.571, 7.691, 6.967, 5.854, 5.042, 4.427, 3.945, 3.56, 3.245, 2.983, 2.142, 1.689, 1.406, 1.212, 0.9602, 0.8037, 0.6963, 0.6178, 0.5577, 0.5102, 0.4716, 0.3527, 0.2913, 0.2538, 0.2285, 0.197, 0.1784, 0.1663, 0.1579, 0.1519, 0.1475, 0.1441, }; 127 AddData(E,G4_WATER_O,factor); 128 129 G4double G4_WATER_F [53]={4.992, 5.575, 6.637, 7.578, 8.418, 9.171, 9.847, 10.45, 11, 12.9, 13.88, 14.35, 14.56, 14.59, 14.4, 14.13, 13.83, 13.51, 13.19, 12.87, 11.44, 10.26, 9.279, 8.463, 7.187, 6.237, 5.506, 4.928, 4.461, 4.076, 3.753, 2.707, 2.137, 1.779, 1.533, 1.215, 1.017, 0.8809, 0.7816, 0.7056, 0.6456, 0.5969, 0.4466, 0.3688, 0.3213, 0.2894, 0.2496, 0.2259, 0.2106, 0.2, 0.1924, 0.1868, 0.1825, }; 130 AddData(E,G4_WATER_F,factor); 131 132 G4double G4_WATER_Ne [53]={5.182, 5.797, 6.931, 7.948, 8.865, 9.693, 10.44, 11.12, 11.74, 13.98, 15.21, 15.85, 16.17, 16.33, 16.21, 15.98, 15.69, 15.38, 15.06, 14.74, 13.24, 11.98, 10.91, 10.01, 8.584, 7.503, 6.66, 5.986, 5.436, 4.979, 4.595, 3.332, 2.635, 2.195, 1.892, 1.499, 1.255, 1.087, 0.9646, 0.8709, 0.7969, 0.7368, 0.5514, 0.4555, 0.3969, 0.3576, 0.3083, 0.2792, 0.2602, 0.2472, 0.2378, 0.2308, 0.2255, }; 133 AddData(E,G4_WATER_Ne,factor); 134 135 G4double G4_WATER_Na [53]={5.352, 5.998, 7.203, 8.3, 9.298, 10.21, 11.04, 11.81, 12.5, 15.13, 16.68, 17.56, 18.04, 18.43, 18.44, 18.29, 18.05, 17.78, 17.48, 17.18, 15.67, 14.32, 13.15, 12.14, 10.5, 9.226, 8.218, 7.401, 6.728, 6.166, 5.69, 4.112, 3.237, 2.686, 2.307, 1.821, 1.521, 1.317, 1.168, 1.054, 0.9644, 0.8917, 0.6674, 0.5514, 0.4806, 0.4329, 0.3734, 0.3381, 0.3152, 0.2993, 0.288, 0.2796, 0.2732, }; 136 AddData(E,G4_WATER_Na,factor); 137 138 G4double G4_WATER_Mg [53]={5.542, 6.193, 7.42, 8.551, 9.59, 10.54, 11.42, 12.23, 12.98, 15.85, 17.62, 18.66, 19.26, 19.76, 19.83, 19.7, 19.47, 19.2, 18.89, 18.58, 17.02, 15.62, 14.41, 13.36, 11.64, 10.3, 9.233, 8.362, 7.64, 7.033, 6.516, 4.777, 3.792, 3.162, 2.725, 2.159, 1.806, 1.565, 1.388, 1.254, 1.147, 1.061, 0.7944, 0.6565, 0.5722, 0.5156, 0.4447, 0.4027, 0.3754, 0.3566, 0.343, 0.333, 0.3254, }; 139 AddData(E,G4_WATER_Mg,factor); 140 141 G4double G4_WATER_Al [53]={5.724, 6.39, 7.649, 8.82, 9.905, 10.91, 11.84, 12.71, 13.51, 16.66, 18.69, 19.93, 20.68, 21.38, 21.56, 21.5, 21.32, 21.07, 20.79, 20.48, 18.91, 17.47, 16.19, 15.08, 13.24, 11.78, 10.61, 9.641, 8.835, 8.153, 7.569, 5.583, 4.444, 3.71, 3.199, 2.534, 2.12, 1.836, 1.629, 1.471, 1.346, 1.245, 0.9325, 0.7707, 0.6719, 0.6054, 0.5223, 0.473, 0.441, 0.4189, 0.403, 0.3912, 0.3823, }; 142 AddData(E,G4_WATER_Al,factor); 143 144 G4double G4_WATER_Si [53]={5.905, 6.583, 7.868, 9.073, 10.2, 11.25, 12.23, 13.14, 13.99, 17.4, 19.66, 21.1, 22.01, 22.91, 23.21, 23.22, 23.09, 22.87, 22.61, 22.32, 20.76, 19.28, 17.95, 16.78, 14.83, 13.26, 11.99, 10.94, 10.06, 9.304, 8.656, 6.43, 5.135, 4.294, 3.705, 2.938, 2.458, 2.129, 1.889, 1.706, 1.561, 1.444, 1.082, 0.8942, 0.7796, 0.7026, 0.6061, 0.549, 0.5119, 0.4862, 0.4678, 0.4542, 0.4438, }; 145 AddData(E,G4_WATER_Si,factor); 146 147 G4double G4_WATER_P [53]={6.12, 6.81, 8.118, 9.352, 10.51, 11.61, 12.63, 13.58, 14.48, 18.13, 20.63, 22.28, 23.34, 24.47, 24.91, 25.02, 24.95, 24.78, 24.55, 24.28, 22.76, 21.26, 19.89, 18.67, 16.59, 14.92, 13.54, 12.39, 11.42, 10.59, 9.867, 7.367, 5.896, 4.935, 4.259, 3.376, 2.824, 2.445, 2.169, 1.959, 1.792, 1.657, 1.242, 1.027, 0.8954, 0.807, 0.6963, 0.6308, 0.5881, 0.5587, 0.5375, 0.5219, 0.51, }; 148 AddData(E,G4_WATER_P,factor); 149 150 G4double G4_WATER_S [53]={6.294, 7, 8.338, 9.604, 10.8, 11.94, 13, 14, 14.94, 18.82, 21.55, 23.41, 24.65, 26.01, 26.6, 26.81, 26.81, 26.69, 26.5, 26.26, 24.79, 23.28, 21.88, 20.61, 18.43, 16.64, 15.16, 13.92, 12.86, 11.95, 11.15, 8.371, 6.715, 5.624, 4.856, 3.847, 3.217, 2.785, 2.47, 2.229, 2.04, 1.886, 1.413, 1.169, 1.019, 0.9187, 0.7929, 0.7183, 0.6697, 0.6362, 0.6122, 0.5944, 0.5808, }; 151 AddData(E,G4_WATER_S,factor); 152 153 G4double G4_WATER_Cl [53]={6.522, 7.237, 8.59, 9.875, 11.1, 12.26, 13.36, 14.39, 15.37, 19.45, 22.4, 24.45, 25.86, 27.46, 28.19, 28.5, 28.57, 28.5, 28.34, 28.13, 26.72, 25.21, 23.78, 22.47, 20.2, 18.32, 16.75, 15.42, 14.28, 13.3, 12.44, 9.392, 7.557, 6.34, 5.479, 4.344, 3.633, 3.145, 2.789, 2.517, 2.303, 2.13, 1.596, 1.32, 1.151, 1.038, 0.8957, 0.8115, 0.7567, 0.7189, 0.6917, 0.6717, 0.6563, }; 154 AddData(E,G4_WATER_Cl,factor); 155 156 G4double G4_WATER_Ar [53]={6.642, 7.369, 8.739, 10.04, 11.28, 12.47, 13.59, 14.66, 15.66, 19.93, 23.08, 25.32, 26.89, 28.72, 29.6, 30.01, 30.15, 30.13, 30, 29.82, 28.46, 26.94, 25.49, 24.15, 21.81, 19.86, 18.22, 16.82, 15.61, 14.57, 13.65, 10.39, 8.394, 7.063, 6.114, 4.859, 4.067, 3.522, 3.125, 2.821, 2.581, 2.387, 1.789, 1.48, 1.291, 1.164, 1.005, 0.9105, 0.8491, 0.8067, 0.7763, 0.7537, 0.7365, }; 157 AddData(E,G4_WATER_Ar,factor); 114 G4double G4_WATER_Li[53]={2.3193,2.5198,2.8539,3.1164,3.3203,3.4756,3.5914,3.6755,3.7347,3.8125,3.7349,3.6134,3.4818,3.2258,2.9949,2.7909,2.611,2.4517,2.3103,2.1841,1.7151,1.4139,1.2053,1.0525,0.84417,0.70862,0.61317,0.54214,0.48708,0.44305,0.40697,0.29312,0.23208,0.19364,0.16706,0.13252,0.11092,0.09608,0.08522,0.076915,0.07035,0.065026,0.048615,0.040137,0.034964,0.03149,0.027148,0.024579,0.022911,0.021761,0.020937,0.020327,0.019862}; 115 AddData(E,G4_WATER_Li,factor); 116 G4double G4_WATER_Be[53]={2.872,3.1439,3.6102,3.9967,4.3169,4.5788,4.7897,4.9568,5.0872,5.387,5.4028,5.3185,5.1971,4.9243,4.6549,4.4036,4.1732,3.9629,3.771,3.5957,2.9117,2.4439,2.1062,1.8518,1.4956,1.2587,1.0901,0.96393,0.86589,0.78742,0.72313,0.52053,0.41214,0.34394,0.2968,0.2355,0.19717,0.17081,0.15152,0.13676,0.1251,0.11564,0.086471,0.071399,0.062202,0.056023,0.048303,0.043734,0.040767,0.038723,0.037256,0.036172,0.035345}; 117 AddData(E,G4_WATER_Be,factor); 118 G4double G4_WATER_B[53]={3.2922,3.6315,4.2226,4.7242,5.1543,5.5219,5.8323,6.0911,6.3043,6.8888,7.0451,7.0309,6.9445,6.6925,6.4129,6.1372,5.8747,5.6283,5.3983,5.1841,4.3121,3.6826,3.2109,2.8459,2.3203,1.9619,1.7028,1.5072,1.3543,1.2315,1.1307,0.81305,0.64344,0.53693,0.46337,0.36777,0.30797,0.26684,0.23674,0.21371,0.1955,0.18072,0.13517,0.11163,0.097256,0.087601,0.075535,0.068395,0.063757,0.060561,0.058268,0.056574,0.05528}; 119 AddData(E,G4_WATER_B,factor); 120 G4double G4_WATER_C[53]={3.6037,4.0035,4.7125,5.3248,5.8601,6.3287,6.7363,7.0875,7.3872,8.2986,8.635,8.7189,8.6879,8.485,8.2162,7.9331,7.6534,7.3839,7.1272,6.884,5.8573,5.0814,4.4808,4.0044,3.3005,2.808,2.4458,2.1691,1.9511,1.7751,1.6302,1.1714,0.9263,0.77269,0.66676,0.52925,0.44328,0.38415,0.34086,0.30773,0.28154,0.26028,0.19473,0.16083,0.14014,0.12624,0.10886,0.098575,0.091894,0.08729,0.083986,0.081545,0.079682}; 121 AddData(E,G4_WATER_C,factor); 122 G4double G4_WATER_N[53]={3.8821,4.3278,5.1314,5.8371,6.4632,7.021,7.5168,7.9545,8.3378,9.5943,10.145,10.356,10.402,10.278,10.041,9.7677,9.4845,9.2035,8.9301,8.6668,7.5173,6.6126,5.8919,5.308,4.4226,3.7883,3.3138,2.9467,2.655,2.4179,2.2217,1.5965,1.2614,1.0516,0.90715,0.71995,0.60305,0.52268,0.46384,0.41882,0.3832,0.35431,0.26515,0.21903,0.19087,0.17194,0.14829,0.13429,0.12519,0.11892,0.11442,0.1111,0.10856}; 123 AddData(E,G4_WATER_N,factor); 124 G4double G4_WATER_O[53]={4.1215,4.6063,5.4945,6.2868,6.9978,7.6391,8.2175,8.7372,9.201,10.808,11.596,11.955,12.096,12.077,11.89,11.639,11.364,11.081,10.799,10.523,9.2787,8.2615,7.4307,6.7431,5.6787,4.8981,4.3045,3.8393,3.4663,3.161,2.9069,2.0903,1.6501,1.3745,1.1851,0.94004,0.78733,0.68244,0.60568,0.54694,0.50048,0.46278,0.34643,0.28622,0.24945,0.22474,0.19384,0.17555,0.16367,0.15547,0.14959,0.14525,0.14194}; 125 AddData(E,G4_WATER_O,factor); 126 G4double G4_WATER_F[53]={4.2951,4.8107,5.7678,6.6342,7.4196,8.1343,8.7858,9.3786,9.9152,11.857,12.89,13.408,13.652,13.749,13.62,13.398,13.136,12.857,12.573,12.291,10.982,9.88,8.9601,8.1871,6.9687,6.0574,5.3535,4.7954,4.3434,3.9704,3.658,2.642,2.0878,1.7393,1.4994,1.1892,0.99601,0.86336,0.76632,0.69207,0.63334,0.58568,0.43857,0.36242,0.3159,0.28462,0.24552,0.22237,0.20733,0.19696,0.18951,0.18401,0.17981}; 127 AddData(E,G4_WATER_F,factor); 128 G4double G4_WATER_Ne[53]={4.4513,4.991,6.004,6.9338,7.7852,8.5662,9.284,9.9435,10.547,12.813,14.099,14.791,15.151,15.382,15.321,15.136,14.895,14.626,14.345,14.061,12.705,11.53,10.532,9.6823,8.3208,7.2846,6.4735,5.8234,5.2919,4.8502,4.4779,3.252,2.5745,2.146,1.8503,1.4675,1.2291,1.0654,0.94575,0.85419,0.78178,0.72301,0.54158,0.44763,0.39022,0.35162,0.30335,0.27477,0.25619,0.24339,0.23419,0.2274,0.22222}; 129 AddData(E,G4_WATER_Ne,factor); 130 G4double G4_WATER_Na[53]={4.5914,5.1553,6.2244,7.2203,8.1435,8.9982,9.7906,10.526,11.206,13.848,15.453,16.384,16.916,17.369,17.442,17.344,17.16,16.93,16.675,16.407,15.045,13.799,12.706,11.753,10.187,8.9672,7.9956,7.2072,6.5562,6.0112,5.5493,4.0154,3.1635,2.6262,2.2573,1.7832,1.4904,1.2907,1.1451,1.0339,0.94615,0.87498,0.65548,0.54186,0.47243,0.42574,0.36734,0.33275,0.31027,0.29478,0.28364,0.27542,0.26915}; 131 AddData(E,G4_WATER_Na,factor); 132 G4double G4_WATER_Mg[53]={4.7537,5.3178,6.3962,7.4137,8.3663,9.2554,10.085,10.859,11.581,14.455,16.279,17.373,18.018,18.598,18.727,18.654,18.479,18.25,17.99,17.716,16.313,15.026,13.895,12.907,11.277,9.9981,8.9727,8.1344,7.4376,6.8507,6.35,4.6625,3.7049,3.0916,2.666,2.1135,1.7695,1.5336,1.3613,1.2296,1.1255,1.041,0.7802,0.64511,0.56252,0.50698,0.43749,0.39634,0.36958,0.35114,0.33789,0.3281,0.32063}; 133 AddData(E,G4_WATER_Mg,factor); 134 G4double G4_WATER_Al[53]={4.911,5.4856,6.5852,7.6302,8.6193,9.551,10.426,11.248,12.018,15.152,17.23,18.531,19.33,20.11,20.354,20.352,20.224,20.024,19.785,19.521,18.12,16.795,15.611,14.565,12.819,11.431,10.305,9.3772,8.6003,7.9413,7.3761,5.45,4.342,3.6272,3.129,2.4808,2.0767,1.7997,1.5975,1.4429,1.3207,1.2216,0.91581,0.75739,0.66053,0.59537,0.51383,0.46554,0.43413,0.41248,0.39693,0.38544,0.37667}; 135 AddData(E,G4_WATER_Al,factor); 136 G4double G4_WATER_Si[53]={5.0697,5.6529,6.77,7.8376,8.8564,9.8229,10.736,11.597,12.409,15.773,18.09,19.594,20.549,21.535,21.9,21.975,21.896,21.73,21.513,21.265,19.878,18.525,17.298,16.202,14.354,12.867,11.651,10.64,9.7876,9.0607,8.4341,6.2756,5.0169,4.1981,3.6247,2.8758,2.4078,2.0868,1.8523,1.6731,1.5315,1.4167,1.0623,0.8787,0.76644,0.69091,0.59637,0.54036,0.50394,0.47883,0.46078,0.44746,0.43728}; 137 AddData(E,G4_WATER_Si,factor); 138 G4double G4_WATER_P[53]={5.2616,5.8538,6.9867,8.074,9.1192,10.117,11.065,11.964,12.815,16.396,18.945,20.657,21.779,22.993,23.502,23.672,23.659,23.54,23.356,23.132,21.792,20.425,19.162,18.02,16.066,14.472,13.155,12.051,11.114,10.311,9.6149,7.1917,5.7614,4.8249,4.1667,3.305,2.7661,2.3966,2.1269,1.9215,1.7603,1.6263,1.2197,1.0091,0.88027,0.79361,0.68512,0.62083,0.57902,0.55019,0.52947,0.51417,0.50248}; 139 AddData(E,G4_WATER_P,factor); 140 G4double G4_WATER_S[53]={5.4129,6.0193,7.1761,8.2871,9.36,10.39,11.373,12.308,13.196,16.986,19.762,21.683,22.976,24.431,25.091,25.363,25.421,25.354,25.208,25.013,23.734,22.366,21.074,19.891,17.84,16.146,14.731,13.536,12.516,11.636,10.869,8.1727,6.5617,5.4997,4.7505,3.7671,3.1514,2.7293,2.4215,2.1866,2.0012,1.8509,1.3881,1.1485,1.002,0.90348,0.78008,0.70695,0.65938,0.62657,0.60299,0.58558,0.57229}; 141 AddData(E,G4_WATER_S,factor); 142 G4double G4_WATER_Cl[53]={5.6171,6.2307,7.3984,8.5209,9.6097,10.661,11.669,12.632,13.551,17.518,20.497,22.615,24.076,25.769,26.58,26.953,27.082,27.066,26.958,26.791,25.579,24.214,22.901,21.685,19.553,17.771,16.271,14.996,13.9,12.949,12.118,9.1688,7.385,6.2003,5.3604,4.2535,3.5588,3.0821,2.7343,2.4689,2.2595,2.0898,1.5673,1.297,1.1318,1.0205,0.88128,0.79874,0.74504,0.708,0.68138,0.66172,0.64671}; 143 AddData(E,G4_WATER_Cl,factor); 144 G4double G4_WATER_Ar[53]={5.7158,6.3394,7.5204,8.6525,9.7528,10.82,11.849,12.836,13.78,17.904,21.07,23.375,24.999,26.928,27.889,28.361,28.559,28.592,28.521,28.381,27.228,25.869,24.541,23.299,21.103,19.255,17.686,16.346,15.19,14.183,13.3,10.137,8.2021,6.9062,5.9819,4.757,3.9841,3.4523,3.0636,2.7668,2.5324,2.3425,1.7573,1.4546,1.2694,1.1448,0.98872,0.8962,0.83601,0.79448,0.76464,0.74259,0.72575}; 145 AddData(E,G4_WATER_Ar,factor); 158 146 159 147 if(corr) { -
trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4WentzelVIModel.cc,v 1. 16 2008/11/19 11:47:50vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4WentzelVIModel.cc,v 1.32 2009/05/10 16:09:29 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 52 52 // 53 53 54 55 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 56 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... … … 60 59 #include "G4LossTableManager.hh" 61 60 #include "G4ParticleChangeForMSC.hh" 62 #include "G4TransportationManager.hh"63 #include "G4SafetyHelper.hh"64 61 #include "G4PhysicsTableHelper.hh" 65 62 #include "G4ElementVector.hh" … … 71 68 72 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 70 71 G4double G4WentzelVIModel::ScreenRSquare[] = {0.0}; 72 G4double G4WentzelVIModel::FormFactor[] = {0.0}; 73 73 74 74 using namespace std; … … 96 96 thePositron = G4Positron::Positron(); 97 97 theProton = G4Proton::Proton(); 98 a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885);98 lowEnergyLimit = 0.1*keV; 99 99 G4double p0 = electron_mass_c2*classic_electr_radius; 100 100 coeff = twopi*p0*p0; 101 constn = 6.937e-6/(MeV*MeV);102 101 tkin = targetZ = mom2 = DBL_MIN; 103 102 ecut = etag = DBL_MAX; … … 106 105 xsecn.resize(nelments); 107 106 prob.resize(nelments); 108 for(size_t j=0; j<100; j++) { 109 FF[j] = 0.0; 110 } 107 108 // Thomas-Fermi screening radii 109 // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282 110 111 if(0.0 == ScreenRSquare[0]) { 112 G4double a0 = electron_mass_c2/0.88534; 113 G4double constn = 6.937e-6/(MeV*MeV); 114 115 ScreenRSquare[0] = alpha2*a0*a0; 116 for(G4int j=1; j<100; j++) { 117 G4double x = a0*fNistManager->GetZ13(j); 118 ScreenRSquare[j] = alpha2*x*x; 119 x = fNistManager->GetA27(j); 120 FormFactor[j] = constn*x*x; 121 } 122 } 111 123 } 112 124 … … 132 144 if(!isInitialized) { 133 145 isInitialized = true; 134 135 if (pParticleChange) 136 fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange); 137 else 138 fParticleChange = new G4ParticleChangeForMSC(); 139 140 safetyHelper = G4TransportationManager::GetTransportationManager() 141 ->GetSafetyHelper(); 142 safetyHelper->InitialiseHelper(); 146 fParticleChange = GetParticleChangeForMSC(); 147 InitialiseSafetyHelper(); 143 148 } 144 149 } … … 156 161 SetupKinematic(ekin, cutEnergy); 157 162 SetupTarget(Z, ekin); 158 G4double xsec = ComputeTransportXSectionPer Volume();163 G4double xsec = ComputeTransportXSectionPerAtom(); 159 164 /* 160 165 G4cout << "CS: e= " << tkin << " cosEl= " << cosTetMaxElec2 … … 167 172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 168 173 169 G4double G4WentzelVIModel::ComputeTransportXSectionPer Volume()174 G4double G4WentzelVIModel::ComputeTransportXSectionPerAtom() 170 175 { 171 176 G4double xSection = 0.0; … … 189 194 y = 0.0; 190 195 } 191 xSection += y/targetZ; 192 } 193 /* 194 G4cout << "G4WentzelVIModel:XS per A " << " Z= " << Z << " e(MeV)= " << kinEnergy/MeV 196 xSection = y; 197 } 198 /* 199 G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ 200 << " e(MeV)= " << tkin/MeV << " XSel= " << xSection 195 201 << " cut(MeV)= " << ecut/MeV 196 202 << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ 197 << " zmaxN= " << (1.0 - cosTetMsxNuc)/screenZ << G4endl; 203 << " zmaxN= " << (1.0 - cosTetMaxNuc2)/screenZ 204 << " costm= " << cosTetMaxNuc2 << G4endl; 198 205 */ 199 200 206 // scattering off nucleus 201 207 if(cosTetMaxNuc2 < 1.0) { 202 208 x = 1.0 - cosTetMaxNuc2; 203 209 x1 = screenZ*formfactA; 204 x2 = 1.0 /(1.0 - x1);210 x2 = 1.0 - x1; 205 211 x3 = x/screenZ; 206 212 x4 = formfactA*x; 207 213 // low-energy limit 208 214 if(x3 < numlimit && x1 < numlimit) { 209 y = 0.5*x3*x3* x2*x2*x2*(1.0 - 1.333333*x3 + 1.5*x3*x3 - 1.5*x1210 + 3.0*x1*x1 + 2.666666*x3*x1);215 y = 0.5*x3*x3*(1.0 - 1.3333333*x3 + 1.5*x3*x3 - 1.5*x1 216 + 3.0*x1*x1 + 2.666666*x3*x1)/(x2*x2*x2); 211 217 // high energy limit 212 } else if( 1.0 < x1) {218 } else if(x2 <= 0.0) { 213 219 x4 = x1*(1.0 + x3); 214 220 y = x3*(1.0 + 0.5*x3 - (2.0 - x1)*(1.0 + x3 + x3*x3/3.0)/x4)/(x4*x4); … … 216 222 } else { 217 223 y = ((1.0 + x1)*x2*log((1. + x3)/(1. + x4)) 218 - x3/(1. + x3) - x4/(1. + x4))*x2*x2; 219 } 224 - x3/(1. + x3) - x4/(1. + x4))/(x2*x2); 225 } 226 //G4cout << "y= " << y << " x1= " <<x1<<" x2= " <<x2 227 // <<" x3= "<<x3<<" x4= " << x4<<G4endl; 220 228 if(y < 0.0) { 221 229 nwarnings++; … … 230 238 y = 0.0; 231 239 } 232 xSection += y; 233 } 234 xSection *= (coeff*targetZ*targetZ*chargeSquare*invbeta2/mom2); 235 // G4cout << " XStotal= " << xSection/barn << " screenZ= " << screenZ 236 // << " formF= " << formfactA << " for " << p->GetParticleName() << G4endl; 240 xSection += y*targetZ; 241 } 242 xSection *= kinFactor; 243 /* 244 G4cout << "Z= " << targetZ << " XStot= " << xSection/barn 245 << " screenZ= " << screenZ << " formF= " << formfactA 246 << " for " << particle->GetParticleName() 247 << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2) 248 << " x= " << x 249 << G4endl; 250 */ 237 251 return xSection; 238 252 } … … 277 291 // i.e. when it is needed for optimization purposes 278 292 if(stepStatus != fGeomBoundary && presafety < tlimitminfix) 279 presafety = safetyHelper->ComputeSafety(sp->GetPosition());293 presafety = ComputeSafety(sp->GetPosition(), tlimit); 280 294 /* 281 295 G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit tlimit= " … … 291 305 G4double rlimit = facrange*lambda0; 292 306 G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1); 293 if(rcut > rlimit) rlimit = std::pow(2.0*rcut*rcut*lambda0,0.33333333); 307 if(rcut > rlimit) rlimit = std::pow(rcut*rcut*rlimit,0.33333333); 308 rlimit = std::min(rlimit, facgeom*currentMaterial->GetRadlen()); 294 309 if(rlimit < tlimit) tlimit = rlimit; 295 310 } … … 404 419 } else { 405 420 406 // define threshold angle as 2 sigma of central value421 // define threshold angle between single and multiple scattering 407 422 cosThetaMin = 1.0 - 3.0*x1; 408 423 … … 570 585 571 586 if(r > tlimitminfix) { 572 G4ThreeVector Position = *(fParticleChange->GetProposedPosition()); 573 G4double fac= 1.; 574 if(r >= safety) { 575 // ******* so safety is computed at boundary too ************ 576 G4double newsafety = 577 safetyHelper->ComputeSafety(Position) - tlimitminfix; 578 if(newsafety <= 0.0) fac = 0.0; 579 else if(r > newsafety) fac = newsafety/r ; 580 //G4cout << "NewSafety= " << newsafety << " fac= " << fac 581 // << " r= " << r << " sint= " << sint << " pos " << Position << G4endl; 582 } 583 584 if(fac > 0.) { 585 // compute new endpoint of the Step 586 G4ThreeVector newPosition = Position + fac*pos; 587 588 // check safety after displacement 589 G4double postsafety = safetyHelper->ComputeSafety(newPosition); 590 591 // displacement to boundary 592 if(postsafety <= 0.0) { 593 safetyHelper->Locate(newPosition, newDirection); 594 595 // not on the boundary 596 } else { 597 safetyHelper->ReLocateWithinVolume(newPosition); 598 // if(fac < 1.0) G4cout << "NewPosition " << newPosition << G4endl; 599 } 600 601 fParticleChange->ProposePosition(newPosition); 602 } 587 pos /= r; 588 ComputeDisplacement(fParticleChange, pos, r, safety); 603 589 } 604 590 } … … 624 610 G4double xs = 0.0; 625 611 626 G4double fac = coeff*chargeSquare*invbeta2/mom2;627 628 612 for (G4int i=0; i<nelm; i++) { 629 613 SetupTarget((*theElementVector)[i]->GetZ(), tkin); … … 635 619 cosTetMaxNuc2 = std::max(cosnm,cosThetaMin); 636 620 cosTetMaxElec2 = std::max(cosem,cosThetaMin); 637 xtsec += ComputeTransportXSectionPer Volume()*density;621 xtsec += ComputeTransportXSectionPerAtom()*density; 638 622 // return limit back 639 623 cosTetMaxElec2 = cosem; … … 643 627 G4double nsec = 0.0; 644 628 G4double x1 = 1.0 - cosThetaMin + screenZ; 645 G4double f = fac*targetZ*density;629 G4double f = kinFactor*density; 646 630 647 631 // scattering off electrons … … 656 640 G4double s = screenZ*formfactA; 657 641 G4double z1 = 1.0 - cosnm + screenZ; 658 G4double d = (1.0 - s)/formfactA; 642 G4double s1 = 1.0 - s; 643 G4double d = s1/formfactA; 659 644 660 645 // check numerical limit … … 667 652 G4double x2 = x1 + d; 668 653 G4double z2 = z1 + d; 669 nsec = (1.0 + 2.0*s)*((cosThetaMin - cosnm)*(1.0/(x1*z1) + 1.0/(x2*z2)) - 670 2.0*log(z1*x2/(z2*x1))/d); 654 nsec = (1.0/x1 - 1.0/z1 + 1.0/x2 - 1.0/z2 - 2.0*log(z1*x2/(z2*x1))/d)/(s1*s1); 671 655 } 672 656 nsec *= f*targetZ; … … 801 785 G4double mom22 = t1*(t1 + 2.0*mass); 802 786 G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22); 803 if(ctm < 1.0) cosTetMaxElec = ctm; 787 if(ctm < 1.0) cosTetMaxElec = ctm; 788 if(ctm < -1.0) cosTetMaxElec = -1.0; 804 789 } 805 790 } … … 808 793 809 794 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 810 811 void G4WentzelVIModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,812 const G4MaterialCutsCouple*,813 const G4DynamicParticle*,814 G4double,815 G4double)816 {}817 818 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -
trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlung.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eBremsstrahlung.cc,v 1.5 5 2008/11/14 19:23:07 vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eBremsstrahlung.cc,v 1.56 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 104 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 105 105 106 G4bool G4eBremsstrahlung::IsApplicable(const G4ParticleDefinition& p) 107 { 108 return (&p == G4Electron::Electron() || &p == G4Positron::Positron()); 109 } 110 111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 112 106 113 void G4eBremsstrahlung::InitialiseEnergyLossProcess( 107 114 const G4ParticleDefinition* p, -
trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eBremsstrahlungModel.cc,v 1.4 3 2008/11/13 19:28:58 vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eBremsstrahlungModel.cc,v 1.44 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 160 160 } 161 161 if(isInitialised) return; 162 163 if(pParticleChange) { 164 fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange); 165 } else { 166 fParticleChange = new G4ParticleChangeForLoss(); 167 } 162 fParticleChange = GetParticleChangeForLoss(); 168 163 isInitialised = true; 169 164 } -
trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eBremsstrahlungRelModel.cc,v 1.1 2 2008/11/13 23:28:27 schaelicExp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eBremsstrahlungRelModel.cc,v 1.14 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 178 178 179 179 if(isInitialised) return; 180 181 if(pParticleChange) { 182 fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange); 183 } else { 184 fParticleChange = new G4ParticleChangeForLoss(); 185 } 180 fParticleChange = GetParticleChangeForLoss(); 186 181 isInitialised = true; 187 182 } -
trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eCoulombScatteringModel.cc,v 1. 59 2008/10/22 18:39:29 vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eCoulombScatteringModel.cc,v 1.69 2009/05/10 16:09:29 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 66 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 67 67 68 G4double G4eCoulombScatteringModel::ScreenRSquare[] = {0.0}; 69 G4double G4eCoulombScatteringModel::FormFactor[] = {0.0}; 70 68 71 using namespace std; 69 72 … … 84 87 currentMaterial = 0; 85 88 currentElement = 0; 86 a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885);89 lowEnergyLimit = keV; 87 90 G4double p0 = electron_mass_c2*classic_electr_radius; 88 91 coeff = twopi*p0*p0; 89 constn = 6.937e-6/(MeV*MeV);90 92 tkin = targetZ = mom2 = DBL_MIN; 91 93 elecXSection = nucXSection = 0.0; 92 recoilThreshold = DBL_MAX;94 recoilThreshold = 100.*keV; 93 95 ecut = DBL_MAX; 94 96 particle = 0; 95 97 currentCouple = 0; 96 for(size_t j=0; j<100; j++) { 97 FF[j] = 0.0; 98 } 98 99 // Thomas-Fermi screening radii 100 // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282 101 102 if(0.0 == ScreenRSquare[0]) { 103 G4double a0 = electron_mass_c2/0.88534; 104 G4double constn = 6.937e-6/(MeV*MeV); 105 106 ScreenRSquare[0] = alpha2*a0*a0; 107 for(G4int j=1; j<100; j++) { 108 G4double x = a0*fNistManager->GetZ13(j); 109 ScreenRSquare[j] = alpha2*x*x; 110 x = fNistManager->GetA27(j); 111 FormFactor[j] = constn*x*x; 112 } 113 } 99 114 } 100 115 … … 121 136 if(!isInitialised) { 122 137 isInitialised = true; 123 124 if(pParticleChange) 125 fParticleChange = 126 reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange); 127 else 128 fParticleChange = new G4ParticleChangeForGamma(); 138 fParticleChange = GetParticleChangeForGamma(); 129 139 } 130 140 if(mass < GeV && particle->GetParticleType() != "nucleus") { … … 157 167 G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22); 158 168 //G4cout << "ctm= " << ctm << G4endl; 159 if(ctm < 1.0) cosTetMaxElec = ctm; 169 if(ctm < 1.0) cosTetMaxElec = ctm; 170 if(ctm < -1.0) cosTetMaxElec = -1.0; 160 171 } 161 172 } … … 196 207 { 197 208 // This method needs initialisation before be called 198 199 G4double fac = coeff*targetZ*chargeSquare*invbeta2/mom2; 209 G4double fac = coeff*targetZ*chargeSquare*kinFactor; 200 210 elecXSection = 0.0; 201 211 nucXSection = 0.0; … … 216 226 G4double s = screenZ*formfactA; 217 227 G4double z1 = 1.0 - cosTetMaxNuc2 + screenZ; 218 G4double d = (1.0 - s)/formfactA; 228 G4double s1 = 1.0 - s; 229 G4double d = s1/formfactA; 219 230 //G4cout <<"x1= "<<x1<<" z1= " <<z1<<" s= "<<s << " d= " <<d <<G4endl; 220 231 if(d < 0.2*x1) { … … 226 237 G4double x2 = x1 + d; 227 238 G4double z2 = z1 + d; 228 x = (1.0 + 2.0*s)*((cosTetMinNuc - cosTetMaxNuc2)*(1.0/(x1*z1) + 1.0/(x2*z2)) - 229 2.0*log(z1*x2/(z2*x1))/d); 239 x = (1.0/x1 - 1.0/z1 + 1.0/x2 - 1.0/z2 - 2.0*log(z1*x2/(z2*x1))/d)/(s1*s1); 230 240 } 231 241 nucXSection += fac*targetZ*x; 232 242 } 233 234 243 //G4cout<<" cross(bn)= "<<nucXSection/barn<<" xsElec(bn)= "<<elecXSection/barn 235 244 // << " Asc= " << screenZ << G4endl; … … 283 292 G4int ia = SelectIsotopeNumber(currentElement); 284 293 G4double Trec = z1*mom2/(amu_c2*G4double(ia)); 285 G4double th = 286 std::min(recoilThreshold, 287 targetZ*currentElement->GetIonisation()->GetMeanExcitationEnergy()); 288 289 if(Trec > th) { 290 G4int iz = G4int(targetZ); 294 295 if(Trec > recoilThreshold) { 291 296 G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz); 292 297 Trec = z1*mom2/ion->GetPDGMass(); -
trunk/source/processes/electromagnetic/standard/src/G4eIonisation.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eIonisation.cc,v 1.5 6 2008/10/20 08:56:41vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eIonisation.cc,v 1.57 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 103 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 104 104 105 G4double G4eIonisation::MinPrimaryEnergy(const G4ParticleDefinition*, 106 const G4Material*, 107 G4double cut) 108 { 109 G4double x = cut; 110 if(isElectron) x += cut; 111 return x; 112 } 113 114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 115 116 G4bool G4eIonisation::IsApplicable(const G4ParticleDefinition& p) 117 { 118 return (&p == G4Electron::Electron() || &p == G4Positron::Positron()); 119 } 120 121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 122 105 123 void G4eIonisation::InitialiseEnergyLossProcess( 106 124 const G4ParticleDefinition* part, -
trunk/source/processes/electromagnetic/standard/src/G4eeToTwoGammaModel.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eeToTwoGammaModel.cc,v 1.1 4 2007/05/23 08:47:35vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eeToTwoGammaModel.cc,v 1.15 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 102 102 { 103 103 if(isInitialised) return; 104 105 if(pParticleChange) 106 fParticleChange = 107 reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange); 108 else 109 fParticleChange = new G4ParticleChangeForGamma(); 110 104 fParticleChange = GetParticleChangeForGamma(); 111 105 isInitialised = true; 112 106 } -
trunk/source/processes/electromagnetic/standard/src/G4eplusAnnihilation.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eplusAnnihilation.cc,v 1. 29 2008/10/15 17:53:44vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eplusAnnihilation.cc,v 1.30 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 79 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 80 80 81 G4bool G4eplusAnnihilation::IsApplicable(const G4ParticleDefinition& p) 82 { 83 return (&p == G4Positron::Positron()); 84 } 85 86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 87 88 G4double G4eplusAnnihilation::AtRestGetPhysicalInteractionLength( 89 const G4Track&, G4ForceCondition* condition) 90 { 91 *condition = NotForced; 92 return 0.0; 93 } 94 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 96 81 97 void G4eplusAnnihilation::InitialiseProcess(const G4ParticleDefinition*) 82 98 { -
trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4hIonisation.cc,v 1.8 1 2008/10/22 16:02:20vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4hIonisation.cc,v 1.82 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 122 122 {} 123 123 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 125 126 G4bool G4hIonisation::IsApplicable(const G4ParticleDefinition& p) 127 { 128 return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV && 129 !p.IsShortLived()); 130 } 131 132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 133 134 G4double G4hIonisation::MinPrimaryEnergy(const G4ParticleDefinition*, 135 const G4Material*, 136 G4double cut) 137 { 138 G4double x = 0.5*cut/electron_mass_c2; 139 G4double g = x*ratio + std::sqrt((1. + x)*(1. + x*ratio*ratio)); 140 return mass*(g - 1.0); 141 } 142 124 143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 125 144 -
trunk/source/processes/electromagnetic/standard/src/G4ionIonisation.cc
r1007 r1055 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4ionIonisation.cc,v 1.6 5 2008/10/15 17:53:44vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4ionIonisation.cc,v 1.66 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 104 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 105 105 106 G4bool G4ionIonisation::IsApplicable(const G4ParticleDefinition& p) 107 { 108 return (p.GetPDGCharge() != 0.0 && !p.IsShortLived() && 109 p.GetParticleType() == "nucleus"); 110 } 111 112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 113 114 G4double G4ionIonisation::MinPrimaryEnergy(const G4ParticleDefinition* p, 115 const G4Material*, 116 G4double cut) 117 { 118 return 119 p->GetPDGMass()*(std::sqrt(1. + 0.5*cut/CLHEP::electron_mass_c2) - 1.0); 120 } 121 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 123 106 124 void G4ionIonisation::InitialiseEnergyLossProcess( 107 125 const G4ParticleDefinition* part,
Note: See TracChangeset
for help on using the changeset viewer.
