BmnFunctionSet.cxx

Go to the documentation of this file.

#include "BmnFunctionSet.h"
 
#include "FairFileHeader.h"
#include "FairLogger.h"
#include "FairRootManager.h"
#include "TChain.h"
#include "TChainElement.h"
#include "TKey.h"
#include "TObjString.h"
#include "TPRegexp.h"
#include "TSystem.h"
#include "TTree.h"
 
#include <TMath.h>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
using namespace std;
using namespace std::chrono;
using namespace std::chrono_literals;
 
// check whether file exists: 1 - exists, 0 - not exists, -1 - exists but cannot access with the mode (default:
// kFileExists)
int BmnFunctionSet::CheckFileExist(TString& fileName, int iVerbose, EAccessMode mode)
{
    gSystem->ExpandPathName(fileName);
    if (gSystem->AccessPathName(fileName.Data(), kFileExists) == true) {
        if (iVerbose)
            cout << "ERROR: no specified file: " << fileName << endl;
        return 0;
    }
    if ((mode != kFileExists) && (gSystem->AccessPathName(fileName.Data(), mode) == true)) {
        if (iVerbose)
            cout << "ERROR: no required permissions to access the file: " << fileName << endl;
        return -1;
    }
 
    return 1;
}
 
// check whether directory exists: 1 - exists, 0 - not exists, -1 - exists but cannot access with the mode (default:
// kFileExists)
int BmnFunctionSet::CheckDirectoryExist(TString& fileName, int iVerbose, EAccessMode mode)
{
    gSystem->ExpandPathName(fileName);
    TString dirName(gSystem->DirName(fileName.Data()));
    if (gSystem->AccessPathName(dirName.Data(), kFileExists) == true) {
        if (iVerbose)
            cout << "ERROR: no specified directory: " << dirName << endl;
        return 0;
    }
    if ((mode != kFileExists) && (gSystem->AccessPathName(dirName.Data(), mode) == true)) {
        if (iVerbose)
            cout << "ERROR: no required permissions to access the directory: " << dirName << endl;
        return -1;
    }
 
    return 1;
}
 
// create directory tree for the file name if not exists: 1 - exists, 0 - created
// error codes: -1 - exists but cannot access with the mode (default: kWritePermission), -2 - cannot be created, -3 -
// empty file name
int BmnFunctionSet::CreateDirectoryTree(TString& fileName, int iVerbose, EAccessMode mode)
{
    if (fileName == "") {
        if (iVerbose)
            cout << "ERROR: the specified file name is empty" << endl;
        return -3;
    }
 
    int status = CheckDirectoryExist(fileName, 0, mode);
    if (status < 0) {
        if (iVerbose)
            cout << "ERROR: the directory of the specified file (\"" << fileName
                 << "\") does not provide required permissions" << endl;
        return status;
    }
    if (status == 0) {
        TString dirName(gSystem->DirName(fileName.Data()));
        int status_parent = CreateDirectoryTree(dirName, 0, mode);
        if (status_parent < 0)
            return status_parent;
        if (gSystem->MakeDirectory(dirName) != 0) {
            if (iVerbose)
                cout << "ERROR: the directory of the specified file (\"" << fileName << "\") could not be created"
                     << endl;
            return -2;
        }
        return 0;
    }
 
    return 1;
}
 
// get the list of input ROOT files for the given input file, which can be a ROOT file with data, a ROOT file with
// a TChain object or a text file (extension is not '.root') with a list of ROOT files
// return value: number of files in the list, < 0 - error codes
int BmnFunctionSet::GetInputFileList(TString inputFileName, vector<TString>& fileList)
{
    fileList.clear();
    if (CheckFileExist(inputFileName) != 1) {
        LOG(error) << "BmnFunctionSet::GetInputFileList: no specified file: " << inputFileName;
        return -1;
    }
 
    if (inputFileName.EndsWith(".root")) {
        TFile* rootFile = TFile::Open(inputFileName.Data());
        if (!rootFile || rootFile->IsZombie()) {
            LOG(error) << "BmnFunctionSet::GetInputFileList: opening the ROOT file failed: " << inputFileName;
            delete rootFile;
            return -2;
        }
 
        // search for a TChain object, otherwise the ROOT file itself contains the data
        TChain* chain = nullptr;
        TIter next(rootFile->GetListOfKeys());
        TKey* key;
        while ((key = (TKey*)next()) != nullptr) {
            if (strcmp(key->GetClassName(), "TChain") == 0) {
                chain = rootFile->Get<TChain>(key->GetName());
                break;
            }
        }
        if (chain == nullptr)
            fileList.push_back(inputFileName);
        else {
            TIter nextElement(chain->GetListOfFiles());
            TChainElement* element;
            while ((element = (TChainElement*)nextElement()) != nullptr)
                fileList.push_back(element->GetTitle());
        }
        rootFile->Close();   // the chain object is deleted together with the file
        delete rootFile;
    } else {
        ifstream listFile(inputFileName.Data());
        if (!listFile.is_open()) {
            LOG(error) << "BmnFunctionSet::GetInputFileList: the file list could not be opened: " << inputFileName;
            return -2;
        }
 
        string line;
        while (getline(listFile, line)) {
            istringstream lineStream(line);
            string word;
            if ((lineStream >> word) && (word[0] != '#'))
                fileList.push_back(word);
        }
    }
 
    // the check also expands the paths, so the files can be given as, e.g., "$VMCWORKDIR/bmndst.root"
    for (auto& fileName : fileList) {
        if (CheckFileExist(fileName) != 1) {
            LOG(error) << "BmnFunctionSet::GetInputFileList: no specified file: " << fileName << " (listed in "
                       << inputFileName << ")";
            fileList.clear();
            return -3;
        }
    }
    if (fileList.empty()) {
        LOG(error) << "BmnFunctionSet::GetInputFileList: no input files found in: " << inputFileName;
        return -4;
    }
 
    return fileList.size();
}
 
int BmnFunctionSet::isSimulationFile(TString fileName)
{
    TFile* fRootFile = nullptr;
    if (fileName == "") {
        FairRootManager* ioman = FairRootManager::Instance();
        if (ioman == nullptr) {
            LOG(error) << "isSimulationFile(): FairRootManager is not instantiated";
            return -1;
        }
        TTree* cur_tree = ioman->GetInTree();
        if (cur_tree == nullptr) {
            LOG(error) << "isSimulationFile(): There is no current tree opened";
            return -2;
        }
 
        fRootFile = cur_tree->GetCurrentFile();
    } else {
        gSystem->ExpandPathName(fileName);
        if (gSystem->AccessPathName(fileName.Data()) == true) {
            LOG(error) << "isSimulationFile(): no specified file: " << fileName;
            return -3;
        }
 
        fRootFile = TFile::Open(fileName.Data());
        if (!fRootFile || fRootFile->IsZombie()) {
            LOG(error) << "isSimulationFile(): opening the input file failed: " << fileName;
            if (fRootFile) {
                fRootFile->Close();
                delete fRootFile;
            }
            return -4;
        }
    }
 
    TObject* branch_list = fRootFile->Get("BranchList");
    if (branch_list != nullptr) {
        TObject* mctrack_find = branch_list->FindObject("MCTrack");
        if (mctrack_find != nullptr) {
            if (fileName != "") {
                fRootFile->Close();
                delete fRootFile;
            }
 
            return 1;
        }
    }
 
    if (fileName != "") {
        fRootFile->Close();
        delete fRootFile;
    }
    return 0;
}
 
int BmnFunctionSet::GetPeriodNumber(TFile* root_file, int method_number)
{
    if (!root_file || root_file->IsZombie()) {
        LOG(error) << "BmnFunctionSet::GetPeriodNumber: ROOT-file is not opened " << root_file->GetName();
        return -2;
    }
 
    int period_number = -1;
    switch (method_number) {
        // method number: 0
        case 0: {
            TIter next(root_file->GetListOfKeys());
            TKey* key;
            while ((key = (TKey*)next())) {
                if (std::string(key->GetClassName()) == "FairRtdbRun") {
                    TString str_period_number = key->GetName();
                    if (str_period_number.IsDigit()) {
                        period_number = str_period_number.Atoi();
                        break;
                    }   // if (str_period_number.IsDigit())
                }   // if "FairRtdbRun" class
            }   // while TKey*)next()
 
            if (period_number == -1) {
                LOG(error) << "BmnFunctionSet::GetPeriodNumber: FairRtdbRun object not found in ROOT-file "
                           << root_file->GetName();
                return -3;
            }
            return period_number;
        }   // method number == 0
        // method number: 1
        case 1: {
            FairFileHeader* header = root_file->Get<FairFileHeader>("FileHeader");
            if (!header) {
                LOG(error) << "BmnFunctionSet::GetPeriodNumber: 'FileHeader' branch not found in ROOT-file "
                           << root_file->GetName();
                return -4;
            }
 
            period_number = (int)header->GetRunId();
            return period_number;
        }
        default:
            return -5;
    }
 
    return -6;
}
 
TString BmnFunctionSet::JoinKeys(const std::vector<std::pair<TString, TString>>& dict)
{
    TString result;
    for (size_t i = 0; i < dict.size(); ++i) {
        if (i > 0)
            result += "|";
        result += dict[i].first;
    }
    return result;
}
 
// get magnetic field value (mV) for the given beam and its energy (e.g. parsed by ParseFileName)
// beam: beam particle name; energy: beam energy (GeV)
// return value: magnetic field value (mV), < 0 - in case of errors
double BmnFunctionSet::GetMagneticFieldByEnergy(const TString& beam, double energy)
{
    double field_value = -1;
    if (beam == "Xe") {
        // comparison tolerance between doubles
        const Double_t kEnergyTolerance = 1e-4;
        // known energy -> field value in mV
        const vector<pair<Double_t, Double_t>> energy_to_field = {{1.6, 51.5}, {2.2, 64.3}, {3.0, 85.0}, {3.8, 103.8}};
 
        if ((energy < 1) || (energy > 6)) {
            LOG(error) << "BmnFunctionSet::GetMagneticFieldByEnergy: defined energy is out of range: " << energy;
            return -3;
        }
 
        // compare with known energies (in tolerance) or use approximate formula if energy not found in the table
        size_t i = 0;
        for (; i < energy_to_field.size(); ++i) {
            const auto& entry = energy_to_field[i];
 
            // check the exact energy value from the known table
            if (TMath::Abs(energy - entry.first) < kEnergyTolerance) {
                field_value = entry.second;
                break;
            }
 
            if (entry.first > energy) {
                LOG(warning) << "BmnFunctionSet::GetMagneticFieldByEnergy: " << energy
                             << " GeV is not table value, using approximate formula for the magnetic field";
                // get nearest energy value from the table
                const auto& ref =
                    (i > 0 && TMath::Abs(energy_to_field[i - 1].first - energy) < TMath::Abs(entry.first - energy))
                        ? energy_to_field[i - 1]
                        : entry;
 
                // formula to calculate magnetic field using base (known) "energy - magnetic field" (ignoring B0)
                field_value =
                    ref.second * TMath::Sqrt(energy * (energy + 1.862)) / TMath::Sqrt(ref.first * (ref.first + 1.862));
                break;
            }
        }
 
        // if given energy greater than maximum value in the table
        if (i == energy_to_field.size()) {
            LOG(warning) << "BmnFunctionSet::GetMagneticFieldByEnergy: " << energy
                         << " GeV is not table value, using approximate formula for the magnetic field";
            const auto& entry = energy_to_field[i - 1];
            // formula to calculate magnetic field using base (known) "energy - magnetic field" (ignoring B0)
            field_value = entry.second * TMath::Sqrt(energy * (energy + 1.862))
                          / TMath::Sqrt(entry.first * (entry.first + 1.862));
        }
    } else
        LOG(warning) << "BmnFunctionSet::GetMagneticFieldByEnergy: field value is not defined for the beam: " << beam;
 
    return field_value;
}
 
// get the correction factor for the magnetic field map scale calculated as 'field voltage / map voltage'
// energy: beam energy (GeV)
// return value: correction factor to multiply the field map scale by, 1 - for the energies without correction
double BmnFunctionSet::GetFieldScaleCorrection(double energy)
{
    // comparison tolerance between doubles
    const Double_t kEnergyTolerance = 1e-4;
    // known energy (GeV) -> correction factor of the field map scale (AZ-070826)
    const vector<pair<Double_t, Double_t>> energy_to_correction = {{1.6, 1.01}, {2.2, 1. / 1.005}, {3.0, 1. / 1.02}};
 
    for (const auto& entry : energy_to_correction)
        if (TMath::Abs(energy - entry.first) < kEnergyTolerance) {
            LOG(debug) << "BmnFunctionSet::GetFieldScaleCorrection: the field map scale is corrected by "
                       << entry.second << " for the beam energy: " << energy << " GeV";
            return entry.second;
        }
 
    // the correction is known for several energies only, the field map scale is not changed for the others
    return 1.;
}
 
// parsing file name to define beam, target and energy of the run
// the file name is parsed with the '_' separator, if it does not fit then with the '.' separator
// file_name: name of the file to parse
// result (beam, target, energy (GeV)) returns via links, energy is < 0 in case of errors
// return value: 0 - success, < 0 - error codes
Int_t BmnFunctionSet::ParseFileName(const TString& file_name, TString& beam, TString& target, Double_t& energy)
{
    Int_t rc = ParseFileNameBySeparator(file_name, "_", beam, target, energy);
    if (rc == 0)
        return 0;
 
    // the dot can be used as a separator in the file name, then it also breaks a float energy value into 2 parts
    LOG(debug) << "BmnFunctionSet::ParseFileName: trying to use '.' separator instead of '_' to parse the file name: "
               << file_name;
    rc = ParseFileNameBySeparator(file_name, ".", beam, target, energy);
    if (rc != 0) {
        LOG(error) << "BmnFunctionSet::ParseFileName: beam, target and energy were not defined by the file name: "
                   << file_name;
        // clean the values, which could be defined by the partial parsing
        beam = "";
        target = "";
    }
 
    return rc;
}
 
// parsing the geometry tag, which can contain the path to an XML file with a local customization of the setup
// geometry after the ';' separator, e.g. "run9_26.09.0;$VMCWORKDIR/my_setup.xml"
Int_t BmnFunctionSet::ParseGeometryTag(TString& geometry_tag, TString& local_settings)
{
    local_settings = "";
    if (!geometry_tag.Contains(';'))
        return 0;
 
    TString setup_tag = "", unexpected_part = "";
    TObjArray* tag_tokens = geometry_tag.Tokenize(";");
    for (Int_t i = 0; i < tag_tokens->GetEntriesFast(); i++) {
        TString token = ((TObjString*)tag_tokens->At(i))->GetString().Strip(TString::kBoth);
        if (token.IsNull())
            continue;
 
        bool is_xml_file = token.EndsWith(".xml", TString::kIgnoreCase);
        if (is_xml_file && local_settings.IsNull())
            local_settings = token;
        else if (!is_xml_file && setup_tag.IsNull())
            setup_tag = token;
        else if (unexpected_part.IsNull())   // e.g. a second XML file or a second setup tag
            unexpected_part = token;
    }
    delete tag_tokens;
 
    // the separator without an XML file after it is considered as a part of the setup tag itself
    if (local_settings.IsNull()) {
        LOG(warning) << "BmnFunctionSet::ParseGeometryTag: no XML file with a geometry customization was found after "
                        "';', the whole value is used as the setup tag: "
                     << geometry_tag;
        return 0;
    }
    if (!unexpected_part.IsNull()) {
        LOG(error) << "BmnFunctionSet::ParseGeometryTag: unexpected part '" << unexpected_part
                   << "' in the geometry tag: " << geometry_tag;
        local_settings = "";
        return -1;
    }
 
    geometry_tag = setup_tag;
    // the check expands a copy of the path, so the file can be given as, e.g., "$VMCWORKDIR/my_setup.xml",
    // while the path is returned as given to be valid on another computer
    TString local_settings_path = local_settings;
    if (CheckFileExist(local_settings_path) != 1) {
        LOG(error) << "BmnFunctionSet::ParseGeometryTag: no XML file with a geometry customization: "
                   << local_settings_path;
        local_settings = "";
        return -2;
    }
 
    return 0;
}
 
// parsing file name with the given separator (used by ParseFileName)
// separator: symbol to split the file name into tokens
// return value: 0 - success, < 0 - error codes
Int_t BmnFunctionSet::ParseFileNameBySeparator(const TString& file_name,
                                               const TString& separator,
                                               TString& beam,
                                               TString& target,
                                               Double_t& energy)
{
    energy = -1;
 
    static const vector<pair<TString, TString>> beam_to_particle = {
        {"bi", "Bi"}, {"pb", "Pb"}, {"au", "Au"}, {"xe", "Xe"}, {"ag", "Ag"},
        {"kr", "Kr"}, {"ar", "Ar"}, {"c", "C"},   {"d", "d"},   {"p", "p"}};
    static const vector<pair<TString, TString>> target_to_particle = {
        {"csi", "CsI"}, {"pb", "Pb"}, {"au", "Au"}, {"cs", "Cs"}, {"xe", "Xe"}, {"sn", "Sn"},
        {"ag", "Ag"},   {"cu", "Cu"}, {"al", "Al"}, {"c", "C"},   {"p", "p"}};
 
    // split file name by separation symbol
    TObjArray* file_tokens = file_name.Tokenize(separator);
    if (!file_tokens) {
        LOG(error) << "BmnFunctionSet::ParseFileNameBySeparator: failed to tokenize file name: " << file_name;
        return -1;
    }
    const Int_t n_tokens = file_tokens->GetEntriesFast();
    Int_t token_num = 0;
 
    // lambda: get token as TString
    auto GetToken = [&](Int_t i) -> TString { return static_cast<TObjString*>(file_tokens->At(i))->GetString(); };
 
    // parse beam and target: group #1 — beam, group #3 — target
    TPRegexp beam_target_re(TString::Format("^(?P<beam>(%s))(?P<target>(%s))$", JoinKeys(beam_to_particle).Data(),
                                            JoinKeys(target_to_particle).Data()));
 
    TString beam_key, target_key;
    Bool_t found = kFALSE;
    while (!found) {
        if (token_num >= n_tokens)
            break;
        TString token_lower = GetToken(token_num);
        token_lower.ToLower();
 
        TObjArray* match = beam_target_re.MatchS(token_lower);
        if (match && match->GetEntriesFast() >= 3) {
            beam_key = static_cast<TObjString*>(match->At(1))->GetString();
            target_key = static_cast<TObjString*>(match->At(3))->GetString();
            found = kTRUE;
        }
 
        delete match;
        ++token_num;
    }
 
    if (!found) {
        LOG(debug) << "BmnFunctionSet::ParseFileNameBySeparator: Beam and Target were not found in the file name: "
                   << file_name << " (separator: " << separator << ")";
        delete file_tokens;
        return -2;
    }
 
    for (const auto& kv : beam_to_particle) {
        if (kv.first == beam_key)
            beam = kv.second;
    }
    for (const auto& kv : target_to_particle) {
        if (kv.first == target_key)
            target = kv.second;
    }
 
    // parse energy: \d+\.?\d*
    TPRegexp energy_re("\\d+\\.?\\d*");
    found = kFALSE;
    bool isUnit = false;
    TString energy_str;
    while (!found) {
        if (token_num >= n_tokens)
            break;
        TString token = GetToken(token_num);
 
        TObjArray* energy_gr = energy_re.MatchS(token);
        if (energy_gr) {
            energy_str = static_cast<TObjString*>(energy_gr->At(0))->GetString();
            found = kTRUE;
            // check if there is GeV or AGeV
            token.ToLower();
            if (token.Contains("gev"))
                isUnit = true;
        }
 
        delete energy_gr;
        ++token_num;
    }
 
    if (!found) {
        LOG(debug) << "BmnFunctionSet::ParseFileNameBySeparator: Energy was not found in the file name: " << file_name
                   << " (separator: " << separator << ")";
        delete file_tokens;
        return -3;
    }
 
    // if dot is separator then it breaks float value of the energy into 2 parts
    if (separator == "." && token_num < n_tokens) {
        TString token = GetToken(token_num);
        TPRegexp digits_re("^\\d+");
        TObjArray* digits_gr = digits_re.MatchS(token);
        if (digits_gr && digits_gr->GetEntriesFast() > 0) {
            energy_str += ".";
            energy_str += static_cast<TObjString*>(digits_gr->At(0))->GetString();
            // check if there is GeV or AGeV
            token.ToLower();
            if (token.Contains("gev"))
                isUnit = true;
 
            delete digits_gr;
            ++token_num;
        }
    }
 
    // if there is no GeV or AGeV label (before or in the next token), then raise error
    if (!isUnit) {
        TString token = (token_num < n_tokens) ? GetToken(token_num) : TString("");
        // check if there is GeV or AGeV
        token.ToLower();
        if (!token.Contains("gev")) {
            LOG(debug) << "BmnFunctionSet::ParseFileNameBySeparator: energy units (gev) were not found in the file "
                          "name: "
                       << file_name << " (separator: " << separator << ")";
            delete file_tokens;
            return -4;
        }
    }
 
    delete file_tokens;
 
    // convert the parsed energy to the number
    const Double_t energy_value = energy_str.Atof();
    if (energy_value <= 0) {
        LOG(debug) << "BmnFunctionSet::ParseFileNameBySeparator: energy value is not correct in the file name: "
                   << file_name << " (separator: " << separator << ")";
        return -5;
    }
    energy = energy_value;
 
    return 0;
}
 
// check whether path is a directory
bool BmnFunctionSet::isDirectory(TString path)
{
    FileStat_t file_stat;
    gSystem->GetPathInfo(path.Data(), file_stat);
    if (R_ISDIR(file_stat.fMode))
        return true;
    else
        return false;
}
 
TTimeStamp BmnFunctionSet::TimePoint2TS(SysPoint p)
{
    Int_t nans = duration_cast<nanoseconds>(p.time_since_epoch()).count() % 1000000000;
    Long64_t sec = duration_cast<nanoseconds>(p.time_since_epoch()).count() / 1000000000;
    return TTimeStamp{time_t(sec), nans};
}
 
SysPoint BmnFunctionSet::TimeStamp2TP(TTimeStamp t)
{
    Int_t nans = t.GetNanoSec();
    Int_t sec = t.GetSec();
    return SysPoint{sec * 1s + nans * 1ns};
}
 
void BmnFunctionSet::PrintTimePoint(SysPoint p)
{
    LOG(info) << TimePoint2String(p);
}
 
std::string BmnFunctionSet::TimePoint2String(SysPoint p)
{
    time_t tt = SysClock::to_time_t(p);
    auto nans = duration_cast<nanoseconds>(p.time_since_epoch()).count() % 1000000000;
    std::stringstream sst;
    sst << std::put_time(std::gmtime(&tt), "%a, %d %b %Y %T %z (%Z)");   // as TTimeStamp
    //    sst << std::put_time(std::localtime(&tt), "%a, %d %b %Y %T %z (%Z)");   // as TTimeStamp
    sst << " +" << std::setw(9) << nans;
    return sst.str();
}
 
SysPoint BmnFunctionSet::TDatime2GMT(TDatime& tdatime, Int_t timezoneOffsetHours)
{
    std::tm tm = {};
    tm.tm_year = tdatime.GetYear() - 1900;
    tm.tm_mon = tdatime.GetMonth() - 1;
    tm.tm_mday = tdatime.GetDay();
    tm.tm_hour = tdatime.GetHour();
    tm.tm_min = tdatime.GetMinute();
    tm.tm_sec = tdatime.GetSecond();
    time_t t = std::time(nullptr);
    std::tm* local = std::localtime(&t);
    int offset = local->tm_gmtoff;
    return (SysClock::from_time_t(std::mktime(&tm) + offset) - 1h * timezoneOffsetHours);
}