Problem of generate random 4pi vectors

[Qi_Zhang]

Please fill out the following information to help in answering your question, and also see tips for posting code snippets. If you don’t provide this information it will take more time to help with your problem!

Geant4 Version:
Operating System:
Compiler/Version:
CMake Version:

—Hi,
I hope to generate random 4pi vectors for the particles in simulation. The code is following:

 G4double  theta_m=G4UniformRand()*1.0*CLHEP::pi;
  G4double  phi_m=2.0*CLHEP::pi*G4UniformRand();
G4ThreeVector vec = G4ThreeVector(0.,0.,1.).rotateY(theta_m);                                                                                                               
  vec = vec.rotateZ(phi_m); 
G4double  pos_num=G4UniformRand();
G4ThreeVector  pos;
 if(pos_num<0.5){ pos=G4ThreeVector(1.5*cm, 1.5*cm, 1.5*cm);
   }else{pos=G4ThreeVector(-1.5*cm,-1.5*cm,-1.5*cm);                                                                                                                              
   }


particleGun->SetParticlePosition(pos);
 particleGun->SetParticleMomentumDirection(vec);

And I placed 4 plan detectors in different positions (0,50,0),(0,-50,0),(0,0,50),(0,0,-50). However I found some results are wrong!

From the histograms, the results in (0,-50,0) and (0,0,50) are wrong! Is this one bug of Geant4?
How can I fix it in my codes? Waiting for your help ! Thanks!

[evc]

A random unit vector uniformly distributed in sphere can be generated as follows:

  G4double z   = 2. * G4UniformRand() - 1.;
  G4double rho = std::sqrt((1. + z) * (1. - z));
  G4double phi = CLHEP::twopi * G4UniformRand();
  G4ThreeVector vector(rho * std::cos(phi), rho * std::sin(phi), z);

Please also see G4RandomDirection.hh.

This header contains two functions, G4RandomDirection() and G4RandomDirection(G4double cosTheta). First function generates a random unit vector uniformly distributed in a sphere. Second function generates a random unit vector uniformly distributed in a cone.

[Qi_Zhang]

I just don’t know why my code is wrong!

[Qi_Zhang]

I have applyed the G4RandomDirection() in my code. However it’s still wrong!

[allison]

Points or vectors uniformly distributed on a sphere are not uniform in theta, but in cos(theta). Think of Earth. Real estate runs out towards the poles.

[evc]

Why do you think that the result is wrong? Do you expect the yellow spot be in the center of the plot? If so, then place the detectors in the origin.

[Qi_Zhang]

I think right result should be similar with (0,0,50)!

[evc]

When you place the detectors at Z=+50,-50, the detectors are parallel to the XY plane, and the images correspond the expectation.

When you place them at Y=+50,-50, they should be parallel to the XZ plane. If not, then this can result in wrong images.

[Qi_Zhang]

I placed the detectors at Y=+50,-50` I get the following

I think it’s wrong!

[evc]

Please check the orientation of the detectors.

[Qi_Zhang]

I placed the detectors at Y=+50,-50` I get the following

I think it’s wrong!

[Qi_Zhang]

The codes of the detector is following:

 double rr=5.*cm;
 G4double phi = 90. * deg;
 G4ThreeVector u = G4ThreeVector(1, 0, 0);
 G4ThreeVector v = G4ThreeVector(0, std::cos(phi), std::sin(phi));
 G4ThreeVector w = G4ThreeVector(0, -1.*std::sin(phi), std::cos(phi));
 G4RotationMatrix *rotm1 = new G4RotationMatrix(u, v, w);
 G4ThreeVector position1 = rr * w;

 G4Box* det_box                                                                                                                      \

   = new G4Box("det_box",5.*cm,5.*cm,0.001*mm);
 G4LogicalVolume* det_log= new G4LogicalVolume(det_box,
                                               Vaccum,"det_log",0,0,0);
ew G4PVPlacement(rotm1, position1,
                    det_log,  // logical volume                                                                                       
                    "det_phys",  // name                                                                                              
                    experimentalHall_log,  // mother volume                                                                           
                    false,  // no boolean operation                                                                                   
                     0);  // copy number   

My code is wrong? please help me!

[evc]

It’s difficult to understand from the code of the detector, how did you place the detector in different positions. When I suggest to check the orientation of the detector, I mean that the detector should be rotated before the move to required position.

I would recommend you to use the constructor for G4PVPlacement with G4Transform3D, instead of the constructor with G4RotationMatrix and G4ThreeVector. The object of the G4Transform3D class can be constructed as a composition of a rotation and then a translation.

I would also recommend you to run your simulation in graphics mode, it will be much easier to detect the source of the problem.

[evc]

For example:

// Position at Z = +50
new G4PVPlacement(G4TranslateZ3D(5.*cm)*G4RotateX3D(0.*deg), ...

// Position at Y = +50
new G4PVPlacement(G4TranslateY3D(5.*cm)*G4RotateX3D(90.*deg), ...

// Position at Z = -50
new G4PVPlacement(G4TranslateZ3D(-5.*cm)*G4RotateX3D(180.*deg), ...

// Position at Y = -50
new G4PVPlacement(G4TranslateY3D(-5.*cm)*G4RotateX3D(270.*deg), ...

In such approach I would expect to see four very similar images.

[Qi_Zhang]

@evc

G4TranslateY3D(-5.*cm)*G4RotateX3D(270.*deg)

I have used “new G4PVPlacement(G4TranslateY3D(-5.*cm)*G4RotateX3D(270.*deg), …” to replace my code.

G4Box* det_box
   = new G4Box("det_box",5.*cm,5.*cm,0.001*mm);
 G4LogicalVolume* det_log= new G4LogicalVolume(det_box,
                                               Vaccum,"det_log",0,0,0);

 new G4PVPlacement(G4TranslateY3D(-5.*cm)*G4RotateX3D(270.*deg),det_log,"det_phys", experimentalHall_log, false, 0);///

However, it’s still wrong!

[jrellin]

If this isn’t some oversight in geometry than I suspect your confusion is, as @allison noted, in confusing uniform in theta with cos(theta). The alternative to normalizing would be to evenly sample on a sphere which is not a trivial problem. You would almost certainly want to use libraries such as healpix.

[Qi_Zhang]

@evc

G4RandomDirection

Now my confusion is that after I used G4RandomDirection() in my codes, the result is till wrong! Why?

particleGun->SetParticleMomentumDirection(G4RandomDirection());

[evc]

I’ve implemented the setup similar to yours. In the attachment you can find
DetectorConstruction.cc (3.3 KB)
PrimaryGeneratorAction.cc (3.3 KB)
for the basic example B1, where the detector has been placed at Y=-50. On the image below you can see the trajectories that intersect the detector. As expected, there are two spots with higher density of intersections points.