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