source: trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc@ 1330

Last change on this file since 1330 was 1315, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 8.8 KB
Line 
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// $Id: G4eCoulombScatteringModel.cc,v 1.89 2010/05/27 14:22:05 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4eCoulombScatteringModel
35//
36// Author: Vladimir Ivanchenko
37//
38// Creation date: 22.08.2005
39//
40// Modifications:
41//
42// 01.08.06 V.Ivanchenko extend upper limit of table to TeV and review the
43// logic of building - only elements from G4ElementTable
44// 08.08.06 V.Ivanchenko build internal table in ekin scale, introduce faclim
45// 19.08.06 V.Ivanchenko add inline function ScreeningParameter
46// 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e-
47// 09.06.08 V.Ivanchenko add SelectIsotope and sampling of the recoil ion
48// 16.06.09 C.Consolandi fixed computation of effective mass
49// 27.05.10 V.Ivanchenko added G4WentzelOKandVIxSection class to
50// compute cross sections and sample scattering angle
51//
52//
53// Class Description:
54//
55// -------------------------------------------------------------------
56//
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
59
60#include "G4eCoulombScatteringModel.hh"
61#include "Randomize.hh"
62#include "G4DataVector.hh"
63#include "G4ElementTable.hh"
64#include "G4ParticleChangeForGamma.hh"
65#include "G4Proton.hh"
66#include "G4ParticleTable.hh"
67#include "G4ProductionCutsTable.hh"
68#include "G4NucleiProperties.hh"
69#include "G4Pow.hh"
70#include "G4LossTableManager.hh"
71#include "G4NistManager.hh"
72
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74
75using namespace std;
76
77G4eCoulombScatteringModel::G4eCoulombScatteringModel(const G4String& nam)
78 : G4VEmModel(nam),
79 cosThetaMin(1.0),
80 cosThetaMax(-1.0),
81 isInitialised(false)
82{
83 fNistManager = G4NistManager::Instance();
84 theParticleTable = G4ParticleTable::GetParticleTable();
85 theProton = G4Proton::Proton();
86 currentMaterial = 0;
87 currentElement = 0;
88 lowEnergyLimit = 1*eV;
89 recoilThreshold = 0.*keV;
90 particle = 0;
91 currentCouple = 0;
92 wokvi = new G4WentzelOKandVIxSection();
93}
94
95//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96
97G4eCoulombScatteringModel::~G4eCoulombScatteringModel()
98{
99 delete wokvi;
100}
101
102//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103
104void G4eCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
105 const G4DataVector& cuts)
106{
107 SetupParticle(p);
108 currentCouple = 0;
109 cosThetaMin = cos(PolarAngleLimit());
110 wokvi->Initialise(p, cosThetaMin);
111 /*
112 G4cout << "G4eCoulombScatteringModel: factorA2(GeV^2) = " << factorA2/(GeV*GeV)
113 << " 1-cos(ThetaLimit)= " << 1 - cosThetaMin
114 << " cos(thetaMax)= " << cosThetaMax
115 << G4endl;
116 */
117 pCuts = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
118 //G4cout << "!!! G4eCoulombScatteringModel::Initialise for "
119 // << p->GetParticleName() << " cos(TetMin)= " << cosThetaMin
120 // << " cos(TetMax)= " << cosThetaMax <<G4endl;
121 // G4cout << "cut0= " << cuts[0] << " cut1= " << cuts[1] << G4endl;
122 if(!isInitialised) {
123 isInitialised = true;
124 fParticleChange = GetParticleChangeForGamma();
125 }
126 if(mass < GeV && particle->GetParticleType() != "nucleus") {
127 InitialiseElementSelectors(p,cuts);
128 }
129}
130
131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132
133G4double G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(
134 const G4ParticleDefinition* p,
135 G4double kinEnergy,
136 G4double Z, G4double,
137 G4double cutEnergy, G4double)
138{
139 //G4cout << "### G4eCoulombScatteringModel::ComputeCrossSectionPerAtom for "
140 // << p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV << G4endl;
141 G4double xsec = 0.0;
142 if(p != particle) { SetupParticle(p); }
143
144 // cross section is set to zero to avoid problems in sample secondary
145 if(kinEnergy < lowEnergyLimit) { return xsec; }
146 DefineMaterial(CurrentCouple());
147 cosTetMinNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
148 if(cosThetaMax < cosTetMinNuc) {
149 G4int iz = G4int(Z);
150 cosTetMinNuc = wokvi->SetupTarget(iz, cutEnergy);
151 cosTetMaxNuc = cosThetaMax;
152 if(iz == 1 && cosTetMaxNuc < 0.0 && particle == theProton) {
153 cosTetMaxNuc = 0.0;
154 }
155 xsec = wokvi->ComputeNuclearCrossSection(cosTetMinNuc, cosTetMaxNuc);
156 elecRatio = wokvi->ComputeElectronCrossSection(cosTetMinNuc, cosThetaMax);
157 xsec += elecRatio;
158 if(xsec > 0.0) { elecRatio /= xsec; }
159 }
160 /*
161 G4cout << "e(MeV)= " << kinEnergy/MeV << " xsec(b)= " << xsec/barn
162 << " 1-cosTetMinNuc= " << 1-cosTetMinNuc
163 << " 1-cosTetMaxNuc2= " << 1-cosTetMaxNuc2
164 << " 1-cosTetMaxElec= " << 1-cosTetMaxElec
165 << " screenZ= " << screenZ
166 << " formfactA= " << formfactA << G4endl;
167 */
168 return xsec;
169}
170
171//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172
173void G4eCoulombScatteringModel::SampleSecondaries(
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 if(kinEnergy < lowEnergyLimit) { return; }
182 SetupParticle(dp->GetDefinition());
183
184 //G4cout << "G4eCoulombScatteringModel::SampleSecondaries e(MeV)= "
185 // << kinEnergy << " " << particle->GetParticleName()
186 // << " cut= " << cutEnergy<< G4endl;
187
188 // Choose nucleus
189 currentElement = SelectRandomAtom(couple,particle,
190 kinEnergy,cutEnergy,kinEnergy);
191
192 G4double Z = currentElement->GetZ();
193
194 if(ComputeCrossSectionPerAtom(particle,kinEnergy, Z,
195 kinEnergy, cutEnergy, kinEnergy) == 0.0)
196 { return; }
197
198 G4int iz = G4int(Z);
199 G4int ia = SelectIsotopeNumber(currentElement);
200 G4double targetMass = G4NucleiProperties::GetNuclearMass(ia, iz);
201
202 G4ThreeVector newDirection =
203 wokvi->SampleSingleScattering(cosTetMinNuc, cosThetaMax, elecRatio);
204 G4double cost = newDirection.z();
205
206 G4ThreeVector direction = dp->GetMomentumDirection();
207 newDirection.rotateUz(direction);
208
209 fParticleChange->ProposeMomentumDirection(newDirection);
210
211 // recoil sampling assuming a small recoil
212 // and first order correction to primary 4-momentum
213 G4double mom2 = wokvi->GetMomentumSquare();
214 G4double trec = mom2*(1.0 - cost)/(targetMass + (mass + kinEnergy)*(1.0 + cost));
215 G4double finalT = kinEnergy - trec;
216 //G4cout<<"G4eCoulombScatteringModel: finalT= "<<finalT<<" Trec= "<<trec<<G4endl;
217 if(finalT <= lowEnergyLimit) {
218 trec = kinEnergy;
219 finalT = 0.0;
220 }
221
222 fParticleChange->SetProposedKineticEnergy(finalT);
223 G4double tcut = recoilThreshold;
224 if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); }
225
226 if(trec > tcut) {
227 G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
228 G4ThreeVector dir = (direction*sqrt(mom2) -
229 newDirection*sqrt(finalT*(2*mass + finalT))).unit();
230 G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, trec);
231 fvect->push_back(newdp);
232 } else {
233 fParticleChange->ProposeLocalEnergyDeposit(trec);
234 fParticleChange->ProposeNonIonizingEnergyDeposit(trec);
235 }
236
237 return;
238}
239
240//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
241
242
Note: See TracBrowser for help on using the repository browser.