Monte Carlo Detector Simulation

1
Monte Carlo Detector Simulation

• Pat Ward

Memorable quote from UA5- A Monte Carlo isnt
the only source of truth B It is in this
experiment
2
Monte Carlo Detector Simulation

• Introduction why, how, random numbers, Monte
  Carlo techniques
• Event Generation
• Geometry and Materials
• Tracking
• Physics Processes
• Hits and Digitisation
• Fast Simulation

3
A Simulated ATLAS event
4
1) Introduction

• Why do we use MC simulations so extensively?
• Design detector
• Understand detector response
• Develop reconstruction and analysis algorithms
• Estimate efficiencies and backgrounds
• These studies need large numbers of simulated
  events often many times the number of data
  events.

5
Full / Fast Simulation

• Full simulation track particles through detailed
  detector geometry, taking account of interactions
  and decays. Secondary particles produced in
  interactions of primaries with detector material
  are tracked.
• Fast simulation parameterize response of
  detector to primary particles.
• Hybrid full tracking in parts of detector (e.g.
  tracking chambers) with parameterization in other
  regions, e.g. calorimeters.
• Will first concentrate on full simulation.

6
All steps in one job or output results at each
stage and input to next stage
Event Generator
Detector Geometry
Digitisation
Tracking
Hits
Primary Event
Digits (Raw data)
Physics Processes
Reconstruction
7
MC Programs

• Most experiments now use GEANT4 toolkit for
  detector simulation (OO, C).
• Contains tracking, physics processes, tools for
  user to build detector geometry, visualisation.
• Developed from earlier Fortran GEANT3, used by
  e.g. LEP experiments.
• Most of following applies to GEANT4.
• S.Agostinelli et al., Nucl. Inst. And Methods
  A506 (2003) 250
• J.Allison et al., IEEE Transactions 53 (2006)
  270.

http//geant4.web.cern.ch/geant4/
8
Random Number Generation

• Random number generation obviously important for
  any Monte Carlo simulations.
• Has caused numerous problems when the random
  number generator was not random enough e.g.
  periodic behaviour.
• A single event uses huge numbers of random
  numbers e.g. a typical atmospheric neutrino
  event (mean energy 1 GeV) in MINOS uses 2.3x105
  .
• Assuming scales with energy, LHC event needs
  3x109 .

9
Types of Random Number Generators

• Truly random from random physical process, e.g.
  radioactive decay.
• Pseudorandom generated from simple numerical
  algorithm, appear random if algorithm not known.
• Quasirandom from numerical algorithm, but
  designed to be as uniform as possible rather than
  to appear random. Used for MC integration.
• Monte Carlo simulations use pseudorandom numbers.

10
Random Number Generators

• Generators for any distribution depend on a basic
  generator producing numbers uniform in 0,1.
• N.B. endpoints of range may or may not be
  included if included beware log(0).
• Set state using initial seed or multiple seeds.
  
• State at any time defined by one or more seeds,
  which can be retrieved for restarting generator.
• All pseudorandom number generators have some
  period after which sequence repeats.

11
Required Properties

• Good distribution i.e. randomness.
• Long period no repeat events.
• Repeatability need to be able to regenerate a
  particular event in middle of run to debug
  problem, for example.
• Long disjoint sequences run many jobs in
  parallel.
• Portability use many different systems.
• Efficiency.

12
Simple Generators

• Simple generators use a single seed.
• With single seed, 32-bit computer -gt maximum
  period 230 (109).
• E.g. multiplicative linear congruential
  generator si1 (a si c) mod m
• a is well-chosen multiplier constant c may
  be 0
• m (slightly lt) largest integer
• Integer seed si converted to floating by
  dividing by m.
• Not adequate for typical MC applications.

13
Some Random Number Generators

• Various techniques exist for improving simple
  generators.
• See F. James, Comp. Phys. Comm. 60 (1990)
  329.
• Some of the commonly used generators, available
  in GEANT4, are
• Ranecu
• From Cernlib RANECU
• Period 1018
• Not easy to make disjoint sequences
• Initialized / restarted by one 32-bit integer

14
Some Random Number Generators

• Ranlux
• From Cernlib RANLUX improvement on RCARRY
• 109 disjoint sequences of length 10161
• Initialized by one 32-bit integer, but need 25
  words to store state
• HepJamesRandom
• From Cernlib RANMAR
• Default in GEANT4
• 109 disjoint sequences of length 1034
• Initialized by one 32-bit integer, but need 100
  words to store state

15
Random Numbers from Common Distributions

• Isotropic distribution in 3D
• Probability density dO d(cos?)df
• cos? uniform (2u1 1)
• F uniform 2pu2
• Gaussian
• z1 sin2pu1v(-2lnu2)
• z2 cos2pu1v(-2lnu2)
• z1 and z2 are independent and Gaussian
  distributed with mean zero and standard deviation
  1
• N.B. there are faster variants

u1, u2 are random numbers uniform in 0,1
16
Inverse Transform Method

• Want to generate random number x from f(x)
• f(x) normalized
• Cumulative distribution F(a) ?af(x)dx
• If a is chosen randomly from f(x), then F(a) is
  uniform in 0,1
• Choose u uniform in 0,1
• Find x F-1(u)
• Good for simple, easily inverted functions

• Useful for obtaining random number according to
  histogram

17
Acceptance-Rejection Method

• Enclose f(x) in shape Ch(x) where h(x) is easily
  generated function
• f, h normalized, C gt 1
• Generate random number x0 from h(x)
• Generate u uniform random number in 0,1
• Accept x0 if uCh(x)f(x)
• Otherwise reject and generate new x0

• Choose Ch(x) as close to f(x) as possible for
  maximum efficiency importance sampling

18
Aside.

• Even with a good basic uniform random number
  generator beware problems in distributions.
• E.g GEANT3 generated Poisson distribution using
  GPOISS which returned integer.
• For mean gt 16, routine used Gaussian approx.
  which calculated a float.
• Result inadvertently rounded down when GPOISS
  returned integer.
• So, e.g. many calls with mean 16.1 would give
  distribution with mean 15.6 .

19
2) Event Generation

• At colliders, generally use stand-alone program
  to generate final-state particles of primary
  event for input to detector simulation.
• Examples are Pythia, Herwig discussed in
  Bryans lectures.
• Need standard, well-defined interface between
  generator and detector simulation program.
• E.g. HEPEVT common block in Fortran or HepMC in
  C .

20
HepMC

• Particles and vertices in tree structure
• Can be input to GEANT4

vertex
Z
Z
vertex
vertex
e
e
q
qbar
21
Particle Codes

• Particles are identified by a unique code often
  called the PDG code.
• Particles positive code, antiparticles negative.
• Quarks d 1, u 2, s 3,
• Leptons e- 11, ?e 12, µ- 13, ?µ 14, .
• Mesons and baryons spin and quark content are
  incorporated in code.
• See Monte Carlo Particle Numbering Scheme in
  PDG.

22
Short-lived Particles

• Event generators usually decay resonances, and
  may optionally decay other short-lived particles.
• Specialized packages may be used e.g. TAUOLA
  for tau decays, EvtGen for B decays.
• At detector simulation stage, these particles
  need to be tracked until decay, but then use the
  decay products already generated.
• This is possible in GEANT4 (G4PreAssignedDecayProd
  ucts)

23
Fixed-Target Event Generation

• For fixed-target experiments the primary event
  generation may be hard to separate from detector
  simulation.
• Simulate primary interaction within detector
  simulation program often with specialized
  generator, e.g. NEUGEN for neutrino interactions.
• Tracking beam particles through detector until
  interact not feasible for neutrinos -
  cross-sections far too small.
• Neutrino beam contains range of energies, and
  detector unlikely to be homogenous. Cross-section
  depends on energy and target nucleus.
• Calculate maximum interaction probability at
  initialisation and reweight all probabilities to
  get efficient generation.

24
Neutrino Event Generation

• Initialize detector geometry, neutrino flux
  generator and neutrino interaction generator.
• Estimate maximum weight (i.e. interaction
  probability).
• In event loop
• Pick neutrino from flux generator.
• Calculate accumulated mass distribution,
  cross-sections, for materials along ? trajectory.
  
• Use total cross-sectionmass to decide whether to
  interact if not, pick another ?.
• Pick material, vertex location and generate
  kinematics.

25
3) Geometry and Materials

• Detector geometry is modelled as a hierarchy of
  volumes (i.e. tree structure).
• A logical volume represents a detector element of
  certain shape
• Made of certain material
• May be a sensitive volume
• May contain other volumes
• A physical volume represents the spatial
  positioning or placement of the logical volume
  within a mother (logical) volume.

26
Example
SCT_Barrel Tube, air
SCTLayer0 Tube, air
SCTInterlink Tube, air
SCTLayer1 SCTLayer2 SCTLayer3
Thermal shield etc.
x2
(13 vols)
Interlink segments, B6 bearings
SCTLayer0Support Tube, air
SCTLayer0Active Tube, air
SCTLayer0Aux Tube, air
Bracket, harness
SCTSki0
SCTClamp
SCTCoolingEnd
x32
x2
x2
SCTFlange
SCTSupportCylinder
Modules, doglegs, cooling blocks
x2
27
Geometry

• A logical volume may be
• A simple shape e.g. box, cylinder, cone.
• Made from boolean combinations of these shapes
• In Geant4 it is also possible to define shapes
  using bounding surfaces (useful for interfacing
  to CAD systems).
• Physical volumes are positioned within the mother
  logical volume using translation and rotation
  relative to the local coordinate system.
• The daughter volumes must not overlap, and must
  not protrude from the mother volume.

28
Example Shape Tube Section
29
Example Shape Polyhedron
30
Materials

• An element can be defined by atomic number and
  atomic mass.
• A material is built from elements or other
  materials need to specify density and fraction
  (by mass) of each element / material.
• From composition and density useful properties
  (radiation length, interaction length, dE/dx etc)
  can be calculated.
• Specifying correct composition is important for
  simulation of energy loss, multiple scattering,
  photon conversion probability etc..

31
Sensitive Volume

• A logical volume may be a sensitive volume one
  for which hits are stored during tracking.
• In GEANT4 a sensitive detector may have a readout
  geometry.
• The readout geometry is completely separate from
  the tracking geometry it is used to determine
  the readout channel for a hit, for example.

32
4) Tracking

• Particles are transported through the detector
  geometry in a series of steps.
• Particle at (x,y,z) with momentum (px,py,pz).
• Calculate step to nearest boundary.
• Exit current volume or enter daughter.
• For each physics process calculate distance to
  interaction or decay.
• Take shortest step.
• Apply dE/dx, multiple scattering and any other
  continuous processes.
• Interact / decay particle if that was what
  limited step.

33
Tracking Optimization

• Tracking of particles is the time-consuming part
  of simulation.
• Various techniques to reduce number of candidate
  volumes for intersection.
• Virtual Divisions (used in Geant3).
• Slice mother volume evenly along one axis.
• Each section stores list of daughters contained.
• Works well in hierarchical geometry.

34
Tracking Optimization

• Smart Voxels used in Geant4.
• One-dimensional virtual division for each mother
  volume.
• Combine divisions containing same volumes.
• If division has too many volumes divide again
  along a different axis.
• Still too many volumes ? divide again along third
  axis.

35
Smart Voxels Example

• Mother volume sliced along horizontal axis.
• Each slice has independent set of vertical
  divisions.

From S. Agostinelli et al., Nucl. Inst. Meth.
A506 (2003) 250
36
Tracking in EM Field

• Charged particle does not follow linear
  trajectory in magnetic field.
• Propagate by integrating equation of motion.
• Uniform field ? helical path, analytic solution.
• Non-uniform field ? non-analytic solution.
• Geant4 uses 4th order Runge-Kutta method by
  default.
• Best method of solution depends on field.
• Near-uniform field start from helical solution
  and use Runge-Kutta iteration.
• Field not smooth use lower order method.

37
Tracking in EM Field

• Curved path broken into linear chord segments.
• Chord segments used to test intersection with
  boundary.
• Need to choose parameters (e.g. miss distance)
  carefully to give desired accuracy, else momentum
  bias.

38
Interactions and Decays

• Distance to interaction or decay characterised by
  mean free path ? (for each process).
• Decay ??vt
• v velocity, t mean lifetime
• Interaction 1/? ?Sixisi/mi
• ? density of material
• xi mass fraction
• si cross-section
• mi mass
• ? varies as particle loses energy or crosses
  boundary.

isotope i
39
Interactions and Decays

• Probability of surviving distance l is
• P(l) exp(-n? ) where n? ?l0 dl/?(l).
• P(n?) is exponential, independent of material and
  energy.
• At particle production point set n? - ln(?)
  where ? is random number uniform in 0,1.
• Different random number for each process.
• Use this to determine point of interaction or
  decay in current material.
• Update remaining n? after each step if particle
  still alive.

40
5) Physics Processes

• Physics processes may happen
• Continuously along a step (e.g. energy loss).
• At the end of a step (decay, interaction).
• At rest (e.g. decay at rest).
• Most important processes (for our simulations)
  are electromagnetic and hadronic interactions.
• N.B. Geant4 has wide range of physics processes
  available, covering wide energy range e.g.
  Compton scattering, photoelectron production,
  synchrotron radiation..

41
Main Physics Processes

• Photon
• Pair production, Compton scattering,
  photoelectric effect
• All charged particles
• Ionization / d-rays, multiple scattering
• Electron / positron
• Bremsstrahlung, annihilation (e)
• Hadron
• Hadronic interactions

42
Range Cuts

• In Geant3 particles below (user-defined) energy
  were stopped and remaining energy deposited ?
  results depend on cut-off used.
• In Geant4 all particles are tracked to end of
  range.
• But number of secondaries increases rapidly as
  energy decreases for d-rays, bremsstrahlung ?
  necessary to have cut-off for production of
  secondary particles.
• Cuts (for e,e-,?) given as range internally
  translated into energy for each material.

43
Range Cuts

• Energy thresholds corresponding to range cuts in
  silicon

Range e e- photon
5 µm 6.7 keV 6.7 keV (990 eV)
10 µm 31 keV 31 keV (990 eV)
50 µm 79 keV 80 keV 1.7 keV
100 µm 119 keV 120 keV 2.3 keV
200 µm 181 keV 185 keV 3.2 keV
1 mm 521 keV 541 keV 6.9 keV
44
Range Cuts

• Different range cuts are appropriate for
  different parts of detector.
• Sensitive detector need range cuts comparable
  to spatial resolution.
• Inert material not interested in detailed
  tracking of particles which do not emerge, use
  higher range cuts to improve speed.
• Example ATLAS SCT (silicon wafers 285 µm thick,
  strip pitch 80 µm). Number of clusters for single
  muon depends on range cut for electrons.

45
Total no of clusters in 1m events
Vary StepLim
46
Ionization Energy Loss

• Important for all charged particles.
• Above electron cut-off kinetic energy Tcut
  simulate discrete d-ray production.
• Below Tcut continuous energy loss along step.
• Imposes limit on step length because
  cross-section energy-dependent.

47
Continuous Energy Loss

• For hadrons typically use Bethe-Bloch restricted
  energy loss

48
Continuous Energy Loss

• Shell correction accounts for interaction of
  atomic electrons with nucleus 10 correction
  for T 2 MeV protons.
• Density effect arises from polarization of
  medium important at high energies.
• Bethe-Bloch good for T gt 2 MeV at lower energy
  use Bragg model.
• For each step, compute mean energy loss and then
  apply fluctuations .

49
Multiple Scattering

• Charged particles traversing medium deflected by
  many small-angle scatters Coulomb scattering
  from nuclei.
• Detailed model simulate all interactions (thin
  foils).
• Condensed model simulate global effects at end
  of track segment (usually used in MC, e.g.
  Geant4).
• Path length correction
• Net displacement
• Change of direction
• For small deflection angles distribution is
  approx. Gaussian with r.m.s. angular deviation in
  a plane

50
Multiple Scattering

• For larger deflections ( gt few ?0) behaves like
  Rutherford scattering longer tails.
• Geant4 uses theory of Lewis (does not assume
  small angle) H. W. Lewis Phys. Rev. 78 (1950)
  526.
• All theoretical distributions have some error
  simulation results depend on step length.
• Older versions of Geant (to 4.7) MS not allowed
  to limit step.
• Newer releases (4.8) MS can limit step (for low
  energy particles). Reduces speed by factor 2
  for ATLAS.

51
Hadronic Interactions

• Wide range of incident particle energies, from
  thermal neutrons to TeV.
• Calorimeters at colliders modelling of pion -
  nucleon interactions fundamental (cross-sections,
  pT distributions).
• Good model of muon production in hadronic
  interactions needed for accurate background in
  muon detectors.
• Different models appropriate for different
  incident particles, energy ranges, applications.

52
Hadronic Interactions

• Model types
• Data driven low energy neutron transport
• Parameterised hadronic showers in calorimeters
  (e.g. GHEISHA)
• Theory based parton string models used for
  extrapolation to high energy
• In Geant4 usually specify the physics list
  set of models to be used for various processes,
  e.g. LHEP, QGSP, QGSP_BERT, QGSP_BIC

53
Physics Lists in Geant4

• LHEP
• Parameterised models (LEP and HEP for low and
  high energies respectively) based on GHEISHA.
• Fast.
• Good shower shape, not so good for detailed
  interactions.
• Used as backbone for other lists (i.e. to fill in
  gap in energy range or particle coverage).

54
Physics Lists in Geant4

• QGSP
• Quark gluon string model for high energy p, n, p,
  K (gt 20 GeV) (uses LHEP for lower energy).
• Creates excited nucleus passed to precompound
  model for nuclear de-excitation.
• Theoretically motivated.
• Problems with shower shape at high energy.
• QGSC
• As QGSP but CHIPS (Chiral Invariant Phase Space)
  model for nuclear de-excitation.

55
Physics Lists in Geant4

• Above may be combined with intra-nuclear
  transport models at low energy.
• QGSP_BERT
• As QGSP but uses Bertini cascade model for p, n,
  p, K lt 10 GeV.
• More secondary protons and neutrons.
• Slower.
• QGSP_BIC
• As QGSP but uses binary cascade model for p, n,
  p, K lt 3 GeV.

56
ATLAS Tile Calorimeter Test Beam
Energy deposited in sensitive elements for
incident pions
Average energy greater with Bertini cascade
Shower longer with Bertini cascade Better
agreement with data
57
6) Hits and Digitisation

• As particles traverse a sensitive volume, record
  information about position, energy loss etc.
• These records are termed hits.
• The process which simulates the detector
  response, i.e. turns the hits into the detector
  signals, is termed digitisation.
• The output of the digitisation process is termed
  digits (or maybe raw data).

58
Hits

• The information stored in a hit depends on the
  type of detector.
• For a tracking detector, typically store step
  start and end points, step length, energy
  deposited, time, particle.
• For a calorimeter such detailed information is
  usually unnecessary, and consumes too much space.
• Store sum of energy deposited in each readout
  cell rather than individual steps.

59
Digitisation

• Details of digitisation very detector dependent.
• Example silicon strip detector (ATLAS SCT)
• Convert energy deposited to charge.
• Electron-hole pair creation energy 3.63
  eV/pair.
• Drift charges to surface of wafer, including
  effects of diffusion (depends on temperature),
  Lorentz angle (depends on bias voltage).
• Sum charges on each strip.
• Add noise (using measured noise occupancy).
• Simulate electronics response.

60
Digitisation

• Example MINOS scintillator strips
• Convert energy deposited to photons.
• Transport photons along wls fibres allowing for
  attenuation, reflection (model efficiency v
  distance).
• Convert photons to electrons at photocathode
  (quantum efficiency).
• Model readout electronics.
• Include noise and cross-talk.

61
7) Fast Simulation

• Full detector simulation is very cpu (and memory)
  intensive for high energy events.
• 10 15 mins / event (or more) for ATLAS
• Electromagnetic shower number of particles
  (hence time) is proportional to energy.
• Usually not feasible to generate large samples of
  events with different generator parameters needed
  for systematic studies, for example.
• Use fast generation techniques instead may be
  adequate for parameter comparisons if not for
  detailed studies of detector performance.

62
Fast Simulation

• Fastest simulation simply smear parameters of
  particles from event generator. 1s / event for
  ATLAS
• May include detailed parameterization of detector
  response as function of production angle (e.g.
  reflecting number of layers of tracking detector
  traversed, thickness and resolution of
  calorimeter etc).
• Hard to include effects of overlapping particles
  shared tracking hits, spatial resolution in
  calorimeters.
• Hard to simulate tails of distributions.
• Hybrid techniques combining full simulation of
  tracking detector with parameterization of
  calorimeter increase speed while still giving
  accurate simulation.

63
EM Shower Simulation

• Most time is taken by detailed simulation of
  showers in calorimeters.
• Methods to reduce this
• Dump energy of low energy particles (lt 10 MeV)
  instead of tracking to end of range.
• Parameterization of showers.
• Frozen shower libraries.
• N.B. need parameterization for electrons and
  positrons only photons soon convert.
• Hadronic showers contain large EM component, so
  hope gains in EM showers help these.

64
EM Shower Parameterization

• EM showers have well-defined (average)
  longitudinal and transverse profile which depend
  on radiation length.
• Longitudinal profile
• Radial profile depends on Molière radius,
  RMX0Es/Ec
• (Ec critical energy, Es 21 MeV)
• Can be parameterized as sum of 2 Gaussians or
  by

65
EM Shower Parameterization

• Determine parameters from full simulation.
• Validate full simulation using test beam
  data.
• During event simulation, when electron enters, or
  is produced in, calorimeter
• Stop tracking.
• Use parameterization to determine mean energy
  deposit in each calorimeter cell.
• Apply fluctuations.
• Parameterization valid above a few GeV.
• Also slow at low energy in fact full
  tracking may be faster.
• At low energy use frozen shower (or continue full
  tracking).

66
Frozen Shower Library

• Generate special samples of electrons incident on
  calorimeter and store hits.
• During routine event generation, retrieve one of
  these pre-generated showers instead of tracking
  particle (matched in energy, angle).
• Map to current shower position.
• Useful for low energies (lt 1 GeV) where
  parameterization not accurate.
• Optimal scheme may combine all techniques.

67
Example Scheme
Gives factor 4 reduction in simulation time in
ATLAS calorimeter for dijet events.
68
Summary

• MC simulations are crucial to design and
  understand detector and extract physics from
  data.
• Full simulation consists of several steps event
  generation tracking through detector geometry
  taking account of energy loss, multiple
  scattering, decays and interactions simulation
  of detector response.
• Most full detector simulations are carried out in
  the GEANT4 framework. http//geant4.web.cern.ch/g
  eant4/
• Full simulation is slow fast or hybrid
  techniques also needed.
• Understanding of these steps is necessary part of
  physics analysis.

69
And Finally

• Monte Carlo simulations are only useful if they
  are validated against (your) data!