| 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 file 28 // GEANT4 Class file
29 // 29 //
30 // 30 //
31 // File name: G4hCoulombScatteringModel 31 // File name: G4hCoulombScatteringModel
32 // 32 //
33 // Author: Vladimir Ivanchenko 33 // Author: Vladimir Ivanchenko
34 // 34 //
35 // Creation date: 08.06.2012 from G4eCoulombSc 35 // Creation date: 08.06.2012 from G4eCoulombScatteringModel
36 // 36 //
37 // Modifications: 37 // Modifications:
38 // 38 //
39 // 39 //
40 // Class Description: 40 // Class Description:
41 // 41 //
42 // ------------------------------------------- 42 // -------------------------------------------------------------------
43 // 43 //
44 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 46
47 #include "G4hCoulombScatteringModel.hh" 47 #include "G4hCoulombScatteringModel.hh"
48 #include "G4PhysicalConstants.hh" 48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh"
50 #include "Randomize.hh" 50 #include "Randomize.hh"
51 #include "G4DataVector.hh" 51 #include "G4DataVector.hh"
52 #include "G4ElementTable.hh" 52 #include "G4ElementTable.hh"
53 #include "G4ParticleChangeForGamma.hh" 53 #include "G4ParticleChangeForGamma.hh"
54 #include "G4Proton.hh" 54 #include "G4Proton.hh"
55 #include "G4ParticleTable.hh" 55 #include "G4ParticleTable.hh"
56 #include "G4IonTable.hh" 56 #include "G4IonTable.hh"
57 #include "G4ProductionCutsTable.hh" 57 #include "G4ProductionCutsTable.hh"
58 #include "G4NucleiProperties.hh" 58 #include "G4NucleiProperties.hh"
59 #include "G4Pow.hh" 59 #include "G4Pow.hh"
60 #include "G4NistManager.hh" 60 #include "G4NistManager.hh"
61 61
62 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 63
64 G4hCoulombScatteringModel::G4hCoulombScatterin 64 G4hCoulombScatteringModel::G4hCoulombScatteringModel(G4bool combined)
65 : G4VEmModel("hCoulombScattering"), 65 : G4VEmModel("hCoulombScattering"),
66 cosThetaMin(1.0), 66 cosThetaMin(1.0),
67 cosThetaMax(-1.0), 67 cosThetaMax(-1.0),
68 isCombined(combined) 68 isCombined(combined)
69 { 69 {
70 fParticleChange = nullptr; 70 fParticleChange = nullptr;
71 fNistManager = G4NistManager::Instance(); 71 fNistManager = G4NistManager::Instance();
72 theIonTable = G4ParticleTable::GetParticleT 72 theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
73 theProton = G4Proton::Proton(); 73 theProton = G4Proton::Proton();
74 currentMaterial = nullptr; 74 currentMaterial = nullptr;
75 fixedCut = -1.0; 75 fixedCut = -1.0;
76 76
77 pCuts = nullptr; 77 pCuts = nullptr;
78 78
79 recoilThreshold = 0.0; // by default does no 79 recoilThreshold = 0.0; // by default does not work
80 80
81 particle = nullptr; 81 particle = nullptr;
82 currentCouple = nullptr; 82 currentCouple = nullptr;
83 wokvi = new G4WentzelVIRelXSection(); 83 wokvi = new G4WentzelVIRelXSection();
84 84
85 currentMaterialIndex = 0; 85 currentMaterialIndex = 0;
86 mass = CLHEP::proton_mass_c2; 86 mass = CLHEP::proton_mass_c2;
87 elecRatio = 0.0; 87 elecRatio = 0.0;
88 } 88 }
89 89
90 //....oooOO0OOooo........oooOO0OOooo........oo 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 91
92 G4hCoulombScatteringModel::~G4hCoulombScatteri 92 G4hCoulombScatteringModel::~G4hCoulombScatteringModel()
93 { 93 {
94 delete wokvi; 94 delete wokvi;
95 } 95 }
96 96
97 //....oooOO0OOooo........oooOO0OOooo........oo 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 98
99 void G4hCoulombScatteringModel::Initialise(con 99 void G4hCoulombScatteringModel::Initialise(const G4ParticleDefinition* part,
100 const G4DataVector& cuts) 100 const G4DataVector& cuts)
101 { 101 {
102 SetupParticle(part); 102 SetupParticle(part);
103 currentCouple = nullptr; 103 currentCouple = nullptr;
104 104
105 // defined theta limit between single and mu 105 // defined theta limit between single and multiple scattering
106 isCombined = true; 106 isCombined = true;
107 G4double tet = PolarAngleLimit(); 107 G4double tet = PolarAngleLimit();
108 108
109 if(tet <= 0.0) { 109 if(tet <= 0.0) {
110 cosThetaMin = 1.0; 110 cosThetaMin = 1.0;
111 isCombined = false; 111 isCombined = false;
112 } else if(tet >= CLHEP::pi) { 112 } else if(tet >= CLHEP::pi) {
113 cosThetaMin = -1.0; 113 cosThetaMin = -1.0;
114 } else { 114 } else {
115 cosThetaMin = std::cos(tet); 115 cosThetaMin = std::cos(tet);
116 } 116 }
117 117
118 wokvi->Initialise(part, cosThetaMin); 118 wokvi->Initialise(part, cosThetaMin);
119 /* 119 /*
120 G4cout << "G4hCoulombScatteringModel: " << p 120 G4cout << "G4hCoulombScatteringModel: " << particle->GetParticleName()
121 << " 1-cos(ThetaLimit)= " << 1 - cos 121 << " 1-cos(ThetaLimit)= " << 1 - cosThetaMin
122 << " cos(thetaMax)= " << cosThetaMax 122 << " cos(thetaMax)= " << cosThetaMax
123 << G4endl; 123 << G4endl;
124 */ 124 */
125 pCuts = &cuts; 125 pCuts = &cuts;
126 //G4ProductionCutsTable::GetProductionCutsTa 126 //G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
127 /* 127 /*
128 G4cout << "!!! G4hCoulombScatteringModel::In 128 G4cout << "!!! G4hCoulombScatteringModel::Initialise for "
129 << part->GetParticleName() << " cos(TetM 129 << part->GetParticleName() << " cos(TetMin)= " << cosThetaMin
130 << " cos(TetMax)= " << cosThetaMax <<G4e 130 << " cos(TetMax)= " << cosThetaMax <<G4endl;
131 G4cout << "cut= " << (*pCuts)[0] << " cut1= 131 G4cout << "cut= " << (*pCuts)[0] << " cut1= " << (*pCuts)[1] << G4endl;
132 */ 132 */
133 if(!fParticleChange) { 133 if(!fParticleChange) {
134 fParticleChange = GetParticleChangeForGamm 134 fParticleChange = GetParticleChangeForGamma();
135 } 135 }
136 if(IsMaster() && mass < CLHEP::GeV && part-> 136 if(IsMaster() && mass < CLHEP::GeV && part->GetParticleName() != "GenericIon") {
137 InitialiseElementSelectors(part, cuts); 137 InitialiseElementSelectors(part, cuts);
138 } 138 }
139 } 139 }
140 140
141 //....oooOO0OOooo........oooOO0OOooo........oo 141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
142 142
143 void G4hCoulombScatteringModel::InitialiseLoca 143 void G4hCoulombScatteringModel::InitialiseLocal(const G4ParticleDefinition*,
144 G4VEmModel* masterModel) 144 G4VEmModel* masterModel)
145 { 145 {
146 SetElementSelectors(masterModel->GetElementS 146 SetElementSelectors(masterModel->GetElementSelectors());
147 } 147 }
148 148
149 //....oooOO0OOooo........oooOO0OOooo........oo 149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150 150
151 G4double 151 G4double
152 G4hCoulombScatteringModel::MinPrimaryEnergy(co 152 G4hCoulombScatteringModel::MinPrimaryEnergy(const G4Material* material,
153 const G4ParticleDefinition* part 153 const G4ParticleDefinition* part,
154 G4double) 154 G4double)
155 { 155 {
156 SetupParticle(part); 156 SetupParticle(part);
157 157
158 // define cut using cuts for proton 158 // define cut using cuts for proton
159 G4double cut = 159 G4double cut =
160 std::max(recoilThreshold, (*pCuts)[Current 160 std::max(recoilThreshold, (*pCuts)[CurrentCouple()->GetIndex()]);
161 161
162 // find out lightest element 162 // find out lightest element
163 const G4ElementVector* theElementVector = ma 163 const G4ElementVector* theElementVector = material->GetElementVector();
164 std::size_t nelm = material->GetNumberOfElem << 164 G4int nelm = material->GetNumberOfElements();
165 165
166 // select lightest element 166 // select lightest element
167 G4int Z = 300; 167 G4int Z = 300;
168 for (std::size_t j=0; j<nelm; ++j) { << 168 for (G4int j=0; j<nelm; ++j) {
169 Z = std::min(Z,(*theElementVector)[j]->Get 169 Z = std::min(Z,(*theElementVector)[j]->GetZasInt());
170 } 170 }
171 G4int A = G4lrint(fNistManager->GetAtomicMas 171 G4int A = G4lrint(fNistManager->GetAtomicMassAmu(Z));
172 G4double targetMass = G4NucleiProperties::Ge 172 G4double targetMass = G4NucleiProperties::GetNuclearMass(A, Z);
173 G4double t = std::max(cut, 0.5*(cut + std::s 173 G4double t = std::max(cut, 0.5*(cut + std::sqrt(2*cut*targetMass)));
174 174
175 return t; 175 return t;
176 } 176 }
177 177
178 //....oooOO0OOooo........oooOO0OOooo........oo 178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
179 179
180 G4double G4hCoulombScatteringModel::ComputeCro 180 G4double G4hCoulombScatteringModel::ComputeCrossSectionPerAtom(
181 const G4ParticleDefinition* p, 181 const G4ParticleDefinition* p,
182 G4double kinEnergy, 182 G4double kinEnergy,
183 G4double Z, G4double, 183 G4double Z, G4double,
184 G4double cutEnergy, G4double) 184 G4double cutEnergy, G4double)
185 { 185 {
186 //G4cout << "### G4hCoulombScatteringModel:: 186 //G4cout << "### G4hCoulombScatteringModel::ComputeCrossSectionPerAtom for "
187 //<< p->GetParticleName()<<" Z= "<<Z<<" e(Me 187 //<< p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV << G4endl;
188 G4double cross = 0.0; 188 G4double cross = 0.0;
189 elecRatio = 0.0; 189 elecRatio = 0.0;
190 if(p != particle) { SetupParticle(p); } 190 if(p != particle) { SetupParticle(p); }
191 191
192 // cross section is set to zero to avoid pro 192 // cross section is set to zero to avoid problems in sample secondary
193 if(kinEnergy <= 0.0) { return cross; } 193 if(kinEnergy <= 0.0) { return cross; }
194 DefineMaterial(CurrentCouple()); 194 DefineMaterial(CurrentCouple());
195 195
196 G4int iz = G4lrint(Z); 196 G4int iz = G4lrint(Z);
197 G4double tmass = (1 == iz) ? proton_mass_c2 197 G4double tmass = (1 == iz) ? proton_mass_c2 :
198 fNistManager->GetAtomicMassAmu(iz)*amu_c2; 198 fNistManager->GetAtomicMassAmu(iz)*amu_c2;
199 wokvi->SetTargetMass(tmass); 199 wokvi->SetTargetMass(tmass);
200 200
201 G4double costmin = 201 G4double costmin =
202 wokvi->SetupKinematic(kinEnergy, currentMa 202 wokvi->SetupKinematic(kinEnergy, currentMaterial);
203 203
204 if(cosThetaMax < costmin) { 204 if(cosThetaMax < costmin) {
205 G4double cut = (0.0 < fixedCut) ? fixedCut 205 G4double cut = (0.0 < fixedCut) ? fixedCut : cutEnergy;
206 costmin = wokvi->SetupTarget(iz, cut); 206 costmin = wokvi->SetupTarget(iz, cut);
207 G4double costmax = 207 G4double costmax =
208 (1 == iz && particle == theProton && cos 208 (1 == iz && particle == theProton && cosThetaMax < 0.0)
209 ? 0.0 : cosThetaMax; 209 ? 0.0 : cosThetaMax;
210 if(costmin > costmax) { 210 if(costmin > costmax) {
211 cross = wokvi->ComputeNuclearCrossSectio 211 cross = wokvi->ComputeNuclearCrossSection(costmin, costmax)
212 + wokvi->ComputeElectronCrossSection(costmin 212 + wokvi->ComputeElectronCrossSection(costmin, costmax);
213 } 213 }
214 /* 214 /*
215 if(p->GetParticleName() == "mu+") 215 if(p->GetParticleName() == "mu+")
216 G4cout << "e(MeV)= " << kinEnergy/MeV << " c 216 G4cout << "e(MeV)= " << kinEnergy/MeV << " cross(b)= " << cross/barn
217 << " 1-costmin= " << 1-costmin 217 << " 1-costmin= " << 1-costmin
218 << " 1-costmax= " << 1-costmax 218 << " 1-costmax= " << 1-costmax
219 << " 1-cosThetaMax= " << 1-cosThetaMax 219 << " 1-cosThetaMax= " << 1-cosThetaMax
220 << G4endl; 220 << G4endl;
221 */ 221 */
222 } 222 }
223 return cross; 223 return cross;
224 } 224 }
225 225
226 //....oooOO0OOooo........oooOO0OOooo........oo 226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227 227
228 void G4hCoulombScatteringModel::SampleSecondar 228 void G4hCoulombScatteringModel::SampleSecondaries(
229 std::vector<G4DynamicParticle* 229 std::vector<G4DynamicParticle*>* fvect,
230 const G4MaterialCutsCouple* couple, 230 const G4MaterialCutsCouple* couple,
231 const G4DynamicParticle* dp, 231 const G4DynamicParticle* dp,
232 G4double cutEnergy, 232 G4double cutEnergy,
233 G4double) 233 G4double)
234 { 234 {
235 G4double kinEnergy = dp->GetKineticEnergy(); 235 G4double kinEnergy = dp->GetKineticEnergy();
236 SetupParticle(dp->GetDefinition()); 236 SetupParticle(dp->GetDefinition());
237 DefineMaterial(couple); 237 DefineMaterial(couple);
238 238
239 // Choose nucleus 239 // Choose nucleus
240 G4double cut = (0.0 < fixedCut) ? fixedCut : 240 G4double cut = (0.0 < fixedCut) ? fixedCut : cutEnergy;
241 241
242 const G4Element* elm = SelectRandomAtom(coup 242 const G4Element* elm = SelectRandomAtom(couple,particle,
243 kinEnergy,cut,kinEnergy); 243 kinEnergy,cut,kinEnergy);
244 244
245 G4int iz = elm->GetZasInt(); 245 G4int iz = elm->GetZasInt();
246 G4int ia = SelectIsotopeNumber(elm); 246 G4int ia = SelectIsotopeNumber(elm);
247 G4double mass2 = G4NucleiProperties::GetNucl 247 G4double mass2 = G4NucleiProperties::GetNuclearMass(ia, iz);
248 248
249 wokvi->SetTargetMass(mass2); 249 wokvi->SetTargetMass(mass2);
250 wokvi->SetupKinematic(kinEnergy, currentMate 250 wokvi->SetupKinematic(kinEnergy, currentMaterial);
251 G4double costmin = wokvi->SetupTarget(iz, cu 251 G4double costmin = wokvi->SetupTarget(iz, cut);
252 G4double costmax = (1 == iz && particle == t 252 G4double costmax = (1 == iz && particle == theProton && cosThetaMax < 0.0)
253 ? 0.0 : cosThetaMax; 253 ? 0.0 : cosThetaMax;
254 if(costmin <= costmax) { return; } 254 if(costmin <= costmax) { return; }
255 255
256 G4double cross = wokvi->ComputeNuclearCrossS 256 G4double cross = wokvi->ComputeNuclearCrossSection(costmin, costmax);
257 G4double ecross = wokvi->ComputeElectronCros 257 G4double ecross = wokvi->ComputeElectronCrossSection(costmin, costmax);
258 G4double ratio = ecross/(cross + ecross); 258 G4double ratio = ecross/(cross + ecross);
259 259
260 G4ThreeVector newDirection = 260 G4ThreeVector newDirection =
261 wokvi->SampleSingleScattering(costmin, cos 261 wokvi->SampleSingleScattering(costmin, costmax, ratio);
262 262
263 // kinematics in the Lab system 263 // kinematics in the Lab system
264 G4double ptot = std::sqrt(kinEnergy*(kinEner 264 G4double ptot = std::sqrt(kinEnergy*(kinEnergy + 2.0*mass));
265 G4double e1 = mass + kinEnergy; 265 G4double e1 = mass + kinEnergy;
266 266
267 // Lab. system kinematics along projectile d 267 // Lab. system kinematics along projectile direction
268 G4LorentzVector v0 = G4LorentzVector(0, 0, p 268 G4LorentzVector v0 = G4LorentzVector(0, 0, ptot, e1+mass2);
269 G4LorentzVector v1 = G4LorentzVector(0, 0, p 269 G4LorentzVector v1 = G4LorentzVector(0, 0, ptot, e1);
270 G4ThreeVector bst = v0.boostVector(); 270 G4ThreeVector bst = v0.boostVector();
271 v1.boost(-bst); 271 v1.boost(-bst);
272 // CM projectile 272 // CM projectile
273 G4double momCM = v1.pz(); 273 G4double momCM = v1.pz();
274 274
275 // Momentum after scattering of incident par 275 // Momentum after scattering of incident particle
276 v1.setX(momCM*newDirection.x()); 276 v1.setX(momCM*newDirection.x());
277 v1.setY(momCM*newDirection.y()); 277 v1.setY(momCM*newDirection.y());
278 v1.setZ(momCM*newDirection.z()); 278 v1.setZ(momCM*newDirection.z());
279 279
280 // CM--->Lab 280 // CM--->Lab
281 v1.boost(bst); 281 v1.boost(bst);
282 282
283 G4ThreeVector dir = dp->GetMomentumDirection 283 G4ThreeVector dir = dp->GetMomentumDirection();
284 newDirection = v1.vect().unit(); 284 newDirection = v1.vect().unit();
285 newDirection.rotateUz(dir); 285 newDirection.rotateUz(dir);
286 286
287 fParticleChange->ProposeMomentumDirection(ne 287 fParticleChange->ProposeMomentumDirection(newDirection);
288 288
289 // recoil 289 // recoil
290 v0 -= v1; 290 v0 -= v1;
291 G4double trec = std::max(v0.e() - mass2, 0.0 291 G4double trec = std::max(v0.e() - mass2, 0.0);
292 G4double edep = 0.0; 292 G4double edep = 0.0;
293 293
294 G4double tcut = recoilThreshold; 294 G4double tcut = recoilThreshold;
295 if(pCuts) { tcut= std::max(tcut,(*pCuts)[cur 295 if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); }
296 296
297 if(trec > tcut) { 297 if(trec > tcut) {
298 G4ParticleDefinition* ion = theIonTable->G 298 G4ParticleDefinition* ion = theIonTable->GetIon(iz, ia, 0);
299 newDirection = v0.vect().unit(); 299 newDirection = v0.vect().unit();
300 newDirection.rotateUz(dir); 300 newDirection.rotateUz(dir);
301 auto newdp = new G4DynamicParticle(ion, ne << 301 G4DynamicParticle* newdp = new G4DynamicParticle(ion, newDirection, trec);
302 fvect->push_back(newdp); 302 fvect->push_back(newdp);
303 } else if(trec > 0.0) { 303 } else if(trec > 0.0) {
304 edep = trec; 304 edep = trec;
305 fParticleChange->ProposeNonIonizingEnergyD 305 fParticleChange->ProposeNonIonizingEnergyDeposit(edep);
306 } 306 }
307 307
308 // finelize primary energy and energy balanc 308 // finelize primary energy and energy balance
309 G4double finalT = v1.e() - mass; 309 G4double finalT = v1.e() - mass;
310 if(finalT < 0.0) { 310 if(finalT < 0.0) {
311 edep += finalT; 311 edep += finalT;
312 finalT = 0.0; 312 finalT = 0.0;
313 } 313 }
314 edep = std::max(edep, 0.0); 314 edep = std::max(edep, 0.0);
315 fParticleChange->SetProposedKineticEnergy(fi 315 fParticleChange->SetProposedKineticEnergy(finalT);
316 fParticleChange->ProposeLocalEnergyDeposit(e 316 fParticleChange->ProposeLocalEnergyDeposit(edep);
317 } 317 }
318 318
319 //....oooOO0OOooo........oooOO0OOooo........oo 319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
320 320