diff --git a/data/distributions/CosmicsCry.root b/data/distributions/CosmicsCry.root
new file mode 100644
index 000000000..7f6dace14
Binary files /dev/null and b/data/distributions/CosmicsCry.root differ
diff --git a/macros/REST_OpenInputFile.C b/macros/REST_OpenInputFile.C
index 7812d8f5f..0f1af5cd6 100644
--- a/macros/REST_OpenInputFile.C
+++ b/macros/REST_OpenInputFile.C
@@ -47,7 +47,7 @@ void REST_OpenInputFile(const std::string& fileName) {
         printf("%s\n", " - dSet->GetDataFrame().GetColumnNames()");
         printf("%s\n\n", " - dSet->GetTree()->GetEntries()");
         printf("%s\n", " - dSet->GetDataFrame().Display(\"\")->Print()");
-        printf("%s\n\n", " - dSet->GetDataFrame().Display({\"colName1,colName2\"})->Print()");
+        printf("%s\n\n", " - dSet->GetDataFrame().Display({\"colName1\", \"colName2\"})->Print()");
         if (dSet) delete dSet;
         dSet = new TRestDataSet();
         dSet->EnableMultiThreading(false);
diff --git a/scripts/generateVersionHeader.py b/scripts/generateVersionHeader.py
index 29cc21712..1d85a9b38 100755
--- a/scripts/generateVersionHeader.py
+++ b/scripts/generateVersionHeader.py
@@ -159,29 +159,29 @@
 )
 print("\n")
 print(
-    "Once your PR has been accepted and merged, you should generate a new Git tag at the master branch.\n"
-)
-print(
-    "-----> git tag -a v"
-    + str(a)
-    + "."
-    + str(b)
-    + "."
-    + str(c)
-    + ' -m "Update to version v'
-    + str(a)
-    + "."
-    + str(b)
-    + "."
-    + str(c)
-    + '"\n'
-)
-print(
-    "And push the changes to repository. You should also push your branch to GitHub if you have not already.\n"
-)
-print("-----> git push origin v" + str(a) + "." + str(b) + "." + str(c) + "\n")
-print(
-    "IMPORTANT. Summarize the changes in the tag generated at the Gitlab repository website.\n"
+    "Once your PR has been accepted and merged, you should generate a new Git tag and RELEASE at the master branch.\n"
 )
+# print(
+#    "-----> git tag -a v"
+#    + str(a)
+#    + "."
+#    + str(b)
+#    + "."
+#    + str(c)
+#    + ' -m "Update to version v'
+#    + str(a)
+#    + "."
+#    + str(b)
+#    + "."
+#    + str(c)
+#    + '"\n'
+# )
+# print(
+#    "And push the changes to repository. You should also push your branch to GitHub if you have not already.\n"
+# )
+# print("-----> git push origin v" + str(a) + "." + str(b) + "." + str(c) + "\n")
+# print(
+#    "IMPORTANT. Summarize the changes in the tag generated at the Gitlab repository website.\n"
+# )
 
 exit(0)
diff --git a/source/framework/analysis/src/TRestDataSetCalibration.cxx b/source/framework/analysis/src/TRestDataSetCalibration.cxx
index d9869d29d..f6794f3a9 100644
--- a/source/framework/analysis/src/TRestDataSetCalibration.cxx
+++ b/source/framework/analysis/src/TRestDataSetCalibration.cxx
@@ -192,7 +192,7 @@ void TRestDataSetCalibration::Calibrate() {
 
     if (fCalibFile.empty()) {
         auto histo = dataSet.GetDataFrame().Histo1D(
-            {"spectrum", "spectrum", fNBins, fCalibRange.X(), fCalibRange.X()}, fCalObservable);
+            {"spectrum", "spectrum", fNBins, fCalibRange.X(), fCalibRange.Y()}, fCalObservable);
         spectrum = std::unique_ptr<TH1F>(static_cast<TH1F*>(histo->DrawClone()));
 
         // Get position of the maximum
diff --git a/source/framework/core/inc/TRestDataSetPlot.h b/source/framework/core/inc/TRestDataSetPlot.h
index fec0c6266..d957631b2 100644
--- a/source/framework/core/inc/TRestDataSetPlot.h
+++ b/source/framework/core/inc/TRestDataSetPlot.h
@@ -95,6 +95,7 @@ class TRestDataSetPlot : public TRestMetadata {
         std::vector<std::pair<std::array<std::string, 3>, TVector2>> variablePos;
         std::vector<std::pair<std::array<std::string, 3>, TVector2>> metadataPos;
         std::vector<std::pair<std::array<std::string, 3>, TVector2>> obsPos;
+        std::vector<std::pair<std::array<std::string, 3>, TVector2>> expPos;
 
         TRestCut* panelCut = nullptr;
 
@@ -182,6 +183,6 @@ class TRestDataSetPlot : public TRestMetadata {
     TRestDataSetPlot(const char* configFilename, std::string name = "");
     ~TRestDataSetPlot();
 
-    ClassDefOverride(TRestDataSetPlot, 1);
+    ClassDefOverride(TRestDataSetPlot, 2);
 };
 #endif
diff --git a/source/framework/core/inc/TRestVersion.h b/source/framework/core/inc/TRestVersion.h
index d7d1c8f66..eb622bc7f 100644
--- a/source/framework/core/inc/TRestVersion.h
+++ b/source/framework/core/inc/TRestVersion.h
@@ -12,12 +12,12 @@
  * #endif
  *
  */
-#define REST_RELEASE "2.4.2"
-#define REST_RELEASE_DATE "Mon Feb 12"
-#define REST_RELEASE_TIME "22:23:31 CET 2024"
-#define REST_RELEASE_NAME "Henry Primakoff"
-#define REST_GIT_COMMIT "d8ec95be"
-#define REST_VERSION_CODE 132098
+#define REST_RELEASE "2.4.3"
+#define REST_RELEASE_DATE "Fri May  3"
+#define REST_RELEASE_TIME "17:36:10 CEST 2024"
+#define REST_RELEASE_NAME "Steven Weinberg"
+#define REST_GIT_COMMIT "bc42645d"
+#define REST_VERSION_CODE 132099
 #define REST_VERSION(a, b, c) (((a) << 16) + ((b) << 8) + (c))
 #define REST_SCHEMA_EVOLUTION "ON"
 #endif
diff --git a/source/framework/core/src/TRestDataSet.cxx b/source/framework/core/src/TRestDataSet.cxx
index 1a3710d86..1465b37dc 100644
--- a/source/framework/core/src/TRestDataSet.cxx
+++ b/source/framework/core/src/TRestDataSet.cxx
@@ -169,6 +169,26 @@
 /// [2] d.GetTree()->GetEntries()
 /// \endcode
 ///
+/// Example 3 Automatically importing a dataset using restRoot
+///
+/// \code
+/// restRoot Dataset_FinalBabyIAXO_XMM_mM_P14.root
+///
+///  REST dataset file identified. It contains a valid TRestDataSet.
+///
+///  Importing dataset /nfs/dust/iaxo/user/jgalan//Dataset_FinalBabyIAXO_XMM_mM_P14.root as `dSet`
+///
+///  The dataset is ready. You may now access the dataset using:
+///
+///   - dSet->PrintMetadata()
+///   - dSet->GetDataFrame().GetColumnNames()
+///   - dSet->GetTree()->GetEntries()
+///
+///   - dSet->GetDataFrame().Display("")->Print()
+///   - dSet->GetDataFrame().Display({"colName1,colName2"})->Print()
+/// [0]
+/// \endcode
+///
 /// ### Relevant quantities
 ///
 /// Sometimes we will be willing that our dataset contains few variables
@@ -241,6 +261,21 @@
 /// where `SolarFlux`,`GeneratorArea` and `Nsim` are the given names of
 /// the relevant quantities inside the dataset.
 ///
+/// ### Adding cuts
+///
+/// It is also possible to add cuts used to filter the data that will
+/// be stored inside the dataset. We can do that including a TRestCut
+/// definition inside the TRestDataSet.
+///
+/// For example, the following cut definition would discard entries
+/// with unexpected values inside the specified column, `process_status`.
+///
+/// \code
+///    <TRestCut>
+///            <cut name="goodData" value="!TMath::IsNaN(process_status)"/>
+///    </TRestCut>
+/// \endcode
+///
 ///----------------------------------------------------------------------
 ///
 /// REST-for-Physics - Software for Rare Event Searches Toolkit
@@ -487,7 +522,7 @@ std::vector<std::string> TRestDataSet::FileSelection() {
 }
 
 ///////////////////////////////////////////////
-/// \brief This function apply a TRestCut to the dataframe
+/// \brief This function applies a TRestCut to the dataframe
 /// and returns a dataframe with the applied cuts. Note that
 /// the cuts are not applied directly to the dataframe on
 /// TRestDataSet, to do so you should do fDataSet = MakeCut(fCut);
diff --git a/source/framework/core/src/TRestDataSetPlot.cxx b/source/framework/core/src/TRestDataSetPlot.cxx
index ec25cf686..d65e118aa 100644
--- a/source/framework/core/src/TRestDataSetPlot.cxx
+++ b/source/framework/core/src/TRestDataSetPlot.cxx
@@ -91,8 +91,12 @@
 /// * **precision**: Precision for the values to be written.
 /// * **value**: If true/ON panel is displayed, otherwise is ignored
 /// Different keys are provided: `metadata` is meant for the metadata info inside the
-/// TRestDataSet, `variable` for a predefined variable e.g. rate and `observable` for an
-/// observable value. All the keys have the same structure which is detailed below:
+/// TRestDataSet (as a RelevantQuantity), `variable` for a predefined variable e.g. rate,
+/// `observable` for an observable value and `expression` for a mathematical expression
+/// that can contain any of the previous. Note that the time-related variables _startTime_,
+/// _endTime_ and _runLength_ (then _meanRate_ too) are obtained from the TRestDataSet and
+/// not the RDataFrame of the TRestDataSet, thus they are not afected by the cuts.
+/// All the keys have the same structure which is detailed below:
 /// * **value**: Name of the metadata, variable or observable value.
 /// * **label**: String of the label that will be written before the observable value.
 /// * **units**: String with the units to be appended to the label.
@@ -107,6 +111,10 @@
 ///          <metadata value="[TRestRun::fRunTag]" label="Run tag" x="0.25" y="0.82" />
 ///          <variable value="[[entries]]" label="Entries" x="0.25" y="0.58" />
 ///          <observable value="alphaTrackAna_angle" label="Mean Angle" units="rad" x="0.25" y="0.01" />
+///          <expression value="cos(alphaTrackAna_angle)^2" label="Cosine of the mean angle" units="" x="0.25"
+///          y="0.12" />
+///          <expression value="[TRestDetector::fDriftField]*[TRestDetector::fPressure]" label="Drift field"
+///          units="V/cm" x="0.25" y="0.24" />
 ///          <addCut name="Fiducial"/>
 ///    </panel>
 /// \endcode
@@ -124,8 +132,10 @@
 /// * **timeDisplay**: If true/ON time display is set in the X axis
 /// * **norm**: Normalization constant in which the plot will be normalized, e.g. use `1` in
 /// case you want to normalize by 1.
-/// * **scale**: Multiply all the histogram bins by a constant given by scale. If you want to
-/// use the size of the bin to normalize you should write down `binSize`
+/// * **scale**: Divide all the histogram bins by a constant given by scale. You may use any
+/// mathematical expression in combination with the special keywords: `binSize`, `entries`,
+/// `runLength` (in hours) and `integral`. Adding a scale will make the histogram to be
+/// drawn with errors. To avoid this, set parameter option to 'hist' in the histogram.
 /// * **value**: If true/ON plot is displayed, otherwise is ignored
 /// * **save**: String with the name of the output file in which the plot will be saved
 /// in a separated file, several formats are supported (root, pdf, eps, jpg,...)
@@ -252,6 +262,8 @@
 /// 2023-04: First implementation of TRestDataSetPlot, based on TRestAnalysisPlot
 /// JuanAn Garcia
 ///
+/// 2024-05: Extend some functionalities, Álvaro Ezquerro
+///
 /// \class TRestDataSetPlot
 /// \author: JuanAn Garcia   e-mail: juanangp@unizar.es
 ///
@@ -389,6 +401,19 @@ void TRestDataSetPlot::ReadPanelInfo() {
 
             observable = GetNextElement(observable);
         }
+        TiXmlElement* expression = GetElement("expression", panelele);
+        while (expression != nullptr) {
+            std::array<std::string, 3> label;
+            label[0] = GetParameter("value", expression, "");
+            label[1] = GetParameter("label", expression, "");
+            label[2] = GetParameter("units", expression, "");
+            double posX = StringToDouble(GetParameter("x", expression, "0.1"));
+            double posY = StringToDouble(GetParameter("y", expression, "0.1"));
+
+            panel.expPos.push_back(std::make_pair(label, TVector2(posX, posY)));
+
+            expression = GetNextElement(expression);
+        }
 
         fPanels.push_back(panel);
         panelele = GetNextElement(panelele);
@@ -505,7 +530,50 @@ void TRestDataSetPlot::GenerateDataSetFromFilePattern(TRestDataSet& dataSet) {
     std::map<std::string, TRestDataSet::RelevantQuantity> quantity;
 
     for (auto& panel : fPanels) {
-        // Add obserbables from panel info, both variables and cuts
+        // Add metadata and observables from expression key from panel info
+        for (auto& [key, posLabel] : panel.expPos) {
+            // look for metadata which are surrounded by [] but not [[]] (variables)
+            auto&& [exp, label, units] = key;
+            std::string text = exp;
+            while (text.find_last_of('[') != std::string::npos) {
+                int squareBracketCorrector = 0;
+                size_t posOpen = text.find_last_of('[');
+                size_t posClose = text.find_first_of(']', posOpen);
+                if (posOpen != 0) {
+                    if (text[posOpen - 1] == '[') {
+                        squareBracketCorrector = 1;
+                    }
+                }
+                std::string varOrMeta = text.substr(posOpen - squareBracketCorrector,
+                                                    posClose + 1 - posOpen + 2 * squareBracketCorrector);
+                if (squareBracketCorrector == 0) {
+                    // Here varOrMeta is a metadata
+                    // Add metadata to relevant quantity from the panel info
+                    TRestDataSet::RelevantQuantity quant;
+                    quant.metadata = varOrMeta;
+                    quant.strategy = "unique";
+                    quantity[label] = quant;
+                }
+                text = Replace(text, varOrMeta, "1");
+            }
+
+            // look for observables (characterized by having a _ in the name)
+            while (text.find("_") != std::string::npos) {
+                size_t pos_ = text.find("_");
+                size_t beginning = text.find_last_of(" -+*/)(^%", pos_) + 1;
+                size_t end = text.find_first_of(" -+*/)(^%", pos_);
+                std::string obs = text.substr(beginning, end - beginning);
+                text = Replace(text, obs, "1");
+                obsList.push_back(obs);
+            }
+
+            if (!(isAExpression(text) || isANumber(text)))
+                RESTWarning << "The expression " << exp
+                            << " has not been correctly parsed into variables, metadata and observables"
+                            << RESTendl;
+        }
+
+        // Add observables from observable key of panel info
         for (auto& [key, posLabel] : panel.obsPos) {
             auto&& [obs, label, units] = key;
             obsList.push_back(obs);
@@ -535,6 +603,7 @@ void TRestDataSetPlot::GenerateDataSetFromFilePattern(TRestDataSet& dataSet) {
 ///
 void TRestDataSetPlot::PlotCombinedCanvas() {
     TRestDataSet dataSet;
+    dataSet.EnableMultiThreading(true);
 
     // Import dataSet
     dataSet.Import(fDataSetName);
@@ -586,6 +655,41 @@ void TRestDataSetPlot::PlotCombinedCanvas() {
         paramMap["[[runLength]]"] = StringWithPrecision(runLength, panel.precision);
         paramMap["[[entries]]"] = StringWithPrecision(entries, panel.precision);
         paramMap["[[meanRate]]"] = StringWithPrecision(meanRate, panel.precision);
+
+        paramMap["[[cutNames]]"] = "";
+        paramMap["[[cuts]]"] = "";
+        if (fCut) {
+            for (const auto& cut : fCut->GetCuts()) {
+                if (paramMap["[[cutNames]]"].empty())
+                    paramMap["[[cutNames]]"] += cut.GetName();
+                else
+                    paramMap["[[cutNames]]"] += "," + (std::string)cut.GetName();
+                if (paramMap["[[cuts]]"].empty())
+                    paramMap["[[cuts]]"] += cut.GetTitle();
+                else
+                    paramMap["[[cuts]]"] += " && " + (std::string)cut.GetTitle();
+            }
+        }
+
+        paramMap["[[panelCutNames]]"] = "";
+        paramMap["[[panelCuts]]"] = "";
+        if (panel.panelCut) {
+            for (const auto& cut : panel.panelCut->GetCuts()) {
+                if (paramMap["[[panelCutNames]]"].empty())
+                    paramMap["[[panelCutNames]]"] += cut.GetName();
+                else
+                    paramMap["[[panelCutNames]]"] += "," + (std::string)cut.GetName();
+                if (paramMap["[[panelCuts]]"].empty())
+                    paramMap["[[panelCuts]]"] += cut.GetTitle();
+                else
+                    paramMap["[[panelCuts]]"] += " && " + (std::string)cut.GetTitle();
+            }
+        }
+
+        RESTInfo << "Global cuts: " << paramMap["[[cuts]]"] << RESTendl;
+        if (!paramMap["[[panelCuts]]"].empty())
+            RESTInfo << "Additional panel cuts: " << paramMap["[[panelCuts]]"] << RESTendl;
+
         // Replace panel variables and generate a TLatex label
         for (const auto& [key, posLabel] : panel.variablePos) {
             auto&& [variable, label, units] = key;
@@ -633,6 +737,35 @@ void TRestDataSetPlot::PlotCombinedCanvas() {
             panel.text.emplace_back(new TLatex(posLabel.X(), posLabel.Y(), lab.c_str()));
         }
 
+        // Replace any expression and generate TLatex label
+        for (const auto& [key, posLabel] : panel.expPos) {
+            auto&& [text, label, units] = key;
+            std::string var = text;
+
+            // replace variables
+            for (const auto& [param, val] : paramMap) {
+                var = Replace(var, param, val);
+            }
+            // replace metadata
+            for (const auto& [name, quant] : quantity) {
+                var = Replace(var, name, quant.value);
+            }
+            // replace observables
+            for (const auto& obs : dataFrame.GetColumnNames()) {
+                if (var.find(obs) == std::string::npos) continue;
+                // here there should be a checking that the mean(obs) can be calculated
+                // (checking obs data type?)
+                double value = *dataFrame.Mean(obs);
+                var = Replace(var, obs, DoubleToString(value));
+            }
+            var = Replace(var, "[", "(");
+            var = Replace(var, "]", ")");
+            var = EvaluateExpression(var);
+
+            std::string lab = label + ": " + var + " " + units;
+            panel.text.emplace_back(new TLatex(posLabel.X(), posLabel.Y(), lab.c_str()));
+        }
+
         // Draw the labels inside the pad
         for (const auto& text : panel.text) {
             text->SetTextColor(1);
@@ -691,13 +824,19 @@ void TRestDataSetPlot::PlotCombinedCanvas() {
             }
             // Scale histos
             if (plots.scale != "") {
-                Double_t scale = 1.;
-                if (plots.scale == "binSize") {
-                    scale = 1. / hist.histo->GetXaxis()->GetBinWidth(1);
-                } else {
-                    scale = StringToDouble(plots.scale);
-                }
-                hist.histo->Scale(scale);
+                std::string inputScale = plots.scale;
+                double binSize = hist.histo->GetXaxis()->GetBinWidth(1);
+                double entries = hist.histo->GetEntries();
+                double runLength = dataSet.GetTotalTimeInSeconds() / 3600.;  // in hours
+                double integral = hist.histo->Integral("width");
+
+                inputScale = Replace(inputScale, "binSize", DoubleToString(binSize));
+                inputScale = Replace(inputScale, "entries", DoubleToString(entries));
+                inputScale = Replace(inputScale, "runLength", DoubleToString(runLength));
+                inputScale = Replace(inputScale, "integral", DoubleToString(integral));
+
+                std::string scale = "1./(" + inputScale + ")";
+                hist.histo->Scale(StringToDouble(EvaluateExpression(scale)));  // -1 if 'scale' isn't valid
             }
 
             // Add histos to the THStack
@@ -881,6 +1020,12 @@ void TRestDataSetPlot::PrintMetadata() {
                          << " Pos (" << posLabel.X() << ", " << posLabel.Y() << ")" << RESTendl;
         }
         RESTMetadata << "****************" << RESTendl;
+        for (auto& [key, posLabel] : panel.expPos) {
+            auto&& [obs, label, units] = key;
+            RESTMetadata << "Label Expression " << obs << ", label " << label << ", units " << units
+                         << " Pos (" << posLabel.X() << ", " << posLabel.Y() << ")" << RESTendl;
+        }
+        RESTMetadata << "****************" << RESTendl;
     }
     RESTMetadata << "-------------------" << RESTendl;
 
diff --git a/source/framework/sensitivity/inc/TRestComponent.h b/source/framework/sensitivity/inc/TRestComponent.h
index 1dfd614dc..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<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;  //<
 
@@ -82,9 +94,6 @@ class TRestComponent : public TRestMetadata {
     /// A canvas for drawing the active node component
     TCanvas* fCanvas = nullptr;  //!
 
-    /// It returns true if any nodes have been defined.
-    Bool_t HasNodes() { return !fParameterizationNodes.empty(); }
-
     Bool_t HasDensity() { return !fNodeDensity.empty(); }
 
     /// It returns true if the node has been properly identified
@@ -104,6 +113,11 @@ 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(); }
+
     virtual void RegenerateActiveNodeDensity() {}
 
     std::string GetNature() const { return fNature; }
@@ -112,7 +126,12 @@ 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; }
     std::vector<TVector2> GetRanges() const { return fRanges; }
@@ -120,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);
 
@@ -127,6 +148,7 @@ class TRestComponent : public TRestMetadata {
 
     void SetPrecision(const Float_t& pr) { fPrecision = pr; }
 
+    Int_t FindActiveNode(Double_t node);
     Int_t SetActiveNode(Double_t node);
     Int_t SetActiveNode(Int_t n) {
         fActiveNode = n;
@@ -169,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<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;
diff --git a/source/framework/sensitivity/inc/TRestExperiment.h b/source/framework/sensitivity/inc/TRestExperiment.h
index 7e25014f2..440ec3512 100644
--- a/source/framework/sensitivity/inc/TRestExperiment.h
+++ b/source/framework/sensitivity/inc/TRestExperiment.h
@@ -42,7 +42,7 @@ class TRestExperiment : public TRestMetadata {
     TRestComponent* fSignal = nullptr;  //<
 
     /// It defines the filename used to load the dataset
-    std::string fExperimentalDataSet = "";
+    std::string fExperimentalDataSet = "";  //<
 
     /// It contains the experimental data (should contain same columns as the components)
     TRestDataSet fExperimentalData;  //<
@@ -53,6 +53,12 @@ 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;  //<
+
     /// Internal process random generator
     TRandom3* fRandom = nullptr;  //!
 
@@ -63,7 +69,8 @@ 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; }
     Bool_t IsDataReady() const { return fDataReady; }
@@ -90,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 5d1737e40..ca0b60e02 100644
--- a/source/framework/sensitivity/inc/TRestExperimentList.h
+++ b/source/framework/sensitivity/inc/TRestExperimentList.h
@@ -50,9 +50,21 @@ class TRestExperimentList : public TRestMetadata {
     /// A vector with a list of experiments includes the background components in this model
     std::vector<TRestExperiment*> fExperiments;  //<
 
+    /// 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 (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;  //<
 
@@ -71,6 +83,10 @@ class TRestExperimentList : public TRestMetadata {
     void SetSignal(TRestComponent* comp) { fSignal = comp; }
     void SetBackground(TRestComponent* comp) { fBackground = comp; }
 
+    void ExtractExperimentParameterizationNodes();
+    std::vector<Double_t> GetParameterizationNodes() { return fParameterizationNodes; }
+    void PrintParameterizationNodes();
+
     std::vector<TRestExperiment*> GetExperiments() { return fExperiments; }
     TRestExperiment* GetExperiment(const size_t& n) {
         if (n >= GetNumberOfExperiments())
@@ -88,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 56c9d0bd6..a6ab4169c 100644
--- a/source/framework/sensitivity/inc/TRestSensitivity.h
+++ b/source/framework/sensitivity/inc/TRestSensitivity.h
@@ -29,22 +29,52 @@
 class TRestSensitivity : public TRestMetadata {
    private:
     /// A list of experimental conditions included to get a final sensitivity plot
-    std::vector<TRestExperiment*> fExperiments;  //<
+    std::vector<TRestExperiment*> fExperiments;  //!
 
+    /// The fusioned list of parameterization nodes found at each experiment signal
+    std::vector<Double_t> fParameterizationNodes;  //<
+
+    /// A vector of calculated sensitivity curves defined as a funtion of the parametric node
+    std::vector<std::vector<Double_t>> fCurves;  //<
+
+    /// 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
+
+    /// It is used to generate a histogram with the signal distribution produced with different signal samples
     TH1D* fSignalTest = nullptr;
 
+    /// A canvas to draw
+    TCanvas* fCanvas = nullptr;  //!
+
    protected:
     void InitFromConfigFile() override;
 
-    Double_t UnbinnedLogLikelihood(const TRestExperiment* experiment, Double_t g4 = 0);
-    Double_t ApproachByFactor(Double_t g4, Double_t chi0, Double_t target, Double_t factor);
+    Double_t UnbinnedLogLikelihood(const TRestExperiment* experiment, Double_t node, Double_t g4 = 0);
+    Double_t ApproachByFactor(Double_t node, Double_t g4, Double_t chi0, Double_t target, Double_t factor);
 
    public:
     void Initialize() override;
 
-    Double_t GetCoupling(Double_t sigma = 2, Double_t precision = 0.01);
+    void ExtractExperimentParameterizationNodes(Bool_t rescan = false);
+    std::vector<Double_t> GetParameterizationNodes() { return fParameterizationNodes; }
+    void PrintParameterizationNodes();
+
+    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);
 
-    TH1D* SignalStatisticalTest(Int_t N);
+    std::vector<Double_t> GetCurve(size_t n = 0);
+    std::vector<Double_t> GetAveragedCurve();
+    std::vector<std::vector<Double_t>> GetLevelCurves(const std::vector<Double_t>& levels);
+
+    void ExportCurve(std::string fname, int n);
+    void ExportAveragedCurve(std::string fname);
+
+    TH1D* SignalStatisticalTest(Double_t node, Int_t N);
+
+    void Freeze() { fFrozen = true; }
 
     std::vector<TRestExperiment*> GetExperiments() { return fExperiments; }
     TRestExperiment* GetExperiment(const size_t& n) {
@@ -55,13 +85,18 @@ class TRestSensitivity : public TRestMetadata {
     }
 
     size_t GetNumberOfExperiments() { return fExperiments.size(); }
+    size_t GetNumberOfCurves() { return fCurves.size(); }
+    size_t GetNumberOfNodes() { return fParameterizationNodes.size(); }
 
     void PrintMetadata() override;
 
+    TCanvas* DrawCurves();
+    TCanvas* DrawLevelCurves();
+
     TRestSensitivity(const char* cfgFileName, const std::string& name = "");
     TRestSensitivity();
     ~TRestSensitivity();
 
-    ClassDefOverride(TRestSensitivity, 1);
+    ClassDefOverride(TRestSensitivity, 2);
 };
 #endif
diff --git a/source/framework/sensitivity/src/TRestComponent.cxx b/source/framework/sensitivity/src/TRestComponent.cxx
index 516bf8104..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<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,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) {
@@ -627,6 +712,25 @@ void TRestComponent::InitFromConfigFile() {
     if (fResponse) fResponse->LoadResponse();
 }
 
+/////////////////////////////////////////////
+/// \brief It returns the position of the fParameterizationNodes
+/// element for the variable name given by argument.
+///
+Int_t TRestComponent::FindActiveNode(Double_t node) {
+    int n = 0;
+    for (const auto& v : fParameterizationNodes) {
+        Double_t pUp = node * (1 + fPrecision / 2);
+        Double_t pDown = node * (1 - fPrecision / 2);
+        if (v > pDown && v < pUp) {
+            fActiveNode = n;
+            return fActiveNode;
+        }
+        n++;
+    }
+
+    return -1;
+}
+
 /////////////////////////////////////////////
 /// \brief It returns the position of the fParameterizationNodes
 /// element for the variable name given by argument.
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<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, "");
diff --git a/source/framework/sensitivity/src/TRestExperiment.cxx b/source/framework/sensitivity/src/TRestExperiment.cxx
index 674c391c9..b3e410756 100644
--- a/source/framework/sensitivity/src/TRestExperiment.cxx
+++ b/source/framework/sensitivity/src/TRestExperiment.cxx
@@ -87,7 +87,9 @@ void TRestExperiment::Initialize() {
     fSeed = fRandom->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,12 +104,16 @@ 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);
 
     fExperimentalData.SetDataFrame(df);
     fExperimentalData.SetTotalTimeInSeconds(fExposureTime * units("s"));
 
+    fExperimentalCounts = *fExperimentalData.GetDataFrame().Count();
+
     fMockData = true;
     fDataReady = true;
 }
@@ -118,6 +124,7 @@ void TRestExperiment::SetExperimentalDataSet(const std::string& filename) {
 
     /// fExposureTime is in standard REST units : us
     fExposureTime = fExperimentalData.GetTotalTimeInSeconds() / units("s");
+    fExperimentalCounts = *fExperimentalData.GetDataFrame().Count();
 
     fMockData = false;
     fDataReady = true;
@@ -149,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");
@@ -196,7 +205,7 @@ void TRestExperiment::InitFromConfigFile() {
     }
 
     if (fExposureTime > 0 && fExperimentalDataSet.empty()) {
-        GenerateMockDataSet();
+        GenerateMockDataSet(fUseAverage);
     } else if (fExposureTime == 0 && !fExperimentalDataSet.empty()) {
         SetExperimentalDataSet(fExperimentalDataSet);
     } else {
@@ -274,16 +283,21 @@ 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;
         }
     }
 
-    RESTMetadata << " - Experimental counts : " << *fExperimentalData.GetDataFrame().Count() << RESTendl;
+    RESTMetadata << " - Experimental counts : " << fExperimentalCounts << RESTendl;
 
     RESTMetadata << "----" << RESTendl;
 }
diff --git a/source/framework/sensitivity/src/TRestExperimentList.cxx b/source/framework/sensitivity/src/TRestExperimentList.cxx
index ab773d024..b370fcb9e 100644
--- a/source/framework/sensitivity/src/TRestExperimentList.cxx
+++ b/source/framework/sensitivity/src/TRestExperimentList.cxx
@@ -177,14 +177,17 @@ void TRestExperimentList::InitFromConfigFile() {
                 }
                 column++;
             } else {
-                experiment->SetExposureInSeconds(fExposureTime * units("s"));
-                // We will generate mock data once we load the background component
-                generateMockData = true;
+                if (ToLower(fExposureStrategy) == "unique") {
+                    experiment->SetExposureInSeconds(fExposureTime * units("s"));
+                    // We will generate mock data once we load the background component
+                    generateMockData = true;
+                }
             }
 
             if (!fSignal) {
                 if (GetComponent(experimentRow[column])) {
                     TRestComponent* sgnl = (TRestComponent*)GetComponent(experimentRow[column])->Clone();
+                    sgnl->SetName((TString)experimentRow[column]);
                     experiment->SetSignal(sgnl);
                 } else {
                     RESTError << "TRestExperimentList. Signal component : " << experimentRow[column]
@@ -209,14 +212,78 @@ void TRestExperimentList::InitFromConfigFile() {
 
             if (generateMockData) {
                 RESTInfo << "TRestExperimentList. Generating mock dataset" << RESTendl;
-                experiment->GenerateMockDataSet();
+                experiment->GenerateMockDataSet(fUseAverage);
+            }
+
+            if (experiment->GetSignal() && experiment->GetBackground()) {
+                experiment->SetName(experiment->GetSignal()->GetName());
+                fExperiments.push_back(experiment);
+            }
+        }
+
+        if (fExposureTime > 0 && ToLower(fExposureStrategy) == "exp") {
+            ExtractExperimentParameterizationNodes();
+
+            Double_t sum = 0;
+            for (size_t n = 0; n < fExperiments.size(); n++) sum += TMath::Exp((double)n * fExposureFactor);
+
+            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);
+            }
+        }
+
+        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);
             }
+        }
 
-            fExperiments.push_back(experiment);
+        if (fExposureTime > 0 && ToLower(fExposureStrategy) == "equal") {
+            ExtractExperimentParameterizationNodes();
+
+            for (size_t n = 0; n < fExperiments.size(); n++) {
+                fExperiments[n]->SetExposureInSeconds(fExposureTime * units("s") / fExperiments.size());
+                fExperiments[n]->GenerateMockDataSet(fUseAverage);
+            }
         }
     }
 }
 
+/////////////////////////////////////////////
+/// \brief It scans all the experiment signals parametric nodes to build a complete list
+/// of nodes used to build a complete sensitivity curve. Some experiments may be
+/// sensitivy to a particular node, while others may be sensitivy to another. If more
+/// than one experiment is sensitivy to a given node, the sensitivity will be combined
+/// later on.
+///
+void TRestExperimentList::ExtractExperimentParameterizationNodes() {
+    fParameterizationNodes.clear();
+
+    for (const auto& experiment : fExperiments) {
+        std::vector<Double_t> nodes = experiment->GetSignal()->GetParameterizationNodes();
+        fParameterizationNodes.insert(fParameterizationNodes.end(), nodes.begin(), nodes.end());
+
+        std::sort(fParameterizationNodes.begin(), fParameterizationNodes.end());
+        auto last = std::unique(fParameterizationNodes.begin(), fParameterizationNodes.end());
+        fParameterizationNodes.erase(last, fParameterizationNodes.end());
+    }
+}
+
+void TRestExperimentList::PrintParameterizationNodes() {
+    std::cout << "Experiment list nodes: ";
+    for (const auto& node : fParameterizationNodes) std::cout << node << "\t";
+    std::cout << std::endl;
+}
+
 TRestComponent* TRestExperimentList::GetComponent(std::string compName) {
     TRestComponent* component = nullptr;
     for (const auto& c : fComponentFiles) {
@@ -244,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 12fe749d6..57d44be79 100644
--- a/source/framework/sensitivity/src/TRestSensitivity.cxx
+++ b/source/framework/sensitivity/src/TRestSensitivity.cxx
@@ -83,7 +83,8 @@ void TRestSensitivity::Initialize() { SetSectionName(this->ClassName()); }
 /// closer to the target value given by argument. The factor will be used to
 /// increase or decrease the coupling, and evaluate the likelihood.
 ///
-Double_t TRestSensitivity::ApproachByFactor(Double_t g4, Double_t chi0, Double_t target, Double_t factor) {
+Double_t TRestSensitivity::ApproachByFactor(Double_t node, Double_t g4, Double_t chi0, Double_t target,
+                                            Double_t factor) {
     if (factor == 1) {
         return 0;
     }
@@ -92,7 +93,7 @@ Double_t TRestSensitivity::ApproachByFactor(Double_t g4, Double_t chi0, Double_t
     Double_t Chi2 = 0;
     do {
         Chi2 = 0;
-        for (const auto& exp : fExperiments) Chi2 += -2 * UnbinnedLogLikelihood(exp, g4);
+        for (const auto& exp : fExperiments) Chi2 += -2 * UnbinnedLogLikelihood(exp, node, g4);
 
         g4 = factor * g4;
     } while (Chi2 - chi0 < target);
@@ -101,7 +102,7 @@ Double_t TRestSensitivity::ApproachByFactor(Double_t g4, Double_t chi0, Double_t
     /// Coarse movement to get to Chi2 below target (/2)
     do {
         Chi2 = 0;
-        for (const auto& exp : fExperiments) Chi2 += -2 * UnbinnedLogLikelihood(exp, g4);
+        for (const auto& exp : fExperiments) Chi2 += -2 * UnbinnedLogLikelihood(exp, node, g4);
 
         g4 = g4 / factor;
     } while (Chi2 - chi0 > target);
@@ -109,21 +110,143 @@ Double_t TRestSensitivity::ApproachByFactor(Double_t g4, Double_t chi0, Double_t
     return g4 * factor;
 }
 
+void TRestSensitivity::GenerateCurve() {
+    ExtractExperimentParameterizationNodes();
+
+    if (GetNumberOfCurves() > 0)
+        for (const auto& exp : fExperiments) {
+            exp->GenerateMockDataSet();
+        }
+
+    RESTInfo << "Generating sensitivity curve" << RESTendl;
+    std::vector<Double_t> curve;
+    for (const auto& node : fParameterizationNodes) {
+        RESTInfo << "Generating node : " << node << RESTendl;
+        curve.push_back(GetCoupling(node));
+    }
+    fCurves.push_back(curve);
+
+    RESTInfo << "Curve has been generated. You may use now TRestSensitivity::ExportCurve( fname.txt )."
+             << RESTendl;
+}
+
+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();
+}
+
+std::vector<Double_t> TRestSensitivity::GetCurve(size_t n) {
+    if (n >= GetNumberOfCurves()) {
+        RESTWarning << "Requested curve number : " << n << " but only " << GetNumberOfCurves() << " generated"
+                    << RESTendl;
+        return std::vector<Double_t>();
+    }
+    return fCurves[n];
+}
+
+std::vector<Double_t> TRestSensitivity::GetAveragedCurve() {
+    if (GetNumberOfCurves() <= 0) return std::vector<Double_t>();
+
+    std::vector<double> averagedCurve(fCurves[0].size(), 0.0);  // Initialize with zeros
+
+    for (const auto& row : fCurves) {
+        for (size_t i = 0; i < row.size(); ++i) {
+            averagedCurve[i] += row[i];
+        }
+    }
+
+    for (double& avg : averagedCurve) {
+        avg /= static_cast<double>(fCurves.size());
+    }
+
+    return averagedCurve;
+}
+
+void TRestSensitivity::ExportAveragedCurve(std::string fname) {
+    std::vector<Double_t> curve = GetAveragedCurve();
+    if (curve.empty()) std::cout << "Curve is empty" << std::endl;
+    if (curve.empty()) return;
+
+    // Open a file for writing
+    std::ofstream outputFile(fname);
+
+    // Check if the file is opened successfully
+    if (!outputFile) {
+        RESTError << "TRestSensitivity::ExportCurve. Error opening file for writing!" << RESTendl;
+        return;
+    }
+
+    if (fParameterizationNodes.size() != curve.size()) {
+        RESTError << "TRestSensitivity::ExportCurve. Curve has not been properly generated" << RESTendl;
+        RESTError << "Parameterization nodes: " << fParameterizationNodes.size() << RESTendl;
+        RESTError << "Try invoking TRestSensitivity::GenerateCurve" << RESTendl;
+        return;
+    }
+
+    int m = 0;
+    for (const auto& node : fParameterizationNodes) {
+        outputFile << node << " " << curve[m] << std::endl;
+        m++;
+    }
+
+    outputFile.close();
+
+    RESTInfo << "TRestSensitivity::ExportCurve. File has been written successfully!" << RESTendl;
+}
+
+void TRestSensitivity::ExportCurve(std::string fname, int n = 0) {
+    std::vector<Double_t> curve = GetCurve(n);
+    if (curve.empty()) return;
+
+    // Open a file for writing
+    std::ofstream outputFile(fname);
+
+    // Check if the file is opened successfully
+    if (!outputFile) {
+        RESTError << "TRestSensitivity::ExportCurve. Error opening file for writing!" << RESTendl;
+        return;
+    }
+
+    if (fParameterizationNodes.size() != curve.size()) {
+        RESTError << "TRestSensitivity::ExportCurve. Curve has not been properly generated" << RESTendl;
+        RESTError << "Try invoking TRestSensitivity::GenerateCurve" << RESTendl;
+        return;
+    }
+
+    int m = 0;
+    for (const auto& node : fParameterizationNodes) {
+        outputFile << node << " " << curve[m] << std::endl;
+        m++;
+    }
+
+    outputFile.close();
+
+    RESTInfo << "TRestSensitivity::ExportCurve. File has been written successfully!" << RESTendl;
+}
+
 ///////////////////////////////////////////////
 /// \brief It will return the coupling value for which Chi=sigma
 ///
-Double_t TRestSensitivity::GetCoupling(Double_t sigma, Double_t precision) {
+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, 0);
+    for (const auto& exp : fExperiments) {
+        Chi2_0 += -2 * UnbinnedLogLikelihood(exp, node, 0);
+    }
 
     Double_t target = sigma * sigma;
 
     Double_t g4 = 0.5;
 
-    g4 = ApproachByFactor(g4, Chi2_0, target, 2);
-    g4 = ApproachByFactor(g4, Chi2_0, target, 1.2);
-    g4 = ApproachByFactor(g4, Chi2_0, target, 1.02);
-    g4 = ApproachByFactor(g4, Chi2_0, target, 1.0002);
+    g4 = ApproachByFactor(node, g4, Chi2_0, target, 2);
+    g4 = ApproachByFactor(node, g4, Chi2_0, target, 1.2);
+    g4 = ApproachByFactor(node, g4, Chi2_0, target, 1.02);
+    g4 = ApproachByFactor(node, g4, Chi2_0, target, 1.0002);
 
     return g4;
 }
@@ -131,7 +254,8 @@ Double_t TRestSensitivity::GetCoupling(Double_t sigma, Double_t precision) {
 ///////////////////////////////////////////////
 /// \brief It returns the Log(L) for the experiment and coupling given by argument.
 ///
-Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experiment, Double_t g4) {
+Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experiment, Double_t node,
+                                                 Double_t g4) {
     Double_t lhood = 0;
     if (!experiment->IsDataReady()) {
         RESTError << "TRestSensitivity::UnbinnedLogLikelihood. Experiment " << experiment->GetName()
@@ -139,10 +263,41 @@ 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 {
+        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
+    /// conditions. For example, higher pressure inside the detector gains in signal sensitivity but
+    /// it will produce also higher background.
+
+    if (experiment->GetBackground()->HasNodes()) {
+        RESTWarning
+            << "TRestSensitivity::UnbinnedLogLikelihood is not ready to have background parameter nodes!"
+            << RESTendl;
+        return 0.0;
+    }
+
     Double_t signal = g4 * experiment->GetSignal()->GetTotalRate() * experiment->GetExposureInSeconds();
 
     lhood = -signal;
 
+    if (experiment->GetExperimentalCounts() == 0) return lhood;
+
     if (ROOT::IsImplicitMTEnabled()) ROOT::DisableImplicitMT();
 
     std::vector<std::vector<Double_t>> trackingData;
@@ -169,12 +324,12 @@ Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experime
 ///////////////////////////////////////////////
 /// \brief
 ///
-TH1D* TRestSensitivity::SignalStatisticalTest(Int_t N) {
+TH1D* TRestSensitivity::SignalStatisticalTest(Double_t node, Int_t N) {
     std::vector<Double_t> couplings;
     for (int n = 0; n < N; n++) {
         for (const auto& exp : fExperiments) exp->GetSignal()->RegenerateActiveNodeDensity();
 
-        Double_t coupling = TMath::Sqrt(TMath::Sqrt(GetCoupling()));
+        Double_t coupling = TMath::Sqrt(TMath::Sqrt(GetCoupling(node)));
         couplings.push_back(coupling);
     }
 
@@ -213,11 +368,254 @@ void TRestSensitivity::InitFromConfigFile() {
     Initialize();
 }
 
+/////////////////////////////////////////////
+/// \brief This method is used to obtain the list of curves that satisfy that each value inside
+/// the curve is placed at a specified level. E.g. if we provide a level 0.5, then the corresponding
+/// curve will be constructed with the central value extracted at each parameter point.
+//
+/// We may then construct the profile of the sensitivity curves at 98%, 95% and 68% C.L. as follows:
+///
+/// \code
+/// TRestSensitivity::GetLevelCurves( {0.01, 0.025, 0.16, 0.84, 0.975, 0.99} );
+/// \endcode
+///
+std::vector<std::vector<Double_t>> TRestSensitivity::GetLevelCurves(const std::vector<Double_t>& levels) {
+    std::vector<std::vector<Double_t>> curves(levels.size());
+
+    for (const auto& l : levels) {
+        if (l >= 1 || l <= 0) {
+            RESTError << "The level value should be between 0 and 1" << RESTendl;
+            return curves;
+        }
+    }
+
+    std::vector<int> intLevels;
+    for (const auto& l : levels) {
+        int val = (int)round(l * fCurves.size());
+        if (val >= (int)fCurves.size()) val = fCurves.size() - 1;
+        if (val < 0) val = 0;
+
+        intLevels.push_back(val);
+    }
+
+    for (size_t m = 0; m < fCurves[0].size(); m++) {
+        std::vector<Double_t> v;
+        for (size_t n = 0; n < fCurves.size(); n++) v.push_back(fCurves[n][m]);
+
+        std::sort(v.begin(), v.begin() + v.size());
+
+        for (size_t n = 0; n < intLevels.size(); n++) curves[n].push_back(v[intLevels[n]]);
+    }
+
+    return curves;
+}
+
+TCanvas* TRestSensitivity::DrawCurves() {
+    if (fCanvas != NULL) {
+        delete fCanvas;
+        fCanvas = NULL;
+    }
+    fCanvas = new TCanvas("canv", "This is the canvas title", 600, 450);
+    fCanvas->Draw();
+
+    TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97);
+    // pad1->Divide(2, 2);
+    pad1->Draw();
+
+    ////// Drawing reflectivity versus angle
+    //   pad1->cd()->SetLogx();
+    pad1->cd()->SetRightMargin(0.09);
+    pad1->cd()->SetLeftMargin(0.15);
+    pad1->cd()->SetBottomMargin(0.15);
+
+    std::vector<TGraph*> graphs;
+
+    for (size_t n = 0; n < 20; n++) {
+        std::string grname = "gr" + IntegerToString(n);
+        TGraph* gr = new TGraph();
+        gr->SetName(grname.c_str());
+        for (size_t m = 0; m < this->GetCurve(n).size(); m++)
+            gr->SetPoint(gr->GetN(), fParameterizationNodes[m],
+                         TMath::Sqrt(TMath::Sqrt(this->GetCurve(n)[m])));
+
+        gr->SetLineColorAlpha(kBlue + n, 0.3);
+        gr->SetLineWidth(1);
+        graphs.push_back(gr);
+    }
+
+    TGraph* avGr = new TGraph();
+    std::vector<Double_t> avCurve = GetAveragedCurve();
+    for (size_t m = 0; m < avCurve.size(); m++)
+        avGr->SetPoint(avGr->GetN(), fParameterizationNodes[m], TMath::Sqrt(TMath::Sqrt(avCurve[m])));
+    avGr->SetLineColor(kBlack);
+    avGr->SetLineWidth(2);
+
+    graphs[0]->GetXaxis()->SetLimits(0, 0.25);
+    //   graphs[0]->GetHistogram()->SetMaximum(1);
+    //   graphs[0]->GetHistogram()->SetMinimum(lowRange);
+
+    graphs[0]->GetXaxis()->SetTitle("mass [eV]");
+    graphs[0]->GetXaxis()->SetTitleSize(0.05);
+    graphs[0]->GetXaxis()->SetLabelSize(0.05);
+    graphs[0]->GetYaxis()->SetTitle("g_{a#gamma} [10^{-10} GeV^{-1}]");
+    graphs[0]->GetYaxis()->SetTitleOffset(1.5);
+    graphs[0]->GetYaxis()->SetTitleSize(0.05);
+    // graphs[0]->GetYaxis()->SetLabelSize(0.05);
+    // graphs[0]->GetYaxis()->SetLabelOffset(0.0);
+    // pad1->cd()->SetLogy();
+    graphs[0]->Draw("AL");
+    for (unsigned int n = 1; n < graphs.size(); n++) graphs[n]->Draw("L");
+    avGr->Draw("L");
+
+    /*
+Double_t lx1 = 0.6, ly1 = 0.75, lx2 = 0.9, ly2 = 0.95;
+if (eLegendCoords.size() > 0) {
+    lx1 = eLegendCoords[0];
+    ly1 = eLegendCoords[1];
+    lx2 = eLegendCoords[2];
+    ly2 = eLegendCoords[3];
+}
+TLegend* legend = new TLegend(lx1, ly1, lx2, ly2);
+
+legend->SetTextSize(0.03);
+legend->SetHeader("Energies", "C");  // option "C" allows to center the header
+for (unsigned int n = 0; n < energies.size(); n++) {
+    std::string lname = "gr" + IntegerToString(n);
+    std::string ltitle = DoubleToString(energies[n]) + " keV";
+
+    legend->AddEntry(lname.c_str(), ltitle.c_str(), "l");
+}
+legend->Draw();
+    */
+
+    return fCanvas;
+}
+
+TCanvas* TRestSensitivity::DrawLevelCurves() {
+    if (fCanvas != NULL) {
+        delete fCanvas;
+        fCanvas = NULL;
+    }
+    fCanvas = new TCanvas("canv", "This is the canvas title", 500, 400);
+    fCanvas->Draw();
+    fCanvas->SetLeftMargin(0.15);
+    fCanvas->SetRightMargin(0.04);
+    fCanvas->SetLogx();
+
+    std::vector<std::vector<Double_t>> levelCurves = GetLevelCurves({0.025, 0.16, 0.375, 0.625, 0.84, 0.975});
+
+    std::vector<TGraph*> graphs;
+    for (size_t n = 0; n < levelCurves.size(); n++) {
+        std::string grname = "gr" + IntegerToString(n);
+        TGraph* gr = new TGraph();
+        gr->SetName(grname.c_str());
+        for (size_t m = 0; m < levelCurves[n].size(); m++)
+            gr->SetPoint(gr->GetN(), fParameterizationNodes[m], TMath::Sqrt(TMath::Sqrt(levelCurves[n][m])));
+
+        gr->SetLineColor(kBlack);
+        gr->SetLineWidth(1);
+        graphs.push_back(gr);
+    }
+
+    TGraph* centralGr = new TGraph();
+    std::vector<Double_t> centralCurve = GetLevelCurves({0.5})[0];
+    for (size_t m = 0; m < centralCurve.size(); m++)
+        centralGr->SetPoint(centralGr->GetN(), fParameterizationNodes[m],
+                            TMath::Sqrt(TMath::Sqrt(centralCurve[m])));
+    centralGr->SetLineColor(kBlack);
+    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]->GetXaxis()->SetTitle("mass [eV]");
+    graphs[0]->GetXaxis()->SetTitleSize(0.04);
+    graphs[0]->GetXaxis()->SetTitleOffset(1.15);
+    graphs[0]->GetXaxis()->SetLabelSize(0.04);
+
+    //   graphs[0]->GetYaxis()->SetLabelFont(43);
+    graphs[0]->GetYaxis()->SetTitle("g_{a#gamma} [10^{-10} GeV^{-1}]");
+    graphs[0]->GetYaxis()->SetTitleOffset(1.5);
+    graphs[0]->GetYaxis()->SetTitleSize(0.04);
+    graphs[0]->GetYaxis()->SetLabelSize(0.04);
+    // graphs[0]->GetYaxis()->SetLabelOffset(0);
+    // graphs[0]->GetYaxis()->SetLabelFont(43);
+    graphs[0]->Draw("AL");
+
+    TGraph* randomGr = new TGraph();
+    std::vector<Double_t> randomCurve = fCurves[13];
+    for (size_t m = 0; m < randomCurve.size(); m++)
+        randomGr->SetPoint(randomGr->GetN(), fParameterizationNodes[m],
+                           TMath::Sqrt(TMath::Sqrt(randomCurve[m])));
+    randomGr->SetLineColor(kBlack);
+    randomGr->SetLineWidth(1);
+    randomGr->SetMarkerSize(0.3);
+    randomGr->SetMarkerStyle(4);
+
+    std::vector<TGraph*> shadeGraphs;
+
+    int M = (int)levelCurves.size();
+    for (int x = 0; x < M / 2; x++) {
+        TGraph* shade = new TGraph();
+        int N = levelCurves[0].size();
+        for (size_t m = 0; m < levelCurves[0].size(); m++)
+            shade->SetPoint(shade->GetN(), fParameterizationNodes[m],
+                            TMath::Sqrt(TMath::Sqrt(levelCurves[x][m])));
+        for (int m = N - 1; m >= 0; --m)
+            shade->SetPoint(shade->GetN(), fParameterizationNodes[m],
+                            TMath::Sqrt(TMath::Sqrt(levelCurves[M - 1 - x][m])));
+        shade->SetFillColorAlpha(kRed, 0.25);
+        shade->Draw("f");
+        shadeGraphs.push_back(shade);
+    }
+
+    for (unsigned int n = 1; n < graphs.size(); n++) graphs[n]->Draw("Lsame");
+    randomGr->Draw("LPsame");
+    // centralGr->Draw("Lsame");
+
+    return fCanvas;
+}
+
+/////////////////////////////////////////////
+/// \brief It scans all the experiment signals parametric nodes to build a complete list
+/// of nodes used to build a complete sensitivity curve. Some experiments may be
+/// sensitivy to a particular node, while others may be sensitivy to another. If more
+/// than one experiment is sensitivy to a given node, the sensitivity will be combined
+/// later on.
+///
+void TRestSensitivity::ExtractExperimentParameterizationNodes(Bool_t rescan) {
+    if (fParameterizationNodes.empty() || rescan) {
+        fParameterizationNodes.clear();
+
+        for (const auto& experiment : fExperiments) {
+            std::vector<Double_t> nodes = experiment->GetSignal()->GetParameterizationNodes();
+            fParameterizationNodes.insert(fParameterizationNodes.end(), nodes.begin(), nodes.end());
+
+            std::sort(fParameterizationNodes.begin(), fParameterizationNodes.end());
+            auto last = std::unique(fParameterizationNodes.begin(), fParameterizationNodes.end());
+            fParameterizationNodes.erase(last, fParameterizationNodes.end());
+        }
+    }
+}
+
+void TRestSensitivity::PrintParameterizationNodes() {
+    std::cout << "Curve sensitivity nodes: ";
+    for (const auto& node : fParameterizationNodes) std::cout << node << "\t";
+    std::cout << std::endl;
+}
+
 /////////////////////////////////////////////
 /// \brief Prints on screen the information about the metadata members of TRestAxionSolarFlux
 ///
 void TRestSensitivity::PrintMetadata() {
     TRestMetadata::PrintMetadata();
 
+    RESTMetadata << " - Number of parameterization nodes : " << GetNumberOfNodes() << RESTendl;
+    RESTMetadata << " - Number of experiments loaded : " << GetNumberOfExperiments() << RESTendl;
+    RESTMetadata << " - Number of sensitivity curves generated : " << GetNumberOfCurves() << RESTendl;
+    RESTMetadata << " " << RESTendl;
+    RESTMetadata << " You may access experiment info using TRestSensitivity::GetExperiment(n)" << RESTendl;
+
     RESTMetadata << "----" << RESTendl;
 }
diff --git a/source/libraries/axion b/source/libraries/axion
index 3800d90f4..3802d108b 160000
--- a/source/libraries/axion
+++ b/source/libraries/axion
@@ -1 +1 @@
-Subproject commit 3800d90f41622ec0c465f7b4d4c5f384e4d9433a
+Subproject commit 3802d108bed118f4777ca142273aabaf368770c0