61 AngDistType =
"planar";
73 UserDistType =
"NULL";
74 UserWRTSurface =
true;
76 IPDFThetaExist =
false;
87 if(atype !=
"iso" && atype !=
"cos" && atype !=
"user" && atype !=
"planar"
88 && atype !=
"beam1d" && atype !=
"beam2d" && atype !=
"focused")
89 G4cout <<
"Error, distribution must be iso, cos, planar, beam1d, beam2d, focused or user" <<
G4endl;
92 if (AngDistType ==
"cos") MaxTheta =
pi/2. ;
93 if (AngDistType ==
"user") {
94 UDefThetaH = IPDFThetaH = ZeroPhysVector ;
95 IPDFThetaExist = false ;
96 UDefPhiH = IPDFPhiH = ZeroPhysVector ;
97 IPDFPhiExist = false ;
103 if(refname ==
"angref1")
104 AngRef1 = ref.
unit();
105 else if(refname ==
"angref2")
106 AngRef2 = ref.
unit();
113 AngRef3 = AngRef1.
cross(AngRef2);
114 AngRef2 = AngRef3.
cross(AngRef1);
116 if(verbosityLevel == 2)
118 G4cout <<
"Angular distribution rotation axes " << AngRef1 <<
" " << AngRef2 <<
" " << AngRef3 <<
G4endl;
159 if(UserDistType ==
"NULL") UserDistType =
"theta";
160 if(UserDistType ==
"phi") UserDistType =
"both";
164 if(verbosityLevel >= 1)
171 if(UserDistType ==
"NULL") UserDistType =
"phi";
172 if(UserDistType ==
"theta") UserDistType =
"both";
176 if(verbosityLevel >= 1)
193 UserWRTSurface = wrtSurf;
200 UserAngRef = userang;
203 void G4SPSAngDistribution::GenerateBeamFlux()
207 if (AngDistType ==
"beam1d")
209 theta = G4RandGauss::shoot(0.0,DR);
214 px = G4RandGauss::shoot(0.0,DX);
215 py = G4RandGauss::shoot(0.0,DY);
216 theta = std::sqrt (px*px + py*py);
218 phi = std::acos(px/theta);
219 if ( py < 0.) phi = -phi;
226 px = -std::sin(theta) * std::cos(phi);
227 py = -std::sin(theta) * std::sin(phi);
228 pz = -std::cos(theta);
230 finx = px, finy =py, finz =pz;
234 finx = (px * AngRef1.
x()) + (py * AngRef2.
x()) + (pz * AngRef3.
x());
235 finy = (px * AngRef1.
y()) + (py * AngRef2.
y()) + (pz * AngRef3.
y());
236 finz = (px * AngRef1.
z()) + (py * AngRef2.
z()) + (pz * AngRef3.
z());
237 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
242 particle_momentum_direction.
setX(finx);
243 particle_momentum_direction.
setY(finy);
244 particle_momentum_direction.
setZ(finz);
247 if(verbosityLevel >= 1)
248 G4cout <<
"Generating beam vector: " << particle_momentum_direction <<
G4endl;
251 void G4SPSAngDistribution::GenerateFocusedFlux()
253 particle_momentum_direction = (FocusPoint - posDist->particle_position).unit();
256 if(verbosityLevel >= 1)
257 G4cout <<
"Generating focused vector: " << particle_momentum_direction <<
G4endl;
260 void G4SPSAngDistribution::GenerateIsotropicFlux()
268 G4double sintheta, sinphi,costheta,cosphi;
270 costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta) - std::cos(MaxTheta));
271 sintheta = std::sqrt(1. - costheta*costheta);
274 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
275 sinphi = std::sin(Phi);
276 cosphi = std::cos(Phi);
278 px = -sintheta * cosphi;
279 py = -sintheta * sinphi;
286 if (posDist->SourcePosType ==
"Point" || posDist->SourcePosType ==
"Volume") {
290 finx = (px * AngRef1.
x()) + (py * AngRef2.
x()) + (pz * AngRef3.
x());
291 finy = (px * AngRef1.
y()) + (py * AngRef2.
y()) + (pz * AngRef3.
y());
292 finz = (px * AngRef1.
z()) + (py * AngRef2.
z()) + (pz * AngRef3.
z());
302 finx = (px * AngRef1.
x()) + (py * AngRef2.
x()) + (pz * AngRef3.
x());
303 finy = (px * AngRef1.
y()) + (py * AngRef2.
y()) + (pz * AngRef3.
y());
304 finz = (px * AngRef1.
z()) + (py * AngRef2.
z()) + (pz * AngRef3.
z());
306 finx = (px*posDist->SideRefVec1.
x()) + (py*posDist->SideRefVec2.
x()) + (pz*posDist->SideRefVec3.
x());
307 finy = (px*posDist->SideRefVec1.
y()) + (py*posDist->SideRefVec2.
y()) + (pz*posDist->SideRefVec3.
y());
308 finz = (px*posDist->SideRefVec1.
z()) + (py*posDist->SideRefVec2.
z()) + (pz*posDist->SideRefVec3.
z());
311 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
316 particle_momentum_direction.
setX(finx);
317 particle_momentum_direction.
setY(finy);
318 particle_momentum_direction.
setZ(finz);
321 if(verbosityLevel >= 1)
322 G4cout <<
"Generating isotropic vector: " << particle_momentum_direction <<
G4endl;
325 void G4SPSAngDistribution::GenerateCosineLawFlux()
331 G4double sintheta, sinphi,costheta,cosphi;
333 sintheta = std::sqrt( rndm * (std::sin(MaxTheta)*std::sin(MaxTheta) - std::sin(MinTheta)*std::sin(MinTheta) )
334 +std::sin(MinTheta)*std::sin(MinTheta) );
335 costheta = std::sqrt(1. -sintheta*sintheta);
338 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
339 sinphi = std::sin(Phi);
340 cosphi = std::cos(Phi);
342 px = -sintheta * cosphi;
343 py = -sintheta * sinphi;
350 if (posDist->SourcePosType ==
"Point" || posDist->SourcePosType ==
"Volume") {
353 finx = (px * AngRef1.
x()) + (py * AngRef2.
x()) + (pz * AngRef3.
x());
354 finy = (px * AngRef1.
y()) + (py * AngRef2.
y()) + (pz * AngRef3.
y());
355 finz = (px * AngRef1.
z()) + (py * AngRef2.
z()) + (pz * AngRef3.
z());
364 finx = (px * AngRef1.
x()) + (py * AngRef2.
x()) + (pz * AngRef3.
x());
365 finy = (px * AngRef1.
y()) + (py * AngRef2.
y()) + (pz * AngRef3.
y());
366 finz = (px * AngRef1.
z()) + (py * AngRef2.
z()) + (pz * AngRef3.
z());
368 finx = (px*posDist->SideRefVec1.
x()) + (py*posDist->SideRefVec2.
x()) + (pz*posDist->SideRefVec3.
x());
369 finy = (px*posDist->SideRefVec1.
y()) + (py*posDist->SideRefVec2.
y()) + (pz*posDist->SideRefVec3.
y());
370 finz = (px*posDist->SideRefVec1.
z()) + (py*posDist->SideRefVec2.
z()) + (pz*posDist->SideRefVec3.
z());
373 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
378 particle_momentum_direction.
setX(finx);
379 particle_momentum_direction.
setY(finy);
380 particle_momentum_direction.
setZ(finz);
383 if(verbosityLevel >= 1)
385 G4cout <<
"Resultant cosine-law unit momentum vector " << particle_momentum_direction <<
G4endl;
389 void G4SPSAngDistribution::GeneratePlanarFlux()
394 if(verbosityLevel >= 1)
396 G4cout <<
"Resultant Planar wave momentum vector " << particle_momentum_direction <<
G4endl;
400 void G4SPSAngDistribution::GenerateUserDefFlux()
404 if(UserDistType ==
"NULL")
406 else if(UserDistType ==
"theta") {
408 while(Theta > MaxTheta || Theta < MinTheta)
409 Theta = GenerateUserDefTheta();
411 while(Phi > MaxPhi || Phi < MinPhi) {
416 else if(UserDistType ==
"phi") {
418 while(Theta > MaxTheta || Theta < MinTheta)
421 Theta = std::acos(1. - (2. * rndm));
424 while(Phi > MaxPhi || Phi < MinPhi)
425 Phi = GenerateUserDefPhi();
427 else if(UserDistType ==
"both")
430 while(Theta > MaxTheta || Theta < MinTheta)
431 Theta = GenerateUserDefTheta();
433 while(Phi > MaxPhi || Phi < MinPhi)
434 Phi = GenerateUserDefPhi();
436 px = -std::sin(Theta) * std::cos(Phi);
437 py = -std::sin(Theta) * std::sin(Phi);
438 pz = -std::cos(Theta);
440 pmag = std::sqrt((px*px) + (py*py) + (pz*pz));
442 if(!UserWRTSurface) {
447 finx = (px * AngRef1.
x()) + (py * AngRef2.
x()) + (pz * AngRef3.
x());
448 finy = (px * AngRef1.
y()) + (py * AngRef2.
y()) + (pz * AngRef3.
y());
449 finz = (px * AngRef1.
z()) + (py * AngRef2.
z()) + (pz * AngRef3.
z());
455 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
460 particle_momentum_direction.
setX(finx);
461 particle_momentum_direction.
setY(finy);
462 particle_momentum_direction.
setZ(finz);
468 if(verbosityLevel > 1) {
469 G4cout <<
"SideRefVecs " <<posDist->SideRefVec1<<posDist->SideRefVec2<<posDist->SideRefVec3<<
G4endl;
470 G4cout <<
"Raw Unit vector "<<pxh<<
","<<pyh<<
","<<pzh<<
G4endl;
472 G4double resultx = (pxh*posDist->SideRefVec1.
x()) + (pyh*posDist->SideRefVec2.
x()) +
473 (pzh*posDist->SideRefVec3.
x());
475 G4double resulty = (pxh*posDist->SideRefVec1.
y()) + (pyh*posDist->SideRefVec2.
y()) +
476 (pzh*posDist->SideRefVec3.
y());
478 G4double resultz = (pxh*posDist->SideRefVec1.
z()) + (pyh*posDist->SideRefVec2.
z()) +
479 (pzh*posDist->SideRefVec3.
z());
481 G4double ResMag = std::sqrt((resultx*resultx) + (resulty*resulty) + (resultz*resultz));
482 resultx = resultx/ResMag;
483 resulty = resulty/ResMag;
484 resultz = resultz/ResMag;
486 particle_momentum_direction.
setX(resultx);
487 particle_momentum_direction.
setY(resulty);
488 particle_momentum_direction.
setZ(resultz);
492 if(verbosityLevel > 0 )
494 G4cout <<
"Final User Defined momentum vector " << particle_momentum_direction <<
G4endl;
498 G4double G4SPSAngDistribution::GenerateUserDefTheta()
502 if(UserDistType ==
"NULL" || UserDistType ==
"phi")
513 if(IPDFThetaExist ==
false)
516 G4double bins[1024],vals[1024], sum;
520 vals[0] = UDefThetaH(
size_t(0));
522 for(ii=1;ii<maxbin;ii++)
525 vals[ii] = UDefThetaH(
size_t(ii)) + vals[ii-1];
526 sum = sum + UDefThetaH(
size_t(ii));
528 for(ii=0;ii<maxbin;ii++)
530 vals[ii] = vals[ii]/sum;
534 IPDFThetaExist =
true;
542 G4double G4SPSAngDistribution::GenerateUserDefPhi()
547 if(UserDistType ==
"NULL" || UserDistType ==
"theta")
558 if(IPDFPhiExist ==
false)
561 G4double bins[1024],vals[1024], sum;
565 vals[0] = UDefPhiH(
size_t(0));
567 for(ii=1;ii<maxbin;ii++)
570 vals[ii] = UDefPhiH(
size_t(ii)) + vals[ii-1];
571 sum = sum + UDefPhiH(
size_t(ii));
574 for(ii=0;ii<maxbin;ii++)
576 vals[ii] = vals[ii]/sum;
590 if (atype ==
"theta") {
591 UDefThetaH = IPDFThetaH = ZeroPhysVector ;
592 IPDFThetaExist = false ;}
593 else if (atype ==
"phi"){
594 UDefPhiH = IPDFPhiH = ZeroPhysVector ;
595 IPDFPhiExist = false ;}
605 if(AngDistType ==
"iso")
606 GenerateIsotropicFlux();
607 else if(AngDistType ==
"cos")
608 GenerateCosineLawFlux();
609 else if(AngDistType ==
"planar")
610 GeneratePlanarFlux();
611 else if(AngDistType ==
"beam1d" || AngDistType ==
"beam2d" )
613 else if(AngDistType ==
"user")
614 GenerateUserDefFlux();
615 else if(AngDistType ==
"focused")
616 GenerateFocusedFlux();
618 G4cout <<
"Error: AngDistType has unusual value" <<
G4endl;
619 return particle_momentum_direction;