Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // >> 23 // $Id: G4IonFluctuations.cc,v 1.2 2004/12/01 19:37:14 vnivanch Exp $ >> 24 // GEANT4 tag $Name: geant4-07-00-patch-01 $ 26 // 25 // 27 // ------------------------------------------- 26 // ------------------------------------------------------------------- 28 // 27 // 29 // GEANT4 Class file 28 // GEANT4 Class file 30 // 29 // 31 // 30 // 32 // File name: G4IonFluctuation 31 // File name: G4IonFluctuation 33 // 32 // 34 // Author: Vladimir Ivanchenko 33 // Author: Vladimir Ivanchenko 35 // 34 // 36 // Creation date: 03.01.2002 35 // Creation date: 03.01.2002 37 // 36 // 38 // Modifications: 37 // Modifications: 39 // 38 // 40 // 28-12-02 add method Dispersion (V.Ivanchenk 39 // 28-12-02 add method Dispersion (V.Ivanchenko) 41 // 07-02-03 change signature (V.Ivanchenko) 40 // 07-02-03 change signature (V.Ivanchenko) 42 // 13-02-03 Add name (V.Ivanchenko) 41 // 13-02-03 Add name (V.Ivanchenko) 43 // 23-05-03 Add control on parthalogical cases 42 // 23-05-03 Add control on parthalogical cases (V.Ivanchenko) 44 // 16-10-03 Changed interface to Initialisatio 43 // 16-10-03 Changed interface to Initialisation (V.Ivanchenko) 45 // 27-09-07 Use FermiEnergy from material, add << 46 // 01-02-08 Add protection for small energies << 47 // 01-06-08 Added initialisation of effective << 48 // 44 // 49 // Class Description: 45 // Class Description: 50 // 46 // 51 // ------------------------------------------- 47 // ------------------------------------------------------------------- 52 // 48 // 53 49 54 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 55 //....oooOO0OOooo........oooOO0OOooo........oo 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 56 52 57 #include "G4IonFluctuations.hh" 53 #include "G4IonFluctuations.hh" 58 #include "G4UniversalFluctuation.hh" << 59 #include "G4PhysicalConstants.hh" << 60 #include "G4SystemOfUnits.hh" << 61 #include "Randomize.hh" 54 #include "Randomize.hh" 62 #include "G4Poisson.hh" 55 #include "G4Poisson.hh" 63 #include "G4MaterialCutsCouple.hh" << 64 #include "G4Material.hh" 56 #include "G4Material.hh" 65 #include "G4DynamicParticle.hh" 57 #include "G4DynamicParticle.hh" 66 #include "G4Pow.hh" << 58 #include "G4ParticleDefinition.hh" 67 #include "G4Log.hh" << 68 59 69 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 70 61 71 using namespace std; 62 using namespace std; 72 63 73 G4IonFluctuations::G4IonFluctuations(const G4S 64 G4IonFluctuations::G4IonFluctuations(const G4String& nam) 74 : G4VEmFluctuationModel(nam), << 65 :G4VEmFluctuationModel(nam), 75 particleMass(CLHEP::proton_mass_c2), << 66 particle(0), 76 parameter(10.0*CLHEP::MeV/CLHEP::proton_ma << 67 minNumberInteractionsBohr(10.0), 77 theBohrBeta2(50.0*keV/CLHEP::proton_mass_c << 68 theBohrBeta2(50.0*keV/proton_mass_c2), 78 minLoss(0.001*CLHEP::eV) << 69 minFraction(0.2), 79 { << 70 xmin(0.2), 80 uniFluct = new G4UniversalFluctuation(); << 71 minLoss(0.001*eV) 81 g4calc = G4Pow::GetInstance(); << 72 {} 82 } << 83 73 84 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 75 86 G4IonFluctuations::~G4IonFluctuations() = defa << 76 G4IonFluctuations::~G4IonFluctuations() >> 77 {} 87 78 88 //....oooOO0OOooo........oooOO0OOooo........oo 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 89 80 90 void G4IonFluctuations::InitialiseMe(const G4P 81 void G4IonFluctuations::InitialiseMe(const G4ParticleDefinition* part) 91 { 82 { 92 particle = part; << 83 particle = part; 93 particleMass = part->GetPDGMass(); << 84 particleMass = part->GetPDGMass(); 94 charge = part->GetPDGCharge()/eplus; << 85 charge = part->GetPDGCharge()/eplus; 95 effChargeSquare = chargeSquare = charge*char << 86 chargeSquare = charge*charge; 96 uniFluct->InitialiseMe(part); << 87 chargeSqRatio = 1.0; 97 } 88 } 98 89 99 //....oooOO0OOooo........oooOO0OOooo........oo 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 100 91 101 G4double << 92 G4double G4IonFluctuations::SampleFluctuations(const G4Material* material, 102 G4IonFluctuations::SampleFluctuations(const G4 << 93 const G4DynamicParticle* dp, 103 const G4 << 94 G4double& tmax, 104 const G4 << 95 G4double& length, 105 const G4 << 96 G4double& meanLoss) 106 const G4 << 107 const G4 << 108 { 97 { 109 // G4cout << "### meanLoss= " << meanLoss < << 110 if(meanLoss <= minLoss) return meanLoss; 98 if(meanLoss <= minLoss) return meanLoss; 111 99 112 //G4cout << "G4IonFluctuations::SampleFluctu << 100 if(dp->GetDefinition() != particle) { 113 // << dp->GetKineticEnergy() << 101 particle = dp->GetDefinition(); 114 // << " Elim(MeV)= " << parameter*c << 102 charge = particle->GetPDGCharge()/eplus; 115 << 116 // Vavilov fluctuations above energy thresho << 117 if(dp->GetKineticEnergy() > parameter*charge << 118 return uniFluct->SampleFluctuations(couple << 119 } 103 } 120 104 121 const G4Material* material = couple->GetMate << 105 G4double siga = Dispersion(material,dp,tmax,length); 122 G4double siga = Dispersion(material,dp,tcut, << 123 G4double loss = meanLoss; 106 G4double loss = meanLoss; 124 107 125 //G4cout << "### siga= " << sqrt(siga) << " << 126 108 127 // Gaussian fluctuation << 109 G4double navr = minNumberInteractionsBohr; 128 110 129 // Increase fluctuations for big fractional << 111 // Gaussian fluctuation 130 //G4cout << "siga= " << siga << G4endl; << 112 G4bool gauss = true; 131 if ( meanLoss > minFraction*kineticEnergy ) << 113 if (meanLoss >= minNumberInteractionsBohr*tmax) { 132 G4double gam = (kineticEnergy - meanLoss)/ << 114 navr = meanLoss*meanLoss/siga; 133 G4double b2 = 1.0 - 1.0/(gam*gam); << 115 if (navr < minNumberInteractionsBohr) gauss = false; 134 if(b2 < xmin*beta2) b2 = xmin*beta2; << 135 G4double x = b2/beta2; << 136 G4double x3 = 1.0/(x*x*x); << 137 siga *= 0.25*(1.0 + x)*(x3 + (1.0/b2 - 0.5 << 138 } 116 } 139 siga = std::sqrt(siga); << 140 G4double sn = meanLoss/siga; << 141 G4double twomeanLoss = meanLoss + meanLoss; << 142 // G4cout << "siga= " << siga << " sn= " < << 143 << 144 CLHEP::HepRandomEngine* rndmEngine = G4Rando << 145 // thick target case << 146 if (sn >= 2.0) { << 147 117 148 do { << 118 if(gauss) { 149 loss = G4RandGauss::shoot(rndmEngine,mea << 119 // Increase fluctuations for big fractional energy loss 150 // Loop checking, 03-Aug-2015, Vladimir << 151 } while (0.0 > loss || twomeanLoss < loss) << 152 120 153 // Gamma distribution << 121 //G4cout << "siga= " << siga << G4endl; 154 } else if(sn > 0.1) { << 122 if ( meanLoss > minFraction*kineticEnergy ) { >> 123 G4double gam = (kineticEnergy - meanLoss)/particleMass + 1.0; >> 124 G4double b2 = 1.0 - 1.0/(gam*gam); >> 125 if(b2 < xmin*beta2) b2 = xmin*beta2; >> 126 G4double x = b2/beta2; >> 127 G4double x3 = 1.0/(x*x*x); >> 128 siga *= 0.25*(1.0 + x)*(x3 + (1.0/b2 - 0.5)/(1.0/beta2 - 0.5) ); >> 129 } >> 130 // G4cout << "siga= " << siga << G4endl; >> 131 siga = sqrt(siga); 155 132 156 G4double neff = sn*sn; << 133 G4double lossmax = meanLoss+meanLoss; 157 loss = meanLoss*G4RandGamma::shoot(rndmEng << 134 do { >> 135 loss = G4RandGauss::shoot(meanLoss,siga); >> 136 } while (0.0 > loss || loss > lossmax); 158 137 159 // uniform distribution for very small ste << 138 // Poisson fluctuations 160 } else { 139 } else { 161 loss = twomeanLoss*rndmEngine->flat(); << 140 >> 141 G4double n = (G4double)(G4Poisson(navr)); >> 142 loss = meanLoss*n/navr; 162 } 143 } 163 144 164 //G4cout << "meanLoss= " << meanLoss << " lo << 145 // G4cout << "meanLoss= " << meanLoss << " loss= " << loss << G4endl; 165 return loss; 146 return loss; 166 } 147 } 167 148 168 //....oooOO0OOooo........oooOO0OOooo........oo 149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 169 150 170 G4double G4IonFluctuations::Dispersion(const G << 151 G4double G4IonFluctuations::Dispersion( 171 const G << 152 const G4Material* material, 172 const G << 153 const G4DynamicParticle* dp, 173 const G << 154 G4double& tmax, 174 const G << 155 G4double& length) 175 { 156 { 176 if(dp->GetDefinition() != particle) { Initia << 157 particleMass = dp->GetMass(); 177 << 158 G4double q = dp->GetCharge()/eplus; 178 const G4double beta = dp->GetBeta(); << 159 chargeSquare = q*q; 179 kineticEnergy = dp->GetKineticEnergy(); << 160 chargeSqRatio = chargeSquare/(charge*charge); 180 beta2 = beta*beta; << 161 181 << 162 G4double electronDensity = material->GetElectronDensity(); 182 G4double siga = (tmax/beta2 - 0.5*tcut)*CLHE << 163 kineticEnergy = dp->GetKineticEnergy(); 183 material->GetElectronDensity()*effChargeSq << 164 // G4cout << "e= " << kineticEnergy << " m= " << particleMass >> 165 // << " tmax= " << tmax << " l= " << length << " q^2= " << chargeSquare << G4endl; >> 166 G4double gam = kineticEnergy/particleMass + 1.0; >> 167 beta2 = 1.0 - 1.0/(gam*gam); >> 168 G4double siga = (1.0/beta2 - 0.5)*tmax*length*electronDensity*twopi_mc2_rcl2*chargeSquare; >> 169 // G4cout << "siga= " << siga << G4endl; 184 170 185 // Low velocity - additional ion charge fluc 171 // Low velocity - additional ion charge fluctuations according to 186 // Q.Yang et al., NIM B61(1991)149-155. 172 // Q.Yang et al., NIM B61(1991)149-155. 187 //G4cout << "sigE= " << sqrt(siga) << " char << 173 G4double zeff = electronDensity/(material->GetTotNbOfAtomsPerVolume()); >> 174 //G4cout << "siga= " << siga << " zeff= " << zeff << G4endl; 188 175 189 G4double Z = material->GetIonisation()->GetZ << 176 if ( beta2 < 3.0*theBohrBeta2*zeff ) { 190 G4double fac = Factor(material, Z); << 191 177 192 // heavy ion correction << 178 G4double a = CoeffitientA (zeff); 193 // G4double f1 = 1.065e-4*chargeSquare; << 179 G4double b = CoeffitientB (material, zeff); 194 // if(beta2 > theBohrBeta2) f1/= beta2; << 180 // G4cout << "a= " << a << " b= " << b << G4endl; 195 // else f1/= theBohrBe << 181 siga *= (a*chargeSqRatio + b); 196 // if(f1 > 2.5) f1 = 2.5; << 182 } else { 197 // fac *= (1.0 + f1); << 183 198 << 184 // H.Geissel et al. NIM B, 195 (2002) 3. 199 // taking into account the cutg << 185 siga *= RelativisticFactor(material, zeff); 200 G4double fac_cut = 1.0 + (fac - 1.0)*2.0*CLH << 201 /(tmax*(1.0 - beta2)); << 202 if(fac_cut > 0.01 && fac > 0.01) { << 203 siga *= fac_cut; << 204 } 186 } 205 /* << 187 206 G4cout << "siga(keV)= " << sqrt(siga)/keV << << 207 << " f1= " << fac_cut << G4endl; << 208 */ << 209 return siga; 188 return siga; 210 } 189 } 211 190 212 //....oooOO0OOooo........oooOO0OOooo........oo 191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 213 192 214 G4double G4IonFluctuations::Factor(const G4Mat << 193 G4double G4IonFluctuations::CoeffitientA(G4double& zeff) 215 { 194 { 216 // The aproximation of energy loss fluctuati 195 // The aproximation of energy loss fluctuations 217 // Q.Yang et al., NIM B61(1991)149-155. 196 // Q.Yang et al., NIM B61(1991)149-155. 218 197 219 // Reduced energy in MeV/AMU 198 // Reduced energy in MeV/AMU 220 G4double energy = kineticEnergy*CLHEP::amu_c << 199 G4double energy = kineticEnergy * amu_c2/(particleMass*MeV) ; 221 << 200 static G4double a[96][4] = { 222 // simple approximation for higher beta2 << 223 G4double s1 = RelativisticFactor(material, Z << 224 << 225 // tabulation for lower beta2 << 226 if( beta2 < 3.0*theBohrBeta2*Z ) { << 227 << 228 static const G4double a[96][4] = { << 229 {-0.3291, -0.8312, 0.2460, -1.0220}, 201 {-0.3291, -0.8312, 0.2460, -1.0220}, 230 {-0.5615, -0.5898, 0.5205, -0.7258}, 202 {-0.5615, -0.5898, 0.5205, -0.7258}, 231 {-0.5280, -0.4981, 0.5519, -0.5865}, 203 {-0.5280, -0.4981, 0.5519, -0.5865}, 232 {-0.5125, -0.4625, 0.5660, -0.5190}, 204 {-0.5125, -0.4625, 0.5660, -0.5190}, 233 {-0.5127, -0.8595, 0.5626, -0.8721}, 205 {-0.5127, -0.8595, 0.5626, -0.8721}, 234 {-0.5174, -1.1930, 0.5565, -1.1980}, 206 {-0.5174, -1.1930, 0.5565, -1.1980}, 235 {-0.5179, -1.1850, 0.5560, -1.2070}, 207 {-0.5179, -1.1850, 0.5560, -1.2070}, 236 {-0.5209, -0.9355, 0.5590, -1.0250}, 208 {-0.5209, -0.9355, 0.5590, -1.0250}, 237 {-0.5255, -0.7766, 0.5720, -0.9412}, 209 {-0.5255, -0.7766, 0.5720, -0.9412}, 238 210 239 {-0.5776, -0.6665, 0.6598, -0.8484}, 211 {-0.5776, -0.6665, 0.6598, -0.8484}, 240 {-0.6013, -0.6045, 0.7321, -0.7671}, 212 {-0.6013, -0.6045, 0.7321, -0.7671}, 241 {-0.5781, -0.5518, 0.7605, -0.6919}, 213 {-0.5781, -0.5518, 0.7605, -0.6919}, 242 {-0.5587, -0.4981, 0.7835, -0.6195}, 214 {-0.5587, -0.4981, 0.7835, -0.6195}, 243 {-0.5466, -0.4656, 0.7978, -0.5771}, 215 {-0.5466, -0.4656, 0.7978, -0.5771}, 244 {-0.5406, -0.4690, 0.8031, -0.5718}, 216 {-0.5406, -0.4690, 0.8031, -0.5718}, 245 {-0.5391, -0.5061, 0.8024, -0.5974}, 217 {-0.5391, -0.5061, 0.8024, -0.5974}, 246 {-0.5380, -0.6483, 0.7962, -0.6970}, 218 {-0.5380, -0.6483, 0.7962, -0.6970}, 247 {-0.5355, -0.7722, 0.7962, -0.7839}, 219 {-0.5355, -0.7722, 0.7962, -0.7839}, 248 {-0.5329, -0.7720, 0.7988, -0.7846}, 220 {-0.5329, -0.7720, 0.7988, -0.7846}, 249 221 250 {-0.5335, -0.7671, 0.7984, -0.7933}, 222 {-0.5335, -0.7671, 0.7984, -0.7933}, 251 {-0.5324, -0.7612, 0.7998, -0.8031}, 223 {-0.5324, -0.7612, 0.7998, -0.8031}, 252 {-0.5305, -0.7300, 0.8031, -0.7990}, 224 {-0.5305, -0.7300, 0.8031, -0.7990}, 253 {-0.5307, -0.7178, 0.8049, -0.8216}, 225 {-0.5307, -0.7178, 0.8049, -0.8216}, 254 {-0.5248, -0.6621, 0.8165, -0.7919}, 226 {-0.5248, -0.6621, 0.8165, -0.7919}, 255 {-0.5180, -0.6502, 0.8266, -0.7986}, 227 {-0.5180, -0.6502, 0.8266, -0.7986}, 256 {-0.5084, -0.6408, 0.8396, -0.8048}, 228 {-0.5084, -0.6408, 0.8396, -0.8048}, 257 {-0.4967, -0.6331, 0.8549, -0.8093}, 229 {-0.4967, -0.6331, 0.8549, -0.8093}, 258 {-0.4861, -0.6508, 0.8712, -0.8432}, 230 {-0.4861, -0.6508, 0.8712, -0.8432}, 259 {-0.4700, -0.6186, 0.8961, -0.8132}, 231 {-0.4700, -0.6186, 0.8961, -0.8132}, 260 232 261 {-0.4545, -0.5720, 0.9227, -0.7710}, 233 {-0.4545, -0.5720, 0.9227, -0.7710}, 262 {-0.4404, -0.5226, 0.9481, -0.7254}, 234 {-0.4404, -0.5226, 0.9481, -0.7254}, 263 {-0.4288, -0.4778, 0.9701, -0.6850}, 235 {-0.4288, -0.4778, 0.9701, -0.6850}, 264 {-0.4199, -0.4425, 0.9874, -0.6539}, 236 {-0.4199, -0.4425, 0.9874, -0.6539}, 265 {-0.4131, -0.4188, 0.9998, -0.6332}, 237 {-0.4131, -0.4188, 0.9998, -0.6332}, 266 {-0.4089, -0.4057, 1.0070, -0.6218}, 238 {-0.4089, -0.4057, 1.0070, -0.6218}, 267 {-0.4039, -0.3913, 1.0150, -0.6107}, 239 {-0.4039, -0.3913, 1.0150, -0.6107}, 268 {-0.3987, -0.3698, 1.0240, -0.5938}, 240 {-0.3987, -0.3698, 1.0240, -0.5938}, 269 {-0.3977, -0.3608, 1.0260, -0.5852}, 241 {-0.3977, -0.3608, 1.0260, -0.5852}, 270 {-0.3972, -0.3600, 1.0260, -0.5842}, 242 {-0.3972, -0.3600, 1.0260, -0.5842}, 271 243 272 {-0.3985, -0.3803, 1.0200, -0.6013}, << 244 {-0.3985, -0.3803, 1.0200, -0.6013}, 273 {-0.3985, -0.3979, 1.0150, -0.6168}, 245 {-0.3985, -0.3979, 1.0150, -0.6168}, 274 {-0.3968, -0.3990, 1.0160, -0.6195}, 246 {-0.3968, -0.3990, 1.0160, -0.6195}, 275 {-0.3971, -0.4432, 1.0050, -0.6591}, 247 {-0.3971, -0.4432, 1.0050, -0.6591}, 276 {-0.3944, -0.4665, 1.0010, -0.6825}, 248 {-0.3944, -0.4665, 1.0010, -0.6825}, 277 {-0.3924, -0.5109, 0.9921, -0.7235}, 249 {-0.3924, -0.5109, 0.9921, -0.7235}, 278 {-0.3882, -0.5158, 0.9947, -0.7343}, 250 {-0.3882, -0.5158, 0.9947, -0.7343}, 279 {-0.3838, -0.5125, 0.9999, -0.7370}, 251 {-0.3838, -0.5125, 0.9999, -0.7370}, 280 {-0.3786, -0.4976, 1.0090, -0.7310}, 252 {-0.3786, -0.4976, 1.0090, -0.7310}, 281 {-0.3741, -0.4738, 1.0200, -0.7155}, 253 {-0.3741, -0.4738, 1.0200, -0.7155}, 282 254 283 {-0.3969, -0.4496, 1.0320, -0.6982}, 255 {-0.3969, -0.4496, 1.0320, -0.6982}, 284 {-0.3663, -0.4297, 1.0430, -0.6828}, 256 {-0.3663, -0.4297, 1.0430, -0.6828}, 285 {-0.3630, -0.4120, 1.0530, -0.6689}, 257 {-0.3630, -0.4120, 1.0530, -0.6689}, 286 {-0.3597, -0.3964, 1.0620, -0.6564}, 258 {-0.3597, -0.3964, 1.0620, -0.6564}, 287 {-0.3555, -0.3809, 1.0720, -0.6454}, 259 {-0.3555, -0.3809, 1.0720, -0.6454}, 288 {-0.3525, -0.3607, 1.0820, -0.6289}, 260 {-0.3525, -0.3607, 1.0820, -0.6289}, 289 {-0.3505, -0.3465, 1.0900, -0.6171}, 261 {-0.3505, -0.3465, 1.0900, -0.6171}, 290 {-0.3397, -0.3570, 1.1020, -0.6384}, 262 {-0.3397, -0.3570, 1.1020, -0.6384}, 291 {-0.3314, -0.3552, 1.1130, -0.6441}, 263 {-0.3314, -0.3552, 1.1130, -0.6441}, 292 {-0.3235, -0.3531, 1.1230, -0.6498}, 264 {-0.3235, -0.3531, 1.1230, -0.6498}, 293 265 294 {-0.3150, -0.3483, 1.1360, -0.6539}, 266 {-0.3150, -0.3483, 1.1360, -0.6539}, 295 {-0.3060, -0.3441, 1.1490, -0.6593}, 267 {-0.3060, -0.3441, 1.1490, -0.6593}, 296 {-0.2968, -0.3396, 1.1630, -0.6649}, 268 {-0.2968, -0.3396, 1.1630, -0.6649}, 297 {-0.2935, -0.3225, 1.1760, -0.6527}, 269 {-0.2935, -0.3225, 1.1760, -0.6527}, 298 {-0.2797, -0.3262, 1.1940, -0.6722}, 270 {-0.2797, -0.3262, 1.1940, -0.6722}, 299 {-0.2704, -0.3202, 1.2100, -0.6770}, 271 {-0.2704, -0.3202, 1.2100, -0.6770}, 300 {-0.2815, -0.3227, 1.2480, -0.6775}, 272 {-0.2815, -0.3227, 1.2480, -0.6775}, 301 {-0.2880, -0.3245, 1.2810, -0.6801}, 273 {-0.2880, -0.3245, 1.2810, -0.6801}, 302 {-0.3034, -0.3263, 1.3270, -0.6778}, 274 {-0.3034, -0.3263, 1.3270, -0.6778}, 303 {-0.2936, -0.3215, 1.3430, -0.6835}, 275 {-0.2936, -0.3215, 1.3430, -0.6835}, 304 276 305 {-0.3282, -0.3200, 1.3980, -0.6650}, 277 {-0.3282, -0.3200, 1.3980, -0.6650}, 306 {-0.3260, -0.3070, 1.4090, -0.6552}, 278 {-0.3260, -0.3070, 1.4090, -0.6552}, 307 {-0.3511, -0.3074, 1.4470, -0.6442}, 279 {-0.3511, -0.3074, 1.4470, -0.6442}, 308 {-0.3501, -0.3064, 1.4500, -0.6442}, 280 {-0.3501, -0.3064, 1.4500, -0.6442}, 309 {-0.3490, -0.3027, 1.4550, -0.6418}, 281 {-0.3490, -0.3027, 1.4550, -0.6418}, 310 {-0.3487, -0.3048, 1.4570, -0.6447}, 282 {-0.3487, -0.3048, 1.4570, -0.6447}, 311 {-0.3478, -0.3074, 1.4600, -0.6483}, 283 {-0.3478, -0.3074, 1.4600, -0.6483}, 312 {-0.3501, -0.3283, 1.4540, -0.6669}, 284 {-0.3501, -0.3283, 1.4540, -0.6669}, 313 {-0.3494, -0.3373, 1.4550, -0.6765}, 285 {-0.3494, -0.3373, 1.4550, -0.6765}, 314 {-0.3485, -0.3373, 1.4570, -0.6774}, 286 {-0.3485, -0.3373, 1.4570, -0.6774}, 315 287 316 {-0.3462, -0.3300, 1.4630, -0.6728}, 288 {-0.3462, -0.3300, 1.4630, -0.6728}, 317 {-0.3462, -0.3225, 1.4690, -0.6662}, 289 {-0.3462, -0.3225, 1.4690, -0.6662}, 318 {-0.3453, -0.3094, 1.4790, -0.6553}, 290 {-0.3453, -0.3094, 1.4790, -0.6553}, 319 {-0.3844, -0.3134, 1.5240, -0.6412}, 291 {-0.3844, -0.3134, 1.5240, -0.6412}, 320 {-0.3848, -0.3018, 1.5310, -0.6303}, 292 {-0.3848, -0.3018, 1.5310, -0.6303}, 321 {-0.3862, -0.2955, 1.5360, -0.6237}, 293 {-0.3862, -0.2955, 1.5360, -0.6237}, 322 {-0.4262, -0.2991, 1.5860, -0.6115}, 294 {-0.4262, -0.2991, 1.5860, -0.6115}, 323 {-0.4278, -0.2910, 1.5900, -0.6029}, 295 {-0.4278, -0.2910, 1.5900, -0.6029}, 324 {-0.4303, -0.2817, 1.5940, -0.5927}, 296 {-0.4303, -0.2817, 1.5940, -0.5927}, 325 {-0.4315, -0.2719, 1.6010, -0.5829}, 297 {-0.4315, -0.2719, 1.6010, -0.5829}, 326 298 327 {-0.4359, -0.2914, 1.6050, -0.6010}, 299 {-0.4359, -0.2914, 1.6050, -0.6010}, 328 {-0.4365, -0.2982, 1.6080, -0.6080}, 300 {-0.4365, -0.2982, 1.6080, -0.6080}, 329 {-0.4253, -0.3037, 1.6120, -0.6150}, 301 {-0.4253, -0.3037, 1.6120, -0.6150}, 330 {-0.4335, -0.3245, 1.6160, -0.6377}, 302 {-0.4335, -0.3245, 1.6160, -0.6377}, 331 {-0.4307, -0.3292, 1.6210, -0.6447}, 303 {-0.4307, -0.3292, 1.6210, -0.6447}, 332 {-0.4284, -0.3204, 1.6290, -0.6380}, 304 {-0.4284, -0.3204, 1.6290, -0.6380}, 333 {-0.4227, -0.3217, 1.6360, -0.6438} 305 {-0.4227, -0.3217, 1.6360, -0.6438} 334 } ; << 306 } ; 335 307 336 G4int iz = G4lrint(Z) - 2; << 308 G4int iz = (G4int)zeff - 2 ; 337 if( 0 > iz ) { iz = 0; } << 309 if( 0 > iz ) iz = 0 ; 338 else if(95 < iz ) { iz = 95; } << 310 if(95 < iz ) iz = 95 ; >> 311 >> 312 G4double q = 1.0 / (1.0 + a[iz][0]*pow(energy,a[iz][1])+ >> 313 + a[iz][2]*pow(energy,a[iz][3])) ; >> 314 >> 315 return q ; >> 316 } >> 317 >> 318 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 319 >> 320 G4double G4IonFluctuations::CoeffitientB(const G4Material* material, G4double& zeff) >> 321 { >> 322 // The aproximation of energy loss fluctuations >> 323 // Q.Yang et al., NIM B61(1991)149-155. >> 324 >> 325 // Reduced energy in MeV/AMU >> 326 G4double energy = kineticEnergy *amu_c2/(particleMass*MeV) ; 339 327 340 const G4double ss = 1.0 + a[iz][0]*g4calc- << 341 + a[iz][2]*g4calc->powA(energy,a[iz][3]) << 342 << 343 // protection for the validity range for l << 344 static const G4double slim = 0.001; << 345 if(ss < slim) { s1 = 1.0/slim; } << 346 // for high value of beta << 347 else if(s1*ss < 1.0) { s1 = 1.0/ss; } << 348 } << 349 G4int i = 0 ; 328 G4int i = 0 ; 350 G4double factor = 1.0 ; 329 G4double factor = 1.0 ; 351 330 352 // The index of set of parameters i = 0 for 331 // The index of set of parameters i = 0 for protons(hadrons) in gases 353 // 1 for 332 // 1 for protons(hadrons) in solids 354 // 2 for 333 // 2 for ions in atomic gases 355 // 3 for 334 // 3 for ions in molecular gases 356 // 4 for 335 // 4 for ions in solids 357 static const G4double b[5][4] = { << 336 static G4double b[5][4] = { 358 {0.1014, 0.3700, 0.9642, 3.987}, 337 {0.1014, 0.3700, 0.9642, 3.987}, 359 {0.1955, 0.6941, 2.522, 1.040}, 338 {0.1955, 0.6941, 2.522, 1.040}, 360 {0.05058, 0.08975, 0.1419, 10.80}, 339 {0.05058, 0.08975, 0.1419, 10.80}, 361 {0.05009, 0.08660, 0.2751, 3.787}, 340 {0.05009, 0.08660, 0.2751, 3.787}, 362 {0.01273, 0.03458, 0.3951, 3.812} 341 {0.01273, 0.03458, 0.3951, 3.812} 363 } ; 342 } ; 364 343 365 // protons (hadrons) 344 // protons (hadrons) 366 if(1.5 > charge) { 345 if(1.5 > charge) { 367 if( kStateGas != material->GetState() ) { << 346 if( kStateGas != material->GetState() ) i = 1 ; 368 347 369 // ions 348 // ions 370 } else { 349 } else { 371 << 350 factor = charge * pow(charge/zeff, 0.3333) ; 372 factor = charge * g4calc->A13(charge/Z); << 373 351 374 if( kStateGas == material->GetState() ) { 352 if( kStateGas == material->GetState() ) { 375 energy /= (charge * std::sqrt(charge)) ; << 353 energy /= (charge * sqrt(charge)) ; 376 354 377 if(1 == (material->GetNumberOfElements() 355 if(1 == (material->GetNumberOfElements())) { 378 i = 2 ; 356 i = 2 ; 379 } else { 357 } else { 380 i = 3 ; 358 i = 3 ; 381 } 359 } 382 360 383 } else { 361 } else { 384 energy /= (charge * std::sqrt(charge*Z)) << 362 energy /= (charge * sqrt(charge*zeff)) ; 385 i = 4 ; 363 i = 4 ; 386 } 364 } 387 } 365 } 388 366 389 G4double x = b[i][2]; << 367 G4double x = b[i][2] * (1.0 - exp( - energy * b[i][3] )) ; 390 G4double y = energy * b[i][3]; << 391 if(y <= 0.2) x *= (y*(1.0 - 0.5*y)); << 392 else x *= (1.0 - g4calc->expA(-y)); << 393 << 394 y = energy - b[i][1]; << 395 << 396 const G4double s2 = factor * x * b[i][0] / ( << 397 /* << 398 G4cout << "s1= " << s1 << " s2= " << s2 << " << 399 << " e= " << energy << G4endl; << 400 */ << 401 return s1*effChargeSquare/chargeSquare + s2; << 402 } << 403 << 404 //....oooOO0OOooo........oooOO0OOooo........oo << 405 << 406 G4double G4IonFluctuations::RelativisticFactor << 407 << 408 { << 409 G4double eF = mat->GetIonisation()->GetFermi << 410 G4double I = mat->GetIonisation()->GetMeanE << 411 << 412 // H.Geissel et al. NIM B, 195 (2002) 3. << 413 G4double bF2= 2.0*eF/CLHEP::electron_mass_c2 << 414 G4double f = 0.4*(1.0 - beta2)/((1.0 - 0.5* << 415 if(beta2 > bF2) f *= G4Log(2.0*CLHEP::electr << 416 else f *= G4Log(4.0*eF/I); << 417 368 418 // G4cout << "f= " << f << " beta2= " << be << 369 G4double q = factor * x * b[i][0] / 419 // << " bf2= " << bF2 << " q^2= " << << 370 ((energy - b[i][1])*(energy - b[i][1]) + x*x) ; 420 371 421 return 1.0 + f; << 372 return q ; 422 } 373 } 423 374 424 //....oooOO0OOooo........oooOO0OOooo........oo 375 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 425 376 426 void G4IonFluctuations::SetParticleAndCharge(c << 377 G4double G4IonFluctuations::RelativisticFactor(const G4Material*, G4double& zeff) 427 G << 428 { 378 { 429 if(part != particle) { << 379 // H.Geissel et al. NIM B, 195 (2002) 3. 430 particle = part; << 380 G4double factor = 1.0 + 0.667*theBohrBeta2*(1.0 - beta2) 431 particleMass = part->GetPDGMass(); << 381 * log(2.0*electron_mass_c2/(5.0*charge*eV)) 432 charge = part->GetPDGCharge()/eplu << 382 / ((1.0 - 0.5*beta2)*beta2*zeff) ; 433 chargeSquare = charge*charge; << 383 factor *= (1.0 + 1.415e-4*chargeSquare/beta2); 434 } << 384 // G4cout << "factor= " << factor << G4endl; 435 effChargeSquare = q2; << 385 return factor; 436 uniFluct->SetParticleAndCharge(part, q2); << 437 } 386 } 438 387 439 //....oooOO0OOooo........oooOO0OOooo........oo 388 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 440 389