44 AngDistType =
"planar";
56 UserDistType =
"NULL";
57 UserWRTSurface =
true;
59 IPDFThetaExist =
false;
74 if(atype !=
"iso" && atype !=
"cos" && atype !=
"user" && atype !=
"planar"
75 && atype !=
"beam1d" && atype !=
"beam2d" && atype !=
"focused")
77 G4cout <<
"Error, distribution must be iso, cos, planar, beam1d, beam2d, focused or user"
84 if (AngDistType ==
"cos") { MaxTheta = pi/2.; }
85 if (AngDistType ==
"user")
87 UDefThetaH = IPDFThetaH = ZeroPhysVector;
88 IPDFThetaExist =
false;
89 UDefPhiH = IPDFPhiH = ZeroPhysVector;
98 if (refname ==
"angref1")
100 else if (refname ==
"angref2")
101 AngRef2 = ref.
unit();
108 AngRef3 = AngRef1.
cross(AngRef2);
109 AngRef2 = AngRef3.
cross(AngRef1);
111 if(verbosityLevel == 2)
113 G4cout <<
"Angular distribution rotation axes " << AngRef1
114 <<
" " << AngRef2 <<
" " << AngRef3 <<
G4endl;
164 particle_momentum_direction = aMomentumDirection.
unit();
188 if(UserDistType ==
"NULL") UserDistType =
"theta";
189 if(UserDistType ==
"phi") UserDistType =
"both";
193 if(verbosityLevel >= 1)
G4cout <<
"In UserDefAngTheta" <<
G4endl;
194 UDefThetaH.InsertValues(thi, val);
230 return particle_momentum_direction;
236 if(UserDistType ==
"NULL") UserDistType =
"phi";
237 if(UserDistType ==
"theta") UserDistType =
"both";
241 if(verbosityLevel >= 1)
G4cout <<
"In UserDefAngPhi" <<
G4endl;
242 UDefPhiH.InsertValues(phhi, val);
260 UserWRTSurface = wrtSurf;
270 UserAngRef = userang;
277 if (AngDistType ==
"beam1d")
279 theta = G4RandGauss::shoot(0.0,DR);
284 px = G4RandGauss::shoot(0.0,DX);
285 py = G4RandGauss::shoot(0.0,DY);
286 theta = std::sqrt (px*px + py*py);
289 phi = std::acos(px/theta);
290 if ( py < 0.) phi = -phi;
297 px = -std::sin(theta) * std::cos(phi);
298 py = -std::sin(theta) * std::sin(phi);
299 pz = -std::cos(theta);
301 finx=px, finy=py, finz=pz;
306 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
307 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
308 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
309 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
320 if(verbosityLevel >= 1)
328 mom = (FocusPoint - posDist->GetParticlePos()).
unit();
332 if(verbosityLevel >= 1)
334 G4cout <<
"Generating focused vector: " << mom <<
G4endl;
346 G4double sintheta, sinphi,costheta,cosphi;
347 rndm = angRndm->GenRandTheta();
348 costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta)
349 - std::cos(MaxTheta));
350 sintheta = std::sqrt(1. - costheta*costheta);
352 rndm2 = angRndm->GenRandPhi();
353 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
354 sinphi = std::sin(Phi);
355 cosphi = std::cos(Phi);
357 px = -sintheta * cosphi;
358 py = -sintheta * sinphi;
366 if (posDist->GetSourcePosType() ==
"Point"
367 || posDist->GetSourcePosType() ==
"Volume")
373 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
374 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
375 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
390 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
391 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
392 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
396 finx = (px*posDist->GetSideRefVec1().x())
397 + (py*posDist->GetSideRefVec2().x())
398 + (pz*posDist->GetSideRefVec3().x());
399 finy = (px*posDist->GetSideRefVec1().y())
400 + (py*posDist->GetSideRefVec2().y())
401 + (pz*posDist->GetSideRefVec3().y());
402 finz = (px*posDist->GetSideRefVec1().z())
403 + (py*posDist->GetSideRefVec2().z())
404 + (pz*posDist->GetSideRefVec3().z());
407 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
418 if(verbosityLevel >= 1)
420 G4cout <<
"Generating isotropic vector: " << mom <<
G4endl;
431 G4double sintheta, sinphi,costheta,cosphi;
432 rndm = angRndm->GenRandTheta();
433 sintheta = std::sqrt( rndm * (std::sin(MaxTheta)*std::sin(MaxTheta)
434 - std::sin(MinTheta)*std::sin(MinTheta) )
435 + std::sin(MinTheta)*std::sin(MinTheta) );
436 costheta = std::sqrt(1. -sintheta*sintheta);
438 rndm2 = angRndm->GenRandPhi();
439 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
440 sinphi = std::sin(Phi);
441 cosphi = std::cos(Phi);
443 px = -sintheta * cosphi;
444 py = -sintheta * sinphi;
452 if (posDist->GetSourcePosType() ==
"Point"
453 || posDist->GetSourcePosType() ==
"Volume")
458 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
459 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
460 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
474 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
475 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
476 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
480 finx = (px*posDist->GetSideRefVec1().x())
481 + (py*posDist->GetSideRefVec2().x())
482 + (pz*posDist->GetSideRefVec3().x());
483 finy = (px*posDist->GetSideRefVec1().y())
484 + (py*posDist->GetSideRefVec2().y())
485 + (pz*posDist->GetSideRefVec3().y());
486 finz = (px*posDist->GetSideRefVec1().z())
487 + (py*posDist->GetSideRefVec2().z())
488 + (pz*posDist->GetSideRefVec3().z());
491 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
502 if(verbosityLevel >= 1)
504 G4cout <<
"Resultant cosine-law unit momentum vector " << mom <<
G4endl;
514 if(verbosityLevel >= 1)
516 G4cout <<
"Resultant Planar wave momentum vector " << mom <<
G4endl;
524 if(UserDistType ==
"NULL")
528 else if(UserDistType ==
"theta")
531 while(Theta > MaxTheta || Theta < MinTheta)
533 Theta = GenerateUserDefTheta();
536 while(Phi > MaxPhi || Phi < MinPhi)
538 rndm = angRndm->GenRandPhi();
542 else if(UserDistType ==
"phi")
545 while(Theta > MaxTheta || Theta < MinTheta)
547 rndm = angRndm->GenRandTheta();
548 Theta = std::acos(1. - (2. * rndm));
551 while(Phi > MaxPhi || Phi < MinPhi)
553 Phi = GenerateUserDefPhi();
556 else if(UserDistType ==
"both")
559 while(Theta > MaxTheta || Theta < MinTheta)
561 Theta = GenerateUserDefTheta();
564 while(Phi > MaxPhi || Phi < MinPhi)
566 Phi = GenerateUserDefPhi();
569 px = -std::sin(Theta) * std::cos(Phi);
570 py = -std::sin(Theta) * std::sin(Phi);
571 pz = -std::cos(Theta);
573 pmag = std::sqrt((px*px) + (py*py) + (pz*pz));
582 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
583 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
584 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
592 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
606 if(verbosityLevel > 1)
608 G4cout <<
"SideRefVecs " << posDist->GetSideRefVec1()
609 << posDist->GetSideRefVec2() << posDist->GetSideRefVec3()
611 G4cout <<
"Raw Unit vector " << pxh
612 <<
"," << pyh <<
"," << pzh <<
G4endl;
614 G4double resultx = (pxh*posDist->GetSideRefVec1().x())
615 + (pyh*posDist->GetSideRefVec2().x())
616 + (pzh*posDist->GetSideRefVec3().x());
618 G4double resulty = (pxh*posDist->GetSideRefVec1().y())
619 + (pyh*posDist->GetSideRefVec2().y())
620 + (pzh*posDist->GetSideRefVec3().y());
622 G4double resultz = (pxh*posDist->GetSideRefVec1().z())
623 + (pyh*posDist->GetSideRefVec2().z())
624 + (pzh*posDist->GetSideRefVec3().z());
626 G4double ResMag = std::sqrt((resultx*resultx)
628 + (resultz*resultz));
629 resultx = resultx/ResMag;
630 resulty = resulty/ResMag;
631 resultz = resultz/ResMag;
640 if(verbosityLevel > 0 )
642 G4cout <<
"Final User Defined momentum vector "
643 << particle_momentum_direction <<
G4endl;
647G4double G4SPSAngDistribution::GenerateUserDefTheta()
652 if(UserDistType ==
"NULL" || UserDistType ==
"phi")
667 G4double bins[1024],vals[1024], sum;
669 auto maxbin =
G4int(UDefThetaH.GetVectorLength());
670 bins[0] = UDefThetaH.GetLowEdgeEnergy(std::size_t(0));
671 vals[0] = UDefThetaH(std::size_t(0));
673 for(ii=1; ii<maxbin; ++ii)
675 bins[ii] = UDefThetaH.GetLowEdgeEnergy(std::size_t(ii));
676 vals[ii] = UDefThetaH(std::size_t(ii)) + vals[ii-1];
677 sum = sum + UDefThetaH(std::size_t(ii));
679 for(ii=0; ii<maxbin; ++ii)
681 vals[ii] = vals[ii]/sum;
682 IPDFThetaH.InsertValues(bins[ii], vals[ii]);
684 IPDFThetaExist =
true;
691 return IPDFThetaH.GetEnergy(rndm);
694G4double G4SPSAngDistribution::GenerateUserDefPhi()
699 if(UserDistType ==
"NULL" || UserDistType ==
"theta")
714 G4double bins[1024],vals[1024], sum;
716 auto maxbin =
G4int(UDefPhiH.GetVectorLength());
717 bins[0] = UDefPhiH.GetLowEdgeEnergy(std::size_t(0));
718 vals[0] = UDefPhiH(std::size_t(0));
720 for(ii=1; ii<maxbin; ++ii)
722 bins[ii] = UDefPhiH.GetLowEdgeEnergy(std::size_t(ii));
723 vals[ii] = UDefPhiH(std::size_t(ii)) + vals[ii-1];
724 sum = sum + UDefPhiH(std::size_t(ii));
726 for(ii=0; ii<maxbin; ++ii)
728 vals[ii] = vals[ii]/sum;
729 IPDFPhiH.InsertValues(bins[ii], vals[ii]);
738 return IPDFPhiH.GetEnergy(rndm);
744 if (atype ==
"theta")
746 UDefThetaH = IPDFThetaH = ZeroPhysVector ;
747 IPDFThetaExist = false ;
749 else if (atype ==
"phi")
751 UDefPhiH = IPDFPhiH = ZeroPhysVector ;
752 IPDFPhiExist = false ;
768 if(AngDistType ==
"iso")
769 GenerateIsotropicFlux(localM);
770 else if(AngDistType ==
"cos")
771 GenerateCosineLawFlux(localM);
772 else if(AngDistType ==
"planar")
773 GeneratePlanarFlux(localM);
774 else if(AngDistType ==
"beam1d" || AngDistType ==
"beam2d" )
775 GenerateBeamFlux(localM);
776 else if(AngDistType ==
"user")
777 GenerateUserDefFlux(localM);
778 else if(AngDistType ==
"focused")
779 GenerateFocusedFlux(localM);
781 G4cout <<
"Error: AngDistType has unusual value" <<
G4endl;
G4TemplateAutoLock< G4Mutex > G4AutoLock
G4ThreeVector G4ParticleMomentum
#define G4MUTEXDESTROY(mutex)
#define G4MUTEXINIT(mutex)
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cout
Hep3Vector cross(const Hep3Vector &) const
void SetBeamSigmaInAngX(G4double)
void SetBeamSigmaInAngR(G4double)
G4ThreeVector GetDirection()
void SetPosDistribution(G4SPSPosDistribution *a)
void UserDefAngTheta(const G4ThreeVector &)
void SetFocusPoint(const G4ThreeVector &)
void SetAngDistType(const G4String &)
void UserDefAngPhi(const G4ThreeVector &)
G4ParticleMomentum GenerateOne()
void SetMaxTheta(G4double)
void SetUseUserAngAxis(G4bool)
void SetMinTheta(G4double)
void SetParticleMomentumDirection(const G4ParticleMomentum &aMomDirection)
void SetBeamSigmaInAngY(G4double)
void ReSetHist(const G4String &)
void SetBiasRndm(G4SPSRandomGenerator *a)
void SetVerbosity(G4int a)
void DefineAngRefAxes(const G4String &, const G4ThreeVector &)
void SetUserWRTSurface(G4bool)
DLL_API const Hep3Vector HepZHat
DLL_API const Hep3Vector HepXHat
DLL_API const Hep3Vector HepYHat