Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4eSingleCoulombScatteringModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 
 56 #include "G4eSingleCoulombScatteringModel.hh"
 57 #include "G4PhysicalConstants.hh"
 58 #include "G4SystemOfUnits.hh"
 59 #include "Randomize.hh"
 60 #include "G4ParticleChangeForGamma.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 
 75 G4eSingleCoulombScatteringModel::G4eSingleCoulombScatteringModel(const G4String& nam)
 76   : G4VEmModel(nam),
 77     cosThetaMin(1.0)
 78 {
 79   fNistManager = G4NistManager::Instance();
 80   theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
 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;
 94   currentMaterialIndex = -1;
 95 
 96   Mottcross = new G4ScreeningMottCrossSection();
 97   //G4cout <<"## G4eSingleCoulombScatteringModel: " << this << "  " << Mottcross << G4endl;
 98 }
 99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 
102 G4eSingleCoulombScatteringModel::~G4eSingleCoulombScatteringModel()
103 {
104   //G4cout <<"## G4eSingleCoulombScatteringModel: delete " << this << "  " << Mottcross << G4endl;
105   delete Mottcross;
106 }
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109 
110 void G4eSingleCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
111              const G4DataVector&  cuts)
112 {
113   G4EmParameters* param = G4EmParameters::Instance();
114 
115   SetupParticle(p);
116   currentCouple = nullptr;
117   currentMaterialIndex = -1;
118   //cosThetaMin = cos(PolarAngleLimit());
119   Mottcross->Initialise(p,cosThetaMin);
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) {
132     fParticleChange = GetParticleChangeForGamma();
133   }
134 
135   if(IsMaster()) {
136     InitialiseElementSelectors(p,cuts);
137   }
138 
139   FormFactor=param->NuclearFormfactorType();
140 
141   //G4cout<<"NUCLEAR FORM FACTOR: "<<FormFactor<<G4endl;
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145 
146 void 
147 G4eSingleCoulombScatteringModel::InitialiseLocal(const G4ParticleDefinition*,
148                                                  G4VEmModel* masterModel)
149 {
150   SetElementSelectors(masterModel->GetElementSelectors());
151 }
152 
153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
154 
155 void G4eSingleCoulombScatteringModel::SetXSectionModel(const G4String& model)
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 
167 G4double G4eSingleCoulombScatteringModel::ComputeCrossSectionPerAtom(
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 
180   DefineMaterial(CurrentCouple());
181 
182   //Total Cross section
183   Mottcross->SetupKinematic(kinEnergy, Z);
184   cross = Mottcross->NuclearCrossSection(FormFactor,XSectionModel);
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 
194 void G4eSingleCoulombScatteringModel::SampleSecondaries(
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);
207   SetupParticle(dp->GetDefinition());
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);
215   G4double mass2 = G4NucleiProperties::GetNuclearMass(ia, iz);
216 
217   //G4cout<<"..Z: "<<Z<<" ..iz: "<<iz<<" ..ia: "<<ia<<" ..mass2: "<<mass2<<G4endl;
218 
219   Mottcross->SetupKinematic(kinEnergy, iz);
220   G4double cross= Mottcross->NuclearCrossSection(FormFactor,XSectionModel);
221   if(cross == 0.0) { return; }
222   //cout<< "Energy: "<<kinEnergy/MeV<<" Z: "<<Z<<"....cross "<<G4BestUnit(cross,"Surface") << " cm2 "<< cross/cm2 <<endl;
223 
224   G4double z1 = Mottcross->GetScatteringAngle(FormFactor,XSectionModel);
225   G4double sint = sin(z1);
226   G4double cost = cos(z1);
227   G4double phi  = twopi* G4UniformRand();
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
250   G4ThreeVector dir = dp->GetMomentumDirection();
251   G4ThreeVector newDirection = v1.vect().unit();
252   newDirection.rotateUz(dir);
253 
254   fParticleChange->ProposeMomentumDirection(newDirection);
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) {
272     G4ParticleDefinition* ion = theIonTable->GetIon(iz, ia, 0);
273     newDirection = v0.vect().unit();
274     newDirection.rotateUz(dir);
275     auto newdp  = new G4DynamicParticle(ion, newDirection, trec);
276     fvect->push_back(newdp);
277   } else if(trec > 0.0) {
278     edep = trec;
279     fParticleChange->ProposeNonIonizingEnergyDeposit(edep);
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);
290   fParticleChange->SetProposedKineticEnergy(finalT);
291   fParticleChange->ProposeLocalEnergyDeposit(edep);
292 
293 }
294 
295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
296