From 29f1a91ff967350cab77f93e12ad91fc3c0f1ff9 Mon Sep 17 00:00:00 2001
From: lintao <lintao51@gmail.com>
Date: Thu, 27 May 2021 23:30:34 +0800
Subject: [PATCH] WIP: calculate the Br/Bz.
---
.../src/FieldMapFileProvider.cpp | 17 +++-
.../src/FieldMapFileProvider.h | 4 +-
.../src/GenericBFieldMapBrBz.cpp | 94 ++++++++++++++++++-
.../MagneticFieldMap/src/IFieldMapProvider.h | 5 +-
4 files changed, 109 insertions(+), 11 deletions(-)
diff --git a/Detector/MagneticFieldMap/src/FieldMapFileProvider.cpp b/Detector/MagneticFieldMap/src/FieldMapFileProvider.cpp
index 2865f6bf..2eb352e4 100644
--- a/Detector/MagneticFieldMap/src/FieldMapFileProvider.cpp
+++ b/Detector/MagneticFieldMap/src/FieldMapFileProvider.cpp
@@ -16,7 +16,7 @@ FieldMapFileProvider::FieldMapFileProvider(const std::string& url_)
init();
}
-int FieldMapFileProvider::rBinIdx(double r) {
+int FieldMapFileProvider::rBinIdx(double r, double& rn) {
//
// | --- | --- | --- |
// ^ ^ ^
@@ -30,20 +30,30 @@ int FieldMapFileProvider::rBinIdx(double r) {
int idx = -1;
+ // std::cout << "FieldMapFileProvider::rBinIdx: "
+ // << " r: " << r
+ // << " rBinMin: " << rBinMin
+ // << " rBinMax: " << rBinMax
+ // << std::endl;
+
if ( rBinMin <= absr && absr < rBinMax) {
idx = (absr - rBinMin) / drBin;
+ double r0 = rBinMin + idx*drBin;
+ rn = (absr - r0)/drBin;
}
return idx;
}
-int FieldMapFileProvider::zBinIdx(double z) {
+int FieldMapFileProvider::zBinIdx(double z, double& zn) {
double absz = std::fabs(z);
int idx = -1;
if ( zBinMin <= absz && absz < zBinMax) {
- idx = (absz - zBinMin) / drBin;
+ idx = (absz - zBinMin) / dzBin;
+ double z0 = zBinMin + idx*dzBin;
+ zn = (absz - z0)/dzBin;
}
return idx;
@@ -55,6 +65,7 @@ void FieldMapFileProvider::access(int rbin, int zbin, double& Br, double& Bz) {
if ((rbin < 0 || rbin >= nr) || (zbin < 0 || zbin >= nr)) {
Br = 0;
Bz = 0;
+ return;
}
// convert to the internal table (with left col and top row)
diff --git a/Detector/MagneticFieldMap/src/FieldMapFileProvider.h b/Detector/MagneticFieldMap/src/FieldMapFileProvider.h
index a6712c66..8a754a7a 100644
--- a/Detector/MagneticFieldMap/src/FieldMapFileProvider.h
+++ b/Detector/MagneticFieldMap/src/FieldMapFileProvider.h
@@ -12,8 +12,8 @@ public:
FieldMapFileProvider(const std::string& url_);
// Meta data about the map
- virtual int rBinIdx(double r);
- virtual int zBinIdx(double z);
+ virtual int rBinIdx(double r, double& rn);
+ virtual int zBinIdx(double z, double& zn);
// The Br and Bz
virtual void access(int rbin, int zbin, double& Br, double& Bz);
diff --git a/Detector/MagneticFieldMap/src/GenericBFieldMapBrBz.cpp b/Detector/MagneticFieldMap/src/GenericBFieldMapBrBz.cpp
index 568d70c8..9b3218b3 100644
--- a/Detector/MagneticFieldMap/src/GenericBFieldMapBrBz.cpp
+++ b/Detector/MagneticFieldMap/src/GenericBFieldMapBrBz.cpp
@@ -3,6 +3,10 @@
#include "FieldMapFileProvider.h"
+#include <cmath>
+
+#include "DD4hep/DD4hepUnits.h"
+
GenericBFieldMapBrBz::GenericBFieldMapBrBz()
: m_provider(nullptr) {
type = dd4hep::CartesianField::MAGNETIC;
@@ -10,11 +14,93 @@ GenericBFieldMapBrBz::GenericBFieldMapBrBz()
}
void GenericBFieldMapBrBz::fieldComponents(const double* pos, double* field) {
- double curfield[3] = {0.0, 0.0, 0.0};
+ if (!m_provider) {
+ std::string error_msg = "[ERROR] GenericBFieldMapBrBz: No provider! ";
+ throw std::runtime_error(error_msg);
+ }
+
+ // convert pos to r/z
+ double x = pos[0] / dd4hep::m; // convert to meter
+ double y = pos[1] / dd4hep::m;
+ double z = pos[2] / dd4hep::m;
+ double r = sqrt(x*x+y*y);
+ double phi = atan2(y, x);
+
+ // std::cout << " r: " << r
+ // << " x: " << x
+ // << " y: " << y
+ // << " z: " << z
+ // << " field Bx: " << field[0]
+ // << " field By: " << field[1]
+ // << " field Bz: " << field[2]
+ // << std::endl;
+
+ // get the point r0/z0, r1/z1
+
+ // -- normalized r/z at the bin
+ double rn = 0.0;
+ double zn = 0.0;
+
+ int ir0 = 0;
+ int iz0 = 0;
+
+ ir0 = m_provider->rBinIdx(r, rn);
+ iz0 = m_provider->zBinIdx(z, zn);
+
+ // not in the valid return
+ if (ir0 < 0 || iz0 < 0) {
+ // std::cout << "SKIP due to "
+ // << " r: " << r
+ // << " x: " << x
+ // << " y: " << y
+ // << " z: " << z
+ // << " ir0: " << ir0
+ // << " iz0: " << iz0
+ // << " rn: " << rn
+ // << " zn: " << zn
+ // << std::endl;
+ return;
+ }
+
+ // std::cout << " r: " << r
+ // << " x: " << x
+ // << " y: " << y
+ // << " z: " << z
+ // << " ir0: " << ir0
+ // << " iz0: " << iz0
+ // << " rn: " << rn
+ // << " zn: " << zn
+ // << std::endl;
+
+ int ir1 = ir0 + 1;
+ int iz1 = iz0 + 1;
+
+ // get the Br/Bz at four points
+ double Br_r0z0 = 0.0; double Bz_r0z0 = 0.0;
+ double Br_r1z0 = 0.0; double Bz_r1z0 = 0.0;
+ double Br_r0z1 = 0.0; double Bz_r0z1 = 0.0;
+ double Br_r1z1 = 0.0; double Bz_r1z1 = 0.0;
+
+ // calculate the field at r/z
+ m_provider->access(ir0, iz0, Br_r0z0, Bz_r0z0);
+ m_provider->access(ir1, iz0, Br_r1z0, Bz_r1z0);
+ m_provider->access(ir0, iz1, Br_r0z1, Bz_r0z1);
+ m_provider->access(ir1, iz1, Br_r1z1, Bz_r1z1);
+
+ double Br = (1.0 - rn) * (1.0 - zn) * Br_r0z0
+ + rn * (1.0 - zn) * Br_r1z0
+ + (1.0 - rn) * zn * Br_r0z1
+ + rn * zn * Br_r1z1;
+
+ double Bz = (1.0 - rn) * (1.0 - zn) * Bz_r0z0
+ + rn * (1.0 - zn) * Bz_r1z0
+ + (1.0 - rn) * zn * Bz_r0z1
+ + rn * zn * Bz_r1z1;
- field[0] += curfield[0];
- field[1] += curfield[1];
- field[1] += curfield[2];
+ // update the global field
+ field[0] += Br*cos(phi);
+ field[1] += Br*sin(phi);
+ field[2] += Bz;
return;
}
diff --git a/Detector/MagneticFieldMap/src/IFieldMapProvider.h b/Detector/MagneticFieldMap/src/IFieldMapProvider.h
index 3d238422..b1e8d948 100644
--- a/Detector/MagneticFieldMap/src/IFieldMapProvider.h
+++ b/Detector/MagneticFieldMap/src/IFieldMapProvider.h
@@ -13,8 +13,9 @@
class IFieldMapProvider {
public:
// Meta data about the map
- virtual int rBinIdx(double r) = 0;
- virtual int zBinIdx(double z) = 0;
+ // return rn/zn is the normalized value in the bin. the value is [0,1]
+ virtual int rBinIdx(double r, double& rn) = 0;
+ virtual int zBinIdx(double z, double& zn) = 0;
// The Br and Bz
virtual void access(int rbin, int zbin, double& Br, double& Bz) = 0;
--
GitLab