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 // $Id: G4eCoulombScatteringModel.cc,v 1.59 2008/10/22 18:39:29 vnivanch Exp $ >> 27 // GEANT4 tag $Name: geant4-09-02-patch-04 $ 26 // 28 // 27 // ------------------------------------------- 29 // ------------------------------------------------------------------- 28 // 30 // 29 // GEANT4 Class file 31 // GEANT4 Class file 30 // 32 // 31 // 33 // 32 // File name: G4eCoulombScatteringModel 34 // File name: G4eCoulombScatteringModel 33 // 35 // 34 // Author: Vladimir Ivanchenko 36 // Author: Vladimir Ivanchenko 35 // 37 // 36 // Creation date: 22.08.2005 38 // Creation date: 22.08.2005 37 // 39 // 38 // Modifications: V.Ivanchenko << 40 // Modifications: 39 // << 41 // 01.08.06 V.Ivanchenko extend upper limit of table to TeV and review the 40 // << 42 // logic of building - only elements from G4ElementTable >> 43 // 08.08.06 V.Ivanchenko build internal table in ekin scale, introduce faclim >> 44 // 19.08.06 V.Ivanchenko add inline function ScreeningParameter >> 45 // 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e- >> 46 // 09.06.08 V.Ivanchenko add SelectIsotope and sampling of the recoil ion 41 // 47 // 42 // Class Description: 48 // Class Description: 43 // 49 // 44 // ------------------------------------------- 50 // ------------------------------------------------------------------- 45 // 51 // 46 //....oooOO0OOooo........oooOO0OOooo........oo 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 47 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 48 54 49 #include "G4eCoulombScatteringModel.hh" 55 #include "G4eCoulombScatteringModel.hh" 50 #include "G4PhysicalConstants.hh" << 51 #include "G4SystemOfUnits.hh" << 52 #include "Randomize.hh" 56 #include "Randomize.hh" 53 #include "G4DataVector.hh" 57 #include "G4DataVector.hh" 54 #include "G4ElementTable.hh" 58 #include "G4ElementTable.hh" >> 59 #include "G4PhysicsLogVector.hh" 55 #include "G4ParticleChangeForGamma.hh" 60 #include "G4ParticleChangeForGamma.hh" >> 61 #include "G4Electron.hh" >> 62 #include "G4Positron.hh" 56 #include "G4Proton.hh" 63 #include "G4Proton.hh" 57 #include "G4ParticleTable.hh" 64 #include "G4ParticleTable.hh" 58 #include "G4IonTable.hh" << 59 #include "G4ProductionCutsTable.hh" << 60 #include "G4NucleiProperties.hh" << 61 #include "G4Pow.hh" << 62 #include "G4NistManager.hh" << 63 65 64 //....oooOO0OOooo........oooOO0OOooo........oo 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 65 67 66 using namespace std; 68 using namespace std; 67 69 68 G4eCoulombScatteringModel::G4eCoulombScatterin << 70 G4eCoulombScatteringModel::G4eCoulombScatteringModel(const G4String& nam) 69 : G4VEmModel("eCoulombScattering"), isCombin << 71 : G4VEmModel(nam), >> 72 cosThetaMin(1.0), >> 73 cosThetaMax(-1.0), >> 74 q2Limit(TeV*TeV), >> 75 alpha2(fine_structure_const*fine_structure_const), >> 76 faclim(100.0), >> 77 isInitialised(false) 70 { 78 { 71 fNistManager = G4NistManager::Instance(); 79 fNistManager = G4NistManager::Instance(); 72 theIonTable = G4ParticleTable::GetParticleT << 80 theParticleTable = G4ParticleTable::GetParticleTable(); 73 theProton = G4Proton::Proton(); << 81 theElectron = G4Electron::Electron(); 74 << 82 thePositron = G4Positron::Positron(); 75 wokvi = new G4WentzelOKandVIxSection(isCombi << 83 theProton = G4Proton::Proton(); 76 << 84 currentMaterial = 0; 77 mass = CLHEP::proton_mass_c2; << 85 currentElement = 0; >> 86 a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885); >> 87 G4double p0 = electron_mass_c2*classic_electr_radius; >> 88 coeff = twopi*p0*p0; >> 89 constn = 6.937e-6/(MeV*MeV); >> 90 tkin = targetZ = mom2 = DBL_MIN; >> 91 elecXSection = nucXSection = 0.0; >> 92 recoilThreshold = DBL_MAX; >> 93 ecut = DBL_MAX; >> 94 particle = 0; >> 95 currentCouple = 0; >> 96 for(size_t j=0; j<100; j++) { >> 97 FF[j] = 0.0; >> 98 } 78 } 99 } 79 100 80 //....oooOO0OOooo........oooOO0OOooo........oo 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 81 102 82 G4eCoulombScatteringModel::~G4eCoulombScatteri 103 G4eCoulombScatteringModel::~G4eCoulombScatteringModel() 83 { << 104 {} 84 delete wokvi; << 85 } << 86 105 87 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 88 107 89 void G4eCoulombScatteringModel::Initialise(con << 108 void G4eCoulombScatteringModel::Initialise(const G4ParticleDefinition* p, 90 const G4DataVector& cuts) 109 const G4DataVector& cuts) 91 { 110 { 92 SetupParticle(part); << 111 SetupParticle(p); 93 currentCouple = nullptr; << 112 currentCouple = 0; 94 << 113 elecXSection = nucXSection = 0.0; 95 G4double tet = PolarAngleLimit(); << 114 tkin = targetZ = mom2 = DBL_MIN; 96 << 115 ecut = etag = DBL_MAX; 97 // defined theta limit between single and mu << 116 cosThetaMin = cos(PolarAngleLimit()); 98 if(isCombined) { << 117 currentCuts = &cuts; 99 if(tet >= CLHEP::pi) { cosThetaMin = -1.0; << 118 //G4cout << "!!! G4eCoulombScatteringModel::Initialise for " 100 else if(tet > 0.0) { cosThetaMin = std::co << 119 // << p->GetParticleName() << " cos(TetMin)= " << cosThetaMin 101 << 120 // << " cos(TetMax)= " << cosThetaMax <<G4endl; 102 // single scattering without multiple << 121 if(!isInitialised) { 103 } else if(tet > 0.0) { << 122 isInitialised = true; 104 cosThetaMin = std::cos(std::min(tet, CLHEP << 123 105 } << 124 if(pParticleChange) 106 << 125 fParticleChange = 107 wokvi->Initialise(part, cosThetaMin); << 126 reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange); 108 pCuts = &cuts; << 127 else 109 /* << 128 fParticleChange = new G4ParticleChangeForGamma(); 110 G4cout << "G4eCoulombScatteringModel::Initia << 111 << part->GetParticleName() << " 1-cos(Tet << 112 << " 1-cos(TetMax)= " << 1. - cosThetaMax << 113 G4cout << "cut[0]= " << (*pCuts)[0] << G4end << 114 */ << 115 if(nullptr == fParticleChange) { << 116 fParticleChange = GetParticleChangeForGamm << 117 } 129 } 118 if(IsMaster() && mass < GeV && part->GetPart << 130 if(mass < GeV && particle->GetParticleType() != "nucleus") { 119 InitialiseElementSelectors(part, cuts); << 131 InitialiseElementSelectors(p,cuts); 120 } 132 } 121 } 133 } 122 134 123 //....oooOO0OOooo........oooOO0OOooo........oo 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 124 136 125 void G4eCoulombScatteringModel::InitialiseLoca << 137 void G4eCoulombScatteringModel::ComputeMaxElectronScattering(G4double cutEnergy) 126 G4VEmModel* masterModel) << 127 { 138 { 128 SetElementSelectors(masterModel->GetElementS << 139 ecut = cutEnergy; 129 } << 140 G4double tmax = tkin; 130 << 141 cosTetMaxElec = 1.0; 131 //....oooOO0OOooo........oooOO0OOooo........oo << 142 if(mass > MeV) { 132 << 143 G4double ratio = electron_mass_c2/mass; 133 G4double << 144 G4double tau = tkin/mass; 134 G4eCoulombScatteringModel::MinPrimaryEnergy(co << 145 tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/ 135 const G4ParticleDefinition* part << 146 (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio); 136 G4double) << 147 cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2; 137 { << 148 } else { 138 SetupParticle(part); << 139 149 140 // define cut using cuts for proton << 150 if(particle == theElectron) tmax *= 0.5; 141 G4double cut = << 151 G4double t = std::min(cutEnergy, tmax); 142 std::max(recoilThreshold, (*pCuts)[Current << 152 G4double mom21 = t*(t + 2.0*electron_mass_c2); 143 << 153 G4double t1 = tkin - t; 144 // find out lightest element << 154 //G4cout << "tkin= " << tkin << " t= " << t << " t1= " << t1 << G4endl; 145 const G4ElementVector* theElementVector = ma << 155 if(t1 > 0.0) { 146 std::size_t nelm = material->GetNumberOfElem << 156 G4double mom22 = t1*(t1 + 2.0*mass); 147 << 157 G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22); 148 // select lightest element << 158 //G4cout << "ctm= " << ctm << G4endl; 149 G4int Z = 300; << 159 if(ctm < 1.0) cosTetMaxElec = ctm; 150 for (std::size_t j=0; j<nelm; ++j) { << 160 } 151 Z = std::min(Z,(*theElementVector)[j]->Get << 152 } 161 } 153 G4int A = G4lrint(fNistManager->GetAtomicMas << 154 G4double targetMass = G4NucleiProperties::Ge << 155 G4double t = std::max(cut, 0.5*(cut + sqrt(2 << 156 << 157 return t; << 158 } 162 } 159 163 160 //....oooOO0OOooo........oooOO0OOooo........oo 164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 161 165 162 G4double G4eCoulombScatteringModel::ComputeCro 166 G4double G4eCoulombScatteringModel::ComputeCrossSectionPerAtom( 163 const G4ParticleDefinition* p, 167 const G4ParticleDefinition* p, 164 G4double kinEnergy, 168 G4double kinEnergy, 165 G4double Z, G4double, 169 G4double Z, G4double, 166 G4double cutEnergy, G4double) 170 G4double cutEnergy, G4double) 167 { 171 { 168 /* << 172 //G4cout << "### G4eCoulombScatteringModel::ComputeCrossSectionPerAtom for " 169 G4cout << "### G4eCoulombScatteringModel::Co << 173 // << p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV << G4endl; 170 << p->GetParticleName()<<" Z= "<<Z<<" e(MeV << 174 G4double xsec = 0.0; 171 << G4endl; << 175 SetupParticle(p); >> 176 G4double ekin = std::max(lowEnergyLimit, kinEnergy); >> 177 SetupKinematic(ekin, cutEnergy); >> 178 if(cosTetMaxNuc < cosTetMinNuc) { >> 179 SetupTarget(Z, ekin); >> 180 xsec = CrossSectionPerAtom(); >> 181 } >> 182 /* >> 183 G4cout << "e(MeV)= " << ekin/MeV << "cosTetMinNuc= " << cosTetMinNuc >> 184 << " cosTetMaxNuc= " << cosTetMaxNuc >> 185 << " cosTetMaxElec= " << cosTetMaxElec >> 186 << " screenZ= " << screenZ >> 187 << " formfactA= " << formfactA >> 188 << " cosTetMaxHad= " << cosTetMaxHad << G4endl; 172 */ 189 */ 173 G4double cross = 0.0; << 190 return xsec; 174 elecRatio = 0.0; << 191 } 175 if(p != particle) { SetupParticle(p); } << 192 176 << 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 177 // cross section is set to zero to avoid pro << 194 178 if(kinEnergy <= 0.0) { return cross; } << 195 G4double G4eCoulombScatteringModel::CrossSectionPerAtom() 179 DefineMaterial(CurrentCouple()); << 196 { 180 G4double costmin = wokvi->SetupKinematic(kin << 197 // This method needs initialisation before be called 181 << 198 182 //G4cout << "cosThetaMax= "<<cosThetaMax<<" << 199 G4double fac = coeff*targetZ*chargeSquare*invbeta2/mom2; 183 << 200 elecXSection = 0.0; 184 if(cosThetaMax < costmin) { << 201 nucXSection = 0.0; 185 G4int iz = G4lrint(Z); << 202 186 G4double cut = (0.0 < fixedCut) ? fixedCut << 203 G4double x = 1.0 - cosTetMinNuc; 187 costmin = wokvi->SetupTarget(iz, cut); << 204 G4double x1 = x + screenZ; 188 //G4cout << "SetupTarget: Z= " << iz << " << 205 189 // << costmin << G4endl; << 206 if(cosTetMaxElec2 < cosTetMinNuc) { 190 G4double costmax = (1 == iz && particle == << 207 elecXSection = fac*(cosTetMinNuc - cosTetMaxElec2)/ 191 ? 0.0 : cosThetaMax; << 208 (x1*(1.0 - cosTetMaxElec2 + screenZ)); 192 if(costmin > costmax) { << 209 nucXSection = elecXSection; 193 cross = wokvi->ComputeNuclearCrossSectio << 210 } 194 + wokvi->ComputeElectronCrossSection(c << 211 >> 212 //G4cout << "XS tkin(MeV)= " << tkin<<" xs= " <<nucXSection >> 213 // << " costmax= " << cosTetMaxNuc2 >> 214 // << " costmin= " << cosTetMinNuc << " Z= " << targetZ <<G4endl; >> 215 if(cosTetMaxNuc2 < cosTetMinNuc) { >> 216 G4double s = screenZ*formfactA; >> 217 G4double z1 = 1.0 - cosTetMaxNuc2 + screenZ; >> 218 G4double d = (1.0 - s)/formfactA; >> 219 //G4cout <<"x1= "<<x1<<" z1= " <<z1<<" s= "<<s << " d= " <<d <<G4endl; >> 220 if(d < 0.2*x1) { >> 221 G4double x2 = x1*x1; >> 222 G4double z2 = z1*z1; >> 223 x = (1.0/(x1*x2) - 1.0/(z1*z2) - d*1.5*(1.0/(x2*x2) - 1.0/(z2*z2)))/ >> 224 (3.0*formfactA*formfactA); >> 225 } else { >> 226 G4double x2 = x1 + d; >> 227 G4double z2 = z1 + d; >> 228 x = (1.0 + 2.0*s)*((cosTetMinNuc - cosTetMaxNuc2)*(1.0/(x1*z1) + 1.0/(x2*z2)) - >> 229 2.0*log(z1*x2/(z2*x1))/d); 195 } 230 } 196 /* << 231 nucXSection += fac*targetZ*x; 197 if(p->GetParticleName() == "e-") << 198 G4cout << "Z= " << Z << " e(MeV)= " << kin << 199 << " cross(b)= " << cross/barn << " 1-cos << 200 << " 1-costmax= " << 1-costmax << 201 << " 1-cosThetaMax= " << 1-cosThetaMax << 202 << " " << currentMaterial->GetName() << 203 << G4endl; << 204 */ << 205 } 232 } 206 //G4cout << "====== cross= " << cross << G4e << 233 207 return cross; << 234 //G4cout<<" cross(bn)= "<<nucXSection/barn<<" xsElec(bn)= "<<elecXSection/barn >> 235 // << " Asc= " << screenZ << G4endl; >> 236 >> 237 return nucXSection; 208 } 238 } 209 239 210 //....oooOO0OOooo........oooOO0OOooo........oo 240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 211 241 212 void G4eCoulombScatteringModel::SampleSecondar 242 void G4eCoulombScatteringModel::SampleSecondaries( 213 std::vector<G4DynamicParticle* 243 std::vector<G4DynamicParticle*>* fvect, 214 const G4MaterialCutsCouple* couple, 244 const G4MaterialCutsCouple* couple, 215 const G4DynamicParticle* dp, 245 const G4DynamicParticle* dp, 216 G4double cutEnergy, 246 G4double cutEnergy, 217 G4double) 247 G4double) 218 { 248 { 219 G4double kinEnergy = dp->GetKineticEnergy(); 249 G4double kinEnergy = dp->GetKineticEnergy(); 220 SetupParticle(dp->GetDefinition()); << 250 if(kinEnergy <= DBL_MIN) return; 221 DefineMaterial(couple); 251 DefineMaterial(couple); 222 /* << 252 SetupParticle(dp->GetDefinition()); 223 G4cout << "G4eCoulombScatteringModel::Sample << 253 G4double ekin = std::max(lowEnergyLimit, kinEnergy); 224 << kinEnergy << " " << particle->GetPart << 254 SetupKinematic(ekin, cutEnergy); 225 << " cut= " << cutEnergy<< G4endl; << 255 //G4cout << "G4eCoulombScatteringModel::SampleSecondaries e(MeV)= " 226 */ << 256 // << kinEnergy << " " << particle->GetParticleName() << G4endl; >> 257 227 // Choose nucleus 258 // Choose nucleus 228 G4double cut = (0.0 < fixedCut) ? fixedCut : << 259 currentElement = SelectRandomAtom(couple,particle,ekin,cutEnergy,ekin); 229 260 230 wokvi->SetupKinematic(kinEnergy, currentMate << 261 SetupTarget(currentElement->GetZ(),ekin); >> 262 >> 263 G4double cost = SampleCosineTheta(); >> 264 G4double z1 = 1.0 - cost; >> 265 if(z1 < 0.0) return; >> 266 >> 267 G4double sint = sqrt(z1*(1.0 + cost)); >> 268 >> 269 //G4cout<<"## Sampled sint= " << sint << " Z= " << targetZ >> 270 // << " screenZ= " << screenZ << " cn= " << formfactA << G4endl; >> 271 >> 272 G4double phi = twopi * G4UniformRand(); 231 273 232 const G4Element* currentElement = SelectTarg << 233 dp->Get << 234 G4int iz = currentElement->GetZasInt(); << 235 << 236 G4double costmin = wokvi->SetupTarget(iz, cu << 237 G4double costmax = (1 == iz && particle == t << 238 ? 0.0 : cosThetaMax; << 239 if(costmin <= costmax) { return; } << 240 << 241 G4double cross = wokvi->ComputeNuclearCrossS << 242 G4double ecross = wokvi->ComputeElectronCros << 243 G4double ratio = ecross/(cross + ecross); << 244 << 245 G4int ia = SelectIsotopeNumber(currentElemen << 246 G4double targetMass = G4NucleiProperties::Ge << 247 wokvi->SetTargetMass(targetMass); << 248 << 249 G4ThreeVector newDirection = << 250 wokvi->SampleSingleScattering(costmin, cos << 251 G4double cost = newDirection.z(); << 252 /* << 253 G4cout << "SampleSec: e(MeV)= " << kinEn << 254 << " 1-costmin= " << 1-costmin << 255 << " 1-costmax= " << 1-costmax << 256 << " 1-cost= " << 1-cost << 257 << " ratio= " << ratio << 258 << G4endl; << 259 */ << 260 G4ThreeVector direction = dp->GetMomentumDir 274 G4ThreeVector direction = dp->GetMomentumDirection(); >> 275 G4ThreeVector newDirection(cos(phi)*sint,sin(phi)*sint,cost); 261 newDirection.rotateUz(direction); 276 newDirection.rotateUz(direction); 262 277 263 fParticleChange->ProposeMomentumDirection(ne 278 fParticleChange->ProposeMomentumDirection(newDirection); 264 279 265 // recoil sampling assuming a small recoil 280 // recoil sampling assuming a small recoil 266 // and first order correction to primary 4-m 281 // and first order correction to primary 4-momentum 267 G4double mom2 = wokvi->GetMomentumSquare(); << 282 if(lowEnergyLimit < kinEnergy) { 268 G4double trec = mom2*(1.0 - cost) << 283 G4int ia = SelectIsotopeNumber(currentElement); 269 /(targetMass + (mass + kinEnergy)*(1.0 - c << 284 G4double Trec = z1*mom2/(amu_c2*G4double(ia)); 270 << 285 G4double th = 271 // the check likely not needed << 286 std::min(recoilThreshold, 272 trec = std::min(trec, kinEnergy); << 287 targetZ*currentElement->GetIonisation()->GetMeanExcitationEnergy()); 273 G4double finalT = kinEnergy - trec; << 288 274 G4double edep = 0.0; << 289 if(Trec > th) { 275 /* << 290 G4int iz = G4int(targetZ); 276 G4cout<<"G4eCoulombScatteringModel: finalT << 291 G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz); 277 <<trec << " Z= " << iz << " A= " << ia << 292 Trec = z1*mom2/ion->GetPDGMass(); 278 << " tcut(keV)= " << (*pCuts)[currentMater << 293 if(Trec < kinEnergy) { 279 */ << 294 G4ThreeVector dir = (direction - newDirection).unit(); 280 G4double tcut = recoilThreshold; << 295 G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, Trec); 281 if(pCuts) { tcut= std::max(tcut,(*pCuts)[cur << 296 fvect->push_back(newdp); 282 << 297 fParticleChange->SetProposedKineticEnergy(kinEnergy - Trec); 283 if(trec > tcut) { << 298 } 284 G4ParticleDefinition* ion = theIonTable->G << 299 } 285 G4ThreeVector dir = (direction*sqrt(mom2) << 300 } 286 newDirection*sqrt(finalT*(2*mass + fina << 301 287 auto newdp = new G4DynamicParticle(ion, di << 302 return; 288 fvect->push_back(newdp); << 303 } 289 } else { << 304 290 edep = trec; << 305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 291 fParticleChange->ProposeNonIonizingEnergyD << 306 >> 307 G4double G4eCoulombScatteringModel::SampleCosineTheta() >> 308 { >> 309 G4double costm = cosTetMaxNuc2; >> 310 G4double formf = formfactA; >> 311 G4double prob = 0.0; >> 312 G4double xs = CrossSectionPerAtom(); >> 313 if(xs > 0.0) prob = elecXSection/xs; >> 314 >> 315 // scattering off e or A? >> 316 if(G4UniformRand() < prob) { >> 317 costm = cosTetMaxElec2; >> 318 formf = 0.0; 292 } 319 } 293 320 294 // finelize primary energy and energy bala << 321 /* 295 // this threshold may be applied only beca << 322 G4cout << "SampleCost: e(MeV)= " << tkin 296 // e+e- msc model is applied << 323 << " ctmin= " << cosThetaMin 297 if(finalT < 0.0) { << 324 << " ctmaxN= " << cosTetMaxNuc 298 edep += finalT; << 325 << " ctmax= " << costm 299 finalT = 0.0; << 326 << " Z= " << targetZ << " A= " << targetA 300 } << 327 << G4endl; 301 edep = std::max(edep, 0.0); << 328 */ 302 fParticleChange->SetProposedKineticEnergy(fi << 329 if(costm >= cosTetMinNuc) return 2.0; 303 fParticleChange->ProposeLocalEnergyDeposit(e << 330 >> 331 G4double x1 = 1. - cosTetMinNuc + screenZ; >> 332 G4double x2 = 1. - costm + screenZ; >> 333 G4double x3 = cosTetMinNuc - costm; >> 334 G4double grej, z1; >> 335 do { >> 336 z1 = x1*x2/(x1 + G4UniformRand()*x3) - screenZ; >> 337 grej = 1.0/(1.0 + formf*z1); >> 338 } while ( G4UniformRand() > grej*grej ); >> 339 >> 340 //G4cout << "z= " << z1 << " cross= " << nucXSection/barn >> 341 // << " crossE= " << elecXSection/barn << G4endl; >> 342 >> 343 return 1.0 - z1; 304 } 344 } 305 345 306 //....oooOO0OOooo........oooOO0OOooo........oo 346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 347 >> 348 307 349