Hidden Photon flux integration

.gitignore

@@ -1,3 +1,4 @@
 build
 *.root
 **/.DS_Store
+*.json

.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

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;  //!

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

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;  //!
 

inc/TRestAxionField.h → 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

inc/TRestAxionSolarFlux.h

@@ -24,8 +24,8 @@
 #define _TRestAxionSolarFlux
 
 #include <TCanvas.h>
-#include <TH1F.h>
-#include <TH2F.h>
+#include <TH1D.h>
+#include <TH2D.h>
 #include <TRandom3.h>
 #include <TRestMetadata.h>
 
@@ -38,6 +38,9 @@ class TRestAxionSolarFlux : public TRestMetadata {
     /// Axion coupling strength
     Double_t fCouplingStrength = 0;  //<
 
+    /// Mass parameter
+    Double_t fMass = 0;  //!
+
     /// Seed used in random generator
     Int_t fSeed = 0;  //<
 
@@ -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<Double_t, Double_t> GetRandomEnergyAndRadius(TVector2 eRange = TVector2(-1, -1),
-                                                                   Double_t mass = 0) = 0;
+    virtual std::pair<Double_t, Double_t> 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();

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 <TRestAxionSolarFlux.h>
+#include <TRestAxionSolarModel.h>
+
+//! 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<TH1D*> fFluxTable;  //!
+
+    /// The tabulated solar flux continuum spectra TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius
+    std::vector<TH1D*> fContinuumTable;  //!
+
+    /// The tabulated resonance width TH1D(200,0,20)keV in eV2 versus solar radius
+    std::vector<TH1D*> fWidthTable;  //!
+
+    /// The solar plasma frequency vector in eV versus solar radius
+    std::vector<TH1D*> fPlasmaFreqTable;  //!
+
+    /// The total solar flux TH1D(200,0,20)keV in cm-2 s-1 keV-1 versus solar radius
+    std::vector<TH1D*> fFullFluxTable;  //!
+
+    /// Accumulative integrated solar flux for each solar ring for continuum spectrum (renormalized to unity)
+    std::vector<Double_t> 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<Double_t, Double_t> 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