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