diff --git a/data b/data index 8a4d3c0f..73e585c0 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 8a4d3c0f544df6e9e80046c336391e10a2a5de2a +Subproject commit 73e585c077745dd0fee21ec05898f08122382ee9 diff --git a/inc/TRestAxionField.h b/inc/TRestAxionField.h index 83cd950c..76be08b6 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,16 @@ class TRestAxionField : public TObject { void Initialize(); /// A pointer to the buffer gas definition - TRestAxionBufferGas* fBufferGas = NULL; //! + TRestAxionBufferGas* fBufferGas = nullptr; //! + + /// A pointer to the magnetic field definition + TRestAxionMagneticField* fMagneticField = nullptr; //! + + std::pair 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; } @@ -62,6 +72,9 @@ class TRestAxionField : public TObject { /// It assigns a gas buffer medium to the calculation void AssignBufferGas(TRestAxionBufferGas* buffGas) { fBufferGas = buffGas; } + /// It assigns a magnetic field to the calculation + void AssignMagneticField(TRestAxionMagneticField* mField) { fMagneticField = mField; } + /// It assigns a gas buffer medium to the calculation void SetBufferGas(TRestAxionBufferGas* buffGas) { fBufferGas = buffGas; } @@ -78,6 +91,21 @@ class TRestAxionField : public TObject { Double_t GammaTransmissionProbability(std::vector 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); + + /// 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", diff --git a/inc/TRestAxionFieldPropagationProcess.h b/inc/TRestAxionFieldPropagationProcess.h index c12e3ac8..8d14dc32 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; //! @@ -63,7 +72,14 @@ class TRestAxionFieldPropagationProcess : public TRestAxionEventProcess { void PrintMetadata() override { BeginPrintProcess(); - RESTMetadata << "Integration step length : " << fIntegrationStep << " mm" << RESTendl; + RESTMetadata << "Integration parameters" << RESTendl; + RESTMetadata << "- Integration accuracy/tolerance : " << fAccuracy << RESTendl; + RESTMetadata << "- Max number of integration intervals : " << fNumIntervals << RESTendl; + RESTMetadata << "- Number of QAWO levels : " << fQawoLevels << RESTendl; + RESTMetadata << " " << RESTendl; + RESTMetadata << "Field re-mapping size : (" << fReMap.X() << ", " << fReMap.Y() << ", " << fReMap.Z() + << ")" << RESTendl; + RESTMetadata << " " << RESTendl; RESTMetadata << "Buffer gas additional length : " << fBufferGasAdditionalLength * units("m") << " m" << RESTendl; diff --git a/inc/TRestAxionMagneticField.h b/inc/TRestAxionMagneticField.h index dfaf80e3..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 + /// 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,9 +68,21 @@ 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; //! + /// The start track position used to parameterize the field along a track + TVector3 fTrackStart = TVector3(0, 0, 0); //! + + /// The track direction used to parameterize the field along a track + TVector3 fTrackDirection = TVector3(0, 0, 0); //! + + /// The total length of the track which defines the limit for field parameterization + Double_t fTrackLength = 0; //! + /// A helper histogram to plot the field TH2D* fHisto; //! @@ -102,6 +115,14 @@ class TRestAxionMagneticField : public TRestMetadata { public: void LoadMagneticVolumes(); + void ReMap(const size_t& n, const TVector3& newMapSize); + + void SetTrack(const TVector3& position, const TVector3& direction); + + Double_t GetTrackLength() const { return fTrackLength; } + TVector3 GetTrackStart() const { return fTrackStart; } + TVector3 GetTrackDirection() const { return fTrackDirection; } + /// It returns true if no magnetic field map was loaded for that volume Bool_t IsFieldConstant(Int_t id) { if (GetMagneticVolume(id)) return GetMagneticVolume(id)->field.size() == 0; @@ -135,6 +156,7 @@ class TRestAxionMagneticField : public TRestMetadata { TVector3 GetVolumeCenter(Int_t id); Double_t GetTransversalComponent(TVector3 position, TVector3 direction); + Double_t GetTransversalComponentInParametricTrack(Double_t x); std::vector GetTransversalComponentAlongPath(TVector3 from, TVector3 to, Double_t dl = 1., Int_t Nmax = 0); diff --git a/macros/REST_Axion_FieldIntegrationTests.C b/macros/REST_Axion_FieldIntegrationTests.C new file mode 100644 index 00000000..6f3788e2 --- /dev/null +++ b/macros/REST_Axion_FieldIntegrationTests.C @@ -0,0 +1,50 @@ +#include "TCanvas.h" +#include "TGraph.h" +#include "TLatex.h" +#include "TLegend.h" +#include "TLine.h" +#include "TRestAxionBufferGas.h" +#include "TRestAxionField.h" +//******************************************************************************************************* +//*** Description: This script will launch the integration of the axion-field with given parameters. +//*** It allows to test different magnetic field cell sizes, for a given mass that can be off-resonance +//*** for dm different from zero, and a given maximum tolerance or error for the integration routine. +//*** +//*** The macro sets the TRestAxionField under debug mode to print the different results on screen. +//*** +//*** -------------- +//*** Usage: restManager FieldIntegrationTests [sX=10] [sX=10] [sZ=10] [dm=0.01] [tolerance=0.1] [Ea=4.2] +//*** -------------- +//*** +//*** Author: Javier Galan +//******************************************************************************************************* +int REST_Axion_FieldIntegrationTests(Double_t sX = 10, Double_t sY = 10, Double_t sZ = 50, Double_t dm = 0.01, + Double_t tolerance = 0.1, Double_t Ea = 4.2) { + /// Setting up magnetic field and track to evaluate + TRestAxionMagneticField magneticField("fields.rml", "babyIAXO_HD"); + + for (size_t n = 0; n < magneticField.GetNumberOfVolumes(); n++) + magneticField.ReMap(n, TVector3(sX, sY, sZ)); + magneticField.SetTrack(TVector3(-110, -110, -11000), TVector3(0.01, 0.01, 1)); + magneticField.PrintMetadata(); + + /// Setting up the gas + TRestAxionBufferGas* gas = new TRestAxionBufferGas(); + gas->SetGasDensity("He", 3.267069078540181e-10); + gas->PrintMetadata(); + Double_t ma = gas->GetPhotonMass(Ea); + + /// Setting up the axion-field and assign gas and magnetic field. + TRestAxionField* ax = new TRestAxionField(); + ax->AssignBufferGas(gas); + ax->AssignMagneticField(&magneticField); + ax->SetAxionEnergy(Ea); + + /// Enable debugging mode in axion-field + ax->SetDebug(true); + + std::pair prob = + ax->GammaTransmissionFieldMapProbability(Ea, ma - dm, tolerance, 100, 25); + + return 0; +} diff --git a/pipeline/ray-tracing/axion-field/Validate.C b/pipeline/ray-tracing/axion-field/Validate.C index 199958c8..e0785e40 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,25 +33,30 @@ Int_t Validate(Double_t prob = 5.18573e-19, Double_t fieldAverage = 1.5764) { return 3; } - run->GetAnalysisTree()->Draw("axionPhoton_fieldAverage", "axionPhoton_fieldAverage"); + /* We do not add the field average observable anymore at TRestAxionFieldPropagationProcess + * We would need to create a new process that does this, since it is computationally + * not negligible + * +run->GetAnalysisTree()->Draw("axionPhoton_fieldAverage", "axionPhoton_fieldAverage"); - TH1D* h2 = (TH1D*)run->GetAnalysisTree()->GetHistogram(); - if (h2 == nullptr) { - std::cout << "Problems generating field average histogram" << std::endl; - return 4; - } +TH1D* h2 = (TH1D*)run->GetAnalysisTree()->GetHistogram(); +if (h2 == nullptr) { + std::cout << "Problems generating field average histogram" << std::endl; + return 4; +} - Double_t field = h2->Integral() / run->GetEntries(); - std::cout << "Average magnetic field: " << field << " T" << std::endl; - std::cout << "Expected field: " << fieldAverage << std::endl; - std::cout << "Tolerance: " << fieldTolerance << std::endl; - std::cout << "Low : " << fieldAverage - fieldTolerance << " high: " << fieldAverage + fieldTolerance - << std::endl; +Double_t field = h2->Integral() / run->GetEntries(); +std::cout << "Average magnetic field: " << field << " T" << std::endl; +std::cout << "Expected field: " << fieldAverage << std::endl; +std::cout << "Tolerance: " << fieldTolerance << std::endl; +std::cout << "Low : " << fieldAverage - fieldTolerance << " high: " << fieldAverage + fieldTolerance + << std::endl; - if (field < fieldAverage - fieldTolerance || field > fieldAverage + fieldTolerance) { - std::cout << "Wrong average field!" << std::endl; - return 5; - } +if (field < fieldAverage - fieldTolerance || field > fieldAverage + fieldTolerance) { + std::cout << "Wrong average field!" << std::endl; + return 5; +} + */ delete run; diff --git a/pipeline/ray-tracing/axion-field/photonConversion.rml b/pipeline/ray-tracing/axion-field/photonConversion.rml index 1d176dc4..f5653ee0 100644 --- a/pipeline/ray-tracing/axion-field/photonConversion.rml +++ b/pipeline/ray-tracing/axion-field/photonConversion.rml @@ -26,7 +26,13 @@ - + + + + + + + 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 @@ + diff --git a/pipeline/ray-tracing/full-chain/ValidateChain.C b/pipeline/ray-tracing/full-chain/ValidateChain.C index b67bb41d..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 < 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; } diff --git a/scripts/csvMentinkToBin.py b/scripts/csvMentinkToBin.py new file mode 100755 index 00000000..eee90c6f --- /dev/null +++ b/scripts/csvMentinkToBin.py @@ -0,0 +1,154 @@ +#!/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 + +# 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.8 +zCenter = 0 + +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()) + +print("---") +print(df[0:5]) +print("---") + +print("Translating to matrix") +xyzBdata = df.values + +print("#####") +print(xyzBdata[0][0:6]) +print(xyzBdata[1][0:6]) +print(xyzBdata[2][0:6]) +print("#####") + +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. + + 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/scripts/xlsToBin.py b/scripts/xlsBikovskiyToBin.py similarity index 100% rename from scripts/xlsToBin.py rename to scripts/xlsBikovskiyToBin.py 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; diff --git a/src/TRestAxionField.cxx b/src/TRestAxionField.cxx index fae12e99..0883e6ce 100644 --- a/src/TRestAxionField.cxx +++ b/src/TRestAxionField.cxx @@ -45,16 +45,14 @@ /// #include "TRestAxionField.h" +#include #include - -#include "TH1F.h" - -#ifdef USE_MPFR -#include "TRestComplex.h" -#endif +#include #include +#include "TH1F.h" + using namespace std; ClassImp(TRestAxionField); @@ -76,10 +74,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 @@ -94,9 +88,10 @@ 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 +106,10 @@ 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; @@ -135,17 +131,9 @@ 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 - mpfr::mpreal axionMass = ma; - mpfr::mpreal cohLength = fLcoh / 1000.; // Default REST units are mm; - - mpfr::mpreal photonMass = mg; + Double_t cohLength = fLcoh * units("m"); // Default REST units are mm; + + Double_t photonMass = mg; if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(fEa); @@ -160,42 +148,41 @@ 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 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 * 100; + Double_t 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; + Double_t 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 } /////////////////////////////////////////////// @@ -233,100 +220,276 @@ 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 - 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 / 1000.; // in m + Double_t cohLength = Lcoh * units("m"); // in m - mpfr::mpreal photonMass = mg; + Double_t photonMass = mg; if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(Ea); 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 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 * 100; + Double_t 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; + Double_t 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) - TRestComplex sum(0, 0); + 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 / 10.; - - TRestComplex qC(0, -q * lStepIneV); - qC = TRestComplex::Exp(qC); + Double_t lStepInCm = ((double)n + 0.5) * deltaL * units("cm"); - TRestComplex gC(0.5 * Gamma * lStepInCm, 0); - gC = TRestComplex::Exp(gC); + TComplex qCgC(0.5 * Gamma * lStepInCm, -q * lStepIneV); + qCgC = TComplex::Exp(qCgC); - TRestComplex integrand = Bmiddle * deltaL * gC * qC; // The integrand is in T by mm + TComplex integrand = Bmiddle * deltaL * qCgC; // The integrand is in T by mm sum += integrand; } - mpfr::mpreal 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))); - */ - - 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::PhMeterIneV * 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; } /////////////////////////////////////////////// @@ -343,23 +506,15 @@ Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bma /// 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::GammaTransmissionProbability will return 0" << RESTendl; - return 0; -#else - mpfr::mpreal axionMass = ma; - mpfr::mpreal cohLength = fLcoh / 1000.; // Default REST units are mm; - - mpfr::mpreal photonMass = mg; + Double_t cohLength = fLcoh * units("m"); // Default REST units are mm; + + Double_t photonMass = mg; if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(fEa); 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 +526,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 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 * 100; + Double_t GammaL = Gamma * cohLength * units("cm") / units("m"); if (fDebug) { RESTDebug << "+------------------------+" << RESTendl; @@ -390,7 +545,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) { @@ -405,7 +560,6 @@ Double_t TRestAxionField::AxionAbsorptionProbability(Double_t ma, Double_t mg, D if (fDebug) RESTDebug << "Axion-photon absorption probability : " << sol << RESTendl; return sol; -#endif } /////////////////////////////////////////////// diff --git a/src/TRestAxionFieldPropagationProcess.cxx b/src/TRestAxionFieldPropagationProcess.cxx index e7de98b3..0de32071 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,13 +165,13 @@ 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); } + for (size_t n = 0; n < fMagneticField->GetNumberOfVolumes(); n++) fMagneticField->ReMap(n, fReMap); + if (!fAxionField) { fAxionField = new TRestAxionField(); @@ -178,10 +180,15 @@ void TRestAxionFieldPropagationProcess::InitProcess() { fAxionField->AssignBufferGas(fBufferGas); else fBufferGasAdditionalLength = 0; + + fAxionField->AssignMagneticField(fMagneticField); + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) + fMagneticField->PrintMetadata(); } RESTDebug << "Axion-field : " << fAxionField << RESTendl; RESTDebug << "Buffer gas : " << fBufferGas << RESTendl; + RESTDebug << "Magnetic field : " << fMagneticField << RESTendl; } TRestEvent* TRestAxionFieldPropagationProcess::ProcessEvent(TRestEvent* evInput) { @@ -189,51 +196,31 @@ 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(); - - 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? diff --git a/src/TRestAxionMagneticField.cxx b/src/TRestAxionMagneticField.cxx index 08473768..fe6ae7fd 100644 --- a/src/TRestAxionMagneticField.cxx +++ b/src/TRestAxionMagneticField.cxx @@ -918,6 +918,10 @@ 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; } @@ -1059,6 +1063,57 @@ 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(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); +} + /////////////////////////////////////////////// /// \brief It returns the corresponding volume index at the given position. If not found it will return /// -1. @@ -1143,6 +1198,38 @@ 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`. @@ -1327,7 +1414,7 @@ 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) { @@ -1400,6 +1487,11 @@ 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++) {