ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eSingleCoulombScatteringModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4eSingleCoulombScatteringModel.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 // G4eSingleCoulombScatteringModel.cc
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class header file
30 //
31 // File name: G4eSingleCoulombScatteringModel
32 //
33 // Author: Cristina Consolandi
34 //
35 // Creation date: 20.10.2012
36 //
37 // Class Description:
38 // Single Scattering model for electron-nuclei interaction.
39 // Suitable for high energy electrons and low scattering angles.
40 //
41 //
42 // Reference:
43 // M.J. Boschini et al. "Non Ionizing Energy Loss induced by Electrons
44 // in the Space Environment" Proc. of the 13th International Conference
45 // on Particle Physics and Advanced Technology
46 //
47 // (13th ICPPAT, Como 3-7/10/2011), World Scientific (Singapore).
48 // Available at: http://arxiv.org/abs/1111.4042v4
49 //
50 //
51 // -------------------------------------------------------------------
52 //
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 
55 
57 #include "G4PhysicalConstants.hh"
58 #include "G4SystemOfUnits.hh"
59 #include "Randomize.hh"
61 #include "G4Proton.hh"
62 #include "G4ProductionCutsTable.hh"
63 #include "G4NucleiProperties.hh"
64 #include "G4NistManager.hh"
65 #include "G4ParticleTable.hh"
66 #include "G4IonTable.hh"
67 
68 #include "G4UnitsTable.hh"
69 #include "G4EmParameters.hh"
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 
73 using namespace std;
74 
76  : G4VEmModel(nam),
77  cosThetaMin(1.0)
78 {
81  fParticleChange = nullptr;
82 
83  pCuts=nullptr;
84  currentMaterial = nullptr;
85  currentElement = nullptr;
86  currentCouple = nullptr;
87 
88  lowEnergyLimit = 0*keV;
89  recoilThreshold = 0.*eV;
90  XSectionModel = 1;
91  FormFactor = 0;
92  particle = nullptr;
93  mass=0.0;
95 
97  //G4cout <<"## G4eSingleCoulombScatteringModel: " << this << " " << Mottcross << G4endl;
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 
103 {
104  //G4cout <<"## G4eSingleCoulombScatteringModel: delete " << this << " " << Mottcross << G4endl;
105  delete Mottcross;
106 }
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109 
111  const G4DataVector& cuts)
112 {
114 
115  SetupParticle(p);
116  currentCouple = nullptr;
118  //cosThetaMin = cos(PolarAngleLimit());
120 
121  pCuts = &cuts;
122  //G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
123 
124  /*
125  G4cout << "!!! G4eSingleCoulombScatteringModel::Initialise for "
126  << part->GetParticleName() << " cos(TetMin)= " << cosThetaMin
127  << " cos(TetMax)= " << cosThetaMax <<G4endl;
128  G4cout << "cut= " << (*pCuts)[0] << " cut1= " << (*pCuts)[1] << G4endl;
129  */
130 
131  if(!fParticleChange) {
133  }
134 
135  if(IsMaster()) {
137  }
138 
140 
141  //G4cout<<"NUCLEAR FORM FACTOR: "<<FormFactor<<G4endl;
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145 
146 void
148  G4VEmModel* masterModel)
149 {
151 }
152 
153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
154 
156 {
157  if(model == "Fast" || model == "fast") { XSectionModel=1; }
158  else if(model == "Precise" || model == "precise") { XSectionModel=0; }
159  else {
160  G4cout<<"G4eSingleCoulombScatteringModel WARNING: "<<model
161  <<" is not a valid model name"<<G4endl;
162  }
163 }
164 
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
166 
168  const G4ParticleDefinition* p,
169  G4double kinEnergy,
170  G4double Z,
171  G4double ,
172  G4double,
173  G4double )
174 {
175  SetupParticle(p);
176 
177  G4double cross =0.0;
178  if(kinEnergy < lowEnergyLimit) return cross;
179 
181 
182  //Total Cross section
183  Mottcross->SetupKinematic(kinEnergy, Z);
185 
186  //cout<< "Compute Cross Section....cross "<<G4BestUnit(cross,"Surface") << " cm2 "<< cross/cm2 <<" Z: "<<Z<<" kinEnergy: "<<kinEnergy<<endl;
187 
188  //G4cout<<"Energy: "<<kinEnergy/MeV<<" Total Cross: "<<cross<<G4endl;
189  return cross;
190 }
191 
192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
193 
195  std::vector<G4DynamicParticle*>* fvect,
196  const G4MaterialCutsCouple* couple,
197  const G4DynamicParticle* dp,
198  G4double cutEnergy,
199  G4double)
200 {
201  G4double kinEnergy = dp->GetKineticEnergy();
202  //cout<<"--- kinEnergy "<<kinEnergy<<endl;
203 
204  if(kinEnergy < lowEnergyLimit) return;
205 
206  DefineMaterial(couple);
208 
209  // Choose nucleus
210  //last two :cutEnergy= min e kinEnergy=max
211  currentElement = SelectTargetAtom(couple, particle, kinEnergy,
212  dp->GetLogKineticEnergy(), cutEnergy, kinEnergy);
213  G4int iz = currentElement->GetZasInt();
214  G4int ia = SelectIsotopeNumber(currentElement);
216 
217  //G4cout<<"..Z: "<<Z<<" ..iz: "<<iz<<" ..ia: "<<ia<<" ..mass2: "<<mass2<<G4endl;
218 
219  Mottcross->SetupKinematic(kinEnergy, iz);
221  if(cross == 0.0) { return; }
222  //cout<< "Energy: "<<kinEnergy/MeV<<" Z: "<<Z<<"....cross "<<G4BestUnit(cross,"Surface") << " cm2 "<< cross/cm2 <<endl;
223 
225  G4double sint = sin(z1);
226  G4double cost = cos(z1);
228 
229  // kinematics in the Lab system
230  G4double ptot = sqrt(kinEnergy*(kinEnergy + 2.0*mass));
231  G4double e1 = mass + kinEnergy;
232 
233  // Lab. system kinematics along projectile direction
234  G4LorentzVector v0 = G4LorentzVector(0, 0, ptot, e1+mass2);
235  G4LorentzVector v1 = G4LorentzVector(0, 0, ptot, e1);
236  G4ThreeVector bst = v0.boostVector();
237  v1.boost(-bst);
238  // CM projectile
239  G4double momCM = v1.pz();
240 
241  // Momentum after scattering of incident particle
242  v1.setX(momCM*sint*cos(phi));
243  v1.setY(momCM*sint*sin(phi));
244  v1.setZ(momCM*cost);
245 
246  // CM--->Lab
247  v1.boost(bst);
248 
249  // Rotate to global system
251  G4ThreeVector newDirection = v1.vect().unit();
252  newDirection.rotateUz(dir);
253 
255 
256  // recoil
257  v0 -= v1;
258  G4double trec = std::max(v0.e() - mass2, 0.0);
259  G4double edep = 0.0;
260 
261  G4double tcut = recoilThreshold;
262 
263  //G4cout<<" Energy Transfered: "<<trec/eV<<G4endl;
264 
265  if(pCuts) {
266  tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]);
267  //G4cout<<"Cuts: "<<(*pCuts)[currentMaterialIndex]/eV<<" eV"<<G4endl;
268  //G4cout<<"Threshold: "<<tcut/eV<<" eV"<<G4endl;
269  }
270 
271  if(trec > tcut) {
273  newDirection = v0.vect().unit();
274  newDirection.rotateUz(dir);
275  G4DynamicParticle* newdp = new G4DynamicParticle(ion, newDirection, trec);
276  fvect->push_back(newdp);
277  } else if(trec > 0.0) {
278  edep = trec;
280  }
281 
282  // finelize primary energy and energy balance
283  G4double finalT = v1.e() - mass;
284  //G4cout<<"Final Energy: "<<finalT/eV<<G4endl;
285  if(finalT <= lowEnergyLimit) {
286  edep += finalT;
287  finalT = 0.0;
288  }
289  edep = std::max(edep, 0.0);
292 
293 }
294 
295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......