GAMOS an easy and flexible framework for GEANT4 simulations Pedro Arce (CIEMAT, Madrid)

1. GAMOS an easy and flexible framework for GEANT4 simulations Pedro Arce (CIEMAT, Madrid) CIEMAT, 5th July 2012

2. Outline • GAMOS plug-in’s • Applications • PET/SPECT • Radiotherapy • Accelerator geometry • Phase space files • Phantom voxelised geometries • Dose in phantom • Analysis and visualization • Optimisationn • Installation/Documentation/User support • Installation • Documentation • Release and user support • Summary • Introduction • MC applications • GAMOS objectives • Defining the input • Geometry from text files • Movements • Primary generator • Physics • Extracting detailed data • GAMOS data • Filtering • Classifying • Scoring • Sensitive detectors • Histogramming • Parameter management • Verbosity • Optimisation

3. Introduction

4. MC applications • Often the use a MC simulation is a difficult task for a non-expert software user • In the case of Geant4 most of the application has to be written in C++ • Several applications try to facilitate the use of MC in a specific field • Providing an scripting language tailored to the field

5. MC applications But users find some problems with these applications • User wants to describe some input not included in the application • A peculiar volume shape, a new primary generator position distribution, some option in the physics,… • User wants to have some kind of output or detailed information to debug or understand some part of the simulation, and the application only provides a limited amount of output possibilities • Dose from the gammas that entered the phantom with small energy, energy lost by particle traversing a volume as a function of the initial energy, … • The available applications cover a limited amount of physics fields, and many users do not find an appropriate one for their needs

6. GAMOS objectives • After almost ten years providing different utilities to make easier the use of Geant4 we created GAMOS with the objectives: • to provide a framework not only easy-to-use but also flexible, that may cover the needs of a wide range of Geant4 users • (a framework should not become an obstacle for an advanced user who wants to exploit all the functionality that Geant4 offers) • Comprehensive scripting language • Very wide range of input possibilities • Allow extraction of detailed information at any moment of the simulation Easily extendible with new user code • Provide a mechanism so that users can easily extend the framework functionalities with their code • Without the need to understand how GAMOS works internally

7. GAMOS history • GAMOS became a framework during first months of 2007 • to help in the simulation of CIEMAT PET detector projects • it was soon extended for radiotherapy applications GAMOS is the fruit of +10 years experience developing Geant4-based frameworks • Member of Geant4 since 1998 • Coordination of the CMS Geant4 framework during its first years • The first LHC Geant4 framework to be fully operative • Coordination the HARP (PS) Geant4 framework (by nomination from S. Giani, first Geant4 spokesman) CIEMAT group had several years experience with PET detectors and radiotherapy

8. Defining the input

9. Geant4 simulation • To do simulation with Geant4 a user has first to define the input: • Geometry • Optionally if a volume is a sensitive detector • Primary generator (initial particles) • Physics • There is no predefined physics in Geant4, user has to choose the one more appropriate for her/his simulation And the output: • Histograms • Scores (counts) • Files with data for later analysis The output is mainly obtained through: • User actions • Sensitive detectors / hits • Scorers

10. Geometry • Three different ways to define it: • C++ code: • The usual GEANT4 way • Add one line to transform your class in a plug-in DEFINE_GAMOS_GEOMETRY (MyGeometry); so that you can select it in your input macro /gamos/geometry MyGeometry Using one of the GAMOS examples for simple cases • Simple PET can be defined through an 8-parameters file (n_crystals, crystal_x/y/z, radius, …) • Simple phantoms with a few user commands Define it in ASCII (text) files • The easiest way to define a geometry • Based on simple tags • Same order of parameters as corresponding GEANT4 classes

11. Geometry from text files MATE Cu 29 63.54 8.9333*g/cm3 :VOLU yoke TUBE  62.*cm 820.*mm 1.27*m Cu :ROTM RM0 0. 0. 0. :PLACE  yoke  1   expHall  RM0   0.   0.   370.*cm • Based on simple tags, with same order of parameters as corresponding GEANT4 classes :MATE Cu 29 63.54 8.9333*g/cm3 :VOLU yoke TUBE  62.*cm 820.*mm 1.27*m Cu :ROTM RM0 0. 0. 0. :PLACE  yoke  1   expHall  RM0   0.   0.   370.*cm MATERIALS: • Isotopes • Elements • Simple materials • Material mixtures by weight, volume or number of atoms • GEANT4 intrinsic materials • + 300 predefined materials common in the medical field SOLIDS: • All GEANT4 “CSG” and “specific” solids • Twisted solids • Tessellated solids • Boolean solids

12. Geometry from text files • ROTATIONMATRICES: • 3 rotation angles around X,Y,Z • 6 theta and phi angles of X,Y,Z axis • 9 matrix values PLACEMENTS: • Simple placements • Divisions • Replicas • Assemblies • Parameterisations • Linear, circular, square or phantom-like • For complicated parameterisations example of how to mix the C++ parameterisation with the ASCII geometry file COLOUR VISUALISATION ON/OFF

13. Geometry from text files • PARAMETERS: • Can be defined to use them later :P InnerR 12. :VOLU yoke :TUBS  Iron   3  $InnerR 400. 300. :VOLU yoke2 :TUBS  Iron   3  $InnerR 200. 200. ARITHMETIC EXPRESSIONS: • Elaborated expressions can be used :SOLID yoke TUBE sin($ANGX)*2+4*exp(1.5)*cm 820. 1270.*mm UNITS: • Default units for each parameter • Each value can be overridden by user INCLUDE OTHER FILES (hierarchical approach): #include mygeom2.txt. :P InnerR 12. :VOLU yoke :TUBS  Iron   3  $InnerR 400. 300. :VOLU yoke2 :TUBS  Iron   3  $InnerR 200. 200.

14. Geometry from text files • User can extend it: add new tags and process it without touching base code • Can mix a C++ geometry with a text geometry • GEANT4 in memory geometry  text files • Install and use it as another GEANT4 library G4VPhysicalVolume* MyDetectorConstruction::Construct(){ G4tgbVolumeMgr* volmgr = G4tgbVolumeMgr::GetInstance(); volmgr->AddTextFile(filename); // Several files can be added return = volmgr->ReadAndConstructDetector(); HISTORY: • In use to build GEANT4 geometries since 10 years ago • An evolving code… • It is part of the GEANT4 release since December ‘08 :ROTM R00 0. 0. 0.:VOLU world BOX 100. 100. 100. G4_AIR:VIS world OFF:VOLU "my tube" TUBE 0. 10. 20. G4_WATER:P POSZ 5 :PLACE "my tube" 1 world R00 0. 0. -$POSZ:VOLU sphere ORB 5. G4_Si :PLACE sphere 1 "my tube" R00 0. 1. $POSZ

15. Movements • User can move any volume with user commands • Several movements of the same volume or different ones can be set at the same time • Displacements, rotations or text file that allows to define any movement • Every N events or every interval of time • N times or forever • Offset can be defined Sinusoidal movement simulated with GAMOS

16. Particle generator • C++ code • The usual GEANT4 way • Add one line to transform your class in a plug-in DEFINE_GAMOS_GENERATOR(MyGenerator); so that you can select it in your input macro /gamos/generator MyGenerator GAMOS generator • Wide range of particles: gamma, e-, e+, proton, neutron, ions, isotopes, … • Combine any number of particles in one event (history) • For each particle or isotope user may select by user commands any combination of time, energy, position and direction distributions • Geant4 radioactive decay can be selected with user command • Decay branching ratios and energy spectra of all isotopes /gamos/generator/timeDist source GmGenerTimeDecay /gamos/generator/directionDist source GmGenerPositionDiscGaussian 1.*mm

17. Generator distributions • Common medical physics distributions are available POSITION: Point Square Rectangle Disc DiscGaussian LineSteps InVolumes (sphere/box/tube_section/ellipsoid) InVolumeSurfaces (sphere/box/tube_section/ellipsoid) InVolumesGeneral (any solid) VoxelPhantomMaterials DIRECTION: Random Const Cone ENERGY: Constant Gaussian RandomFlat BetaDecay ConstantIsotopeDecay TIME: Constant Decay POSITION & DIRECTION: InVolumeSurfaceTowardsCentre TowardsBox • User can create its own distribution and select it with a user command (plug-in) F18 decay energy

18. Generator distributions EXAMPLE: Select as primary particle a gamma of 5.75 MeV: /gamos/generator/addSingleParticleSource MySource gamma 5.75*MeV Select as position uniformly distributed into the volume called ‘source’: /gamos/generator/positionDist MySource GmGenerDistPositionInG4Volumes source Select as direction a cone with axis along X and aperture 3 degrees: /gamos/generator/directionDist MySource GmGenerDistDirectionCone 1. 0. 0. 3*deg

19. Physics GAMOS physics list • Selection of all options of GEANT4 electromagnetic physics with user commands • standard / low energy / Penelope • different multiple scattering models • gamma/electron – nuclear interactions • atomic deexcitation • X-ray reflection (new development by Zhentian Wang, Dept. of Physics Engineering, Tsinghua University, China) • All physics lists in Geant4 code can be selected through a user command • Production cuts: can be set different for each volume region /gamos/physics/setCuts collimator 1.*mm 1*mm /gamos/physics/setCuts MLCs 10.*mm 1*mm • User limits: (min Ekin, min range, max track length, max step length, max time of flight) can be set different for each volume Utility to automatically optimize production cuts and user limits (see User’s Guide)

20. Extracting • detailed data

21. Obtaining detailed data • Most users are researchers, it is not enough to provide some final results, like dose distribution, or PET event classification table • Want to have a deep understanding of what happens in the simulation • Want to have the capability to evaluate the reliability of the results • Want to choose the best physics configuration For example: • Which is the dose contributed by the electrons that reach the phantom in a gamma accelerator treatment? • How many gammas traverse completely the jaws? How much energy they lose? • What is the length travelled by electrons produced by Compton interactions in a crystal? • … • Something we GAMOS developers have never imagined…

22. GAMOS data • Almost 200 types of data can be obtained at each step, track, event or run • event ID, trackID, particle, process, kinetic energy, energy lost, energy deposited, step length, number of secondary's, position (X/Y/Z/R2/R3/phi/theta), direction, change in position, change in angle, energy secondary/energy primary, angle secondary – primary, etc. • Also composed data, using arithmetic expressionssqrt(2.*FinalLocalPosX*FinalLocalPosY) • With a few user commands data can be used to fill histograms (1D, 2D, profile), be dumped in text or binary files, or be printed on the screen • EXAMPLE: Plot at each track step the logarithm of the initial kinetic energy and also the final X position vs final Y position /gamos/userAction GmStepDataHistosUA /gamos/setParam GmStepDataHistosUA:DataList log10(InitialKineticEnergy) FinalPosX.vs.FinalPosY

23. Filtering A filter is an algorithm that accepts or reject a step or track ● > 100 filters in GAMOS Starts/Ends/Enters/Exits/Traverses/In LogicalVolume/PhysicalVolume/Touchable/Region, or any of their children, particle type, process type, primary or secondary, charged or neutral, etc. Composed filters: ● AND, OR, XOR, Inverse, OnSecondary (apply filter to secondary tracks created in a step), with history (pass filter in any previous track step or any ancestor previous track step), etc. Filter on numeric data:if the data is inside an interval /gamos/filter energyF GmNumericDataFilter InitialKineticEnergy 100.*keV 1.*MeV Filter on string data:if the data is included in a list /gamos/filter scatteringF GmStringDataFilter Process compt rayleigh

24. Filtering example Dump in a text file the energy lost by the tracks in all kind of jaws # only when track steps happen in JAWS_X or JAWS_Y /gamos/filter jawsF GmInLogicalVolumeFilter MY_JAWS* # Select what to dump into text file /gamos/setParam GmTrackDataTextFileUA_jawsF:DataList EventID TrackID Particle AccumulatedEnergyLost # Include a header line in text file /gamos/setParam GmTrackDataTextFileUA_jawsF:WriteHeader 1 # Dump data in text file /gamos/userAction GmTrackDataTextFileUAjawsF HEADER: 4,"EventID","TrackID","Particle","AccumulatedEnergyLost" 2,2,"e-",0.0198374 40,2,"e-",0.151526 63,2,"e-",0.188402 119,52,"e-",0.494147 119,49,"gamma",0.514096 119,46,"gamma",4.02731

25. Classifying A classifier is an algorithm that classifies a step or track information, returning an index ● ByParticle, ByKineticEnergy, ByNAncestors, ByLogicalVolume, ByPhysicalVolume, ByTouchable, ByRegion, ByProcess, ByMaterial, etc. Compound classifiers:return an index combining two or more classifiers /gamos/classifier particleAndProcessC GmClassifierByParticle GmClassifierByProcess Classifier on numeric data:different index for each interval /gamos/classifier energyC GmClassifierByNumericData InitialKineticEnergy 0. 10.*MeV 1*MeV Classifier on string data:different index for each data value /gamos/classifier processC GmClassifierByStringData CreatorProcess

26. Classifying example Plot for each Z layer the X vs Y position only when an interaction occurred in the “phantom” volume # only when track steps happen in phantom /gamos/filter InPhantomF GmInLogicalVolumeFilter phantom # classify by the Z position (a different histogram for each Z interval) /gamos/classifier ZposC GmClassifierByNumericData FinalPosZ -100. 100. 200./10 # select what to plot /gamos/setParam GmStepDataHistosUA_InPhantomF_ZposC:DataList FinalPosX.vs.FinalPosY # build histograms at each step using filter and classifier gamos/userAction GmStepDataHistosUA InPhantomF ZposC -100 < Z < -80 -80 < Z < -60 -60 < Z < -40

27. Scoring • Scoring may be an important part of a simulation powerful and flexible framework developed, fully based on user commands: • Many possible quantities can be scored in one or several volumes (based on Geant4 scorers) ● Dose ● Deposited energy ● Flux (in/out/passage) ● Current (in/out/passage) ● Charge ● Step length ● Number of particles ● Number of interactions ● Number of 2ary particles ● Number of steps ● Minimum kinetic energy ● Kerma • For each scored quantity one of several filters can be used • only electrons, only particles with energy in a given interval, … • Several ways to classify the different scores • One different score for each volume copy, or volume name, or energy bin, … • Results can be printed in one or several formats for each scored quantity • Standard output, text/binary file, histograms • Scoring can be made in real or in parallel worlds • All scored quantities can be calculated with/without errors • All scored quantities can be calculated per event or per run • Taking into account correlations from particles from same event

28. Scoring • Everything managed with user commands For example: score the number of interactions (collisions) of gammas in each volume copy - Select all volumes for scoring /gamos/scoring/createMFDetector nInterDet * - Score number of interactions /gamos/scoring/addScorer2MFD nInterScorer GmG4PSNofCollision nInterDet - One score per each physical volume /gamos/scoring/assignClassifier2Scorer GmClassifierByPhysicalVolume nInterScorer - Only score for gammas /gamos/scoring/addFilter2Scorer GmGammaFilter nInterScorer OUTPUT: MultiFunctionalDet: nInterDet PrimitiveScorer: nInterScorer Number of entries= 6 index: target:1 = 2.111 +-(REL) 0.024593555 index: primary collimator:1 = 1.066 +-(REL) 0.049482843 index: flattening filter:1 = 0.028 +-(REL) 0.28856012 index: jaws_X:1 = 0.003 +-(REL) 1 index: jaws_X:2 = 0.006 +-(REL) 0.6231189 index: expHall:0 = 0.002 +-(REL) 0.70675279

29. Sensitive Detectors • To produce hits in GEANT4 a user has to (with C++): • Define a class inheriting from G4VSensitiveDetector • Associate it to a G4LogicalVolume • Create hits in the ProcessHits method • Clean the hits at EndOfEvent In GAMOS you can do all this with a user command /gamos/assocSD2LogVol SD_CLASS SD_TYPE LOGVOL_NAME • SD_CLASS: the G4VSensitiveDetector class • SD_TYPE: an identifier string, so that different SD/hits can have different treatment Users can create her/his own SD class and make it a plug-in DEFINE_GAMOS_SENSDET(MySD);

30. Hits Bi fluorescence peak 511 keV photoelectric (PE) peak backscattering peak escape peak E > 511 keV photoelectric peaks • A GAMOS hit has the following information: • ID of the sensitive volume copy • event ID • energy • time of the first E deposit • time of the last E deposit • position • list of all tracks that contributed • list of all ‘primary´tracks that contributed • list of all deposited energies • SD type • Users can create his/her own hit class • Hits can be stored in a job and read in another job • Format compatible with STIR software (PET reconstruction) /gamos/userAction GmHitsHistosUA

31. Detector effects • Measuring time • - A detector is not able to separate signals from different events if they come close in time • Dead time • When a detector is triggered, this detector (or even the whole group it belongs to) is not able to take data during some time • Paralizable (if a new event happens while it is dead the dead time restarts) • Non-paralizable (dead time is fixed, a new event has no effect) • Both can be set by the user in the input macro • A different time for each SD_TYPE /gamos/setParam SD:MeasuringTime:MySD_1 1000.*ns /gamos/setParam SD:MeasuringTime:MySD_2 250.*ns

32. Hits digitization and reconstruction • Framework to transform simulated hits digital signals reconstructed hits (energy/time/position) • Digitization and hits reconstruction are very detector specific : it is not possible to provide a general solution • GAMOS provides a few simple digitizers and hits reconstructor’s • 1 hit  1 digit • Merge hits close enough • Same set of sensitive volumes • Closer than a given distance • User may implement detector specific effects • Trigger • Pulse simulation • Sampling • Noise • … -- hits -- clusters of hits /gamos/userAction GmHitsHistosUA /gamos/userAction GmRecHitsHistosUA

33. Histograms • By default histograms are written in ROOT format • But many medical users have never heard of ROOT, while they are experts on Matlab, Origin, … • Own format: can write histograms in text files (CSV) without any external dependency • Read them in Matlab, Origin, Octave, GNUplot,…, even MS Excel /gamos/userAction GmGenerHistosUA /gamos/analysis/fileFormat CSV MS Excel graphics of energy of e+ from F18 decay

34. Parameter management • Many classes change their behaviour depending on the value of some parameters, that the user can set with user commands • GAMOS code use parameters extensively to maximise the flexibility • Always a default value in case the user does not care • Parameters can be numbers, strings, list of numbers or list of strings • Units and mathematical expressions can be used 10*cos(30*deg) (3-0.5*0.6)*becquerel • Wildcards can be used crystal_* *layer* • GAMOS provides a unique command to manage all these parameters • A parameter is defined in the input macro /gamos/setParam SD:EnergyResol 0.1*mm • Any class can use this value in any part of the code float enerResol = GmParameterMgr::GetInstance() ->GetNumericValue(“SD:EnergyResol”,0.); Over 300 parameters (all clearly documented)

35. Verbosity When a user finds a problem, or to better understand the simulation, she/he wants to get a lot of detailed output • But only for a few events, not always! Framework should provide a simple mechanism to activate/deactivate verbosity • Per event and/or track • Per simulation component (geometry, physics, PET, radiotherapy, …)

36. GAMOS verbosity • Verbosity controlled through user commands • Different verbosities for different simulation parts (geometry, physics, generator, PET, radiotherapy, …) • 6 levels of verbosity (silent, error, warning, info, debug, test) • Controlled per event • Selectable by user commands /gamos/verbosity GenerVerbosity info /gamos/verbosity PETVerbosity debug Verbosity managers are plug-ins: user can easily create its own one • GEANT4 tracking verbosity can be controlled by event and by track /gamos/setParam GmTrackingVerboseUA:EventMin 14632 /gamos/setParam GmTrackingVerboseUA:EventMin 14635 /gamos/setparam GmTrackingVerboseUA:TrackMin 1 /gamos/setparam GmTrackingVerboseUA:TrackMax 20

37. Optimisation

38. Optimisation: time studies • User commands to get CPU time study • By particle, energy bins, volume, region (or combination of them) Just add user commands: /gamos/classifier ParticleAndEnergyClassifier GmCompoundClassifier GmClassifierByParticle GmClassifierByEnergy /gamos/userAction GmTimeStudyUA ParticleAndEnergyClassifier gamma/0.001-0.01: User=0.01 Real=0 Sys=0 gamma/0.01-0.1: User=2.01 Real=2.45 Sys=0.27 gamma/0.1-1: User=19.12 Real=22.05 Sys=1.51 gamma/1-10: User=4.25 Real=5.4 Sys=0.3 e-/0.0001-0.001: User=0.07 Real=0.1 Sys=0 e-/0.001-0.01: User=0.54 Real=0.69 Sys=0.06 e-/0.01-0.1: User=4.71 Real=5.41 Sys=0.38 e-/0.1-1: User=15.59 Real=18.19 Sys=1.79 e-/1-10: User=82.83 Real=98.62 Sys=7.45

39. Optimisation: physics cuts • Two kind of physics cuts in Geant4: • Production thresholds • Limit secondary production for ionisation and bremsstrahlung (and e+e-) • User Limits: • Limit step size: better done tuning EM physics parameters • Limit track length • Limit time of flight • Kinetic energy/range minimum (stop particle if kinE/range below a limit): useful to kill particles when energy not enough to reach target • Better use range because more uniform for different regions • Better use kinetic energy because more efficient calculation • GAMOS solution: use kinetic energy inside but user provides it as range through a command • Different physics cuts per particle and region (group of volumes) can be defined through user commands

40. Optimisation: setting best cuts • Normally people want that optimising the cuts produces only few % (or less) difference in results • need to run very big statistics • Cuts can be very different for different particle/regions, and a cut in one particle/region may affect another • need to run many sets of cut values RESULT: usually people think it is too complicated • Do not optimise cuts • Just copy the cuts from someone, which usually are not the best one for them In GAMOS we propose an Automatic determination ofbest production cuts or user limits in one job

41. Extending the scripting language: • plug-in’s

42. GAMOS plug-in’s If I want some functionality that GAMOS does not have? • Best solution for biggest flexibility: plug-in’s What’s is a plug-in? It is the same in software that USB in hardware: The easiest way to add a new device (class), without touching the operative system (framework): no need to install a driver (modify framework classes) How it works in GAMOS: • If you want to use, for example, your own physics list instead of one of the GAMOS ones • Add one line in user’s code DEFINE_GAMOS_PHYSICS(MyPhysicsList); • Code is transformed into a plug-in • Automatically it may be selected with a user command /gamos/physics MyPhysicsList

43. GAMOS plug-in’s Advantages of plug-in’s: • No need to understand how GAMOS works internally (how GAMOS would invoke my code?) or modify GAMOS code • No need to recompile each time I want to alternate between the GAMOS component and my own one • GAMOS has no predefined components: user has full freedom in choosing components • Any user written code (geometry, primary generator, physics list, sensitive detector, user actions, …) can substitute any GAMOS component while still using the rest of GAMOS utilities • If you have a working application, you may still use it, while you take profit of the part of GAMOS you like • No restrictions on the way to do things: all Geant4 functionality is available to GAMOS users

44. Applications

45. PET/SPECT

46. PET/SPECT Hits energy LaBr3 • Any PET/SPECT detector can be simulated with simple text format • Utility to create simple PET geometries from a few parameters (number of crystals per block, number of block per ring, radius, ...) • Several sensitive detectors types selectable by user commands • User can create it own one and make it a plug-in • Automatic creation of hits with detailed information /gamos/assocSD2LogVol GmSDSimple SD_TYPE VOLUME_NAME • Writing hits into text or binary file • Retrieving them in next job (for doing parametric studies) • Format compatible with STIR software • Several algorithms for identifying first interaction for detectors with DOI info

47. PET/SPECT • Detector effects • Energy resolution • Time resolution • Dead time (paralizable / non-paralizable) • Measuring time • Coincidence time • Digitizer and reconstructed hits framework and examples ClearPET sensitivity profile optical photons in BrainPET

48. PET/SPECT • Table of PET/SPECT classification • True / scattered / random coincidences / scattered&random • PET/SPECT line far or close to vertex • Optionally merge hits if there is more than one because of Compton interactions • SPECT: gamma passed through collimator or not • Histograms of hits or reconstructed hits with a user command • Histograms of positron history • Histograms of PET/SPECT information PixelPET sensitivity of different event types ClearPET derenzo reconstruction

49. Radiotherapy

50. Radiotherapy • GAMOS framework covers the needs of an external beam radiotherapy simulation user: • Easy way to build accelerator geometries • Reading DICOM patient geometries • Common primary generator distributions • Writing/reading phase space files • Dose deposition in phantoms • Histograms and binary files • CPU Optimisation • Optimising cuts • Variance reduction techniques • Debugging tools