packages/PHField/PHField3DCartesian.cc

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