Sensitivity and masks classes updates

source/framework/masks/inc/TRestSpiderMask.h

@@ -28,7 +28,6 @@
 /// A class used to define and generate a spider structure mask
 class TRestSpiderMask : public TRestPatternMask {
    private:
-    void Initialize() override;
     /// The angle between two consecutive spider arms measured in radians.
     Double_t fArmsSeparationAngle = 0;  //<
 
@@ -48,6 +47,8 @@ class TRestSpiderMask : public TRestPatternMask {
     std::vector<std::pair<Double_t, Double_t>> fNegativeRanges;  //!
 
    public:
+    void Initialize() override;
+
     void GenerateSpider();
 
     virtual Int_t GetRegion(Double_t& x, Double_t& y) override;
@@ -58,6 +59,12 @@ class TRestSpiderMask : public TRestPatternMask {
     /// It returns the angular width of each spider arm in radians
     Double_t GetArmsWidth() { return fArmsWidth; }
 
+    /// It returns the number of arms in the spider structure
+    size_t GetNumberOfArms() {
+        if (fArmsSeparationAngle == 0) return 0;
+        return (size_t)(2 * TMath::Pi() / fArmsSeparationAngle);
+    }
+
     /// It returns the inner ring radius that defines the inner start of the spider structure
     Double_t GetInitialRadius() { return fInitialRadius; }
 

source/framework/masks/src/TRestPatternMask.cxx

@@ -168,6 +168,7 @@ TCanvas* TRestPatternMask::DrawMonteCarlo(Int_t nSamples) {
         fCanvas = NULL;
     }
     fCanvas = new TCanvas("canv", "This is the canvas title", 1200, 1200);
+    fCanvas->SetRealAspectRatio();
     fCanvas->Draw();
 
     TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97);
@@ -213,12 +214,12 @@ TCanvas* TRestPatternMask::DrawMonteCarlo(Int_t nSamples) {
         if (nGraphs == 0) {
             gr->SetLineColor(kBlack);
             gr->SetMarkerColor(kBlack);
-            gr->SetMarkerSize(0.6);
+            gr->SetMarkerSize(0.3);
             gr->SetLineWidth(2);
         } else {
             gr->SetMarkerColor((nGraphs * 3) % 29 + 20);
             gr->SetLineColor((nGraphs * 3) % 29 + 20);
-            gr->SetMarkerSize(0.6);
+            gr->SetMarkerSize(0.3);
             gr->SetLineWidth(2);
         }
 

source/framework/masks/src/TRestSpiderMask.cxx

@@ -139,11 +139,11 @@ TRestSpiderMask::TRestSpiderMask() : TRestPatternMask() { Initialize(); }
 /// corresponding TRestSpiderMask section inside the RML.
 ///
 TRestSpiderMask::TRestSpiderMask(const char* cfgFileName, std::string name) : TRestPatternMask(cfgFileName) {
-    Initialize();
-
     LoadConfigFromFile(fConfigFileName, name);
 
     if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata();
+
+    Initialize();
 }
 
 ///////////////////////////////////////////////

source/framework/sensitivity/inc/TRestComponent.h

@@ -64,6 +64,9 @@ class TRestComponent : public TRestMetadata {
     /// It defines the increasing step for automatic parameter list generation
     Double_t fStepParameterValue = 0;  //<
 
+    /// It keeps track of total rates for quick access
+    std::vector<Double_t> fCachedRates;  //<
+
     /// It true the parametric values automatically generated will grow exponentially
     Bool_t fExponential = false;  //<
 
@@ -117,6 +120,7 @@ class TRestComponent : public TRestMetadata {
 
     /// It returns true if any nodes have been defined.
     Bool_t HasNodes() { return !fParameterizationNodes.empty(); }
+    size_t GetNumberOfParameterizationNodes() { return fParameterizationNodes.size(); }
 
     virtual void RegenerateActiveNodeDensity() {}
 
@@ -151,7 +155,7 @@ class TRestComponent : public TRestMetadata {
     Int_t FindActiveNode(Double_t node);
     Int_t SetActiveNode(Double_t node);
     Int_t SetActiveNode(Int_t n) {
-        fActiveNode = n;
+        fActiveNode = n % fParameterizationNodes.size();
         return fActiveNode;
     }
 
@@ -191,6 +195,6 @@ class TRestComponent : public TRestMetadata {
     TRestComponent();
     ~TRestComponent();
 
-    ClassDefOverride(TRestComponent, 6);
+    ClassDefOverride(TRestComponent, 7);
 };
 #endif

source/framework/sensitivity/inc/TRestResponse.h

@@ -64,17 +64,31 @@ class TRestResponse : public TRestMetadata {
     std::string GetVariable() const { return fVariable; }
 
     TVector2 GetInputRange() const {
+        if (fResponseMatrix.empty()) {
+            RESTError
+                << " TRestResponse::GetInputRange. Response matrix not loaded yet, try using LoadResponse"
+                << RESTendl;
+            return 0;
+        }
         return TVector2(fOrigin.X(), fOrigin.X() + fResponseMatrix[0].size() * fBinSize);
     }
 
     TVector2 GetOutputRange() const {
+        if (fResponseMatrix.empty()) {
+            RESTError
+                << " TRestResponse::GetOutputRange. Response matrix not loaded yet, try using LoadResponse"
+                << RESTendl;
+            return 0;
+        }
         return TVector2(fOrigin.Y(), fOrigin.Y() + fResponseMatrix.size() * fBinSize);
     }
 
     void Initialize() override;
 
     void LoadResponse(Bool_t transpose = true);
 
+    Double_t GetOverallEfficiency(Double_t from, Double_t to);
+
     std::vector<std::pair<Double_t, Double_t>> GetResponse(Double_t input);
 
     void PrintResponseMatrix(Int_t fromRow, Int_t toRow);

source/framework/sensitivity/inc/TRestSensitivity.h

@@ -29,14 +29,17 @@
 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;  //<
 
+    /// If disabled the experiments will not be saved to disk
+    Bool_t fSaveExperiments = false;  //<
+
     /// 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
 
@@ -66,6 +69,7 @@ class TRestSensitivity : public TRestMetadata {
     void GenerateCurves(Int_t N);
 
     std::vector<Double_t> GetCurve(size_t n = 0);
+    std::vector<Double_t> ExportCurve(size_t n = 0) { return GetCurve(n); }
     std::vector<Double_t> GetAveragedCurve();
     std::vector<std::vector<Double_t>> GetLevelCurves(const std::vector<Double_t>& levels);
 
@@ -74,7 +78,8 @@ class TRestSensitivity : public TRestMetadata {
 
     TH1D* SignalStatisticalTest(Double_t node, Int_t N);
 
-    void Freeze() { fFrozen = true; }
+    void Freeze(Bool_t freeze = true) { fFrozen = freeze; }
+    void SaveExperiments(Bool_t save = true) { fSaveExperiments = save; }
 
     std::vector<TRestExperiment*> GetExperiments() { return fExperiments; }
     TRestExperiment* GetExperiment(const size_t& n) {
@@ -90,13 +95,15 @@ class TRestSensitivity : public TRestMetadata {
 
     void PrintMetadata() override;
 
+    virtual Int_t Write(const char* name = nullptr, Int_t option = 0, Int_t bufsize = 0) override;
+
     TCanvas* DrawCurves();
     TCanvas* DrawLevelCurves();
 
     TRestSensitivity(const char* cfgFileName, const std::string& name = "");
     TRestSensitivity();
     ~TRestSensitivity();
 
-    ClassDefOverride(TRestSensitivity, 2);
+    ClassDefOverride(TRestSensitivity, 4);
 };
 #endif

source/framework/sensitivity/src/TRestComponent.cxx

@@ -62,7 +62,7 @@ TRestComponent::TRestComponent() {}
 /// The default behaviour is that the config file must be specified with
 /// full path, absolute or relative.
 ///
-/// \param cfgFileName A const char* giving the path to an RML file.
+// \param cfgFileName A const char* giving the path to an RML file.
 /// \param name The name of the specific metadata. It will be used to find the
 /// corresponding TRestAxionMagneticField section inside the RML.
 ///
@@ -103,6 +103,9 @@ void TRestComponent::Initialize() {
     } else {
         if (!fParameterizationNodes.empty()) FillHistograms();
     }
+
+    fCachedRates.clear();
+    if (fNodeDensity.size() > 0) fCachedRates.resize(fNodeDensity.size(), 0.0);
 }
 
 /////////////////////////////////////////////
@@ -325,6 +328,10 @@ Double_t TRestComponent::GetRawRate(std::vector<Double_t> point) {
 /// The result will be returned in s-1.
 ///
 Double_t TRestComponent::GetTotalRate() {
+    if (fCachedRates.size() == 0 && fNodeDensity.size() > 0) fCachedRates.resize(fNodeDensity.size(), 0.0);
+    if (fActiveNode >= 0 && (Int_t)fCachedRates.size() > fActiveNode && fCachedRates[fActiveNode] > 0)
+        return fCachedRates[fActiveNode];
+
     THnD* dHist = GetDensityForActiveNode();
     if (!dHist) return 0;
 
@@ -343,6 +350,8 @@ Double_t TRestComponent::GetTotalRate() {
         if (!skip) integral += GetRate(point);
     }
 
+    if ((Int_t)fCachedRates.size() > fActiveNode) fCachedRates[fActiveNode] = integral;
+
     return integral;
 }
 

source/framework/sensitivity/src/TRestComponentDataSet.cxx

@@ -126,11 +126,11 @@ TRestComponentDataSet::TRestComponentDataSet(const char* cfgFileName, const std:
 /// (or observables) that have been defined by the user.
 ///
 void TRestComponentDataSet::Initialize() {
-    TRestComponent::Initialize();
-
     SetSectionName(this->ClassName());
 
     LoadDataSets();
+
+    TRestComponent::Initialize();
 }
 
 /////////////////////////////////////////////

source/framework/sensitivity/src/TRestExperiment.cxx

@@ -173,7 +173,7 @@ void TRestExperiment::InitFromConfigFile() {
 
     auto ele = GetElement("addComponent");
     if (ele != nullptr) {
-        std::string filename = GetParameter("filename", ele, "");
+        std::string filename = SearchFile(GetParameter("filename", ele, ""));
         std::string component = GetParameter("component", ele, "");
 
         if (filename.empty())
@@ -221,7 +221,7 @@ void TRestExperiment::InitFromConfigFile() {
 
     if (!fSignal) {
         RESTError << "TRestExperiment : " << GetName() << RESTendl;
-        RESTError << "There was a problem initiazing the signal component!" << RESTendl;
+        RESTError << "There was a problem initializing the signal component!" << RESTendl;
         return;
     }
 

source/framework/sensitivity/src/TRestResponse.cxx

@@ -129,6 +129,28 @@ void TRestResponse::LoadResponse(Bool_t transpose) {
     RESTError << "Extension format - " << extension << " - not recognized!" << RESTendl;
 }
 
+///////////////////////////////////////////////
+/// \brief It is used to integrate the response matrix and obtain the overall efficiency
+/// inside an energy range given by argument.
+///
+Double_t TRestResponse::GetOverallEfficiency(Double_t from, Double_t to) {
+    if (fResponseMatrix.empty()) {
+        RESTError
+            << " TRestResponse::GetOverallEfficiency. Response matrix not loaded yet, try using LoadResponse"
+            << RESTendl;
+        return 0;
+    }
+
+    Int_t nBins = 0;
+    Double_t eff = 0;
+    for (double en = from; en < to; en += fBinSize) {
+        std::vector<std::pair<Double_t, Double_t>> effs = GetResponse(en);
+        for (const auto& ef : effs) eff += ef.second;
+        nBins++;
+    }
+    return eff / nBins;
+}
+
 /////////////////////////////////////////////
 /// \brief This method will return a vector of std::pair, each pair will contain the
 /// output energy together with the corresponding response (or efficiency), for the

source/framework/sensitivity/src/TRestSensitivity.cxx

@@ -76,7 +76,10 @@ TRestSensitivity::TRestSensitivity(const char* cfgFileName, const std::string& n
 /// \brief It will initialize the data frame with the filelist and column names
 /// (or observables) that have been defined by the user.
 ///
-void TRestSensitivity::Initialize() { SetSectionName(this->ClassName()); }
+void TRestSensitivity::Initialize() {
+    SetSectionName(this->ClassName());
+    ExtractExperimentParameterizationNodes();
+}
 
 ///////////////////////////////////////////////
 /// \brief It will return a value of the coupling, g4, such that (chi-chi0) gets
@@ -275,8 +278,6 @@ Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experime
     if (nd >= 0)
         experiment->GetSignal()->SetActiveNode(nd);
     else {
-        RESTWarning << "Node : " << node << " not found in signal : " << experiment->GetSignal()->GetName()
-                    << RESTendl;
         return 0.0;
     }
 
@@ -326,8 +327,15 @@ Double_t TRestSensitivity::UnbinnedLogLikelihood(const TRestExperiment* experime
 ///
 TH1D* TRestSensitivity::SignalStatisticalTest(Double_t node, Int_t N) {
     std::vector<Double_t> couplings;
+    Int_t nodeCheck = 0;
+    if (fExperiments.size() > 0) nodeCheck = fExperiments[0]->GetSignal()->GetActiveNode();
     for (int n = 0; n < N; n++) {
-        for (const auto& exp : fExperiments) exp->GetSignal()->RegenerateActiveNodeDensity();
+        for (const auto& exp : fExperiments) {
+            if (exp->GetSignal()->GetActiveNode() != nodeCheck)
+                RESTError << "TRestSensitivity::SignalStatisticalTest. Problem" << RESTendl;
+            exp->GetSignal()->SetActiveNode(exp->GetSignal()->GetActiveNode() + 1);
+            exp->GetSignal()->RegenerateActiveNodeDensity();
+        }
 
         Double_t coupling = TMath::Sqrt(TMath::Sqrt(GetCoupling(node)));
         couplings.push_back(coupling);
@@ -338,7 +346,7 @@ TH1D* TRestSensitivity::SignalStatisticalTest(Double_t node, Int_t N) {
     double max_value = *std::max_element(couplings.begin(), couplings.end());
 
     if (fSignalTest) delete fSignalTest;
-    fSignalTest = new TH1D("SignalTest", "A signal test", 100, 0.9 * min_value, 1.1 * max_value);
+    fSignalTest = new TH1D("SignalTest", "A signal test", 1000, 0.9 * min_value, 1.1 * max_value);
     for (const auto& coup : couplings) fSignalTest->Fill(coup);
 
     return fSignalTest;
@@ -496,11 +504,12 @@ TCanvas* TRestSensitivity::DrawLevelCurves() {
         delete fCanvas;
         fCanvas = NULL;
     }
-    fCanvas = new TCanvas("canv", "This is the canvas title", 500, 400);
+    fCanvas = new TCanvas("canv", "This is the canvas title", 600, 500);
     fCanvas->Draw();
     fCanvas->SetLeftMargin(0.15);
     fCanvas->SetRightMargin(0.04);
     fCanvas->SetLogx();
+    fCanvas->SetLogy();
 
     std::vector<std::vector<Double_t>> levelCurves = GetLevelCurves({0.025, 0.16, 0.375, 0.625, 0.84, 0.975});
 
@@ -526,9 +535,9 @@ TCanvas* TRestSensitivity::DrawLevelCurves() {
     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]->GetYaxis()->SetRangeUser(0, 0.31);
+    graphs[0]->GetXaxis()->SetRangeUser(0.001, 0.31);
+    graphs[0]->GetXaxis()->SetLimits(0.0001, 0.31);
     graphs[0]->GetXaxis()->SetTitle("mass [eV]");
     graphs[0]->GetXaxis()->SetTitleSize(0.04);
     graphs[0]->GetXaxis()->SetTitleOffset(1.15);
@@ -539,12 +548,13 @@ TCanvas* TRestSensitivity::DrawLevelCurves() {
     graphs[0]->GetYaxis()->SetTitleOffset(1.5);
     graphs[0]->GetYaxis()->SetTitleSize(0.04);
     graphs[0]->GetYaxis()->SetLabelSize(0.04);
+    graphs[0]->GetYaxis()->SetRangeUser(0.1, 30);
     // 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];
+    std::vector<Double_t> randomCurve = fCurves[GetNumberOfCurves() / 2];
     for (size_t m = 0; m < randomCurve.size(); m++)
         randomGr->SetPoint(randomGr->GetN(), fParameterizationNodes[m],
                            TMath::Sqrt(TMath::Sqrt(randomCurve[m])));
@@ -572,7 +582,7 @@ TCanvas* TRestSensitivity::DrawLevelCurves() {
 
     for (unsigned int n = 1; n < graphs.size(); n++) graphs[n]->Draw("Lsame");
     randomGr->Draw("LPsame");
-    // centralGr->Draw("Lsame");
+    //   centralGr->Draw("Lsame");
 
     return fCanvas;
 }
@@ -619,3 +629,13 @@ void TRestSensitivity::PrintMetadata() {
 
     RESTMetadata << "----" << RESTendl;
 }
+
+Int_t TRestSensitivity::Write(const char* name, Int_t option, Int_t bufsize) {
+    if (!fSaveExperiments) {
+        RESTInfo << "TRestSensitivity::Write. Removing experiments before writting to disk." << RESTendl;
+        RESTInfo << "Use TRestSensitivity::SaveExperiments( true ) to change this behaviour" << RESTendl;
+        fExperiments.clear();
+    }
+
+    return TRestMetadata::Write(name, option, bufsize);
+}