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 // G4SPSAngDistribution class implementation 26 // G4SPSAngDistribution 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: Andrea Dotti, SLAC 30 // Revisions: Andrea Dotti, SLAC 31 // ------------------------------------------- 31 // -------------------------------------------------------------------- 32 32 33 #include "G4SPSAngDistribution.hh" 33 #include "G4SPSAngDistribution.hh" 34 34 35 #include "Randomize.hh" 35 #include "Randomize.hh" 36 #include "G4PhysicalConstants.hh" 36 #include "G4PhysicalConstants.hh" 37 37 38 G4SPSAngDistribution::G4SPSAngDistribution() 38 G4SPSAngDistribution::G4SPSAngDistribution() 39 { 39 { 40 // Angular distribution Variables 40 // Angular distribution Variables 41 G4ThreeVector zero; 41 G4ThreeVector zero; 42 particle_momentum_direction = G4ParticleMome 42 particle_momentum_direction = G4ParticleMomentum(0,0,-1); 43 43 44 AngDistType = "planar"; 44 AngDistType = "planar"; 45 AngRef1 = CLHEP::HepXHat; 45 AngRef1 = CLHEP::HepXHat; 46 AngRef2 = CLHEP::HepYHat; 46 AngRef2 = CLHEP::HepYHat; 47 AngRef3 = CLHEP::HepZHat; 47 AngRef3 = CLHEP::HepZHat; 48 MinTheta = 0.; 48 MinTheta = 0.; 49 MaxTheta = pi; 49 MaxTheta = pi; 50 MinPhi = 0.; 50 MinPhi = 0.; 51 MaxPhi = twopi; 51 MaxPhi = twopi; 52 DR = 0.; 52 DR = 0.; 53 DX = 0.; 53 DX = 0.; 54 DY = 0.; 54 DY = 0.; 55 FocusPoint = G4ThreeVector(0., 0., 0.); 55 FocusPoint = G4ThreeVector(0., 0., 0.); 56 UserDistType = "NULL"; 56 UserDistType = "NULL"; 57 UserWRTSurface = true; 57 UserWRTSurface = true; 58 UserAngRef = false; 58 UserAngRef = false; 59 IPDFThetaExist = false; 59 IPDFThetaExist = false; 60 IPDFPhiExist = false; 60 IPDFPhiExist = false; 61 verbosityLevel = 0; 61 verbosityLevel = 0; 62 62 63 G4MUTEXINIT(mutex); 63 G4MUTEXINIT(mutex); 64 } 64 } 65 65 66 G4SPSAngDistribution::~G4SPSAngDistribution() 66 G4SPSAngDistribution::~G4SPSAngDistribution() 67 { 67 { 68 G4MUTEXDESTROY(mutex); 68 G4MUTEXDESTROY(mutex); 69 } 69 } 70 70 71 void G4SPSAngDistribution::SetAngDistType(cons 71 void G4SPSAngDistribution::SetAngDistType(const G4String& atype) 72 { 72 { 73 G4AutoLock l(&mutex); 73 G4AutoLock l(&mutex); 74 if(atype != "iso" && atype != "cos" && atype 74 if(atype != "iso" && atype != "cos" && atype != "user" && atype != "planar" 75 && atype != "beam1d" && atype != "beam2d" 75 && atype != "beam1d" && atype != "beam2d" && atype != "focused") 76 { 76 { 77 G4cout << "Error, distribution must be iso 77 G4cout << "Error, distribution must be iso, cos, planar, beam1d, beam2d, focused or user" 78 << G4endl; 78 << G4endl; 79 } 79 } 80 else 80 else 81 { 81 { 82 AngDistType = atype; 82 AngDistType = atype; 83 } 83 } 84 if (AngDistType == "cos") { MaxTheta = pi/2 84 if (AngDistType == "cos") { MaxTheta = pi/2.; } 85 if (AngDistType == "user") 85 if (AngDistType == "user") 86 { 86 { 87 UDefThetaH = IPDFThetaH = ZeroPhysVector; 87 UDefThetaH = IPDFThetaH = ZeroPhysVector; 88 IPDFThetaExist = false; 88 IPDFThetaExist = false; 89 UDefPhiH = IPDFPhiH = ZeroPhysVector; 89 UDefPhiH = IPDFPhiH = ZeroPhysVector; 90 IPDFPhiExist = false; 90 IPDFPhiExist = false; 91 } 91 } 92 } 92 } 93 93 94 void G4SPSAngDistribution::DefineAngRefAxes(co 94 void G4SPSAngDistribution::DefineAngRefAxes(const G4String& refname, 95 co 95 const G4ThreeVector& ref) 96 { 96 { 97 G4AutoLock l(&mutex); 97 G4AutoLock l(&mutex); 98 if (refname == "angref1") 98 if (refname == "angref1") 99 AngRef1 = ref.unit(); // x' 99 AngRef1 = ref.unit(); // x' 100 else if (refname == "angref2") 100 else if (refname == "angref2") 101 AngRef2 = ref.unit(); // vector in x'y' pl 101 AngRef2 = ref.unit(); // vector in x'y' plane 102 102 103 // User defines x' (AngRef1) and a vector in 103 // User defines x' (AngRef1) and a vector in the x'y' 104 // plane (AngRef2). Then, AngRef1 x AngRef2 104 // plane (AngRef2). Then, AngRef1 x AngRef2 = AngRef3 105 // the z' vector. Then, AngRef3 x AngRef1 = 105 // the z' vector. Then, AngRef3 x AngRef1 = AngRef2 106 // which will now be y'. 106 // which will now be y'. 107 107 108 AngRef3 = AngRef1.cross(AngRef2); // z' 108 AngRef3 = AngRef1.cross(AngRef2); // z' 109 AngRef2 = AngRef3.cross(AngRef1); // y' 109 AngRef2 = AngRef3.cross(AngRef1); // y' 110 UserAngRef = true ; 110 UserAngRef = true ; 111 if(verbosityLevel == 2) 111 if(verbosityLevel == 2) 112 { 112 { 113 G4cout << "Angular distribution rotation a 113 G4cout << "Angular distribution rotation axes " << AngRef1 114 << " " << AngRef2 << " " << AngRef3 114 << " " << AngRef2 << " " << AngRef3 << G4endl; 115 } 115 } 116 } 116 } 117 117 118 void G4SPSAngDistribution::SetMinTheta(G4doubl 118 void G4SPSAngDistribution::SetMinTheta(G4double mint) 119 { 119 { 120 G4AutoLock l(&mutex); 120 G4AutoLock l(&mutex); 121 MinTheta = mint; 121 MinTheta = mint; 122 } 122 } 123 123 124 void G4SPSAngDistribution::SetMinPhi(G4double 124 void G4SPSAngDistribution::SetMinPhi(G4double minp) 125 { 125 { 126 G4AutoLock l(&mutex); 126 G4AutoLock l(&mutex); 127 MinPhi = minp; 127 MinPhi = minp; 128 } 128 } 129 129 130 void G4SPSAngDistribution::SetMaxTheta(G4doubl 130 void G4SPSAngDistribution::SetMaxTheta(G4double maxt) 131 { 131 { 132 G4AutoLock l(&mutex); 132 G4AutoLock l(&mutex); 133 MaxTheta = maxt; 133 MaxTheta = maxt; 134 } 134 } 135 135 136 void G4SPSAngDistribution::SetMaxPhi(G4double 136 void G4SPSAngDistribution::SetMaxPhi(G4double maxp) 137 { 137 { 138 G4AutoLock l(&mutex); 138 G4AutoLock l(&mutex); 139 MaxPhi = maxp; 139 MaxPhi = maxp; 140 } 140 } 141 141 142 void G4SPSAngDistribution::SetBeamSigmaInAngR( 142 void G4SPSAngDistribution::SetBeamSigmaInAngR(G4double r) 143 { 143 { 144 G4AutoLock l(&mutex); 144 G4AutoLock l(&mutex); 145 DR = r; 145 DR = r; 146 } 146 } 147 147 148 void G4SPSAngDistribution::SetBeamSigmaInAngX( 148 void G4SPSAngDistribution::SetBeamSigmaInAngX(G4double r) 149 { 149 { 150 G4AutoLock l(&mutex); 150 G4AutoLock l(&mutex); 151 DX = r; 151 DX = r; 152 } 152 } 153 153 154 void G4SPSAngDistribution::SetBeamSigmaInAngY( 154 void G4SPSAngDistribution::SetBeamSigmaInAngY(G4double r) 155 { 155 { 156 G4AutoLock l(&mutex); 156 G4AutoLock l(&mutex); 157 DY = r; 157 DY = r; 158 } 158 } 159 159 160 void G4SPSAngDistribution:: 160 void G4SPSAngDistribution:: 161 SetParticleMomentumDirection(const G4ParticleM 161 SetParticleMomentumDirection(const G4ParticleMomentum& aMomentumDirection) 162 { 162 { 163 G4AutoLock l(&mutex); 163 G4AutoLock l(&mutex); 164 particle_momentum_direction = aMomentumDirec 164 particle_momentum_direction = aMomentumDirection.unit(); 165 } 165 } 166 166 167 void G4SPSAngDistribution::SetPosDistribution( 167 void G4SPSAngDistribution::SetPosDistribution(G4SPSPosDistribution* a) 168 { 168 { 169 G4AutoLock l(&mutex); 169 G4AutoLock l(&mutex); 170 posDist = a; 170 posDist = a; 171 } 171 } 172 172 173 void G4SPSAngDistribution::SetBiasRndm(G4SPSRa 173 void G4SPSAngDistribution::SetBiasRndm(G4SPSRandomGenerator* a) 174 { 174 { 175 G4AutoLock l(&mutex); 175 G4AutoLock l(&mutex); 176 angRndm = a; 176 angRndm = a; 177 } 177 } 178 178 179 void G4SPSAngDistribution::SetVerbosity(G4int 179 void G4SPSAngDistribution::SetVerbosity(G4int a) 180 { 180 { 181 G4AutoLock l(&mutex); 181 G4AutoLock l(&mutex); 182 verbosityLevel = a; 182 verbosityLevel = a; 183 } 183 } 184 184 185 void G4SPSAngDistribution::UserDefAngTheta(con 185 void G4SPSAngDistribution::UserDefAngTheta(const G4ThreeVector& input) 186 { 186 { 187 G4AutoLock l(&mutex); 187 G4AutoLock l(&mutex); 188 if(UserDistType == "NULL") UserDistType = "t 188 if(UserDistType == "NULL") UserDistType = "theta"; 189 if(UserDistType == "phi") UserDistType = "bo 189 if(UserDistType == "phi") UserDistType = "both"; 190 G4double thi, val; 190 G4double thi, val; 191 thi = input.x(); 191 thi = input.x(); 192 val = input.y(); 192 val = input.y(); 193 if(verbosityLevel >= 1) G4cout << "In UserDe 193 if(verbosityLevel >= 1) G4cout << "In UserDefAngTheta" << G4endl; 194 UDefThetaH.InsertValues(thi, val); 194 UDefThetaH.InsertValues(thi, val); 195 } 195 } 196 196 197 G4String G4SPSAngDistribution::GetDistType() 197 G4String G4SPSAngDistribution::GetDistType() 198 { 198 { 199 G4AutoLock l(&mutex); 199 G4AutoLock l(&mutex); 200 return AngDistType; 200 return AngDistType; 201 } 201 } 202 202 203 G4double G4SPSAngDistribution::GetMinTheta() 203 G4double G4SPSAngDistribution::GetMinTheta() 204 { 204 { 205 G4AutoLock l(&mutex); 205 G4AutoLock l(&mutex); 206 return MinTheta; 206 return MinTheta; 207 } 207 } 208 208 209 G4double G4SPSAngDistribution::GetMaxTheta() 209 G4double G4SPSAngDistribution::GetMaxTheta() 210 { 210 { 211 G4AutoLock l(&mutex); 211 G4AutoLock l(&mutex); 212 return MaxTheta; 212 return MaxTheta; 213 } 213 } 214 214 215 G4double G4SPSAngDistribution::GetMinPhi() 215 G4double G4SPSAngDistribution::GetMinPhi() 216 { 216 { 217 G4AutoLock l(&mutex); 217 G4AutoLock l(&mutex); 218 return MinPhi; 218 return MinPhi; 219 } 219 } 220 220 221 G4double G4SPSAngDistribution::GetMaxPhi() 221 G4double G4SPSAngDistribution::GetMaxPhi() 222 { 222 { 223 G4AutoLock l(&mutex); 223 G4AutoLock l(&mutex); 224 return MaxPhi; 224 return MaxPhi; 225 } 225 } 226 226 227 G4ThreeVector G4SPSAngDistribution::GetDirecti 227 G4ThreeVector G4SPSAngDistribution::GetDirection() 228 { 228 { 229 G4AutoLock l(&mutex); 229 G4AutoLock l(&mutex); 230 return particle_momentum_direction; 230 return particle_momentum_direction; 231 } 231 } 232 232 233 void G4SPSAngDistribution::UserDefAngPhi(const 233 void G4SPSAngDistribution::UserDefAngPhi(const G4ThreeVector& input) 234 { 234 { 235 G4AutoLock l(&mutex); 235 G4AutoLock l(&mutex); 236 if(UserDistType == "NULL") UserDistType = "p 236 if(UserDistType == "NULL") UserDistType = "phi"; 237 if(UserDistType == "theta") UserDistType = " 237 if(UserDistType == "theta") UserDistType = "both"; 238 G4double phhi, val; 238 G4double phhi, val; 239 phhi = input.x(); 239 phhi = input.x(); 240 val = input.y(); 240 val = input.y(); 241 if(verbosityLevel >= 1) G4cout << "In UserDe 241 if(verbosityLevel >= 1) G4cout << "In UserDefAngPhi" << G4endl; 242 UDefPhiH.InsertValues(phhi, val); 242 UDefPhiH.InsertValues(phhi, val); 243 } 243 } 244 244 245 void G4SPSAngDistribution::SetFocusPoint(const 245 void G4SPSAngDistribution::SetFocusPoint(const G4ThreeVector& input) 246 { 246 { 247 G4AutoLock l(&mutex); 247 G4AutoLock l(&mutex); 248 FocusPoint = input; 248 FocusPoint = input; 249 } 249 } 250 250 251 void G4SPSAngDistribution::SetUserWRTSurface(G 251 void G4SPSAngDistribution::SetUserWRTSurface(G4bool wrtSurf) 252 { 252 { 253 G4AutoLock l(&mutex); 253 G4AutoLock l(&mutex); 254 254 255 // if UserWRTSurface = true then the user wa 255 // if UserWRTSurface = true then the user wants momenta with respect 256 // to the surface normals. 256 // to the surface normals. 257 // When doing this theta has to be 0-90 only 257 // When doing this theta has to be 0-90 only otherwise there will be 258 // errors, which currently are flagged anywh 258 // errors, which currently are flagged anywhere. 259 // 259 // 260 UserWRTSurface = wrtSurf; 260 UserWRTSurface = wrtSurf; 261 } 261 } 262 262 263 void G4SPSAngDistribution::SetUseUserAngAxis(G 263 void G4SPSAngDistribution::SetUseUserAngAxis(G4bool userang) 264 { 264 { 265 G4AutoLock l(&mutex); 265 G4AutoLock l(&mutex); 266 266 267 // if UserAngRef = true the angular distrib 267 // if UserAngRef = true the angular distribution is defined wrt 268 // the user defined coordinates 268 // the user defined coordinates 269 // 269 // 270 UserAngRef = userang; 270 UserAngRef = userang; 271 } 271 } 272 272 273 void G4SPSAngDistribution::GenerateBeamFlux(G4 273 void G4SPSAngDistribution::GenerateBeamFlux(G4ParticleMomentum& mom) 274 { 274 { 275 G4double theta, phi; 275 G4double theta, phi; 276 G4double px, py, pz; 276 G4double px, py, pz; 277 if (AngDistType == "beam1d") 277 if (AngDistType == "beam1d") 278 { 278 { 279 theta = G4RandGauss::shoot(0.0,DR); 279 theta = G4RandGauss::shoot(0.0,DR); 280 phi = twopi * G4UniformRand(); 280 phi = twopi * G4UniformRand(); 281 } 281 } 282 else 282 else 283 { 283 { 284 px = G4RandGauss::shoot(0.0,DX); 284 px = G4RandGauss::shoot(0.0,DX); 285 py = G4RandGauss::shoot(0.0,DY); 285 py = G4RandGauss::shoot(0.0,DY); 286 theta = std::sqrt (px*px + py*py); 286 theta = std::sqrt (px*px + py*py); 287 if (theta != 0.) 287 if (theta != 0.) 288 { 288 { 289 phi = std::acos(px/theta); 289 phi = std::acos(px/theta); 290 if ( py < 0.) phi = -phi; 290 if ( py < 0.) phi = -phi; 291 } 291 } 292 else 292 else 293 { 293 { 294 phi = 0.0; 294 phi = 0.0; 295 } 295 } 296 } 296 } 297 px = -std::sin(theta) * std::cos(phi); 297 px = -std::sin(theta) * std::cos(phi); 298 py = -std::sin(theta) * std::sin(phi); 298 py = -std::sin(theta) * std::sin(phi); 299 pz = -std::cos(theta); 299 pz = -std::cos(theta); 300 G4double finx, finy, finz; 300 G4double finx, finy, finz; 301 finx=px, finy=py, finz=pz; 301 finx=px, finy=py, finz=pz; 302 if (UserAngRef) 302 if (UserAngRef) 303 { 303 { 304 // Apply Angular Rotation Matrix 304 // Apply Angular Rotation Matrix 305 // x * AngRef1, y * AngRef2 and z * AngRef 305 // x * AngRef1, y * AngRef2 and z * AngRef3 306 finx = (px * AngRef1.x()) + (py * AngRef2. 306 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); 307 finy = (px * AngRef1.y()) + (py * AngRef2. 307 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); 308 finz = (px * AngRef1.z()) + (py * AngRef2. 308 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); 309 G4double ResMag = std::sqrt((finx*finx) + 309 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz)); 310 finx = finx/ResMag; 310 finx = finx/ResMag; 311 finy = finy/ResMag; 311 finy = finy/ResMag; 312 finz = finz/ResMag; 312 finz = finz/ResMag; 313 } 313 } 314 mom.setX(finx); 314 mom.setX(finx); 315 mom.setY(finy); 315 mom.setY(finy); 316 mom.setZ(finz); 316 mom.setZ(finz); 317 317 318 // particle_momentum_direction now holds uni 318 // particle_momentum_direction now holds unit momentum vector 319 319 320 if(verbosityLevel >= 1) 320 if(verbosityLevel >= 1) 321 { 321 { 322 G4cout << "Generating beam vector: " << mo 322 G4cout << "Generating beam vector: " << mom << G4endl; 323 } 323 } 324 } 324 } 325 325 326 void G4SPSAngDistribution::GenerateFocusedFlux 326 void G4SPSAngDistribution::GenerateFocusedFlux(G4ParticleMomentum& mom) 327 { 327 { 328 mom = (FocusPoint - posDist->GetParticlePos( 328 mom = (FocusPoint - posDist->GetParticlePos()).unit(); 329 329 330 // particle_momentum_direction now holds uni 330 // particle_momentum_direction now holds unit momentum vector. 331 331 332 if(verbosityLevel >= 1) 332 if(verbosityLevel >= 1) 333 { 333 { 334 G4cout << "Generating focused vector: " << 334 G4cout << "Generating focused vector: " << mom << G4endl; 335 } 335 } 336 } 336 } 337 337 338 void G4SPSAngDistribution::GenerateIsotropicFl 338 void G4SPSAngDistribution::GenerateIsotropicFlux(G4ParticleMomentum& mom) 339 { 339 { 340 // generates isotropic flux. 340 // generates isotropic flux. 341 // No vectors are needed. 341 // No vectors are needed. 342 342 343 G4double rndm, rndm2; 343 G4double rndm, rndm2; 344 G4double px, py, pz; 344 G4double px, py, pz; 345 345 346 G4double sintheta, sinphi,costheta,cosphi; 346 G4double sintheta, sinphi,costheta,cosphi; 347 rndm = angRndm->GenRandTheta(); 347 rndm = angRndm->GenRandTheta(); 348 costheta = std::cos(MinTheta) - rndm * (std: 348 costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta) 349 - std: 349 - std::cos(MaxTheta)); 350 sintheta = std::sqrt(1. - costheta*costheta) 350 sintheta = std::sqrt(1. - costheta*costheta); 351 351 352 rndm2 = angRndm->GenRandPhi(); 352 rndm2 = angRndm->GenRandPhi(); 353 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; 353 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; 354 sinphi = std::sin(Phi); 354 sinphi = std::sin(Phi); 355 cosphi = std::cos(Phi); 355 cosphi = std::cos(Phi); 356 356 357 px = -sintheta * cosphi; 357 px = -sintheta * cosphi; 358 py = -sintheta * sinphi; 358 py = -sintheta * sinphi; 359 pz = -costheta; 359 pz = -costheta; 360 360 361 // For volume and point source use mother or 361 // For volume and point source use mother or user defined coordinates 362 // for plane and surface source user surface 362 // for plane and surface source user surface-normal or user-defined 363 // coordinates 363 // coordinates 364 // 364 // 365 G4double finx, finy, finz; 365 G4double finx, finy, finz; 366 if (posDist->GetSourcePosType() == "Point" 366 if (posDist->GetSourcePosType() == "Point" 367 || posDist->GetSourcePosType() == "Volume") 367 || posDist->GetSourcePosType() == "Volume") 368 { 368 { 369 if (UserAngRef) 369 if (UserAngRef) 370 { 370 { 371 // Apply Rotation Matrix 371 // Apply Rotation Matrix 372 // x * AngRef1, y * AngRef2 and z * AngR 372 // x * AngRef1, y * AngRef2 and z * AngRef3 373 finx = (px * AngRef1.x()) + (py * AngRef 373 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); 374 finy = (px * AngRef1.y()) + (py * AngRef 374 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); 375 finz = (px * AngRef1.z()) + (py * AngRef 375 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); 376 } 376 } 377 else 377 else 378 { 378 { 379 finx = px; 379 finx = px; 380 finy = py; 380 finy = py; 381 finz = pz; 381 finz = pz; 382 } 382 } 383 } 383 } 384 else 384 else 385 { // for plane and surface source 385 { // for plane and surface source 386 if (UserAngRef) 386 if (UserAngRef) 387 { 387 { 388 // Apply Rotation Matrix 388 // Apply Rotation Matrix 389 // x * AngRef1, y * AngRef2 and z * AngR 389 // x * AngRef1, y * AngRef2 and z * AngRef3 390 finx = (px * AngRef1.x()) + (py * AngRef 390 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); 391 finy = (px * AngRef1.y()) + (py * AngRef 391 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); 392 finz = (px * AngRef1.z()) + (py * AngRef 392 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); 393 } 393 } 394 else 394 else 395 { 395 { 396 finx = (px*posDist->GetSideRefVec1().x() 396 finx = (px*posDist->GetSideRefVec1().x()) 397 + (py*posDist->GetSideRefVec2().x() 397 + (py*posDist->GetSideRefVec2().x()) 398 + (pz*posDist->GetSideRefVec3().x() 398 + (pz*posDist->GetSideRefVec3().x()); 399 finy = (px*posDist->GetSideRefVec1().y() 399 finy = (px*posDist->GetSideRefVec1().y()) 400 + (py*posDist->GetSideRefVec2().y() 400 + (py*posDist->GetSideRefVec2().y()) 401 + (pz*posDist->GetSideRefVec3().y() 401 + (pz*posDist->GetSideRefVec3().y()); 402 finz = (px*posDist->GetSideRefVec1().z() 402 finz = (px*posDist->GetSideRefVec1().z()) 403 + (py*posDist->GetSideRefVec2().z() 403 + (py*posDist->GetSideRefVec2().z()) 404 + (pz*posDist->GetSideRefVec3().z() 404 + (pz*posDist->GetSideRefVec3().z()); 405 } 405 } 406 } 406 } 407 G4double ResMag = std::sqrt((finx*finx) + (f 407 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz)); 408 finx = finx/ResMag; 408 finx = finx/ResMag; 409 finy = finy/ResMag; 409 finy = finy/ResMag; 410 finz = finz/ResMag; 410 finz = finz/ResMag; 411 411 412 mom.setX(finx); 412 mom.setX(finx); 413 mom.setY(finy); 413 mom.setY(finy); 414 mom.setZ(finz); 414 mom.setZ(finz); 415 415 416 // particle_momentum_direction now holds uni 416 // particle_momentum_direction now holds unit momentum vector. 417 417 418 if(verbosityLevel >= 1) 418 if(verbosityLevel >= 1) 419 { 419 { 420 G4cout << "Generating isotropic vector: " 420 G4cout << "Generating isotropic vector: " << mom << G4endl; 421 } 421 } 422 } 422 } 423 423 424 void G4SPSAngDistribution::GenerateCosineLawFl 424 void G4SPSAngDistribution::GenerateCosineLawFlux(G4ParticleMomentum& mom) 425 { 425 { 426 // Method to generate flux distributed with 426 // Method to generate flux distributed with a cosine law 427 427 428 G4double px, py, pz; 428 G4double px, py, pz; 429 G4double rndm, rndm2; 429 G4double rndm, rndm2; 430 430 431 G4double sintheta, sinphi,costheta,cosphi; 431 G4double sintheta, sinphi,costheta,cosphi; 432 rndm = angRndm->GenRandTheta(); 432 rndm = angRndm->GenRandTheta(); 433 sintheta = std::sqrt( rndm * (std::sin(MaxTh 433 sintheta = std::sqrt( rndm * (std::sin(MaxTheta)*std::sin(MaxTheta) 434 - std::sin(MinTh 434 - std::sin(MinTheta)*std::sin(MinTheta) ) 435 + std::sin(MinTheta)*std 435 + std::sin(MinTheta)*std::sin(MinTheta) ); 436 costheta = std::sqrt(1. -sintheta*sintheta); 436 costheta = std::sqrt(1. -sintheta*sintheta); 437 437 438 rndm2 = angRndm->GenRandPhi(); 438 rndm2 = angRndm->GenRandPhi(); 439 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; 439 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; 440 sinphi = std::sin(Phi); 440 sinphi = std::sin(Phi); 441 cosphi = std::cos(Phi); 441 cosphi = std::cos(Phi); 442 442 443 px = -sintheta * cosphi; 443 px = -sintheta * cosphi; 444 py = -sintheta * sinphi; 444 py = -sintheta * sinphi; 445 pz = -costheta; 445 pz = -costheta; 446 446 447 // for volume and point source use mother or 447 // for volume and point source use mother or user defined coordinates 448 // for plane and surface source user surface 448 // for plane and surface source user surface-normal or userdefined 449 // coordinates 449 // coordinates 450 // 450 // 451 G4double finx, finy, finz; 451 G4double finx, finy, finz; 452 if (posDist->GetSourcePosType() == "Point" 452 if (posDist->GetSourcePosType() == "Point" 453 || posDist->GetSourcePosType() == "Volume") 453 || posDist->GetSourcePosType() == "Volume") 454 { 454 { 455 if (UserAngRef) 455 if (UserAngRef) 456 { 456 { 457 // Apply Rotation Matrix 457 // Apply Rotation Matrix 458 finx = (px * AngRef1.x()) + (py * AngRef 458 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); 459 finy = (px * AngRef1.y()) + (py * AngRef 459 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); 460 finz = (px * AngRef1.z()) + (py * AngRef 460 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); 461 } 461 } 462 else 462 else 463 { 463 { 464 finx = px; 464 finx = px; 465 finy = py; 465 finy = py; 466 finz = pz; 466 finz = pz; 467 } 467 } 468 } 468 } 469 else 469 else 470 { // for plane and surface source 470 { // for plane and surface source 471 if (UserAngRef) 471 if (UserAngRef) 472 { 472 { 473 // Apply Rotation Matrix 473 // Apply Rotation Matrix 474 finx = (px * AngRef1.x()) + (py * AngRef 474 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); 475 finy = (px * AngRef1.y()) + (py * AngRef 475 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); 476 finz = (px * AngRef1.z()) + (py * AngRef 476 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); 477 } 477 } 478 else 478 else 479 { 479 { 480 finx = (px*posDist->GetSideRefVec1().x() 480 finx = (px*posDist->GetSideRefVec1().x()) 481 + (py*posDist->GetSideRefVec2().x() 481 + (py*posDist->GetSideRefVec2().x()) 482 + (pz*posDist->GetSideRefVec3().x() 482 + (pz*posDist->GetSideRefVec3().x()); 483 finy = (px*posDist->GetSideRefVec1().y() 483 finy = (px*posDist->GetSideRefVec1().y()) 484 + (py*posDist->GetSideRefVec2().y() 484 + (py*posDist->GetSideRefVec2().y()) 485 + (pz*posDist->GetSideRefVec3().y() 485 + (pz*posDist->GetSideRefVec3().y()); 486 finz = (px*posDist->GetSideRefVec1().z() 486 finz = (px*posDist->GetSideRefVec1().z()) 487 + (py*posDist->GetSideRefVec2().z() 487 + (py*posDist->GetSideRefVec2().z()) 488 + (pz*posDist->GetSideRefVec3().z() 488 + (pz*posDist->GetSideRefVec3().z()); 489 } 489 } 490 } 490 } 491 G4double ResMag = std::sqrt((finx*finx) + (f 491 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz)); 492 finx = finx/ResMag; 492 finx = finx/ResMag; 493 finy = finy/ResMag; 493 finy = finy/ResMag; 494 finz = finz/ResMag; 494 finz = finz/ResMag; 495 495 496 mom.setX(finx); 496 mom.setX(finx); 497 mom.setY(finy); 497 mom.setY(finy); 498 mom.setZ(finz); 498 mom.setZ(finz); 499 499 500 // particle_momentum_direction now contains 500 // particle_momentum_direction now contains unit momentum vector. 501 501 502 if(verbosityLevel >= 1) 502 if(verbosityLevel >= 1) 503 { 503 { 504 G4cout << "Resultant cosine-law unit momen 504 G4cout << "Resultant cosine-law unit momentum vector " << mom << G4endl; 505 } 505 } 506 } 506 } 507 507 508 void G4SPSAngDistribution::GeneratePlanarFlux( 508 void G4SPSAngDistribution::GeneratePlanarFlux(G4ParticleMomentum& mom) 509 { 509 { 510 // particle_momentum_direction now contains 510 // particle_momentum_direction now contains unit momentum vector. 511 // nothing need be done here as the m-direct 511 // nothing need be done here as the m-directions have been set directly 512 // under this option 512 // under this option 513 513 514 if(verbosityLevel >= 1) 514 if(verbosityLevel >= 1) 515 { 515 { 516 G4cout << "Resultant Planar wave momentum 516 G4cout << "Resultant Planar wave momentum vector " << mom << G4endl; 517 } 517 } 518 } 518 } 519 519 520 void G4SPSAngDistribution::GenerateUserDefFlux 520 void G4SPSAngDistribution::GenerateUserDefFlux(G4ParticleMomentum& mom) 521 { 521 { 522 G4double rndm, px, py, pz, pmag; 522 G4double rndm, px, py, pz, pmag; 523 523 524 if(UserDistType == "NULL") 524 if(UserDistType == "NULL") 525 { 525 { 526 G4cout << "Error: UserDistType undefined" 526 G4cout << "Error: UserDistType undefined" << G4endl; 527 } 527 } 528 else if(UserDistType == "theta") 528 else if(UserDistType == "theta") 529 { 529 { 530 Theta = 10.; 530 Theta = 10.; 531 while(Theta > MaxTheta || Theta < MinTheta 531 while(Theta > MaxTheta || Theta < MinTheta) 532 { 532 { 533 Theta = GenerateUserDefTheta(); 533 Theta = GenerateUserDefTheta(); 534 } 534 } 535 Phi = 10.; 535 Phi = 10.; 536 while(Phi > MaxPhi || Phi < MinPhi) 536 while(Phi > MaxPhi || Phi < MinPhi) 537 { 537 { 538 rndm = angRndm->GenRandPhi(); 538 rndm = angRndm->GenRandPhi(); 539 Phi = twopi * rndm; 539 Phi = twopi * rndm; 540 } 540 } 541 } 541 } 542 else if(UserDistType == "phi") 542 else if(UserDistType == "phi") 543 { 543 { 544 Theta = 10.; 544 Theta = 10.; 545 while(Theta > MaxTheta || Theta < MinTheta 545 while(Theta > MaxTheta || Theta < MinTheta) 546 { 546 { 547 rndm = angRndm->GenRandTheta(); 547 rndm = angRndm->GenRandTheta(); 548 Theta = std::acos(1. - (2. * rndm)); 548 Theta = std::acos(1. - (2. * rndm)); 549 } 549 } 550 Phi = 10.; 550 Phi = 10.; 551 while(Phi > MaxPhi || Phi < MinPhi) 551 while(Phi > MaxPhi || Phi < MinPhi) 552 { 552 { 553 Phi = GenerateUserDefPhi(); 553 Phi = GenerateUserDefPhi(); 554 } 554 } 555 } 555 } 556 else if(UserDistType == "both") 556 else if(UserDistType == "both") 557 { 557 { 558 Theta = 10.; 558 Theta = 10.; 559 while(Theta > MaxTheta || Theta < MinTheta 559 while(Theta > MaxTheta || Theta < MinTheta) 560 { 560 { 561 Theta = GenerateUserDefTheta(); 561 Theta = GenerateUserDefTheta(); 562 } 562 } 563 Phi = 10.; 563 Phi = 10.; 564 while(Phi > MaxPhi || Phi < MinPhi) 564 while(Phi > MaxPhi || Phi < MinPhi) 565 { 565 { 566 Phi = GenerateUserDefPhi(); 566 Phi = GenerateUserDefPhi(); 567 } 567 } 568 } 568 } 569 px = -std::sin(Theta) * std::cos(Phi); 569 px = -std::sin(Theta) * std::cos(Phi); 570 py = -std::sin(Theta) * std::sin(Phi); 570 py = -std::sin(Theta) * std::sin(Phi); 571 pz = -std::cos(Theta); 571 pz = -std::cos(Theta); 572 572 573 pmag = std::sqrt((px*px) + (py*py) + (pz*pz) 573 pmag = std::sqrt((px*px) + (py*py) + (pz*pz)); 574 574 575 if(!UserWRTSurface) 575 if(!UserWRTSurface) 576 { 576 { 577 G4double finx, finy, finz; 577 G4double finx, finy, finz; 578 if (UserAngRef) 578 if (UserAngRef) 579 { 579 { 580 // Apply Rotation Matrix 580 // Apply Rotation Matrix 581 // x * AngRef1, y * AngRef2 and z * AngR 581 // x * AngRef1, y * AngRef2 and z * AngRef3 582 finx = (px * AngRef1.x()) + (py * AngRef 582 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); 583 finy = (px * AngRef1.y()) + (py * AngRef 583 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); 584 finz = (px * AngRef1.z()) + (py * AngRef 584 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); 585 } 585 } 586 else // use mother coordinates 586 else // use mother coordinates 587 { 587 { 588 finx = px; 588 finx = px; 589 finy = py; 589 finy = py; 590 finz = pz; 590 finz = pz; 591 } 591 } 592 G4double ResMag = std::sqrt((finx*finx) + 592 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz)); 593 finx = finx/ResMag; 593 finx = finx/ResMag; 594 finy = finy/ResMag; 594 finy = finy/ResMag; 595 finz = finz/ResMag; 595 finz = finz/ResMag; 596 596 597 mom.setX(finx); 597 mom.setX(finx); 598 mom.setY(finy); 598 mom.setY(finy); 599 mom.setZ(finz); 599 mom.setZ(finz); 600 } 600 } 601 else // UserWRTSurface = true 601 else // UserWRTSurface = true 602 { 602 { 603 G4double pxh = px/pmag; 603 G4double pxh = px/pmag; 604 G4double pyh = py/pmag; 604 G4double pyh = py/pmag; 605 G4double pzh = pz/pmag; 605 G4double pzh = pz/pmag; 606 if(verbosityLevel > 1) 606 if(verbosityLevel > 1) 607 { 607 { 608 G4cout << "SideRefVecs " << posDist->Get 608 G4cout << "SideRefVecs " << posDist->GetSideRefVec1() 609 << posDist->GetSideRefVec2() << p 609 << posDist->GetSideRefVec2() << posDist->GetSideRefVec3() 610 << G4endl; 610 << G4endl; 611 G4cout << "Raw Unit vector " << pxh 611 G4cout << "Raw Unit vector " << pxh 612 << "," << pyh << "," << pzh << G4 612 << "," << pyh << "," << pzh << G4endl; 613 } 613 } 614 G4double resultx = (pxh*posDist->GetSideRe 614 G4double resultx = (pxh*posDist->GetSideRefVec1().x()) 615 + (pyh*posDist->GetSideRe 615 + (pyh*posDist->GetSideRefVec2().x()) 616 + (pzh*posDist->GetSideRe 616 + (pzh*posDist->GetSideRefVec3().x()); 617 617 618 G4double resulty = (pxh*posDist->GetSideRe 618 G4double resulty = (pxh*posDist->GetSideRefVec1().y()) 619 + (pyh*posDist->GetSideRe 619 + (pyh*posDist->GetSideRefVec2().y()) 620 + (pzh*posDist->GetSideRe 620 + (pzh*posDist->GetSideRefVec3().y()); 621 621 622 G4double resultz = (pxh*posDist->GetSideRe 622 G4double resultz = (pxh*posDist->GetSideRefVec1().z()) 623 + (pyh*posDist->GetSideRe 623 + (pyh*posDist->GetSideRefVec2().z()) 624 + (pzh*posDist->GetSideRe 624 + (pzh*posDist->GetSideRefVec3().z()); 625 625 626 G4double ResMag = std::sqrt((resultx*resul 626 G4double ResMag = std::sqrt((resultx*resultx) 627 + (resulty*resul 627 + (resulty*resulty) 628 + (resultz*resul 628 + (resultz*resultz)); 629 resultx = resultx/ResMag; 629 resultx = resultx/ResMag; 630 resulty = resulty/ResMag; 630 resulty = resulty/ResMag; 631 resultz = resultz/ResMag; 631 resultz = resultz/ResMag; 632 632 633 mom.setX(resultx); 633 mom.setX(resultx); 634 mom.setY(resulty); 634 mom.setY(resulty); 635 mom.setZ(resultz); 635 mom.setZ(resultz); 636 } 636 } 637 637 638 // particle_momentum_direction now contains 638 // particle_momentum_direction now contains unit momentum vector. 639 639 640 if(verbosityLevel > 0 ) 640 if(verbosityLevel > 0 ) 641 { 641 { 642 G4cout << "Final User Defined momentum vec 642 G4cout << "Final User Defined momentum vector " 643 << particle_momentum_direction << G 643 << particle_momentum_direction << G4endl; 644 } 644 } 645 } 645 } 646 646 647 G4double G4SPSAngDistribution::GenerateUserDef 647 G4double G4SPSAngDistribution::GenerateUserDefTheta() 648 { 648 { 649 // Create cumulative histogram if not alread 649 // Create cumulative histogram if not already done so. 650 // Then use RandFlat::shoot to generate the 650 // Then use RandFlat::shoot to generate the output Theta value. 651 651 652 if(UserDistType == "NULL" || UserDistType == 652 if(UserDistType == "NULL" || UserDistType == "phi") 653 { 653 { 654 // No user defined theta distribution 654 // No user defined theta distribution 655 G4cout << "Error ***********************" 655 G4cout << "Error ***********************" << G4endl; 656 G4cout << "UserDistType = " << UserDistTyp 656 G4cout << "UserDistType = " << UserDistType << G4endl; 657 return (0.); 657 return (0.); 658 } 658 } 659 659 660 // UserDistType = theta or both and so a the 660 // UserDistType = theta or both and so a theta distribution 661 // is defined. This should be integrated if 661 // is defined. This should be integrated if not already done. 662 G4AutoLock l(&mutex); 662 G4AutoLock l(&mutex); 663 if(!IPDFThetaExist) 663 if(!IPDFThetaExist) 664 { 664 { 665 // IPDF has not been created, so create it 665 // IPDF has not been created, so create it 666 // 666 // 667 G4double bins[1024],vals[1024], sum; 667 G4double bins[1024],vals[1024], sum; 668 G4int ii; 668 G4int ii; 669 auto maxbin = G4int(UDefThetaH.GetVectorL << 669 G4int maxbin = G4int(UDefThetaH.GetVectorLength()); 670 bins[0] = UDefThetaH.GetLowEdgeEnergy(std: 670 bins[0] = UDefThetaH.GetLowEdgeEnergy(std::size_t(0)); 671 vals[0] = UDefThetaH(std::size_t(0)); 671 vals[0] = UDefThetaH(std::size_t(0)); 672 sum = vals[0]; 672 sum = vals[0]; 673 for(ii=1; ii<maxbin; ++ii) 673 for(ii=1; ii<maxbin; ++ii) 674 { 674 { 675 bins[ii] = UDefThetaH.GetLowEdgeEnergy(s 675 bins[ii] = UDefThetaH.GetLowEdgeEnergy(std::size_t(ii)); 676 vals[ii] = UDefThetaH(std::size_t(ii)) + 676 vals[ii] = UDefThetaH(std::size_t(ii)) + vals[ii-1]; 677 sum = sum + UDefThetaH(std::size_t(ii)); 677 sum = sum + UDefThetaH(std::size_t(ii)); 678 } 678 } 679 for(ii=0; ii<maxbin; ++ii) 679 for(ii=0; ii<maxbin; ++ii) 680 { 680 { 681 vals[ii] = vals[ii]/sum; 681 vals[ii] = vals[ii]/sum; 682 IPDFThetaH.InsertValues(bins[ii], vals[i 682 IPDFThetaH.InsertValues(bins[ii], vals[ii]); 683 } 683 } 684 IPDFThetaExist = true; 684 IPDFThetaExist = true; 685 } 685 } 686 l.unlock(); 686 l.unlock(); 687 687 688 // IPDF has been created so carry on 688 // IPDF has been created so carry on 689 // 689 // 690 G4double rndm = G4UniformRand(); 690 G4double rndm = G4UniformRand(); 691 return IPDFThetaH.GetEnergy(rndm); 691 return IPDFThetaH.GetEnergy(rndm); 692 } 692 } 693 693 694 G4double G4SPSAngDistribution::GenerateUserDef 694 G4double G4SPSAngDistribution::GenerateUserDefPhi() 695 { 695 { 696 // Create cumulative histogram if not alread 696 // Create cumulative histogram if not already done so. 697 // Then use RandFlat::shoot to generate the 697 // Then use RandFlat::shoot to generate the output Theta value. 698 698 699 if(UserDistType == "NULL" || UserDistType == 699 if(UserDistType == "NULL" || UserDistType == "theta") 700 { 700 { 701 // No user defined phi distribution 701 // No user defined phi distribution 702 G4cout << "Error ***********************" 702 G4cout << "Error ***********************" << G4endl; 703 G4cout << "UserDistType = " << UserDistTyp 703 G4cout << "UserDistType = " << UserDistType << G4endl; 704 return(0.); 704 return(0.); 705 } 705 } 706 706 707 // UserDistType = phi or both and so a phi d 707 // UserDistType = phi or both and so a phi distribution 708 // is defined. This should be integrated if 708 // is defined. This should be integrated if not already done. 709 G4AutoLock l(&mutex); 709 G4AutoLock l(&mutex); 710 if(!IPDFPhiExist) 710 if(!IPDFPhiExist) 711 { 711 { 712 // IPDF has not been created, so create it 712 // IPDF has not been created, so create it 713 // 713 // 714 G4double bins[1024],vals[1024], sum; 714 G4double bins[1024],vals[1024], sum; 715 G4int ii; 715 G4int ii; 716 auto maxbin = G4int(UDefPhiH.GetVectorLen << 716 G4int maxbin = G4int(UDefPhiH.GetVectorLength()); 717 bins[0] = UDefPhiH.GetLowEdgeEnergy(std::s 717 bins[0] = UDefPhiH.GetLowEdgeEnergy(std::size_t(0)); 718 vals[0] = UDefPhiH(std::size_t(0)); 718 vals[0] = UDefPhiH(std::size_t(0)); 719 sum = vals[0]; 719 sum = vals[0]; 720 for(ii=1; ii<maxbin; ++ii) 720 for(ii=1; ii<maxbin; ++ii) 721 { 721 { 722 bins[ii] = UDefPhiH.GetLowEdgeEnergy(std 722 bins[ii] = UDefPhiH.GetLowEdgeEnergy(std::size_t(ii)); 723 vals[ii] = UDefPhiH(std::size_t(ii)) + v 723 vals[ii] = UDefPhiH(std::size_t(ii)) + vals[ii-1]; 724 sum = sum + UDefPhiH(std::size_t(ii)); 724 sum = sum + UDefPhiH(std::size_t(ii)); 725 } 725 } 726 for(ii=0; ii<maxbin; ++ii) 726 for(ii=0; ii<maxbin; ++ii) 727 { 727 { 728 vals[ii] = vals[ii]/sum; 728 vals[ii] = vals[ii]/sum; 729 IPDFPhiH.InsertValues(bins[ii], vals[ii] 729 IPDFPhiH.InsertValues(bins[ii], vals[ii]); 730 } 730 } 731 IPDFPhiExist = true; 731 IPDFPhiExist = true; 732 } 732 } 733 l.unlock(); 733 l.unlock(); 734 734 735 // IPDF has been create so carry on 735 // IPDF has been create so carry on 736 // 736 // 737 G4double rndm = G4UniformRand(); 737 G4double rndm = G4UniformRand(); 738 return IPDFPhiH.GetEnergy(rndm); 738 return IPDFPhiH.GetEnergy(rndm); 739 } 739 } 740 740 741 void G4SPSAngDistribution::ReSetHist(const G4S 741 void G4SPSAngDistribution::ReSetHist(const G4String& atype) 742 { 742 { 743 G4AutoLock l(&mutex); 743 G4AutoLock l(&mutex); 744 if (atype == "theta") 744 if (atype == "theta") 745 { 745 { 746 UDefThetaH = IPDFThetaH = ZeroPhysVector ; 746 UDefThetaH = IPDFThetaH = ZeroPhysVector ; 747 IPDFThetaExist = false ; 747 IPDFThetaExist = false ; 748 } 748 } 749 else if (atype == "phi") 749 else if (atype == "phi") 750 { 750 { 751 UDefPhiH = IPDFPhiH = ZeroPhysVector ; 751 UDefPhiH = IPDFPhiH = ZeroPhysVector ; 752 IPDFPhiExist = false ; 752 IPDFPhiExist = false ; 753 } 753 } 754 else 754 else 755 { 755 { 756 G4cout << "Error, histtype not accepted " 756 G4cout << "Error, histtype not accepted " << G4endl; 757 } 757 } 758 } 758 } 759 759 760 G4ParticleMomentum G4SPSAngDistribution::Gener 760 G4ParticleMomentum G4SPSAngDistribution::GenerateOne() 761 { 761 { 762 // Local copy for thread safety 762 // Local copy for thread safety 763 // 763 // 764 G4ParticleMomentum localM = particle_momentu 764 G4ParticleMomentum localM = particle_momentum_direction; 765 765 766 // Angular stuff 766 // Angular stuff 767 // 767 // 768 if(AngDistType == "iso") 768 if(AngDistType == "iso") 769 GenerateIsotropicFlux(localM); 769 GenerateIsotropicFlux(localM); 770 else if(AngDistType == "cos") 770 else if(AngDistType == "cos") 771 GenerateCosineLawFlux(localM); 771 GenerateCosineLawFlux(localM); 772 else if(AngDistType == "planar") 772 else if(AngDistType == "planar") 773 GeneratePlanarFlux(localM); 773 GeneratePlanarFlux(localM); 774 else if(AngDistType == "beam1d" || AngDistTy 774 else if(AngDistType == "beam1d" || AngDistType == "beam2d" ) 775 GenerateBeamFlux(localM); 775 GenerateBeamFlux(localM); 776 else if(AngDistType == "user") 776 else if(AngDistType == "user") 777 GenerateUserDefFlux(localM); 777 GenerateUserDefFlux(localM); 778 else if(AngDistType == "focused") 778 else if(AngDistType == "focused") 779 GenerateFocusedFlux(localM); 779 GenerateFocusedFlux(localM); 780 else 780 else 781 G4cout << "Error: AngDistType has unusual 781 G4cout << "Error: AngDistType has unusual value" << G4endl; 782 return localM; 782 return localM; 783 } 783 } 784 784