Merge pull request #63 from rest-for-physics/solar_flux_update

TRestAxionSolarFlux now is an abstract class that allows to encapsulate different solar flux types

Author: jgalan

.github/pr-badge.yml

@@ -13,7 +13,5 @@
   message: "Ok: $additions"
   color: "green"
   when: "$additions < 100"
-- imageUrl: "https://gitlab.cern.ch/rest-for-physics/axionlib/badges/$branchName/pipeline.svg"
-  url: "https://gitlab.cern.ch/rest-for-physics/axionlib/-/commits/$branchName"
-- imageUrl: "https://github.com/rest-for-physics/axionlib/actions/workflows/validation.yml/badge.svg?branch=$branchName"
+- imageUrl: "https://github.com/rest-for-physics/axionlib/actions/workflows/frameworkValidation.yml/badge.svg?branch=$branchName"
   url: "https://github.com/rest-for-physics/axionlib/commits/$branchName"

.github/workflows/validation.yml

@@ -187,7 +187,7 @@ jobs:
         wget https://rest-for-physics.github.io/axionlib-data/solarFlux/Primakoff_Gianotti_201904.dat
         wget https://rest-for-physics.github.io/axionlib-data/solarFlux/Primakoff_LennertHoof_202203.dat
         python solarTests.py
-        python solarPlot.py
+        python solarPlotQCD.py
         python compare.py
 
   Ray-tracing:

CMakeLists.txt

@@ -4,7 +4,7 @@ cmake_minimum_required(VERSION 3.5 FATAL_ERROR)
 include(ExternalProject)
 set(CMAKE_MACOSX_RPATH 1)
 
-set(LibraryVersion "1.1")
+set(LibraryVersion "1.2")
 add_definitions(-DLIBRARY_VERSION="${LibraryVersion}")
 
 if (${REST_MPFR} MATCHES "ON")

data

@@ -47,7 +47,7 @@ class TRestAxionEventProcess : public TRestEventProcess {
     /// The rotation angle around CenterSetup with respect to X-axis
     Double_t fExternalPitch = 0;  //<
 
-    /// If enabled it will skip the end rotation that recovers the original axion trajectory direction
+    /// If enabled it skips the end rotation that recovers the original axion trajectory direction
     Bool_t fLocalAxis = true;  //<
 
    protected:

inc/TRestAxionEventProcess.h

@@ -27,18 +27,11 @@
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TRandom3.h>
-#include <TRestAxionSolarModel.h>
 #include <TRestMetadata.h>
 
 //! A metadata class to load tabulated solar axion fluxes
 class TRestAxionSolarFlux : public TRestMetadata {
    private:
-    /// The filename containning the solar flux table with continuum spectrum
-    std::string fFluxDataFile = "";  //<
-
-    /// The filename containning the solar flux spectra for monochromatic spectrum
-    std::string fFluxSptFile = "";  //<
-
     /// Axion coupling. Defines coupling type and strength.
     std::string fCouplingType = "";  //<
 
@@ -48,101 +41,55 @@ class TRestAxionSolarFlux : public TRestMetadata {
     /// Seed used in random generator
     Int_t fSeed = 0;  //<
 
-    /// It will be used when loading `.flux` files to define the input file energy binsize in eV.
-    Double_t fBinSize = 0;  //<
-
-    /// It will be used when loading `.flux` files to define the threshold for peak identification
-    Double_t fPeakSigma = 0;  //<
-
-    /// The tabulated solar flux continuum spectra TH1F(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius
-    std::vector<TH1F*> fFluxTable;  //!
-
-    /// The tabulated solar flux in cm-2 s-1 for a number of monochromatic energies versus solar radius
-    std::map<Double_t, TH1F*> fFluxLines;  //!
-
-    /// Accumulative integrated solar flux for each solar ring for continuum spectrum (renormalized to unity)
-    std::vector<Double_t> fFluxTableIntegrals;  //!
-
-    /// Accumulative integrated solar flux for each monochromatic energy (renormalized to unity)
-    std::vector<Double_t> fFluxLineIntegrals;  //!
-
-    /// Total solar flux for monochromatic contributions
-    Double_t fTotalMonochromaticFlux = 0;  //!
-
-    /// Total solar flux for monochromatic contributions
-    Double_t fTotalContinuumFlux = 0;  //!
-
-    /// The ratio between monochromatic and total flux
-    Double_t fFluxRatio = 0;  //!
-
-    /// Random number generator
-    TRandom3* fRandom = nullptr;  //!
+    /// A metadata member to control if this class has been initialized
+    Bool_t fTablesLoaded = false;  //!
 
+   protected:
     /// A canvas pointer for drawing
     TCanvas* fCanvas = nullptr;  //!
 
-    /// A pointer to the continuum spectrum histogram
-    TH1F* fContinuumHist = nullptr;  //!
-
-    /// A pointer to the monochromatic spectrum histogram
-    TH1F* fMonoHist = nullptr;  //!
-
-    /// A pointer to the superposed monochromatic and continuum spectrum histogram
-    TH1F* fTotalHist = nullptr;  //!
-
-    /// A metadata member to control if the tables have been loaded
-    Bool_t fTablesLoaded = false;  //!
+    /// Random number generator
+    TRandom3* fRandom = nullptr;  //!
 
-    void Initialize();
+    TRestAxionSolarFlux();
+    TRestAxionSolarFlux(const char* cfgFileName, std::string name = "");
 
-    void ReadFluxFile();
-    void LoadContinuumFluxTable();
-    void LoadMonoChromaticFluxTable();
-    void IntegrateSolarFluxes();
+    /// It defines how to read the solar tables at the inhereted class
+    virtual Bool_t LoadTables() = 0;
 
    public:
-    /// It returns true if continuum flux spectra was loaded
-    Bool_t isSolarTableLoaded() { return fFluxTable.size() > 0; }
-
-    /// It returns true if monochromatic flux spectra was loaded
-    Bool_t isSolarSpectrumLoaded() { return fFluxLines.size() > 0; }
+    /// It is required in order to load solar flux tables into memory
+    void Initialize();
 
     /// It returns the integrated flux at earth in cm-2 s-1 for the given energy range
-    Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1));
+    virtual Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1), Double_t mass = 0) = 0;
 
     /// It returns the total integrated flux at earth in cm-2 s-1
-    Double_t GetTotalFlux() { return fTotalContinuumFlux + fTotalMonochromaticFlux; }
+    virtual Double_t GetTotalFlux(Double_t mass = 0) = 0;
 
-    std::pair<Double_t, Double_t> GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1));
+    /// It defines how to generate Monte Carlo energy and radius values to reproduce the flux
+    virtual std::pair<Double_t, Double_t> GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1),
+                                                                   Double_t mass = 0) = 0;
 
-    void LoadTables();
+    /// It returns an energy integrated spectrum in cm-2 s-1 keV-1
+    virtual TH1F* GetEnergySpectrum(Double_t m = 0) = 0;
 
-    TH1F* GetContinuumSpectrum();
-    TH1F* GetMonochromaticSpectrum();
-    TH1F* GetTotalSpectrum();
+    virtual TCanvas* DrawSolarFlux();
 
-    TH1F* GetFluxHistogram(std::string fname, Double_t binSize = 0.01);
-    TCanvas* DrawFluxFile(std::string fname, Double_t binSize = 0.01);
-    TCanvas* DrawSolarFlux();
-
-    /// Tables might be loaded using a solar model description by TRestAxionSolarModel
-    void InitializeSolarTable(TRestAxionSolarModel* model) {
-        // TOBE implemented
-        // This method should initialize the tables fFluxTable and fFluxLines
+    virtual void ExportTables(Bool_t ascii = false) {
+        RESTWarning << "TRestAxionSolarFlux::ExportTables must be re-implemented in the inherited class"
+                    << RESTendl;
     }
 
-    void ExportTables(Bool_t ascii = false);
+    Bool_t AreTablesLoaded() { return fTablesLoaded; }
 
-    void PrintMetadata();
+    TH1F* GetFluxHistogram(std::string fname, Double_t binSize = 0.01);
+    TCanvas* DrawFluxFile(std::string fname, Double_t binSize = 0.01);
 
-    void PrintContinuumSolarTable();
-    void PrintIntegratedRingFlux();
-    void PrintMonoChromaticFlux();
+    void PrintMetadata();
 
-    TRestAxionSolarFlux();
-    TRestAxionSolarFlux(const char* cfgFileName, std::string name = "");
     ~TRestAxionSolarFlux();
 
-    ClassDef(TRestAxionSolarFlux, 1);
+    ClassDef(TRestAxionSolarFlux, 2);
 };
 #endif

inc/TRestAxionSolarFlux.h

@@ -0,0 +1,130 @@
+/*************************************************************************
+ * This file is part of the REST software framework.                     *
+ *                                                                       *
+ * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
+ * For more information see http://gifna.unizar.es/trex                  *
+ *                                                                       *
+ * REST is free software: you can redistribute it and/or modify          *
+ * it under the terms of the GNU General Public License as published by  *
+ * the Free Software Foundation, either version 3 of the License, or     *
+ * (at your option) any later version.                                   *
+ *                                                                       *
+ * REST is distributed in the hope that it will be useful,               *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
+ * GNU General Public License for more details.                          *
+ *                                                                       *
+ * You should have a copy of the GNU General Public License along with   *
+ * REST in $REST_PATH/LICENSE.                                           *
+ * If not, see http://www.gnu.org/licenses/.                             *
+ * For the list of contributors see $REST_PATH/CREDITS.                  *
+ *************************************************************************/
+
+#ifndef _TRestAxionSolarQCDFlux
+#define _TRestAxionSolarQCDFlux
+
+#include <TRestAxionSolarFlux.h>
+#include <TRestAxionSolarModel.h>
+
+//! A metadata class to load tabulated solar axion fluxes. Mass independent.
+class TRestAxionSolarQCDFlux : public TRestAxionSolarFlux {
+   private:
+    /// The filename containning the solar flux table with continuum spectrum
+    std::string fFluxDataFile = "";  //<
+
+    /// The filename containning the solar flux spectra for monochromatic spectrum
+    std::string fFluxSptFile = "";  //<
+
+    /// It will be used when loading `.flux` files to define the input file energy binsize in eV.
+    Double_t fBinSize = 0;  //<
+
+    /// It will be used when loading `.flux` files to define the threshold for peak identification
+    Double_t fPeakSigma = 0;  //<
+
+    /// The tabulated solar flux continuum spectra TH1F(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius
+    std::vector<TH1F*> fFluxTable;  //!
+
+    /// The tabulated solar flux in cm-2 s-1 for a number of monochromatic energies versus solar radius
+    std::map<Double_t, TH1F*> fFluxLines;  //!
+
+    /// Accumulative integrated solar flux for each solar ring for continuum spectrum (renormalized to unity)
+    std::vector<Double_t> fFluxTableIntegrals;  //!
+
+    /// Accumulative integrated solar flux for each monochromatic energy (renormalized to unity)
+    std::vector<Double_t> fFluxLineIntegrals;  //!
+
+    /// Total solar flux for monochromatic contributions
+    Double_t fTotalMonochromaticFlux = 0;  //!
+
+    /// Total solar flux for monochromatic contributions
+    Double_t fTotalContinuumFlux = 0;  //!
+
+    /// The ratio between monochromatic and total flux
+    Double_t fFluxRatio = 0;  //!
+
+    /// A pointer to the continuum spectrum histogram
+    TH1F* fContinuumHist = nullptr;  //!
+
+    /// A pointer to the monochromatic spectrum histogram
+    TH1F* fMonoHist = nullptr;  //!
+
+    /// A pointer to the superposed monochromatic and continuum spectrum histogram
+    TH1F* fTotalHist = nullptr;  //!
+
+    void ReadFluxFile();
+    void LoadContinuumFluxTable();
+    void LoadMonoChromaticFluxTable();
+    void IntegrateSolarFluxes();
+
+   public:
+    /// It returns true if continuum flux spectra was loaded
+    Bool_t isSolarTableLoaded() { return fFluxTable.size() > 0; }
+
+    /// It returns true if monochromatic flux spectra was loaded
+    Bool_t isSolarSpectrumLoaded() { return fFluxLines.size() > 0; }
+
+    /// It returns the integrated flux at earth in cm-2 s-1 for the given energy range
+    Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1), Double_t mass = 0) override;
+
+    /// It defines how to generate Monte Carlo energy and radius values to reproduce the flux
+    std::pair<Double_t, Double_t> GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1),
+                                                           Double_t mass = 0) override;
+
+    /// It defines how to read the solar tables at the inhereted class
+    Bool_t LoadTables() override;
+
+    /// It returns the total integrated flux at earth in cm-2 s-1
+    Double_t GetTotalFlux(Double_t mass = 0) override {
+        return fTotalContinuumFlux + fTotalMonochromaticFlux;
+    }
+
+    /// It returns an energy integrated spectrum in cm-2 s-1 keV-1
+    TH1F* GetEnergySpectrum(Double_t m = 0) override { return GetTotalSpectrum(); }
+
+    TH1F* GetContinuumSpectrum();
+    TH1F* GetMonochromaticSpectrum();
+    TH1F* GetTotalSpectrum();
+
+    virtual TCanvas* DrawSolarFlux() override;
+
+    /// Tables might be loaded using a solar model description by TRestAxionSolarModel
+    void InitializeSolarTable(TRestAxionSolarModel* model) {
+        // TOBE implemented
+        // This method should initialize the tables fFluxTable and fFluxLines
+    }
+
+    void ExportTables(Bool_t ascii = false) override;
+
+    void PrintMetadata() override;
+
+    void PrintContinuumSolarTable();
+    void PrintIntegratedRingFlux();
+    void PrintMonoChromaticFlux();
+
+    TRestAxionSolarQCDFlux();
+    TRestAxionSolarQCDFlux(const char* cfgFileName, std::string name = "");
+    ~TRestAxionSolarQCDFlux();
+
+    ClassDefOverride(TRestAxionSolarQCDFlux, 1);
+};
+#endif

inc/TRestAxionSolarQCDFlux.h

@@ -40,8 +40,8 @@
 pad1.Divide(2, 1)
 pad1.Draw()
 
-primakoffLH = ROOT.TRestAxionSolarFlux("fluxes.rml", "LennertHoofPrimakoff")
-primakoffG = ROOT.TRestAxionSolarFlux("fluxes.rml", "Gianotti")
+primakoffLH = ROOT.TRestAxionSolarQCDFlux("fluxes.rml", "LennertHoofPrimakoff")
+primakoffG = ROOT.TRestAxionSolarQCDFlux("fluxes.rml", "Gianotti")
 
 if primakoffLH.GetError():
     print(primakoffLH.GetErrorMessage())

pipeline/metadata/solarFlux/compare.py

@@ -63,7 +63,7 @@
 if args.samples != None:
     samples = args.samples
 
-validation = False
+validation = True
 if (
     rmlfile == "fluxes.rml"
     and fluxname == "combined"
@@ -83,8 +83,8 @@
 pad1.Divide(2, 2)
 pad1.Draw()
 
-combinedFlux = ROOT.TRestAxionSolarFlux(rmlfile, fluxname)
-combinedFlux.LoadTables()
+combinedFlux = ROOT.TRestAxionSolarQCDFlux(rmlfile, fluxname)
+combinedFlux.Initialize()
 combinedFlux.PrintMetadata()
 
 if combinedFlux.GetError():
@@ -166,14 +166,14 @@
 print("\nMaximum energy bin is " + str(enSpt.GetMaximumBin()))
 if validation:
     if enSpt.GetMaximumBin() != 8001:
-        print("\nMaximum Bin is not the expected one! Exit code : 1")
+        print("\nMaximum Bin is not the expected one (8001)! Exit code : 1")
         exit(1)
 
 print("\nMaximum radius bin is " + str(rSpt.GetMaximumBin()))
 
 if validation:
     if rSpt.GetMaximumBin() != 25:
-        print("\nMaximum Bin is not the expected one! Exit code : 2")
+        print("\nMaximum Bin is not the expected one (25)! Exit code : 2")
         exit(2)
 
 exit(0)

pipeline/metadata/solarFlux/solarPlot.py

@@ -29,8 +29,8 @@
 ROOT.gSystem.Load("libRestFramework.so")
 ROOT.gSystem.Load("libRestAxion.so")
 
-monoFlux = ROOT.TRestAxionSolarFlux("fluxes.rml", "mono")
-monoFlux.LoadTables()
+monoFlux = ROOT.TRestAxionSolarQCDFlux("fluxes.rml", "mono")
+monoFlux.Initialize()
 monoFlux.PrintMetadata()
 
 if monoFlux.GetError():
@@ -47,7 +47,7 @@
 
 print("[\033[92m OK \x1b[0m]")
 
-continuumFlux = ROOT.TRestAxionSolarFlux("fluxes.rml", "Gianotti")
+continuumFlux = ROOT.TRestAxionSolarQCDFlux("fluxes.rml", "Gianotti")
 continuumFlux.LoadTables()
 continuumFlux.PrintMetadata()