TRestExperiment.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 "TRestExperiment.h"
 
ClassImp(TRestExperiment);
 
TRestExperiment::TRestExperiment() { Initialize(); }
 
TRestExperiment::~TRestExperiment() {}
 
TRestExperiment::TRestExperiment(const char* cfgFileName, const std::string& name)
    : TRestMetadata(cfgFileName) {
    LoadConfigFromFile(fConfigFileName, name);
}
 
void TRestExperiment::Initialize() {
    SetSectionName(this->ClassName());
 
    if (!fRandom) {
        delete fRandom;
        fRandom = nullptr;
    }
 
    fRandom = new TRandom3(fSeed);
    fSeed = fRandom->TRandom::GetSeed();
}
 
void TRestExperiment::GenerateMockDataSet(Bool_t useAverage) {
    fUseAverage = useAverage;
 
    if (!fBackground) {
        RESTError << "TRestExperiment::GenerateMockData. Background component was not initialized!"
                  << RESTendl;
        return;
    }
 
    if (fExposureTime <= 0) {
        RESTError << "The experimental exposure time has not been defined" << RESTendl;
        RESTError << "This time is required to create the mock dataset" << RESTendl;
    }
 
    Double_t meanCounts = GetBackground()->GetTotalRate() * fExposureTime * units("s");
 
    Int_t N = fRandom->Poisson(meanCounts);
    if (fUseAverage) N = (Int_t)meanCounts;
    RESTInfo << "Experiment: " << GetName() << " Generating mock dataset. Counts: " << N << RESTendl;
 
    ROOT::RDF::RNode df = fBackground->GetMonteCarloDataFrame(N);
 
    fExperimentalData.SetDataFrame(df);
    fExperimentalData.SetTotalTimeInSeconds(fExposureTime * units("s"));
 
    fExperimentalCounts = *fExperimentalData.GetDataFrame().Count();
 
    fMockData = true;
    fDataReady = true;
}
 
void TRestExperiment::SetExperimentalDataSet(const std::string& filename) {
    fExperimentalDataSet = SearchFile(filename);
    fExperimentalData.Import(fExperimentalDataSet);
 
    fExposureTime = fExperimentalData.GetTotalTimeInSeconds() / units("s");
    fExperimentalCounts = *fExperimentalData.GetDataFrame().Count();
 
    fMockData = false;
    fDataReady = true;
 
    if (!fSignal || !fBackground) {
        RESTWarning << "TRestExperiment::SetExperimentalDataSet. Signal and background components must "
                       "be available before atempt to load experimental data"
                    << RESTendl;
        fDataReady = false;
        return;
    }
 
    std::vector<std::string> columns = fExperimentalData.GetDataFrame().GetColumnNames();
    for (const auto& v : fSignal->GetVariables()) {
        if (std::find(columns.begin(), columns.end(), v) == columns.end()) {
            RESTError << "TRestExperiment::SetExperimentalDataSetFile a component variable was not found in "
                         "the dataset!"
                      << RESTendl;
            fDataReady = false;
            return;
        }
    }
}
 
void TRestExperiment::InitFromConfigFile() {
    TRestMetadata::InitFromConfigFile();
 
    int cont = 0;
    TRestMetadata* md = (TRestMetadata*)this->InstantiateChildMetadata(cont);
    while (md != nullptr) {
        if (md->InheritsFrom("TRestComponent")) {
            TRestComponent* comp = (TRestComponent*)md;
            if (ToLower(comp->GetNature()) == "background")
                fBackground = comp;
            else if (ToLower(comp->GetNature()) == "signal")
                fSignal = comp;
            else
                RESTWarning << "TRestExperiment::InitFromConfigFile. Unknown component!" << RESTendl;
        }
        cont++;
        md = (TRestMetadata*)this->InstantiateChildMetadata(cont);
    }
 
    auto ele = GetElement("addComponent");
    if (ele != nullptr) {
        std::string filename = GetParameter("filename", ele, "");
        std::string component = GetParameter("component", ele, "");
 
        if (filename.empty())
            RESTWarning
                << "TRestExperiment. There is a problem with `filename` definition inside <addComponent."
                << RESTendl;
        if (component.empty())
            RESTWarning
                << "TRestExperiment. There is a problem with `component` definition inside <addComponent."
                << RESTendl;
 
        if (TRestTools::fileExists(filename) && TRestTools::isRootFile(filename)) {
            TFile* file = TFile::Open(filename.c_str(), "READ");
            if (file != nullptr) {
                TRestComponent* comp = file->Get<TRestComponent>(component.c_str());
                if (comp) {
                    if (ToLower(comp->GetNature()) == "signal")
                        fSignal = comp;
                    else if (ToLower(comp->GetNature()) == "background")
                        fBackground = comp;
                    else
                        RESTError << "TRestExperiment::InitFromConfigFile. Component : " << component
                                  << ". Nature unknown!" << RESTendl;
                } else
                    RESTError << "TRestExperiment::InitFromConfigFile. Component : " << component
                              << " not found! File : " << filename << RESTendl;
            }
        }
    }
 
    if (fExposureTime > 0 && fExperimentalDataSet.empty()) {
        GenerateMockDataSet(fUseAverage);
    } else if (fExposureTime == 0 && !fExperimentalDataSet.empty()) {
        SetExperimentalDataSet(fExperimentalDataSet);
    } else {
        RESTWarning << "The exposure time is not zero but a experimental data filename was defined!"
                    << RESTendl;
        RESTWarning << "The exposure time will be recovered from the dataset duration. To avoid confusion is"
                    << RESTendl;
        RESTWarning << "better that you set exposure time to zero inside the RML definition," << RESTendl;
        RESTError
            << " or do not define a dataset if you wish to generate mock data using the exposure time given"
            << RESTendl;
    }
 
    if (!fSignal) {
        RESTError << "TRestExperiment : " << GetName() << RESTendl;
        RESTError << "There was a problem initiazing the signal component!" << RESTendl;
        return;
    }
 
    if (!fBackground) {
        RESTError << "TRestExperiment : " << GetName() << RESTendl;
        RESTError << "There was a problem initiazing the background component!" << RESTendl;
        return;
    }
 
    if (fSignal->GetVariables() != fBackground->GetVariables()) {
        RESTError << "TRestExperiment : " << GetName() << RESTendl;
        RESTError << "Background and signal components got different variable names or variable ordering!"
                  << RESTendl;
        RESTError << "This will lead to undesired results during Likelihood calculation!" << RESTendl;
        return;
    }
 
    if (fSignal->GetNbins() != fBackground->GetNbins()) {
        RESTError << "TRestExperiment : " << GetName() << RESTendl;
        RESTError << "Background and signal components got different binning values!" << RESTendl;
        RESTError << "This will lead to undesired results during Likelihood calculation!" << RESTendl;
        return;
    }
 
    cont = 0;
    std::vector<TVector2> bckRanges = fBackground->GetRanges();
    std::vector<TVector2> sgnlRanges = fSignal->GetRanges();
    for (const TVector2& sRange : sgnlRanges) {
        if (sRange.X() != bckRanges[cont].X() || sRange.Y() != bckRanges[cont].Y()) {
            RESTError << "TRestExperiment : " << GetName() << RESTendl;
            RESTError << "Background and signal components got different range definitions!" << RESTendl;
            RESTError << "This will lead to undesired results during Likelihood calculation!" << RESTendl;
            return;
        }
        cont++;
    }
 
    Initialize();
}
 
void TRestExperiment::PrintMetadata() {
    TRestMetadata::PrintMetadata();
 
    RESTMetadata << "Random seed : " << fSeed << RESTendl;
    RESTMetadata << " " << RESTendl;
    if (fExposureTime > 0) {
        RESTMetadata << " - Exposure time : " << fExposureTime * units("s") << " seconds" << RESTendl;
        RESTMetadata << " - Exposure time : " << fExposureTime * units("hr") << " hours" << RESTendl;
        RESTMetadata << " - Exposure time : " << fExposureTime * units("day") << " days" << RESTendl;
    }
 
    if (fSignal) RESTMetadata << " - Signal component : " << fSignal->GetName() << RESTendl;
 
    if (fBackground) RESTMetadata << " - Background component : " << fBackground->GetName() << RESTendl;
 
    if (fMockData) {
        RESTMetadata << " " << RESTendl;
        if (fMockData) {
            RESTMetadata << "The dataset was MC-generated" << RESTendl;
            if (fUseAverage)
                RESTMetadata
                    << " - The number of counts in dataset was generated with the mean background counts"
                    << RESTendl;
 
        } else {
            RESTMetadata << "The dataset was loaded from an existing dataset file" << RESTendl;
            if (!fExperimentalDataSet.empty())
                RESTMetadata << " - Experimental dataset file : " << fExperimentalDataSet << RESTendl;
        }
    }
 
    RESTMetadata << " - Experimental counts : " << fExperimentalCounts << RESTendl;
 
    RESTMetadata << "----" << RESTendl;
}