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