diff --git a/inc/TRestAxionSolarHiddenPhotonFlux.h b/inc/TRestAxionSolarHiddenPhotonFlux.h index 2df02f96..bdc44040 100644 --- a/inc/TRestAxionSolarHiddenPhotonFlux.h +++ b/inc/TRestAxionSolarHiddenPhotonFlux.h @@ -45,13 +45,13 @@ class TRestAxionSolarHiddenPhotonFlux : public TRestAxionSolarFlux { 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; //! + std::vector fContinuumTable; //! /// The tabulated resonance width TH1D(200,0,20)keV in eV2 versus solar radius - std::vector> fWidthTable; //! + std::vector fWidthTable; //! /// The solar plasma frequency vector in eV versus solar radius - std::vector> fPlasmaFreqTable; //! + 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; //! diff --git a/inc/TRestAxionSolarHiddenPhotonFlux.h.new b/inc/TRestAxionSolarHiddenPhotonFlux.h.new deleted file mode 100644 index 2df02f96..00000000 --- a/inc/TRestAxionSolarHiddenPhotonFlux.h.new +++ /dev/null @@ -1,131 +0,0 @@ -/************************************************************************* - * 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 deleted file mode 100644 index bdc44040..00000000 --- a/inc/TRestAxionSolarHiddenPhotonFlux.h.old +++ /dev/null @@ -1,131 +0,0 @@ -/************************************************************************* - * 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/src/TRestAxionSolarHiddenPhotonFlux.cxx b/src/TRestAxionSolarHiddenPhotonFlux.cxx index 01f443fa..0c213f20 100644 --- a/src/TRestAxionSolarHiddenPhotonFlux.cxx +++ b/src/TRestAxionSolarHiddenPhotonFlux.cxx @@ -230,18 +230,26 @@ void TRestAxionSolarHiddenPhotonFlux::LoadContinuumFluxTable() { RESTDebug << "Loading table from file : " << RESTendl; RESTDebug << "File : " << fullPathName << RESTendl; + std::vector> fluxTable; + if (!TRestTools::IsBinaryFile(fFluxDataFile)) { - fContinuumTable.clear(); + fluxTable.clear(); RESTError << "File is not in binary format!" << RESTendl; } - TRestTools::ReadBinaryTable(fullPathName, fContinuumTable); + TRestTools::ReadBinaryTable(fullPathName, fluxTable); - if (fContinuumTable.size() != 1000 || fContinuumTable[0].size() != 200) { + 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; - fContinuumTable.clear(); + } + + 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); } } @@ -261,18 +269,27 @@ void TRestAxionSolarHiddenPhotonFlux::LoadWidthTable() { RESTDebug << "Loading table from file : " << RESTendl; RESTDebug << "File : " << fullPathName << RESTendl; + std::vector> fluxTable; + if (!TRestTools::IsBinaryFile(fWidthDataFile)) { - fWidthTable.clear(); + fluxTable.clear(); RESTError << "File is not in binary format!" << RESTendl; } - TRestTools::ReadBinaryTable(fullPathName, fWidthTable); + TRestTools::ReadBinaryTable(fullPathName, fluxTable); + //RESTMetadata << "Width table rows / columns: " << fluxTable.size() << " " << fluxTable[0].size() << RESTendl; - if (fWidthTable.size() != 1000 || fWidthTable[0].size() != 200) { + 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; - fWidthTable.clear(); + } + + 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); } } @@ -293,18 +310,26 @@ 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, fPlasmaFreqTable); + TRestTools::ReadBinaryTable(fullPathName, fluxTable); - if (fPlasmaFreqTable.size() != 1000 || fPlasmaFreqTable[0].size() != 1) { + 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; - fPlasmaFreqTable.clear(); + } + + 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); } } @@ -331,22 +356,27 @@ void TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux() { } 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][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); + 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 = 0; 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 ) - 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); } } diff --git a/src/TRestAxionSolarHiddenPhotonFlux.cxx.new b/src/TRestAxionSolarHiddenPhotonFlux.cxx.new deleted file mode 100644 index 26c782c1..00000000 --- a/src/TRestAxionSolarHiddenPhotonFlux.cxx.new +++ /dev/null @@ -1,583 +0,0 @@ -/******************** 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 deleted file mode 100644 index e61f0630..00000000 --- a/src/TRestAxionSolarHiddenPhotonFlux.cxx.old +++ /dev/null @@ -1,613 +0,0 @@ -/******************** 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; -} -