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