GRAPPA is a Geant4 application that simulates the interaction of an incoming particle beam (typically electron or photons) with an High-Z material. The basic GRAPPA geometry consists of an incoming beam and a square target of given thickness (foil) that can be rotated around the vertical axis. Detectors are placed onto the world boundary; the world is a sphere that surrounds the initial particles and the target. Final positions, momenta and proper time are registered by the detectors, that are triggered by the passage of
- Primary particles
- Electrons
- Positrons
- Pions
- Muons
- Photons
The output of a GRAPPA simulation consists of a histogram plot file with pre-analyzed particle statistics, a histogram ROOT file and an Ntuple ROOT file. In Geant4, an Ntuple represents a file where each column stores particular data (for example x coordinate or y momentum) of a given particle.
GRAPPA is a Geant4 (G4) application. A comprehensive guide on how to install G4 can be found on the G4 installation guide.
GRAPPA requires:
Geant4C++compiler (standard required C++-17)cmake
Additional features:
- Mutithreading support: if G4 is built with
-DGEANT4_BUILD_MULTITHREADED:BOOL=ON, GRAPPA will automatically execute using multiple threads, unless specified differently in the simulation parameters (e.g., in the macro file). - Visualization support: In order to produce visual outputs, G4 provided visualization drivers must be installed (e.g. OpenGL, RayTracer, QT6, etc...).
GRAPPA is built using cmake. From the source folder, GRAPPA is built and installed using the commands
cmake -S . -B build -DCMAKE_INSTALL_PREFIX=. -DGeant4_DIR={G4 cmake directory on the system} -DCMAKE_BUILD_TYPE="Release"
cmake --build build --target install
Notes:
- The
Geant4Config.cmakefile of the G4 installation should be specified via theDGeant4_DIRflag. For Geant4 versions < 11.1 it is located in${G4_base_install_dir}/lib/Geant4-versionotherwise in${G4_base_install_dir}/lib/cmake/Geant4. - The default installation directory is system-dependent. Please specify the installation directory via the flag
CMAKE_INSTALL_PREFIXor make sure you have the necessary privileges to install it in a system folder. The example build command, as shown before, installs GRAPPA in thebindirectory of the source code folder.
GRAPPA can either run using an interactive user interface (UI) or using a provided macro script. Before running, the G4 initialization script, provided by G4, must be sourced. It is located in the bin directory of the G4 installation directory and it should be sourced with
source ${G4_base_install_dir}/bin/geant4.sh
Launch GRAPPA with no extra argument using
./GRAPPA
and you will be prompted to the interactive session.
From there, you can control the simulation via the G4 UI commands.
GRAPPA can only be run in interactive mode if some visualization driver is provided.
In order to execute a macro, you should provide it when running the program as
./GRAPPA /path/to/macro
Examples of macros can be found in the script/run folder and they will be installed in the bin/run folder.
Macro files in GRAPPA act as input files, providing commands to initialize simulations, change parameters, set up the in-situ data analysis, change the physics, etc...
GRAPPA accepts all the basic Geant4 commands and introduces some specific ones. In this guide we only discuss some simple use cases, for a more complete manual please refer to the G4 application developers guide.
In a typical macro file, one wants to address some aspects of the simulation:
- Define the initial particle source
- Assemble the target geometry
- Modify the physics package
- Personalize the data analysis
- Set up the visualization
- Start the Run
- Others
Macro files are saved using the .mac extension and interpret
everything after a # as comments.
GRAPPA can shoot particle either from a General Particle Source (that is the default behavior), where the beam source statistics must be specified, or from a file of pre-generated list of positions and momenta. Choice between GPS and particles from file must be made before the run initialization. An example of GPS settings is
# =/=/=/=/= Particle source section =/=/=/=/=
# Set the particle source: must be done after run initialization.
# Particle type can be chosen between e- and gamma (for our use)
/gps/verbose 0 # verbose level of the GPS
/gps/particle e- #electron chosen as source
# Particle source spatial information
/gps/pos/type Beam # This is the source type, it represents a beam in the focal position
/gps/pos/centre 0 0 5 mm
/gps/pos/shape Circle
# The radius parameter defines a central flat plateau
/gps/pos/radius 0 mm
# The sigma_r parameter adds a gaussian uncertainty around the center
/gps/pos/sigma_r 10 um
# Angular distribution of the beam
/gps/ang/type beam1d
/gps/ang/sigma_r 0.002 rad # Beam divergence
# Energy distribution of the beam
/gps/ene/mono 10 GeV # Central beam energy
/gps/ene/type Gauss
/gps/ene/sigma 10 MeV # Beam energy spread
In case you want to use a particle source extracted from a file, you should provide an ASCII file with the phase space coordinates of the particles. The command to extract particles from a given ASCII file is
# Run this commands before initializing the run
/particlesource/particlesFromFile true
# ... (Other operations)
# Run initialization
/run/initialize
# ... (Other operations)
The structure of an ascii file is a six-column list of all the particles, in the form:
| x [mm] | y [mm] | z [mm] | px [MeV] | py [MeV] | pz [MeV] |
(separated by spaces or tabs). Other commands are available to modify the input from file and they can be called both before and after the initialization and even between runs
# Select the particle type that is being read from the file.
# If not specified, a geantino is created
/particlesource/setParticle gamma
# Provide the filename with the particles phase space
/particlesource/setFileName filename.dat
# Command to print the current setup
/particlesource/print
⚠️ If you request particles from a file, the particle type must be changed via the command/particlesource/setParticle <particlename>and not via the/guncommand. The/guncommand is not aware of the request to read particles from a file and therefore cannot update the particle container consistently.
If the file contains fewer particles than the ones requested for a run, GRAPPA will keep looping through the file relying on the fact that the random seed will be constantly updated.
GRAPPA supports two ways of handling simulation geometry. Geometry can be specified either via command line or by passing a geometry description file (GDML). The user chooses which geometry to use via macro commands
/geometry/useGDML true|false
Hardcoded geometry was the default system adopted by the code since its creation. GRAPPA defines a spherical world and a simple target, plus an thin "absorber" layer that coincides with the world. These three elements can be modified using macro commands. Commands to modify the geometry are internally defined and, as such, you will not find them in the G4 official guide.
They must be invoked before the run is initialized in order to avoid the loss of the geometry1.
# Modify geometry (custom functions)
# Modify the world material
# /geometry/world/material G4_AIR
# Modify the world radius
/geometry/world/radius 120 cm
# Modify the thickness of the absorber that encompasses the world
/geometry/absorber/thickness 0.2 mm
# Square target side length
/geometry/foil/size 10 mm
# Where the center of the target is placed
/geometry/foil/center 0 0 0 mm
# Target thickness
/geometry/foil/thickness 3 mm
# Target rotation around the y axis
/geometry/foil/rotation 45 deg
# Target material
/geometry/foil/material G4_W
The world is a sphere that surrounds everything else, so please ensure that its radius is large enough. Sensitive detectors are positioned on the world boundary (i.e. they are spherical), such that when a particle reaches them, it is registered before being killed.
GRAPPA also implements a more general system to handle the geometry based on GDML files. Using these files the user can change the geometry of the simulation entirely without needing to program C++ files and recompile the code from scratch. GRAPPA is sensitive to two particular auxiliary variables that can be specified within the geometry:
-
SensitiveDetectorTypewhich carries the value of a specific sensitive detector that should be attached to that volume. Refer to the section on sensitive detectors to check the supported types of detectors. The auxiliary value associated with the variable is the detector name. The detector name will generate ntuples with different names. Otherwise, one can associate multiple detectors of the same type with the same name, in which case all their ntuples will be merged in one. A single volume cannot have multiple sensitive detectors associated to it. This would require implementing a different detector class, not currently supported. -
Biasingwhich specifies viatrueorfalseif the volume is biased according the the biasing operation explained later.
An auxiliary line typically looks like: <auxiliary auxtype="SensitiveDetector" auxvalue="FinalAbsorber"/>.
The geometry is fully specified in the .gdml file, that is passed via
/geometry/gdmlFile filename.gdml
In addition to this, the user can choose to have GRAPPA generate an .obj file per volume defined in the GDML file via
/geometry/dumpObjGeometry true|false
/geometry/objDumpDirectory directory
where directory indicates the name of the directory where to store the list of files. This is a useful choice in case the user wants to visualize the object using external softwares (such as Paraview) or simply to load them into the Geant4 internal visualization system.
Important
The GDML geometry is only enabled if supported by the Geant4 installation. Make sure to compile Geant4 with -DGEANT4_USE_GDML=ON. If the GDML support is not available in Geant4, GRAPPA will default on the standard geometry definition that was used in the past.
Analysis commands allow to personalize both the in-situ histograms that GRAPPA generates at the end of the simulation and the output files for the raw data.
The complete list of generated histogram is the following:
# Histogram list
# 1D
# - 0 Final Energy of primaries
# - 1 Final Energy of positrons
# - 2 Final Energy of electrons
# - 3 Final Energy of photons
# - 4 Final Energy of pions
# - 5 Final Energy of muons
# - 6 Initial Energy of primaries
# - 7 Final theta of primaries
# - 8 Final theta of electrons
# - 9 Final theta of positrons
# - 10 Final theta of photons
# - 11 Final theta of pions
# - 12 Final theta of muons
# - 13 Final phi of primaries
# - 14 Final phi of electrons
# - 15 Final phi of positrons
# - 16 Final phi of photons
# - 17 Final phi of pions
# - 18 Final phi of muons
#
# 2D
# - 0 Initial primary transverse distribution
# - 1 Initial primary longitudinal distribution
# - 2 Initial primary angle distribution
# - 3 Initial X phase space of primaries
# - 4 Initial Y phase space of primaries
# - 5 Final primary transverse distribution
# - 6 Final electron transverse distribution
# - 7 Final positron transverse distribution
# - 8 Final photon transverse distribution
# - 9 Final pions transverse distribution
# - 10 Final muons transverse distribution
# - 11 Final primary angle distribution
# - 12 Final electron angle distribution
# - 13 Final positron angle distribution
# - 14 Final photon angle distribution
# - 15 Final pions angle distribution
# - 16 Final muons angle distribution
When passing an histogram command, the corresponding histogram id must be provided. Some example is
# Modifying the histogram 1D number 3 (Final energy of photons)
# Setting 200 bins, from 1 to 100 MeV
# `none` refers to the function that should be applied to the bin value
# `linear` refers to the binning scheme of the histogram
/analysis/h1/set 3 200 1 100 MeV none linear
# Modifying the histogram 2D number 16 (Final muons angle distribution)
# The parameters are analogous to the 1D case, applied to the x and y axis respectively
/analysis/h2/set 16 60 -100 100 mrad none linear 60 -100 100 mrad none linear
The complete list of generated ntuples is the following:
# - 0 Initial primary particles
# - 1 Final primary particles
# - 2 Final positrons
# - 3 Final electrons
# - 4 Final photon
# - 5 Final pions
# - 6 Final muons
# - 7 Number of primaries evolved
# - 8 Final particles (all particles not captured by the previous ntuples)
The name of the file containing the output can be changed too, by using the command
/analysis/setFileName Filename
File extension must be omitted as it's introduced automatically depending on the file type (ROOT is the default type for GRAPPA).
Plotting style "ROOT" can only be used if
Freetypewas available at compile time. The corresponding command is/analysis/plot/setStyle ROOT_default, which must be omitted in caseFreetypeis missing.
Each or all of the histograms and Ntuples can be activated and deactivated for plotting (histogram only) and being generated at all. Additional histograms can also be created on the fly when setting up the simulation.
/ntuplecontrol/list
that lists all the Ntuples, and
/ntuplecontrol/setDump ntupleID ifdump
that sets the activation status of the ntupleID ntuple.
From version v0.5.0, GRAPPA uses a custom physics list based on a combinaiton of default lists with the "basic" electromagnetic option (option 0).
In addition to the default options, the custom physics list replaces the G4StoppingPhysics approach to muon stopping with the dedicated G4MuonMinusAtomicCapture.
G4MuonMinusAtomicCapture is deactivated in GRAPPA in favor of the default G4MuonMinusCapture as in the last few runs we experienced some code crashes related to muonic atoms definitions.
The latter includes a special description of muonic atoms and of the mu- capturing process.
Moreover, the code already includes muon production via Bethe-Heitler.
For more information read the Physics Reference Manual.
G4EmExtraPhysics list could be controlled by commands such as
/physics_lists/em/GammaToMuons true
/physics_lists/em/PositronToMuons true
This is not the case anymore and these commands won't work as the G4EmExtraPhysics list was removed.
Of this list, we kept only the muon production via pair decay.
Now, all the commands have been replaced by custom ones:
/GRAPPA/physicsList/list
/GRAPPA/physicsList/setMuonPairCrossSectionFactor
/GRAPPA/physicsList/getMuonPairCrossSectionFactor
In future releases, we may consider a more detailed approach to enable fine-grained control on other physics list components.
GRAPPA implements two kinds of sensitive detectors.
-
FinalAbsorberrecords particles the first time they touch it. It records their position, momentum, energy, particleID, weight, interaction time, and creator process. In the current implementation, a few "major" particles are saved in dedicated NTuples, while the rest of the particles are saved in a separated one. -
StoppingDetectorrecords all the particles stopping inside it, where stopping is defined as particles that are killed (a particle that leaves the detector is never registered). The detector saves the stopping position, the particle energy when it was killed the particle ID and weight, the stopping time and the stopping process.
StoppingDetector defines a dedicated messenger that is accessed via /GRAPPA/stoppingDetector/detectorName after initialization and that is used to define which particles and processes the detector us sensitive to. Use all to make the detector sensitive to all particles/processes.
The code provides biasing capabilities for muon production. In particular, we bias the muon pair production and pion decay by splitting the created particles in certain volumes. The user can provide independent biasing factors for the two processes using dedicated macro commands. When a process is biased, the weight of the products is reduced by a factor 1/splittingFactor to maintain a correct physical description on average.
Important
Since it is possible to increase the cross section of muon pair production, the biasing takes that into account and reduces the muon weight by the corresponding factor. This finally returns correct results on average and solves the need to manually adjust the muon weight in post-processing. However, since the weighting is now left to the system, we have to make sure that all the volumes where a significant number of muon pairs are generated are included in the biasing process. Otherwise, in some of this volumes the muon weight is not adjusted corrispondingly and the results can become unphysical.
We introduced a few new commands:
/GRAPPA/biasing/pionDecay/setSplittingFactor factor Sets the splitting factor for pion decay biasing
/GRAPPA/biasing/pionDecay/getSplittingFactor Prints the current splitting factor for pion decay biasing
/GRAPPA/biasing/muonPair/setSplittingFactor factor Sets the splitting factor for muon pair production biasing
/GRAPPA/biasing/muonPair/getSplittingFactor Prints the current splitting factor for muon pair production biasing
/GRAPPA/biasing/muonPair/setScaleMuonWeigthWithCrossSection true/false Set if the weight of the muons produced via pair production must be scaled with the artificial increase in cross section. This will violate the total energy conservation but it will return particle numbers proportional to the initial number of primaries. Warning, the scaling is only performed in biased volumes.
/GRAPPA/biasing/muonPair/getScaleMuonWeigthWithCrossSection Returns if the muon weights are scaled with the artificial cross section
Biasing is enabled by default, but it can be disabled (and/or enabled again) using
/GRAPPA/physicsList/setBiasing true|false, which is only accepted in the pre-initialization phase.
GRAPPA can easily generate plots of the system being simulated, both in interactive and in batch mode. Visual aspects can be modified with some commands:
# Use this open statement to create an OpenGL view:
/vis/open OGL 1280x720-0+0
# Disable auto refresh and quieten vis messages whilst scene and
# trajectories are established:
/vis/viewer/set/autoRefresh false
/vis/verbose errors
# Draw geometry:
/vis/drawVolume
# Set background color
/vis/viewer/set/background 0 0 0
# Draw smooth trajectories at end of event, showing trajectory points
# as markers 5 pixels wide:
/vis/scene/add/trajectories
/vis/modeling/trajectories/create/drawByCharge
/vis/modeling/trajectories/drawByCharge-0/default/setDrawStepPts true
/vis/modeling/trajectories/drawByCharge-0/default/setStepPtsSize 5
# Filter on particle type (show only electrons and positrons)
/vis/filtering/trajectories/create/particleFilter
/vis/filtering/trajectories/particleFilter-0/add e-
/vis/filtering/trajectories/particleFilter-0/add e+
#/vis/filtering/trajectories/particleFilter-0/add gamma
/vis/viewer/set/viewpointThetaPhi
/vis/viewer/zoomTo 0.9
/vis/viewer/flush
/random/setSeeds 2 8 0
/run/beamOn 3
/vis/ogl/export snapshot_1.pdf
The image is updated after every run.
Notes
- Advanced histogram plotting is available with the
Freetypelibrary: G4 must be built with-DGEANT4_USE_FREETYPE:BOOL=ON. - For the analysis of the
ROOT(default) output file, ROOT or the python package uproot must be installed.
The Run object is the Geant4 component that manages every program run. In G4, a run is the collection of some number of particle events, each from the generation to the end of each particle trajectory. The run commands specify the number of threads used, how often to print debug messages, the initial random seed, etc... After the run has been initialized, many aspects of the simulation such as the geometry, the numbers of threads or the particle source type are fixed and cannot be changed. Many others can still be modified afterwards. The command to initialize the run is
/run/initialize
which will print the summary information about the run being generated. After the initialization, data analysis, visualization and incoming particle statistics can still be modified. A run is started with
/run/beamOn N
where N is an integer number determining the number of extractions from the particle source. After the run is completed, the analysis is automatically performed and closed and the plots (if available) generated. Ending the run does not end the simulation. In fact, one can have multiple runs in the same simulation. As discussed, analysis, visualization and particle sources can be modified between runs. Re-runnning overwrites the analysis files previously generated, so it is a good practice to change the analysis file name between runs.
The simulation ends after all the lines of a macro file have been executed, then the codes cleans up the memory and return a successful message.
Tip: if you don't remember the syntax for some GRAPPA commands, run a simulation in interactive mode. There, you can find a full list of commands with their usage cases.
We provide some ROOT scripts to facilitate the analysis of the final GRAPPA data.
Those are automatically installed in the bin/analysis folder as scriptname.C.
You can copy them in your output folder and execute them in ROOT as
>> root scriptname.C
For more insights and to learn how to customize the scripts, visit the ROOT manual.
Important
We will slowly drop support for pure ROOT scripts. We found that python scripts based on the uproot package are much more agile to write and easier to support.
Previous scripts may become out-of-date, not supporting newer changes in the ntuple structure and they are left as a base for future record.
We will slowly add new python scripts over time.