00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00027
00028
00029
00030
00031
00032
00033
00034
00035
00037
00038
00039
00040
00041
00042
00043
00044
00045
00047
00048
00049 #include "G4SPSAngDistribution.hh"
00050
00051 #include "Randomize.hh"
00052 #include "G4PhysicalConstants.hh"
00053
00054 G4SPSAngDistribution::G4SPSAngDistribution()
00055 : Theta(0.), Phi(0.)
00056 {
00057
00058 G4ThreeVector zero;
00059 particle_momentum_direction = G4ParticleMomentum(0,0,-1);
00060
00061 AngDistType = "planar";
00062 AngRef1 = CLHEP::HepXHat;
00063 AngRef2 = CLHEP::HepYHat;
00064 AngRef3 = CLHEP::HepZHat;
00065 MinTheta = 0.;
00066 MaxTheta = pi;
00067 MinPhi = 0.;
00068 MaxPhi = twopi;
00069 DR = 0.;
00070 DX = 0.;
00071 DY = 0.;
00072 FocusPoint = G4ThreeVector(0., 0., 0.);
00073 UserDistType = "NULL";
00074 UserWRTSurface = true;
00075 UserAngRef = false;
00076 IPDFThetaExist = false;
00077 IPDFPhiExist = false;
00078 verbosityLevel = 0 ;
00079 }
00080
00081 G4SPSAngDistribution::~G4SPSAngDistribution()
00082 {}
00083
00084
00085 void G4SPSAngDistribution::SetAngDistType(G4String atype)
00086 {
00087 if(atype != "iso" && atype != "cos" && atype != "user" && atype != "planar"
00088 && atype != "beam1d" && atype != "beam2d" && atype != "focused")
00089 G4cout << "Error, distribution must be iso, cos, planar, beam1d, beam2d, focused or user" << G4endl;
00090 else
00091 AngDistType = atype;
00092 if (AngDistType == "cos") MaxTheta = pi/2. ;
00093 if (AngDistType == "user") {
00094 UDefThetaH = IPDFThetaH = ZeroPhysVector ;
00095 IPDFThetaExist = false ;
00096 UDefPhiH = IPDFPhiH = ZeroPhysVector ;
00097 IPDFPhiExist = false ;
00098 }
00099 }
00100
00101 void G4SPSAngDistribution::DefineAngRefAxes(G4String refname, G4ThreeVector ref)
00102 {
00103 if(refname == "angref1")
00104 AngRef1 = ref.unit();
00105 else if(refname == "angref2")
00106 AngRef2 = ref.unit();
00107
00108
00109
00110
00111
00112
00113 AngRef3 = AngRef1.cross(AngRef2);
00114 AngRef2 = AngRef3.cross(AngRef1);
00115 UserAngRef = true ;
00116 if(verbosityLevel == 2)
00117 {
00118 G4cout << "Angular distribution rotation axes " << AngRef1 << " " << AngRef2 << " " << AngRef3 << G4endl;
00119 }
00120 }
00121
00122 void G4SPSAngDistribution::SetMinTheta(G4double mint)
00123 {
00124 MinTheta = mint;
00125 }
00126
00127 void G4SPSAngDistribution::SetMinPhi(G4double minp)
00128 {
00129 MinPhi = minp;
00130 }
00131
00132 void G4SPSAngDistribution::SetMaxTheta(G4double maxt)
00133 {
00134 MaxTheta = maxt;
00135 }
00136
00137 void G4SPSAngDistribution::SetMaxPhi(G4double maxp)
00138 {
00139 MaxPhi = maxp;
00140 }
00141
00142 void G4SPSAngDistribution::SetBeamSigmaInAngR(G4double r)
00143 {
00144 DR = r;
00145 }
00146
00147 void G4SPSAngDistribution::SetBeamSigmaInAngX(G4double r)
00148 {
00149 DX = r;
00150 }
00151
00152 void G4SPSAngDistribution::SetBeamSigmaInAngY(G4double r)
00153 {
00154 DY = r;
00155 }
00156
00157 void G4SPSAngDistribution::UserDefAngTheta(G4ThreeVector input)
00158 {
00159 if(UserDistType == "NULL") UserDistType = "theta";
00160 if(UserDistType == "phi") UserDistType = "both";
00161 G4double thi, val;
00162 thi = input.x();
00163 val = input.y();
00164 if(verbosityLevel >= 1)
00165 G4cout << "In UserDefAngTheta" << G4endl;
00166 UDefThetaH.InsertValues(thi, val);
00167 }
00168
00169 void G4SPSAngDistribution::UserDefAngPhi(G4ThreeVector input)
00170 {
00171 if(UserDistType == "NULL") UserDistType = "phi";
00172 if(UserDistType == "theta") UserDistType = "both";
00173 G4double phhi, val;
00174 phhi = input.x();
00175 val = input.y();
00176 if(verbosityLevel >= 1)
00177 G4cout << "In UserDefAngPhi" << G4endl;
00178 UDefPhiH.InsertValues(phhi, val);
00179 }
00180
00181 void G4SPSAngDistribution::SetFocusPoint(G4ThreeVector input)
00182 {
00183 FocusPoint = input;
00184 }
00185
00186 void G4SPSAngDistribution::SetUserWRTSurface(G4bool wrtSurf)
00187 {
00188
00189
00190
00191
00192
00193 UserWRTSurface = wrtSurf;
00194 }
00195
00196 void G4SPSAngDistribution::SetUseUserAngAxis(G4bool userang)
00197 {
00198
00199
00200 UserAngRef = userang;
00201 }
00202
00203 void G4SPSAngDistribution::GenerateBeamFlux()
00204 {
00205 G4double theta, phi;
00206 G4double px, py, pz;
00207 if (AngDistType == "beam1d")
00208 {
00209 theta = G4RandGauss::shoot(0.0,DR);
00210 phi = twopi * G4UniformRand();
00211 }
00212 else
00213 {
00214 px = G4RandGauss::shoot(0.0,DX);
00215 py = G4RandGauss::shoot(0.0,DY);
00216 theta = std::sqrt (px*px + py*py);
00217 if (theta != 0.) {
00218 phi = std::acos(px/theta);
00219 if ( py < 0.) phi = -phi;
00220 }
00221 else
00222 {
00223 phi = 0.0;
00224 }
00225 }
00226 px = -std::sin(theta) * std::cos(phi);
00227 py = -std::sin(theta) * std::sin(phi);
00228 pz = -std::cos(theta);
00229 G4double finx, finy, finz ;
00230 finx = px, finy =py, finz =pz;
00231 if (UserAngRef){
00232
00233
00234 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
00235 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
00236 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
00237 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
00238 finx = finx/ResMag;
00239 finy = finy/ResMag;
00240 finz = finz/ResMag;
00241 }
00242 particle_momentum_direction.setX(finx);
00243 particle_momentum_direction.setY(finy);
00244 particle_momentum_direction.setZ(finz);
00245
00246
00247 if(verbosityLevel >= 1)
00248 G4cout << "Generating beam vector: " << particle_momentum_direction << G4endl;
00249 }
00250
00251 void G4SPSAngDistribution::GenerateFocusedFlux()
00252 {
00253 particle_momentum_direction = (FocusPoint - posDist->particle_position).unit();
00254
00255
00256 if(verbosityLevel >= 1)
00257 G4cout << "Generating focused vector: " << particle_momentum_direction << G4endl;
00258 }
00259
00260 void G4SPSAngDistribution::GenerateIsotropicFlux()
00261 {
00262
00263
00264 G4double rndm, rndm2;
00265 G4double px, py, pz;
00266
00267
00268 G4double sintheta, sinphi,costheta,cosphi;
00269 rndm = angRndm->GenRandTheta();
00270 costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta) - std::cos(MaxTheta));
00271 sintheta = std::sqrt(1. - costheta*costheta);
00272
00273 rndm2 = angRndm->GenRandPhi();
00274 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
00275 sinphi = std::sin(Phi);
00276 cosphi = std::cos(Phi);
00277
00278 px = -sintheta * cosphi;
00279 py = -sintheta * sinphi;
00280 pz = -costheta;
00281
00282
00283
00284
00285 G4double finx, finy, finz;
00286 if (posDist->SourcePosType == "Point" || posDist->SourcePosType == "Volume") {
00287 if (UserAngRef){
00288
00289
00290 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
00291 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
00292 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
00293 } else {
00294 finx = px;
00295 finy = py;
00296 finz = pz;
00297 }
00298 } else {
00299 if (UserAngRef){
00300
00301
00302 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
00303 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
00304 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
00305 } else {
00306 finx = (px*posDist->SideRefVec1.x()) + (py*posDist->SideRefVec2.x()) + (pz*posDist->SideRefVec3.x());
00307 finy = (px*posDist->SideRefVec1.y()) + (py*posDist->SideRefVec2.y()) + (pz*posDist->SideRefVec3.y());
00308 finz = (px*posDist->SideRefVec1.z()) + (py*posDist->SideRefVec2.z()) + (pz*posDist->SideRefVec3.z());
00309 }
00310 }
00311 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
00312 finx = finx/ResMag;
00313 finy = finy/ResMag;
00314 finz = finz/ResMag;
00315
00316 particle_momentum_direction.setX(finx);
00317 particle_momentum_direction.setY(finy);
00318 particle_momentum_direction.setZ(finz);
00319
00320
00321 if(verbosityLevel >= 1)
00322 G4cout << "Generating isotropic vector: " << particle_momentum_direction << G4endl;
00323 }
00324
00325 void G4SPSAngDistribution::GenerateCosineLawFlux()
00326 {
00327
00328 G4double px, py, pz;
00329 G4double rndm, rndm2;
00330
00331 G4double sintheta, sinphi,costheta,cosphi;
00332 rndm = angRndm->GenRandTheta();
00333 sintheta = std::sqrt( rndm * (std::sin(MaxTheta)*std::sin(MaxTheta) - std::sin(MinTheta)*std::sin(MinTheta) )
00334 +std::sin(MinTheta)*std::sin(MinTheta) );
00335 costheta = std::sqrt(1. -sintheta*sintheta);
00336
00337 rndm2 = angRndm->GenRandPhi();
00338 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
00339 sinphi = std::sin(Phi);
00340 cosphi = std::cos(Phi);
00341
00342 px = -sintheta * cosphi;
00343 py = -sintheta * sinphi;
00344 pz = -costheta;
00345
00346
00347
00348
00349 G4double finx, finy, finz;
00350 if (posDist->SourcePosType == "Point" || posDist->SourcePosType == "Volume") {
00351 if (UserAngRef){
00352
00353 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
00354 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
00355 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
00356 } else {
00357 finx = px;
00358 finy = py;
00359 finz = pz;
00360 }
00361 } else {
00362 if (UserAngRef){
00363
00364 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
00365 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
00366 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
00367 } else {
00368 finx = (px*posDist->SideRefVec1.x()) + (py*posDist->SideRefVec2.x()) + (pz*posDist->SideRefVec3.x());
00369 finy = (px*posDist->SideRefVec1.y()) + (py*posDist->SideRefVec2.y()) + (pz*posDist->SideRefVec3.y());
00370 finz = (px*posDist->SideRefVec1.z()) + (py*posDist->SideRefVec2.z()) + (pz*posDist->SideRefVec3.z());
00371 }
00372 }
00373 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
00374 finx = finx/ResMag;
00375 finy = finy/ResMag;
00376 finz = finz/ResMag;
00377
00378 particle_momentum_direction.setX(finx);
00379 particle_momentum_direction.setY(finy);
00380 particle_momentum_direction.setZ(finz);
00381
00382
00383 if(verbosityLevel >= 1)
00384 {
00385 G4cout << "Resultant cosine-law unit momentum vector " << particle_momentum_direction << G4endl;
00386 }
00387 }
00388
00389 void G4SPSAngDistribution::GeneratePlanarFlux()
00390 {
00391
00392
00393
00394 if(verbosityLevel >= 1)
00395 {
00396 G4cout << "Resultant Planar wave momentum vector " << particle_momentum_direction << G4endl;
00397 }
00398 }
00399
00400 void G4SPSAngDistribution::GenerateUserDefFlux()
00401 {
00402 G4double rndm, px, py, pz, pmag;
00403
00404 if(UserDistType == "NULL")
00405 G4cout << "Error: UserDistType undefined" << G4endl;
00406 else if(UserDistType == "theta") {
00407 Theta = 10.;
00408 while(Theta > MaxTheta || Theta < MinTheta)
00409 Theta = GenerateUserDefTheta();
00410 Phi = 10.;
00411 while(Phi > MaxPhi || Phi < MinPhi) {
00412 rndm = angRndm->GenRandPhi();
00413 Phi = twopi * rndm;
00414 }
00415 }
00416 else if(UserDistType == "phi") {
00417 Theta = 10.;
00418 while(Theta > MaxTheta || Theta < MinTheta)
00419 {
00420 rndm = angRndm->GenRandTheta();
00421 Theta = std::acos(1. - (2. * rndm));
00422 }
00423 Phi = 10.;
00424 while(Phi > MaxPhi || Phi < MinPhi)
00425 Phi = GenerateUserDefPhi();
00426 }
00427 else if(UserDistType == "both")
00428 {
00429 Theta = 10.;
00430 while(Theta > MaxTheta || Theta < MinTheta)
00431 Theta = GenerateUserDefTheta();
00432 Phi = 10.;
00433 while(Phi > MaxPhi || Phi < MinPhi)
00434 Phi = GenerateUserDefPhi();
00435 }
00436 px = -std::sin(Theta) * std::cos(Phi);
00437 py = -std::sin(Theta) * std::sin(Phi);
00438 pz = -std::cos(Theta);
00439
00440 pmag = std::sqrt((px*px) + (py*py) + (pz*pz));
00441
00442 if(!UserWRTSurface) {
00443 G4double finx, finy, finz;
00444 if (UserAngRef) {
00445
00446
00447 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
00448 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
00449 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
00450 } else {
00451 finx = px;
00452 finy = py;
00453 finz = pz;
00454 }
00455 G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
00456 finx = finx/ResMag;
00457 finy = finy/ResMag;
00458 finz = finz/ResMag;
00459
00460 particle_momentum_direction.setX(finx);
00461 particle_momentum_direction.setY(finy);
00462 particle_momentum_direction.setZ(finz);
00463 }
00464 else {
00465 G4double pxh = px/pmag;
00466 G4double pyh = py/pmag;
00467 G4double pzh = pz/pmag;
00468 if(verbosityLevel > 1) {
00469 G4cout <<"SideRefVecs " <<posDist->SideRefVec1<<posDist->SideRefVec2<<posDist->SideRefVec3<<G4endl;
00470 G4cout <<"Raw Unit vector "<<pxh<<","<<pyh<<","<<pzh<<G4endl;
00471 }
00472 G4double resultx = (pxh*posDist->SideRefVec1.x()) + (pyh*posDist->SideRefVec2.x()) +
00473 (pzh*posDist->SideRefVec3.x());
00474
00475 G4double resulty = (pxh*posDist->SideRefVec1.y()) + (pyh*posDist->SideRefVec2.y()) +
00476 (pzh*posDist->SideRefVec3.y());
00477
00478 G4double resultz = (pxh*posDist->SideRefVec1.z()) + (pyh*posDist->SideRefVec2.z()) +
00479 (pzh*posDist->SideRefVec3.z());
00480
00481 G4double ResMag = std::sqrt((resultx*resultx) + (resulty*resulty) + (resultz*resultz));
00482 resultx = resultx/ResMag;
00483 resulty = resulty/ResMag;
00484 resultz = resultz/ResMag;
00485
00486 particle_momentum_direction.setX(resultx);
00487 particle_momentum_direction.setY(resulty);
00488 particle_momentum_direction.setZ(resultz);
00489 }
00490
00491
00492 if(verbosityLevel > 0 )
00493 {
00494 G4cout << "Final User Defined momentum vector " << particle_momentum_direction << G4endl;
00495 }
00496 }
00497
00498 G4double G4SPSAngDistribution::GenerateUserDefTheta()
00499 {
00500
00501
00502 if(UserDistType == "NULL" || UserDistType == "phi")
00503 {
00504
00505 G4cout << "Error ***********************" << G4endl;
00506 G4cout << "UserDistType = " << UserDistType << G4endl;
00507 return (0.);
00508 }
00509 else
00510 {
00511
00512
00513 if(IPDFThetaExist == false)
00514 {
00515
00516 G4double bins[1024],vals[1024], sum;
00517 G4int ii;
00518 G4int maxbin = G4int(UDefThetaH.GetVectorLength());
00519 bins[0] = UDefThetaH.GetLowEdgeEnergy(size_t(0));
00520 vals[0] = UDefThetaH(size_t(0));
00521 sum = vals[0];
00522 for(ii=1;ii<maxbin;ii++)
00523 {
00524 bins[ii] = UDefThetaH.GetLowEdgeEnergy(size_t(ii));
00525 vals[ii] = UDefThetaH(size_t(ii)) + vals[ii-1];
00526 sum = sum + UDefThetaH(size_t(ii));
00527 }
00528 for(ii=0;ii<maxbin;ii++)
00529 {
00530 vals[ii] = vals[ii]/sum;
00531 IPDFThetaH.InsertValues(bins[ii], vals[ii]);
00532 }
00533
00534 IPDFThetaExist = true;
00535 }
00536
00537 G4double rndm = G4UniformRand();
00538 return(IPDFThetaH.GetEnergy(rndm));
00539 }
00540 }
00541
00542 G4double G4SPSAngDistribution::GenerateUserDefPhi()
00543 {
00544
00545
00546
00547 if(UserDistType == "NULL" || UserDistType == "theta")
00548 {
00549
00550 G4cout << "Error ***********************" << G4endl;
00551 G4cout << "UserDistType = " << UserDistType << G4endl;
00552 return(0.);
00553 }
00554 else
00555 {
00556
00557
00558 if(IPDFPhiExist == false)
00559 {
00560
00561 G4double bins[1024],vals[1024], sum;
00562 G4int ii;
00563 G4int maxbin = G4int(UDefPhiH.GetVectorLength());
00564 bins[0] = UDefPhiH.GetLowEdgeEnergy(size_t(0));
00565 vals[0] = UDefPhiH(size_t(0));
00566 sum = vals[0];
00567 for(ii=1;ii<maxbin;ii++)
00568 {
00569 bins[ii] = UDefPhiH.GetLowEdgeEnergy(size_t(ii));
00570 vals[ii] = UDefPhiH(size_t(ii)) + vals[ii-1];
00571 sum = sum + UDefPhiH(size_t(ii));
00572 }
00573
00574 for(ii=0;ii<maxbin;ii++)
00575 {
00576 vals[ii] = vals[ii]/sum;
00577 IPDFPhiH.InsertValues(bins[ii], vals[ii]);
00578 }
00579
00580 IPDFPhiExist = true;
00581 }
00582
00583 G4double rndm = G4UniformRand();
00584 return(IPDFPhiH.GetEnergy(rndm));
00585 }
00586 }
00587
00588 void G4SPSAngDistribution::ReSetHist(G4String atype)
00589 {
00590 if (atype == "theta") {
00591 UDefThetaH = IPDFThetaH = ZeroPhysVector ;
00592 IPDFThetaExist = false ;}
00593 else if (atype == "phi"){
00594 UDefPhiH = IPDFPhiH = ZeroPhysVector ;
00595 IPDFPhiExist = false ;}
00596 else {
00597 G4cout << "Error, histtype not accepted " << G4endl;
00598 }
00599 }
00600
00601
00602 G4ParticleMomentum G4SPSAngDistribution::GenerateOne()
00603 {
00604
00605 if(AngDistType == "iso")
00606 GenerateIsotropicFlux();
00607 else if(AngDistType == "cos")
00608 GenerateCosineLawFlux();
00609 else if(AngDistType == "planar")
00610 GeneratePlanarFlux();
00611 else if(AngDistType == "beam1d" || AngDistType == "beam2d" )
00612 GenerateBeamFlux();
00613 else if(AngDistType == "user")
00614 GenerateUserDefFlux();
00615 else if(AngDistType == "focused")
00616 GenerateFocusedFlux();
00617 else
00618 G4cout << "Error: AngDistType has unusual value" << G4endl;
00619 return particle_momentum_direction;
00620 }
00621
00622
00623
00624
00625
00626
00627
00628
00629