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