| 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
>> 148 //cout<< "..........cross "<<G4BestUnit(cross,"Surface") <<endl;
147 return cross; 149 return cross;
148 } 150 }
149 151
150 //....oooOO0OOooo........oooOO0OOooo........oo 152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 153
152 void G4IonCoulombScatteringModel::SampleSecond 154 void G4IonCoulombScatteringModel::SampleSecondaries(
153 std::vector<G4DynamicParticle*>* 155 std::vector<G4DynamicParticle*>* fvect,
154 const G4MaterialCutsCouple* coupl 156 const G4MaterialCutsCouple* couple,
155 const G4DynamicParticle* dp, 157 const G4DynamicParticle* dp,
156 G4double, G4double) 158 G4double, G4double)
157 { 159 {
158 G4double kinEnergy = dp->GetKineticEnergy(); 160 G4double kinEnergy = dp->GetKineticEnergy();
>> 161
159 DefineMaterial(couple); 162 DefineMaterial(couple);
>> 163
160 SetupParticle(dp->GetDefinition()); 164 SetupParticle(dp->GetDefinition());
161 165
162 // Choose nucleus 166 // Choose nucleus
163 currentElement = SelectTargetAtom(couple, pa << 167 currentElement = SelectRandomAtom(couple, particle, kinEnergy);
164 dp->GetLog <<
165 168
166 G4int iz = currentElement->GetZasInt(); 169 G4int iz = currentElement->GetZasInt();
167 G4int ia = SelectIsotopeNumber(currentElemen 170 G4int ia = SelectIsotopeNumber(currentElement);
168 G4double mass2 = G4NucleiProperties::GetNucl 171 G4double mass2 = G4NucleiProperties::GetNuclearMass(ia, iz);
169 172
170 ioncross->SetupKinematic(kinEnergy, mass2); 173 ioncross->SetupKinematic(kinEnergy, mass2);
>> 174
171 ioncross->SetupTarget(currentElement->GetZ() 175 ioncross->SetupTarget(currentElement->GetZ(), kinEnergy, heavycorr);
172 176
173 //scattering angle, z1 == (1-cost) 177 //scattering angle, z1 == (1-cost)
174 G4double z1 = ioncross->SampleCosineTheta(); 178 G4double z1 = ioncross->SampleCosineTheta();
175 if(z1 > 2.0) { z1 = 2.0; } 179 if(z1 > 2.0) { z1 = 2.0; }
176 else if(z1 < 0.0) { z1 = 0.0; } 180 else if(z1 < 0.0) { z1 = 0.0; }
177 /* << 181
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; 182 G4double cost = 1.0 - z1;
185 G4double sint = sqrt(z1*(1.0 + cost)); 183 G4double sint = sqrt(z1*(1.0 + cost));
186 G4double phi = twopi * G4UniformRand(); 184 G4double phi = twopi * G4UniformRand();
187 185
188 // kinematics in the Lab system 186 // kinematics in the Lab system
189 G4double ptot = sqrt(kinEnergy*(kinEnergy + << 187 G4double ptot = dp->GetTotalMomentum();
190 G4double e1 = mass + kinEnergy; << 188 G4double e1 = dp->GetTotalEnergy();
191 189
192 // Lab. system kinematics along projectile d 190 // Lab. system kinematics along projectile direction
193 G4LorentzVector v0 = G4LorentzVector(0, 0, p << 191 G4LorentzVector v0 = G4LorentzVector(0, 0, ptot, e1);
194 G4LorentzVector v1 = G4LorentzVector(0, 0, p << 192 G4double bet = ptot/(e1 + mass2);
195 G4ThreeVector bst = v0.boostVector(); << 193 G4double gam = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
196 v1.boost(-bst); << 194
197 // CM projectile 195 // CM projectile
198 G4double momCM = v1.pz(); << 196 G4double momCM = gam*(ptot - bet*e1);
>> 197 G4double eCM = gam*(e1 - bet*ptot);
199 198
200 // Momentum after scattering of incident par 199 // Momentum after scattering of incident particle
201 v1.setX(momCM*sint*cos(phi)); << 200 G4double pxCM = momCM*sint*cos(phi);
202 v1.setY(momCM*sint*sin(phi)); << 201 G4double pyCM = momCM*sint*sin(phi);
203 v1.setZ(momCM*cost); << 202 G4double pzCM = momCM*cost;
204 203
205 // CM--->Lab 204 // CM--->Lab
206 v1.boost(bst); << 205 G4LorentzVector v1(pxCM , pyCM, gam*(pzCM + bet*eCM), gam*(eCM + bet*pzCM));
207 206
208 // Rotate to global system 207 // Rotate to global system
209 G4ThreeVector dir = dp->GetMomentumDirection 208 G4ThreeVector dir = dp->GetMomentumDirection();
210 G4ThreeVector newDirection = v1.vect().unit( 209 G4ThreeVector newDirection = v1.vect().unit();
211 newDirection.rotateUz(dir); 210 newDirection.rotateUz(dir);
212 211
213 fParticleChange->ProposeMomentumDirection(ne 212 fParticleChange->ProposeMomentumDirection(newDirection);
214 213
215 // recoil v0 energy is kinetic 214 // recoil v0 energy is kinetic
216 v0 -= v1; 215 v0 -= v1;
217 G4double trec = std::max(v0.e() - mass2, 0.0 << 216 G4double trec = v0.e();
218 G4double edep = 0.0; 217 G4double edep = 0.0;
219 218
220 G4double tcut = recoilThreshold; 219 G4double tcut = recoilThreshold;
221 if(pCuts) { 220 if(pCuts) {
222 tcut= std::max(tcut,(*pCuts)[currentMateri 221 tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]);
223 //G4cout<<" tcut eV "<<tcut/eV<<endl; 222 //G4cout<<" tcut eV "<<tcut/eV<<endl;
224 } 223 }
225 224
226 // Recoil 225 // Recoil
227 if(trec > tcut) { 226 if(trec > tcut) {
228 G4ParticleDefinition* ion = theIonTable->G 227 G4ParticleDefinition* ion = theIonTable->GetIon(iz, ia, 0);
229 newDirection = v0.vect().unit(); 228 newDirection = v0.vect().unit();
230 newDirection.rotateUz(dir); 229 newDirection.rotateUz(dir);
231 auto newdp = new G4DynamicParticle(ion, ne << 230 G4DynamicParticle* newdp = new G4DynamicParticle(ion, newDirection, trec);
232 fvect->push_back(newdp); 231 fvect->push_back(newdp);
233 } else if(trec > 0.0) { 232 } else if(trec > 0.0) {
234 edep = trec; 233 edep = trec;
235 fParticleChange->ProposeNonIonizingEnergyD 234 fParticleChange->ProposeNonIonizingEnergyDeposit(edep);
236 } 235 }
237 236
238 // finelize primary energy and energy balanc 237 // finelize primary energy and energy balance
239 G4double finalT = v1.e() - mass; 238 G4double finalT = v1.e() - mass;
240 if(finalT < 0.0) { 239 if(finalT < 0.0) {
241 edep += finalT; 240 edep += finalT;
242 finalT = 0.0; 241 finalT = 0.0;
243 } 242 }
244 edep = std::max(edep, 0.0); 243 edep = std::max(edep, 0.0);
245 //G4cout << "Efinal(MeV)= " << finalT << " E <<
246 // << " Trec(MeV)= " << trec << G4endl; <<
247 fParticleChange->SetProposedKineticEnergy(fi 244 fParticleChange->SetProposedKineticEnergy(finalT);
248 fParticleChange->ProposeLocalEnergyDeposit(e 245 fParticleChange->ProposeLocalEnergyDeposit(edep);
249 } 246 }
250 247
251 //....oooOO0OOooo........oooOO0OOooo........oo 248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
252 249
253 250