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 // ------------------------------------------------------------------- 27 // 27 // 28 // GEANT4 Class header file 28 // GEANT4 Class header file 29 // 29 // 30 // 30 // 31 // File name: G4LindhardSorensenIonModel 31 // File name: G4LindhardSorensenIonModel 32 // 32 // 33 // Author: Alexander Bagulya & Vladimir 33 // Author: Alexander Bagulya & Vladimir Ivanchenko 34 // 34 // 35 // Creation date: 16.04.2018 35 // Creation date: 16.04.2018 36 // 36 // 37 // 37 // 38 // ------------------------------------------- 38 // ------------------------------------------------------------------- 39 // 39 // 40 40 41 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 42 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 43 43 44 #include "G4LindhardSorensenIonModel.hh" 44 #include "G4LindhardSorensenIonModel.hh" 45 #include "Randomize.hh" 45 #include "Randomize.hh" 46 #include "G4PhysicalConstants.hh" 46 #include "G4PhysicalConstants.hh" 47 #include "G4SystemOfUnits.hh" 47 #include "G4SystemOfUnits.hh" 48 #include "G4Electron.hh" 48 #include "G4Electron.hh" 49 #include "G4LossTableManager.hh" 49 #include "G4LossTableManager.hh" 50 #include "G4EmCorrections.hh" 50 #include "G4EmCorrections.hh" 51 #include "G4ParticleChangeForLoss.hh" 51 #include "G4ParticleChangeForLoss.hh" 52 #include "G4Log.hh" 52 #include "G4Log.hh" 53 #include "G4DeltaAngle.hh" 53 #include "G4DeltaAngle.hh" 54 #include "G4LindhardSorensenData.hh" 54 #include "G4LindhardSorensenData.hh" 55 #include "G4BraggModel.hh" << 55 #include "G4BraggIonModel.hh" 56 #include "G4BetheBlochModel.hh" 56 #include "G4BetheBlochModel.hh" 57 #include "G4IonICRU73Data.hh" 57 #include "G4IonICRU73Data.hh" 58 #include "G4AutoLock.hh" << 59 58 60 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 61 60 62 using namespace std; 61 using namespace std; 63 62 64 G4LindhardSorensenData* G4LindhardSorensenIonM 63 G4LindhardSorensenData* G4LindhardSorensenIonModel::lsdata = nullptr; 65 G4IonICRU73Data* G4LindhardSorensenIonModel::f 64 G4IonICRU73Data* G4LindhardSorensenIonModel::fIonData = nullptr; 66 << 65 std::vector<G4float>* G4LindhardSorensenIonModel::fact[] = {nullptr}; 67 namespace << 68 { << 69 G4Mutex ionXSMutex = G4MUTEX_INITIALIZER; << 70 } << 71 66 72 G4LindhardSorensenIonModel::G4LindhardSorensen 67 G4LindhardSorensenIonModel::G4LindhardSorensenIonModel(const G4ParticleDefinition*, 73 68 const G4String& nam) 74 : G4VEmModel(nam), 69 : G4VEmModel(nam), >> 70 particle(nullptr), 75 twoln10(2.0*G4Log(10.0)) 71 twoln10(2.0*G4Log(10.0)) 76 { 72 { >> 73 fParticleChange = nullptr; 77 theElectron = G4Electron::Electron(); 74 theElectron = G4Electron::Electron(); 78 corr = G4LossTableManager::Instance()->EmCor 75 corr = G4LossTableManager::Instance()->EmCorrections(); 79 nist = G4NistManager::Instance(); 76 nist = G4NistManager::Instance(); 80 fBraggModel = new G4BraggModel(); << 77 fBraggModel = new G4BraggIonModel(); 81 fBBModel = new G4BetheBlochModel(); 78 fBBModel = new G4BetheBlochModel(); 82 fElimit = 2.0*CLHEP::MeV; 79 fElimit = 2.0*CLHEP::MeV; 83 } 80 } 84 81 85 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 86 83 87 G4LindhardSorensenIonModel::~G4LindhardSorense << 84 G4LindhardSorensenIonModel::~G4LindhardSorensenIonModel() = default; 88 if(isFirst) { << 89 delete lsdata; << 90 delete fIonData; << 91 lsdata = nullptr; << 92 fIonData = nullptr; << 93 } << 94 } << 95 85 96 //....oooOO0OOooo........oooOO0OOooo........oo 86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 97 87 98 void G4LindhardSorensenIonModel::Initialise(co 88 void G4LindhardSorensenIonModel::Initialise(const G4ParticleDefinition* p, 99 co 89 const G4DataVector& ptr) 100 { 90 { 101 fBraggModel->Initialise(p, ptr); 91 fBraggModel->Initialise(p, ptr); 102 fBBModel->Initialise(p, ptr); 92 fBBModel->Initialise(p, ptr); 103 SetParticle(p); 93 SetParticle(p); >> 94 //G4cout << "G4LindhardSorensenIonModel::Initialise for " >> 95 // << p->GetParticleName() << G4endl; 104 96 105 // always false before the run 97 // always false before the run 106 SetDeexcitationFlag(false); 98 SetDeexcitationFlag(false); 107 99 108 if(nullptr == fParticleChange) { 100 if(nullptr == fParticleChange) { 109 fParticleChange = GetParticleChangeForLoss 101 fParticleChange = GetParticleChangeForLoss(); 110 if(UseAngularGeneratorFlag() && nullptr == 102 if(UseAngularGeneratorFlag() && nullptr == GetAngularDistribution()) { 111 SetAngularDistribution(new G4DeltaAngle( 103 SetAngularDistribution(new G4DeltaAngle()); 112 } 104 } 113 } 105 } 114 if(nullptr == lsdata) { << 106 if(IsMaster()) { 115 G4AutoLock l(&ionXSMutex); << 116 if(nullptr == lsdata) { 107 if(nullptr == lsdata) { 117 isFirst = true; << 118 lsdata = new G4LindhardSorensenData(); 108 lsdata = new G4LindhardSorensenData(); >> 109 } >> 110 if(nullptr == fIonData) { 119 fIonData = new G4IonICRU73Data(); 111 fIonData = new G4IonICRU73Data(); 120 } 112 } 121 l.unlock(); << 122 } << 123 if(isFirst) { << 124 fIonData->Initialise(); 113 fIonData->Initialise(); 125 } 114 } 126 } 115 } 127 116 128 //....oooOO0OOooo........oooOO0OOooo........oo 117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 129 118 130 G4double 119 G4double 131 G4LindhardSorensenIonModel::GetChargeSquareRat 120 G4LindhardSorensenIonModel::GetChargeSquareRatio(const G4ParticleDefinition* p, 132 121 const G4Material* mat, 133 122 G4double kinEnergy) 134 { 123 { 135 chargeSquare = corr->EffectiveChargeSquareRa << 124 chargeSquare = corr->EffectiveChargeSquareRatio(p,mat,kinEnergy); 136 return chargeSquare; 125 return chargeSquare; 137 } 126 } 138 127 139 //....oooOO0OOooo........oooOO0OOooo........oo 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 140 129 141 G4double 130 G4double 142 G4LindhardSorensenIonModel::GetParticleCharge( 131 G4LindhardSorensenIonModel::GetParticleCharge(const G4ParticleDefinition* p, 143 132 const G4Material* mat, 144 133 G4double kinEnergy) 145 { 134 { 146 return corr->GetParticleCharge(p,mat,kinEner 135 return corr->GetParticleCharge(p,mat,kinEnergy); 147 } 136 } 148 137 149 //....oooOO0OOooo........oooOO0OOooo........oo 138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 150 139 151 void G4LindhardSorensenIonModel::SetupParamete 140 void G4LindhardSorensenIonModel::SetupParameters() 152 { 141 { 153 mass = particle->GetPDGMass(); 142 mass = particle->GetPDGMass(); 154 spin = particle->GetPDGSpin(); 143 spin = particle->GetPDGSpin(); 155 charge = particle->GetPDGCharge()*inveplus; 144 charge = particle->GetPDGCharge()*inveplus; 156 Zin = G4lrint(std::abs(charge)); 145 Zin = G4lrint(std::abs(charge)); 157 chargeSquare = charge*charge; 146 chargeSquare = charge*charge; 158 eRatio = CLHEP::electron_mass_c2/mass; 147 eRatio = CLHEP::electron_mass_c2/mass; 159 pRatio = CLHEP::proton_mass_c2/mass; 148 pRatio = CLHEP::proton_mass_c2/mass; 160 const G4double aMag = 149 const G4double aMag = 161 1./(0.5*CLHEP::eplus*CLHEP::hbar_Planck*CL 150 1./(0.5*CLHEP::eplus*CLHEP::hbar_Planck*CLHEP::c_squared); 162 G4double magmom = particle->GetPDGMagneticMo 151 G4double magmom = particle->GetPDGMagneticMoment()*mass*aMag; 163 magMoment2 = magmom*magmom - 1.0; 152 magMoment2 = magmom*magmom - 1.0; 164 G4double x = 0.8426*CLHEP::GeV; 153 G4double x = 0.8426*CLHEP::GeV; 165 if(spin == 0.0 && mass < GeV) { x = 0.736*CL 154 if(spin == 0.0 && mass < GeV) { x = 0.736*CLHEP::GeV; } 166 else if (Zin > 1) { x /= nist->GetA27(Zin); 155 else if (Zin > 1) { x /= nist->GetA27(Zin); } 167 156 168 formfact = 2.0*CLHEP::electron_mass_c2/(x*x) 157 formfact = 2.0*CLHEP::electron_mass_c2/(x*x); 169 tlimit = 2.0/formfact; 158 tlimit = 2.0/formfact; 170 } 159 } 171 160 172 //....oooOO0OOooo........oooOO0OOooo........oo 161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 173 162 174 G4double G4LindhardSorensenIonModel::MinEnergy 163 G4double G4LindhardSorensenIonModel::MinEnergyCut(const G4ParticleDefinition*, 175 const 164 const G4MaterialCutsCouple* couple) 176 { 165 { 177 return couple->GetMaterial()->GetIonisation( 166 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy(); 178 } 167 } 179 168 180 //....oooOO0OOooo........oooOO0OOooo........oo 169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 181 170 182 G4double 171 G4double 183 G4LindhardSorensenIonModel::ComputeCrossSectio 172 G4LindhardSorensenIonModel::ComputeCrossSectionPerElectron( 184 173 const G4ParticleDefinition* p, 185 174 G4double kinEnergy, 186 175 G4double cutEnergy, 187 176 G4double maxKinEnergy) 188 { 177 { 189 // take into account formfactor 178 // take into account formfactor 190 G4double tmax = MaxSecondaryEnergy(p, kinEne 179 G4double tmax = MaxSecondaryEnergy(p, kinEnergy); 191 G4double emax = std::min(tmax, maxKinEnergy) 180 G4double emax = std::min(tmax, maxKinEnergy); 192 G4double escaled = kinEnergy*pRatio; 181 G4double escaled = kinEnergy*pRatio; 193 G4double cross = (escaled <= fElimit) 182 G4double cross = (escaled <= fElimit) 194 ? fBraggModel->ComputeCrossSectionPerElect 183 ? fBraggModel->ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,emax) 195 : fBBModel->ComputeCrossSectionPerElectron 184 : fBBModel->ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,emax); 196 // G4cout << "LS: e= " << kinEnergy << " tmi 185 // G4cout << "LS: e= " << kinEnergy << " tmin= " << cutEnergy 197 // << " tmax= " << maxEnergy << " cro 186 // << " tmax= " << maxEnergy << " cross= " << cross << G4endl; 198 return cross; 187 return cross; 199 } 188 } 200 189 201 //....oooOO0OOooo........oooOO0OOooo........oo 190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 202 191 203 G4double G4LindhardSorensenIonModel::ComputeCr 192 G4double G4LindhardSorensenIonModel::ComputeCrossSectionPerAtom( 204 con 193 const G4ParticleDefinition* p, 205 194 G4double kineticEnergy, 206 195 G4double Z, G4double, 207 196 G4double cutEnergy, 208 197 G4double maxEnergy) 209 { 198 { 210 return Z*ComputeCrossSectionPerElectron(p,ki 199 return Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy); 211 } 200 } 212 201 213 //....oooOO0OOooo........oooOO0OOooo........oo 202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 214 203 215 G4double G4LindhardSorensenIonModel::CrossSect 204 G4double G4LindhardSorensenIonModel::CrossSectionPerVolume( 216 con 205 const G4Material* material, 217 con 206 const G4ParticleDefinition* p, 218 207 G4double kineticEnergy, 219 208 G4double cutEnergy, 220 209 G4double maxEnergy) 221 { 210 { 222 return material->GetElectronDensity() 211 return material->GetElectronDensity() 223 *ComputeCrossSectionPerElectron(p,kineticE 212 *ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy); 224 } 213 } 225 214 226 //....oooOO0OOooo........oooOO0OOooo........oo 215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 227 216 228 G4double 217 G4double 229 G4LindhardSorensenIonModel::ComputeDEDXPerVolu 218 G4LindhardSorensenIonModel::ComputeDEDXPerVolume(const G4Material* mat, 230 219 const G4ParticleDefinition* p, 231 220 G4double kinEnergy, 232 221 G4double cut) 233 { 222 { 234 // formfactor is taken into account in Corre 223 // formfactor is taken into account in CorrectionsAlongStep(..) 235 G4double tmax = MaxSecondaryEnergy(p, k 224 G4double tmax = MaxSecondaryEnergy(p, kinEnergy); 236 G4double cutEnergy = std::min(std::min(cut,t 225 G4double cutEnergy = std::min(std::min(cut,tmax), tlimit); 237 226 238 G4double escaled = kinEnergy*pRatio; 227 G4double escaled = kinEnergy*pRatio; 239 G4double dedx = (escaled <= fElimit) 228 G4double dedx = (escaled <= fElimit) 240 ? fBraggModel->ComputeDEDXPerVolume(mat, p 229 ? fBraggModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cutEnergy) 241 : fBBModel->ComputeDEDXPerVolume(mat, p, k 230 : fBBModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cutEnergy); 242 231 243 //G4cout << "E(MeV)=" << kinEnergy/MeV << " 232 //G4cout << "E(MeV)=" << kinEnergy/MeV << " dedx=" << dedx 244 // << " " << mat->GetName() << " Ecut(MeV) << 233 // << " " << material->GetName() << " Ecut(MeV)=" << cutEnergy << G4endl; 245 return dedx; 234 return dedx; 246 } 235 } 247 236 248 //....oooOO0OOooo........oooOO0OOooo........oo 237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 249 238 250 void G4LindhardSorensenIonModel::CorrectionsAl 239 void G4LindhardSorensenIonModel::CorrectionsAlongStep( 251 const G4Mater 240 const G4MaterialCutsCouple* couple, 252 const G4Dynam 241 const G4DynamicParticle* dp, 253 const G4doubl 242 const G4double& length, 254 G4doubl 243 G4double& eloss) 255 { 244 { 256 // no correction at the last step 245 // no correction at the last step 257 const G4double preKinEnergy = dp->GetKinetic 246 const G4double preKinEnergy = dp->GetKineticEnergy(); 258 if(eloss >= preKinEnergy) { return; } 247 if(eloss >= preKinEnergy) { return; } 259 248 260 const G4ParticleDefinition* p = dp->GetDefin 249 const G4ParticleDefinition* p = dp->GetDefinition(); 261 SetParticle(p); 250 SetParticle(p); 262 const G4Material* mat = couple->GetMaterial( 251 const G4Material* mat = couple->GetMaterial(); 263 const G4double eDensity = mat->GetElectronDe 252 const G4double eDensity = mat->GetElectronDensity(); 264 const G4double e = std::max(preKinEnergy - e 253 const G4double e = std::max(preKinEnergy - eloss*0.5, preKinEnergy*0.5); 265 const G4double tmax = MaxSecondaryEnergy(p, 254 const G4double tmax = MaxSecondaryEnergy(p, e); 266 const G4double escaled = e*pRatio; 255 const G4double escaled = e*pRatio; 267 const G4double tau = e/mass; 256 const G4double tau = e/mass; >> 257 >> 258 const G4double q20 = corr->EffectiveChargeSquareRatio(p, mat, preKinEnergy); 268 const G4double q2 = corr->EffectiveChargeSqu 259 const G4double q2 = corr->EffectiveChargeSquareRatio(p, mat, e); 269 const G4int Z = p->GetAtomicNumber(); 260 const G4int Z = p->GetAtomicNumber(); 270 261 271 G4double res = 0.0; 262 G4double res = 0.0; 272 if(escaled <= fElimit) { 263 if(escaled <= fElimit) { 273 // data from ICRU73 or ICRU90 264 // data from ICRU73 or ICRU90 274 if(Z > 2 && Z <= 80) { 265 if(Z > 2 && Z <= 80) { 275 res = fIonData->GetDEDX(mat, Z, escaled, 266 res = fIonData->GetDEDX(mat, Z, escaled, G4Log(escaled)); 276 /* 267 /* 277 G4cout << " GetDEDX for Z=" << Z << " in " << 268 G4cout << "GetDEDX for Z=" << Z << " in " << mat->GetName() 278 << " Escaled=" << escaled << " E=" 269 << " Escaled=" << escaled << " E=" 279 << e << " dEdx=" << res << G4endl; 270 << e << " dEdx=" << res << G4endl; 280 */ 271 */ 281 } 272 } 282 if(res > 0.0) { 273 if(res > 0.0) { 283 auto pcuts = couple->GetProductionCuts() 274 auto pcuts = couple->GetProductionCuts(); 284 G4double cut = (nullptr == pcuts) ? tmax 275 G4double cut = (nullptr == pcuts) ? tmax : pcuts->GetProductionCut(1); 285 if(cut < tmax) { 276 if(cut < tmax) { 286 const G4double x = cut/tmax; 277 const G4double x = cut/tmax; 287 res += (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau 278 res += (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau * (tau + 2.0)) + 1.0 - x) 288 *q2*CLHEP::twopi_mc2_rcl2*eDensity; 279 *q2*CLHEP::twopi_mc2_rcl2*eDensity; 289 } 280 } 290 res *= length; 281 res *= length; 291 } else { 282 } else { 292 // simplified correction 283 // simplified correction 293 res = eloss*q2/chargeSquare; << 284 res = eloss*q2/q20; 294 } 285 } 295 } else { 286 } else { 296 // Lindhard-Sorensen model 287 // Lindhard-Sorensen model 297 const G4double gam = tau + 1.0; 288 const G4double gam = tau + 1.0; 298 const G4double beta2 = tau * (tau+2.0)/(ga 289 const G4double beta2 = tau * (tau+2.0)/(gam*gam); 299 G4double deltaL0 = 2.0*corr->BarkasCorrect 290 G4double deltaL0 = 2.0*corr->BarkasCorrection(p, mat, e)*(charge-1.)/charge; 300 G4double deltaL = lsdata->GetDeltaL(Zin, g 291 G4double deltaL = lsdata->GetDeltaL(Zin, gam); 301 292 302 res = eloss + 293 res = eloss + 303 CLHEP::twopi_mc2_rcl2*q2*eDensity*(delta 294 CLHEP::twopi_mc2_rcl2*q2*eDensity*(deltaL+deltaL0)*length/beta2; 304 /* << 295 /* 305 G4cout << " E(GeV)=" << preKinEnergy/GeV << 296 G4cout << "G4LindhardSorensenIonModel::CorrectionsAlongStep: E(GeV)= " 306 << " L= " << eloss*beta2/(twopi_mc2_rcl2* << 297 << preKinEnergy/GeV << " eloss(MeV)= " << eloss 307 << " dL0= " << deltaL0 << 298 << " L= " << eloss*beta2/(twopi_mc2_rcl2*chargeSquare*eDensity*length) 308 << " dL= " << deltaL << " dE(MeV)=" << re << 299 << " dL0= " << deltaL0 309 */ << 300 << " dL= " << deltaL << G4endl; >> 301 */ 310 } 302 } 311 if(res > preKinEnergy || 2*res < eloss) { re << 303 if(res > preKinEnergy) { res = preKinEnergy; } >> 304 else if(res < 0.0) { res = eloss; } 312 /* 305 /* 313 G4cout << " G4LindhardSorensenIonModel::Cor << 306 G4cout << "G4LindhardSorensenIonModel::CorrectionsAlongStep: E(GeV)=" 314 << preKinEnergy/GeV << " eloss(MeV)=" 307 << preKinEnergy/GeV << " eloss(MeV)=" << eloss 315 << " res(MeV)=" << res << G4endl; 308 << " res(MeV)=" << res << G4endl; 316 */ 309 */ 317 eloss = res; 310 eloss = res; 318 } 311 } 319 312 320 //....oooOO0OOooo........oooOO0OOooo........oo 313 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 321 314 322 void G4LindhardSorensenIonModel::SampleSeconda 315 void G4LindhardSorensenIonModel::SampleSecondaries( 323 std::vector<G4DynamicParticle*>* v 316 std::vector<G4DynamicParticle*>* vdp, 324 cons 317 const G4MaterialCutsCouple* couple, 325 cons 318 const G4DynamicParticle* dp, 326 G4do 319 G4double cut, 327 G4do 320 G4double maxEnergy) 328 { 321 { 329 G4double kineticEnergy = dp->GetKineticEnerg 322 G4double kineticEnergy = dp->GetKineticEnergy(); 330 // take into account formfactor 323 // take into account formfactor 331 const G4double tmax = MaxSecondaryEnergy(dp- << 324 G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy); 332 const G4double minKinEnergy = std::min(cut, << 325 G4double minKinEnergy = std::min(cut, tmax); 333 const G4double maxKinEnergy = std::min(maxEn << 326 G4double maxKinEnergy = std::min(maxEnergy,tmax); 334 if(minKinEnergy >= maxKinEnergy) { return; } 327 if(minKinEnergy >= maxKinEnergy) { return; } 335 328 336 //G4cout << "G4LindhardSorensenIonModel::Sam 329 //G4cout << "G4LindhardSorensenIonModel::SampleSecondaries Emin= " 337 // << minKinEnergy << " Emax= " << maxKinEne 330 // << minKinEnergy << " Emax= " << maxKinEnergy << G4endl; 338 331 339 G4double totEnergy = kineticEnergy + mass; << 332 G4double totEnergy = kineticEnergy + mass; 340 G4double etot2 = totEnergy*totEnergy; << 333 G4double etot2 = totEnergy*totEnergy; 341 G4double beta2 = kineticEnergy*(kineticEnerg << 334 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2; 342 335 343 G4double deltaKinEnergy, f; 336 G4double deltaKinEnergy, f; 344 G4double f1 = 0.0; 337 G4double f1 = 0.0; 345 G4double fmax = 1.0; 338 G4double fmax = 1.0; 346 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy* 339 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; } 347 340 348 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra 341 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine(); 349 G4double rndm[2]; 342 G4double rndm[2]; 350 343 351 // sampling without nuclear size effect 344 // sampling without nuclear size effect 352 do { 345 do { 353 rndmEngineMod->flatArray(2, rndm); 346 rndmEngineMod->flatArray(2, rndm); 354 deltaKinEnergy = minKinEnergy*maxKinEnergy 347 deltaKinEnergy = minKinEnergy*maxKinEnergy 355 /(minKinEnergy*(1.0 - rndm 348 /(minKinEnergy*(1.0 - rndm[0]) + maxKinEnergy*rndm[0]); 356 349 357 f = 1.0 - beta2*deltaKinEnergy/tmax; 350 f = 1.0 - beta2*deltaKinEnergy/tmax; 358 if( 0.0 < spin ) { 351 if( 0.0 < spin ) { 359 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/e 352 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2; 360 f += f1; 353 f += f1; 361 } 354 } 362 355 363 // Loop checking, 03-Aug-2015, Vladimir Iv 356 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 364 } while( fmax*rndm[1] > f); 357 } while( fmax*rndm[1] > f); 365 358 366 // projectile formfactor - suppresion of hig 359 // projectile formfactor - suppresion of high energy 367 // delta-electron production at high energy 360 // delta-electron production at high energy 368 361 369 G4double x = formfact*deltaKinEnergy; 362 G4double x = formfact*deltaKinEnergy; 370 if(x > 1.e-6) { 363 if(x > 1.e-6) { 371 364 372 G4double x1 = 1.0 + x; 365 G4double x1 = 1.0 + x; 373 G4double grej = 1.0/(x1*x1); 366 G4double grej = 1.0/(x1*x1); 374 if( 0.0 < spin ) { 367 if( 0.0 < spin ) { 375 G4double x2 = 0.5*electron_mass_c2*delta 368 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass); 376 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1 369 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2)); 377 } 370 } 378 if(grej > 1.1) { 371 if(grej > 1.1) { 379 G4cout << "### G4LindhardSorensenIonMode 372 G4cout << "### G4LindhardSorensenIonModel WARNING: grej= " << grej 380 << " " << dp->GetDefinition()->G 373 << " " << dp->GetDefinition()->GetParticleName() 381 << " Ekin(MeV)= " << kineticEner 374 << " Ekin(MeV)= " << kineticEnergy 382 << " delEkin(MeV)= " << deltaKinE 375 << " delEkin(MeV)= " << deltaKinEnergy 383 << G4endl; 376 << G4endl; 384 } 377 } 385 if(rndmEngineMod->flat() > grej) { return; 378 if(rndmEngineMod->flat() > grej) { return; } 386 } 379 } 387 380 388 G4ThreeVector deltaDirection; 381 G4ThreeVector deltaDirection; 389 382 390 if(UseAngularGeneratorFlag()) { 383 if(UseAngularGeneratorFlag()) { 391 384 392 const G4Material* mat = couple->GetMateri 385 const G4Material* mat = couple->GetMaterial(); 393 G4int Z = SelectRandomAtomNumber(mat); 386 G4int Z = SelectRandomAtomNumber(mat); 394 387 395 deltaDirection = 388 deltaDirection = 396 GetAngularDistribution()->SampleDirectio 389 GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat); 397 390 398 } else { 391 } else { 399 392 400 G4double deltaMomentum = 393 G4double deltaMomentum = 401 std::sqrt(deltaKinEnergy * (deltaKinEner 394 std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2)); 402 G4double cost = deltaKinEnergy * (totEnerg 395 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) / 403 (deltaMomentum * dp->GetTotalMomentum()) 396 (deltaMomentum * dp->GetTotalMomentum()); 404 cost = std::min(cost, 1.0); 397 cost = std::min(cost, 1.0); 405 G4double sint = std::sqrt((1.0 - cost)*(1. 398 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost)); 406 399 407 G4double phi = CLHEP::twopi*rndmEngineMod- 400 G4double phi = CLHEP::twopi*rndmEngineMod->flat(); 408 401 409 deltaDirection.set(sint*std::cos(phi),sint 402 deltaDirection.set(sint*std::cos(phi),sint*std::sin(phi), cost) ; 410 deltaDirection.rotateUz(dp->GetMomentumDir 403 deltaDirection.rotateUz(dp->GetMomentumDirection()); 411 } 404 } 412 /* 405 /* 413 G4cout << "### G4LindhardSorensenIonModel 406 G4cout << "### G4LindhardSorensenIonModel " 414 << dp->GetDefinition()->GetParticle 407 << dp->GetDefinition()->GetParticleName() 415 << " Ekin(MeV)= " << kineticEnergy 408 << " Ekin(MeV)= " << kineticEnergy 416 << " delEkin(MeV)= " << deltaKinEne 409 << " delEkin(MeV)= " << deltaKinEnergy 417 << " tmin(MeV)= " << minKinEnergy 410 << " tmin(MeV)= " << minKinEnergy 418 << " tmax(MeV)= " << maxKinEnergy 411 << " tmax(MeV)= " << maxKinEnergy 419 << " dir= " << dp->GetMomentumDirec 412 << " dir= " << dp->GetMomentumDirection() 420 << " dirDelta= " << deltaDirection 413 << " dirDelta= " << deltaDirection 421 << G4endl; 414 << G4endl; 422 */ 415 */ 423 // create G4DynamicParticle object for delta 416 // create G4DynamicParticle object for delta ray 424 auto delta = new G4DynamicParticle(theElectr 417 auto delta = new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy); 425 418 426 vdp->push_back(delta); 419 vdp->push_back(delta); 427 420 428 // Change kinematics of primary particle 421 // Change kinematics of primary particle 429 kineticEnergy -= deltaKinEnergy; 422 kineticEnergy -= deltaKinEnergy; 430 G4ThreeVector finalP = dp->GetMomentum() - d 423 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum(); 431 finalP = finalP.unit(); << 424 finalP = finalP.unit(); 432 425 433 fParticleChange->SetProposedKineticEnergy(ki 426 fParticleChange->SetProposedKineticEnergy(kineticEnergy); 434 fParticleChange->SetProposedMomentumDirectio 427 fParticleChange->SetProposedMomentumDirection(finalP); 435 } 428 } 436 429 437 //....oooOO0OOooo........oooOO0OOooo........oo 430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 438 431 439 G4double 432 G4double 440 G4LindhardSorensenIonModel::MaxSecondaryEnergy 433 G4LindhardSorensenIonModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd, 441 434 G4double kinEnergy) 442 { 435 { 443 // here particle type is checked for any met 436 // here particle type is checked for any method 444 SetParticle(pd); 437 SetParticle(pd); 445 G4double tau = kinEnergy/mass; 438 G4double tau = kinEnergy/mass; 446 return 2.0*CLHEP::electron_mass_c2*tau*(tau 439 return 2.0*CLHEP::electron_mass_c2*tau*(tau + 2.) / 447 (1. + 2.0*(tau + 1.)*eRatio + eRatio*eRati 440 (1. + 2.0*(tau + 1.)*eRatio + eRatio*eRatio); 448 } 441 } 449 442 450 //....oooOO0OOooo........oooOO0OOooo........oo 443 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 451 444