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