From 2faa963ffb26e7f6462ef396a1ff4dd41ba4ec17 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 19 Mar 2024 12:50:07 +0100 Subject: [PATCH 1/7] scripts/xlsBikovskiyToBin.py updating to a more modern python style --- scripts/xlsBikovskiyToBin.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/scripts/xlsBikovskiyToBin.py b/scripts/xlsBikovskiyToBin.py index a8e8dd23..07c89cff 100755 --- a/scripts/xlsBikovskiyToBin.py +++ b/scripts/xlsBikovskiyToBin.py @@ -22,13 +22,14 @@ print("Starting to read") # loading the data into a matrix (xyzBdata) -file = r"../../../data/magneticField/Bykovskiy_201906.xls" +file = r"../data/magneticField/Bykovskiy_202004.xls" df = pd.read_excel(file) print(df[1:5]) print("Translating to matrix") -xyzBdata = df.as_matrix(columns=df.columns[0:]) +#xyzBdata = df.values(columns=df.columns[0:]) +xyzBdata = df[df.columns[0:]].values print(xyzBdata[0][0:6]) From 20d251bf6c40046f5d141d14e515fec929f1c24d Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 19 Mar 2024 12:50:46 +0100 Subject: [PATCH 2/7] Uploading xlsToBin scripts --- scripts/xlsBikovskiyToBin.py | 1 - scripts/xlsMentinkToBin.py | 153 +++++++++++++++++++++++++++++++++++ 2 files changed, 153 insertions(+), 1 deletion(-) create mode 100755 scripts/xlsMentinkToBin.py diff --git a/scripts/xlsBikovskiyToBin.py b/scripts/xlsBikovskiyToBin.py index 07c89cff..09a3a99d 100755 --- a/scripts/xlsBikovskiyToBin.py +++ b/scripts/xlsBikovskiyToBin.py @@ -28,7 +28,6 @@ print(df[1:5]) print("Translating to matrix") -#xyzBdata = df.values(columns=df.columns[0:]) xyzBdata = df[df.columns[0:]].values print(xyzBdata[0][0:6]) diff --git a/scripts/xlsMentinkToBin.py b/scripts/xlsMentinkToBin.py new file mode 100755 index 00000000..26ca25fd --- /dev/null +++ b/scripts/xlsMentinkToBin.py @@ -0,0 +1,153 @@ +#!/usr/bin/python + +# This basic script serves to create a plain binary file from the XLS files provided by Bykovsky +# for the Baby-IAXO magnet. We will only upload the resulting bin file, if the original XLS file +# would be required at some point, do not hesitate to contact me at javier.galan@unizar.es. +# +# I upload this script in case it is necessary to convert future field maps. Some revision will be +# needed in order to adapt to other input/output schemes. +# + +import sys +import struct +from pandas import DataFrame, read_csv +import matplotlib.pyplot as plt +import pandas as pd + +# Field map is provided only for the top magnetic bore. We will recenter the field map. +# Since this is a condition for TRestAxionMagneticField. +xCenter = 0 +yCenter = 0.8 +zCenter = 0 + +print("Starting to read") +# loading the data into a matrix (xyzBdata) +file = r"../data/magneticField/Mentink_202401.txt" +df = pd.read_excel(file) + +print(df[1:5]) + +print("Translating to matrix") +#xyzBdata = df.values(columns=df.columns[0:]) +xyzBdata = df[df.columns[0:]].values + +print(xyzBdata[0][0:6]) + +print(len(xyzBdata)) + +# Printing out some info about data +xMax = -100000 +xMin = 100000 +yMax = -100000 +yMin = 100000 +zMax = -100000 +zMin = 100000 + +xBMax = -100000 +xBMin = 100000 +yBMax = -100000 +yBMin = 100000 +zBMax = -100000 +zBMin = 100000 + +f = open("output.txt", "wt") + +for x in xyzBdata: + f.write( + str(x[0]) + + "\t" + + str(x[1]) + + "\t" + + str(x[2]) + + "\t" + + str(x[3]) + + "\t" + + str(x[4]) + + "\t" + + str(x[5]) + + "\n" + ) + if xMax < x[0]: + xMax = x[0] + if yMax < x[1]: + yMax = x[1] + if zMax < x[2]: + zMax = x[2] + + if xMin > x[0]: + xMin = x[0] + if yMin > x[1]: + yMin = x[1] + if zMin > x[2]: + zMin = x[2] + + if xBMax < x[3]: + xBMax = x[3] + if yBMax < x[4]: + yBMax = x[4] + if zBMax < x[5]: + zBMax = x[5] + + if xBMin > x[3]: + xBMin = x[3] + if yBMin > x[4]: + yBMin = x[4] + if zBMin > x[5]: + zBMin = x[5] + +f.close() + +print("X max : " + str(xMax)) +print("Y max : " + str(yMax)) +print("Z max : " + str(zMax) + "\n") + +print("X min : " + str(xMin)) +print("Y min : " + str(yMin)) +print("Z min : " + str(zMin) + "\n") + +print("BX max : " + str(xBMax)) +print("BY max : " + str(yBMax)) +print("BZ max : " + str(zBMax) + "\n") + +print("BX min : " + str(xBMin)) +print("BY min : " + str(yBMin)) +print("BZ min : " + str(zBMin) + "\n") + +# Creating the binary file +fbin = open("output.bin", "wb") +count = 0 +for x in xyzBdata: + # We recenter the volume and redefine axis (x becomes z, y becomes x and z becomes y) + # XLS file distances are expressed in m. We translate to mm. + y = [ + 1000 * (x[0] - xCenter), + 1000 * (x[1] - yCenter), + 1000 * (x[2] - zCenter), + x[3], + x[4], + x[5], + ] + + # Baby-IAXO is symmetric respect to z (direction along axis) and x (direction parallel to the ground). + # I.e. x and y in the previous axis definition. + # The y-axis symmetry (in the new axis definition) affects the second bore that is not described in this map. + # We apply the corresponding symmetries to define the full map of one magnetic bore in the output binary file. + + count = count + 1 + fbin.write(struct.pack("<%df" % len(y), *y)) + if count < 6: + print(len(y)) + print(x[0:6]) + print(y[0:6]) + # The original file was only missing the x-axis, when we change the sign of x-axis we must change the sign of Bx. + if y[0] < 0: + y[0] = -y[0] + y[3] = -y[3] + count = count + 1 + fbin.write(struct.pack("<%df" % len(y), *y)) + if count < 6: + print(len(y)) + print(x[0:6]) + print(y[0:6]) +fbin.close() +print("Lines writen : " + str(count)) From 8c0c7200de8f0b27425ca563c381484c2822cdb1 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 19 Mar 2024 12:58:37 +0100 Subject: [PATCH 3/7] TRestAxionBufferGas. Adding extra example in documentation --- src/TRestAxionBufferGas.cxx | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/TRestAxionBufferGas.cxx b/src/TRestAxionBufferGas.cxx index 5631950e..23acf354 100644 --- a/src/TRestAxionBufferGas.cxx +++ b/src/TRestAxionBufferGas.cxx @@ -51,6 +51,7 @@ /// gas->SetGasMixture( "He+Xe", "2.6e-6g/dm^3+5.6mg/m^3" ); /// \endcode /// +/// /// The corresponding RML section for initialization through a configuration /// file would be as follows. /// @@ -61,6 +62,15 @@ /// /// \endcode /// +/// Example 3: +/// +/// We may also use a definition found inside the official REST-for-Physics +/// data directory. +/// +/// \code +/// TRestAxionBufferGas *gas = new TRestAxionBufferGas("bufferGases.rml", "helium"); +/// \endcode +/// ///-------------------------------------------------------------------------- /// /// RESTsoft - Software for Rare Event Searches with TPCs From 5dc99bc4548c6364137fed5c1199bfc337111655 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 19 Mar 2024 19:24:53 +0100 Subject: [PATCH 4/7] TRestAxionField. Adding documentation for the integration of the field --- src/TRestAxionField.cxx | 179 +++++++++++++++++++++++++++++++++++----- 1 file changed, 160 insertions(+), 19 deletions(-) diff --git a/src/TRestAxionField.cxx b/src/TRestAxionField.cxx index 0883e6ce..a34ef0f5 100644 --- a/src/TRestAxionField.cxx +++ b/src/TRestAxionField.cxx @@ -21,13 +21,152 @@ *************************************************************************/ ////////////////////////////////////////////////////////////////////////// +//// /// TRestAxionField 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. +/// and determine the probability of the particle being in a photon +/// state after propagating inside a magnetic field. +/// +/// The class allows to assign a buffer gas using a TRestAxionBufferGas and +/// a magnetic field map loaded through TRestAxionMagneticField. If these +/// objects have been assigned, the methods implemented in this class +/// may consider those - if needed - inside the calculation. +/// +/// In practice this class provides different versions of the method +/// GammaTransmissionProbability that allows to calculate the axion-photon +/// probability using different strategies. +/// +/// \section constant Calculating axion-photon probability in a constant field. +/// +/// ### In vacuum +/// +/// For calculations inside a constant magnetic field one may simply +/// invoke the following code which will launch the calculation in *vacuum*, +/// for a field of 2T, a coherence length of 10000mm, and an axion energy +/// and mass of 4.2keV and 0.1eV, respectively. +/// +/// \code +/// TRestAxionField field; +/// field->GammaTransmissionProbability( 2, 10000, 4.2, 0.1 ); +/// \endcode +/// +/// It is possible to reduce the number of arguments we are giving to this +/// function by assigning few data member values present inside TRestAxionField +/// as follows: +/// +/// \code +/// TRestAxionField field; +/// field.SetMagneticField(2); +/// field.SetCoherenceLength(10000); +/// field.SetAxionEnergy(4.2); +/// field.GammaTransmissionProbability( 0.1 ); +/// \endcode +/// +/// Indeed, the complete version with all arguments will also update the +/// data member values, so that next calls to this method will use the same +/// magnetic field, coherence length and energy. In the following example +/// the second call to GammaTransmissionProbability will use again the same +/// values for magnetic field, coherence length and energy from the +/// preceding call. +/// +/// \code +/// TRestAxionField field; +/// field.GammaTransmissionProbability( 2, 10000, 4.2, 0.1 ); +/// field.GammaTransmissionProbability( 0.01 ); +/// \endcode +/// +/// ### In a buffer gas medium +/// +/// The axion-photon probability can also be calculated in the sin of a +/// gaseous medium. For that we need to assign a buffer gas instance containing +/// the relevant gas properties in the form of a TRestAxionBufferGas. +/// +/// The following piece of code will assign the gas properties to the field +/// calculation as defined inside the `bufferGases.rml` official file. The gas +/// provides the absorption and equivalent photon mass as a funtion of the +/// energy. +/// +/// \code +/// TRestAxionField field; +/// TRestAxionBufferGas gas("bufferGases.rml", "helium"); +/// field.AssignBufferGas( &gas ); +/// field.SetMagneticField(2); +/// field.SetCoherenceLength(10000); +/// field.SetAxionEnergy(4.2); +/// field.GammaTransmissionProbability( 0.1 ); +/// \endcode +/// +/// \section non-constant Calculating axion-photon probability in an unhomogeneous field. +/// +/// There are two main strategies presently implemented inside this class to integrate +/// the axion field along an unhomogeneous magnetic field. The first one uses a simple +/// trapezoidal integration method with the advantage that the computing time is controlled +/// by length of the integration step, but no error calculation is given. The second +/// method advantage is that in principle it allows to fix the accuracy of our calculation +/// and an error of the integration is returned by the GSL integration method used. The +/// main disadvantage is that by fixing the accuracy we lose control on the computation +/// time required to reach such accuracy. Here we describe both methods. +/// +/// ### 1. Using external field profile +/// +/// Using a external field profile requires to prepare the field profile and pass it +/// along to TRestAxionField::GammaTransmissionProbability method performing the +/// integration. The profile will be placed in a `std::vector` that describes +/// the magnetic field profile with a fixed length step. +/// The method TRestAxionMagneticField::GetTransversalComponentAlongPath can be used +/// to generate such profile and calculate the field using the first method, as in +/// the following example: +/// +/// \code +/// Double_t deltaL = 100; // mm +/// Double_t Ea = 4.2; // keV +/// Double_t ma = 0.001; // eV +/// +/// // Recovering the B-field profile along the axis +/// TRestAxionMagneticField mField("fields.rml", "babyIAXO_2024"); +/// std::vector bField = mField.GetTransversalComponentAlongPath( TVector3(0,0,-6000), +/// TVector3(0,0,6000), deltaL ); +/// +/// TRestAxionField aField; +/// aField.GammaTransmissionProbability( bField, deltaL, Ea, ma ); +/// \endcode +/// +/// In a similar way to the constant field case, it would suffice to assign a buffer gas to +/// the TRestAxionField instance to have the calculation performed in a buffer gas. +/// +/// ### 2. Using a field map instance +/// +/// A new method has been added in which the integration method is given a pointer to the +/// function that will allow evaluation of the magnetic field along the B-profile +/// parameterized track. Thus, it is the integration method the one requesting the points +/// on the parametric curve that will be necessary to perform the integration. The required +/// field evaluation is encoded inside TRestAxionField::Integrand which accesses the +/// magnetic field instance TRestAxionField::fMagneticField in order to recover a parametric +/// track along the path. +/// +/// The method name using this technique requires thus direct access to a +/// TRestAxionMagneticField instance that must be assigned beforehand. For computing +/// the integral following this strategy one needs to write the following code where we +/// consider the most generic case including a gas buffer medium. +/// +/// \code +/// TRestAxionField aField; +/// +/// TRestAxionBufferGas gas("bufferGases.rml", "helium"); +/// aField.AssignBufferGas( &gas ); +/// +/// TRestAxionMagneticField mField("fields.rml", "babyIAXO_2024"); +/// aField.AssignMagneticField( &mField ); +/// +/// Double_t ma = 0.001; +/// Double_t Ea = 4.2; +/// Double_t accuracy = 0.1; +/// // The track that will be integrated must be predefined before inside the magnetic field class +/// mField.SetTrack( TVector3(0,0,-6000), TVector3(0,0,1) ); +/// std::pair prob = aField.GammaTransmissionFieldMapProbability(Ea, ma, accuracy); +/// \endcode +/// +/// Where the later call will return a `std::pair` with the probability and its +/// corresponding error. /// ///-------------------------------------------------------------------------- /// @@ -584,8 +723,8 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t Bmag, Double_t Lco /// Default value for the parameters are: Bmag = 2.5 T, Lcoh = 10000 mm, Ea = 4 keV /// eV, /// -/// 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. +/// 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 maxMass = 10; // 10eV is the maximum mass (exit condition) @@ -599,9 +738,9 @@ Double_t TRestAxionField::GammaTransmissionFWHM(Double_t step) { while (Pmax / 2 < GammaTransmissionProbability(scanMass)) { scanMass += step; if (scanMass > maxMass) { - RESTError - << "TRestAxionField::GammaTransmissionProbability. Something went wrong when calculating FWHM" - << RESTendl; + RESTError << "TRestAxionField::GammaTransmissionProbability. Something went wrong when " + "calculating FWHM" + << RESTendl; return maxMass; } } @@ -615,23 +754,25 @@ Double_t TRestAxionField::GammaTransmissionFWHM(Double_t step) { } /////////////////////////////////////////////// -/// \brief This method determines the proper densities to be used in an axion helioscope experiment in order -/// to achieve a continuous axion mass scan. +/// \brief This method determines the proper densities to be used in an axion helioscope experiment in +/// order to achieve a continuous axion mass scan. /// /// The first scanning density is placed where the axion-photon vacuum probability reaches half the value, /// `P_ag(max)/2`. Once the first density, or step, has been obtained, the method calculates the FWHM -/// resonance probability for each density/mass and moves the next scanning axion mass by a step of amplitude -/// `FWHM/factor`, where the factor is a value that moves from 2 to 1 as the mass increases, and it falls down -/// with a velocity related with the argument `rampDown`, where the factor follows this formula: +/// resonance probability for each density/mass and moves the next scanning axion mass by a step of +/// amplitude `FWHM/factor`, where the factor is a value that moves from 2 to 1 as the mass increases, and +/// it falls down with a velocity related with the argument `rampDown`, where the factor follows this +/// formula: /// /// factor = TMath::Exp( - ma*rampDown ) + 1; /// -/// The method stops when the axion mass is bigger than the maximum axion mass given as an argument, `maMax`. +/// The method stops when the axion mass is bigger than the maximum axion mass given as an argument, +/// `maMax`. /// /// Default arguments: gasName="He", ma_max=0.15 eV, rampDown=5 /// -/// \return It returns a vector of pair with the values for the scan, the first one is the axion mass and the -/// second one is the density. +/// \return It returns a vector of pair with the values for the scan, the first one is the axion mass and +/// the second one is the density. /// /// For additional info see PR: https://github.com/rest-for-physics/axionlib/pull/78 /// From fe7e72f29e34e1b90c822691acc7b32c4557d84f Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 19 Mar 2024 19:25:30 +0100 Subject: [PATCH 5/7] TRestAxionMagneticField. Adding documentation for ReMap --- src/TRestAxionMagneticField.cxx | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/src/TRestAxionMagneticField.cxx b/src/TRestAxionMagneticField.cxx index fe6ae7fd..763819de 100644 --- a/src/TRestAxionMagneticField.cxx +++ b/src/TRestAxionMagneticField.cxx @@ -161,8 +161,8 @@ /// /// There are also geometric functions that allow to identify the boundaries of /// the magnetic volume. A test particle will penetrate in the bounding box and identify -/// the moment where the field changes to a value different from (0,0,0) in order -/// to identify the entrance and exit point. +/// the moment where the field changes to a value different from (0,0,0) in order("fields.rml", +/// "babyIAXO_2024"); to identify the entrance and exit point. /// /// The following code will return two TVector3 with the magnetic volume entrance /// and exit coordinates, for a test particle being placed at x=10cm, y=10cm and @@ -177,6 +177,22 @@ /// In the other hand, TRestAxionMagneticField::GetVolumeBoundaries will return the /// bounding box containing that magnetic field. /// +/// ### Remapping the field +/// +/// It is possible to remap the magnetic field in order to increase the size of the grid +/// cells in order to gain in performance. We may only increase the size of the cells +/// by a multiple of the original mesh size. +/// +/// Given the original mesh granularity by `10mmx10mmx50mm`, we might increase the +/// granularity of the first volume, index `0` to `50mmx50mmx200mm` as follows: +/// +/// \code +/// TRestAxionMagneticField *mag = new TRestAxionMagneticField("fields.rml", "babyIAXO" ); +/// mag->PrintMetadata(); +/// mag->ReMap( 0, TVector3(50, 50, 200) ); +/// mag->PrintMetadata(); +/// \endcode +/// /// ### Visualizing the magnetic field /// /// You may visualize the magnetic field profile along tracks towards a vanishing point From b680476d252b3eca6e42d71724952b7fc7bafa56 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 19 Mar 2024 18:27:31 +0000 Subject: [PATCH 6/7] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- scripts/xlsMentinkToBin.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/xlsMentinkToBin.py b/scripts/xlsMentinkToBin.py index 26ca25fd..ee8de25c 100755 --- a/scripts/xlsMentinkToBin.py +++ b/scripts/xlsMentinkToBin.py @@ -28,7 +28,7 @@ print(df[1:5]) print("Translating to matrix") -#xyzBdata = df.values(columns=df.columns[0:]) +# xyzBdata = df.values(columns=df.columns[0:]) xyzBdata = df[df.columns[0:]].values print(xyzBdata[0][0:6]) From 867a4aff1b3ae3cd75fbce94727d5d6e50741982 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Wed, 20 Mar 2024 11:07:36 +0100 Subject: [PATCH 7/7] TRestAxionField updating documentation --- src/TRestAxionField.cxx | 35 ++++++++++++++++++++++++++++------- 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/src/TRestAxionField.cxx b/src/TRestAxionField.cxx index a34ef0f5..a0706be4 100644 --- a/src/TRestAxionField.cxx +++ b/src/TRestAxionField.cxx @@ -35,9 +35,9 @@ /// GammaTransmissionProbability that allows to calculate the axion-photon /// probability using different strategies. /// -/// \section constant Calculating axion-photon probability in a constant field. +/// ## Calculating axion-photon probability in a constant field. /// -/// ### In vacuum +/// ### 1. In vacuum /// /// For calculations inside a constant magnetic field one may simply /// invoke the following code which will launch the calculation in *vacuum*, @@ -74,7 +74,7 @@ /// field.GammaTransmissionProbability( 0.01 ); /// \endcode /// -/// ### In a buffer gas medium +/// ### 2. In a buffer gas medium /// /// The axion-photon probability can also be calculated in the sin of a /// gaseous medium. For that we need to assign a buffer gas instance containing @@ -95,7 +95,14 @@ /// field.GammaTransmissionProbability( 0.1 ); /// \endcode /// -/// \section non-constant Calculating axion-photon probability in an unhomogeneous field. +/// Once we have assigned a buffer gas to the class we may return back to the +/// vacuum state by just assigning a nullptr to the gas. +/// +/// \code +/// field.AssignBufferGas(nullptr); +/// \endcode +/// +/// ## Calculating axion-photon probability in an unhomogeneous field. /// /// There are two main strategies presently implemented inside this class to integrate /// the axion field along an unhomogeneous magnetic field. The first one uses a simple @@ -168,6 +175,20 @@ /// Where the later call will return a `std::pair` with the probability and its /// corresponding error. /// +/// ## Determining density steps for continuous scanning +/// +/// A method may be used to determine the masses and gas densities to achieve +/// a continuous scanning TRestAxionField::GetMassDensityScanning. This method +/// will place a new mass or gas density setting at FWHM/2 till it reaches +/// the maximum axion mass specified in the argument. +/// +/// The following code recovers the axion masses and density settings required +/// for a continuous scan till 0.2eV. +/// +/// \code +/// std::pair aField.GetMassDensityScanning( "He", 0.2 ); +/// \endcode +/// ///-------------------------------------------------------------------------- /// /// RESTsoft - Software for Rare Event Searches with TPCs @@ -176,9 +197,12 @@ /// /// 2019-March: First concept and implementation of TRestAxionField class. /// Javier Galan +/// 2023-December: Implementing methods to recover the axion mass scanning nodes +/// Fran Candon /// /// \class TRestAxionField /// \author Javier Galan +/// \author Fran Candon /// ///
/// @@ -209,9 +233,6 @@ TRestAxionField::~TRestAxionField() {} /////////////////////////////////////////////// /// \brief Initialization of TRestAxionField 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() { fBufferGas = NULL;