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