source: trunk/source/processes/electromagnetic/standard/src/G4IonCoulombScatteringModel.cc@ 1357

Last change on this file since 1357 was 1350, checked in by garnier, 15 years ago

update to last version 4.9.4
File size: 8.1 KB
RevLine 
[1350]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 "Randomize.hh"
60//#include "G4DataVector.hh"
61#include "G4ParticleChangeForGamma.hh"
62#include "G4Proton.hh"
63#include "G4ProductionCutsTable.hh"
64#include "G4NucleiProperties.hh"
65
66#include "G4UnitsTable.hh"
67
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69
70using namespace std;
71
72G4IonCoulombScatteringModel::G4IonCoulombScatteringModel(const G4String& nam)
73 : G4VEmModel(nam),
74
75 cosThetaMin(1.0),
76 isInitialised(false)
77{
78 fNistManager = G4NistManager::Instance();
79 theParticleTable = G4ParticleTable::GetParticleTable();
80 theProton = G4Proton::Proton();
81
82 pCuts=0;
83 currentMaterial = 0;
84 currentElement = 0;
85 currentCouple = 0;
86
87 lowEnergyLimit = 100*eV;
88 recoilThreshold = 0.*eV;
89 heavycorr =0;
90 particle = 0;
91 mass=0;
92
93 ioncross = new G4IonCoulombCrossSection();
94
95}
96
97
98//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99
100G4IonCoulombScatteringModel::~G4IonCoulombScatteringModel()
101{ delete ioncross;}
102
103//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104
105void G4IonCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
106 const G4DataVector& )
107{
108 SetupParticle(p);
109 currentCouple = 0;
110 cosThetaMin = cos(PolarAngleLimit());
111 ioncross->Initialise(p,cosThetaMin);
112
113 pCuts = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
114
115
116 if(!isInitialised) {
117 isInitialised = true;
118 fParticleChange = GetParticleChangeForGamma();
119 }
120}
121
122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123
124G4double G4IonCoulombScatteringModel::ComputeCrossSectionPerAtom(
125 const G4ParticleDefinition* p,
126 G4double kinEnergy,
127 G4double Z,
128 G4double,
129 G4double cutEnergy,
130 G4double)
131{
132
133 SetupParticle(p);
134
135 G4double xsec =0.0;
136 if(kinEnergy < lowEnergyLimit) return xsec;
137
138 DefineMaterial(CurrentCouple());
139
140 G4int iz = G4int(Z);
141
142 //from lab to pCM & mu_rel of effective particle
143 ioncross->SetupKinematic(kinEnergy, cutEnergy,iz);
144
145
146 ioncross->SetupTarget(Z, kinEnergy, heavycorr);
147
148
149 xsec = ioncross->NuclearCrossSection();
150
151//cout<< "..........xsec "<<G4BestUnit(xsec,"Surface") <<endl;
152 return xsec;
153}
154
155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156
157void G4IonCoulombScatteringModel::SampleSecondaries(
158 std::vector<G4DynamicParticle*>* fvect,
159 const G4MaterialCutsCouple* couple,
160 const G4DynamicParticle* dp,
161 G4double cutEnergy,
162 G4double)
163{
164 G4double kinEnergy = dp->GetKineticEnergy();
165
166 if(kinEnergy < lowEnergyLimit) return;
167
168 DefineMaterial(couple);
169
170 SetupParticle(dp->GetDefinition());
171
172 // Choose nucleus
173 currentElement = SelectRandomAtom(couple,particle,
174 kinEnergy,cutEnergy,kinEnergy);
175
176 G4double Z = currentElement->GetZ();
177 G4int iz = G4int(Z);
178 G4int ia = SelectIsotopeNumber(currentElement);
179 G4double m2 = G4NucleiProperties::GetNuclearMass(ia, iz);
180
181
182
183 G4double xsec= ComputeCrossSectionPerAtom(particle,kinEnergy, Z,
184 kinEnergy, cutEnergy, kinEnergy) ;
185 if(xsec == 0.0)return;
186
187 //scattering angle, z1 == (1-cost)
188 G4double z1 = ioncross->SampleCosineTheta();
189
190 if(z1 <= 0.0) { return; }
191
192 G4double cost = 1.0 - z1;
193 G4double sint = sqrt(z1*(1.0 + cost));
194 G4double phi = twopi * G4UniformRand();
195
196
197 // kinematics in the Lab system
198 G4double etot = kinEnergy + mass;
199 G4double mom2= kinEnergy*(kinEnergy+2.0*mass);
200 G4double ptot = sqrt(mom2);
201
202 //CM particle 1
203 G4double bet = ptot/(etot + m2);
204 G4double gam = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
205
206 //CM
207 G4double momCM2= ioncross->GetMomentum2();
208 G4double momCM =std::sqrt(momCM2);
209 //energy & momentum after scattering of incident particle
210 G4double pxCM = momCM*sint*cos(phi);
211 G4double pyCM = momCM*sint*sin(phi);
212 G4double pzCM = momCM*cost;
213 G4double eCM = sqrt(momCM2 + mass*mass);
214
215 //CM--->Lab
216 G4ThreeVector v1(pxCM , pyCM, gam*(pzCM + bet*eCM));
217 G4ThreeVector dir = dp->GetMomentumDirection();
218
219 G4ThreeVector newDirection = v1.unit();
220 newDirection.rotateUz(dir);
221
222 fParticleChange->ProposeMomentumDirection(newDirection);
223
224 // recoil.......................................
225 G4double trec =(1.0 - cost)* m2*(etot*etot - mass*mass )/
226 (mass*mass + m2*m2+ 2.*m2*etot);
227
228 G4double finalT = kinEnergy - trec;
229
230
231 if(finalT <= lowEnergyLimit) {
232 trec = kinEnergy;
233 finalT = 0.0;
234 }
235
236 fParticleChange->SetProposedKineticEnergy(finalT);
237
238 G4double tcut = recoilThreshold;
239 if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]);
240
241 //G4cout<<" tcut eV "<<tcut/eV<<endl;
242 }
243
244 if(trec > tcut) {
245 G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
246 G4double plab = sqrt(finalT*(finalT + 2.0*mass));
247 G4ThreeVector p2 = (ptot*dir - plab*newDirection).unit();
248 G4DynamicParticle* newdp = new G4DynamicParticle(ion, p2, trec);
249 fvect->push_back(newdp);
250 } else if(trec > 0.0) {
251 fParticleChange->ProposeLocalEnergyDeposit(trec);
252 fParticleChange->ProposeNonIonizingEnergyDeposit(trec);
253 }
254
255
256}
257
258//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
259
Note: See TracBrowser for help on using the repository browser.