Geant4  10.02.p01
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();
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;
89  G4DataVector p;
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 
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 
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 
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 }
static const double MeV
Definition: G4SIunits.hh:211
G4double Function(G4double x, const G4DataVector &p) const
G4String name
Definition: TRTMaterials.hh:40
G4double AverageValue(G4double xMin, G4double xMax, const G4DataVector &p) const
int G4int
Definition: G4Types.hh:78
G4double Probability(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell=0, const G4ParticleDefinition *pd=0) const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double ParameterC(G4int index) const
G4double Parameter(G4int parameterIndex, G4int Z, G4double energy) const
G4double Excitation(G4int Z, G4double kineticEnergy) const
G4double SampleEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell=0, const G4ParticleDefinition *pd=0) const
static const double eV
Definition: G4SIunits.hh:212
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4double x[NPOINTSGL]
G4double AverageEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell=0, const G4ParticleDefinition *pd=0) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
#define G4endl
Definition: G4ios.hh:61
G4double IntSpectrum(G4double xMin, G4double xMax, const G4DataVector &p) const
double G4double
Definition: G4Types.hh:76
G4RDBremsstrahlungParameters * theBRparam
G4RDeBremsstrahlungSpectrum(const G4DataVector &bins, const G4String &name)