diff --git a/Detector/CMakeLists.txt b/Detector/CMakeLists.txt
index 57e6f2ce1389dc6690edf864c20b8862b81313f5..d628d7ab5b402422bff9135b662a12688c56841f 100644
--- a/Detector/CMakeLists.txt
+++ b/Detector/CMakeLists.txt
@@ -7,3 +7,4 @@ add_subdirectory(DetInterface)
add_subdirectory(DetSegmentation)
add_subdirectory(GeomSvc)
add_subdirectory(Identifier)
+add_subdirectory(MagneticFieldMap)
diff --git a/Detector/MagneticFieldMap/CMakeLists.txt b/Detector/MagneticFieldMap/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7c65b50318a2d0040922057a2fcc36a84db3ec50
--- /dev/null
+++ b/Detector/MagneticFieldMap/CMakeLists.txt
@@ -0,0 +1,22 @@
+#################################################################################
+##Package : MagneticFieldMap
+#################################################################################
+
+gaudi_add_module(MagneticFieldMap
+ SOURCES src/GenericBFieldMapBrBz.cpp
+ src/GenericBFieldMapBrBzFactory.cpp
+ src/FieldMapFileProvider.cpp
+ LINK Gaudi::GaudiKernel
+ ${DD4hep_COMPONENT_LIBRARIES}
+ ${ROOT_LIBRARIES}
+)
+
+set(LIBRARY_OUTPUT_PATH ${CMAKE_LIBRARY_OUTPUT_DIRECTORY})
+message(STATUS "LIBRARY_OUTPUT_PATH -> ${LIBRARY_OUTPUT_PATH}")
+dd4hep_generate_rootmap(MagneticFieldMap)
+
+install(TARGETS MagneticFieldMap
+ EXPORT CEPCSWTargets
+ RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+ LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+ COMPONENT dev)
diff --git a/Detector/MagneticFieldMap/src/FieldMapFileProvider.cpp b/Detector/MagneticFieldMap/src/FieldMapFileProvider.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2eb352e4912a7d63eb797104beab25f59ec706ec
--- /dev/null
+++ b/Detector/MagneticFieldMap/src/FieldMapFileProvider.cpp
@@ -0,0 +1,283 @@
+#include "FieldMapFileProvider.h"
+
+#include <cmath>
+#include <string>
+#include <map>
+#include <vector>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+
+FieldMapFileProvider::FieldMapFileProvider(const std::string& url_)
+ : m_url(url_), nr(-1), nz(-1),
+ rBinMin(-1), zBinMin(-1),
+ rBinMax(-1), zBinMax(-1),
+ drBin(-1), dzBin(-1) {
+ init();
+}
+
+int FieldMapFileProvider::rBinIdx(double r, double& rn) {
+ //
+ // | --- | --- | --- |
+ // ^ ^ ^
+ // rmin | rmax
+ // r
+ // |
+ // return
+
+
+ double absr = std::fabs(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, double& zn) {
+ double absz = std::fabs(z);
+
+ int idx = -1;
+
+ if ( zBinMin <= absz && absz < zBinMax) {
+ idx = (absz - zBinMin) / dzBin;
+ double z0 = zBinMin + idx*dzBin;
+ zn = (absz - z0)/dzBin;
+ }
+
+ return idx;
+}
+
+void FieldMapFileProvider::access(int rbin, int zbin, double& Br, double& Bz) {
+ // the valid bin should between [0, n)
+ // if the point is not in the valid region, return 0
+ 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)
+ int ridx = rbin+1;
+ int zidx = zbin+1;
+
+ int globalidx = ridx + zidx*(nr+1);
+
+ Br = Brdata[globalidx];
+ Bz = Bzdata[globalidx];
+}
+
+// ======================
+// Below are private impl
+// ======================
+
+void FieldMapFileProvider::init() {
+ // parse the url
+ // Br=/tmp/lint/CEPCSW/Br.csv;Bz=/tmp/lint/CEPCSW/Bz.csv
+
+ std::map<std::string, std::string> url_parsed;
+
+ size_t idx_begin = 0;
+ size_t idx_end = m_url.find(";");
+
+ while (true) {
+ std::string keyval = m_url.substr(idx_begin, idx_end-idx_begin);
+
+ if (keyval.size()) {
+ std::cout << "---> keyval: " << keyval << std::endl;
+ // separate by '='
+ size_t idx_sep = keyval.find("=");
+ if (idx_sep == std::string::npos) {
+ std::string error_msg = "[ERROR] FieldMapFileProvider: Pleaes specify label=path. ";
+ throw std::runtime_error(error_msg);
+
+ }
+ std::string key = keyval.substr(0, idx_sep);
+ std::string val = keyval.substr(idx_sep+1);
+ std::cout << "-----> key: " << key << std::endl;
+ std::cout << "-----> val: " << val << std::endl;
+ // insert into map
+
+ // duplicated
+ if (url_parsed.count(key)) {
+ std::string error_msg = "[ERROR] FieldMapFileProvider: duplicated key '" + key + "'";
+ throw std::runtime_error(error_msg);
+ }
+
+ url_parsed[key] = val;
+ }
+
+ if (idx_end == std::string::npos) {
+ break;
+ }
+
+ idx_begin = idx_end+1;
+ idx_end = m_url.find(';', idx_begin);
+ }
+
+ // access the map 'url_parsed'
+ std::string key;
+
+ // load Br
+ key = "Br";
+ if (url_parsed.count(key) == 0) {
+ std::string error_msg = "[ERROR] FieldMapFileProvider: missing key '" + key + "'";
+ throw std::runtime_error(error_msg);
+ }
+ loadBr(url_parsed[key]);
+
+ // load Bz
+ key = "Bz";
+ if (url_parsed.count(key) == 0) {
+ std::string error_msg = "[ERROR] FieldMapFileProvider: missing key '" + key + "'";
+ throw std::runtime_error(error_msg);
+ }
+ loadBz(url_parsed[key]);
+
+}
+
+
+
+bool FieldMapFileProvider::loadBr(const std::string& fn) {
+
+ // ----> r
+ // |
+ // v
+ // z
+ // the shape
+ int ncol = -1;
+ int nrow = -1;
+
+ // load data
+ bool st = loadCSV(fn, Brdata, ncol, nrow);
+
+ // update the metadata
+ // for (int i = 0; i < ncol; ++i) {
+ // std::cout << "i: " << i << " data[i]: " << data[i] << std::endl;
+ // }
+ nr = ncol-1; // skip the top row and left col
+ nz = nrow-1;
+
+ rBinMin = Brdata[1];
+ zBinMin = Brdata[ncol];
+
+ rBinMax = Brdata[ncol-1];
+ zBinMax = Brdata[ncol*(nrow-1)];
+
+ drBin = (rBinMax-rBinMin) / (nr-1);
+ dzBin = (zBinMax-zBinMin) / (nz-1);
+
+ std::cout << "nr: " << nr << std::endl;
+ std::cout << "nz: " << nz << std::endl;
+ std::cout << "rBinMin: " << rBinMin << std::endl;
+ std::cout << "zBinMin: " << zBinMin << std::endl;
+ std::cout << "rBinMax: " << rBinMax << std::endl;
+ std::cout << "zBinMax: " << zBinMax << std::endl;
+ std::cout << "drBin: " << drBin << std::endl;
+ std::cout << "dzBin: " << dzBin << std::endl;
+
+
+
+ return true;
+}
+
+bool FieldMapFileProvider::loadBz(const std::string& fn) {
+
+ // the shape
+ int ncol = -1;
+ int nrow = -1;
+ bool st = loadCSV(fn, Bzdata, ncol, nrow);
+
+ return true;
+}
+
+bool FieldMapFileProvider::loadCSV(const std::string& fn,
+ std::vector<double>& data,
+ int& ncol, int& nrow) {
+
+ std::ifstream input(fn);
+ std::string tmpline;
+
+ ncol = 0;
+ nrow = 0;
+
+
+ while(std::getline(input, tmpline)) {
+ std::cout << "------> " << tmpline << std::endl;
+ ++nrow;
+
+ // split by ,
+ int col_per_line = 0;
+
+ size_t idx_begin = 0;
+ size_t idx_end = tmpline.find(",");
+
+ while (true) {
+ std::cout << "----------------> idx: "
+ << idx_begin << "->" << idx_end << std::endl;
+ std::string val = tmpline.substr(idx_begin, idx_end-idx_begin);
+ double v = 0.0;
+ // parse val
+ std::stringstream ss;
+ ss << val;
+ ss >> v;
+ // if (not ss.good()) {
+ // v = 0;
+ // ss.clear();
+ // }
+
+ std::cout << "----------------> v: " << v << std::endl;
+ // if (val.size()==0 || ) { // if it is empty or not number, use 0
+ // data.push_back(0.0);
+ // } else {
+ // }
+ data.push_back(v);
+
+ ++col_per_line;
+
+ // break
+ if (idx_end == std::string::npos) {
+ break;
+ }
+
+ // goto next val
+ idx_begin = idx_end+1;
+ idx_end = tmpline.find(',', idx_begin);
+
+ }
+
+ // if it is the first line, we need to store the column sizes
+ if (ncol == 0) {
+ ncol = col_per_line;
+ }
+
+ std::cout << "--> ncol: " << ncol
+ << " col_per_line: " << col_per_line
+ << std::endl;
+
+ if (ncol != col_per_line) {
+ std::string error_msg = "[ERROR] FieldMapFileProvider: Mismatch columns. ";
+ throw std::runtime_error(error_msg);
+ }
+
+ }
+
+ std::cout << "SUMMARY: "
+ << " size of the data: " << data.size()
+ << " shape (ncol: " << ncol
+ << " , nrow: " << nrow
+ << ")" << std::endl;
+
+ return true;
+}
diff --git a/Detector/MagneticFieldMap/src/FieldMapFileProvider.h b/Detector/MagneticFieldMap/src/FieldMapFileProvider.h
new file mode 100644
index 0000000000000000000000000000000000000000..8a754a7a2a1ae41864ca7c339adc968a12fb288a
--- /dev/null
+++ b/Detector/MagneticFieldMap/src/FieldMapFileProvider.h
@@ -0,0 +1,52 @@
+#ifndef FieldMapFileProvider_h
+#define FieldMapFileProvider_h
+
+#include "IFieldMapProvider.h"
+
+#include <string>
+#include <vector>
+
+class FieldMapFileProvider: public IFieldMapProvider {
+
+public:
+ FieldMapFileProvider(const std::string& url_);
+
+ // Meta data about the map
+ 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);
+
+private:
+ void init();
+
+ bool loadBr(const std::string& fn);
+ bool loadBz(const std::string& fn);
+
+ bool loadCSV(const std::string& fn, std::vector<double>& data, int& ncol, int& nrow);
+private:
+ std::string m_url; // the url passed from GenericBFieldMapBrBz
+
+ std::string m_filename_Br;
+ std::string m_filename_Bz;
+
+ // the data will include the left col (z) and top row (r)
+ std::vector<double> Brdata;
+ std::vector<double> Bzdata;
+
+ // in this impl, the rBin/zBin should be consistent between Br and Bz.
+ int nr; // include the two endpoints, so nr = nrBin + 1
+ int nz;
+
+ double rBinMin;
+ double zBinMin;
+
+ double rBinMax;
+ double zBinMax;
+
+ double drBin;
+ double dzBin;
+};
+
+#endif
diff --git a/Detector/MagneticFieldMap/src/GenericBFieldMapBrBz.cpp b/Detector/MagneticFieldMap/src/GenericBFieldMapBrBz.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9b3218b3c2f47521c892e854e2732728cb7e7a72
--- /dev/null
+++ b/Detector/MagneticFieldMap/src/GenericBFieldMapBrBz.cpp
@@ -0,0 +1,118 @@
+
+#include "GenericBFieldMapBrBz.h"
+
+#include "FieldMapFileProvider.h"
+
+#include <cmath>
+
+#include "DD4hep/DD4hepUnits.h"
+
+GenericBFieldMapBrBz::GenericBFieldMapBrBz()
+ : m_provider(nullptr) {
+ type = dd4hep::CartesianField::MAGNETIC;
+
+}
+
+void GenericBFieldMapBrBz::fieldComponents(const double* pos, double* field) {
+ 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;
+
+ // update the global field
+ field[0] += Br*cos(phi);
+ field[1] += Br*sin(phi);
+ field[2] += Bz;
+
+ return;
+}
+
+void GenericBFieldMapBrBz::init_provider(const std::string& provider, const std::string& url) {
+ if (provider == "file") {
+ std::cout << "Initialize provider with file. " << std::endl;
+ m_provider = new FieldMapFileProvider(url);
+ } else if (provider == "db") {
+ std::cout << "Initialize provider with file. " << std::endl;
+ } else {
+ std::string error_msg = "[ERROR] GenericBFieldMapBrBz: Unknown provider: " + provider;
+ throw std::runtime_error(error_msg);
+ }
+}
diff --git a/Detector/MagneticFieldMap/src/GenericBFieldMapBrBz.h b/Detector/MagneticFieldMap/src/GenericBFieldMapBrBz.h
new file mode 100644
index 0000000000000000000000000000000000000000..93abfcfbbb5650d5e5945b774280e6586b1c17c3
--- /dev/null
+++ b/Detector/MagneticFieldMap/src/GenericBFieldMapBrBz.h
@@ -0,0 +1,35 @@
+#ifndef GenericBFieldMapBrBz_h
+#define GenericBFieldMapBrBz_h
+
+/*
+ * GenericBFieldMapBrBz is an extension of Cartesian Field in DD4hep.
+ * - It enables the DD4hep to access the Magnetic Service from Gaudi.
+ * - It also enables the calculation of Br/Bz at position X.
+ * - It will get the map from an abstract class IFieldMapProvider.
+ *
+ * -- Tao Lin <lintao AT ihep.ac.cn>
+ */
+
+#include <DD4hep/FieldTypes.h>
+
+#include "IFieldMapProvider.h"
+
+
+class GenericBFieldMapBrBz: public dd4hep::CartesianField::Object {
+public:
+
+ GenericBFieldMapBrBz();
+
+ virtual void fieldComponents(const double* pos, double* field);
+
+public:
+ // following are interfaces to configure this field map
+ void init_provider(const std::string& provider, const std::string& url);
+
+private:
+
+ IFieldMapProvider* m_provider;
+};
+
+#endif
+
diff --git a/Detector/MagneticFieldMap/src/GenericBFieldMapBrBzFactory.cpp b/Detector/MagneticFieldMap/src/GenericBFieldMapBrBzFactory.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d7d9fe302d75e214e290f52ebf2382441df305c0
--- /dev/null
+++ b/Detector/MagneticFieldMap/src/GenericBFieldMapBrBzFactory.cpp
@@ -0,0 +1,79 @@
+/*
+ * In this file, the xml is parsed and the GenericBFieldMapBrBz object is created and configured.
+ *
+ * The properties for the GenericBFieldMapBrBz
+ * - provider (attribute)
+ * - [file, db]
+ * - source (tag)
+ * - the attributes include:
+ * - url
+ * - file path for the 'file' mode.
+ * - DB instance ... for the 'db' mode.
+ * - rhoMin, rhoMax, zMin, zMax
+ *
+ * -- Tao Lin <lintao AT ihep.ac.cn>
+ */
+#include "GenericBFieldMapBrBz.h"
+
+
+#include <DD4hep/Version.h>
+#if DD4HEP_VERSION_GE(0,24)
+#include <DD4hep/detail/Handle.inl>
+#else
+#include <DD4hep/Handle.inl>
+#endif
+
+#include <DD4hep/FieldTypes.h>
+
+
+#include <DD4hep/DetFactoryHelper.h>
+#include <XML/Utilities.h>
+
+
+static dd4hep::Ref_t create_GenericBFieldMapBrBz(dd4hep::Detector& ,
+ dd4hep::xml::Handle_t handle ) {
+
+ // 1. retrieve the parameters from xml
+ dd4hep::xml::Component xmlParameter(handle);
+
+ // - provider
+ bool hasProvider = xmlParameter.hasAttr(_Unicode(provider));
+ if (!hasProvider) {
+ std::string error_msg = "[ERROR] GenericBFieldMapBrBz: Must specify the 'provider'. ";
+ throw std::runtime_error(error_msg);
+ }
+
+ std::string provider = xmlParameter.attr<std::string>(_Unicode(provider));
+
+ // - source
+ bool hasSource = xmlParameter.hasChild(_Unicode(source));
+ if (!hasSource) {
+ std::string error_msg = "[ERROR] GenericBFieldMapBrBz: Must specify the 'source' tag. ";
+ throw std::runtime_error(error_msg);
+
+ }
+
+ dd4hep::xml::Component source(xmlParameter.child(_Unicode(source)));
+
+ // get the url
+ bool hasUrl = source.hasAttr(_Unicode(url));
+ if (!hasUrl) {
+ std::string error_msg = "[ERROR] GenericBFieldMapBrBz: Must specify the 'url' in 'source' tag. ";
+ throw std::runtime_error(error_msg);
+ }
+
+ std::string url = source.attr<std::string>(_Unicode(url));
+
+ // 2. create the CartesianField
+ dd4hep::CartesianField obj;
+ GenericBFieldMapBrBz* ptr = new GenericBFieldMapBrBz();
+
+ ptr->init_provider(provider, url);
+
+ obj.assign(ptr, xmlParameter.nameStr(), xmlParameter.typeStr());
+
+ return obj;
+
+}
+
+DECLARE_XMLELEMENT(GenericBFieldMapBrBz,create_GenericBFieldMapBrBz)
diff --git a/Detector/MagneticFieldMap/src/IFieldMapProvider.h b/Detector/MagneticFieldMap/src/IFieldMapProvider.h
new file mode 100644
index 0000000000000000000000000000000000000000..b1e8d948ef3797926f43a892fe21ef61eafdb207
--- /dev/null
+++ b/Detector/MagneticFieldMap/src/IFieldMapProvider.h
@@ -0,0 +1,26 @@
+#ifndef IFieldMapProvider_h
+#define IFieldMapProvider_h
+
+/*
+ * IFieldMapProvider will provide the B-Field at a 2D grid (r-z).
+ * The provider will return the bin index for r/z. Then the GenericBFieldMapBrBz
+ * will use this bin index to get the values at (r0,z0), (r1,z0), (r0,z1), (r1,z1).
+ * The interpolation is not computated IFieldMapProvider. Please see GenericBFieldMapBrBz.
+ *
+ * -- Tao Lin <lintao AT ihep.ac.cn>
+ */
+
+class IFieldMapProvider {
+public:
+ // Meta data about the map
+ // 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;
+
+};
+
+#endif
+