diff --git a/CMakeLists.txt b/CMakeLists.txt index 6d4e7928f..f58e2e85e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -188,6 +188,27 @@ else () set(REST_MPFR OFF) endif (${REST_MPFR} MATCHES "ON") +if (((${REST_ALL_LIBS} MATCHES "ON") AND (NOT DEFINED RESTLIB_AXION)) + OR (${RESTLIB_AXION} MATCHES "ON")) + # GSL #### Find GSL + find_package(GSL REQUIRED) + + # Include GSL directories + set(external_include_dirs ${external_include_dirs} ${GSL_INCLUDE_DIRS}) + + # Link GSL libraries + set(external_libs ${external_libs};${GSL_LIBRARIES}) + + message(STATUS "Found GSL libraries : ${GSL_LIBRARIES}") + message(STATUS "GSL headers : ${GSL_INCLUDE_DIRS}") + + set(REST_GSL ON) +else () + set(REST_GSL OFF) + message(STATUS "GSL libraries not required") +endif (((${REST_ALL_LIBS} MATCHES "ON") AND (NOT DEFINED RESTLIB_AXION)) + OR (${RESTLIB_AXION} MATCHES "ON")) + # CURL ##### find_library(CURL_LIB curl) if (NOT ${CURL_LIB} STREQUAL "CURL_LIB-NOTFOUND") diff --git a/doc/doxygen/images/component_hitmap.png b/doc/doxygen/images/component_hitmap.png new file mode 100644 index 000000000..ba18c0a1c Binary files /dev/null and b/doc/doxygen/images/component_hitmap.png differ diff --git a/doc/doxygen/images/component_spectra.png b/doc/doxygen/images/component_spectra.png new file mode 100644 index 000000000..d206ced4e Binary files /dev/null and b/doc/doxygen/images/component_spectra.png differ diff --git a/macros/REST_AddComponentDataSet.C b/macros/REST_AddComponentDataSet.C new file mode 100644 index 000000000..f0c60044c --- /dev/null +++ b/macros/REST_AddComponentDataSet.C @@ -0,0 +1,46 @@ +#include "TRestComponent.h" +#include "TRestTask.h" + +#ifndef RestTask_AddComponent +#define RestTask_AddComponent + +//******************************************************************************************************* +//*** Description: This macro will load from an RML the component chosen in the arguments and it +//*** will write it inside the file given as outputFile +//*** +//*** -------------- +//*** Usage: restManager AddComponentDataSet components.rml sectionName [outputFile] [componentName] [update] +//*** +//*** Arguments description: +//*** +//*** - cfgFile: The RML configuration file where the component definition can be found. +//*** - sectionName: The section name used to select a component inside the RML file. +//*** - outputFile: The file where the component is written, by default is components.root. +//*** - componentName: This argument allows to change the component name stored in the output file. +//*** By default it will take the same value as section name. +//*** - update: If disabled it will create a new file erasing any other previously added components. +//*** It is enabled by default. +//*** +//******************************************************************************************************* + +Int_t REST_AddComponentDataSet(std::string cfgFile, std::string sectionName, + std::string outputFile = "components.root", std::string componentName = "", + Bool_t update = true) { + TRestComponentDataSet comp(cfgFile.c_str(), sectionName.c_str()); + comp.Initialize(); + + TFile* f; + if (update) + f = TFile::Open(outputFile.c_str(), "UPDATE"); + else + f = TFile::Open(outputFile.c_str(), "RECREATE"); + + if (componentName == "") componentName = sectionName; + + comp.Write(componentName.c_str()); + + f->Close(); + + return 0; +} +#endif diff --git a/macros/REST_AddComponentFormula.C b/macros/REST_AddComponentFormula.C new file mode 100644 index 000000000..940866c37 --- /dev/null +++ b/macros/REST_AddComponentFormula.C @@ -0,0 +1,46 @@ +#include "TRestComponent.h" +#include "TRestTask.h" + +#ifndef RestTask_AddComponentFormula +#define RestTask_AddComponentFormula + +//******************************************************************************************************* +//*** Description: This macro will load from an RML the component chosen in the arguments and it +//*** will write it inside the file given as outputFile +//*** +//*** -------------- +//*** Usage: restManager AddComponentFormula components.rml sectionName [outputFile] [componentName] [update] +//*** +//*** Arguments description: +//*** +//*** - cfgFile: The RML configuration file where the component definition can be found. +//*** - sectionName: The section name used to select a component inside the RML file. +//*** - outputFile: The file where the component is written, by default is components.root. +//*** - componentName: This argument allows to change the component name stored in the output file. +//*** By default it will take the same value as section name. +//*** - update: If disabled it will create a new file erasing any other previously added components. +//*** It is enabled by default. +//*** +//******************************************************************************************************* + +Int_t REST_AddComponentFormula(std::string cfgFile, std::string sectionName, + std::string outputFile = "components.root", std::string componentName = "", + Bool_t update = true) { + TRestComponentFormula comp(cfgFile.c_str(), sectionName.c_str()); + comp.Initialize(); + + TFile* f; + if (update) + f = TFile::Open(outputFile.c_str(), "UPDATE"); + else + f = TFile::Open(outputFile.c_str(), "RECREATE"); + + if (componentName == "") componentName = sectionName; + + comp.Write(componentName.c_str()); + + f->Close(); + + return 0; +} +#endif diff --git a/macros/REST_CheckValidRuns.C b/macros/REST_CheckValidRuns.C index 2031a7eb9..18592be6e 100644 --- a/macros/REST_CheckValidRuns.C +++ b/macros/REST_CheckValidRuns.C @@ -37,15 +37,14 @@ namespace fs = std::filesystem; //*** CAUTION: Be aware that any non-REST file in the list will be removed if you use purge=1 //*** //******************************************************************************************************* -Int_t REST_CheckValidRuns(TString namePattern, Bool_t purge = false) { +Int_t REST_CheckValidRuns(std::string namePattern, Bool_t purge = false) { TGeoManager::SetVerboseLevel(0); vector filesNotWellClosed; TRestStringOutput RESTLog; - string a = TRestTools::Execute((string)("ls -d -1 " + namePattern)); - vector b = Split(a, "\n"); + std::vector b = TRestTools::GetFilesMatchingPattern(namePattern, true); Double_t totalTime = 0; int cont = 0; @@ -76,11 +75,12 @@ Int_t REST_CheckValidRuns(TString namePattern, Bool_t purge = false) { } RESTLog << "Run time (hours) : " << run->GetRunLength() / 3600. << RESTendl; - if (run->GetRunLength() > 0) totalTime += run->GetRunLength() / 3600.; + RESTLog << "Entries : " << run->GetEntries() << RESTendl; - if (run->GetEndTimestamp() == 0 || run->GetRunLength() < 0) { + if (run->GetEndTimestamp() == 0 || run->GetRunLength() < 0 || run->GetEntries() == 0) { filesNotWellClosed.push_back(filename); - } + } else if (run->GetRunLength() > 0) + totalTime += run->GetRunLength() / 3600.; delete run; @@ -100,7 +100,6 @@ Int_t REST_CheckValidRuns(TString namePattern, Bool_t purge = false) { if (purge) { RESTLog << "---------------------------" << RESTendl; RESTLog << "The files have been removed" << RESTendl; - RESTLog << "---------------------------" << RESTendl; } } diff --git a/macros/REST_OpenInputFile.C b/macros/REST_OpenInputFile.C index 3087791c2..7812d8f5f 100644 --- a/macros/REST_OpenInputFile.C +++ b/macros/REST_OpenInputFile.C @@ -46,6 +46,7 @@ void REST_OpenInputFile(const std::string& fileName) { printf("\n%s\n", " - dSet->PrintMetadata()"); printf("%s\n", " - dSet->GetDataFrame().GetColumnNames()"); printf("%s\n\n", " - dSet->GetTree()->GetEntries()"); + printf("%s\n", " - dSet->GetDataFrame().Display(\"\")->Print()"); printf("%s\n\n", " - dSet->GetDataFrame().Display({\"colName1,colName2\"})->Print()"); if (dSet) delete dSet; dSet = new TRestDataSet(); diff --git a/pipeline/panda-x/P3AutoChain.rml b/pipeline/panda-x/P3AutoChain.rml index 75327d268..e4b504391 100644 --- a/pipeline/panda-x/P3AutoChain.rml +++ b/pipeline/panda-x/P3AutoChain.rml @@ -1,5 +1,4 @@ - - + - - - - //change this value to "auto" to enable database - - - + //change this value to "auto" to enable database + - - - - + /// We define all observables except MinValue because it is not yet in validation chain @@ -61,31 +52,20 @@ - - + - - - - - - - - - - + diff --git a/source/bin/restRoot.cxx b/source/bin/restRoot.cxx index 17e2a1e0a..758516f47 100644 --- a/source/bin/restRoot.cxx +++ b/source/bin/restRoot.cxx @@ -61,7 +61,7 @@ int main(int argc, char* argv[]) { printf("\n"); printf(" restRoot --m [0,1]\n"); printf("\n"); - printf(" Option 0 will disable macro loading. Option 1 is the default.\n"); + printf(" Option 0 will disable macro loading. Option 0 is the default.\n"); printf("\n"); exit(0); } @@ -81,7 +81,18 @@ int main(int argc, char* argv[]) { gROOT->ProcessLine("#include "); gROOT->ProcessLine("#include "); if (loadMacros) { - if (!silent) printf("= Loading macros ...\n"); + if (!silent) { + printf("= Loading macros ...\n"); + printf( + "= HINT. Loading all macros may require a long time till you get access to the interactive " + "shell\n"); + printf("= HINT. You may use `restListMacros` outside `restRoot` to check the available macros\n"); + printf( + "= HINT. Then, you may execute the macro externally by using: `restManager MacroName " + "[ARGUMENTS]\n"); + printf("= HINT. `MacroName` is the name of the macro after removing the macro name header\n"); + printf("= HINT. E.g. REST_Detector_XYZ(arg1,arg2) may be called as: restManager XYZ arg1 arg2\n"); + } vector macroFiles; const vector patterns = { REST_PATH + "/macros/REST_*.C", // framework diff --git a/source/framework/CMakeLists.txt b/source/framework/CMakeLists.txt index 731c59cd7..b6440e4dc 100644 --- a/source/framework/CMakeLists.txt +++ b/source/framework/CMakeLists.txt @@ -4,7 +4,7 @@ link_libraries("-fPIC") add_subdirectory(external) -set(contents external/tinyxml tools core analysis masks) +set(contents external/tinyxml tools core analysis masks sensitivity) file(GLOB_RECURSE addon_src "tiny*cpp" "startup.cpp") diff --git a/source/framework/core/inc/TRestDataSet.h b/source/framework/core/inc/TRestDataSet.h index 4efe5699e..483351510 100644 --- a/source/framework/core/inc/TRestDataSet.h +++ b/source/framework/core/inc/TRestDataSet.h @@ -82,7 +82,7 @@ class TRestDataSet : public TRestMetadata { std::map fQuantity; //< /// Parameter cuts over the selected dataset - TRestCut* fCut = nullptr; + TRestCut* fCut = nullptr; //< /// The total integrated run time of selected files Double_t fTotalDuration = 0; //< @@ -103,12 +103,14 @@ class TRestDataSet : public TRestMetadata { std::vector fImportedFiles; //< /// A list of new columns together with its corresponding expressions added to the dataset - std::vector> fColumnNameExpressions; + std::vector> fColumnNameExpressions; //< /// A flag to enable Multithreading during dataframe generation Bool_t fMT = false; //< - inline auto GetAddedColumns() const { return fColumnNameExpressions; } + // If the dataframe was defined externally it will be true + Bool_t fExternal = false; //< + /// The resulting RDF::RNode object after initialization ROOT::RDF::RNode fDataSet = ROOT::RDataFrame(0); //! @@ -123,20 +125,28 @@ class TRestDataSet : public TRestMetadata { public: /// Gives access to the RDataFrame ROOT::RDF::RNode GetDataFrame() const { - if (fTree == nullptr) RESTWarning << "DataFrame has not been yet initialized" << RESTendl; + if (!fExternal && fTree == nullptr) + RESTWarning << "DataFrame has not been yet initialized" << RESTendl; return fDataSet; } - void SetDataFrame(const ROOT::RDF::RNode& dS) { fDataSet = dS; } - void EnableMultiThreading(Bool_t enable = true) { fMT = enable; } /// Gives access to the tree TTree* GetTree() const { + if (fTree == nullptr && fExternal) { + RESTInfo << "The tree is not accessible. Only GetDataFrame can be used in an externally " + "generated dataset" + << RESTendl; + RESTInfo << "You may write a tree using GetDataFrame()->Snapshot(\"MyTree\", \"output.root\");" + << RESTendl; + return fTree; + } + if (fTree == nullptr) { RESTError << "Tree has not been yet initialized" << RESTendl; - RESTError << "You should invoke TRestDataSet::GenerateDataSet() before trying to access the tree" - << RESTendl; + RESTError << "You should invoke TRestDataSet::GenerateDataSet() or " << RESTendl; + RESTError << "TRestDataSet::Import( fname ) before trying to access the tree" << RESTendl; } return fTree; } @@ -167,6 +177,7 @@ class TRestDataSet : public TRestMetadata { inline auto GetFilterLowerThan() const { return fFilterLowerThan; } inline auto GetFilterEqualsTo() const { return fFilterEqualsTo; } inline auto GetQuantity() const { return fQuantity; } + inline auto GetAddedColumns() const { return fColumnNameExpressions; } inline auto GetCut() const { return fCut; } inline auto IsMergedDataSet() const { return fMergedDataset; } @@ -174,14 +185,19 @@ class TRestDataSet : public TRestMetadata { inline void SetFilePattern(const std::string& pattern) { fFilePattern = pattern; } inline void SetQuantity(const std::map& quantity) { fQuantity = quantity; } + void SetTotalTimeInSeconds(Double_t seconds) { fTotalDuration = seconds; } + void SetDataFrame(const ROOT::RDF::RNode& dS) { + fDataSet = dS; + fExternal = true; + } + TRestDataSet& operator=(TRestDataSet& dS); Bool_t Merge(const TRestDataSet& dS); void Import(const std::string& fileName); void Import(std::vector fileNames); - void Export(const std::string& filename); + void Export(const std::string& filename, std::vector excludeColumns = {}); ROOT::RDF::RNode MakeCut(const TRestCut* cut); - ROOT::RDF::RNode DefineColumn(const std::string& columnName, const std::string& formula); void PrintMetadata() override; @@ -193,6 +209,6 @@ class TRestDataSet : public TRestMetadata { TRestDataSet(const char* cfgFileName, const std::string& name = ""); ~TRestDataSet(); - ClassDefOverride(TRestDataSet, 5); + ClassDefOverride(TRestDataSet, 7); }; #endif diff --git a/source/framework/core/inc/TRestMetadata.h b/source/framework/core/inc/TRestMetadata.h index d87a60139..561256f99 100644 --- a/source/framework/core/inc/TRestMetadata.h +++ b/source/framework/core/inc/TRestMetadata.h @@ -332,8 +332,8 @@ class TRestMetadata : public TNamed { ~TRestMetadata(); // Making class constructors protected to keep this class abstract - TRestMetadata& operator=(const TRestMetadata&) = delete; - TRestMetadata(const TRestMetadata&) = delete; + TRestMetadata& operator=(const TRestMetadata&); + TRestMetadata(const TRestMetadata&); /// Call CINT to generate streamers for this class ClassDef(TRestMetadata, 9); diff --git a/source/framework/core/inc/TRestSystemOfUnits.h b/source/framework/core/inc/TRestSystemOfUnits.h index 782979d52..5d21fe2dd 100644 --- a/source/framework/core/inc/TRestSystemOfUnits.h +++ b/source/framework/core/inc/TRestSystemOfUnits.h @@ -1,15 +1,24 @@ -///______________________________________________________________________________ -///______________________________________________________________________________ -///______________________________________________________________________________ -/// -/// -/// RESTSoft : Software for Rare Event Searches with TPCs -/// -/// TRestSystemOfUnits.h -/// -/// Feb 2016: First concept -/// author : Javier Galan -///_______________________________________________________________________________ +/************************************************************************* + * This file is part of the REST software framework. * + * * + * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza) * + * For more information see http://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 http://www.gnu.org/licenses/. * + * For the list of contributors see $REST_PATH/CREDITS. * + *************************************************************************/ #ifndef RestCore_TRestSystemOfUnits #define RestCore_TRestSystemOfUnits @@ -105,11 +114,16 @@ AddUnit(us, REST_Units::Time, 1.); AddUnit(ms, REST_Units::Time, 1.e-3); AddUnit(s, REST_Units::Time, 1.e-6); AddUnit(Hz, REST_Units::Time, 1.e6); -AddUnit(minu, REST_Units::Time, 1.67e-8); -AddUnit(hr, REST_Units::Time, 2.78e-10); -AddUnit(day, REST_Units::Time, 1.16e-11); -AddUnit(mon, REST_Units::Time, 3.85e-13); -AddUnit(yr, REST_Units::Time, 3.17e-14); +AddUnit(minu, REST_Units::Time, 1.e-6 / 60.); +AddUnit(minutes, REST_Units::Time, 1.e-6 / 60.); +AddUnit(hr, REST_Units::Time, 1e-6 / 3600.); +AddUnit(hours, REST_Units::Time, 1e-6 / 3600.); +AddUnit(day, REST_Units::Time, 1e-6 / 3600. / 24.); +AddUnit(days, REST_Units::Time, 1e-6 / 3600. / 24.); +AddUnit(mon, REST_Units::Time, 1e-6 / 3600. / 24. / 30); +AddUnit(months, REST_Units::Time, 1e-6 / 3600. / 24. / 30); +AddUnit(yr, REST_Units::Time, 1e-6 / 3600. / 24. / 365.25); +AddUnit(years, REST_Units::Time, 1e-6 / 3600. / 24. / 365.25); // length unit multiplier AddUnit(nm, REST_Units::Length, 1e6); diff --git a/source/framework/core/inc/TRestVersion.h b/source/framework/core/inc/TRestVersion.h index 3cfdf91aa..d7d1c8f66 100644 --- a/source/framework/core/inc/TRestVersion.h +++ b/source/framework/core/inc/TRestVersion.h @@ -12,12 +12,12 @@ * #endif * */ -#define REST_RELEASE "2.4.1" -#define REST_RELEASE_DATE "Sat Dec 16" -#define REST_RELEASE_TIME "11:14:21 CET 2023" -#define REST_RELEASE_NAME "Igor G. Irastorza" -#define REST_GIT_COMMIT "6b9d0650" -#define REST_VERSION_CODE 132097 +#define REST_RELEASE "2.4.2" +#define REST_RELEASE_DATE "Mon Feb 12" +#define REST_RELEASE_TIME "22:23:31 CET 2024" +#define REST_RELEASE_NAME "Henry Primakoff" +#define REST_GIT_COMMIT "d8ec95be" +#define REST_VERSION_CODE 132098 #define REST_VERSION(a, b, c) (((a) << 16) + ((b) << 8) + (c)) #define REST_SCHEMA_EVOLUTION "ON" #endif diff --git a/source/framework/core/src/TRestAnalysisTree.cxx b/source/framework/core/src/TRestAnalysisTree.cxx index 01d2c6622..d52950b24 100644 --- a/source/framework/core/src/TRestAnalysisTree.cxx +++ b/source/framework/core/src/TRestAnalysisTree.cxx @@ -623,6 +623,7 @@ Double_t TRestAnalysisTree::GetDblObservableValue(Int_t n) { } if (GetObservableType(n) == "int") return GetObservableValue(n); + if (GetObservableType(n) == "float") return GetObservableValue(n); if (GetObservableType(n) == "double") return GetObservableValue(n); cout << "TRestAnalysisTree::GetDblObservableValue. Type " << GetObservableType(n) diff --git a/source/framework/core/src/TRestDataSet.cxx b/source/framework/core/src/TRestDataSet.cxx index a3d13b94d..56011722a 100644 --- a/source/framework/core/src/TRestDataSet.cxx +++ b/source/framework/core/src/TRestDataSet.cxx @@ -391,15 +391,12 @@ std::vector TRestDataSet::FileSelection() { std::vector fileNames = TRestTools::GetFilesMatchingPattern(fFilePattern); - if (!fileNames.empty()) { - RESTInfo << "TRestDataSet::FileSelection. Starting file selection." << RESTendl; - RESTInfo << "Total files : " << fileNames.size() << RESTendl; - RESTInfo << "This process may take long computation time in case there are many files." << RESTendl; - } + RESTInfo << "TRestDataSet::FileSelection. Starting file selection." << RESTendl; + RESTInfo << "Total files : " << fileNames.size() << RESTendl; + RESTInfo << "This process may take long computation time in case there are many files." << RESTendl; fTotalDuration = 0; - std::cout << "Total files : " << fileNames.size() << std::endl; - std::cout << "Processing file selection ."; + std::cout << "Processing file selection."; int cnt = 1; for (const auto& file : fileNames) { if (cnt % 100 == 0) { @@ -413,6 +410,7 @@ std::vector TRestDataSet::FileSelection() { double runEnd = run.GetEndTimestamp(); if (runStart < time_stamp_start || runEnd > time_stamp_end) { + RESTInfo << "Rejecting file out of date range: " << file << RESTendl; continue; } @@ -783,10 +781,42 @@ void TRestDataSet::InitFromConfigFile() { /// Snapshot of the current dataset, i.e. in standard TTree format, together with a copy /// of the TRestDataSet instance that contains the conditions used to generate the dataset. /// -void TRestDataSet::Export(const std::string& filename) { +void TRestDataSet::Export(const std::string& filename, std::vector excludeColumns) { RESTInfo << "Exporting dataset" << RESTendl; + + std::vector columns = fDataSet.GetColumnNames(); + if (!excludeColumns.empty()) { + columns.erase(std::remove_if(columns.begin(), columns.end(), + [&excludeColumns](std::string elem) { + return std::find(excludeColumns.begin(), excludeColumns.end(), + elem) != excludeColumns.end(); + }), + columns.end()); + + RESTInfo << "Re-Generating snapshot." << RESTendl; + std::string user = getenv("USER"); + std::string fOutName = "/tmp/rest_output_" + user + ".root"; + fDataSet.Snapshot("AnalysisTree", fOutName, columns); + + RESTInfo << "Re-importing analysis tree." << RESTendl; + fDataSet = ROOT::RDataFrame("AnalysisTree", fOutName); + + TFile* f = TFile::Open(fOutName.c_str()); + fTree = (TChain*)f->Get("AnalysisTree"); + } + if (TRestTools::GetFileNameExtension(filename) == "txt" || TRestTools::GetFileNameExtension(filename) == "csv") { + if (excludeColumns.empty()) { + RESTInfo << "Re-Generating snapshot." << RESTendl; + std::string user = getenv("USER"); + std::string fOutName = "/tmp/rest_output_" + user + ".root"; + fDataSet.Snapshot("AnalysisTree", fOutName); + + TFile* f = TFile::Open(fOutName.c_str()); + fTree = (TChain*)f->Get("AnalysisTree"); + } + std::vector dataTypes; for (int n = 0; n < fTree->GetListOfBranches()->GetEntries(); n++) { std::string bName = fTree->GetListOfBranches()->At(n)->GetName(); @@ -876,12 +906,15 @@ void TRestDataSet::Export(const std::string& filename) { fDataSet.Snapshot("AnalysisTree", filename); TFile* f = TFile::Open(filename.c_str(), "UPDATE"); - this->Write(); + std::string name = this->GetName(); + if (name.empty()) name = "mock"; + this->Write(name.c_str()); f->Close(); } else { RESTWarning << "TRestDataSet::Export. Extension " << TRestTools::GetFileNameExtension(filename) << " not recognized" << RESTendl; } + RESTInfo << "Dataset generated: " << filename << RESTendl; } /////////////////////////////////////////////// @@ -903,6 +936,7 @@ TRestDataSet& TRestDataSet::operator=(TRestDataSet& dS) { fFilterLowerThan = dS.GetFilterLowerThan(); fFilterEqualsTo = dS.GetFilterEqualsTo(); fQuantity = dS.GetQuantity(); + fColumnNameExpressions = dS.GetAddedColumns(); fTotalDuration = dS.GetTotalTimeInSeconds(); fCut = dS.GetCut(); diff --git a/source/framework/core/src/TRestMetadata.cxx b/source/framework/core/src/TRestMetadata.cxx index 718f62daa..e5e598b41 100644 --- a/source/framework/core/src/TRestMetadata.cxx +++ b/source/framework/core/src/TRestMetadata.cxx @@ -512,6 +512,28 @@ TRestMetadata::TRestMetadata() : RESTendl(this) { #endif } +TRestMetadata::TRestMetadata(const TRestMetadata&) : RESTendl(this) { + fStore = true; + fElementGlobal = nullptr; + fElement = nullptr; + fVerboseLevel = gVerbose; + fVariables.clear(); + fConstants.clear(); + fHostmgr = nullptr; + + fConfigFileName = "null"; + configBuffer = ""; + RESTMetadata.setlength(100); + +#ifdef WIN32 + fOfficialRelease = true; + fCleanState = true; +#else + if (TRestTools::Execute("rest-config --release") == "Yes") fOfficialRelease = true; + if (TRestTools::Execute("rest-config --clean") == "Yes") fCleanState = true; +#endif +} + /////////////////////////////////////////////// /// \brief constructor /// @@ -2549,6 +2571,8 @@ void TRestMetadata::ReadOneParameter(string name, string value) { Double_t valueY = REST_Units::ConvertValueToRESTUnits(value.Y(), unit); Double_t valueZ = REST_Units::ConvertValueToRESTUnits(value.Z(), unit); *(TVector3*)datamember = TVector3(valueX, valueY, valueZ); + } else if (datamember.type == "string") { + // We just ignore this case } else { RESTWarning << this->ClassName() << " find unit definition in parameter: " << name << ", but the corresponding data member doesn't support it. Data " diff --git a/source/framework/core/src/TRestProcessRunner.cxx b/source/framework/core/src/TRestProcessRunner.cxx index 85a7efbd0..75758a91a 100644 --- a/source/framework/core/src/TRestProcessRunner.cxx +++ b/source/framework/core/src/TRestProcessRunner.cxx @@ -1099,9 +1099,7 @@ void TRestProcessRunner::PrintProcessedEvents(Int_t rateE) { // Nevents is unknown, reading root file { prog = fRunInfo->GetCurrentEntry() / (double)inputtreeentries * 100; - } - - else { + } else { prog = fProcessedEvents / (double)fEventsToProcess * 100; } diff --git a/source/framework/core/src/TRestRun.cxx b/source/framework/core/src/TRestRun.cxx index 31f4da4e4..af011ab61 100644 --- a/source/framework/core/src/TRestRun.cxx +++ b/source/framework/core/src/TRestRun.cxx @@ -515,9 +515,8 @@ void TRestRun::ReadInputFileTrees() { if (fNFilesSplit > 0) { // fNFilesSplit=1: split to 1 additional file RESTEssential << "Linking analysis tree from split data files" << RESTendl; - fAnalysisTree = - (TRestAnalysisTree*) - fAnalysisTree->Clone(); // we must make a copy to have TBrowser correctly browsed. + fAnalysisTree = (TRestAnalysisTree*)fAnalysisTree->Clone(); // we must make a copy to have + // TBrowser correctly browsed. for (int i = 1; i <= fNFilesSplit; i++) { string filename = fInputFile->GetName() + (string) "." + ToString(i); RESTInfo << filename << " --> "; @@ -1368,7 +1367,8 @@ Long64_t TRestRun::GetEntries() const { if (fAnalysisTree != nullptr) { return fAnalysisTree->GetEntries(); } - return REST_MAXIMUM_EVENTS; + return fEntriesSaved; + // return REST_MAXIMUM_EVENTS; } // Getters diff --git a/source/framework/core/src/TRestSystemOfUnits.cxx b/source/framework/core/src/TRestSystemOfUnits.cxx index ef01c59fc..ad8a58ea4 100644 --- a/source/framework/core/src/TRestSystemOfUnits.cxx +++ b/source/framework/core/src/TRestSystemOfUnits.cxx @@ -1,3 +1,25 @@ +/************************************************************************* + * This file is part of the REST software framework. * + * * + * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza) * + * For more information see http://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 http://www.gnu.org/licenses/. * + * For the list of contributors see $REST_PATH/CREDITS. * + *************************************************************************/ + #include #include @@ -45,10 +67,16 @@ map> __ListOfRESTUnits; // name, {type, scale} /// /// History of developments: /// -/// 2017-Nov: First concept and implementation of REST_Units namespace. -/// \author Javier Galan +/// 2016-Feb: First concept of REST standard system of units +/// Javier Galán +/// +/// 2017-Aug: Major upgrades +/// Kaixiang Ni /// +/// \class TRestSystemOfUnits /// \namespace REST_Units +/// \author Javier Galan +/// \author Kaixiang Ni /// ///
namespace REST_Units { @@ -280,7 +308,7 @@ TRestSystemOfUnits::TRestSystemOfUnits(string unitsStr) { orderprefix = -1; } else if (unitsStr[pos1 - 1] == '-' || unitsStr[pos1 - 1] == '*') { } else { - RESTWarning << "illegeal unit combiner \"" << unitsStr[pos1 - 1] << "\"" << RESTendl; + RESTWarning << "illegal unit combiner \"" << unitsStr[pos1 - 1] << "\"" << RESTendl; } } diff --git a/source/framework/sensitivity/inc/TRestComponent.h b/source/framework/sensitivity/inc/TRestComponent.h new file mode 100644 index 000000000..1dfd614dc --- /dev/null +++ b/source/framework/sensitivity/inc/TRestComponent.h @@ -0,0 +1,174 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef REST_TRestComponent +#define REST_TRestComponent + +#include +#include +#include + +#include +#include + +#include "TRestDataSet.h" +#include "TRestMetadata.h" +#include "TRestResponse.h" + +/// It defines a background/signal model distribution in a given parameter space (tipically x,y,en) +class TRestComponent : public TRestMetadata { + protected: + /// It defines the component type (unknown/signal/background) + std::string fNature = "unknown"; //< + + /// A list with the branches that will be used to create the distribution space + std::vector fVariables; //< + + /// The range of each of the variables used to create the PDF distribution + std::vector fRanges; //< + + /// The number of bins in which we should divide each variable + std::vector fNbins; //< + + /// It is used to parameterize a set of distribution densities (e.g. WIMP or axion mass) + std::string fParameter = ""; //< + + /// It defines the nodes of the parameterization (Initialized by the dataset) + std::vector fParameterizationNodes; //< + + /// It is used to define the node that will be accessed for rate retrieval + Int_t fActiveNode = -1; //< + + /// The generated N-dimensional variable space density for a given node + std::vector fNodeDensity; //< + + /// It introduces a fixed number of samples (if 0 it will take all available samples) + Int_t fSamples = 0; //< + + /// Enables or disables the interpolation at TRestComponentDataSet::GetRawRate + Bool_t fInterpolation = true; //< + + /// A pointer to the detector response + TRestResponse* fResponse = nullptr; //< + + /// A precision used to select the node value with a given range defined as a fraction of the value + Float_t fPrecision = 0.01; //< + + /// Internal process random generator + TRandom3* fRandom = nullptr; //! + + /// Seed used in random generator + UInt_t fSeed = 0; //< + + /// A canvas for drawing the active node component + TCanvas* fCanvas = nullptr; //! + + /// It returns true if any nodes have been defined. + Bool_t HasNodes() { return !fParameterizationNodes.empty(); } + + Bool_t HasDensity() { return !fNodeDensity.empty(); } + + /// It returns true if the node has been properly identified + Bool_t ValidNode(Double_t node) { + SetActiveNode(node); + if (GetActiveNode() >= 0) return true; + return false; + } + + Int_t GetVariableIndex(std::string varName); + + void InitFromConfigFile() override; + + virtual void FillHistograms() = 0; + + public: + void Initialize() override; + void RegenerateHistograms(UInt_t seed = 0); + + virtual void RegenerateActiveNodeDensity() {} + + std::string GetNature() const { return fNature; } + TRestResponse* GetResponse() const { return fResponse; } + Float_t GetPrecision() { return fPrecision; } + size_t GetDimensions() { return fVariables.size(); } + Int_t GetSamples() { return fSamples; } + Int_t GetActiveNode() { return fActiveNode; } + Double_t GetActiveNodeValue() { return fParameterizationNodes[fActiveNode]; } + + std::vector GetVariables() const { return fVariables; } + std::vector GetRanges() const { return fRanges; } + std::vector GetNbins() const { return fNbins; } + + Double_t GetRawRate(std::vector point); + Double_t GetTotalRate(); + Double_t GetNormalizedRate(std::vector point); + Double_t GetRate(std::vector point); + + Double_t GetBinCenter(Int_t nDim, const Int_t bin); + + void SetPrecision(const Float_t& pr) { fPrecision = pr; } + + Int_t SetActiveNode(Double_t node); + Int_t SetActiveNode(Int_t n) { + fActiveNode = n; + return fActiveNode; + } + + void SetSamples(Int_t samples) { fSamples = samples; } + + Bool_t Interpolation() { return fInterpolation; } + void EnableInterpolation() { fInterpolation = true; } + void DisableInterpolation() { fInterpolation = false; } + + THnD* GetDensityForNode(Double_t value); + THnD* GetDensityForActiveNode(); + THnD* GetDensity() { return GetDensityForActiveNode(); } + + TH1D* GetHistogram(Double_t node, std::string varName); + TH2D* GetHistogram(Double_t node, std::string varName1, std::string varName2); + TH3D* GetHistogram(Double_t node, std::string varName1, std::string varName2, std::string varName3); + + TH1D* GetHistogram(std::string varName); + TH2D* GetHistogram(std::string varName1, std::string varName2); + TH3D* GetHistogram(std::string varName1, std::string varName2, std::string varName3); + + ROOT::RVecD GetRandom(); + + ROOT::RDF::RNode GetMonteCarloDataFrame(Int_t N = 100); + + TCanvas* DrawComponent(std::vector drawVariables, std::vector scanVariables, + Int_t binScanSize = 1, TString drawOption = ""); + + void LoadResponse(const TRestResponse& resp); + + void PrintMetadata() override; + + void PrintStatistics(); + void PrintNodes(); + + TRestComponent(const char* cfgFileName, const std::string& name = ""); + TRestComponent(); + ~TRestComponent(); + + ClassDefOverride(TRestComponent, 5); +}; +#endif diff --git a/source/framework/sensitivity/inc/TRestComponentDataSet.h b/source/framework/sensitivity/inc/TRestComponentDataSet.h new file mode 100644 index 000000000..9a2154f05 --- /dev/null +++ b/source/framework/sensitivity/inc/TRestComponentDataSet.h @@ -0,0 +1,89 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef REST_TRestComponentDataSet +#define REST_TRestComponentDataSet + +#include + +#include "TRestComponent.h" +#include "TRestDataSet.h" + +/// It defines a background/signal model distribution in a given parameter space (tipically x,y,en) +class TRestComponentDataSet : public TRestComponent { + private: + /// A list with the dataset columns used to weight the distribution density and define rate + std::vector fWeights; //< + + /// It defines the number of entries in the sample for each parameterization node (Initialized by the + /// dataset) + std::vector fNSimPerNode; //< + + /// It defines the total number of entries for each parameterization node (Initialized by the dataset) + std::vector fTotalSamples; //< + + /// The filename of the dataset used + std::vector fDataSetFileNames; //< + + /// TODO we need to define multiple datasets and weigth. The weight will be used + /// to create a model, such as weighting different background contaminations or + /// different signal coupling contributions. + + /// TODO Then we probably need here a `std::vector ` and another vector + /// with the weights (isotope activity/flux/etc). + + /// The dataset used to initialize the distribution + TRestDataSet fDataSet; //! + + /// It is true of the dataset was loaded without issues + Bool_t fDataSetLoaded = false; //! + + Bool_t ValidDataSet(); + + protected: + void RegenerateActiveNodeDensity() override; + + std::vector ExtractParameterizationNodes(); + std::vector ExtractNodeStatistics(); + void FillHistograms() override; + + Bool_t VariablesOk(); + Bool_t WeightsOk(); + + Bool_t LoadDataSets(); + + public: + Bool_t IsDataSetLoaded() { return fDataSetLoaded; } + + void PrintStatistics(); + + void PrintMetadata() override; + void Initialize() override; + void InitFromConfigFile() override; + + TRestComponentDataSet(); + TRestComponentDataSet(const char* cfgFileName, const std::string& name); + ~TRestComponentDataSet(); + + ClassDefOverride(TRestComponentDataSet, 3); +}; +#endif diff --git a/source/framework/sensitivity/inc/TRestComponentFormula.h b/source/framework/sensitivity/inc/TRestComponentFormula.h new file mode 100644 index 000000000..9f9046824 --- /dev/null +++ b/source/framework/sensitivity/inc/TRestComponentFormula.h @@ -0,0 +1,57 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef REST_TRestComponentFormula +#define REST_TRestComponentFormula + +#include + +#include "TRestComponent.h" +#include "TRestDataSet.h" + +/// It defines an analytical component model distribution in a given parameter space (tipically x,y,en) +class TRestComponentFormula : public TRestComponent { + private: + /// A vector of formulas that will be added up to integrate a given rate + std::vector fFormulas; + + /// The formulas should be expressed in the following units + std::string fUnits = "cm^-2*keV^-1"; //< + + protected: + void InitFromConfigFile() override; + + void FillHistograms() override; + + public: + Double_t GetFormulaRate(std::vector point); + + void PrintMetadata() override; + + void Initialize() override; + TRestComponentFormula(const char* cfgFileName, const std::string& name); + TRestComponentFormula(); + ~TRestComponentFormula(); + + ClassDefOverride(TRestComponentFormula, 1); +}; +#endif diff --git a/source/framework/sensitivity/inc/TRestExperiment.h b/source/framework/sensitivity/inc/TRestExperiment.h new file mode 100644 index 000000000..7e25014f2 --- /dev/null +++ b/source/framework/sensitivity/inc/TRestExperiment.h @@ -0,0 +1,95 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef REST_TRestExperiment +#define REST_TRestExperiment + +#include + +#include "TRestComponent.h" +#include "TRestDataSet.h" +#include "TRestMetadata.h" + +/// It includes a model definition and experimental data used to obtain a final experimental sensitivity +class TRestExperiment : public TRestMetadata { + private: + /// The exposure time. If 0 it will be extracted from the tracking dataset (In us, standard REST unit) + Double_t fExposureTime = 0; //< + + /// A pointer to the background component + TRestComponent* fBackground = nullptr; //< + + /// A pointer to the signal component + TRestComponent* fSignal = nullptr; //< + + /// It defines the filename used to load the dataset + std::string fExperimentalDataSet = ""; + + /// It contains the experimental data (should contain same columns as the components) + TRestDataSet fExperimentalData; //< + + /// If enabled it means that the experimental data was MC-generated + Bool_t fMockData = false; //< + + /// Only if it is true we will be able to calculate the LogLikelihood + Bool_t fDataReady = false; //< + + /// Internal process random generator + TRandom3* fRandom = nullptr; //! + + /// Seed used in random generator + UInt_t fSeed = 0; //< + + protected: + void InitFromConfigFile() override; + + public: + void GenerateMockDataSet(); + + Bool_t IsMockData() const { return fMockData; } + Bool_t IsDataReady() const { return fDataReady; } + + void SetExposureInSeconds(const Double_t exposure) { fExposureTime = exposure / units("s"); } + void SetSignal(TRestComponent* comp) { fSignal = comp; } + void SetBackground(TRestComponent* comp) { fBackground = comp; } + + void SetExperimentalDataSet(const std::string& filename); + + Double_t GetExposureInSeconds() const { return fExposureTime * units("s"); } + TRestComponent* GetBackground() const { return fBackground; } + TRestComponent* GetSignal() const { return fSignal; } + TRestDataSet GetExperimentalDataSet() const { return fExperimentalData; } + ROOT::RDF::RNode GetExperimentalDataFrame() const { return fExperimentalData.GetDataFrame(); } + + void PrintExperimentalData() { GetExperimentalDataFrame().Display("")->Print(); } + + void Initialize() override; + + void PrintMetadata() override; + + TRestExperiment(const char* cfgFileName, const std::string& name = ""); + TRestExperiment(); + ~TRestExperiment(); + + ClassDefOverride(TRestExperiment, 1); +}; +#endif diff --git a/source/framework/sensitivity/inc/TRestExperimentList.h b/source/framework/sensitivity/inc/TRestExperimentList.h new file mode 100644 index 000000000..5d1737e40 --- /dev/null +++ b/source/framework/sensitivity/inc/TRestExperimentList.h @@ -0,0 +1,93 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef REST_TRestExperimentList +#define REST_TRestExperimentList + +#include "TRestExperiment.h" +#include "TRestMetadata.h" + +/// A helper metadata class to create a list of TRestExperiment instances +class TRestExperimentList : public TRestMetadata { + private: + /// A fullpath filename pattern helping to initialize the component files vector + std::string fComponentPattern = ""; //< + + /// A vector with filenames containing the components + std::vector fComponentFiles; //< + + /// A fullpath filename pattern helping to initialize the dataset files vector + std::string fDataSetPattern = ""; //< + + /// A vector with filenames containing the datasets with experimental data + std::vector fDataSetFilenames; //< + + /// A file where we define experiment components, exposureTime, and tracking data of each experiment + std::string fExperimentsFile = ""; //< Exposure/TrackingData - SignalComponent - BackgroundComponent + + /// A table with the experiment file information + std::vector > fExperimentsTable; //< + + /// A vector with a list of experiments includes the background components in this model + std::vector fExperiments; //< + + /// If not zero this will be the common exposure time in micro-seconds (standard REST units) + Double_t fExposureTime = 0; + + /// If not null this will be the common signal used in each experiment + TRestComponent* fSignal = nullptr; //< + + /// If not null this will be the common signal used in each experiment + TRestComponent* fBackground = nullptr; //< + + protected: + TRestComponent* GetComponent(std::string compName); + + void InitFromConfigFile() override; + + public: + void Initialize() override; + + void SetExposure(const Double_t& exposure) { fExposureTime = exposure; } + void SetSignal(TRestComponent* comp) { fSignal = comp; } + void SetBackground(TRestComponent* comp) { fBackground = comp; } + + std::vector GetExperiments() { return fExperiments; } + TRestExperiment* GetExperiment(const size_t& n) { + if (n >= GetNumberOfExperiments()) + return nullptr; + else + return fExperiments[n]; + } + + size_t GetNumberOfExperiments() { return fExperiments.size(); } + + void PrintMetadata() override; + + TRestExperimentList(const char* cfgFileName, const std::string& name); + + TRestExperimentList(); + ~TRestExperimentList(); + + ClassDefOverride(TRestExperimentList, 1); +}; +#endif diff --git a/source/framework/sensitivity/inc/TRestResponse.h b/source/framework/sensitivity/inc/TRestResponse.h new file mode 100644 index 000000000..25bdb4ed1 --- /dev/null +++ b/source/framework/sensitivity/inc/TRestResponse.h @@ -0,0 +1,92 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef REST_TRestResponse +#define REST_TRestResponse + +#include "TRestMetadata.h" + +/// A response matrix that might be applied to a given component inside a TRestComponent +class TRestResponse : public TRestMetadata { + private: + /// The filename used to import the response matrix + std::string fFilename = ""; + + /// It defines the variable name for which the response should be applied to + std::string fVariable = ""; + + /// First element of the response matrix (input/incident, output/detected) + TVector2 fOrigin = TVector2(0, 0); + + /// The resolution of the response matrix (binning) + Double_t fBinSize = 0.1; //< + + /// The response matrix + std::vector> fResponseMatrix; //< + + /// Determines if the response matrix has been transposed + Bool_t fTransposed = false; //< + + /// It allows to decide if the returned response should be interpolated (default:false) + Bool_t fInterpolation = false; //! + + public: + void SetBinSize(Double_t bSize) { fBinSize = bSize; } + void SetResponseFilename(std::string responseFile) { fFilename = responseFile; } + void SetOrigin(const TVector2& v) { fOrigin = v; } + void SetVariable(const std::string& var) { fVariable = var; } + + Bool_t ApplyInterpolation() { return fInterpolation; } + void Interpolate(Bool_t interpolate = true) { fInterpolation = interpolate; } + + Double_t GetBinSize() const { return fBinSize; } + std::string GetResponseFilename() const { return fFilename; } + TVector2 GetOrigin() const { return fOrigin; } + std::string GetVariable() const { return fVariable; } + + TVector2 GetInputRange() const { + return TVector2(fOrigin.X(), fOrigin.X() + fResponseMatrix[0].size() * fBinSize); + } + + TVector2 GetOutputRange() const { + return TVector2(fOrigin.Y(), fOrigin.Y() + fResponseMatrix.size() * fBinSize); + } + + void Initialize() override; + + void LoadResponse(Bool_t transpose = true); + + std::vector> GetResponse(Double_t input); + + void PrintResponseMatrix(Int_t fromRow, Int_t toRow); + + void PrintMetadata() override; + + std::vector> GetMatrix() const { return fResponseMatrix; } + + TRestResponse(const char* cfgFileName, const std::string& name = ""); + TRestResponse(); + ~TRestResponse(); + + ClassDefOverride(TRestResponse, 1); +}; +#endif diff --git a/source/framework/sensitivity/inc/TRestSensitivity.h b/source/framework/sensitivity/inc/TRestSensitivity.h new file mode 100644 index 000000000..56c9d0bd6 --- /dev/null +++ b/source/framework/sensitivity/inc/TRestSensitivity.h @@ -0,0 +1,67 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef REST_TRestSensitivity +#define REST_TRestSensitivity + +#include "TRestExperiment.h" + +/// It combines a number of experimental conditions allowing to calculate a combined experimental sensitivity +class TRestSensitivity : public TRestMetadata { + private: + /// A list of experimental conditions included to get a final sensitivity plot + std::vector fExperiments; //< + + TH1D* fSignalTest = nullptr; + + protected: + void InitFromConfigFile() override; + + Double_t UnbinnedLogLikelihood(const TRestExperiment* experiment, Double_t g4 = 0); + Double_t ApproachByFactor(Double_t g4, Double_t chi0, Double_t target, Double_t factor); + + public: + void Initialize() override; + + Double_t GetCoupling(Double_t sigma = 2, Double_t precision = 0.01); + + TH1D* SignalStatisticalTest(Int_t N); + + std::vector GetExperiments() { return fExperiments; } + TRestExperiment* GetExperiment(const size_t& n) { + if (n >= GetNumberOfExperiments()) + return nullptr; + else + return fExperiments[n]; + } + + size_t GetNumberOfExperiments() { return fExperiments.size(); } + + void PrintMetadata() override; + + TRestSensitivity(const char* cfgFileName, const std::string& name = ""); + TRestSensitivity(); + ~TRestSensitivity(); + + ClassDefOverride(TRestSensitivity, 1); +}; +#endif diff --git a/source/framework/sensitivity/src/TRestComponent.cxx b/source/framework/sensitivity/src/TRestComponent.cxx new file mode 100644 index 000000000..516bf8104 --- /dev/null +++ b/source/framework/sensitivity/src/TRestComponent.cxx @@ -0,0 +1,757 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +///////////////////////////////////////////////////////////////////////// +/// This class allows to ... +/// +/// +///---------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2023-December: First implementation of TRestComponent +/// Javier Galan +/// +/// \class TRestComponent +/// \author: Javier Galan (javier.galan.lacarra@cern.ch) +/// +///
+/// +#include "TRestComponent.h" + +#include +#include + +#include + +ClassImp(TRestComponent); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestComponent::TRestComponent() {} + +///////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// \param name The name of the specific metadata. It will be used to find the +/// corresponding TRestAxionMagneticField section inside the RML. +/// +TRestComponent::TRestComponent(const char* cfgFileName, const std::string& name) + : TRestMetadata(cfgFileName) { + RESTDebug << "Entering TRestComponent constructor( cfgFileName, name )" << RESTendl; + RESTDebug << "File: " << cfgFileName << " Name: " << name << RESTendl; +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestComponent::~TRestComponent() {} + +/////////////////////////////////////////////// +/// \brief It initializes the random number. We avoid to define the section name +/// here since we will never define a TRestComponent section in our RML file, +/// since this class is pure virtual. It will be the inherited class the +/// responsible to define the section name. +/// +void TRestComponent::Initialize() { + // SetSectionName(this->ClassName()); + + if (!fRandom) { + delete fRandom; + fRandom = nullptr; + } + + fRandom = new TRandom3(fSeed); + fSeed = fRandom->TRandom::GetSeed(); +} + +///////////////////////////////////////////// +/// \brief It will produce a histogram with the distribution defined using the +/// variables and the weights for each of the parameter nodes. +/// +/// fPrecision is used to define the active node +/// +void TRestComponent::RegenerateHistograms(UInt_t seed) { + fNodeDensity.clear(); + + fSeed = seed; + TRestComponent::Initialize(); + + FillHistograms(); +} + +/////////////////////////////////////////// +/// \brief It returns the position of the fVariable element for the variable +/// name given by argument. +/// +Int_t TRestComponent::GetVariableIndex(std::string varName) { + int n = 0; + for (const auto& v : fVariables) { + if (v == varName) return n; + n++; + } + + return -1; +} + +/////////////////////////////////////////////// +/// \brief It returns the intensity/rate (in seconds) corresponding to the +/// generated distribution or formula evaluated at the position of the parameter +/// space given by point. +/// +/// The response matrix (if defined) will be used to convolute the expected rate. +/// The TRestResponse metadata class defines the variable where the response will +/// be applied. +/// +Double_t TRestComponent::GetRate(std::vector point) { + if (!fResponse) { + return GetRawRate(point); + } + + std::string responseVariable = fResponse->GetVariable(); + Int_t respVarIndex = GetVariableIndex(responseVariable); + + if (respVarIndex == -1) { + RESTError << "The response variable `" << responseVariable << "`, defined inside TRestResponse," + << RESTendl; + RESTError << "could not be found inside the component." << RESTendl; + RESTError << "Please, check the component variable names." << RESTendl; + return 0; + } + + std::vector > response = fResponse->GetResponse(point[respVarIndex]); + + Double_t rate = 0; + for (const auto& resp : response) { + std::vector newPoint = point; + newPoint[respVarIndex] = resp.first; + rate += resp.second * GetRawRate(newPoint); + } + + return rate; +} + +/////////////////////////////////////////////// +/// \brief It returns the intensity/rate (in seconds) corresponding to the +/// generated distribution or formula evaluated at the position of the parameter +/// space given by point. +/// +/// The rate returned by the TRestComponent::GetRawRate method will be normalized +/// to the corresponding parameter space. Thus, if the parameter consists of +/// 2-spatial dimensions and 1-energy dimension, the returned rate will be +/// expressed in standard REST units as, s-1 mm-2 keV-1. +/// +/// The returned value may be recovered back with the desired units using +/// the REST_Units namespace. +/// +/// \code +/// component->GetNormalizedRate( {0,0,0} ) * units("cm^-2*keV^-1") +/// \endcode +/// +/// The response matrix (if defined) will be used to convolute the expected rate. +/// The TRestResponse metadata class defines the variable where the response will +/// be applied. +/// +Double_t TRestComponent::GetNormalizedRate(std::vector point) { + Double_t normFactor = 1; + for (size_t n = 0; n < GetDimensions(); n++) normFactor *= fNbins[n] / (fRanges[n].Y() - fRanges[n].X()); + + return normFactor * GetRate(point); +} + +/////////////////////////////////////////////// +/// \brief It returns the intensity/rate (in seconds) corresponding to the +/// generated distribution or formula evaluated at the position of the parameter +/// space given by point. The returned rate is integrated to the granularity +/// of the parameter space (cell size). To get a normalized rate use +/// TRestComponent::GetNormalizedRate. +/// +/// The size of the point vector must have the same dimension as the dimensions +/// of the distribution. +/// +/// If interpolation is enabled (which is disabled by default) the rate will be +/// evaluated using interpolation with neighbour histogram cells. +/// +/// Interpolation technique extracted from: +/// https://math.stackexchange.com/questions/1342364/formula-for-n-dimensional-linear-interpolation +/// +/// 𝑓(𝑥0,𝑥1,𝑥2)=𝐴000(1−𝑥0)(1−𝑥1)(1−𝑥2)+𝐴001𝑥0(1−𝑥1)(1−𝑥2)+𝐴010(1−𝑥0)𝑥1(1−𝑥2)⋯+𝐴111𝑥0𝑥1𝑥 +/// +Double_t TRestComponent::GetRawRate(std::vector point) { + if (point.size() != GetDimensions()) { + RESTError << "The size of the point given is : " << point.size() << RESTendl; + RESTError << "The density distribution dimensions are : " << GetDimensions() << RESTendl; + RESTError << "Point must have the same dimension as the distribution" << RESTendl; + return 0; + } + + if (!HasNodes() && !HasDensity()) { + RESTError << "TRestComponent::GetRawRate. The component has no nodes!" << RESTendl; + RESTError << "Try calling TRestComponent::Initialize" << RESTendl; + + RESTInfo << "Trying to initialize" << RESTendl; + Initialize(); + if (HasNodes()) + RESTInfo << "Sucess!" << RESTendl; + else + return 0; + } + + if (HasNodes() && fActiveNode == -1) { + RESTError << "TRestComponent::GetRawRate. Active node has not been defined" << RESTendl; + return 0; + } + + Int_t centerBin[GetDimensions()]; + Double_t centralDensity = GetDensity()->GetBinContent(GetDensity()->GetBin(point.data()), centerBin); + if (!Interpolation()) return centralDensity; + + std::vector direction; + std::vector nDist; + for (size_t dim = 0; dim < GetDimensions(); dim++) { + Double_t x1 = GetBinCenter(dim, centerBin[dim] - 1); + Double_t x2 = GetBinCenter(dim, centerBin[dim] + 1); + + if (centerBin[dim] == 1 || centerBin[dim] == fNbins[dim]) { + direction.push_back(0); + nDist.push_back(0); + } else if (x2 - point[dim] > point[dim] - x1) { + // we chose left bin (x1) since it is closer than right bin + direction.push_back(-1); + nDist.push_back(1 - 2 * (point[dim] - x1) / (x2 - x1)); + } else { + direction.push_back(1); + nDist.push_back(1 - 2 * (x2 - point[dim]) / (x2 - x1)); + } + } + + // In 3-dimensions we got 8 points to interpolate + // In 4-dimensions we would get 16 points to interpolate + // ... + Int_t nPoints = (Int_t)TMath::Power(2, (Int_t)GetDimensions()); + + Double_t sum = 0; + for (int n = 0; n < nPoints; n++) { + std::vector cell = REST_StringHelper::IntegerToBinary(n, GetDimensions()); + + Double_t weightDistance = 1; + int cont = 0; + for (const auto& c : cell) { + if (c == 0) + weightDistance *= (1 - nDist[cont]); + else + weightDistance *= nDist[cont]; + cont++; + } + + for (size_t k = 0; k < cell.size(); k++) cell[k] = cell[k] * direction[k] + centerBin[k]; + + Double_t density = GetDensity()->GetBinContent(cell.data()); + sum += density * weightDistance; + } + + return sum; +} + +/////////////////////////////////////////////// +/// \brief This method integrates the rate to all the parameter space defined in the density function. +/// The result will be returned in s-1. +/// +Double_t TRestComponent::GetTotalRate() { + THnD* dHist = GetDensityForActiveNode(); + + Double_t integral = 0; + if (dHist != nullptr) { + TH1D* h1 = dHist->Projection(0); + integral = h1->Integral(); + delete h1; + } + + return integral; +} + +/////////////////////////////////////////////// +/// \brief It returns the bin center of the given component dimension. +/// +/// It required implementation since I did not find a method inside THnD. Surprising. +/// +Double_t TRestComponent::GetBinCenter(Int_t nDim, const Int_t bin) { + return fRanges[nDim].X() + (fRanges[nDim].Y() - fRanges[nDim].X()) * ((double)bin - 0.5) / fNbins[nDim]; +} + +ROOT::RVecD TRestComponent::GetRandom() { + Double_t* tuple = new Double_t[GetDimensions()]; + + ROOT::RVecD result; + if (!GetDensity()) { + for (size_t n = 0; n < GetDimensions(); n++) result.push_back(0); + RESTWarning << "TRestComponent::GetRandom. Component might not be initialized! Use " + "TRestComponent::Initialize()." + << RESTendl; + return result; + } + + GetDensity()->GetRandom(tuple); + + for (size_t n = 0; n < GetDimensions(); n++) result.push_back(tuple[n]); + return result; +} + +ROOT::RDF::RNode TRestComponent::GetMonteCarloDataFrame(Int_t N) { + ROOT::RDF::RNode df = ROOT::RDataFrame(N); + + // Function to fill the RDataFrame using GetRandom method + auto fillRndm = [&]() { + ROOT::RVecD randomValues = GetRandom(); + return randomValues; + }; + df = df.Define("Rndm", fillRndm); + + // Creating dedicated columns for each GetRandom component + for (size_t i = 0; i < fVariables.size(); ++i) { + auto varName = fVariables[i]; + auto FillRand = [i](const ROOT::RVecD& randomValues) { return randomValues[i]; }; + df = df.Define(varName, FillRand, {"Rndm"}); + } + + /* Excluding Rndm from df */ + std::vector columns = df.GetColumnNames(); + std::vector excludeColumns = {"Rndm"}; + columns.erase(std::remove_if(columns.begin(), columns.end(), + [&excludeColumns](std::string elem) { + return std::find(excludeColumns.begin(), excludeColumns.end(), elem) != + excludeColumns.end(); + }), + columns.end()); + + std::string user = getenv("USER"); + std::string fOutName = "/tmp/rest_output_" + user + ".root"; + df.Snapshot("AnalysisTree", fOutName, columns); + + df = ROOT::RDataFrame("AnalysisTree", fOutName); + + return df; +} + +/////////////////////////////////////////////// +/// \brief A method allowing to draw a series of plots representing the density distributions. +/// +/// The method will produce 1- or 2-dimensional histograms of the `drawVariables` given in the +/// argument. A third scan variable must be provided in order to show the distribution slices +/// along the scan variable. +/// +/// The binScanSize argument can be used to define the binSize of the scanning variables. +/// +TCanvas* TRestComponent::DrawComponent(std::vector drawVariables, + std::vector scanVariables, Int_t binScanSize, + TString drawOption) { + if (drawVariables.size() > 2 || drawVariables.size() == 0) { + RESTError << "TRestComponent::DrawComponent. The number of variables to be drawn must " + "be 1 or 2!" + << RESTendl; + return fCanvas; + } + + if (scanVariables.size() > 2 || scanVariables.size() == 0) { + RESTError << "TRestComponent::DrawComponent. The number of variables to be scanned must " + "be 1 or 2!" + << RESTendl; + return fCanvas; + } + + //// Checking the number of plots to be generated + std::vector scanIndexes; + for (const auto& x : scanVariables) scanIndexes.push_back(GetVariableIndex(x)); + + Int_t nPlots = 1; + size_t n = 0; + for (const auto& x : scanIndexes) { + if (fNbins[x] % binScanSize != 0) { + RESTWarning << "The variable " << scanVariables[n] << " contains " << fNbins[x] + << " bins and it doesnt match with a bin size " << binScanSize << RESTendl; + RESTWarning << "The bin size must be a multiple of the number of bins." << RESTendl; + RESTWarning << "Redefining bin size to 1." << RESTendl; + binScanSize = 1; + } + nPlots *= fNbins[x] / binScanSize; + n++; + } + + /// Finding canvas division scheme + Int_t nPlotsX = 0; + Int_t nPlotsY = 0; + + if (scanIndexes.size() == 2) { + nPlotsX = fNbins[scanIndexes[0]] / binScanSize; + nPlotsY = fNbins[scanIndexes[1]] / binScanSize; + } else { + nPlotsX = TRestTools::CanvasDivisions(nPlots)[1]; + nPlotsY = TRestTools::CanvasDivisions(nPlots)[0]; + } + + RESTInfo << "Number of plots to be generated: " << nPlots << RESTendl; + RESTInfo << "Canvas size : " << nPlotsX << " x " << nPlotsY << RESTendl; + + //// Setting up the canvas with the appropriate number of divisions + if (fCanvas != nullptr) { + delete fCanvas; + fCanvas = nullptr; + } + + fCanvas = new TCanvas(this->GetName(), this->GetName(), 0, 0, nPlotsX * 640, nPlotsY * 480); + fCanvas->Divide(nPlotsX, nPlotsY, 0.01, 0.01); + + std::vector variableIndexes; + for (const auto& x : drawVariables) variableIndexes.push_back(GetVariableIndex(x)); + + for (int n = 0; n < nPlotsX; n++) + for (int m = 0; m < nPlotsY; m++) { + TPad* pad = (TPad*)fCanvas->cd(n * nPlotsY + m + 1); + pad->SetFixedAspectRatio(true); + + THnD* hnd = GetDensity(); + + int binXo = binScanSize * n + 1; + int binXf = binScanSize * n + binScanSize; + int binYo = binScanSize * m + 1; + int binYf = binScanSize * m + binScanSize; + + if (scanVariables.size() == 2) { + hnd->GetAxis(scanIndexes[0])->SetRange(binXo, binXf); + hnd->GetAxis(scanIndexes[1])->SetRange(binYo, binYf); + } else if (scanVariables.size() == 1) { + binXo = binScanSize * nPlotsY * n + binScanSize * m + 1; + binXf = binScanSize * nPlotsY * n + binScanSize * m + binScanSize; + hnd->GetAxis(scanIndexes[0])->SetRange(binXo, binXf); + } + + if (variableIndexes.size() == 1) { + TH1D* h1 = hnd->Projection(variableIndexes[0]); + std::string hName; + + if (scanIndexes.size() == 2) + hName = scanVariables[0] + "(" + IntegerToString(binXo) + ", " + IntegerToString(binXf) + + ") " + scanVariables[1] + "(" + IntegerToString(binYo) + ", " + + IntegerToString(binYf) + ") "; + + if (scanIndexes.size() == 1) + hName = scanVariables[0] + "(" + IntegerToString(binXo) + ", " + IntegerToString(binXf) + + ") "; + + TH1D* newh = (TH1D*)h1->Clone(hName.c_str()); + newh->SetTitle(hName.c_str()); + newh->SetStats(false); + newh->GetXaxis()->SetTitle((TString)drawVariables[0]); + newh->SetMarkerStyle(kFullCircle); + newh->Draw("PLC PMC"); + + TString entriesStr = "Entries: " + IntegerToString(newh->GetEntries()); + TLatex* textLatex = new TLatex(0.62, 0.825, entriesStr); + textLatex->SetNDC(); + textLatex->SetTextColor(1); + textLatex->SetTextSize(0.05); + textLatex->Draw("same"); + delete h1; + } + + if (variableIndexes.size() == 2) { + TH2D* h2 = hnd->Projection(variableIndexes[0], variableIndexes[1]); + + std::string hName; + if (scanIndexes.size() == 2) + hName = scanVariables[0] + "(" + IntegerToString(binXo) + ", " + IntegerToString(binXf) + + ") " + scanVariables[1] + "(" + IntegerToString(binYo) + ", " + + IntegerToString(binYf) + ") "; + + if (scanIndexes.size() == 1) + hName = scanVariables[0] + "(" + IntegerToString(binXo) + ", " + IntegerToString(binXf) + + ") "; + + TH2D* newh = (TH2D*)h2->Clone(hName.c_str()); + newh->SetStats(false); + newh->GetXaxis()->SetTitle((TString)drawVariables[0]); + newh->GetYaxis()->SetTitle((TString)drawVariables[1]); + newh->SetTitle(hName.c_str()); + newh->Draw(drawOption); + + TString entriesStr = "Entries: " + IntegerToString(newh->GetEntries()); + TLatex* textLatex = new TLatex(0.62, 0.825, entriesStr); + textLatex->SetNDC(); + textLatex->SetTextColor(1); + textLatex->SetTextSize(0.05); + textLatex->Draw("same"); + delete h2; + } + } + + return fCanvas; +} + +/////////////////////////////////////////////// +/// \brief +/// +void TRestComponent::LoadResponse(const TRestResponse& resp) { + if (fResponse) { + delete fResponse; + fResponse = nullptr; + } + + fResponse = (TRestResponse*)resp.Clone("response"); + if (fResponse) fResponse->LoadResponse(); + + fResponse->PrintMetadata(); +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux +/// +void TRestComponent::PrintMetadata() { + TRestMetadata::PrintMetadata(); + + RESTMetadata << "Component nature : " << fNature << RESTendl; + RESTMetadata << " " << RESTendl; + + RESTMetadata << "Random seed : " << fSeed << RESTendl; + if (fSamples) RESTMetadata << "Samples : " << fSamples << RESTendl; + RESTMetadata << " " << RESTendl; + + if (fVariables.size() != fRanges.size()) + RESTWarning << "The number of variables does not match with the number of defined ranges!" + << RESTendl; + + else if (fVariables.size() != fNbins.size()) + RESTWarning << "The number of variables does not match with the number of defined bins!" << RESTendl; + else { + int n = 0; + RESTMetadata << " === Variables === " << RESTendl; + for (const auto& varName : fVariables) { + RESTMetadata << " - Name: " << varName << " Range: (" << fRanges[n].X() << ", " << fRanges[n].Y() + << ") bins: " << fNbins[n] << RESTendl; + n++; + } + } + + if (!fParameter.empty()) { + RESTMetadata << " " << RESTendl; + RESTMetadata << " === Parameterization === " << RESTendl; + RESTMetadata << "- Parameter : " << fParameter << RESTendl; + + RESTMetadata << " - Number of parametric nodes : " << fParameterizationNodes.size() << RESTendl; + RESTMetadata << " " << RESTendl; + RESTMetadata << " Use : PrintNodes() for additional info" << RESTendl; + } + + if (fResponse) { + RESTMetadata << " " << RESTendl; + RESTMetadata << "A response matrix was loaded inside the component" << RESTendl; + RESTMetadata << "You may get more details using TRestComponent::GetResponse()->PrintMetadata()" + << RESTendl; + } + + RESTMetadata << "----" << RESTendl; +} + +///////////////////////////////////////////// +/// \brief It prints out on screen the values of the parametric node +/// +void TRestComponent::PrintNodes() { + std::cout << " - Number of nodes : " << fParameterizationNodes.size() << std::endl; + for (const auto& x : fParameterizationNodes) std::cout << x << " "; + std::cout << std::endl; +} + +///////////////////////////////////////////// +/// \brief It customizes the retrieval of XML data values of this class +/// +void TRestComponent::InitFromConfigFile() { + TRestMetadata::InitFromConfigFile(); + + auto ele = GetElement("cVariable"); + while (ele != nullptr) { + std::string name = GetParameter("name", ele, ""); + TVector2 v = Get2DVectorParameterWithUnits("range", ele, TVector2(-1, -1)); + Int_t bins = StringToInteger(GetParameter("bins", ele, "0")); + + if (name.empty() || (v.X() == -1 && v.Y() == -1) || bins == 0) { + RESTWarning << "TRestComponentFormula::fVariable. Problem with definition." << RESTendl; + RESTWarning << "Name: " << name << " range: (" << v.X() << ", " << v.Y() << ") bins: " << bins + << RESTendl; + } else { + fVariables.push_back(name); + fRanges.push_back(v); + fNbins.push_back(bins); + } + + ele = GetNextElement(ele); + } + + if (fNbins.size() == 0) + RESTError << "TRestComponent::InitFromConfigFile. No cVariables where found!" << RESTendl; + + if (fResponse) { + delete fResponse; + fResponse = nullptr; + } + + fResponse = (TRestResponse*)this->InstantiateChildMetadata("Response"); + if (fResponse) fResponse->LoadResponse(); +} + +///////////////////////////////////////////// +/// \brief It returns the position of the fParameterizationNodes +/// element for the variable name given by argument. +/// +Int_t TRestComponent::SetActiveNode(Double_t node) { + int n = 0; + for (const auto& v : fParameterizationNodes) { + Double_t pUp = node * (1 + fPrecision / 2); + Double_t pDown = node * (1 - fPrecision / 2); + if (v > pDown && v < pUp) { + fActiveNode = n; + return fActiveNode; + } + n++; + } + + RESTError << "Parametric node : " << node << " was not found in component" << RESTendl; + RESTError << "Keeping previous node as active : " << fParameterizationNodes[fActiveNode] << RESTendl; + PrintNodes(); + + return fActiveNode; +} + +///////////////////////////////////////////// +/// \brief +/// +THnD* TRestComponent::GetDensityForNode(Double_t node) { + int n = 0; + for (const auto& x : fParameterizationNodes) { + if (x == node) { + return fNodeDensity[n]; + } + n++; + } + + RESTError << "Parametric node : " << node << " was not found in component" << RESTendl; + PrintNodes(); + return nullptr; +} + +///////////////////////////////////////////// +/// \brief +/// +THnD* TRestComponent::GetDensityForActiveNode() { + if (fActiveNode >= 0) return fNodeDensity[fActiveNode]; + + RESTError << "The active node is invalid" << RESTendl; + PrintNodes(); + return nullptr; +} + +///////////////////////////////////////////// +/// \brief It returns a 1-dimensional projected histogram for the variable names +/// provided in the argument +/// +TH1D* TRestComponent::GetHistogram(Double_t node, std::string varName) { + SetActiveNode(node); + return GetHistogram(varName); +} + +///////////////////////////////////////////// +/// \brief It returns a 1-dimensional projected histogram for the variable names +/// provided in the argument. It will recover the histogram corresponding to +/// the active node. +/// +TH1D* TRestComponent::GetHistogram(std::string varName) { + if (fActiveNode < 0) return nullptr; + + Int_t v1 = GetVariableIndex(varName); + + if (v1 >= 0 && GetDensityForActiveNode()) return GetDensityForActiveNode()->Projection(v1); + + return nullptr; +} + +///////////////////////////////////////////// +/// \brief It returns the 2-dimensional projected histogram for the variable names +/// provided in the argument +/// +TH2D* TRestComponent::GetHistogram(Double_t node, std::string varName1, std::string varName2) { + SetActiveNode(node); + return GetHistogram(varName1, varName2); +} + +///////////////////////////////////////////// +/// \brief It returns a 2-dimensional projected histogram for the variable names +/// provided in the argument. It will recover the histogram corresponding to +/// the active node. +/// +TH2D* TRestComponent::GetHistogram(std::string varName1, std::string varName2) { + if (fActiveNode < 0) return nullptr; + + Int_t v1 = GetVariableIndex(varName1); + Int_t v2 = GetVariableIndex(varName2); + + if (v1 >= 0 && v2 >= 0) + if (GetDensityForActiveNode()) return GetDensityForActiveNode()->Projection(v1, v2); + + return nullptr; +} + +///////////////////////////////////////////// +/// \brief It returns the 3-dimensional projected histogram for the variable names +/// provided in the argument +/// +TH3D* TRestComponent::GetHistogram(Double_t node, std::string varName1, std::string varName2, + std::string varName3) { + SetActiveNode(node); + return GetHistogram(varName1, varName2, varName3); +} + +///////////////////////////////////////////// +/// \brief It returns a 3-dimensional projected histogram for the variable names +/// provided in the argument. It will recover the histogram corresponding to +/// the active node. +/// +TH3D* TRestComponent::GetHistogram(std::string varName1, std::string varName2, std::string varName3) { + if (fActiveNode < 0) return nullptr; + + Int_t v1 = GetVariableIndex(varName1); + Int_t v2 = GetVariableIndex(varName2); + Int_t v3 = GetVariableIndex(varName3); + + if (v1 >= 0 && v2 >= 0 && v3 >= 0) + if (GetDensityForActiveNode()) return GetDensityForActiveNode()->Projection(v1, v2, v3); + + return nullptr; +} diff --git a/source/framework/sensitivity/src/TRestComponentDataSet.cxx b/source/framework/sensitivity/src/TRestComponentDataSet.cxx new file mode 100644 index 000000000..951018acc --- /dev/null +++ b/source/framework/sensitivity/src/TRestComponentDataSet.cxx @@ -0,0 +1,550 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +///////////////////////////////////////////////////////////////////////// +/// This class ... +/// +/// \code +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// \endcode +/// +/// \code +/// restRoot +/// [0] TRestComponentDataSet comp("components.rml", "agSignal_vacuum"); +/// [1] comp.LoadDataSets() +/// [2] TFile *f = TFile::Open("vacuumComponent.root", "RECREATE"); +/// [3] comp.Write("agSignal_vacuum"); +/// \endcode +/// +/// \code +/// restRoot vacuumComponents.root +/// [0] TCanvas *c = agVacuum->DrawComponent( { "final_posX", "final_posY"}, {"final_energy"}, 2); +/// [1] c->Print("component_hitmaps.png"); +/// \endcode +/// +/// \htmlonly \endhtmlonly +/// ![A 2-dimensional histogram scan versus the `final_energy` observable, generated by the DrawComponent +/// method](component_hitmap.png) +/// +/// \code +/// restRoot vacuumComponents.root +/// [0] TCanvas *c = agVacuum->DrawComponent( { "final_energy"}, {"final_posX", "final_posY"}, 2); +/// [1] c->Print("component_hitmaps.png"); +/// \endcode +/// +/// In both cases each plot will regroup 2 bins. +/// +/// \htmlonly \endhtmlonly +/// ![A 1-dimensional histogram scan versus the `final_posX` and `final_posY` observables, generated by the +/// DrawComponent method](component_spectra.png) +/// +///---------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2023-December: First implementation of TRestComponentDataSet +/// Javier Galan +/// +/// \class TRestComponentDataSet +/// \author: Javier Galan (javier.galan.lacarra@cern.ch) +/// +///
+/// +#include "TRestComponentDataSet.h" + +#include + +#include + +ClassImp(TRestComponentDataSet); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestComponentDataSet::TRestComponentDataSet() {} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestComponentDataSet::~TRestComponentDataSet() {} + +///////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// \param name The name of the specific metadata. It will be used to find the +/// corresponding TRestAxionMagneticField section inside the RML. +/// +TRestComponentDataSet::TRestComponentDataSet(const char* cfgFileName, const std::string& name) + : TRestComponent(cfgFileName) { + LoadConfigFromFile(fConfigFileName, name); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); +} + +/////////////////////////////////////////////// +/// \brief It will initialize the data frame with the filelist and column names +/// (or observables) that have been defined by the user. +/// +void TRestComponentDataSet::Initialize() { + TRestComponent::Initialize(); + + SetSectionName(this->ClassName()); + + LoadDataSets(); +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux +/// +void TRestComponentDataSet::PrintMetadata() { + TRestComponent::PrintMetadata(); + + if (!fDataSetFileNames.empty()) { + RESTMetadata << " " << RESTendl; + RESTMetadata << " == Dataset filenames ==" << RESTendl; + + for (const auto& x : fDataSetFileNames) RESTMetadata << "- " << x << RESTendl; + + RESTMetadata << " " << RESTendl; + } + + if (!fParameter.empty() && fParameterizationNodes.empty()) { + RESTMetadata << "This component has no nodes!" << RESTendl; + RESTMetadata << " Use: LoadDataSets() to initialize the nodes" << RESTendl; + } + + if (!fWeights.empty()) { + RESTMetadata << " " << RESTendl; + RESTMetadata << " == Weights ==" << RESTendl; + + for (const auto& x : fWeights) RESTMetadata << "- " << x << RESTendl; + + RESTMetadata << " " << RESTendl; + } + + RESTMetadata << " Use : PrintStatistics() to check node statistics" << RESTendl; + RESTMetadata << "----" << RESTendl; +} + +///////////////////////////////////////////// +/// \brief It prints out the statistics available for each parametric node +/// +void TRestComponentDataSet::PrintStatistics() { + if (fNSimPerNode.empty() && IsDataSetLoaded()) fNSimPerNode = ExtractNodeStatistics(); + + if (!HasNodes() && !IsDataSetLoaded()) { + RESTWarning << "TRestComponentDataSet::PrintStatistics. Empty nodes and no dataset loaded!" + << RESTendl; + RESTWarning << "Invoking TRestComponentDataSet::Initialize() might solve the problem" << RESTendl; + return; + } + + auto result = std::accumulate(fNSimPerNode.begin(), fNSimPerNode.end(), 0); + RESTInfo << "Total counts : " << result << RESTendl; + std::cout << std::endl; + + RESTInfo << " Parameter node statistics (" << fParameter << ")" << RESTendl; + int n = 0; + for (const auto& p : fParameterizationNodes) { + RESTInfo << " - Value : " << p << " Counts: " << fNSimPerNode[n] << RESTendl; + n++; + } +} + +///////////////////////////////////////////// +/// \brief It customizes the retrieval of XML data values of this class +/// +void TRestComponentDataSet::InitFromConfigFile() { + TRestComponent::InitFromConfigFile(); + + auto ele = GetElement("dataset"); + while (ele != nullptr) { + fDataSetFileNames.push_back(GetParameter("filename", ele, "")); + ele = GetNextElement(ele); + } + + if (!fDataSetFileNames.empty()) Initialize(); +} + +///////////////////////////////////////////// +/// \brief It will produce a histogram with the distribution defined using the +/// variables and the weights for each of the parameter nodes. +/// +/// fPrecision is used to define the active node +/// +void TRestComponentDataSet::FillHistograms() { + if (!fNodeDensity.empty()) return; + + if (fNbins.size() == 0) { + RESTError + << "TRestComponentDataSet::FillHistograms. Trying to fill histograms but no variables found!" + << RESTendl; + return; + } + + fNSimPerNode = ExtractNodeStatistics(); + + if (!IsDataSetLoaded()) { + RESTError << "TRestComponentDataSet::FillHistograms. Dataset has not been initialized!" << RESTendl; + return; + } + + if (fParameterizationNodes.empty()) { + RESTWarning << "Nodes have not been defined" << RESTendl; + RESTWarning << "The full dataset will be used to generate the density distribution" << RESTendl; + fParameterizationNodes.push_back(-137); + } + + RESTInfo << "Generating N-dim histograms" << RESTendl; + int nIndex = 0; + for (const auto& node : fParameterizationNodes) { + Int_t from = 0; + Int_t to = 0; + if (fSamples > 0 && fTotalSamples[nIndex] - fSamples > 0) { + from = fRandom->Integer(fTotalSamples[nIndex] - fSamples); + to = from + fSamples; + fNSimPerNode[nIndex] = fSamples; + } + + ROOT::RDF::RNode df = ROOT::RDataFrame(0); + //// Yet not tested in the case when we want to define a unique node without filters + //// Needs to be improved + if (fParameterizationNodes.size() == 1 && node == -137) { + RESTInfo << "Creating component with no parameters (full dataset used)" << RESTendl; + df = fDataSet.GetDataFrame().Range(from, to); + fParameterizationNodes.clear(); + } else { + RESTInfo << "Creating THnD for parameter " << fParameter << ": " << DoubleToString(node) + << RESTendl; + Double_t pUp = node * (1 + fPrecision / 2); + Double_t pDown = node * (1 - fPrecision / 2); + std::string filter = fParameter + " < " + DoubleToString(pUp) + " && " + fParameter + " > " + + DoubleToString(pDown); + df = fDataSet.GetDataFrame().Filter(filter).Range(from, to); + } + + Int_t* bins = new Int_t[fNbins.size()]; + Double_t* xmin = new Double_t[fNbins.size()]; + Double_t* xmax = new Double_t[fNbins.size()]; + + for (size_t n = 0; n < fNbins.size(); n++) { + bins[n] = fNbins[n]; + xmin[n] = fRanges[n].X(); + xmax[n] = fRanges[n].Y(); + } + + TString hName = fParameter + "_" + DoubleToString(node); + if (fParameterizationNodes.empty()) hName = "full"; + + std::vector varsAndWeight = fVariables; + + if (!fWeights.empty()) { + std::string weightsStr = ""; + for (size_t n = 0; n < fWeights.size(); n++) { + if (n > 0) weightsStr += "*"; + + weightsStr += fWeights[n]; + } + df = df.Define("componentWeight", weightsStr); + varsAndWeight.push_back("componentWeight"); + } + + auto hn = df.HistoND({hName, hName, (int)fNbins.size(), bins, xmin, xmax}, varsAndWeight); + THnD* hNd = new THnD(*hn); + hNd->Scale(1. / fNSimPerNode[nIndex]); + + fNodeDensity.push_back(hNd); + fActiveNode = nIndex; + nIndex++; + } +} + +///////////////////////////////////////////// +/// \brief It will regenerate the density histogram for the active node. It is +/// practical in the case when the number of samples fSamples is lower than the total +/// number of samples. The density distribution will be then re-generated with a +/// different random sample. +/// +void TRestComponentDataSet::RegenerateActiveNodeDensity() { + if (fActiveNode >= 0 && fNodeDensity[fActiveNode]) { + delete fNodeDensity[fActiveNode]; + } else { + RESTError << "TRestComponentDataSet::RegenerateActiveNode. Active node undefined!" << RESTendl; + return; + } + + Int_t from = 0; + Int_t to = 0; + if (fSamples > 0 && fTotalSamples[fActiveNode] - fSamples > 0) { + from = fRandom->Integer(fTotalSamples[fActiveNode] - fSamples); + to = from + fSamples; + fNSimPerNode[fActiveNode] = fSamples; + } + + Double_t node = GetActiveNodeValue(); + RESTInfo << "Creating THnD for parameter " << fParameter << ": " << DoubleToString(node) << RESTendl; + + ROOT::RDF::RNode df = ROOT::RDataFrame(0); + Double_t pUp = node * (1 + fPrecision / 2); + Double_t pDown = node * (1 - fPrecision / 2); + std::string filter = + fParameter + " < " + DoubleToString(pUp) + " && " + fParameter + " > " + DoubleToString(pDown); + df = fDataSet.GetDataFrame().Filter(filter).Range(from, to); + + Int_t* bins = new Int_t[fNbins.size()]; + Double_t* xmin = new Double_t[fNbins.size()]; + Double_t* xmax = new Double_t[fNbins.size()]; + + for (size_t n = 0; n < fNbins.size(); n++) { + bins[n] = fNbins[n]; + xmin[n] = fRanges[n].X(); + xmax[n] = fRanges[n].Y(); + } + + TString hName = fParameter + "_" + DoubleToString(node); + if (fParameterizationNodes.empty()) hName = "full"; + + std::vector varsAndWeight = fVariables; + + if (!fWeights.empty()) { + std::string weightsStr = ""; + for (size_t n = 0; n < fWeights.size(); n++) { + if (n > 0) weightsStr += "*"; + + weightsStr += fWeights[n]; + } + df = df.Define("componentWeight", weightsStr); + varsAndWeight.push_back("componentWeight"); + } + + auto hn = df.HistoND({hName, hName, (int)fNbins.size(), bins, xmin, xmax}, varsAndWeight); + THnD* hNd = new THnD(*hn); + hNd->Scale(1. / fNSimPerNode[fActiveNode]); + + fNodeDensity[fActiveNode] = hNd; +} + +///////////////////////////////////////////// +/// \brief It returns a vector with all the different values found on +/// the dataset column for the user given parameterization variable. +/// +/// If fParameterizationNodes has already been initialized it will +/// directly return its value. +/// +std::vector TRestComponentDataSet::ExtractParameterizationNodes() { + if (!fParameterizationNodes.empty()) return fParameterizationNodes; + + RESTInfo << "Extracting parameterization nodes" << RESTendl; + + std::vector vs; + if (!IsDataSetLoaded()) { + RESTError << "TRestComponentDataSet::ExtractParameterizationNodes. Dataset has not been initialized!" + << RESTendl; + return vs; + } + + auto parValues = fDataSet.GetDataFrame().Take(fParameter); + for (const auto v : parValues) vs.push_back(v); + + std::vector::iterator ip; + ip = std::unique(vs.begin(), vs.begin() + vs.size()); + vs.resize(std::distance(vs.begin(), ip)); + std::sort(vs.begin(), vs.end()); + ip = std::unique(vs.begin(), vs.end()); + vs.resize(std::distance(vs.begin(), ip)); + + return vs; +} + +///////////////////////////////////////////// +/// \brief It returns a vector with the number of entries found for each +/// parameterization node. +/// +/// If fNSimPerNode has already been initialized it will directly return its value. +/// +/// fPrecision will be used to include a thin range where to select +/// the node values. The value defines the range with a fraction proportional to +/// the parameter value. +/// +std::vector TRestComponentDataSet::ExtractNodeStatistics() { + if (!fNSimPerNode.empty()) return fNSimPerNode; + + fTotalSamples.clear(); + + std::vector stats; + if (!IsDataSetLoaded()) { + RESTError << "TRestComponentDataSet::ExtractNodeStatistics. Dataset has not been initialized!" + << RESTendl; + return stats; + } + + RESTInfo << "Counting statistics for each node ..." << RESTendl; + RESTInfo << "Number of nodes : " << fParameterizationNodes.size() << RESTendl; + for (const auto& p : fParameterizationNodes) { + Double_t pUp = p * (1 + fPrecision / 2); + Double_t pDown = p * (1 - fPrecision / 2); + std::string filter = + fParameter + " < " + DoubleToString(pUp) + " && " + fParameter + " > " + DoubleToString(pDown); + RESTInfo << "Counting stats for : " << fParameter << " = " << p << RESTendl; + auto nEv = fDataSet.GetDataFrame().Filter(filter).Count(); + fTotalSamples.push_back(*nEv); + RESTInfo << "Total entries for " << fParameter << ":" << p << " = " << *nEv << RESTendl; + if (fSamples != 0) { + nEv = fDataSet.GetDataFrame().Filter(filter).Range(fSamples).Count(); + } + + if ((Int_t)*nEv < fSamples) { + RESTWarning << "The number of requested samples (" << fSamples + << ") is higher than the number of dataset entries (" << *nEv << ")" << RESTendl; + } + RESTInfo << "Samples to be used for " << fParameter << ":" << p << " = " << *nEv << RESTendl; + stats.push_back(*nEv); + } + return stats; +} + +///////////////////////////////////////////// +/// \brief A method responsible to import a list of TRestDataSet into fDataSet +/// and check that the variables and weights defined by the user can be found +/// inside the dataset. +/// +Bool_t TRestComponentDataSet::LoadDataSets() { + if (fDataSetFileNames.empty()) { + fDataSetLoaded = false; + return fDataSetLoaded; + } + + RESTInfo << "Loading datasets" << RESTendl; + + std::vector fullFileNames; + for (const auto& name : fDataSetFileNames) { + // TODO we get here a list of files. However, we will need to weight each dataset with a factor + // to consider the contribution of each background component. + // Of course, we could previously take a factor into account already in the dataset, through the + // definition of a new column. But being this way would allow us to play around with the + // background model without having to regenerate the dataset. + std::string fileName = SearchFile(name); + if (fileName.empty()) { + RESTError << "TRestComponentDataSet::LoadDataSet. Error loading file : " << name << RESTendl; + RESTError << "Does the file exist?" << RESTendl; + RESTError << "You may use ` = TRestStringOutput::REST_Verbose_Level::REST_Info) fDataSet.PrintMetadata(); + + if (VariablesOk() && WeightsOk()) { + fParameterizationNodes = ExtractParameterizationNodes(); + FillHistograms(); + return fDataSetLoaded; + } + + return fDataSetLoaded; +} + +///////////////////////////////////////////// +/// \brief It returns true if all variables have been found inside TRestDataSet +/// +Bool_t TRestComponentDataSet::VariablesOk() { + Bool_t ok = true; + std::vector cNames = fDataSet.GetDataFrame().GetColumnNames(); + + for (const auto& var : fVariables) + if (std::count(cNames.begin(), cNames.end(), var) == 0) { + RESTError << "Variable ---> " << var << " <--- NOT found on dataset" << RESTendl; + ok = false; + } + return ok; +} + +///////////////////////////////////////////// +/// \brief It returns true if all weights have been found inside TRestDataSet +/// +Bool_t TRestComponentDataSet::WeightsOk() { + Bool_t ok = true; + std::vector cNames = fDataSet.GetDataFrame().GetColumnNames(); + + for (const auto& var : fWeights) + if (std::count(cNames.begin(), cNames.end(), var) == 0) { + RESTError << "Weight ---> " << var << " <--- NOT found on dataset" << RESTendl; + ok = false; + } + return ok; +} + +///////////////////////////////////////////// +/// \brief Takes care of initializing datasets if have not been initialized. +/// On sucess it returns true. +/// +Bool_t TRestComponentDataSet::ValidDataSet() { + if (!IsDataSetLoaded()) { + RESTWarning << "TRestComponentDataSet::ValidDataSet. Dataset has not been loaded" << RESTendl; + RESTWarning << "Try calling TRestComponentDataSet::Initialize()" << RESTendl; + + RESTInfo << "Trying to load datasets" << RESTendl; + LoadDataSets(); + if (IsDataSetLoaded()) { + RESTInfo << "Sucess!" << RESTendl; + } else { + RESTError << "Failed loading datasets" << RESTendl; + return false; + } + } + + if (HasNodes() && fActiveNode == -1) { + RESTError << "TRestComponentDataSet::ValidDataSet. Active node has not been defined" << RESTendl; + return false; + } + return true; +} diff --git a/source/framework/sensitivity/src/TRestComponentFormula.cxx b/source/framework/sensitivity/src/TRestComponentFormula.cxx new file mode 100644 index 000000000..06c8ef3d6 --- /dev/null +++ b/source/framework/sensitivity/src/TRestComponentFormula.cxx @@ -0,0 +1,260 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +///////////////////////////////////////////////////////////////////////// +/// This class allows to ... +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// // A spatial gaussian component contribution (energy normalized) +/// +/// // A flat contribution +/// +/// +/// +///---------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2023-December: First implementation of TRestComponentFormula +/// Javier Galan +/// +/// \class TRestComponentFormula +/// \author: Javier Galan (javier.galan.lacarra@cern.ch) +/// +///
+/// +#include "TRestComponentFormula.h" + +#include + +#include "TKey.h" + +ClassImp(TRestComponentFormula); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestComponentFormula::TRestComponentFormula() { Initialize(); } + +///////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// \param name The name of the specific metadata. It will be used to find the +/// corresponding TRestAxionMagneticField section inside the RML. +/// +TRestComponentFormula::TRestComponentFormula(const char* cfgFileName, const std::string& name) + : TRestComponent(cfgFileName) { + Initialize(); + + LoadConfigFromFile(fConfigFileName, name); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestComponentFormula::~TRestComponentFormula() {} + +/////////////////////////////////////////////// +/// \brief It will initialize the data frame with the filelist and column names +/// (or observables) that have been defined by the user. +/// +void TRestComponentFormula::Initialize() { + TRestComponent::Initialize(); + + SetSectionName(this->ClassName()); + + FillHistograms(); +} + +/////////////////////////////////////////////// +/// \brief It returns the intensity/rate (in seconds) corresponding to the +/// formula evaluated at the position of the parameter space given by point +/// and integrated to the parameter space cell volume. +/// +/// The size of the point vector must have the same dimension as the dimensions +/// of the variables of the distribution. +/// +Double_t TRestComponentFormula::GetFormulaRate(std::vector point) { + if (fVariables.size() != point.size()) { + RESTError << "Point should have same dimensions as number of variables!" << RESTendl; + return 0; + } + + Double_t result = 0; + for (auto& formula : fFormulas) { + for (size_t n = 0; n < fVariables.size(); n++) formula.SetParameter(fVariables[n].c_str(), point[n]); + + result += formula.EvalPar(nullptr); + } + + Double_t normFactor = 1; + for (size_t n = 0; n < GetDimensions(); n++) { + normFactor *= (fRanges[n].Y() - fRanges[n].X()) / fNbins[n]; + } + + return normFactor * result / units(fUnits); +} + +///////////////////////////////////////////// +/// \brief It will produce a histogram with the distribution using the formula +/// contributions. +/// +/// For the moment this method will just fill one node (without fParameter). But +/// if the component expression depends on the node parameter it might require +/// further development. +/// +/// TODO: The histogram is filled just by evaluating the formula, but it would +/// be more realistic that we fill the histograms with a number N of entries +/// that mimic a MC generation scheme similar to TRestComponentDataSet. +/// +void TRestComponentFormula::FillHistograms() { + if (fFormulas.empty()) return; + + RESTInfo << "Generating N-dim histogram for " << GetName() << RESTendl; + + TString hName = "formula"; + + Int_t* bins = new Int_t[fNbins.size()]; + Double_t* xlow = new Double_t[fNbins.size()]; + Double_t* xhigh = new Double_t[fNbins.size()]; + + for (size_t n = 0; n < fNbins.size(); n++) { + bins[n] = fNbins[n]; + xlow[n] = fRanges[n].X(); + xhigh[n] = fRanges[n].Y(); + } + + THnD* hNd = new THnD(hName, hName, fNbins.size(), bins, xlow, xhigh); + + // Calculate the bin width in each dimension + std::vector binWidths; + for (size_t i = 0; i < fNbins.size(); ++i) { + double width = static_cast(xhigh[i] - xlow[i]) / bins[i]; + binWidths.push_back(width); + } + + // Nested loop to iterate over each bin and print its center + std::vector binIndices(fNbins.size(), 0); // Initialize bin indices to 0 in each dimension + + bool carry = false; + while (!carry) { + // Calculate the center of the current bin in each dimension + std::vector binCenter; + for (size_t i = 0; i < fNbins.size(); ++i) + binCenter.push_back(xlow[i] + (binIndices[i] + 0.5) * binWidths[i]); + + hNd->Fill(binCenter.data(), GetFormulaRate(binCenter)); + + // Update bin indices for the next iteration + carry = true; + for (size_t i = 0; i < fNbins.size(); ++i) { + binIndices[i]++; + if (binIndices[i] < bins[i]) { + carry = false; + break; + } + binIndices[i] = 0; + } + } + + fNodeDensity.clear(); + fNodeDensity.push_back(hNd); + fActiveNode = 0; // For the moment only 1-node! +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux +/// +void TRestComponentFormula::PrintMetadata() { + TRestComponent::PrintMetadata(); + + RESTMetadata << " " << RESTendl; + RESTMetadata << "Formula units: " << fUnits << RESTendl; + + if (!fFormulas.empty()) { + RESTMetadata << " " << RESTendl; + RESTMetadata << " == Contributions implemented inside the component ==" << RESTendl; + + for (const auto& x : fFormulas) + RESTMetadata << "- " << x.GetName() << " = " << x.GetExpFormula() << RESTendl; + + RESTMetadata << " " << RESTendl; + } + + RESTMetadata << "----" << RESTendl; +} + +///////////////////////////////////////////// +/// \brief It customizes the retrieval of XML data values of this class +/// +void TRestComponentFormula::InitFromConfigFile() { + TRestComponent::InitFromConfigFile(); + + if (!fFormulas.empty()) return; + + auto ele = GetElement("formula"); + while (ele != nullptr) { + std::string name = GetParameter("name", ele, ""); + std::string expression = GetParameter("expression", ele, ""); + + if (expression.empty()) { + RESTWarning << "TRestComponentFormula::InitFromConfigFile. Invalid formula" << RESTendl; + } else { + TFormula formula(name.c_str(), expression.c_str()); + + for (Int_t n = 0; n < formula.GetNpar(); n++) { + if (std::find(fVariables.begin(), fVariables.end(), formula.GetParName(n)) == + fVariables.end()) { + RESTError << "Variable : " << formula.GetParName(n) << " not found in component! " + << RESTendl; + RESTError << "TRestComponentFormula evaluation will lead to wrong results!" << RESTendl; + } + } + + for (const auto& varName : fVariables) formula.SetParameter(varName.c_str(), 0.0); + + fFormulas.push_back(formula); + } + + ele = GetNextElement(ele); + } +} diff --git a/source/framework/sensitivity/src/TRestExperiment.cxx b/source/framework/sensitivity/src/TRestExperiment.cxx new file mode 100644 index 000000000..674c391c9 --- /dev/null +++ b/source/framework/sensitivity/src/TRestExperiment.cxx @@ -0,0 +1,289 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +///////////////////////////////////////////////////////////////////////// +/// Documentation TOBE written +/// +/// +///---------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2022-December: First implementation of TRestExperiment +/// Javier Galan +/// +/// \class TRestExperiment +/// \author: Javier Galan (javier.galan.lacarra@cern.ch) +/// +///
+/// +#include "TRestExperiment.h" + +ClassImp(TRestExperiment); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestExperiment::TRestExperiment() { Initialize(); } + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestExperiment::~TRestExperiment() {} + +///////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// \param name The name of the specific metadata. It will be used to find the +/// corresponding TRestAxionMagneticField section inside the RML. +/// +TRestExperiment::TRestExperiment(const char* cfgFileName, const std::string& name) + : TRestMetadata(cfgFileName) { + LoadConfigFromFile(fConfigFileName, name); +} + +/////////////////////////////////////////////// +/// \brief It will initialize the data frame with the filelist and column names +/// (or observables) that have been defined by the user. +/// +void TRestExperiment::Initialize() { + SetSectionName(this->ClassName()); + + if (!fRandom) { + delete fRandom; + fRandom = nullptr; + } + + fRandom = new TRandom3(fSeed); + fSeed = fRandom->TRandom::GetSeed(); +} + +void TRestExperiment::GenerateMockDataSet() { + 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); + + ROOT::RDF::RNode df = fBackground->GetMonteCarloDataFrame(N); + + fExperimentalData.SetDataFrame(df); + fExperimentalData.SetTotalTimeInSeconds(fExposureTime * units("s")); + + fMockData = true; + fDataReady = true; +} + +void TRestExperiment::SetExperimentalDataSet(const std::string& filename) { + fExperimentalDataSet = SearchFile(filename); + fExperimentalData.Import(fExperimentalDataSet); + + /// fExposureTime is in standard REST units : us + fExposureTime = fExperimentalData.GetTotalTimeInSeconds() / units("s"); + + 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 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; + } + } +} + +///////////////////////////////////////////// +/// \brief It customizes the retrieval of XML data values of this class +/// +void TRestExperiment::InitFromConfigFile() { + TRestMetadata::InitFromConfigFile(); + + 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 << "TRestExperiment::InitFromConfigFile. Unknown component!" << RESTendl; + + cont++; + comp = (TRestComponent*)this->InstantiateChildMetadata(cont, "Component"); + } + + 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 Get(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(); + } 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; + } + + /// Checking that signal/background/tracking got the same variable names and ranges + 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 bckRanges = fBackground->GetRanges(); + std::vector 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(); +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux +/// +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; + else { + RESTMetadata << "The dataset was loaded from an existing dataset file" << RESTendl; + if (!fExperimentalDataSet.empty()) + RESTMetadata << " - Experimental dataset file : " << fExperimentalDataSet << RESTendl; + } + } + + RESTMetadata << " - Experimental counts : " << *fExperimentalData.GetDataFrame().Count() << RESTendl; + + RESTMetadata << "----" << RESTendl; +} diff --git a/source/framework/sensitivity/src/TRestExperimentList.cxx b/source/framework/sensitivity/src/TRestExperimentList.cxx new file mode 100644 index 000000000..ab773d024 --- /dev/null +++ b/source/framework/sensitivity/src/TRestExperimentList.cxx @@ -0,0 +1,248 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +///////////////////////////////////////////////////////////////////////// +/// Documentation TOBE written +/// +/// +///---------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2022-December: First implementation of TRestExperimentList +/// Javier Galan +/// +/// \class TRestExperimentList +/// \author: Javier Galan (javier.galan.lacarra@cern.ch) +/// +///
+/// +#include "TRestExperimentList.h" + +ClassImp(TRestExperimentList); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestExperimentList::TRestExperimentList() { Initialize(); } + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestExperimentList::~TRestExperimentList() {} + +///////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// \param name The name of the specific metadata. +/// +TRestExperimentList::TRestExperimentList(const char* cfgFileName, const std::string& name) + : TRestMetadata(cfgFileName) { + LoadConfigFromFile(fConfigFileName, name); +} + +/////////////////////////////////////////////// +/// \brief It will initialize the data frame with the filelist and column names +/// (or observables) that have been defined by the user. +/// +void TRestExperimentList::Initialize() { SetSectionName(this->ClassName()); } + +///////////////////////////////////////////// +/// \brief It customizes the retrieval of XML data values of this class +/// +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 { + 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(); + 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(); + } + + fExperiments.push_back(experiment); + } + } +} + +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; +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux +/// +void TRestExperimentList::PrintMetadata() { + TRestMetadata::PrintMetadata(); + + RESTMetadata << "Number of experiments loaded: " << fExperiments.size() << RESTendl; + + RESTMetadata << "----" << RESTendl; +} diff --git a/source/framework/sensitivity/src/TRestResponse.cxx b/source/framework/sensitivity/src/TRestResponse.cxx new file mode 100644 index 000000000..1c9c769bf --- /dev/null +++ b/source/framework/sensitivity/src/TRestResponse.cxx @@ -0,0 +1,248 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +///////////////////////////////////////////////////////////////////////// +/// Documentation TOBE written +/// +/// +///---------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2023-December: First implementation of TRestResponse +/// Javier Galan +/// +/// \class TRestResponse +/// \author: Javier Galan (javier.galan.lacarra@cern.ch) +/// +///
+/// +#include "TRestResponse.h" + +ClassImp(TRestResponse); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestResponse::TRestResponse() { Initialize(); } + +///////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// \param name The name of the specific metadata. It will be used to find the +/// corresponding TRestAxionMagneticField section inside the RML. +/// +TRestResponse::TRestResponse(const char* cfgFileName, const std::string& name) : TRestMetadata(cfgFileName) { + Initialize(); + + LoadConfigFromFile(fConfigFileName, name); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestResponse::~TRestResponse() {} + +/////////////////////////////////////////////// +/// \brief It will initialize the data frame with the filelist and column names +/// (or observables) that have been defined by the user. +/// +void TRestResponse::Initialize() { SetSectionName(this->ClassName()); } + +/////////////////////////////////////////////// +/// \brief It loads into the fResponseMatrix data member the response from a file. +/// +/// For the moment only binary data files with format .N???f have been implemented. +/// +/// The response file should be arranged in a way that the more internal `std::vector` +/// (row) inside the `std::vector ` table corresponds to a specific +/// transformed energy (output). For example, if we have the incident particle energy +/// (input) and the expected detected energy response (output) for that energy, then +/// each row in the matrix corresponds to a detected energy and each element of that +/// row defines the fraction of incident energies that contributed to that detected +/// energy. +/// +/// Thats why we may need to transpose the matrix, so that when we can extract a row +/// (detected energy) from the matrix directly, such as fMatrix[n], and we just get +/// the vector that should convolute with the Signal(Ei) that is a vector of signals +/// as a function of incident energy. The resulting scalar product will give the +/// expected signal at the detection energy. +/// +void TRestResponse::LoadResponse(Bool_t transpose) { + if (fFilename == "") { + RESTError << "TRestResponse::LoadResponse. The response filename was not defined" << RESTendl; + return; + } + + std::string fullFilename = SearchFile(fFilename); + if (fullFilename.empty()) { + RESTError << "TRestResponse::LoadResponse. The response filename was not found!" << RESTendl; + RESTError << "Filename : " << fFilename << RESTendl; + RESTError << "You may want to define definition" << RESTendl; + return; + } + + std::string extension = TRestTools::GetFileNameExtension(fFilename); + if (!extension.empty() && extension[0] == 'N' && extension.back() == 'f') { + TRestTools::ReadBinaryTable(fullFilename, fResponseMatrix); + + fTransposed = false; + if (transpose) { + fTransposed = transpose; + TRestTools::TransposeTable(fResponseMatrix); + } + + return; + } + + RESTError << "Extension format - " << extension << " - not recognized!" << RESTendl; +} + +///////////////////////////////////////////// +/// \brief This method will return a vector of std::pair, each pair will contain the +/// output energy together with the corresponding response (or efficiency), for the +/// given input energy. +/// +/// The output value will be mapped following the binning and the origin given on the +/// metadata members. +/// +std::vector> TRestResponse::GetResponse(Double_t input) { + std::vector> response; + + if (fResponseMatrix.empty()) { + RESTError << "TRestResponse::GetResponse. Response matrix has not been loaded yet!" << RESTendl; + return response; + } + + if (input < GetInputRange().X() || input > GetInputRange().Y()) { + RESTError << "TRestResponse::GetResponse. The input value " << input << " is outside range!" + << RESTendl; + return response; + } + + if (!fInterpolation) { + Int_t bin = (Int_t)((input - fOrigin.X()) / fBinSize); + + for (std::size_t n = 0; n < fResponseMatrix[bin].size(); n++) { + Double_t output = fOrigin.Y() + ((double)n + 0.5) * fBinSize; + Double_t value = fResponseMatrix[bin][n]; + + std::pair outp{output, value}; + + response.push_back(outp); + } + + return response; + } + + Int_t binLeft = (Int_t)((input - fBinSize / 2. - fOrigin.X()) / fBinSize); + Int_t binRight = binLeft + 1; + + Double_t distLeft = (input - fBinSize / 2. + fOrigin.X()) - binLeft * fBinSize; + + if (input <= GetInputRange().X() + fBinSize / 2. || input >= GetInputRange().Y() - fBinSize / 2.) + binRight = binLeft; + + /* + std::cout << "Top : " << GetInputRange().Y() - fBinSize/2. << std::endl; + std::cout << "binLeft : " << binLeft << std::endl; + std::cout << "binRight : " << binRight << std::endl; + std::cout << "dLeft : " << distLeft << std::endl; + std::cout << "dLeft/fBinSize : " << distLeft/fBinSize << std::endl; + std::cout << "1 - distLeft/fBinSize : " << 1 - distLeft/fBinSize << std::endl; + */ + + for (std::size_t n = 0; n < fResponseMatrix[binLeft].size(); n++) { + Double_t output = fOrigin.Y() + ((double)n + 0.5) * fBinSize; + + Double_t value = fResponseMatrix[binLeft][n] * (1 - distLeft / fBinSize) + + fResponseMatrix[binRight][n] * distLeft / fBinSize; + + std::pair outp{output, value}; + + response.push_back(outp); + + /* + std::cout << "n: " << n << " output : " << output << std::endl; + std::cout << "response: " << response << std::endl; + */ + } + + return response; +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux +/// +void TRestResponse::PrintResponseMatrix(Int_t fromRow = 0, Int_t toRow = 0) { + TRestTools::PrintTable(fResponseMatrix, fromRow, toRow); +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux +/// +void TRestResponse::PrintMetadata() { + TRestMetadata::PrintMetadata(); + + RESTMetadata << "Response file : " << fFilename << RESTendl; + RESTMetadata << "Variable : " << fVariable << RESTendl; + RESTMetadata << "Bin size : " << fBinSize << RESTendl; + RESTMetadata << " " << RESTendl; + + if (!fResponseMatrix.empty()) { + RESTMetadata << "Response matrix has been loaded" << RESTendl; + RESTMetadata << " - Number of columns: " << fResponseMatrix[0].size() << RESTendl; + RESTMetadata << " - Number of rows : " << fResponseMatrix.size() << RESTendl; + RESTMetadata << " - Input range : " << GetInputRange().X() << " - " << GetInputRange().Y() + << RESTendl; + RESTMetadata << " - Output range : " << GetOutputRange().X() << " - " << GetOutputRange().Y() + << RESTendl; + + if (fTransposed) { + RESTMetadata << " " << RESTendl; + RESTMetadata << "Original matrix was transposed" << RESTendl; + } + } else { + RESTMetadata << "Response matrix has NOT been loaded" << RESTendl; + RESTMetadata << "Try calling TRestResponse::LoadResponse()" << RESTendl; + } + if (fInterpolation) { + RESTMetadata << " " << RESTendl; + RESTMetadata << "Interpolation is enabled" << RESTendl; + } else { + RESTMetadata << " " << RESTendl; + RESTMetadata << "Interpolation is disabled" << RESTendl; + } + RESTMetadata << "----" << RESTendl; +} diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx new file mode 100644 index 000000000..12fe749d6 --- /dev/null +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -0,0 +1,223 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +///////////////////////////////////////////////////////////////////////// +/// Documentation TOBE written +/// +/// +///---------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2022-December: First implementation of TRestSensitivity +/// Javier Galan +/// +/// \class TRestSensitivity +/// \author: Javier Galan (javier.galan.lacarra@cern.ch) +/// +///
+/// +#include +#include + +ClassImp(TRestSensitivity); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestSensitivity::TRestSensitivity() { Initialize(); } + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestSensitivity::~TRestSensitivity() {} + +///////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// \param name The name of the specific metadata. It will be used to find the +/// corresponding TRestAxionMagneticField section inside the RML. +/// +TRestSensitivity::TRestSensitivity(const char* cfgFileName, const std::string& name) + : TRestMetadata(cfgFileName) { + LoadConfigFromFile(fConfigFileName, name); +} + +/////////////////////////////////////////////// +/// \brief It will initialize the data frame with the filelist and column names +/// (or observables) that have been defined by the user. +/// +void TRestSensitivity::Initialize() { SetSectionName(this->ClassName()); } + +/////////////////////////////////////////////// +/// \brief It will return a value of the coupling, g4, such that (chi-chi0) gets +/// closer to the target value given by argument. The factor will be used to +/// increase or decrease the coupling, and evaluate the likelihood. +/// +Double_t TRestSensitivity::ApproachByFactor(Double_t g4, Double_t chi0, Double_t target, Double_t factor) { + if (factor == 1) { + return 0; + } + + /// Coarse movement to get to Chi2 above target + Double_t Chi2 = 0; + do { + Chi2 = 0; + for (const auto& exp : fExperiments) Chi2 += -2 * UnbinnedLogLikelihood(exp, g4); + + g4 = factor * g4; + } while (Chi2 - chi0 < target); + g4 = g4 / factor; + + /// Coarse movement to get to Chi2 below target (/2) + do { + Chi2 = 0; + for (const auto& exp : fExperiments) Chi2 += -2 * UnbinnedLogLikelihood(exp, g4); + + g4 = g4 / factor; + } while (Chi2 - chi0 > target); + + return g4 * factor; +} + +/////////////////////////////////////////////// +/// \brief It will return the coupling value for which Chi=sigma +/// +Double_t TRestSensitivity::GetCoupling(Double_t sigma, Double_t precision) { + Double_t Chi2_0 = 0; + for (const auto& exp : fExperiments) Chi2_0 += -2 * UnbinnedLogLikelihood(exp, 0); + + Double_t target = sigma * sigma; + + Double_t g4 = 0.5; + + g4 = ApproachByFactor(g4, Chi2_0, target, 2); + g4 = ApproachByFactor(g4, Chi2_0, target, 1.2); + g4 = ApproachByFactor(g4, Chi2_0, target, 1.02); + g4 = ApproachByFactor(g4, Chi2_0, target, 1.0002); + + return g4; +} + +/////////////////////////////////////////////// +/// \brief It returns the Log(L) for the experiment and coupling given by argument. +/// +Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experiment, Double_t g4) { + Double_t lhood = 0; + if (!experiment->IsDataReady()) { + RESTError << "TRestSensitivity::UnbinnedLogLikelihood. Experiment " << experiment->GetName() + << " is not ready!" << RESTendl; + return lhood; + } + + Double_t signal = g4 * experiment->GetSignal()->GetTotalRate() * experiment->GetExposureInSeconds(); + + lhood = -signal; + + if (ROOT::IsImplicitMTEnabled()) ROOT::DisableImplicitMT(); + + std::vector> trackingData; + for (const auto& var : experiment->GetSignal()->GetVariables()) { + auto values = experiment->GetExperimentalDataFrame().Take(var); + std::vector vDbl = std::move(*values); + trackingData.push_back(vDbl); + } + + for (size_t n = 0; n < trackingData[0].size(); n++) { + std::vector point; + for (size_t m = 0; m < trackingData.size(); m++) point.push_back(trackingData[m][n]); + + Double_t bckRate = experiment->GetBackground()->GetRate(point); + Double_t sgnlRate = experiment->GetSignal()->GetRate(point); + + Double_t expectedRate = bckRate + g4 * sgnlRate; + lhood += TMath::Log(expectedRate); + } + + return lhood; +} + +/////////////////////////////////////////////// +/// \brief +/// +TH1D* TRestSensitivity::SignalStatisticalTest(Int_t N) { + std::vector couplings; + for (int n = 0; n < N; n++) { + for (const auto& exp : fExperiments) exp->GetSignal()->RegenerateActiveNodeDensity(); + + Double_t coupling = TMath::Sqrt(TMath::Sqrt(GetCoupling())); + couplings.push_back(coupling); + } + + // Directly assign the minimum and maximum values + double min_value = *std::min_element(couplings.begin(), couplings.end()); + double max_value = *std::max_element(couplings.begin(), couplings.end()); + + if (fSignalTest) delete fSignalTest; + fSignalTest = new TH1D("SignalTest", "A signal test", 100, 0.9 * min_value, 1.1 * max_value); + for (const auto& coup : couplings) fSignalTest->Fill(coup); + + return fSignalTest; +} + +///////////////////////////////////////////// +/// \brief It customizes the retrieval of XML data values of this class +/// +void TRestSensitivity::InitFromConfigFile() { + TRestMetadata::InitFromConfigFile(); + + int cont = 0; + TRestMetadata* metadata = (TRestMetadata*)this->InstantiateChildMetadata(cont, "Experiment"); + while (metadata != nullptr) { + cont++; + if (metadata->InheritsFrom("TRestExperimentList")) { + TRestExperimentList* experimentsList = (TRestExperimentList*)metadata; + std::vector exList = experimentsList->GetExperiments(); + fExperiments.insert(fExperiments.end(), exList.begin(), exList.end()); + } else if (metadata->InheritsFrom("TRestExperiment")) { + fExperiments.push_back((TRestExperiment*)metadata); + } + + metadata = (TRestMetadata*)this->InstantiateChildMetadata(cont, "Experiment"); + } + + Initialize(); +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux +/// +void TRestSensitivity::PrintMetadata() { + TRestMetadata::PrintMetadata(); + + RESTMetadata << "----" << RESTendl; +} diff --git a/source/framework/tools/inc/TRestPhysics.h b/source/framework/tools/inc/TRestPhysics.h index dca1f0637..e48018eda 100644 --- a/source/framework/tools/inc/TRestPhysics.h +++ b/source/framework/tools/inc/TRestPhysics.h @@ -49,7 +49,7 @@ constexpr double kBoltzman = 1.380E-23; constexpr double hPlanck = 1.054E-34; /// A meter in eV -constexpr double PhMeterIneV = 5067731.236453719; // 8.0655447654281218E5;// 506.773123645372; +constexpr double PhMeterIneV = 5067731.236453719; /// A second in eV (using natural units) constexpr double secondIneV = 1519225802531030.2; /// Electron charge in natural units diff --git a/source/framework/tools/inc/TRestStringHelper.h b/source/framework/tools/inc/TRestStringHelper.h index 879e5ee43..c34c8cbcf 100644 --- a/source/framework/tools/inc/TRestStringHelper.h +++ b/source/framework/tools/inc/TRestStringHelper.h @@ -35,6 +35,7 @@ Float_t StringToFloat(std::string in); Double_t StringToDouble(std::string in); Int_t StringToInteger(std::string in); std::string IntegerToString(Int_t n, std::string format = "%d"); +std::vector IntegerToBinary(int number, size_t dimension = 0); std::string DoubleToString(Double_t d, std::string format = "%8.6e"); Bool_t StringToBool(std::string booleanString); Long64_t StringToLong(std::string in); @@ -45,6 +46,7 @@ std::vector Split(std::string in, std::string separator, bool allow std::vector StringToElements(std::string in, std::string separator); std::vector StringToElements(std::string in, std::string headChar, std::string separator, std::string tailChar); +std::string RemoveDelimiters(std::string in); std::string RemoveWhiteSpaces(std::string in); std::string Replace(std::string in, std::string thisString, std::string byThisString, size_t fromPosition = 0, Int_t N = 0); diff --git a/source/framework/tools/inc/TRestTools.h b/source/framework/tools/inc/TRestTools.h index 58836b4bf..34e1e18ab 100644 --- a/source/framework/tools/inc/TRestTools.h +++ b/source/framework/tools/inc/TRestTools.h @@ -63,6 +63,8 @@ class TRestTools { Int_t skipLines = 0, std::string separator = "\t"); static int ReadASCIITable(std::string fName, std::vector>& data, Int_t skipLines = 0, std::string separator = "\t"); + static int ReadASCIITable(std::string fName, std::vector>& data, + Int_t skipLines = 0, std::string separator = "\t"); static int ReadCSVFile(std::string fName, std::vector>& data, Int_t skipLines = 0); static int ReadCSVFile(std::string fName, std::vector>& data, Int_t skipLines = 0); @@ -119,7 +121,7 @@ class TRestTools { static std::string GetPureFileName(const std::string& fullPathFileName); static std::string SearchFileInPath(std::vector path, std::string filename); static bool CheckFileIsAccessible(const std::string&); - static std::vector GetFilesMatchingPattern(std::string pattern); + static std::vector GetFilesMatchingPattern(std::string pattern, bool unlimited = false); static int ConvertVersionCode(std::string in); static std::istream& GetLine(std::istream& is, std::string& t); @@ -131,6 +133,8 @@ class TRestTools { static std::string POSTRequest(const std::string& url, const std::map& keys); static void ChangeDirectory(const std::string& toDirectory); + + static std::vector CanvasDivisions(int n); }; namespace REST_InitTools { diff --git a/source/framework/tools/src/TRestStringHelper.cxx b/source/framework/tools/src/TRestStringHelper.cxx index 0d96327b6..101e56a99 100644 --- a/source/framework/tools/src/TRestStringHelper.cxx +++ b/source/framework/tools/src/TRestStringHelper.cxx @@ -256,11 +256,13 @@ vector REST_StringHelper::Split(string in, string separator, bool allowB /////////////////////////////////////////////// /// \brief Convert the input string into a vector of double elements /// +/// The method will remove any delimiters found in the string (), [] or {}. +/// /// e.g. Input: "1,2,3,4", Output: {1.,2.,3.,4.} /// vector REST_StringHelper::StringToElements(string in, string separator) { vector result; - vector vec_str = REST_StringHelper::Split(in, separator); + vector vec_str = REST_StringHelper::Split(RemoveDelimiters(in), separator); for (unsigned int i = 0; i < vec_str.size(); i++) { double temp = REST_StringHelper::StringToDouble(vec_str[i]); result.push_back(temp); @@ -294,6 +296,19 @@ vector REST_StringHelper::StringToElements(string in, string headChar, s return result; } +/////////////////////////////////////////////// +/// \brief Returns the input string removing any delimiters ({[]}) +/// +string REST_StringHelper::RemoveDelimiters(string in) { + string out = in; + size_t pos = out.find_first_of("+-*/e^%"); + while ((pos = out.find_first_of("({[]})")) != string::npos) { + out.erase(pos, 1); + } + + return out; +} + /////////////////////////////////////////////// /// \brief Returns the input string removing all white spaces. /// @@ -464,7 +479,7 @@ Int_t REST_StringHelper::DiffString(const string& source, const string& target) } /////////////////////////////////////////////// -/// \brief Replace every occurences of **thisSring** by **byThisString** inside +/// \brief Replace any occurences of **thisSring** by **byThisString** inside /// string **in**. /// string REST_StringHelper::Replace(string in, string thisString, string byThisString, size_t fromPosition, @@ -618,6 +633,37 @@ Int_t REST_StringHelper::StringToInteger(string in) { /// string REST_StringHelper::IntegerToString(Int_t n, string format) { return Form(format.c_str(), n); } +/////////////////////////////////////////////// +/// \brief It returns an integer vector with the binary digits decomposed. +/// +/// Example: IntegerToBinary(7) will return { 1, 1, 1 }. +/// +/// Optionally we can fix the minimum number of digits to be returned, so that +/// it will be filled with zeros to the left. +/// +/// Example: IntegerToBinary(9,8) will return { 0, 0, 0, 0, 1, 0, 0, 1 }. +/// +std::vector REST_StringHelper::IntegerToBinary(int number, size_t dimension) { + std::vector binaryNumber; + + if (number == 0) { + binaryNumber.insert(binaryNumber.begin(), dimension, 0); + if (binaryNumber.empty()) binaryNumber.push_back(0); + return binaryNumber; + } + + while (number > 0) { + int digit = number % 2; + binaryNumber.insert(binaryNumber.begin(), digit); + number /= 2; + } + + if (dimension > binaryNumber.size()) + binaryNumber.insert(binaryNumber.begin(), dimension - binaryNumber.size(), 0); + + return binaryNumber; +} + /////////////////////////////////////////////// /// \brief Gets a string from a double /// diff --git a/source/framework/tools/src/TRestTools.cxx b/source/framework/tools/src/TRestTools.cxx index 0af3237fa..9502be51a 100644 --- a/source/framework/tools/src/TRestTools.cxx +++ b/source/framework/tools/src/TRestTools.cxx @@ -174,6 +174,8 @@ template int TRestTools::PrintTable(std::vector> data, template int TRestTools::PrintTable(std::vector> data, Int_t start, Int_t end); template int TRestTools::PrintTable(std::vector> data, Int_t start, Int_t end); +template int TRestTools::PrintTable(std::vector> data, Int_t start, + Int_t end); /////////////////////////////////////////////// /// \brief Writes the contents of the vector table given as argument to `fname`. @@ -507,6 +509,58 @@ template std::vector TRestTools::GetColumnFromTable( template std::vector TRestTools::GetColumnFromTable( const std::vector>& data, unsigned int column); +template std::vector TRestTools::GetColumnFromTable( + const std::vector>& data, unsigned int column); + +/////////////////////////////////////////////// +/// \brief Reads an ASCII file containing a table with values +/// +/// This method will open the file fName. This file should contain a tabulated +/// ASCII table containing any format values. The values on the table will be +/// loaded in the matrix provided through the argument `data`. The content of +/// `data` will be cleared in this method. +/// +/// If any header in the file is present, it should be skipped using the argument +///`skipLines` or preceding any line inside the header using `#`. +/// +/// This table will just split the ASCII elements inside a std::string matrix +/// +int TRestTools::ReadASCIITable(string fName, std::vector>& data, Int_t skipLines, + std::string separator) { + if (!TRestTools::isValidFile((string)fName)) { + cout << "TRestTools::ReadASCIITable. Error" << endl; + cout << "Cannot open file : " << fName << endl; + return 0; + } + + data.clear(); + + std::ifstream fin(fName); + + // First we create a table with string values + std::vector> values; + + for (string line; std::getline(fin, line);) { + if (skipLines > 0) { + skipLines--; + continue; + } + + if (line.find("#") == string::npos) { + std::istringstream in(line); + + std::string token; + std::vector tokens; + while (std::getline(in, token, (char)separator[0])) { + tokens.push_back(token); + } + data.push_back(tokens); + } + } + + return 1; +} + /////////////////////////////////////////////// /// \brief Reads an ASCII file containing a table with values /// @@ -915,7 +969,11 @@ bool TRestTools::CheckFileIsAccessible(const std::string& filename) { /// \brief Returns a list of files whose name match the pattern string. Key word /// is "*". e.g. abc00*.root /// -vector TRestTools::GetFilesMatchingPattern(string pattern) { +/// Argument unlimited will fix an issue with the number of files being to high. +/// However, it causes issues when searching/listing the macros. +/// The default value for unlimited is `false`. +/// +vector TRestTools::GetFilesMatchingPattern(string pattern, bool unlimited) { std::vector outputFileNames; if (pattern != "") { vector items = Split(pattern, "\n"); @@ -944,11 +1002,25 @@ vector TRestTools::GetFilesMatchingPattern(string pattern) { } } #else - string a = Execute("find " + item); - auto b = Split(a, "\n"); + auto path_name = SeparatePathAndName(item); + if (unlimited) { + std::string currentDir = filesystem::current_path(); + ChangeDirectory(path_name.first); + string a = Execute("find -type f -name \'" + path_name.second + "\'"); + ChangeDirectory(currentDir); + auto b = Split(a, "\n"); + + for (unsigned int i = 0; i < b.size(); i++) { + outputFileNames.push_back(path_name.first + "/" + b[i]); + } + + } else { + string a = Execute("find " + item); + auto b = Split(a, "\n"); - for (unsigned int i = 0; i < b.size(); i++) { - outputFileNames.push_back(b[i]); + for (unsigned int i = 0; i < b.size(); i++) { + outputFileNames.push_back(b[i]); + } } #endif @@ -1249,6 +1321,58 @@ int TRestTools::UploadToServer(string localFile, string remoteFile, string metho void TRestTools::ChangeDirectory(const string& toDirectory) { filesystem::current_path(toDirectory); } +/////////////////////////////////////////////// +/// \brief It returns a vector with 2 components {a,b}, the components satisfy that `a x b = n`, +/// being the ratio a/b as close to 1 as possible. +/// +/// This method can be used to help dividing a canvas that will contain a number `n` of plots. +/// +/// If `n` is a prime number, then the pair generated will be `n x 1`. +/// +std::vector TRestTools::CanvasDivisions(int n) { + std::vector r; + for (int i = 2; i * i <= n; i += 1 + (i > 2)) { + while ((n % i) == 0) { + r.push_back(i); + n /= i; + } + } + if (n != 1) r.push_back(n); + + while (r.size() > 2) { + // We multiply the 2 lowest elements and + // replace the elements in the vector by the result + auto min1 = std::min_element(r.begin(), r.end()); + int low1 = *min1; + + // Remove the first element equal to min1 (efficient way) + auto it = std::find(r.begin(), r.end(), low1); + if (it != r.end()) { + std::iter_swap(it, r.end() - 1); + r.erase(r.end() - 1); + } + + auto min2 = std::min_element(r.begin(), r.end()); + int low2 = *min2; + + // Remove the first element equal to min2 (efficient way) + it = std::find(r.begin(), r.end(), low2); + if (it != r.end()) { + std::iter_swap(it, r.end() - 1); + r.erase(r.end() - 1); + } + + int resultado = low1 * low2; + r.push_back(resultado); + } + + std::sort(r.begin(), r.end()); + + if (r.size() == 1) r.push_back(1); + + return r; +} + string ValueWithQuantity::ToString() const { string unit; auto value = fValue; diff --git a/source/libraries/axion b/source/libraries/axion index 21375385f..3800d90f4 160000 --- a/source/libraries/axion +++ b/source/libraries/axion @@ -1 +1 @@ -Subproject commit 21375385fc4cad92bda60ba29f641b9afad9339e +Subproject commit 3800d90f41622ec0c465f7b4d4c5f384e4d9433a diff --git a/source/libraries/legacy b/source/libraries/legacy index 831424307..caa4e2540 160000 --- a/source/libraries/legacy +++ b/source/libraries/legacy @@ -1 +1 @@ -Subproject commit 83142430760cb35c2b3c597b24c1d884a01c0fff +Subproject commit caa4e2540f87b33071b1ae8e21f277aa69e60861