Merge pull request #508 from rest-for-physics/jgalan_sensitivity

Sensitivity classes updated, axionlib upgrade, new version header will be generated 2.4.3

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);

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)

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

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
@@ -483,7 +518,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);

source/framework/sensitivity/inc/TRestComponent.h

@@ -82,9 +82,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 +101,9 @@ class TRestComponent : public TRestMetadata {
     void Initialize() override;
     void RegenerateHistograms(UInt_t seed = 0);
 
+    /// It returns true if any nodes have been defined.
+    Bool_t HasNodes() { return !fParameterizationNodes.empty(); }
+
     virtual void RegenerateActiveNodeDensity() {}
 
     std::string GetNature() const { return fNature; }
@@ -113,6 +113,7 @@ class TRestComponent : public TRestMetadata {
     Int_t GetSamples() { return fSamples; }
     Int_t GetActiveNode() { return fActiveNode; }
     Double_t GetActiveNodeValue() { return fParameterizationNodes[fActiveNode]; }
+    std::vector<Double_t> GetParameterizationNodes() { return fParameterizationNodes; }
 
     std::vector<std::string> GetVariables() const { return fVariables; }
     std::vector<TVector2> GetRanges() const { return fRanges; }
@@ -127,6 +128,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;

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,9 @@ class TRestExperiment : public TRestMetadata {
     /// Only if it is true we will be able to calculate the LogLikelihood
     Bool_t fDataReady = false;  //<
 
+    /// It keeps track on the number of counts inside the dataset
+    Int_t fExperimentalCounts = 0;  //<
+
     /// Internal process random generator
     TRandom3* fRandom = nullptr;  //!
 
@@ -64,6 +67,7 @@ class TRestExperiment : public TRestMetadata {
 
    public:
     void GenerateMockDataSet();
+    Int_t GetExperimentalCounts() const { return fExperimentalCounts; }
 
     Bool_t IsMockData() const { return fMockData; }
     Bool_t IsDataReady() const { return fDataReady; }

source/framework/sensitivity/inc/TRestExperimentList.h

@@ -50,9 +50,15 @@ 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;  //<
+
     /// 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";
+
     /// If not null this will be the common signal used in each experiment
     TRestComponent* fSignal = nullptr;  //<
 
@@ -71,6 +77,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())

source/framework/sensitivity/inc/TRestSensitivity.h

@@ -29,22 +29,51 @@
 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 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 +84,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

source/framework/sensitivity/src/TRestComponent.cxx

@@ -627,6 +627,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.

source/framework/sensitivity/src/TRestExperiment.cxx

@@ -102,12 +102,15 @@ void TRestExperiment::GenerateMockDataSet() {
     Double_t meanCounts = GetBackground()->GetTotalRate() * fExposureTime * units("s");
 
     Int_t N = fRandom->Poisson(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 +121,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;
@@ -283,7 +287,7 @@ void TRestExperiment::PrintMetadata() {
         }
     }
 
-    RESTMetadata << " - Experimental counts : " << *fExperimentalData.GetDataFrame().Count() << RESTendl;
+    RESTMetadata << " - Experimental counts : " << fExperimentalCounts << RESTendl;
 
     RESTMetadata << "----" << RESTendl;
 }