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/.gitmodules b/.gitmodules index 29be5f3e..7fe95128 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,6 +1,6 @@ [submodule "data"] path = data - url = https://github.com/rest-for-physics/axionlib-data.git + url = https://github.com/rest-for-physics/axionlib-data [submodule "external/solarAxionFlux"] path = external/solarAxionFlux url = https://github.com/sebhoof/SolarAxionFlux.git diff --git a/images/SolarHiddenPhotonFlux/flux0.001.png b/images/SolarHiddenPhotonFlux/flux0.001.png new file mode 100644 index 00000000..87ca31e2 Binary files /dev/null and b/images/SolarHiddenPhotonFlux/flux0.001.png differ diff --git a/images/SolarHiddenPhotonFlux/flux0.01.png b/images/SolarHiddenPhotonFlux/flux0.01.png new file mode 100644 index 00000000..d90face0 Binary files /dev/null and b/images/SolarHiddenPhotonFlux/flux0.01.png differ diff --git a/images/SolarHiddenPhotonFlux/flux0.1.png b/images/SolarHiddenPhotonFlux/flux0.1.png new file mode 100644 index 00000000..f76ea5ec Binary files /dev/null and b/images/SolarHiddenPhotonFlux/flux0.1.png differ diff --git a/images/SolarHiddenPhotonFlux/flux1.png b/images/SolarHiddenPhotonFlux/flux1.png new file mode 100644 index 00000000..60cb71b3 Binary files /dev/null and b/images/SolarHiddenPhotonFlux/flux1.png differ diff --git a/images/SolarHiddenPhotonFlux/flux10.png b/images/SolarHiddenPhotonFlux/flux10.png new file mode 100644 index 00000000..dbf8a3c3 Binary files /dev/null and b/images/SolarHiddenPhotonFlux/flux10.png differ diff --git a/images/SolarHiddenPhotonFlux/flux100.png b/images/SolarHiddenPhotonFlux/flux100.png new file mode 100644 index 00000000..02b84648 Binary files /dev/null and b/images/SolarHiddenPhotonFlux/flux100.png differ diff --git a/images/SolarHiddenPhotonFlux/flux1000.png b/images/SolarHiddenPhotonFlux/flux1000.png new file mode 100644 index 00000000..639aa3eb Binary files /dev/null and b/images/SolarHiddenPhotonFlux/flux1000.png differ diff --git a/images/SolarHiddenPhotonFlux/flux20.png b/images/SolarHiddenPhotonFlux/flux20.png new file mode 100644 index 00000000..8f8641a1 Binary files /dev/null and b/images/SolarHiddenPhotonFlux/flux20.png differ diff --git a/images/SolarHiddenPhotonFlux/flux30.png b/images/SolarHiddenPhotonFlux/flux30.png new file mode 100644 index 00000000..b8e12ead Binary files /dev/null and b/images/SolarHiddenPhotonFlux/flux30.png differ diff --git a/images/SolarHiddenPhotonFlux/flux50.png b/images/SolarHiddenPhotonFlux/flux50.png new file mode 100644 index 00000000..76c4c5b7 Binary files /dev/null and b/images/SolarHiddenPhotonFlux/flux50.png differ diff --git a/inc/TRestAxionFieldPropagationProcess.h b/inc/TRestAxionFieldPropagationProcess.h index c12e3ac8..8313bf50 100644 --- a/inc/TRestAxionFieldPropagationProcess.h +++ b/inc/TRestAxionFieldPropagationProcess.h @@ -25,8 +25,8 @@ #include "TRestAxionEvent.h" #include "TRestAxionEventProcess.h" -#include "TRestAxionField.h" #include "TRestAxionMagneticField.h" +#include "TRestAxionQCDField.h" #include "TRestPhysics.h" #include "TVector3.h" #include "TVectorD.h" @@ -44,7 +44,7 @@ class TRestAxionFieldPropagationProcess : public TRestAxionEventProcess { TRestAxionMagneticField* fMagneticField = nullptr; //! /// A pointer to TRestAxionField that implements probability calculations - TRestAxionField* fAxionField = nullptr; //! + TRestAxionQCDField* fAxionField = nullptr; //! /// A pointer to TRestBufferGas given to TRestAxionField to perform calculations in a particular gas TRestAxionBufferGas* fBufferGas = nullptr; //! diff --git a/inc/TRestAxionHiddenPhotonField.h b/inc/TRestAxionHiddenPhotonField.h new file mode 100644 index 00000000..0db3630e --- /dev/null +++ b/inc/TRestAxionHiddenPhotonField.h @@ -0,0 +1,65 @@ +/************************************************************************* + * 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 _TRestAxionHiddenPhotonField +#define _TRestAxionHiddenPhotonField + +#include "TRestAxionBufferGas.h" + +//! A basic class to define analytical axion-photon conversion calculations for axion helioscopes +class TRestAxionHiddenPhotonField : public TObject { + private: + Bool_t fDebug = false; //! + + void Initialize(); + + /// A pointer to the buffer gas definition + TRestAxionBufferGas* fBufferGas = NULL; //! + + public: + /// momentum difference q + Double_t momentumTransfer(Double_t Ea, Double_t ma, Double_t mg); + + /// vacuum conversion probability + Double_t VacuumConversion(Double_t Lcoh, Double_t Ea, Double_t ma); + + /// It enables/disables debug mode + void SetDebug(Bool_t v) { fDebug = v; } + + /// It assigns a gas buffer medium to the calculation + void AssignBufferGas(TRestAxionBufferGas* buffGas) { fBufferGas = buffGas; } + + /// It assigns a gas buffer medium to the calculation + void SetBufferGas(TRestAxionBufferGas* buffGas) { fBufferGas = buffGas; } + + Double_t GammaTransmissionProbability(Double_t Lcoh, Double_t Ea, Double_t ma, Double_t mg = 0, + Double_t absLength = 0); + + // Double_t PhotonAbsorptionProbability(Double_t Lcoh, Double_t Ea, Double_t ma, + // Double_t mg = 0, Double_t absLength = 0); + + TRestAxionHiddenPhotonField(); + ~TRestAxionHiddenPhotonField(); + + ClassDef(TRestAxionHiddenPhotonField, 2); +}; +#endif diff --git a/inc/TRestAxionLikelihood.h b/inc/TRestAxionLikelihood.h index 055bb819..a542d1c9 100644 --- a/inc/TRestAxionLikelihood.h +++ b/inc/TRestAxionLikelihood.h @@ -27,7 +27,7 @@ #include "TRandom3.h" #include "TRestAxionBufferGas.h" -#include "TRestAxionField.h" +#include "TRestAxionQCDField.h" #include "TRestAxionSolarModel.h" #include "TRestAxionSpectrum.h" @@ -62,7 +62,7 @@ class TRestAxionLikelihood : public TRestMetadata { Double_t fLastStepDensity = 0.; //-> - TRestAxionField* fAxionField; //! + TRestAxionQCDField* fAxionField; //! TRestAxionBufferGas* fBufferGas; //! TRestAxionSpectrum* fAxionSpectrum; //! diff --git a/inc/TRestAxionField.h b/inc/TRestAxionQCDField.h similarity index 95% rename from inc/TRestAxionField.h rename to inc/TRestAxionQCDField.h index 83cd950c..91f25d16 100644 --- a/inc/TRestAxionField.h +++ b/inc/TRestAxionQCDField.h @@ -20,13 +20,13 @@ * For the list of contributors see $REST_PATH/CREDITS. * *************************************************************************/ -#ifndef _TRestAxionField -#define _TRestAxionField +#ifndef _TRestAxionQCDField +#define _TRestAxionQCDField #include "TRestAxionBufferGas.h" //! A basic class to define analytical axion-photon conversion calculations for axion helioscopes -class TRestAxionField : public TObject { +class TRestAxionQCDField : public TObject { private: Bool_t fDebug = false; //! @@ -84,9 +84,9 @@ class TRestAxionField : public TObject { Double_t maMax = 0.15, Double_t rampDown = 5.0); - TRestAxionField(); - ~TRestAxionField(); + TRestAxionQCDField(); + ~TRestAxionQCDField(); - ClassDef(TRestAxionField, 2); + ClassDef(TRestAxionQCDField, 2); }; #endif diff --git a/inc/TRestAxionSolarFlux.h b/inc/TRestAxionSolarFlux.h index 73bcc00d..6c740852 100644 --- a/inc/TRestAxionSolarFlux.h +++ b/inc/TRestAxionSolarFlux.h @@ -24,8 +24,8 @@ #define _TRestAxionSolarFlux #include -#include -#include +#include +#include #include #include @@ -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; //< @@ -61,18 +64,26 @@ class TRestAxionSolarFlux : public TRestMetadata { /// It is required in order to load solar flux tables into memory void Initialize(); + /// It is required in order to load solar flux tables into memory for specific mass + void InitializeMass(Double_t mass) { + SetMass(mass); + Initialize(); + } + + /// Set mass and reinitialise + void SetMass(Double_t m) { fMass = m; } // Initialize(); } + /// It returns the integrated flux at earth in cm-2 s-1 for the given energy range - virtual Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1), Double_t mass = 0) = 0; + virtual Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1)) = 0; /// It returns the total integrated flux at earth in cm-2 s-1 - virtual Double_t GetTotalFlux(Double_t mass = 0) = 0; + virtual Double_t GetTotalFlux() = 0; /// It defines how to generate Monte Carlo energy and radius values to reproduce the flux - virtual std::pair GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1), - Double_t mass = 0) = 0; + virtual std::pair GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1)) = 0; /// It returns an energy integrated spectrum in cm-2 s-1 keV-1 - virtual TH1F* GetEnergySpectrum(Double_t m = 0) = 0; + virtual TH1D* GetEnergySpectrum() = 0; virtual TCanvas* DrawSolarFlux(); @@ -83,7 +94,9 @@ class TRestAxionSolarFlux : public TRestMetadata { Bool_t AreTablesLoaded() { return fTablesLoaded; } - TH1F* GetFluxHistogram(std::string fname, Double_t binSize = 0.01); + Double_t GetMass() { return fMass; } + + TH1D* GetFluxHistogram(std::string fname, Double_t binSize = 0.01); TCanvas* DrawFluxFile(std::string fname, Double_t binSize = 0.01); void PrintMetadata(); diff --git a/inc/TRestAxionSolarHiddenPhotonFlux.h b/inc/TRestAxionSolarHiddenPhotonFlux.h new file mode 100644 index 00000000..5a6b0398 --- /dev/null +++ b/inc/TRestAxionSolarHiddenPhotonFlux.h @@ -0,0 +1,131 @@ +/************************************************************************* + * This file is part of the REST software framework. * + * * + * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza) * + * For more information see http://gifna.unizar.es/trex * + * * + * REST is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * REST is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have a copy of the GNU General Public License along with * + * REST in $REST_PATH/LICENSE. * + * If not, see http://www.gnu.org/licenses/. * + * For the list of contributors see $REST_PATH/CREDITS. * + *************************************************************************/ + +#ifndef _TRestAxionSolarHiddenPhotonFlux +#define _TRestAxionSolarHiddenPhotonFlux + +#include +#include + +//! A metadata class to load tabulated solar hidden photon fluxes. Kinetic mixing set to 1. +class TRestAxionSolarHiddenPhotonFlux : public TRestAxionSolarFlux { + private: + /// The filename containing the solar flux table with continuum spectrum + std::string fFluxDataFile = ""; //< + + /// The filename containing the resonance width (wGamma) + std::string fWidthDataFile = ""; //< + + /// The filename containing the plasma frequency (wp) table + std::string fPlasmaFreqDataFile = ""; //< + + /// It will be used when loading `.flux` files to define the input file energy binsize in eV. + Double_t fBinSize = 0; //< + + /// The tabulated solar flux continuum spectra TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fFluxTable; //! + + /// The tabulated solar flux continuum spectra TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fContinuumTable; //! + + /// The tabulated resonance width TH1D(200,0,20)keV in eV2 versus solar radius + std::vector fWidthTable; //! + + /// The solar plasma frequency vector in eV versus solar radius + std::vector fPlasmaFreqTable; //! + + /// The total solar flux TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fFullFluxTable; //! + + /// Accumulative integrated solar flux for each solar ring for continuum spectrum (renormalized to unity) + std::vector fFluxTableIntegrals; //! + + /// Total solar flux for monochromatic contributions + Double_t fTotalContinuumFlux = 0; //! + + /// A pointer to the continuum spectrum histogram + TH1D* fContinuumHist = nullptr; //! + + /// A pointer to the superposed monochromatic and continuum spectrum histogram + TH1D* fTotalHist = nullptr; //! + + void ReadFluxFile(); + void LoadContinuumFluxTable(); + void LoadMonoChromaticFluxTable(); + void IntegrateSolarFluxes(); + + public: + /// It returns true if continuum flux spectra was loaded + Bool_t isSolarTableLoaded() { return fFluxTable.size() > 0; } + + /// It returns true if width table was loaded + Bool_t isWidthTableLoaded() { return fWidthTable.size() > 0; } + + /// It returns true if plasma frequency table was loaded + Bool_t isPlasmaFreqLoaded() { return fPlasmaFreqTable.size() > 0; } + + /// It returns the integrated flux at earth in cm-2 s-1 for the given energy range + Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1)) override; + + /// It defines how to generate Monte Carlo energy and radius values to reproduce the flux + std::pair GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1)) override; + + /// It defines how to read the solar tables at the inhereted class for a given mass in eV + Bool_t LoadTables() override; + + void LoadContinuumTable(); + void LoadWidthTable(); + void LoadPlasmaFreqTable(); + + // calculate solar HP flux from the 3 tables and mass + void CalculateSolarFlux(); + /// It returns the total integrated flux at earth in cm-2 s-1 + Double_t GetTotalFlux() override { return fTotalContinuumFlux; } + + /// It returns an energy integrated spectrum in cm-2 s-1 keV-1 + TH1D* GetEnergySpectrum() override { return GetTotalSpectrum(); } + + TH1D* GetContinuumSpectrum(); + TH1D* GetTotalSpectrum(); + + virtual TCanvas* DrawSolarFlux() override; + + /// Tables might be loaded using a solar model description by TRestAxionSolarModel + void InitializeSolarTable(TRestAxionSolarModel* model) { + // TOBE implemented + // This method should initialize the tables fFluxTable and fFluxLines + } + + void ExportTables(Bool_t ascii = false) override; + + void PrintMetadata() override; + + void PrintContinuumSolarTable(); + void PrintIntegratedRingFlux(); + + TRestAxionSolarHiddenPhotonFlux(); + TRestAxionSolarHiddenPhotonFlux(const char* cfgFileName, std::string name = ""); + ~TRestAxionSolarHiddenPhotonFlux(); + + ClassDefOverride(TRestAxionSolarHiddenPhotonFlux, 1); +}; +#endif diff --git a/inc/TRestAxionSolarQCDFlux.h b/inc/TRestAxionSolarQCDFlux.h index 373eab33..fb99ebb0 100644 --- a/inc/TRestAxionSolarQCDFlux.h +++ b/inc/TRestAxionSolarQCDFlux.h @@ -41,11 +41,11 @@ class TRestAxionSolarQCDFlux : public TRestAxionSolarFlux { /// It will be used when loading `.flux` files to define the threshold for peak identification Double_t fPeakSigma = 0; //< - /// The tabulated solar flux continuum spectra TH1F(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius - std::vector fFluxTable; //! + /// The tabulated solar flux continuum spectra TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius + std::vector fFluxTable; //! /// The tabulated solar flux in cm-2 s-1 for a number of monochromatic energies versus solar radius - std::map fFluxLines; //! + std::map fFluxLines; //! /// Accumulative integrated solar flux for each solar ring for continuum spectrum (renormalized to unity) std::vector fFluxTableIntegrals; //! @@ -63,13 +63,13 @@ class TRestAxionSolarQCDFlux : public TRestAxionSolarFlux { Double_t fFluxRatio = 0; //! /// A pointer to the continuum spectrum histogram - TH1F* fContinuumHist = nullptr; //! + TH1D* fContinuumHist = nullptr; //! /// A pointer to the monochromatic spectrum histogram - TH1F* fMonoHist = nullptr; //! + TH1D* fMonoHist = nullptr; //! /// A pointer to the superposed monochromatic and continuum spectrum histogram - TH1F* fTotalHist = nullptr; //! + TH1D* fTotalHist = nullptr; //! void ReadFluxFile(); void LoadContinuumFluxTable(); @@ -84,26 +84,23 @@ class TRestAxionSolarQCDFlux : public TRestAxionSolarFlux { Bool_t isSolarSpectrumLoaded() { return fFluxLines.size() > 0; } /// It returns the integrated flux at earth in cm-2 s-1 for the given energy range - Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1), Double_t mass = 0) override; + Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1)) override; /// It defines how to generate Monte Carlo energy and radius values to reproduce the flux - std::pair GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1), - Double_t mass = 0) override; + std::pair GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1)) override; /// It defines how to read the solar tables at the inhereted class Bool_t LoadTables() override; /// It returns the total integrated flux at earth in cm-2 s-1 - Double_t GetTotalFlux(Double_t mass = 0) override { - return fTotalContinuumFlux + fTotalMonochromaticFlux; - } + Double_t GetTotalFlux() override { return fTotalContinuumFlux + fTotalMonochromaticFlux; } /// It returns an energy integrated spectrum in cm-2 s-1 keV-1 - TH1F* GetEnergySpectrum(Double_t m = 0) override { return GetTotalSpectrum(); } + TH1D* GetEnergySpectrum() override { return GetTotalSpectrum(); } - TH1F* GetContinuumSpectrum(); - TH1F* GetMonochromaticSpectrum(); - TH1F* GetTotalSpectrum(); + TH1D* GetContinuumSpectrum(); + TH1D* GetMonochromaticSpectrum(); + TH1D* GetTotalSpectrum(); virtual TCanvas* DrawSolarFlux() override; diff --git a/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py b/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py new file mode 100755 index 00000000..5b83b758 --- /dev/null +++ b/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py @@ -0,0 +1,165 @@ +#!/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("--m", dest="mass", type=float, help="Hidden photon mass [eV]") +args = parser.parse_args() + +mass = 10.0 # 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 = False + +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.SetMass(mass) +print(combinedFlux.GetMass()) +combinedFlux.Initialize() +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 i in range(samples): + x = combinedFlux.GetRandomEnergyAndRadius((-1, -1)) + # print(x) + 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) + +# exit(0) diff --git a/src/TRestAxionFieldPropagationProcess.cxx b/src/TRestAxionFieldPropagationProcess.cxx index e7de98b3..03021654 100644 --- a/src/TRestAxionFieldPropagationProcess.cxx +++ b/src/TRestAxionFieldPropagationProcess.cxx @@ -171,7 +171,7 @@ void TRestAxionFieldPropagationProcess::InitProcess() { } if (!fAxionField) { - fAxionField = new TRestAxionField(); + fAxionField = new TRestAxionQCDField(); fBufferGas = (TRestAxionBufferGas*)this->GetMetadata("TRestAxionBufferGas"); if (fBufferGas) diff --git a/src/TRestAxionHiddenPhotonField.cxx b/src/TRestAxionHiddenPhotonField.cxx new file mode 100644 index 00000000..7c0d9560 --- /dev/null +++ b/src/TRestAxionHiddenPhotonField.cxx @@ -0,0 +1,325 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// TRestAxionHiddenPhotonField is a class used to calculate the axion-photon mixing +/// and determine the probability of the particle being in an axion or photon +/// state. +/// +/// A peculiarity from this class is that it encapsulates internally the high +/// precision calculations using the real precisions types using TRestComplex. +/// It is known that double precision is not good enough in some scenarios. +/// +///-------------------------------------------------------------------------- +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2019-March: First concept and implementation of TRestAxionHiddenPhotonField class. +/// Javier Galan +/// 2023-September: Separation of QCD and HiddenPhoton classes. +/// Tom O'Shea +/// +/// \class TRestAxionHiddenPhotonField +/// \author Javier Galan +/// +///
+/// +#include "TRestAxionHiddenPhotonField.h" + +#include + +#include "TH1F.h" + +#ifdef USE_MPFR +#include "TRestComplex.h" +#endif + +#include + +using namespace std; + +ClassImp(TRestAxionHiddenPhotonField); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestAxionHiddenPhotonField::TRestAxionHiddenPhotonField() { Initialize(); } + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestAxionHiddenPhotonField::~TRestAxionHiddenPhotonField() {} + +/////////////////////////////////////////////// +/// \brief Initialization of TRestAxionHiddenPhotonField class +/// +/// It sets the default real precision to be used with mpfr types. Now it is 30 digits. +/// So that we can still calculate numbers such as : 1.0 - 1.e-30 +/// +void TRestAxionHiddenPhotonField::Initialize() { +#ifdef USE_MPFR + TRestComplex::SetPrecision(30); +#endif + + fBufferGas = NULL; + + /// MOVED TO TRestAxionHiddenPhotonFieldPropagationProcess class + /// faxion = SetComplexReal(1, 0); + /// fAem = SetComplexReal(0, 0); +} + +/////////////////////////////////////////////// +/// momentum difference q in eV, useful in all calculations +/// Ea in keV, ma in eV, mg in eV +Double_t TRestAxionHiddenPhotonField::momentumTransfer(Double_t Ea, Double_t ma, Double_t mg) { + Ea *= 1000; // [eV] Ea is given in keV + return sqrt(Ea * Ea - mg * mg) - sqrt(Ea * Ea - ma * ma); //[eV] +} + +/////////////////////////////////////////////// +/// vacuum conversion probability +/// ie. when photonMass = 0 +/// Ea in keV, ma in eV, Lcoh in mm +Double_t TRestAxionHiddenPhotonField::VacuumConversion(Double_t Lcoh, Double_t Ea, Double_t ma) { + Lcoh *= REST_Physics::PhMeterIneV / 1000; // [eV-1] + Double_t q = momentumTransfer(Ea, ma, 0.); // [eV] + return pow(2 * sin(q * Lcoh / 2), 2); +} + +/////////////////////////////////////////////// +/// \brief Performs the calculation of hidden photon - photon conversion probability. +/// +/// If m_gamma (mg) is not given as an argument, i.e. it is equal to zero, then m_gamma +/// will be obtainned from the buffer gas definition. If no buffer gas has been assigned +/// then the medium will be assumed to be vacuum. +/// +/// Ea in keV, ma in eV, Lcoh in mm, mg in eV, absLength in cm-1 +/// +/// The returned value is given for chi = 1 +/// +Double_t TRestAxionHiddenPhotonField::GammaTransmissionProbability(Double_t Lcoh, Double_t Ea, Double_t ma, + Double_t mg, Double_t absLength) { +#ifndef USE_MPFR + RESTWarning + << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFR=ON to your REST compilation" + << RESTendl; + RESTWarning << "TRestAxionHiddenPhotonField::GammaTransmissionProbability will return 0" << RESTendl; + return 0; +#else + mpfr::mpreal hiddenPhotonMass = ma; // [eV] + mpfr::mpreal cohLength = Lcoh / 1000.; // [m] Default REST units are mm + + mpfr::mpreal photonMass = mg; // [eV] + + if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(Ea); + + RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; + RESTDebug << " TRestAxionHiddenPhotonField::GammaTransmissionProbability. Parameter summary" << RESTendl; + RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; + RESTDebug << " Hidden photon mass : " << ma << " eV" << RESTendl; + RESTDebug << " Hidden photon energy : " << Ea << " keV" << RESTendl; + RESTDebug << " Lcoh : " << Lcoh << " mm" << RESTendl; + RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; + + if (photonMass == 0.0) return VacuumConversion(Lcoh, Ea, ma); + + mpfr::mpreal q = momentumTransfer(Ea, ma, mg); // eV + mpfr::mpreal l = cohLength * REST_Physics::PhMeterIneV; // eV-1 + + mpfr::mpreal Gamma = absLength; // cm-1 + if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 + mpfr::mpreal GammaL = Gamma * cohLength * 100; + Gamma *= 100 / REST_Physics::PhMeterIneV; // eV + + if (fDebug) { + RESTDebug << "+------------------------+" << RESTendl; + RESTDebug << " Intermediate calculations" << RESTendl; + RESTDebug << " q : " << q << " eV" << RESTendl; + RESTDebug << " l : " << l << " eV-1" << RESTendl; + RESTDebug << " ql : " << q * l << RESTendl; + RESTDebug << "Gamma : " << Gamma << " eV" << RESTendl; + RESTDebug << "GammaL : " << GammaL << RESTendl; + RESTDebug << "+------------------------+" << RESTendl; + } + + mpfr::mpreal MFactor = + pow(ma, 4) / (pow(ma * ma - mg * mg, 2) + pow(Ea * Gamma, 2)); // all should be in eV + + if (fDebug) { + RESTDebug << "Mfactor : " << MFactor << RESTendl; + RESTDebug << "cos(ql) : " << cos(q * l) << RESTendl; + RESTDebug << "Exp(-GammaL) : " << exp(-GammaL) << RESTendl; + } + + double sol = (double)(MFactor * (1 + exp(-GammaL) - 2 * exp(-GammaL / 2) * cos(q * l))); + + RESTDebug << "Axion-photon transmission probability : " << sol << RESTendl; + + return sol; +#endif +} + +/* +/////////////////////////////////////////////// +/// \brief Performs the calculation of the photon absorbtion probability in the buffer gas. +/// +/// If not provided as argument, m_g it will be attempted to obtain the value from the buffer gas +/// definition. If no buffer gas has been assigned the medium will be assumed to be vacuum. +/// +/// Ea in keV, ma in eV, Lcoh in mm, mg in eV, absLength in cm-1 +/// +/// The returned value is given for chi = 1 +/// +Double_t TRestAxionHiddenPhotonField::AxionAbsorptionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, +Double_t ma, Double_t mg, Double_t absLength) { #ifndef USE_MPFR RESTWarning + << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFr=ON to your REST compilation" + << RESTendl; + RESTWarning << "TRestAxionHiddenPhotonField::GammaTransmissionProbability will return 0" << RESTendl; + return 0; +#else + mpfr::mpreal axionMass = ma; + mpfr::mpreal cohLength = Lcoh / 1000.; // Default REST units are mm; + + mpfr::mpreal photonMass = mg; + if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(Ea); + + if (fDebug) { + RESTDebug << "+--------------------------------------------------------------------------+" + << RESTendl; + RESTDebug << " TRestAxionHiddenPhotonField::GammaTransmissionProbability. Parameter summary" << +RESTendl; RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; RESTDebug << " Axion mass : " << +ma << " eV" << RESTendl; RESTDebug << " Axion energy : " << Ea << " keV" << RESTendl; RESTDebug << " Lcoh : " +<< Lcoh << " mm" << RESTendl; RESTDebug << " Bmag : " << Bmag << " T" << RESTendl; RESTDebug << +"+--------------------------------------------------------------------------+" + << RESTendl; + } + + if (ma == 0.0 && photonMass == 0.0) return BLHalfSquared(Bmag, Lcoh); + + mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / Ea / 1000.0; + mpfr::mpreal l = cohLength * REST_Physics::PhMeterIneV; + mpfr::mpreal phi = q * l; + + mpfr::mpreal Gamma = absLength; + if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 + mpfr::mpreal GammaL = Gamma * cohLength * 100; + + if (fDebug) { + RESTDebug << "+------------------------+" << RESTendl; + RESTDebug << " Intermediate calculations" << RESTendl; + RESTDebug << " q : " << q << " eV" << RESTendl; + RESTDebug << " l : " << l << " eV-1" << RESTendl; + RESTDebug << " phi : " << phi << RESTendl; + RESTDebug << "Gamma : " << Gamma << RESTendl; + RESTDebug << "GammaL : " << GammaL << RESTendl; + RESTDebug << "+------------------------+" << RESTendl; + } + + mpfr::mpreal MFactor = phi * phi + GammaL * GammaL / 4.0; + MFactor = 1.0 / MFactor; + + if (fDebug) { + RESTDebug << "Mfactor : " << MFactor << RESTendl; + RESTDebug << "(BL/2)^2 : " << BLHalfSquared(Bmag, Lcoh) << RESTendl; + RESTDebug << "cos(phi) : " << cos(phi) << RESTendl; + RESTDebug << "Exp(-GammaL) : " << exp(-GammaL) << RESTendl; + } + + double sol = (double)(MFactor * BLHalfSquared(Bmag, Lcoh) * GammaL); + + if (fDebug) RESTDebug << "Axion-photon absorption probability : " << sol << RESTendl; + + return sol; +#endif +} +*/ + +/// Commented because it uses ComplexReal structure that is moved to +/// TRestAxionHiddenPhotonFieldPropagationProcess class +/* +void TRestAxionHiddenPhotonField::PropagateAxion(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, + Double_t mg, Double_t absLength) { + mpfr::mpreal axionMass = ma; + mpfr::mpreal cohLength = Lcoh / 1000.; // Default REST units are mm; + + mpfr::mpreal photonMass = mg; + if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(Ea); + + if (fDebug) { + RESTDebug << "+--------------------------------------------------------------------------+" << +RESTendl; RESTDebug << " TRestAxionHiddenPhotonField::GammaTransmissionProbability. Parameter summary" << +RESTendl; RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; RESTDebug << " Axion mass : " << +ma << " eV" << RESTendl; RESTDebug << " Axion energy : " << Ea << " keV" << RESTendl; RESTDebug << " Lcoh : " +<< Lcoh << " mm" << RESTendl; RESTDebug << " Bmag : " << Bmag << " T" << RESTendl; RESTDebug << +"+--------------------------------------------------------------------------+" << RESTendl; + } + + mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / Ea / 1000.0; + mpfr::mpreal l = cohLength * PhMeterIneV; + mpfr::mpreal phi = q * l; + + mpfr::mpreal Gamma = absLength; + if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 + mpfr::mpreal GammaL = Gamma * cohLength * 100; + + if (fDebug) { + RESTDebug << "+------------------------+" << RESTendl; + RESTDebug << " Intermediate calculations" << RESTendl; + RESTDebug << " q : " << q << " eV" << RESTendl; + RESTDebug << " l : " << l << " eV-1" << RESTendl; + RESTDebug << " phi : " << phi << RESTendl; + RESTDebug << "Gamma : " << Gamma << RESTendl; + RESTDebug << "GammaL : " << GammaL << RESTendl; + RESTDebug << "+------------------------+" << RESTendl; + } + + mpfr::mpreal bl = BL(Bmag, Lcoh); + + /// We have now calculated the main quantities BL, QL, and GammaL + + ComplexReal I = SetComplexReal(0, 1); + ComplexReal ExpPhi = SetComplexReal(cos(-phi), sin(-phi)); + + ComplexReal Bcomplex = SetComplexReal(BL(Bmag, Lcoh), 0); + ComplexReal Qcomplex = SetComplexReal(phi, -GammaL / 2); + + ComplexReal Bterm = ComplexCocient(Bcomplex, Qcomplex); + Bterm = ComplexProduct(I, Bterm); + + mpfr::mpreal ExpGamma = exp(-GammaL / 2.); + Double_t ExpGammaDouble = TMath::Exp((Double_t)-GammaL / 2.); + + cout.precision(30); + + if (fDebug) { + cout << "ExpGamma : " << ExpGamma << RESTendl; + cout << "ExpGammaDouble : " << ExpGammaDouble << RESTendl; + RESTDebug << "(BL/2)^2 : " << BLHalfSquared(Bmag, Lcoh) << RESTendl; + RESTDebug << "cos(phi) : " << cos(phi) << RESTendl; + RESTDebug << "Exp(-GammaL) : " << exp(-GammaL) << RESTendl; + } + + // if (fDebug) RESTDebug << "Axion-photon absorption probability : " << sol << RESTendl; +} +*/ diff --git a/src/TRestAxionLikelihood.cxx b/src/TRestAxionLikelihood.cxx index f5d26c17..4f9b9e74 100644 --- a/src/TRestAxionLikelihood.cxx +++ b/src/TRestAxionLikelihood.cxx @@ -84,7 +84,7 @@ void TRestAxionLikelihood::Initialize() { fBufferGas = new TRestAxionBufferGas(); // Conversion probabilities as defined by van Bibber paper - fAxionField = new TRestAxionField(); + fAxionField = new TRestAxionQCDField(); fAxionField->AssignBufferGas(fBufferGas); diff --git a/src/TRestAxionField.cxx b/src/TRestAxionQCDField.cxx similarity index 88% rename from src/TRestAxionField.cxx rename to src/TRestAxionQCDField.cxx index fae12e99..6a5c5f3a 100644 --- a/src/TRestAxionField.cxx +++ b/src/TRestAxionQCDField.cxx @@ -21,7 +21,7 @@ *************************************************************************/ ////////////////////////////////////////////////////////////////////////// -/// TRestAxionField is a class used to calculate the axion-photon mixing +/// TRestAxionQCDField is a class used to calculate the axion-photon mixing /// and determine the probability of the particle being in an axion or photon /// state. /// @@ -35,15 +35,17 @@ /// /// History of developments: /// -/// 2019-March: First concept and implementation of TRestAxionField class. +/// 2019-March: First concept and implementation of TRestAxionQCDField class. /// Javier Galan +/// 2023-September: Separation of QCD and HiddenPhoton classes. +/// Tom O'Shea /// -/// \class TRestAxionField +/// \class TRestAxionQCDField /// \author Javier Galan /// ///
/// -#include "TRestAxionField.h" +#include "TRestAxionQCDField.h" #include @@ -57,32 +59,32 @@ using namespace std; -ClassImp(TRestAxionField); +ClassImp(TRestAxionQCDField); /////////////////////////////////////////////// /// \brief Default constructor /// -TRestAxionField::TRestAxionField() { Initialize(); } +TRestAxionQCDField::TRestAxionQCDField() { Initialize(); } /////////////////////////////////////////////// /// \brief Default destructor /// -TRestAxionField::~TRestAxionField() {} +TRestAxionQCDField::~TRestAxionQCDField() {} /////////////////////////////////////////////// -/// \brief Initialization of TRestAxionField class +/// \brief Initialization of TRestAxionQCDField class /// /// It sets the default real precision to be used with mpfr types. Now it is 30 digits. /// So that we can still calculate numbers such as : 1.0 - 1.e-30 /// -void TRestAxionField::Initialize() { +void TRestAxionQCDField::Initialize() { #ifdef USE_MPFR TRestComplex::SetPrecision(30); #endif fBufferGas = NULL; - /// MOVED TO TRestAxionFieldPropagationProcess class + /// MOVED TO TRestAxionQCDFieldPropagationProcess class /// faxion = SetComplexReal(1, 0); /// fAem = SetComplexReal(0, 0); } @@ -93,7 +95,7 @@ void TRestAxionField::Initialize() { /// `Lcoh` should be expressed in `mm`, and `Bmag` in `T`. /// The result will be given for an axion-photon coupling of 10^{-10} GeV^{-1} /// -double TRestAxionField::BL(Double_t Bmag, Double_t Lcoh) { +double TRestAxionQCDField::BL(Double_t Bmag, Double_t Lcoh) { Double_t lengthInMeters = Lcoh / 1000.; Double_t tm = REST_Physics::lightSpeed / REST_Physics::naturalElectron * 1.0e-9; // GeV @@ -109,7 +111,7 @@ double TRestAxionField::BL(Double_t Bmag, Double_t Lcoh) { /// `Lcoh` should be expressed in `mm`, and `Bmag` in `T`. /// The result will be given for an axion-photon coupling of 10^{-10} GeV^{-1} /// -double TRestAxionField::BLHalfSquared(Double_t Bmag, Double_t Lcoh) // (BL/2)**2 +double TRestAxionQCDField::BLHalfSquared(Double_t Bmag, Double_t Lcoh) // (BL/2)**2 { Double_t lengthInMeters = Lcoh / 1000.; @@ -134,12 +136,12 @@ double TRestAxionField::BLHalfSquared(Double_t Bmag, Double_t Lcoh) // (BL/2)** /// /// The returned value is given for g_ag = 10^-10 GeV-1 /// -Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, Double_t absLength) { +Double_t TRestAxionQCDField::GammaTransmissionProbability(Double_t ma, Double_t mg, Double_t absLength) { #ifndef USE_MPFR RESTWarning << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFR=ON to your REST compilation" << RESTendl; - RESTWarning << "TRestAxionField::GammaTransmissionProbability will return 0" << RESTendl; + RESTWarning << "TRestAxionQCDField::GammaTransmissionProbability will return 0" << RESTendl; return 0; #else mpfr::mpreal axionMass = ma; @@ -150,7 +152,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(fEa); RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << RESTendl; + RESTDebug << " TRestAxionQCDField::GammaTransmissionProbability. Parameter summary" << RESTendl; RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; RESTDebug << " Axion mass : " << ma << " eV" << RESTendl; RESTDebug << " Axion energy : " << fEa << " keV" << RESTendl; @@ -202,8 +204,8 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, /// \brief On top of calculating the gamma transmission probability it will assign new values /// for the magnetic field (Bmag/T), coherence length (Lcoh/mm) and axion energy (Ea/keV). /// -Double_t TRestAxionField::GammaTransmissionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, - Double_t mg, Double_t absLength) { +Double_t TRestAxionQCDField::GammaTransmissionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, + Double_t ma, Double_t mg, Double_t absLength) { fBmag = Bmag; fLcoh = Lcoh; fEa = Ea; @@ -230,14 +232,14 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t Bmag, Double_t L /// to solve the problem with a density profile. TOBE implemented in a new method if needed, where /// Gamma is not constant and \integral{q(z)} is integrated at each step. /// -Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bmag, Double_t deltaL, - Double_t Ea, Double_t ma, Double_t mg, - Double_t absLength) { +Double_t TRestAxionQCDField::GammaTransmissionProbability(std::vector Bmag, Double_t deltaL, + Double_t Ea, Double_t ma, Double_t mg, + Double_t absLength) { #ifndef USE_MPFR RESTWarning << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFR=ON to your REST compilation" << RESTendl; - RESTWarning << "TRestAxionField::GammaTransmissionProbability will return 0" << RESTendl; + RESTWarning << "TRestAxionQCDField::GammaTransmissionProbability will return 0" << RESTendl; return 0; #else mpfr::mpreal axionMass = ma; @@ -254,7 +256,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma if (Bmag.size() > 0) fieldAverage = std::accumulate(Bmag.begin(), Bmag.end(), 0.0) / Bmag.size(); RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << RESTendl; + RESTDebug << " TRestAxionQCDField::GammaTransmissionProbability. Parameter summary" << RESTendl; RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; RESTDebug << " Axion mass : " << ma << " eV" << RESTendl; RESTDebug << " Axion energy : " << Ea << " keV" << RESTendl; @@ -342,12 +344,12 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma /// /// The returned value is given for g_ag = 10^-10 GeV-1 /// -Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, Double_t absLength) { +Double_t TRestAxionQCDField::AxionAbsorptionProbability(Double_t ma, Double_t mg, Double_t absLength) { #ifndef USE_MPFR RESTWarning - << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFr=ON to your REST compilation" + << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFR=ON to your REST compilation" << RESTendl; - RESTWarning << "TRestAxionField::GammaTransmissionProbability will return 0" << RESTendl; + RESTWarning << "TRestAxionQCDField::AxionAbsorptionProbability will return 0" << RESTendl; return 0; #else mpfr::mpreal axionMass = ma; @@ -359,7 +361,7 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, D if (fDebug) { RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << RESTendl; + RESTDebug << " TRestAxionQCDField::GammaTransmissionProbability. Parameter summary" << RESTendl; RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; RESTDebug << " Axion mass : " << ma << " eV" << RESTendl; RESTDebug << " Axion energy : " << fEa << " keV" << RESTendl; @@ -412,8 +414,8 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, D /// \brief On top of calculating the axion absorption probability it will assign new values /// for the magnetic field (Bmag/T), coherence length (Lcoh/mm) and axion energy (Ea/keV). /// -Double_t TRestAxionField::AxionAbsorptionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, - Double_t mg, Double_t absLength) { +Double_t TRestAxionQCDField::AxionAbsorptionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, + Double_t ma, Double_t mg, Double_t absLength) { fBmag = Bmag; fLcoh = Lcoh; fEa = Ea; @@ -433,7 +435,7 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t Bmag, Double_t Lco /// IMPORTANT: In the case that the buffer gas is not defined, this method will return the mass at which the /// probability reaches half of the maximum **vacuum** probability. /// -Double_t TRestAxionField::GammaTransmissionFWHM(Double_t step) { +Double_t TRestAxionQCDField::GammaTransmissionFWHM(Double_t step) { Double_t maxMass = 10; // 10eV is the maximum mass (exit condition) Double_t resonanceMass = 0; @@ -481,9 +483,9 @@ Double_t TRestAxionField::GammaTransmissionFWHM(Double_t step) { /// /// For additional info see PR: https://github.com/rest-for-physics/axionlib/pull/78 /// -std::vector> TRestAxionField::GetMassDensityScanning(std::string gasName, - double maMax, - double rampDown) { +std::vector> TRestAxionQCDField::GetMassDensityScanning(std::string gasName, + double maMax, + double rampDown) { std::vector> massDensityPairs; // Storing the gas pointer, if there was one diff --git a/src/TRestAxionSolarFlux.cxx b/src/TRestAxionSolarFlux.cxx index a2e0ccd6..62549922 100644 --- a/src/TRestAxionSolarFlux.cxx +++ b/src/TRestAxionSolarFlux.cxx @@ -64,7 +64,7 @@ /// on the flux distribution initialized. /// - **LoadTables**: It is called by TRestAxionSolarFlux::Initialize to allow the inherited /// class to load all the necessary tables in memory. -/// - **GetEnergySpectrum**: It should return a TH1F pointer with a energy spectrum histogram. +/// - **GetEnergySpectrum**: It should return a TH1D pointer with a energy spectrum histogram. /// ///-------------------------------------------------------------------------- /// @@ -125,12 +125,7 @@ void TRestAxionSolarFlux::Initialize() { fTablesLoaded = false; if (LoadTables()) fTablesLoaded = true; - if (!fRandom) { - delete fRandom; - fRandom = nullptr; - } - - if (fRandom != nullptr) { + if (fRandom) { delete fRandom; fRandom = nullptr; } @@ -139,11 +134,17 @@ void TRestAxionSolarFlux::Initialize() { fSeed = fRandom->TRandom::GetSeed(); } +/////////////////////////////////////////////// +/// \brief Initialization of TRestAxionSolarFlux members with specific mass +/// +// void TRestAxionSolarFlux::InitializeMass( Double_t mass ) { SetMass(mass); RESTMetadata << GetMass() << +// RESTendl; } // SetMass calls Initialize + /////////////////////////////////////////////// /// \brief It builds a histogram using the contents of the .flux file given /// in the argument. /// -TH1F* TRestAxionSolarFlux::GetFluxHistogram(string fname, Double_t binSize) { +TH1D* TRestAxionSolarFlux::GetFluxHistogram(string fname, Double_t binSize) { string fullPathName = SearchFile(fname); std::vector> fluxData; @@ -160,7 +161,7 @@ TH1F* TRestAxionSolarFlux::GetFluxHistogram(string fname, Double_t binSize) { originalHist->Fill(r, en, flux); } - return (TH1F*)originalHist->ProjectionY(); + return (TH1D*)originalHist->ProjectionY(); } /////////////////////////////////////////////// @@ -206,7 +207,7 @@ TCanvas* TRestAxionSolarFlux::DrawSolarFlux() { pad1->SetLeftMargin(0.15); pad1->SetBottomMargin(0.15); - TH1F* ht = GetEnergySpectrum(); + TH1D* ht = GetEnergySpectrum(); ht->SetLineColor(kBlack); ht->SetFillStyle(4050); ht->SetFillColor(kBlue - 10); diff --git a/src/TRestAxionSolarHiddenPhotonFlux.cxx b/src/TRestAxionSolarHiddenPhotonFlux.cxx new file mode 100644 index 00000000..66ffa0ec --- /dev/null +++ b/src/TRestAxionSolarHiddenPhotonFlux.cxx @@ -0,0 +1,629 @@ +/******************** 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 = 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 ) + + // m4 * wG * flux / ( (m2 - wp2)^2 + (w G)^2) + // 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. +/// +TH1D* TRestAxionSolarHiddenPhotonFlux::GetContinuumSpectrum() { + if (fContinuumHist != nullptr) { + delete fContinuumHist; + fContinuumHist = nullptr; + } + + fContinuumHist = new TH1D("ContinuumHist", "", 200, 0, 20); + for (const auto& x : fFluxTable) { + fContinuumHist->Add(x); + } + + fContinuumHist->SetStats(0); + fContinuumHist->GetXaxis()->SetTitle("Energy [keV]"); + fContinuumHist->GetXaxis()->SetTitleSize(0.05); + fContinuumHist->GetXaxis()->SetLabelSize(0.05); + fContinuumHist->GetYaxis()->SetTitle("Flux [cm-2 s-1 keV-1]"); + fContinuumHist->GetYaxis()->SetTitleSize(0.05); + fContinuumHist->GetYaxis()->SetLabelSize(0.05); + + return fContinuumHist; +} + +/////////////////////////////////////////////// +/// \brief Same as GetContinuumSpectrum, the flux will be +/// expressed in cm-2 s-1 keV-1. Binned in 1eV steps. +/// +TH1D* TRestAxionSolarHiddenPhotonFlux::GetTotalSpectrum() { + TH1D* hc = GetContinuumSpectrum(); + + if (fTotalHist != nullptr) { + delete fTotalHist; + fTotalHist = nullptr; + } + + fTotalHist = new TH1D("fTotalHist", "", 20000, 0, 20); + for (int n = 0; n < hc->GetNbinsX(); n++) { + for (int m = 0; m < 100; m++) { + fTotalHist->SetBinContent(n * 100 + 1 + m, hc->GetBinContent(n + 1)); + } + } + + fTotalHist->SetStats(0); + fTotalHist->GetXaxis()->SetTitle("Energy [keV]"); + fTotalHist->GetXaxis()->SetTitleSize(0.05); + fTotalHist->GetXaxis()->SetLabelSize(0.05); + fTotalHist->GetYaxis()->SetTitle("Flux [cm-2 s-1 keV-1]"); + fTotalHist->GetYaxis()->SetTitleSize(0.05); + fTotalHist->GetYaxis()->SetLabelSize(0.05); + + return fTotalHist; +} + +/////////////////////////////////////////////// +/// \brief A helper method to initialize the internal class data members with the +/// integrated flux for each solar ring. It will be called by TRestAxionSolarHiddenPhotonFlux::Initialize. +/// +void TRestAxionSolarHiddenPhotonFlux::IntegrateSolarFluxes() { + fTotalContinuumFlux = 0.0; + for (unsigned int n = 0; n < fFluxTable.size(); n++) { + fTotalContinuumFlux += fFluxTable[n]->Integral() * 0.1; // We integrate in 100eV steps + fFluxTableIntegrals.push_back(fTotalContinuumFlux); + } + + for (unsigned int n = 0; n < fFluxTableIntegrals.size(); n++) + fFluxTableIntegrals[n] /= fTotalContinuumFlux; +} + +/////////////////////////////////////////////// +/// \brief It returns the integrated flux at earth in cm-2 s-1 for the given energy range +/// +Double_t TRestAxionSolarHiddenPhotonFlux::IntegrateFluxInRange(TVector2 eRange) { + if (eRange.X() == -1 && eRange.Y() == -1) { + if (GetTotalFlux() == 0) IntegrateSolarFluxes(); + return GetTotalFlux(); + } + + Double_t flux = 0; + fTotalContinuumFlux = 0.0; + for (unsigned int n = 0; n < fFluxTable.size(); n++) { + flux += fFluxTable[n]->Integral(fFluxTable[n]->FindFixBin(eRange.X()), + fFluxTable[n]->FindFixBin(eRange.Y())) * + 0.1; // We integrate in 100eV steps + } + + return flux; +} + +/////////////////////////////////////////////// +/// \brief It returns a random solar radius position and energy according to the +/// flux distributions defined inside the solar tables loaded in the class +/// +std::pair TRestAxionSolarHiddenPhotonFlux::GetRandomEnergyAndRadius(TVector2 eRange) { + std::pair result = {0, 0}; + if (!AreTablesLoaded()) { + RESTWarning << "Tables not loaded!!" << RESTendl; + return result; + } + Double_t rnd = fRandom->Rndm(); + for (unsigned int r = 0; r < fFluxTableIntegrals.size(); r++) { + if (rnd < fFluxTableIntegrals[r]) { + Double_t energy = fFluxTable[r]->GetRandom(); + if (eRange.X() != -1 && eRange.Y() != -1) { + if (energy < eRange.X() || energy > eRange.Y()) return GetRandomEnergyAndRadius(eRange); + } + Double_t radius = ((Double_t)r + fRandom->Rndm()) * 0.001; + std::pair p = {energy, radius}; + return p; + } + } + return result; +} + +/////////////////////////////////////////////// +/// \brief It prints on screen the flux table that has been read from file +/// +void TRestAxionSolarHiddenPhotonFlux::PrintContinuumSolarTable() { + cout << "Continuum solar flux table: " << endl; + cout << "--------------------------- " << endl; + for (unsigned int n = 0; n < fFluxTable.size(); n++) { + for (int m = 0; m < fFluxTable[n]->GetNbinsX(); m++) + cout << fFluxTable[n]->GetBinContent(m + 1) << "\t"; + cout << endl; + cout << endl; + } + cout << endl; +} + +/////////////////////////////////////////////// +/// \brief It prints on screen the integrated solar flux per solar ring +/// +void TRestAxionSolarHiddenPhotonFlux::PrintIntegratedRingFlux() { + cout << "Integrated solar flux per solar ring: " << endl; + cout << "--------------------------- " << endl; + /* + for (int n = 0; n < fFluxPerRadius.size(); n++) + cout << "n : " << n << " flux : " << fFluxPerRadius[n] << endl; + cout << endl; + */ +} + +/////////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestAxionSolarHiddenPhotonFlux +/// +void TRestAxionSolarHiddenPhotonFlux::PrintMetadata() { + TRestAxionSolarFlux::PrintMetadata(); + + RESTMetadata << " - Solar hidden photon datafile (flux) : " << fFluxDataFile << RESTendl; + RESTMetadata << " - Solar hidden photon datafile (width) : " << fWidthDataFile << RESTendl; + RESTMetadata << " - Solar hidden photon datafile (plasma freq) : " << fPlasmaFreqDataFile << RESTendl; + RESTMetadata << "-------" << RESTendl; + RESTMetadata << " - Total continuum flux : " << fTotalContinuumFlux << " cm-2 s-1" << RESTendl; + RESTMetadata << "++++++++++++++++++" << RESTendl; + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { + PrintContinuumSolarTable(); + PrintIntegratedRingFlux(); + } +} + +/////////////////////////////////////////////// +/// \brief It will create files with spectra to be used +/// in a later ocasion. +/// +void TRestAxionSolarHiddenPhotonFlux::ExportTables(Bool_t ascii) { + string rootFilename = TRestTools::GetFileNameRoot(fFluxDataFile); + + string path = REST_USER_PATH + "/export/"; + + if (!TRestTools::fileExists(path)) { + std::cout << "Creating path: " << path << std::endl; + int z = system(("mkdir -p " + path).c_str()); + if (z != 0) RESTError << "Could not create directory " << path << RESTendl; + } + + if (fFluxTable.size() > 0) { + std::vector> table; + for (const auto& x : fFluxTable) { + std::vector row; + for (int n = 0; n < x->GetNbinsX(); n++) row.push_back(x->GetBinContent(n + 1)); + + table.push_back(row); + } + + if (!ascii) + TRestTools::ExportBinaryTable(path + "/" + rootFilename + ".N200f", table); + else + TRestTools::ExportASCIITable(path + "/" + rootFilename + ".dat", table); + } +} + +/////////////////////////////////////////////// +/// \brief It draws the contents of a .flux file. This method just receives the +/// +TCanvas* TRestAxionSolarHiddenPhotonFlux::DrawSolarFlux() { + if (fCanvas != nullptr) { + delete fCanvas; + fCanvas = nullptr; + } + fCanvas = new TCanvas("canv", "This is the canvas title", 1200, 500); + fCanvas->Draw(); + + TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97); + pad1->Divide(2, 1); + pad1->Draw(); + + pad1->cd(1); + pad1->cd(1)->SetLogy(); + pad1->cd(1)->SetRightMargin(0.09); + pad1->cd(1)->SetLeftMargin(0.15); + pad1->cd(1)->SetBottomMargin(0.15); + + TH1D* ht = GetTotalSpectrum(); + ht->SetLineColor(kBlack); + ht->SetFillStyle(4050); + ht->SetFillColor(kBlue - 10); + + ht->Draw("hist"); + + pad1->cd(2); + pad1->cd(2)->SetRightMargin(0.09); + pad1->cd(2)->SetLeftMargin(0.15); + pad1->cd(2)->SetBottomMargin(0.15); + + ht->Draw("hist"); + + return fCanvas; +} diff --git a/src/TRestAxionSolarQCDFlux.cxx b/src/TRestAxionSolarQCDFlux.cxx index 5cdd2d46..79b5452e 100644 --- a/src/TRestAxionSolarQCDFlux.cxx +++ b/src/TRestAxionSolarQCDFlux.cxx @@ -265,7 +265,7 @@ Bool_t TRestAxionSolarQCDFlux::LoadTables() { } /////////////////////////////////////////////// -/// \brief A helper method to load the data file containning continuum spectra as a +/// \brief A helper method to load the data file containing continuum spectra as a /// function of the solar radius. It will be called by TRestAxionSolarQCDFlux::Initialize. /// void TRestAxionSolarQCDFlux::LoadContinuumFluxTable() { @@ -280,7 +280,7 @@ void TRestAxionSolarQCDFlux::LoadContinuumFluxTable() { RESTDebug << "Loading table from file : " << RESTendl; RESTDebug << "File : " << fullPathName << RESTendl; - std::vector> fluxTable; + std::vector> fluxTable; if (TRestTools::GetFileNameExtension(fFluxDataFile) == "dat") { std::vector> doubleTable; if (!TRestTools::ReadASCIITable(fullPathName, doubleTable)) { @@ -289,8 +289,8 @@ void TRestAxionSolarQCDFlux::LoadContinuumFluxTable() { return; } for (const auto& row : doubleTable) { - std::vector floatVec(row.begin(), row.end()); - fluxTable.push_back(floatVec); + std::vector doubleVec(row.begin(), row.end()); + fluxTable.push_back(doubleVec); } } else if (TRestTools::IsBinaryFile(fFluxDataFile)) TRestTools::ReadBinaryTable(fullPathName, fluxTable); @@ -309,14 +309,14 @@ void TRestAxionSolarQCDFlux::LoadContinuumFluxTable() { } for (unsigned int n = 0; n < fluxTable.size(); n++) { - TH1F* h = new TH1F(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), "", 200, 0, 20); + TH1D* h = new TH1D(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), "", 200, 0, 20); for (unsigned int m = 0; m < fluxTable[n].size(); m++) h->SetBinContent(m + 1, fluxTable[n][m]); fFluxTable.push_back(h); } } /////////////////////////////////////////////// -/// \brief A helper method to load the data file containning monochromatic spectral +/// \brief A helper method to load the data file containing monochromatic spectral /// lines as a function of the solar radius. It will be called by TRestAxionSolarQCDFlux::Initialize. /// void TRestAxionSolarQCDFlux::LoadMonoChromaticFluxTable() { @@ -339,10 +339,10 @@ void TRestAxionSolarQCDFlux::LoadMonoChromaticFluxTable() { return; } - std::vector> asciiTable; + std::vector> asciiTable; for (const auto& row : doubleTable) { - std::vector floatVec(row.begin(), row.end()); - asciiTable.push_back(floatVec); + std::vector doubleVec(row.begin(), row.end()); + asciiTable.push_back(doubleVec); } fFluxLines.clear(); @@ -355,8 +355,8 @@ void TRestAxionSolarQCDFlux::LoadMonoChromaticFluxTable() { } for (unsigned int en = 0; en < asciiTable[0].size(); en++) { - Float_t energy = asciiTable[0][en]; - TH1F* h = new TH1F(Form("%s_MonochromeFluxAtEnergy%6.4lf", GetName(), energy), "", 100, 0, 1); + Double_t energy = asciiTable[0][en]; + TH1D* h = new TH1D(Form("%s_MonochromeFluxAtEnergy%6.4lf", GetName(), energy), "", 100, 0, 1); for (unsigned int r = 1; r < asciiTable.size(); r++) h->SetBinContent(r, asciiTable[r][en]); fFluxLines[energy] = h; } @@ -399,12 +399,12 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { Double_t fluxBinSize = TRestTools::GetLowestIncreaseFromTable(fluxData, 1); for (const auto& data : fluxData) { - Float_t r = 0.005 + data[0]; - Float_t en = data[1] - 0.005; - Float_t flux = data[2] * fluxBinSize; // flux in cm-2 s-1 bin-1 + Double_t r = 0.005 + data[0]; + Double_t en = data[1] - 0.005; + Double_t flux = data[2] * fluxBinSize; // flux in cm-2 s-1 bin-1 - originalHist->Fill(r, en, (Float_t)flux); - continuumHist->Fill(r, en, (Float_t)flux); + originalHist->Fill(r, en, (Double_t)flux); + continuumHist->Fill(r, en, (Double_t)flux); } RESTDebug << "Histograms filled" << RESTendl; @@ -415,8 +415,8 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { const int smearPoints = (Int_t)(5 / (fBinSize * 100)); const int excludePoints = smearPoints / 5; for (const auto& data : fluxData) { - Float_t r = 0.005 + data[0]; - Float_t en = data[1] - 0.005; + Double_t r = 0.005 + data[0]; + Double_t en = data[1] - 0.005; Int_t binR = continuumHist->GetXaxis()->FindBin(r); Int_t binE = continuumHist->GetYaxis()->FindBin(en); @@ -430,8 +430,8 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { } avgFlux /= n; - Float_t targetBinFlux = continuumHist->GetBinContent(binR, binE); - Float_t thrFlux = avgFlux + fPeakSigma * TMath::Sqrt(avgFlux); + Double_t targetBinFlux = continuumHist->GetBinContent(binR, binE); + Double_t thrFlux = avgFlux + fPeakSigma * TMath::Sqrt(avgFlux); if (targetBinFlux > 0 && targetBinFlux > thrFlux) { continuumHist->SetBinContent(binR, binE, avgFlux); peaks++; @@ -441,8 +441,8 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { for (int n = 0; n < originalHist->GetNbinsX(); n++) for (int m = 0; m < originalHist->GetNbinsY(); m++) { - Float_t orig = originalHist->GetBinContent(n + 1, m + 1); - Float_t cont = continuumHist->GetBinContent(n + 1, m + 1); + Double_t orig = originalHist->GetBinContent(n + 1, m + 1); + Double_t cont = continuumHist->GetBinContent(n + 1, m + 1); spectrumHist->SetBinContent(n + 1, m + 1, orig - cont); } @@ -453,8 +453,8 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { fFluxTable.clear(); for (int n = 0; n < continuumHist->GetNbinsX(); n++) { - TH1F* hc = - (TH1F*)continuumHist->ProjectionY(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), n + 1, n + 1); + TH1D* hc = + (TH1D*)continuumHist->ProjectionY(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), n + 1, n + 1); fFluxTable.push_back(hc); } @@ -462,7 +462,7 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { for (int n = 0; n < spectrumHist->GetNbinsY(); n++) { if (spectrumHist->ProjectionX("", n + 1, n + 1)->Integral() > 0) { Double_t energy = spectrumHist->ProjectionY()->GetBinCenter(n + 1); - TH1F* hm = (TH1F*)spectrumHist->ProjectionX( + TH1D* hm = (TH1D*)spectrumHist->ProjectionX( Form("%s_MonochromeFluxAtEnergy%5.3lf", GetName(), energy), n + 1, n + 1); fFluxLines[energy] = hm; } @@ -475,13 +475,13 @@ void TRestAxionSolarQCDFlux::ReadFluxFile() { /// \brief It builds a histogram with the continuum spectrum component. /// The flux will be expressed in cm-2 s-1 keV-1. Binned in 100eV steps. /// -TH1F* TRestAxionSolarQCDFlux::GetContinuumSpectrum() { +TH1D* TRestAxionSolarQCDFlux::GetContinuumSpectrum() { if (fContinuumHist != nullptr) { delete fContinuumHist; fContinuumHist = nullptr; } - fContinuumHist = new TH1F("ContinuumHist", "", 200, 0, 20); + fContinuumHist = new TH1D("ContinuumHist", "", 200, 0, 20); for (const auto& x : fFluxTable) { fContinuumHist->Add(x); } @@ -501,13 +501,13 @@ TH1F* TRestAxionSolarQCDFlux::GetContinuumSpectrum() { /// \brief It builds a histogram with the monochromatic spectrum component. /// The flux will be expressed in cm-2 s-1 eV-1. Binned in 1eV steps. /// -TH1F* TRestAxionSolarQCDFlux::GetMonochromaticSpectrum() { +TH1D* TRestAxionSolarQCDFlux::GetMonochromaticSpectrum() { if (fMonoHist != nullptr) { delete fMonoHist; fMonoHist = nullptr; } - fMonoHist = new TH1F("MonochromaticHist", "", 20000, 0, 20); + fMonoHist = new TH1D("MonochromaticHist", "", 20000, 0, 20); for (const auto& x : fFluxLines) { fMonoHist->Fill(x.first, x.second->Integral()); // cm-2 s-1 eV-1 } @@ -528,16 +528,16 @@ TH1F* TRestAxionSolarQCDFlux::GetMonochromaticSpectrum() { /// spectrum component. The flux will be expressed in cm-2 s-1 keV-1. /// Binned in 1eV steps. /// -TH1F* TRestAxionSolarQCDFlux::GetTotalSpectrum() { - TH1F* hm = GetMonochromaticSpectrum(); - TH1F* hc = GetContinuumSpectrum(); +TH1D* TRestAxionSolarQCDFlux::GetTotalSpectrum() { + TH1D* hm = GetMonochromaticSpectrum(); + TH1D* hc = GetContinuumSpectrum(); if (fTotalHist != nullptr) { delete fTotalHist; fTotalHist = nullptr; } - fTotalHist = new TH1F("fTotalHist", "", 20000, 0, 20); + fTotalHist = new TH1D("fTotalHist", "", 20000, 0, 20); for (int n = 0; n < hc->GetNbinsX(); n++) { for (int m = 0; m < 100; m++) { fTotalHist->SetBinContent(n * 100 + 1 + m, hc->GetBinContent(n + 1)); @@ -580,12 +580,12 @@ TCanvas* TRestAxionSolarQCDFlux::DrawSolarFlux() { pad1->cd(1)->SetLeftMargin(0.15); pad1->cd(1)->SetBottomMargin(0.15); - TH1F* ht = GetTotalSpectrum(); + TH1D* ht = GetTotalSpectrum(); ht->SetLineColor(kBlack); ht->SetFillStyle(4050); ht->SetFillColor(kBlue - 10); - TH1F* hm = GetMonochromaticSpectrum(); + TH1D* hm = GetMonochromaticSpectrum(); hm->SetLineColor(kBlack); hm->Scale(100); // renormalizing per 100eV-1 @@ -634,7 +634,7 @@ void TRestAxionSolarQCDFlux::IntegrateSolarFluxes() { /////////////////////////////////////////////// /// \brief It returns the integrated flux at earth in cm-2 s-1 for the given energy range /// -Double_t TRestAxionSolarQCDFlux::IntegrateFluxInRange(TVector2 eRange, Double_t mass) { +Double_t TRestAxionSolarQCDFlux::IntegrateFluxInRange(TVector2 eRange) { if (eRange.X() == -1 && eRange.Y() == -1) { if (GetTotalFlux() == 0) IntegrateSolarFluxes(); return GetTotalFlux(); @@ -658,8 +658,7 @@ Double_t TRestAxionSolarQCDFlux::IntegrateFluxInRange(TVector2 eRange, Double_t /// \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 TRestAxionSolarQCDFlux::GetRandomEnergyAndRadius(TVector2 eRange, - Double_t mass) { +std::pair TRestAxionSolarQCDFlux::GetRandomEnergyAndRadius(TVector2 eRange) { std::pair result = {0, 0}; if (!AreTablesLoaded()) return result; Double_t rnd = fRandom->Rndm(); @@ -771,9 +770,9 @@ void TRestAxionSolarQCDFlux::ExportTables(Bool_t ascii) { } if (fFluxTable.size() > 0) { - std::vector> table; + std::vector> table; for (const auto& x : fFluxTable) { - std::vector row; + std::vector row; for (int n = 0; n < x->GetNbinsX(); n++) row.push_back(x->GetBinContent(n + 1)); table.push_back(row); @@ -786,9 +785,9 @@ void TRestAxionSolarQCDFlux::ExportTables(Bool_t ascii) { } if (fFluxLines.size() > 0) { - std::vector> table; + std::vector> table; for (const auto& x : fFluxLines) { - std::vector row; + std::vector row; row.push_back(x.first); for (int n = 0; n < x.second->GetNbinsX(); n++) row.push_back(x.second->GetBinContent(n + 1));