Isotropic vector generation

[nkx111]

restG4/src/PrimaryGeneratorAction.cc

Line 671 in 8e6c7ef

G4ThreeVector PrimaryGeneratorAction::GetIsotropicVector() {

Generating vector in this way may not be isotropic. Shouldn't we use θ and Φ to randomize direction?

[jgalan]

I believe the present implementation should be isotropic. Using angles would be more complex implementation

Here is a good example
http://corysimon.github.io/articles/uniformdistn-on-sphere/

Do you have any clues pointing to a problem on isotropy with present vector?

[jgalan]

Our present implementation is closer to the "Alternative method 1", however, we do not protect normalization by very small vectors, as they do, this could cause problems for us:

They require that the norm is larger than 0.001, we could improve that in our method.

[lobis]

I think we should really use the internal features of geant4 for this, using the G4GeneralParticleSource. it will also work with the random seed internally. I will update this function as an example to see the difference it makes, I think it will greatly simplify our code.

We can only use the angular distribution from the generator, we do not need to change from G4ParticleGun.

I will do this later this morning.

[nkx111]

The currecnt method is by throwing points randomly in a cube, and calculate their directions. I guess there will be more events near the 8 corners of the cube

[nkx111]

OK, I see we rejected the events with n >1. Then this is right. Isotropic indeed. Just be less elegant.

[jgalan]

Yes, you always get a uniform distribution within a sphere of radius 1. Indeed, efficiency goes with 4/3 * PI * (D/2)^3/D^3. Volume of the sphere versus volume of the cube.

However, the calculation of cosines and other mathematical functions required in other methods believe is more expensive computationally.

[jgalan]

Ok, finally I believe there is a bug there, because the random number is running between -0.5 and 0.5, so the radius is 0.5 and not 1!

[jgalan]

No no, it is ok, it is ok, it is divided by 0.5 so it is runnning between -1 and 1. :)

[lobis]

This generator is improved in commit (3826fd1) of the branch lobis-patch. There you can see this new approach, which I think is much better.

I will do the pull request once it is better organized, feel free to discuss the approach, I think it could be ported to all distributions we are currently using.

[jgalan]

It is probably more time consuming than the simple implementation we had. On top of that it loses transparency, before we were at least able to go to the method and read clearly what it was doing.

[jgalan]

void G4SPSAngDistribution::GenerateIsotropicFlux(G4ParticleMomentum& mom)
  311 {
  312   // generates isotropic flux.
  313   // No vectors are needed.
  314   G4double rndm, rndm2;
  315   G4double px, py, pz;
  316 
  317   //
  318   G4double sintheta, sinphi,costheta,cosphi;
  319   rndm = angRndm->GenRandTheta();
  320   costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta) - std::cos(MaxTheta));
  321   sintheta = std::sqrt(1. - costheta*costheta);
  322   
  323   rndm2 = angRndm->GenRandPhi();
  324   Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; 
  325   sinphi = std::sin(Phi);
  326   cosphi = std::cos(Phi);
  327 
  328   px = -sintheta * cosphi;
  329   py = -sintheta * sinphi;
  330   pz = -costheta;
  331 
  332   // for volume and ponit source use mother or user defined co-ordinates
  333   // for plane and surface source user surface-normal or userdefined co-ordinates
  334   //
  335   G4double finx, finy, finz;
  336   if (posDist->GetSourcePosType() == "Point" || posDist->GetSourcePosType() == "Volume") {
  337     if (UserAngRef){
  338       // Apply Rotation Matrix
  339       // x * AngRef1, y * AngRef2 and z * AngRef3
  340       finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
  341       finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
  342       finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
  343     } else {
  344       finx = px;
  345       finy = py;
  346       finz = pz;
  347     }
  348   } else {    // for plane and surface source   
  349     if (UserAngRef){
  350       // Apply Rotation Matrix
  351       // x * AngRef1, y * AngRef2 and z * AngRef3
  352       finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
  353       finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
  354       finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
  355     } else {
  356       finx = (px*posDist->GetSideRefVec1().x()) + (py*posDist->GetSideRefVec2().x()) + (pz*posDist->GetSideRefVec3().x());
  357       finy = (px*posDist->GetSideRefVec1().y()) + (py*posDist->GetSideRefVec2().y()) + (pz*posDist->GetSideRefVec3().y());
  358       finz = (px*posDist->GetSideRefVec1().z()) + (py*posDist->GetSideRefVec2().z()) + (pz*posDist->GetSideRefVec3().z());
  359     }
  360   }
  361   G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
  362   finx = finx/ResMag;
  363   finy = finy/ResMag;
  364   finz = finz/ResMag;
  365 
  366   mom.setX(finx);
  367   mom.setY(finy);
  368   mom.setZ(finz);
  369 
  370   // particle_momentum_direction now holds unit momentum vector.
  371   if(verbosityLevel >= 1)
  372     G4cout << "Generating isotropic vector: " << mom << G4endl;
  373 }

[lobis]

It is probably more time consuming than the simple implementation we had. On top of that it loses transparency, before we were at least able to go to the method and read clearly what it was doing.

Well, I don't think the performance is a bottleneck here, as the time generating primaries is most of the time neglectable compared to the tracking. Maybe in the worst case scenario you could save a few minutes.

Regarding the "transparency", I guess this is a matter of opinion. I am of the opinion that the less code, the better, as more code means more maintenance. Personally I think the tradeoff if justified here, but it is just an opinion.