Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eSingleCoulombScatteringModel.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 // 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 {
79  fNistManager = G4NistManager::Instance();
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  FormFactor = 0;
91  particle = nullptr;
92  mass=0.0;
93  currentMaterialIndex = -1;
94 
95  Mottcross = new G4ScreeningMottCrossSection();
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99 
101 {
102  delete Mottcross;
103 }
104 
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
106 
108  const G4DataVector& cuts)
109 {
111 
112  SetupParticle(p);
113  currentCouple = nullptr;
114  currentMaterialIndex = -1;
115  //cosThetaMin = cos(PolarAngleLimit());
116  Mottcross->Initialise(p,cosThetaMin);
117 
118  pCuts = &cuts;
119  //G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
120 
121  /*
122  G4cout << "!!! G4eSingleCoulombScatteringModel::Initialise for "
123  << part->GetParticleName() << " cos(TetMin)= " << cosThetaMin
124  << " cos(TetMax)= " << cosThetaMax <<G4endl;
125  G4cout << "cut= " << (*pCuts)[0] << " cut1= " << (*pCuts)[1] << G4endl;
126  */
127 
128  if(!fParticleChange) {
129  fParticleChange = GetParticleChangeForGamma();
130  }
131 
132  if(IsMaster()) {
134  }
135 
136  FormFactor=param->NuclearFormfactorType();
137 
138  //G4cout<<"NUCLEAR FORM FACTOR: "<<FormFactor<<G4endl;
139 }
140 
142  G4VEmModel* masterModel)
143 {
145 }
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147 
149  const G4ParticleDefinition* p,
150  G4double kinEnergy,
151  G4double Z,
152  G4double ,
153  G4double,
154  G4double )
155 {
156  SetupParticle(p);
157 
158  G4double cross =0.0;
159  if(kinEnergy < lowEnergyLimit) return cross;
160 
161  DefineMaterial(CurrentCouple());
162 
163  //Total Cross section
164  Mottcross->SetupKinematic(kinEnergy, Z);
165  cross = Mottcross->NuclearCrossSection(FormFactor);
166 
167  //cout<< "Compute Cross Section....cross "<<G4BestUnit(cross,"Surface") << " cm2 "<< cross/cm2 <<" Z: "<<Z<<" kinEnergy: "<<kinEnergy<<endl;
168  return cross;
169 }
170 
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172 
174  std::vector<G4DynamicParticle*>* fvect,
175  const G4MaterialCutsCouple* couple,
176  const G4DynamicParticle* dp,
177  G4double cutEnergy,
178  G4double)
179 {
180  G4double kinEnergy = dp->GetKineticEnergy();
181  //cout<<"--- kinEnergy "<<kinEnergy<<endl;
182 
183  if(kinEnergy < lowEnergyLimit) return;
184 
185  DefineMaterial(couple);
186  SetupParticle(dp->GetDefinition());
187 
188  // Choose nucleus
189  //last two :cutEnergy= min e kinEnergy=max
190  currentElement = SelectRandomAtom(couple,particle,
191  kinEnergy,cutEnergy,kinEnergy);
192 
193  G4double Z = currentElement->GetZ();
194  G4int iz = G4int(Z);
195  G4int ia = SelectIsotopeNumber(currentElement);
197 
198  //G4cout<<"..Z: "<<Z<<" ..iz: "<<iz<<" ..ia: "<<ia<<" ..mass2: "<<mass2<<G4endl;
199 
200  Mottcross->SetupKinematic(kinEnergy, Z);
201  G4double cross= Mottcross->NuclearCrossSection(FormFactor); //MODIFY TO LOAD TABLE
202  if(cross == 0.0) { return; }
203  //cout<< "Energy: "<<kinEnergy/MeV<<" Z: "<<Z<<"....cross "<<G4BestUnit(cross,"Surface") << " cm2 "<< cross/cm2 <<endl;
204 
205  G4double z1 = Mottcross->GetScatteringAngle();
206  G4double sint = sin(z1);
207  G4double cost = sqrt(1.0 - sint*sint);
208  G4double phi = twopi* G4UniformRand();
209 
210  // kinematics in the Lab system
211  G4double ptot = dp->GetTotalMomentum();
212  G4double e1 = dp->GetTotalEnergy();
213  // Lab. system kinematics along projectile direction
214  G4LorentzVector v0 = G4LorentzVector(0, 0, ptot, e1);
215  G4double bet = ptot/(v0.e() + mass2);
216  G4double gam = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
217 
218  //CM Projectile
219  G4double momCM = gam*(ptot - bet*e1);
220  G4double eCM = gam*(e1 - bet*ptot);
221  //energy & momentum after scattering of incident particle
222  G4double pxCM = momCM*sint*cos(phi);
223  G4double pyCM = momCM*sint*sin(phi);
224  G4double pzCM = momCM*cost;
225 
226  //CM--->Lab
227  G4LorentzVector v1(pxCM , pyCM, gam*(pzCM + bet*eCM), gam*(eCM + bet*pzCM));
228 
229  // Rotate to global system
231  G4ThreeVector newDirection = v1.vect().unit();
232  newDirection.rotateUz(dir);
233 
234  fParticleChange->ProposeMomentumDirection(newDirection);
235 
236  // recoil
237  v0 -= v1;
238  G4double trec = v0.e();
239  G4double edep = 0.0;
240 
241  G4double tcut = recoilThreshold;
242 
243  //G4cout<<" Energy Transfered: "<<trec/eV<<G4endl;
244 
245  if(pCuts) {
246  tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]);
247  }
248 
249  if(trec > tcut) {
250  G4ParticleDefinition* ion = theIonTable->GetIon(iz, ia, 0);
251  newDirection = v0.vect().unit();
252  newDirection.rotateUz(dir);
253  G4DynamicParticle* newdp = new G4DynamicParticle(ion, newDirection, trec);
254  fvect->push_back(newdp);
255  } else if(trec > 0.0) {
256  edep = trec;
257  fParticleChange->ProposeNonIonizingEnergyDeposit(edep);
258  }
259 
260  // finelize primary energy and energy balance
261  G4double finalT = v1.e() - mass;
262  //G4cout<<"Final Energy: "<<finalT/eV<<G4endl;
263  if(finalT <= lowEnergyLimit) {
264  edep += finalT;
265  finalT = 0.0;
266  }
267  edep = std::max(edep, 0.0);
268  fParticleChange->SetProposedKineticEnergy(finalT);
269  fParticleChange->ProposeLocalEnergyDeposit(edep);
270 
271 }
272 
273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) final
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetKineticEnergy() const
void Initialise(const G4ParticleDefinition *, G4double cosThetaLim)
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
G4double GetTotalEnergy() const
const char * p
Definition: xmltok.h:285
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
G4bool IsMaster() const
Definition: G4VEmModel.hh:725
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) final
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double GetTotalMomentum() const
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double cut, G4double emax) final
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 constexpr double eV
Definition: G4SIunits.hh:215
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) final
G4int SelectIsotopeNumber(const G4Element *)
Definition: G4VEmModel.hh:587
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:809
G4NuclearFormfactorType NuclearFormfactorType() const
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:817
G4eSingleCoulombScatteringModel(const G4String &nam="eSingleCoulombScat")
static G4ParticleTable * GetParticleTable()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Hep3Vector unit() const
static G4EmParameters * Instance()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
double G4double
Definition: G4Types.hh:76
void SetupKinematic(G4double kinEnergy, G4double Z)
static constexpr double keV
Definition: G4SIunits.hh:216
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