| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4SPSEneDistribution class implementation << 26 ///////////////////////////////////////////////////////////////////////////////
27 // 27 //
28 // Author: Fan Lei, QinetiQ ltd - 05/02/2004 << 28 // MODULE: G4SPSEneDistribution.cc
29 // Customer: ESA/ESTEC << 29 //
30 // Revisions: Andrew Green, Andrea Dotti << 30 // Version: 1.0
31 // ------------------------------------------- << 31 // Date: 5/02/04
>> 32 // Author: Fan Lei
>> 33 // Organisation: QinetiQ ltd.
>> 34 // Customer: ESA/ESTEC
>> 35 //
>> 36 ///////////////////////////////////////////////////////////////////////////////
>> 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 //
>> 45 ///////////////////////////////////////////////////////////////////////////////
>> 46 //
>> 47
32 #include "G4SPSEneDistribution.hh" 48 #include "G4SPSEneDistribution.hh"
33 49
34 #include "G4Exp.hh" <<
35 #include "G4SystemOfUnits.hh" 50 #include "G4SystemOfUnits.hh"
36 #include "G4UnitsTable.hh" <<
37 #include "Randomize.hh" 51 #include "Randomize.hh"
38 #include "G4AutoLock.hh" <<
39 #include "G4Threading.hh" <<
40 52
41 G4SPSEneDistribution::G4SPSEneDistribution() 53 G4SPSEneDistribution::G4SPSEneDistribution()
>> 54 : particle_definition(0), eneRndm(0), Splinetemp(0)
42 { 55 {
43 G4MUTEXINIT(mutex); << 56 //
44 << 57 // Initialise all variables
45 // Initialise all variables << 58 particle_energy = 1.0 * MeV;
46 << 59
47 particle_energy = 1.0 * MeV; << 60 EnergyDisType = "Mono";
48 EnergyDisType = "Mono"; << 61 weight = 1.;
49 weight=1.; << 62 MonoEnergy = 1 * MeV;
50 MonoEnergy = 1 * MeV; << 63 Emin = 0.;
51 Emin = 0.; << 64 Emax = 1.e30;
52 Emax = 1.e30; << 65 alpha = 0.;
53 alpha = 0.; << 66 biasalpha = 0.;
54 biasalpha = 0.; << 67 prob_norm = 1.0;
55 prob_norm = 1.0; << 68 Ezero = 0.;
56 Ezero = 0.; << 69 SE = 0.;
57 SE = 0.; << 70 Temp = 0.;
58 Temp = 0.; << 71 grad = 0.;
59 grad = 0.; << 72 cept = 0.;
60 cept = 0.; << 73 Biased = false; // not biased
61 IntType = "NULL"; // Interpolation type << 74 EnergySpec = true; // true - energy spectra, false - momentum spectra
62 << 75 DiffSpec = true; // true - differential spec, false integral spec
63 ArbEmin = 0.; << 76 IntType = "NULL"; // Interpolation type
64 ArbEmax = 1.e30; << 77 IPDFEnergyExist = false;
65 << 78 IPDFArbExist = false;
66 verbosityLevel = 0; << 79
67 << 80 ArbEmin = 0.;
68 threadLocal_t& data = threadLocalData.Get(); << 81 ArbEmax = 1.e30;
69 data.Emax = Emax; << 82
70 data.Emin = Emin; << 83 verbosityLevel = 0;
71 data.alpha =alpha; << 84
72 data.cept = cept; << 85 }
73 data.Ezero = Ezero; << 86
74 data.grad = grad; << 87 G4SPSEneDistribution::~G4SPSEneDistribution() {
75 data.particle_energy = 0.; << 88 }
76 data.particle_definition = nullptr; << 89
77 data.weight = weight; << 90 void G4SPSEneDistribution::SetEnergyDisType(G4String DisType) {
78 } << 91 EnergyDisType = DisType;
79 << 92 if (EnergyDisType == "User") {
80 G4SPSEneDistribution::~G4SPSEneDistribution() << 93 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
81 { << 94 IPDFEnergyExist = false;
82 G4MUTEXDESTROY(mutex); << 95 } else if (EnergyDisType == "Arb") {
83 if(Arb_grad_cept_flag) << 96 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
84 { << 97 IPDFArbExist = false;
85 delete [] Arb_grad; << 98 } else if (EnergyDisType == "Epn") {
86 delete [] Arb_cept; << 99 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
87 } << 100 IPDFEnergyExist = false;
88 << 101 EpnEnergyH = ZeroPhysVector;
89 if(Arb_alpha_Const_flag) << 102 }
90 { << 103 }
91 delete [] Arb_alpha; << 104
92 delete [] Arb_Const; << 105 void G4SPSEneDistribution::SetEmin(G4double emi) {
93 } << 106 Emin = emi;
94 << 107 }
95 if(Arb_ezero_flag) << 108
96 { << 109 void G4SPSEneDistribution::SetEmax(G4double ema) {
97 delete [] Arb_ezero; << 110 Emax = ema;
98 } << 111 }
99 delete Bbody_x; << 112
100 delete BBHist; << 113 void G4SPSEneDistribution::SetMonoEnergy(G4double menergy) {
101 delete CP_x; << 114 MonoEnergy = menergy;
102 delete CPHist; << 115 }
103 for (auto & it : SplineInt) << 116
104 { << 117 void G4SPSEneDistribution::SetBeamSigmaInE(G4double e) {
105 delete it; << 118 SE = e;
106 it = nullptr; << 119 }
107 } << 120 void G4SPSEneDistribution::SetAlpha(G4double alp) {
108 SplineInt.clear(); << 121 alpha = alp;
109 } << 122 }
110 << 123
111 void G4SPSEneDistribution::SetEnergyDisType(co << 124 void G4SPSEneDistribution::SetBiasAlpha(G4double alp) {
112 { << 125 biasalpha = alp;
113 G4AutoLock l(&mutex); << 126 Biased = true;
114 EnergyDisType = DisType; << 127 }
115 if (EnergyDisType == "User") << 128
116 { << 129 void G4SPSEneDistribution::SetTemp(G4double tem) {
117 UDefEnergyH = IPDFEnergyH = ZeroPhysVector << 130 Temp = tem;
118 IPDFEnergyExist = false; << 131 }
119 } << 132
120 else if (EnergyDisType == "Arb") << 133 void G4SPSEneDistribution::SetEzero(G4double eze) {
121 { << 134 Ezero = eze;
122 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVect << 135 }
123 IPDFArbExist = false; << 136
124 } << 137 void G4SPSEneDistribution::SetGradient(G4double gr) {
125 else if (EnergyDisType == "Epn") << 138 grad = gr;
126 { << 139 }
127 UDefEnergyH = IPDFEnergyH = ZeroPhysVector << 140
128 IPDFEnergyExist = false; << 141 void G4SPSEneDistribution::SetInterCept(G4double c) {
129 EpnEnergyH = ZeroPhysVector; << 142 cept = c;
130 } << 143 }
131 } << 144
132 << 145 void G4SPSEneDistribution::UserEnergyHisto(G4ThreeVector input) {
133 const G4String& G4SPSEneDistribution::GetEnerg << 146 G4double ehi, val;
134 { << 147 ehi = input.x();
135 G4AutoLock l(&mutex); << 148 val = input.y();
136 return EnergyDisType; << 149 if (verbosityLevel > 1) {
137 } << 150 G4cout << "In UserEnergyHisto" << G4endl;
138 << 151 G4cout << " " << ehi << " " << val << G4endl;
139 void G4SPSEneDistribution::SetEmin(G4double em << 152 }
140 { << 153 UDefEnergyH.InsertValues(ehi, val);
141 G4AutoLock l(&mutex); << 154 Emax = ehi;
142 Emin = emi; << 155 }
143 threadLocalData.Get().Emin = Emin; << 156
144 } << 157 void G4SPSEneDistribution::ArbEnergyHisto(G4ThreeVector input) {
145 << 158 G4double ehi, val;
146 G4double G4SPSEneDistribution::GetEmin() const << 159 ehi = input.x();
147 { << 160 val = input.y();
148 return threadLocalData.Get().Emin; << 161 if (verbosityLevel > 1) {
149 } << 162 G4cout << "In ArbEnergyHisto" << G4endl;
150 << 163 G4cout << " " << ehi << " " << val << G4endl;
151 G4double G4SPSEneDistribution::GetArbEmin() << 164 }
152 { << 165 ArbEnergyH.InsertValues(ehi, val);
153 G4AutoLock l(&mutex); << 166 }
154 return ArbEmin; << 167
155 } << 168 void G4SPSEneDistribution::ArbEnergyHistoFile(G4String filename) {
156 << 169 std::ifstream infile(filename, std::ios::in);
157 G4double G4SPSEneDistribution::GetArbEmax() << 170 if (!infile)
158 { << 171 G4Exception("G4SPSEneDistribution::ArbEnergyHistoFile",
159 G4AutoLock l(&mutex); << 172 "Event0301",FatalException,
160 return ArbEmax; << 173 "Unable to open the histo ASCII file");
161 } << 174 G4double ehi, val;
162 << 175 while (infile >> ehi >> val) {
163 void G4SPSEneDistribution::SetEmax(G4double em << 176 ArbEnergyH.InsertValues(ehi, val);
164 { << 177 }
165 G4AutoLock l(&mutex); << 178 }
166 Emax = ema; << 179
167 threadLocalData.Get().Emax = Emax; << 180 void G4SPSEneDistribution::EpnEnergyHisto(G4ThreeVector input) {
168 } << 181 G4double ehi, val;
169 << 182 ehi = input.x();
170 G4double G4SPSEneDistribution::GetEmax() const << 183 val = input.y();
171 { << 184 if (verbosityLevel > 1) {
172 return threadLocalData.Get().Emax; << 185 G4cout << "In EpnEnergyHisto" << G4endl;
173 } << 186 G4cout << " " << ehi << " " << val << G4endl;
174 << 187 }
175 void G4SPSEneDistribution::SetMonoEnergy(G4dou << 188 EpnEnergyH.InsertValues(ehi, val);
176 { << 189 Emax = ehi;
177 G4AutoLock l(&mutex); << 190 Epnflag = true;
178 MonoEnergy = menergy; << 191 }
179 } << 192
180 << 193 void G4SPSEneDistribution::Calculate() {
181 void G4SPSEneDistribution::SetBeamSigmaInE(G4d << 194 if (EnergyDisType == "Cdg")
182 { << 195 CalculateCdgSpectrum();
183 G4AutoLock l(&mutex); << 196 else if (EnergyDisType == "Bbody")
184 SE = e; << 197 CalculateBbodySpectrum();
185 } << 198 }
186 void G4SPSEneDistribution::SetAlpha(G4double a << 199
187 { << 200 void G4SPSEneDistribution::CalculateCdgSpectrum() {
188 G4AutoLock l(&mutex); << 201 // This uses the spectrum from The INTEGRAL Mass Model (TIMM)
189 alpha = alp; << 202 // to generate a Cosmic Diffuse X/gamma ray spectrum.
190 threadLocalData.Get().alpha = alpha; << 203 G4double pfact[2] = { 8.5, 112 };
191 } << 204 G4double spind[2] = { 1.4, 2.3 };
192 << 205 G4double ene_line[3] = { 1. * keV, 18. * keV, 1E6 * keV };
193 void G4SPSEneDistribution::SetBiasAlpha(G4doub << 206 G4int n_par;
194 { << 207
195 G4AutoLock l(&mutex); << 208 ene_line[0] = Emin;
196 biasalpha = alp; << 209 if (Emin < 18 * keV) {
197 Biased = true; << 210 n_par = 2;
198 } << 211 ene_line[2] = Emax;
199 << 212 if (Emax < 18 * keV) {
200 void G4SPSEneDistribution::SetTemp(G4double te << 213 n_par = 1;
201 { << 214 ene_line[1] = Emax;
202 G4AutoLock l(&mutex); << 215 }
203 Temp = tem; << 216 } else {
204 } << 217 n_par = 1;
205 << 218 pfact[0] = 112.;
206 void G4SPSEneDistribution::SetEzero(G4double e << 219 spind[0] = 2.3;
207 { << 220 ene_line[1] = Emax;
208 G4AutoLock l(&mutex); << 221 }
209 Ezero = eze; << 222
210 threadLocalData.Get().Ezero = Ezero; << 223 // Create a cumulative histogram.
211 } << 224 CDGhist[0] = 0.;
212 << 225 G4double omalpha;
213 void G4SPSEneDistribution::SetGradient(G4doubl << 226 G4int i = 0;
214 { << 227
215 G4AutoLock l(&mutex); << 228 while (i < n_par) {
216 grad = gr; << 229 omalpha = 1. - spind[i];
217 threadLocalData.Get().grad = grad; << 230 CDGhist[i + 1] = CDGhist[i] + (pfact[i] / omalpha) * (std::pow(
218 } << 231 ene_line[i + 1] / keV, omalpha) - std::pow(ene_line[i] / keV,
219 << 232 omalpha));
220 void G4SPSEneDistribution::SetInterCept(G4doub << 233 i++;
221 { << 234 }
222 G4AutoLock l(&mutex); << 235
223 cept = c; << 236 // Normalise histo and divide by 1000 to make MeV.
224 threadLocalData.Get().cept = cept; << 237 i = 0;
225 } << 238 while (i < n_par) {
226 << 239 CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[n_par];
227 const G4String& G4SPSEneDistribution::GetIntTy << 240 // G4cout << CDGhist[i] << CDGhist[n_par] << G4endl;
228 { << 241 i++;
229 G4AutoLock l(&mutex); << 242 }
230 return IntType; << 243 }
231 } << 244
232 << 245 void G4SPSEneDistribution::CalculateBbodySpectrum() {
233 void G4SPSEneDistribution::SetBiasRndm(G4SPSRa << 246 // create bbody spectrum
234 { << 247 // Proved very hard to integrate indefinitely, so different
235 G4AutoLock l(&mutex); << 248 // method. User inputs emin, emax and T. These are used to
236 eneRndm = a; << 249 // create a 10,000 bin histogram.
237 } << 250 // Use photon density spectrum = 2 nu**2/c**2 * (std::exp(h nu/kT)-1)
238 << 251 // = 2 E**2/h**2c**2 times the exponential
239 void G4SPSEneDistribution::SetVerbosity(G4int << 252 G4double erange = Emax - Emin;
240 { << 253 G4double steps = erange / 10000.;
241 G4AutoLock l(&mutex); << 254 G4double Bbody_y[10000];
242 verbosityLevel = a; << 255 G4double k = 8.6181e-11; //Boltzmann const in MeV/K
243 } << 256 G4double h = 4.1362e-21; // Plancks const in MeV s
244 << 257 G4double c = 3e8; // Speed of light
245 G4double G4SPSEneDistribution::GetWeight() con << 258 G4double h2 = h * h;
246 { << 259 G4double c2 = c * c;
247 return threadLocalData.Get().weight; << 260 G4int count = 0;
248 } << 261 G4double sum = 0.;
249 << 262 BBHist[0] = 0.;
250 G4double G4SPSEneDistribution::GetMonoEnergy() << 263 while (count < 10000) {
251 { << 264 Bbody_x[count] = Emin + G4double(count * steps);
252 G4AutoLock l(&mutex); << 265 Bbody_y[count] = (2. * std::pow(Bbody_x[count], 2.)) / (h2 * c2
253 return MonoEnergy; << 266 * (std::exp(Bbody_x[count] / (k * Temp)) - 1.));
254 } << 267 sum = sum + Bbody_y[count];
255 << 268 BBHist[count + 1] = BBHist[count] + Bbody_y[count];
256 G4double G4SPSEneDistribution::GetSE() << 269 count++;
257 { << 270 }
258 G4AutoLock l(&mutex); << 271
259 return SE; << 272 Bbody_x[10000] = Emax;
260 } << 273 // Normalise cumulative histo.
261 << 274 count = 0;
262 G4double G4SPSEneDistribution::Getalpha() cons << 275 while (count < 10001) {
263 { << 276 BBHist[count] = BBHist[count] / sum;
264 return threadLocalData.Get().alpha; << 277 count++;
265 } << 278 }
266 << 279 }
267 G4double G4SPSEneDistribution::GetEzero() cons << 280
268 { << 281 void G4SPSEneDistribution::InputEnergySpectra(G4bool value) {
269 return threadLocalData.Get().Ezero; << 282 // Allows user to specifiy spectrum is momentum
270 } << 283 EnergySpec = value; // false if momentum
271 << 284 if (verbosityLevel > 1)
272 G4double G4SPSEneDistribution::GetTemp() << 285 G4cout << "EnergySpec has value " << EnergySpec << G4endl;
273 { << 286 }
274 G4AutoLock l(&mutex); << 287
275 return Temp; << 288 void G4SPSEneDistribution::InputDifferentialSpectra(G4bool value) {
276 } << 289 // Allows user to specify integral or differential spectra
277 << 290 DiffSpec = value; // true = differential, false = integral
278 G4double G4SPSEneDistribution::Getgrad() const << 291 if (verbosityLevel > 1)
279 { << 292 G4cout << "Diffspec has value " << DiffSpec << G4endl;
280 return threadLocalData.Get().grad; << 293 }
281 } << 294
282 << 295 void G4SPSEneDistribution::ArbInterpolate(G4String IType) {
283 G4double G4SPSEneDistribution::Getcept() const << 296 if (EnergyDisType != "Arb")
284 { << 297 G4cout << "Error: this is for arbitrary distributions" << G4endl;
285 return threadLocalData.Get().cept; << 298 IntType = IType;
286 } << 299 ArbEmax = ArbEnergyH.GetMaxLowEdgeEnergy();
287 << 300 ArbEmin = ArbEnergyH.GetMinLowEdgeEnergy();
288 G4PhysicsFreeVector G4SPSEneDistribution::GetU << 301
289 { << 302 // Now interpolate points
290 G4AutoLock l(&mutex); << 303 if (IntType == "Lin")
291 return UDefEnergyH; << 304 LinearInterpolation();
292 } << 305 if (IntType == "Log")
293 << 306 LogInterpolation();
294 G4PhysicsFreeVector G4SPSEneDistribution::GetA << 307 if (IntType == "Exp")
295 { << 308 ExpInterpolation();
296 G4AutoLock l(&mutex); << 309 if (IntType == "Spline")
297 return ArbEnergyH; << 310 SplineInterpolation();
298 } << 311 }
299 << 312
300 void G4SPSEneDistribution::UserEnergyHisto(con << 313 void G4SPSEneDistribution::LinearInterpolation() {
301 { << 314 // Method to do linear interpolation on the Arb points
302 G4AutoLock l(&mutex); << 315 // Calculate equation of each line segment, max 1024.
303 G4double ehi = input.x(), << 316 // Calculate Area under each segment
304 val = input.y(); << 317 // Create a cumulative array which is then normalised Arb_Cum_Area
305 if (verbosityLevel > 1) << 318
306 { << 319 G4double Area_seg[1024]; // Stores area under each segment
307 G4cout << "In UserEnergyHisto" << G4endl; << 320 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
308 G4cout << " " << ehi << " " << val << G4en << 321 G4int i, count;
309 } << 322 G4int maxi = ArbEnergyH.GetVectorLength();
310 UDefEnergyH.InsertValues(ehi, val); << 323 for (i = 0; i < maxi; i++) {
311 Emax = ehi; << 324 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
312 threadLocalData.Get().Emax = Emax; << 325 Arb_y[i] = ArbEnergyH(size_t(i));
313 } << 326 }
314 << 327 // Points are now in x,y arrays. If the spectrum is integral it has to be
315 void G4SPSEneDistribution::ArbEnergyHisto(cons << 328 // made differential and if momentum it has to be made energy.
316 { << 329 if (DiffSpec == false) {
317 G4AutoLock l(&mutex); << 330 // Converts integral point-wise spectra to Differential
318 G4double ehi = input.x(), << 331 for (count = 0; count < maxi - 1; count++) {
319 val = input.y(); << 332 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
320 if (verbosityLevel > 1) << 333 / (Arb_x[count + 1] - Arb_x[count]);
321 { << 334 }
322 G4cout << "In ArbEnergyHisto" << G4endl; << 335 maxi--;
323 G4cout << " " << ehi << " " << val << G4en << 336 }
324 } << 337 //
325 ArbEnergyH.InsertValues(ehi, val); << 338 if (EnergySpec == false) {
326 } << 339 // change currently stored values (emin etc) which are actually momenta
327 << 340 // to energies.
328 void G4SPSEneDistribution::ArbEnergyHistoFile( << 341 if (particle_definition == NULL)
329 { << 342 G4cout << "Error: particle not defined" << G4endl;
330 G4AutoLock l(&mutex); << 343 else {
331 std::ifstream infile(filename, std::ios::in) << 344 // Apply Energy**2 = p**2c**2 + m0**2c**4
332 if (!infile) << 345 // p should be entered as E/c i.e. without the division by c
333 { << 346 // being done - energy equivalent.
334 G4Exception("G4SPSEneDistribution::ArbEner << 347 G4double mass = particle_definition->GetPDGMass();
335 FatalException, "Unable to ope << 348 // convert point to energy unit and its value to per energy unit
336 } << 349 G4double total_energy;
337 G4double ehi, val; << 350 for (count = 0; count < maxi; count++) {
338 while (infile >> ehi >> val) << 351 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
339 { << 352 * mass)); // total energy
340 ArbEnergyH.InsertValues(ehi, val); << 353
341 } << 354 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
342 } << 355 Arb_x[count] = total_energy - mass; // kinetic energy
343 << 356 }
344 void G4SPSEneDistribution::EpnEnergyHisto(cons << 357 }
345 { << 358 }
346 G4AutoLock l(&mutex); << 359 //
347 G4double ehi = input.x(), << 360 i = 1;
348 val = input.y(); << 361 Arb_grad[0] = 0.;
349 if (verbosityLevel > 1) << 362 Arb_cept[0] = 0.;
350 { << 363 Area_seg[0] = 0.;
351 G4cout << "In EpnEnergyHisto" << G4endl; << 364 Arb_Cum_Area[0] = 0.;
352 G4cout << " " << ehi << " " << val << G4en << 365 while (i < maxi) {
353 } << 366 // calc gradient and intercept for each segment
354 EpnEnergyH.InsertValues(ehi, val); << 367 Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[i - 1]);
355 Emax = ehi; << 368 if (verbosityLevel == 2)
356 threadLocalData.Get().Emax = Emax; << 369 G4cout << Arb_grad[i] << G4endl;
357 Epnflag = true; << 370 if (Arb_grad[i] > 0.) {
358 } << 371 if (verbosityLevel == 2)
359 << 372 G4cout << "Arb_grad is positive" << G4endl;
360 void G4SPSEneDistribution::Calculate() << 373 Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]);
361 { << 374 } else if (Arb_grad[i] < 0.) {
362 G4AutoLock l(&mutex); << 375 if (verbosityLevel == 2)
363 if (EnergyDisType == "Cdg") << 376 G4cout << "Arb_grad is negative" << G4endl;
364 { << 377 Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]);
365 CalculateCdgSpectrum(); << 378 } else {
366 } << 379 if (verbosityLevel == 2)
367 else if (EnergyDisType == "Bbody") << 380 G4cout << "Arb_grad is 0." << G4endl;
368 { << 381 Arb_cept[i] = Arb_y[i];
369 if(!BBhistInit) << 382 }
370 { << 383
371 BBInitHists(); << 384 Area_seg[i] = ((Arb_grad[i] / 2) * (Arb_x[i] * Arb_x[i] - Arb_x[i - 1]
372 } << 385 * Arb_x[i - 1]) + Arb_cept[i] * (Arb_x[i] - Arb_x[i - 1]));
373 CalculateBbodySpectrum(); << 386 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
374 } << 387 sum = sum + Area_seg[i];
375 else if (EnergyDisType == "CPow") << 388 if (verbosityLevel == 2)
376 { << 389 G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i]
377 if(!CPhistInit) << 390 << G4endl;
378 { << 391 i++;
379 CPInitHists(); << 392 }
380 } << 393
381 CalculateCPowSpectrum(); << 394 i = 0;
382 } << 395 while (i < maxi) {
383 } << 396 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation
384 << 397 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
385 void G4SPSEneDistribution::BBInitHists() // M << 398 i++;
386 { << 399 }
387 BBHist = new std::vector<G4double>(10001, 0. << 400
388 Bbody_x = new std::vector<G4double>(10001, 0 << 401 // now scale the ArbEnergyH, needed by Probability()
389 BBhistInit = true; << 402 ArbEnergyH.ScaleVector(1., 1./sum);
390 } << 403
391 << 404 if (verbosityLevel >= 1) {
392 void G4SPSEneDistribution::CPInitHists() // M << 405 G4cout << "Leaving LinearInterpolation" << G4endl;
393 { << 406 ArbEnergyH.DumpValues();
394 CPHist = new std::vector<G4double>(10001, 0. << 407 IPDFArbEnergyH.DumpValues();
395 CP_x = new std::vector<G4double>(10001, 0.0) << 408 }
396 CPhistInit = true; << 409 }
397 } << 410
398 << 411 void G4SPSEneDistribution::LogInterpolation() {
399 void G4SPSEneDistribution::CalculateCdgSpectru << 412 // Interpolation based on Logarithmic equations
400 { << 413 // Generate equations of line segments
401 // This uses the spectrum from the INTEGRAL << 414 // y = Ax**alpha => log y = alpha*logx + logA
402 // to generate a Cosmic Diffuse X/gamma ray << 415 // Find area under line segments
403 << 416 // create normalised, cumulative array Arb_Cum_Area
404 G4double pfact[2] = { 8.5, 112 }; << 417 G4double Area_seg[1024]; // Stores area under each segment
405 G4double spind[2] = { 1.4, 2.3 }; << 418 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
406 G4double ene_line[3] = { 1. * keV, 18. * keV << 419 G4int i, count;
407 G4int n_par; << 420 G4int maxi = ArbEnergyH.GetVectorLength();
408 << 421 for (i = 0; i < maxi; i++) {
409 ene_line[0] = threadLocalData.Get().Emin; << 422 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
410 if (threadLocalData.Get().Emin < 18 * keV) << 423 Arb_y[i] = ArbEnergyH(size_t(i));
411 { << 424 }
412 n_par = 2; << 425 // Points are now in x,y arrays. If the spectrum is integral it has to be
413 ene_line[2] = threadLocalData.Get().Emax; << 426 // made differential and if momentum it has to be made energy.
414 if (threadLocalData.Get().Emax < 18 * keV) << 427 if (DiffSpec == false) {
415 { << 428 // Converts integral point-wise spectra to Differential
416 n_par = 1; << 429 for (count = 0; count < maxi - 1; count++) {
417 ene_line[1] = threadLocalData.Get().Emax << 430 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
418 } << 431 / (Arb_x[count + 1] - Arb_x[count]);
419 } << 432 }
420 else << 433 maxi--;
421 { << 434 }
422 n_par = 1; << 435 //
423 pfact[0] = 112.; << 436 if (EnergySpec == false) {
424 spind[0] = 2.3; << 437 // change currently stored values (emin etc) which are actually momenta
425 ene_line[1] = threadLocalData.Get().Emax; << 438 // to energies.
426 } << 439 if (particle_definition == NULL)
427 << 440 G4cout << "Error: particle not defined" << G4endl;
428 // Create a cumulative histogram << 441 else {
429 // << 442 // Apply Energy**2 = p**2c**2 + m0**2c**4
430 CDGhist[0] = 0.; << 443 // p should be entered as E/c i.e. without the division by c
431 G4double omalpha; << 444 // being done - energy equivalent.
432 G4int i = 0; << 445 G4double mass = particle_definition->GetPDGMass();
433 while (i < n_par) << 446 // convert point to energy unit and its value to per energy unit
434 { << 447 G4double total_energy;
435 omalpha = 1. - spind[i]; << 448 for (count = 0; count < maxi; count++) {
436 CDGhist[i + 1] = CDGhist[i] + (pfact[i] / << 449 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
437 * (std::pow(en << 450 * mass)); // total energy
438 - std::pow(ene << 451
439 ++i; << 452 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
440 } << 453 Arb_x[count] = total_energy - mass; // kinetic energy
441 << 454 }
442 // Normalise histo and divide by 1000 to mak << 455 }
443 // << 456 }
444 i = 0; << 457 //
445 while (i < n_par) << 458 i = 1;
446 { << 459 Arb_alpha[0] = 0.;
447 CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[ << 460 Arb_Const[0] = 0.;
448 ++i; << 461 Area_seg[0] = 0.;
449 } << 462 Arb_Cum_Area[0] = 0.;
450 } << 463 if (Arb_x[0] <= 0. || Arb_y[0] <= 0.) {
451 << 464 G4cout << "You should not use log interpolation with points <= 0."
452 void G4SPSEneDistribution::CalculateBbodySpect << 465 << G4endl;
453 { << 466 G4cout << "These will be changed to 1e-20, which may cause problems"
454 // Create bbody spectrum << 467 << G4endl;
455 // Proved very hard to integrate indefinitel << 468 if (Arb_x[0] <= 0.)
456 // User inputs emin, emax and T. These are u << 469 Arb_x[0] = 1e-20;
457 // bin histogram. << 470 if (Arb_y[0] <= 0.)
458 // Use photon density spectrum = 2 nu**2/c** << 471 Arb_y[0] = 1e-20;
459 // = 2 E**2/h**2c**2 times the exponential << 472 }
460 << 473
461 G4double erange = threadLocalData.Get().Emax << 474 G4double alp;
462 G4double steps = erange / 10000.; << 475 while (i < maxi) {
463 << 476 // Incase points are negative or zero
464 const G4double k = 8.6181e-11; //Boltzmann c << 477 if (Arb_x[i] <= 0. || Arb_y[i] <= 0.) {
465 const G4double h = 4.1362e-21; // Plancks co << 478 G4cout << "You should not use log interpolation with points <= 0."
466 const G4double c = 3e8; // Speed of light << 479 << G4endl;
467 const G4double h2 = h * h; << 480 G4cout
468 const G4double c2 = c * c; << 481 << "These will be changed to 1e-20, which may cause problems"
469 G4int count = 0; << 482 << G4endl;
470 G4double sum = 0.; << 483 if (Arb_x[i] <= 0.)
471 BBHist->at(0) = 0.; << 484 Arb_x[i] = 1e-20;
472 << 485 if (Arb_y[i] <= 0.)
473 while (count < 10000) << 486 Arb_y[i] = 1e-20;
474 { << 487 }
475 Bbody_x->at(count) = threadLocalData.Get() << 488
476 G4double Bbody_y = (2. * std::pow(Bbody_x- << 489 Arb_alpha[i] = (std::log10(Arb_y[i]) - std::log10(Arb_y[i - 1]))
477 / (h2*c2*(std::exp(Bbody_ << 490 / (std::log10(Arb_x[i]) - std::log10(Arb_x[i - 1]));
478 sum = sum + Bbody_y; << 491 Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[i], Arb_alpha[i]));
479 BBHist->at(count + 1) = BBHist->at(count) << 492 alp = Arb_alpha[i] + 1;
480 ++count; << 493 if (alp == 0.) {
481 } << 494 Area_seg[i] = Arb_Const[i] * (std::log(Arb_x[i]) - std::log(Arb_x[i - 1]));
482 << 495 } else {
483 Bbody_x->at(10000) = threadLocalData.Get().E << 496 Area_seg[i] = (Arb_Const[i] / alp) * (std::pow(Arb_x[i], alp)
484 << 497 - std::pow(Arb_x[i - 1], alp));
485 // Normalise cumulative histo << 498 }
486 // << 499 sum = sum + Area_seg[i];
487 count = 0; << 500 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
488 while (count < 10001) << 501 if (verbosityLevel == 2)
489 { << 502 G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl;
490 BBHist->at(count) = BBHist->at(count) / su << 503 i++;
491 ++count; << 504 }
492 } << 505
493 } << 506 i = 0;
494 << 507 while (i < maxi) {
495 void G4SPSEneDistribution::CalculateCPowSpectr << 508 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
496 { << 509 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
497 // Create cutoff power-law spectrum, x^a exp << 510 i++;
498 // The integral of this function is an incom << 511 }
499 // is only available in the Boost library. << 512
500 // << 513 // now scale the ArbEnergyH, needed by Probability()
501 // User inputs are emin, emax and alpha and << 514 ArbEnergyH.ScaleVector(1., 1./sum);
502 // create a 10,000 bin histogram. << 515
503 << 516 if (verbosityLevel >= 1)
504 G4double erange = threadLocalData.Get().Emax << 517 G4cout << "Leaving LogInterpolation " << G4endl;
505 G4double steps = erange / 10000.; << 518 }
506 alpha = threadLocalData.Get().alpha ; << 519
507 Ezero = threadLocalData.Get().Ezero ; << 520 void G4SPSEneDistribution::ExpInterpolation() {
508 << 521 // Interpolation based on Exponential equations
509 G4int count = 0; << 522 // Generate equations of line segments
510 G4double sum = 0.; << 523 // y = Ae**-(x/e0) => ln y = -x/e0 + lnA
511 CPHist->at(0) = 0.; << 524 // Find area under line segments
512 << 525 // create normalised, cumulative array Arb_Cum_Area
513 while (count < 10000) << 526 G4double Area_seg[1024]; // Stores area under each segment
514 { << 527 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
515 CP_x->at(count) = threadLocalData.Get().Em << 528 G4int i, count;
516 G4double CP_y = std::pow(CP_x->at(count), << 529 G4int maxi = ArbEnergyH.GetVectorLength();
517 * std::exp(-CP_x->at(count) << 530 for (i = 0; i < maxi; i++) {
518 sum = sum + CP_y; << 531 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
519 CPHist->at(count + 1) = CPHist->at(count) << 532 Arb_y[i] = ArbEnergyH(size_t(i));
520 ++count; << 533 }
521 } << 534 // Points are now in x,y arrays. If the spectrum is integral it has to be
522 << 535 // made differential and if momentum it has to be made energy.
523 CP_x->at(10000) = threadLocalData.Get().Emax << 536 if (DiffSpec == false) {
524 << 537 // Converts integral point-wise spectra to Differential
525 // Normalise cumulative histo << 538 for (count = 0; count < maxi - 1; count++) {
526 // << 539 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
527 count = 0; << 540 / (Arb_x[count + 1] - Arb_x[count]);
528 while (count < 10001) << 541 }
529 { << 542 maxi--;
530 CPHist->at(count) = CPHist->at(count) / su << 543 }
531 ++count; << 544 //
532 } << 545 if (EnergySpec == false) {
533 } << 546 // change currently stored values (emin etc) which are actually momenta
534 << 547 // to energies.
535 void G4SPSEneDistribution::InputEnergySpectra( << 548 if (particle_definition == NULL)
536 { << 549 G4cout << "Error: particle not defined" << G4endl;
537 G4AutoLock l(&mutex); << 550 else {
538 << 551 // Apply Energy**2 = p**2c**2 + m0**2c**4
539 // Allows user to specify spectrum is moment << 552 // p should be entered as E/c i.e. without the division by c
540 // << 553 // being done - energy equivalent.
541 EnergySpec = value; // false if momentum << 554 G4double mass = particle_definition->GetPDGMass();
542 if (verbosityLevel > 1) << 555 // convert point to energy unit and its value to per energy unit
543 { << 556 G4double total_energy;
544 G4cout << "EnergySpec has value " << Energ << 557 for (count = 0; count < maxi; count++) {
545 } << 558 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
546 } << 559 * mass)); // total energy
547 << 560
548 void G4SPSEneDistribution::InputDifferentialSp << 561 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
549 { << 562 Arb_x[count] = total_energy - mass; // kinetic energy
550 G4AutoLock l(&mutex); << 563 }
551 << 564 }
552 // Allows user to specify integral or differ << 565 }
553 // << 566 //
554 DiffSpec = value; // true = differential, fa << 567 i = 1;
555 if (verbosityLevel > 1) << 568 Arb_ezero[0] = 0.;
556 { << 569 Arb_Const[0] = 0.;
557 G4cout << "Diffspec has value " << DiffSpe << 570 Area_seg[0] = 0.;
558 } << 571 Arb_Cum_Area[0] = 0.;
559 } << 572 while (i < maxi) {
560 << 573 G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i - 1]);
561 void G4SPSEneDistribution::ArbInterpolate(cons << 574 if (test > 0. || test < 0.) {
562 { << 575 Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1]) / (std::log(Arb_y[i])
563 G4AutoLock l(&mutex); << 576 - std::log(Arb_y[i - 1]));
564 << 577 Arb_Const[i] = Arb_y[i] / (std::exp(-Arb_x[i] / Arb_ezero[i]));
565 IntType = IType; << 578 Area_seg[i] = -(Arb_Const[i] * Arb_ezero[i]) * (std::exp(-Arb_x[i]
566 ArbEmax = ArbEnergyH.GetMaxEnergy(); << 579 / Arb_ezero[i]) - std::exp(-Arb_x[i - 1] / Arb_ezero[i]));
567 ArbEmin = ArbEnergyH.Energy(0); << 580 } else {
568 << 581 G4cout << "Flat line segment: problem" << G4endl;
569 // Now interpolate points << 582 Arb_ezero[i] = 0.;
570 << 583 Arb_Const[i] = 0.;
571 if (IntType == "Lin") LinearInterpolation(); << 584 Area_seg[i] = 0.;
572 if (IntType == "Log") LogInterpolation(); << 585 }
573 if (IntType == "Exp") ExpInterpolation(); << 586 sum = sum + Area_seg[i];
574 if (IntType == "Spline") SplineInterpolation << 587 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
575 } << 588 if (verbosityLevel == 2)
576 << 589 G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl;
577 void G4SPSEneDistribution::LinearInterpolation << 590 i++;
578 { << 591 }
579 // Method to do linear interpolation on the << 592
580 // Calculate equation of each line segment, << 593 i = 0;
581 // Calculate Area under each segment << 594 while (i < maxi) {
582 // Create a cumulative array which is then n << 595 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
583 << 596 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
584 G4double Area_seg[1024]; // Stores area unde << 597 i++;
585 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 598 }
586 std::size_t i, count; << 599
587 std::size_t maxi = ArbEnergyH.GetVectorLengt << 600 // now scale the ArbEnergyH, needed by Probability()
588 for (i = 0; i < maxi; ++i) << 601 ArbEnergyH.ScaleVector(1., 1./sum);
589 { << 602
590 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 603 if (verbosityLevel >= 1)
591 Arb_y[i] = ArbEnergyH(i); << 604 G4cout << "Leaving ExpInterpolation " << G4endl;
592 } << 605 }
593 << 606
594 // Points are now in x,y arrays. If the spec << 607 void G4SPSEneDistribution::SplineInterpolation() {
595 // made differential and if momentum it has << 608 // Interpolation using Splines.
596 << 609 // Create Normalised arrays, make x 0->1 and y hold
597 if (!DiffSpec) << 610 // the function (Energy)
598 { << 611 //
599 // Converts integral point-wise spectra to << 612 // Current method based on the above will not work in all cases.
600 // << 613 // new method is implemented below.
601 for (count = 0; count < maxi-1; ++count) <<
602 { <<
603 Arb_y[count] = (Arb_y[count] - Arb_y[cou <<
604 / (Arb_x[count + 1] - Arb_x <<
605 } <<
606 --maxi; <<
607 } <<
608 <<
609 if (!EnergySpec) <<
610 { <<
611 // change currently stored values (emin et <<
612 // to energies <<
613 // <<
614 G4ParticleDefinition* pdef = threadLocalDa <<
615 if (pdef == nullptr) <<
616 { <<
617 G4Exception("G4SPSEneDistribution::Linea <<
618 "Event0302", FatalException, <<
619 "Error: particle not defined <<
620 } <<
621 else <<
622 { <<
623 // Apply Energy**2 = p**2c**2 + m0**2c** <<
624 // p should be entered as E/c i.e. witho <<
625 // being done - energy equivalent <<
626 <<
627 G4double mass = pdef->GetPDGMass(); <<
628 <<
629 // Convert point to energy unit and its <<
630 // <<
631 G4double total_energy; <<
632 for (count = 0; count < maxi; ++count) <<
633 { <<
634 total_energy = std::sqrt((Arb_x[count] <<
635 + (mass * mass)); // tota <<
636 Arb_y[count] = Arb_y[count] * Arb_x[co <<
637 Arb_x[count] = total_energy - mass; // <<
638 } <<
639 } <<
640 } <<
641 <<
642 i = 1; <<
643 <<
644 Arb_grad = new G4double [1024]; <<
645 Arb_cept = new G4double [1024]; <<
646 Arb_grad_cept_flag = true; <<
647 <<
648 Arb_grad[0] = 0.; <<
649 Arb_cept[0] = 0.; <<
650 Area_seg[0] = 0.; <<
651 Arb_Cum_Area[0] = 0.; <<
652 while (i < maxi) <<
653 { <<
654 // calculate gradient and intercept for ea <<
655 // <<
656 Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / <<
657 if (verbosityLevel == 2) <<
658 { <<
659 G4cout << Arb_grad[i] << G4endl; <<
660 } <<
661 if (Arb_grad[i] > 0.) <<
662 { <<
663 if (verbosityLevel == 2) <<
664 { <<
665 G4cout << "Arb_grad is positive" << G <<
666 } <<
667 Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * <<
668 } <<
669 else if (Arb_grad[i] < 0.) <<
670 { <<
671 if (verbosityLevel == 2) <<
672 { <<
673 G4cout << "Arb_grad is negative" << G4 <<
674 } <<
675 Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * <<
676 } <<
677 else <<
678 { <<
679 if (verbosityLevel == 2) <<
680 { <<
681 G4cout << "Arb_grad is 0." << G4endl; <<
682 } <<
683 Arb_cept[i] = Arb_y[i]; <<
684 } <<
685 <<
686 Area_seg[i] = ((Arb_grad[i] / 2) <<
687 * (Arb_x[i] * Arb_x[i] - Arb_x <<
688 + Arb_cept[i] * (Arb_x[i] - Ar <<
689 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Ar <<
690 sum = sum + Area_seg[i]; <<
691 if (verbosityLevel == 2) <<
692 { <<
693 G4cout << Arb_x[i] << Arb_y[i] << Area_s <<
694 << Arb_grad[i] << G4endl; <<
695 } <<
696 ++i; <<
697 } <<
698 <<
699 i = 0; <<
700 while (i < maxi) <<
701 { <<
702 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; / <<
703 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_ <<
704 ++i; <<
705 } <<
706 <<
707 // now scale the ArbEnergyH, needed by Proba <<
708 // <<
709 ArbEnergyH.ScaleVector(1., 1./sum); <<
710 <<
711 if (verbosityLevel >= 1) <<
712 { <<
713 G4cout << "Leaving LinearInterpolation" << <<
714 ArbEnergyH.DumpValues(); <<
715 IPDFArbEnergyH.DumpValues(); <<
716 } <<
717 } <<
718 <<
719 void G4SPSEneDistribution::LogInterpolation() <<
720 { <<
721 // Interpolation based on Logarithmic equati <<
722 // Generate equations of line segments <<
723 // y = Ax**alpha => log y = alpha*logx + log <<
724 // Find area under line segments <<
725 // Create normalised, cumulative array Arb_C <<
726 <<
727 G4double Area_seg[1024]; // Stores area unde <<
728 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 <<
729 std::size_t i, count; <<
730 std::size_t maxi = ArbEnergyH.GetVectorLengt <<
731 for (i = 0; i < maxi; ++i) <<
732 { <<
733 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); <<
734 Arb_y[i] = ArbEnergyH(i); <<
735 } <<
736 <<
737 // Points are now in x,y arrays. If the spec <<
738 // made differential and if momentum it has <<
739 <<
740 if (!DiffSpec) <<
741 { <<
742 // Converts integral point-wise spectra to <<
743 // <<
744 for (count = 0; count < maxi-1; ++count) <<
745 { <<
746 Arb_y[count] = (Arb_y[count] - Arb_y[cou <<
747 / (Arb_x[count + 1] - Arb_x <<
748 } <<
749 --maxi; <<
750 } <<
751 <<
752 if (!EnergySpec) <<
753 { <<
754 // Change currently stored values (emin et <<
755 // to energies <<
756 <<
757 G4ParticleDefinition* pdef = threadLocalDa <<
758 if (pdef == nullptr) <<
759 { <<
760 G4Exception("G4SPSEneDistribution::LogIn <<
761 "Event0302", FatalException, <<
762 "Error: particle not defined <<
763 } <<
764 else <<
765 { <<
766 // Apply Energy**2 = p**2c**2 + m0**2c** <<
767 // p should be entered as E/c i.e. witho <<
768 // being done - energy equivalent <<
769 <<
770 G4double mass = pdef->GetPDGMass(); <<
771 <<
772 // Convert point to energy unit and its <<
773 // <<
774 G4double total_energy; <<
775 for (count = 0; count < maxi; ++count) <<
776 { <<
777 total_energy = std::sqrt((Arb_x[count] <<
778 Arb_y[count] = Arb_y[count] * Arb_x[co <<
779 Arb_x[count] = total_energy - mass; // <<
780 } <<
781 } <<
782 } <<
783 <<
784 i = 1; <<
785 <<
786 if ( Arb_ezero != nullptr ) { delete [] Arb_ <<
787 if ( Arb_Const != nullptr ) { delete [] Arb_ <<
788 Arb_alpha = new G4double [1024]; <<
789 Arb_Const = new G4double [1024]; <<
790 Arb_alpha_Const_flag = true; <<
791 <<
792 Arb_alpha[0] = 0.; <<
793 Arb_Const[0] = 0.; <<
794 Area_seg[0] = 0.; <<
795 Arb_Cum_Area[0] = 0.; <<
796 if (Arb_x[0] <= 0. || Arb_y[0] <= 0.) <<
797 { <<
798 G4cout << "You should not use log interpol <<
799 << G4endl; <<
800 G4cout << "These will be changed to 1e-20, <<
801 << G4endl; <<
802 if (Arb_x[0] <= 0.) <<
803 { <<
804 Arb_x[0] = 1e-20; <<
805 } <<
806 if (Arb_y[0] <= 0.) <<
807 { <<
808 Arb_y[0] = 1e-20; <<
809 } <<
810 } <<
811 <<
812 G4double alp; <<
813 while (i < maxi) <<
814 { <<
815 // In case points are negative or zero <<
816 // <<
817 if (Arb_x[i] <= 0. || Arb_y[i] <= 0.) <<
818 { <<
819 G4cout << "You should not use log interp <<
820 << G4endl; <<
821 G4cout << "These will be changed to 1e-2 <<
822 << G4endl; <<
823 if (Arb_x[i] <= 0.) <<
824 { <<
825 Arb_x[i] = 1e-20; <<
826 } <<
827 if (Arb_y[i] <= 0.) <<
828 { <<
829 Arb_y[i] = 1e-20; <<
830 } <<
831 } <<
832 <<
833 Arb_alpha[i] = (std::log10(Arb_y[i]) - std <<
834 / (std::log10(Arb_x[i]) - std <<
835 Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[ <<
836 alp = Arb_alpha[i] + 1; <<
837 if (alp == 0.) <<
838 { <<
839 Area_seg[i] = Arb_Const[i] <<
840 * (std::log(Arb_x[i]) - std: <<
841 } <<
842 else <<
843 { <<
844 Area_seg[i] = (Arb_Const[i] / alp) <<
845 * (std::pow(Arb_x[i], alp) - <<
846 } <<
847 sum = sum + Area_seg[i]; <<
848 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Ar <<
849 if (verbosityLevel == 2) <<
850 { <<
851 G4cout << Arb_alpha[i] << Arb_Const[i] < <<
852 } <<
853 ++i; <<
854 } <<
855 <<
856 i = 0; <<
857 while (i < maxi) <<
858 { <<
859 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; <<
860 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_ <<
861 ++i; <<
862 } <<
863 <<
864 // Now scale the ArbEnergyH, needed by Proba <<
865 // <<
866 ArbEnergyH.ScaleVector(1., 1./sum); <<
867 <<
868 if (verbosityLevel >= 1) <<
869 { <<
870 G4cout << "Leaving LogInterpolation " << G <<
871 } <<
872 } <<
873 <<
874 void G4SPSEneDistribution::ExpInterpolation() <<
875 { <<
876 // Interpolation based on Exponential equati <<
877 // Generate equations of line segments <<
878 // y = Ae**-(x/e0) => ln y = -x/e0 + lnA <<
879 // Find area under line segments <<
880 // Create normalised, cumulative array Arb_C <<
881 <<
882 G4double Area_seg[1024]; // Stores area unde <<
883 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 <<
884 std::size_t i, count; <<
885 std::size_t maxi = ArbEnergyH.GetVectorLengt <<
886 for (i = 0; i < maxi; ++i) <<
887 { <<
888 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); <<
889 Arb_y[i] = ArbEnergyH(i); <<
890 } <<
891 <<
892 // Points are now in x,y arrays. If the spec <<
893 // made differential and if momentum it has <<
894 <<
895 if (!DiffSpec) <<
896 { <<
897 // Converts integral point-wise spectra to <<
898 // <<
899 for (count = 0; count < maxi - 1; ++count) <<
900 { <<
901 Arb_y[count] = (Arb_y[count] - Arb_y[cou <<
902 / (Arb_x[count + 1] - Arb_x <<
903 } <<
904 --maxi; <<
905 } <<
906 <<
907 if (!EnergySpec) <<
908 { <<
909 // Change currently stored values (emin et <<
910 // to energies <<
911 // <<
912 G4ParticleDefinition* pdef = threadLocalDa <<
913 if (pdef == nullptr) <<
914 { <<
915 G4Exception("G4SPSEneDistribution::ExpIn <<
916 "Event0302", FatalException, <<
917 "Error: particle not defined <<
918 } <<
919 else <<
920 { <<
921 // Apply Energy**2 = p**2c**2 + m0**2c** <<
922 // p should be entered as E/c i.e. witho <<
923 // being done - energy equivalent <<
924 <<
925 G4double mass = pdef->GetPDGMass(); <<
926 <<
927 // Convert point to energy unit and its <<
928 // <<
929 G4double total_energy; <<
930 for (count = 0; count < maxi; ++count) <<
931 { <<
932 total_energy = std::sqrt((Arb_x[count] <<
933 Arb_y[count] = Arb_y[count] * Arb_x[co <<
934 Arb_x[count] = total_energy - mass; // <<
935 } <<
936 } <<
937 } <<
938 <<
939 i = 1; <<
940 <<
941 if ( Arb_ezero != nullptr ) { delete[] Arb_e <<
942 if ( Arb_Const != nullptr ) { delete[] Arb_C <<
943 Arb_ezero = new G4double [1024]; <<
944 Arb_Const = new G4double [1024]; <<
945 Arb_ezero_flag = true; <<
946 <<
947 Arb_ezero[0] = 0.; <<
948 Arb_Const[0] = 0.; <<
949 Area_seg[0] = 0.; <<
950 Arb_Cum_Area[0] = 0.; <<
951 while (i < maxi) <<
952 { <<
953 G4double test = std::log(Arb_y[i]) - std:: <<
954 if (test > 0. || test < 0.) <<
955 { <<
956 Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1] <<
957 / (std::log(Arb_y[i]) - std <<
958 Arb_Const[i] = Arb_y[i] / (std::exp(-Arb <<
959 Area_seg[i] = -(Arb_Const[i] * Arb_ezero <<
960 * (std::exp(-Arb_x[i] / Arb_ <<
961 - std::exp(-Arb_x[i - 1] / <<
962 } <<
963 else <<
964 { <<
965 G4Exception("G4SPSEneDistribution::ExpIn <<
966 "Event0302", JustWarning, <<
967 "Flat line segment: problem, <<
968 G4cout << "Flat line segment: problem" < <<
969 Arb_ezero[i] = 0.; <<
970 Arb_Const[i] = 0.; <<
971 Area_seg[i] = 0.; <<
972 } <<
973 sum = sum + Area_seg[i]; <<
974 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Ar <<
975 if (verbosityLevel == 2) <<
976 { <<
977 G4cout << Arb_ezero[i] << Arb_Const[i] < <<
978 } <<
979 ++i; <<
980 } <<
981 <<
982 i = 0; <<
983 while (i < maxi) <<
984 { <<
985 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; <<
986 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_ <<
987 ++i; <<
988 } <<
989 <<
990 // Now scale the ArbEnergyH, needed by Proba <<
991 // <<
992 ArbEnergyH.ScaleVector(1., 1./sum); <<
993 <<
994 if (verbosityLevel >= 1) <<
995 { <<
996 G4cout << "Leaving ExpInterpolation " << G <<
997 } <<
998 } <<
999 <<
1000 void G4SPSEneDistribution::SplineInterpolatio <<
1001 { <<
1002 // Interpolation using Splines. <<
1003 // Create Normalised arrays, make x 0->1 an <<
1004 // <<
1005 // Current method based on the above will n <<
1006 // New method is implemented below. <<
1007 614
1008 G4double sum, Arb_x[1024]={0.}, Arb_y[1024] << 615 G4double sum, Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
1009 std::size_t i, count; << 616 G4int i, count;
1010 std::size_t maxi = ArbEnergyH.GetVectorLeng <<
1011 <<
1012 for (i = 0; i < maxi; ++i) <<
1013 { <<
1014 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i) <<
1015 Arb_y[i] = ArbEnergyH(i); <<
1016 } <<
1017 <<
1018 // Points are now in x,y arrays. If the spe <<
1019 // made differential and if momentum it has <<
1020 <<
1021 if (!DiffSpec) <<
1022 { <<
1023 // Converts integral point-wise spectra t <<
1024 // <<
1025 for (count = 0; count < maxi - 1; ++count <<
1026 { <<
1027 Arb_y[count] = (Arb_y[count] - Arb_y[co <<
1028 / (Arb_x[count + 1] - Arb_ <<
1029 } <<
1030 --maxi; <<
1031 } <<
1032 <<
1033 if (!EnergySpec) <<
1034 { <<
1035 // Change currently stored values (emin e <<
1036 // to energies <<
1037 // <<
1038 G4ParticleDefinition* pdef = threadLocalD <<
1039 if (pdef == nullptr) <<
1040 { <<
1041 G4Exception("G4SPSEneDistribution::Spli <<
1042 "Event0302", FatalException <<
1043 "Error: particle not define <<
1044 } <<
1045 else <<
1046 { <<
1047 // Apply Energy**2 = p**2c**2 + m0**2c* <<
1048 // p should be entered as E/c i.e. with <<
1049 // being done - energy equivalent <<
1050 <<
1051 G4double mass = pdef->GetPDGMass(); <<
1052 <<
1053 // Convert point to energy unit and its <<
1054 // <<
1055 G4double total_energy; <<
1056 for (count = 0; count < maxi; ++count) <<
1057 { <<
1058 total_energy = std::sqrt((Arb_x[count <<
1059 Arb_y[count] = Arb_y[count] * Arb_x[c <<
1060 Arb_x[count] = total_energy - mass; / <<
1061 } <<
1062 } <<
1063 } <<
1064 <<
1065 i = 1; <<
1066 Arb_Cum_Area[0] = 0.; <<
1067 sum = 0.; <<
1068 Splinetemp = new G4DataInterpolation(Arb_x, <<
1069 G4double ei[101], prob[101]; <<
1070 for (auto & it : SplineInt) <<
1071 { <<
1072 delete it; <<
1073 it = nullptr; <<
1074 } <<
1075 SplineInt.clear(); <<
1076 SplineInt.resize(1024,nullptr); <<
1077 while (i < maxi) <<
1078 { <<
1079 // 100 step per segment for the integrati <<
1080 <<
1081 G4double de = (Arb_x[i] - Arb_x[i - 1])/1 <<
1082 G4double area = 0.; <<
1083 <<
1084 for (count = 0; count < 101; ++count) <<
1085 { <<
1086 ei[count] = Arb_x[i - 1] + de*count ; <<
1087 prob[count] = Splinetemp->CubicSplineI <<
1088 if (prob[count] < 0.) <<
1089 { <<
1090 G4ExceptionDescription ED; <<
1091 ED << "Warning: G4DataInterpolation r <<
1092 << " " << ei[count] << G4endl; <<
1093 G4Exception("G4SPSEneDistribution::Sp <<
1094 FatalException, ED); <<
1095 } <<
1096 area += prob[count]*de; <<
1097 } <<
1098 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + a <<
1099 sum += area; <<
1100 <<
1101 prob[0] = prob[0]/(area/de); <<
1102 for (count = 1; count < 100; ++count) <<
1103 { <<
1104 prob[count] = prob[count-1] + prob[coun <<
1105 } <<
1106 <<
1107 SplineInt[i] = new G4DataInterpolation(pr <<
1108 <<
1109 // NOTE: i starts from 1! <<
1110 // <<
1111 ++i; <<
1112 } <<
1113 <<
1114 i = 0; <<
1115 while (i < maxi) <<
1116 { <<
1117 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; <<
1118 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb <<
1119 ++i; <<
1120 } <<
1121 <<
1122 // Now scale the ArbEnergyH, needed by Prob <<
1123 // <<
1124 ArbEnergyH.ScaleVector(1., 1./sum); <<
1125 <<
1126 if (verbosityLevel > 0) <<
1127 { <<
1128 G4cout << "Leaving SplineInterpolation " <<
1129 } <<
1130 } <<
1131 <<
1132 void G4SPSEneDistribution::GenerateMonoEnerge <<
1133 { <<
1134 // Method to generate MonoEnergetic particl <<
1135 <<
1136 threadLocalData.Get().particle_energy = Mon <<
1137 } <<
1138 617
1139 void G4SPSEneDistribution::GenerateGaussEnerg << 618 G4int maxi = ArbEnergyH.GetVectorLength();
1140 { << 619 for (i = 0; i < maxi; i++) {
1141 // Method to generate Gaussian particles << 620 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
1142 << 621 Arb_y[i] = ArbEnergyH(size_t(i));
1143 G4double ene = G4RandGauss::shoot(MonoEnerg << 622 }
1144 if (ene < 0) ene = 0.; << 623 // Points are now in x,y arrays. If the spectrum is integral it has to be
1145 threadLocalData.Get().particle_energy = ene << 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.) {
>> 673 G4ExceptionDescription ED;
>> 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 }
1146 } 1196 }
1147 1197
1148 void G4SPSEneDistribution::GenerateLinearEner << 1198 //
1149 { << 1199 void G4SPSEneDistribution::ReSetHist(G4String atype) {
1150 G4double rndm; << 1200 if (atype == "energy") {
1151 threadLocal_t& params = threadLocalData.Get << 1201 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
1152 G4double emaxsq = std::pow(params.Emax, 2.) << 1202 IPDFEnergyExist = false;
1153 G4double eminsq = std::pow(params.Emin, 2.) << 1203 Emin = 0.;
1154 G4double intersq = std::pow(params.cept, 2. << 1204 Emax = 1e30;
1155 << 1205 } else if (atype == "arb") {
1156 if (bArb) rndm = G4UniformRand(); << 1206 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
1157 else rndm = eneRndm->GenRandEnergy(); << 1207 IPDFArbExist = false;
1158 << 1208 } else if (atype == "epn") {
1159 G4double bracket = ((params.grad / 2.) << 1209 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
1160 * (emaxsq - eminsq) << 1210 IPDFEnergyExist = false;
1161 + params.cept * (params.Em << 1211 EpnEnergyH = ZeroPhysVector;
1162 bracket = bracket * rndm; << 1212 } else {
1163 bracket = bracket + (params.grad / 2.) * em << 1213 G4cout << "Error, histtype not accepted " << G4endl;
1164 << 1214 }
1165 // Now have a quad of form m/2 E**2 + cE - << 1215 }
1166 // << 1216
1167 bracket = -bracket; << 1217 G4double G4SPSEneDistribution::GenerateOne(G4ParticleDefinition* a) {
1168 << 1218 particle_definition = a;
1169 if (params.grad != 0.) << 1219 particle_energy = -1.;
1170 { << 1220
1171 G4double sqbrack = (intersq - 4 * (params << 1221 while ((EnergyDisType == "Arb") ? (particle_energy < ArbEmin
1172 sqbrack = std::sqrt(sqbrack); << 1222 || particle_energy > ArbEmax) : (particle_energy < Emin
1173 G4double root1 = -params.cept + sqbrack; << 1223 || particle_energy > Emax)) {
1174 root1 = root1 / (2. * (params.grad / 2.)) << 1224 if (Biased) {
1175 << 1225 GenerateBiasPowEnergies();
1176 G4double root2 = -params.cept - sqbrack; << 1226 } else {
1177 root2 = root2 / (2. * (params.grad / 2.)) << 1227 if (EnergyDisType == "Mono")
1178 << 1228 GenerateMonoEnergetic();
1179 if (root1 > params.Emin && root1 < params << 1229 else if (EnergyDisType == "Lin")
1180 { << 1230 GenerateLinearEnergies();
1181 params.particle_energy = root1; << 1231 else if (EnergyDisType == "Pow")
1182 } << 1232 GeneratePowEnergies();
1183 if (root2 > params.Emin && root2 < params << 1233 else if (EnergyDisType == "Exp")
1184 { << 1234 GenerateExpEnergies();
1185 params.particle_energy = root2; << 1235 else if (EnergyDisType == "Gauss")
1186 } << 1236 GenerateGaussEnergies();
1187 } << 1237 else if (EnergyDisType == "Brem")
1188 else if (params.grad == 0.) << 1238 GenerateBremEnergies();
1189 { << 1239 else if (EnergyDisType == "Bbody")
1190 // have equation of form cE - bracket =0 << 1240 GenerateBbodyEnergies();
1191 // << 1241 else if (EnergyDisType == "Cdg")
1192 params.particle_energy = bracket / params << 1242 GenerateCdgEnergies();
1193 } << 1243 else if (EnergyDisType == "User")
1194 << 1244 GenUserHistEnergies();
1195 if (params.particle_energy < 0.) << 1245 else if (EnergyDisType == "Arb")
1196 { << 1246 GenArbPointEnergies();
1197 params.particle_energy = -params.particle << 1247 else if (EnergyDisType == "Epn")
1198 } << 1248 GenEpnHistEnergies();
1199 << 1249 else
1200 if (verbosityLevel >= 1) << 1250 G4cout << "Error: EnergyDisType has unusual value" << G4endl;
1201 { << 1251 }
1202 G4cout << "Energy is " << params.particle << 1252 }
1203 } << 1253 return particle_energy;
1204 } << 1254 }
1205 << 1255
1206 void G4SPSEneDistribution::GeneratePowEnergie << 1256 G4double G4SPSEneDistribution::GetProbability(G4double ene) {
1207 { << 1257 G4double prob = 1.;
1208 // Method to generate particle energies dis << 1258
1209 << 1259 if (EnergyDisType == "Lin") {
1210 G4double rndm; << 1260 if (prob_norm == 1.) {
1211 G4double emina, emaxa; << 1261 prob_norm = 0.5*grad*Emax*Emax + cept*Emax - 0.5*grad*Emin*Emin - cept*Emin;
1212 << 1262 }
1213 threadLocal_t& params = threadLocalData.Get << 1263 prob = cept + grad * ene;
1214 << 1264 prob /= prob_norm;
1215 emina = std::pow(params.Emin, params.alpha << 1265 }
1216 emaxa = std::pow(params.Emax, params.alpha << 1266 else if (EnergyDisType == "Pow") {
1217 << 1267 if (prob_norm == 1.) {
1218 if (bArb) rndm = G4UniformRand(); << 1268 if (alpha != -1.) {
1219 else rndm = eneRndm->GenRandEnergy(); << 1269 G4double emina = std::pow(Emin, alpha + 1);
1220 << 1270 G4double emaxa = std::pow(Emax, alpha + 1);
1221 if (params.alpha != -1.) << 1271 prob_norm = 1./(1.+alpha) * (emaxa - emina);
1222 { << 1272 } else {
1223 G4double ene = ((rndm * (emaxa - emina)) << 1273 prob_norm = std::log(Emax) - std::log(Emin) ;
1224 ene = std::pow(ene, (1. / (params.alpha + << 1274 }
1225 params.particle_energy = ene; << 1275 }
1226 } << 1276 prob = std::pow(ene, alpha)/prob_norm;
1227 else << 1277 }
1228 { << 1278 else if (EnergyDisType == "Exp"){
1229 G4double ene = (std::log(params.Emin) << 1279 if (prob_norm == 1.) {
1230 + rndm * (std::log(params.Em << 1280 prob_norm = -Ezero*(std::exp(-Emax/Ezero) - std::exp(Emin/Ezero));
1231 params.particle_energy = std::exp(ene); << 1281 }
1232 } << 1282 prob = std::exp(-ene / Ezero);
1233 if (verbosityLevel >= 1) << 1283 prob /= prob_norm;
1234 { << 1284 }
1235 G4cout << "Energy is " << params.particle << 1285 else if (EnergyDisType == "Arb") {
1236 } << 1286 prob = ArbEnergyH.Value(ene);
1237 } << 1287 // prob = ArbEInt->CubicSplineInterpolation(ene);
1238 << 1288 //G4double deltaY;
1239 void G4SPSEneDistribution::GenerateCPowEnergi << 1289 //prob = ArbEInt->PolynomInterpolation(ene, deltaY);
1240 { << 1290 if (prob <= 0.) {
1241 // Method to generate particle energies dis << 1291 //G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << " " <<deltaY<< G4endl;
1242 // cutoff power-law distribution << 1292 G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << G4endl;
1243 // << 1293 prob = 1e-30;
1244 // CP_x holds Energies, and CPHist holds th << 1294 }
1245 // binary search to find correct bin then l << 1295 // already normalised
1246 // Use the earlier defined histogram + Rand << 1296 }
1247 // random numbers following the histos dist << 1297 else
1248 << 1298 G4cout << "Error: EnergyDisType not supported" << G4endl;
1249 G4double rndm = eneRndm->GenRandEnergy(); <<
1250 G4int nabove = 10001, nbelow = 0, middle; <<
1251 <<
1252 G4AutoLock l(&mutex); <<
1253 G4bool done = CPhistCalcd; <<
1254 l.unlock(); <<
1255 <<
1256 if(!done) <<
1257 { <<
1258 Calculate(); //This is has a lock inside, <<
1259 l.lock(); <<
1260 CPhistCalcd = true; <<
1261 l.unlock(); <<
1262 } <<
1263 <<
1264 // Binary search to find bin that rndm is i <<
1265 // <<
1266 while (nabove - nbelow > 1) <<
1267 { <<
1268 middle = (nabove + nbelow) / 2; <<
1269 if (rndm == CPHist->at(middle)) <<
1270 { <<
1271 break; <<
1272 } <<
1273 if (rndm < CPHist->at(middle)) <<
1274 { <<
1275 nabove = middle; <<
1276 } <<
1277 else <<
1278 { <<
1279 nbelow = middle; <<
1280 } <<
1281 } <<
1282 <<
1283 // Now interpolate in that bin to find the <<
1284 // <<
1285 G4double x1, x2, y1, y2, t, q; <<
1286 x1 = CP_x->at(nbelow); <<
1287 if(nbelow+1 == static_cast<G4int>(CP_x->siz <<
1288 { <<
1289 x2 = CP_x->back(); <<
1290 } <<
1291 else <<
1292 { <<
1293 x2 = CP_x->at(nbelow + 1); <<
1294 } <<
1295 y1 = CPHist->at(nbelow); <<
1296 if(nbelow+1 == static_cast<G4int>(CPHist->s <<
1297 { <<
1298 G4cout << CPHist->back() << G4endl; <<
1299 y2 = CPHist->back(); <<
1300 } <<
1301 else <<
1302 { <<
1303 y2 = CPHist->at(nbelow + 1); <<
1304 } <<
1305 t = (y2 - y1) / (x2 - x1); <<
1306 q = y1 - t * x1; <<
1307 <<
1308 threadLocalData.Get().particle_energy = (rn <<
1309 <<
1310 if (verbosityLevel >= 1) <<
1311 { <<
1312 G4cout << "Energy is " << threadLocalData <<
1313 } <<
1314 } <<
1315 <<
1316 void G4SPSEneDistribution::GenerateBiasPowEne <<
1317 { <<
1318 // Method to generate particle energies dis <<
1319 // in biased power-law and calculate its we <<
1320 <<
1321 threadLocal_t& params = threadLocalData.Get <<
1322 <<
1323 G4double rndm; <<
1324 G4double emina, emaxa, emin, emax; <<
1325 <<
1326 G4double normal = 1.; <<
1327 <<
1328 emin = params.Emin; <<
1329 emax = params.Emax; <<
1330 <<
1331 rndm = eneRndm->GenRandEnergy(); <<
1332 <<
1333 if (biasalpha != -1.) <<
1334 { <<
1335 emina = std::pow(emin, biasalpha + 1); <<
1336 emaxa = std::pow(emax, biasalpha + 1); <<
1337 G4double ee = ((rndm * (emaxa - emina)) + <<
1338 params.particle_energy = std::pow(ee, (1. <<
1339 normal = 1./(1+biasalpha) * (emaxa - emin <<
1340 } <<
1341 else <<
1342 { <<
1343 G4double ee = (std::log(emin) + rndm * (s <<
1344 params.particle_energy = std::exp(ee); <<
1345 normal = std::log(emax) - std::log(emin); <<
1346 } <<
1347 params.weight = GetProbability(params.parti <<
1348 / (std::pow(params.particle_e <<
1349 <<
1350 if (verbosityLevel >= 1) <<
1351 { <<
1352 G4cout << "Energy is " << params.particle <<
1353 } <<
1354 } <<
1355 <<
1356 void G4SPSEneDistribution::GenerateExpEnergie <<
1357 { <<
1358 // Method to generate particle energies dis <<
1359 // to an exponential curve <<
1360 <<
1361 G4double rndm; <<
1362 <<
1363 if (bArb) rndm = G4UniformRand(); <<
1364 else rndm = eneRndm->GenRandEnergy(); <<
1365 <<
1366 threadLocal_t& params = threadLocalData.Get <<
1367 params.particle_energy = -params.Ezero <<
1368 * (std::log(rndm * ( <<
1369 <<
1370 - <<
1371 <<
1372 + std::exp <<
1373 if (verbosityLevel >= 1) <<
1374 { <<
1375 G4cout << "Energy is " << params.particle <<
1376 } <<
1377 } <<
1378 <<
1379 void G4SPSEneDistribution::GenerateBremEnergi <<
1380 { <<
1381 // Method to generate particle energies dis <<
1382 // to a Bremstrahlung equation of the form <<
1383 // I = const*((kT)**1/2)*E*(e**(-E/kT)) <<
1384 <<
1385 G4double rndm = eneRndm->GenRandEnergy(); <<
1386 G4double expmax, expmin, k; <<
1387 <<
1388 k = 8.6181e-11; // Boltzmann's const in MeV <<
1389 G4double ksq = std::pow(k, 2.); // k square <<
1390 G4double Tsq = std::pow(Temp, 2.); // Temp <<
1391 <<
1392 threadLocal_t& params = threadLocalData.Get <<
1393 <<
1394 expmax = std::exp(-params.Emax / (k * Temp) <<
1395 expmin = std::exp(-params.Emin / (k * Temp) <<
1396 <<
1397 // If either expmax or expmin are zero then <<
1398 // Most probably this will be because T is <<
1399 <<
1400 if (expmax == 0.) <<
1401 { <<
1402 G4Exception("G4SPSEneDistribution::Genera <<
1403 "Event0302", FatalException, <<
1404 "*****EXPMAX=0. Choose differ <<
1405 } <<
1406 if (expmin == 0.) <<
1407 { <<
1408 G4Exception("G4SPSEneDistribution::Genera <<
1409 "Event0302", FatalException, <<
1410 "*****EXPMIN=0. Choose differ <<
1411 } <<
1412 <<
1413 G4double tempvar = rndm * ((-k) * Temp * (p <<
1414 - p <<
1415 - (ksq * Tsq * (exp <<
1416 <<
1417 G4double bigc = (tempvar - k * Temp * param <<
1418 - ksq * Tsq * expmin) / (-k <<
1419 <<
1420 // This gives an equation of form: Ee(-E/kT <<
1421 // Solve this iteratively, step from Emin t <<
1422 // and take the best solution. <<
1423 <<
1424 G4double erange = params.Emax - params.Emin <<
1425 G4double steps = erange / 1000.; <<
1426 G4int i; <<
1427 G4double etest, diff, err = 100000.; <<
1428 <<
1429 for (i = 1; i < 1000; ++i) <<
1430 { <<
1431 etest = params.Emin + (i - 1) * steps; <<
1432 diff = etest * (std::exp(-etest / (k * Te <<
1433 + k * Temp * (std::exp(-etest / (k * <<
1434 <<
1435 if (diff < 0.) <<
1436 { <<
1437 diff = -diff; <<
1438 } <<
1439 <<
1440 if (diff < err) <<
1441 { <<
1442 err = diff; <<
1443 params.particle_energy = etest; <<
1444 } <<
1445 } <<
1446 if (verbosityLevel >= 1) <<
1447 { <<
1448 G4cout << "Energy is " << params.particle <<
1449 } <<
1450 } <<
1451 <<
1452 void G4SPSEneDistribution::GenerateBbodyEnerg <<
1453 { <<
1454 // BBody_x holds Energies, and BBHist holds <<
1455 // Binary search to find correct bin then l <<
1456 // Use the earlier defined histogram + Rand <<
1457 // random numbers following the histos dist <<
1458 <<
1459 G4double rndm = eneRndm->GenRandEnergy(); <<
1460 G4int nabove = 10001, nbelow = 0, middle; <<
1461 <<
1462 G4AutoLock l(&mutex); <<
1463 G4bool done = BBhistCalcd; <<
1464 l.unlock(); <<
1465 <<
1466 if(!done) <<
1467 { <<
1468 Calculate(); //This is has a lock inside, <<
1469 l.lock(); <<
1470 BBhistCalcd = true; <<
1471 l.unlock(); <<
1472 } <<
1473 <<
1474 // Binary search to find bin that rndm is i <<
1475 // <<
1476 while (nabove - nbelow > 1) <<
1477 { <<
1478 middle = (nabove + nbelow) / 2; <<
1479 if (rndm == BBHist->at(middle)) <<
1480 { <<
1481 break; <<
1482 } <<
1483 if (rndm < BBHist->at(middle)) <<
1484 { <<
1485 nabove = middle; <<
1486 } <<
1487 else <<
1488 { <<
1489 nbelow = middle; <<
1490 } <<
1491 } <<
1492 <<
1493 // Now interpolate in that bin to find the <<
1494 // <<
1495 G4double x1, x2, y1, y2, t, q; <<
1496 x1 = Bbody_x->at(nbelow); <<
1497 <<
1498 if(nbelow+1 == static_cast<G4int>(Bbody_x-> <<
1499 { <<
1500 x2 = Bbody_x->back(); <<
1501 } <<
1502 else <<
1503 { <<
1504 x2 = Bbody_x->at(nbelow + 1); <<
1505 } <<
1506 y1 = BBHist->at(nbelow); <<
1507 if(nbelow+1 == static_cast<G4int>(BBHist->s <<
1508 { <<
1509 G4cout << BBHist->back() << G4endl; <<
1510 y2 = BBHist->back(); <<
1511 } <<
1512 else <<
1513 { <<
1514 y2 = BBHist->at(nbelow + 1); <<
1515 } <<
1516 t = (y2 - y1) / (x2 - x1); <<
1517 q = y1 - t * x1; <<
1518 <<
1519 threadLocalData.Get().particle_energy = (rn <<
1520 <<
1521 if (verbosityLevel >= 1) <<
1522 { <<
1523 G4cout << "Energy is " << threadLocalData <<
1524 } <<
1525 } <<
1526 <<
1527 void G4SPSEneDistribution::GenerateCdgEnergie <<
1528 { <<
1529 // Generate random numbers, compare with va <<
1530 // to find appropriate part of spectrum and <<
1531 // generate energy in the usual inversion w <<
1532 <<
1533 G4double rndm, rndm2; <<
1534 G4double ene_line[3]={0,0,0}; <<
1535 G4double omalpha[2]={0,0}; <<
1536 threadLocal_t& params = threadLocalData.Get <<
1537 if (params.Emin < 18 * keV && params.Emax < <<
1538 { <<
1539 omalpha[0] = 1. - 1.4; <<
1540 ene_line[0] = params.Emin; <<
1541 ene_line[1] = params.Emax; <<
1542 } <<
1543 if (params.Emin < 18 * keV && params.Emax > <<
1544 { <<
1545 omalpha[0] = 1. - 1.4; <<
1546 omalpha[1] = 1. - 2.3; <<
1547 ene_line[0] = params.Emin; <<
1548 ene_line[1] = 18. * keV; <<
1549 ene_line[2] = params.Emax; <<
1550 } <<
1551 if (params.Emin > 18 * keV) <<
1552 { <<
1553 omalpha[0] = 1. - 2.3; <<
1554 ene_line[0] = params.Emin; <<
1555 ene_line[1] = params.Emax; <<
1556 } <<
1557 rndm = eneRndm->GenRandEnergy(); <<
1558 rndm2 = eneRndm->GenRandEnergy(); <<
1559 <<
1560 G4int i = 0; <<
1561 while (rndm >= CDGhist[i]) <<
1562 { <<
1563 ++i; <<
1564 } <<
1565 <<
1566 // Generate final energy <<
1567 // <<
1568 G4double ene = (std::pow(ene_line[i - 1], o <<
1569 + (std::pow(ene_line[i], omalp <<
1570 - std::pow(ene_line[i - 1], o <<
1571 params.particle_energy = std::pow(ene, (1. <<
1572 <<
1573 if (verbosityLevel >= 1) <<
1574 { <<
1575 G4cout << "Energy is " << params.particle <<
1576 } <<
1577 } <<
1578 <<
1579 void G4SPSEneDistribution::GenUserHistEnergie <<
1580 { <<
1581 // Histograms are DIFFERENTIAL <<
1582 <<
1583 G4AutoLock l(&mutex); <<
1584 <<
1585 if (!IPDFEnergyExist) <<
1586 { <<
1587 std::size_t ii; <<
1588 std::size_t maxbin = UDefEnergyH.GetVecto <<
1589 G4double bins[1024], vals[1024], sum; <<
1590 for ( ii = 0 ; ii<1024 ; ++ii ) { bins[ii <<
1591 sum = 0.; <<
1592 <<
1593 if ( (!EnergySpec) <<
1594 && (threadLocalData.Get().particle_defi <<
1595 { <<
1596 G4Exception("G4SPSEneDistribution::GenU <<
1597 "Event0302", FatalException <<
1598 "Error: particle definition <<
1599 } <<
1600 <<
1601 if (maxbin > 1024) <<
1602 { <<
1603 G4Exception("G4SPSEneDistribution::GenU <<
1604 "Event0302", JustWarning, <<
1605 "Maxbin>1024\n Setting maxbi <<
1606 maxbin = 1024; <<
1607 } <<
1608 <<
1609 if (!DiffSpec) <<
1610 { <<
1611 G4cout << "Histograms are Differential! <<
1612 } <<
1613 else <<
1614 { <<
1615 bins[0] = UDefEnergyH.GetLowEdgeEnergy( <<
1616 vals[0] = UDefEnergyH(0); <<
1617 sum = vals[0]; <<
1618 for (ii = 1; ii < maxbin; ++ii) <<
1619 { <<
1620 bins[ii] = UDefEnergyH.GetLowEdgeEner <<
1621 vals[ii] = UDefEnergyH(ii) + vals[ii <<
1622 sum = sum + UDefEnergyH(ii); <<
1623 } <<
1624 } <<
1625 <<
1626 if (!EnergySpec) <<
1627 { <<
1628 G4double mass = threadLocalData.Get().p <<
1629 <<
1630 // Multiply the function (vals) up by t <<
1631 // to make the function counts/s (i.e. <<
1632 <<
1633 for (ii = 1; ii < maxbin; ++ii) <<
1634 { <<
1635 vals[ii] = vals[ii] * (bins[ii] - bin <<
1636 } <<
1637 <<
1638 // Put energy bins into new histo, plus <<
1639 // to make evals counts/s/energy <<
1640 // <<
1641 for (ii = 0; ii < maxbin; ++ii) <<
1642 { <<
1643 // kinetic energy <<
1644 // <<
1645 bins[ii] = std::sqrt((bins[ii]*bins[i <<
1646 } <<
1647 for (ii = 1; ii < maxbin; ++ii) <<
1648 { <<
1649 vals[ii] = vals[ii] / (bins[ii] - bin <<
1650 } <<
1651 sum = vals[maxbin - 1]; <<
1652 vals[0] = 0.; <<
1653 } <<
1654 for (ii = 0; ii < maxbin; ++ii) <<
1655 { <<
1656 vals[ii] = vals[ii] / sum; <<
1657 IPDFEnergyH.InsertValues(bins[ii], vals <<
1658 } <<
1659 <<
1660 IPDFEnergyExist = true; <<
1661 if (verbosityLevel > 1) <<
1662 { <<
1663 IPDFEnergyH.DumpValues(); <<
1664 } <<
1665 } <<
1666 l.unlock(); <<
1667 <<
1668 // IPDF has been create so carry on <<
1669 // <<
1670 G4double rndm = eneRndm->GenRandEnergy(); <<
1671 threadLocalData.Get().particle_energy= IPDF <<
1672 <<
1673 if (verbosityLevel >= 1) <<
1674 { <<
1675 G4cout << "Energy is " << particle_energy <<
1676 } <<
1677 } <<
1678 <<
1679 G4double G4SPSEneDistribution::GetArbEneWeigh <<
1680 { <<
1681 auto nbelow = IPDFArbEnergyH.FindBin(ene,(I <<
1682 G4double wei = 0.; <<
1683 if(IntType=="Lin") <<
1684 { <<
1685 // note: grad[i] and cept[i] are calculat <<
1686 auto gr = Arb_grad[nbelow + 1]; <<
1687 auto ce = Arb_cept[nbelow + 1]; <<
1688 wei = ene*gr + ce; <<
1689 } <<
1690 else if(IntType=="Log") <<
1691 { <<
1692 auto alp = Arb_alpha[nbelow + 1]; <<
1693 auto cns = Arb_Const[nbelow + 1]; <<
1694 wei = cns * std::pow(ene,alp); <<
1695 } <<
1696 else if(IntType=="Exp") <<
1697 { <<
1698 auto e0 = Arb_ezero[nbelow + 1]; <<
1699 auto cns = Arb_Const[nbelow + 1]; <<
1700 wei = cns * std::exp(-ene/e0); <<
1701 } <<
1702 else if(IntType=="Spline") <<
1703 { <<
1704 wei = SplineInt[nbelow+1]->CubicSplineInt <<
1705 } <<
1706 return wei; <<
1707 } <<
1708 <<
1709 void G4SPSEneDistribution::GenArbPointEnergie <<
1710 { <<
1711 if (verbosityLevel > 0) <<
1712 { <<
1713 G4cout << "In GenArbPointEnergies" << G4e <<
1714 } <<
1715 <<
1716 G4double rndm = eneRndm->GenRandEnergy(); <<
1717 <<
1718 // Find the Bin, have x, y, no of points, a <<
1719 // <<
1720 std::size_t nabove = IPDFArbEnergyH.GetVect <<
1721 <<
1722 // Binary search to find bin that rndm is i <<
1723 // <<
1724 while (nabove - nbelow > 1) <<
1725 { <<
1726 middle = (nabove + nbelow) / 2; <<
1727 if (rndm == IPDFArbEnergyH(middle)) <<
1728 { <<
1729 break; <<
1730 } <<
1731 if (rndm < IPDFArbEnergyH(middle)) <<
1732 { <<
1733 nabove = middle; <<
1734 } <<
1735 else <<
1736 { <<
1737 nbelow = middle; <<
1738 } <<
1739 } <<
1740 threadLocal_t& params = threadLocalData.Get <<
1741 if (IntType == "Lin") <<
1742 { <<
1743 // Update thread-local copy of parameters <<
1744 // <<
1745 params.Emax = IPDFArbEnergyH.GetLowEdgeEn <<
1746 params.Emin = IPDFArbEnergyH.GetLowEdgeEn <<
1747 params.grad = Arb_grad[nbelow + 1]; <<
1748 params.cept = Arb_cept[nbelow + 1]; <<
1749 GenerateLinearEnergies(true); <<
1750 } <<
1751 else if (IntType == "Log") <<
1752 { <<
1753 params.Emax = IPDFArbEnergyH.GetLowEdgeEn <<
1754 params.Emin = IPDFArbEnergyH.GetLowEdgeEn <<
1755 params.alpha = Arb_alpha[nbelow + 1]; <<
1756 GeneratePowEnergies(true); <<
1757 } <<
1758 else if (IntType == "Exp") <<
1759 { <<
1760 params.Emax = IPDFArbEnergyH.GetLowEdgeEn <<
1761 params.Emin = IPDFArbEnergyH.GetLowEdgeEn <<
1762 params.Ezero = Arb_ezero[nbelow + 1]; <<
1763 GenerateExpEnergies(true); <<
1764 } <<
1765 else if (IntType == "Spline") <<
1766 { <<
1767 params.Emax = IPDFArbEnergyH.GetLowEdgeEn <<
1768 params.Emin = IPDFArbEnergyH.GetLowEdgeEn <<
1769 params.particle_energy = -1e100; <<
1770 rndm = eneRndm->GenRandEnergy(); <<
1771 while (params.particle_energy < params.Em <<
1772 || params.particle_energy > params.Em <<
1773 { <<
1774 params.particle_energy = <<
1775 SplineInt[nbelow+1]->CubicSplineInter <<
1776 rndm = eneRndm->GenRandEnergy(); <<
1777 } <<
1778 if (verbosityLevel >= 1) <<
1779 { <<
1780 G4cout << "Energy is " << params.partic <<
1781 } <<
1782 } <<
1783 else <<
1784 { <<
1785 G4Exception("G4SPSEneDistribution::GenArb <<
1786 FatalException, "Error: IntTy <<
1787 } <<
1788 } <<
1789 <<
1790 void G4SPSEneDistribution::GenEpnHistEnergies <<
1791 { <<
1792 // Firstly convert to energy if not already <<
1793 <<
1794 G4AutoLock l(&mutex); <<
1795 <<
1796 if (Epnflag) // true means spectrum is epn <<
1797 { <<
1798 // Convert to energy by multiplying by A <<
1799 // <<
1800 ConvertEPNToEnergy(); <<
1801 } <<
1802 if (!IPDFEnergyExist) <<
1803 { <<
1804 // IPDF has not been created, so create i <<
1805 // <<
1806 G4double bins[1024], vals[1024], sum; <<
1807 std::size_t ii; <<
1808 std::size_t maxbin = UDefEnergyH.GetVecto <<
1809 bins[0] = UDefEnergyH.GetLowEdgeEnergy(0) <<
1810 vals[0] = UDefEnergyH(0); <<
1811 sum = vals[0]; <<
1812 for (ii = 1; ii < maxbin; ++ii) <<
1813 { <<
1814 bins[ii] = UDefEnergyH.GetLowEdgeEnergy <<
1815 vals[ii] = UDefEnergyH(ii) + vals[ii - <<
1816 sum = sum + UDefEnergyH(ii); <<
1817 } <<
1818 <<
1819 l.lock(); <<
1820 for (ii = 0; ii < maxbin; ++ii) <<
1821 { <<
1822 vals[ii] = vals[ii] / sum; <<
1823 IPDFEnergyH.InsertValues(bins[ii], vals <<
1824 } <<
1825 IPDFEnergyExist = true; <<
1826 1299
1827 } << 1300 return prob;
1828 l.unlock(); <<
1829 <<
1830 // IPDF has been create so carry on <<
1831 // <<
1832 G4double rndm = eneRndm->GenRandEnergy(); <<
1833 threadLocalData.Get().particle_energy = IPD <<
1834 <<
1835 if (verbosityLevel >= 1) <<
1836 { <<
1837 G4cout << "Energy is " << threadLocalData <<
1838 } <<
1839 } <<
1840 <<
1841 void G4SPSEneDistribution::ConvertEPNToEnergy <<
1842 { <<
1843 // Use this before particle generation to c <<
1844 // currently stored histogram from energy/n <<
1845 <<
1846 threadLocal_t& params = threadLocalData.Get <<
1847 if (params.particle_definition == nullptr) <<
1848 { <<
1849 G4cout << "Error: particle not defined" < <<
1850 } <<
1851 else <<
1852 { <<
1853 // Need to multiply histogram by the numb <<
1854 // Baryon Number looks to hold the no. of <<
1855 // <<
1856 G4int Bary = params.particle_definition-> <<
1857 <<
1858 // Change values in histogram, Read it ou <<
1859 // <<
1860 std::size_t count, maxcount; <<
1861 maxcount = EpnEnergyH.GetVectorLength(); <<
1862 G4double ebins[1024], evals[1024]; <<
1863 if (maxcount > 1024) <<
1864 { <<
1865 G4Exception("G4SPSEneDistribution::Conv <<
1866 "gps001", JustWarning, <<
1867 "Histogram contains more th <<
1868 Those above 1024 will be i <<
1869 maxcount = 1024; <<
1870 } <<
1871 if (maxcount < 1) <<
1872 { <<
1873 G4Exception("G4SPSEneDistribution::Conv <<
1874 "gps001", FatalException, <<
1875 "Histogram contains less tha <<
1876 return; <<
1877 } <<
1878 for (count = 0; count < maxcount; ++count <<
1879 { <<
1880 // Read out <<
1881 ebins[count] = EpnEnergyH.GetLowEdgeEne <<
1882 evals[count] = EpnEnergyH(count); <<
1883 } <<
1884 <<
1885 // Multiply the channels by the nucleon n <<
1886 // <<
1887 for (count = 0; count < maxcount; ++count <<
1888 { <<
1889 ebins[count] = ebins[count] * Bary; <<
1890 } <<
1891 <<
1892 // Set Emin and Emax <<
1893 // <<
1894 params.Emin = ebins[0]; <<
1895 if (maxcount > 1) <<
1896 { <<
1897 params.Emax = ebins[maxcount - 1]; <<
1898 } <<
1899 else <<
1900 { <<
1901 params.Emax = ebins[0]; <<
1902 } <<
1903 <<
1904 // Put energy bins into new histogram - U <<
1905 // <<
1906 for (count = 0; count < maxcount; ++count <<
1907 { <<
1908 UDefEnergyH.InsertValues(ebins[count], <<
1909 } <<
1910 Epnflag = false; // so that you dont repe <<
1911 } <<
1912 } <<
1913 <<
1914 void G4SPSEneDistribution::ReSetHist(const G4 <<
1915 { <<
1916 G4AutoLock l(&mutex); <<
1917 if (atype == "energy") <<
1918 { <<
1919 UDefEnergyH = IPDFEnergyH = ZeroPhysVecto <<
1920 IPDFEnergyExist = false; <<
1921 Emin = 0.; <<
1922 Emax = 1e30; <<
1923 } <<
1924 else if (atype == "arb") <<
1925 { <<
1926 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVec <<
1927 IPDFArbExist = false; <<
1928 } <<
1929 else if (atype == "epn") <<
1930 { <<
1931 UDefEnergyH = IPDFEnergyH = ZeroPhysVecto <<
1932 IPDFEnergyExist = false; <<
1933 EpnEnergyH = ZeroPhysVector; <<
1934 } <<
1935 else <<
1936 { <<
1937 G4cout << "Error, histtype not accepted " <<
1938 } <<
1939 } <<
1940 <<
1941 G4double G4SPSEneDistribution::GenerateOne(G4 <<
1942 { <<
1943 // Copy global shared status to thread-loca <<
1944 // <<
1945 threadLocal_t& params = threadLocalData.Get <<
1946 params.particle_definition=a; <<
1947 params.particle_energy=-1; <<
1948 if(applyEvergyWeight) <<
1949 { <<
1950 params.Emax = ArbEmax; <<
1951 params.Emin = ArbEmin; <<
1952 } <<
1953 else <<
1954 { <<
1955 params.Emax = Emax; <<
1956 params.Emin = Emin; <<
1957 } <<
1958 params.alpha = alpha; <<
1959 params.Ezero = Ezero; <<
1960 params.grad = grad; <<
1961 params.cept = cept; <<
1962 params.weight = weight; <<
1963 // particle_energy = -1.; <<
1964 <<
1965 if((EnergyDisType == "Mono") && ((MonoEnerg <<
1966 { <<
1967 G4ExceptionDescription ed; <<
1968 ed << "MonoEnergy " << G4BestUnit(MonoEne <<
1969 << " is outside of [Emin,Emax] = [" <<
1970 << G4BestUnit(Emin,"Energy") << ", " <<
1971 << G4BestUnit(Emax,"Energy") << ". Mon <<
1972 G4Exception("G4SPSEneDistribution::Genera <<
1973 "GPS0001", JustWarning, ed); <<
1974 params.particle_energy=MonoEnergy; <<
1975 return params.particle_energy; <<
1976 } <<
1977 while ( (EnergyDisType == "Arb") <<
1978 ? (params.particle_energy < ArbEmin <<
1979 || params.particle_energy > ArbEmax <<
1980 : (params.particle_energy < params. <<
1981 || params.particle_energy > params. <<
1982 { <<
1983 if (Biased) <<
1984 { <<
1985 GenerateBiasPowEnergies(); <<
1986 } <<
1987 else <<
1988 { <<
1989 if (EnergyDisType == "Mono") <<
1990 { <<
1991 GenerateMonoEnergetic(); <<
1992 } <<
1993 else if (EnergyDisType == "Lin") <<
1994 { <<
1995 GenerateLinearEnergies(false); <<
1996 } <<
1997 else if (EnergyDisType == "Pow") <<
1998 { <<
1999 GeneratePowEnergies(false); <<
2000 } <<
2001 else if (EnergyDisType == "CPow") <<
2002 { <<
2003 GenerateCPowEnergies(); <<
2004 } <<
2005 else if (EnergyDisType == "Exp") <<
2006 { <<
2007 GenerateExpEnergies(false); <<
2008 } <<
2009 else if (EnergyDisType == "Gauss") <<
2010 { <<
2011 GenerateGaussEnergies(); <<
2012 } <<
2013 else if (EnergyDisType == "Brem") <<
2014 { <<
2015 GenerateBremEnergies(); <<
2016 } <<
2017 else if (EnergyDisType == "Bbody") <<
2018 { <<
2019 GenerateBbodyEnergies(); <<
2020 } <<
2021 else if (EnergyDisType == "Cdg") <<
2022 { <<
2023 GenerateCdgEnergies(); <<
2024 } <<
2025 else if (EnergyDisType == "User") <<
2026 { <<
2027 GenUserHistEnergies(); <<
2028 } <<
2029 else if (EnergyDisType == "Arb") <<
2030 { <<
2031 GenArbPointEnergies(); <<
2032 } <<
2033 else if (EnergyDisType == "Epn") <<
2034 { <<
2035 GenEpnHistEnergies(); <<
2036 } <<
2037 else <<
2038 { <<
2039 G4cout << "Error: EnergyDisType has u <<
2040 } <<
2041 } <<
2042 } <<
2043 return params.particle_energy; <<
2044 } <<
2045 <<
2046 G4double G4SPSEneDistribution::GetProbability <<
2047 { <<
2048 G4double prob = 1.; <<
2049 <<
2050 threadLocal_t& params = threadLocalData.Get <<
2051 if (EnergyDisType == "Lin") <<
2052 { <<
2053 if (prob_norm == 1.) <<
2054 { <<
2055 prob_norm = 0.5*params.grad*params.Emax <<
2056 + params.cept*params.Emax <<
2057 - 0.5*params.grad*params.Emin <<
2058 - params.cept*params.Emin; <<
2059 } <<
2060 prob = params.cept + params.grad * ene; <<
2061 prob /= prob_norm; <<
2062 } <<
2063 else if (EnergyDisType == "Pow") <<
2064 { <<
2065 if (prob_norm == 1.) <<
2066 { <<
2067 if (alpha != -1.) <<
2068 { <<
2069 G4double emina = std::pow(params.Emin <<
2070 G4double emaxa = std::pow(params.Emax <<
2071 prob_norm = 1./(1.+alpha) * (emaxa - <<
2072 } <<
2073 else <<
2074 { <<
2075 prob_norm = std::log(params.Emax) - s <<
2076 } <<
2077 } <<
2078 prob = std::pow(ene, params.alpha)/prob_n <<
2079 } <<
2080 else if (EnergyDisType == "Exp") <<
2081 { <<
2082 if (prob_norm == 1.) <<
2083 { <<
2084 prob_norm = -params.Ezero*(std::exp(-pa <<
2085 - std::exp(par <<
2086 } <<
2087 prob = std::exp(-ene / params.Ezero); <<
2088 prob /= prob_norm; <<
2089 } <<
2090 else if (EnergyDisType == "Arb") <<
2091 { <<
2092 prob = ArbEnergyH.Value(ene); <<
2093 <<
2094 if (prob <= 0.) <<
2095 { <<
2096 G4cout << " Warning:G4SPSEneDistributio <<
2097 << prob << " " << ene << G4endl; <<
2098 prob = 1e-30; <<
2099 } <<
2100 } <<
2101 else <<
2102 { <<
2103 G4cout << "Error: EnergyDisType not suppo <<
2104 } <<
2105 <<
2106 return prob; <<
2107 } 1301 }
2108 1302