source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlungContinuous.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: 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.12 2009/06/10 13:32:36 mantero Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
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   }
92 G4String pathString(path);
93 G4String filename = "br-pen-cont-";
94 std::ostringstream ost; 
95 ost << filename << Zmat << ".dat";
96 G4String name(ost.str());
97 G4String dirFile = pathString + "/penelope/" + name;
98 std::ifstream file(dirFile);
99 if (!file.is_open())
100     {
101      G4String excep = "G4PenelopeBremsstrahlungContinuous - data file " + name + " not found!";
102      G4Exception(excep);
103     }
104 G4double a1 = -1.;
105 for (size_t i=0;i<NumberofEPoints;i++)
106   {
107     file >> a1;
108     //1) reads energy in MeV
109     Energies[i]=a1*MeV;
110     //2) read 32 cross sections in cm2
111     for (size_t j=0;j<NumberofKPoints;j++){
112       file >> a1;
113       ReducedCS[i][j]=a1*cm2; 
114     }
115     //3) read the total cross section, in cm2
116     file >> a1;
117     TotalCS[i]=a1*cm2; 
118     // Check closing item
119     file >> a1;
120     if (a1 != ((G4double) -1))
121       {
122         G4String excep = "G4PenelopeBremsstrahlungContinuous - Check the bremms data file " 
123           + name;
124         G4Exception(excep);
125       }
126   }
127 file.close();
128}
129
130//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131
132void G4PenelopeBremsstrahlungContinuous::PrepareInterpolationTable()
133{
134  //The energy-loss spectrum is re-normalized to reproduce the
135  //total scaled cross-section of Berger and Seltzer
136
137  G4double pK[NumberofKPoints] = {1.0e-12,0.05,0.075,0.1,0.125,0.15,0.2,0.25,
138                                  0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,
139                                  0.75,0.8,0.85,0.9,0.925,0.95,0.97,0.99,
140                                  0.995,0.999,0.9995,0.9999,0.99995,0.99999,1.0};
141  G4double pY[NumberofKPoints];
142
143  size_t i=0,j=0;
144  for (i=0;i<NumberofEPoints;i++){
145    for (j=0;j<NumberofKPoints;j++){
146      pY[j]= ReducedCS[i][j];
147    }
148    G4PenelopeInterpolator* interpolator = new G4PenelopeInterpolator(pK,pY,NumberofKPoints);
149    G4double Rsum = interpolator->CalculateMomentum(1.0,0);
150   
151    delete interpolator;
152    G4double Fact = (Energies[i]+electron_mass_c2)*(1.0/fine_structure_const)*millibarn/
153      (classic_electr_radius*classic_electr_radius*(Energies[i]+2.0*electron_mass_c2));
154    G4double Normalization = TotalCS[i]/(Rsum*Fact);
155    G4double TST = std::abs(Normalization-1.0);
156    if (TST > 1.0) {
157      G4String excep = "G4PenelopeBremsstrahlungContinuous - Check the bremms data file";
158      G4Exception(excep);
159    }
160    for (j=0;j<NumberofKPoints;j++){
161      ReducedCS[i][j] = ReducedCS[i][j]*Normalization;
162    }
163  }
164
165  //Compute the scaled energy loss distribution and sampling parameters
166  // for the energies in the simulation grid
167 
168  // Interpolation in E
169  G4double pX[NumberofEPoints];
170  G4double pYY[NumberofEPoints];
171  for (i=0;i<NumberofEPoints;i++){
172    pX[i] = std::log(Energies[i]);
173  }
174 
175  for (j=0;j<NumberofKPoints;j++){
176    for (i=0;i<NumberofEPoints;i++){
177      pYY[i] = std::log(ReducedCS[i][j]);
178    }
179    G4PenelopeInterpolator* interpolator2 = new G4PenelopeInterpolator(pX,pYY,NumberofEPoints);
180    for (i=0;i<NumberofExtendedEGrid;i++){
181      G4double ELL = ExtendedLogEnergy[i];
182      if (ELL >= pX[0]) {
183        p0[i][j] = std::exp(interpolator2->CubicSplineInterpolation(ELL));
184      }
185      else
186        {
187          G4double F1=interpolator2->CubicSplineInterpolation(pX[0]);
188          G4double FP1 = interpolator2->FirstDerivative(pX[0]);
189          p0[i][j] = std::exp(F1+FP1*(ELL-pX[0]));
190        }
191    }
192    delete interpolator2;
193  }
194 
195  //All this stuff might be useful later on
196  /*
197  G4double PDF[NumberofKPoints];
198 
199  for (i=0;i<NumberofExtendedEGrid;i++){
200    for (j=0;j<NumberofKPoints;j++){
201      PDF[j]=p0[i][j];
202    }
203    G4double Xc=0;
204    if (i<(NumberofExtendedEGrid-1)){
205      Xc=tCut/std::exp(ExtendedLogEnergy[i+1]);
206    }
207    else
208      {
209        Xc=tCut/std::exp(ExtendedLogEnergy[NumberofExtendedEGrid-1]);
210      }
211   
212    G4PenelopeInterpolator* interpolator3 = new G4PenelopeInterpolator(pK,PDF,NumberofKPoints);
213    Pbcut[i]=interpolator3->CalculateMomentum(Xc,-1);
214    delete interpolator3;
215  }
216  */
217 
218}
219
220//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
221                       
222G4double G4PenelopeBremsstrahlungContinuous::CalculateStopping(G4double energy)
223  //Stopping power expressed in MeV/mm*2
224{
225  G4double Xel=std::max(std::log(energy),ExtendedLogEnergy[0]);
226  G4double Xe=1.0+(Xel-ExtendedLogEnergy[0])*DLFC;
227  G4int Ke = (G4int) Xe; 
228  G4double Xek = Xe-Ke; 
229 
230  //Global x-section factor
231  G4double Fact=Zmat*Zmat*(energy+electron_mass_c2)*(energy+electron_mass_c2)/
232    (energy*(energy+2.0*electron_mass_c2));
233  Fact *= PositronCorrection(energy);
234
235  //Moments of the scaled bremss x-section
236  G4double wcre = tCut/energy;
237  G4double pY[NumberofKPoints];
238  G4double pK[NumberofKPoints] = {1.0e-12,0.05,0.075,0.1,0.125,0.15,0.2,0.25,
239                                  0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,
240                                  0.75,0.8,0.85,0.9,0.925,0.95,0.97,0.99,
241                                  0.995,0.999,0.9995,0.9999,0.99995,0.99999,1.0};
242
243  for (size_t i=0;i<NumberofKPoints;i++){
244    pY[i] = p0[Ke][i];
245  }
246  G4PenelopeInterpolator* interpolator1 = new G4PenelopeInterpolator(pK,pY,NumberofKPoints);
247  G4double XS1A = interpolator1->CalculateMomentum(wcre,0);
248  for (size_t k=0;k<NumberofKPoints;k++){
249    pY[k] = p0[std::min(Ke+1,(G4int) NumberofExtendedEGrid-1)][k];
250  }
251  G4PenelopeInterpolator* interpolator2 = new G4PenelopeInterpolator (pK,pY,NumberofKPoints);
252  G4double XS1B = interpolator2->CalculateMomentum(wcre,0);
253 
254  G4double XS1 = ((1.0-Xek)*XS1A+Xek*XS1B)*Fact*energy; //weighted mean between the energy bin of the grid
255
256  //This is the 2nd momentum (straggling cross section). It might be useful later on
257  /*
258    G4double XS2A = interpolator1->CalculateMomentum(wcre,1);
259    G4double XS2 = ((1.0-Xek)*XS2A+Xek*XS2B)*Fact*energy*energy;
260    G4double XS2B = interpolator2->CalculateMomentum(wcre,1);
261  */
262
263  delete interpolator1;
264  delete interpolator2;
265
266  return XS1;
267}
268
269//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270
271G4double G4PenelopeBremsstrahlungContinuous::PositronCorrection(G4double energy)
272{
273  const G4double Coeff[7]={-1.2359e-01,6.1274e-2,-3.1516e-2,7.7446e-3,-1.0595e-3,
274                           7.0568e-5,-1.8080e-6};
275  G4double correct=0;
276  if (partName == "e-") {
277    return 1.0; //no correction for electrons
278  }
279  else if (partName == "e+"){
280    G4double T=std::log(1+((1e6*energy)/(Zmat*Zmat*electron_mass_c2)));
281    for (G4int i=0;i<7;i++){
282      correct += Coeff[i]*std::pow(T,i+1);
283    }
284    correct = 1.0-std::exp(correct);
285    return correct;
286  }
287  else //neither electrons nor positrons...exception
288    {
289      G4String excep = "G4PenelopeBremmstrahlungContinuous: the particle is not e- nor e+!";
290      G4Exception(excep);
291      return 0;
292    }
293}
Note: See TracBrowser for help on using the repository browser.