source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlungContinuous.cc

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

geant4 tag 9.4
File size: 10.2 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: G4PenelopeBremsstrahlungContinuous.cc,v 1.13 2010/11/25 09:43:47 pandola Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29// --------------------------------------------------------------
30//
31// File name:     G4PenelopeBremsstrahlungContinuous
32//
33// Author:        Luciano Pandola
34//
35// Creation date: February 2003
36// History:
37// -----------
38// 20 Feb 2003  L. Pandola       1st implementation
39// 17 Mar 2003  L. Pandola       Added the correction for positrons
40// 19 Mar 2003  L. Pandola       Bugs fixed
41// 17 Mar 2004  L. Pandola       Removed unnecessary calls to std::pow(a,b)
42// 09 Dec 2008  L. Pandola       Update ReadFile() in a way to make explicit use of the units of
43//                               measurement. Also some cosmetics to improve readibility.
44//----------------------------------------------------------------
45
46#include "G4PenelopeBremsstrahlungContinuous.hh"
47#include "G4PenelopeInterpolator.hh"
48#include "G4Electron.hh"
49#include "G4Positron.hh"
50#include "Randomize.hh"
51#include "G4ParticleDefinition.hh"
52#include "globals.hh"
53#include <sstream>
54
55G4PenelopeBremsstrahlungContinuous::G4PenelopeBremsstrahlungContinuous (G4int Zed,G4double cut,
56                                                                        G4double e1,
57                                                                        G4double e2,
58                                                                        const G4String name)  : 
59  Zmat(Zed),tCut(cut),MinE(e1),MaxE(e2),partName(name)
60{
61  //Construct extended energy table     
62  //200 bins between MinE and MaxE (logarithmic)
63  G4double EL=0.99999*MinE;
64  G4double EU=1.00001*MaxE;
65  DLFC=std::log(EU/EL)/((G4double) (NumberofExtendedEGrid-1));
66  ExtendedLogEnergy[0]=std::log(EL);
67  for (size_t i=1;i<NumberofExtendedEGrid;i++){
68    ExtendedLogEnergy[i]=ExtendedLogEnergy[i-1]+DLFC;
69  }
70  DLFC=1.0/DLFC;
71
72  LoadFromFile();
73  PrepareInterpolationTable();
74}
75
76//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78G4PenelopeBremsstrahlungContinuous::~G4PenelopeBremsstrahlungContinuous()
79{;}
80
81//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82
83void G4PenelopeBremsstrahlungContinuous::LoadFromFile()
84{
85  //Read information from DataBase File
86 char* path = getenv("G4LEDATA");
87 if (!path)
88   {
89     G4String excep = "G4PenelopeBremsstrahlungContinuous - G4LEDATA environment variable not set!";
90     G4Exception(excep);
91     return;
92   }
93 G4String pathString(path);
94 G4String filename = "br-pen-cont-";
95 std::ostringstream ost; 
96 ost << filename << Zmat << ".dat";
97 G4String name(ost.str());
98 G4String dirFile = pathString + "/penelope/" + name;
99 std::ifstream file(dirFile);
100 if (!file.is_open())
101     {
102      G4String excep = "G4PenelopeBremsstrahlungContinuous - data file " + name + " not found!";
103      G4Exception(excep);
104     }
105 G4double a1 = -1.;
106 for (size_t i=0;i<NumberofEPoints;i++)
107   {
108     file >> a1;
109     //1) reads energy in MeV
110     Energies[i]=a1*MeV;
111     //2) read 32 cross sections in cm2
112     for (size_t j=0;j<NumberofKPoints;j++){
113       file >> a1;
114       ReducedCS[i][j]=a1*cm2; 
115     }
116     //3) read the total cross section, in cm2
117     file >> a1;
118     TotalCS[i]=a1*cm2; 
119     // Check closing item
120     file >> a1;
121     if (a1 != ((G4double) -1))
122       {
123         G4String excep = "G4PenelopeBremsstrahlungContinuous - Check the bremms data file " 
124           + name;
125         G4Exception(excep);
126       }
127   }
128 file.close();
129}
130
131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132
133void G4PenelopeBremsstrahlungContinuous::PrepareInterpolationTable()
134{
135  //The energy-loss spectrum is re-normalized to reproduce the
136  //total scaled cross-section of Berger and Seltzer
137
138  G4double pK[NumberofKPoints] = {1.0e-12,0.05,0.075,0.1,0.125,0.15,0.2,0.25,
139                                  0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,
140                                  0.75,0.8,0.85,0.9,0.925,0.95,0.97,0.99,
141                                  0.995,0.999,0.9995,0.9999,0.99995,0.99999,1.0};
142  G4double pY[NumberofKPoints];
143
144  size_t i=0,j=0;
145  for (i=0;i<NumberofEPoints;i++){
146    for (j=0;j<NumberofKPoints;j++){
147      pY[j]= ReducedCS[i][j];
148    }
149    G4PenelopeInterpolator* interpolator = new G4PenelopeInterpolator(pK,pY,NumberofKPoints);
150    G4double Rsum = interpolator->CalculateMomentum(1.0,0);
151   
152    delete interpolator;
153    G4double Fact = (Energies[i]+electron_mass_c2)*(1.0/fine_structure_const)*millibarn/
154      (classic_electr_radius*classic_electr_radius*(Energies[i]+2.0*electron_mass_c2));
155    G4double Normalization = TotalCS[i]/(Rsum*Fact);
156    G4double TST = std::abs(Normalization-1.0);
157    if (TST > 1.0) {
158      G4String excep = "G4PenelopeBremsstrahlungContinuous - Check the bremms data file";
159      G4Exception(excep);
160    }
161    for (j=0;j<NumberofKPoints;j++){
162      ReducedCS[i][j] = ReducedCS[i][j]*Normalization;
163    }
164  }
165
166  //Compute the scaled energy loss distribution and sampling parameters
167  // for the energies in the simulation grid
168 
169  // Interpolation in E
170  G4double pX[NumberofEPoints];
171  G4double pYY[NumberofEPoints];
172  for (i=0;i<NumberofEPoints;i++){
173    pX[i] = std::log(Energies[i]);
174  }
175 
176  for (j=0;j<NumberofKPoints;j++){
177    for (i=0;i<NumberofEPoints;i++){
178      pYY[i] = std::log(ReducedCS[i][j]);
179    }
180    G4PenelopeInterpolator* interpolator2 = new G4PenelopeInterpolator(pX,pYY,NumberofEPoints);
181    for (i=0;i<NumberofExtendedEGrid;i++){
182      G4double ELL = ExtendedLogEnergy[i];
183      if (ELL >= pX[0]) {
184        p0[i][j] = std::exp(interpolator2->CubicSplineInterpolation(ELL));
185      }
186      else
187        {
188          G4double F1=interpolator2->CubicSplineInterpolation(pX[0]);
189          G4double FP1 = interpolator2->FirstDerivative(pX[0]);
190          p0[i][j] = std::exp(F1+FP1*(ELL-pX[0]));
191        }
192    }
193    delete interpolator2;
194  }
195 
196  //All this stuff might be useful later on
197  /*
198  G4double PDF[NumberofKPoints];
199 
200  for (i=0;i<NumberofExtendedEGrid;i++){
201    for (j=0;j<NumberofKPoints;j++){
202      PDF[j]=p0[i][j];
203    }
204    G4double Xc=0;
205    if (i<(NumberofExtendedEGrid-1)){
206      Xc=tCut/std::exp(ExtendedLogEnergy[i+1]);
207    }
208    else
209      {
210        Xc=tCut/std::exp(ExtendedLogEnergy[NumberofExtendedEGrid-1]);
211      }
212   
213    G4PenelopeInterpolator* interpolator3 = new G4PenelopeInterpolator(pK,PDF,NumberofKPoints);
214    Pbcut[i]=interpolator3->CalculateMomentum(Xc,-1);
215    delete interpolator3;
216  }
217  */
218 
219}
220
221//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222                       
223G4double G4PenelopeBremsstrahlungContinuous::CalculateStopping(G4double energy)
224  //Stopping power expressed in MeV/mm*2
225{
226  G4double Xel=std::max(std::log(energy),ExtendedLogEnergy[0]);
227  G4double Xe=1.0+(Xel-ExtendedLogEnergy[0])*DLFC;
228  G4int Ke = (G4int) Xe; 
229  G4double Xek = Xe-Ke; 
230 
231  //Global x-section factor
232  G4double Fact=Zmat*Zmat*(energy+electron_mass_c2)*(energy+electron_mass_c2)/
233    (energy*(energy+2.0*electron_mass_c2));
234  Fact *= PositronCorrection(energy);
235
236  //Moments of the scaled bremss x-section
237  G4double wcre = tCut/energy;
238  G4double pY[NumberofKPoints];
239  G4double pK[NumberofKPoints] = {1.0e-12,0.05,0.075,0.1,0.125,0.15,0.2,0.25,
240                                  0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,
241                                  0.75,0.8,0.85,0.9,0.925,0.95,0.97,0.99,
242                                  0.995,0.999,0.9995,0.9999,0.99995,0.99999,1.0};
243
244  for (size_t i=0;i<NumberofKPoints;i++){
245    pY[i] = p0[Ke][i];
246  }
247  G4PenelopeInterpolator* interpolator1 = new G4PenelopeInterpolator(pK,pY,NumberofKPoints);
248  G4double XS1A = interpolator1->CalculateMomentum(wcre,0);
249  for (size_t k=0;k<NumberofKPoints;k++){
250    pY[k] = p0[std::min(Ke+1,(G4int) NumberofExtendedEGrid-1)][k];
251  }
252  G4PenelopeInterpolator* interpolator2 = new G4PenelopeInterpolator (pK,pY,NumberofKPoints);
253  G4double XS1B = interpolator2->CalculateMomentum(wcre,0);
254 
255  G4double XS1 = ((1.0-Xek)*XS1A+Xek*XS1B)*Fact*energy; //weighted mean between the energy bin of the grid
256
257  //This is the 2nd momentum (straggling cross section). It might be useful later on
258  /*
259    G4double XS2A = interpolator1->CalculateMomentum(wcre,1);
260    G4double XS2 = ((1.0-Xek)*XS2A+Xek*XS2B)*Fact*energy*energy;
261    G4double XS2B = interpolator2->CalculateMomentum(wcre,1);
262  */
263
264  delete interpolator1;
265  delete interpolator2;
266
267  return XS1;
268}
269
270//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
271
272G4double G4PenelopeBremsstrahlungContinuous::PositronCorrection(G4double energy)
273{
274  const G4double Coeff[7]={-1.2359e-01,6.1274e-2,-3.1516e-2,7.7446e-3,-1.0595e-3,
275                           7.0568e-5,-1.8080e-6};
276  G4double correct=0;
277  if (partName == "e-") {
278    return 1.0; //no correction for electrons
279  }
280  else if (partName == "e+"){
281    G4double T=std::log(1+((1e6*energy)/(Zmat*Zmat*electron_mass_c2)));
282    for (G4int i=0;i<7;i++){
283      correct += Coeff[i]*std::pow(T,i+1);
284    }
285    correct = 1.0-std::exp(correct);
286    return correct;
287  }
288  else //neither electrons nor positrons...exception
289    {
290      G4String excep = "G4PenelopeBremmstrahlungContinuous: the particle is not e- nor e+!";
291      G4Exception(excep);
292      return 0;
293    }
294}
Note: See TracBrowser for help on using the repository browser.