ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4IonCoulombScatteringModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4IonCoulombScatteringModel.cc
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 {
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;
94 
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;
115 
116  pCuts = &cuts;
117  // G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
118  if(!fParticleChange) {
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 
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  return cross;
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 
153  std::vector<G4DynamicParticle*>* fvect,
154  const G4MaterialCutsCouple* couple,
155  const G4DynamicParticle* dp,
157 {
158  G4double kinEnergy = dp->GetKineticEnergy();
159  DefineMaterial(couple);
161 
162  // Choose nucleus
163  currentElement = SelectTargetAtom(couple, particle, kinEnergy,
164  dp->GetLogKineticEnergy());
165 
166  G4int iz = currentElement->GetZasInt();
167  G4int ia = SelectIsotopeNumber(currentElement);
169 
170  ioncross->SetupKinematic(kinEnergy, mass2);
171  ioncross->SetupTarget(currentElement->GetZ(), kinEnergy, heavycorr);
172 
173  //scattering angle, z1 == (1-cost)
175  if(z1 > 2.0) { z1 = 2.0; }
176  else if(z1 < 0.0) { z1 = 0.0; }
177  /*
178  G4cout << "Sample: " << particle->GetParticleName()
179  << " mass(GeV)= " << mass/GeV
180  << " Ekin(MeV)= " << kinEnergy << " cost= " << 1. - z1 << G4endl;
181  G4cout << " Z= " << iz << " A= " << ia
182  << " mass(GeV)= " << mass2/GeV << G4endl;
183  */
184  G4double cost = 1.0 - z1;
185  G4double sint = sqrt(z1*(1.0 + cost));
187 
188  // kinematics in the Lab system
189  G4double ptot = sqrt(kinEnergy*(kinEnergy + 2.0*mass));
190  G4double e1 = mass + kinEnergy;
191 
192  // Lab. system kinematics along projectile direction
193  G4LorentzVector v0 = G4LorentzVector(0, 0, ptot, e1+mass2);
194  G4LorentzVector v1 = G4LorentzVector(0, 0, ptot, e1);
195  G4ThreeVector bst = v0.boostVector();
196  v1.boost(-bst);
197  // CM projectile
198  G4double momCM = v1.pz();
199 
200  // Momentum after scattering of incident particle
201  v1.setX(momCM*sint*cos(phi));
202  v1.setY(momCM*sint*sin(phi));
203  v1.setZ(momCM*cost);
204 
205  // CM--->Lab
206  v1.boost(bst);
207 
208  // Rotate to global system
210  G4ThreeVector newDirection = v1.vect().unit();
211  newDirection.rotateUz(dir);
212 
214 
215  // recoil v0 energy is kinetic
216  v0 -= v1;
217  G4double trec = std::max(v0.e() - mass2, 0.0);
218  G4double edep = 0.0;
219 
220  G4double tcut = recoilThreshold;
221  if(pCuts) {
222  tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]);
223  //G4cout<<" tcut eV "<<tcut/eV<<endl;
224  }
225 
226  // Recoil
227  if(trec > tcut) {
229  newDirection = v0.vect().unit();
230  newDirection.rotateUz(dir);
231  G4DynamicParticle* newdp = new G4DynamicParticle(ion, newDirection, trec);
232  fvect->push_back(newdp);
233  } else if(trec > 0.0) {
234  edep = trec;
236  }
237 
238  // finelize primary energy and energy balance
239  G4double finalT = v1.e() - mass;
240  if(finalT < 0.0) {
241  edep += finalT;
242  finalT = 0.0;
243  }
244  edep = std::max(edep, 0.0);
245  //G4cout << "Efinal(MeV)= " << finalT << " Edep(MeV)= " << edep
246  // << " Trec(MeV)= " << trec << G4endl;
249 }
250 
251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
252