Single Particle Simulations with `ddsim`

Overview

Teaching: 30 min
Exercises: 20 min

Questions

How can I simulate single particles for detetor studies?

Objectives

Know where to find the available options for ddsim.

Understand the differences between steering files and command line options.

Know what the key output file collections are.

In this first episode we will go through the running of single particle events with ddsim, using the built-in event generator of ddsim. This is the quickest way so run some straightforward tests of the geometry and produce output hits in the detectors for further analysis.

Passing options to `ddsim`

The program ddsim is part of the DD4hep installation when it is compiled with Geant4 support. In the EIC standard environment eic-shell it is available and used for many simulations in the suite of continuous integration and benchmarking checks for the geometry. Simply entering ddsim --help will show the large wealth of options that can be passed to ddsim on the command line.

It is also possible to use a python steering file to control the simulations. This is often a more convenient approach for sharing settings with others or archiving them for later reference or to iteratively add functionality. However, it should be noted that this can also lead to divergence in ‘default’ running conditions when these are not propagated to other environments.

Exercise:

Consult the avialable options of ddsim by passing it the --help flag.

Find the appropriate option to specify the output file when using ddsim.

Use the --dumpSteeringFile flag to print out a default steering file, and redirect it into a file steering.py.

Available options in `ddsim`

There are several option blocks that can be used with ddsim. Those are best looked at in the steering.py file, where additional documentation is added.

SIM.action.* or --action.* options can be used to tune sensitive detector actions,
SIM.field.* or --field.* options affect the magnetic field steppers,
SIM.filter.* or --filter.* options can add filters to sensitive detectors,
SIM.gun.* or --gun.* options can set single particle gun settings,
SIM.physics.* or --physics.* options allow setting the physics list,
SIM.random.* or --random.* options can be used to fix the random seed.

Some options, such as SIM.physics.setupUserPhysics, take as argument a python function, so they can only be used inside the python steering files. We will come back to this later.

When using the steering file approach, it is often useful to remove all options which you will not change (this allows you to take advantage of updates to the ddsim command itself without being stuck on old default settings). In this case, you would simply start from a steering file that only contains:

from DDSim.DD4hepSimulation import DD4hepSimulation
from g4units import mm, GeV, MeV
SIM = DD4hepSimulation()

and which can be passed to ddsim with the --steeringFile flag:

$ ddsim --steeringFile steering.py

This steering file merely sets up the simulation object that can then be configured with settings that deviate from the default.

Exercise:

Compose a ‘minimal’ steering file that only contains the required header line.

Attempt to run this steering file and note what ddsim claims is missing (we will specify this on the command line next).

Running a first single-particle simulation

When we used the minimal steering file, ddsim pointed out that we did not specify the ‘geometry compact file’, nor the number of events and source of those events. In this sections we’ll specify the geometry and tell ddsim to use a particle gun.

The compact file is the entry point of our geometry, for which we must load the geometry environment first

$ source /opt/detector/setup.sh
$ ddsim --steeringFile steering.py --compactFile $DETECTOR_PATH/$DETECTOR_CONFIG.xml

Next, we will specify that we want ddsim to use the DD4hep particle gun, with --enableGun (or -G), and that we want 10 events, with --numberOfEvents 10 (or -N 10):

$ ddsim --steeringFile steering.py --compactFile $DETECTOR_PATH/$DETECTOR_CONFIG.xml -G -N 10

Note: If you are running eic-shell directly on cvmfs, this may take a little while the first time. All large Geant4 data files are accessed for the first time and need to be retrieved.

When we run ddsim, by default it prints out some information for each generated event:

GenerationInit   INFO  +++ Initializing event 1. Within run:0 event 1.
Gun              INFO  Particle [0] mu-          Mom:10.000 GeV vertex:( 0.000  0.000  0.000)[mm] direction:( 0.000  0.000  1.000)
Gun              INFO  Shoot [0] 10.000 GeV mu- pos:(0.000 0.000 0.000)[mm] dir:( 0.000  0.000  1.000)
Gun              INFO  +-> Interaction [0] 10.000 GeV mu- pos:(0.000 0.000 0.000)[mm]
Gun              INFO  +++   +-> ID:  0 mu-          status:00000002 PDG:    13 Vtx:(+0.00e+00,+0.00e+00,+0.00e+00)[mm] time: +0.00e+00 [ns] #Dau:  0 #Par:0      
PrimaryHandler   INFO  +++++ G4PrimaryVertex at (+0.00e+00,+0.00e+00,+0.00e+00) [mm] +0.00e+00 [ns]
ParticleHandler  INFO  +++ Event 0 Begin event action. Access event related information.

You will notice that the particle gun has reverted to a default particle in a default direction: a 10 GeV muon in the positive z direction. That is, of course, not going to result in many hits in our ePIC detector… We will now take a closer look at some of the particle gun options.

  --gun.energy GUN.ENERGY
  --gun.particle GUN.PARTICLE
  --gun.multiplicity GUN.MULTIPLICITY
  --gun.phiMin GUN.PHIMIN
                        Minimal azimuthal angle for random distribution
  --gun.phiMax GUN.PHIMAX
  --gun.thetaMin GUN.THETAMIN
  --gun.thetaMax GUN.THETAMAX
  --gun.momentumMin GUN.MOMENTUMMIN
                        Minimal momentum when using distribution (default = 0.0)
  --gun.momentumMax GUN.MOMENTUMMAX
  --gun.direction GUN.DIRECTION
                         direction of the particle gun, 3 vector 
  --gun.distribution {uniform,cos(theta),eta,pseudorapidity,ffbar}
                        choose the distribution of the random direction for theta
                        
                            Options for random distributions:
                        
                            'uniform' is the default distribution, flat in theta
                            'cos(theta)' is flat in cos(theta)
                            'eta', or 'pseudorapidity' is flat in pseudorapity
                            'ffbar' is distributed according to 1+cos^2(theta)
                        
                            Setting a distribution will set isotrop = True
                            
  --gun.isotrop GUN.ISOTROP
                         isotropic distribution for the particle gun
                        
                            use the options phiMin, phiMax, thetaMin, and thetaMax to limit the range of randomly distributed directions
                            if one of these options is not None the random distribution will be set to True and cannot be turned off!
                            
  --gun.position GUN.POSITION
                         position of the particle gun, 3 vector 

While many of the options have straighforward names, others may be more confusing. The gun.distribution option is particularly relevant when we want to distribute single particle events over a range of angles. Depending on your needs, you may prefer one over the other, but in this tutorial we will simply use cos(theta) to throw uniformly on the unit sphere in the forward direction:

$ ddsim --steeringFile steering.py --compactFile $DETECTOR_PATH/$DETECTOR_CONFIG.xml -G -N 10 --gun.thetaMin "3*deg" --gun.thetaMax "45*deg" --gun.distribution "cos(theta)" --gun.momentumMin "1*GeV" --gun.momentumMax "10*GeV" --gun.particle "pi+"

Note: Avoid the use of the gun.energy option, since it is inherently more ambiguous than gun.momentum for massive particles (in the context of EIC).

Note how we pass arguments with units in this example. The double quotes are in many cases necessary on the command line to avoid having the * be expanded by your shell. In the python steering file they can be ommitted and one can simply write, e.g., SIM.gun.thetaMin = 3*deg.

When running this simulation above, note the output on the command line

GenerationInit   INFO  +++ Initializing event 1. Within run:0 event 1.
Gun              INFO  Particle [0] pi+          Mom:4.924 GeV vertex:( 0.000  0.000  0.000)[mm] direction:( 0.017  0.403  0.915)
Gun              INFO  Shoot [0] 10.000 GeV pi+ pos:(0.000 0.000 0.000)[mm] dir:( 0.000  0.000  1.000)
Gun              INFO  +-> Interaction [0] 10.000 GeV pi+ pos:(0.000 0.000 0.000)[mm]
Gun              INFO  +++   +-> ID:  0 pi+          status:00000002 PDG:   211 Vtx:(+0.00e+00,+0.00e+00,+0.00e+00)[mm] time: +0.00e+00 [ns] #Dau:  0 #Par:0      
PrimaryHandler   INFO  +++++ G4PrimaryVertex at (+0.00e+00,+0.00e+00,+0.00e+00) [mm] +0.00e+00 [ns]
ParticleHandler  INFO  +++ Event 0 Begin event action. Access event related information.

For historical reasons (which are being addressed), the relevant line is the one with Particle [0].

Exercise:

Simulate 10 events in the negative endcap, using an angle distribution that is uniform in cos(theta), and with an electron multiplicity of 2.

Add all gun options (but not the number of events) to the minimal steering file you constructed earlier and rename the steering file to ee_1GeV_10GeV_EndcapN.py and run the simulation again.

Verify that the momentum ranges and angular ranges are correct in the output.

Output files

Until now we have not bothered to check the output files (in case you were wondering and went exploring, you may have noticed that output went into dummyOutput.slcio). In this section we’ll explore how to write output files.

The command line option to use to specify the output file is the --outputFile option (or SIM.outputFile in the steering file). Depending on the extension of the output file, a specific output file format is chosen. The default output is in the slcio format, but we have standardized on the EDM4hep data model inside ROOT files. To choose this output file format, use a file extension .edm4hep.root. We could, for example, run the following command:

$ ddsim --steeringFile ee_1GeV_10GeV_EndcapN.py --compactFile $DETECTOR_PATH/$DETECTOR_CONFIG.xml --numberOfEvents 10 --outputFile ee_1GeV_10GeV_EndcapN_1e1.edm4hep.root

Note: The extension .edm4hep.root of the output file is important since ddsim infers the output file type from the extension.

Let’s take a look at the output file in ROOT (you can do this inside or outside the container, depending on your system and facility with opening ROOT browsers).

$ root -l ee_1GeV_10GeV_EndcapN_1e1.edm4hep.root
root [1] .ls
TFile**         ee_1GeV_10GeV_EndcapN_1e1.edm4hep.root  data file
 TFile*         ee_1GeV_10GeV_EndcapN_1e1.edm4hep.root  data file
  KEY: TTree    events;1        Events tree
  KEY: TTree    metadata;1      Metadata tree
  KEY: TTree    run_metadata;1  Run metadata tree
  KEY: TTree    evt_metadata;1  Event metadata tree
  KEY: TTree    col_metadata;1  Collection metadata tree

All EDM4hep files will have the same structure. We focus here on the events tree (if you see EVENT, in capital letters, you will need to ensure that you are indeed using the .edm4hep.root extension).

Note: In a future version of ddsim the output file format will change slightly, and the various metatadata trees will not be there anymore.

`MCParticles`

Let’s first focus on the MCParticles branch. This contains ‘truth’ information and exact Geant4 step output for selected tracks. It is filled independent of sensitive detectors defined in the geometry. It is most useful to analyze the initial state of the simulation (i.e. the final state of the event generator). In our case, the number of generated particles from the particle gun is always 2, which is what we expect to recover by looking at the @MCParticles.size() in the events tree:

root [2] events->Draw("@MCParticles.size()")

Wait a minute. What happened here? The MCParticles branch contains more than just the generated events. To include only the generated particles that were ‘thrown’ into the detector, we can select those with generatorStatus equal to 1.

root [3] events->Draw("MCParticles.PDG","MCParticles.generatorStatus==1")

Another approach which is sometimes useful to reduce the number of ‘other’ particles that are written is to pass the option --part.minimalKineticEnergy "1*TeV". This restricts the particles that are written to those that have an energy higher than 1 TeV. Initial state event generator particles are exempt from this requirement and are always written.

`*Hits`

The second set of branches that are important, and that are used as input to the event reconstruction framework, are the *Hits branches. These branches include the hits in sensitive detectors: either tracker hits or calorimeter hits. You will notice that these two different types of hits contain slightly different information.

Exercise:

Open the file ee_1GeV_10GeV_EndcapN_1e1.edm4hep.root which you just created in a local ROOT installation or online at https://root.cern/js/latest.

Plot the z component of the momentum for generated particles only and verify that this is indeed negative for particles going towards the negative endcap.

Verify that the total number of entries is consistent with the multiplicity and number of events you have simulated.

Plot the deposited energy of the hits in the endcap Ecal and compare the number of hits in the positive and negative endcaps.

Key Points

ddsim can be used with a particle gun to generate single particle simulations.

Running physics simulations with `ddsim`

Overview

Teaching: 30 min
Exercises: 20 min

Questions

How can I simulate events from physics event generators?

Objectives

ddsim can simulate HepMC3 events in DD4hep geometries.

npsim can be used as an alternative for simulations with optical photons.

We now move on to running simulations on HepMC3 event input files from (in this case) Pythia8.

Using centrally produced input files

The large input files for simulation campaigns are stored on S3, but the eic-shell environment contains the Minio cient mc to access them. We first need to set up our access credentials, though:

mc config host add S3 https://eics3.sdcc.bnl.gov:9000 $S3_ACCESS_KEY $S3_SECRET_KEY
mc ls S3/eictest/EPIC/Tutorials/pythia8NCDIS_10x100_minQ2=1_beamEffects_xAngle=-0.025_hiDiv.hepmc

Note: The values for $S3_ACCESS_KEY and $S3_SECRET_KEY will be provided in the tutorial, but are not provided here.

You will notice that the input file is large (GBs). For this tutorial we only need the first few thousand lines:

mc head -n 20000 S3/eictest/EPIC/Tutorials/pythia8NCDIS_10x100_minQ2=1_beamEffects_xAngle=-0.025_hiDiv.hepmc > pythia8NCDIS_10x100.hepmc

We can now specify this HepMC3 input file as input to ddsim:

ddsim --compactFile $DETECTOR_PATH/$DETECTOR_CONFIG.xml --numberOfEvents 10 --inputFiles pythia8NCDIS_10x100.hepmc --outputFile pythia8NCDIS_10x100.edm4hep.root

Instead of downloading files, we can also request events on-demand from the publicly accessible EIC XRootD server, but in this case we must use the hepmc3.tree.root input file extension:

ddsim --compactFile $DETECTOR_PATH/$DETECTOR_CONFIG.xml --numberOfEvents 10 --inputFiles root://dtn-eic.jlab.org//work/eic2/EPIC/Tutorials/pythia8NCDIS_10x100_minQ2=1_beamEffects_xAngle=-0.025_hiDiv.hepmc3.tree.root --outputFile pythia8NCDIS_10x100.edm4hep.root

Note: Many files on S3 under the S3/eictest/EPIC location are mirrored on XRootD under the root://dtn-eic.jlab.org//work/eic2/EPIC location. Note the use of the double slash in this URI!

Creating your own input files

Rather than relying on the centrally produced events, we can also create events ourselves. This gives us flexibility to run on and off certain event generator effects.

Head-on versus rotated collision frames: the “afterburner”

In this exercise, we will use Pythia8 to generate DIS neutral current interactions, but we could use other event generators as well. However, we have to pay attention to the reference frames in which interactions are generated. Most event generators are set up to generate events in the head-on collision frame of reference. This is not the reference frame in which beams collide at the EIC: the beams have a crossing angle of -0.025 mrad. In addition, there are beam energy smearing effects that cause the beam energies to deviated from the ‘exact’ values indicated in the settings: a 10 GeV electron beam contains in reality electrons with energies distributed around 10 GeV. To correct for crossing angle and beam energy smearing, we can modify the event generator or we can apply an “afterburner” which rotates and boosts events from the head-on fram into the correct frame. The afterburner is beyond the scope of this episode, but can be found at https://github.com/eic/afterburner.

Using Pythia8 with crossing angle and beam energy corrections

Rather than relying on the afterburner, we have modified Pythia8 to include the required corrections directly upon event generation. The steering code and input files can be found at https://github.com/eic/eicSimuBeamEffects, so we start with using git to obtain this code.

git clone https://github.com/eic/eicSimuBeamEffects

We can compile the code inside the eic-shell environment (which includes the Pythia8 event generator libraries that are used by this simulation):

cd eicSimuBeamEffects/Pythia8
make

After compilation, we can use the executable runBeamShapeHepMC.exe to generate events, but we need to provide some arguments:

$ ./runBeamShapeHepMC.exe
Wrong number of arguments
program.exe steer configuration hadronE leptonE xangle out.hist.root out.hepmc

The various steering files in steerFiles contain various beam conitions. Here we will use the 10 GeV electron on 100 GeV proton conditions in the high beam divergence setting (hiDiv), or the steering file dis_eicBeam_hiDiv_10x100. The hiDiv setting requires the configuration flag value 1 (as explained in the README.md file).

./runBeamShapeHepMC.exe steerFiles/dis_eicBeam_hiDiv_10x100 1 100 10 -0.025 \ 
  pythia8NCDIS_10x100_minQ2=1_beamEffects_xAngle=-0.025_hiDiv.hist.root \
  pythia8NCDIS_10x100_minQ2=1_beamEffects_xAngle=-0.025_hiDiv.hepmc

We can now run the output files through the ddsim simulation as before.

NPSim as an alternative to `ddsim`

Since some of the options that we pass to ddsim can only be provided through a steering file (such as python functions), or are otherwise cumbersome to provide on the command line, we provide npsim as a layer on top of ddsim that has these options pre-configured. This is as if you would take your steering file options and contribute them back to a central location for others to use them.

npsim can be easily interpreted (since it has sections that look exactly like the steering file). We can look at its python source code, located at /usr/local/bin/npsim.py in the eic-shell environment.

Currently npsim has the following additional options:

Cerenkov and optical photon physics are added through a python setup function,
an optical photon filter is created and added to the DRICH,
a modified tracking detector action which absorbs optical photon,
the minimal energy deposition in the tracker is set to zero.

You can run npsim exactly as you would run ddsim.

Exercise:

Rerun the previous Pythia8 simulation with npsim, and notice any difference in running time. Because of the addition of optical photon physics, the simulation will run more slowly.

Open the output file and verify that more hits (from optical photons with PDG code -22) are stored in the hits branches for the relevant RICH detector.

Key Points

ddsim or npsim are both able to simulate physics events.

EIC Tutorial: Simulations using ddsim and geant4

Single Particle Simulations with `ddsim`

Overview

Passing options to ddsim

Available options in ddsim

Running a first single-particle simulation

Output files

MCParticles

*Hits

Key Points

Running physics simulations with `ddsim`

Overview

Using centrally produced input files

Creating your own input files

Head-on versus rotated collision frames: the “afterburner”

Using Pythia8 with crossing angle and beam energy corrections

NPSim as an alternative to ddsim

Key Points

Passing options to `ddsim`

Available options in `ddsim`

`MCParticles`

`*Hits`

NPSim as an alternative to `ddsim`