TRestExperimentList.cxx

/*************************************************************************
 * This file is part of the REST software framework.                     *
 *                                                                       *
 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
 * For more information see https://gifna.unizar.es/trex                 *
 *                                                                       *
 * REST is free software: you can redistribute it and/or modify          *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * REST is distributed in the hope that it will be useful,               *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have a copy of the GNU General Public License along with   *
 * REST in $REST_PATH/LICENSE.                                           *
 * If not, see https://www.gnu.org/licenses/.                            *
 * For the list of contributors see $REST_PATH/CREDITS.                  *
 *************************************************************************/
 
#include "TRestExperimentList.h"
 
ClassImp(TRestExperimentList);
 
TRestExperimentList::TRestExperimentList() { Initialize(); }
 
TRestExperimentList::~TRestExperimentList() {}
 
TRestExperimentList::TRestExperimentList(const char* cfgFileName, const std::string& name)
    : TRestMetadata(cfgFileName) {
    LoadConfigFromFile(fConfigFileName, name);
}
 
void TRestExperimentList::Initialize() { SetSectionName(this->ClassName()); }
 
void TRestExperimentList::InitFromConfigFile() {
    TRestMetadata::InitFromConfigFile();
 
    if (!fExperimentsFile.empty() && fExperiments.empty()) {
        TRestTools::ReadASCIITable(fExperimentsFile, fExperimentsTable);
 
        for (auto& row : fExperimentsTable)
            for (auto& el : row) el = REST_StringHelper::ReplaceMathematicalExpressions(el);
 
        if (fExperimentsTable.empty()) {
            RESTError << "TRestExperimentList::InitFromConfigFile. The experiments table is empty!"
                      << RESTendl;
            return;
        }
 
        Int_t nTableColumns = fExperimentsTable[0].size();
 
        int cont = 0;
        TRestComponent* comp = (TRestComponent*)this->InstantiateChildMetadata(cont, "Component");
        while (comp != nullptr) {
            if (ToLower(comp->GetNature()) == "background")
                fBackground = comp;
            else if (ToLower(comp->GetNature()) == "signal")
                fSignal = comp;
            else
                RESTWarning << "TRestExperimentList::InitFromConfigFile. Unknown component!" << RESTendl;
 
            cont++;
            comp = (TRestComponent*)this->InstantiateChildMetadata(cont, "Component");
        }
 
        Int_t nExpectedColumns = 3;
        if (fSignal) nExpectedColumns--;
        if (fBackground) nExpectedColumns--;
        if (fExposureTime > 0) nExpectedColumns--;
 
        if (nExpectedColumns == 0) {
            RESTError << "TRestExperimentList::InitFromConfigFile. At least one free parameter required! "
                         "(Exposure/Background/Signal)"
                      << RESTendl;
            return;
        }
 
        if (nExpectedColumns != nTableColumns) {
            std::string expectedColumns = "";
            if (fExposureTime == 0) expectedColumns += "exposure";
            if (!fSignal) {
                if (fExposureTime == 0) expectedColumns += "/";
                expectedColumns += "signal";
            }
            if (!fBackground) {
                if (fExposureTime == 0 || !fSignal) expectedColumns += "/";
                expectedColumns += "background";
            }
 
            RESTError << "TRestExperimentList::InitFromConfigFile. Number of expected columns does not match "
                         "the number of table columns"
                      << RESTendl;
            RESTError << "Number of table columns : " << nTableColumns << RESTendl;
            RESTError << "Number of expected columns : " << nExpectedColumns << RESTendl;
            RESTError << "Expected columns : " << expectedColumns << RESTendl;
            return;
        }
 
        fComponentFiles = TRestTools::GetFilesMatchingPattern(fComponentPattern);
 
        Bool_t generateMockData = false;
        for (const auto& experimentRow : fExperimentsTable) {
            TRestExperiment* experiment = new TRestExperiment();
 
            std::string rowStr = "";
            for (const auto& el : experimentRow) {
                rowStr += el + " ";
            }
 
            RESTInfo << "TRestExperimentList. Loading experiment: " << rowStr << RESTendl;
 
            int column = 0;
            if (fExposureTime == 0) {
                if (REST_StringHelper::isANumber(experimentRow[column])) {
                    experiment->SetExposureInSeconds(
                        REST_StringHelper::StringToDouble(experimentRow[column]));
                    // We will generate mock data once we load the background component
                    generateMockData = true;
                } else if (TRestTools::isRootFile(experimentRow[column])) {
                    // We load the file with the dataset into the experimental data
                    std::string fname = SearchFile(experimentRow[column]);
                    experiment->SetExperimentalDataSet(fname);
                    RESTWarning << "Loading experimental data havent been tested yet!" << RESTendl;
                    RESTWarning
                        << "It might require further development. Remove these lines once it works smooth!"
                        << RESTendl;
                } else {
                    RESTError << experimentRow[column] << " is not a exposure time or an experimental dataset"
                              << RESTendl;
                }
                column++;
            } else {
                if (ToLower(fExposureStrategy) == "unique") {
                    experiment->SetExposureInSeconds(fExposureTime * units("s"));
                    // We will generate mock data once we load the background component
                    generateMockData = true;
                }
            }
 
            if (!fSignal) {
                if (GetComponent(experimentRow[column])) {
                    TRestComponent* sgnl = (TRestComponent*)GetComponent(experimentRow[column])->Clone();
                    sgnl->SetName((TString)experimentRow[column]);
                    experiment->SetSignal(sgnl);
                } else {
                    RESTError << "TRestExperimentList. Signal component : " << experimentRow[column]
                              << " not found!" << RESTendl;
                }
                column++;
            } else {
                experiment->SetSignal(fSignal);
            }
 
            if (!fBackground) {
                if (GetComponent(experimentRow[column])) {
                    TRestComponent* bck = (TRestComponent*)GetComponent(experimentRow[column])->Clone();
                    experiment->SetBackground(bck);
                } else {
                    RESTError << "TRestExperimentList. Background component : " << experimentRow[column]
                              << " not found!" << RESTendl;
                }
            } else {
                experiment->SetBackground(fBackground);
            }
 
            if (generateMockData) {
                RESTInfo << "TRestExperimentList. Generating mock dataset" << RESTendl;
                experiment->GenerateMockDataSet();
            }
 
            if (experiment->GetSignal() && experiment->GetBackground()) {
                experiment->SetName(experiment->GetSignal()->GetName());
                fExperiments.push_back(experiment);
            }
        }
 
        if (fExposureTime > 0 && ToLower(fExposureStrategy) == "ksvz") {
            ExtractExperimentParameterizationNodes();
 
            Double_t Coefficient = 0;
            for (const auto& node : fParameterizationNodes) {
                Double_t expectedRate = 0;
                for (const auto& experiment : fExperiments) {
                    Int_t nd = experiment->GetSignal()->FindActiveNode(node);
                    if (nd >= 0) {
                        experiment->GetSignal()->SetActiveNode(nd);
                        expectedRate += experiment->GetSignal()->GetTotalRate();
                    }
                }
 
                Coefficient += 1. / node / expectedRate;
            }
 
            Double_t expectedRate = 0;
            for (const auto& experiment : fExperiments) {
                // We consider the contribution to each node for a given experiment
                for (const auto& node : fParameterizationNodes) {
                    Int_t nd = experiment->GetSignal()->FindActiveNode(node);
                    if (nd >= 0) {
                        experiment->GetSignal()->SetActiveNode(nd);
                        expectedRate += node * experiment->GetSignal()->GetTotalRate();
                    }
                }
                Double_t experimentTime = fExposureTime / Coefficient / expectedRate;
                experiment->SetExposureInSeconds(experimentTime * units("s"));
            }
 
            Double_t totalExp = 0;
            for (const auto& experiment : fExperiments) totalExp += experiment->GetExposureInSeconds();
 
            for (const auto& experiment : fExperiments) {
                Double_t xp = experiment->GetExposureInSeconds();
                experiment->SetExposureInSeconds(xp * fExposureTime * units("s") / totalExp);
                experiment->GenerateMockDataSet();
            }
        }
    }
}
 
void TRestExperimentList::ExtractExperimentParameterizationNodes() {
    fParameterizationNodes.clear();
 
    for (const auto& experiment : fExperiments) {
        std::vector<Double_t> nodes = experiment->GetSignal()->GetParameterizationNodes();
        fParameterizationNodes.insert(fParameterizationNodes.end(), nodes.begin(), nodes.end());
 
        std::sort(fParameterizationNodes.begin(), fParameterizationNodes.end());
        auto last = std::unique(fParameterizationNodes.begin(), fParameterizationNodes.end());
        fParameterizationNodes.erase(last, fParameterizationNodes.end());
    }
}
 
void TRestExperimentList::PrintParameterizationNodes() {
    std::cout << "Experiment list nodes: ";
    for (const auto& node : fParameterizationNodes) std::cout << node << "\t";
    std::cout << std::endl;
}
 
TRestComponent* TRestExperimentList::GetComponent(std::string compName) {
    TRestComponent* component = nullptr;
    for (const auto& c : fComponentFiles) {
        TFile* f = TFile::Open(c.c_str(), "READ");
        TObject* obj = f->Get((TString)compName);
 
        if (!obj) continue;
 
        if (obj->InheritsFrom("TRestComponent")) {
            return (TRestComponent*)obj;
        } else {
            RESTError << "An object named : " << compName
                      << " exists inside the file, but it does not inherit from TRestComponent" << RESTendl;
        }
    }
 
    return component;
}
 
void TRestExperimentList::PrintMetadata() {
    TRestMetadata::PrintMetadata();
 
    RESTMetadata << "Number of experiments loaded: " << fExperiments.size() << RESTendl;
 
    RESTMetadata << "----" << RESTendl;
}