00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00027
00028
00029
00030
00031
00032
00033
00034
00035
00037
00038
00039
00040
00041
00042
00043
00044
00046
00047
00048 #include "G4SPSEneDistribution.hh"
00049
00050 #include "G4SystemOfUnits.hh"
00051 #include "Randomize.hh"
00052
00053 G4SPSEneDistribution::G4SPSEneDistribution()
00054 : particle_definition(0), eneRndm(0), Splinetemp(0)
00055 {
00056
00057
00058 particle_energy = 1.0 * MeV;
00059
00060 EnergyDisType = "Mono";
00061 weight = 1.;
00062 MonoEnergy = 1 * MeV;
00063 Emin = 0.;
00064 Emax = 1.e30;
00065 alpha = 0.;
00066 biasalpha = 0.;
00067 prob_norm = 1.0;
00068 Ezero = 0.;
00069 SE = 0.;
00070 Temp = 0.;
00071 grad = 0.;
00072 cept = 0.;
00073 Biased = false;
00074 EnergySpec = true;
00075 DiffSpec = true;
00076 IntType = "NULL";
00077 IPDFEnergyExist = false;
00078 IPDFArbExist = false;
00079
00080 ArbEmin = 0.;
00081 ArbEmax = 1.e30;
00082
00083 verbosityLevel = 0;
00084
00085 }
00086
00087 G4SPSEneDistribution::~G4SPSEneDistribution() {
00088 }
00089
00090 void G4SPSEneDistribution::SetEnergyDisType(G4String DisType) {
00091 EnergyDisType = DisType;
00092 if (EnergyDisType == "User") {
00093 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
00094 IPDFEnergyExist = false;
00095 } else if (EnergyDisType == "Arb") {
00096 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
00097 IPDFArbExist = false;
00098 } else if (EnergyDisType == "Epn") {
00099 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
00100 IPDFEnergyExist = false;
00101 EpnEnergyH = ZeroPhysVector;
00102 }
00103 }
00104
00105 void G4SPSEneDistribution::SetEmin(G4double emi) {
00106 Emin = emi;
00107 }
00108
00109 void G4SPSEneDistribution::SetEmax(G4double ema) {
00110 Emax = ema;
00111 }
00112
00113 void G4SPSEneDistribution::SetMonoEnergy(G4double menergy) {
00114 MonoEnergy = menergy;
00115 }
00116
00117 void G4SPSEneDistribution::SetBeamSigmaInE(G4double e) {
00118 SE = e;
00119 }
00120 void G4SPSEneDistribution::SetAlpha(G4double alp) {
00121 alpha = alp;
00122 }
00123
00124 void G4SPSEneDistribution::SetBiasAlpha(G4double alp) {
00125 biasalpha = alp;
00126 Biased = true;
00127 }
00128
00129 void G4SPSEneDistribution::SetTemp(G4double tem) {
00130 Temp = tem;
00131 }
00132
00133 void G4SPSEneDistribution::SetEzero(G4double eze) {
00134 Ezero = eze;
00135 }
00136
00137 void G4SPSEneDistribution::SetGradient(G4double gr) {
00138 grad = gr;
00139 }
00140
00141 void G4SPSEneDistribution::SetInterCept(G4double c) {
00142 cept = c;
00143 }
00144
00145 void G4SPSEneDistribution::UserEnergyHisto(G4ThreeVector input) {
00146 G4double ehi, val;
00147 ehi = input.x();
00148 val = input.y();
00149 if (verbosityLevel > 1) {
00150 G4cout << "In UserEnergyHisto" << G4endl;
00151 G4cout << " " << ehi << " " << val << G4endl;
00152 }
00153 UDefEnergyH.InsertValues(ehi, val);
00154 Emax = ehi;
00155 }
00156
00157 void G4SPSEneDistribution::ArbEnergyHisto(G4ThreeVector input) {
00158 G4double ehi, val;
00159 ehi = input.x();
00160 val = input.y();
00161 if (verbosityLevel > 1) {
00162 G4cout << "In ArbEnergyHisto" << G4endl;
00163 G4cout << " " << ehi << " " << val << G4endl;
00164 }
00165 ArbEnergyH.InsertValues(ehi, val);
00166 }
00167
00168 void G4SPSEneDistribution::ArbEnergyHistoFile(G4String filename) {
00169 std::ifstream infile(filename, std::ios::in);
00170 if (!infile)
00171 G4Exception("G4SPSEneDistribution::ArbEnergyHistoFile",
00172 "Event0301",FatalException,
00173 "Unable to open the histo ASCII file");
00174 G4double ehi, val;
00175 while (infile >> ehi >> val) {
00176 ArbEnergyH.InsertValues(ehi, val);
00177 }
00178 }
00179
00180 void G4SPSEneDistribution::EpnEnergyHisto(G4ThreeVector input) {
00181 G4double ehi, val;
00182 ehi = input.x();
00183 val = input.y();
00184 if (verbosityLevel > 1) {
00185 G4cout << "In EpnEnergyHisto" << G4endl;
00186 G4cout << " " << ehi << " " << val << G4endl;
00187 }
00188 EpnEnergyH.InsertValues(ehi, val);
00189 Emax = ehi;
00190 Epnflag = true;
00191 }
00192
00193 void G4SPSEneDistribution::Calculate() {
00194 if (EnergyDisType == "Cdg")
00195 CalculateCdgSpectrum();
00196 else if (EnergyDisType == "Bbody")
00197 CalculateBbodySpectrum();
00198 }
00199
00200 void G4SPSEneDistribution::CalculateCdgSpectrum() {
00201
00202
00203 G4double pfact[2] = { 8.5, 112 };
00204 G4double spind[2] = { 1.4, 2.3 };
00205 G4double ene_line[3] = { 1. * keV, 18. * keV, 1E6 * keV };
00206 G4int n_par;
00207
00208 ene_line[0] = Emin;
00209 if (Emin < 18 * keV) {
00210 n_par = 2;
00211 ene_line[2] = Emax;
00212 if (Emax < 18 * keV) {
00213 n_par = 1;
00214 ene_line[1] = Emax;
00215 }
00216 } else {
00217 n_par = 1;
00218 pfact[0] = 112.;
00219 spind[0] = 2.3;
00220 ene_line[1] = Emax;
00221 }
00222
00223
00224 CDGhist[0] = 0.;
00225 G4double omalpha;
00226 G4int i = 0;
00227
00228 while (i < n_par) {
00229 omalpha = 1. - spind[i];
00230 CDGhist[i + 1] = CDGhist[i] + (pfact[i] / omalpha) * (std::pow(
00231 ene_line[i + 1] / keV, omalpha) - std::pow(ene_line[i] / keV,
00232 omalpha));
00233 i++;
00234 }
00235
00236
00237 i = 0;
00238 while (i < n_par) {
00239 CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[n_par];
00240
00241 i++;
00242 }
00243 }
00244
00245 void G4SPSEneDistribution::CalculateBbodySpectrum() {
00246
00247
00248
00249
00250
00251
00252 G4double erange = Emax - Emin;
00253 G4double steps = erange / 10000.;
00254 G4double Bbody_y[10000];
00255 G4double k = 8.6181e-11;
00256 G4double h = 4.1362e-21;
00257 G4double c = 3e8;
00258 G4double h2 = h * h;
00259 G4double c2 = c * c;
00260 G4int count = 0;
00261 G4double sum = 0.;
00262 BBHist[0] = 0.;
00263 while (count < 10000) {
00264 Bbody_x[count] = Emin + G4double(count * steps);
00265 Bbody_y[count] = (2. * std::pow(Bbody_x[count], 2.)) / (h2 * c2
00266 * (std::exp(Bbody_x[count] / (k * Temp)) - 1.));
00267 sum = sum + Bbody_y[count];
00268 BBHist[count + 1] = BBHist[count] + Bbody_y[count];
00269 count++;
00270 }
00271
00272 Bbody_x[10000] = Emax;
00273
00274 count = 0;
00275 while (count < 10001) {
00276 BBHist[count] = BBHist[count] / sum;
00277 count++;
00278 }
00279 }
00280
00281 void G4SPSEneDistribution::InputEnergySpectra(G4bool value) {
00282
00283 EnergySpec = value;
00284 if (verbosityLevel > 1)
00285 G4cout << "EnergySpec has value " << EnergySpec << G4endl;
00286 }
00287
00288 void G4SPSEneDistribution::InputDifferentialSpectra(G4bool value) {
00289
00290 DiffSpec = value;
00291 if (verbosityLevel > 1)
00292 G4cout << "Diffspec has value " << DiffSpec << G4endl;
00293 }
00294
00295 void G4SPSEneDistribution::ArbInterpolate(G4String IType) {
00296 if (EnergyDisType != "Arb")
00297 G4cout << "Error: this is for arbitrary distributions" << G4endl;
00298 IntType = IType;
00299 ArbEmax = ArbEnergyH.GetMaxLowEdgeEnergy();
00300 ArbEmin = ArbEnergyH.GetMinLowEdgeEnergy();
00301
00302
00303 if (IntType == "Lin")
00304 LinearInterpolation();
00305 if (IntType == "Log")
00306 LogInterpolation();
00307 if (IntType == "Exp")
00308 ExpInterpolation();
00309 if (IntType == "Spline")
00310 SplineInterpolation();
00311 }
00312
00313 void G4SPSEneDistribution::LinearInterpolation() {
00314
00315
00316
00317
00318
00319 G4double Area_seg[1024];
00320 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
00321 G4int i, count;
00322 G4int maxi = ArbEnergyH.GetVectorLength();
00323 for (i = 0; i < maxi; i++) {
00324 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
00325 Arb_y[i] = ArbEnergyH(size_t(i));
00326 }
00327
00328
00329 if (DiffSpec == false) {
00330
00331 for (count = 0; count < maxi - 1; count++) {
00332 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
00333 / (Arb_x[count + 1] - Arb_x[count]);
00334 }
00335 maxi--;
00336 }
00337
00338 if (EnergySpec == false) {
00339
00340
00341 if (particle_definition == NULL)
00342 G4cout << "Error: particle not defined" << G4endl;
00343 else {
00344
00345
00346
00347 G4double mass = particle_definition->GetPDGMass();
00348
00349 G4double total_energy;
00350 for (count = 0; count < maxi; count++) {
00351 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
00352 * mass));
00353
00354 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
00355 Arb_x[count] = total_energy - mass;
00356 }
00357 }
00358 }
00359
00360 i = 1;
00361 Arb_grad[0] = 0.;
00362 Arb_cept[0] = 0.;
00363 Area_seg[0] = 0.;
00364 Arb_Cum_Area[0] = 0.;
00365 while (i < maxi) {
00366
00367 Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[i - 1]);
00368 if (verbosityLevel == 2)
00369 G4cout << Arb_grad[i] << G4endl;
00370 if (Arb_grad[i] > 0.) {
00371 if (verbosityLevel == 2)
00372 G4cout << "Arb_grad is positive" << G4endl;
00373 Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]);
00374 } else if (Arb_grad[i] < 0.) {
00375 if (verbosityLevel == 2)
00376 G4cout << "Arb_grad is negative" << G4endl;
00377 Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]);
00378 } else {
00379 if (verbosityLevel == 2)
00380 G4cout << "Arb_grad is 0." << G4endl;
00381 Arb_cept[i] = Arb_y[i];
00382 }
00383
00384 Area_seg[i] = ((Arb_grad[i] / 2) * (Arb_x[i] * Arb_x[i] - Arb_x[i - 1]
00385 * Arb_x[i - 1]) + Arb_cept[i] * (Arb_x[i] - Arb_x[i - 1]));
00386 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
00387 sum = sum + Area_seg[i];
00388 if (verbosityLevel == 2)
00389 G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i]
00390 << G4endl;
00391 i++;
00392 }
00393
00394 i = 0;
00395 while (i < maxi) {
00396 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
00397 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
00398 i++;
00399 }
00400
00401
00402 ArbEnergyH.ScaleVector(1., 1./sum);
00403
00404 if (verbosityLevel >= 1) {
00405 G4cout << "Leaving LinearInterpolation" << G4endl;
00406 ArbEnergyH.DumpValues();
00407 IPDFArbEnergyH.DumpValues();
00408 }
00409 }
00410
00411 void G4SPSEneDistribution::LogInterpolation() {
00412
00413
00414
00415
00416
00417 G4double Area_seg[1024];
00418 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
00419 G4int i, count;
00420 G4int maxi = ArbEnergyH.GetVectorLength();
00421 for (i = 0; i < maxi; i++) {
00422 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
00423 Arb_y[i] = ArbEnergyH(size_t(i));
00424 }
00425
00426
00427 if (DiffSpec == false) {
00428
00429 for (count = 0; count < maxi - 1; count++) {
00430 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
00431 / (Arb_x[count + 1] - Arb_x[count]);
00432 }
00433 maxi--;
00434 }
00435
00436 if (EnergySpec == false) {
00437
00438
00439 if (particle_definition == NULL)
00440 G4cout << "Error: particle not defined" << G4endl;
00441 else {
00442
00443
00444
00445 G4double mass = particle_definition->GetPDGMass();
00446
00447 G4double total_energy;
00448 for (count = 0; count < maxi; count++) {
00449 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
00450 * mass));
00451
00452 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
00453 Arb_x[count] = total_energy - mass;
00454 }
00455 }
00456 }
00457
00458 i = 1;
00459 Arb_alpha[0] = 0.;
00460 Arb_Const[0] = 0.;
00461 Area_seg[0] = 0.;
00462 Arb_Cum_Area[0] = 0.;
00463 if (Arb_x[0] <= 0. || Arb_y[0] <= 0.) {
00464 G4cout << "You should not use log interpolation with points <= 0."
00465 << G4endl;
00466 G4cout << "These will be changed to 1e-20, which may cause problems"
00467 << G4endl;
00468 if (Arb_x[0] <= 0.)
00469 Arb_x[0] = 1e-20;
00470 if (Arb_y[0] <= 0.)
00471 Arb_y[0] = 1e-20;
00472 }
00473
00474 G4double alp;
00475 while (i < maxi) {
00476
00477 if (Arb_x[i] <= 0. || Arb_y[i] <= 0.) {
00478 G4cout << "You should not use log interpolation with points <= 0."
00479 << G4endl;
00480 G4cout
00481 << "These will be changed to 1e-20, which may cause problems"
00482 << G4endl;
00483 if (Arb_x[i] <= 0.)
00484 Arb_x[i] = 1e-20;
00485 if (Arb_y[i] <= 0.)
00486 Arb_y[i] = 1e-20;
00487 }
00488
00489 Arb_alpha[i] = (std::log10(Arb_y[i]) - std::log10(Arb_y[i - 1]))
00490 / (std::log10(Arb_x[i]) - std::log10(Arb_x[i - 1]));
00491 Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[i], Arb_alpha[i]));
00492 alp = Arb_alpha[i] + 1;
00493 if (alp == 0.) {
00494 Area_seg[i] = Arb_Const[i] * (std::log(Arb_x[i]) - std::log(Arb_x[i - 1]));
00495 } else {
00496 Area_seg[i] = (Arb_Const[i] / alp) * (std::pow(Arb_x[i], alp)
00497 - std::pow(Arb_x[i - 1], alp));
00498 }
00499 sum = sum + Area_seg[i];
00500 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
00501 if (verbosityLevel == 2)
00502 G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl;
00503 i++;
00504 }
00505
00506 i = 0;
00507 while (i < maxi) {
00508 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
00509 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
00510 i++;
00511 }
00512
00513
00514 ArbEnergyH.ScaleVector(1., 1./sum);
00515
00516 if (verbosityLevel >= 1)
00517 G4cout << "Leaving LogInterpolation " << G4endl;
00518 }
00519
00520 void G4SPSEneDistribution::ExpInterpolation() {
00521
00522
00523
00524
00525
00526 G4double Area_seg[1024];
00527 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
00528 G4int i, count;
00529 G4int maxi = ArbEnergyH.GetVectorLength();
00530 for (i = 0; i < maxi; i++) {
00531 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
00532 Arb_y[i] = ArbEnergyH(size_t(i));
00533 }
00534
00535
00536 if (DiffSpec == false) {
00537
00538 for (count = 0; count < maxi - 1; count++) {
00539 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
00540 / (Arb_x[count + 1] - Arb_x[count]);
00541 }
00542 maxi--;
00543 }
00544
00545 if (EnergySpec == false) {
00546
00547
00548 if (particle_definition == NULL)
00549 G4cout << "Error: particle not defined" << G4endl;
00550 else {
00551
00552
00553
00554 G4double mass = particle_definition->GetPDGMass();
00555
00556 G4double total_energy;
00557 for (count = 0; count < maxi; count++) {
00558 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
00559 * mass));
00560
00561 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
00562 Arb_x[count] = total_energy - mass;
00563 }
00564 }
00565 }
00566
00567 i = 1;
00568 Arb_ezero[0] = 0.;
00569 Arb_Const[0] = 0.;
00570 Area_seg[0] = 0.;
00571 Arb_Cum_Area[0] = 0.;
00572 while (i < maxi) {
00573 G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i - 1]);
00574 if (test > 0. || test < 0.) {
00575 Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1]) / (std::log(Arb_y[i])
00576 - std::log(Arb_y[i - 1]));
00577 Arb_Const[i] = Arb_y[i] / (std::exp(-Arb_x[i] / Arb_ezero[i]));
00578 Area_seg[i] = -(Arb_Const[i] * Arb_ezero[i]) * (std::exp(-Arb_x[i]
00579 / Arb_ezero[i]) - std::exp(-Arb_x[i - 1] / Arb_ezero[i]));
00580 } else {
00581 G4cout << "Flat line segment: problem" << G4endl;
00582 Arb_ezero[i] = 0.;
00583 Arb_Const[i] = 0.;
00584 Area_seg[i] = 0.;
00585 }
00586 sum = sum + Area_seg[i];
00587 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
00588 if (verbosityLevel == 2)
00589 G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl;
00590 i++;
00591 }
00592
00593 i = 0;
00594 while (i < maxi) {
00595 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
00596 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
00597 i++;
00598 }
00599
00600
00601 ArbEnergyH.ScaleVector(1., 1./sum);
00602
00603 if (verbosityLevel >= 1)
00604 G4cout << "Leaving ExpInterpolation " << G4endl;
00605 }
00606
00607 void G4SPSEneDistribution::SplineInterpolation() {
00608
00609
00610
00611
00612
00613
00614
00615 G4double sum, Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
00616 G4int i, count;
00617
00618 G4int maxi = ArbEnergyH.GetVectorLength();
00619 for (i = 0; i < maxi; i++) {
00620 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
00621 Arb_y[i] = ArbEnergyH(size_t(i));
00622 }
00623
00624
00625 if (DiffSpec == false) {
00626
00627 for (count = 0; count < maxi - 1; count++) {
00628 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
00629 / (Arb_x[count + 1] - Arb_x[count]);
00630 }
00631 maxi--;
00632 }
00633
00634 if (EnergySpec == false) {
00635
00636
00637 if (particle_definition == NULL)
00638 G4Exception("G4SPSEneDistribution::SplineInterpolation",
00639 "Event0302",FatalException,
00640 "Error: particle not defined");
00641 else {
00642
00643
00644
00645 G4double mass = particle_definition->GetPDGMass();
00646
00647 G4double total_energy;
00648 for (count = 0; count < maxi; count++) {
00649 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
00650 * mass));
00651
00652 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
00653 Arb_x[count] = total_energy - mass;
00654 }
00655 }
00656 }
00657
00658
00659 i = 1;
00660 Arb_Cum_Area[0] = 0.;
00661 sum = 0.;
00662 Splinetemp = new G4DataInterpolation(Arb_x, Arb_y, maxi, 0., 0.);
00663 G4double ei[101],prob[101];
00664 while (i < maxi) {
00665
00666 G4double de = (Arb_x[i] - Arb_x[i - 1])/100.;
00667 G4double area = 0.;
00668
00669 for (count = 0; count < 101; count++) {
00670 ei[count] = Arb_x[i - 1] + de*count ;
00671 prob[count] = Splinetemp->CubicSplineInterpolation(ei[count]);
00672 if (prob[count] < 0.) {
00673 G4ExceptionDescription ED;
00674 ED << "Warning: G4DataInterpolation returns value < 0 " << prob[count] <<" "<<ei[count]<< G4endl;
00675 G4Exception("G4SPSEneDistribution::SplineInterpolation","Event0303",
00676 FatalException,ED);
00677 }
00678 area += prob[count]*de;
00679 }
00680 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + area;
00681 sum += area;
00682
00683 prob[0] = prob[0]/(area/de);
00684 for (count = 1; count < 100; count++)
00685 prob[count] = prob[count-1] + prob[count]/(area/de);
00686
00687 SplineInt[i] = new G4DataInterpolation(prob, ei, 101, 0., 0.);
00688
00689 i++;
00690 }
00691 i = 0;
00692 while (i < maxi) {
00693 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
00694 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
00695 i++;
00696 }
00697
00698 ArbEnergyH.ScaleVector(1., 1./sum);
00699
00700 if (verbosityLevel > 0)
00701 G4cout << "Leaving SplineInterpolation " << G4endl;
00702 }
00703
00704 void G4SPSEneDistribution::GenerateMonoEnergetic() {
00705
00706 particle_energy = MonoEnergy;
00707 }
00708
00709 void G4SPSEneDistribution::GenerateGaussEnergies() {
00710
00711 particle_energy = G4RandGauss::shoot(MonoEnergy,SE);
00712 if (particle_energy < 0) particle_energy = 0.;
00713 }
00714
00715 void G4SPSEneDistribution::GenerateLinearEnergies(G4bool bArb = false) {
00716 G4double rndm;
00717 G4double emaxsq = std::pow(Emax, 2.);
00718 G4double eminsq = std::pow(Emin, 2.);
00719 G4double intersq = std::pow(cept, 2.);
00720
00721 if (bArb)
00722 rndm = G4UniformRand();
00723 else
00724 rndm = eneRndm->GenRandEnergy();
00725
00726 G4double bracket = ((grad / 2.) * (emaxsq - eminsq) + cept * (Emax - Emin));
00727 bracket = bracket * rndm;
00728 bracket = bracket + (grad / 2.) * eminsq + cept * Emin;
00729
00730 bracket = -bracket;
00731
00732 if (grad != 0.) {
00733 G4double sqbrack = (intersq - 4 * (grad / 2.) * (bracket));
00734
00735 sqbrack = std::sqrt(sqbrack);
00736 G4double root1 = -cept + sqbrack;
00737 root1 = root1 / (2. * (grad / 2.));
00738
00739 G4double root2 = -cept - sqbrack;
00740 root2 = root2 / (2. * (grad / 2.));
00741
00742
00743
00744 if (root1 > Emin && root1 < Emax)
00745 particle_energy = root1;
00746 if (root2 > Emin && root2 < Emax)
00747 particle_energy = root2;
00748 } else if (grad == 0.)
00749
00750 particle_energy = bracket / cept;
00751
00752 if (particle_energy < 0.)
00753 particle_energy = -particle_energy;
00754
00755 if (verbosityLevel >= 1)
00756 G4cout << "Energy is " << particle_energy << G4endl;
00757 }
00758
00759 void G4SPSEneDistribution::GeneratePowEnergies(G4bool bArb = false) {
00760
00761
00762
00763 G4double rndm;
00764 G4double emina, emaxa;
00765
00766 emina = std::pow(Emin, alpha + 1);
00767 emaxa = std::pow(Emax, alpha + 1);
00768
00769 if (bArb)
00770 rndm = G4UniformRand();
00771 else
00772 rndm = eneRndm->GenRandEnergy();
00773
00774 if (alpha != -1.) {
00775 particle_energy = ((rndm * (emaxa - emina)) + emina);
00776 particle_energy = std::pow(particle_energy, (1. / (alpha + 1.)));
00777 } else {
00778 particle_energy = (std::log(Emin) + rndm * (std::log(Emax) - std::log(
00779 Emin)));
00780 particle_energy = std::exp(particle_energy);
00781 }
00782 if (verbosityLevel >= 1)
00783 G4cout << "Energy is " << particle_energy << G4endl;
00784 }
00785
00786 void G4SPSEneDistribution::GenerateBiasPowEnergies() {
00787
00788
00789
00790 G4double rndm;
00791 G4double emina, emaxa, emin, emax;
00792
00793 G4double normal = 1. ;
00794
00795 emin = Emin;
00796 emax = Emax;
00797
00798
00799
00800
00801
00802 rndm = eneRndm->GenRandEnergy();
00803
00804 if (biasalpha != -1.) {
00805 emina = std::pow(emin, biasalpha + 1);
00806 emaxa = std::pow(emax, biasalpha + 1);
00807 particle_energy = ((rndm * (emaxa - emina)) + emina);
00808 particle_energy = std::pow(particle_energy, (1. / (biasalpha + 1.)));
00809 normal = 1./(1+biasalpha) * (emaxa - emina);
00810 } else {
00811 particle_energy = (std::log(emin) + rndm * (std::log(emax) - std::log(
00812 emin)));
00813 particle_energy = std::exp(particle_energy);
00814 normal = std::log(emax) - std::log(emin) ;
00815 }
00816 weight = GetProbability(particle_energy) / (std::pow(particle_energy,biasalpha)/normal);
00817
00818 if (verbosityLevel >= 1)
00819 G4cout << "Energy is " << particle_energy << G4endl;
00820 }
00821
00822 void G4SPSEneDistribution::GenerateExpEnergies(G4bool bArb = false) {
00823
00824
00825 G4double rndm;
00826
00827 if (bArb)
00828 rndm = G4UniformRand();
00829 else
00830 rndm = eneRndm->GenRandEnergy();
00831
00832 particle_energy = -Ezero * (std::log(rndm * (std::exp(-Emax / Ezero)
00833 - std::exp(-Emin / Ezero)) + std::exp(-Emin / Ezero)));
00834 if (verbosityLevel >= 1)
00835 G4cout << "Energy is " << particle_energy << G4endl;
00836 }
00837
00838 void G4SPSEneDistribution::GenerateBremEnergies() {
00839
00840
00841
00842
00843 G4double rndm;
00844 rndm = eneRndm->GenRandEnergy();
00845 G4double expmax, expmin, k;
00846
00847 k = 8.6181e-11;
00848 G4double ksq = std::pow(k, 2.);
00849 G4double Tsq = std::pow(Temp, 2.);
00850
00851 expmax = std::exp(-Emax / (k * Temp));
00852 expmin = std::exp(-Emin / (k * Temp));
00853
00854
00855
00856
00857 if (expmax == 0.)
00858 G4cout << "*****EXPMAX=0. Choose different E's or Temp" << G4endl;
00859 if (expmin == 0.)
00860 G4cout << "*****EXPMIN=0. Choose different E's or Temp" << G4endl;
00861
00862 G4double tempvar = rndm * ((-k) * Temp * (Emax * expmax - Emin * expmin)
00863 - (ksq * Tsq * (expmax - expmin)));
00864
00865 G4double bigc = (tempvar - k * Temp * Emin * expmin - ksq * Tsq * expmin)
00866 / (-k * Temp);
00867
00868
00869
00870
00871
00872 G4double erange = Emax - Emin;
00873 G4double steps = erange / 1000.;
00874 G4int i;
00875 G4double etest, diff, err;
00876
00877 err = 100000.;
00878
00879 for (i = 1; i < 1000; i++) {
00880 etest = Emin + (i - 1) * steps;
00881
00882 diff = etest * (std::exp(-etest / (k * Temp))) + k * Temp * (std::exp(
00883 -etest / (k * Temp))) - bigc;
00884
00885 if (diff < 0.)
00886 diff = -diff;
00887
00888 if (diff < err) {
00889 err = diff;
00890 particle_energy = etest;
00891 }
00892 }
00893 if (verbosityLevel >= 1)
00894 G4cout << "Energy is " << particle_energy << G4endl;
00895 }
00896
00897 void G4SPSEneDistribution::GenerateBbodyEnergies() {
00898
00899
00900
00901
00902 G4double rndm;
00903 G4int nabove, nbelow = 0, middle;
00904 nabove = 10001;
00905 rndm = eneRndm->GenRandEnergy();
00906
00907
00908 while (nabove - nbelow > 1) {
00909 middle = (nabove + nbelow) / 2;
00910 if (rndm == BBHist[middle])
00911 break;
00912 if (rndm < BBHist[middle])
00913 nabove = middle;
00914 else
00915 nbelow = middle;
00916 }
00917
00918
00919 G4double x1, x2, y1, y2, t, q;
00920 x1 = Bbody_x[nbelow];
00921 x2 = Bbody_x[nbelow + 1];
00922 y1 = BBHist[nbelow];
00923 y2 = BBHist[nbelow + 1];
00924 t = (y2 - y1) / (x2 - x1);
00925 q = y1 - t * x1;
00926
00927 particle_energy = (rndm - q) / t;
00928
00929 if (verbosityLevel >= 1) {
00930 G4cout << "Energy is " << particle_energy << G4endl;
00931 }
00932 }
00933
00934 void G4SPSEneDistribution::GenerateCdgEnergies() {
00935
00936
00937
00938
00939
00940
00941 G4double rndm, rndm2;
00942 G4double ene_line[3];
00943 G4double omalpha[2];
00944 if (Emin < 18 * keV && Emax < 18 * keV) {
00945 omalpha[0] = 1. - 1.4;
00946 ene_line[0] = Emin;
00947 ene_line[1] = Emax;
00948 }
00949 if (Emin < 18 * keV && Emax > 18 * keV) {
00950 omalpha[0] = 1. - 1.4;
00951 omalpha[1] = 1. - 2.3;
00952 ene_line[0] = Emin;
00953 ene_line[1] = 18. * keV;
00954 ene_line[2] = Emax;
00955 }
00956 if (Emin > 18 * keV) {
00957 omalpha[0] = 1. - 2.3;
00958 ene_line[0] = Emin;
00959 ene_line[1] = Emax;
00960 }
00961 rndm = eneRndm->GenRandEnergy();
00962 rndm2 = eneRndm->GenRandEnergy();
00963
00964 G4int i = 0;
00965 while (rndm >= CDGhist[i]) {
00966 i++;
00967 }
00968
00969 particle_energy = (std::pow(ene_line[i - 1], omalpha[i - 1]) + (std::pow(
00970 ene_line[i], omalpha[i - 1]) - std::pow(ene_line[i - 1], omalpha[i
00971 - 1])) * rndm2);
00972 particle_energy = std::pow(particle_energy, (1. / omalpha[i - 1]));
00973
00974 if (verbosityLevel >= 1)
00975 G4cout << "Energy is " << particle_energy << G4endl;
00976 }
00977
00978 void G4SPSEneDistribution::GenUserHistEnergies() {
00979
00980
00981 if (IPDFEnergyExist == false) {
00982 G4int ii;
00983 G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
00984 G4double bins[1024], vals[1024], sum;
00985 sum = 0.;
00986
00987 if ((EnergySpec == false) && (particle_definition == NULL))
00988 G4cout << "Error: particle definition is NULL" << G4endl;
00989
00990 if (maxbin > 1024) {
00991 G4cout << "Maxbin > 1024" << G4endl;
00992 G4cout << "Setting maxbin to 1024, other bins are lost" << G4endl;
00993 }
00994
00995 if (DiffSpec == false)
00996 G4cout << "Histograms are Differential!!! " << G4endl;
00997 else {
00998 bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
00999 vals[0] = UDefEnergyH(size_t(0));
01000 sum = vals[0];
01001 for (ii = 1; ii < maxbin; ii++) {
01002 bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
01003 vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1];
01004 sum = sum + UDefEnergyH(size_t(ii));
01005 }
01006 }
01007
01008 if (EnergySpec == false) {
01009 G4double mass = particle_definition->GetPDGMass();
01010
01011
01012
01013 for (ii = 1; ii < maxbin; ii++) {
01014 vals[ii] = vals[ii] * (bins[ii] - bins[ii - 1]);
01015 }
01016
01017
01018 for (ii = 0; ii < maxbin; ii++) {
01019 bins[ii] = std::sqrt((bins[ii] * bins[ii]) + (mass * mass))
01020 - mass;
01021 }
01022 for (ii = 1; ii < maxbin; ii++) {
01023 vals[ii] = vals[ii] / (bins[ii] - bins[ii - 1]);
01024 }
01025 sum = vals[maxbin - 1];
01026 vals[0] = 0.;
01027 }
01028 for (ii = 0; ii < maxbin; ii++) {
01029 vals[ii] = vals[ii] / sum;
01030 IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
01031 }
01032
01033
01034 IPDFEnergyExist = true;
01035 if (verbosityLevel > 1)
01036 IPDFEnergyH.DumpValues();
01037 }
01038
01039
01040 G4double rndm = eneRndm->GenRandEnergy();
01041 particle_energy = IPDFEnergyH.GetEnergy(rndm);
01042
01043 if (verbosityLevel >= 1)
01044 G4cout << "Energy is " << particle_energy << G4endl;
01045 }
01046
01047 void G4SPSEneDistribution::GenArbPointEnergies() {
01048 if (verbosityLevel > 0)
01049 G4cout << "In GenArbPointEnergies" << G4endl;
01050 G4double rndm;
01051 rndm = eneRndm->GenRandEnergy();
01052
01053
01054
01055 G4int nabove, nbelow = 0, middle;
01056 nabove = IPDFArbEnergyH.GetVectorLength();
01057
01058
01059 while (nabove - nbelow > 1) {
01060 middle = (nabove + nbelow) / 2;
01061 if (rndm == IPDFArbEnergyH(size_t(middle)))
01062 break;
01063 if (rndm < IPDFArbEnergyH(size_t(middle)))
01064 nabove = middle;
01065 else
01066 nbelow = middle;
01067 }
01068 if (IntType == "Lin") {
01069 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
01070 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
01071 grad = Arb_grad[nbelow + 1];
01072 cept = Arb_cept[nbelow + 1];
01073
01074 GenerateLinearEnergies(true);
01075 } else if (IntType == "Log") {
01076 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
01077 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
01078 alpha = Arb_alpha[nbelow + 1];
01079
01080 GeneratePowEnergies(true);
01081 } else if (IntType == "Exp") {
01082 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
01083 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
01084 Ezero = Arb_ezero[nbelow + 1];
01085
01086 GenerateExpEnergies(true);
01087 } else if (IntType == "Spline") {
01088 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
01089 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
01090 particle_energy = -1e100;
01091 rndm = eneRndm->GenRandEnergy();
01092 while (particle_energy < Emin || particle_energy > Emax) {
01093 particle_energy = SplineInt[nbelow+1]->CubicSplineInterpolation(rndm);
01094 rndm = eneRndm->GenRandEnergy();
01095 }
01096 if (verbosityLevel >= 1)
01097 G4cout << "Energy is " << particle_energy << G4endl;
01098 } else
01099 G4cout << "Error: IntType unknown type" << G4endl;
01100 }
01101
01102 void G4SPSEneDistribution::GenEpnHistEnergies() {
01103
01104
01105
01106 if (Epnflag == true)
01107
01108 {
01109
01110 ConvertEPNToEnergy();
01111
01112
01113 }
01114
01115
01116 if (IPDFEnergyExist == false) {
01117
01118 G4double bins[1024], vals[1024], sum;
01119 G4int ii;
01120 G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
01121 bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
01122 vals[0] = UDefEnergyH(size_t(0));
01123 sum = vals[0];
01124 for (ii = 1; ii < maxbin; ii++) {
01125 bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
01126 vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1];
01127 sum = sum + UDefEnergyH(size_t(ii));
01128 }
01129
01130 for (ii = 0; ii < maxbin; ii++) {
01131 vals[ii] = vals[ii] / sum;
01132 IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
01133 }
01134
01135 IPDFEnergyExist = true;
01136 }
01137
01138
01139 G4double rndm = eneRndm->GenRandEnergy();
01140 particle_energy = IPDFEnergyH.GetEnergy(rndm);
01141
01142 if (verbosityLevel >= 1)
01143 G4cout << "Energy is " << particle_energy << G4endl;
01144 }
01145
01146 void G4SPSEneDistribution::ConvertEPNToEnergy() {
01147
01148
01149
01150
01151 if (particle_definition == NULL)
01152 G4cout << "Error: particle not defined" << G4endl;
01153 else {
01154
01155
01156 G4int Bary = particle_definition->GetBaryonNumber();
01157
01158
01159 G4int count, maxcount;
01160 maxcount = G4int(EpnEnergyH.GetVectorLength());
01161
01162 G4double ebins[1024], evals[1024];
01163 if (maxcount > 1024) {
01164 G4cout << "Histogram contains more than 1024 bins!" << G4endl;
01165 G4cout << "Those above 1024 will be ignored" << G4endl;
01166 maxcount = 1024;
01167 }
01168 if (maxcount < 1) {
01169 G4cout << "Histogram contains less than 1 bin!" << G4endl;
01170 G4cout << "Redefine the histogram" << G4endl;
01171 return;
01172 }
01173 for (count = 0; count < maxcount; count++) {
01174
01175 ebins[count] = EpnEnergyH.GetLowEdgeEnergy(size_t(count));
01176 evals[count] = EpnEnergyH(size_t(count));
01177 }
01178
01179
01180 for (count = 0; count < maxcount; count++) {
01181 ebins[count] = ebins[count] * Bary;
01182 }
01183
01184
01185 Emin = ebins[0];
01186 if (maxcount > 1)
01187 Emax = ebins[maxcount - 1];
01188 else
01189 Emax = ebins[0];
01190
01191 for (count = 0; count < maxcount; count++) {
01192 UDefEnergyH.InsertValues(ebins[count], evals[count]);
01193 }
01194 Epnflag = false;
01195 }
01196 }
01197
01198
01199 void G4SPSEneDistribution::ReSetHist(G4String atype) {
01200 if (atype == "energy") {
01201 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
01202 IPDFEnergyExist = false;
01203 Emin = 0.;
01204 Emax = 1e30;
01205 } else if (atype == "arb") {
01206 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
01207 IPDFArbExist = false;
01208 } else if (atype == "epn") {
01209 UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
01210 IPDFEnergyExist = false;
01211 EpnEnergyH = ZeroPhysVector;
01212 } else {
01213 G4cout << "Error, histtype not accepted " << G4endl;
01214 }
01215 }
01216
01217 G4double G4SPSEneDistribution::GenerateOne(G4ParticleDefinition* a) {
01218 particle_definition = a;
01219 particle_energy = -1.;
01220
01221 while ((EnergyDisType == "Arb") ? (particle_energy < ArbEmin
01222 || particle_energy > ArbEmax) : (particle_energy < Emin
01223 || particle_energy > Emax)) {
01224 if (Biased) {
01225 GenerateBiasPowEnergies();
01226 } else {
01227 if (EnergyDisType == "Mono")
01228 GenerateMonoEnergetic();
01229 else if (EnergyDisType == "Lin")
01230 GenerateLinearEnergies();
01231 else if (EnergyDisType == "Pow")
01232 GeneratePowEnergies();
01233 else if (EnergyDisType == "Exp")
01234 GenerateExpEnergies();
01235 else if (EnergyDisType == "Gauss")
01236 GenerateGaussEnergies();
01237 else if (EnergyDisType == "Brem")
01238 GenerateBremEnergies();
01239 else if (EnergyDisType == "Bbody")
01240 GenerateBbodyEnergies();
01241 else if (EnergyDisType == "Cdg")
01242 GenerateCdgEnergies();
01243 else if (EnergyDisType == "User")
01244 GenUserHistEnergies();
01245 else if (EnergyDisType == "Arb")
01246 GenArbPointEnergies();
01247 else if (EnergyDisType == "Epn")
01248 GenEpnHistEnergies();
01249 else
01250 G4cout << "Error: EnergyDisType has unusual value" << G4endl;
01251 }
01252 }
01253 return particle_energy;
01254 }
01255
01256 G4double G4SPSEneDistribution::GetProbability(G4double ene) {
01257 G4double prob = 1.;
01258
01259 if (EnergyDisType == "Lin") {
01260 if (prob_norm == 1.) {
01261 prob_norm = 0.5*grad*Emax*Emax + cept*Emax - 0.5*grad*Emin*Emin - cept*Emin;
01262 }
01263 prob = cept + grad * ene;
01264 prob /= prob_norm;
01265 }
01266 else if (EnergyDisType == "Pow") {
01267 if (prob_norm == 1.) {
01268 if (alpha != -1.) {
01269 G4double emina = std::pow(Emin, alpha + 1);
01270 G4double emaxa = std::pow(Emax, alpha + 1);
01271 prob_norm = 1./(1.+alpha) * (emaxa - emina);
01272 } else {
01273 prob_norm = std::log(Emax) - std::log(Emin) ;
01274 }
01275 }
01276 prob = std::pow(ene, alpha)/prob_norm;
01277 }
01278 else if (EnergyDisType == "Exp"){
01279 if (prob_norm == 1.) {
01280 prob_norm = -Ezero*(std::exp(-Emax/Ezero) - std::exp(Emin/Ezero));
01281 }
01282 prob = std::exp(-ene / Ezero);
01283 prob /= prob_norm;
01284 }
01285 else if (EnergyDisType == "Arb") {
01286 prob = ArbEnergyH.Value(ene);
01287
01288
01289
01290 if (prob <= 0.) {
01291
01292 G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << G4endl;
01293 prob = 1e-30;
01294 }
01295
01296 }
01297 else
01298 G4cout << "Error: EnergyDisType not supported" << G4endl;
01299
01300 return prob;
01301 }