# MURE 2: SMURE, Serpent-MCNP Utility for Reactor Evolution User Guide -Version 1 Main Contributors

O. Méplan, Jan Hajnrych, A. Bidaud, S. David, N. Capellan, B. Leniau, A. Nuttin, Frantisek Havluj, Radim Vocka, J. Wilson, et al.

## ► To cite this version:

O. Méplan, Jan Hajnrych, A. Bidaud, S. David, N. Capellan, et al. MURE 2: SMURE, Serpent-MCNP Utility for Reactor Evolution User Guide -Version 1 Main Contributors. [Technical Report] LPSC17002, Laboratoire de Physique Subatomique et de Cosmologie. 2017. in2p3-01491116

HAL Id: in2p3-01491116

https://hal.in2p3.fr/in2p3-01491116

Submitted on 16 Mar 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - NonCommercial - ShareAlike 4.0 International License
MURE 2 : SMURE,
Serpent-MCNP Utility for Reactor Evolution

User Guide - Version 1

<figure><img src="image_2.png" /></figure>

February 2017

**Main Contributors :**

* O. Méplan (alias PTO), LPSC Grenoble

* Jan Hajnrych, Ecole des Mines de Nantes, and Warsaw University of Technology (Poland)

* A. Bidaud (alias le bid), LPSC Grenoble

* S. David (alias GTS), IPN Orsay

* N. Capellan (alias Nico la Star), LPSC Grenoble

* B. Leniau (alias BLG), Subatech Nantes

* A. Nuttin (alias Nut), LPSC Grenoble

* Frantisek Havluj, UJV, Czech Republic

* Radim Vocka, UJV, Czech Republic

* and in the past:

- J. Wilson (alias JW), IPN Orsay

- R. Chambon (alias Le caribou), LPSC, Grenoble, now back to Canada
- F. Michel-Sendis (alias FMS), IPN Orsay, Now @NEA

- F. Perdu (alias WEC), LPSC, Grenoble, Now @CEA

- L. Perrot, IPN Orsay

See also the FAQ
# Contents

<table><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Introduction</td><td>8</td></tr><tr><td>1.1</td><td>Installation</td><td>10</td></tr><tr><td>1.1.1</td><td>Compilation</td><td>10</td></tr><tr><td>1.1.2</td><td>Remarks on MCNP/makxsf compilation</td><td>11</td></tr><tr><td>1.1.3</td><td>Remarks on MURE with Serpent2</td><td>11</td></tr><tr><td>1.1.4</td><td>Building files for evolution</td><td>11</td></tr><tr><td>1.1.5</td><td>Running some examples</td><td>12</td></tr><tr><td>1.1.5.1</td><td>basic MURE possibilities</td><td>12</td></tr><tr><td>1.1.5.2</td><td>Fuel Evolution in a Sphere</td><td>13</td></tr><tr><td>1.2</td><td>MURE Package structure</td><td>14</td></tr><tr><td>1.3</td><td>MURE class</td><td>15</td></tr><tr><td>1.4</td><td>MURE basic files</td><td>16</td></tr><tr><td>1.4.1</td><td>Files of "MURE/documentation"</td><td>16</td></tr><tr><td>1.4.2</td><td>Files in "MURE/data"</td><td>16</td></tr><tr><td>1.4.3</td><td>Files of "MURE/utils"</td><td>17</td></tr><tr><td>1.4.4</td><td>Graphical User Interface "MURE/gui"</td><td>18</td></tr><tr><td>1.4.5</td><td>MURE Source files "MURE/source/include" and "MURE/source/src"</td><td>18</td></tr><tr><td>1.4.6</td><td>Other files</td><td>19</td></tr><tr><td>2</td><td>What's new in MURE</td><td>20</td></tr><tr><td>3</td><td>From MURE to SMURE : Choosing the Monte-Carlo Transport Code</td><td>25</td></tr><tr><td>3.1</td><td>Switching from MCNP to Serpent in a MURE input file</td><td>25</td></tr><tr><td>3.1.1</td><td>Principle of the implementation</td><td>25</td></tr><tr><td>3.1.2</td><td>A example</td><td>26</td></tr><tr><td>3.2</td><td>How migrate my old MURE V1.x file to the MURE V2.x file</td><td>26</td></tr><tr><td>3.2.1</td><td>What has definitely change</td><td>26</td></tr><tr><td>3.2.2</td><td>What is still existing but can be done in a more elegant way</td><td>28</td></tr><tr><td>4</td><td>Geometry Definition</td><td>31</td></tr><tr><td>4.1</td><td>Introduction</td><td>31</td></tr><tr><td>4.2</td><td>Definition of geometrical shapes</td><td>31</td></tr><tr><td>4.2.1</td><td>Units</td><td>32</td></tr></tbody></table>
<table>
  <tbody>
    <tr>
      <td>4.2.2 Examples of simple shapes</td>
      <td>32</td>
    </tr>
    <tr>
      <td>4.2.3 Examples of simple Nodes</td>
      <td>33</td>
    </tr>
    <tr>
      <td>4.2.4 Moving a Shape</td>
      <td>33</td>
    </tr>
    <tr>
      <td>4.3 The "put in" operator &gt;&gt;</td>
      <td>33</td>
    </tr>
    <tr>
      <td>4.4 Clone of a Shape</td>
      <td>34</td>
    </tr>
    <tr>
      <td>4.5 Boundary conditions</td>
      <td>35</td>
    </tr>
    <tr>
      <td>4.6 Recommendations</td>
      <td>35</td>
    </tr>
    <tr>
      <td>4.7 Definition of MC cells</td>
      <td>36</td>
    </tr>
    <tr>
      <td>4.7.1 Cell and clone shapes</td>
      <td>37</td>
    </tr>
    <tr>
      <td>4.8 Lattice</td>
      <td>37</td>
    </tr>
    <tr>
      <td>4.8.1 An implicit lattice example (SimpleLattice.cxx &amp; SimpleLattice_serpent.cxx)</td>
      <td>38</td>
    </tr>
    <tr>
      <td>4.8.2 An explicit lattice example with different zones (SimpleLattice2.cxx &amp; SimpleLattice2_serpent.cxx)</td>
      <td>40</td>
    </tr>
    <tr>
      <td>4.8.3 A lattice with more than one simple shape (Stadium.cxx &amp; Stadium_serpent.cxx)</td>
      <td>41</td>
    </tr>
    <tr>
      <td>4.8.4 Lattice of a Lattice (LatticeOfLattice.cxx &amp; LatticeOfLattice_serpent.cxx)</td>
      <td>43</td>
    </tr>
    <tr>
      <td><b>5 Materials, Sources, Tallies, ...</b></td>
      <td><b>47</b></td>
    </tr>
    <tr>
      <td>5.1 Definition of material</td>
      <td>47</td>
    </tr>
    <tr>
      <td>5.1.1 Clone</td>
      <td>47</td>
    </tr>
    <tr>
      <td>5.1.2 Mix (Mix.cxx)</td>
      <td>48</td>
    </tr>
    <tr>
      <td>5.1.3 Units</td>
      <td>48</td>
    </tr>
    <tr>
      <td>5.1.4 Material extension for MC</td>
      <td>49</td>
    </tr>
    <tr>
      <td>5.1.5 Pseudo Materials</td>
      <td>50</td>
    </tr>
    <tr>
      <td>5.1.6 MC material Printing</td>
      <td>51</td>
    </tr>
    <tr>
      <td>5.1.7 Examples</td>
      <td>51</td>
    </tr>
    <tr>
      <td>5.1.8 Automatic extension finding</td>
      <td>51</td>
    </tr>
    <tr>
      <td>5.1.9 Automatic XSDIR construction</td>
      <td>52</td>
    </tr>
    <tr>
      <td>5.2 Particle Source</td>
      <td>52</td>
    </tr>
    <tr>
      <td>5.2.1 Setting Source for MURE</td>
      <td>53</td>
    </tr>
    <tr>
      <td>5.2.2 Criticallity calculation: Kcode mode</td>
      <td>53</td>
    </tr>
    <tr>
      <td>5.2.3 More elaborated sources (MCNP only)</td>
      <td>53</td>
    </tr>
    <tr>
      <td>5.2.3.1 Source define via Spectrum class (MCNP only)</td>
      <td>54</td>
    </tr>
    <tr>
      <td>5.2.3.2 Tube Source (MCNP only)</td>
      <td>54</td>
    </tr>
    <tr>
      <td>5.2.3.3 Define MCNP Source from the result of an evolution (MCNP only)</td>
      <td>54</td>
    </tr>
    <tr>
      <td>5.3 Tally class</td>
      <td>55</td>
    </tr>
    <tr>
      <td>5.3.1 Fluence to dose conversion (MCNP Only)</td>
      <td>56</td>
    </tr>
    <tr>
      <td><b>6 Nuclei Tree</b></td>
      <td><b>58</b></td>
    </tr>
    <tr>
      <td>6.1 Nuclei Trees: general considerations</td>
      <td>58</td>
    </tr>
    <tr>
      <td>6.1.1 How the tree information is stored</td>
      <td>58</td>
    </tr>
    <tr>
      <td>6.1.2 Cutting or simplifying the tree</td>
      <td>58</td>
    </tr>
    <tr>
      <td>6.1.3 The difference/similarity of material trees</td>
      <td>60</td>
    </tr>
    <tr>
      <td>6.1.4 How many nuclei need to be considered in the evolution?</td>
      <td>63</td>
    </tr>
    <tr>
      <td>6.1.5 Fission Products</td>
      <td>63</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>6.1.6</td>
      <td>Isomer Production from (n,gamma) or (n,2n) reactions</td>
      <td>67</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Nuclei Tree: the implementation</td>
      <td>67</td>
    </tr>
    <tr>
      <td>6.2.1</td>
      <td>Important files</td>
      <td>67</td>
    </tr>
    <tr>
      <td>6.2.2</td>
      <td>Reaction Auto-detection</td>
      <td>68</td>
    </tr>
    <tr>
      <td>6.2.3</td>
      <td>The recursion depth cut</td>
      <td>70</td>
    </tr>
    <tr>
      <td><b>7</b></td>
      <td><b>Evolution in MURE</b></td>
      <td><b>72</b></td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>Preliminary Remark</td>
      <td>72</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>General Considerations</td>
      <td>72</td>
    </tr>
    <tr>
      <td>7.2.1</td>
      <td>Time discretization</td>
      <td>73</td>
    </tr>
    <tr>
      <td>7.2.2</td>
      <td>The Bateman's Equation</td>
      <td>73</td>
    </tr>
    <tr>
      <td>7.2.2.1</td>
      <td>The evolution matrix</td>
      <td>74</td>
    </tr>
    <tr>
      <td>7.2.2.2</td>
      <td>Multi-Threading parallelization</td>
      <td>74</td>
    </tr>
    <tr>
      <td>7.2.2.3</td>
      <td>Cooling period</td>
      <td>74</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Different way of evolution</td>
      <td>75</td>
    </tr>
    <tr>
      <td>7.3.1</td>
      <td>Cross-section Evolution</td>
      <td>76</td>
    </tr>
    <tr>
      <td>7.3.1.1</td>
      <td>How does it work?</td>
      <td>76</td>
    </tr>
    <tr>
      <td>7.4</td>
      <td>Reactivity Control</td>
      <td>77</td>
    </tr>
    <tr>
      <td>7.4.1</td>
      <td>A simple example using standard EvolutionControl</td>
      <td>78</td>
    </tr>
    <tr>
      <td>7.4.1.1</td>
      <td>Poison material declaration</td>
      <td>78</td>
    </tr>
    <tr>
      <td>7.4.1.2</td>
      <td>Escape calculation</td>
      <td>78</td>
    </tr>
    <tr>
      <td>7.4.1.3</td>
      <td>Evolution definition</td>
      <td>78</td>
    </tr>
    <tr>
      <td>7.4.2</td>
      <td>Using your own EvolutionControl</td>
      <td>79</td>
    </tr>
    <tr>
      <td>7.5</td>
      <td>Evolution of a MCNP user defined geometry</td>
      <td>79</td>
    </tr>
    <tr>
      <td>7.6</td>
      <td>More complex evolution conditions</td>
      <td>80</td>
    </tr>
    <tr>
      <td>7.7</td>
      <td>Equilibrium of Xe-135</td>
      <td>80</td>
    </tr>
    <tr>
      <td><b>8</b></td>
      <td><b>Looking at results from the evolution</b></td>
      <td><b>82</b></td>
    </tr>
    <tr>
      <td>8.1</td>
      <td>The MURE data files</td>
      <td>82</td>
    </tr>
    <tr>
      <td>8.2</td>
      <td>Reading the data files with Tcl/Tk graphical interface scripts</td>
      <td>82</td>
    </tr>
    <tr>
      <td>8.3</td>
      <td>Reading the data files with ROOT graphical interface</td>
      <td>83</td>
    </tr>
    <tr>
      <td>8.3.1</td>
      <td>The 8 tabs of the main window</td>
      <td>84</td>
    </tr>
    <tr>
      <td>8.3.1.1</td>
      <td>Out flux radiotoxicity evolution</td>
      <td>91</td>
    </tr>
    <tr>
      <td>8.3.2</td>
      <td>Frame "Cell Spatial Var" in the Main window</td>
      <td>91</td>
    </tr>
    <tr>
      <td>8.3.3</td>
      <td>Frame "Time Interpolation" in the main window</td>
      <td>92</td>
    </tr>
    <tr>
      <td>8.3.4</td>
      <td>Frame "Error Bars" in the main window</td>
      <td>92</td>
    </tr>
    <tr>
      <td>8.3.5</td>
      <td>The Plot/Save/Quit buttons in the main window</td>
      <td>93</td>
    </tr>
    <tr>
      <td>8.3.6</td>
      <td>More about Radiotoxicity Tab</td>
      <td>93</td>
    </tr>
    <tr>
      <td>8.3.6.1</td>
      <td>Radiotoxicity of a user input composition (cooling only)</td>
      <td>93</td>
    </tr>
    <tr>
      <td>8.3.6.2</td>
      <td>Gamma, Beta, Alpha, and Neutron spectra of an evolving material</td>
      <td>94</td>
    </tr>
    <tr>
      <td>8.3.6.3</td>
      <td>Radiotoxicity of a sample</td>
      <td>95</td>
    </tr>
    <tr>
      <td>8.4</td>
      <td>Spectrum Classes (MCNP only)</td>
      <td>95</td>
    </tr>
    <tr>
      <td>8.4.1</td>
      <td>Spectrum Class</td>
      <td>97</td>
    </tr>
  </tbody>
</table>
<table><tr><td>8.4.2 GammaSpectrum class</td><td>97</td></tr><tr><td>8.4.3 AlphaSpectrum class</td><td>98</td></tr><tr><td>8.4.4 BetaSpectrum class</td><td>98</td></tr><tr><td>8.4.5 NeutronSpectrum class</td><td>98</td></tr><tr><td>8.4.5.1 Neutron from spontaneous fission</td><td>98</td></tr><tr><td>8.4.5.2 Neutron from (α, n) reactions</td><td>99</td></tr><tr><td>8.4.5.3 Neutron from (β⁻n) decays</td><td>99</td></tr><tr><td>8.4.6 Define MCNP Source with Spectrum object</td><td>99</td></tr></table>

9 Thermal-hydraulics/neutronics coupling 100

<table><tr><td>9.1 Preliminary Remarks</td><td>100</td></tr><tr><td>9.2 General Overview</td><td>100</td></tr><tr><td>9.3 Description of methods</td><td>101</td></tr><tr><td>9.3.1 Introduction</td><td>101</td></tr><tr><td>9.3.2 ReactorAssembly class</td><td>103</td></tr><tr><td>9.3.2.1 Using a ReactorAssembly to fill a core lattice</td><td>104</td></tr><tr><td>9.3.3 The COBRA_EN class: coupled neutronics/thermal-hydraulics calculations</td><td>106</td></tr><tr><td>9.3.3.1 Coupling Thermo-hydraulics and neutronics</td><td>108</td></tr><tr><td>9.3.3.2 Input/Output file</td><td>108</td></tr><tr><td>9.3.4 The BATH class (NOT YET TESTED IN SMURE-MURE v2.0)</td><td>110</td></tr><tr><td>9.3.4.1 Capabilities</td><td>110</td></tr><tr><td>9.3.4.2 What is solved</td><td>110</td></tr><tr><td>9.3.4.3 How to add data</td><td>111</td></tr><tr><td>9.3.4.4 Use the code</td><td>111</td></tr><tr><td>9.3.4.5 Output data</td><td>111</td></tr></table>

A Basic of C++ to understand MURE 113

<table><tr><td>A.1 Introduction</td><td>113</td></tr><tr><td>A.2 Class & Object</td><td>113</td></tr><tr><td>A.2.1 Header and implementation files</td><td>115</td></tr><tr><td>A.3 Default arguments</td><td>116</td></tr><tr><td>A.4 Pointers</td><td>117</td></tr><tr><td>A.5 Inheritance</td><td>118</td></tr><tr><td>A.6 Namespace</td><td>118</td></tr></table>

B Node tree simplification 120

C Stochastic volume calculation 123

<table><tr><td>C.1 Important Note</td><td>124</td></tr></table>

D Back-Stage processes: everything you always wanted to know about MURE* (*But were afraid to ask) 125

<table><tr><td>D.1 Steady-state power: Tally Normalization</td><td>125</td></tr><tr><td>D.1.1 Energy Released from MURE fissions</td><td>125</td></tr></table>
<table><tr><td>D.1.2 How the normalization factor is calculated</td><td>126</td></tr><tr><td>D.1.3 An alternative approach (for information only)</td><td>127</td></tr><tr><td><b>Bibliography</b></td><td><b>128</b></td></tr></table>
# Chapter 1

# Introduction

In the following, we refer either to **MURE** (MCNP Utility for Reactor Evolution), **MURE 2** or **SMURE** (Serpent/MCNP Utility for Reactor Evolution)<fnref n="1" />; they represent the same package (MURE is the version 1 of SMURE/MURE 2). The main aim of the *MURE*[1, 2]/*SMURE* package is to perform nuclear reactor time-evolution using the widely-used particle transport code *MCNP*[3] (a Monte Carlo code which is mostly written in FORTRAN) or *Serpent*[2, 4, 5]. Many depletion codes exist for determining time-dependent fuel composition and reaction rates. These codes are either based on solving Boltzman equation using deterministic methods or based on Monte-Carlo method for neutron transport. Among them, one has to cite MCNPX/CINDER 90[6], MONTEBURN[7], KENO/ORIGEN[8], MOCUP[9], MCB[10], VESTA/MORET[17, 18], TRIPOLI-4D[19], Serpent[4], ... which provide neutron transport and depletion capabilities. However, the way to control (or interact with) the evolution are either limited to specific procedure and/or difficult to implement.

In (S)MURE, due to the Object-oriented programming, any user can define his own way to interact with evolution. From an academic point of view, it is also good to have lots of M-C evolution codes to compare and benchmark them to understand physics approximations of each one. Moreover, SMURE provides a simple graphical interface to visualize the results. It also provides a way to couple the neutronics (with or without fuel burn-up) and thermohydraulics using either an open source simple code developed in SMURE (BATH, Basic Approach of Thermal Hydraulics) or a sub-channel 3D code, COBRA-EN[12, 13]. But SMURE can also be used just as an interface to MCNP or Serpent to build geometries (e.g. for neutronics experiments simulation).

SMURE is based on C++ objects allowing a great flexibility in the use<fnref n="2" />. There are 4 main parts in this library:

1. Definition of the geometry, materials, neutron source, tallies, ...

2. Construction of the nuclear tree, the network of links between neighbouring nuclei via radioactive decays and nuclear reactions.

3. Evolution of some materials, by solving the corresponding Bateman's equations.

4. Thermal-hydraulics: it couples neutronics, thermal-hydraulics and, if needed, fuel evolution.

<footnote n="1">¹In french, a mûre is a blackberry; the blue fish with legs of the MURE's logo is the Darwin evolution symbol (a fish leaving the sea), with a small wrench that symbolize that MURE is a tool. The fish climb on a tree with blackberries; they represents the nuclei tree need for the evolution, each drupelet of the berries being the nucleons. The small serpent is of course link with the biblical tree of knowledge (but here for the neutronics knowledge) as well as the Serpent Monte-Carlo code used together with MCNP.</footnote>

<footnote n="2">²Basic knowledge of C/C++ may help to understand MURE. Nevertheless, careful reading of examples is sufficient to understand and use the package. In appendix A, a very short introduction to C++ terminology is given in order to facilitate the use of the User Guide and examples.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1.1: Principle of fuel evolution in MURE.</figcaption></figure>

* Part 1 can be used independently of the 2 others; it allows “easy” generation of *Serpent/MCNP* input files by providing a set of classes for describing complex geometries. The ability to make quick global changes to reactor component dimensions and the ability to create large lattices of similar components are two important features that can be implemented by the C++ interface. It should be noted that some knowledge of *MCNP* or *Serpent* is very useful in understanding the geometry generation philosophy.

* Part 2 builds the specific nuclear tree from an initial material composition (list of nuclei). The tree of each “evolving”<fnref n="3" /> nucleus is created by following the links between neighbours via radioactive decay and/or reactions until a self-consistent set of linked nuclei is extracted. Nuclei with half-lives very much shorter than the evolution time steps, could be removed from the tree; mothers and daughters of these removed nuclei are re-linked in the correct way. Part 2 can also be used independently of the other two parts to process cross-sections for *MCNP/Serpent* at the desired temperature.

* Part 3 simulates the evolution of the fuel within a given reactor over a time period of up to several years, by successive steps of *Serpent/MCNP* calculation and numerical integration of Bateman’s equations. Each time the *MC* code is called, the reactor fuel composition will change due to the fission/capture/decay process occurring inside. Changes in geometry, temperature, external feeding or extraction during the evolution can also be taken into account. Obviously this part is not independent of the 2 others<fnref n="4" /> (see figure 1.1).

* Part 4 consists of coupling the Oak Ridge National Laboratory code **COBRA-EN** (COolant Boiling in Rod Arrays) with **MURE**. **COBRA** is a sub-channel code that allows steady-state and transient analysis of the coolant in rod arrays. The simulation of flow is based on a three or four partial differential equations: conservation of mass, energy and momentum vector for the water liquid/vapor mixture (optionally a fourth equation can be added which tracks the vapor mass separately). The heat transfer model is featured by a full boiling curve, comprising the basic heat transfer regimes: single phase forced convection, sub-cooled nucleate boiling, saturated nucleate boiling, transition and film boiling. Heat conduction in the fuel and the cladding is calculated using the balance equation.

The use of this package requires the following installation:

1. a C++ compiler (mandatory). *SMURE* is developed using *gcc* (all version between 2.96 and 6.3.1 are known to work perfectly).

2. *Serpent2*, *MCNP* or *MCNPX* (mandatory, available at the NEA DataBank & RSICC). These codes will be refered as “*MC*” in the following.

<footnote n="3">³For reactor physics, it is not generally necessary to take into account the evolution of every material (such as reflectors, vessels, ...).</footnote>

<footnote n="4">⁴It is, nevertheless, relatively easy to make evolution of user defined MCNP file.</footnote>
(a) Nuclear data files process for MCNP (ACE format) and/or NJOY to process ENDF files into the ACE format (Nuclear data in ACE format are mandatory).

(b) **Serpent 1** can not be used with MURE because it does not support unions.

3. **COBRA-EN** for the thermal-hydraulic part (available the atNEA DataBank & RSICC)

4. The **ROOT** graphical tools developed at CERN (http://root.cern.ch) is necessary if user wants to use the post treatment GUI tools. In the GUI, radiotoxicity of fuel/waste can be calculated and plotted if the **LAPACK** library is installed (Linear Algebra PACKage, see http://www.netlib.org/lapack/index.html or standard LINUX repositories (Redhat, Fedora, Ubuntu, ...) for pre-installed versions). If you are using ROOT-5 or before, you should modify the beginning of the Makefile of MureGui.

At present, SMURE has only been compiled and tested on LINUX/UNIX platforms.

# 1.1 Installation

TO BE NOTICED: for some evident reasons<fnref n="5" />, upgrading for MURE version 1 to MURE version 2 breaks some backward compatibility. Nevertheless, for “old” MURE user changes will be minors; we have try to keep the essential of MURE, but given more flexibility. A chapter is devoted to the migration of old MURE users.

## 1.1.1 Compilation

* Uncompress the archive file where you want to install MURE source.

```
tar zxvf MURE_XX.tgz
```

or

```
gunzip MURE_XX.tgz ; tar xvf MURE_XX.tar
```

This will create the "data", "documentation", "examples", "gui", "lib", "source" and "utils" directories in the MURE directory.

* Configure and Compile MURE

- A bash script is given in the MURE directory

```
./install.sh -h : full flag option
```

* main options are

- `-MCNP-version=Num` where `Num` is the MCNP version (or its ACE data format) ; `Num` can be either 4 (for MCNP 4 (MCNPX <2.5) compiled with 32 bit data) or 5 (for MCNP ≥5 (or MCNPX ≥2.5) compiled with 64 bit data. MCNP Version 5 is the default.

- `-ENSDF-path=path` where `path` is a valid directory path (absolute or relative) to ENSDF main directory files ; these optional data can be obtained at http://ie.lbl.gov/databases/ensdfserve.html, download the "Complete ENSDF database" and unzip it in `/path_to_ENSDF_data/ENSDF`.

<footnote n="5">⁵MURE v1 was designed only for MCNP. It was closely linked to MCNP. The first phase of MURE v2 has been to decouple MURE from MCNP ; the second phase was to couple MURE to a MC code in a “generic” way”. This has caused the break in backward compatibility.</footnote>
* Quick install (with the default MCNP 5 version for binary ACE data format 64 bits)

./install.sh

This script will create a *config* directory that contains *Makefile.config* and *config.hxx* necessary for all *Makefiles* and *MureHeader.hxx* respectively. Then, it compiles "external" libraries (*ValErr* and *mctal*) and the *MUREpkg* library and put them in "MURE/lib". These libraries are shared libraries.

* Set the LD_LIBRARY_PATH (or SHLIB_PATH on HP-UX or LIBPATH on AIX)
  in tcsh/csh

setenv LD_LIBRARY_PATH "${LD_LIBRARY_PATH}:/path_to_MURE/MURE/lib"

or in bash/sh

export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:/path_to_MURE/MURE/lib"

### 1.1.2 Remarks on MCNP/makxsf compilation

In *MURE*, the access to ACE *MCNP* cross-sections is done for different purposes (find if cross-sections exist, find the total cross-sections in order to write in *MC* input file only "significant nucleus" (see 5.1.6), use of "multi-group flux" evolution method (see item 5 of section 7.3), ...). The use of binary ACE files improves the reading time (and also the disk space necessary to store the cross-sections). In C/C++, reading/writing in unformatted files (binary) is done by reading records of 1 byte long. But some FORTRAN compilers such as IFORT of *Intel* read/write in unformatted file with 4 bytes records. Thus if you want to use such a compiler for MCNP/makxsf you have to force the record length in unformatted files to be 1 byte. In the *Intel* IFORT compiler this is done via the *-assume byterecl* flag of *ifort* command line (see the *System.gcf* file of the *config* dir of MCNP distribution where *System* is either Linux, AIX, ...).

### 1.1.3 Remarks on MURE with Serpent2

At the present time, *Serpent2* seems only to handle ASCII file for ACE cross-section. Thus one has to use this ASCII data if one use *Serpent*. Conversion from the MCNP xsdir to Serpent xsdata is done either with the code provided with Serpent, either directly in *MURE* if one use the *MURE:SetAutoXSDIR()* method.

### 1.1.4 Building files for evolution

In this section, it is explained how to build the 2 necessary files *BaseSummary.dat* and *AvailableReactionChart.dat* (these files are not provided in the distribution). One supposes that a user has a standard *MCNP5* ENDF B6 base provided with MCNP5. The provided *xsdir* file and the nuclear data are located in */path_2_ace_files* directory (Verify that in the xsdir the absolute path to data files is present).

* Then go to *MUREutils/datadir*.

* Compile the *ExtractXsdir.cxx* file (compilation line is at the end of the file).

* Then run it with "ExtractXsdir"

- Answer *ENDFB* to the first question
- 6.8 to the second one

- /path_2_ace_files/ xsdir to the third one

- and STD to the last one.

* After a while the **BaseSummary.dat** is created and contains the following

<table><tr><td>1</td><td>1</td><td>0 ENDFB</td><td>6.1</td><td>293.62</td><td>.24c</td><td>STD</td><td>1001.24c 0.999170 /path_2_ace_files/la150n2 0 2 1 10106 4096 512 2.5301E-08</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1</td><td>1</td><td>0 ENDFB</td><td>6.1</td><td>293.62</td><td>.66c</td><td>STD</td><td>1001.66c 0.999170 /path_2_ace_files/ endf66a2 0 2 1 10128 4096 512 2.5301E-08</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>92</td><td>235</td><td>0 ENDFB</td><td>6.1</td><td>293.62</td><td>.66c</td><td>STD</td><td>92235.66c 233.025000 /path_2_ace_files/ endf66c2 0 2 6899 722105 4096 512 2.5301E-08 ptable</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></table>

* Then compile (compilation line at the end of the file) and run **CheckReaction**

* copy **AvailableReactionChart.dat** and **BaseSummary.dat** in MURE/data

* Now you can run any examples of MURE.

### 1.1.5 Running some examples

The directory MURE/examples contains commented examples on how to use MURE. When this is possible, the same example exist for both MCNP and Serpent2; as it will be emphasis in the following, the differences in these 2 examples is generally limited to 3 or 4 lines to change, whatever the complexity of the system. A README file suggests an order to execute examples. For each example, a compilation line is put at the end of the file. Here, 2 examples are described; other examples are described in this User Guide (see for example, § 4.8.1, 4.8.2, 4.8.3 and 4.5).

In examples directory, one finds static examples (i.e. without burn-up calculations) and their outputs<fnref n="6" /> (in MURE/examples/Output/ directory). To run some of these examples you must provide an MCNP xsdir file, whereas for other ones user has to build the **BaseSummary.dat** and **AvailableReactionChart.dat** files. The “evolution” examples are in MURE/examples/Evolution/ directory. Again you will find an “Output” directory that contains directory output of the examples (also gzip). In these result directories, one has suppressed the MCNP “o” and “m” files but left the MCNP input files as well as the MURE results files.

#### 1.1.5.1 basic MURE possibilities

This example (MURE/examples/Putin_example.cxx and MURE/examples/Putin_example_sepent.cxx) shows basic geometrical methods.

* The 2nd and 3rd line of these files says that the resulting files will be either MCNP or Serpent input.

* Then a Connector plugin is built; this connector links MURE object with Serpent or MCNP ones. It is MANDATORY.

* One first builds 2 Materials (Graphite & Fuel).

* Then pure geometrical Shapes are defined: in this example, one uses only Sphere.

<footnote n="6">⁶Output files are mainly MCNP/Serpent input files named inp.example_name. All output files have been compressed; so you have to make a "gunzip *" inside MURE/examples/Output.</footnote>
* The “put in” operator (“>>”) is used to put 3 spheres like Russian dolls: *Smallest* and *Small* are called **Inside Shapes** of *Medium*.

* Then, this *Medium* "matrioshka" is duplicated with the `Shape::Clone()` method: the *Medium2* clone is an exact copy of *Medium*: it contains a small sphere which contains a smaller one.

* *Medium* and *Medium2* shapes are then translated (with their inside shapes) and then put in the *Big* sphere.

* Then, **Cells** are built (*Serpent/MCNP* cells): it associates a **Shape** with a **Material**. The exterior (*outBig* cell) is void and one has to specify a zero neutron importance because neutrons will not be followed (default importance is 1). For *Serpent*, this will correspond to the “outside” material. One has to build all **Cells** for all **Shapes**:

- for *Medium*, *Small* and *Smallest* shapes, it is easy because they have been explicitly declared

- *Medium2* being a clone of *Medium*, it has also 2 inside shapes: thus one has to define cells for these inside shapes: one access to `Shape::GetOriginalInsideShape()`; by convention, the most inner inside shape has the index 0, then the next one has the index 1, and so on.

* And finally, the *Serpent/MCNP* file is built; name is `minp` for MCNP input file and `sinp` for *Serpent* input file.

### 1.1.5.2 Fuel Evolution in a Sphere

The aim of this example<fnref n="7" /> of evolution (*MURE/examples/Evolution/EvolvingSphere.cxx* and *MURE/examples/Evolution/EvolvingSphere.cxx* and *MURE/examples/Evolution/EvolvingSphere.cxx*) is to show basic burn-up calculation in *MURE*. The example describes a 50cm radius sphere with an inner sphere (R=30cm) of metal uranium, a middle ring (thickness=10cm) of UOx fuel and an external water reflector (thickness=10cm). The water will not evolve with time; the thermal power of this “reactor” is 100 MW and it is kept constant for the whole evolution performed from *t* = 0 to *t* = 3 years.

* After defining the *Connector*<fnref n="8" />, the file begins with some general *MURE* settings: DATADIR, Temperature precision, Fission Product selection, ...

* Then 3 *Materials* are built: the water will not evolve whereas the metal and oxide uranium are evolving *Materials*.

* Then *Shapes* are defined (the 3 spheres and the exterior of the “reactor”).

* The *Cell* part associates *Shapes* and *Materials*.

* Then one defines an isotropic, mono-energetic neutron source (500 n/batch) to run in critical mode (KCODE) with 30 active cycles, 10 inactive cycles and an estimated initial $k_{eff} = 1$.

* The directory for the evolution is set (where *MC* files and evolution results will be written).

* One important function concerning reactor evolution is to be called here: it is the thermal power at which evolution is to be simulated. This value is set using the `gMURE->SetPower()` where the argument must be in watts. (c.f. Steady-state Power Normalization in § D.1).

<footnote n="7">⁷See comments in the file.</footnote>

<footnote n="8">⁸The “r” files of *MCNP*, in the *MCNP* version, will be removed (they are in general not useful and very large).</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 1.2: Main class structure of MURE. Main attributes are shown in blue. A Shape_ptr is only a typedef of Reference_ptr for a Shape. Idem for Nucleus_ptr.</figcaption></figure>

**NOTE:** For other applications or if the user already knows which neutron flux he wants to simulate, no power needs to be entered here but a tally normalization factor must then be specified. This **Tally Normalization Factor** can be set directly through the `gMURE->SetTallyNormalizationFactor()` method.

* Then one defines the 6 time steps at which a *MC* run is done, in a vector. The first *MC* step is always at time t=0 ; if this step is not defined, it will be insert in the vector. The evolution is performed up to the end but the last *MC* run is just the antepenultimate time step. This is important to know because this means the last values for *k*<sub>eff</sub>, fluxes, cross-sections are not *MCNP* values: for *k*<sub>eff</sub> the last value considered is the antepenultimate one.

## 1.2 MURE Package structure

The distribution package contains:

* The **PDF version** of this user guide.

* A complete **useful description of each class** (headers)

The general structure of main MURE's classes is shown in figure 1.2.

As it can be seen in figure 1.2, the link between classes is assured by the main class *MURE*. The *Shape_ptr* and *Nucleus_ptr* of the figure are smart pointer on *Shape* and *Nucleus* objects. The figure 1.3 shows the *Shape* inheritance graph using the 2 name spaces associated to the MC transport code used (MCNP or Serpent).
<figure><img src="image_5.png" /><figcaption>Figure 1.3: Shape inheritance scheme using the 2 namespaces.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 1.4: MureTally class.</figcaption></figure>

These name spaces are very important (see §3.1) ; one must use 1 (and only one) of them in any MURE input file.
Results from *MC* run (Serpent detectors or MCNP tallies) are stored in the MureTally class which contains MureBin
(cell, surface, ... bins) and TallyMultiplicator for reaction rates (see fig. 1.4 and 1.5).

1.3 MURE class

The *MURE* class is some kind of super class that handles connections between all the other classes. It controls all
flows: *MC* input file writing, nuclei trees building, evolution, ...

A global pointer (gMURE) on this class is defined to allow interaction between classes. For example, when the
geometry has been defined, as well as the *MC* source, tallies (or detectors), ..., the *MC* file will be written using

gMURE->BuildMCFile("myfile");

<figure><img src="image_7.png" /><figcaption>Figure 1.5: MureTallyBin class.</figcaption></figure>
this will generate an *MC* input file called "*myfile*" (if no argument is given to the previous method, the default input file name is "*inp*").

The name of *MC* exec command is given by

```
gMURE->SetMCExec("mcnp");
or
gMURE->SetMCExec("sss2");
```

By default, **MCNP::Connector** set the default value for the *MC* executable to "*mcnp5*" and **Serpent::Connector** set it to "*sss2*".

**MCNP Only:** You can specify the particle running mode by either `gMURE->SetModeN()` (the default) to run *MCNP* with neutron transport, `gMURE->SetModeNP()` to transport neutrons and photons and so on. This is not taken into account in the Serpent version.

Many methods exist and are grouped in sections to help the users; therefore examination of *MURE* class header is **strongly** recommended (see **Doxygen class description**).

## 1.4 MURE basic files

Different files are used or built in *MURE*. It is important to know what these files are to understand the behavior of *MURE*.

### 1.4.1 Files of "MURE/documentation"

This directory contains the most updated user guide and class documentation.

* the "*MURE/documentation/pdf*" directory contains the PDF files

* the "*MURE/documentation/html*" directory contains the HTML files, in particular an HTML version of the user guide (which is a translation of the PDF one). HTML links are easier to use as well as copy/paste for examples. Moreover the "*MURE/documentation/html/doxygen*" directory contains the automatic Doxygen description of classes.

### 1.4.2 Files in "MURE/data"

In the root *MURE* tree, there is a very important directory called "**data**" (known in *MURE* as the *DATADIR*); this directory can be changed by using either the environment variable **DATADIR**, or using the *MURE::SetDATADIR()* method<fnref n="9" />. In this directory, you will find:

* **chart.JEF3T**: this (old) file contains, for each nucleus, the half-life time and the decay modes. It is a "hand-made" file from nuclear data book and JEFF 2.0 library, from "*Nuclides and Isotopes*", [14] and "*Table of Isotopes*", [15].

* **chart.jeff3.1.1**: this file contains, for each nucleus, the half-life time and the decay modes. It is the updated version of the previous one. Data are extracted from NUBASE-2003, ENSDF, LNHB, UKHEDD-2.x and UKADD-6.x (see [16]).

<footnote n="9">⁹One can use either, `gMURE->SetDATADIR(path)` or the environment variable `DATADIR` (`setenv DATADIR path` in csh and `export DATADIR=path` in bash). Note that if you are using `gMURE->SetDATADIR(path)` you must call it before any Material definition.</footnote>
* **Mass.dat** and **NaturalIsotopeMass.dat** contain respectively the atomic mass for all nuclei and some natural compositions.

* **FPavailable.dat** and **FPyield.bin** are fission product yield files: the first one is an ASCII file containing the Z,A of fissiles, and the address of the fission yield (for the available energies) in the binary *FPyield.bin* file. These files are built from ENDF/B6 file. Source programs that build these files are in *MURE/utilson/fp*. These files are part of the MURE package and have been built for 32 and 64 bits computers ; if you want to use other FP yield, you need to rebuild these files using *MURE/utilson/fp/GenerateFPYield.cxx* (the compilation line is at the end of the file).

* **x DIRPreQU.dat** is the "header" of a general *x DIR* (from the top to the "directory" keyword). It is used in the automatic XDIR construction for *MCNP*.

* **BaseSummary.dat** is the file that contains all available nuclei, their temperature and the corresponding *x DIR* line.

  It is built from *MURE/utilson/datadir/ExtractTree.cxx* and/or *MURE/utilson/datadir/ExtractXsdir.cxx*. The first one extracts information from a directory tree and the second from an existing *x DIR*. This file is the only way for *MURE* to use automatic MC nucleus code. The *xsdata* file for *Serpent* is also built from this file, using a conversion technique close to the perl script of *Serpent*.

* **AvailableReactionChart.dat** contains, for all nuclei of the **chart.jeff3.1.1** (or the old **chart.JEF3T**) and **BaseSummary.dat**, whether or not reactions are available for **MC** codes. The aim of this file is to speed up the NucleiTree construction. It is built from *MURE/utilson/datadir/CheckReaction.cxx*.

* **IsomerProduction.dat** contains necessary information to allow the $(n, \gamma)$ isomer productions for some important nuclei such as $^{241}$Am, $^{109}$Ag, ...

**NOTE:** When you add new nucleus or temperature in nuclear cross-sections files, these modifications will be taken into account ONLY if *BaseSummary.dat* and *AvailableReactionChart.dat* have been updated. You also need to remove the local *ReactionList* directory (see after § 1.4.6).

### 1.4.3 Files of "MURE/utilson"

This directory contains codes that can be very helpful to users.

* "endf2ace" directory contains a small interface to **NJOY** (see documentation in the directory)

* "fp" contains a utility to rebuild **FPyield.bin** and **FPavailable.dat** of *MURE/data* directory. Main file is *GenerateFPYield.cxx* which extracts from an ASCII ENDF file (containing the FP yields) the fission product yield in a format used by MURE.

* "datadir" contains programs that allows users to build or complete the *BaseSummary.dat* and *AvailableReactionChart.dat*.

- *ExtractXsdir.cxx* allows users to build/complete a *BaseSummary.dat* from an existing *MCNP* *x DIR*.

- *ExtractTree.cxx* allows users to build/complete a *BaseSummary.dat* from a directory tree of cross-sections ; this tree can be built by the "endf2ace" utility. The form of the directory tree is
base_name/base_version/Z/AAA/Isomeric_state/temperature/bin/

where, for example, base_name=ENDFB, base_version=6.8, Z=92, AAA=235, Isomeric_state=0 (for ground state), temperature=600 (in K). In the directory "bin" one find the binary cross-section of the desired nucleus in the file "ace" and an xsdir line for this nucleus in a file named "dir" (for example the line looks like '92235.10c 233.025000 ace 0 2 1 621851 4096 512 5.1704E-08 ptable')

In both cases, the compilation line is at the end of the source code.

NOTE: If an existing *BaseSummary.dat* is already present in *MURE/utilsexdir*, this file is completed ; otherwise it is created.

- *CheckReaction.cxx* must be used after building/modifying the *BaseSummary.dat* file in order to build the *AvailableReactionChart.dat*.

NOTE: After running *ExtractXsdir/ExtractTree* and *CheckReaction*, you have to copy the *BaseSummary.dat* and *AvailableReactionChart.dat* in your *MURE/data*.

### 1.4.4 Graphical User Interface “MURE/gui”

The GUI of MURE is located in the *MURE/gui* directory. The executable name is “MureGui” and it reads the result of a *MURE* or *Dragon* evolution. To use it, user must install *ROOT* (free download at https://root.cern.ch). Normally, the *install.sh* script of *MURE*, defined all available libraries in the *MURE/config* directory. If the Lapack package is needed, (see http://www.netlib.org/lapack/index.html)<fnref n="10" />. The compilation is obtained by a "make" in the *MURE/gui* directory. Using *MureGui* without argument gives a short description of the code used (see also § 8.3).

### 1.4.5 MURE Source files “MURE/source/include” and “MURE/source/src”

The MURE source files are located in "MURE/source/include" for the header files (`.hxx`) and in "MURE/source/src" for implementation (`.cxx`). When user is modifying the MURE source, he has to recompile the MURE package (by a "make" in *MURE/source/src*). In *MURE/source/external*, one finds 2 auxiliary libraries, **ValErr** and **Mctal** that are used in MURE. The **ValErr** defines a class to handle numbers and their errors, and **Mctal** is dedicated to reading/writing MCNP "m" files.

In any program using MURE, you have to include at least:

```c
#include <MureHeaders.hxx>
```

and either one of the following

```c
#include <MureMCNPHeaders.hxx>
```

or

```c
#include <MureSerpentHeaders.hxx>
```

If needed on can add also :

```c
#include <libmctal/TMTally.hxx>
```

```c
#include <libmctal/TMComment.hxx>
```

```c
#include <libmctal/TMctal.hxx>
```

<footnote n="10">¹⁰Most of LINUX distributions provide LAPACK package ; a 'locate liblapack.so' shows if it is already installed ; if only liblapack.so.version_number is found, you either have to install the "lapack-devel" package or to make a 'ln -s /path_to_lapacklib/liblapack.so.version liblapack.so' in the *MURE/lib* directory.</footnote>
1.4.6 Other files

Each time a new evolution is run, a directory **ReactionList** is locally built (if it is not already there) with the available
reaction of user's nuclear data base (one binary file for each Z). A list of suppressed reactions for some of the nuclei
(because they are lower than a given threshold) is also written in **ReactionList/SuppressReaction.dat**. The aim
of these files is to save time when **using the same nuclear data base AND the same reaction threshold,**
**life-time cutting** and so on. Thus if you modify one of these, **DO NOT FORGET** to remove the **ReactionList**
directory.
# Chapter 2

## What's new in MURE

1. November 2016:

(a) Add Z revolution axis Torus for MCNP and Serpent; see *MathZTorus*, *MCNPZTorus* and *SerpentZTorus* classes.

(b) Update "chart.jef3T" file that contains periods and decay modes for all the chart; this file (JEFF 2 lib.) is kept unchanged, but a new file has been created "chart.jeff3.1.1" (from ref. [16]); it is the new default.

2. November 2015:

(a) *ReactorMesh* and *ReactorChannel* class have been removed; a new *GenericReactorAssembly*, and its concrete versions *MCNP::ReactorAssembly* and *Serpent::ReactorAssembly* replace the *ReactorMesh* class. It is a complete rewriting of that class, but the spirit is kept (see § 9.1 and 9.3.2). It adds lot of new possibilities (add duct to an assembly, filling core with a *ReactorAssembly*, ...).

(b) *COBRA* class is renamed *COBRA_EN* and rewritten; again backward compatibility is broken (it does not inherit from *GenericReactorAssembly*) but the spirit is not changed.

3. October 2015:

(a) Correct a bug when reading the "standard" `xsdir` of MCNP5; isomeric cross-section libraries are now in accordance with the MCNP5 documentation.

4. July 2015:

* **New MURE version 2.0 (SMURE).** This is a major release of MURE. It couples now both MCNP & Serpent. Due to technical, philosophical difficulty and also to beauty in programming (dry codes), backward compatibility is broken. But the changes to user input files between MURE v1 and MURE v2 are small. See 3.

5. November 2014:

(a) Add the possibility to use standard MCNP tallies $^{238}U$ and multi-group tallies for other nuclei (see 5).

6. October 2014:
(a) Add an optional equilibrium treatment for $^{135}$Xe (see 7.7)

7. April 2013:

(a) Version of base and isomeric states (metastable and ground state) for $^{242}$Am in *BaseSummary.dat* file for standard library (distributed with **MCNP** or available@NEA) are now more conform to the reality; you have to regenerate your *BaseSummary.dat/AvailableReaction.dat* files with the exec of *MURE/utilsexadadir*. **DON'T FORGET TO REMOVE YOUR ReactionList DIRECTORY IN ORDER TO TAKE INTO AC- COUNT THESE MODIFICATION.**

(b) *BasePriority* has been debug and it seems to work as expected...

8. March 2013:

(a) Add a generalization of the *MURE::SetMode()* method to allow any type of particle transport (mainly for *MCNPX*)

(b) Cell importance may have different value for each transported particles (by successive call to *Cell::AddParticle* method)

(c) extend the *MCNPSource* possibilities (particle distributions are now allowed for *MCNPX*)

9. July 2012:

(a) Generalization of isomer production by $(n, \gamma)$ for special cases (such as $^{241}$Am, $^{109}$Ag, ...), and bug correction in $(242\text{Cm}, 242\text{Pu})$ production, see 6.1.6.

(b) Improve **MurGui** Spectrum radiotoxicity tab

(c) Allow to plot **neutron balance** in the "Reaction Rate" tab.

10. June 2012:

(a) Add functions in **MureGui** (see MureGui's Radiotoxicity tab)

i. Nuclei extraction is now possible after a given cooling time.

ii. $\gamma,\beta, \alpha$ and neutron spectra for evolving materials can be computed.

(b) New **MURE** install procedure (more simple, more robust): **any fresh install or update NEED to run first the install.sh script** (require bash shell).

(c) Add the possibility to used multi-threading (*OpenMP*) during the evolution and for updating multigroup $\sigma\phi$ when the *gcc* version support this option (see *INSTALL* procedure and *install.sh* script and section 7.2.2.2).

(d) Add some new kind of *MCNPSource* (see section 5.2.3)

(e) Add fluence to dose conversion for tallies

11. May 2011:

(a) Add a new class *GammaSpectrum* class (section 8.4.2)
(b) New option of *MureGui*: use the "radiotoxicity tab" of *MureGui* on a dumped Material created with *MURE* (see 8.3.6 and "MURE/example/GammaSpectrumExample.cxx")

(c) New methods `MCNPSource::SetAXS` and `MCNPSource::SetVec` allows user to define collimated sources.

(d) In `Tally::Tally(int type,const char *particle)`: if `type < 0`, it's change the Tally units.

i. For example `Tally *t=new Tally(-4,"P")` define a tally in `MeV/cm²` rather than `Particles/cm²`.

(e) A new method `MURE::SetModeP()` can be used to allow only the photon transport.

12. July 2009:

(a) Improve the English in User Guide: **many thanks to Erica Agostinho for this painful work.**

(b) Implementation of *PseudoMaterial*: in order to take into account temperature effects, one can process the nuclear data at the desired temperature (using *NJOY*) or use an interpolation between 2 already existing temperatures; this later method is used in the *PseudoMaterial* techniques.

13. Avril 2009: MURE is available at NEA DataBank

14. January 2009:

(a) Improve reading time of *BaseSummary.dat* file: it is now greatly recommended that this file is ordered by Z,A,I (this is the case if it is generated by *ExtractedXsdir.cxx* and *ExtractTree.cxx*).

(b) Rewrite long parts of the documentation (*User Guide* and *HTML* class description)

(c) Modify examples directory: it now contains documented examples and output.

15. November 2008: Improve radiotoxicity post treatment in *MureGui*.

16. September 2008: Implement "multigroup" calculation for reaction rates: a very narrow energy binning for flux calculation is put in each *MCNP* run; then reaction rates are obtained by reading *ACE MCNP* files after each *MCNP* run. The result of such method saves a considerable amount of CPU time for *MCNP* (at least a factor 30) with only a low percentage of discrepancy (~1 to 3%) in result compared to the standard calculation (reaction rates are tallied in *MCNP*).

17. April 2008:

(a) implementation of Predictor-Corrector method in the evolution

(b) possibility to read ASCII nuclear data file (ACE format): this avoids problems due to binary compatibility (size of real (float or double?), little or big endian, ...). BUT it is much longer to (1) build the *ReactionList* directory and (2) run a *MCNP*.

18. December 2007:

(a) Disable the `σϕ` extrapolation by default. It has been shown, but not really understood, that this treatment introduces a larger dispersion in the result, after N identical evolutions, than doing nothing (no `σϕ extrapolation`).

(b) Modify the evolution using a *MCNP* User Input geometry file:
i. "like but" cells can be used (but not evolved)

ii. "MCNP Transformation" cards can be used

iii. the 3rd block of the MCNP file is read and copied except the materials (that must be defined in the *MURE* file). Thus the source as well as all other cards of this block can be used without defining them in the *MURE* file (this is also true for user defined tallies).

19. October 2007:

(a) Switch from `ccdoc` to `Doxygen` for class documentation

(b) **Change completely the MURE directory trees**

(c) **NON BACKWARD COMPATIBILITY**

i. Material definition has changed: now only 2 constructors must be used:

* `Material()`: for standard material (you have to specify density, ... with the Set methods)

* `Material(int)`: for using materials from an MCNP input file geometry.

ii. THE COPY CONSTRUCTOR HAS NOT TO BE CALL: CALL only the `Material::Copy()`.

iii. The `Material::Mix` has been modified (number of arguments, units required). see `Material.hxx`

iv. The `Material::AddNucleus` has been modified (number of arguments). By default the proportion units are "kpMOL" (i.e. molar proportion). But the unit must be specified if you use a moderator (MT card of MCNP).

v. The proportion units (both for Density and Proportion) must be used for any `Material::GetProportion()` and `Material::GetDensity()` methods

* for Proportion the only valid units are: `kpMOL` (molar prop), `kpMASS` (mass prop), `kpATCM3` (at/cm3), `kpATBCM` (at/barn.cm)

* for density the only valid units are: `kdGCM3` (g/cm3), `kdATCM3` (at/cm3), `kdATBCM` (at/barn.cm)

vi. A new material has been defined : `ControlMaterial` (public of Material). This class is used for Poison, Fissile or other control of reactivity (e.g. `poison.cxx` in *MURE/examples*)

vii. All `Print()` methods now return a **string** instead of a void: to use them: `cout<<Mat->Print()`<<endl;` for example.

viii. `EvolutionControl` class has been cleaned (as well as *MURE* class). If you want to use special control you have to write your own derivative class. Examples using `PoisonControl`, `FissileControl` & `HNControl` (but Adrien you have to rewrite them and look carefully at `TMSR.cxx` in *MURE/examples*.) and Rod control are defined in `source/src`. You can use them as they are or defined your own using these examples.

(d) Almost all `cout/cerr` have been removed from classes; used instead `LOG_DEBUG`, `LOG_INFO`, ... ; `LOG_INFO` is now independent of `LogLevelMessage`; it is always printed. If `MURE::SetMessageLevel` is set to `LOG_LEVEL_DEBUG`, all `LOG_DEBUG` are printed. But if `MURE::SetSilentDebug` is used, only the `LOG_DEBUG`s of methods where a `"int DODEBUG=1"` is inserted are printed. Thus, using `LOG_DEBUG`, avoids to comment all `"cout"` when no debugging is desired.

(e) ALL EXAMPLES HAVE BEEN UPDATED TO TAKE INTO ACCOUNT THESE MODIFICATIONS: please READ THEM!!!!!!!!!
20. September 2007:

(a) Correct an important bug in evolution using an *MCNP* User Input geometry file: number of Materials were not correct (Thanks to Jan Frybort).

21. Juin 2007:

(a) Rename the MURE header file `Shapes.hxx` into `MureHeaders.hxx`: this is more logical...but you have to change your MURE files....

(b) Add a new class `EvolutionWrapper` to simplify and extend `EvolutionControl` capabilities

(c) Suppress the writing of `BDATA_ xxx` and `DATA_ xxx` ; now, by default, only `BDATA_ xxx` are written.
This can be changed using the `MURE::SetWriteBinaryData()` and `MURE::SetWriteASCIIData()` methods.

(d) One can start an evolution from a given step : suppose that the evolution stops at step *i* ; an evolution can be started from the step *i+1* using `MURE::Evolution(T, i+1)`. Warning: it is probably not correct for `OutCoreEvolutiveSystemVector`...you must do the evolution from the first step as before.
# Chapter 3

## From MURE to SMURE: Choosing the Monte-Carlo Transport Code

As already said, coupling *MURE* with a other *MC* code had 2 consequences:

* decouple *MURE* from *MCNP*, which was a big job because they were closely linked.

* couple *MURE* to a *MC* code (*Serpent* or *MCNP*), which is also a big job because the *Serpent* an *MCNP* philosophy are sometimes completely different.

The aim is with a minor change to input file, it will be able to generate output for *MCNP* or for *Serpent 2*. In this chapter, it is explained first how to switch from an *MCNP* input file to a *Serpent* input file (or vice-versa) and then how to migrate from *MURE 1.0* to *SMURE (MURE 2.0)*.

### 3.1 Switching from MCNP to Serpent in a MURE input file

*SMURE* has been designed to make the change from *Serpent* input to *MCNP* input (or vice-versa) as simple as possible. Of course, *Serpent* is not *MCNP* and conversely. Thus each of these *MC* codes have their own possibility, not implemented in the other. Moreover, the maturity of the *MURE* coupling is much greater (few years of work for *MCNP* versus 4 months of work for *Serpent*). Thus this section deals only with things that are possible with both *MC* codes and has been implemented in *MURE*. Nevertheless, it seems to us that it is enough to make a standard reactor study (including evolution).

#### 3.1.1 Principle of the implementation

In *MURE*, the choice of *MC* output is performed via a **ConnectorPlugin**; this class has 2 sibling, a **MCNP::Connector** and a **Serpent::Connector** (see Fig. 3.1); thus the main trick is the use of a specific namespace (*MCNP* or *Serpent*). Geometrical shapes are defined in *MURE* as **MathShape** (i.e. a “pure” mathematical Shape, e.g. a **MathSphere**). Then, in each namespace, are defined concrete (or real) shapes that essentially provide a *print* method derivating from the **MathShapes** (see Fig. 1.3). For example, a **Serpent::Sphere** and **MCNP::Sphere** are daughters of **MathSphere** and implement a **Serpent::Sphere::Print()** or a **MCNP::Sphere::Print()** that will produce respectively the right output for the corresponding MC code. The same trick applies for tallies (or detectors in *Serpent*).
<figure><img src="image_8.png" /><figcaption>Figure 3.1: The ConnectorPlugin class.</figcaption></figure>

Thus, using a **Connector**, a **Sphere**, ... in *MURE* is just “resolved” by the namespace (that should not be forgotten). In order it works without problems, appropriate header files must be included. First header file, is the “pure” *MURE* header : *MureHeaders.hxx* which includes all pure *MURE* classes. It is thus included in all *MURE* input files, whatever will be the output. The second header file is linked to the output ; *MureMCNPHeaders.hxx* and *MureSerpentHeaders.hxx* includes all classes that couple *MURE* with *MCNP* and *Serpent* respectively. They include particularly **MCNP::Connector** (*Serpent::Connector*), **MCNP::Sphere** (*Serpent::Sphere*), ...

Of course connectors, because of some profound differences between *MCNP* & *Serpent*, have specific methods existing only in one type of connector. The same applies also for the **Tally** class which is more complete in *MCNP*. The last case is for particle sources which are very different from *MCNP* to *Serpent*. This is why, sibling of **MCSource** class keep the *MC* code prefix in there name (*MCNPSource* and *SerpentSource*) in order that the user don't forget its specificity. Another cause of problems, switching from one code to the other is the **MURE::AddSpecialCard()** which allows to use specific *Serpent* or *MCNP* cards that have not been implemented in *MURE*. Thus a careful check of **Connector**, **MCSource** and **MURE::AddSpecialCard()** should be done when switching from one code to the other.

### 3.1.2 A example

Switching from one *MC* code to the other can be performed very fast : edit the file (let say a *MCNP* version) and replace all **“MCNP”** string by **“Serpent”** string. This should change 4 occurrences: the header, the namespace and the source name (See example of Tab. 3.1). Then remove specific command for *MCNP* and add specific command for *Serpent*. To illustrate this, let us consider the “evolving sphere” example. The “red” parts of Table 3.1 correspond to a “replace all” “*MCNP*” string by “*Serpent*” string. The magenta words are “output” directory name, that here, are different to keep a trace of both evolution. The blue lines is an example of a pure specific method of one code (here *MCNP*) of one of the 3 classes **Connector**, (*MCNP/Serpent*)Source or **Tally**. And then, in green, a **MURE::AddSpecialCard()** that in this case is only for *Serpent*.

As Shown on this example, the number of line to modify is very limited ; moreover, dealing with a more complex case, this number will not increase significantly (probably less than a factor 2).

## 3.2 How migrate my old MURE V1.x file to the MURE V2.x file

### 3.2.1 What has definitely change

In this part, only the coupling with *MCNP* is considered (as in the old *MURE* v1). Main changes are due to

* Header file has change (row 1 of the Tab. 3.2) ; one must use a namespace *MCNP*.

* One has to define a **ConnectorPlugin** of type **MCNP::Connector** (see row 1 of Tab. 3.2) just after the main declaration. The connector is responsible for all pure *MCNP* methods. Row 3 of Tab. 3.2 give an ex-
<caption>Table 3.1: Switching from MCNP (left) to Serpent (right) output files.</caption>

```cpp
#include <>iostream>
using namespace std;
#include <MureHeaders.hxx>
#include <MureMCNPHeaders.hxx>
using namespace MCNP;
int main(int argc, char **argv)
{
    Connector *plugin=new Connector();
    //remove the MCNP "r" file after each evo step
    plugin->SetRemove_r_files();
    gMURE->SetConnectorPlugin coupler);
    // Nuclei tree construction & nuclear data library
    //
    gMURE->SetMessageLevel(3);
    gMURE->SetDATADIR("../../data/");
    gMURE->SetAutoXSDIR();
    gMURE->SetShortestHalfLife(3600);
    gMURE->GetTemperatureMap()->SetDeltaTPrecision(1500);
    gMURE->KeepOnlyFissionProductSelection();
    //
    // Materials
    //
    Material *H20= new Material();
    H20->SetDensity(1.0); // density in g/cm3
    H20->AddNucleus(1, 1, 0, 2.); // H of H20
    H20->AddNucleus(8, 16, 0, 1.); // O of H20
    Material *Uranium_metal= new Material();
    Uranium_metal->SetDensity(19.0);
    Uranium_metal->SetTemperature(600.0); //Temperature in K
    Uranium_metal->AddNucleus(92,235,0.20); //U-235 in mol
    Uranium_metal->AddNucleus(92,238,0.80); //U-238 in mol
    Uranium_metal->SetEvolution(); //this is an evolving material
    //
    Material *UOx= new Material();
    UOx->SetDensity(10.2);
    UOx->AddNucleus(92,235,0, 0.05); // 5%(in mol) of U-235
    UOx->AddNucleus(92,238,0, 0.95); // 95%(in mol) of U-238
    UOx->AddNucleus(8,16,0, 2.); // 2 O for 1 UO2
    UOx->SetEvolution(); //this is an evolving material
    //
    // Shapes
    //
    Shape_ptr Shape_MySphere(new Sphere(0.5));
    Shape_ptr Shape_Exterior(!Shape_MySphere);
    Shape_ptr Shape_MySphere2(new Sphere(0.4));
    Shape_ptr Shape_MySphere3(new Sphere(0.3));
    Shape_MySphere3>>Shape_MySphere2>>Shape_MySphere;
    //
    // Cells
    //
    Cell* CellSphere = new Cell(Shape_MySphere, H20);
    Cell* CellSphere2 = new Cell(Shape_MySphere2, UOx);
    Cell* CellSphere3 = new Cell(Shape_MySphere3, Uranium_metal);
    Cell* CellExterior = new Cell(Shape_Exterior, 0, 0);
    //
    // General Cards
    //
    //define the neutron source (KCODE)
    MCNPSource *s=new MCNPSource(500);
    s->SetKcode(60,30,1.);
    s->SetEnergy(2e6);
    gMURE->SetSource(s);
    // MC option
    gMURE->SetComment("Evolution of a Sphere");
    gMURE->SetMCInputFileName("inp");
    gMURE->SetMCRunDirectory("sphere","keep");
    //
    // Evolution
    //
    gMURE->SetPower(1e8); //give the total power: 100 MW
double TEnd=3*365.25*3600.*24;
vector double T;
int Nt=6;
double dt=TEnd/Nt;
for(int i=0; i<Nt;i++)
    T.push_back(i*dt);
gMURE->Evolution(T);
}
```

```cpp
#include <>iostream>
using namespace std;
#include <MureHeaders.hxx>
#include <MureSerpentHeaders.hxx>
using namespace Serpent;
int main(int argc, char **argv)
{
    Connector *plugin=new Connector();
    gMURE->SetConnectorPlugin coupler);
    // Nuclei tree construction & nuclear data library
    //
    gMURE->SetMessageLevel(3);
    gMURE->SetDATADIR("../../data/");
    gMURE->SetAutoXSDIR();
    gMURE->SetShortestHalfLife(3600);
    gMURE->GetTemperatureMap()->SetDeltaTPrecision(1500);
    gMURE->KeepOnlyFissionProductSelection();
    //
    // Materials
    //
    Material *H20= new Material();
    H20->SetDensity(1.0); // density in g/cm3
    H20->AddNucleus(1, 1, 0, 2.); // H of H20
    H20->AddNucleus(8, 16, 0, 1.); // O of H20
    Material *Uranium_metal= new Material();
    Uranium_metal->SetDensity(19.0);
    Uranium_metal->SetTemperature(600.0); //Temperature in K
    Uranium_metal->AddNucleus(92,235,0.20); //U-235 in mol
    Uranium_metal->AddNucleus(92,238,0.80); //U-238 in mol
    Uranium_metal->SetEvolution(); //this is an evolving material
    //
    Material *UOx= new Material();
    UOx->SetDensity(10.2);
    UOx->AddNucleus(92,235,0, 0.05); // 5%(in mol) of U-235
    UOx->AddNucleus(92,238,0, 0.95); // 95%(in mol) of U-238
    UOx->AddNucleus(8,16,0, 2.); // 2 O for 1 UO2
    UOx->SetEvolution(); //this is an evolving material
    //
    // Shapes
    //
    Shape_ptr Shape_MySphere(new Sphere(0.5));
    Shape_ptr Shape_Exterior(!Shape_MySphere);
    Shape_ptr Shape_MySphere2(new Sphere(0.4));
    Shape_ptr Shape_MySphere3(new Sphere(0.3));
    Shape_MySphere3>>Shape_MySphere2>>Shape_MySphere;
    //
    // Cells
    //
    Cell* CellSphere = new Cell(Shape_MySphere, H20);
    Cell* CellSphere2 = new Cell(Shape_MySphere2, UOx);
    Cell* CellSphere3 = new Cell(Shape_MySphere3, Uranium_metal);
    Cell* CellExterior = new Cell(Shape_Exterior, 0, 0);
    //
    // General Cards
    //
    //define the neutron source (KCODE)
    SerpentSource *s=new SerpentSource(500);
    s->SetKcode(60,30,1.);
    s->SetEnergy(2e6);
    gMURE->SetSource(s);
    // MC option
    gMURE->SetComment("Evolution of a Sphere");
    gMURE->SetMCInputFileName("inp");
    gMURE->SetMCRunDirectory("sphere_serpent","keep");
    gMURE->AddSpecialCard("plot 3 2000 2000");
    //
    // Evolution
    //
    gMURE->SetPower(1e8); //give the total power: 100 MW
double TEnd=3*365.25*3600.*24;
vector double T;
int Nt=6;
double dt=TEnd/Nt;
for(int i=0; i<Nt;i++)
    T.push_back(i*dt);
gMURE->Evolution(T);
}
```
<table><caption>Table 3.2: Summary of changes from MURE to SMURE.</caption><thead><tr><th>MURE version 1</th><th>MURE version 2 : SMURE</th></tr></thead><tbody><tr><td><pre>#include &lt;MureHeaders.hxx&gt;</pre></td><td><pre>#include "MureHeaders.hxx"<br>#include "MureMCNPHeaders.hxx"<br>using namespace MCNP;</pre></td></tr><tr><td><pre>int main(...)<br>{</pre></td><td><pre>int main(...)<br>{<br> Connector *plugin=new Connector();<br> gMURE-&gt;SetConnectorPlugin plu<code>gMURE-&gt;SetRemove_r_files();</code></td></tr><tr><td><pre>gMURE-&gt;BuildMCNPFile();</pre></td><td><pre>plugin-&gt;SetRemove_r_files();</pre></td></tr><tr><td><pre>gMURE-&gt;SetMCNPInputFileName("inp");</pre></td><td><pre>gMURE-&gt;BuildMCFile();</pre></td></tr><tr><td><pre>gMURE-&gt;SetMCNPRunDirectory("sphere","keep")</pre></td><td><pre>gMURE-&gt;SetMCInputFileName("inp");</pre></td></tr><tr><td><pre>gMURE-&gt;SetMCNPExec("mcnp5");</pre></td><td><pre>gMURE-&gt;SetMCRunDirectory("sphere","keep")</pre></td></tr><tr><td><pre>gMURE-&gt;SetWriteASCIIData();</pre></td><td><pre>gMURE-&gt;SetMCExec("mcnp5");</pre></td></tr><tr><td><pre>gMURE-&gt;SetWriteBinaryData();</pre></td><td><pre>gMURE-&gt;GetEvolutionSolver()-&gt;SetWriteASCIIData();</pre></td></tr><tr><td><pre>gMURE-&gt;UseMultiGroupTallies();</pre></td><td><pre>gMURE-&gt;GetEvolutionSolver()-&gt;SetWriteBinaryData();</pre></td></tr><tr><td><pre>gMURE-&gt;SetUseNewDBCN();</pre></td><td><pre>gMURE-&gt;GetEvolutionSolver()-&gt;UseMultiGroupTallies();</pre></td></tr><tr><td><pre>Cell *MyLattice=new Cell(AGenerator);<br>MyLattice-&gt;Lattice(1,-rx,rx,-ry,ry,-rz,rz);</pre></td><td><pre>LatticeCell *MyLattice=new LatticeCell(AGenerator);<br>MyLattice-&gt;SetLatticeRange(-rx,rx,-ry,ry,-rz,rz);</pre></td></tr></tbody></table>

ample for removing MCNP “r” files; but the same apply for `SetPRDMP()`, `SetMCNP4B()` which are now in `MCNP::Connector` and no more in `MURE`.

* Method names of `MURE` containing “MCNP” are (in general) now named with “MC” ; examples are given from row 4 to 7.

* Evolution is no more in `MURE` class but in the `EvolutionSolver` class ; example of change are given in row 8 to 10. Almost all `MURE` evolution related method are now in `EvolutionSolver` excepted `MURE;SetPower()` and `MURE::Evolution()`.

* Lattice implementation has changed :

- Replace the declaration from a simple `Cell` to a `LatticeCell` (row 11)

- the `Cell::Lattice()` is replaced by `LatticeCell::SetLatticeRange()` and the lattice type is no more necessary (see row 12)

* `ReactorChannel` & `ReactorMesh` classes: the first one is, since a while, obsolete ; it is removed in MURE V2.0. The second one has been deeply rewritten and rename as `ReactorAssembly` but the spirit has not change ; most of work (and thus time spent to migrate from MURE V1.x to MURE V2.0 will be in the change of this `ReactorMesh` to `ReactorAssembly` class but we think that the improvements of `ReactorAssembly` justify this “extra-work”. For more details on the `ReactorAssembly` class see section 9.3.2.

The time spent to change a very big file from version 1 to version 2 is around 5 minutes maximum.

### 3.2.2 What is still existing but can be done in a more elegant way

In most of real case fro reactor physics, one has to use lattice, fuel pin, ... A new improvement of *MURE version 2* is to implement a **PinCell** which in fact copy the elegant *Serpent's* pin definition. To illustrate the advantage of such
a **PinCell** class, the main part of the lattice example in the “old” way is presented and then rewrite using **PinCell**.
For sake of clarity, only *Shape* and *Cell* part are presented. All sizes have been defined, as well as 3 Materials : H2O,
Iron, and UOx. As shown on Tab. 3.3, the use of **PinCell** makes the code shorter (the 8 blue lines of the shape
declaration have been suppressed), and easier to define and to read. The **PinCell** declaration make cell part slightly
longer (4 "magenta" lines versus 6 "green" lines). To be noticed, the use of a **PinCell** without a layer allows to define
the whole universe. And a last comment on the *FillLattice* methods : in the first version, a Shape is used as argument
(red variables on the left) where as in the new version this is a **PinCell** (red variables on the right).

Thus user is strongly encourage to use these PinCells.
# Chapter 4

## Geometry Definition

### 4.1 Introduction

C++ logic for building geometries is slightly different for the *MCNP* or *Serpent* one; therefore, each time a new geometry is built you should check it with *MCNP/Serpent* before using it.

There are two base classes to build a geometry: *Shape* and *Cell*. *Shape* describes only geometrical shapes, and *Cell* corresponds to an *MC* code cell (i.e., it has a material, importance, etc.).

Shape objects correspond to simple geometrical shapes (sphere, plane, ...) as well as more complex ones resulting from the intersection and/or union of simple shapes (Intersections/unions are defined by the *Node* class). A *Node* is a “tree” of intersections/unions of Shapes. For fast calculations, a node tree has to be as simple as possible. Special methods are available for simplifying the node trees which can (in general) determine whether or not a Shape is disjointed (or included) of (in) another Shape (see example in figures B.1 to B.4 in Appendix B).

These simplifications may result in the deletion of some Shapes. But, because one must not destroy a Shape in a Node if that Shape belongs to another Node, a special way to handle Shape creation/destruction has been implemented (via *Reference_ptr*, which is a C++ smart pointer).

**In conclusion: the user must only use *Shape_ptr* (a “reference Shape”) and not Shape. Shape_ptr is a pointer to Shape with Reference_ptr.** In that way, the *Shape_ptr* will be destructed only if it is no longer referenced, otherwise, its “deletion” leads to a decrement of the number of references.

### 4.2 Definition of geometrical shapes

* There are 4 available base *Shapes* :

- *Plane* (infinite),

- *Cylinder* (infinite),

- *Sphere*,

- *Brick* (finite or infinite)

* Then one can define *Node* (unions or intersections of Shapes) with 2 already defined Nodes:
- Tube

- Hexagon (finite or infinite)

A *Brick* is a rectangular parallelepiped (and can be infinite). A *Tube* is a finite cylinder with an optional inner radius that defines a “tube”. *Hexagons* can have finite or infinite height. In general, a user will only need to use *Spheres*, *Bricks*, *Tubes* and *Hexagons*.

One can define the interior or the exterior of each *Shape*. The complement of a *Shape* can be defined by the `Not()` method as well as by the `!` operator (see examples given later on).

Other shapes available in MCNP(X) or Serpent may be added by any user. This is not too hard even if it requires some work. The best way to do it, is to read the implementation of existing shapes to have an idea.

### 4.2.1 Units

**WARNING:** In MURE, the length unit is the meter (whereas in MCNP and Serpent it is the `cm`). In particular, volumes in MURE are in $m^3$.

<table><tr><td>Length</td><td>m</td></tr><tr><td>Energy</td><td>eV</td></tr><tr><td>Temperature</td><td>K</td></tr><tr><td>Density</td><td>g/cm³, atom/cm³ or atom/(barn.cm)</td></tr><tr><td>Proportion</td><td>%mol, %mass, atom/cm³ or atom/(barn.cm)</td></tr></table>

### 4.2.2 Examples of simple shapes

To define the interior of a origin center sphere of radius $R$:

```
Shape_ptr S(new Sphere(R));
```

To define the outer part of $S$ one can do

```
Shape_ptr Ext_S(!S);
```

or

```
Shape_ptr Ext_S(S->Not());
```

or

```
Shape_ptr Ext_S(new Sphere(R,0,0,0,1));
```

where 3 zeros correspond to the sphere center (the origin) and +1 to the exterior of the *Sphere* (default=-1).

For *Hexagons* and *Bricks*, two versions exist: finite or infinite *Shapes*. For finite bricks and hexagons, you should define it as:

```
Shape_ptr B(new Brick(HalfX, HalfY, HalfZ, Signe));
```

```
Shape_ptr H(new Hexagon(HalfHeight, Side, Signe));
```

where `HalfX` (resp. `Y` and `Z`) is the half length of the brick in the `X` (resp. `Y` and `Z`) direction, `HalfHeight` the half height (!) of the hexagon, `Side` its side, and `Signe` the sign defining whether it is the inner or the outer part, just as for the sphere. For infinite ones, to avoid conflicts between the definitions, a string at the beginning is necessary (and, of course, `HalfHeight` and `HalfZ` are irrelevant):
Shape_ptr B(new Brick("any string you want", HalfX, HalfY, Signe));
Shape_ptr H(new Hexagon("any string you want", Side, Signe));

The infinite shapes obtained are parallel to the Z axis, but they may be rotated afterwards.

### 4.2.3 Examples of simple Nodes

A Node is the Union (+1) or the Intersection (-1) of Shapes (i.e., simple Shapes or Nodes).

To define the intersection of a sphere centered at $(x_0, y_0, z_0)$ of radius $R$ with a square brick of side $a$ centered at the origin, one may use:

Shape_ptr S(new Sphere(R,x0,y0,z0));
Shape_ptr B(new Brick(a/2,a/2,a/2));
Shape_ptr Inter=S & B;

whereas the union of the sphere and the brick may be defined as

Shape_ptr Union=S | B;

### 4.2.4 Moving a Shape

It is possible to move (translation/rotation) a Shape (or a Node) via the `Shape::Translate` and `Shape::Rotate` methods.
For example, to translate the `Shape_ptr` B of $(dx,dy,dz)$,

B->Translate(dx,dy,dz);

and to rotate clockwise B around $(x_0, y_0, z_0)$ of $\phi, \theta$ and $\psi$ around the $z, y$ and $x$ axis respectively:

double center[3]={x0,y0,z0};
B->Rotate(phi,theta,psi,center);

Note that the translation/rotation of a Shape translates/rotates also the inside shapes (see next section).

**To be noticed:** Angles are in radians.

## 4.3 The “put in” operator >>

It is possible to put a `Shape_ptr` A inside another B one, via the "put in"<fnref n="1" /> operator `A>>B`; this operator works differently depending on A:

* if A is a normal `Shape_ptr`: a `A>>B` modifies both A and B; A becomes $A \cap B$ and B becomes $B \cap !A$ (see figure 4.1)

* if A is a `Shape_ptr` with a universe number (e.g. after a call to `Shape_ptr::SetUniverse()` for a lattice): `A>>B` does not modify neither A nor B; a MCNP Fill card is just added to B in order to fill B with universe of A.

Example:

Shape_ptr S(new Sphere(R));
Shape_ptr B(new Brick(a/2,a/2,a/2));
B>>S;

<footnote n="1">¹In C++ this operator is known as the bitwise right shift operator</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 4.1: Operator >> (put in). (a) A and B before the action of the operator ; (b) A and B after the action of the operator.</figcaption></figure>

In this example, the *Brick B* is put inside S.

Note that now S and B are new *Shape_ptr*: S=S & !B (i.e. the sphere without the brick) and B=S & B (i.e, the intersection of the brick and the sphere) as already mentioned above (see fig. 4.1).

**Note:**

* The brick B is an *inside shape* of S: this means that if one moves S, one moves also B.

* It is possible to link inclusion:

A>>B>>C;

means that A is put in B and the result is also put in C. A is an inside shape of B and C; B is an inside shape of C.

* if after the previous example, one makes

C>>D;

Inside shapes of C are cleared and D has A, B and C as Inside Shapes.

## 4.4 Clone of a Shape

It may be sometimes very useful to clone a *Shape*, i.e., to create a new *Shape* with the same properties as the original one. Here is an example:

Shape_ptr B1(new Tube(1,0.5));
Shape_ptr T1(new Tube(1,0.4));
Shape_ptr T2(new Tube(1,0.3));
T2>>T1>>B1;
Shape_ptr B2=B1->Clone();
B2->Translate(1,0,0);

B2 has exactly the same aspect as B1 (with 2 tubes inside) but it is translated by 1m in the *x* direction, whereas B1 is centered at the origin.
## 4.5 Boundary conditions

In MCNP, it is possible to use special boundary conditions: mirror, white or periodic boundary conditions. In Serpent, boundary conditions can be applied only on the shapes associated to cells with the “outside” Serpent’s material. The only valid possibilities for Serpent are black (totally absorbing, the default), mirror or periodic.

**To be noticed also, the “outside” cells of Serpent don’t support unions. For any geometry, if the outside cell (Serpent namespace only) is not a “simple” shape (simple in the Serpent’s sense, it includes Spheres, unrotated Tubes, Bricks or Hexagons), then the user define outside shape (the one which is in the 0 importance cell) is put into a Sphere.**

In Serpent, boundary conditions only apply to Bricks (square, cuboid, ...) or Hexagons (hexxc, hexxprism, ...). Mirror and Black boundary conditions have been implemented in MURE for unrotated Bricks or Hexagons in Serpent. In any case, the use of boundary conditions in MURE could be used but have not been fully tested...so it is advised to use them carefully. In all cases, the exterior of the Shape_ptr with boundary conditions must have a zero importance. You can apply only one type of boundary conditions to a Shape_ptr. Here is a brief description of these boundary conditions:

* Mirror conditions correspond to a standard reflection on a shape ; the method to define such conditions is

```cpp
Shape_ptr A(...);
A->SetMirrorBoundary();
```

* For some shapes, (Hexagons and Bricks) mix boundary conditions can be used (e.g. x-y surfaces are mirrors whereas top and bottom z-planes are open (black)).

* White conditions correspond to a reflection with a cosine direction distribution on the surface ; the method to define such conditions is

```cpp
Shape_ptr A(...);
A->SetWhiteBoundary();
```

* Periodic condition: when a particle leaves a given plane it re-enters through another one. This method could only be applied to Bricks or Hexagons, with the restriction that the Top and Bottom planes (before any reflections) are either Mirror boundaries, White boundaries or Infinite boundaries. The method to define such conditions is

```cpp
Shape_ptr A(new Brick(1,1,1));
A->SetPeriodicBoundary(true, "mirror");
```
or

```cpp
Shape_ptr A(new Brick(1,1,1));
A->SetPeriodicBoundary(true, "white");
```

## 4.6 Recommendations

* It is generally simpler to define a geometry from the inner part to the outer part : for example,
  - define first shape A, then B and C and D.
- Put shapes inside each other like

A>>B>>D;
C>>D;

- After doing that, you can't move and/or rotate A, B or C. But if you move/rotate D, then it will move/rotate also the inner shapes (A, B and D). If A, B or C have to be rotated or translated to a given position in D, then you must do it BEFORE putting these shapes inside D.

* Avoid complex shapes: it is more efficient (in term of CPU time at least for MCNP) to divide a complex in several simpler shapes.

## 4.7 Definition of MC cells

A **Cell** is defined by a **Shape_ptr**, and if needed a **Material**.

MC cells are defined via the *Cell* class; one has to give

* the Shape_ptr corresponding to the geometric shape of the Cell,

* the *Material* (see section 5.1) if exists (default=0 for void),

* the importance of particles<fnref n="2" /> (default=1),

For example to construct a cell composed of a full tube of radius *R* and height *H*, made of $B_4C$

Shape_ptr C(new Tube(H/2,R));
Cell* c=new Cell(C,B4C);

where $B_4C$ is a $(Material*)^3$. In this example, one can define the exterior cell as

Shape_ptr Exterior(!C);
Cell* exterior=new Cell(Exterior,0,0);

The first zero means that the material is void and the second one means that neutron importance is set to 0 in the exterior cell ("outside" material for *Serpent*).

**NOTE:** If not only neutrons are transported in MCNP, it is useful to define the desired mode (i.e., NP, NE, NPE for Neutron and Photon, Neutron and Electron and Neutron, Photon and Electron) before any Cell definition; indeed, in this case, the given cell importance applies to all particle types. For example:

gMURE->SetModeNP();
Cell *c=new Cell(C,B4C);
Cell *exterior=new Cell(Exterior,0,0);

will set the transport mode to neutrons and photons and the cell *c* will have an *IMP:N=1* and *IMP:P=1* whereas the cell *exterior* will have an *IMP:N=0* and *IMP:P=0*.

<footnote n="2">²The importance given in the *Cell* constructor will be the same for all transported particles (define by one of the *MURE::SetMode* methods); if user want to give specific importances for each particle, one needs to call for ALL transported particles the *Cell::AddParticle("particle_Name", importance_for _that_part)*.</footnote>

<footnote n="3">³i.e. a pointer on a *Material*, see Appendix A</footnote>
### 4.7.1 Cell and clone shapes

Suppose we have defined 2 tubes B1 and B2 as follows:

```
Shape_ptr B1(new Tube(1,0.5));
Shape_ptr T1(new Tube(1,0.4));
Shape_ptr T2(new Tube(1,0.3));
T2>>T1>>B1;
Shape_ptr B2=B1->Clone();
B2->Translate(1,0,0);
```

the cell definition for B1 is not a problem:

```
Cell *b1=new Cell(B1,Graphite);
Cell *t1=new Cell(T1,Iron);
Cell *t2=new Cell(T2);
```

where Graphite and Iron are two *Material*. Here, the result will be a tube of graphite containing a tube of iron containing a void tube. To define analog cells for B2 we have to do:

```
Cell *b2=new Cell(B2,Graphite);
Cell *tt1=new Cell(B2->GetOriginalInsideShape(1),Iron);
Cell *tt2=new Cell(B2->GetOriginalInsideShape(0));
```

Indeed, the “Inside shapes” of B2 are ordered from the most inner (`B2->GetOriginalInsideShape(0)` clone of T2), to the most outer (`B2->GetOriginalInsideShape(1)` clone of T1).

## 4.8 Lattice

Lattices are used in MC to fill cells of repeated structure. In *MCNP*, there are 2 types of lattice, hexahedra (type=1) or hexagonal (type=2). In *Serpent*, circular lattice are also defined but they are yet not implemented in *MURE*.

In this section, four examples of lattices are presented (from the most simple to the most complex case). Of course the lattice type does not change the declaration (except that the lattice generator is a *Brick* for hexahedra lattice whereas it is an *Hexagon* for the hexagonal one). Lattice philosophy is somehow different from *MCNP* and *Serpent*; the *MURE* implementation is closer to the *MCNP*'s one.

The general philosophy for lattice declaration is the following:

1. In the Shape section

(a) define the *Shape_ptr C* that will be filled by the lattice (like core),

(b) define the *Shape_ptr G* that is used as lattice generator (a **Brick** or an **Hexagon**) and give to that generator a universe number (via **Shape::SetUniverse**()),

(c) put G in C with $G>>C$

(d) define all *Shape_ptr B*<sub>*i*</sub> that shall be used in the lattice and assign them a universe number (via **Shape::SetUniverse**()) or use **PinCell** (see (a) of Cell Section)

2. In the Cell section

(a) define **PinCell** $P_i$ (and/or use (d) of previous part)
<figure><img src="image_10.png" /><figcaption>Figure 4.2: A simple *implicit* lattice.</figcaption></figure>

(b) build the **Cell** for *Shape_ptr C*

(c) build the **LatticeCell** with the generator *Shape_ptr G*

(d) [optional] indicate to **LatticeCell** the extension of the lattice (via *LatticeCell::SetLatticeRange()* or *LatticeCell::SetNumberOfLatticeElement()*)

(e) fill the lattice with all the desired *Shape_ptr Bᵢ* and/or the **PinCell** *P*ᵢ (via *LatticeCell::FillLattice()*)

### 4.8.1 An implicit lattice example (SimpleLattice.cxx & SimpleLattice_serpent.cxx)

This example illustrates *implicit lattices*, that is to say a lattice filled by a unique universe. A simple tube (a Vessel) is filled by an hexagonal lattice. Each hexagon of the lattice is made of water and filled with a small cylinder of UOx fuel (a pin) inside an Iron cladding. For this example, 2 different implementation are proposed to illustrate the only use of shape (like in the 1.d of the beginning of this section) or the use of **PinCell** (like in 2.a).

As it can be seen, the version using **PinCell** is more compact and more easy to understand (Blue part disappear and magenta part is replaced by green one ; the lattice is thus filled by the **PinCell**). These two implementations give the same result which is shown in Figure 4.2. Note that in the first implementation, outside of the cladding (*OutClad Shape*) has the same universe than the Clad Shape because it is defined as the complementary Shape of Clad (the ! operator). Note also that *OutClad* and *Exterior shapes* must be defined before the “put-in” operator (>>>) because this latter modify the definition of Vessel or Clad. Few words for *Serpent* users: the “lat” card of Serpent need an origin, and a pitch. These quantities are deduced from the generator shape ; translating it will thus change “lat” origin. The sizes of the generator shape will determine the pitch (number of values as well as values depending on the dimension (2D or 3D) of the lattice).
<caption>Table 4.1: SimpleLattice example. Implementation using only shapes (left) and using PinCell (right)</caption>

//
// define materials (see § 5.1)
//
Material *H20= new Material();
H20->SetDensity(1.);
H20->AddNucleus(1, 1, 2.);
H20->AddNucleus(8, 16, 1.);

Material *Iron=new Material();
Iron->SetDensity(7.87);
Iron->AddNucleus(26,56,1.);

Material *UOx= new Material();
UOx->SetDensity(10.2);
UOx->AddNucleus(92, 235, 0.1);
UOx->AddNucleus(92, 238, 0.9);
UOx->AddNucleus(8, 16, 2.);
//
// Geometrical Data Size
//
double VesselH=1;
double VesselR=1;

double HexaH=VesselH;
double HexaSide=0.2;

double PinR=0.12;
double CladR=PinR+0.05;
//
//Shapes
//
//the Vessel : a full Tube
Shape_ptr Vessel(new Tube(VesselH/2,VesselR));
//the Exterior
Shape_ptr Exterior(!Vessel);
//the lattice generator
Shape_ptr LatGene(new Hexagon(HexaH/2,HexaSide));
LatGene->SetUniverse();

// put the lattice generator in the vessel
LatGene >> Vessel;

// a small cylinder that will contain UOx
Shape_ptr Pin(new Cylinder(PinR));
Pin->SetUniverse();

// a bigger one that define fuel clad
Shape_ptr Clad(new Cylinder(CladR));
Clad->SetUniverse(Pin->GetUniverse());

//Outside of the Clad
Shape_ptr OutClad=!Clad;
Pin>>Clad;
//
//Cells
//
Cell *pin=new Cell(Pin,UOx);
Cell *clad=new Cell(Clad,Iron);
Cell *surrounding_clad=new Cell(OutClad,H20);

Cell *exterior=new Cell(Exterior,0,0);
exterior->SetComment("Exterior of the cylinder");
//vessel
Cell *vessel=new Cell(Vessel);
vessel->SetComment("The vessel");
//the lattice
//---------------- Line A -----------------
LatticeCell *Pavage=new LatticeCell(LatGene);
//fill the lattice the Pin shape
Pavage->FillLattice(Pin);
//---------------- Line B -----------------

//
// define materials (see § 5.1)
//
Material *H20= new Material();
H20->SetDensity(1.);
H20->AddNucleus(1, 1, 2.);
H20->AddNucleus(8, 16, 1.);

Material *Iron=new Material();
Iron->SetDensity(7.87);
Iron->AddNucleus(26,56,1.);

Material *UOx= new Material();
UOx->SetDensity(10.2);
UOx->AddNucleus(92, 235, 0.1);
UOx->AddNucleus(92, 238, 0.9);
UOx->AddNucleus(8, 16, 2.);
//
// Geometrical Data Size
//
double VesselH=1;
double VesselR=1;

double HexaH=VesselH;
double HexaSide=0.2;

double PinR=0.12;
double CladR=PinR+0.05;
//
//Shapes
//
//the Vessel : a full Tube
Shape_ptr Vessel(new Tube(VesselH/2,VesselR));
//the Exterior
Shape_ptr Exterior(!Vessel);
//the lattice generator
Shape_ptr LatGene(new Hexagon(HexaH/2,HexaSide));
LatGene->SetUniverse();

// put the lattice generator in the vessel
LatGene >> Vessel;

// Cells
//
PinCell* pinUox=new PinCell;
pinUox->AddLayer(Uox,PinR);
pinUox->AddLayer(Iron,CladR);
pinUox->SetSurroundingMaterial(H20);

Cell *exterior=new Cell(Exterior,0,0);
exterior->SetComment("Exterior of the cylinder");
//vessel
Cell *vessel=new Cell(Vessel);
vessel->SetComment("The vessel");

//the lattice
//---------------- Line A -----------------
LatticeCell *Pavage=new LatticeCell(LatGene);
//fill the lattice with the pincell
Pavage->FillLattice(pinUox);
//---------------- Line B -----------------
### 4.8.2 An explicit lattice example with different zones (SimpleLattice2.cxx & Simple-Lattice2_serpent.cxx)

Explicit lattices are filled with different universes, allowing to put different structures at some given position. To illustrate this, let us modify the previous example slightly: we want to fill the Vessel with 2 types of hexagons: the first ones are full water hexagons (the reflector) in the outer part of the Vessel and the second ones are the previous hexagons with the small pin rod inside. For shortness, we chose to fill the Vessel with second hexagon type if it is inside a "virtual" cylinder of radius 90cm. Both implementations using Shape only or PinCell has been proposed in Tab. 3.3. Here, a focus is made on the PinCell implementation. The principal changes take place between the lines named "Line A" and "Line B". We have to replace this part by:

```cpp
//---------------- Line A ------------------
PinCell* whole=new PinCell;
whole->SetSurroundingMaterial(H2O);
LatticeCell *Pavage=new LatticeCell(LatticeGenerator);
int range=int(VesselR/(1.5 * HexaSide))+1;
Pavage->SetLatticeRange(-range,range,-range,range); //define an explicit lattice
for(int i=-range; i<=range; i++)
    for(int j=-range; j<=range; j++)
    {
        int pos[3]=[i,j,0]; //the lattice index
        double xt = HexaSide*sqrt(3.) * (i + j*0.5 ) ;
        double yt = 1.5 * HexaSide * j;
        double X[2]=[xt,yt]; //the center of each hexagons
        if (IsHexagonInTube(X,LatticeGenerator,0.9)) //true if the hexa center at X is
            Pavage->FillLattice.pinUox(pos); //inside a tube of radius 90cm
        else
            Pavage->FillLattice(whole pos); //inside a tube of radius 90cm
    }
//---------------- Line B ------------------
```

* A new **PinCell** (*whole*) is defined; because just a surrounding material is defined, this **PinCell** represents the whole space (it is a "pseudo-PinCell").

* Then the lattice range is defined as $(2range+1) \times (2range+1)$ matrix (in *x* and *y* direction and infinite in *z*). The "range" has been defined to cover at least the entire vessel.

* Then, each position of hexagons in the lattice is computed; if a hexagon is inside the virtual cylinder of radius 0.9m, a hexagon with the pin rod is used (**PinCell** *pinUox* in the *Cell::FillLattice* method) at the given position *pos* else the "whole" **PinCell** is used (and as always with lattices, it is automatically cut by the lattice meshes).

The result is shown in Figure 4.3.

Note: if one prefers use only shapes (and not **PinCell**), one can use the following trick (**but it is not implemented for Serpent!!**)

* In the *Shape* part, define a *Shape_ptr* corresponding to the *Whole* space (intersection of nothing<fnref n="4" />):

<footnote n="4">⁴Because we are in a lattice, this whole space is cut by the lattice generator. It is also possible to give a bigger volume to obtain the same kind of result; however, MCNP has to check all surfaces of this volume for a neutron in the cell so it will result in a longer MCNP run...</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 4.3: A lattice with different zones.</figcaption></figure>

```cpp
Shape_ptr Whole(new Node(-1)); // -1=intersection, here, of nothing
Whole->SetUniverse();
```

* and in the Cell part, add

```cpp
Cell *whole=new Cell(Whole,H2O);
```

### 4.8.3 A lattice with more than one simple shape (Stadium.cxx & Stadium_serpent.cxx)

In this example, we consider a hexahedra lattice. Each brick of the lattice will have parts of “stadiums” on each (x,y) sides (see Figure 4.4, this is the base geometry of the MSRE reactor).

```cpp
//
//materials
//
Material *Salt=new Material(); // this is the fuel LiF-ThF4/UF4
Salt->SetDensity(3.1);
Salt->AddNucleus(3,7,0.7);
Salt->AddNucleus(9,19,1.9);
Salt->AddNucleus(90,232,0.15);
Salt->AddNucleus(92,233,0.15);
Material *Graphite=new Material();
Graphite->SetDensity(1.86);
Graphite->AddNucleus(6,0,1);
//
// Data
//
double VesselH=0.8; // the core H and R
double VesselR=0.4;
double SquareH=VesselH; // stadium height
double SquareS=0.2; // stadium straight line length
```
double BendR=0.02; //stadium bend curve radius

double BendH=VesselH;

double BrickW=0.04;

double BrickL=0.08;

//
//Shapes
//
//the Vessel
Shape_ptr Vessel(new Tube(VesselH/2,VesselR));
//the Exterior
Shape_ptr Exterior(!Vessel);
//the lattice generator
Shape_ptr LatticeGenerator(new Brick(SquareS/2,SquareS/2,SquareH/2));
LatticeGenerator->SetUniverse();
//put the lattice generator inside the core
LatticeGenerator>>Vessel;
//Definition of the stadiums (horizontal & vertical)
Shape_ptr Circle0(new Tube(BendH/2,BendR));
Circle0->SetUniverse();
//clone circles for horizontal and vertical stadiums
Shape_ptr Circle1=Circle0->Clone();
Shape_ptr Circle2=Circle0->Clone();
Shape_ptr Circle3=Circle0->Clone();
//translate circle
Circle0->Translate(-0.1,+0.04,0);
Circle1->Translate(-0.1,-0.04,0);
Circle2->Translate(-0.04,0.1,0);
Circle3->Translate(+0.04,0.1,0);
//2 bricks for horizontal and vertical stadiums
Shape_ptr Brick0(new Brick(BrickW/2,BrickL/2,SquareH/2));
Brick0->SetUniverse(Circle0->GetUniverse()); //belong to the same universe than Circle0
Brick0->Translate(-0.1,0,0);
Shape_ptr Brick1(new Brick(BrickL/2,BrickW/2,SquareH/2));
Brick1->SetUniverse(Circle0->GetUniverse()); //belong to the same universe than Circle0
Brick1->Translate(0,0.1,0);
//2 vertical stadiums
Shape_ptr Stade0=Circle0 | Brick0 | Circle1; //union of circles and brick
Shape_ptr Stade2=Stade0->Clone();
Stade2->Translate(0.2,0,0);
//2 horizontal stadiums
Shape_ptr Stade1=Circle2 | Brick1 | Circle3;
Shape_ptr Stade3=Stade1->Clone();
Stade3->Translate(0,-0.2,0);
//interior of all stadiums
Shape_ptr AllStades=Stade0 | Stade1 | Stade2 | Stade3; //to define all stadium exterior

//
// Cells
//
Cell *exterior=new Cell(Exterior,0,0);
exterior->SetComment("The exterior of the reactor");
//vessel
Cell *vessel=new Cell(Vessel);
<figure><img src="image_12.png" /><figcaption>Figure 4.4: A lattice with more than one shape.</figcaption></figure>

```cpp
vessel->SetComment("The reactor vessel");

//Lattice
LatticeCell *Pavage=new LatticeCell(LatticeGenerator);
Pavage->FillLattice(AllStades);
Cell *stade0=new Cell(Stade0,Salt);
Cell *stade1=new Cell(Stade1,Salt);
Cell *stade2=new Cell(Stade2,Salt);
Cell *stade3=new Cell(Stade3,Salt);
Cell *exterior_stade=new Cell(!AllStades,Graphite);
```

* In the Shape part:

- A *Tube* (Circle0) is created at the origin; then it is cloned in 3 other *Tubes* (universe number of Circle0 is given to its clones). Then these *Tubes* are translated to the right place.

- 2 *Bricks* are created (with the same universe as Circle0) and translated to the left and the top of the square. All elements of all stadiums have, here, the same universe.

- Then 2 stadiums are constructed as the union of 2 tubes and a brick. The right and bottom stadiums are cloned from the 2 previous ones, and a union of all stadiums is constructed.

* In the Cell part:

- The lattice cell is defined, and filled by the universe of *AllStades*.

- then, stadiums cells and exterior of these stadium are built.

### 4.8.4 Lattice of a Lattice (LatticeOfLattice.cxx & LatticeOfLattice_serpent.cxx)

The purpose of this example is to illustrate a hexagonal Lattice of big structural hexagons, themselves filled with a hexagonal lattice of fuel rods (like in a SFR core)
//
// Material
//
Material *Iron=new Material();
  Iron->SetDensity(7.87); //density in g/cm3
  Iron->AddNucleus(26,56,1.);

  Material *Na=new Material();
  Na->SetTemperature(600);
  Na->SetDensity(0.97);
  Na->AddNucleus(11,23,1);
  Material *Na_cold=Na->Clone(300); //just to see in Serpent all hexagons
  Material *Fuel=new Material();
  Fuel->SetTemperature(900);
  Fuel->SetDensity(19.);
  Fuel->AddNucleus(92,238,0.8);
  Fuel->AddNucleus(94,239,0.2);
//
// data size
//
double VesselH=1.1; //the core size
double VesselR=1.1;
double FuelHexaSide=0.05; // (small) fuel hexagons side
double PinRadius=0.02; //fuel pin radius
double StructHexaSide=8./sqrt(3.)*FuelHexaSide; // (big) structure hexa
double StructThick=0.01; //the thickness of the structure (inner
// part of this structure is filled with
// FuelHexa.
//
// Shapes
//
//the Core
Shape_ptr Vessel(new Tube(VesselH,VesselR));
//the Exterior
Shape_ptr Exterior(!Vessel);
//the outside structure hexagon
Shape_ptr StructHexa(new Hexagon(VesselH,StructHexaSide+StructThick));
StructHexa->SetUniverse();
//the inside structure hexagon
Shape_ptr InnerStructHexa(new Hexagon(VesselH,StructHexaSide));
InnerStructHexa->SetUniverse();
// put the structure hexagons in the vessel
StructHexa>>Vessel;
//
//a "fuel hexagon"
//
Shape_ptr FuelHexa(new Hexagon(VesselH,FuelHexaSide));
FuelHexa->Rotate(-Pi/2);
FuelHexa->SetUniverse();
//put the fuel hexagons in the structure (inner part)
FuelHexa>>InnerStructHexa;
//
// Cells
//
PinCell *Reflector=new PinCell(); //Reflector=beetwen Core and big
Reflector->SetSurroundingMaterial(Na_cold); //hexagons. This is "cold Na" mat
PinCell *FuelHexaBorder=new PinCell(); //Fill the inner part of big
FuelHexaBorder->SetSurroundingMaterial(Na_cold); //hexagons with cold Na"
PinCell *Carandache=new PinCell(); //This is a fuel pin with
Carandache->AddLayer(Fuel,PinRadius); //a cylinder of fuel
Carandache->SetSurroundingMaterial(Na); //surrounded by Na
Cell *exterior=new Cell(Exterior,0,0); //out side of the core
exterior->SetComment("The exterior of the reactor");
Cell *vessel=new Cell(Vessel); //the core
vessel->SetComment("The reactor vessel");
Cell *innerstructhexa=new Cell(InnerStructHexa);
innerstructhexa->SetComment("The Inside struct Hexa");
Cell *exterior_innerstructhexa=new Cell(!InnerStructHexa,Iron); //space beetwen InnerStructHexa
exterior_innerstructhexa->SetComment("The Big Hexa claddings");//and StructHexa fill with Iron
LatticeCell *structhexa=new LatticeCell(StructHexa);
structhexa->SetComment("The Ouside struct of Big Hexa (lattice)");
int range1=int(VesselR/(1.5 * StructHexaSide)+1);
structhexa->SetLatticeRange(-range1,range1,-range1,range1);
double StructHexaWidth=sqrt(3.)*StructHexaSide;
for(int i=-range1; i<=range1; i++)
for(int j=-range1; j<=range1; j++)
{
    int pos[3]=[i,j,0];
    double xt = StructHexaWidth * (i + j*0.5 );
    double yt = 1.5 * StructHexaSide * j;
    double X[2]=[xt,yt];
    if (IsHexagonInTube(X,StructHexa,VesselR))
        structhexa->FillLattice(InnerStructHexa pos);
    else
        structhexa->FillLattice(Reflector pos);
}
LatticeCell *fuel_hexa=new LatticeCell(FuelHexa);
fuel_hexa->SetComment("The fuel hexagonal lattice");
int range=int(StructHexaSide/(1.5 * FuelHexaSide)+1);
fuel_hexa->SetLatticeRange(-range,range,-range,range);
double FuelHexaWidth=sqrt(3.)*FuelHexaSide;
for(int i=-range; i<=range; i++)
for(int j=-range; j<=range; j++)
{
    int pos[3]=[i,j,0];
    double xt = 1.5 * FuelHexaSide * i;
    double yt = FuelHexaWidth * ( j + i*0.5 );
    double X[2]=[xt,yt];
    if (IsHexagonInHexagon(X,FuelHexa,InnerStructHexa))
        fuel_hexa->FillLattice(Carandache pos);
    else
        fuel_hexa->FillLattice(FuelHexaBorder pos);
}
<figure><img src="image_13.png" /><figcaption>Figure 4.5: Lattice of a Lattice.</figcaption></figure>
Chapter 5

Materials, Sources, Tallies, ...

**5.1 Definition of material**

A **Material** is a set of nuclei with their proportions (default units is %mol), it has a density and a temperature.

For example to define $B_4C$ of density 2.52 g/cm$^3$ at 296 K:

Material *B4C=new Material();
B4C->SetDensity(2.52); //default units are g/cm3
B4C->AddNucleus(5,10,0.199); //Boron is composed 19.9% of $^{10}_{5}B$
B4C->AddNucleus(5,11,0.801); // and 80.1% of $^{11}_{5}B$
B4C->AddNucleus(6,0,0.25); // This is natural C (Z=6, A=0) ; for 1 mol of B one has 0.25 mol of C

Just before writing Material in MC file, proportions are normalized to 1. The *Material::AddNucleus()* method has the general form

Material::AddNucleus(Z,A,I,prop);

or

Material::AddNucleus(Z,A,prop);

where Z is an integer corresponding to the proton number of the nucleus, A, an integer for the number of nucleons, I, an integer for the isomeric state and prop, a double for the proportion of the nucleus in the *Material*. For **natural isotopes A=0**; for **ground state I=0** (this is the default value when not specified). See § 5.1.3 for units.

* A *Material* can contain one nucleus which is a “moderator”: this means that the special treatment $S(\alpha, \beta)$ is applied for low energy in *MC* run for this nucleus. In order to say that a nucleus is a moderator, one has to add a string for the last parameter of *Material::AddNucleus* (see 5.1.4 for a full description) ; defining thermal treatment for hydrogen in light water leads to:

H2O->AddNucleus(1,1,2.,kpMOL,"H2O"); //here the unit of proportion must be given explicitly

H2O->AddNucleus(8,16,1.);

**5.1.1 Clone**

A *Material* can be cloned via the *Material::Clone(double Temperature)* method. There are 2 types of clones depending on:
* whether the material does not evolve and the Temperature argument of `Material::Clone()` is negative (the default) (see example of `B4C_Light` in § 5.1.7): the clone is only used because it makes possible to define a material with the same number in `MCNP` but also present at different densities in different cells. A call to `Material::SetTemperature()` will have no effect on such a clone!

* or whether the material will evolve and/or the Temperature argument is positive or null (see example of `Fuel2` and `B4C_Hot` in § 5.1.7): the clone is a "true" clone, i.e., all compositions and proportions are duplicated (2 different material numbers in `MCNP`). Of course the density, temperature,... of the clone could be changed after cloning, as well as new nuclei could be added. This kind of clone is the correct way to handle 2 evolving cells starting with the same composition but evolving independently. Note that by default the density of the clone is the one of the original material.

### 5.1.2 Mix (Mix.cxx)

Two Materials of different densities can be mixed together to create a third material. This is useful in many situations, e.g., the case of reactor fuels with a given enrichment of fissile isotopes. The principle is the following:

M12=M1->Mix(M2,part,kpMOL);

M12 is a new material and is created from the molar part mixture of materials M1 and M2 (for example Uranium and Plutonium Oxides). M1 and M2 can not be used in the problem, because the `Material::Mix` method forces them to be fictitious materials by using the `Material::ForbidPrint()` method: **THIS IS DONE IN THE Material::Mix() METHOD.** The density of the new material M12, is calculated as

$$ \rho_{12} = \frac{1}{\omega_1/\rho_1 + \omega_2/\rho_2} $$

where $\omega_i$ and $\rho_i$ are the mass fraction and density of material i. The temperature of the new material is the same as that of material M1.

### 5.1.3 Units

#### Density

Default units (k=constant, d=density) for density is `kdGCM3` (g/cm³); but, to avoid confusion, you must set the density unit in `Material::GetDensity()` (either `kdGCM3` (g/cm³), `kdATCM3` (at/cm³) or `kdATBCM` (at/barn.cm)). User may give the density in `atom/barn.cm` by using `Material::SetDensity(density, kdATBCM)`.

**NOTE:** If the density is not given (i.e. set to the default value 0), it could be automatically found if the proportions are given in `atoms/barn.cm` or in `atoms/cm³`.

#### Proportions

Default units for proportions are in molar percentage. To give proportion in other units, use

`Material::AddNucleus(Z,A,I,prop,UNITS)`

where `UNITS` stands for (k=constant, p=proportion): `kpMOL` (`"%mol"` (default)) or `kpMASS` (`"%mass"`) or `kpATBCM` (`"at/barn.cm"`) or `kpATCM3` (`"at/cm3"`).
### 5.1.4 Material extension for MC

In *MCNP* or *Serpent*, material codes have extensions (like ".60c"). These extensions are very important because they distinguish isotopes at different temperatures, bases, versions, ... In practice, if no extension is specified in *MCNP* file, the first isotope encountered in the *xsdir* is chosen regardless the temperature imposed in the cell where the material is (but this is not the case in *Serpent*). In *MURE* user has different ways to specify an extension.

* This can be done manually for each nucleus of the material composition:

```
mat->GetNucleus(Z,A,I)->SetXSExtension(".49c")
```

* This can be done manually for all nuclei of the material composition:

```
mat->SetDefaultXSExtension(".49c")
```

then each nucleus of the Material *mat* will have this ":.49c" extension.

* If the *MCNP* code for a particular Nucleus (Z,A,I) of a given Material *mat* is known, it can be specified by using:

```
mat->GetNucleus(Z,A,I)->SetXSExtension(".49c")
```

* This can be done automatically (this is the recommended way): the specific file, the so-called *BaseSummary.dat* file, which contains *xsdir/xsdata* information written in a different way for search efficiency, is used. The extension is chosen according to the *Material* temperature. In fact the closest temperature to the one of the *Material* is chosen with the help of classes *BasePriority* and *TemperatureMap*. (see below)

**NOTE:**

* If evolution is needed, the automatic procedure is certainly required... but I am not sure (PTO).

* If one specifies the extension value (either via `Nucleus::SetXSExtension()` or `Material::SetDefaultXSExtension()` or `Nucleus::SetModeratorName()`) then you must provide in the *xsdir* the corresponding data files.

#### S(α, β) Treatment and MT card

When the special treatment *S*(α, β) is desired for low energy, *MCNP MT* card arguments or *Serpent moder* and *therm* cards have to be specified. Consider the water example:

```
H2O->AddNucleus(1,1,2.,kpMOL,"H2O");
```

```
H2O->AddNucleus(8,16,1.);
```

* this can be done manually: the string "H2O" must not be empty; then

```
H2O->GetNucleus(1,1)->SetModeratorName("lwtr.01t");
```

* this can be done automatically: the *BaseSummary.dat* file is read to find the code of *S*(α, β) treatment file; this time, because *S*(α, β) treatment is material dependent, the string has to refer to a given "Category" to specify which kind of effects are taken into account; in the above example, hydrogen has different *S*(α, β) treatment for light water, polyethylene, benzene, ... The complete list of available Categories is given in **Nucleus** in the *Nucleus::SetModeratorCategory()*:.
<table><tr><td>H2O</td><td>for H in light water</td></tr><tr><td>H/Zr</td><td>for H in ZrH</td></tr><tr><td>poly</td><td>for H in polyethylene</td></tr><tr><td>D2O</td><td>for D in heavy water</td></tr><tr><td>BeO</td><td>for Be in Be oxide</td></tr><tr><td>Be</td><td>for Be in Be metal</td></tr><tr><td>Gr</td><td>for graphite</td></tr><tr><td>Zr/H</td><td>for Zr in ZrH</td></tr></table>

The moderator category can also be any real MT file name such as “lwtr” or “u/o2”, ...

### 5.1.5 Pseudo Materials

PseudoMaterials are particularly useful for thermal-hydraulics problems. Normally, the material temperature dictates which extension (e.g. “.60c”) for cross-sections will be used from the evaluated data libraries. The default behaviour in MURE is that the nearest available temperature is used, the maximum tolerable difference between this temperature and the desired material temperature being defined in the TemperatureMap class. However, there are occasions where it is useful to perform a stochastic interpolation at the desired temperature by adding two identical nuclei at different temperatures in the material in varying proportions. These types of materials are called PseudoMaterials, because, in principle, they behave as a material at a precise temperature, T, even though the cross-sections in the evaluated data libraries at this precise temperature do not exist. To declare a Material (Mat) to be a PseudoMaterial requires the following simple line:

Mat->SetPseudoMaterial();

From this point onward, MURE will automatically create two identical nuclei with different extensions/temperatures and in the correct proportions for, every nucleus which is added to this material. No other commands are necessary. However, in the case of an evolving material (a reactor fuel rod for example), this can cause the addition of many, many extra nuclei to the material, thus slowing down MCNP calculations by perhaps as much as 20%. In the case that the user only wants to declare certain nuclei to be PseudoNuclei (for example, the most important ones, such as the fissile and fertile nuclei) the following kind of lines can be added after the creation of a Material:

UO2Fuel->AddNucleus(92,235,0.04);

UO2Fuel->AddNucleus(92,238,0.96);

UO2Fuel->AddNucleus(8,16,2.0);

UO2Fuel->GetNucleus(92,235)->SetPseudoNucleus();

UO2Fuel->GetNucleus(92,238)->SetPseudoNucleus();

Now, only the selected nuclei will be treated as pseudo nuclei and automatically inserted twice into the MCNP file at the two nearest data base temperatures in the correct proportions. The proportions, $\omega_1$ and $\omega_2$, for the interpolation are calculated in the following way:

$$ \omega_2 = \frac{\sqrt{T} - \sqrt{T_2}}{\sqrt{T_2} - \sqrt{T_1}} $$

$$ \omega_1 = 1 - \omega_2 $$
where $T_1$ and $T_2$ are the nearest two temperatures in the evaluated data base files straddling the desired temperature $T$.

### 5.1.6 MC material Printing

This section concerns only evolving Materials. When evolution is required, at the beginning, many nuclei are present in the material but in a very small proportion. Each time a neutron enters in a cell, MC code tries to find the nucleus on which the potential reaction could happen. Thus, the more complex is a material, the longer is the MC run. The idea is to suppress only for the MC print, the nuclei which do not contribute "significantly" to neutron transport to save CPU time. The criterion is based on the evaluation of total cross-section $\sigma_T$ and nucleus proportion $p$. If a nucleus of a material has a $p \times \sigma_T < \varepsilon$, it is not printed. The value of $\varepsilon$ is by default set to $0.1\% \times 1\,\text{mb} = 10^{-6}$. This can be changed via `MURE::SetMCNPNucleusThreshold()`. To be noticed that perturbative tallies of all evolving nuclei are written (in order to evaluate the cross-section...). The effect of this procedure is thus to accelerate MC runs at the beginning of the evolution where time steps between MC are generally smaller, and add more and more nuclei as far as the evolution is going on (with larger time steps between MC runs).

### 5.1.7 Examples

Please take a look at `Material_test.cxx` and `Material_test2.cxx` and the corresponding `Serpent` version for the second one. The use of `Nucleus::SetXSExtension` or `Material::SetDefaultXSExtension` supposed that the user provide a correct `xsdir` file with all necessary information.

### 5.1.8 Automatic extension finding

Here is a short description of the automatic nucleus extension finding.

#### Making the `BaseSummary.dat` file

In order to work, a specific file, `BaseSummary.dat`, must be present in the data directory. The aim of this file is to facilitate the search for a given Z,A,I and related `xsdir/xsdata` information (see also § 1.4.2 and 1.4.3).

#### Extracting the MC nucleus code from the `BaseSummary.dat` file

Suppose the desired ZAI exists in the `BaseSummary.dat` file. Then the choice of the desired nucleus is based on 2 criteria: the first one is related to the preference in a nuclear base (ENDB, JEFF, JENDL) and the version of this base, and the second criterion is related to the nucleus temperature. A score is calculated for each criterion, and the best score is chosen. Two classes are in charge of these criteria: `BasePriority` and `TemperatureMap`.

* `BasePriority` class: a priority base is organized in 3 levels (base names) and maximum 3 sub-levels (base versions). Level 1 has a greater priority than level 2 and so on. By default, the 3 levels are: (1) ENDFB (USA), (2) JEFF (European) and (3) JENDL (Japan). Only 1 sub-level is defined for each base (6.1 for ENDFB, 3.0 for JEFF and 3.3 for JENDL). If a user defines for each level the same base (e.g., (1) ENDFB, (2) ENDB, (3) ENDB), then this imposes the base choice (i.e., no other choice is possible).

For example

```
BasePriority *BP=new BasePriority();
BP->SetBasePriority(0,"JEFF"); //1st prefered base
```
BP->SetVersionPriority(0,"3.1",0); //1st prefered version for 1st prefered base
BP->SetVersionPriority(1,"2.2",0); //2nd prefered version for 1st prefered base
BP->SetBasePriority(1,"ENDFB"); //2nd prefered base
BP->SetBasePriority(2,"JENDL"); //2nd prefered base
gMURE->SetBasePriority (BP);

will choose nuclei first in JEFF 3.1, the if not found, in JEFF 2.2, then in JEFF (what ever is the version), then
in ENDFB (default first preferred version VI.8) and then in JENDL (default first preferred version 3.3).

If one want to impose nuclei in ENDFB VII.0, then user should do

BasePriority *BP=new BasePriority();
for(int i=0; i<3; i++)
{
  BP->SetBasePriority(i,"ENDFB");
  for(int j=0; j<3; j++)
    BP->SetVersionPriority(j,"7.0",i);
}
gMURE->SetBasePriority (BP);

In that case, only nuclei given in the wanted nuclear base are used.

*   **TemperatureMap class:** this class defines a temperature "binning" $T[i]$ and a precision $\Delta T$; an input temperature $T_{in}$ can be searched in this binning with respect to the precision. For example, if $T[0] = 300$ K, $T[1] = 500$ K and $\Delta T = 50$ K, then all $T_{in}$ in $[250 \text{ K}, 350 \text{ K}]$ will have a score of $\frac{|T_{in} - T[0]|}{\Delta T}$, all $T_{in}$ in $[350 \text{ K}, 450 \text{ K}]$ will have a 0 score. The precision $\Delta T$ can be set via `gMURE->GetTemperatureMap()->SetDeltaT(TemperatureMap)`. A zero score in temperature stops the MURE execution: you may either increase the $\Delta T$ precision or build the nucleus data file (ACE format) at the desired temperature.

### 5.1.9 Automatic XSDIR construction

Using the extension finding, MURE can build the `XSDIR/XSDATA` file that is used by `MC` code to get the cross-sections used in the problem. To use this option, you just have to set :

gMURE->SetAutoXSDIR();

When getting the code for `MC` file, `MURE` will then save the best line it finds in `BaseSummary.dat`. Therefore, right before building the `MC` file (the `XSDIR` is built at the first step of an evolution), `MURE` will build the `XSDIR/XSDATA` file by using all these best lines it saved.

**IMPORTANT**: If an existing `XSDIR/XSDATA` is already present in the `MC` run directory, then `MURE` will **NOT** build a new one ; the existing one will be used by `MC` code. Thus either the `xsdir/xsdata` is correct, either you have to remove it.

## 5.2 Particle Source

Because the difference between particle source definition in `MCNP` and `Serpent` is great, `MURE` source definition uses 2 different classes, `MCNPSource` and `SerpentSource` to avoid confusion for the user ; they are sibling of the `MCSource` class.
The *MCNPSource SerpentSource* classes enable the definition of particle sources (up to now only neutron sources may be defined for *Serpent*). Nevertheless, some source definition methods are common to both classes. Two type of source can be defined, source for criticallity calculation (called *kcode mode*) and fix sources for non fissile or subcritical system only. For *Serpent*, source possibilities are much more simpler than for *MCNP*; in *MURE*, the source implementation for *Serpent* is also much less developed. Mainly position and energy of punctual monoenergetic isotrope source can be defined. Criticallity is defined also by the “*SerpentSource::SetKcode*” method. In the following section, the *MCNPSource* is used because of its greater complexity, but looking to *SerpentSource* methods and their analog counterpart in *MCNPSource* is self understandable.

The source may be defined with *MCNP* card or as an external source. In the latter case (*MCNPSource::SetExtern()*), external source will be read in a special file. This is valid in subcritical mode as well as in critical (*KCODE*) mode (not true for *Serpent*). One chooses to run in subcritical mode (default) or in *KCODE* mode via the *MC-NPSource::SetKcode()* method.

### 5.2.1 Setting Source for MURE

Once an *MCSource* has been defined, one has to give this source to *MURE* class:

```
MCNPSource *mysource=new MCNPSource(1000); //1000 source particles
gMURE->SetSource(mysource);
```

### 5.2.2 Criticallity calculation: Kcode mode

To be noticed that, for *Kcode*, the number of source particles defined in the constructor of *MCSource* is in fact the number of neutrons per cycle.

To initiate a *KCODE*, a *SDEF* (or *KSCR*) card is needed (if no external source is given!). By default in *MCNP*, the position of this initial source is (0,0,0). If you want to change the source position, use the *MCNPSource::SetPosition()*. In *Serpent* the default positions are chosen in cells containing fissile materials.

```
MCNPSource *mysource=new MCNPSource(1000); //1000 source particles/cycle
//Kcode with first expected keff=1, 10 inactive cycles and 80 active ones
mysource->SetKcode(80,10,1);
gMURE->SetSource(mysource);
```

If one want to use a source from a previous *KCODE*:

```
gMURE->UsePreviousRunSource("olds");
```

where the source file *olds* is the result of a previous run. In case of evolution, if you do not specify any argument, the source will be calculated as specified by inactive cycles, and then for next *MCNP* runs, the source is taken from the previous cycle.

### 5.2.3 More elaborated sources (MCNP only)

The *MCNPSource* gives more details on *MCNPSource* class possibilities. Two simple examples are presented: First
```
MCNPSource *source=new MCNPSource(1000); //1000 particles per cycle
source->SetKcode(500,20,1); //critical mode (500 active cycles, 20 inactive, supposed keff=1)
source->SetKSRC(); //use fission spectrum for netron energy
source->AddPosition(0,0,0); //add one position for the source at (0,0,0)
source->AddPosition(0.1,0.1,0.1); //add an other one at (10cm,10cm,10cm)
source->AddPosition(0,0,0.2); //add an other one at (0,0,20cm)
source->AddPosition(0,0,-0.2); //add an other one at (0,0,-20cm)
```

and the second

```
MCNPSource *source=new MCNPSource(); //SDEF source with 10000 particles
source->SetEnergie("D1"); //energy is defined by the source probability 1
source->SetDistribution("RAD=D2 EXT=D3"); // radial and axial extension are
                                        // defined by source information 1 and 2
source->AddBias("SP1 -5 2.0"); //p(E)=E*exp(E/2.MeV)
source->AddBias("SI2 1. 2."); // radial extension from r=1cm to r=2cm
source->AddBias("SI3 -10. 10."); // axial extension from z=-10cm to z=10cm
```

### 5.2.3.1 Source define via Spectrum class (MCNP only)

One can define energy distribution source using the `Spectrum` class and the `MCNPSource::UseThisEnergyDistribution()` method. Any call to that method will dump in `MCNPSource_YourEnergyDistribution.dat`, the ASCII spectrum of this source (for user own use) as well as the total number of particles emitted (for normalization purpose).

**To be noticed: this method will change the source particle type according to the Spectrum type (neutrons for NeutronSpectrum, γ for GammaSpectrum, ...).**

This method only change energy distribution of the source (using `MCNP` distribution source cards (D, SI, ...)) numbered by default from 800).

### 5.2.3.2 Tube Source (MCNP only)

It is possible to define volumic cylindrical sources with the help of the `Tube` class. Two methods are available: `MCNPSource::AddTubeSource(Tube *tube, string energy)` and `MCNPSource::AddTubeSource(Tube *tube, Spectrum* spectrum)`. In both cases, the `tube` is used to define the shape of the source. The first one allows to use either mono-energetic sources or using source distributions "D" card. **In the latter case, the user need to provide himself the distribution description via the `MCNPSource::AddBias()` method (see TubeSource.cxx example).** The method using `Spectrum` class argument use `Spectrum` to define both particle type of the source (neutrons for NeutronSpectrum, γ for GammaSpectrum, ...) and energy spectrum (see TubeSource2.cxx example).

These 2 methods use internally source distribution cards (D, SP, SI, ...). The starting number of these distribution is 900. Thus don't use number above 900 to define your own source distribution and biasing.

Of course you can add more than one `Tube` source; but **it is only possible to add more Tube sources with same calling method (either MCNPSource::AddTubeSource(Tube *tube, string energy) or MCNPSource::AddTubeSource(Tube *tube, Spectrum* spectrum). One cannot mix them.**

### 5.2.3.3 Define MCNP Source from the result of an evolution (MCNP only)

If one want to use the energy distribution of a spent fuel which has cooled for a time `t`, use `MureGui` and dump the source in a `MURE` format ("Save Data" button, then select the "MURE INPUT Spectrum" radio widget). The generated file can be inserted in a `MURE` "cxx" code (c.f. `MureGui` section 8.3.6).
## 5.3 Tally class

The *Tally* class enables tallies for *MCNP/Serpent* to be defined.

**WARNING: NEVER USE: "Tally *f=new Tally[N]" BUT USE INSTEAD "vector<Tally*> f(N); for (int i=0; i<N; i++) f[i]=new Tally;". The reason is that when tallies are deleted, one need to use either "delete []" or "delete"; in MURE, the second case is chosen which is incompatible with the former declaration.**

*MCNP* tallies or equivalent *Serpent* detectors are used to store quantities from a run; *Serpent* detectors and *MCNP* tallies have lots of common points but the *MCNP* version has more possibilities (Group tallies, tally normalization, ...). The implementation of tallies in *MURE* is (at least for historical reasons) much more developed for *MCNP*. Nevertheless, **surface tallies (type 1 and 2 in *MCNP*) are not fully implemented (Surface calculation, ...) even if one can use them with special precaution/verification.**

There are 2 main tally types: Surface tallies (type 1 and 2 in *MCNP*, same convention is kept for *Serpent*) and Cell tallies (4,...). At present, only *Cell* tally of type 1 (current through a surface), 2 (flux through a surface), 4 (flux in cell), 6 (energy deposition in a cell) and 7 (fission energy deposition in a cell) are implemented in *MURE*. **Surface tallies have not been tested at all and stochastic surface calculation is not implemented in *MURE*.**

A *Tally* consists of

* a tally type (what is the desired quantity, e.g. flux in a cell),

* a list of *TallyBin*: this may be *SimpleBin* (a cell or surface number), *GroupBin* (a group of cells or surfaces) or a *LatticeBin* which allows to obtain tallies in cells with universes. **The GroupBin does not exist in *Serpent*; the LatticeBin is not implemented in *MURE* for *Serpent*.**

* a Time and/or Energy binning (linear, log or arbitrary binning)

* a list of *TallyMultiplicator* (*Tally* Multiplicator) that allows to obtain cross-sections, .... A *TallyMultiplicator* is composed of a multiplicative constant, a material number and a *Reaction*. The multiplicative constant is not used in *Serpent*.

Of course *Cell* tally bin types could not be mixed with surface tally bin type. Tally bins are added with the *Tally::Add()* method (for cell/surface bins) whereas energy and/or time bins are added with *Tally::AddEnergy()* or *Tally::AddTime()*.

*Tally* multiplicator may be added with *Tally::AddMultiplicator()* to obtain cross sections, reaction rates and so on. Example of Tallies:

```
Tally *t1=new Tally(); //default tally type is 4: flux in cell
t1->Add(cell_1);
t1->Add(cell_2);
LatticeBin *lpb=new LatticeBin(rod); //rod is a cell associated to an element of a Lattice
lpb1->AddContainer(core); //core is a cell fills by "rods"
t1->Add(lpb);
t1->AddEnergy(1e-2,1.e7); //log energy binning from 0.01 eV to 10 MeV with 10 bin/decade (default)
t1->AddMultiplicator mat_1,new Reaction(102)); //capture reaction rate
t1->AddMultiplicator mat_1,-6); //fission reaction rate
```
When in LatticeBins, you can precise the positions of the cell in the lattice, use a string which gives the position in the MCNP way, that is :

```cpp
string pos="[0 0 0]";
LatticeBin *lpb=new LatticeBin(rod); //rod is a cell associated to an element of a Lattice
lpb->AddContainer(core,pos); //core is a cell fills by "rod" (e.g. a lattice)
```

This means that you are interested by the flux in the cell rod which is in the cell core at the lattice position [0,0,0]. The result in MCNP file will be (given that the number of cell rod is 1 and cell core is 2):

```
F4:N (1<2 [0 0 0])
```

Finally, you can use the universe shorthand (u=1 for instance), simply by giving the universe number, either directly if you know it or by using `Shape::GetUniverse()`:

```cpp
Tally *t1=new Tally();
t1->Add(MainShape->GetUniverse());
```

Please note that, in MCNP, the universe shorthand doesn't work for lattice elements filled with the universe number, that is to say, when a lattice is filled by position, for instance when in MURE we use `Cell::FillLattice(Universe Number, Position)`, the universe is not detected by MCNP. (see example `Tally_ Test.cxx`)

As this shorthand is replaced by **all** the cells filled with the universe, we strongly recommend that you use it only if you are perfectly sure of the number of cells involved in the tally defined in that way.

Each tally must have a definite volume (or surface). The volume is set via the `Tally::SetBinVolume()` or for a given bin, via the `TallyBin::SetBinVolume()` or `TallyBin::SetVolume()`. This may be done manually or automatically: if no bin volume has been provided, the volume is automatically calculated by a stochastic method (see also the note C.1 concerning the volume calculation). For lattices and universe shorthand, as multiple volumes are needed for each bin, the `SetBinVolume()` method precises a partial number.

```cpp
Tally *t1=new Tally(4);
t1->Add(Cell1);
t1->SetBinVolume(0,1.e-6,0); // Set the volume of Cell 1 to 1 cm3
LatticeBin *lpb2=new LatticeBin(rod1);
lpb2->AddContainer(rod2);
lpb2->AddContainer(core);
Tally *t2=new Tally(4);
t2->Add(lpb2); // A Lattice Bin giving somthing like (1 2 < 3)
t2->SetBinVolume(0,1.e-6,0); // Set the volume of Cell 1 in cell 3
t2->SetBinVolume(0,1.e-6,1); // Set the volume of Cell 2 in cell 3
```

The partial number follows the order of the output. (to know exactly the order for lattice bin volumes, see MCNP doc or example `Tally_test.cxx`).

### 5.3.1 Fluence to dose conversion (MCNP Only)

The `Tally::AddFluenceToDoseConversion()` method converts fluence (only for F4 type tally) to dose rate (rem/h) using fluence to dose conversion factors. It can be used only for photons or neutrons. The conversion factor are stored in `/path_to_mure/MURE/data/GammaFluenceToDose.dat` and `NeutronFluenceToDose.dat` (they are taken from [22]).
Interpolation is done for energies and for factors. This interpolation is done by default in log way. The 2 boolean
parameters of `Tally::AddFluenceToDoseConversion()` allows to specify linear interpolation for energy (the 1st one)
and conversion factors (the 2nd one).

```cpp
Tally *MyTally=new Tally(4,ParticleType); //fluence to dose works only for ParticleType='P' or 'N'
MyTally->Add(...);
MyTally->AddEnergy(...);
MyTally->AddFluenceToDoseConversion();
```
# Chapter 6

## Nuclei Tree

### 6.1 Nuclei Trees: general considerations

To evolve a given material, it is essential to know in advance the ensemble of nuclei which it can produce via successive nuclear reactions and decays. This information comes from the *Nuclei Tree* object, which is defined by the **NucleiTree** class, and contains the linking information between all the nuclei that exist in the chart of nuclides. Nuclei can be transformed into other nuclei by nuclear decays (β, α, electron capture, etc.) and nuclear reactions ((n, γ), (n, 2n), fission, etc.). In general a particular nucleus has daughters (those nuclei it can transform into) and parents (those nuclei that can produce it). To calculate how much of a particular nucleus is produced during the evolution it is essential to know the ways in which it can be produced and the ways it can be destroyed.

We cannot only define the total global nuclear tree (the connections between all 3834 nuclei), but also the local or particular tree for a given nucleus. The particular tree for a certain nucleus is, in general, a subset of nuclei and their linking from the global tree. To extract an individual tree from the global tree, the *NucleiTree::ExtractZAI* method is used. The local tree is explored recursively by starting at the given nucleus (see figure 6.1) and then finding each successive reaction and decay daughters.

#### 6.1.1 How the tree information is stored

* The global nuclear tree has a vector of all the physically possible ZAI objects, representing the different isotopes, including long-lived isomers of the same isotope, as an attribute. Each ZAI object contains vectors of pointers to other ZAI which are either its reaction daughters or its decay daughters. When a particular tree is extracted, the parental relations are also determined (see figure 6.2). It is then possible to follow the tree by jumping from ZAI to ZAI, and to understand for each ZAI object all the ways it can be created and destroyed.

#### 6.1.2 Cutting or simplifying the tree

Isotopes in the chart of nuclides have half-lives ranging over ~40 orders of magnitude. In a reactor core neutron fluxes are much lower (~10¹⁵ neutron/cm²/s) than in supernova explosions (up to ~10²⁴ neutrons/cm²/s) and thus decay rates are usually much higher than reactions rates for most short-lived nuclei. Furthermore, the reactor fuel composition only changes significantly after some years, so short-lived isotopes play no important role in the evolution. We can therefore assume that nuclei with short half-lives decay instantaneously which allows us to simplify the tree
<figure><img src="image_14.png" /><figcaption>Figure 6.1: The chart of nuclides used in MURE containing 3834 Nuclei. Nuclei colored in red are those present in the Nuclear Tree of $^{232}$Th. To build this tree, simple rules are used for determining which reactions are possible: if $Z < 50$, $(n, \gamma)$ reactions are allowed, if $Z < 90$, $(n, \gamma)$ and $(n, 2n)$ reactions are allowed, and if $Z \ge 90$, $(n, \gamma)$ and $(n, 2n)$ fission reactions are allowed.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 6.2: The general case of a ZAI object X in a tree with multiple reaction/decay parents and daughters.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 6.3: Tree simplification if nucleus X is produced via radioactive decay</figcaption></figure>

and reduce the number of nuclei in the calculation and hence the calculation time. To simplify a nuclear tree requires these short-lived nuclei to be cut out of the tree and their reaction/decay parents and daughters re-linked into the tree correctly. There are two types of simplification for a particular tree branch, depending on whether the nucleus X in question is produced via decay (see figure 6.3) or nuclear reaction (see figure 6.4).

Note that this means that in some cases a single type of reaction can produce multiple daughters with different probabilities (branching ratios), since the reaction daughter had more than one decay daughter itself, and was cut out of the tree because its half-life was too short (see figure 6.4). After the simplification process, these reaction daughters will have weights which are proportional to the decay branching ratios of the nucleus that was removed from the tree.

For a particular tree many nuclei towards the neutron drip-line end up getting cut out due to their progressively shorter half-lives. In figure 6.5 we can see the nuclei explored and then removed from the tree right up to the edge of the nuclear chart.

### 6.1.3 The difference/similarity of material trees

Two nuclei X, and Y will have identical trees if there is a pathway of decays and reactions that produces Y from X and also a pathway of decays and/or reactions that produces X from Y. Reactions (n, γ) and (n, 2n) move up/down in mass, and decays only move down in mass. Once the mass is sufficiently low that (n, 2n) reactions are not permitted, then decays will only take nuclei as far as the valley of stability and a stable nucleus are reached. At this point the tree branch ends. However, successive neutron captures are always possible, therefore fuels containing lighter elements, such as the Oxygen in Uranium Oxide (UOX) fuel, will have a different tree from a pure fuel's one such as the <sup>232</sup>Th (see figure 6.6).
<figure><img src="image_17.png" /><figcaption>b) After Simplification</figcaption></figure>

<figure><img src="image_18.png" /></figure>

Figure 6.4: Tree simplification if nucleus X is produced via a reaction. The decay daughters of X become reaction daughters of A.
<figure><img src="image_19.png" /><figcaption>Figure 6.5: The simplified nuclear tree for $^{232}$Th. Nuclei in yellow are those which were cut out during the tree exploration process because their half-lives were too short (less than 1 hour). Nuclei in red remain in the tree. Reactions follow the same rules that those defined in fig. 6.1.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 6.6: The UOX tree for $T_{1/2} \ge 1$ hour. Reactions follow the same rules that those defined in fig. 6.1.</figcaption></figure>

### 6.1.4 How many nuclei need to be considered in the evolution?

In all the previous chart figures, the reactions taken into account are: if $Z < 50$, $(n, \gamma)$, if $Z < 90$, $(n, \gamma)$ and $(n, 2n)$, and if $Z \ge 90$, $(n, \gamma)$ and $(n, 2n)$ are fission reactions. In reality, of course, nuclear cross sections are only available for a subset of all nuclei. Figures 6.7 and 6.8 show the nuclei with available reaction data found in the default MCNP4 cross-section (from ENDFB base). All the available reactions are plotted in figure 6.9.

For a realistic half-life cut, the actinide nuclear tree still contains several hundred nuclei. However, simplification has reduced considerably the complexity of the problem(see figure 6.10)<fnref n="1" />.

### 6.1.5 Fission Products

Fission products are built according to the fission yields available in nuclear data libraries (ENDF B 6.8). They are the major contributor to nuclei increased in the reaction tree. Thus, the more FP in nuclear libraries, the longer the *MCNP/Serpent* execution time will be due to the increase in tally number. Thus, in order to obtain a first approximation, it may be convenient to run only a selection of fission products that may be considered as “fundamental”. We provide a method, *NucleiTree::FissionProductSelection()* to choose the most important FP. All other FP will not evolve (no resolution of Bateman equations) but instead are added in the real materials of the *MC* input file in order to have a good flux description. A default selection has been entered with the nuclei of Table 6.1; nevertheless the user can provide a

<footnote n="1">¹In this figure, we also show the MURE cross-section data base built at IPNO/LPSC; this base stores nuclei from ENDF, JENDL and JEFF. The relatively small difference in nuclei number between MCNP and MURE bases is due to the fact that this tree has been built at 300K. When choosing higher temperatures (as in standard reactors), the difference will increase because the MURE base has been processed for different temperatures (500K, 600K, 800K, 1000K, 1200K, 1500K) for minor actinides and main fission products. This base is not part of the MURE package.</footnote>
<figure><img src="image_21.png" /><figcaption>Figure 6.7: <sup>232</sup>Th tree with MCNP default cross-sections and T<sub>1/2</sub> ≥ 1 hour.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 6.8: The UOX tree with MCNP default cross-sections for T<sub>1/2</sub> ≥ 1 hour.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 6.9: MCNP default available reaction data are shown in light blue.</figcaption></figure>

file with its own selection. To use this selection you must use `gMURE->KeepOnlyFissionProductSelection()`; to give your own selection, use `gMURE->KeepOnlyFissionProductSelection("MyFPSelection.dat")` where the file `MyFPSelection.dat` is an ASCII file with each FP you prefer to keep for evolution: for example

```
Z1 A1
Z2 A2
...
```

<table><caption>Table 6.1: Default 58 fission product selection in `NucleiTree::FissionProductSelection()`.</caption><tr><td>36-Kr- 83</td><td>46-Pd-107</td><td>54-Xe-134</td><td>60-Nd-145</td><td>62-Sm-152</td><td>66-Dy-160</td></tr><tr><td>40-Zr- 93</td><td>46-Pd-108</td><td>54-Xe-135</td><td>61-Pm-147</td><td>63-Eu-151</td><td>66-Dy-161</td></tr><tr><td>42-Mo- 95</td><td>47-Ag-109</td><td>54-Xe-136</td><td>61-Pm-148</td><td>63-Eu-152</td><td>66-Dy-162</td></tr><tr><td>43-Tc- 99</td><td>48-Cd-113</td><td>55-Cs-133</td><td>61-Pm-149</td><td>63-Eu-153</td><td>66-Dy-163</td></tr><tr><td>44-Ru-101</td><td>49-In-115</td><td>55-Cs-134</td><td>62-Sm-147</td><td>63-Eu-154</td><td>66-Dy-164</td></tr><tr><td>44-Ru-103</td><td>51-Sb-125</td><td>55-Cs-137</td><td>62-Sm-148</td><td>63-Eu-155</td><td>67-Ho-165</td></tr><tr><td>44-Ru-106</td><td>52-Te-127</td><td>55-Cs-135</td><td>62-Sm-149</td><td>64-Gd-155</td><td>68-Er-166</td></tr><tr><td>45-Rh-103</td><td>53-I-127</td><td>56-Ba-138</td><td>62-Sm-150</td><td>64-Gd-156</td><td>68-Er-167</td></tr><tr><td>45-Rh-105</td><td>53-I-135</td><td>59-Pr-141</td><td>62-Sm-151</td><td>64-Gd-157</td><td></td></tr><tr><td>46-Pd-105</td><td>54-Xe-131</td><td>60-Nd-143</td><td>64-Gd-154</td><td>64-Gd-158</td><td></td></tr></table>
<figure><img src="image_24.png" /><figcaption>Figure 6.10: Number of nuclei in the <sup>232</sup>Th tree as a function of the minimum half life allowed, T<sub>1/2min</sub>. The blue curve shows the number of nuclei in the tree if the simple rules for allowed reactions are used (see figure. 6.1). The black curve shows the number of nuclei in the tree if allowed reactions are taken from the *MCNP* default reaction data base. The green curve corresponds to the *MURE* data base built at IPNO/LPSC (not part of the *MURE* package). The red dot shows the number of nuclei in the tree for the default cut (1 hour).</figcaption></figure>
### 6.1.6 Isomer Production from (n,gamma) or (n,2n) reactions

**TO BE NOTICED:** Before July 2012, MURE only takes into account the case of $^{241}$Am. In the previous treatment, all ground state production was replaced by $^{242}$Cm; this was of course wrong and lead to an over estimation of $^{242}$Cm production and an under estimation of $^{242}$Pu. This has been corrected since July 2012.

In some cases, $(n, \gamma)$ (or $(n, 2n)$) reactions can lead to either ground state and isomeric state of a nucleus. For example, let take the case of $^{241}$Am. Neutron capture on this nucleus leads to $^{242}$Am($T_{1/2} \sim 16h$) and $^{242m}$Am(141y). Ground state decays in $^{242}$Pu(17.3%) by electronic capture and in $^{242}$Cm(82.7%) by $\beta$-decay. A file (*IsomerProduction.dat*) is provided to give the branching ratio (for thermal and fast spectrum) to produced its ground and isomeric states of some nuclei. The format of this file is

```
Z A ReactionCode BRth BRf halflife CutFlag
...
```

where Z, A, ReactionCode, BRth and BRf are respectively the Z (95 for $^{241}$Am), A (241 for $^{241}$Am), the reaction type (allowed type are $n\_gamma$ or $n\_2n$) and the thermal (87.33% for $^{241}$Am) and fast (85% for $^{241}$Am) branching ratio of the $(n, \gamma)$ reaction that leads to ground state (of $^{242}$Am). The “halflife” parameter is the half life of the produced ground state ($T_{1/2} \sim 16h$ for $^{242}$Am) follows by its unit (here “h”) ; this parameter is not used, it is just an indication that can help the user to decide if the last parameter (CutFlag) is a “X” or a “M”. “M” means that MURE decide what to do according to MURE::GetShortestHalfLife() method. “X” forced the replacement of the ground state (e.g. $^{242}$Am) by its decay daughters. To be noticed, if “M” is chosen, wrong results can be obtained in $^{242}$Pu and $^{242}$Cm for example, if no cross section is available for the ground state ($^{242}$Am).

**Branching Ratio** Evaluated branching ratio data for isomer production can be obtained at NNDC-ENDF, in the “advanced” or “extended” retrieval tab, by choosing “n,g” in the Reaction line and “MRNP” ( multiplicities for production of radioactive elements, MF=9) in the Quantity line. The branching ratio given in MURE/data/IsomerProduction.dat are effective branching ratio computed on a typical thermal and fast spectrum as explained in [25]:

$$BR_{eff} = \frac{\int_{0}^{\infty} BR(E)\sigma_{react}(E)\phi(E)dE}{\int_{0}^{\infty} \sigma_{react}(E)\phi(E)dE}$$

The NNDC-ENDF site gives either the branching ratio to the ground state or to the isomeric state(s); but in MURE the branching ratio given is always to the ground state.

## 6.2 Nuclei Tree: the implementation

### 6.2.1 Important files

All provided data files are located in MURE/data. One can use either, gMURE->SetDATADIR(path) or the environment variable DATADIR (setenv DATADIR path in csh).

* When building trees, nuclear mass data are also incorporated: the two files used for the input are Mass.dat (a modified version of the "The 2003 Atomic Mass Evaluation" data from the National Nuclear Data Center (BNL)) and NaturalIsotopeMass.dat that is used for non evolving natural isotopes (i.e. with A=0); this file has been built from the NNCD nuclear wallet card. The names of these files can be defined by gMURE->SetMassDataFileName() and gMURE->SetNaturalIsotopeMassFileName() respectively.
*   **Decay data**: The decay data for the ~ 4000 nuclei in the nuclear chart comes from the *chart.jeff3.1.1* file which has been constructed from the ref. [16]. The name of this file is set through the command `gMURE->SetNucleiChartFileName()`.

*   **Reaction data**: Allowed reactions are defined in **ReactionList** class. This class allows different kinds of possibilities for allowed reaction initialization. The reactions allowed by using the simple rules (see Fig. 6.1) are defined in **ReactionList::InitJon()**; the automatic detection from available cross-sections is defined in **ReactionList::InitPTO()** (see next paragraph for more details). It uses cross-sections in **ACE format**. The **AvailableReactionChart.dat** helps to save times (this file informs if cross-section are available or not for each nucleus of *chart.jeff3.1.1* file). This file could be made by using *MURE/util/fp/CheckReaction*. Other initialization functions are free to be defined by the user. This is particularly useful for prohibiting all reactions except a few of particular interest to the given problem.

One can choose a given method to generate *ReactionList* with

```R
gMURE->SetReactionListInitMethod(new TSpecificFunctor<ReactionList>(0, &ReactionList::InitMethod));
```

where `InitMethod` could be `InitJon`, `InitPTO`, ..., `ReactionList::InitPTO()` being the default.

*   **Isomer production**: a file *IsomerProduction.dat* is provided in *MURE/data* that contains information to produce isomer state and ground state from (n, γ) reactions. This file contains Z,A, thermal and fast branching ratio for a mother nucleus (e.g. $^{241}$Am) as well as half-life of the produced ground state (e.g. $^{242}$Am) and a string flag (either “M” or “X”) to know if user want to replace the ground state by its daughters (e.g. electronic capture $^{242}$Pu and β⁻-decay $^{242}$Cm daughters). “M” means that *MURE* decide what to do according to *MURE::GetShortestHalfLife()* method. “X” forced the replacement of the ground state (e.g. $^{242}$Am) by its decay daughters. **If a file *IsomerProduction.dat* exists in the local run directory, this file is used instead of the *MURE/data/IsomerProduction.dat* one.**

*   **Fission**: Two files are needed for the inclusion of spontaneous fission decays and neutron induced fission reactions: *FPavailable.dat* and *FPyield.bin*, which have been built by *MURE/util/fp/GenerateFPYield*. The first file is in ASCII format and contains the fission energies for each available ZAI, and the position (records) of the fission product yields in the binary file *FPyield.bin* (see for example Fig. 6.11 and 6.12). Even if these files are provided with *MURE*, you can regenerate them by providing the appropriate ENDF fission yield file.

### 6.2.2 Reaction Auto-detection

Reaction auto-detection is based on information found in the *ACE* binary<fnref n="2" /> <fnref n="3" /> format files. Nuclei, for which reaction cross-section database information is available, can be found. The auto-detection process relies on the **ReactionList** class.

<footnote n="2">²One can make ACE binary format from ASCII ACE format using MCNP *makxsf* exec. Even if NJOY is theoretically able to create such files (ACE binary files), our experience shows that only ASCII files made by NJOY are usable. So you can process ACE ASCII files from ENDF format with NJOY and then, use *makxsf* to convert them to binary.</footnote>

<footnote n="3">³An ASCII file may be used but since the run time increases significantly, we recommend to convert the ASCII ACE file to a binary ACE file with *makxsf*.</footnote>
<figure><img src="image_25.png" /><figcaption>Figure 6.11: Fission yield for <sup>235</sup>U.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 6.12: Fission yield for <sup>235</sup>U with T<sub>1/2</sub> > 1 hour.</figcaption></figure>
### How does it work?

* First, try to find out if auto-detection has been performed for the given ZAI by reading the binary file "ReactionList/Z.bin" where Z in the proton number of the ZAI. This file contains an array of Boolean for each reaction.

* If ZAI is not found in "ReactionList/Z.bin", start auto-detection:

- Try to find in BaseSummary.dat file the ZAI (or an isomeric state of the ZAI)

- If it is found, read the ACE file and, among all reactions found in this file, find those which are significant according to a given threshold (see next paragraph)

- Allow significant reactions and disable all others in ReactionList and add the ZAI to the "ReactionList/Z.bin" file with its reactions.

**PLEASE NOTE:** if the BaseSummary.dat file is modified (add new nucleus, or new base) or if the threshold used to find significant reactions is changed, **THE ReactionList DIRECTORY MUST BE REMOVED**. Otherwise, your modifications will not be taken into account.

### What are significant reactions?

A quantity $\Xi_i$ is calculated as

$$ \Xi_i = \frac{\int \sigma_i d \log E}{\int d \log E} $$

where i stands for reaction i. Then if and only if $\Xi_i \ge \sigma_{min}$ where $\sigma_{min}$ is a minimum threshold, the reaction is taken into account; The default threshold has been set to 0.01 barn$^4$. This can be modified by MURE::SetReactionThreshold() called before any MURE::SetReactionListInitMethod() or Material definition. In order to save time, the procedure to allow or disable a reaction for a given nucleus is performed only the first time an evolution for that nucleus is asked; then the result for that nucleus (i.e., a ReactionList object) is stored in a file (one per Z) in a specific local directory named ReactionList.

**NOTE** Any change of the MURE::SetReactionListInitMethod() or MURE::SetReactionThreshold() after a first run will be taken into account only if the nucleus has not been already processed: thus you must destroy the ReactionList directory in order to have these changes taken into account.

### 6.2.3 The recursion depth cut

It is possible to put a limit on the recursion depth of the reactions when calling the NucleiTree::ExtractZAI method. For example prohibiting 30 or more successive neutron captures limits the nuclei included in the tree to the top of the actinide region (e.g Fermium, Californium, etc.) where experimental cross-section information does not exist, and many successive neutron captures are the only way in which these nuclei can be produced. The recursion depth parameter should be used with care. Figure 6.13 shows the effect of the recursion depth. The default value is 10000...

<footnote>$^4$This threshold has been chosen in order to suppress "exotic" reactions, such as (n,α), as well as other insignificant reactions and to keep "essential" reactions; a file, SuppressReaction.dat in the ReactionList directory contains the removed reactions. Don't forget to verify it!!!</footnote>
<figure><img src="image_27.png" /><figcaption>Figure 6.13: Effect from limiting the recursion depth for successive reactions in the <sup>232</sup>Th tree, where nuclei with T<sub>1/2</sub> < 1 hour have been removed.</figcaption></figure>
# Chapter 7

## Evolution in MURE

### 7.1 Preliminary Remark

Before doing any evolution you must provide a system (geometry based on Cells, all Materials, running MC code mode (critical or subcritical), ...). This system definition is normally provided by the MURE geometry input file (a C++ code). But if you have a very large MCNP file already defined, it is possible to perform an evolution using this file without redefining the whole MCNP geometry in MURE style. The way of doing this is explained in section 7.5. This has not been yet implemeted for Serpent.

### 7.2 General Considerations

Simulating time evolution involves solving the well-known Bateman equations for every given nucleus in every given cell for which evolution is desired, and carrying out this integration over a discrete number of time steps. Simulating the fuel burn-up may also involve adjusting the neutron flux value at every step of the integration process so as to model constant reactor power, adjusting poison concentration to keep a given value for $k_{eff}$, ...

Implementation of evolution related methods are handled via the **EvolutionSolver** class. An object of that class is automatically built in **MURE** constructor, thus can access to it via `gMURE->GetEvolutionSolver()`.

The general scheme for the evolution is shown in figure 7.1:

* The nuclei tree is built once and for all at the first MC run; it is controlled via the `Material::SetEvolution()` method which must be called for each evolving *Material* before any cell definitions are made. Initial compositions of all materials are entered by the user and these will evolve automatically.

<figure><img src="image_28.png" /><figcaption>Figure 7.1: Evolution principle.</figcaption></figure>
* All the necessary tallies for calculating mean neutron fluxes and cross-sections in evolving cells are automatically built.

* The MC input file with the composition at a given time $t_i$ is built and a MC run is performed.

* Bateman equations are solved for (Runge-Kutta method) using fluxes and cross-sections of the MC run over a given $\Delta t_i$

* A new MC file with the composition at $t_{i+1} = t_i + \Delta t_i$ is built, and so on.

At the present time:

* User must provide all $t_i$ in a vector of double where the ith element is the value of instant $t_i$ at which the Monte-Carlo transport is to be performed.

* EvolutionControl is used to define evolution conditions (constant power, $k_{eff}$ adjustments, ...). Using inheritance, you may develop your own EvolutionControl class (see for example the PoisonEvolutionControl class and § 7.4).

The results of an evolution consist of a number of files stored in an “evolution” directory; the name of this directory should be defined by the user via the MURE::SetMCRunDirectory function. Two general sorts of files can be found in this directory: On the one hand, there are the MC input and output files (“o”, “m”, “s” and so on for MCNP and “_res.m”, “_det0.m”, ... for Serpent) and on the other, the MURE-generated summary files (DATA_i, BDATA_i, ...) that can be read by a ROOT (see ROOT home page) C++ interface: MureGui.

### 7.2.1 Time discretization

There are 3 levels of time discretization:

* The first level is the MC step, i.e., the number of $\Delta t_i$ or MC runs. These steps are user-defined and may be unregular.

* The second level is the discretization within a $\Delta t_i$. Each $\Delta t_i$ is divided in $N_{RK}$ equally spaced Runge-Kutta $\delta t$ steps (which we'll also refer to as RK steps). At each $t_k = k\delta t$, Bateman equations are built with given cross-sections. Special methods controlling the evolution could be called at these times (e.g. the flux is renormalized to keep the constant power, ...).

* Then the last discretization step is performed automatically by the adaptive step size RK method. During this step (discretization of $\delta t$ in $dt_{rk}$) cross-sections as well as fluxes are kept constant.

### 7.2.2 The Bateman's Equation

Here is the general form of the first-order differential equation that has to be solved for every nucleus in every evolving cell.

$$ \frac{\partial N_i}{\partial t} = -\lambda_i N_i + \underbrace{\sum_j \lambda_j^{j \to i} N_j}_{\text{Decay}} + \underbrace{\sum_j N_j \sigma_j^{j \to i} \langle \phi \rangle - N_i \sum_{\forall r} \sigma_i^{(r)} \langle \phi \rangle}_{\text{Reaction}} $$
The writing and solving of these equations are systematically carried out by using internal methods of the *EvolutiveSystem* class. Only the general approach will be specified here. This approach involves the definition of what is called the *evolution matrix E* for every evolving cell. The user does not need to know in detail how it works and it is only being described here for the record.

### 7.2.2.1 The evolution matrix

Internally, every evolving cell knows which nuclei it contains (initially) or will eventually contain due to the evolution of its initial composition, thanks to the prior construction of the nuclei tree (c.f. *NucleiTree* in chapter 6). The approach taken in *MURE* is to assign unique indexes to every nucleus in the cell's composition so as to build its evolution matrix *E*. It is a square matrix of *n* dimensions, *n* being the total number of nuclei in the composition of the cell.

After completion of the matrix, the ith nucleus of a given cell will evolve according to its own Bateman's equation which will then be simply written as the sum over all nuclei:

$$ \frac{\partial N_i}{\partial t} = \sum_j E_{ij} \times N_j $$

where $E_{ij}$ is the corresponding element of matrix *E* and represents the contribution of the jth nucleus in the cell's composition to the rate of change in time of the ith nucleus.

The actual writing of the matrix follows the natural structure of Bateman's equations and is done in two separate parts : first by taking into account all parent/daughter decay relations between nuclei considering proper decay branching ratios. After this, the parent/reaction daughters term is added, where, in the case of fission, fission yields are properly considered.

The evolution matrix is rebuilt in this way at every Runge-Kutta step (c.f : Time Discretization) for all evolving cells. The solving for the material's composition new proportion vector $\vec{N} = \{N_0, N_1, ..., N_n\}$ is done automatically using a fourth-order Runge-Kutta-type integration method.

### 7.2.2.2 Multi-Threading parallelization

Since *gcc* version 4.1, *OpenMp* multi-threading is available with *gcc* compilation. Depending on the distribution, *libgomp* library (*omp.h* and *libgomp*) need to be installed. The multi-threading behavior is controlled by the *OMP_NUM env*ironment variable. Setting it to 1 is equivalent to the single processor mode. Setting it to the maximum number of cores available, will run on all that cores (which is the default when this variable is not set). This variable only controls the behavior of the evolution NOT the *MC* run itself (which is still controlled by the *MURE::SetOMP()* method).

### 7.2.2.3 Cooling period

Whereas in general the evolution is useful when the reactor runs at its given power, it is also possible to introduce some cooling period, where the power (and thus the flux) is set to 0. This is done by using an array of booleans of same size than the time vector at which *MC* runs are normally performed. Each time this array contains a "true", the flux is imposed to 0 for the considered step, else the evolution is a standard one. The array is passed to *MURE* via *gMURE->SetCooling()* method.
## 7.3 Different way of evolution

We suppose that the evolution is done at a constant power. Evolution can be performed in different ways that can influence the result.

1. **Evolution with a "constant" reaction rate:** this is the default. Between 2 *MC* steps, reaction rates are kept "constant"; a new reaction rate evaluation is performed at the next *MC* run for the next step. In fact, they are not really constant: due to renormalization of the flux to keep constant power, reaction rates used in Bateman equations are in fact, $\alpha\sigma\phi$ where $\alpha$ is the flux renormalization coefficient. Thus, periodically within a *MC* step, these reaction rates are readjusted (this explains the use of the term "constant" in quotes). An example of such evolution can be found in *MURE/examples/Evolution/EvolvingSphere.cxx* or *MURE/examples/Evolution/EvolvingSphere_serpent.cxx*

2. **Cross-section evolution:** this is no longer the default. Previous *MC* runs are used to linearly extrapolate the reaction rates for the next step.

3. **Window average evolution:** a mean ($\sigma\phi$) is calculated from the previous *MC* runs and the actual one; this mean value is used to solve Bateman's equations. The window to compute average is based on *Evolution-Solver::GetFitRangeNumber()* previous runs (default is 4).

4. **Predictor-corrector methods.** There are 3 Predictor-corrector (PC) methods implemented. This is actually new and not fully tested (despite the fact that it currently works). Suppose that the *MC* run at time $T_N$ has been done (leading to reaction rates ($\sigma\phi)_N$). The evolution has to be performed up to $T_{N+1}$.

(a) **PC at middle point:** the predictor run is performed with "constant" reaction rates ($\sigma\phi)_N$ up to $T_C = \frac{T_N+T_{N+1}}{2}$; then a corrector *MC* run is performed and leads to ($\sigma\phi)_C$. Starting with the composition of $T_N$, the evolution is done again from $T_N$ to $T_{N+1}$ with the "constant" reaction rates ($\sigma\phi)_C$ evaluated at the middle point.

(b) **PC at end-point:** the predictor run is performed with "constant" reaction rates ($\sigma\phi)_N$ up to $T_C = T_{N+1}$ leading a composition $C_{N+1}^P$; then a corrector *MC* run is performed and leads to ($\sigma\phi)_{N+1}$. Starting with the composition of $T_N$, the evolution is done again from $T_N$ to $T_{N+1}$ with the "constant" reaction rates ($\sigma\phi)_{N+1}$ leading to a composition $C_{N+1}^C$. Then, the composition used for the next predictor run (for $T_{N+1}$ to $T_{N+2}$) is $C_{N+1} = \frac{C_{N+1}^F+C_{N+1}^C}{2}$.

(c) **Kind of PC at end-point:** the predictor run is done like in b) but for the corrector run, we used a linear interpolation of reaction rates from ($\sigma\phi)_N$ to ($\sigma\phi)_{N+1}$. The next composition for $T_{N+1}$ to $T_{N+2}$ is $C_{N+1} = C_{N+1}^C$.

First tests seem to give an advantage to method b) for the quality of its results (compositions, dispersion of composition for *M* identical run with different random seeds, ...).

5. **Multi-Group evolution**[17]. If this method is used, not that reaction rate tallies are not built in the *MC* input file; instead, a thin energy binning flux will be built for each evolving cell. Then, after the *MC* run, reaction rates are calculated as $\int \phi(E)\sigma_i^r(E)dE$ where the sum is done over each energy group, $\sigma_i^r$ is the pointwise cross-section of reaction *r* for nucleus *i* read in the ACE file used for *MC* run (if ASCII
version is used, please note that this will take longer...). The main advantage of this method is a considerable CPU time saving (more than 30 times faster<fnref n="1" /> than the “standard” evolution). The accuracy of the result are quite good (~1% for thermal systems and 5% for fast ones<fnref n="2" />). To use this method, one has to define `"gMURE->GetEvolutionSolver()->UseMultiGroupTallies();"` at the beginning of the input file (see `MURE/examples/Evolution/EvolvingSphere2.cxx` example). The default energy binning for flux is (17900 groups)

<table><thead><tr><th>Energy range</th><th>10<sup>-4</sup>eV - 1eV</th><th>1eV - 10eV</th><th>10eV - 10keV</th><th>10keV - 0.1MeV</th><th>0.1MeV - 20MeV</th></tr></thead><tbody><tr><td>Number of bins/decade</td><td>100</td><td>500</td><td>5000</td><td>1000</td><td>500</td></tr></tbody></table>

The main problem of multi-group evolution is the ²³⁸U. Its capture cross section has been chosen to “optimize” the group structure (number of bins/energy interval). In order to take the advantage of the multi-group CPU time saving, it has been implemented the possibility to use multi-group runs for all nuclei but ²³⁸U: in that case, the “standard” method (real MC tally) is used for ²³⁸U, and multi-group tallies for the other nuclei. This is done by setting the `"StdTallyFor238U"` to `"true"` in the `"EvolutionSolver::UseMultiGroupTallies()"` method.

### 7.3.1 Cross-section Evolution

Cross-sections are evolving with time because of the flux shape evolution. Between 2 MCNP steps, cross-sections are either constant, either evolving. Taking the cross-section evolution during the “evolving step” should enable a better precision and increase evolution step duration, reducing the CPU time. It is based on linear extrapolation of previous MC runs.

!!!! After a study of error propagation, it appears that fitting $\sigma\phi$ leads to a bigger dispersion of cross-section evolution in the result (with a constant power). Thus this treatment is no more used by default. To use it you must specify it by a `MURE::FitSigmaPhi(true)`.

#### 7.3.1.1 How does it work?

The treatment is applied for each cross-section of each nucleus in each cell. Thus, for the sake of simplicity, we will only consider a reaction cross-section $\sigma_i$ of a nucleus in a cell.

Since for evolution, reaction rates are used, we should rather consider the reaction rate $\sigma_i\phi$ and not the cross-section $\sigma_i$. Suppose that we have performed an evolution up to step $k$ (i.e., $k$th MC run is just finished). To know the reaction rate evolution up to step $k+1$ (i.e., just before a new MC run), a linear fit (least square method from Numerical Recipes) is applied to `EvolutionSolver::GetFitRangeNumber()` previous MC reaction rates<fnref n="3" /> (that is to say, MC runs $k-3, k-2, k-1$, and $k$). Then, for the evolution over step $k \to k+1$, the reaction rate is extrapolated with $(\sigma_i\phi)(t) = at + b$, where $a$ and $b$ are the results of the linear regression. Note that for the first `EvolutionSolver::GetFitRangeNumber()` MC runs, reaction rates are kept constant.

As already mentioned, integration of Bateman equations is in fact done in $N_{rk}$ Runge-Kutta steps which are automatically divided in other shorter time steps (adaptive step RK method); the reaction rate evolution is only

<footnote n="1">¹When using parallel calculation for MCNP (omp,mpi) this factor is reduced: for example, running MCNP in omp on a 8 processor machine leads to only 16 times faster for a simple problem.</footnote>

<footnote n="2">²This is of the same order of magnitude than the deviation observed in a series of N identical evolution starting with a different random seed. Using this multi-group method, unresolved probability tables are not taken into account for reaction rates evaluations; this explains that for fast systems, the discrepancy is larger than for thermal ones.</footnote>

<footnote n="3">³The default is 4, and it can be changed with `MURE::SetFitRangeNumber()`.</footnote>
<figure><img src="image_29.png" /><figcaption>Figure 7.2: Illustration of a reaction rate evolution. MCNP results are in black dots. The reaction rate used during the evolution from step $k$ to step $k + 1$ is in green whereas the red squares indicate the time where the extrapolation is performed.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Figure 7.3: EvolutionControl class and siblings.</figcaption></figure>

performed at each RK step level and it is kept constant within the current RK step. This constant value is evaluated in the middle of the RK step to be more accurate (see fig. 7.2).

In order to better consider the rapid variation of a reaction rate slope, we have slightly improved the method as follows (here $N_{fit} = EvolutionSolver:: GetFitRangeNumber()$):

* a linear regression is performed on $N_{fit}$ to obtain $a'$ and $b'$ for the reaction rate evolution

* a linear regression is performed on $N_{fit} - 1$ to obtain $a''$ and $b''$

* then if the 2 slopes $a'$ and $a''$ are of opposite sign, the slope taken for the evolution is the mean $a = \frac{a'+a''}{2}$ and $b = \frac{b'+b''}{2}$

* or else, $a = a'$ and $b = b'$ if $|a'| < |a''|$ or $a = a''$ and $b = b''$ if $|a''| < |a'|$

## 7.4 Reactivity Control

NOTE: only few tests have been performed and therefore many bugs may be present in the Evolution Control.

As already mentioned, the control of the whole evolution is done via the *EvolutionControl* class. This class already includes some useful methods, but the best thing to do in order to implement your own control is to write a derived class. In a first example, the reactivity control uses the original *EvolutionControl* class; then the way of implementing your own control will be given. The EvolutionControl class inheritance tree is shown on Fig. 7.3
### 7.4.1 A simple example using standard EvolutionControl

The example file **poison.cxx** (MURE/example/poison.cxx) corresponds to a very simple reactor (cylinder of borated water with a hexahedral fuel pin lattice). The reactivity control is done by adjusting the concentration of a poison (the boron) in the water. We assume that a large reserve of boron is available and thus, when necessary, we can adjust the boron concentration at a desired value. The method principle is to calculate an estimation of $k_{eff}$ using the reaction rate extrapolations:

$$k_{eff} = \frac{\nu N_f \sigma_f \phi}{N_a \sigma_{abs w/o poison} \phi + N_p \sigma_{poison} \phi - N \sigma_{n,2n} \phi - 2N \sigma_{n,3n} \phi + escape}$$

(n,xn), "escape" and $\sigma_{poison}$ are supposed to be constant between 2 MC steps; the poison concentration $N_p$ is thus determined using the desired $k_{eff}$ in the formula.

#### 7.4.1.1 Poison material declaration

```cpp
ControlMaterial *BoratedWtr=new Material (0.7187023, 673);
BoratedWtr->SetDensity(0.719);
BoratedWtr->AddTheControlNucleus(5, 10, 0.00012);
BoratedWtr->AddTheControlNucleus(5,11,0.00048);
BoratedWtr->AddNucleus(8,16,0.333);
BoratedWtr->AddNucleus(1,1,0.667, "H2O");
```

A *ControlMaterial* is a *Material* which is used to control the evolution, and in this case, a poison to control the reactivity. The 3th and 4th lines informs us that boron 10 (and 11) are the nuclei used to control the reactivity. See the "poison.cxx" example.

This material must not be an evolving material, because we'll consider that we can adjust boron concentration every time.

#### 7.4.1.2 Escape calculation

In the $k_{eff}$ formula, neutron leakage is needed; in order to calculate these escapes, one has to define a F2 tally (flux through surface) through the most outer Shape of the geometry; in the example this leads to:

```cpp
gMURE->SetOutermostShape(Vessel);
```

because "Vessel" is the outer cylinder of this "reactor".

#### 7.4.1.3 Evolution definition

You have to declare the Evolution Control object:

```cpp
PoisonEvolutionControl *ECon=new PoisonEvolutionControl();
ECon->AddReactions(BoratedWtr);
gMURE->SetEvolutionControl(ECon);
```

The 2nd line gives the specific EvolutionControl command to add the required reactions useful for the control (poison absorption, ...)
### 7.4.2 Using your own EvolutionControl

The aim of this section is to help user implement its own `EvolutionControl` using inheritance mechanism. Then you will be able to write your own specific methods that will be called for the evolution control. Supposed we wish to create a `MyEvolutionControl` class: the `hxx` and `cxx` files are in `MURE/example/`.

NOTE: Don't forget to redefine the `EvolutionControl* Clone()` as:

```c
EvolutionControl* Clone() {return new MyEvolutionControl (*this);}
```

Each method and variable of `EvolutionControl` are known in `MyEvolutionControl`; but because of the virtuality of some methods, the `ControlAtEachMCNPStep()`, `ControlAtEachRKStep()` and `ControlAfterEndOfRKIntegration()` methods of `MyEvolutionControl` will be called instead of the ones of `EvolutionControl`. This means that these original methods of `EvolutionControl` are no more called, and thus you have to read them carefully to implement what is needed or necessary for your evolution control.

Now to use this new class in `Myprog.cxx`:

```cpp
#include <iostream>
#include <cmath>
using namespace std;
#include "MureHeaders.hxx"
#include "MyEvolutionControl.hxx"
...
int main(int argc, char** argv)
{
...
MyEvolutionControl *NewEC=new MyEvolutionControl();
gMURE->SetEvolutionControl(NewEC);
}
```

and to compile the Makefile for the compilation will be

```
g++ -g -O -Wall -fPIC -I../..source/include -I../..source/external -c MyEvolutionControl.cxx
g++ -g -O -o Myprog Myprog.cxx MyEvolutionControl.o -I../..source/include -I../..source/external/\ -L../..lib -lMUREpkg -lvalerr -lmctal
```

## 7.5 Evolution of a MCNP user defined geometry

This section is devoted to users who have already defined a complex MCNP file and would rather not use the MURE geometry capability when making evolution. Even if this task requires some work, the procedure is simple. User has to redefine in MURE style, **ALL materials** of its geometry and **all EVOLVING cells** (but not the shapes associated to the cells). Of course, user will need all files required by a standard evolution (BaseSummary.dat, AvailableReactionChart.dat, ...)

The principle is the following: a *line by line* of a MCNP user input file are copied, but

* if the cell number found corresponds to one of the "virtual" evolving cells defined, then the density of the evolving material at the given evolution step is taken into account. Be sure that 2 evolving cells must have 2 different materials (because they usually evolve differently even if at the beginning they are the same); thus, it is at present not possible to make the evolution of a cell defined by "like but" because there is no material on this cell line.
* The general MCNP block will be read, but all the material description ("M" and "MT" cards) will be skipped and replaced by the MURE material definition. Note that if one wishes to use source MURE definitions, particle MODE definitions or PRDMP MURE definitions, any "SDEF", "KSRC", and "KCODE" (for source) or "MODE" or "PRDMP" must be removed from the original file. If you have defined tallies, they will be copied (line by line) and the necessary tallies for evolution will be constructed starting from your highest tally number : for example, if you defined 2 tallies F34 and F44, the automatic tallies for evolution will start at F54 ; tally numbers 4,14 and 24, are not used.

For material declaration, you must use the special *Material* constructor that takes as a single argument the material number ; this number is the one you have written in the MCNP input file. Don't forget to declare ALL materials of the geometry (evolving or non evolving). For non evolving materials, you must ask *MURE* to include them in the *MCNP* output file by using *Material::AddToGlobalNucleiVector()* after the last input nucleus.

For cell declaration, you must use the special *Cell* constructor ; it takes the cell number (the one you have written in the *MCNP* input file), an evolving *Material* pointer and the cell volume (in $m^3$). **ONLY evolving cells should be declared.**

The Original MCNP input name is defined by `MURE::SetUserGeometryInputFile()` method. You can also defined all desired cards of the general data block (source, physical card, ...). The rest of the evolution is completely analog to what is explained in the above sections. Don't forget to provide a mean to calculate the tally normalization factor (you may either give a predefined power or use the `MURE::SetTallyNormalizationFactor()`).

An example is shown in `UserGeo.cxx` with the `UserMCNPGeo` MCNP input file.

## 7.6 More complex evolution conditions

Sometimes, the evolution is performed with power variations, some parameters such as boron concentration, temperature, ... are changing. This can be achieved through the **EvolutionWrapper** class (see *EvolutionWrapper* class). Such conditions are illustrated on the following example:

```
EvolutionWrapper *EC=new EvolutionWrapper(); // create EvolutionWrapper
EC->AddPhase(100,3,1); // create a new 100 days phase with three MCNP steps, (1=logarithmic discretization)
EC->SetPowerLinear(0,4e4); // raise power from 0 to 40 kW
EC->SetMaterialBoronLinear(Cell_of_BoratedWater->GetMaterial(),1000e-6,500e-6); //boron depletion from 1000 to 500 ppm
EC->AddPhase(100,1,0); // a new single-step 100 days cooling phase
EC->SetPowerCooling();
EC->AddPhase(100,3,0); // a new 100 days phase with three equal(0=linear discretization) MCNP steps
EC->SetPowerConstant(4e4); // power is kept constant at 40kW in that phase
EC->Evolve(); // run the evolution (DO NOT CALL the MURE:Evolution method ; it is done here)
```

## 7.7 Equilibrium of Xe-135

In some circumstances, flux oscillations in weak coupled systems create numerical (unphysical) Xe oscillations. These Xe oscillations enhanced the flux oscillations causing numerical instabilities. To reduce these instabilities, a special treatment has been implemented. To use this treatment, one has to call the `EvolutionSolver::SetXe135Equilibrium(double teq)` method, where $t_{eq}$ is the time from which the treatment is really applied (default time value is $3 \times T_{1/2}^{135}X_e \sim 27h$).

Here is the description of what is done:
• if $t < t_{eq}$, then the standard evolution is carried out (i.e. no special treatment)

• if $t \ge t_{eq}$, then

– the $^{135}$Xe Bateman equation is imposed to be at equilibrium : $\frac{dN_{Xe}}{dt} = 0$

– the $^{135}$Xe is calculated as

$$N_{Xe} = \frac{y_{Xe} \Sigma_{fission} \tilde{\phi} + \lambda_{135I} N_{135I}}{\lambda_{135Xe} + \sigma_{n,\gamma}^{Xe} \tilde{\phi}} \quad (7.1)$$

where $y_{Xe}$ is the fission yield for $^{135}$Xe, $\Sigma_{fission} \tilde{\phi}$ is the number of fissions, $\lambda_i$ is the decay constant for nuclei i and $\tilde{\phi}$ is a “mean” flux in the cell calculated as

$$\tilde{\phi} = \frac{\phi^{\text{previous MCNP run}} + \phi^{\text{current MCNP run}}}{2}$$

if there is a factor greater than 1.5 between $\phi^{\text{previous MCNP run}}$ and $\phi^{\text{current MCNP run}}$, else $\tilde{\phi} = \phi^{\text{current MCNP run}}$

– the $^{135}$I used in (7.1) is calculated in the same way :

$$N_I = \frac{y_I \Sigma_{fission} \tilde{\phi} + \lambda_{135Te} N_{135Te}}{\lambda_{135I} + \sigma_{n,\gamma}^{I} \tilde{\phi}}$$

but this value is only used in for the (7.1). To be noticed, that in general, due to its very short period ($T_{1/2} \sim 19s$), the $^{135}$Te is absent of the evolution ; it is then taken into account in the iodine fission yield $y_I$ as explained in the Nuclei Tree section (c.f. 6).
Chapter 8

Looking at results from the evolution

8.1 The MURE data files

At present, results from the evolution are stored in the *MC* Run directory. There are 4 types of *MURE* files generated
in addition to the *MC* files for each evolution step:

*   **DATA\_xxx** files are written for each step number *xxx*, and contain not only the Material composition information, but also average cross-sections for every possible reactions for all evolving cells in the geometry. In addition, average flux for each cell, cell spatial variables (if defined), *k*<sub>eff</sub>, and time are also written out. They are written just before RK steps ONLY IF *EvolutionSolver::SetWriteASCIIData()* has been called.

*   **BDATA\_xxx** files contain the same information as the above files but in binary format. Initially it was expected that these files would be much shorter in length than the ASCII files, but it appears that not much space is saved if binary files are used to store *MURE* results. Nevertheless it is much faster to use these files for post-treatment. They are written just before RK steps by default. On can disable this writing by using *EvolutionSolver::SetWriteBinaryData(false)*.

*   **KDATA** contains all the *k*<sub>eff</sub> information (with statistical errors) for all evolution steps in one file. These data are also contained in [B]**DATA\_xxx**.

*   **FDATA** contains all the data pertaining to the fission rates of various nuclei. For all the 50 or so nuclei which can fission, the rate (in fissions per second) integrated over all the evolving cells, is written out.

When using reprocessing feature **DATA\_bxxx** or/and **BDATA\_bxxx** files are also written just at the end of
RK steps. When using *PoisonEvolutionControl*, additional files are also written such as **POISON\_PROPS** which
contains the poison proportion as a function of time.

8.2 Reading the data files with Tcl/Tk graphical interface scripts

This graphical interface is now deprecated. Except for ExamTree.tcl

*   *ExamTree.tcl*
(one for each tab).

Whereas MureGui is graphical interface, the “Terminal” window used to run this program is often used for some important warnings/information/output: don’t forget to look at it!

### 8.3.1 The 8 tabs of the main window

* In the Inventory Tab, colors of nuclei (see Fig. 8.1) correspond to the total mass scale importance (red: large mass->black: small mass). Nuclei are sorted in 4 categories (Actinides, Fission Products, Gas, and Miscellaneous) but a nucleus could belong to more than one category (a gaseous fission product is in FP and in Gas). “FP” and “AM” check boxes allows respectively the plotting of total fission products and total minor actinides (Np and Am, Cm, ...) . The total mass of checked nuclei can also be automatically performed using the “Sum of Selected” check boxes<fnref n="3" />. Finally, nuclei inventories are given in “g” or “Mol”. When several physical cells i are selected (see later with Spatial Variable combo box), the total inventory of one nucleus is simply given by the sum of inventories in each region. Thus we have:

$$N_{tot} = \sum_{i} N_{i}$$

* In the Cross-section Tab, colors only help users to see the isotopes (all reactions for a same isotope are on the same color). The “FP” check box allows the plotting of the mean fission products total cross-section. Similarly with the Inventory tab, the total cross-section can be obtained with the “Sum of Selected” check box (but FP and AM check boxes are not taken into account in the sum). Microscopic and macroscopic cross-sections are available. When several physical cells i are selected (see later with Spatial Variable combo box), an equivalent cross-section $\bar{\Sigma}$ is computed for all reactions to maintain their total reaction rate:

$$\bar{\Sigma} = \frac{\sum_{i} n_{i} \sigma_{i} \phi_{i} V_{i}}{\bar{\phi} V_{tot}}$$

where the nuclei density is $n_i = N_i/V_i$, the mean flux is $\bar{\phi} = \frac{\sum_i \phi_i V_i}{\sum_i V_i}$ and the total volume of all cells is $V_{tot} = \sum_i V_i$. The numerator is the total number of reactions in all i-cells. Thus, one can define an equivalent microscopic cross-section:

$$\bar{\sigma} = \frac{\bar{\Sigma}}{\bar{n}} = \frac{\bar{\Sigma}}{\frac{\sum_{i} n_{i} V_{i}}{\sum_{i} V_{i}}} = \frac{\sum_{i} n_{i} \sigma_{i} \phi_{i} V_{i}}{N_{tot} \bar{\phi}}.$$

**Note:**

- when the density of one nuclei is null, a very small and uniform density is used to compute the equivalent microscopic cross-section ($n = 10^{-20} \text{at.cm}^{-3}$).

- This mean cross-section is useful to compute “simplified” Bateman equations. But one has to be very attentive to its physical signification when lot of cells with different fluxes and concentrations are involved;

<footnote n="3">³But if FP and/or AM check boxes are selected, they are not part of the sum</footnote>
<figure><img src="image_31.png" /><figcaption>Figure 8.1: MureGui main window. Inventory Tab.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 8.2: Neutron Balance window. In this example, the reaction rates less than 3 per mil are grouped in the “Others” category.</figcaption></figure>

then this mean value may be greater than all individual cross-section of each cell due to an under estimation of the denominator.

- Sometimes, it may be useful to compute the mean microscopic cross-section over the spectrum:

$$ \langle \sigma \rangle = \frac{\int \sigma(E) \phi(E) dE}{\int \phi(E) dE} $$

Then one can compute a average cross-section over the cells as $\langle \sigma \rangle = \frac{\sum_i \langle \sigma_i \rangle \phi_i V_i}{\sum_i \phi_i V_i}$. In MureGui, one can plot this king of microscopic cross-section by selecting the “Micro (spectrum)” radio widget.

* In the Flux Tab, the “All cell” and “Sum of Selected” check boxes allow users to plot a mean flux over all or selected evolving cells i. The mean flux is given by:

$$ \bar{\phi} = \frac{\sum_i \phi_i V_i}{\sum_i V_i} $$

* In the Reaction Rate Tab, the “FP” and “Sum of Selected” check boxes allow users to plot total reaction rate of fission products and the sum of selected reaction rates respectively. The generated power can also be obtained using the “Power” check box. Power results are returned in the terminal when the graph is refreshed. The “(Un)Select Fissile” button selects (or unselects) all available fission reactions automatically. Note that released energy per fission values (Q_values) are automatically selected according to the transport code previously used (DRAGON, MCNP or Serpent). In this Tab, it is also possible to plot Neutron Balance normalized per 1000 fissions. Once the user has selected the wanted reaction rates (or all reaction rates for the fissile (“(Un)Select Fissile” button) or for all nuclei (“(Un)Select All” button), the neutron balance is plotted using histogram bars. The box “Threshold to Others” allow user to group all reaction rates less than the threshold (default 0.1 per mil) in a category named “Others” (see fig. 8.2 and fig 8.3). The “Dir” number is the number of the Dirname argument of MureGui and the number following the “@” character corresponds to the time at which the MC run is done (thus figure 8.2 corresponds to first Dirname at time t=1004.44 days where a MCNP run was performed)

* In the Keff Tab, the “Keff” check box allows user to plot the multiplication factor $k_{eff}$ of the assembly with evolution time. The “Average Keff” check box computes and plots the average neutron multiplication factor $\overline{k_{eff}}$ of an assembly along its life. It is particularly useful in CANDU reactors at the refueling equilibrium state. The average value $\overline{k_{eff}}$ is given by:

$$ \overline{k_{eff}}(t) = \frac{1}{t} \int_{0}^{t} k_{eff}(t') dt' $$
<figure><img src="image_33.png" /><figcaption>Figure 8.3: Neutron balance plot.</figcaption></figure>
* In the Breeding Ratio Tab, users have to select some nuclei among all heavy elements in order to define them as fissile for the calculation of Breeding Ratio (BR) and Fissile Inventory Ratio (FIR). The (instantaneous) BR is computed at a given time t as the ratio of the total production rate of the so-defined fissile nuclei at time t over the total disappearance rate of the same nuclei at the same time. Up to now, only fissile nuclei among $^{233}$U, $^{235}$U, $^{239}$Pu and $^{241}$Pu can be defined as fissile and taken into account for BR calculation, with a special treatment for $^{233}$U (which production rate is $^{233}$Pa decay rate) and $^{241}$Pu (which disappearance rate includes its decay rate).

The FIR is computed as the ratio of fissile nuclei between actual and initial time :

$$FIR = \frac{M_{\text{FN}}(t)}{M_{\text{FN}}(0)}$$

Note that if the initial fissile nuclei mass is 0, no plot is returned. In the case of $^{233}$U selected as fissile, selecting $^{233}$Pa as fissile too gives access to the effective FIR values, with each calculation time as a possible End Of Cycle. Indeed this accounts for $^{233}$Pa out-of-flux decay, producing $^{233}$U without loss. In this case, a warning is printed to inform the user that $^{233}$Pa is used this way only for FIR (BR is still computed without Pa-233 defined as fissile, and using its decay rate as $^{233}$U production rate).

* In the Radiotoxicity Tab (see also section 8.3.6), shown in Fig. 8.4, (only if Lapack library has been used), decay calculations can be performed for the final inventory after an evolution has been done. Determination of total and partial radiotoxicities, heat released and activities as a function of time can be computed for a given system. Note that if, for radiotoxicity post-treatment, MureGui is used from in a directory where the MURE "data" directory is not present just above the current path, you must define the DATADIR variable.

- Users should click on the red "Build Mat" button before taking any further action. This will compute the Material and the matrices. The material is based on the material defined in the DATA or BDATA files where all the nuclei are copied, and all their decay children. No cut-off is used to select the children. By default, spontaneous fissions are not taken into account; in order to allow it, check the "with SF" button BEFORE clicking the "Built Mat" button. The matrices for decay calculations A are also computed once and for all. The variation of the isotopic vector N is given by:

$$\frac{dN(t_c)}{dt_c} = -A.N(t_c)$$

where $t_c$ represents the cooling time (after irradiation). The solution of this equation is simply given by:

$$N(t_c) = N(0).exp(-A.t_c)$$

This equation represents a set of coupled differential equations for the decays which is solved by a matrix diagonalization procedure based on the LAPACK standard libraries (resolution is performed at the beginning). Inventory composition over time is then directly obtained from the initial composition.

- The "Nuclei Extraction" frame allows user to specify what percentage of nuclei of the current composition is taken into account for activity, radiotoxicity, ... calculation. When one choose "Partial fuel extraction" a new window appears to choose the per cent of category you keep: 100% means that all that category is kept. For example, putting 100% for all categories except 1% for U and Pu mean that the radiotoxicity is
<figure><img src="image_34.png" /><figcaption>Figure 8.4: MureGui : the radiotoxicity Tab.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Figure 8.5: “Irradiation” frame. Here, the first Dirname at time t=0, i.e., the first MCNP step, is used. The fuel cooling will be plotted from 1 year to 100 million of years with 50 log steps.</figcaption></figure>

calculated for the initial composition where only 1% of the U and Pu nuclei are kept, all other nuclei being completely taken into account. User can choose a cooling time before the extraction occurs. In that case, the plot window will start for t = 0 after the core discharge to the cooling time value^4; then, the partial extraction is done and the result is plotted. A dash red line is drawn to show the time of extraction. A “Chars” logo button allows user to plot energy spectra of γ, β, α and neutron emission of “spent fuel” (see section 8.3.6).

- In the “Miscellaneous” frame, one can also select the value to be plotted: “Activity”, “Radiotoxicity”, “Heat” released and “Inventory”. By selecting the “By mother” check button, the radiotoxicity is computed by mother nucleus (i.e. present at the starting time) in order to see the contribution to the activity of a given nuclei (the daughters of this nucleus are part of its activity). Note, the “By mother” calculation takes much more CPU time! ”FP”, “MA” and “Total” check boxes allow one to plot the total value for fission products, minor actinides, and all nuclei respectively. A combo box can be used to plot “selected” nuclei, “selected nuclei & their sum” or “only the sum” of selected nuclei (FP and MA check boxes are not taken into account in the sum).

- In the “Irradiation” frame, see Fig 8.5, one can choose the burn-up of the spent fuel (i.e., the *MCNP* step given at a time). The “Dir # num” allows to choose which evolution directory to consider in the case or more than one evolution directory is given to *MurGui* arguments. The evolution (cooling) of the spent fuel is specified by the evolution parameter: the cooling starts at *t* = 0, but the plot start at the first value of the “Evolution Period” frame ; the final time, number of time-steps and whether the steps are logarithmic or linear are then given. Note that the graph are on a log scale, thus, initial time has to be strictly positive (minimum value is 1 minute). The number of steps does not affect the results, but only the smoothness of the curve.

* The “Spatial Var” tab allows the plotting of spatial flux (or power deposit) dependence if you have defined spatial variables (see the *Cell* class and next paragraph). **It is present only if Spatial variables have been defined.** The “Plotting Time” frame allows user to select the time at which the plot is done. The “Variable choice” frame allows user to select what to plot (in the “Plotting Variable” combo box) and with, eventually, a condition on spatial variables (“Conditional Variable” combo box). For example, if a core has been divided in N radial zone and M axial zone, you can plot the mean flux as a function of Z for every radial zone (see figure 8.6). One can also plot this flux but only for a given radial zone by selecting the desired “ring” in the “Conditional Variable” combo box. The “Same Plot” check button, allow users to overlap plots.

<footnote>^4The input cooling time is then chosen as the lower bound of the time binning defined in the “Evolution period” frame ; see the exact value in the terminal window.</footnote>
<figure><img src="image_36.png" /><figcaption>Figure 8.6: Spatial Var Tab - This will plot the flux as a function of Z for t = 0 whatever the spatial variables (no condition).</figcaption></figure>

### 8.3.1.1 Out flux radiotoxicity evolution

It is possible to evaluate the time evolution of radiotoxicity (or activity, ...) of an user input composition without making any MURE/MCNP evolution. This is of course an out flux evolution (only decays are taken into account). User should create a file name `DATA_000` in a given directory (let say `TestDir`). This file must contain

```
V -1
UNIT
n
Z1 A1 I1 M1 P1
Z2 A2 I2 M2 P2
...
Zn An In Mn Pn
```

where `UNIT` is a string among `At`, `Mol`, `Wgt(g)`, `Wgt(kg)`, `Wgt(t)` giving the unit of each `proportion Pi`; `n` is the number of different isotopes (Zi,Ai,Ii) and `Mi` is the molar mass of the isotope. For example, to see the evolution of 2.975 kg of $^{238}$U the `DATA_000` file is

```
V -1
Wgt(kg)
1
92 238 0 238.05 2.975
```

Then launch MureGui as :

```
MureGui -type a TestDir
```

and go in the `Radiotoxicity` tab ; this gives a total radiotoxicity of about 1mCi for $^{238}$U at t = 0 ; then it reaches the secular equilibrium after 2 million years to $14 \times 1mCi$ (13 radioactive daughters of $^{238}U+^{238}U$).

### 8.3.2 Frame "Cell Spatial Var" in the Main window

Users can define Spatial variables for cells (see the *Cell* class). The aim of these variables is only to extract results in a simpler way (this means that you can add spatial variables after the evolution is done and redo the evolution keeping
<figure><img src="image_37.png" /><figcaption>Figure 8.9: MureGui - Interpolation and Error Bars.</figcaption></figure>

MCNP files without any problems). If such variables have been defined you can plot quantities according to spatial variables.

Suppose we have done the evolution of 3 rings of a cylindrical core divided in 3 parts along the cylinder axis. If, for each of these 9 evolving cells, 2 spatial variables named “R” and “Z” with value corresponding to the typical cell locations have been defined, then the combo box of the “Cell Spatial Var” will allow users to plot for example $^{235}$U quantity for all cells with a particular radius $R_i$ (i.e. whatever other Z value, see Fig 8.7) or to plot this quantity only for a given $R_i$ and $Z_i$ (see Fig 8.8)

<figure><img src="image_38.png" /><figcaption>Figure 8.7: MureGui - Spatial Variable combo box selected for all cells of $R = 1$.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Figure 8.8: MureGui - Spatial Variable combo box selected for the cell located at $R = 1$ and $Z = 0$.</figcaption></figure>

### 8.3.3 Frame “Time Interpolation” in the main window

A linear interpolation is performed on the plotted data for the selected Tab at the corresponding time (and not the selected data of this Tab). Results are returned in the terminal. Note that no time unit is specified, thus it must be coherent with the plotted data.

### 8.3.4 Frame “Error Bars” in the main window

When statistical errors are available, error bars may be added to the corresponding graphs. Note that MCNP provides errors on fluxes, reaction rates and $k_{eff}$ values. Thus, the Error Bars Button has no effect on the Inventory and Conversion Rate tabs. In [B]DATA_xxx files, are written $k_{eff}$, fluxes (with their MCNP absolute error), compositions of each nuclei (in atoms), and cross-sections with their error. This error is not really correct because it has been calculated as $\frac{\Delta\sigma}{\sigma} = \sqrt{\left(\frac{\Delta(\sigma\phi)}{(\sigma\phi)}\right)^2 + \left(\frac{\Delta\phi}{\phi}\right)^2}$. This is correct only if $\phi$ and $(\sigma\phi)$ are not correlated...which is obviously not the case. Thus only error for $k_{eff}$ and fluxes are really correct.

For the average flux $\bar{\phi}$ on several individual cells (with an individual flux of $\phi_i$), we have:

$$ \Delta\bar{\phi} = \frac{1}{V_{tot}} \sqrt{\sum_i (V_i \Delta\phi_i)^2} $$

If the sum of cross-sections is computed, the associated errors is given by:
$$ \Delta \overline{\sigma}_{sum} = \frac{1}{N_{tot}} \sqrt{\sum_{j} (N_{j,tot} \Delta \overline{\sigma}_{j})^2} $$

For the average reaction rate $N_{j,tot} \overline{\sigma} \overline{\phi}$ of one specific reaction j (linked to the corresponding nucleus) on several individual cells i, we have:

$$ \Delta(N_{j,tot} \overline{\sigma}_j \overline{\phi}) = \sqrt{\sum_i \Delta(N_{j,i} \sigma_{j,i} \phi_i)^2} = \sqrt{\sum_i (N_{j,i} \sigma_{j,i} \phi_i)^2 \left[ \left( \frac{\Delta \sigma_{j,i}}{\sigma_{j,i}} \right)^2 + \left( \frac{\Delta \phi_i}{\phi_i} \right)^2 \right]} $$

For the error associated with the power computation, we can used similar formula as for reaction rates, i.e.:

$$ \Delta(\sum_j P_j) = \sqrt{\sum_j \sum_i (Q_j N_{j,i} \sigma_{j,i} \phi_i)^2 \left[ \left( \frac{\Delta \sigma_{j,i}}{\sigma_{j,i}} \right)^2 + \left( \frac{\Delta \phi_i}{\phi_i} \right)^2 \right]} $$

### 8.3.5 The Plot/Save/Quit buttons in the main window

* **“Plot (All)” button:** Except for the “Cell Spatial Var” frame in the main window where by selecting a variable, the plot is automatically done, you must push the “Plot(All)” button to plot a graph; The plot window is selected according to the Tab in the main window.

* **“Save Data” button:** This button is used to save plotted graph data from the selected plot window in an ASCII file.

* **“Exec Macro” button:** This button is used to apply a ROOT/C++ interpreted (no compilation required) function to selected plot window. See for instance the **MyMacro.cxx** example in the **MURE/gui** directory. To use this macro (which does a linear fit of the ith graph), just click on the “Exec Macro” button and type **MyMacro.cxx(i)** where i is the number of the graph desired.

### 8.3.6 More about Radiotoxicity Tab

#### 8.3.6.1 Radiotoxicity of a user input composition (cooling only)

It is possible to obtain activity, inventory, heat of a user input mixture composed of n nuclei as a function of time in a “cooling” phase (i.e. outside of a flux) in a very easy way. User has to create a directory (let’s call it “tox”) containing a file named **DATA_000**. This file has the following form:

```
V -1
units
n
Z0 A0 IO MO XO
Z0 A0 IO MO XO
...
Zn An In Mn Xn
```

where *units* is one of “At”, “Mol”, “Wgt(g)”, “Wgt(kg)” or “Wgt(t)” (At=atomes, Mol=moles, Wgt(g)= mass in g, Wgt(kg)= mass in kg and Wgt(t)=mass in metric tons) that corresponds to the unit of the $X_i$ values input for each nuclei defined by its $Z_i$ (proton number), $A_i$ (atomic number), $I_i$ (isomeric state, 0 for ground state, 1 for 1st isomeric state, ...). $M_i$ is the molar mass in g/mol of the $i^{th}$ nuclei. For example, one can see how 2.975 kg of $^{238}$U reaches the secular equilibrium with