150 for ( std::vector<G4DataInterpolation*>::iterator it =
SplineInt.begin() ; it!=
SplineInt.end() ; ++it)
371 G4Exception(
"G4SPSEneDistribution::ArbEnergyHistoFile",
"Event0301",
FatalException,
"Unable to open the histo ASCII file");
374 while (infile >> ehi >> val)
425 BBHist =
new std::vector<G4double>(10001, 0.0);
426 Bbody_x =
new std::vector<G4double>(10001, 0.0);
433 CPHist =
new std::vector<G4double>(10001, 0.0);
434 CP_x =
new std::vector<G4double>(10001, 0.0);
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));
512 while (count < 10000) {
523 while (count < 10001) {
548 while (count < 10000) {
559 while (count < 10001) {
584 G4Exception(
"G4SPSEneDistribution::ArbInterpolate",
586 "Error: this is for arbitrary distributions");
610 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
613 for (i = 0; i <
maxi; i++)
623 for (count = 0; count < maxi - 1; count++)
625 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
626 / (Arb_x[count + 1] - Arb_x[count]);
638 G4Exception(
"G4SPSEneDistribution::LinearInterpolation",
640 "Error: particle not defined");
650 for (count = 0; count <
maxi; count++)
652 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
655 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
656 Arb_x[count] = total_energy - mass;
672 Arb_Cum_Area[0] = 0.;
676 Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[i - 1]);
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];
719 Arb_Cum_Area[i] = Arb_Cum_Area[i] /
sum;
744 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
747 for (i = 0; i <
maxi; i++)
757 for (count = 0; count < maxi - 1; count++)
759 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) / (Arb_x[count + 1] - Arb_x[count]);
771 G4Exception(
"G4SPSEneDistribution::LogInterpolation",
773 "Error: particle not defined");
783 for (count = 0; count <
maxi; count++)
785 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
788 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
789 Arb_x[count] = total_energy - mass;
807 Arb_Cum_Area[0] = 0.;
808 if (Arb_x[0] <= 0. || Arb_y[0] <= 0.)
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;
826 if (Arb_x[i] <= 0. || Arb_y[i] <= 0.)
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;
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]));
845 Area_seg[i] =
Arb_Const[i] * (std::log(Arb_x[i]) - std::log(Arb_x[i - 1]));
849 Area_seg[i] = (
Arb_Const[i] / alp) * (std::pow(Arb_x[i], alp) - std::pow(Arb_x[i - 1], alp));
851 sum = sum + Area_seg[i];
852 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
863 Arb_Cum_Area[i] = Arb_Cum_Area[i] /
sum;
886 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
889 for (i = 0; i <
maxi; i++)
899 for (count = 0; count < maxi - 1; count++)
901 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
902 / (Arb_x[count + 1] - Arb_x[count]);
914 G4Exception(
"G4SPSEneDistribution::ExpInterpolation",
916 "Error: particle not defined");
926 for (count = 0; count <
maxi; count++)
928 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass * mass));
929 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
930 Arb_x[count] = total_energy - mass;
947 Arb_Cum_Area[0] = 0.;
950 G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i - 1]);
951 if (test > 0. || test < 0.)
953 Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1]) / (std::log(Arb_y[i]) - std::log(Arb_y[i - 1]));
959 G4Exception(
"G4SPSEneDistribution::ExpInterpolation",
961 "Flat line segment: problem, setting to zero parameters.");
967 sum = sum + Area_seg[i];
968 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
979 Arb_Cum_Area[i] = Arb_Cum_Area[i] /
sum;
1003 G4double sum, Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
1007 for (i = 0; i <
maxi; i++) {
1015 for (count = 0; count < maxi - 1; count++)
1017 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) / (Arb_x[count + 1] - Arb_x[count]);
1027 G4Exception(
"G4SPSEneDistribution::SplineInterpolation",
1029 "Error: particle not defined");
1037 for (count = 0; count <
maxi; count++) {
1038 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
1041 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
1042 Arb_x[count] = total_energy - mass;
1049 Arb_Cum_Area[0] = 0.;
1053 for ( std::vector<G4DataInterpolation*>::iterator it =
SplineInt.begin() ; it!=
SplineInt.end() ; ++it)
1062 G4double de = (Arb_x[i] - Arb_x[i - 1])/100.;
1065 for (count = 0; count < 101; count++) {
1066 ei[count] = Arb_x[i - 1] + de*count ;
1068 if (prob[count] < 0.) {
1070 ED <<
"Warning: G4DataInterpolation returns value < 0 " << prob[count] <<
" "<<ei[count]<<
G4endl;
1071 G4Exception(
"G4SPSEneDistribution::SplineInterpolation",
"Event0303",
1074 area += prob[count]*de;
1076 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + area;
1079 prob[0] = prob[0]/(area/de);
1080 for (count = 1; count < 100; count++)
1081 prob[count] = prob[count-1] + prob[count]/(area/de);
1089 Arb_Cum_Area[i] = Arb_Cum_Area[i] /
sum;
1108 if (ene < 0) ene = 0.;
1125 bracket = bracket * rndm;
1126 bracket = bracket + (params.
grad / 2.) * eminsq + params.
cept * params.
Emin;
1130 if (params.
grad != 0.)
1132 G4double sqbrack = (intersq - 4 * (params.
grad / 2.) * (bracket));
1134 sqbrack = std::sqrt(sqbrack);
1136 root1 = root1 / (2. * (params.
grad / 2.));
1139 root2 = root2 / (2. * (params.
grad / 2.));
1143 if (root1 > params.
Emin && root1 < params.
Emax)
1147 if (root2 > params.
Emin && root2 < params.
Emax)
1152 else if (params.
grad == 0.)
1179 emina = std::pow(params.
Emin, params.
alpha + 1);
1180 emaxa = std::pow(params.
Emax, params.
alpha + 1);
1191 if (params.
alpha != -1.)
1193 G4double ene = ((rndm * (emaxa - emina)) + emina);
1194 ene = std::pow(ene, (1. / (params.
alpha + 1.)));
1199 G4double ene = (std::log(params.
Emin) + rndm * (std::log(params.
Emax) - std::log(params.
Emin)));
1218 G4int nabove, nbelow = 0, middle;
1233 while (nabove - nbelow > 1)
1235 middle = (nabove + nbelow) / 2;
1236 if (rndm ==
CPHist->at(middle))
1240 if (rndm < CPHist->at(middle))
1252 x1 =
CP_x->at(nbelow);
1253 if(nbelow+1 == static_cast<G4int>(
CP_x->size()))
1259 x2 =
CP_x->at(nbelow + 1);
1262 if(nbelow+1 == static_cast<G4int>(
CPHist->size()))
1269 y2 =
CPHist->at(nbelow + 1);
1271 t = (y2 -
y1) / (x2 - x1);
1306 G4double ee = ((rndm * (emaxa - emina)) + emina);
1308 normal = 1./(1+
biasalpha) * (emaxa - emina);
1312 G4double ee = (std::log(emin) + rndm * (std::log(emax) - std::log(emin)));
1314 normal = std::log(emax) - std::log(emin) ;
1364 expmax = std::exp(-params.
Emax / (k *
Temp));
1365 expmin = std::exp(-params.
Emin / (k * Temp));
1372 G4Exception(
"G4SPSEneDistribution::GenerateBremEnergies",
1374 "*****EXPMAX=0. Choose different E's or Temp");
1378 G4Exception(
"G4SPSEneDistribution::GenerateBremEnergies",
1380 "*****EXPMIN=0. Choose different E's or Temp");
1383 G4double tempvar = rndm * ((-k) * Temp * (params.
Emax * expmax - params.
Emin * expmin) - (ksq * Tsq * (expmax - expmin)));
1385 G4double bigc = (tempvar - k * Temp * params.
Emin * expmin - ksq * Tsq * expmin) / (-k * Temp);
1398 for (i = 1; i < 1000; i++)
1400 etest = params.
Emin + (i - 1) * steps;
1402 diff = etest * (std::exp(-etest / (k * Temp))) + k * Temp * (std::exp(-etest / (k * Temp))) - bigc;
1428 G4int nabove, nbelow = 0, middle;
1443 while (nabove - nbelow > 1)
1445 middle = (nabove + nbelow) / 2;
1446 if (rndm ==
BBHist->at(middle))
1450 if (rndm < BBHist->at(middle))
1463 if(nbelow+1 == static_cast<G4int>(
Bbody_x->size()))
1472 if(nbelow+1 == static_cast<G4int>(
BBHist->size()))
1479 y2 =
BBHist->at(nbelow + 1);
1481 t = (y2 -
y1) / (x2 - x1);
1504 omalpha[0] = 1. - 1.4;
1505 ene_line[0] = params.
Emin;
1506 ene_line[1] = params.
Emax;
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;
1518 omalpha[0] = 1. - 2.3;
1519 ene_line[0] = params.
Emin;
1520 ene_line[1] = params.
Emax;
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);
1550 for ( ii = 0 ; ii<1024 ; ++ii ) { bins[ii]=0; vals[ii]=0; }
1555 G4Exception(
"G4SPSEneDistribution::GenUserHistEnergies",
1557 "Error: particle definition is NULL");
1562 G4Exception(
"G4SPSEneDistribution::GenUserHistEnergies",
1563 "Event0302",
JustWarning,
"Maxbin>1024\n Setting maxbin to 1024, other bins are lost");
1576 for (ii = 1; ii < maxbin; ii++)
1579 vals[ii] =
UDefEnergyH(
size_t(ii)) + vals[ii - 1];
1590 for (ii = 1; ii < maxbin; ii++)
1592 vals[ii] = vals[ii] * (bins[ii] - bins[ii - 1]);
1596 for (ii = 0; ii < maxbin; ii++)
1598 bins[ii] = std::sqrt((bins[ii] * bins[ii]) + (mass * mass))- mass;
1600 for (ii = 1; ii < maxbin; ii++)
1602 vals[ii] = vals[ii] / (bins[ii] - bins[ii - 1]);
1604 sum = vals[maxbin - 1];
1607 for (ii = 0; ii < maxbin; ii++)
1609 vals[ii] = vals[ii] /
sum;
1643 G4int nabove, nbelow = 0, middle;
1647 while (nabove - nbelow > 1)
1649 middle = (nabove + nbelow) / 2;
1708 G4Exception(
"G4SPSEneDistribution::GenArbPointEnergies",
"Event0302",
1735 for (ii = 1; ii < maxbin; ii++)
1738 vals[ii] =
UDefEnergyH(
size_t(ii)) + vals[ii - 1];
1743 for (ii = 0; ii < maxbin; ii++)
1745 vals[ii] = vals[ii] /
sum;
1783 G4int count, maxcount;
1787 if (maxcount > 1024)
1790 "Histogram contains more than 1024 bins!\nThose above 1024 will be ignored");
1795 G4Exception(
"G4SPSEneDistribution::ConvertEPNToEnergy()",
"gps001",
FatalException,
"Histogram contains less than 1 bin!\nRedefine the histogram");
1798 for (count = 0; count < maxcount; count++)
1806 for (count = 0; count < maxcount; count++)
1808 ebins[count] = ebins[count] * Bary;
1812 params.
Emin = ebins[0];
1815 params.
Emax = ebins[maxcount - 1];
1819 params.
Emax = ebins[0];
1822 for (count = 0; count < maxcount; count++)
1834 if (atype ==
"energy")
1841 else if (atype ==
"arb")
1846 else if (atype ==
"epn")
1930 G4cout <<
"Error: EnergyDisType has unusual value" <<
G4endl;
1948 prob = params.
cept + params.
grad * ene;
1974 prob = std::exp(-ene / params.
Ezero);
1986 G4cout <<
" Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<
" "<<ene <<
G4endl;
1993 G4cout <<
"Error: EnergyDisType not supported" <<
G4endl;
G4double GetProbability(G4double)
#define G4MUTEXDESTROY(mutex)
std::vector< G4double > * CPHist
G4PhysicsOrderedFreeVector IPDFEnergyH
std::ostringstream G4ExceptionDescription
std::vector< ExP01TrackerHit * > a
void SetGradient(G4double)
void ArbEnergyHisto(G4ThreeVector)
void SetEnergyDisType(G4String)
void GenUserHistEnergies()
static constexpr double MeV
G4double GetLowEdgeEnergy(size_t binNumber) const
Float_t y1[n_points_granero]
void InsertValues(G4double energy, G4double value)
void GenerateCdgEnergies()
static constexpr double keV
void EpnEnergyHisto(G4ThreeVector)
void GenArbPointEnergies()
void GenerateMonoEnergetic()
void ArbInterpolate(G4String)
Float_t x1[n_points_granero]
void SetBiasAlpha(G4double)
virtual void ScaleVector(G4double factorE, G4double factorV)
G4double GetEnergy(G4double aValue)
void GeneratePowEnergies(G4bool)
void ArbEnergyHistoFile(G4String)
std::vector< G4DataInterpolation * > SplineInt
G4double GenerateOne(G4ParticleDefinition *)
void SetBiasRndm(G4SPSRandomGenerator *a)
void LinearInterpolation()
void GenerateCPowEnergies()
static double normal(HepRandomEngine *eptr)
std::vector< G4double > * CP_x
std::vector< G4double > * BBHist
G4PhysicsOrderedFreeVector UDefEnergyH
static const G4double emax
G4double Value(G4double theEnergy, size_t &lastidx) const
G4bool Arb_grad_cept_flag
G4double GetPDGMass() const
void GenerateExpEnergies(G4bool)
G4double GetMaxLowEdgeEnergy()
void GenerateGaussEnergies()
Float_t y2[n_points_geant4]
void SetBeamSigmaInE(G4double)
const XML_Char const XML_Char * data
G4PhysicsOrderedFreeVector ZeroPhysVector
void DumpValues(G4double unitE=1.0, G4double unitV=1.0) const
G4double GetMinLowEdgeEnergy()
void ConvertEPNToEnergy()
#define G4MUTEXINIT(mutex)
const XML_Char int const XML_Char * value
std::vector< G4double > * Bbody_x
G4int GetBaryonNumber() const
G4PhysicsOrderedFreeVector IPDFArbEnergyH
G4PhysicsOrderedFreeVector EpnEnergyH
void SetVerbosity(G4int a)
void CalculateCdgSpectrum()
ThreeVector shoot(const G4int Ap, const G4int Af)
void SplineInterpolation()
void GenerateBiasPowEnergies()
void InputEnergySpectra(G4bool)
void UserEnergyHisto(G4ThreeVector)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
void GenerateBbodyEnergies()
void GenEpnHistEnergies()
G4double CubicSplineInterpolation(G4double pX) const
void CalculateBbodySpectrum()
G4PhysicsOrderedFreeVector ArbEnergyH
void SetInterCept(G4double)
void SetMonoEnergy(G4double)
G4ParticleDefinition * particle_definition
G4bool Arb_alpha_Const_flag
G4SPSRandomGenerator * eneRndm
G4Cache< threadLocal_t > threadLocalData
G4GLOB_DLL std::ostream G4cout
void CalculateCPowSpectrum()
void GenerateBremEnergies()
G4String GetEnergyDisType()
G4PhysicsOrderedFreeVector GetUserDefinedEnergyHisto()
Float_t x2[n_points_geant4]
G4DataInterpolation * Splinetemp
void InputDifferentialSpectra(G4bool)
void GenerateLinearEnergies(G4bool)
G4PhysicsOrderedFreeVector GetArbEnergyHisto()
size_t GetVectorLength() const