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

Last change on this file since 968 was 961, checked in by garnier, 17 years ago

update processes
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.11 2008/12/15 09:23:06 pandola Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
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.