Simulation
This page is a draft, please help complete it!
This page page needs expert review!
Introduction
The simulation of particles in the detector are used for many purposes. Standard model physics is simulated to estimate the backgrounds we can expect. Conversion electrons are simulated to study the efficiency for detecting new physics processes. Other simulations might tell us if a detector component will receive a survivable radiation dose. Still others might be used to optimize detector design or reconstruction algorithms.
Simulations start with a generator step. This step creates the initial particles to simulate, for example, protons approaching the primary target or cosmic rays approaching the top of the detector. Our generators include protons hitting the primary target, cosmic rays, and particle guns. Once the initial particles are set in motion, the details of the particle interactions with magnetic fields and materials are simulated by the geant package, which has been used widely in HEP for years and is the industry standard. Geant steps particles forward in time and at each steps allows for particles to scatter, range out, decay or interact with the material and produce more particles to trace. The processes are drawn from physics models using random numbers to to determine quantities which are only known statistically, such as scattering angles, which of several processes might occur, and decay times. The careful use random number seeds insures that the simulated events are statistically independent. In this way, the code builds up each event consisting of many particles traced through the material in as much detail as is needed by the user's analysis.
The Mu2eG4 module, which is often labeled "g4run", controls geant and the process of creating and tracing particles, including energy deposits. After this step, modules for each detector are run to simulate the measurement of the energy deposit, including uncertainty and noise. The output are the digi products, ADC's and TDC's for detector channels. At this point, the use of random numbers, and the simulation process, is done. The digis are the products that the real data will contain, and the reconstruction process follows.
The following sections are high-level overviews of topics in simulation. Due to the complexity of simulation configuration, the myriad possibilities, and how code evolves, we only present an overview of concepts here, with some example fcl and code snippets. Specific examples of full configurations can be found in the JobConfig Offline directory. Some specific topics have separate documentation. Usually, an existing standard configuration can be modified for a custom purpose. Sample jobs can be run and the resulting events can be printed to make sure they have the right contents. For large production runs, an expert should review the final fcl. It also might make sense to make a new code release or tag the code.
Practical techniques
Detailed simulation can be computationally expensive. The collaboration may need to run grid jobs for months to create enough statistics; it all depends on the needs of the project, which may be quite varied. Several techniques have been developed to create and use simulation efficiently, especially minimizing CPU, redundant processing, disk space, and risk.
- Staging
- simulation may be interrupted at a certain point, saved, and started again later, for various practical reasons
- Examples: stop beam simulation just outside the detectors, write it out. Simulate cosmic rays to the outside of the tracker, stop and write the result out.
- Mixing
- simulate different parts of an event, for example a conversion electron and everything else, simulated up to a certain point, and mix these simulations together to one event with an electron and all the standard model noise.
- Example: simulate all standard model physics in a event and a conversion electron separately, then mix them
- Time simulation
- The decay of particles may be turned off, or delayed
- Example: stop simulation when muons stop in the target - decay/interact them in later jobs
- Resampling
- Take the output of a stage and simulate it many times with different random number seeds.
- Example, a fairly small set of stopped muons can be used over and over again. This works as long as the variations in the events are significant. In the example of stopped muons, each resampling decays a muon with the electron launched into a different direction. This leads to very different result for hits in the detector even though the same stopped muon was used.
- Filtering
- Only write events which pass some criteria
- Examples: require high-momentum e-, or some activity in tracker before saving a cosmic ray event
- Dropping
- Do not write out all data products
- Example: drop MC digi truth information, drop the digis after tracks are made
- Multiple streams
- Write some tracks from events to one output file, some to another
- Example: events w/stopped muons to one file, the rest of the particles to another
- Note: The way the streams are defined, one event can trigger multiple streams.
- Example: one primary proton from stage 1 can produce a neutron and a muon that end up in dsregion and mubeam streams respectively. These tracks might be later mixed in two different micro-bunches, breaking the correlation between daughter particles from the original primary proton.
- Variants
- Run one simulation stage with variations in the detector geometry or physics parameters
- Example: repeating background measurements with different calorimeter designs
- Compression
- Compress particle lists, removing unneeded entries
- Example: Compress large EM showers into smaller summary products; remove particles that do not leave any trace in a detector
- Error streams
- flag events as having errors or producing extremely long particle lists
- Example: cosmic ray events with more than 100K particles are stopped and written to the "truncated" output stream
See also practical grid job planning.
Products
The output of simulation includes high-level truth products which describe the event, and low-level products that describe each step of each relevant particle.
GenParticle
This product includes information about the type of generator and the particles produced by the generator. For primary protons, this will include the position and momentum of the primary particle.
SimParticle
Simparticle product includes the initial and final time, position, momentum, and type of all particles geant traced. It records the physics process that created or destroyed the particle. If a particle decays or interacts to produce more particles, then those new particles will point to the original particle as their parent.
If a simulation is stopped, when it is restarted, the code will find SimParticles that are still viable, and start tracing them again. Each particle will have a unique number. Typically, each simulation stage N will number its particles starting with (N-1)*100,000. If a dataset is produced through mixing, the SimParticle numbering can get quite complex, see the mixing page
The SimParticle product can become very large, especially in events with EM showers, such that it can create memory or storage issues. The SimParticle product is often purged of any particles that don't leave energy in a detector, or produce particles that leave energy in a detector.
StepPointMC
This product records the time, positions and energy loss of a geant step. It is only recorded while the particle is in a sensitive detector: a straw, a calorimeter cell, or CRV bar. These steps are used to create the digis, the products in the events from real data detector readout. A StepPointMC is also used as a way to record where a particle crossed a boundary of a virtual detector, for studies or staging, such as the virtual volume enclosing the DS, or the tracker.
MCTrajectory
This product will contain every geant step (position and time) that a particle takes. It can be large, so it is not usually created, except for cosmic ray particles (which can have complex paths) and special studies. What particles should have trajectories saved is specified in the stanza:
physics.producers.g4run.TrajectoryControl
StatusG4
The Mu2eG4 module, which runs geant, produces an art product called StatusG4. This contains the size and CPU time for the event, and error flags. FilterStatusG4 is a filter module which reads StatusG4 product to detect and flag event with errors or overflowing particle lists. The filter is used two ways:
- as a filter to collect events with errors into a special output stream
- as a veto to remove events with defects from an output stream
StatusG4Analyzer is an analysis module to histogram simulation operations quantities such as time and event sizes.
Geometry
For geant to run efficiently, it requires that all surfaces and volumes are described in its geometry hierarchy. The geometry is constructed using a combination of a SimpleConfig file and code. The geometry config file is determined by the fcl services entry:
services : {
GeometryService : { inputFile : "JobConfig/cd3/geom_baseline.txt" }
}
Many of the features of the geometry, such as materials and positions of detector components can be specified or overridden in the config file, however, any changes must be consistent with the capabilities of the code. Some examples are changing the materials in a collimator, of the exact z position of target foils. An example of change which requires a change to the code is adding a new foil.
Once the geometry is constructed, it is constant for the duration of the art exe run. If a simulation is staged, the geometry can change between stages. This can be helpful or disastrous. If a geometry is modified incorrectly, it might lead to volumes overlapping, which will usually cause geant to fail and crash. There are simple procedures to search points on a surface to see if there are in overlapping volumes. See the geometry pages for how to specify, modify, and check geometries.
Geometries may be visualized.
sensitive detectors
Sensitive detectors are a list geometry objects that we are particularly interested in. The basic sensitive detectors are places where energy is measured: straws, calorimeter cells and CRV bars. Each sensitive detector can produce an art product that contains StepPointMC's that occur in that sensitive detector. These StepPointMC's are used to create digi's which are the energies as read out by the detector channels. If the staging of a simulation is designed so that a detector is not simulated that sensitive detector should be left disabled.
A typical sensitive detector configuration
physics.producers.g4run.SDConfig.enableSD : [ tracker, calorimeter, calorimeterRO, CRV, virtualdetector ]
A typical set of StepPointMC products
Friendly Class Name Module Label Instance Name Process Name Product ID mu2e::StepPointMCs g4run CRV AllPatRecReco 1:40 mu2e::StepPointMCs g4run tracker AllPatRecReco 1:43 mu2e::StepPointMCs g4run calorimeter AllPatRecReco 1:41 mu2e::StepPointMCs g4run calorimeterRO AllPatRecReco 1:42 mu2e::StepPointMCs g4run virtualdetector AllPatRecReco 1:44
sensitive volumes
Hi Bob, I am not sure how you would use the 'sensitive volumes’ method of adding the many Production Target Core Sensitive Detectors you want to add, unless you explicitly put all of the distinct names in the list of sensitive volumes in the fcl file. The code in SDHelper needs that list of distinct names.
You need to have a generic name for your new set of SDs, eg productiontargetcore which you enable in the fcl file in this way: physics.producers.g4run.SDConfig.enableSD : [productiontargetcore]
Then add into Mu2eWorld::instantiateSensitiveDetectors() the equivalent of the following code:
if(sdHelper_->enabled(StepInstanceName::tracker)) {
StrawSD* strawSD =
new StrawSD( SensitiveDetectorName::TrackerGas(), _config );
SDman->AddNewDetector(strawSD);
//loop over all of the LV names and find the ones we need
//set the SensitiveDetectors for these
for(G4LogicalVolumeStore::iterator pos=store->begin(); pos!=store->end(); pos++){
G4String LVname = (*pos)->GetName();
//from ConstructTrackerTDR and constructTrackerv3
if (LVname.find("TrackerStrawGas_") != std::string::npos) {
(*pos)->SetSensitiveDetector(strawSD);
}
}//for
}//if tracker
I am assuming you have some Sensitive Detector class, something akin to StrawSD, for these new volumes.
virtual detectors
Virtual detectors are geometry volumes that do not correspond to real material, i.e. are virtual. They are used to create boundaries that are useful for staging or physics studies. Some examples include a volume just inside the DS, or a volume just outside the tracker. Whether a virtual detector is created depends on whether the physics detector is created. For example, the virtual detector volume around the tracker is created if the tracker is created. Available virtual detectors can be seen in Mu2eG4/src/constructVirtualDetectors.cc
Here is an example from a cosmics job fcl (for Mu2eG4 module) which stops tracing particles on the surface of the tracker and calorimeter and writes the StepPointMC's to a collection called "crvStage1".
Mu2eG4CommonCut: {
type: inVolume
pars: [ TrackerMother, CalorimeterMother ]
write: crvStage1
}
time virtual detectors
The time virtual detector is used to record the current steps of all SimParticles at a given time. A typical use might be to locate the global state of the event at the start of the readout time (usually about 700 ns from the proton bunch). It is controlled by the stanza:
physics.producers.g4run.SDConfig.timeVD : [ 700, 750 ]
and writes a StepPointMC collection labeled "timeVD".
persistence
SimParticles save an index which indicates what geometry volume the particle started, stopped, or interacted in. If the same geometry is recreated, then the index can be converted to a geant volume, and the volume can be asked its name and material, etc, which may be a critical part of studies. But since geometry can can customized, may evolve quickly, and may be different at different stages of a simulation, it not always easy to decode the index.
A partial solution to this issue is currently implemented. For each stage of simulation, the code looks into the SimParticle list and makes a list of indices that are used on interesting particles. A CompressPhysicalVolumes module (often labeled CompressPV) writes an art product which includes the index and the volume name. This product, called PhysicalVolumeInfo, is written in the SubRun record.
Generators
A generator produces the first particle in a simulation chain and write the GenParticle art product.
Generators include:
- Protons - produce a proton in the beamline approaching the proton target. This is the basis of many simulations, including the signal, and all the standard model backgrounds in the detector.
- Cosmics - generate a cosmic ray approaching the detector
- Stopped particles - simulate the decay of a stopped muon or pion. The positions of the stopped particles are taken from the proton simulation.
- Particle guns - produce a partcile of any type at an position with any momentum. This is useful for studying detector materials.
Please see the Generators page for details.
Mu2eG4
The module that controls geant is Mu2eG4/src/Mu2eG4_module.cc. It controls the configuration of the physics, stepping (moving a particle forward in time) and stacking (creating a particle and adding it to the list), when to stop simulation, and some features of the output products.
Configuration
The first part of the Mu2eG4 module stanza defines the physics of the Geant4 simulation. Many definitions and defaults are in Mu2eG4/fcl/prolog.fcl.
module_type: Mu2eG4 physics: @local::mu2eg4DefaultPhysics # use standard physics processes "ShieldingM" ResourceLimits: @local::mu2eg4DefaultResourceLimits # define how big particles lists can be TrajectoryControl: @local::mu2eg4NoTrajectories # do not save trajectories debug: @local::mu2eg4DefaultDebug # do not print visualization: @local::mu2eg4NoVisualization # do not run any graphics
Some useful options are to turn on trajectories (uses a lot of disk space).
TrajectoryControl: @local::mu2eg4DefaultTrajectories
Or request a different physics list. These lists determine how detailed the simulation is at each step and which physics models are used
physics.producers.g4run.physics.physicsListName: "ShieldingM"
More information on Geant4 physics lists can be found at: http://geant4-userdoc.web.cern.ch/geant4-userdoc/UsersGuides/PhysicsListGuide/html/index.html
The debug stanza has the ability to request to print various steps as they happen.
Diagnostic printouts
Using various parameters of the physics.producers.g4run.debug pset one can affect what gets printed and at what detail. Changing the following parameters will result in an increased printout
diagLevel : 0
trackingVerbosityLevel : 0
steppingVerbosityLevel : 0
The following will put the Geant4 geometry navigator in the so called check mode:
physics.producers.g4run.debug.navigatorCheckMode : true
and can be used to debug geometry construction when Geant4 tracking exceptions are encountered.
Step Limits
Step limits determine how big a step Geant4 takes. We may ask for smaller steps in sensitive detectors when want to have a very good model of how much energy is deposited or for visualization purposes in the vacuum, e.g.:
physics.producers.g4run.physics.bfieldMaxStep : 20. // mm; for step limiter in vacua; impacts tracking accuracy as well physics.producers.g4run.physics.strawGasMaxStep : -1.0 // mm; for straw step limiter; impacts tracking accuracy as well (set negative to disable)
See Mu2eG4/src/Mu2eWorld::constructStepLimiters() for more details.
Range Cuts
Range, or production cuts, determine energy below which no secondary electrons, positrons and gammas are generated in some electromagnetic processes. That energy is expressed in mm and is the range which a particle would have traveled before ranging out in a given material. Instead of creating the particles, their energy is deposited at the given position. The above cuts can also be imposed selectively for a given, so called, Geant4 geometrical region, i.e. a root volume in the volume hierarchy. The production cut assigned to protons is the production threshold of nuclei for hadron elastic processes. See below for examples.
physics.producers.g4run.physics.minRangeCut : 1. // mm
physics.producers.g4run.physics.protonProductionCut: 1.0
physics.producers.g4run.physics.minRangeRegionCuts: {
CalorimeterMother : 0.1
TrackerMother : 0.001 }
// {RegionName : mm }
The range cuts are described in more detail in Geant4 Book For Application Developers, section 2.4.2 and chapter 5.4 and sections mentioned therein, see http://geant4.web.cern.ch/support/user_documentation
We have not seen much difference between protonProductionCut 0 and 1 mm, though 0 produces more particles, especially excited nuclei.
Cuts and Collections
Cuts may be specified to indicate when to stop simulation of a particle. For example, the simulation might be stopped on the outside of a detector, as particles cross a plane, or satisfy kinematic cuts. Cuts may combined with a "union" (OR) or "intersection" (AND). A kineticEnergy cut passes when the particle has energy less than the cut. A cut on a plane passes when the particle moves into the region pointed to by the normal.
For a set of particles passing a cut, you can request to write to a specific StepPointMC collection. These features are fundamental to staging simulation.
Here is an example from the second stage of a protons-on-target (beam) job:
Mu2eG4CommonCut: { # stop when cuts are satisfied
type: union # there is a one cut here, which is the union (OR) of this list of two following cuts
pars: [
@local::mu2eg4CutDeltaElectrons, # first cut : particle=electron and E<1
{ # start of second cut
type: inVolume # passes when step ends in the following volumes
pars: [TS5Vacuum, DS1Vacuum, DS2Vacuum, DS3Vacuum ]
write: DSVacuum # write the inVolume particle StepointMC's to a collection with this name
}
]
}
Here is an example of stopping simulation on a pair of intersecting planes and write it out. The normal of the plane points to the region where the cut passes when the particle enters. For example, particles tend to move towards +z, so an x-y stopping plane would usually have a normal 0,0,1.
type: intersection
pars: [
{ type: plane
normal: [ 0, 0, 1 ]
point : [ 0, 0, -4851 ]
},
{ type: plane
normal: [ -1, 0, 0 ]
point : [ 1612., 0, 0 ]
}
]
write outplane
Staging Inputs
When simulation is stopped, written to a file and restarted, the input to the second stage has to be specified in the second stage job, which is done with the MultiStageParameters stanza. For example:
MultiStageParameters : {
simParticleNumberOffset: 100000 // safe b/c of g4.particlesSizeLimit in stage1
genInputHits: [ "muonBeamFilter:mubeam" ]
inputSimParticles: "muonBeamFilter"
inputMCTrajectories: ""
inputPhysVolumeMultiInfo: "compressPVmuonBeam"
}
Each stage will produce some new SimParticles. The number offset makes sure they have a SimParticle number which indicates which stage it produced it. The rest of the parameters specify the collection writing in the previous stage, to be continued in this stage.
Multi-threaded Mode
As of early June 2019, the Mu2e Offline code supports multi-threaded running only in the Geant4 module, Mu2eG4/src/Mu2eG4_module.cc. Work is currently being done to integrate the multithreading capabilities of art3 into the Offline code, so that the Generator, Mu2eG4, and Filter modules can run in MT mode.
The following generators have been modified to work with Mu2eG4 running in MT mode:
- StoppedParticleReactionGunN_module.cc (an example fcl file is here: Mu2eG4/fcl/g4test_03MT.fcl)
- CRYEventGeneratorN_module.cc (an example fcl file is here: JobConfig/cosmic/cosmic_s1_general_cryMT.fcl)
- PrimaryProtonGunN_module.cc (an example fcl file is here: JobConfig/beam/beam_g4s1MT.fcl)
The Generator module must produce a 'stash' of generated particles to be used by the threads of G4. If you want to use a Generator module not listed above, you need to modify the module code to make it compatible with Mu2eG4 running in MT mode. Reference any of the above modules for details.
To run Mu2eG4 in multithreaded mode, your fcl file must include these fcl parameters
stashSize : <an int> runinMTMode : true numberOfThreads : <an int>
The 'stashSize' is needed for both the Generator and Mu2eG4 module.
Filtering
FilterStatusG4
The Mu2eG4 module, which runs geant, produces an art product called StatusG4. This contains the size and CPU time for the event, and error flags. FilterStatusG4 is a filter module which reads StatusG4 product to detect and flag event with errors or overflowing particle lists. The filter is used two ways:
- as a filter to collect events with errors into a special output stream (sometimes called
g4status) - as a veto to remove events with defects from an output stream (sometimes called
g4consistent)
StatusG4Analyzer is an analysis module to histogram simulation operations quantities such as time and event sizes.
FilterG4Out
This module is used to prepare geant collections for output. It can combine, veto and compress collections of SimParticle, StepPointMC and MCTrajectory collections.
Special categories of particles may be flagged by modules, typically by adding them to a StepPointMC collection. These lists may be particles that passed a momentum cut, or passed through a plane, or entered the volume surrounding a detector. These collections are typically produced by specifying a set of cuts to the Mu2eG4 module, and writing the particles passing cuts to a collection. This collection serves as the primary "good" collection for the FilterG4Out module. The primary function of the module is then to find the SimParticles related to the flagged particles, and then remove all the rest from the collections. The point is that a user will be interested in the particles that are in the tracker and doesn't need to save all the little photons that showered in the shielding.
An option extraHitInputs allows the user to specify collections that should also be saved after filtered for the interesting particles as defined by the main collections. Here a typical use would be to save the virtualdetector StepPointMC related to the interesting particles.
An option vetoParticles is available to remove particles and their daughters that are flagged in a special SimParticle collection. A related option vetoDaughters removes only the daughters of the flagged particles, but not the particles themselves. This is often used to stop simulation after a muon has stopped in material.
The module can also be given a simple list of SimParticles to keep.
How are collections combined and named?
Output Modules
There will be one output module, with module type "RootOutput", for each art file output stream. There are three main configuration points.
Here is an example.
outputs: {
filteredOutput : {
module_type : RootOutput
SelectEvents: ["trigFilter"]
outputCommands: [
"drop *_*_*_*",
"keep mu2e::GenParticles_*_*_*",
"keep *_cosmicFilter_*_*",
"keep *_compressPV_*_*"
]
}
}
outputs.filteredOutput.fileName : "sim.owner.cd3-cosmic-g4s1-general.version.sequencer.art"
- how to select which events to save. Here is determined by the filtering of the trigFilter result.
- what art products are kept, determined by the outputCommands. The names go by the 4-paramter names as described here
- the output file name
See also details of output modules.