diff --git a/.gitignore b/.gitignore index f57a1172..9632e928 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ build *.root **/.DS_Store +*.json diff --git a/data b/data deleted file mode 160000 index 957d9fc5..00000000 --- a/data +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 957d9fc5e24d2e5aef97cfd1f99982f2c3a07a98 diff --git a/data b/data new file mode 120000 index 00000000..95fc0d20 --- /dev/null +++ b/data @@ -0,0 +1 @@ +../axionlib-data \ No newline at end of file diff --git a/flux.png b/flux.png new file mode 100644 index 00000000..26f7d1e5 Binary files /dev/null and b/flux.png differ diff --git a/inc/TRestAxionSolarFlux.h b/inc/TRestAxionSolarFlux.h index 73bcc00d..1d068863 100644 --- a/inc/TRestAxionSolarFlux.h +++ b/inc/TRestAxionSolarFlux.h @@ -38,6 +38,9 @@ class TRestAxionSolarFlux : public TRestMetadata { /// Axion coupling strength Double_t fCouplingStrength = 0; //< + /// Mass parameter + Double_t fMass = 0; //! + /// Seed used in random generator Int_t fSeed = 0; //< @@ -55,7 +58,7 @@ 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() = 0; + virtual Bool_t LoadTables(Double_t mass = 0) = 0; public: /// It is required in order to load solar flux tables into memory @@ -83,6 +86,9 @@ 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); TCanvas* DrawFluxFile(std::string fname, Double_t binSize = 0.01); diff --git a/inc/TRestAxionSolarHiddenPhotonFlux.h b/inc/TRestAxionSolarHiddenPhotonFlux.h new file mode 100644 index 00000000..502983af --- /dev/null +++ b/inc/TRestAxionSolarHiddenPhotonFlux.h @@ -0,0 +1,132 @@ +/************************************************************************* + * 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 TH1F(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 resonance width TH1F(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 TH1F(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 + TH1F* fContinuumHist = nullptr; //! + + /// A pointer to the superposed monochromatic and continuum spectrum histogram + TH1F* fTotalHist = nullptr; //! + + void ReadFluxFile(); + void LoadContinuumFluxTable(); + void LoadMonoChromaticFluxTable(); + void IntegrateSolarFluxes(); + + public: + /// It returns true if continuum flux spectra was loaded + Bool_t isSolarTableLoaded() { return fFluxTable.size() > 0; } + + /// It returns true if 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), Double_t mass = 0) override; + + /// It defines how to generate Monte Carlo energy and radius values to reproduce the flux + std::pair GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1), + Double_t mass = 0) 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; + + 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(Double_t mass = 0) 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(); } + + TH1F* GetContinuumSpectrum(); + TH1F* GetTotalSpectrum(); + + virtual TCanvas* DrawSolarFlux() override; + + /// Tables might be loaded using a solar model description by TRestAxionSolarModel + void InitializeSolarTable(TRestAxionSolarModel* model) { + // TOBE implemented + // This method should initialize the tables fFluxTable and fFluxLines + } + + void ExportTables(Bool_t ascii = false) override; + + void PrintMetadata() override; + + void PrintContinuumSolarTable(); + void PrintIntegratedRingFlux(); + + TRestAxionSolarHiddenPhotonFlux(); + TRestAxionSolarHiddenPhotonFlux(const char* cfgFileName, std::string name = ""); + ~TRestAxionSolarHiddenPhotonFlux(); + + ClassDefOverride(TRestAxionSolarHiddenPhotonFlux, 1); +}; +#endif diff --git a/inc/TRestAxionSolarQCDFlux.h b/inc/TRestAxionSolarQCDFlux.h index 373eab33..240fd395 100644 --- a/inc/TRestAxionSolarQCDFlux.h +++ b/inc/TRestAxionSolarQCDFlux.h @@ -91,7 +91,7 @@ class TRestAxionSolarQCDFlux : public TRestAxionSolarFlux { Double_t mass = 0) override; /// It defines how to read the solar tables at the inhereted class - Bool_t LoadTables() override; + Bool_t LoadTables(Double_t mass = 0) override; /// It returns the total integrated flux at earth in cm-2 s-1 Double_t GetTotalFlux(Double_t mass = 0) override { diff --git a/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py b/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py new file mode 100755 index 00000000..d24dfcc1 --- /dev/null +++ b/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py @@ -0,0 +1,176 @@ +#!/usr/bin/python3 + +import math +import ROOT +from ROOT import ( + TChain, + TFile, + TTree, + TCanvas, + TPad, + TRandom3, + TH1D, + TH2D, + TH3D, + TProfile, + TProfile2D, + TProfile3D, + TGraph, + TGraph2D, + TF1, + TF2, + TF3, + TFormula, + TLorentzVector, + TVector3, +) + +import argparse + +parser = argparse.ArgumentParser() +parser.add_argument("--rml", dest="rmlfile", type=str, help="The input RML file .rml") +parser.add_argument( + "--out", + dest="outfname", + type=str, + help="The output filename in png,pdf,C or any other ROOT accepted format", +) +parser.add_argument( + "--fluxname", + dest="fluxname", + type=str, + help="The name of the flux definition to be chosen from the RML", +) +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]") + +args = parser.parse_args() + +mass = 10 # eV +if args.mass != None: + mass = args.mass + +rmlfile = "fluxes.rml" +if args.rmlfile != None: + rmlfile = args.rmlfile + +outfname = "flux.png" +if args.outfname != None: + outfname = args.outfname + +fluxname = "HiddenPhoton" +if args.fluxname != None: + fluxname = args.fluxname + +samples = 20000 +if args.samples != None: + samples = args.samples + +validation = True + +ROOT.gSystem.Load("libRestFramework.so") +ROOT.gSystem.Load("libRestAxion.so") + +c1 = TCanvas("c1", "My canvas", 800, 700) +c1.GetFrame().SetBorderSize(6) +c1.GetFrame().SetBorderMode(-1) + +pad1 = TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97) +pad1.Divide(2, 2) +pad1.Draw() + +combinedFlux = ROOT.TRestAxionSolarHiddenPhotonFlux(rmlfile, fluxname) +combinedFlux.Initialize(mass) +combinedFlux.PrintMetadata() + +if combinedFlux.GetError(): + print(combinedFlux.GetErrorMessage()) + print("\nSolar flux initialization failed! Exit code : 101") + exit(101) + +comb_spt = TH2D("comb_spt", "Energy versus solar radius", 20000, 0, 20, 100, 0, 1) +for x in range(samples): + x = combinedFlux.GetRandomEnergyAndRadius((-1, -1)) + comb_spt.Fill(x[0], x[1]) + +rnd = TRandom3(0) +solarDisk = TH2D("solar_disk", "SolarDisk", 120, -1.2, 1.2, 120, -1.2, 1.2) +for x in range(samples): + angle = rnd.Rndm() * 2 * math.pi + x = combinedFlux.GetRandomEnergyAndRadius((-1, -1)) + solarDisk.Fill(x[1] * math.cos(angle), x[1] * math.sin(angle)) + +pad1.cd(1) +pad1.cd(1).SetRightMargin(0.09) +pad1.cd(1).SetLeftMargin(0.15) +pad1.cd(1).SetBottomMargin(0.15) + +comb_spt.SetStats(0) +comb_spt.GetXaxis().SetTitle("Energy [keV]") +comb_spt.GetXaxis().SetTitleSize(0.05) +comb_spt.GetXaxis().SetLabelSize(0.05) +comb_spt.GetYaxis().SetTitle("Solar radius") +comb_spt.GetYaxis().SetTitleSize(0.05) +comb_spt.GetYaxis().SetLabelSize(0.05) +comb_spt.Draw("colz") + +pad1.cd(2) +pad1.cd(2).SetLogy() +pad1.cd(2).SetRightMargin(0.09) +pad1.cd(2).SetLeftMargin(0.15) +pad1.cd(2).SetBottomMargin(0.15) +enSpt = comb_spt.ProjectionX() +enSpt.SetTitle("Energy spectrum") +enSpt.GetYaxis().SetTitleSize(0.05) +enSpt.SetStats(0) +enSpt.SetFillStyle(4050) +enSpt.SetFillColor(ROOT.kBlue - 9) +enSpt.SetLineColor(ROOT.kBlack) +enSpt.Draw() + +pad1.cd(3) +pad1.cd(3).SetRightMargin(0.09) +pad1.cd(3).SetLeftMargin(0.15) +pad1.cd(3).SetBottomMargin(0.15) +rSpt = comb_spt.ProjectionY() +rSpt.SetTitle("Radial distribution") +rSpt.GetYaxis().SetTitleSize(0.05) +rSpt.SetStats(0) +rSpt.SetFillStyle(4050) +rSpt.SetFillColor(ROOT.kBlue - 9) +rSpt.SetLineColor(ROOT.kBlack) +rSpt.Draw() + +pad1.cd(4) +pad1.cd(4).SetRightMargin(0.09) +pad1.cd(4).SetLeftMargin(0.15) +pad1.cd(4).SetBottomMargin(0.15) +solarDisk.SetStats(0) +solarDisk.GetXaxis().SetTitle("X") +solarDisk.GetXaxis().SetTitleSize(0.05) +solarDisk.GetXaxis().SetLabelSize(0.05) +solarDisk.GetYaxis().SetTitle("Y") +solarDisk.GetYaxis().SetTitleOffset(1) +solarDisk.GetYaxis().SetTitleSize(0.05) +solarDisk.GetYaxis().SetLabelSize(0.05) +solarDisk.Draw("colz") + +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) diff --git a/src/TRestAxionSolarHiddenPhotonFlux.cxx b/src/TRestAxionSolarHiddenPhotonFlux.cxx new file mode 100644 index 00000000..fe285b04 --- /dev/null +++ b/src/TRestAxionSolarHiddenPhotonFlux.cxx @@ -0,0 +1,588 @@ +/******************** 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(Double_t mass) { + if (fFluxDataFile == "" || fWidthDataFile == "" || fPlasmaFreqDataFile == "") return false; + + SetMass(mass); + + 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::LoadContinuumTable() { + if (fFluxDataFile == "") { + RESTDebug + << "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++) { + 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); + } +} + +/////////////////////////////////////////////// +/// \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 == "") { + RESTDebug << "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); + + 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++) { + 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); + } +} + +/////////////////////////////////////////////// +/// \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 == "") { + RESTDebug + << "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 << "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_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; + } + 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)); + + 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. +/// +TH1F* TRestAxionSolarHiddenPhotonFlux::GetContinuumSpectrum() { + if (fContinuumHist != nullptr) { + delete fContinuumHist; + fContinuumHist = nullptr; + } + + fContinuumHist = new TH1F("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. +/// +TH1F* TRestAxionSolarHiddenPhotonFlux::GetTotalSpectrum() { + TH1F* hc = GetContinuumSpectrum(); + + if (fTotalHist != nullptr) { + delete fTotalHist; + fTotalHist = nullptr; + } + + fTotalHist = new TH1F("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, Double_t mass) { + 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, + Double_t HPmass) { + LoadTables(HPmass); // load tables with specified hidden photon mass + std::pair result = {0, 0}; + if (!AreTablesLoaded()) 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.01; + std::pair p = {energy, radius}; + return p; + } + } + return result; +} + +/////////////////////////////////////////////// +/// \brief It prints on screen the table that has been loaded in memory +/// +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 << "-------" << 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); + + TH1F* 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 5cdd2d46..c64969df 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() { +Bool_t TRestAxionSolarQCDFlux::LoadTables(Double_t mass) { if (fFluxDataFile == "" && fFluxSptFile == "") return false; if (TRestTools::GetFileNameExtension(fFluxDataFile) == "flux") {