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