Monte Carlo Particle Lists: MCPL
Thomas Kittelmann, Esben Klinkby, Erik B Knudsen, Peter Willendrup, Xiao Xiao Cai, Kalliopi Kanaki
MMonte Carlo Particle Lists: MCPL
T Kittelmann a, ∗ , E Klinkby b , E B Knudsen c , P Willendrup c,a , X X Cai a,b ,K Kanaki a a European Spallation Source ERIC, Sweden b DTU Nutech, Technical University of Denmark, Denmark c DTU Physics, Technical University of Denmark, Denmark
Abstract
A binary format with lists of particle state information, for interchanging parti-cles between various Monte Carlo simulation applications, is presented. Portable C code for file manipulation is made available to the scientific community, alongwith converters and plugins for several popular simulation packages. PROGRAM SUMMARY
Manuscript Title:
Monte Carlo Particle Lists: MCPL
Authors:
T Kittelmann, E Klinkby, E B Knudsen, P Willendrup, X X Cai andK Kanaki
Program Title:
MCPL
Journal Reference:Catalogue identifier:Licensing provisions:
CC0 for core MCPL, see LICENSE file for details.
Programming language: C and C++
Operating system:
Linux, OSX, Windows
Keywords:
MCPL , Monte Carlo, particles, data storage, simulation, interchange format,
C++ , C , Geant4 , MCNP , McStas , McXtrace
Classification:
4, 9, 11, 17
External routines/libraries:
Geant4 , MCNP , McStas , McXtrace
Nature of problem: ∗ Corresponding author.
Email address: [email protected]
Preprint submitted to Computer Physics Communications March 31, 2017 a r X i v : . [ phy s i c s . c o m p - ph ] M a r ersistification of particle states in Monte Carlo simulations, for interchange betweensimulation packages or for reuse within a single package. Solution method:
Binary interchange format with associated code written in portable C along with toolsand interfaces for relevant simulation packages.
1. Introduction
The usage of Monte Carlo simulations to study the transport and interac-tion of particles and radiation is a powerful and popular technique, finding usethroughout a wide range of fields – including but not limited to both high en-ergy and nuclear physics, as well as space and medical sciences [1]. Naturally,a plethora of different frameworks and applications exist for carrying out thesesimulations (cf. section 3 for examples), with implementations in different lan-guages and domains ranging from general purpose to highly specialised field-and application-specific.A common principle used in the implementation of these applications is therepresentation of particles by a set of state parameters – usually including atleast particle type, time coordinate, position and velocity or momentum vectors– and a suitable representation of the geometry of the problem (either via de-scriptions of actual surfaces and volumes in a virtual three-dimensional space, orthrough suitable parameterisations). In the simplest scenario where no variance-reduction techniques are employed, simulations are typically carried out by pro-ceeding iteratively in steps from an initial set of particles states, with the stateinformation being updated along the way as a result of the pseudo-random ordeterministic modelling of processes affecting the particle. The modelling canrepresent particle self-interactions, interactions with the material of the simu-lated geometry, or simply its forward transport through the geometry, usingeither straight-forward ray-tracing techniques or more complicated trajectorycalculations as appropriate. In addition to a simple update of state parameters,the modelling can result in termination of the simulation for the given particle2r in the creation of new secondary particle states, which will in turn undergosimulation themselves.Occasionally, use-cases arise in which it would be beneficial to be able tocapture a certain subset of particle states present in a given simulation, in orderto continue their simulation at a later point in either the same or a differentframework. Such capabilities have typically been implemented using customapplication-specific means of data exchange, often involving the tedious writingof custom input and output hooks for the specific frameworks and use-cases inquestion. Here is instead presented a standard format for exchange of particlestate data,
Monte Carlo Particle Lists ( MCPL ), which is intended to replace theplethora of custom converters with a more convenient scenario in which expertsof each framework implement converters to the common format, as a one-timeeffort. The idea being that users of the various frameworks then gain the abilityto simply activate those pre-existing and validated converters in order to carryout their work.The present work originated in the needs for simulations at neutron scat-tering facilities, where a multitude of simulation frameworks are typically usedto describe the various components from neutron production to detection, buthistorically other conceptually similar formats have been and are used in highenergy physics to communicate particle states between event generators anddetector simulations [2, 3, 4]. However, these formats were developed for some-what different purposes than the one presented here, keeping simulation histo-ries, focusing on the description of intermediate unphysical or bound particles,existing primarily in-memory rather than on-disk, or implemented in languagesnot readily accessible to applications based on different technologies. For in-stance, [2] is defined as an in-memory
FORTRAN common block, [3] provides a
C++ infrastructure for in-memory data with customisable persistification, and[4] defines a text-based format focused on descriptions of intermediate particlesand lacking particle positions. These existing solutions were thus deemed unfitfor the goals of the work presented here: a compact yet flexible on-disk binaryformat for particle state information, portable, well-defined and able to accom-3odate a wide range of use-cases with close to optimal storage requirements.The accompanying code with which to access and manipulate the files shouldbe small, efficient and easily integrated into existing codes and build systems.Consequently, it was chosen to implement the format through a set of C func-tions declared in a single header file, mcpl.h , and implemented in a single file, mcpl.c . These two files will here be referred to as the core MCPL code, and aremade freely available under the CC0 1.0 Universal Creative Commons license.Along with associated code examples, documentation, configuration files (cf.section 2.5) and application-specific interface code which is not embedded inthe relevant upstream projects (cf. sections 3.1 and 3.2), these files constitutethe
MCPL distribution. The present text concerns the second public release of
MCPL , version 1.1.0. Future updates to the distribution will be made availableat the project website [5].
2. The
MCPL format
MCPL is a binary file format in which a header section, with configurationand meta-data, is followed by a data section, where the state information ofthe contained particles is kept. Data compression is available but optional (cf.section 2.4). The uncompressed storage size of a particle entry in the datasection is determined by overall settings in the header section, and dependson what exact information is stored for the particles in a given file, as willbe discussed shortly. Within a given file, all particle entries will always be ofequal length, allowing for trivial calculation of the absolute data location for aparticle at a given index in the file – and thus for efficient seeking and skippingbetween particles if desired. It is expected and recommended that
MCPL fileswill be manipulated, directly or indirectly, by calls to the functions in mcpl.h (cf. section 2.2), but for reference a complete specification of the binary layoutof data in the files is provided in Appendix A.4 ile header information
Field Description
File type magic number 0x4d43504c ("MCPL") All
MCPL files start with this 4-byte word.Version File format version.Endianness Whether numbers in file are in little- or big-endian format.Number of particles in file 64 bit integer.Flag : Particles have polarisation info If false, all loaded particles will have polarisation vectors (0,0,0).Flag : Particles have user-flags field If false, all loaded particles will have user-flags 0x00000000.Flag : Particle info use double-precision If true, floating point storage use double-precision.Global PDG code If this 32 bit integer is non-zero, all loaded particles will havethis PDG code.Global weight If this double-precision floating point number is non-zero, allloaded particles will have this weight.Source name String indicating the application which created the
MCPL file.Comments A variable number of comments (strings) added at file creation.Binary blobs A variable number of binary data blobs, indexed by keys(strings). This allows arbitrary custom data to be embedded.
Table 1: Information available in the header section of
MCPL files.
The information available in the file header is indicated in Table 1: a unique4-byte magic number identifying the format always starts all files, and is followedby the format version, the endianness ( little or big ) in which numbers in the fileare stored, and the number of particles in the file. The versioning providesa clear path for future updates to the format, without losing the ability toread files created with previous versions of the MCPL code, and the endiannessinformation prevents interpretation errors on different machines (although atpresent, most consumer platforms are little-endian). Next come five optionsindicating what data is stored per-particle, which will be discussed in the nextparagraph. Finally, the header contains several options for embedding customfree-form information: first of all, the source name, in the form of a single stringcontaining the name and perhaps version of the application which created the In the current implementation, reading a little-endian
MCPL file on a big-endian machineor vice versa triggers an error message. It is envisioned that a future version of the
MCPL codecould instead transparently correct the endianness at load time. article state information Field Description Bytes of storage used per entry (FP = 4 or 8 bytes)
PDG code 32 bit integer indicating particle type. 0 or 4Position Vector, values in centimetres. 3FPDirection Unit vector along the particle momentum. 2FPKinetic energy Value in MeV. 1FPTime Value in milliseconds. 1FPWeight Weight or intensity. 0 or 1FPPolarisation Vector. 0 or 3FPUser-flags 32 bit integer with custom info. 0 or 4
Table 2: Particle state information available and uncompressed storage requirements for eachentry in the data section of
MCPL files. file. Secondly, any number of strings can be added as human readable comments,and, thirdly, any number of binary data blobs can be added, each identified by astring key. The
MCPL format itself provides no restrictions on what data, if any,can be stored in these binary blobs, but useful content could for instance be acopy of configuration data used by the source application when the given filewas produced, kept for later reference. Also note that, for reasons of security,no code in the
MCPL distribution ever attempts to interpret contents stored insuch binary data blobs.Table 2 shows the state information available per-particle in
MCPL files, alongwith the storage requirements of each field. Particle position, direction, kineticenergy and time are always stored. Polarisation vectors and so-called user-flags in the form of unsigned 32 bit integers are only stored when relevant flagsin the header are enabled and weights are only stored explicitly in each entrywhen no global common value was set in the header. Likewise, the particle typeinformation in the form of so-called PDG codes is only stored when a global Note that a valid alternative to storing the directional unit vector along with the kineticenergy would have been the momentum vector. However, the choice here is consistent with thevariables used in interfaces of both
MCNP and
Geant4 , and means that the mcpl2ssw converterdiscussed in section 3.2 can be implemented without access to an unwieldy database of particleand isotope masses.
MCPL format is thus designed to beflexible enough to handle use-cases requiring a high level of detail in the parti-cle state information, without imposing excessive storage requirements on lessdemanding scenarios.Note that while the units for position, energy and time indicated in Table 2of course must be respected, the choices themselves are somewhat arbitrary andshould in no way be taken to indicate the suitability of the
MCPL format for agiven simulation task. In particular, note that within the dynamic range of agiven floating point representation, the relative numerical precision is essentiallyindependent of the magnitude of the numbers involved and is determined by thenumber of bits allocated for the significand [7]. Thus, it is important to realisethat usage of the
MCPL format to deal with a simulation task whose natural unitsare many orders of magnitude different than the ones in Table 2 does not implyany detrimental impact on numerical precision.Packing of the three-dimensional unit directional vector into just two floatingpoint numbers of storage is carried out via a new packing algorithm, tentativelynamed
Adaptive Projection Packing , discussed in detail in Appendix B. Unlikeother popular packing strategies considered, the chosen algorithm provides whatis for all practical purposes flawless performance, with a precision comparableto the one existing absent any packing (i.e. direct storage of all coordinates intothree floating point numbers). It does so without suffering from domain validityissues, and the implemented code is not significantly slower to execute than thealternatives. 7 .2. Accessing or creating
MCPL files programmatically
While a complete documentation of the programming API provided by theimplementation of
MCPL in mcpl.h and mcpl.c can be found in Appendix C,the present discussion will restrict itself to a more digestible overview.The main feature provided by the API is naturally the ability to create new MCPL files and access the contents of existing ones, using a set of dedicatedfunctions. No matter which settings were chosen when a given
MCPL file wascreated, the interface for accessing the header and particle state informationwithin it is the same, as can be seen in Listing 1: after obtaining a file handlevia mcpl_open_file , a pointer to an mcpl_particle_t struct , whose fieldscontain the state information available for a given particle, is returned by calling mcpl_read . This also advances the position in the file, and returns a null-pointerwhen there are no more particles in the file, ending the loop. If a file was createdwith either polarisation vectors or user-flags disabled, the corresponding fields onthe particle will contain zeroes (thus representing polarisation information withnull-vectors and user-flags with an integer with no bits enabled). All floatingpoint fields on mcpl_particle_t are represented with a double-precision type,but the actual precision of the numbers will obviously be limited to that storedin the input file. In addition to the interface illustrated by Listing 1, functionscan be found in mcpl.h for accessing any information available in the file header(see Table 1), or for seeking and skipping to particles at specific positions in thefile, rather than simply iterating through the full file.Code creating
MCPL files is typically slightly more involved, as the creationprocess also involves deciding on the values of the various header flags andfilling of free-form information like source name and comments. An exampleproducing a file with 1000 particles is shown in Listing 2. The first part of theprocedure is to obtain a file handle through a call to mcpl_create_outfile ,configure the header and overall flags, and prepare a zero-initialised instance of mcpl_particle_t . Next comes the loop filling the particles into the file, whichhappens by updating the state information on the mcpl_particle_t instanceas needed, and passing it to mcpl_add_particle each time. At the end, a call8 isting 1: Simple example for looping over all particles in an existing
MCPL file. to mcpl_close_outfile finishes up by flushing all internal buffers to disk andupdating the field containing the number of particles at the beginning of thefile.Should the program abort before the call to mcpl_close_outfile , particlesalready written into the output file are normally recoverable: upon openingsuch an incomplete file, the MCPL code detects that the actual size of the file isinconsistent with the value of the field in the header containing the number ofparticles. Thus, it emits a warning message and calculates a more appropriatevalue for the field, ignoring any partially written particle state entry at the endof the file. This ability to transparently correct incomplete files upon load alsomeans that it is possible to inspect (with the mcpltool command discussedin section 2.3) or analyse files that are still being created. To avoid seeing awarning each time a file left over from an aborted job is opened, mcpl.h alsoprovides the function mcpl_repair which can be used to permanently correctthe header of the file. 9 isting 2: Simple example for creating an
MCPL file with 1000 particles. isting 3: Example extracting low-energy neutrons (PDG code 2112) from an MCPL file.
Likewise, mcpl.h also provides the function mcpl_merge_files which canbe used to merge a list of compatible
MCPL files into a new one, which mighttypically be useful when gathering up the output of simulations carried outvia parallel processing techniques. Compatibility here means that the files musthave essentially identical header sections, except for the field holding the numberof particles. Finally, the function mcpl_transfer_metadata can be used toeasily implement custom extraction of particle subsets from existing
MCPL filesinto new (smaller) ones. An example of this is illustrated in Listing 3.
MCPL files from the command line
Compared with simpler text-based formats (e.g. ASCII files with data for-matted in columns), one potential disadvantage of a binary data format like
MCPL is the lack of an easy way for users to quickly inspect a file and inves-tigate its contents. To alleviate this, mcpl.h provides a function which, in astraight-forward manner, can be used to build a generic mcpltool command-line executable: int mcpl_tool(int argc,char** argv) , for which full usageinstructions can be found in Appendix D or by invoking it with the --help flag.11imply running this command on an
MCPL file without specifying other argu-ments, results in a short summary of the file content being printed to standardoutput, which includes a listing of the first 10 contained particles. An exam-ple of such a summary is provided in Listing 4: it is clear from the displayedmeta-data that the particles in the given file represent a transmission spectrumresulting from illumination of a block of lead by a
10 GeV proton beam in a
Geant4 [8, 9] simulation. The displayed header information and data columnsshould be mostly self-explanatory, noting that ( x , y , z ) indicates the particleposition, ( ux , uy , uz ) its normalised direction, and that the pdgcode column in-deed shows particle types typical in a hadronic shower: π + (211), γ (22), protons(2212), π − ( − ) and neutrons (2112). If the file had user-flags or polarisationvectors enabled, appropriate columns for those would be shown as well. Finally,note that the
36 bytes/particle refers to uncompressed storage, and that inthis particular case the file actually has a compression ratio of approximately70%, meaning that about 25 bytes of on-disk storage is used per particle (cf.section 2.4).By providing suitable arguments (cf. Appendix D) to mcpltool , it is pos-sible to modify what information from the file is displayed. This includes thepossibility to change what particles from the file, if any, should be listed, as wellas the option to extract the contents of a given binary data blob to standardoutput. The latter might be particularly handy when entire configuration fileshave been embedded (cf. sections 3.2 and 3.3). Finally, the mcpltool commandalso allows file merging and repairing, as discussed in section 2.2, and providesfunctionality for selecting a subset of particles from a given file and extractingthem into a new smaller file.Advanced functionality such as graphics display and interactive GUI-basedinvestigation or manipulation of the contents of
MCPL files is not provided bythe mcpltool , since those would imply additional unwanted dependencies tothe core
MCPL code, which is required by design to be light-weight and widelyportable. However, it is the hope that the existence of a standard format like
MCPL will encourage development of such tools, and indeed some already exist in12 i s t i n g4 : E x a m p l e o u t pu t o f r unn i n g mcpltool w i t hn oa r g u m e n t s o n a s p ec i fi c MCPL fi l e . OpenedMCPLfilemyoutput.mcpl.gz: BasicinfoFormat:MCPL-3No.ofparticles:1106933Headerstorage:140bytesDatastorage:39849588bytes CustommetadataSource:"G4MCPLWriter[G4MCPLWriter]"Numberofcomments:1->comment0:"Transmissionspectrumfrom10GeVprotonbeamon20cmlead"Numberofblobs:0 ParticledataformatUserflags:noPolarisationinfo:noFixedpart.type:noFixedpart.weight:noFPprecision:singleEndianness:littleStorage:36bytes/particle indexpdgcodeekin[MeV]x[cm]y[cm]z[cm]uxuyuztime[ms]weight0211487.02-0.58981.83520-0.0924070.204910.974417.3346e-0711221.53261.063511.351200.0804410.660260.746721.0882e-0612223.9526-0.439078.647320-0.566160.505580.651041.0286e-0613220.825911.74449.7622200.0920990.795970.598291.0378e-0614221.19582.18068.6416200.219970.664350.714321.0124e-0615221.25253.09497.7366200.489030.307890.816121.013e-0616222.62473.9485.681200.625030.642210.443749.1152e-07172212824.28-1.8797-2.512420-0.3077-0.404960.8617.6539e-0718-2113459.8-0.795210.9148120-0.134410.144380.980357.0618e-071921120.3055354.47133.386200.48620.0119580.873770.000164421
MCPL to includerelevant parts of these tools as a separate and optional component.
The utilisation of data compression in a format like
MCPL is potentially animportant feature, since on-disk storage size could be a concern for some appli-cations. Aiming to maximise flexibility, transparency and portability, optionalcompression of
MCPL files is simply provided by allowing whole-file compressioninto the widespread
GZIP format [13] (changing the file extension from .mcpl to .mcpl.gz in the process). This utilises the DEFLATE compression algorithm [14]which offers a good performance compromise with a reasonable compressionratio and an excellent speed of compression and decompression.Relying on a standard format such as
GZIP means that, if needed, users canavail themselves of existing tools (like the gzip and gunzip commands availableon most
UNIX platforms) to change the compression state of an existing
MCPL file.However, when the code in mcpl.c is linked with the ubiquitous
ZLIB [15, 16](cf. section 2.5), compressed
MCPL files can be read directly. For convenience, mcpl.h additionally provides a function mcpl_closeandgzip_outfile , whichcan be used instead of mcpl_close_outfile (cf. Listing 2) to ensure that newlycreated
MCPL files are automatically compressed if possible (either through acall to an external gzip command or through custom
ZLIB -dependent code,depending on availability).
It is the hope that eventually
MCPL capabilities will be included upstreamin many applications, and that users of those consequently won’t have to doanything extra to start using it. As will be discussed in section 3, this is atpresent the case for users of recent versions of
McStas [17, 18] and
McXtrace [19],and is additionally the case for users of the in-house
Geant4 -based frameworkof the ESS Detector Group [10]. 14y design, it is expected that most developers wishing to add
MCPL supportto their application will simply place copies of mcpl.h and mcpl.c into theirexisting build system and include mcpl.h from either C or C++ code. In order tomake the resulting binary code able to manipulate compressed files directly (cf.section 2.4), the code in mcpl.c must usually be compiled against and linkedwith an installation of
ZLIB (see detailed instructions regarding build flags atthe top of mcpl.c ). Alternatively, the
MCPL distribution presented here containsa “fat” auto-generated drop-in replacement for mcpl.c named mcpl_fat.c , inwhich the source code of
ZLIB has been included in its entirety. Using thissomewhat larger file enables
ZLIB -dependent code in
MCPL even in situationswhere
ZLIB might not be otherwise available.In addition to the core
MCPL code, the
MCPL distribution also contains asmall file providing the mcpltool executable,
C++ files implementing the
Geant4 classes discussed in section 3.1, C files for the mcpl2ssw and ssw2mcpl executa-bles discussed in section 3.2, and a few examples show-casing how user codemight look.Building of all of these parts should be straight-forward using standard tools,but a configuration file for CMake [20] which builds and installs everything isnonetheless provided for reference and convenience. Additionally, “fat” single-fileversions of all command line utilities ( mcpltool , mcpl2ssw and ssw2mcpl ) arealso provided, containing both MCPL and
ZLIB code within as appropriate. Thus,any of these single-file versions can be compiled directly into the correspondingcommand line executable, without any other dependencies than a C compiler.For more details about how to build and deploy, refer to the INSTALL file shippedwith the
MCPL distribution. Compilation of mcpl.c can happen with any of the following standards:
C99 , C11 , C++98 , C++11 , C++14 , or later. In addition to those, mcpl.h is also
C89 compatible. Note that onplatforms where the standard C math function sqrt is provided in a separate library, thatlibrary must be available at link-time. Note that all
ZLIB symbols have been prefixed, to guard against potential run-time clasheswhere a separate
ZLIB is nonetheless loaded. . Application-specific converters and plugins While the examples in section 2.2 show how it is possible to manipulate
MCPL files directly from C or C++ code, it is not envisioned that most users willhave to write such code themselves. Rather, in addition to using available tools(such as the mcpltool described in section 2.3) to access the contents of filesas needed, users would ideally simply use pre-existing plugins and converterswritten by application-specific experts, to load particles from
MCPL files intotheir given Monte Carlo applications, or extract particles from those into
MCPL files. At the time of this initial public release of
MCPL , four such applicationsare already
MCPL -aware in this manner:
Geant4 , MCNP , McStas and
McXtrace ,and the details of the corresponding converters and plugins are discussed in thefollowing sub-sections, after a few general pieces of advice for other implementersin the next paragraphs.In order for
MCPL files to be as widely exchangeable as possible, code loadingparticles from
MCPL files into a given Monte Carlo application should preferablybe as accepting as possible. In particular, this means that warnings ratherthan errors should result if the input file contains PDG codes corresponding toparticle types that can not be handled by the application in question. As anexample, a detailed
MCNP or Geant4 simulation of a moderated neutron sourcewill typically produce files containing not only neutrons, but also gammas andother particles. It should certainly be possible to load such a file into a neutron-only simulation application like
McStas , resulting in simulation of the containedneutrons (preferably with a warning or informative message drawing attentionto some particles being ignored).Applications employing parallel processing techniques, must always pay par-ticular attention when implementing file-based I/O, and this is naturally alsothe case when creating
MCPL -aware plugins for them. However, the availablefunctionality for merging of
MCPL files makes the scenario of file creation partic-ularly simple to implement: each sub-task can simply write its own file, with thesubsequent merging into a single file taking place during post-processing. For16eading of particles in existing
MCPL files, it is recommended that each sub-taskperforms a separate call to mcpl_open_file , and use the skipping and seek-ing functionality to load just a subset of the particles within, as required. Inthe case of a multi-threading application, it is of course also possible to handleconcurrent input or output directly through a single file handle. In this case,however, calls to mcpl_add_particle and mcpl_read must be protected againstconcurrent invocations with a suitable lock or mutex.The following three sub-sections are dedicated to discussions of presentlyavailable
MCPL interfaces for specific Monte Carlo applications. The discussionswill in each case presuppose familiarity with the application in question.
Geant4 interface
In the most typical mode of working with the
Geant4 [8, 9] toolkit, userscreate custom
C++ classes, sub-classing appropriate abstract interfaces, in or-der to set up geometry, particle generation, custom data readout and physicsmodelling. At run-time, those classes are then instantiated and registered withthe framework. Accordingly, the
MCPL – Geant4 integration takes the form of twosuch sub-classes of
Geant4 interface classes, which can be either directly in-stantiated or further sub-classed themselves as needed:
G4MCPLGenerator and
G4MCPLWriter . They are believed to be compatible with any recent version of
Geant4 and were explicitly tested with versions 10.00.p03 and 10.02.p02.First, the
G4MCPLGenerator , the relevant parts of which are shown in List-ing 5, implements a
Geant4 generator by sub-classing the
G4VUserPrimary-GeneratorAction interface class. The constructor of
G4MCPLGenerator mustbe provided with the path to an
MCPL file, which will then be read one particleat a time whenever
Geant4 calls the
GeneratePrimaries method, in order togenerate
Geant4 events with a single primary particle in each. If the file runsout of particles before the
Geant4 simulation is ended for other reasons, the
G4MCPLGenerator graciously requests the
G4RunManager to abort the simula-tion. Thus, a convenient way in which to use the entire input file for simulationis to launch the simulation with a very high number of events requested, as is17 isting 5: The
G4MCPLGenerator class. class G4MCPLGenerator : public G4VUserPrimaryGeneratorAction{ public :G4MCPLGenerator ( const G4String & inputFile );virtual ~ G4MCPLGenerator ();virtual void GeneratePrimaries ( G4Event *);protected :// Reimplement this to filter input particles ( default// implementation accepts all particles ):virtual bool UseParticle ( const mcpl_particle_t *) const ;// Reimplement this to change coordinates or weights of// input particles before using them ( does nothing by// default ):virtual void ModifyParticle ( G4ThreeVector & position ,G4ThreeVector & direction ,G4ThreeVector & polarisation ,G4double & time ,G4double & weight ) const ;private :// ..}; done in the example in Listing 6. In case the user wishes to use only certain particles from the input file for sim-ulation, the
G4MCPLGenerator class must be sub-classed and the
UseParticle method reimplemented, returning false for particles which should be skipped.Likewise, if it is desired to perform coordinate transformations or reweighingbefore using the loaded particles, the
ModifyParticle method must be reim-plemented.The
G4MCPLWriter class, the relevant parts of which are shown in Listing 7,is a
G4VSensitiveDetector which in the default configuration “consumes” allparticles which, during a simulation, enter any geometrical volume(s) to whichit is attached by the user and stores them into the specified
MCPL file. At thesame time it asks
Geant4 to end further simulation of those particles (“killing” Unfortunately, due to a limitation in the
G4RunManager interface, this number will belimited by the highest number representable with a
G4int , which on most modern platformsis 2147483647. isting 6: Example showing how to load particles from an MCPL file into a
Geant4 simulation. them). This strategy of killing particles stored into the file was chosen as a sen-sible default behaviour, as it prevents potential double-counting in the scenarioswhere a particle (or its induced secondary particles) would otherwise be able toenter a given volume multiple times. If it is desired to modify this strategy, theuser must sub-class
G4MCPLWriter and reimplement the
ProcessHits method,using calls to
StorePreStep , StorePostStep and
Kill , as appropriate. Forreference, code responsible for the default implementation is shown in Listing 8.Likewise, to add
MCPL user-flags into the file, the
UserFlagsDescription and
UserFlags methods must simply be reimplemented - the description naturallyending up as a comment in the output file.In Listing 9 is shown how the
G4MCPLWriter will typically be configured andattached to logical volume(s) of the geometry.
MCNP interface
Most users of
MCNP are currently employing one of three distinct flavours:
MCNPX [21, 22],
MCNP5 [23] or
MCNP6 [24]. In the most typical mode of workingwith any of these software packages, users edit and launch
MCNP through the useof text-based configuration files (so-called input decks ), in order to set up detailsof the simulation including geometry, particle generation, and data extraction.The latter typically results in the creation of data files containing simulation19 isting 7: The
G4MCPLWriter class. class G4MCPLWriter : public G4VSensitiveDetector{ public :// Basic interface :G4MCPLWriter ( const G4String & outputFile ,const G4String & name = " G4MCPLWriter " );virtual ~ G4MCPLWriter ();void AddComment ( const G4String & );void AddData ( const G4String & data_key ,size_t data_length ,const char * data );void EnableDoublePrecision ();void EnablePolarisation ();void EnableUniversalWeight ( G4double );// Optional reimplement this to change default//" store -and - kill at entry " strategy :virtual G4bool ProcessHits ( G4Step * step ,G4TouchableHistory * );// Optional reimplement these to add MCPL userflags :virtual G4String UserFlagsDescription () const { return ""; }virtual uint32_t UserFlags ( const G4Step *) const { return 0; }protected :// Methods that can be used if reimplementing ProcessHits ():void StorePreStep ( const G4Step *);void StorePostStep ( const G4Step *);void Kill ( G4Step *);private :// ...}; isting 8: The default ProcessHits implementation in the
G4MCPLWriter class.
G4bool G4MCPLWriter :: ProcessHits ( G4Step * step , G4TouchableHistory *){ // Only consider particle steps originating at the boundary// of the monitored volume :if ( step -> GetPreStepPoint ()-> GetStepStatus () != fGeomBoundary )return false ;// Store the state at the beginning of the step , but avoid// particles taking their very first step ( this would double -// count secondary particles created at the volume edge ):if ( step -> GetTrack ()-> GetCurrentStepNumber () > 1 )StorePreStep ( step );// Tell Geant4 to stop further tracking of the particle :Kill ( step );return true ;}
Listing 9: Example showing how to produce an
MCPL file from a
Geant4 simulation. // Provide output filename when creating G4MCPLWriter instance :G4MCPLWriter * mcplwriter = new G4MCPLWriter (" myoutput . mcpl ");// Optional calls which add meta - data or change flags :mcplwriter -> AddComment (" Some useful description here ");mcplwriter -> AddData ( ... );mcplwriter -> EnableDoublePrecision ();mcplwriter -> EnablePolarisation ();mcplwriter -> EnableUniversalWeight (1.0);// Register with G4SDManager and on one or more logical// volumes to activate :G4SDManager :: GetSDMpointer ()-> AddNewDetector ( mcplwriter );alogvol -> SetSensitiveDetector ( mcplwriter );
FORTRAN -compatible
MCPL hooks for
MCNP , such an approach would require users to undertake someform of compilation and linking procedure. This would likely impose a change inworking mode for the majority of
MCNP users, in addition to possibly requiring aspecial license for source-level access to
MCNP . Instead, the
MCNP – MCPL interfacepresented here exploits the existing
MCNP capability to stop and subsequentlyrestart simulations at a user-defined set of surfaces through the
Surface SourceWrite/Read ( SSW / SSR ) functionality. As the name suggests, the state parame-ters of simulated particles crossing a given surface are stored on disk in dedicatedfiles, with the intentional use as a surface source in subsequent simulations withthe same
MCNP setup. Presumably, these files (henceforth denoted “
SSW files” inthe present text) are intended for this internal intermediate usage only, sincetheir format differs between different flavours of
MCNP , and little effort has beenmade to document the format in publicly available manuals. Despite these ob-stacles, the
SSW format is stable enough that several existing
MCNP -aware tools(e.g. [25, 26, 27]) have chosen to provide converters for this format, with variouslevels of functionality, and it was thus deemed suitable also for the needs of the
MCPL project.Thus, the
MCPL distribution presented here includes dependency-free C codefor two standalone executables, mcpl2ssw and ssw2mcpl , which users can invokefrom the command-line in order to convert between MCPL and
SSW files. Theusage of these two executables will be discussed here, while users are referred tothe relevant
MCNP manuals for details of how to set up their input decks to enable
SSW input or output in their
MCNP simulations: [28, Ch. II.3.7], [29, Ch. 5.5.5]and [30, Ch. 3.3.4.7]. Note that through usage of ssw2mcpl and mcpl2ssw , itis even possible to transfer particles between different flavours and versions of
MCNP , which is otherwise not possible with
SSW files. Prior work in [25, 26] served as valuable input when developing code for interpreting datasections in
SSW files. isting 10: Usage instructions for the ssw2mcpl command. Usage:ssw2mcpl [ options ] input.ssw [ output .mcpl]Converts the Monte Carlo particles in the input.ssw file (MCNP SurfaceSource Write format ) to MCPL format and stores in the designated outputfile ( defaults to " output .mcpl ").Options :-h, --help : Show this usage information .-d, --double : Enable double - precision storage of floating point values .-s, --surf : Store SSW surface IDs in the MCPL userflags .-n, --nogzip : Do not attempt to gzip output file.-c FILE : Embed entire configuration FILE (the input deck)used to produce input.ssw in the MCPL header .
First, the ssw2mcpl command, for which the full usage instructions areshown in Listing 10, is in its most base invocation straight-forward to use.Simply provide it with the name of an existing
SSW file to run on, and it willresult in the creation of a new (compressed)
MCPL file, output.mcpl.gz , con-taining a copy of all particles found in the
SSW file. The
MCNP flavour responsiblefor creating the
SSW file is automatically detected, the resulting differences inthe file format are taken into account behind the scenes, and the detected
MCNP version is documented as a comment in the header of the resulting
MCPL file.The only relevant piece of information which is by default not transferredfrom the
SSW particle state into the
MCPL file is the numerical ID of the surfacewhere the particle was registered in the
MCNP simulation. By supplying the -s option, ssw2mcpl will transfer those to the MCPL user-flags field, and documentthis in the
MCPL header. Additionally, while floating point numbers in the
SSW file are always stored in double-precision, the transfer to
MCPL will by defaultconvert them to single-precision. This was chosen as the default behaviour tokeep usual storage requirements low, as single-precision is arguably sufficientfor most studies. By supplying the -d option, ssw2mcpl will keep the numbersin double-precision in the MCPL file as well. Depending on compression andthe applied flags, the on-disk size of the resulting
MCPL file will typically besomewhere between 20% and 80% of the on-disk size of the
SSW file from whichit was converted. 23inally it is possible, via the -c FILE flag, to point the ssw2mcpl commandto the input deck file used when producing the provided
SSW file. Doing so willresult in a complete copy of that file being stored in the
MCPL header as a binarydata blob under the string key "mcnp_input_deck" , thus providing users witha convenient snapshot in the
MCPL file of the
MCNP setup used. Unfortunately,it was not possible to automate this procedure completely, and it thus relieson the user to provide the correct input deck for a given
SSW file. But the ssw2mcpl command does at least check that the specified file is a text-file andthat it contains somewhere the correct value of the so-called problem title : acustom free-form string which is specified by the user in the input deck andembedded in the
SSW file by
MCNP . The input deck embedded in a given
MCPL file can later be inspected from the command line by invoking the command“ mcpltool -bmcnp_input_deck
MCPL file, the user must also supply a reference
SSW file in a format suitable for the
MCNP setup in which the resulting
SSW file is subsequently intended to be used asinput. The need for this added complexity stems from the constraint that the
SSW format is merely intended as an internal format in which it is possible to stopand restart particles while remaining within a given setup of an
MCNP simulation– meaning at the very least that the
MCNP version and the configuration of thegeometrical surfaces involved in the
Surface Source Write/Read proceduremust be unchanged. Thus, for maximal robustness, the user must supply areference
SSW file which was produced by the setup in which the
SSW file createdwith mcpl2ssw is to be used (it does not matter how many particles the referencefile contains). What will actually happen is that in addition to the particle statedata itself, the newly created
SSW file will contain the exact same header as theone in the reference
SSW file, apart from the fields related to the number ofparticles in the file.Additionally, the user must consider carefully which
MCNP surface IDs theparticles from the
MCPL file should be associated with, once transferred to the24 isting 11: Usage instructions for the mcpl2ssw command.
Usage:mcpl2ssw [ options ]
SSW file. By default it will assume that the
MCPL user-flags field contains exactlythis ID, but more often than not, users will have to specify a global surface ID forall of the particles through the -s
SSW files do not contain polarisation information, andany such polarisation information in the input
MCPL file will consequently bediscarded in the translation. Likewise, in cases where the input
MCPL file containsone or more particles whose type does not have a representation in the targetedflavour of
MCNP , they will be ignored with suitable warnings.
McStas and
McXtrace interfaces
Recent releases of the neutron ray tracing software package
McStas [17, 18](version 2.3 and later) and its X-ray sibling package
McXtrace [19] (version 1.4and later) include
MCPL -interfaces. Although
McStas and
McXtrace are twodistinct software packages, they are implemented upon a common technologi-cal platform,
McCode , and the discussions here will for simplicity use the term25 isting 12: Code enabling
MCPL input in its simplest form.
COMPONENT vin = MCPL_input ( filename =" myfile . mcpl " )AT (0 ,0 ,0) RELATIVE Origin
McCode where the instructions are otherwise identical for users of the two pack-ages.The particle model adopted in
McCode is directly compatible with
MCPL . Inessence, apart from simple unit conversions, particles are read from or writtento
MCPL files at one or more predefined logical points defined in the
McCode con-figuration files (so-called instrument files ). Specifically, two new components,
MCPL_input and
MCPL_output , are provided, which users can activate by addingentries at relevant points in their instrument files as is usual when working with
McCode .First, when using the
MCPL_input component, particles are directly readfrom an
MCPL input file and injected into the simulation at the desired point,thus playing the role of a source. In Listing 12 is shown how, in its simplestform, users would insert an
MCPL_input component in their instrument file.This will result in the
MCPL file being read in its entirety, and all found neutrons(for
McStas ) or gamma particles (for
McXtrace ) traced through the
McCode simulation. Listing 13 indicates how the user can additionally impose an al-lowed energy range when loading particles by supplying the
Emin and
Emax parameters. The units are meV and keV respectively for
McStas and
McXtrace .Thus, the code in Listing 13 would select –
100 meV neutrons in
McStas and –
100 keV gammas in
McXtrace . A particle from the
MCPL file is injected atthe position indicated by its
MCPL coordinates relative to the position of the
MCPL_input component in the
McCode instrument. Thus, a user can imposecoordinate transformations by altering the positioning of
MCPL_input as shownin Listing 14, which would shift the initial position of the particles by ( X, Y, Z ) and rotate their initial velocities around the x , y and z axes (in that order) byrespectively Rx , Ry and Rz degrees. Furthermore, Listing 14 shows a way tointroduce a time shift of to all particles, using an EXTEND code block.26 isting 13: Code enabling
MCPL input, selecting particles in a given energy range.
COMPONENT vin = MCPL_input ( filename =" myfile . mcpl ",Emin =12 , Emax =100 )AT (0 ,0 ,0) RELATIVE Origin
Listing 14: Code enabling
MCPL input, applying spatial and temporal transformations.
COMPONENT vin = MCPL_input ( filename =" myfile . mcpl " )AT(X,Y,Z) RELATIVE OriginROTATED (Rx ,Ry ,Rz) RELATIVE OriginEXTEND%{ t=t +2;%}
For technical reasons, the number of particles to be simulated in
McCode must be fixed at initialisation time. Thus, the number of particles will be set tothe total number of particles in the input file, as this is provided through thecorresponding
MCPL header field. If and when a particle is encountered whichcan not be used (due to having a wrong particle type or energy), it will leadto an empty event in which no particles leave the source. At the end of therun, the number of particles skipped over will be summarised for the user. Thisapproach obviates the need for running twice over the input file and avoids thepotential introduction of statistical bias from reading a partial file.Note that if running
McCode in parallel processing mode using MPI [31],each process will operate on all particles in the entire file, but the particleswill get their statistical weights reduced accordingly upon load. This behaviouris not specific to the
MCPL_input component, but is a general feature of howmultiprocessing is implemented in
McCode .When adding an
MCPL_output component to a
McCode instrument file, thecomplete state of all particles reaching that component is written to the re-quested output file. In Listing 15 is shown how, in its simplest form, userswould insert such a component in their instrument file, and get particles writ-ten with coordinates relative to the component preceding it, into the out-put file (replace
RELATIVE PREVIOUS with
RELATIVE ABSOLUTE to write ab-27 isting 15: Code enabling
MCPL output in its simplest form.
COMPONENT mcplout = MCPL_output ( filename =" myoutput . mcpl " )AT (0 ,0 ,0) RELATIVE PREVIOUS
Listing 16: Code enabling
MCPL output with polarisation and double-precision numbers.
COMPONENT mcplout = MCPL_output ( filename =" myoutput . mcpl ",polarisationuse =1 ,doubleprec =1 )AT (0 ,0 ,0) RELATIVE PREVIOUS solute coordinates instead). For reference, a copy of the complete instrumentfile is stored in the
MCPL header as a binary data blob under the string key "mccode_instr_file" . This feature provides users with a convenient snapshotof the generating setup. The instrument file embedded in a given
MCPL filecan be inspected from the command line by invoking the command “ mcpltool-bmccode_instr_file
McCode in parallel processing mode using MPI, each process willcreate a separate output file named after the pattern myoutput.node_idx.mcpl where idx is the process number (assuming filename="myoutput.mcpl" as inListing 15), and those files will be automatically merged during post-processinginto a single file. To avoid generating unnecessarily large files, the
MCPL_output componentstores particle state data using the global PDG code feature (cf. section 2.1),uses single-precision floating point numbers, and does not by default store po-larisation vectors. The two latter settings may be changed by the user throughthe polarisationuse and doubleprec parameters respectively, as shown inlisting 16.Finally, if desired, custom information might be stored per-particle into the
MCPL user-flags field for later reference. This could be any property, such as This automatic merging only happens in
McStas version 2.4 or later and
McXtrace version1.4 or later, and can be disabled by setting the parameter merge_mpi=0 . Users of earlierversions will have to use the mcpltool command to perform the merging manually, if desired. isting 17: Code enabling MCPL output with custom user-flags information. /* some upstream component setting a variable ( customvar ) */COMPONENT some_comp = Some_Component ( /* some parameters here */ )AT (0 ,0 ,0) ABSOLUTEEXTEND %{customvar =( uint32_t ) mcget_run_num ();%}/* MCPL output capturing customvar into MCPL user - flags */COMPONENT vout = MCPL_output ( filename =" myoutput . mcpl ",userflag = customvar ,userflagcomment =" Particle Id" )AT (0 ,0 ,0) RELATIVE PREVIOUS for instance the number of reflections along a neutron guide, or the type ofscattering process in a crystal, etc. Listing 17 shows a simple example of thiswhere the particle ID, in the form of its
McCode ray number (returned from the
McCode library function mcget_run_num ), is stored into the user-flags field. Astring, userflagcomment , is required in order to describe the significance of theextra data, and will end up as a comment in the resulting
MCPL file.
4. Example scientific use cases
The possible uses for
MCPL are envisioned to be many and varied, facilitatingboth straight-forward transfers of particle data between different simulations,as well as data reuse and cross-code comparisons. Actual scientific studies arealready being performed with the help of
MCPL , demonstrating the suitability ofthe format “in the field”. By way of example, it will be discussed in the followinghow
MCPL is used in two such ongoing studies.
The ongoing construction of the European Spallation Source (ESS) [11, 12]has initiated significant development of novel neutronic technologies in the past5 years. The performance requirements for neutron instruments at the ESS,in particular those resulting from the unprecedented cold and thermal neutronbrightness, are at or beyond the capabilities of detector technologies currently29vailable [32]. Additionally, shortage of He [33, 34], upon which the vast major-ity of previous detectors were based, augments the need for development of newefficient and cost-effective detectors based on other isotopes with high neutronicconversion cross sections.A typical approach to instrument design and optimisation at ESS involvesthe development of a
McStas -based simulation of the instrument. Such a sim-ulation includes an appropriate neutron source description and detailed modelsof the major instrument components, such as benders, neutron guides, choppersystems, collimators, sample environment and sample. See [35] for an introduc-tion to the role of the various instrument components. Detector components in
McStas are, however, typically not implemented with any detailed modelling,and are simply registering all neutrons as they arrive. Thus, while the setupin
McStas allows for an efficient and precise optimisation of most of the instru-ment parameters, detailed detector optimisation studies must out of necessitybe carried out in a separate simulation package, such as
Geant4 .As the detector development progresses in parallel with the general instru-ment design, it is crucial to be able to optimise the detector setup for the exactinstrument conditions under investigation in
McStas . The
MCPL format, alongwith the interfaces discussed in sections 3.1 and 3.3, facilitates this by allowingfor easy transfer of neutron states from the
McStas instrument simulation into
Geant4 simulations with detailed setups of proposed detector designs.Technically, this is done by placing the
MCPL_output component just afterthe relevant sample component in the
McStas instrument file. Additionally,using the procedure for creation and storage of custom
MCPL user-flags also dis-cussed in section 3.3, it is possible to differentiate neutrons that scattered onthe sample from those which continued undisturbed, and to carry this informa-tion into the
Geant4 simulations. This information is needed to understand theimpact of the direct beam on the low angle measurements, in order to study therequirements for a so-called zero-angle detector.For example, in order to optimise the detector technology that the LoKIinstrument [36, 37, 38] might adopt, a series of
McStas simulations of the in-30 igure 1: Layout of the
McStas model of the LoKI instrument. Neutrons originate at thesource located at z = 0 and progress through the various instrument components toward thesample at z = 22 . . strument components and the interactions in realistic samples [39] are performed(see Figure 1 for a view of the instrument in McStas ). The parameters of theinstrument and the samples in the
McStas model are chosen in such a way,that various aspects of the detector performance can be investigated, includingrate capability and spatial resolution. The neutrons emerging from the samplein
McStas are then transferred via
MCPL to the detector simulation in
Geant4 ,where a detailed detector geometry and appropriate materials are implemented(see Figure 2 for a visualisation of the
Geant4 model).Neutrons traversing the detector geometry in
Geant4 undergo interactionswith the materials they pass on their flight-path, according to the physics pro-cesses and respective cross sections available in the setup. Special attentionis needed when configuring the
Geant4 physics modelling, to ensure that allprocesses relevant for neutron detection are taken into account and handled31 igure 2:
Geant4 model of a potential detector geometry for the LoKI instrument. Neutronsfrom the sample hitting the active detector area appear in red. correctly. Specifically, the setup utilises the high-precision neutron models in
Geant4 extended with [40], and is implemented in [10]. In the solid-converterbased detectors under consideration, a neutron absorption results in emissionof charged products which then travel a certain range inside the detector anddeposit energy in a counting gas. It is possible to extract position and timeinformation from the energy deposition profile and use these space-time coordi-nates for further analysis, in the same way that measurements in a real detectorwould be treated. This way it becomes possible to reproduce the distributionsof observable quantities relevant for Small Angle Neutron Scattering (SANS)analysis [41, 42].One such observable quantity is the Q distribution [35, Ch. 2.3.3], where Q is defined as the momentum change of the neutron as it scatters on the sample,divided by ¯ h : Q ≡ | ∆ (cid:126)p | / ¯ h . Figure 3 demonstrates such a distribution, based onthe simulated output of the middle detector bank of LoKI (cf. Figure 2), for acertain instrument setup – including a sample modelled as consisting of sphereswith radii of Å. The raw Q distribution is calculated both based on the neu-tron states as they emerge from the sample in McStas , and from the simulatedmeasurements in
Geant4 . With such a procedure, resolution-smearing effectscan be correctly attributed to their sources, geometrical acceptance and detec-32 igure 3: Raw Q distribution for a subset of the LoKI detectors (middle detector bank ofFigure 2). The McStas post-sample output appears in blue, while the distribution calculatedfrom the simulated measurements in
Geant4 appears in red. tor efficiency can be studied in detail, and the impact of engineering featuressuch as dead space can be accurately considered.
The use of radionuclides produced in-situ by cosmic rays for dating purposeshas, in the last two decades, revolutionised the earth surface sciences [43]. Theprecise determination of the production rate of such isotopes, like Be and Al, poses the key challenge for this technique and relies on a folding of cos-mic fluxes with energy dependent production cross sections [44]. The presentdiscussion will focus on the evaluation of the neutron flux induced by cosmicradiation, and in particular on how
MCPL can be exploited both to facilitate thereuse of computationally intensive simulations, and as a means for cross-codecomparisons.At sea level, neutrons constitute the most abundant hadronic componentof cosmic ray induced showers, and possess relatively high cross sections for33roduction of isotopes relevant for radionuclide dating. Thus, it is the dominantcontributor to the relevant isotopic production in the first few meters below thesurface [45]. Extending further below the surface, the neutron flux decreasesrapidly, and as a consequence the isotopic production rate induced by cosmicmuons eventually becomes the most significant factor [46, 47]. At a depthof approximately below the surface, the production rate due to muons iscomparable with the rate from neutrons [45]. Considering non-erosive surfacesand samples at depths significantly less than , the production rates can thusbe estimated by considering just the flux of neutrons. Thus, given known crosssections for neutronic production of Be or Al, properties such as the cosmicirradiation time of a given sample can be directly inferred from its isotopiccontent – providing information about geological activity. In the present study,Monte Carlo methods are used to simulate atmospheric cosmic rays [48, 49] andsubsequently estimate the neutron flux spectra as a function of depth under thesurface of the Earth.Primary cosmic rays constantly bombard the solar system and initiate cosmicray showers in the Earth’s atmosphere, leading to the production of atmosphericneutrons. Figure 4 shows the trajectories of a simulated air shower induced bya single
100 GeV proton in
Geant4 : very large numbers of secondary particlesare generated in each shower, all of which must themselves undergo simulation.Full scale simulation of such showers is therefore relatively time consuming. Onthe other hand, simulations of the propagation of sea level neutrons in a fewmeters of solid material are relatively fast. In the present work of estimatingneutron spectra for different underground materials,
MCPL is used to recordparticle information at sea level. Using the recorded data as input, subsequentsimulations are dedicated to the neutron transport in different undergroundmaterials. In this way, repetition of the time consuming parts of the simulationis avoided.
Geant4 is used to simulate the air shower in this work, while both
Geant4 and
MCNPX are used to simulate neutron spectra underground.In the
Geant4 simulation of the Earth’s atmosphere, the geometry is imple-mented as a
100 km thick shell with an inner radius of , sub-divided34 igure 4: Cosmic shower simulated in
Geant4 . The incident proton energy is
100 GeV andthe length of the x -axis is . The straight grey trajectories are neutrinos. The yellow andgreen trajectories are photons and neutrons, respectively. into 50 equally thick layers, the effective temperatures and densities of whichare calculated using the “U.S. standard atmosphere, 1976” model [50]. Usingthe plugins described in section 3.1, the simulation of any particle reaching theinner surface of the atmosphere is ended and its state stored in an MCPL file.To compare the simulated and measured [51] spectra at New York city, a lowercutoff of E c = 2 .
08 GeV on the kinetic energy of the primary proton is applied,to take the geomagnetic field shielding effect at this location into account. Therelationship between the number of simulated primary protons, N , and the realworld time-span, δt , to which such a sample-size corresponds, is given by thefollowing equation: δt = N ∞ (cid:82) E c J ( E ) dE × π × πr Here, r is the outer radius of the simulated atmosphere and J the differentialspectrum of Usoskin’s model [52] using the parameterisation in [53].In the simulation of . × p rimary protons, the resulting integral neu-35ronic flux above
20 MeV at sea level was found to be . × − cm − , corre-sponding to an absolute surface flux at New York city of . × − cm − s − .Integrating the measured reference neutron spectrum tabulated in [51] above
20 MeV , an integral flux of . × − cm − s − is obtained. The simulationthus overestimates the measured flux by 34%, which is a level of disagreementcompatible with the variation in the predicted value of the integral flux be-tween different models of the local interstellar spectrum [53]. Therefore, theperformance of the atmospheric simulation is concluded to be satisfactory.In the subsequent underground simulations presented here, the Earth is forsimplicity modelled as consisting entirely of quartz (SiO ), which is a samplematerial widely used in cosmogenic dating applications [43], as both Al and Be are produced within when subjected to neutron radiation – normally viaspallation. The
MCPL files generated by the computationally expensive atmo-spheric shower simulation described above, is input to the underground simu-lations implemented in both
Geant4 and
MCNPX , using the interfaces describedin sections 3.1 and 3.2. The geometries in both cases are defined as
20 cm thickspherical shells consisting of pure quartz. As the threshold energies of the re-lated spallation reactions are well above
20 MeV , only spectra above this energyare compared in this study. The simulated volume spectra in a few layers arecompared in Figure 5. Good agreement between
Geant4 and
MCNPX is observed.In conclusion, a useful method for disentangling the resource intensive simu-lation of cosmic showers from subsequent faster simulations of neutron transportin the Earth crust has been demonstrated using
MCPL as an intermediate steppingstone. The simulation strategy thus employed eases the use of computationalresources, and provides a means for cross-comparison between simulation codes.Given reliable energy dependent cross sections, many of the key parameters forcosmogenic dating applications can be provided based on the work described inthis section. 36 energy, MeV10 -32 -31 -30 -29 -28 -27 -26 -25 -24 -23 d i ff e r e n t i a l f l u x , c m − s − M e V − Figure 5: Comparisons of simulated neutron spectra in underground quartz.
5. Summary and outlook
The
MCPL format provides flexible yet efficient storage of particle-state in-formation, aimed at simplifying and standardising interchange of such data be-tween applications and processes. The core parts of
MCPL are implemented inportable and legally unencumbered C code. This is intended to facilitate adop-tion into existing packages and build systems, and the creation of application-specific converters and plugins.In connection with the initial release presented here, MCPL interfaces werecreated for several popular Monte Carlo particle simulation packages:
Geant4 , MCNP , McStas and
McXtrace . It is the intention and hope that the number ofsuch
MCPL -aware applications will increase going forward. A website [5] has beenset up for the
MCPL project, on which users will be able to locate future updatesto the
MCPL distribution, as well as relevant documentation.37 cknowledgements
This work was supported in part by the European Union’s Horizon 2020 re-search and innovation programme under grant agreement No 676548 (the Bright-nESS project) and under grant agreement No 654000 (the SINE2020 project).
ReferencesReferences [1] D. A. Kling, et al., Advanced Monte Carlo for Radiation Physics, ParticleTransport Simulation and Applications, Proceedings of the Monte Carlo2000 Conference, Lisbon, 23-26 October 2000, Springer, 2000.[2] T. Sjöstrand, et al., in: G. Altarelli, et al. (Eds.), Z physics at LEP 1,volume 3, CERN, Geneva, 1989, pp. 327–330. doi: .[3] M. Dobbs, J. B. Hansen, Computer Physics Communications 134 (2001)41 – 46. doi: .[4] J. Alwall, et al., Computer Physics Communications 176 (2007) 300 – 304.doi: .[5] T Kittelmann and others, Website of MCPL: Monte Carlo Particle Lists,2016. URL: https://mctools.github.io/mcpl/ .[6] K A Olive, et al. (Particle Data Group), Chinese Physics C 38 (2014)090001. doi: .[7] IEEE Task P754, IEEE 754-2008, Standard for Floating-Point Arithmetic,pub-IEEE-STD, 2008. doi: .[8] S. Agostinelli, et al. (GEANT4), Nucl. Instrum. Meth. A506 (2003) 250–303. doi: .389] J. Allison, et al., IEEE Trans. Nucl. Sci. 53 (2006) 270. doi: .[10] T. Kittelmann, et al., J. Phys: Conf. Ser. 513 (2014) 022017. doi: .[11] S. Peggs, et al., ESS Conceptual Design Report, ESS 2012-001, http://europeanspallationsource.se/scientific-technical-documentation , 2012.[12] S. Peggs, et al., ESS Technical Design Report, ESS 2013-001, http://europeanspallationsource.se/scientific-technical-documentation , 2013.[13] L. P. Deutsch, GZIP file format specification version 4.3, RFC 1952, IETF,1996. doi: .[14] L. P. Deutsch, DEFLATE Compressed Data Format Specification version1.3, RFC 1951, IETF, 1996. doi: .[15] Jean-Loup Gailly and Mark Adler, ZLIB Compression library, 2013. URL: http://zlib.net/ , v1.2.8.[16] L. P. Deutsch, J.-L. Gailly, ZLIB Compressed Data Format Specificationversion 3.3, RFC 1950, IETF, 1996. doi: .[17] K. Lefmann, K. Nielsen, Neutron News 10 (1999) 20–23. doi: .[18] P. Willendrup, E. Farhi, K. Lefmann, Physica B: Condensed Matter 350(2004) E735 – E737. doi: .[19] E. Bergbäck Knudsen, et al., Journal of Applied Crystallography 46 (2013)679–696. doi: .[20] K. Martin, B. Hoffman, Mastering CMake: A Cross-Platform Build System,Kitware Inc, 2015. 3921] L. S. Waters, et al., AIP Conf. Proc. 896 (2007) 81–90. doi: .[22] D. B. Pelowitz, et al., MCNPX 2.7.0 Extensions, Technical Report LA-UR-11-02295, Los Alamos National Laboratory, 2011.[23] B. K. Forrest Brown, J. Bull, MCNP5-1.60 Release Notes, Technical ReportLA-UR-10-06235, Los Alamos National Laboratory, 2010.[24] J. Goorley, et al., Initial MCNP6 Release Overview - MCNP6 version 1.0,Technical Report LA-UR-13-22934, Los Alamos National Laboratory, 2013.[25] E. Klinkby, et al., Nucl. Instrum. Meth. A700 (2013) 106 – 110. doi: .[26] K. Batkov, et al., mctools - Some Monte Carlo tools for MCNP(X), PHITSand FLUKA, https://github.com/kbat/mc-tools , 2015.[27] The PyNE Development Team, PyNE: The Nuclear Engineering Toolkit,2014. URL: http://pyne.io .[28] X-5 Monte Carlo Team, MCNP - A General N-Particle Transport Code,Version 5, Technical Report, Los Alamos National Laboratory, 2003. Vol-ume II: User’s Guide, LA-CP-03-0245.[29] D. B. Pelowitz, MCNPX User’s Manual, Version 2.6.0, Technical ReportLA-CP-07-1473, Los Alamos National Laboratory, 2008.[30] D. B. Pelowitz, et al., MCNP6 Users Manual, Version 1.0, Technical ReportLA-CP-13-00634, Los Alamos National Laboratory, 2013.[31] M. P. I. Forum, MPI: A Message-Passing Interface Standard, Version 3.1,High Performance Computing Center Stuttgart (HLRS), 2015.[32] O. Kirstein, et al., Proceedings of Science (Vertex 2014) (2014) 029.doi: arXiv:1411.6194 .[33] A. Cho, Science 326 (2009) 778–779. doi: .4034] T. M. Persons, G. Aloise, Technology Assessment: Neutron Detectors: Al-ternatives to Using Helium-3, , Technical Report, U.S. General AccountingOffice, 2011. GAO-11-753.[35] B. Willis, C. Carlile, Experimental Neutron Scattering, Oxford UniversityPress, 2009.[36] A. Jackson, et al., Proceedings ICANS XXI (2015). doi: .[37] K. Kanaki, et al., Journal of Applied Crystallography 46 (2013) 1031–1037.doi: doi:10.1107/S0021889813011862 .[38] K. Kanaki, et al., Journal of Applied Crystallography 46 (2013) 1528.doi: doi:10.1107/S0021889813022826 .[39] SasView, SasView Model Functions, 2014. URL: , v3.0.0.[40] T. Kittelmann, M. Boin, Computer Physics Communications 189 (2015)114 – 118. doi: .[41] L. Feigin, D. Svergun, Structure Analysis by Small-Angle X-Ray and Neu-tron Scattering, Springer, 1987.[42] T. Imae, et al., Neutrons in Soft Matter, Wiley, 2011.[43] T. J. Dunai, Cosmogenic Nuclides: Principles, Concepts and Applicationsin the Earth Surface Sciences, Cambridge University Press, 2010.[44] R. Reedy, Nucl. Instrum. Meth. B 293 (2013) 470 – 474. doi: doi:10.1016/j.nimb.2011.08.034 .[45] J. Gosse, F. Phillips, Quat. Sci. Rev. 20 (2001) 1475 – 1560.[46] B. Heisinger, et al., Earth Planet. Sci. Lett. 200 (2002) 345 – 355.[47] B. Heisinger, et al., Earth Planet. Sci. Lett. 200 (2002) 357 – 369.4148] J. Masarik, J. Beer, J. Geophys. Res. 104 (1999) 12099 – 12111.[49] J. Masarik, J. Beer, J. Geophys. Res. 114 (2009) D11103. doi: .[50] A. Jursa (Ed.), Handbook of Geophysics and the Space Environment, 1985.[51] M. Gordon, et al., IEEE Trans. Nucl. Sci. 51 (2004) 3427 – 3434. doi: .[52] I. Usoskin, et al., J. Geophys. Res. 110 (2005) A12108. doi: .[53] K. Herbst, et al., J. Geophys. Res. 115 (2010) D00I20. doi: .[54] Z. H. Cigolle, S. Donow, D. Evangelakos, M. Mara, M. McGuire, Q. Meyer,Journal of Computer Graphics Techniques (JCGT) 3 (2014) 1–30. URL: http://jcgt.org/published/0003/02/01/ .[55] T. Engelhardt, C. Dachsbacher, Proceedings of the Vision, Modeling, andVisualization Conference 2008, VMV 2008, Konstanz, Germany (2008)383–388.[56] F. Johansson, et al., mpmath: a Python library for arbitrary-precisionfloating-point arithmetic (version 0.19), 2014. http://mpmath.org/ .42 ppendix A. Detailed layout of
MCPL files
It is recommended to manipulate
MCPL files through (direct or indirect)calls to functions in mcpl.h . This not only ensures consistency and simplic-ity, but also allows painless format evolution where clients simply need to ob-tain updated versions of mcpl.h and mcpl.c in order to support new versionsof the
MCPL file format, while typically retaining full backwards compatibility.Nonetheless, this appendix provides reference information concerning the exactbinary layout of data in
MCPL files, specifically those in the third version of theformat (written by
MCPL starting from version 1.1.0). The file format version ofa given
MCPL file is encoded in the first few bytes of a file (see below), is printedupon inspection with the mcpltool command line utility, and is available pro-grammatically with the function mcpl_hdr_version and macro
MCPL_FORMATVERSION (cf. Appendix C). Note that all floating point numbers in
MCPL files must ad-here to the
IEEE-754 floating point standard [7] for single- (32 bit) or double-precision (64 bit) as relevant. All signed integers in the file must follow theubiquitous “two’s complement” representation.The first 48 bytes of an
MCPL file follow a fixed layout, as indicated in Ta-ble A.1, providing flags and values needed to read and interpret both the particledata section and the remainder of the header section. Note that abbreviationsused for type information in tables in this appendix are UINT32 and UINT64for 32 and 64 bit unsigned integers respectively, INT32 for 32 bit signed integersand FP64 for double-precision 64 bit floating point numbers.The layout of the second part of the header section is indicated in Table A.2,and includes both optional and repeated entries and entries with flexible length.Entries with type listed as “DATA ARRAY” are arbitrary length byte-arrays inwhich the first four bytes are unsigned 32 bit integers indicating the byte lengthof the data payload which follows just after. Note that strings are stored likeany other data, with the only twist being that terminating
NULL characters are not stored.Next, after the header section, the remainder of the file consists of the particle43 eader layout (first part)
Position Size Type Description
MCPL file. Value is always0x4d43504c ("MCPL" in
ASCII ).4 3B - File format version encoded as 3 digit zero-prefixed
ASCII num-ber (e.g. "003").7 1B - Endianness of numbers in file, 0x4c (
ASCII “L”) for little or 0x42(
ASCII “B”) for big .8 8B UINT64 Number of particles in file.16 4B UINT32 Number of custom comments in file,
NCMTS .20 4B UINT32 Number of custom binary data blobs in file,
NBLOBS .24 4B UINT32 Flag signifying whether user-flags are enabled (0x1) or disabled(0x0).28 4B UINT32 Flag signifying whether polarisation vectors are enabled (0x1)or disabled (0x0).32 4B UINT32 Flag signifying whether floating point numbers in the particledata section are in double- (0x0) or single-precision (0x1).36 4B INT32 Value of universal PDG code. A value of 0x0 means that parti-cles in the file all have their own PDG code field.40 4B UINT32 Data length per particle (redundant information, as it can beinferred from the other flags and values).44 4B UINT32 Flag signifying whether a universal weight is present in theheader (0x1) or if particles in the file all have their own weightfield (0x0).
Table A.1: Detailed layout of the first part of the header section of an
MCPL file. eader layout (second part) Size Type Description
0B or 8B FP64 Value of the universal weight, if enabled.4B+ DATA ARRAY Data is a string holding the user provided “Source name”. Thisis always present, but might be empty.
NCMTS × NBLOBS × NBLOBS × Table A.2: Detailed layout of the second part of the header section of an
MCPL file. Thepresence and count of entries here depends on values found in the first part of the headersection (cf. Table A.1). data section. It contains one entry for each particle in the file, with detailedlayout of each as indicated in Table A.3. Here, FP is either single- (32 bit)or double-precision (64 bit) numbers, depending on the relevant flag in the fileheader. Concerning the 3 floating point numbers used to represent the packeddirection vector and kinetic energy, the scheme is as discussed in Appendix B:the first and second of the 3 floating point numbers are respectively FP1 andFP2 from Table B.1, while the third is a number whose magnitude is given bythe particle’s kinetic energy and whose sign bit is used to store the last bit ofinformation needed for the direction vector, indicated in the “+1 bit” columnof Table B.1. Specifically, the sign bit is set when the number indicated in the“+1 bit” column is negative. 45 article data layout
Presence Count & type Description
OPTIONAL 3 × FP Polarisation vector (if enabled in file).ALWAYS 3 × FP Position vectorALWAYS 3 × FP Packed direction vector and kinetic energy.ALWAYS 1 × FP Time.OPTIONAL 1 × FP Weight (if file does not have universal weight).OPTIONAL 1 × INT32 PDG code (if file does not have universal PDG code).OPTIONAL 1 × UINT32 User-flags (if enabled in file).
Table A.3: Detailed layout of the data associated with each particle in an
MCPL file. ppendix B. Unit vector packing From a purely mathematical perspective, it is trivial to “pack” unit vectorsspecified in three-dimensional Cartesian coordinates into a two-dimensional co-ordinate space, and it can for instance be achieved by the standard transforma-tion between Cartesian and Spherical coordinates. However, when consideringfloating point numbers rather than the ideal mathematical abstraction of thecomplete set of real numbers, issues of numerical imprecision during the packingand subsequent unpacking transformations become crucial. Where it matters,the discussion in this appendix will, like the
MCPL format in general, assumefloating point numbers to adhere to the relevant
IEEE floating point standard [7].A few different packing algorithms will be compared in this appendix. Firstof all, the spherical transformation between Cartesian unit vectors ( u x , u y , u z ) and spherical coordinates ( θ, φ ) is investigated, using the C math library func-tions acos , atan2 , sin and cos in an obvious manner when implementing thetransformations. In addition to the significant computational overhead involvedwhen evaluating trigonometric functions, numerical uncertainties also tend toblow up when a Cartesian coordinate is very small. For instance, consider aunit vector with u z = (cid:15) for some | (cid:15) | (cid:28) . Then, θ = arccos( (cid:15) ) which is to lowestorder equal to π/ − (cid:15) , a subtraction which out of necessity will cause a lossof precision when it is stored as a floating point number: if (cid:15) is N orders ofmagnitude smaller than π/ , then the stored result will be insensitive to the N least significant digits of information in the storage of (cid:15) . The calculation of φ suffers from similar problems.The next packing algorithm to be considered is what will be denoted the Static Projection method in the following. It represents the straight-forwardand widespread solution of storing two Cartesian components, u x and u y , di-rectly and recovering the magnitude of the third by the expression | u z | = (cid:113) − u x − u y . This method, which incidentally is the one used internally inthe SSW files produced by
MCNP (cf. section 3.2), requires a single bit of addi-tional storage to be available, in order to recover the sign of u z as well as its47agnitude. This requirement of an extra bit of storage is seen in several packingschemes and is not necessarily a problem, as will be discussed later. What isproblematic, however, is that the calculation of (cid:113) − u x − u y will result in largenumerical uncertainties when the magnitude of u z is small, as it implies the sub-traction of large and nearly equal quantities and a resulting loss of significantdigits in the result.Finally, inspired by a recent survey of unit vector packing techniques [54],a packing scheme using a so-called Octahedral Projection [55] was also inves-tigated. It is a variant of the Static Projection method in which the originalpoint on the unit sphere is first projected onto an octahedral surface beforethe resulting x and y coordinates are stored. Thus, the stored variables are ( o x , o y ) = ( u x /n, u y /n ) where n = | u x | + | u y | + | u z | . Unpacking proceeds byfirst recovering the point on the octahedron as ( o x , o y , − | o x | − | o y | ) , and thenprojecting back out to the unit sphere with a simple normalisation, requiringthe evaluation of a square root. Although the algorithm has improved perfor-mance over the methods already discussed, it once again suffers from numericalprecision issues when applied to certain unit vectors. Consider for instance theunit vector ( √ − (cid:15) , , (cid:15) ) for a very small but positive epsilon. Packing willresult in o y = 0 and o x = 1 / ( √ − (cid:15) + (cid:15) ) ≈ − (cid:15) , storage of which discardsthe N least significant digits of (cid:15) , when (cid:15) is N orders of magnitude smaller thanunity. Finally, it should be mentioned that like the Static Projection method,the Octahedral Projection method also needs the sign of u z stored in a bit else-where. The present discussion thus ignores the method described in [54, 55]of encoding the sign of u z by folding ( o x , o y ) into a separate part of the planewhen u z < , as this operation introduces additional numerical imprecision andan undesired asymmetry into the algorithm.Although the three pre-existing packing methods discussed so far all rep-resent potentially useful approaches to unit vector packing depending on theparticular needs of a given use-case, they were nonetheless deemed undesirablefor a general purpose scientific format like MCPL . This is because
MCPL is intendedfor usage in a wide variety of simulation scenarios, including those where some48omponents of the directional unit vectors of stored particles could be truly mi-nuscule but non-zero in magnitude. In order to better fulfil the requirements, anew unit vector packing algorithm was devised for
MCPL . It is tentatively named
Adaptive Projection Packing , and will be presented in the following. The algo-rithm is based upon the simple observation that although the Static Projectionmethod has issues when | u z | (cid:28) , it provides very precise results when u z is thecomponent with the largest magnitude. Thus, rather than always storing u x and u y , performance can be significantly improved by letting the packing algo-rithm select the two components with the smallest magnitudes and store those.That leaves the issue that the unpacking algorithm must be able to recognisewhich components were stored. The solution to that is based upon the fact thatall three components can have at most unit magnitude, and that | u z | ≤ / √ whenever it must be stored, implying that /u z will have a magnitude largerthan . Thus, by storing /u z in place of the one of u x and u y which is greater inmagnitude, the unpacking algorithm can easily tell which components are storedin which of the two packed numbers, merely by looking at how their magnitudescompare to unity (the code should obviously store the floating point represen-tation of ∞ when u z = 0 ). The resulting encoding scheme is summarised inTable B.1: one of three storage scenarios is picked by the packing code, de-pending on the magnitudes of the three unit vector components. The uniquepacked signature of each scenario, listed in the last column, enables them to beeasily distinguished by unpacking code. Although obviously more complicatedthan the Static Projection method, the added computational cost of using theAdaptive Projection method is at most a few branches and a division, which isfound in the context of MCPL to be comparable to the cost of the OctahedralProjection method and much faster than the Spherical coordinate method dueto the expensive trigonometric function calls. In any case, the overhead is foundto be insignificant given the intended usage in
MCPL .As indicated in Table B.1, the issue of needing one extra bit of storage spaceto hold the sign of the component which was projected away, has naturally beeninherited from the Static Projection method. Although finding an unused bit49 daptive Projection Packing
Scenario FP1 FP2 + | u x | largest /u z u y sign( u x ) | FP1 | > , | FP2 | < | u y | largest u x /u z sign( u y ) | FP1 | < , | FP2 | > | u z | largest u x u y sign( u z ) | FP1 | < , | FP2 | < Table B.1: Breakdown of the Adaptive Projection Packing method, in which a unit vector, ( u x , u y , u z ) is stored into two floating point numbers, FP1 and FP2, and one extra bit ofinformation. could certainly present a challenge for some applications of unit vector packing,it is not actually an issue for the MCPL format where an extra bit of informationcan be encoded into the sign of the floating point number otherwise used to storethe kinetic energy of the particle (cf. Table 2). This bit is otherwise unused,since kinetic energy is a non-negative quantity by definition, and is availablefor usage even if particle states are specified with a kinetic energy of exactlyzero. The latter follows from the
IEEE floating point standard [7], which ensuresavailability of the sign bit even for zero (i.e. it provides a signed zero floatingpoint implementation with distinct bit patterns for − and +0 ). Although it isthus not needed for MCPL , it should be noted for completeness that it is straight-forward to implement a variant of the Adaptive Projection method which simplyallocates the “extra” bit of storage internally in one of the two resulting floatingpoint numbers. This can be achieved by encoding the information into the leastsignificant bit of the significand of either FP1 or FP2, thus sacrificing a small –but hopefully inconsequential – amount of precision in the process.In order to compare the performance of the considered packing algorithmsquantitatively, the packing precision when a unit vector u is transformed by thepacking and unpacking into u p is defined in the following. First, the packingprecision of a single component, i ∈ { x, y, z } , is defined as δ i ≡ min(1 , | u pi /u i |− ,except when u i = 0 in which case it is defined to be if u pi = u i , and otherwise.The packing precision for the entire vector is then defined as being the worst50recision of any component, δ ≡ max i δ i . The resulting values will all lie in theunit interval and if for instance δ = 10 − , the packing algorithm in questioncan be said to have preserved the vector with at least 16 significant digits inall components. The reason for not picking a somewhat simpler figure of merit,such as the angular separation between u and u p , is that it would be insensitiveto components which are very small in magnitude.To be able to provide meaningful results for all considered packing methods,the analysis presented in the following was carried out using software [56] capableof arbitrary precision floating-point arithmetic. Additionally it should be notedthat, as the focus of the present analysis is on precision and storage consumptionfor a format like MCPL , the algorithms packing numbers into single-precision (32bit) floating point storage are actually implemented using double-precision (64bit) floating point code for all intermediate calculations. For other particularuse-cases, such as graphics rendering with surface normals on a particular GPU,one could of course imagine also implementing the packing code itself usingsingle-precision everywhere.Motivated by the fact that the considered packing methods all have par-ticular issues or strategies when u z is small in magnitude, Figure B.1 illus-trates their performance when applied to vectors with specific (possibly tiny)values of u z . For each given value of u z , a test set of unit vectors isformed by sampling an azimuthal angle φ uniformly in [0 , π ) , and letting ( u x , u y ) = (cid:112) − u z · (cos φ, sin φ ) . Each packing method is then applied tothe test set, and the average value of the resulting packing precisions is plotted.The plot clearly confirms that while the Octahedral Projection method outper-forms the Spherical coordinate and Static Projection methods, they all degradein performance as the magnitude of u z tends to zero. This is qualitatively differ-ent from the Adaptive Projection method, which shows constant performanceover the entire range, at a level which is practically indistinguishable from thecase of not using any packing – i.e. storing the Cartesian coordinates directlyinto 3 dedicated floating point numbers. Note that for clarity the plot is onlyshown here for positive u z larger than − , but it was verified in a full analysis51 − − − − − − − − − − Third coordinate of unit vector − − − − − − − − − A v e r ag e p a c k i n g p r ec i s i o n b i t F P b i t F P Unit vector packing scheme
No packing (3FP)Adaptive Projection (2FP+1bit)Adaptive Projection (2FP)Spherical coordinates (2FP)Octahedral Projection (2FP+1bit)Static Projection (2FP+1bit)
Figure B.1: Average packing precision of the considered unit vector packing methods, asdiscussed in the text, for both double-precision (64 bit, full curves) and single-precision (32bit, dashed opaque curves). The Adaptive Projection (blue, green) and No packing (black)curves are almost coinciding. that, as expected, the flat level of the Adaptive Projection curves is independentof the sign of u z , and continues unchanged over the entire dynamic range of nor-malised floating point numbers, down to about − and − respectivelyfor single- and double-precision.As a natural figure of merit, Table B.2 shows the average packing preci-sion of the considered methods when applied to a sample of unit vectorssampled at random from an isotropic distribution. Although all methods pro-vide some level of precision, it is clear that only Adaptive Projection Packingprovides performance comparable to not using packing at all. In fact, for single-precision storage, Adaptive Projection Packing in the variant relevantfor MCPL , even seems to be outperforming the case of not using any packing. Thissomewhat counter-intuitive result is understood to arise from the fact that the52 verage packing precision in isotropic sample
32 bit FP 64 bit FP
Static projection (2FP+1bit) . · − . · − Spherical coordinates (2FP) . · − . · − Octahedral packing (2FP+1bit) . · − . · − Adaptive projection (2FP) . · − . · − Adaptive projection (2FP+1bit) . · − . · − No packing (3FP) . · − . · − Table B.2: The average observed packing precision of each considered method when tested on unit vectors sampled at random from an isotropic distribution, using either single- (32 bit)or double-precision (64 bit) floating point numbers for storage of the packed representations. unpacking code is implemented in full double-precision, allowing the implicitusage of the condition u x + u y + u z = 1 for recovery of a component magnitudeto inject a small amount of added precision into the final result. For double-precision storage, the precision of Adaptive Projection Packing is also seen tobe practically indistinguishable from the no packing scenario.Potentially just as crucial for users of MCPL as the average performance,Table B.3 shows the worst packing precision of any of the sampled vectors.Here, the pre-existing packing methods, listed in the first three rows, all exhibitvalues far from their average precision in Table B.2. This indicates the existenceof vectors in the sample for which the performance of the algorithms break down– a feature which was anticipated from the preceding discussion of flaws in thosealgorithms. On the other hand, the worst encountered packing precision of theconsidered Adaptive Projection Packing variants are all reasonably close to thecorresponding average values listed in Table B.2. This once again indicates therobustness of those packing algorithms, as they indeed do not seem to sufferfrom breakdown on certain domains.In summary, the presented Adaptive Projection Packing method and in par-ticular the variant adopted for
MCPL , provides a packing precision53 orst packing precision in isotropic sample
32 bit FP 64 bit FP
Static projection (2FP+1bit)
Spherical coordinates (2FP) . · − Octahedral packing (2FP+1bit) . · − Adaptive projection (2FP) . · − . · − Adaptive projection (2FP+1bit) . · − . · − No packing (3FP) . · − . · − Table B.3: The worst observed packing precision of each considered method when tested on unit vectors sampled at random from an isotropic distribution, using either single- (32 bit)or double-precision (64 bit) floating point numbers for storage of the packed representations. which is for all practical purposes comparable to that given by three floatingpoint numbers. The difference in performance is so small that it is hard to imag-ine use-cases which would be satisfied with the latter but not by the former.For that reason, it was decided to not add a no-packing option to MCPL at thispoint. 54 ppendix C. C -API for MCPL files
For reference, this appendix documents the available functions, data struc-tures and constants in the API exposed in mcpl.h as of
MCPL version 1.1.0. Thefile and its contents can be directly included and used from code compiled withany official standard of either C or C++ . Appendix C.1. Data structures and constants
Pre-processor macros providing information about the version of MCPL. They should hope-fully be self-explanatory except perhaps
MCPL_VERSION which encode the version into a singleinteger, 10000*MAJOR+100*MINOR+PATCH, allowing easy comparison of versions numbers.
Pre-processor macro providing the file format version of files written by the installation of
MCPL .Note that this file format version (currently 3) is not the same as the version number of the
MCPL distribution (currently 1.1.0), and the latter is expected to be updated more frequently.The current third version of the file format is described in detail in Appendix A, and the fileformat of a given
MCPL file can be queried through the function mcpl_hdr_version describedin Appendix C.3. typedef struct {double ekin; /* kinetic energy [MeV] */double polarisation [3]; /* polarisation vector */double position [3]; /* position [cm] */double direction [3]; /* momentum direction (unit vector ) */double time; /* time -stamp [ millisecond ] */double weight ; /* weight or intensity */int32_t pdgcode ;uint32_t userflags ;} mcpl_particle_t ;
Data structure representing a particle. Pointers to mcpl_particle_t instances are respectivelyreturned from and passed to mcpl_read (cf. Appendix C.3) and mcpl_add_particle (cf. Ap-pendix C.2) when extracting or adding particles. Refer to the description of those functionsfor further details. typedef struct { void * internal ; } mcpl_file_ttypedef struct { void * internal ; } mcpl_outfile_t
Data structures representing file handles to
MCPL files. Typically file handles are created andreturned by a call to either mcpl_open_file or mcpl_create_outfile as appropriate, andfunctions performing subsequent operations on the file all take the file handle as one of thearguments. Note that the size of these structures is identical to that of a pointer, and they aretherefore appropriate to pass or return by value from functions with no special overhead. The oid * internal pointer is, as the name implies, for purely internal usage in mcpl.c . Clientcode should never access or modify this pointer in any way, except perhaps to set it to NULLin order to mark a file handle as uninitialised or invalid. Refer to Appendix C.2 and AppendixC.3 for further details of how these file handles can be used. Appendix C.2. Functions for file creation
Creation of new
MCPL files always begins with a call to mcpl_create_outfile , re-turning a file handle of the type mcpl_outfile_t (cf. Appendix C.1). This file han-dle is then passed in to various functions in order to first set flags and meta-dataand then add particles with mcpl_add_particle . Finally, the file must be properlyclosed with a call to either mcpl_close_outfile or mcpl_closeandgzip_outfile . mcpl_outfile_t mcpl_create_outfile (const char * filename ) Function used to start creation of a new
MCPL file. It attempts to open the indicated file forwriting (overriding any existing file) and returns a handle to the caller. The latter must subse-quently be passed in as a parameter to other functions below, in order to configure the file, addparticles to it, and ultimately close it. Note that if the provided file name does not end with theextension “.mcpl” , it will be automatically appended to it (see also mcpl_outfile_filename below). const char * mcpl_outfile_filename ( mcpl_outfile_t )
Filename being written to. If the filename passed to mcpl_create_outfile ended with “.mcpl”,it will be identical to what is returned. Otherwise, a postfix of “.mcpl” will have been appended. void mcpl_hdr_set_srcname ( mcpl_outfile_t , const char *)
Optionally set name of the generating application which will be stored in the file header. Ifnot called, a string with content “unknown” will be stored instead. This function must becalled before any calls to mcpl_add_particle . void mcpl_hdr_add_comment ( mcpl_outfile_t ,const char *) Add one or more human-readable comments to the file. This function must be called beforeany calls to mcpl_add_particle . void mcpl_hdr_add_data ( mcpl_outfile_t , const char * key ,uint32_t ldata , const char * data) Add a binary data “blob” and associate it with a given key, which must be unique to the file.This function must be called before any calls to mcpl_add_particle . void mcpl_enable_userflags ( mcpl_outfile_t ) Enable per-particle user-flags to be written in the file. If not called, any non-zero value in the userflags field of added particles will be ignored by mcpl_add_particle . This functionmust be called before any calls to mcpl_add_particle . oid mcpl_enable_polarisation ( mcpl_outfile_t ) Enable per-particle polarisation vectors to be written in the file. If not called, any non-zerovalues in the polarisation field of added particles will be ignored by mcpl_add_particle .This function must be called before any calls to mcpl_add_particle . void mcpl_enable_doubleprec ( mcpl_outfile_t ) Enables double-precision (64 bit) storage of floating point numbers in particle data. De-fault is otherwise single-precision (32 bit). This function must be called before any callsto mcpl_add_particle . void mcpl_enable_universal_pdgcode ( mcpl_outfile_t , int32_t pdgcode ) Prevent per-particle PDG codes from being written in the file, letting instead all particleshave the same common code. This means that values in the pdgcode field of added particleswill be ignored by mcpl_add_particle . This function must be called before any calls to mcpl_add_particle . void mcpl_enable_universal_weight ( mcpl_outfile_t , double uw) Prevent per-particle weights from being written in the file, letting instead all particles havethe same common weight. This means that values in the weight field of added particles willbe ignored by mcpl_add_particle . Typical usage is to save one floating point of per-particlestorage, when simulation output is known to not be weighted, by setting a universal weight of . This function must be called before any calls to mcpl_add_particle . void mcpl_add_particle ( mcpl_outfile_t ,const mcpl_particle_t *) Add the particle state represented by the values on the provided mcpl_particle_t instance, tothe file. Note that some fields on the mcpl_particle_t instance might be ignored, dependingon the settings of the file (cf. the preceding descriptions of functions altering the file settings).Note that after this function has been called, the header is flushed to disk and file settings canconsequently no longer be modified – attempting to do so will result in an error. Also notethat the mcpl_get_empty_particle function described below provides a convenient mannerto obtain a mcpl_particle_t instance to modify and pass to this function, without having toworry about memory management. void mcpl_close_outfile ( mcpl_outfile_t )
Close file and flush all pending data to disk, thus representing the last but non-optional stepin the creation of a new
MCPL file. It is undefined behaviour to attempt to use a file handle foranything after passing it to this function. int mcpl_closeandgzip_outfile ( mcpl_outfile_t )
A convenience function which first calls mcpl_close_outfile and then mcpl_gzip_file (cf. Appendix C.4), forwarding the return value of the latter. This means that, on platformswhere this is possible, the function will return 1 to indicate that the newly created
MCPL filewill have been compressed with
GZIP and the final output file will have had “.gz” appended toits name. If the function returns 0, the
MCPL file will still have been created, but it will nothave been compressed with
GZIP . cpl_particle_t * mcpl_get_empty_particle ( mcpl_outfile_t ) Convenience function which returns a pointer to a nulled-out particle struct which can be usedto edit and pass to mcpl_add_particle . It can be reused and will be automatically deallocatedwhen the file is closed.
Appendix C.3. Functions for reading file data
Access to data in existing
MCPL files always begins with a call to mcpl_open_file ,returning a file handle of the type mcpl_file_t (cf. Appendix C.1). This file han-dle is then passed in to various functions in order to either access file settingsand meta-data, or particle data with the mcpl_read function. The latter functionreads the particle at the current position and skips forward to the next one,but functions of course also exist which facilitate skipping and seeking to anyposition in the file. Finally, to ensure release of all resources, the file handleshould be released with a call to mcpl_close_file . mcpl_file_t mcpl_open_file (const char * filename ); Open indicated file and load header information into memory, skip to the first (if any) availableparticle. When
ZLIB support is enabled (cf. section 2.5), this function can read compressed
MCPL files with the .mcpl.gz extension directly. unsigned mcpl_hdr_version ( mcpl_file_t )
File format version of the
MCPL file. See the description of the
MCPL_FORMATVERSION macro inAppendix C.1 for more details. uint64_t mcpl_hdr_nparticles ( mcpl_file_t )
Number of particles stored in the file. const char* mcpl_hdr_srcname ( mcpl_file_t )
Name of the generating application. unsigned mcpl_hdr_ncomments ( mcpl_file_t )
Number of comments stored in the file. const char * mcpl_hdr_comment (mcpl_file_t , unsigned i)
Access i ’th comment stored in the file. If i is not a number smaller than the number ofcomments in the file (cf. mcpl_hdr_ncomments ), an error results. int mcpl_hdr_nblobs ( mcpl_file_t ); Number of binary data “blobs” stored in the file. onst char ** mcpl_hdr_blobkeys ( mcpl_file_t ) Returns a list of the binary data keys available in the file. The function returns NULL if thereare no binary data “blobs” in the file. int mcpl_hdr_blob (mcpl_file_t , const char* key ,uint32_t * ldata , const char ** data)
Access the binary data array associated with a given key. The function returns 0 if the keydoes not exist in the file (cf. mcpl_hdr_blobkeys ). Otherwise it returns 1 and ldata and data will have been modified to contain respectively the length of the data and the address of thedata. int mcpl_hdr_has_userflags ( mcpl_file_t );
Returns 1 if per-particle user-flags are stored in the file. If not, 0 is returned and all particlesread from the file will have a userflags field value of . int mcpl_hdr_has_polarisation ( mcpl_file_t ); Returns 1 if per-particle polarisation vectors are stored in the file. If not, 0 is returned and allparticles read from the file will have a null vector in the polarisation field. int mcpl_hdr_has_doubleprec ( mcpl_file_t );
Returns 1 if floating point numbers in the particle data in the file is stored using double-precision (64 bit) as opposed to single-precision (32 bit) numbers. uint64_t mcpl_hdr_header_size ( mcpl_file_t );
Returns the number of bytes consumed by the header on disk (uncompressed). int mcpl_hdr_particle_size ( mcpl_file_t );
Returns the number of bytes consumed by each particle on disk (uncompressed). int32_t mcpl_hdr_universal_pdgcode ( mcpl_file_t );
Returns zero if per-particle PDG codes are stored in the file. If not, the returned value is thecommon value which all particles read from the file will have in the pdgcode field. double mcpl_hdr_universal_weight ( mcpl_file_t );
Returns zero if per-particle weights are stored in the file. If not, the returned value is thecommon value which all particles read from the file will have in the weight field. int mcpl_hdr_little_endian ( mcpl_file_t );
Returns 1 if the numbers in the file are stored in little-endian form, and 0 if they are stored inbig-endian form. const mcpl_particle_t * mcpl_read ( mcpl_file_t ); ttempts to read the particle at the current position in file and skips forward to the nextparticle. Return value will be NULL in case there was no particle at the current location(normally due to end-of-file), otherwise it will be a pointer to an mcpl_particle_t instancerepresenting the particle just read. Note that the returned pointer is invalidated if the file isclosed or mcpl_read is called again. uint64_t mcpl_currentposition ( mcpl_file_t ); Returns current position in the file, which is a number less than or equal to the number ofparticles in the file, N (cf. mcpl_hdr_nparticles ). If N is returned, this indicates an end-of-file condition where no more particles can be read with mcpl_read . int mcpl_skipforward (mcpl_file_t , uint64_t n); Advance position in file n steps. Returns 0 if this causes the end-of-file to be reached and thereis no particle to be read at the new position. Otherwise returns 1. int mcpl_rewind ( mcpl_file_t ); Rewinds position in file to 0. Returns 0 if this causes the end-of-file to be reached and there isno particle to be read at the new position (under normal conditions, this should only happenfor empty files with no particles). Otherwise returns 1. int mcpl_seek (mcpl_file_t , uint64_t ipos );
Seek directly to specified position in file. Returns 0 if this causes the end-of-file to be reachedand there is no particle to be read at the new position. Otherwise returns 1. void mcpl_close_file ( mcpl_file_t );
Deallocate memory and release file-handle. It is undefined behaviour to attempt to use a filehandle for anything after passing it to this function.
Appendix C.4. Other functions void mcpl_dump (const char * file , int parts , uint64_t nskip , uint64_t nlimit )
Prints information about the specified
MCPL file to standard output. This is similar to what isprinted by the mcpltool at the command line. The parameter parts can be used to controlwhether or not to print information from just the file header ( parts=1 ), just the particle statedata ( parts=2 ) or both ( parts=0 ). If particle state data is listed, nskip and nlimit canbe used to control which of the contained particles to list: nlimit is an upper limit on thenumber of particles printed ( nlimit=0 means no limit), and nskip can be used to skip thatmany positions into the file before starting the printouts. mcpl_outfile_t mcpl_merge_files (const char* file_output ,unsigned nfiles , const char ** files) erge contents of a list of files by concatenating all particle contents into a new output file.This results in an error unless all meta-data and settings in the files are identical. For safety,this fails if file_output already exists. The function returns a handle to the output file whichhas had all particles added to it, but has not yet been closed. Note that if any file is specifiedmore than once in the input list, a warning will be printed to standard output. int mcpl_can_merge (const char * file1 , const char* file2) Test if files could be merged by mcpl_merge_files . This returns 1 if all meta-data and settingsin the files are identical, otherwise 0. void mcpl_merge_inplace (const char * file1 , const char* file2 );
Similar to mcpl_merge_files , but merges two files by appending all particles in file2 to thelist in file1 (thus file1 grows while file2 stays untouched). Note that this requires similarversion of the file format of the two files, in addition to the other checks in mcpl_can_merge .Careful usage of this function can be more efficient than mcpl_merge_files , but it is poten-tially also less safe as file1 is left modified. Note that if file1 and file2 are the same file,a warning will be printed to standard output. void mcpl_repair (const char * file1)
Attempts to repair a broken file which was never properly closed. This is intended for recoveryof contents in files produced in long jobs which were interrupted for some reason, and thusnever had the number of particles field in the file header updated correctly. It works by usingthe file size to calculate the number of complete particle entries in the file, and then updatingthe header. int mcpl_tool (int argc , char ** argv)
This function implements the command line mcpltool command, and should be wrapped insidea standard main function of a C application, which should then be compiled into an executablenamed mcpltool or similar. int mcpl_gzip_file (const char * filename ); Attempts to compress file with
GZIP , appending “.gz” to its name in the process. On platformswhere this is possible, the function will return 1 to indicate success. If the function insteadreturns 0, the file will have been left untouched. void mcpl_transfer_metadata ( mcpl_file_t source , mcpl_outfile_t target );
Convenience function which transfers all settings, blobs and comments to outfile. Intended tomake it easy to filter files via custom C code. void mcpl_set_error_handler (void (* handler )( const char *));
Override the default
MCPL error handler, which is a function that will get called with a stringdescribing the error and which should never return to the calling code. If no handler is set,errors will get printed to standard output and the process terminated. ppendix D. Usage instructions for mcpltool Usage instructions of the mcpltool command are available from the com-mand line by invoking it with the --help flag. For reference the output isrepeated here in Listing D.1.
Listing D.1: Usage instructions for the mcpltool command (output of “ mcpltool --help ”).