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 // 44 << 45 // Initialise all variables 55 // Initialise all variables >> 56 particle_energy = 1.0*MeV; 46 57 47 particle_energy = 1.0 * MeV; << 48 EnergyDisType = "Mono"; 58 EnergyDisType = "Mono"; 49 weight=1.; << 59 MonoEnergy = 1*MeV; 50 MonoEnergy = 1 * MeV; << 51 Emin = 0.; 60 Emin = 0.; 52 Emax = 1.e30; 61 Emax = 1.e30; 53 alpha = 0.; 62 alpha = 0.; 54 biasalpha = 0.; << 55 prob_norm = 1.0; << 56 Ezero = 0.; 63 Ezero = 0.; 57 SE = 0.; 64 SE = 0.; 58 Temp = 0.; 65 Temp = 0.; 59 grad = 0.; 66 grad = 0.; 60 cept = 0.; 67 cept = 0.; >> 68 EnergySpec = true; // true - energy spectra, false - momentum spectra >> 69 DiffSpec = true; // true - differential spec, false integral spec 61 IntType = "NULL"; // Interpolation type 70 IntType = "NULL"; // Interpolation type >> 71 IPDFEnergyExist = false; >> 72 IPDFArbExist = false; 62 73 63 ArbEmin = 0.; 74 ArbEmin = 0.; 64 ArbEmax = 1.e30; 75 ArbEmax = 1.e30; 65 76 66 verbosityLevel = 0; << 77 verbosityLevel = 0 ; 67 << 78 68 threadLocal_t& data = threadLocalData.Get(); << 69 data.Emax = Emax; << 70 data.Emin = Emin; << 71 data.alpha =alpha; << 72 data.cept = cept; << 73 data.Ezero = Ezero; << 74 data.grad = grad; << 75 data.particle_energy = 0.; << 76 data.particle_definition = nullptr; << 77 data.weight = weight; << 78 } 79 } 79 80 80 G4SPSEneDistribution::~G4SPSEneDistribution() 81 G4SPSEneDistribution::~G4SPSEneDistribution() 81 { << 82 {} 82 G4MUTEXDESTROY(mutex); << 83 if(Arb_grad_cept_flag) << 84 { << 85 delete [] Arb_grad; << 86 delete [] Arb_cept; << 87 } << 88 << 89 if(Arb_alpha_Const_flag) << 90 { << 91 delete [] Arb_alpha; << 92 delete [] Arb_Const; << 93 } << 94 << 95 if(Arb_ezero_flag) << 96 { << 97 delete [] Arb_ezero; << 98 } << 99 delete Bbody_x; << 100 delete BBHist; << 101 delete CP_x; << 102 delete CPHist; << 103 for (auto & it : SplineInt) << 104 { << 105 delete it; << 106 it = nullptr; << 107 } << 108 SplineInt.clear(); << 109 } << 110 83 111 void G4SPSEneDistribution::SetEnergyDisType(co << 84 void G4SPSEneDistribution::SetEnergyDisType(G4String DisType) 112 { 85 { 113 G4AutoLock l(&mutex); << 114 EnergyDisType = DisType; 86 EnergyDisType = DisType; 115 if (EnergyDisType == "User") << 87 if (EnergyDisType == "User"){ 116 { << 88 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ; 117 UDefEnergyH = IPDFEnergyH = ZeroPhysVector << 89 IPDFEnergyExist = false ; 118 IPDFEnergyExist = false; << 90 } else if ( EnergyDisType == "Arb"){ 119 } << 91 ArbEnergyH =IPDFArbEnergyH = ZeroPhysVector ; 120 else if (EnergyDisType == "Arb") << 121 { << 122 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVect << 123 IPDFArbExist = false; 92 IPDFArbExist = false; 124 } << 93 } else if (EnergyDisType == "Epn"){ 125 else if (EnergyDisType == "Epn") << 94 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ; 126 { << 95 IPDFEnergyExist = false ; 127 UDefEnergyH = IPDFEnergyH = ZeroPhysVector << 96 EpnEnergyH = ZeroPhysVector ; 128 IPDFEnergyExist = false; << 129 EpnEnergyH = ZeroPhysVector; << 130 } 97 } 131 } 98 } 132 99 133 const G4String& G4SPSEneDistribution::GetEnerg << 134 { << 135 G4AutoLock l(&mutex); << 136 return EnergyDisType; << 137 } << 138 << 139 void G4SPSEneDistribution::SetEmin(G4double em 100 void G4SPSEneDistribution::SetEmin(G4double emi) 140 { 101 { 141 G4AutoLock l(&mutex); << 142 Emin = emi; 102 Emin = emi; 143 threadLocalData.Get().Emin = Emin; << 144 } << 145 << 146 G4double G4SPSEneDistribution::GetEmin() const << 147 { << 148 return threadLocalData.Get().Emin; << 149 } << 150 << 151 G4double G4SPSEneDistribution::GetArbEmin() << 152 { << 153 G4AutoLock l(&mutex); << 154 return ArbEmin; << 155 } << 156 << 157 G4double G4SPSEneDistribution::GetArbEmax() << 158 { << 159 G4AutoLock l(&mutex); << 160 return ArbEmax; << 161 } 103 } 162 104 163 void G4SPSEneDistribution::SetEmax(G4double em 105 void G4SPSEneDistribution::SetEmax(G4double ema) 164 { 106 { 165 G4AutoLock l(&mutex); << 166 Emax = ema; 107 Emax = ema; 167 threadLocalData.Get().Emax = Emax; << 168 } << 169 << 170 G4double G4SPSEneDistribution::GetEmax() const << 171 { << 172 return threadLocalData.Get().Emax; << 173 } 108 } 174 109 175 void G4SPSEneDistribution::SetMonoEnergy(G4dou 110 void G4SPSEneDistribution::SetMonoEnergy(G4double menergy) 176 { 111 { 177 G4AutoLock l(&mutex); << 178 MonoEnergy = menergy; 112 MonoEnergy = menergy; 179 } 113 } 180 114 181 void G4SPSEneDistribution::SetBeamSigmaInE(G4d 115 void G4SPSEneDistribution::SetBeamSigmaInE(G4double e) 182 { 116 { 183 G4AutoLock l(&mutex); << 184 SE = e; 117 SE = e; 185 } 118 } 186 void G4SPSEneDistribution::SetAlpha(G4double a 119 void G4SPSEneDistribution::SetAlpha(G4double alp) 187 { 120 { 188 G4AutoLock l(&mutex); << 189 alpha = alp; 121 alpha = alp; 190 threadLocalData.Get().alpha = alpha; << 191 } << 192 << 193 void G4SPSEneDistribution::SetBiasAlpha(G4doub << 194 { << 195 G4AutoLock l(&mutex); << 196 biasalpha = alp; << 197 Biased = true; << 198 } 122 } 199 123 200 void G4SPSEneDistribution::SetTemp(G4double te 124 void G4SPSEneDistribution::SetTemp(G4double tem) 201 { 125 { 202 G4AutoLock l(&mutex); << 203 Temp = tem; 126 Temp = tem; 204 } 127 } 205 128 206 void G4SPSEneDistribution::SetEzero(G4double e 129 void G4SPSEneDistribution::SetEzero(G4double eze) 207 { 130 { 208 G4AutoLock l(&mutex); << 209 Ezero = eze; 131 Ezero = eze; 210 threadLocalData.Get().Ezero = Ezero; << 211 } 132 } 212 133 213 void G4SPSEneDistribution::SetGradient(G4doubl 134 void G4SPSEneDistribution::SetGradient(G4double gr) 214 { 135 { 215 G4AutoLock l(&mutex); << 216 grad = gr; 136 grad = gr; 217 threadLocalData.Get().grad = grad; << 218 } 137 } 219 138 220 void G4SPSEneDistribution::SetInterCept(G4doub 139 void G4SPSEneDistribution::SetInterCept(G4double c) 221 { 140 { 222 G4AutoLock l(&mutex); << 223 cept = c; 141 cept = c; 224 threadLocalData.Get().cept = cept; << 225 } << 226 << 227 const G4String& G4SPSEneDistribution::GetIntTy << 228 { << 229 G4AutoLock l(&mutex); << 230 return IntType; << 231 } << 232 << 233 void G4SPSEneDistribution::SetBiasRndm(G4SPSRa << 234 { << 235 G4AutoLock l(&mutex); << 236 eneRndm = a; << 237 } << 238 << 239 void G4SPSEneDistribution::SetVerbosity(G4int << 240 { << 241 G4AutoLock l(&mutex); << 242 verbosityLevel = a; << 243 } << 244 << 245 G4double G4SPSEneDistribution::GetWeight() con << 246 { << 247 return threadLocalData.Get().weight; << 248 } << 249 << 250 G4double G4SPSEneDistribution::GetMonoEnergy() << 251 { << 252 G4AutoLock l(&mutex); << 253 return MonoEnergy; << 254 } << 255 << 256 G4double G4SPSEneDistribution::GetSE() << 257 { << 258 G4AutoLock l(&mutex); << 259 return SE; << 260 } << 261 << 262 G4double G4SPSEneDistribution::Getalpha() cons << 263 { << 264 return threadLocalData.Get().alpha; << 265 } << 266 << 267 G4double G4SPSEneDistribution::GetEzero() cons << 268 { << 269 return threadLocalData.Get().Ezero; << 270 } << 271 << 272 G4double G4SPSEneDistribution::GetTemp() << 273 { << 274 G4AutoLock l(&mutex); << 275 return Temp; << 276 } << 277 << 278 G4double G4SPSEneDistribution::Getgrad() const << 279 { << 280 return threadLocalData.Get().grad; << 281 } 142 } 282 143 283 G4double G4SPSEneDistribution::Getcept() const << 144 void G4SPSEneDistribution::UserEnergyHisto(G4ThreeVector input) 284 { 145 { 285 return threadLocalData.Get().cept; << 146 G4double ehi, val; 286 } << 147 ehi = input.x(); 287 << 148 val = input.y(); 288 G4PhysicsFreeVector G4SPSEneDistribution::GetU << 149 if(verbosityLevel > 1) { 289 { << 290 G4AutoLock l(&mutex); << 291 return UDefEnergyH; << 292 } << 293 << 294 G4PhysicsFreeVector G4SPSEneDistribution::GetA << 295 { << 296 G4AutoLock l(&mutex); << 297 return ArbEnergyH; << 298 } << 299 << 300 void G4SPSEneDistribution::UserEnergyHisto(con << 301 { << 302 G4AutoLock l(&mutex); << 303 G4double ehi = input.x(), << 304 val = input.y(); << 305 if (verbosityLevel > 1) << 306 { << 307 G4cout << "In UserEnergyHisto" << G4endl; 150 G4cout << "In UserEnergyHisto" << G4endl; 308 G4cout << " " << ehi << " " << val << G4en 151 G4cout << " " << ehi << " " << val << G4endl; 309 } 152 } 310 UDefEnergyH.InsertValues(ehi, val); 153 UDefEnergyH.InsertValues(ehi, val); 311 Emax = ehi; 154 Emax = ehi; 312 threadLocalData.Get().Emax = Emax; << 313 } 155 } 314 156 315 void G4SPSEneDistribution::ArbEnergyHisto(cons << 157 void G4SPSEneDistribution::ArbEnergyHisto(G4ThreeVector input) 316 { 158 { 317 G4AutoLock l(&mutex); << 159 G4double ehi, val; 318 G4double ehi = input.x(), << 160 ehi = input.x(); 319 val = input.y(); << 161 val = input.y(); 320 if (verbosityLevel > 1) << 162 if(verbosityLevel >1 ) { 321 { << 322 G4cout << "In ArbEnergyHisto" << G4endl; 163 G4cout << "In ArbEnergyHisto" << G4endl; 323 G4cout << " " << ehi << " " << val << G4en 164 G4cout << " " << ehi << " " << val << G4endl; 324 } 165 } 325 ArbEnergyH.InsertValues(ehi, val); 166 ArbEnergyH.InsertValues(ehi, val); 326 } 167 } 327 168 328 void G4SPSEneDistribution::ArbEnergyHistoFile( << 169 void G4SPSEneDistribution::EpnEnergyHisto(G4ThreeVector input) 329 { 170 { 330 G4AutoLock l(&mutex); << 331 std::ifstream infile(filename, std::ios::in) << 332 if (!infile) << 333 { << 334 G4Exception("G4SPSEneDistribution::ArbEner << 335 FatalException, "Unable to ope << 336 } << 337 G4double ehi, val; 171 G4double ehi, val; 338 while (infile >> ehi >> val) << 172 ehi = input.x(); 339 { << 173 val = input.y(); 340 ArbEnergyH.InsertValues(ehi, val); << 174 if(verbosityLevel > 1) { 341 } << 342 } << 343 << 344 void G4SPSEneDistribution::EpnEnergyHisto(cons << 345 { << 346 G4AutoLock l(&mutex); << 347 G4double ehi = input.x(), << 348 val = input.y(); << 349 if (verbosityLevel > 1) << 350 { << 351 G4cout << "In EpnEnergyHisto" << G4endl; 175 G4cout << "In EpnEnergyHisto" << G4endl; 352 G4cout << " " << ehi << " " << val << G4en 176 G4cout << " " << ehi << " " << val << G4endl; 353 } 177 } 354 EpnEnergyH.InsertValues(ehi, val); 178 EpnEnergyH.InsertValues(ehi, val); 355 Emax = ehi; 179 Emax = ehi; 356 threadLocalData.Get().Emax = Emax; << 357 Epnflag = true; 180 Epnflag = true; 358 } 181 } 359 182 360 void G4SPSEneDistribution::Calculate() 183 void G4SPSEneDistribution::Calculate() 361 { 184 { 362 G4AutoLock l(&mutex); << 185 if(EnergyDisType == "Cdg") 363 if (EnergyDisType == "Cdg") << 364 { << 365 CalculateCdgSpectrum(); 186 CalculateCdgSpectrum(); 366 } << 187 else if(EnergyDisType == "Bbody") 367 else if (EnergyDisType == "Bbody") << 368 { << 369 if(!BBhistInit) << 370 { << 371 BBInitHists(); << 372 } << 373 CalculateBbodySpectrum(); 188 CalculateBbodySpectrum(); 374 } << 375 else if (EnergyDisType == "CPow") << 376 { << 377 if(!CPhistInit) << 378 { << 379 CPInitHists(); << 380 } << 381 CalculateCPowSpectrum(); << 382 } << 383 } 189 } 384 190 385 void G4SPSEneDistribution::BBInitHists() // M << 191 void G4SPSEneDistribution::CalculateCdgSpectrum() 386 { 192 { 387 BBHist = new std::vector<G4double>(10001, 0. << 193 // This uses the spectrum from The INTEGRAL Mass Model (TIMM) 388 Bbody_x = new std::vector<G4double>(10001, 0 << 389 BBhistInit = true; << 390 } << 391 << 392 void G4SPSEneDistribution::CPInitHists() // M << 393 { << 394 CPHist = new std::vector<G4double>(10001, 0. << 395 CP_x = new std::vector<G4double>(10001, 0.0) << 396 CPhistInit = true; << 397 } << 398 << 399 void G4SPSEneDistribution::CalculateCdgSpectru << 400 { << 401 // This uses the spectrum from the INTEGRAL << 402 // to generate a Cosmic Diffuse X/gamma ray 194 // to generate a Cosmic Diffuse X/gamma ray spectrum. 403 << 195 G4double pfact[2] = {8.5, 112}; 404 G4double pfact[2] = { 8.5, 112 }; << 196 G4double spind[2] = {1.4, 2.3}; 405 G4double spind[2] = { 1.4, 2.3 }; << 197 G4double ene_line[3] = {1.*keV, 18.*keV, 1E6*keV}; 406 G4double ene_line[3] = { 1. * keV, 18. * keV << 407 G4int n_par; 198 G4int n_par; 408 199 409 ene_line[0] = threadLocalData.Get().Emin; << 200 ene_line[0] = Emin; 410 if (threadLocalData.Get().Emin < 18 * keV) << 201 if(Emin < 18*keV) 411 { << 412 n_par = 2; << 413 ene_line[2] = threadLocalData.Get().Emax; << 414 if (threadLocalData.Get().Emax < 18 * keV) << 415 { 202 { 416 n_par = 1; << 203 n_par = 2; 417 ene_line[1] = threadLocalData.Get().Emax << 204 ene_line[2] = Emax; >> 205 if(Emax < 18*keV) >> 206 { >> 207 n_par = 1; >> 208 ene_line[1] = Emax; >> 209 } 418 } 210 } 419 } << 420 else 211 else 421 { << 212 { 422 n_par = 1; << 213 n_par = 1; 423 pfact[0] = 112.; << 214 pfact[0] = 112.; 424 spind[0] = 2.3; << 215 spind[0] = 2.3; 425 ene_line[1] = threadLocalData.Get().Emax; << 216 ene_line[1] = Emax; 426 } << 217 } 427 << 218 428 // Create a cumulative histogram << 219 // Create a cumulative histogram. 429 // << 430 CDGhist[0] = 0.; 220 CDGhist[0] = 0.; 431 G4double omalpha; 221 G4double omalpha; 432 G4int i = 0; 222 G4int i = 0; 433 while (i < n_par) << 434 { << 435 omalpha = 1. - spind[i]; << 436 CDGhist[i + 1] = CDGhist[i] + (pfact[i] / << 437 * (std::pow(en << 438 - std::pow(ene << 439 ++i; << 440 } << 441 223 442 // Normalise histo and divide by 1000 to mak << 224 while(i < n_par) 443 // << 225 { >> 226 omalpha = 1. - spind[i]; >> 227 CDGhist[i+1] = CDGhist[i] + (pfact[i]/omalpha)* >> 228 (std::pow(ene_line[i+1],omalpha)-std::pow(ene_line[i],omalpha)); >> 229 i++; >> 230 } >> 231 >> 232 // Normalise histo and divide by 1000 to make MeV. 444 i = 0; 233 i = 0; 445 while (i < n_par) << 234 while(i < n_par) 446 { << 235 { 447 CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[ << 236 CDGhist[i+1] = CDGhist[i+1]/CDGhist[n_par]; 448 ++i; << 237 // G4cout << CDGhist[i] << CDGhist[n_par] << G4endl; 449 } << 238 i++; >> 239 } 450 } 240 } 451 241 452 void G4SPSEneDistribution::CalculateBbodySpect << 242 void G4SPSEneDistribution::CalculateBbodySpectrum() 453 { 243 { 454 // Create bbody spectrum << 244 // create bbody spectrum 455 // Proved very hard to integrate indefinitel << 245 // Proved very hard to integrate indefinitely, so different 456 // User inputs emin, emax and T. These are u << 246 // method. User inputs emin, emax and T. These are used to 457 // bin histogram. << 247 // create a 10,000 bin histogram. 458 // Use photon density spectrum = 2 nu**2/c** 248 // Use photon density spectrum = 2 nu**2/c**2 * (std::exp(h nu/kT)-1) 459 // = 2 E**2/h**2c**2 times the exponential 249 // = 2 E**2/h**2c**2 times the exponential 460 << 250 G4double erange = Emax - Emin; 461 G4double erange = threadLocalData.Get().Emax << 251 G4double steps = erange/10000.; 462 G4double steps = erange / 10000.; << 252 G4double Bbody_y[10000]; 463 << 253 G4double k = 8.6181e-11; //Boltzmann const in MeV/K 464 const G4double k = 8.6181e-11; //Boltzmann c << 254 G4double h = 4.1362e-21; // Plancks const in MeV s 465 const G4double h = 4.1362e-21; // Plancks co << 255 G4double c = 3e8; // Speed of light 466 const G4double c = 3e8; // Speed of light << 256 G4double h2 = h*h; 467 const G4double h2 = h * h; << 257 G4double c2 = c*c; 468 const G4double c2 = c * c; << 469 G4int count = 0; << 470 G4double sum = 0.; << 471 BBHist->at(0) = 0.; << 472 << 473 while (count < 10000) << 474 { << 475 Bbody_x->at(count) = threadLocalData.Get() << 476 G4double Bbody_y = (2. * std::pow(Bbody_x- << 477 / (h2*c2*(std::exp(Bbody_ << 478 sum = sum + Bbody_y; << 479 BBHist->at(count + 1) = BBHist->at(count) << 480 ++count; << 481 } << 482 << 483 Bbody_x->at(10000) = threadLocalData.Get().E << 484 << 485 // Normalise cumulative histo << 486 // << 487 count = 0; << 488 while (count < 10001) << 489 { << 490 BBHist->at(count) = BBHist->at(count) / su << 491 ++count; << 492 } << 493 } << 494 << 495 void G4SPSEneDistribution::CalculateCPowSpectr << 496 { << 497 // Create cutoff power-law spectrum, x^a exp << 498 // The integral of this function is an incom << 499 // is only available in the Boost library. << 500 // << 501 // User inputs are emin, emax and alpha and << 502 // create a 10,000 bin histogram. << 503 << 504 G4double erange = threadLocalData.Get().Emax << 505 G4double steps = erange / 10000.; << 506 alpha = threadLocalData.Get().alpha ; << 507 Ezero = threadLocalData.Get().Ezero ; << 508 << 509 G4int count = 0; 258 G4int count = 0; 510 G4double sum = 0.; 259 G4double sum = 0.; 511 CPHist->at(0) = 0.; << 260 BBHist[0] = 0.; 512 << 261 while(count < 10000) 513 while (count < 10000) << 262 { 514 { << 263 Bbody_x[count] = Emin + G4double(count*steps); 515 CP_x->at(count) = threadLocalData.Get().Em << 264 Bbody_y[count] = (2.*std::pow(Bbody_x[count],2.))/ 516 G4double CP_y = std::pow(CP_x->at(count), << 265 (h2*c2*(std::exp(Bbody_x[count]/(k*Temp)) - 1.)); 517 * std::exp(-CP_x->at(count) << 266 sum = sum + Bbody_y[count]; 518 sum = sum + CP_y; << 267 BBHist[count+1] = BBHist[count] + Bbody_y[count]; 519 CPHist->at(count + 1) = CPHist->at(count) << 268 count++; 520 ++count; << 269 } 521 } << 522 << 523 CP_x->at(10000) = threadLocalData.Get().Emax << 524 270 525 // Normalise cumulative histo << 271 Bbody_x[10000] = Emax; 526 // << 272 // Normalise cumulative histo. 527 count = 0; 273 count = 0; 528 while (count < 10001) << 274 while(count<10001) 529 { << 275 { 530 CPHist->at(count) = CPHist->at(count) / su << 276 BBHist[count] = BBHist[count]/sum; 531 ++count; << 277 count++; 532 } << 278 } 533 } 279 } 534 280 535 void G4SPSEneDistribution::InputEnergySpectra( 281 void G4SPSEneDistribution::InputEnergySpectra(G4bool value) 536 { 282 { 537 G4AutoLock l(&mutex); << 283 // Allows user to specifiy spectrum is momentum 538 << 539 // Allows user to specify spectrum is moment << 540 // << 541 EnergySpec = value; // false if momentum 284 EnergySpec = value; // false if momentum 542 if (verbosityLevel > 1) << 285 if(verbosityLevel > 1) 543 { << 544 G4cout << "EnergySpec has value " << Energ 286 G4cout << "EnergySpec has value " << EnergySpec << G4endl; 545 } << 546 } 287 } 547 288 548 void G4SPSEneDistribution::InputDifferentialSp 289 void G4SPSEneDistribution::InputDifferentialSpectra(G4bool value) 549 { 290 { 550 G4AutoLock l(&mutex); << 551 << 552 // Allows user to specify integral or differ 291 // Allows user to specify integral or differential spectra 553 // << 554 DiffSpec = value; // true = differential, fa 292 DiffSpec = value; // true = differential, false = integral 555 if (verbosityLevel > 1) << 293 if(verbosityLevel > 1) 556 { << 557 G4cout << "Diffspec has value " << DiffSpe 294 G4cout << "Diffspec has value " << DiffSpec << G4endl; 558 } << 559 } 295 } 560 296 561 void G4SPSEneDistribution::ArbInterpolate(cons << 297 void G4SPSEneDistribution::ArbInterpolate(G4String IType) 562 { 298 { 563 G4AutoLock l(&mutex); << 299 if(EnergyDisType != "Arb") 564 << 300 G4cout << "Error: this is for arbitrary distributions" << G4endl; 565 IntType = IType; 301 IntType = IType; 566 ArbEmax = ArbEnergyH.GetMaxEnergy(); << 302 ArbEmax = Emax; 567 ArbEmin = ArbEnergyH.Energy(0); << 303 ArbEmin = Emin; 568 304 569 // Now interpolate points 305 // Now interpolate points 570 << 306 if(IntType == "Lin") 571 if (IntType == "Lin") LinearInterpolation(); << 307 LinearInterpolation(); 572 if (IntType == "Log") LogInterpolation(); << 308 if(IntType == "Log") 573 if (IntType == "Exp") ExpInterpolation(); << 309 LogInterpolation(); 574 if (IntType == "Spline") SplineInterpolation << 310 if(IntType == "Exp") >> 311 ExpInterpolation(); >> 312 if(IntType == "Spline") >> 313 SplineInterpolation(); 575 } 314 } 576 315 577 void G4SPSEneDistribution::LinearInterpolation << 316 void G4SPSEneDistribution::LinearInterpolation() 578 { 317 { 579 // Method to do linear interpolation on the 318 // Method to do linear interpolation on the Arb points 580 // Calculate equation of each line segment, 319 // Calculate equation of each line segment, max 1024. 581 // Calculate Area under each segment 320 // Calculate Area under each segment 582 // Create a cumulative array which is then n 321 // Create a cumulative array which is then normalised Arb_Cum_Area 583 322 584 G4double Area_seg[1024]; // Stores area unde 323 G4double Area_seg[1024]; // Stores area under each segment 585 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 324 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024]; 586 std::size_t i, count; << 325 G4int i, count; 587 std::size_t maxi = ArbEnergyH.GetVectorLengt << 326 G4int maxi = ArbEnergyH.GetVectorLength(); 588 for (i = 0; i < maxi; ++i) << 327 for(i=0;i<maxi;i++) { 589 { << 328 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); 590 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 329 Arb_y[i] = ArbEnergyH(size_t(i)); 591 Arb_y[i] = ArbEnergyH(i); << 592 } 330 } 593 << 594 // Points are now in x,y arrays. If the spec 331 // Points are now in x,y arrays. If the spectrum is integral it has to be 595 // made differential and if momentum it has << 332 // made differential and if momentum it has to be made energy. 596 << 333 if(DiffSpec == false) { 597 if (!DiffSpec) << 598 { << 599 // Converts integral point-wise spectra to 334 // Converts integral point-wise spectra to Differential 600 // << 335 for( count=0;count < maxi-1;count++) { 601 for (count = 0; count < maxi-1; ++count) << 336 Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]); 602 { << 603 Arb_y[count] = (Arb_y[count] - Arb_y[cou << 604 / (Arb_x[count + 1] - Arb_x << 605 } 337 } 606 --maxi; << 338 maxi--; 607 } 339 } 608 << 340 // 609 if (!EnergySpec) << 341 if(EnergySpec == false) { 610 { << 611 // change currently stored values (emin et 342 // change currently stored values (emin etc) which are actually momenta 612 // to energies << 343 // to energies. 613 // << 344 if(particle_definition == NULL) 614 G4ParticleDefinition* pdef = threadLocalDa << 345 G4cout << "Error: particle not defined" << G4endl; 615 if (pdef == nullptr) << 346 else { 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** 347 // Apply Energy**2 = p**2c**2 + m0**2c**4 624 // p should be entered as E/c i.e. witho 348 // p should be entered as E/c i.e. without the division by c 625 // being done - energy equivalent << 349 // being done - energy equivalent. 626 << 350 G4double mass = particle_definition->GetPDGMass(); 627 G4double mass = pdef->GetPDGMass(); << 351 // convert point to energy unit and its value to per energy unit 628 << 629 // Convert point to energy unit and its << 630 // << 631 G4double total_energy; 352 G4double total_energy; 632 for (count = 0; count < maxi; ++count) << 353 for(count=0;count<maxi;count++) { 633 { << 354 total_energy = std::sqrt((Arb_x[count]*Arb_x[count]) 634 total_energy = std::sqrt((Arb_x[count] << 355 + (mass*mass)); // total energy 635 + (mass * mass)); // tota << 356 636 Arb_y[count] = Arb_y[count] * Arb_x[co << 357 Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy; 637 Arb_x[count] = total_energy - mass; // << 358 Arb_x[count] = total_energy - mass ; // kinetic energy 638 } 359 } 639 } 360 } 640 } 361 } 641 << 362 // 642 i = 1; << 363 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.; 364 Arb_grad[0] = 0.; 649 Arb_cept[0] = 0.; 365 Arb_cept[0] = 0.; 650 Area_seg[0] = 0.; 366 Area_seg[0] = 0.; 651 Arb_Cum_Area[0] = 0.; 367 Arb_Cum_Area[0] = 0.; 652 while (i < maxi) << 368 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 { 369 { 671 if (verbosityLevel == 2) << 370 // calc gradient and intercept for each segment 672 { << 371 Arb_grad[i] = (Arb_y[i] - Arb_y[i-1]) / (Arb_x[i] - Arb_x[i-1]); 673 G4cout << "Arb_grad is negative" << G4 << 372 if(verbosityLevel == 2) 674 } << 373 G4cout << Arb_grad[i] << G4endl; 675 Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * << 374 if(Arb_grad[i] > 0.) >> 375 { >> 376 if(verbosityLevel == 2) >> 377 G4cout << "Arb_grad is positive" << G4endl; >> 378 Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]); >> 379 } >> 380 else if(Arb_grad[i] < 0.) >> 381 { >> 382 if(verbosityLevel == 2) >> 383 G4cout << "Arb_grad is negative" << G4endl; >> 384 Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]); >> 385 } >> 386 else >> 387 { >> 388 if(verbosityLevel == 2) >> 389 G4cout << "Arb_grad is 0." << G4endl; >> 390 Arb_cept[i] = Arb_y[i]; >> 391 } >> 392 >> 393 Area_seg[i] = ((Arb_grad[i]/2)*(Arb_x[i]*Arb_x[i] - Arb_x[i-1]*Arb_x[i-1]) + Arb_cept[i]*(Arb_x[i] - Arb_x[i-1])); >> 394 Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i]; >> 395 sum = sum + Area_seg[i]; >> 396 if(verbosityLevel == 2) >> 397 G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i] << G4endl; >> 398 i++; 676 } 399 } 677 else << 400 >> 401 i=0; >> 402 while(i < maxi) 678 { 403 { 679 if (verbosityLevel == 2) << 404 Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum; // normalisation 680 { << 405 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); 681 G4cout << "Arb_grad is 0." << G4endl; << 406 i++; 682 } << 683 Arb_cept[i] = Arb_y[i]; << 684 } 407 } 685 408 686 Area_seg[i] = ((Arb_grad[i] / 2) << 409 if(verbosityLevel >= 1) 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 { 410 { 693 G4cout << Arb_x[i] << Arb_y[i] << Area_s << 411 G4cout << "Leaving LinearInterpolation" << G4endl; 694 << Arb_grad[i] << G4endl; << 412 ArbEnergyH.DumpValues(); >> 413 IPDFArbEnergyH.DumpValues(); 695 } 414 } 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 } 415 } 718 416 719 void G4SPSEneDistribution::LogInterpolation() << 417 void G4SPSEneDistribution::LogInterpolation() 720 { 418 { 721 // Interpolation based on Logarithmic equati 419 // Interpolation based on Logarithmic equations 722 // Generate equations of line segments 420 // Generate equations of line segments 723 // y = Ax**alpha => log y = alpha*logx + log 421 // y = Ax**alpha => log y = alpha*logx + logA 724 // Find area under line segments 422 // Find area under line segments 725 // Create normalised, cumulative array Arb_C << 423 // create normalised, cumulative array Arb_Cum_Area 726 << 727 G4double Area_seg[1024]; // Stores area unde 424 G4double Area_seg[1024]; // Stores area under each segment 728 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 425 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024]; 729 std::size_t i, count; << 426 G4int i, count; 730 std::size_t maxi = ArbEnergyH.GetVectorLengt << 427 G4int maxi = ArbEnergyH.GetVectorLength(); 731 for (i = 0; i < maxi; ++i) << 428 for(i=0;i<maxi;i++) { 732 { << 429 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); 733 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 430 Arb_y[i] = ArbEnergyH(size_t(i)); 734 Arb_y[i] = ArbEnergyH(i); << 735 } 431 } 736 << 737 // Points are now in x,y arrays. If the spec 432 // Points are now in x,y arrays. If the spectrum is integral it has to be 738 // made differential and if momentum it has << 433 // made differential and if momentum it has to be made energy. 739 << 434 if(DiffSpec == false) { 740 if (!DiffSpec) << 741 { << 742 // Converts integral point-wise spectra to 435 // Converts integral point-wise spectra to Differential 743 // << 436 for( count=0;count<maxi-1;count++) { 744 for (count = 0; count < maxi-1; ++count) << 437 Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]); 745 { << 746 Arb_y[count] = (Arb_y[count] - Arb_y[cou << 747 / (Arb_x[count + 1] - Arb_x << 748 } 438 } 749 --maxi; << 439 maxi--; 750 } 440 } 751 << 441 // 752 if (!EnergySpec) << 442 if(EnergySpec == false) { 753 { << 443 // change currently stored values (emin etc) which are actually momenta 754 // Change currently stored values (emin et << 444 // to energies. 755 // to energies << 445 if(particle_definition == NULL) 756 << 446 G4cout << "Error: particle not defined" << G4endl; 757 G4ParticleDefinition* pdef = threadLocalDa << 447 else { 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** 448 // Apply Energy**2 = p**2c**2 + m0**2c**4 767 // p should be entered as E/c i.e. witho 449 // p should be entered as E/c i.e. without the division by c 768 // being done - energy equivalent << 450 // being done - energy equivalent. 769 << 451 G4double mass = particle_definition->GetPDGMass(); 770 G4double mass = pdef->GetPDGMass(); << 452 // convert point to energy unit and its value to per energy unit 771 << 772 // Convert point to energy unit and its << 773 // << 774 G4double total_energy; 453 G4double total_energy; 775 for (count = 0; count < maxi; ++count) << 454 for(count=0;count<maxi;count++) { 776 { << 455 total_energy = std::sqrt((Arb_x[count]*Arb_x[count]) 777 total_energy = std::sqrt((Arb_x[count] << 456 + (mass*mass)); // total energy 778 Arb_y[count] = Arb_y[count] * Arb_x[co << 457 779 Arb_x[count] = total_energy - mass; // << 458 Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy; >> 459 Arb_x[count] = total_energy - mass ; // kinetic energy 780 } 460 } 781 } 461 } 782 } 462 } 783 << 463 // 784 i = 1; << 464 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.; 465 Arb_alpha[0] = 0.; 793 Arb_Const[0] = 0.; 466 Arb_Const[0] = 0.; 794 Area_seg[0] = 0.; 467 Area_seg[0] = 0.; 795 Arb_Cum_Area[0] = 0.; << 468 if(Arb_x[0] <= 0. || Arb_y[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 { 469 { 808 Arb_y[0] = 1e-20; << 470 G4cout << "You should not use log interpolation with points <= 0." << G4endl; >> 471 G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl; >> 472 if(Arb_x[0] <= 0.) >> 473 Arb_x[0] = 1e-20; >> 474 if(Arb_y[0] <= 0.) >> 475 Arb_y[0] = 1e-20; >> 476 } >> 477 >> 478 G4double alp; >> 479 while(i <maxi) >> 480 { >> 481 // Incase points are negative or zero >> 482 if(Arb_x[i] <= 0. || Arb_y[i] <= 0.) >> 483 { >> 484 G4cout << "You should not use log interpolation with points <= 0." << G4endl; >> 485 G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl; >> 486 if(Arb_x[i] <= 0.) >> 487 Arb_x[i] = 1e-20; >> 488 if(Arb_y[i] <= 0.) >> 489 Arb_y[i] = 1e-20; >> 490 } >> 491 >> 492 Arb_alpha[i] = (std::log10(Arb_y[i])-std::log10(Arb_y[i-1]))/(std::log10(Arb_x[i])-std::log10(Arb_x[i-1])); >> 493 Arb_Const[i] = Arb_y[i]/(std::pow(Arb_x[i],Arb_alpha[i])); >> 494 alp = Arb_alpha[i] + 1; >> 495 Area_seg[i] = (Arb_Const[i]/alp) * (std::pow(Arb_x[i],alp) - std::pow(Arb_x[i-1],alp)); >> 496 sum = sum + Area_seg[i]; >> 497 Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i]; >> 498 if(verbosityLevel == 2) >> 499 G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl; >> 500 i++; 809 } 501 } 810 } << 502 811 << 503 i=0; 812 G4double alp; << 504 while(i<maxi) 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 { 505 { 851 G4cout << Arb_alpha[i] << Arb_Const[i] < << 506 Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum; >> 507 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); >> 508 i++; 852 } 509 } 853 ++i; << 510 if(verbosityLevel >= 1) 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 511 G4cout << "Leaving LogInterpolation " << G4endl; 871 } << 872 } 512 } 873 513 874 void G4SPSEneDistribution::ExpInterpolation() << 514 void G4SPSEneDistribution::ExpInterpolation() 875 { 515 { 876 // Interpolation based on Exponential equati 516 // Interpolation based on Exponential equations 877 // Generate equations of line segments 517 // Generate equations of line segments 878 // y = Ae**-(x/e0) => ln y = -x/e0 + lnA 518 // y = Ae**-(x/e0) => ln y = -x/e0 + lnA 879 // Find area under line segments 519 // Find area under line segments 880 // Create normalised, cumulative array Arb_C << 520 // create normalised, cumulative array Arb_Cum_Area 881 << 882 G4double Area_seg[1024]; // Stores area unde 521 G4double Area_seg[1024]; // Stores area under each segment 883 G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 522 G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024]; 884 std::size_t i, count; << 523 G4int i, count; 885 std::size_t maxi = ArbEnergyH.GetVectorLengt << 524 G4int maxi = ArbEnergyH.GetVectorLength(); 886 for (i = 0; i < maxi; ++i) << 525 for(i=0;i<maxi;i++) { 887 { << 526 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); 888 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 527 Arb_y[i] = ArbEnergyH(size_t(i)); 889 Arb_y[i] = ArbEnergyH(i); << 890 } 528 } 891 << 892 // Points are now in x,y arrays. If the spec 529 // Points are now in x,y arrays. If the spectrum is integral it has to be 893 // made differential and if momentum it has << 530 // made differential and if momentum it has to be made energy. 894 << 531 if(DiffSpec == false) { 895 if (!DiffSpec) << 896 { << 897 // Converts integral point-wise spectra to 532 // Converts integral point-wise spectra to Differential 898 // << 533 for( count=0;count< maxi-1;count++) { 899 for (count = 0; count < maxi - 1; ++count) << 534 Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]); 900 { << 901 Arb_y[count] = (Arb_y[count] - Arb_y[cou << 902 / (Arb_x[count + 1] - Arb_x << 903 } 535 } 904 --maxi; << 536 maxi--; 905 } 537 } 906 << 538 // 907 if (!EnergySpec) << 539 if(EnergySpec == false) { 908 { << 540 // change currently stored values (emin etc) which are actually momenta 909 // Change currently stored values (emin et << 541 // to energies. 910 // to energies << 542 if(particle_definition == NULL) 911 // << 543 G4cout << "Error: particle not defined" << G4endl; 912 G4ParticleDefinition* pdef = threadLocalDa << 544 else { 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** 545 // Apply Energy**2 = p**2c**2 + m0**2c**4 922 // p should be entered as E/c i.e. witho 546 // p should be entered as E/c i.e. without the division by c 923 // being done - energy equivalent << 547 // being done - energy equivalent. 924 << 548 G4double mass = particle_definition->GetPDGMass(); 925 G4double mass = pdef->GetPDGMass(); << 549 // convert point to energy unit and its value to per energy unit 926 << 927 // Convert point to energy unit and its << 928 // << 929 G4double total_energy; 550 G4double total_energy; 930 for (count = 0; count < maxi; ++count) << 551 for(count=0;count<maxi;count++) { 931 { << 552 total_energy = std::sqrt((Arb_x[count]*Arb_x[count]) 932 total_energy = std::sqrt((Arb_x[count] << 553 + (mass*mass)); // total energy 933 Arb_y[count] = Arb_y[count] * Arb_x[co << 554 934 Arb_x[count] = total_energy - mass; // << 555 Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy; >> 556 Arb_x[count] = total_energy - mass ; // kinetic energy 935 } 557 } 936 } 558 } 937 } 559 } 938 << 560 // 939 i = 1; << 561 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.; 562 Arb_ezero[0] = 0.; 948 Arb_Const[0] = 0.; 563 Arb_Const[0] = 0.; 949 Area_seg[0] = 0.; 564 Area_seg[0] = 0.; 950 Arb_Cum_Area[0] = 0.; 565 Arb_Cum_Area[0] = 0.; 951 while (i < maxi) << 566 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 { 567 { 965 G4Exception("G4SPSEneDistribution::ExpIn << 568 G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i-1]); 966 "Event0302", JustWarning, << 569 if(test > 0. || test < 0.) 967 "Flat line segment: problem, << 570 { 968 G4cout << "Flat line segment: problem" < << 571 Arb_ezero[i] = -(Arb_x[i] - Arb_x[i-1])/(std::log(Arb_y[i]) - std::log(Arb_y[i-1])); 969 Arb_ezero[i] = 0.; << 572 Arb_Const[i] = Arb_y[i]/(std::exp(-Arb_x[i]/Arb_ezero[i])); 970 Arb_Const[i] = 0.; << 573 Area_seg[i]=-(Arb_Const[i]*Arb_ezero[i])*(std::exp(-Arb_x[i]/Arb_ezero[i]) 971 Area_seg[i] = 0.; << 574 -std::exp(-Arb_x[i-1]/Arb_ezero[i])); 972 } << 575 } 973 sum = sum + Area_seg[i]; << 576 else 974 Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Ar << 577 { 975 if (verbosityLevel == 2) << 578 G4cout << "Flat line segment: problem" << G4endl; >> 579 Arb_ezero[i] = 0.; >> 580 Arb_Const[i] = 0.; >> 581 Area_seg[i] = 0.; >> 582 } >> 583 sum = sum + Area_seg[i]; >> 584 Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i]; >> 585 if(verbosityLevel == 2) >> 586 G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl; >> 587 i++; >> 588 } >> 589 >> 590 i=0; >> 591 while(i<maxi) 976 { 592 { 977 G4cout << Arb_ezero[i] << Arb_Const[i] < << 593 Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum; >> 594 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); >> 595 i++; 978 } 596 } 979 ++i; << 597 if(verbosityLevel >= 1) 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 598 G4cout << "Leaving ExpInterpolation " << G4endl; 997 } << 998 } 599 } 999 600 1000 void G4SPSEneDistribution::SplineInterpolatio << 601 void G4SPSEneDistribution::SplineInterpolation() 1001 { 602 { 1002 // Interpolation using Splines. 603 // Interpolation using Splines. 1003 // Create Normalised arrays, make x 0->1 an << 604 // Create Normalised arrays, make x 0->1 and y hold 1004 // << 605 // the function (Energy) 1005 // Current method based on the above will n << 606 G4double Arb_x[1024], Arb_y[1024]; 1006 // New method is implemented below. << 607 G4int i, count; 1007 << 608 G4int maxi = ArbEnergyH.GetVectorLength(); 1008 G4double sum, Arb_x[1024]={0.}, Arb_y[1024] << 609 for(i=0;i<maxi;i++) { 1009 std::size_t i, count; << 610 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); 1010 std::size_t maxi = ArbEnergyH.GetVectorLeng << 611 Arb_y[i] = ArbEnergyH(size_t(i)); 1011 << 1012 for (i = 0; i < maxi; ++i) << 1013 { << 1014 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i) << 1015 Arb_y[i] = ArbEnergyH(i); << 1016 } 612 } 1017 << 1018 // Points are now in x,y arrays. If the spe 613 // Points are now in x,y arrays. If the spectrum is integral it has to be 1019 // made differential and if momentum it has << 614 // made differential and if momentum it has to be made energy. 1020 << 615 if(DiffSpec == false) { 1021 if (!DiffSpec) << 1022 { << 1023 // Converts integral point-wise spectra t 616 // Converts integral point-wise spectra to Differential 1024 // << 617 for( count=0;count< maxi-1;count++) { 1025 for (count = 0; count < maxi - 1; ++count << 618 Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]); 1026 { << 1027 Arb_y[count] = (Arb_y[count] - Arb_y[co << 1028 / (Arb_x[count + 1] - Arb_ << 1029 } 619 } 1030 --maxi; << 620 maxi--; 1031 } 621 } 1032 << 622 // 1033 if (!EnergySpec) << 623 if(EnergySpec == false) { 1034 { << 624 // change currently stored values (emin etc) which are actually momenta 1035 // Change currently stored values (emin e << 625 // to energies. 1036 // to energies << 626 if(particle_definition == NULL) 1037 // << 627 G4cout << "Error: particle not defined" << G4endl; 1038 G4ParticleDefinition* pdef = threadLocalD << 628 else { 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* 629 // Apply Energy**2 = p**2c**2 + m0**2c**4 1048 // p should be entered as E/c i.e. with 630 // p should be entered as E/c i.e. without the division by c 1049 // being done - energy equivalent << 631 // being done - energy equivalent. 1050 << 632 G4double mass = particle_definition->GetPDGMass(); 1051 G4double mass = pdef->GetPDGMass(); << 633 // convert point to energy unit and its value to per energy unit 1052 << 1053 // Convert point to energy unit and its << 1054 // << 1055 G4double total_energy; 634 G4double total_energy; 1056 for (count = 0; count < maxi; ++count) << 635 for(count=0;count<maxi;count++) { 1057 { << 636 total_energy = std::sqrt((Arb_x[count]*Arb_x[count]) 1058 total_energy = std::sqrt((Arb_x[count << 637 + (mass*mass)); // total energy 1059 Arb_y[count] = Arb_y[count] * Arb_x[c << 638 1060 Arb_x[count] = total_energy - mass; / << 639 Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy; >> 640 Arb_x[count] = total_energy - mass ; // kinetic energy 1061 } 641 } 1062 } 642 } 1063 } 643 } 1064 << 644 // 1065 i = 1; << 645 for(i=1;i<maxi;i++) 1066 Arb_Cum_Area[0] = 0.; << 646 Arb_y[i] += Arb_y[i-1]; 1067 sum = 0.; << 647 1068 Splinetemp = new G4DataInterpolation(Arb_x, << 648 for(i=0;i<maxi;i++) 1069 G4double ei[101], prob[101]; << 649 Arb_y[i] /= Arb_y[maxi-1]; 1070 for (auto & it : SplineInt) << 650 // now Arb_y is accumulated normalised probabilities 1071 { << 651 /* for(i=0; i<maxi;i++) { 1072 delete it; << 652 if(verbosityLevel >1) 1073 it = nullptr; << 653 G4cout << i <<" "<< Arb_x[i] << " " << Arb_y[i] << G4endl; 1074 } << 654 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_y[i]); 1075 SplineInt.clear(); << 655 } 1076 SplineInt.resize(1024,nullptr); << 656 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(IPDFArbEnergyH.GetVectorLength()-1); 1077 while (i < maxi) << 657 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(0); 1078 { << 658 */ 1079 // 100 step per segment for the integrati << 659 // Should now have normalised cumulative probabilities in Arb_y 1080 << 660 // and energy values in Arb_x. 1081 G4double de = (Arb_x[i] - Arb_x[i - 1])/1 << 661 // maxi = maxi + 1; 1082 G4double area = 0.; << 662 // Put y into x and x into y. The spline interpolation will then 1083 << 663 // go through x-axis to find where to interpolate (cum probability) 1084 for (count = 0; count < 101; ++count) << 664 // then generate a y (which will now be energy). 1085 { << 665 SplineInt = new G4DataInterpolation(Arb_y,Arb_x,maxi,1e30,1e30); 1086 ei[count] = Arb_x[i - 1] + de*count ; << 666 if(verbosityLevel >1 ) 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 { 667 { 1104 prob[count] = prob[count-1] + prob[coun << 668 G4cout << SplineInt << G4endl; >> 669 G4cout << SplineInt->LocateArgument(1.0) << G4endl; 1105 } 670 } 1106 << 671 if(verbosityLevel > 0 ) 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 " 672 G4cout << "Leaving SplineInterpolation " << G4endl; 1129 } << 1130 } 673 } 1131 674 1132 void G4SPSEneDistribution::GenerateMonoEnerge 675 void G4SPSEneDistribution::GenerateMonoEnergetic() 1133 { 676 { 1134 // Method to generate MonoEnergetic particl << 677 // Method to generate MonoEnergetic particles. 1135 << 678 particle_energy = MonoEnergy; 1136 threadLocalData.Get().particle_energy = Mon << 1137 } 679 } 1138 680 1139 void G4SPSEneDistribution::GenerateGaussEnerg 681 void G4SPSEneDistribution::GenerateGaussEnergies() 1140 { 682 { 1141 // Method to generate Gaussian particles << 683 // Method to generate Gaussian particles. 1142 << 684 particle_energy = G4RandGauss::shoot(MonoEnergy,SE); 1143 G4double ene = G4RandGauss::shoot(MonoEnerg << 685 if (particle_energy < 0) particle_energy = 0.; 1144 if (ene < 0) ene = 0.; << 1145 threadLocalData.Get().particle_energy = ene << 1146 } 686 } 1147 687 1148 void G4SPSEneDistribution::GenerateLinearEner 688 void G4SPSEneDistribution::GenerateLinearEnergies(G4bool bArb = false) 1149 { 689 { 1150 G4double rndm; 690 G4double rndm; 1151 threadLocal_t& params = threadLocalData.Get << 691 G4double emaxsq = std::pow(Emax,2.); //Emax squared 1152 G4double emaxsq = std::pow(params.Emax, 2.) << 692 G4double eminsq = std::pow(Emin,2.); //Emin squared 1153 G4double eminsq = std::pow(params.Emin, 2.) << 693 G4double intersq = std::pow(cept,2.); //cept squared 1154 G4double intersq = std::pow(params.cept, 2. << 1155 694 1156 if (bArb) rndm = G4UniformRand(); 695 if (bArb) rndm = G4UniformRand(); 1157 else rndm = eneRndm->GenRandEnergy(); << 696 else rndm = eneRndm->GenRandEnergy(); 1158 << 697 1159 G4double bracket = ((params.grad / 2.) << 698 G4double bracket = ((grad/2.)*(emaxsq - eminsq) + cept*(Emax-Emin)); 1160 * (emaxsq - eminsq) << 1161 + params.cept * (params.Em << 1162 bracket = bracket * rndm; 699 bracket = bracket * rndm; 1163 bracket = bracket + (params.grad / 2.) * em << 700 bracket = bracket + (grad/2.)*eminsq + cept*Emin; 1164 << 1165 // Now have a quad of form m/2 E**2 + cE - 701 // Now have a quad of form m/2 E**2 + cE - bracket = 0 1166 // << 1167 bracket = -bracket; 702 bracket = -bracket; 1168 << 703 // G4cout << "BRACKET" << bracket << G4endl; 1169 if (params.grad != 0.) << 704 if(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 { 705 { 1181 params.particle_energy = root1; << 706 G4double sqbrack = (intersq - 4*(grad/2.)*(bracket)); >> 707 // G4cout << "SQBRACK" << sqbrack << G4endl; >> 708 sqbrack = std::sqrt(sqbrack); >> 709 G4double root1 = -cept + sqbrack; >> 710 root1 = root1/(2.*(grad/2.)); >> 711 >> 712 G4double root2 = -cept - sqbrack; >> 713 root2 = root2/(2.*(grad/2.)); >> 714 >> 715 // G4cout << root1 << " roots " << root2 << G4endl; >> 716 >> 717 if(root1 > Emin && root1 < Emax) >> 718 particle_energy = root1; >> 719 if(root2 > Emin && root2 < Emax) >> 720 particle_energy = root2; 1182 } 721 } 1183 if (root2 > params.Emin && root2 < params << 722 else if(grad == 0.) 1184 { << 1185 params.particle_energy = root2; << 1186 } << 1187 } << 1188 else if (params.grad == 0.) << 1189 { << 1190 // have equation of form cE - bracket =0 723 // have equation of form cE - bracket =0 1191 // << 724 particle_energy = bracket/cept; 1192 params.particle_energy = bracket / params << 1193 } << 1194 << 1195 if (params.particle_energy < 0.) << 1196 { << 1197 params.particle_energy = -params.particle << 1198 } << 1199 725 1200 if (verbosityLevel >= 1) << 726 if(particle_energy < 0.) 1201 { << 727 particle_energy = -particle_energy; 1202 G4cout << "Energy is " << params.particle << 728 1203 } << 729 if(verbosityLevel >= 1) >> 730 G4cout << "Energy is " << particle_energy << G4endl; 1204 } 731 } 1205 732 1206 void G4SPSEneDistribution::GeneratePowEnergie 733 void G4SPSEneDistribution::GeneratePowEnergies(G4bool bArb = false) 1207 { 734 { 1208 // Method to generate particle energies dis << 735 // Method to generate particle energies distributed as >> 736 // a powerlaw 1209 737 1210 G4double rndm; 738 G4double rndm; 1211 G4double emina, emaxa; 739 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 740 1218 if (bArb) rndm = G4UniformRand(); << 741 emina = std::pow(Emin,alpha+1); 1219 else rndm = eneRndm->GenRandEnergy(); << 742 emaxa = std::pow(Emax,alpha+1); 1220 743 1221 if (params.alpha != -1.) << 744 if (bArb) rndm = G4UniformRand(); 1222 { << 745 else rndm = eneRndm->GenRandEnergy(); 1223 G4double ene = ((rndm * (emaxa - emina)) << 746 1224 ene = std::pow(ene, (1. / (params.alpha + << 747 if(alpha != -1.) 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 { 748 { 1275 nabove = middle; << 749 particle_energy = ((rndm*(emaxa - emina)) + emina); >> 750 particle_energy = std::pow(particle_energy,(1./(alpha+1.))); 1276 } 751 } 1277 else << 752 else if(alpha == -1.) 1278 { 753 { 1279 nbelow = middle; << 754 particle_energy = (std::log(Emin) + rndm*(std::log(Emax) - std::log(Emin))); >> 755 particle_energy = std::exp(particle_energy); 1280 } 756 } 1281 } << 757 if(verbosityLevel >= 1) 1282 << 758 G4cout << "Energy is " << particle_energy << G4endl; 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 } 759 } 1355 760 1356 void G4SPSEneDistribution::GenerateExpEnergie 761 void G4SPSEneDistribution::GenerateExpEnergies(G4bool bArb = false) 1357 { 762 { 1358 // Method to generate particle energies dis 763 // Method to generate particle energies distributed according 1359 // to an exponential curve << 764 // to an exponential curve. 1360 << 1361 G4double rndm; 765 G4double rndm; 1362 << 766 1363 if (bArb) rndm = G4UniformRand(); 767 if (bArb) rndm = G4UniformRand(); 1364 else rndm = eneRndm->GenRandEnergy(); << 768 else rndm = eneRndm->GenRandEnergy(); 1365 769 1366 threadLocal_t& params = threadLocalData.Get << 770 particle_energy = -Ezero*(std::log(rndm*(std::exp(-Emax/Ezero) - std::exp(-Emin/Ezero)) + 1367 params.particle_energy = -params.Ezero << 771 std::exp(-Emin/Ezero))); 1368 * (std::log(rndm * ( << 772 if(verbosityLevel >= 1) 1369 << 773 G4cout << "Energy is " << particle_energy << G4endl; 1370 - << 1371 << 1372 + std::exp << 1373 if (verbosityLevel >= 1) << 1374 { << 1375 G4cout << "Energy is " << params.particle << 1376 } << 1377 } 774 } 1378 775 1379 void G4SPSEneDistribution::GenerateBremEnergi 776 void G4SPSEneDistribution::GenerateBremEnergies() 1380 { 777 { 1381 // Method to generate particle energies dis << 778 // Method to generate particle energies distributed according 1382 // to a Bremstrahlung equation of the form << 779 // to a Bremstrahlung equation of 1383 // I = const*((kT)**1/2)*E*(e**(-E/kT)) << 780 // form I = const*((kT)**1/2)*E*(e**(-E/kT)) 1384 << 781 1385 G4double rndm = eneRndm->GenRandEnergy(); << 782 G4double rndm; >> 783 rndm = eneRndm->GenRandEnergy(); 1386 G4double expmax, expmin, k; 784 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 785 1392 threadLocal_t& params = threadLocalData.Get << 786 k = 8.6181e-11; // Boltzmann's const in MeV/K >> 787 G4double ksq = std::pow(k,2.); // k squared >> 788 G4double Tsq = std::pow(Temp,2.); // Temp squared 1393 789 1394 expmax = std::exp(-params.Emax / (k * Temp) << 790 expmax = std::exp(-Emax/(k*Temp)); 1395 expmin = std::exp(-params.Emin / (k * Temp) << 791 expmin = std::exp(-Emin/(k*Temp)); 1396 792 1397 // If either expmax or expmin are zero then 793 // If either expmax or expmin are zero then this will cause problems 1398 // Most probably this will be because T is 794 // Most probably this will be because T is too low or E is too high 1399 795 1400 if (expmax == 0.) << 796 if(expmax == 0.) 1401 { << 797 G4cout << "*****EXPMAX=0. Choose different E's or Temp" << G4endl; 1402 G4Exception("G4SPSEneDistribution::Genera << 798 if(expmin == 0.) 1403 "Event0302", FatalException, << 799 G4cout << "*****EXPMIN=0. Choose different E's or Temp" << G4endl; 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 800 1413 G4double tempvar = rndm * ((-k) * Temp * (p << 801 G4double tempvar = rndm *((-k)*Temp*(Emax*expmax - Emin*expmin) - 1414 - p << 802 (ksq*Tsq*(expmax-expmin))); 1415 - (ksq * Tsq * (exp << 1416 803 1417 G4double bigc = (tempvar - k * Temp * param << 804 G4double bigc = (tempvar - k*Temp*Emin*expmin - ksq*Tsq*expmin)/(-k*Temp); 1418 - ksq * Tsq * expmin) / (-k << 1419 805 1420 // This gives an equation of form: Ee(-E/kT 806 // This gives an equation of form: Ee(-E/kT) + kTe(-E/kT) - C =0 1421 // Solve this iteratively, step from Emin t 807 // Solve this iteratively, step from Emin to Emax in 1000 steps 1422 // and take the best solution. 808 // and take the best solution. 1423 809 1424 G4double erange = params.Emax - params.Emin << 810 G4double erange = Emax - Emin; 1425 G4double steps = erange / 1000.; << 811 G4double steps = erange/1000.; 1426 G4int i; 812 G4int i; 1427 G4double etest, diff, err = 100000.; << 813 G4double etest, diff, err; 1428 << 814 1429 for (i = 1; i < 1000; ++i) << 815 err = 100000.; 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 816 1440 if (diff < err) << 817 for(i=1; i<1000; i++) 1441 { 818 { 1442 err = diff; << 819 etest = Emin + (i-1)*steps; 1443 params.particle_energy = etest; << 820 1444 } << 821 diff = etest*(std::exp(-etest/(k*Temp))) + k*Temp*(std::exp(-etest/(k*Temp))) - bigc; 1445 } << 822 1446 if (verbosityLevel >= 1) << 823 if(diff < 0.) 1447 { << 824 diff = -diff; 1448 G4cout << "Energy is " << params.particle << 825 1449 } << 826 if(diff < err) >> 827 { >> 828 err = diff; >> 829 particle_energy = etest; >> 830 } >> 831 } >> 832 if(verbosityLevel >= 1) >> 833 G4cout << "Energy is " << particle_energy << G4endl; 1450 } 834 } 1451 835 1452 void G4SPSEneDistribution::GenerateBbodyEnerg 836 void G4SPSEneDistribution::GenerateBbodyEnergies() 1453 { 837 { 1454 // BBody_x holds Energies, and BBHist holds 838 // BBody_x holds Energies, and BBHist holds the cumulative histo. 1455 // Binary search to find correct bin then l << 839 // binary search to find correct bin then lin interpolation. 1456 // Use the earlier defined histogram + Rand 840 // Use the earlier defined histogram + RandGeneral method to generate 1457 // random numbers following the histos dist << 841 // random numbers following the histos distribution. 1458 << 842 G4double rndm; 1459 G4double rndm = eneRndm->GenRandEnergy(); << 843 G4int nabove, nbelow = 0, middle; 1460 G4int nabove = 10001, nbelow = 0, middle; << 844 nabove = 10001; 1461 << 845 rndm = eneRndm->GenRandEnergy(); 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 846 1474 // Binary search to find bin that rndm is i 847 // Binary search to find bin that rndm is in 1475 // << 848 while(nabove-nbelow > 1) 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 { 849 { 1485 nabove = middle; << 850 middle = (nabove + nbelow)/2; >> 851 if(rndm == BBHist[middle]) break; >> 852 if(rndm < BBHist[middle]) nabove = middle; >> 853 else nbelow = middle; 1486 } 854 } 1487 else << 855 >> 856 // Now interpolate in that bin to find the correct output value. >> 857 G4double x1, x2, y1, y2, m, q; >> 858 x1 = Bbody_x[nbelow]; >> 859 x2 = Bbody_x[nbelow+1]; >> 860 y1 = BBHist[nbelow]; >> 861 y2 = BBHist[nbelow+1]; >> 862 m = (y2-y1)/(x2-x1); >> 863 q = y1 - m*x1; >> 864 >> 865 particle_energy = (rndm - q)/m; >> 866 >> 867 if(verbosityLevel >= 1) 1488 { 868 { 1489 nbelow = middle; << 869 G4cout << "Energy is " << particle_energy << G4endl; 1490 } 870 } 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 } 871 } 1526 872 1527 void G4SPSEneDistribution::GenerateCdgEnergie 873 void G4SPSEneDistribution::GenerateCdgEnergies() 1528 { 874 { 1529 // Generate random numbers, compare with va << 875 // Gen random numbers, compare with values in cumhist 1530 // to find appropriate part of spectrum and << 876 // to find appropriate part of spectrum and then 1531 // generate energy in the usual inversion w << 877 // generate energy in the usual inversion way. 1532 << 878 // G4double pfact[2] = {8.5, 112}; >> 879 // G4double spind[2] = {1.4, 2.3}; >> 880 // G4double ene_line[3] = {1., 18., 1E6}; 1533 G4double rndm, rndm2; 881 G4double rndm, rndm2; 1534 G4double ene_line[3]={0,0,0}; << 882 G4double ene_line[3]; 1535 G4double omalpha[2]={0,0}; << 883 G4double omalpha[2]; 1536 threadLocal_t& params = threadLocalData.Get << 884 if(Emin < 18*keV && Emax < 18*keV) 1537 if (params.Emin < 18 * keV && params.Emax < << 885 { 1538 { << 886 omalpha[0] = 1. - 1.4; 1539 omalpha[0] = 1. - 1.4; << 887 ene_line[0] = Emin; 1540 ene_line[0] = params.Emin; << 888 ene_line[1] = Emax; 1541 ene_line[1] = params.Emax; << 889 } 1542 } << 890 if(Emin < 18*keV && Emax > 18*keV) 1543 if (params.Emin < 18 * keV && params.Emax > << 891 { 1544 { << 892 omalpha[0] = 1. - 1.4; 1545 omalpha[0] = 1. - 1.4; << 893 omalpha[1] = 1. - 2.3; 1546 omalpha[1] = 1. - 2.3; << 894 ene_line[0] = Emin; 1547 ene_line[0] = params.Emin; << 895 ene_line[1] = 18.; 1548 ene_line[1] = 18. * keV; << 896 ene_line[2] = Emax; 1549 ene_line[2] = params.Emax; << 897 } 1550 } << 898 if(Emin > 18*keV) 1551 if (params.Emin > 18 * keV) << 899 { 1552 { << 900 omalpha[0] = 1. - 2.3; 1553 omalpha[0] = 1. - 2.3; << 901 ene_line[0] = Emin; 1554 ene_line[0] = params.Emin; << 902 ene_line[1] = Emax; 1555 ene_line[1] = params.Emax; << 903 } 1556 } << 1557 rndm = eneRndm->GenRandEnergy(); 904 rndm = eneRndm->GenRandEnergy(); 1558 rndm2 = eneRndm->GenRandEnergy(); 905 rndm2 = eneRndm->GenRandEnergy(); 1559 906 1560 G4int i = 0; 907 G4int i = 0; 1561 while (rndm >= CDGhist[i]) << 908 while( rndm >= CDGhist[i]) 1562 { << 909 { 1563 ++i; << 910 i++; 1564 } << 911 } >> 912 // Generate final energy. >> 913 particle_energy = (std::pow(ene_line[i-1],omalpha[i-1]) + (std::pow(ene_line[i],omalpha[i-1]) >> 914 - std::pow(ene_line[i-1],omalpha[i-1]))*rndm2); >> 915 particle_energy = std::pow(particle_energy,(1./omalpha[i-1])); 1565 916 1566 // Generate final energy << 917 if(verbosityLevel >= 1) 1567 // << 918 G4cout << "Energy is " << particle_energy << G4endl; 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 } 919 } 1578 920 1579 void G4SPSEneDistribution::GenUserHistEnergie 921 void G4SPSEneDistribution::GenUserHistEnergies() 1580 { 922 { 1581 // Histograms are DIFFERENTIAL << 923 // Histograms are DIFFERENTIAL. 1582 << 924 // G4cout << "In GenUserHistEnergies " << G4endl; 1583 G4AutoLock l(&mutex); << 925 if(IPDFEnergyExist == false) 1584 << 926 { 1585 if (!IPDFEnergyExist) << 927 G4int ii; 1586 { << 928 G4int maxbin = G4int(UDefEnergyH.GetVectorLength()); 1587 std::size_t ii; << 929 G4double bins[1024], vals[1024], sum; 1588 std::size_t maxbin = UDefEnergyH.GetVecto << 930 sum=0.; 1589 G4double bins[1024], vals[1024], sum; << 931 1590 for ( ii = 0 ; ii<1024 ; ++ii ) { bins[ii << 932 if((EnergySpec == false) && (particle_definition == NULL)) 1591 sum = 0.; << 933 G4cout << "Error: particle definition is NULL" << G4endl; 1592 << 934 1593 if ( (!EnergySpec) << 935 if(maxbin > 1024) 1594 && (threadLocalData.Get().particle_defi << 936 { 1595 { << 937 G4cout << "Maxbin > 1024" << G4endl; 1596 G4Exception("G4SPSEneDistribution::GenU << 938 G4cout << "Setting maxbin to 1024, other bins are lost" << G4endl; 1597 "Event0302", FatalException << 939 } 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 940 1660 IPDFEnergyExist = true; << 941 if(DiffSpec == false) 1661 if (verbosityLevel > 1) << 942 G4cout << "Histograms are Differential!!! " << G4endl; 1662 { << 943 else 1663 IPDFEnergyH.DumpValues(); << 944 { >> 945 bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0)); >> 946 vals[0] = UDefEnergyH(size_t(0)); >> 947 sum = vals[0]; >> 948 for(ii=1;ii<maxbin;ii++) >> 949 { >> 950 bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii)); >> 951 vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii-1]; >> 952 sum = sum + UDefEnergyH(size_t(ii)); >> 953 } >> 954 } >> 955 >> 956 if(EnergySpec == false) >> 957 { >> 958 G4double mass = particle_definition->GetPDGMass(); >> 959 // multiply the function (vals) up by the bin width >> 960 // to make the function counts/s (i.e. get rid of momentum >> 961 // dependence). >> 962 for(ii=1;ii<maxbin;ii++) >> 963 { >> 964 vals[ii] = vals[ii] * (bins[ii] - bins[ii-1]); >> 965 } >> 966 // Put energy bins into new histo, plus divide by energy bin width >> 967 // to make evals counts/s/energy >> 968 for(ii=0;ii<maxbin;ii++) >> 969 { >> 970 bins[ii] = std::sqrt((bins[ii]*bins[ii]) + (mass*mass)) - mass; //kinetic energy >> 971 } >> 972 for(ii=1;ii<maxbin;ii++) >> 973 { >> 974 vals[ii] = vals[ii]/(bins[ii] - bins[ii-1]); >> 975 } >> 976 sum = vals[maxbin-1]; >> 977 vals[0] = 0.; >> 978 } >> 979 for(ii=0;ii<maxbin;ii++) >> 980 { >> 981 vals[ii] = vals[ii]/sum; >> 982 IPDFEnergyH.InsertValues(bins[ii], vals[ii]); >> 983 } >> 984 >> 985 // Make IPDFEnergyExist = true >> 986 IPDFEnergyExist = true; >> 987 if(verbosityLevel > 1) >> 988 IPDFEnergyH.DumpValues(); 1664 } 989 } 1665 } << 990 1666 l.unlock(); << 1667 << 1668 // IPDF has been create so carry on 991 // IPDF has been create so carry on 1669 // << 1670 G4double rndm = eneRndm->GenRandEnergy(); 992 G4double rndm = eneRndm->GenRandEnergy(); 1671 threadLocalData.Get().particle_energy= IPDF << 993 particle_energy = IPDFEnergyH.GetEnergy(rndm); 1672 << 994 1673 if (verbosityLevel >= 1) << 995 if(verbosityLevel >= 1) 1674 { << 1675 G4cout << "Energy is " << particle_energy 996 G4cout << "Energy is " << particle_energy << G4endl; 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 } 997 } 1708 998 1709 void G4SPSEneDistribution::GenArbPointEnergie 999 void G4SPSEneDistribution::GenArbPointEnergies() 1710 { 1000 { 1711 if (verbosityLevel > 0) << 1001 if(verbosityLevel > 0) 1712 { << 1713 G4cout << "In GenArbPointEnergies" << G4e 1002 G4cout << "In GenArbPointEnergies" << G4endl; 1714 } << 1003 G4double rndm; 1715 << 1004 rndm = eneRndm->GenRandEnergy(); 1716 G4double rndm = eneRndm->GenRandEnergy(); << 1005 if(IntType != "Spline") 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 { 1006 { 1780 G4cout << "Energy is " << params.partic << 1007 // IPDFArbEnergyH.DumpValues(); >> 1008 // Find the Bin >> 1009 // have x, y, no of points, and cumulative area distribution >> 1010 G4int nabove, nbelow = 0, middle; >> 1011 nabove = IPDFArbEnergyH.GetVectorLength(); >> 1012 // G4cout << nabove << G4endl; >> 1013 // Binary search to find bin that rndm is in >> 1014 while(nabove-nbelow > 1) >> 1015 { >> 1016 middle = (nabove + nbelow)/2; >> 1017 if(rndm == IPDFArbEnergyH(size_t(middle))) break; >> 1018 if(rndm < IPDFArbEnergyH(size_t(middle))) nabove = middle; >> 1019 else nbelow = middle; >> 1020 } >> 1021 if(IntType == "Lin") >> 1022 { >> 1023 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1)); >> 1024 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); >> 1025 grad = Arb_grad[nbelow+1]; >> 1026 cept = Arb_cept[nbelow+1]; >> 1027 // G4cout << rndm << " " << Emax << " " << Emin << " " << grad << " " << cept << G4endl; >> 1028 GenerateLinearEnergies(true); >> 1029 } >> 1030 else if(IntType == "Log") >> 1031 { >> 1032 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1)); >> 1033 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); >> 1034 alpha = Arb_alpha[nbelow+1]; >> 1035 // G4cout << rndm << " " << Emax << " " << Emin << " " << alpha << G4endl; >> 1036 GeneratePowEnergies(true); >> 1037 } >> 1038 else if(IntType == "Exp") >> 1039 { >> 1040 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1)); >> 1041 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); >> 1042 Ezero = Arb_ezero[nbelow+1]; >> 1043 // G4cout << rndm << " " << Emax << " " << Emin << " " << Ezero << G4endl; >> 1044 GenerateExpEnergies(true); >> 1045 } >> 1046 } >> 1047 else if(IntType == "Spline") >> 1048 { >> 1049 if(verbosityLevel > 1) >> 1050 G4cout << "IntType = Spline " << rndm << G4endl; >> 1051 // in SplineInterpolation created SplineInt >> 1052 // Now generate a random number put it into CubicSplineInterpolation >> 1053 // and you should get out an energy!?! >> 1054 particle_energy = -1e100; >> 1055 while (particle_energy < Emin || particle_energy > Emax ) { >> 1056 particle_energy = SplineInt->CubicSplineInterpolation(rndm); >> 1057 rndm = eneRndm->GenRandEnergy(); >> 1058 } >> 1059 if(verbosityLevel >= 1) >> 1060 G4cout << "Energy is " << particle_energy << G4endl; 1781 } 1061 } 1782 } << 1783 else 1062 else 1784 { << 1063 G4cout << "Error: IntType unknown type" << G4endl; 1785 G4Exception("G4SPSEneDistribution::GenArb << 1786 FatalException, "Error: IntTy << 1787 } << 1788 } 1064 } 1789 1065 1790 void G4SPSEneDistribution::GenEpnHistEnergies 1066 void G4SPSEneDistribution::GenEpnHistEnergies() 1791 { 1067 { 1792 // Firstly convert to energy if not already << 1068 // G4cout << "In GenEpnHistEnergies " << Epnflag << G4endl; 1793 << 1069 1794 G4AutoLock l(&mutex); << 1070 // Firstly convert to energy if not already done. 1795 << 1071 if(Epnflag == true) 1796 if (Epnflag) // true means spectrum is epn << 1072 // epnflag = true means spectrum is epn, false means e. 1797 { << 1073 { 1798 // Convert to energy by multiplying by A << 1074 // convert to energy by multiplying by A number 1799 // << 1075 ConvertEPNToEnergy(); 1800 ConvertEPNToEnergy(); << 1076 // EpnEnergyH will be replace by UDefEnergyH. 1801 } << 1077 // UDefEnergyH.DumpValues(); 1802 if (!IPDFEnergyExist) << 1078 } 1803 { << 1079 1804 // IPDF has not been created, so create i << 1080 // G4cout << "Creating IPDFEnergy if not already done so" << G4endl; 1805 // << 1081 if(IPDFEnergyExist == false) 1806 G4double bins[1024], vals[1024], sum; << 1082 { 1807 std::size_t ii; << 1083 // IPDF has not been created, so create it 1808 std::size_t maxbin = UDefEnergyH.GetVecto << 1084 G4double bins[1024],vals[1024], sum; 1809 bins[0] = UDefEnergyH.GetLowEdgeEnergy(0) << 1085 G4int ii; 1810 vals[0] = UDefEnergyH(0); << 1086 G4int maxbin = G4int(UDefEnergyH.GetVectorLength()); 1811 sum = vals[0]; << 1087 bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0)); 1812 for (ii = 1; ii < maxbin; ++ii) << 1088 vals[0] = UDefEnergyH(size_t(0)); 1813 { << 1089 sum = vals[0]; 1814 bins[ii] = UDefEnergyH.GetLowEdgeEnergy << 1090 for(ii=1;ii<maxbin;ii++) 1815 vals[ii] = UDefEnergyH(ii) + vals[ii - << 1091 { 1816 sum = sum + UDefEnergyH(ii); << 1092 bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii)); 1817 } << 1093 vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii-1]; 1818 << 1094 sum = sum + UDefEnergyH(size_t(ii)); 1819 l.lock(); << 1095 } 1820 for (ii = 0; ii < maxbin; ++ii) << 1096 1821 { << 1097 for(ii=0;ii<maxbin;ii++) 1822 vals[ii] = vals[ii] / sum; << 1098 { 1823 IPDFEnergyH.InsertValues(bins[ii], vals << 1099 vals[ii] = vals[ii]/sum; >> 1100 IPDFEnergyH.InsertValues(bins[ii], vals[ii]); >> 1101 } >> 1102 // Make IPDFEpnExist = true >> 1103 IPDFEnergyExist = true; 1824 } 1104 } 1825 IPDFEnergyExist = true; << 1105 // IPDFEnergyH.DumpValues(); 1826 << 1827 } << 1828 l.unlock(); << 1829 << 1830 // IPDF has been create so carry on 1106 // IPDF has been create so carry on 1831 // << 1832 G4double rndm = eneRndm->GenRandEnergy(); 1107 G4double rndm = eneRndm->GenRandEnergy(); 1833 threadLocalData.Get().particle_energy = IPD << 1108 particle_energy = IPDFEnergyH.GetEnergy(rndm); 1834 1109 1835 if (verbosityLevel >= 1) << 1110 if(verbosityLevel >= 1) 1836 { << 1111 G4cout << "Energy is " << particle_energy << G4endl; 1837 G4cout << "Energy is " << threadLocalData << 1838 } << 1839 } 1112 } 1840 1113 1841 void G4SPSEneDistribution::ConvertEPNToEnergy << 1114 void G4SPSEneDistribution::ConvertEPNToEnergy() 1842 { 1115 { 1843 // Use this before particle generation to c << 1116 // Use this before particle generation to convert the 1844 // currently stored histogram from energy/n << 1117 // currently stored histogram from energy/nucleon 1845 << 1118 // to energy. 1846 threadLocal_t& params = threadLocalData.Get << 1119 // G4cout << "In ConvertEpntoEnergy " << G4endl; 1847 if (params.particle_definition == nullptr) << 1120 if(particle_definition==NULL) 1848 { << 1849 G4cout << "Error: particle not defined" < 1121 G4cout << "Error: particle not defined" << G4endl; 1850 } << 1851 else 1122 else 1852 { << 1853 // Need to multiply histogram by the numb << 1854 // Baryon Number looks to hold the no. of << 1855 // << 1856 G4int Bary = params.particle_definition-> << 1857 << 1858 // Change values in histogram, Read it ou << 1859 // << 1860 std::size_t count, maxcount; << 1861 maxcount = EpnEnergyH.GetVectorLength(); << 1862 G4double ebins[1024], evals[1024]; << 1863 if (maxcount > 1024) << 1864 { << 1865 G4Exception("G4SPSEneDistribution::Conv << 1866 "gps001", JustWarning, << 1867 "Histogram contains more th << 1868 Those above 1024 will be i << 1869 maxcount = 1024; << 1870 } << 1871 if (maxcount < 1) << 1872 { 1123 { 1873 G4Exception("G4SPSEneDistribution::Conv << 1124 // Need to multiply histogram by the number of nucleons. 1874 "gps001", FatalException, << 1125 // Baryon Number looks to hold the no. of nucleons. 1875 "Histogram contains less tha << 1126 G4int Bary = particle_definition->GetBaryonNumber(); 1876 return; << 1127 // G4cout << "Baryon No. " << Bary << G4endl; 1877 } << 1128 // Change values in histogram, Read it out, delete it, re-create it 1878 for (count = 0; count < maxcount; ++count << 1129 G4int count, maxcount; 1879 { << 1130 maxcount = G4int(EpnEnergyH.GetVectorLength()); 1880 // Read out << 1131 // G4cout << maxcount << G4endl; 1881 ebins[count] = EpnEnergyH.GetLowEdgeEne << 1132 G4double ebins[1024],evals[1024]; 1882 evals[count] = EpnEnergyH(count); << 1133 if(maxcount > 1024) 1883 } << 1134 { 1884 << 1135 G4cout << "Histogram contains more than 1024 bins!" << G4endl; 1885 // Multiply the channels by the nucleon n << 1136 G4cout << "Those above 1024 will be ignored" << G4endl; 1886 // << 1137 maxcount = 1024; 1887 for (count = 0; count < maxcount; ++count << 1138 } 1888 { << 1139 for(count=0;count<maxcount;count++) 1889 ebins[count] = ebins[count] * Bary; << 1140 { 1890 } << 1141 // Read out 1891 << 1142 ebins[count] = EpnEnergyH.GetLowEdgeEnergy(size_t(count)); 1892 // Set Emin and Emax << 1143 evals[count] = EpnEnergyH(size_t(count)); 1893 // << 1144 } 1894 params.Emin = ebins[0]; << 1145 1895 if (maxcount > 1) << 1146 // Multiply the channels by the nucleon number to give energies 1896 { << 1147 for(count=0;count<maxcount;count++) 1897 params.Emax = ebins[maxcount - 1]; << 1148 { 1898 } << 1149 ebins[count] = ebins[count] * Bary; 1899 else << 1150 } 1900 { << 1151 1901 params.Emax = ebins[0]; << 1152 // Set Emin and Emax 1902 } << 1153 Emin = ebins[0]; 1903 << 1154 Emax = ebins[maxcount-1]; 1904 // Put energy bins into new histogram - U << 1155 1905 // << 1156 // Put energy bins into new histogram - UDefEnergyH. 1906 for (count = 0; count < maxcount; ++count << 1157 for(count=0;count<maxcount;count++) 1907 { << 1158 { 1908 UDefEnergyH.InsertValues(ebins[count], << 1159 UDefEnergyH.InsertValues(ebins[count], evals[count]); 1909 } << 1160 } 1910 Epnflag = false; // so that you dont repe << 1161 Epnflag = false; //so that you dont repeat this method. 1911 } << 1162 } 1912 } 1163 } 1913 1164 1914 void G4SPSEneDistribution::ReSetHist(const G4 << 1165 // >> 1166 void G4SPSEneDistribution::ReSetHist(G4String atype) 1915 { 1167 { 1916 G4AutoLock l(&mutex); << 1168 if (atype == "energy"){ 1917 if (atype == "energy") << 1169 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ; 1918 { << 1170 IPDFEnergyExist = false ; 1919 UDefEnergyH = IPDFEnergyH = ZeroPhysVecto << 1920 IPDFEnergyExist = false; << 1921 Emin = 0.; 1171 Emin = 0.; 1922 Emax = 1e30; << 1172 Emax = 1e30;} 1923 } << 1173 else if ( atype == "arb"){ 1924 else if (atype == "arb") << 1174 ArbEnergyH =IPDFArbEnergyH = ZeroPhysVector ; 1925 { << 1175 IPDFArbExist = false;} 1926 ArbEnergyH = IPDFArbEnergyH = ZeroPhysVec << 1176 else if ( atype == "epn"){ 1927 IPDFArbExist = false; << 1177 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ; 1928 } << 1178 IPDFEnergyExist = false ; 1929 else if (atype == "epn") << 1179 EpnEnergyH = ZeroPhysVector ;} 1930 { << 1180 else { 1931 UDefEnergyH = IPDFEnergyH = ZeroPhysVecto << 1932 IPDFEnergyExist = false; << 1933 EpnEnergyH = ZeroPhysVector; << 1934 } << 1935 else << 1936 { << 1937 G4cout << "Error, histtype not accepted " 1181 G4cout << "Error, histtype not accepted " << G4endl; 1938 } 1182 } 1939 } 1183 } 1940 1184 1941 G4double G4SPSEneDistribution::GenerateOne(G4 1185 G4double G4SPSEneDistribution::GenerateOne(G4ParticleDefinition* a) 1942 { 1186 { 1943 // Copy global shared status to thread-loca << 1187 particle_definition = a; 1944 // << 1188 particle_energy = -1.; 1945 threadLocal_t& params = threadLocalData.Get << 1189 while ( (EnergyDisType == "Arb")? (particle_energy < ArbEmin || particle_energy > ArbEmax) 1946 params.particle_definition=a; << 1190 : (particle_energy < Emin || particle_energy > Emax) ) { 1947 params.particle_energy=-1; << 1191 if(EnergyDisType == "Mono") 1948 if(applyEvergyWeight) << 1192 GenerateMonoEnergetic(); 1949 { << 1193 else if(EnergyDisType == "Lin") 1950 params.Emax = ArbEmax; << 1194 GenerateLinearEnergies(); 1951 params.Emin = ArbEmin; << 1195 else if(EnergyDisType == "Pow") 1952 } << 1196 GeneratePowEnergies(); 1953 else << 1197 else if(EnergyDisType == "Exp") 1954 { << 1198 GenerateExpEnergies(); 1955 params.Emax = Emax; << 1199 else if(EnergyDisType == "Gauss") 1956 params.Emin = Emin; << 1200 GenerateGaussEnergies(); 1957 } << 1201 else if(EnergyDisType == "Brem") 1958 params.alpha = alpha; << 1202 GenerateBremEnergies(); 1959 params.Ezero = Ezero; << 1203 else if(EnergyDisType == "Bbody") 1960 params.grad = grad; << 1204 GenerateBbodyEnergies(); 1961 params.cept = cept; << 1205 else if(EnergyDisType == "Cdg") 1962 params.weight = weight; << 1206 GenerateCdgEnergies(); 1963 // particle_energy = -1.; << 1207 else if(EnergyDisType == "User") 1964 << 1208 GenUserHistEnergies(); 1965 if((EnergyDisType == "Mono") && ((MonoEnerg << 1209 else if(EnergyDisType == "Arb") 1966 { << 1210 GenArbPointEnergies(); 1967 G4ExceptionDescription ed; << 1211 else if(EnergyDisType == "Epn") 1968 ed << "MonoEnergy " << G4BestUnit(MonoEne << 1212 GenEpnHistEnergies(); 1969 << " is outside of [Emin,Emax] = [" << 1970 << G4BestUnit(Emin,"Energy") << ", " << 1971 << G4BestUnit(Emax,"Energy") << ". Mon << 1972 G4Exception("G4SPSEneDistribution::Genera << 1973 "GPS0001", JustWarning, ed); << 1974 params.particle_energy=MonoEnergy; << 1975 return params.particle_energy; << 1976 } << 1977 while ( (EnergyDisType == "Arb") << 1978 ? (params.particle_energy < ArbEmin << 1979 || params.particle_energy > ArbEmax << 1980 : (params.particle_energy < params. << 1981 || params.particle_energy > params. << 1982 { << 1983 if (Biased) << 1984 { << 1985 GenerateBiasPowEnergies(); << 1986 } << 1987 else 1213 else 1988 { << 1214 G4cout << "Error: EnergyDisType has unusual value" << G4endl; 1989 if (EnergyDisType == "Mono") << 1990 { << 1991 GenerateMonoEnergetic(); << 1992 } << 1993 else if (EnergyDisType == "Lin") << 1994 { << 1995 GenerateLinearEnergies(false); << 1996 } << 1997 else if (EnergyDisType == "Pow") << 1998 { << 1999 GeneratePowEnergies(false); << 2000 } << 2001 else if (EnergyDisType == "CPow") << 2002 { << 2003 GenerateCPowEnergies(); << 2004 } << 2005 else if (EnergyDisType == "Exp") << 2006 { << 2007 GenerateExpEnergies(false); << 2008 } << 2009 else if (EnergyDisType == "Gauss") << 2010 { << 2011 GenerateGaussEnergies(); << 2012 } << 2013 else if (EnergyDisType == "Brem") << 2014 { << 2015 GenerateBremEnergies(); << 2016 } << 2017 else if (EnergyDisType == "Bbody") << 2018 { << 2019 GenerateBbodyEnergies(); << 2020 } << 2021 else if (EnergyDisType == "Cdg") << 2022 { << 2023 GenerateCdgEnergies(); << 2024 } << 2025 else if (EnergyDisType == "User") << 2026 { << 2027 GenUserHistEnergies(); << 2028 } << 2029 else if (EnergyDisType == "Arb") << 2030 { << 2031 GenArbPointEnergies(); << 2032 } << 2033 else if (EnergyDisType == "Epn") << 2034 { << 2035 GenEpnHistEnergies(); << 2036 } << 2037 else << 2038 { << 2039 G4cout << "Error: EnergyDisType has u << 2040 } << 2041 } << 2042 } 1215 } 2043 return params.particle_energy; << 1216 return particle_energy; 2044 } 1217 } 2045 1218 2046 G4double G4SPSEneDistribution::GetProbability << 2047 { << 2048 G4double prob = 1.; << 2049 1219 2050 threadLocal_t& params = threadLocalData.Get << 2051 if (EnergyDisType == "Lin") << 2052 { << 2053 if (prob_norm == 1.) << 2054 { << 2055 prob_norm = 0.5*params.grad*params.Emax << 2056 + params.cept*params.Emax << 2057 - 0.5*params.grad*params.Emin << 2058 - params.cept*params.Emin; << 2059 } << 2060 prob = params.cept + params.grad * ene; << 2061 prob /= prob_norm; << 2062 } << 2063 else if (EnergyDisType == "Pow") << 2064 { << 2065 if (prob_norm == 1.) << 2066 { << 2067 if (alpha != -1.) << 2068 { << 2069 G4double emina = std::pow(params.Emin << 2070 G4double emaxa = std::pow(params.Emax << 2071 prob_norm = 1./(1.+alpha) * (emaxa - << 2072 } << 2073 else << 2074 { << 2075 prob_norm = std::log(params.Emax) - s << 2076 } << 2077 } << 2078 prob = std::pow(ene, params.alpha)/prob_n << 2079 } << 2080 else if (EnergyDisType == "Exp") << 2081 { << 2082 if (prob_norm == 1.) << 2083 { << 2084 prob_norm = -params.Ezero*(std::exp(-pa << 2085 - std::exp(par << 2086 } << 2087 prob = std::exp(-ene / params.Ezero); << 2088 prob /= prob_norm; << 2089 } << 2090 else if (EnergyDisType == "Arb") << 2091 { << 2092 prob = ArbEnergyH.Value(ene); << 2093 << 2094 if (prob <= 0.) << 2095 { << 2096 G4cout << " Warning:G4SPSEneDistributio << 2097 << prob << " " << ene << G4endl; << 2098 prob = 1e-30; << 2099 } << 2100 } << 2101 else << 2102 { << 2103 G4cout << "Error: EnergyDisType not suppo << 2104 } << 2105 1220 2106 return prob; << 1221 2107 } << 1222 >> 1223 >> 1224 >> 1225 >> 1226 2108 1227