| 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 //
>> 28 // MODULE: G4SPSEneDistribution.cc
>> 29 //
>> 30 // Version: 1.0
>> 31 // Date: 5/02/04
>> 32 // Author: Fan Lei
>> 33 // Organisation: QinetiQ ltd.
>> 34 // Customer: ESA/ESTEC
>> 35 //
>> 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 ///////////////////////////////////////////////////////////////////////////////
27 // 46 //
28 // Author: Fan Lei, QinetiQ ltd - 05/02/2004 <<
29 // Customer: ESA/ESTEC <<
30 // Revisions: Andrew Green, Andrea Dotti <<
31 // ------------------------------------------- <<
32 #include "G4SPSEneDistribution.hh" <<
33 <<
34 #include "G4Exp.hh" <<
35 #include "G4SystemOfUnits.hh" <<
36 #include "G4UnitsTable.hh" <<
37 #include "Randomize.hh" 47 #include "Randomize.hh"
38 #include "G4AutoLock.hh" << 48 //#include <cmath>
39 #include "G4Threading.hh" << 49
>> 50 #include "G4SPSEneDistribution.hh"
40 51
41 G4SPSEneDistribution::G4SPSEneDistribution() 52 G4SPSEneDistribution::G4SPSEneDistribution()
42 { 53 {
43 G4MUTEXINIT(mutex); << 54 //
44 <<
45 // Initialise all variables 55 // Initialise all variables
>> 56 particle_energy = 1.0*MeV;
46 57
47 particle_energy = 1.0 * MeV; <<
48 EnergyDisType = "Mono"; 58 EnergyDisType = "Mono";
49 weight=1.; << 59 MonoEnergy = 1*MeV;
50 MonoEnergy = 1 * MeV; <<
51 Emin = 0.; 60 Emin = 0.;
52 Emax = 1.e30; 61 Emax = 1.e30;
53 alpha = 0.; 62 alpha = 0.;
54 biasalpha = 0.; <<
55 prob_norm = 1.0; <<
56 Ezero = 0.; 63 Ezero = 0.;
57 SE = 0.; 64 SE = 0.;
58 Temp = 0.; 65 Temp = 0.;
59 grad = 0.; 66 grad = 0.;
60 cept = 0.; 67 cept = 0.;
>> 68 EnergySpec = true; // true - energy spectra, false - momentum spectra
>> 69 DiffSpec = true; // true - differential spec, false integral spec
61 IntType = "NULL"; // Interpolation type 70 IntType = "NULL"; // Interpolation type
>> 71 IPDFEnergyExist = false;
>> 72 IPDFArbExist = false;
62 73
63 ArbEmin = 0.; 74 ArbEmin = 0.;
64 ArbEmax = 1.e30; 75 ArbEmax = 1.e30;
65 76
66 verbosityLevel = 0; << 77 verbosityLevel = 0 ;
67 << 78
68 threadLocal_t& data = threadLocalData.Get(); <<
69 data.Emax = Emax; <<
70 data.Emin = Emin; <<
71 data.alpha =alpha; <<
72 data.cept = cept; <<
73 data.Ezero = Ezero; <<
74 data.grad = grad; <<
75 data.particle_energy = 0.; <<
76 data.particle_definition = nullptr; <<
77 data.weight = weight; <<
78 } 79 }
79 80
80 G4SPSEneDistribution::~G4SPSEneDistribution() 81 G4SPSEneDistribution::~G4SPSEneDistribution()
81 { << 82 {}
82 G4MUTEXDESTROY(mutex); <<
83 if(Arb_grad_cept_flag) <<
84 { <<
85 delete [] Arb_grad; <<
86 delete [] Arb_cept; <<
87 } <<
88 <<
89 if(Arb_alpha_Const_flag) <<
90 { <<
91 delete [] Arb_alpha; <<
92 delete [] Arb_Const; <<
93 } <<
94 <<
95 if(Arb_ezero_flag) <<
96 { <<
97 delete [] Arb_ezero; <<
98 } <<
99 delete Bbody_x; <<
100 delete BBHist; <<
101 delete CP_x; <<
102 delete CPHist; <<
103 for (auto & it : SplineInt) <<
104 { <<
105 delete it; <<
106 it = nullptr; <<
107 } <<
108 SplineInt.clear(); <<
109 } <<
110 83
111 void G4SPSEneDistribution::SetEnergyDisType(co << 84 void G4SPSEneDistribution::SetEnergyDisType(G4String DisType)
112 { 85 {
113 G4AutoLock l(&mutex); <<
114 EnergyDisType = DisType; 86 EnergyDisType = DisType;
115 if (EnergyDisType == "User") << 87 if (EnergyDisType == "User"){
116 { << 88 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ;
117 UDefEnergyH = IPDFEnergyH = ZeroPhysVector << 89 IPDFEnergyExist = false ;
118 IPDFEnergyExist = false; << 90 } else if ( EnergyDisType == "Arb"){
119 } << 91 ArbEnergyH =IPDFArbEnergyH = ZeroPhysVector ;
120 else if (EnergyDisType == "Arb") <<
121 { <<
122 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVect <<
123 IPDFArbExist = false; 92 IPDFArbExist = false;
124 } << 93 } else if (EnergyDisType == "Epn"){
125 else if (EnergyDisType == "Epn") << 94 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ;
126 { << 95 IPDFEnergyExist = false ;
127 UDefEnergyH = IPDFEnergyH = ZeroPhysVector << 96 EpnEnergyH = ZeroPhysVector ;
128 IPDFEnergyExist = false; <<
129 EpnEnergyH = ZeroPhysVector; <<
130 } 97 }
131 } 98 }
132 99
133 const G4String& G4SPSEneDistribution::GetEnerg <<
134 { <<
135 G4AutoLock l(&mutex); <<
136 return EnergyDisType; <<
137 } <<
138 <<
139 void G4SPSEneDistribution::SetEmin(G4double em 100 void G4SPSEneDistribution::SetEmin(G4double emi)
140 { 101 {
141 G4AutoLock l(&mutex); <<
142 Emin = emi; 102 Emin = emi;
143 threadLocalData.Get().Emin = Emin; <<
144 } <<
145 <<
146 G4double G4SPSEneDistribution::GetEmin() const <<
147 { <<
148 return threadLocalData.Get().Emin; <<
149 } <<
150 <<
151 G4double G4SPSEneDistribution::GetArbEmin() <<
152 { <<
153 G4AutoLock l(&mutex); <<
154 return ArbEmin; <<
155 } <<
156 <<
157 G4double G4SPSEneDistribution::GetArbEmax() <<
158 { <<
159 G4AutoLock l(&mutex); <<
160 return ArbEmax; <<
161 } 103 }
162 104
163 void G4SPSEneDistribution::SetEmax(G4double em 105 void G4SPSEneDistribution::SetEmax(G4double ema)
164 { 106 {
165 G4AutoLock l(&mutex); <<
166 Emax = ema; 107 Emax = ema;
167 threadLocalData.Get().Emax = Emax; <<
168 } <<
169 <<
170 G4double G4SPSEneDistribution::GetEmax() const <<
171 { <<
172 return threadLocalData.Get().Emax; <<
173 } 108 }
174 109
175 void G4SPSEneDistribution::SetMonoEnergy(G4dou 110 void G4SPSEneDistribution::SetMonoEnergy(G4double menergy)
176 { 111 {
177 G4AutoLock l(&mutex); <<
178 MonoEnergy = menergy; 112 MonoEnergy = menergy;
179 } 113 }
180 114
181 void G4SPSEneDistribution::SetBeamSigmaInE(G4d 115 void G4SPSEneDistribution::SetBeamSigmaInE(G4double e)
182 { 116 {
183 G4AutoLock l(&mutex); <<
184 SE = e; 117 SE = e;
185 } 118 }
186 void G4SPSEneDistribution::SetAlpha(G4double a 119 void G4SPSEneDistribution::SetAlpha(G4double alp)
187 { 120 {
188 G4AutoLock l(&mutex); <<
189 alpha = alp; 121 alpha = alp;
190 threadLocalData.Get().alpha = alpha; <<
191 } <<
192 <<
193 void G4SPSEneDistribution::SetBiasAlpha(G4doub <<
194 { <<
195 G4AutoLock l(&mutex); <<
196 biasalpha = alp; <<
197 Biased = true; <<
198 } 122 }
199 123
200 void G4SPSEneDistribution::SetTemp(G4double te 124 void G4SPSEneDistribution::SetTemp(G4double tem)
201 { 125 {
202 G4AutoLock l(&mutex); <<
203 Temp = tem; 126 Temp = tem;
204 } 127 }
205 128
206 void G4SPSEneDistribution::SetEzero(G4double e 129 void G4SPSEneDistribution::SetEzero(G4double eze)
207 { 130 {
208 G4AutoLock l(&mutex); <<
209 Ezero = eze; 131 Ezero = eze;
210 threadLocalData.Get().Ezero = Ezero; <<
211 } 132 }
212 133
213 void G4SPSEneDistribution::SetGradient(G4doubl 134 void G4SPSEneDistribution::SetGradient(G4double gr)
214 { 135 {
215 G4AutoLock l(&mutex); <<
216 grad = gr; 136 grad = gr;
217 threadLocalData.Get().grad = grad; <<
218 } 137 }
219 138
220 void G4SPSEneDistribution::SetInterCept(G4doub 139 void G4SPSEneDistribution::SetInterCept(G4double c)
221 { 140 {
222 G4AutoLock l(&mutex); <<
223 cept = c; 141 cept = c;
224 threadLocalData.Get().cept = cept; <<
225 } <<
226 <<
227 const G4String& G4SPSEneDistribution::GetIntTy <<
228 { <<
229 G4AutoLock l(&mutex); <<
230 return IntType; <<
231 } <<
232 <<
233 void G4SPSEneDistribution::SetBiasRndm(G4SPSRa <<
234 { <<
235 G4AutoLock l(&mutex); <<
236 eneRndm = a; <<
237 } <<
238 <<
239 void G4SPSEneDistribution::SetVerbosity(G4int <<
240 { <<
241 G4AutoLock l(&mutex); <<
242 verbosityLevel = a; <<
243 } <<
244 <<
245 G4double G4SPSEneDistribution::GetWeight() con <<
246 { <<
247 return threadLocalData.Get().weight; <<
248 } <<
249 <<
250 G4double G4SPSEneDistribution::GetMonoEnergy() <<
251 { <<
252 G4AutoLock l(&mutex); <<
253 return MonoEnergy; <<
254 } <<
255 <<
256 G4double G4SPSEneDistribution::GetSE() <<
257 { <<
258 G4AutoLock l(&mutex); <<
259 return SE; <<
260 } <<
261 <<
262 G4double G4SPSEneDistribution::Getalpha() cons <<
263 { <<
264 return threadLocalData.Get().alpha; <<
265 } <<
266 <<
267 G4double G4SPSEneDistribution::GetEzero() cons <<
268 { <<
269 return threadLocalData.Get().Ezero; <<
270 } <<
271 <<
272 G4double G4SPSEneDistribution::GetTemp() <<
273 { <<
274 G4AutoLock l(&mutex); <<
275 return Temp; <<
276 } <<
277 <<
278 G4double G4SPSEneDistribution::Getgrad() const <<
279 { <<
280 return threadLocalData.Get().grad; <<
281 } 142 }
282 143
283 G4double G4SPSEneDistribution::Getcept() const << 144 void G4SPSEneDistribution::UserEnergyHisto(G4ThreeVector input)
284 { 145 {
285 return threadLocalData.Get().cept; << 146 G4double ehi, val;
286 } << 147 ehi = input.x();
287 << 148 val = input.y();
288 G4PhysicsFreeVector G4SPSEneDistribution::GetU << 149 if(verbosityLevel > 1) {
289 { <<
290 G4AutoLock l(&mutex); <<
291 return UDefEnergyH; <<
292 } <<
293 <<
294 G4PhysicsFreeVector G4SPSEneDistribution::GetA <<
295 { <<
296 G4AutoLock l(&mutex); <<
297 return ArbEnergyH; <<
298 } <<
299 <<
300 void G4SPSEneDistribution::UserEnergyHisto(con <<
301 { <<
302 G4AutoLock l(&mutex); <<
303 G4double ehi = input.x(), <<
304 val = input.y(); <<
305 if (verbosityLevel > 1) <<
306 { <<
307 G4cout << "In UserEnergyHisto" << G4endl; 150 G4cout << "In UserEnergyHisto" << G4endl;
308 G4cout << " " << ehi << " " << val << G4en 151 G4cout << " " << ehi << " " << val << G4endl;
309 } 152 }
310 UDefEnergyH.InsertValues(ehi, val); 153 UDefEnergyH.InsertValues(ehi, val);
311 Emax = ehi; 154 Emax = ehi;
312 threadLocalData.Get().Emax = Emax; <<
313 } 155 }
314 156
315 void G4SPSEneDistribution::ArbEnergyHisto(cons << 157 void G4SPSEneDistribution::ArbEnergyHisto(G4ThreeVector input)
316 { 158 {
317 G4AutoLock l(&mutex); << 159 G4double ehi, val;
318 G4double ehi = input.x(), << 160 ehi = input.x();
319 val = input.y(); << 161 val = input.y();
320 if (verbosityLevel > 1) << 162 if(verbosityLevel >1 ) {
321 { <<
322 G4cout << "In ArbEnergyHisto" << G4endl; 163 G4cout << "In ArbEnergyHisto" << G4endl;
323 G4cout << " " << ehi << " " << val << G4en 164 G4cout << " " << ehi << " " << val << G4endl;
324 } 165 }
325 ArbEnergyH.InsertValues(ehi, val); 166 ArbEnergyH.InsertValues(ehi, val);
326 } 167 }
327 168
328 void G4SPSEneDistribution::ArbEnergyHistoFile( << 169 void G4SPSEneDistribution::EpnEnergyHisto(G4ThreeVector input)
329 { 170 {
330 G4AutoLock l(&mutex); <<
331 std::ifstream infile(filename, std::ios::in) <<
332 if (!infile) <<
333 { <<
334 G4Exception("G4SPSEneDistribution::ArbEner <<
335 FatalException, "Unable to ope <<
336 } <<
337 G4double ehi, val; 171 G4double ehi, val;
338 while (infile >> ehi >> val) << 172 ehi = input.x();
339 { << 173 val = input.y();
340 ArbEnergyH.InsertValues(ehi, val); << 174 if(verbosityLevel > 1) {
341 } <<
342 } <<
343 <<
344 void G4SPSEneDistribution::EpnEnergyHisto(cons <<
345 { <<
346 G4AutoLock l(&mutex); <<
347 G4double ehi = input.x(), <<
348 val = input.y(); <<
349 if (verbosityLevel > 1) <<
350 { <<
351 G4cout << "In EpnEnergyHisto" << G4endl; 175 G4cout << "In EpnEnergyHisto" << G4endl;
352 G4cout << " " << ehi << " " << val << G4en 176 G4cout << " " << ehi << " " << val << G4endl;
353 } 177 }
354 EpnEnergyH.InsertValues(ehi, val); 178 EpnEnergyH.InsertValues(ehi, val);
355 Emax = ehi; 179 Emax = ehi;
356 threadLocalData.Get().Emax = Emax; <<
357 Epnflag = true; 180 Epnflag = true;
358 } 181 }
359 182
360 void G4SPSEneDistribution::Calculate() 183 void G4SPSEneDistribution::Calculate()
361 { 184 {
362 G4AutoLock l(&mutex); << 185 if(EnergyDisType == "Cdg")
363 if (EnergyDisType == "Cdg") <<
364 { <<
365 CalculateCdgSpectrum(); 186 CalculateCdgSpectrum();
366 } << 187 else if(EnergyDisType == "Bbody")
367 else if (EnergyDisType == "Bbody") <<
368 { <<
369 if(!BBhistInit) <<
370 { <<
371 BBInitHists(); <<
372 } <<
373 CalculateBbodySpectrum(); 188 CalculateBbodySpectrum();
374 } <<
375 else if (EnergyDisType == "CPow") <<
376 { <<
377 if(!CPhistInit) <<
378 { <<
379 CPInitHists(); <<
380 } <<
381 CalculateCPowSpectrum(); <<
382 } <<
383 } 189 }
384 190
385 void G4SPSEneDistribution::BBInitHists() // M << 191 void G4SPSEneDistribution::CalculateCdgSpectrum()
386 { 192 {
387 BBHist = new std::vector<G4double>(10001, 0. << 193 // This uses the spectrum from The INTEGRAL Mass Model (TIMM)
388 Bbody_x = new std::vector<G4double>(10001, 0 <<
389 BBhistInit = true; <<
390 } <<
391 <<
392 void G4SPSEneDistribution::CPInitHists() // M <<
393 { <<
394 CPHist = new std::vector<G4double>(10001, 0. <<
395 CP_x = new std::vector<G4double>(10001, 0.0) <<
396 CPhistInit = true; <<
397 } <<
398 <<
399 void G4SPSEneDistribution::CalculateCdgSpectru <<
400 { <<
401 // This uses the spectrum from the INTEGRAL <<
402 // to generate a Cosmic Diffuse X/gamma ray 194 // to generate a Cosmic Diffuse X/gamma ray spectrum.
403 << 195 G4double pfact[2] = {8.5, 112};
404 G4double pfact[2] = { 8.5, 112 }; << 196 G4double spind[2] = {1.4, 2.3};
405 G4double spind[2] = { 1.4, 2.3 }; << 197 G4double ene_line[3] = {1.*keV, 18.*keV, 1E6*keV};
406 G4double ene_line[3] = { 1. * keV, 18. * keV <<
407 G4int n_par; 198 G4int n_par;
408 199
409 ene_line[0] = threadLocalData.Get().Emin; << 200 ene_line[0] = Emin;
410 if (threadLocalData.Get().Emin < 18 * keV) << 201 if(Emin < 18*keV)
411 { <<
412 n_par = 2; <<
413 ene_line[2] = threadLocalData.Get().Emax; <<
414 if (threadLocalData.Get().Emax < 18 * keV) <<
415 { 202 {
416 n_par = 1; << 203 n_par = 2;
417 ene_line[1] = threadLocalData.Get().Emax << 204 ene_line[2] = Emax;
>> 205 if(Emax < 18*keV)
>> 206 {
>> 207 n_par = 1;
>> 208 ene_line[1] = Emax;
>> 209 }
418 } 210 }
419 } <<
420 else 211 else
421 { << 212 {
422 n_par = 1; << 213 n_par = 1;
423 pfact[0] = 112.; << 214 pfact[0] = 112.;
424 spind[0] = 2.3; << 215 spind[0] = 2.3;
425 ene_line[1] = threadLocalData.Get().Emax; << 216 ene_line[1] = Emax;
426 } << 217 }
427 << 218
428 // Create a cumulative histogram << 219 // Create a cumulative histogram.
429 // <<
430 CDGhist[0] = 0.; 220 CDGhist[0] = 0.;
431 G4double omalpha; 221 G4double omalpha;
432 G4int i = 0; 222 G4int i = 0;
433 while (i < n_par) <<
434 { <<
435 omalpha = 1. - spind[i]; <<
436 CDGhist[i + 1] = CDGhist[i] + (pfact[i] / <<
437 * (std::pow(en <<
438 - std::pow(ene <<
439 ++i; <<
440 } <<
441 223
442 // Normalise histo and divide by 1000 to mak << 224 while(i < n_par)
443 // << 225 {
>> 226 omalpha = 1. - spind[i];
>> 227 CDGhist[i+1] = CDGhist[i] + (pfact[i]/omalpha)*
>> 228 (std::pow(ene_line[i+1]/keV,omalpha)-std::pow(ene_line[i]/keV,omalpha));
>> 229 i++;
>> 230 }
>> 231
>> 232 // Normalise histo and divide by 1000 to make MeV.
444 i = 0; 233 i = 0;
445 while (i < n_par) << 234 while(i < n_par)
446 { << 235 {
447 CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[ << 236 CDGhist[i+1] = CDGhist[i+1]/CDGhist[n_par];
448 ++i; << 237 // G4cout << CDGhist[i] << CDGhist[n_par] << G4endl;
449 } << 238 i++;
>> 239 }
450 } 240 }
451 241
452 void G4SPSEneDistribution::CalculateBbodySpect << 242 void G4SPSEneDistribution::CalculateBbodySpectrum()
453 { 243 {
454 // Create bbody spectrum << 244 // create bbody spectrum
455 // Proved very hard to integrate indefinitel << 245 // Proved very hard to integrate indefinitely, so different
456 // User inputs emin, emax and T. These are u << 246 // method. User inputs emin, emax and T. These are used to
457 // bin histogram. << 247 // create a 10,000 bin histogram.
458 // Use photon density spectrum = 2 nu**2/c** 248 // Use photon density spectrum = 2 nu**2/c**2 * (std::exp(h nu/kT)-1)
459 // = 2 E**2/h**2c**2 times the exponential 249 // = 2 E**2/h**2c**2 times the exponential
460 << 250 G4double erange = Emax - Emin;
461 G4double erange = threadLocalData.Get().Emax << 251 G4double steps = erange/10000.;
462 G4double steps = erange / 10000.; << 252 G4double Bbody_y[10000];
463 << 253 G4double k = 8.6181e-11; //Boltzmann const in MeV/K
464 const G4double k = 8.6181e-11; //Boltzmann c << 254 G4double h = 4.1362e-21; // Plancks const in MeV s
465 const G4double h = 4.1362e-21; // Plancks co << 255 G4double c = 3e8; // Speed of light
466 const G4double c = 3e8; // Speed of light << 256 G4double h2 = h*h;
467 const G4double h2 = h * h; << 257 G4double c2 = c*c;
468 const G4double c2 = c * c; <<
469 G4int count = 0; <<
470 G4double sum = 0.; <<
471 BBHist->at(0) = 0.; <<
472 <<
473 while (count < 10000) <<
474 { <<
475 Bbody_x->at(count) = threadLocalData.Get() <<
476 G4double Bbody_y = (2. * std::pow(Bbody_x- <<
477 / (h2*c2*(std::exp(Bbody_ <<
478 sum = sum + Bbody_y; <<
479 BBHist->at(count + 1) = BBHist->at(count) <<
480 ++count; <<
481 } <<
482 <<
483 Bbody_x->at(10000) = threadLocalData.Get().E <<
484 <<
485 // Normalise cumulative histo <<
486 // <<
487 count = 0; <<
488 while (count < 10001) <<
489 { <<
490 BBHist->at(count) = BBHist->at(count) / su <<
491 ++count; <<
492 } <<
493 } <<
494 <<
495 void G4SPSEneDistribution::CalculateCPowSpectr <<
496 { <<
497 // Create cutoff power-law spectrum, x^a exp <<
498 // The integral of this function is an incom <<
499 // is only available in the Boost library. <<
500 // <<
501 // User inputs are emin, emax and alpha and <<
502 // create a 10,000 bin histogram. <<
503 <<
504 G4double erange = threadLocalData.Get().Emax <<
505 G4double steps = erange / 10000.; <<
506 alpha = threadLocalData.Get().alpha ; <<
507 Ezero = threadLocalData.Get().Ezero ; <<
508 <<
509 G4int count = 0; 258 G4int count = 0;
510 G4double sum = 0.; 259 G4double sum = 0.;
511 CPHist->at(0) = 0.; << 260 BBHist[0] = 0.;
512 << 261 while(count < 10000)
513 while (count < 10000) << 262 {
514 { << 263 Bbody_x[count] = Emin + G4double(count*steps);
515 CP_x->at(count) = threadLocalData.Get().Em << 264 Bbody_y[count] = (2.*std::pow(Bbody_x[count],2.))/
516 G4double CP_y = std::pow(CP_x->at(count), << 265 (h2*c2*(std::exp(Bbody_x[count]/(k*Temp)) - 1.));
517 * std::exp(-CP_x->at(count) << 266 sum = sum + Bbody_y[count];
518 sum = sum + CP_y; << 267 BBHist[count+1] = BBHist[count] + Bbody_y[count];
519 CPHist->at(count + 1) = CPHist->at(count) << 268 count++;
520 ++count; << 269 }
521 } <<
522 <<
523 CP_x->at(10000) = threadLocalData.Get().Emax <<
524 270
525 // Normalise cumulative histo << 271 Bbody_x[10000] = Emax;
526 // << 272 // Normalise cumulative histo.
527 count = 0; 273 count = 0;
528 while (count < 10001) << 274 while(count<10001)
529 { << 275 {
530 CPHist->at(count) = CPHist->at(count) / su << 276 BBHist[count] = BBHist[count]/sum;
531 ++count; << 277 count++;
532 } << 278 }
533 } 279 }
534 280
535 void G4SPSEneDistribution::InputEnergySpectra( 281 void G4SPSEneDistribution::InputEnergySpectra(G4bool value)
536 { 282 {
537 G4AutoLock l(&mutex); << 283 // Allows user to specifiy spectrum is momentum
538 <<
539 // Allows user to specify spectrum is moment <<
540 // <<
541 EnergySpec = value; // false if momentum 284 EnergySpec = value; // false if momentum
542 if (verbosityLevel > 1) << 285 if(verbosityLevel > 1)
543 { <<
544 G4cout << "EnergySpec has value " << Energ 286 G4cout << "EnergySpec has value " << EnergySpec << G4endl;
545 } <<
546 } 287 }
547 288
548 void G4SPSEneDistribution::InputDifferentialSp 289 void G4SPSEneDistribution::InputDifferentialSpectra(G4bool value)
549 { 290 {
550 G4AutoLock l(&mutex); <<
551 <<
552 // Allows user to specify integral or differ 291 // Allows user to specify integral or differential spectra
553 // <<
554 DiffSpec = value; // true = differential, fa 292 DiffSpec = value; // true = differential, false = integral
555 if (verbosityLevel > 1) << 293 if(verbosityLevel > 1)
556 { <<
557 G4cout << "Diffspec has value " << DiffSpe 294 G4cout << "Diffspec has value " << DiffSpec << G4endl;
558 } <<
559 } 295 }
560 296
561 void G4SPSEneDistribution::ArbInterpolate(cons << 297 void G4SPSEneDistribution::ArbInterpolate(G4String IType)
562 { 298 {
563 G4AutoLock l(&mutex); << 299 if(EnergyDisType != "Arb")
564 << 300 G4cout << "Error: this is for arbitrary distributions" << G4endl;
565 IntType = IType; 301 IntType = IType;
566 ArbEmax = ArbEnergyH.GetMaxEnergy(); << 302 ArbEmax = Emax;
567 ArbEmin = ArbEnergyH.Energy(0); << 303 ArbEmin = Emin;
568 304
569 // Now interpolate points 305 // Now interpolate points
570 << 306 if(IntType == "Lin")
571 if (IntType == "Lin") LinearInterpolation(); << 307 LinearInterpolation();
572 if (IntType == "Log") LogInterpolation(); << 308 if(IntType == "Log")
573 if (IntType == "Exp") ExpInterpolation(); << 309 LogInterpolation();
574 if (IntType == "Spline") SplineInterpolation << 310 if(IntType == "Exp")
>> 311 ExpInterpolation();
>> 312 if(IntType == "Spline")
>> 313 SplineInterpolation();
575 } 314 }
576 315
577 void G4SPSEneDistribution::LinearInterpolation << 316 void G4SPSEneDistribution::LinearInterpolation()
578 { 317 {
579 // Method to do linear interpolation on the 318 // Method to do linear interpolation on the Arb points
580 // Calculate equation of each line segment, 319 // Calculate equation of each line segment, max 1024.
581 // Calculate Area under each segment 320 // Calculate Area under each segment
582 // Create a cumulative array which is then n 321 // Create a cumulative array which is then normalised Arb_Cum_Area
583 322
584 G4double Area_seg[1024]; // Stores area unde 323 G4double Area_seg[1024]; // Stores area under each segment
585 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 324 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
586 std::size_t i, count; << 325 G4int i, count;
587 std::size_t maxi = ArbEnergyH.GetVectorLengt << 326 G4int maxi = ArbEnergyH.GetVectorLength();
588 for (i = 0; i < maxi; ++i) << 327 for(i=0;i<maxi;i++) {
589 { << 328 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
590 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 329 Arb_y[i] = ArbEnergyH(size_t(i));
591 Arb_y[i] = ArbEnergyH(i); <<
592 } 330 }
593 <<
594 // Points are now in x,y arrays. If the spec 331 // Points are now in x,y arrays. If the spectrum is integral it has to be
595 // made differential and if momentum it has << 332 // made differential and if momentum it has to be made energy.
596 << 333 if(DiffSpec == false) {
597 if (!DiffSpec) <<
598 { <<
599 // Converts integral point-wise spectra to 334 // Converts integral point-wise spectra to Differential
600 // << 335 for( count=0;count < maxi-1;count++) {
601 for (count = 0; count < maxi-1; ++count) << 336 Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]);
602 { <<
603 Arb_y[count] = (Arb_y[count] - Arb_y[cou <<
604 / (Arb_x[count + 1] - Arb_x <<
605 } 337 }
606 --maxi; << 338 maxi--;
607 } 339 }
608 << 340 //
609 if (!EnergySpec) << 341 if(EnergySpec == false) {
610 { <<
611 // change currently stored values (emin et 342 // change currently stored values (emin etc) which are actually momenta
612 // to energies << 343 // to energies.
613 // << 344 if(particle_definition == NULL)
614 G4ParticleDefinition* pdef = threadLocalDa << 345 G4cout << "Error: particle not defined" << G4endl;
615 if (pdef == nullptr) << 346 else {
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** 347 // Apply Energy**2 = p**2c**2 + m0**2c**4
624 // p should be entered as E/c i.e. witho 348 // p should be entered as E/c i.e. without the division by c
625 // being done - energy equivalent << 349 // being done - energy equivalent.
626 << 350 G4double mass = particle_definition->GetPDGMass();
627 G4double mass = pdef->GetPDGMass(); << 351 // convert point to energy unit and its value to per energy unit
628 <<
629 // Convert point to energy unit and its <<
630 // <<
631 G4double total_energy; 352 G4double total_energy;
632 for (count = 0; count < maxi; ++count) << 353 for(count=0;count<maxi;count++) {
633 { << 354 total_energy = std::sqrt((Arb_x[count]*Arb_x[count])
634 total_energy = std::sqrt((Arb_x[count] << 355 + (mass*mass)); // total energy
635 + (mass * mass)); // tota << 356
636 Arb_y[count] = Arb_y[count] * Arb_x[co << 357 Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy;
637 Arb_x[count] = total_energy - mass; // << 358 Arb_x[count] = total_energy - mass ; // kinetic energy
638 } 359 }
639 } 360 }
640 } 361 }
641 << 362 //
642 i = 1; << 363 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.; 364 Arb_grad[0] = 0.;
649 Arb_cept[0] = 0.; 365 Arb_cept[0] = 0.;
650 Area_seg[0] = 0.; 366 Area_seg[0] = 0.;
651 Arb_Cum_Area[0] = 0.; 367 Arb_Cum_Area[0] = 0.;
652 while (i < maxi) << 368 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 { 369 {
671 if (verbosityLevel == 2) << 370 // calc gradient and intercept for each segment
672 { << 371 Arb_grad[i] = (Arb_y[i] - Arb_y[i-1]) / (Arb_x[i] - Arb_x[i-1]);
673 G4cout << "Arb_grad is negative" << G4 << 372 if(verbosityLevel == 2)
674 } << 373 G4cout << Arb_grad[i] << G4endl;
675 Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * << 374 if(Arb_grad[i] > 0.)
>> 375 {
>> 376 if(verbosityLevel == 2)
>> 377 G4cout << "Arb_grad is positive" << G4endl;
>> 378 Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]);
>> 379 }
>> 380 else if(Arb_grad[i] < 0.)
>> 381 {
>> 382 if(verbosityLevel == 2)
>> 383 G4cout << "Arb_grad is negative" << G4endl;
>> 384 Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]);
>> 385 }
>> 386 else
>> 387 {
>> 388 if(verbosityLevel == 2)
>> 389 G4cout << "Arb_grad is 0." << G4endl;
>> 390 Arb_cept[i] = Arb_y[i];
>> 391 }
>> 392
>> 393 Area_seg[i] = ((Arb_grad[i]/2)*(Arb_x[i]*Arb_x[i] - Arb_x[i-1]*Arb_x[i-1]) + Arb_cept[i]*(Arb_x[i] - Arb_x[i-1]));
>> 394 Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i];
>> 395 sum = sum + Area_seg[i];
>> 396 if(verbosityLevel == 2)
>> 397 G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i] << G4endl;
>> 398 i++;
676 } 399 }
677 else << 400
>> 401 i=0;
>> 402 while(i < maxi)
678 { 403 {
679 if (verbosityLevel == 2) << 404 Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum; // normalisation
680 { << 405 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
681 G4cout << "Arb_grad is 0." << G4endl; << 406 i++;
682 } <<
683 Arb_cept[i] = Arb_y[i]; <<
684 } 407 }
685 408
686 Area_seg[i] = ((Arb_grad[i] / 2) << 409 if(verbosityLevel >= 1)
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 { 410 {
693 G4cout << Arb_x[i] << Arb_y[i] << Area_s << 411 G4cout << "Leaving LinearInterpolation" << G4endl;
694 << Arb_grad[i] << G4endl; << 412 ArbEnergyH.DumpValues();
>> 413 IPDFArbEnergyH.DumpValues();
695 } 414 }
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 } 415 }
718 416
719 void G4SPSEneDistribution::LogInterpolation() << 417 void G4SPSEneDistribution::LogInterpolation()
720 { 418 {
721 // Interpolation based on Logarithmic equati 419 // Interpolation based on Logarithmic equations
722 // Generate equations of line segments 420 // Generate equations of line segments
723 // y = Ax**alpha => log y = alpha*logx + log 421 // y = Ax**alpha => log y = alpha*logx + logA
724 // Find area under line segments 422 // Find area under line segments
725 // Create normalised, cumulative array Arb_C << 423 // create normalised, cumulative array Arb_Cum_Area
726 <<
727 G4double Area_seg[1024]; // Stores area unde 424 G4double Area_seg[1024]; // Stores area under each segment
728 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 425 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
729 std::size_t i, count; << 426 G4int i, count;
730 std::size_t maxi = ArbEnergyH.GetVectorLengt << 427 G4int maxi = ArbEnergyH.GetVectorLength();
731 for (i = 0; i < maxi; ++i) << 428 for(i=0;i<maxi;i++) {
732 { << 429 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
733 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 430 Arb_y[i] = ArbEnergyH(size_t(i));
734 Arb_y[i] = ArbEnergyH(i); <<
735 } 431 }
736 <<
737 // Points are now in x,y arrays. If the spec 432 // Points are now in x,y arrays. If the spectrum is integral it has to be
738 // made differential and if momentum it has << 433 // made differential and if momentum it has to be made energy.
739 << 434 if(DiffSpec == false) {
740 if (!DiffSpec) <<
741 { <<
742 // Converts integral point-wise spectra to 435 // Converts integral point-wise spectra to Differential
743 // << 436 for( count=0;count<maxi-1;count++) {
744 for (count = 0; count < maxi-1; ++count) << 437 Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]);
745 { <<
746 Arb_y[count] = (Arb_y[count] - Arb_y[cou <<
747 / (Arb_x[count + 1] - Arb_x <<
748 } 438 }
749 --maxi; << 439 maxi--;
750 } 440 }
751 << 441 //
752 if (!EnergySpec) << 442 if(EnergySpec == false) {
753 { << 443 // change currently stored values (emin etc) which are actually momenta
754 // Change currently stored values (emin et << 444 // to energies.
755 // to energies << 445 if(particle_definition == NULL)
756 << 446 G4cout << "Error: particle not defined" << G4endl;
757 G4ParticleDefinition* pdef = threadLocalDa << 447 else {
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** 448 // Apply Energy**2 = p**2c**2 + m0**2c**4
767 // p should be entered as E/c i.e. witho 449 // p should be entered as E/c i.e. without the division by c
768 // being done - energy equivalent << 450 // being done - energy equivalent.
769 << 451 G4double mass = particle_definition->GetPDGMass();
770 G4double mass = pdef->GetPDGMass(); << 452 // convert point to energy unit and its value to per energy unit
771 <<
772 // Convert point to energy unit and its <<
773 // <<
774 G4double total_energy; 453 G4double total_energy;
775 for (count = 0; count < maxi; ++count) << 454 for(count=0;count<maxi;count++) {
776 { << 455 total_energy = std::sqrt((Arb_x[count]*Arb_x[count])
777 total_energy = std::sqrt((Arb_x[count] << 456 + (mass*mass)); // total energy
778 Arb_y[count] = Arb_y[count] * Arb_x[co << 457
779 Arb_x[count] = total_energy - mass; // << 458 Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy;
>> 459 Arb_x[count] = total_energy - mass ; // kinetic energy
780 } 460 }
781 } 461 }
782 } 462 }
783 << 463 //
784 i = 1; << 464 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.; 465 Arb_alpha[0] = 0.;
793 Arb_Const[0] = 0.; 466 Arb_Const[0] = 0.;
794 Area_seg[0] = 0.; 467 Area_seg[0] = 0.;
795 Arb_Cum_Area[0] = 0.; << 468 Arb_Cum_Area[0]=0. ;
796 if (Arb_x[0] <= 0. || Arb_y[0] <= 0.) << 469 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 { 470 {
808 Arb_y[0] = 1e-20; << 471 G4cout << "You should not use log interpolation with points <= 0." << G4endl;
>> 472 G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl;
>> 473 if(Arb_x[0] <= 0.)
>> 474 Arb_x[0] = 1e-20;
>> 475 if(Arb_y[0] <= 0.)
>> 476 Arb_y[0] = 1e-20;
>> 477 }
>> 478
>> 479 G4double alp;
>> 480 while(i <maxi)
>> 481 {
>> 482 // Incase points are negative or zero
>> 483 if(Arb_x[i] <= 0. || Arb_y[i] <= 0.)
>> 484 {
>> 485 G4cout << "You should not use log interpolation with points <= 0." << G4endl;
>> 486 G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl;
>> 487 if(Arb_x[i] <= 0.)
>> 488 Arb_x[i] = 1e-20;
>> 489 if(Arb_y[i] <= 0.)
>> 490 Arb_y[i] = 1e-20;
>> 491 }
>> 492
>> 493 Arb_alpha[i] = (std::log10(Arb_y[i])-std::log10(Arb_y[i-1]))/(std::log10(Arb_x[i])-std::log10(Arb_x[i-1]));
>> 494 Arb_Const[i] = Arb_y[i]/(std::pow(Arb_x[i],Arb_alpha[i]));
>> 495 alp = Arb_alpha[i] + 1;
>> 496 Area_seg[i] = (Arb_Const[i]/alp) * (std::pow(Arb_x[i],alp) - std::pow(Arb_x[i-1],alp));
>> 497 sum = sum + Area_seg[i];
>> 498 Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i];
>> 499 if(verbosityLevel == 2)
>> 500 G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl;
>> 501 i++;
809 } 502 }
810 } << 503
811 << 504 i=0;
812 G4double alp; << 505 while(i<maxi)
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 { 506 {
851 G4cout << Arb_alpha[i] << Arb_Const[i] < << 507 Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum;
>> 508 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
>> 509 i++;
852 } 510 }
853 ++i; << 511 if(verbosityLevel >= 1)
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 512 G4cout << "Leaving LogInterpolation " << G4endl;
871 } <<
872 } 513 }
873 514
874 void G4SPSEneDistribution::ExpInterpolation() << 515 void G4SPSEneDistribution::ExpInterpolation()
875 { 516 {
876 // Interpolation based on Exponential equati 517 // Interpolation based on Exponential equations
877 // Generate equations of line segments 518 // Generate equations of line segments
878 // y = Ae**-(x/e0) => ln y = -x/e0 + lnA 519 // y = Ae**-(x/e0) => ln y = -x/e0 + lnA
879 // Find area under line segments 520 // Find area under line segments
880 // Create normalised, cumulative array Arb_C << 521 // create normalised, cumulative array Arb_Cum_Area
881 <<
882 G4double Area_seg[1024]; // Stores area unde 522 G4double Area_seg[1024]; // Stores area under each segment
883 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 523 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
884 std::size_t i, count; << 524 G4int i, count;
885 std::size_t maxi = ArbEnergyH.GetVectorLengt << 525 G4int maxi = ArbEnergyH.GetVectorLength();
886 for (i = 0; i < maxi; ++i) << 526 for(i=0;i<maxi;i++) {
887 { << 527 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
888 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 528 Arb_y[i] = ArbEnergyH(size_t(i));
889 Arb_y[i] = ArbEnergyH(i); <<
890 } 529 }
891 <<
892 // Points are now in x,y arrays. If the spec 530 // Points are now in x,y arrays. If the spectrum is integral it has to be
893 // made differential and if momentum it has << 531 // made differential and if momentum it has to be made energy.
894 << 532 if(DiffSpec == false) {
895 if (!DiffSpec) <<
896 { <<
897 // Converts integral point-wise spectra to 533 // Converts integral point-wise spectra to Differential
898 // << 534 for( count=0;count< maxi-1;count++) {
899 for (count = 0; count < maxi - 1; ++count) << 535 Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]);
900 { <<
901 Arb_y[count] = (Arb_y[count] - Arb_y[cou <<
902 / (Arb_x[count + 1] - Arb_x <<
903 } 536 }
904 --maxi; << 537 maxi--;
905 } 538 }
906 << 539 //
907 if (!EnergySpec) << 540 if(EnergySpec == false) {
908 { << 541 // change currently stored values (emin etc) which are actually momenta
909 // Change currently stored values (emin et << 542 // to energies.
910 // to energies << 543 if(particle_definition == NULL)
911 // << 544 G4cout << "Error: particle not defined" << G4endl;
912 G4ParticleDefinition* pdef = threadLocalDa << 545 else {
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** 546 // Apply Energy**2 = p**2c**2 + m0**2c**4
922 // p should be entered as E/c i.e. witho 547 // p should be entered as E/c i.e. without the division by c
923 // being done - energy equivalent << 548 // being done - energy equivalent.
924 << 549 G4double mass = particle_definition->GetPDGMass();
925 G4double mass = pdef->GetPDGMass(); << 550 // convert point to energy unit and its value to per energy unit
926 <<
927 // Convert point to energy unit and its <<
928 // <<
929 G4double total_energy; 551 G4double total_energy;
930 for (count = 0; count < maxi; ++count) << 552 for(count=0;count<maxi;count++) {
931 { << 553 total_energy = std::sqrt((Arb_x[count]*Arb_x[count])
932 total_energy = std::sqrt((Arb_x[count] << 554 + (mass*mass)); // total energy
933 Arb_y[count] = Arb_y[count] * Arb_x[co << 555
934 Arb_x[count] = total_energy - mass; // << 556 Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy;
>> 557 Arb_x[count] = total_energy - mass ; // kinetic energy
935 } 558 }
936 } 559 }
937 } 560 }
938 << 561 //
939 i = 1; << 562 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.; 563 Arb_ezero[0] = 0.;
948 Arb_Const[0] = 0.; 564 Arb_Const[0] = 0.;
949 Area_seg[0] = 0.; 565 Area_seg[0] = 0.;
950 Arb_Cum_Area[0] = 0.; 566 Arb_Cum_Area[0] = 0.;
951 while (i < maxi) << 567 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 { 568 {
965 G4Exception("G4SPSEneDistribution::ExpIn << 569 G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i-1]);
966 "Event0302", JustWarning, << 570 if(test > 0. || test < 0.)
967 "Flat line segment: problem, << 571 {
968 G4cout << "Flat line segment: problem" < << 572 Arb_ezero[i] = -(Arb_x[i] - Arb_x[i-1])/(std::log(Arb_y[i]) - std::log(Arb_y[i-1]));
969 Arb_ezero[i] = 0.; << 573 Arb_Const[i] = Arb_y[i]/(std::exp(-Arb_x[i]/Arb_ezero[i]));
970 Arb_Const[i] = 0.; << 574 Area_seg[i]=-(Arb_Const[i]*Arb_ezero[i])*(std::exp(-Arb_x[i]/Arb_ezero[i])
971 Area_seg[i] = 0.; << 575 -std::exp(-Arb_x[i-1]/Arb_ezero[i]));
972 } << 576 }
973 sum = sum + Area_seg[i]; << 577 else
974 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Ar << 578 {
975 if (verbosityLevel == 2) << 579 G4cout << "Flat line segment: problem" << G4endl;
>> 580 Arb_ezero[i] = 0.;
>> 581 Arb_Const[i] = 0.;
>> 582 Area_seg[i] = 0.;
>> 583 }
>> 584 sum = sum + Area_seg[i];
>> 585 Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i];
>> 586 if(verbosityLevel == 2)
>> 587 G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl;
>> 588 i++;
>> 589 }
>> 590
>> 591 i=0;
>> 592 while(i<maxi)
976 { 593 {
977 G4cout << Arb_ezero[i] << Arb_Const[i] < << 594 Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum;
>> 595 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
>> 596 i++;
978 } 597 }
979 ++i; << 598 if(verbosityLevel >= 1)
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 599 G4cout << "Leaving ExpInterpolation " << G4endl;
997 } <<
998 } 600 }
999 601
1000 void G4SPSEneDistribution::SplineInterpolatio << 602 void G4SPSEneDistribution::SplineInterpolation()
1001 { 603 {
1002 // Interpolation using Splines. 604 // Interpolation using Splines.
1003 // Create Normalised arrays, make x 0->1 an << 605 // Create Normalised arrays, make x 0->1 and y hold
1004 // << 606 // the function (Energy)
1005 // Current method based on the above will n << 607 G4double Arb_x[1024], Arb_y[1024];
1006 // New method is implemented below. << 608 G4int i, count;
1007 << 609 G4int maxi = ArbEnergyH.GetVectorLength();
1008 G4double sum, Arb_x[1024]={0.}, Arb_y[1024] << 610 for(i=0;i<maxi;i++) {
1009 std::size_t i, count; << 611 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
1010 std::size_t maxi = ArbEnergyH.GetVectorLeng << 612 Arb_y[i] = ArbEnergyH(size_t(i));
1011 <<
1012 for (i = 0; i < maxi; ++i) <<
1013 { <<
1014 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i) <<
1015 Arb_y[i] = ArbEnergyH(i); <<
1016 } 613 }
1017 <<
1018 // Points are now in x,y arrays. If the spe 614 // Points are now in x,y arrays. If the spectrum is integral it has to be
1019 // made differential and if momentum it has << 615 // made differential and if momentum it has to be made energy.
1020 << 616 if(DiffSpec == false) {
1021 if (!DiffSpec) <<
1022 { <<
1023 // Converts integral point-wise spectra t 617 // Converts integral point-wise spectra to Differential
1024 // << 618 for( count=0;count< maxi-1;count++) {
1025 for (count = 0; count < maxi - 1; ++count << 619 Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]);
1026 { <<
1027 Arb_y[count] = (Arb_y[count] - Arb_y[co <<
1028 / (Arb_x[count + 1] - Arb_ <<
1029 } 620 }
1030 --maxi; << 621 maxi--;
1031 } 622 }
1032 << 623 //
1033 if (!EnergySpec) << 624 if(EnergySpec == false) {
1034 { << 625 // change currently stored values (emin etc) which are actually momenta
1035 // Change currently stored values (emin e << 626 // to energies.
1036 // to energies << 627 if(particle_definition == NULL)
1037 // << 628 G4cout << "Error: particle not defined" << G4endl;
1038 G4ParticleDefinition* pdef = threadLocalD << 629 else {
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* 630 // Apply Energy**2 = p**2c**2 + m0**2c**4
1048 // p should be entered as E/c i.e. with 631 // p should be entered as E/c i.e. without the division by c
1049 // being done - energy equivalent << 632 // being done - energy equivalent.
1050 << 633 G4double mass = particle_definition->GetPDGMass();
1051 G4double mass = pdef->GetPDGMass(); << 634 // convert point to energy unit and its value to per energy unit
1052 <<
1053 // Convert point to energy unit and its <<
1054 // <<
1055 G4double total_energy; 635 G4double total_energy;
1056 for (count = 0; count < maxi; ++count) << 636 for(count=0;count<maxi;count++) {
1057 { << 637 total_energy = std::sqrt((Arb_x[count]*Arb_x[count])
1058 total_energy = std::sqrt((Arb_x[count << 638 + (mass*mass)); // total energy
1059 Arb_y[count] = Arb_y[count] * Arb_x[c << 639
1060 Arb_x[count] = total_energy - mass; / << 640 Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy;
>> 641 Arb_x[count] = total_energy - mass ; // kinetic energy
1061 } 642 }
1062 } 643 }
1063 } 644 }
1064 << 645 //
1065 i = 1; << 646 for(i=1;i<maxi;i++)
1066 Arb_Cum_Area[0] = 0.; << 647 Arb_y[i] += Arb_y[i-1];
1067 sum = 0.; << 648
1068 Splinetemp = new G4DataInterpolation(Arb_x, << 649 for(i=0;i<maxi;i++)
1069 G4double ei[101], prob[101]; << 650 Arb_y[i] /= Arb_y[maxi-1];
1070 for (auto & it : SplineInt) << 651 // now Arb_y is accumulated normalised probabilities
1071 { << 652 /* for(i=0; i<maxi;i++) {
1072 delete it; << 653 if(verbosityLevel >1)
1073 it = nullptr; << 654 G4cout << i <<" "<< Arb_x[i] << " " << Arb_y[i] << G4endl;
1074 } << 655 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_y[i]);
1075 SplineInt.clear(); << 656 }
1076 SplineInt.resize(1024,nullptr); << 657 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(IPDFArbEnergyH.GetVectorLength()-1);
1077 while (i < maxi) << 658 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(0);
1078 { << 659 */
1079 // 100 step per segment for the integrati << 660 // Should now have normalised cumulative probabilities in Arb_y
1080 << 661 // and energy values in Arb_x.
1081 G4double de = (Arb_x[i] - Arb_x[i - 1])/1 << 662 // maxi = maxi + 1;
1082 G4double area = 0.; << 663 // Put y into x and x into y. The spline interpolation will then
1083 << 664 // go through x-axis to find where to interpolate (cum probability)
1084 for (count = 0; count < 101; ++count) << 665 // then generate a y (which will now be energy).
1085 { << 666 SplineInt = new G4DataInterpolation(Arb_y,Arb_x,maxi,1e30,1e30);
1086 ei[count] = Arb_x[i - 1] + de*count ; << 667 if(verbosityLevel >1 )
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 { 668 {
1104 prob[count] = prob[count-1] + prob[coun << 669 G4cout << SplineInt << G4endl;
>> 670 G4cout << SplineInt->LocateArgument(1.0) << G4endl;
1105 } 671 }
1106 << 672 if(verbosityLevel > 0 )
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 " 673 G4cout << "Leaving SplineInterpolation " << G4endl;
1129 } <<
1130 } 674 }
1131 675
1132 void G4SPSEneDistribution::GenerateMonoEnerge 676 void G4SPSEneDistribution::GenerateMonoEnergetic()
1133 { 677 {
1134 // Method to generate MonoEnergetic particl << 678 // Method to generate MonoEnergetic particles.
1135 << 679 particle_energy = MonoEnergy;
1136 threadLocalData.Get().particle_energy = Mon <<
1137 } 680 }
1138 681
1139 void G4SPSEneDistribution::GenerateGaussEnerg 682 void G4SPSEneDistribution::GenerateGaussEnergies()
1140 { 683 {
1141 // Method to generate Gaussian particles << 684 // Method to generate Gaussian particles.
1142 << 685 particle_energy = G4RandGauss::shoot(MonoEnergy,SE);
1143 G4double ene = G4RandGauss::shoot(MonoEnerg << 686 if (particle_energy < 0) particle_energy = 0.;
1144 if (ene < 0) ene = 0.; <<
1145 threadLocalData.Get().particle_energy = ene <<
1146 } 687 }
1147 688
1148 void G4SPSEneDistribution::GenerateLinearEner 689 void G4SPSEneDistribution::GenerateLinearEnergies(G4bool bArb = false)
1149 { 690 {
1150 G4double rndm; 691 G4double rndm;
1151 threadLocal_t& params = threadLocalData.Get << 692 G4double emaxsq = std::pow(Emax,2.); //Emax squared
1152 G4double emaxsq = std::pow(params.Emax, 2.) << 693 G4double eminsq = std::pow(Emin,2.); //Emin squared
1153 G4double eminsq = std::pow(params.Emin, 2.) << 694 G4double intersq = std::pow(cept,2.); //cept squared
1154 G4double intersq = std::pow(params.cept, 2. <<
1155 695
1156 if (bArb) rndm = G4UniformRand(); 696 if (bArb) rndm = G4UniformRand();
1157 else rndm = eneRndm->GenRandEnergy(); << 697 else rndm = eneRndm->GenRandEnergy();
1158 << 698
1159 G4double bracket = ((params.grad / 2.) << 699 G4double bracket = ((grad/2.)*(emaxsq - eminsq) + cept*(Emax-Emin));
1160 * (emaxsq - eminsq) <<
1161 + params.cept * (params.Em <<
1162 bracket = bracket * rndm; 700 bracket = bracket * rndm;
1163 bracket = bracket + (params.grad / 2.) * em << 701 bracket = bracket + (grad/2.)*eminsq + cept*Emin;
1164 <<
1165 // Now have a quad of form m/2 E**2 + cE - 702 // Now have a quad of form m/2 E**2 + cE - bracket = 0
1166 // <<
1167 bracket = -bracket; 703 bracket = -bracket;
1168 << 704 // G4cout << "BRACKET" << bracket << G4endl;
1169 if (params.grad != 0.) << 705 if(grad != 0.)
1170 { <<
1171 G4double sqbrack = (intersq - 4 * (params <<
1172 sqbrack = std::sqrt(sqbrack); <<
1173 G4double root1 = -params.cept + sqbrack; <<
1174 root1 = root1 / (2. * (params.grad / 2.)) <<
1175 <<
1176 G4double root2 = -params.cept - sqbrack; <<
1177 root2 = root2 / (2. * (params.grad / 2.)) <<
1178 <<
1179 if (root1 > params.Emin && root1 < params <<
1180 { 706 {
1181 params.particle_energy = root1; << 707 G4double sqbrack = (intersq - 4*(grad/2.)*(bracket));
>> 708 // G4cout << "SQBRACK" << sqbrack << G4endl;
>> 709 sqbrack = std::sqrt(sqbrack);
>> 710 G4double root1 = -cept + sqbrack;
>> 711 root1 = root1/(2.*(grad/2.));
>> 712
>> 713 G4double root2 = -cept - sqbrack;
>> 714 root2 = root2/(2.*(grad/2.));
>> 715
>> 716 // G4cout << root1 << " roots " << root2 << G4endl;
>> 717
>> 718 if(root1 > Emin && root1 < Emax)
>> 719 particle_energy = root1;
>> 720 if(root2 > Emin && root2 < Emax)
>> 721 particle_energy = root2;
1182 } 722 }
1183 if (root2 > params.Emin && root2 < params << 723 else if(grad == 0.)
1184 { <<
1185 params.particle_energy = root2; <<
1186 } <<
1187 } <<
1188 else if (params.grad == 0.) <<
1189 { <<
1190 // have equation of form cE - bracket =0 724 // have equation of form cE - bracket =0
1191 // << 725 particle_energy = bracket/cept;
1192 params.particle_energy = bracket / params <<
1193 } <<
1194 <<
1195 if (params.particle_energy < 0.) <<
1196 { <<
1197 params.particle_energy = -params.particle <<
1198 } <<
1199 726
1200 if (verbosityLevel >= 1) << 727 if(particle_energy < 0.)
1201 { << 728 particle_energy = -particle_energy;
1202 G4cout << "Energy is " << params.particle << 729
1203 } << 730 if(verbosityLevel >= 1)
>> 731 G4cout << "Energy is " << particle_energy << G4endl;
1204 } 732 }
1205 733
1206 void G4SPSEneDistribution::GeneratePowEnergie 734 void G4SPSEneDistribution::GeneratePowEnergies(G4bool bArb = false)
1207 { 735 {
1208 // Method to generate particle energies dis << 736 // Method to generate particle energies distributed as
>> 737 // a powerlaw
1209 738
1210 G4double rndm; 739 G4double rndm;
1211 G4double emina, emaxa; 740 G4double emina, emaxa;
1212 <<
1213 threadLocal_t& params = threadLocalData.Get <<
1214 <<
1215 emina = std::pow(params.Emin, params.alpha <<
1216 emaxa = std::pow(params.Emax, params.alpha <<
1217 741
1218 if (bArb) rndm = G4UniformRand(); << 742 emina = std::pow(Emin,alpha+1);
1219 else rndm = eneRndm->GenRandEnergy(); << 743 emaxa = std::pow(Emax,alpha+1);
1220 744
1221 if (params.alpha != -1.) << 745 if (bArb) rndm = G4UniformRand();
1222 { << 746 else rndm = eneRndm->GenRandEnergy();
1223 G4double ene = ((rndm * (emaxa - emina)) << 747
1224 ene = std::pow(ene, (1. / (params.alpha + << 748 if(alpha != -1.)
1225 params.particle_energy = ene; <<
1226 } <<
1227 else <<
1228 { <<
1229 G4double ene = (std::log(params.Emin) <<
1230 + rndm * (std::log(params.Em <<
1231 params.particle_energy = std::exp(ene); <<
1232 } <<
1233 if (verbosityLevel >= 1) <<
1234 { <<
1235 G4cout << "Energy is " << params.particle <<
1236 } <<
1237 } <<
1238 <<
1239 void G4SPSEneDistribution::GenerateCPowEnergi <<
1240 { <<
1241 // Method to generate particle energies dis <<
1242 // cutoff power-law distribution <<
1243 // <<
1244 // CP_x holds Energies, and CPHist holds th <<
1245 // binary search to find correct bin then l <<
1246 // Use the earlier defined histogram + Rand <<
1247 // random numbers following the histos dist <<
1248 <<
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 { 749 {
1275 nabove = middle; << 750 particle_energy = ((rndm*(emaxa - emina)) + emina);
>> 751 particle_energy = std::pow(particle_energy,(1./(alpha+1.)));
1276 } 752 }
1277 else << 753 else if(alpha == -1.)
1278 { 754 {
1279 nbelow = middle; << 755 particle_energy = (std::log(Emin) + rndm*(std::log(Emax) - std::log(Emin)));
>> 756 particle_energy = std::exp(particle_energy);
1280 } 757 }
1281 } << 758 if(verbosityLevel >= 1)
1282 << 759 G4cout << "Energy is " << particle_energy << G4endl;
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 } 760 }
1355 761
1356 void G4SPSEneDistribution::GenerateExpEnergie 762 void G4SPSEneDistribution::GenerateExpEnergies(G4bool bArb = false)
1357 { 763 {
1358 // Method to generate particle energies dis 764 // Method to generate particle energies distributed according
1359 // to an exponential curve << 765 // to an exponential curve.
1360 <<
1361 G4double rndm; 766 G4double rndm;
1362 << 767
1363 if (bArb) rndm = G4UniformRand(); 768 if (bArb) rndm = G4UniformRand();
1364 else rndm = eneRndm->GenRandEnergy(); << 769 else rndm = eneRndm->GenRandEnergy();
1365 770
1366 threadLocal_t& params = threadLocalData.Get << 771 particle_energy = -Ezero*(std::log(rndm*(std::exp(-Emax/Ezero) - std::exp(-Emin/Ezero)) +
1367 params.particle_energy = -params.Ezero << 772 std::exp(-Emin/Ezero)));
1368 * (std::log(rndm * ( << 773 if(verbosityLevel >= 1)
1369 << 774 G4cout << "Energy is " << particle_energy << G4endl;
1370 - <<
1371 <<
1372 + std::exp <<
1373 if (verbosityLevel >= 1) <<
1374 { <<
1375 G4cout << "Energy is " << params.particle <<
1376 } <<
1377 } 775 }
1378 776
1379 void G4SPSEneDistribution::GenerateBremEnergi 777 void G4SPSEneDistribution::GenerateBremEnergies()
1380 { 778 {
1381 // Method to generate particle energies dis << 779 // Method to generate particle energies distributed according
1382 // to a Bremstrahlung equation of the form << 780 // to a Bremstrahlung equation of
1383 // I = const*((kT)**1/2)*E*(e**(-E/kT)) << 781 // form I = const*((kT)**1/2)*E*(e**(-E/kT))
1384 << 782
1385 G4double rndm = eneRndm->GenRandEnergy(); << 783 G4double rndm;
>> 784 rndm = eneRndm->GenRandEnergy();
1386 G4double expmax, expmin, k; 785 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 786
1392 threadLocal_t& params = threadLocalData.Get << 787 k = 8.6181e-11; // Boltzmann's const in MeV/K
>> 788 G4double ksq = std::pow(k,2.); // k squared
>> 789 G4double Tsq = std::pow(Temp,2.); // Temp squared
1393 790
1394 expmax = std::exp(-params.Emax / (k * Temp) << 791 expmax = std::exp(-Emax/(k*Temp));
1395 expmin = std::exp(-params.Emin / (k * Temp) << 792 expmin = std::exp(-Emin/(k*Temp));
1396 793
1397 // If either expmax or expmin are zero then 794 // If either expmax or expmin are zero then this will cause problems
1398 // Most probably this will be because T is 795 // Most probably this will be because T is too low or E is too high
1399 796
1400 if (expmax == 0.) << 797 if(expmax == 0.)
1401 { << 798 G4cout << "*****EXPMAX=0. Choose different E's or Temp" << G4endl;
1402 G4Exception("G4SPSEneDistribution::Genera << 799 if(expmin == 0.)
1403 "Event0302", FatalException, << 800 G4cout << "*****EXPMIN=0. Choose different E's or Temp" << G4endl;
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 801
1413 G4double tempvar = rndm * ((-k) * Temp * (p << 802 G4double tempvar = rndm *((-k)*Temp*(Emax*expmax - Emin*expmin) -
1414 - p << 803 (ksq*Tsq*(expmax-expmin)));
1415 - (ksq * Tsq * (exp <<
1416 804
1417 G4double bigc = (tempvar - k * Temp * param << 805 G4double bigc = (tempvar - k*Temp*Emin*expmin - ksq*Tsq*expmin)/(-k*Temp);
1418 - ksq * Tsq * expmin) / (-k <<
1419 806
1420 // This gives an equation of form: Ee(-E/kT 807 // This gives an equation of form: Ee(-E/kT) + kTe(-E/kT) - C =0
1421 // Solve this iteratively, step from Emin t 808 // Solve this iteratively, step from Emin to Emax in 1000 steps
1422 // and take the best solution. 809 // and take the best solution.
1423 810
1424 G4double erange = params.Emax - params.Emin << 811 G4double erange = Emax - Emin;
1425 G4double steps = erange / 1000.; << 812 G4double steps = erange/1000.;
1426 G4int i; 813 G4int i;
1427 G4double etest, diff, err = 100000.; << 814 G4double etest, diff, err;
1428 << 815
1429 for (i = 1; i < 1000; ++i) << 816 err = 100000.;
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 817
1440 if (diff < err) << 818 for(i=1; i<1000; i++)
1441 { 819 {
1442 err = diff; << 820 etest = Emin + (i-1)*steps;
1443 params.particle_energy = etest; << 821
1444 } << 822 diff = etest*(std::exp(-etest/(k*Temp))) + k*Temp*(std::exp(-etest/(k*Temp))) - bigc;
1445 } << 823
1446 if (verbosityLevel >= 1) << 824 if(diff < 0.)
1447 { << 825 diff = -diff;
1448 G4cout << "Energy is " << params.particle << 826
1449 } << 827 if(diff < err)
>> 828 {
>> 829 err = diff;
>> 830 particle_energy = etest;
>> 831 }
>> 832 }
>> 833 if(verbosityLevel >= 1)
>> 834 G4cout << "Energy is " << particle_energy << G4endl;
1450 } 835 }
1451 836
1452 void G4SPSEneDistribution::GenerateBbodyEnerg 837 void G4SPSEneDistribution::GenerateBbodyEnergies()
1453 { 838 {
1454 // BBody_x holds Energies, and BBHist holds 839 // BBody_x holds Energies, and BBHist holds the cumulative histo.
1455 // Binary search to find correct bin then l << 840 // binary search to find correct bin then lin interpolation.
1456 // Use the earlier defined histogram + Rand 841 // Use the earlier defined histogram + RandGeneral method to generate
1457 // random numbers following the histos dist << 842 // random numbers following the histos distribution.
1458 << 843 G4double rndm;
1459 G4double rndm = eneRndm->GenRandEnergy(); << 844 G4int nabove, nbelow = 0, middle;
1460 G4int nabove = 10001, nbelow = 0, middle; << 845 nabove = 10001;
1461 << 846 rndm = eneRndm->GenRandEnergy();
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 847
1474 // Binary search to find bin that rndm is i 848 // Binary search to find bin that rndm is in
1475 // << 849 while(nabove-nbelow > 1)
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 { 850 {
1485 nabove = middle; << 851 middle = (nabove + nbelow)/2;
>> 852 if(rndm == BBHist[middle]) break;
>> 853 if(rndm < BBHist[middle]) nabove = middle;
>> 854 else nbelow = middle;
1486 } 855 }
1487 else << 856
>> 857 // Now interpolate in that bin to find the correct output value.
>> 858 G4double x1, x2, y1, y2, m, q;
>> 859 x1 = Bbody_x[nbelow];
>> 860 x2 = Bbody_x[nbelow+1];
>> 861 y1 = BBHist[nbelow];
>> 862 y2 = BBHist[nbelow+1];
>> 863 m = (y2-y1)/(x2-x1);
>> 864 q = y1 - m*x1;
>> 865
>> 866 particle_energy = (rndm - q)/m;
>> 867
>> 868 if(verbosityLevel >= 1)
1488 { 869 {
1489 nbelow = middle; << 870 G4cout << "Energy is " << particle_energy << G4endl;
1490 } 871 }
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 } 872 }
1526 873
1527 void G4SPSEneDistribution::GenerateCdgEnergie 874 void G4SPSEneDistribution::GenerateCdgEnergies()
1528 { 875 {
1529 // Generate random numbers, compare with va << 876 // Gen random numbers, compare with values in cumhist
1530 // to find appropriate part of spectrum and << 877 // to find appropriate part of spectrum and then
1531 // generate energy in the usual inversion w << 878 // generate energy in the usual inversion way.
1532 << 879 // G4double pfact[2] = {8.5, 112};
>> 880 // G4double spind[2] = {1.4, 2.3};
>> 881 // G4double ene_line[3] = {1., 18., 1E6};
1533 G4double rndm, rndm2; 882 G4double rndm, rndm2;
1534 G4double ene_line[3]={0,0,0}; << 883 G4double ene_line[3];
1535 G4double omalpha[2]={0,0}; << 884 G4double omalpha[2];
1536 threadLocal_t& params = threadLocalData.Get << 885 if(Emin < 18*keV && Emax < 18*keV)
1537 if (params.Emin < 18 * keV && params.Emax < << 886 {
1538 { << 887 omalpha[0] = 1. - 1.4;
1539 omalpha[0] = 1. - 1.4; << 888 ene_line[0] = Emin;
1540 ene_line[0] = params.Emin; << 889 ene_line[1] = Emax;
1541 ene_line[1] = params.Emax; << 890 }
1542 } << 891 if(Emin < 18*keV && Emax > 18*keV)
1543 if (params.Emin < 18 * keV && params.Emax > << 892 {
1544 { << 893 omalpha[0] = 1. - 1.4;
1545 omalpha[0] = 1. - 1.4; << 894 omalpha[1] = 1. - 2.3;
1546 omalpha[1] = 1. - 2.3; << 895 ene_line[0] = Emin;
1547 ene_line[0] = params.Emin; << 896 ene_line[1] = 18.*keV;
1548 ene_line[1] = 18. * keV; << 897 ene_line[2] = Emax;
1549 ene_line[2] = params.Emax; << 898 }
1550 } << 899 if(Emin > 18*keV)
1551 if (params.Emin > 18 * keV) << 900 {
1552 { << 901 omalpha[0] = 1. - 2.3;
1553 omalpha[0] = 1. - 2.3; << 902 ene_line[0] = Emin;
1554 ene_line[0] = params.Emin; << 903 ene_line[1] = Emax;
1555 ene_line[1] = params.Emax; << 904 }
1556 } <<
1557 rndm = eneRndm->GenRandEnergy(); 905 rndm = eneRndm->GenRandEnergy();
1558 rndm2 = eneRndm->GenRandEnergy(); 906 rndm2 = eneRndm->GenRandEnergy();
1559 907
1560 G4int i = 0; 908 G4int i = 0;
1561 while (rndm >= CDGhist[i]) << 909 while( rndm >= CDGhist[i])
1562 { << 910 {
1563 ++i; << 911 i++;
1564 } << 912 }
>> 913 // Generate final energy.
>> 914 particle_energy = (std::pow(ene_line[i-1],omalpha[i-1]) + (std::pow(ene_line[i],omalpha[i-1])
>> 915 - std::pow(ene_line[i-1],omalpha[i-1]))*rndm2);
>> 916 particle_energy = std::pow(particle_energy,(1./omalpha[i-1]));
1565 917
1566 // Generate final energy << 918 if(verbosityLevel >= 1)
1567 // << 919 G4cout << "Energy is " << particle_energy << G4endl;
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 } 920 }
1578 921
1579 void G4SPSEneDistribution::GenUserHistEnergie 922 void G4SPSEneDistribution::GenUserHistEnergies()
1580 { 923 {
1581 // Histograms are DIFFERENTIAL << 924 // Histograms are DIFFERENTIAL.
1582 << 925 // G4cout << "In GenUserHistEnergies " << G4endl;
1583 G4AutoLock l(&mutex); << 926 if(IPDFEnergyExist == false)
1584 << 927 {
1585 if (!IPDFEnergyExist) << 928 G4int ii;
1586 { << 929 G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
1587 std::size_t ii; << 930 G4double bins[1024], vals[1024], sum;
1588 std::size_t maxbin = UDefEnergyH.GetVecto << 931 sum=0.;
1589 G4double bins[1024], vals[1024], sum; << 932
1590 for ( ii = 0 ; ii<1024 ; ++ii ) { bins[ii << 933 if((EnergySpec == false) && (particle_definition == NULL))
1591 sum = 0.; << 934 G4cout << "Error: particle definition is NULL" << G4endl;
1592 << 935
1593 if ( (!EnergySpec) << 936 if(maxbin > 1024)
1594 && (threadLocalData.Get().particle_defi << 937 {
1595 { << 938 G4cout << "Maxbin > 1024" << G4endl;
1596 G4Exception("G4SPSEneDistribution::GenU << 939 G4cout << "Setting maxbin to 1024, other bins are lost" << G4endl;
1597 "Event0302", FatalException << 940 }
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 941
1660 IPDFEnergyExist = true; << 942 if(DiffSpec == false)
1661 if (verbosityLevel > 1) << 943 G4cout << "Histograms are Differential!!! " << G4endl;
1662 { << 944 else
1663 IPDFEnergyH.DumpValues(); << 945 {
>> 946 bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
>> 947 vals[0] = UDefEnergyH(size_t(0));
>> 948 sum = vals[0];
>> 949 for(ii=1;ii<maxbin;ii++)
>> 950 {
>> 951 bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
>> 952 vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii-1];
>> 953 sum = sum + UDefEnergyH(size_t(ii));
>> 954 }
>> 955 }
>> 956
>> 957 if(EnergySpec == false)
>> 958 {
>> 959 G4double mass = particle_definition->GetPDGMass();
>> 960 // multiply the function (vals) up by the bin width
>> 961 // to make the function counts/s (i.e. get rid of momentum
>> 962 // dependence).
>> 963 for(ii=1;ii<maxbin;ii++)
>> 964 {
>> 965 vals[ii] = vals[ii] * (bins[ii] - bins[ii-1]);
>> 966 }
>> 967 // Put energy bins into new histo, plus divide by energy bin width
>> 968 // to make evals counts/s/energy
>> 969 for(ii=0;ii<maxbin;ii++)
>> 970 {
>> 971 bins[ii] = std::sqrt((bins[ii]*bins[ii]) + (mass*mass)) - mass; //kinetic energy
>> 972 }
>> 973 for(ii=1;ii<maxbin;ii++)
>> 974 {
>> 975 vals[ii] = vals[ii]/(bins[ii] - bins[ii-1]);
>> 976 }
>> 977 sum = vals[maxbin-1];
>> 978 vals[0] = 0.;
>> 979 }
>> 980 for(ii=0;ii<maxbin;ii++)
>> 981 {
>> 982 vals[ii] = vals[ii]/sum;
>> 983 IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
>> 984 }
>> 985
>> 986 // Make IPDFEnergyExist = true
>> 987 IPDFEnergyExist = true;
>> 988 if(verbosityLevel > 1)
>> 989 IPDFEnergyH.DumpValues();
1664 } 990 }
1665 } << 991
1666 l.unlock(); <<
1667 <<
1668 // IPDF has been create so carry on 992 // IPDF has been create so carry on
1669 // <<
1670 G4double rndm = eneRndm->GenRandEnergy(); 993 G4double rndm = eneRndm->GenRandEnergy();
1671 threadLocalData.Get().particle_energy= IPDF << 994 particle_energy = IPDFEnergyH.GetEnergy(rndm);
1672 << 995
1673 if (verbosityLevel >= 1) << 996 if(verbosityLevel >= 1)
1674 { <<
1675 G4cout << "Energy is " << particle_energy 997 G4cout << "Energy is " << particle_energy << G4endl;
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 } 998 }
1708 999
1709 void G4SPSEneDistribution::GenArbPointEnergie 1000 void G4SPSEneDistribution::GenArbPointEnergies()
1710 { 1001 {
1711 if (verbosityLevel > 0) << 1002 if(verbosityLevel > 0)
1712 { <<
1713 G4cout << "In GenArbPointEnergies" << G4e 1003 G4cout << "In GenArbPointEnergies" << G4endl;
1714 } << 1004 G4double rndm;
1715 << 1005 rndm = eneRndm->GenRandEnergy();
1716 G4double rndm = eneRndm->GenRandEnergy(); << 1006 if(IntType != "Spline")
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 { 1007 {
1780 G4cout << "Energy is " << params.partic << 1008 // IPDFArbEnergyH.DumpValues();
>> 1009 // Find the Bin
>> 1010 // have x, y, no of points, and cumulative area distribution
>> 1011 G4int nabove, nbelow = 0, middle;
>> 1012 nabove = IPDFArbEnergyH.GetVectorLength();
>> 1013 // G4cout << nabove << G4endl;
>> 1014 // Binary search to find bin that rndm is in
>> 1015 while(nabove-nbelow > 1)
>> 1016 {
>> 1017 middle = (nabove + nbelow)/2;
>> 1018 if(rndm == IPDFArbEnergyH(size_t(middle))) break;
>> 1019 if(rndm < IPDFArbEnergyH(size_t(middle))) nabove = middle;
>> 1020 else nbelow = middle;
>> 1021 }
>> 1022 if(IntType == "Lin")
>> 1023 {
>> 1024 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1));
>> 1025 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
>> 1026 grad = Arb_grad[nbelow+1];
>> 1027 cept = Arb_cept[nbelow+1];
>> 1028 // G4cout << rndm << " " << Emax << " " << Emin << " " << grad << " " << cept << G4endl;
>> 1029 GenerateLinearEnergies(true);
>> 1030 }
>> 1031 else if(IntType == "Log")
>> 1032 {
>> 1033 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1));
>> 1034 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
>> 1035 alpha = Arb_alpha[nbelow+1];
>> 1036 // G4cout << rndm << " " << Emax << " " << Emin << " " << alpha << G4endl;
>> 1037 GeneratePowEnergies(true);
>> 1038 }
>> 1039 else if(IntType == "Exp")
>> 1040 {
>> 1041 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1));
>> 1042 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
>> 1043 Ezero = Arb_ezero[nbelow+1];
>> 1044 // G4cout << rndm << " " << Emax << " " << Emin << " " << Ezero << G4endl;
>> 1045 GenerateExpEnergies(true);
>> 1046 }
>> 1047 }
>> 1048 else if(IntType == "Spline")
>> 1049 {
>> 1050 if(verbosityLevel > 1)
>> 1051 G4cout << "IntType = Spline " << rndm << G4endl;
>> 1052 // in SplineInterpolation created SplineInt
>> 1053 // Now generate a random number put it into CubicSplineInterpolation
>> 1054 // and you should get out an energy!?!
>> 1055 particle_energy = -1e100;
>> 1056 while (particle_energy < Emin || particle_energy > Emax ) {
>> 1057 particle_energy = SplineInt->CubicSplineInterpolation(rndm);
>> 1058 rndm = eneRndm->GenRandEnergy();
>> 1059 }
>> 1060 if(verbosityLevel >= 1)
>> 1061 G4cout << "Energy is " << particle_energy << G4endl;
1781 } 1062 }
1782 } <<
1783 else 1063 else
1784 { << 1064 G4cout << "Error: IntType unknown type" << G4endl;
1785 G4Exception("G4SPSEneDistribution::GenArb <<
1786 FatalException, "Error: IntTy <<
1787 } <<
1788 } 1065 }
1789 1066
1790 void G4SPSEneDistribution::GenEpnHistEnergies 1067 void G4SPSEneDistribution::GenEpnHistEnergies()
1791 { 1068 {
1792 // Firstly convert to energy if not already << 1069 // G4cout << "In GenEpnHistEnergies " << Epnflag << G4endl;
1793 << 1070
1794 G4AutoLock l(&mutex); << 1071 // Firstly convert to energy if not already done.
1795 << 1072 if(Epnflag == true)
1796 if (Epnflag) // true means spectrum is epn << 1073 // epnflag = true means spectrum is epn, false means e.
1797 { << 1074 {
1798 // Convert to energy by multiplying by A << 1075 // convert to energy by multiplying by A number
1799 // << 1076 ConvertEPNToEnergy();
1800 ConvertEPNToEnergy(); << 1077 // EpnEnergyH will be replace by UDefEnergyH.
1801 } << 1078 // UDefEnergyH.DumpValues();
1802 if (!IPDFEnergyExist) << 1079 }
1803 { << 1080
1804 // IPDF has not been created, so create i << 1081 // G4cout << "Creating IPDFEnergy if not already done so" << G4endl;
1805 // << 1082 if(IPDFEnergyExist == false)
1806 G4double bins[1024], vals[1024], sum; << 1083 {
1807 std::size_t ii; << 1084 // IPDF has not been created, so create it
1808 std::size_t maxbin = UDefEnergyH.GetVecto << 1085 G4double bins[1024],vals[1024], sum;
1809 bins[0] = UDefEnergyH.GetLowEdgeEnergy(0) << 1086 G4int ii;
1810 vals[0] = UDefEnergyH(0); << 1087 G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
1811 sum = vals[0]; << 1088 bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
1812 for (ii = 1; ii < maxbin; ++ii) << 1089 vals[0] = UDefEnergyH(size_t(0));
1813 { << 1090 sum = vals[0];
1814 bins[ii] = UDefEnergyH.GetLowEdgeEnergy << 1091 for(ii=1;ii<maxbin;ii++)
1815 vals[ii] = UDefEnergyH(ii) + vals[ii - << 1092 {
1816 sum = sum + UDefEnergyH(ii); << 1093 bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
1817 } << 1094 vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii-1];
1818 << 1095 sum = sum + UDefEnergyH(size_t(ii));
1819 l.lock(); << 1096 }
1820 for (ii = 0; ii < maxbin; ++ii) << 1097
1821 { << 1098 for(ii=0;ii<maxbin;ii++)
1822 vals[ii] = vals[ii] / sum; << 1099 {
1823 IPDFEnergyH.InsertValues(bins[ii], vals << 1100 vals[ii] = vals[ii]/sum;
>> 1101 IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
>> 1102 }
>> 1103 // Make IPDFEpnExist = true
>> 1104 IPDFEnergyExist = true;
1824 } 1105 }
1825 IPDFEnergyExist = true; << 1106 // IPDFEnergyH.DumpValues();
1826 <<
1827 } <<
1828 l.unlock(); <<
1829 <<
1830 // IPDF has been create so carry on 1107 // IPDF has been create so carry on
1831 // <<
1832 G4double rndm = eneRndm->GenRandEnergy(); 1108 G4double rndm = eneRndm->GenRandEnergy();
1833 threadLocalData.Get().particle_energy = IPD << 1109 particle_energy = IPDFEnergyH.GetEnergy(rndm);
1834 1110
1835 if (verbosityLevel >= 1) << 1111 if(verbosityLevel >= 1)
1836 { << 1112 G4cout << "Energy is " << particle_energy << G4endl;
1837 G4cout << "Energy is " << threadLocalData <<
1838 } <<
1839 } 1113 }
1840 1114
1841 void G4SPSEneDistribution::ConvertEPNToEnergy << 1115 void G4SPSEneDistribution::ConvertEPNToEnergy()
1842 { 1116 {
1843 // Use this before particle generation to c << 1117 // Use this before particle generation to convert the
1844 // currently stored histogram from energy/n << 1118 // currently stored histogram from energy/nucleon
1845 << 1119 // to energy.
1846 threadLocal_t& params = threadLocalData.Get << 1120 // G4cout << "In ConvertEpntoEnergy " << G4endl;
1847 if (params.particle_definition == nullptr) << 1121 if(particle_definition==NULL)
1848 { <<
1849 G4cout << "Error: particle not defined" < 1122 G4cout << "Error: particle not defined" << G4endl;
1850 } <<
1851 else 1123 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 { 1124 {
1873 G4Exception("G4SPSEneDistribution::Conv << 1125 // Need to multiply histogram by the number of nucleons.
1874 "gps001", FatalException, << 1126 // Baryon Number looks to hold the no. of nucleons.
1875 "Histogram contains less tha << 1127 G4int Bary = particle_definition->GetBaryonNumber();
1876 return; << 1128 // G4cout << "Baryon No. " << Bary << G4endl;
1877 } << 1129 // Change values in histogram, Read it out, delete it, re-create it
1878 for (count = 0; count < maxcount; ++count << 1130 G4int count, maxcount;
1879 { << 1131 maxcount = G4int(EpnEnergyH.GetVectorLength());
1880 // Read out << 1132 // G4cout << maxcount << G4endl;
1881 ebins[count] = EpnEnergyH.GetLowEdgeEne << 1133 G4double ebins[1024],evals[1024];
1882 evals[count] = EpnEnergyH(count); << 1134 if(maxcount > 1024)
1883 } << 1135 {
1884 << 1136 G4cout << "Histogram contains more than 1024 bins!" << G4endl;
1885 // Multiply the channels by the nucleon n << 1137 G4cout << "Those above 1024 will be ignored" << G4endl;
1886 // << 1138 maxcount = 1024;
1887 for (count = 0; count < maxcount; ++count << 1139 }
1888 { << 1140 for(count=0;count<maxcount;count++)
1889 ebins[count] = ebins[count] * Bary; << 1141 {
1890 } << 1142 // Read out
1891 << 1143 ebins[count] = EpnEnergyH.GetLowEdgeEnergy(size_t(count));
1892 // Set Emin and Emax << 1144 evals[count] = EpnEnergyH(size_t(count));
1893 // << 1145 }
1894 params.Emin = ebins[0]; << 1146
1895 if (maxcount > 1) << 1147 // Multiply the channels by the nucleon number to give energies
1896 { << 1148 for(count=0;count<maxcount;count++)
1897 params.Emax = ebins[maxcount - 1]; << 1149 {
1898 } << 1150 ebins[count] = ebins[count] * Bary;
1899 else << 1151 }
1900 { << 1152
1901 params.Emax = ebins[0]; << 1153 // Set Emin and Emax
1902 } << 1154 Emin = ebins[0];
1903 << 1155 Emax = ebins[maxcount-1];
1904 // Put energy bins into new histogram - U << 1156
1905 // << 1157 // Put energy bins into new histogram - UDefEnergyH.
1906 for (count = 0; count < maxcount; ++count << 1158 for(count=0;count<maxcount;count++)
1907 { << 1159 {
1908 UDefEnergyH.InsertValues(ebins[count], << 1160 UDefEnergyH.InsertValues(ebins[count], evals[count]);
1909 } << 1161 }
1910 Epnflag = false; // so that you dont repe << 1162 Epnflag = false; //so that you dont repeat this method.
1911 } << 1163 }
1912 } 1164 }
1913 1165
1914 void G4SPSEneDistribution::ReSetHist(const G4 << 1166 //
>> 1167 void G4SPSEneDistribution::ReSetHist(G4String atype)
1915 { 1168 {
1916 G4AutoLock l(&mutex); << 1169 if (atype == "energy"){
1917 if (atype == "energy") << 1170 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ;
1918 { << 1171 IPDFEnergyExist = false ;
1919 UDefEnergyH = IPDFEnergyH = ZeroPhysVecto <<
1920 IPDFEnergyExist = false; <<
1921 Emin = 0.; 1172 Emin = 0.;
1922 Emax = 1e30; << 1173 Emax = 1e30;}
1923 } << 1174 else if ( atype == "arb"){
1924 else if (atype == "arb") << 1175 ArbEnergyH =IPDFArbEnergyH = ZeroPhysVector ;
1925 { << 1176 IPDFArbExist = false;}
1926 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVec << 1177 else if ( atype == "epn"){
1927 IPDFArbExist = false; << 1178 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ;
1928 } << 1179 IPDFEnergyExist = false ;
1929 else if (atype == "epn") << 1180 EpnEnergyH = ZeroPhysVector ;}
1930 { << 1181 else {
1931 UDefEnergyH = IPDFEnergyH = ZeroPhysVecto <<
1932 IPDFEnergyExist = false; <<
1933 EpnEnergyH = ZeroPhysVector; <<
1934 } <<
1935 else <<
1936 { <<
1937 G4cout << "Error, histtype not accepted " 1182 G4cout << "Error, histtype not accepted " << G4endl;
1938 } 1183 }
1939 } 1184 }
1940 1185
1941 G4double G4SPSEneDistribution::GenerateOne(G4 1186 G4double G4SPSEneDistribution::GenerateOne(G4ParticleDefinition* a)
1942 { 1187 {
1943 // Copy global shared status to thread-loca << 1188 particle_definition = a;
1944 // << 1189 particle_energy = -1.;
1945 threadLocal_t& params = threadLocalData.Get << 1190 while ( (EnergyDisType == "Arb")? (particle_energy < ArbEmin || particle_energy > ArbEmax)
1946 params.particle_definition=a; << 1191 : (particle_energy < Emin || particle_energy > Emax) ) {
1947 params.particle_energy=-1; << 1192 if(EnergyDisType == "Mono")
1948 if(applyEvergyWeight) << 1193 GenerateMonoEnergetic();
1949 { << 1194 else if(EnergyDisType == "Lin")
1950 params.Emax = ArbEmax; << 1195 GenerateLinearEnergies();
1951 params.Emin = ArbEmin; << 1196 else if(EnergyDisType == "Pow")
1952 } << 1197 GeneratePowEnergies();
1953 else << 1198 else if(EnergyDisType == "Exp")
1954 { << 1199 GenerateExpEnergies();
1955 params.Emax = Emax; << 1200 else if(EnergyDisType == "Gauss")
1956 params.Emin = Emin; << 1201 GenerateGaussEnergies();
1957 } << 1202 else if(EnergyDisType == "Brem")
1958 params.alpha = alpha; << 1203 GenerateBremEnergies();
1959 params.Ezero = Ezero; << 1204 else if(EnergyDisType == "Bbody")
1960 params.grad = grad; << 1205 GenerateBbodyEnergies();
1961 params.cept = cept; << 1206 else if(EnergyDisType == "Cdg")
1962 params.weight = weight; << 1207 GenerateCdgEnergies();
1963 // particle_energy = -1.; << 1208 else if(EnergyDisType == "User")
1964 << 1209 GenUserHistEnergies();
1965 if((EnergyDisType == "Mono") && ((MonoEnerg << 1210 else if(EnergyDisType == "Arb")
1966 { << 1211 GenArbPointEnergies();
1967 G4ExceptionDescription ed; << 1212 else if(EnergyDisType == "Epn")
1968 ed << "MonoEnergy " << G4BestUnit(MonoEne << 1213 GenEpnHistEnergies();
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 1214 else
1988 { << 1215 G4cout << "Error: EnergyDisType has unusual value" << G4endl;
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 } 1216 }
2043 return params.particle_energy; << 1217 return particle_energy;
2044 } 1218 }
2045 1219
2046 G4double G4SPSEneDistribution::GetProbability <<
2047 { <<
2048 G4double prob = 1.; <<
2049 1220
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 1221
2106 return prob; << 1222
2107 } << 1223
>> 1224
>> 1225
>> 1226
>> 1227
2108 1228