How to make pencil beam radius

Particles, Track, Event, Run and Biasing
1

Hello experts!

How to make pencil beam radius?

How to correct my code ?
Help please!

void SimpleParticleSource::GeneratePrimaries(G4Event* anEvent)
{
G4double radius = 0.4*cm;
G4double phi = 360.*deg;

    G4double rRan = radius*(G4UniformRand());
    G4double phiRan = phi*(G4UniformRand());
    G4double z0 = -6.8*cm;

fParticleGun->SetParticlePosition(G4ThreeVector(rRancos(phiRan)-9.18cm,rRan*sin(phiRan),z0));
if(flag) fParticleGun->GeneratePrimaryVertex(anEvent);
}

Best regards,

2

the way your code generates coordinates inside the disk oversamples positions close to the center.

I would recommend to use general particle source gps instead, there you can configure the shape to be a circle of desired radius.

3

Do you mean using macro file ? I don’t understand

4

yes, using macro is what I had in mind. There are probably several examples shipped with geant4, and also this: Geant4 GPS Examples

regarding the statement of oversampling at the center: math - Generate a random point within a circle (uniformly) - Stack Overflow

5

I don’t use macro file.

6 
void SimpleParticleSource::GeneratePrimaries(G4Event* anEvent)
{
    G4double radius = 0.4*cm;
    G4double phi =  360.*deg;

    G4double u = G4UniformRand() + G4UniformRand();
    G4double rRan = radius*(( u>1) ? 2-u : u);

    G4double phiRan = phi*G4UniformRand();
    G4double z0 = -6.8*cm;

    fParticleGun->SetParticlePosition(G4ThreeVector(rRan*cos(phiRan)-9.18*cm,rRan*sin(phiRan),z0));
    fParticleGun->GeneratePrimaryVertex(anEvent);
}
1 Like
7

@weller Why do you set u = G4UniformRand()+G4UniformRand() rather than simply u = 2*G4UniformRand()? I don’t think it fundamentally matters for the distribution, but I’m curious.

8

honestly — no clue. found the code snipped in there and switched off brain:

edit: they do a vivid discussion about the differences though!

9

Thanks, @weller ! From the StackOverflow comments:

Try using both schemes to generate random numbers and see what you get. 2*random() gives numbers uniformly distributed in the range 0 to 2. random()+random() gives you numbers in the range 0 to 2 but there will (usually) be more numbers near 1.0 than near 0.0 or 2.0. It’s like how rolling two dice and summing is more likely to give 7 than any other number.

I need to think about this. We do something entirely different, generating differential area elements:

  G4double rho = Rmax * sqrt(G4UniformRand());
  G4double phi = twopi * G4UniformRand();
  G4TwoVector point; point.setPolar(rho, phi);
10

the two approaches seem equally valid in result, yours is also mentioned in that StackOverflow discussion: Disk Point Picking -- from Wolfram MathWorld

maybe the performance (Sqrt vs. 2nd random throw) is the difference then?
edit: depends on the generator, for the ones available in matlab sqrt is faster.

11

Ah, you may be right. Naively, I would have expected a random generator to be more expensive than sqrt(), but possibly not.

12

to finally come back to the topic: if you want to avoid introducing macro files to your simulation, you can also issue commands in source code like this:

  G4UImanager* UImanager = G4UImanager::GetUIpointer();
  UImanager->ApplyCommand("/gps/pos/type Plane");
  UImanager->ApplyCommand("/gps/pos/shape Circle");
  UImanager->ApplyCommand("/gps/direction 0 0 1");
  UImanager->ApplyCommand("/gps/pos/centre -9.18 0. -6.8. cm");
  UImanager->ApplyCommand("/gps/pos/radius 0.4 cm");
1 Like
13

Thank You so much! It works

14

Hello!

How to rotate the circular plane source 30 degrees to axes X and 60 degrees to axes Z ?
(I mean the particles must interact with material by the form of circle not a elliptic which we see on the picture).

    G4double radius = 1.5*cm;
    G4double phi = 360.*deg;
    
    G4double u = G4UniformRand() + G4UniformRand();
    G4double rRan = radius*(( u>1) ? 2-u : u);
    
    G4double phiRan = phi*(G4UniformRand());
    G4double z0 = -6.8*cm;
    G4double x0 = -9.18*cm;

fParticleGun->SetParticlePosition(G4ThreeVector(x0,rRansin(phiRan),rRancos(phiRan)+z0));
fParticleGun->SetParticleMomentumDirection(G4ParticleMomentum(cos(30.0deg),cos(90.0deg),cos(60.0*deg)));

Screenshot from 2023-01-27 15-40-39
Screenshot from 2023-01-27 15-40-39750×439 43.2 KB

15

one way is to calculate the rotated base vectors and issue the commands /gps/pos/rot1 and /gps/pos/rot2, see
https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/GettingStarted/generalParticleSource.html#source-direction-and-angular-distribution

1 Like
16

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.