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