Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "PHField3DCylindrical.h"
0002 
0003 #include <TDirectory.h>  // for TDirectory, gDirectory
0004 #include <TFile.h>
0005 #include <TNtuple.h>
0006 
0007 #include <Geant4/G4SystemOfUnits.hh>
0008 
0009 #include <algorithm>
0010 #include <cassert>
0011 #include <cmath>
0012 #include <cstdlib>
0013 #include <iostream>
0014 #include <iterator>
0015 #include <set>
0016 #include <utility>
0017 
0018 PHField3DCylindrical::PHField3DCylindrical(const std::string &filename, const int verb, const float magfield_rescale)
0019   : PHField(verb)
0020 {
0021   std::cout << "\n================ Begin Construct Mag Field =====================" << std::endl;
0022   std::cout << "\n-----------------------------------------------------------"
0023             << "\n      Magnetic field Module - Verbosity:" << Verbosity()
0024             << "\n-----------------------------------------------------------";
0025 
0026   // open file
0027   TFile *rootinput = TFile::Open(filename.c_str());
0028   if (!rootinput)
0029   {
0030     std::cout << "\n could not open " << filename << " exiting now" << std::endl;
0031     exit(1);
0032   }
0033   std::cout << "\n ---> "
0034                "Reading the field grid from "
0035             << filename << " ... " << std::endl;
0036   rootinput->cd();
0037 
0038   //  get root NTuple objects
0039   TNtuple *field_map = (TNtuple *) gDirectory->Get("map");
0040   Float_t ROOT_Z;
0041   Float_t ROOT_R;
0042   Float_t ROOT_PHI;
0043   Float_t ROOT_BZ;
0044   Float_t ROOT_BR;
0045   Float_t ROOT_BPHI;
0046   field_map->SetBranchAddress("z", &ROOT_Z);
0047   field_map->SetBranchAddress("r", &ROOT_R);
0048   field_map->SetBranchAddress("phi", &ROOT_PHI);
0049   field_map->SetBranchAddress("bz", &ROOT_BZ);
0050   field_map->SetBranchAddress("br", &ROOT_BR);
0051   field_map->SetBranchAddress("bphi", &ROOT_BPHI);
0052 
0053   // get the number of entries in the tree
0054   int nz;
0055   int nr;
0056   int nphi;
0057   nz = field_map->GetEntries("z>-1e6");
0058   nr = field_map->GetEntries("r>-1e6");
0059   nphi = field_map->GetEntries("phi>-1e6");
0060   static const int NENTRIES = field_map->GetEntries();
0061 
0062   // run checks on entries
0063   std::cout << " ---> The field grid contained " << NENTRIES << " entries" << std::endl;
0064   if (Verbosity() > 0)
0065   {
0066     std::cout << "\n  NENTRIES should be the same as the following values:"
0067               << "\n  [ Number of values r,phi,z: "
0068               << nr << " " << nphi << " " << nz << " ]! " << std::endl;
0069   }
0070 
0071   if (nz != nr || nz != nphi || nr != nphi)
0072   {
0073     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
0074               << "\n The file you entered is not a \"table\" of values"
0075               << "\n Something very likely went oh so wrong"
0076               << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0077   }
0078 
0079   // Keep track of the unique z, r, phi values in the grid using sets
0080   std::set<float> z_set;
0081   std::set<float> r_set;
0082   std::set<float> phi_set;
0083 
0084   // Sort the entries to get rid of any stupid ordering problems...
0085 
0086   // We copy the TNtuple into a std::map (which is always sorted)
0087   // using a 3-tuple of (z, r, phi) so it is sorted in z, then r, then
0088   // phi.
0089   if (Verbosity() > 0)
0090   {
0091     std::cout << "  --> Sorting Entries..." << std::endl;
0092   }
0093   std::map<trio, trio> sorted_map;
0094   for (int i = 0; i < field_map->GetEntries(); i++)
0095   {
0096     field_map->GetEntry(i);
0097     trio coord_key(ROOT_Z, ROOT_R, ROOT_PHI);
0098     trio field_val(ROOT_BZ, ROOT_BR, ROOT_BPHI);
0099     sorted_map[coord_key] = field_val;
0100 
0101     z_set.insert(ROOT_Z * cm);
0102     r_set.insert(ROOT_R * cm);
0103     phi_set.insert(ROOT_PHI * deg);
0104   }
0105 
0106   // std::couts for assurance
0107   if (Verbosity() > 4)
0108   {
0109     std::map<trio, trio>::iterator it = sorted_map.begin();
0110     print_map(it);
0111     float last_z = std::get<0>(it->first);
0112     for (it = sorted_map.begin(); it != sorted_map.end(); ++it)
0113     {
0114       if (std::get<0>(it->first) != last_z)
0115       {
0116         last_z = std::get<0>(it->first);
0117         print_map(it);
0118       }
0119     }
0120   }
0121 
0122   if (Verbosity() > 0)
0123   {
0124     std::cout << "  --> Putting entries into containers... " << std::endl;
0125   }
0126 
0127   // grab the minimum and maximum z values
0128   minz_ = *(z_set.begin());
0129   std::set<float>::iterator ziter = z_set.end();
0130   --ziter;
0131   maxz_ = *ziter;
0132 
0133   // initialize maps
0134   nz = z_set.size();
0135   nr = r_set.size();
0136   nphi = phi_set.size();
0137   r_map_.resize(nr, 0);
0138   z_map_.resize(nz, 0);
0139   phi_map_.resize(nphi, 0);
0140 
0141   std::copy(z_set.begin(), z_set.end(), z_map_.begin());
0142   std::copy(phi_set.begin(), phi_set.end(), phi_map_.begin());
0143   std::copy(r_set.begin(), r_set.end(), r_map_.begin());
0144 
0145   // initialize the field map vectors to the correct sizes
0146   BFieldR_.resize(nz, std::vector<std::vector<float> >(nr, std::vector<float>(nphi, 0)));
0147   BFieldPHI_.resize(nz, std::vector<std::vector<float> >(nr, std::vector<float>(nphi, 0)));
0148   BFieldZ_.resize(nz, std::vector<std::vector<float> >(nr, std::vector<float>(nphi, 0)));
0149 
0150   // all of this assumes that  z_prev < z , i.e. the table is ordered (as of right now)
0151   unsigned int ir = 0;
0152   unsigned int iphi = 0;
0153   unsigned int iz = 0;  // useful indexes to keep track of
0154   std::map<trio, trio>::iterator iter = sorted_map.begin();
0155   for (; iter != sorted_map.end(); ++iter)
0156   {
0157     // equivalent to ->GetEntry(iter)
0158     float z = std::get<0>(iter->first) * cm;
0159     float r = std::get<1>(iter->first) * cm;
0160     float phi = std::get<2>(iter->first) * deg;
0161     float Bz = std::get<0>(iter->second) * gauss;
0162     float Br = std::get<1>(iter->second) * gauss;
0163     float Bphi = std::get<2>(iter->second) * gauss;
0164 
0165     maxz_ = std::max(z, maxz_);
0166     minz_ = std::min(z, minz_);
0167 
0168     // check for change in z value, when z changes we have a ton of updates to do
0169     if (z != z_map_[iz])
0170     {
0171       ++iz;
0172       ir = 0;
0173       iphi = 0;  // reset indices
0174     }
0175     else if (r != r_map_[ir])
0176     {  // check for change in r value
0177       ++ir;
0178       iphi = 0;
0179     }
0180     else if (phi != phi_map_[iphi])
0181     {  // change in phi value? (should be every time)
0182       ++iphi;
0183     }
0184 
0185     // shouldn't happen
0186     if (iz > 0 && z < z_map_[iz - 1])
0187     {
0188       std::cout << "!!!!!!!!! Your map isn't ordered.... z: " << z << " zprev: " << z_map_[iz - 1] << std::endl;
0189     }
0190 
0191     BFieldR_[iz][ir][iphi] = Br * magfield_rescale;
0192     BFieldPHI_[iz][ir][iphi] = Bphi * magfield_rescale;
0193     BFieldZ_[iz][ir][iphi] = Bz * magfield_rescale;
0194 
0195     // you can change this to check table values for correctness
0196     // print_map prints the values in the root table, and the
0197     // std::couts print the values entered into the vectors
0198     if (std::fabs(z) < 10 && ir < 10 /*&& iphi==2*/ && Verbosity() > 3)
0199     {
0200       print_map(iter);
0201 
0202       std::cout << " B("
0203                 << r_map_[ir] << ", "
0204                 << phi_map_[iphi] << ", "
0205                 << z_map_[iz] << "):  ("
0206                 << BFieldR_[iz][ir][iphi] << ", "
0207                 << BFieldPHI_[iz][ir][iphi] << ", "
0208                 << BFieldZ_[iz][ir][iphi] << ")" << std::endl;
0209     }
0210 
0211   }  // end loop over root field map file
0212 
0213   rootinput->Close();
0214 
0215   std::cout << "\n ---> ... read file successfully "
0216             << "\n ---> Z Boundaries ~ zlow, zhigh: "
0217             << minz_ / cm << "," << maxz_ / cm << " cm " << std::endl;
0218 
0219   std::cout << "\n================= End Construct Mag Field ======================\n"
0220             << std::endl;
0221 }
0222 
0223 void PHField3DCylindrical::GetFieldValue(const double point[4], double *Bfield) const
0224 {
0225   if (Verbosity() > 2)
0226   {
0227     std::cout << "\nPHField3DCylindrical::GetFieldValue" << std::endl;
0228   }
0229   double x = point[0];
0230   double y = point[1];
0231   double z = point[2];
0232   double r = sqrt(x * x + y * y);
0233   double phi;
0234   if (x == 0)
0235   {
0236     phi = atan2(y, 0.00000000001);
0237   }
0238   else
0239   {
0240     phi = atan2(y, x);
0241   }
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
0266   {  // x,y,z is outside of z range of the field map
0267 
0268     Bfield[0] = 0.0;
0269     Bfield[1] = 0.0;
0270     Bfield[2] = 0.0;
0271     if (Verbosity() > 2)
0272     {
0273       std::cout << "!!!!!!!!!! Field point not in defined region (outside of z bounds)" << std::endl;
0274     }
0275   }
0276 
0277   if (Verbosity() > 2)
0278   {
0279     std::cout << "END PHField3DCylindrical::GetFieldValue\n"
0280               << "  --->  {Bx, By, Bz} : "
0281               << "< " << Bfield[0] << ", " << Bfield[1] << ", " << Bfield[2] << " >" << std::endl;
0282   }
0283 
0284   return;
0285 }
0286 
0287 void PHField3DCylindrical::GetFieldCyl(const double CylPoint[4], double *BfieldCyl) const
0288 {
0289   float z = CylPoint[0];
0290   float r = CylPoint[1];
0291   float phi = CylPoint[2];
0292 
0293   BfieldCyl[0] = 0.0;
0294   BfieldCyl[1] = 0.0;
0295   BfieldCyl[2] = 0.0;
0296 
0297   if (Verbosity() > 2)
0298   {
0299     std::cout << "GetFieldCyl@ <z,r,phi>: {" << z << "," << r << "," << phi << "}" << std::endl;
0300   }
0301 
0302   if (z <= z_map_[0] || z >= z_map_[z_map_.size() - 1])
0303   {
0304     if (Verbosity() > 2)
0305     {
0306       std::cout << "!!!! Point not in defined region (|z| too large)" << std::endl;
0307     }
0308     return;
0309   }
0310   if (r < r_map_[0])
0311   {
0312     r = r_map_[0];
0313     if (Verbosity() > 2)
0314     {
0315       std::cout << "!!!! Point not in defined region (radius too small in specific z-plane). Use min radius" << std::endl;
0316     }
0317     //    return;
0318   }
0319   if (r > r_map_[r_map_.size() - 1])
0320   {
0321     if (Verbosity() > 2)
0322     {
0323       std::cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
0324     }
0325     return;
0326   }
0327 
0328   std::vector<float>::const_iterator ziter = upper_bound(z_map_.begin(), z_map_.end(), z);
0329   int z_index0 = distance(z_map_.begin(), ziter) - 1;
0330   int z_index1 = z_index0 + 1;
0331 
0332   assert(z_index0 >= 0);
0333   assert(z_index1 >= 0);
0334   assert(z_index0 < (int) z_map_.size());
0335   assert(z_index1 < (int) z_map_.size());
0336 
0337   std::vector<float>::const_iterator riter = upper_bound(r_map_.begin(), r_map_.end(), r);
0338   int r_index0 = distance(r_map_.begin(), riter) - 1;
0339   if (r_index0 >= (int) r_map_.size())
0340   {
0341     if (Verbosity() > 2)
0342     {
0343       std::cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
0344     }
0345     return;
0346   }
0347 
0348   int r_index1 = r_index0 + 1;
0349   if (r_index1 >= (int) r_map_.size())
0350   {
0351     if (Verbosity() > 2)
0352     {
0353       std::cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
0354     }
0355     return;
0356   }
0357 
0358   assert(r_index0 >= 0);
0359   assert(r_index1 >= 0);
0360 
0361   std::vector<float>::const_iterator phiiter = upper_bound(phi_map_.begin(), phi_map_.end(), phi);
0362   int phi_index0 = distance(phi_map_.begin(), phiiter) - 1;
0363   int phi_index1 = phi_index0 + 1;
0364   if (phi_index1 >= (int) phi_map_.size())
0365   {
0366     phi_index1 = 0;
0367   }
0368 
0369   assert(phi_index0 >= 0);
0370   assert(phi_index0 < (int) phi_map_.size());
0371   assert(phi_index1 >= 0);
0372 
0373   double Br000 = BFieldR_[z_index0][r_index0][phi_index0];
0374   double Br001 = BFieldR_[z_index0][r_index0][phi_index1];
0375   double Br010 = BFieldR_[z_index0][r_index1][phi_index0];
0376   double Br011 = BFieldR_[z_index0][r_index1][phi_index1];
0377   double Br100 = BFieldR_[z_index1][r_index0][phi_index0];
0378   double Br101 = BFieldR_[z_index1][r_index0][phi_index1];
0379   double Br110 = BFieldR_[z_index1][r_index1][phi_index0];
0380   double Br111 = BFieldR_[z_index1][r_index1][phi_index1];
0381 
0382   double Bphi000 = BFieldPHI_[z_index0][r_index0][phi_index0];
0383   double Bphi001 = BFieldPHI_[z_index0][r_index0][phi_index1];
0384   double Bphi010 = BFieldPHI_[z_index0][r_index1][phi_index0];
0385   double Bphi011 = BFieldPHI_[z_index0][r_index1][phi_index1];
0386   double Bphi100 = BFieldPHI_[z_index1][r_index0][phi_index0];
0387   double Bphi101 = BFieldPHI_[z_index1][r_index0][phi_index1];
0388   double Bphi110 = BFieldPHI_[z_index1][r_index1][phi_index0];
0389   double Bphi111 = BFieldPHI_[z_index1][r_index1][phi_index1];
0390 
0391   double Bz000 = BFieldZ_[z_index0][r_index0][phi_index0];
0392   double Bz001 = BFieldZ_[z_index0][r_index0][phi_index1];
0393   double Bz100 = BFieldZ_[z_index1][r_index0][phi_index0];
0394   double Bz101 = BFieldZ_[z_index1][r_index0][phi_index1];
0395   double Bz010 = BFieldZ_[z_index0][r_index1][phi_index0];
0396   double Bz110 = BFieldZ_[z_index1][r_index1][phi_index0];
0397   double Bz011 = BFieldZ_[z_index0][r_index1][phi_index1];
0398   double Bz111 = BFieldZ_[z_index1][r_index1][phi_index1];
0399 
0400   double zweight = z - z_map_[z_index0];
0401   double zspacing = z_map_[z_index1] - z_map_[z_index0];
0402   zweight /= zspacing;
0403 
0404   double rweight = r - r_map_[r_index0];
0405   double rspacing = r_map_[r_index1] - r_map_[r_index0];
0406   rweight /= rspacing;
0407 
0408   double phiweight = phi - phi_map_[phi_index0];
0409   double phispacing = phi_map_[phi_index1] - phi_map_[phi_index0];
0410   if (phi_index1 == 0)
0411   {
0412     phispacing += 2 * M_PI;
0413   }
0414   phiweight /= phispacing;
0415 
0416   // Z direction of B-field
0417   BfieldCyl[0] =
0418       (1 - zweight) * ((1 - rweight) * ((1 - phiweight) * Bz000 + phiweight * Bz001) +
0419                        rweight * ((1 - phiweight) * Bz010 + phiweight * Bz011)) +
0420       zweight * ((1 - rweight) * ((1 - phiweight) * Bz100 + phiweight * Bz101) +
0421                  rweight * ((1 - phiweight) * Bz110 + phiweight * Bz111));
0422 
0423   // R direction of B-field
0424   BfieldCyl[1] =
0425       (1 - zweight) * ((1 - rweight) * ((1 - phiweight) * Br000 + phiweight * Br001) +
0426                        rweight * ((1 - phiweight) * Br010 + phiweight * Br011)) +
0427       zweight * ((1 - rweight) * ((1 - phiweight) * Br100 + phiweight * Br101) +
0428                  rweight * ((1 - phiweight) * Br110 + phiweight * Br111));
0429 
0430   // PHI Direction of B-field
0431   BfieldCyl[2] =
0432       (1 - zweight) * ((1 - rweight) * ((1 - phiweight) * Bphi000 + phiweight * Bphi001) +
0433                        rweight * ((1 - phiweight) * Bphi010 + phiweight * Bphi011)) +
0434       zweight * ((1 - rweight) * ((1 - phiweight) * Bphi100 + phiweight * Bphi101) +
0435                  rweight * ((1 - phiweight) * Bphi110 + phiweight * Bphi111));
0436 
0437   //     std::cout << "wr: " << rweight << " wz: " << zweight << " wphi: " << phiweight << std::endl;
0438   //     std::cout << "Bz000: " << Bz000 << std::endl
0439   //          << "Bz001: " << Bz001 << std::endl
0440   //          << "Bz010: " << Bz010 << std::endl
0441   //          << "Bz011: " << Bz011 << std::endl
0442   //          << "Bz100: " << Bz100 << std::endl
0443   //          << "Bz101: " << Bz101 << std::endl
0444   //          << "Bz110: " << Bz110 << std::endl
0445   //          << "Bz111: " << Bz111 << std::endl
0446   //          << "Bz:    " << BfieldCyl[0] << std::endl << std::endl;
0447 
0448   if (Verbosity() > 2)
0449   {
0450     std::cout << "End GFCyl Call: <bz,br,bphi> : {"
0451               << BfieldCyl[0] / gauss << "," << BfieldCyl[1] / gauss << "," << BfieldCyl[2] / gauss << "}"
0452               << std::endl;
0453   }
0454 
0455   return;
0456 }
0457 
0458 // a binary search algorithm that puts the location that "key" would be, into index...
0459 // it returns true if key was found, and false if not.
0460 bool PHField3DCylindrical::bin_search(const std::vector<float> &vec, unsigned start, unsigned end, const float &key, unsigned &index) const
0461 {
0462   // Termination condition: start index greater than end index
0463   if (start > end)
0464   {
0465     index = start;
0466     return false;
0467   }
0468 
0469   // Find the middle element of the vector and use that for splitting
0470   // the array into two pieces.
0471   unsigned middle = start + ((end - start) / 2);
0472 
0473   if (vec[middle] == key)
0474   {
0475     index = middle;
0476     return true;
0477   }
0478   if (vec[middle] > key)
0479   {
0480     return bin_search(vec, start, middle - 1, key, index);
0481   }
0482   return bin_search(vec, middle + 1, end, key, index);
0483 }
0484 
0485 // debug function to print key/value pairs in map
0486 void PHField3DCylindrical::print_map(std::map<trio, trio>::iterator &it) const
0487 {
0488   std::cout << "    Key: <"
0489             << std::get<0>(it->first) * cm << ","
0490             << std::get<1>(it->first) * cm << ","
0491             << std::get<2>(it->first) * deg << ">"
0492 
0493             << " Value: <"
0494             << std::get<0>(it->second) * gauss << ","
0495             << std::get<1>(it->second) * gauss << ","
0496             << std::get<2>(it->second) * gauss << ">" << std::endl;
0497 }
[ 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!