Skip to content

Commit 8b97d34

Browse files
authored
Merge pull request #521 from rest-for-physics/jgalan_bIAXO_example
Adding IAXO sensitivity examples and updating sensitivity classes when necessary
2 parents 41c2174 + 513988e commit 8b97d34

File tree

10 files changed

+209
-66
lines changed

10 files changed

+209
-66
lines changed

source/framework/sensitivity/inc/TRestComponent.h

Lines changed: 22 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -55,6 +55,18 @@ class TRestComponent : public TRestMetadata {
5555
/// It defines the nodes of the parameterization (Initialized by the dataset)
5656
std::vector<Double_t> fParameterizationNodes; //<
5757

58+
/// It defines the first parametric node value in case of automatic parameter generation
59+
Double_t fFirstParameterValue = 0; //<
60+
61+
/// It defines the upper limit for the automatic parametric node values generation
62+
Double_t fLastParameterValue = 0; //<
63+
64+
/// It defines the increasing step for automatic parameter list generation
65+
Double_t fStepParameterValue = 0; //<
66+
67+
/// It true the parametric values automatically generated will grow exponentially
68+
Bool_t fExponential = false; //<
69+
5870
/// It is used to define the node that will be accessed for rate retrieval
5971
Int_t fActiveNode = -1; //<
6072

@@ -101,6 +113,8 @@ class TRestComponent : public TRestMetadata {
101113
void Initialize() override;
102114
void RegenerateHistograms(UInt_t seed = 0);
103115

116+
void RegenerateParametricNodes(Double_t from, Double_t to, Double_t step, Bool_t expIncrease = false);
117+
104118
/// It returns true if any nodes have been defined.
105119
Bool_t HasNodes() { return !fParameterizationNodes.empty(); }
106120

@@ -112,7 +126,11 @@ class TRestComponent : public TRestMetadata {
112126
size_t GetDimensions() { return fVariables.size(); }
113127
Int_t GetSamples() { return fSamples; }
114128
Int_t GetActiveNode() { return fActiveNode; }
115-
Double_t GetActiveNodeValue() { return fParameterizationNodes[fActiveNode]; }
129+
Double_t GetActiveNodeValue() {
130+
if (fActiveNode >= 0 && fActiveNode < (Int_t)fParameterizationNodes.size())
131+
return fParameterizationNodes[fActiveNode];
132+
return 0;
133+
}
116134
std::vector<Double_t> GetParameterizationNodes() { return fParameterizationNodes; }
117135

118136
std::vector<std::string> GetVariables() const { return fVariables; }
@@ -121,6 +139,8 @@ class TRestComponent : public TRestMetadata {
121139

122140
Double_t GetRawRate(std::vector<Double_t> point);
123141
Double_t GetTotalRate();
142+
Double_t GetMaxRate();
143+
Double_t GetAllNodesIntegratedRate();
124144
Double_t GetNormalizedRate(std::vector<Double_t> point);
125145
Double_t GetRate(std::vector<Double_t> point);
126146

@@ -171,6 +191,6 @@ class TRestComponent : public TRestMetadata {
171191
TRestComponent();
172192
~TRestComponent();
173193

174-
ClassDefOverride(TRestComponent, 5);
194+
ClassDefOverride(TRestComponent, 6);
175195
};
176196
#endif

source/framework/sensitivity/inc/TRestComponentFormula.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@ class TRestComponentFormula : public TRestComponent {
3535
std::vector<TFormula> fFormulas;
3636

3737
/// The formulas should be expressed in the following units
38-
std::string fUnits = "cm^-2*keV^-1"; //<
38+
std::string fFormulaUnits = "cm^-2*keV^-1"; //<
3939

4040
protected:
4141
void InitFromConfigFile() override;

source/framework/sensitivity/inc/TRestExperiment.h

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,9 @@ class TRestExperiment : public TRestMetadata {
5353
/// Only if it is true we will be able to calculate the LogLikelihood
5454
Bool_t fDataReady = false; //<
5555

56+
/// The mock dataset will be generated using the mean counts instead of a real MonteCarlo
57+
Bool_t fUseAverage = false; //<
58+
5659
/// It keeps track on the number of counts inside the dataset
5760
Int_t fExperimentalCounts = 0; //<
5861

@@ -66,7 +69,7 @@ class TRestExperiment : public TRestMetadata {
6669
void InitFromConfigFile() override;
6770

6871
public:
69-
void GenerateMockDataSet();
72+
void GenerateMockDataSet(Bool_t useAverage = false);
7073
Int_t GetExperimentalCounts() const { return fExperimentalCounts; }
7174

7275
Bool_t IsMockData() const { return fMockData; }
@@ -94,6 +97,6 @@ class TRestExperiment : public TRestMetadata {
9497
TRestExperiment();
9598
~TRestExperiment();
9699

97-
ClassDefOverride(TRestExperiment, 1);
100+
ClassDefOverride(TRestExperiment, 2);
98101
};
99102
#endif

source/framework/sensitivity/inc/TRestExperimentList.h

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -53,11 +53,17 @@ class TRestExperimentList : public TRestMetadata {
5353
/// The fusioned list of parameterization nodes found at each experiment signal
5454
std::vector<Double_t> fParameterizationNodes; //<
5555

56+
/// The mock dataset will be generated using the mean counts instead of a real MonteCarlo
57+
Bool_t fUseAverage = false; //<
58+
5659
/// If not zero this will be the common exposure time in micro-seconds (standard REST units)
5760
Double_t fExposureTime = 0;
5861

59-
/// In case an exposure time is given it defines how to assign time to each experiment (unique/ksvz).
60-
std::string fExposureStrategy = "unique";
62+
/// In case an exposure time is given it defines how to assign time to each experiment (equal/ksvz).
63+
std::string fExposureStrategy = "equal";
64+
65+
/// The factor used on the exponential exposure time as a function of the experiment number
66+
Double_t fExposureFactor = 0;
6167

6268
/// If not null this will be the common signal used in each experiment
6369
TRestComponent* fSignal = nullptr; //<
@@ -98,6 +104,6 @@ class TRestExperimentList : public TRestMetadata {
98104
TRestExperimentList();
99105
~TRestExperimentList();
100106

101-
ClassDefOverride(TRestExperimentList, 1);
107+
ClassDefOverride(TRestExperimentList, 2);
102108
};
103109
#endif

source/framework/sensitivity/inc/TRestSensitivity.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,7 @@ class TRestSensitivity : public TRestMetadata {
6161

6262
Double_t GetCoupling(Double_t node, Double_t sigma = 2, Double_t precision = 0.01);
6363
void AddCurve(const std::vector<Double_t>& curve) { fCurves.push_back(curve); }
64+
void ImportCurve(const std::vector<Double_t>& curve) { AddCurve(curve); }
6465
void GenerateCurve();
6566
void GenerateCurves(Int_t N);
6667

source/framework/sensitivity/src/TRestComponent.cxx

Lines changed: 92 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -86,13 +86,23 @@ TRestComponent::~TRestComponent() {}
8686
void TRestComponent::Initialize() {
8787
// SetSectionName(this->ClassName());
8888

89+
/// Avoiding double initialization
90+
if (!fNodeDensity.empty() && fRandom) return;
91+
8992
if (!fRandom) {
9093
delete fRandom;
9194
fRandom = nullptr;
9295
}
9396

9497
fRandom = new TRandom3(fSeed);
9598
fSeed = fRandom->TRandom::GetSeed();
99+
100+
if (fStepParameterValue > 0) {
101+
RegenerateParametricNodes(fFirstParameterValue, fLastParameterValue, fStepParameterValue,
102+
fExponential);
103+
} else {
104+
if (!fParameterizationNodes.empty()) FillHistograms();
105+
}
96106
}
97107

98108
/////////////////////////////////////////////
@@ -105,11 +115,33 @@ void TRestComponent::RegenerateHistograms(UInt_t seed) {
105115
fNodeDensity.clear();
106116

107117
fSeed = seed;
108-
TRestComponent::Initialize();
109-
110118
FillHistograms();
111119
}
112120

121+
/////////////////////////////////////////////
122+
/// \brief It allows to produce a parameter nodes list providing the initial
123+
/// value, the final value and the step. We might chose the step growing in
124+
/// linear increase steps or exponential. Linear is the default value.
125+
///
126+
void TRestComponent::RegenerateParametricNodes(Double_t from, Double_t to, Double_t step,
127+
Bool_t expIncrease) {
128+
fStepParameterValue = step;
129+
fFirstParameterValue = from;
130+
fLastParameterValue = to;
131+
fExponential = expIncrease;
132+
133+
fParameterizationNodes.clear();
134+
135+
if (expIncrease) {
136+
for (double p = from; p < to; p *= step) fParameterizationNodes.push_back(p);
137+
} else {
138+
for (double p = from; p < to; p += step) fParameterizationNodes.push_back(p);
139+
}
140+
141+
if (fParameterizationNodes.empty()) return;
142+
RegenerateHistograms(fSeed);
143+
}
144+
113145
///////////////////////////////////////////
114146
/// \brief It returns the position of the fVariable element for the variable
115147
/// name given by argument.
@@ -232,6 +264,11 @@ Double_t TRestComponent::GetRawRate(std::vector<Double_t> point) {
232264
return 0;
233265
}
234266

267+
for (size_t n = 0; n < point.size(); n++) {
268+
// The point is outside boundaries
269+
if (point[n] < fRanges[n].X() || point[n] > fRanges[n].Y()) return 0;
270+
}
271+
235272
Int_t centerBin[GetDimensions()];
236273
Double_t centralDensity = GetDensity()->GetBinContent(GetDensity()->GetBin(point.data()), centerBin);
237274
if (!Interpolation()) return centralDensity;
@@ -289,17 +326,53 @@ Double_t TRestComponent::GetRawRate(std::vector<Double_t> point) {
289326
///
290327
Double_t TRestComponent::GetTotalRate() {
291328
THnD* dHist = GetDensityForActiveNode();
329+
if (!dHist) return 0;
292330

293331
Double_t integral = 0;
294-
if (dHist != nullptr) {
295-
TH1D* h1 = dHist->Projection(0);
296-
integral = h1->Integral();
297-
delete h1;
332+
for (Int_t n = 0; n < dHist->GetNbins(); ++n) {
333+
Int_t centerBin[GetDimensions()];
334+
std::vector<Double_t> point;
335+
336+
dHist->GetBinContent(n, centerBin);
337+
for (size_t d = 0; d < GetDimensions(); ++d) point.push_back(GetBinCenter(d, centerBin[d]));
338+
339+
Bool_t skip = false;
340+
for (size_t d = 0; d < GetDimensions(); ++d) {
341+
if (point[d] < fRanges[d].X() || point[d] > fRanges[d].Y()) skip = true;
342+
}
343+
if (!skip) integral += GetRate(point);
298344
}
299345

300346
return integral;
301347
}
302348

349+
///////////////////////////////////////////////
350+
/// \brief This method returns the total rate for the node that has the highest contribution
351+
/// The result will be returned in s-1.
352+
///
353+
Double_t TRestComponent::GetMaxRate() {
354+
Double_t maxRate = 0;
355+
for (size_t n = 0; n < fParameterizationNodes.size(); n++) {
356+
SetActiveNode((Int_t)n);
357+
Double_t rate = GetTotalRate();
358+
if (rate > maxRate) maxRate = rate;
359+
}
360+
return maxRate;
361+
}
362+
363+
///////////////////////////////////////////////
364+
/// \brief This method returns the integrated total rate for all the nodes
365+
/// The result will be returned in s-1.
366+
///
367+
Double_t TRestComponent::GetAllNodesIntegratedRate() {
368+
Double_t rate = 0;
369+
for (size_t n = 0; n < fParameterizationNodes.size(); n++) {
370+
SetActiveNode((Int_t)n);
371+
rate += GetTotalRate();
372+
}
373+
return rate;
374+
}
375+
303376
///////////////////////////////////////////////
304377
/// \brief It returns the bin center of the given component dimension.
305378
///
@@ -561,14 +634,26 @@ void TRestComponent::PrintMetadata() {
561634
}
562635
}
563636

564-
if (!fParameter.empty()) {
637+
if (!fParameterizationNodes.empty()) {
565638
RESTMetadata << " " << RESTendl;
566639
RESTMetadata << " === Parameterization === " << RESTendl;
567640
RESTMetadata << "- Parameter : " << fParameter << RESTendl;
568641

569642
RESTMetadata << " - Number of parametric nodes : " << fParameterizationNodes.size() << RESTendl;
570643
RESTMetadata << " " << RESTendl;
571644
RESTMetadata << " Use : PrintNodes() for additional info" << RESTendl;
645+
646+
if (fStepParameterValue > 0) {
647+
RESTMetadata << " " << RESTendl;
648+
RESTMetadata << " Nodes were automatically generated using these parameters" << RESTendl;
649+
RESTMetadata << " - First node : " << fFirstParameterValue << RESTendl;
650+
RESTMetadata << " - Upper limit node : " << fLastParameterValue << RESTendl;
651+
RESTMetadata << " - Increasing step : " << fStepParameterValue << RESTendl;
652+
if (fExponential)
653+
RESTMetadata << " - Increases exponentially" << RESTendl;
654+
else
655+
RESTMetadata << " - Increases linearly" << RESTendl;
656+
}
572657
}
573658

574659
if (fResponse) {

source/framework/sensitivity/src/TRestComponentFormula.cxx

Lines changed: 11 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -131,7 +131,7 @@ Double_t TRestComponentFormula::GetFormulaRate(std::vector<Double_t> point) {
131131
normFactor *= (fRanges[n].Y() - fRanges[n].X()) / fNbins[n];
132132
}
133133

134-
return normFactor * result / units(fUnits);
134+
return normFactor * result / units(fFormulaUnits);
135135
}
136136

137137
/////////////////////////////////////////////
@@ -149,6 +149,12 @@ Double_t TRestComponentFormula::GetFormulaRate(std::vector<Double_t> point) {
149149
void TRestComponentFormula::FillHistograms() {
150150
if (fFormulas.empty()) return;
151151

152+
if (GetDimensions() == 0) {
153+
RESTError << "TRestComponentFormula::FillHistograms. No variables defined!" << RESTendl;
154+
RESTError << "Did you add a <cVariable entry?" << RESTendl;
155+
return;
156+
}
157+
152158
RESTInfo << "Generating N-dim histogram for " << GetName() << RESTendl;
153159

154160
TString hName = "formula";
@@ -208,7 +214,7 @@ void TRestComponentFormula::PrintMetadata() {
208214
TRestComponent::PrintMetadata();
209215

210216
RESTMetadata << " " << RESTendl;
211-
RESTMetadata << "Formula units: " << fUnits << RESTendl;
217+
RESTMetadata << "Formula units: " << fFormulaUnits << RESTendl;
212218

213219
if (!fFormulas.empty()) {
214220
RESTMetadata << " " << RESTendl;
@@ -231,6 +237,9 @@ void TRestComponentFormula::InitFromConfigFile() {
231237

232238
if (!fFormulas.empty()) return;
233239

240+
/// For some reason I need to do this manually. Dont understand why!
241+
fFormulaUnits = GetParameter("formulaUnits");
242+
234243
auto ele = GetElement("formula");
235244
while (ele != nullptr) {
236245
std::string name = GetParameter("name", ele, "");

0 commit comments

Comments
 (0)