diff --git a/inc/TRestAxionSolarFlux.h b/inc/TRestAxionSolarFlux.h index 30c49340..9a917211 100644 --- a/inc/TRestAxionSolarFlux.h +++ b/inc/TRestAxionSolarFlux.h @@ -24,8 +24,8 @@ #define _TRestAxionSolarFlux #include -#include -#include +#include +#include #include #include @@ -58,23 +58,29 @@ class TRestAxionSolarFlux : public TRestMetadata { TRestAxionSolarFlux(const char* cfgFileName, std::string name = ""); /// It defines how to read the solar tables at the inhereted class - virtual Bool_t LoadTables(Double_t mass = 0) = 0; + virtual Bool_t LoadTables() = 0; public: /// It is required in order to load solar flux tables into memory void Initialize(); + /// It is required in order to load solar flux tables into memory for specific mass + void InitializeMass( Double_t mass ) { SetMass(mass); Initialize(); } + + /// Set mass and reinitialise + void SetMass( Double_t m ) { fMass = m; } //Initialize(); } + /// It returns the integrated flux at earth in cm-2 s-1 for the given energy range - virtual Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1), Double_t mass = 0) = 0; + virtual Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1)) = 0; /// It returns the total integrated flux at earth in cm-2 s-1 - virtual Double_t GetTotalFlux(Double_t mass = 0) = 0; + virtual Double_t GetTotalFlux() = 0; /// It defines how to generate Monte Carlo energy and radius values to reproduce the flux virtual std::pair GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1)) = 0; /// It returns an energy integrated spectrum in cm-2 s-1 keV-1 - virtual TH1F* GetEnergySpectrum(Double_t m = 0) = 0; + virtual TH1D* GetEnergySpectrum() = 0; virtual TCanvas* DrawSolarFlux(); @@ -86,9 +92,8 @@ class TRestAxionSolarFlux : public TRestMetadata { Bool_t AreTablesLoaded() { return fTablesLoaded; } Double_t GetMass() { return fMass; } - void SetMass(const Double_t& m) { fMass = m; } - TH1F* GetFluxHistogram(std::string fname, Double_t binSize = 0.01); + TH1D* GetFluxHistogram(std::string fname, Double_t binSize = 0.01); TCanvas* DrawFluxFile(std::string fname, Double_t binSize = 0.01); void PrintMetadata(); diff --git a/inc/TRestAxionSolarHiddenPhotonFlux.h b/inc/TRestAxionSolarHiddenPhotonFlux.h index 876a29ef..2df02f96 100644 --- a/inc/TRestAxionSolarHiddenPhotonFlux.h +++ b/inc/TRestAxionSolarHiddenPhotonFlux.h @@ -41,20 +41,20 @@ class TRestAxionSolarHiddenPhotonFlux : public TRestAxionSolarFlux { /// It will be used when loading `.flux` files to define the input file energy binsize in eV. Double_t fBinSize = 0; //< - /// The tabulated solar flux continuum spectra TH1F(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius - std::vector fFluxTable; //! + /// The tabulated solar flux continuum spectra TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fFluxTable; //! - /// The tabulated solar flux continuum spectra TH1F(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius - std::vector fContinuumTable; //! + /// The tabulated solar flux continuum spectra TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector> fContinuumTable; //! - /// The tabulated resonance width TH1F(200,0,20)keV in eV2 versus solar radius - std::vector fWidthTable; //! + /// The tabulated resonance width TH1D(200,0,20)keV in eV2 versus solar radius + std::vector> fWidthTable; //! /// The solar plasma frequency vector in eV versus solar radius - std::vector fPlasmaFreqTable; //! + std::vector> fPlasmaFreqTable; //! - /// The total solar flux TH1F(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius - std::vector fFullFluxTable; //! + /// The total solar flux TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fFullFluxTable; //! /// Accumulative integrated solar flux for each solar ring for continuum spectrum (renormalized to unity) std::vector fFluxTableIntegrals; //! @@ -63,17 +63,17 @@ class TRestAxionSolarHiddenPhotonFlux : public TRestAxionSolarFlux { Double_t fTotalContinuumFlux = 0; //! /// A pointer to the continuum spectrum histogram - TH1F* fContinuumHist = nullptr; //! + TH1D* fContinuumHist = nullptr; //! /// A pointer to the superposed monochromatic and continuum spectrum histogram - TH1F* fTotalHist = nullptr; //! + TH1D* fTotalHist = nullptr; //! void ReadFluxFile(); void LoadContinuumFluxTable(); void LoadMonoChromaticFluxTable(); void IntegrateSolarFluxes(); - public: + public: /// It returns true if continuum flux spectra was loaded Bool_t isSolarTableLoaded() { return fFluxTable.size() > 0; } @@ -84,13 +84,13 @@ class TRestAxionSolarHiddenPhotonFlux : public TRestAxionSolarFlux { Bool_t isPlasmaFreqLoaded() { return fPlasmaFreqTable.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; + Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1)) override; /// It defines how to generate Monte Carlo energy and radius values to reproduce the flux std::pair GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1)) override; /// It defines how to read the solar tables at the inhereted class for a given mass in eV - Bool_t LoadTables(Double_t mass) override; + Bool_t LoadTables() override; void LoadContinuumTable(); void LoadWidthTable(); @@ -99,13 +99,13 @@ class TRestAxionSolarHiddenPhotonFlux : public TRestAxionSolarFlux { // calculate solar HP flux from the 3 tables and mass void CalculateSolarFlux(); /// It returns the total integrated flux at earth in cm-2 s-1 - Double_t GetTotalFlux(Double_t mass = 0) override { return fTotalContinuumFlux; } + Double_t GetTotalFlux() override { return fTotalContinuumFlux; } /// It returns an energy integrated spectrum in cm-2 s-1 keV-1 - TH1F* GetEnergySpectrum(Double_t m = 0) override { return GetTotalSpectrum(); } + TH1D* GetEnergySpectrum() override { return GetTotalSpectrum(); } - TH1F* GetContinuumSpectrum(); - TH1F* GetTotalSpectrum(); + TH1D* GetContinuumSpectrum(); + TH1D* GetTotalSpectrum(); virtual TCanvas* DrawSolarFlux() override; diff --git a/inc/TRestAxionSolarHiddenPhotonFlux.h.new b/inc/TRestAxionSolarHiddenPhotonFlux.h.new new file mode 100644 index 00000000..2df02f96 --- /dev/null +++ b/inc/TRestAxionSolarHiddenPhotonFlux.h.new @@ -0,0 +1,131 @@ +/************************************************************************* + * 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 _TRestAxionSolarHiddenPhotonFlux +#define _TRestAxionSolarHiddenPhotonFlux + +#include +#include + +//! A metadata class to load tabulated solar hidden photon fluxes. Kinetic mixing set to 1. +class TRestAxionSolarHiddenPhotonFlux : public TRestAxionSolarFlux { + private: + /// The filename containing the solar flux table with continuum spectrum + std::string fFluxDataFile = ""; //< + + /// The filename containing the resonance width (wGamma) + std::string fWidthDataFile = ""; //< + + /// The filename containing the plasma frequency (wp) table + std::string fPlasmaFreqDataFile = ""; //< + + /// It will be used when loading `.flux` files to define the input file energy binsize in eV. + Double_t fBinSize = 0; //< + + /// The tabulated solar flux continuum spectra TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fFluxTable; //! + + /// The tabulated solar flux continuum spectra TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector> fContinuumTable; //! + + /// The tabulated resonance width TH1D(200,0,20)keV in eV2 versus solar radius + std::vector> fWidthTable; //! + + /// The solar plasma frequency vector in eV versus solar radius + std::vector> fPlasmaFreqTable; //! + + /// The total solar flux TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fFullFluxTable; //! + + /// Accumulative integrated solar flux for each solar ring for continuum spectrum (renormalized to unity) + std::vector fFluxTableIntegrals; //! + + /// Total solar flux for monochromatic contributions + Double_t fTotalContinuumFlux = 0; //! + + /// A pointer to the continuum spectrum histogram + TH1D* fContinuumHist = nullptr; //! + + /// A pointer to the superposed monochromatic and continuum spectrum histogram + TH1D* 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 width table was loaded + Bool_t isWidthTableLoaded() { return fWidthTable.size() > 0; } + + /// It returns true if plasma frequency table was loaded + Bool_t isPlasmaFreqLoaded() { return fPlasmaFreqTable.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)) override; + + /// It defines how to generate Monte Carlo energy and radius values to reproduce the flux + std::pair GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1)) override; + + /// It defines how to read the solar tables at the inhereted class for a given mass in eV + Bool_t LoadTables() override; + + void LoadContinuumTable(); + void LoadWidthTable(); + void LoadPlasmaFreqTable(); + + // calculate solar HP flux from the 3 tables and mass + void CalculateSolarFlux(); + /// It returns the total integrated flux at earth in cm-2 s-1 + Double_t GetTotalFlux() override { return fTotalContinuumFlux; } + + /// It returns an energy integrated spectrum in cm-2 s-1 keV-1 + TH1D* GetEnergySpectrum() override { return GetTotalSpectrum(); } + + TH1D* GetContinuumSpectrum(); + TH1D* 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(); + + TRestAxionSolarHiddenPhotonFlux(); + TRestAxionSolarHiddenPhotonFlux(const char* cfgFileName, std::string name = ""); + ~TRestAxionSolarHiddenPhotonFlux(); + + ClassDefOverride(TRestAxionSolarHiddenPhotonFlux, 1); +}; +#endif diff --git a/inc/TRestAxionSolarHiddenPhotonFlux.h.old b/inc/TRestAxionSolarHiddenPhotonFlux.h.old new file mode 100644 index 00000000..bdc44040 --- /dev/null +++ b/inc/TRestAxionSolarHiddenPhotonFlux.h.old @@ -0,0 +1,131 @@ +/************************************************************************* + * 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 _TRestAxionSolarHiddenPhotonFlux +#define _TRestAxionSolarHiddenPhotonFlux + +#include +#include + +//! A metadata class to load tabulated solar hidden photon fluxes. Kinetic mixing set to 1. +class TRestAxionSolarHiddenPhotonFlux : public TRestAxionSolarFlux { + private: + /// The filename containing the solar flux table with continuum spectrum + std::string fFluxDataFile = ""; //< + + /// The filename containing the resonance width (wGamma) + std::string fWidthDataFile = ""; //< + + /// The filename containing the plasma frequency (wp) table + std::string fPlasmaFreqDataFile = ""; //< + + /// It will be used when loading `.flux` files to define the input file energy binsize in eV. + Double_t fBinSize = 0; //< + + /// The tabulated solar flux continuum spectra TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fFluxTable; //! + + /// The tabulated solar flux continuum spectra TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fContinuumTable; //! + + /// The tabulated resonance width TH1D(200,0,20)keV in eV2 versus solar radius + std::vector fWidthTable; //! + + /// The solar plasma frequency vector in eV versus solar radius + std::vector fPlasmaFreqTable; //! + + /// The total solar flux TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fFullFluxTable; //! + + /// Accumulative integrated solar flux for each solar ring for continuum spectrum (renormalized to unity) + std::vector fFluxTableIntegrals; //! + + /// Total solar flux for monochromatic contributions + Double_t fTotalContinuumFlux = 0; //! + + /// A pointer to the continuum spectrum histogram + TH1D* fContinuumHist = nullptr; //! + + /// A pointer to the superposed monochromatic and continuum spectrum histogram + TH1D* 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 width table was loaded + Bool_t isWidthTableLoaded() { return fWidthTable.size() > 0; } + + /// It returns true if plasma frequency table was loaded + Bool_t isPlasmaFreqLoaded() { return fPlasmaFreqTable.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)) override; + + /// It defines how to generate Monte Carlo energy and radius values to reproduce the flux + std::pair GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1)) override; + + /// It defines how to read the solar tables at the inhereted class for a given mass in eV + Bool_t LoadTables() override; + + void LoadContinuumTable(); + void LoadWidthTable(); + void LoadPlasmaFreqTable(); + + // calculate solar HP flux from the 3 tables and mass + void CalculateSolarFlux(); + /// It returns the total integrated flux at earth in cm-2 s-1 + Double_t GetTotalFlux() override { return fTotalContinuumFlux; } + + /// It returns an energy integrated spectrum in cm-2 s-1 keV-1 + TH1D* GetEnergySpectrum() override { return GetTotalSpectrum(); } + + TH1D* GetContinuumSpectrum(); + TH1D* 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(); + + TRestAxionSolarHiddenPhotonFlux(); + TRestAxionSolarHiddenPhotonFlux(const char* cfgFileName, std::string name = ""); + ~TRestAxionSolarHiddenPhotonFlux(); + + ClassDefOverride(TRestAxionSolarHiddenPhotonFlux, 1); +}; +#endif diff --git a/inc/TRestAxionSolarQCDFlux.h b/inc/TRestAxionSolarQCDFlux.h index 13246023..10db402e 100644 --- a/inc/TRestAxionSolarQCDFlux.h +++ b/inc/TRestAxionSolarQCDFlux.h @@ -41,11 +41,11 @@ class TRestAxionSolarQCDFlux : public TRestAxionSolarFlux { /// 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 fFluxTable; //! + /// The tabulated solar flux continuum spectra TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fFluxTable; //! /// The tabulated solar flux in cm-2 s-1 for a number of monochromatic energies versus solar radius - std::map fFluxLines; //! + std::map fFluxLines; //! /// Accumulative integrated solar flux for each solar ring for continuum spectrum (renormalized to unity) std::vector fFluxTableIntegrals; //! @@ -63,13 +63,13 @@ class TRestAxionSolarQCDFlux : public TRestAxionSolarFlux { Double_t fFluxRatio = 0; //! /// A pointer to the continuum spectrum histogram - TH1F* fContinuumHist = nullptr; //! + TH1D* fContinuumHist = nullptr; //! /// A pointer to the monochromatic spectrum histogram - TH1F* fMonoHist = nullptr; //! + TH1D* fMonoHist = nullptr; //! /// A pointer to the superposed monochromatic and continuum spectrum histogram - TH1F* fTotalHist = nullptr; //! + TH1D* fTotalHist = nullptr; //! void ReadFluxFile(); void LoadContinuumFluxTable(); @@ -84,25 +84,25 @@ class TRestAxionSolarQCDFlux : public TRestAxionSolarFlux { 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; + Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1)) override; /// It defines how to generate Monte Carlo energy and radius values to reproduce the flux std::pair GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1)) override; /// It defines how to read the solar tables at the inhereted class - Bool_t LoadTables(Double_t mass = 0) override; + 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 { + Double_t GetTotalFlux() 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(); } + TH1D* GetEnergySpectrum() override { return GetTotalSpectrum(); } - TH1F* GetContinuumSpectrum(); - TH1F* GetMonochromaticSpectrum(); - TH1F* GetTotalSpectrum(); + TH1D* GetContinuumSpectrum(); + TH1D* GetMonochromaticSpectrum(); + TH1D* GetTotalSpectrum(); virtual TCanvas* DrawSolarFlux() override; diff --git a/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py b/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py index d24dfcc1..c07ac580 100755 --- a/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py +++ b/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py @@ -44,11 +44,12 @@ parser.add_argument( "--N", dest="samples", type=int, help="The number of generated particles" ) -parser.add_argument("--mass", dest="mass", type=float, help="Hidden photon mass [eV]") - +parser.add_argument( + "--m", dest="mass", type=float, help="Hidden photon mass [eV]" +) args = parser.parse_args() -mass = 10 # eV +mass = 10. # eV if args.mass != None: mass = args.mass @@ -68,7 +69,7 @@ if args.samples != None: samples = args.samples -validation = True +validation = False ROOT.gSystem.Load("libRestFramework.so") ROOT.gSystem.Load("libRestAxion.so") @@ -82,7 +83,9 @@ pad1.Draw() combinedFlux = ROOT.TRestAxionSolarHiddenPhotonFlux(rmlfile, fluxname) -combinedFlux.Initialize(mass) +combinedFlux.SetMass(mass) +print(combinedFlux.GetMass()) +combinedFlux.Initialize() combinedFlux.PrintMetadata() if combinedFlux.GetError(): @@ -91,8 +94,9 @@ exit(101) comb_spt = TH2D("comb_spt", "Energy versus solar radius", 20000, 0, 20, 100, 0, 1) -for x in range(samples): +for i in range(samples): x = combinedFlux.GetRandomEnergyAndRadius((-1, -1)) + #print(x) comb_spt.Fill(x[0], x[1]) rnd = TRandom3(0) @@ -160,17 +164,4 @@ c1.Print(outfname) print("Generated file : " + outfname) -print("\nMaximum energy bin is " + str(enSpt.GetMaximumBin())) -if validation: - if enSpt.GetMaximumBin() != 8001: - 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 (25)! Exit code : 2") - exit(2) - -exit(0) +#exit(0) diff --git a/src/TRestAxionSolarFlux.cxx b/src/TRestAxionSolarFlux.cxx index a2e0ccd6..f82d9e34 100644 --- a/src/TRestAxionSolarFlux.cxx +++ b/src/TRestAxionSolarFlux.cxx @@ -64,7 +64,7 @@ /// on the flux distribution initialized. /// - **LoadTables**: It is called by TRestAxionSolarFlux::Initialize to allow the inherited /// class to load all the necessary tables in memory. -/// - **GetEnergySpectrum**: It should return a TH1F pointer with a energy spectrum histogram. +/// - **GetEnergySpectrum**: It should return a TH1D pointer with a energy spectrum histogram. /// ///-------------------------------------------------------------------------- /// @@ -125,12 +125,7 @@ void TRestAxionSolarFlux::Initialize() { fTablesLoaded = false; if (LoadTables()) fTablesLoaded = true; - if (!fRandom) { - delete fRandom; - fRandom = nullptr; - } - - if (fRandom != nullptr) { + if (fRandom) { delete fRandom; fRandom = nullptr; } @@ -139,11 +134,18 @@ void TRestAxionSolarFlux::Initialize() { fSeed = fRandom->TRandom::GetSeed(); } + +/////////////////////////////////////////////// +/// \brief Initialization of TRestAxionSolarFlux members with specific mass +/// +//void TRestAxionSolarFlux::InitializeMass( Double_t mass ) { SetMass(mass); RESTMetadata << GetMass() << RESTendl; } // SetMass calls Initialize + + /////////////////////////////////////////////// /// \brief It builds a histogram using the contents of the .flux file given /// in the argument. /// -TH1F* TRestAxionSolarFlux::GetFluxHistogram(string fname, Double_t binSize) { +TH1D* TRestAxionSolarFlux::GetFluxHistogram(string fname, Double_t binSize) { string fullPathName = SearchFile(fname); std::vector> fluxData; @@ -160,7 +162,7 @@ TH1F* TRestAxionSolarFlux::GetFluxHistogram(string fname, Double_t binSize) { originalHist->Fill(r, en, flux); } - return (TH1F*)originalHist->ProjectionY(); + return (TH1D*)originalHist->ProjectionY(); } /////////////////////////////////////////////// @@ -206,7 +208,7 @@ TCanvas* TRestAxionSolarFlux::DrawSolarFlux() { pad1->SetLeftMargin(0.15); pad1->SetBottomMargin(0.15); - TH1F* ht = GetEnergySpectrum(); + TH1D* ht = GetEnergySpectrum(); ht->SetLineColor(kBlack); ht->SetFillStyle(4050); ht->SetFillColor(kBlue - 10); diff --git a/src/TRestAxionSolarHiddenPhotonFlux.cxx b/src/TRestAxionSolarHiddenPhotonFlux.cxx index 77abd9b8..01f443fa 100644 --- a/src/TRestAxionSolarHiddenPhotonFlux.cxx +++ b/src/TRestAxionSolarHiddenPhotonFlux.cxx @@ -183,10 +183,25 @@ TRestAxionSolarHiddenPhotonFlux::~TRestAxionSolarHiddenPhotonFlux() {} /// \brief It will load the tables in memory by using the filename information provided /// inside the metadata members, and calculate the solar flux for a given m. /// -Bool_t TRestAxionSolarHiddenPhotonFlux::LoadTables(Double_t mass) { - if (fFluxDataFile == "" || fWidthDataFile == "" || fPlasmaFreqDataFile == "") return false; +Bool_t TRestAxionSolarHiddenPhotonFlux::LoadTables() { - SetMass(mass); + if (GetMass() <= 0 ) { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - hidden photon mass not yet defined" << RESTendl; + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - mass given as " << GetMass() << " eV" << RESTendl; + return false; + } + if ( fFluxDataFile == "" ){ + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - flux table not found!!\n " << fFluxDataFile << RESTendl; + return false; + } + if ( fWidthDataFile == "") { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - width table not found!!\n " << fWidthDataFile << RESTendl; + return false; + } + if ( fPlasmaFreqDataFile == "" ) { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - plasma frequency table not found!!\n " << fPlasmaFreqDataFile << RESTendl; + return false; + } LoadContinuumFluxTable(); LoadWidthTable(); @@ -195,7 +210,6 @@ Bool_t TRestAxionSolarHiddenPhotonFlux::LoadTables(Double_t mass) { CalculateSolarFlux(); IntegrateSolarFluxes(); - return true; } @@ -203,9 +217,9 @@ Bool_t TRestAxionSolarHiddenPhotonFlux::LoadTables(Double_t mass) { /// \brief A helper method to load the data file containing continuum spectra as a /// function of the solar radius. It will be called by TRestAxionSolarHiddenPhotonFlux::Initialize. /// -void TRestAxionSolarHiddenPhotonFlux::LoadContinuumTable() { +void TRestAxionSolarHiddenPhotonFlux::LoadContinuumFluxTable() { if (fFluxDataFile == "") { - RESTDebug + RESTError << "TRestAxionSolarHiddenPhotonFlux::LoadContinuumFluxTable. No solar flux table was defined" << RESTendl; return; @@ -216,26 +230,18 @@ void TRestAxionSolarHiddenPhotonFlux::LoadContinuumTable() { RESTDebug << "Loading table from file : " << RESTendl; RESTDebug << "File : " << fullPathName << RESTendl; - std::vector> fluxTable; - if (!TRestTools::IsBinaryFile(fFluxDataFile)) { - fluxTable.clear(); + fContinuumTable.clear(); RESTError << "File is not in binary format!" << RESTendl; } - TRestTools::ReadBinaryTable(fullPathName, fluxTable); + TRestTools::ReadBinaryTable(fullPathName, fContinuumTable); - if (fluxTable.size() != 1000 && fluxTable[0].size() != 200) { - fluxTable.clear(); + if (fContinuumTable.size() != 1000 || fContinuumTable[0].size() != 200) { RESTError << "LoadContinuumFluxTable. The table does not contain the right number of rows or columns" << RESTendl; RESTError << "Table will not be populated" << RESTendl; - } - - for (unsigned int n = 0; n < fluxTable.size(); n++) { - TH1F* h = new TH1F(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), "", 200, 0, 20); - for (unsigned int m = 0; m < fluxTable[n].size(); m++) h->SetBinContent(m + 1, fluxTable[n][m]); - fContinuumTable.push_back(h); + fContinuumTable.clear(); } } @@ -245,7 +251,7 @@ void TRestAxionSolarHiddenPhotonFlux::LoadContinuumTable() { /// void TRestAxionSolarHiddenPhotonFlux::LoadWidthTable() { if (fFluxDataFile == "") { - RESTDebug << "TRestAxionSolarHiddenPhotonFlux::LoadWidthTable. No width table was defined" + RESTError << "TRestAxionSolarHiddenPhotonFlux::LoadWidthTable. No width table was defined" << RESTendl; return; } @@ -255,26 +261,18 @@ void TRestAxionSolarHiddenPhotonFlux::LoadWidthTable() { RESTDebug << "Loading table from file : " << RESTendl; RESTDebug << "File : " << fullPathName << RESTendl; - std::vector> fluxTable; - if (!TRestTools::IsBinaryFile(fWidthDataFile)) { - fluxTable.clear(); + fWidthTable.clear(); RESTError << "File is not in binary format!" << RESTendl; } - TRestTools::ReadBinaryTable(fullPathName, fluxTable); - - if (fluxTable.size() != 1000 && fluxTable[0].size() != 200) { - fluxTable.clear(); + TRestTools::ReadBinaryTable(fullPathName, fWidthTable); + + if (fWidthTable.size() != 1000 || fWidthTable[0].size() != 200) { RESTError << "LoadWidthTable. The table does not contain the right number of rows or columns" << RESTendl; RESTError << "Table will not be populated" << RESTendl; - } - - for (unsigned int n = 0; n < fluxTable.size(); n++) { - TH1F* h = new TH1F(Form("%s_ResonanceWidthAtRadius%d", GetName(), n), "", 200, 0, 20); - for (unsigned int m = 0; m < fluxTable[n].size(); m++) h->SetBinContent(m + 1, fluxTable[n][m]); - fWidthTable.push_back(h); + fWidthTable.clear(); } } @@ -284,7 +282,7 @@ void TRestAxionSolarHiddenPhotonFlux::LoadWidthTable() { /// void TRestAxionSolarHiddenPhotonFlux::LoadPlasmaFreqTable() { if (fFluxDataFile == "") { - RESTDebug + RESTError << "TRestAxionSolarHiddenPhotonFlux::LoadPlasmaFreqTable. No plasma frequency table was defined" << RESTendl; return; @@ -295,26 +293,18 @@ void TRestAxionSolarHiddenPhotonFlux::LoadPlasmaFreqTable() { RESTDebug << "Loading table from file : " << RESTendl; RESTDebug << "File : " << fullPathName << RESTendl; - std::vector> fluxTable; - if (!TRestTools::IsBinaryFile(fWidthDataFile)) { - fluxTable.clear(); RESTError << "File is not in binary format!" << RESTendl; + fPlasmaFreqTable.clear(); } - TRestTools::ReadBinaryTable(fullPathName, fluxTable); - - if (fluxTable.size() != 1000 && fluxTable[0].size() != 1) { - fluxTable.clear(); - RESTError << "LoadWidthTable. The table does not contain the right number of rows or columns" + TRestTools::ReadBinaryTable(fullPathName, fPlasmaFreqTable); + + if (fPlasmaFreqTable.size() != 1000 || fPlasmaFreqTable[0].size() != 1) { + RESTError << "LoadPlasmaFreqTable. The table does not contain the right number of rows or columns" << RESTendl; RESTError << "Table will not be populated" << RESTendl; - } - - for (unsigned int n = 0; n < fluxTable.size(); n++) { - TH1F* h = new TH1F(Form("%s_PlasmaFreqAtRadius%d", GetName(), n), "", 1, 0, 20); - for (unsigned int m = 0; m < fluxTable[n].size(); m++) h->SetBinContent(m + 1, fluxTable[n][m]); - fPlasmaFreqTable.push_back(h); + fPlasmaFreqTable.clear(); } } @@ -326,35 +316,37 @@ void TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux() { if (GetMass() == 0) { RESTError << "CalculateSolarFlux. The hidden photon mass is set to zero!" << RESTendl; return; + } + if (fContinuumTable.size() == 0 ) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty flux table!" << RESTendl; + return; } + if (fPlasmaFreqTable.size() == 0 ) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty plasma freq table!" << RESTendl; + return; + } + if (fWidthTable.size() == 0 ) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty width table!" << RESTendl; + return; + } + Double_t mass = GetMass(); for (unsigned int n = 0; n < fContinuumTable.size(); n++) { // m4 * chi2 * wG * flux / ( (m2 - wp2)^2 + (w G)^2 ) - vector mass2Vector(200, pow(mass, 2)); - double wp = fPlasmaFreqTable[n]->GetBinContent(1); - vector wp2Vector(200, pow(wp, 2)); - vector weights(200, 1); - - TH1F* hMass = new TH1F("hMass", "hMass", 200, 0, 20); - TH1F* hWp = new TH1F("hWp", "hWp", 200, 0, 20); - TH1F* hWg2 = (TH1F*)fWidthTable[n]->Clone(); - hWg2->Multiply(hWg2); // (w G)^2 - - // hMass->FillN(200, mass2Vector); // m^2 hist - hMass->FillN(200, mass2Vector.data(), weights.data()); // m^2 hist - // hWp->FillN(200, wp2Vector); // wp^2 hist - hWp->FillN(200, wp2Vector.data(), weights.data()); // wp^2 hist - - hMass->Add(hWp, -1); // (m2 - wp2) - hMass->Multiply(hMass); // (m2 - wp2)^2 - hMass->Add(hWg2); // (m2 - wp2)^2 + (w G)^2 - - TH1F* h = (TH1F*)fWidthTable[n]->Clone(); - h->Multiply(fContinuumTable[n]); - h->Divide(hMass); - h->Scale(pow(mass, 4)); - + Double_t wp = fPlasmaFreqTable[n][0]; + vector wG = fWidthTable[n]; + vector flux = fContinuumTable[n]; + vector v; + for( unsigned int c; c < wG.size(); c++ ) { + Double_t d1 = ( wG[c] * flux[c] * pow(mass,4) ); // m4 * wG * flux + Double_t d2 = ( pow( pow(mass,2) - pow(wp,2) , 2 ) + pow(wG[c],2) ); // m4 * wG * flux / ( (m2 - wp2)^2 + (w G)^2 ) + v.push_back(d1/d2); + } + + TH1D* h = new TH1D(Form("%s_TotalFluxTable%d", GetName(), n), "", 200, 0, 20); + for (unsigned int c = 0; c < v.size(); c++) { h->SetBinContent(c + 1, v[c]); } + //for (unsigned int c = 0; c < v.size(); c++) { h->Fill(v[c]); } fFluxTable.push_back(h); } } @@ -363,13 +355,13 @@ void TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux() { /// \brief It builds a histogram with the continuum spectrum. /// The flux will be expressed in cm-2 s-1 keV-1. Binned in 100eV steps. /// -TH1F* TRestAxionSolarHiddenPhotonFlux::GetContinuumSpectrum() { +TH1D* TRestAxionSolarHiddenPhotonFlux::GetContinuumSpectrum() { if (fContinuumHist != nullptr) { delete fContinuumHist; fContinuumHist = nullptr; } - fContinuumHist = new TH1F("ContinuumHist", "", 200, 0, 20); + fContinuumHist = new TH1D("ContinuumHist", "", 200, 0, 20); for (const auto& x : fFluxTable) { fContinuumHist->Add(x); } @@ -389,15 +381,15 @@ TH1F* TRestAxionSolarHiddenPhotonFlux::GetContinuumSpectrum() { /// \brief Same as GetContinuumSpectrum, the flux will be /// expressed in cm-2 s-1 keV-1. Binned in 1eV steps. /// -TH1F* TRestAxionSolarHiddenPhotonFlux::GetTotalSpectrum() { - TH1F* hc = GetContinuumSpectrum(); +TH1D* TRestAxionSolarHiddenPhotonFlux::GetTotalSpectrum() { + TH1D* hc = GetContinuumSpectrum(); if (fTotalHist != nullptr) { delete fTotalHist; fTotalHist = nullptr; } - fTotalHist = new TH1F("fTotalHist", "", 20000, 0, 20); + fTotalHist = new TH1D("fTotalHist", "", 20000, 0, 20); for (int n = 0; n < hc->GetNbinsX(); n++) { for (int m = 0; m < 100; m++) { fTotalHist->SetBinContent(n * 100 + 1 + m, hc->GetBinContent(n + 1)); @@ -433,7 +425,7 @@ void TRestAxionSolarHiddenPhotonFlux::IntegrateSolarFluxes() { /////////////////////////////////////////////// /// \brief It returns the integrated flux at earth in cm-2 s-1 for the given energy range /// -Double_t TRestAxionSolarHiddenPhotonFlux::IntegrateFluxInRange(TVector2 eRange, Double_t mass) { +Double_t TRestAxionSolarHiddenPhotonFlux::IntegrateFluxInRange(TVector2 eRange) { if (eRange.X() == -1 && eRange.Y() == -1) { if (GetTotalFlux() == 0) IntegrateSolarFluxes(); return GetTotalFlux(); @@ -455,8 +447,9 @@ Double_t TRestAxionSolarHiddenPhotonFlux::IntegrateFluxInRange(TVector2 eRange, /// flux distributions defined inside the solar tables loaded in the class /// std::pair TRestAxionSolarHiddenPhotonFlux::GetRandomEnergyAndRadius(TVector2 eRange) { + std::pair result = {0, 0}; - if (!AreTablesLoaded()) return result; + if (!AreTablesLoaded()) { RESTWarning << "Tables not loaded!!" << RESTendl; return result; } Double_t rnd = fRandom->Rndm(); for (unsigned int r = 0; r < fFluxTableIntegrals.size(); r++) { if (rnd < fFluxTableIntegrals[r]) { @@ -464,7 +457,7 @@ std::pair TRestAxionSolarHiddenPhotonFlux::GetRandomEnergyAn if (eRange.X() != -1 && eRange.Y() != -1) { if (energy < eRange.X() || energy > eRange.Y()) return GetRandomEnergyAndRadius(eRange); } - Double_t radius = ((Double_t)r + fRandom->Rndm()) * 0.01; + Double_t radius = ((Double_t)r + fRandom->Rndm()) * 0.001; std::pair p = {energy, radius}; return p; } @@ -473,7 +466,7 @@ std::pair TRestAxionSolarHiddenPhotonFlux::GetRandomEnergyAn } /////////////////////////////////////////////// -/// \brief It prints on screen the table that has been loaded in memory +/// \brief It prints on screen the flux table that has been read from file /// void TRestAxionSolarHiddenPhotonFlux::PrintContinuumSolarTable() { cout << "Continuum solar flux table: " << endl; @@ -506,6 +499,9 @@ void TRestAxionSolarHiddenPhotonFlux::PrintIntegratedRingFlux() { void TRestAxionSolarHiddenPhotonFlux::PrintMetadata() { TRestAxionSolarFlux::PrintMetadata(); + RESTMetadata << " - Solar hidden photon datafile (flux) : " << fFluxDataFile << RESTendl; + RESTMetadata << " - Solar hidden photon datafile (width) : " << fWidthDataFile << RESTendl; + RESTMetadata << " - Solar hidden photon datafile (plasma freq) : " << fPlasmaFreqDataFile << RESTendl; RESTMetadata << "-------" << RESTendl; RESTMetadata << " - Total continuum flux : " << fTotalContinuumFlux << " cm-2 s-1" << RESTendl; RESTMetadata << "++++++++++++++++++" << RESTendl; @@ -532,9 +528,9 @@ void TRestAxionSolarHiddenPhotonFlux::ExportTables(Bool_t ascii) { } if (fFluxTable.size() > 0) { - std::vector> table; + std::vector> table; for (const auto& x : fFluxTable) { - std::vector row; + std::vector row; for (int n = 0; n < x->GetNbinsX(); n++) row.push_back(x->GetBinContent(n + 1)); table.push_back(row); @@ -568,7 +564,7 @@ TCanvas* TRestAxionSolarHiddenPhotonFlux::DrawSolarFlux() { pad1->cd(1)->SetLeftMargin(0.15); pad1->cd(1)->SetBottomMargin(0.15); - TH1F* ht = GetTotalSpectrum(); + TH1D* ht = GetTotalSpectrum(); ht->SetLineColor(kBlack); ht->SetFillStyle(4050); ht->SetFillColor(kBlue - 10); @@ -584,3 +580,4 @@ TCanvas* TRestAxionSolarHiddenPhotonFlux::DrawSolarFlux() { return fCanvas; } + diff --git a/src/TRestAxionSolarHiddenPhotonFlux.cxx.new b/src/TRestAxionSolarHiddenPhotonFlux.cxx.new new file mode 100644 index 00000000..26c782c1 --- /dev/null +++ b/src/TRestAxionSolarHiddenPhotonFlux.cxx.new @@ -0,0 +1,583 @@ +/******************** REST disclaimer *********************************** + * 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. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// TRestAxionSolarHiddenPhotonFlux will use a file in binary format to initialize +/// a solar flux table that will describe the solar hidden photon flux spectrum as a function +/// of the solar radius. +/// +/// This class loads the hidden photon flux that depends on the mass and kinetic mixing parameter. +/// For axion-like particle prodution independent of mass there is the class +/// TRestAxionSolarQCDFlux. Both classes are prototyped by a pure base class TRestAxionSolarFlux +/// that defines common methods used to evaluate the flux, and generate Monte-Carlo events inside +/// TRestAxionGeneratorProcess. +/// +/// ### Basic use +/// +/// Once the class has been initialized, the main use of this class will be provided +/// by the method TRestAxionSolarHiddenPhotonFlux::GetRandomEnergyAndRadius. This method will +/// return a random axion energy and position inside the solar radius following the +/// distributions given by the solar flux tables. +/// +/// Description of the specific parameters accepted by this metadata class. +/// - *fluxDataFile:* A table with 1000 rows representing the solar ring flux from the +/// center to the corona, and 200 columns representing the flux, measured in cm-2 s-1 keV-1, +/// for the range (0,20)keV in steps of 100eV. The table is provided as a binary table using +/// `.N200f` extension. +/// - *widthDataFile:* A table with 1000 rows representing the width of the hidden photon +/// resonant production (wG) for each solar ring from the center to the corona, and 200 columns +/// representing the width, measured in eV2, for the range (0,20)keV in steps of 100eV. The +/// table is provided as a binary table using `.N200f` extension. +/// - *plasmaFreqDataFile:* A table with 1000 rows and only 1 column representing the solar +/// plasma frequency (wp) for each solar ring from the center to the corona, measured in eV. The +/// table is provided as a binary table using `.N1f` extension. +/// +/// Pre-generated solar axion flux tables will be available at the +/// [axionlib-data](https://github.com/rest-for-physics/axionlib-data/tree/master) +/// repository. The different RML flux definitions used to load those tables +/// will be found at the +/// [fluxes.rml](https://github.com/rest-for-physics/axionlib-data/blob/master/solarFlux/fluxes.rml) +/// file found at the axionlib-data repository. +/// +/// Inside a local REST installation, the `fluxes.rml` file will be found at the REST +/// installation directory, and it will be located automatically by the +/// TRestMetadata::SearchFile method. +/// +/// ### A basic RML definition +/// +/// The following definition integrates an axion-photon component with a continuum +/// spectrum using a Primakoff production model, and a dummy spectrum file that +/// includes two monocrhomatic lines at different solar disk radius positions. +/// +/// \code +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// \endcode +/// +/// \warning When the flux is loaded manually inside the `restRoot` interactive +/// shell, or inside a macro or script, after metadata initialization, it is necessary +/// to call the method TRestAxionSolarHiddenPhotonFlux::LoadTables(mass) to trigger the tables +/// initialization. +/// +/// ### Performing MonteCarlo tests using pre-loaded tables +/// +/// In order to test the response of different solar flux definitions we may use the script +/// `solarPlot.py` found at `pipeline/metadata/solarFlux/`. This script will generate a +/// number of particles and it will assign to each particle an energy and solar disk +/// location with the help of the method TRestAxionSolarHiddenPhotonFlux::GetRandomEnergyAndRadius. +/// +/// \code +/// python3 solarPlotHiddenPhoton.py --fluxname HiddenPhoton --N 1000000 +/// \endcode +/// +/// By default, it will load the flux definition found at `fluxes.rml` from the +/// `axionlib-data` repository, and generate a `png` image with the resuts from the +/// Monte Carlo execution. +/// +/// \htmlonly \endhtmlonly +/// +/// ![Solar flux distributions MC-generated with TRestAxionSolarQCDFlux.](ABC_flux_MC.png) +/// +/// ### Exporting the solar flux tables +/// +/// On top of that, we will be able to export those tables to the TRestAxionSolarHiddenPhotonFlux +/// standard format to be used in later occasions. +/// +/// \code +/// TRestAxionSolarHiddenPhotonFlux *sFlux = new TRestAxionSolarHiddenPhotonFlux("fluxes.rml", +/// "HiddenPhoton") sFlux->Initialize() sFlux->ExportTables() +/// \endcode +/// +/// which will produce a binary table `.N200f` with the continuum flux. The filename root will be +/// extracted from the original `.flux` file. Optionally we may export the continuum flux to an +/// ASCII file by indicating it at the TRestAxionSolarHiddenPhotonFlux::ExportTables method call. +/// The files will be placed at the REST user space, at `$HOME/.rest/export/` directory. +/// +/// TODO Implement the method TRestAxionSolarQCDFlux::InitializeSolarTable using +/// a solar model description by TRestAxionSolarModel. +/// +/// TODO Perhaps it would be interesting to replace fFluxTable for a TH2D +/// +///-------------------------------------------------------------------------- +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2023-May: Specific methods extracted from TRestAxionSolarFlux +/// Javier Galan +/// 2023-June: TRestAxionSolarHiddenPhotonFlux created by editing TRestAxionSolarQCDFlux +/// Tomas O'Shea +/// +/// \class TRestAxionSolarHiddenPhotonFlux +/// \author Javier Galan +/// +///
+/// + +#include "TRestAxionSolarHiddenPhotonFlux.h" +using namespace std; + +ClassImp(TRestAxionSolarHiddenPhotonFlux); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestAxionSolarHiddenPhotonFlux::TRestAxionSolarHiddenPhotonFlux() : TRestAxionSolarFlux() {} + +/////////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// 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 name The name of the specific metadata. It will be used to find the +/// corresponding TRestAxionMagneticField section inside the RML. +/// +TRestAxionSolarHiddenPhotonFlux::TRestAxionSolarHiddenPhotonFlux(const char* cfgFileName, string name) + : TRestAxionSolarFlux(cfgFileName) { + LoadConfigFromFile(fConfigFileName, name); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestAxionSolarHiddenPhotonFlux::~TRestAxionSolarHiddenPhotonFlux() {} + +/////////////////////////////////////////////// +/// \brief It will load the tables in memory by using the filename information provided +/// inside the metadata members, and calculate the solar flux for a given m. +/// +Bool_t TRestAxionSolarHiddenPhotonFlux::LoadTables() { + + if (GetMass() <= 0 ) { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - hidden photon mass not yet defined" << RESTendl; + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - mass given as " << GetMass() << " eV" << RESTendl; + return false; + } + if ( fFluxDataFile == "" ){ + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - flux table not found!!\n " << fFluxDataFile << RESTendl; + return false; + } + if ( fWidthDataFile == "") { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - width table not found!!\n " << fWidthDataFile << RESTendl; + return false; + } + if ( fPlasmaFreqDataFile == "" ) { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - plasma frequency table not found!!\n " << fPlasmaFreqDataFile << RESTendl; + return false; + } + + LoadContinuumFluxTable(); + LoadWidthTable(); + LoadPlasmaFreqTable(); + + CalculateSolarFlux(); + + IntegrateSolarFluxes(); + return true; +} + +/////////////////////////////////////////////// +/// \brief A helper method to load the data file containing continuum spectra as a +/// function of the solar radius. It will be called by TRestAxionSolarHiddenPhotonFlux::Initialize. +/// +void TRestAxionSolarHiddenPhotonFlux::LoadContinuumFluxTable() { + if (fFluxDataFile == "") { + RESTError + << "TRestAxionSolarHiddenPhotonFlux::LoadContinuumFluxTable. No solar flux table was defined" + << RESTendl; + return; + } + + string fullPathName = SearchFile((string)fFluxDataFile); + + RESTDebug << "Loading table from file : " << RESTendl; + RESTDebug << "File : " << fullPathName << RESTendl; + + if (!TRestTools::IsBinaryFile(fFluxDataFile)) { + fContinuumTable.clear(); + RESTError << "File is not in binary format!" << RESTendl; + } + + TRestTools::ReadBinaryTable(fullPathName, fContinuumTable); + + if (fContinuumTable.size() != 1000 || fContinuumTable[0].size() != 200) { + RESTError << "LoadContinuumFluxTable. The table does not contain the right number of rows or columns" + << RESTendl; + RESTError << "Table will not be populated" << RESTendl; + fContinuumTable.clear(); + } +} + +/////////////////////////////////////////////// +/// \brief A helper method to load the data file containing resonance width as a +/// function of the solar radius. It will be called by TRestAxionSolarHiddenPhotonFlux::Initialize. +/// +void TRestAxionSolarHiddenPhotonFlux::LoadWidthTable() { + if (fFluxDataFile == "") { + RESTError << "TRestAxionSolarHiddenPhotonFlux::LoadWidthTable. No width table was defined" + << RESTendl; + return; + } + + string fullPathName = SearchFile((string)fWidthDataFile); + + RESTDebug << "Loading table from file : " << RESTendl; + RESTDebug << "File : " << fullPathName << RESTendl; + + if (!TRestTools::IsBinaryFile(fWidthDataFile)) { + fWidthTable.clear(); + RESTError << "File is not in binary format!" << RESTendl; + } + + TRestTools::ReadBinaryTable(fullPathName, fWidthTable); + + if (fWidthTable.size() != 1000 || fWidthTable[0].size() != 200) { + RESTError << "LoadWidthTable. The table does not contain the right number of rows or columns" + << RESTendl; + RESTError << "Table will not be populated" << RESTendl; + fWidthTable.clear(); + } +} + +/////////////////////////////////////////////// +/// \brief A helper method to load the data file containing resonance width as a +/// function of the solar radius. It will be called by TRestAxionSolarHiddenPhotonFlux::Initialize. +/// +void TRestAxionSolarHiddenPhotonFlux::LoadPlasmaFreqTable() { + if (fFluxDataFile == "") { + RESTError + << "TRestAxionSolarHiddenPhotonFlux::LoadPlasmaFreqTable. No plasma frequency table was defined" + << RESTendl; + return; + } + + string fullPathName = SearchFile((string)fPlasmaFreqDataFile); + + RESTDebug << "Loading table from file : " << RESTendl; + RESTDebug << "File : " << fullPathName << RESTendl; + + if (!TRestTools::IsBinaryFile(fWidthDataFile)) { + RESTError << "File is not in binary format!" << RESTendl; + fPlasmaFreqTable.clear(); + } + + TRestTools::ReadBinaryTable(fullPathName, fPlasmaFreqTable); + + if (fPlasmaFreqTable.size() != 1000 || fPlasmaFreqTable[0].size() != 1) { + RESTError << "LoadPlasmaFreqTable. The table does not contain the right number of rows or columns" + << RESTendl; + RESTError << "Table will not be populated" << RESTendl; + fPlasmaFreqTable.clear(); + } +} + +/////////////////////////////////////////////// +/// \brief A helper method to calculate the real solar flux spectrum from the 3 tables, the +/// and the hidden photon mass for chi=1. +/// +void TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux() { + if (GetMass() == 0) { + RESTError << "CalculateSolarFlux. The hidden photon mass is set to zero!" << RESTendl; + return; + } + if (fContinuumTable.size() == 0 ) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty flux table!" << RESTendl; + return; + } + if (fPlasmaFreqTable.size() == 0 ) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty plasma freq table!" << RESTendl; + return; + } + if (fWidthTable.size() == 0 ) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty width table!" << RESTendl; + return; + } + fFluxTable.clear(); + Double_t mass = GetMass(); + for (unsigned int n = 0; n < fContinuumTable.size(); n++) { + // m4 * chi2 * wG * flux / ( (m2 - wp2)^2 + (w G)^2 ) + + Double_t wp = fPlasmaFreqTable[n][0]; + vector wG = fWidthTable[n]; + vector flux = fContinuumTable[n]; + vector v; + for( unsigned int c; c < wG.size(); c++ ) { + Double_t d1 = ( wG[c] * flux[c] * pow(mass,4) ); // m4 * wG * flux + Double_t d2 = ( pow( pow(mass,2) - pow(wp,2) , 2 ) + pow(wG[c],2) ); // m4 * wG * flux / ( (m2 - wp2)^2 + (w G)^2 ) + v.push_back(d1/d2); + } + + TH1D* h = new TH1D(Form("%s_TotalFluxTable%d", GetName(), n), "", 200, 0, 20); + for (unsigned int c = 0; c < v.size(); c++) { h->SetBinContent(c + 1, v[c]); } + //for (unsigned int c = 0; c < v.size(); c++) { h->Fill(v[c]); } + fFluxTable.push_back(h); + } +} + +/////////////////////////////////////////////// +/// \brief It builds a histogram with the continuum spectrum. +/// The flux will be expressed in cm-2 s-1 keV-1. Binned in 100eV steps. +/// +TH1D* TRestAxionSolarHiddenPhotonFlux::GetContinuumSpectrum() { + if (fContinuumHist != nullptr) { + delete fContinuumHist; + fContinuumHist = nullptr; + } + + fContinuumHist = new TH1D("ContinuumHist", "", 200, 0, 20); + for (const auto& x : fFluxTable) { + fContinuumHist->Add(x); + } + + fContinuumHist->SetStats(0); + fContinuumHist->GetXaxis()->SetTitle("Energy [keV]"); + fContinuumHist->GetXaxis()->SetTitleSize(0.05); + fContinuumHist->GetXaxis()->SetLabelSize(0.05); + fContinuumHist->GetYaxis()->SetTitle("Flux [cm-2 s-1 keV-1]"); + fContinuumHist->GetYaxis()->SetTitleSize(0.05); + fContinuumHist->GetYaxis()->SetLabelSize(0.05); + + return fContinuumHist; +} + +/////////////////////////////////////////////// +/// \brief Same as GetContinuumSpectrum, the flux will be +/// expressed in cm-2 s-1 keV-1. Binned in 1eV steps. +/// +TH1D* TRestAxionSolarHiddenPhotonFlux::GetTotalSpectrum() { + TH1D* hc = GetContinuumSpectrum(); + + if (fTotalHist != nullptr) { + delete fTotalHist; + fTotalHist = nullptr; + } + + fTotalHist = new TH1D("fTotalHist", "", 20000, 0, 20); + for (int n = 0; n < hc->GetNbinsX(); n++) { + for (int m = 0; m < 100; m++) { + fTotalHist->SetBinContent(n * 100 + 1 + m, hc->GetBinContent(n + 1)); + } + } + + fTotalHist->SetStats(0); + fTotalHist->GetXaxis()->SetTitle("Energy [keV]"); + fTotalHist->GetXaxis()->SetTitleSize(0.05); + fTotalHist->GetXaxis()->SetLabelSize(0.05); + fTotalHist->GetYaxis()->SetTitle("Flux [cm-2 s-1 keV-1]"); + fTotalHist->GetYaxis()->SetTitleSize(0.05); + fTotalHist->GetYaxis()->SetLabelSize(0.05); + + return fTotalHist; +} + +/////////////////////////////////////////////// +/// \brief A helper method to initialize the internal class data members with the +/// integrated flux for each solar ring. It will be called by TRestAxionSolarHiddenPhotonFlux::Initialize. +/// +void TRestAxionSolarHiddenPhotonFlux::IntegrateSolarFluxes() { + fTotalContinuumFlux = 0.0; + for (unsigned int n = 0; n < fFluxTable.size(); n++) { + fTotalContinuumFlux += fFluxTable[n]->Integral() * 0.1; // We integrate in 100eV steps + fFluxTableIntegrals.push_back(fTotalContinuumFlux); + } + + for (unsigned int n = 0; n < fFluxTableIntegrals.size(); n++) + fFluxTableIntegrals[n] /= fTotalContinuumFlux; +} + +/////////////////////////////////////////////// +/// \brief It returns the integrated flux at earth in cm-2 s-1 for the given energy range +/// +Double_t TRestAxionSolarHiddenPhotonFlux::IntegrateFluxInRange(TVector2 eRange) { + if (eRange.X() == -1 && eRange.Y() == -1) { + if (GetTotalFlux() == 0) IntegrateSolarFluxes(); + return GetTotalFlux(); + } + + Double_t flux = 0; + fTotalContinuumFlux = 0.0; + for (unsigned int n = 0; n < fFluxTable.size(); n++) { + flux += fFluxTable[n]->Integral(fFluxTable[n]->FindFixBin(eRange.X()), + fFluxTable[n]->FindFixBin(eRange.Y())) * + 0.1; // We integrate in 100eV steps + } + + return flux; +} + +/////////////////////////////////////////////// +/// \brief It returns a random solar radius position and energy according to the +/// flux distributions defined inside the solar tables loaded in the class +/// +std::pair TRestAxionSolarHiddenPhotonFlux::GetRandomEnergyAndRadius(TVector2 eRange) { + + std::pair result = {0, 0}; + if (!AreTablesLoaded()) { RESTWarning << "Tables not loaded!!" << RESTendl; return result; } + Double_t rnd = fRandom->Rndm(); + for (unsigned int r = 0; r < fFluxTableIntegrals.size(); r++) { + if (rnd < fFluxTableIntegrals[r]) { + Double_t energy = fFluxTable[r]->GetRandom(); + if (eRange.X() != -1 && eRange.Y() != -1) { + if (energy < eRange.X() || energy > eRange.Y()) return GetRandomEnergyAndRadius(eRange); + } + Double_t radius = ((Double_t)r + fRandom->Rndm()) * 0.001; + std::pair p = {energy, radius}; + return p; + } + } + return result; +} + +/////////////////////////////////////////////// +/// \brief It prints on screen the flux table that has been read from file +/// +void TRestAxionSolarHiddenPhotonFlux::PrintContinuumSolarTable() { + cout << "Continuum solar flux table: " << endl; + cout << "--------------------------- " << endl; + for (unsigned int n = 0; n < fFluxTable.size(); n++) { + for (int m = 0; m < fFluxTable[n]->GetNbinsX(); m++) + cout << fFluxTable[n]->GetBinContent(m + 1) << "\t"; + cout << endl; + cout << endl; + } + cout << endl; +} + +/////////////////////////////////////////////// +/// \brief It prints on screen the integrated solar flux per solar ring +/// +void TRestAxionSolarHiddenPhotonFlux::PrintIntegratedRingFlux() { + cout << "Integrated solar flux per solar ring: " << endl; + cout << "--------------------------- " << endl; + /* + for (int n = 0; n < fFluxPerRadius.size(); n++) + cout << "n : " << n << " flux : " << fFluxPerRadius[n] << endl; + cout << endl; + */ +} + +/////////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestAxionSolarHiddenPhotonFlux +/// +void TRestAxionSolarHiddenPhotonFlux::PrintMetadata() { + TRestAxionSolarFlux::PrintMetadata(); + + RESTMetadata << " - Solar hidden photon datafile (flux) : " << fFluxDataFile << RESTendl; + RESTMetadata << " - Solar hidden photon datafile (width) : " << fWidthDataFile << RESTendl; + RESTMetadata << " - Solar hidden photon datafile (plasma freq) : " << fPlasmaFreqDataFile << RESTendl; + RESTMetadata << "-------" << RESTendl; + RESTMetadata << " - Total continuum flux : " << fTotalContinuumFlux << " cm-2 s-1" << RESTendl; + RESTMetadata << "++++++++++++++++++" << RESTendl; + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { + PrintContinuumSolarTable(); + PrintIntegratedRingFlux(); + } +} + +/////////////////////////////////////////////// +/// \brief It will create files with spectra to be used +/// in a later ocasion. +/// +void TRestAxionSolarHiddenPhotonFlux::ExportTables(Bool_t ascii) { + string rootFilename = TRestTools::GetFileNameRoot(fFluxDataFile); + + string path = REST_USER_PATH + "/export/"; + + if (!TRestTools::fileExists(path)) { + std::cout << "Creating path: " << path << std::endl; + int z = system(("mkdir -p " + path).c_str()); + if (z != 0) RESTError << "Could not create directory " << path << RESTendl; + } + + if (fFluxTable.size() > 0) { + std::vector> table; + for (const auto& x : fFluxTable) { + std::vector row; + for (int n = 0; n < x->GetNbinsX(); n++) row.push_back(x->GetBinContent(n + 1)); + + table.push_back(row); + } + + if (!ascii) + TRestTools::ExportBinaryTable(path + "/" + rootFilename + ".N200f", table); + else + TRestTools::ExportASCIITable(path + "/" + rootFilename + ".dat", table); + } +} + +/////////////////////////////////////////////// +/// \brief It draws the contents of a .flux file. This method just receives the +/// +TCanvas* TRestAxionSolarHiddenPhotonFlux::DrawSolarFlux() { + if (fCanvas != nullptr) { + delete fCanvas; + fCanvas = nullptr; + } + fCanvas = new TCanvas("canv", "This is the canvas title", 1200, 500); + fCanvas->Draw(); + + TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97); + pad1->Divide(2, 1); + pad1->Draw(); + + pad1->cd(1); + pad1->cd(1)->SetLogy(); + pad1->cd(1)->SetRightMargin(0.09); + pad1->cd(1)->SetLeftMargin(0.15); + pad1->cd(1)->SetBottomMargin(0.15); + + TH1D* ht = GetTotalSpectrum(); + ht->SetLineColor(kBlack); + ht->SetFillStyle(4050); + ht->SetFillColor(kBlue - 10); + + ht->Draw("hist"); + + pad1->cd(2); + pad1->cd(2)->SetRightMargin(0.09); + pad1->cd(2)->SetLeftMargin(0.15); + pad1->cd(2)->SetBottomMargin(0.15); + + ht->Draw("hist"); + + return fCanvas; +} + diff --git a/src/TRestAxionSolarHiddenPhotonFlux.cxx.old b/src/TRestAxionSolarHiddenPhotonFlux.cxx.old new file mode 100644 index 00000000..e61f0630 --- /dev/null +++ b/src/TRestAxionSolarHiddenPhotonFlux.cxx.old @@ -0,0 +1,613 @@ +/******************** REST disclaimer *********************************** + * 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. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// TRestAxionSolarHiddenPhotonFlux will use a file in binary format to initialize +/// a solar flux table that will describe the solar hidden photon flux spectrum as a function +/// of the solar radius. +/// +/// This class loads the hidden photon flux that depends on the mass and kinetic mixing parameter. +/// For axion-like particle prodution independent of mass there is the class +/// TRestAxionSolarQCDFlux. Both classes are prototyped by a pure base class TRestAxionSolarFlux +/// that defines common methods used to evaluate the flux, and generate Monte-Carlo events inside +/// TRestAxionGeneratorProcess. +/// +/// ### Basic use +/// +/// Once the class has been initialized, the main use of this class will be provided +/// by the method TRestAxionSolarHiddenPhotonFlux::GetRandomEnergyAndRadius. This method will +/// return a random axion energy and position inside the solar radius following the +/// distributions given by the solar flux tables. +/// +/// Description of the specific parameters accepted by this metadata class. +/// - *fluxDataFile:* A table with 1000 rows representing the solar ring flux from the +/// center to the corona, and 200 columns representing the flux, measured in cm-2 s-1 keV-1, +/// for the range (0,20)keV in steps of 100eV. The table is provided as a binary table using +/// `.N200f` extension. +/// - *widthDataFile:* A table with 1000 rows representing the width of the hidden photon +/// resonant production (wG) for each solar ring from the center to the corona, and 200 columns +/// representing the width, measured in eV2, for the range (0,20)keV in steps of 100eV. The +/// table is provided as a binary table using `.N200f` extension. +/// - *plasmaFreqDataFile:* A table with 1000 rows and only 1 column representing the solar +/// plasma frequency (wp) for each solar ring from the center to the corona, measured in eV. The +/// table is provided as a binary table using `.N1f` extension. +/// +/// Pre-generated solar axion flux tables will be available at the +/// [axionlib-data](https://github.com/rest-for-physics/axionlib-data/tree/master) +/// repository. The different RML flux definitions used to load those tables +/// will be found at the +/// [fluxes.rml](https://github.com/rest-for-physics/axionlib-data/blob/master/solarFlux/fluxes.rml) +/// file found at the axionlib-data repository. +/// +/// Inside a local REST installation, the `fluxes.rml` file will be found at the REST +/// installation directory, and it will be located automatically by the +/// TRestMetadata::SearchFile method. +/// +/// ### A basic RML definition +/// +/// The following definition integrates an axion-photon component with a continuum +/// spectrum using a Primakoff production model, and a dummy spectrum file that +/// includes two monocrhomatic lines at different solar disk radius positions. +/// +/// \code +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// \endcode +/// +/// \warning When the flux is loaded manually inside the `restRoot` interactive +/// shell, or inside a macro or script, after metadata initialization, it is necessary +/// to call the method TRestAxionSolarHiddenPhotonFlux::LoadTables(mass) to trigger the tables +/// initialization. +/// +/// ### Performing MonteCarlo tests using pre-loaded tables +/// +/// In order to test the response of different solar flux definitions we may use the script +/// `solarPlot.py` found at `pipeline/metadata/solarFlux/`. This script will generate a +/// number of particles and it will assign to each particle an energy and solar disk +/// location with the help of the method TRestAxionSolarHiddenPhotonFlux::GetRandomEnergyAndRadius. +/// +/// \code +/// python3 solarPlotHiddenPhoton.py --fluxname HiddenPhoton --N 1000000 +/// \endcode +/// +/// By default, it will load the flux definition found at `fluxes.rml` from the +/// `axionlib-data` repository, and generate a `png` image with the resuts from the +/// Monte Carlo execution. +/// +/// \htmlonly \endhtmlonly +/// +/// ![Solar flux distributions MC-generated with TRestAxionSolarQCDFlux.](ABC_flux_MC.png) +/// +/// ### Exporting the solar flux tables +/// +/// On top of that, we will be able to export those tables to the TRestAxionSolarHiddenPhotonFlux +/// standard format to be used in later occasions. +/// +/// \code +/// TRestAxionSolarHiddenPhotonFlux *sFlux = new TRestAxionSolarHiddenPhotonFlux("fluxes.rml", +/// "HiddenPhoton") sFlux->Initialize() sFlux->ExportTables() +/// \endcode +/// +/// which will produce a binary table `.N200f` with the continuum flux. The filename root will be +/// extracted from the original `.flux` file. Optionally we may export the continuum flux to an +/// ASCII file by indicating it at the TRestAxionSolarHiddenPhotonFlux::ExportTables method call. +/// The files will be placed at the REST user space, at `$HOME/.rest/export/` directory. +/// +/// TODO Implement the method TRestAxionSolarQCDFlux::InitializeSolarTable using +/// a solar model description by TRestAxionSolarModel. +/// +/// TODO Perhaps it would be interesting to replace fFluxTable for a TH2D +/// +///-------------------------------------------------------------------------- +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2023-May: Specific methods extracted from TRestAxionSolarFlux +/// Javier Galan +/// 2023-June: TRestAxionSolarHiddenPhotonFlux created by editing TRestAxionSolarQCDFlux +/// Tomas O'Shea +/// +/// \class TRestAxionSolarHiddenPhotonFlux +/// \author Javier Galan +/// +///
+/// + +#include "TRestAxionSolarHiddenPhotonFlux.h" +using namespace std; + +ClassImp(TRestAxionSolarHiddenPhotonFlux); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestAxionSolarHiddenPhotonFlux::TRestAxionSolarHiddenPhotonFlux() : TRestAxionSolarFlux() {} + +/////////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// 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 name The name of the specific metadata. It will be used to find the +/// corresponding TRestAxionMagneticField section inside the RML. +/// +TRestAxionSolarHiddenPhotonFlux::TRestAxionSolarHiddenPhotonFlux(const char* cfgFileName, string name) + : TRestAxionSolarFlux(cfgFileName) { + LoadConfigFromFile(fConfigFileName, name); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestAxionSolarHiddenPhotonFlux::~TRestAxionSolarHiddenPhotonFlux() {} + +/////////////////////////////////////////////// +/// \brief It will load the tables in memory by using the filename information provided +/// inside the metadata members, and calculate the solar flux for a given m. +/// +Bool_t TRestAxionSolarHiddenPhotonFlux::LoadTables() { + + if (GetMass() <= 0 ) { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - hidden photon mass not yet defined" << RESTendl; + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - mass given as " << GetMass() << " eV" << RESTendl; + return false; + } + if ( fFluxDataFile == "" ){ + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - flux table not found!!\n " << fFluxDataFile << RESTendl; + return false; + } + if ( fWidthDataFile == "") { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - width table not found!!\n " << fWidthDataFile << RESTendl; + return false; + } + if ( fPlasmaFreqDataFile == "" ) { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - plasma frequency table not found!!\n " << fPlasmaFreqDataFile << RESTendl; + return false; + } + + LoadContinuumFluxTable(); + LoadWidthTable(); + LoadPlasmaFreqTable(); + + CalculateSolarFlux(); + + IntegrateSolarFluxes(); + return true; +} + +/////////////////////////////////////////////// +/// \brief A helper method to load the data file containing continuum spectra as a +/// function of the solar radius. It will be called by TRestAxionSolarHiddenPhotonFlux::Initialize. +/// +void TRestAxionSolarHiddenPhotonFlux::LoadContinuumFluxTable() { + if (fFluxDataFile == "") { + RESTError + << "TRestAxionSolarHiddenPhotonFlux::LoadContinuumFluxTable. No solar flux table was defined" + << RESTendl; + return; + } + + string fullPathName = SearchFile((string)fFluxDataFile); + + RESTDebug << "Loading table from file : " << RESTendl; + RESTDebug << "File : " << fullPathName << RESTendl; + + std::vector> fluxTable; + + if (!TRestTools::IsBinaryFile(fFluxDataFile)) { + fluxTable.clear(); + RESTError << "File is not in binary format!" << RESTendl; + } + + TRestTools::ReadBinaryTable(fullPathName, fluxTable); + + if (fluxTable.size() != 1000 || fluxTable[0].size() != 200) { + fluxTable.clear(); + RESTError << "LoadContinuumFluxTable. The table does not contain the right number of rows or columns" + << RESTendl; + RESTError << "Table will not be populated" << RESTendl; + } + + for (unsigned int n = 0; n < fluxTable.size(); n++) { + TH1D* h = new TH1D(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), "", 200, 0, 20); + for (unsigned int m = 0; m < fluxTable[n].size(); m++) { h->SetBinContent(m + 1, fluxTable[n][m]); } + fContinuumTable.push_back(h); + } +} + +/////////////////////////////////////////////// +/// \brief A helper method to load the data file containing resonance width as a +/// function of the solar radius. It will be called by TRestAxionSolarHiddenPhotonFlux::Initialize. +/// +void TRestAxionSolarHiddenPhotonFlux::LoadWidthTable() { + if (fFluxDataFile == "") { + RESTError << "TRestAxionSolarHiddenPhotonFlux::LoadWidthTable. No width table was defined" + << RESTendl; + return; + } + + string fullPathName = SearchFile((string)fWidthDataFile); + + RESTDebug << "Loading table from file : " << RESTendl; + RESTDebug << "File : " << fullPathName << RESTendl; + + std::vector> fluxTable; + + if (!TRestTools::IsBinaryFile(fWidthDataFile)) { + fluxTable.clear(); + RESTError << "File is not in binary format!" << RESTendl; + } + + TRestTools::ReadBinaryTable(fullPathName, fluxTable); + //RESTMetadata << "Width table rows / columns: " << fluxTable.size() << " " << fluxTable[0].size() << RESTendl; + + if (fluxTable.size() != 1000 || fluxTable[0].size() != 200) { + fluxTable.clear(); + RESTError << "LoadWidthTable. The table does not contain the right number of rows or columns" + << RESTendl; + RESTError << "Table will not be populated" << RESTendl; + } + + for (unsigned int n = 0; n < fluxTable.size(); n++) { + TH1D* h = new TH1D(Form("%s_ResonanceWidthAtRadius%d", GetName(), n), "", 200, 0, 20); + for (unsigned int m = 0; m < fluxTable[n].size(); m++) { h->SetBinContent(m + 1, fluxTable[n][m]); } + fWidthTable.push_back(h); + } +} + +/////////////////////////////////////////////// +/// \brief A helper method to load the data file containing resonance width as a +/// function of the solar radius. It will be called by TRestAxionSolarHiddenPhotonFlux::Initialize. +/// +void TRestAxionSolarHiddenPhotonFlux::LoadPlasmaFreqTable() { + if (fFluxDataFile == "") { + RESTError + << "TRestAxionSolarHiddenPhotonFlux::LoadPlasmaFreqTable. No plasma frequency table was defined" + << RESTendl; + return; + } + + string fullPathName = SearchFile((string)fPlasmaFreqDataFile); + + RESTDebug << "Loading table from file : " << RESTendl; + RESTDebug << "File : " << fullPathName << RESTendl; + + std::vector> fluxTable; + + if (!TRestTools::IsBinaryFile(fWidthDataFile)) { + fluxTable.clear(); + RESTError << "File is not in binary format!" << RESTendl; + } + + TRestTools::ReadBinaryTable(fullPathName, fluxTable); + + if (fluxTable.size() != 1000 || fluxTable[0].size() != 1) { + fluxTable.clear(); + RESTError << "LoadPlasmaFreqTable. The table does not contain the right number of rows or columns" + << RESTendl; + RESTError << "Table will not be populated" << RESTendl; + } + + for (unsigned int n = 0; n < fluxTable.size(); n++) { + TH1D* h = new TH1D(Form("%s_PlasmaFreqAtRadius%d", GetName(), n), "", 1, 0, 20); + for (unsigned int m = 0; m < fluxTable[n].size(); m++) { h->SetBinContent(m + 1, fluxTable[n][m]); } + fPlasmaFreqTable.push_back(h); + } +} + +/////////////////////////////////////////////// +/// \brief A helper method to calculate the real solar flux spectrum from the 3 tables, the +/// and the hidden photon mass for chi=1. +/// +void TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux() { + if (GetMass() == 0) { + RESTError << "CalculateSolarFlux. The hidden photon mass is set to zero!" << RESTendl; + return; + } + if (fContinuumTable.size() == 0 ) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty flux table!" << RESTendl; + return; + } + if (fPlasmaFreqTable.size() == 0 ) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty plasma freq table!" << RESTendl; + return; + } + if (fWidthTable.size() == 0 ) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty width table!" << RESTendl; + return; + } + + Double_t mass = GetMass(); + cout << mass << endl; + for (unsigned int n = 0; n < fContinuumTable.size(); n++) { + // m4 * chi2 * wG * flux / ( (m2 - wp2)^2 + (w G)^2 ) + + Double_t wp = fPlasmaFreqTable[n]->GetBinContent(1); + TH1D* hMass = new TH1D(Form("%s_hMass%d", GetName(), n), "hMass", 200, 0, 20); + TH1D* hWg2 = (TH1D*)fWidthTable[n]->Clone(); + hWg2->Multiply(hWg2); // (w G)^2 + + for ( unsigned int c; c < 200; c++ ) { + Double_t wG = fWidthTable[n]->GetBinContent(c+1); + hMass->SetBinContent( c+1, pow(mass,-4) * ( pow( pow(mass,2) - pow(wp,2) , 2 )));// + pow(wG,2) ) ); // m2 + } + + hMass->Add(hWg2); // (m2 - wp2)^2 + (w G)^2 + + TH1D* h = (TH1D*)fWidthTable[n]->Clone(); // wG + h->Multiply(fContinuumTable[n]); // wG * flux + h->Divide(hMass); // wG * flux / ( (m2 - wp2)^2 + (w G)^2 ) + //h->Scale(pow(mass, 4)); // m4 * wG * flux / ( (m2 - wp2)^2 + (w G)^2 ) + + fFluxTable.push_back(h); + } +} + +/////////////////////////////////////////////// +/// \brief It builds a histogram with the continuum spectrum. +/// The flux will be expressed in cm-2 s-1 keV-1. Binned in 100eV steps. +/// +TH1D* TRestAxionSolarHiddenPhotonFlux::GetContinuumSpectrum() { + if (fContinuumHist != nullptr) { + delete fContinuumHist; + fContinuumHist = nullptr; + } + + fContinuumHist = new TH1D("ContinuumHist", "", 200, 0, 20); + for (const auto& x : fFluxTable) { + fContinuumHist->Add(x); + } + + fContinuumHist->SetStats(0); + fContinuumHist->GetXaxis()->SetTitle("Energy [keV]"); + fContinuumHist->GetXaxis()->SetTitleSize(0.05); + fContinuumHist->GetXaxis()->SetLabelSize(0.05); + fContinuumHist->GetYaxis()->SetTitle("Flux [cm-2 s-1 keV-1]"); + fContinuumHist->GetYaxis()->SetTitleSize(0.05); + fContinuumHist->GetYaxis()->SetLabelSize(0.05); + + return fContinuumHist; +} + +/////////////////////////////////////////////// +/// \brief Same as GetContinuumSpectrum, the flux will be +/// expressed in cm-2 s-1 keV-1. Binned in 1eV steps. +/// +TH1D* TRestAxionSolarHiddenPhotonFlux::GetTotalSpectrum() { + TH1D* hc = GetContinuumSpectrum(); + + if (fTotalHist != nullptr) { + delete fTotalHist; + fTotalHist = nullptr; + } + + fTotalHist = new TH1D("fTotalHist", "", 20000, 0, 20); + for (int n = 0; n < hc->GetNbinsX(); n++) { + for (int m = 0; m < 100; m++) { + fTotalHist->SetBinContent(n * 100 + 1 + m, hc->GetBinContent(n + 1)); + } + } + + fTotalHist->SetStats(0); + fTotalHist->GetXaxis()->SetTitle("Energy [keV]"); + fTotalHist->GetXaxis()->SetTitleSize(0.05); + fTotalHist->GetXaxis()->SetLabelSize(0.05); + fTotalHist->GetYaxis()->SetTitle("Flux [cm-2 s-1 keV-1]"); + fTotalHist->GetYaxis()->SetTitleSize(0.05); + fTotalHist->GetYaxis()->SetLabelSize(0.05); + + return fTotalHist; +} + +/////////////////////////////////////////////// +/// \brief A helper method to initialize the internal class data members with the +/// integrated flux for each solar ring. It will be called by TRestAxionSolarHiddenPhotonFlux::Initialize. +/// +void TRestAxionSolarHiddenPhotonFlux::IntegrateSolarFluxes() { + fTotalContinuumFlux = 0.0; + for (unsigned int n = 0; n < fFluxTable.size(); n++) { + fTotalContinuumFlux += fFluxTable[n]->Integral() * 0.1; // We integrate in 100eV steps + fFluxTableIntegrals.push_back(fTotalContinuumFlux); + } + + for (unsigned int n = 0; n < fFluxTableIntegrals.size(); n++) + fFluxTableIntegrals[n] /= fTotalContinuumFlux; +} + +/////////////////////////////////////////////// +/// \brief It returns the integrated flux at earth in cm-2 s-1 for the given energy range +/// +Double_t TRestAxionSolarHiddenPhotonFlux::IntegrateFluxInRange(TVector2 eRange) { + if (eRange.X() == -1 && eRange.Y() == -1) { + if (GetTotalFlux() == 0) IntegrateSolarFluxes(); + return GetTotalFlux(); + } + + Double_t flux = 0; + fTotalContinuumFlux = 0.0; + for (unsigned int n = 0; n < fFluxTable.size(); n++) { + flux += fFluxTable[n]->Integral(fFluxTable[n]->FindFixBin(eRange.X()), + fFluxTable[n]->FindFixBin(eRange.Y())) * + 0.1; // We integrate in 100eV steps + } + + return flux; +} + +/////////////////////////////////////////////// +/// \brief It returns a random solar radius position and energy according to the +/// flux distributions defined inside the solar tables loaded in the class +/// +std::pair TRestAxionSolarHiddenPhotonFlux::GetRandomEnergyAndRadius(TVector2 eRange) { + + std::pair result = {0, 0}; + if (!AreTablesLoaded()) { RESTWarning << "Tables not loaded!!" << RESTendl; return result; } + Double_t rnd = fRandom->Rndm(); + for (unsigned int r = 0; r < fFluxTableIntegrals.size(); r++) { + if (rnd < fFluxTableIntegrals[r]) { + Double_t energy = fFluxTable[r]->GetRandom(); + if (eRange.X() != -1 && eRange.Y() != -1) { + if (energy < eRange.X() || energy > eRange.Y()) return GetRandomEnergyAndRadius(eRange); + } + Double_t radius = ((Double_t)r + fRandom->Rndm()) * 0.001; + std::pair p = {energy, radius}; + return p; + } + } + return result; +} + +/////////////////////////////////////////////// +/// \brief It prints on screen the flux table that has been read from file +/// +void TRestAxionSolarHiddenPhotonFlux::PrintContinuumSolarTable() { + cout << "Continuum solar flux table: " << endl; + cout << "--------------------------- " << endl; + for (unsigned int n = 0; n < fFluxTable.size(); n++) { + for (int m = 0; m < fFluxTable[n]->GetNbinsX(); m++) + cout << fFluxTable[n]->GetBinContent(m + 1) << "\t"; + cout << endl; + cout << endl; + } + cout << endl; +} + +/////////////////////////////////////////////// +/// \brief It prints on screen the integrated solar flux per solar ring +/// +void TRestAxionSolarHiddenPhotonFlux::PrintIntegratedRingFlux() { + cout << "Integrated solar flux per solar ring: " << endl; + cout << "--------------------------- " << endl; + /* + for (int n = 0; n < fFluxPerRadius.size(); n++) + cout << "n : " << n << " flux : " << fFluxPerRadius[n] << endl; + cout << endl; + */ +} + +/////////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestAxionSolarHiddenPhotonFlux +/// +void TRestAxionSolarHiddenPhotonFlux::PrintMetadata() { + TRestAxionSolarFlux::PrintMetadata(); + + RESTMetadata << " - Solar hidden photon datafile (flux) : " << fFluxDataFile << RESTendl; + RESTMetadata << " - Solar hidden photon datafile (width) : " << fWidthDataFile << RESTendl; + RESTMetadata << " - Solar hidden photon datafile (plasma freq) : " << fPlasmaFreqDataFile << RESTendl; + RESTMetadata << "-------" << RESTendl; + RESTMetadata << " - Total continuum flux : " << fTotalContinuumFlux << " cm-2 s-1" << RESTendl; + RESTMetadata << "++++++++++++++++++" << RESTendl; + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { + PrintContinuumSolarTable(); + PrintIntegratedRingFlux(); + } +} + +/////////////////////////////////////////////// +/// \brief It will create files with spectra to be used +/// in a later ocasion. +/// +void TRestAxionSolarHiddenPhotonFlux::ExportTables(Bool_t ascii) { + string rootFilename = TRestTools::GetFileNameRoot(fFluxDataFile); + + string path = REST_USER_PATH + "/export/"; + + if (!TRestTools::fileExists(path)) { + std::cout << "Creating path: " << path << std::endl; + int z = system(("mkdir -p " + path).c_str()); + if (z != 0) RESTError << "Could not create directory " << path << RESTendl; + } + + if (fFluxTable.size() > 0) { + std::vector> table; + for (const auto& x : fFluxTable) { + std::vector row; + for (int n = 0; n < x->GetNbinsX(); n++) row.push_back(x->GetBinContent(n + 1)); + + table.push_back(row); + } + + if (!ascii) + TRestTools::ExportBinaryTable(path + "/" + rootFilename + ".N200f", table); + else + TRestTools::ExportASCIITable(path + "/" + rootFilename + ".dat", table); + } +} + +/////////////////////////////////////////////// +/// \brief It draws the contents of a .flux file. This method just receives the +/// +TCanvas* TRestAxionSolarHiddenPhotonFlux::DrawSolarFlux() { + if (fCanvas != nullptr) { + delete fCanvas; + fCanvas = nullptr; + } + fCanvas = new TCanvas("canv", "This is the canvas title", 1200, 500); + fCanvas->Draw(); + + TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97); + pad1->Divide(2, 1); + pad1->Draw(); + + pad1->cd(1); + pad1->cd(1)->SetLogy(); + pad1->cd(1)->SetRightMargin(0.09); + pad1->cd(1)->SetLeftMargin(0.15); + pad1->cd(1)->SetBottomMargin(0.15); + + TH1D* ht = GetTotalSpectrum(); + ht->SetLineColor(kBlack); + ht->SetFillStyle(4050); + ht->SetFillColor(kBlue - 10); + + ht->Draw("hist"); + + pad1->cd(2); + pad1->cd(2)->SetRightMargin(0.09); + pad1->cd(2)->SetLeftMargin(0.15); + pad1->cd(2)->SetBottomMargin(0.15); + + ht->Draw("hist"); + + return fCanvas; +} + diff --git a/src/TRestAxionSolarQCDFlux.cxx b/src/TRestAxionSolarQCDFlux.cxx index 457bd31f..79b5452e 100644 --- a/src/TRestAxionSolarQCDFlux.cxx +++ b/src/TRestAxionSolarQCDFlux.cxx @@ -249,7 +249,7 @@ TRestAxionSolarQCDFlux::~TRestAxionSolarQCDFlux() {} /// \brief It will load the tables in memory by using the filename information provided /// inside the metadata members. /// -Bool_t TRestAxionSolarQCDFlux::LoadTables(Double_t mass) { +Bool_t TRestAxionSolarQCDFlux::LoadTables() { if (fFluxDataFile == "" && fFluxSptFile == "") return false; if (TRestTools::GetFileNameExtension(fFluxDataFile) == "flux") { @@ -265,7 +265,7 @@ Bool_t TRestAxionSolarQCDFlux::LoadTables(Double_t mass) { } /////////////////////////////////////////////// -/// \brief A helper method to load the data file containning continuum spectra as a +/// \brief A helper method to load the data file containing continuum spectra as a /// function of the solar radius. It will be called by TRestAxionSolarQCDFlux::Initialize. /// void TRestAxionSolarQCDFlux::LoadContinuumFluxTable() { @@ -280,7 +280,7 @@ void TRestAxionSolarQCDFlux::LoadContinuumFluxTable() { RESTDebug << "Loading table from file : " << RESTendl; RESTDebug << "File : " << fullPathName << RESTendl; - std::vector> fluxTable; + std::vector> fluxTable; if (TRestTools::GetFileNameExtension(fFluxDataFile) == "dat") { std::vector> doubleTable; if (!TRestTools::ReadASCIITable(fullPathName, doubleTable)) { @@ -289,8 +289,8 @@ void TRestAxionSolarQCDFlux::LoadContinuumFluxTable() { return; } for (const auto& row : doubleTable) { - std::vector floatVec(row.begin(), row.end()); - fluxTable.push_back(floatVec); + std::vector doubleVec(row.begin(), row.end()); + fluxTable.push_back(doubleVec); } } else if (TRestTools::IsBinaryFile(fFluxDataFile)) TRestTools::ReadBinaryTable(fullPathName, fluxTable); @@ -309,14 +309,14 @@ void TRestAxionSolarQCDFlux::LoadContinuumFluxTable() { } for (unsigned int n = 0; n < fluxTable.size(); n++) { - TH1F* h = new TH1F(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), "", 200, 0, 20); + TH1D* h = new TH1D(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), "", 200, 0, 20); for (unsigned int m = 0; m < fluxTable[n].size(); m++) h->SetBinContent(m + 1, fluxTable[n][m]); fFluxTable.push_back(h); } } /////////////////////////////////////////////// -/// \brief A helper method to load the data file containning monochromatic spectral +/// \brief A helper method to load the data file containing monochromatic spectral /// lines as a function of the solar radius. It will be called by TRestAxionSolarQCDFlux::Initialize. /// void TRestAxionSolarQCDFlux::LoadMonoChromaticFluxTable() { @@ -339,10 +339,10 @@ void TRestAxionSolarQCDFlux::LoadMonoChromaticFluxTable() { return; } - std::vector> asciiTable; + std::vector> asciiTable; for (const auto& row : doubleTable) { - std::vector floatVec(row.begin(), row.end()); - asciiTable.push_back(floatVec); + std::vector doubleVec(row.begin(), row.end()); + asciiTable.push_back(doubleVec); } fFluxLines.clear(); @@ -355,8 +355,8 @@ void TRestAxionSolarQCDFlux::LoadMonoChromaticFluxTable() { } for (unsigned int en = 0; en < asciiTable[0].size(); en++) { - Float_t energy = asciiTable[0][en]; - TH1F* h = new TH1F(Form("%s_MonochromeFluxAtEnergy%6.4lf", GetName(), energy), "", 100, 0, 1); + Double_t energy = asciiTable[0][en]; + TH1D* h = new TH1D(Form("%s_MonochromeFluxAtEnergy%6.4lf", GetName(), energy), "", 100, 0, 1); for (unsigned int r = 1; r < asciiTable.size(); r++) h->SetBinContent(r, asciiTable[r][en]); fFluxLines[energy] = h; } @@ -399,12 +399,12 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { Double_t fluxBinSize = TRestTools::GetLowestIncreaseFromTable(fluxData, 1); for (const auto& data : fluxData) { - Float_t r = 0.005 + data[0]; - Float_t en = data[1] - 0.005; - Float_t flux = data[2] * fluxBinSize; // flux in cm-2 s-1 bin-1 + Double_t r = 0.005 + data[0]; + Double_t en = data[1] - 0.005; + Double_t flux = data[2] * fluxBinSize; // flux in cm-2 s-1 bin-1 - originalHist->Fill(r, en, (Float_t)flux); - continuumHist->Fill(r, en, (Float_t)flux); + originalHist->Fill(r, en, (Double_t)flux); + continuumHist->Fill(r, en, (Double_t)flux); } RESTDebug << "Histograms filled" << RESTendl; @@ -415,8 +415,8 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { const int smearPoints = (Int_t)(5 / (fBinSize * 100)); const int excludePoints = smearPoints / 5; for (const auto& data : fluxData) { - Float_t r = 0.005 + data[0]; - Float_t en = data[1] - 0.005; + Double_t r = 0.005 + data[0]; + Double_t en = data[1] - 0.005; Int_t binR = continuumHist->GetXaxis()->FindBin(r); Int_t binE = continuumHist->GetYaxis()->FindBin(en); @@ -430,8 +430,8 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { } avgFlux /= n; - Float_t targetBinFlux = continuumHist->GetBinContent(binR, binE); - Float_t thrFlux = avgFlux + fPeakSigma * TMath::Sqrt(avgFlux); + Double_t targetBinFlux = continuumHist->GetBinContent(binR, binE); + Double_t thrFlux = avgFlux + fPeakSigma * TMath::Sqrt(avgFlux); if (targetBinFlux > 0 && targetBinFlux > thrFlux) { continuumHist->SetBinContent(binR, binE, avgFlux); peaks++; @@ -441,8 +441,8 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { for (int n = 0; n < originalHist->GetNbinsX(); n++) for (int m = 0; m < originalHist->GetNbinsY(); m++) { - Float_t orig = originalHist->GetBinContent(n + 1, m + 1); - Float_t cont = continuumHist->GetBinContent(n + 1, m + 1); + Double_t orig = originalHist->GetBinContent(n + 1, m + 1); + Double_t cont = continuumHist->GetBinContent(n + 1, m + 1); spectrumHist->SetBinContent(n + 1, m + 1, orig - cont); } @@ -453,8 +453,8 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { fFluxTable.clear(); for (int n = 0; n < continuumHist->GetNbinsX(); n++) { - TH1F* hc = - (TH1F*)continuumHist->ProjectionY(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), n + 1, n + 1); + TH1D* hc = + (TH1D*)continuumHist->ProjectionY(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), n + 1, n + 1); fFluxTable.push_back(hc); } @@ -462,7 +462,7 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { for (int n = 0; n < spectrumHist->GetNbinsY(); n++) { if (spectrumHist->ProjectionX("", n + 1, n + 1)->Integral() > 0) { Double_t energy = spectrumHist->ProjectionY()->GetBinCenter(n + 1); - TH1F* hm = (TH1F*)spectrumHist->ProjectionX( + TH1D* hm = (TH1D*)spectrumHist->ProjectionX( Form("%s_MonochromeFluxAtEnergy%5.3lf", GetName(), energy), n + 1, n + 1); fFluxLines[energy] = hm; } @@ -475,13 +475,13 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { /// \brief It builds a histogram with the continuum spectrum component. /// The flux will be expressed in cm-2 s-1 keV-1. Binned in 100eV steps. /// -TH1F* TRestAxionSolarQCDFlux::GetContinuumSpectrum() { +TH1D* TRestAxionSolarQCDFlux::GetContinuumSpectrum() { if (fContinuumHist != nullptr) { delete fContinuumHist; fContinuumHist = nullptr; } - fContinuumHist = new TH1F("ContinuumHist", "", 200, 0, 20); + fContinuumHist = new TH1D("ContinuumHist", "", 200, 0, 20); for (const auto& x : fFluxTable) { fContinuumHist->Add(x); } @@ -501,13 +501,13 @@ TH1F* TRestAxionSolarQCDFlux::GetContinuumSpectrum() { /// \brief It builds a histogram with the monochromatic spectrum component. /// The flux will be expressed in cm-2 s-1 eV-1. Binned in 1eV steps. /// -TH1F* TRestAxionSolarQCDFlux::GetMonochromaticSpectrum() { +TH1D* TRestAxionSolarQCDFlux::GetMonochromaticSpectrum() { if (fMonoHist != nullptr) { delete fMonoHist; fMonoHist = nullptr; } - fMonoHist = new TH1F("MonochromaticHist", "", 20000, 0, 20); + fMonoHist = new TH1D("MonochromaticHist", "", 20000, 0, 20); for (const auto& x : fFluxLines) { fMonoHist->Fill(x.first, x.second->Integral()); // cm-2 s-1 eV-1 } @@ -528,16 +528,16 @@ TH1F* TRestAxionSolarQCDFlux::GetMonochromaticSpectrum() { /// spectrum component. The flux will be expressed in cm-2 s-1 keV-1. /// Binned in 1eV steps. /// -TH1F* TRestAxionSolarQCDFlux::GetTotalSpectrum() { - TH1F* hm = GetMonochromaticSpectrum(); - TH1F* hc = GetContinuumSpectrum(); +TH1D* TRestAxionSolarQCDFlux::GetTotalSpectrum() { + TH1D* hm = GetMonochromaticSpectrum(); + TH1D* hc = GetContinuumSpectrum(); if (fTotalHist != nullptr) { delete fTotalHist; fTotalHist = nullptr; } - fTotalHist = new TH1F("fTotalHist", "", 20000, 0, 20); + fTotalHist = new TH1D("fTotalHist", "", 20000, 0, 20); for (int n = 0; n < hc->GetNbinsX(); n++) { for (int m = 0; m < 100; m++) { fTotalHist->SetBinContent(n * 100 + 1 + m, hc->GetBinContent(n + 1)); @@ -580,12 +580,12 @@ TCanvas* TRestAxionSolarQCDFlux::DrawSolarFlux() { pad1->cd(1)->SetLeftMargin(0.15); pad1->cd(1)->SetBottomMargin(0.15); - TH1F* ht = GetTotalSpectrum(); + TH1D* ht = GetTotalSpectrum(); ht->SetLineColor(kBlack); ht->SetFillStyle(4050); ht->SetFillColor(kBlue - 10); - TH1F* hm = GetMonochromaticSpectrum(); + TH1D* hm = GetMonochromaticSpectrum(); hm->SetLineColor(kBlack); hm->Scale(100); // renormalizing per 100eV-1 @@ -634,7 +634,7 @@ void TRestAxionSolarQCDFlux::IntegrateSolarFluxes() { /////////////////////////////////////////////// /// \brief It returns the integrated flux at earth in cm-2 s-1 for the given energy range /// -Double_t TRestAxionSolarQCDFlux::IntegrateFluxInRange(TVector2 eRange, Double_t mass) { +Double_t TRestAxionSolarQCDFlux::IntegrateFluxInRange(TVector2 eRange) { if (eRange.X() == -1 && eRange.Y() == -1) { if (GetTotalFlux() == 0) IntegrateSolarFluxes(); return GetTotalFlux(); @@ -770,9 +770,9 @@ void TRestAxionSolarQCDFlux::ExportTables(Bool_t ascii) { } if (fFluxTable.size() > 0) { - std::vector> table; + std::vector> table; for (const auto& x : fFluxTable) { - std::vector row; + std::vector row; for (int n = 0; n < x->GetNbinsX(); n++) row.push_back(x->GetBinContent(n + 1)); table.push_back(row); @@ -785,9 +785,9 @@ void TRestAxionSolarQCDFlux::ExportTables(Bool_t ascii) { } if (fFluxLines.size() > 0) { - std::vector> table; + std::vector> table; for (const auto& x : fFluxLines) { - std::vector row; + std::vector row; row.push_back(x.first); for (int n = 0; n < x.second->GetNbinsX(); n++) row.push_back(x.second->GetBinContent(n + 1));