diff --git a/.github/workflows/validation.yml b/.github/workflows/validation.yml index ff9f80245..4070ac1f5 100644 --- a/.github/workflows/validation.yml +++ b/.github/workflows/validation.yml @@ -462,7 +462,7 @@ jobs: source ${{ env.REST_PATH }}/thisREST.sh cd ${{ env.REST_PATH }}/examples/restG4/07.FullChainDecay/ restG4 fullChain.rml -o Run00001_U238_FullChainDecay.root - restRoot -b -q Validate.C'("Run00001_U238_FullChainDecay.root", 15)' + restRoot -b -q Validate.C'("Run00001_U238_FullChainDecay.root", 16)' restG4 singleDecay.rml -o Run00002_U238_SingleChainDecay.root restRoot -b -q Validate.C'("Run00002_U238_SingleChainDecay.root", 1)' export REST_ISOTOPE=Be7 diff --git a/CMakeLists.txt b/CMakeLists.txt index f58e2e85e..557218d49 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -66,6 +66,11 @@ endif () set(external_include_dirs) set(external_libs) +if (CMAKE_SYSTEM_NAME MATCHES "Darwin") + set(external_include_dirs ${external_include_dirs} "/usr/local/include/") + message(STATUS "Adding BREW include directory: /usr/local/include") +endif (CMAKE_SYSTEM_NAME MATCHES "Darwin") + set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake ${DECAY_PATH}/cmake ${CMAKE_MODULE_PATH}) set(CMAKE_MACOSX_RPATH 1) @@ -356,7 +361,7 @@ if (CMAKE_SYSTEM_NAME MATCHES "Darwin") # we must call library install here in endif () # Copy pcm files -if (CMAKE_SYSTEM_NAME MATCHES "Darwin" OR CMAKE_SYSTEM_NAME MATCHES "Windows") +if (CMAKE_SYSTEM_NAME MATCHES "Windows") # Copy pcm files to bin install( diff --git a/README.md b/README.md index 42d99934d..8183f1c82 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # The REST Framework [![DOI](https://zenodo.org/badge/324291710.svg)](http://doi.org/10.5281/zenodo.4528985) -[![pipeline status](https://gitlab.cern.ch/rest-for-physics/framework/badges/master/pipeline.svg)](https://gitlab.cern.ch/rest-for-physics/framework/-/commits/master) +[![pipeline status](https://github.com/rest-for-physics/framework/actions/workflows/validation.yml/badge.svg)](https://github.com/rest-for-physics/framework) [![website](https://img.shields.io/badge/user-guide-E8B6FF.svg)](https://rest-for-physics.github.io) [![api](https://img.shields.io/badge/user-API-FFCA78.svg)](https://rest-for-physics.github.io/framework/) [![forum](https://img.shields.io/badge/user-forum-AAFF90.svg)](https://rest-forum.unizar.es/) @@ -77,6 +77,8 @@ A major change at 2.3 will prevent from backwards compatibility, since class nam - PandaX-III: Searching for neutrinoless double beta decay with high pressure 136Xe gas time projection chambers. [X. Chen et al., Science China Physics, Mechanics & Astronomy 60, 061011 (2017)](https://doi.org/10.1007/s11433-017-9028-0) [arXiv:1610.08883](https://arxiv.org/abs/1610.08883). ## Presentations +- Background model studies for IAXO using the REST-for-Physics framework, Luis Obis, [2024-Apr, VIEnna Workshop on Simulations 2024 (VIEWS24), Vienna](https://indico.cern.ch/event/1275551/contributions/5858816/). +- REST-for-Physics framework for Geant4 simulations and data analysis, Álvaro Ezquerro, [2024-Apr, VIEnna Workshop on Simulations 2024 (VIEWS24), Vienna](https://indico.cern.ch/event/1275551/contributions/5858814/). - REST-for-Physics, Luis Obis, [2022-May, ROOT Users Workshop, FermiLab](https://indico.fnal.gov/event/23628/contributions/240755/). - REST v2.0 : A data analysis and simulation framework for micro-patterned readout detectors., Javier Galan, [2016-Dec, 8th Symposium on Large TPCs for low-energy rare event detection, Paris](https://indico.cern.ch/event/473362/contributions/2334838/). diff --git a/cmake/cmake_uninstall.cmake.in b/cmake/cmake_uninstall.cmake.in index e925304f8..dc36259c0 100644 --- a/cmake/cmake_uninstall.cmake.in +++ b/cmake/cmake_uninstall.cmake.in @@ -8,7 +8,7 @@ string(REGEX REPLACE "\n" ";" files "${files}") foreach(file ${files}) message(STATUS "Uninstalling $ENV{DESTDIR}${file}") if(IS_SYMLINK "$ENV{DESTDIR}${file}" OR EXISTS "$ENV{DESTDIR}${file}") - exec_program( + execute_process( "${CMAKE_COMMAND}" ARGS "-E remove \"$ENV{DESTDIR}${file}\"" OUTPUT_VARIABLE rm_out RETURN_VALUE rm_retval diff --git a/data/distributions/CosmicsCry.root b/data/distributions/CosmicsCry.root new file mode 100644 index 000000000..7f6dace14 Binary files /dev/null and b/data/distributions/CosmicsCry.root differ diff --git a/doc/doxygen/images/radialstrippedmask.png b/doc/doxygen/images/radialstrippedmask.png new file mode 100644 index 000000000..b4c928a32 Binary files /dev/null and b/doc/doxygen/images/radialstrippedmask.png differ diff --git a/examples/masks.rml b/examples/masks.rml index 15e8a9c2b..da33120f0 100644 --- a/examples/masks.rml +++ b/examples/masks.rml @@ -14,6 +14,17 @@ + + + + + + + + + + + diff --git a/macros/REST_ListMacros.C b/macros/REST_ListMacros.C index 3059259af..4333dead8 100644 --- a/macros/REST_ListMacros.C +++ b/macros/REST_ListMacros.C @@ -13,40 +13,21 @@ using namespace std; //*** Lists all the official macros together with its documentation //*** //******************************************************************************************************* -Int_t REST_ListMacros() { +Int_t REST_ListMacros(int details = 0, std::string onlyThisMacro = "") { string macrosPath = (string)getenv("REST_PATH") + "/macros"; vector directories = TRestTools::GetSubdirectories(macrosPath); - cout << "Directory : " << macrosPath << endl; + cout << endl; + cout << "Entering directory : " << macrosPath << endl; vector main_files = TRestTools::GetFilesMatchingPattern(macrosPath + "/REST_*.C"); - for (int m = 0; m < main_files.size(); m++) { - cout << " ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl; - cout << " ++ " - << " Macro : " << TRestTools::SeparatePathAndName(main_files[m]).second << endl; - std::ifstream t(main_files[m]); - std::string str((std::istreambuf_iterator(t)), std::istreambuf_iterator()); - - std::vector lines = REST_StringHelper::Split(str, "\n"); - - cout << " ++ " << endl; - for (int l = 0; l < lines.size(); l++) - if (lines[l].find("//*** ") != string::npos) - cout << " ++ " << lines[l].substr(6, -1) << endl; - cout << " ++ " << endl; - } - cout << " ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl; - - for (int n = 0; n < directories.size(); n++) { - cout << "Directory : " << directories[n] << endl; - if (directories[n].find("pipeline") != string::npos) continue; - if (directories[n].find("/macros/mac/") != string::npos) continue; - vector files = TRestTools::GetFilesMatchingPattern(directories[n] + "/REST_*.C"); - - for (int m = 0; m < files.size(); m++) { + if (details) { + for (int m = 0; m < main_files.size(); m++) { + std::string macro = TRestTools::SeparatePathAndName(main_files[m]).second; + if (!onlyThisMacro.empty() && onlyThisMacro != macro) continue; cout << " ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl; cout << " ++ " - << " Macro : " << TRestTools::SeparatePathAndName(files[m]).second << endl; - std::ifstream t(files[m]); + << " Macro : " << macro << endl; + std::ifstream t(main_files[m]); std::string str((std::istreambuf_iterator(t)), std::istreambuf_iterator()); std::vector lines = REST_StringHelper::Split(str, "\n"); @@ -58,6 +39,53 @@ Int_t REST_ListMacros() { cout << " ++ " << endl; } cout << " ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl; + } else { + for (int m = 0; m < main_files.size(); m++) { + cout << " + " << TRestTools::SeparatePathAndName(main_files[m]).second << endl; + } + } + + for (int n = 0; n < directories.size(); n++) { + cout << endl; + cout << "Entering directory : " << directories[n] << endl; + if (directories[n].find("pipeline") != string::npos) continue; + if (directories[n].find("/macros/mac/") != string::npos) continue; + vector files = TRestTools::GetFilesMatchingPattern(directories[n] + "/REST_*.C"); + + if (details) { + for (int m = 0; m < files.size(); m++) { + std::string macro = TRestTools::SeparatePathAndName(files[m]).second; + if (!onlyThisMacro.empty() && onlyThisMacro != macro) continue; + cout << " ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" + << endl; + cout << " ++ Macro : " << macro << endl; + std::ifstream t(files[m]); + std::string str((std::istreambuf_iterator(t)), std::istreambuf_iterator()); + + std::vector lines = REST_StringHelper::Split(str, "\n"); + + cout << " ++ " << endl; + for (int l = 0; l < lines.size(); l++) + if (lines[l].find("//*** ") != string::npos) + cout << " ++ " << lines[l].substr(6, -1) << endl; + cout << " ++ " << endl; + } + cout << " ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl; + } else { + for (int m = 0; m < files.size(); m++) { + cout << " + " << TRestTools::SeparatePathAndName(files[m]).second << endl; + } + std::cout << std::endl; + std::cout << " ------- " << std::endl; + std::cout << "IMPORTANT. To get a more detailed macro documentation enable the details argument:" + << std::endl; + std::cout << "Execute: restListMacros 1" << std::endl; + std::cout << "OR if you want to show only the documentation of a given macro, add the macro.C " + "filename as argument" + << std::endl; + std::cout << "Execute: restListMacros 1 REST_CheckValidRuns.C" << std::endl; + std::cout << " ------- " << std::endl; + } } return 0; } diff --git a/macros/REST_OpenInputFile.C b/macros/REST_OpenInputFile.C index 7812d8f5f..0f1af5cd6 100644 --- a/macros/REST_OpenInputFile.C +++ b/macros/REST_OpenInputFile.C @@ -47,7 +47,7 @@ void REST_OpenInputFile(const std::string& fileName) { 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()"); + printf("%s\n\n", " - dSet->GetDataFrame().Display({\"colName1\", \"colName2\"})->Print()"); if (dSet) delete dSet; dSet = new TRestDataSet(); dSet->EnableMultiThreading(false); diff --git a/scripts/generateVersionHeader.py b/scripts/generateVersionHeader.py index 29cc21712..1d85a9b38 100755 --- a/scripts/generateVersionHeader.py +++ b/scripts/generateVersionHeader.py @@ -159,29 +159,29 @@ ) print("\n") print( - "Once your PR has been accepted and merged, you should generate a new Git tag at the master branch.\n" -) -print( - "-----> git tag -a v" - + str(a) - + "." - + str(b) - + "." - + str(c) - + ' -m "Update to version v' - + str(a) - + "." - + str(b) - + "." - + str(c) - + '"\n' -) -print( - "And push the changes to repository. You should also push your branch to GitHub if you have not already.\n" -) -print("-----> git push origin v" + str(a) + "." + str(b) + "." + str(c) + "\n") -print( - "IMPORTANT. Summarize the changes in the tag generated at the Gitlab repository website.\n" + "Once your PR has been accepted and merged, you should generate a new Git tag and RELEASE at the master branch.\n" ) +# print( +# "-----> git tag -a v" +# + str(a) +# + "." +# + str(b) +# + "." +# + str(c) +# + ' -m "Update to version v' +# + str(a) +# + "." +# + str(b) +# + "." +# + str(c) +# + '"\n' +# ) +# print( +# "And push the changes to repository. You should also push your branch to GitHub if you have not already.\n" +# ) +# print("-----> git push origin v" + str(a) + "." + str(b) + "." + str(c) + "\n") +# print( +# "IMPORTANT. Summarize the changes in the tag generated at the Gitlab repository website.\n" +# ) exit(0) diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h index 50000a01b..617c5f0aa 100644 --- a/source/framework/analysis/inc/TRestDataSetGainMap.h +++ b/source/framework/analysis/inc/TRestDataSetGainMap.h @@ -33,6 +33,7 @@ #include #include +#include "TRestCut.h" #include "TRestDataSet.h" #include "TRestMetadata.h" @@ -42,13 +43,26 @@ class TRestDataSetGainMap : public TRestMetadata { class Module; private: - std::string fObservable = ""; //"rawAna_ThresholdIntegral"; //< - std::string fSpatialObservableX = ""; //"hitsAna_xMean"; //< - std::string fSpatialObservableY = ""; //"hitsAna_yMean"; //< + /// Observable that will be used to calculate the gain map + std::string fObservable = ""; //< - std::vector fModulesCal = {}; - std::string fCalibFileName = ""; - std::string fOutputFileName = ""; + /// Observable that will be used to segmentize the gain map in the x direction + std::string fSpatialObservableX = ""; //< + + /// Observable that will be used to segmentize the gain map in the y direction + std::string fSpatialObservableY = ""; //< + + /// List of modules + std::vector fModulesCal = {}; //< + + /// Name of the file that contains the calibration data + std::string fCalibFileName = ""; //< + + /// Name of the file where the gain map was (or will be) exported + std::string fOutputFileName = ""; //< + + /// Cut to be applied to the calibration data + TRestCut* fCut = nullptr; //< void Initialize() override; void InitFromConfigFile() override; @@ -69,14 +83,18 @@ class TRestDataSetGainMap : public TRestMetadata { std::string GetObservable() const { return fObservable; } std::string GetSpatialObservableX() const { return fSpatialObservableX; } std::string GetSpatialObservableY() const { return fSpatialObservableY; } + TRestCut* GetCut() const { return fCut; } + Module* GetModule(const size_t index = 0); Module* GetModule(const int planeID, const int moduleID); double GetSlopeParameter(const int planeID, const int moduleID, const double x, const double y); double GetInterceptParameter(const int planeID, const int moduleID, const double x, const double y); + double GetSlopeParameterFullSpc(const int planeID, const int moduleID); + double GetInterceptParameterFullSpc(const int planeID, const int moduleID); void SetCalibrationFileName(const std::string& fileName) { fCalibFileName = fileName; } void SetOutputFileName(const std::string& fileName) { fOutputFileName = fileName; } - void SetModuleCalibration(const Module& moduleCal); + void SetModule(const Module& moduleCal); void SetObservable(const std::string& observable) { fObservable = observable; } void SetSpatialObservableX(const std::string& spatialObservableX) { fSpatialObservableX = spatialObservableX; @@ -84,52 +102,106 @@ class TRestDataSetGainMap : public TRestMetadata { void SetSpatialObservableY(const std::string& spatialObservableY) { fSpatialObservableY = spatialObservableY; } + void SetCut(TRestCut* cut) { fCut = cut; } void Import(const std::string& fileName); void Export(const std::string& fileName = ""); TRestDataSetGainMap& operator=(TRestDataSetGainMap& src); - public: void PrintMetadata() override; void GenerateGainMap(); - void CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName = ""); + void CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName = "", + std::vector excludeColumns = {}); TRestDataSetGainMap(); TRestDataSetGainMap(const char* configFilename, std::string name = ""); ~TRestDataSetGainMap(); - ClassDefOverride(TRestDataSetGainMap, 1); + ClassDefOverride(TRestDataSetGainMap, 2); class Module { private: - const TRestDataSetGainMap* p = nullptr; // fEnergyPeaks = {}; - std::vector fRangePeaks = {}; //{TVector2(230000, 650000), TVector2(40000, 230000)}; - TVector2 fCalibRange = TVector2(0, 0); //< // Calibration range - Int_t fNBins = 100; //< // Number of bins for the spectrum histograms - std::string fDefinitionCut = ""; //"TREXsides_tagId == 2"; //< - - Int_t fNumberOfSegmentsX = 1; //< - Int_t fNumberOfSegmentsY = 1; //< - TVector2 fReadoutRange = TVector2(-1, 246.24); //< // Readout dimensions - std::set fSplitX = {}; //< - std::set fSplitY = {}; //< - - std::string fDataSetFileName = ""; //< // File name for the calibration dataset - - std::vector> fSlope = {}; //< + /// Pointer to the parent class + const TRestDataSetGainMap* p = nullptr; // FitPeaks(TH1F* hSeg, TGraph* gr); + std::pair UpdateCalibrationFits(TH1F* hSeg, TGraph* gr); + + public: + /// Plane ID (unique identifier). Although it is not linked to any TRestDetectorReadout it is + /// recommended to use the same. + Int_t fPlaneId = -1; //< + + // Module ID (unique identifier). Although it is not linked to any TRestDetectorReadout it is + // recommended to use the same. + Int_t fModuleId = -1; //< + + /// Energy of the peaks to be used for the calibration. + std::vector fEnergyPeaks = {}; //< + + /// Range of the peaks to be used for the calibration. If empty it will be automatically calculated. + std::vector fRangePeaks = {}; //< + + /// Calibration range. If fCalibRange.X()>=fCalibRange.Y() the range will be automatically calculated. + TVector2 fCalibRange = TVector2(0, 0); //< + + /// Number of bins for the spectrum histograms. + Int_t fNBins = 100; //< + + /// Cut that defines which events are from this module. + std::string fDefinitionCut = ""; //< + + /// Number of segments in the x direction. + Int_t fNumberOfSegmentsX = 1; //< + + /// Number of segments in the y direction. + Int_t fNumberOfSegmentsY = 1; //< + + /// Readout dimensions + TVector2 fReadoutRange = TVector2(-1, 246.24); //< + + /// Split points in the x direction. + std::set fSplitX = {}; //< + + /// Split points in the y direction. + std::set fSplitY = {}; //< + + /// Name of the file that contains the calibration data. If empty, it will use its parent + /// TRestDataSetGainMap::fCalibFileName. + std::string fDataSetFileName = ""; //< + + /// Array containing the slope of the linear fit for each segment. + std::vector> fSlope = {}; //< + + /// Slope of the calibration linear fit of whole module + double fFullSlope = 0; //< + + /// Intercept of the calibration linear fit of whole module + double fFullIntercept = 0; //< + + /// Array containing the intercept of the linear fit for each segment. std::vector> fIntercept = {}; //< - bool fZeroPoint = false; //< Zero point will be automatically added if there are less than 2 peaks - bool fAutoRangePeaks = true; //< Automatic range peaks - std::vector> fSegSpectra = {}; - std::vector> fSegLinearFit = {}; + /// Add zero point to the calibration linear fit. Zero point will be automatically added if there are + /// less than 2 points. + bool fZeroPoint = false; //< + + /// Automatic range for the peaks fitting. See GenerateGainMap() for more information of the logic. + bool fAutoRangePeaks = true; //< + + /// Array containing the observable spectrum for each segment. + std::vector> fSegSpectra = {}; //< + + /// Array containing the calibration linear fit for each segment. + std::vector> fSegLinearFit = {}; //< + + /// Spectrum of the observable for the whole module. + TH1F* fFullSpectrum = nullptr; //< + + /// Calibration linear fit for the whole module. + TGraph* fFullLinearFit = nullptr; //< public: void AddPeak(const double& energyPeak, const TVector2& rangePeak = TVector2(0, 0)) { @@ -139,9 +211,12 @@ class TRestDataSetGainMap : public TRestMetadata { void LoadConfigFromTiXmlElement(const TiXmlElement* module); + TRestDataSetGainMap* GetParent() const { return const_cast(p); } std::pair GetIndexMatrix(const double x, const double y) const; double GetSlope(const double x, const double y) const; double GetIntercept(const double x, const double y) const; + double GetSlopeFullSpc() const { return fFullSlope; }; + double GetInterceptFullSpc() const { return fFullIntercept; }; Int_t GetPlaneId() const { return fPlaneId; } Int_t GetModuleId() const { return fModuleId; } @@ -155,24 +230,27 @@ class TRestDataSetGainMap : public TRestMetadata { std::set GetSplitX() const { return fSplitX; } std::set GetSplitY() const { return fSplitY; } std::string GetDataSetFileName() const { return fDataSetFileName; } - TVector2 GetReadoutRangeVar() const { return fReadoutRange; } + TVector2 GetReadoutRange() const { return fReadoutRange; } void DrawSpectrum(const bool drawFits = true, const int color = -1, TCanvas* c = nullptr); void DrawSpectrum(const TVector2& position, bool drawFits = true, int color = -1, TCanvas* c = nullptr); void DrawSpectrum(const int index_x, const int index_y, bool drawFits = true, int color = -1, TCanvas* c = nullptr); - void DrawFullSpectrum(); + void DrawFullSpectrum(const bool drawFits = true, const int color = -1, TCanvas* c = nullptr); - void DrawLinearFit(); + void DrawLinearFit(TCanvas* c = nullptr); void DrawLinearFit(const TVector2& position, TCanvas* c = nullptr); void DrawLinearFit(const int index_x, const int index_y, TCanvas* c = nullptr); - void DrawGainMap(const int peakNumber = 0); + void DrawGainMap(const int peakNumber = 0, const bool fullModuleAsRef = true); void Refit(const TVector2& position, const double energy, const TVector2& range); void Refit(const size_t x, const size_t y, const size_t peakNumber, const TVector2& range); + void RefitFullSpc(const double energy, const TVector2& range); + void RefitFullSpc(const size_t peakNumber, const TVector2& range); void UpdateCalibrationFits(const size_t x, const size_t y); + void UpdateCalibrationFitsFullSpc(); void SetPlaneId(const Int_t& planeId) { fPlaneId = planeId; } void SetModuleId(const Int_t& moduleId) { fModuleId = moduleId; } diff --git a/source/framework/analysis/src/TRestDataSetCalibration.cxx b/source/framework/analysis/src/TRestDataSetCalibration.cxx index d9869d29d..f6794f3a9 100644 --- a/source/framework/analysis/src/TRestDataSetCalibration.cxx +++ b/source/framework/analysis/src/TRestDataSetCalibration.cxx @@ -192,7 +192,7 @@ void TRestDataSetCalibration::Calibrate() { if (fCalibFile.empty()) { auto histo = dataSet.GetDataFrame().Histo1D( - {"spectrum", "spectrum", fNBins, fCalibRange.X(), fCalibRange.X()}, fCalObservable); + {"spectrum", "spectrum", fNBins, fCalibRange.X(), fCalibRange.Y()}, fCalObservable); spectrum = std::unique_ptr(static_cast(histo->DrawClone())); // Get position of the maximum diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index 53a5fa0a0..eab96f424 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -186,6 +186,8 @@ void TRestDataSetGainMap::InitFromConfigFile() { moduleDefinition = GetNextElement(moduleDefinition); } + fCut = (TRestCut*)InstantiateChildMetadata("TRestCut"); + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) PrintMetadata(); } @@ -194,30 +196,38 @@ void TRestDataSetGainMap::InitFromConfigFile() { /// void TRestDataSetGainMap::GenerateGainMap() { for (auto& mod : fModulesCal) { - RESTInfo << "Calibrating plane " << mod.GetPlaneId() << " module " << mod.GetModuleId() << RESTendl; + RESTInfo << "Generating gain map of plane " << mod.GetPlaneId() << " module " << mod.GetModuleId() + << RESTendl; mod.GenerateGainMap(); if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) { mod.DrawSpectrum(); + mod.DrawFullSpectrum(); mod.DrawGainMap(); } } } ///////////////////////////////////////////// -/// \brief Function to calibrate a dataset. +/// \brief Function to calibrate a dataset with this gain map. /// /// \param dataSetFileName the name of the root file where the TRestDataSet to be calibrated is stored. /// \param outputFileName the name of the output (root) file where the calibrated TRestDataSet will be -/// exported. If empty, the output file will be named as the input file with the suffix "_cc". E.g. -/// "data/myDataSet.root" -> "data/myDataSet_cc.root". +/// exported. If empty, the output file will be named as the input file plus the name of the +/// TRestDataSetGainMap. E.g. "data/myDataSet.root" -> "data/myDataSet_.root". +/// \param excludeColumns a vector of strings with the names of the columns to be excluded from the +/// output file. If empty, all columns will be included. If "all" is in the list, all columns will be +/// excluded except the calibrated observable, the calibrated observable with no segmentation and +/// the plane-module identifier (pmID). /// -void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName) { +void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName, + std::vector excludeColumns) { if (fModulesCal.empty()) { RESTError << "TRestDataSetGainMap::CalibrateDataSet: No modules defined." << RESTendl; return; } TRestDataSet dataSet; + dataSet.EnableMultiThreading(true); dataSet.Import(dataSetFileName); auto dataFrame = dataSet.GetDataFrame(); @@ -253,18 +263,54 @@ void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, s dataFrame = dataFrame.Define(calibObsName, calibrate, {fObservable, fSpatialObservableX, fSpatialObservableY, pmIDname}); + // Define a new column with the calibrated observable for the whole module calibration + auto calibrateFullSpc = [this](double val, int pmID) { + for (auto& m : fModulesCal) { + if (pmID == m.GetPlaneId() * 10 + m.GetModuleId()) + return m.GetSlopeFullSpc() * val + m.GetInterceptFullSpc(); + } + // RESTError << "TRestDataSetGainMap::CalibrateDataSet: Module not found for pmID " << pmID << + // RESTendl; + return std::numeric_limits::quiet_NaN(); + }; + std::string calibObsNameFullSpc = (std::string)GetName() + "_"; + calibObsNameFullSpc += + GetObservable().erase(0, GetObservable().find("_") + 1); // remove the "rawAna_" part + calibObsNameFullSpc += "_NoSegmentation"; + dataFrame = dataFrame.Define(calibObsNameFullSpc, calibrateFullSpc, {fObservable, pmIDname}); + dataSet.SetDataFrame(dataFrame); // Format the output file name and export the dataSet if (outputFileName.empty()) outputFileName = dataSetFileName; - if (outputFileName == dataSetFileName) { // TRestDataSet cannot be overwritten - outputFileName = outputFileName.substr(0, outputFileName.find_last_of(".")) + "_cc."; + std::string gmName = GetName(); + outputFileName = outputFileName.substr(0, outputFileName.find_last_of(".")); // remove extension + outputFileName += "_" + gmName; outputFileName += TRestTools::GetFileNameExtension(dataSetFileName); } - RESTInfo << "Exporting calibrated dataSet to " << outputFileName << RESTendl; - dataSet.Export(outputFileName); + // Export dataset. Exclude columns if requested. + std::set excludeCol; // vector with the explicit column names to be excluded + auto columns = dataSet.GetDataFrame().GetColumnNames(); + // Get the columns to be excluded from the list of columns. It accepts wildcards "*" and "?" + for (auto& eC : excludeColumns) { + if (eC.find("*") != std::string::npos || eC.find("?") != std::string::npos) { + for (auto& c : columns) + if (MatchString(c, eC)) excludeCol.insert(c); + } else if (std::find(columns.begin(), columns.end(), eC) != columns.end()) + excludeCol.insert(eC); + } + // Remove the calibObsName, calibObsNameFullSpc and pmIDname from the list of columns to be excluded + excludeCol.erase(calibObsName); + excludeCol.erase(calibObsNameFullSpc); + excludeCol.erase(pmIDname); + + RESTDebug << "Excluding columns: "; + for (auto& c : excludeCol) RESTDebug << c << ", "; + RESTDebug << RESTendl; + + dataSet.Export(outputFileName, std::vector(excludeCol.begin(), excludeCol.end())); // Add this TRestDataSetGainMap metadata to the output file TFile* f = TFile::Open(outputFileName.c_str(), "UPDATE"); @@ -273,6 +319,21 @@ void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, s delete f; } +///////////////////////////////////////////// +/// \brief Function to retrieve the module calibration by index. Default is 0. +/// +/// +TRestDataSetGainMap::Module* TRestDataSetGainMap::GetModule(const size_t index) { + if (index < fModulesCal.size()) return &fModulesCal[index]; + + RESTError << "No ModuleCalibration with index " << index; + if (fModulesCal.empty()) + RESTError << ". There are no modules defined." << RESTendl; + else + RESTError << ". Max index is " << fModulesCal.size() - 1 << RESTendl; + return nullptr; +} + ///////////////////////////////////////////// /// \brief Function to retrieve the module calibration with planeID and moduleID /// @@ -297,6 +358,16 @@ double TRestDataSetGainMap::GetSlopeParameter(const int planeID, const int modul return moduleCal->GetSlope(x, y); } +///////////////////////////////////////////// +/// \brief Function to get the slope parameter of the whole module with +/// planeID and moduleID. +/// +double TRestDataSetGainMap::GetSlopeParameterFullSpc(const int planeID, const int moduleID) { + Module* moduleCal = GetModule(planeID, moduleID); + if (moduleCal == nullptr) return 0; // return numeric_limits::quiet_NaN() + return moduleCal->GetSlopeFullSpc(); +} + ///////////////////////////////////////////// /// \brief Function to get the intercept parameter of the module with /// planeID and moduleID at physical position (x,y) @@ -309,6 +380,16 @@ double TRestDataSetGainMap::GetInterceptParameter(const int planeID, const int m return moduleCal->GetIntercept(x, y); } +///////////////////////////////////////////// +/// \brief Function to get the intercept parameter of the whole module with +/// planeID and moduleID. +/// +double TRestDataSetGainMap::GetInterceptParameterFullSpc(const int planeID, const int moduleID) { + Module* moduleCal = GetModule(planeID, moduleID); + if (moduleCal == nullptr) return 0; // return numeric_limits::quiet_NaN() + return moduleCal->GetInterceptFullSpc(); +} + ///////////////////////////////////////////// /// \brief Function to get a list (set) of the plane IDs /// @@ -346,6 +427,7 @@ TRestDataSetGainMap& TRestDataSetGainMap::operator=(TRestDataSetGainMap& src) { fObservable = src.GetObservable(); fSpatialObservableX = src.GetSpatialObservableX(); fSpatialObservableY = src.GetSpatialObservableY(); + fCut = src.GetCut(); fModulesCal.clear(); for (auto pID : src.GetPlaneIDs()) for (auto mID : src.GetModuleIDs(pID)) fModulesCal.push_back(*src.GetModule(pID, mID)); @@ -356,7 +438,7 @@ TRestDataSetGainMap& TRestDataSetGainMap::operator=(TRestDataSetGainMap& src) { /// \brief Function to set a module calibration. If the module calibration /// already exists (same planeId and moduleId), it will be replaced. /// -void TRestDataSetGainMap::SetModuleCalibration(const Module& moduleCal) { +void TRestDataSetGainMap::SetModule(const Module& moduleCal) { for (auto& i : fModulesCal) { if (i.GetPlaneId() == moduleCal.GetPlaneId() && i.GetModuleId() == moduleCal.GetModuleId()) { i = moduleCal; @@ -431,7 +513,16 @@ void TRestDataSetGainMap::Export(const std::string& fileName) { void TRestDataSetGainMap::PrintMetadata() { TRestMetadata::PrintMetadata(); RESTMetadata << " Calibration dataset: " << fCalibFileName << RESTendl; + if (fCut) { + RESTMetadata << " Cuts applied: "; + /* print only cutStrings and paramCut because + TRestDataSet::MakeCut() uses fCut->GetCutStrings() and fCut->GetParamCut() */ + for (const auto& cut : fCut->GetCutStrings()) RESTMetadata << cut << ", " << RESTendl; + for (const auto& cut : fCut->GetParamCut()) + RESTMetadata << cut.first << " " << cut.second << ", " << RESTendl; + } RESTMetadata << " Output file: " << fOutputFileName << RESTendl; + RESTMetadata << RESTendl; RESTMetadata << " Number of planes: " << GetNumberOfPlanes() << RESTendl; RESTMetadata << " Number of modules: " << GetNumberOfModules() << RESTendl; RESTMetadata << " Calibration observable: " << fObservable << RESTendl; @@ -628,10 +719,11 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { } if (!TRestTools::isDataSet(dsFileName)) RESTWarning << dsFileName << " is not a dataset." << p->RESTendl; TRestDataSet dataSet; + dataSet.EnableMultiThreading(true); dataSet.Import(dsFileName); fDataSetFileName = dsFileName; - SetSplits(); + dataSet.SetDataFrame(dataSet.MakeCut(p->GetCut())); if (fSplitX.empty()) SetSplitX(); if (fSplitY.empty()) SetSplitY(); @@ -663,20 +755,28 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { << p->RESTendl; } + // --- Definition of histogram whole module --- + std::string hModuleName = "hSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); + delete fFullSpectrum; + fFullSpectrum = new TH1F(hModuleName.c_str(), "", fNBins, fCalibRange.X(), fCalibRange.Y()); + + // build the spectrum for the whole module + std::string cut = fDefinitionCut; + if (cut.empty()) cut = "1"; + auto histoMod = dataSet.GetDataFrame().Filter(cut).Histo1D( + {"tempMod", "", fNBins, fCalibRange.X(), fCalibRange.Y()}, GetObservable()); + std::unique_ptr hpuntMod = std::unique_ptr(static_cast(histoMod->Clone())); + fFullSpectrum->Add(hpuntMod.get()); + //--- Definition of histogram matrix --- std::vector> h(fNumberOfSegmentsX, std::vector(fNumberOfSegmentsY, nullptr)); for (size_t i = 0; i < h.size(); i++) { for (size_t j = 0; j < h.at(0).size(); j++) { - std::string name = "hSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId) + "_" + - std::to_string(i) + "_" + std::to_string(j); + std::string name = hModuleName + "_" + std::to_string(i) + "_" + std::to_string(j); h[i][j] = new TH1F(name.c_str(), "", fNBins, fCalibRange.X(), fCalibRange.Y()); // h[column][row] equivalent to h[x][y] } } - std::vector> calParSlope(fNumberOfSegmentsX, - std::vector(fNumberOfSegmentsY, 0)); - std::vector> calParIntercept(fNumberOfSegmentsX, - std::vector(fNumberOfSegmentsY, 0)); // build the spectrum for each segment auto itX = fSplitX.begin(); @@ -713,111 +813,130 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { } //--- Fit every peak energy for every segment --- + std::vector> calParSlope(fNumberOfSegmentsX, + std::vector(fNumberOfSegmentsY, 0)); + std::vector> calParIntercept(fNumberOfSegmentsX, + std::vector(fNumberOfSegmentsY, 0)); fSegLinearFit = std::vector(h.size(), std::vector(h.at(0).size(), nullptr)); - for (size_t i = 0; i < h.size(); i++) { - for (size_t j = 0; j < h.at(0).size(); j++) { - RESTExtreme << "Segment[" << i << "][" << j << "]" << p->RESTendl; - // Search for peaks --> peakPos - std::unique_ptr s(new TSpectrum(2 * fEnergyPeaks.size() + 1)); - std::vector peakPos; - s->Search(h[i][j], 2, "goff", 0.1); - for (int k = 0; k < s->GetNPeaks(); k++) peakPos.push_back(s->GetPositionX()[k]); - std::sort(peakPos.begin(), peakPos.end(), std::greater()); - const double ratio = peakPos.size() == 0 - ? 1 - : peakPos.front() / fEnergyPeaks.front(); // to estimate peak position - - // Initialize graph for linear fit - std::shared_ptr gr; - gr = std::shared_ptr(new TGraph()); - gr->SetName("grFit"); - gr->SetTitle((";" + GetObservable() + ";energy").c_str()); - - // Fit every energy peak - int c = 0; - double mu = 0; - for (const auto& energy : fEnergyPeaks) { - RESTExtreme << "\t fitting energy " << DoubleToString(energy, "%g") << p->RESTendl; - // estimation of the peak position is between start and end - double pos = energy * ratio; - double start = pos * 0.8; - double end = pos * 1.2; - if (fRangePeaks.at(c).X() < fRangePeaks.at(c).Y()) { // if range is defined use it - start = fRangePeaks.at(c).X(); - end = fRangePeaks.at(c).Y(); - } + for (size_t i = 0; i < h.size(); i++) + for (int j = 0; j < (int)h.at(0).size(); j++) { + fSegLinearFit[i][j] = new TGraph(); + auto [intercept, slope] = FitPeaks(h[i][j], fSegLinearFit[i][j]); + calParSlope[i][j] = slope; + calParIntercept[i][j] = intercept; + } + fSlope = calParSlope; + fIntercept = calParIntercept; + fSegSpectra = h; + + //--- Fit every peak energy for the whole module --- + delete fFullLinearFit; + fFullLinearFit = new TGraph(); + auto [intercept, slope] = FitPeaks(fFullSpectrum, fFullLinearFit); + fFullSlope = slope; + fFullIntercept = intercept; +} + +std::pair TRestDataSetGainMap::Module::FitPeaks(TH1F* hSeg, TGraph* gr) { + if (!hSeg) { + RESTError << "No histogram for fitting" << p->RESTendl; + return std::make_pair(0, 0); + } + if (hSeg->Integral() == 0) { + RESTError << "Empty spectrum " << hSeg->GetName() << p->RESTendl; + return std::make_pair(0, 0); + } + std::shared_ptr graph = std::shared_ptr(new TGraph()); + RESTExtreme << "Fitting peaks for " << hSeg->GetName() << p->RESTendl; + + // Search for peaks --> peakPos + std::unique_ptr s(new TSpectrum(2 * fEnergyPeaks.size() + 1)); + std::vector peakPos; + s->Search(hSeg, 2, "goff", 0.1); + for (int k = 0; k < s->GetNPeaks(); k++) peakPos.push_back(s->GetPositionX()[k]); + std::sort(peakPos.begin(), peakPos.end(), std::greater()); + const double ratio = + peakPos.size() == 0 ? 1 : peakPos.front() / fEnergyPeaks.front(); // to estimate peak position + + // Initialize graph for linear fit + graph->SetName("grFit"); + graph->SetTitle((";" + GetObservable() + ";energy").c_str()); + + // Fit every energy peak + int c = 0; + double mu = 0; + for (const auto& energy : fEnergyPeaks) { + RESTExtreme << "\t fitting energy " << DoubleToString(energy, "%g") << p->RESTendl; + // estimation of the peak position is between start and end + double pos = energy * ratio; + double start = pos * 0.8; + double end = pos * 1.2; + if (fRangePeaks.at(c).X() < fRangePeaks.at(c).Y()) { // if range is defined use it + start = fRangePeaks.at(c).X(); + end = fRangePeaks.at(c).Y(); + } - do { - if (fAutoRangePeaks) { - if (peakPos.size() > 0) { - // Find the peak position that is between start and end + do { + if (fAutoRangePeaks) { + if (peakPos.size() > 0) { + // Find the peak position that is between start and end + pos = peakPos.at(0); + while (!(start < pos && pos < end)) { + // if none of the peak position is + // between start and end, use the greater one. + if (pos == peakPos.back()) { pos = peakPos.at(0); - while (!(start < pos && pos < end)) { - // if none of the peak position is - // between start and end, use the greater one. - if (pos == peakPos.back()) { - pos = peakPos.at(0); - break; - } - pos = *std::next(std::find(peakPos.begin(), peakPos.end(), - pos)); // get next peak position - } - peakPos.erase(std::find(peakPos.begin(), peakPos.end(), - pos)); // remove this peak position from the list - // final estimation of the peak range (idem fitting window) with this peak - // position pos - start = pos * 0.8; - end = pos * 1.2; - const double relDist = peakPos.size() > 0 ? (pos - peakPos.front()) / pos : 999; - if (relDist < 0.2) { // if the next peak is too close reduce the window width - start = pos * (1 - relDist / 2); - end = pos * (1 + relDist / 2); - } + break; } + pos = *std::next(std::find(peakPos.begin(), peakPos.end(), + pos)); // get next peak position } - - std::string name = "g" + std::to_string(c); - TF1* g = new TF1(name.c_str(), "gaus", start, end); - RESTExtreme << "\t\tat " << DoubleToString(pos, "%.3g") << ". Range(" - << DoubleToString(start, "%.3g") << ", " << DoubleToString(end, "%.3g") << ")" - << p->RESTendl; - - if (h[i][j]->GetFunction(name.c_str())) // remove previous fit - h[i][j]->GetListOfFunctions()->Remove(h[i][j]->GetFunction(name.c_str())); - - h[i][j]->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it - mu = g->GetParameter(1); - RESTExtreme << "\t\tgaus mean " << DoubleToString(mu, "%g") << p->RESTendl; - } while (fAutoRangePeaks && peakPos.size() > 0 && - !(start < mu && mu < end)); // avoid small peaks on main peak tail - gr->SetPoint(c++, mu, energy); - } - s.reset(); // delete s; - - if (fZeroPoint) gr->SetPoint(c++, 0, 0); - while (gr->GetN() < 2) { // minimun 2 points needed for linear fit - gr->SetPoint(c++, 0, 0); - SetZeroPoint(true); - RESTDebug << "Not enough points for linear fit. Adding and setting zero point to true" - << p->RESTendl; + peakPos.erase(std::find(peakPos.begin(), peakPos.end(), + pos)); // remove this peak position from the list + // final estimation of the peak range (idem fitting window) with this peak + // position pos + start = pos * 0.8; + end = pos * 1.2; + const double relDist = peakPos.size() > 0 ? (pos - peakPos.front()) / pos : 999; + if (relDist < 0.2) { // if the next peak is too close reduce the window width + start = pos * (1 - relDist / 2); + end = pos * (1 + relDist / 2); + } + } } - // Linear fit - std::unique_ptr linearFit; - linearFit = std::unique_ptr(new TF1("linearFit", "pol1")); - gr->Fit("linearFit", "SQ"); // Q for quiet mode + std::string name = "g" + std::to_string(c); + TF1* g = new TF1(name.c_str(), "gaus", start, end); + RESTExtreme << "\t\tat " << DoubleToString(pos, "%.3g") << ". Range(" + << DoubleToString(start, "%.3g") << ", " << DoubleToString(end, "%.3g") << ")" + << p->RESTendl; - fSegLinearFit.at(i).at(j) = (TGraph*)gr->Clone(); + if (hSeg->GetFunction(name.c_str())) // remove previous fit + hSeg->GetListOfFunctions()->Remove(hSeg->GetFunction(name.c_str())); - const double slope = linearFit->GetParameter(1); - const double intercept = linearFit->GetParameter(0); - calParSlope.at(i).at(j) = slope; - calParIntercept.at(i).at(j) = intercept; - } + hSeg->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it + mu = g->GetParameter(1); + RESTExtreme << "\t\tgaus mean " << DoubleToString(mu, "%g") << p->RESTendl; + } while (fAutoRangePeaks && peakPos.size() > 0 && + !(start < mu && mu < end)); // avoid small peaks on main peak tail + graph->SetPoint(c++, mu, energy); } - fSegSpectra = h; - fSlope = calParSlope; - fIntercept = calParIntercept; + s.reset(); // delete s; + + if (fZeroPoint) graph->SetPoint(c++, 0, 0); + while (graph->GetN() < 2) { // minimun 2 points needed for linear fit + graph->SetPoint(c++, 0, 0); + SetZeroPoint(true); + RESTDebug << "Not enough points for linear fit. Adding and setting zero point to true" << p->RESTendl; + } + + // Linear fit + std::unique_ptr linearFit; + linearFit = std::unique_ptr(new TF1("linearFit", "pol1")); + graph->Fit("linearFit", "SQ"); // Q for quiet mode + + if (gr) *gr = *(TGraph*)graph->Clone(); // if nullptr is passed, do not copy the graph + return std::make_pair(linearFit->GetParameter(0), linearFit->GetParameter(1)); } ///////////////////////////////////////////// @@ -879,6 +998,55 @@ void TRestDataSetGainMap::Module::Refit(const size_t x, const size_t y, const si UpdateCalibrationFits(x, y); } +///////////////////////////////////////////// +/// \brief Function to fit again manually a peak for the whole module spectrum. The +/// calibration curve is updated after the fit. +/// +/// \param energyPeak The energy of the peak to be fitted (in physical units). +/// \param range The range for the fitting of the peak (in the observables corresponding units). +/// +void TRestDataSetGainMap::Module::RefitFullSpc(const double energyPeak, const TVector2& range) { + int peakNumber = -1; + for (size_t i = 0; i < fEnergyPeaks.size(); i++) + if (fEnergyPeaks.at(i) == energyPeak) { + peakNumber = i; + break; + } + if (peakNumber == -1) { + RESTError << "Energy " << energyPeak << " not found in the list of energy peaks" << p->RESTendl; + return; + } + RefitFullSpc((size_t)peakNumber, range); +} + +///////////////////////////////////////////// +/// \brief Function to fit again manually a peak for the whole module spectrum. The +/// calibration curve is updated after the fit. +/// +/// \param peakNumber The index of the peak to be fitted. +/// \param range The range for the fitting of the peak (in the observables corresponding units). +/// +void TRestDataSetGainMap::Module::RefitFullSpc(const size_t peakNumber, const TVector2& range) { + if (!fFullSpectrum) { + RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl; + return; + } + if (peakNumber >= fEnergyPeaks.size()) { + RESTError << "Peak with index " << peakNumber << " not found" << p->RESTendl; + return; + } + + // Refit the desired peak + std::string name = "g" + std::to_string(peakNumber); + TF1* g = new TF1(name.c_str(), "gaus", range.X(), range.Y()); + while (fFullSpectrum->GetFunction(name.c_str())) // clear previous fits for this peakNumber + fFullSpectrum->GetListOfFunctions()->Remove(fFullSpectrum->GetFunction(name.c_str())); + fFullSpectrum->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it + + // Change the point of the graph + UpdateCalibrationFitsFullSpc(); +} + ///////////////////////////////////////////// /// \brief Function to update the calibration curve for a given segment of the module. The calibration /// curve is cleared and then the means of the gaussian fits for each energy peak are added. If there are @@ -900,6 +1068,25 @@ void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const si TH1F* h = fSegSpectra.at(x).at(y); TGraph* gr = fSegLinearFit.at(x).at(y); + auto [intercept, slope] = UpdateCalibrationFits(h, gr); + fSlope[x][y] = slope; + fIntercept[x][y] = intercept; +} + +std::pair TRestDataSetGainMap::Module::UpdateCalibrationFits(TH1F* h, TGraph* gr) { + if (!h) { + RESTError << "No histogram for updating fits" << p->RESTendl; + return std::make_pair(0, 0); + } + if (!gr) { + RESTError << "No graph for updating fits" << p->RESTendl; + return std::make_pair(0, 0); + } + if (h->Integral() == 0) { + RESTError << "Empty spectrum " << h->GetName() << p->RESTendl; + return std::make_pair(0, 0); + } + // Clear the points of the graph for (size_t i = 0; i < fEnergyPeaks.size(); i++) gr->RemovePoint(i); // Add the new points to the graph @@ -909,7 +1096,7 @@ void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const si TF1* g = h->GetFunction(fitName.c_str()); if (!g) { RESTWarning << "No fit ( " << fitName << " ) found for energy peak " << fEnergyPeaks[i] - << " in segment " << x << "," << y << p->RESTendl; + << " in histogram " << h->GetName() << p->RESTendl; continue; } gr->SetPoint(c++, g->GetParameter(1), fEnergyPeaks[i]); @@ -918,8 +1105,6 @@ void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const si // Add zero points if needed (if there are less than 2 points) while (gr->GetN() < 2) { gr->SetPoint(c++, 0, 0); - RESTWarning << "Not enough points for linear fit at segment (" << x << ", " << y - << "). Adding zero point." << p->RESTendl; } // Refit the calibration curve @@ -929,8 +1114,19 @@ void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const si else lf = new TF1("linearFit", "pol1"); gr->Fit(lf, "SQ"); // Q for quiet mode - fSlope.at(x).at(y) = lf->GetParameter(1); - fIntercept.at(x).at(y) = lf->GetParameter(0); + + return std::make_pair(lf->GetParameter(0), lf->GetParameter(1)); +} + +///////////////////////////////////////////// +/// \brief Function to update the calibration curve for the whole module. The calibration +/// curve is cleared and then the means of the gaussian fits for each energy peak are added. If there are +/// less than 2 fits, zero points are added. Then, the calibration curve is refitted (linearFit). +/// +void TRestDataSetGainMap::Module::UpdateCalibrationFitsFullSpc() { + auto [intercept, slope] = UpdateCalibrationFits(fFullSpectrum, fFullLinearFit); + fFullSlope = slope; + fFullIntercept = intercept; } ///////////////////////////////////////////// @@ -1021,6 +1217,10 @@ void TRestDataSetGainMap::Module::DrawSpectrum(const int index_x, const int inde RESTError << "Index out of range." << p->RESTendl; return; } + if (!fSegSpectra[index_x][index_y]) { + RESTError << "No Spectrum for segment (" << index_x << ", " << index_y << ")." << p->RESTendl; + return; + } if (!c) { std::string t = "spectrum_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId) + "_" + @@ -1044,7 +1244,8 @@ void TRestDataSetGainMap::Module::DrawSpectrum(const int index_x, const int inde if (drawFits) for (size_t c = 0; c < fEnergyPeaks.size(); c++) { auto fit = fSegSpectra[index_x][index_y]->GetFunction(("g" + std::to_string(c)).c_str()); - if (!fit) RESTWarning << "Fit for energy peak" << fEnergyPeaks[c] << " not found." << p->RESTendl; + if (!fit) + RESTWarning << "Fit for energy peak " << fEnergyPeaks[c] << " not found." << p->RESTendl; if (!fit) continue; fit->SetLineColor(c + 2 != colorT ? c + 2 : ++colorT); /* does not work with kRed, kBlue, etc. as they are not defined with the same number as the first 10 basic colors. See @@ -1095,30 +1296,41 @@ void TRestDataSetGainMap::Module::DrawSpectrum(const bool drawFits, const int co for (size_t i = 0; i < fSegSpectra.size(); i++) { for (size_t j = 0; j < fSegSpectra[i].size(); j++) { - c->cd(i + 1 + fSegSpectra[i].size() * j); - DrawSpectrum(i, fSegSpectra[i].size() - 1 - j, drawFits, color, c); + int pad = fSegSpectra.size() * (fSegSpectra[i].size() - 1) + 1 + i - fSegSpectra.size() * j; + c->cd(pad); + DrawSpectrum(i, j, drawFits, color, c); } } } -void TRestDataSetGainMap::Module::DrawFullSpectrum() { - if (fSegSpectra.size() == 0) { - RESTError << "Spectra matrix is empty." << p->RESTendl; +void TRestDataSetGainMap::Module::DrawFullSpectrum(const bool drawFits, const int color, TCanvas* c) { + if (!fFullSpectrum) { + RESTError << "Spectrum is empty." << p->RESTendl; return; } - TH1F* sumHist = - new TH1F("fullSpc", "", fSegSpectra[0][0]->GetNbinsX(), fSegSpectra[0][0]->GetXaxis()->GetXmin(), - fSegSpectra[0][0]->GetXaxis()->GetXmax()); - sumHist->SetTitle(("Full spectrum;" + GetObservable() + ";counts").c_str()); - for (size_t i = 0; i < fSegSpectra.size(); i++) { - for (size_t j = 0; j < fSegSpectra.at(0).size(); j++) { - sumHist->Add(fSegSpectra[i][j]); - } + if (!c) { + std::string t = "fullSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); + c = new TCanvas(t.c_str(), t.c_str()); } - std::string t = "fullSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); - TCanvas* c = new TCanvas(t.c_str(), t.c_str()); c->cd(); - sumHist->Draw(); + + fFullSpectrum->SetTitle(("Full spectrum;" + GetObservable() + ";counts").c_str()); + + if (color > 0) fFullSpectrum->SetLineColor(color); + size_t colorT = fFullSpectrum->GetLineColor(); + fFullSpectrum->Draw("same"); + + if (drawFits) + for (size_t c = 0; c < fEnergyPeaks.size(); c++) { + auto fit = fFullSpectrum->GetFunction(("g" + std::to_string(c)).c_str()); + if (!fit) RESTWarning << "Fit for energy peak" << fEnergyPeaks[c] << " not found." << p->RESTendl; + if (!fit) continue; + fit->SetLineColor(c + 2 != colorT ? c + 2 : ++colorT); /* does not work with kRed, kBlue, etc. + as they are not defined with the same number as the first 10 basic colors. See + https://root.cern.ch/doc/master/classTColor.html#C01 and + https://root.cern.ch/doc/master/classTColor.html#C02 */ + fit->Draw("same"); + } } void TRestDataSetGainMap::Module::DrawLinearFit(const TVector2& position, TCanvas* c) { @@ -1136,6 +1348,11 @@ void TRestDataSetGainMap::Module::DrawLinearFit(const int index_x, const int ind RESTError << "Index out of range." << p->RESTendl; return; } + if (!fSegLinearFit[index_x][index_y]) { + RESTError << "No linear fit for segment (" << index_x << ", " << index_y << ")." << p->RESTendl; + return; + } + if (!c) { std::string t = "linearFit_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId) + "_" + std::to_string(index_x) + "_" + std::to_string(index_y); @@ -1152,25 +1369,58 @@ void TRestDataSetGainMap::Module::DrawLinearFit(const int index_x, const int ind fSegLinearFit[index_x][index_y]->Draw("AL*"); } -void TRestDataSetGainMap::Module::DrawLinearFit() { - std::string t = "linearFits_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); - TCanvas* myCanvas = new TCanvas(t.c_str(), t.c_str()); - myCanvas->Divide(fSegLinearFit.size(), fSegLinearFit.at(0).size()); +void TRestDataSetGainMap::Module::DrawLinearFit(TCanvas* c) { + if (fSegLinearFit.size() == 0) { + RESTError << "Spectra matrix is empty." << p->RESTendl; + return; + } + if (!c) { + std::string t = "linearFits_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); + c = new TCanvas(t.c_str(), t.c_str()); + } + + size_t nPads = 0; + for (const auto& object : *c->GetListOfPrimitives()) + if (object->InheritsFrom(TVirtualPad::Class())) ++nPads; + if (nPads != 0 && nPads != fSegLinearFit.size() * fSegLinearFit.at(0).size()) { + RESTError << "Canvas " << c->GetName() << " has " << nPads << " pads, but " + << fSegLinearFit.size() * fSegLinearFit.at(0).size() << " are needed." << p->RESTendl; + return; + } else if (nPads == 0) + c->Divide(fSegLinearFit.size(), fSegLinearFit.at(0).size()); + for (size_t i = 0; i < fSegLinearFit.size(); i++) { for (size_t j = 0; j < fSegLinearFit[i].size(); j++) { - myCanvas->cd(i + 1 + fSegLinearFit[i].size() * j); - // fSegLinearFit[i][j]->Draw("AL*"); - DrawLinearFit(i, fSegSpectra[i].size() - 1 - j, myCanvas); + int pad = fSegLinearFit.size() * (fSegLinearFit[i].size() - 1) + 1 + i - fSegLinearFit.size() * j; + c->cd(pad); + DrawLinearFit(i, j, c); } } } -void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { +///////////////////////////////////////////// +/// \brief Function to draw the relative gain map for a given energy peak of the module. +/// +/// \param peakNumber The index of the peak to be drawn (remember they are in descending order). +/// \param fullModuleAsRef If true, it will use the peak position at the full module spectrum +/// as reference for the gain map. If false, it will use the centered segment of the module +/// as reference. +/// +void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber, bool fullModuleAsRef) { if (peakNumber < 0 || peakNumber >= (int)fEnergyPeaks.size()) { RESTError << "Peak number out of range (peakNumber should be between 0 and " << fEnergyPeaks.size() - 1 << " )" << p->RESTendl; return; } + if (fSegLinearFit.size() == 0) { + RESTError << "Linear fit matrix is empty." << p->RESTendl; + return; + } + if (!fFullLinearFit) { + RESTError << "Full linear fit is empty." << p->RESTendl; + return; + } + double peakEnergy = fEnergyPeaks[peakNumber]; std::string title = "Gain map for energy " + DoubleToString(peakEnergy, "%g") + ";" + GetSpatialObservableX() + ";" + GetSpatialObservableY(); // + " keV"; @@ -1178,11 +1428,15 @@ void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { std::to_string(fModuleId); TCanvas* gainMap = new TCanvas(t.c_str(), t.c_str()); gainMap->cd(); - TH2F* hGainMap = new TH2F(("h" + t).c_str(), title.c_str(), fNumberOfSegmentsY, fReadoutRange.X(), - fReadoutRange.Y(), fNumberOfSegmentsX, fReadoutRange.X(), fReadoutRange.Y()); - - const double peakPosRef = - fSegLinearFit[(fNumberOfSegmentsX - 1) / 2][(fNumberOfSegmentsY - 1) / 2]->GetPointX(peakNumber); + TH2F* hGainMap = new TH2F(("h" + t).c_str(), title.c_str(), fNumberOfSegmentsX, fReadoutRange.X(), + fReadoutRange.Y(), fNumberOfSegmentsY, fReadoutRange.X(), fReadoutRange.Y()); + + double peakPosRef = fFullLinearFit->GetPointX(peakNumber); + if (!fullModuleAsRef) { + int index_x = fNumberOfSegmentsX > 0 ? (fNumberOfSegmentsX - 1) / 2 : 0; + int index_y = fNumberOfSegmentsY > 0 ? (fNumberOfSegmentsY - 1) / 2 : 0; + peakPosRef = fSegLinearFit[index_x][index_y]->GetPointX(peakNumber); + } auto itX = fSplitX.begin(); for (size_t i = 0; i < fSegLinearFit.size(); i++) { @@ -1192,8 +1446,11 @@ void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { auto xUpper = *std::next(itX); auto yLower = *itY; auto yUpper = *std::next(itY); - hGainMap->Fill((xUpper + xLower) / 2., (yUpper + yLower) / 2., - fSegLinearFit[i][j]->GetPointX(peakNumber) / peakPosRef); + float xMean = (xUpper + xLower) / 2.; + float yMean = (yUpper + yLower) / 2.; + auto [index_x, index_y] = GetIndexMatrix(xMean, yMean); + if (!fSegLinearFit[index_x][index_y]) continue; + hGainMap->Fill(xMean, yMean, fSegLinearFit[index_x][index_y]->GetPointX(peakNumber) / peakPosRef); itY++; } itX++; @@ -1281,5 +1538,9 @@ void TRestDataSetGainMap::Module::Print() const { } RESTMetadata << p->RESTendl; } + RESTMetadata << p->RESTendl; + RESTMetadata << " Full slope: " << DoubleToString(fFullSlope, "%.3e") << p->RESTendl; + RESTMetadata << " Full intercept: " << DoubleToString(fFullIntercept, "%+.3e") << p->RESTendl; + RESTMetadata << "-----------------------------------------------" << p->RESTendl; } diff --git a/source/framework/core/inc/TRestDataSet.h b/source/framework/core/inc/TRestDataSet.h index 483351510..c70eae202 100644 --- a/source/framework/core/inc/TRestDataSet.h +++ b/source/framework/core/inc/TRestDataSet.h @@ -112,7 +112,7 @@ class TRestDataSet : public TRestMetadata { Bool_t fExternal = false; //< /// The resulting RDF::RNode object after initialization - ROOT::RDF::RNode fDataSet = ROOT::RDataFrame(0); //! + ROOT::RDF::RNode fDataFrame = ROOT::RDataFrame(0); //! /// A pointer to the generated tree TChain* fTree = nullptr; //! @@ -122,12 +122,14 @@ class TRestDataSet : public TRestMetadata { protected: virtual std::vector FileSelection(); + void RegenerateTree(std::vector finalList = {}); + public: /// Gives access to the RDataFrame ROOT::RDF::RNode GetDataFrame() const { if (!fExternal && fTree == nullptr) RESTWarning << "DataFrame has not been yet initialized" << RESTendl; - return fDataSet; + return fDataFrame; } void EnableMultiThreading(Bool_t enable = true) { fMT = enable; } @@ -152,7 +154,7 @@ class TRestDataSet : public TRestMetadata { } /// Number of variables (or observables) - size_t GetNumberOfColumns() { return fDataSet.GetColumnNames().size(); } + size_t GetNumberOfColumns() { return fDataFrame.GetColumnNames().size(); } /// Number of variables (or observables) size_t GetNumberOfBranches() { return GetNumberOfColumns(); } @@ -187,7 +189,7 @@ class TRestDataSet : public TRestMetadata { void SetTotalTimeInSeconds(Double_t seconds) { fTotalDuration = seconds; } void SetDataFrame(const ROOT::RDF::RNode& dS) { - fDataSet = dS; + fDataFrame = dS; fExternal = true; } @@ -198,8 +200,12 @@ class TRestDataSet : public TRestMetadata { void Export(const std::string& filename, std::vector excludeColumns = {}); ROOT::RDF::RNode MakeCut(const TRestCut* cut); + ROOT::RDF::RNode ApplyRange(size_t from, size_t to); + ROOT::RDF::RNode Range(size_t from, size_t to); ROOT::RDF::RNode DefineColumn(const std::string& columnName, const std::string& formula); + size_t GetEntries(); + void PrintMetadata() override; void Initialize() override; @@ -209,6 +215,6 @@ class TRestDataSet : public TRestMetadata { TRestDataSet(const char* cfgFileName, const std::string& name = ""); ~TRestDataSet(); - ClassDefOverride(TRestDataSet, 7); + ClassDefOverride(TRestDataSet, 8); }; #endif diff --git a/source/framework/core/inc/TRestDataSetPlot.h b/source/framework/core/inc/TRestDataSetPlot.h index fec0c6266..d957631b2 100644 --- a/source/framework/core/inc/TRestDataSetPlot.h +++ b/source/framework/core/inc/TRestDataSetPlot.h @@ -95,6 +95,7 @@ class TRestDataSetPlot : public TRestMetadata { std::vector, TVector2>> variablePos; std::vector, TVector2>> metadataPos; std::vector, TVector2>> obsPos; + std::vector, TVector2>> expPos; TRestCut* panelCut = nullptr; @@ -182,6 +183,6 @@ class TRestDataSetPlot : public TRestMetadata { TRestDataSetPlot(const char* configFilename, std::string name = ""); ~TRestDataSetPlot(); - ClassDefOverride(TRestDataSetPlot, 1); + ClassDefOverride(TRestDataSetPlot, 2); }; #endif diff --git a/source/framework/core/inc/TRestMetadata.h b/source/framework/core/inc/TRestMetadata.h index 561256f99..e0e51018d 100644 --- a/source/framework/core/inc/TRestMetadata.h +++ b/source/framework/core/inc/TRestMetadata.h @@ -286,6 +286,10 @@ class TRestMetadata : public TNamed { inline Bool_t isCleanState() const { return fCleanState; } + UInt_t GetVersionMajor() const; + UInt_t GetVersionMinor() const; + UInt_t GetVersionPatch() const; + Int_t GetVersionCode(); /// Returns a std::string with the path used for data storage inline TString GetDataPath() { diff --git a/source/framework/core/inc/TRestVersion.h b/source/framework/core/inc/TRestVersion.h index d7d1c8f66..eb622bc7f 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.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_RELEASE "2.4.3" +#define REST_RELEASE_DATE "Fri May 3" +#define REST_RELEASE_TIME "17:36:10 CEST 2024" +#define REST_RELEASE_NAME "Steven Weinberg" +#define REST_GIT_COMMIT "bc42645d" +#define REST_VERSION_CODE 132099 #define REST_VERSION(a, b, c) (((a) << 16) + ((b) << 8) + (c)) #define REST_SCHEMA_EVOLUTION "ON" #endif diff --git a/source/framework/core/src/TRestDataSet.cxx b/source/framework/core/src/TRestDataSet.cxx index 56011722a..09c2eaa95 100644 --- a/source/framework/core/src/TRestDataSet.cxx +++ b/source/framework/core/src/TRestDataSet.cxx @@ -169,6 +169,26 @@ /// [2] d.GetTree()->GetEntries() /// \endcode /// +/// Example 3 Automatically importing a dataset using restRoot +/// +/// \code +/// restRoot Dataset_FinalBabyIAXO_XMM_mM_P14.root +/// +/// REST dataset file identified. It contains a valid TRestDataSet. +/// +/// Importing dataset /nfs/dust/iaxo/user/jgalan//Dataset_FinalBabyIAXO_XMM_mM_P14.root as `dSet` +/// +/// The dataset is ready. You may now access the dataset using: +/// +/// - dSet->PrintMetadata() +/// - dSet->GetDataFrame().GetColumnNames() +/// - dSet->GetTree()->GetEntries() +/// +/// - dSet->GetDataFrame().Display("")->Print() +/// - dSet->GetDataFrame().Display({"colName1,colName2"})->Print() +/// [0] +/// \endcode +/// /// ### Relevant quantities /// /// Sometimes we will be willing that our dataset contains few variables @@ -241,6 +261,21 @@ /// where `SolarFlux`,`GeneratorArea` and `Nsim` are the given names of /// the relevant quantities inside the dataset. /// +/// ### Adding cuts +/// +/// It is also possible to add cuts used to filter the data that will +/// be stored inside the dataset. We can do that including a TRestCut +/// definition inside the TRestDataSet. +/// +/// For example, the following cut definition would discard entries +/// with unexpected values inside the specified column, `process_status`. +/// +/// \code +/// +/// +/// +/// \endcode +/// ///---------------------------------------------------------------------- /// /// REST-for-Physics - Software for Rare Event Searches Toolkit @@ -347,30 +382,40 @@ void TRestDataSet::GenerateDataSet() { ROOT::DisableImplicitMT(); RESTInfo << "Initializing dataset" << RESTendl; - fDataSet = ROOT::RDataFrame("AnalysisTree", fFileSelection); + fDataFrame = ROOT::RDataFrame("AnalysisTree", fFileSelection); RESTInfo << "Making cuts" << RESTendl; - fDataSet = MakeCut(fCut); + fDataFrame = MakeCut(fCut); // Adding new user columns added to the dataset for (const auto& [cName, cExpression] : fColumnNameExpressions) { RESTInfo << "Adding column to dataset: " << cName << RESTendl; finalList.emplace_back(cName); - fDataSet = DefineColumn(cName, cExpression); + fDataFrame = DefineColumn(cName, cExpression); } + RegenerateTree(finalList); + + RESTInfo << " - Dataset generated!" << RESTendl; +} + +/////////////////////////////////////////////// +/// \brief It regenerates the tree so that it is an exact copy of the present DataFrame +/// +void TRestDataSet::RegenerateTree(std::vector finalList) { RESTInfo << "Generating snapshot." << RESTendl; std::string user = getenv("USER"); std::string fOutName = "/tmp/rest_output_" + user + ".root"; - fDataSet.Snapshot("AnalysisTree", fOutName, finalList); + if (!finalList.empty()) + fDataFrame.Snapshot("AnalysisTree", fOutName, finalList); + else + fDataFrame.Snapshot("AnalysisTree", fOutName); RESTInfo << "Re-importing analysis tree." << RESTendl; - fDataSet = ROOT::RDataFrame("AnalysisTree", fOutName); + fDataFrame = ROOT::RDataFrame("AnalysisTree", fOutName); TFile* f = TFile::Open(fOutName.c_str()); fTree = (TChain*)f->Get("AnalysisTree"); - - RESTInfo << " - Dataset generated!" << RESTendl; } /////////////////////////////////////////////// @@ -483,13 +528,31 @@ std::vector TRestDataSet::FileSelection() { } /////////////////////////////////////////////// -/// \brief This function apply a TRestCut to the dataframe +/// \brief This method returns a RDataFrame node with the number of +/// samples inside the dataset by selecting a range. It will not +/// modify internally the dataset. See ApplyRange to modify internally +/// the dataset. +/// +ROOT::RDF::RNode TRestDataSet::Range(size_t from, size_t to) { return fDataFrame.Range(from, to); } + +/////////////////////////////////////////////// +/// \brief This method reduces the number of samples inside the +/// dataset by selecting a range. +/// +ROOT::RDF::RNode TRestDataSet::ApplyRange(size_t from, size_t to) { + fDataFrame = fDataFrame.Range(from, to); + RegenerateTree(); + return fDataFrame; +} + +/////////////////////////////////////////////// +/// \brief This function applies a TRestCut to the dataframe /// and returns a dataframe with the applied cuts. Note that /// the cuts are not applied directly to the dataframe on -/// TRestDataSet, to do so you should do fDataSet = MakeCut(fCut); +/// TRestDataSet, to do so you should do fDataFrame = MakeCut(fCut); /// ROOT::RDF::RNode TRestDataSet::MakeCut(const TRestCut* cut) { - auto df = fDataSet; + auto df = fDataFrame; if (cut == nullptr) return df; @@ -526,6 +589,20 @@ ROOT::RDF::RNode TRestDataSet::MakeCut(const TRestCut* cut) { return df; } +/////////////////////////////////////////////// +/// \brief It returns the number of entries found inside fDataFrame +/// and prints out a warning if the number of entries inside the +/// tree is not the same. +/// +size_t TRestDataSet::GetEntries() { + auto nEntries = fDataFrame.Count(); + if (*nEntries == (long long unsigned int)GetTree()->GetEntries()) return *nEntries; + RESTWarning << "TRestDataSet::GetEntries. Number of tree entries is not the same as RDataFrame entries." + << RESTendl; + RESTWarning << "Returning RDataFrame entries" << RESTendl; + return *nEntries; +} + /////////////////////////////////////////////// /// \brief This function will add a new column to the RDataFrame using /// the same scheme as the usual RDF::Define method, but it will on top of @@ -539,7 +616,7 @@ ROOT::RDF::RNode TRestDataSet::MakeCut(const TRestCut* cut) { /// \endcode /// ROOT::RDF::RNode TRestDataSet::DefineColumn(const std::string& columnName, const std::string& formula) { - auto df = fDataSet; + auto df = fDataFrame; std::string evalFormula = formula; for (auto const& [name, properties] : fQuantity) @@ -784,7 +861,7 @@ void TRestDataSet::InitFromConfigFile() { void TRestDataSet::Export(const std::string& filename, std::vector excludeColumns) { RESTInfo << "Exporting dataset" << RESTendl; - std::vector columns = fDataSet.GetColumnNames(); + std::vector columns = fDataFrame.GetColumnNames(); if (!excludeColumns.empty()) { columns.erase(std::remove_if(columns.begin(), columns.end(), [&excludeColumns](std::string elem) { @@ -796,10 +873,10 @@ void TRestDataSet::Export(const std::string& filename, std::vector RESTInfo << "Re-Generating snapshot." << RESTendl; std::string user = getenv("USER"); std::string fOutName = "/tmp/rest_output_" + user + ".root"; - fDataSet.Snapshot("AnalysisTree", fOutName, columns); + fDataFrame.Snapshot("AnalysisTree", fOutName, columns); RESTInfo << "Re-importing analysis tree." << RESTendl; - fDataSet = ROOT::RDataFrame("AnalysisTree", fOutName); + fDataFrame = ROOT::RDataFrame("AnalysisTree", fOutName); TFile* f = TFile::Open(fOutName.c_str()); fTree = (TChain*)f->Get("AnalysisTree"); @@ -811,7 +888,7 @@ void TRestDataSet::Export(const std::string& filename, std::vector RESTInfo << "Re-Generating snapshot." << RESTendl; std::string user = getenv("USER"); std::string fOutName = "/tmp/rest_output_" + user + ".root"; - fDataSet.Snapshot("AnalysisTree", fOutName); + fDataFrame.Snapshot("AnalysisTree", fOutName); TFile* f = TFile::Open(fOutName.c_str()); fTree = (TChain*)f->Get("AnalysisTree"); @@ -875,7 +952,7 @@ void TRestDataSet::Export(const std::string& filename, std::vector fprintf(f, "###\n"); fprintf(f, "### Data starts here\n"); - auto obsNames = fDataSet.GetColumnNames(); + auto obsNames = fDataFrame.GetColumnNames(); std::string obsListStr = ""; for (const auto& l : obsNames) { if (!obsListStr.empty()) obsListStr += ":"; @@ -903,7 +980,7 @@ void TRestDataSet::Export(const std::string& filename, std::vector return; } else if (TRestTools::GetFileNameExtension(filename) == "root") { - fDataSet.Snapshot("AnalysisTree", filename); + fDataFrame.Snapshot("AnalysisTree", filename); TFile* f = TFile::Open(filename.c_str(), "UPDATE"); std::string name = this->GetName(); @@ -1003,7 +1080,7 @@ void TRestDataSet::Import(const std::string& fileName) { else ROOT::DisableImplicitMT(); - fDataSet = ROOT::RDataFrame("AnalysisTree", fileName); + fDataFrame = ROOT::RDataFrame("AnalysisTree", fileName); fTree = (TChain*)file->Get("AnalysisTree"); } @@ -1069,7 +1146,7 @@ void TRestDataSet::Import(std::vector fileNames) { } RESTInfo << "Opening list of files. First file: " << fileNames[0] << RESTendl; - fDataSet = ROOT::RDataFrame("AnalysisTree", fileNames); + fDataFrame = ROOT::RDataFrame("AnalysisTree", fileNames); if (fTree != nullptr) { delete fTree; diff --git a/source/framework/core/src/TRestDataSetPlot.cxx b/source/framework/core/src/TRestDataSetPlot.cxx index ec25cf686..d65e118aa 100644 --- a/source/framework/core/src/TRestDataSetPlot.cxx +++ b/source/framework/core/src/TRestDataSetPlot.cxx @@ -91,8 +91,12 @@ /// * **precision**: Precision for the values to be written. /// * **value**: If true/ON panel is displayed, otherwise is ignored /// Different keys are provided: `metadata` is meant for the metadata info inside the -/// TRestDataSet, `variable` for a predefined variable e.g. rate and `observable` for an -/// observable value. All the keys have the same structure which is detailed below: +/// TRestDataSet (as a RelevantQuantity), `variable` for a predefined variable e.g. rate, +/// `observable` for an observable value and `expression` for a mathematical expression +/// that can contain any of the previous. Note that the time-related variables _startTime_, +/// _endTime_ and _runLength_ (then _meanRate_ too) are obtained from the TRestDataSet and +/// not the RDataFrame of the TRestDataSet, thus they are not afected by the cuts. +/// All the keys have the same structure which is detailed below: /// * **value**: Name of the metadata, variable or observable value. /// * **label**: String of the label that will be written before the observable value. /// * **units**: String with the units to be appended to the label. @@ -107,6 +111,10 @@ /// /// /// +/// +/// /// /// /// \endcode @@ -124,8 +132,10 @@ /// * **timeDisplay**: If true/ON time display is set in the X axis /// * **norm**: Normalization constant in which the plot will be normalized, e.g. use `1` in /// case you want to normalize by 1. -/// * **scale**: Multiply all the histogram bins by a constant given by scale. If you want to -/// use the size of the bin to normalize you should write down `binSize` +/// * **scale**: Divide all the histogram bins by a constant given by scale. You may use any +/// mathematical expression in combination with the special keywords: `binSize`, `entries`, +/// `runLength` (in hours) and `integral`. Adding a scale will make the histogram to be +/// drawn with errors. To avoid this, set parameter option to 'hist' in the histogram. /// * **value**: If true/ON plot is displayed, otherwise is ignored /// * **save**: String with the name of the output file in which the plot will be saved /// in a separated file, several formats are supported (root, pdf, eps, jpg,...) @@ -252,6 +262,8 @@ /// 2023-04: First implementation of TRestDataSetPlot, based on TRestAnalysisPlot /// JuanAn Garcia /// +/// 2024-05: Extend some functionalities, Álvaro Ezquerro +/// /// \class TRestDataSetPlot /// \author: JuanAn Garcia e-mail: juanangp@unizar.es /// @@ -389,6 +401,19 @@ void TRestDataSetPlot::ReadPanelInfo() { observable = GetNextElement(observable); } + TiXmlElement* expression = GetElement("expression", panelele); + while (expression != nullptr) { + std::array label; + label[0] = GetParameter("value", expression, ""); + label[1] = GetParameter("label", expression, ""); + label[2] = GetParameter("units", expression, ""); + double posX = StringToDouble(GetParameter("x", expression, "0.1")); + double posY = StringToDouble(GetParameter("y", expression, "0.1")); + + panel.expPos.push_back(std::make_pair(label, TVector2(posX, posY))); + + expression = GetNextElement(expression); + } fPanels.push_back(panel); panelele = GetNextElement(panelele); @@ -505,7 +530,50 @@ void TRestDataSetPlot::GenerateDataSetFromFilePattern(TRestDataSet& dataSet) { std::map quantity; for (auto& panel : fPanels) { - // Add obserbables from panel info, both variables and cuts + // Add metadata and observables from expression key from panel info + for (auto& [key, posLabel] : panel.expPos) { + // look for metadata which are surrounded by [] but not [[]] (variables) + auto&& [exp, label, units] = key; + std::string text = exp; + while (text.find_last_of('[') != std::string::npos) { + int squareBracketCorrector = 0; + size_t posOpen = text.find_last_of('['); + size_t posClose = text.find_first_of(']', posOpen); + if (posOpen != 0) { + if (text[posOpen - 1] == '[') { + squareBracketCorrector = 1; + } + } + std::string varOrMeta = text.substr(posOpen - squareBracketCorrector, + posClose + 1 - posOpen + 2 * squareBracketCorrector); + if (squareBracketCorrector == 0) { + // Here varOrMeta is a metadata + // Add metadata to relevant quantity from the panel info + TRestDataSet::RelevantQuantity quant; + quant.metadata = varOrMeta; + quant.strategy = "unique"; + quantity[label] = quant; + } + text = Replace(text, varOrMeta, "1"); + } + + // look for observables (characterized by having a _ in the name) + while (text.find("_") != std::string::npos) { + size_t pos_ = text.find("_"); + size_t beginning = text.find_last_of(" -+*/)(^%", pos_) + 1; + size_t end = text.find_first_of(" -+*/)(^%", pos_); + std::string obs = text.substr(beginning, end - beginning); + text = Replace(text, obs, "1"); + obsList.push_back(obs); + } + + if (!(isAExpression(text) || isANumber(text))) + RESTWarning << "The expression " << exp + << " has not been correctly parsed into variables, metadata and observables" + << RESTendl; + } + + // Add observables from observable key of panel info for (auto& [key, posLabel] : panel.obsPos) { auto&& [obs, label, units] = key; obsList.push_back(obs); @@ -535,6 +603,7 @@ void TRestDataSetPlot::GenerateDataSetFromFilePattern(TRestDataSet& dataSet) { /// void TRestDataSetPlot::PlotCombinedCanvas() { TRestDataSet dataSet; + dataSet.EnableMultiThreading(true); // Import dataSet dataSet.Import(fDataSetName); @@ -586,6 +655,41 @@ void TRestDataSetPlot::PlotCombinedCanvas() { paramMap["[[runLength]]"] = StringWithPrecision(runLength, panel.precision); paramMap["[[entries]]"] = StringWithPrecision(entries, panel.precision); paramMap["[[meanRate]]"] = StringWithPrecision(meanRate, panel.precision); + + paramMap["[[cutNames]]"] = ""; + paramMap["[[cuts]]"] = ""; + if (fCut) { + for (const auto& cut : fCut->GetCuts()) { + if (paramMap["[[cutNames]]"].empty()) + paramMap["[[cutNames]]"] += cut.GetName(); + else + paramMap["[[cutNames]]"] += "," + (std::string)cut.GetName(); + if (paramMap["[[cuts]]"].empty()) + paramMap["[[cuts]]"] += cut.GetTitle(); + else + paramMap["[[cuts]]"] += " && " + (std::string)cut.GetTitle(); + } + } + + paramMap["[[panelCutNames]]"] = ""; + paramMap["[[panelCuts]]"] = ""; + if (panel.panelCut) { + for (const auto& cut : panel.panelCut->GetCuts()) { + if (paramMap["[[panelCutNames]]"].empty()) + paramMap["[[panelCutNames]]"] += cut.GetName(); + else + paramMap["[[panelCutNames]]"] += "," + (std::string)cut.GetName(); + if (paramMap["[[panelCuts]]"].empty()) + paramMap["[[panelCuts]]"] += cut.GetTitle(); + else + paramMap["[[panelCuts]]"] += " && " + (std::string)cut.GetTitle(); + } + } + + RESTInfo << "Global cuts: " << paramMap["[[cuts]]"] << RESTendl; + if (!paramMap["[[panelCuts]]"].empty()) + RESTInfo << "Additional panel cuts: " << paramMap["[[panelCuts]]"] << RESTendl; + // Replace panel variables and generate a TLatex label for (const auto& [key, posLabel] : panel.variablePos) { auto&& [variable, label, units] = key; @@ -633,6 +737,35 @@ void TRestDataSetPlot::PlotCombinedCanvas() { panel.text.emplace_back(new TLatex(posLabel.X(), posLabel.Y(), lab.c_str())); } + // Replace any expression and generate TLatex label + for (const auto& [key, posLabel] : panel.expPos) { + auto&& [text, label, units] = key; + std::string var = text; + + // replace variables + for (const auto& [param, val] : paramMap) { + var = Replace(var, param, val); + } + // replace metadata + for (const auto& [name, quant] : quantity) { + var = Replace(var, name, quant.value); + } + // replace observables + for (const auto& obs : dataFrame.GetColumnNames()) { + if (var.find(obs) == std::string::npos) continue; + // here there should be a checking that the mean(obs) can be calculated + // (checking obs data type?) + double value = *dataFrame.Mean(obs); + var = Replace(var, obs, DoubleToString(value)); + } + var = Replace(var, "[", "("); + var = Replace(var, "]", ")"); + var = EvaluateExpression(var); + + std::string lab = label + ": " + var + " " + units; + panel.text.emplace_back(new TLatex(posLabel.X(), posLabel.Y(), lab.c_str())); + } + // Draw the labels inside the pad for (const auto& text : panel.text) { text->SetTextColor(1); @@ -691,13 +824,19 @@ void TRestDataSetPlot::PlotCombinedCanvas() { } // Scale histos if (plots.scale != "") { - Double_t scale = 1.; - if (plots.scale == "binSize") { - scale = 1. / hist.histo->GetXaxis()->GetBinWidth(1); - } else { - scale = StringToDouble(plots.scale); - } - hist.histo->Scale(scale); + std::string inputScale = plots.scale; + double binSize = hist.histo->GetXaxis()->GetBinWidth(1); + double entries = hist.histo->GetEntries(); + double runLength = dataSet.GetTotalTimeInSeconds() / 3600.; // in hours + double integral = hist.histo->Integral("width"); + + inputScale = Replace(inputScale, "binSize", DoubleToString(binSize)); + inputScale = Replace(inputScale, "entries", DoubleToString(entries)); + inputScale = Replace(inputScale, "runLength", DoubleToString(runLength)); + inputScale = Replace(inputScale, "integral", DoubleToString(integral)); + + std::string scale = "1./(" + inputScale + ")"; + hist.histo->Scale(StringToDouble(EvaluateExpression(scale))); // -1 if 'scale' isn't valid } // Add histos to the THStack @@ -881,6 +1020,12 @@ void TRestDataSetPlot::PrintMetadata() { << " Pos (" << posLabel.X() << ", " << posLabel.Y() << ")" << RESTendl; } RESTMetadata << "****************" << RESTendl; + for (auto& [key, posLabel] : panel.expPos) { + auto&& [obs, label, units] = key; + RESTMetadata << "Label Expression " << obs << ", label " << label << ", units " << units + << " Pos (" << posLabel.X() << ", " << posLabel.Y() << ")" << RESTendl; + } + RESTMetadata << "****************" << RESTendl; } RESTMetadata << "-------------------" << RESTendl; diff --git a/source/framework/core/src/TRestEvent.cxx b/source/framework/core/src/TRestEvent.cxx index 7acb1831f..191dda2f6 100644 --- a/source/framework/core/src/TRestEvent.cxx +++ b/source/framework/core/src/TRestEvent.cxx @@ -49,6 +49,8 @@ #include "TRestEvent.h" +#include "TRestRun.h" + using namespace std; ClassImp(TRestEvent); @@ -170,7 +172,11 @@ void TRestEvent::RestartPad(Int_t nElements) { void TRestEvent::InitializeWithMetadata(TRestRun* run) { Initialize(); } -void TRestEvent::InitializeReferences(TRestRun* run) { fRun = run; } +void TRestEvent::InitializeReferences(TRestRun* run) { + fRun = run; + SetRunOrigin(fRun->GetRunNumber()); + SetSubRunOrigin(fRun->GetSubRunNumber()); +} ////////////////////////////////////////////////////////////////////////// /// Run to print event data info on console diff --git a/source/framework/core/src/TRestMetadata.cxx b/source/framework/core/src/TRestMetadata.cxx index e5e598b41..8a416e438 100644 --- a/source/framework/core/src/TRestMetadata.cxx +++ b/source/framework/core/src/TRestMetadata.cxx @@ -1284,7 +1284,7 @@ void TRestMetadata::ExpandIncludeFile(TiXmlElement* e) { url = _filename; } - filename = TRestTools::DownloadRemoteFile(url); + filename = TRestTools::DownloadRemoteFile(url, true); } else { filename = SearchFile(_filename); } @@ -2702,3 +2702,18 @@ void TRestMetadata::Merge(const TRestMetadata& metadata) { fName = metadata.GetName(); } } + +UInt_t TRestMetadata::GetVersionMajor() const { + TString major = fVersion(0, fVersion.First('.')); + return major.Atoi(); +} + +UInt_t TRestMetadata::GetVersionMinor() const { + TString minor = fVersion(fVersion.First('.') + 1, fVersion.Last('.')); + return minor.Atoi(); +} + +UInt_t TRestMetadata::GetVersionPatch() const { + TString patch = fVersion(fVersion.Last('.') + 1, fVersion.Length()); + return patch.Atoi(); +} diff --git a/source/framework/masks/inc/TRestPatternMask.h b/source/framework/masks/inc/TRestPatternMask.h index 12979a560..c16c6a89a 100644 --- a/source/framework/masks/inc/TRestPatternMask.h +++ b/source/framework/masks/inc/TRestPatternMask.h @@ -29,7 +29,7 @@ /// An abstract class used to encapsulate different mask pattern class definitions. class TRestPatternMask : public TRestMetadata { private: - /// It is used to introduce an offset on the pattern + /// It is used to introduce an offset on the pattern (not the mask, mask is always centered) TVector2 fOffset = TVector2(0, 0); //< /// An angle (in radians) used to introduce a rotation to the pattern diff --git a/source/framework/masks/inc/TRestRadialStrippedMask.h b/source/framework/masks/inc/TRestRadialStrippedMask.h new file mode 100644 index 000000000..b8d95caf2 --- /dev/null +++ b/source/framework/masks/inc/TRestRadialStrippedMask.h @@ -0,0 +1,70 @@ +/************************************************************************* + * 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_TRestRadialStrippedMask +#define REST_TRestRadialStrippedMask + +#include + +/// A class used to define a stripped mask pattern +class TRestRadialStrippedMask : public TRestPatternMask { + private: + void Initialize() override; + + /// The periodity of the stripped structure in radians + Double_t fStripsAngle = TMath::Pi() / 3; //< + + /// The width of the stripped structure in mm + Double_t fStripsThickness = 0.5; //< + + /// The spacers structure will be effective from this radius, in mm. Default is from 20 mm. + Double_t fInitialRadius = 20.; //< + + /// Radius of an internal circular region defined inside the fInitialRadius. If 0, there will be no region + Double_t fInternalRegionRadius = 0.; //< + + /// It defines the maximum number of cells/regions in each axis + Int_t fModulus = 10; + + public: + virtual Int_t GetRegion(Double_t& x, Double_t& y) override; + + /// It returns the gap/periodicity of the strips in degrees + Double_t GetStripsAngle() { return fStripsAngle * units("degrees"); } + + /// It returns the thickness of the strips in mm + Double_t GetStripsThickness() { return fStripsThickness; } + + /// It returns the modulus used to define a finite set of ids + Int_t GetModulus() { return fModulus; } + + void PrintMetadata() override; + void PrintMaskMembers() override; + void PrintMask() override; + + TRestRadialStrippedMask(); + TRestRadialStrippedMask(const char* cfgFileName, std::string name = ""); + ~TRestRadialStrippedMask(); + + ClassDefOverride(TRestRadialStrippedMask, 1); +}; +#endif diff --git a/source/framework/masks/src/TRestPatternMask.cxx b/source/framework/masks/src/TRestPatternMask.cxx index 74fa9143c..5388888b3 100644 --- a/source/framework/masks/src/TRestPatternMask.cxx +++ b/source/framework/masks/src/TRestPatternMask.cxx @@ -167,7 +167,7 @@ TCanvas* TRestPatternMask::DrawMonteCarlo(Int_t nSamples) { delete fCanvas; fCanvas = NULL; } - fCanvas = new TCanvas("canv", "This is the canvas title", 1400, 1200); + fCanvas = new TCanvas("canv", "This is the canvas title", 1200, 1200); fCanvas->Draw(); TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97); @@ -183,17 +183,19 @@ TCanvas* TRestPatternMask::DrawMonteCarlo(Int_t nSamples) { TRandom3* rnd = new TRandom3(0); for (int n = 0; n < nSamples; n++) { - Double_t x = 2.5 * (rnd->Rndm() - 0.5) * fMaskRadius; - Double_t y = 2.5 * (rnd->Rndm() - 0.5) * fMaskRadius; + Double_t xO = 2.5 * (rnd->Rndm() - 0.5) * fMaskRadius; + Double_t yO = 2.5 * (rnd->Rndm() - 0.5) * fMaskRadius; + Double_t x = xO; + Double_t y = yO; Int_t id = GetRegion(x, y); if (points.count(id) == 0) { std::vector a; - a.push_back(TVector2(x, y)); + a.push_back(TVector2(xO, yO)); points[id] = a; } else { - points[id].push_back(TVector2(x, y)); + points[id].push_back(TVector2(xO, yO)); } } diff --git a/source/framework/masks/src/TRestRadialStrippedMask.cxx b/source/framework/masks/src/TRestRadialStrippedMask.cxx new file mode 100644 index 000000000..6985e92c3 --- /dev/null +++ b/source/framework/masks/src/TRestRadialStrippedMask.cxx @@ -0,0 +1,221 @@ +/************************************************************************* + * 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 defines a stripped pattern. It defines a periodicity +/// and a thickness for the strips. The method TRestRadialStrippedMask::GetRegion +/// will return a unique id for each region in between strips. +/// +/// The stripped structure is centered in (0,0) and it can be shifted using +/// the offset defined inside TRestPatternMask. The pattern will be only +/// delimited by the limits imposed inside TRestPatternMask. +/// +/// ### Specific stripped metadata parameters +/// +/// * **stripsAngle**: This parameter defines the strips angular periodicity. +/// * **stripsThickness**: The thickness of the strips. +/// * **modulus**: A number that defines the range of ids used to identify +/// the different regions inside the stripped pattern. If modulus is 10, +/// then we will only be able to identify up to 10 unique regions. If a +/// larger amount of regions is found, it will happen that two regions will +/// be assigned the same id. +/// +/// ### Common pattern metadata parameters +/// +/// On top of the metadata class parameters, we may define common pattern +/// parameters to induce an offset and rotation to the pattern. +/// +/// * **offset**: A parameter to shift the pattern window mask. +/// * **rotationAngle**: An angle given in radians to rotate the pattern. +/// * **maskRadius**: A radius defining the limits of the circular mask. +/// +/// ### Examples +/// +/// Mask pattern RML definitions can be found inside the file +/// `REST_PATH/examples/masks.rml`. +/// +/// The following definition ilustrates a complete RML implementation of a +/// TRestRadialStrippedMask. +/// +/// \code +/// +/// +/// +/// +/// +/// +/// +/// +/// \endcode +/// +/// The basic use of this class is provided by the TRestRadialStrippedMask::GetRegion +/// method. For example: +/// +/// \code +/// TRestRadialStrippedMask mask("masks.rml", "radialStrips"); +/// Int_t id = mask.GetRegion( 12.5, 4.3 ); +/// std::cout << "Region id is : " << id << endl; +/// \endcode +/// +/// The following figure may be generated using the TRestPatternMask::DrawMonteCarlo +/// method. +/// +/// \code +/// TRestRadialStrippedMask mask("masks.rml", "radialStrips"); +/// TCanvas *c = mask.DrawMonteCarlo(30000); +/// c->Draw(); +/// c->Print("radialstrippedmask.png"); +/// \endcode +/// +/// \htmlonly \endhtmlonly +/// ![An illustration of the montecarlo mask test using DrawMonteCarlo](strippedmask.png) +/// +///---------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2022-05: First implementation of TRestRadialStrippedMask +/// Javier Galan +/// +/// \class TRestRadialStrippedMask +/// \author: Javier Galan - javier.galan@unizar.es +/// +///
+/// + +#include "TRestRadialStrippedMask.h" + +#include "TRandom3.h" + +ClassImp(TRestRadialStrippedMask); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestRadialStrippedMask::TRestRadialStrippedMask() : TRestPatternMask() { 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 TRestRadialStrippedMask section inside the RML. +/// +TRestRadialStrippedMask::TRestRadialStrippedMask(const char* cfgFileName, std::string name) + : TRestPatternMask(cfgFileName) { + Initialize(); + + LoadConfigFromFile(fConfigFileName, name); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestRadialStrippedMask::~TRestRadialStrippedMask() {} + +/////////////////////////////////////////////// +/// \brief Function to initialize input/output event members and define +/// the section name +/// +void TRestRadialStrippedMask::Initialize() { + SetSectionName(this->ClassName()); + SetType("RadialStripped"); +} + +/////////////////////////////////////////////// +/// \brief It returns a number identifying the region where the particle +/// with coordinates (x,y) felt in. The method returns 0 if the particle +/// hits the pattern. +/// +/// The particle will be counter-rotated to emulate the mask rotation +/// using the method TRestPatternMask::ApplyCommonMaskTransformation +/// +Int_t TRestRadialStrippedMask::GetRegion(Double_t& x, Double_t& y) { + if (TRestPatternMask::GetRegion(x, y)) return 0; + + Double_t d = TMath::Sqrt(x * x + y * y); + + if (d < fInitialRadius) { + if (fInternalRegionRadius > 0 && d < fInternalRegionRadius) return 1; + + return 0; + } + + TVector2 point(x, y); + Double_t phi = point.Phi(); + + /// phi determines the region where the point is found + Int_t region = (Int_t)(phi / fStripsAngle); + region = 2 + region % fMaxRegions; + + Double_t angle = 0; + /// Checking if we hit an arm + while (angle < 2 * TMath::Pi()) { + if (point.Y() < fStripsThickness / 2. && point.Y() > -fStripsThickness / 2. && point.X() >= 0) + return 0; + + point = point.Rotate(fStripsAngle); + angle += fStripsAngle; + } + + return 1 + region % fModulus; +} + +///////////////////////////////////////////// +/// \brief Prints on screen the complete information about the metadata members from this class +/// +void TRestRadialStrippedMask::PrintMetadata() { + TRestPatternMask::PrintMetadata(); + + PrintMaskMembers(); + RESTMetadata << "++++" << RESTendl; +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestRingsMask, +/// including common pattern headers, but without common metadata headers. +/// +void TRestRadialStrippedMask::PrintMask() { + PrintCommonPatternMembers(); + RESTMetadata << "----" << RESTendl; + PrintMaskMembers(); +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestRingsMask, +/// excluding common metadata headers. +/// +void TRestRadialStrippedMask::PrintMaskMembers() { + RESTMetadata << " - Strips angle : " << fStripsAngle * units("degrees") << " degrees" << RESTendl; + RESTMetadata << " - Strips thickness : " << fStripsThickness << " mm" << RESTendl; +} diff --git a/source/framework/sensitivity/inc/TRestComponent.h b/source/framework/sensitivity/inc/TRestComponent.h index 1dfd614dc..1e22722a6 100644 --- a/source/framework/sensitivity/inc/TRestComponent.h +++ b/source/framework/sensitivity/inc/TRestComponent.h @@ -55,6 +55,18 @@ class TRestComponent : public TRestMetadata { /// It defines the nodes of the parameterization (Initialized by the dataset) std::vector fParameterizationNodes; //< + /// It defines the first parametric node value in case of automatic parameter generation + Double_t fFirstParameterValue = 0; //< + + /// It defines the upper limit for the automatic parametric node values generation + Double_t fLastParameterValue = 0; //< + + /// It defines the increasing step for automatic parameter list generation + Double_t fStepParameterValue = 0; //< + + /// It true the parametric values automatically generated will grow exponentially + Bool_t fExponential = false; //< + /// It is used to define the node that will be accessed for rate retrieval Int_t fActiveNode = -1; //< @@ -82,9 +94,6 @@ class TRestComponent : public TRestMetadata { /// 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 @@ -104,6 +113,11 @@ class TRestComponent : public TRestMetadata { void Initialize() override; void RegenerateHistograms(UInt_t seed = 0); + void RegenerateParametricNodes(Double_t from, Double_t to, Double_t step, Bool_t expIncrease = false); + + /// It returns true if any nodes have been defined. + Bool_t HasNodes() { return !fParameterizationNodes.empty(); } + virtual void RegenerateActiveNodeDensity() {} std::string GetNature() const { return fNature; } @@ -112,7 +126,12 @@ class TRestComponent : public TRestMetadata { size_t GetDimensions() { return fVariables.size(); } Int_t GetSamples() { return fSamples; } Int_t GetActiveNode() { return fActiveNode; } - Double_t GetActiveNodeValue() { return fParameterizationNodes[fActiveNode]; } + Double_t GetActiveNodeValue() { + if (fActiveNode >= 0 && fActiveNode < (Int_t)fParameterizationNodes.size()) + return fParameterizationNodes[fActiveNode]; + return 0; + } + std::vector GetParameterizationNodes() { return fParameterizationNodes; } std::vector GetVariables() const { return fVariables; } std::vector GetRanges() const { return fRanges; } @@ -120,6 +139,8 @@ class TRestComponent : public TRestMetadata { Double_t GetRawRate(std::vector point); Double_t GetTotalRate(); + Double_t GetMaxRate(); + Double_t GetAllNodesIntegratedRate(); Double_t GetNormalizedRate(std::vector point); Double_t GetRate(std::vector point); @@ -127,6 +148,7 @@ class TRestComponent : public TRestMetadata { void SetPrecision(const Float_t& pr) { fPrecision = pr; } + Int_t FindActiveNode(Double_t node); Int_t SetActiveNode(Double_t node); Int_t SetActiveNode(Int_t n) { fActiveNode = n; @@ -169,6 +191,6 @@ class TRestComponent : public TRestMetadata { TRestComponent(); ~TRestComponent(); - ClassDefOverride(TRestComponent, 5); + ClassDefOverride(TRestComponent, 6); }; #endif diff --git a/source/framework/sensitivity/inc/TRestComponentDataSet.h b/source/framework/sensitivity/inc/TRestComponentDataSet.h index 9a2154f05..59e276268 100644 --- a/source/framework/sensitivity/inc/TRestComponentDataSet.h +++ b/source/framework/sensitivity/inc/TRestComponentDataSet.h @@ -54,6 +54,12 @@ class TRestComponentDataSet : public TRestComponent { /// The dataset used to initialize the distribution TRestDataSet fDataSet; //! + /// It helps to split large datasets when extracting the parameterization nodes + long long unsigned int fSplitEntries = 600000000; + + /// It creates a sample subset using a range definition + TVector2 fDFRange = TVector2(0, 0); + /// It is true of the dataset was loaded without issues Bool_t fDataSetLoaded = false; //! @@ -84,6 +90,6 @@ class TRestComponentDataSet : public TRestComponent { TRestComponentDataSet(const char* cfgFileName, const std::string& name); ~TRestComponentDataSet(); - ClassDefOverride(TRestComponentDataSet, 3); + ClassDefOverride(TRestComponentDataSet, 4); }; #endif diff --git a/source/framework/sensitivity/inc/TRestComponentFormula.h b/source/framework/sensitivity/inc/TRestComponentFormula.h index 9f9046824..1913c5f3f 100644 --- a/source/framework/sensitivity/inc/TRestComponentFormula.h +++ b/source/framework/sensitivity/inc/TRestComponentFormula.h @@ -35,7 +35,7 @@ class TRestComponentFormula : public TRestComponent { std::vector fFormulas; /// The formulas should be expressed in the following units - std::string fUnits = "cm^-2*keV^-1"; //< + std::string fFormulaUnits = "cm^-2*keV^-1"; //< protected: void InitFromConfigFile() override; diff --git a/source/framework/sensitivity/inc/TRestExperiment.h b/source/framework/sensitivity/inc/TRestExperiment.h index 7e25014f2..440ec3512 100644 --- a/source/framework/sensitivity/inc/TRestExperiment.h +++ b/source/framework/sensitivity/inc/TRestExperiment.h @@ -42,7 +42,7 @@ class TRestExperiment : public TRestMetadata { TRestComponent* fSignal = nullptr; //< /// It defines the filename used to load the dataset - std::string fExperimentalDataSet = ""; + std::string fExperimentalDataSet = ""; //< /// It contains the experimental data (should contain same columns as the components) TRestDataSet fExperimentalData; //< @@ -53,6 +53,12 @@ class TRestExperiment : public TRestMetadata { /// Only if it is true we will be able to calculate the LogLikelihood Bool_t fDataReady = false; //< + /// The mock dataset will be generated using the mean counts instead of a real MonteCarlo + Bool_t fUseAverage = false; //< + + /// It keeps track on the number of counts inside the dataset + Int_t fExperimentalCounts = 0; //< + /// Internal process random generator TRandom3* fRandom = nullptr; //! @@ -63,7 +69,8 @@ class TRestExperiment : public TRestMetadata { void InitFromConfigFile() override; public: - void GenerateMockDataSet(); + void GenerateMockDataSet(Bool_t useAverage = false); + Int_t GetExperimentalCounts() const { return fExperimentalCounts; } Bool_t IsMockData() const { return fMockData; } Bool_t IsDataReady() const { return fDataReady; } @@ -90,6 +97,6 @@ class TRestExperiment : public TRestMetadata { TRestExperiment(); ~TRestExperiment(); - ClassDefOverride(TRestExperiment, 1); + ClassDefOverride(TRestExperiment, 2); }; #endif diff --git a/source/framework/sensitivity/inc/TRestExperimentList.h b/source/framework/sensitivity/inc/TRestExperimentList.h index 5d1737e40..ca0b60e02 100644 --- a/source/framework/sensitivity/inc/TRestExperimentList.h +++ b/source/framework/sensitivity/inc/TRestExperimentList.h @@ -50,9 +50,21 @@ class TRestExperimentList : public TRestMetadata { /// A vector with a list of experiments includes the background components in this model std::vector fExperiments; //< + /// The fusioned list of parameterization nodes found at each experiment signal + std::vector fParameterizationNodes; //< + + /// The mock dataset will be generated using the mean counts instead of a real MonteCarlo + Bool_t fUseAverage = false; //< + /// If not zero this will be the common exposure time in micro-seconds (standard REST units) Double_t fExposureTime = 0; + /// In case an exposure time is given it defines how to assign time to each experiment (equal/ksvz). + std::string fExposureStrategy = "equal"; + + /// The factor used on the exponential exposure time as a function of the experiment number + Double_t fExposureFactor = 0; + /// If not null this will be the common signal used in each experiment TRestComponent* fSignal = nullptr; //< @@ -71,6 +83,10 @@ class TRestExperimentList : public TRestMetadata { void SetSignal(TRestComponent* comp) { fSignal = comp; } void SetBackground(TRestComponent* comp) { fBackground = comp; } + void ExtractExperimentParameterizationNodes(); + std::vector GetParameterizationNodes() { return fParameterizationNodes; } + void PrintParameterizationNodes(); + std::vector GetExperiments() { return fExperiments; } TRestExperiment* GetExperiment(const size_t& n) { if (n >= GetNumberOfExperiments()) @@ -88,6 +104,6 @@ class TRestExperimentList : public TRestMetadata { TRestExperimentList(); ~TRestExperimentList(); - ClassDefOverride(TRestExperimentList, 1); + ClassDefOverride(TRestExperimentList, 2); }; #endif diff --git a/source/framework/sensitivity/inc/TRestSensitivity.h b/source/framework/sensitivity/inc/TRestSensitivity.h index 56c9d0bd6..c60350f9e 100644 --- a/source/framework/sensitivity/inc/TRestSensitivity.h +++ b/source/framework/sensitivity/inc/TRestSensitivity.h @@ -29,22 +29,52 @@ class TRestSensitivity : public TRestMetadata { private: /// A list of experimental conditions included to get a final sensitivity plot - std::vector fExperiments; //< + std::vector fExperiments; //! - TH1D* fSignalTest = nullptr; + /// The fusioned list of parameterization nodes found at each experiment signal + std::vector fParameterizationNodes; //< + + /// A vector of calculated sensitivity curves defined as a funtion of the parametric node + std::vector> fCurves; //< + + /// A flag that will frozen adding more experiments in the future. + Bool_t fFrozen = false; //< Only needed if we add experiments by other means than RML + + /// It is used to generate a histogram with the signal distribution produced with different signal samples + TH1D* fSignalTest = nullptr; //< + + /// A canvas to draw + TCanvas* fCanvas = 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); + Double_t UnbinnedLogLikelihood(const TRestExperiment* experiment, Double_t node, Double_t g4 = 0); + Double_t ApproachByFactor(Double_t node, 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); + void ExtractExperimentParameterizationNodes(Bool_t rescan = false); + std::vector GetParameterizationNodes() { return fParameterizationNodes; } + void PrintParameterizationNodes(); + + Double_t GetCoupling(Double_t node, Double_t sigma = 2, Double_t precision = 0.01); + void AddCurve(const std::vector& curve) { fCurves.push_back(curve); } + void ImportCurve(const std::vector& curve) { AddCurve(curve); } + void GenerateCurve(); + void GenerateCurves(Int_t N); - TH1D* SignalStatisticalTest(Int_t N); + std::vector GetCurve(size_t n = 0); + std::vector GetAveragedCurve(); + std::vector> GetLevelCurves(const std::vector& levels); + + void ExportCurve(std::string fname, Double_t factor = 1.e-10, Double_t power = 0.25, int n = 0); + void ExportAveragedCurve(std::string fname, Double_t factor = 1.e-10, Double_t power = 0.25); + + TH1D* SignalStatisticalTest(Double_t node, Int_t N); + + void Freeze() { fFrozen = true; } std::vector GetExperiments() { return fExperiments; } TRestExperiment* GetExperiment(const size_t& n) { @@ -55,13 +85,18 @@ class TRestSensitivity : public TRestMetadata { } size_t GetNumberOfExperiments() { return fExperiments.size(); } + size_t GetNumberOfCurves() { return fCurves.size(); } + size_t GetNumberOfNodes() { return fParameterizationNodes.size(); } void PrintMetadata() override; + TCanvas* DrawCurves(); + TCanvas* DrawLevelCurves(); + TRestSensitivity(const char* cfgFileName, const std::string& name = ""); TRestSensitivity(); ~TRestSensitivity(); - ClassDefOverride(TRestSensitivity, 1); + ClassDefOverride(TRestSensitivity, 2); }; #endif diff --git a/source/framework/sensitivity/src/TRestComponent.cxx b/source/framework/sensitivity/src/TRestComponent.cxx index 516bf8104..c723bd374 100644 --- a/source/framework/sensitivity/src/TRestComponent.cxx +++ b/source/framework/sensitivity/src/TRestComponent.cxx @@ -86,6 +86,9 @@ TRestComponent::~TRestComponent() {} void TRestComponent::Initialize() { // SetSectionName(this->ClassName()); + /// Avoiding double initialization + if (!fNodeDensity.empty() && fRandom) return; + if (!fRandom) { delete fRandom; fRandom = nullptr; @@ -93,6 +96,13 @@ void TRestComponent::Initialize() { fRandom = new TRandom3(fSeed); fSeed = fRandom->TRandom::GetSeed(); + + if (fStepParameterValue > 0) { + RegenerateParametricNodes(fFirstParameterValue, fLastParameterValue, fStepParameterValue, + fExponential); + } else { + if (!fParameterizationNodes.empty()) FillHistograms(); + } } ///////////////////////////////////////////// @@ -105,11 +115,33 @@ void TRestComponent::RegenerateHistograms(UInt_t seed) { fNodeDensity.clear(); fSeed = seed; - TRestComponent::Initialize(); - FillHistograms(); } +///////////////////////////////////////////// +/// \brief It allows to produce a parameter nodes list providing the initial +/// value, the final value and the step. We might chose the step growing in +/// linear increase steps or exponential. Linear is the default value. +/// +void TRestComponent::RegenerateParametricNodes(Double_t from, Double_t to, Double_t step, + Bool_t expIncrease) { + fStepParameterValue = step; + fFirstParameterValue = from; + fLastParameterValue = to; + fExponential = expIncrease; + + fParameterizationNodes.clear(); + + if (expIncrease) { + for (double p = from; p < to; p *= step) fParameterizationNodes.push_back(p); + } else { + for (double p = from; p < to; p += step) fParameterizationNodes.push_back(p); + } + + if (fParameterizationNodes.empty()) return; + RegenerateHistograms(fSeed); +} + /////////////////////////////////////////// /// \brief It returns the position of the fVariable element for the variable /// name given by argument. @@ -232,6 +264,11 @@ Double_t TRestComponent::GetRawRate(std::vector point) { return 0; } + for (size_t n = 0; n < point.size(); n++) { + // The point is outside boundaries + if (point[n] < fRanges[n].X() || point[n] > fRanges[n].Y()) return 0; + } + Int_t centerBin[GetDimensions()]; Double_t centralDensity = GetDensity()->GetBinContent(GetDensity()->GetBin(point.data()), centerBin); if (!Interpolation()) return centralDensity; @@ -289,17 +326,53 @@ Double_t TRestComponent::GetRawRate(std::vector point) { /// Double_t TRestComponent::GetTotalRate() { THnD* dHist = GetDensityForActiveNode(); + if (!dHist) return 0; Double_t integral = 0; - if (dHist != nullptr) { - TH1D* h1 = dHist->Projection(0); - integral = h1->Integral(); - delete h1; + for (Int_t n = 0; n < dHist->GetNbins(); ++n) { + Int_t centerBin[GetDimensions()]; + std::vector point; + + dHist->GetBinContent(n, centerBin); + for (size_t d = 0; d < GetDimensions(); ++d) point.push_back(GetBinCenter(d, centerBin[d])); + + Bool_t skip = false; + for (size_t d = 0; d < GetDimensions(); ++d) { + if (point[d] < fRanges[d].X() || point[d] > fRanges[d].Y()) skip = true; + } + if (!skip) integral += GetRate(point); } return integral; } +/////////////////////////////////////////////// +/// \brief This method returns the total rate for the node that has the highest contribution +/// The result will be returned in s-1. +/// +Double_t TRestComponent::GetMaxRate() { + Double_t maxRate = 0; + for (size_t n = 0; n < fParameterizationNodes.size(); n++) { + SetActiveNode((Int_t)n); + Double_t rate = GetTotalRate(); + if (rate > maxRate) maxRate = rate; + } + return maxRate; +} + +/////////////////////////////////////////////// +/// \brief This method returns the integrated total rate for all the nodes +/// The result will be returned in s-1. +/// +Double_t TRestComponent::GetAllNodesIntegratedRate() { + Double_t rate = 0; + for (size_t n = 0; n < fParameterizationNodes.size(); n++) { + SetActiveNode((Int_t)n); + rate += GetTotalRate(); + } + return rate; +} + /////////////////////////////////////////////// /// \brief It returns the bin center of the given component dimension. /// @@ -561,7 +634,7 @@ void TRestComponent::PrintMetadata() { } } - if (!fParameter.empty()) { + if (!fParameterizationNodes.empty()) { RESTMetadata << " " << RESTendl; RESTMetadata << " === Parameterization === " << RESTendl; RESTMetadata << "- Parameter : " << fParameter << RESTendl; @@ -569,6 +642,18 @@ void TRestComponent::PrintMetadata() { RESTMetadata << " - Number of parametric nodes : " << fParameterizationNodes.size() << RESTendl; RESTMetadata << " " << RESTendl; RESTMetadata << " Use : PrintNodes() for additional info" << RESTendl; + + if (fStepParameterValue > 0) { + RESTMetadata << " " << RESTendl; + RESTMetadata << " Nodes were automatically generated using these parameters" << RESTendl; + RESTMetadata << " - First node : " << fFirstParameterValue << RESTendl; + RESTMetadata << " - Upper limit node : " << fLastParameterValue << RESTendl; + RESTMetadata << " - Increasing step : " << fStepParameterValue << RESTendl; + if (fExponential) + RESTMetadata << " - Increases exponentially" << RESTendl; + else + RESTMetadata << " - Increases linearly" << RESTendl; + } } if (fResponse) { @@ -627,6 +712,25 @@ void TRestComponent::InitFromConfigFile() { if (fResponse) fResponse->LoadResponse(); } +///////////////////////////////////////////// +/// \brief It returns the position of the fParameterizationNodes +/// element for the variable name given by argument. +/// +Int_t TRestComponent::FindActiveNode(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++; + } + + return -1; +} + ///////////////////////////////////////////// /// \brief It returns the position of the fParameterizationNodes /// element for the variable name given by argument. diff --git a/source/framework/sensitivity/src/TRestComponentDataSet.cxx b/source/framework/sensitivity/src/TRestComponentDataSet.cxx index 951018acc..94d1ef7bc 100644 --- a/source/framework/sensitivity/src/TRestComponentDataSet.cxx +++ b/source/framework/sensitivity/src/TRestComponentDataSet.cxx @@ -148,6 +148,11 @@ void TRestComponentDataSet::PrintMetadata() { RESTMetadata << " " << RESTendl; } + if (fDFRange.X() != 0 || fDFRange.Y() != 0) { + RESTMetadata << " DataFrame range: ( " << fDFRange.X() << ", " << fDFRange.Y() << ")" << RESTendl; + RESTMetadata << " " << RESTendl; + } + if (!fParameter.empty() && fParameterizationNodes.empty()) { RESTMetadata << "This component has no nodes!" << RESTendl; RESTMetadata << " Use: LoadDataSets() to initialize the nodes" << RESTendl; @@ -383,15 +388,17 @@ std::vector TRestComponentDataSet::ExtractParameterizationNodes() { return vs; } - auto parValues = fDataSet.GetDataFrame().Take(fParameter); - for (const auto v : parValues) vs.push_back(v); + auto GetUniqueElements = [](const std::vector& vec) { + std::set uniqueSet(vec.begin(), vec.end()); + return std::vector(uniqueSet.begin(), uniqueSet.end()); + }; - 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)); + for (size_t n = 0; n < 1 + fDataSet.GetEntries() / fSplitEntries; n++) { + auto nEn = fDataSet.Range(n * fSplitEntries, (n + 1) * fSplitEntries).Count(); + auto parValues = fDataSet.Range(n * fSplitEntries, (n + 1) * fSplitEntries).Take(fParameter); + std::vector uniqueVec = GetUniqueElements(*parValues); + vs.insert(vs.end(), uniqueVec.begin(), uniqueVec.end()); + } return vs; } @@ -476,6 +483,9 @@ Bool_t TRestComponentDataSet::LoadDataSets() { fDataSet.Import(fullFileNames); fDataSetLoaded = true; + if (fDFRange.X() != 0 || fDFRange.Y() != 0) + fDataSet.ApplyRange((size_t)fDFRange.X(), (size_t)fDFRange.Y()); + if (fDataSet.GetTree() == nullptr) { RESTError << "Problem loading dataset from file list :" << RESTendl; for (const auto& f : fDataSetFileNames) RESTError << " - " << f << RESTendl; @@ -486,6 +496,7 @@ Bool_t TRestComponentDataSet::LoadDataSets() { if (VariablesOk() && WeightsOk()) { fParameterizationNodes = ExtractParameterizationNodes(); + RESTInfo << "Filling histograms" << RESTendl; FillHistograms(); return fDataSetLoaded; } @@ -515,11 +526,12 @@ 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) { + for (const auto& var : fWeights) { + if (!isANumber(var) && std::count(cNames.begin(), cNames.end(), var) == 0) { RESTError << "Weight ---> " << var << " <--- NOT found on dataset" << RESTendl; ok = false; } + } return ok; } diff --git a/source/framework/sensitivity/src/TRestComponentFormula.cxx b/source/framework/sensitivity/src/TRestComponentFormula.cxx index 06c8ef3d6..db5d518fd 100644 --- a/source/framework/sensitivity/src/TRestComponentFormula.cxx +++ b/source/framework/sensitivity/src/TRestComponentFormula.cxx @@ -131,7 +131,7 @@ Double_t TRestComponentFormula::GetFormulaRate(std::vector point) { normFactor *= (fRanges[n].Y() - fRanges[n].X()) / fNbins[n]; } - return normFactor * result / units(fUnits); + return normFactor * result / units(fFormulaUnits); } ///////////////////////////////////////////// @@ -149,6 +149,12 @@ Double_t TRestComponentFormula::GetFormulaRate(std::vector point) { void TRestComponentFormula::FillHistograms() { if (fFormulas.empty()) return; + if (GetDimensions() == 0) { + RESTError << "TRestComponentFormula::FillHistograms. No variables defined!" << RESTendl; + RESTError << "Did you add a TRandom::GetSeed(); } -void TRestExperiment::GenerateMockDataSet() { +void TRestExperiment::GenerateMockDataSet(Bool_t useAverage) { + fUseAverage = useAverage; + if (!fBackground) { RESTError << "TRestExperiment::GenerateMockData. Background component was not initialized!" << RESTendl; @@ -102,12 +104,16 @@ void TRestExperiment::GenerateMockDataSet() { Double_t meanCounts = GetBackground()->GetTotalRate() * fExposureTime * units("s"); Int_t N = fRandom->Poisson(meanCounts); + if (fUseAverage) N = (Int_t)meanCounts; + RESTInfo << "Experiment: " << GetName() << " Generating mock dataset. Counts: " << N << RESTendl; ROOT::RDF::RNode df = fBackground->GetMonteCarloDataFrame(N); fExperimentalData.SetDataFrame(df); fExperimentalData.SetTotalTimeInSeconds(fExposureTime * units("s")); + fExperimentalCounts = *fExperimentalData.GetDataFrame().Count(); + fMockData = true; fDataReady = true; } @@ -118,6 +124,7 @@ void TRestExperiment::SetExperimentalDataSet(const std::string& filename) { /// fExposureTime is in standard REST units : us fExposureTime = fExperimentalData.GetTotalTimeInSeconds() / units("s"); + fExperimentalCounts = *fExperimentalData.GetDataFrame().Count(); fMockData = false; fDataReady = true; @@ -149,17 +156,19 @@ 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; - + TRestMetadata* md = (TRestMetadata*)this->InstantiateChildMetadata(cont); + while (md != nullptr) { + if (md->InheritsFrom("TRestComponent")) { + TRestComponent* comp = (TRestComponent*)md; + if (ToLower(comp->GetNature()) == "background") + fBackground = comp; + else if (ToLower(comp->GetNature()) == "signal") + fSignal = comp; + else + RESTWarning << "TRestExperiment::InitFromConfigFile. Unknown component!" << RESTendl; + } cont++; - comp = (TRestComponent*)this->InstantiateChildMetadata(cont, "Component"); + md = (TRestMetadata*)this->InstantiateChildMetadata(cont); } auto ele = GetElement("addComponent"); @@ -196,7 +205,7 @@ void TRestExperiment::InitFromConfigFile() { } if (fExposureTime > 0 && fExperimentalDataSet.empty()) { - GenerateMockDataSet(); + GenerateMockDataSet(fUseAverage); } else if (fExposureTime == 0 && !fExperimentalDataSet.empty()) { SetExperimentalDataSet(fExperimentalDataSet); } else { @@ -274,16 +283,21 @@ void TRestExperiment::PrintMetadata() { if (fMockData) { RESTMetadata << " " << RESTendl; - if (fMockData) + if (fMockData) { RESTMetadata << "The dataset was MC-generated" << RESTendl; - else { + if (fUseAverage) + RESTMetadata + << " - The number of counts in dataset was generated with the mean background counts" + << RESTendl; + + } else { RESTMetadata << "The dataset was loaded from an existing dataset file" << RESTendl; if (!fExperimentalDataSet.empty()) RESTMetadata << " - Experimental dataset file : " << fExperimentalDataSet << RESTendl; } } - RESTMetadata << " - Experimental counts : " << *fExperimentalData.GetDataFrame().Count() << RESTendl; + RESTMetadata << " - Experimental counts : " << fExperimentalCounts << RESTendl; RESTMetadata << "----" << RESTendl; } diff --git a/source/framework/sensitivity/src/TRestExperimentList.cxx b/source/framework/sensitivity/src/TRestExperimentList.cxx index ab773d024..b370fcb9e 100644 --- a/source/framework/sensitivity/src/TRestExperimentList.cxx +++ b/source/framework/sensitivity/src/TRestExperimentList.cxx @@ -177,14 +177,17 @@ void TRestExperimentList::InitFromConfigFile() { } column++; } else { - experiment->SetExposureInSeconds(fExposureTime * units("s")); - // We will generate mock data once we load the background component - generateMockData = true; + if (ToLower(fExposureStrategy) == "unique") { + experiment->SetExposureInSeconds(fExposureTime * units("s")); + // We will generate mock data once we load the background component + generateMockData = true; + } } if (!fSignal) { if (GetComponent(experimentRow[column])) { TRestComponent* sgnl = (TRestComponent*)GetComponent(experimentRow[column])->Clone(); + sgnl->SetName((TString)experimentRow[column]); experiment->SetSignal(sgnl); } else { RESTError << "TRestExperimentList. Signal component : " << experimentRow[column] @@ -209,14 +212,78 @@ void TRestExperimentList::InitFromConfigFile() { if (generateMockData) { RESTInfo << "TRestExperimentList. Generating mock dataset" << RESTendl; - experiment->GenerateMockDataSet(); + experiment->GenerateMockDataSet(fUseAverage); + } + + if (experiment->GetSignal() && experiment->GetBackground()) { + experiment->SetName(experiment->GetSignal()->GetName()); + fExperiments.push_back(experiment); + } + } + + if (fExposureTime > 0 && ToLower(fExposureStrategy) == "exp") { + ExtractExperimentParameterizationNodes(); + + Double_t sum = 0; + for (size_t n = 0; n < fExperiments.size(); n++) sum += TMath::Exp((double)n * fExposureFactor); + + Double_t A = fExposureTime * units("s") / sum; + for (size_t n = 0; n < fExperiments.size(); n++) { + fExperiments[n]->SetExposureInSeconds(A * TMath::Exp((double)n * fExposureFactor)); + fExperiments[n]->GenerateMockDataSet(fUseAverage); + } + } + + if (fExposureTime > 0 && ToLower(fExposureStrategy) == "power") { + ExtractExperimentParameterizationNodes(); + + Double_t sum = 0; + for (size_t n = 0; n < fExperiments.size(); n++) sum += TMath::Power((double)n, fExposureFactor); + + Double_t A = fExposureTime * units("s") / sum; + for (size_t n = 0; n < fExperiments.size(); n++) { + fExperiments[n]->SetExposureInSeconds(A * TMath::Power((double)n, fExposureFactor)); + fExperiments[n]->GenerateMockDataSet(fUseAverage); } + } - fExperiments.push_back(experiment); + if (fExposureTime > 0 && ToLower(fExposureStrategy) == "equal") { + ExtractExperimentParameterizationNodes(); + + for (size_t n = 0; n < fExperiments.size(); n++) { + fExperiments[n]->SetExposureInSeconds(fExposureTime * units("s") / fExperiments.size()); + fExperiments[n]->GenerateMockDataSet(fUseAverage); + } } } } +///////////////////////////////////////////// +/// \brief It scans all the experiment signals parametric nodes to build a complete list +/// of nodes used to build a complete sensitivity curve. Some experiments may be +/// sensitivy to a particular node, while others may be sensitivy to another. If more +/// than one experiment is sensitivy to a given node, the sensitivity will be combined +/// later on. +/// +void TRestExperimentList::ExtractExperimentParameterizationNodes() { + fParameterizationNodes.clear(); + + for (const auto& experiment : fExperiments) { + std::vector nodes = experiment->GetSignal()->GetParameterizationNodes(); + fParameterizationNodes.insert(fParameterizationNodes.end(), nodes.begin(), nodes.end()); + + std::sort(fParameterizationNodes.begin(), fParameterizationNodes.end()); + auto last = std::unique(fParameterizationNodes.begin(), fParameterizationNodes.end()); + fParameterizationNodes.erase(last, fParameterizationNodes.end()); + } +} + +void TRestExperimentList::PrintParameterizationNodes() { + std::cout << "Experiment list nodes: "; + for (const auto& node : fParameterizationNodes) std::cout << node << "\t"; + std::cout << std::endl; +} + TRestComponent* TRestExperimentList::GetComponent(std::string compName) { TRestComponent* component = nullptr; for (const auto& c : fComponentFiles) { @@ -244,5 +311,7 @@ void TRestExperimentList::PrintMetadata() { RESTMetadata << "Number of experiments loaded: " << fExperiments.size() << RESTendl; + if (fUseAverage) RESTMetadata << "Average MonteCarlo counts generation was enabled" << RESTendl; + RESTMetadata << "----" << RESTendl; } diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index 12fe749d6..b9dc66893 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -83,7 +83,8 @@ void TRestSensitivity::Initialize() { SetSectionName(this->ClassName()); } /// 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) { +Double_t TRestSensitivity::ApproachByFactor(Double_t node, Double_t g4, Double_t chi0, Double_t target, + Double_t factor) { if (factor == 1) { return 0; } @@ -92,7 +93,7 @@ Double_t TRestSensitivity::ApproachByFactor(Double_t g4, Double_t chi0, Double_t Double_t Chi2 = 0; do { Chi2 = 0; - for (const auto& exp : fExperiments) Chi2 += -2 * UnbinnedLogLikelihood(exp, g4); + for (const auto& exp : fExperiments) Chi2 += -2 * UnbinnedLogLikelihood(exp, node, g4); g4 = factor * g4; } while (Chi2 - chi0 < target); @@ -101,7 +102,7 @@ Double_t TRestSensitivity::ApproachByFactor(Double_t g4, Double_t chi0, Double_t /// Coarse movement to get to Chi2 below target (/2) do { Chi2 = 0; - for (const auto& exp : fExperiments) Chi2 += -2 * UnbinnedLogLikelihood(exp, g4); + for (const auto& exp : fExperiments) Chi2 += -2 * UnbinnedLogLikelihood(exp, node, g4); g4 = g4 / factor; } while (Chi2 - chi0 > target); @@ -109,21 +110,143 @@ Double_t TRestSensitivity::ApproachByFactor(Double_t g4, Double_t chi0, Double_t return g4 * factor; } +void TRestSensitivity::GenerateCurve() { + ExtractExperimentParameterizationNodes(); + + if (GetNumberOfCurves() > 0) + for (const auto& exp : fExperiments) { + exp->GenerateMockDataSet(); + } + + RESTInfo << "Generating sensitivity curve" << RESTendl; + std::vector curve; + for (const auto& node : fParameterizationNodes) { + RESTInfo << "Generating node : " << node << RESTendl; + curve.push_back(GetCoupling(node)); + } + fCurves.push_back(curve); + + RESTInfo << "Curve has been generated. You may use now TRestSensitivity::ExportCurve( fname.txt )." + << RESTendl; +} + +void TRestSensitivity::GenerateCurves(Int_t N) { + /* + std::cout << "TRestSensitivity::GenerateCurves." << std::endl; + std::cout << "There is a potential memory leak when generating several curves." << std::endl; + std::cout << "This code needs to be reviewed" << std::endl; + return; + */ + + for (int n = 0; n < N; n++) GenerateCurve(); +} + +std::vector TRestSensitivity::GetCurve(size_t n) { + if (n >= GetNumberOfCurves()) { + RESTWarning << "Requested curve number : " << n << " but only " << GetNumberOfCurves() << " generated" + << RESTendl; + return std::vector(); + } + return fCurves[n]; +} + +std::vector TRestSensitivity::GetAveragedCurve() { + if (GetNumberOfCurves() <= 0) return std::vector(); + + std::vector averagedCurve(fCurves[0].size(), 0.0); // Initialize with zeros + + for (const auto& row : fCurves) { + for (size_t i = 0; i < row.size(); ++i) { + averagedCurve[i] += row[i]; + } + } + + for (double& avg : averagedCurve) { + avg /= static_cast(fCurves.size()); + } + + return averagedCurve; +} + +void TRestSensitivity::ExportAveragedCurve(std::string fname, Double_t factor, Double_t power) { + std::vector curve = GetAveragedCurve(); + if (curve.empty()) std::cout << "Curve is empty" << std::endl; + if (curve.empty()) return; + + // Open a file for writing + std::ofstream outputFile(fname); + + // Check if the file is opened successfully + if (!outputFile) { + RESTError << "TRestSensitivity::ExportCurve. Error opening file for writing!" << RESTendl; + return; + } + + if (fParameterizationNodes.size() != curve.size()) { + RESTError << "TRestSensitivity::ExportCurve. Curve has not been properly generated" << RESTendl; + RESTError << "Parameterization nodes: " << fParameterizationNodes.size() << RESTendl; + RESTError << "Try invoking TRestSensitivity::GenerateCurve" << RESTendl; + return; + } + + int m = 0; + for (const auto& node : fParameterizationNodes) { + outputFile << node << " " << factor * TMath::Power(curve[m], power) << std::endl; + m++; + } + + outputFile.close(); + + RESTInfo << "TRestSensitivity::ExportCurve. File has been written successfully!" << RESTendl; +} + +void TRestSensitivity::ExportCurve(std::string fname, Double_t factor, Double_t power, int n) { + std::vector curve = GetCurve(n); + if (curve.empty()) return; + + // Open a file for writing + std::ofstream outputFile(fname); + + // Check if the file is opened successfully + if (!outputFile) { + RESTError << "TRestSensitivity::ExportCurve. Error opening file for writing!" << RESTendl; + return; + } + + if (fParameterizationNodes.size() != curve.size()) { + RESTError << "TRestSensitivity::ExportCurve. Curve has not been properly generated" << RESTendl; + RESTError << "Try invoking TRestSensitivity::GenerateCurve" << RESTendl; + return; + } + + int m = 0; + for (const auto& node : fParameterizationNodes) { + outputFile << node << " " << factor * TMath::Power(curve[m], power) << std::endl; + m++; + } + + outputFile.close(); + + RESTInfo << "TRestSensitivity::ExportCurve. File has been written successfully!" << RESTendl; +} + /////////////////////////////////////////////// /// \brief It will return the coupling value for which Chi=sigma /// -Double_t TRestSensitivity::GetCoupling(Double_t sigma, Double_t precision) { +Double_t TRestSensitivity::GetCoupling(Double_t node, Double_t sigma, Double_t precision) { Double_t Chi2_0 = 0; - for (const auto& exp : fExperiments) Chi2_0 += -2 * UnbinnedLogLikelihood(exp, 0); + for (const auto& exp : fExperiments) { + Chi2_0 += -2 * UnbinnedLogLikelihood(exp, node, 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); + g4 = ApproachByFactor(node, g4, Chi2_0, target, 2); + g4 = ApproachByFactor(node, g4, Chi2_0, target, 1.2); + g4 = ApproachByFactor(node, g4, Chi2_0, target, 1.02); + g4 = ApproachByFactor(node, g4, Chi2_0, target, 1.0002); return g4; } @@ -131,7 +254,8 @@ Double_t TRestSensitivity::GetCoupling(Double_t sigma, Double_t precision) { /////////////////////////////////////////////// /// \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 TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experiment, Double_t node, + Double_t g4) { Double_t lhood = 0; if (!experiment->IsDataReady()) { RESTError << "TRestSensitivity::UnbinnedLogLikelihood. Experiment " << experiment->GetName() @@ -139,10 +263,41 @@ Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experime return lhood; } + if (!experiment->GetSignal()->HasNodes()) { + RESTError << "Experiment signal : " << experiment->GetSignal()->GetName() << " has no nodes!" + << RESTendl; + return lhood; + } + + /// We check if the signal component is sensitive to that particular node + /// If not, this experiment will not contribute to that node + Int_t nd = experiment->GetSignal()->FindActiveNode(node); + if (nd >= 0) + experiment->GetSignal()->SetActiveNode(nd); + else { + RESTWarning << "Node : " << node << " not found in signal : " << experiment->GetSignal()->GetName() + << RESTendl; + return 0.0; + } + + /// We could check if background has also components, but for the moment we do not have a background + /// for each node, although it could be the case, if for example the background depends on the detector + /// conditions. For example, higher pressure inside the detector gains in signal sensitivity but + /// it will produce also higher background. + + if (experiment->GetBackground()->HasNodes()) { + RESTWarning + << "TRestSensitivity::UnbinnedLogLikelihood is not ready to have background parameter nodes!" + << RESTendl; + return 0.0; + } + Double_t signal = g4 * experiment->GetSignal()->GetTotalRate() * experiment->GetExposureInSeconds(); lhood = -signal; + if (experiment->GetExperimentalCounts() == 0) return lhood; + if (ROOT::IsImplicitMTEnabled()) ROOT::DisableImplicitMT(); std::vector> trackingData; @@ -169,12 +324,12 @@ Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experime /////////////////////////////////////////////// /// \brief /// -TH1D* TRestSensitivity::SignalStatisticalTest(Int_t N) { +TH1D* TRestSensitivity::SignalStatisticalTest(Double_t node, 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())); + Double_t coupling = TMath::Sqrt(TMath::Sqrt(GetCoupling(node))); couplings.push_back(coupling); } @@ -213,11 +368,254 @@ void TRestSensitivity::InitFromConfigFile() { Initialize(); } +///////////////////////////////////////////// +/// \brief This method is used to obtain the list of curves that satisfy that each value inside +/// the curve is placed at a specified level. E.g. if we provide a level 0.5, then the corresponding +/// curve will be constructed with the central value extracted at each parameter point. +// +/// We may then construct the profile of the sensitivity curves at 98%, 95% and 68% C.L. as follows: +/// +/// \code +/// TRestSensitivity::GetLevelCurves( {0.01, 0.025, 0.16, 0.84, 0.975, 0.99} ); +/// \endcode +/// +std::vector> TRestSensitivity::GetLevelCurves(const std::vector& levels) { + std::vector> curves(levels.size()); + + for (const auto& l : levels) { + if (l >= 1 || l <= 0) { + RESTError << "The level value should be between 0 and 1" << RESTendl; + return curves; + } + } + + std::vector intLevels; + for (const auto& l : levels) { + int val = (int)round(l * fCurves.size()); + if (val >= (int)fCurves.size()) val = fCurves.size() - 1; + if (val < 0) val = 0; + + intLevels.push_back(val); + } + + for (size_t m = 0; m < fCurves[0].size(); m++) { + std::vector v; + for (size_t n = 0; n < fCurves.size(); n++) v.push_back(fCurves[n][m]); + + std::sort(v.begin(), v.begin() + v.size()); + + for (size_t n = 0; n < intLevels.size(); n++) curves[n].push_back(v[intLevels[n]]); + } + + return curves; +} + +TCanvas* TRestSensitivity::DrawCurves() { + if (fCanvas != NULL) { + delete fCanvas; + fCanvas = NULL; + } + fCanvas = new TCanvas("canv", "This is the canvas title", 600, 450); + fCanvas->Draw(); + + TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97); + // pad1->Divide(2, 2); + pad1->Draw(); + + ////// Drawing reflectivity versus angle + // pad1->cd()->SetLogx(); + pad1->cd()->SetRightMargin(0.09); + pad1->cd()->SetLeftMargin(0.15); + pad1->cd()->SetBottomMargin(0.15); + + std::vector graphs; + + for (size_t n = 0; n < 20; n++) { + std::string grname = "gr" + IntegerToString(n); + TGraph* gr = new TGraph(); + gr->SetName(grname.c_str()); + for (size_t m = 0; m < this->GetCurve(n).size(); m++) + gr->SetPoint(gr->GetN(), fParameterizationNodes[m], + TMath::Sqrt(TMath::Sqrt(this->GetCurve(n)[m]))); + + gr->SetLineColorAlpha(kBlue + n, 0.3); + gr->SetLineWidth(1); + graphs.push_back(gr); + } + + TGraph* avGr = new TGraph(); + std::vector avCurve = GetAveragedCurve(); + for (size_t m = 0; m < avCurve.size(); m++) + avGr->SetPoint(avGr->GetN(), fParameterizationNodes[m], TMath::Sqrt(TMath::Sqrt(avCurve[m]))); + avGr->SetLineColor(kBlack); + avGr->SetLineWidth(2); + + graphs[0]->GetXaxis()->SetLimits(0, 0.25); + // graphs[0]->GetHistogram()->SetMaximum(1); + // graphs[0]->GetHistogram()->SetMinimum(lowRange); + + graphs[0]->GetXaxis()->SetTitle("mass [eV]"); + graphs[0]->GetXaxis()->SetTitleSize(0.05); + graphs[0]->GetXaxis()->SetLabelSize(0.05); + graphs[0]->GetYaxis()->SetTitle("g_{a#gamma} [10^{-10} GeV^{-1}]"); + graphs[0]->GetYaxis()->SetTitleOffset(1.5); + graphs[0]->GetYaxis()->SetTitleSize(0.05); + // graphs[0]->GetYaxis()->SetLabelSize(0.05); + // graphs[0]->GetYaxis()->SetLabelOffset(0.0); + // pad1->cd()->SetLogy(); + graphs[0]->Draw("AL"); + for (unsigned int n = 1; n < graphs.size(); n++) graphs[n]->Draw("L"); + avGr->Draw("L"); + + /* +Double_t lx1 = 0.6, ly1 = 0.75, lx2 = 0.9, ly2 = 0.95; +if (eLegendCoords.size() > 0) { + lx1 = eLegendCoords[0]; + ly1 = eLegendCoords[1]; + lx2 = eLegendCoords[2]; + ly2 = eLegendCoords[3]; +} +TLegend* legend = new TLegend(lx1, ly1, lx2, ly2); + +legend->SetTextSize(0.03); +legend->SetHeader("Energies", "C"); // option "C" allows to center the header +for (unsigned int n = 0; n < energies.size(); n++) { + std::string lname = "gr" + IntegerToString(n); + std::string ltitle = DoubleToString(energies[n]) + " keV"; + + legend->AddEntry(lname.c_str(), ltitle.c_str(), "l"); +} +legend->Draw(); + */ + + return fCanvas; +} + +TCanvas* TRestSensitivity::DrawLevelCurves() { + if (fCanvas != NULL) { + delete fCanvas; + fCanvas = NULL; + } + fCanvas = new TCanvas("canv", "This is the canvas title", 500, 400); + fCanvas->Draw(); + fCanvas->SetLeftMargin(0.15); + fCanvas->SetRightMargin(0.04); + fCanvas->SetLogx(); + + std::vector> levelCurves = GetLevelCurves({0.025, 0.16, 0.375, 0.625, 0.84, 0.975}); + + std::vector graphs; + for (size_t n = 0; n < levelCurves.size(); n++) { + std::string grname = "gr" + IntegerToString(n); + TGraph* gr = new TGraph(); + gr->SetName(grname.c_str()); + for (size_t m = 0; m < levelCurves[n].size(); m++) + gr->SetPoint(gr->GetN(), fParameterizationNodes[m], TMath::Sqrt(TMath::Sqrt(levelCurves[n][m]))); + + gr->SetLineColor(kBlack); + gr->SetLineWidth(1); + graphs.push_back(gr); + } + + TGraph* centralGr = new TGraph(); + std::vector centralCurve = GetLevelCurves({0.5})[0]; + for (size_t m = 0; m < centralCurve.size(); m++) + centralGr->SetPoint(centralGr->GetN(), fParameterizationNodes[m], + TMath::Sqrt(TMath::Sqrt(centralCurve[m]))); + centralGr->SetLineColor(kBlack); + centralGr->SetLineWidth(2); + centralGr->SetMarkerSize(0.1); + + graphs[0]->GetYaxis()->SetRangeUser(0, 0.5); + graphs[0]->GetXaxis()->SetRangeUser(0.001, 0.25); + graphs[0]->GetXaxis()->SetLimits(0.0001, 0.25); + graphs[0]->GetXaxis()->SetTitle("mass [eV]"); + graphs[0]->GetXaxis()->SetTitleSize(0.04); + graphs[0]->GetXaxis()->SetTitleOffset(1.15); + graphs[0]->GetXaxis()->SetLabelSize(0.04); + + // graphs[0]->GetYaxis()->SetLabelFont(43); + graphs[0]->GetYaxis()->SetTitle("g_{a#gamma} [10^{-10} GeV^{-1}]"); + graphs[0]->GetYaxis()->SetTitleOffset(1.5); + graphs[0]->GetYaxis()->SetTitleSize(0.04); + graphs[0]->GetYaxis()->SetLabelSize(0.04); + // graphs[0]->GetYaxis()->SetLabelOffset(0); + // graphs[0]->GetYaxis()->SetLabelFont(43); + graphs[0]->Draw("AL"); + + TGraph* randomGr = new TGraph(); + std::vector randomCurve = fCurves[13]; + for (size_t m = 0; m < randomCurve.size(); m++) + randomGr->SetPoint(randomGr->GetN(), fParameterizationNodes[m], + TMath::Sqrt(TMath::Sqrt(randomCurve[m]))); + randomGr->SetLineColor(kBlack); + randomGr->SetLineWidth(1); + randomGr->SetMarkerSize(0.3); + randomGr->SetMarkerStyle(4); + + std::vector shadeGraphs; + + int M = (int)levelCurves.size(); + for (int x = 0; x < M / 2; x++) { + TGraph* shade = new TGraph(); + int N = levelCurves[0].size(); + for (size_t m = 0; m < levelCurves[0].size(); m++) + shade->SetPoint(shade->GetN(), fParameterizationNodes[m], + TMath::Sqrt(TMath::Sqrt(levelCurves[x][m]))); + for (int m = N - 1; m >= 0; --m) + shade->SetPoint(shade->GetN(), fParameterizationNodes[m], + TMath::Sqrt(TMath::Sqrt(levelCurves[M - 1 - x][m]))); + shade->SetFillColorAlpha(kRed, 0.25); + shade->Draw("f"); + shadeGraphs.push_back(shade); + } + + for (unsigned int n = 1; n < graphs.size(); n++) graphs[n]->Draw("Lsame"); + randomGr->Draw("LPsame"); + // centralGr->Draw("Lsame"); + + return fCanvas; +} + +///////////////////////////////////////////// +/// \brief It scans all the experiment signals parametric nodes to build a complete list +/// of nodes used to build a complete sensitivity curve. Some experiments may be +/// sensitivy to a particular node, while others may be sensitivy to another. If more +/// than one experiment is sensitivy to a given node, the sensitivity will be combined +/// later on. +/// +void TRestSensitivity::ExtractExperimentParameterizationNodes(Bool_t rescan) { + if (fParameterizationNodes.empty() || rescan) { + fParameterizationNodes.clear(); + + for (const auto& experiment : fExperiments) { + std::vector nodes = experiment->GetSignal()->GetParameterizationNodes(); + fParameterizationNodes.insert(fParameterizationNodes.end(), nodes.begin(), nodes.end()); + + std::sort(fParameterizationNodes.begin(), fParameterizationNodes.end()); + auto last = std::unique(fParameterizationNodes.begin(), fParameterizationNodes.end()); + fParameterizationNodes.erase(last, fParameterizationNodes.end()); + } + } +} + +void TRestSensitivity::PrintParameterizationNodes() { + std::cout << "Curve sensitivity nodes: "; + for (const auto& node : fParameterizationNodes) std::cout << node << "\t"; + std::cout << std::endl; +} + ///////////////////////////////////////////// /// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux /// void TRestSensitivity::PrintMetadata() { TRestMetadata::PrintMetadata(); + RESTMetadata << " - Number of parameterization nodes : " << GetNumberOfNodes() << RESTendl; + RESTMetadata << " - Number of experiments loaded : " << GetNumberOfExperiments() << RESTendl; + RESTMetadata << " - Number of sensitivity curves generated : " << GetNumberOfCurves() << RESTendl; + RESTMetadata << " " << RESTendl; + RESTMetadata << " You may access experiment info using TRestSensitivity::GetExperiment(n)" << RESTendl; + RESTMetadata << "----" << RESTendl; } diff --git a/source/framework/tools/inc/TRestTools.h b/source/framework/tools/inc/TRestTools.h index 34e1e18ab..81069d187 100644 --- a/source/framework/tools/inc/TRestTools.h +++ b/source/framework/tools/inc/TRestTools.h @@ -127,7 +127,7 @@ class TRestTools { static std::string Execute(std::string cmd); - static std::string DownloadRemoteFile(const std::string& remoteFile); + static std::string DownloadRemoteFile(const std::string& remoteFile, bool pidPrefix = false); static int DownloadRemoteFile(std::string remoteFile, std::string localFile); static int UploadToServer(std::string localFile, std::string remoteFile, std::string methodUrl = ""); diff --git a/source/framework/tools/src/TRestTools.cxx b/source/framework/tools/src/TRestTools.cxx index 9502be51a..5c868261c 100644 --- a/source/framework/tools/src/TRestTools.cxx +++ b/source/framework/tools/src/TRestTools.cxx @@ -1121,7 +1121,7 @@ std::istream& TRestTools::GetLine(std::istream& is, std::string& t) { /// will be used, including scp, wget. Downloads to REST_USER_PATH + "/download/" + filename /// by default. /// -std::string TRestTools::DownloadRemoteFile(const string& url) { +std::string TRestTools::DownloadRemoteFile(const string& url, bool pidPrefix) { string pureName = TRestTools::GetPureFileName(url); if (pureName.empty()) { RESTWarning << "error! (TRestTools::DownloadRemoteFile): url is not a file!" << RESTendl; @@ -1133,7 +1133,8 @@ std::string TRestTools::DownloadRemoteFile(const string& url) { if (url.find("local:") == 0) { return Replace(url, "local:", ""); } else { - string fullpath = REST_USER_PATH + "/download/" + pureName; + string fullpath = + REST_USER_PATH + "/download/" + (pidPrefix ? "PID_" + ToString(getpid()) + "_" : "") + pureName; int out; int attempts = 10; do { @@ -1153,7 +1154,6 @@ std::string TRestTools::DownloadRemoteFile(const string& url) { return ""; } } - return ""; } /////////////////////////////////////////////// diff --git a/source/libraries/axion b/source/libraries/axion index 3800d90f4..3802d108b 160000 --- a/source/libraries/axion +++ b/source/libraries/axion @@ -1 +1 @@ -Subproject commit 3800d90f41622ec0c465f7b4d4c5f384e4d9433a +Subproject commit 3802d108bed118f4777ca142273aabaf368770c0