| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4SPSEneDistribution class implementation
27 //
28 // Author: Fan Lei, QinetiQ ltd - 05/02/2004
29 // Customer: ESA/ESTEC
30 // Revisions: Andrew Green, Andrea Dotti
31 // --------------------------------------------------------------------
32 #include "G4SPSEneDistribution.hh"
33
34 #include "G4Exp.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4UnitsTable.hh"
37 #include "Randomize.hh"
38 #include "G4AutoLock.hh"
39 #include "G4Threading.hh"
40
41 G4SPSEneDistribution::G4SPSEneDistribution()
42 {
43 G4MUTEXINIT(mutex);
44
45 // Initialise all variables
46
47 particle_energy = 1.0 * MeV;
48 EnergyDisType = "Mono";
49 weight=1.;
50 MonoEnergy = 1 * MeV;
51 Emin = 0.;
52 Emax = 1.e30;
53 alpha = 0.;
54 biasalpha = 0.;
55 prob_norm = 1.0;
56 Ezero = 0.;
57 SE = 0.;
58 Temp = 0.;
59 grad = 0.;
60 cept = 0.;
61 IntType = "NULL"; // Interpolation type
62
63 ArbEmin = 0.;
64 ArbEmax = 1.e30;
65
66 verbosityLevel = 0;
67
68 threadLocal_t& data = threadLocalData.Get();
69 data.Emax = Emax;
70 data.Emin = Emin;
71 data.alpha =alpha;
72 data.cept = cept;
73 data.Ezero = Ezero;
74 data.grad = grad;
75 data.particle_energy = 0.;
76 data.particle_definition = nullptr;
77 data.weight = weight;
78 }
79
80 G4SPSEneDistribution::~G4SPSEneDistribution()
81 {
82 G4MUTEXDESTROY(mutex);
83 if(Arb_grad_cept_flag)
84 {
85 delete [] Arb_grad;
86 delete [] Arb_cept;
87 }
88
89 if(Arb_alpha_Const_flag)
90 {
91 delete [] Arb_alpha;
92 delete [] Arb_Const;
93 }
94
95 if(Arb_ezero_flag)
96 {
97 delete [] Arb_ezero;
98 }
99 delete Bbody_x;
100 delete BBHist;
101 delete CP_x;
102 delete CPHist;
103 for (auto & it : SplineInt)
104 {
105 delete it;
106 it = nullptr;
107 }
108 SplineInt.clear();
109 }
110
111 void G4SPSEneDistribution::SetEnergyDisType(const G4String& DisType)
112 {
113 G4AutoLock l(&mutex);
114 EnergyDisType = DisType;
115 if (EnergyDisType == "User")
116 {
117 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
118 IPDFEnergyExist = false;
119 }
120 else if (EnergyDisType == "Arb")
121 {
122 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
123 IPDFArbExist = false;
124 }
125 else if (EnergyDisType == "Epn")
126 {
127 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
128 IPDFEnergyExist = false;
129 EpnEnergyH = ZeroPhysVector;
130 }
131 }
132
133 const G4String& G4SPSEneDistribution::GetEnergyDisType()
134 {
135 G4AutoLock l(&mutex);
136 return EnergyDisType;
137 }
138
139 void G4SPSEneDistribution::SetEmin(G4double emi)
140 {
141 G4AutoLock l(&mutex);
142 Emin = emi;
143 threadLocalData.Get().Emin = Emin;
144 }
145
146 G4double G4SPSEneDistribution::GetEmin() const
147 {
148 return threadLocalData.Get().Emin;
149 }
150
151 G4double G4SPSEneDistribution::GetArbEmin()
152 {
153 G4AutoLock l(&mutex);
154 return ArbEmin;
155 }
156
157 G4double G4SPSEneDistribution::GetArbEmax()
158 {
159 G4AutoLock l(&mutex);
160 return ArbEmax;
161 }
162
163 void G4SPSEneDistribution::SetEmax(G4double ema)
164 {
165 G4AutoLock l(&mutex);
166 Emax = ema;
167 threadLocalData.Get().Emax = Emax;
168 }
169
170 G4double G4SPSEneDistribution::GetEmax() const
171 {
172 return threadLocalData.Get().Emax;
173 }
174
175 void G4SPSEneDistribution::SetMonoEnergy(G4double menergy)
176 {
177 G4AutoLock l(&mutex);
178 MonoEnergy = menergy;
179 }
180
181 void G4SPSEneDistribution::SetBeamSigmaInE(G4double e)
182 {
183 G4AutoLock l(&mutex);
184 SE = e;
185 }
186 void G4SPSEneDistribution::SetAlpha(G4double alp)
187 {
188 G4AutoLock l(&mutex);
189 alpha = alp;
190 threadLocalData.Get().alpha = alpha;
191 }
192
193 void G4SPSEneDistribution::SetBiasAlpha(G4double alp)
194 {
195 G4AutoLock l(&mutex);
196 biasalpha = alp;
197 Biased = true;
198 }
199
200 void G4SPSEneDistribution::SetTemp(G4double tem)
201 {
202 G4AutoLock l(&mutex);
203 Temp = tem;
204 }
205
206 void G4SPSEneDistribution::SetEzero(G4double eze)
207 {
208 G4AutoLock l(&mutex);
209 Ezero = eze;
210 threadLocalData.Get().Ezero = Ezero;
211 }
212
213 void G4SPSEneDistribution::SetGradient(G4double gr)
214 {
215 G4AutoLock l(&mutex);
216 grad = gr;
217 threadLocalData.Get().grad = grad;
218 }
219
220 void G4SPSEneDistribution::SetInterCept(G4double c)
221 {
222 G4AutoLock l(&mutex);
223 cept = c;
224 threadLocalData.Get().cept = cept;
225 }
226
227 const G4String& G4SPSEneDistribution::GetIntType()
228 {
229 G4AutoLock l(&mutex);
230 return IntType;
231 }
232
233 void G4SPSEneDistribution::SetBiasRndm(G4SPSRandomGenerator* a)
234 {
235 G4AutoLock l(&mutex);
236 eneRndm = a;
237 }
238
239 void G4SPSEneDistribution::SetVerbosity(G4int a)
240 {
241 G4AutoLock l(&mutex);
242 verbosityLevel = a;
243 }
244
245 G4double G4SPSEneDistribution::GetWeight() const
246 {
247 return threadLocalData.Get().weight;
248 }
249
250 G4double G4SPSEneDistribution::GetMonoEnergy()
251 {
252 G4AutoLock l(&mutex);
253 return MonoEnergy;
254 }
255
256 G4double G4SPSEneDistribution::GetSE()
257 {
258 G4AutoLock l(&mutex);
259 return SE;
260 }
261
262 G4double G4SPSEneDistribution::Getalpha() const
263 {
264 return threadLocalData.Get().alpha;
265 }
266
267 G4double G4SPSEneDistribution::GetEzero() const
268 {
269 return threadLocalData.Get().Ezero;
270 }
271
272 G4double G4SPSEneDistribution::GetTemp()
273 {
274 G4AutoLock l(&mutex);
275 return Temp;
276 }
277
278 G4double G4SPSEneDistribution::Getgrad() const
279 {
280 return threadLocalData.Get().grad;
281 }
282
283 G4double G4SPSEneDistribution::Getcept() const
284 {
285 return threadLocalData.Get().cept;
286 }
287
288 G4PhysicsFreeVector G4SPSEneDistribution::GetUserDefinedEnergyHisto()
289 {
290 G4AutoLock l(&mutex);
291 return UDefEnergyH;
292 }
293
294 G4PhysicsFreeVector G4SPSEneDistribution::GetArbEnergyHisto()
295 {
296 G4AutoLock l(&mutex);
297 return ArbEnergyH;
298 }
299
300 void G4SPSEneDistribution::UserEnergyHisto(const G4ThreeVector& input)
301 {
302 G4AutoLock l(&mutex);
303 G4double ehi = input.x(),
304 val = input.y();
305 if (verbosityLevel > 1)
306 {
307 G4cout << "In UserEnergyHisto" << G4endl;
308 G4cout << " " << ehi << " " << val << G4endl;
309 }
310 UDefEnergyH.InsertValues(ehi, val);
311 Emax = ehi;
312 threadLocalData.Get().Emax = Emax;
313 }
314
315 void G4SPSEneDistribution::ArbEnergyHisto(const G4ThreeVector& input)
316 {
317 G4AutoLock l(&mutex);
318 G4double ehi = input.x(),
319 val = input.y();
320 if (verbosityLevel > 1)
321 {
322 G4cout << "In ArbEnergyHisto" << G4endl;
323 G4cout << " " << ehi << " " << val << G4endl;
324 }
325 ArbEnergyH.InsertValues(ehi, val);
326 }
327
328 void G4SPSEneDistribution::ArbEnergyHistoFile(const G4String& filename)
329 {
330 G4AutoLock l(&mutex);
331 std::ifstream infile(filename, std::ios::in);
332 if (!infile)
333 {
334 G4Exception("G4SPSEneDistribution::ArbEnergyHistoFile", "Event0301",
335 FatalException, "Unable to open the histo ASCII file");
336 }
337 G4double ehi, val;
338 while (infile >> ehi >> val)
339 {
340 ArbEnergyH.InsertValues(ehi, val);
341 }
342 }
343
344 void G4SPSEneDistribution::EpnEnergyHisto(const G4ThreeVector& input)
345 {
346 G4AutoLock l(&mutex);
347 G4double ehi = input.x(),
348 val = input.y();
349 if (verbosityLevel > 1)
350 {
351 G4cout << "In EpnEnergyHisto" << G4endl;
352 G4cout << " " << ehi << " " << val << G4endl;
353 }
354 EpnEnergyH.InsertValues(ehi, val);
355 Emax = ehi;
356 threadLocalData.Get().Emax = Emax;
357 Epnflag = true;
358 }
359
360 void G4SPSEneDistribution::Calculate()
361 {
362 G4AutoLock l(&mutex);
363 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() // MT: Lock in caller
386 {
387 BBHist = new std::vector<G4double>(10001, 0.0);
388 Bbody_x = new std::vector<G4double>(10001, 0.0);
389 BBhistInit = true;
390 }
391
392 void G4SPSEneDistribution::CPInitHists() // MT: Lock in caller
393 {
394 CPHist = new std::vector<G4double>(10001, 0.0);
395 CP_x = new std::vector<G4double>(10001, 0.0);
396 CPhistInit = true;
397 }
398
399 void G4SPSEneDistribution::CalculateCdgSpectrum() // MT: Lock in caller
400 {
401 // This uses the spectrum from the INTEGRAL Mass Model (TIMM)
402 // to generate a Cosmic Diffuse X/gamma ray spectrum.
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, 1E6 * 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] / omalpha)
437 * (std::pow(ene_line[i + 1] / keV, omalpha)
438 - std::pow(ene_line[i] / keV,omalpha));
439 ++i;
440 }
441
442 // Normalise histo and divide by 1000 to make MeV
443 //
444 i = 0;
445 while (i < n_par)
446 {
447 CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[n_par];
448 ++i;
449 }
450 }
451
452 void G4SPSEneDistribution::CalculateBbodySpectrum() // MT: Lock in caller
453 {
454 // Create bbody spectrum
455 // Proved very hard to integrate indefinitely, so different method.
456 // User inputs emin, emax and T. These are used to create a 10,000
457 // bin histogram.
458 // Use photon density spectrum = 2 nu**2/c**2 * (std::exp(h nu/kT)-1)
459 // = 2 E**2/h**2c**2 times the exponential
460
461 G4double erange = threadLocalData.Get().Emax - threadLocalData.Get().Emin;
462 G4double steps = erange / 10000.;
463
464 const G4double k = 8.6181e-11; //Boltzmann const in MeV/K
465 const G4double h = 4.1362e-21; // Plancks const in MeV s
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().Emin + G4double(count * steps);
476 G4double Bbody_y = (2. * std::pow(Bbody_x->at(count), 2.))
477 / (h2*c2*(std::exp(Bbody_x->at(count) / (k*Temp)) - 1.));
478 sum = sum + Bbody_y;
479 BBHist->at(count + 1) = BBHist->at(count) + Bbody_y;
480 ++count;
481 }
482
483 Bbody_x->at(10000) = threadLocalData.Get().Emax;
484
485 // Normalise cumulative histo
486 //
487 count = 0;
488 while (count < 10001)
489 {
490 BBHist->at(count) = BBHist->at(count) / sum;
491 ++count;
492 }
493 }
494
495 void G4SPSEneDistribution::CalculateCPowSpectrum() // MT: Lock in caller
496 {
497 // Create cutoff power-law spectrum, x^a exp(-x/b)
498 // The integral of this function is an incomplete gamma function, which
499 // is only available in the Boost library.
500 //
501 // User inputs are emin, emax and alpha and Ezero. These are used to
502 // create a 10,000 bin histogram.
503
504 G4double erange = threadLocalData.Get().Emax - threadLocalData.Get().Emin;
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().Emin + G4double(count * steps);
516 G4double CP_y = std::pow(CP_x->at(count), alpha)
517 * std::exp(-CP_x->at(count) / Ezero);
518 sum = sum + CP_y;
519 CPHist->at(count + 1) = CPHist->at(count) + CP_y;
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) / sum;
531 ++count;
532 }
533 }
534
535 void G4SPSEneDistribution::InputEnergySpectra(G4bool value)
536 {
537 G4AutoLock l(&mutex);
538
539 // Allows user to specify spectrum is momentum
540 //
541 EnergySpec = value; // false if momentum
542 if (verbosityLevel > 1)
543 {
544 G4cout << "EnergySpec has value " << EnergySpec << G4endl;
545 }
546 }
547
548 void G4SPSEneDistribution::InputDifferentialSpectra(G4bool value)
549 {
550 G4AutoLock l(&mutex);
551
552 // Allows user to specify integral or differential spectra
553 //
554 DiffSpec = value; // true = differential, false = integral
555 if (verbosityLevel > 1)
556 {
557 G4cout << "Diffspec has value " << DiffSpec << G4endl;
558 }
559 }
560
561 void G4SPSEneDistribution::ArbInterpolate(const G4String& IType)
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() // MT: Lock in caller
578 {
579 // Method to do linear interpolation on the Arb points
580 // Calculate equation of each line segment, max 1024.
581 // Calculate Area under each segment
582 // Create a cumulative array which is then normalised Arb_Cum_Area
583
584 G4double Area_seg[1024]; // Stores area under each segment
585 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1024]={0.}, Arb_Cum_Area[1024]={0.};
586 std::size_t i, count;
587 std::size_t maxi = ArbEnergyH.GetVectorLength();
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 spectrum is integral it has to be
595 // made differential and if momentum it has to be made energy
596
597 if (!DiffSpec)
598 {
599 // Converts integral point-wise spectra to Differential
600 //
601 for (count = 0; count < maxi-1; ++count)
602 {
603 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
604 / (Arb_x[count + 1] - Arb_x[count]);
605 }
606 --maxi;
607 }
608
609 if (!EnergySpec)
610 {
611 // change currently stored values (emin etc) which are actually momenta
612 // to energies
613 //
614 G4ParticleDefinition* pdef = threadLocalData.Get().particle_definition;
615 if (pdef == nullptr)
616 {
617 G4Exception("G4SPSEneDistribution::LinearInterpolation",
618 "Event0302", FatalException,
619 "Error: particle not defined");
620 }
621 else
622 {
623 // Apply Energy**2 = p**2c**2 + m0**2c**4
624 // p should be entered as E/c i.e. without the division by c
625 // being done - energy equivalent
626
627 G4double mass = pdef->GetPDGMass();
628
629 // Convert point to energy unit and its value to per energy unit
630 //
631 G4double total_energy;
632 for (count = 0; count < maxi; ++count)
633 {
634 total_energy = std::sqrt((Arb_x[count] * Arb_x[count])
635 + (mass * mass)); // total energy
636 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
637 Arb_x[count] = total_energy - mass; // kinetic energy
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 each segment
655 //
656 Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[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" << G4endl;
666 }
667 Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]);
668 }
669 else if (Arb_grad[i] < 0.)
670 {
671 if (verbosityLevel == 2)
672 {
673 G4cout << "Arb_grad is negative" << G4endl;
674 }
675 Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]);
676 }
677 else
678 {
679 if (verbosityLevel == 2)
680 {
681 G4cout << "Arb_grad is 0." << G4endl;
682 }
683 Arb_cept[i] = Arb_y[i];
684 }
685
686 Area_seg[i] = ((Arb_grad[i] / 2)
687 * (Arb_x[i] * Arb_x[i] - Arb_x[i - 1] * Arb_x[i - 1])
688 + Arb_cept[i] * (Arb_x[i] - Arb_x[i - 1]));
689 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
690 sum = sum + Area_seg[i];
691 if (verbosityLevel == 2)
692 {
693 G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum
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; // normalisation
703 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
704 ++i;
705 }
706
707 // now scale the ArbEnergyH, needed by Probability()
708 //
709 ArbEnergyH.ScaleVector(1., 1./sum);
710
711 if (verbosityLevel >= 1)
712 {
713 G4cout << "Leaving LinearInterpolation" << G4endl;
714 ArbEnergyH.DumpValues();
715 IPDFArbEnergyH.DumpValues();
716 }
717 }
718
719 void G4SPSEneDistribution::LogInterpolation() // MT: Lock in caller
720 {
721 // Interpolation based on Logarithmic equations
722 // Generate equations of line segments
723 // y = Ax**alpha => log y = alpha*logx + logA
724 // Find area under line segments
725 // Create normalised, cumulative array Arb_Cum_Area
726
727 G4double Area_seg[1024]; // Stores area under each segment
728 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1024]={0.}, Arb_Cum_Area[1024]={0.};
729 std::size_t i, count;
730 std::size_t maxi = ArbEnergyH.GetVectorLength();
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 spectrum is integral it has to be
738 // made differential and if momentum it has to be made energy
739
740 if (!DiffSpec)
741 {
742 // Converts integral point-wise spectra to Differential
743 //
744 for (count = 0; count < maxi-1; ++count)
745 {
746 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
747 / (Arb_x[count + 1] - Arb_x[count]);
748 }
749 --maxi;
750 }
751
752 if (!EnergySpec)
753 {
754 // Change currently stored values (emin etc) which are actually momenta
755 // to energies
756
757 G4ParticleDefinition* pdef = threadLocalData.Get().particle_definition;
758 if (pdef == nullptr)
759 {
760 G4Exception("G4SPSEneDistribution::LogInterpolation",
761 "Event0302", FatalException,
762 "Error: particle not defined");
763 }
764 else
765 {
766 // Apply Energy**2 = p**2c**2 + m0**2c**4
767 // p should be entered as E/c i.e. without the division by c
768 // being done - energy equivalent
769
770 G4double mass = pdef->GetPDGMass();
771
772 // Convert point to energy unit and its value to per energy unit
773 //
774 G4double total_energy;
775 for (count = 0; count < maxi; ++count)
776 {
777 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass * mass));
778 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
779 Arb_x[count] = total_energy - mass; // kinetic energy
780 }
781 }
782 }
783
784 i = 1;
785
786 if ( Arb_ezero != nullptr ) { delete [] Arb_ezero; Arb_ezero = nullptr; }
787 if ( Arb_Const != nullptr ) { delete [] Arb_Const; Arb_Const = nullptr; }
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 interpolation with points <= 0."
799 << G4endl;
800 G4cout << "These will be changed to 1e-20, which may cause problems"
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 interpolation with points <= 0."
820 << G4endl;
821 G4cout << "These will be changed to 1e-20, which may cause problems"
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::log10(Arb_y[i - 1]))
834 / (std::log10(Arb_x[i]) - std::log10(Arb_x[i - 1]));
835 Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[i], Arb_alpha[i]));
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::log(Arb_x[i - 1]));
841 }
842 else
843 {
844 Area_seg[i] = (Arb_Const[i] / alp)
845 * (std::pow(Arb_x[i], alp) - std::pow(Arb_x[i - 1], alp));
846 }
847 sum = sum + Area_seg[i];
848 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
849 if (verbosityLevel == 2)
850 {
851 G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl;
852 }
853 ++i;
854 }
855
856 i = 0;
857 while (i < maxi)
858 {
859 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
860 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
861 ++i;
862 }
863
864 // Now scale the ArbEnergyH, needed by Probability()
865 //
866 ArbEnergyH.ScaleVector(1., 1./sum);
867
868 if (verbosityLevel >= 1)
869 {
870 G4cout << "Leaving LogInterpolation " << G4endl;
871 }
872 }
873
874 void G4SPSEneDistribution::ExpInterpolation() // MT: Lock in caller
875 {
876 // Interpolation based on Exponential equations
877 // Generate equations of line segments
878 // y = Ae**-(x/e0) => ln y = -x/e0 + lnA
879 // Find area under line segments
880 // Create normalised, cumulative array Arb_Cum_Area
881
882 G4double Area_seg[1024]; // Stores area under each segment
883 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1024]={0.}, Arb_Cum_Area[1024]={0.};
884 std::size_t i, count;
885 std::size_t maxi = ArbEnergyH.GetVectorLength();
886 for (i = 0; i < maxi; ++i)
887 {
888 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i);
889 Arb_y[i] = ArbEnergyH(i);
890 }
891
892 // Points are now in x,y arrays. If the spectrum is integral it has to be
893 // made differential and if momentum it has to be made energy
894
895 if (!DiffSpec)
896 {
897 // Converts integral point-wise spectra to Differential
898 //
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)
908 {
909 // Change currently stored values (emin etc) which are actually momenta
910 // to energies
911 //
912 G4ParticleDefinition* pdef = threadLocalData.Get().particle_definition;
913 if (pdef == nullptr)
914 {
915 G4Exception("G4SPSEneDistribution::ExpInterpolation",
916 "Event0302", FatalException,
917 "Error: particle not defined");
918 }
919 else
920 {
921 // Apply Energy**2 = p**2c**2 + m0**2c**4
922 // p should be entered as E/c i.e. without the division by c
923 // being done - energy equivalent
924
925 G4double mass = pdef->GetPDGMass();
926
927 // Convert point to energy unit and its value to per energy unit
928 //
929 G4double total_energy;
930 for (count = 0; count < maxi; ++count)
931 {
932 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass * mass));
933 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
934 Arb_x[count] = total_energy - mass; // kinetic energy
935 }
936 }
937 }
938
939 i = 1;
940
941 if ( Arb_ezero != nullptr ) { delete[] Arb_ezero; Arb_ezero = nullptr; }
942 if ( Arb_Const != nullptr ) { delete[] Arb_Const; Arb_Const = nullptr; }
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::log(Arb_y[i - 1]);
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::log(Arb_y[i - 1]));
958 Arb_Const[i] = Arb_y[i] / (std::exp(-Arb_x[i] / Arb_ezero[i]));
959 Area_seg[i] = -(Arb_Const[i] * Arb_ezero[i])
960 * (std::exp(-Arb_x[i] / Arb_ezero[i])
961 - std::exp(-Arb_x[i - 1] / Arb_ezero[i]));
962 }
963 else
964 {
965 G4Exception("G4SPSEneDistribution::ExpInterpolation",
966 "Event0302", JustWarning,
967 "Flat line segment: problem, setting to zero parameters.");
968 G4cout << "Flat line segment: problem" << G4endl;
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] + Area_seg[i];
975 if (verbosityLevel == 2)
976 {
977 G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl;
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_Cum_Area[i]);
987 ++i;
988 }
989
990 // Now scale the ArbEnergyH, needed by Probability()
991 //
992 ArbEnergyH.ScaleVector(1., 1./sum);
993
994 if (verbosityLevel >= 1)
995 {
996 G4cout << "Leaving ExpInterpolation " << G4endl;
997 }
998 }
999
1000 void G4SPSEneDistribution::SplineInterpolation() // MT: Lock in caller
1001 {
1002 // Interpolation using Splines.
1003 // Create Normalised arrays, make x 0->1 and y hold the function (Energy)
1004 //
1005 // Current method based on the above will not work in all cases.
1006 // New method is implemented below.
1007
1008 G4double sum, Arb_x[1024]={0.}, Arb_y[1024]={0.}, Arb_Cum_Area[1024]={0.};
1009 std::size_t i, count;
1010 std::size_t maxi = ArbEnergyH.GetVectorLength();
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 spectrum is integral it has to be
1019 // made differential and if momentum it has to be made energy
1020
1021 if (!DiffSpec)
1022 {
1023 // Converts integral point-wise spectra to Differential
1024 //
1025 for (count = 0; count < maxi - 1; ++count)
1026 {
1027 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
1028 / (Arb_x[count + 1] - Arb_x[count]);
1029 }
1030 --maxi;
1031 }
1032
1033 if (!EnergySpec)
1034 {
1035 // Change currently stored values (emin etc) which are actually momenta
1036 // to energies
1037 //
1038 G4ParticleDefinition* pdef = threadLocalData.Get().particle_definition;
1039 if (pdef == nullptr)
1040 {
1041 G4Exception("G4SPSEneDistribution::SplineInterpolation",
1042 "Event0302", FatalException,
1043 "Error: particle not defined");
1044 }
1045 else
1046 {
1047 // Apply Energy**2 = p**2c**2 + m0**2c**4
1048 // p should be entered as E/c i.e. without the division by c
1049 // being done - energy equivalent
1050
1051 G4double mass = pdef->GetPDGMass();
1052
1053 // Convert point to energy unit and its value to per energy unit
1054 //
1055 G4double total_energy;
1056 for (count = 0; count < maxi; ++count)
1057 {
1058 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass * mass));
1059 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
1060 Arb_x[count] = total_energy - mass; // kinetic energy
1061 }
1062 }
1063 }
1064
1065 i = 1;
1066 Arb_Cum_Area[0] = 0.;
1067 sum = 0.;
1068 Splinetemp = new G4DataInterpolation(Arb_x, Arb_y, (G4int)maxi, 0., 0.);
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 integration of area
1080
1081 G4double de = (Arb_x[i] - Arb_x[i - 1])/100.;
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->CubicSplineInterpolation(ei[count]);
1088 if (prob[count] < 0.)
1089 {
1090 G4ExceptionDescription ED;
1091 ED << "Warning: G4DataInterpolation returns value < 0 " << prob[count]
1092 << " " << ei[count] << G4endl;
1093 G4Exception("G4SPSEneDistribution::SplineInterpolation", "Event0303",
1094 FatalException, ED);
1095 }
1096 area += prob[count]*de;
1097 }
1098 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + area;
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[count]/(area/de);
1105 }
1106
1107 SplineInt[i] = new G4DataInterpolation(prob, ei, 101, 0., 0.);
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; // normalisation
1118 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
1119 ++i;
1120 }
1121
1122 // Now scale the ArbEnergyH, needed by Probability()
1123 //
1124 ArbEnergyH.ScaleVector(1., 1./sum);
1125
1126 if (verbosityLevel > 0)
1127 {
1128 G4cout << "Leaving SplineInterpolation " << G4endl;
1129 }
1130 }
1131
1132 void G4SPSEneDistribution::GenerateMonoEnergetic()
1133 {
1134 // Method to generate MonoEnergetic particles
1135
1136 threadLocalData.Get().particle_energy = MonoEnergy;
1137 }
1138
1139 void G4SPSEneDistribution::GenerateGaussEnergies()
1140 {
1141 // Method to generate Gaussian particles
1142
1143 G4double ene = G4RandGauss::shoot(MonoEnergy,SE);
1144 if (ene < 0) ene = 0.;
1145 threadLocalData.Get().particle_energy = ene;
1146 }
1147
1148 void G4SPSEneDistribution::GenerateLinearEnergies(G4bool bArb = false)
1149 {
1150 G4double rndm;
1151 threadLocal_t& params = threadLocalData.Get();
1152 G4double emaxsq = std::pow(params.Emax, 2.); // Emax squared
1153 G4double eminsq = std::pow(params.Emin, 2.); // Emin squared
1154 G4double intersq = std::pow(params.cept, 2.); // cept squared
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.Emax - params.Emin));
1162 bracket = bracket * rndm;
1163 bracket = bracket + (params.grad / 2.) * eminsq + params.cept * params.Emin;
1164
1165 // Now have a quad of form m/2 E**2 + cE - bracket = 0
1166 //
1167 bracket = -bracket;
1168
1169 if (params.grad != 0.)
1170 {
1171 G4double sqbrack = (intersq - 4 * (params.grad / 2.) * (bracket));
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.Emax)
1180 {
1181 params.particle_energy = root1;
1182 }
1183 if (root2 > params.Emin && root2 < params.Emax)
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.cept;
1193 }
1194
1195 if (params.particle_energy < 0.)
1196 {
1197 params.particle_energy = -params.particle_energy;
1198 }
1199
1200 if (verbosityLevel >= 1)
1201 {
1202 G4cout << "Energy is " << params.particle_energy << G4endl;
1203 }
1204 }
1205
1206 void G4SPSEneDistribution::GeneratePowEnergies(G4bool bArb = false)
1207 {
1208 // Method to generate particle energies distributed as a power-law
1209
1210 G4double rndm;
1211 G4double emina, emaxa;
1212
1213 threadLocal_t& params = threadLocalData.Get();
1214
1215 emina = std::pow(params.Emin, params.alpha + 1);
1216 emaxa = std::pow(params.Emax, params.alpha + 1);
1217
1218 if (bArb) rndm = G4UniformRand();
1219 else rndm = eneRndm->GenRandEnergy();
1220
1221 if (params.alpha != -1.)
1222 {
1223 G4double ene = ((rndm * (emaxa - emina)) + emina);
1224 ene = std::pow(ene, (1. / (params.alpha + 1.)));
1225 params.particle_energy = ene;
1226 }
1227 else
1228 {
1229 G4double ene = (std::log(params.Emin)
1230 + rndm * (std::log(params.Emax) - std::log(params.Emin)));
1231 params.particle_energy = std::exp(ene);
1232 }
1233 if (verbosityLevel >= 1)
1234 {
1235 G4cout << "Energy is " << params.particle_energy << G4endl;
1236 }
1237 }
1238
1239 void G4SPSEneDistribution::GenerateCPowEnergies()
1240 {
1241 // Method to generate particle energies distributed in
1242 // cutoff power-law distribution
1243 //
1244 // CP_x holds Energies, and CPHist holds the cumulative histo.
1245 // binary search to find correct bin then lin interpolation.
1246 // Use the earlier defined histogram + RandGeneral method to generate
1247 // random numbers following the histos distribution
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, risk is to do it twice
1259 l.lock();
1260 CPhistCalcd = true;
1261 l.unlock();
1262 }
1263
1264 // Binary search to find bin that rndm is in
1265 //
1266 while (nabove - nbelow > 1)
1267 {
1268 middle = (nabove + nbelow) / 2;
1269 if (rndm == CPHist->at(middle))
1270 {
1271 break;
1272 }
1273 if (rndm < CPHist->at(middle))
1274 {
1275 nabove = middle;
1276 }
1277 else
1278 {
1279 nbelow = middle;
1280 }
1281 }
1282
1283 // Now interpolate in that bin to find the correct output value
1284 //
1285 G4double x1, x2, y1, y2, t, q;
1286 x1 = CP_x->at(nbelow);
1287 if(nbelow+1 == static_cast<G4int>(CP_x->size()))
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->size()))
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 = (rndm - q) / t;
1309
1310 if (verbosityLevel >= 1)
1311 {
1312 G4cout << "Energy is " << threadLocalData.Get().particle_energy << G4endl;
1313 }
1314 }
1315
1316 void G4SPSEneDistribution::GenerateBiasPowEnergies()
1317 {
1318 // Method to generate particle energies distributed as
1319 // in biased power-law and calculate its weight
1320
1321 threadLocal_t& params = threadLocalData.Get();
1322
1323 G4double rndm;
1324 G4double emina, emaxa, emin, emax;
1325
1326 G4double normal = 1.;
1327
1328 emin = params.Emin;
1329 emax = params.Emax;
1330
1331 rndm = eneRndm->GenRandEnergy();
1332
1333 if (biasalpha != -1.)
1334 {
1335 emina = std::pow(emin, biasalpha + 1);
1336 emaxa = std::pow(emax, biasalpha + 1);
1337 G4double ee = ((rndm * (emaxa - emina)) + emina);
1338 params.particle_energy = std::pow(ee, (1. / (biasalpha + 1.)));
1339 normal = 1./(1+biasalpha) * (emaxa - emina);
1340 }
1341 else
1342 {
1343 G4double ee = (std::log(emin) + rndm * (std::log(emax) - std::log(emin)));
1344 params.particle_energy = std::exp(ee);
1345 normal = std::log(emax) - std::log(emin);
1346 }
1347 params.weight = GetProbability(params.particle_energy)
1348 / (std::pow(params.particle_energy,biasalpha)/normal);
1349
1350 if (verbosityLevel >= 1)
1351 {
1352 G4cout << "Energy is " << params.particle_energy << G4endl;
1353 }
1354 }
1355
1356 void G4SPSEneDistribution::GenerateExpEnergies(G4bool bArb = false)
1357 {
1358 // Method to generate particle energies distributed according
1359 // to an exponential curve
1360
1361 G4double rndm;
1362
1363 if (bArb) rndm = G4UniformRand();
1364 else rndm = eneRndm->GenRandEnergy();
1365
1366 threadLocal_t& params = threadLocalData.Get();
1367 params.particle_energy = -params.Ezero
1368 * (std::log(rndm * (std::exp(-params.Emax
1369 / params.Ezero)
1370 - std::exp(-params.Emin
1371 / params.Ezero))
1372 + std::exp(-params.Emin / params.Ezero)));
1373 if (verbosityLevel >= 1)
1374 {
1375 G4cout << "Energy is " << params.particle_energy << G4endl;
1376 }
1377 }
1378
1379 void G4SPSEneDistribution::GenerateBremEnergies()
1380 {
1381 // Method to generate particle energies distributed according
1382 // to a Bremstrahlung equation of the form
1383 // I = const*((kT)**1/2)*E*(e**(-E/kT))
1384
1385 G4double rndm = eneRndm->GenRandEnergy();
1386 G4double expmax, expmin, k;
1387
1388 k = 8.6181e-11; // Boltzmann's const in MeV/K
1389 G4double ksq = std::pow(k, 2.); // k squared
1390 G4double Tsq = std::pow(Temp, 2.); // Temp squared
1391
1392 threadLocal_t& params = threadLocalData.Get();
1393
1394 expmax = std::exp(-params.Emax / (k * Temp));
1395 expmin = std::exp(-params.Emin / (k * Temp));
1396
1397 // If either expmax or expmin are zero then this will cause problems
1398 // Most probably this will be because T is too low or E is too high
1399
1400 if (expmax == 0.)
1401 {
1402 G4Exception("G4SPSEneDistribution::GenerateBremEnergies",
1403 "Event0302", FatalException,
1404 "*****EXPMAX=0. Choose different E's or Temp");
1405 }
1406 if (expmin == 0.)
1407 {
1408 G4Exception("G4SPSEneDistribution::GenerateBremEnergies",
1409 "Event0302", FatalException,
1410 "*****EXPMIN=0. Choose different E's or Temp");
1411 }
1412
1413 G4double tempvar = rndm * ((-k) * Temp * (params.Emax * expmax
1414 - params.Emin * expmin)
1415 - (ksq * Tsq * (expmax - expmin)));
1416
1417 G4double bigc = (tempvar - k * Temp * params.Emin * expmin
1418 - ksq * Tsq * expmin) / (-k * Temp);
1419
1420 // This gives an equation of form: Ee(-E/kT) + kTe(-E/kT) - C =0
1421 // Solve this iteratively, step from Emin to Emax in 1000 steps
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 * Temp)))
1433 + k * Temp * (std::exp(-etest / (k * Temp))) - bigc;
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_energy << G4endl;
1449 }
1450 }
1451
1452 void G4SPSEneDistribution::GenerateBbodyEnergies()
1453 {
1454 // BBody_x holds Energies, and BBHist holds the cumulative histo.
1455 // Binary search to find correct bin then lin interpolation.
1456 // Use the earlier defined histogram + RandGeneral method to generate
1457 // random numbers following the histos distribution
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, risk is to do it twice
1469 l.lock();
1470 BBhistCalcd = true;
1471 l.unlock();
1472 }
1473
1474 // Binary search to find bin that rndm is in
1475 //
1476 while (nabove - nbelow > 1)
1477 {
1478 middle = (nabove + nbelow) / 2;
1479 if (rndm == BBHist->at(middle))
1480 {
1481 break;
1482 }
1483 if (rndm < BBHist->at(middle))
1484 {
1485 nabove = middle;
1486 }
1487 else
1488 {
1489 nbelow = middle;
1490 }
1491 }
1492
1493 // Now interpolate in that bin to find the correct output value
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->size()))
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->size()))
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 = (rndm - q) / t;
1520
1521 if (verbosityLevel >= 1)
1522 {
1523 G4cout << "Energy is " << threadLocalData.Get().particle_energy << G4endl;
1524 }
1525 }
1526
1527 void G4SPSEneDistribution::GenerateCdgEnergies()
1528 {
1529 // Generate random numbers, compare with values in cumhist
1530 // to find appropriate part of spectrum and then
1531 // generate energy in the usual inversion way
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 < 18 * keV)
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 > 18 * keV)
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], omalpha[i - 1])
1569 + (std::pow(ene_line[i], omalpha[i - 1])
1570 - std::pow(ene_line[i - 1], omalpha[i- 1])) * rndm2);
1571 params.particle_energy = std::pow(ene, (1. / omalpha[i - 1]));
1572
1573 if (verbosityLevel >= 1)
1574 {
1575 G4cout << "Energy is " << params.particle_energy << G4endl;
1576 }
1577 }
1578
1579 void G4SPSEneDistribution::GenUserHistEnergies()
1580 {
1581 // Histograms are DIFFERENTIAL
1582
1583 G4AutoLock l(&mutex);
1584
1585 if (!IPDFEnergyExist)
1586 {
1587 std::size_t ii;
1588 std::size_t maxbin = UDefEnergyH.GetVectorLength();
1589 G4double bins[1024], vals[1024], sum;
1590 for ( ii = 0 ; ii<1024 ; ++ii ) { bins[ii]=0; vals[ii]=0; }
1591 sum = 0.;
1592
1593 if ( (!EnergySpec)
1594 && (threadLocalData.Get().particle_definition == nullptr))
1595 {
1596 G4Exception("G4SPSEneDistribution::GenUserHistEnergies",
1597 "Event0302", FatalException,
1598 "Error: particle definition is NULL");
1599 }
1600
1601 if (maxbin > 1024)
1602 {
1603 G4Exception("G4SPSEneDistribution::GenUserHistEnergies",
1604 "Event0302", JustWarning,
1605 "Maxbin>1024\n Setting maxbin to 1024, other bins are lost");
1606 maxbin = 1024;
1607 }
1608
1609 if (!DiffSpec)
1610 {
1611 G4cout << "Histograms are Differential!!! " << G4endl;
1612 }
1613 else
1614 {
1615 bins[0] = UDefEnergyH.GetLowEdgeEnergy(0);
1616 vals[0] = UDefEnergyH(0);
1617 sum = vals[0];
1618 for (ii = 1; ii < maxbin; ++ii)
1619 {
1620 bins[ii] = UDefEnergyH.GetLowEdgeEnergy(ii);
1621 vals[ii] = UDefEnergyH(ii) + vals[ii - 1];
1622 sum = sum + UDefEnergyH(ii);
1623 }
1624 }
1625
1626 if (!EnergySpec)
1627 {
1628 G4double mass = threadLocalData.Get().particle_definition->GetPDGMass();
1629
1630 // Multiply the function (vals) up by the bin width
1631 // to make the function counts/s (i.e. get rid of momentum dependence)
1632
1633 for (ii = 1; ii < maxbin; ++ii)
1634 {
1635 vals[ii] = vals[ii] * (bins[ii] - bins[ii - 1]);
1636 }
1637
1638 // Put energy bins into new histo, plus divide by energy bin width
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[ii])+(mass*mass))-mass;
1646 }
1647 for (ii = 1; ii < maxbin; ++ii)
1648 {
1649 vals[ii] = vals[ii] / (bins[ii] - bins[ii - 1]);
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[ii]);
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= IPDFEnergyH.GetEnergy(rndm);
1672
1673 if (verbosityLevel >= 1)
1674 {
1675 G4cout << "Energy is " << particle_energy << G4endl;
1676 }
1677 }
1678
1679 G4double G4SPSEneDistribution::GetArbEneWeight(G4double ene)
1680 {
1681 auto nbelow = IPDFArbEnergyH.FindBin(ene,(IPDFArbEnergyH.GetVectorLength())/2);
1682 G4double wei = 0.;
1683 if(IntType=="Lin")
1684 {
1685 // note: grad[i] and cept[i] are calculated with x[i-1] and x[i]
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]->CubicSplineInterpolation(ene);
1705 }
1706 return wei;
1707 }
1708
1709 void G4SPSEneDistribution::GenArbPointEnergies()
1710 {
1711 if (verbosityLevel > 0)
1712 {
1713 G4cout << "In GenArbPointEnergies" << G4endl;
1714 }
1715
1716 G4double rndm = eneRndm->GenRandEnergy();
1717
1718 // Find the Bin, have x, y, no of points, and cumulative area distribution
1719 //
1720 std::size_t nabove = IPDFArbEnergyH.GetVectorLength(), nbelow = 0, middle;
1721
1722 // Binary search to find bin that rndm is in
1723 //
1724 while (nabove - nbelow > 1)
1725 {
1726 middle = (nabove + nbelow) / 2;
1727 if (rndm == IPDFArbEnergyH(middle))
1728 {
1729 break;
1730 }
1731 if (rndm < IPDFArbEnergyH(middle))
1732 {
1733 nabove = middle;
1734 }
1735 else
1736 {
1737 nbelow = middle;
1738 }
1739 }
1740 threadLocal_t& params = threadLocalData.Get();
1741 if (IntType == "Lin")
1742 {
1743 // Update thread-local copy of parameters
1744 //
1745 params.Emax = IPDFArbEnergyH.GetLowEdgeEnergy(nbelow + 1);
1746 params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(nbelow);
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.GetLowEdgeEnergy(nbelow + 1);
1754 params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(nbelow);
1755 params.alpha = Arb_alpha[nbelow + 1];
1756 GeneratePowEnergies(true);
1757 }
1758 else if (IntType == "Exp")
1759 {
1760 params.Emax = IPDFArbEnergyH.GetLowEdgeEnergy(nbelow + 1);
1761 params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(nbelow);
1762 params.Ezero = Arb_ezero[nbelow + 1];
1763 GenerateExpEnergies(true);
1764 }
1765 else if (IntType == "Spline")
1766 {
1767 params.Emax = IPDFArbEnergyH.GetLowEdgeEnergy(nbelow + 1);
1768 params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(nbelow);
1769 params.particle_energy = -1e100;
1770 rndm = eneRndm->GenRandEnergy();
1771 while (params.particle_energy < params.Emin
1772 || params.particle_energy > params.Emax)
1773 {
1774 params.particle_energy =
1775 SplineInt[nbelow+1]->CubicSplineInterpolation(rndm);
1776 rndm = eneRndm->GenRandEnergy();
1777 }
1778 if (verbosityLevel >= 1)
1779 {
1780 G4cout << "Energy is " << params.particle_energy << G4endl;
1781 }
1782 }
1783 else
1784 {
1785 G4Exception("G4SPSEneDistribution::GenArbPointEnergies", "Event0302",
1786 FatalException, "Error: IntType unknown type");
1787 }
1788 }
1789
1790 void G4SPSEneDistribution::GenEpnHistEnergies()
1791 {
1792 // Firstly convert to energy if not already done
1793
1794 G4AutoLock l(&mutex);
1795
1796 if (Epnflag) // true means spectrum is epn, false means e
1797 {
1798 // Convert to energy by multiplying by A number
1799 //
1800 ConvertEPNToEnergy();
1801 }
1802 if (!IPDFEnergyExist)
1803 {
1804 // IPDF has not been created, so create it
1805 //
1806 G4double bins[1024], vals[1024], sum;
1807 std::size_t ii;
1808 std::size_t maxbin = UDefEnergyH.GetVectorLength();
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(ii);
1815 vals[ii] = UDefEnergyH(ii) + vals[ii - 1];
1816 sum = sum + UDefEnergyH(ii);
1817 }
1818
1819 l.lock();
1820 for (ii = 0; ii < maxbin; ++ii)
1821 {
1822 vals[ii] = vals[ii] / sum;
1823 IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
1824 }
1825 IPDFEnergyExist = true;
1826
1827 }
1828 l.unlock();
1829
1830 // IPDF has been create so carry on
1831 //
1832 G4double rndm = eneRndm->GenRandEnergy();
1833 threadLocalData.Get().particle_energy = IPDFEnergyH.GetEnergy(rndm);
1834
1835 if (verbosityLevel >= 1)
1836 {
1837 G4cout << "Energy is " << threadLocalData.Get().particle_energy << G4endl;
1838 }
1839 }
1840
1841 void G4SPSEneDistribution::ConvertEPNToEnergy() // MT: lock in caller
1842 {
1843 // Use this before particle generation to convert the
1844 // currently stored histogram from energy/nucleon to energy.
1845
1846 threadLocal_t& params = threadLocalData.Get();
1847 if (params.particle_definition == nullptr)
1848 {
1849 G4cout << "Error: particle not defined" << G4endl;
1850 }
1851 else
1852 {
1853 // Need to multiply histogram by the number of nucleons.
1854 // Baryon Number looks to hold the no. of nucleons
1855 //
1856 G4int Bary = params.particle_definition->GetBaryonNumber();
1857
1858 // Change values in histogram, Read it out, delete it, re-create it
1859 //
1860 std::size_t count, maxcount;
1861 maxcount = EpnEnergyH.GetVectorLength();
1862 G4double ebins[1024], evals[1024];
1863 if (maxcount > 1024)
1864 {
1865 G4Exception("G4SPSEneDistribution::ConvertEPNToEnergy()",
1866 "gps001", JustWarning,
1867 "Histogram contains more than 1024 bins!\n\
1868 Those above 1024 will be ignored");
1869 maxcount = 1024;
1870 }
1871 if (maxcount < 1)
1872 {
1873 G4Exception("G4SPSEneDistribution::ConvertEPNToEnergy()",
1874 "gps001", FatalException,
1875 "Histogram contains less than 1 bin!\nRedefine the histogram");
1876 return;
1877 }
1878 for (count = 0; count < maxcount; ++count)
1879 {
1880 // Read out
1881 ebins[count] = EpnEnergyH.GetLowEdgeEnergy(count);
1882 evals[count] = EpnEnergyH(count);
1883 }
1884
1885 // Multiply the channels by the nucleon number to give energies
1886 //
1887 for (count = 0; count < maxcount; ++count)
1888 {
1889 ebins[count] = ebins[count] * Bary;
1890 }
1891
1892 // Set Emin and Emax
1893 //
1894 params.Emin = ebins[0];
1895 if (maxcount > 1)
1896 {
1897 params.Emax = ebins[maxcount - 1];
1898 }
1899 else
1900 {
1901 params.Emax = ebins[0];
1902 }
1903
1904 // Put energy bins into new histogram - UDefEnergyH
1905 //
1906 for (count = 0; count < maxcount; ++count)
1907 {
1908 UDefEnergyH.InsertValues(ebins[count], evals[count]);
1909 }
1910 Epnflag = false; // so that you dont repeat this method
1911 }
1912 }
1913
1914 void G4SPSEneDistribution::ReSetHist(const G4String& atype)
1915 {
1916 G4AutoLock l(&mutex);
1917 if (atype == "energy")
1918 {
1919 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
1920 IPDFEnergyExist = false;
1921 Emin = 0.;
1922 Emax = 1e30;
1923 }
1924 else if (atype == "arb")
1925 {
1926 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
1927 IPDFArbExist = false;
1928 }
1929 else if (atype == "epn")
1930 {
1931 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
1932 IPDFEnergyExist = false;
1933 EpnEnergyH = ZeroPhysVector;
1934 }
1935 else
1936 {
1937 G4cout << "Error, histtype not accepted " << G4endl;
1938 }
1939 }
1940
1941 G4double G4SPSEneDistribution::GenerateOne(G4ParticleDefinition* a)
1942 {
1943 // Copy global shared status to thread-local one
1944 //
1945 threadLocal_t& params = threadLocalData.Get();
1946 params.particle_definition=a;
1947 params.particle_energy=-1;
1948 if(applyEvergyWeight)
1949 {
1950 params.Emax = ArbEmax;
1951 params.Emin = ArbEmin;
1952 }
1953 else
1954 {
1955 params.Emax = Emax;
1956 params.Emin = Emin;
1957 }
1958 params.alpha = alpha;
1959 params.Ezero = Ezero;
1960 params.grad = grad;
1961 params.cept = cept;
1962 params.weight = weight;
1963 // particle_energy = -1.;
1964
1965 if((EnergyDisType == "Mono") && ((MonoEnergy>Emax)||(MonoEnergy<Emin)))
1966 {
1967 G4ExceptionDescription ed;
1968 ed << "MonoEnergy " << G4BestUnit(MonoEnergy,"Energy")
1969 << " is outside of [Emin,Emax] = ["
1970 << G4BestUnit(Emin,"Energy") << ", "
1971 << G4BestUnit(Emax,"Energy") << ". MonoEnergy is used anyway.";
1972 G4Exception("G4SPSEneDistribution::GenerateOne()",
1973 "GPS0001", JustWarning, ed);
1974 params.particle_energy=MonoEnergy;
1975 return params.particle_energy;
1976 }
1977 while ( (EnergyDisType == "Arb")
1978 ? (params.particle_energy < ArbEmin
1979 || params.particle_energy > ArbEmax)
1980 : (params.particle_energy < params.Emin
1981 || params.particle_energy > params.Emax) )
1982 {
1983 if (Biased)
1984 {
1985 GenerateBiasPowEnergies();
1986 }
1987 else
1988 {
1989 if (EnergyDisType == "Mono")
1990 {
1991 GenerateMonoEnergetic();
1992 }
1993 else if (EnergyDisType == "Lin")
1994 {
1995 GenerateLinearEnergies(false);
1996 }
1997 else if (EnergyDisType == "Pow")
1998 {
1999 GeneratePowEnergies(false);
2000 }
2001 else if (EnergyDisType == "CPow")
2002 {
2003 GenerateCPowEnergies();
2004 }
2005 else if (EnergyDisType == "Exp")
2006 {
2007 GenerateExpEnergies(false);
2008 }
2009 else if (EnergyDisType == "Gauss")
2010 {
2011 GenerateGaussEnergies();
2012 }
2013 else if (EnergyDisType == "Brem")
2014 {
2015 GenerateBremEnergies();
2016 }
2017 else if (EnergyDisType == "Bbody")
2018 {
2019 GenerateBbodyEnergies();
2020 }
2021 else if (EnergyDisType == "Cdg")
2022 {
2023 GenerateCdgEnergies();
2024 }
2025 else if (EnergyDisType == "User")
2026 {
2027 GenUserHistEnergies();
2028 }
2029 else if (EnergyDisType == "Arb")
2030 {
2031 GenArbPointEnergies();
2032 }
2033 else if (EnergyDisType == "Epn")
2034 {
2035 GenEpnHistEnergies();
2036 }
2037 else
2038 {
2039 G4cout << "Error: EnergyDisType has unusual value" << G4endl;
2040 }
2041 }
2042 }
2043 return params.particle_energy;
2044 }
2045
2046 G4double G4SPSEneDistribution::GetProbability(G4double ene)
2047 {
2048 G4double prob = 1.;
2049
2050 threadLocal_t& params = threadLocalData.Get();
2051 if (EnergyDisType == "Lin")
2052 {
2053 if (prob_norm == 1.)
2054 {
2055 prob_norm = 0.5*params.grad*params.Emax*params.Emax
2056 + params.cept*params.Emax
2057 - 0.5*params.grad*params.Emin*params.Emin
2058 - params.cept*params.Emin;
2059 }
2060 prob = params.cept + params.grad * ene;
2061 prob /= prob_norm;
2062 }
2063 else if (EnergyDisType == "Pow")
2064 {
2065 if (prob_norm == 1.)
2066 {
2067 if (alpha != -1.)
2068 {
2069 G4double emina = std::pow(params.Emin, params.alpha + 1);
2070 G4double emaxa = std::pow(params.Emax, params.alpha + 1);
2071 prob_norm = 1./(1.+alpha) * (emaxa - emina);
2072 }
2073 else
2074 {
2075 prob_norm = std::log(params.Emax) - std::log(params.Emin) ;
2076 }
2077 }
2078 prob = std::pow(ene, params.alpha)/prob_norm;
2079 }
2080 else if (EnergyDisType == "Exp")
2081 {
2082 if (prob_norm == 1.)
2083 {
2084 prob_norm = -params.Ezero*(std::exp(-params.Emax/params.Ezero)
2085 - std::exp(params.Emin/params.Ezero));
2086 }
2087 prob = std::exp(-ene / params.Ezero);
2088 prob /= prob_norm;
2089 }
2090 else if (EnergyDisType == "Arb")
2091 {
2092 prob = ArbEnergyH.Value(ene);
2093
2094 if (prob <= 0.)
2095 {
2096 G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "
2097 << prob << " " << ene << G4endl;
2098 prob = 1e-30;
2099 }
2100 }
2101 else
2102 {
2103 G4cout << "Error: EnergyDisType not supported" << G4endl;
2104 }
2105
2106 return prob;
2107 }
2108