Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eBremsstrahlungSpectrum.cc
Go to the documentation of this file.
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: G4eBremsstrahlungSpectrum.cc 95950 2016-03-03 10:42:48Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4eBremsstrahlungSpectrum
34 //
35 // Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
36 //
37 // Creation date: 29 September 2001
38 //
39 // Modifications:
40 // 10.10.01 MGP Revision to improve code quality and consistency with design
41 // 15.11.01 VI Update spectrum model Bethe-Haitler spectrum at high energy
42 // 30.05.02 VI Update interpolation between 2 last energy points in the
43 // parametrisation
44 // 21.02.03 V.Ivanchenko Energy bins are defined in the constructor
45 // 28.02.03 V.Ivanchenko Filename is defined in the constructor
46 //
47 // -------------------------------------------------------------------
48 
51 #include "Randomize.hh"
52 #include "G4SystemOfUnits.hh"
53 #include "G4Exp.hh"
54 
57  lowestE(0.1*eV),
58  xp(bins)
59 {
60  length = xp.size();
61  theBRparam = new G4BremsstrahlungParameters(name,length+1);
62 
63  verbose = 0;
64 }
65 
66 
68 {
69  delete theBRparam;
70 }
71 
72 
74  G4double tmin,
75  G4double tmax,
76  G4double e,
77  G4int,
78  const G4ParticleDefinition*) const
79 {
80  G4double tm = std::min(tmax, e);
81  G4double t0 = std::max(tmin, lowestE);
82  if(t0 >= tm) return 0.0;
83 
84  t0 /= e;
85  tm /= e;
86 
87  G4double z0 = lowestE/e;
89 
90  // Access parameters
91  for (size_t i=0; i<=length; i++) {
92  p.push_back(theBRparam->Parameter(i, Z, e));
93  }
94 
95  G4double x = IntSpectrum(t0, tm, p);
96  G4double y = IntSpectrum(z0, 1.0, p);
97 
98 
99  if(1 < verbose) {
100  G4cout << "tcut(MeV)= " << tmin/MeV
101  << "; tMax(MeV)= " << tmax/MeV
102  << "; t0= " << t0
103  << "; tm= " << tm
104  << "; xp[0]= " << xp[0]
105  << "; z= " << z0
106  << "; val= " << x
107  << "; nor= " << y
108  << G4endl;
109  }
110  p.clear();
111 
112  if(y > 0.0) x /= y;
113  else x = 0.0;
114  // if(x < 0.0) x = 0.0;
115 
116  return x;
117 }
118 
119 
121  G4double tmin,
122  G4double tmax,
123  G4double e,
124  G4int,
125  const G4ParticleDefinition*) const
126 {
127  G4double tm = std::min(tmax, e);
128  G4double t0 = std::max(tmin, lowestE);
129  if(t0 >= tm) return 0.0;
130 
131  t0 /= e;
132  tm /= e;
133 
134  G4double z0 = lowestE/e;
135 
136  G4DataVector p;
137 
138  // Access parameters
139  for (size_t i=0; i<=length; i++) {
140  p.push_back(theBRparam->Parameter(i, Z, e));
141  }
142 
143  G4double x = AverageValue(t0, tm, p);
144  G4double y = IntSpectrum(z0, 1.0, p);
145 
146  // Add integrant over lowest energies
147  G4double zmin = tmin/e;
148  if(zmin < t0) {
149  G4double c = std::sqrt(theBRparam->ParameterC(Z));
150  x += p[0]*(t0 - zmin - c*(std::atan(t0/c) - std::atan(zmin/c)));
151  }
152  x *= e;
153 
154  if(1 < verbose) {
155  G4cout << "tcut(MeV)= " << tmin/MeV
156  << "; tMax(MeV)= " << tmax/MeV
157  << "; e(MeV)= " << e/MeV
158  << "; t0= " << t0
159  << "; tm= " << tm
160  << "; y= " << y
161  << "; x= " << x
162  << G4endl;
163  }
164  p.clear();
165 
166  if(y > 0.0) x /= y;
167  else x = 0.0;
168  // if(x < 0.0) x = 0.0;
169 
170  return x;
171 }
172 
173 
175  G4double tmin,
176  G4double tmax,
177  G4double e,
178  G4int,
179  const G4ParticleDefinition*) const
180 {
181  G4double tm = std::min(tmax, e);
182  G4double t0 = std::max(tmin, lowestE);
183  if(t0 >= tm) return 0.0;
184 
185  t0 /= e;
186  tm /= e;
187 
188  G4DataVector p;
189 
190  for (size_t i=0; i<=length; i++) {
191  p.push_back(theBRparam->Parameter(i, Z, e));
192  }
193  G4double amaj = std::max(p[length], 1. - (p[1] - p[0])*xp[0]/(xp[1] - xp[0]) );
194 
195  G4double amax = std::log(tm);
196  G4double amin = std::log(t0);
197  G4double tgam, q, fun;
198 
199  do {
200  G4double x = amin + G4UniformRand()*(amax - amin);
201  tgam = G4Exp(x);
202  fun = Function(tgam, p);
203 
204  if(fun > amaj) {
205  G4cout << "WARNING in G4eBremsstrahlungSpectrum::SampleEnergy:"
206  << " Majoranta " << amaj
207  << " < " << fun
208  << G4endl;
209  }
210 
211  q = amaj * G4UniformRand();
212  } while (q > fun);
213 
214  tgam *= e;
215 
216  p.clear();
217 
218  return tgam;
219 }
220 
221 G4double G4eBremsstrahlungSpectrum::IntSpectrum(G4double xMin,
222  G4double xMax,
223  const G4DataVector& p) const
224 {
225  G4double x1 = std::min(xMin, xp[0]);
226  G4double x2 = std::min(xMax, xp[0]);
227  G4double sum = 0.0;
228 
229  if(x1 < x2) {
230  G4double k = (p[1] - p[0])/(xp[1] - xp[0]);
231  sum += (1. - k*xp[0])*std::log(x2/x1) + k*(x2 - x1);
232  }
233 
234  for (size_t i=0; i<length-1; i++) {
235  x1 = std::max(xMin, xp[i]);
236  x2 = std::min(xMax, xp[i+1]);
237  if(x1 < x2) {
238  G4double z1 = p[i];
239  G4double z2 = p[i+1];
240  sum += z2 - z1 + std::log(x2/x1)*(z1*x2 - z2*x1)/(x2 - x1);
241  }
242  }
243  if(sum < 0.0) sum = 0.0;
244  return sum;
245 }
246 
247 G4double G4eBremsstrahlungSpectrum::AverageValue(G4double xMin,
248  G4double xMax,
249  const G4DataVector& p) const
250 {
251  G4double x1 = std::min(xMin, xp[0]);
252  G4double x2 = std::min(xMax, xp[0]);
253  G4double z1 = x1;
254  G4double z2 = x2;
255  G4double sum = 0.0;
256 
257  if(x1 < x2) {
258  G4double k = (p[1] - p[0])/(xp[1] - xp[0]);
259  sum += (z2 - z1)*(1. - k*xp[0]);
260  z1 *= x1;
261  z2 *= x2;
262  sum += 0.5*k*(z2 - z1);
263  }
264 
265  for (size_t i=0; i<length-1; i++) {
266  x1 = std::max(xMin, xp[i]);
267  x2 = std::min(xMax, xp[i+1]);
268  if(x1 < x2) {
269  z1 = p[i];
270  z2 = p[i+1];
271  sum += 0.5*(z2 - z1)*(x2 + x1) + z1*x2 - z2*x1;
272  }
273  }
274  if(sum < 0.0) sum = 0.0;
275  return sum;
276 }
277 
278 G4double G4eBremsstrahlungSpectrum::Function(G4double x,
279  const G4DataVector& p) const
280 {
281  G4double f = 0.0;
282 
283  if(x <= xp[0]) {
284  f = p[0] + (p[1] - p[0])*(x - xp[0])/(xp[1] - xp[0]);
285 
286  } else {
287 
288  for (size_t i=0; i<length-1; i++) {
289 
290  if(x <= xp[i+1]) {
291  f = p[i] + (p[i+1] - p[i])*(x - xp[i])/(xp[i+1] - xp[i]);
292  break;
293  }
294  }
295  }
296  return f;
297 }
298 
300 { theBRparam->PrintData(); }
301 
303 {
304  return 0.0;
305 }
306 
308  G4int, // Z,
309  const G4ParticleDefinition*) const
310 {
311  return kineticEnergy;
312 }
const XML_Char * name
Definition: expat.h:151
const char * p
Definition: xmltok.h:285
G4double Excitation(G4int Z, G4double kineticEnergy) const
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
G4eBremsstrahlungSpectrum(const G4DataVector &bins, const G4String &name)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
static constexpr double eV
Definition: G4SIunits.hh:215
G4double Parameter(G4int parameterIndex, G4int Z, G4double energy) const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double SampleEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell=0, const G4ParticleDefinition *pd=0) const
G4double AverageEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell=0, const G4ParticleDefinition *pd=0) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double Probability(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell=0, const G4ParticleDefinition *pd=0) const
double G4double
Definition: G4Types.hh:76
G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
tuple c
Definition: test.py:13
G4double ParameterC(G4int index) const