source: trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc @ 846

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 9.1 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: G4CoulombScatteringModel.cc,v 1.29 2007/11/09 11:45:45 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:     G4CoulombScatteringModel
35//
36// Author:        Vladimir Ivanchenko
37//
38// Creation date: 22.08.2005
39//
40// Modifications:
41// 01.08.06 V.Ivanchenko extend upper limit of table to TeV and review the
42//          logic of building - only elements from G4ElementTable
43// 08.08.06 V.Ivanchenko build internal table in ekin scale, introduce faclim
44// 19.10.06 V.Ivanchenko use inheritance from G4eCoulombScatteringModel
45// 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e-
46//
47// Class Description:
48//
49// -------------------------------------------------------------------
50//
51//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53
54#include "G4CoulombScatteringModel.hh"
55#include "Randomize.hh"
56#include "G4ParticleChangeForGamma.hh"
57#include "G4NistManager.hh"
58#include "G4ParticleTable.hh"
59#include "G4IonTable.hh"
60#include "G4Proton.hh"
61
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63
64using namespace std;
65
66G4CoulombScatteringModel::G4CoulombScatteringModel(
67  G4double thetaMin, G4double thetaMax, G4bool build, 
68  G4double tlim, const G4String& nam)
69  : G4eCoulombScatteringModel(thetaMin,thetaMax,build,tlim,nam)
70{
71  theMatManager    = G4NistManager::Instance();
72  theParticleTable = G4ParticleTable::GetParticleTable();
73}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77G4CoulombScatteringModel::~G4CoulombScatteringModel()
78{}
79
80//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81
82G4double G4CoulombScatteringModel::ComputeCrossSectionPerAtom(
83                                const G4ParticleDefinition* p,
84                                G4double kinEnergy, 
85                                G4double Z, 
86                                G4double A, 
87                                G4double cutEnergy,
88                                G4double)
89{
90  if(p == particle && kinEnergy == tkin && Z == targetZ &&
91     A == targetA && cutEnergy == ecut) return nucXSection;
92
93  // Lab system
94  G4double ekin = std::max(keV, kinEnergy);
95  nucXSection = ComputeElectronXSectionPerAtom(p,ekin,Z,A,cutEnergy);
96
97  // CM system
98  G4int iz      = G4int(Z);
99  G4double m1   = theMatManager->GetAtomicMassAmu(iz)*amu_c2;
100  G4double etot = tkin + mass;
101  G4double ptot = sqrt(mom2);
102  G4double bet  = ptot/(etot + m1);
103  G4double gam  = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
104  G4double momCM= gam*(ptot - bet*etot);
105
106  //  G4cout << "ptot= " << ptot << " etot= " << etot << " beta= "
107  //     << bet << " gam= " << gam << " Z= " << Z << " A= " << A << G4endl;
108  // G4cout << " CM. mom= " << momCM  << " m= " << m
109  // << " m1= " << m1 << " iz= " << iz <<G4endl;
110
111  G4double momCM2 = momCM*momCM;
112  cosTetMaxNuc = std::max(cosThetaMax, 1.0 - 0.5*q2Limit/momCM2);
113  if(1.5 > targetA && p == theProton && cosTetMaxNuc < 0.0) cosTetMaxNuc = 0.0;
114  //G4cout << " ctmax= " << cosTetMaxNuc
115  //<< " ctmin= " << cosThetaMin << G4endl; 
116
117  // Cross section in CM system
118  if(cosTetMaxNuc < cosThetaMin) {
119    G4double effmass = mass*m1/(mass + m1);
120    G4double x1 = 1.0 - cosThetaMin;
121    G4double x2 = 1.0 - cosTetMaxNuc;
122    G4double z1 = x1 + screenZ;
123    G4double z2 = x2 + screenZ;
124    G4double d  = 1.0/formfactA;
125    G4double zn1= x1 + d;
126    G4double zn2= x2 + d;
127    nucXSection += coeff*Z*Z*chargeSquare*(1.0 +  effmass*effmass/momCM2)
128      *(1./z1 - 1./z2 + 1./zn1 - 1./zn2 + 
129        2.0*formfactA*std::log(z1*zn2/(z2*zn1)))/momCM2;
130    //G4cout << "XS: x1= " << x1 << " x2= " << x2
131    //<< " cross= " << cross << G4endl;
132    //G4cout << "momCM2= " << momCM2 << " invbeta2= " << invbeta2
133    //       << " coeff= " << coeff << G4endl;
134  }
135  if(nucXSection < 0.0) nucXSection = 0.0;
136  return nucXSection;
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140
141G4double G4CoulombScatteringModel::SelectIsotope(const G4Element* elm)
142{
143  G4double N = elm->GetN(); 
144  G4int ni   = elm->GetNumberOfIsotopes();
145  if(ni > 0) {
146    G4double* ab = elm->GetRelativeAbundanceVector();
147    G4double x = G4UniformRand();
148    G4int idx;
149    for(idx=0; idx<ni; idx++) { 
150      x -= ab[idx];
151      if (x <= 0.0) break;
152    }
153    if(idx >= ni) {
154      G4cout << "G4CoulombScatteringModel::SelectIsotope WARNING: "
155             << "abandance vector for"
156             << elm->GetName() << " is not normalised to unit" << G4endl;
157    } else {
158      N = G4double(elm->GetIsotope(idx)->GetN());
159    }
160  }
161  return N;
162}
163
164//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
165
166void G4CoulombScatteringModel::SampleSecondaries(
167                               std::vector<G4DynamicParticle*>* fvect,
168                               const G4MaterialCutsCouple* couple,
169                               const G4DynamicParticle* dp,
170                               G4double cutEnergy, 
171                               G4double maxEnergy)
172{
173  const G4Material* aMaterial = couple->GetMaterial();
174  const G4ParticleDefinition* p = dp->GetDefinition();
175  G4double kinEnergy = dp->GetKineticEnergy();
176
177  // Select isotope and setup
178  SetupParticle(p);
179  const G4Element* elm = 
180    SelectRandomAtom(aMaterial,p,kinEnergy,cutEnergy,maxEnergy);
181  G4double Z  = elm->GetZ();
182  G4double A  = SelectIsotope(elm);
183  G4int iz    = G4int(Z);
184  G4int ia    = G4int(A + 0.5);
185
186  G4double cross = 
187    ComputeCrossSectionPerAtom(p,kinEnergy,Z,A,cutEnergy,maxEnergy);
188
189  G4double costm = cosTetMaxNuc;
190  G4double formf = formfactA;
191  if(G4UniformRand()*cross < elecXSection) {
192    costm = cosTetMaxElec;
193    formf = 0.0;
194  }
195
196  //  G4cout << "SampleSec: Ekin= " << kinEnergy << " m1= " << m1
197  // << " Z= "<< Z << " A= " <<A<< G4endl;
198
199  if(costm >= cosThetaMin) return; 
200
201  // kinematics in CM system
202  G4double m1   = theParticleTable->GetIonTable()->GetNucleusMass(iz, ia);
203  G4double etot = kinEnergy + mass;
204  G4double ptot = sqrt(mom2);
205  G4double bet  = ptot/(etot + m1);
206  G4double gam  = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
207  G4double pCM  = gam*(ptot - bet*etot);
208  G4double eCM  = gam*(etot - bet*ptot);
209
210  G4double x1 = 1. - cosThetaMin + screenZ;
211  G4double x2 = 1. - costm;
212  G4double x3 = cosThetaMin - costm;
213
214  G4double grej,  z, z1; 
215  do {
216    z  = G4UniformRand()*x3;
217    z1 = (x1*x2 - screenZ*z)/(x1 + z);
218    if(z1 < 0.0) z1 = 0.0;
219    else if(z1 > 2.0) z1 = 2.0;
220    grej = 1.0/(1.0 + formf*z1);
221  } while ( G4UniformRand() > grej*grej ); 
222 
223  G4double cost = 1.0 - z1;
224  G4double sint= sqrt(z1*(2.0 - z1));
225
226  G4double phi = twopi * G4UniformRand();
227
228  // projectile after scattering
229  G4double pzCM = pCM*cost;
230  G4ThreeVector v1(pCM*cos(phi)*sint,pCM*sin(phi)*sint,gam*(pzCM + bet*eCM));
231  G4ThreeVector dir = dp->GetMomentumDirection(); 
232  G4ThreeVector newDirection = v1.unit();
233  newDirection.rotateUz(dir);   
234  fParticleChange->ProposeMomentumDirection(newDirection);   
235  G4double elab = gam*(eCM + bet*pzCM);
236  G4double ekin = elab - mass;
237  if(ekin < 0.0) ekin = 0.0;
238  G4double plab = sqrt(ekin*(ekin + 2.0*mass));
239  fParticleChange->SetProposedKineticEnergy(ekin);
240
241  // recoil
242  G4double erec = kinEnergy - ekin;
243  if(erec > Z*aMaterial->GetIonisation()->GetMeanExcitationEnergy()) {
244    G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
245    G4ThreeVector p2 = (ptot*dir - plab*newDirection).unit();
246    G4DynamicParticle* newdp  = new G4DynamicParticle(ion, p2, erec);
247    fvect->push_back(newdp);
248  } else if(erec > 0.0) {
249    fParticleChange->ProposeLocalEnergyDeposit(erec);
250    fParticleChange->ProposeNonIonizingEnergyDeposit(erec);
251  }
252}
253
254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255
256
Note: See TracBrowser for help on using the repository browser.