GSL Integration implementation, TRestAxionMagneticField ReMap feature and HOT bug fix.

inc/TRestAxionField.h

@@ -24,16 +24,17 @@
 #define _TRestAxionField
 
 #include "TRestAxionBufferGas.h"
+#include "TRestAxionMagneticField.h"
 
 //! A basic class to define analytical axion-photon conversion calculations for axion helioscopes
 class TRestAxionField : public TObject {
    private:
     Bool_t fDebug = false;  //!
 
-    /// The magnetic field in Teslas
+    /// The magnetic field in Teslas (used for constant field formulas)
     Double_t fBmag = 2.5;
 
-    /// The coherence lenght (in mm) where the magnetic field is defined
+    /// The coherence lenght (in mm) where the magnetic field is defined (for constant field)
     Double_t fLcoh = 10000;
 
     /// The energy of the axion in keV
@@ -42,7 +43,16 @@ class TRestAxionField : public TObject {
     void Initialize();
 
     /// A pointer to the buffer gas definition
-    TRestAxionBufferGas* fBufferGas = NULL;  //!
+    TRestAxionBufferGas* fBufferGas = nullptr;  //!
+
+    /// A pointer to the magnetic field definition
+    TRestAxionMagneticField* fMagneticField = nullptr;  //!
+
+    std::pair<Double_t, Double_t> ComputeOffResonanceIntegral(Double_t q, Double_t Gamma, Double_t accuracy,
+                                                              Int_t num_intervals, Int_t qawo_levels);
+
+    std::pair<Double_t, Double_t> ComputeResonanceIntegral(Double_t Gamma, Double_t accuracy,
+                                                           Int_t num_intervals);
 
    public:
     void SetMagneticField(Double_t b) { fBmag = b; }
@@ -62,6 +72,9 @@ class TRestAxionField : public TObject {
     /// 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 AssignMagneticField(TRestAxionMagneticField* mField) { fMagneticField = mField; }
+
     /// It assigns a gas buffer medium to the calculation
     void SetBufferGas(TRestAxionBufferGas* buffGas) { fBufferGas = buffGas; }
 
@@ -78,6 +91,21 @@ class TRestAxionField : public TObject {
     Double_t GammaTransmissionProbability(std::vector<Double_t> Bmag, Double_t deltaL, Double_t Ea,
                                           Double_t ma, Double_t mg = 0, Double_t absLength = 0);
 
+    std::pair<Double_t, Double_t> GammaTransmissionFieldMapProbability(Double_t Ea, Double_t ma,
+                                                                       Double_t accuracy = 1.e-1,
+                                                                       Int_t num_intervals = 100,
+                                                                       Int_t qawo_levels = 20);
+
+    /// Integrand used for axion-photon probability integration
+    static double Integrand(double x, void* params) {
+        auto* data = reinterpret_cast<std::pair<TRestAxionMagneticField*, double>*>(params);
+
+        TRestAxionMagneticField* field = data->first;
+        double gamma = data->second;
+
+        return exp(0.5 * gamma * x) * field->GetTransversalComponentInParametricTrack(x);
+    }
+
     Double_t GammaTransmissionFWHM(Double_t step = 0.00001);
 
     std::vector<std::pair<Double_t, Double_t>> GetMassDensityScanning(std::string gasName = "He",

inc/TRestAxionFieldPropagationProcess.h

@@ -34,12 +34,21 @@
 //! A process to introduce the magnetic field profile integration along the track
 class TRestAxionFieldPropagationProcess : public TRestAxionEventProcess {
    private:
-    /// The differential length in mm used for the field integration
-    Double_t fIntegrationStep = 50;  //<
-
     /// The additional length in mm that the converted photon propagates without magnetic field
     Double_t fBufferGasAdditionalLength = 0;  //<
 
+    /// The tolerance or accuracy used inside the GSL integrator
+    Double_t fAccuracy = 0.1;
+
+    /// Number of intervales used by the GSL integrator
+    Int_t fNumIntervals = 100;
+
+    /// Number of levels used by the GSL integrator to parameterize the cosine integral
+    Int_t fQawoLevels = 20;
+
+    /// It will re-size the cells in the magnetic field (affecting the time required for integral computation)
+    TVector3 fReMap = TVector3(30, 30, 100);
+
     /// A pointer to the magnetic field description stored in TRestRun
     TRestAxionMagneticField* fMagneticField = nullptr;  //!
 
@@ -63,7 +72,14 @@ class TRestAxionFieldPropagationProcess : public TRestAxionEventProcess {
     void PrintMetadata() override {
         BeginPrintProcess();
 
-        RESTMetadata << "Integration step length : " << fIntegrationStep << " mm" << RESTendl;
+        RESTMetadata << "Integration parameters" << RESTendl;
+        RESTMetadata << "- Integration accuracy/tolerance : " << fAccuracy << RESTendl;
+        RESTMetadata << "- Max number of integration intervals : " << fNumIntervals << RESTendl;
+        RESTMetadata << "- Number of QAWO levels : " << fQawoLevels << RESTendl;
+        RESTMetadata << " " << RESTendl;
+        RESTMetadata << "Field re-mapping size : (" << fReMap.X() << ", " << fReMap.Y() << ", " << fReMap.Z()
+                     << ")" << RESTendl;
+        RESTMetadata << " " << RESTendl;
         RESTMetadata << "Buffer gas additional length : " << fBufferGasAdditionalLength * units("m") << " m"
                      << RESTendl;
 

inc/TRestAxionMagneticField.h

@@ -58,7 +58,8 @@ class TRestAxionMagneticField : public TRestMetadata {
     /// A constant field component that will be added to the field map
     std::vector<TVector3> fConstantField;  //<
 
-    /// The size of a grid element from the mesh in mm
+    /// The size of a grid element from the mesh in mm. Initially, it must be the same as the binary input
+    /// data
     std::vector<TVector3> fMeshSize;  //<
 
     /// The type of the mesh used (default is cylindrical)
@@ -67,9 +68,21 @@ class TRestAxionMagneticField : public TRestMetadata {
     /// A vector to store the maximum bounding box values
     std::vector<TVector3> fBoundMax;  //<
 
+    /// A vector that defines the new mesh cell volume. It will re-scale the original fMeshSize.
+    TVector3 fReMap = TVector3(0, 0, 0);  //<
+
     /// A magnetic field volume structure to store field data and mesh.
     std::vector<MagneticFieldVolume> fMagneticFieldVolumes;  //!
 
+    /// The start track position used to parameterize the field along a track
+    TVector3 fTrackStart = TVector3(0, 0, 0);  //!
+
+    /// The track direction used to parameterize the field along a track
+    TVector3 fTrackDirection = TVector3(0, 0, 0);  //!
+
+    /// The total length of the track which defines the limit for field parameterization
+    Double_t fTrackLength = 0;  //!
+
     /// A helper histogram to plot the field
     TH2D* fHisto;  //!
 
@@ -102,6 +115,14 @@ class TRestAxionMagneticField : public TRestMetadata {
    public:
     void LoadMagneticVolumes();
 
+    void ReMap(const size_t& n, const TVector3& newMapSize);
+
+    void SetTrack(const TVector3& position, const TVector3& direction);
+
+    Double_t GetTrackLength() const { return fTrackLength; }
+    TVector3 GetTrackStart() const { return fTrackStart; }
+    TVector3 GetTrackDirection() const { return fTrackDirection; }
+
     /// It returns true if no magnetic field map was loaded for that volume
     Bool_t IsFieldConstant(Int_t id) {
         if (GetMagneticVolume(id)) return GetMagneticVolume(id)->field.size() == 0;
@@ -135,6 +156,7 @@ class TRestAxionMagneticField : public TRestMetadata {
     TVector3 GetVolumeCenter(Int_t id);
 
     Double_t GetTransversalComponent(TVector3 position, TVector3 direction);
+    Double_t GetTransversalComponentInParametricTrack(Double_t x);
 
     std::vector<Double_t> GetTransversalComponentAlongPath(TVector3 from, TVector3 to, Double_t dl = 1.,
                                                            Int_t Nmax = 0);

macros/REST_Axion_FieldIntegrationTests.C

@@ -0,0 +1,50 @@
+#include "TCanvas.h"
+#include "TGraph.h"
+#include "TLatex.h"
+#include "TLegend.h"
+#include "TLine.h"
+#include "TRestAxionBufferGas.h"
+#include "TRestAxionField.h"
+//*******************************************************************************************************
+//*** Description: This script will launch the integration of the axion-field with given parameters.
+//*** It allows to test different magnetic field cell sizes, for a given mass that can be off-resonance
+//*** for dm different from zero, and a given maximum tolerance or error for the integration routine.
+//***
+//*** The macro sets the TRestAxionField under debug mode to print the different results on screen.
+//***
+//*** --------------
+//*** Usage: restManager FieldIntegrationTests [sX=10] [sX=10] [sZ=10] [dm=0.01] [tolerance=0.1] [Ea=4.2]
+//*** --------------
+//***
+//*** Author: Javier Galan
+//*******************************************************************************************************
+int REST_Axion_FieldIntegrationTests(Double_t sX = 10, Double_t sY = 10, Double_t sZ = 50, Double_t dm = 0.01,
+                                     Double_t tolerance = 0.1, Double_t Ea = 4.2) {
+    /// Setting up magnetic field and track to evaluate
+    TRestAxionMagneticField magneticField("fields.rml", "babyIAXO_HD");
+
+    for (size_t n = 0; n < magneticField.GetNumberOfVolumes(); n++)
+        magneticField.ReMap(n, TVector3(sX, sY, sZ));
+    magneticField.SetTrack(TVector3(-110, -110, -11000), TVector3(0.01, 0.01, 1));
+    magneticField.PrintMetadata();
+
+    /// Setting up the gas
+    TRestAxionBufferGas* gas = new TRestAxionBufferGas();
+    gas->SetGasDensity("He", 3.267069078540181e-10);
+    gas->PrintMetadata();
+    Double_t ma = gas->GetPhotonMass(Ea);
+
+    /// Setting up the axion-field and assign gas and magnetic field.
+    TRestAxionField* ax = new TRestAxionField();
+    ax->AssignBufferGas(gas);
+    ax->AssignMagneticField(&magneticField);
+    ax->SetAxionEnergy(Ea);
+
+    /// Enable debugging mode in axion-field
+    ax->SetDebug(true);
+
+    std::pair<double, double> prob =
+        ax->GammaTransmissionFieldMapProbability(Ea, ma - dm, tolerance, 100, 25);
+
+    return 0;
+}

pipeline/ray-tracing/axion-field/Validate.C

@@ -1,10 +1,9 @@
 #include <TH1D.h>
 #include <TRestRun.h>
 
-Double_t probTolerance = 1.e-21;
-Double_t fieldTolerance = 0.01;
+Double_t probTolerance = 1.e-22;
 
-Int_t Validate(Double_t prob = 5.18573e-19, Double_t fieldAverage = 1.5764) {
+Int_t Validate(Double_t prob = 2.34295e-21) {
     TRestRun* run = new TRestRun("AxionPhotonProbability.root");
 
     if (run->GetEntries() != 100) {
@@ -34,25 +33,30 @@ Int_t Validate(Double_t prob = 5.18573e-19, Double_t fieldAverage = 1.5764) {
         return 3;
     }
 
-    run->GetAnalysisTree()->Draw("axionPhoton_fieldAverage", "axionPhoton_fieldAverage");
+    /* We do not add the field average observable anymore at TRestAxionFieldPropagationProcess
+     * We would need to create a new process that does this, since it is computationally
+     * not negligible
+     *
+run->GetAnalysisTree()->Draw("axionPhoton_fieldAverage", "axionPhoton_fieldAverage");
 
-    TH1D* h2 = (TH1D*)run->GetAnalysisTree()->GetHistogram();
-    if (h2 == nullptr) {
-        std::cout << "Problems generating field average histogram" << std::endl;
-        return 4;
-    }
+TH1D* h2 = (TH1D*)run->GetAnalysisTree()->GetHistogram();
+if (h2 == nullptr) {
+    std::cout << "Problems generating field average histogram" << std::endl;
+    return 4;
+}
 
-    Double_t field = h2->Integral() / run->GetEntries();
-    std::cout << "Average magnetic field: " << field << " T" << std::endl;
-    std::cout << "Expected field: " << fieldAverage << std::endl;
-    std::cout << "Tolerance: " << fieldTolerance << std::endl;
-    std::cout << "Low : " << fieldAverage - fieldTolerance << " high: " << fieldAverage + fieldTolerance
-              << std::endl;
+Double_t field = h2->Integral() / run->GetEntries();
+std::cout << "Average magnetic field: " << field << " T" << std::endl;
+std::cout << "Expected field: " << fieldAverage << std::endl;
+std::cout << "Tolerance: " << fieldTolerance << std::endl;
+std::cout << "Low : " << fieldAverage - fieldTolerance << " high: " << fieldAverage + fieldTolerance
+          << std::endl;
 
-    if (field < fieldAverage - fieldTolerance || field > fieldAverage + fieldTolerance) {
-        std::cout << "Wrong average field!" << std::endl;
-        return 5;
-    }
+if (field < fieldAverage - fieldTolerance || field > fieldAverage + fieldTolerance) {
+    std::cout << "Wrong average field!" << std::endl;
+    return 5;
+}
+    */
 
     delete run;
 

pipeline/ray-tracing/axion-field/photonConversion.rml

@@ -26,7 +26,13 @@
         </addProcess>
         <addProcess type="TRestAxionTransportProcess" zPosition="-5m" name="to5m"/>
         <addProcess type="TRestAxionDeviationProcess" name="dev" devAngle="2deg" seed="137"/>
-        <addProcess type="TRestAxionFieldPropagationProcess" name="axionPhoton" integrationStep="5cm" bufferGasAdditionalLength="3m" observables="all" verboseLevel="info"/>
+
+        <addProcess type="TRestAxionFieldPropagationProcess" name="axionPhoton" bufferGasAdditionalLength="3m" observables="all" verboseLevel="info">
+			<parameter name="accuracy" value="0.1"/>
+			<parameter name="numIntervals" value="100" />
+			<parameter name="qawoLevels" value="20"/>
+		</addProcess>
+
         <addProcess type="TRestAxionAnalysisProcess" name="at5m" observables="all"/>
     </TRestProcessRunner>
     <addTask command="EventProcess-&gt;RunProcess()" value="ON"/>

pipeline/ray-tracing/axion-field/plots.rml

@@ -1,6 +1,7 @@
 <?xml version="1.0" encoding="UTF-8" standalone="no"?>
 <TRestManager name="SpecPlot" title="Example" verboseLevel="info">
     <TRestRun name="DummyRun"/>
+	<!--- Some of these observables are only available up to release v2.4.1. After release v2.4.2 these observables have changed! -->
     <TRestAnalysisPlot name="restplot" title="Optics Plots" verboseLevel="warning">
         <parameter name="previewPlot" value="false"/>
         <canvas size="(3600,2400)" divide="(3,3)" save="axionFieldPlots.png"/>

pipeline/ray-tracing/full-chain/ValidateChain.C

@@ -46,7 +46,7 @@ Int_t ValidateChain(std::string fname) {
     Double_t integral2 = g->Integral() / run->GetEntries();
     std::cout << "Average window transmission : " << integral2 << std::endl;
 
-    if (integral < 9.71489e-24 || integral > 9.91489e-24) {
+    if (integral > 2.01953e-22 || integral < 1.99953e-22) {
         std::cout << "Axion-photon probability is not within the expected range!" << std::endl;
         return 3;
     }