Line data    Source code

       1             : #include "crpropa/ParticleMass.h"
       2             : #include "crpropa/ParticleID.h"
       3             : #include "crpropa/Common.h"
       4             : #include "crpropa/Units.h"
       5             : 
       6             : #include "kiss/convert.h"
       7             : #include "kiss/logger.h"
       8             : 
       9             : #include <vector>
      10             : #include <fstream>
      11             : #include <stdexcept>
      12             : #include <limits>
      13             : 
      14             : namespace crpropa {
      15             : 
      16             : struct NuclearMassTable {
      17             :         bool initialized;
      18             :         std::vector<double> table;
      19             : 
      20             :         NuclearMassTable() {
      21             :                 initialized = false;
      22             :         }
      23             : 
      24          12 :         void init() {
      25          24 :                 std::string filename = getDataPath("nuclear_mass.txt");
      26          12 :                 std::ifstream infile(filename.c_str());
      27             : 
      28          12 :                 if (!infile.good())
      29           0 :                         throw std::runtime_error("crpropa: could not open file " + filename);
      30             : 
      31             :                 int Z, N;
      32             :                 double mass;
      33       10092 :                 while (infile.good()) {
      34       10080 :                         if (infile.peek() != '#') {
      35       10056 :                                 infile >> Z >> N >> mass;
      36       10056 :                                 table.push_back(mass);
      37             :                         }
      38       10080 :                         infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
      39             :                 }
      40             : 
      41          12 :                 infile.close();
      42          12 :                 initialized = true;
      43          24 :         }
      44             : 
      45      211320 :         double getMass(std::size_t idx) {
      46      211320 :                 if (!initialized) {
      47          12 : #pragma omp critical(init)
      48          12 :                         init();
      49             :                 }
      50      211320 :                 return table[idx];
      51             :         }
      52             : };
      53             : 
      54             : static NuclearMassTable nuclearMassTable;
      55             : 
      56      211296 : double nuclearMass(int id) {
      57      211296 :         int A = massNumber(id);
      58      211296 :         int Z = chargeNumber(id);
      59      211296 :         return nuclearMass(A, Z);
      60             : }
      61             : 
      62      211321 : double nuclearMass(int A, int Z) {
      63      211321 :         if ((A < 1) or (A > 56) or (Z < 0) or (Z > 26) or (Z > A)) {
      64           2 :                 KISS_LOG_WARNING <<
      65           1 :                 "nuclearMass: nuclear mass not found in the mass table for " <<
      66           1 :                 "A = " << A << ", Z = " << Z << ". " <<
      67           1 :                 "Approximated value used A * amu - Z * m_e instead.";
      68           1 :                 return A * amu - Z * mass_electron;
      69             :         }
      70      211320 :         int N = A - Z;
      71      211320 :         return nuclearMassTable.getMass(Z * 31 + N);
      72             : }
      73             : 
      74             : } // namespace crpropa