==== Technical design description ==== 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. === What is an SDID? === 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 [[atris: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. === What is an DIID? === 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 === What is a PDG code? === 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[(http://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf)]. Note: Antiparticles are designated with a negative sign. === How is the ICRU sphere water phantom response data read-in and stored in RAM? === Please read the [[AtRIS: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. == The energy index $k$ == 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 === How does code determine if a given secondary particle is to be used for the dose-rate estimation ? === The author has performed an analysis of the secondary particles population[(Note: plots will be made after publication available)]. He has determined that only the following particles contribute to the dose rate: ^ particle: ^ GPS name: ^ PDG code: ^ Frequency: ^ Dose contribution: ^ | $\gamma$ | photon | 22 | very frequent | normal, lowest | | $\pi^+,\pi^-$ | pions | 221,-221 | infrequent | enhanced through decay | | $K^+,K^-1,K^0$ | kaons | 321,-321,130 | very rare | enhanced through decay | | $e$ | e- | 11 | frequent | normal, low | | $\bar{e}$ | e+ | -11 | frequent | enhanced through annihilation | | $\mu$ | mu- | 13 | frequent | enhanced through decay | | $\bar{\mu}$ | mu+ | -13 | frequent | enhanced through decay | | $\bar{p}$ | proton | 2212 | frequent | normal, low | | $\bar{p}$ | anti_proton | -2212 | infrequent | enhanced through annihilation | | $n,\bar{n}$ | neutron | 2112, | frequent | enhanced through decay | | $\bar{n}$ | anti_neutron | -2112 | infrequent | enhanced through decay, annihilation | | $d,t$ | deuteron,triton | 1000010020,1000010030 | very rare | normal, low | | $He^3,\alpha$ | Helium3, alpha | 1000020030,1000020040 | very rare | normal, low | 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 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 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. == Structure of log data? == AtRIS produces 3 different log files: - **.log** is the main log file. It contains lots of information. Check it out! - **.out** is a copy of the stdout - **.error** is a redirection of the stderr. Check this if the simulation crushes. Here are a couple of reasons: - debugging - debugging - reproducability - documentation - diminishing performance impact since we are not writing out to stdout, but to the hard drive - ~~REFNOTES~~ ~~DISCUSSION|Ask a question or suggest an improvement~~