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

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

import all except CVS
File size: 11.5 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.40 2008/01/07 08:32:01 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
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// 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.08.06 V.Ivanchenko add inline function ScreeningParameter
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 "G4eCoulombScatteringModel.hh"
55#include "Randomize.hh"
56#include "G4DataVector.hh"
57#include "G4ElementTable.hh"
58#include "G4PhysicsLogVector.hh"
59#include "G4ParticleChangeForGamma.hh"
60#include "G4Electron.hh"
61#include "G4Positron.hh"
62#include "G4Proton.hh"
63
64//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65
66using namespace std;
67
68G4eCoulombScatteringModel::G4eCoulombScatteringModel(
69  G4double thetaMin, G4double thetaMax, G4bool build, 
70  G4double tlim, const G4String& nam)
71  : G4VEmModel(nam),
72    cosThetaMin(cos(thetaMin)),
73    cosThetaMax(cos(thetaMax)),
74    q2Limit(tlim),
75    theCrossSectionTable(0),
76    lowKEnergy(keV),
77    highKEnergy(TeV),
78    alpha2(fine_structure_const*fine_structure_const),
79    faclim(100.0),
80    nbins(12),
81    nmax(100),
82    buildTable(build),
83    isInitialised(false)
84{
85  fNistManager = G4NistManager::Instance();
86  theElectron = G4Electron::Electron();
87  thePositron = G4Positron::Positron();
88  theProton   = G4Proton::Proton();
89  a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885);
90  G4double p0 = electron_mass_c2*classic_electr_radius;
91  coeff  = twopi*p0*p0;
92  constn = 6.937e-6/(MeV*MeV);
93  tkin = targetZ = targetA = mom2 = DBL_MIN;
94  elecXSection = nucXSection = 0.0;
95  ecut = DBL_MAX;
96  particle = 0;
97  for(size_t j=0; j<100; j++) {index[j] = -1;} 
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101
102G4eCoulombScatteringModel::~G4eCoulombScatteringModel()
103{
104  if(theCrossSectionTable) {
105    theCrossSectionTable->clearAndDestroy();
106    delete theCrossSectionTable;
107  }
108}
109
110//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111
112void G4eCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
113                                           const G4DataVector&)
114{
115  //  G4cout << "!!! G4eCoulombScatteringModel::Initialise for "
116  // << p->GetParticleName() << "  cos(TetMin)= " << cosThetaMin
117  // << "  cos(TetMax)= " << cosThetaMax <<G4endl;
118  if(!isInitialised) {
119    isInitialised = true;
120
121    if(pParticleChange)
122      fParticleChange = 
123        reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
124    else
125      fParticleChange = new G4ParticleChangeForGamma();
126  } else {
127    return;
128  }
129
130  if(p->GetParticleType() == "nucleus") buildTable = false;
131  if(!buildTable) return;
132
133  // Compute log cross section table per atom
134  if(!theCrossSectionTable) theCrossSectionTable = new G4PhysicsTable();
135 
136  nbins = 2*G4int(log10(highKEnergy/lowKEnergy));
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140
141G4double G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(
142                const G4ParticleDefinition* p,
143                G4double kinEnergy,
144                G4double Z, G4double A,
145                G4double cutEnergy, G4double)
146{
147  if(p == particle && kinEnergy == tkin && Z == targetZ &&
148     A == targetA && cutEnergy == ecut) return nucXSection;
149
150  //G4cout << "### G4eCoulombScatteringModel::ComputeCrossSectionPerAtom  for "
151  //     << p->GetParticleName() << " Z= " << Z << " A= " << A
152  //     << " e= " << kinEnergy << G4endl;
153
154  nucXSection = ComputeElectronXSectionPerAtom(p,kinEnergy,Z,A,cutEnergy);
155
156  // nuclear cross section
157  if(theCrossSectionTable) {
158    G4bool b;
159    G4int iz  = G4int(Z);
160    G4int idx = index[iz];
161
162    // compute table for given Z
163    if(-1 == idx) {
164      idx = theCrossSectionTable->size();
165      index[iz] = idx;
166      G4PhysicsLogVector* ptrVector
167        = new G4PhysicsLogVector(lowKEnergy, highKEnergy, nbins);
168      //  G4cout << "New vector Z= " << iz << " A= " << A << " idx= " << idx << G4endl;
169      G4double e, value;
170      for(G4int i=0; i<=nbins; i++) {
171        e     = ptrVector->GetLowEdgeEnergy( i ) ;
172        value = CalculateCrossSectionPerAtom(p, e, Z, A); 
173        ptrVector->PutValue( i, log(value) );
174      }
175      theCrossSectionTable->push_back(ptrVector);
176    }
177
178      // take value from the table
179    nucXSection += 
180      std::exp((((*theCrossSectionTable)[idx]))->GetValue(kinEnergy, b));
181
182    // compute value from scratch
183  } else nucXSection += CalculateCrossSectionPerAtom(p, kinEnergy, Z, A);
184 
185  //  G4cout << " cross(bn)= " << nucXSection/barn << G4endl;
186 
187  if(nucXSection < 0.0) nucXSection = 0.0;
188  return nucXSection;
189}
190
191//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
192
193G4double G4eCoulombScatteringModel::ComputeElectronXSectionPerAtom(
194                const G4ParticleDefinition* p,
195                G4double kinEnergy,
196                G4double Z,
197                G4double A,
198                G4double cutEnergy)
199{
200  if(p == particle && kinEnergy == tkin && Z == targetZ &&
201     cutEnergy == ecut) return elecXSection;
202  ecut = cutEnergy;
203  elecXSection = 0.0;
204  SetupParticle(p);
205  G4double ekin = std::max(keV, kinEnergy);
206  //G4double ekin = kinEnergy;
207  SetupTarget(Z, A, ekin);
208
209  G4double tmax = tkin;
210  if(p == theElectron) tmax *= 0.5;
211  else if(p != thePositron) {
212    G4double ratio = electron_mass_c2/mass;
213    G4double tau = tkin/mass;
214    tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
215      (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio); 
216  }
217
218  cosTetMaxElec = cosTetMaxNuc;
219  G4double t = std::min(cutEnergy, tmax);
220  G4double mom21 = t*(t + 2.0*electron_mass_c2);
221  G4double t1 = tkin - t;
222  if(t1 > 0.0) {
223    G4double mom22 = t1*(t1 + 2.0*mass);
224    G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
225    if(ctm > cosTetMaxElec && ctm <= 1.0) cosTetMaxElec = ctm;
226  }
227
228  if(cosTetMaxElec < cosThetaMin) {
229    G4double x1 = 1.0 - cosThetaMin  + screenZ;
230    G4double x2 = 1.0 - cosTetMaxElec + screenZ;
231    elecXSection = coeff*Z*chargeSquare*invbeta2*
232      (cosThetaMin - cosTetMaxElec)/(x1*x2*mom2);
233  }
234  //  G4cout << "cut= " << ecut << " e= " << tkin
235  // << " croosE(barn)= " << elecXSection/barn
236  // << " cosEl= " << cosTetMaxElec << " costmin= " << cosThetaMin << G4endl;
237  return elecXSection;
238}
239
240//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
241
242G4double G4eCoulombScatteringModel::CalculateCrossSectionPerAtom(
243                             const G4ParticleDefinition* p,
244                             G4double kinEnergy, 
245                             G4double Z, G4double A)
246{
247  G4double cross = 0.0;
248  SetupParticle(p);
249  G4double ekin = std::max(keV, kinEnergy);
250  //G4double ekin = kinEnergy;
251  SetupTarget(Z, A, ekin);
252
253  if(cosTetMaxNuc < cosThetaMin) {
254    G4double x1 = 1.0 - cosThetaMin;
255    G4double x2 = 1.0 - cosTetMaxNuc;
256    G4double x3 = cosThetaMin - cosTetMaxNuc;
257    G4double z1 = x1 + screenZ;
258    G4double z2 = x2 + screenZ;
259    G4double d  = 1.0/formfactA - screenZ;
260    G4double d1 = 1.0 - formfactA*screenZ;
261    G4double zn1= x1 + d;
262    G4double zn2= x2 + d;
263    cross = coeff*Z*Z*chargeSquare*invbeta2
264      *(x3/(z1*z2) + x3/(zn1*zn2) + 
265        2.0*std::log(z1*zn2/(z2*zn1))/d) / (mom2*d1*d1);
266  }
267 
268  //  G4cout << "CalculateCrossSectionPerAtom: e(MeV)= " << tkin
269  // << " cross(b)= " << cross/barn << " ctmin= " << cosThetaMin
270  // << " ctmax= " << cosTetMaxNuc << G4endl;
271 
272  return cross;
273}
274
275//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
276
277void G4eCoulombScatteringModel::SampleSecondaries(
278                std::vector<G4DynamicParticle*>*,
279                const G4MaterialCutsCouple* couple,
280                const G4DynamicParticle* dp,
281                G4double cutEnergy,
282                G4double maxEnergy)
283{
284  const G4Material* aMaterial = couple->GetMaterial();
285  const G4ParticleDefinition* p = dp->GetDefinition();
286  G4double kinEnergy = dp->GetKineticEnergy();
287
288  // Select atom and setup
289  SetupParticle(p);
290  const G4Element* elm = 
291    SelectRandomAtom(aMaterial,p,kinEnergy,cutEnergy,maxEnergy);
292  G4double Z  = elm->GetZ();
293  G4double A  = elm->GetN();
294
295  G4double cross = 
296    ComputeCrossSectionPerAtom(p,kinEnergy,Z,A,cutEnergy,maxEnergy);
297
298  G4double costm = cosTetMaxNuc;
299  G4double formf = formfactA;
300  if(G4UniformRand()*cross < elecXSection) {
301    costm = cosTetMaxElec;
302    formf = 0.0;
303  }
304  /*
305  G4cout << "G4eCoul...SampleSecondaries: e(MeV)= " << tkin
306         << " ctmin= " << cosThetaMin
307         << " ctmaxN= " << cosTetMaxNuc
308         << " ctmax= " << costm
309         << " Z= " << Z << " A= " << A
310         << " cross= " << cross/barn << " crossE= " << elecXSection/barn
311         << G4endl;
312  */
313  if(costm >= cosThetaMin) return; 
314
315  G4double x1 = 1. - cosThetaMin + screenZ;
316  G4double x2 = 1. - costm;
317  G4double x3 = cosThetaMin - costm;
318  G4double grej,  z, z1; 
319  do {
320    z  = G4UniformRand()*x3;
321    z1 = (x1*x2 - screenZ*z)/(x1 + z);
322    if(z1 < 0.0) z1 = 0.0;
323    else if(z1 > 2.0) z1 = 2.0;
324    grej = 1.0/(1.0 + formf*z1);
325  } while ( G4UniformRand() > grej*grej ); 
326 
327  G4double cost = 1.0 - z1;
328  G4double sint= sqrt(z1*(2.0 - z1));
329  /*
330  if(sint > 0.1)
331    G4cout<<"## SampleSecondaries: e(MeV)= " << kinEnergy
332          << " sint= " << sint << "  Z= " << Z << "  screenZ= " << screenZ
333          << " cn= " << formf
334          << G4endl;
335  */
336  G4double phi  = twopi * G4UniformRand();
337
338  G4ThreeVector direction = dp->GetMomentumDirection(); 
339  G4ThreeVector newDirection(cos(phi)*sint,sin(phi)*sint,cost);
340  newDirection.rotateUz(direction);   
341
342  fParticleChange->ProposeMomentumDirection(newDirection);   
343 
344  return;
345}
346
347//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
348
349
Note: See TracBrowser for help on using the repository browser.