Geant4  10.03.p03
 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 #include "G4ParticleTable.hh"
67 #include "G4IonTable.hh"
68 
69 #include "G4UnitsTable.hh"
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 
73 using namespace std;
74 
76  : G4VEmModel(nam),
77  cosThetaMin(1.0)
78 {
79  fNistManager = G4NistManager::Instance();
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 = 0;
92  mass=0;
93  currentMaterialIndex = -1;
94 
95  ioncross = new G4IonCoulombCrossSection();
96 }
97 
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
102 {
103  delete ioncross;
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 
109  const G4DataVector& cuts)
110 {
111  SetupParticle(p);
112  currentCouple = 0;
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 
126  const G4ParticleDefinition* p,
127  G4double kinEnergy,
128  G4double Z,
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 
148  //cout<< "..........cross "<<G4BestUnit(cross,"Surface") <<endl;
149  return cross;
150 }
151 
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
153 
155  std::vector<G4DynamicParticle*>* fvect,
156  const G4MaterialCutsCouple* couple,
157  const G4DynamicParticle* dp,
159 {
160  G4double kinEnergy = dp->GetKineticEnergy();
161 
162  DefineMaterial(couple);
163 
164  SetupParticle(dp->GetDefinition());
165 
166  // Choose nucleus
167  currentElement = SelectRandomAtom(couple, particle, kinEnergy);
168 
169  G4int iz = currentElement->GetZasInt();
170  G4int ia = SelectIsotopeNumber(currentElement);
172 
173  ioncross->SetupKinematic(kinEnergy, mass2);
174 
175  ioncross->SetupTarget(currentElement->GetZ(), kinEnergy, heavycorr);
176 
177  //scattering angle, z1 == (1-cost)
178  G4double z1 = ioncross->SampleCosineTheta();
179  if(z1 > 2.0) { z1 = 2.0; }
180  else if(z1 < 0.0) { z1 = 0.0; }
181 
182  G4double cost = 1.0 - z1;
183  G4double sint = sqrt(z1*(1.0 + cost));
184  G4double phi = twopi * G4UniformRand();
185 
186  // kinematics in the Lab system
187  G4double ptot = dp->GetTotalMomentum();
188  G4double e1 = dp->GetTotalEnergy();
189 
190  // Lab. system kinematics along projectile direction
191  G4LorentzVector v0 = G4LorentzVector(0, 0, ptot, e1);
192  G4double bet = ptot/(e1 + mass2);
193  G4double gam = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
194 
195  // CM projectile
196  G4double momCM = gam*(ptot - bet*e1);
197  G4double eCM = gam*(e1 - bet*ptot);
198 
199  // Momentum after scattering of incident particle
200  G4double pxCM = momCM*sint*cos(phi);
201  G4double pyCM = momCM*sint*sin(phi);
202  G4double pzCM = momCM*cost;
203 
204  // CM--->Lab
205  G4LorentzVector v1(pxCM , pyCM, gam*(pzCM + bet*eCM), gam*(eCM + bet*pzCM));
206 
207  // Rotate to global system
208  G4ThreeVector dir = dp->GetMomentumDirection();
209  G4ThreeVector newDirection = v1.vect().unit();
210  newDirection.rotateUz(dir);
211 
212  fParticleChange->ProposeMomentumDirection(newDirection);
213 
214  // recoil v0 energy is kinetic
215  v0 -= v1;
216  G4double trec = v0.e();
217  G4double edep = 0.0;
218 
219  G4double tcut = recoilThreshold;
220  if(pCuts) {
221  tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]);
222  //G4cout<<" tcut eV "<<tcut/eV<<endl;
223  }
224 
225  // Recoil
226  if(trec > tcut) {
227  G4ParticleDefinition* ion = theIonTable->GetIon(iz, ia, 0);
228  newDirection = v0.vect().unit();
229  newDirection.rotateUz(dir);
230  G4DynamicParticle* newdp = new G4DynamicParticle(ion, newDirection, trec);
231  fvect->push_back(newdp);
232  } else if(trec > 0.0) {
233  edep = trec;
234  fParticleChange->ProposeNonIonizingEnergyDeposit(edep);
235  }
236 
237  // finelize primary energy and energy balance
238  G4double finalT = v1.e() - mass;
239  if(finalT < 0.0) {
240  edep += finalT;
241  finalT = 0.0;
242  }
243  edep = std::max(edep, 0.0);
244  fParticleChange->SetProposedKineticEnergy(finalT);
245  fParticleChange->ProposeLocalEnergyDeposit(edep);
246 }
247 
248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249 
void SetupTarget(G4double Z, G4double kinEnergy, G4int heavycorr)
static G4double GetNuclearMass(const G4double A, const G4double Z)
void SetupKinematic(G4double kinEnergy, G4double tmass)
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
const char * p
Definition: xmltok.h:285
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
Definition: G4SIunits.hh:76
G4IonCoulombScatteringModel(const G4String &nam="IonCoulombScattering")
G4double GetTotalMomentum() const
G4IonTable * GetIonTable() const
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:454
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
void ProposeNonIonizingEnergyDeposit(G4double anEnergyPart)
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
static G4Proton * Proton()
Definition: G4Proton.cc:93
static constexpr double eV
Definition: G4SIunits.hh:215
float proton_mass_c2
Definition: hepunit.py:275
G4int SelectIsotopeNumber(const G4Element *)
Definition: G4VEmModel.hh:587
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double cut, G4double emax) final
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) final
static G4ParticleTable * GetParticleTable()
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Hep3Vector unit() const
G4double GetAtomicMassAmu(const G4String &symb) const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) final
double G4double
Definition: G4Types.hh:76
float amu_c2
Definition: hepunit.py:277
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
CLHEP::HepLorentzVector G4LorentzVector
void Initialise(const G4ParticleDefinition *, G4double cosThetaLim)