tkensdf_gamma_rec.cpp

/********************************************************************************
 *   Copyright (c) : Université de Lyon 1, CNRS/IN2P3, UMR5822,                 *
 *                   IP2I, F-69622 Villeurbanne Cedex, France                   *
 *                   Normandie Université, ENSICAEN, UNICAEN, CNRS/IN2P3,       *
 *                   LPC Caen, F-14000 Caen, France                             *
 *   Contibutor(s) :                                                            *
 *      Jérémie Dudouet jeremie.dudouet@cnrs.fr [2020]                          *
 *      Diego Gruyer    diego.gruyer@cnrs.fr    [2020]                          *
 *                                                                              *
 *    This software is governed by the CeCILL-B license under French law and    *
 *    abiding by the  rules of distribution of free  software.  You can use,    *
 *    modify  and/ or  redistribute  the  software under  the  terms of  the    *
 *    CeCILL-B license as circulated by CEA, CNRS and INRIA at the following    *
 *    URL \"http://www.cecill.info\".                                           *
 *                                                                              *
 *    As a counterpart to the access  to the source code and rights to copy,    *
 *    modify  and redistribute granted  by the  license, users  are provided    *
 *    only with a limited warranty  and the software's author, the holder of    *
 *    the economic  rights, and the  successive licensors have  only limited    *
 *    liability.                                                                *
 *                                                                              *
 *    In this respect, the user's attention is drawn to the risks associated    *
 *    with loading,  using, modifying  and/or developing or  reproducing the    *
 *    software by the user in light of its specific status of free software,    *
 *    that  may mean that  it is  complicated to  manipulate, and  that also    *
 *    therefore  means that it  is reserved  for developers  and experienced    *
 *    professionals having in-depth  computer knowledge. Users are therefore    *
 *    encouraged  to load  and test  the software's  suitability  as regards    *
 *    their  requirements  in  conditions  enabling the  security  of  their    *
 *    systems  and/or data to  be ensured  and, more  generally, to  use and    *
 *    operate it in the same conditions as regards security.                    *
 *                                                                              *
 *    The fact that  you are presently reading this means  that you have had    *
 *    knowledge of the CeCILL-B license and that you accept its terms.          *
 ********************************************************************************/
 
#include "tkensdf_gamma_rec.h"
 
#include <cmath>
 
#include "tklog.h"
 
namespace tkn {
}
 
using namespace  tkn;
using namespace  std;
 
tkensdf_gamma_rec::tkensdf_gamma_rec()
{
    ftype = kgamma;
}
 
tkensdf_gamma_rec::tkensdf_gamma_rec(const tkensdf_gamma_rec &_gamma_record) : tkensdf_record(_gamma_record)
{
    fIs_final_level_set = _gamma_record.fIs_final_level_set;
    fFinal_level_energy = _gamma_record.fFinal_level_energy;
}
 
bool tkensdf_gamma_rec::set_record(const tkstring &_record)
{
    frecord = _record;
    return true;
}
 
void tkensdf_gamma_rec::analyse_record()
{
    // check if the level is uncertain
    fUncertain_record = frecord.substr(79,1);
 
    // decode gamma energy
    decode_energy(get_energy_str(frecord));
 
    // decode relative intensity
    set_relative_intensity(frecord);
 
    // decode multipolarity
    set_multipolarity(frecord);
 
    // decode mixing ratio
    set_mixing_ratio(frecord);
 
    // decode conversion coeff
    set_conv_coeff(frecord);
 
    analyse_continuation_record();
}
 
void tkensdf_gamma_rec::set_relative_intensity(const tkstring &_record)
{
    auto RI_str  = _record.substr(21,8);
    auto DRI_str = _record.substr(29,2);
 
    if(RI_str.contains("(")) {
        frelative_intensity.info += ";VAL=()";
        RI_str.remove_all("(").remove_all(")");
    }
    frelative_intensity.value = RI_str.atof();
    if(frelative_intensity.value>0) frelative_intensity.filled = true;
    else return;
 
    frelative_intensity.unit = "%";
    frelative_intensity.err = DRI_str.atof();
}
 
void tkensdf_gamma_rec::set_multipolarity(const tkstring &_record)
{
    auto str  = _record.substr(31,10);
    if(!str.is_empty()) fmultipolarity = str.remove_all_extra_white_space();
}
 
void tkensdf_gamma_rec::set_mixing_ratio(const tkstring &_record)
{
    auto MR_str  = _record.substr(41,8).strip_all_extra_white_space();
    auto DMR_str = _record.substr(49,6).strip_all_extra_white_space();
 
    fmixing_ratio.value = MR_str.atof();
    if(fmixing_ratio.value==0) return;
 
    if(!DMR_str.is_empty()) {
 
        if(DMR_str.is_alpha()) {
            fmixing_ratio.info += ";ERR=" + DMR_str;
        }
        else {
            // define precision according to the value string
            double precision = tkstring::get_precision(MR_str);
 
            // check if asym errors
            if(DMR_str.contains("+") && DMR_str.contains("-")) {
                bool inv = DMR_str.index("-") <  DMR_str.index("+");
                DMR_str.replace_all("+"," ");
                DMR_str.replace_all("-"," ");
                auto tokens_lt_err = DMR_str.tokenize(" ");
                if(tokens_lt_err.size()!=2) {
                    glog << error << " bad asymetric error decoding in record " << frecord << do_endl;
                    return;
                }
                tkstring err_low_str = tokens_lt_err.back();
                tkstring err_high_str = tokens_lt_err.front();
                if(inv) {
                    err_low_str = tokens_lt_err.front();
                    err_high_str  = tokens_lt_err.back();
                }
 
                fmixing_ratio.err_low = err_low_str.atof()*precision;
                fmixing_ratio.err_high = err_high_str.atof()*precision;
            }
            else {
                fmixing_ratio.err = DMR_str.atof()*precision;
            }
        }
    }
    else
        fmixing_ratio.info += ";ERR=?";
 
    fmixing_ratio.filled = true;
}
 
void tkensdf_gamma_rec::set_conv_coeff(const tkstring &_record)
{
    auto CC_str  = _record.substr(55,7);
    auto DCC_str = _record.substr(62,2);
 
    fconv_coeff.value = CC_str.atof();
    if(fconv_coeff.value>0) fconv_coeff.filled = true;
    else return;
 
    fconv_coeff.unit = "";
    fconv_coeff.err = DCC_str.atof();
}
 
void tkensdf_gamma_rec::set_red_trans_prob(const tkstring &_record)
{
    auto tmp_rec = _record;
    tmp_rec.strip_all_extra_white_space();
 
    tkstring trans_name{};
 
    bool is_elec=false;
    bool is_mag=false;
    if(tmp_rec.begins_with("BM")) {
        is_mag=true;
        trans_name += "BM";
    }
    if(tmp_rec.begins_with("BE")) {
        is_elec=true;
        trans_name += "BE";
    }
    if(!(is_elec || is_mag)) return;
 
    trans_name += tmp_rec[2];
    int mult = tmp_rec.substr(2,1).atoi();
 
    bool wu = (tmp_rec.substr(3,1) == "W");
    if(wu) trans_name += "W";
 
    tkdb_table::measure_data_struct *the_struct=nullptr;
    if(is_elec && wu) the_struct = &fBEW;
    else if(is_elec && !wu) the_struct = &fBE;
    else if(is_mag && wu) the_struct = &fBMW;
    else if(is_mag && !wu) the_struct = &fBM;
    if(!the_struct) return;
    if(the_struct->filled) return;
 
    if(!wu) {
        if(is_elec) {
            the_struct->value *= pow(100,mult);
            the_struct->err *= pow(100,mult);
            the_struct->err_low *= pow(100,mult);
            the_struct->err_high *= pow(100,mult);
            the_struct->unit = tkstring::form("e^2 x (fm)^%d",mult*2);
        }
        else {
            the_struct->value *= pow(100,mult-1);
            the_struct->err *= pow(100,mult-1);
            the_struct->err_low *= pow(100,mult-1);
            the_struct->err_high *= pow(100,mult-1);
            if(mult==1) the_struct->unit = "mu^2";
            else the_struct->unit = tkstring::form("mu^2 x (fm)^%d",2*(mult-1));
        }
    }
    else the_struct->unit = "Weisskopf";
 
    the_struct->info += ";TYPE=" + trans_name;
 
    tmp_rec.remove_all(trans_name);
    // no info on what are the parenthesis, so we remove it
    tmp_rec.remove_all("(").remove_all(")");
 
    // check if the value is known
    if(tmp_rec.substr(0,1) == "=") {
        tmp_rec.remove_all("=");
        auto toks = tmp_rec.tokenize(" ");
        tkstring val = toks.at(0);
        tkstring err = "";
        if(toks.size()>1) err = toks.at(1);
 
        the_struct->value = val.atof();
 
        if(!err.is_empty()) {
            // define precision according to the value string
            double precision = tkstring::get_precision(val);
 
            // check if asym errors
            if(err.contains("+") && err.contains("-")) {
                bool inv = err.index("-") <  err.index("+");
                err.replace_all("+"," ");
                err.replace_all("-"," ");
                auto tokens_lt_err = err.tokenize(" ");
                if(tokens_lt_err.size()!=2) {
                    glog << error << " bad asymetric error decoding in record " << frecord << do_endl;
                    return;
                }
                tkstring err_low_str = tokens_lt_err.back();
                tkstring err_high_str = tokens_lt_err.front();
                if(inv) {
                    err_low_str = tokens_lt_err.front();
                    err_high_str  = tokens_lt_err.back();
                }
 
                the_struct->err_low = err_low_str.atof()*precision;
                the_struct->err_high = err_high_str.atof()*precision;
            }
            else {
                the_struct->err = err.atof()*precision;
            }
        }
        else
            the_struct->info += ";ERR=?";
    }
    else {
        tmp_rec.replace_all(">", " GT ").replace_all("<", " LT ");
 
        auto toks = tmp_rec.tokenize(" ");
        tkstring lim = toks.at(0);
        tkstring value = toks.at(1);
 
        the_struct->value = value.atof();
        the_struct->info += ";ERR=" + lim;
    }
 
    the_struct->filled = true;
}
 
void tkensdf_gamma_rec::analyse_continuation_record()
{
    auto tokens = fcontinuation_record.tokenize("$");
    for(auto &token: tokens) {
        // check if final level energy is defined
        if(token.contains("FL=")) {
            fIs_final_level_set = true;
            token.remove_all_extra_white_space().remove_all("FL=");
            if(token.is_float()) fFinal_level_energy = token.atof();
            else fFinal_level_energy = -1;
        }
        if(token.contains("BE") || token.contains("BM")) {
            set_red_trans_prob(token);
        }
    }
 
    check_red_trans_prob();
}
 
void tkensdf_gamma_rec::check_red_trans_prob()
{
    if((fBE.filled && !fBEW.filled) || (!fBE.filled && fBEW.filled) || (fBM.filled && !fBMW.filled) || (!fBM.filled && fBMW.filled) ) {
        double mass = frecord.substr(0,3).atof();
        if(!fBE.filled && fBEW.filled) {
            tkstring transname = fBEW.info.tokenize(";").front().remove_all("TYPE=");
            int mult = transname.substr(2,1).atoi();
            double conv = weisskopf_conversion(true,mass,mult);
            fBE.value = fBEW.value*conv;
            fBE.err = fBEW.err*conv;
            fBE.err_low = fBEW.err_low*conv;
            fBE.err_high = fBEW.err_high*conv;
            fBE.unit = tkstring::form("e^2 x (fm)^%d",2*mult);
            fBE.info = fBEW.info.copy().replace_all(transname,transname.copy().remove_all("W"));
            fBE.info += ";ERR=CONV";
            fBE.filled = true;
        }
        if(fBE.filled && !fBEW.filled) {
            tkstring transname = fBE.info.tokenize(";").front().remove_all("TYPE=");
            int mult = transname.substr(2,1).atoi();
            double conv = 1./weisskopf_conversion(true,mass,mult);
            fBEW.value = fBE.value*conv;
            fBEW.err = fBE.err*conv;
            fBEW.err_low = fBE.err_low*conv;
            fBEW.err_high = fBE.err_high*conv;
            fBEW.unit = "Weisskopf";
            fBEW.info = fBE.info.copy().replace_all(transname,transname.copy().append("W"));
            fBEW.info += ";ERR=CONV";
            fBEW.filled = true;
        }
        if(!fBM.filled && fBMW.filled) {
            tkstring transname = fBMW.info.tokenize(";").front().remove_all("TYPE=");
            int mult = transname.substr(2,1).atoi();
            double conv = weisskopf_conversion(false,mass,mult);
            fBM.value = fBMW.value*conv;
            fBM.err = fBMW.err*conv;
            fBM.err_low = fBMW.err_low*conv;
            fBM.err_high = fBMW.err_high*conv;
            if(mult==1) fBM.unit = "mu^2";
            else fBM.unit = tkstring::form("mu^2 x (fm)^%d",2*(mult-1));
            fBM.info = fBMW.info.copy().replace_all(transname,transname.copy().remove_all("W"));
            fBM.info += ";ERR=CONV";
            fBM.filled = true;
        }
        if(fBM.filled && !fBMW.filled) {
            tkstring transname = fBM.info.tokenize(";").front().remove_all("TYPE=");
            int mult = transname.substr(2,1).atoi();
            double conv = 1./weisskopf_conversion(false,mass,mult);
            fBMW.value = fBM.value*conv;
            fBMW.err = fBM.err*conv;
            fBMW.err_low = fBM.err_low*conv;
            fBMW.err_high = fBM.err_high*conv;
            fBMW.unit = "Weisskopf";
            fBMW.info = fBM.info.copy().replace_all(transname,transname.copy().append("W"));
            fBMW.info += ";ERR=CONV";
            fBMW.filled = true;
        }
    }
}
 
double tkensdf_gamma_rec::weisskopf_conversion(bool electric, double mass, double mult)
{
    if(electric) return pow(1.2,2.*mult)/(4.*M_PI)*3./(mult+3.)*3./(mult+3.)*pow(mass,2.*mult/3.);
    return 10./M_PI*pow(1.2,2.*mult-2.)*3./(mult+3.)*3./(mult+3.)*pow(mass,(2.*mult-2.)/3.);
}
 
void tkensdf_gamma_rec::clear()
{
    tkensdf_record::clear();
 
    fIs_final_level_set = false;
    fFinal_level_energy = -1.;
    frelative_intensity.clear();
    fmultipolarity.clear();
    fmixing_ratio.clear();
    fconv_coeff.clear();
    fBE.clear();
    fBEW.clear();
    fBM.clear();
    fBMW.clear();
 
}
 
void tkensdf_gamma_rec::print(std::ostream &/*out*/) const
{
 
}