Merge pull request #521 from rest-for-physics/jgalan_bIAXO_example

Adding IAXO sensitivity examples and updating sensitivity classes when necessary

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<Double_t> 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<Double_t> GetParameterizationNodes() { return fParameterizationNodes; }
 
     std::vector<std::string> GetVariables() const { return fVariables; }
@@ -121,6 +139,8 @@ class TRestComponent : public TRestMetadata {
 
     Double_t GetRawRate(std::vector<Double_t> point);
     Double_t GetTotalRate();
+    Double_t GetMaxRate();
+    Double_t GetAllNodesIntegratedRate();
     Double_t GetNormalizedRate(std::vector<Double_t> point);
     Double_t GetRate(std::vector<Double_t> point);
 
@@ -171,6 +191,6 @@ class TRestComponent : public TRestMetadata {
     TRestComponent();
     ~TRestComponent();
 
-    ClassDefOverride(TRestComponent, 5);
+    ClassDefOverride(TRestComponent, 6);
 };
 #endif

source/framework/sensitivity/inc/TRestComponentFormula.h

@@ -35,7 +35,7 @@ class TRestComponentFormula : public TRestComponent {
     std::vector<TFormula> 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;

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

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<Double_t> 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

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<Double_t>& curve) { fCurves.push_back(curve); }
+    void ImportCurve(const std::vector<Double_t>& curve) { AddCurve(curve); }
     void GenerateCurve();
     void GenerateCurves(Int_t N);
 

source/framework/sensitivity/src/TRestComponent.cxx

@@ -86,13 +86,23 @@ TRestComponent::~TRestComponent() {}
 void TRestComponent::Initialize() {
     //   SetSectionName(this->ClassName());
 
+    /// Avoiding double initialization
+    if (!fNodeDensity.empty() && fRandom) return;
+
     if (!fRandom) {
         delete fRandom;
         fRandom = nullptr;
     }
 
     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<Double_t> 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<Double_t> 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<Double_t> 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,14 +634,26 @@ void TRestComponent::PrintMetadata() {
         }
     }
 
-    if (!fParameter.empty()) {
+    if (!fParameterizationNodes.empty()) {
         RESTMetadata << " " << RESTendl;
         RESTMetadata << " === Parameterization === " << RESTendl;
         RESTMetadata << "- Parameter : " << fParameter << RESTendl;
 
         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) {

source/framework/sensitivity/src/TRestComponentFormula.cxx

@@ -131,7 +131,7 @@ Double_t TRestComponentFormula::GetFormulaRate(std::vector<Double_t> 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<Double_t> point) {
 void TRestComponentFormula::FillHistograms() {
     if (fFormulas.empty()) return;
 
+    if (GetDimensions() == 0) {
+        RESTError << "TRestComponentFormula::FillHistograms. No variables defined!" << RESTendl;
+        RESTError << "Did you add a <cVariable entry?" << RESTendl;
+        return;
+    }
+
     RESTInfo << "Generating N-dim histogram for " << GetName() << RESTendl;
 
     TString hName = "formula";
@@ -208,7 +214,7 @@ void TRestComponentFormula::PrintMetadata() {
     TRestComponent::PrintMetadata();
 
     RESTMetadata << " " << RESTendl;
-    RESTMetadata << "Formula units: " << fUnits << RESTendl;
+    RESTMetadata << "Formula units: " << fFormulaUnits << RESTendl;
 
     if (!fFormulas.empty()) {
         RESTMetadata << " " << RESTendl;
@@ -231,6 +237,9 @@ void TRestComponentFormula::InitFromConfigFile() {
 
     if (!fFormulas.empty()) return;
 
+    /// For some reason I need to do this manually. Dont understand why!
+    fFormulaUnits = GetParameter("formulaUnits");
+
     auto ele = GetElement("formula");
     while (ele != nullptr) {
         std::string name = GetParameter("name", ele, "");