| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4IonCoulombScatteringModel.cc
27 // -------------------------------------------
28 //
29 // GEANT4 Class header file
30 //
31 // File name: G4IonCoulombScatteringModel
32 //
33 // Author: Cristina Consolandi
34 //
35 // Creation date: 05.10.2010 from G4eCoulombSc
36 // & G4CoulombSc
37 //
38 // Class Description:
39 // Single Scattering Model for
40 // for protons, alpha and heavy Ions
41 //
42 // Reference:
43 // M.J. Boschini et al. "Nuclear and Non-
44 // for Coulomb ScatteredParticles from Lo
45 // Regime in Space Radiation Environment"
46 // Accepted for publication in the Procee
47 // on Cosmic Rays for Particle and Astrop
48 // October, 2010, to be published by Wor
49 //
50 // Available for downloading at:
51 // http://arxiv.org/abs/1011.4822
52 //
53 // -------------------------------------------
54 //
55 //....oooOO0OOooo........oooOO0OOooo........oo
56
57
58 #include "G4IonCoulombScatteringModel.hh"
59 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh"
61 #include "Randomize.hh"
62 #include "G4ParticleChangeForGamma.hh"
63 #include "G4Proton.hh"
64 #include "G4ProductionCutsTable.hh"
65 #include "G4NucleiProperties.hh"
66 #include "G4ParticleTable.hh"
67 #include "G4IonTable.hh"
68
69 #include "G4UnitsTable.hh"
70
71 //....oooOO0OOooo........oooOO0OOooo........oo
72
73 using namespace std;
74
75 G4IonCoulombScatteringModel::G4IonCoulombScatt
76 : G4VEmModel(nam),
77 cosThetaMin(1.0)
78 {
79 fNistManager = G4NistManager::Instance();
80 theIonTable = G4ParticleTable::GetParticleT
81 theProton = G4Proton::Proton();
82
83 pCuts = nullptr;
84 currentMaterial = nullptr;
85 currentElement = nullptr;
86 currentCouple = nullptr;
87 fParticleChange = nullptr;
88
89 recoilThreshold = 0.*eV;
90 heavycorr =0;
91 particle = nullptr;
92 mass=0;
93 currentMaterialIndex = -1;
94
95 ioncross = new G4IonCoulombCrossSection();
96 }
97
98
99 //....oooOO0OOooo........oooOO0OOooo........oo
100
101 G4IonCoulombScatteringModel::~G4IonCoulombScat
102 {
103 delete ioncross;
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oo
107
108 void G4IonCoulombScatteringModel::Initialise(c
109 const G4DataVector& cuts)
110 {
111 SetupParticle(p);
112 currentCouple = nullptr;
113 currentMaterialIndex = -1;
114 ioncross->Initialise(p,cosThetaMin);
115
116 pCuts = &cuts;
117 // G4ProductionCutsTable::GetProductionCuts
118 if(!fParticleChange) {
119 fParticleChange = GetParticleChangeForGamm
120 }
121 }
122
123 //....oooOO0OOooo........oooOO0OOooo........oo
124
125 G4double G4IonCoulombScatteringModel::ComputeC
126 const G4Partic
127 G4double kinEnergy,
128 G4double Z,
129 G4double, G4double, G4double)
130 {
131 SetupParticle(p);
132
133 G4double cross = 0.0;
134
135 DefineMaterial(CurrentCouple());
136
137 G4int iz = G4lrint(Z);
138
139 //from lab to pCM & mu_rel of effective part
140 G4double tmass = proton_mass_c2;
141 if(1 < iz) {
142 tmass = fNistManager->GetAtomicMassAmu(iz)
143 }
144 ioncross->SetupKinematic(kinEnergy, tmass);
145 ioncross->SetupTarget(Z, kinEnergy, heavycor
146 cross = ioncross->NuclearCrossSection();
147 return cross;
148 }
149
150 //....oooOO0OOooo........oooOO0OOooo........oo
151
152 void G4IonCoulombScatteringModel::SampleSecond
153 std::vector<G4DynamicParticle*>*
154 const G4MaterialCutsCouple* coupl
155 const G4DynamicParticle* dp,
156 G4double, G4double)
157 {
158 G4double kinEnergy = dp->GetKineticEnergy();
159 DefineMaterial(couple);
160 SetupParticle(dp->GetDefinition());
161
162 // Choose nucleus
163 currentElement = SelectTargetAtom(couple, pa
164 dp->GetLog
165
166 G4int iz = currentElement->GetZasInt();
167 G4int ia = SelectIsotopeNumber(currentElemen
168 G4double mass2 = G4NucleiProperties::GetNucl
169
170 ioncross->SetupKinematic(kinEnergy, mass2);
171 ioncross->SetupTarget(currentElement->GetZ()
172
173 //scattering angle, z1 == (1-cost)
174 G4double z1 = ioncross->SampleCosineTheta();
175 if(z1 > 2.0) { z1 = 2.0; }
176 else if(z1 < 0.0) { z1 = 0.0; }
177 /*
178 G4cout << "Sample: " << particle->GetParticl
179 << " mass(GeV)= " << mass/GeV
180 << " Ekin(MeV)= " << kinEnergy << " cost= "
181 G4cout << " Z= " << iz << " A= " << ia
182 << " mass(GeV)= " << mass2/GeV << G4endl;
183 */
184 G4double cost = 1.0 - z1;
185 G4double sint = sqrt(z1*(1.0 + cost));
186 G4double phi = twopi * G4UniformRand();
187
188 // kinematics in the Lab system
189 G4double ptot = sqrt(kinEnergy*(kinEnergy +
190 G4double e1 = mass + kinEnergy;
191
192 // Lab. system kinematics along projectile d
193 G4LorentzVector v0 = G4LorentzVector(0, 0, p
194 G4LorentzVector v1 = G4LorentzVector(0, 0, p
195 G4ThreeVector bst = v0.boostVector();
196 v1.boost(-bst);
197 // CM projectile
198 G4double momCM = v1.pz();
199
200 // Momentum after scattering of incident par
201 v1.setX(momCM*sint*cos(phi));
202 v1.setY(momCM*sint*sin(phi));
203 v1.setZ(momCM*cost);
204
205 // CM--->Lab
206 v1.boost(bst);
207
208 // Rotate to global system
209 G4ThreeVector dir = dp->GetMomentumDirection
210 G4ThreeVector newDirection = v1.vect().unit(
211 newDirection.rotateUz(dir);
212
213 fParticleChange->ProposeMomentumDirection(ne
214
215 // recoil v0 energy is kinetic
216 v0 -= v1;
217 G4double trec = std::max(v0.e() - mass2, 0.0
218 G4double edep = 0.0;
219
220 G4double tcut = recoilThreshold;
221 if(pCuts) {
222 tcut= std::max(tcut,(*pCuts)[currentMateri
223 //G4cout<<" tcut eV "<<tcut/eV<<endl;
224 }
225
226 // Recoil
227 if(trec > tcut) {
228 G4ParticleDefinition* ion = theIonTable->G
229 newDirection = v0.vect().unit();
230 newDirection.rotateUz(dir);
231 auto newdp = new G4DynamicParticle(ion, ne
232 fvect->push_back(newdp);
233 } else if(trec > 0.0) {
234 edep = trec;
235 fParticleChange->ProposeNonIonizingEnergyD
236 }
237
238 // finelize primary energy and energy balanc
239 G4double finalT = v1.e() - mass;
240 if(finalT < 0.0) {
241 edep += finalT;
242 finalT = 0.0;
243 }
244 edep = std::max(edep, 0.0);
245 //G4cout << "Efinal(MeV)= " << finalT << " E
246 // << " Trec(MeV)= " << trec << G4endl;
247 fParticleChange->SetProposedKineticEnergy(fi
248 fParticleChange->ProposeLocalEnergyDeposit(e
249 }
250
251 //....oooOO0OOooo........oooOO0OOooo........oo
252
253