Beta Spectrum Generator程序说明

来自Jimmy's Wiki
跳到导航 跳到搜索

BSG_EXEC running process

1. Initialize instance

    bsg::BSGOptionContainer::GetInstance(argc, argv);

Get params via ParseOptions Library (PO)

All params are stored in a variables map.

Define all options in PO

  • transitionOptions
  • configOptions
    • spectrumOptions
    • constants
  • genericOptions

Parse options from terminal input

Parse config options

Parse input options

Initialize NME

2. Construct the environment

    bsg::Generator* gen = new bsg::Generator();

Initialize loggers

Output logs.

Initialize constants

Load constants.

Z, A, and R

R is the nuclear radius in natural units (HBAR=c=m_e=1).

A is mass number.

Z is the proton number of the daughter nucleus.

  • Z is defined(continue to use) by its daughter Z’. A is the same with Z.
    Z = GetBSGOpt(int, Daughter.Z);
    A = GetBSGOpt(int, Daughter.A);
  • A must equal to its mother. Z must equal to its mother + betatype.
    if (A != GetBSGOpt(int, Mother.A)) {
        consoleLogger->error("Mother and daughter mass numbers are not the same.");
    }
    if (Z != GetBSGOpt(int, Mother.Z)+betaType) {
        consoleLogger->error("Mother and daughter cannot be obtained through {} process", process);
    }
  • However, R have judgements.
    R = GetBSGOpt(double, Daughter.Radius) * 1e-15 / NATURAL_LENGTH * std::sqrt(5. / 3.);
    if (R == 0.0) {
        debugFileLogger->debug("Radius not found. Using standard formula.");
        R = 1.2 * std::pow(A, 1. / 3.) * 1e-15 / NATURAL_LENGTH;
    }

NATURAL_LENGTH is the beta decay natural length scale using HBAR=c=m_e=1. 3.861592643659598e-13 m

These constants get from previous config unconditionally

  • motherBeta2
  • daughterBeta2
  • motherSpinParity
  • daughterSpinParity
  • motherExcitationEn
  • daughterExcitationEn
  • gA
  • gP
  • gM
  • QValue
  • atomicEnergyDeficit
  • Beta type
  • Process(B+ or B-)
  • Transition type(Fermi, Gamow-Teller, or mixed)

QValue is Q value of the particular decay.

motherBeta2 is the quadrupole deofrmation of the mother nucleus.

daughterBeta2 is the quadrupole deformation of the daughter nucleus.

motherSpinParity is the double of the spin parity of the mother state.

daughterSpinParity is the double of the spin parity of the daughter state.

gA, gP and gM are coupling constants of the weak Hamiltonian.

atomicEnergyDeficit is explicit energy difference between Q value and endpoint energy due to incomplete atomic overlap (atomic excitations).

Beta type is internal state of the beta type, 1 in 2 values: BETA_PLUS(-1), BETA_MINUS(1).

Process is the decay process, 1 in 2 valid values: B+, B-.

Transition type is the type of beta decay, 1 in 3 valid values: FERMI, GAMOW_TELLER, MIXED.

W0

W0 is the total electron energy in units of its rest mass.

W0 is determined by Process.

    if (betaType == BETA_MINUS) {
        W0 = (QValue - atomicEnergyDeficit + motherExcitationEn - daughterExcitationEn) / ELECTRON_MASS_KEV + 1.;
    } else {
        W0 = (QValue - atomicEnergyDeficit + motherExcitationEn - daughterExcitationEn) / ELECTRON_MASS_KEV - 1.;
    }
    W0 = W0 - (W0 * W0 - 1) / 2. / A / (NUCLEON_MASS_KEV / ELECTRON_MASS_KEV);

ELECTRON_MASS_KEV is the electron rest mass in keV. 510.998902 keV

NEUTRON_MASS_KEV is the neutron rest mass in keV. 939565.378 keV

PROTON_MASS_KEV is the proton rest mass in keV. 938271.998 keV

NUCLEON_MASS_KEV is the average nucleus mass in keV. (NEUTRON_MASS_KEV + PROTON_MASS_KEV) / 2

Initialize shape parameters

Gauss fit

Detect Gauss Fit.

  • If Gauss Fit is disabled, fit the charge distribution as constructed using Harmonic Oscillator functions to a Modified Gaussian distribution.
  • If Gauss Fit is enabled, use previous options hoFit.
    if (!BSGOptExists(Spectrum.ModGaussFit)) {
        hoFit = CD::FitHODist(Z, R * std::sqrt(3. / 5.));
    } else {
        hoFit = GetBSGOpt(double, Spectrum.ModGaussFit);
    }

hoFit is the fit value obtained after fitting the nuclear charge distribution with a Modified Gaussian.

v and v’

v, aka vOld, is the power expansion in r for old electrostatic potential.

v', aka vNew, is the power expansion in r for new electrostatic potential.

ESShape is the name denoting the base shape to use for the U correction.

NSShape is the name denoting the base shape to use for the C correction.

Detect ESShape and NSShape. NSShape uses the params before. While ESShape can determine v and v’.

  • If ESShape is Modified_Gaussian, use preset v and v’.
  • If v and v’ are defined before, use previous options.
  • As long as one of v and v’ is not defined before and now, throw error.
    if (ESShape == "Modified_Gaussian") {
        debugFileLogger->debug("Found Modified_Gaussian shape");
        vOld[0] = 3./2.;
        vOld[1] = -1./2.;
        vNew[0] = std::sqrt(5./2.)*4.*(1.+hoFit)*std::sqrt(2.+5.*hoFit)/std::sqrt(M_PI)*std::pow(2.+3.*hoFit, 3./2.);
        vNew[1] = -4./3./(3.*hoFit+2)/std::sqrt(M_PI)*std::pow(5.*(2.+5.*hoFit)/2./(2.+3.*hoFit), 3./2.);
        vNew[2] = (2.-7.*hoFit)/5./(3.*hoFit+2)/std::sqrt(M_PI)*std::pow(5.*(2.+5.*hoFit)/2./(2.+3.*hoFit), 5./3.);
    } else {
        if (BSGOptExists(Spectrum.vold) && BSGOptExists(Spectrum.vnew)) {
            debugFileLogger->debug("Found v and v'");
            vOld = GetBSGOpt(std::vector<double>, Spectrum.vold);
            vNew = GetBSGOpt(std::vector<double>, Spectrum.vnew);
        } else if (BSGOptExists(vold) || BSGOptExists(vnew)) {
            consoleLogger->error("ERROR: Both old and new potential expansions must be given.");
        }
    }

Initialize L0 constants

Calculate aNeg and aPos

aNeg is an array containing Wilkinson’s fit coefficients for the L0 correction for beta- decay.

aPos is an array containing Wilkinson’s fit coefficients for the L0 correction for beta+ decay.

  • Define bNeg and bPos.
bNeg[0] bNeg[1] bNeg[2] bNeg[3] bNeg[4] bNeg[5] bNeg[6]
0.115 -0.00062 0.02482 -0.14038 0.008152 1.2145 -1.5632
-1.8123 0.007165 -0.5975 3.64953 -1.15664 -23.9931 33.4192
8.2498 0.01841 4.84199 -38.8143 49.9663 149.9718 -255.1333
-11.223 -0.53736 -15.3374 172.1368 -273.711 -471.2985 938.5297
-14.854 1.2691 23.9774 -346.708 657.6292 662.1909 -1641.2845
32.086 -1.5467 -12.6534 288.7873 -603.7033 -305.6804 1095.358
bPos[0] bPos[1] bPos[2] bPos[3] bPos[4] bPos[5] bPos[6]
0.0701 -0.002308 0.07936 -0.93832 4.276181 -8.2135 5.4583
-2.572 0.066463 -2.09284 22.02513 -96.82411 179.0862 -115.8922
27.5971 -0.6407 18.45462 -197.00221 835.26505 -1492.1295 940.8305
-128.658 2.63606 -80.9375 807.1878 -3355.8441 5872.5362 -3633.9181
272.264 -5.6317 160.8384 -1566.6077 6411.3255 -11038.7299 6727.6296
-214.925 4.0011 -124.8927 1156.3287 -4681.573 7963.4701 -4795.0481
    double bNeg[7][6];
    //...
    double bPos[7][6];
    //...
  • Calc aNeg and aPos using codes below.
    for (int i = 0; i < 7; i++) {
        aPos[i] = 0;
        aNeg[i] = 0;
        for (int j = 0; j < 6; j++) {
            aNeg[i] += bNeg[i][j] * std::pow(ALPHA * Z, j + 1);
            aPos[i] += bPos[i][j] * std::pow(ALPHA * Z, j + 1);
        }
    }

ALPHA is the fine-structure constant. 0.0072973525664

Load exchange parameters (optional)

If Atomic Exchange is enabled, this process will be gone through.

Open exchange data file

Select the right exchange data

exPars is an array for the fit coefficients of the atomic exchange correction.

  1. Exchange Data file contains a serial number and nine values per line, so the first step is to get each value.
    • z is related to beta type.
  2. If beta type corresponds with z, exPars will be set to the correct value in the Exchange Data file.
    • Generally speaking, only one set(line) of values should be accepted in exPars.
    if (paramStream.is_open()) {
        while (getline(paramStream, line)) {
            double z;
            double a, b, c, d, e, f, g, h, i;
            std::istringstream iss(line);
            iss >> z >> a >> b >> c >> d >> e >> f >> g >> h >> i;
            if (z == Z - betaType) {
                exPars[0] = a;
                exPars[1] = b;
                exPars[2] = c;
                exPars[3] = d;
                exPars[4] = e;
                exPars[5] = f;
                exPars[6] = g;
                exPars[7] = h;
                exPars[8] = i;
            }
        }
    } else {
        consoleLogger->error("ERROR: Can't find Exchange parameters file at {}.", exParamFile);
    }

Initialize NSM info

Initialize Nuclear Structure Manager

    nsm = new NS::NuclearStructureManager();

Nuclear Structure Manager is a array in NME that contains elements related to matrix.

More details will be discussed in NME_EXEC(module) running process.

Initialize ESP state in NME module

    if (BSGOptExists(connect)) {
        int dKi, dKf;
        nsm->GetESPStates(spsi, spsf, dKi, dKf);
    }

spsi, spsf are single particle states calculated from the NME library and used in the C_I correction when turned on.

3. Calculate spectrum

Define a 2-dimensional array for data storage

It is in line with the raw output which is a 2-dimensional array called spectrum.

    spectrum = new std::vector<std::vector<double> >();

Get start, end points and step of the calculation

Get beginEn and endEn and convert them into beginW and endW

beginEn is the starting energy in keV from which to start the spectrum calculation.

endEn is the last energy in keV for which to calculate the spectrum.

beginEn and endEn are both from config.

    double beginEn = GetBSGOpt(double, Spectrum.Begin);
    double endEn = GetBSGOpt(double, Spectrum.End);

Calc beginW and endW. Guess that W means “relative”. Same below.

  • beginW = beginEn / ELECTRON_MASS_KEV + 1.
  • endW = endEn / ELECTRON_MASS_KEV + 1.
    • But when endEn equals to 0, endW will equal to W0.
    double beginW = beginEn / ELECTRON_MASS_KEV + 1.;
    double endW = endEn / ELECTRON_MASS_KEV + 1.;
    if (endEn == 0.0) {
        endW = W0;
    }

Get stepW

Steps is the number of steps in the total spectrum. StepSize is the stepsize in keV.

If you specify the Steps in config, stepW will use (endW - beginW) / Steps first. Otherwise, stepW will equal to StepSize / ELECTRON_MASS_KEV.

    double stepW = GetBSGOpt(double, Spectrum.StepSize) / ELECTRON_MASS_KEV;
    if (BSGOptExists(Spectrum.Steps)) {
        stepW = (endW-beginW)/GetBSGOpt(int, Spectrum.Steps);
    }

Start calculation

Loop for each entry(datum)

  • The loop pointer currentW starts from beginW, ends in endW. Every cycle, currentW adds stepW once.
  • Each entry calls CalculateDecayRate() once, and push the return value into the 2-dimensional array spectrum mentioned above.
  • Raw input currentW give CalculateDecayRate() to get return value from this function.
    double currentW = beginW;
    while (currentW <= endW) {
        auto result = CalculateDecayRate(currentW);
        std::vector<double> entry = {currentW, std::get<0>(result), std::get<1>(result)};
        spectrum->push_back(entry);
        currentW += stepW;
    }

So here comes the result: all corrections and calculations are in CalculateDecayRate() and based on raw input currentW and config.

Calculate decay rate for each discrete data

This process is the logic of CalculateDecayRate().

The function accepts one parameter: W. currentW passes to this function as W.

Variables
  • result
    • Default value: 1
  • neutrinoResult
    • Default value: 1
    double result = 1;
    double neutrinoResult = 1;
Constants
  • Wv = W0 - W + 1

W0 is the total electron energy in units of its rest mass.

    double Wv = W0 - W + 1;
Corrections

There are several corrections.

There may be some judgment conditions in some of the correction processes, but in general, every correction has 2 steps:

  • Calculate factor (use W), then multiply the current calculation result result by this factor.
  • Calculate factor (use Wv), then multiply the current calculation neutrino result neutrinoResult by this factor.

Here are the corrections:

  • Phasespace: The phase space factor.
  • Fermi: The Fermi Function.
  • C: C correction.
    • This correction affects by NME connection.
  • Isovector: The isovector correction to C.
  • Relativistic: Relativistic corrections.
  • ES Deformation: Deformation corrections.
  • ES Finite Size: L0 correction.
  • U: U correction.
  • Coulomb Recoil: Q correction.
  • Radiative: Radiative corrections.
  • Recoil: Kinematic recoil corrections.
  • Screening: Atomic screening.
  • Exchange: Atomic exchange.
    • It enables only when betatype is B-.
  • Atomic Mismatch: The atomic mismatch correction.
    • It enables only when atomicEnergyDeficit equals to 0.

atomicEnergyDeficit is the explicit energy difference between Q value and endpoint energy due to incomplete atomic overlap (atomic excitations).

Return

If result and neutrinoResult are greater than 0, then return this number. Otherwise, return 0.

Prepare output file