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