Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:20:18 

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