source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeGammaConversionModel.cc @ 1315

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

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 17.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: G4PenelopeGammaConversionModel.cc,v 1.6 2009/06/11 15:47:08 mantero Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 06 Oct 2008   L Pandola    Migration from process to model
34// 17 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
35//                  - apply internal high-energy limit only in constructor
36//                  - do not apply low-energy limit (default is 0)
37//                  - do not apply production threshold on level of the model
38// 19 May 2009   L Pandola    Explicitely set to zero pointers deleted in
39//                            Initialise(), since they might be checked later on
40//
41
42#include "G4PenelopeGammaConversionModel.hh"
43#include "G4ParticleDefinition.hh"
44#include "G4MaterialCutsCouple.hh"
45#include "G4ProductionCutsTable.hh"
46#include "G4DynamicParticle.hh"
47#include "G4Element.hh"
48#include "G4Gamma.hh"
49#include "G4Electron.hh"
50#include "G4Positron.hh"
51#include "G4CrossSectionHandler.hh"
52
53//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54
55
56G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel(const G4ParticleDefinition*,
57                                             const G4String& nam)
58  :G4VEmModel(nam),fTheScreeningRadii(0),crossSectionHandler(0),isInitialised(false)
59{
60  fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2;
61  fIntrinsicHighEnergyLimit = 100.0*GeV;
62  fSmallEnergy = 1.1*MeV;
63
64  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
65  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
66  //
67  verboseLevel= 0;
68  // Verbosity scale:
69  // 0 = nothing
70  // 1 = warning for energy non-conservation
71  // 2 = details of energy budget
72  // 3 = calculation of cross sections, file openings, sampling of atoms
73  // 4 = entering in methods
74}
75
76//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel()
79{ 
80  if (crossSectionHandler) delete crossSectionHandler;
81  if (fTheScreeningRadii) delete fTheScreeningRadii;
82}
83
84//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85
86void G4PenelopeGammaConversionModel::Initialise(const G4ParticleDefinition*,
87                                                const G4DataVector& )
88{
89  if (verboseLevel > 3)
90    G4cout << "Calling  G4PenelopeGammaConversionModel::Initialise()" << G4endl;
91
92  //Delete the old cross section handler, if necessary
93  if (crossSectionHandler)
94    {
95      crossSectionHandler->Clear();
96      delete crossSectionHandler;
97      crossSectionHandler = 0;
98    }
99 
100  //Re-initialize cross section handler
101  crossSectionHandler = new G4CrossSectionHandler();
102  crossSectionHandler->Initialise(0,fIntrinsicLowEnergyLimit,HighEnergyLimit(),400);
103  crossSectionHandler->Clear();
104  G4String crossSectionFile = "penelope/pp-cs-pen-";
105  crossSectionHandler->LoadData(crossSectionFile);
106  //This is used to retrieve cross section values later on
107  crossSectionHandler->BuildMeanFreePathForMaterials();
108
109  if (verboseLevel > 2) 
110    G4cout << "Loaded cross section files for PenelopeGammaConversion" << G4endl;
111
112  if (verboseLevel > 0) { 
113    G4cout << "Penelope Gamma Conversion model is initialized " << G4endl
114           << "Energy range: "
115           << LowEnergyLimit() / MeV << " MeV - "
116           << HighEnergyLimit() / GeV << " GeV"
117           << G4endl;
118  }
119
120  if(isInitialised) return;
121  fParticleChange = GetParticleChangeForGamma();
122  isInitialised = true;
123}
124
125//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
126
127G4double G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom(
128                                       const G4ParticleDefinition*,
129                                             G4double energy,
130                                             G4double Z, G4double,
131                                             G4double, G4double)
132{
133  //
134  // Penelope model.
135  // Cross section (including triplet production) read from database and managed
136  // through the G4CrossSectionHandler utility. Cross section data are from
137  // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598
138  //
139 
140  if (verboseLevel > 3)
141    G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl;
142
143  G4int iZ = (G4int) Z;
144  //  if (!crossSectionHandler) //VI: should not be checked in run time
145  //  {
146  //    G4cout << "G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom" << G4endl;
147  //    G4cout << "The cross section handler is not correctly initialized" << G4endl;
148  //    G4Exception();
149  //  }
150  G4double cs = crossSectionHandler->FindValue(iZ,energy);
151
152  if (verboseLevel > 2)
153    G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z << 
154      " = " << cs/barn << " barn" << G4endl;
155  return cs;
156}
157
158//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
159
160void 
161G4PenelopeGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
162                                                  const G4MaterialCutsCouple* couple,
163                                                  const G4DynamicParticle* aDynamicGamma,
164                                                  G4double,
165                                                  G4double)
166{
167  //
168  // Penelope model.
169  // Final state is sampled according to the Bethe-Heitler model with Coulomb
170  // corrections, according to the semi-empirical model of
171  //  J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531.
172  //
173  // The model uses the high energy Coulomb correction from
174  //  H. Davies et al., Phys. Rev. 93 (1954) 788
175  // and atomic screening radii tabulated from
176  //  J.H. Hubbel et al., J. Phys. Chem. Ref. Data 9 (1980) 1023
177  // for Z= 1 to 92. This managed in this model by the method
178  // GetScreeningRadius().
179  //
180  if (verboseLevel > 3)
181    G4cout << "Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" << G4endl;
182
183  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
184
185  // Always kill primary
186  fParticleChange->ProposeTrackStatus(fStopAndKill);
187  fParticleChange->SetProposedKineticEnergy(0.);
188
189  if (photonEnergy <= fIntrinsicLowEnergyLimit)
190    {
191      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
192      return ;
193    }
194
195  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
196
197  G4double eps ;
198  G4double eki = electron_mass_c2 / photonEnergy ;
199
200  // Do it fast if photon energy < 1.1 MeV
201  if (photonEnergy < fSmallEnergy )
202    {
203      eps = eki + (1-2*eki) * G4UniformRand();
204    }
205  else
206    {
207      // Select randomly one element in the current material
208      if (verboseLevel > 2)
209        G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
210      //use crossSectionHandler instead of G4EmElementSelector because in this case
211      //the dimension of the table is equal to the dimension of the database
212      //(less interpolation errors)
213      G4int Z_int = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
214      if (verboseLevel > 2)
215        G4cout << "Selected Z = " << Z_int << G4endl;
216     
217      //Low energy and Coulomb corrections
218      G4double Z=(G4double) Z_int;
219      G4double ZAlpha = Z*fine_structure_const;
220      G4double ScreenRadius = GetScreeningRadius(Z);
221      G4double funct1=0,g0=0;
222      G4double g1min=0,g2min=0;
223      funct1 = 4.0*std::log(ScreenRadius);
224      g0 = funct1-4*CoulombCorrection(ZAlpha)+LowEnergyCorrection(ZAlpha,eki); 
225      G4double bmin = 2*eki*ScreenRadius;
226      std::vector<G4double> ScreenFunctionValues = ScreenFunction(bmin);
227      if (ScreenFunctionValues.size() != 2)
228        {
229          G4cout << "G4PenelopeGammaConversionModel::SampleSecondaries" << G4endl;
230          G4cout << "ScreenFunction did not return 2 values! Something wrong! " << G4endl;
231          G4Exception();
232        }
233      g1min=g0+ScreenFunctionValues[0];
234      g2min=g0+ScreenFunctionValues[1];
235      G4double xr,a1,p1;
236      xr=0.5-eki;
237      a1=(2.0/3.0)*g1min*xr*xr;
238      p1=a1/(a1+g2min);
239
240      //Random sampling of eps
241      G4double rand1,rand2,rand3,b;
242      G4double g1;
243   
244      do{
245        rand1 = G4UniformRand();
246        if (rand1 < p1) {
247          rand2 = 2.0*G4UniformRand()-1.0;
248          if (rand2 < 0) {
249            eps = 0.5 - xr*std::pow(std::abs(rand2),(1./3.));
250          }
251          else
252            {
253              eps = 0.5 + xr*std::pow(rand2,(1./3.));
254            }
255          b = (eki*ScreenRadius)/(2*eps*(1.0-eps));
256          std::vector<G4double> ScreenFunctionSampling = ScreenFunction(b);
257          g1 = g0+ScreenFunctionSampling[0];
258          if (g1 < 0) g1=0;
259          rand3 = G4UniformRand()*g1min;
260        }
261        else
262          {
263            eps = eki+2.0*xr*G4UniformRand();
264            b = (eki*ScreenRadius)/(2*eps*(1.0-eps));
265            std::vector<G4double> ScreenFunctionSampling = ScreenFunction(b);
266            g1 = g0+ScreenFunctionSampling[1];
267            if (g1 < 0) g1=0; 
268            rand3 = G4UniformRand()*g2min;
269          }       
270      } while (rand3>g1);
271    } //End of eps sampling
272
273  G4double electronTotEnergy = eps*photonEnergy;
274  G4double positronTotEnergy = (1.0-eps)*photonEnergy;
275 
276  // Scattered electron (positron) angles. ( Z - axis along the parent photon)
277
278  //electron kinematics
279  G4double costheta_el,costheta_po;
280  G4double phi_el,phi_po;
281  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ; 
282  costheta_el = G4UniformRand()*2.0-1.0;
283  G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
284  costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
285  phi_el  = twopi * G4UniformRand() ;
286  G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
287  G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
288  G4double dirZ_el = costheta_el;
289
290  //positron kinematics
291  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
292  costheta_po = G4UniformRand()*2.0-1.0;
293  kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
294  costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
295  phi_po  = twopi * G4UniformRand() ;
296  G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
297  G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
298  G4double dirZ_po = costheta_po;
299
300  // Kinematics of the created pair:
301  // the electron and positron are assumed to have a symetric angular
302  // distribution with respect to the Z axis along the parent photon
303  G4double localEnergyDeposit = 0. ;
304
305  //Generate explicitely the electron in the pair, only if it is > threshold
306  //VI: applying cut here provides inconsistency
307
308  if (electronKineEnergy > 0.0)
309    {
310      G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
311      electronDirection.rotateUz(photonDirection);
312      G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
313                                                           electronDirection,
314                                                           electronKineEnergy);
315      fvect->push_back(electron);
316    }
317  else
318    {
319      localEnergyDeposit += electronKineEnergy;
320      electronKineEnergy = 0;
321    }
322
323  //Generate the positron. Real particle in any case, because it will annihilate. If below
324  //threshold, produce it at rest
325  // VI: here there was a bug - positron and electron cuts are different
326  if (positronKineEnergy < 0.0)
327    {
328      localEnergyDeposit += positronKineEnergy;
329      positronKineEnergy = 0; //produce it at rest
330    }
331  G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po);
332  positronDirection.rotateUz(photonDirection);
333  G4DynamicParticle* positron = new G4DynamicParticle(G4Positron::Positron(),
334                                                      positronDirection, positronKineEnergy);
335  fvect->push_back(positron);
336
337  //Add rest of energy to the local energy deposit
338  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
339 
340  if (verboseLevel > 1)
341    {
342      G4cout << "-----------------------------------------------------------" << G4endl;
343      G4cout << "Energy balance from G4PenelopeGammaConversion" << G4endl;
344      G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
345      G4cout << "-----------------------------------------------------------" << G4endl;
346      if (electronKineEnergy)
347        G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV" 
348               << G4endl;
349      if (positronKineEnergy)
350        G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl;
351      G4cout << "Rest masses of e+/- " << 2.0*electron_mass_c2/keV << " keV" << G4endl;
352      if (localEnergyDeposit)
353        G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
354      G4cout << "Total final state: " << (electronKineEnergy+positronKineEnergy+
355                                          localEnergyDeposit+2.0*electron_mass_c2)/keV <<
356        " keV" << G4endl;
357      G4cout << "-----------------------------------------------------------" << G4endl;
358    }
359 if (verboseLevel > 0)
360    {
361      G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
362                                      localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
363      if (energyDiff > 0.05*keV)
364        G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: " 
365               << (electronKineEnergy+positronKineEnergy+
366                   localEnergyDeposit+2.0*electron_mass_c2)/keV
367               << " keV (final) vs. " << photonEnergy/keV << " keV (initial)" << G4endl;
368    } 
369}
370
371//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
372
373std::vector<G4double> G4PenelopeGammaConversionModel::ScreenFunction(G4double b)
374{
375  std::vector<G4double> result;
376  result.clear();
377  G4double bsquare=b*b;
378  G4double a0,f1,f2;
379  f1=2.0-2*std::log(1+bsquare);
380  f2=f1-(2.0/3.0);
381  if (b < 1.0e-10)
382    {
383      f1=f1-twopi*b;
384    }
385  else
386    {
387      a0 = 4*b*std::atan(1.0/b);
388      f1 = f1 - a0;
389      f2 = f2+2*bsquare*(4.0-a0-3*std::log((1+bsquare)/bsquare));
390    }
391  result.push_back(0.5*(3*f1-f2));
392  result.push_back(0.25*(3*f1+f2));
393  return result;
394}
395
396//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
397
398G4double G4PenelopeGammaConversionModel::GetScreeningRadius(G4double Z)
399{
400  G4double result = 0;
401  G4bool foundElement = false;
402  G4int iZ = (G4int) Z;
403  if (!fTheScreeningRadii)
404    fTheScreeningRadii = new std::map<G4int,G4double>;
405 
406  if (fTheScreeningRadii->count(iZ))
407    {
408      //The element is already loaded: just return it
409      result = fTheScreeningRadii->find(iZ)->second;
410      return result;
411    }
412  else //retrieve all from file
413    {
414      char* path = getenv("G4LEDATA");
415      if (!path)
416        {
417          G4String excep = "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!";
418          G4Exception(excep);
419        }
420      G4String pathString(path);
421      G4String pathFile = pathString + "/penelope/pp-pen.dat";
422      std::ifstream file(pathFile);
423     
424      if (!(file.is_open()))
425        {
426          G4String excep = "G4PenelopeGammaConversionModel - data file " + pathFile + "not found!";
427          G4Exception(excep);
428        }
429      G4int k;
430      G4double a1,a2;
431      while(!file.eof()) {
432        file >> k >> a1 >> a2;
433        fTheScreeningRadii->insert(std::make_pair(k,a1));
434        if ((G4double) k == Z)
435          {
436            result = a1;
437            foundElement = true;
438          }
439      }
440      file.close();
441      if (verboseLevel > 2)
442        G4cout << "Read file pp-pen.dat" << G4endl;
443      if (foundElement)
444        return result;
445      else
446        {
447          G4String excep = "G4PenelopeGammaConversionModel - Screening Radius for not found in the data file";
448          G4Exception(excep);
449          return 0;
450        }
451    }
452}
453
454//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
455
456G4double G4PenelopeGammaConversionModel::CoulombCorrection(G4double a)
457{
458  G4double fc=0;
459  G4double b[7] = {0.202059,-0.03693,0.00835,-0.00201,0.00049,-0.00012,0.00003};
460  G4double aSquared = a*a;
461  G4double aFourth = aSquared*aSquared;
462  G4double aEighth = aFourth*aFourth;
463
464  fc = ((1.0/(1.0+a*a))+b[0]+b[1]*aSquared+b[2]*aFourth+b[3]*(aSquared*aFourth)+
465        b[4]*aEighth+b[5]*(aEighth*aSquared)+b[6]*(aEighth*aFourth));
466  fc=aSquared*fc;
467  return fc;
468}
469
470//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
471
472G4double G4PenelopeGammaConversionModel::LowEnergyCorrection(G4double a,G4double eki)
473{
474  G4double f0=0,t=0;
475  G4double b[12] = {-1.744,-12.10,11.18,8.523,73.26,-41.41,-13.52,-121.1,94.41,8.946,62.05,-63.41};
476  t=std::sqrt(2.0*eki);
477  G4double tSq = t*t;
478  f0=(b[0]+b[1]*a+b[2]*a*a)*t+(b[3]+b[4]*a+b[5]*a*a)*(tSq)+(b[6]+b[7]*a+b[8]*a*a)*(tSq*t)+
479    (b[9]+b[10]*a+b[11]*a*a)*(tSq*tSq);
480  return f0;
481
482}
Note: See TracBrowser for help on using the repository browser.