Monte Carlo Detector Simulation

1. Monte Carlo Detector Simulation Pat Ward Memorable quote from UA5:- A: “Monte Carlo isn’t the only source of truth” B: “It is in this experiment” C.P. Ward

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 C.P. Ward

3. A Simulated ATLAS event C.P. Ward

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. C.P. Ward

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. C.P. Ward

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 C.P. Ward

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/ C.P. Ward

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 . C.P. Ward

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. C.P. Ward

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. C.P. Ward

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. C.P. Ward

12. Simple Generators • Simple generators use a single seed. • With single seed, 32-bit computer -> maximum period ~230 (~109). • E.g. multiplicative linear congruential generator: si+1 = (a si + c) mod m a is well-chosen multiplier; constant c may be 0; m (slightly <) largest integer Integer seed siconverted to floating by dividing by m. • Not adequate for typical MC applications. C.P. Ward

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 C.P. Ward

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 C.P. Ward

15. Random Numbers from Common Distributions • Isotropic distribution in 3D • Probability density ~ dΩ ~ d(cosθ)dφ • cosθ uniform (2u1 – 1) • Φ uniform 2πu2 • Gaussian • z1 = sin2πu1√(-2lnu2) • z2 = cos2πu1√(-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] C.P. Ward

16. 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 Inverse Transform Method C.P. Ward

17. Enclose f(x) in shape C*h(x) where h(x) is easily generated function f, h normalized, C > 1 Generate random number x0 from h(x) Generate u = uniform random number in [0,1] Accept x0 if u*C*h(x)≤f(x) Otherwise reject and generate new x0 Choose C*h(x) as close to f(x) as possible for maximum efficiency: importance sampling Acceptance-Rejection Method C.P. Ward

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 > 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 . C.P. Ward

19. Another Aside…… • Reproducibility problem with Gaussian random numbers generated in pairs. Generator in CLHEP: • First entry – generate pair, return one of them and set flag • Second entry – return second number without calling uniform generator again and reset flag • Hence state not fully defined by uniform generator state – not reproducible C.P. Ward

20. 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 Bryan’s lectures. • Need standard, well-defined interface between generator and detector simulation program. • E.g. HEPEVT common block in Fortran or HepMC in C++ . C.P. Ward

21. HepMC • Particles and vertices in tree structure • Can be input to GEANT4 vertex Z Z vertex vertex e e q qbar C.P. Ward

22. 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. C.P. Ward

23. 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 (G4PreAssignedDecayProducts) C.P. Ward

24. 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. C.P. Ward

25. 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-section*mass to decide whether to interact – if not, pick another ν. • Pickmaterial, vertex location and generate kinematics. C.P. Ward

26. 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 certainmaterial • May be asensitive volume • May containother volumes • A physical volume represents the spatial positioning or placement of the logical volume within a mother (logical) volume. C.P. Ward

27. Example SCT_Barrel Tube, air SCT:Layer0 Tube, air SCT:Layer1 SCT:Layer2 SCT:Layer3 SCT:Interlink Tube, air Thermal shield etc. x2 (13 vols) Interlink segments, B6 bearings SCT:Layer0Active Tube, air SCT:Layer0Aux Tube, air SCT:Layer0Support Tube, air Bracket, harness SCT:Ski0 SCT:Clamp SCT:CoolingEnd x32 x2 x2 SCT:SupportCylinder SCT:Flange Modules, doglegs, cooling blocks x2 C.P. Ward

28. 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. C.P. Ward

29. Example Shape: Tube Section C.P. Ward

30. Example Shape: Polyhedron C.P. Ward

31. 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.. C.P. Ward

32. 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. C.P. Ward

33. 4) Tracking • Particles are transported through the detector geometry in a series of steps. • Particleat (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. C.P. Ward

34. 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. C.P. Ward

35. 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. C.P. Ward

36. 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 C.P. Ward

37. 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. C.P. Ward

38. 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. C.P. Ward

39. Interactions and Decays • Distance to interaction or decay characterised by mean free path λ(for each process). • Decay: λ=γvτ • v = velocity, τ = mean lifetime • Interaction: 1/λ = ρΣi{xiσi/mi} • ρ = density of material • xi = mass fraction • σi = cross-section • mi = mass • λ varies as particle loses energy or crosses boundary. } isotope i C.P. Ward

40. Interactions and Decays • Probability of surviving distance l is P(l) = exp(-nλ) wherenλ = ∫l0dl/λ(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. C.P. Ward

41. 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…….. C.P. Ward

42. Main Physics Processes • Photon: • Pair production, Compton scattering, photoelectric effect • All charged particles: • Ionization / δ-rays, multiple scattering • Electron / positron • Bremsstrahlung, annihilation (e+) • Hadron: • Hadronic interactions C.P. Ward

43. 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 δ-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. C.P. Ward

44. Range Cuts • Energy thresholds corresponding to range cuts in silicon C.P. Ward

45. 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. C.P. Ward

46. Total no of clusters in 1m events Vary StepLim C.P. Ward

47. Ionization Energy Loss • Important for all charged particles. • Above electron cut-off kinetic energy Tcut simulate discrete δ-ray production. • Below Tcut continuous energy loss along step. • Imposes limit on step length because cross-section energy-dependent. C.P. Ward

48. Continuous Energy Loss • For hadrons typically use Bethe-Bloch restricted energy loss: C.P. Ward

49. 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 > 2 MeV; at lower energy use Bragg model. • For each step, compute mean energy loss and then apply fluctuations . C.P. Ward

50. 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 C.P. Ward