Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // G4SPSEneDistribution class implementation << 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 /////////////////////////////////////////////////////////////////////////////// 27 // 46 // 28 // Author: Fan Lei, QinetiQ ltd - 05/02/2004 << 29 // Customer: ESA/ESTEC << 30 // Revisions: Andrew Green, Andrea Dotti << 31 // ------------------------------------------- << 32 #include "G4SPSEneDistribution.hh" << 33 << 34 #include "G4Exp.hh" << 35 #include "G4SystemOfUnits.hh" << 36 #include "G4UnitsTable.hh" << 37 #include "Randomize.hh" 47 #include "Randomize.hh" 38 #include "G4AutoLock.hh" << 48 //#include <cmath> 39 #include "G4Threading.hh" << 49 >> 50 #include "G4SPSEneDistribution.hh" 40 51 41 G4SPSEneDistribution::G4SPSEneDistribution() << 52 G4SPSEneDistribution::G4SPSEneDistribution() { 42 { << 53 // 43 G4MUTEXINIT(mutex); << 54 // Initialise all variables 44 << 55 particle_energy = 1.0 * MeV; 45 // Initialise all variables << 56 46 << 57 EnergyDisType = "Mono"; 47 particle_energy = 1.0 * MeV; << 58 weight = 1.; 48 EnergyDisType = "Mono"; << 59 MonoEnergy = 1 * MeV; 49 weight=1.; << 60 Emin = 0.; 50 MonoEnergy = 1 * MeV; << 61 Emax = 1.e30; 51 Emin = 0.; << 62 alpha = 0.; 52 Emax = 1.e30; << 63 biasalpha = 0.; 53 alpha = 0.; << 64 prob_norm = 1.0; 54 biasalpha = 0.; << 65 Ezero = 0.; 55 prob_norm = 1.0; << 66 SE = 0.; 56 Ezero = 0.; << 67 Temp = 0.; 57 SE = 0.; << 68 grad = 0.; 58 Temp = 0.; << 69 cept = 0.; 59 grad = 0.; << 70 Biased = false; // not biased 60 cept = 0.; << 71 EnergySpec = true; // true - energy spectra, false - momentum spectra 61 IntType = "NULL"; // Interpolation type << 72 DiffSpec = true; // true - differential spec, false integral spec 62 << 73 IntType = "NULL"; // Interpolation type 63 ArbEmin = 0.; << 74 IPDFEnergyExist = false; 64 ArbEmax = 1.e30; << 75 IPDFArbExist = false; 65 << 76 66 verbosityLevel = 0; << 77 ArbEmin = 0.; 67 << 78 ArbEmax = 1.e30; 68 threadLocal_t& data = threadLocalData.Get(); << 79 69 data.Emax = Emax; << 80 verbosityLevel = 0; 70 data.Emin = Emin; << 81 71 data.alpha =alpha; << 82 } 72 data.cept = cept; << 83 73 data.Ezero = Ezero; << 84 G4SPSEneDistribution::~G4SPSEneDistribution() { 74 data.grad = grad; << 85 } 75 data.particle_energy = 0.; << 86 76 data.particle_definition = nullptr; << 87 void G4SPSEneDistribution::SetEnergyDisType(G4String DisType) { 77 data.weight = weight; << 88 EnergyDisType = DisType; 78 } << 89 if (EnergyDisType == "User") { 79 << 90 UDefEnergyH = IPDFEnergyH = ZeroPhysVector; 80 G4SPSEneDistribution::~G4SPSEneDistribution() << 91 IPDFEnergyExist = false; 81 { << 92 } else if (EnergyDisType == "Arb") { 82 G4MUTEXDESTROY(mutex); << 93 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector; 83 if(Arb_grad_cept_flag) << 94 IPDFArbExist = false; 84 { << 95 } else if (EnergyDisType == "Epn") { 85 delete [] Arb_grad; << 96 UDefEnergyH = IPDFEnergyH = ZeroPhysVector; 86 delete [] Arb_cept; << 97 IPDFEnergyExist = false; 87 } << 98 EpnEnergyH = ZeroPhysVector; 88 << 99 } 89 if(Arb_alpha_Const_flag) << 100 } 90 { << 101 91 delete [] Arb_alpha; << 102 void G4SPSEneDistribution::SetEmin(G4double emi) { 92 delete [] Arb_Const; << 103 Emin = emi; 93 } << 104 } 94 << 105 95 if(Arb_ezero_flag) << 106 void G4SPSEneDistribution::SetEmax(G4double ema) { 96 { << 107 Emax = ema; 97 delete [] Arb_ezero; << 108 } 98 } << 109 99 delete Bbody_x; << 110 void G4SPSEneDistribution::SetMonoEnergy(G4double menergy) { 100 delete BBHist; << 111 MonoEnergy = menergy; 101 delete CP_x; << 112 } 102 delete CPHist; << 113 103 for (auto & it : SplineInt) << 114 void G4SPSEneDistribution::SetBeamSigmaInE(G4double e) { 104 { << 115 SE = e; 105 delete it; << 116 } 106 it = nullptr; << 117 void G4SPSEneDistribution::SetAlpha(G4double alp) { 107 } << 118 alpha = alp; 108 SplineInt.clear(); << 119 } 109 } << 120 110 << 121 void G4SPSEneDistribution::SetBiasAlpha(G4double alp) { 111 void G4SPSEneDistribution::SetEnergyDisType(co << 122 biasalpha = alp; 112 { << 123 Biased = true; 113 G4AutoLock l(&mutex); << 124 } 114 EnergyDisType = DisType; << 125 115 if (EnergyDisType == "User") << 126 void G4SPSEneDistribution::SetTemp(G4double tem) { 116 { << 127 Temp = tem; 117 UDefEnergyH = IPDFEnergyH = ZeroPhysVector << 128 } 118 IPDFEnergyExist = false; << 129 119 } << 130 void G4SPSEneDistribution::SetEzero(G4double eze) { 120 else if (EnergyDisType == "Arb") << 131 Ezero = eze; 121 { << 132 } 122 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVect << 133 123 IPDFArbExist = false; << 134 void G4SPSEneDistribution::SetGradient(G4double gr) { 124 } << 135 grad = gr; 125 else if (EnergyDisType == "Epn") << 136 } 126 { << 137 127 UDefEnergyH = IPDFEnergyH = ZeroPhysVector << 138 void G4SPSEneDistribution::SetInterCept(G4double c) { 128 IPDFEnergyExist = false; << 139 cept = c; 129 EpnEnergyH = ZeroPhysVector; << 140 } 130 } << 141 131 } << 142 void G4SPSEneDistribution::UserEnergyHisto(G4ThreeVector input) { 132 << 143 G4double ehi, val; 133 const G4String& G4SPSEneDistribution::GetEnerg << 144 ehi = input.x(); 134 { << 145 val = input.y(); 135 G4AutoLock l(&mutex); << 146 if (verbosityLevel > 1) { 136 return EnergyDisType; << 147 G4cout << "In UserEnergyHisto" << G4endl; 137 } << 148 G4cout << " " << ehi << " " << val << G4endl; 138 << 149 } 139 void G4SPSEneDistribution::SetEmin(G4double em << 150 UDefEnergyH.InsertValues(ehi, val); 140 { << 151 Emax = ehi; 141 G4AutoLock l(&mutex); << 152 } 142 Emin = emi; << 153 143 threadLocalData.Get().Emin = Emin; << 154 void G4SPSEneDistribution::ArbEnergyHisto(G4ThreeVector input) { 144 } << 155 G4double ehi, val; 145 << 156 ehi = input.x(); 146 G4double G4SPSEneDistribution::GetEmin() const << 157 val = input.y(); 147 { << 158 if (verbosityLevel > 1) { 148 return threadLocalData.Get().Emin; << 159 G4cout << "In ArbEnergyHisto" << G4endl; 149 } << 160 G4cout << " " << ehi << " " << val << G4endl; 150 << 161 } 151 G4double G4SPSEneDistribution::GetArbEmin() << 162 ArbEnergyH.InsertValues(ehi, val); 152 { << 163 } 153 G4AutoLock l(&mutex); << 164 154 return ArbEmin; << 165 void G4SPSEneDistribution::ArbEnergyHistoFile(G4String filename) { 155 } << 166 std::ifstream infile(filename, std::ios::in); 156 << 167 if (!infile) 157 G4double G4SPSEneDistribution::GetArbEmax() << 168 G4Exception("Unable to open the histo ASCII file"); 158 { << 169 G4double ehi, val; 159 G4AutoLock l(&mutex); << 170 while (infile >> ehi >> val) { 160 return ArbEmax; << 171 ArbEnergyH.InsertValues(ehi, val); 161 } << 172 } 162 << 173 } 163 void G4SPSEneDistribution::SetEmax(G4double em << 174 164 { << 175 void G4SPSEneDistribution::EpnEnergyHisto(G4ThreeVector input) { 165 G4AutoLock l(&mutex); << 176 G4double ehi, val; 166 Emax = ema; << 177 ehi = input.x(); 167 threadLocalData.Get().Emax = Emax; << 178 val = input.y(); 168 } << 179 if (verbosityLevel > 1) { 169 << 180 G4cout << "In EpnEnergyHisto" << G4endl; 170 G4double G4SPSEneDistribution::GetEmax() const << 181 G4cout << " " << ehi << " " << val << G4endl; 171 { << 182 } 172 return threadLocalData.Get().Emax; << 183 EpnEnergyH.InsertValues(ehi, val); 173 } << 184 Emax = ehi; 174 << 185 Epnflag = true; 175 void G4SPSEneDistribution::SetMonoEnergy(G4dou << 186 } 176 { << 187 177 G4AutoLock l(&mutex); << 188 void G4SPSEneDistribution::Calculate() { 178 MonoEnergy = menergy; << 189 if (EnergyDisType == "Cdg") 179 } << 190 CalculateCdgSpectrum(); 180 << 191 else if (EnergyDisType == "Bbody") 181 void G4SPSEneDistribution::SetBeamSigmaInE(G4d << 192 CalculateBbodySpectrum(); 182 { << 193 } 183 G4AutoLock l(&mutex); << 194 184 SE = e; << 195 void G4SPSEneDistribution::CalculateCdgSpectrum() { 185 } << 196 // This uses the spectrum from The INTEGRAL Mass Model (TIMM) 186 void G4SPSEneDistribution::SetAlpha(G4double a << 197 // to generate a Cosmic Diffuse X/gamma ray spectrum. 187 { << 198 G4double pfact[2] = { 8.5, 112 }; 188 G4AutoLock l(&mutex); << 199 G4double spind[2] = { 1.4, 2.3 }; 189 alpha = alp; << 200 G4double ene_line[3] = { 1. * keV, 18. * keV, 1E6 * keV }; 190 threadLocalData.Get().alpha = alpha; << 201 G4int n_par; 191 } << 202 192 << 203 ene_line[0] = Emin; 193 void G4SPSEneDistribution::SetBiasAlpha(G4doub << 204 if (Emin < 18 * keV) { 194 { << 205 n_par = 2; 195 G4AutoLock l(&mutex); << 206 ene_line[2] = Emax; 196 biasalpha = alp; << 207 if (Emax < 18 * keV) { 197 Biased = true; << 208 n_par = 1; 198 } << 209 ene_line[1] = Emax; 199 << 210 } 200 void G4SPSEneDistribution::SetTemp(G4double te << 211 } else { 201 { << 212 n_par = 1; 202 G4AutoLock l(&mutex); << 213 pfact[0] = 112.; 203 Temp = tem; << 214 spind[0] = 2.3; 204 } << 215 ene_line[1] = Emax; 205 << 216 } 206 void G4SPSEneDistribution::SetEzero(G4double e << 217 207 { << 218 // Create a cumulative histogram. 208 G4AutoLock l(&mutex); << 219 CDGhist[0] = 0.; 209 Ezero = eze; << 220 G4double omalpha; 210 threadLocalData.Get().Ezero = Ezero; << 221 G4int i = 0; 211 } << 222 212 << 223 while (i < n_par) { 213 void G4SPSEneDistribution::SetGradient(G4doubl << 224 omalpha = 1. - spind[i]; 214 { << 225 CDGhist[i + 1] = CDGhist[i] + (pfact[i] / omalpha) * (std::pow( 215 G4AutoLock l(&mutex); << 226 ene_line[i + 1] / keV, omalpha) - std::pow(ene_line[i] / keV, 216 grad = gr; << 227 omalpha)); 217 threadLocalData.Get().grad = grad; << 228 i++; 218 } << 229 } 219 << 230 220 void G4SPSEneDistribution::SetInterCept(G4doub << 231 // Normalise histo and divide by 1000 to make MeV. 221 { << 232 i = 0; 222 G4AutoLock l(&mutex); << 233 while (i < n_par) { 223 cept = c; << 234 CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[n_par]; 224 threadLocalData.Get().cept = cept; << 235 // G4cout << CDGhist[i] << CDGhist[n_par] << G4endl; 225 } << 236 i++; 226 << 237 } 227 const G4String& G4SPSEneDistribution::GetIntTy << 238 } 228 { << 239 229 G4AutoLock l(&mutex); << 240 void G4SPSEneDistribution::CalculateBbodySpectrum() { 230 return IntType; << 241 // create bbody spectrum 231 } << 242 // Proved very hard to integrate indefinitely, so different 232 << 243 // method. User inputs emin, emax and T. These are used to 233 void G4SPSEneDistribution::SetBiasRndm(G4SPSRa << 244 // create a 10,000 bin histogram. 234 { << 245 // Use photon density spectrum = 2 nu**2/c**2 * (std::exp(h nu/kT)-1) 235 G4AutoLock l(&mutex); << 246 // = 2 E**2/h**2c**2 times the exponential 236 eneRndm = a; << 247 G4double erange = Emax - Emin; 237 } << 248 G4double steps = erange / 10000.; 238 << 249 G4double Bbody_y[10000]; 239 void G4SPSEneDistribution::SetVerbosity(G4int << 250 G4double k = 8.6181e-11; //Boltzmann const in MeV/K 240 { << 251 G4double h = 4.1362e-21; // Plancks const in MeV s 241 G4AutoLock l(&mutex); << 252 G4double c = 3e8; // Speed of light 242 verbosityLevel = a; << 253 G4double h2 = h * h; 243 } << 254 G4double c2 = c * c; 244 << 255 G4int count = 0; 245 G4double G4SPSEneDistribution::GetWeight() con << 256 G4double sum = 0.; 246 { << 257 BBHist[0] = 0.; 247 return threadLocalData.Get().weight; << 258 while (count < 10000) { 248 } << 259 Bbody_x[count] = Emin + G4double(count * steps); 249 << 260 Bbody_y[count] = (2. * std::pow(Bbody_x[count], 2.)) / (h2 * c2 250 G4double G4SPSEneDistribution::GetMonoEnergy() << 261 * (std::exp(Bbody_x[count] / (k * Temp)) - 1.)); 251 { << 262 sum = sum + Bbody_y[count]; 252 G4AutoLock l(&mutex); << 263 BBHist[count + 1] = BBHist[count] + Bbody_y[count]; 253 return MonoEnergy; << 264 count++; 254 } << 265 } 255 << 266 256 G4double G4SPSEneDistribution::GetSE() << 267 Bbody_x[10000] = Emax; 257 { << 268 // Normalise cumulative histo. 258 G4AutoLock l(&mutex); << 269 count = 0; 259 return SE; << 270 while (count < 10001) { 260 } << 271 BBHist[count] = BBHist[count] / sum; 261 << 272 count++; 262 G4double G4SPSEneDistribution::Getalpha() cons << 273 } 263 { << 274 } 264 return threadLocalData.Get().alpha; << 275 265 } << 276 void G4SPSEneDistribution::InputEnergySpectra(G4bool value) { 266 << 277 // Allows user to specifiy spectrum is momentum 267 G4double G4SPSEneDistribution::GetEzero() cons << 278 EnergySpec = value; // false if momentum 268 { << 279 if (verbosityLevel > 1) 269 return threadLocalData.Get().Ezero; << 280 G4cout << "EnergySpec has value " << EnergySpec << G4endl; 270 } << 281 } 271 << 282 272 G4double G4SPSEneDistribution::GetTemp() << 283 void G4SPSEneDistribution::InputDifferentialSpectra(G4bool value) { 273 { << 284 // Allows user to specify integral or differential spectra 274 G4AutoLock l(&mutex); << 285 DiffSpec = value; // true = differential, false = integral 275 return Temp; << 286 if (verbosityLevel > 1) 276 } << 287 G4cout << "Diffspec has value " << DiffSpec << G4endl; 277 << 288 } 278 G4double G4SPSEneDistribution::Getgrad() const << 289 279 { << 290 void G4SPSEneDistribution::ArbInterpolate(G4String IType) { 280 return threadLocalData.Get().grad; << 291 if (EnergyDisType != "Arb") 281 } << 292 G4cout << "Error: this is for arbitrary distributions" << G4endl; 282 << 293 IntType = IType; 283 G4double G4SPSEneDistribution::Getcept() const << 294 ArbEmax = ArbEnergyH.GetMaxLowEdgeEnergy(); 284 { << 295 ArbEmin = ArbEnergyH.GetMinLowEdgeEnergy(); 285 return threadLocalData.Get().cept; << 296 286 } << 297 // Now interpolate points 287 << 298 if (IntType == "Lin") 288 G4PhysicsFreeVector G4SPSEneDistribution::GetU << 299 LinearInterpolation(); 289 { << 300 if (IntType == "Log") 290 G4AutoLock l(&mutex); << 301 LogInterpolation(); 291 return UDefEnergyH; << 302 if (IntType == "Exp") 292 } << 303 ExpInterpolation(); 293 << 304 if (IntType == "Spline") 294 G4PhysicsFreeVector G4SPSEneDistribution::GetA << 305 SplineInterpolation(); 295 { << 306 } 296 G4AutoLock l(&mutex); << 307 297 return ArbEnergyH; << 308 void G4SPSEneDistribution::LinearInterpolation() { 298 } << 309 // Method to do linear interpolation on the Arb points 299 << 310 // Calculate equation of each line segment, max 1024. 300 void G4SPSEneDistribution::UserEnergyHisto(con << 311 // Calculate Area under each segment 301 { << 312 // Create a cumulative array which is then normalised Arb_Cum_Area 302 G4AutoLock l(&mutex); << 313 303 G4double ehi = input.x(), << 314 G4double Area_seg[1024]; // Stores area under each segment 304 val = input.y(); << 315 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024]; 305 if (verbosityLevel > 1) << 316 G4int i, count; 306 { << 317 G4int maxi = ArbEnergyH.GetVectorLength(); 307 G4cout << "In UserEnergyHisto" << G4endl; << 318 for (i = 0; i < maxi; i++) { 308 G4cout << " " << ehi << " " << val << G4en << 319 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); 309 } << 320 Arb_y[i] = ArbEnergyH(size_t(i)); 310 UDefEnergyH.InsertValues(ehi, val); << 321 } 311 Emax = ehi; << 322 // Points are now in x,y arrays. If the spectrum is integral it has to be 312 threadLocalData.Get().Emax = Emax; << 323 // made differential and if momentum it has to be made energy. 313 } << 324 if (DiffSpec == false) { 314 << 325 // Converts integral point-wise spectra to Differential 315 void G4SPSEneDistribution::ArbEnergyHisto(cons << 326 for (count = 0; count < maxi - 1; count++) { 316 { << 327 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) 317 G4AutoLock l(&mutex); << 328 / (Arb_x[count + 1] - Arb_x[count]); 318 G4double ehi = input.x(), << 329 } 319 val = input.y(); << 330 maxi--; 320 if (verbosityLevel > 1) << 331 } 321 { << 332 // 322 G4cout << "In ArbEnergyHisto" << G4endl; << 333 if (EnergySpec == false) { 323 G4cout << " " << ehi << " " << val << G4en << 334 // change currently stored values (emin etc) which are actually momenta 324 } << 335 // to energies. 325 ArbEnergyH.InsertValues(ehi, val); << 336 if (particle_definition == NULL) 326 } << 337 G4cout << "Error: particle not defined" << G4endl; 327 << 338 else { 328 void G4SPSEneDistribution::ArbEnergyHistoFile( << 339 // Apply Energy**2 = p**2c**2 + m0**2c**4 329 { << 340 // p should be entered as E/c i.e. without the division by c 330 G4AutoLock l(&mutex); << 341 // being done - energy equivalent. 331 std::ifstream infile(filename, std::ios::in) << 342 G4double mass = particle_definition->GetPDGMass(); 332 if (!infile) << 343 // convert point to energy unit and its value to per energy unit 333 { << 344 G4double total_energy; 334 G4Exception("G4SPSEneDistribution::ArbEner << 345 for (count = 0; count < maxi; count++) { 335 FatalException, "Unable to ope << 346 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass 336 } << 347 * mass)); // total energy 337 G4double ehi, val; << 348 338 while (infile >> ehi >> val) << 349 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy; 339 { << 350 Arb_x[count] = total_energy - mass; // kinetic energy 340 ArbEnergyH.InsertValues(ehi, val); << 351 } 341 } << 352 } 342 } << 353 } 343 << 354 // 344 void G4SPSEneDistribution::EpnEnergyHisto(cons << 355 i = 1; 345 { << 356 Arb_grad[0] = 0.; 346 G4AutoLock l(&mutex); << 357 Arb_cept[0] = 0.; 347 G4double ehi = input.x(), << 358 Area_seg[0] = 0.; 348 val = input.y(); << 359 Arb_Cum_Area[0] = 0.; 349 if (verbosityLevel > 1) << 360 while (i < maxi) { 350 { << 361 // calc gradient and intercept for each segment 351 G4cout << "In EpnEnergyHisto" << G4endl; << 362 Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[i - 1]); 352 G4cout << " " << ehi << " " << val << G4en << 363 if (verbosityLevel == 2) 353 } << 364 G4cout << Arb_grad[i] << G4endl; 354 EpnEnergyH.InsertValues(ehi, val); << 365 if (Arb_grad[i] > 0.) { 355 Emax = ehi; << 366 if (verbosityLevel == 2) 356 threadLocalData.Get().Emax = Emax; << 367 G4cout << "Arb_grad is positive" << G4endl; 357 Epnflag = true; << 368 Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]); 358 } << 369 } else if (Arb_grad[i] < 0.) { 359 << 370 if (verbosityLevel == 2) 360 void G4SPSEneDistribution::Calculate() << 371 G4cout << "Arb_grad is negative" << G4endl; 361 { << 372 Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]); 362 G4AutoLock l(&mutex); << 373 } else { 363 if (EnergyDisType == "Cdg") << 374 if (verbosityLevel == 2) 364 { << 375 G4cout << "Arb_grad is 0." << G4endl; 365 CalculateCdgSpectrum(); << 376 Arb_cept[i] = Arb_y[i]; 366 } << 377 } 367 else if (EnergyDisType == "Bbody") << 378 368 { << 379 Area_seg[i] = ((Arb_grad[i] / 2) * (Arb_x[i] * Arb_x[i] - Arb_x[i - 1] 369 if(!BBhistInit) << 380 * Arb_x[i - 1]) + Arb_cept[i] * (Arb_x[i] - Arb_x[i - 1])); 370 { << 381 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i]; 371 BBInitHists(); << 382 sum = sum + Area_seg[i]; 372 } << 383 if (verbosityLevel == 2) 373 CalculateBbodySpectrum(); << 384 G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i] 374 } << 385 << G4endl; 375 else if (EnergyDisType == "CPow") << 386 i++; 376 { << 387 } 377 if(!CPhistInit) << 388 378 { << 389 i = 0; 379 CPInitHists(); << 390 while (i < maxi) { 380 } << 391 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation 381 CalculateCPowSpectrum(); << 392 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); 382 } << 393 i++; 383 } << 394 } 384 << 395 385 void G4SPSEneDistribution::BBInitHists() // M << 396 // now scale the ArbEnergyH, needed by Probability() 386 { << 397 ArbEnergyH.ScaleVector(1., 1./sum); 387 BBHist = new std::vector<G4double>(10001, 0. << 398 388 Bbody_x = new std::vector<G4double>(10001, 0 << 399 if (verbosityLevel >= 1) { 389 BBhistInit = true; << 400 G4cout << "Leaving LinearInterpolation" << G4endl; 390 } << 401 ArbEnergyH.DumpValues(); 391 << 402 IPDFArbEnergyH.DumpValues(); 392 void G4SPSEneDistribution::CPInitHists() // M << 403 } 393 { << 404 } 394 CPHist = new std::vector<G4double>(10001, 0. << 405 395 CP_x = new std::vector<G4double>(10001, 0.0) << 406 void G4SPSEneDistribution::LogInterpolation() { 396 CPhistInit = true; << 407 // Interpolation based on Logarithmic equations 397 } << 408 // Generate equations of line segments 398 << 409 // y = Ax**alpha => log y = alpha*logx + logA 399 void G4SPSEneDistribution::CalculateCdgSpectru << 410 // Find area under line segments 400 { << 411 // create normalised, cumulative array Arb_Cum_Area 401 // This uses the spectrum from the INTEGRAL << 412 G4double Area_seg[1024]; // Stores area under each segment 402 // to generate a Cosmic Diffuse X/gamma ray << 413 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024]; 403 << 414 G4int i, count; 404 G4double pfact[2] = { 8.5, 112 }; << 415 G4int maxi = ArbEnergyH.GetVectorLength(); 405 G4double spind[2] = { 1.4, 2.3 }; << 416 for (i = 0; i < maxi; i++) { 406 G4double ene_line[3] = { 1. * keV, 18. * keV << 417 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); 407 G4int n_par; << 418 Arb_y[i] = ArbEnergyH(size_t(i)); 408 << 419 } 409 ene_line[0] = threadLocalData.Get().Emin; << 420 // Points are now in x,y arrays. If the spectrum is integral it has to be 410 if (threadLocalData.Get().Emin < 18 * keV) << 421 // made differential and if momentum it has to be made energy. 411 { << 422 if (DiffSpec == false) { 412 n_par = 2; << 423 // Converts integral point-wise spectra to Differential 413 ene_line[2] = threadLocalData.Get().Emax; << 424 for (count = 0; count < maxi - 1; count++) { 414 if (threadLocalData.Get().Emax < 18 * keV) << 425 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) 415 { << 426 / (Arb_x[count + 1] - Arb_x[count]); 416 n_par = 1; << 427 } 417 ene_line[1] = threadLocalData.Get().Emax << 428 maxi--; 418 } << 429 } 419 } << 430 // 420 else << 431 if (EnergySpec == false) { 421 { << 432 // change currently stored values (emin etc) which are actually momenta 422 n_par = 1; << 433 // to energies. 423 pfact[0] = 112.; << 434 if (particle_definition == NULL) 424 spind[0] = 2.3; << 435 G4cout << "Error: particle not defined" << G4endl; 425 ene_line[1] = threadLocalData.Get().Emax; << 436 else { 426 } << 437 // Apply Energy**2 = p**2c**2 + m0**2c**4 427 << 438 // p should be entered as E/c i.e. without the division by c 428 // Create a cumulative histogram << 439 // being done - energy equivalent. 429 // << 440 G4double mass = particle_definition->GetPDGMass(); 430 CDGhist[0] = 0.; << 441 // convert point to energy unit and its value to per energy unit 431 G4double omalpha; << 442 G4double total_energy; 432 G4int i = 0; << 443 for (count = 0; count < maxi; count++) { 433 while (i < n_par) << 444 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass 434 { << 445 * mass)); // total energy 435 omalpha = 1. - spind[i]; << 446 436 CDGhist[i + 1] = CDGhist[i] + (pfact[i] / << 447 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy; 437 * (std::pow(en << 448 Arb_x[count] = total_energy - mass; // kinetic energy 438 - std::pow(ene << 449 } 439 ++i; << 450 } 440 } << 451 } 441 << 452 // 442 // Normalise histo and divide by 1000 to mak << 453 i = 1; 443 // << 454 Arb_alpha[0] = 0.; 444 i = 0; << 455 Arb_Const[0] = 0.; 445 while (i < n_par) << 456 Area_seg[0] = 0.; 446 { << 457 Arb_Cum_Area[0] = 0.; 447 CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[ << 458 if (Arb_x[0] <= 0. || Arb_y[0] <= 0.) { 448 ++i; << 459 G4cout << "You should not use log interpolation with points <= 0." 449 } << 460 << G4endl; 450 } << 461 G4cout << "These will be changed to 1e-20, which may cause problems" 451 << 462 << G4endl; 452 void G4SPSEneDistribution::CalculateBbodySpect << 463 if (Arb_x[0] <= 0.) 453 { << 464 Arb_x[0] = 1e-20; 454 // Create bbody spectrum << 465 if (Arb_y[0] <= 0.) 455 // Proved very hard to integrate indefinitel << 466 Arb_y[0] = 1e-20; 456 // User inputs emin, emax and T. These are u << 467 } 457 // bin histogram. << 468 458 // Use photon density spectrum = 2 nu**2/c** << 469 G4double alp; 459 // = 2 E**2/h**2c**2 times the exponential << 470 while (i < maxi) { 460 << 471 // Incase points are negative or zero 461 G4double erange = threadLocalData.Get().Emax << 472 if (Arb_x[i] <= 0. || Arb_y[i] <= 0.) { 462 G4double steps = erange / 10000.; << 473 G4cout << "You should not use log interpolation with points <= 0." 463 << 474 << G4endl; 464 const G4double k = 8.6181e-11; //Boltzmann c << 475 G4cout 465 const G4double h = 4.1362e-21; // Plancks co << 476 << "These will be changed to 1e-20, which may cause problems" 466 const G4double c = 3e8; // Speed of light << 477 << G4endl; 467 const G4double h2 = h * h; << 478 if (Arb_x[i] <= 0.) 468 const G4double c2 = c * c; << 479 Arb_x[i] = 1e-20; 469 G4int count = 0; << 480 if (Arb_y[i] <= 0.) 470 G4double sum = 0.; << 481 Arb_y[i] = 1e-20; 471 BBHist->at(0) = 0.; << 482 } 472 << 483 473 while (count < 10000) << 484 Arb_alpha[i] = (std::log10(Arb_y[i]) - std::log10(Arb_y[i - 1])) 474 { << 485 / (std::log10(Arb_x[i]) - std::log10(Arb_x[i - 1])); 475 Bbody_x->at(count) = threadLocalData.Get() << 486 Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[i], Arb_alpha[i])); 476 G4double Bbody_y = (2. * std::pow(Bbody_x- << 487 alp = Arb_alpha[i] + 1; 477 / (h2*c2*(std::exp(Bbody_ << 488 if (alp == 0.) { 478 sum = sum + Bbody_y; << 489 Area_seg[i] = Arb_Const[i] * (std::log(Arb_x[i]) - std::log(Arb_x[i - 1])); 479 BBHist->at(count + 1) = BBHist->at(count) << 490 } else { 480 ++count; << 491 Area_seg[i] = (Arb_Const[i] / alp) * (std::pow(Arb_x[i], alp) 481 } << 492 - std::pow(Arb_x[i - 1], alp)); 482 << 493 } 483 Bbody_x->at(10000) = threadLocalData.Get().E << 494 sum = sum + Area_seg[i]; 484 << 495 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i]; 485 // Normalise cumulative histo << 496 if (verbosityLevel == 2) 486 // << 497 G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl; 487 count = 0; << 498 i++; 488 while (count < 10001) << 499 } 489 { << 500 490 BBHist->at(count) = BBHist->at(count) / su << 501 i = 0; 491 ++count; << 502 while (i < maxi) { 492 } << 503 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; 493 } << 504 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); 494 << 505 i++; 495 void G4SPSEneDistribution::CalculateCPowSpectr << 506 } 496 { << 507 497 // Create cutoff power-law spectrum, x^a exp << 508 // now scale the ArbEnergyH, needed by Probability() 498 // The integral of this function is an incom << 509 ArbEnergyH.ScaleVector(1., 1./sum); 499 // is only available in the Boost library. << 510 500 // << 511 if (verbosityLevel >= 1) 501 // User inputs are emin, emax and alpha and << 512 G4cout << "Leaving LogInterpolation " << G4endl; 502 // create a 10,000 bin histogram. << 513 } 503 << 514 504 G4double erange = threadLocalData.Get().Emax << 515 void G4SPSEneDistribution::ExpInterpolation() { 505 G4double steps = erange / 10000.; << 516 // Interpolation based on Exponential equations 506 alpha = threadLocalData.Get().alpha ; << 517 // Generate equations of line segments 507 Ezero = threadLocalData.Get().Ezero ; << 518 // y = Ae**-(x/e0) => ln y = -x/e0 + lnA 508 << 519 // Find area under line segments 509 G4int count = 0; << 520 // create normalised, cumulative array Arb_Cum_Area 510 G4double sum = 0.; << 521 G4double Area_seg[1024]; // Stores area under each segment 511 CPHist->at(0) = 0.; << 522 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024]; 512 << 523 G4int i, count; 513 while (count < 10000) << 524 G4int maxi = ArbEnergyH.GetVectorLength(); 514 { << 525 for (i = 0; i < maxi; i++) { 515 CP_x->at(count) = threadLocalData.Get().Em << 526 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); 516 G4double CP_y = std::pow(CP_x->at(count), << 527 Arb_y[i] = ArbEnergyH(size_t(i)); 517 * std::exp(-CP_x->at(count) << 528 } 518 sum = sum + CP_y; << 529 // Points are now in x,y arrays. If the spectrum is integral it has to be 519 CPHist->at(count + 1) = CPHist->at(count) << 530 // made differential and if momentum it has to be made energy. 520 ++count; << 531 if (DiffSpec == false) { 521 } << 532 // Converts integral point-wise spectra to Differential 522 << 533 for (count = 0; count < maxi - 1; count++) { 523 CP_x->at(10000) = threadLocalData.Get().Emax << 534 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) 524 << 535 / (Arb_x[count + 1] - Arb_x[count]); 525 // Normalise cumulative histo << 536 } 526 // << 537 maxi--; 527 count = 0; << 538 } 528 while (count < 10001) << 539 // 529 { << 540 if (EnergySpec == false) { 530 CPHist->at(count) = CPHist->at(count) / su << 541 // change currently stored values (emin etc) which are actually momenta 531 ++count; << 542 // to energies. 532 } << 543 if (particle_definition == NULL) 533 } << 544 G4cout << "Error: particle not defined" << G4endl; 534 << 545 else { 535 void G4SPSEneDistribution::InputEnergySpectra( << 546 // Apply Energy**2 = p**2c**2 + m0**2c**4 536 { << 547 // p should be entered as E/c i.e. without the division by c 537 G4AutoLock l(&mutex); << 548 // being done - energy equivalent. 538 << 549 G4double mass = particle_definition->GetPDGMass(); 539 // Allows user to specify spectrum is moment << 550 // convert point to energy unit and its value to per energy unit 540 // << 551 G4double total_energy; 541 EnergySpec = value; // false if momentum << 552 for (count = 0; count < maxi; count++) { 542 if (verbosityLevel > 1) << 553 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass 543 { << 554 * mass)); // total energy 544 G4cout << "EnergySpec has value " << Energ << 555 545 } << 556 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy; 546 } << 557 Arb_x[count] = total_energy - mass; // kinetic energy 547 << 558 } 548 void G4SPSEneDistribution::InputDifferentialSp << 559 } 549 { << 560 } 550 G4AutoLock l(&mutex); << 561 // 551 << 562 i = 1; 552 // Allows user to specify integral or differ << 563 Arb_ezero[0] = 0.; 553 // << 564 Arb_Const[0] = 0.; 554 DiffSpec = value; // true = differential, fa << 565 Area_seg[0] = 0.; 555 if (verbosityLevel > 1) << 566 Arb_Cum_Area[0] = 0.; 556 { << 567 while (i < maxi) { 557 G4cout << "Diffspec has value " << DiffSpe << 568 G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i - 1]); 558 } << 569 if (test > 0. || test < 0.) { 559 } << 570 Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1]) / (std::log(Arb_y[i]) 560 << 571 - std::log(Arb_y[i - 1])); 561 void G4SPSEneDistribution::ArbInterpolate(cons << 572 Arb_Const[i] = Arb_y[i] / (std::exp(-Arb_x[i] / Arb_ezero[i])); 562 { << 573 Area_seg[i] = -(Arb_Const[i] * Arb_ezero[i]) * (std::exp(-Arb_x[i] 563 G4AutoLock l(&mutex); << 574 / Arb_ezero[i]) - std::exp(-Arb_x[i - 1] / Arb_ezero[i])); 564 << 575 } else { 565 IntType = IType; << 576 G4cout << "Flat line segment: problem" << G4endl; 566 ArbEmax = ArbEnergyH.GetMaxEnergy(); << 577 Arb_ezero[i] = 0.; 567 ArbEmin = ArbEnergyH.Energy(0); << 578 Arb_Const[i] = 0.; 568 << 579 Area_seg[i] = 0.; 569 // Now interpolate points << 580 } 570 << 581 sum = sum + Area_seg[i]; 571 if (IntType == "Lin") LinearInterpolation(); << 582 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i]; 572 if (IntType == "Log") LogInterpolation(); << 583 if (verbosityLevel == 2) 573 if (IntType == "Exp") ExpInterpolation(); << 584 G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl; 574 if (IntType == "Spline") SplineInterpolation << 585 i++; 575 } << 586 } 576 << 587 577 void G4SPSEneDistribution::LinearInterpolation << 588 i = 0; 578 { << 589 while (i < maxi) { 579 // Method to do linear interpolation on the << 590 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; 580 // Calculate equation of each line segment, << 591 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); 581 // Calculate Area under each segment << 592 i++; 582 // Create a cumulative array which is then n << 593 } 583 << 594 584 G4double Area_seg[1024]; // Stores area unde << 595 // now scale the ArbEnergyH, needed by Probability() 585 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 596 ArbEnergyH.ScaleVector(1., 1./sum); 586 std::size_t i, count; << 597 587 std::size_t maxi = ArbEnergyH.GetVectorLengt << 598 if (verbosityLevel >= 1) 588 for (i = 0; i < maxi; ++i) << 599 G4cout << "Leaving ExpInterpolation " << G4endl; 589 { << 600 } 590 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 601 591 Arb_y[i] = ArbEnergyH(i); << 602 void G4SPSEneDistribution::SplineInterpolation() { 592 } << 603 // Interpolation using Splines. 593 << 604 // Create Normalised arrays, make x 0->1 and y hold 594 // Points are now in x,y arrays. If the spec << 605 // the function (Energy) 595 // made differential and if momentum it has << 606 // 596 << 607 // Current method based on the above will not work in all cases. 597 if (!DiffSpec) << 608 // new method is implemented below. 598 { << 599 // Converts integral point-wise spectra to << 600 // << 601 for (count = 0; count < maxi-1; ++count) << 602 { << 603 Arb_y[count] = (Arb_y[count] - Arb_y[cou << 604 / (Arb_x[count + 1] - Arb_x << 605 } << 606 --maxi; << 607 } << 608 << 609 if (!EnergySpec) << 610 { << 611 // change currently stored values (emin et << 612 // to energies << 613 // << 614 G4ParticleDefinition* pdef = threadLocalDa << 615 if (pdef == nullptr) << 616 { << 617 G4Exception("G4SPSEneDistribution::Linea << 618 "Event0302", FatalException, << 619 "Error: particle not defined << 620 } << 621 else << 622 { << 623 // Apply Energy**2 = p**2c**2 + m0**2c** << 624 // p should be entered as E/c i.e. witho << 625 // being done - energy equivalent << 626 << 627 G4double mass = pdef->GetPDGMass(); << 628 << 629 // Convert point to energy unit and its << 630 // << 631 G4double total_energy; << 632 for (count = 0; count < maxi; ++count) << 633 { << 634 total_energy = std::sqrt((Arb_x[count] << 635 + (mass * mass)); // tota << 636 Arb_y[count] = Arb_y[count] * Arb_x[co << 637 Arb_x[count] = total_energy - mass; // << 638 } << 639 } << 640 } << 641 << 642 i = 1; << 643 << 644 Arb_grad = new G4double [1024]; << 645 Arb_cept = new G4double [1024]; << 646 Arb_grad_cept_flag = true; << 647 << 648 Arb_grad[0] = 0.; << 649 Arb_cept[0] = 0.; << 650 Area_seg[0] = 0.; << 651 Arb_Cum_Area[0] = 0.; << 652 while (i < maxi) << 653 { << 654 // calculate gradient and intercept for ea << 655 // << 656 Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / << 657 if (verbosityLevel == 2) << 658 { << 659 G4cout << Arb_grad[i] << G4endl; << 660 } << 661 if (Arb_grad[i] > 0.) << 662 { << 663 if (verbosityLevel == 2) << 664 { << 665 G4cout << "Arb_grad is positive" << G << 666 } << 667 Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * << 668 } << 669 else if (Arb_grad[i] < 0.) << 670 { << 671 if (verbosityLevel == 2) << 672 { << 673 G4cout << "Arb_grad is negative" << G4 << 674 } << 675 Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * << 676 } << 677 else << 678 { << 679 if (verbosityLevel == 2) << 680 { << 681 G4cout << "Arb_grad is 0." << G4endl; << 682 } << 683 Arb_cept[i] = Arb_y[i]; << 684 } << 685 << 686 Area_seg[i] = ((Arb_grad[i] / 2) << 687 * (Arb_x[i] * Arb_x[i] - Arb_x << 688 + Arb_cept[i] * (Arb_x[i] - Ar << 689 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Ar << 690 sum = sum + Area_seg[i]; << 691 if (verbosityLevel == 2) << 692 { << 693 G4cout << Arb_x[i] << Arb_y[i] << Area_s << 694 << Arb_grad[i] << G4endl; << 695 } << 696 ++i; << 697 } << 698 << 699 i = 0; << 700 while (i < maxi) << 701 { << 702 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; / << 703 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_ << 704 ++i; << 705 } << 706 << 707 // now scale the ArbEnergyH, needed by Proba << 708 // << 709 ArbEnergyH.ScaleVector(1., 1./sum); << 710 << 711 if (verbosityLevel >= 1) << 712 { << 713 G4cout << "Leaving LinearInterpolation" << << 714 ArbEnergyH.DumpValues(); << 715 IPDFArbEnergyH.DumpValues(); << 716 } << 717 } << 718 << 719 void G4SPSEneDistribution::LogInterpolation() << 720 { << 721 // Interpolation based on Logarithmic equati << 722 // Generate equations of line segments << 723 // y = Ax**alpha => log y = alpha*logx + log << 724 // Find area under line segments << 725 // Create normalised, cumulative array Arb_C << 726 << 727 G4double Area_seg[1024]; // Stores area unde << 728 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 729 std::size_t i, count; << 730 std::size_t maxi = ArbEnergyH.GetVectorLengt << 731 for (i = 0; i < maxi; ++i) << 732 { << 733 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 734 Arb_y[i] = ArbEnergyH(i); << 735 } << 736 << 737 // Points are now in x,y arrays. If the spec << 738 // made differential and if momentum it has << 739 << 740 if (!DiffSpec) << 741 { << 742 // Converts integral point-wise spectra to << 743 // << 744 for (count = 0; count < maxi-1; ++count) << 745 { << 746 Arb_y[count] = (Arb_y[count] - Arb_y[cou << 747 / (Arb_x[count + 1] - Arb_x << 748 } << 749 --maxi; << 750 } << 751 << 752 if (!EnergySpec) << 753 { << 754 // Change currently stored values (emin et << 755 // to energies << 756 << 757 G4ParticleDefinition* pdef = threadLocalDa << 758 if (pdef == nullptr) << 759 { << 760 G4Exception("G4SPSEneDistribution::LogIn << 761 "Event0302", FatalException, << 762 "Error: particle not defined << 763 } << 764 else << 765 { << 766 // Apply Energy**2 = p**2c**2 + m0**2c** << 767 // p should be entered as E/c i.e. witho << 768 // being done - energy equivalent << 769 << 770 G4double mass = pdef->GetPDGMass(); << 771 << 772 // Convert point to energy unit and its << 773 // << 774 G4double total_energy; << 775 for (count = 0; count < maxi; ++count) << 776 { << 777 total_energy = std::sqrt((Arb_x[count] << 778 Arb_y[count] = Arb_y[count] * Arb_x[co << 779 Arb_x[count] = total_energy - mass; // << 780 } << 781 } << 782 } << 783 << 784 i = 1; << 785 << 786 if ( Arb_ezero != nullptr ) { delete [] Arb_ << 787 if ( Arb_Const != nullptr ) { delete [] Arb_ << 788 Arb_alpha = new G4double [1024]; << 789 Arb_Const = new G4double [1024]; << 790 Arb_alpha_Const_flag = true; << 791 << 792 Arb_alpha[0] = 0.; << 793 Arb_Const[0] = 0.; << 794 Area_seg[0] = 0.; << 795 Arb_Cum_Area[0] = 0.; << 796 if (Arb_x[0] <= 0. || Arb_y[0] <= 0.) << 797 { << 798 G4cout << "You should not use log interpol << 799 << G4endl; << 800 G4cout << "These will be changed to 1e-20, << 801 << G4endl; << 802 if (Arb_x[0] <= 0.) << 803 { << 804 Arb_x[0] = 1e-20; << 805 } << 806 if (Arb_y[0] <= 0.) << 807 { << 808 Arb_y[0] = 1e-20; << 809 } << 810 } << 811 << 812 G4double alp; << 813 while (i < maxi) << 814 { << 815 // In case points are negative or zero << 816 // << 817 if (Arb_x[i] <= 0. || Arb_y[i] <= 0.) << 818 { << 819 G4cout << "You should not use log interp << 820 << G4endl; << 821 G4cout << "These will be changed to 1e-2 << 822 << G4endl; << 823 if (Arb_x[i] <= 0.) << 824 { << 825 Arb_x[i] = 1e-20; << 826 } << 827 if (Arb_y[i] <= 0.) << 828 { << 829 Arb_y[i] = 1e-20; << 830 } << 831 } << 832 << 833 Arb_alpha[i] = (std::log10(Arb_y[i]) - std << 834 / (std::log10(Arb_x[i]) - std << 835 Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[ << 836 alp = Arb_alpha[i] + 1; << 837 if (alp == 0.) << 838 { << 839 Area_seg[i] = Arb_Const[i] << 840 * (std::log(Arb_x[i]) - std: << 841 } << 842 else << 843 { << 844 Area_seg[i] = (Arb_Const[i] / alp) << 845 * (std::pow(Arb_x[i], alp) - << 846 } << 847 sum = sum + Area_seg[i]; << 848 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Ar << 849 if (verbosityLevel == 2) << 850 { << 851 G4cout << Arb_alpha[i] << Arb_Const[i] < << 852 } << 853 ++i; << 854 } << 855 << 856 i = 0; << 857 while (i < maxi) << 858 { << 859 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; << 860 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_ << 861 ++i; << 862 } << 863 << 864 // Now scale the ArbEnergyH, needed by Proba << 865 // << 866 ArbEnergyH.ScaleVector(1., 1./sum); << 867 << 868 if (verbosityLevel >= 1) << 869 { << 870 G4cout << "Leaving LogInterpolation " << G << 871 } << 872 } << 873 << 874 void G4SPSEneDistribution::ExpInterpolation() << 875 { << 876 // Interpolation based on Exponential equati << 877 // Generate equations of line segments << 878 // y = Ae**-(x/e0) => ln y = -x/e0 + lnA << 879 // Find area under line segments << 880 // Create normalised, cumulative array Arb_C << 881 << 882 G4double Area_seg[1024]; // Stores area unde << 883 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 884 std::size_t i, count; << 885 std::size_t maxi = ArbEnergyH.GetVectorLengt << 886 for (i = 0; i < maxi; ++i) << 887 { << 888 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 889 Arb_y[i] = ArbEnergyH(i); << 890 } << 891 << 892 // Points are now in x,y arrays. If the spec << 893 // made differential and if momentum it has << 894 << 895 if (!DiffSpec) << 896 { << 897 // Converts integral point-wise spectra to << 898 // << 899 for (count = 0; count < maxi - 1; ++count) << 900 { << 901 Arb_y[count] = (Arb_y[count] - Arb_y[cou << 902 / (Arb_x[count + 1] - Arb_x << 903 } << 904 --maxi; << 905 } << 906 << 907 if (!EnergySpec) << 908 { << 909 // Change currently stored values (emin et << 910 // to energies << 911 // << 912 G4ParticleDefinition* pdef = threadLocalDa << 913 if (pdef == nullptr) << 914 { << 915 G4Exception("G4SPSEneDistribution::ExpIn << 916 "Event0302", FatalException, << 917 "Error: particle not defined << 918 } << 919 else << 920 { << 921 // Apply Energy**2 = p**2c**2 + m0**2c** << 922 // p should be entered as E/c i.e. witho << 923 // being done - energy equivalent << 924 << 925 G4double mass = pdef->GetPDGMass(); << 926 << 927 // Convert point to energy unit and its << 928 // << 929 G4double total_energy; << 930 for (count = 0; count < maxi; ++count) << 931 { << 932 total_energy = std::sqrt((Arb_x[count] << 933 Arb_y[count] = Arb_y[count] * Arb_x[co << 934 Arb_x[count] = total_energy - mass; // << 935 } << 936 } << 937 } << 938 << 939 i = 1; << 940 << 941 if ( Arb_ezero != nullptr ) { delete[] Arb_e << 942 if ( Arb_Const != nullptr ) { delete[] Arb_C << 943 Arb_ezero = new G4double [1024]; << 944 Arb_Const = new G4double [1024]; << 945 Arb_ezero_flag = true; << 946 << 947 Arb_ezero[0] = 0.; << 948 Arb_Const[0] = 0.; << 949 Area_seg[0] = 0.; << 950 Arb_Cum_Area[0] = 0.; << 951 while (i < maxi) << 952 { << 953 G4double test = std::log(Arb_y[i]) - std:: << 954 if (test > 0. || test < 0.) << 955 { << 956 Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1] << 957 / (std::log(Arb_y[i]) - std << 958 Arb_Const[i] = Arb_y[i] / (std::exp(-Arb << 959 Area_seg[i] = -(Arb_Const[i] * Arb_ezero << 960 * (std::exp(-Arb_x[i] / Arb_ << 961 - std::exp(-Arb_x[i - 1] / << 962 } << 963 else << 964 { << 965 G4Exception("G4SPSEneDistribution::ExpIn << 966 "Event0302", JustWarning, << 967 "Flat line segment: problem, << 968 G4cout << "Flat line segment: problem" < << 969 Arb_ezero[i] = 0.; << 970 Arb_Const[i] = 0.; << 971 Area_seg[i] = 0.; << 972 } << 973 sum = sum + Area_seg[i]; << 974 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Ar << 975 if (verbosityLevel == 2) << 976 { << 977 G4cout << Arb_ezero[i] << Arb_Const[i] < << 978 } << 979 ++i; << 980 } << 981 << 982 i = 0; << 983 while (i < maxi) << 984 { << 985 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; << 986 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_ << 987 ++i; << 988 } << 989 << 990 // Now scale the ArbEnergyH, needed by Proba << 991 // << 992 ArbEnergyH.ScaleVector(1., 1./sum); << 993 << 994 if (verbosityLevel >= 1) << 995 { << 996 G4cout << "Leaving ExpInterpolation " << G << 997 } << 998 } << 999 << 1000 void G4SPSEneDistribution::SplineInterpolatio << 1001 { << 1002 // Interpolation using Splines. << 1003 // Create Normalised arrays, make x 0->1 an << 1004 // << 1005 // Current method based on the above will n << 1006 // New method is implemented below. << 1007 609 1008 G4double sum, Arb_x[1024]={0.}, Arb_y[1024] << 610 G4double sum, Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024]; 1009 std::size_t i, count; << 611 G4int i, count; 1010 std::size_t maxi = ArbEnergyH.GetVectorLeng << 1011 << 1012 for (i = 0; i < maxi; ++i) << 1013 { << 1014 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i) << 1015 Arb_y[i] = ArbEnergyH(i); << 1016 } << 1017 << 1018 // Points are now in x,y arrays. If the spe << 1019 // made differential and if momentum it has << 1020 << 1021 if (!DiffSpec) << 1022 { << 1023 // Converts integral point-wise spectra t << 1024 // << 1025 for (count = 0; count < maxi - 1; ++count << 1026 { << 1027 Arb_y[count] = (Arb_y[count] - Arb_y[co << 1028 / (Arb_x[count + 1] - Arb_ << 1029 } << 1030 --maxi; << 1031 } << 1032 << 1033 if (!EnergySpec) << 1034 { << 1035 // Change currently stored values (emin e << 1036 // to energies << 1037 // << 1038 G4ParticleDefinition* pdef = threadLocalD << 1039 if (pdef == nullptr) << 1040 { << 1041 G4Exception("G4SPSEneDistribution::Spli << 1042 "Event0302", FatalException << 1043 "Error: particle not define << 1044 } << 1045 else << 1046 { << 1047 // Apply Energy**2 = p**2c**2 + m0**2c* << 1048 // p should be entered as E/c i.e. with << 1049 // being done - energy equivalent << 1050 << 1051 G4double mass = pdef->GetPDGMass(); << 1052 << 1053 // Convert point to energy unit and its << 1054 // << 1055 G4double total_energy; << 1056 for (count = 0; count < maxi; ++count) << 1057 { << 1058 total_energy = std::sqrt((Arb_x[count << 1059 Arb_y[count] = Arb_y[count] * Arb_x[c << 1060 Arb_x[count] = total_energy - mass; / << 1061 } << 1062 } << 1063 } << 1064 << 1065 i = 1; << 1066 Arb_Cum_Area[0] = 0.; << 1067 sum = 0.; << 1068 Splinetemp = new G4DataInterpolation(Arb_x, << 1069 G4double ei[101], prob[101]; << 1070 for (auto & it : SplineInt) << 1071 { << 1072 delete it; << 1073 it = nullptr; << 1074 } << 1075 SplineInt.clear(); << 1076 SplineInt.resize(1024,nullptr); << 1077 while (i < maxi) << 1078 { << 1079 // 100 step per segment for the integrati << 1080 << 1081 G4double de = (Arb_x[i] - Arb_x[i - 1])/1 << 1082 G4double area = 0.; << 1083 << 1084 for (count = 0; count < 101; ++count) << 1085 { << 1086 ei[count] = Arb_x[i - 1] + de*count ; << 1087 prob[count] = Splinetemp->CubicSplineI << 1088 if (prob[count] < 0.) << 1089 { << 1090 G4ExceptionDescription ED; << 1091 ED << "Warning: G4DataInterpolation r << 1092 << " " << ei[count] << G4endl; << 1093 G4Exception("G4SPSEneDistribution::Sp << 1094 FatalException, ED); << 1095 } << 1096 area += prob[count]*de; << 1097 } << 1098 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + a << 1099 sum += area; << 1100 << 1101 prob[0] = prob[0]/(area/de); << 1102 for (count = 1; count < 100; ++count) << 1103 { << 1104 prob[count] = prob[count-1] + prob[coun << 1105 } << 1106 << 1107 SplineInt[i] = new G4DataInterpolation(pr << 1108 << 1109 // NOTE: i starts from 1! << 1110 // << 1111 ++i; << 1112 } << 1113 << 1114 i = 0; << 1115 while (i < maxi) << 1116 { << 1117 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; << 1118 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb << 1119 ++i; << 1120 } << 1121 << 1122 // Now scale the ArbEnergyH, needed by Prob << 1123 // << 1124 ArbEnergyH.ScaleVector(1., 1./sum); << 1125 << 1126 if (verbosityLevel > 0) << 1127 { << 1128 G4cout << "Leaving SplineInterpolation " << 1129 } << 1130 } << 1131 << 1132 void G4SPSEneDistribution::GenerateMonoEnerge << 1133 { << 1134 // Method to generate MonoEnergetic particl << 1135 << 1136 threadLocalData.Get().particle_energy = Mon << 1137 } << 1138 << 1139 void G4SPSEneDistribution::GenerateGaussEnerg << 1140 { << 1141 // Method to generate Gaussian particles << 1142 << 1143 G4double ene = G4RandGauss::shoot(MonoEnerg << 1144 if (ene < 0) ene = 0.; << 1145 threadLocalData.Get().particle_energy = ene << 1146 } << 1147 << 1148 void G4SPSEneDistribution::GenerateLinearEner << 1149 { << 1150 G4double rndm; << 1151 threadLocal_t& params = threadLocalData.Get << 1152 G4double emaxsq = std::pow(params.Emax, 2.) << 1153 G4double eminsq = std::pow(params.Emin, 2.) << 1154 G4double intersq = std::pow(params.cept, 2. << 1155 << 1156 if (bArb) rndm = G4UniformRand(); << 1157 else rndm = eneRndm->GenRandEnergy(); << 1158 << 1159 G4double bracket = ((params.grad / 2.) << 1160 * (emaxsq - eminsq) << 1161 + params.cept * (params.Em << 1162 bracket = bracket * rndm; << 1163 bracket = bracket + (params.grad / 2.) * em << 1164 << 1165 // Now have a quad of form m/2 E**2 + cE - << 1166 // << 1167 bracket = -bracket; << 1168 << 1169 if (params.grad != 0.) << 1170 { << 1171 G4double sqbrack = (intersq - 4 * (params << 1172 sqbrack = std::sqrt(sqbrack); << 1173 G4double root1 = -params.cept + sqbrack; << 1174 root1 = root1 / (2. * (params.grad / 2.)) << 1175 << 1176 G4double root2 = -params.cept - sqbrack; << 1177 root2 = root2 / (2. * (params.grad / 2.)) << 1178 << 1179 if (root1 > params.Emin && root1 < params << 1180 { << 1181 params.particle_energy = root1; << 1182 } << 1183 if (root2 > params.Emin && root2 < params << 1184 { << 1185 params.particle_energy = root2; << 1186 } << 1187 } << 1188 else if (params.grad == 0.) << 1189 { << 1190 // have equation of form cE - bracket =0 << 1191 // << 1192 params.particle_energy = bracket / params << 1193 } << 1194 << 1195 if (params.particle_energy < 0.) << 1196 { << 1197 params.particle_energy = -params.particle << 1198 } << 1199 << 1200 if (verbosityLevel >= 1) << 1201 { << 1202 G4cout << "Energy is " << params.particle << 1203 } << 1204 } << 1205 << 1206 void G4SPSEneDistribution::GeneratePowEnergie << 1207 { << 1208 // Method to generate particle energies dis << 1209 << 1210 G4double rndm; << 1211 G4double emina, emaxa; << 1212 << 1213 threadLocal_t& params = threadLocalData.Get << 1214 << 1215 emina = std::pow(params.Emin, params.alpha << 1216 emaxa = std::pow(params.Emax, params.alpha << 1217 << 1218 if (bArb) rndm = G4UniformRand(); << 1219 else rndm = eneRndm->GenRandEnergy(); << 1220 << 1221 if (params.alpha != -1.) << 1222 { << 1223 G4double ene = ((rndm * (emaxa - emina)) << 1224 ene = std::pow(ene, (1. / (params.alpha + << 1225 params.particle_energy = ene; << 1226 } << 1227 else << 1228 { << 1229 G4double ene = (std::log(params.Emin) << 1230 + rndm * (std::log(params.Em << 1231 params.particle_energy = std::exp(ene); << 1232 } << 1233 if (verbosityLevel >= 1) << 1234 { << 1235 G4cout << "Energy is " << params.particle << 1236 } << 1237 } << 1238 << 1239 void G4SPSEneDistribution::GenerateCPowEnergi << 1240 { << 1241 // Method to generate particle energies dis << 1242 // cutoff power-law distribution << 1243 // << 1244 // CP_x holds Energies, and CPHist holds th << 1245 // binary search to find correct bin then l << 1246 // Use the earlier defined histogram + Rand << 1247 // random numbers following the histos dist << 1248 << 1249 G4double rndm = eneRndm->GenRandEnergy(); << 1250 G4int nabove = 10001, nbelow = 0, middle; << 1251 << 1252 G4AutoLock l(&mutex); << 1253 G4bool done = CPhistCalcd; << 1254 l.unlock(); << 1255 << 1256 if(!done) << 1257 { << 1258 Calculate(); //This is has a lock inside, << 1259 l.lock(); << 1260 CPhistCalcd = true; << 1261 l.unlock(); << 1262 } << 1263 << 1264 // Binary search to find bin that rndm is i << 1265 // << 1266 while (nabove - nbelow > 1) << 1267 { << 1268 middle = (nabove + nbelow) / 2; << 1269 if (rndm == CPHist->at(middle)) << 1270 { << 1271 break; << 1272 } << 1273 if (rndm < CPHist->at(middle)) << 1274 { << 1275 nabove = middle; << 1276 } << 1277 else << 1278 { << 1279 nbelow = middle; << 1280 } << 1281 } << 1282 << 1283 // Now interpolate in that bin to find the << 1284 // << 1285 G4double x1, x2, y1, y2, t, q; << 1286 x1 = CP_x->at(nbelow); << 1287 if(nbelow+1 == static_cast<G4int>(CP_x->siz << 1288 { << 1289 x2 = CP_x->back(); << 1290 } << 1291 else << 1292 { << 1293 x2 = CP_x->at(nbelow + 1); << 1294 } << 1295 y1 = CPHist->at(nbelow); << 1296 if(nbelow+1 == static_cast<G4int>(CPHist->s << 1297 { << 1298 G4cout << CPHist->back() << G4endl; << 1299 y2 = CPHist->back(); << 1300 } << 1301 else << 1302 { << 1303 y2 = CPHist->at(nbelow + 1); << 1304 } << 1305 t = (y2 - y1) / (x2 - x1); << 1306 q = y1 - t * x1; << 1307 << 1308 threadLocalData.Get().particle_energy = (rn << 1309 << 1310 if (verbosityLevel >= 1) << 1311 { << 1312 G4cout << "Energy is " << threadLocalData << 1313 } << 1314 } << 1315 << 1316 void G4SPSEneDistribution::GenerateBiasPowEne << 1317 { << 1318 // Method to generate particle energies dis << 1319 // in biased power-law and calculate its we << 1320 << 1321 threadLocal_t& params = threadLocalData.Get << 1322 << 1323 G4double rndm; << 1324 G4double emina, emaxa, emin, emax; << 1325 << 1326 G4double normal = 1.; << 1327 << 1328 emin = params.Emin; << 1329 emax = params.Emax; << 1330 << 1331 rndm = eneRndm->GenRandEnergy(); << 1332 << 1333 if (biasalpha != -1.) << 1334 { << 1335 emina = std::pow(emin, biasalpha + 1); << 1336 emaxa = std::pow(emax, biasalpha + 1); << 1337 G4double ee = ((rndm * (emaxa - emina)) + << 1338 params.particle_energy = std::pow(ee, (1. << 1339 normal = 1./(1+biasalpha) * (emaxa - emin << 1340 } << 1341 else << 1342 { << 1343 G4double ee = (std::log(emin) + rndm * (s << 1344 params.particle_energy = std::exp(ee); << 1345 normal = std::log(emax) - std::log(emin); << 1346 } << 1347 params.weight = GetProbability(params.parti << 1348 / (std::pow(params.particle_e << 1349 << 1350 if (verbosityLevel >= 1) << 1351 { << 1352 G4cout << "Energy is " << params.particle << 1353 } << 1354 } << 1355 << 1356 void G4SPSEneDistribution::GenerateExpEnergie << 1357 { << 1358 // Method to generate particle energies dis << 1359 // to an exponential curve << 1360 << 1361 G4double rndm; << 1362 << 1363 if (bArb) rndm = G4UniformRand(); << 1364 else rndm = eneRndm->GenRandEnergy(); << 1365 << 1366 threadLocal_t& params = threadLocalData.Get << 1367 params.particle_energy = -params.Ezero << 1368 * (std::log(rndm * ( << 1369 << 1370 - << 1371 << 1372 + std::exp << 1373 if (verbosityLevel >= 1) << 1374 { << 1375 G4cout << "Energy is " << params.particle << 1376 } << 1377 } << 1378 << 1379 void G4SPSEneDistribution::GenerateBremEnergi << 1380 { << 1381 // Method to generate particle energies dis << 1382 // to a Bremstrahlung equation of the form << 1383 // I = const*((kT)**1/2)*E*(e**(-E/kT)) << 1384 << 1385 G4double rndm = eneRndm->GenRandEnergy(); << 1386 G4double expmax, expmin, k; << 1387 << 1388 k = 8.6181e-11; // Boltzmann's const in MeV << 1389 G4double ksq = std::pow(k, 2.); // k square << 1390 G4double Tsq = std::pow(Temp, 2.); // Temp << 1391 << 1392 threadLocal_t& params = threadLocalData.Get << 1393 << 1394 expmax = std::exp(-params.Emax / (k * Temp) << 1395 expmin = std::exp(-params.Emin / (k * Temp) << 1396 << 1397 // If either expmax or expmin are zero then << 1398 // Most probably this will be because T is << 1399 << 1400 if (expmax == 0.) << 1401 { << 1402 G4Exception("G4SPSEneDistribution::Genera << 1403 "Event0302", FatalException, << 1404 "*****EXPMAX=0. Choose differ << 1405 } << 1406 if (expmin == 0.) << 1407 { << 1408 G4Exception("G4SPSEneDistribution::Genera << 1409 "Event0302", FatalException, << 1410 "*****EXPMIN=0. Choose differ << 1411 } << 1412 << 1413 G4double tempvar = rndm * ((-k) * Temp * (p << 1414 - p << 1415 - (ksq * Tsq * (exp << 1416 << 1417 G4double bigc = (tempvar - k * Temp * param << 1418 - ksq * Tsq * expmin) / (-k << 1419 << 1420 // This gives an equation of form: Ee(-E/kT << 1421 // Solve this iteratively, step from Emin t << 1422 // and take the best solution. << 1423 << 1424 G4double erange = params.Emax - params.Emin << 1425 G4double steps = erange / 1000.; << 1426 G4int i; << 1427 G4double etest, diff, err = 100000.; << 1428 << 1429 for (i = 1; i < 1000; ++i) << 1430 { << 1431 etest = params.Emin + (i - 1) * steps; << 1432 diff = etest * (std::exp(-etest / (k * Te << 1433 + k * Temp * (std::exp(-etest / (k * << 1434 << 1435 if (diff < 0.) << 1436 { << 1437 diff = -diff; << 1438 } << 1439 << 1440 if (diff < err) << 1441 { << 1442 err = diff; << 1443 params.particle_energy = etest; << 1444 } << 1445 } << 1446 if (verbosityLevel >= 1) << 1447 { << 1448 G4cout << "Energy is " << params.particle << 1449 } << 1450 } << 1451 << 1452 void G4SPSEneDistribution::GenerateBbodyEnerg << 1453 { << 1454 // BBody_x holds Energies, and BBHist holds << 1455 // Binary search to find correct bin then l << 1456 // Use the earlier defined histogram + Rand << 1457 // random numbers following the histos dist << 1458 << 1459 G4double rndm = eneRndm->GenRandEnergy(); << 1460 G4int nabove = 10001, nbelow = 0, middle; << 1461 << 1462 G4AutoLock l(&mutex); << 1463 G4bool done = BBhistCalcd; << 1464 l.unlock(); << 1465 << 1466 if(!done) << 1467 { << 1468 Calculate(); //This is has a lock inside, << 1469 l.lock(); << 1470 BBhistCalcd = true; << 1471 l.unlock(); << 1472 } << 1473 << 1474 // Binary search to find bin that rndm is i << 1475 // << 1476 while (nabove - nbelow > 1) << 1477 { << 1478 middle = (nabove + nbelow) / 2; << 1479 if (rndm == BBHist->at(middle)) << 1480 { << 1481 break; << 1482 } << 1483 if (rndm < BBHist->at(middle)) << 1484 { << 1485 nabove = middle; << 1486 } << 1487 else << 1488 { << 1489 nbelow = middle; << 1490 } << 1491 } << 1492 << 1493 // Now interpolate in that bin to find the << 1494 // << 1495 G4double x1, x2, y1, y2, t, q; << 1496 x1 = Bbody_x->at(nbelow); << 1497 << 1498 if(nbelow+1 == static_cast<G4int>(Bbody_x-> << 1499 { << 1500 x2 = Bbody_x->back(); << 1501 } << 1502 else << 1503 { << 1504 x2 = Bbody_x->at(nbelow + 1); << 1505 } << 1506 y1 = BBHist->at(nbelow); << 1507 if(nbelow+1 == static_cast<G4int>(BBHist->s << 1508 { << 1509 G4cout << BBHist->back() << G4endl; << 1510 y2 = BBHist->back(); << 1511 } << 1512 else << 1513 { << 1514 y2 = BBHist->at(nbelow + 1); << 1515 } << 1516 t = (y2 - y1) / (x2 - x1); << 1517 q = y1 - t * x1; << 1518 << 1519 threadLocalData.Get().particle_energy = (rn << 1520 << 1521 if (verbosityLevel >= 1) << 1522 { << 1523 G4cout << "Energy is " << threadLocalData << 1524 } << 1525 } << 1526 << 1527 void G4SPSEneDistribution::GenerateCdgEnergie << 1528 { << 1529 // Generate random numbers, compare with va << 1530 // to find appropriate part of spectrum and << 1531 // generate energy in the usual inversion w << 1532 << 1533 G4double rndm, rndm2; << 1534 G4double ene_line[3]={0,0,0}; << 1535 G4double omalpha[2]={0,0}; << 1536 threadLocal_t& params = threadLocalData.Get << 1537 if (params.Emin < 18 * keV && params.Emax < << 1538 { << 1539 omalpha[0] = 1. - 1.4; << 1540 ene_line[0] = params.Emin; << 1541 ene_line[1] = params.Emax; << 1542 } << 1543 if (params.Emin < 18 * keV && params.Emax > << 1544 { << 1545 omalpha[0] = 1. - 1.4; << 1546 omalpha[1] = 1. - 2.3; << 1547 ene_line[0] = params.Emin; << 1548 ene_line[1] = 18. * keV; << 1549 ene_line[2] = params.Emax; << 1550 } << 1551 if (params.Emin > 18 * keV) << 1552 { << 1553 omalpha[0] = 1. - 2.3; << 1554 ene_line[0] = params.Emin; << 1555 ene_line[1] = params.Emax; << 1556 } << 1557 rndm = eneRndm->GenRandEnergy(); << 1558 rndm2 = eneRndm->GenRandEnergy(); << 1559 << 1560 G4int i = 0; << 1561 while (rndm >= CDGhist[i]) << 1562 { << 1563 ++i; << 1564 } << 1565 << 1566 // Generate final energy << 1567 // << 1568 G4double ene = (std::pow(ene_line[i - 1], o << 1569 + (std::pow(ene_line[i], omalp << 1570 - std::pow(ene_line[i - 1], o << 1571 params.particle_energy = std::pow(ene, (1. << 1572 << 1573 if (verbosityLevel >= 1) << 1574 { << 1575 G4cout << "Energy is " << params.particle << 1576 } << 1577 } << 1578 << 1579 void G4SPSEneDistribution::GenUserHistEnergie << 1580 { << 1581 // Histograms are DIFFERENTIAL << 1582 << 1583 G4AutoLock l(&mutex); << 1584 << 1585 if (!IPDFEnergyExist) << 1586 { << 1587 std::size_t ii; << 1588 std::size_t maxbin = UDefEnergyH.GetVecto << 1589 G4double bins[1024], vals[1024], sum; << 1590 for ( ii = 0 ; ii<1024 ; ++ii ) { bins[ii << 1591 sum = 0.; << 1592 << 1593 if ( (!EnergySpec) << 1594 && (threadLocalData.Get().particle_defi << 1595 { << 1596 G4Exception("G4SPSEneDistribution::GenU << 1597 "Event0302", FatalException << 1598 "Error: particle definition << 1599 } << 1600 << 1601 if (maxbin > 1024) << 1602 { << 1603 G4Exception("G4SPSEneDistribution::GenU << 1604 "Event0302", JustWarning, << 1605 "Maxbin>1024\n Setting maxbi << 1606 maxbin = 1024; << 1607 } << 1608 << 1609 if (!DiffSpec) << 1610 { << 1611 G4cout << "Histograms are Differential! << 1612 } << 1613 else << 1614 { << 1615 bins[0] = UDefEnergyH.GetLowEdgeEnergy( << 1616 vals[0] = UDefEnergyH(0); << 1617 sum = vals[0]; << 1618 for (ii = 1; ii < maxbin; ++ii) << 1619 { << 1620 bins[ii] = UDefEnergyH.GetLowEdgeEner << 1621 vals[ii] = UDefEnergyH(ii) + vals[ii << 1622 sum = sum + UDefEnergyH(ii); << 1623 } << 1624 } << 1625 << 1626 if (!EnergySpec) << 1627 { << 1628 G4double mass = threadLocalData.Get().p << 1629 << 1630 // Multiply the function (vals) up by t << 1631 // to make the function counts/s (i.e. << 1632 << 1633 for (ii = 1; ii < maxbin; ++ii) << 1634 { << 1635 vals[ii] = vals[ii] * (bins[ii] - bin << 1636 } << 1637 << 1638 // Put energy bins into new histo, plus << 1639 // to make evals counts/s/energy << 1640 // << 1641 for (ii = 0; ii < maxbin; ++ii) << 1642 { << 1643 // kinetic energy << 1644 // << 1645 bins[ii] = std::sqrt((bins[ii]*bins[i << 1646 } << 1647 for (ii = 1; ii < maxbin; ++ii) << 1648 { << 1649 vals[ii] = vals[ii] / (bins[ii] - bin << 1650 } << 1651 sum = vals[maxbin - 1]; << 1652 vals[0] = 0.; << 1653 } << 1654 for (ii = 0; ii < maxbin; ++ii) << 1655 { << 1656 vals[ii] = vals[ii] / sum; << 1657 IPDFEnergyH.InsertValues(bins[ii], vals << 1658 } << 1659 << 1660 IPDFEnergyExist = true; << 1661 if (verbosityLevel > 1) << 1662 { << 1663 IPDFEnergyH.DumpValues(); << 1664 } << 1665 } << 1666 l.unlock(); << 1667 << 1668 // IPDF has been create so carry on << 1669 // << 1670 G4double rndm = eneRndm->GenRandEnergy(); << 1671 threadLocalData.Get().particle_energy= IPDF << 1672 << 1673 if (verbosityLevel >= 1) << 1674 { << 1675 G4cout << "Energy is " << particle_energy << 1676 } << 1677 } << 1678 << 1679 G4double G4SPSEneDistribution::GetArbEneWeigh << 1680 { << 1681 auto nbelow = IPDFArbEnergyH.FindBin(ene,(I << 1682 G4double wei = 0.; << 1683 if(IntType=="Lin") << 1684 { << 1685 // note: grad[i] and cept[i] are calculat << 1686 auto gr = Arb_grad[nbelow + 1]; << 1687 auto ce = Arb_cept[nbelow + 1]; << 1688 wei = ene*gr + ce; << 1689 } << 1690 else if(IntType=="Log") << 1691 { << 1692 auto alp = Arb_alpha[nbelow + 1]; << 1693 auto cns = Arb_Const[nbelow + 1]; << 1694 wei = cns * std::pow(ene,alp); << 1695 } << 1696 else if(IntType=="Exp") << 1697 { << 1698 auto e0 = Arb_ezero[nbelow + 1]; << 1699 auto cns = Arb_Const[nbelow + 1]; << 1700 wei = cns * std::exp(-ene/e0); << 1701 } << 1702 else if(IntType=="Spline") << 1703 { << 1704 wei = SplineInt[nbelow+1]->CubicSplineInt << 1705 } << 1706 return wei; << 1707 } << 1708 << 1709 void G4SPSEneDistribution::GenArbPointEnergie << 1710 { << 1711 if (verbosityLevel > 0) << 1712 { << 1713 G4cout << "In GenArbPointEnergies" << G4e << 1714 } << 1715 << 1716 G4double rndm = eneRndm->GenRandEnergy(); << 1717 << 1718 // Find the Bin, have x, y, no of points, a << 1719 // << 1720 std::size_t nabove = IPDFArbEnergyH.GetVect << 1721 << 1722 // Binary search to find bin that rndm is i << 1723 // << 1724 while (nabove - nbelow > 1) << 1725 { << 1726 middle = (nabove + nbelow) / 2; << 1727 if (rndm == IPDFArbEnergyH(middle)) << 1728 { << 1729 break; << 1730 } << 1731 if (rndm < IPDFArbEnergyH(middle)) << 1732 { << 1733 nabove = middle; << 1734 } << 1735 else << 1736 { << 1737 nbelow = middle; << 1738 } << 1739 } << 1740 threadLocal_t& params = threadLocalData.Get << 1741 if (IntType == "Lin") << 1742 { << 1743 // Update thread-local copy of parameters << 1744 // << 1745 params.Emax = IPDFArbEnergyH.GetLowEdgeEn << 1746 params.Emin = IPDFArbEnergyH.GetLowEdgeEn << 1747 params.grad = Arb_grad[nbelow + 1]; << 1748 params.cept = Arb_cept[nbelow + 1]; << 1749 GenerateLinearEnergies(true); << 1750 } << 1751 else if (IntType == "Log") << 1752 { << 1753 params.Emax = IPDFArbEnergyH.GetLowEdgeEn << 1754 params.Emin = IPDFArbEnergyH.GetLowEdgeEn << 1755 params.alpha = Arb_alpha[nbelow + 1]; << 1756 GeneratePowEnergies(true); << 1757 } << 1758 else if (IntType == "Exp") << 1759 { << 1760 params.Emax = IPDFArbEnergyH.GetLowEdgeEn << 1761 params.Emin = IPDFArbEnergyH.GetLowEdgeEn << 1762 params.Ezero = Arb_ezero[nbelow + 1]; << 1763 GenerateExpEnergies(true); << 1764 } << 1765 else if (IntType == "Spline") << 1766 { << 1767 params.Emax = IPDFArbEnergyH.GetLowEdgeEn << 1768 params.Emin = IPDFArbEnergyH.GetLowEdgeEn << 1769 params.particle_energy = -1e100; << 1770 rndm = eneRndm->GenRandEnergy(); << 1771 while (params.particle_energy < params.Em << 1772 || params.particle_energy > params.Em << 1773 { << 1774 params.particle_energy = << 1775 SplineInt[nbelow+1]->CubicSplineInter << 1776 rndm = eneRndm->GenRandEnergy(); << 1777 } << 1778 if (verbosityLevel >= 1) << 1779 { << 1780 G4cout << "Energy is " << params.partic << 1781 } << 1782 } << 1783 else << 1784 { << 1785 G4Exception("G4SPSEneDistribution::GenArb << 1786 FatalException, "Error: IntTy << 1787 } << 1788 } << 1789 << 1790 void G4SPSEneDistribution::GenEpnHistEnergies << 1791 { << 1792 // Firstly convert to energy if not already << 1793 << 1794 G4AutoLock l(&mutex); << 1795 << 1796 if (Epnflag) // true means spectrum is epn << 1797 { << 1798 // Convert to energy by multiplying by A << 1799 // << 1800 ConvertEPNToEnergy(); << 1801 } << 1802 if (!IPDFEnergyExist) << 1803 { << 1804 // IPDF has not been created, so create i << 1805 // << 1806 G4double bins[1024], vals[1024], sum; << 1807 std::size_t ii; << 1808 std::size_t maxbin = UDefEnergyH.GetVecto << 1809 bins[0] = UDefEnergyH.GetLowEdgeEnergy(0) << 1810 vals[0] = UDefEnergyH(0); << 1811 sum = vals[0]; << 1812 for (ii = 1; ii < maxbin; ++ii) << 1813 { << 1814 bins[ii] = UDefEnergyH.GetLowEdgeEnergy << 1815 vals[ii] = UDefEnergyH(ii) + vals[ii - << 1816 sum = sum + UDefEnergyH(ii); << 1817 } << 1818 << 1819 l.lock(); << 1820 for (ii = 0; ii < maxbin; ++ii) << 1821 { << 1822 vals[ii] = vals[ii] / sum; << 1823 IPDFEnergyH.InsertValues(bins[ii], vals << 1824 } << 1825 IPDFEnergyExist = true; << 1826 << 1827 } << 1828 l.unlock(); << 1829 612 1830 // IPDF has been create so carry on << 613 G4int maxi = ArbEnergyH.GetVectorLength(); 1831 // << 614 for (i = 0; i < maxi; i++) { 1832 G4double rndm = eneRndm->GenRandEnergy(); << 615 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); 1833 threadLocalData.Get().particle_energy = IPD << 616 Arb_y[i] = ArbEnergyH(size_t(i)); 1834 << 617 } 1835 if (verbosityLevel >= 1) << 618 // Points are now in x,y arrays. If the spectrum is integral it has to be 1836 { << 619 // made differential and if momentum it has to be made energy. 1837 G4cout << "Energy is " << threadLocalData << 620 if (DiffSpec == false) { 1838 } << 621 // Converts integral point-wise spectra to Differential 1839 } << 622 for (count = 0; count < maxi - 1; count++) { 1840 << 623 Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) 1841 void G4SPSEneDistribution::ConvertEPNToEnergy << 624 / (Arb_x[count + 1] - Arb_x[count]); 1842 { << 625 } 1843 // Use this before particle generation to c << 626 maxi--; 1844 // currently stored histogram from energy/n << 627 } 1845 << 628 // 1846 threadLocal_t& params = threadLocalData.Get << 629 if (EnergySpec == false) { 1847 if (params.particle_definition == nullptr) << 630 // change currently stored values (emin etc) which are actually momenta 1848 { << 631 // to energies. 1849 G4cout << "Error: particle not defined" < << 632 if (particle_definition == NULL) 1850 } << 633 G4cout << "Error: particle not defined" << G4endl; 1851 else << 634 else { 1852 { << 635 // Apply Energy**2 = p**2c**2 + m0**2c**4 1853 // Need to multiply histogram by the numb << 636 // p should be entered as E/c i.e. without the division by c 1854 // Baryon Number looks to hold the no. of << 637 // being done - energy equivalent. 1855 // << 638 G4double mass = particle_definition->GetPDGMass(); 1856 G4int Bary = params.particle_definition-> << 639 // convert point to energy unit and its value to per energy unit 1857 << 640 G4double total_energy; 1858 // Change values in histogram, Read it ou << 641 for (count = 0; count < maxi; count++) { 1859 // << 642 total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass 1860 std::size_t count, maxcount; << 643 * mass)); // total energy 1861 maxcount = EpnEnergyH.GetVectorLength(); << 644 1862 G4double ebins[1024], evals[1024]; << 645 Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy; 1863 if (maxcount > 1024) << 646 Arb_x[count] = total_energy - mass; // kinetic energy 1864 { << 647 } 1865 G4Exception("G4SPSEneDistribution::Conv << 648 } 1866 "gps001", JustWarning, << 649 } 1867 "Histogram contains more th << 650 1868 Those above 1024 will be i << 651 // 1869 maxcount = 1024; << 652 i = 1; 1870 } << 653 Arb_Cum_Area[0] = 0.; 1871 if (maxcount < 1) << 654 sum = 0.; 1872 { << 655 Splinetemp = new G4DataInterpolation(Arb_x, Arb_y, maxi, 0., 0.); 1873 G4Exception("G4SPSEneDistribution::Conv << 656 G4double ei[101],prob[101]; 1874 "gps001", FatalException, << 657 while (i < maxi) { 1875 "Histogram contains less tha << 658 // 100 step per segment for the integration of area 1876 return; << 659 G4double de = (Arb_x[i] - Arb_x[i - 1])/100.; 1877 } << 660 G4double area = 0.; 1878 for (count = 0; count < maxcount; ++count << 661 1879 { << 662 for (count = 0; count < 101; count++) { 1880 // Read out << 663 ei[count] = Arb_x[i - 1] + de*count ; 1881 ebins[count] = EpnEnergyH.GetLowEdgeEne << 664 prob[count] = Splinetemp->CubicSplineInterpolation(ei[count]); 1882 evals[count] = EpnEnergyH(count); << 665 if (prob[count] < 0.) { 1883 } << 666 G4cout << "Warning: G4DataInterpolation returns value < 0 " << prob[count] <<" "<<ei[count]<< G4endl; 1884 << 667 G4Exception(" Please use an alternative method, e.g. Lin, for interpolation"); 1885 // Multiply the channels by the nucleon n << 668 } 1886 // << 669 area += prob[count]*de; 1887 for (count = 0; count < maxcount; ++count << 670 } 1888 { << 671 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + area; 1889 ebins[count] = ebins[count] * Bary; << 672 sum += area; 1890 } << 673 1891 << 674 prob[0] = prob[0]/(area/de); 1892 // Set Emin and Emax << 675 for (count = 1; count < 100; count++) 1893 // << 676 prob[count] = prob[count-1] + prob[count]/(area/de); 1894 params.Emin = ebins[0]; << 677 1895 if (maxcount > 1) << 678 SplineInt[i] = new G4DataInterpolation(prob, ei, 101, 0., 0.); 1896 { << 679 // note i start from 1! 1897 params.Emax = ebins[maxcount - 1]; << 680 i++; 1898 } << 681 } 1899 else << 682 i = 0; 1900 { << 683 while (i < maxi) { 1901 params.Emax = ebins[0]; << 684 Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation 1902 } << 685 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); 1903 << 686 i++; 1904 // Put energy bins into new histogram - U << 687 } 1905 // << 688 // now scale the ArbEnergyH, needed by Probability() 1906 for (count = 0; count < maxcount; ++count << 689 ArbEnergyH.ScaleVector(1., 1./sum); 1907 { << 690 1908 UDefEnergyH.InsertValues(ebins[count], << 691 if (verbosityLevel > 0) 1909 } << 692 G4cout << "Leaving SplineInterpolation " << G4endl; 1910 Epnflag = false; // so that you dont repe << 693 } 1911 } << 694 1912 } << 695 void G4SPSEneDistribution::GenerateMonoEnergetic() { 1913 << 696 // Method to generate MonoEnergetic particles. 1914 void G4SPSEneDistribution::ReSetHist(const G4 << 697 particle_energy = MonoEnergy; 1915 { << 698 } 1916 G4AutoLock l(&mutex); << 699 1917 if (atype == "energy") << 700 void G4SPSEneDistribution::GenerateGaussEnergies() { 1918 { << 701 // Method to generate Gaussian particles. 1919 UDefEnergyH = IPDFEnergyH = ZeroPhysVecto << 702 particle_energy = G4RandGauss::shoot(MonoEnergy,SE); 1920 IPDFEnergyExist = false; << 703 if (particle_energy < 0) particle_energy = 0.; 1921 Emin = 0.; << 704 } 1922 Emax = 1e30; << 705 1923 } << 706 void G4SPSEneDistribution::GenerateLinearEnergies(G4bool bArb = false) { 1924 else if (atype == "arb") << 707 G4double rndm; 1925 { << 708 G4double emaxsq = std::pow(Emax, 2.); //Emax squared 1926 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVec << 709 G4double eminsq = std::pow(Emin, 2.); //Emin squared 1927 IPDFArbExist = false; << 710 G4double intersq = std::pow(cept, 2.); //cept squared 1928 } << 711 1929 else if (atype == "epn") << 712 if (bArb) 1930 { << 713 rndm = G4UniformRand(); 1931 UDefEnergyH = IPDFEnergyH = ZeroPhysVecto << 714 else 1932 IPDFEnergyExist = false; << 715 rndm = eneRndm->GenRandEnergy(); 1933 EpnEnergyH = ZeroPhysVector; << 716 1934 } << 717 G4double bracket = ((grad / 2.) * (emaxsq - eminsq) + cept * (Emax - Emin)); 1935 else << 718 bracket = bracket * rndm; 1936 { << 719 bracket = bracket + (grad / 2.) * eminsq + cept * Emin; 1937 G4cout << "Error, histtype not accepted " << 720 // Now have a quad of form m/2 E**2 + cE - bracket = 0 1938 } << 721 bracket = -bracket; 1939 } << 722 // G4cout << "BRACKET" << bracket << G4endl; 1940 << 723 if (grad != 0.) { 1941 G4double G4SPSEneDistribution::GenerateOne(G4 << 724 G4double sqbrack = (intersq - 4 * (grad / 2.) * (bracket)); 1942 { << 725 // G4cout << "SQBRACK" << sqbrack << G4endl; 1943 // Copy global shared status to thread-loca << 726 sqbrack = std::sqrt(sqbrack); 1944 // << 727 G4double root1 = -cept + sqbrack; 1945 threadLocal_t& params = threadLocalData.Get << 728 root1 = root1 / (2. * (grad / 2.)); 1946 params.particle_definition=a; << 729 1947 params.particle_energy=-1; << 730 G4double root2 = -cept - sqbrack; 1948 if(applyEvergyWeight) << 731 root2 = root2 / (2. * (grad / 2.)); 1949 { << 732 1950 params.Emax = ArbEmax; << 733 // G4cout << root1 << " roots " << root2 << G4endl; 1951 params.Emin = ArbEmin; << 734 1952 } << 735 if (root1 > Emin && root1 < Emax) 1953 else << 736 particle_energy = root1; 1954 { << 737 if (root2 > Emin && root2 < Emax) 1955 params.Emax = Emax; << 738 particle_energy = root2; 1956 params.Emin = Emin; << 739 } else if (grad == 0.) 1957 } << 740 // have equation of form cE - bracket =0 1958 params.alpha = alpha; << 741 particle_energy = bracket / cept; 1959 params.Ezero = Ezero; << 742 1960 params.grad = grad; << 743 if (particle_energy < 0.) 1961 params.cept = cept; << 744 particle_energy = -particle_energy; 1962 params.weight = weight; << 745 1963 // particle_energy = -1.; << 746 if (verbosityLevel >= 1) 1964 << 747 G4cout << "Energy is " << particle_energy << G4endl; 1965 if((EnergyDisType == "Mono") && ((MonoEnerg << 748 } 1966 { << 749 1967 G4ExceptionDescription ed; << 750 void G4SPSEneDistribution::GeneratePowEnergies(G4bool bArb = false) { 1968 ed << "MonoEnergy " << G4BestUnit(MonoEne << 751 // Method to generate particle energies distributed as 1969 << " is outside of [Emin,Emax] = [" << 752 // a power-law 1970 << G4BestUnit(Emin,"Energy") << ", " << 753 1971 << G4BestUnit(Emax,"Energy") << ". Mon << 754 G4double rndm; 1972 G4Exception("G4SPSEneDistribution::Genera << 755 G4double emina, emaxa; 1973 "GPS0001", JustWarning, ed); << 756 1974 params.particle_energy=MonoEnergy; << 757 emina = std::pow(Emin, alpha + 1); 1975 return params.particle_energy; << 758 emaxa = std::pow(Emax, alpha + 1); 1976 } << 759 1977 while ( (EnergyDisType == "Arb") << 760 if (bArb) 1978 ? (params.particle_energy < ArbEmin << 761 rndm = G4UniformRand(); 1979 || params.particle_energy > ArbEmax << 762 else 1980 : (params.particle_energy < params. << 763 rndm = eneRndm->GenRandEnergy(); 1981 || params.particle_energy > params. << 764 1982 { << 765 if (alpha != -1.) { 1983 if (Biased) << 766 particle_energy = ((rndm * (emaxa - emina)) + emina); 1984 { << 767 particle_energy = std::pow(particle_energy, (1. / (alpha + 1.))); 1985 GenerateBiasPowEnergies(); << 768 } else { 1986 } << 769 particle_energy = (std::log(Emin) + rndm * (std::log(Emax) - std::log( 1987 else << 770 Emin))); 1988 { << 771 particle_energy = std::exp(particle_energy); 1989 if (EnergyDisType == "Mono") << 772 } 1990 { << 773 if (verbosityLevel >= 1) 1991 GenerateMonoEnergetic(); << 774 G4cout << "Energy is " << particle_energy << G4endl; 1992 } << 775 } 1993 else if (EnergyDisType == "Lin") << 776 1994 { << 777 void G4SPSEneDistribution::GenerateBiasPowEnergies() { 1995 GenerateLinearEnergies(false); << 778 // Method to generate particle energies distributed as 1996 } << 779 // in biased power-law and calculate its weight 1997 else if (EnergyDisType == "Pow") << 780 1998 { << 781 G4double rndm; 1999 GeneratePowEnergies(false); << 782 G4double emina, emaxa, emin, emax; 2000 } << 783 2001 else if (EnergyDisType == "CPow") << 784 G4double normal = 1. ; 2002 { << 785 2003 GenerateCPowEnergies(); << 786 emin = Emin; 2004 } << 787 emax = Emax; 2005 else if (EnergyDisType == "Exp") << 788 // if (EnergyDisType == "Arb") { 2006 { << 789 // emin = ArbEmin; 2007 GenerateExpEnergies(false); << 790 // emax = ArbEmax; 2008 } << 791 //} 2009 else if (EnergyDisType == "Gauss") << 792 2010 { << 793 rndm = eneRndm->GenRandEnergy(); 2011 GenerateGaussEnergies(); << 794 2012 } << 795 if (biasalpha != -1.) { 2013 else if (EnergyDisType == "Brem") << 796 emina = std::pow(emin, biasalpha + 1); 2014 { << 797 emaxa = std::pow(emax, biasalpha + 1); 2015 GenerateBremEnergies(); << 798 particle_energy = ((rndm * (emaxa - emina)) + emina); 2016 } << 799 particle_energy = std::pow(particle_energy, (1. / (biasalpha + 1.))); 2017 else if (EnergyDisType == "Bbody") << 800 normal = 1./(1+biasalpha) * (emaxa - emina); 2018 { << 801 } else { 2019 GenerateBbodyEnergies(); << 802 particle_energy = (std::log(emin) + rndm * (std::log(emax) - std::log( 2020 } << 803 emin))); 2021 else if (EnergyDisType == "Cdg") << 804 particle_energy = std::exp(particle_energy); 2022 { << 805 normal = std::log(emax) - std::log(emin) ; 2023 GenerateCdgEnergies(); << 806 } 2024 } << 807 weight = GetProbability(particle_energy) / (std::pow(particle_energy,biasalpha)/normal); 2025 else if (EnergyDisType == "User") << 808 2026 { << 809 if (verbosityLevel >= 1) 2027 GenUserHistEnergies(); << 810 G4cout << "Energy is " << particle_energy << G4endl; 2028 } << 811 } 2029 else if (EnergyDisType == "Arb") << 812 2030 { << 813 void G4SPSEneDistribution::GenerateExpEnergies(G4bool bArb = false) { 2031 GenArbPointEnergies(); << 814 // Method to generate particle energies distributed according 2032 } << 815 // to an exponential curve. 2033 else if (EnergyDisType == "Epn") << 816 G4double rndm; 2034 { << 817 2035 GenEpnHistEnergies(); << 818 if (bArb) 2036 } << 819 rndm = G4UniformRand(); 2037 else << 820 else 2038 { << 821 rndm = eneRndm->GenRandEnergy(); 2039 G4cout << "Error: EnergyDisType has u << 822 2040 } << 823 particle_energy = -Ezero * (std::log(rndm * (std::exp(-Emax / Ezero) 2041 } << 824 - std::exp(-Emin / Ezero)) + std::exp(-Emin / Ezero))); 2042 } << 825 if (verbosityLevel >= 1) 2043 return params.particle_energy; << 826 G4cout << "Energy is " << particle_energy << G4endl; 2044 } << 827 } 2045 << 828 2046 G4double G4SPSEneDistribution::GetProbability << 829 void G4SPSEneDistribution::GenerateBremEnergies() { 2047 { << 830 // Method to generate particle energies distributed according 2048 G4double prob = 1.; << 831 // to a Bremstrahlung equation of 2049 << 832 // form I = const*((kT)**1/2)*E*(e**(-E/kT)) 2050 threadLocal_t& params = threadLocalData.Get << 833 2051 if (EnergyDisType == "Lin") << 834 G4double rndm; 2052 { << 835 rndm = eneRndm->GenRandEnergy(); 2053 if (prob_norm == 1.) << 836 G4double expmax, expmin, k; 2054 { << 837 2055 prob_norm = 0.5*params.grad*params.Emax << 838 k = 8.6181e-11; // Boltzmann's const in MeV/K 2056 + params.cept*params.Emax << 839 G4double ksq = std::pow(k, 2.); // k squared 2057 - 0.5*params.grad*params.Emin << 840 G4double Tsq = std::pow(Temp, 2.); // Temp squared 2058 - params.cept*params.Emin; << 841 2059 } << 842 expmax = std::exp(-Emax / (k * Temp)); 2060 prob = params.cept + params.grad * ene; << 843 expmin = std::exp(-Emin / (k * Temp)); 2061 prob /= prob_norm; << 844 2062 } << 845 // If either expmax or expmin are zero then this will cause problems 2063 else if (EnergyDisType == "Pow") << 846 // Most probably this will be because T is too low or E is too high 2064 { << 847 2065 if (prob_norm == 1.) << 848 if (expmax == 0.) 2066 { << 849 G4cout << "*****EXPMAX=0. Choose different E's or Temp" << G4endl; 2067 if (alpha != -1.) << 850 if (expmin == 0.) 2068 { << 851 G4cout << "*****EXPMIN=0. Choose different E's or Temp" << G4endl; 2069 G4double emina = std::pow(params.Emin << 852 2070 G4double emaxa = std::pow(params.Emax << 853 G4double tempvar = rndm * ((-k) * Temp * (Emax * expmax - Emin * expmin) 2071 prob_norm = 1./(1.+alpha) * (emaxa - << 854 - (ksq * Tsq * (expmax - expmin))); 2072 } << 855 2073 else << 856 G4double bigc = (tempvar - k * Temp * Emin * expmin - ksq * Tsq * expmin) 2074 { << 857 / (-k * Temp); 2075 prob_norm = std::log(params.Emax) - s << 858 2076 } << 859 // This gives an equation of form: Ee(-E/kT) + kTe(-E/kT) - C =0 2077 } << 860 // Solve this iteratively, step from Emin to Emax in 1000 steps 2078 prob = std::pow(ene, params.alpha)/prob_n << 861 // and take the best solution. 2079 } << 862 2080 else if (EnergyDisType == "Exp") << 863 G4double erange = Emax - Emin; 2081 { << 864 G4double steps = erange / 1000.; 2082 if (prob_norm == 1.) << 865 G4int i; 2083 { << 866 G4double etest, diff, err; 2084 prob_norm = -params.Ezero*(std::exp(-pa << 867 2085 - std::exp(par << 868 err = 100000.; 2086 } << 869 2087 prob = std::exp(-ene / params.Ezero); << 870 for (i = 1; i < 1000; i++) { 2088 prob /= prob_norm; << 871 etest = Emin + (i - 1) * steps; 2089 } << 872 2090 else if (EnergyDisType == "Arb") << 873 diff = etest * (std::exp(-etest / (k * Temp))) + k * Temp * (std::exp( 2091 { << 874 -etest / (k * Temp))) - bigc; 2092 prob = ArbEnergyH.Value(ene); << 875 2093 << 876 if (diff < 0.) 2094 if (prob <= 0.) << 877 diff = -diff; 2095 { << 878 2096 G4cout << " Warning:G4SPSEneDistributio << 879 if (diff < err) { 2097 << prob << " " << ene << G4endl; << 880 err = diff; 2098 prob = 1e-30; << 881 particle_energy = etest; 2099 } << 882 } 2100 } << 883 } 2101 else << 884 if (verbosityLevel >= 1) 2102 { << 885 G4cout << "Energy is " << particle_energy << G4endl; 2103 G4cout << "Error: EnergyDisType not suppo << 886 } 2104 } << 887 >> 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(); >> 897 >> 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; >> 907 } >> 908 >> 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; >> 922 } >> 923 } >> 924 >> 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; >> 939 } >> 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; >> 946 } >> 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; >> 967 } >> 968 >> 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(); >> 1028 } >> 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)); >> 1061 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); >> 1062 grad = Arb_grad[nbelow + 1]; >> 1063 cept = Arb_cept[nbelow + 1]; >> 1064 // G4cout << rndm << " " << Emax << " " << Emin << " " << grad << " " << cept << G4endl; >> 1065 GenerateLinearEnergies(true); >> 1066 } else if (IntType == "Log") { >> 1067 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1)); >> 1068 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); >> 1069 alpha = Arb_alpha[nbelow + 1]; >> 1070 // G4cout << rndm << " " << Emax << " " << Emin << " " << alpha << G4endl; >> 1071 GeneratePowEnergies(true); >> 1072 } else if (IntType == "Exp") { >> 1073 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1)); >> 1074 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); >> 1075 Ezero = Arb_ezero[nbelow + 1]; >> 1076 // G4cout << rndm << " " << Emax << " " << Emin << " " << Ezero << G4endl; >> 1077 GenerateExpEnergies(true); >> 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; >> 1091 } >> 1092 >> 1093 void G4SPSEneDistribution::GenEpnHistEnergies() { >> 1094 // G4cout << "In GenEpnHistEnergies " << Epnflag << G4endl; >> 1095 >> 1096 // Firstly convert to energy if not already done. >> 1097 if (Epnflag == true) >> 1098 // epnflag = true means spectrum is epn, false means e. >> 1099 { >> 1100 // convert to energy by multiplying by A number >> 1101 ConvertEPNToEnergy(); >> 1102 // EpnEnergyH will be replace by UDefEnergyH. >> 1103 // UDefEnergyH.DumpValues(); >> 1104 } >> 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; >> 1127 } >> 1128 // IPDFEnergyH.DumpValues(); >> 1129 // IPDF has been create so carry on >> 1130 G4double rndm = eneRndm->GenRandEnergy(); >> 1131 particle_energy = IPDFEnergyH.GetEnergy(rndm); >> 1132 >> 1133 if (verbosityLevel >= 1) >> 1134 G4cout << "Energy is " << particle_energy << G4endl; >> 1135 } >> 1136 >> 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 } >> 1169 >> 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 } >> 1174 >> 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. >> 1186 } >> 1187 } 2105 1188 2106 return prob; << 1189 // >> 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 } >> 1206 } >> 1207 >> 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; >> 1245 } >> 1246 >> 1247 G4double G4SPSEneDistribution::GetProbability(G4double ene) { >> 1248 G4double prob = 1.; >> 1249 >> 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; 2107 } 1292 } 2108 1293