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           13 :         void init() {
      25           13 :                 if(initialized)
      26            0 :                         return;
      27           26 :                 std::string filename = getDataPath("nuclear_mass.txt");
      28           13 :                 std::ifstream infile(filename.c_str());
      29              : 
      30           13 :                 if (!infile.good())
      31            0 :                         throw std::runtime_error("crpropa: could not open file " + filename);
      32              : 
      33           13 :                 table.assign((NUCLEAR_ZMAX + 1) * NUCLEAR_NSTRIDE, 0.0);
      34              : 
      35              :                 int Z, N;
      36              :                 double mass;
      37       143572 :                 while (infile.good()) {
      38       143559 :                         if (infile.peek() != '#') {
      39       143520 :                                 infile >> Z >> N >> mass;
      40       143520 :                                 if (Z <= NUCLEAR_ZMAX && N <= NUCLEAR_NMAX)
      41       143520 :                                         table[Z * NUCLEAR_NSTRIDE + N] = mass;
      42              :                         }
      43       143559 :                         infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
      44              :                 }
      45              : 
      46           13 :                 infile.close();
      47           13 :                 initialized = true;
      48           13 :         }
      49              : 
      50       180895 :         double getMass(std::size_t idx) {
      51       180895 :                 if (!initialized) {
      52           13 : #pragma omp critical(init)
      53           13 :                         init();
      54              :                 }
      55       180895 :                 if (table[idx] == 0.0)
      56            0 :                         return 0.0; // triggers approximation in nuclearMass()
      57              :                 return table[idx];
      58              :         }
      59              : };
      60              : 
      61              : static NuclearMassTable nuclearMassTable;
      62              : 
      63     20386883 : double particleMass(int id) {
      64     20386883 :         if (isNucleus(id))
      65       178553 :                 return nuclearMass(id);
      66     20208330 :         if (abs(id) == 11)
      67        15173 :                 return mass_electron;
      68              :         return 0.0;
      69              : }
      70              : 
      71       180879 : double nuclearMass(int id) {
      72       180879 :         int A = massNumber(id);
      73       180879 :         int Z = chargeNumber(id);
      74       180879 :         return nuclearMass(A, Z);
      75              : }
      76              : 
      77       180906 : double nuclearMass(int A, int Z) {
      78       180906 :         int N = A - Z;
      79       180906 :         if ((A < 1) or (Z < 0) or (Z > A) or (Z > NUCLEAR_ZMAX) or (N > NUCLEAR_NMAX)) {
      80           22 :                 KISS_LOG_WARNING <<
      81           11 :                 "nuclearMass: nuclear mass not found in the mass table for " <<
      82           11 :                 "A = " << A << ", Z = " << Z << ". " <<
      83           11 :                 "Approximated value used A * amu - Z * m_e instead.";
      84           11 :                 return A * amu - Z * mass_electron;
      85              :         }
      86       180895 :         double m = nuclearMassTable.getMass(Z * NUCLEAR_NSTRIDE + N);
      87       180895 :         if (m == 0.0) {
      88            0 :                 KISS_LOG_WARNING <<
      89            0 :                 "nuclearMass: nuclear mass not found in the mass table for " <<
      90            0 :                 "A = " << A << ", Z = " << Z << ". " <<
      91            0 :                 "Approximated value used A * amu - Z * m_e instead.";
      92            0 :                 return A * amu - Z * mass_electron;
      93              :         }
      94              :         return m;
      95              : }
      96              : 
      97              : } // namespace crpropa