The author of AtRIS is determined to provide enough documentation so that the project does not get abandoned because of lack of thereof (Other reasons are OK). Still, in comparison to data analysis, writing documentation can be a rather tedious and unrewarding task. This here is an effort to describe some of the design decisions that the author has made to make things work as efficiently as possible. If you see points which could be further optimized to significantly improve performance, please leave a ticket in the discussion section which can be found at the bottom of the page.

Sensitive Detector ID is the integer designation of atmospheric volume elements (sensitive detectors). The SDID is used throughout the code to as an address for data storage, for example when filling in the ionization histograms. Therefore, the SDID must be properly specified for each volume element, as has been described in input. Note, to improve performance, particles are killed when they enter the core volume (SDID=0) or the loss volume (largest SDID). Still, even though no particle tracking is performed in these two special volumes, they enter in to the histogram generation. The user needs to be aware of this behavior.

Detector Interface ID designates the interface between the volumes with SDID=DIID and SDID=DIID+1. In the code, AtRIS looks for particles which have just entered a new volume and it knows the SDID of the new volume. In order to determine the DIID, we need to know if the particle is moving in a downward or in an upward direction:

 
G4ThreeVector position  = pre->GetPosition().unit();   // the position vector normalized to 1
G4ThreeVector direction = pre->GetMomentumDirection(); // a normalized direction vector
G4double inclination    = (position*direction);        // a scalar product of two unit vectors give the enclosed angle
bool upward             = inclination < halfpi;        // flag which is true if the particle is upward
G4int diid              = sdid - (1*upward);           // to get the diid from the sdid, reduce by 1 if the particle is upward

For Geant4 it was necessary to develop an integer naming scheme which can be used to specify all particles. If numerically we want to check the type of a particle, it is much easier to consider an integer designation than parse a string. The Particle Data Group was tasked with creating such a numbering scheme. The PDG naming scheme covers all possible particles, ranging from quarks, diquarks, etc., to all the known nuclei. More information can be found here¹⁾.

Note: Antiparticles are designated with a negative sign.

Please read the model_description section on the ICRU sphere water phantom. The response step-functions for the different relevant particles (see below) are stored as a an ascii table with a single line header. The header, starting with “#” specifies the $$\rm{PDG}\quad\Leftrightarrow\quad j,$$ correspondence, where $j$ is the column number. To parse this data, the procedure is as follows:

• open the file icru_response.txt.
• parse the header and extract all the PDG codes, count how many you have (nrpdg).
• store them in given order in to an G4int array of length nrpdg named order, while keeping the order information ($j$).
• construct an array named icru of doubles with length 650 (nr of energy bins) x nrpdg (number of relevant particles). While this is a sequential array, the user can imagine it to be a 2d array with 650 rows and nrpdg columns. We access the individual entries with a single index = 650*j + k instead of with a pair of indices. Here $k$ designates the energy index. We will discuss latter how to determine in simulation what is the correct $k$.
• The number of relevant particles is not hardcoded and can be modified by replacing the icru_response.txt file.

For the simulation of $j$, we have used 650 energy bins ranging from $1$ eV to $10$ TeV, equidistant on a $\log_{10}$ scale:

import numpy as np                  # not MeV is the default energy unit. Thus 1e0=1MeV
exponents = np.linspace(-6,7,651)   # generate 651 binedges with equidistant 0.02 spacing
binedges  = 10**exponents           # raise 10 to the calculated exponents
np.savetxt("new_icru.bins",binedges)# save as an txt file 

This energy binning is hard-coded within the code of AtRIS. Namely, the $k$ index selection is hardcoded. The number of energy bins is hardcoded in form of a global variable nricrubins=650. Here is the code which:

• finds the kinetic energy at boundary $E$
• checks in to which bin the particle fits (determines $k$)

G4double    step = 0.02;                          // untransformed step size
G4StepPoint* pre = step->GetPreStepPoint();       // get step point at the beginning of the step
G4double    ekin = pre->GetKineticEnergy();       // retrieve kinetic energy
G4int          k = (G4int) (log10(ekin/eV)/0.02); // we must divide by eV in order to get positive indices

The author has performed an analysis of the secondary particles population²⁾. He has determined that only the following particles contribute to the dose rate:

Lets call the above particles relevant.

The dose rate is estimated each time a secondary passes a volume interface. The procedure goes as follows:

• Was the particle at the beginning of the step (prestep-point) at an volume boundary? If not, stop.
• Get the PDG code.
• At this point we need to check if the PDG code corresponds to one of the above particles. Here are some options:
  • trick: If the particle is in, so is the anti-particle!
  • Create a dictionary of relevant particles and check if the current particle is one of those. (dictionary iteration is expensive!)
  • Create an array of bools (zeros) with as much elements as possible PDG codes and then change those entries whose location corresponds to a PDG of one of the relevant particles to one. A particle is to be used if array[PDG] is true. (this would consume even with bool storing using bit packing more then one GB of ram per process).
  • trick: Except for the couple of heavy ions ($d$,$t$,$He^3$,$\alpha$), all relevant particles have a PDG code that is in magnitude less than 2213.
  • If the PDG is less than 2213, Use the above strategy. Otherwise if the we have a heavy ion (PDG above 1000010001), use dictionary lookup.
  • Problem: The class std::vector<bool> is a horror to use. Typically, the smallest block of memory which can be addressed is one byte. That is, the CPU addresses memory by bytes. As a consequence, when you store a bool, even though it requires only one bit of storage, it occupies 1 bytes of memory. To avoid this, the C++ committee has implement bit packing in to the std::vector<bool> class, and It is generally a good strategy to avoid using it at all costs.
  • We thus store the relevant particles flags in an array named useornot of ints with a length of 2213*2+1. You can do if (useornot[pdg+2213]) and this will return you True, if useornot[pdg+2213] is nonzero. The offset 2213 is necessary to accommodate for negative PDG values (antiparticles). Furthermore, the actual value of useornot[pdg] is the index $j$ (see above), which specifies the column in the icru array which corresponds to the particle with the specified PDG code. Thus now, if we determine the energy index $k$, we can access the factor $f$ from the icru array.
  • The above scheme involves only addition and a single array element access and is quite fast.
• Now we Determine the DIID, as described above.

AtRIS produces 3 different log files:

1. .log is the main log file. It contains lots of information. Check it out!
2. .out is a copy of the stdout
3. .error is a redirection of the stderr. Check this if the simulation crushes.

Here are a couple of reasons:

1. debugging
2. debugging
3. reproducability
4. documentation
5. diminishing performance impact since we are not writing out to stdout, but to the hard drive