diff --git a/source/framework/masks/inc/TRestSpiderMask.h b/source/framework/masks/inc/TRestSpiderMask.h index 83dd9c295..ee8ca76de 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; @@ -58,6 +59,12 @@ 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 (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; } diff --git a/source/framework/masks/src/TRestPatternMask.cxx b/source/framework/masks/src/TRestPatternMask.cxx index 5388888b3..1eedb76f6 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); } diff --git a/source/framework/masks/src/TRestSpiderMask.cxx b/source/framework/masks/src/TRestSpiderMask.cxx index 877fed7a3..ec096ccca 100644 --- a/source/framework/masks/src/TRestSpiderMask.cxx +++ b/source/framework/masks/src/TRestSpiderMask.cxx @@ -139,11 +139,11 @@ 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(); } /////////////////////////////////////////////// diff --git a/source/framework/sensitivity/inc/TRestComponent.h b/source/framework/sensitivity/inc/TRestComponent.h index 1e22722a6..a122f66fd 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; //< @@ -117,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() {} @@ -151,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; } @@ -191,6 +195,6 @@ class TRestComponent : public TRestMetadata { TRestComponent(); ~TRestComponent(); - ClassDefOverride(TRestComponent, 6); + ClassDefOverride(TRestComponent, 7); }; #endif diff --git a/source/framework/sensitivity/inc/TRestResponse.h b/source/framework/sensitivity/inc/TRestResponse.h index 25bdb4ed1..1260cb186 100644 --- a/source/framework/sensitivity/inc/TRestResponse.h +++ b/source/framework/sensitivity/inc/TRestResponse.h @@ -64,10 +64,22 @@ 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); } @@ -75,6 +87,8 @@ class TRestResponse : public TRestMetadata { 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/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/TRestComponent.cxx b/source/framework/sensitivity/src/TRestComponent.cxx index c723bd374..b53bab64f 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; } diff --git a/source/framework/sensitivity/src/TRestComponentDataSet.cxx b/source/framework/sensitivity/src/TRestComponentDataSet.cxx index 94d1ef7bc..7d53df62a 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(); } ///////////////////////////////////////////// diff --git a/source/framework/sensitivity/src/TRestExperiment.cxx b/source/framework/sensitivity/src/TRestExperiment.cxx index b3e410756..982b48e4a 100644 --- a/source/framework/sensitivity/src/TRestExperiment.cxx +++ b/source/framework/sensitivity/src/TRestExperiment.cxx @@ -173,7 +173,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 +221,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; } diff --git a/source/framework/sensitivity/src/TRestResponse.cxx b/source/framework/sensitivity/src/TRestResponse.cxx index 1c9c769bf..d915602e9 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 diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index b9dc66893..d70c53a40 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 @@ -275,8 +278,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; } @@ -326,8 +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(); 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); @@ -338,7 +346,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; @@ -496,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}); @@ -526,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); @@ -539,12 +548,13 @@ 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"); 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]))); @@ -572,7 +582,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; } @@ -619,3 +629,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); +}