Simulation: Difference between revisions
No edit summary |
|||
(79 intermediate revisions by 7 users not shown) | |||
Line 4: | Line 4: | ||
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. | 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. Once the initial particles are set in motion, the details of the particle interactions with | 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 [http://geant4.web.cern.ch/geant4/ 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 [[RandomNumbersBasic|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|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 <code>JobConfig</code> Offline directory. Some specific topics have separate [[Code|documentation]]. Usually, an existing standard configuration can be modified for a custom purpose. Sample jobs can be run and the resulting events can be [[Printing|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. | |||
The remainder of this page presumes familiarity with the FHiCL run time configuration language and how it is used by Mu2e; this is documented on the [[FclIntro]] page. | |||
==Practical techniques== | ==Practical techniques== | ||
Detailed simulation can be computationally expensive. The collaboration may need to run grid jobs for months to create enough statistics | 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]]''' | * '''[[Staging]]''' | ||
** simulation may be | ** simulation may be interrupted at a certain point, saved, and started again later, for various practical reasons | ||
** Examples: stop simulation just outside the detectors, write it out | ** 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]]'' | <span id="SimMixing"></span> | ||
* '''[[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. | **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 | **Example: simulate all standard model physics in a event and a conversion electron separately, then mix them | ||
* '''[[TimeSim| | * '''[[TimeSim|Time simulation]]''' | ||
**The decay of particles may be turned off, or delayed | **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 | ** Example: stop simulation when muons stop in the target - decay/interact them in later jobs | ||
* '''Resampling''' | * '''Resampling''' | ||
**Take the output of a stage and simulate it many times with different random seeds. | **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. | **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 | * '''Filtering''' | ||
**Only write events which pass some criteria | **Only write events which pass some criteria | ||
**Examples: require high-momentum e-, or some activity in tracker before saving a cosmic ray event | **Examples: require high-momentum e-, or some activity in tracker before saving a cosmic ray event | ||
*Dropping | *'''Dropping''' | ||
**Do not write out all data products | **Do not write out all data products | ||
**Example: drop MC digi truth information, drop the digis after tracks are made | **Example: drop MC digi truth information, drop the digis after tracks are made | ||
*Multiple streams | * '''Multiple streams''' | ||
**Write some events to one output file, some to another | **Write some tracks from events to one output file, some to another | ||
**Example: events w/stopped muons to one file, rest to another | **Example: events w/stopped muons to one file, the rest of the particles to another | ||
*Variants | **Note: The way the streams are defined, one event can trigger multiple streams. | ||
**Run one stage with variations in the parameters | ***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. | ||
**Example: repeating background measurements with different | * '''Variants''' | ||
* Compression | **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 | **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 | **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 [[JobPlan|grid job planning]] and [[SimulationFCL]] for examples of how to make standard physics processes. | |||
==Products== | ==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=== | ===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=== | ||
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 [[#SimMixing|mixing]], the SimParticle numbering can get quite complex, see the [[Mixing|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=== | ===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=== | ===StatusG4=== | ||
The <code>Mu2eG4</code> module, which runs geant, produces an art product called <code>StatusG4</code>. This contains the size and CPU time for the event, and error flags. <code>FilterStatusG4</code> is a filter module which reads <code>StatusG4</code> 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 | |||
<code>StatusG4Analyzer</code> 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: | |||
<pre> | |||
services : { | |||
GeometryService : { inputFile : "JobConfig/cd3/geom_baseline.txt" } | |||
} | |||
</pre> | |||
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|geometry]] pages for how to specify, modify, and check geometries. | |||
Geometries may [[TutorialForGeometryVisualization|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 | |||
<pre> | |||
physics.producers.g4run.SDConfig.enableSD : [ tracker, calorimeter, calorimeterRO, CRV, virtualdetector ] | |||
</pre> | |||
A typical set of StepPointMC products | |||
<pre style="font-size:90%"> | |||
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 | |||
</pre> | |||
===sensitive volumes=== | |||
At times it is useful, for example, to count radiation dose in a geometry volume that is not already a sensitive detector. In this case, the arbitrary volume can be turned into a sensitive detector and a StepPointMC collection can be written for that volume. | |||
This requires a few internal code changes. Here is a recipe: | |||
My advice is to stick with whatever bouncing cap scheme you choose so you don't waste time making names match. | |||
1) You need to have a generic name for your new set of SDs, ''e.g.'' ProductionTargetCoreSection. | |||
Then add into Mu2eWorld::instantiateSensitiveDetectors() the equivalent of the following code: | |||
<pre> | |||
if(sdHelper_->enabled(StepInstanceName::ProductionTargetCoreSection)) { | |||
Mu2eSensitiveDetector* ProductionTargetCoreSectionSD = | |||
new Mu2eSensitiveDetector( SensitiveDetectorName::ProductionTargetCoreSection(), _config ); | |||
SDman->AddNewDetector(ProductionTargetCoreSectionSD); | |||
for(G4LogicalVolumeStore::iterator pos=store->begin(); pos!=store->end(); pos++){ | |||
G4String LVname = (*pos)->GetName(); | |||
if (LVname.find("ProductionTargetCoreSection") != std::string::npos) { | |||
(*pos)->SetSensitiveDetector(ProductionTargetCoreSectionSD); | |||
} | |||
} | |||
} | |||
</pre> | |||
but you're not done. | |||
2) go to Mu2eG4/inc/SensitiveDetectorName.hh and add | |||
<pre> | |||
static char const * ProductionTargetCoreSection(){ | |||
return StepInstanceName::name(StepInstanceName::ProductionTargetCoreSection).c_str(); | |||
} | |||
</pre> | |||
3) go to MCDataProducts/inc/StepInstanceName.hh and add "ProductionTargetCoreSection" to the enum and the list below, making sure to append your new one to the end and not changing the order. | |||
4) then in the fcl, | |||
<pre> | |||
physics.producers.g4run.SDConfig.enableSD : [ProductionTargetCoreSection] | |||
</pre> | |||
and danger, danger there is a different way of doing this that is still around in some fcl files. The sensitiveVolumes line will cause an exception to be thrown in Mu2eG4_module.cc that you can't catch with gdb. | |||
<pre> | |||
SDConfig: { | |||
sensitiveVolumes: [ProductionTargetCoreSection] | |||
TimeVD: { times: [] } | |||
} | |||
</pre> | |||
Finally, in an Analyzer module you might set up the following in your fcl: | |||
<pre> | |||
prodCore: { | |||
module_type: PrimaryProtonEnergyDumper TimeOffsets: {} | |||
hitsInputTag: "g4run:ProductionTargetCoreSection" | |||
} | |||
</pre> | |||
to get at the hits. | |||
===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 <code>Mu2eG4/src/constructVirtualDetectors.cc</code> | |||
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". | |||
<pre> | |||
Mu2eG4CommonCut: { | |||
type: inVolume | |||
pars: [ TrackerMother, CalorimeterMother ] | |||
write: crvStage1 | |||
} | |||
</pre> | |||
===notes on G4 Geometry=== | |||
Building geometries in G4 is notoriously non-trivial. A number of useful routines like nestTubs have been written to make your life simpler, but for anything past the simplest constructs you need to know more. | |||
I give a worked-out example for a piece of the Production Target. One big problem here is that the target is rotated at 14^o about the y-axis, but the mother volume is not, and for good reason. There are many ways to handle this, all dissatisfactory. | |||
The target has a core section, consisting of disks along the internal z-axis of the target. Those are very simple to handle and nestTubs works fine. The problem arises when the object is not along the z-axis, so displaced in x or y or both. The guts of G4 are described in their standard talks, e.g. | |||
<pre> | |||
geant4.web.cern.ch/sites/geant4.web.cern.ch/files/geant4/collaboration/workshops/users2002/talks/lectures/geobasics.pdf | |||
</pre> | |||
or | |||
<pre> | |||
geant4-userdoc.web.cern.ch/geant4-userdoc/UsersGuides/ForApplicationDeveloper/html/Detector/Geometry/geomSolids.html?highlight=g4multiunion | |||
</pre> | |||
There are two tricky things you need to keep in mind: | |||
1) CLHEP and G4 use active vs passive rotations. If you look in ProductionTargetGeom/inc/ProductionTarget.hh you will see: | |||
<pre> | |||
// "passive" rotation, used for placing the production target. | |||
// This is the inverse of protonBeamRotation. | |||
const CLHEP::HepRotation& productionTargetRotation() const { return _protonBeamInverseRotation; } | |||
</pre> | |||
which nicely sums up how confusing this can be. | |||
2) The main sequence is this: | |||
<pre> | |||
G4LogicalVolume* finWithCutoutLogical = new G4LogicalVolume(finWithCutoutSolid,prodTargetFinMaterial,name); | |||
G4VPhysicalVolume* finWithCutout = new G4PVPlacement(rotFinG4,finTranslation,finWithCutoutLogical, | |||
name, prodTargetMotherInfo.logical,0,finCopyNumber,false); | |||
// see below the comment on using Mu2e finishNesting instead of G4PVPlacement though | |||
</pre> | |||
and let's take a look at the penultimate line. rotFinG4 is the G4 style rotation matrix for an internal rotation of the object. Follow Euler angle rules. Since you passed Goldstein, this is tedious but not difficult. The finTranslation piece is the problematic part. It is the location of the object in the mother volume, which means you must calculate where that is in the mother yourself. For our Production Target, since the mother is along the magnetic-- not internal target -- z-axis, I had to apply the following line: | |||
<pre> | |||
CLHEP::Hep3VectoroffsetRelativeToCore(distanceFromCenter*TMath::Cos(currentFinAngle), | |||
distanceFromCenter*TMath::Sin(currentFinAngle),0.); | |||
offsetRelativeToCore *= tgt->productionTargetRotation().inverse(); | |||
</pre> | |||
and from 1) you know why I needed the inverse, since this offset is a CLHEP 3-vector and the production target rotation was inverted. | |||
Once you know this, it's not impossible to keep track and there's nothing that a few cout's won't cure. | |||
( When writing new Mu2e code do not use G4PVPlacement. Use Mu2e finishNesting instead. This will assure automatic invocation of Mu2eG4Helper call to add the volumes to the Mu2e VolumeInfo collection and set visibility attributes and the flag controlling the overlap/surface checks. Regarding "lookup token" used e.g., in finishNesting and Mu2e nest... functions, it is described in: GeometryService/inc/G4GeometryOptions.hh with most of the related global options set in Mu2eG4/geom/g4_visOptions.txt; In many cases there will be a Mu2e nest function in Mu2eG4/inc, use it together with objects in GeomPrimitives/inc. When a logical volume needs to be placed multiple times, set "placePV" to false when using e.g. nestBox function. ) | |||
Boolean volumes are also annoying. Make an object, and relative to the geometry OF THAT OBJECT do the subtraction, etc. Now here you must cheat a little, for two reasons: | |||
(1) you don't want the same edge to belong to two different objects. My problem was subtracting a cylindrical segment from the bottom of a box, cylinder from 0 to Pi, both the same thickness, and the bottom edge exactly aligned with the bottom of the box. Which volume does the planar piece at the bottom of the box belong to? Or the sides of the box? You would like to make G4 pick one. | |||
(2) the visualizer | |||
(say GDML) doesn't know if the z=0 plane (for example) belongs to one object or the other. So the view on your screen now depends on the angle at which you view the geometry. I just make the subtraction volume (for me) a tiny amount bigger than it needs to be and that solves both problems. | |||
<pre> | |||
G4Box* finTopBox = new G4Box("finTopBox",tgt->haymanFinThickness()/2.,finTopHalfHeight,tgt->thicknessOfGapPerSection().at(ithSection)/2.); | |||
G4Tubs* finTopCutout = new G4Tubs("finTopCutout",0.,gapRadius+.0001,tgt->haymanFinThickness()/2. + .0001,0.,2.*M_PI); | |||
</pre> | |||
I found a super-useful debugging trick with subtraction volumes is to change it to a union volume and displace the one you're subtracting so you can see them both. After you make sure everything is lined up except for the displacement, get rid of the displacement and switch it back to a subtraction volume. | |||
===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: | |||
<pre> | |||
physics.producers.g4run.SDConfig.timeVD : [ 700, 750 ] | |||
</pre> | |||
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 <code>CompressPhysicalVolumes</code> module (often labeled <code>CompressPV</code>) writes an art product which includes the index and the volume name. This product, called <code>PhysicalVolumeInfo</code>, is written in the SubRun record. | |||
==Generators== | |||
A generator produces the first particle in a simulation chain and write the <code>GenParticle</code> 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== | ==Mu2eG4== | ||
===Cuts=== | |||
===Collections=== | The module that controls geant is <code>Mu2eG4/src/Mu2eG4_module.cc</code>. 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. | ||
===Filtering=== | |||
===Configuration=== | |||
The first part of the Mu2eG4 module stanza defines the physics of the Geant4 simulation. Many definitions and defaults are in <code>Mu2eG4/fcl/prolog.fcl</code>. | |||
The authoritative description of configuration parameters for this module is not in a wiki, but can be printed by | |||
mu2e --print-description Mu2eG4 | |||
A typical fcl file should contain | |||
g4run: @local::g4run | |||
and override individual parameters as needed. | |||
Some of the settings inside @local::g4run are | |||
<pre> | |||
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 | |||
</pre> | |||
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 | |||
<pre> | |||
diagLevel : 0 | |||
trackingVerbosityLevel : 0 | |||
steppingVerbosityLevel : 0 | |||
</pre> | |||
The following will put the Geant4 geometry navigator in the so called check mode: | |||
<pre> | |||
physics.producers.g4run.debug.navigatorCheckMode : true | |||
</pre> | |||
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.: | |||
<pre> | |||
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) | |||
</pre> | |||
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. | |||
<pre> | |||
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 } | |||
</pre> | |||
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|staging]] simulation. | |||
Here is an example from the second stage of a protons-on-target (beam) job: | |||
<pre> | |||
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 | |||
} | |||
] | |||
} | |||
</pre> | |||
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. | |||
<pre> | |||
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 | |||
</pre> | |||
===Staging Inputs=== | |||
When simulation is stopped, written to a file and [[Staging|restarted]], the input to the second stage has to be specified in the second stage job, which is done with the <code>inputs</code> stanza. For example: | |||
<pre> | |||
inputs : { | |||
primaryType: "StepPoints" | |||
primaryTag: "muonBeamFilter:mubeam" | |||
simParticleNumberOffset: 100000 // safe b/c of g4.particlesSizeLimit in stage1 | |||
inputSimParticles: "muonBeamFilter" | |||
inputMCTrajectories: "" | |||
inputPhysVolumeMultiInfo: "compressPVmuonBeam" | |||
} | |||
</pre> | |||
Each stage will produce some new SimParticles. The number offset makes sure they have a unique SimParticle number. The rest of the parameters specify the collection written in the previous stage, to be continued in this stage. | |||
====Legacy staging config==== | |||
Older versions of Offline use <code>MultiStageParameters</code> instead of <code>inputs</code>. For example: | |||
<pre> | |||
MultiStageParameters : { | |||
simParticleNumberOffset: 100000 // safe b/c of g4.particlesSizeLimit in stage1 | |||
genInputHits: [ "muonBeamFilter:mubeam" ] | |||
inputSimParticles: "muonBeamFilter" | |||
inputMCTrajectories: "" | |||
inputPhysVolumeMultiInfo: "compressPVmuonBeam" | |||
} | |||
</pre> | |||
===Multi-threaded Mode=== | |||
The Mu2e Offline code supports multi-threaded running in the Geant4 module, <code>Mu2eG4/src/Mu2eG4MT_module.cc</code>. Other modules such as Generator and Filter modules can run parallel threads through art-supported multithreading. | |||
To run Mu2eG4 in multithreaded mode, the recommended configuration setup is to base your MT fcl file on a non-MT fcl file, | |||
include the non-MT fcl file, and add in the additional necessary parameters. Here is an example fcl file, <code>Mu2eG4/fcl/g4test_03MT.fcl</code>: | |||
<pre> | |||
# Configuration file for G4Test03MT | |||
#include "Mu2eG4/fcl/g4test_03.fcl" | |||
# Note that all parameters except the following should mostly be | |||
# set in the non MT file included above (which includes prolog files) | |||
physics.producers.g4run.module_type : "Mu2eG4MT" | |||
# These set the number of threads used in MT mode. | |||
# The number and threads and number of schedules | |||
# should be the same. | |||
services.scheduler.num_schedules : 5 | |||
services.scheduler.num_threads : 5 | |||
</pre> | |||
==Filtering== | |||
===FilterStatusG4=== | |||
The <code>Mu2eG4</code> module, which runs geant, produces an art product called <code>StatusG4</code>. This contains the size and CPU time for the event, and error flags. <code>FilterStatusG4</code> is a filter module which reads <code>StatusG4</code> 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 <code>g4status</code>) | |||
* as a veto to remove events with defects from an output stream (sometimes called <code>g4consistent</code>) | |||
<code>StatusG4Analyzer</code> 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 <code>extraHitInputs</code> 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 <code>virtualdetector</code> <code>StepPointMC</code> related to the interesting particles. | |||
An option <code>vetoParticles</code> is available to remove particles and their daughters that are flagged in a special SimParticle collection. A related option <code>vetoDaughters</code> 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. | |||
<span style="color:red">How are collections combined and named?</span> | |||
==Output Modules== | ==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: ["<font style="color:red">trigFilter</font>"] | |||
<font style="color:blue">outputCommands:</font> [ | |||
"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 <font style="color:red">trigFilter</font> result. | |||
* what art products are kept, determined by the <font style="color:blue">outputCommands</font>. The names go by the 4-paramter names as [[ReadProducts|described here]] | |||
* the output file name | |||
See also details of [[FclIOModules|output modules]]. | |||
[[Category:Computing]] | |||
[[Category:Code]] |
Latest revision as of 16:27, 15 October 2021
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.
The remainder of this page presumes familiarity with the FHiCL run time configuration language and how it is used by Mu2e; this is documented on the FclIntro page.
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 and SimulationFCL for examples of how to make standard physics processes.
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
At times it is useful, for example, to count radiation dose in a geometry volume that is not already a sensitive detector. In this case, the arbitrary volume can be turned into a sensitive detector and a StepPointMC collection can be written for that volume.
This requires a few internal code changes. Here is a recipe:
My advice is to stick with whatever bouncing cap scheme you choose so you don't waste time making names match.
1) You need to have a generic name for your new set of SDs, e.g. ProductionTargetCoreSection.
Then add into Mu2eWorld::instantiateSensitiveDetectors() the equivalent of the following code:
if(sdHelper_->enabled(StepInstanceName::ProductionTargetCoreSection)) { Mu2eSensitiveDetector* ProductionTargetCoreSectionSD = new Mu2eSensitiveDetector( SensitiveDetectorName::ProductionTargetCoreSection(), _config ); SDman->AddNewDetector(ProductionTargetCoreSectionSD); for(G4LogicalVolumeStore::iterator pos=store->begin(); pos!=store->end(); pos++){ G4String LVname = (*pos)->GetName(); if (LVname.find("ProductionTargetCoreSection") != std::string::npos) { (*pos)->SetSensitiveDetector(ProductionTargetCoreSectionSD); } } }
but you're not done.
2) go to Mu2eG4/inc/SensitiveDetectorName.hh and add
static char const * ProductionTargetCoreSection(){ return StepInstanceName::name(StepInstanceName::ProductionTargetCoreSection).c_str(); }
3) go to MCDataProducts/inc/StepInstanceName.hh and add "ProductionTargetCoreSection" to the enum and the list below, making sure to append your new one to the end and not changing the order.
4) then in the fcl,
physics.producers.g4run.SDConfig.enableSD : [ProductionTargetCoreSection]
and danger, danger there is a different way of doing this that is still around in some fcl files. The sensitiveVolumes line will cause an exception to be thrown in Mu2eG4_module.cc that you can't catch with gdb.
SDConfig: { sensitiveVolumes: [ProductionTargetCoreSection] TimeVD: { times: [] } }
Finally, in an Analyzer module you might set up the following in your fcl:
prodCore: { module_type: PrimaryProtonEnergyDumper TimeOffsets: {} hitsInputTag: "g4run:ProductionTargetCoreSection" }
to get at the hits.
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 }
notes on G4 Geometry
Building geometries in G4 is notoriously non-trivial. A number of useful routines like nestTubs have been written to make your life simpler, but for anything past the simplest constructs you need to know more.
I give a worked-out example for a piece of the Production Target. One big problem here is that the target is rotated at 14^o about the y-axis, but the mother volume is not, and for good reason. There are many ways to handle this, all dissatisfactory.
The target has a core section, consisting of disks along the internal z-axis of the target. Those are very simple to handle and nestTubs works fine. The problem arises when the object is not along the z-axis, so displaced in x or y or both. The guts of G4 are described in their standard talks, e.g.
geant4.web.cern.ch/sites/geant4.web.cern.ch/files/geant4/collaboration/workshops/users2002/talks/lectures/geobasics.pdf
or
geant4-userdoc.web.cern.ch/geant4-userdoc/UsersGuides/ForApplicationDeveloper/html/Detector/Geometry/geomSolids.html?highlight=g4multiunion
There are two tricky things you need to keep in mind:
1) CLHEP and G4 use active vs passive rotations. If you look in ProductionTargetGeom/inc/ProductionTarget.hh you will see:
// "passive" rotation, used for placing the production target. // This is the inverse of protonBeamRotation. const CLHEP::HepRotation& productionTargetRotation() const { return _protonBeamInverseRotation; }
which nicely sums up how confusing this can be.
2) The main sequence is this:
G4LogicalVolume* finWithCutoutLogical = new G4LogicalVolume(finWithCutoutSolid,prodTargetFinMaterial,name); G4VPhysicalVolume* finWithCutout = new G4PVPlacement(rotFinG4,finTranslation,finWithCutoutLogical, name, prodTargetMotherInfo.logical,0,finCopyNumber,false); // see below the comment on using Mu2e finishNesting instead of G4PVPlacement though
and let's take a look at the penultimate line. rotFinG4 is the G4 style rotation matrix for an internal rotation of the object. Follow Euler angle rules. Since you passed Goldstein, this is tedious but not difficult. The finTranslation piece is the problematic part. It is the location of the object in the mother volume, which means you must calculate where that is in the mother yourself. For our Production Target, since the mother is along the magnetic-- not internal target -- z-axis, I had to apply the following line:
CLHEP::Hep3VectoroffsetRelativeToCore(distanceFromCenter*TMath::Cos(currentFinAngle), distanceFromCenter*TMath::Sin(currentFinAngle),0.); offsetRelativeToCore *= tgt->productionTargetRotation().inverse();
and from 1) you know why I needed the inverse, since this offset is a CLHEP 3-vector and the production target rotation was inverted. Once you know this, it's not impossible to keep track and there's nothing that a few cout's won't cure.
( When writing new Mu2e code do not use G4PVPlacement. Use Mu2e finishNesting instead. This will assure automatic invocation of Mu2eG4Helper call to add the volumes to the Mu2e VolumeInfo collection and set visibility attributes and the flag controlling the overlap/surface checks. Regarding "lookup token" used e.g., in finishNesting and Mu2e nest... functions, it is described in: GeometryService/inc/G4GeometryOptions.hh with most of the related global options set in Mu2eG4/geom/g4_visOptions.txt; In many cases there will be a Mu2e nest function in Mu2eG4/inc, use it together with objects in GeomPrimitives/inc. When a logical volume needs to be placed multiple times, set "placePV" to false when using e.g. nestBox function. )
Boolean volumes are also annoying. Make an object, and relative to the geometry OF THAT OBJECT do the subtraction, etc. Now here you must cheat a little, for two reasons:
(1) you don't want the same edge to belong to two different objects. My problem was subtracting a cylindrical segment from the bottom of a box, cylinder from 0 to Pi, both the same thickness, and the bottom edge exactly aligned with the bottom of the box. Which volume does the planar piece at the bottom of the box belong to? Or the sides of the box? You would like to make G4 pick one.
(2) the visualizer (say GDML) doesn't know if the z=0 plane (for example) belongs to one object or the other. So the view on your screen now depends on the angle at which you view the geometry. I just make the subtraction volume (for me) a tiny amount bigger than it needs to be and that solves both problems.
G4Box* finTopBox = new G4Box("finTopBox",tgt->haymanFinThickness()/2.,finTopHalfHeight,tgt->thicknessOfGapPerSection().at(ithSection)/2.); G4Tubs* finTopCutout = new G4Tubs("finTopCutout",0.,gapRadius+.0001,tgt->haymanFinThickness()/2. + .0001,0.,2.*M_PI);
I found a super-useful debugging trick with subtraction volumes is to change it to a union volume and displace the one you're subtracting so you can see them both. After you make sure everything is lined up except for the displacement, get rid of the displacement and switch it back to a subtraction volume.
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
.
The authoritative description of configuration parameters for this module is not in a wiki, but can be printed by
mu2e --print-description Mu2eG4
A typical fcl file should contain
g4run: @local::g4run
and override individual parameters as needed. Some of the settings inside @local::g4run are
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 inputs
stanza. For example:
inputs : { primaryType: "StepPoints" primaryTag: "muonBeamFilter:mubeam" simParticleNumberOffset: 100000 // safe b/c of g4.particlesSizeLimit in stage1 inputSimParticles: "muonBeamFilter" inputMCTrajectories: "" inputPhysVolumeMultiInfo: "compressPVmuonBeam" }
Each stage will produce some new SimParticles. The number offset makes sure they have a unique SimParticle number. The rest of the parameters specify the collection written in the previous stage, to be continued in this stage.
Legacy staging config
Older versions of Offline use MultiStageParameters
instead of inputs
. For example:
MultiStageParameters : { simParticleNumberOffset: 100000 // safe b/c of g4.particlesSizeLimit in stage1 genInputHits: [ "muonBeamFilter:mubeam" ] inputSimParticles: "muonBeamFilter" inputMCTrajectories: "" inputPhysVolumeMultiInfo: "compressPVmuonBeam" }
Multi-threaded Mode
The Mu2e Offline code supports multi-threaded running in the Geant4 module, Mu2eG4/src/Mu2eG4MT_module.cc
. Other modules such as Generator and Filter modules can run parallel threads through art-supported multithreading.
To run Mu2eG4 in multithreaded mode, the recommended configuration setup is to base your MT fcl file on a non-MT fcl file,
include the non-MT fcl file, and add in the additional necessary parameters. Here is an example fcl file, Mu2eG4/fcl/g4test_03MT.fcl
:
# Configuration file for G4Test03MT #include "Mu2eG4/fcl/g4test_03.fcl" # Note that all parameters except the following should mostly be # set in the non MT file included above (which includes prolog files) physics.producers.g4run.module_type : "Mu2eG4MT" # These set the number of threads used in MT mode. # The number and threads and number of schedules # should be the same. services.scheduler.num_schedules : 5 services.scheduler.num_threads : 5
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.