MCPL

MCPL (Monte Carlo Particle Lists) files contain lists of particle state information, and allows for easy storage and interchange of particles between various Monte Carlo simulation applications like Geant4, McStas, MCNP, etc. It is described in great detail in a published paper and recipes for how to use it with other programmes can be found on the MCPL website. Rather than being comprehensive, the present page will focus on MCPL tools and recipes that are most likely to be of interest for users of dgcode. Since the paper, the number of applications supporting MCPL has only continued to grow, but the basic concepts are still relevant.

Inspecting MCPL files

Run mcpltool myfile.mcpl.gz to print contents of files, including both file-level meta-data and the particle data itself. It can also be used to merge or repair files, or perform simple data extraction. Run mcpltool --help for more instructions.
Use pymcpltool --stats myfile.mcpl.gz to get printed a quick summary of the particle contents of a file.
Use sb_mcplextra_browse myfile.mcpl.gz to view a standard set of histograms of the content of the file (using the SimpleHists browser).
You can also add a filter expression (see below) and only get certain particles collected in your histograms. Example: sb_mcplextra_browse myfile.mcpl.gz where "is_neutron && ekin>1keV". Run sb_mcplextra_browse --help for more instructions.
If you need to, you can create custom C/C++ or Python code which opens an MCPL file and loops over the particles inside. You can find examples at the MCPL website here and here. In order to use MCPL inside dgcode, you should add USEEXT MCPL to your package dependencies in your pkg.info file.

MCPL files as simulation input

As documented here, MCPL files are easily used as input to your Geant4 simulation by using the G4MCPLPlugins.MCPLGen module in your sim script:

import G4MCPLPlugins.MCPLGen as Gen
gen = Gen.create()
gen.input_file = "mydata.mcpl.gz"

See the dedicated documentation for more options.

If the input file has a stat:sum:initial_ray_count value, all weights are divided by it, so that the weights in the Geant4 simulation are per initial ray of the simulation which wrote the file, and do not depend on its number of rays. This is for instance the case for files written by the MCPL_output component of McStas 3.10 and later (with weight_mode=1, the default), which store the weights multiplied by the number of rays. The value of the whole file is always used, also when only some of its particles are used (e.g. with skip_events or a limited number of events), so the results of jobs which each process a part of a file can be added up. Files without the value are used with their weights as stored. To always use the weights as stored, set gen.normalise_by_initial_ray_count = False.

MCPL files as simulation output

Normally this requires no code changes. Simply use the --mcpl flag to your sim script to capture all particles entering volumes with a certain name and write them into a newly created MCPL file. So if you have a volume called "FwdDet" and you wish to capture all particles reaching it into a newly created MCPL file, simply run your sim script like (assuming your package name is TriCorder):

$> sb_tricorder_sim [other options] --mcpl=FwdDet

A more advanced example, capturing any high energy gamma entering any volume called either "Detector" or "Sample" (see more about filter expressions further down on this page):

$> sb_mystuff_sim [other options] --mcpl="Detector,Sample where trk.pdgcode==22 and step.pre.ekin > 1MeV"

There are many other configuration possibilities, including attaching custom per-particle info to the stored particles (MCPL “user flags”). Run with the flag --mcpl=help for complete instructions:

$> sb_tricorder_sim --mcpl=help
Capture and write simulated particles to standardised MCPL files
by supplying one or more flags according to the following syntax
(note the usage of quotation marks)

    --mcpl="KEYWORDS [where CONDITION] [to FILENAME] [withflags FLAGS]"

KEYWORDS is a comma separated list containing both names of the
Geant4 volume in which G4Step's of particles will be considered
by the MCPLWriter instance being configured, and/or the follow-
ing list of special keywords:

   on_vol_entry  : Only treat G4Step's at volume entry [default]
   on_vol_exit   : Instead treat G4Step's at volume exit
   kill_always   : Tracks always killed upon treatment [default]
   kill_never    : Tracks never killed
   kill_filtered : Tracks killed only if filter passed
   all_vols      : Consider all volumes in active geometry
   opt_dp        : Store double-precision numbers in output.
   opt_pol       : Store polarisation info in output.
   opt_userflags : Store special user flags in output.
   opt_irc       : Store stat:sum:initial_ray_count of the MCPL input
                   file, with weights in the same convention as in
                   that file (only if each particle can be treated
                   as an independent ray).
   grabsrc       : Short-cut to grab just particles from the particle
                   generator. Implies kill_never and all_vols plus
                   prepends "trk.is_primary && step.stepnbr==1" to
                   CONDITION

CONDITION is an optional filter expression based on the G4ExprParser.

FILENAME is the name of the output file (defaults to "particles" if
not provided).

FLAGS is an expression generating a value for the 32bit integer user
flag field in the output file.

Examples:

1) Capture (and halt further simulation of) any particle entering
   volumes named "Target" to a file named "particles.mcpl.gz":

    --mcpl=Target

2) Capture generated neutrons to srcneutrons.mcpl.gz:

    --mcpl="grabsrc where trk.is_neutron to srcneutrons.mcpl.gz"

3) Capture high energy gammas entering volumes named "Detector"
   or "Sample":

    --mcpl="Detector,Sample where trk.pdgcode==22 and step.pre.ekin > 1MeV"

4) Capture any particle leaving the "World" volume with double-
   precision storage:

    --mcpl="World,opt_dp,on_vol_exit"

5) Capture any particle entering volumes named "Detector" and embed copy
   number information about the Detector volume and it's mother volume
   into the user flag field of the mcpl file:

    --mcpl="World withflags step.volcopyno+1000*step.volcopyno_1"

Finally, it is of course also possible to configure MCPL writing using Python code inside your sim script, rather than using command-line options, by using the G4MCPLPlugins.MCPLWriter module. See an example here for how to use it.

The MCPL files written by Geant4 contain stat:sum values which describe the simulation as a whole (independently of how many particles ended up in the file), and are added up when files are merged (e.g. with mcpltool --merge, or when merging the output of multiple processes):

stat:sum:geant4_nevts: The number of simulated events.
stat:sum:initial_geant4_weights: The sum of the weights of the primary particles of the simulated events, in the same convention as the weights in the file. For instance, the sum of the weights of the particles in the file divided by this value is an efficiency.
stat:sum:original_initial_ray_count: The stat:sum:initial_ray_count value of the input file read with MCPLGen (if it has one). This is for reference only: the weights in the file are the weights used in Geant4.

The stat:sum:initial_ray_count convention of instrument simulations like McStas (see above) assumes that each particle in a file is an independent ray. This is usually not the case for files written by Geant4, where a single event can produce several (correlated) particles. Therefore Geant4 output files do not claim the convention by default. If it is appropriate for your use case, the opt_irc keyword of --mcpl (or setWriteInitialRayCount(True) on an MCPLWriter) instead records the value of the input file as stat:sum:initial_ray_count, and writes the weights in the same convention as in the input file (i.e. multiplied by the value if MCPLGen divided the weights by it). Reading such a file with MCPLGen again then gives the same weights as in the original simulation.

When the output of multiple processes is merged, each process writes its share of the values describing the whole job (like initial_ray_count), so the merged file gets the value.

Custom stat:sum values describing the whole simulation (for instance the duration of a beam pulse) can be added to all MCPL files written by Geant4 from the sim script, before the simulation is started:

import G4MCPLPlugins.MCPLWriter as Writer
Writer.addStatSum( 'my_beam_time', 3600.0 )

Such values are only added once when the output of multiple processes is merged. C++ code can add values computed at the end of the simulation via the functions in G4MCPL/OutputStatSums.hh, which is also how MCPLGen provides the stat:sum:initial_ray_count of its input file.

Modify MCPL files

The sb_mcplextra_filterfile command can be used to produce a new MCPL file from an existing one, containing only a subset of the original particles. The following example shows how to extract just neutrons with long wavelength (see below for more info on the filter expressions used here):

sb_mcplextra_filterfile infile.mcpl outfile.mcpl "is_neutron && neutron_wl > 0.5Aa"

And here is an example extracting all particles starting within \(1\mathrm{cm}\) of the \(Z\)-axis:

sb_mcplextra_filterfile infile.mcpl outfile.mcpl "sqrt(x^2+y^2)<1cm"

Or, getting the first 1000 particles from a file (perhaps to get a smaller representative file which can be more easily shared for debugging purposes):

sb_mcplextra_filterfile infile.mcpl outfile.mcpl -l1000

Any stat:sum values of the original file (like stat:sum:initial_ray_count from McStas) are kept when filtering particles, since the new file is simply a narrower view of the same simulation. But when only a part of the particles is used due to the limit (-l), the values are marked as unavailable, since the file was then split rather than filtered.

Run sb_mcplextra_filterfile --help for full usage instructions. Note that there also exists a sb_mcplextra_filterview command, which merely prints the selected particles, rather than producing a new file.

Filter expressions

Several of the commands or Python modules above accept flexible filter expressions. This functionality is implemented via the custom ExprParser package inside dgcode, and which variables one can use depends on the context.

Always available:
- Mathematical functions (sin, exp, acos, etc.). See full list in ASTStdMath.cc (look in the create_std_math_function function).
- Physics units (mm, cm, keV, MeV, etc.). See the full list in Units/libinc/Units.hh or the equivalent Python file.
- Mathematical and physics constants (pi, c_light, h_Planck, etc.). Again, see the full list in Units/libinc/Units.hh or the equivalent Python file.
Available when filtering based on content of MCPL files (i.e. when using MCPL file as input):
- Position: x, y, z (aliases: posx, posy, posz)
- Direction: ux, uy, uz (aliases: dirx, diry, dirz)
- Particle type, as PDG code (2112=neutron, 22=gamma, etc.): pdgcode
  - For convenience, the following derived boolean variables are provided: is_neutron, is_gamma (alias is_photon), is_neutrino, is_ion
- Kinetic energy: ekin
  - For convenience, the following derived variable is provided for neutrons, to make it easy to filter on the wavelength: neutron_wl. Note that this variable is only defined for neutrons, so it should be protected in expressions with is_neutron, i.e. "is_neutron && neutron_wl>1Aa”.
- Polarisation: polx, poly, polz
- Time stamp: time
- Weight/intensity: weight
- MCPL userflag value, as integer: userflag
Available when filtering on Geant4 simulations in progress (i.e. when writing MCPL files from Geant4 simulations):
- G4Track-level parameters (refer to Geant4 documentation for details): trk.trkid, trk.parentid, trk.is_primary, trk.is_secondary, trk.creatorprocess, trk.weight, trk.mass, trk.width, trk.charge, trk.lifetime, trk.name, trk.type, trk.subtype, trk.atomicnumber, trk.atomicmass, trk.magneticmoment, trk.spin, trk.is_stable, trk.is_shortlived, trk.pdgcode, trk.is_neutron, trk.is_neutrino, trk.is_opticalphoton, trk.is_photon (optical or standard photon), trk.is_gamma (same as trk.is_photon), trk.is_ion and trk.mcplflag (the MCPL userflag from the input file, in case the G4 event is generated from an MCPL file).
- G4Step-level parameters (again, refer to Geant4 documentation for details): step.edep, step.edep_nonion, step.edep_ion, step.delta_e, step.delta_mom, step.delta_pos, step.delta_time, step.steplength, step.stepnbr, step.volname, step.volname_1 (name of mother volume), step.volname_2 (name of grandmother volume), step.volcopyno, step.volcopyno_1 (copy number of mother volume), step.volcopyno_2 (copy number of grandmother volume), step.mat_name, step.mat_state, step.mat_density, step.mat_temperature, step.mat_pressure, step.process_defined_step.
- G4StepPoint-level parameters (again, refer to Geant4 documentation for details). Must be accessed either as post-step or pre-step, i.e. step.pre.time or step.post.time. Here they are just listed in the step.pre.xxx form: step.pre.stepstatus, step.pre.time, step.pre.globaltime, step.pre.localtime, step.pre.propertime, step.pre.mat_name, step.pre.mat_state, step.pre.mat_density, step.pre.mat_temperature, step.pre.mat_pressure, step.pre.ekin, step.pre.neutron_wl, step.pre.etot, step.pre.velocity, step.pre.beta, step.pre.gamma, step.pre.at_voledge, step.pre.x, step.pre.y, step.pre.z, step.pre.px, step.pre.py, step.pre.pz, step.pre.polx, step.pre.poly, step.pre.polz, step.pre.mom, step.pre.mom2, step.pre.weight, step.pre.process_defined_step, step.pre.volname, step.pre.volname_1, step.pre.volname_2, step.pre.volcopyno, step.pre.volcopyno_1, step.pre.volcopyno_2, step.pre.local_x, step.pre.local_y, step.pre.local_z and step.pre.exists (boolean flag, false if the step-point info is not available - i.e. if step.post.exists is false, one should not access other step.post.xxx variables - this is usually only a concern for the post-step, not the pre-step).