Simulation: Difference between revisions

From Mu2eWiki
Jump to navigation Jump to search
 
(23 intermediate revisions by 6 users not shown)
Line 11: Line 11:
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 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; 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.
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]]'''  
Line 49: Line 50:
** Example: cosmic ray events with more than 100K particles are stopped and written to the "truncated" output stream
** 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]].
See also practical [[JobPlan|grid job planning]] and [[SimulationFCL]] for examples of how to make standard physics processes.


==Products==
==Products==
Line 109: Line 110:


===sensitive volumes===
===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:
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.
physics.producers.g4run.SDConfig.enableSD : [productiontargetcore]
 
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:
Then add into Mu2eWorld::instantiateSensitiveDetectors() the equivalent of the following code:
<pre>
<pre>
    if(sdHelper_->enabled(StepInstanceName::tracker)) {
    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>


        StrawSD* strawSD =
but you're not done.
        new StrawSD( SensitiveDetectorName::TrackerGas(), _config );
        SDman->AddNewDetector(strawSD);


        //loop over all of the LV names and find the ones we need
2) go to Mu2eG4/inc/SensitiveDetectorName.hh and add
        //set the SensitiveDetectors for these
<pre>
        for(G4LogicalVolumeStore::iterator pos=store->begin(); pos!=store->end(); pos++){
    static char const * ProductionTargetCoreSection(){
             G4String LVname = (*pos)->GetName();
      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>


            //from ConstructTrackerTDR and constructTrackerv3
Finally, in an Analyzer module you might set up the following in your fcl:
            if (LVname.find("TrackerStrawGas_") != std::string::npos) {
<pre>
              (*pos)->SetSensitiveDetector(strawSD);
        prodCore: {
            }
module_type: PrimaryProtonEnergyDumper TimeOffsets: {}
        }//for
hitsInputTag: "g4run:ProductionTargetCoreSection"
    }//if tracker
      }
</pre>
</pre>
I am assuming you have some Sensitive Detector class, something akin to StrawSD, for these new volumes.
to get at the hits.


===virtual detectors===
===virtual detectors===
Line 148: Line 178:
   }
   }
</pre>
</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===
===time virtual detectors===
Line 179: Line 269:


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 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>
<pre>
module_type: Mu2eG4
module_type: Mu2eG4
Line 290: Line 391:


===Staging Inputs===
===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>MultiStageParameters</code> stanza.  For example:
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>
<pre>
MultiStageParameters : {
MultiStageParameters : {
Line 300: Line 415:
}
}
</pre>
</pre>
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===
===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 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"


The following generators have been modified to work with Mu2eG4 running in MT mode:
# Note that all parameters except the following should mostly be
*StoppedParticleReactionGunN_module.cc (an example fcl file is here: Mu2eG4/fcl/g4test_03MT.fcl)
# set in the non MT file included above (which includes prolog files)
*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.
physics.producers.g4run.module_type : "Mu2eG4MT"


To run Mu2eG4 in multithreaded mode, your fcl file must include these fcl parameters
# These set the number of threads used in MT mode.
<pre>
# The number and threads and number of schedules
stashSize    : <an int>
# should be the same.
runinMTMode    : true
services.scheduler.num_schedules : 5
numberOfThreads : <an int>
services.scheduler.num_threads  : 5
</pre>
</pre>


The 'stashSize' is needed for both the Generator and Mu2eG4 module.


==Filtering==
==Filtering==

Latest revision as of 16:27, 15 October 2021

Construction.jpeg This page is a draft, please help complete it!

Expert.jpeg 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.