From 47258743e126cde667fb200fbe04d4804a65e93f Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 23 Jan 2024 16:04:40 +0100 Subject: [PATCH 01/32] TRestAxionEvent::PrintEvent. Adding mass output --- src/TRestAxionEvent.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/TRestAxionEvent.cxx b/src/TRestAxionEvent.cxx index c5e2c2c0..4e38922e 100644 --- a/src/TRestAxionEvent.cxx +++ b/src/TRestAxionEvent.cxx @@ -159,8 +159,9 @@ void TRestAxionEvent::Translate(const TVector3& delta) { fPosition += delta; } void TRestAxionEvent::PrintEvent() { TRestEvent::PrintEvent(); - cout << "Energy : " << GetEnergy() << endl; - cout << "Position : ( " << fPosition.X() << ", " << fPosition.Y() << ", " << fPosition.Z() << " )" + cout << "Mass : " << GetMass() * units("eV") << " eV" << endl; + cout << "Energy : " << GetEnergy() << " keV" << endl; + cout << "Position : ( " << fPosition.X() << ", " << fPosition.Y() << ", " << fPosition.Z() << " ) mm" << endl; cout << "Direction : ( " << fDirection.X() << ", " << fDirection.Y() << ", " << fDirection.Z() << " )" << endl; From 62a44b3464393658a99507742ef0bc35ec65a58a Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 23 Jan 2024 16:08:55 +0100 Subject: [PATCH 02/32] TRestAxionFieldPropagationProcess. Hot fix --- src/TRestAxionFieldPropagationProcess.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TRestAxionFieldPropagationProcess.cxx b/src/TRestAxionFieldPropagationProcess.cxx index e7de98b3..eaa9ae7a 100644 --- a/src/TRestAxionFieldPropagationProcess.cxx +++ b/src/TRestAxionFieldPropagationProcess.cxx @@ -208,7 +208,7 @@ TRestEvent* TRestAxionFieldPropagationProcess::ProcessEvent(TRestEvent* evInput) fieldAverage = std::accumulate(bProfile.begin(), bProfile.end(), 0.0) / bProfile.size(); Double_t Ea = fAxionEvent->GetEnergy(); - Double_t ma = fAxionEvent->GetMass(); + Double_t ma = fAxionEvent->GetMass() * units("eV"); prob = fAxionField->GammaTransmissionProbability(bProfile, fIntegrationStep, Ea, ma); From 6d8149774cec83f3eb2c8b8645e044e7a342d017 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Wed, 24 Jan 2024 12:05:15 +0100 Subject: [PATCH 03/32] Reviewing scripts to produce magnetic field maps --- scripts/{xlsToBin.py => xlsBikovskiyToBin.py} | 0 scripts/xlsMentinkToBin.py | 152 ++++++++++++++++++ 2 files changed, 152 insertions(+) rename scripts/{xlsToBin.py => xlsBikovskiyToBin.py} (100%) create mode 100755 scripts/xlsMentinkToBin.py diff --git a/scripts/xlsToBin.py b/scripts/xlsBikovskiyToBin.py similarity index 100% rename from scripts/xlsToBin.py rename to scripts/xlsBikovskiyToBin.py diff --git a/scripts/xlsMentinkToBin.py b/scripts/xlsMentinkToBin.py new file mode 100755 index 00000000..8e6a8d1e --- /dev/null +++ b/scripts/xlsMentinkToBin.py @@ -0,0 +1,152 @@ +#!/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 +zCenter = 0.7975 + +print("Starting to read") +# loading the data into a matrix (xyzBdata) +file = r"../../../data/magneticField/Mentink_202401.xlsx" +df = pd.read_excel(file) + +print(df[1:5]) + +print("Translating to matrix") +xyzBdata = df.as_matrix(columns=df.columns[0:]) + +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[1] - yCenter), + 1000 * (x[2] - zCenter), + 1000 * (x[0] - xCenter), + x[4], + x[5], + x[3], + ] + + # 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 z-axis, when we change the sign of z-axis we must change the sign of Bz. + if y[2] > 0: + y[2] = -y[2] + y[5] = -y[5] + 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 9edb0b3367d7fdc55b2265487304917478087af8 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Mon, 5 Feb 2024 16:58:46 +0100 Subject: [PATCH 04/32] TRestAxionMagneticField::ReMap method added --- inc/TRestAxionMagneticField.h | 2 ++ src/TRestAxionMagneticField.cxx | 49 +++++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/inc/TRestAxionMagneticField.h b/inc/TRestAxionMagneticField.h index dfaf80e3..05f703ec 100644 --- a/inc/TRestAxionMagneticField.h +++ b/inc/TRestAxionMagneticField.h @@ -102,6 +102,8 @@ class TRestAxionMagneticField : public TRestMetadata { public: void LoadMagneticVolumes(); + void ReMap( Double_t sX, Double_t sY, Double_t sZ ); + /// 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; diff --git a/src/TRestAxionMagneticField.cxx b/src/TRestAxionMagneticField.cxx index 08473768..d21c405e 100644 --- a/src/TRestAxionMagneticField.cxx +++ b/src/TRestAxionMagneticField.cxx @@ -1059,6 +1059,55 @@ TVector3 TRestAxionMagneticField::GetMagneticField(TVector3 pos, Bool_t showWarn } } +/////////////////////////////////////////////// +/// \brief It allows to remap the magnetic field to a larger mesh size. The new mesh +/// size granularity must be provided by argument, and each dimension must be a factor of the +/// present mesh size. +/// +void TRestAxionMagneticField::ReMap( Double_t sX, Double_t sY, Double_t sZ ) +{ + if( sX == 0 || sY == 0 || sZ == 0 ) + { + RESTError << "TRestAxionMagneticField::ReMap. The mesh granularity cannot be 0" << RESTendl; + RESTError << "Remapping will not take effect" << RESTendl; + return; + } + + Double_t remainder = std::fmod(sX, fMeshSize[0].X()) + std::fmod(sY, fMeshSize[0].Y()) + std::fmod(sZ, fMeshSize[0].Z()); + if( remainder != 0 ) + { + RESTError << "TRestAxionMagneticField::ReMap. The field cannot be remapped."<< RESTendl; + RESTError << "The new mesh granularity must be a multiple of the existing granularity." << RESTendl; + RESTError << "Present mesh size : (" << fMeshSize[0].X() << ", " << fMeshSize[0].Y() << ", " << fMeshSize[0].Z() << ")" << RESTendl; + RESTError << "Requested mesh size : (" << sX << ", " << sY << ", " << sZ << ")" << RESTendl; + RESTError << "Remapping will not take effect" << RESTendl; + return; + } + + Int_t scaleX = (Int_t) (sX/fMeshSize[0].X()); + Int_t scaleY = (Int_t) (sY/fMeshSize[0].Y()); + Int_t scaleZ = (Int_t) (sZ/fMeshSize[0].Z()); + + Int_t newNodesX = (fMagneticFieldVolumes[0].mesh.GetNodesX()-1)/scaleX + 1; + Int_t newNodesY = (fMagneticFieldVolumes[0].mesh.GetNodesY()-1)/scaleY + 1; + Int_t newNodesZ = (fMagneticFieldVolumes[0].mesh.GetNodesZ()-1)/scaleZ + 1; + + for( Int_t nx = 0; nx < newNodesX; nx++ ) + for( Int_t ny = 0; ny < newNodesY; ny++ ) + for( Int_t nz = 0; nz < newNodesZ; nz++ ) + fMagneticFieldVolumes[0].field[nx][ny][nz] = fMagneticFieldVolumes[0].field[nx*scaleX][ny*scaleY][nz*scaleZ]; + + fMagneticFieldVolumes[0].mesh.SetNodes(newNodesX, newNodesY, newNodesZ); + fMagneticFieldVolumes[0].field.resize(fMagneticFieldVolumes[0].mesh.GetNodesX()); + for (unsigned int n = 0; n < fMagneticFieldVolumes[0].field.size(); n++) { + fMagneticFieldVolumes[0].field[n].resize(fMagneticFieldVolumes[0].mesh.GetNodesY()); + for (unsigned int m = 0; m < fMagneticFieldVolumes[0].field[n].size(); m++) + fMagneticFieldVolumes[0].field[n][m].resize(fMagneticFieldVolumes[0].mesh.GetNodesZ()); + } + + fMeshSize[0] = TVector3( fMeshSize[0].X() * scaleX, fMeshSize[0].Y() * scaleY, fMeshSize[0].Z() * scaleZ ); +} + /////////////////////////////////////////////// /// \brief It returns the corresponding volume index at the given position. If not found it will return /// -1. From a98a5bdd65abb04aeaf2bc8e98c75532757fb947 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Fri, 9 Feb 2024 12:49:36 +0100 Subject: [PATCH 05/32] TRestAxionField. Adding fMagneticField pointer --- inc/TRestAxionField.h | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/inc/TRestAxionField.h b/inc/TRestAxionField.h index 83cd950c..21692d1b 100644 --- a/inc/TRestAxionField.h +++ b/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,10 @@ 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; //! public: void SetMagneticField(Double_t b) { fBmag = b; } @@ -62,6 +66,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 +85,18 @@ class TRestAxionField : public TObject { Double_t GammaTransmissionProbability(std::vector Bmag, Double_t deltaL, Double_t Ea, Double_t ma, Double_t mg = 0, Double_t absLength = 0); + Double_t GammaTransmissionFieldMapProbability( Double_t Ea, Double_t ma ); + + /// Integrand used for axion-photon probability integration + static double Integrand(double x, void *params) { + auto *data = reinterpret_cast*>(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> GetMassDensityScanning(std::string gasName = "He", From ecefc75710b00a391d38940198187ea7584c7213 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sat, 10 Feb 2024 20:00:44 +0100 Subject: [PATCH 06/32] TRestAxionField. Implemented integration routines using GSL libraries --- inc/TRestAxionField.h | 6 +- src/TRestAxionField.cxx | 316 ++++++++++++++++++++++++++++++++-------- 2 files changed, 258 insertions(+), 64 deletions(-) diff --git a/inc/TRestAxionField.h b/inc/TRestAxionField.h index 21692d1b..48777a8a 100644 --- a/inc/TRestAxionField.h +++ b/inc/TRestAxionField.h @@ -48,6 +48,10 @@ class TRestAxionField : public TObject { /// A pointer to the magnetic field definition TRestAxionMagneticField* fMagneticField = nullptr; //! + std::pair ComputeOffResonanceIntegral(Double_t q, Double_t Gamma, Double_t accuracy, Int_t num_intervals, Int_t qawo_levels ); + + std::pair ComputeResonanceIntegral( Double_t Gamma, Double_t accuracy, Int_t num_intervals); + public: void SetMagneticField(Double_t b) { fBmag = b; } void SetCoherenceLength(Double_t l) { fLcoh = l; } @@ -85,7 +89,7 @@ class TRestAxionField : public TObject { Double_t GammaTransmissionProbability(std::vector Bmag, Double_t deltaL, Double_t Ea, Double_t ma, Double_t mg = 0, Double_t absLength = 0); - Double_t GammaTransmissionFieldMapProbability( Double_t Ea, Double_t ma ); + std::pair 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) { diff --git a/src/TRestAxionField.cxx b/src/TRestAxionField.cxx index fae12e99..18d0160d 100644 --- a/src/TRestAxionField.cxx +++ b/src/TRestAxionField.cxx @@ -54,6 +54,7 @@ #endif #include +#include using namespace std; @@ -94,9 +95,9 @@ void TRestAxionField::Initialize() { /// The result will be given for an axion-photon coupling of 10^{-10} GeV^{-1} /// double TRestAxionField::BL(Double_t Bmag, Double_t Lcoh) { - Double_t lengthInMeters = Lcoh / 1000.; + Double_t lengthInMeters = Lcoh * units("m"); - Double_t tm = REST_Physics::lightSpeed / REST_Physics::naturalElectron * 1.0e-9; // GeV + Double_t tm = REST_Physics::lightSpeed / REST_Physics::naturalElectron * units("GeV")/units("eV"); // eV --> GeV Double_t sol = lengthInMeters * Bmag * tm; sol = sol * 1.0e-10; @@ -111,9 +112,9 @@ double TRestAxionField::BL(Double_t Bmag, Double_t Lcoh) { /// double TRestAxionField::BLHalfSquared(Double_t Bmag, Double_t Lcoh) // (BL/2)**2 { - Double_t lengthInMeters = Lcoh / 1000.; + Double_t lengthInMeters = Lcoh * units("m"); - Double_t tm = REST_Physics::lightSpeed / REST_Physics::naturalElectron * 1.0e-9; // gev + Double_t tm = REST_Physics::lightSpeed / REST_Physics::naturalElectron * units("GeV") / units( "eV"); // eV --> GeV Double_t sol = lengthInMeters * Bmag * tm / 2; sol = sol * sol * 1.0e-20; @@ -143,7 +144,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, return 0; #else mpfr::mpreal axionMass = ma; - mpfr::mpreal cohLength = fLcoh / 1000.; // Default REST units are mm; + mpfr::mpreal cohLength = fLcoh * units("m"); // Default REST units are mm; mpfr::mpreal photonMass = mg; @@ -160,39 +161,40 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, if (ma == 0.0 && photonMass == 0.0) return BLHalfSquared(fBmag, fLcoh); - mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / fEa / 1000.0; + mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / (fEa * units("eV")); mpfr::mpreal l = cohLength * REST_Physics::PhMeterIneV; mpfr::mpreal phi = q * l; mpfr::mpreal Gamma = absLength; if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(fEa); // cm-1 - mpfr::mpreal GammaL = Gamma * cohLength * 100; + mpfr::mpreal GammaL = Gamma * cohLength * units("cm")/units("m"); // m --> cm if (fDebug) { - RESTDebug << "+------------------------+" << RESTendl; - RESTDebug << " Intermediate calculations" << RESTendl; - RESTDebug << " q : " << q << " eV" << RESTendl; - RESTDebug << " l : " << l << " eV-1" << RESTendl; - RESTDebug << " phi : " << phi << RESTendl; - RESTDebug << "Gamma : " << Gamma << RESTendl; - RESTDebug << "GammaL : " << GammaL << RESTendl; - RESTDebug << "+------------------------+" << RESTendl; + std::cout << "+------------------------+" << std::endl; + std::cout << " Intermediate calculations" << std::endl; + std::cout << " q : " << q << " eV" << std::endl; + std::cout << " l : " << l << " eV-1" << std::endl; + std::cout << " phi : " << phi << std::endl; + std::cout << "Gamma : " << Gamma << std::endl; + std::cout << "GammaL : " << GammaL << std::endl; + std::cout << "+------------------------+" << std::endl; } mpfr::mpreal MFactor = phi * phi + GammaL * GammaL / 4.0; MFactor = 1.0 / MFactor; if (fDebug) { - RESTDebug << "Mfactor : " << MFactor << RESTendl; - RESTDebug << "(BL/2)^2 : " << BLHalfSquared(fBmag, fLcoh) << RESTendl; - RESTDebug << "cos(phi) : " << cos(phi) << RESTendl; - RESTDebug << "Exp(-GammaL) : " << exp(-GammaL) << RESTendl; + std::cout << "Mfactor : " << MFactor << std::endl; + std::cout << "(BL/2)^2 : " << BLHalfSquared(fBmag, fLcoh) << std::endl; + std::cout << "cos(phi) : " << cos(phi) << std::endl; + std::cout << "Exp(-GammaL) : " << exp(-GammaL) << std::endl; } double sol = (double)(MFactor * BLHalfSquared(fBmag, fLcoh) * (1 + exp(-GammaL) - 2 * exp(-GammaL / 2) * cos(phi))); - RESTDebug << "Axion-photon transmission probability : " << sol << RESTendl; + if( fDebug ) + std::cout << "Axion-photon transmission probability : " << sol << std::endl; return sol; #endif @@ -244,7 +246,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma // Default REST units are mm. We express cohLength in m. Double_t Lcoh = (Bmag.size() - 1) * deltaL; // in mm - Double_t cohLength = Lcoh / 1000.; // in m + Double_t cohLength = Lcoh * units("m"); // in m mpfr::mpreal photonMass = mg; @@ -253,69 +255,72 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma Double_t fieldAverage = 0; if (Bmag.size() > 0) fieldAverage = std::accumulate(Bmag.begin(), Bmag.end(), 0.0) / Bmag.size(); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << RESTendl; - RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; - RESTDebug << " Axion mass : " << ma << " eV" << RESTendl; - RESTDebug << " Axion energy : " << Ea << " keV" << RESTendl; - RESTDebug << " Lcoh : " << cohLength << " mm" << RESTendl; - RESTDebug << " Bmag average : " << fieldAverage << " T" << RESTendl; - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; + if( fDebug ) + { + std::cout << "+--------------------------------------------------------------------------+" << std::endl; + std::cout << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << std::endl; + std::cout << " Photon mass : " << photonMass << " eV" << std::endl; + std::cout << " Axion mass : " << ma << " eV" << std::endl; + std::cout << " Axion energy : " << Ea << " keV" << std::endl; + std::cout << " Lcoh : " << cohLength << " m" << std::endl; + std::cout << " Bmag average : " << fieldAverage << " T" << std::endl; + std::cout << "+--------------------------------------------------------------------------+" << std::endl; + } // In vacuum if (ma == 0.0 && photonMass == 0.0) return BLHalfSquared(fieldAverage, Lcoh); - mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / Ea / 1000.0; + mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / (Ea * units("eV")); mpfr::mpreal l = cohLength * REST_Physics::PhMeterIneV; mpfr::mpreal phi = q * l; mpfr::mpreal Gamma = absLength; if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 - mpfr::mpreal GammaL = Gamma * cohLength * 100; + mpfr::mpreal GammaL = Gamma * cohLength * units("cm")/units("m"); // m --> cm if (fDebug) { - RESTDebug << "+------------------------+" << RESTendl; - RESTDebug << " Intermediate calculations" << RESTendl; - RESTDebug << " q : " << q << " eV" << RESTendl; - RESTDebug << " l : " << l << " eV-1" << RESTendl; - RESTDebug << " phi : " << phi << RESTendl; - RESTDebug << "Gamma : " << Gamma << RESTendl; - RESTDebug << "GammaL : " << GammaL << RESTendl; - RESTDebug << "+------------------------+" << RESTendl; + std::cout << "+------------------------+" << std::endl; + std::cout << " Intermediate calculations" << std::endl; + std::cout << " q : " << q << " eV" << std::endl; + std::cout << " l : " << l << " eV-1" << std::endl; + std::cout << " phi : " << phi << std::endl; + std::cout << "Gamma : " << Gamma << std::endl; + std::cout << "GammaL : " << GammaL << std::endl; + std::cout << "+------------------------+" << std::endl; } mpfr::mpreal MFactor = phi * phi + GammaL * GammaL / 4.0; MFactor = 1.0 / MFactor; if (fDebug) { - RESTDebug << "Mfactor : " << MFactor << RESTendl; - RESTDebug << "(BL/2)^2 : " << BLHalfSquared(fieldAverage, Lcoh) << RESTendl; - RESTDebug << "cos(phi) : " << cos(phi) << RESTendl; - RESTDebug << "Exp(-GammaL) : " << exp(-GammaL) << RESTendl; + std::cout << "Mfactor : " << MFactor << std::endl; + std::cout << "(BL/2)^2 : " << BLHalfSquared(fieldAverage, Lcoh) << std::endl; + std::cout << "cos(phi) : " << cos(phi) << std::endl; + std::cout << "Exp(-GammaL) : " << exp(-GammaL) << std::endl; } - Double_t deltaIneV = deltaL / 1000. * REST_Physics::PhMeterIneV; /// We integrate following the Midpoint rule method. (Other potential options : Trapezoidal, Simpsons) + /* + Double_t deltaIneV = deltaL * units("m") * REST_Physics::PhMeterIneV; TRestComplex sum(0, 0); - for (unsigned int n = 0; n < Bmag.size() - 1; n++) { - Double_t Bmiddle = 0.5 * (Bmag[n] + Bmag[n + 1]); + for (unsigned int n = 0; n < Bmag.size() - 1; n++) { + Double_t Bmiddle = 0.5 * (Bmag[n] + Bmag[n + 1]); - Double_t lStepIneV = ((double)n + 0.5) * deltaIneV; - Double_t lStepInCm = ((double)n + 0.5) * deltaL / 10.; + Double_t lStepIneV = ((double)n + 0.5) * deltaIneV; + Double_t lStepInCm = ((double)n + 0.5) * deltaL * units("cm"); - TRestComplex qC(0, -q * lStepIneV); - qC = TRestComplex::Exp(qC); + TRestComplex qCgC(0.5 * Gamma * lStepInCm, -q * lStepIneV); + qCgC = TRestComplex::Exp(qCgC); - TRestComplex gC(0.5 * Gamma * lStepInCm, 0); - gC = TRestComplex::Exp(gC); + TRestComplex integrand = Bmiddle * deltaL * qCgC; // The integrand is in T by mm - TRestComplex integrand = Bmiddle * deltaL * gC * qC; // The integrand is in T by mm + sum += integrand; + } + */ - sum += integrand; - } - - mpfr::mpreal sol = exp(-GammaL) * sum.Rho2() * BLHalfSquared(1, 1); + mpfr::mpreal sol = 0; + // mpfr::mpreal sol = exp(-GammaL) * sum.Rho2() * BLHalfSquared(1, 1); // Now T and mm have been recalculated in natural units using BLHalfSquared(1,1). /* @@ -323,12 +328,197 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma (double)(MFactor * BLHalfSquared(Bmag, Lcoh) * (1 + exp(-GammaL) - 2 * exp(-GammaL / 2) * cos(phi))); */ - RESTDebug << "Axion-photon transmission probability : " << sol << RESTendl; + if( fDebug ) + std::cout << "Axion-photon transmission probability : " << sol << std::endl; return (Double_t)sol; #endif } +/////////////////////////////////////////////// +/// \brief Performs the calculation of axion-photon conversion probability using directly +/// equation (28) from J. Redondo and A. Ringwald, Light shinning through walls. +/// https://arxiv.org/pdf/1011.3741.pdf +/// +/// m_gamma will be obtainned from the buffer gas definition. If no buffer gas has been +/// assigned then the medium will be assumed to be vacuum. +/// +/// Ea in keV, ma in eV +/// +/// The mgamma and absorption lengths are the ones defined by the gas. +/// +/// The returned value is given for g_ag = 10^-10 GeV-1 +/// +/// \note The density is for the moment homogeneous. We would need to implemnent a double integral +/// to solve the problem with a density profile. +/// +std::pair TRestAxionField::GammaTransmissionFieldMapProbability( Double_t Ea, Double_t ma, Double_t accuracy, Int_t num_intervals, Int_t qawo_levels ) { + + if( !fMagneticField ) + { + RESTError << "TRestAxionField::GammaTransmissionFieldMapProbability requires a magnetic field map!" << RESTendl; + RESTError << "Use TRestAxionField::AssignMagneticField method to assign one" << RESTendl; + return {0.0,0.0}; + } + + double photonMass = 0; // Vacuum + if (fBufferGas) photonMass = fBufferGas->GetPhotonMass(Ea); + + if( fDebug ) + { + std::cout << "+--------------------------------------------------------------------------+" << std::endl; + std::cout << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << std::endl; + std::cout << " Photon mass : " << photonMass << " eV" << std::endl; + std::cout << " Axion mass : " << ma << " eV" << std::endl; + std::cout << " Axion energy : " << Ea << " keV" << std::endl; + std::cout << "+--------------------------------------------------------------------------+" << std::endl; + } + + double q = (ma * ma - photonMass * photonMass) / 2. / (Ea * units("eV") ); + q = q / REST_Physics::eVInPhMeter*units("m")/units("mm"); // mm-1 + + double Gamma = 0; + if (fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea) * units("cm") / units("mm"); // mm-1 + + if (fDebug) { + std::cout << "+------------------------+" << std::endl; + std::cout << " Intermediate calculations" << std::endl; + std::cout << " q : " << q << " eV" << std::endl; + std::cout << "Gamma : " << Gamma << std::endl; + std::cout << "+------------------------+" << std::endl; + } + + if( q == 0 ) return ComputeResonanceIntegral( Gamma, accuracy, num_intervals ); + else return ComputeOffResonanceIntegral( q, Gamma, accuracy, num_intervals, qawo_levels); + + return {0.0,0.0}; +} + +/////////////////////////////////////////////// +/// \brief Performs the calculation of axion-photon conversion probability using directly +/// equation (28) from J. Redondo and A. Ringwald, Light shinning through walls. +/// https://arxiv.org/pdf/1011.3741.pdf +/// +/// It integrates the Integrand function defined on the header of TRestAxionField. +/// +/// This method uses the GSL QAG integration method. We use this method when the cosine +/// function vanishes when q = 0. +/// +/// See https://www.gnu.org/software/gsl/doc/html/integration.html for more details on +/// the integration parameters. +/// +/// The Gamma function should be expressed in mm-1 since the field map accessed in the integrand +/// is evaluated using mm. +/// +std::pair TRestAxionField::ComputeResonanceIntegral( Double_t Gamma, Double_t accuracy, Int_t num_intervals ) +{ + double reprob,rerr; + + std::pair params = {fMagneticField, Gamma}; + + gsl_integration_workspace *workspace = gsl_integration_workspace_alloc(num_intervals); + + gsl_function F; + F.function = &Integrand; + F.params = ¶ms; + + auto start = std::chrono::system_clock::now(); + + gsl_integration_qag( &F, 0, fMagneticField->GetTrackLength(), accuracy, accuracy, num_intervals, GSL_INTEG_GAUSS61, workspace, &reprob, &rerr); + + auto end = std::chrono::system_clock::now(); + auto seconds = std::chrono::duration_cast(end-start); + + double GammaL = Gamma * fMagneticField->GetTrackLength(); + double C = exp(-GammaL) * BLHalfSquared(1, 1); + + double prob = C * reprob*reprob; + double proberr = 2 * C * reprob * rerr; + + if( fDebug ) + { + std::cout << " ---- TRestAxionField::ComputeResonanceIntegral (QAG) ----" << std::endl; + std::cout << "Gamma: " << Gamma << " mm-1" << std::endl; + std::cout << "accuracy: " << accuracy << std::endl; + std::cout << "num_intervals: " << num_intervals << std::endl; + std::cout << " -------" << std::endl; + std::cout << "Probability: " << prob << std::endl; + std::cout << "Probability error: " << proberr << std::endl; + std::cout << "Computing time: " << seconds.count() << std::endl; + std::cout << " -------" << std::endl; + } + + std::pair sol = { prob, proberr }; + return sol; +} + +/////////////////////////////////////////////// +/// \brief Performs the calculation of axion-photon conversion probability using directly +/// equation (28) from J. Redondo and A. Ringwald, Light shinning through walls. +/// https://arxiv.org/pdf/1011.3741.pdf +/// +/// It integrates the Integrand function defined on the header of TRestAxionField. +/// +/// This method uses the GSL QAWO method for oscillatory functions. That is the case when +/// the q is not zero, in the off-resonance case. +/// +/// See https://www.gnu.org/software/gsl/doc/html/integration.html for more details on +/// the integration parameters. +/// +/// The Gamma function should be expressed in mm-1 since the field map accessed in the integrand +/// is evaluated using mm. +/// +std::pair TRestAxionField::ComputeOffResonanceIntegral(Double_t q, Double_t Gamma, Double_t accuracy, Int_t num_intervals, Int_t qawo_levels ) +{ + double reprob,rerr; + double improb,imerr; + + std::pair params = {fMagneticField, Gamma}; + + gsl_integration_workspace *workspace = gsl_integration_workspace_alloc(num_intervals); + + gsl_function F; + F.function = &Integrand; + F.params = ¶ms; + + auto start = std::chrono::system_clock::now(); + + gsl_integration_qawo_table *wf = gsl_integration_qawo_table_alloc(q, fMagneticField->GetTrackLength(), GSL_INTEG_COSINE, qawo_levels); + gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &reprob, &rerr); + + gsl_integration_qawo_table_set( wf, q, fMagneticField->GetTrackLength(), GSL_INTEG_SINE ); + gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &improb, &imerr); + + gsl_integration_qawo_table_free(wf); + + auto end = std::chrono::system_clock::now(); + auto seconds = std::chrono::duration_cast(end-start); + + double GammaL = Gamma * fMagneticField->GetTrackLength(); + double C = exp(-GammaL) * BLHalfSquared(1, 1); + + double prob = C * (reprob*reprob + improb*improb); + double proberr = 2 * C * TMath::Sqrt( reprob*reprob * rerr*rerr + improb*improb * imerr*imerr ); + + if( fDebug ) + { + std::cout << " ---- TRestAxionField::ComputeOffResonanceIntegral (QAWO) ----" << std::endl; + std::cout << "Gamma: " << Gamma << " mm-1" << std::endl; + std::cout << "q: " << q << "mm-1" << std::endl; + std::cout << "accuracy: " << accuracy << std::endl; + std::cout << "num_intervals: " << num_intervals << std::endl; + std::cout << "qawo_levels: " << qawo_levels << std::endl; + std::cout << " -------" << std::endl; + std::cout << "Probability: " << prob << std::endl; + std::cout << "Probability error: " << proberr << std::endl; + std::cout << "Computing time: " << seconds.count() << std::endl; + std::cout << " -------" << std::endl; + } + + std::pair sol = { prob, proberr }; + return sol; +} + /////////////////////////////////////////////// /// \brief Performs the calculation of axion-photon absorption probability using directly /// equation (18) from van Bibber, Phys Rev D Part Fields. 1989. @@ -347,11 +537,11 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, D RESTWarning << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFr=ON to your REST compilation" << RESTendl; - RESTWarning << "TRestAxionField::GammaTransmissionProbability will return 0" << RESTendl; + RESTWarning << "TRestAxionField::AxionAbsorptionProbability will return 0" << RESTendl; return 0; #else mpfr::mpreal axionMass = ma; - mpfr::mpreal cohLength = fLcoh / 1000.; // Default REST units are mm; + mpfr::mpreal cohLength = fLcoh * units("m"); // Default REST units are mm; mpfr::mpreal photonMass = mg; if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(fEa); @@ -359,7 +549,7 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, D if (fDebug) { RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << RESTendl; + RESTDebug << " TRestAxionField::AxionAbsorptionProbability. Parameter summary" << RESTendl; RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; RESTDebug << " Axion mass : " << ma << " eV" << RESTendl; RESTDebug << " Axion energy : " << fEa << " keV" << RESTendl; @@ -371,13 +561,13 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, D if (ma == 0.0 && photonMass == 0.0) return BLHalfSquared(fBmag, fLcoh); - mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / fEa / 1000.0; + mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / (fEa * units("eV")); mpfr::mpreal l = cohLength * REST_Physics::PhMeterIneV; mpfr::mpreal phi = q * l; mpfr::mpreal Gamma = absLength; if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(fEa); // cm-1 - mpfr::mpreal GammaL = Gamma * cohLength * 100; + mpfr::mpreal GammaL = Gamma * cohLength * units("cm")/units("m"); if (fDebug) { RESTDebug << "+------------------------+" << RESTendl; From 62d187a80db28b787af88c31b6ea7dcad8cc558e Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sat, 10 Feb 2024 20:02:04 +0100 Subject: [PATCH 07/32] TRestAxionMagneticField::SetTrack/GetTransversalComponentInParametricTrack(Double_t x) implemented to parameterize the field component along a track --- inc/TRestAxionMagneticField.h | 16 +++++++++++++++ src/TRestAxionMagneticField.cxx | 35 +++++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/inc/TRestAxionMagneticField.h b/inc/TRestAxionMagneticField.h index 05f703ec..c29f638f 100644 --- a/inc/TRestAxionMagneticField.h +++ b/inc/TRestAxionMagneticField.h @@ -70,6 +70,15 @@ class TRestAxionMagneticField : public TRestMetadata { /// A magnetic field volume structure to store field data and mesh. std::vector 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; //! @@ -104,6 +113,12 @@ class TRestAxionMagneticField : public TRestMetadata { void ReMap( Double_t sX, Double_t sY, Double_t sZ ); + 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; @@ -137,6 +152,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 GetTransversalComponentAlongPath(TVector3 from, TVector3 to, Double_t dl = 1., Int_t Nmax = 0); diff --git a/src/TRestAxionMagneticField.cxx b/src/TRestAxionMagneticField.cxx index d21c405e..51af0918 100644 --- a/src/TRestAxionMagneticField.cxx +++ b/src/TRestAxionMagneticField.cxx @@ -1192,6 +1192,41 @@ std::vector TRestAxionMagneticField::GetTransversalComponentAlongPath( return Bt; } +/////////////////////////////////////////////// +/// \brief It initializes the field track boundaries data members of this class using a +/// track position and direction so that these values can be used later on in parameterization. +/// +void TRestAxionMagneticField::SetTrack( const TVector3 &position, const TVector3 &direction ) +{ + std::vector trackBounds = GetFieldBoundaries(position, direction); + + if( trackBounds.size() != 2 ) + { + fTrackStart = TVector3(0,0,0); + fTrackDirection = TVector3(0,0,0); + fTrackLength = 0; + } + + fTrackStart = trackBounds[0]; + fTrackLength = (trackBounds[1] - trackBounds[0]).Mag() - 1; + fTrackDirection = (trackBounds[1] - trackBounds[0]).Unit(); +} + +/////////////////////////////////////////////// +/// \brief It will return the transversal magnetic field component evaluated at a parametric +/// distance `x` (given by argument) for the track defined inside the class. The track will +/// be defined by the data members fStartTrack and fEndTrack which should be initialized by +/// external means by invoking the method SetTrack( position, direction ); +/// +/// This method will be used by the integration method +/// +Double_t TRestAxionMagneticField::GetTransversalComponentInParametricTrack( Double_t x ) { + if( x < 0 || x > fTrackLength ) + return 0; + + return GetTransversalComponent(fTrackStart + x * fTrackDirection, fTrackDirection); +} + /////////////////////////////////////////////// /// \brief It returns the average of the transversal magnetic field intensity between the 3-d coordinates /// `from` and `to`. From 2bac74a52da4f11627b9eb0dd2f655b12222261c Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sat, 10 Feb 2024 20:09:46 +0100 Subject: [PATCH 08/32] TRestAxionMagneticField::ReMap. Adapting to any number of magnetic volumes --- src/TRestAxionMagneticField.cxx | 65 +++++++++++++++++---------------- 1 file changed, 34 insertions(+), 31 deletions(-) diff --git a/src/TRestAxionMagneticField.cxx b/src/TRestAxionMagneticField.cxx index 51af0918..644f6106 100644 --- a/src/TRestAxionMagneticField.cxx +++ b/src/TRestAxionMagneticField.cxx @@ -1073,39 +1073,42 @@ void TRestAxionMagneticField::ReMap( Double_t sX, Double_t sY, Double_t sZ ) return; } - Double_t remainder = std::fmod(sX, fMeshSize[0].X()) + std::fmod(sY, fMeshSize[0].Y()) + std::fmod(sZ, fMeshSize[0].Z()); - if( remainder != 0 ) + for( size_t n = 0; n < fMagneticFieldVolumes.size(); n++ ) { - RESTError << "TRestAxionMagneticField::ReMap. The field cannot be remapped."<< RESTendl; - RESTError << "The new mesh granularity must be a multiple of the existing granularity." << RESTendl; - RESTError << "Present mesh size : (" << fMeshSize[0].X() << ", " << fMeshSize[0].Y() << ", " << fMeshSize[0].Z() << ")" << RESTendl; - RESTError << "Requested mesh size : (" << sX << ", " << sY << ", " << sZ << ")" << RESTendl; - RESTError << "Remapping will not take effect" << RESTendl; - return; + Double_t remainder = std::fmod(sX, fMeshSize[n].X()) + std::fmod(sY, fMeshSize[n].Y()) + std::fmod(sZ, fMeshSize[n].Z()); + if( remainder != 0 ) + { + RESTError << "TRestAxionMagneticField::ReMap. The field cannot be remapped."<< RESTendl; + RESTError << "The new mesh granularity must be a multiple of the existing granularity." << RESTendl; + RESTError << "Present mesh size : (" << fMeshSize[n].X() << ", " << fMeshSize[n].Y() << ", " << fMeshSize[n].Z() << ")" << RESTendl; + RESTError << "Requested mesh size : (" << sX << ", " << sY << ", " << sZ << ")" << RESTendl; + RESTError << "Remapping will not take effect" << RESTendl; + return; + } + + Int_t scaleX = (Int_t) (sX/fMeshSize[n].X()); + Int_t scaleY = (Int_t) (sY/fMeshSize[n].Y()); + Int_t scaleZ = (Int_t) (sZ/fMeshSize[n].Z()); + + Int_t newNodesX = (fMagneticFieldVolumes[n].mesh.GetNodesX()-1)/scaleX + 1; + Int_t newNodesY = (fMagneticFieldVolumes[n].mesh.GetNodesY()-1)/scaleY + 1; + Int_t newNodesZ = (fMagneticFieldVolumes[n].mesh.GetNodesZ()-1)/scaleZ + 1; + + for( Int_t nx = 0; nx < newNodesX; nx++ ) + for( Int_t ny = 0; ny < newNodesY; ny++ ) + for( Int_t nz = 0; nz < newNodesZ; nz++ ) + fMagneticFieldVolumes[n].field[nx][ny][nz] = fMagneticFieldVolumes[n].field[nx*scaleX][ny*scaleY][nz*scaleZ]; + + fMagneticFieldVolumes[n].mesh.SetNodes(newNodesX, newNodesY, newNodesZ); + fMagneticFieldVolumes[n].field.resize(fMagneticFieldVolumes[n].mesh.GetNodesX()); + for (unsigned int i = 0; i < fMagneticFieldVolumes[n].field.size(); n++) { + fMagneticFieldVolumes[n].field[n].resize(fMagneticFieldVolumes[n].mesh.GetNodesY()); + for (unsigned int j = 0; j < fMagneticFieldVolumes[n].field[i].size(); m++) + fMagneticFieldVolumes[n].field[i][j].resize(fMagneticFieldVolumes[n].mesh.GetNodesZ()); + } + + fMeshSize[n] = TVector3( fMeshSize[n].X() * scaleX, fMeshSize[n].Y() * scaleY, fMeshSize[n].Z() * scaleZ ); } - - Int_t scaleX = (Int_t) (sX/fMeshSize[0].X()); - Int_t scaleY = (Int_t) (sY/fMeshSize[0].Y()); - Int_t scaleZ = (Int_t) (sZ/fMeshSize[0].Z()); - - Int_t newNodesX = (fMagneticFieldVolumes[0].mesh.GetNodesX()-1)/scaleX + 1; - Int_t newNodesY = (fMagneticFieldVolumes[0].mesh.GetNodesY()-1)/scaleY + 1; - Int_t newNodesZ = (fMagneticFieldVolumes[0].mesh.GetNodesZ()-1)/scaleZ + 1; - - for( Int_t nx = 0; nx < newNodesX; nx++ ) - for( Int_t ny = 0; ny < newNodesY; ny++ ) - for( Int_t nz = 0; nz < newNodesZ; nz++ ) - fMagneticFieldVolumes[0].field[nx][ny][nz] = fMagneticFieldVolumes[0].field[nx*scaleX][ny*scaleY][nz*scaleZ]; - - fMagneticFieldVolumes[0].mesh.SetNodes(newNodesX, newNodesY, newNodesZ); - fMagneticFieldVolumes[0].field.resize(fMagneticFieldVolumes[0].mesh.GetNodesX()); - for (unsigned int n = 0; n < fMagneticFieldVolumes[0].field.size(); n++) { - fMagneticFieldVolumes[0].field[n].resize(fMagneticFieldVolumes[0].mesh.GetNodesY()); - for (unsigned int m = 0; m < fMagneticFieldVolumes[0].field[n].size(); m++) - fMagneticFieldVolumes[0].field[n][m].resize(fMagneticFieldVolumes[0].mesh.GetNodesZ()); - } - - fMeshSize[0] = TVector3( fMeshSize[0].X() * scaleX, fMeshSize[0].Y() * scaleY, fMeshSize[0].Z() * scaleZ ); } /////////////////////////////////////////////// From 187ac3c999419316340189d5c7a1e838d9e51daf Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sat, 10 Feb 2024 22:39:21 +0100 Subject: [PATCH 09/32] TRestAxionField. mpfr to double type conversion --- src/TRestAxionField.cxx | 53 +++++++++++++++++++---------------------- 1 file changed, 25 insertions(+), 28 deletions(-) diff --git a/src/TRestAxionField.cxx b/src/TRestAxionField.cxx index 18d0160d..76b8bc38 100644 --- a/src/TRestAxionField.cxx +++ b/src/TRestAxionField.cxx @@ -143,10 +143,9 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, RESTWarning << "TRestAxionField::GammaTransmissionProbability will return 0" << RESTendl; return 0; #else - mpfr::mpreal axionMass = ma; - mpfr::mpreal cohLength = fLcoh * units("m"); // Default REST units are mm; + Double_t cohLength = fLcoh * units("m"); // Default REST units are mm; - mpfr::mpreal photonMass = mg; + Double_t photonMass = mg; if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(fEa); @@ -161,13 +160,13 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, if (ma == 0.0 && photonMass == 0.0) return BLHalfSquared(fBmag, fLcoh); - mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / (fEa * units("eV")); - mpfr::mpreal l = cohLength * REST_Physics::PhMeterIneV; - mpfr::mpreal phi = q * l; + Double_t q = (ma * ma - photonMass * photonMass) / 2. / (fEa * units("eV")); + Double_t l = cohLength * REST_Physics::PhMeterIneV; + Double_t phi = q * l; - mpfr::mpreal Gamma = absLength; + Double_t Gamma = absLength; if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(fEa); // cm-1 - mpfr::mpreal GammaL = Gamma * cohLength * units("cm")/units("m"); // m --> cm + Double_t GammaL = Gamma * cohLength * units("cm")/units("m"); // m --> cm if (fDebug) { std::cout << "+------------------------+" << std::endl; @@ -180,7 +179,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, std::cout << "+------------------------+" << std::endl; } - mpfr::mpreal MFactor = phi * phi + GammaL * GammaL / 4.0; + Double_t MFactor = phi * phi + GammaL * GammaL / 4.0; MFactor = 1.0 / MFactor; if (fDebug) { @@ -242,13 +241,12 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma RESTWarning << "TRestAxionField::GammaTransmissionProbability will return 0" << RESTendl; return 0; #else - mpfr::mpreal axionMass = ma; // Default REST units are mm. We express cohLength in m. Double_t Lcoh = (Bmag.size() - 1) * deltaL; // in mm Double_t cohLength = Lcoh * units("m"); // in m - mpfr::mpreal photonMass = mg; + Double_t photonMass = mg; if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(Ea); @@ -270,13 +268,13 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma // In vacuum if (ma == 0.0 && photonMass == 0.0) return BLHalfSquared(fieldAverage, Lcoh); - mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / (Ea * units("eV")); - mpfr::mpreal l = cohLength * REST_Physics::PhMeterIneV; - mpfr::mpreal phi = q * l; + Double_t q = (ma * ma - photonMass * photonMass) / 2. / (Ea * units("eV")); + Double_t l = cohLength * REST_Physics::PhMeterIneV; + Double_t phi = q * l; - mpfr::mpreal Gamma = absLength; + Double_t Gamma = absLength; if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 - mpfr::mpreal GammaL = Gamma * cohLength * units("cm")/units("m"); // m --> cm + Double_t GammaL = Gamma * cohLength * units("cm")/units("m"); // m --> cm if (fDebug) { std::cout << "+------------------------+" << std::endl; @@ -289,7 +287,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma std::cout << "+------------------------+" << std::endl; } - mpfr::mpreal MFactor = phi * phi + GammaL * GammaL / 4.0; + Double_t MFactor = phi * phi + GammaL * GammaL / 4.0; MFactor = 1.0 / MFactor; if (fDebug) { @@ -319,8 +317,8 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma } */ - mpfr::mpreal sol = 0; - // mpfr::mpreal sol = exp(-GammaL) * sum.Rho2() * BLHalfSquared(1, 1); + Double_t sol = 0; + // Double_t sol = exp(-GammaL) * sum.Rho2() * BLHalfSquared(1, 1); // Now T and mm have been recalculated in natural units using BLHalfSquared(1,1). /* @@ -540,10 +538,9 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, D RESTWarning << "TRestAxionField::AxionAbsorptionProbability will return 0" << RESTendl; return 0; #else - mpfr::mpreal axionMass = ma; - mpfr::mpreal cohLength = fLcoh * units("m"); // Default REST units are mm; + Double_t cohLength = fLcoh * units("m"); // Default REST units are mm; - mpfr::mpreal photonMass = mg; + Double_t photonMass = mg; if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(fEa); if (fDebug) { @@ -561,13 +558,13 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, D if (ma == 0.0 && photonMass == 0.0) return BLHalfSquared(fBmag, fLcoh); - mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / (fEa * units("eV")); - mpfr::mpreal l = cohLength * REST_Physics::PhMeterIneV; - mpfr::mpreal phi = q * l; + Double_t q = (ma * ma - photonMass * photonMass) / 2. / (fEa * units("eV")); + Double_t l = cohLength * REST_Physics::PhMeterIneV; + Double_t phi = q * l; - mpfr::mpreal Gamma = absLength; + Double_t Gamma = absLength; if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(fEa); // cm-1 - mpfr::mpreal GammaL = Gamma * cohLength * units("cm")/units("m"); + Double_t GammaL = Gamma * cohLength * units("cm")/units("m"); if (fDebug) { RESTDebug << "+------------------------+" << RESTendl; @@ -580,7 +577,7 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, D RESTDebug << "+------------------------+" << RESTendl; } - mpfr::mpreal MFactor = phi * phi + GammaL * GammaL / 4.0; + Double_t MFactor = phi * phi + GammaL * GammaL / 4.0; MFactor = 1.0 / MFactor; if (fDebug) { From 6c5032187e1ae45bdea1d9591f217a034069dc2c Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sat, 10 Feb 2024 22:41:18 +0100 Subject: [PATCH 10/32] TRestAxionMagneticField. Adding fReMap to allow remapping from RML config file --- inc/TRestAxionMagneticField.h | 7 ++- src/TRestAxionMagneticField.cxx | 83 ++++++++++++++++++--------------- 2 files changed, 51 insertions(+), 39 deletions(-) diff --git a/inc/TRestAxionMagneticField.h b/inc/TRestAxionMagneticField.h index c29f638f..de43ded3 100644 --- a/inc/TRestAxionMagneticField.h +++ b/inc/TRestAxionMagneticField.h @@ -58,7 +58,7 @@ class TRestAxionMagneticField : public TRestMetadata { /// A constant field component that will be added to the field map std::vector 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 fMeshSize; //< /// The type of the mesh used (default is cylindrical) @@ -67,6 +67,9 @@ class TRestAxionMagneticField : public TRestMetadata { /// A vector to store the maximum bounding box values std::vector 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 fMagneticFieldVolumes; //! @@ -111,7 +114,7 @@ class TRestAxionMagneticField : public TRestMetadata { public: void LoadMagneticVolumes(); - void ReMap( Double_t sX, Double_t sY, Double_t sZ ); + void ReMap( const size_t &n, const TVector3 &newMapSize ); void SetTrack( const TVector3 &position, const TVector3 &direction ); diff --git a/src/TRestAxionMagneticField.cxx b/src/TRestAxionMagneticField.cxx index 644f6106..ae5ce638 100644 --- a/src/TRestAxionMagneticField.cxx +++ b/src/TRestAxionMagneticField.cxx @@ -918,6 +918,11 @@ void TRestAxionMagneticField::LoadMagneticVolumes() { RESTError << "TRestAxionMagneticField::LoadMagneticVolumes. Volumes overlap!" << RESTendl; exit(1); } + + for( size_t n = 0; n < fMagneticFieldVolumes.size(); n++ ) + if( fReMap != TVector3(0,0,0) ) + ReMap( n, fReMap ); + RESTDebug << "Finished loading magnetic volumes" << RESTendl; } @@ -1064,51 +1069,48 @@ TVector3 TRestAxionMagneticField::GetMagneticField(TVector3 pos, Bool_t showWarn /// size granularity must be provided by argument, and each dimension must be a factor of the /// present mesh size. /// -void TRestAxionMagneticField::ReMap( Double_t sX, Double_t sY, Double_t sZ ) +void TRestAxionMagneticField::ReMap( const size_t &n, const TVector3 &newMapSize ) { - if( sX == 0 || sY == 0 || sZ == 0 ) + if( newMapSize.X() == 0 || newMapSize.Y() == 0 || newMapSize.Z() == 0 ) { RESTError << "TRestAxionMagneticField::ReMap. The mesh granularity cannot be 0" << RESTendl; RESTError << "Remapping will not take effect" << RESTendl; return; } - for( size_t n = 0; n < fMagneticFieldVolumes.size(); n++ ) + Double_t remainder = std::fmod(newMapSize.X(), fMeshSize[n].X()) + std::fmod(newMapSize.Y(), fMeshSize[n].Y()) + std::fmod(newMapSize.Z(), fMeshSize[n].Z()); + if( remainder != 0 ) { - Double_t remainder = std::fmod(sX, fMeshSize[n].X()) + std::fmod(sY, fMeshSize[n].Y()) + std::fmod(sZ, fMeshSize[n].Z()); - if( remainder != 0 ) - { - RESTError << "TRestAxionMagneticField::ReMap. The field cannot be remapped."<< RESTendl; - RESTError << "The new mesh granularity must be a multiple of the existing granularity." << RESTendl; - RESTError << "Present mesh size : (" << fMeshSize[n].X() << ", " << fMeshSize[n].Y() << ", " << fMeshSize[n].Z() << ")" << RESTendl; - RESTError << "Requested mesh size : (" << sX << ", " << sY << ", " << sZ << ")" << RESTendl; - RESTError << "Remapping will not take effect" << RESTendl; - return; - } - - Int_t scaleX = (Int_t) (sX/fMeshSize[n].X()); - Int_t scaleY = (Int_t) (sY/fMeshSize[n].Y()); - Int_t scaleZ = (Int_t) (sZ/fMeshSize[n].Z()); - - Int_t newNodesX = (fMagneticFieldVolumes[n].mesh.GetNodesX()-1)/scaleX + 1; - Int_t newNodesY = (fMagneticFieldVolumes[n].mesh.GetNodesY()-1)/scaleY + 1; - Int_t newNodesZ = (fMagneticFieldVolumes[n].mesh.GetNodesZ()-1)/scaleZ + 1; - - for( Int_t nx = 0; nx < newNodesX; nx++ ) - for( Int_t ny = 0; ny < newNodesY; ny++ ) - for( Int_t nz = 0; nz < newNodesZ; nz++ ) - fMagneticFieldVolumes[n].field[nx][ny][nz] = fMagneticFieldVolumes[n].field[nx*scaleX][ny*scaleY][nz*scaleZ]; - - fMagneticFieldVolumes[n].mesh.SetNodes(newNodesX, newNodesY, newNodesZ); - fMagneticFieldVolumes[n].field.resize(fMagneticFieldVolumes[n].mesh.GetNodesX()); - for (unsigned int i = 0; i < fMagneticFieldVolumes[n].field.size(); n++) { - fMagneticFieldVolumes[n].field[n].resize(fMagneticFieldVolumes[n].mesh.GetNodesY()); - for (unsigned int j = 0; j < fMagneticFieldVolumes[n].field[i].size(); m++) - fMagneticFieldVolumes[n].field[i][j].resize(fMagneticFieldVolumes[n].mesh.GetNodesZ()); - } - - fMeshSize[n] = TVector3( fMeshSize[n].X() * scaleX, fMeshSize[n].Y() * scaleY, fMeshSize[n].Z() * scaleZ ); + RESTError << "TRestAxionMagneticField::ReMap. The field cannot be remapped."<< RESTendl; + RESTError << "The new mesh granularity must be a multiple of the existing granularity." << RESTendl; + RESTError << "Present mesh size : (" << fMeshSize[n].X() << ", " << fMeshSize[n].Y() << ", " << fMeshSize[n].Z() << ")" << RESTendl; + RESTError << "Requested mesh size : (" << newMapSize.X() << ", " << newMapSize.Y() << ", " << newMapSize.Z() << ")" << RESTendl; + RESTError << "Remapping will not take effect" << RESTendl; + return; } + + Int_t scaleX = (Int_t) (newMapSize.X()/fMeshSize[n].X()); + Int_t scaleY = (Int_t) (newMapSize.Y()/fMeshSize[n].Y()); + Int_t scaleZ = (Int_t) (newMapSize.Z()/fMeshSize[n].Z()); + + Int_t newNodesX = (fMagneticFieldVolumes[n].mesh.GetNodesX()-1)/scaleX + 1; + Int_t newNodesY = (fMagneticFieldVolumes[n].mesh.GetNodesY()-1)/scaleY + 1; + Int_t newNodesZ = (fMagneticFieldVolumes[n].mesh.GetNodesZ()-1)/scaleZ + 1; + + for( Int_t nx = 0; nx < newNodesX; nx++ ) + for( Int_t ny = 0; ny < newNodesY; ny++ ) + for( Int_t nz = 0; nz < newNodesZ; nz++ ) + fMagneticFieldVolumes[n].field[nx][ny][nz] = fMagneticFieldVolumes[n].field[nx*scaleX][ny*scaleY][nz*scaleZ]; + + fMagneticFieldVolumes[n].mesh.SetNodes(newNodesX, newNodesY, newNodesZ); + fMagneticFieldVolumes[n].field.resize(fMagneticFieldVolumes[n].mesh.GetNodesX()); + for (unsigned int i = 0; i < fMagneticFieldVolumes[n].field.size(); i++) { + fMagneticFieldVolumes[n].field[n].resize(fMagneticFieldVolumes[n].mesh.GetNodesY()); + for (unsigned int j = 0; j < fMagneticFieldVolumes[n].field[i].size(); j++) + fMagneticFieldVolumes[n].field[i][j].resize(fMagneticFieldVolumes[n].mesh.GetNodesZ()); + } + + fMeshSize[n] = TVector3( fMeshSize[n].X() * scaleX, fMeshSize[n].Y() * scaleY, fMeshSize[n].Z() * scaleZ ); } /////////////////////////////////////////////// @@ -1414,7 +1416,8 @@ std::vector TRestAxionMagneticField::GetFieldBoundaries(TVector3 pos, /// \brief Initialization of TRestAxionMagnetic field members through a RML file /// void TRestAxionMagneticField::InitFromConfigFile() { - this->Initialize(); + + TRestMetadata::InitFromConfigFile(); auto magVolumeDef = GetElement("addMagneticVolume"); while (magVolumeDef) { @@ -1487,6 +1490,12 @@ void TRestAxionMagneticField::InitFromConfigFile() { void TRestAxionMagneticField::PrintMetadata() { TRestMetadata::PrintMetadata(); + if( fReMap != TVector3(0,0,0) ) + { + RESTMetadata << "Fields original mesh size has been remapped" << RESTendl; + RESTMetadata << " " << RESTendl; + } + RESTMetadata << " - Number of magnetic volumes : " << GetNumberOfVolumes() << RESTendl; RESTMetadata << " ------------------------------------------------ " << RESTendl; for (unsigned int p = 0; p < GetNumberOfVolumes(); p++) { From 509dd315a92ac74f786e3874d0ffb6cc49b34fb5 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sat, 10 Feb 2024 23:11:44 +0100 Subject: [PATCH 11/32] TRestAxionFieldPropagationProcess. Removing fIntegrationStep --- inc/TRestAxionFieldPropagationProcess.h | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/inc/TRestAxionFieldPropagationProcess.h b/inc/TRestAxionFieldPropagationProcess.h index c12e3ac8..4b41bb38 100644 --- a/inc/TRestAxionFieldPropagationProcess.h +++ b/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; //! From fbf5da1108926da0195a181fcde509d99decb3b8 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sat, 10 Feb 2024 23:12:25 +0100 Subject: [PATCH 12/32] TRestAxionFieldPropagationProcess. Implementing the new GSL integration methods --- src/TRestAxionFieldPropagationProcess.cxx | 72 +++++++++-------------- 1 file changed, 27 insertions(+), 45 deletions(-) diff --git a/src/TRestAxionFieldPropagationProcess.cxx b/src/TRestAxionFieldPropagationProcess.cxx index eaa9ae7a..43e6f62b 100644 --- a/src/TRestAxionFieldPropagationProcess.cxx +++ b/src/TRestAxionFieldPropagationProcess.cxx @@ -56,18 +56,20 @@ /// /// Two metadata parameters can be defined inside this process: /// -/// * **integrationStep** (Default:50mm): The integration length used for field integration along the particle -/// track. /// * **bufferGasAdditionalLength** (Default:0mm): In the case we use a buffer gas we may include an /// additional length that the particle needs to travel inside that buffer gas but outside the magnetic field /// volume, the result will be written as an independent efficiency at the `transmission` observable inside /// the analysis tree. +/// * **accuracy** (Default:0.1): The integration method will try to get an error lower than 10%. +/// * **numIntervals** (Default:100): The max number of intervals used by the integration method. +/// * **qawoLevels** (Default:20): The size of pre-generated table for cosine integrals. +/// * **reMap** (Default: (30,30,100)mm): The mesh cell size of the magnetic field will be re-mapped to +/// the value provided in this parameter. /// /// List of observables generated by this process: /// -/// * **fieldAverage**: The averaged magnetic field calculated along the particle track. /// * **probability**: The final axion-photon conversion probability. -/// * **coherenceLength**: The lentgh of magnetic field region traversed by the particle. +/// * **error**: The error on the final axion-photon conversion probability. /// * **transmission**: This observable will register the photon transmission produced by an additional /// buffer gas length at the end of the magnetic region. /// @@ -163,8 +165,6 @@ void TRestAxionFieldPropagationProcess::InitProcess() { fMagneticField = (TRestAxionMagneticField*)this->GetMetadata("TRestAxionMagneticField"); - RESTDebug << "Magnetic field : " << fMagneticField << RESTendl; - if (!fMagneticField) { RESTError << "TRestAxionFieldPropagationprocess. Magnetic Field was not defined!" << RESTendl; exit(0); @@ -178,10 +178,13 @@ void TRestAxionFieldPropagationProcess::InitProcess() { fAxionField->AssignBufferGas(fBufferGas); else fBufferGasAdditionalLength = 0; + + fAxionField->AssignMagneticField(fMagneticField); } RESTDebug << "Axion-field : " << fAxionField << RESTendl; RESTDebug << "Buffer gas : " << fBufferGas << RESTendl; + RESTDebug << "Magnetic field : " << fMagneticField << RESTendl; } TRestEvent* TRestAxionFieldPropagationProcess::ProcessEvent(TRestEvent* evInput) { @@ -189,51 +192,30 @@ TRestEvent* TRestAxionFieldPropagationProcess::ProcessEvent(TRestEvent* evInput) RESTDebug << "TRestAxionFieldPropagationProcess::ProcessEvent : " << fAxionEvent->GetID() << RESTendl; - std::vector trackBounds = - fMagneticField->GetFieldBoundaries(fAxionEvent->GetPosition(), fAxionEvent->GetDirection()); + fMagneticField->SetTrack( fAxionEvent->GetPosition(), fAxionEvent->GetDirection() ); + + Double_t Ea = fAxionEvent->GetEnergy(); + Double_t ma = fAxionEvent->GetMass() * units("eV"); + + std::pair prob = fAxionField->GammaTransmissionFieldMapProbability( Ea, ma, fAccuracy, fNumIntervals, fQawoLevels ); + + SetObservableValue("probability", prob.first); + SetObservableValue("error", prob.second); - Double_t prob = 0; - Double_t lCoh = 0; Double_t transmission = 1; - Double_t fieldAverage = 0; - if (trackBounds.size() == 2) { - RESTDebug << "-- Track bounds" << RESTendl; - RESTDebug << "X1:" << trackBounds[0].X() << " Y1: " << trackBounds[0].Y() - << " Z1: " << trackBounds[0].Z() << RESTendl; - RESTDebug << "X2:" << trackBounds[1].X() << " Y2: " << trackBounds[1].Y() - << " Z2: " << trackBounds[1].Z() << RESTendl; - std::vector bProfile = fMagneticField->GetTransversalComponentAlongPath( - trackBounds[0], trackBounds[1], fIntegrationStep); - - fieldAverage = std::accumulate(bProfile.begin(), bProfile.end(), 0.0) / bProfile.size(); - - Double_t Ea = fAxionEvent->GetEnergy(); - Double_t ma = fAxionEvent->GetMass() * units("eV"); - - prob = fAxionField->GammaTransmissionProbability(bProfile, fIntegrationStep, Ea, ma); - - lCoh = (bProfile.size() - 1) * fIntegrationStep; - - if (fBufferGas && fBufferGasAdditionalLength > 0) { - Double_t Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 - Double_t GammaL = Gamma * lCoh * units("cm"); - transmission = exp(-GammaL); - } - } else { - SetWarning("TRestAxionFieldPropagationProcess. Track does not cross the field volume!", false); - } + if (fBufferGas && fBufferGasAdditionalLength > 0) { + Double_t Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 + Double_t GammaL = Gamma * fBufferGasAdditionalLength * units("cm"); + transmission = exp(-GammaL); + } + + SetObservableValue("transmission", transmission); RESTDebug << " --- Process observables: " << RESTendl; - RESTDebug << "Field average: " << fieldAverage << " T" << RESTendl; - RESTDebug << "Probability: " << prob << " T" << RESTendl; - RESTDebug << "Coherence length: " << lCoh << " mm" << RESTendl; + RESTDebug << "Probability: " << prob.first << " T" << RESTendl; + RESTDebug << "Error: " << prob.second << " T" << RESTendl; RESTDebug << "Transmission: " << transmission << " mm" << RESTendl; - SetObservableValue("fieldAverage", fieldAverage); - SetObservableValue("probability", prob); - SetObservableValue("coherenceLength", lCoh); - SetObservableValue("transmission", transmission); - if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) fAxionEvent->PrintEvent(); /// Missing to propagate the axion to the end of magnet bore? From d3402d1d6953362927e4f6e71796df12336f6584 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 10 Feb 2024 22:13:59 +0000 Subject: [PATCH 13/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- inc/TRestAxionField.h | 25 +- inc/TRestAxionFieldPropagationProcess.h | 16 +- inc/TRestAxionMagneticField.h | 33 +-- src/TRestAxionField.cxx | 270 +++++++++++----------- src/TRestAxionFieldPropagationProcess.cxx | 21 +- src/TRestAxionMagneticField.cxx | 138 ++++++----- 6 files changed, 257 insertions(+), 246 deletions(-) diff --git a/inc/TRestAxionField.h b/inc/TRestAxionField.h index 48777a8a..aa73017a 100644 --- a/inc/TRestAxionField.h +++ b/inc/TRestAxionField.h @@ -48,9 +48,11 @@ class TRestAxionField : public TObject { /// A pointer to the magnetic field definition TRestAxionMagneticField* fMagneticField = nullptr; //! - std::pair ComputeOffResonanceIntegral(Double_t q, Double_t Gamma, Double_t accuracy, Int_t num_intervals, Int_t qawo_levels ); + std::pair ComputeOffResonanceIntegral(Double_t q, Double_t Gamma, Double_t accuracy, + Int_t num_intervals, Int_t qawo_levels); - std::pair ComputeResonanceIntegral( Double_t Gamma, Double_t accuracy, Int_t num_intervals); + std::pair ComputeResonanceIntegral(Double_t Gamma, Double_t accuracy, + Int_t num_intervals); public: void SetMagneticField(Double_t b) { fBmag = b; } @@ -89,17 +91,20 @@ class TRestAxionField : public TObject { Double_t GammaTransmissionProbability(std::vector Bmag, Double_t deltaL, Double_t Ea, Double_t ma, Double_t mg = 0, Double_t absLength = 0); - std::pair GammaTransmissionFieldMapProbability( Double_t Ea, Double_t ma, Double_t accuracy = 1.e-1, Int_t num_intervals = 100, Int_t qawo_levels = 20 ); + std::pair 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*>(params); + /// Integrand used for axion-photon probability integration + static double Integrand(double x, void* params) { + auto* data = reinterpret_cast*>(params); - TRestAxionMagneticField* field = data->first; - double gamma = data->second; + TRestAxionMagneticField* field = data->first; + double gamma = data->second; - return exp(0.5 * gamma * x) * field->GetTransversalComponentInParametricTrack(x); - } + return exp(0.5 * gamma * x) * field->GetTransversalComponentInParametricTrack(x); + } Double_t GammaTransmissionFWHM(Double_t step = 0.00001); diff --git a/inc/TRestAxionFieldPropagationProcess.h b/inc/TRestAxionFieldPropagationProcess.h index 4b41bb38..bb8267e7 100644 --- a/inc/TRestAxionFieldPropagationProcess.h +++ b/inc/TRestAxionFieldPropagationProcess.h @@ -37,17 +37,17 @@ class TRestAxionFieldPropagationProcess : public TRestAxionEventProcess { /// 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; + /// 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 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; + /// 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); + /// 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; //! diff --git a/inc/TRestAxionMagneticField.h b/inc/TRestAxionMagneticField.h index de43ded3..bc33354f 100644 --- a/inc/TRestAxionMagneticField.h +++ b/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 fConstantField; //< - /// The size of a grid element from the mesh in mm. Initially, it must be the same as the binary input data + /// The size of a grid element from the mesh in mm. Initially, it must be the same as the binary input + /// data std::vector fMeshSize; //< /// The type of the mesh used (default is cylindrical) @@ -68,19 +69,19 @@ class TRestAxionMagneticField : public TRestMetadata { std::vector fBoundMax; //< /// A vector that defines the new mesh cell volume. It will re-scale the original fMeshSize. - TVector3 fReMap = TVector3(0,0,0); //< + TVector3 fReMap = TVector3(0, 0, 0); //< /// A magnetic field volume structure to store field data and mesh. std::vector 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 start track position used to parameterize the field along a track + TVector3 fTrackStart = TVector3(0, 0, 0); //! - /// The total length of the track which defines the limit for field parameterization - Double_t fTrackLength = 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; //! @@ -114,13 +115,13 @@ class TRestAxionMagneticField : public TRestMetadata { public: void LoadMagneticVolumes(); - void ReMap( const size_t &n, const TVector3 &newMapSize ); + void ReMap(const size_t& n, const TVector3& newMapSize); + + void SetTrack(const TVector3& position, const TVector3& direction); - void SetTrack( const TVector3 &position, const TVector3 &direction ); - - Double_t GetTrackLength() const { return fTrackLength; } - TVector3 GetTrackStart() const { return fTrackStart; } - TVector3 GetTrackDirection() const { return fTrackDirection; } + 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) { @@ -155,7 +156,7 @@ class TRestAxionMagneticField : public TRestMetadata { TVector3 GetVolumeCenter(Int_t id); Double_t GetTransversalComponent(TVector3 position, TVector3 direction); - Double_t GetTransversalComponentInParametricTrack( Double_t x ); + Double_t GetTransversalComponentInParametricTrack(Double_t x); std::vector GetTransversalComponentAlongPath(TVector3 from, TVector3 to, Double_t dl = 1., Int_t Nmax = 0); diff --git a/src/TRestAxionField.cxx b/src/TRestAxionField.cxx index 76b8bc38..27020039 100644 --- a/src/TRestAxionField.cxx +++ b/src/TRestAxionField.cxx @@ -53,9 +53,10 @@ #include "TRestComplex.h" #endif -#include #include +#include + using namespace std; ClassImp(TRestAxionField); @@ -97,7 +98,8 @@ void TRestAxionField::Initialize() { double TRestAxionField::BL(Double_t Bmag, Double_t Lcoh) { Double_t lengthInMeters = Lcoh * units("m"); - Double_t tm = REST_Physics::lightSpeed / REST_Physics::naturalElectron * units("GeV")/units("eV"); // eV --> GeV + Double_t tm = + REST_Physics::lightSpeed / REST_Physics::naturalElectron * units("GeV") / units("eV"); // eV --> GeV Double_t sol = lengthInMeters * Bmag * tm; sol = sol * 1.0e-10; @@ -114,7 +116,8 @@ double TRestAxionField::BLHalfSquared(Double_t Bmag, Double_t Lcoh) // (BL/2)** { Double_t lengthInMeters = Lcoh * units("m"); - Double_t tm = REST_Physics::lightSpeed / REST_Physics::naturalElectron * units("GeV") / units( "eV"); // eV --> GeV + Double_t tm = + REST_Physics::lightSpeed / REST_Physics::naturalElectron * units("GeV") / units("eV"); // eV --> GeV Double_t sol = lengthInMeters * Bmag * tm / 2; sol = sol * sol * 1.0e-20; @@ -166,7 +169,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, Double_t Gamma = absLength; if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(fEa); // cm-1 - Double_t GammaL = Gamma * cohLength * units("cm")/units("m"); // m --> cm + Double_t GammaL = Gamma * cohLength * units("cm") / units("m"); // m --> cm if (fDebug) { std::cout << "+------------------------+" << std::endl; @@ -192,8 +195,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, double sol = (double)(MFactor * BLHalfSquared(fBmag, fLcoh) * (1 + exp(-GammaL) - 2 * exp(-GammaL / 2) * cos(phi))); - if( fDebug ) - std::cout << "Axion-photon transmission probability : " << sol << std::endl; + if (fDebug) std::cout << "Axion-photon transmission probability : " << sol << std::endl; return sol; #endif @@ -244,7 +246,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma // Default REST units are mm. We express cohLength in m. Double_t Lcoh = (Bmag.size() - 1) * deltaL; // in mm - Double_t cohLength = Lcoh * units("m"); // in m + Double_t cohLength = Lcoh * units("m"); // in m Double_t photonMass = mg; @@ -253,17 +255,18 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma Double_t fieldAverage = 0; if (Bmag.size() > 0) fieldAverage = std::accumulate(Bmag.begin(), Bmag.end(), 0.0) / Bmag.size(); - if( fDebug ) - { - std::cout << "+--------------------------------------------------------------------------+" << std::endl; - std::cout << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << std::endl; - std::cout << " Photon mass : " << photonMass << " eV" << std::endl; - std::cout << " Axion mass : " << ma << " eV" << std::endl; - std::cout << " Axion energy : " << Ea << " keV" << std::endl; - std::cout << " Lcoh : " << cohLength << " m" << std::endl; - std::cout << " Bmag average : " << fieldAverage << " T" << std::endl; - std::cout << "+--------------------------------------------------------------------------+" << std::endl; - } + if (fDebug) { + std::cout << "+--------------------------------------------------------------------------+" + << std::endl; + std::cout << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << std::endl; + std::cout << " Photon mass : " << photonMass << " eV" << std::endl; + std::cout << " Axion mass : " << ma << " eV" << std::endl; + std::cout << " Axion energy : " << Ea << " keV" << std::endl; + std::cout << " Lcoh : " << cohLength << " m" << std::endl; + std::cout << " Bmag average : " << fieldAverage << " T" << std::endl; + std::cout << "+--------------------------------------------------------------------------+" + << std::endl; + } // In vacuum if (ma == 0.0 && photonMass == 0.0) return BLHalfSquared(fieldAverage, Lcoh); @@ -274,7 +277,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma Double_t Gamma = absLength; if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 - Double_t GammaL = Gamma * cohLength * units("cm")/units("m"); // m --> cm + Double_t GammaL = Gamma * cohLength * units("cm") / units("m"); // m --> cm if (fDebug) { std::cout << "+------------------------+" << std::endl; @@ -297,28 +300,27 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma std::cout << "Exp(-GammaL) : " << exp(-GammaL) << std::endl; } - /// We integrate following the Midpoint rule method. (Other potential options : Trapezoidal, Simpsons) - /* - Double_t deltaIneV = deltaL * units("m") * REST_Physics::PhMeterIneV; - TRestComplex sum(0, 0); - for (unsigned int n = 0; n < Bmag.size() - 1; n++) { - Double_t Bmiddle = 0.5 * (Bmag[n] + Bmag[n + 1]); + /* +Double_t deltaIneV = deltaL * units("m") * REST_Physics::PhMeterIneV; +TRestComplex sum(0, 0); + for (unsigned int n = 0; n < Bmag.size() - 1; n++) { + Double_t Bmiddle = 0.5 * (Bmag[n] + Bmag[n + 1]); - Double_t lStepIneV = ((double)n + 0.5) * deltaIneV; - Double_t lStepInCm = ((double)n + 0.5) * deltaL * units("cm"); + Double_t lStepIneV = ((double)n + 0.5) * deltaIneV; + Double_t lStepInCm = ((double)n + 0.5) * deltaL * units("cm"); - TRestComplex qCgC(0.5 * Gamma * lStepInCm, -q * lStepIneV); - qCgC = TRestComplex::Exp(qCgC); + TRestComplex qCgC(0.5 * Gamma * lStepInCm, -q * lStepIneV); + qCgC = TRestComplex::Exp(qCgC); - TRestComplex integrand = Bmiddle * deltaL * qCgC; // The integrand is in T by mm + TRestComplex integrand = Bmiddle * deltaL * qCgC; // The integrand is in T by mm - sum += integrand; - } - */ + sum += integrand; + } + */ - Double_t sol = 0; - // Double_t sol = exp(-GammaL) * sum.Rho2() * BLHalfSquared(1, 1); + Double_t sol = 0; + // Double_t sol = exp(-GammaL) * sum.Rho2() * BLHalfSquared(1, 1); // Now T and mm have been recalculated in natural units using BLHalfSquared(1,1). /* @@ -326,8 +328,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma (double)(MFactor * BLHalfSquared(Bmag, Lcoh) * (1 + exp(-GammaL) - 2 * exp(-GammaL / 2) * cos(phi))); */ - if( fDebug ) - std::cout << "Axion-photon transmission probability : " << sol << std::endl; + if (fDebug) std::cout << "Axion-photon transmission probability : " << sol << std::endl; return (Double_t)sol; #endif @@ -350,30 +351,33 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma /// \note The density is for the moment homogeneous. We would need to implemnent a double integral /// to solve the problem with a density profile. /// -std::pair TRestAxionField::GammaTransmissionFieldMapProbability( Double_t Ea, Double_t ma, Double_t accuracy, Int_t num_intervals, Int_t qawo_levels ) { - - if( !fMagneticField ) - { - RESTError << "TRestAxionField::GammaTransmissionFieldMapProbability requires a magnetic field map!" << RESTendl; - RESTError << "Use TRestAxionField::AssignMagneticField method to assign one" << RESTendl; - return {0.0,0.0}; - } +std::pair TRestAxionField::GammaTransmissionFieldMapProbability(Double_t Ea, Double_t ma, + Double_t accuracy, + Int_t num_intervals, + Int_t qawo_levels) { + if (!fMagneticField) { + RESTError << "TRestAxionField::GammaTransmissionFieldMapProbability requires a magnetic field map!" + << RESTendl; + RESTError << "Use TRestAxionField::AssignMagneticField method to assign one" << RESTendl; + return {0.0, 0.0}; + } - double photonMass = 0; // Vacuum + double photonMass = 0; // Vacuum if (fBufferGas) photonMass = fBufferGas->GetPhotonMass(Ea); - if( fDebug ) - { - std::cout << "+--------------------------------------------------------------------------+" << std::endl; - std::cout << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << std::endl; - std::cout << " Photon mass : " << photonMass << " eV" << std::endl; - std::cout << " Axion mass : " << ma << " eV" << std::endl; - std::cout << " Axion energy : " << Ea << " keV" << std::endl; - std::cout << "+--------------------------------------------------------------------------+" << std::endl; - } + if (fDebug) { + std::cout << "+--------------------------------------------------------------------------+" + << std::endl; + std::cout << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << std::endl; + std::cout << " Photon mass : " << photonMass << " eV" << std::endl; + std::cout << " Axion mass : " << ma << " eV" << std::endl; + std::cout << " Axion energy : " << Ea << " keV" << std::endl; + std::cout << "+--------------------------------------------------------------------------+" + << std::endl; + } - double q = (ma * ma - photonMass * photonMass) / 2. / (Ea * units("eV") ); - q = q / REST_Physics::eVInPhMeter*units("m")/units("mm"); // mm-1 + double q = (ma * ma - photonMass * photonMass) / 2. / (Ea * units("eV")); + q = q / REST_Physics::eVInPhMeter * units("m") / units("mm"); // mm-1 double Gamma = 0; if (fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea) * units("cm") / units("mm"); // mm-1 @@ -384,12 +388,14 @@ std::pair TRestAxionField::GammaTransmissionFieldMapProbabili std::cout << " q : " << q << " eV" << std::endl; std::cout << "Gamma : " << Gamma << std::endl; std::cout << "+------------------------+" << std::endl; - } + } - if( q == 0 ) return ComputeResonanceIntegral( Gamma, accuracy, num_intervals ); - else return ComputeOffResonanceIntegral( q, Gamma, accuracy, num_intervals, qawo_levels); + if (q == 0) + return ComputeResonanceIntegral(Gamma, accuracy, num_intervals); + else + return ComputeOffResonanceIntegral(q, Gamma, accuracy, num_intervals, qawo_levels); - return {0.0,0.0}; + return {0.0, 0.0}; } /////////////////////////////////////////////// @@ -397,7 +403,7 @@ std::pair TRestAxionField::GammaTransmissionFieldMapProbabili /// equation (28) from J. Redondo and A. Ringwald, Light shinning through walls. /// https://arxiv.org/pdf/1011.3741.pdf /// -/// It integrates the Integrand function defined on the header of TRestAxionField. +/// It integrates the Integrand function defined on the header of TRestAxionField. /// /// This method uses the GSL QAG integration method. We use this method when the cosine /// function vanishes when q = 0. @@ -408,46 +414,46 @@ std::pair TRestAxionField::GammaTransmissionFieldMapProbabili /// The Gamma function should be expressed in mm-1 since the field map accessed in the integrand /// is evaluated using mm. /// -std::pair TRestAxionField::ComputeResonanceIntegral( Double_t Gamma, Double_t accuracy, Int_t num_intervals ) -{ - double reprob,rerr; +std::pair TRestAxionField::ComputeResonanceIntegral(Double_t Gamma, Double_t accuracy, + Int_t num_intervals) { + double reprob, rerr; - std::pair params = {fMagneticField, Gamma}; + std::pair params = {fMagneticField, Gamma}; - gsl_integration_workspace *workspace = gsl_integration_workspace_alloc(num_intervals); + gsl_integration_workspace* workspace = gsl_integration_workspace_alloc(num_intervals); - gsl_function F; - F.function = &Integrand; - F.params = ¶ms; + gsl_function F; + F.function = &Integrand; + F.params = ¶ms; - auto start = std::chrono::system_clock::now(); + auto start = std::chrono::system_clock::now(); - gsl_integration_qag( &F, 0, fMagneticField->GetTrackLength(), accuracy, accuracy, num_intervals, GSL_INTEG_GAUSS61, workspace, &reprob, &rerr); + gsl_integration_qag(&F, 0, fMagneticField->GetTrackLength(), accuracy, accuracy, num_intervals, + GSL_INTEG_GAUSS61, workspace, &reprob, &rerr); - auto end = std::chrono::system_clock::now(); - auto seconds = std::chrono::duration_cast(end-start); + auto end = std::chrono::system_clock::now(); + auto seconds = std::chrono::duration_cast(end - start); - double GammaL = Gamma * fMagneticField->GetTrackLength(); - double C = exp(-GammaL) * BLHalfSquared(1, 1); + double GammaL = Gamma * fMagneticField->GetTrackLength(); + double C = exp(-GammaL) * BLHalfSquared(1, 1); - double prob = C * reprob*reprob; - double proberr = 2 * C * reprob * rerr; + double prob = C * reprob * reprob; + double proberr = 2 * C * reprob * rerr; - if( fDebug ) - { - std::cout << " ---- TRestAxionField::ComputeResonanceIntegral (QAG) ----" << std::endl; - std::cout << "Gamma: " << Gamma << " mm-1" << std::endl; - std::cout << "accuracy: " << accuracy << std::endl; - std::cout << "num_intervals: " << num_intervals << std::endl; - std::cout << " -------" << std::endl; - std::cout << "Probability: " << prob << std::endl; - std::cout << "Probability error: " << proberr << std::endl; - std::cout << "Computing time: " << seconds.count() << std::endl; - std::cout << " -------" << std::endl; - } + if (fDebug) { + std::cout << " ---- TRestAxionField::ComputeResonanceIntegral (QAG) ----" << std::endl; + std::cout << "Gamma: " << Gamma << " mm-1" << std::endl; + std::cout << "accuracy: " << accuracy << std::endl; + std::cout << "num_intervals: " << num_intervals << std::endl; + std::cout << " -------" << std::endl; + std::cout << "Probability: " << prob << std::endl; + std::cout << "Probability error: " << proberr << std::endl; + std::cout << "Computing time: " << seconds.count() << std::endl; + std::cout << " -------" << std::endl; + } - std::pair sol = { prob, proberr }; - return sol; + std::pair sol = {prob, proberr}; + return sol; } /////////////////////////////////////////////// @@ -455,7 +461,7 @@ std::pair TRestAxionField::ComputeResonanceIntegral( Double_t /// equation (28) from J. Redondo and A. Ringwald, Light shinning through walls. /// https://arxiv.org/pdf/1011.3741.pdf /// -/// It integrates the Integrand function defined on the header of TRestAxionField. +/// It integrates the Integrand function defined on the header of TRestAxionField. /// /// This method uses the GSL QAWO method for oscillatory functions. That is the case when /// the q is not zero, in the off-resonance case. @@ -466,55 +472,57 @@ std::pair TRestAxionField::ComputeResonanceIntegral( Double_t /// The Gamma function should be expressed in mm-1 since the field map accessed in the integrand /// is evaluated using mm. /// -std::pair TRestAxionField::ComputeOffResonanceIntegral(Double_t q, Double_t Gamma, Double_t accuracy, Int_t num_intervals, Int_t qawo_levels ) -{ - double reprob,rerr; - double improb,imerr; +std::pair TRestAxionField::ComputeOffResonanceIntegral(Double_t q, Double_t Gamma, + Double_t accuracy, + Int_t num_intervals, + Int_t qawo_levels) { + double reprob, rerr; + double improb, imerr; - std::pair params = {fMagneticField, Gamma}; + std::pair params = {fMagneticField, Gamma}; - gsl_integration_workspace *workspace = gsl_integration_workspace_alloc(num_intervals); + gsl_integration_workspace* workspace = gsl_integration_workspace_alloc(num_intervals); - gsl_function F; - F.function = &Integrand; - F.params = ¶ms; + gsl_function F; + F.function = &Integrand; + F.params = ¶ms; - auto start = std::chrono::system_clock::now(); + auto start = std::chrono::system_clock::now(); - gsl_integration_qawo_table *wf = gsl_integration_qawo_table_alloc(q, fMagneticField->GetTrackLength(), GSL_INTEG_COSINE, qawo_levels); - gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &reprob, &rerr); + gsl_integration_qawo_table* wf = + gsl_integration_qawo_table_alloc(q, fMagneticField->GetTrackLength(), GSL_INTEG_COSINE, qawo_levels); + gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &reprob, &rerr); - gsl_integration_qawo_table_set( wf, q, fMagneticField->GetTrackLength(), GSL_INTEG_SINE ); - gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &improb, &imerr); + gsl_integration_qawo_table_set(wf, q, fMagneticField->GetTrackLength(), GSL_INTEG_SINE); + gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &improb, &imerr); - gsl_integration_qawo_table_free(wf); + gsl_integration_qawo_table_free(wf); - auto end = std::chrono::system_clock::now(); - auto seconds = std::chrono::duration_cast(end-start); + auto end = std::chrono::system_clock::now(); + auto seconds = std::chrono::duration_cast(end - start); - double GammaL = Gamma * fMagneticField->GetTrackLength(); - double C = exp(-GammaL) * BLHalfSquared(1, 1); + double GammaL = Gamma * fMagneticField->GetTrackLength(); + double C = exp(-GammaL) * BLHalfSquared(1, 1); - double prob = C * (reprob*reprob + improb*improb); - double proberr = 2 * C * TMath::Sqrt( reprob*reprob * rerr*rerr + improb*improb * imerr*imerr ); + double prob = C * (reprob * reprob + improb * improb); + double proberr = 2 * C * TMath::Sqrt(reprob * reprob * rerr * rerr + improb * improb * imerr * imerr); - if( fDebug ) - { - std::cout << " ---- TRestAxionField::ComputeOffResonanceIntegral (QAWO) ----" << std::endl; - std::cout << "Gamma: " << Gamma << " mm-1" << std::endl; - std::cout << "q: " << q << "mm-1" << std::endl; - std::cout << "accuracy: " << accuracy << std::endl; - std::cout << "num_intervals: " << num_intervals << std::endl; - std::cout << "qawo_levels: " << qawo_levels << std::endl; - std::cout << " -------" << std::endl; - std::cout << "Probability: " << prob << std::endl; - std::cout << "Probability error: " << proberr << std::endl; - std::cout << "Computing time: " << seconds.count() << std::endl; - std::cout << " -------" << std::endl; - } + if (fDebug) { + std::cout << " ---- TRestAxionField::ComputeOffResonanceIntegral (QAWO) ----" << std::endl; + std::cout << "Gamma: " << Gamma << " mm-1" << std::endl; + std::cout << "q: " << q << "mm-1" << std::endl; + std::cout << "accuracy: " << accuracy << std::endl; + std::cout << "num_intervals: " << num_intervals << std::endl; + std::cout << "qawo_levels: " << qawo_levels << std::endl; + std::cout << " -------" << std::endl; + std::cout << "Probability: " << prob << std::endl; + std::cout << "Probability error: " << proberr << std::endl; + std::cout << "Computing time: " << seconds.count() << std::endl; + std::cout << " -------" << std::endl; + } - std::pair sol = { prob, proberr }; - return sol; + std::pair sol = {prob, proberr}; + return sol; } /////////////////////////////////////////////// @@ -564,7 +572,7 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, D Double_t Gamma = absLength; if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(fEa); // cm-1 - Double_t GammaL = Gamma * cohLength * units("cm")/units("m"); + Double_t GammaL = Gamma * cohLength * units("cm") / units("m"); if (fDebug) { RESTDebug << "+------------------------+" << RESTendl; diff --git a/src/TRestAxionFieldPropagationProcess.cxx b/src/TRestAxionFieldPropagationProcess.cxx index 43e6f62b..5b25ab05 100644 --- a/src/TRestAxionFieldPropagationProcess.cxx +++ b/src/TRestAxionFieldPropagationProcess.cxx @@ -179,7 +179,7 @@ void TRestAxionFieldPropagationProcess::InitProcess() { else fBufferGasAdditionalLength = 0; - fAxionField->AssignMagneticField(fMagneticField); + fAxionField->AssignMagneticField(fMagneticField); } RESTDebug << "Axion-field : " << fAxionField << RESTendl; @@ -192,22 +192,23 @@ TRestEvent* TRestAxionFieldPropagationProcess::ProcessEvent(TRestEvent* evInput) RESTDebug << "TRestAxionFieldPropagationProcess::ProcessEvent : " << fAxionEvent->GetID() << RESTendl; - fMagneticField->SetTrack( fAxionEvent->GetPosition(), fAxionEvent->GetDirection() ); + fMagneticField->SetTrack(fAxionEvent->GetPosition(), fAxionEvent->GetDirection()); - Double_t Ea = fAxionEvent->GetEnergy(); - Double_t ma = fAxionEvent->GetMass() * units("eV"); + Double_t Ea = fAxionEvent->GetEnergy(); + Double_t ma = fAxionEvent->GetMass() * units("eV"); - std::pair prob = fAxionField->GammaTransmissionFieldMapProbability( Ea, ma, fAccuracy, fNumIntervals, fQawoLevels ); + std::pair prob = + fAxionField->GammaTransmissionFieldMapProbability(Ea, ma, fAccuracy, fNumIntervals, fQawoLevels); SetObservableValue("probability", prob.first); SetObservableValue("error", prob.second); Double_t transmission = 1; - if (fBufferGas && fBufferGasAdditionalLength > 0) { - Double_t Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 - Double_t GammaL = Gamma * fBufferGasAdditionalLength * units("cm"); - transmission = exp(-GammaL); - } + if (fBufferGas && fBufferGasAdditionalLength > 0) { + Double_t Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 + Double_t GammaL = Gamma * fBufferGasAdditionalLength * units("cm"); + transmission = exp(-GammaL); + } SetObservableValue("transmission", transmission); diff --git a/src/TRestAxionMagneticField.cxx b/src/TRestAxionMagneticField.cxx index ae5ce638..fe6ae7fd 100644 --- a/src/TRestAxionMagneticField.cxx +++ b/src/TRestAxionMagneticField.cxx @@ -919,9 +919,8 @@ void TRestAxionMagneticField::LoadMagneticVolumes() { exit(1); } - for( size_t n = 0; n < fMagneticFieldVolumes.size(); n++ ) - if( fReMap != TVector3(0,0,0) ) - ReMap( n, fReMap ); + for (size_t n = 0; n < fMagneticFieldVolumes.size(); n++) + if (fReMap != TVector3(0, 0, 0)) ReMap(n, fReMap); RESTDebug << "Finished loading magnetic volumes" << RESTendl; } @@ -1066,51 +1065,53 @@ TVector3 TRestAxionMagneticField::GetMagneticField(TVector3 pos, Bool_t showWarn /////////////////////////////////////////////// /// \brief It allows to remap the magnetic field to a larger mesh size. The new mesh -/// size granularity must be provided by argument, and each dimension must be a factor of the +/// size granularity must be provided by argument, and each dimension must be a factor of the /// present mesh size. /// -void TRestAxionMagneticField::ReMap( const size_t &n, const TVector3 &newMapSize ) -{ - if( newMapSize.X() == 0 || newMapSize.Y() == 0 || newMapSize.Z() == 0 ) - { - RESTError << "TRestAxionMagneticField::ReMap. The mesh granularity cannot be 0" << RESTendl; - RESTError << "Remapping will not take effect" << RESTendl; - return; - } - - Double_t remainder = std::fmod(newMapSize.X(), fMeshSize[n].X()) + std::fmod(newMapSize.Y(), fMeshSize[n].Y()) + std::fmod(newMapSize.Z(), fMeshSize[n].Z()); - if( remainder != 0 ) - { - RESTError << "TRestAxionMagneticField::ReMap. The field cannot be remapped."<< RESTendl; - RESTError << "The new mesh granularity must be a multiple of the existing granularity." << RESTendl; - RESTError << "Present mesh size : (" << fMeshSize[n].X() << ", " << fMeshSize[n].Y() << ", " << fMeshSize[n].Z() << ")" << RESTendl; - RESTError << "Requested mesh size : (" << newMapSize.X() << ", " << newMapSize.Y() << ", " << newMapSize.Z() << ")" << RESTendl; - RESTError << "Remapping will not take effect" << RESTendl; - return; - } - - Int_t scaleX = (Int_t) (newMapSize.X()/fMeshSize[n].X()); - Int_t scaleY = (Int_t) (newMapSize.Y()/fMeshSize[n].Y()); - Int_t scaleZ = (Int_t) (newMapSize.Z()/fMeshSize[n].Z()); - - Int_t newNodesX = (fMagneticFieldVolumes[n].mesh.GetNodesX()-1)/scaleX + 1; - Int_t newNodesY = (fMagneticFieldVolumes[n].mesh.GetNodesY()-1)/scaleY + 1; - Int_t newNodesZ = (fMagneticFieldVolumes[n].mesh.GetNodesZ()-1)/scaleZ + 1; - - for( Int_t nx = 0; nx < newNodesX; nx++ ) - for( Int_t ny = 0; ny < newNodesY; ny++ ) - for( Int_t nz = 0; nz < newNodesZ; nz++ ) - fMagneticFieldVolumes[n].field[nx][ny][nz] = fMagneticFieldVolumes[n].field[nx*scaleX][ny*scaleY][nz*scaleZ]; - - fMagneticFieldVolumes[n].mesh.SetNodes(newNodesX, newNodesY, newNodesZ); - fMagneticFieldVolumes[n].field.resize(fMagneticFieldVolumes[n].mesh.GetNodesX()); - for (unsigned int i = 0; i < fMagneticFieldVolumes[n].field.size(); i++) { - fMagneticFieldVolumes[n].field[n].resize(fMagneticFieldVolumes[n].mesh.GetNodesY()); - for (unsigned int j = 0; j < fMagneticFieldVolumes[n].field[i].size(); j++) - fMagneticFieldVolumes[n].field[i][j].resize(fMagneticFieldVolumes[n].mesh.GetNodesZ()); - } - - fMeshSize[n] = TVector3( fMeshSize[n].X() * scaleX, fMeshSize[n].Y() * scaleY, fMeshSize[n].Z() * scaleZ ); +void TRestAxionMagneticField::ReMap(const size_t& n, const TVector3& newMapSize) { + if (newMapSize.X() == 0 || newMapSize.Y() == 0 || newMapSize.Z() == 0) { + RESTError << "TRestAxionMagneticField::ReMap. The mesh granularity cannot be 0" << RESTendl; + RESTError << "Remapping will not take effect" << RESTendl; + return; + } + + Double_t remainder = std::fmod(newMapSize.X(), fMeshSize[n].X()) + + std::fmod(newMapSize.Y(), fMeshSize[n].Y()) + + std::fmod(newMapSize.Z(), fMeshSize[n].Z()); + if (remainder != 0) { + RESTError << "TRestAxionMagneticField::ReMap. The field cannot be remapped." << RESTendl; + RESTError << "The new mesh granularity must be a multiple of the existing granularity." << RESTendl; + RESTError << "Present mesh size : (" << fMeshSize[n].X() << ", " << fMeshSize[n].Y() << ", " + << fMeshSize[n].Z() << ")" << RESTendl; + RESTError << "Requested mesh size : (" << newMapSize.X() << ", " << newMapSize.Y() << ", " + << newMapSize.Z() << ")" << RESTendl; + RESTError << "Remapping will not take effect" << RESTendl; + return; + } + + Int_t scaleX = (Int_t)(newMapSize.X() / fMeshSize[n].X()); + Int_t scaleY = (Int_t)(newMapSize.Y() / fMeshSize[n].Y()); + Int_t scaleZ = (Int_t)(newMapSize.Z() / fMeshSize[n].Z()); + + Int_t newNodesX = (fMagneticFieldVolumes[n].mesh.GetNodesX() - 1) / scaleX + 1; + Int_t newNodesY = (fMagneticFieldVolumes[n].mesh.GetNodesY() - 1) / scaleY + 1; + Int_t newNodesZ = (fMagneticFieldVolumes[n].mesh.GetNodesZ() - 1) / scaleZ + 1; + + for (Int_t nx = 0; nx < newNodesX; nx++) + for (Int_t ny = 0; ny < newNodesY; ny++) + for (Int_t nz = 0; nz < newNodesZ; nz++) + fMagneticFieldVolumes[n].field[nx][ny][nz] = + fMagneticFieldVolumes[n].field[nx * scaleX][ny * scaleY][nz * scaleZ]; + + fMagneticFieldVolumes[n].mesh.SetNodes(newNodesX, newNodesY, newNodesZ); + fMagneticFieldVolumes[n].field.resize(fMagneticFieldVolumes[n].mesh.GetNodesX()); + for (unsigned int i = 0; i < fMagneticFieldVolumes[n].field.size(); i++) { + fMagneticFieldVolumes[n].field[n].resize(fMagneticFieldVolumes[n].mesh.GetNodesY()); + for (unsigned int j = 0; j < fMagneticFieldVolumes[n].field[i].size(); j++) + fMagneticFieldVolumes[n].field[i][j].resize(fMagneticFieldVolumes[n].mesh.GetNodesZ()); + } + + fMeshSize[n] = TVector3(fMeshSize[n].X() * scaleX, fMeshSize[n].Y() * scaleY, fMeshSize[n].Z() * scaleZ); } /////////////////////////////////////////////// @@ -1198,38 +1199,35 @@ std::vector TRestAxionMagneticField::GetTransversalComponentAlongPath( } /////////////////////////////////////////////// -/// \brief It initializes the field track boundaries data members of this class using a +/// \brief It initializes the field track boundaries data members of this class using a /// track position and direction so that these values can be used later on in parameterization. /// -void TRestAxionMagneticField::SetTrack( const TVector3 &position, const TVector3 &direction ) -{ - std::vector trackBounds = GetFieldBoundaries(position, direction); +void TRestAxionMagneticField::SetTrack(const TVector3& position, const TVector3& direction) { + std::vector trackBounds = GetFieldBoundaries(position, direction); - if( trackBounds.size() != 2 ) - { - fTrackStart = TVector3(0,0,0); - fTrackDirection = TVector3(0,0,0); - fTrackLength = 0; - } + if (trackBounds.size() != 2) { + fTrackStart = TVector3(0, 0, 0); + fTrackDirection = TVector3(0, 0, 0); + fTrackLength = 0; + } - fTrackStart = trackBounds[0]; - fTrackLength = (trackBounds[1] - trackBounds[0]).Mag() - 1; - fTrackDirection = (trackBounds[1] - trackBounds[0]).Unit(); + fTrackStart = trackBounds[0]; + fTrackLength = (trackBounds[1] - trackBounds[0]).Mag() - 1; + fTrackDirection = (trackBounds[1] - trackBounds[0]).Unit(); } /////////////////////////////////////////////// /// \brief It will return the transversal magnetic field component evaluated at a parametric -/// distance `x` (given by argument) for the track defined inside the class. The track will +/// distance `x` (given by argument) for the track defined inside the class. The track will /// be defined by the data members fStartTrack and fEndTrack which should be initialized by /// external means by invoking the method SetTrack( position, direction ); /// /// This method will be used by the integration method /// -Double_t TRestAxionMagneticField::GetTransversalComponentInParametricTrack( Double_t x ) { - if( x < 0 || x > fTrackLength ) - return 0; +Double_t TRestAxionMagneticField::GetTransversalComponentInParametricTrack(Double_t x) { + if (x < 0 || x > fTrackLength) return 0; - return GetTransversalComponent(fTrackStart + x * fTrackDirection, fTrackDirection); + return GetTransversalComponent(fTrackStart + x * fTrackDirection, fTrackDirection); } /////////////////////////////////////////////// @@ -1416,8 +1414,7 @@ std::vector TRestAxionMagneticField::GetFieldBoundaries(TVector3 pos, /// \brief Initialization of TRestAxionMagnetic field members through a RML file /// void TRestAxionMagneticField::InitFromConfigFile() { - - TRestMetadata::InitFromConfigFile(); + TRestMetadata::InitFromConfigFile(); auto magVolumeDef = GetElement("addMagneticVolume"); while (magVolumeDef) { @@ -1490,11 +1487,10 @@ void TRestAxionMagneticField::InitFromConfigFile() { void TRestAxionMagneticField::PrintMetadata() { TRestMetadata::PrintMetadata(); - if( fReMap != TVector3(0,0,0) ) - { - RESTMetadata << "Fields original mesh size has been remapped" << RESTendl; - RESTMetadata << " " << RESTendl; - } + if (fReMap != TVector3(0, 0, 0)) { + RESTMetadata << "Fields original mesh size has been remapped" << RESTendl; + RESTMetadata << " " << RESTendl; + } RESTMetadata << " - Number of magnetic volumes : " << GetNumberOfVolumes() << RESTendl; RESTMetadata << " ------------------------------------------------ " << RESTendl; From 062f6ce449e3c189f3facf7fa2fd5cda62ae5f3c Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sun, 11 Feb 2024 09:22:16 +0100 Subject: [PATCH 14/32] Adding warning on pipeline scripts --- pipeline/ray-tracing/axion-field/photonConversion.rml | 1 + pipeline/ray-tracing/axion-field/plots.rml | 1 + 2 files changed, 2 insertions(+) diff --git a/pipeline/ray-tracing/axion-field/photonConversion.rml b/pipeline/ray-tracing/axion-field/photonConversion.rml index 1d176dc4..9fbf6930 100644 --- a/pipeline/ray-tracing/axion-field/photonConversion.rml +++ b/pipeline/ray-tracing/axion-field/photonConversion.rml @@ -26,6 +26,7 @@ + diff --git a/pipeline/ray-tracing/axion-field/plots.rml b/pipeline/ray-tracing/axion-field/plots.rml index bdc704ff..c092fcb2 100644 --- a/pipeline/ray-tracing/axion-field/plots.rml +++ b/pipeline/ray-tracing/axion-field/plots.rml @@ -1,6 +1,7 @@ + From 6c4e496e6658964db247e09718a960b595b5f008 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sun, 11 Feb 2024 09:32:50 +0100 Subject: [PATCH 15/32] TRestAxionFieldPropagationProcess::PrintMetadata. Updating --- inc/TRestAxionFieldPropagationProcess.h | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/inc/TRestAxionFieldPropagationProcess.h b/inc/TRestAxionFieldPropagationProcess.h index bb8267e7..2acd3063 100644 --- a/inc/TRestAxionFieldPropagationProcess.h +++ b/inc/TRestAxionFieldPropagationProcess.h @@ -72,7 +72,13 @@ 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; From 78461746c86d11b4fe91d354d2856664e7c1353c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 11 Feb 2024 08:33:01 +0000 Subject: [PATCH 16/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- inc/TRestAxionFieldPropagationProcess.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/inc/TRestAxionFieldPropagationProcess.h b/inc/TRestAxionFieldPropagationProcess.h index 2acd3063..8d14dc32 100644 --- a/inc/TRestAxionFieldPropagationProcess.h +++ b/inc/TRestAxionFieldPropagationProcess.h @@ -72,13 +72,14 @@ class TRestAxionFieldPropagationProcess : public TRestAxionEventProcess { void PrintMetadata() override { BeginPrintProcess(); - RESTMetadata << "Integration parameters" << 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 << " " << 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; From bce8a7a8ebf80f1c455db39719d30fd92f68cf41 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sun, 11 Feb 2024 15:43:05 +0100 Subject: [PATCH 17/32] TRestAxionField. Recovering classic integration routine and removing MPFR --- src/TRestAxionField.cxx | 57 +++++++++-------------------------------- 1 file changed, 12 insertions(+), 45 deletions(-) diff --git a/src/TRestAxionField.cxx b/src/TRestAxionField.cxx index 27020039..6f08ae56 100644 --- a/src/TRestAxionField.cxx +++ b/src/TRestAxionField.cxx @@ -139,13 +139,6 @@ double TRestAxionField::BLHalfSquared(Double_t Bmag, Double_t Lcoh) // (BL/2)** /// The returned value is given for g_ag = 10^-10 GeV-1 /// Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, Double_t absLength) { -#ifndef USE_MPFR - RESTWarning - << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFR=ON to your REST compilation" - << RESTendl; - RESTWarning << "TRestAxionField::GammaTransmissionProbability will return 0" << RESTendl; - return 0; -#else Double_t cohLength = fLcoh * units("m"); // Default REST units are mm; Double_t photonMass = mg; @@ -198,7 +191,6 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t ma, Double_t mg, if (fDebug) std::cout << "Axion-photon transmission probability : " << sol << std::endl; return sol; -#endif } /////////////////////////////////////////////// @@ -236,14 +228,6 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t Bmag, Double_t L Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bmag, Double_t deltaL, Double_t Ea, Double_t ma, Double_t mg, Double_t absLength) { -#ifndef USE_MPFR - RESTWarning - << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFR=ON to your REST compilation" - << RESTendl; - RESTWarning << "TRestAxionField::GammaTransmissionProbability will return 0" << RESTendl; - return 0; -#else - // Default REST units are mm. We express cohLength in m. Double_t Lcoh = (Bmag.size() - 1) * deltaL; // in mm Double_t cohLength = Lcoh * units("m"); // in m @@ -301,37 +285,28 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma } /// We integrate following the Midpoint rule method. (Other potential options : Trapezoidal, Simpsons) - /* -Double_t deltaIneV = deltaL * units("m") * REST_Physics::PhMeterIneV; -TRestComplex sum(0, 0); - for (unsigned int n = 0; n < Bmag.size() - 1; n++) { - Double_t Bmiddle = 0.5 * (Bmag[n] + Bmag[n + 1]); + Double_t deltaIneV = deltaL * units("m") * REST_Physics::PhMeterIneV; + TRestComplex sum(0, 0); + for (unsigned int n = 0; n < Bmag.size() - 1; n++) { + Double_t Bmiddle = 0.5 * (Bmag[n] + Bmag[n + 1]); - Double_t lStepIneV = ((double)n + 0.5) * deltaIneV; - Double_t lStepInCm = ((double)n + 0.5) * deltaL * units("cm"); + Double_t lStepIneV = ((double)n + 0.5) * deltaIneV; + Double_t lStepInCm = ((double)n + 0.5) * deltaL * units("cm"); - TRestComplex qCgC(0.5 * Gamma * lStepInCm, -q * lStepIneV); - qCgC = TRestComplex::Exp(qCgC); + TRestComplex qCgC(0.5 * Gamma * lStepInCm, -q * lStepIneV); + qCgC = TRestComplex::Exp(qCgC); - TRestComplex integrand = Bmiddle * deltaL * qCgC; // The integrand is in T by mm + TRestComplex integrand = Bmiddle * deltaL * qCgC; // The integrand is in T by mm - sum += integrand; - } - */ + sum += integrand; + } - Double_t sol = 0; - // Double_t sol = exp(-GammaL) * sum.Rho2() * BLHalfSquared(1, 1); + Double_t sol = exp(-GammaL) * sum.Rho2() * BLHalfSquared(1, 1); // Now T and mm have been recalculated in natural units using BLHalfSquared(1,1). - /* - double sol = - (double)(MFactor * BLHalfSquared(Bmag, Lcoh) * (1 + exp(-GammaL) - 2 * exp(-GammaL / 2) * cos(phi))); - */ - if (fDebug) std::cout << "Axion-photon transmission probability : " << sol << std::endl; return (Double_t)sol; -#endif } /////////////////////////////////////////////// @@ -539,13 +514,6 @@ std::pair TRestAxionField::ComputeOffResonanceIntegral(Doubl /// The returned value is given for g_ag = 10^-10 GeV-1 /// Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, Double_t absLength) { -#ifndef USE_MPFR - RESTWarning - << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFr=ON to your REST compilation" - << RESTendl; - RESTWarning << "TRestAxionField::AxionAbsorptionProbability will return 0" << RESTendl; - return 0; -#else Double_t cohLength = fLcoh * units("m"); // Default REST units are mm; Double_t photonMass = mg; @@ -600,7 +568,6 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, D if (fDebug) RESTDebug << "Axion-photon absorption probability : " << sol << RESTendl; return sol; -#endif } /////////////////////////////////////////////// From c02a4bd3a1053ff7b738060e28bfd715afc33db3 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sun, 11 Feb 2024 17:17:00 +0100 Subject: [PATCH 18/32] TRestAxionField. Replacing TRestComplex (based on MPFR) by TComplex --- src/TRestAxionField.cxx | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/src/TRestAxionField.cxx b/src/TRestAxionField.cxx index 6f08ae56..1ad68e95 100644 --- a/src/TRestAxionField.cxx +++ b/src/TRestAxionField.cxx @@ -49,10 +49,8 @@ #include "TH1F.h" -#ifdef USE_MPFR -#include "TRestComplex.h" -#endif +#include #include #include @@ -78,10 +76,6 @@ TRestAxionField::~TRestAxionField() {} /// So that we can still calculate numbers such as : 1.0 - 1.e-30 /// void TRestAxionField::Initialize() { -#ifdef USE_MPFR - TRestComplex::SetPrecision(30); -#endif - fBufferGas = NULL; /// MOVED TO TRestAxionFieldPropagationProcess class @@ -286,17 +280,17 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma /// We integrate following the Midpoint rule method. (Other potential options : Trapezoidal, Simpsons) Double_t deltaIneV = deltaL * units("m") * REST_Physics::PhMeterIneV; - TRestComplex sum(0, 0); + TComplex sum(0, 0); for (unsigned int n = 0; n < Bmag.size() - 1; n++) { Double_t Bmiddle = 0.5 * (Bmag[n] + Bmag[n + 1]); Double_t lStepIneV = ((double)n + 0.5) * deltaIneV; Double_t lStepInCm = ((double)n + 0.5) * deltaL * units("cm"); - TRestComplex qCgC(0.5 * Gamma * lStepInCm, -q * lStepIneV); - qCgC = TRestComplex::Exp(qCgC); + TComplex qCgC(0.5 * Gamma * lStepInCm, -q * lStepIneV); + qCgC = TComplex::Exp(qCgC); - TRestComplex integrand = Bmiddle * deltaL * qCgC; // The integrand is in T by mm + TComplex integrand = Bmiddle * deltaL * qCgC; // The integrand is in T by mm sum += integrand; } @@ -352,7 +346,7 @@ std::pair TRestAxionField::GammaTransmissionFieldMapProbabil } double q = (ma * ma - photonMass * photonMass) / 2. / (Ea * units("eV")); - q = q / REST_Physics::eVInPhMeter * units("m") / units("mm"); // mm-1 + q = q * REST_Physics::PhMeterIneV * units("m") / units("mm"); // mm-1 double Gamma = 0; if (fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea) * units("cm") / units("mm"); // mm-1 From b54d36a344762ec408ed36a1c4df90e614bd1635 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 11 Feb 2024 16:18:35 +0000 Subject: [PATCH 19/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/TRestAxionField.cxx | 30 ++++++++++++++---------------- 1 file changed, 14 insertions(+), 16 deletions(-) diff --git a/src/TRestAxionField.cxx b/src/TRestAxionField.cxx index 1ad68e95..0883e6ce 100644 --- a/src/TRestAxionField.cxx +++ b/src/TRestAxionField.cxx @@ -45,16 +45,14 @@ /// #include "TRestAxionField.h" -#include - -#include "TH1F.h" - - #include +#include #include #include +#include "TH1F.h" + using namespace std; ClassImp(TRestAxionField); @@ -279,21 +277,21 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma } /// We integrate following the Midpoint rule method. (Other potential options : Trapezoidal, Simpsons) - Double_t deltaIneV = deltaL * units("m") * REST_Physics::PhMeterIneV; - TComplex sum(0, 0); - for (unsigned int n = 0; n < Bmag.size() - 1; n++) { - Double_t Bmiddle = 0.5 * (Bmag[n] + Bmag[n + 1]); + Double_t deltaIneV = deltaL * units("m") * REST_Physics::PhMeterIneV; + TComplex sum(0, 0); + for (unsigned int n = 0; n < Bmag.size() - 1; n++) { + Double_t Bmiddle = 0.5 * (Bmag[n] + Bmag[n + 1]); - Double_t lStepIneV = ((double)n + 0.5) * deltaIneV; - Double_t lStepInCm = ((double)n + 0.5) * deltaL * units("cm"); + Double_t lStepIneV = ((double)n + 0.5) * deltaIneV; + Double_t lStepInCm = ((double)n + 0.5) * deltaL * units("cm"); - TComplex qCgC(0.5 * Gamma * lStepInCm, -q * lStepIneV); - qCgC = TComplex::Exp(qCgC); + TComplex qCgC(0.5 * Gamma * lStepInCm, -q * lStepIneV); + qCgC = TComplex::Exp(qCgC); - TComplex integrand = Bmiddle * deltaL * qCgC; // The integrand is in T by mm + TComplex integrand = Bmiddle * deltaL * qCgC; // The integrand is in T by mm - sum += integrand; - } + sum += integrand; + } Double_t sol = exp(-GammaL) * sum.Rho2() * BLHalfSquared(1, 1); // Now T and mm have been recalculated in natural units using BLHalfSquared(1,1). From 3e8052f185f0ca9ee808880e5d86fb19126e9c5f Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sun, 11 Feb 2024 18:11:47 +0100 Subject: [PATCH 20/32] Adding REST_Axion_FieldIntegrationTests.C macro --- macros/REST_Axion_FieldIntegrationTests.C | 48 +++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 macros/REST_Axion_FieldIntegrationTests.C diff --git a/macros/REST_Axion_FieldIntegrationTests.C b/macros/REST_Axion_FieldIntegrationTests.C new file mode 100644 index 00000000..345402a7 --- /dev/null +++ b/macros/REST_Axion_FieldIntegrationTests.C @@ -0,0 +1,48 @@ +#include "TCanvas.h" +#include "TGraph.h" +#include "TLatex.h" +#include "TLegend.h" +#include "TLine.h" +#include "TRestAxionBufferGas.h" +#include "TRestAxionField.h" +//******************************************************************************************************* +//*** Description: This script plots the transmission probability as a function of the axion mass for a given +//*** +//*** -------------- +//*** Usage: restManager PlotResonances [ma_max=0.1] [ma_min=0] [Ea=4.2] [Bmag=2.5] [Lcoh=10000] [gasName=He] +//*** [vacuum=true] [n_ma=10000] +//*** +//*** Being all of them optional arguments. +//*** -------------- +//*** 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 prob = ax->GammaTransmissionFieldMapProbability( Ea, ma - dm, tolerance, 100, 25 ); + + return 0; +} From d279801685b1e1fd3cababe6c4cc8c5361d3fddb Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 11 Feb 2024 17:12:24 +0000 Subject: [PATCH 21/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- macros/REST_Axion_FieldIntegrationTests.C | 57 ++++++++++++----------- 1 file changed, 29 insertions(+), 28 deletions(-) diff --git a/macros/REST_Axion_FieldIntegrationTests.C b/macros/REST_Axion_FieldIntegrationTests.C index 345402a7..d7cfc0cf 100644 --- a/macros/REST_Axion_FieldIntegrationTests.C +++ b/macros/REST_Axion_FieldIntegrationTests.C @@ -17,32 +17,33 @@ //*** 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 prob = ax->GammaTransmissionFieldMapProbability( Ea, ma - dm, tolerance, 100, 25 ); - - return 0; +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 prob = + ax->GammaTransmissionFieldMapProbability(Ea, ma - dm, tolerance, 100, 25); + + return 0; } From 9b80b469821bfcc8b3ac35c2c5e694472ffd327a Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sun, 11 Feb 2024 18:29:49 +0100 Subject: [PATCH 22/32] Updating axion-field probability validation --- pipeline/ray-tracing/axion-field/Validate.C | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/pipeline/ray-tracing/axion-field/Validate.C b/pipeline/ray-tracing/axion-field/Validate.C index 199958c8..8381f676 100644 --- a/pipeline/ray-tracing/axion-field/Validate.C +++ b/pipeline/ray-tracing/axion-field/Validate.C @@ -1,10 +1,9 @@ #include #include -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,6 +33,10 @@ Int_t Validate(Double_t prob = 5.18573e-19, Double_t fieldAverage = 1.5764) { return 3; } + /* 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(); @@ -53,6 +56,7 @@ Int_t Validate(Double_t prob = 5.18573e-19, Double_t fieldAverage = 1.5764) { std::cout << "Wrong average field!" << std::endl; return 5; } + */ delete run; From 59582684db7ee20d5717ffcf0eaabb7db1f508a7 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sun, 11 Feb 2024 18:33:01 +0100 Subject: [PATCH 23/32] pipeline/photonConversion.rml updating to the new TRestAxionFieldPropagationProcess --- pipeline/ray-tracing/axion-field/photonConversion.rml | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/pipeline/ray-tracing/axion-field/photonConversion.rml b/pipeline/ray-tracing/axion-field/photonConversion.rml index 9fbf6930..f5653ee0 100644 --- a/pipeline/ray-tracing/axion-field/photonConversion.rml +++ b/pipeline/ray-tracing/axion-field/photonConversion.rml @@ -26,8 +26,13 @@ - - + + + + + + + From 1cb4cb2f494cbc670e296b3b6ec7fb4d1dd2423d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 11 Feb 2024 17:34:28 +0000 Subject: [PATCH 24/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- pipeline/ray-tracing/axion-field/Validate.C | 42 ++++++++++----------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/pipeline/ray-tracing/axion-field/Validate.C b/pipeline/ray-tracing/axion-field/Validate.C index 8381f676..e0785e40 100644 --- a/pipeline/ray-tracing/axion-field/Validate.C +++ b/pipeline/ray-tracing/axion-field/Validate.C @@ -33,30 +33,30 @@ Int_t Validate(Double_t prob = 2.34295e-21) { return 3; } - /* 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"); + /* 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; From 05c1cc0d0bcb07149a82e667d6c004662b826107 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sun, 11 Feb 2024 19:00:13 +0100 Subject: [PATCH 25/32] pipeline/ray-tracing/full-chain/ValidateChain.C updating --- pipeline/ray-tracing/full-chain/ValidateChain.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pipeline/ray-tracing/full-chain/ValidateChain.C b/pipeline/ray-tracing/full-chain/ValidateChain.C index b67bb41d..264d3e05 100644 --- a/pipeline/ray-tracing/full-chain/ValidateChain.C +++ b/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; } From 15451d072cceab32a20783215b8837b57efe4c33 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sun, 11 Feb 2024 19:38:43 +0100 Subject: [PATCH 26/32] pipeline/ray-tracing/full-chain/ValidateChain.C updating conditions --- pipeline/ray-tracing/full-chain/ValidateChain.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pipeline/ray-tracing/full-chain/ValidateChain.C b/pipeline/ray-tracing/full-chain/ValidateChain.C index 264d3e05..1769fad8 100644 --- a/pipeline/ray-tracing/full-chain/ValidateChain.C +++ b/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 < 2.01953e-22 || integral > 1.99953e-22) { + if (integral > 2.01953e-22 || integral < 1.99953e-22) { std::cout << "Axion-photon probability is not within the expected range!" << std::endl; return 3; } From 5675a82d7ee65812e1ec4a0fa31bf88074bd6fcc Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Sun, 11 Feb 2024 21:00:36 +0100 Subject: [PATCH 27/32] REST_Axion_FieldIntegrationTests.C Updating documentation --- macros/REST_Axion_FieldIntegrationTests.C | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/macros/REST_Axion_FieldIntegrationTests.C b/macros/REST_Axion_FieldIntegrationTests.C index d7cfc0cf..6f3788e2 100644 --- a/macros/REST_Axion_FieldIntegrationTests.C +++ b/macros/REST_Axion_FieldIntegrationTests.C @@ -6,17 +6,18 @@ #include "TRestAxionBufferGas.h" #include "TRestAxionField.h" //******************************************************************************************************* -//*** Description: This script plots the transmission probability as a function of the axion mass for a given +//*** 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. //*** -//*** -------------- -//*** Usage: restManager PlotResonances [ma_max=0.1] [ma_min=0] [Ea=4.2] [Bmag=2.5] [Lcoh=10000] [gasName=He] -//*** [vacuum=true] [n_ma=10000] +//*** The macro sets the TRestAxionField under debug mode to print the different results on screen. //*** -//*** Being all of them optional arguments. //*** -------------- +//*** 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 From 4c74ea91655cc3905984f654e1cd5ff08f9c8343 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Mon, 12 Feb 2024 10:38:22 +0100 Subject: [PATCH 28/32] TRestAxionFieldPropagationProcess. It will now remap the magnetic field --- scripts/{xlsMentinkToBin.py => csvMentinkToBin.py} | 0 src/TRestAxionFieldPropagationProcess.cxx | 5 +++++ 2 files changed, 5 insertions(+) rename scripts/{xlsMentinkToBin.py => csvMentinkToBin.py} (100%) diff --git a/scripts/xlsMentinkToBin.py b/scripts/csvMentinkToBin.py similarity index 100% rename from scripts/xlsMentinkToBin.py rename to scripts/csvMentinkToBin.py diff --git a/src/TRestAxionFieldPropagationProcess.cxx b/src/TRestAxionFieldPropagationProcess.cxx index 5b25ab05..6ca499e4 100644 --- a/src/TRestAxionFieldPropagationProcess.cxx +++ b/src/TRestAxionFieldPropagationProcess.cxx @@ -170,6 +170,10 @@ void TRestAxionFieldPropagationProcess::InitProcess() { exit(0); } + for( size_t n = 0; n < fMagneticField->GetNumberOfVolumes(); n++ ) + fMagneticField->ReMap( n, fReMap ); + + if (!fAxionField) { fAxionField = new TRestAxionField(); @@ -180,6 +184,7 @@ void TRestAxionFieldPropagationProcess::InitProcess() { fBufferGasAdditionalLength = 0; fAxionField->AssignMagneticField(fMagneticField); + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) fMagneticField->PrintMetadata(); } RESTDebug << "Axion-field : " << fAxionField << RESTendl; From 3951cb45cda15a0bb9909d7a7f38d738cbe97d8c Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Mon, 12 Feb 2024 10:39:13 +0100 Subject: [PATCH 29/32] Uploading csvMentingToBin.py used to produce binary magnetic field table --- scripts/csvMentinkToBin.py | 89 +++++++++++++++++--------------------- 1 file changed, 39 insertions(+), 50 deletions(-) diff --git a/scripts/csvMentinkToBin.py b/scripts/csvMentinkToBin.py index 8e6a8d1e..c580338d 100755 --- a/scripts/csvMentinkToBin.py +++ b/scripts/csvMentinkToBin.py @@ -14,23 +14,34 @@ import matplotlib.pyplot as plt import pandas as pd +# We insert a cut-off on the field z-axis (in meters) +zCutOffLow = -6 +zCutOffHigh = 6 + # 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 -zCenter = 0.7975 +yCenter = 0.8 +zCenter = 0 print("Starting to read") # loading the data into a matrix (xyzBdata) -file = r"../../../data/magneticField/Mentink_202401.xlsx" -df = pd.read_excel(file) +file = r"../data/magneticField/Mentink_20240124.txt" +df = pd.read_csv(file, comment="%", sep='\t',header=None) +print( df.head() ) -print(df[1:5]) +print( "---" ) +print(df[0:5]) +print( "---" ) print("Translating to matrix") -xyzBdata = df.as_matrix(columns=df.columns[0:]) +xyzBdata = df.values +print( "#####" ) print(xyzBdata[0][0:6]) +print(xyzBdata[1][0:6]) +print(xyzBdata[2][0:6]) +print( "#####" ) print(len(xyzBdata)) @@ -52,20 +63,7 @@ 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" - ) + 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]: @@ -118,35 +116,26 @@ 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[1] - yCenter), - 1000 * (x[2] - zCenter), - 1000 * (x[0] - xCenter), - x[4], - x[5], - x[3], - ] - - # 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 z-axis, when we change the sign of z-axis we must change the sign of Bz. - if y[2] > 0: - y[2] = -y[2] - y[5] = -y[5] - 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]) + + if( x[2] <= zCutOffHigh and x[2] > zCutOffLow): + 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]) fbin.close() print("Lines writen : " + str(count)) From 027b54074b86866c6e5de4bce85440fe0316b24b Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 12 Feb 2024 09:49:05 +0000 Subject: [PATCH 30/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- scripts/csvMentinkToBin.py | 67 ++++++++++++++--------- src/TRestAxionFieldPropagationProcess.cxx | 7 +-- 2 files changed, 43 insertions(+), 31 deletions(-) diff --git a/scripts/csvMentinkToBin.py b/scripts/csvMentinkToBin.py index c580338d..eee90c6f 100755 --- a/scripts/csvMentinkToBin.py +++ b/scripts/csvMentinkToBin.py @@ -27,21 +27,21 @@ print("Starting to read") # loading the data into a matrix (xyzBdata) file = r"../data/magneticField/Mentink_20240124.txt" -df = pd.read_csv(file, comment="%", sep='\t',header=None) -print( df.head() ) +df = pd.read_csv(file, comment="%", sep="\t", header=None) +print(df.head()) -print( "---" ) +print("---") print(df[0:5]) -print( "---" ) +print("---") print("Translating to matrix") xyzBdata = df.values -print( "#####" ) +print("#####") print(xyzBdata[0][0:6]) print(xyzBdata[1][0:6]) print(xyzBdata[2][0:6]) -print( "#####" ) +print("#####") print(len(xyzBdata)) @@ -63,7 +63,20 @@ 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") + 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]: @@ -117,25 +130,25 @@ # 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. - if( x[2] <= zCutOffHigh and x[2] > zCutOffLow): - 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]) + if x[2] <= zCutOffHigh and x[2] > zCutOffLow: + 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]) fbin.close() print("Lines writen : " + str(count)) diff --git a/src/TRestAxionFieldPropagationProcess.cxx b/src/TRestAxionFieldPropagationProcess.cxx index 6ca499e4..0de32071 100644 --- a/src/TRestAxionFieldPropagationProcess.cxx +++ b/src/TRestAxionFieldPropagationProcess.cxx @@ -170,9 +170,7 @@ void TRestAxionFieldPropagationProcess::InitProcess() { exit(0); } - for( size_t n = 0; n < fMagneticField->GetNumberOfVolumes(); n++ ) - fMagneticField->ReMap( n, fReMap ); - + for (size_t n = 0; n < fMagneticField->GetNumberOfVolumes(); n++) fMagneticField->ReMap(n, fReMap); if (!fAxionField) { fAxionField = new TRestAxionField(); @@ -184,7 +182,8 @@ void TRestAxionFieldPropagationProcess::InitProcess() { fBufferGasAdditionalLength = 0; fAxionField->AssignMagneticField(fMagneticField); - if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) fMagneticField->PrintMetadata(); + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) + fMagneticField->PrintMetadata(); } RESTDebug << "Axion-field : " << fAxionField << RESTendl; From 42c3d8c97d3a3a61193222e2237ce2e5764aa9cc Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Mon, 12 Feb 2024 13:53:43 +0100 Subject: [PATCH 31/32] Updating data submodule --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 8a4d3c0f..73e585c0 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 8a4d3c0f544df6e9e80046c336391e10a2a5de2a +Subproject commit 73e585c077745dd0fee21ec05898f08122382ee9 From a6c5786026c704d83558094c3639e710e8479ae8 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Mon, 12 Feb 2024 14:10:45 +0100 Subject: [PATCH 32/32] Update inc/TRestAxionField.h MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Francisco Rodríguez Candón --- inc/TRestAxionField.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inc/TRestAxionField.h b/inc/TRestAxionField.h index aa73017a..76be08b6 100644 --- a/inc/TRestAxionField.h +++ b/inc/TRestAxionField.h @@ -72,7 +72,7 @@ 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 + /// It assigns a magnetic field to the calculation void AssignMagneticField(TRestAxionMagneticField* mField) { fMagneticField = mField; } /// It assigns a gas buffer medium to the calculation