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       147202 :         double getMass(std::size_t idx) {
      46       147202 :                 if (!initialized) {
      47           12 : #pragma omp critical(init)
      48           12 :                         init();
      49              :                 }
      50       147202 :                 return table[idx];
      51              :         }
      52              : };
      53              : 
      54              : static NuclearMassTable nuclearMassTable;
      55              : 
      56     20338042 : double particleMass(int id) {
      57     20338042 :         if (isNucleus(id))
      58       146180 :                 return nuclearMass(id);
      59     20191862 :         if (abs(id) == 11)
      60        14340 :                 return mass_electron;
      61              :         return 0.0;
      62              : }
      63              : 
      64       147178 : double nuclearMass(int id) {
      65       147178 :         int A = massNumber(id);
      66       147178 :         int Z = chargeNumber(id);
      67       147178 :         return nuclearMass(A, Z);
      68              : }
      69              : 
      70       147203 : double nuclearMass(int A, int Z) {
      71       147203 :         if ((A < 1) or (A > 56) or (Z < 0) or (Z > 26) or (Z > A)) {
      72            2 :                 KISS_LOG_WARNING <<
      73            1 :                 "nuclearMass: nuclear mass not found in the mass table for " <<
      74            1 :                 "A = " << A << ", Z = " << Z << ". " <<
      75            1 :                 "Approximated value used A * amu - Z * m_e instead.";
      76            1 :                 return A * amu - Z * mass_electron;
      77              :         }
      78       147202 :         int N = A - Z;
      79       147202 :         return nuclearMassTable.getMass(Z * 31 + N);
      80              : }
      81              : 
      82              : } // namespace crpropa