Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SPSEneDistribution.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 //
27 //
28 // MODULE: G4SPSEneDistribution.cc
29 //
30 // Version: 1.0
31 // Date: 5/02/04
32 // Author: Fan Lei
33 // Organisation: QinetiQ ltd.
34 // Customer: ESA/ESTEC
35 //
37 //
38 // CHANGE HISTORY
39 // --------------
40 //
41 //
42 // Version 1.0, 05/02/2004, Fan Lei, Created.
43 // Based on the G4GeneralParticleSource class in Geant4 v6.0
44 //
46 //
47 
48 #include "G4SPSEneDistribution.hh"
49 
50 #include "G4SystemOfUnits.hh"
51 #include "Randomize.hh"
52 
54  : particle_definition(0), eneRndm(0), Splinetemp(0)
55 {
56  //
57  // Initialise all variables
58  particle_energy = 1.0 * MeV;
59 
60  EnergyDisType = "Mono";
61  weight = 1.;
62  MonoEnergy = 1 * MeV;
63  Emin = 0.;
64  Emax = 1.e30;
65  alpha = 0.;
66  biasalpha = 0.;
67  prob_norm = 1.0;
68  Ezero = 0.;
69  SE = 0.;
70  Temp = 0.;
71  grad = 0.;
72  cept = 0.;
73  Biased = false; // not biased
74  EnergySpec = true; // true - energy spectra, false - momentum spectra
75  DiffSpec = true; // true - differential spec, false integral spec
76  IntType = "NULL"; // Interpolation type
77  IPDFEnergyExist = false;
78  IPDFArbExist = false;
79 
80  ArbEmin = 0.;
81  ArbEmax = 1.e30;
82 
83  verbosityLevel = 0;
84 
85 }
86 
88 }
89 
91  EnergyDisType = DisType;
92  if (EnergyDisType == "User") {
93  UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
94  IPDFEnergyExist = false;
95  } else if (EnergyDisType == "Arb") {
96  ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
97  IPDFArbExist = false;
98  } else if (EnergyDisType == "Epn") {
99  UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
100  IPDFEnergyExist = false;
101  EpnEnergyH = ZeroPhysVector;
102  }
103 }
104 
106  Emin = emi;
107 }
108 
110  Emax = ema;
111 }
112 
114  MonoEnergy = menergy;
115 }
116 
118  SE = e;
119 }
121  alpha = alp;
122 }
123 
125  biasalpha = alp;
126  Biased = true;
127 }
128 
130  Temp = tem;
131 }
132 
134  Ezero = eze;
135 }
136 
138  grad = gr;
139 }
140 
142  cept = c;
143 }
144 
146  G4double ehi, val;
147  ehi = input.x();
148  val = input.y();
149  if (verbosityLevel > 1) {
150  G4cout << "In UserEnergyHisto" << G4endl;
151  G4cout << " " << ehi << " " << val << G4endl;
152  }
153  UDefEnergyH.InsertValues(ehi, val);
154  Emax = ehi;
155 }
156 
158  G4double ehi, val;
159  ehi = input.x();
160  val = input.y();
161  if (verbosityLevel > 1) {
162  G4cout << "In ArbEnergyHisto" << G4endl;
163  G4cout << " " << ehi << " " << val << G4endl;
164  }
165  ArbEnergyH.InsertValues(ehi, val);
166 }
167 
169  std::ifstream infile(filename, std::ios::in);
170  if (!infile)
171  G4Exception("G4SPSEneDistribution::ArbEnergyHistoFile",
172  "Event0301",FatalException,
173  "Unable to open the histo ASCII file");
174  G4double ehi, val;
175  while (infile >> ehi >> val) {
176  ArbEnergyH.InsertValues(ehi, val);
177  }
178 }
179 
181  G4double ehi, val;
182  ehi = input.x();
183  val = input.y();
184  if (verbosityLevel > 1) {
185  G4cout << "In EpnEnergyHisto" << G4endl;
186  G4cout << " " << ehi << " " << val << G4endl;
187  }
188  EpnEnergyH.InsertValues(ehi, val);
189  Emax = ehi;
190  Epnflag = true;
191 }
192 
194  if (EnergyDisType == "Cdg")
195  CalculateCdgSpectrum();
196  else if (EnergyDisType == "Bbody")
197  CalculateBbodySpectrum();
198 }
199 
200 void G4SPSEneDistribution::CalculateCdgSpectrum() {
201  // This uses the spectrum from The INTEGRAL Mass Model (TIMM)
202  // to generate a Cosmic Diffuse X/gamma ray spectrum.
203  G4double pfact[2] = { 8.5, 112 };
204  G4double spind[2] = { 1.4, 2.3 };
205  G4double ene_line[3] = { 1. * keV, 18. * keV, 1E6 * keV };
206  G4int n_par;
207 
208  ene_line[0] = Emin;
209  if (Emin < 18 * keV) {
210  n_par = 2;
211  ene_line[2] = Emax;
212  if (Emax < 18 * keV) {
213  n_par = 1;
214  ene_line[1] = Emax;
215  }
216  } else {
217  n_par = 1;
218  pfact[0] = 112.;
219  spind[0] = 2.3;
220  ene_line[1] = Emax;
221  }
222 
223  // Create a cumulative histogram.
224  CDGhist[0] = 0.;
225  G4double omalpha;
226  G4int i = 0;
227 
228  while (i < n_par) {
229  omalpha = 1. - spind[i];
230  CDGhist[i + 1] = CDGhist[i] + (pfact[i] / omalpha) * (std::pow(
231  ene_line[i + 1] / keV, omalpha) - std::pow(ene_line[i] / keV,
232  omalpha));
233  i++;
234  }
235 
236  // Normalise histo and divide by 1000 to make MeV.
237  i = 0;
238  while (i < n_par) {
239  CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[n_par];
240  // G4cout << CDGhist[i] << CDGhist[n_par] << G4endl;
241  i++;
242  }
243 }
244 
245 void G4SPSEneDistribution::CalculateBbodySpectrum() {
246  // create bbody spectrum
247  // Proved very hard to integrate indefinitely, so different
248  // method. User inputs emin, emax and T. These are used to
249  // create a 10,000 bin histogram.
250  // Use photon density spectrum = 2 nu**2/c**2 * (std::exp(h nu/kT)-1)
251  // = 2 E**2/h**2c**2 times the exponential
252  G4double erange = Emax - Emin;
253  G4double steps = erange / 10000.;
254  G4double Bbody_y[10000];
255  G4double k = 8.6181e-11; //Boltzmann const in MeV/K
256  G4double h = 4.1362e-21; // Plancks const in MeV s
257  G4double c = 3e8; // Speed of light
258  G4double h2 = h * h;
259  G4double c2 = c * c;
260  G4int count = 0;
261  G4double sum = 0.;
262  BBHist[0] = 0.;
263  while (count < 10000) {
264  Bbody_x[count] = Emin + G4double(count * steps);
265  Bbody_y[count] = (2. * std::pow(Bbody_x[count], 2.)) / (h2 * c2
266  * (std::exp(Bbody_x[count] / (k * Temp)) - 1.));
267  sum = sum + Bbody_y[count];
268  BBHist[count + 1] = BBHist[count] + Bbody_y[count];
269  count++;
270  }
271 
272  Bbody_x[10000] = Emax;
273  // Normalise cumulative histo.
274  count = 0;
275  while (count < 10001) {
276  BBHist[count] = BBHist[count] / sum;
277  count++;
278  }
279 }
280 
282  // Allows user to specifiy spectrum is momentum
283  EnergySpec = value; // false if momentum
284  if (verbosityLevel > 1)
285  G4cout << "EnergySpec has value " << EnergySpec << G4endl;
286 }
287 
289  // Allows user to specify integral or differential spectra
290  DiffSpec = value; // true = differential, false = integral
291  if (verbosityLevel > 1)
292  G4cout << "Diffspec has value " << DiffSpec << G4endl;
293 }
294 
296  if (EnergyDisType != "Arb")
297  G4cout << "Error: this is for arbitrary distributions" << G4endl;
298  IntType = IType;
299  ArbEmax = ArbEnergyH.GetMaxLowEdgeEnergy();
300  ArbEmin = ArbEnergyH.GetMinLowEdgeEnergy();
301 
302  // Now interpolate points
303  if (IntType == "Lin")
304  LinearInterpolation();
305  if (IntType == "Log")
306  LogInterpolation();
307  if (IntType == "Exp")
308  ExpInterpolation();
309  if (IntType == "Spline")
310  SplineInterpolation();
311 }
312 
313 void G4SPSEneDistribution::LinearInterpolation() {
314  // Method to do linear interpolation on the Arb points
315  // Calculate equation of each line segment, max 1024.
316  // Calculate Area under each segment
317  // Create a cumulative array which is then normalised Arb_Cum_Area
318 
319  G4double Area_seg[1024]; // Stores area under each segment
320  G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
321  G4int i, count;
322  G4int maxi = ArbEnergyH.GetVectorLength();
323  for (i = 0; i < maxi; i++) {
324  Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
325  Arb_y[i] = ArbEnergyH(size_t(i));
326  }
327  // Points are now in x,y arrays. If the spectrum is integral it has to be
328  // made differential and if momentum it has to be made energy.
329  if (DiffSpec == false) {
330  // Converts integral point-wise spectra to Differential
331  for (count = 0; count < maxi - 1; count++) {
332  Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
333  / (Arb_x[count + 1] - Arb_x[count]);
334  }
335  maxi--;
336  }
337  //
338  if (EnergySpec == false) {
339  // change currently stored values (emin etc) which are actually momenta
340  // to energies.
341  if (particle_definition == NULL)
342  G4cout << "Error: particle not defined" << G4endl;
343  else {
344  // Apply Energy**2 = p**2c**2 + m0**2c**4
345  // p should be entered as E/c i.e. without the division by c
346  // being done - energy equivalent.
347  G4double mass = particle_definition->GetPDGMass();
348  // convert point to energy unit and its value to per energy unit
349  G4double total_energy;
350  for (count = 0; count < maxi; count++) {
351  total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
352  * mass)); // total energy
353 
354  Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
355  Arb_x[count] = total_energy - mass; // kinetic energy
356  }
357  }
358  }
359  //
360  i = 1;
361  Arb_grad[0] = 0.;
362  Arb_cept[0] = 0.;
363  Area_seg[0] = 0.;
364  Arb_Cum_Area[0] = 0.;
365  while (i < maxi) {
366  // calc gradient and intercept for each segment
367  Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[i - 1]);
368  if (verbosityLevel == 2)
369  G4cout << Arb_grad[i] << G4endl;
370  if (Arb_grad[i] > 0.) {
371  if (verbosityLevel == 2)
372  G4cout << "Arb_grad is positive" << G4endl;
373  Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]);
374  } else if (Arb_grad[i] < 0.) {
375  if (verbosityLevel == 2)
376  G4cout << "Arb_grad is negative" << G4endl;
377  Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]);
378  } else {
379  if (verbosityLevel == 2)
380  G4cout << "Arb_grad is 0." << G4endl;
381  Arb_cept[i] = Arb_y[i];
382  }
383 
384  Area_seg[i] = ((Arb_grad[i] / 2) * (Arb_x[i] * Arb_x[i] - Arb_x[i - 1]
385  * Arb_x[i - 1]) + Arb_cept[i] * (Arb_x[i] - Arb_x[i - 1]));
386  Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
387  sum = sum + Area_seg[i];
388  if (verbosityLevel == 2)
389  G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i]
390  << G4endl;
391  i++;
392  }
393 
394  i = 0;
395  while (i < maxi) {
396  Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation
397  IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
398  i++;
399  }
400 
401  // now scale the ArbEnergyH, needed by Probability()
402  ArbEnergyH.ScaleVector(1., 1./sum);
403 
404  if (verbosityLevel >= 1) {
405  G4cout << "Leaving LinearInterpolation" << G4endl;
406  ArbEnergyH.DumpValues();
407  IPDFArbEnergyH.DumpValues();
408  }
409 }
410 
411 void G4SPSEneDistribution::LogInterpolation() {
412  // Interpolation based on Logarithmic equations
413  // Generate equations of line segments
414  // y = Ax**alpha => log y = alpha*logx + logA
415  // Find area under line segments
416  // create normalised, cumulative array Arb_Cum_Area
417  G4double Area_seg[1024]; // Stores area under each segment
418  G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
419  G4int i, count;
420  G4int maxi = ArbEnergyH.GetVectorLength();
421  for (i = 0; i < maxi; i++) {
422  Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
423  Arb_y[i] = ArbEnergyH(size_t(i));
424  }
425  // Points are now in x,y arrays. If the spectrum is integral it has to be
426  // made differential and if momentum it has to be made energy.
427  if (DiffSpec == false) {
428  // Converts integral point-wise spectra to Differential
429  for (count = 0; count < maxi - 1; count++) {
430  Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
431  / (Arb_x[count + 1] - Arb_x[count]);
432  }
433  maxi--;
434  }
435  //
436  if (EnergySpec == false) {
437  // change currently stored values (emin etc) which are actually momenta
438  // to energies.
439  if (particle_definition == NULL)
440  G4cout << "Error: particle not defined" << G4endl;
441  else {
442  // Apply Energy**2 = p**2c**2 + m0**2c**4
443  // p should be entered as E/c i.e. without the division by c
444  // being done - energy equivalent.
445  G4double mass = particle_definition->GetPDGMass();
446  // convert point to energy unit and its value to per energy unit
447  G4double total_energy;
448  for (count = 0; count < maxi; count++) {
449  total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
450  * mass)); // total energy
451 
452  Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
453  Arb_x[count] = total_energy - mass; // kinetic energy
454  }
455  }
456  }
457  //
458  i = 1;
459  Arb_alpha[0] = 0.;
460  Arb_Const[0] = 0.;
461  Area_seg[0] = 0.;
462  Arb_Cum_Area[0] = 0.;
463  if (Arb_x[0] <= 0. || Arb_y[0] <= 0.) {
464  G4cout << "You should not use log interpolation with points <= 0."
465  << G4endl;
466  G4cout << "These will be changed to 1e-20, which may cause problems"
467  << G4endl;
468  if (Arb_x[0] <= 0.)
469  Arb_x[0] = 1e-20;
470  if (Arb_y[0] <= 0.)
471  Arb_y[0] = 1e-20;
472  }
473 
474  G4double alp;
475  while (i < maxi) {
476  // Incase points are negative or zero
477  if (Arb_x[i] <= 0. || Arb_y[i] <= 0.) {
478  G4cout << "You should not use log interpolation with points <= 0."
479  << G4endl;
480  G4cout
481  << "These will be changed to 1e-20, which may cause problems"
482  << G4endl;
483  if (Arb_x[i] <= 0.)
484  Arb_x[i] = 1e-20;
485  if (Arb_y[i] <= 0.)
486  Arb_y[i] = 1e-20;
487  }
488 
489  Arb_alpha[i] = (std::log10(Arb_y[i]) - std::log10(Arb_y[i - 1]))
490  / (std::log10(Arb_x[i]) - std::log10(Arb_x[i - 1]));
491  Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[i], Arb_alpha[i]));
492  alp = Arb_alpha[i] + 1;
493  if (alp == 0.) {
494  Area_seg[i] = Arb_Const[i] * (std::log(Arb_x[i]) - std::log(Arb_x[i - 1]));
495  } else {
496  Area_seg[i] = (Arb_Const[i] / alp) * (std::pow(Arb_x[i], alp)
497  - std::pow(Arb_x[i - 1], alp));
498  }
499  sum = sum + Area_seg[i];
500  Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
501  if (verbosityLevel == 2)
502  G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl;
503  i++;
504  }
505 
506  i = 0;
507  while (i < maxi) {
508  Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
509  IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
510  i++;
511  }
512 
513  // now scale the ArbEnergyH, needed by Probability()
514  ArbEnergyH.ScaleVector(1., 1./sum);
515 
516  if (verbosityLevel >= 1)
517  G4cout << "Leaving LogInterpolation " << G4endl;
518 }
519 
520 void G4SPSEneDistribution::ExpInterpolation() {
521  // Interpolation based on Exponential equations
522  // Generate equations of line segments
523  // y = Ae**-(x/e0) => ln y = -x/e0 + lnA
524  // Find area under line segments
525  // create normalised, cumulative array Arb_Cum_Area
526  G4double Area_seg[1024]; // Stores area under each segment
527  G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
528  G4int i, count;
529  G4int maxi = ArbEnergyH.GetVectorLength();
530  for (i = 0; i < maxi; i++) {
531  Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
532  Arb_y[i] = ArbEnergyH(size_t(i));
533  }
534  // Points are now in x,y arrays. If the spectrum is integral it has to be
535  // made differential and if momentum it has to be made energy.
536  if (DiffSpec == false) {
537  // Converts integral point-wise spectra to Differential
538  for (count = 0; count < maxi - 1; count++) {
539  Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
540  / (Arb_x[count + 1] - Arb_x[count]);
541  }
542  maxi--;
543  }
544  //
545  if (EnergySpec == false) {
546  // change currently stored values (emin etc) which are actually momenta
547  // to energies.
548  if (particle_definition == NULL)
549  G4cout << "Error: particle not defined" << G4endl;
550  else {
551  // Apply Energy**2 = p**2c**2 + m0**2c**4
552  // p should be entered as E/c i.e. without the division by c
553  // being done - energy equivalent.
554  G4double mass = particle_definition->GetPDGMass();
555  // convert point to energy unit and its value to per energy unit
556  G4double total_energy;
557  for (count = 0; count < maxi; count++) {
558  total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
559  * mass)); // total energy
560 
561  Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
562  Arb_x[count] = total_energy - mass; // kinetic energy
563  }
564  }
565  }
566  //
567  i = 1;
568  Arb_ezero[0] = 0.;
569  Arb_Const[0] = 0.;
570  Area_seg[0] = 0.;
571  Arb_Cum_Area[0] = 0.;
572  while (i < maxi) {
573  G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i - 1]);
574  if (test > 0. || test < 0.) {
575  Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1]) / (std::log(Arb_y[i])
576  - std::log(Arb_y[i - 1]));
577  Arb_Const[i] = Arb_y[i] / (std::exp(-Arb_x[i] / Arb_ezero[i]));
578  Area_seg[i] = -(Arb_Const[i] * Arb_ezero[i]) * (std::exp(-Arb_x[i]
579  / Arb_ezero[i]) - std::exp(-Arb_x[i - 1] / Arb_ezero[i]));
580  } else {
581  G4cout << "Flat line segment: problem" << G4endl;
582  Arb_ezero[i] = 0.;
583  Arb_Const[i] = 0.;
584  Area_seg[i] = 0.;
585  }
586  sum = sum + Area_seg[i];
587  Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
588  if (verbosityLevel == 2)
589  G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl;
590  i++;
591  }
592 
593  i = 0;
594  while (i < maxi) {
595  Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
596  IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
597  i++;
598  }
599 
600  // now scale the ArbEnergyH, needed by Probability()
601  ArbEnergyH.ScaleVector(1., 1./sum);
602 
603  if (verbosityLevel >= 1)
604  G4cout << "Leaving ExpInterpolation " << G4endl;
605 }
606 
607 void G4SPSEneDistribution::SplineInterpolation() {
608  // Interpolation using Splines.
609  // Create Normalised arrays, make x 0->1 and y hold
610  // the function (Energy)
611  //
612  // Current method based on the above will not work in all cases.
613  // new method is implemented below.
614 
615  G4double sum, Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
616  G4int i, count;
617 
618  G4int maxi = ArbEnergyH.GetVectorLength();
619  for (i = 0; i < maxi; i++) {
620  Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
621  Arb_y[i] = ArbEnergyH(size_t(i));
622  }
623  // Points are now in x,y arrays. If the spectrum is integral it has to be
624  // made differential and if momentum it has to be made energy.
625  if (DiffSpec == false) {
626  // Converts integral point-wise spectra to Differential
627  for (count = 0; count < maxi - 1; count++) {
628  Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
629  / (Arb_x[count + 1] - Arb_x[count]);
630  }
631  maxi--;
632  }
633  //
634  if (EnergySpec == false) {
635  // change currently stored values (emin etc) which are actually momenta
636  // to energies.
637  if (particle_definition == NULL)
638  G4Exception("G4SPSEneDistribution::SplineInterpolation",
639  "Event0302",FatalException,
640  "Error: particle not defined");
641  else {
642  // Apply Energy**2 = p**2c**2 + m0**2c**4
643  // p should be entered as E/c i.e. without the division by c
644  // being done - energy equivalent.
645  G4double mass = particle_definition->GetPDGMass();
646  // convert point to energy unit and its value to per energy unit
647  G4double total_energy;
648  for (count = 0; count < maxi; count++) {
649  total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
650  * mass)); // total energy
651 
652  Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
653  Arb_x[count] = total_energy - mass; // kinetic energy
654  }
655  }
656  }
657 
658  //
659  i = 1;
660  Arb_Cum_Area[0] = 0.;
661  sum = 0.;
662  Splinetemp = new G4DataInterpolation(Arb_x, Arb_y, maxi, 0., 0.);
663  G4double ei[101],prob[101];
664  while (i < maxi) {
665  // 100 step per segment for the integration of area
666  G4double de = (Arb_x[i] - Arb_x[i - 1])/100.;
667  G4double area = 0.;
668 
669  for (count = 0; count < 101; count++) {
670  ei[count] = Arb_x[i - 1] + de*count ;
671  prob[count] = Splinetemp->CubicSplineInterpolation(ei[count]);
672  if (prob[count] < 0.) {
674  ED << "Warning: G4DataInterpolation returns value < 0 " << prob[count] <<" "<<ei[count]<< G4endl;
675  G4Exception("G4SPSEneDistribution::SplineInterpolation","Event0303",
676  FatalException,ED);
677  }
678  area += prob[count]*de;
679  }
680  Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + area;
681  sum += area;
682 
683  prob[0] = prob[0]/(area/de);
684  for (count = 1; count < 100; count++)
685  prob[count] = prob[count-1] + prob[count]/(area/de);
686 
687  SplineInt[i] = new G4DataInterpolation(prob, ei, 101, 0., 0.);
688  // note i start from 1!
689  i++;
690  }
691  i = 0;
692  while (i < maxi) {
693  Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation
694  IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
695  i++;
696  }
697  // now scale the ArbEnergyH, needed by Probability()
698  ArbEnergyH.ScaleVector(1., 1./sum);
699 
700  if (verbosityLevel > 0)
701  G4cout << "Leaving SplineInterpolation " << G4endl;
702 }
703 
704 void G4SPSEneDistribution::GenerateMonoEnergetic() {
705  // Method to generate MonoEnergetic particles.
706  particle_energy = MonoEnergy;
707 }
708 
709 void G4SPSEneDistribution::GenerateGaussEnergies() {
710  // Method to generate Gaussian particles.
711  particle_energy = G4RandGauss::shoot(MonoEnergy,SE);
712  if (particle_energy < 0) particle_energy = 0.;
713 }
714 
715 void G4SPSEneDistribution::GenerateLinearEnergies(G4bool bArb = false) {
716  G4double rndm;
717  G4double emaxsq = std::pow(Emax, 2.); //Emax squared
718  G4double eminsq = std::pow(Emin, 2.); //Emin squared
719  G4double intersq = std::pow(cept, 2.); //cept squared
720 
721  if (bArb)
722  rndm = G4UniformRand();
723  else
724  rndm = eneRndm->GenRandEnergy();
725 
726  G4double bracket = ((grad / 2.) * (emaxsq - eminsq) + cept * (Emax - Emin));
727  bracket = bracket * rndm;
728  bracket = bracket + (grad / 2.) * eminsq + cept * Emin;
729  // Now have a quad of form m/2 E**2 + cE - bracket = 0
730  bracket = -bracket;
731  // G4cout << "BRACKET" << bracket << G4endl;
732  if (grad != 0.) {
733  G4double sqbrack = (intersq - 4 * (grad / 2.) * (bracket));
734  // G4cout << "SQBRACK" << sqbrack << G4endl;
735  sqbrack = std::sqrt(sqbrack);
736  G4double root1 = -cept + sqbrack;
737  root1 = root1 / (2. * (grad / 2.));
738 
739  G4double root2 = -cept - sqbrack;
740  root2 = root2 / (2. * (grad / 2.));
741 
742  // G4cout << root1 << " roots " << root2 << G4endl;
743 
744  if (root1 > Emin && root1 < Emax)
745  particle_energy = root1;
746  if (root2 > Emin && root2 < Emax)
747  particle_energy = root2;
748  } else if (grad == 0.)
749  // have equation of form cE - bracket =0
750  particle_energy = bracket / cept;
751 
752  if (particle_energy < 0.)
753  particle_energy = -particle_energy;
754 
755  if (verbosityLevel >= 1)
756  G4cout << "Energy is " << particle_energy << G4endl;
757 }
758 
759 void G4SPSEneDistribution::GeneratePowEnergies(G4bool bArb = false) {
760  // Method to generate particle energies distributed as
761  // a power-law
762 
763  G4double rndm;
764  G4double emina, emaxa;
765 
766  emina = std::pow(Emin, alpha + 1);
767  emaxa = std::pow(Emax, alpha + 1);
768 
769  if (bArb)
770  rndm = G4UniformRand();
771  else
772  rndm = eneRndm->GenRandEnergy();
773 
774  if (alpha != -1.) {
775  particle_energy = ((rndm * (emaxa - emina)) + emina);
776  particle_energy = std::pow(particle_energy, (1. / (alpha + 1.)));
777  } else {
778  particle_energy = (std::log(Emin) + rndm * (std::log(Emax) - std::log(
779  Emin)));
780  particle_energy = std::exp(particle_energy);
781  }
782  if (verbosityLevel >= 1)
783  G4cout << "Energy is " << particle_energy << G4endl;
784 }
785 
786 void G4SPSEneDistribution::GenerateBiasPowEnergies() {
787  // Method to generate particle energies distributed as
788  // in biased power-law and calculate its weight
789 
790  G4double rndm;
791  G4double emina, emaxa, emin, emax;
792 
793  G4double normal = 1. ;
794 
795  emin = Emin;
796  emax = Emax;
797  // if (EnergyDisType == "Arb") {
798  // emin = ArbEmin;
799  // emax = ArbEmax;
800  //}
801 
802  rndm = eneRndm->GenRandEnergy();
803 
804  if (biasalpha != -1.) {
805  emina = std::pow(emin, biasalpha + 1);
806  emaxa = std::pow(emax, biasalpha + 1);
807  particle_energy = ((rndm * (emaxa - emina)) + emina);
808  particle_energy = std::pow(particle_energy, (1. / (biasalpha + 1.)));
809  normal = 1./(1+biasalpha) * (emaxa - emina);
810  } else {
811  particle_energy = (std::log(emin) + rndm * (std::log(emax) - std::log(
812  emin)));
813  particle_energy = std::exp(particle_energy);
814  normal = std::log(emax) - std::log(emin) ;
815  }
816  weight = GetProbability(particle_energy) / (std::pow(particle_energy,biasalpha)/normal);
817 
818  if (verbosityLevel >= 1)
819  G4cout << "Energy is " << particle_energy << G4endl;
820 }
821 
822 void G4SPSEneDistribution::GenerateExpEnergies(G4bool bArb = false) {
823  // Method to generate particle energies distributed according
824  // to an exponential curve.
825  G4double rndm;
826 
827  if (bArb)
828  rndm = G4UniformRand();
829  else
830  rndm = eneRndm->GenRandEnergy();
831 
832  particle_energy = -Ezero * (std::log(rndm * (std::exp(-Emax / Ezero)
833  - std::exp(-Emin / Ezero)) + std::exp(-Emin / Ezero)));
834  if (verbosityLevel >= 1)
835  G4cout << "Energy is " << particle_energy << G4endl;
836 }
837 
838 void G4SPSEneDistribution::GenerateBremEnergies() {
839  // Method to generate particle energies distributed according
840  // to a Bremstrahlung equation of
841  // form I = const*((kT)**1/2)*E*(e**(-E/kT))
842 
843  G4double rndm;
844  rndm = eneRndm->GenRandEnergy();
845  G4double expmax, expmin, k;
846 
847  k = 8.6181e-11; // Boltzmann's const in MeV/K
848  G4double ksq = std::pow(k, 2.); // k squared
849  G4double Tsq = std::pow(Temp, 2.); // Temp squared
850 
851  expmax = std::exp(-Emax / (k * Temp));
852  expmin = std::exp(-Emin / (k * Temp));
853 
854  // If either expmax or expmin are zero then this will cause problems
855  // Most probably this will be because T is too low or E is too high
856 
857  if (expmax == 0.)
858  G4cout << "*****EXPMAX=0. Choose different E's or Temp" << G4endl;
859  if (expmin == 0.)
860  G4cout << "*****EXPMIN=0. Choose different E's or Temp" << G4endl;
861 
862  G4double tempvar = rndm * ((-k) * Temp * (Emax * expmax - Emin * expmin)
863  - (ksq * Tsq * (expmax - expmin)));
864 
865  G4double bigc = (tempvar - k * Temp * Emin * expmin - ksq * Tsq * expmin)
866  / (-k * Temp);
867 
868  // This gives an equation of form: Ee(-E/kT) + kTe(-E/kT) - C =0
869  // Solve this iteratively, step from Emin to Emax in 1000 steps
870  // and take the best solution.
871 
872  G4double erange = Emax - Emin;
873  G4double steps = erange / 1000.;
874  G4int i;
875  G4double etest, diff, err;
876 
877  err = 100000.;
878 
879  for (i = 1; i < 1000; i++) {
880  etest = Emin + (i - 1) * steps;
881 
882  diff = etest * (std::exp(-etest / (k * Temp))) + k * Temp * (std::exp(
883  -etest / (k * Temp))) - bigc;
884 
885  if (diff < 0.)
886  diff = -diff;
887 
888  if (diff < err) {
889  err = diff;
890  particle_energy = etest;
891  }
892  }
893  if (verbosityLevel >= 1)
894  G4cout << "Energy is " << particle_energy << G4endl;
895 }
896 
897 void G4SPSEneDistribution::GenerateBbodyEnergies() {
898  // BBody_x holds Energies, and BBHist holds the cumulative histo.
899  // binary search to find correct bin then lin interpolation.
900  // Use the earlier defined histogram + RandGeneral method to generate
901  // random numbers following the histos distribution.
902  G4double rndm;
903  G4int nabove, nbelow = 0, middle;
904  nabove = 10001;
905  rndm = eneRndm->GenRandEnergy();
906 
907  // Binary search to find bin that rndm is in
908  while (nabove - nbelow > 1) {
909  middle = (nabove + nbelow) / 2;
910  if (rndm == BBHist[middle])
911  break;
912  if (rndm < BBHist[middle])
913  nabove = middle;
914  else
915  nbelow = middle;
916  }
917 
918  // Now interpolate in that bin to find the correct output value.
919  G4double x1, x2, y1, y2, t, q;
920  x1 = Bbody_x[nbelow];
921  x2 = Bbody_x[nbelow + 1];
922  y1 = BBHist[nbelow];
923  y2 = BBHist[nbelow + 1];
924  t = (y2 - y1) / (x2 - x1);
925  q = y1 - t * x1;
926 
927  particle_energy = (rndm - q) / t;
928 
929  if (verbosityLevel >= 1) {
930  G4cout << "Energy is " << particle_energy << G4endl;
931  }
932 }
933 
934 void G4SPSEneDistribution::GenerateCdgEnergies() {
935  // Gen random numbers, compare with values in cumhist
936  // to find appropriate part of spectrum and then
937  // generate energy in the usual inversion way.
938  // G4double pfact[2] = {8.5, 112};
939  // G4double spind[2] = {1.4, 2.3};
940  // G4double ene_line[3] = {1., 18., 1E6};
941  G4double rndm, rndm2;
942  G4double ene_line[3];
943  G4double omalpha[2];
944  if (Emin < 18 * keV && Emax < 18 * keV) {
945  omalpha[0] = 1. - 1.4;
946  ene_line[0] = Emin;
947  ene_line[1] = Emax;
948  }
949  if (Emin < 18 * keV && Emax > 18 * keV) {
950  omalpha[0] = 1. - 1.4;
951  omalpha[1] = 1. - 2.3;
952  ene_line[0] = Emin;
953  ene_line[1] = 18. * keV;
954  ene_line[2] = Emax;
955  }
956  if (Emin > 18 * keV) {
957  omalpha[0] = 1. - 2.3;
958  ene_line[0] = Emin;
959  ene_line[1] = Emax;
960  }
961  rndm = eneRndm->GenRandEnergy();
962  rndm2 = eneRndm->GenRandEnergy();
963 
964  G4int i = 0;
965  while (rndm >= CDGhist[i]) {
966  i++;
967  }
968  // Generate final energy.
969  particle_energy = (std::pow(ene_line[i - 1], omalpha[i - 1]) + (std::pow(
970  ene_line[i], omalpha[i - 1]) - std::pow(ene_line[i - 1], omalpha[i
971  - 1])) * rndm2);
972  particle_energy = std::pow(particle_energy, (1. / omalpha[i - 1]));
973 
974  if (verbosityLevel >= 1)
975  G4cout << "Energy is " << particle_energy << G4endl;
976 }
977 
978 void G4SPSEneDistribution::GenUserHistEnergies() {
979  // Histograms are DIFFERENTIAL.
980  // G4cout << "In GenUserHistEnergies " << G4endl;
981  if (IPDFEnergyExist == false) {
982  G4int ii;
983  G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
984  G4double bins[1024], vals[1024], sum;
985  sum = 0.;
986 
987  if ((EnergySpec == false) && (particle_definition == NULL))
988  G4cout << "Error: particle definition is NULL" << G4endl;
989 
990  if (maxbin > 1024) {
991  G4cout << "Maxbin > 1024" << G4endl;
992  G4cout << "Setting maxbin to 1024, other bins are lost" << G4endl;
993  }
994 
995  if (DiffSpec == false)
996  G4cout << "Histograms are Differential!!! " << G4endl;
997  else {
998  bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
999  vals[0] = UDefEnergyH(size_t(0));
1000  sum = vals[0];
1001  for (ii = 1; ii < maxbin; ii++) {
1002  bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
1003  vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1];
1004  sum = sum + UDefEnergyH(size_t(ii));
1005  }
1006  }
1007 
1008  if (EnergySpec == false) {
1009  G4double mass = particle_definition->GetPDGMass();
1010  // multiply the function (vals) up by the bin width
1011  // to make the function counts/s (i.e. get rid of momentum
1012  // dependence).
1013  for (ii = 1; ii < maxbin; ii++) {
1014  vals[ii] = vals[ii] * (bins[ii] - bins[ii - 1]);
1015  }
1016  // Put energy bins into new histo, plus divide by energy bin width
1017  // to make evals counts/s/energy
1018  for (ii = 0; ii < maxbin; ii++) {
1019  bins[ii] = std::sqrt((bins[ii] * bins[ii]) + (mass * mass))
1020  - mass; //kinetic energy
1021  }
1022  for (ii = 1; ii < maxbin; ii++) {
1023  vals[ii] = vals[ii] / (bins[ii] - bins[ii - 1]);
1024  }
1025  sum = vals[maxbin - 1];
1026  vals[0] = 0.;
1027  }
1028  for (ii = 0; ii < maxbin; ii++) {
1029  vals[ii] = vals[ii] / sum;
1030  IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
1031  }
1032 
1033  // Make IPDFEnergyExist = true
1034  IPDFEnergyExist = true;
1035  if (verbosityLevel > 1)
1036  IPDFEnergyH.DumpValues();
1037  }
1038 
1039  // IPDF has been create so carry on
1040  G4double rndm = eneRndm->GenRandEnergy();
1041  particle_energy = IPDFEnergyH.GetEnergy(rndm);
1042 
1043  if (verbosityLevel >= 1)
1044  G4cout << "Energy is " << particle_energy << G4endl;
1045 }
1046 
1047 void G4SPSEneDistribution::GenArbPointEnergies() {
1048  if (verbosityLevel > 0)
1049  G4cout << "In GenArbPointEnergies" << G4endl;
1050  G4double rndm;
1051  rndm = eneRndm->GenRandEnergy();
1052  // IPDFArbEnergyH.DumpValues();
1053  // Find the Bin
1054  // have x, y, no of points, and cumulative area distribution
1055  G4int nabove, nbelow = 0, middle;
1056  nabove = IPDFArbEnergyH.GetVectorLength();
1057  // G4cout << nabove << G4endl;
1058  // Binary search to find bin that rndm is in
1059  while (nabove - nbelow > 1) {
1060  middle = (nabove + nbelow) / 2;
1061  if (rndm == IPDFArbEnergyH(size_t(middle)))
1062  break;
1063  if (rndm < IPDFArbEnergyH(size_t(middle)))
1064  nabove = middle;
1065  else
1066  nbelow = middle;
1067  }
1068  if (IntType == "Lin") {
1069  Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
1070  Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
1071  grad = Arb_grad[nbelow + 1];
1072  cept = Arb_cept[nbelow + 1];
1073  // G4cout << rndm << " " << Emax << " " << Emin << " " << grad << " " << cept << G4endl;
1074  GenerateLinearEnergies(true);
1075  } else if (IntType == "Log") {
1076  Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
1077  Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
1078  alpha = Arb_alpha[nbelow + 1];
1079  // G4cout << rndm << " " << Emax << " " << Emin << " " << alpha << G4endl;
1080  GeneratePowEnergies(true);
1081  } else if (IntType == "Exp") {
1082  Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
1083  Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
1084  Ezero = Arb_ezero[nbelow + 1];
1085  // G4cout << rndm << " " << Emax << " " << Emin << " " << Ezero << G4endl;
1086  GenerateExpEnergies(true);
1087  } else if (IntType == "Spline") {
1088  Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
1089  Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
1090  particle_energy = -1e100;
1091  rndm = eneRndm->GenRandEnergy();
1092  while (particle_energy < Emin || particle_energy > Emax) {
1093  particle_energy = SplineInt[nbelow+1]->CubicSplineInterpolation(rndm);
1094  rndm = eneRndm->GenRandEnergy();
1095  }
1096  if (verbosityLevel >= 1)
1097  G4cout << "Energy is " << particle_energy << G4endl;
1098  } else
1099  G4cout << "Error: IntType unknown type" << G4endl;
1100 }
1101 
1102 void G4SPSEneDistribution::GenEpnHistEnergies() {
1103  // G4cout << "In GenEpnHistEnergies " << Epnflag << G4endl;
1104 
1105  // Firstly convert to energy if not already done.
1106  if (Epnflag == true)
1107  // epnflag = true means spectrum is epn, false means e.
1108  {
1109  // convert to energy by multiplying by A number
1110  ConvertEPNToEnergy();
1111  // EpnEnergyH will be replace by UDefEnergyH.
1112  // UDefEnergyH.DumpValues();
1113  }
1114 
1115  // G4cout << "Creating IPDFEnergy if not already done so" << G4endl;
1116  if (IPDFEnergyExist == false) {
1117  // IPDF has not been created, so create it
1118  G4double bins[1024], vals[1024], sum;
1119  G4int ii;
1120  G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
1121  bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
1122  vals[0] = UDefEnergyH(size_t(0));
1123  sum = vals[0];
1124  for (ii = 1; ii < maxbin; ii++) {
1125  bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
1126  vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1];
1127  sum = sum + UDefEnergyH(size_t(ii));
1128  }
1129 
1130  for (ii = 0; ii < maxbin; ii++) {
1131  vals[ii] = vals[ii] / sum;
1132  IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
1133  }
1134  // Make IPDFEpnExist = true
1135  IPDFEnergyExist = true;
1136  }
1137  // IPDFEnergyH.DumpValues();
1138  // IPDF has been create so carry on
1139  G4double rndm = eneRndm->GenRandEnergy();
1140  particle_energy = IPDFEnergyH.GetEnergy(rndm);
1141 
1142  if (verbosityLevel >= 1)
1143  G4cout << "Energy is " << particle_energy << G4endl;
1144 }
1145 
1146 void G4SPSEneDistribution::ConvertEPNToEnergy() {
1147  // Use this before particle generation to convert the
1148  // currently stored histogram from energy/nucleon
1149  // to energy.
1150  // G4cout << "In ConvertEpntoEnergy " << G4endl;
1151  if (particle_definition == NULL)
1152  G4cout << "Error: particle not defined" << G4endl;
1153  else {
1154  // Need to multiply histogram by the number of nucleons.
1155  // Baryon Number looks to hold the no. of nucleons.
1156  G4int Bary = particle_definition->GetBaryonNumber();
1157  // G4cout << "Baryon No. " << Bary << G4endl;
1158  // Change values in histogram, Read it out, delete it, re-create it
1159  G4int count, maxcount;
1160  maxcount = G4int(EpnEnergyH.GetVectorLength());
1161  // G4cout << maxcount << G4endl;
1162  G4double ebins[1024], evals[1024];
1163  if (maxcount > 1024) {
1164  G4cout << "Histogram contains more than 1024 bins!" << G4endl;
1165  G4cout << "Those above 1024 will be ignored" << G4endl;
1166  maxcount = 1024;
1167  }
1168  if (maxcount < 1) {
1169  G4cout << "Histogram contains less than 1 bin!" << G4endl;
1170  G4cout << "Redefine the histogram" << G4endl;
1171  return;
1172  }
1173  for (count = 0; count < maxcount; count++) {
1174  // Read out
1175  ebins[count] = EpnEnergyH.GetLowEdgeEnergy(size_t(count));
1176  evals[count] = EpnEnergyH(size_t(count));
1177  }
1178 
1179  // Multiply the channels by the nucleon number to give energies
1180  for (count = 0; count < maxcount; count++) {
1181  ebins[count] = ebins[count] * Bary;
1182  }
1183 
1184  // Set Emin and Emax
1185  Emin = ebins[0];
1186  if (maxcount > 1)
1187  Emax = ebins[maxcount - 1];
1188  else
1189  Emax = ebins[0];
1190  // Put energy bins into new histogram - UDefEnergyH.
1191  for (count = 0; count < maxcount; count++) {
1192  UDefEnergyH.InsertValues(ebins[count], evals[count]);
1193  }
1194  Epnflag = false; //so that you dont repeat this method.
1195  }
1196 }
1197 
1198 //
1200  if (atype == "energy") {
1201  UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
1202  IPDFEnergyExist = false;
1203  Emin = 0.;
1204  Emax = 1e30;
1205  } else if (atype == "arb") {
1206  ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
1207  IPDFArbExist = false;
1208  } else if (atype == "epn") {
1209  UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
1210  IPDFEnergyExist = false;
1211  EpnEnergyH = ZeroPhysVector;
1212  } else {
1213  G4cout << "Error, histtype not accepted " << G4endl;
1214  }
1215 }
1216 
1218  particle_definition = a;
1219  particle_energy = -1.;
1220 
1221  while ((EnergyDisType == "Arb") ? (particle_energy < ArbEmin
1222  || particle_energy > ArbEmax) : (particle_energy < Emin
1223  || particle_energy > Emax)) {
1224  if (Biased) {
1225  GenerateBiasPowEnergies();
1226  } else {
1227  if (EnergyDisType == "Mono")
1228  GenerateMonoEnergetic();
1229  else if (EnergyDisType == "Lin")
1230  GenerateLinearEnergies();
1231  else if (EnergyDisType == "Pow")
1232  GeneratePowEnergies();
1233  else if (EnergyDisType == "Exp")
1234  GenerateExpEnergies();
1235  else if (EnergyDisType == "Gauss")
1236  GenerateGaussEnergies();
1237  else if (EnergyDisType == "Brem")
1238  GenerateBremEnergies();
1239  else if (EnergyDisType == "Bbody")
1240  GenerateBbodyEnergies();
1241  else if (EnergyDisType == "Cdg")
1242  GenerateCdgEnergies();
1243  else if (EnergyDisType == "User")
1244  GenUserHistEnergies();
1245  else if (EnergyDisType == "Arb")
1246  GenArbPointEnergies();
1247  else if (EnergyDisType == "Epn")
1248  GenEpnHistEnergies();
1249  else
1250  G4cout << "Error: EnergyDisType has unusual value" << G4endl;
1251  }
1252  }
1253  return particle_energy;
1254 }
1255 
1257  G4double prob = 1.;
1258 
1259  if (EnergyDisType == "Lin") {
1260  if (prob_norm == 1.) {
1261  prob_norm = 0.5*grad*Emax*Emax + cept*Emax - 0.5*grad*Emin*Emin - cept*Emin;
1262  }
1263  prob = cept + grad * ene;
1264  prob /= prob_norm;
1265  }
1266  else if (EnergyDisType == "Pow") {
1267  if (prob_norm == 1.) {
1268  if (alpha != -1.) {
1269  G4double emina = std::pow(Emin, alpha + 1);
1270  G4double emaxa = std::pow(Emax, alpha + 1);
1271  prob_norm = 1./(1.+alpha) * (emaxa - emina);
1272  } else {
1273  prob_norm = std::log(Emax) - std::log(Emin) ;
1274  }
1275  }
1276  prob = std::pow(ene, alpha)/prob_norm;
1277  }
1278  else if (EnergyDisType == "Exp"){
1279  if (prob_norm == 1.) {
1280  prob_norm = -Ezero*(std::exp(-Emax/Ezero) - std::exp(Emin/Ezero));
1281  }
1282  prob = std::exp(-ene / Ezero);
1283  prob /= prob_norm;
1284  }
1285  else if (EnergyDisType == "Arb") {
1286  prob = ArbEnergyH.Value(ene);
1287  // prob = ArbEInt->CubicSplineInterpolation(ene);
1288  //G4double deltaY;
1289  //prob = ArbEInt->PolynomInterpolation(ene, deltaY);
1290  if (prob <= 0.) {
1291  //G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << " " <<deltaY<< G4endl;
1292  G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << G4endl;
1293  prob = 1e-30;
1294  }
1295  // already normalised
1296  }
1297  else
1298  G4cout << "Error: EnergyDisType not supported" << G4endl;
1299 
1300  return prob;
1301 }