From 5cf58ae5f8f318e1c323139a8d6f4d72ec146c53 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Fri, 14 Jun 2024 19:14:43 +0200 Subject: [PATCH 01/16] TRestSpiderMask::GetNumberOfArms method added --- source/framework/masks/inc/TRestSpiderMask.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/source/framework/masks/inc/TRestSpiderMask.h b/source/framework/masks/inc/TRestSpiderMask.h index 83dd9c295..87b86bad9 100644 --- a/source/framework/masks/inc/TRestSpiderMask.h +++ b/source/framework/masks/inc/TRestSpiderMask.h @@ -58,6 +58,9 @@ class TRestSpiderMask : public TRestPatternMask { /// It returns the angular width of each spider arm in radians Double_t GetArmsWidth() { return fArmsWidth; } + /// It returns the number of arms in the spider structure + size_t GetNumberOfArms() { if( fArmsSeparationAngle == 0 ) return 0; return (Int_t) (2*TMath::Pi()/fArmsSeparationAngle); } + /// It returns the inner ring radius that defines the inner start of the spider structure Double_t GetInitialRadius() { return fInitialRadius; } From be0271dc252b683906e166f0052976ad75d6d683 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 18 Jun 2024 09:29:43 +0200 Subject: [PATCH 02/16] TRestSpiderMask. Now spider generation occurs after RML loading --- source/framework/masks/inc/TRestSpiderMask.h | 3 ++- source/framework/masks/src/TRestSpiderMask.cxx | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/source/framework/masks/inc/TRestSpiderMask.h b/source/framework/masks/inc/TRestSpiderMask.h index 87b86bad9..c8ea03b5e 100644 --- a/source/framework/masks/inc/TRestSpiderMask.h +++ b/source/framework/masks/inc/TRestSpiderMask.h @@ -28,7 +28,6 @@ /// A class used to define and generate a spider structure mask class TRestSpiderMask : public TRestPatternMask { private: - void Initialize() override; /// The angle between two consecutive spider arms measured in radians. Double_t fArmsSeparationAngle = 0; //< @@ -48,6 +47,8 @@ class TRestSpiderMask : public TRestPatternMask { std::vector> fNegativeRanges; //! public: + void Initialize() override; + void GenerateSpider(); virtual Int_t GetRegion(Double_t& x, Double_t& y) override; diff --git a/source/framework/masks/src/TRestSpiderMask.cxx b/source/framework/masks/src/TRestSpiderMask.cxx index 877fed7a3..8f1bebef2 100644 --- a/source/framework/masks/src/TRestSpiderMask.cxx +++ b/source/framework/masks/src/TRestSpiderMask.cxx @@ -139,11 +139,12 @@ TRestSpiderMask::TRestSpiderMask() : TRestPatternMask() { Initialize(); } /// corresponding TRestSpiderMask section inside the RML. /// TRestSpiderMask::TRestSpiderMask(const char* cfgFileName, std::string name) : TRestPatternMask(cfgFileName) { - Initialize(); LoadConfigFromFile(fConfigFileName, name); if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); + + Initialize(); } /////////////////////////////////////////////// From 62e13497d3e0a0a9ff199ba2c854b1348b3d6247 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 18 Jun 2024 09:31:49 +0200 Subject: [PATCH 03/16] TRestPatternMask::DrawMonteCarlo updated points size and aspect ratio --- source/framework/masks/src/TRestPatternMask.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/source/framework/masks/src/TRestPatternMask.cxx b/source/framework/masks/src/TRestPatternMask.cxx index 5388888b3..51868bbc5 100644 --- a/source/framework/masks/src/TRestPatternMask.cxx +++ b/source/framework/masks/src/TRestPatternMask.cxx @@ -168,6 +168,7 @@ TCanvas* TRestPatternMask::DrawMonteCarlo(Int_t nSamples) { fCanvas = NULL; } fCanvas = new TCanvas("canv", "This is the canvas title", 1200, 1200); + fCanvas->SetRealAspectRatio(); fCanvas->Draw(); TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97); @@ -213,12 +214,12 @@ TCanvas* TRestPatternMask::DrawMonteCarlo(Int_t nSamples) { if (nGraphs == 0) { gr->SetLineColor(kBlack); gr->SetMarkerColor(kBlack); - gr->SetMarkerSize(0.6); + gr->SetMarkerSize(0.3); gr->SetLineWidth(2); } else { gr->SetMarkerColor((nGraphs * 3) % 29 + 20); gr->SetLineColor((nGraphs * 3) % 29 + 20); - gr->SetMarkerSize(0.6); + gr->SetMarkerSize(0.3); gr->SetLineWidth(2); } From e37b6473af1d02cb595aeeb7a21cac06b861b5b8 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 18 Jun 2024 09:33:53 +0200 Subject: [PATCH 04/16] TRestSpiderMask::GetNumberOfArms minor update --- source/framework/masks/inc/TRestSpiderMask.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/source/framework/masks/inc/TRestSpiderMask.h b/source/framework/masks/inc/TRestSpiderMask.h index c8ea03b5e..ee8ca76de 100644 --- a/source/framework/masks/inc/TRestSpiderMask.h +++ b/source/framework/masks/inc/TRestSpiderMask.h @@ -59,8 +59,11 @@ class TRestSpiderMask : public TRestPatternMask { /// It returns the angular width of each spider arm in radians Double_t GetArmsWidth() { return fArmsWidth; } - /// It returns the number of arms in the spider structure - size_t GetNumberOfArms() { if( fArmsSeparationAngle == 0 ) return 0; return (Int_t) (2*TMath::Pi()/fArmsSeparationAngle); } + /// It returns the number of arms in the spider structure + size_t GetNumberOfArms() { + if (fArmsSeparationAngle == 0) return 0; + return (size_t)(2 * TMath::Pi() / fArmsSeparationAngle); + } /// It returns the inner ring radius that defines the inner start of the spider structure Double_t GetInitialRadius() { return fInitialRadius; } From 8a98a91847910a89611582421ecc7cfb24841e06 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 18 Jun 2024 09:35:01 +0200 Subject: [PATCH 05/16] TRestExperiment. Fixing typos --- .../sensitivity/src/TRestExperiment.cxx | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/source/framework/sensitivity/src/TRestExperiment.cxx b/source/framework/sensitivity/src/TRestExperiment.cxx index b3e410756..ac5dbd799 100644 --- a/source/framework/sensitivity/src/TRestExperiment.cxx +++ b/source/framework/sensitivity/src/TRestExperiment.cxx @@ -131,8 +131,7 @@ 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; } @@ -141,8 +140,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; } @@ -173,7 +171,7 @@ void TRestExperiment::InitFromConfigFile() { auto ele = GetElement("addComponent"); if (ele != nullptr) { - std::string filename = GetParameter("filename", ele, ""); + std::string filename = SearchFile(GetParameter("filename", ele, "")); std::string component = GetParameter("component", ele, ""); if (filename.empty()) @@ -221,7 +219,7 @@ void TRestExperiment::InitFromConfigFile() { if (!fSignal) { RESTError << "TRestExperiment : " << GetName() << RESTendl; - RESTError << "There was a problem initiazing the signal component!" << RESTendl; + RESTError << "There was a problem initializing the signal component!" << RESTendl; return; } @@ -272,9 +270,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; From ba09216a3d2aa673eccbe1660f5dacaf2529561a Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 18 Jun 2024 09:37:37 +0200 Subject: [PATCH 06/16] TRestSensitivity::UnbinnedLogLikelihood avoiding unnecessary warning --- source/framework/sensitivity/src/TRestSensitivity.cxx | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index b9dc66893..32d0dde1b 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -191,7 +191,7 @@ void TRestSensitivity::ExportAveragedCurve(std::string fname, Double_t factor, D int m = 0; for (const auto& node : fParameterizationNodes) { - outputFile << node << " " << factor * TMath::Power(curve[m], power) << std::endl; + outputFile << node << " " << factor* TMath::Power(curve[m], power) << std::endl; m++; } @@ -221,7 +221,7 @@ void TRestSensitivity::ExportCurve(std::string fname, Double_t factor, Double_t int m = 0; for (const auto& node : fParameterizationNodes) { - outputFile << node << " " << factor * TMath::Power(curve[m], power) << std::endl; + outputFile << node << " " << factor* TMath::Power(curve[m], power) << std::endl; m++; } @@ -236,6 +236,7 @@ void TRestSensitivity::ExportCurve(std::string fname, Double_t factor, Double_t 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, node, 0); } @@ -275,8 +276,6 @@ Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experime if (nd >= 0) experiment->GetSignal()->SetActiveNode(nd); else { - RESTWarning << "Node : " << node << " not found in signal : " << experiment->GetSignal()->GetName() - << RESTendl; return 0.0; } From 4b3c0c34884d7d1ec4878e0c8e4f0c8a4187b9fb Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Wed, 19 Jun 2024 11:21:22 +0200 Subject: [PATCH 07/16] TRestComponent::fCachedRates added to speed up rate evaluation --- .../sensitivity/inc/TRestComponent.h | 5 +++- .../sensitivity/src/TRestComponent.cxx | 29 +++++++++++-------- 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/source/framework/sensitivity/inc/TRestComponent.h b/source/framework/sensitivity/inc/TRestComponent.h index 1e22722a6..2065fc222 100644 --- a/source/framework/sensitivity/inc/TRestComponent.h +++ b/source/framework/sensitivity/inc/TRestComponent.h @@ -64,6 +64,9 @@ class TRestComponent : public TRestMetadata { /// It defines the increasing step for automatic parameter list generation Double_t fStepParameterValue = 0; //< + /// It keeps track of total rates for quick access + std::vector fCachedRates; //< + /// It true the parametric values automatically generated will grow exponentially Bool_t fExponential = false; //< @@ -191,6 +194,6 @@ class TRestComponent : public TRestMetadata { TRestComponent(); ~TRestComponent(); - ClassDefOverride(TRestComponent, 6); + ClassDefOverride(TRestComponent, 7); }; #endif diff --git a/source/framework/sensitivity/src/TRestComponent.cxx b/source/framework/sensitivity/src/TRestComponent.cxx index c723bd374..dc81ae90a 100644 --- a/source/framework/sensitivity/src/TRestComponent.cxx +++ b/source/framework/sensitivity/src/TRestComponent.cxx @@ -62,7 +62,7 @@ TRestComponent::TRestComponent() {} /// 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 cfgFileName A const char* giving the path to an RML file. /// \param name The name of the specific metadata. It will be used to find the /// corresponding TRestAxionMagneticField section inside the RML. /// @@ -103,6 +103,9 @@ void TRestComponent::Initialize() { } else { if (!fParameterizationNodes.empty()) FillHistograms(); } + + fCachedRates.clear(); + if (fNodeDensity.size() > 0) fCachedRates.resize(fNodeDensity.size(), 0.0); } ///////////////////////////////////////////// @@ -325,6 +328,10 @@ Double_t TRestComponent::GetRawRate(std::vector point) { /// The result will be returned in s-1. /// Double_t TRestComponent::GetTotalRate() { + if (fCachedRates.size() == 0 && fNodeDensity.size() > 0) fCachedRates.resize(fNodeDensity.size(), 0.0); + if (fActiveNode >= 0 && (Int_t)fCachedRates.size() > fActiveNode && fCachedRates[fActiveNode] > 0) + return fCachedRates[fActiveNode]; + THnD* dHist = GetDensityForActiveNode(); if (!dHist) return 0; @@ -343,6 +350,8 @@ Double_t TRestComponent::GetTotalRate() { if (!skip) integral += GetRate(point); } + if ((Int_t)fCachedRates.size() > fActiveNode) fCachedRates[fActiveNode] = integral; + return integral; } @@ -389,8 +398,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; } @@ -420,11 +428,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"); @@ -450,15 +457,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; } From 15ec422b6adb07abb04c671a9540bac4dbfcc62f Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sat, 22 Jun 2024 21:10:14 +0200 Subject: [PATCH 08/16] TRestSensitivity. Adding Write with sSaveExperiments --- .../sensitivity/inc/TRestSensitivity.h | 13 ++++++++++--- .../sensitivity/src/TRestSensitivity.cxx | 19 ++++++++++++++++--- 2 files changed, 26 insertions(+), 6 deletions(-) diff --git a/source/framework/sensitivity/inc/TRestSensitivity.h b/source/framework/sensitivity/inc/TRestSensitivity.h index c60350f9e..c2bf4901e 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; //< @@ -37,6 +37,9 @@ class TRestSensitivity : public TRestMetadata { /// A vector of calculated sensitivity curves defined as a funtion of the parametric node std::vector> fCurves; //< + /// If disabled the experiments will not be saved to disk + Bool_t fSaveExperiments = false; //< + /// 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 @@ -66,6 +69,7 @@ class TRestSensitivity : public TRestMetadata { void GenerateCurves(Int_t N); std::vector GetCurve(size_t n = 0); + std::vector ExportCurve(size_t n = 0) { return GetCurve(n); } std::vector GetAveragedCurve(); std::vector> GetLevelCurves(const std::vector& levels); @@ -74,7 +78,8 @@ class TRestSensitivity : public TRestMetadata { TH1D* SignalStatisticalTest(Double_t node, Int_t N); - void Freeze() { fFrozen = true; } + void Freeze(Bool_t freeze = true) { fFrozen = freeze; } + void SaveExperiments(Bool_t save = true) { fSaveExperiments = save; } std::vector GetExperiments() { return fExperiments; } TRestExperiment* GetExperiment(const size_t& n) { @@ -90,6 +95,8 @@ class TRestSensitivity : public TRestMetadata { void PrintMetadata() override; + virtual Int_t Write(const char* name = nullptr, Int_t option = 0, Int_t bufsize = 0) override; + TCanvas* DrawCurves(); TCanvas* DrawLevelCurves(); @@ -97,6 +104,6 @@ class TRestSensitivity : public TRestMetadata { TRestSensitivity(); ~TRestSensitivity(); - ClassDefOverride(TRestSensitivity, 2); + ClassDefOverride(TRestSensitivity, 4); }; #endif diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index 32d0dde1b..127920ff3 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -76,7 +76,10 @@ 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 @@ -543,7 +546,7 @@ TCanvas* TRestSensitivity::DrawLevelCurves() { graphs[0]->Draw("AL"); TGraph* randomGr = new TGraph(); - std::vector randomCurve = fCurves[13]; + std::vector randomCurve = fCurves[GetNumberOfCurves() / 2]; for (size_t m = 0; m < randomCurve.size(); m++) randomGr->SetPoint(randomGr->GetN(), fParameterizationNodes[m], TMath::Sqrt(TMath::Sqrt(randomCurve[m]))); @@ -571,7 +574,7 @@ TCanvas* TRestSensitivity::DrawLevelCurves() { for (unsigned int n = 1; n < graphs.size(); n++) graphs[n]->Draw("Lsame"); randomGr->Draw("LPsame"); - // centralGr->Draw("Lsame"); + // centralGr->Draw("Lsame"); return fCanvas; } @@ -618,3 +621,13 @@ void TRestSensitivity::PrintMetadata() { RESTMetadata << "----" << RESTendl; } + +Int_t TRestSensitivity::Write(const char* name, Int_t option, Int_t bufsize) { + if (!fSaveExperiments) { + RESTInfo << "TRestSensitivity::Write. Removing experiments before writting to disk." << RESTendl; + RESTInfo << "Use TRestSensitivity::SaveExperiments( true ) to change this behaviour" << RESTendl; + fExperiments.clear(); + } + + return TRestMetadata::Write(name, option, bufsize); +} From 479da14deb19d96e922468a79184f0a4b7eda621 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Wed, 3 Jul 2024 11:59:16 +0200 Subject: [PATCH 09/16] TRestResponse::GetOverallEfficiency method added --- .../framework/sensitivity/inc/TRestResponse.h | 11 ++++++++++ .../sensitivity/src/TRestResponse.cxx | 22 +++++++++++++++++++ 2 files changed, 33 insertions(+) diff --git a/source/framework/sensitivity/inc/TRestResponse.h b/source/framework/sensitivity/inc/TRestResponse.h index 25bdb4ed1..bf509f677 100644 --- a/source/framework/sensitivity/inc/TRestResponse.h +++ b/source/framework/sensitivity/inc/TRestResponse.h @@ -64,17 +64,28 @@ class TRestResponse : public TRestMetadata { std::string GetVariable() const { return fVariable; } TVector2 GetInputRange() const { + if( fResponseMatrix.empty() ) { + RESTError << " TRestResponse::GetInputRange. Response matrix not loaded yet, try using LoadResponse" << RESTendl; + return 0; + } return TVector2(fOrigin.X(), fOrigin.X() + fResponseMatrix[0].size() * fBinSize); } TVector2 GetOutputRange() const { + if( fResponseMatrix.empty() ) { + RESTError << " TRestResponse::GetOutputRange. Response matrix not loaded yet, try using LoadResponse" << RESTendl; + return 0; + } return TVector2(fOrigin.Y(), fOrigin.Y() + fResponseMatrix.size() * fBinSize); } + void Initialize() override; void LoadResponse(Bool_t transpose = true); + Double_t GetOverallEfficiency(Double_t from, Double_t to); + std::vector> GetResponse(Double_t input); void PrintResponseMatrix(Int_t fromRow, Int_t toRow); diff --git a/source/framework/sensitivity/src/TRestResponse.cxx b/source/framework/sensitivity/src/TRestResponse.cxx index 1c9c769bf..8fab5af2e 100644 --- a/source/framework/sensitivity/src/TRestResponse.cxx +++ b/source/framework/sensitivity/src/TRestResponse.cxx @@ -129,6 +129,28 @@ void TRestResponse::LoadResponse(Bool_t transpose) { RESTError << "Extension format - " << extension << " - not recognized!" << RESTendl; } +/////////////////////////////////////////////// +/// \brief It is used to integrate the response matrix and obtain the overall efficiency +/// inside an energy range given by argument. +/// +Double_t TRestResponse::GetOverallEfficiency(Double_t from, Double_t to) { + if( fResponseMatrix.empty() ) { + RESTError << " TRestResponse::GetOverallEfficiency. Response matrix not loaded yet, try using LoadResponse" << RESTendl; + return 0; + } + + Int_t nBins = 0; + Double_t eff = 0; + for( double en = from; en < to; en += fBinSize ) + { + std::vector> effs = GetResponse( en ); + for( const auto &ef: effs ) + eff += ef.second; + nBins++; + } + return eff/nBins; +} + ///////////////////////////////////////////// /// \brief This method will return a vector of std::pair, each pair will contain the /// output energy together with the corresponding response (or efficiency), for the From 56dce397dbf5458bc28f8cf4ceb51f310063da20 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 13 Jul 2024 16:53:08 +0000 Subject: [PATCH 10/16] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../framework/masks/src/TRestPatternMask.cxx | 2 +- .../framework/masks/src/TRestSpiderMask.cxx | 1 - .../framework/sensitivity/inc/TRestResponse.h | 23 +++++++------- .../sensitivity/src/TRestComponent.cxx | 18 ++++++----- .../sensitivity/src/TRestExperiment.cxx | 12 ++++---- .../sensitivity/src/TRestResponse.cxx | 30 +++++++++---------- .../sensitivity/src/TRestSensitivity.cxx | 5 ++-- 7 files changed, 49 insertions(+), 42 deletions(-) diff --git a/source/framework/masks/src/TRestPatternMask.cxx b/source/framework/masks/src/TRestPatternMask.cxx index 51868bbc5..1eedb76f6 100644 --- a/source/framework/masks/src/TRestPatternMask.cxx +++ b/source/framework/masks/src/TRestPatternMask.cxx @@ -168,7 +168,7 @@ TCanvas* TRestPatternMask::DrawMonteCarlo(Int_t nSamples) { fCanvas = NULL; } fCanvas = new TCanvas("canv", "This is the canvas title", 1200, 1200); - fCanvas->SetRealAspectRatio(); + fCanvas->SetRealAspectRatio(); fCanvas->Draw(); TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97); diff --git a/source/framework/masks/src/TRestSpiderMask.cxx b/source/framework/masks/src/TRestSpiderMask.cxx index 8f1bebef2..ec096ccca 100644 --- a/source/framework/masks/src/TRestSpiderMask.cxx +++ b/source/framework/masks/src/TRestSpiderMask.cxx @@ -139,7 +139,6 @@ TRestSpiderMask::TRestSpiderMask() : TRestPatternMask() { Initialize(); } /// corresponding TRestSpiderMask section inside the RML. /// TRestSpiderMask::TRestSpiderMask(const char* cfgFileName, std::string name) : TRestPatternMask(cfgFileName) { - LoadConfigFromFile(fConfigFileName, name); if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); diff --git a/source/framework/sensitivity/inc/TRestResponse.h b/source/framework/sensitivity/inc/TRestResponse.h index bf509f677..1260cb186 100644 --- a/source/framework/sensitivity/inc/TRestResponse.h +++ b/source/framework/sensitivity/inc/TRestResponse.h @@ -64,27 +64,30 @@ class TRestResponse : public TRestMetadata { std::string GetVariable() const { return fVariable; } TVector2 GetInputRange() const { - if( fResponseMatrix.empty() ) { - RESTError << " TRestResponse::GetInputRange. Response matrix not loaded yet, try using LoadResponse" << RESTendl; - return 0; - } + if (fResponseMatrix.empty()) { + RESTError + << " TRestResponse::GetInputRange. Response matrix not loaded yet, try using LoadResponse" + << RESTendl; + return 0; + } return TVector2(fOrigin.X(), fOrigin.X() + fResponseMatrix[0].size() * fBinSize); } TVector2 GetOutputRange() const { - if( fResponseMatrix.empty() ) { - RESTError << " TRestResponse::GetOutputRange. Response matrix not loaded yet, try using LoadResponse" << RESTendl; - return 0; - } + if (fResponseMatrix.empty()) { + RESTError + << " TRestResponse::GetOutputRange. Response matrix not loaded yet, try using LoadResponse" + << RESTendl; + return 0; + } return TVector2(fOrigin.Y(), fOrigin.Y() + fResponseMatrix.size() * fBinSize); } - void Initialize() override; void LoadResponse(Bool_t transpose = true); - Double_t GetOverallEfficiency(Double_t from, Double_t to); + Double_t GetOverallEfficiency(Double_t from, Double_t to); std::vector> GetResponse(Double_t input); diff --git a/source/framework/sensitivity/src/TRestComponent.cxx b/source/framework/sensitivity/src/TRestComponent.cxx index dc81ae90a..b53bab64f 100644 --- a/source/framework/sensitivity/src/TRestComponent.cxx +++ b/source/framework/sensitivity/src/TRestComponent.cxx @@ -398,7 +398,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; } @@ -428,10 +429,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"); @@ -457,13 +459,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 ac5dbd799..982b48e4a 100644 --- a/source/framework/sensitivity/src/TRestExperiment.cxx +++ b/source/framework/sensitivity/src/TRestExperiment.cxx @@ -131,7 +131,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; } @@ -140,7 +141,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; } @@ -270,9 +272,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/TRestResponse.cxx b/source/framework/sensitivity/src/TRestResponse.cxx index 8fab5af2e..d915602e9 100644 --- a/source/framework/sensitivity/src/TRestResponse.cxx +++ b/source/framework/sensitivity/src/TRestResponse.cxx @@ -134,21 +134,21 @@ void TRestResponse::LoadResponse(Bool_t transpose) { /// inside an energy range given by argument. /// Double_t TRestResponse::GetOverallEfficiency(Double_t from, Double_t to) { - if( fResponseMatrix.empty() ) { - RESTError << " TRestResponse::GetOverallEfficiency. Response matrix not loaded yet, try using LoadResponse" << RESTendl; - return 0; - } - - Int_t nBins = 0; - Double_t eff = 0; - for( double en = from; en < to; en += fBinSize ) - { - std::vector> effs = GetResponse( en ); - for( const auto &ef: effs ) - eff += ef.second; - nBins++; - } - return eff/nBins; + if (fResponseMatrix.empty()) { + RESTError + << " TRestResponse::GetOverallEfficiency. Response matrix not loaded yet, try using LoadResponse" + << RESTendl; + return 0; + } + + Int_t nBins = 0; + Double_t eff = 0; + for (double en = from; en < to; en += fBinSize) { + std::vector> effs = GetResponse(en); + for (const auto& ef : effs) eff += ef.second; + nBins++; + } + return eff / nBins; } ///////////////////////////////////////////// diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index 127920ff3..cee6f4755 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -194,7 +194,7 @@ void TRestSensitivity::ExportAveragedCurve(std::string fname, Double_t factor, D int m = 0; for (const auto& node : fParameterizationNodes) { - outputFile << node << " " << factor* TMath::Power(curve[m], power) << std::endl; + outputFile << node << " " << factor * TMath::Power(curve[m], power) << std::endl; m++; } @@ -224,7 +224,7 @@ void TRestSensitivity::ExportCurve(std::string fname, Double_t factor, Double_t int m = 0; for (const auto& node : fParameterizationNodes) { - outputFile << node << " " << factor* TMath::Power(curve[m], power) << std::endl; + outputFile << node << " " << factor * TMath::Power(curve[m], power) << std::endl; m++; } @@ -239,7 +239,6 @@ void TRestSensitivity::ExportCurve(std::string fname, Double_t factor, Double_t 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, node, 0); } From 5ce4d3ba21f03a2888ed9c701b9b9c2c2fff91fd Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 16 Jul 2024 21:03:34 +0200 Subject: [PATCH 11/16] TRestComponentDataSet::Initialize. Updating initialization order --- source/framework/sensitivity/src/TRestComponentDataSet.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/framework/sensitivity/src/TRestComponentDataSet.cxx b/source/framework/sensitivity/src/TRestComponentDataSet.cxx index 94d1ef7bc..bd382193d 100644 --- a/source/framework/sensitivity/src/TRestComponentDataSet.cxx +++ b/source/framework/sensitivity/src/TRestComponentDataSet.cxx @@ -126,11 +126,11 @@ TRestComponentDataSet::TRestComponentDataSet(const char* cfgFileName, const std: /// (or observables) that have been defined by the user. /// void TRestComponentDataSet::Initialize() { - TRestComponent::Initialize(); - SetSectionName(this->ClassName()); LoadDataSets(); + + TRestComponent::Initialize(); } ///////////////////////////////////////////// From c3453762734f0005346465b6c952d0a1c9724ed2 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 16 Jul 2024 19:07:18 +0000 Subject: [PATCH 12/16] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- source/framework/sensitivity/src/TRestComponentDataSet.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/framework/sensitivity/src/TRestComponentDataSet.cxx b/source/framework/sensitivity/src/TRestComponentDataSet.cxx index bd382193d..7d53df62a 100644 --- a/source/framework/sensitivity/src/TRestComponentDataSet.cxx +++ b/source/framework/sensitivity/src/TRestComponentDataSet.cxx @@ -129,7 +129,7 @@ void TRestComponentDataSet::Initialize() { SetSectionName(this->ClassName()); LoadDataSets(); - + TRestComponent::Initialize(); } From a876b4935a35f5f3215dfb5e507781a82295e803 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sat, 27 Jul 2024 16:50:31 +0200 Subject: [PATCH 13/16] TRestComponent::SetActiveNode. Protecting --- source/framework/sensitivity/inc/TRestComponent.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/framework/sensitivity/inc/TRestComponent.h b/source/framework/sensitivity/inc/TRestComponent.h index 2065fc222..14628cf60 100644 --- a/source/framework/sensitivity/inc/TRestComponent.h +++ b/source/framework/sensitivity/inc/TRestComponent.h @@ -120,6 +120,7 @@ class TRestComponent : public TRestMetadata { /// It returns true if any nodes have been defined. Bool_t HasNodes() { return !fParameterizationNodes.empty(); } + size_t GetNumberOfParameterizationNodes() { return fParameterizationNodes.size();} virtual void RegenerateActiveNodeDensity() {} @@ -154,7 +155,7 @@ class TRestComponent : public TRestMetadata { Int_t FindActiveNode(Double_t node); Int_t SetActiveNode(Double_t node); Int_t SetActiveNode(Int_t n) { - fActiveNode = n; + fActiveNode = n%fParameterizationNodes.size(); return fActiveNode; } From 000f1d79c319c14de615d7a8667610ebc51c79c5 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sat, 27 Jul 2024 16:54:17 +0200 Subject: [PATCH 14/16] TRestSensitivity::SignalStatisticalTest now loops over all nodes --- .../framework/sensitivity/src/TRestSensitivity.cxx | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index cee6f4755..eea41cdfe 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -327,8 +327,17 @@ Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experime /// TH1D* TRestSensitivity::SignalStatisticalTest(Double_t node, Int_t N) { std::vector couplings; + Int_t nodeCheck = 0; + if( fExperiments.size() > 0 ) + nodeCheck = fExperiments[0]->GetSignal()->GetActiveNode(); for (int n = 0; n < N; n++) { - for (const auto& exp : fExperiments) exp->GetSignal()->RegenerateActiveNodeDensity(); + for (const auto& exp : fExperiments) + { + if( exp->GetSignal()->GetActiveNode() != nodeCheck ) + RESTError << "TRestSensitivity::SignalStatisticalTest. Problem" << RESTendl; + exp->GetSignal()->SetActiveNode(exp->GetSignal()->GetActiveNode()+1); + exp->GetSignal()->RegenerateActiveNodeDensity(); + } Double_t coupling = TMath::Sqrt(TMath::Sqrt(GetCoupling(node))); couplings.push_back(coupling); @@ -339,7 +348,7 @@ TH1D* TRestSensitivity::SignalStatisticalTest(Double_t node, Int_t N) { double max_value = *std::max_element(couplings.begin(), couplings.end()); if (fSignalTest) delete fSignalTest; - fSignalTest = new TH1D("SignalTest", "A signal test", 100, 0.9 * min_value, 1.1 * max_value); + fSignalTest = new TH1D("SignalTest", "A signal test", 1000, 0.9 * min_value, 1.1 * max_value); for (const auto& coup : couplings) fSignalTest->Fill(coup); return fSignalTest; From e177ae0aade136a4607ce278df5c25ce96e303cf Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 27 Jul 2024 14:54:59 +0000 Subject: [PATCH 15/16] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../framework/sensitivity/inc/TRestComponent.h | 4 ++-- .../sensitivity/src/TRestSensitivity.cxx | 18 ++++++++---------- 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/source/framework/sensitivity/inc/TRestComponent.h b/source/framework/sensitivity/inc/TRestComponent.h index 14628cf60..a122f66fd 100644 --- a/source/framework/sensitivity/inc/TRestComponent.h +++ b/source/framework/sensitivity/inc/TRestComponent.h @@ -120,7 +120,7 @@ class TRestComponent : public TRestMetadata { /// It returns true if any nodes have been defined. Bool_t HasNodes() { return !fParameterizationNodes.empty(); } - size_t GetNumberOfParameterizationNodes() { return fParameterizationNodes.size();} + size_t GetNumberOfParameterizationNodes() { return fParameterizationNodes.size(); } virtual void RegenerateActiveNodeDensity() {} @@ -155,7 +155,7 @@ class TRestComponent : public TRestMetadata { Int_t FindActiveNode(Double_t node); Int_t SetActiveNode(Double_t node); Int_t SetActiveNode(Int_t n) { - fActiveNode = n%fParameterizationNodes.size(); + fActiveNode = n % fParameterizationNodes.size(); return fActiveNode; } diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index eea41cdfe..7aeccaa35 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -327,17 +327,15 @@ Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experime /// TH1D* TRestSensitivity::SignalStatisticalTest(Double_t node, Int_t N) { std::vector couplings; - Int_t nodeCheck = 0; - if( fExperiments.size() > 0 ) - nodeCheck = fExperiments[0]->GetSignal()->GetActiveNode(); + Int_t nodeCheck = 0; + if (fExperiments.size() > 0) nodeCheck = fExperiments[0]->GetSignal()->GetActiveNode(); for (int n = 0; n < N; n++) { - for (const auto& exp : fExperiments) - { - if( exp->GetSignal()->GetActiveNode() != nodeCheck ) - RESTError << "TRestSensitivity::SignalStatisticalTest. Problem" << RESTendl; - exp->GetSignal()->SetActiveNode(exp->GetSignal()->GetActiveNode()+1); - exp->GetSignal()->RegenerateActiveNodeDensity(); - } + for (const auto& exp : fExperiments) { + if (exp->GetSignal()->GetActiveNode() != nodeCheck) + RESTError << "TRestSensitivity::SignalStatisticalTest. Problem" << RESTendl; + exp->GetSignal()->SetActiveNode(exp->GetSignal()->GetActiveNode() + 1); + exp->GetSignal()->RegenerateActiveNodeDensity(); + } Double_t coupling = TMath::Sqrt(TMath::Sqrt(GetCoupling(node))); couplings.push_back(coupling); From ce6538a320e75a14d5bcf7273dc8ab81eb8ee669 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sat, 27 Jul 2024 17:12:46 +0200 Subject: [PATCH 16/16] TRestSensitivity::DrawLevelCurves esthetical updates --- source/framework/sensitivity/src/TRestSensitivity.cxx | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index 7aeccaa35..d70c53a40 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -504,11 +504,12 @@ TCanvas* TRestSensitivity::DrawLevelCurves() { delete fCanvas; fCanvas = NULL; } - fCanvas = new TCanvas("canv", "This is the canvas title", 500, 400); + fCanvas = new TCanvas("canv", "This is the canvas title", 600, 500); fCanvas->Draw(); fCanvas->SetLeftMargin(0.15); fCanvas->SetRightMargin(0.04); fCanvas->SetLogx(); + fCanvas->SetLogy(); std::vector> levelCurves = GetLevelCurves({0.025, 0.16, 0.375, 0.625, 0.84, 0.975}); @@ -534,9 +535,9 @@ TCanvas* TRestSensitivity::DrawLevelCurves() { 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]->GetYaxis()->SetRangeUser(0, 0.31); + graphs[0]->GetXaxis()->SetRangeUser(0.001, 0.31); + graphs[0]->GetXaxis()->SetLimits(0.0001, 0.31); graphs[0]->GetXaxis()->SetTitle("mass [eV]"); graphs[0]->GetXaxis()->SetTitleSize(0.04); graphs[0]->GetXaxis()->SetTitleOffset(1.15); @@ -547,6 +548,7 @@ TCanvas* TRestSensitivity::DrawLevelCurves() { graphs[0]->GetYaxis()->SetTitleOffset(1.5); graphs[0]->GetYaxis()->SetTitleSize(0.04); graphs[0]->GetYaxis()->SetLabelSize(0.04); + graphs[0]->GetYaxis()->SetRangeUser(0.1, 30); // graphs[0]->GetYaxis()->SetLabelOffset(0); // graphs[0]->GetYaxis()->SetLabelFont(43); graphs[0]->Draw("AL");