From 73637e0c5bea11bf4da15a07533c836d30eef27e Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Wed, 14 Feb 2024 10:20:45 +0100 Subject: [PATCH 01/22] TRestComponent. Adding methods to recover the parameterization nodes from the component --- .../sensitivity/inc/TRestComponent.h | 8 ++-- .../sensitivity/src/TRestComponent.cxx | 37 +++++++++++++------ 2 files changed, 31 insertions(+), 14 deletions(-) diff --git a/source/framework/sensitivity/inc/TRestComponent.h b/source/framework/sensitivity/inc/TRestComponent.h index 1dfd614dc..a4d8d213e 100644 --- a/source/framework/sensitivity/inc/TRestComponent.h +++ b/source/framework/sensitivity/inc/TRestComponent.h @@ -82,9 +82,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 +101,9 @@ class TRestComponent : public TRestMetadata { void Initialize() override; void RegenerateHistograms(UInt_t seed = 0); + /// It returns true if any nodes have been defined. + Bool_t HasNodes() { return !fParameterizationNodes.empty(); } + virtual void RegenerateActiveNodeDensity() {} std::string GetNature() const { return fNature; } @@ -113,6 +113,7 @@ class TRestComponent : public TRestMetadata { Int_t GetSamples() { return fSamples; } Int_t GetActiveNode() { return fActiveNode; } Double_t GetActiveNodeValue() { return fParameterizationNodes[fActiveNode]; } + std::vector GetParameterizationNodes() { return fParameterizationNodes; } std::vector GetVariables() const { return fVariables; } std::vector GetRanges() const { return fRanges; } @@ -127,6 +128,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; diff --git a/source/framework/sensitivity/src/TRestComponent.cxx b/source/framework/sensitivity/src/TRestComponent.cxx index 516bf8104..c435ac6a7 100644 --- a/source/framework/sensitivity/src/TRestComponent.cxx +++ b/source/framework/sensitivity/src/TRestComponent.cxx @@ -316,8 +316,7 @@ ROOT::RVecD TRestComponent::GetRandom() { if (!GetDensity()) { for (size_t n = 0; n < GetDimensions(); n++) result.push_back(0); RESTWarning << "TRestComponent::GetRandom. Component might not be initialized! Use " - "TRestComponent::Initialize()." - << RESTendl; + "TRestComponent::Initialize()." << RESTendl; return result; } @@ -347,11 +346,10 @@ ROOT::RDF::RNode TRestComponent::GetMonteCarloDataFrame(Int_t N) { /* Excluding Rndm from df */ std::vector columns = df.GetColumnNames(); std::vector excludeColumns = {"Rndm"}; - columns.erase(std::remove_if(columns.begin(), columns.end(), - [&excludeColumns](std::string elem) { - return std::find(excludeColumns.begin(), excludeColumns.end(), elem) != - excludeColumns.end(); - }), + columns.erase(std::remove_if(columns.begin(), columns.end(), [&excludeColumns](std::string elem) { + return std::find(excludeColumns.begin(), excludeColumns.end(), elem) != + excludeColumns.end(); + }), columns.end()); std::string user = getenv("USER"); @@ -377,15 +375,13 @@ TCanvas* TRestComponent::DrawComponent(std::vector drawVariables, TString drawOption) { if (drawVariables.size() > 2 || drawVariables.size() == 0) { RESTError << "TRestComponent::DrawComponent. The number of variables to be drawn must " - "be 1 or 2!" - << RESTendl; + "be 1 or 2!" << RESTendl; return fCanvas; } if (scanVariables.size() > 2 || scanVariables.size() == 0) { RESTError << "TRestComponent::DrawComponent. The number of variables to be scanned must " - "be 1 or 2!" - << RESTendl; + "be 1 or 2!" << RESTendl; return fCanvas; } @@ -627,6 +623,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. From df6104f45b3c70e0e7ce40334b84ad09525e0676 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Wed, 14 Feb 2024 10:29:26 +0100 Subject: [PATCH 02/22] TRestSensitivity::GenerateCurve method and other necessary updates added --- .../sensitivity/inc/TRestSensitivity.h | 21 ++- .../sensitivity/src/TRestSensitivity.cxx | 123 ++++++++++++++++-- 2 files changed, 127 insertions(+), 17 deletions(-) diff --git a/source/framework/sensitivity/inc/TRestSensitivity.h b/source/framework/sensitivity/inc/TRestSensitivity.h index 56c9d0bd6..7e0b51080 100644 --- a/source/framework/sensitivity/inc/TRestSensitivity.h +++ b/source/framework/sensitivity/inc/TRestSensitivity.h @@ -31,20 +31,33 @@ class TRestSensitivity : public TRestMetadata { /// A list of experimental conditions included to get a final sensitivity plot std::vector fExperiments; //< + /// The fusioned list of parameterization nodes found at each experiment signal + std::vector fParameterizationNodes; //< + + /// The calculated coupling for each parametric node + std::vector fCouplingForNode; //< + + /// It is used to generate a histogram with the signal distribution produced with different signal samples TH1D* fSignalTest = nullptr; protected: void InitFromConfigFile() override; - Double_t UnbinnedLogLikelihood(const TRestExperiment* experiment, Double_t g4 = 0); - Double_t ApproachByFactor(Double_t g4, Double_t chi0, Double_t target, Double_t factor); + 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(); + std::vector GetParameterizationNodes() { return fParameterizationNodes; } + void PrintParameterizationNodes(); + + Double_t GetCoupling(Double_t node, Double_t sigma = 2, Double_t precision = 0.01); + void GenerateCurve(); + void ExportCurve(std::string fname); - TH1D* SignalStatisticalTest(Int_t N); + TH1D* SignalStatisticalTest(Double_t node, Int_t N); std::vector GetExperiments() { return fExperiments; } TRestExperiment* GetExperiment(const size_t& n) { diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index 12fe749d6..5533bff2f 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -76,14 +76,19 @@ TRestSensitivity::TRestSensitivity(const char* cfgFileName, const std::string& n /// \brief It will initialize the data frame with the filelist and column names /// (or observables) that have been defined by the user. /// -void TRestSensitivity::Initialize() { SetSectionName(this->ClassName()); } +void TRestSensitivity::Initialize() { + SetSectionName(this->ClassName()); + + ExtractExperimentParameterizationNodes(); +} /////////////////////////////////////////////// /// \brief It will return a value of the coupling, g4, such that (chi-chi0) gets /// closer to the target value given by argument. The factor will be used to /// increase or decrease the coupling, and evaluate the likelihood. /// -Double_t TRestSensitivity::ApproachByFactor(Double_t g4, Double_t chi0, Double_t target, Double_t factor) { +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 +97,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 +106,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 +114,64 @@ Double_t TRestSensitivity::ApproachByFactor(Double_t g4, Double_t chi0, Double_t return g4 * factor; } +void TRestSensitivity::GenerateCurve() { + ExtractExperimentParameterizationNodes(); + + RESTInfo << "Generating sensitivity curve" << RESTendl; + fCouplingForNode.clear(); + for (const auto& node : fParameterizationNodes) { + RESTInfo << "Generating node : " << node << RESTendl; + fCouplingForNode.push_back(GetCoupling(node)); + } + + RESTInfo << "Curve has been generated. You may use now TRestSensitivity::ExportCurve( fname.txt )." + << RESTendl; +} + +void TRestSensitivity::ExportCurve(std::string fname) { + // 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() != fCouplingForNode.size()) { + RESTError << "TRestSensitivity::ExportCurve. Curve has not been generated" << RESTendl; + RESTError << "Try invoking TRestSensitivity::ExportCurve" << RESTendl; + return; + } + + int n = 0; + for (const auto& node : fParameterizationNodes) { + outputFile << node << " " << fCouplingForNode[n] << std::endl; + n++; + } + + 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 +179,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 +188,32 @@ Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experime 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 + 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 +240,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,6 +284,32 @@ void TRestSensitivity::InitFromConfigFile() { Initialize(); } +///////////////////////////////////////////// +/// \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() { + 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 /// From 24e5bf535c5b33b46c85cc8e778c5ca82d6a0022 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Wed, 14 Feb 2024 10:30:14 +0100 Subject: [PATCH 03/22] TRestExperiment. Now the experimental counts are registed as a metadata member --- .../sensitivity/inc/TRestExperiment.h | 6 +++++- .../sensitivity/src/TRestExperiment.cxx | 18 ++++++++++-------- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/source/framework/sensitivity/inc/TRestExperiment.h b/source/framework/sensitivity/inc/TRestExperiment.h index 7e25014f2..ee9efd923 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,9 @@ class TRestExperiment : public TRestMetadata { /// Only if it is true we will be able to calculate the LogLikelihood Bool_t fDataReady = false; //< + /// It keeps track on the number of counts inside the dataset + Int_t fExperimentalCounts = 0; //< + /// Internal process random generator TRandom3* fRandom = nullptr; //! @@ -64,6 +67,7 @@ class TRestExperiment : public TRestMetadata { public: void GenerateMockDataSet(); + Int_t GetExperimentalCounts() const { return fExperimentalCounts; } Bool_t IsMockData() const { return fMockData; } Bool_t IsDataReady() const { return fDataReady; } diff --git a/source/framework/sensitivity/src/TRestExperiment.cxx b/source/framework/sensitivity/src/TRestExperiment.cxx index 674c391c9..8cf526b00 100644 --- a/source/framework/sensitivity/src/TRestExperiment.cxx +++ b/source/framework/sensitivity/src/TRestExperiment.cxx @@ -102,12 +102,15 @@ void TRestExperiment::GenerateMockDataSet() { Double_t meanCounts = GetBackground()->GetTotalRate() * fExposureTime * units("s"); Int_t N = fRandom->Poisson(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,14 +121,14 @@ 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; if (!fSignal || !fBackground) { RESTWarning << "TRestExperiment::SetExperimentalDataSet. Signal and background components must " - "be available before atempt to load experimental data" - << RESTendl; + "be available before atempt to load experimental data" << RESTendl; fDataReady = false; return; } @@ -134,8 +137,7 @@ void TRestExperiment::SetExperimentalDataSet(const std::string& filename) { for (const auto& v : fSignal->GetVariables()) { if (std::find(columns.begin(), columns.end(), v) == columns.end()) { RESTError << "TRestExperiment::SetExperimentalDataSetFile a component variable was not found in " - "the dataset!" - << RESTendl; + "the dataset!" << RESTendl; fDataReady = false; return; } @@ -263,9 +265,9 @@ void TRestExperiment::PrintMetadata() { RESTMetadata << "Random seed : " << fSeed << RESTendl; RESTMetadata << " " << RESTendl; if (fExposureTime > 0) { - RESTMetadata << " - Exposure time : " << fExposureTime * units("s") << " seconds" << RESTendl; - RESTMetadata << " - Exposure time : " << fExposureTime * units("hr") << " hours" << RESTendl; - RESTMetadata << " - Exposure time : " << fExposureTime * units("day") << " days" << RESTendl; + RESTMetadata << " - Exposure time : " << fExposureTime* units("s") << " seconds" << RESTendl; + RESTMetadata << " - Exposure time : " << fExposureTime* units("hr") << " hours" << RESTendl; + RESTMetadata << " - Exposure time : " << fExposureTime* units("day") << " days" << RESTendl; } if (fSignal) RESTMetadata << " - Signal component : " << fSignal->GetName() << RESTendl; @@ -283,7 +285,7 @@ void TRestExperiment::PrintMetadata() { } } - RESTMetadata << " - Experimental counts : " << *fExperimentalData.GetDataFrame().Count() << RESTendl; + RESTMetadata << " - Experimental counts : " << fExperimentalCounts << RESTendl; RESTMetadata << "----" << RESTendl; } From 2832781fdd77a3debd835102fd77ee522c3f1f4a Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Wed, 14 Feb 2024 10:31:33 +0100 Subject: [PATCH 04/22] TRestExperimentList. Implementing KSVZ exposure strategy --- .../sensitivity/inc/TRestExperimentList.h | 10 +++ .../sensitivity/src/TRestExperimentList.cxx | 90 +++++++++++++++++-- 2 files changed, 92 insertions(+), 8 deletions(-) diff --git a/source/framework/sensitivity/inc/TRestExperimentList.h b/source/framework/sensitivity/inc/TRestExperimentList.h index 5d1737e40..2c99e57a4 100644 --- a/source/framework/sensitivity/inc/TRestExperimentList.h +++ b/source/framework/sensitivity/inc/TRestExperimentList.h @@ -50,9 +50,15 @@ 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; //< + /// 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 (unique/ksvz). + std::string fExposureStrategy = "unique"; + /// If not null this will be the common signal used in each experiment TRestComponent* fSignal = nullptr; //< @@ -71,6 +77,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()) diff --git a/source/framework/sensitivity/src/TRestExperimentList.cxx b/source/framework/sensitivity/src/TRestExperimentList.cxx index ab773d024..aa6f9ad2e 100644 --- a/source/framework/sensitivity/src/TRestExperimentList.cxx +++ b/source/framework/sensitivity/src/TRestExperimentList.cxx @@ -117,8 +117,7 @@ void TRestExperimentList::InitFromConfigFile() { if (nExpectedColumns == 0) { RESTError << "TRestExperimentList::InitFromConfigFile. At least one free parameter required! " - "(Exposure/Background/Signal)" - << RESTendl; + "(Exposure/Background/Signal)" << RESTendl; return; } @@ -135,8 +134,7 @@ void TRestExperimentList::InitFromConfigFile() { } RESTError << "TRestExperimentList::InitFromConfigFile. Number of expected columns does not match " - "the number of table columns" - << RESTendl; + "the number of table columns" << RESTendl; RESTError << "Number of table columns : " << nTableColumns << RESTendl; RESTError << "Number of expected columns : " << nExpectedColumns << RESTendl; RESTError << "Expected columns : " << expectedColumns << RESTendl; @@ -177,14 +175,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] @@ -212,11 +213,84 @@ void TRestExperimentList::InitFromConfigFile() { experiment->GenerateMockDataSet(); } - fExperiments.push_back(experiment); + if (experiment->GetSignal() && experiment->GetBackground()) { + experiment->SetName(experiment->GetSignal()->GetName()); + fExperiments.push_back(experiment); + } + } + + /// TODO. I am being lazy here. It would be more appropriate that + /// I create a TRestAxionExperimentList::TRestExperimentList + /// that implements this axion dedicated piece of code + if (fExposureTime > 0 && ToLower(fExposureStrategy) == "ksvz") { + ExtractExperimentParameterizationNodes(); + + Double_t Coefficient = 0; + for (const auto& node : fParameterizationNodes) { + Double_t expectedRate = 0; + for (const auto& experiment : fExperiments) { + Int_t nd = experiment->GetSignal()->FindActiveNode(node); + if (nd >= 0) { + experiment->GetSignal()->SetActiveNode(nd); + expectedRate += experiment->GetSignal()->GetTotalRate(); + } + } + + Coefficient += 1. / node / expectedRate; + } + + Double_t expectedRate = 0; + for (const auto& experiment : fExperiments) { + // We consider the contribution to each node for a given experiment + for (const auto& node : fParameterizationNodes) { + Int_t nd = experiment->GetSignal()->FindActiveNode(node); + if (nd >= 0) { + experiment->GetSignal()->SetActiveNode(nd); + expectedRate += node * experiment->GetSignal()->GetTotalRate(); + } + } + Double_t experimentTime = fExposureTime / Coefficient / expectedRate; + experiment->SetExposureInSeconds(experimentTime * units("s")); + } + + Double_t totalExp = 0; + for (const auto& experiment : fExperiments) totalExp += experiment->GetExposureInSeconds(); + + for (const auto& experiment : fExperiments) { + Double_t xp = experiment->GetExposureInSeconds(); + experiment->SetExposureInSeconds(xp * fExposureTime * units("s") / totalExp); + experiment->GenerateMockDataSet(); + } } } } +///////////////////////////////////////////// +/// \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) { From 54b941f5ce8d507ffc111eb8436102289da2fed8 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 14 Feb 2024 09:33:00 +0000 Subject: [PATCH 05/22] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../sensitivity/src/TRestComponent.cxx | 18 +++++++++++------- .../sensitivity/src/TRestExperiment.cxx | 12 +++++++----- .../sensitivity/src/TRestExperimentList.cxx | 6 ++++-- 3 files changed, 22 insertions(+), 14 deletions(-) diff --git a/source/framework/sensitivity/src/TRestComponent.cxx b/source/framework/sensitivity/src/TRestComponent.cxx index c435ac6a7..75906a898 100644 --- a/source/framework/sensitivity/src/TRestComponent.cxx +++ b/source/framework/sensitivity/src/TRestComponent.cxx @@ -316,7 +316,8 @@ ROOT::RVecD TRestComponent::GetRandom() { if (!GetDensity()) { for (size_t n = 0; n < GetDimensions(); n++) result.push_back(0); RESTWarning << "TRestComponent::GetRandom. Component might not be initialized! Use " - "TRestComponent::Initialize()." << RESTendl; + "TRestComponent::Initialize()." + << RESTendl; return result; } @@ -346,10 +347,11 @@ ROOT::RDF::RNode TRestComponent::GetMonteCarloDataFrame(Int_t N) { /* Excluding Rndm from df */ std::vector columns = df.GetColumnNames(); std::vector excludeColumns = {"Rndm"}; - columns.erase(std::remove_if(columns.begin(), columns.end(), [&excludeColumns](std::string elem) { - return std::find(excludeColumns.begin(), excludeColumns.end(), elem) != - excludeColumns.end(); - }), + columns.erase(std::remove_if(columns.begin(), columns.end(), + [&excludeColumns](std::string elem) { + return std::find(excludeColumns.begin(), excludeColumns.end(), elem) != + excludeColumns.end(); + }), columns.end()); std::string user = getenv("USER"); @@ -375,13 +377,15 @@ TCanvas* TRestComponent::DrawComponent(std::vector drawVariables, TString drawOption) { if (drawVariables.size() > 2 || drawVariables.size() == 0) { RESTError << "TRestComponent::DrawComponent. The number of variables to be drawn must " - "be 1 or 2!" << RESTendl; + "be 1 or 2!" + << RESTendl; return fCanvas; } if (scanVariables.size() > 2 || scanVariables.size() == 0) { RESTError << "TRestComponent::DrawComponent. The number of variables to be scanned must " - "be 1 or 2!" << RESTendl; + "be 1 or 2!" + << RESTendl; return fCanvas; } diff --git a/source/framework/sensitivity/src/TRestExperiment.cxx b/source/framework/sensitivity/src/TRestExperiment.cxx index 8cf526b00..65e331608 100644 --- a/source/framework/sensitivity/src/TRestExperiment.cxx +++ b/source/framework/sensitivity/src/TRestExperiment.cxx @@ -128,7 +128,8 @@ void TRestExperiment::SetExperimentalDataSet(const std::string& filename) { if (!fSignal || !fBackground) { RESTWarning << "TRestExperiment::SetExperimentalDataSet. Signal and background components must " - "be available before atempt to load experimental data" << RESTendl; + "be available before atempt to load experimental data" + << RESTendl; fDataReady = false; return; } @@ -137,7 +138,8 @@ void TRestExperiment::SetExperimentalDataSet(const std::string& filename) { for (const auto& v : fSignal->GetVariables()) { if (std::find(columns.begin(), columns.end(), v) == columns.end()) { RESTError << "TRestExperiment::SetExperimentalDataSetFile a component variable was not found in " - "the dataset!" << RESTendl; + "the dataset!" + << RESTendl; fDataReady = false; return; } @@ -265,9 +267,9 @@ void TRestExperiment::PrintMetadata() { RESTMetadata << "Random seed : " << fSeed << RESTendl; RESTMetadata << " " << RESTendl; if (fExposureTime > 0) { - RESTMetadata << " - Exposure time : " << fExposureTime* units("s") << " seconds" << RESTendl; - RESTMetadata << " - Exposure time : " << fExposureTime* units("hr") << " hours" << RESTendl; - RESTMetadata << " - Exposure time : " << fExposureTime* units("day") << " days" << RESTendl; + RESTMetadata << " - Exposure time : " << fExposureTime * units("s") << " seconds" << RESTendl; + RESTMetadata << " - Exposure time : " << fExposureTime * units("hr") << " hours" << RESTendl; + RESTMetadata << " - Exposure time : " << fExposureTime * units("day") << " days" << RESTendl; } if (fSignal) RESTMetadata << " - Signal component : " << fSignal->GetName() << RESTendl; diff --git a/source/framework/sensitivity/src/TRestExperimentList.cxx b/source/framework/sensitivity/src/TRestExperimentList.cxx index aa6f9ad2e..18050b92a 100644 --- a/source/framework/sensitivity/src/TRestExperimentList.cxx +++ b/source/framework/sensitivity/src/TRestExperimentList.cxx @@ -117,7 +117,8 @@ void TRestExperimentList::InitFromConfigFile() { if (nExpectedColumns == 0) { RESTError << "TRestExperimentList::InitFromConfigFile. At least one free parameter required! " - "(Exposure/Background/Signal)" << RESTendl; + "(Exposure/Background/Signal)" + << RESTendl; return; } @@ -134,7 +135,8 @@ void TRestExperimentList::InitFromConfigFile() { } RESTError << "TRestExperimentList::InitFromConfigFile. Number of expected columns does not match " - "the number of table columns" << RESTendl; + "the number of table columns" + << RESTendl; RESTError << "Number of table columns : " << nTableColumns << RESTendl; RESTError << "Number of expected columns : " << nExpectedColumns << RESTendl; RESTError << "Expected columns : " << expectedColumns << RESTendl; From 1a3e2a39830cf4e3b650258727db2a4250729e5d Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 20 Feb 2024 20:29:42 +0100 Subject: [PATCH 06/22] TRestDataSet. Adding cuts documentation --- source/framework/core/src/TRestDataSet.cxx | 46 ++++++++++++++-------- 1 file changed, 29 insertions(+), 17 deletions(-) diff --git a/source/framework/core/src/TRestDataSet.cxx b/source/framework/core/src/TRestDataSet.cxx index 56011722a..3a35ca4ee 100644 --- a/source/framework/core/src/TRestDataSet.cxx +++ b/source/framework/core/src/TRestDataSet.cxx @@ -241,6 +241,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 @@ -353,7 +368,7 @@ void TRestDataSet::GenerateDataSet() { fDataSet = MakeCut(fCut); // Adding new user columns added to the dataset - for (const auto& [cName, cExpression] : fColumnNameExpressions) { + for (const auto & [ cName, cExpression ] : fColumnNameExpressions) { RESTInfo << "Adding column to dataset: " << cName << RESTendl; finalList.emplace_back(cName); fDataSet = DefineColumn(cName, cExpression); @@ -384,8 +399,7 @@ std::vector TRestDataSet::FileSelection() { if (!time_stamp_end || !time_stamp_start) { RESTError << "TRestDataSet::FileSelect. Start or end dates not properly formed. Please, check " - "REST_StringHelper::StringToTimeStamp documentation for valid formats" - << RESTendl; + "REST_StringHelper::StringToTimeStamp documentation for valid formats" << RESTendl; return fFileSelection; } @@ -438,7 +452,7 @@ std::vector TRestDataSet::FileSelection() { if (!accept) continue; Double_t acc = 0; - for (auto& [name, properties] : fQuantity) { + for (auto & [ name, properties ] : fQuantity) { std::string value = run.ReplaceMetadataMembers(properties.metadata); const Double_t val = REST_StringHelper::StringToDouble(value); @@ -483,7 +497,7 @@ std::vector TRestDataSet::FileSelection() { } /////////////////////////////////////////////// -/// \brief This function apply a TRestCut to the dataframe +/// \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); @@ -495,7 +509,7 @@ ROOT::RDF::RNode TRestDataSet::MakeCut(const TRestCut* cut) { auto paramCut = cut->GetParamCut(); auto obsList = df.GetColumnNames(); - for (const auto& [param, condition] : paramCut) { + for (const auto & [ param, condition ] : paramCut) { if (std::find(obsList.begin(), obsList.end(), param) != obsList.end()) { std::string pCut = param + condition; RESTDebug << "Applying cut " << pCut << RESTendl; @@ -542,7 +556,7 @@ ROOT::RDF::RNode TRestDataSet::DefineColumn(const std::string& columnName, const auto df = fDataSet; std::string evalFormula = formula; - for (auto const& [name, properties] : fQuantity) + for (auto const & [ name, properties ] : fQuantity) evalFormula = REST_StringHelper::Replace(evalFormula, name, properties.value); df = df.Define(columnName, evalFormula); @@ -608,7 +622,7 @@ void TRestDataSet::PrintMetadata() { RESTMetadata << " Relevant quantities: " << RESTendl; RESTMetadata << " -------------------- " << RESTendl; - for (auto const& [name, properties] : fQuantity) { + for (auto const & [ name, properties ] : fQuantity) { RESTMetadata << " - Name : " << name << ". Value : " << properties.value << ". Strategy: " << properties.strategy << RESTendl; RESTMetadata << " - Metadata: " << properties.metadata << RESTendl; @@ -620,7 +634,7 @@ void TRestDataSet::PrintMetadata() { if (!fColumnNameExpressions.empty()) { RESTMetadata << " New columns added to generated dataframe: " << RESTendl; RESTMetadata << " ---------------------------------------- " << RESTendl; - for (const auto& [cName, cExpression] : fColumnNameExpressions) { + for (const auto & [ cName, cExpression ] : fColumnNameExpressions) { RESTMetadata << " - Name : " << cName << RESTendl; RESTMetadata << " - Expression: " << cExpression << RESTendl; RESTMetadata << " " << RESTendl; @@ -786,11 +800,10 @@ void TRestDataSet::Export(const std::string& filename, std::vector std::vector columns = fDataSet.GetColumnNames(); if (!excludeColumns.empty()) { - columns.erase(std::remove_if(columns.begin(), columns.end(), - [&excludeColumns](std::string elem) { - return std::find(excludeColumns.begin(), excludeColumns.end(), - elem) != excludeColumns.end(); - }), + columns.erase(std::remove_if(columns.begin(), columns.end(), [&excludeColumns](std::string elem) { + return std::find(excludeColumns.begin(), excludeColumns.end(), elem) != + excludeColumns.end(); + }), columns.end()); RESTInfo << "Re-Generating snapshot." << RESTendl; @@ -825,8 +838,7 @@ void TRestDataSet::Export(const std::string& filename, std::vector if (type != "Double_t" && type != "Int_t") { RESTError << "Branch name : " << bName << " is type : " << type << RESTendl; RESTError << "Only Int_t and Double_t types are allowed for " - "exporting to ASCII table" - << RESTendl; + "exporting to ASCII table" << RESTendl; RESTError << "File will not be generated" << RESTendl; return; } @@ -861,7 +873,7 @@ void TRestDataSet::Export(const std::string& filename, std::vector } fprintf(f, "###\n"); fprintf(f, "### Relevant quantities: \n"); - for (auto& [name, properties] : fQuantity) { + for (auto & [ name, properties ] : fQuantity) { fprintf(f, "### - %s : %s - %s\n", name.c_str(), properties.value.c_str(), properties.description.c_str()); } From 26047138a2b55eff3cbd88499f6ad69f6158cfda Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 20 Feb 2024 19:31:23 +0000 Subject: [PATCH 07/22] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- source/framework/core/src/TRestDataSet.cxx | 29 ++++++++++++---------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/source/framework/core/src/TRestDataSet.cxx b/source/framework/core/src/TRestDataSet.cxx index 3a35ca4ee..a67092a52 100644 --- a/source/framework/core/src/TRestDataSet.cxx +++ b/source/framework/core/src/TRestDataSet.cxx @@ -368,7 +368,7 @@ void TRestDataSet::GenerateDataSet() { fDataSet = MakeCut(fCut); // Adding new user columns added to the dataset - for (const auto & [ cName, cExpression ] : fColumnNameExpressions) { + for (const auto& [cName, cExpression] : fColumnNameExpressions) { RESTInfo << "Adding column to dataset: " << cName << RESTendl; finalList.emplace_back(cName); fDataSet = DefineColumn(cName, cExpression); @@ -399,7 +399,8 @@ std::vector TRestDataSet::FileSelection() { if (!time_stamp_end || !time_stamp_start) { RESTError << "TRestDataSet::FileSelect. Start or end dates not properly formed. Please, check " - "REST_StringHelper::StringToTimeStamp documentation for valid formats" << RESTendl; + "REST_StringHelper::StringToTimeStamp documentation for valid formats" + << RESTendl; return fFileSelection; } @@ -452,7 +453,7 @@ std::vector TRestDataSet::FileSelection() { if (!accept) continue; Double_t acc = 0; - for (auto & [ name, properties ] : fQuantity) { + for (auto& [name, properties] : fQuantity) { std::string value = run.ReplaceMetadataMembers(properties.metadata); const Double_t val = REST_StringHelper::StringToDouble(value); @@ -509,7 +510,7 @@ ROOT::RDF::RNode TRestDataSet::MakeCut(const TRestCut* cut) { auto paramCut = cut->GetParamCut(); auto obsList = df.GetColumnNames(); - for (const auto & [ param, condition ] : paramCut) { + for (const auto& [param, condition] : paramCut) { if (std::find(obsList.begin(), obsList.end(), param) != obsList.end()) { std::string pCut = param + condition; RESTDebug << "Applying cut " << pCut << RESTendl; @@ -556,7 +557,7 @@ ROOT::RDF::RNode TRestDataSet::DefineColumn(const std::string& columnName, const auto df = fDataSet; std::string evalFormula = formula; - for (auto const & [ name, properties ] : fQuantity) + for (auto const& [name, properties] : fQuantity) evalFormula = REST_StringHelper::Replace(evalFormula, name, properties.value); df = df.Define(columnName, evalFormula); @@ -622,7 +623,7 @@ void TRestDataSet::PrintMetadata() { RESTMetadata << " Relevant quantities: " << RESTendl; RESTMetadata << " -------------------- " << RESTendl; - for (auto const & [ name, properties ] : fQuantity) { + for (auto const& [name, properties] : fQuantity) { RESTMetadata << " - Name : " << name << ". Value : " << properties.value << ". Strategy: " << properties.strategy << RESTendl; RESTMetadata << " - Metadata: " << properties.metadata << RESTendl; @@ -634,7 +635,7 @@ void TRestDataSet::PrintMetadata() { if (!fColumnNameExpressions.empty()) { RESTMetadata << " New columns added to generated dataframe: " << RESTendl; RESTMetadata << " ---------------------------------------- " << RESTendl; - for (const auto & [ cName, cExpression ] : fColumnNameExpressions) { + for (const auto& [cName, cExpression] : fColumnNameExpressions) { RESTMetadata << " - Name : " << cName << RESTendl; RESTMetadata << " - Expression: " << cExpression << RESTendl; RESTMetadata << " " << RESTendl; @@ -800,10 +801,11 @@ void TRestDataSet::Export(const std::string& filename, std::vector std::vector columns = fDataSet.GetColumnNames(); if (!excludeColumns.empty()) { - columns.erase(std::remove_if(columns.begin(), columns.end(), [&excludeColumns](std::string elem) { - return std::find(excludeColumns.begin(), excludeColumns.end(), elem) != - excludeColumns.end(); - }), + columns.erase(std::remove_if(columns.begin(), columns.end(), + [&excludeColumns](std::string elem) { + return std::find(excludeColumns.begin(), excludeColumns.end(), + elem) != excludeColumns.end(); + }), columns.end()); RESTInfo << "Re-Generating snapshot." << RESTendl; @@ -838,7 +840,8 @@ void TRestDataSet::Export(const std::string& filename, std::vector if (type != "Double_t" && type != "Int_t") { RESTError << "Branch name : " << bName << " is type : " << type << RESTendl; RESTError << "Only Int_t and Double_t types are allowed for " - "exporting to ASCII table" << RESTendl; + "exporting to ASCII table" + << RESTendl; RESTError << "File will not be generated" << RESTendl; return; } @@ -873,7 +876,7 @@ void TRestDataSet::Export(const std::string& filename, std::vector } fprintf(f, "###\n"); fprintf(f, "### Relevant quantities: \n"); - for (auto & [ name, properties ] : fQuantity) { + for (auto& [name, properties] : fQuantity) { fprintf(f, "### - %s : %s - %s\n", name.c_str(), properties.value.c_str(), properties.description.c_str()); } From bd94cf7e58a13d466e98d531ff0fc057032e4ca4 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 20 Feb 2024 20:37:25 +0100 Subject: [PATCH 08/22] TRestDataSet. Adding documentation for opening dataset --- source/framework/core/src/TRestDataSet.cxx | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/source/framework/core/src/TRestDataSet.cxx b/source/framework/core/src/TRestDataSet.cxx index 3a35ca4ee..fc6961cff 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 From e9ee1da81051116828226d2eac59b90967c41926 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Wed, 21 Feb 2024 12:08:57 +0100 Subject: [PATCH 09/22] TRestSensitivity. Adding capability to generate and store several experimental curves. Method to obtain an average added too --- .../sensitivity/inc/TRestSensitivity.h | 14 ++- .../sensitivity/src/TRestSensitivity.cxx | 117 ++++++++++++++---- 2 files changed, 106 insertions(+), 25 deletions(-) diff --git a/source/framework/sensitivity/inc/TRestSensitivity.h b/source/framework/sensitivity/inc/TRestSensitivity.h index 7e0b51080..201ea41d8 100644 --- a/source/framework/sensitivity/inc/TRestSensitivity.h +++ b/source/framework/sensitivity/inc/TRestSensitivity.h @@ -35,7 +35,7 @@ class TRestSensitivity : public TRestMetadata { std::vector fParameterizationNodes; //< /// The calculated coupling for each parametric node - std::vector fCouplingForNode; //< + std::vector> fCouplingsForNode; //< /// It is used to generate a histogram with the signal distribution produced with different signal samples TH1D* fSignalTest = nullptr; @@ -49,13 +49,21 @@ class TRestSensitivity : public TRestMetadata { public: void Initialize() override; - void ExtractExperimentParameterizationNodes(); + 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 GenerateCurve(); - void ExportCurve(std::string fname); + void GenerateCurves(Int_t N); + + std::vector GetAveragedCurve(); + + void ExportCurve(std::string fname, int n); + void ExportAveragedCurve(std::string fname); + + size_t GetNumberOfCurves() { return fCouplingsForNode.size(); } + std::vector GetSensitivityCurve(size_t n = 0); TH1D* SignalStatisticalTest(Double_t node, Int_t N); diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index 5533bff2f..3d652b7fc 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -76,11 +76,7 @@ TRestSensitivity::TRestSensitivity(const char* cfgFileName, const std::string& n /// \brief It will initialize the data frame with the filelist and column names /// (or observables) that have been defined by the user. /// -void TRestSensitivity::Initialize() { - SetSectionName(this->ClassName()); - - ExtractExperimentParameterizationNodes(); -} +void TRestSensitivity::Initialize() { SetSectionName(this->ClassName()); } /////////////////////////////////////////////// /// \brief It will return a value of the coupling, g4, such that (chi-chi0) gets @@ -117,18 +113,92 @@ Double_t TRestSensitivity::ApproachByFactor(Double_t node, Double_t g4, Double_t void TRestSensitivity::GenerateCurve() { ExtractExperimentParameterizationNodes(); + if (GetNumberOfCurves() > 0) + for (const auto& exp : fExperiments) { + exp->GenerateMockDataSet(); + } + RESTInfo << "Generating sensitivity curve" << RESTendl; - fCouplingForNode.clear(); + std::vector curve; for (const auto& node : fParameterizationNodes) { RESTInfo << "Generating node : " << node << RESTendl; - fCouplingForNode.push_back(GetCoupling(node)); + curve.push_back(GetCoupling(node)); } + fCouplingsForNode.push_back(curve); RESTInfo << "Curve has been generated. You may use now TRestSensitivity::ExportCurve( fname.txt )." << RESTendl; } -void TRestSensitivity::ExportCurve(std::string fname) { +void TRestSensitivity::GenerateCurves(Int_t N) { + for (int n = 0; n < N; n++) GenerateCurve(); +} + +std::vector TRestSensitivity::GetSensitivityCurve(size_t n) { + if (n >= GetNumberOfCurves()) { + RESTWarning << "Requested curve number : " << n << " but only " << GetNumberOfCurves() << " generated" + << RESTendl; + return std::vector(); + } + return fCouplingsForNode[n]; +} + +std::vector TRestSensitivity::GetAveragedCurve() { + if (GetNumberOfCurves() <= 0) return std::vector(); + + std::cout << "Points : " << fCouplingsForNode[0].size() << std::endl; + std::vector averagedCurve(fCouplingsForNode[0].size(), 0.0); // Initialize with zeros + + for (const auto& row : fCouplingsForNode) { + for (size_t i = 0; i < row.size(); ++i) { + averagedCurve[i] += row[i]; + } + } + + for (double& avg : averagedCurve) { + avg /= static_cast(fCouplingsForNode.size()); + } + + return averagedCurve; +} + +void TRestSensitivity::ExportAveragedCurve(std::string fname) { + 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 << " " << curve[m] << std::endl; + m++; + } + + outputFile.close(); + + RESTInfo << "TRestSensitivity::ExportCurve. File has been written successfully!" << RESTendl; +} + +void TRestSensitivity::ExportCurve(std::string fname, int n = 0) { + + std::vector curve = GetSensitivityCurve(n); + if (curve.empty()) return; + // Open a file for writing std::ofstream outputFile(fname); @@ -138,16 +208,16 @@ void TRestSensitivity::ExportCurve(std::string fname) { return; } - if (fParameterizationNodes.size() != fCouplingForNode.size()) { - RESTError << "TRestSensitivity::ExportCurve. Curve has not been generated" << RESTendl; - RESTError << "Try invoking TRestSensitivity::ExportCurve" << RESTendl; + if (fParameterizationNodes.size() != curve.size()) { + RESTError << "TRestSensitivity::ExportCurve. Curve has not been properly generated" << RESTendl; + RESTError << "Try invoking TRestSensitivity::GenerateCurve" << RESTendl; return; } - int n = 0; + int m = 0; for (const auto& node : fParameterizationNodes) { - outputFile << node << " " << fCouplingForNode[n] << std::endl; - n++; + outputFile << node << " " << curve[m] << std::endl; + m++; } outputFile.close(); @@ -291,16 +361,19 @@ void TRestSensitivity::InitFromConfigFile() { /// than one experiment is sensitivy to a given node, the sensitivity will be combined /// later on. /// -void TRestSensitivity::ExtractExperimentParameterizationNodes() { - fParameterizationNodes.clear(); +void TRestSensitivity::ExtractExperimentParameterizationNodes(Bool_t rescan) { - for (const auto& experiment : fExperiments) { - std::vector nodes = experiment->GetSignal()->GetParameterizationNodes(); - fParameterizationNodes.insert(fParameterizationNodes.end(), nodes.begin(), nodes.end()); + if (fParameterizationNodes.empty() || rescan) { + fParameterizationNodes.clear(); - std::sort(fParameterizationNodes.begin(), fParameterizationNodes.end()); - auto last = std::unique(fParameterizationNodes.begin(), fParameterizationNodes.end()); - fParameterizationNodes.erase(last, fParameterizationNodes.end()); + 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()); + } } } From da5b3e4c2d8d3385b99f2bfedd71a07b3af72c90 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 21 Feb 2024 11:09:54 +0000 Subject: [PATCH 10/22] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- source/framework/sensitivity/src/TRestSensitivity.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index 3d652b7fc..64e2c9d18 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -195,7 +195,6 @@ void TRestSensitivity::ExportAveragedCurve(std::string fname) { } void TRestSensitivity::ExportCurve(std::string fname, int n = 0) { - std::vector curve = GetSensitivityCurve(n); if (curve.empty()) return; @@ -362,7 +361,6 @@ void TRestSensitivity::InitFromConfigFile() { /// later on. /// void TRestSensitivity::ExtractExperimentParameterizationNodes(Bool_t rescan) { - if (fParameterizationNodes.empty() || rescan) { fParameterizationNodes.clear(); From a16b7b1c3307e3c67b5c0413e4cf56c62552c06a Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Wed, 21 Feb 2024 18:43:15 +0100 Subject: [PATCH 11/22] TRestSensitivity::PrintMetadata implemented --- .../sensitivity/inc/TRestSensitivity.h | 14 ++++++++++---- .../sensitivity/src/TRestSensitivity.cxx | 17 +++++++++++------ 2 files changed, 21 insertions(+), 10 deletions(-) diff --git a/source/framework/sensitivity/inc/TRestSensitivity.h b/source/framework/sensitivity/inc/TRestSensitivity.h index 201ea41d8..000d88388 100644 --- a/source/framework/sensitivity/inc/TRestSensitivity.h +++ b/source/framework/sensitivity/inc/TRestSensitivity.h @@ -34,8 +34,11 @@ class TRestSensitivity : public TRestMetadata { /// The fusioned list of parameterization nodes found at each experiment signal std::vector fParameterizationNodes; //< - /// The calculated coupling for each parametric node - std::vector> fCouplingsForNode; //< + /// 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; @@ -62,11 +65,12 @@ class TRestSensitivity : public TRestMetadata { void ExportCurve(std::string fname, int n); void ExportAveragedCurve(std::string fname); - size_t GetNumberOfCurves() { return fCouplingsForNode.size(); } std::vector GetSensitivityCurve(size_t n = 0); TH1D* SignalStatisticalTest(Double_t node, Int_t N); + void Freeze() { fFrozen = true; } + std::vector GetExperiments() { return fExperiments; } TRestExperiment* GetExperiment(const size_t& n) { if (n >= GetNumberOfExperiments()) @@ -76,6 +80,8 @@ 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; @@ -83,6 +89,6 @@ class TRestSensitivity : public TRestMetadata { TRestSensitivity(); ~TRestSensitivity(); - ClassDefOverride(TRestSensitivity, 1); + ClassDefOverride(TRestSensitivity, 2); }; #endif diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index 3d652b7fc..688f59dab 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -124,7 +124,7 @@ void TRestSensitivity::GenerateCurve() { RESTInfo << "Generating node : " << node << RESTendl; curve.push_back(GetCoupling(node)); } - fCouplingsForNode.push_back(curve); + fCurves.push_back(curve); RESTInfo << "Curve has been generated. You may use now TRestSensitivity::ExportCurve( fname.txt )." << RESTendl; @@ -140,23 +140,22 @@ std::vector TRestSensitivity::GetSensitivityCurve(size_t n) { << RESTendl; return std::vector(); } - return fCouplingsForNode[n]; + return fCurves[n]; } std::vector TRestSensitivity::GetAveragedCurve() { if (GetNumberOfCurves() <= 0) return std::vector(); - std::cout << "Points : " << fCouplingsForNode[0].size() << std::endl; - std::vector averagedCurve(fCouplingsForNode[0].size(), 0.0); // Initialize with zeros + std::vector averagedCurve(fCurves[0].size(), 0.0); // Initialize with zeros - for (const auto& row : fCouplingsForNode) { + 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(fCouplingsForNode.size()); + avg /= static_cast(fCurves.size()); } return averagedCurve; @@ -389,5 +388,11 @@ void TRestSensitivity::PrintParameterizationNodes() { 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; } From d6d24752b0dc0b35e3d849d60d6db7a51b388602 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Thu, 22 Feb 2024 20:45:26 +0100 Subject: [PATCH 12/22] REST_OpenInputFile.C Fixing typo in RDataFrame helper --- macros/REST_OpenInputFile.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/macros/REST_OpenInputFile.C b/macros/REST_OpenInputFile.C index 7812d8f5f..9037d6402 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); From f40b7ded1be717eb4e208b92b60af5b3ca3256f3 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Mon, 11 Mar 2024 09:15:51 +0100 Subject: [PATCH 13/22] REST_OpenInputFile.C fixed a typo --- macros/REST_OpenInputFile.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/macros/REST_OpenInputFile.C b/macros/REST_OpenInputFile.C index 9037d6402..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); From 672446e5770bf291a50c9336ff0da42a0a5fade1 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Mon, 18 Mar 2024 14:30:45 +0100 Subject: [PATCH 14/22] TRestSenstivity. Adding DrawCurves --- .../sensitivity/inc/TRestSensitivity.h | 11 ++- .../sensitivity/src/TRestSensitivity.cxx | 85 ++++++++++++++++++- 2 files changed, 91 insertions(+), 5 deletions(-) diff --git a/source/framework/sensitivity/inc/TRestSensitivity.h b/source/framework/sensitivity/inc/TRestSensitivity.h index 000d88388..0c79c404e 100644 --- a/source/framework/sensitivity/inc/TRestSensitivity.h +++ b/source/framework/sensitivity/inc/TRestSensitivity.h @@ -29,7 +29,7 @@ class TRestSensitivity : public TRestMetadata { private: /// A list of experimental conditions included to get a final sensitivity plot - std::vector fExperiments; //< + std::vector fExperiments; //! /// The fusioned list of parameterization nodes found at each experiment signal std::vector fParameterizationNodes; //< @@ -43,6 +43,9 @@ class TRestSensitivity : public TRestMetadata { /// 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; @@ -57,16 +60,16 @@ class TRestSensitivity : public TRestMetadata { 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 GenerateCurve(); void GenerateCurves(Int_t N); + std::vector GetCurve(size_t n = 0); std::vector GetAveragedCurve(); void ExportCurve(std::string fname, int n); void ExportAveragedCurve(std::string fname); - std::vector GetSensitivityCurve(size_t n = 0); - TH1D* SignalStatisticalTest(Double_t node, Int_t N); void Freeze() { fFrozen = true; } @@ -85,6 +88,8 @@ class TRestSensitivity : public TRestMetadata { void PrintMetadata() override; + TCanvas* DrawCurves(); + TRestSensitivity(const char* cfgFileName, const std::string& name = ""); TRestSensitivity(); ~TRestSensitivity(); diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index 688f59dab..036c3ea7e 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -134,7 +134,7 @@ void TRestSensitivity::GenerateCurves(Int_t N) { for (int n = 0; n < N; n++) GenerateCurve(); } -std::vector TRestSensitivity::GetSensitivityCurve(size_t n) { +std::vector TRestSensitivity::GetCurve(size_t n) { if (n >= GetNumberOfCurves()) { RESTWarning << "Requested curve number : " << n << " but only " << GetNumberOfCurves() << " generated" << RESTendl; @@ -195,7 +195,7 @@ void TRestSensitivity::ExportAveragedCurve(std::string fname) { void TRestSensitivity::ExportCurve(std::string fname, int n = 0) { - std::vector curve = GetSensitivityCurve(n); + std::vector curve = GetCurve(n); if (curve.empty()) return; // Open a file for writing @@ -353,6 +353,87 @@ void TRestSensitivity::InitFromConfigFile() { Initialize(); } +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); + // 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; +} + ///////////////////////////////////////////// /// \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 From 24d772e21ce5d8cd1a1cb6620327f2e1e3a9b9c9 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sun, 7 Apr 2024 10:49:04 +0200 Subject: [PATCH 15/22] TRestSensitivity::GetLevelCurves method added --- .../sensitivity/inc/TRestSensitivity.h | 1 + .../sensitivity/src/TRestSensitivity.cxx | 42 +++++++++++++++++++ 2 files changed, 43 insertions(+) diff --git a/source/framework/sensitivity/inc/TRestSensitivity.h b/source/framework/sensitivity/inc/TRestSensitivity.h index 0c79c404e..db689cf01 100644 --- a/source/framework/sensitivity/inc/TRestSensitivity.h +++ b/source/framework/sensitivity/inc/TRestSensitivity.h @@ -66,6 +66,7 @@ class TRestSensitivity : public TRestMetadata { std::vector GetCurve(size_t n = 0); std::vector GetAveragedCurve(); + std::vector> GetLevelCurves(const std::vector& levels); void ExportCurve(std::string fname, int n); void ExportAveragedCurve(std::string fname); diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index 036c3ea7e..8d8c42755 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -353,6 +353,48 @@ 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) { From 749881ebc0ac8ded64e97bc2ad13d0dd9baddc20 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Mon, 8 Apr 2024 11:11:13 +0200 Subject: [PATCH 16/22] TRestSensitivity::DrawLevelCurves method added --- .../sensitivity/inc/TRestSensitivity.h | 1 + .../sensitivity/src/TRestSensitivity.cxx | 90 ++++++++++++++++++- 2 files changed, 90 insertions(+), 1 deletion(-) diff --git a/source/framework/sensitivity/inc/TRestSensitivity.h b/source/framework/sensitivity/inc/TRestSensitivity.h index db689cf01..38b1845a5 100644 --- a/source/framework/sensitivity/inc/TRestSensitivity.h +++ b/source/framework/sensitivity/inc/TRestSensitivity.h @@ -90,6 +90,7 @@ class TRestSensitivity : public TRestMetadata { void PrintMetadata() override; TCanvas* DrawCurves(); + TCanvas* DrawLevelCurves(); TRestSensitivity(const char* cfgFileName, const std::string& name = ""); TRestSensitivity(); diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index 03ce2a53a..daecdac5e 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -446,7 +446,8 @@ TCanvas* TRestSensitivity::DrawCurves() { 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()->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"); @@ -476,6 +477,93 @@ 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 From b592589bf01825379f65a2e6bdb1f7f9687b5670 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 8 Apr 2024 09:12:05 +0000 Subject: [PATCH 17/22] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- source/framework/sensitivity/src/TRestSensitivity.cxx | 3 --- 1 file changed, 3 deletions(-) diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index daecdac5e..b7b548729 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -194,7 +194,6 @@ void TRestSensitivity::ExportAveragedCurve(std::string fname) { } void TRestSensitivity::ExportCurve(std::string fname, int n = 0) { - std::vector curve = GetCurve(n); if (curve.empty()) return; @@ -396,7 +395,6 @@ std::vector> TRestSensitivity::GetLevelCurves(const std::v } TCanvas* TRestSensitivity::DrawCurves() { - if (fCanvas != NULL) { delete fCanvas; fCanvas = NULL; @@ -478,7 +476,6 @@ legend->Draw(); } TCanvas* TRestSensitivity::DrawLevelCurves() { - if (fCanvas != NULL) { delete fCanvas; fCanvas = NULL; From bc42645da03cadba75321a68743cfc39b35ec440 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Fri, 3 May 2024 17:31:36 +0200 Subject: [PATCH 18/22] Updating axionlib to version 2.3 --- source/libraries/axion | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From b8c414ec104d9f26b43e326dab05e3dfd4c13a14 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Fri, 3 May 2024 17:36:25 +0200 Subject: [PATCH 19/22] Updated TRestVersion.h to v2.4.3 --- source/framework/core/inc/TRestVersion.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) 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 From b6fd87f0e0e3dd125a0fa8e8dfbbce410d376fb9 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Fri, 3 May 2024 17:39:03 +0200 Subject: [PATCH 20/22] Updated generateVersionHeader.py info --- scripts/generateVersionHeader.py | 44 ++++++++++++++++---------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/scripts/generateVersionHeader.py b/scripts/generateVersionHeader.py index 29cc21712..96d0c774c 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" + "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) +#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" -) +#) +#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) From 6209c634317dadf9c4afc1a007ed1ca61ee1165f Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Fri, 3 May 2024 17:44:16 +0200 Subject: [PATCH 21/22] generateVersionHeader.py fixing typo --- scripts/generateVersionHeader.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/generateVersionHeader.py b/scripts/generateVersionHeader.py index 96d0c774c..ab65d4cae 100755 --- a/scripts/generateVersionHeader.py +++ b/scripts/generateVersionHeader.py @@ -174,7 +174,7 @@ # + str(b) # + "." # + str(c) - + '"\n' +# + '"\n' #) #print( # "And push the changes to repository. You should also push your branch to GitHub if you have not already.\n" From b310063baf37098f5b89ac7f635f39b219a97dc1 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 3 May 2024 15:44:46 +0000 Subject: [PATCH 22/22] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- scripts/generateVersionHeader.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/scripts/generateVersionHeader.py b/scripts/generateVersionHeader.py index ab65d4cae..1d85a9b38 100755 --- a/scripts/generateVersionHeader.py +++ b/scripts/generateVersionHeader.py @@ -159,9 +159,9 @@ ) print("\n") print( - "Once your PR has been accepted and merged, you should generate a new Git tag and RELEASE at the master branch.\n" + "Once your PR has been accepted and merged, you should generate a new Git tag and RELEASE at the master branch.\n" ) -#print( +# print( # "-----> git tag -a v" # + str(a) # + "." @@ -175,13 +175,13 @@ # + "." # + str(c) # + '"\n' -#) -#print( +# ) +# 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( +# ) +# 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)