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