Geant4 Cross Reference |
1 // << 1 // 2 // ******************************************* << 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * License and Disclaimer * 4 // * << 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided << 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License << 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ << 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. << 9 // * include a list of copyright holders. * 10 // * << 10 // * * 11 // * Neither the authors of this software syst << 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin << 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran << 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum << 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio << 16 // * for the full disclaimer and the limitation of liability. * 17 // * << 17 // * * 18 // * This code implementation is the result << 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio << 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag << 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati << 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof << 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* << 24 // ******************************************************************** 25 // << 25 // 26 // << 26 // $Id: G4EmCorrections.cc,v 1.51.2.1 2010/01/26 14:33:54 gcosmo Exp $ 27 // ------------------------------------------- << 27 // GEANT4 tag $Name: geant4-09-02-patch-04 $ 28 // << 28 // 29 // GEANT4 Class << 29 // ------------------------------------------------------------------- 30 // << 30 // 31 // File name: G4EmCorrections << 31 // GEANT4 Class 32 // << 32 // 33 // Author: Vladimir Ivanchenko << 33 // File name: G4EmCorrections 34 // << 34 // 35 // Creation date: 13.01.2005 << 35 // Author: Vladimir Ivanchenko 36 // << 36 // 37 // Modifications: << 37 // Creation date: 13.01.2005 38 // 05.05.2005 V.Ivanchenko Fix misprint in Mot << 38 // 39 // 26.11.2005 V.Ivanchenko Fix effective charg << 39 // Modifications: 40 // 28.04.2006 V.Ivanchenko General cleanup, ad << 40 // 05.05.2005 V.Ivanchenko Fix misprint in Mott term 41 // 13.05.2006 V.Ivanchenko Add corrections for << 41 // 26.11.2005 V.Ivanchenko Fix effective charge for heavy ions using original paper 42 // 08.05.2007 V.Ivanchenko Use G4IonTable for << 42 // 28.04.2006 V.Ivanchenko General cleanup, add finite size corrections 43 // division by zero << 43 // 13.05.2006 V.Ivanchenko Add corrections for ion stopping 44 // 29.02.2008 V.Ivanchenko use expantions for << 44 // 08.05.2007 V.Ivanchenko Use G4IonTable for ion mass instead of NistTable to avoid 45 // 21.04.2008 Updated computations for ions (V << 45 // division by zero 46 // 20.05.2008 Removed Finite Size correction ( << 46 // 29.02.2008 V.Ivanchenko use expantions for log and power function 47 // 19.04.2012 Fix reproducibility problem (A.R << 47 // 21.04.2008 Updated computations for ions (V.Ivanchenko) 48 // << 48 // 20.05.2008 Removed Finite Size correction (V.Ivanchenko) 49 // << 49 // 50 // Class Description: << 50 // 51 // << 51 // Class Description: 52 // This class provides calculation of EM corre << 52 // 53 // << 53 // This class provides calculation of EM corrections to ionisation 54 << 54 // 55 // ------------------------------------------- << 55 56 // << 56 // ------------------------------------------------------------------- 57 << 57 // 58 #include "G4EmCorrections.hh" << 58 59 #include "G4PhysicalConstants.hh" << 59 #include "G4EmCorrections.hh" 60 #include "G4SystemOfUnits.hh" << 60 #include "Randomize.hh" 61 #include "G4ParticleTable.hh" << 61 #include "G4ParticleTable.hh" 62 #include "G4IonTable.hh" << 62 #include "G4IonTable.hh" 63 #include "G4VEmModel.hh" << 63 #include "G4VEmModel.hh" 64 #include "G4Proton.hh" << 64 #include "G4Proton.hh" 65 #include "G4GenericIon.hh" << 65 #include "G4GenericIon.hh" 66 #include "G4PhysicsLogVector.hh" << 66 #include "G4LPhysicsFreeVector.hh" 67 #include "G4ProductionCutsTable.hh" << 67 #include "G4PhysicsLogVector.hh" 68 #include "G4MaterialCutsCouple.hh" << 68 #include "G4ProductionCutsTable.hh" 69 #include "G4AtomicShells.hh" << 69 #include "G4MaterialCutsCouple.hh" 70 #include "G4Log.hh" << 70 71 #include "G4Exp.hh" << 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 72 #include "G4Pow.hh" << 72 73 << 73 G4EmCorrections::G4EmCorrections() 74 //....oooOO0OOooo........oooOO0OOooo........oo << 74 { 75 << 75 particle = 0; 76 namespace << 76 curParticle= 0; 77 { << 77 material = 0; 78 constexpr G4double inveplus = 1.0/CLHEP::epl << 78 curMaterial= 0; 79 constexpr G4double alpha2 = CLHEP::fine_stru << 79 curVector = 0; 80 } << 80 kinEnergy = 0.0; 81 const G4double G4EmCorrections::ZD[11] = << 81 ionLEModel = 0; 82 {0., 0., 0., 1.72, 2.09, 2.48, 2.82, 3.16, << 82 ionHEModel = 0; 83 const G4double G4EmCorrections::UK[20] = {1.99 << 83 nIons = 0; 84 2.0817, 2.0945, 2.0 << 84 verbose = 1; 85 2.1197, 2.1246, 2.1 << 85 ncouples = 0; 86 2.1310, 2.1310, 2.1 << 86 massFactor = 1.0; 87 const G4double G4EmCorrections::VK[20] = {8.34 << 87 eth = 2.0*MeV; 88 8.3219, 8.3201, 8.3 << 88 nbinCorr = 20; 89 8.3191, 8.3199, 8.3 << 89 eCorrMin = 25.*keV; 90 8.3244, 8.3264, 8.3 << 90 eCorrMax = 250.*MeV; 91 G4double G4EmCorrections::ZK[] = {0.0}; << 91 nist = G4NistManager::Instance(); 92 const G4double G4EmCorrections::Eta[29] = {0.0 << 92 ionTable = G4ParticleTable::GetParticleTable()->GetIonTable(); 93 0.03, 0.04, 0.05 << 93 Initialise(); 94 0.1, 0.15, 0.2, << 94 } 95 0.5, 0.6, 0.7, << 95 96 1.2, 1.4, 1.5, << 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 97 G4double G4EmCorrections::CK[20][29]; << 97 98 G4double G4EmCorrections::CL[26][28]; << 98 G4EmCorrections::~G4EmCorrections() 99 const G4double G4EmCorrections::UL[] = {0.1215 << 99 { 100 1.4379, 1.5032, 1.5 << 100 for(G4int i=0; i<nIons; i++) {delete stopData[i];} 101 1.8036, 1.8543, 1.8 << 101 } 102 1.9508, 1.9696, 1.9 << 102 103 2.0001, 2.0039, 2.0 << 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 104 G4double G4EmCorrections::VL[] = {0.0}; << 104 105 G4double G4EmCorrections::sWmaxBarkas = 10.0; << 105 G4double G4EmCorrections::HighOrderCorrections(const G4ParticleDefinition* p, 106 << 106 const G4Material* mat, 107 G4PhysicsFreeVector* G4EmCorrections::sBarkasC << 107 G4double e, G4double) 108 G4PhysicsFreeVector* G4EmCorrections::sThetaK << 108 { 109 G4PhysicsFreeVector* G4EmCorrections::sThetaL << 109 // . Z^3 Barkas effect in the stopping power of matter for charged particles 110 << 110 // J.C Ashley and R.H.Ritchie 111 G4EmCorrections::G4EmCorrections(G4int verb) << 111 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397 112 : verbose(verb) << 112 // and ICRU49 report 113 { << 113 // valid for kineticEnergy < 0.5 MeV 114 eth = 2.0*CLHEP::MeV; << 114 // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980 115 eCorrMin = 25.*CLHEP::keV; << 115 116 eCorrMax = 1.*CLHEP::GeV; << 116 SetupKinematics(p, mat, e); 117 << 117 if(tau <= 0.0) return 0.0; 118 ionTable = G4ParticleTable::GetParticleTable << 118 119 g4calc = G4Pow::GetInstance(); << 119 G4double Barkas = BarkasCorrection (p, mat, e); 120 << 120 G4double Bloch = BlochCorrection (p, mat, e); 121 // fill vectors << 121 G4double Mott = MottCorrection (p, mat, e); 122 if (nullptr == sBarkasCorr) { << 122 123 Initialise(); << 123 G4double sum = (2.0*(Barkas + Bloch) + Mott); 124 isInitializer = true; << 124 125 } << 125 if(verbose > 1) 126 } << 126 G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas 127 << 127 << " Bloch= " << Bloch << " Mott= " << Mott 128 //....oooOO0OOooo........oooOO0OOooo........oo << 128 << " Sum= " << sum << G4endl; 129 << 129 130 G4EmCorrections::~G4EmCorrections() << 130 sum *= material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2; 131 { << 131 return sum; 132 for (G4int i=0; i<nIons; ++i) { delete stopD << 132 } 133 if (isInitializer) { << 133 134 delete sBarkasCorr; << 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 135 delete sThetaK; << 135 136 delete sThetaL; << 136 G4double G4EmCorrections::IonBarkasCorrection(const G4ParticleDefinition* p, 137 sBarkasCorr = sThetaK = sThetaL = nullptr; << 137 const G4Material* mat, 138 } << 138 G4double e) 139 } << 139 { 140 << 140 // . Z^3 Barkas effect in the stopping power of matter for charged particles 141 void G4EmCorrections::SetupKinematics(const G4 << 141 // J.C Ashley and R.H.Ritchie 142 const G4Material* mat, << 142 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397 143 const G4double kineticEnergy) << 143 // and ICRU49 report 144 { << 144 // valid for kineticEnergy < 0.5 MeV 145 if(kineticEnergy != kinEnergy || p != partic << 145 146 particle = p; << 146 SetupKinematics(p, mat, e); 147 kinEnergy = kineticEnergy; << 147 G4double res = 0.0; 148 mass = p->GetPDGMass(); << 148 if(tau > 0.0) 149 tau = kineticEnergy / mass; << 149 res = 2.0*BarkasCorrection(p, mat, e)* 150 gamma = 1.0 + tau; << 150 material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2; 151 bg2 = tau * (tau+2.0); << 151 return res; 152 beta2 = bg2/(gamma*gamma); << 152 } 153 beta = std::sqrt(beta2); << 153 154 ba2 = beta2/alpha2; << 154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 155 const G4double ratio = CLHEP::electron_mas << 155 156 tmax = 2.0*CLHEP::electron_mass_c2*bg2 << 156 G4double G4EmCorrections::ComputeIonCorrections(const G4ParticleDefinition* p, 157 /(1. + 2.0*gamma*ratio + ratio*ratio); << 157 const G4Material* mat, 158 charge = p->GetPDGCharge()*inveplus; << 158 G4double e) 159 if(charge > 1.5) { charge = effCharge.Effe << 159 { 160 q2 = charge*charge; << 160 // . Z^3 Barkas effect in the stopping power of matter for charged particles 161 } << 161 // J.C Ashley and R.H.Ritchie 162 if(mat != material) { << 162 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397 163 material = mat; << 163 // and ICRU49 report 164 theElementVector = material->GetElementVec << 164 // valid for kineticEnergy < 0.5 MeV 165 atomDensity = material->GetAtomicNumDensi << 165 // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980 166 numberOfElements = (G4int)material->GetNum << 166 SetupKinematics(p, mat, e); 167 } << 167 if(tau <= 0.0) return 0.0; 168 } << 168 169 << 169 G4double Barkas = BarkasCorrection (p, mat, e); 170 //....oooOO0OOooo........oooOO0OOooo........oo << 170 G4double Bloch = BlochCorrection (p, mat, e); 171 << 171 G4double Mott = MottCorrection (p, mat, e); 172 G4double G4EmCorrections::HighOrderCorrections << 172 173 << 173 G4double sum = 2.0*(Barkas*(charge - 1.0)/charge + Bloch) + Mott; 174 << 174 175 { << 175 if(verbose > 1) { 176 // . Z^3 Barkas effect in the stopping power << 176 G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas 177 // J.C Ashley and R.H.Ritchie << 177 << " Bloch= " << Bloch << " Mott= " << Mott 178 // Physical review B Vol.5 No.7 1 April 19 << 178 << " Sum= " << sum << G4endl; 179 // and ICRU49 report << 179 } 180 // valid for kineticEnergy < 0.5 MeV << 180 sum *= material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2; 181 // Other corrections from S.P.Ahlen Rev. M << 181 182 << 182 if(verbose > 1) G4cout << " Sum= " << sum << G4endl; 183 SetupKinematics(p, mat, e); << 183 return sum; 184 if(tau <= 0.0) { return 0.0; } << 184 } 185 << 185 186 const G4double Barkas = BarkasCorrection(p, << 186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 187 const G4double Bloch = BlochCorrection(p, m << 187 188 const G4double Mott = MottCorrection(p, mat, << 188 G4double G4EmCorrections::IonHighOrderCorrections(const G4ParticleDefinition* p, 189 << 189 const G4MaterialCutsCouple* couple, 190 G4double sum = 2.0*(Barkas + Bloch) + Mott; << 190 G4double e) 191 << 191 { 192 if(verbose > 1) { << 192 // . Z^3 Barkas effect in the stopping power of matter for charged particles 193 G4cout << "EmCorrections: E(MeV)= " << e/M << 193 // J.C Ashley and R.H.Ritchie 194 << " Bloch= " << Bloch << " Mott= " << 194 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397 195 << " Sum= " << sum << " q2= " << q2 << 195 // and ICRU49 report 196 G4cout << " ShellCorrection: " << ShellCor << 196 // valid for kineticEnergy < 0.5 MeV 197 << " Kshell= " << KShellCorrection( << 197 // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980 198 << " Lshell= " << LShellCorrection( << 198 199 << " " << mat->GetName() << G4end << 199 G4double sum = 0.0; 200 } << 200 201 sum *= material->GetElectronDensity()*q2*CLH << 201 if(ionHEModel) { 202 return sum; << 202 G4int Z = G4int(p->GetPDGCharge()/eplus + 0.5); 203 } << 203 if(Z >= 100) Z = 99; 204 << 204 else if(Z < 1) Z = 1; 205 //....oooOO0OOooo........oooOO0OOooo........oo << 205 206 << 206 // fill vector 207 G4double G4EmCorrections::IonBarkasCorrection( << 207 if(thcorr[Z].size() == 0) { 208 << 208 thcorr[Z].resize(ncouples); 209 << 209 G4double ethscaled = eth*p->GetPDGMass()/proton_mass_c2; 210 { << 210 211 SetupKinematics(p, mat, e); << 211 for(size_t i=0; i<ncouples; i++) { 212 return 2.0*BarkasCorrection(p, mat, e, true) << 212 (thcorr[Z])[i] = ethscaled*ComputeIonCorrections(p, currmat[i], ethscaled); 213 material->GetElectronDensity() * q2 * CLHE << 213 //G4cout << i << ". ethscaled= " << ethscaled 214 } << 214 //<< " corr= " << (thcorr[Z])[i]/ethscaled << G4endl; 215 << 215 } 216 //....oooOO0OOooo........oooOO0OOooo........oo << 216 } 217 << 217 G4double rest = (thcorr[Z])[couple->GetIndex()]; 218 G4double G4EmCorrections::ComputeIonCorrection << 218 219 << 219 sum = ComputeIonCorrections(p,couple->GetMaterial(),e) - rest/e; 220 << 220 221 { << 221 if(verbose > 1) G4cout << " Sum= " << sum << " dSum= " << rest/e << G4endl; 222 // . Z^3 Barkas effect in the stopping power << 222 } 223 // J.C Ashley and R.H.Ritchie << 223 return sum; 224 // Physical review B Vol.5 No.7 1 April 19 << 224 } 225 // and ICRU49 report << 225 226 // valid for kineticEnergy < 0.5 MeV << 226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 227 // Other corrections from S.P.Ahlen Rev. M << 227 228 SetupKinematics(p, mat, e); << 228 G4double G4EmCorrections::Bethe(const G4ParticleDefinition* p, 229 if(tau <= 0.0) { return 0.0; } << 229 const G4Material* mat, 230 << 230 G4double e) 231 const G4double Barkas = BarkasCorrection (p, << 231 { 232 const G4double Bloch = BlochCorrection (p, << 232 SetupKinematics(p, mat, e); 233 const G4double Mott = MottCorrection (p, mat << 233 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy(); 234 << 234 G4double eexc2 = eexc*eexc; 235 G4double sum = 2.0*(Barkas*(charge - 1.0)/ch << 235 G4double dedx = 0.5*std::log(2.0*electron_mass_c2*bg2*tmax/eexc2)-beta2; 236 << 236 return dedx; 237 if(verbose > 1) { << 237 } 238 G4cout << "EmCorrections: E(MeV)= " << e/M << 238 239 << " Bloch= " << Bloch << " Mott= " << 239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 240 << " Sum= " << sum << G4endl; << 240 241 } << 241 G4double G4EmCorrections::SpinCorrection(const G4ParticleDefinition* p, 242 sum *= material->GetElectronDensity() * q2 * << 242 const G4Material* mat, 243 << 243 G4double e) 244 if(verbose > 1) { G4cout << " Sum= " << sum << 244 { 245 return sum; << 245 SetupKinematics(p, mat, e); 246 } << 246 G4double dedx = 0.5*tmax/(kinEnergy + mass); 247 << 247 return 0.5*dedx*dedx; 248 //....oooOO0OOooo........oooOO0OOooo........oo << 248 } 249 << 249 250 G4double G4EmCorrections::IonHighOrderCorrecti << 250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 251 << 251 252 << 252 G4double G4EmCorrections:: KShellCorrection(const G4ParticleDefinition* p, 253 { << 253 const G4Material* mat, 254 // . Z^3 Barkas effect in the stopping power << 254 G4double e) 255 // J.C Ashley and R.H.Ritchie << 255 { 256 // Physical review B Vol.5 No.7 1 April 19 << 256 SetupKinematics(p, mat, e); 257 // and ICRU49 report << 257 G4double term = 0.0; 258 // valid for kineticEnergy < 0.5 MeV << 258 for (G4int i = 0; i<numberOfElements; i++) { 259 // Other corrections from S.P.Ahlen Rev. M << 259 260 << 260 G4double Z = (*theElementVector)[i]->GetZ(); 261 G4double sum = 0.0; << 261 G4int iz = G4int(Z); 262 << 262 G4double f = 1.0; 263 if (nullptr != ionHEModel) { << 263 G4double Z2= (Z-0.3)*(Z-0.3); 264 G4int Z = G4lrint(p->GetPDGCharge()*invepl << 264 if(1 == iz) { 265 Z = std::max(std::min(Z, 99), 1); << 265 f = 0.5; 266 << 266 Z2 = 1.0; 267 const G4double ethscaled = eth*p->GetPDGMa << 267 } 268 const G4int ionPDG = p->GetPDGEncoding(); << 268 G4double e0= 13.6*eV*Z2; 269 auto iter = thcorr.find(ionPDG); << 269 term += f*atomDensity[i]*KShell(shells.GetBindingEnergy(iz,0)/e0,ba2/Z2)/Z; 270 if (iter == thcorr.end()) { // Not found: << 270 } 271 std::vector<G4double> v; << 271 272 for(std::size_t i=0; i<ncouples; ++i){ << 272 term /= material->GetTotNbOfAtomsPerVolume(); 273 v.push_back(ethscaled*ComputeIonCorrec << 273 274 } << 274 return term; 275 thcorr.insert(std::pair< G4int, std::vec << 275 } 276 } << 276 277 G4double rest = 0.0; << 277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 278 iter = thcorr.find(ionPDG); << 278 279 if (iter != thcorr.end()) { rest = (iter-> << 279 G4double G4EmCorrections:: LShellCorrection(const G4ParticleDefinition* p, 280 << 280 const G4Material* mat, 281 sum = ComputeIonCorrections(p,couple->GetM << 281 G4double e) 282 << 282 { 283 if(verbose > 1) { << 283 SetupKinematics(p, mat, e); 284 G4cout << " Sum= " << sum << " dSum= " < << 284 G4double term = 0.0; 285 } << 285 for (G4int i = 0; i<numberOfElements; i++) { 286 } << 286 287 return sum; << 287 G4double Z = (*theElementVector)[i]->GetZ(); 288 } << 288 G4int iz = G4int(Z); 289 << 289 if(2 < iz) { 290 //....oooOO0OOooo........oooOO0OOooo........oo << 290 G4double Zeff = Z - ZD[10]; 291 << 291 if(iz < 10) Zeff = Z - ZD[iz]; 292 G4double G4EmCorrections::Bethe(const G4Partic << 292 G4double Z2= Zeff*Zeff; 293 const G4Materi << 293 G4double e0= 13.6*eV*Z2*0.25; 294 const G4double << 294 G4double f = 0.125; 295 { << 295 G4int nmax = std::min(4,shells.GetNumberOfShells(iz)); 296 SetupKinematics(p, mat, e); << 296 for(G4int j=1; j<nmax; j++) { 297 const G4double eexc = material->GetIonisati << 297 G4double ne = G4double(shells.GetNumberOfElectrons(iz,j)); 298 const G4double eexc2 = eexc*eexc; << 298 G4double e1 = shells.GetBindingEnergy(iz,j); 299 return 0.5*G4Log(2.0*electron_mass_c2*bg2*tm << 299 // G4cout << "LShell: j= " << j << " ne= " << ne << " e(eV)= " << e/eV 300 } << 300 // << " e0(eV)= " << e0/eV << G4endl; 301 << 301 term += f*ne*atomDensity[i]*LShell(e1/e0,ba2/Z2)/Z; 302 //....oooOO0OOooo........oooOO0OOooo........oo << 302 } 303 << 303 } 304 G4double G4EmCorrections::SpinCorrection(const << 304 } 305 const << 305 306 const << 306 term /= material->GetTotNbOfAtomsPerVolume(); 307 { << 307 308 SetupKinematics(p, mat, e); << 308 return term; 309 const G4double dedx = 0.5*tmax/(kinEnergy + << 309 } 310 return 0.5*dedx*dedx; << 310 311 } << 311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 312 << 312 313 //....oooOO0OOooo........oooOO0OOooo........oo << 313 G4double G4EmCorrections::KShell(G4double tet, G4double eta) 314 << 314 { 315 G4double G4EmCorrections:: KShellCorrection(co << 315 G4double corr = 0.0; 316 co << 316 317 co << 317 G4double x = tet; 318 { << 318 if(tet < TheK[0]) x = TheK[0]; 319 SetupKinematics(p, mat, e); << 319 G4int itet = Index(x, TheK, nK); 320 G4double term = 0.0; << 320 321 for (G4int i = 0; i<numberOfElements; ++i) { << 321 // assimptotic case 322 << 322 if(eta >= Eta[nEtaK-1]) { 323 G4double Z = (*theElementVector)[i]->GetZ( << 323 corr = (Value(x, TheK[itet], TheK[itet+1], UK[itet], UK[itet+1]) + 324 G4int iz = (*theElementVector)[i]->GetZa << 324 Value(x, TheK[itet], TheK[itet+1], VK[itet], VK[itet+1])/eta + 325 G4double f = 1.0; << 325 Value(x, TheK[itet], TheK[itet+1], ZK[itet], ZK[itet+1])/(eta*eta))/eta; 326 G4double Z2= (Z-0.3)*(Z-0.3); << 326 } else { 327 if(1 == iz) { << 327 G4double y = eta; 328 f = 0.5; << 328 if(eta < Eta[0]) y = Eta[0]; 329 Z2 = 1.0; << 329 G4int ieta = Index(y, Eta, nEtaK); 330 } << 330 corr = Value2(x, y, TheK[itet], TheK[itet+1], Eta[ieta], Eta[ieta+1], 331 const G4double eta = ba2/Z2; << 331 CK[itet][ieta], CK[itet+1][ieta], 332 const G4double tet = (11 < iz) ? sThetaK-> << 332 CK[itet][ieta+1], CK[itet+1][ieta+1]); 333 term += f*atomDensity[i]*KShell(tet,eta)/Z << 333 //G4cout << " x= " <<x<<" y= "<<y<<" tet= " <<TheK[itet] 334 } << 334 //<<" "<< TheK[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1] 335 << 335 // <<" CK= " << CK[itet][ieta]<<" "<< CK[itet+1][ieta] 336 term /= material->GetTotNbOfAtomsPerVolume() << 336 //<<" "<< CK[itet][ieta+1]<<" "<< CK[itet+1][ieta+1]<<G4endl; 337 << 337 } 338 return term; << 338 //G4cout << "Kshell: tet= " << tet << " eta= " << eta << " C= " << corr 339 } << 339 //<< " itet,ieta= " << itet <<G4endl; 340 << 340 return corr; 341 //....oooOO0OOooo........oooOO0OOooo........oo << 341 } 342 << 342 343 G4double G4EmCorrections:: LShellCorrection(co << 343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 344 co << 344 345 co << 345 G4double G4EmCorrections::LShell(G4double tet, G4double eta) 346 { << 346 { 347 SetupKinematics(p, mat, e); << 347 G4double corr = 0.0; 348 G4double term = 0.0; << 348 349 for (G4int i = 0; i<numberOfElements; ++i) { << 349 G4double x = tet; 350 << 350 if(tet < TheL[0]) x = TheL[0]; 351 const G4double Z = (*theElementVector)[i]- << 351 G4int itet = Index(x, TheL, nL); 352 const G4int iz = (*theElementVector)[i]->G << 352 353 if(2 < iz) { << 353 // assimptotic case 354 const G4double Zeff = (iz < 10) ? Z - ZD << 354 if(eta >= Eta[nEtaL-1]) { 355 const G4double Z2= Zeff*Zeff; << 355 corr = (Value(x, TheL[itet], TheL[itet+1], UL[itet], UL[itet+1]) 356 const G4double eta = ba2/Z2; << 356 + Value(x, TheL[itet], TheL[itet+1], VL[itet], VL[itet+1])/eta 357 G4double tet = sThetaL->Value(Z); << 357 )/eta; 358 G4int nmax = std::min(4, G4AtomicShells: << 358 } else { 359 for (G4int j=1; j<nmax; ++j) { << 359 G4double y = eta; 360 G4int ne = G4AtomicShells::GetNumberOf << 360 if(eta < Eta[0]) y = Eta[0]; 361 if (15 >= iz) { << 361 G4int ieta = Index(y, Eta, nEtaL); 362 tet = (3 > j) ? 0.25*Z2*(1.0 + 5*Z2* << 362 corr = Value2(x, y, TheL[itet], TheL[itet+1], Eta[ieta], Eta[ieta+1], 363 0.25*Z2*(1.0 + Z2*alpha2/16.); << 363 CL[itet][ieta], CL[itet+1][ieta], CL[itet][ieta+1], CL[itet+1][ieta+1]); 364 } << 364 //G4cout << " x= " <<x<<" y= "<<y<<" tet= " <<TheL[itet] 365 //G4cout << " LShell: j= " << j << " n << 365 //<<" "<< TheL[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1] 366 // << " ThetaL= " << tet << G4en << 366 // <<" CL= " << CL[itet][ieta]<<" "<< CL[itet+1][ieta] 367 term += 0.125*ne*atomDensity[i]*LShell << 367 //<<" "<< CL[itet][ieta+1]<<" "<< CL[itet+1][ieta+1]<<G4endl; 368 } << 368 } 369 } << 369 // G4cout << "Lshell: tet= " << tet << " eta= " << eta << " C= " << corr << " itet= " << itet <<G4endl; 370 } << 370 return corr; 371 << 371 } 372 term /= material->GetTotNbOfAtomsPerVolume() << 372 373 << 373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 374 return term; << 374 375 } << 375 G4double G4EmCorrections::ShellCorrectionSTD(const G4ParticleDefinition* p, 376 << 376 const G4Material* mat, 377 //....oooOO0OOooo........oooOO0OOooo........oo << 377 G4double e) 378 << 378 { 379 G4double G4EmCorrections::KShell(const G4doubl << 379 SetupKinematics(p, mat, e); 380 { << 380 G4double taulim= 8.0*MeV/mass; 381 G4double corr = 0.0; << 381 G4double bg2lim= taulim * (taulim+2.0); 382 << 382 383 static const G4double TheK[20] = << 383 G4double* shellCorrectionVector = 384 {0.64, 0.65, 0.66, 0.68, 0.70, 0.72, 0.74, << 384 material->GetIonisation()->GetShellCorrectionVector(); 385 0.80, 0.82, 0.84, 0.85, 0.86, 0.88, 0.90, << 385 G4double sh = 0.0; 386 << 386 G4double x = 1.0; 387 << 387 G4double taul = material->GetIonisation()->GetTaul(); 388 G4double x = tet; << 388 389 G4int itet = 0; << 389 if ( bg2 >= bg2lim ) { 390 G4int ieta = 0; << 390 for (G4int k=0; k<3; k++) { 391 if(tet < TheK[0]) { << 391 x *= bg2 ; 392 x = TheK[0]; << 392 sh += shellCorrectionVector[k]/x; 393 } else if(tet > TheK[nK-1]) { << 393 } 394 x = TheK[nK-1]; << 394 395 itet = nK-2; << 395 } else { 396 } else { << 396 for (G4int k=0; k<3; k++) { 397 itet = Index(x, TheK, nK); << 397 x *= bg2lim ; 398 } << 398 sh += shellCorrectionVector[k]/x; 399 // asymptotic case << 399 } 400 if(eta >= Eta[nEtaK-1]) { << 400 sh *= std::log(tau/taul)/std::log(taulim/taul); 401 corr = << 401 } 402 (Value(x, TheK[itet], TheK[itet+1], UK[i << 402 sh *= 0.5; 403 Value(x, TheK[itet], TheK[itet+1], VK[i << 403 return sh; 404 Value(x, TheK[itet], TheK[itet+1], ZK[i << 404 } 405 } else { << 405 406 G4double y = eta; << 406 407 if(eta < Eta[0]) { << 407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 408 y = Eta[0]; << 408 409 } else { << 409 G4double G4EmCorrections::ShellCorrection(const G4ParticleDefinition* p, 410 ieta = Index(y, Eta, nEtaK); << 410 const G4Material* mat, 411 } << 411 G4double e) 412 corr = Value2(x, y, TheK[itet], TheK[itet+ << 412 { 413 CK[itet][ieta], CK[itet+1][i << 413 SetupKinematics(p, mat, e); 414 CK[itet][ieta+1], CK[itet+1] << 414 415 //G4cout << " x= " <<x<<" y= "<<y<<" tet << 415 G4double term = 0.0; 416 // <<" "<< TheK[itet+1]<<" eta= << 416 417 // <<" CK= " << CK[itet][ieta]<< << 417 for (G4int i = 0; i<numberOfElements; i++) { 418 // <<" "<< CK[itet][ieta+1]<<" " << 418 419 } << 419 G4double Z = (*theElementVector)[i]->GetZ(); 420 //G4cout << "Kshell: tet= " << tet << " eta << 420 G4int iz = G4int(Z); 421 // << " itet= " << itet << " ieta= << 421 G4double Z2= (Z-0.3)*(Z-0.3); 422 return corr; << 422 G4double f = 1.0; 423 } << 423 if(1 == iz) { 424 << 424 f = 0.5; 425 //....oooOO0OOooo........oooOO0OOooo........oo << 425 Z2 = 1.0; 426 << 426 } 427 G4double G4EmCorrections::LShell(const G4doubl << 427 G4double e0= 13.6*eV*Z2; 428 { << 428 term += f*atomDensity[i]*KShell(shells.GetBindingEnergy(iz,0)/e0,ba2/Z2)/Z; 429 G4double corr = 0.0; << 429 if(2 < iz) { 430 << 430 G4double Zeff = Z - ZD[10]; 431 static const G4double TheL[26] = << 431 if(iz < 10) Zeff = Z - ZD[iz]; 432 {0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.3 << 432 Z2= Zeff*Zeff; 433 0.42, 0.44, 0.45, 0.46, 0.48, 0.50, 0.5 << 433 e0= 13.6*eV*Z2*0.25; 434 0.58, 0.60, 0.62, 0.64, 0.65, 0.66}; << 434 f = 0.125; 435 << 435 G4double eta = ba2/Z2; 436 G4double x = tet; << 436 G4int ntot = shells.GetNumberOfShells(iz); 437 G4int itet = 0; << 437 G4int nmax = std::min(4, ntot); 438 G4int ieta = 0; << 438 G4double norm = 0.0; 439 if(tet < TheL[0]) { << 439 G4double eshell = 0.0; 440 x = TheL[0]; << 440 for(G4int j=1; j<nmax; j++) { 441 } else if(tet > TheL[nL-1]) { << 441 G4double x = G4double(shells.GetNumberOfElectrons(iz,j)); 442 x = TheL[nL-1]; << 442 G4double e1 = shells.GetBindingEnergy(iz,j); 443 itet = nL-2; << 443 norm += x; 444 } else { << 444 eshell += e1*x; 445 itet = Index(x, TheL, nL); << 445 term += f*x*atomDensity[i]*LShell(e1/e0,eta)/Z; 446 } << 446 } 447 << 447 if(10 < iz) { 448 // asymptotic case << 448 eshell /= norm; 449 if(eta >= Eta[nEtaL-1]) { << 449 G4double eeff = eshell*eta; 450 corr = (Value(x, TheL[itet], TheL[itet+1], << 450 for(G4int k=nmax; k<ntot; k++) { 451 + Value(x, TheL[itet], TheL[itet << 451 G4double x = G4double(shells.GetNumberOfElectrons(iz,k)); 452 )/eta; << 452 G4double e1 = shells.GetBindingEnergy(iz,k); 453 } else { << 453 term += f*x*atomDensity[i]*LShell(e1/e0,eeff/e1)/Z; 454 G4double y = eta; << 454 // term += f*x*atomDensity[i]*LShell(eshell/e0,eeff/e1)/Z; 455 if(eta < Eta[0]) { << 455 } 456 y = Eta[0]; << 456 /* 457 } else { << 457 if(28 >= iz) { 458 ieta = Index(y, Eta, nEtaL); << 458 term += f*(Z - 10.)*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z; 459 } << 459 } else if(32 >= iz) { 460 corr = Value2(x, y, TheL[itet], TheL[itet+ << 460 term += f*18.0*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z; 461 CL[itet][ieta], CL[itet+1][i << 461 } else if(60 >= iz) { 462 CL[itet][ieta+1], CL[itet+1] << 462 term += f*18.0*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z; 463 //G4cout << " x= " <<x<<" y= "<<y<<" tet << 463 term += f*(Z - 28.)*atomDensity[i]*LShell(eshell,HN[iz-33]*eta)/Z; 464 // <<" "<< TheL[itet+1]<<" eta= << 464 } else { 465 // <<" CL= " << CL[itet][ieta]<< << 465 term += f*18.0*atomDensity[i]*LShell(eshell,HM[53]*eta)/Z; 466 // <<" "<< CL[itet][ieta+1]<<" " << 466 term += f*32.0*atomDensity[i]*LShell(eshell,HN[30]*eta)/Z; 467 } << 467 term += f*(Z - 60.)*atomDensity[i]*LShell(eshell,150.*eta)/Z; 468 //G4cout<<"Lshell: tet= "<<tet<<" eta= "<<e << 468 } 469 // <<" ieta= "<<ieta<<" Corr= "<<cor << 469 */ 470 return corr; << 470 } 471 } << 471 } 472 << 472 //term += atomDensity[i]*MSH[iz]/(ba2*ba2); 473 //....oooOO0OOooo........oooOO0OOooo........oo << 473 } 474 << 474 475 G4double G4EmCorrections::ShellCorrectionSTD(c << 475 term /= material->GetTotNbOfAtomsPerVolume(); 476 c << 476 return term; 477 c << 477 } 478 { << 478 479 SetupKinematics(p, mat, e); << 479 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 480 G4double taulim= 8.0*MeV/mass; << 480 481 G4double bg2lim= taulim * (taulim+2.0); << 481 G4double G4EmCorrections::DensityCorrection(const G4ParticleDefinition* p, 482 << 482 const G4Material* mat, 483 G4double* shellCorrectionVector = << 483 G4double e) 484 material->GetIonisation()->GetShel << 484 { 485 G4double sh = 0.0; << 485 SetupKinematics(p, mat, e); 486 G4double x = 1.0; << 486 487 G4double taul = material->GetIonisation()-> << 487 G4double cden = material->GetIonisation()->GetCdensity(); 488 << 488 G4double mden = material->GetIonisation()->GetMdensity(); 489 if ( bg2 >= bg2lim ) { << 489 G4double aden = material->GetIonisation()->GetAdensity(); 490 for (G4int k=0; k<3; ++k) { << 490 G4double x0den = material->GetIonisation()->GetX0density(); 491 x *= bg2 ; << 491 G4double x1den = material->GetIonisation()->GetX1density(); 492 sh += shellCorrectionVector[k]/x; << 492 493 } << 493 G4double twoln10 = 2.0*std::log(10.0); 494 << 494 G4double dedx = 0.0; 495 } else { << 495 496 for (G4int k=0; k<3; ++k) { << 496 // density correction 497 x *= bg2lim ; << 497 G4double x = std::log(bg2)/twoln10; 498 sh += shellCorrectionVector[k]/x; << 498 if ( x >= x0den ) { 499 } << 499 dedx = twoln10*x - cden ; 500 sh *= G4Log(tau/taul)/G4Log(taulim/taul); << 500 if ( x < x1den ) dedx += aden*std::pow((x1den-x),mden) ; 501 } << 501 } 502 sh *= 0.5; << 502 503 return sh; << 503 return dedx; 504 } << 504 } 505 << 505 506 << 506 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 507 //....oooOO0OOooo........oooOO0OOooo........oo << 507 508 << 508 G4double G4EmCorrections::BarkasCorrection(const G4ParticleDefinition* p, 509 G4double G4EmCorrections::ShellCorrection(cons << 509 const G4Material* mat, 510 cons << 510 G4double e) 511 cons << 511 { 512 { << 512 // . Z^3 Barkas effect in the stopping power of matter for charged particles 513 SetupKinematics(p, mat, ekin); << 513 // J.C Ashley and R.H.Ritchie 514 G4double term = 0.0; << 514 // Physical review B Vol.5 No.7 1 April 1972 pp. 2393-2397 515 //G4cout << "### G4EmCorrections::ShellCorre << 515 // valid for kineticEnergy > 0.5 MeV 516 // << " " << ekin/MeV << " MeV " << << 516 517 for (G4int i = 0; i<numberOfElements; ++i) { << 517 SetupKinematics(p, mat, e); 518 << 518 G4double BarkasTerm = 0.0; 519 G4double res = 0.0; << 519 520 G4double res0 = 0.0; << 520 for (G4int i = 0; i<numberOfElements; i++) { 521 const G4double Z = (*theElementVector)[i]- << 521 522 const G4int iz = (*theElementVector)[i]->G << 522 G4double Z = (*theElementVector)[i]->GetZ(); 523 G4double Z2= (Z-0.3)*(Z-0.3); << 523 G4int iz = G4int(Z); 524 G4double f = 1.0; << 524 525 if(1 == iz) { << 525 G4double X = ba2 / Z; 526 f = 0.5; << 526 G4double b = 1.3; 527 Z2 = 1.0; << 527 if(1 == iz) { 528 } << 528 if(material->GetName() == "G4_lH2") b = 0.6; 529 G4double eta = ba2/Z2; << 529 else b = 1.8; 530 G4double tet = (11 < iz) ? sThetaK->Value( << 530 } 531 res0 = f*KShell(tet,eta); << 531 else if(2 == iz) b = 0.6; 532 res += res0; << 532 else if(10 >= iz) b = 1.8; 533 //G4cout << " Z= " << iz << " Shell 0" << << 533 else if(17 >= iz) b = 1.4; 534 // << " eta= " << eta << " resK= " << 534 else if(18 == iz) b = 1.8; 535 << 535 else if(25 >= iz) b = 1.4; 536 if(2 < iz) { << 536 else if(50 >= iz) b = 1.35; 537 const G4double Zeff = (iz < 10) ? Z - ZD << 537 538 Z2= Zeff*Zeff; << 538 G4double W = b/std::sqrt(X); 539 eta = ba2/Z2; << 539 540 tet = sThetaL->Value(Z); << 540 G4double val; 541 f = 0.125; << 541 if(W <= engBarkas[0]) val = corBarkas[0]; 542 const G4int ntot = G4AtomicShells::GetNu << 542 else if(W >= engBarkas[46]) val = corBarkas[46]*engBarkas[46]/W; 543 const G4int nmax = std::min(4, ntot); << 543 else { 544 G4double norm = 0.0; << 544 G4int iw = Index(W, engBarkas, 47); 545 G4double eshell = 0.0; << 545 val = Value(W, engBarkas[iw], engBarkas[iw+1], 546 for(G4int j=1; j<nmax; ++j) { << 546 corBarkas[iw], corBarkas[iw+1]); 547 G4int ne = G4AtomicShells::GetNumberOf << 547 } 548 if(15 >= iz) { << 548 // G4cout << "i= " << i << " b= " << b << " W= " << W 549 tet = (3 > j) ? 0.25*Z2*(1.0 + 5*Z2* << 549 // << " Z= " << Z << " X= " << X << " val= " << val<< G4endl; 550 0.25*Z2*(1.0 + Z2*alpha2/16.); << 550 BarkasTerm += val*atomDensity[i] / (std::sqrt(Z*X)*X); 551 } << 551 } 552 norm += ne; << 552 553 eshell += tet*ne; << 553 BarkasTerm *= 1.29*charge/material->GetTotNbOfAtomsPerVolume(); 554 res0 = f*ne*LShell(tet,eta); << 554 return BarkasTerm; 555 res += res0; << 555 } 556 //G4cout << " Zeff= " << Zeff << " She << 556 557 // << " tet= " << tet << " eta= << 557 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 558 // << " resL= " << res0 << G4en << 558 559 } << 559 G4double G4EmCorrections::BlochCorrection(const G4ParticleDefinition* p, 560 if(ntot > nmax) { << 560 const G4Material* mat, 561 if (norm > 0.0) { norm = 1.0/norm; } << 561 G4double e) 562 eshell *= norm; << 562 { 563 << 563 SetupKinematics(p, mat, e); 564 static const G4double HM[53] = { << 564 565 12.0, 12.0, 12.0, 12.0, 11.9, 11.7, 11.5, << 565 G4double y2 = q2/ba2; 566 10.0, 9.51, 8.97, 8.52, 8.03, 7.46, 6.95, << 566 567 5.61, 5.39, 5.19, 5.01, 4.86, 4.72, 4.62, << 567 G4double term = 1.0/(1.0 + y2); 568 4.32, 4.26, 4.20, 4.15, 4.1, 4.04, 4.00, << 568 G4double del; 569 3.90, 3.89, 3.89, 3.88, 3.88, 3.88, 3.88, << 569 G4double j = 1.0; 570 3.90, 3.92, 3.93 }; << 570 do { 571 static const G4double HN[31] = { << 571 j += 1.0; 572 75.5, 61.9, 52.2, 45.1, 39.6, 35.4, 31.9, << 572 del = 1.0/(j* (j*j + y2)); 573 24.5, 23.6, 22.7, 22.0, 21.4, 20.9, 20.5, << 573 term += del; 574 19.5, 19.3, 19.2, 19.1, 18.4, 18.8, 18.7, << 574 } while (del > 0.01*term); 575 18.2}; << 575 576 << 576 return -y2*term; 577 // Add M-shell << 577 } 578 if(28 > iz) { << 578 579 res += f*(iz - 10)*LShell(eshell,HM[ << 579 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 580 } else if(63 > iz) { << 580 581 res += f*18*LShell(eshell,HM[iz-11]* << 581 G4double G4EmCorrections::MottCorrection(const G4ParticleDefinition* p, 582 } else { << 582 const G4Material* mat, 583 res += f*18*LShell(eshell,HM[52]*eta << 583 G4double e) 584 } << 584 { 585 // Add N-shell << 585 SetupKinematics(p, mat, e); 586 if(32 < iz) { << 586 G4double mterm = CLHEP::pi*fine_structure_const*beta*charge; 587 if(60 > iz) { << 587 return mterm; 588 res += f*(iz - 28)*LShell(eshell,H << 588 } 589 } else if(63 > iz) { << 589 590 res += 4*LShell(eshell,HN[iz-33]*e << 590 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 591 } else { << 591 592 res += 4*LShell(eshell,HN[30]*eta) << 592 G4double G4EmCorrections::NuclearDEDX(const G4ParticleDefinition* p, 593 } << 593 const G4Material* mat, 594 // Add O-P-shells << 594 G4double e, 595 if(60 < iz) { << 595 G4bool fluct) 596 res += f*(iz - 60)*LShell(eshell,1 << 596 { 597 } << 597 G4double nloss = 0.0; 598 } << 598 if(e <= 0.0) return nloss; 599 } << 599 SetupKinematics(p, mat, e); 600 } << 600 601 term += res*atomDensity[i]/Z; << 601 lossFlucFlag = fluct; 602 } << 602 603 << 603 // Projectile nucleus 604 term /= material->GetTotNbOfAtomsPerVolume() << 604 G4double z1 = std::abs(particle->GetPDGCharge()/eplus); 605 //G4cout << "##Shell Correction=" << term << << 605 G4double m1 = mass/amu_c2; 606 return term; << 606 607 } << 607 // loop for the elements in the material 608 << 608 for (G4int iel=0; iel<numberOfElements; iel++) { 609 //....oooOO0OOooo........oooOO0OOooo........oo << 609 const G4Element* element = (*theElementVector)[iel] ; 610 << 610 G4double z2 = element->GetZ(); 611 G4double G4EmCorrections::DensityCorrection(co << 611 G4double m2 = element->GetA()*mole/g ; 612 co << 612 nloss += (NuclearStoppingPower(kinEnergy, z1, z2, m1, m2)) 613 co << 613 * atomDensity[iel] ; 614 { << 614 } 615 SetupKinematics(p, mat, e); << 615 nloss *= theZieglerFactor; 616 << 616 return nloss; 617 G4double cden = material->GetIonisation()-> << 617 } 618 G4double mden = material->GetIonisation()-> << 618 619 G4double aden = material->GetIonisation()-> << 619 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 620 G4double x0den = material->GetIonisation()-> << 620 621 G4double x1den = material->GetIonisation()-> << 621 G4double G4EmCorrections::NuclearStoppingPower(G4double kineticEnergy, 622 << 622 G4double z1, G4double z2, 623 G4double dedx = 0.0; << 623 G4double m1, G4double m2) 624 << 624 { 625 // density correction << 625 G4double energy = kineticEnergy/keV ; // energy in keV 626 static const G4double twoln10 = 2.0*G4Log(10 << 626 G4double nloss = 0.0; 627 G4double x = G4Log(bg2)/twoln10; << 627 628 if ( x >= x0den ) { << 628 G4double rm; 629 dedx = twoln10*x - cden ; << 629 if(z1 > 1.5) rm = (m1 + m2) * ( Z23[G4int(z1)] + Z23[G4int(z2)] ) ; 630 if ( x < x1den ) dedx += aden*G4Exp(G4Log( << 630 else rm = (m1 + m2) * nist->GetZ13(G4int(z2)); 631 } << 631 632 << 632 G4double er = 32.536 * m2 * energy / ( z1 * z2 * rm ) ; // reduced energy 633 return dedx; << 633 634 } << 634 if (er >= ed[0]) nloss = a[0]; 635 << 635 else { 636 //....oooOO0OOooo........oooOO0OOooo........oo << 636 // the table is inverse in energy 637 << 637 for (G4int i=102; i>=0; i--) 638 G4double G4EmCorrections::BarkasCorrection(con << 638 { 639 con << 639 if (er <= ed[i]) { 640 con << 640 nloss = (a[i] - a[i+1])*(er - ed[i+1])/(ed[i] - ed[i+1]) + a[i+1]; 641 con << 641 break; 642 { << 642 } 643 // . Z^3 Barkas effect in the stopping power << 643 } 644 // J.C Ashley and R.H.Ritchie << 644 } 645 // Physical review B Vol.5 No.7 1 April 19 << 645 646 // valid for kineticEnergy > 0.5 MeV << 646 // Stragling 647 << 647 if(lossFlucFlag) { 648 if (!isInitialized) { SetupKinematics(p, mat << 648 // G4double sig = 4.0 * m1 * m2 / ((m1 + m2)*(m1 + m2)* 649 G4double BarkasTerm = 0.0; << 649 // (4.0 + 0.197*std::pow(er,-1.6991)+6.584*std::pow(er,-1.0494))) ; 650 << 650 G4double sig = 4.0 * m1 * m2 / ((m1 + m2)*(m1 + m2)* 651 for (G4int i = 0; i<numberOfElements; ++i) { << 651 (4.0 + 0.197/(er*er) + 6.584/er)); 652 << 652 653 const G4int iz = (*theElementVector)[i]->G << 653 nloss *= G4RandGauss::shoot(1.0,sig) ; 654 if(iz == 47) { << 654 } 655 BarkasTerm += atomDensity[i]*0.006812*G4 << 655 656 } else if(iz >= 64) { << 656 nloss *= 8.462 * z1 * z2 * m1 / rm ; // Return to [ev/(10^15 atoms/cm^2] 657 BarkasTerm += atomDensity[i]*0.002833*G4 << 657 658 } else { << 658 if ( nloss < 0.0) nloss = 0.0 ; 659 << 659 660 const G4double Z = (*theElementVector)[i << 660 return nloss; 661 const G4double X = ba2 / Z; << 661 } 662 G4double b = 1.3; << 662 663 if(1 == iz) { b = (material->GetName() = << 663 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 664 else if(2 == iz) { b = 0.6; } << 664 665 else if(10 >= iz) { b = 1.8; } << 665 G4double G4EmCorrections::EffectiveChargeCorrection(const G4ParticleDefinition* p, 666 else if(17 >= iz) { b = 1.4; } << 666 const G4Material* mat, 667 else if(18 == iz) { b = 1.8; } << 667 G4double ekin) 668 else if(25 >= iz) { b = 1.4; } << 668 { 669 else if(50 >= iz) { b = 1.35;} << 669 G4double factor = 1.0; 670 << 670 if(p->GetPDGCharge() <= 2.5*eplus || nIons <= 0) return factor; 671 const G4double W = b/std::sqrt(X); << 671 /* 672 << 672 if(verbose > 1) { 673 G4double val = sBarkasCorr->Value(W, idx << 673 G4cout << "EffectiveChargeCorrection: " << p->GetParticleName() 674 if (W > sWmaxBarkas) { val *= (sWmaxBark << 674 << " in " << mat->GetName() 675 // G4cout << "i= " << i << " b= " << << 675 << " ekin(MeV)= " << ekin/MeV << G4endl; 676 // << " Z= " << Z << " X= " << X << " va << 676 } 677 BarkasTerm += val*atomDensity[i] / (std: << 677 */ 678 } << 678 if(p != curParticle || mat != curMaterial) { 679 } << 679 curParticle = p; 680 << 680 curMaterial = mat; 681 BarkasTerm *= 1.29*charge/material->GetTotNb << 681 curVector = 0; 682 << 682 currentZ = p->GetAtomicNumber(); 683 return BarkasTerm; << 683 if(verbose > 1) { 684 } << 684 G4cout << "G4EmCorrections::EffectiveChargeCorrection: Zion= " 685 << 685 << currentZ << " Aion= " << p->GetPDGMass()/amu_c2 << G4endl; 686 //....oooOO0OOooo........oooOO0OOooo........oo << 686 } 687 << 687 massFactor = proton_mass_c2/p->GetPDGMass(); 688 G4double G4EmCorrections::BlochCorrection(cons << 688 idx = -1; 689 cons << 689 690 cons << 690 for(G4int i=0; i<nIons; i++) { 691 cons << 691 if(materialList[i] == mat && currentZ == Zion[i]) { 692 { << 692 idx = i; 693 if (!isInitialized) { SetupKinematics(p, mat << 693 break; 694 << 694 } 695 G4double y2 = q2/ba2; << 695 } 696 << 696 // G4cout << " idx= " << idx << " dz= " << dz << G4endl; 697 G4double term = 1.0/(1.0 + y2); << 697 if(idx >= 0) { 698 G4double del; << 698 if(!ionList[idx]) BuildCorrectionVector(); 699 G4double j = 1.0; << 699 if(ionList[idx]) curVector = stopData[idx]; 700 do { << 700 } else { return factor; } 701 j += 1.0; << 701 } 702 del = 1.0/(j* (j*j + y2)); << 702 if(curVector) { 703 term += del; << 703 G4bool b; 704 // Loop checking, 03-Aug-2015, Vladimir Iv << 704 factor = curVector->GetValue(ekin*massFactor,b); 705 } while (del > 0.01*term); << 705 if(verbose > 1) { 706 << 706 G4cout << "E= " << ekin << " factor= " << factor << " massfactor= " 707 return -y2*term; << 707 << massFactor << G4endl; 708 } << 708 } 709 << 709 } 710 //....oooOO0OOooo........oooOO0OOooo........oo << 710 return factor; 711 << 711 } 712 G4double G4EmCorrections::MottCorrection(const << 712 713 const << 713 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 714 const << 714 715 const << 715 void G4EmCorrections::AddStoppingData(G4int Z, G4int A, 716 { << 716 const G4String& mname, 717 if (!isInitialized) { SetupKinematics(p, mat << 717 G4PhysicsVector* dVector) 718 return CLHEP::pi*CLHEP::fine_structure_const << 718 { 719 } << 719 G4int i = 0; 720 << 720 for(; i<nIons; i++) { 721 //....oooOO0OOooo........oooOO0OOooo........oo << 721 if(Z == Zion[i] && A == Aion[i] && mname == materialName[i]) break; 722 << 722 } 723 G4double << 723 if(i == nIons) { 724 G4EmCorrections::EffectiveChargeCorrection(con << 724 Zion.push_back(Z); 725 con << 725 Aion.push_back(A); 726 con << 726 materialName.push_back(mname); 727 { << 727 materialList.push_back(0); 728 G4double factor = 1.0; << 728 ionList.push_back(0); 729 if(p->GetPDGCharge() <= 2.5*CLHEP::eplus || << 729 stopData.push_back(dVector); 730 << 730 nIons++; 731 if(verbose > 1) { << 731 if(verbose>1) { 732 G4cout << "EffectiveChargeCorrection: " << << 732 G4cout << "AddStoppingData Z= " << Z << " A= " << A << " " << mname 733 << " in " << mat->GetName() << 733 << " idx= " << i << G4endl; 734 << " ekin(MeV)= " << ekin/MeV << G4 << 734 } 735 } << 735 } 736 << 736 } 737 if(p != curParticle || mat != curMaterial) { << 737 738 curParticle = p; << 738 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 739 curMaterial = mat; << 739 740 curVector = nullptr; << 740 void G4EmCorrections::BuildCorrectionVector() 741 currentZ = p->GetAtomicNumber(); << 741 { 742 if(verbose > 1) { << 742 if(!ionLEModel || !ionHEModel) { 743 G4cout << "G4EmCorrections::EffectiveCha << 743 return; 744 << currentZ << " Aion= " << p->Ge << 744 } 745 } << 745 746 massFactor = CLHEP::proton_mass_c2/p->GetP << 746 const G4ParticleDefinition* ion = curParticle; 747 idx = -1; << 747 G4int Z = Zion[idx]; 748 << 748 if(currentZ != Z) { 749 for(G4int i=0; i<nIons; ++i) { << 749 ion = G4ParticleTable::GetParticleTable()->FindIon(Z, Aion[idx], 0, Z); 750 if(materialList[i] == mat && currentZ == << 750 } 751 idx = i; << 751 //G4cout << "BuildCorrectionVector: idx= " << idx << " Z= " << Z 752 break; << 752 // << " curZ= " << currentZ << G4endl; 753 } << 753 754 } << 754 // G4double A = nist->GetAtomicMassAmu(Z); 755 //G4cout << " idx= " << idx << " dz= " << << 755 G4double A = G4double(ion->GetBaryonNumber()); 756 if(idx >= 0) { << 756 G4PhysicsVector* v = stopData[idx]; 757 if(nullptr == ionList[idx]) { BuildCorre << 757 758 curVector = stopData[idx]; << 758 const G4ParticleDefinition* p = G4GenericIon::GenericIon(); 759 } else { << 759 G4double massRatio = proton_mass_c2/ion->GetPDGMass(); 760 return factor; << 760 761 } << 761 if(verbose>1) { 762 } << 762 G4cout << "BuildCorrectionVector: Stopping for " 763 if(nullptr != curVector) { << 763 << curParticle->GetParticleName() << " in " 764 factor = curVector->Value(ekin*massFactor) << 764 << materialName[idx] << " Ion Z= " << Z << " A= " << A 765 if(verbose > 1) { << 765 << " massRatio= " << massRatio << G4endl; 766 G4cout << "E= " << ekin << " factor= " < << 766 } 767 << massFactor << G4endl; << 767 G4bool b; 768 } << 768 769 } << 769 G4PhysicsLogVector* vv = 770 return factor; << 770 new G4PhysicsLogVector(eCorrMin,eCorrMax,nbinCorr); 771 } << 771 vv->SetSpline(true); 772 << 772 G4double e, eion, dedx, dedx1; 773 //....oooOO0OOooo........oooOO0OOooo........oo << 773 G4double eth0 = v->GetLowEdgeEnergy(0); 774 << 774 G4double escal = eth/massRatio; 775 void G4EmCorrections::AddStoppingData(const G4 << 775 G4double qe = 776 const G4 << 776 effCharge.EffectiveChargeSquareRatio(ion, curMaterial, escal); 777 G4Physic << 777 G4double dedxt = 778 { << 778 ionLEModel->ComputeDEDXPerVolume(curMaterial, p, eth, eth)*qe; 779 G4int i = 0; << 779 G4double dedx1t = 780 for(; i<nIons; ++i) { << 780 ionHEModel->ComputeDEDXPerVolume(curMaterial, p, eth, eth)*qe 781 if(Z == Zion[i] && A == Aion[i] && mname = << 781 + ComputeIonCorrections(curParticle, curMaterial, escal); 782 } << 782 G4double rest = escal*(dedxt - dedx1t); 783 if(i == nIons) { << 783 //G4cout << "Escal(MeV)= "<<escal<<" dedxt0= " <<dedxt 784 Zion.push_back(Z); << 784 // << " dedxt1= " << dedx1t << G4endl; 785 Aion.push_back(A); << 785 786 materialName.push_back(mname); << 786 for(G4int i=0; i<nbinCorr; i++) { 787 materialList.push_back(nullptr); << 787 e = vv->GetLowEdgeEnergy(i); 788 ionList.push_back(nullptr); << 788 escal = e/massRatio; 789 stopData.push_back(dVector); << 789 eion = escal/A; 790 nIons++; << 790 if(eion <= eth0) { 791 if(verbose > 1) { << 791 dedx = v->GetValue(eth0, b)*std::sqrt(eion/eth0); 792 G4cout << "AddStoppingData Z= " << Z << << 792 } else { 793 << " idx= " << i << G4endl; << 793 dedx = v->GetValue(eion, b); 794 } << 794 } 795 } << 795 qe = effCharge.EffectiveChargeSquareRatio(curParticle,curMaterial,escal); 796 } << 796 if(e <= eth) { 797 << 797 dedx1 = ionLEModel->ComputeDEDXPerVolume(curMaterial, p, e, e)*qe; 798 //....oooOO0OOooo........oooOO0OOooo........oo << 798 } else { 799 << 799 dedx1 = ionHEModel->ComputeDEDXPerVolume(curMaterial, p, e, e)*qe + 800 void G4EmCorrections::BuildCorrectionVector() << 800 ComputeIonCorrections(curParticle, curMaterial, escal) + rest/escal; 801 { << 801 } 802 if(nullptr == ionLEModel || nullptr == ionHE << 802 vv->PutValue(i, dedx/dedx1); 803 return; << 803 if(verbose>1) { 804 } << 804 G4cout << " E(meV)= " << e/MeV << " Correction= " << dedx/dedx1 805 << 805 << " " << dedx << " " << dedx1 806 const G4ParticleDefinition* ion = curParticl << 806 << " massF= " << massFactor << G4endl; 807 const G4ParticleDefinition* gion = G4Generic << 807 } 808 G4int Z = Zion[idx]; << 808 } 809 G4double A = Aion[idx]; << 809 delete v; 810 G4PhysicsVector* v = stopData[idx]; << 810 ionList[idx] = ion; 811 << 811 stopData[idx] = vv; 812 if(verbose > 1) { << 812 if(verbose>1) G4cout << "End data set " << G4endl; 813 G4cout << "BuildCorrectionVector: Stopping << 813 } 814 << curParticle->GetParticleName() < << 814 815 << materialName[idx] << " Ion Z= " << 815 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 816 << " massFactor= " << massFactor << << 816 817 G4cout << " Nbins=" << nbinCorr << " Em << 817 void G4EmCorrections::InitialiseForNewRun() 818 << " Emax(MeV)=" << eCorrMax << " ion: " << 818 { 819 << ion->GetParticleName() << G4endl << 819 G4ProductionCutsTable* tb = G4ProductionCutsTable::GetProductionCutsTable(); 820 } << 820 ncouples = tb->GetTableSize(); 821 << 821 if(currmat.size() != ncouples) { 822 auto vv = new G4PhysicsLogVector(eCorrMin,eC << 822 currmat.resize(ncouples); 823 const G4double eth0 = v->Energy(0); << 823 size_t i; 824 const G4double escal = eth/massFactor; << 824 for(i=0; i<100; i++) {thcorr[i].clear();} 825 G4double qe = << 825 for(i=0; i<ncouples; i++) { 826 effCharge.EffectiveChargeSquareRatio(curPa << 826 currmat[i] = tb->GetMaterialCutsCouple(i)->GetMaterial(); 827 const G4double dedxt = << 827 G4String nam = currmat[i]->GetName(); 828 ionLEModel->ComputeDEDXPerVolume(curMateri << 828 for(G4int j=0; j<nIons; j++) { 829 const G4double dedx1t = << 829 if(nam == materialName[j]) { materialList[j] = currmat[i]; } 830 ionHEModel->ComputeDEDXPerVolume(curMateri << 830 } 831 + ComputeIonCorrections(curParticle, curMa << 831 } 832 const G4double rest = escal*(dedxt - dedx1t) << 832 } 833 if(verbose > 1) { << 833 } 834 G4cout << "Escal(MeV)= " << escal << " qe= << 834 835 << " dedxt= " << dedxt << " dedx1t= << 835 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 836 } << 836 837 for(G4int i=0; i<=nbinCorr; ++i) { << 837 void G4EmCorrections::Initialise() 838 // energy in the local table (GenericIon) << 838 { 839 G4double e = vv->Energy(i); << 839 // . Z^3 Barkas effect in the stopping power of matter for charged particles 840 // energy of the real ion << 840 // J.C Ashley and R.H.Ritchie 841 G4double eion = e/massFactor; << 841 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397 842 // energy in the imput stopping data vecto << 842 G4int i, j; 843 G4double e1 = eion/A; << 843 const G4double fTable[47][2] = { 844 G4double dedx = (e1 < eth0) << 844 { 0.02, 21.5}, 845 ? v->Value(eth0)*std::sqrt(e1/eth0) : v- << 845 { 0.03, 20.0}, 846 qe = effCharge.EffectiveChargeSquareRatio( << 846 { 0.04, 18.0}, 847 G4double dedx1 = (e <= eth) << 847 { 0.05, 15.6}, 848 ? ionLEModel->ComputeDEDXPerVolume(curMa << 848 { 0.06, 15.0}, 849 : ionHEModel->ComputeDEDXPerVolume(curMa << 849 { 0.07, 14.0}, 850 ComputeIonCorrections(curParticle, curMa << 850 { 0.08, 13.5}, 851 vv->PutValue(i, dedx/dedx1); << 851 { 0.09, 13.}, 852 if(verbose > 1) { << 852 { 0.1, 12.2}, 853 G4cout << "E(MeV)=" << e/CLHEP::MeV << " << 853 { 0.2, 9.25}, 854 << " e1=" << e1 << " dedxRatio= " << de << 854 { 0.3, 7.0}, 855 << " dedx=" << dedx << " dedx1=" << 855 { 0.4, 6.0}, 856 << " qe=" << qe << " rest/eion=" << 856 { 0.5, 4.5}, 857 } << 857 { 0.6, 3.5}, 858 } << 858 { 0.7, 3.0}, 859 delete v; << 859 { 0.8, 2.5}, 860 ionList[idx] = ion; << 860 { 0.9, 2.0}, 861 stopData[idx] = vv; << 861 { 1.0, 1.7}, 862 if(verbose > 1) { G4cout << "End data set " << 862 { 1.2, 1.2}, 863 } << 863 { 1.3, 1.0}, 864 << 864 { 1.4, 0.86}, 865 //....oooOO0OOooo........oooOO0OOooo........oo << 865 { 1.5, 0.7}, 866 << 866 { 1.6, 0.61}, 867 void G4EmCorrections::InitialiseForNewRun() << 867 { 1.7, 0.52}, 868 { << 868 { 1.8, 0.5}, 869 G4ProductionCutsTable* tb = G4ProductionCuts << 869 { 1.9, 0.43}, 870 ncouples = tb->GetTableSize(); << 870 { 2.0, 0.42}, 871 if(currmat.size() != ncouples) { << 871 { 2.1, 0.3}, 872 currmat.resize(ncouples); << 872 { 2.4, 0.2}, 873 for(auto it = thcorr.begin(); it != thcorr << 873 { 3.0, 0.13}, 874 (it->second).clear(); << 874 { 3.08, 0.1}, 875 } << 875 { 3.1, 0.09}, 876 thcorr.clear(); << 876 { 3.3, 0.08}, 877 for(std::size_t i=0; i<ncouples; ++i) { << 877 { 3.5, 0.07}, 878 currmat[i] = tb->GetMaterialCutsCouple(( << 878 { 3.8, 0.06}, 879 G4String nam = currmat[i]->GetName(); << 879 { 4.0, 0.051}, 880 for(G4int j=0; j<nIons; ++j) { << 880 { 4.1, 0.04}, 881 if(nam == materialName[j]) { materialL << 881 { 4.8, 0.03}, 882 } << 882 { 5.0, 0.024}, 883 } << 883 { 5.1, 0.02}, 884 } << 884 { 6.0, 0.013}, 885 } << 885 { 6.5, 0.01}, 886 << 886 { 7.0, 0.009}, 887 //....oooOO0OOooo........oooOO0OOooo........oo << 887 { 7.1, 0.008}, 888 << 888 { 8.0, 0.006}, 889 void G4EmCorrections::Initialise() << 889 { 9.0, 0.0032}, 890 { << 890 { 10.0, 0.0025} }; 891 // Z^3 Barkas effect in the stopping power o << 891 892 // J.C Ashley and R.H.Ritchie << 892 for(i=0; i<47; i++) { 893 // Physical review B Vol.5 No.7 1 April 1972 << 893 engBarkas[i] = fTable[i][0]; 894 // G.S. Khandelwal Nucl. Phys. A116(1968)97 << 894 corBarkas[i] = fTable[i][1]; 895 // "Shell corrections for K- and L- electron << 895 } 896 << 896 897 G4int i, j; << 897 const G4double nuca[104][2] = { 898 static const G4double fTable[47][2] = { << 898 { 1.0E+8, 5.831E-8}, 899 { 0.02, 21.5}, << 899 { 8.0E+7, 7.288E-8}, 900 { 0.03, 20.0}, << 900 { 6.0E+7, 9.719E-8}, 901 { 0.04, 18.0}, << 901 { 5.0E+7, 1.166E-7}, 902 { 0.05, 15.6}, << 902 { 4.0E+7, 1.457E-7}, 903 { 0.06, 15.0}, << 903 { 3.0E+7, 1.942E-7}, 904 { 0.07, 14.0}, << 904 { 2.0E+7, 2.916E-7}, 905 { 0.08, 13.5}, << 905 { 1.5E+7, 3.887E-7}, 906 { 0.09, 13.}, << 906 907 { 0.1, 12.2}, << 907 { 1.0E+7, 5.833E-7}, 908 { 0.2, 9.25}, << 908 { 8.0E+6, 7.287E-7}, 909 { 0.3, 7.0}, << 909 { 6.0E+6, 9.712E-7}, 910 { 0.4, 6.0}, << 910 { 5.0E+6, 1.166E-6}, 911 { 0.5, 4.5}, << 911 { 4.0E+6, 1.457E-6}, 912 { 0.6, 3.5}, << 912 { 3.0E+6, 1.941E-6}, 913 { 0.7, 3.0}, << 913 { 2.0E+6, 2.911E-6}, 914 { 0.8, 2.5}, << 914 { 1.5E+6, 3.878E-6}, 915 { 0.9, 2.0}, << 915 916 { 1.0, 1.7}, << 916 { 1.0E+6, 5.810E-6}, 917 { 1.2, 1.2}, << 917 { 8.0E+5, 7.262E-6}, 918 { 1.3, 1.0}, << 918 { 6.0E+5, 9.663E-6}, 919 { 1.4, 0.86}, << 919 { 5.0E+5, 1.157E-5}, 920 { 1.5, 0.7}, << 920 { 4.0E+5, 1.442E-5}, 921 { 1.6, 0.61}, << 921 { 3.0E+5, 1.913E-5}, 922 { 1.7, 0.52}, << 922 { 2.0E+5, 2.845E-5}, 923 { 1.8, 0.5}, << 923 { 1.5E+5, 3.762E-5}, 924 { 1.9, 0.43}, << 924 925 { 2.0, 0.42}, << 925 { 1.0E+5, 5.554E-5}, 926 { 2.1, 0.3}, << 926 { 8.0E+4, 6.866E-5}, 927 { 2.4, 0.2}, << 927 { 6.0E+4, 9.020E-5}, 928 { 3.0, 0.13}, << 928 { 5.0E+4, 1.070E-4}, 929 { 3.08, 0.1}, << 929 { 4.0E+4, 1.319E-4}, 930 { 3.1, 0.09}, << 930 { 3.0E+4, 1.722E-4}, 931 { 3.3, 0.08}, << 931 { 2.0E+4, 2.499E-4}, 932 { 3.5, 0.07}, << 932 { 1.5E+4, 3.248E-4}, 933 { 3.8, 0.06}, << 933 934 { 4.0, 0.051}, << 934 { 1.0E+4, 4.688E-4}, 935 { 4.1, 0.04}, << 935 { 8.0E+3, 5.729E-4}, 936 { 4.8, 0.03}, << 936 { 6.0E+3, 7.411E-4}, 937 { 5.0, 0.024}, << 937 { 5.0E+3, 8.718E-4}, 938 { 5.1, 0.02}, << 938 { 4.0E+3, 1.063E-3}, 939 { 6.0, 0.013}, << 939 { 3.0E+3, 1.370E-3}, 940 { 6.5, 0.01}, << 940 { 2.0E+3, 1.955E-3}, 941 { 7.0, 0.009}, << 941 { 1.5E+3, 2.511E-3}, 942 { 7.1, 0.008}, << 942 943 { 8.0, 0.006}, << 943 { 1.0E+3, 3.563E-3}, 944 { 9.0, 0.0032}, << 944 { 8.0E+2, 4.314E-3}, 945 { 10.0, 0.0025} }; << 945 { 6.0E+2, 5.511E-3}, 946 << 946 { 5.0E+2, 6.430E-3}, 947 sBarkasCorr = new G4PhysicsFreeVector(47, fa << 947 { 4.0E+2, 7.756E-3}, 948 for(i=0; i<47; ++i) { sBarkasCorr->PutValues << 948 { 3.0E+2, 9.855E-3}, 949 << 949 { 2.0E+2, 1.375E-2}, 950 const G4double SK[20] = {1.9477, 1.9232, 1.8 << 950 { 1.5E+2, 1.736E-2}, 951 1.7754, 1.7396, 1.7 << 951 952 1.6461, 1.6189, 1.5 << 952 { 1.0E+2, 2.395E-2}, 953 1.5467, 1.5254, 1.5 << 953 { 8.0E+1, 2.850E-2}, 954 const G4double TK[20] = {2.5222, 2.5125, 2.5 << 954 { 6.0E+1, 3.552E-2}, 955 2.4388, 2.4163, 2.4 << 955 { 5.0E+1, 4.073E-2}, 956 2.3466, 2.3229, 2.2 << 956 { 4.0E+1, 4.802E-2}, 957 2.2515, 2.2277, 2.2 << 957 { 3.0E+1, 5.904E-2}, 958 << 958 { 1.5E+1, 9.426E-2}, 959 const G4double SL[26] = {15.3343, 13.9389, 1 << 959 960 10.3424, 10.0371, << 960 { 1.0E+1, 1.210E-1}, 961 8.4114, 8.0683, << 961 { 8.0E+0, 1.377E-1}, 962 7.2506, 7.0327, << 962 { 6.0E+0, 1.611E-1}, 963 6.4969, 6.3498, << 963 { 5.0E+0, 1.768E-1}, 964 const G4double TL[26] = {35.0669, 33.4344, 3 << 964 { 4.0E+0, 1.968E-1}, 965 28.6128, 28.1449, 2 << 965 { 3.0E+0, 2.235E-1}, 966 25.4058, 24.7587, 2 << 966 { 2.0E+0, 2.613E-1}, 967 23.0771, 22.5880, 2 << 967 { 1.5E+0, 2.871E-1}, 968 21.2872, 20.9006, 2 << 968 969 << 969 { 1.0E+0, 3.199E-1}, 970 const G4double bk1[29][11] = { << 970 { 8.0E-1, 3.354E-1}, 971 {0.005, 1.34782E-8, 1.46132E-8, 1.72179E-8, << 971 { 6.0E-1, 3.523E-1}, 972 {0.007, 6.87555E-8, 7.44373E-8, 8.74397E-8, << 972 { 5.0E-1, 3.609E-1}, 973 {0.01, 3.78413E-7, 4.08831E-7, 4.78154E-7, 5 << 973 { 4.0E-1, 3.693E-1}, 974 {0.015, 2.53200E-6, 2.72664E-6, 3.16738E-6, << 974 { 3.0E-1, 3.766E-1}, 975 {0.02, 9.43891E-6, 1.01339E-5, 1.16984E-5, 1 << 975 { 2.0E-1, 3.803E-1}, 976 {0.03, 5.67088E-5, 6.05576E-5, 6.91311E-5, 7 << 976 { 1.5E-1, 3.788E-1}, 977 {0.04, 1.91576E-4, 2.03626E-4, 2.30230E-4, 2 << 977 978 {0.05, 4.74226E-4, 5.02006E-4, 5.62872E-4, 6 << 978 { 1.0E-1, 3.711E-1}, 979 {0.06, 9.67285E-4, 1.02029E-3, 1.13566E-3, 1 << 979 { 8.0E-2, 3.644E-1}, 980 {0.08, 2.81537E-3, 2.95200E-3, 3.24599E-3, 3 << 980 { 6.0E-2, 3.530E-1}, 981 {0.1, 6.14216E-3, 6.40864E-3, 6.97750E-3, 7. << 981 { 5.0E-2, 3.444E-1}, 982 {0.15, 2.23096E-2, 2.30790E-2, 2.46978E-2, 2 << 982 { 4.0E-2, 3.323E-1}, 983 {0.2, 5.04022E-2, 5.18492E-2, 5.48682E-2, 5. << 983 { 3.0E-2, 3.144E-1}, 984 {0.3, 1.38018E-1, 1.41026E-1, 1.47244E-1, 1. << 984 { 2.0E-2, 2.854E-1}, 985 {0.4, 2.56001E-1, 2.60576E-1, 2.69992E-1, 2. << 985 { 1.5E-2, 2.629E-1}, 986 {0.5, 3.92181E-1, 3.98213E-1, 4.10597E-1, 4. << 986 987 {0.6, 5.37913E-1, 5.45268E-1, 5.60350E-1, 5. << 987 { 1.0E-2, 2.298E-1}, 988 {0.7, 6.87470E-1, 6.96021E-1, 7.13543E-1, 7. << 988 { 8.0E-3, 2.115E-1}, 989 {0.8, 8.37159E-1, 8.46790E-1, 8.66525E-1, 8. << 989 { 6.0E-3, 1.883E-1}, 990 {1.0, 1.12850, 1.14002, 1.16362, 1.18799, 1. << 990 { 5.0E-3, 1.741E-1}, 991 {1.2, 1.40232, 1.41545, 1.44232, 1.47007, 1. << 991 { 4.0E-3, 1.574E-1}, 992 {1.4, 1.65584, 1.67034, 1.70004, 1.73072, 1. << 992 { 3.0E-3, 1.372E-1}, 993 {1.5, 1.77496, 1.79009, 1.82107, 1.85307, 1. << 993 { 2.0E-3, 1.116E-1}, 994 {1.7, 1.99863, 2.01490, 2.04822, 2.08265, 2. << 994 { 1.5E-3, 9.559E-2}, 995 {2.0, 2.29325, 2.31100, 2.34738, 2.38501, 2. << 995 996 {3.0, 3.08887, 3.11036, 3.15445, 3.20013, 3. << 996 { 1.0E-3, 7.601E-2}, 997 {5.0, 4.07599, 4.10219, 4.15606, 4.21199, 4. << 997 { 8.0E-4, 6.668E-2}, 998 {7.0, 4.69647, 4.72577, 4.78608, 4.84877, 4. << 998 { 6.0E-4, 5.605E-2}, 999 {10.0, 5.32590, 5.35848, 5.42560, 5.49547, 5 << 999 { 5.0E-4, 5.008E-2}, 1000 }; << 1000 { 4.0E-4, 4.352E-2}, 1001 << 1001 { 3.0E-4, 3.617E-2}, 1002 const G4double bk2[29][11] = { << 1002 { 2.0E-4, 2.768E-2}, 1003 {0.005, 5.98040E-8, 7.25636E-8, 8.00602E-8, << 1003 { 1.5E-4, 2.279E-2}, 1004 {0.007, 2.95345E-7, 3.56497E-7, 3.92247E-7, << 1004 1005 {0.01, 1.55232E-6, 1.86011E-6, 2.03881E-6, << 1005 { 1.0E-4, 1.723E-2}, 1006 {0.015, 9.67802E-6, 1.14707E-5, 1.25008E-5, << 1006 { 8.0E-5, 1.473E-2}, 1007 {0.02, 3.38425E-5, 3.97259E-5, 4.30763E-5, << 1007 { 6.0E-5, 1.200E-2}, 1008 {0.03, 1.81920E-4, 2.10106E-4, 2.25918E-4, << 1008 { 5.0E-5, 1.052E-2}, 1009 {0.04, 5.59802E-4, 6.38103E-4, 6.81511E-4, << 1009 { 4.0E-5, 8.950E-3}, 1010 {0.05, 1.28002E-3, 1.44336E-3, 1.53305E-3, << 1010 { 3.0E-5, 7.246E-3}, 1011 {0.06, 2.42872E-3, 2.72510E-3, 2.88111E-3, << 1011 { 2.0E-5, 5.358E-3}, 1012 {0.08, 6.35222E-3, 6.99730E-3, 7.34446E-3, << 1012 { 1.5E-5, 4.313E-3}, 1013 {0.1, 1.26929E-2, 1.38803E-2, 1.44371E-2, 1 << 1013 { 0.0, 3.166E-3} 1014 {0.15, 3.96872E-2, 4.24699E-2, 4.39340E-2, << 1014 }; 1015 {0.2, 8.15290E-2, 8.62830E-2, 8.87650E-2, 9 << 1015 1016 {0.3, 1.99528E-1, 2.08471E-1, 2.13103E-1, 2 << 1016 for(i=0; i<104; i++) { 1017 {0.4, 3.47383E-1, 3.60369E-1, 3.67073E-1, 3 << 1017 ed[i] = nuca[i][0]; 1018 {0.5, 5.11214E-1, 5.27935E-1, 5.36554E-1, 5 << 1018 a[i] = nuca[i][1]; 1019 {0.6, 6.82122E-1, 7.02260E-1, 7.12631E-1, 7 << 1019 } 1020 {0.7, 8.54544E-1, 8.77814E-1, 8.89791E-1, 9 << 1020 1021 {0.8, 1.02508, 1.05121, 1.06466, 1.07838, 1 << 1021 // Constants 1022 {1.0, 1.35307, 1.38429, 1.40036, 1.41676, 1 << 1022 theZieglerFactor = eV*cm2*1.0e-15 ; 1023 {1.2, 1.65823, 1.69385, 1.71220, 1.73092, 1 << 1023 alpha2 = fine_structure_const*fine_structure_const; 1024 {1.4, 1.93902, 1.97852, 1.99887, 2.01964, 2 << 1024 lossFlucFlag = true; 1025 {1.5, 2.07055, 2.11182, 2.13309, 2.15480, 2 << 1025 1026 {1.7, 2.31700, 2.36154, 2.38451, 2.40798, 2 << 1026 // G.S. Khandelwal Nucl. Phys. A116(1968)97 - 111. 1027 {2.0, 2.64162, 2.69053, 2.71576, 2.74154, 2 << 1027 // "Shell corrections for K- and L- electrons 1028 {3.0, 3.51376, 3.57394, 3.60503, 3.63684, 3 << 1028 nK = 20; 1029 {5.0, 4.59935, 4.67433, 4.71316, 4.75293, 4 << 1029 nL = 26; 1030 {7.0, 5.28542, 5.37040, 5.41445, 5.45962, 5 << 1030 nEtaK = 29; 1031 {10.0, 5.98474, 6.08046, 6.13015, 6.18112, << 1031 nEtaL = 28; 1032 }; << 1032 1033 << 1033 const G4double d[11] = {0., 0., 0., 1.72, 2.09, 2.48, 2.82, 3.16, 3.53, 3.84, 4.15}; 1034 const G4double bls1[28][10] = { << 1034 const G4double thek[20] = {0.64, 0.65, 0.66, 0.68, 0.70, 0.72, 0.74, 0.75, 0.76, 0.78, 1035 {0.005, 2.4111E-4, 2.5612E-4, 2.7202E-4, 3. << 1035 0.80, 0.82, 0.84, 0.85, 0.86, 0.88, 0.90, 0.92, 0.94, 0.95}; 1036 {0.007, 6.3947E-4, 6.7058E-4, 7.0295E-4, 7. << 1036 const G4double sk[20] = {1.9477, 1.9232, 1.8996, 1.8550, 1.8137, 1037 {0.01, 1.5469E-3, 1.6036E-3, 1.6622E-3, 1.7 << 1037 1.7754, 1.7396, 1.7223, 1.7063, 1.6752, 1038 {0.015, 3.7221E-3, 3.8404E-3, 3.9650E-3, 4. << 1038 1.6461, 1.6189, 1.5933, 1.5811, 1.5693, 1039 {0.02, 6.9449E-3, 7.1910E-3, 7.4535E-3, 8.0 << 1039 1.5467, 1.5254, 1.5053, 1.4863, 1.4772}; 1040 {0.03, 1.7411E-2, 1.8180E-2, 1.8997E-2, 2.0 << 1040 const G4double tk[20] = {2.5222, 2.5125, 2.5026, 2.4821, 2.4608, 1041 {0.04, 3.8474E-2, 4.0056E-2, 4.1718E-2, 4.5 << 1041 2.4388, 2.4163, 2.4044, 2.3933, 2.3701, 1042 {0.05, 6.7371E-2, 6.9928E-2, 7.2596E-2, 7.8 << 1042 2.3466, 2.3229, 2.2992, 2.2872, 2.2753, 1043 {0.06, 1.0418E-1, 1.0778E-1, 1.1152E-1, 1.1 << 1043 2.2515, 2.2277, 2.2040, 2.1804, 2.1686}; 1044 {0.08, 1.9647E-1, 2.0217E-1, 2.0805E-1, 2.2 << 1044 const G4double uk[20] = {1.9999, 2.0134, 2.0258, 2.0478, 2.0662, 1045 {0.1, 3.0594E-1, 3.1361E-1, 3.2150E-1, 3.37 << 1045 2.0817, 2.0945, 2.0999, 2.1049, 2.1132, 1046 {0.15, 6.1411E-1, 6.2597E-1, 6.3811E-1, 6.6 << 1046 2.1197, 2.1246, 2.1280, 2.1292, 2.1301, 1047 {0.2, 9.3100E-1, 9.5614E-1, 9.7162E-1, 1.00 << 1047 2.1310, 2.1310, 2.1300, 2.1283, 2.1271}; 1048 {0.3, 1.5172, 1.5372, 1.5576, 1.5998, 1.643 << 1048 const G4double vk[20] = {8.3410, 8.3373, 8.3340, 8.3287, 8.3247, 1049 {0.4, 2.0173, 2.0408, 2.0647, 2.1142, 2.165 << 1049 8.3219, 8.3201, 8.3194, 8.3191, 8.3188, 1050 {0.5, 2.3932, 2.4193, 2.4460, 2.5011, 2.558 << 1050 8.3191, 8.3199, 8.3211, 8.3218, 8.3226, 1051 {0.6, 2.7091, 2.7374, 2.7663, 2.8260, 2.888 << 1051 8.3244, 8.3264, 8.3285, 8.3308, 8.3320}; 1052 {0.7, 2.9742, 3.0044, 3.0352, 3.0988, 3.165 << 1052 1053 {0.8, 3.2222, 3.2539, 3.2863, 3.3532, 3.423 << 1053 for(i=0; i<11; i++) { ZD[i] = d[i];} 1054 {1.0, 3.6690, 3.7033, 3.7384, 3.8108, 3.886 << 1054 1055 {1.2, 3.9819, 4.0183, 4.0555, 4.1324, 4.213 << 1055 for(i=0; i<nK; i++) { 1056 {1.4, 4.2745, 4.3127, 4.3517, 4.4324, 4.517 << 1056 TheK[i] = thek[i]; 1057 {1.5, 4.4047, 4.4436, 4.4834, 4.5658, 4.652 << 1057 SK[i] = sk[i]; 1058 {1.7, 4.6383, 4.6787, 4.7200, 4.8054, 4.894 << 1058 TK[i] = tk[i]; 1059 {2.0, 4.9369, 4.9791, 5.0223, 5.1116, 5.205 << 1059 UK[i] = uk[i]; 1060 {3.0, 5.6514, 5.6981, 5.7459, 5.8450, 5.948 << 1060 VK[i] = vk[i]; 1061 {5.0, 6.4665, 6.5189, 6.5724, 6.6835, 6.800 << 1061 } 1062 {7.0, 6.8634, 6.9194, 6.9767, 7.0957, 7.220 << 1062 1063 }; << 1063 const G4double thel[26] = {0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.35, 0.36, 0.38, 0.40, 1064 << 1064 0.42, 0.44, 0.45, 0.46, 0.48, 0.50, 0.52, 0.54, 0.55, 0.56, 1065 const G4double bls2[28][10] = { << 1065 0.58, 0.60, 0.62, 0.64, 0.65, 0.66}; 1066 {0.005, 5.4561E-4, 6.0905E-4, 6.7863E-4, 7. << 1066 const G4double sl[26] = {15.3343, 13.9389, 12.7909, 11.8343, 11.0283, 1067 {0.007, 1.2068E-3, 1.3170E-3, 1.4377E-3, 1. << 1067 10.3424, 10.0371, 9.7537, 9.2443, 8.8005, 1068 {0.01, 2.6832E-3, 2.9017E-3, 3.1534E-3, 3.4 << 1068 8.4114, 8.0683, 7.9117, 7.7641, 7.4931, 1069 {0.015, 6.2775E-3, 6.9077E-3, 7.6525E-3, 8. << 1069 7.2506, 7.0327, 6.8362, 6.7452, 6.6584, 1070 {0.02, 1.2561E-2, 1.3943E-2, 1.5553E-2, 1.7 << 1070 6.4969, 6.3498, 6.2154, 6.0923, 6.0345, 5.9792}; 1071 {0.03, 3.3750E-2, 3.7470E-2, 4.1528E-2, 4.6 << 1071 const G4double tl[26] = {35.0669, 33.4344, 32.0073, 30.7466, 29.6226, 1072 {0.04, 6.9619E-2, 7.6098E-2, 8.3249E-2, 9.1 << 1072 28.6128, 28.4149, 27.6991, 26.8674, 26.1061, 1073 {0.05, 1.1522E-1, 1.2470E-1, 1.3504E-1, 1.4 << 1073 25.4058, 24.7587, 24.4531, 24.1583, 23.5992, 1074 {0.06, 1.6931E-1, 1.8179E-1, 1.9530E-1, 2.0 << 1074 23.0771, 22.5880, 22.1285, 21.9090, 21.6958, 1075 {0.08, 2.9540E-1, 3.1361E-1, 3.3312E-1, 3.5 << 1075 21.2872, 20.9006, 20.5341, 20.1859, 20.0183, 19.8546}; 1076 {0.1, 4.3613E-1, 4.5956E-1, 4.852E-1, 5.1 << 1076 const G4double ul[26] = {0.1215, 0.5265, 0.8411, 1.0878, 1.2828, 1077 {0.15, 8.1014E-1, 8.4453E-1, 8.8093E-1, 9.1 << 1077 1.4379, 1.5032, 1.5617, 1.6608, 1.7401, 1078 {0.2, 1.1888, 1.2319, 1.2774, 1.3255, 1.350 << 1078 1.8036, 1.8543, 1.8756, 1.8945, 1.9262, 1079 {0.3, 1.8422, 1.8983, 1.9575, 2.0201, 2.052 << 1079 1.9508, 1.9696, 1.9836, 1.9890, 1.9935, 1080 {0.4, 2.3984, 2.4642, 2.5336, 2.6070, 2.645 << 1080 2.0001, 2.0039, 2.0053, 2.0049, 2.0040, 2.0028}; 1081 {0.5, 2.8181, 2.8915, 2.9690, 3.0509, 3.093 << 1081 for(i=0; i<nL; i++) { 1082 {0.6, 3.1698, 3.2494, 3.3336, 3.4226, 3.469 << 1082 TheL[i] = thel[i]; 1083 {0.7, 3.4652, 3.5502, 3.6400, 3.7351, 3.784 << 1083 SL[i] = sl[i]; 1084 {0.8, 3.7392, 3.8289, 3.9236, 4.0239, 4.076 << 1084 TL[i] = tl[i]; 1085 {1.0, 4.2295, 4.3269, 4.4299, 4.5391, 4.596 << 1085 UL[i] = ul[i]; 1086 {1.2, 4.5777, 4.6814, 4.7912, 4.9076, 4.968 << 1086 } 1087 {1.4, 4.9001, 5.0092, 5.1247, 5.2473, 5.311 << 1087 1088 {1.5, 5.0434, 5.1550, 5.2731, 5.3984, 5.464 << 1088 const G4double eta[29] = {0.005, 0.007, 0.01, 0.015, 0.02, 1089 {1.7, 5.3011, 5.4170, 5.5398, 5.6701, 5.737 << 1089 0.03, 0.04, 0.05, 0.06, 0.08, 1090 {2.0, 5.6308, 5.7523, 5.8811, 6.0174, 6.089 << 1090 0.1, 0.15, 0.2, 0.3, 0.4, 1091 {3.0, 6.4224, 6.5580, 6.7019, 6.8549, 6.935 << 1091 0.5, 0.6, 0.7, 0.8, 1.0, 1092 {5.0, 7.3338, 7.4872, 7.6500, 7.8235, 7.914 << 1092 1.2, 1.4, 1.5, 1.7, 2.0, 3.0, 5.0, 7.0, 10.}; 1093 {7.0, 7.7938, 7.9588, 8.1342, 8.3211, 8.419 << 1093 1094 }; << 1094 const G4double bk1[29][11] = { 1095 << 1095 {0.005, 1.34782E-8, 1.46132E-8, 1.72179E-8, 2.03521E-8, 2.41370E-8, 2.87247E-8, 3.13778E-8, 3.43072E-8, 4.11274E-8, 4.94946E-8}, 1096 const G4double bls3[28][9] = { << 1096 {0.007, 6.87555E-8, 7.44373E-8, 8.74397E-8, 1.03022E-7, 1.21760E-7, 1.44370E-7, 1.57398E-7, 1.71747E-7, 2.05023E-7, 2.45620E-7}, 1097 {0.005, 1.2895E-3, 1.3670E-3, 1.4524E-3, 1. << 1097 {0.01, 3.78413E-7, 4.08831E-7, 4.78154E-7, 5.60760E-7, 6.59478E-7, 7.77847E-7, 8.45709E-7, 9.20187E-7, 1.09192E-6, 1.29981E-6}, 1098 {0.007, 2.6467E-3, 2.8242E-3, 3.0238E-3, 3. << 1098 {0.015, 2.53200E-6, 2.72664E-6, 3.16738E-6, 3.68791E-6, 4.30423E-6, 5.03578E-6, 5.45200E-6, 5.90633E-6, 6.94501E-6, 8.18757E-6}, 1099 {0.01, 6.1472E-3, 6.6086E-3, 7.1246E-3, 8.3 << 1099 {0.02, 9.43891E-6, 1.01339E-5, 1.16984E-5, 1.35313E-5, 1.56829E-5, 1.82138E-5, 1.96439E-5, 2.11973E-5, 2.47216E-5, 2.88935E-5}, 1100 {0.015, 1.63316E-2, 1.7572E-2, 1.8932E-2, 2 << 1100 {0.03, 5.67088E-5, 6.05576E-5, 6.91311E-5, 7.90324E-5, 9.04832E-5, 1.03744E-4, 1.11147E-4, 1.19122E-4, 1.36980E-4, 1.57744E-4}, 1101 {0.02, 3.2634E-2, 3.4900E-2, 3.7348E-2, 4.2 << 1101 {0.04, 1.91576E-4, 2.03626E-4, 2.30230E-4, 2.60584E-4, 2.95248E-4, 3.34870E-4, 3.56771E-4, 3.80200E-4, 4.32104E-4, 4.91584E-4}, 1102 {0.03, 7.9907E-2, 8.4544E-2, 8.9486E-2, 1.0 << 1102 {0.05, 4.74226E-4, 5.02006E-4, 5.62872E-4, 6.31597E-4, 7.09244E-4, 7.97023E-4, 8.45134E-4, 8.96410E-4, 1.00867E-3, 1.13590E-3}, 1103 {0.04, 1.4523E-1, 1.5237E-1, 1.5985E-1, 1.7 << 1103 {0.06, 9.67285E-4, 1.02029E-3, 1.13566E-3, 1.26476E-3, 1.46928E-3, 1.57113E-3, 1.65921E-3, 1.75244E-3, 1.95562E-3, 2.18336E-3}, 1104 {0.05, 2.2082E-1, 2.3036E-1, 2.4038E-1, 2.6 << 1104 {0.08, 2.81537E-3, 2.95200E-3, 3.24599E-3, 3.57027E-3, 3.92793E-3, 4.32246E-3, 4.53473E-3, 4.75768E-3, 5.23785E-3, 5.76765E-3}, 1105 {0.06, 3.0423E-1, 3.1611E-1, 3.2854E-1, 3.5 << 1105 {0.1, 6.14216E-3, 6.40864E-3, 6.97750E-3, 7.59781E-3, 8.27424E-3, 9.01187E-3, 9.40534E-3, 9.81623E-3, 1.06934E-2, 1.16498E-2}, 1106 {0.08, 4.8536E-1, 5.0156E-1, 5.1846E-1, 5.5 << 1106 {0.15, 2.23096E-2, 2.30790E-2, 2.46978E-2, 2.64300E-2, 2.82832E-2, 3.02661E-2, 3.13090E-2, 3.23878E-2, 3.46580E-2, 3.70873E-2}, 1107 {0.1, 6.7586E-1, 6.9596E-1, 7.1688E-1, 7.61 << 1107 {0.2, 5.04022E-2, 5.18492E-2, 5.48682E-2, 5.80617E-2, 6.14403E-2, 6.50152E-2, 6.68798E-2, 6.87981E-2, 7.28020E-2, 7.70407E-2}, 1108 {0.15, 1.1544, 1.1828, 1.2122, 1.2746, 1.34 << 1108 {0.3, 1.38018E-1, 1.41026E-1, 1.47244E-1, 1.53743E-1, 1.60536E-1, 1.67641E-1, 1.71315E-1, 1.75074E-1, 1.82852E-1, 1.90997E-1}, 1109 {0.2, 1.6167, 1.6517, 1.6880, 1.7650, 1.848 << 1109 {0.4, 2.56001E-1, 2.60576E-1, 2.69992E-1, 2.79776E-1, 2.89946E-1, 3.00525E-1, 3.05974E-1, 3.11533E-1, 3.22994E-1, 3.34935E-1}, 1110 {0.3, 2.3979, 2.4432, 2.4902, 2.5899, 2.698 << 1110 {0.5, 3.92181E-1, 3.98213E-1, 4.10597E-1, 4.23427E-1, 4.36726E-1, 4.50519E-1, 4.57610E-1, 4.64834E-1, 4.79700E-1, 4.95148E-1}, 1111 {0.4, 3.0502, 3.1034, 3.1586, 3.2758, 3.403 << 1111 {0.6, 5.37913E-1, 5.45268E-1, 5.60350E-1, 5.75948E-1, 5.92092E-1, 6.08811E-1, 6.17396E-1, 6.26136E-1, 6.44104E-1, 6.62752E-1}, 1112 {0.5, 3.5466, 3.6062, 3.6681, 3.7994, 3.942 << 1112 {0.7, 6.87470E-1, 6.96021E-1, 7.13543E-1, 7.31650E-1, 7.50373E-1, 7.69748E-1, 7.79591E-1, 7.89811E-1, 8.10602E-1, 8.32167E-1}, 1113 {0.6, 3.9620, 4.0270, 4.0945, 4.2378, 4.393 << 1113 {0.8, 8.37159E-1, 8.46790E-1, 8.66525E-1, 8.86910E-1, 9.07979E-1, 9.29772E-1, 9.40953E-1, 9.52331E-1, 9.75701E-1, 9.99934E-1}, 1114 {0.7, 4.3020, 4.3715, 4.4438, 4.5974, 4.764 << 1114 {1.0, 1.12850, 1.14002, 1.16362, 1.18799, 1.21317, 1.23921, 1.25257, 1.26616, 1.29408, 1.32303}, 1115 {0.8, 4.6336, 4.7072, 4.7837, 4.9463, 5.123 << 1115 {1.2, 1.40232, 1.41545, 1.44232, 1.47007, 1.49875, 1.52842, 1.54364, 1.55913, 1.59095, 1.62396}, 1116 {1.0, 5.2041, 5.2845, 5.3682, 5.5462, 5.740 << 1116 {1.4, 1.65584, 1.67034, 1.70004, 1.73072, 1.76244, 1.79526, 1.81210, 1.82925, 1.86448, 1.90104}, 1117 {1.2, 5.6182, 5.7044, 5.7940, 5.9848, 6.192 << 1117 {1.5, 1.77496, 1.79009, 1.82107, 1.85307, 1.88617, 1.92042, 1.93800, 1.95590, 1.99269, 2.03087}, 1118 {1.4, 5.9967, 6.0876, 6.1823, 6.3839, 6.603 << 1118 {1.7, 1.99863, 2.01490, 2.04822, 2.08265, 2.11827, 2.15555, 2.17409, 2.19337, 2.23302, 2.27419}, 1119 {1.5, 6.1652, 6.2583, 6.3553, 6.5618, 6.787 << 1119 {2.0, 2.29325, 2.31100, 2.34738, 2.38501, 2.42395, 2.46429, 2.48401, 2.50612, 2.54955, 2.59468}, 1120 {1.7, 6.4686, 6.5657, 6.6668, 6.8823, 7.117 << 1120 {3.0, 3.08887, 3.11036, 3.15445, 3.20013, 3.24748, 3.29664, 3.32192, 3.34770, 3.40081, 3.45611}, 1121 {2.0, 6.8577, 6.9600, 7.0666, 7.2937, 7.541 << 1121 {5.0, 4.07599, 4.10219, 4.15606, 4.21199, 4.27010, 4.33056, 4.36172, 4.39353, 4.45918, 4.52772}, 1122 {3.0, 7.7981, 7.9134, 8.0336, 8.2901, 8.570 << 1122 {7.0, 4.69647, 4.72577, 4.78608, 4.84877, 4.91402, 4.98200, 5.01707, 5.05290, 5.12695, 5.20436}, 1123 {5.0, 8.8978, 9.0297, 9.1673, 9.4612, 9.783 << 1123 {10.0, 5.32590, 5.35848, 5.42560, 5.49547, 5.56830, 5.64429, 5.68353, 5.72366, 5.80666, 5.89359} 1124 {7.0, 9.4819, 9.6248, 9.7739, 10.0926, 10.4 << 1124 }; 1125 }; << 1125 1126 << 1126 const G4double bk2[29][11] = { 1127 const G4double bll1[28][10] = { << 1127 {0.005, 5.98040E-8, 7.25636E-8, 8.00602E-8, 8.84294E-8, 1.08253E-7, 1.33148E-7, 1.64573E-7, 2.04459E-7, 2.28346E-7, 2.55370E-7}, 1128 {0.005, 3.6324E-5, 4.0609E-5, 4.5430E-5, 5. << 1128 {0.007, 2.95345E-7, 3.56497E-7, 3.92247E-7, 4.32017E-7, 5.25688E-7, 6.42391E-7, 7.88464E-7, 9.72171E-7, 1.08140E-6, 1.20435E-6}, 1129 {0.007, 1.8110E-4, 2.0001E-4, 2.2099E-4, 2. << 1129 {0.01, 1.55232E-6, 1.86011E-6, 2.03881E-6, 2.23662E-6, 2.69889E-6, 3.26860E-6, 3.26860E-6, 4.84882E-6, 5.36428E-6, 5.94048E-6}, 1130 {0.01, 8.6524E-4, 9.4223E-4, 1.0262E-3, 1.2 << 1130 {0.015, 9.67802E-6, 1.14707E-5, 1.25008E-5, 1.36329E-5, 1.62480E-5, 1.94200E-5, 2.32783E-5, 2.79850E-5, 3.07181E-5, 3.37432E-5}, 1131 {0.015, 4.2293E-3, 4.5314E-3, 4.8551E-3, 5. << 1131 {0.02, 3.38425E-5, 3.97259E-5, 4.30763E-5, 4.67351E-5, 5.51033E-5, 6.51154E-5, 7.71154E-5, 9.15431E-5, 9.98212E-5, 1.08909E-4}, 1132 {0.02, 1.1485E-2, 1.2172E-2, 1.2900E-2, 1.4 << 1132 {0.03, 1.81920E-4, 2.10106E-4, 2.25918E-4, 2.43007E-4, 2.81460E-4, 3.26458E-4, 3.79175E-4, 4.41006E-4, 4.75845E-4, 5.13606E-4}, 1133 {0.03, 3.9471E-2, 4.1270E-2, 4.3149E-2, 4.7 << 1133 {0.04, 5.59802E-4, 6.38103E-4, 6.81511E-4, 7.28042E-4, 8.31425E-4, 9.50341E-4, 1.08721E-3, 1.24485E-3, 1.33245E-3, 1.42650E-3}, 1134 {0.04, 8.4454E-2, 8.7599E-2, 9.0860E-2, 9.7 << 1134 {0.05, 1.28002E-3, 1.44336E-3, 1.53305E-3, 1.62855E-3, 1.83861E-3, 2.07396E-3, 2.34750E-3, 2.65469E-3, 2.82358E-3, 3.00358E-3}, 1135 {0.05, 1.4339E-1, 1.4795E-1, 1.5266E-1, 1.6 << 1135 {0.06, 2.42872E-3, 2.72510E-3, 2.88111E-3, 3.04636E-3, 3.40681E-3, 3.81132E-3, 4.26536E-3, 4.77507E-3, 5.05294E-3, 5.34739E-3}, 1136 {0.06, 2.1304E-1, 2.1899E-1, 2.2512E-1, 2.3 << 1136 {0.08, 6.35222E-3, 6.99730E-3, 7.34446E-3, 7.70916E-3, 8.49478E-3, 9.36187E-3, 1.03189E-2, 1.13754E-2, 1.19441E-2, 1.25417E-2}, 1137 {0.08, 3.7382E-1, 3.8241E-1, 3.9122E-1, 4.0 << 1137 {0.1, 1.26929E-2, 1.38803E-2, 1.44371E-2, 1.50707E-2, 1.64235E-2, 1.78989E-2, 1.95083E-2, 2.12640E-2, 2.22009E-2, 2.31795E-2}, 1138 {0.1, 5.5056E-1, 5.6151E-1, 5.7273E-1, 5.96 << 1138 {0.15, 3.96872E-2, 4.24699E-2, 4.39340E-2, 4.54488E-2, 4.86383E-2, 5.20542E-2, 5.57135E-2, 5.96350E-2, 6.17003E-2, 6.38389E-2}, 1139 {0.15, 1.0066, 1.0224, 1.0386, 1.0721, 1.10 << 1139 {0.2, 8.15290E-2, 8.62830E-2, 8.87650E-2, 9.13200E-2, 9.66589E-2, 1.02320E-1, 1.08326E-1, 1.14701E-1, 1.18035E-1, 1.21472E-1}, 1140 {0.2, 1.4376, 1.4572, 1.4773, 1.5188, 1.562 << 1140 {0.3, 1.99528E-1, 2.08471E-1, 2.13103E-1, 2.17848E-1, 2.27689E-1, 2.38022E-1, 2.48882E-1, 2.60304E-1, 2.66239E-1, 2.72329E-1}, 1141 {0.3, 2.1712, 2.1964, 2.2223, 2.2758, 2.332 << 1141 {0.4, 3.47383E-1, 3.60369E-1, 3.67073E-1, 3.73925E-1, 3.88089E-1, 4.02900E-1, 4.18404E-1, 4.34647E-1, 4.43063E-1, 4.51685E-1}, 1142 {0.4, 2.7500, 2.7793, 2.8094, 2.8719, 2.937 << 1142 {0.5, 5.11214E-1, 5.27935E-1, 5.36554E-1, 5.45354E-1, 5.63515E-1, 5.82470E-1, 6.02273E-1, 6.22986E-1, 6.33705E-1, 6.44677E-1}, 1143 {0.5, 3.2033, 3.2359, 3.2693, 3.3389, 3.412 << 1143 {0.6, 6.82122E-1, 7.02260E-1, 7.12631E-1, 7.23214E-1, 7.45041E-1, 7.67800E-1, 7.91559E-1, 8.16391E-1, 8.29235E-1, 8.42380E-1}, 1144 {0.6, 3.6038, 3.6391, 3.6753, 3.7506, 3.830 << 1144 {0.7, 8.54544E-1, 8.77814E-1, 8.89791E-1, 9.02008E-1, 9.27198E-1, 9.53454E-1, 9.80856E-1, 1.00949, 1.02430, 1.03945}, 1145 {0.7, 3.9106, 3.9482, 3.9867, 4.0670, 4.152 << 1145 {0.8, 1.02508, 1.05121, 1.06466, 1.07838, 1.10667, 1.13615, 1.16692, 1.19907, 1.21570, 1.23272}, 1146 {0.8, 4.1790, 4.2185, 4.2590, 4.3437, 4.433 << 1146 {1.0, 1.35307, 1.38429, 1.40036, 1.41676, 1.45057, 1.48582, 1.52263, 1.56111, 1.58102, 1.60140}, 1147 {1.0, 4.6344, 4.6772, 4.7212, 4.8131, 4.910 << 1147 {1.2, 1.65823, 1.69385, 1.71220, 1.73092, 1.76954, 1.80983, 1.85192, 1.89596, 1.91876, 1.94211}, 1148 {1.2, 4.9787, 5.0242, 5.0711, 5.1689, 5.272 << 1148 {1.4, 1.93902, 1.97852, 1.99887, 2.01964, 2.06251, 2.10727, 2.15406, 2.20304, 2.22842, 2.25442}, 1149 {1.4, 5.2688, 5.3166, 5.3658, 5.4687, 5.578 << 1149 {1.5, 2.07055, 2.11182, 2.13309, 2.15480, 2.19963, 2.24644, 2.29539, 2.34666, 2.37323, 2.40045}, 1150 {1.5, 5.3966, 5.4454, 5.4957, 5.6009, 5.712 << 1150 {1.7, 2.31700, 2.36154, 2.38451, 2.40798, 2.45641, 2.50703, 2.56000, 2.61552, 2.64430, 2.67381}, 1151 {1.7, 5.6255, 5.6762, 5.7284, 5.8377, 5.954 << 1151 {2.0, 2.64162, 2.69053, 2.71576, 2.74154, 2.79481, 2.85052, 2.90887, 2.97009, 3.00185, 3.03442}, 1152 {2.0, 5.9170, 5.9701, 6.0248, 6.1395, 6.261 << 1152 {3.0, 3.51376, 3.57394, 3.60503, 3.63684, 3.70268, 3.77170, 3.84418, 3.92040, 3.96003, 4.00073}, 1153 {3.0, 6.6210, 6.6801, 6.7411, 6.8692, 7.006 << 1153 {5.0, 4.59935, 4.67433, 4.71316, 4.75293, 4.83543, 4.92217, 5.01353, 5.10992, 5.16014, 5.21181}, 1154 {5.0, 7.4620, 7.5288, 7.5977, 7.7428, 7.898 << 1154 {7.0, 5.28542, 5.37040, 5.41445, 5.45962, 5.55344, 5.65226, 5.79496, 5.90517, 5.96269, 6.02191}, 1155 {7.0, 7.7362, 7.8079, 7.8821, 8.0383, 8.206 << 1155 {10.0, 5.98474, 6.08046, 6.13015, 6.18112, 6.28715, 6.39903, 6.51728, 6.64249, 6.70792, 6.77535} 1156 }; << 1156 }; 1157 << 1157 1158 const G4double bll2[28][10] = { << 1158 const G4double bls1[28][10] = { 1159 {0.005, 1.8339E-4, 2.3330E-4, 2.9738E-4, 3. << 1159 {0.005, 2.4111E-4, 2.5612E-4, 2.7202E-4, 3.0658E-4, 3.4511E-4, 3.8795E-4, 4.3542E-4, 4.6100E-4, 4.8786E-4}, 1160 {0.007, 7.5042E-4, 9.2355E-4, 1.1375E-3, 1. << 1160 {0.007, 6.3947E-4, 6.7058E-4, 7.0295E-4, 7.7167E-4, 8.4592E-4, 9.2605E-4, 1.0125E-3, 1.0583E-3, 1.1058E-3}, 1161 {0.01, 2.8829E-3, 3.4275E-3, 4.0758E-3, 4.8 << 1161 {0.01, 1.5469E-3, 1.6036E-3, 1.6622E-3, 1.7856E-3, 1.9181E-3, 2.1615E-3, 2.3178E-3, 2.4019E-3, 2.4904E-3}, 1162 {0.015, 1.1087E-2, 1.2716E-2, 1.4581E-2, 1. << 1162 {0.015, 3.7221E-3, 3.8404E-3, 3.9650E-3, 4.2354E-3, 4.5396E-3, 4.8853E-3, 5.2820E-3, 5.5031E-3, 5.7414E-3}, 1163 {0.02, 2.5786E-2, 2.8922E-2, 3.2435E-2, 3.6 << 1163 {0.02, 6.9449E-3, 7.1910E-3, 7.4535E-3, 8.0336E-3, 8.6984E-3, 9.4638E-3, 1.0348E-2, 1.0841E-2, 1.1372E-2}, 1164 {0.03, 7.3461E-2, 8.0264E-2, 8.7705E-2, 9.5 << 1164 {0.03, 1.7411E-2, 1.8180E-2, 1.8997E-2, 2.0784E-2, 2.2797E-2, 2.5066E-2, 2.7622E-2, 2.9020E-2, 3.0503E-2}, 1165 {0.04, 1.4094E-1, 1.5172E-1, 1.6336E-1, 1.7 << 1165 {0.04, 3.8474E-2, 4.0056E-2, 4.1718E-2, 4.5300E-2, 4.9254E-2, 5.3619E-2, 5.8436E-2, 6.1028E-2, 6.3752E-2}, 1166 {0.05, 2.2289E-1, 2.3762E-1, 2.5344E-1, 2.7 << 1166 {0.05, 6.7371E-2, 6.9928E-2, 7.2596E-2, 7.8282E-2, 8.4470E-2, 9.1206E-2, 9.8538E-2, 1.0244E-1, 1.0652E-1}, 1167 {0.06, 3.1503E-1, 3.3361E-1, 3.5346E-1, 3.7 << 1167 {0.06, 1.0418E-1, 1.0778E-1, 1.1152E-1, 1.1943E-1, 1.2796E-1, 1.3715E-1, 1.4706E-1, 1.5231E-1, 1.5776E-1}, 1168 {0.08, 5.1793E-1, 5.4368E-1, 5.7109E-1, 6.0 << 1168 {0.08, 1.9647E-1, 2.0217E-1, 2.0805E-1, 2.2038E-1, 2.3351E-1, 2.4751E-1, 2.6244E-1, 2.7027E-1, 2.7837E-1}, 1169 {0.1, 7.3258E-1, 7.6481E-1, 7.9907E-1, 8.35 << 1169 {0.1, 3.0594E-1, 3.1361E-1, 3.2150E-1, 3.3796E-1, 3.5537E-1, 3.7381E-1, 3.9336E-1, 4.0357E-1, 4.1410E-1}, 1170 {0.15, 1.2677, 1.3137, 1.3626, 1.4147, 1.44 << 1170 {0.15, 6.1411E-1, 6.2597E-1, 6.3811E-1, 6.6330E-1, 6.8974E-1, 7.1753E-1, 7.4678E-1, 7.6199E-1, 7.7761E-1}, 1171 {0.2, 1.7615, 1.8188, 1.8797, 1.9446, 1.978 << 1171 {0.2, 9.3100E-1, 9.5614E-1, 9.7162E-1, 1.0037, 1.0372, 1.0723, 1.1092, 1.1284, 1.1480}, 1172 {0.3, 2.5909, 2.6658, 2.7457, 2.8312, 2.876 << 1172 {0.3, 1.5172, 1.5372, 1.5576, 1.5998, 1.6438, 1.6899, 1.7382, 1.7632, 1.7889}, 1173 {0.4, 3.2417, 3.3302, 3.4249, 3.5266, 3.580 << 1173 {0.4, 2.0173, 2.0408, 2.0647, 2.1142, 2.1659, 2.2199, 2.2765, 2.3059, 2.3360}, 1174 {0.5, 3.7527, 3.8523, 3.9591, 4.0741, 4.135 << 1174 {0.5, 2.3932, 2.4193, 2.4460, 2.5011, 2.5587, 2.6190, 2.6821, 2.7148, 2.7484}, 1175 {0.6, 4.2013, 4.3103, 4.4274, 4.5537, 4.620 << 1175 {0.6, 2.7091, 2.7374, 2.7663, 2.8260, 2.8884, 2.9538, 3.0222, 3.0577, 3.0941}, 1176 {0.7, 4.5493, 4.6664, 4.7925, 4.9286, 5.000 << 1176 {0.7, 2.9742, 3.0044, 3.0352, 3.0988, 3.1652, 3.2349, 3.3079, 3.3457, 3.3845}, 1177 {0.8, 4.8537, 4.9780, 5.1119, 5.2568, 5.333 << 1177 {0.8, 3.2222, 3.2539, 3.2863, 3.3532, 3.4232, 3.4965, 3.5734, 3.6133, 3.6542}, 1178 {1.0, 5.3701, 5.5066, 5.6540, 5.8138, 5.898 << 1178 {1.0, 3.6690, 3.7033, 3.7384, 3.8108, 3.8866, 3.9661, 4.0495, 4.0928, 4.1371}, 1179 {1.2, 5.7648, 5.9114, 6.0701, 6.2424, 6.334 << 1179 {1.2, 3.9819, 4.0183, 4.0555, 4.1324, 4.2130, 4.2974, 4.3861, 4.4321, 4.4794}, 1180 {1.4, 6.0976, 6.2530, 6.4214, 6.6044, 6.702 << 1180 {1.4, 4.2745, 4.3127, 4.3517, 4.4324, 4.5170, 4.6056, 4.6988, 4.7471, 4.7968}, 1181 {1.5, 6.2447, 6.4041, 6.5768, 6.7647, 6.865 << 1181 {1.5, 4.4047, 4.4436, 4.4834, 4.5658, 4.6521, 4.7426, 4.8378, 4.8872, 4.9379}, 1182 {1.7, 6.5087, 6.6752, 6.8558, 7.0526, 7.157 << 1182 {1.7, 4.6383, 4.6787, 4.7200, 4.8054, 4.8949, 4.9888, 5.0876, 5.1388, 5.1915}, 1183 {2.0, 6.8458, 7.0218, 7.2129, 7.4213, 7.532 << 1183 {2.0, 4.9369, 4.9791, 5.0223, 5.1116, 5.2053, 5.3036, 5.4070, 5.4607, 5.5159}, 1184 {3.0, 7.6647, 7.8644, 8.0819, 8.3189, 8.447 << 1184 {3.0, 5.6514, 5.6981, 5.7459, 5.8450, 5.9489, 6.0581, 6.1730, 6.2328, 6.2943}, 1185 {5.0, 8.6515, 8.8816, 9.1330, 9.4090, 9.557 << 1185 {5.0, 6.4665, 6.5189, 6.5724, 6.6835, 6.8003, 6.9231, 7.0525, 7.1199, 7.1892}, 1186 {7.0, 9.0221, 9.2724, 9.5464, 9.8477, 10.00 << 1186 {7.0, 6.8634, 6.9194, 6.9767, 7.0957, 7.2208, 7.3526, 7.4915, 7.5639, 7.6384} 1187 }; << 1187 }; 1188 << 1188 1189 const G4double bll3[28][9] = { << 1189 1190 {0.005, 1.3190E-3, 1.4961E-3, 1.6974E-3, 2. << 1190 const G4double bls2[28][10] = { 1191 {0.007, 4.0158E-3, 4.4623E-3, 4.9592E-3, 6. << 1191 {0.005, 5.4561E-4, 6.0905E-4, 6.7863E-4, 7.5494E-4, 7.9587E-4, 8.3883E-4, 9.3160E-4, 1.0352E-3, 1.1529E-3}, 1192 {0.01, 1.1509E-2, 1.2548E-2, 1.3681E-2, 1.6 << 1192 {0.007, 1.2068E-3, 1.3170E-3, 1.4377E-3, 1.5719E-3, 1.6451E-3, 1.7231E-3, 1.8969E-3, 2.1009E-3, 2.3459E-3}, 1193 {0.015, 3.3070E-2, 3.5408E-2, 3.7914E-2, 4. << 1193 {0.01, 2.6832E-3, 2.9017E-3, 3.1534E-3, 3.4479E-3, 3.6149E-3, 3.7976E-3, 4.2187E-3, 4.7320E-3, 5.3636E-3}, 1194 {0.02, 6.4555E-2, 6.8394E-2, 7.2472E-2, 8.1 << 1194 {0.015, 6.2775E-3, 6.9077E-3, 7.6525E-3, 8.5370E-2, 9.0407E-3, 9.5910E-3, 1.0850E-2, 1.2358E-2, 1.4165E-2}, 1195 {0.03, 1.5030E-1, 1.5101E-1, 1.5844E-1, 1.7 << 1195 {0.02, 1.2561E-2, 1.3943E-2, 1.5553E-2, 1.7327E-2, 1.8478E-2, 1.9612E-2, 2.2160E-2, 2.5130E-2, 2.8594E-2}, 1196 {0.04, 2.5743E-1, 2.6774E-1, 2.7855E-1, 3.0 << 1196 {0.03, 3.3750E-2, 3.7470E-2, 4.1528E-2, 4.6170E-2, 4.8708E-2, 5.1401E-2, 5.7297E-2, 6.3943E-2, 7.1441E-2}, 1197 {0.05, 3.7846E-1, 3.9195E-1, 4.0607E-1, 4.3 << 1197 {0.04, 6.9619E-2, 7.6098E-2, 8.3249E-2, 9.1150E-2, 9.5406E-2, 9.9881E-2, 1.0954E-1, 1.2023E-1, 1.3208E-1}, 1198 {0.06, 5.0839E-1, 5.2497E-1, 5.4230E-1, 5.7 << 1198 {0.05, 1.1522E-1, 1.2470E-1, 1.3504E-1, 1.4632E-1, 1.5234E-1, 1.5864E-1, 1.7211E-1, 1.8686E-1, 2.0304E-1}, 1199 {0.08, 7.8230E-1, 8.0474E-1, 8.2818E-1, 8.7 << 1199 {0.06, 1.6931E-1, 1.8179E-1, 1.9530E-1, 2.0991E-1, 2.1767E-1, 2.2576E-1, 2.4295E-1, 2.6165E-1, 2.8201E-1}, 1200 {0.1, 1.0621, 1.0900, 1.1192, 1.1816, 1.250 << 1200 {0.08, 2.9540E-1, 3.1361E-1, 3.3312E-1, 3.5403E-1, 3.6506E-1, 3.7650E-1, 4.0067E-1, 4.2673E-1, 4.5488E-1}, 1201 {0.15, 1.7389, 1.7790, 1.8210, 1.9112, 2.01 << 1201 {0.1, 4.3613E-1, 4.5956E-1, 4.852E-1, 5.1115E-1, 5.2514E-1, 5.3961E-1, 5.7008E-1, 6.0277E-1, 6.3793E-1}, 1202 {0.2, 2.3516, 2.4024, 2.4556, 2.5701, 2.697 << 1202 {0.15, 8.1014E-1, 8.4453E-1, 8.8093E-1, 9.1954E-1, 9.3973E-1, 9.6056E-1, 1.0043, 1.0509, 1.1008}, 1203 {0.3, 3.3741, 3.4427, 3.5148, 3.6706, 3.844 << 1203 {0.2, 1.1888, 1.2319, 1.2774, 1.3255, 1.3506, 1.3765, 1.4308, 1.4886, 1.5504}, 1204 {0.4, 4.1788, 4.2620, 4.3496, 4.5398, 4.753 << 1204 {0.3, 1.8422, 1.8983, 1.9575, 2.0201, 2.0528, 2.0864, 2.1569, 2.2319, 2.3120}, 1205 {0.5, 4.8180, 4.9137, 5.0146, 5.2341, 5.481 << 1205 {0.4, 2.3984, 2.4642, 2.5336, 2.6070, 2.6452, 2.6847, 2.7674, 2.8554, 2.9494}, 1206 {0.6, 5.3765, 5.4830, 5.5954, 5.8406, 6.117 << 1206 {0.5, 2.8181, 2.8915, 2.9690, 3.0509, 3.0937, 3.1378, 3.2301, 3.3285, 3.4337}, 1207 {0.7, 5.8208, 5.9369, 6.0596, 6.3275, 6.630 << 1207 {0.6, 3.1698, 3.2494, 3.3336, 3.4226, 3.4691, 3.5171, 3.6175, 3.7246, 3.8391}, 1208 {0.8, 6.2109, 6.3355, 6.4674, 6.7558, 7.082 << 1208 {0.7, 3.4652, 3.5502, 3.6400, 3.7351, 3.7848, 3.8360, 3.9433, 4.0578, 4.1804}, 1209 {1.0, 6.8747, 7.0142, 7.1621, 7.4861, 7.854 << 1209 {0.8, 3.7392, 3.8289, 3.9236, 4.0239, 4.0764, 4.1304, 4.2438, 4.3648, 4.4944}, 1210 {1.2, 7.3933, 7.5454, 7.7068, 8.0612, 8.465 << 1210 {1.0, 4.2295, 4.3269, 4.4299, 4.5391, 4.5962, 4.6551, 4.7786, 4.9106, 5.0520}, 1211 {1.4, 7.8331, 7.9967, 8.1694, 8.5502, 8.985 << 1211 {1.2, 4.5777, 4.6814, 4.7912, 4.9076, 4.9685, 5.0314, 5.1633, 5.3043, 5.4555}, 1212 {1.5, 8.0286, 8.1967, 8.3753, 8.7681, 9.218 << 1212 {1.4, 4.9001, 5.0092, 5.1247, 5.2473, 5.3114, 5.3776, 5.5166, 5.6653, 5.8249}, 1213 {1.7, 8.3813, 8.5585, 8.7469, 9.1618, 9.636 << 1213 {1.5, 5.0434, 5.1550, 5.2731, 5.3984, 5.4640, 5.5317, 5.6739, 5.8260, 5.9893}, 1214 {2.0, 8.8352, 9.0245, 9.2260, 9.6701, 10.17 << 1214 {1.7, 5.3011, 5.4170, 5.5398, 5.6701, 5.7373, 5.8088, 5.9568, 6.1152, 6.2853}, 1215 {3.0, 9.9511, 10.1714, 10.4062, 10.9254, 11 << 1215 {2.0, 5.6308, 5.7523, 5.8811, 6.0174, 6.0896, 6.1636, 6.3192, 6.4857, 6.6647}, 1216 {5.0, 11.3211, 11.5818, 11.8601, 12.4771, 1 << 1216 {3.0, 6.4224, 6.5580, 6.7019, 6.8549, 6.9351, 7.0180, 7.1925, 7.3795, 7.5808}, 1217 {7.0, 11.9480, 12.2357, 12.5432, 13.2260, 1 << 1217 {5.0, 7.3338, 7.4872, 7.6500, 7.8235, 7.9146, 8.0087, 8.2071, 8.4200, 8.6496}, 1218 }; << 1218 {7.0, 7.7938, 7.9588, 8.1342, 8.3211, 8.4193, 8.5209, 8.7350, 8.9651, 9.2133} 1219 << 1219 }; 1220 G4double b, bs; << 1220 1221 for(i=0; i<nEtaK; ++i) { << 1221 const G4double bls3[28][9] = { 1222 << 1222 {0.005, 1.2895E-3, 1.3670E-3, 1.4524E-3, 1.6524E-3, 1.9078E-3, 2.2414E-3, 2.6889E-3, 3.3006E-3}, 1223 G4double et = Eta[i]; << 1223 {0.007, 2.6467E-3, 2.8242E-3, 3.0238E-3, 3.5045E-3, 4.1260E-3, 4.9376E-3, 6.0050E-3, 7.4152E-3}, 1224 G4double loget = G4Log(et); << 1224 {0.01, 6.1472E-3, 6.6086E-3, 7.1246E-3, 8.3491E-3, 9.8871E-3, 1.1822E-2, 1.4261E-2, 1.7335E-2}, 1225 << 1225 {0.015, 1.63316E-2, 1.7572E-2, 1.8932E-2, 2.2053E-2, 2.5803E-2, 3.0308E-2, 3.5728E-2, 4.2258E-2}, 1226 for(j=0; j<nK; ++j) { << 1226 {0.02, 3.2634E-2, 3.4900E-2, 3.7348E-2, 4.2850E-2, 4.9278E-2, 5.6798E-2, 6.5610E-2, 7.5963E-2}, 1227 << 1227 {0.03, 7.9907E-2, 8.4544E-2, 8.9486E-2, 1.0032E-1, 1.1260E-1, 1.2656E-1, 1.4248E-1, 1.6071E-1}, 1228 if(j < 10) { b = bk2[i][10-j]; } << 1228 {0.04, 1.4523E-1, 1.5237E-1, 1.5985E-1, 1.7614E-1, 1.9434E-1, 2.1473E-1, 2.3766E-1, 2.6357E-1}, 1229 else { b = bk1[i][20-j]; } << 1229 {0.05, 2.2082E-1, 2.3036E-1, 2.4038E-1, 2.6199E-1, 2.8590E-1, 3.1248E-1, 3.4212E-1, 3.7536E-1}, 1230 << 1230 {0.06, 3.0423E-1, 3.1611E-1, 3.2854E-1, 3.5522E-1, 3.8459E-1, 4.1704E-1, 4.5306E-1, 4.9326E-1}, 1231 CK[j][i] = SK[j]*loget + TK[j] - b; << 1231 {0.08, 4.8536E-1, 5.0156E-1, 5.1846E-1, 5.5453E-1, 5.9397E-1, 6.3728E-1, 6.8507E-1, 7.3810E-1}, 1232 << 1232 {0.1, 6.7586E-1, 6.9596E-1, 7.1688E-1, 7.6141E-1, 8.0992E-1, 8.6301E-1, 9.2142E-1, 9.8604E-1}, 1233 if(i == nEtaK-1) { << 1233 {0.15, 1.1544, 1.1828, 1.2122, 1.2746, 1.3424, 1.4163, 1.4974, 1.5868}, 1234 ZK[j] = et*(et*et*CK[j][i] - et*UK[j] << 1234 {0.2, 1.6167, 1.6517, 1.6880, 1.7650, 1.8484, 1.9394, 2.0390, 2.1489}, 1235 //G4cout << "i= " << i << " j= " << j << 1235 {0.3, 2.3979, 2.4432, 2.4902, 2.5899, 2.6980, 2.8159, 2.9451, 3.0876}, 1236 // << " CK[j][i]= " << CK[j][i << 1236 {0.4, 3.0502, 3.1034, 3.1586, 3.2758, 3.4030, 3.5416, 3.6938, 3.8620}, 1237 // << " ZK[j]= " << ZK[j] << " << 1237 {0.5, 3.5466, 3.6062, 3.6681, 3.7994, 3.9421, 4.0978, 4.2688, 4.4580}, 1238 } << 1238 {0.6, 3.9620, 4.0270, 4.0945, 4.2378, 4.3935, 4.5636, 4.7506, 4.9576}, 1239 } << 1239 {0.7, 4.3020, 4.3715, 4.4438, 4.5974, 4.7644, 4.9470, 5.1478, 5.3703}, 1240 //G4cout << G4endl; << 1240 {0.8, 4.6336, 4.7072, 4.7837, 4.9463, 5.1233, 5.3169, 5.5300, 5.7661}, 1241 if(i < nEtaL) { << 1241 {1.0, 5.2041, 5.2845, 5.3682, 5.5462, 5.7400, 5.9523, 6.1863, 6.4458}, 1242 //G4cout << " LShell:" <<G4endl; << 1242 {1.2, 5.6182, 5.7044, 5.7940, 5.9848, 6.1927, 6.4206, 6.6719, 6.9510}, 1243 for(j=0; j<nL; ++j) { << 1243 {1.4, 5.9967, 6.0876, 6.1823, 6.3839, 6.6038, 6.8451, 7.1113, 7.4071}, 1244 << 1244 {1.5, 6.1652, 6.2583, 6.3553, 6.5618, 6.7871, 7.0344, 7.3073, 7.6107}, 1245 if(j < 8) { << 1245 {1.7, 6.4686, 6.5657, 6.6668, 6.8823, 7.1175, 7.3757, 7.6609, 7.9782}, 1246 bs = bls3[i][8-j]; << 1246 {2.0, 6.8577, 6.9600, 7.0666, 7.2937, 7.5418, 7.8143, 8.1156, 8.4510}, 1247 b = bll3[i][8-j]; << 1247 {3.0, 7.7981, 7.9134, 8.0336, 8.2901, 8.5708, 8.8796, 9.2214, 9.6027}, 1248 } else if(j < 17) { << 1248 {5.0, 8.8978, 9.0297, 9.1673, 9.4612, 9.7834, 10.1384, 10.5323, 10.9722}, 1249 bs = bls2[i][17-j]; << 1249 {7.0, 9.4819, 9.6248, 9.7739, 10.0926, 10.4423, 10.8282, 11.2565, 11.7356} 1250 b = bll2[i][17-j]; << 1250 }; 1251 } else { << 1251 1252 bs = bls1[i][26-j]; << 1252 const G4double bll1[28][10] = { 1253 b = bll1[i][26-j]; << 1253 {0.005, 3.6324E-5, 4.0609E-5, 4.5430E-5, 5.6969E-5, 7.1625E-5, 9.0279E-5, 1.1407E-4, 1.2834E-4, 1.4447E-4}, 1254 } << 1254 {0.007, 1.8110E-4, 2.0001E-4, 2.2099E-4, 2.7006E-4, 3.3049E-4, 4.0498E-4, 4.9688E-4, 5.5061E-4, 6.1032E-4}, 1255 G4double c = SL[j]*loget + TL[j]; << 1255 {0.01, 8.6524E-4, 9.4223E-4, 1.0262E-3, 1.2178E-3, 1.4459E-3, 1.7174E-3, 2.0405E-3, 2.2245E-3, 2.4252E-3}, 1256 CL[j][i] = c - bs - 3.0*b; << 1256 {0.015, 4.2293E-3, 4.5314E-3, 4.8551E-3, 5.5731E-3, 6.3968E-3, 7.3414E-3, 8.4242E-3, 9.0236E-3, 9.6652E-3}, 1257 if(i == nEtaL-1) { << 1257 {0.02, 1.1485E-2, 1.2172E-2, 1.2900E-2, 1.4486E-2, 1.6264E-2, 1.8256E-2, 2.0487E-2, 2.1702E-2, 2.2989E-2}, 1258 VL[j] = et*(et*CL[j][i] - UL[j]); << 1258 {0.03, 3.9471E-2, 4.1270E-2, 4.3149E-2, 4.7163E-2, 5.1543E-2, 5.6423E-2, 6.1540E-2, 6.4326E-2, 6.7237E-2}, 1259 //G4cout << "i= " << i << " j= " << << 1259 {0.04, 8.4454E-2, 8.7599E-2, 9.0860E-2, 9.7747E-2, 1.0516E-1, 1.1313E-1, 1.2171E-1, 1.2625E-1, 1.3096E-1}, 1260 // << " CL[j][i]= " << CL[ << 1260 {0.05, 1.4339E-1, 1.4795E-1, 1.5266E-1, 1.6253E-1, 1.7306E-1, 1.8430E-1, 1.9630E-1, 2.0261E-1, 2.0924E-1}, 1261 // << " VL[j]= " << VL[j] < << 1261 {0.06, 2.1304E-1, 2.1899E-1, 2.2512E-1, 2.3794E-1, 2.5153E-1, 2.6596E-1, 2.8130E-1, 2.8934E-1, 2.9763E-1}, 1262 // << " et= " << et << G4en << 1262 {0.08, 3.7382E-1, 3.8241E-1, 3.9122E-1, 4.0955E-1, 4.2888E-1, 4.4928E-1, 4.7086E-1, 4.8212E-1, 4.9371E-1}, 1263 //" UL= " << UL[j] << " TL= " << TL << 1263 {0.1, 5.5056E-1, 5.6151E-1, 5.7273E-1, 5.9601E-1, 6.2049E-1, 6.4627E-1, 6.7346E-1, 6.8762E-1, 7.0218E-1}, 1264 } << 1264 {0.15, 1.0066, 1.0224, 1.0386, 1.0721, 1.1073, 1.1443, 1.1832, 1.2035, 1.2243}, 1265 } << 1265 {0.2, 1.4376, 1.4572, 1.4773, 1.5188, 1.5624, 1.6083, 1.6566, 1.6817, 1.7076}, 1266 //G4cout << G4endl; << 1266 {0.3, 2.1712, 2.1964, 2.2223, 2.2758, 2.3322, 2.3915, 2.4542, 2.4869, 2.5205}, 1267 } << 1267 {0.4, 2.7500, 2.7793, 2.8094, 2.8719, 2.9377, 3.0072, 3.0807, 3.1192, 3.1587}, 1268 } << 1268 {0.5, 3.2033, 3.2359, 3.2693, 3.3389, 3.4122, 3.4898, 3.5721, 3.6151, 3.6595}, 1269 << 1269 {0.6, 3.6038, 3.6391, 3.6753, 3.7506, 3.8303, 3.9146, 4.0042, 4.0511, 4.0995}, 1270 const G4double xzk[34] = { 11.7711, << 1270 {0.7, 3.9106, 3.9482, 3.9867, 4.0670, 4.1520, 4.2421, 4.3380, 4.3882, 4.4401}, 1271 13.3669, 15.5762, 17.1715, 18.7667, 20.85 << 1271 {0.8, 4.1790, 4.2185, 4.2590, 4.3437, 4.4333, 4.5285, 4.6298, 4.6830, 4.7380}, 1272 34.3457, 37.4119, 40.3555, 42.3177, 44.77 << 1272 {1.0, 4.6344, 4.6772, 4.7212, 4.8131, 4.9106, 5.0144, 5.1250, 5.1831, 5.2432}, 1273 60.9586, 63.6567, 66.5998, 68.807, 71.872 << 1273 {1.2, 4.9787, 5.0242, 5.0711, 5.1689, 5.2729, 5.3837, 5.5050, 5.5642, 5.6287}, 1274 93.4549, 96.2753 << 1274 {1.4, 5.2688, 5.3166, 5.3658, 5.4687, 5.5782, 5.6950, 5.8198, 5.8855, 5.9554}, 1275 const G4double yzk[34] = { 0.70663, << 1275 {1.5, 5.3966, 5.4454, 5.4957, 5.6009, 5.7128, 5.8323, 5.9601, 6.0274, 6.0972}, 1276 0.72033, 0.73651, 0.74647, 0.75518, 0.763 << 1276 {1.7, 5.6255, 5.6762, 5.7284, 5.8377, 5.9541, 6.0785, 6.2116, 6.2818, 6.3546}, 1277 0.80611, 0.8123, 0.8185, 0.82097, 0.82467 << 1277 {2.0, 5.9170, 5.9701, 6.0248, 6.1395, 6.2618, 6.3925, 6.5327, 6.6066, 6.6833}, 1278 0.84565, 0.84936, 0.85181, 0.85303, 0.855 << 1278 {3.0, 6.6210, 6.6801, 6.7411, 6.8692, 7.0062, 7.1529, 7.3107, 7.3941, 7.4807}, 1279 0.86891, 0.87011 << 1279 {5.0, 7.4620, 7.5288, 7.5977, 7.7428, 7.8982, 8.0653, 8.2454, 8.3409, 8.4402}, 1280 << 1280 {7.0, 7.7362, 7.8079, 7.8821, 8.0383, 8.2061, 8.3866, 8.5816, 8.6850, 8.7927} 1281 const G4double xzl[36] = { 15.5102, << 1281 }; 1282 16.7347, 17.9592, 19.551, 21.0204, 22.612 << 1282 1283 34.4898, 36.2041, 38.4082, 40.3674, 42.57 << 1283 const G4double bll2[28][10] = { 1284 59.3469, 61.9184, 64.6122, 67.4286, 71.46 << 1284 {0.005, 1.8339E-4, 2.3330E-4, 2.9738E-4, 3.7977E-4, 4.2945E-4, 4.8582E-4, 6.2244E-4, 7.9858E-4, 1.0258E-3}, 1285 91.551, 94.2449, << 1285 {0.007, 7.5042E-4, 9.2355E-4, 1.1375E-3, 1.4021E-3, 1.5570E-3, 1.7292E-3, 2.1335E-3, 2.6335E-3, 3.2515E-3}, 1286 const G4double yzl[36] = { 0.29875, << 1286 {0.01, 2.8829E-3, 3.4275E-3, 4.0758E-3, 4.8457E-3, 5.2839E-3, 5.7617E-3, 6.8504E-3, 8.1442E-3, 9.6816E-3}, 1287 0.31746, 0.33368, 0.35239, 0.36985, 0.387 << 1287 {0.015, 1.1087E-2, 1.2716E-2, 1.4581E-2, 1.6717E-2, 1.7898E-2, 1.9163E-2, 2.1964E-2, 2.5173E-2, 2.8851E-2}, 1288 0.4896, 0.50083, 0.51331, 0.52328, 0.5307 << 1288 {0.02, 2.5786E-2, 2.8922E-2, 3.2435E-2, 3.6371E-2, 3.8514E-2, 4.0784E-2, 4.5733E-2, 5.1288E-2, 5.7531E-2}, 1289 0.58191, 0.5869, 0.59189, 0.60062, 0.6068 << 1289 {0.03, 7.3461E-2, 8.0264E-2, 8.7705E-2, 9.5852E-2, 1.0021E-1, 1.0478E-1, 1.1458E-1, 1.2535E-1, 1.3721E-1}, 1290 0.6368, 0.64054 << 1290 {0.04, 1.4094E-1, 1.5172E-1, 1.6336E-1, 1.7596E-1, 1.8265E-1, 1.8962E-1, 2.0445E-1, 2.2058E-1, 2.3818E-1}, 1291 << 1291 {0.05, 2.2289E-1, 2.3762E-1, 2.5344E-1, 2.7045E-1, 2.7944E-1, 2.8877E-1, 3.0855E-1, 3.2995E-1, 3.5318E-1}, 1292 sThetaK = new G4PhysicsFreeVector(34, false << 1292 {0.06, 3.1503E-1, 3.3361E-1, 3.5346E-1, 3.7473E-1, 3.8594E-1, 3.9756E-1, 4.2212E-1, 4.4861E-1, 4.7727E-1}, 1293 for(i=0; i<34; ++i) { sThetaK->PutValues(i, << 1293 {0.08, 5.1793E-1, 5.4368E-1, 5.7109E-1, 6.0032E-1, 6.1569E-1, 6.3159E-1, 6.6512E-1, 7.0119E-1, 7.4012E-1}, 1294 sThetaL = new G4PhysicsFreeVector(36, false << 1294 {0.1, 7.3258E-1, 7.6481E-1, 7.9907E-1, 8.3556E-1, 8.5472E-1, 8.7454E-1, 9.1630E-1, 9.6119E-1, 1.0096}, 1295 for(i=0; i<36; ++i) { sThetaL->PutValues(i, << 1295 {0.15, 1.2677, 1.3137, 1.3626, 1.4147, 1.4421, 1.4703, 1.5300, 1.5942, 1.6636}, 1296 } << 1296 {0.2, 1.7615, 1.8188, 1.8797, 1.9446, 1.9788, 2.0142, 2.0889, 2.1694, 2.2567}, 1297 << 1297 {0.3, 2.5909, 2.6658, 2.7457, 2.8312, 2.8762, 2.9231, 3.0222, 3.1295, 3.2463}, 1298 //....oooOO0OOooo........oooOO0OOooo........o << 1298 {0.4, 3.2417, 3.3302, 3.4249, 3.5266, 3.5803, 3.6361, 3.7546, 3.8835, 4.0242}, 1299 << 1299 {0.5, 3.7527, 3.8523, 3.9591, 4.0741, 4.1350, 4.1983, 4.3330, 4.4799, 4.6408}, >> 1300 {0.6, 4.2013, 4.3103, 4.4274, 4.5537, 4.6206, 4.6904, 4.8390, 5.0013, 5.1796}, >> 1301 {0.7, 4.5493, 4.6664, 4.7925, 4.9286, 5.0009, 5.0762, 5.2370, 5.4129, 5.6066}, >> 1302 {0.8, 4.8537, 4.9780, 5.1119, 5.2568, 5.3338, 5.4141, 5.5857, 5.7738, 5.9811}, >> 1303 {1.0, 5.3701, 5.5066, 5.6540, 5.8138, 5.8989, 5.9878, 6.1780, 6.3870, 6.6179}, >> 1304 {1.2, 5.7648, 5.9114, 6.0701, 6.2424, 6.3343, 6.4303, 6.6361, 6.8626, 7.1137}, >> 1305 {1.4, 6.0976, 6.2530, 6.4214, 6.6044, 6.7021, 6.8043, 7.0237, 7.2655, 7.5338}, >> 1306 {1.5, 6.2447, 6.4041, 6.5768, 6.7647, 6.8650, 6.9700, 7.1954, 7.4442, 7.7203}, >> 1307 {1.7, 6.5087, 6.6752, 6.8558, 7.0526, 7.1578, 7.2679, 7.5045, 7.7660, 8.0565}, >> 1308 {2.0, 6.8458, 7.0218, 7.2129, 7.4213, 7.5328, 7.6496, 7.9010, 8.1791, 8.4886}, >> 1309 {3.0, 7.6647, 7.8644, 8.0819, 8.3189, 8.4475, 8.5814, 8.8702, 9.1908, 9.5488}, >> 1310 {5.0, 8.6515, 8.8816, 9.1330, 9.4090, 9.5572, 9.7132, 10.0504, 10.4259, 10.8466}, >> 1311 {7.0, 9.0221, 9.2724, 9.5464, 9.8477, 10.0099, 10.1805, 10.5499, 10.9622, 11.4250} >> 1312 }; >> 1313 >> 1314 const G4double bll3[28][9] = { >> 1315 {0.005, 1.3190E-3, 1.4961E-3, 1.6974E-3, 2.1858E-3, 2.8163E-3, 3.6302E-3, 4.6814E-3, 6.0395E-3}, >> 1316 {0.007, 4.0158E-3, 4.4623E-3, 4.9592E-3, 6.1257E-3, 7.5675E-3, 9.3502E-3, 1.1556E-2, 1.4290E-2}, >> 1317 {0.01, 1.1509E-2, 1.2548E-2, 1.3681E-2, 1.6263E-2, 1.9336E-2, 2.2999E-2, 2.7370E-2, 3.2603E-2}, >> 1318 {0.015, 3.3070E-2, 3.5408E-2, 3.7914E-2, 4.3483E-2, 4.9898E-2, 5.7304E-2, 6.5884E-2, 7.5861E-2}, >> 1319 {0.02, 6.4555E-2, 6.8394E-2, 7.2472E-2, 8.1413E-2, 9.1539E-2, 1.0304E-1, 1.1617E-1, 1.3121E-1}, >> 1320 {0.03, 1.5030E-1, 1.5101E-1, 1.5844E-1, 1.7451E-1, 1.9244E-1, 2.1244E-1, 2.3496E-1, 2.6044E-1}, >> 1321 {0.04, 2.5743E-1, 2.6774E-1, 2.7855E-1, 3.0180E-1, 3.2751E-1, 3.5608E-1, 3.8803E-1, 4.2401E-1}, >> 1322 {0.05, 3.7846E-1, 3.9195E-1, 4.0607E-1, 4.3635E-1, 4.6973E-1, 5.0672E-1, 5.4798E-1, 5.9436E-1}, >> 1323 {0.06, 5.0839E-1, 5.2497E-1, 5.4230E-1, 5.7943E-1, 6.2028E-1, 6.6549E-1, 7.1589E-1, 7.7252E-1}, >> 1324 {0.08, 7.8230E-1, 8.0474E-1, 8.2818E-1, 8.7836E-1, 9.3355E-1, 9.9462E-1, 1.0627, 1.1394}, >> 1325 {0.1, 1.0621, 1.0900, 1.1192, 1.1816, 1.2503, 1.3265, 1.4116, 1.5076}, >> 1326 {0.15, 1.7389, 1.7790, 1.8210, 1.9112, 2.0108, 2.1217, 2.2462, 2.3876}, >> 1327 {0.2, 2.3516, 2.4024, 2.4556, 2.5701, 2.6971, 2.8391, 2.9994, 3.1822}, >> 1328 {0.3, 3.3741, 3.4427, 3.5148, 3.6706, 3.8445, 4.0404, 4.2631, 4.5193}, >> 1329 {0.4, 4.1788, 4.2620, 4.3496, 4.5398, 4.7530, 4.9944, 5.2703, 5.5895}, >> 1330 {0.5, 4.8180, 4.9137, 5.0146, 5.2341, 5.4811, 5.7618, 6.0840, 6.4583}, >> 1331 {0.6, 5.3765, 5.4830, 5.5954, 5.8406, 6.1173, 6.4326, 6.7958, 7.2191}, >> 1332 {0.7, 5.8208, 5.9369, 6.0596, 6.3275, 6.6306, 6.9769, 7.3767, 7.8440}, >> 1333 {0.8, 6.2109, 6.3355, 6.4674, 6.7558, 7.0827, 7.4570, 7.8900, 8.3972}, >> 1334 {1.0, 6.8747, 7.0142, 7.1621, 7.4861, 7.8546, 8.2778, 8.7690, 9.3464}, >> 1335 {1.2, 7.3933, 7.5454, 7.7068, 8.0612, 8.4652, 8.9302, 9.4713, 10.1090}, >> 1336 {1.4, 7.8331, 7.9967, 8.1694, 8.5502, 8.9851, 9.4866, 10.0713, 10.7619}, >> 1337 {1.5, 8.0286, 8.1967, 8.3753, 8.7681, 9.2181, 9.7352, 10.3399, 11.0546}, >> 1338 {1.7, 8.3813, 8.5585, 8.7469, 9.1618, 9.6367, 10.1856, 10.8270, 11.5863}, >> 1339 {2.0, 8.8352, 9.0245, 9.2260, 9.6701, 10.1793, 10.7688, 11.4590, 12.2775}, >> 1340 {3.0, 9.9511, 10.1714, 10.4062, 10.9254, 11.5229, 12.2172, 13.0332, 14.0048}, >> 1341 {5.0, 11.3211, 11.5818, 11.8601, 12.4771, 13.1898, 14.0213, 15.0024, 16.1752}, >> 1342 {7.0, 11.9480, 12.2357, 12.5432, 13.2260, 14.0164, 14.9404, 16.0330, 17.3420} >> 1343 }; >> 1344 >> 1345 >> 1346 G4double b, bs; >> 1347 for(i=0; i<nEtaK; i++) { >> 1348 >> 1349 G4double et = eta[i]; >> 1350 G4double loget = std::log(et); >> 1351 Eta[i] = et; >> 1352 // G4cout << "### eta[" << i << "]= " << et << " KShell: tet= " << TheK[0]<<" - " <<TheK[nK-1]<< G4endl; >> 1353 >> 1354 for(j=0; j<nK; j++) { >> 1355 >> 1356 if(j < 10) b = bk2[i][10-j]; >> 1357 else b = bk1[i][20-j]; >> 1358 >> 1359 CK[j][i] = SK[j]*loget + TK[j] - b; >> 1360 //G4cout << " " << CK[j][i]; >> 1361 if(i == nEtaK-1) ZK[j] = et*(et*et*CK[j][i] - et*UK[j] - VK[j]); >> 1362 } >> 1363 //G4cout << G4endl; >> 1364 if(i < nEtaL) { >> 1365 //G4cout << " LShell:" <<G4endl; >> 1366 for(j=0; j<nL; j++) { >> 1367 >> 1368 if(j < 8) { >> 1369 bs = bls3[i][8-j]; >> 1370 b = bll3[i][8-j]; >> 1371 } else if(j < 17) { >> 1372 bs = bls2[i][17-j]; >> 1373 b = bll2[i][17-j]; >> 1374 } else { >> 1375 bs = bls1[i][26-j]; >> 1376 b = bll1[i][26-j]; >> 1377 } >> 1378 G4double c = SL[j]*loget + TL[j]; >> 1379 CL[j][i] = c - bs - 3.0*b; >> 1380 //G4cout << " " << CL[j][i]; >> 1381 if(i == nEtaL-1) VL[j] = et*(et*CL[j][i] - UL[j]); >> 1382 } >> 1383 //G4cout << G4endl; >> 1384 } >> 1385 } >> 1386 const G4double hm[53] = { >> 1387 12.0, 12.0, 12.0, 12.0, 11.9, 11.7, 11.5, 11.2, 10.8, 10.4, >> 1388 10.0, 9.51, 8.97, 8.52, 8.03, 7.46, 6.95, 6.53, 6.18, 5.87, >> 1389 5.61, 5.39, 5.19, 5.01, 4.86, 4.72, 4.62, 4.53, 4.44, 4.38, >> 1390 4.32, 4.26, 4.20, 4.15, 4.1, 4.04, 4.00, 3.95, 3.93, 3.91, >> 1391 3.90, 3.89, 3.89, 3.88, 3.88, 3.88, 3.88, 3.88, 3.89, 3.89, >> 1392 3.90, 3.92, 3.93 }; >> 1393 const G4double hn[31] = { >> 1394 75.5, 61.9, 52.2, 45.1, 39.6, 35.4, 31.9, 29.1, 27.2, 25.8, >> 1395 24.5, 23.6, 22.7, 22.0, 21.4, 20.9, 20.5, 20.2, 19.9, 19.7, >> 1396 19.5, 19.3, 19.2, 19.1, 18.4, 18.8, 18.7, 18.6, 18.5, 18.4, >> 1397 18.2 >> 1398 }; >> 1399 for(i=0; i<53; i++) {HM[i] = hm[i];} >> 1400 for(i=0; i<31; i++) {HN[i] = hn[i];} >> 1401 >> 1402 const G4double mm[93] = { >> 1403 0.0, >> 1404 /* >> 1405 -0.0001577, 5.396e-05, 0.0004194, -0.0001969, -0.0003447, -0.0001904, 9.612e-05, 4.16e-05, 0.0001899, 0.0001322, >> 1406 0.0001485, 4.698e-05, -7.726e-05, -6.448e-05, 3.269e-05, -0.0001293, 0.0001483, -8.559e-05, 4.018e-05, 6.239e-05, >> 1407 4.447e-05, 2.212e-05, -5.568e-06, 3.782e-05, 4.398e-05, 8.942e-05, 0.0002202, 0.0001357, 0.0001195, 0.0001812, >> 1408 0.0002365, 4.508e-05, 0.0001629, 7.75e-05, 0.0001886, 0.0001601, 0.00015, 5.679e-05, 5.956e-06, -7.309e-05, >> 1409 -0.0001144, -8.585e-05, -0.0001051, -8.209e-05, -5.871e-05, -3.744e-05, -6.707e-05, -9.199e-06, 1.181e-05, 2.252e-05, >> 1410 1.051e-05, 7.63e-05, 1.71e-05, 1.705e-05, 1.561e-05, 8.987e-06, 4.957e-06, -1.473e-05, -1.902e-05, -1.491e-05, >> 1411 -9.891e-08, 2.584e-05, 3.696e-05, 8.651e-06, -3.601e-06, 1.647e-05, 4.536e-05, 9.954e-05, 8.922e-05, 0.0001026, >> 1412 4.094e-05, 3.788e-05, 1.578e-05, 2.731e-05, -6.035e-05, -9.112e-05, -8.846e-05, -0.0001361, -0.0001959, -0.0001883, >> 1413 -0.0002854, -0.0004129, -0.0004611, -0.0005523, -0.0005869, -0.0005008, -0.000659, -0.0007045, -0.0007322, -0.0008374, >> 1414 -0.0008731, -0.0009327, >> 1415 >> 1416 -0.0003427, -4.685e-05, 0.0007304, -0.0003158, -0.0004104, -0.001133, -0.0007987, -0.0001528, 8.976e-05, 8.63e-05, >> 1417 1.764e-05, -0.000136, -0.0002895, -0.0002618, -0.0002878, -0.0003813, 1.417e-05, -0.0004216, -0.0003859, -0.0001028, >> 1418 -0.0001617, -0.0002141, -0.0002178, 5.59e-05, 9.964e-07, 7.528e-05, 0.0002261, 4.241e-05, 0.0005255, 0.0002139, >> 1419 -5.821e-05, -2.661e-05, 0.0001039, 3.359e-05, 0.0004229, -0.0004872, 0.0001374, -2.79e-05, -1.319e-05, -0.0001077, >> 1420 -0.0001923, -0.0001189, -0.0001091, -8.486e-05, 7.506e-05, 0.0001053, 0.0004596, 0.000177, 4.805e-05, 3.97e-05, >> 1421 0.0004008, 0.0002194, 0.0001639, 0.0003285, 0.000219, 0.0001339, 6.323e-05, 2.013e-05, 2.568e-05, 4.346e-05, >> 1422 8.455e-05, 8.785e-05, 0.0001393, 9.907e-05, 0.0001003, 0.0001508, 0.0002189, 0.0003885, 0.0003603, 0.0003839, >> 1423 0.0003321, 0.0002207, 0.0001627, 0.0002282, 0.0002177, 0.0002156, 0.0002855, 0.0002361, 0.0007735, 0.0001157, >> 1424 -0.0001214, -0.0002944, -0.0003193, -0.0004025, -0.0002436, -7.573e-05, -0.0003741, -0.0003678, -0.0004839, -0.0003934, >> 1425 -0.0005817, -0.000726, >> 1426 >> 1427 -118, 40.37, 313.8, -147.3, -257.9, -142.4, 71.91, 31.12, 142.1, 98.92, >> 1428 112.3, 41.36, -47.15, -17.84, 46.46, -66.91, 141.3, -19.96, 84.66, 110.2, >> 1429 113.2, 108.7, 92.89, 127.3, 139.3, 187, 288.5, 233.4, 225.3, 284.6, >> 1430 344.5, 246.2, 359.7, 324.1, 455, 457.1, 481, 459.8, 447.9, 407.5, >> 1431 399.7, 457.7, 470.4, 511.8, 545.9, 581, 574.1, 638.6, 685.6, 719.6, >> 1432 740.8, 817.2, 808.3, 839.7, 866.2, 873.2, 879.3, 886.3, 902, 912.7, >> 1433 927.7, 942.1, 967.1, 962.2, 961.8, 993.4, 1006, 1060, 1046, 1061, >> 1434 1036, 1060, 1062, 1084, 1048, 1063, 1098, 1093, 1089, 1117, >> 1435 1085, 1044, 1033, 1027, 1051, 1150, 1064, 1078, 1090, 1140, >> 1436 1119, 1097, >> 1437 */ >> 1438 >> 1439 -511.3, -69.91, 1090, -471.2, -612.4, -1690, -1192, -228, 133.9, 128.8, >> 1440 28.71, -190.6, -408.7, -326.6, -378.5, -507.3, 93.63, -531, -459.2, -21.83, >> 1441 -74.93, -121.7, -122.8, 285.2, 217.4, 366.3, 598.7, 341.8, 1078, 634.2, >> 1442 280.5, 430.7, 663.5, 596.4, 1307, 6.814, 992.6, 887, 967.6, 868, >> 1443 775.9, 1034, 1050, 1204, 1458, 1541, 2186, 1806, 1676, 1718, >> 1444 2416, 2204, 2197, 2496, 2371, 2370, 2330, 2318, 2379, 2434, >> 1445 2513, 2533, 2650, 2650, 2684, 2788, 2868, 3157, 3150, 3196, >> 1446 3162, 3055, 3008, 3133, 3176, 3212, 3356, 3324, 4170, 3320, >> 1447 3109, 2941, 2973, 2858, 3088, 3447, 3134, 3121, 3086, 3199, >> 1448 3115, 2979, >> 1449 >> 1450 }; >> 1451 >> 1452 const G4double ntau[93] = { >> 1453 0.0, >> 1454 -251.1, 512.8, 1724, 649, 658.8, 594.1, 946.6, 969.2, 1154, 997.8, >> 1455 939.3, 775.5, 619.9, 848.3, 973.4, 1035, 1252, 698.8, 1043, 1080, >> 1456 957, 885.6, 772.1, 847.1, 799.7, 901.7, 1042, 880.8, 857.7, 968, >> 1457 850.4, 482.8, 781.4, 685, 1093, 847.4, 824.4, 751.5, 638, 489.2, >> 1458 266.5, 371, 319.2, 356.6, 424.9, 405, 390.8, 610.6, 404.3, 444.6, >> 1459 439, 579, 466.4, 520.5, 551.1, 441.6, 328.8, 264.3, 260.9, 274, >> 1460 264.8, 235.1, 290, 196, 178.2, 248.4, 224.3, 268.2, 239.5, 253.9, >> 1461 219.2, 172.7, 169.5, 199, 53.35, 65.78, 120.2, 77.13, 32.76, 45.53, >> 1462 -148, -331.4, -386.4, -430.8, -361.1, -108.3, -382.4, -400.1, -458.4, -430.8, >> 1463 -536.9, -627.6, >> 1464 }; >> 1465 >> 1466 for(i=0; i<93; i++) { >> 1467 MSH[i] = mm[i]; >> 1468 TAU[i] = ntau[i]; >> 1469 } >> 1470 const G4double coseb[14] = {0.0,0.05,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8, >> 1471 1.0,1.2,1.5,2.0}; >> 1472 const G4double cosxi[14] = {1.0000, 0.9905, 0.9631, 0.9208, 0.8680, >> 1473 0.7478, 0.6303, 0.5290, 0.4471, 0.3323, >> 1474 0.2610, 0.2145, 0.1696, 0.1261}; >> 1475 for(i=0; i<14; i++) { >> 1476 COSEB[i] = coseb[i]; >> 1477 COSXI[i] = cosxi[i]; >> 1478 } >> 1479 >> 1480 for(i=1; i<100; i++) { >> 1481 Z23[i] = std::pow(G4double(i),0.23); >> 1482 } >> 1483 } >> 1484 >> 1485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 1486 1300 1487