[816] | 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 | // |
---|
| 26 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 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 | // |
---|
| 36 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 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 | // |
---|
| 45 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 46 | // |
---|
| 47 | #include "Randomize.hh" |
---|
| 48 | //#include <cmath> |
---|
| 49 | |
---|
| 50 | #include "G4SPSEneDistribution.hh" |
---|
| 51 | |
---|
[1347] | 52 | G4SPSEneDistribution::G4SPSEneDistribution() { |
---|
| 53 | // |
---|
| 54 | // Initialise all variables |
---|
| 55 | particle_energy = 1.0 * MeV; |
---|
[816] | 56 | |
---|
[1347] | 57 | EnergyDisType = "Mono"; |
---|
| 58 | weight = 1.; |
---|
| 59 | MonoEnergy = 1 * MeV; |
---|
| 60 | Emin = 0.; |
---|
| 61 | Emax = 1.e30; |
---|
| 62 | alpha = 0.; |
---|
| 63 | biasalpha = 0.; |
---|
| 64 | prob_norm = 1.0; |
---|
| 65 | Ezero = 0.; |
---|
| 66 | SE = 0.; |
---|
| 67 | Temp = 0.; |
---|
| 68 | grad = 0.; |
---|
| 69 | cept = 0.; |
---|
| 70 | Biased = false; // not biased |
---|
| 71 | EnergySpec = true; // true - energy spectra, false - momentum spectra |
---|
| 72 | DiffSpec = true; // true - differential spec, false integral spec |
---|
| 73 | IntType = "NULL"; // Interpolation type |
---|
| 74 | IPDFEnergyExist = false; |
---|
| 75 | IPDFArbExist = false; |
---|
[816] | 76 | |
---|
[1347] | 77 | ArbEmin = 0.; |
---|
| 78 | ArbEmax = 1.e30; |
---|
[816] | 79 | |
---|
[1347] | 80 | verbosityLevel = 0; |
---|
[816] | 81 | |
---|
| 82 | } |
---|
| 83 | |
---|
[1347] | 84 | G4SPSEneDistribution::~G4SPSEneDistribution() { |
---|
| 85 | } |
---|
[816] | 86 | |
---|
[1347] | 87 | void G4SPSEneDistribution::SetEnergyDisType(G4String DisType) { |
---|
| 88 | EnergyDisType = DisType; |
---|
| 89 | if (EnergyDisType == "User") { |
---|
| 90 | UDefEnergyH = IPDFEnergyH = ZeroPhysVector; |
---|
| 91 | IPDFEnergyExist = false; |
---|
| 92 | } else if (EnergyDisType == "Arb") { |
---|
| 93 | ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector; |
---|
| 94 | IPDFArbExist = false; |
---|
| 95 | } else if (EnergyDisType == "Epn") { |
---|
| 96 | UDefEnergyH = IPDFEnergyH = ZeroPhysVector; |
---|
| 97 | IPDFEnergyExist = false; |
---|
| 98 | EpnEnergyH = ZeroPhysVector; |
---|
| 99 | } |
---|
[816] | 100 | } |
---|
| 101 | |
---|
[1347] | 102 | void G4SPSEneDistribution::SetEmin(G4double emi) { |
---|
| 103 | Emin = emi; |
---|
[816] | 104 | } |
---|
| 105 | |
---|
[1347] | 106 | void G4SPSEneDistribution::SetEmax(G4double ema) { |
---|
| 107 | Emax = ema; |
---|
[816] | 108 | } |
---|
| 109 | |
---|
[1347] | 110 | void G4SPSEneDistribution::SetMonoEnergy(G4double menergy) { |
---|
| 111 | MonoEnergy = menergy; |
---|
[816] | 112 | } |
---|
| 113 | |
---|
[1347] | 114 | void G4SPSEneDistribution::SetBeamSigmaInE(G4double e) { |
---|
| 115 | SE = e; |
---|
[816] | 116 | } |
---|
[1347] | 117 | void G4SPSEneDistribution::SetAlpha(G4double alp) { |
---|
| 118 | alpha = alp; |
---|
[816] | 119 | } |
---|
| 120 | |
---|
[1347] | 121 | void G4SPSEneDistribution::SetBiasAlpha(G4double alp) { |
---|
| 122 | biasalpha = alp; |
---|
| 123 | Biased = true; |
---|
[816] | 124 | } |
---|
| 125 | |
---|
[1347] | 126 | void G4SPSEneDistribution::SetTemp(G4double tem) { |
---|
| 127 | Temp = tem; |
---|
[816] | 128 | } |
---|
| 129 | |
---|
[1347] | 130 | void G4SPSEneDistribution::SetEzero(G4double eze) { |
---|
| 131 | Ezero = eze; |
---|
[816] | 132 | } |
---|
| 133 | |
---|
[1347] | 134 | void G4SPSEneDistribution::SetGradient(G4double gr) { |
---|
| 135 | grad = gr; |
---|
[816] | 136 | } |
---|
| 137 | |
---|
[1347] | 138 | void G4SPSEneDistribution::SetInterCept(G4double c) { |
---|
| 139 | cept = c; |
---|
[816] | 140 | } |
---|
| 141 | |
---|
[1347] | 142 | void G4SPSEneDistribution::UserEnergyHisto(G4ThreeVector input) { |
---|
| 143 | G4double ehi, val; |
---|
| 144 | ehi = input.x(); |
---|
| 145 | val = input.y(); |
---|
| 146 | if (verbosityLevel > 1) { |
---|
| 147 | G4cout << "In UserEnergyHisto" << G4endl; |
---|
| 148 | G4cout << " " << ehi << " " << val << G4endl; |
---|
| 149 | } |
---|
| 150 | UDefEnergyH.InsertValues(ehi, val); |
---|
| 151 | Emax = ehi; |
---|
[816] | 152 | } |
---|
| 153 | |
---|
[1347] | 154 | void G4SPSEneDistribution::ArbEnergyHisto(G4ThreeVector input) { |
---|
| 155 | G4double ehi, val; |
---|
| 156 | ehi = input.x(); |
---|
| 157 | val = input.y(); |
---|
| 158 | if (verbosityLevel > 1) { |
---|
| 159 | G4cout << "In ArbEnergyHisto" << G4endl; |
---|
| 160 | G4cout << " " << ehi << " " << val << G4endl; |
---|
| 161 | } |
---|
| 162 | ArbEnergyH.InsertValues(ehi, val); |
---|
[816] | 163 | } |
---|
| 164 | |
---|
[1347] | 165 | void G4SPSEneDistribution::ArbEnergyHistoFile(G4String filename) { |
---|
| 166 | std::ifstream infile(filename, std::ios::in); |
---|
| 167 | if (!infile) |
---|
| 168 | G4Exception("Unable to open the histo ASCII file"); |
---|
| 169 | G4double ehi, val; |
---|
| 170 | while (infile >> ehi >> val) { |
---|
| 171 | ArbEnergyH.InsertValues(ehi, val); |
---|
| 172 | } |
---|
[816] | 173 | } |
---|
| 174 | |
---|
[1347] | 175 | void G4SPSEneDistribution::EpnEnergyHisto(G4ThreeVector input) { |
---|
| 176 | G4double ehi, val; |
---|
| 177 | ehi = input.x(); |
---|
| 178 | val = input.y(); |
---|
| 179 | if (verbosityLevel > 1) { |
---|
| 180 | G4cout << "In EpnEnergyHisto" << G4endl; |
---|
| 181 | G4cout << " " << ehi << " " << val << G4endl; |
---|
| 182 | } |
---|
| 183 | EpnEnergyH.InsertValues(ehi, val); |
---|
| 184 | Emax = ehi; |
---|
| 185 | Epnflag = true; |
---|
| 186 | } |
---|
[816] | 187 | |
---|
[1347] | 188 | void G4SPSEneDistribution::Calculate() { |
---|
| 189 | if (EnergyDisType == "Cdg") |
---|
| 190 | CalculateCdgSpectrum(); |
---|
| 191 | else if (EnergyDisType == "Bbody") |
---|
| 192 | CalculateBbodySpectrum(); |
---|
| 193 | } |
---|
| 194 | |
---|
| 195 | void G4SPSEneDistribution::CalculateCdgSpectrum() { |
---|
| 196 | // This uses the spectrum from The INTEGRAL Mass Model (TIMM) |
---|
| 197 | // to generate a Cosmic Diffuse X/gamma ray spectrum. |
---|
| 198 | G4double pfact[2] = { 8.5, 112 }; |
---|
| 199 | G4double spind[2] = { 1.4, 2.3 }; |
---|
| 200 | G4double ene_line[3] = { 1. * keV, 18. * keV, 1E6 * keV }; |
---|
| 201 | G4int n_par; |
---|
| 202 | |
---|
| 203 | ene_line[0] = Emin; |
---|
| 204 | if (Emin < 18 * keV) { |
---|
| 205 | n_par = 2; |
---|
| 206 | ene_line[2] = Emax; |
---|
| 207 | if (Emax < 18 * keV) { |
---|
| 208 | n_par = 1; |
---|
| 209 | ene_line[1] = Emax; |
---|
| 210 | } |
---|
| 211 | } else { |
---|
| 212 | n_par = 1; |
---|
| 213 | pfact[0] = 112.; |
---|
| 214 | spind[0] = 2.3; |
---|
| 215 | ene_line[1] = Emax; |
---|
[816] | 216 | } |
---|
| 217 | |
---|
[1347] | 218 | // Create a cumulative histogram. |
---|
| 219 | CDGhist[0] = 0.; |
---|
| 220 | G4double omalpha; |
---|
| 221 | G4int i = 0; |
---|
| 222 | |
---|
| 223 | while (i < n_par) { |
---|
| 224 | omalpha = 1. - spind[i]; |
---|
| 225 | CDGhist[i + 1] = CDGhist[i] + (pfact[i] / omalpha) * (std::pow( |
---|
| 226 | ene_line[i + 1] / keV, omalpha) - std::pow(ene_line[i] / keV, |
---|
| 227 | omalpha)); |
---|
| 228 | i++; |
---|
| 229 | } |
---|
| 230 | |
---|
| 231 | // Normalise histo and divide by 1000 to make MeV. |
---|
| 232 | i = 0; |
---|
| 233 | while (i < n_par) { |
---|
| 234 | CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[n_par]; |
---|
| 235 | // G4cout << CDGhist[i] << CDGhist[n_par] << G4endl; |
---|
| 236 | i++; |
---|
| 237 | } |
---|
[816] | 238 | } |
---|
| 239 | |
---|
[1347] | 240 | void G4SPSEneDistribution::CalculateBbodySpectrum() { |
---|
| 241 | // create bbody spectrum |
---|
| 242 | // Proved very hard to integrate indefinitely, so different |
---|
| 243 | // method. User inputs emin, emax and T. These are used to |
---|
| 244 | // create a 10,000 bin histogram. |
---|
| 245 | // Use photon density spectrum = 2 nu**2/c**2 * (std::exp(h nu/kT)-1) |
---|
| 246 | // = 2 E**2/h**2c**2 times the exponential |
---|
| 247 | G4double erange = Emax - Emin; |
---|
| 248 | G4double steps = erange / 10000.; |
---|
| 249 | G4double Bbody_y[10000]; |
---|
| 250 | G4double k = 8.6181e-11; //Boltzmann const in MeV/K |
---|
| 251 | G4double h = 4.1362e-21; // Plancks const in MeV s |
---|
| 252 | G4double c = 3e8; // Speed of light |
---|
| 253 | G4double h2 = h * h; |
---|
| 254 | G4double c2 = c * c; |
---|
| 255 | G4int count = 0; |
---|
| 256 | G4double sum = 0.; |
---|
| 257 | BBHist[0] = 0.; |
---|
| 258 | while (count < 10000) { |
---|
| 259 | Bbody_x[count] = Emin + G4double(count * steps); |
---|
| 260 | Bbody_y[count] = (2. * std::pow(Bbody_x[count], 2.)) / (h2 * c2 |
---|
| 261 | * (std::exp(Bbody_x[count] / (k * Temp)) - 1.)); |
---|
| 262 | sum = sum + Bbody_y[count]; |
---|
| 263 | BBHist[count + 1] = BBHist[count] + Bbody_y[count]; |
---|
| 264 | count++; |
---|
| 265 | } |
---|
[816] | 266 | |
---|
[1347] | 267 | Bbody_x[10000] = Emax; |
---|
| 268 | // Normalise cumulative histo. |
---|
| 269 | count = 0; |
---|
| 270 | while (count < 10001) { |
---|
| 271 | BBHist[count] = BBHist[count] / sum; |
---|
| 272 | count++; |
---|
| 273 | } |
---|
[816] | 274 | } |
---|
| 275 | |
---|
[1347] | 276 | void G4SPSEneDistribution::InputEnergySpectra(G4bool value) { |
---|
| 277 | // Allows user to specifiy spectrum is momentum |
---|
| 278 | EnergySpec = value; // false if momentum |
---|
| 279 | if (verbosityLevel > 1) |
---|
| 280 | G4cout << "EnergySpec has value " << EnergySpec << G4endl; |
---|
[816] | 281 | } |
---|
| 282 | |
---|
[1347] | 283 | void G4SPSEneDistribution::InputDifferentialSpectra(G4bool value) { |
---|
| 284 | // Allows user to specify integral or differential spectra |
---|
| 285 | DiffSpec = value; // true = differential, false = integral |
---|
| 286 | if (verbosityLevel > 1) |
---|
| 287 | G4cout << "Diffspec has value " << DiffSpec << G4endl; |
---|
[816] | 288 | } |
---|
| 289 | |
---|
[1347] | 290 | void G4SPSEneDistribution::ArbInterpolate(G4String IType) { |
---|
| 291 | if (EnergyDisType != "Arb") |
---|
| 292 | G4cout << "Error: this is for arbitrary distributions" << G4endl; |
---|
| 293 | IntType = IType; |
---|
| 294 | ArbEmax = ArbEnergyH.GetMaxLowEdgeEnergy(); |
---|
| 295 | ArbEmin = ArbEnergyH.GetMinLowEdgeEnergy(); |
---|
[816] | 296 | |
---|
[1347] | 297 | // Now interpolate points |
---|
| 298 | if (IntType == "Lin") |
---|
| 299 | LinearInterpolation(); |
---|
| 300 | if (IntType == "Log") |
---|
| 301 | LogInterpolation(); |
---|
| 302 | if (IntType == "Exp") |
---|
| 303 | ExpInterpolation(); |
---|
| 304 | if (IntType == "Spline") |
---|
| 305 | SplineInterpolation(); |
---|
[816] | 306 | } |
---|
| 307 | |
---|
[1347] | 308 | void G4SPSEneDistribution::LinearInterpolation() { |
---|
| 309 | // Method to do linear interpolation on the Arb points |
---|
| 310 | // Calculate equation of each line segment, max 1024. |
---|
| 311 | // Calculate Area under each segment |
---|
| 312 | // Create a cumulative array which is then normalised Arb_Cum_Area |
---|
[816] | 313 | |
---|
[1347] | 314 | G4double Area_seg[1024]; // Stores area under each segment |
---|
| 315 | G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024]; |
---|
| 316 | G4int i, count; |
---|
| 317 | G4int maxi = ArbEnergyH.GetVectorLength(); |
---|
| 318 | for (i = 0; i < maxi; i++) { |
---|
| 319 | Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); |
---|
| 320 | Arb_y[i] = ArbEnergyH(size_t(i)); |
---|
[816] | 321 | } |
---|
[1347] | 322 | // Points are now in x,y arrays. If the spectrum is integral it has to be |
---|
| 323 | // made differential and if momentum it has to be made energy. |
---|
| 324 | if (DiffSpec == false) { |
---|
| 325 | // Converts integral point-wise spectra to Differential |
---|
| 326 | for (count = 0; count < maxi - 1; count++) { |
---|
| 327 | Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) |
---|
| 328 | / (Arb_x[count + 1] - Arb_x[count]); |
---|
| 329 | } |
---|
| 330 | maxi--; |
---|
[816] | 331 | } |
---|
[1347] | 332 | // |
---|
| 333 | if (EnergySpec == false) { |
---|
| 334 | // change currently stored values (emin etc) which are actually momenta |
---|
| 335 | // to energies. |
---|
| 336 | if (particle_definition == NULL) |
---|
| 337 | G4cout << "Error: particle not defined" << G4endl; |
---|
| 338 | else { |
---|
| 339 | // Apply Energy**2 = p**2c**2 + m0**2c**4 |
---|
| 340 | // p should be entered as E/c i.e. without the division by c |
---|
| 341 | // being done - energy equivalent. |
---|
| 342 | G4double mass = particle_definition->GetPDGMass(); |
---|
| 343 | // convert point to energy unit and its value to per energy unit |
---|
| 344 | G4double total_energy; |
---|
| 345 | for (count = 0; count < maxi; count++) { |
---|
| 346 | total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass |
---|
| 347 | * mass)); // total energy |
---|
| 348 | |
---|
| 349 | Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy; |
---|
| 350 | Arb_x[count] = total_energy - mass; // kinetic energy |
---|
| 351 | } |
---|
| 352 | } |
---|
[816] | 353 | } |
---|
[1347] | 354 | // |
---|
| 355 | i = 1; |
---|
| 356 | Arb_grad[0] = 0.; |
---|
| 357 | Arb_cept[0] = 0.; |
---|
| 358 | Area_seg[0] = 0.; |
---|
| 359 | Arb_Cum_Area[0] = 0.; |
---|
| 360 | while (i < maxi) { |
---|
| 361 | // calc gradient and intercept for each segment |
---|
| 362 | Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[i - 1]); |
---|
| 363 | if (verbosityLevel == 2) |
---|
| 364 | G4cout << Arb_grad[i] << G4endl; |
---|
| 365 | if (Arb_grad[i] > 0.) { |
---|
| 366 | if (verbosityLevel == 2) |
---|
| 367 | G4cout << "Arb_grad is positive" << G4endl; |
---|
| 368 | Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]); |
---|
| 369 | } else if (Arb_grad[i] < 0.) { |
---|
| 370 | if (verbosityLevel == 2) |
---|
| 371 | G4cout << "Arb_grad is negative" << G4endl; |
---|
| 372 | Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]); |
---|
| 373 | } else { |
---|
| 374 | if (verbosityLevel == 2) |
---|
| 375 | G4cout << "Arb_grad is 0." << G4endl; |
---|
| 376 | Arb_cept[i] = Arb_y[i]; |
---|
| 377 | } |
---|
[816] | 378 | |
---|
[1347] | 379 | Area_seg[i] = ((Arb_grad[i] / 2) * (Arb_x[i] * Arb_x[i] - Arb_x[i - 1] |
---|
| 380 | * Arb_x[i - 1]) + Arb_cept[i] * (Arb_x[i] - Arb_x[i - 1])); |
---|
| 381 | Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i]; |
---|
| 382 | sum = sum + Area_seg[i]; |
---|
| 383 | if (verbosityLevel == 2) |
---|
| 384 | G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i] |
---|
| 385 | << G4endl; |
---|
| 386 | i++; |
---|
| 387 | } |
---|
[816] | 388 | |
---|
[1347] | 389 | i = 0; |
---|
| 390 | while (i < maxi) { |
---|
| 391 | Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation |
---|
| 392 | IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); |
---|
| 393 | i++; |
---|
| 394 | } |
---|
[816] | 395 | |
---|
[1347] | 396 | // now scale the ArbEnergyH, needed by Probability() |
---|
| 397 | ArbEnergyH.ScaleVector(1., 1./sum); |
---|
| 398 | |
---|
| 399 | if (verbosityLevel >= 1) { |
---|
| 400 | G4cout << "Leaving LinearInterpolation" << G4endl; |
---|
| 401 | ArbEnergyH.DumpValues(); |
---|
| 402 | IPDFArbEnergyH.DumpValues(); |
---|
| 403 | } |
---|
[816] | 404 | } |
---|
| 405 | |
---|
[1347] | 406 | void G4SPSEneDistribution::LogInterpolation() { |
---|
| 407 | // Interpolation based on Logarithmic equations |
---|
| 408 | // Generate equations of line segments |
---|
| 409 | // y = Ax**alpha => log y = alpha*logx + logA |
---|
| 410 | // Find area under line segments |
---|
| 411 | // create normalised, cumulative array Arb_Cum_Area |
---|
| 412 | G4double Area_seg[1024]; // Stores area under each segment |
---|
| 413 | G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024]; |
---|
| 414 | G4int i, count; |
---|
| 415 | G4int maxi = ArbEnergyH.GetVectorLength(); |
---|
| 416 | for (i = 0; i < maxi; i++) { |
---|
| 417 | Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); |
---|
| 418 | Arb_y[i] = ArbEnergyH(size_t(i)); |
---|
| 419 | } |
---|
| 420 | // Points are now in x,y arrays. If the spectrum is integral it has to be |
---|
| 421 | // made differential and if momentum it has to be made energy. |
---|
| 422 | if (DiffSpec == false) { |
---|
| 423 | // Converts integral point-wise spectra to Differential |
---|
| 424 | for (count = 0; count < maxi - 1; count++) { |
---|
| 425 | Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) |
---|
| 426 | / (Arb_x[count + 1] - Arb_x[count]); |
---|
| 427 | } |
---|
| 428 | maxi--; |
---|
| 429 | } |
---|
| 430 | // |
---|
| 431 | if (EnergySpec == false) { |
---|
| 432 | // change currently stored values (emin etc) which are actually momenta |
---|
| 433 | // to energies. |
---|
| 434 | if (particle_definition == NULL) |
---|
| 435 | G4cout << "Error: particle not defined" << G4endl; |
---|
| 436 | else { |
---|
| 437 | // Apply Energy**2 = p**2c**2 + m0**2c**4 |
---|
| 438 | // p should be entered as E/c i.e. without the division by c |
---|
| 439 | // being done - energy equivalent. |
---|
| 440 | G4double mass = particle_definition->GetPDGMass(); |
---|
| 441 | // convert point to energy unit and its value to per energy unit |
---|
| 442 | G4double total_energy; |
---|
| 443 | for (count = 0; count < maxi; count++) { |
---|
| 444 | total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass |
---|
| 445 | * mass)); // total energy |
---|
[816] | 446 | |
---|
[1347] | 447 | Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy; |
---|
| 448 | Arb_x[count] = total_energy - mass; // kinetic energy |
---|
| 449 | } |
---|
| 450 | } |
---|
[816] | 451 | } |
---|
[1347] | 452 | // |
---|
| 453 | i = 1; |
---|
| 454 | Arb_alpha[0] = 0.; |
---|
| 455 | Arb_Const[0] = 0.; |
---|
| 456 | Area_seg[0] = 0.; |
---|
| 457 | Arb_Cum_Area[0] = 0.; |
---|
| 458 | if (Arb_x[0] <= 0. || Arb_y[0] <= 0.) { |
---|
| 459 | G4cout << "You should not use log interpolation with points <= 0." |
---|
| 460 | << G4endl; |
---|
| 461 | G4cout << "These will be changed to 1e-20, which may cause problems" |
---|
| 462 | << G4endl; |
---|
| 463 | if (Arb_x[0] <= 0.) |
---|
| 464 | Arb_x[0] = 1e-20; |
---|
| 465 | if (Arb_y[0] <= 0.) |
---|
| 466 | Arb_y[0] = 1e-20; |
---|
| 467 | } |
---|
[816] | 468 | |
---|
[1347] | 469 | G4double alp; |
---|
| 470 | while (i < maxi) { |
---|
| 471 | // Incase points are negative or zero |
---|
| 472 | if (Arb_x[i] <= 0. || Arb_y[i] <= 0.) { |
---|
| 473 | G4cout << "You should not use log interpolation with points <= 0." |
---|
| 474 | << G4endl; |
---|
| 475 | G4cout |
---|
| 476 | << "These will be changed to 1e-20, which may cause problems" |
---|
| 477 | << G4endl; |
---|
| 478 | if (Arb_x[i] <= 0.) |
---|
| 479 | Arb_x[i] = 1e-20; |
---|
| 480 | if (Arb_y[i] <= 0.) |
---|
| 481 | Arb_y[i] = 1e-20; |
---|
| 482 | } |
---|
| 483 | |
---|
| 484 | Arb_alpha[i] = (std::log10(Arb_y[i]) - std::log10(Arb_y[i - 1])) |
---|
| 485 | / (std::log10(Arb_x[i]) - std::log10(Arb_x[i - 1])); |
---|
| 486 | Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[i], Arb_alpha[i])); |
---|
| 487 | alp = Arb_alpha[i] + 1; |
---|
| 488 | if (alp == 0.) { |
---|
| 489 | Area_seg[i] = Arb_Const[i] * (std::log(Arb_x[i]) - std::log(Arb_x[i - 1])); |
---|
| 490 | } else { |
---|
| 491 | Area_seg[i] = (Arb_Const[i] / alp) * (std::pow(Arb_x[i], alp) |
---|
| 492 | - std::pow(Arb_x[i - 1], alp)); |
---|
| 493 | } |
---|
| 494 | sum = sum + Area_seg[i]; |
---|
| 495 | Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i]; |
---|
| 496 | if (verbosityLevel == 2) |
---|
| 497 | G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl; |
---|
| 498 | i++; |
---|
| 499 | } |
---|
| 500 | |
---|
| 501 | i = 0; |
---|
| 502 | while (i < maxi) { |
---|
| 503 | Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; |
---|
| 504 | IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); |
---|
| 505 | i++; |
---|
| 506 | } |
---|
| 507 | |
---|
| 508 | // now scale the ArbEnergyH, needed by Probability() |
---|
| 509 | ArbEnergyH.ScaleVector(1., 1./sum); |
---|
| 510 | |
---|
| 511 | if (verbosityLevel >= 1) |
---|
| 512 | G4cout << "Leaving LogInterpolation " << G4endl; |
---|
[816] | 513 | } |
---|
| 514 | |
---|
[1347] | 515 | void G4SPSEneDistribution::ExpInterpolation() { |
---|
| 516 | // Interpolation based on Exponential equations |
---|
| 517 | // Generate equations of line segments |
---|
| 518 | // y = Ae**-(x/e0) => ln y = -x/e0 + lnA |
---|
| 519 | // Find area under line segments |
---|
| 520 | // create normalised, cumulative array Arb_Cum_Area |
---|
| 521 | G4double Area_seg[1024]; // Stores area under each segment |
---|
| 522 | G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024]; |
---|
| 523 | G4int i, count; |
---|
| 524 | G4int maxi = ArbEnergyH.GetVectorLength(); |
---|
| 525 | for (i = 0; i < maxi; i++) { |
---|
| 526 | Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); |
---|
| 527 | Arb_y[i] = ArbEnergyH(size_t(i)); |
---|
[816] | 528 | } |
---|
[1347] | 529 | // Points are now in x,y arrays. If the spectrum is integral it has to be |
---|
| 530 | // made differential and if momentum it has to be made energy. |
---|
| 531 | if (DiffSpec == false) { |
---|
| 532 | // Converts integral point-wise spectra to Differential |
---|
| 533 | for (count = 0; count < maxi - 1; count++) { |
---|
| 534 | Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) |
---|
| 535 | / (Arb_x[count + 1] - Arb_x[count]); |
---|
| 536 | } |
---|
| 537 | maxi--; |
---|
[816] | 538 | } |
---|
[1347] | 539 | // |
---|
| 540 | if (EnergySpec == false) { |
---|
| 541 | // change currently stored values (emin etc) which are actually momenta |
---|
| 542 | // to energies. |
---|
| 543 | if (particle_definition == NULL) |
---|
| 544 | G4cout << "Error: particle not defined" << G4endl; |
---|
| 545 | else { |
---|
| 546 | // Apply Energy**2 = p**2c**2 + m0**2c**4 |
---|
| 547 | // p should be entered as E/c i.e. without the division by c |
---|
| 548 | // being done - energy equivalent. |
---|
| 549 | G4double mass = particle_definition->GetPDGMass(); |
---|
| 550 | // convert point to energy unit and its value to per energy unit |
---|
| 551 | G4double total_energy; |
---|
| 552 | for (count = 0; count < maxi; count++) { |
---|
| 553 | total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass |
---|
| 554 | * mass)); // total energy |
---|
| 555 | |
---|
| 556 | Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy; |
---|
| 557 | Arb_x[count] = total_energy - mass; // kinetic energy |
---|
| 558 | } |
---|
| 559 | } |
---|
| 560 | } |
---|
| 561 | // |
---|
| 562 | i = 1; |
---|
| 563 | Arb_ezero[0] = 0.; |
---|
| 564 | Arb_Const[0] = 0.; |
---|
| 565 | Area_seg[0] = 0.; |
---|
| 566 | Arb_Cum_Area[0] = 0.; |
---|
| 567 | while (i < maxi) { |
---|
| 568 | G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i - 1]); |
---|
| 569 | if (test > 0. || test < 0.) { |
---|
| 570 | Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1]) / (std::log(Arb_y[i]) |
---|
| 571 | - std::log(Arb_y[i - 1])); |
---|
| 572 | Arb_Const[i] = Arb_y[i] / (std::exp(-Arb_x[i] / Arb_ezero[i])); |
---|
| 573 | Area_seg[i] = -(Arb_Const[i] * Arb_ezero[i]) * (std::exp(-Arb_x[i] |
---|
| 574 | / Arb_ezero[i]) - std::exp(-Arb_x[i - 1] / Arb_ezero[i])); |
---|
| 575 | } else { |
---|
| 576 | G4cout << "Flat line segment: problem" << G4endl; |
---|
| 577 | Arb_ezero[i] = 0.; |
---|
| 578 | Arb_Const[i] = 0.; |
---|
| 579 | Area_seg[i] = 0.; |
---|
| 580 | } |
---|
| 581 | sum = sum + Area_seg[i]; |
---|
| 582 | Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i]; |
---|
| 583 | if (verbosityLevel == 2) |
---|
| 584 | G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl; |
---|
| 585 | i++; |
---|
| 586 | } |
---|
| 587 | |
---|
| 588 | i = 0; |
---|
| 589 | while (i < maxi) { |
---|
| 590 | Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; |
---|
| 591 | IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); |
---|
| 592 | i++; |
---|
| 593 | } |
---|
| 594 | |
---|
| 595 | // now scale the ArbEnergyH, needed by Probability() |
---|
| 596 | ArbEnergyH.ScaleVector(1., 1./sum); |
---|
| 597 | |
---|
| 598 | if (verbosityLevel >= 1) |
---|
| 599 | G4cout << "Leaving ExpInterpolation " << G4endl; |
---|
[816] | 600 | } |
---|
| 601 | |
---|
[1347] | 602 | void G4SPSEneDistribution::SplineInterpolation() { |
---|
| 603 | // Interpolation using Splines. |
---|
| 604 | // Create Normalised arrays, make x 0->1 and y hold |
---|
| 605 | // the function (Energy) |
---|
| 606 | // |
---|
| 607 | // Current method based on the above will not work in all cases. |
---|
| 608 | // new method is implemented below. |
---|
| 609 | |
---|
| 610 | G4double sum, Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024]; |
---|
| 611 | G4int i, count; |
---|
[816] | 612 | |
---|
[1347] | 613 | G4int maxi = ArbEnergyH.GetVectorLength(); |
---|
| 614 | for (i = 0; i < maxi; i++) { |
---|
| 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 | // Converts integral point-wise spectra to Differential |
---|
| 622 | for (count = 0; count < maxi - 1; count++) { |
---|
| 623 | Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) |
---|
| 624 | / (Arb_x[count + 1] - Arb_x[count]); |
---|
| 625 | } |
---|
| 626 | maxi--; |
---|
| 627 | } |
---|
| 628 | // |
---|
| 629 | if (EnergySpec == false) { |
---|
| 630 | // change currently stored values (emin etc) which are actually momenta |
---|
| 631 | // to energies. |
---|
| 632 | if (particle_definition == NULL) |
---|
| 633 | G4cout << "Error: particle not defined" << G4endl; |
---|
| 634 | else { |
---|
| 635 | // Apply Energy**2 = p**2c**2 + m0**2c**4 |
---|
| 636 | // p should be entered as E/c i.e. without the division by c |
---|
| 637 | // being done - energy equivalent. |
---|
| 638 | G4double mass = particle_definition->GetPDGMass(); |
---|
| 639 | // convert point to energy unit and its value to per energy unit |
---|
| 640 | G4double total_energy; |
---|
| 641 | for (count = 0; count < maxi; count++) { |
---|
| 642 | total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass |
---|
| 643 | * mass)); // total energy |
---|
| 644 | |
---|
| 645 | Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy; |
---|
| 646 | Arb_x[count] = total_energy - mass; // kinetic energy |
---|
| 647 | } |
---|
| 648 | } |
---|
| 649 | } |
---|
| 650 | |
---|
| 651 | // |
---|
| 652 | i = 1; |
---|
| 653 | Arb_Cum_Area[0] = 0.; |
---|
| 654 | sum = 0.; |
---|
| 655 | Splinetemp = new G4DataInterpolation(Arb_x, Arb_y, maxi, 0., 0.); |
---|
| 656 | G4double ei[101],prob[101]; |
---|
| 657 | while (i < maxi) { |
---|
| 658 | // 100 step per segment for the integration of area |
---|
| 659 | G4double de = (Arb_x[i] - Arb_x[i - 1])/100.; |
---|
| 660 | G4double area = 0.; |
---|
| 661 | |
---|
| 662 | for (count = 0; count < 101; count++) { |
---|
| 663 | ei[count] = Arb_x[i - 1] + de*count ; |
---|
| 664 | prob[count] = Splinetemp->CubicSplineInterpolation(ei[count]); |
---|
| 665 | if (prob[count] < 0.) { |
---|
| 666 | G4cout << "Warning: G4DataInterpolation returns value < 0 " << prob[count] <<" "<<ei[count]<< G4endl; |
---|
| 667 | G4Exception(" Please use an alternative method, e.g. Lin, for interpolation"); |
---|
| 668 | } |
---|
| 669 | area += prob[count]*de; |
---|
| 670 | } |
---|
| 671 | Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + area; |
---|
| 672 | sum += area; |
---|
| 673 | |
---|
| 674 | prob[0] = prob[0]/(area/de); |
---|
| 675 | for (count = 1; count < 100; count++) |
---|
| 676 | prob[count] = prob[count-1] + prob[count]/(area/de); |
---|
| 677 | |
---|
| 678 | SplineInt[i] = new G4DataInterpolation(prob, ei, 101, 0., 0.); |
---|
| 679 | // note i start from 1! |
---|
| 680 | i++; |
---|
| 681 | } |
---|
| 682 | i = 0; |
---|
| 683 | while (i < maxi) { |
---|
| 684 | Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation |
---|
| 685 | IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); |
---|
| 686 | i++; |
---|
| 687 | } |
---|
| 688 | // now scale the ArbEnergyH, needed by Probability() |
---|
| 689 | ArbEnergyH.ScaleVector(1., 1./sum); |
---|
| 690 | |
---|
| 691 | if (verbosityLevel > 0) |
---|
| 692 | G4cout << "Leaving SplineInterpolation " << G4endl; |
---|
[816] | 693 | } |
---|
| 694 | |
---|
[1347] | 695 | void G4SPSEneDistribution::GenerateMonoEnergetic() { |
---|
| 696 | // Method to generate MonoEnergetic particles. |
---|
| 697 | particle_energy = MonoEnergy; |
---|
[816] | 698 | } |
---|
| 699 | |
---|
[1347] | 700 | void G4SPSEneDistribution::GenerateGaussEnergies() { |
---|
| 701 | // Method to generate Gaussian particles. |
---|
| 702 | particle_energy = G4RandGauss::shoot(MonoEnergy,SE); |
---|
| 703 | if (particle_energy < 0) particle_energy = 0.; |
---|
[816] | 704 | } |
---|
| 705 | |
---|
[1347] | 706 | void G4SPSEneDistribution::GenerateLinearEnergies(G4bool bArb = false) { |
---|
| 707 | G4double rndm; |
---|
| 708 | G4double emaxsq = std::pow(Emax, 2.); //Emax squared |
---|
| 709 | G4double eminsq = std::pow(Emin, 2.); //Emin squared |
---|
| 710 | G4double intersq = std::pow(cept, 2.); //cept squared |
---|
[816] | 711 | |
---|
[1347] | 712 | if (bArb) |
---|
| 713 | rndm = G4UniformRand(); |
---|
| 714 | else |
---|
| 715 | rndm = eneRndm->GenRandEnergy(); |
---|
[816] | 716 | |
---|
[1347] | 717 | G4double bracket = ((grad / 2.) * (emaxsq - eminsq) + cept * (Emax - Emin)); |
---|
| 718 | bracket = bracket * rndm; |
---|
| 719 | bracket = bracket + (grad / 2.) * eminsq + cept * Emin; |
---|
| 720 | // Now have a quad of form m/2 E**2 + cE - bracket = 0 |
---|
| 721 | bracket = -bracket; |
---|
| 722 | // G4cout << "BRACKET" << bracket << G4endl; |
---|
| 723 | if (grad != 0.) { |
---|
| 724 | G4double sqbrack = (intersq - 4 * (grad / 2.) * (bracket)); |
---|
| 725 | // G4cout << "SQBRACK" << sqbrack << G4endl; |
---|
| 726 | sqbrack = std::sqrt(sqbrack); |
---|
| 727 | G4double root1 = -cept + sqbrack; |
---|
| 728 | root1 = root1 / (2. * (grad / 2.)); |
---|
[816] | 729 | |
---|
[1347] | 730 | G4double root2 = -cept - sqbrack; |
---|
| 731 | root2 = root2 / (2. * (grad / 2.)); |
---|
[816] | 732 | |
---|
[1347] | 733 | // G4cout << root1 << " roots " << root2 << G4endl; |
---|
| 734 | |
---|
| 735 | if (root1 > Emin && root1 < Emax) |
---|
| 736 | particle_energy = root1; |
---|
| 737 | if (root2 > Emin && root2 < Emax) |
---|
| 738 | particle_energy = root2; |
---|
| 739 | } else if (grad == 0.) |
---|
| 740 | // have equation of form cE - bracket =0 |
---|
| 741 | particle_energy = bracket / cept; |
---|
| 742 | |
---|
| 743 | if (particle_energy < 0.) |
---|
| 744 | particle_energy = -particle_energy; |
---|
| 745 | |
---|
| 746 | if (verbosityLevel >= 1) |
---|
| 747 | G4cout << "Energy is " << particle_energy << G4endl; |
---|
[816] | 748 | } |
---|
| 749 | |
---|
[1347] | 750 | void G4SPSEneDistribution::GeneratePowEnergies(G4bool bArb = false) { |
---|
| 751 | // Method to generate particle energies distributed as |
---|
| 752 | // a power-law |
---|
[816] | 753 | |
---|
[1347] | 754 | G4double rndm; |
---|
| 755 | G4double emina, emaxa; |
---|
[816] | 756 | |
---|
[1347] | 757 | emina = std::pow(Emin, alpha + 1); |
---|
| 758 | emaxa = std::pow(Emax, alpha + 1); |
---|
[816] | 759 | |
---|
[1347] | 760 | if (bArb) |
---|
| 761 | rndm = G4UniformRand(); |
---|
| 762 | else |
---|
| 763 | rndm = eneRndm->GenRandEnergy(); |
---|
| 764 | |
---|
| 765 | if (alpha != -1.) { |
---|
| 766 | particle_energy = ((rndm * (emaxa - emina)) + emina); |
---|
| 767 | particle_energy = std::pow(particle_energy, (1. / (alpha + 1.))); |
---|
| 768 | } else { |
---|
| 769 | particle_energy = (std::log(Emin) + rndm * (std::log(Emax) - std::log( |
---|
| 770 | Emin))); |
---|
| 771 | particle_energy = std::exp(particle_energy); |
---|
| 772 | } |
---|
| 773 | if (verbosityLevel >= 1) |
---|
| 774 | G4cout << "Energy is " << particle_energy << G4endl; |
---|
[816] | 775 | } |
---|
| 776 | |
---|
[1347] | 777 | void G4SPSEneDistribution::GenerateBiasPowEnergies() { |
---|
| 778 | // Method to generate particle energies distributed as |
---|
| 779 | // in biased power-law and calculate its weight |
---|
[816] | 780 | |
---|
[1347] | 781 | G4double rndm; |
---|
| 782 | G4double emina, emaxa, emin, emax; |
---|
| 783 | |
---|
| 784 | G4double normal = 1. ; |
---|
| 785 | |
---|
| 786 | emin = Emin; |
---|
| 787 | emax = Emax; |
---|
| 788 | // if (EnergyDisType == "Arb") { |
---|
| 789 | // emin = ArbEmin; |
---|
| 790 | // emax = ArbEmax; |
---|
| 791 | //} |
---|
| 792 | |
---|
| 793 | rndm = eneRndm->GenRandEnergy(); |
---|
| 794 | |
---|
| 795 | if (biasalpha != -1.) { |
---|
| 796 | emina = std::pow(emin, biasalpha + 1); |
---|
| 797 | emaxa = std::pow(emax, biasalpha + 1); |
---|
| 798 | particle_energy = ((rndm * (emaxa - emina)) + emina); |
---|
| 799 | particle_energy = std::pow(particle_energy, (1. / (biasalpha + 1.))); |
---|
| 800 | normal = 1./(1+biasalpha) * (emaxa - emina); |
---|
| 801 | } else { |
---|
| 802 | particle_energy = (std::log(emin) + rndm * (std::log(emax) - std::log( |
---|
| 803 | emin))); |
---|
| 804 | particle_energy = std::exp(particle_energy); |
---|
| 805 | normal = std::log(emax) - std::log(emin) ; |
---|
| 806 | } |
---|
| 807 | weight = GetProbability(particle_energy) / (std::pow(particle_energy,biasalpha)/normal); |
---|
| 808 | |
---|
| 809 | if (verbosityLevel >= 1) |
---|
| 810 | G4cout << "Energy is " << particle_energy << G4endl; |
---|
[816] | 811 | } |
---|
| 812 | |
---|
[1347] | 813 | void G4SPSEneDistribution::GenerateExpEnergies(G4bool bArb = false) { |
---|
| 814 | // Method to generate particle energies distributed according |
---|
| 815 | // to an exponential curve. |
---|
| 816 | G4double rndm; |
---|
[816] | 817 | |
---|
[1347] | 818 | if (bArb) |
---|
| 819 | rndm = G4UniformRand(); |
---|
| 820 | else |
---|
| 821 | rndm = eneRndm->GenRandEnergy(); |
---|
[816] | 822 | |
---|
[1347] | 823 | particle_energy = -Ezero * (std::log(rndm * (std::exp(-Emax / Ezero) |
---|
| 824 | - std::exp(-Emin / Ezero)) + std::exp(-Emin / Ezero))); |
---|
| 825 | if (verbosityLevel >= 1) |
---|
| 826 | G4cout << "Energy is " << particle_energy << G4endl; |
---|
| 827 | } |
---|
[816] | 828 | |
---|
[1347] | 829 | void G4SPSEneDistribution::GenerateBremEnergies() { |
---|
| 830 | // Method to generate particle energies distributed according |
---|
| 831 | // to a Bremstrahlung equation of |
---|
| 832 | // form I = const*((kT)**1/2)*E*(e**(-E/kT)) |
---|
[816] | 833 | |
---|
[1347] | 834 | G4double rndm; |
---|
| 835 | rndm = eneRndm->GenRandEnergy(); |
---|
| 836 | G4double expmax, expmin, k; |
---|
[816] | 837 | |
---|
[1347] | 838 | k = 8.6181e-11; // Boltzmann's const in MeV/K |
---|
| 839 | G4double ksq = std::pow(k, 2.); // k squared |
---|
| 840 | G4double Tsq = std::pow(Temp, 2.); // Temp squared |
---|
[816] | 841 | |
---|
[1347] | 842 | expmax = std::exp(-Emax / (k * Temp)); |
---|
| 843 | expmin = std::exp(-Emin / (k * Temp)); |
---|
[816] | 844 | |
---|
[1347] | 845 | // If either expmax or expmin are zero then this will cause problems |
---|
| 846 | // Most probably this will be because T is too low or E is too high |
---|
[816] | 847 | |
---|
[1347] | 848 | if (expmax == 0.) |
---|
| 849 | G4cout << "*****EXPMAX=0. Choose different E's or Temp" << G4endl; |
---|
| 850 | if (expmin == 0.) |
---|
| 851 | G4cout << "*****EXPMIN=0. Choose different E's or Temp" << G4endl; |
---|
[816] | 852 | |
---|
[1347] | 853 | G4double tempvar = rndm * ((-k) * Temp * (Emax * expmax - Emin * expmin) |
---|
| 854 | - (ksq * Tsq * (expmax - expmin))); |
---|
[816] | 855 | |
---|
[1347] | 856 | G4double bigc = (tempvar - k * Temp * Emin * expmin - ksq * Tsq * expmin) |
---|
| 857 | / (-k * Temp); |
---|
[816] | 858 | |
---|
[1347] | 859 | // This gives an equation of form: Ee(-E/kT) + kTe(-E/kT) - C =0 |
---|
| 860 | // Solve this iteratively, step from Emin to Emax in 1000 steps |
---|
| 861 | // and take the best solution. |
---|
[816] | 862 | |
---|
[1347] | 863 | G4double erange = Emax - Emin; |
---|
| 864 | G4double steps = erange / 1000.; |
---|
| 865 | G4int i; |
---|
| 866 | G4double etest, diff, err; |
---|
[816] | 867 | |
---|
[1347] | 868 | err = 100000.; |
---|
[816] | 869 | |
---|
[1347] | 870 | for (i = 1; i < 1000; i++) { |
---|
| 871 | etest = Emin + (i - 1) * steps; |
---|
[816] | 872 | |
---|
[1347] | 873 | diff = etest * (std::exp(-etest / (k * Temp))) + k * Temp * (std::exp( |
---|
| 874 | -etest / (k * Temp))) - bigc; |
---|
[816] | 875 | |
---|
[1347] | 876 | if (diff < 0.) |
---|
| 877 | diff = -diff; |
---|
[816] | 878 | |
---|
[1347] | 879 | if (diff < err) { |
---|
| 880 | err = diff; |
---|
| 881 | particle_energy = etest; |
---|
| 882 | } |
---|
| 883 | } |
---|
| 884 | if (verbosityLevel >= 1) |
---|
| 885 | G4cout << "Energy is " << particle_energy << G4endl; |
---|
[816] | 886 | } |
---|
| 887 | |
---|
[1347] | 888 | void G4SPSEneDistribution::GenerateBbodyEnergies() { |
---|
| 889 | // BBody_x holds Energies, and BBHist holds the cumulative histo. |
---|
| 890 | // binary search to find correct bin then lin interpolation. |
---|
| 891 | // Use the earlier defined histogram + RandGeneral method to generate |
---|
| 892 | // random numbers following the histos distribution. |
---|
| 893 | G4double rndm; |
---|
| 894 | G4int nabove, nbelow = 0, middle; |
---|
| 895 | nabove = 10001; |
---|
| 896 | rndm = eneRndm->GenRandEnergy(); |
---|
[816] | 897 | |
---|
[1347] | 898 | // Binary search to find bin that rndm is in |
---|
| 899 | while (nabove - nbelow > 1) { |
---|
| 900 | middle = (nabove + nbelow) / 2; |
---|
| 901 | if (rndm == BBHist[middle]) |
---|
| 902 | break; |
---|
| 903 | if (rndm < BBHist[middle]) |
---|
| 904 | nabove = middle; |
---|
| 905 | else |
---|
| 906 | nbelow = middle; |
---|
[816] | 907 | } |
---|
| 908 | |
---|
[1347] | 909 | // Now interpolate in that bin to find the correct output value. |
---|
| 910 | G4double x1, x2, y1, y2, m, q; |
---|
| 911 | x1 = Bbody_x[nbelow]; |
---|
| 912 | x2 = Bbody_x[nbelow + 1]; |
---|
| 913 | y1 = BBHist[nbelow]; |
---|
| 914 | y2 = BBHist[nbelow + 1]; |
---|
| 915 | m = (y2 - y1) / (x2 - x1); |
---|
| 916 | q = y1 - m * x1; |
---|
| 917 | |
---|
| 918 | particle_energy = (rndm - q) / m; |
---|
| 919 | |
---|
| 920 | if (verbosityLevel >= 1) { |
---|
| 921 | G4cout << "Energy is " << particle_energy << G4endl; |
---|
[816] | 922 | } |
---|
[1347] | 923 | } |
---|
[816] | 924 | |
---|
[1347] | 925 | void G4SPSEneDistribution::GenerateCdgEnergies() { |
---|
| 926 | // Gen random numbers, compare with values in cumhist |
---|
| 927 | // to find appropriate part of spectrum and then |
---|
| 928 | // generate energy in the usual inversion way. |
---|
| 929 | // G4double pfact[2] = {8.5, 112}; |
---|
| 930 | // G4double spind[2] = {1.4, 2.3}; |
---|
| 931 | // G4double ene_line[3] = {1., 18., 1E6}; |
---|
| 932 | G4double rndm, rndm2; |
---|
| 933 | G4double ene_line[3]; |
---|
| 934 | G4double omalpha[2]; |
---|
| 935 | if (Emin < 18 * keV && Emax < 18 * keV) { |
---|
| 936 | omalpha[0] = 1. - 1.4; |
---|
| 937 | ene_line[0] = Emin; |
---|
| 938 | ene_line[1] = Emax; |
---|
[816] | 939 | } |
---|
[1347] | 940 | if (Emin < 18 * keV && Emax > 18 * keV) { |
---|
| 941 | omalpha[0] = 1. - 1.4; |
---|
| 942 | omalpha[1] = 1. - 2.3; |
---|
| 943 | ene_line[0] = Emin; |
---|
| 944 | ene_line[1] = 18. * keV; |
---|
| 945 | ene_line[2] = Emax; |
---|
[816] | 946 | } |
---|
[1347] | 947 | if (Emin > 18 * keV) { |
---|
| 948 | omalpha[0] = 1. - 2.3; |
---|
| 949 | ene_line[0] = Emin; |
---|
| 950 | ene_line[1] = Emax; |
---|
| 951 | } |
---|
| 952 | rndm = eneRndm->GenRandEnergy(); |
---|
| 953 | rndm2 = eneRndm->GenRandEnergy(); |
---|
| 954 | |
---|
| 955 | G4int i = 0; |
---|
| 956 | while (rndm >= CDGhist[i]) { |
---|
| 957 | i++; |
---|
| 958 | } |
---|
| 959 | // Generate final energy. |
---|
| 960 | particle_energy = (std::pow(ene_line[i - 1], omalpha[i - 1]) + (std::pow( |
---|
| 961 | ene_line[i], omalpha[i - 1]) - std::pow(ene_line[i - 1], omalpha[i |
---|
| 962 | - 1])) * rndm2); |
---|
| 963 | particle_energy = std::pow(particle_energy, (1. / omalpha[i - 1])); |
---|
| 964 | |
---|
| 965 | if (verbosityLevel >= 1) |
---|
| 966 | G4cout << "Energy is " << particle_energy << G4endl; |
---|
[816] | 967 | } |
---|
| 968 | |
---|
[1347] | 969 | void G4SPSEneDistribution::GenUserHistEnergies() { |
---|
| 970 | // Histograms are DIFFERENTIAL. |
---|
| 971 | // G4cout << "In GenUserHistEnergies " << G4endl; |
---|
| 972 | if (IPDFEnergyExist == false) { |
---|
| 973 | G4int ii; |
---|
| 974 | G4int maxbin = G4int(UDefEnergyH.GetVectorLength()); |
---|
| 975 | G4double bins[1024], vals[1024], sum; |
---|
| 976 | sum = 0.; |
---|
| 977 | |
---|
| 978 | if ((EnergySpec == false) && (particle_definition == NULL)) |
---|
| 979 | G4cout << "Error: particle definition is NULL" << G4endl; |
---|
| 980 | |
---|
| 981 | if (maxbin > 1024) { |
---|
| 982 | G4cout << "Maxbin > 1024" << G4endl; |
---|
| 983 | G4cout << "Setting maxbin to 1024, other bins are lost" << G4endl; |
---|
| 984 | } |
---|
| 985 | |
---|
| 986 | if (DiffSpec == false) |
---|
| 987 | G4cout << "Histograms are Differential!!! " << G4endl; |
---|
| 988 | else { |
---|
| 989 | bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0)); |
---|
| 990 | vals[0] = UDefEnergyH(size_t(0)); |
---|
| 991 | sum = vals[0]; |
---|
| 992 | for (ii = 1; ii < maxbin; ii++) { |
---|
| 993 | bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii)); |
---|
| 994 | vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1]; |
---|
| 995 | sum = sum + UDefEnergyH(size_t(ii)); |
---|
| 996 | } |
---|
| 997 | } |
---|
| 998 | |
---|
| 999 | if (EnergySpec == false) { |
---|
| 1000 | G4double mass = particle_definition->GetPDGMass(); |
---|
| 1001 | // multiply the function (vals) up by the bin width |
---|
| 1002 | // to make the function counts/s (i.e. get rid of momentum |
---|
| 1003 | // dependence). |
---|
| 1004 | for (ii = 1; ii < maxbin; ii++) { |
---|
| 1005 | vals[ii] = vals[ii] * (bins[ii] - bins[ii - 1]); |
---|
| 1006 | } |
---|
| 1007 | // Put energy bins into new histo, plus divide by energy bin width |
---|
| 1008 | // to make evals counts/s/energy |
---|
| 1009 | for (ii = 0; ii < maxbin; ii++) { |
---|
| 1010 | bins[ii] = std::sqrt((bins[ii] * bins[ii]) + (mass * mass)) |
---|
| 1011 | - mass; //kinetic energy |
---|
| 1012 | } |
---|
| 1013 | for (ii = 1; ii < maxbin; ii++) { |
---|
| 1014 | vals[ii] = vals[ii] / (bins[ii] - bins[ii - 1]); |
---|
| 1015 | } |
---|
| 1016 | sum = vals[maxbin - 1]; |
---|
| 1017 | vals[0] = 0.; |
---|
| 1018 | } |
---|
| 1019 | for (ii = 0; ii < maxbin; ii++) { |
---|
| 1020 | vals[ii] = vals[ii] / sum; |
---|
| 1021 | IPDFEnergyH.InsertValues(bins[ii], vals[ii]); |
---|
| 1022 | } |
---|
| 1023 | |
---|
| 1024 | // Make IPDFEnergyExist = true |
---|
| 1025 | IPDFEnergyExist = true; |
---|
| 1026 | if (verbosityLevel > 1) |
---|
| 1027 | IPDFEnergyH.DumpValues(); |
---|
[816] | 1028 | } |
---|
[1347] | 1029 | |
---|
| 1030 | // IPDF has been create so carry on |
---|
| 1031 | G4double rndm = eneRndm->GenRandEnergy(); |
---|
| 1032 | particle_energy = IPDFEnergyH.GetEnergy(rndm); |
---|
| 1033 | |
---|
| 1034 | if (verbosityLevel >= 1) |
---|
| 1035 | G4cout << "Energy is " << particle_energy << G4endl; |
---|
| 1036 | } |
---|
| 1037 | |
---|
| 1038 | void G4SPSEneDistribution::GenArbPointEnergies() { |
---|
| 1039 | if (verbosityLevel > 0) |
---|
| 1040 | G4cout << "In GenArbPointEnergies" << G4endl; |
---|
| 1041 | G4double rndm; |
---|
| 1042 | rndm = eneRndm->GenRandEnergy(); |
---|
| 1043 | // IPDFArbEnergyH.DumpValues(); |
---|
| 1044 | // Find the Bin |
---|
| 1045 | // have x, y, no of points, and cumulative area distribution |
---|
| 1046 | G4int nabove, nbelow = 0, middle; |
---|
| 1047 | nabove = IPDFArbEnergyH.GetVectorLength(); |
---|
| 1048 | // G4cout << nabove << G4endl; |
---|
| 1049 | // Binary search to find bin that rndm is in |
---|
| 1050 | while (nabove - nbelow > 1) { |
---|
| 1051 | middle = (nabove + nbelow) / 2; |
---|
| 1052 | if (rndm == IPDFArbEnergyH(size_t(middle))) |
---|
| 1053 | break; |
---|
| 1054 | if (rndm < IPDFArbEnergyH(size_t(middle))) |
---|
| 1055 | nabove = middle; |
---|
| 1056 | else |
---|
| 1057 | nbelow = middle; |
---|
| 1058 | } |
---|
| 1059 | if (IntType == "Lin") { |
---|
| 1060 | Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1)); |
---|
[816] | 1061 | Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); |
---|
[1347] | 1062 | grad = Arb_grad[nbelow + 1]; |
---|
| 1063 | cept = Arb_cept[nbelow + 1]; |
---|
[816] | 1064 | // G4cout << rndm << " " << Emax << " " << Emin << " " << grad << " " << cept << G4endl; |
---|
| 1065 | GenerateLinearEnergies(true); |
---|
[1347] | 1066 | } else if (IntType == "Log") { |
---|
| 1067 | Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1)); |
---|
[816] | 1068 | Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); |
---|
[1347] | 1069 | alpha = Arb_alpha[nbelow + 1]; |
---|
[816] | 1070 | // G4cout << rndm << " " << Emax << " " << Emin << " " << alpha << G4endl; |
---|
| 1071 | GeneratePowEnergies(true); |
---|
[1347] | 1072 | } else if (IntType == "Exp") { |
---|
| 1073 | Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1)); |
---|
[816] | 1074 | Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); |
---|
[1347] | 1075 | Ezero = Arb_ezero[nbelow + 1]; |
---|
[816] | 1076 | // G4cout << rndm << " " << Emax << " " << Emin << " " << Ezero << G4endl; |
---|
| 1077 | GenerateExpEnergies(true); |
---|
[1347] | 1078 | } else if (IntType == "Spline") { |
---|
| 1079 | Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1)); |
---|
| 1080 | Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); |
---|
| 1081 | particle_energy = -1e100; |
---|
| 1082 | rndm = eneRndm->GenRandEnergy(); |
---|
| 1083 | while (particle_energy < Emin || particle_energy > Emax) { |
---|
| 1084 | particle_energy = SplineInt[nbelow+1]->CubicSplineInterpolation(rndm); |
---|
| 1085 | rndm = eneRndm->GenRandEnergy(); |
---|
| 1086 | } |
---|
| 1087 | if (verbosityLevel >= 1) |
---|
| 1088 | G4cout << "Energy is " << particle_energy << G4endl; |
---|
| 1089 | } else |
---|
| 1090 | G4cout << "Error: IntType unknown type" << G4endl; |
---|
[816] | 1091 | } |
---|
| 1092 | |
---|
[1347] | 1093 | void G4SPSEneDistribution::GenEpnHistEnergies() { |
---|
| 1094 | // G4cout << "In GenEpnHistEnergies " << Epnflag << G4endl; |
---|
[816] | 1095 | |
---|
[1347] | 1096 | // Firstly convert to energy if not already done. |
---|
| 1097 | if (Epnflag == true) |
---|
| 1098 | // epnflag = true means spectrum is epn, false means e. |
---|
[816] | 1099 | { |
---|
[1347] | 1100 | // convert to energy by multiplying by A number |
---|
| 1101 | ConvertEPNToEnergy(); |
---|
| 1102 | // EpnEnergyH will be replace by UDefEnergyH. |
---|
| 1103 | // UDefEnergyH.DumpValues(); |
---|
[816] | 1104 | } |
---|
[1347] | 1105 | |
---|
| 1106 | // G4cout << "Creating IPDFEnergy if not already done so" << G4endl; |
---|
| 1107 | if (IPDFEnergyExist == false) { |
---|
| 1108 | // IPDF has not been created, so create it |
---|
| 1109 | G4double bins[1024], vals[1024], sum; |
---|
| 1110 | G4int ii; |
---|
| 1111 | G4int maxbin = G4int(UDefEnergyH.GetVectorLength()); |
---|
| 1112 | bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0)); |
---|
| 1113 | vals[0] = UDefEnergyH(size_t(0)); |
---|
| 1114 | sum = vals[0]; |
---|
| 1115 | for (ii = 1; ii < maxbin; ii++) { |
---|
| 1116 | bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii)); |
---|
| 1117 | vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1]; |
---|
| 1118 | sum = sum + UDefEnergyH(size_t(ii)); |
---|
| 1119 | } |
---|
| 1120 | |
---|
| 1121 | for (ii = 0; ii < maxbin; ii++) { |
---|
| 1122 | vals[ii] = vals[ii] / sum; |
---|
| 1123 | IPDFEnergyH.InsertValues(bins[ii], vals[ii]); |
---|
| 1124 | } |
---|
| 1125 | // Make IPDFEpnExist = true |
---|
| 1126 | IPDFEnergyExist = true; |
---|
[816] | 1127 | } |
---|
[1347] | 1128 | // IPDFEnergyH.DumpValues(); |
---|
| 1129 | // IPDF has been create so carry on |
---|
| 1130 | G4double rndm = eneRndm->GenRandEnergy(); |
---|
| 1131 | particle_energy = IPDFEnergyH.GetEnergy(rndm); |
---|
[816] | 1132 | |
---|
[1347] | 1133 | if (verbosityLevel >= 1) |
---|
| 1134 | G4cout << "Energy is " << particle_energy << G4endl; |
---|
[816] | 1135 | } |
---|
| 1136 | |
---|
[1347] | 1137 | void G4SPSEneDistribution::ConvertEPNToEnergy() { |
---|
| 1138 | // Use this before particle generation to convert the |
---|
| 1139 | // currently stored histogram from energy/nucleon |
---|
| 1140 | // to energy. |
---|
| 1141 | // G4cout << "In ConvertEpntoEnergy " << G4endl; |
---|
| 1142 | if (particle_definition == NULL) |
---|
| 1143 | G4cout << "Error: particle not defined" << G4endl; |
---|
| 1144 | else { |
---|
| 1145 | // Need to multiply histogram by the number of nucleons. |
---|
| 1146 | // Baryon Number looks to hold the no. of nucleons. |
---|
| 1147 | G4int Bary = particle_definition->GetBaryonNumber(); |
---|
| 1148 | // G4cout << "Baryon No. " << Bary << G4endl; |
---|
| 1149 | // Change values in histogram, Read it out, delete it, re-create it |
---|
| 1150 | G4int count, maxcount; |
---|
| 1151 | maxcount = G4int(EpnEnergyH.GetVectorLength()); |
---|
| 1152 | // G4cout << maxcount << G4endl; |
---|
| 1153 | G4double ebins[1024], evals[1024]; |
---|
| 1154 | if (maxcount > 1024) { |
---|
| 1155 | G4cout << "Histogram contains more than 1024 bins!" << G4endl; |
---|
| 1156 | G4cout << "Those above 1024 will be ignored" << G4endl; |
---|
| 1157 | maxcount = 1024; |
---|
| 1158 | } |
---|
| 1159 | if (maxcount < 1) { |
---|
| 1160 | G4cout << "Histogram contains less than 1 bin!" << G4endl; |
---|
| 1161 | G4cout << "Redefine the histogram" << G4endl; |
---|
| 1162 | return; |
---|
| 1163 | } |
---|
| 1164 | for (count = 0; count < maxcount; count++) { |
---|
| 1165 | // Read out |
---|
| 1166 | ebins[count] = EpnEnergyH.GetLowEdgeEnergy(size_t(count)); |
---|
| 1167 | evals[count] = EpnEnergyH(size_t(count)); |
---|
| 1168 | } |
---|
[816] | 1169 | |
---|
[1347] | 1170 | // Multiply the channels by the nucleon number to give energies |
---|
| 1171 | for (count = 0; count < maxcount; count++) { |
---|
| 1172 | ebins[count] = ebins[count] * Bary; |
---|
| 1173 | } |
---|
[816] | 1174 | |
---|
[1347] | 1175 | // Set Emin and Emax |
---|
| 1176 | Emin = ebins[0]; |
---|
| 1177 | if (maxcount > 1) |
---|
| 1178 | Emax = ebins[maxcount - 1]; |
---|
| 1179 | else |
---|
| 1180 | Emax = ebins[0]; |
---|
| 1181 | // Put energy bins into new histogram - UDefEnergyH. |
---|
| 1182 | for (count = 0; count < maxcount; count++) { |
---|
| 1183 | UDefEnergyH.InsertValues(ebins[count], evals[count]); |
---|
| 1184 | } |
---|
| 1185 | Epnflag = false; //so that you dont repeat this method. |
---|
[816] | 1186 | } |
---|
| 1187 | } |
---|
| 1188 | |
---|
| 1189 | // |
---|
[1347] | 1190 | void G4SPSEneDistribution::ReSetHist(G4String atype) { |
---|
| 1191 | if (atype == "energy") { |
---|
| 1192 | UDefEnergyH = IPDFEnergyH = ZeroPhysVector; |
---|
| 1193 | IPDFEnergyExist = false; |
---|
| 1194 | Emin = 0.; |
---|
| 1195 | Emax = 1e30; |
---|
| 1196 | } else if (atype == "arb") { |
---|
| 1197 | ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector; |
---|
| 1198 | IPDFArbExist = false; |
---|
| 1199 | } else if (atype == "epn") { |
---|
| 1200 | UDefEnergyH = IPDFEnergyH = ZeroPhysVector; |
---|
| 1201 | IPDFEnergyExist = false; |
---|
| 1202 | EpnEnergyH = ZeroPhysVector; |
---|
| 1203 | } else { |
---|
| 1204 | G4cout << "Error, histtype not accepted " << G4endl; |
---|
| 1205 | } |
---|
[816] | 1206 | } |
---|
| 1207 | |
---|
[1347] | 1208 | G4double G4SPSEneDistribution::GenerateOne(G4ParticleDefinition* a) { |
---|
| 1209 | particle_definition = a; |
---|
| 1210 | particle_energy = -1.; |
---|
| 1211 | |
---|
| 1212 | while ((EnergyDisType == "Arb") ? (particle_energy < ArbEmin |
---|
| 1213 | || particle_energy > ArbEmax) : (particle_energy < Emin |
---|
| 1214 | || particle_energy > Emax)) { |
---|
| 1215 | if (Biased) { |
---|
| 1216 | GenerateBiasPowEnergies(); |
---|
| 1217 | } else { |
---|
| 1218 | if (EnergyDisType == "Mono") |
---|
| 1219 | GenerateMonoEnergetic(); |
---|
| 1220 | else if (EnergyDisType == "Lin") |
---|
| 1221 | GenerateLinearEnergies(); |
---|
| 1222 | else if (EnergyDisType == "Pow") |
---|
| 1223 | GeneratePowEnergies(); |
---|
| 1224 | else if (EnergyDisType == "Exp") |
---|
| 1225 | GenerateExpEnergies(); |
---|
| 1226 | else if (EnergyDisType == "Gauss") |
---|
| 1227 | GenerateGaussEnergies(); |
---|
| 1228 | else if (EnergyDisType == "Brem") |
---|
| 1229 | GenerateBremEnergies(); |
---|
| 1230 | else if (EnergyDisType == "Bbody") |
---|
| 1231 | GenerateBbodyEnergies(); |
---|
| 1232 | else if (EnergyDisType == "Cdg") |
---|
| 1233 | GenerateCdgEnergies(); |
---|
| 1234 | else if (EnergyDisType == "User") |
---|
| 1235 | GenUserHistEnergies(); |
---|
| 1236 | else if (EnergyDisType == "Arb") |
---|
| 1237 | GenArbPointEnergies(); |
---|
| 1238 | else if (EnergyDisType == "Epn") |
---|
| 1239 | GenEpnHistEnergies(); |
---|
| 1240 | else |
---|
| 1241 | G4cout << "Error: EnergyDisType has unusual value" << G4endl; |
---|
| 1242 | } |
---|
| 1243 | } |
---|
| 1244 | return particle_energy; |
---|
[816] | 1245 | } |
---|
| 1246 | |
---|
[1347] | 1247 | G4double G4SPSEneDistribution::GetProbability(G4double ene) { |
---|
| 1248 | G4double prob = 1.; |
---|
[816] | 1249 | |
---|
[1347] | 1250 | if (EnergyDisType == "Lin") { |
---|
| 1251 | if (prob_norm == 1.) { |
---|
| 1252 | prob_norm = 0.5*grad*Emax*Emax + cept*Emax - 0.5*grad*Emin*Emin - cept*Emin; |
---|
| 1253 | } |
---|
| 1254 | prob = cept + grad * ene; |
---|
| 1255 | prob /= prob_norm; |
---|
| 1256 | } |
---|
| 1257 | else if (EnergyDisType == "Pow") { |
---|
| 1258 | if (prob_norm == 1.) { |
---|
| 1259 | if (alpha != -1.) { |
---|
| 1260 | G4double emina = std::pow(Emin, alpha + 1); |
---|
| 1261 | G4double emaxa = std::pow(Emax, alpha + 1); |
---|
| 1262 | prob_norm = 1./(1.+alpha) * (emaxa - emina); |
---|
| 1263 | } else { |
---|
| 1264 | prob_norm = std::log(Emax) - std::log(Emin) ; |
---|
| 1265 | } |
---|
| 1266 | } |
---|
| 1267 | prob = std::pow(ene, alpha)/prob_norm; |
---|
| 1268 | } |
---|
| 1269 | else if (EnergyDisType == "Exp"){ |
---|
| 1270 | if (prob_norm == 1.) { |
---|
| 1271 | prob_norm = -Ezero*(std::exp(-Emax/Ezero) - std::exp(Emin/Ezero)); |
---|
| 1272 | } |
---|
| 1273 | prob = std::exp(-ene / Ezero); |
---|
| 1274 | prob /= prob_norm; |
---|
| 1275 | } |
---|
| 1276 | else if (EnergyDisType == "Arb") { |
---|
| 1277 | prob = ArbEnergyH.Value(ene); |
---|
| 1278 | // prob = ArbEInt->CubicSplineInterpolation(ene); |
---|
| 1279 | //G4double deltaY; |
---|
| 1280 | //prob = ArbEInt->PolynomInterpolation(ene, deltaY); |
---|
| 1281 | if (prob <= 0.) { |
---|
| 1282 | //G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << " " <<deltaY<< G4endl; |
---|
| 1283 | G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << G4endl; |
---|
| 1284 | prob = 1e-30; |
---|
| 1285 | } |
---|
| 1286 | // already normalised |
---|
| 1287 | } |
---|
| 1288 | else |
---|
| 1289 | G4cout << "Error: EnergyDisType not supported" << G4endl; |
---|
| 1290 | |
---|
| 1291 | return prob; |
---|
| 1292 | } |
---|