Merge pull request #92 from rest-for-physics/jgalan_updates

Fixing argument not being passed by reference. Strongly enhances field computation time!

Author: LouisHelary

inc/TRestAxionMagneticField.h

@@ -95,7 +95,7 @@ class TRestAxionMagneticField : public TRestMetadata {
 
     void LoadMagneticFieldData(MagneticFieldVolume& mVol, std::vector<std::vector<Float_t>> data);
 
-    TVector3 GetMagneticVolumeNode(MagneticFieldVolume mVol, TVector3 pos);
+    TVector3 GetMagneticVolumeNode(size_t id, TVector3 pos);
 
     /// \brief This private method returns true if the magnetic field volumes loaded are the same as
     /// the volumes defined.
@@ -161,6 +161,9 @@ class TRestAxionMagneticField : public TRestMetadata {
     std::vector<Double_t> GetTransversalComponentAlongPath(TVector3 from, TVector3 to, Double_t dl = 1.,
                                                            Int_t Nmax = 0);
 
+    std::vector<Double_t> GetComponentAlongPath(Int_t axis, TVector3 from, TVector3 to, Double_t dl = 1.,
+                                                Int_t Nmax = 0);
+
     Double_t GetTransversalFieldAverage(TVector3 from, TVector3 to, Double_t dl = 1., Int_t Nmax = 0);
 
     TVector3 GetFieldAverageTransverseVector(TVector3 from, TVector3 to, Double_t dl = 10., Int_t Nmax = 0);
@@ -170,6 +173,8 @@ class TRestAxionMagneticField : public TRestMetadata {
 
     TCanvas* DrawTracks(TVector3 vanishingPoint, Int_t divisions, Int_t volId = 0);
 
+    TCanvas* DrawComponents(Int_t volId = 0);
+
     void PrintMetadata();
 
     TRestAxionMagneticField();

macros/REST_Axion_FieldIntegrationTests.C

@@ -21,7 +21,7 @@
 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");
+    TRestAxionMagneticField magneticField("fields.rml", "babyIAXO_2024");
 
     for (size_t n = 0; n < magneticField.GetNumberOfVolumes(); n++)
         magneticField.ReMap(n, TVector3(sX, sY, sZ));

src/TRestAxionField.cxx

@@ -210,6 +210,7 @@
 
 #include <TComplex.h>
 #include <TVectorD.h>
+#include <gsl/gsl_errno.h>
 #include <gsl/gsl_integration.h>
 
 #include <numeric>
@@ -482,13 +483,17 @@ std::pair<Double_t, Double_t> TRestAxionField::GammaTransmissionFieldMapProbabil
                                                                                     Double_t accuracy,
                                                                                     Int_t num_intervals,
                                                                                     Int_t qawo_levels) {
+    gsl_set_error_handler_off();
+
     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};
     }
 
+    if (fMagneticField->GetTrackLength() <= 0) return {0.0, 0.0};
+
     double photonMass = 0;  // Vacuum
     if (fBufferGas) photonMass = fBufferGas->GetPhotonMass(Ea);
 
@@ -555,8 +560,10 @@ std::pair<Double_t, Double_t> TRestAxionField::ComputeResonanceIntegral(Double_t
 
     auto start = std::chrono::system_clock::now();
 
-    gsl_integration_qag(&F, 0, fMagneticField->GetTrackLength(), accuracy, accuracy, num_intervals,
-                        GSL_INTEG_GAUSS61, workspace, &reprob, &rerr);
+    int status = gsl_integration_qag(&F, 0, fMagneticField->GetTrackLength(), accuracy, accuracy,
+                                     num_intervals, GSL_INTEG_GAUSS61, workspace, &reprob, &rerr);
+
+    if (status > 0) return {0, status};
 
     auto end = std::chrono::system_clock::now();
     auto seconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
@@ -618,12 +625,19 @@ std::pair<Double_t, Double_t> TRestAxionField::ComputeOffResonanceIntegral(Doubl
 
     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);
+    int status =
+        gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &reprob, &rerr);
+    if (status > 0) {
+        gsl_integration_qawo_table_free(wf);
+        return {0, status};
+    }
 
     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);
+    status = gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &improb, &imerr);
+    if (status > 0) {
+        gsl_integration_qawo_table_free(wf);
+        return {0, status};
+    }
 
     auto end = std::chrono::system_clock::now();
     auto seconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

src/TRestAxionMagneticField.cxx

@@ -584,6 +584,65 @@ TCanvas* TRestAxionMagneticField::DrawHistogram(TString projection, TString Bcom
     return fCanvas;
 }
 
+///////////////////////////////////////////////
+/// \brief A method that creates a canvas where tracks traversing the magnetic volume
+/// are drawm together with their corresponding field intensity profile along the Z-axis.
+///
+TCanvas* TRestAxionMagneticField::DrawComponents(Int_t volId) {
+    Int_t divisions = 100;
+    if (fCanvas != NULL) {
+        delete fCanvas;
+        fCanvas = NULL;
+    }
+    fCanvas = new TCanvas("fCanvas", "", 1600, 600);
+
+    TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97);
+    //   pad1->Divide(2, 1);
+    pad1->Draw();
+    pad1->cd();
+
+    Int_t n = 0;
+    Double_t genPositionZ = fPositions[volId][2] - fBoundMax[volId].Z() - 2000;
+    Double_t finalPositionZ = fPositions[volId][2] + fBoundMax[volId].Z() + 2000;
+    TVector3 position(0, 0, genPositionZ);
+    TVector3 direction(0, 0, 1);
+
+    RESTDebug << RESTendl;
+    RESTDebug << "Start moving along" << RESTendl;
+    RESTDebug << "++++++++++++++++++" << RESTendl;
+
+    TGraph* fieldGr = new TGraph();
+    Double_t posZ = fPositions[volId][2] - fBoundMax[volId].Z() - 10;
+    Double_t delta = fBoundMax[volId][2] * 2. / divisions;
+
+    while (posZ <= fPositions[volId][2] + fBoundMax[volId].Z()) {
+        TVector3 posAlongAxis = TVector3(fPositions[volId][0], fPositions[volId][1], posZ);
+
+        position = MoveToPlane(position, direction, TVector3(0, 0, 1), posAlongAxis);
+        Double_t fieldY = this->GetMagneticField(position).Y();
+
+        fieldGr->SetPoint(fieldGr->GetN(), posZ, fieldY);
+
+        posZ += delta;
+    }
+
+    fieldGr->SetLineWidth(3);
+    fieldGr->SetLineColor(38 + n);
+    fieldGr->GetXaxis()->SetLimits(genPositionZ - 500, finalPositionZ + 500);
+    fieldGr->GetHistogram()->SetMaximum(2.5);
+    fieldGr->GetHistogram()->SetMinimum(0);
+    fieldGr->GetXaxis()->SetTitle("Z [mm]");
+    fieldGr->GetXaxis()->SetTitleSize(0.05);
+    fieldGr->GetXaxis()->SetLabelSize(0.05);
+    fieldGr->GetYaxis()->SetTitle("B [T]");
+    fieldGr->GetYaxis()->SetTitleOffset(1.3);
+    fieldGr->GetYaxis()->SetTitleSize(0.05);
+    fieldGr->GetYaxis()->SetLabelSize(0.05);
+    fieldGr->Draw("AL");
+
+    return fCanvas;
+}
+
 ///////////////////////////////////////////////
 /// \brief A method that creates a canvas where tracks traversing the magnetic volume
 /// are drawm together with their corresponding field intensity profile along the Z-axis.
@@ -966,7 +1025,7 @@ TVector3 TRestAxionMagneticField::GetMagneticField(TVector3 pos, Bool_t showWarn
         return TVector3(0, 0, 0);
     } else {
         if (IsFieldConstant(id)) return fConstantField[id];
-        TVector3 node = GetMagneticVolumeNode(fMagneticFieldVolumes[id], pos);
+        TVector3 node = GetMagneticVolumeNode((size_t)id, pos);
         Int_t nX = node.X();
         Int_t nY = node.Y();
         Int_t nZ = node.Z();
@@ -1214,6 +1273,49 @@ std::vector<Double_t> TRestAxionMagneticField::GetTransversalComponentAlongPath(
     return Bt;
 }
 
+///////////////////////////////////////////////
+/// \brief It returns a vector describing the magnetic field profile component requested using the `axis`
+/// argument which may take values 0,1,2 for X,Y,Z components. The field profile will be constructed in the
+/// track defined between `from` and `to` positions given by argument.
+///
+/// The differential element `dl` is by default 1mm, but it can be modified through the third argument of
+/// this function.
+///
+/// The maximum number of divisions (unlimited by default) of the output vector can be fixed by the forth
+/// argument. In that case, the differential element `dl` length might be increased to fullfil such
+/// condition.
+///
+std::vector<Double_t> TRestAxionMagneticField::GetComponentAlongPath(Int_t axis, TVector3 from, TVector3 to,
+                                                                     Double_t dl, Int_t Nmax) {
+    std::vector<Double_t> Bt;
+    if (axis != 0 && axis != 1 && axis != 2) {
+        RESTError << "TRestAxionMagneticField::GetComponentAlongPath. Axis must take values 0,1 or 2"
+                  << RESTendl;
+        return Bt;
+    }
+    Double_t length = (to - from).Mag();
+
+    Double_t diff = dl;
+    if (Nmax > 0) {
+        if (length / dl > Nmax) {
+            diff = length / Nmax;
+            RESTWarning << "TRestAxionMagneticField::GetComponentAlongPath. Nmax reached!" << RESTendl;
+            RESTWarning << "Nmax = " << Nmax << RESTendl;
+            RESTWarning << "Adjusting differential step to : " << diff << " mm" << RESTendl;
+        }
+    }
+
+    TVector3 direction = (to - from).Unit();
+
+    for (Double_t d = 0; d < length; d += diff) {
+        if (axis == 0) Bt.push_back(GetMagneticField(from + d * direction).X());
+        if (axis == 1) Bt.push_back(GetMagneticField(from + d * direction).Y());
+        if (axis == 2) Bt.push_back(GetMagneticField(from + d * direction).Z());
+    }
+
+    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.
@@ -1225,6 +1327,7 @@ void TRestAxionMagneticField::SetTrack(const TVector3& position, const TVector3&
         fTrackStart = TVector3(0, 0, 0);
         fTrackDirection = TVector3(0, 0, 0);
         fTrackLength = 0;
+        return;
     }
 
     fTrackStart = trackBounds[0];
@@ -1327,10 +1430,10 @@ TVector3 TRestAxionMagneticField::GetFieldAverageTransverseVector(TVector3 from,
 ///
 /// This method will be made private, no reason to use it outside this class.
 ///
-TVector3 TRestAxionMagneticField::GetMagneticVolumeNode(MagneticFieldVolume mVol, TVector3 pos) {
-    Int_t nx = mVol.mesh.GetNodeX(pos.X());
-    Int_t ny = mVol.mesh.GetNodeY(pos.Y());
-    Int_t nz = mVol.mesh.GetNodeZ(pos.Z());
+TVector3 TRestAxionMagneticField::GetMagneticVolumeNode(size_t id, TVector3 pos) {
+    Int_t nx = fMagneticFieldVolumes[id].mesh.GetNodeX(pos.X());
+    Int_t ny = fMagneticFieldVolumes[id].mesh.GetNodeY(pos.Y());
+    Int_t nz = fMagneticFieldVolumes[id].mesh.GetNodeZ(pos.Z());
     return TVector3(nx, ny, nz);
 }