Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4IonCoulombScatteringModel.cc
Go to the documentation of this file.
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 
59 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh"
61 #include "Randomize.hh"
63 #include "G4Proton.hh"
64 #include "G4ProductionCutsTable.hh"
65 #include "G4NucleiProperties.hh"
66 
67 #include "G4UnitsTable.hh"
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70 
71 using namespace std;
72 
74  : G4VEmModel(nam),
75 
76  cosThetaMin(1.0),
77  isInitialised(false)
78 {
82 
83  pCuts=0;
84  currentMaterial = 0;
85  currentElement = 0;
86  currentCouple = 0;
87 
88  lowEnergyLimit = 100*eV;
89  recoilThreshold = 0.*eV;
90  heavycorr =0;
91  particle = 0;
92  mass=0;
94 
96 
97 }
98 
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 
103 { delete ioncross;}
104 
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
106 
108  const G4DataVector& )
109 {
110  SetupParticle(p);
111  currentCouple = 0;
113  cosThetaMin = cos(PolarAngleLimit());
115 
117 
118 
119  if(!isInitialised) {
120  isInitialised = true;
122  }
123 }
124 
125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
126 
128  const G4ParticleDefinition* p,
129  G4double kinEnergy,
130  G4double Z,
131  G4double,
132  G4double cutEnergy,
133  G4double)
134 {
135 
136  SetupParticle(p);
137 
138  G4double cross =0.0;
139  if(kinEnergy < lowEnergyLimit) return cross;
140 
142 
143  G4int iz = G4int(Z);
144 
145  //from lab to pCM & mu_rel of effective particle
146  ioncross->SetupKinematic(kinEnergy, cutEnergy,iz);
147 
148 
149  ioncross->SetupTarget(Z, kinEnergy, heavycorr);
150 
151  cross = ioncross->NuclearCrossSection();
152 
153 //cout<< "..........cross "<<G4BestUnit(cross,"Surface") <<endl;
154  return cross;
155 }
156 
157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158 
160  std::vector<G4DynamicParticle*>* fvect,
161  const G4MaterialCutsCouple* couple,
162  const G4DynamicParticle* dp,
163  G4double cutEnergy,
164  G4double)
165 {
166  G4double kinEnergy = dp->GetKineticEnergy();
167 
168  if(kinEnergy < lowEnergyLimit) return;
169 
170  DefineMaterial(couple);
171 
173 
174  // Choose nucleus
176  kinEnergy,cutEnergy,kinEnergy);
177 
178  G4double Z = currentElement->GetZ();
179  G4int iz = G4int(Z);
180  G4int ia = SelectIsotopeNumber(currentElement);
182 
183 
184 
185  G4double cross= ComputeCrossSectionPerAtom(particle,kinEnergy, Z,
186  kinEnergy, cutEnergy, kinEnergy) ;
187  if(cross == 0.0) { return; }
188 
189  //scattering angle, z1 == (1-cost)
191  if(z1 > 2.0) { z1 = 2.0; }
192  else if(z1 < 0.0) { z1 = 0.0; }
193 
194  G4double cost = 1.0 - z1;
195  G4double sint = sqrt(z1*(1.0 + cost));
196  G4double phi = twopi * G4UniformRand();
197 
198 
199  // kinematics in the Lab system
200  G4double etot = kinEnergy + mass;
201  G4double mom2= kinEnergy*(kinEnergy+2.0*mass);
202  G4double ptot = sqrt(mom2);
203 
204  //CM particle 1
205  G4double bet = ptot/(etot + mass2);
206  G4double gam = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
207 
208  //CM
209  G4double momCM2= ioncross->GetMomentum2();
210  G4double momCM =std::sqrt(momCM2);
211  //energy & momentum after scattering of incident particle
212  G4double pxCM = momCM*sint*cos(phi);
213  G4double pyCM = momCM*sint*sin(phi);
214  G4double pzCM = momCM*cost;
215  G4double eCM = sqrt(momCM2 + mass*mass);
216 
217  //CM--->Lab
218  G4ThreeVector v1(pxCM , pyCM, gam*(pzCM + bet*eCM));
220 
221  G4ThreeVector newDirection = v1.unit();
222  newDirection.rotateUz(dir);
223 
225 
226  // V.Ivanchenko fix of final energies after scattering
227  // recoil.......................................
228  //G4double trec =(1.0 - cost)* mass2*(etot*etot - mass*mass )/
229  // (mass*mass + mass2*mass2+ 2.*mass2*etot);
230  //G4double finalT = kinEnergy - trec;
231 
232  // new computation
233  G4double finalT = gam*(eCM + bet*pzCM) - mass;
234  G4double trec = kinEnergy - finalT;
235 
236  if(finalT <= lowEnergyLimit) {
237  trec = kinEnergy;
238  finalT = 0.0;
239  } else if(trec < 0.0) {
240  trec = 0.0;
241  finalT = kinEnergy;
242  }
243 
245 
246  G4double tcut = recoilThreshold;
247  if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]);
248 
249  //G4cout<<" tcut eV "<<tcut/eV<<endl;
250  }
251 
252  if(trec > tcut) {
254  G4double plab = sqrt(finalT*(finalT + 2.0*mass));
255  G4ThreeVector p2 = (ptot*dir - plab*newDirection).unit();
256  G4DynamicParticle* newdp = new G4DynamicParticle(ion, p2, trec);
257  fvect->push_back(newdp);
258  } else if(trec > 0.0) {
261  }
262 
263 
264 }
265 
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
267