Weihe Zeng
Tsinghua University
Dec. 4, 2017

Content

Random number generators and sampling

random number generator: Psudorandom, True random.
random sampling based on uniform RN for various distributions.
Random number utilities in GEANT4
Examples: position, direction, energy.

Primary Particles

Particles in GEANT4
Abstract Classes: G4VUserPrimaryGeneratorAction and G4VPrimaryGenerator
PariticleGun: set methods
GeneralParticleSource: UI commands

Pseudo-random Number Generator

A pseudorandom number generator (PRNG), also known as a deterministic random bit generator (DRBG), is an algorithm for generating a sequence of numbers whose properties approximate the properties of sequences of random numbers.

Select from possible values with equal probability (Uniform)
Each generation must be independent.
Long period for simulation.
~~Cryptographically secure against attacker (CSPRNG).~~

Linear Congruential Generator

$$ X{n+1} = (aXn + c) \text{ mod } m $$ $0 \leq X_0 \lt m$ is the start value or seed

Other generators implemented in GEANT4/CLHEP: James, MIXMAX, Ranlux, Ranecu.

True Random Number Generator

A hardware random number generator (TRNG) is a device that generates random numbers from a physical process, rather than a computer program.

Nuclear decay process
Thermal noise, atmospheric noise.

30%

Check out random device on your Operating System (Unix-like): /dev/random and /dev/urandom.

Comparison of PRNG and TRNG

| Characteristic | PRNG | TRNG | | -------------- | ------------ | ---------------- | | Efficiency | Excellent | Poor | | Determinism | Determinstic | Nondeterministic | | Periodicity | Periodic | Aperiodic |

| Application | Most Suitable Generator | | ------------------------ | ----------------------- | | Lotteries and Draws | TRNG | | Games and Gambling | TRNG | | Simulation and Modelling | PRNG | | Security | TRNG | | The Arts | Varies |

Statistical test of RNG

Visualization of random number in a 2 dimensional scatter plot.

140%

Ref: http://boallen.com/random-numbers.html

Random Sampling for a Probability Distribution

Problem Description:

Given a sequence of random numbers $u1, u2, u3,...$ uniformly distributed in [0,1], determine a sequence $x1,x2,x3,...$ distributed according to the p.d.f. $f(x)$.

Inverse transformation method

Generate a random number $u$ uniformly distributed in [0,1].
Compute the value $x$ such that $FX(x)=u$, i.e. $x = F^{-1}X(u)$.
Take $x$ to be the random number drawn from the distribution described by $F_X(x)$.

Example:

p.d.f. $f(x) = e^{-x},, x$ in $[0, \infty),, \text{c.d.f. } F(x) = 1-e^{-x}$

Inversion of c.d.f. is $F^{-1}(y) = -\ln(1-y)$

Acceptance-rejection method

Consider a p.d.f $f(x)$ with $x$ in [$x{min}, x{max}$] and a maximum value not greater than $f_{max}$.

Generate a random number x, uniformly distributed in [$x{min}, x{max}$] from $u1 \sim U[0,1]$, i.e. $x = x{min} + u1 (x{max} - x_{min})$.
Generate a new random number $y$ uniformly distributed in $[0,f{max}]$ from $u2$, i.e. $y = u2 f{max}$.
Accept $x$ if $y < f(x)$. If not, reject $x$ and repeat from step 1.

70%

Random Module in GEANT4/CLHEP

Textbook: Chapter 3.2.2, page 51-54.

HepRandom module in CLHEP (A Class Library for High Energy Physics) package.

Change random engines and set/get seed:

// default engine: HepJamesRandom
CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine); 

CLHEP::HepRandom::setTheSeed(seed); // set seed(G4long)  
CLHEP::HepRandom::getTheSeed();     // get current seed

Some built-in distributions:

RandFlat, RandGauss, RandExponential, RandBreitWigner,  
RandPoison, RandLandau, RandGamma ...

GEANT4 Wrapper for HepRandom

Defined in Randomize.hh, here are some simple usages.

#include "Randomize.hh"

// THE MOST IMPORTANT ONE
// uniform random number in [0,1]
G4double u = G4UniformRand();

// other distributions
G4double uGauss = G4RandGauss::shoot(mu,sigma);  
G4double uExpo = G4RandExponential::shoot(mean);  
G4double uFlat = G4RandFlat::shoot(a, b);

Examples of Random Sampling in Physics Simulation

NOTE: $u1, u2, u_3, ...$ are random numbers uniformly distributed in [0,1], e.g. from G4UniformRand().

Position:

Uniformly distributed inside a 3D box with $a, b, c$ as length, width and height.

$$ x = au1,, y = bu2, z = cu_3 $$

Uniformly distributed inside a circle plane centered at (0,0,0) with radius of $r$.

$$ x = r(2u1 - 1), y=r(2u2 -1), z=0$$

If $x^2 + y^2 > r^2$, repeat last step, else accept $x, y$.

Uniformly distributed on the lateral surface of a cylinder paralleled to z-axis and centered at (0,0,0) with radius $r$ and height $h$.

$$ \phi = 2\pi u1, z=hu2 - h/2,, x = r\cos(\phi), y=r\sin(\phi) $$

Direction (Momentum)

Isotropic distribution:

$$ \cos \theta = 2u1 - 1,, \phi = 2\pi u2 $$

$$px = \sin \theta \cos \phi,, py = \sin\theta \sin\phi,, p_z = \cos\theta$$

80%

$$ \mathrm{d} S = 2\pi a (r \mathrm{d} \theta) = 2\pi (r \sin \theta) ( r\mathrm{d} \theta) = -2\pi r^2 \mathrm{d} \cos\theta$$

Let $t=\cos\theta$, we can calculate surface of sphere cap as $S{cap} = \int{0}^{\theta0} \mathrm{d} S = \int{t=1}^{\cos\theta0} -2\pi r^2 \mathrm{d} t = 2\pi r^2 (1-\cos\theta0)$

Energy

Given discrete energy spectrum$(E(i), h(i))$, where $i=1,2,3,...,n$, $E(i)$ is the right boundary of i-th energy bin and $(i)h$ is the probability of energy in i-th bin.

60%

Calculate cumulated function $H(i)$ from $h(i)$.
Generate $u1$ and find $i$ for which $H(i-1) \leq u1 \lt H(i)$.
Sample $e$ uniformly within energy bin [$E(i-1), E(i)$].

Primary Particle Generators

Code references:

examples/basic/B1  
examples/extended/eventgenerator/  
source/event/[include|src]/G4SPS[Ang|Ene|Pos]Distritbution.[hh|cc]

Textbook:

Particles: Chapter 5.3, page 216-220  
G4ParticleGun: Chapter 2.6, page 15-17  
G4GeneralParticleSource: Chapter 2.7, page 17-27

Online source:

Geant4 Tutorial Course

Particles in GEANT4

There are a large number of elementary particles and nuclei. Geant4 provides the G4ParticleDefinition class to represent particles, and various particles, such as the electron, proton, and gamma have their own classes derived from G4ParticleDefinition. And particles in Geant4 have properties such as

PDG encoding, mass, width, electric charge
spin, isospin and parity, life time and decay modes

// code to get pre-defined particle
G4ParticleTable* particleTable =  
        G4ParticleTable::GetParticleTable();
G4ParticleDefinition* fGamma =  
        particleTable->FindParticle("gamma"),
fElec = particleTable->FindParticle("e-"),  
fNeut = particleTable->FindParticle("neutron");

Abstract Classes for Primary Generation

G4VUserPrimaryGeneratorAction is the abstract base class to control the primary particle generation.
G4VPrimaryGenerator is the abstract class that do the actual generation work.

Users must implement a concrete class of G4VUserPrimaryGeneratorAction and instantiate it in the Build() method in G4VUserActionInitialization.

// source/run/include/G4VUserPrimaryGeneratorAction.hh
class G4VUserPrimaryGeneratorAction{  
  public:
    G4VUserPrimaryGeneratorAction();
    virtual ~G4VUserPrimaryGeneratorAction();
  public:
    // Invoked at the beginning of each event
    virtual void GeneratePrimaries(G4Event* anEvent) = 0;
};

There are some built-in concrete implementations of G4VPrimaryGenerator.

We should use these classes to generate primary particles inside our implementation of G4VUserPrimaryGeneratorAction::GeneratePrimaries.

G4ParticleGun

I want “particle shotgun”, “particle machinegun”.

Instead of implementing such a fancy weapon, in the implementation of G4UserPrimaryGeneratorAction:

Shoot random numbers in arbitrary distribution
Use set methods of G4ParticleGun
Use G4ParticleGun as many times as you want

// set methods of G4ParticleGun
void SetParticleDefinition(G4ParticleDefinition*)  
void SetParticleMomentum(G4ParticleMomentum)  
void SetParticleMomentumDirection(G4ThreeVector)  
void SetParticleEnergy(G4double)  
void SetParticleTime(G4double)  
void SetParticlePosition(G4ThreeVector)  
void SetParticlePolarization(G4ThreeVector)  
void SetNumberOfParticles(G4int)

// generate method of G4ParticleGun
void GeneratePrimaryVertex(G4Event*)

Using G4ParticleGun in PrimaryGeneratorAction

// constructor of PrimaryGeneratorAction invoke only once
B1PrimaryGeneratorAction::B1PrimaryGeneratorAction()  
: G4VUserPrimaryGeneratorAction(),
  fParticleGun(0), 
  fEnvelopeBox(0)
{
  G4int n_particle = 1;
  fParticleGun  = new G4ParticleGun(n_particle);

  // default particle definition
  G4ParticleTable* particleTable = 
      G4ParticleTable::GetParticleTable();
  G4ParticleDefinition* particle = 
      particleTable->FindParticle("gamma");

  // set default particle properties  
  fParticleGun->SetParticleDefinition(particle);
  fParticleGun->SetParticleMomentumDirection(
    G4ThreeVector(0.,0.,1.));
  fParticleGun->SetParticleEnergy(6.*MeV);
}

Notes:

fParticleGun is declared as private member of B1PrimaryGeneratorAction
n_particle and SetNumberOfParticles method: generate n_particle each invoke of GeneratePrimaryVertex.
These default values can also be modified/randomized in GeneratePrimaries of PrimaryGeneratorAction.

void B1PrimaryGeneratorAction::GeneratePrimaries(  
                                    G4Event* anEvent)
{
  //this function is called at the begining of ecah event
  //  
  G4double envSizeXY = ...; // from logical volume
  G4double envSizeZ = ...;

  G4double size = 0.8; 
  G4double x0 = size * envSizeXY * (G4UniformRand()-0.5);
  G4double y0 = size * envSizeXY * (G4UniformRand()-0.5);
  G4double z0 = -0.5 * envSizeZ;

  fParticleGun->SetParticlePosition(
                             G4ThreeVector(x0,y0,z0));
  fParticleGun->GeneratePrimaryVertex(anEvent);

  /* you can repeat the steps and generate a second vertex 
  x0 = -x0;  
  fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
  fParticleGun->GeneratePrimaryVertex(anEvent);
  */
}

G4GeneralParticleSource (GPS)

compPG

Using G4ParticleGun in PrimaryGeneratorAction

examples/extended/eventgenerator/exgps

class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction{  
  public:
    PrimaryGeneratorAction();
   ~PrimaryGeneratorAction();
    virtual void GeneratePrimaries(G4Event*);
  private:
    G4GeneralParticleSource* fParticleGun;
};
PrimaryGeneratorAction::PrimaryGeneratorAction()  
 : G4VUserPrimaryGeneratorAction() {
   fParticleGun = new G4GeneralParticleSource();
}
PrimaryGeneratorAction::~PrimaryGeneratorAction(){  
  delete fParticleGun;
}
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent){  
  fParticleGun->GeneratePrimaryVertex(anEvent);
}

Controlling GPS with UI Commands

# spherical volume source, isotropic radiation, linear energy
/run/initialize

/gps/particle geantino
# source
/gps/pos/type Volume
/gps/pos/shape Sphere
/gps/pos/centre 2. 2. 2. cm
/gps/pos/radius 2.5 cm

# angular distribution
/gps/ang/type iso

# energy distribution
/gps/ene/type Lin
/gps/ene/min 2. MeV
/gps/ene/max 10. MeV
/gps/ene/gradient 1.
/gps/ene/intercept 0.

See more examples in examples/extended/eventgenerator/exgps/macros.

Visualization of Particles

Use /run/beamOn <number-of-events> after your usual geometry visualization (try exampleB1).

vis

Homework

Deadline: 18:00, Dec. 11, 2017

Read textbook chapter 2.6-2.7, page 15-27
Add primary particle source to the B1 example with your human model.
- B1PrimaryGeneratorAction should be modified. You can choose either G4ParticleGun or GPS to accomplish this work.
- Primary particle position uniformly distributed on the entire surface of a cylinder ($r=1m, h=2m$). If your model is a Giant, feel free to enlarge this cylinder.
- Initial directions of particles are isotropically distributed.
- Primary particles are $\gamma$ photon with 1 MeV energy.

Notes on examples/basic/B1 modification:

You should assign fScoringVolume private member of B1DetectorConstruction inside Construct() method to one of your logical volume, e.g.fScoringVolume = logicHead;
If you want to use GPS as primary generator, you need to comment out the following lines in src/B1RunAction.cc

const B1PrimaryGeneratorAction* generatorAction  
 = static_cast<const B1PrimaryGeneratorAction*>
   (G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction());
G4String runCondition;  
if (generatorAction)  
{
  const G4ParticleGun* particleGun = generatorAction->GetParticleGun();
  runCondition += particleGun->GetParticleDefinition()->GetParticleName();
  runCondition += " of ";
  G4double particleEnergy = particleGun->GetParticleEnergy();
  runCondition += G4BestUnit(particleEnergy,"Energy");
}
...
<< " The run consists of " << nofEvents << " "<< runCondition

Backup Slides

GPS macro

macro文件（*.mac）是geant4的输入文件，可以不修改源码，对模拟进行控制。使用方法是：program_name.exe macro_name.mac，main函数里这几行解释了为什么可以这样做：

    G4String command = "/control/execute ";
    G4String fileName = argv[1];
    UImanager->ApplyCommand(command+fileName);

GPS的用法在Official document 有详细的描述，但是还是例子更直观（？存疑）

例子见：examples/extended/eventgenerator/exgps

这里简单说明几种常见的源项设定，可以在exgps/macros里找到，说明见exgps/macros/README：

注1：这些指令都要在/run/initialize之后使用！ 注2：/analysis/ commands are independent of gps

位置抽样

球型体源

# test10.mac
/gps/pos/type Volume
/gps/pos/shape Sphere
/gps/pos/centre 2. 2. 2. cm
/gps/pos/radius 2.5 cm

方形面源

# test28.mac
/gps/pos/type Plane
/gps/pos/shape Square
/gps/pos/centre 1. 2. 1. cm
/gps/pos/halfx 2. cm
/gps/pos/halfy 2. cm

长方体源（可以不是直角的，parallelepiped，需要设置角度）

/gps/pos/type Volume
/gps/pos/shape Para
/gps/pos/centre 30. 0. 321.5 mm
/gps/pos/halfx 10. mm
/gps/pos/halfy 10. mm
/gps/pos/halfz 1. mm

额外条件，如果你的源项不规则，可以用上面的方法定义一个较大的源项，然后限定它：

/gps/pos/confine Electronics_Box_0 # name of the physical volume

角度抽样

默认是各向同性：/gps/ang/type iso
还有cos分布：/gps/ang/type cos
聚焦的源等

能量抽样

单能的源（默认）

/gps/energy 1.0 MeV

各种函数（linear, gaussian, power, exponential ...)
自定义的能量分布，中间插值

# test19.mac
# energy distribution
/gps/ene/type Arb
/gps/hist/type arb
/gps/hist/point 1.0 2.  # energy weight
/gps/hist/point 2.0 5.
/gps/hist/point 7.0 1.
/gps/hist/point 10. 1.
/gps/hist/inter Lin # linear interpolation

外部文件输入的能量分布

# test38.mac and spectrum.dat
# energy distribution
/gps/ene/type Arb
/gps/hist/file ./macros/spectrum.dat
/gps/hist/inter Exp

粒子类型

/gps/particle/particle_name

可以选的粒子很多：gamma，e-，mu+等，详见：geant4 particle

如果物理过程中定义了G4RadioactiveDecayPhysics，如果是不稳定的ion，会发生相应的衰变，例如：

# Uranium-238 source (assuming the nuclei is still)
/gps/particle ion
/gps/ion 92 238 0 0.
/gps/energy 0.0 MeV

完整的macro文件

/run/initialize

# source 
/gps/particle ion
/gps/ion 92 238 0 0.
/gps/energy 0.0 MeV

# in a box 
/gps/pos/type Volume
/gps/pos/shape Para
/gps/pos/centre 0. 0. 0. mm
/gps/pos/halfx 20. mm
/gps/pos/halfy 20. mm
/gps/pos/halfz 10. mm

/run/beamOn 100000

Menu

[Tutorial] GEANT4 Introduction: Random Sampling & Primary Particles

Content

Random number generators and sampling

Primary Particles

Pseudo-random Number Generator

True Random Number Generator

Comparison of PRNG and TRNG

Statistical test of RNG

Random Sampling for a Probability Distribution

Inverse transformation method

Acceptance-rejection method

Random Module in GEANT4/CLHEP

GEANT4 Wrapper for HepRandom

Examples of Random Sampling in Physics Simulation

Position:

Direction (Momentum)

Energy

Primary Particle Generators

Particles in GEANT4

Abstract Classes for Primary Generation

G4ParticleGun

Using G4ParticleGun in PrimaryGeneratorAction

G4GeneralParticleSource (GPS)

Using G4ParticleGun in PrimaryGeneratorAction

Controlling GPS with UI Commands

Visualization of Particles

Homework

Backup Slides

GPS macro

位置抽样

角度抽样

能量抽样

粒子类型

完整的macro文件