diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index e8ed7a56..11e22758 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -1,5 +1,5 @@ Please, REMOVE this text! This is just a template to serve as a reminder. - + Please, give as much detail as possible at your PR and create the simplest possible PR to facilitate reviewing. If your PR fixes a given issue, please specify the number followed by # diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index bb3810c6..53dd9fbc 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -30,7 +30,7 @@ validateLibrary: - $CRONJOB == "YES" build: - type: build + stage: build script: - echo "**${CI_PROJECT_DIR}**" - rm -rf ${CI_PROJECT_DIR}/install @@ -38,13 +38,9 @@ build: - cd framework - ./scripts/checkoutRemoteBranch.sh ${CI_COMMIT_BRANCH} - git submodule init source/libraries/axion - # This will download the data from axionlib. - # erhaps we should make the data accessible throught https and then open the remote files from the pipeline? - git submodule update source/libraries/axion - cd source/libraries/axion/ - git checkout ${CI_COMMIT_BRANCH} - - git submodule init data - - git submodule update data - cd ../../../ - mkdir build - cd build @@ -63,7 +59,7 @@ build: expire_in: 1 day loadRESTLibs: - type: loadRESTLibs + stage: loadRESTLibs script: - . ${CI_PROJECT_DIR}/install/thisREST.sh # - . ${CI_PROJECT_DIR}/framework/source/libraries/axion/external/solarAxionFlux/bin/thisSolarAxionFluxLib.sh @@ -73,14 +69,16 @@ loadRESTLibs: - $CRONJOB magneticField: - type: metadata + stage: metadata script: # - . ${CI_PROJECT_DIR}/framework/source/libraries/axion/external/solarAxionFlux/bin/thisSolarAxionFluxLib.sh - . ${CI_PROJECT_DIR}/install/thisREST.sh - cd ${CI_PROJECT_DIR}/pipeline/metadata/magneticField/ - - ls $REST_PATH/data/axion/magneticField/ + - wget https://sultan.unizar.es/axionlib-data/magneticField/fields.rml + - wget https://sultan.unizar.es/axionlib-data/magneticField/Bykovskiy_201906.dat - python magneticField.py - cd trilinear + - wget https://sultan.unizar.es/axionlib-data/magneticField/fields.rml - restRoot -b -q GetMagneticField_test.C - cd ${CI_PROJECT_DIR}/pipeline/metadata/magneticField/boundary/ - restRoot -b -q Boundaries_test.C @@ -89,10 +87,16 @@ magneticField: - $CRONJOB optics: - type: metadata + stage: metadata script: - . ${CI_PROJECT_DIR}/install/thisREST.sh - cd ${CI_PROJECT_DIR}/pipeline/metadata/optics/ + - wget https://sultan.unizar.es/axionlib-data/optics/mirrors.rml + - wget https://sultan.unizar.es/axionlib-data/optics/optics.rml + - wget https://sultan.unizar.es/axionlib-data/optics/Reflectivity_Single_C_30_SiO2_0.N901f + - wget https://sultan.unizar.es/axionlib-data/optics/Transmission_Single_C_30_SiO2_0.N901f + - wget https://sultan.unizar.es/axionlib-data/optics/Reflectivity_Single_Au_250_Ni_0.4.N901f + - wget https://sultan.unizar.es/axionlib-data/optics/Transmission_Single_Au_250_Ni_0.4.N901f - python mirrors.py - python optics.py - python basic.py @@ -100,21 +104,29 @@ optics: variables: - $CRONJOB - xRayTransmission: - type: metadata + stage: metadata script: - . ${CI_PROJECT_DIR}/install/thisREST.sh - cd ${CI_PROJECT_DIR}/pipeline/metadata/transmission/ + - wget https://sultan.unizar.es/axionlib-data/transmission/windows.rml + - wget https://sultan.unizar.es/axionlib-data/transmission/Al.sol + - wget https://sultan.unizar.es/axionlib-data/transmission/Si.sol + - wget https://sultan.unizar.es/axionlib-data/transmission/Si3N4.sol - export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${CI_PROJECT_DIR}/mpfr-4.0.2/install/lib - python windowPlot.py solarFlux: - type: metadata + stage: metadata script: - . ${CI_PROJECT_DIR}/install/thisREST.sh - cd ${CI_PROJECT_DIR}/pipeline/metadata/solarFlux/ - - python solarFlux.py + - wget https://sultan.unizar.es/axionlib-data/solarFlux/fluxes.rml + - wget https://sultan.unizar.es/axionlib-data/solarFlux/Dummy_Galan_202202.spt + - wget https://sultan.unizar.es/axionlib-data/solarFlux/Primakoff_Gianotti_201904.dat + - wget https://sultan.unizar.es/axionlib-data/solarFlux/Primakoff_LennertHoof_202203.dat + - python solarTests.py + - python solarPlot.py - python compare.py except: variables: diff --git a/data b/data index 227fb74b..22cf5a4d 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 227fb74b9bc0ca440b0896e24cf6e786a5488eea +Subproject commit 22cf5a4dd9f3b07e0d0017c767c62839988c99c7 diff --git a/images/ABC_FluxTable.png b/images/ABC_FluxTable.png new file mode 100644 index 00000000..1453af39 Binary files /dev/null and b/images/ABC_FluxTable.png differ diff --git a/images/ABC_flux_MC.png b/images/ABC_flux_MC.png new file mode 100644 index 00000000..ac7e002f Binary files /dev/null and b/images/ABC_flux_MC.png differ diff --git a/inc/TRestAxionSolarFlux.h b/inc/TRestAxionSolarFlux.h index e435e940..0e88c17b 100644 --- a/inc/TRestAxionSolarFlux.h +++ b/inc/TRestAxionSolarFlux.h @@ -23,7 +23,9 @@ #ifndef _TRestAxionSolarFlux #define _TRestAxionSolarFlux -#include +#include +#include +#include #include #include @@ -35,8 +37,6 @@ class TRestAxionSolarFlux : public TRestMetadata { private: void Initialize(); - void InitFromConfigFile(); - /// The filename containning the solar flux table with continuum spectrum std::string fFluxDataFile = ""; //< @@ -52,11 +52,17 @@ class TRestAxionSolarFlux : public TRestMetadata { /// Seed used in random generator Int_t fSeed = 0; //< - /// The tabulated solar flux continuum spectra TH1D(100,0,20)keV in cm-2 s-1 keV-1 versus solar radius - std::vector fFluxTable; //! + /// It will be used when loading `.flux` files to define the input file energy binsize in eV. + Double_t fBinSize = 0; //< + + /// It will be used when loading `.flux` files to define the threshold for peak identification + Double_t fPeakSigma = 0; //< + + /// The tabulated solar flux continuum spectra TH1F(100,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; //! @@ -76,6 +82,22 @@ class TRestAxionSolarFlux : public TRestMetadata { /// Random number generator TRandom3* fRandom = nullptr; //! + /// A canvas pointer for drawing + TCanvas* fCanvas = nullptr; //! + + /// A pointer to the continuum spectrum histogram + TH1F* fContinuumHist = nullptr; //! + + /// A pointer to the monochromatic spectrum histogram + TH1F* fMonoHist = nullptr; //! + + /// A pointer to the superposed monochromatic and continuum spectrum histogram + TH1F* fTotalHist = nullptr; //! + + /// A metadata member to control if the tables have been loaded + Bool_t fTablesLoaded = false; //! + + void ReadFluxFile(); void LoadContinuumFluxTable(); void LoadMonoChromaticFluxTable(); void IntegrateSolarFluxes(); @@ -92,19 +114,23 @@ class TRestAxionSolarFlux : public TRestMetadata { std::pair GetRandomEnergyAndRadius(); + void LoadTables(); + + TH1F* GetContinuumSpectrum(); + TH1F* GetMonochromaticSpectrum(); + TH1F* GetTotalSpectrum(); + + TH1F* GetFluxHistogram(std::string fname, Double_t binSize = 0.01); + TCanvas* DrawFluxFile(std::string fname, Double_t binSize = 0.01); + TCanvas* DrawSolarFlux(); + /// Tables might be loaded using a solar model description by TRestAxionSolarModel void InitializeSolarTable(TRestAxionSolarModel* model) { // TOBE implemented // This method should initialize the tables fFluxTable and fFluxLines } - void ExportTables(std::string fname) { - // TOBE implemented. Creates fname.dat and fname.spt - // If we have external methods to initialize solar flux tables this method - // might be used to generate the tables that can be used later on directly - // - // Check data/solarFlux/README.md for data format and file naming conventions - } + void ExportTables(Bool_t ascii = false); void PrintMetadata(); diff --git a/pipeline/.gitignore b/pipeline/.gitignore new file mode 100644 index 00000000..e33609d2 --- /dev/null +++ b/pipeline/.gitignore @@ -0,0 +1 @@ +*.png diff --git a/pipeline/metadata/magneticField/fields.rml b/pipeline/metadata/magneticField/fields.rml deleted file mode 100644 index a58a7623..00000000 --- a/pipeline/metadata/magneticField/fields.rml +++ /dev/null @@ -1,40 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/pipeline/metadata/magneticField/trilinear/GetMagneticField_test.C b/pipeline/metadata/magneticField/trilinear/GetMagneticField_test.C index b19db3c0..2ddf65b1 100644 --- a/pipeline/metadata/magneticField/trilinear/GetMagneticField_test.C +++ b/pipeline/metadata/magneticField/trilinear/GetMagneticField_test.C @@ -12,7 +12,7 @@ Int_t CheckPoint(TVector3 point); Int_t GetMagneticField_test() { TVector3 coordinates; - myField = new TRestAxionMagneticField("../fields.rml", "bField"); + myField = new TRestAxionMagneticField("fields.rml", "bField"); // changing x, while y and z are constant coordinates = TVector3(-349.8, -325.0, -4975.0); diff --git a/pipeline/metadata/optics/basic.py b/pipeline/metadata/optics/basic.py index bc43c111..ce0006f3 100755 --- a/pipeline/metadata/optics/basic.py +++ b/pipeline/metadata/optics/basic.py @@ -27,8 +27,8 @@ totalSamples = 10000 genSize = 80 -basicOptics = ROOT.TRestAxionOptics("basic.rml", "basic") -spiderOptics = ROOT.TRestAxionOptics("basic.rml", "basic_spider") +basicOptics = ROOT.TRestAxionOptics("optics.rml", "basic") +spiderOptics = ROOT.TRestAxionOptics("optics.rml", "basic_spider") rings = basicOptics.GetNumberOfRings() print( "Number of rings (no-spider): " + str(rings) ) diff --git a/pipeline/metadata/optics/basic.rml b/pipeline/metadata/optics/basic.rml deleted file mode 100644 index c0d8e262..00000000 --- a/pipeline/metadata/optics/basic.rml +++ /dev/null @@ -1,36 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/pipeline/metadata/optics/mirrors.rml b/pipeline/metadata/optics/mirrors.rml deleted file mode 100644 index b2403114..00000000 --- a/pipeline/metadata/optics/mirrors.rml +++ /dev/null @@ -1,35 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/pipeline/metadata/optics/optics.py b/pipeline/metadata/optics/optics.py index 151f63ba..7ab66e49 100755 --- a/pipeline/metadata/optics/optics.py +++ b/pipeline/metadata/optics/optics.py @@ -5,7 +5,7 @@ ROOT.gSystem.Load("libRestFramework.so") ROOT.gSystem.Load("libRestAxion.so") -mcplOptics_1 = ROOT.TRestAxionMCPLOptics("setups.rml", "mcpl") +mcplOptics_1 = ROOT.TRestAxionMCPLOptics("optics.rml", "mcpl") mcplOptics_1.PrintMetadata() rings = mcplOptics_1.GetNumberOfRings() diff --git a/pipeline/metadata/optics/setups.rml b/pipeline/metadata/optics/setups.rml deleted file mode 100644 index 51489cec..00000000 --- a/pipeline/metadata/optics/setups.rml +++ /dev/null @@ -1,19 +0,0 @@ - - - - - - - - - - - - - - - - - - - diff --git a/pipeline/metadata/solarFlux/compare.py b/pipeline/metadata/solarFlux/compare.py index 00da2726..6b1d45b8 100755 --- a/pipeline/metadata/solarFlux/compare.py +++ b/pipeline/metadata/solarFlux/compare.py @@ -25,7 +25,7 @@ pad1.Divide(2,1) pad1.Draw() -primakoffLH = ROOT.TRestAxionSolarFlux("fluxes.rml", "LennertHoof") +primakoffLH = ROOT.TRestAxionSolarFlux("fluxes.rml", "LennertHoofPrimakoff") primakoffG = ROOT.TRestAxionSolarFlux("fluxes.rml", "Gianotti") if primakoffLH.GetError(): diff --git a/pipeline/metadata/solarFlux/fluxes.rml b/pipeline/metadata/solarFlux/fluxes.rml deleted file mode 100644 index edcc41b2..00000000 --- a/pipeline/metadata/solarFlux/fluxes.rml +++ /dev/null @@ -1,43 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/pipeline/metadata/solarFlux/solarFlux.py b/pipeline/metadata/solarFlux/solarFlux.py deleted file mode 100755 index 88f4488c..00000000 --- a/pipeline/metadata/solarFlux/solarFlux.py +++ /dev/null @@ -1,119 +0,0 @@ -#!/usr/bin/python3 - -totalSamples = 20000 - -outfname = "flux.png" - -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) - -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.TRestAxionSolarFlux("fluxes.rml", "combined") - -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", 200, 0, 20, 100, 0, 1 ) -for x in range(totalSamples): - x = combinedFlux.GetRandomEnergyAndRadius() - 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(totalSamples): - angle = rnd.Rndm() * 2 * math.pi - x = combinedFlux.GetRandomEnergyAndRadius() - - solarDisk.Fill( x[1]*math.cos(angle), x[1]*math.sin(angle) ) - -pad1.cd(1) -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("box") - -pad1.cd(2) -enSpt = comb_spt.ProjectionX() -enSpt.SetTitle("Energy spectrum") -enSpt.GetYaxis().SetTitleSize(0.05); -enSpt.SetStats(0) -enSpt.Draw() -if ( enSpt.GetMaximumBin() != 41 ): - exit(1) - -pad1.cd(3) -rSpt = comb_spt.ProjectionY() -rSpt.SetTitle("Radial distribution") -rSpt.GetYaxis().SetTitleSize(0.05); -rSpt.SetStats(0) -rSpt.Draw() -if ( rSpt.GetMaximumBin() != 25 ): - exit(2) - -pad1.cd(4) -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() - -c1.Print(outfname) - -monoFlux = ROOT.TRestAxionSolarFlux("fluxes.rml", "mono") - -if monoFlux.GetError(): - print ( monoFlux.GetErrorMessage() ) - print ( "\nSolar flux initialization failed! Exit code : 101" ) - exit(101) - -x = monoFlux.GetRandomEnergyAndRadius() -if ( int( x[0]*100 ) != 800 or int( x[1]*100 ) != 58 ): - print ( "\nMonochromatic flux values seem to be wrong! Exit code : 201" ) - exit(201) - -print ("[\033[92m OK \x1b[0m]") - -continuumFlux = ROOT.TRestAxionSolarFlux("fluxes.rml", "Gianotti") - -if continuumFlux.GetError(): - print ( continuumFlux.GetErrorMessage() ) - print ( "\nSolar flux initialization failed! Exit code : 101" ) - exit(101) - -x = monoFlux.GetRandomEnergyAndRadius() -if ( int( x[0]*100 ) != 400 or int( x[1]*100 ) != 24 ): - print ( "\nMonochromatic flux values seem to be wrong! Exit code : 202" ) - exit(202) - -print ("All tests passed! [\033[92m OK \x1b[0m]") - -print ("") - -exit(0) diff --git a/pipeline/metadata/solarFlux/solarPlot.py b/pipeline/metadata/solarFlux/solarPlot.py new file mode 100755 index 00000000..78e16b87 --- /dev/null +++ b/pipeline/metadata/solarFlux/solarPlot.py @@ -0,0 +1,147 @@ +#!/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') +args = parser.parse_args() + +rmlfile = "fluxes.rml" +if args.rmlfile != None: + rmlfile = args.rmlfile + +outfname = "flux.png" +if args.outfname != None: + outfname = args.outfname + +fluxname = "combined" +if args.fluxname != None: + fluxname = args.fluxname + + +samples = 20000 +if args.samples != None: + samples = args.samples + +validation = False +if rmlfile == "fluxes.rml" and fluxname == "combined" and outfname == "flux.png" and samples == 20000: + 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.TRestAxionSolarFlux(rmlfile, fluxname) +combinedFlux.LoadTables() +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() + 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() + + 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! Exit code : 1" ) + exit(1) + +print ( "\nMaximum radius bin is " + str(rSpt.GetMaximumBin() ) ) + +if validation: + if ( rSpt.GetMaximumBin() != 25 ): + print ( "\nMaximum Bin is not the expected one! Exit code : 2" ) + exit(2) + +exit(0) diff --git a/pipeline/metadata/solarFlux/solarTests.py b/pipeline/metadata/solarFlux/solarTests.py new file mode 100755 index 00000000..80926871 --- /dev/null +++ b/pipeline/metadata/solarFlux/solarTests.py @@ -0,0 +1,53 @@ +#!/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) + + +ROOT.gSystem.Load("libRestFramework.so") +ROOT.gSystem.Load("libRestAxion.so") + +monoFlux = ROOT.TRestAxionSolarFlux("fluxes.rml", "mono") +monoFlux.LoadTables() +monoFlux.PrintMetadata() + +if monoFlux.GetError(): + print ( monoFlux.GetErrorMessage() ) + print ( "\nSolar flux initialization failed! Exit code : 101" ) + exit(101) + +x = monoFlux.GetRandomEnergyAndRadius() +print( x[0] ) +print( x[1] ) +if ( int( x[0]*100 ) != 800 or int( x[1]*100 ) != 83 ): + print ( "\nMonochromatic flux values seem to be wrong! Exit code : 201" ) + exit(201) + +print ("[\033[92m OK \x1b[0m]") + +continuumFlux = ROOT.TRestAxionSolarFlux("fluxes.rml", "Gianotti") +continuumFlux.LoadTables() +continuumFlux.PrintMetadata() + +if continuumFlux.GetError(): + print ( continuumFlux.GetErrorMessage() ) + print ( "\nSolar flux initialization failed! Exit code : 101" ) + exit(101) + +x = monoFlux.GetRandomEnergyAndRadius() +if ( int( x[0]*100 ) != 400 or int( x[1]*100 ) != 24 ): + print ( "\nContinuum flux values seem to be wrong! Exit code : 202" ) + exit(202) + +print ("All tests passed! [\033[92m OK \x1b[0m]") + +print ("") + +exit(0) diff --git a/pipeline/metadata/transmission/windows.rml b/pipeline/metadata/transmission/windows.rml deleted file mode 100644 index 290081e9..00000000 --- a/pipeline/metadata/transmission/windows.rml +++ /dev/null @@ -1,32 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/scripts/fluxToTables.py b/scripts/fluxToTables.py index acaf8105..ae63d4c7 100755 --- a/scripts/fluxToTables.py +++ b/scripts/fluxToTables.py @@ -10,10 +10,10 @@ import sys, os import ROOT +from array import array from ROOT import TH1D, TH2D import pandas as pd - import argparse parser = argparse.ArgumentParser() @@ -41,11 +41,13 @@ print (args.fname.find(".")) outfname = os.path.basename(args.fname) -outfname = outfname[0:outfname.find(".")] + ".dat" +outfname = outfname[0:outfname.find(".")] + ".Nxxf" if args.outfname != None: - outfname = args.outfname + ".dat" + outfname = args.outfname + "Nxxf" print (outfname) +extension = outfname[outfname.find(".")+1:] +print (extension) data = pd.read_csv(args.fname, sep=" ", skiprows=skiprows, header=None) @@ -65,13 +67,24 @@ #print( "R: " + str(r) + " En: " + str(en) + " Flux: " + str(flux) + "\n" ) print (outfname) -outf = open( outfname, "w") -print( str(fluxTable.GetNbinsX()) + " " + str(fluxTable.GetNbinsY()) ) -for n in range(0,fluxTable.GetNbinsX()): - for m in range(0,fluxTable.GetNbinsY()): - outf.write( str(fluxTable.GetBinContent( n+1, m+1 )) ) - if m == fluxTable.GetNbinsY() - 1: - outf.write("\n") - else: - outf.write("\t") -outf.close() + +if extension == "Nxxf": + print("Writting binary" ) + output_file = open(outfname, 'wb') + for n in range(0,fluxTable.GetNbinsX()): + float_array = array('f', []) + for m in range(0,fluxTable.GetNbinsY()): + float_array.append(float(fluxTable.GetBinContent( n+1, m+1 ))) + float_array.tofile(output_file) + output_file.close() +else: + outf = open( outfname, "w") + print( str(fluxTable.GetNbinsX()) + " " + str(fluxTable.GetNbinsY()) ) + for n in range(0,fluxTable.GetNbinsX()): + for m in range(0,fluxTable.GetNbinsY()): + outf.write( str(fluxTable.GetBinContent( n+1, m+1 )) ) + if m == fluxTable.GetNbinsY() - 1: + outf.write("\n") + else: + outf.write("\t") + outf.close() diff --git a/src/TRestAxionOpticsMirror.cxx b/src/TRestAxionOpticsMirror.cxx index 7e995001..f48dab20 100644 --- a/src/TRestAxionOpticsMirror.cxx +++ b/src/TRestAxionOpticsMirror.cxx @@ -285,7 +285,7 @@ void TRestAxionOpticsMirror::LoadTables() { /// members /// std::string TRestAxionOpticsMirror::GetReflectivityFilename() { - string fnameR = "Reflectivy_" + fMirrorType + "_" + fLayerTop + "_" + fLayerThicknessTop + "_" + + string fnameR = "Reflectivity_" + fMirrorType + "_" + fLayerTop + "_" + fLayerThicknessTop + "_" + fSubstrate + "_" + fSigmaTop + ".N901f"; if (fMirrorType == "Bilayer") fnameR = "Reflectivity_" + fMirrorType + "_" + fLayerTop + "_" + fLayerThicknessTop + "_" + diff --git a/src/TRestAxionSolarFlux.cxx b/src/TRestAxionSolarFlux.cxx index fe1b8fb2..22df80cb 100644 --- a/src/TRestAxionSolarFlux.cxx +++ b/src/TRestAxionSolarFlux.cxx @@ -42,13 +42,48 @@ /// given in the solar flux tables. /// - *fluxDataFile:* A table with 100 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. +/// for the range (0,20)keV in steps of 100eV. The table could be provided in ASCII format, +/// using `.dat` extension, or it might be a binary table using `.N200f` extension. /// - *fluxSptFile:* A table where each column represents a monochromatic energy. The /// column contains 101 rows, the first element is the energy of the monochromatic line /// while the next 100 elements contain the flux, measured in cm-2 s-1, integrated to /// each solar ring, being the second element the ring in the center of the sun. /// -/// The following code shows how to define this class inside a RML file. +/// Additionally this class will be able to read `.flux` files that are the original files +/// produced in 3-columns format (inner radius [solar units] / energy [keV] / +/// flux [cm-2 s-1 keV-1]). The `.flux` files may contain the full information, +/// continuum and spectral components. Those components will be splited into two independent +/// contributions by TRestAxionSolarFlux::ReadFluxFile to be managed internally. Two +/// additional parameters *will be required* to translate the `.flux` files into the tables +/// that are understood by this class. +/// - *binSize:* The energy binning used on the `.flux` file and inside the histogram used +/// for monochromatic lines identification. +/// - *peakSigma:* The ratio between the flux provided at the `.flux` file and the +/// average flux calculated in the peak surroundings. If the flux ratio is higher than +/// this value, the flux at that particular bin will be considered a peak. +/// +/// Optionally, if we want to consider a different binning on the monochromatic/continuum +/// histogram used internally for the calculation we may specify optionally a new parameter. +/// In that case, `fBinSize` will be the binning of the internal histogram, while the new +/// parameter will be the binning given inside the `.flux` file. +/// - *fluxBinSize:* The bin size used on the `.flux` table. +/// +/// 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 /// @@ -59,16 +94,97 @@ /// /// \endcode /// -/// The previous definition was used to generate the following figure using the script -/// pipeline/solarFlux/solarFlux.py. +/// \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 TRestAxionSolarFlux::LoadTables to trigger the tables +/// initialization. /// -/// \htmlonly \endhtmlonly +/// ### Performing MonteCarlo tests using pre-loaded tables /// -/// ![Solar flux distributions generated with TRestAxionSolarFlux.](AxionSolarFlux.png) +/// 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 TRestAxionSolarFlux::GetRandomEnergyAndRadius. +/// +/// \code +/// python3 solarPlot.py --fluxname LennertHoofABC --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 TRestAxionSolarFlux.](ABC_flux_MC.png) +/// +/// ### Reading solar flux tables from `.flux` files +/// +/// In a similar way we may initialize the class using a `.flux` file. The `.flux` described +/// previously will contain a high definition flux table measured in `cm-2 s-1 keV-1` as a +/// function of the energy (from 0 to 20keV) and the inner solar radius (from 0 to 1). The +/// binning of this table will be typically between 1eV and 10eV. +/// +/// The class TRestAxionSolarFlux will be initialized directly using the `.flux` file +/// provided under the `fluxDataFile` parameter. During the initialization, the flux will be +/// splitted into two independent flux components. One smooth continuum component integrated +/// in 100 eV steps, and a monochromatic peak components. +/// +/// In order to help with the identification of peaks we need to define the `binSize` used in +/// the `.flux` table and the `peakSigma` defining the number of sigmas over the average for +/// a bin to be considered a peak. +/// +/// \code +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// \endcode +/// +/// We will be able to load this file as usual, using the following recipe inside +/// `restRoot`, +/// +/// \code +/// TRestAxionSolarFlux *sFlux = new TRestAxionSolarFlux("fluxes.rml", "LennertHoofABC") +/// sFlux->LoadTables() +/// TCanvas *c = sFlux->DrawSolarFluxes() +/// c->Print("ABC_FluxTable.png" ) +/// \endcode +/// +/// will generate the following figure. +/// +/// \htmlonly \endhtmlonly +/// +/// ![Solar flux distributions generated with TRestAxionSolarFlux::DrawSolarFlux.](ABC_FluxTable.png) +/// +/// ### Exporting the solar flux tables +/// +/// On top of that, we will be able to export those tables to the TRestAxionSolarFlux standard +/// format to be used in later occasions. +/// +/// \code +/// TRestAxionSolarFlux *sFlux = new TRestAxionSolarFlux("fluxes.rml", "LennertHoofABC") +/// sFlux->LoadTables() +/// sFlux->ExportTables() +/// \endcode +/// +/// which will produce two files, a binary table `.N200f` with the continuum flux, and an ASCII +/// table containning the `.spt` monochromatic lines. 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 TRestAxionSolarFlux::ExportTables method call. The files will be placed +/// at the REST user space, at `$HOME/.rest/export/` directory. /// /// TODO Implement the method TRestAxionSolarFlux::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 @@ -109,8 +225,6 @@ TRestAxionSolarFlux::TRestAxionSolarFlux() : TRestMetadata() {} /// corresponding TRestAxionMagneticField section inside the RML. /// TRestAxionSolarFlux::TRestAxionSolarFlux(const char* cfgFileName, string name) : TRestMetadata(cfgFileName) { - cout << "Entering TRestAxionSolarFlux constructor( cfgFileName, name )" << endl; - LoadConfigFromFile(fConfigFileName, name); if (GetVerboseLevel() >= REST_Info) PrintMetadata(); @@ -128,9 +242,23 @@ void TRestAxionSolarFlux::Initialize() { SetSectionName(this->ClassName()); SetLibraryVersion(LIBRARY_VERSION); - LoadContinuumFluxTable(); + fTablesLoaded = false; + LoadTables(); +} - LoadMonoChromaticFluxTable(); +/////////////////////////////////////////////// +/// \brief It will load the tables in memory by using the filename information provided +/// inside the metadata members. +/// +void TRestAxionSolarFlux::LoadTables() { + if (fFluxDataFile == "" && fFluxSptFile == "") return; + + if (TRestTools::GetFileNameExtension(fFluxDataFile) == "flux") { + ReadFluxFile(); + } else { + LoadContinuumFluxTable(); + LoadMonoChromaticFluxTable(); + } IntegrateSolarFluxes(); @@ -141,6 +269,8 @@ void TRestAxionSolarFlux::Initialize() { fRandom = new TRandom3(fSeed); if (fSeed == 0) fSeed = fRandom->GetSeed(); + + fTablesLoaded = true; } /////////////////////////////////////////////// @@ -149,7 +279,9 @@ void TRestAxionSolarFlux::Initialize() { /// void TRestAxionSolarFlux::LoadContinuumFluxTable() { if (fFluxDataFile == "") { - debug << "TRestAxionSolarflux::LoadContinuumFluxTable. No solar flux table was defined" << endl; + debug << "TRestAxionSolarflux::LoadContinuumFluxTable. No solar flux continuum table was " + "defined" + << endl; return; } @@ -158,18 +290,31 @@ void TRestAxionSolarFlux::LoadContinuumFluxTable() { debug << "Loading table from file : " << endl; debug << "File : " << fullPathName << endl; - std::vector> fluxTable; - TRestTools::ReadASCIITable(fullPathName, fluxTable); + std::vector> fluxTable; + if (TRestTools::GetFileNameExtension(fFluxDataFile) == "dat") { + std::vector> doubleTable; + TRestTools::ReadASCIITable(fullPathName, doubleTable); + for (const auto& row : doubleTable) { + std::vector floatVec(row.begin(), row.end()); + fluxTable.push_back(floatVec); + } + } else if (TRestTools::IsBinaryFile(fFluxDataFile)) + TRestTools::ReadBinaryTable(fullPathName, fluxTable); + else { + fluxTable.clear(); + ferr << "Filename extension was not recognized!" << endl; + ferr << "Solar flux table will not be populated" << endl; + } if (fluxTable.size() != 100 && fluxTable[0].size() != 200) { fluxTable.clear(); ferr << "LoadContinuumFluxTable. The table does not contain the right number of rows or columns" << endl; - ferr << "Table will not be populated" << endl; + ferr << "Solar flux table will not be populated" << endl; } for (int n = 0; n < fluxTable.size(); n++) { - TH1D* h = new TH1D(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), "", 200, 0, 20); + TH1F* h = new TH1F(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), "", 200, 0, 20); for (int m = 0; m < fluxTable[n].size(); m++) h->SetBinContent(m + 1, fluxTable[n][m]); fFluxTable.push_back(h); } @@ -192,8 +337,14 @@ void TRestAxionSolarFlux::LoadMonoChromaticFluxTable() { debug << "Loading monochromatic lines from file : " << endl; debug << "File : " << fullPathName << endl; - std::vector> asciiTable; - TRestTools::ReadASCIITable(fullPathName, asciiTable); + std::vector> doubleTable; + TRestTools::ReadASCIITable(fullPathName, doubleTable); + + std::vector> asciiTable; + for (const auto& row : doubleTable) { + std::vector floatVec(row.begin(), row.end()); + asciiTable.push_back(floatVec); + } fFluxLines.clear(); @@ -204,14 +355,301 @@ void TRestAxionSolarFlux::LoadMonoChromaticFluxTable() { } for (int en = 0; en < asciiTable[0].size(); en++) { - Double_t energy = asciiTable[0][en]; - std::vector profile; - TH1D* h = new TH1D(Form("%s_MonochromeFluxAtEnergy%4.2lf", GetName(), energy), "", 100, 0, 1); + Float_t energy = asciiTable[0][en]; + TH1F* h = new TH1F(Form("%s_MonochromeFluxAtEnergy%6.4lf", GetName(), energy), "", 100, 0, 1); for (int r = 1; r < asciiTable.size(); r++) h->SetBinContent(r, asciiTable[r][en]); fFluxLines[energy] = h; } } +/////////////////////////////////////////////// +/// \brief It loads a .flux file. It will split continuum and monochromatic peaks, loading +/// both internal flux tables. +/// +void TRestAxionSolarFlux::ReadFluxFile() { + if (fBinSize <= 0) { + ferr << "TRestAxionSolarflux::ReadFluxFile. Energy bin size of .flux file must be specified." << endl; + ferr << "Please, define binSize parameter in eV." << endl; + return; + } + + if (fPeakSigma <= 0) { + warning << "TRestAxionSolarflux::ReadFluxFile. Peak sigma must be specified to generate " + "monochromatic spectrum." + << endl; + warning + << "Only continuum table will be generated. If this was intentional, please, ignore this warning." + << endl; + return; + } + + string fullPathName = SearchFile((string)fFluxDataFile); + + debug << "Loading flux table ... " << endl; + debug << "File : " << fullPathName << endl; + std::vector> fluxData; + TRestTools::ReadASCIITable(fullPathName, fluxData, 3); + + debug << "Table loaded" << endl; + TH2F* originalHist = new TH2F("FullTable", "", 100, 0., 1., (Int_t)(20. / fBinSize), 0., 20.); + TH2F* continuumHist = new TH2F("ContinuumTable", "", 100, 0., 1., (Int_t)(20. / fBinSize), 0., 20.); + TH2F* spectrumHist = new TH2F("LinesTable", "", 100, 0., 1., (Int_t)(20. / fBinSize), 0., 20.); + + 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 + + originalHist->Fill(r, en, (Float_t)flux); + continuumHist->Fill(r, en, (Float_t)flux); + } + debug << "Histograms filled" << endl; + + Int_t peaks = 0; + do { + peaks = 0; + // We just identify pronounced peaks, not smoothed gaussians! + 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; + Float_t flux = data[2]; // flux per cm-2 s-1 keV-1 + + Int_t binR = continuumHist->GetXaxis()->FindBin(r); + Int_t binE = continuumHist->GetYaxis()->FindBin(en); + + Double_t avgFlux = 0; + Int_t n = 0; + for (int e = binE - smearPoints; e <= binE + smearPoints; e++) { + if (e < 1 || (e < binE + excludePoints && e > binE - excludePoints)) continue; + n++; + avgFlux += continuumHist->GetBinContent(binR, e); + } + avgFlux /= n; + + Float_t targetBinFlux = continuumHist->GetBinContent(binR, binE); + Float_t thrFlux = avgFlux + fPeakSigma * TMath::Sqrt(avgFlux); + if (targetBinFlux > 0 && targetBinFlux > thrFlux) { + continuumHist->SetBinContent(binR, binE, avgFlux); + peaks++; + } + } + } while (peaks > 0); + + 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); + + spectrumHist->SetBinContent(n + 1, m + 1, orig - cont); + } + + continuumHist->Rebin2D(1, (Int_t)(0.1 / fBinSize)); // cm-2 s-1 (100eV)-1 + continuumHist->Scale(10); // cm-2 s-1 keV-1 + // It could be over here if we would use directly a TH2D + + 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); + fFluxTable.push_back(hc); + } + + fFluxLines.clear(); + 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( + Form("%s_MonochromeFluxAtEnergy%5.3lf", GetName(), energy), n + 1, n + 1); + fFluxLines[energy] = hm; + } + } + + cout << "Number of peaks identified: " << fFluxLines.size() << endl; +} + +/////////////////////////////////////////////// +/// \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* TRestAxionSolarFlux::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 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* TRestAxionSolarFlux::GetMonochromaticSpectrum() { + if (fMonoHist != nullptr) { + delete fMonoHist; + fMonoHist = nullptr; + } + + fMonoHist = new TH1F("MonochromaticHist", "", 20000, 0, 20); + for (const auto& x : fFluxLines) { + fMonoHist->Fill(x.first, x.second->Integral()); // cm-2 s-1 eV-1 + } + + fMonoHist->SetStats(0); + fMonoHist->GetXaxis()->SetTitle("Energy [keV]"); + fMonoHist->GetXaxis()->SetTitleSize(0.05); + fMonoHist->GetXaxis()->SetLabelSize(0.05); + fMonoHist->GetYaxis()->SetTitle("Flux [cm-2 s-1 eV-1]"); + fMonoHist->GetYaxis()->SetTitleSize(0.05); + fMonoHist->GetYaxis()->SetLabelSize(0.05); + + return fMonoHist; +} + +/////////////////////////////////////////////// +/// \brief It builds a histogram adding the continuum and the monochromatic +/// spectrum component. The flux will be expressed in cm-2 s-1 keV-1. +/// Binned in 1eV steps. +/// +TH1F* TRestAxionSolarFlux::GetTotalSpectrum() { + TH1F* hm = GetMonochromaticSpectrum(); + 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)); + } + } + + for (int n = 0; n < hm->GetNbinsX(); n++) + // 1e-2 is the renormalization from 20000 bins to 200 bins + fTotalHist->SetBinContent(n + 1, fTotalHist->GetBinContent(n + 1) + 100 * hm->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 It builds a histogram using the contents of the .flux file given +/// in the argument. +/// +TH1F* TRestAxionSolarFlux::GetFluxHistogram(string fname, Double_t binSize) { + string fullPathName = SearchFile(fname); + + std::vector> fluxData; + TRestTools::ReadASCIITable(fullPathName, fluxData, 3); + + TH2F* originalHist = + new TH2F(Form("FluxTable", GetName()), "", 100, 0., 1., (Int_t)(20. / binSize), 0., 20.); + + for (const auto& data : fluxData) { + Double_t r = 0.005 + data[0]; + Double_t en = data[1] - 0.005; + Double_t flux = data[2] * binSize; // flux in cm-2 s-1 bin-1 + + originalHist->Fill(r, en, flux); + } + + return (TH1F*)originalHist->ProjectionY(); +} + +/////////////////////////////////////////////// +/// \brief It draws the contents of a .flux file. This method just receives the +/// name of the .flux file and it works stand-alone. +/// +TCanvas* TRestAxionSolarFlux::DrawFluxFile(string fname, Double_t binSize) { + if (fCanvas != nullptr) { + delete fCanvas; + fCanvas = nullptr; + } + fCanvas = new TCanvas("canv", "This is the canvas title", 1400, 1200); + fCanvas->Draw(); + + TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97); + pad1->Draw(); + + fCanvas->cd(); + pad1->cd(); + + GetFluxHistogram(fname, binSize)->Draw("hist"); + + return fCanvas; +} + +/////////////////////////////////////////////// +/// \brief It draws the contents of a .flux file. This method just receives the +/// +TCanvas* TRestAxionSolarFlux::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); + + TH1F* hm = GetMonochromaticSpectrum(); + hm->SetLineColor(kBlack); + hm->Scale(100); // renormalizing per 100eV-1 + + ht->Draw("hist"); + hm->Draw("hist same"); + + pad1->cd(2); + pad1->cd(2)->SetRightMargin(0.09); + pad1->cd(2)->SetLeftMargin(0.15); + pad1->cd(2)->SetBottomMargin(0.15); + + ht->Draw("hist"); + hm->Draw("hist same"); + + return fCanvas; +} + /////////////////////////////////////////////// /// \brief A helper method to initialize the internal class data members with the /// integrated flux for each solar ring. It will be called by TRestAxionSolarFlux::Initialize. @@ -219,6 +657,7 @@ void TRestAxionSolarFlux::LoadMonoChromaticFluxTable() { void TRestAxionSolarFlux::IntegrateSolarFluxes() { fFluxLineIntegrals.clear(); fTotalMonochromaticFlux = 0; + for (const auto& line : fFluxLines) { fTotalMonochromaticFlux += line.second->Integral(); fFluxLineIntegrals.push_back(fTotalMonochromaticFlux); @@ -237,27 +676,13 @@ void TRestAxionSolarFlux::IntegrateSolarFluxes() { fFluxRatio = fTotalMonochromaticFlux / (fTotalContinuumFlux + fTotalMonochromaticFlux); } -/////////////////////////////////////////////// -/// \brief Initialization of TRestAxionSolarFlux metadata members through a RML file -/// -void TRestAxionSolarFlux::InitFromConfigFile() { - debug << "Entering TRestAxionSolarFlux::InitFromConfigFile" << endl; - - fFluxDataFile = GetParameter("fluxDataFile", ""); - fFluxSptFile = GetParameter("fluxSptFile", ""); - fCouplingType = GetParameter("couplingType", "g_ag"); - fCouplingStrength = StringToDouble(GetParameter("couplingStrength", "1.e-10")); - fSeed = StringToInteger(GetParameter("seed", "0")); - - this->Initialize(); -} - /////////////////////////////////////////////// /// \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 TRestAxionSolarFlux::GetRandomEnergyAndRadius() { std::pair result = {0, 0}; + if (!fTablesLoaded) return result; Double_t rnd = fRandom->Rndm(); if (fTotalMonochromaticFlux == 0 || fRandom->Rndm() > fFluxRatio) { // Continuum @@ -305,9 +730,9 @@ void TRestAxionSolarFlux::PrintIntegratedRingFlux() { cout << "Integrated solar flux per solar ring: " << endl; cout << "--------------------------- " << endl; /* -for (int n = 0; n < fFluxPerRadius.size(); n++) + for (int n = 0; n < fFluxPerRadius.size(); n++) cout << "n : " << n << " flux : " << fFluxPerRadius[n] << endl; -cout << endl; + cout << endl; */ } @@ -343,6 +768,8 @@ void TRestAxionSolarFlux::PrintMetadata() { metadata << " - Total continuum flux : " << fTotalContinuumFlux << " cm-2 s-1" << endl; metadata << "--------" << endl; metadata << " - Random seed : " << fSeed << endl; + if (fBinSize > 0) metadata << " - Energy bin size : " << fBinSize * units("eV") << " eV" << endl; + if (fPeakSigma > 0) metadata << " - Peak signal-to-noise in sigmas : " << fPeakSigma << endl; metadata << "++++++++++++++++++" << endl; if (GetVerboseLevel() >= REST_Debug) { @@ -351,3 +778,48 @@ void TRestAxionSolarFlux::PrintMetadata() { PrintIntegratedRingFlux(); } } + +/////////////////////////////////////////////// +/// \brief It will create files with the continuum and spectral flux components to be used +/// in a later ocasion. +/// +void TRestAxionSolarFlux::ExportTables(Bool_t ascii) { + string rootFilename = TRestTools::GetFileNameRoot(fFluxDataFile); + + string path = REST_USER_PATH + "/export/"; + + if (!TRestTools::fileExists(path)) { + cout << "Creating path: " << path << endl; + system(("mkdir -p " + path).c_str()); + } + + 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); + } + + if (fFluxLines.size() > 0) { + std::vector> table; + for (const auto& x : fFluxLines) { + 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)); + + table.push_back(row); + } + + TRestTools::TransposeTable(table); + + TRestTools::ExportASCIITable(path + "/" + rootFilename + ".spt", table); + } +}