Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:17:57 

0001 #include "PHField2D.h"
0002 
0003 // root framework
0004 #include <TDirectory.h>
0005 #include <TFile.h>
0006 #include <TNtuple.h>
0007 #include <TSystem.h>
0008 
0009 #include <Geant4/G4SystemOfUnits.hh>
0010 
0011 #include <algorithm>
0012 #include <cmath>
0013 #include <cstdlib>
0014 #include <iostream>
0015 #include <iterator>
0016 #include <set>
0017 #include <utility>
0018 
0019 PHField2D::PHField2D(const std::string &filename, const int verb, const float magfield_rescale)
0020   : PHField(verb)
0021   , r_index0_cache(0)
0022   , r_index1_cache(0)
0023   , z_index0_cache(0)
0024   , z_index1_cache(0)
0025 {
0026   if (Verbosity() > 0)
0027   {
0028     std::cout << " ------------- PHField2D::PHField2D() ------------------" << std::endl;
0029   }
0030   // open file
0031   TFile *rootinput = TFile::Open(filename.c_str());
0032   if (!rootinput)
0033   {
0034     std::cout << " could not open " << filename << " exiting now" << std::endl;
0035     gSystem->Exit(1);
0036     exit(1);
0037   }
0038   if (Verbosity() > 0)
0039   {
0040     std::cout << "  Field grid file: " << filename << std::endl;
0041   }
0042   rootinput->cd();
0043 
0044   Float_t ROOT_Z, ROOT_R;
0045   Float_t ROOT_BZ, ROOT_BR;
0046   //  get root NTuple objects
0047   TNtuple *field_map = (TNtuple *) gDirectory->Get("fieldmap");
0048   if (field_map)
0049   {
0050     magfield_unit = tesla;
0051   }
0052   else
0053   {
0054     field_map = (TNtuple *) gDirectory->Get("map");
0055     if (!field_map)
0056     {
0057       std::cout << "PHField2D: could not locate ntuple of name map or fieldmap, exiting now" << std::endl;
0058       exit(1);
0059     }
0060     magfield_unit = gauss;
0061   }
0062   field_map->SetBranchAddress("z", &ROOT_Z);
0063   field_map->SetBranchAddress("r", &ROOT_R);
0064   field_map->SetBranchAddress("bz", &ROOT_BZ);
0065   field_map->SetBranchAddress("br", &ROOT_BR);
0066 
0067   // get the number of entries in the tree
0068   int nz, nr;
0069   nz = field_map->GetEntries("z>-1e6");
0070   nr = field_map->GetEntries("r>-1e6");
0071   static const int NENTRIES = field_map->GetEntries();
0072 
0073   // run checks on entries
0074   if (Verbosity() > 0)
0075   {
0076     std::cout << "  The field grid contained " << NENTRIES << " entries" << std::endl;
0077   }
0078   if (Verbosity() > 1)
0079   {
0080     std::cout << "\n  NENTRIES should be the same as the following values:"
0081               << "\n  [ Number of values r,z: "
0082               << nr << " " << nz << " ]! " << std::endl;
0083   }
0084 
0085   if (nz != nr)
0086   {
0087     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
0088               << "\n The file you entered is not a \"table\" of values"
0089               << "\n Something very likely went oh so wrong"
0090               << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0091   }
0092 
0093   // Keep track of the unique z, r, phi values in the grid using sets
0094   std::set<float> z_set, r_set;
0095 
0096   // Sort the entries to get rid of any stupid ordering problems...
0097 
0098   // We copy the TNtuple into a std::map (which is always sorted)
0099   // using a 3-tuple of (z, r, phi) so it is sorted in z, then r, then
0100   // phi.
0101   if (Verbosity() > 1)
0102   {
0103     std::cout << "  --> Sorting Entries..." << std::endl;
0104   }
0105   std::map<trio, trio> sorted_map;
0106   for (int i = 0; i < field_map->GetEntries(); i++)
0107   {
0108     field_map->GetEntry(i);
0109     trio coord_key(ROOT_Z, ROOT_R);
0110     trio field_val(ROOT_BZ, ROOT_BR);
0111     sorted_map[coord_key] = field_val;
0112 
0113     z_set.insert(ROOT_Z * cm);
0114     r_set.insert(ROOT_R * cm);
0115   }
0116 
0117   // std::couts for assurance
0118   if (Verbosity() > 4)
0119   {
0120     std::map<trio, trio>::iterator it = sorted_map.begin();
0121     print_map(it);
0122     float last_z = std::get<0>(it->first);
0123     for (it = sorted_map.begin(); it != sorted_map.end(); ++it)
0124     {
0125       if (std::get<0>(it->first) != last_z)
0126       {
0127         last_z = std::get<0>(it->first);
0128         print_map(it);
0129       }
0130     }
0131   }
0132 
0133   if (Verbosity() > 1)
0134   {
0135     std::cout << "  --> Putting entries into containers... " << std::endl;
0136   }
0137 
0138   // grab the minimum and maximum z values
0139   minz_ = *(z_set.begin());
0140   std::set<float>::iterator ziter = z_set.end();
0141   --ziter;
0142   maxz_ = *ziter;
0143 
0144   // initialize maps
0145   nz = z_set.size();
0146   nr = r_set.size();
0147   r_map_.resize(nr, 0);
0148   z_map_.resize(nz, 0);
0149 
0150   std::copy(z_set.begin(), z_set.end(), z_map_.begin());
0151   std::copy(r_set.begin(), r_set.end(), r_map_.begin());
0152 
0153   // initialize the field map vectors to the correct sizes
0154   BFieldR_.resize(nz, std::vector<float>(nr, 0));
0155   BFieldZ_.resize(nz, std::vector<float>(nr, 0));
0156 
0157   // all of this assumes that  z_prev < z , i.e. the table is ordered (as of right now)
0158   unsigned int ir = 0, iz = 0;  // useful indexes to keep track of
0159   std::map<trio, trio>::iterator iter = sorted_map.begin();
0160   for (; iter != sorted_map.end(); ++iter)
0161   {
0162     // equivalent to ->GetEntry(iter)
0163     float z = std::get<0>(iter->first) * cm;
0164     float r = std::get<1>(iter->first) * cm;
0165     float Bz = std::get<0>(iter->second) * magfield_unit;
0166     float Br = std::get<1>(iter->second) * magfield_unit;
0167 
0168     if (z > maxz_)
0169     {
0170       maxz_ = z;
0171     }
0172     if (z < minz_)
0173     {
0174       minz_ = z;
0175     }
0176 
0177     // check for change in z value, when z changes we have a ton of updates to do
0178     if (z != z_map_[iz])
0179     {
0180       ++iz;
0181       ir = 0;
0182     }
0183     else if (r != r_map_[ir])  // check for change in r value
0184     {
0185       ++ir;
0186     }
0187 
0188     // shouldn't happen
0189     if (iz > 0 && z < z_map_[iz - 1])
0190     {
0191       std::cout << "!!!!!!!!! Your map isn't ordered.... z: " << z << " zprev: " << z_map_[iz - 1] << std::endl;
0192     }
0193 
0194     BFieldR_[iz][ir] = Br * magfield_rescale;
0195     BFieldZ_[iz][ir] = Bz * magfield_rescale;
0196 
0197     // you can change this to check table values for correctness
0198     // print_map prints the values in the root table, and the
0199     // std::couts print the values entered into the vectors
0200     if (std::fabs(z) < 10 && ir < 10 /*&& iphi==2*/ && Verbosity() > 3)
0201     {
0202       print_map(iter);
0203 
0204       std::cout << " B("
0205                 << r_map_[ir] << ", "
0206                 << z_map_[iz] << "):  ("
0207                 << BFieldR_[iz][ir] << ", "
0208                 << BFieldZ_[iz][ir] << ")" << std::endl;
0209     }
0210 
0211   }  // end loop over root field map file
0212 
0213   rootinput->Close();
0214 
0215   if (Verbosity() > 0)
0216   {
0217     std::cout << "  Mag field z boundaries (min,max): (" << minz_ / cm << ", " << maxz_ / cm << ") cm" << std::endl;
0218   }
0219   if (Verbosity() > 0)
0220   {
0221     std::cout << "  Mag field r max boundary: " << r_map_.back() / cm << " cm" << std::endl;
0222   }
0223 
0224   if (Verbosity() > 0)
0225   {
0226     std::cout << " -----------------------------------------------------------" << std::endl;
0227   }
0228 }
0229 
0230 void PHField2D::GetFieldValue(const double point[4], double *Bfield) const
0231 {
0232   if (Verbosity() > 2)
0233   {
0234     std::cout << "\nPHField2D::GetFieldValue" << std::endl;
0235   }
0236   double x = point[0];
0237   double y = point[1];
0238   double z = point[2];
0239   double r = sqrt(x * x + y * y);
0240   double phi;
0241   phi = atan2(y, x);
0242   if (phi < 0)
0243   {
0244     phi += 2 * M_PI;  // normalize phi to be over the range [0,2*pi]
0245   }
0246 
0247   // Check that the point is within the defined z region (check r in a second)
0248   if ((z >= minz_) && (z <= maxz_))
0249   {
0250     double BFieldCyl[3];
0251     double cylpoint[4] = {z, r, phi, 0};
0252 
0253     // take <z,r,phi> location and return a vector of <Bz, Br, Bphi>
0254     GetFieldCyl(cylpoint, BFieldCyl);
0255 
0256     // X direction of B-field ( Bx = Br*cos(phi) - Bphi*sin(phi)
0257     Bfield[0] = cos(phi) * BFieldCyl[1] - sin(phi) * BFieldCyl[2];  // unit vector transformations
0258 
0259     // Y direction of B-field ( By = Br*sin(phi) + Bphi*cos(phi)
0260     Bfield[1] = sin(phi) * BFieldCyl[1] + cos(phi) * BFieldCyl[2];
0261 
0262     // Z direction of B-field
0263     Bfield[2] = BFieldCyl[0];
0264   }
0265   else  // x,y,z is outside of z range of the field map
0266   {
0267     Bfield[0] = 0.0;
0268     Bfield[1] = 0.0;
0269     Bfield[2] = 0.0;
0270     if (Verbosity() > 2)
0271     {
0272       std::cout << "!!!!!!!!!! Field point not in defined region (outside of z bounds)" << std::endl;
0273     }
0274   }
0275 
0276   if (Verbosity() > 2)
0277   {
0278     std::cout << "END PHField2D::GetFieldValue\n"
0279               << "  --->  {Bx, By, Bz} : "
0280               << "< " << Bfield[0] << ", " << Bfield[1] << ", " << Bfield[2] << " >" << std::endl;
0281   }
0282 
0283   return;
0284 }
0285 
0286 void PHField2D::GetFieldValue_nocache(const double point[4], double *Bfield) const
0287 {
0288   if (Verbosity() > 2)
0289   {
0290     std::cout << "\nPHField2D::GetFieldValue" << std::endl;
0291   }
0292   double x = point[0];
0293   double y = point[1];
0294   double z = point[2];
0295   double r = sqrt(x * x + y * y);
0296   double phi;
0297   phi = atan2(y, x);
0298   if (phi < 0)
0299   {
0300     phi += 2 * M_PI;  // normalize phi to be over the range [0,2*pi]
0301   }
0302 
0303   // Check that the point is within the defined z region (check r in a second)
0304   if ((z >= minz_) && (z <= maxz_))
0305   {
0306     double BFieldCyl[3];
0307     double cylpoint[4] = {z, r, phi, 0};
0308 
0309     // take <z,r,phi> location and return a vector of <Bz, Br, Bphi>
0310     GetFieldCyl_nocache(cylpoint, BFieldCyl);
0311 
0312     // X direction of B-field ( Bx = Br*cos(phi) - Bphi*sin(phi)
0313     Bfield[0] = cos(phi) * BFieldCyl[1] - sin(phi) * BFieldCyl[2];  // unit vector transformations
0314 
0315     // Y direction of B-field ( By = Br*sin(phi) + Bphi*cos(phi)
0316     Bfield[1] = sin(phi) * BFieldCyl[1] + cos(phi) * BFieldCyl[2];
0317 
0318     // Z direction of B-field
0319     Bfield[2] = BFieldCyl[0];
0320   }
0321   else  // x,y,z is outside of z range of the field map
0322   {
0323     Bfield[0] = 0.0;
0324     Bfield[1] = 0.0;
0325     Bfield[2] = 0.0;
0326     if (Verbosity() > 2)
0327     {
0328       std::cout << "!!!!!!!!!! Field point not in defined region (outside of z bounds)" << std::endl;
0329     }
0330   }
0331 
0332   if (Verbosity() > 2)
0333   {
0334     std::cout << "END PHField2D::GetFieldValue\n"
0335               << "  --->  {Bx, By, Bz} : "
0336               << "< " << Bfield[0] << ", " << Bfield[1] << ", " << Bfield[2] << " >" << std::endl;
0337   }
0338 
0339   return;
0340 }
0341 
0342 void PHField2D::GetFieldCyl(const double CylPoint[4], double *BfieldCyl) const
0343 {
0344   float z = CylPoint[0];
0345   float r = CylPoint[1];
0346 
0347   BfieldCyl[0] = 0.0;
0348   BfieldCyl[1] = 0.0;
0349   BfieldCyl[2] = 0.0;
0350 
0351   if (Verbosity() > 2)
0352   {
0353     std::cout << "GetFieldCyl@ <z,r>: {" << z << "," << r << "}" << std::endl;
0354   }
0355 
0356   if (z < z_map_[0] || z > z_map_[z_map_.size() - 1])
0357   {
0358     if (Verbosity() > 2)
0359     {
0360       std::cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
0361     }
0362     return;
0363   }
0364 
0365   // since GEANT4 looks up the field ~95% of the time in the same voxel
0366   // between subsequent calls, we can save on the expense of the upper_bound
0367   // lookup (~10-15% of central event run time) with some caching between calls
0368   unsigned int r_index0 = r_index0_cache;
0369   unsigned int r_index1 = r_index1_cache;
0370 
0371   if (!((r > r_map_[r_index0]) && (r < r_map_[r_index1])))
0372   {
0373     // if miss cached r values, search through the lookup table
0374     std::vector<float>::const_iterator riter = upper_bound(r_map_.begin(), r_map_.end(), r);
0375     r_index0 = distance(r_map_.begin(), riter) - 1;
0376     if (r_index0 >= r_map_.size())
0377     {
0378       if (Verbosity() > 2)
0379       {
0380         std::cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
0381       }
0382       return;
0383     }
0384 
0385     r_index1 = r_index0 + 1;
0386     if (r_index1 >= r_map_.size())
0387     {
0388       if (Verbosity() > 2)
0389       {
0390         std::cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
0391       }
0392       return;
0393     }
0394 
0395     // update cache
0396     r_index0_cache = r_index0;
0397     r_index1_cache = r_index1;
0398   }
0399 
0400   unsigned int z_index0 = z_index0_cache;
0401   unsigned int z_index1 = z_index1_cache;
0402 
0403   if (!((z > z_map_[z_index0]) && (z < z_map_[z_index1])))
0404   {
0405     // if miss cached z values, search through the lookup table
0406     std::vector<float>::const_iterator ziter = upper_bound(z_map_.begin(), z_map_.end(), z);
0407     z_index0 = distance(z_map_.begin(), ziter) - 1;
0408     z_index1 = z_index0 + 1;
0409     if (z_index1 >= z_map_.size())
0410     {
0411       if (Verbosity() > 2)
0412       {
0413         std::cout << "!!!! Point not in defined region (z too large in specific r-plane)" << std::endl;
0414       }
0415       return;
0416     }
0417 
0418     // update cache
0419     z_index0_cache = z_index0;
0420     z_index1_cache = z_index1;
0421   }
0422 
0423   double Br000 = BFieldR_[z_index0][r_index0];
0424   double Br010 = BFieldR_[z_index0][r_index1];
0425   double Br100 = BFieldR_[z_index1][r_index0];
0426   double Br110 = BFieldR_[z_index1][r_index1];
0427 
0428   double Bz000 = BFieldZ_[z_index0][r_index0];
0429   double Bz100 = BFieldZ_[z_index1][r_index0];
0430   double Bz010 = BFieldZ_[z_index0][r_index1];
0431   double Bz110 = BFieldZ_[z_index1][r_index1];
0432 
0433   double zweight = z - z_map_[z_index0];
0434   double zspacing = z_map_[z_index1] - z_map_[z_index0];
0435   zweight /= zspacing;
0436 
0437   double rweight = r - r_map_[r_index0];
0438   double rspacing = r_map_[r_index1] - r_map_[r_index0];
0439   rweight /= rspacing;
0440 
0441   // Z direction of B-field
0442   BfieldCyl[0] =
0443       (1 - zweight) * ((1 - rweight) * Bz000 +
0444                        rweight * Bz010) +
0445       zweight * ((1 - rweight) * Bz100 +
0446                  rweight * Bz110);
0447 
0448   // R direction of B-field
0449   BfieldCyl[1] =
0450       (1 - zweight) * ((1 - rweight) * Br000 +
0451                        rweight * Br010) +
0452       zweight * ((1 - rweight) * Br100 +
0453                  rweight * Br110);
0454 
0455   // PHI Direction of B-field
0456   BfieldCyl[2] = 0;
0457 
0458   if (Verbosity() > 2)
0459   {
0460     std::cout << "End GFCyl Call: <bz,br,bphi> : {"
0461               << BfieldCyl[0] / gauss << "," << BfieldCyl[1] / gauss << "," << BfieldCyl[2] / gauss << "}"
0462               << std::endl;
0463   }
0464 
0465   return;
0466 }
0467 
0468 
0469 void PHField2D::GetFieldCyl_nocache(const double CylPoint[4], double *BfieldCyl) const
0470 {
0471   float z = CylPoint[0];
0472   float r = CylPoint[1];
0473 
0474   BfieldCyl[0] = 0.0;
0475   BfieldCyl[1] = 0.0;
0476   BfieldCyl[2] = 0.0;
0477 
0478   if (Verbosity() > 2)
0479   {
0480     std::cout << "GetFieldCyl@ <z,r>: {" << z << "," << r << "}" << std::endl;
0481   }
0482 
0483   if (z < z_map_[0] || z > z_map_[z_map_.size() - 1])
0484   {
0485     if (Verbosity() > 2)
0486     {
0487       std::cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
0488     }
0489     return;
0490   }
0491 
0492   // since GEANT4 looks up the field ~95% of the time in the same voxel
0493   // between subsequent calls, we can save on the expense of the upper_bound
0494   // lookup (~10-15% of central event run time) with some caching between calls
0495 
0496   // if miss cached r values, search through the lookup table
0497   std::vector<float>::const_iterator riter = upper_bound(r_map_.begin(), r_map_.end(), r);
0498   const unsigned int r_index0 = distance(r_map_.begin(), riter) - 1;
0499   if (r_index0 >= r_map_.size())
0500   {
0501     if (Verbosity() > 2)
0502     {
0503       std::cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
0504     }
0505     return;
0506   }
0507 
0508   const unsigned int r_index1 = r_index0 + 1;
0509   if (r_index1 >= r_map_.size())
0510   {
0511     if (Verbosity() > 2)
0512     {
0513       std::cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
0514     }
0515     return;
0516   }
0517 
0518   // if miss cached z values, search through the lookup table
0519   std::vector<float>::const_iterator ziter = upper_bound(z_map_.begin(), z_map_.end(), z);
0520   const unsigned int z_index0 = distance(z_map_.begin(), ziter) - 1;
0521   const unsigned int z_index1 = z_index0 + 1;
0522   if (z_index1 >= z_map_.size())
0523   {
0524     if (Verbosity() > 2)
0525     {
0526       std::cout << "!!!! Point not in defined region (z too large in specific r-plane)" << std::endl;
0527     }
0528     return;
0529   }
0530 
0531   double Br000 = BFieldR_[z_index0][r_index0];
0532   double Br010 = BFieldR_[z_index0][r_index1];
0533   double Br100 = BFieldR_[z_index1][r_index0];
0534   double Br110 = BFieldR_[z_index1][r_index1];
0535 
0536   double Bz000 = BFieldZ_[z_index0][r_index0];
0537   double Bz100 = BFieldZ_[z_index1][r_index0];
0538   double Bz010 = BFieldZ_[z_index0][r_index1];
0539   double Bz110 = BFieldZ_[z_index1][r_index1];
0540 
0541   double zweight = z - z_map_[z_index0];
0542   double zspacing = z_map_[z_index1] - z_map_[z_index0];
0543   zweight /= zspacing;
0544 
0545   double rweight = r - r_map_[r_index0];
0546   double rspacing = r_map_[r_index1] - r_map_[r_index0];
0547   rweight /= rspacing;
0548 
0549   // Z direction of B-field
0550   BfieldCyl[0] =
0551       (1 - zweight) * ((1 - rweight) * Bz000 +
0552                        rweight * Bz010) +
0553       zweight * ((1 - rweight) * Bz100 +
0554                  rweight * Bz110);
0555 
0556   // R direction of B-field
0557   BfieldCyl[1] =
0558       (1 - zweight) * ((1 - rweight) * Br000 +
0559                        rweight * Br010) +
0560       zweight * ((1 - rweight) * Br100 +
0561                  rweight * Br110);
0562 
0563   // PHI Direction of B-field
0564   BfieldCyl[2] = 0;
0565 
0566   if (Verbosity() > 2)
0567   {
0568     std::cout << "End GFCyl Call: <bz,br,bphi> : {"
0569               << BfieldCyl[0] / gauss << "," << BfieldCyl[1] / gauss << "," << BfieldCyl[2] / gauss << "}"
0570               << std::endl;
0571   }
0572 
0573   return;
0574 }
0575 
0576 // debug function to print key/value pairs in map
0577 void PHField2D::print_map(std::map<trio, trio>::iterator &it) const
0578 {
0579   std::cout << "    Key: <"
0580             << std::get<0>(it->first) / cm << ","
0581             << std::get<1>(it->first) / cm << ">"
0582 
0583             << " Value: <"
0584             << std::get<0>(it->second) / magfield_unit << ","
0585             << std::get<1>(it->second) / magfield_unit << ">" << std::endl;
0586 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!