| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
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 G4eCoulombScatteringModel
36 // & G4CoulombScatteringModel
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-Ionizing Energy-Loss
44 // for Coulomb ScatteredParticles from Low Energy up to Relativistic
45 // Regime in Space Radiation Environment"
46 // Accepted for publication in the Proceedings of the ICATPP Conference
47 // on Cosmic Rays for Particle and Astroparticle Physics, Villa Olmo, 7-8
48 // October, 2010, to be published by World Scientific (Singapore).
49 //
50 // Available for downloading at:
51 // http://arxiv.org/abs/1011.4822
52 //
53 // -------------------------------------------------------------------
54 //
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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........oooOO0OOooo........oooOO0OOooo....
72
73 using namespace std;
74
75 G4IonCoulombScatteringModel::G4IonCoulombScatteringModel(const G4String& nam)
76 : G4VEmModel(nam),
77 cosThetaMin(1.0)
78 {
79 fNistManager = G4NistManager::Instance();
80 theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
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........oooOO0OOooo........oooOO0OOooo....
100
101 G4IonCoulombScatteringModel::~G4IonCoulombScatteringModel()
102 {
103 delete ioncross;
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107
108 void G4IonCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
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::GetProductionCutsTable()->GetEnergyCutsVector(3);
118 if(!fParticleChange) {
119 fParticleChange = GetParticleChangeForGamma();
120 }
121 }
122
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124
125 G4double G4IonCoulombScatteringModel::ComputeCrossSectionPerAtom(
126 const G4ParticleDefinition* p,
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 particle
140 G4double tmass = proton_mass_c2;
141 if(1 < iz) {
142 tmass = fNistManager->GetAtomicMassAmu(iz)*amu_c2;
143 }
144 ioncross->SetupKinematic(kinEnergy, tmass);
145 ioncross->SetupTarget(Z, kinEnergy, heavycorr);
146 cross = ioncross->NuclearCrossSection();
147 return cross;
148 }
149
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151
152 void G4IonCoulombScatteringModel::SampleSecondaries(
153 std::vector<G4DynamicParticle*>* fvect,
154 const G4MaterialCutsCouple* couple,
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, particle, kinEnergy,
164 dp->GetLogKineticEnergy());
165
166 G4int iz = currentElement->GetZasInt();
167 G4int ia = SelectIsotopeNumber(currentElement);
168 G4double mass2 = G4NucleiProperties::GetNuclearMass(ia, iz);
169
170 ioncross->SetupKinematic(kinEnergy, mass2);
171 ioncross->SetupTarget(currentElement->GetZ(), kinEnergy, heavycorr);
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->GetParticleName()
179 << " mass(GeV)= " << mass/GeV
180 << " Ekin(MeV)= " << kinEnergy << " cost= " << 1. - z1 << G4endl;
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 + 2.0*mass));
190 G4double e1 = mass + kinEnergy;
191
192 // Lab. system kinematics along projectile direction
193 G4LorentzVector v0 = G4LorentzVector(0, 0, ptot, e1+mass2);
194 G4LorentzVector v1 = G4LorentzVector(0, 0, ptot, e1);
195 G4ThreeVector bst = v0.boostVector();
196 v1.boost(-bst);
197 // CM projectile
198 G4double momCM = v1.pz();
199
200 // Momentum after scattering of incident particle
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(newDirection);
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)[currentMaterialIndex]);
223 //G4cout<<" tcut eV "<<tcut/eV<<endl;
224 }
225
226 // Recoil
227 if(trec > tcut) {
228 G4ParticleDefinition* ion = theIonTable->GetIon(iz, ia, 0);
229 newDirection = v0.vect().unit();
230 newDirection.rotateUz(dir);
231 auto newdp = new G4DynamicParticle(ion, newDirection, trec);
232 fvect->push_back(newdp);
233 } else if(trec > 0.0) {
234 edep = trec;
235 fParticleChange->ProposeNonIonizingEnergyDeposit(edep);
236 }
237
238 // finelize primary energy and energy balance
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 << " Edep(MeV)= " << edep
246 // << " Trec(MeV)= " << trec << G4endl;
247 fParticleChange->SetProposedKineticEnergy(finalT);
248 fParticleChange->ProposeLocalEnergyDeposit(edep);
249 }
250
251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
252
253