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.

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.

## 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

1. Generate a random number $u$ uniformly distributed in [0,1].
2. Compute the value $x$ such that $FX(x)=u$, i.e. $x = F^{-1}X(u)$.
3. 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}$.

1. 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})$.
2. Generate a new random number $y$ uniformly distributed in $[0,f{max}]$ from $u2$, i.e. $y = u2 f{max}$.
3. Accept $x$ if $y < f(x)$. If not, reject $x$ and repeat from step 1.

## 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:

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

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

1. 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$.

1. 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$$

$$\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.

1. Calculate cumulated function $H(i)$ from $h(i)$.
2. Generate $u1$ and find $i$ for which $H(i-1) \leq u1 \lt H(i)$.
3. 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:

## 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)

### 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

# angular distribution
/gps/ang/type iso

# energy distribution
/gps/ene/type Lin
/gps/ene/min 2. MeV
/gps/ene/max 10. MeV
/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).

# Homework

1. Read textbook chapter 2.6-2.7, page 15-27
2. 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 有详细的描述，但是还是例子更直观（？存疑）

#### 位置抽样

1. 球型体源
# test10.mac
/gps/pos/type Volume
/gps/pos/shape Sphere
/gps/pos/centre 2. 2. 2. cm

1. 方形面源
# 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


1. 长方体源 （可以不是直角的，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

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

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

#### 角度抽样

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

#### 能量抽样

1. 单能的源（默认）
/gps/energy 1.0 MeV

1. 各种函数（linear, gaussian, power, exponential ...)

2. 自定义的能量分布，中间插值

# 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


1. 外部文件输入的能量分布
# 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

# 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