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

Last change on this file since 1348 was 1340, checked in by garnier, 14 years ago

update ti head
File size: 8.9 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.90 2010/10/26 10:35:22 vnivanch Exp $
27// GEANT4 tag $Name: emstand-V09-03-24 $
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  currentMaterialIndex = 0;
95
96  cosTetMinNuc = 1.0;
97  cosTetMaxNuc = -1.0;
98  elecRatio = 0.0;
99  mass = proton_mass_c2;
100}
101
102//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103
104G4eCoulombScatteringModel::~G4eCoulombScatteringModel()
105{
106  delete wokvi;
107}
108
109//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110
111void G4eCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
112                                           const G4DataVector& cuts)
113{
114  SetupParticle(p);
115  currentCouple = 0;
116  cosThetaMin = cos(PolarAngleLimit());
117  wokvi->Initialise(p, cosThetaMin);
118  /*
119  G4cout << "G4eCoulombScatteringModel: factorA2(GeV^2) = " << factorA2/(GeV*GeV)
120         << "  1-cos(ThetaLimit)= " << 1 - cosThetaMin
121         << "  cos(thetaMax)= " <<  cosThetaMax
122         << G4endl;
123  */
124  pCuts = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
125  //G4cout << "!!! G4eCoulombScatteringModel::Initialise for "
126  //     << p->GetParticleName() << "  cos(TetMin)= " << cosThetaMin
127  //     << "  cos(TetMax)= " << cosThetaMax <<G4endl;
128  // G4cout << "cut0= " << cuts[0] << "  cut1= " << cuts[1] << G4endl;
129  if(!isInitialised) {
130    isInitialised = true;
131    fParticleChange = GetParticleChangeForGamma();
132  }
133  if(mass < GeV && particle->GetParticleType() != "nucleus") {
134    InitialiseElementSelectors(p,cuts);
135  }
136}
137
138//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
139
140G4double G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(
141                const G4ParticleDefinition* p,
142                G4double kinEnergy,
143                G4double Z, G4double,
144                G4double cutEnergy, G4double)
145{
146  //G4cout << "### G4eCoulombScatteringModel::ComputeCrossSectionPerAtom  for "
147  //  << p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV << G4endl;
148  G4double xsec = 0.0;
149  if(p != particle) { SetupParticle(p); }
150
151  // cross section is set to zero to avoid problems in sample secondary
152  if(kinEnergy < lowEnergyLimit) { return xsec; }
153  DefineMaterial(CurrentCouple());
154  cosTetMinNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
155  if(cosThetaMax < cosTetMinNuc) {
156    G4int iz = G4int(Z);
157    cosTetMinNuc = wokvi->SetupTarget(iz, cutEnergy);
158    cosTetMaxNuc = cosThetaMax; 
159    if(iz == 1 && cosTetMaxNuc < 0.0 && particle == theProton) { 
160      cosTetMaxNuc = 0.0; 
161    }
162    xsec =  wokvi->ComputeNuclearCrossSection(cosTetMinNuc, cosTetMaxNuc);
163    elecRatio = wokvi->ComputeElectronCrossSection(cosTetMinNuc, cosThetaMax);
164    xsec += elecRatio;
165    if(xsec > 0.0) { elecRatio /= xsec; } 
166  }
167  /*
168  G4cout << "e(MeV)= " << kinEnergy/MeV << " xsec(b)= " << xsec/barn 
169         << " 1-cosTetMinNuc= " << 1-cosTetMinNuc
170         << " 1-cosTetMaxNuc2= " << 1-cosTetMaxNuc2
171         << " 1-cosTetMaxElec= " << 1-cosTetMaxElec
172         << " screenZ= " << screenZ
173         << " formfactA= " << formfactA << G4endl;
174  */
175  return xsec; 
176}
177
178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179
180void G4eCoulombScatteringModel::SampleSecondaries(
181                std::vector<G4DynamicParticle*>* fvect,
182                const G4MaterialCutsCouple* couple,
183                const G4DynamicParticle* dp,
184                G4double cutEnergy,
185                G4double)
186{
187  G4double kinEnergy = dp->GetKineticEnergy();
188  if(kinEnergy < lowEnergyLimit) { return; }
189  SetupParticle(dp->GetDefinition());
190
191  //G4cout << "G4eCoulombScatteringModel::SampleSecondaries e(MeV)= "
192  //     << kinEnergy << "  " << particle->GetParticleName()
193  //     << " cut= " << cutEnergy<< G4endl;
194 
195  // Choose nucleus
196  currentElement = SelectRandomAtom(couple,particle,
197                                    kinEnergy,cutEnergy,kinEnergy);
198
199  G4double Z = currentElement->GetZ();
200 
201  if(ComputeCrossSectionPerAtom(particle,kinEnergy, Z,
202                                kinEnergy, cutEnergy, kinEnergy) == 0.0) 
203    { return; }
204
205  G4int iz = G4int(Z);
206  G4int ia = SelectIsotopeNumber(currentElement);
207  G4double targetMass = G4NucleiProperties::GetNuclearMass(ia, iz);
208
209  G4ThreeVector newDirection = 
210    wokvi->SampleSingleScattering(cosTetMinNuc, cosThetaMax, elecRatio);
211  G4double cost = newDirection.z();
212
213  G4ThreeVector direction = dp->GetMomentumDirection(); 
214  newDirection.rotateUz(direction);   
215
216  fParticleChange->ProposeMomentumDirection(newDirection);   
217
218  // recoil sampling assuming a small recoil
219  // and first order correction to primary 4-momentum
220  G4double mom2 = wokvi->GetMomentumSquare();
221  G4double trec = mom2*(1.0 - cost)/(targetMass + (mass + kinEnergy)*(1.0 + cost));
222  G4double finalT = kinEnergy - trec; 
223  //G4cout<<"G4eCoulombScatteringModel: finalT= "<<finalT<<" Trec= "<<trec<<G4endl;
224  if(finalT <= lowEnergyLimit) { 
225    trec = kinEnergy; 
226    finalT = 0.0;
227  } 
228
229  fParticleChange->SetProposedKineticEnergy(finalT);
230  G4double tcut = recoilThreshold;
231  if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); }
232
233  if(trec > tcut) {
234    G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
235    G4ThreeVector dir = (direction*sqrt(mom2) - 
236                         newDirection*sqrt(finalT*(2*mass + finalT))).unit();
237    G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, trec);
238    fvect->push_back(newdp);
239  } else {
240    fParticleChange->ProposeLocalEnergyDeposit(trec);
241    fParticleChange->ProposeNonIonizingEnergyDeposit(trec);
242  }
243 
244  return;
245}
246
247//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
248
249
Note: See TracBrowser for help on using the repository browser.