diff --git a/source/framework/sensitivity/inc/TRestComponent.h b/source/framework/sensitivity/inc/TRestComponent.h index a4d8d213e..1e22722a6 100644 --- a/source/framework/sensitivity/inc/TRestComponent.h +++ b/source/framework/sensitivity/inc/TRestComponent.h @@ -55,6 +55,18 @@ class TRestComponent : public TRestMetadata { /// It defines the nodes of the parameterization (Initialized by the dataset) std::vector fParameterizationNodes; //< + /// It defines the first parametric node value in case of automatic parameter generation + Double_t fFirstParameterValue = 0; //< + + /// It defines the upper limit for the automatic parametric node values generation + Double_t fLastParameterValue = 0; //< + + /// It defines the increasing step for automatic parameter list generation + Double_t fStepParameterValue = 0; //< + + /// It true the parametric values automatically generated will grow exponentially + Bool_t fExponential = false; //< + /// It is used to define the node that will be accessed for rate retrieval Int_t fActiveNode = -1; //< @@ -101,6 +113,8 @@ class TRestComponent : public TRestMetadata { void Initialize() override; void RegenerateHistograms(UInt_t seed = 0); + void RegenerateParametricNodes(Double_t from, Double_t to, Double_t step, Bool_t expIncrease = false); + /// It returns true if any nodes have been defined. Bool_t HasNodes() { return !fParameterizationNodes.empty(); } @@ -112,7 +126,11 @@ class TRestComponent : public TRestMetadata { size_t GetDimensions() { return fVariables.size(); } Int_t GetSamples() { return fSamples; } Int_t GetActiveNode() { return fActiveNode; } - Double_t GetActiveNodeValue() { return fParameterizationNodes[fActiveNode]; } + Double_t GetActiveNodeValue() { + if (fActiveNode >= 0 && fActiveNode < (Int_t)fParameterizationNodes.size()) + return fParameterizationNodes[fActiveNode]; + return 0; + } std::vector GetParameterizationNodes() { return fParameterizationNodes; } std::vector GetVariables() const { return fVariables; } @@ -121,6 +139,8 @@ class TRestComponent : public TRestMetadata { Double_t GetRawRate(std::vector point); Double_t GetTotalRate(); + Double_t GetMaxRate(); + Double_t GetAllNodesIntegratedRate(); Double_t GetNormalizedRate(std::vector point); Double_t GetRate(std::vector point); @@ -171,6 +191,6 @@ class TRestComponent : public TRestMetadata { TRestComponent(); ~TRestComponent(); - ClassDefOverride(TRestComponent, 5); + ClassDefOverride(TRestComponent, 6); }; #endif diff --git a/source/framework/sensitivity/inc/TRestComponentFormula.h b/source/framework/sensitivity/inc/TRestComponentFormula.h index 9f9046824..1913c5f3f 100644 --- a/source/framework/sensitivity/inc/TRestComponentFormula.h +++ b/source/framework/sensitivity/inc/TRestComponentFormula.h @@ -35,7 +35,7 @@ class TRestComponentFormula : public TRestComponent { std::vector fFormulas; /// The formulas should be expressed in the following units - std::string fUnits = "cm^-2*keV^-1"; //< + std::string fFormulaUnits = "cm^-2*keV^-1"; //< protected: void InitFromConfigFile() override; diff --git a/source/framework/sensitivity/inc/TRestExperiment.h b/source/framework/sensitivity/inc/TRestExperiment.h index ee9efd923..440ec3512 100644 --- a/source/framework/sensitivity/inc/TRestExperiment.h +++ b/source/framework/sensitivity/inc/TRestExperiment.h @@ -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; //< + /// The mock dataset will be generated using the mean counts instead of a real MonteCarlo + Bool_t fUseAverage = false; //< + /// It keeps track on the number of counts inside the dataset Int_t fExperimentalCounts = 0; //< @@ -66,7 +69,7 @@ class TRestExperiment : public TRestMetadata { void InitFromConfigFile() override; public: - void GenerateMockDataSet(); + void GenerateMockDataSet(Bool_t useAverage = false); Int_t GetExperimentalCounts() const { return fExperimentalCounts; } Bool_t IsMockData() const { return fMockData; } @@ -94,6 +97,6 @@ class TRestExperiment : public TRestMetadata { TRestExperiment(); ~TRestExperiment(); - ClassDefOverride(TRestExperiment, 1); + ClassDefOverride(TRestExperiment, 2); }; #endif diff --git a/source/framework/sensitivity/inc/TRestExperimentList.h b/source/framework/sensitivity/inc/TRestExperimentList.h index 2c99e57a4..ca0b60e02 100644 --- a/source/framework/sensitivity/inc/TRestExperimentList.h +++ b/source/framework/sensitivity/inc/TRestExperimentList.h @@ -53,11 +53,17 @@ class TRestExperimentList : public TRestMetadata { /// The fusioned list of parameterization nodes found at each experiment signal std::vector fParameterizationNodes; //< + /// The mock dataset will be generated using the mean counts instead of a real MonteCarlo + Bool_t fUseAverage = false; //< + /// If not zero this will be the common exposure time in micro-seconds (standard REST units) Double_t fExposureTime = 0; - /// In case an exposure time is given it defines how to assign time to each experiment (unique/ksvz). - std::string fExposureStrategy = "unique"; + /// In case an exposure time is given it defines how to assign time to each experiment (equal/ksvz). + std::string fExposureStrategy = "equal"; + + /// The factor used on the exponential exposure time as a function of the experiment number + Double_t fExposureFactor = 0; /// If not null this will be the common signal used in each experiment TRestComponent* fSignal = nullptr; //< @@ -98,6 +104,6 @@ class TRestExperimentList : public TRestMetadata { TRestExperimentList(); ~TRestExperimentList(); - ClassDefOverride(TRestExperimentList, 1); + ClassDefOverride(TRestExperimentList, 2); }; #endif diff --git a/source/framework/sensitivity/inc/TRestSensitivity.h b/source/framework/sensitivity/inc/TRestSensitivity.h index 38b1845a5..a6ab4169c 100644 --- a/source/framework/sensitivity/inc/TRestSensitivity.h +++ b/source/framework/sensitivity/inc/TRestSensitivity.h @@ -61,6 +61,7 @@ class TRestSensitivity : public TRestMetadata { Double_t GetCoupling(Double_t node, Double_t sigma = 2, Double_t precision = 0.01); void AddCurve(const std::vector& curve) { fCurves.push_back(curve); } + void ImportCurve(const std::vector& curve) { AddCurve(curve); } void GenerateCurve(); void GenerateCurves(Int_t N); diff --git a/source/framework/sensitivity/src/TRestComponent.cxx b/source/framework/sensitivity/src/TRestComponent.cxx index 75906a898..c723bd374 100644 --- a/source/framework/sensitivity/src/TRestComponent.cxx +++ b/source/framework/sensitivity/src/TRestComponent.cxx @@ -86,6 +86,9 @@ TRestComponent::~TRestComponent() {} void TRestComponent::Initialize() { // SetSectionName(this->ClassName()); + /// Avoiding double initialization + if (!fNodeDensity.empty() && fRandom) return; + if (!fRandom) { delete fRandom; fRandom = nullptr; @@ -93,6 +96,13 @@ void TRestComponent::Initialize() { fRandom = new TRandom3(fSeed); fSeed = fRandom->TRandom::GetSeed(); + + if (fStepParameterValue > 0) { + RegenerateParametricNodes(fFirstParameterValue, fLastParameterValue, fStepParameterValue, + fExponential); + } else { + if (!fParameterizationNodes.empty()) FillHistograms(); + } } ///////////////////////////////////////////// @@ -105,11 +115,33 @@ void TRestComponent::RegenerateHistograms(UInt_t seed) { fNodeDensity.clear(); fSeed = seed; - TRestComponent::Initialize(); - FillHistograms(); } +///////////////////////////////////////////// +/// \brief It allows to produce a parameter nodes list providing the initial +/// value, the final value and the step. We might chose the step growing in +/// linear increase steps or exponential. Linear is the default value. +/// +void TRestComponent::RegenerateParametricNodes(Double_t from, Double_t to, Double_t step, + Bool_t expIncrease) { + fStepParameterValue = step; + fFirstParameterValue = from; + fLastParameterValue = to; + fExponential = expIncrease; + + fParameterizationNodes.clear(); + + if (expIncrease) { + for (double p = from; p < to; p *= step) fParameterizationNodes.push_back(p); + } else { + for (double p = from; p < to; p += step) fParameterizationNodes.push_back(p); + } + + if (fParameterizationNodes.empty()) return; + RegenerateHistograms(fSeed); +} + /////////////////////////////////////////// /// \brief It returns the position of the fVariable element for the variable /// name given by argument. @@ -232,6 +264,11 @@ Double_t TRestComponent::GetRawRate(std::vector point) { return 0; } + for (size_t n = 0; n < point.size(); n++) { + // The point is outside boundaries + if (point[n] < fRanges[n].X() || point[n] > fRanges[n].Y()) return 0; + } + Int_t centerBin[GetDimensions()]; Double_t centralDensity = GetDensity()->GetBinContent(GetDensity()->GetBin(point.data()), centerBin); if (!Interpolation()) return centralDensity; @@ -289,17 +326,53 @@ Double_t TRestComponent::GetRawRate(std::vector point) { /// Double_t TRestComponent::GetTotalRate() { THnD* dHist = GetDensityForActiveNode(); + if (!dHist) return 0; Double_t integral = 0; - if (dHist != nullptr) { - TH1D* h1 = dHist->Projection(0); - integral = h1->Integral(); - delete h1; + for (Int_t n = 0; n < dHist->GetNbins(); ++n) { + Int_t centerBin[GetDimensions()]; + std::vector point; + + dHist->GetBinContent(n, centerBin); + for (size_t d = 0; d < GetDimensions(); ++d) point.push_back(GetBinCenter(d, centerBin[d])); + + Bool_t skip = false; + for (size_t d = 0; d < GetDimensions(); ++d) { + if (point[d] < fRanges[d].X() || point[d] > fRanges[d].Y()) skip = true; + } + if (!skip) integral += GetRate(point); } return integral; } +/////////////////////////////////////////////// +/// \brief This method returns the total rate for the node that has the highest contribution +/// The result will be returned in s-1. +/// +Double_t TRestComponent::GetMaxRate() { + Double_t maxRate = 0; + for (size_t n = 0; n < fParameterizationNodes.size(); n++) { + SetActiveNode((Int_t)n); + Double_t rate = GetTotalRate(); + if (rate > maxRate) maxRate = rate; + } + return maxRate; +} + +/////////////////////////////////////////////// +/// \brief This method returns the integrated total rate for all the nodes +/// The result will be returned in s-1. +/// +Double_t TRestComponent::GetAllNodesIntegratedRate() { + Double_t rate = 0; + for (size_t n = 0; n < fParameterizationNodes.size(); n++) { + SetActiveNode((Int_t)n); + rate += GetTotalRate(); + } + return rate; +} + /////////////////////////////////////////////// /// \brief It returns the bin center of the given component dimension. /// @@ -561,7 +634,7 @@ void TRestComponent::PrintMetadata() { } } - if (!fParameter.empty()) { + if (!fParameterizationNodes.empty()) { RESTMetadata << " " << RESTendl; RESTMetadata << " === Parameterization === " << RESTendl; RESTMetadata << "- Parameter : " << fParameter << RESTendl; @@ -569,6 +642,18 @@ void TRestComponent::PrintMetadata() { RESTMetadata << " - Number of parametric nodes : " << fParameterizationNodes.size() << RESTendl; RESTMetadata << " " << RESTendl; RESTMetadata << " Use : PrintNodes() for additional info" << RESTendl; + + if (fStepParameterValue > 0) { + RESTMetadata << " " << RESTendl; + RESTMetadata << " Nodes were automatically generated using these parameters" << RESTendl; + RESTMetadata << " - First node : " << fFirstParameterValue << RESTendl; + RESTMetadata << " - Upper limit node : " << fLastParameterValue << RESTendl; + RESTMetadata << " - Increasing step : " << fStepParameterValue << RESTendl; + if (fExponential) + RESTMetadata << " - Increases exponentially" << RESTendl; + else + RESTMetadata << " - Increases linearly" << RESTendl; + } } if (fResponse) { diff --git a/source/framework/sensitivity/src/TRestComponentFormula.cxx b/source/framework/sensitivity/src/TRestComponentFormula.cxx index 06c8ef3d6..db5d518fd 100644 --- a/source/framework/sensitivity/src/TRestComponentFormula.cxx +++ b/source/framework/sensitivity/src/TRestComponentFormula.cxx @@ -131,7 +131,7 @@ Double_t TRestComponentFormula::GetFormulaRate(std::vector point) { normFactor *= (fRanges[n].Y() - fRanges[n].X()) / fNbins[n]; } - return normFactor * result / units(fUnits); + return normFactor * result / units(fFormulaUnits); } ///////////////////////////////////////////// @@ -149,6 +149,12 @@ Double_t TRestComponentFormula::GetFormulaRate(std::vector point) { void TRestComponentFormula::FillHistograms() { if (fFormulas.empty()) return; + if (GetDimensions() == 0) { + RESTError << "TRestComponentFormula::FillHistograms. No variables defined!" << RESTendl; + RESTError << "Did you add a TRandom::GetSeed(); } -void TRestExperiment::GenerateMockDataSet() { +void TRestExperiment::GenerateMockDataSet(Bool_t useAverage) { + fUseAverage = useAverage; + if (!fBackground) { RESTError << "TRestExperiment::GenerateMockData. Background component was not initialized!" << RESTendl; @@ -102,6 +104,7 @@ void TRestExperiment::GenerateMockDataSet() { Double_t meanCounts = GetBackground()->GetTotalRate() * fExposureTime * units("s"); Int_t N = fRandom->Poisson(meanCounts); + if (fUseAverage) N = (Int_t)meanCounts; RESTInfo << "Experiment: " << GetName() << " Generating mock dataset. Counts: " << N << RESTendl; ROOT::RDF::RNode df = fBackground->GetMonteCarloDataFrame(N); @@ -153,17 +156,19 @@ void TRestExperiment::InitFromConfigFile() { TRestMetadata::InitFromConfigFile(); int cont = 0; - TRestComponent* comp = (TRestComponent*)this->InstantiateChildMetadata(cont, "Component"); - while (comp != nullptr) { - if (ToLower(comp->GetNature()) == "background") - fBackground = comp; - else if (ToLower(comp->GetNature()) == "signal") - fSignal = comp; - else - RESTWarning << "TRestExperiment::InitFromConfigFile. Unknown component!" << RESTendl; - + TRestMetadata* md = (TRestMetadata*)this->InstantiateChildMetadata(cont); + while (md != nullptr) { + if (md->InheritsFrom("TRestComponent")) { + TRestComponent* comp = (TRestComponent*)md; + if (ToLower(comp->GetNature()) == "background") + fBackground = comp; + else if (ToLower(comp->GetNature()) == "signal") + fSignal = comp; + else + RESTWarning << "TRestExperiment::InitFromConfigFile. Unknown component!" << RESTendl; + } cont++; - comp = (TRestComponent*)this->InstantiateChildMetadata(cont, "Component"); + md = (TRestMetadata*)this->InstantiateChildMetadata(cont); } auto ele = GetElement("addComponent"); @@ -200,7 +205,7 @@ void TRestExperiment::InitFromConfigFile() { } if (fExposureTime > 0 && fExperimentalDataSet.empty()) { - GenerateMockDataSet(); + GenerateMockDataSet(fUseAverage); } else if (fExposureTime == 0 && !fExperimentalDataSet.empty()) { SetExperimentalDataSet(fExperimentalDataSet); } else { @@ -278,9 +283,14 @@ void TRestExperiment::PrintMetadata() { if (fMockData) { RESTMetadata << " " << RESTendl; - if (fMockData) + if (fMockData) { RESTMetadata << "The dataset was MC-generated" << RESTendl; - else { + if (fUseAverage) + RESTMetadata + << " - The number of counts in dataset was generated with the mean background counts" + << RESTendl; + + } else { RESTMetadata << "The dataset was loaded from an existing dataset file" << RESTendl; if (!fExperimentalDataSet.empty()) RESTMetadata << " - Experimental dataset file : " << fExperimentalDataSet << RESTendl; diff --git a/source/framework/sensitivity/src/TRestExperimentList.cxx b/source/framework/sensitivity/src/TRestExperimentList.cxx index 18050b92a..b370fcb9e 100644 --- a/source/framework/sensitivity/src/TRestExperimentList.cxx +++ b/source/framework/sensitivity/src/TRestExperimentList.cxx @@ -212,7 +212,7 @@ void TRestExperimentList::InitFromConfigFile() { if (generateMockData) { RESTInfo << "TRestExperimentList. Generating mock dataset" << RESTendl; - experiment->GenerateMockDataSet(); + experiment->GenerateMockDataSet(fUseAverage); } if (experiment->GetSignal() && experiment->GetBackground()) { @@ -221,47 +221,38 @@ void TRestExperimentList::InitFromConfigFile() { } } - /// 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") { + if (fExposureTime > 0 && ToLower(fExposureStrategy) == "exp") { 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(); - } - } + Double_t sum = 0; + for (size_t n = 0; n < fExperiments.size(); n++) sum += TMath::Exp((double)n * fExposureFactor); - Coefficient += 1. / node / expectedRate; + Double_t A = fExposureTime * units("s") / sum; + for (size_t n = 0; n < fExperiments.size(); n++) { + fExperiments[n]->SetExposureInSeconds(A * TMath::Exp((double)n * fExposureFactor)); + fExperiments[n]->GenerateMockDataSet(fUseAverage); } + } - 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")); + if (fExposureTime > 0 && ToLower(fExposureStrategy) == "power") { + ExtractExperimentParameterizationNodes(); + + Double_t sum = 0; + for (size_t n = 0; n < fExperiments.size(); n++) sum += TMath::Power((double)n, fExposureFactor); + + Double_t A = fExposureTime * units("s") / sum; + for (size_t n = 0; n < fExperiments.size(); n++) { + fExperiments[n]->SetExposureInSeconds(A * TMath::Power((double)n, fExposureFactor)); + fExperiments[n]->GenerateMockDataSet(fUseAverage); } + } - Double_t totalExp = 0; - for (const auto& experiment : fExperiments) totalExp += experiment->GetExposureInSeconds(); + if (fExposureTime > 0 && ToLower(fExposureStrategy) == "equal") { + ExtractExperimentParameterizationNodes(); - for (const auto& experiment : fExperiments) { - Double_t xp = experiment->GetExposureInSeconds(); - experiment->SetExposureInSeconds(xp * fExposureTime * units("s") / totalExp); - experiment->GenerateMockDataSet(); + for (size_t n = 0; n < fExperiments.size(); n++) { + fExperiments[n]->SetExposureInSeconds(fExposureTime * units("s") / fExperiments.size()); + fExperiments[n]->GenerateMockDataSet(fUseAverage); } } } @@ -320,5 +311,7 @@ void TRestExperimentList::PrintMetadata() { RESTMetadata << "Number of experiments loaded: " << fExperiments.size() << RESTendl; + if (fUseAverage) RESTMetadata << "Average MonteCarlo counts generation was enabled" << RESTendl; + RESTMetadata << "----" << RESTendl; } diff --git a/source/framework/sensitivity/src/TRestSensitivity.cxx b/source/framework/sensitivity/src/TRestSensitivity.cxx index b7b548729..57d44be79 100644 --- a/source/framework/sensitivity/src/TRestSensitivity.cxx +++ b/source/framework/sensitivity/src/TRestSensitivity.cxx @@ -131,6 +131,13 @@ void TRestSensitivity::GenerateCurve() { } void TRestSensitivity::GenerateCurves(Int_t N) { + /* + std::cout << "TRestSensitivity::GenerateCurves." << std::endl; + std::cout << "There is a potential memory leak when generating several curves." << std::endl; + std::cout << "This code needs to be reviewed" << std::endl; + return; + */ + for (int n = 0; n < N; n++) GenerateCurve(); } @@ -256,13 +263,22 @@ Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experime return lhood; } + if (!experiment->GetSignal()->HasNodes()) { + RESTError << "Experiment signal : " << experiment->GetSignal()->GetName() << " has no nodes!" + << RESTendl; + return lhood; + } + /// We check if the signal component is sensitive to that particular node /// If not, this experiment will not contribute to that node Int_t nd = experiment->GetSignal()->FindActiveNode(node); if (nd >= 0) experiment->GetSignal()->SetActiveNode(nd); - else + else { + RESTWarning << "Node : " << node << " not found in signal : " << experiment->GetSignal()->GetName() + << RESTendl; return 0.0; + } /// We could check if background has also components, but for the moment we do not have a background /// for each node, although it could be the case, if for example the background depends on the detector