File indexing completed on 2025-08-06 08:17:58
0001 #include "PHField3DCartesian.h"
0002
0003 #include <phool/phool.h>
0004
0005 #include <TDirectory.h> // for TDirectory, gDirectory
0006 #include <TFile.h>
0007 #include <TNtuple.h>
0008 #include <TSystem.h>
0009
0010 #include <Geant4/G4SystemOfUnits.hh>
0011
0012
0013 #pragma GCC diagnostic push
0014 #pragma GCC diagnostic ignored "-Wshadow"
0015 #include <boost/stacktrace.hpp>
0016 #pragma GCC diagnostic pop
0017
0018 #include <cassert>
0019 #include <cmath>
0020 #include <cstdlib>
0021 #include <iostream>
0022 #include <iterator>
0023 #include <set>
0024 #include <utility>
0025
0026 PHField3DCartesian::PHField3DCartesian(const std::string &fname, const float magfield_rescale, const float innerradius, const float outerradius, const float size_z)
0027 : filename(fname)
0028 {
0029 std::cout << "PHField3DCartesian::PHField3DCartesian" << std::endl;
0030
0031 for (int i = 0; i < 2; i++)
0032 {
0033 for (int j = 0; j < 2; j++)
0034 {
0035 for (int k = 0; k < 2; k++)
0036 {
0037 for (int l = 0; l < 3; l++)
0038 {
0039 bf[i][j][k][l] = NAN;
0040 }
0041 }
0042 }
0043 }
0044 std::cout << "\n================ Begin Construct Mag Field =====================" << std::endl;
0045 std::cout << "\n-----------------------------------------------------------"
0046 << "\n Magnetic field Module - Verbosity:"
0047 << "\n-----------------------------------------------------------";
0048
0049
0050 TFile *rootinput = TFile::Open(filename.c_str());
0051 if (!rootinput)
0052 {
0053 std::cout << "\n could not open " << filename << " exiting now" << std::endl;
0054 gSystem->Exit(1);
0055 exit(1);
0056 }
0057 std::cout << "\n ---> "
0058 "Reading the field grid from "
0059 << filename << " ... " << std::endl;
0060
0061
0062 TNtuple *field_map = nullptr;
0063 rootinput->GetObject("fieldmap", field_map);
0064 if (field_map == nullptr)
0065 {
0066 std::cout << PHWHERE << " Could not load fieldmap ntuple from "
0067 << filename << " exiting now" << std::endl;
0068 gSystem->Exit(1);
0069 exit(1);
0070 }
0071 Float_t ROOT_X, ROOT_Y, ROOT_Z;
0072 Float_t ROOT_BX, ROOT_BY, ROOT_BZ;
0073 field_map->SetBranchAddress("x", &ROOT_X);
0074 field_map->SetBranchAddress("y", &ROOT_Y);
0075 field_map->SetBranchAddress("z", &ROOT_Z);
0076 field_map->SetBranchAddress("bx", &ROOT_BX);
0077 field_map->SetBranchAddress("by", &ROOT_BY);
0078 field_map->SetBranchAddress("bz", &ROOT_BZ);
0079 for (int i = 0; i < field_map->GetEntries(); i++)
0080 {
0081 field_map->GetEntry(i);
0082 trio coord_key(ROOT_X * cm, ROOT_Y * cm, ROOT_Z * cm);
0083 trio field_val(ROOT_BX * tesla * magfield_rescale, ROOT_BY * tesla * magfield_rescale, ROOT_BZ * tesla * magfield_rescale);
0084 xvals.insert(ROOT_X * cm);
0085 yvals.insert(ROOT_Y * cm);
0086 zvals.insert(ROOT_Z * cm);
0087 if ((std::sqrt(ROOT_X * cm * ROOT_X * cm + ROOT_Y * cm * ROOT_Y * cm) >= innerradius &&
0088 std::sqrt(ROOT_X * cm * ROOT_X * cm + ROOT_Y * cm * ROOT_Y * cm) <= outerradius) ||
0089 std::abs(ROOT_Z * cm) > size_z)
0090 {
0091 fieldmap[coord_key] = field_val;
0092 }
0093 }
0094 xmin = *(xvals.begin());
0095 xmax = *(xvals.rbegin());
0096
0097 ymin = *(yvals.begin());
0098 ymax = *(yvals.rbegin());
0099 if (ymin != xmin || ymax != xmax)
0100 {
0101 std::cout << "PHField3DCartesian: Compiler bug!!!!!!!! Do not use inlining!!!!!!" << std::endl;
0102 std::cout << "exiting now - recompile with -fno-inline" << std::endl;
0103 exit(1);
0104 }
0105
0106 zmin = *(zvals.begin());
0107 zmax = *(zvals.rbegin());
0108
0109 xstepsize = (xmax - xmin) / (xvals.size() - 1);
0110 ystepsize = (ymax - ymin) / (yvals.size() - 1);
0111 zstepsize = (zmax - zmin) / (zvals.size() - 1);
0112
0113 delete field_map;
0114 delete rootinput;
0115 std::cout << "\n================= End Construct Mag Field ======================\n"
0116 << std::endl;
0117 }
0118
0119 PHField3DCartesian::~PHField3DCartesian()
0120 {
0121 std::cout << "PHField3DCartesian::~PHField3DCartesian" << std::endl;
0122 if (Verbosity() > 0)
0123 {
0124 std::cout << "PHField3DCartesian: cache hits: " << cache_hits
0125 << " cache misses: " << cache_misses
0126 << std::endl;
0127 }
0128 }
0129
0130 void PHField3DCartesian::GetFieldValue(const double point[4], double *Bfield) const
0131 {
0132
0133 static double xsav = -1000000.;
0134 static double ysav = -1000000.;
0135 static double zsav = -1000000.;
0136
0137 const double& x = point[0];
0138 const double& y = point[1];
0139 const double& z = point[2];
0140
0141 Bfield[0] = 0.0;
0142 Bfield[1] = 0.0;
0143 Bfield[2] = 0.0;
0144 if (!std::isfinite(x) || !std::isfinite(y) || !std::isfinite(z))
0145 {
0146 static int ifirst = 0;
0147 if (ifirst < 10)
0148 {
0149 std::cout << "PHField3DCartesian::GetFieldValue: "
0150 << "Invalid coordinates: "
0151 << "x: " << x / cm
0152 << ", y: " << y / cm
0153 << ", z: " << z / cm
0154 << " bailing out returning zero bfield"
0155 << std::endl;
0156 std::cout << "previous point: "
0157 << "x: " << xsav / cm
0158 << ", y: " << ysav / cm
0159 << ", z: " << zsav / cm
0160 << std::endl;
0161 std::cout << "Here is the stacktrace: " << std::endl;
0162 std::cout << boost::stacktrace::stacktrace();
0163 std::cout << "This is not a segfault. Check the stacktrace for the guilty party (typically #2)" << std::endl;
0164 ifirst++;
0165 }
0166 return;
0167 }
0168
0169 xsav = x;
0170 ysav = y;
0171 zsav = z;
0172
0173 if (point[0] < xmin || point[0] > xmax ||
0174 point[1] < ymin || point[1] > ymax ||
0175 point[2] < zmin || point[2] > zmax)
0176 {
0177 return;
0178 }
0179
0180 double xkey[2];
0181 std::set<float>::const_iterator it = xvals.lower_bound(x);
0182 xkey[0] = *it;
0183 if (it == xvals.begin())
0184 {
0185 xkey[1] = *it;
0186 if (x < xkey[0])
0187 {
0188 std::cout << PHWHERE << ": This should not happen! x too small - outside range: " << x / cm << std::endl;
0189 return;
0190 }
0191 }
0192 else
0193 {
0194 --it;
0195 xkey[1] = *it;
0196 }
0197
0198 double ykey[2];
0199 it = yvals.lower_bound(y);
0200 ykey[0] = *it;
0201 if (it == yvals.begin())
0202 {
0203 ykey[1] = *it;
0204 if (y < ykey[0])
0205 {
0206 std::cout << PHWHERE << ": This should not happen! y too small - outside range: " << y / cm << std::endl;
0207 return;
0208 }
0209 }
0210 else
0211 {
0212 --it;
0213 ykey[1] = *it;
0214 }
0215 double zkey[2];
0216 it = zvals.lower_bound(z);
0217 zkey[0] = *it;
0218 if (it == zvals.begin())
0219 {
0220 zkey[1] = *it;
0221 if (z < zkey[0])
0222 {
0223 std::cout << PHWHERE << ": This should not happen! z too small - outside range: " << z / cm << std::endl;
0224 return;
0225 }
0226 }
0227 else
0228 {
0229 --it;
0230 zkey[1] = *it;
0231 }
0232
0233 if (xkey_save != xkey[0] ||
0234 ykey_save != ykey[0] ||
0235 zkey_save != zkey[0])
0236 {
0237 cache_misses++;
0238 xkey_save = xkey[0];
0239 ykey_save = ykey[0];
0240 zkey_save = zkey[0];
0241
0242 std::map<std::tuple<float, float, float>, std::tuple<float, float, float> >::const_iterator magval;
0243 trio key;
0244 for (int i = 0; i < 2; i++)
0245 {
0246 for (int j = 0; j < 2; j++)
0247 {
0248 for (int k = 0; k < 2; k++)
0249 {
0250 key = std::make_tuple(xkey[i], ykey[j], zkey[k]);
0251 magval = fieldmap.find(key);
0252 if (magval == fieldmap.end())
0253 {
0254 std::cout << PHWHERE << " could not locate key in " << filename
0255 << " value: x: " << xkey[i] / cm
0256 << ", y: " << ykey[j] / cm
0257 << ", z: " << zkey[k] / cm << std::endl;
0258 return;
0259 }
0260 bf[i][j][k][0] = std::get<0>(magval->second);
0261 bf[i][j][k][1] = std::get<1>(magval->second);
0262 bf[i][j][k][2] = std::get<2>(magval->second);
0263 if (Verbosity() > 0)
0264 {
0265 const double x_loc = std::get<0>(magval->first);
0266 const double y_loc = std::get<1>(magval->first);
0267 const double z_loc = std::get<2>(magval->first);
0268
0269 std::cout << "read x/y/z: " << x_loc / cm << "/"
0270 << y_loc / cm << "/"
0271 << z_loc / cm << " bx/by/bz: "
0272 << bf[i][j][k][0] / tesla << "/"
0273 << bf[i][j][k][1] / tesla << "/"
0274 << bf[i][j][k][2] / tesla << std::endl;
0275 }
0276 }
0277 }
0278 }
0279 }
0280 else
0281 {
0282 cache_hits++;
0283 }
0284
0285
0286 double xinblock = point[0] - xkey[1];
0287 double yinblock = point[1] - ykey[1];
0288 double zinblock = point[2] - zkey[1];
0289
0290 double fractionx = xinblock / xstepsize;
0291 double fractiony = yinblock / ystepsize;
0292 double fractionz = zinblock / zstepsize;
0293 if (Verbosity() > 0)
0294 {
0295 std::cout << "x/y/z stepsize: " << xstepsize / cm << "/" << ystepsize / cm << "/" << zstepsize / cm << std::endl;
0296 std::cout << "x/y/z inblock: " << xinblock / cm << "/" << yinblock / cm << "/" << zinblock / cm << std::endl;
0297 std::cout << "x/y/z fraction: " << fractionx << "/" << fractiony << "/" << fractionz << std::endl;
0298 }
0299
0300
0301
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312 for (int i = 0; i < 3; i++)
0313 {
0314 Bfield[i] = bf[0][0][0][i] * fractionx * fractiony * fractionz +
0315 bf[1][0][0][i] * (1. - fractionx) * fractiony * fractionz +
0316 bf[0][1][0][i] * fractionx * (1. - fractiony) * fractionz +
0317 bf[0][0][1][i] * fractionx * fractiony * (1. - fractionz) +
0318 bf[1][0][1][i] * (1. - fractionx) * fractiony * (1. - fractionz) +
0319 bf[0][1][1][i] * fractionx * (1. - fractiony) * (1. - fractionz) +
0320 bf[1][1][0][i] * (1. - fractionx) * (1. - fractiony) * fractionz +
0321 bf[1][1][1][i] * (1. - fractionx) * (1. - fractiony) * (1. - fractionz);
0322 }
0323
0324 return;
0325 }
0326
0327
0328 void PHField3DCartesian::GetFieldValue_nocache(const double point[4], double *Bfield) const
0329 {
0330
0331 const double& x = point[0];
0332 const double& y = point[1];
0333 const double& z = point[2];
0334
0335 Bfield[0] = 0.0;
0336 Bfield[1] = 0.0;
0337 Bfield[2] = 0.0;
0338 if (!std::isfinite(x) || !std::isfinite(y) || !std::isfinite(z))
0339 { return; }
0340
0341 if (point[0] < xmin || point[0] > xmax ||
0342 point[1] < ymin || point[1] > ymax ||
0343 point[2] < zmin || point[2] > zmax)
0344 { return; }
0345
0346 double xkey[2];
0347 std::set<float>::const_iterator it = xvals.lower_bound(x);
0348 xkey[0] = *it;
0349 if (it == xvals.begin())
0350 {
0351 xkey[1] = *it;
0352 if (x < xkey[0])
0353 {
0354 std::cout << PHWHERE << ": This should not happen! x too small - outside range: " << x / cm << std::endl;
0355 return;
0356 }
0357 }
0358 else
0359 {
0360 --it;
0361 xkey[1] = *it;
0362 }
0363
0364 double ykey[2];
0365 it = yvals.lower_bound(y);
0366 ykey[0] = *it;
0367 if (it == yvals.begin())
0368 {
0369 ykey[1] = *it;
0370 if (y < ykey[0])
0371 {
0372 std::cout << PHWHERE << ": This should not happen! y too small - outside range: " << y / cm << std::endl;
0373 return;
0374 }
0375 }
0376 else
0377 {
0378 --it;
0379 ykey[1] = *it;
0380 }
0381 double zkey[2];
0382 it = zvals.lower_bound(z);
0383 zkey[0] = *it;
0384 if (it == zvals.begin())
0385 {
0386 zkey[1] = *it;
0387 if (z < zkey[0])
0388 {
0389 std::cout << PHWHERE << ": This should not happen! z too small - outside range: " << z / cm << std::endl;
0390 return;
0391 }
0392 }
0393 else
0394 {
0395 --it;
0396 zkey[1] = *it;
0397 }
0398
0399
0400 double bf_loc[2][2][2][3]{};
0401
0402 std::map<std::tuple<float, float, float>, std::tuple<float, float, float> >::const_iterator magval;
0403 trio key;
0404 for (int i = 0; i < 2; i++)
0405 {
0406 for (int j = 0; j < 2; j++)
0407 {
0408 for (int k = 0; k < 2; k++)
0409 {
0410 key = std::make_tuple(xkey[i], ykey[j], zkey[k]);
0411 magval = fieldmap.find(key);
0412 if (magval == fieldmap.end())
0413 {
0414 std::cout << PHWHERE << " could not locate key in " << filename
0415 << " value: x: " << xkey[i] / cm
0416 << ", y: " << ykey[j] / cm
0417 << ", z: " << zkey[k] / cm << std::endl;
0418 return;
0419 }
0420
0421 bf_loc[i][j][k][0] = std::get<0>(magval->second);
0422 bf_loc[i][j][k][1] = std::get<1>(magval->second);
0423 bf_loc[i][j][k][2] = std::get<2>(magval->second);
0424 if (Verbosity() > 0)
0425 {
0426
0427 const double x_loc = std::get<0>(magval->first);
0428 const double y_loc = std::get<1>(magval->first);
0429 const double z_loc = std::get<2>(magval->first);
0430
0431 std::cout << "read x/y/z: " << x_loc / cm << "/"
0432 << y_loc / cm << "/"
0433 << z_loc / cm << " bx/by/bz: "
0434 << bf_loc[i][j][k][0] / tesla << "/"
0435 << bf_loc[i][j][k][1] / tesla << "/"
0436 << bf_loc[i][j][k][2] / tesla << std::endl;
0437 }
0438 }
0439 }
0440 }
0441
0442
0443 double xinblock = point[0] - xkey[1];
0444 double yinblock = point[1] - ykey[1];
0445 double zinblock = point[2] - zkey[1];
0446
0447 double fractionx = xinblock / xstepsize;
0448 double fractiony = yinblock / ystepsize;
0449 double fractionz = zinblock / zstepsize;
0450 if (Verbosity() > 0)
0451 {
0452 std::cout << "x/y/z stepsize: " << xstepsize / cm << "/" << ystepsize / cm << "/" << zstepsize / cm << std::endl;
0453 std::cout << "x/y/z inblock: " << xinblock / cm << "/" << yinblock / cm << "/" << zinblock / cm << std::endl;
0454 std::cout << "x/y/z fraction: " << fractionx << "/" << fractiony << "/" << fractionz << std::endl;
0455 }
0456
0457
0458
0459
0460
0461
0462
0463
0464
0465
0466
0467
0468
0469 for (int i = 0; i < 3; i++)
0470 {
0471 Bfield[i] = bf_loc[0][0][0][i] * fractionx * fractiony * fractionz +
0472 bf_loc[1][0][0][i] * (1. - fractionx) * fractiony * fractionz +
0473 bf_loc[0][1][0][i] * fractionx * (1. - fractiony) * fractionz +
0474 bf_loc[0][0][1][i] * fractionx * fractiony * (1. - fractionz) +
0475 bf_loc[1][0][1][i] * (1. - fractionx) * fractiony * (1. - fractionz) +
0476 bf_loc[0][1][1][i] * fractionx * (1. - fractiony) * (1. - fractionz) +
0477 bf_loc[1][1][0][i] * (1. - fractionx) * (1. - fractiony) * fractionz +
0478 bf_loc[1][1][1][i] * (1. - fractionx) * (1. - fractiony) * (1. - fractionz);
0479 }
0480
0481 return;
0482 }
