372 G4Exception(
"G4SPSEneDistribution::ArbEnergyHistoFile",
"Event0301",
FatalException,
"Unable to open the histo ASCII file");
375 while (infile >> ehi >> val)
426 BBHist =
new std::vector<G4double>(10001, 0.0);
427 Bbody_x =
new std::vector<G4double>(10001, 0.0);
434 CPHist =
new std::vector<G4double>(10001, 0.0);
435 CP_x =
new std::vector<G4double>(10001, 0.0);
476 omalpha = 1. - spind[i];
477 CDGhist[i + 1] =
CDGhist[i] + (pfact[i] / omalpha) * (std::pow(ene_line[i + 1] /
keV, omalpha) - std::pow(ene_line[i] /
keV,omalpha));
513 while (count < 10000) {
524 while (count < 10001) {
549 while (count < 10000) {
560 while (count < 10001) {
585 G4Exception(
"G4SPSEneDistribution::ArbInterpolate",
587 "Error: this is for arbitrary distributions");
611 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
614 for (i = 0; i <
maxi; i++)
624 for (count = 0; count < maxi - 1; count++)
626 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
627 / (Arb_x[count + 1] - Arb_x[count]);
639 G4Exception(
"G4SPSEneDistribution::LinearInterpolation",
641 "Error: particle not defined");
651 for (count = 0; count <
maxi; count++)
653 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
656 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
657 Arb_x[count] = total_energy -
mass;
673 Arb_Cum_Area[0] = 0.;
677 Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[i - 1]);
707 Area_seg[i] = ((
Arb_grad[i] / 2) * (Arb_x[i] * Arb_x[i] - Arb_x[i - 1] * Arb_x[i - 1]) +
Arb_cept[i] * (Arb_x[i] - Arb_x[i - 1]));
708 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
709 sum = sum + Area_seg[i];
720 Arb_Cum_Area[i] = Arb_Cum_Area[i] /
sum;
745 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
748 for (i = 0; i <
maxi; i++)
758 for (count = 0; count < maxi - 1; count++)
760 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) / (Arb_x[count + 1] - Arb_x[count]);
772 G4Exception(
"G4SPSEneDistribution::LogInterpolation",
774 "Error: particle not defined");
784 for (count = 0; count <
maxi; count++)
786 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
789 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
790 Arb_x[count] = total_energy -
mass;
808 Arb_Cum_Area[0] = 0.;
809 if (Arb_x[0] <= 0. || Arb_y[0] <= 0.)
811 G4cout <<
"You should not use log interpolation with points <= 0." <<
G4endl;
812 G4cout <<
"These will be changed to 1e-20, which may cause problems" <<
G4endl;
827 if (Arb_x[i] <= 0. || Arb_y[i] <= 0.)
829 G4cout <<
"You should not use log interpolation with points <= 0." <<
G4endl;
830 G4cout <<
"These will be changed to 1e-20, which may cause problems" <<
G4endl;
841 Arb_alpha[i] = (std::log10(Arb_y[i]) - std::log10(Arb_y[i - 1])) / (std::log10(Arb_x[i]) - std::log10(Arb_x[i - 1]));
846 Area_seg[i] =
Arb_Const[i] * (std::log(Arb_x[i]) - std::log(Arb_x[i - 1]));
850 Area_seg[i] = (
Arb_Const[i] / alp) * (std::pow(Arb_x[i], alp) - std::pow(Arb_x[i - 1], alp));
852 sum = sum + Area_seg[i];
853 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
864 Arb_Cum_Area[i] = Arb_Cum_Area[i] /
sum;
887 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
890 for (i = 0; i <
maxi; i++)
900 for (count = 0; count < maxi - 1; count++)
902 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
903 / (Arb_x[count + 1] - Arb_x[count]);
915 G4Exception(
"G4SPSEneDistribution::ExpInterpolation",
917 "Error: particle not defined");
927 for (count = 0; count <
maxi; count++)
929 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass * mass));
930 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
931 Arb_x[count] = total_energy -
mass;
948 Arb_Cum_Area[0] = 0.;
951 G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i - 1]);
952 if (test > 0. || test < 0.)
954 Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1]) / (std::log(Arb_y[i]) - std::log(Arb_y[i - 1]));
960 G4Exception(
"G4SPSEneDistribution::ExpInterpolation",
962 "Flat line segment: problem, setting to zero parameters.");
968 sum = sum + Area_seg[i];
969 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
980 Arb_Cum_Area[i] = Arb_Cum_Area[i] /
sum;
1004 G4double sum, Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
1008 for (i = 0; i <
maxi; i++) {
1016 for (count = 0; count < maxi - 1; count++)
1018 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) / (Arb_x[count + 1] - Arb_x[count]);
1028 G4Exception(
"G4SPSEneDistribution::SplineInterpolation",
1030 "Error: particle not defined");
1038 for (count = 0; count <
maxi; count++) {
1039 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
1042 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
1043 Arb_x[count] = total_energy -
mass;
1050 Arb_Cum_Area[0] = 0.;
1063 G4double de = (Arb_x[i] - Arb_x[i - 1])/100.;
1066 for (count = 0; count < 101; count++) {
1067 ei[count] = Arb_x[i - 1] + de*count ;
1069 if (prob[count] < 0.) {
1071 ED <<
"Warning: G4DataInterpolation returns value < 0 " << prob[count] <<
" "<<ei[count]<<
G4endl;
1072 G4Exception(
"G4SPSEneDistribution::SplineInterpolation",
"Event0303",
1075 area += prob[count]*de;
1077 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + area;
1080 prob[0] = prob[0]/(area/de);
1081 for (count = 1; count < 100; count++)
1082 prob[count] = prob[count-1] + prob[count]/(area/de);
1090 Arb_Cum_Area[i] = Arb_Cum_Area[i] /
sum;
1109 if (ene < 0) ene = 0.;
1126 bracket = bracket * rndm;
1127 bracket = bracket + (params.
grad / 2.) * eminsq + params.
cept * params.
Emin;
1131 if (params.
grad != 0.)
1133 G4double sqbrack = (intersq - 4 * (params.
grad / 2.) * (bracket));
1135 sqbrack = std::sqrt(sqbrack);
1137 root1 = root1 / (2. * (params.
grad / 2.));
1140 root2 = root2 / (2. * (params.
grad / 2.));
1144 if (root1 > params.
Emin && root1 < params.
Emax)
1148 if (root2 > params.
Emin && root2 < params.
Emax)
1153 else if (params.
grad == 0.)
1180 emina = std::pow(params.
Emin, params.
alpha + 1);
1181 emaxa = std::pow(params.
Emax, params.
alpha + 1);
1192 if (params.
alpha != -1.)
1194 G4double ene = ((rndm * (emaxa - emina)) + emina);
1195 ene = std::pow(ene, (1. / (params.
alpha + 1.)));
1200 G4double ene = (std::log(params.
Emin) + rndm * (std::log(params.
Emax) - std::log(params.
Emin)));
1219 G4int nabove, nbelow = 0, middle;
1234 while (nabove - nbelow > 1)
1236 middle = (nabove + nbelow) / 2;
1237 if (rndm ==
CPHist->at(middle))
1241 if (rndm < CPHist->
at(middle))
1253 x1 =
CP_x->at(nbelow);
1254 if(nbelow+1 == static_cast<G4int>(
CP_x->size()))
1260 x2 =
CP_x->at(nbelow + 1);
1263 if(nbelow+1 == static_cast<G4int>(
CPHist->size()))
1270 y2 =
CPHist->at(nbelow + 1);
1272 t = (y2 -
y1) / (x2 - x1);
1307 G4double ee = ((rndm * (emaxa - emina)) + emina);
1309 normal = 1./(1+
biasalpha) * (emaxa - emina);
1313 G4double ee = (std::log(emin) + rndm * (std::log(emax) - std::log(emin)));
1315 normal = std::log(emax) - std::log(emin) ;
1365 expmax = std::exp(-params.
Emax / (k *
Temp));
1366 expmin = std::exp(-params.
Emin / (k * Temp));
1373 G4Exception(
"G4SPSEneDistribution::GenerateBremEnergies",
1375 "*****EXPMAX=0. Choose different E's or Temp");
1379 G4Exception(
"G4SPSEneDistribution::GenerateBremEnergies",
1381 "*****EXPMIN=0. Choose different E's or Temp");
1384 G4double tempvar = rndm * ((-
k) * Temp * (params.
Emax * expmax - params.
Emin * expmin) - (ksq * Tsq * (expmax - expmin)));
1386 G4double bigc = (tempvar - k * Temp * params.
Emin * expmin - ksq * Tsq * expmin) / (-k * Temp);
1399 for (i = 1; i < 1000; i++)
1401 etest = params.
Emin + (i - 1) * steps;
1403 diff = etest * (std::exp(-etest / (k * Temp))) + k * Temp * (std::exp(-etest / (k * Temp))) - bigc;
1429 G4int nabove, nbelow = 0, middle;
1444 while (nabove - nbelow > 1)
1446 middle = (nabove + nbelow) / 2;
1447 if (rndm ==
BBHist->at(middle))
1451 if (rndm < BBHist->
at(middle))
1464 if(nbelow+1 == static_cast<G4int>(
Bbody_x->size()))
1473 if(nbelow+1 == static_cast<G4int>(
BBHist->size()))
1480 y2 =
BBHist->at(nbelow + 1);
1482 t = (y2 -
y1) / (x2 - x1);
1505 omalpha[0] = 1. - 1.4;
1506 ene_line[0] = params.
Emin;
1507 ene_line[1] = params.
Emax;
1511 omalpha[0] = 1. - 1.4;
1512 omalpha[1] = 1. - 2.3;
1513 ene_line[0] = params.
Emin;
1514 ene_line[1] = 18. *
keV;
1515 ene_line[2] = params.
Emax;
1519 omalpha[0] = 1. - 2.3;
1520 ene_line[0] = params.
Emin;
1521 ene_line[1] = params.
Emax;
1532 G4double ene = (std::pow(ene_line[i - 1], omalpha[i - 1]) + (std::pow(ene_line[i], omalpha[i - 1]) - std::pow(ene_line[i - 1], omalpha[i- 1])) * rndm2);
1551 for ( ii = 0 ; ii<1024 ; ++ii ) { bins[ii]=0; vals[ii]=0; }
1556 G4Exception(
"G4SPSEneDistribution::GenUserHistEnergies",
1558 "Error: particle definition is NULL");
1563 G4Exception(
"G4SPSEneDistribution::GenUserHistEnergies",
1564 "Event0302",
JustWarning,
"Maxbin>1024\n Setting maxbin to 1024, other bins are lost");
1577 for (ii = 1; ii < maxbin; ii++)
1580 vals[ii] =
UDefEnergyH(
size_t(ii)) + vals[ii - 1];
1591 for (ii = 1; ii < maxbin; ii++)
1593 vals[ii] = vals[ii] * (bins[ii] - bins[ii - 1]);
1597 for (ii = 0; ii < maxbin; ii++)
1599 bins[ii] = std::sqrt((bins[ii] * bins[ii]) + (mass * mass))-
mass;
1601 for (ii = 1; ii < maxbin; ii++)
1603 vals[ii] = vals[ii] / (bins[ii] - bins[ii - 1]);
1605 sum = vals[maxbin - 1];
1608 for (ii = 0; ii < maxbin; ii++)
1610 vals[ii] = vals[ii] /
sum;
1644 G4int nabove, nbelow = 0, middle;
1648 while (nabove - nbelow > 1)
1650 middle = (nabove + nbelow) / 2;
1709 G4Exception(
"G4SPSEneDistribution::GenArbPointEnergies",
"Event0302",
1736 for (ii = 1; ii < maxbin; ii++)
1739 vals[ii] =
UDefEnergyH(
size_t(ii)) + vals[ii - 1];
1744 for (ii = 0; ii < maxbin; ii++)
1746 vals[ii] = vals[ii] /
sum;
1784 G4int count, maxcount;
1788 if (maxcount > 1024)
1791 "Histogram contains more than 1024 bins!\nThose above 1024 will be ignored");
1796 G4Exception(
"G4SPSEneDistribution::ConvertEPNToEnergy()",
"gps001",
FatalException,
"Histogram contains less than 1 bin!\nRedefine the histogram");
1799 for (count = 0; count < maxcount; count++)
1807 for (count = 0; count < maxcount; count++)
1809 ebins[count] = ebins[count] * Bary;
1813 params.
Emin = ebins[0];
1816 params.
Emax = ebins[maxcount - 1];
1820 params.
Emax = ebins[0];
1823 for (count = 0; count < maxcount; count++)
1835 if (atype ==
"energy")
1842 else if (atype ==
"arb")
1847 else if (atype ==
"epn")
1940 G4cout <<
"Error: EnergyDisType has unusual value" <<
G4endl;
1958 prob = params.
cept + params.
grad * ene;
1984 prob = std::exp(-ene / params.
Ezero);
1996 G4cout <<
" Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<
" "<<ene <<
G4endl;
2003 G4cout <<
"Error: EnergyDisType not supported" <<
G4endl;