Beta Spectrum Generator程序说明
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 usingHBAR=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+
orB-
)Transition type
(Fermi
,Gamow-Teller
, ormixed
)
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
andgM
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 optionshoFit
.
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
isModified_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
andbPos
.
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
andaPos
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.
- 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.
- 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
.
- Generally speaking, only one set(line) of values should be accepted in
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 toW0
.
- But when
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 frombeginW
, ends inendW
. Every cycle,currentW
addsstepW
once. - Each entry calls
CalculateDecayRate()
once, and push the return value into the 2-dimensional arrayspectrum
mentioned above. - Raw input
currentW
giveCalculateDecayRate()
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 resultresult
by this factor. - Calculate factor (use
Wv
), then multiply the current calculation neutrino resultneutrinoResult
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
isB-
.
- It enables only when
- Atomic Mismatch: The atomic mismatch correction.
- It enables only when
atomicEnergyDeficit
equals to 0.
- It enables only when
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.