54 : particle_definition(0), eneRndm(0), Splinetemp(0)
169 std::ifstream infile(filename, std::ios::in);
171 G4Exception(
"G4SPSEneDistribution::ArbEnergyHistoFile",
173 "Unable to open the histo ASCII file");
175 while (infile >> ehi >> val) {
229 omalpha = 1. - spind[i];
231 ene_line[i + 1] /
keV, omalpha) - std::pow(ene_line[i] /
keV,
263 while (count < 10000) {
265 Bbody_y[count] = (2. * std::pow(
Bbody_x[count], 2.)) / (h2 * c2
267 sum = sum + Bbody_y[count];
275 while (count < 10001) {
297 G4cout <<
"Error: this is for arbitrary distributions" <<
G4endl;
320 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
323 for (i = 0; i < maxi; i++) {
331 for (count = 0; count < maxi - 1; count++) {
332 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
333 / (Arb_x[count + 1] - Arb_x[count]);
350 for (count = 0; count < maxi; count++) {
351 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
354 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
355 Arb_x[count] = total_energy - mass;
364 Arb_Cum_Area[0] = 0.;
367 Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[i - 1]);
384 Area_seg[i] = ((
Arb_grad[i] / 2) * (Arb_x[i] * Arb_x[i] - Arb_x[i - 1]
385 * Arb_x[i - 1]) +
Arb_cept[i] * (Arb_x[i] - Arb_x[i - 1]));
386 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
387 sum = sum + Area_seg[i];
389 G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum <<
Arb_grad[i]
396 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
418 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
421 for (i = 0; i < maxi; i++) {
429 for (count = 0; count < maxi - 1; count++) {
430 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
431 / (Arb_x[count + 1] - Arb_x[count]);
448 for (count = 0; count < maxi; count++) {
449 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
452 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
453 Arb_x[count] = total_energy - mass;
462 Arb_Cum_Area[0] = 0.;
463 if (Arb_x[0] <= 0. || Arb_y[0] <= 0.) {
464 G4cout <<
"You should not use log interpolation with points <= 0."
466 G4cout <<
"These will be changed to 1e-20, which may cause problems"
477 if (Arb_x[i] <= 0. || Arb_y[i] <= 0.) {
478 G4cout <<
"You should not use log interpolation with points <= 0."
481 <<
"These will be changed to 1e-20, which may cause problems"
489 Arb_alpha[i] = (std::log10(Arb_y[i]) - std::log10(Arb_y[i - 1]))
490 / (std::log10(Arb_x[i]) - std::log10(Arb_x[i - 1]));
494 Area_seg[i] =
Arb_Const[i] * (std::log(Arb_x[i]) - std::log(Arb_x[i - 1]));
496 Area_seg[i] = (
Arb_Const[i] / alp) * (std::pow(Arb_x[i], alp)
497 - std::pow(Arb_x[i - 1], alp));
499 sum = sum + Area_seg[i];
500 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
508 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
527 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
530 for (i = 0; i < maxi; i++) {
538 for (count = 0; count < maxi - 1; count++) {
539 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
540 / (Arb_x[count + 1] - Arb_x[count]);
557 for (count = 0; count < maxi; count++) {
558 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
561 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
562 Arb_x[count] = total_energy - mass;
571 Arb_Cum_Area[0] = 0.;
573 G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i - 1]);
574 if (test > 0. || test < 0.) {
575 Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1]) / (std::log(Arb_y[i])
576 - std::log(Arb_y[i - 1]));
586 sum = sum + Area_seg[i];
587 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
595 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
615 G4double sum, Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
619 for (i = 0; i < maxi; i++) {
627 for (count = 0; count < maxi - 1; count++) {
628 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
629 / (Arb_x[count + 1] - Arb_x[count]);
638 G4Exception(
"G4SPSEneDistribution::SplineInterpolation",
640 "Error: particle not defined");
648 for (count = 0; count < maxi; count++) {
649 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
652 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
653 Arb_x[count] = total_energy - mass;
660 Arb_Cum_Area[0] = 0.;
666 G4double de = (Arb_x[i] - Arb_x[i - 1])/100.;
669 for (count = 0; count < 101; count++) {
670 ei[count] = Arb_x[i - 1] + de*count ;
672 if (prob[count] < 0.) {
674 ED <<
"Warning: G4DataInterpolation returns value < 0 " << prob[count] <<
" "<<ei[count]<<
G4endl;
675 G4Exception(
"G4SPSEneDistribution::SplineInterpolation",
"Event0303",
678 area += prob[count]*de;
680 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + area;
683 prob[0] = prob[0]/(area/de);
684 for (count = 1; count < 100; count++)
685 prob[count] = prob[count-1] + prob[count]/(area/de);
693 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
727 bracket = bracket * rndm;
733 G4double sqbrack = (intersq - 4 * (
grad / 2.) * (bracket));
735 sqbrack = std::sqrt(sqbrack);
737 root1 = root1 / (2. * (
grad / 2.));
740 root2 = root2 / (2. * (
grad / 2.));
744 if (root1 > Emin && root1 <
Emax)
746 if (root2 > Emin && root2 <
Emax)
748 }
else if (
grad == 0.)
809 normal = 1./(1+
biasalpha) * (emaxa - emina);
814 normal = std::log(emax) - std::log(emin) ;
851 expmax = std::exp(-
Emax / (k *
Temp));
852 expmin = std::exp(-
Emin / (k * Temp));
858 G4cout <<
"*****EXPMAX=0. Choose different E's or Temp" <<
G4endl;
860 G4cout <<
"*****EXPMIN=0. Choose different E's or Temp" <<
G4endl;
863 - (ksq * Tsq * (expmax - expmin)));
865 G4double bigc = (tempvar - k * Temp *
Emin * expmin - ksq * Tsq * expmin)
879 for (i = 1; i < 1000; i++) {
880 etest = Emin + (i - 1) * steps;
882 diff = etest * (std::exp(-etest / (k * Temp))) + k * Temp * (std::exp(
883 -etest / (k * Temp))) - bigc;
903 G4int nabove, nbelow = 0, middle;
908 while (nabove - nbelow > 1) {
909 middle = (nabove + nbelow) / 2;
910 if (rndm ==
BBHist[middle])
912 if (rndm <
BBHist[middle])
924 t = (y2 - y1) / (x2 - x1);
945 omalpha[0] = 1. - 1.4;
949 if (Emin < 18 * keV && Emax > 18 *
keV) {
950 omalpha[0] = 1. - 1.4;
951 omalpha[1] = 1. - 2.3;
953 ene_line[1] = 18. *
keV;
956 if (
Emin > 18 * keV) {
957 omalpha[0] = 1. - 2.3;
969 particle_energy = (std::pow(ene_line[i - 1], omalpha[i - 1]) + (std::pow(
970 ene_line[i], omalpha[i - 1]) - std::pow(ene_line[i - 1], omalpha[i
984 G4double bins[1024], vals[1024], sum;
988 G4cout <<
"Error: particle definition is NULL" <<
G4endl;
992 G4cout <<
"Setting maxbin to 1024, other bins are lost" <<
G4endl;
1001 for (ii = 1; ii < maxbin; ii++) {
1003 vals[ii] =
UDefEnergyH(
size_t(ii)) + vals[ii - 1];
1013 for (ii = 1; ii < maxbin; ii++) {
1014 vals[ii] = vals[ii] * (bins[ii] - bins[ii - 1]);
1018 for (ii = 0; ii < maxbin; ii++) {
1019 bins[ii] = std::sqrt((bins[ii] * bins[ii]) + (mass * mass))
1022 for (ii = 1; ii < maxbin; ii++) {
1023 vals[ii] = vals[ii] / (bins[ii] - bins[ii - 1]);
1025 sum = vals[maxbin - 1];
1028 for (ii = 0; ii < maxbin; ii++) {
1029 vals[ii] = vals[ii] / sum;
1055 G4int nabove, nbelow = 0, middle;
1059 while (nabove - nbelow > 1) {
1060 middle = (nabove + nbelow) / 2;
1075 }
else if (
IntType ==
"Log") {
1081 }
else if (
IntType ==
"Exp") {
1087 }
else if (
IntType ==
"Spline") {
1092 while (particle_energy < Emin || particle_energy >
Emax) {
1118 G4double bins[1024], vals[1024], sum;
1124 for (ii = 1; ii < maxbin; ii++) {
1126 vals[ii] =
UDefEnergyH(
size_t(ii)) + vals[ii - 1];
1130 for (ii = 0; ii < maxbin; ii++) {
1131 vals[ii] = vals[ii] / sum;
1159 G4int count, maxcount;
1163 if (maxcount > 1024) {
1164 G4cout <<
"Histogram contains more than 1024 bins!" <<
G4endl;
1165 G4cout <<
"Those above 1024 will be ignored" <<
G4endl;
1169 G4cout <<
"Histogram contains less than 1 bin!" <<
G4endl;
1173 for (count = 0; count < maxcount; count++) {
1180 for (count = 0; count < maxcount; count++) {
1181 ebins[count] = ebins[count] * Bary;
1187 Emax = ebins[maxcount - 1];
1191 for (count = 0; count < maxcount; count++) {
1200 if (atype ==
"energy") {
1205 }
else if (atype ==
"arb") {
1208 }
else if (atype ==
"epn") {
1250 G4cout <<
"Error: EnergyDisType has unusual value" <<
G4endl;
1282 prob = std::exp(-ene /
Ezero);
1292 G4cout <<
" Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<
" "<<ene <<
G4endl;
1298 G4cout <<
"Error: EnergyDisType not supported" <<
G4endl;
void SetBeamSigmaInE(G4double)
ThreeVector shoot(const G4int Ap, const G4int Af)
static c2_factory< G4double > c2
void CalculateBbodySpectrum()
void LinearInterpolation()
void GenEpnHistEnergies()
void CalculateCdgSpectrum()
std::ostringstream G4ExceptionDescription
CLHEP::Hep3Vector G4ThreeVector
void ArbEnergyHisto(G4ThreeVector)
void ConvertEPNToEnergy()
void GenerateCdgEnergies()
void InsertValues(G4double energy, G4double value)
void GenArbPointEnergies()
size_t GetVectorLength() const
G4double GetLowEdgeEnergy(size_t binNumber) const
G4double CubicSplineInterpolation(G4double pX) const
void GenerateGaussEnergies()
G4PhysicsOrderedFreeVector IPDFEnergyH
G4double GetMaxLowEdgeEnergy()
static double normal(HepRandomEngine *eptr)
G4ParticleDefinition * particle_definition
G4GLOB_DLL std::ostream G4cout
void SetEnergyDisType(G4String)
void GenerateBremEnergies()
G4double GetProbability(G4double)
G4DataInterpolation * Splinetemp
void GenerateMonoEnergetic()
G4double GenerateOne(G4ParticleDefinition *)
virtual void ScaleVector(G4double factorE, G4double factorV)
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)
void ArbInterpolate(G4String)
void EpnEnergyHisto(G4ThreeVector)
G4PhysicsOrderedFreeVector ArbEnergyH
void InputEnergySpectra(G4bool)
G4SPSRandomGenerator * eneRndm
G4double GetPDGMass() const
G4PhysicsOrderedFreeVector UDefEnergyH
G4DataInterpolation * SplineInt[1024]
void SetGradient(G4double)
void UserEnergyHisto(G4ThreeVector)
void GenUserHistEnergies()
void GenerateBbodyEnergies()
void GenerateExpEnergies(G4bool)
void SetInterCept(G4double)
void InputDifferentialSpectra(G4bool)
G4PhysicsOrderedFreeVector EpnEnergyH
void GeneratePowEnergies(G4bool)
void SetBiasAlpha(G4double)
void ArbEnergyHistoFile(G4String)
G4double GetMinLowEdgeEnergy()
G4PhysicsOrderedFreeVector IPDFArbEnergyH
void GenerateBiasPowEnergies()
void GenerateLinearEnergies(G4bool)
void SplineInterpolation()
void SetMonoEnergy(G4double)
G4int GetBaryonNumber() const
G4PhysicsOrderedFreeVector ZeroPhysVector