Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DMagneticFieldMapPS2DMap.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: DMagneticFieldMapPS2DMap.cc
4 
5 #include <sys/stat.h>
6 #include <cmath>
7 using namespace std;
8 #ifdef HAVE_EVIO
9 #include <evioFileChannel.hxx>
10 #include <evioUtil.hxx>
11 using namespace evio;
12 #endif
13 
15 
16 //class DApplication;
17 #include <HDGEOMETRY/DGeometry.h>
18 #include <DANA/DApplication.h>
19 
20 //---------------------------------
21 // DMagneticFieldMapPS2DMap (Constructor)
22 //---------------------------------
23 DMagneticFieldMapPS2DMap::DMagneticFieldMapPS2DMap(JApplication *japp, int32_t runnumber, string namepath)
24 {
25  jcalib = japp->GetJCalibration(runnumber);
26  jresman = japp->GetJResourceManager(runnumber);
27  DApplication *dapp = dynamic_cast<DApplication *>(japp);
28  geom = dapp->GetDGeometry(runnumber);
29 
30  JParameterManager *jparms = japp->GetJParameterManager();
31  jparms->SetDefaultParameter("PSBFIELD_MAP", namepath);
32 
33  int Npoints = ReadMap(namepath, runnumber);
34  if(Npoints==0){
35  _DBG_<<"Error getting JCalibration object for magnetic field!"<<endl;
36  japp->Quit();
37  }
38 }
39 
40 
41 //---------------------------------
42 // DMagneticFieldMapPS2DMap (Constructor)
43 //---------------------------------
44 DMagneticFieldMapPS2DMap::DMagneticFieldMapPS2DMap(JCalibration *jcalib, string namepath)
45 {
46  this->jcalib = jcalib;
47  this->geom = NULL;
48  this->jresman = NULL;
49  if(ReadMap(namepath)==0){
50  _DBG_<<"Error getting JCalibration object for magnetic field!"<<endl;
51  exit(1);
52  }
53 }
54 
55 //---------------------------------
56 // ~DMagneticFieldMapPS2DMap (Destructor)
57 //---------------------------------
59 {
60 
61 }
62 
63 //---------------------------------
64 // ReadMap
65 //---------------------------------
66 int DMagneticFieldMapPS2DMap::ReadMap(string namepath, int32_t runnumber, string context)
67 {
68  /// Read the magnetic field map in from the calibration database.
69  /// This will read in the map and figure out the number of grid
70  /// points in each direction (x,y, and z) and the range in each.
71  /// The gradiant of the field is calculated for all but the most
72  /// exterior points and saved to use in later calls to GetField(...).
73 
74  // Read in map from calibration database. This should really be
75  // built into a run-dependent geometry framework, but for now
76  // we do it this way.
77  if(!jcalib){
78  jerr << "ERROR: jcalib pointer is NULL in DMagneticFieldMapPS2DMap::ReadMap() !" << endl;
79  return 0;
80  }
81  if(!jresman){
82  jerr << "ERROR: jresman pointer is NULL in DMagneticFieldMapPS2DMap::ReadMap() !" << endl;
83  return 0;
84  }
85  if(!geom){
86  jerr << "ERROR: geom pointer is NULL in DMagneticFieldMapPS2DMap::ReadMap() !" << endl;
87  return 0;
88  }
89 
90  jout<<endl;
91  jout<<"Reading Pair Spectrometer Magnetic field map from "<<namepath<<" ..."<<endl;
92  vector< vector<float> > Bmap;
93 
94  // Try getting it as a JANA resource.
95  jresman->Get(namepath, Bmap);
96 
97  jout<<Bmap.size()<<" entries found (";
98  // The map should be on a grid with equal spacing in x,y, and z.
99  // Here we want to determine the number of points in each of these
100  // dimensions and the range. Note that all of or maps right now are
101  // 2-dimensional with extent only in x and z.
102  // The easiest way to do this is to use a map<float, int> to make a
103  // histogram of the entries by using the key to hold the extent
104  // so that the number of entries will be equal to the number of
105  // different values.
106  map<float, int> xvals;
107  map<float, int> yvals;
108  map<float, int> zvals;
109  xmin = ymin = zmin = 1.0E6;
110  xmax = ymax = zmax = -1.0E6;
111  for(unsigned int i=0; i<Bmap.size(); i++){
112  vector<float> &a = Bmap[i];
113  float &x = a[0];
114  float &y = a[1];
115  float &z = a[2];
116 
117  xvals[x] = 1;
118  yvals[y] = 1;
119  zvals[z] = 1;
120  if(x<xmin)xmin=x;
121  if(y<ymin)ymin=y;
122  if(z<zmin)zmin=z;
123  if(x>xmax)xmax=x;
124  if(y>ymax)ymax=y;
125  if(z>zmax)zmax=z;
126  }
127  Nx = xvals.size();
128  Ny = yvals.size();
129  Nz = zvals.size();
130  cout<<" Nx="<<Nx;
131  cout<<" Ny="<<Ny;
132  cout<<" Nz="<<Nz;
133  cout<<" ) at 0x"<<hex<<(unsigned long)this<<dec<<endl;
134 
135  // Create 3D vector so we can index the values by [x][y][z]
136  vector<DBfieldPoint_t> zvec(Nz);
137  vector< vector<DBfieldPoint_t> > yvec;
138  for(int i=0; i<Ny; i++)yvec.push_back(zvec);
139  for(int i=0; i<Nx; i++)Btable.push_back(yvec);
140 
141  // Distance between map points for r and z
142  dx = (xmax-xmin)/(double)(Nx-1);
143  dy = (ymax-ymin)/(double)((Ny < 2)? 1 : Ny-1);
144  dz = (zmax-zmin)/(double)(Nz-1);
145 
146  one_over_dx=1./dx;
147  one_over_dz=1./dz;
148 
149  // Copy values into Btable
150  for(unsigned int i=0; i<Bmap.size(); i++){
151  vector<float> &a = Bmap[i];
152  int xindex = (int)floor((a[0]-xmin+dx/2.0)/dx); // the +dx/2.0 guarantees against round-off errors
153  int yindex = (int)(Ny<2 ? 0:floor((a[1]-ymin+dy/2.0)/dy));
154  int zindex = (int)floor((a[2]-zmin+dz/2.0)/dz);
155  DBfieldPoint_t *b = &Btable[xindex][yindex][zindex];
156  b->x = a[0];
157  b->y = a[1];
158  b->z = a[2];
159  // convert to T
160  b->Bx = a[3]/10000.;
161  b->By = a[4]/10000.;
162  b->Bz = a[5]/10000.;
163  //b->Bx = a[3];
164  //b->By = a[4];
165  //b->Bz = a[5];
166  }
167 
168  // Calculate gradient at every point in the map.
169  // The derivatives are calculated wrt to the *index* of the
170  // field map, not the physical unit. This is because
171  // the fractions used in interpolation are fractions
172  // between indices.
173  for(int index_x=0; index_x<Nx; index_x++){
174  int index_x0 = index_x - (index_x>0 ? 1:0);
175  int index_x1 = index_x + (index_x<Nx-1 ? 1:0);
176  for(int index_z=1; index_z<Nz-1; index_z++){
177  int index_z0 = index_z - (index_z>0 ? 1:0);
178  int index_z1 = index_z + (index_z<Nz-1 ? 1:0);
179 
180  // Right now, all of our maps are 2-D with no y extent
181  // so we comment out the following for-loop and hardwire
182  // the index_y values to zero
183  //for(int index_y=1; index_y<Ny-1; index_y++){
184  // int index_y0 = index_y-1;
185  // int index_y1 = index_y+1;
186  if(Ny>1){
187  _DBG_<<"Field map appears to be 3 dimensional. Code is currently"<<endl;
188  _DBG_<<"unable to handle this. Exiting..."<<endl;
189  exit(-1);
190  }else{
191  int index_y=0;
192  /*
193  int index_y0, index_y1;
194  index_y=index_y0=index_y1=0;
195  */
196  double d_index_x=double(index_x1-index_x0);
197  double d_index_z=double(index_z1-index_z0);
198 
199  DBfieldPoint_t *Bx0 = &Btable[index_x0][index_y][index_z];
200  DBfieldPoint_t *Bx1 = &Btable[index_x1][index_y][index_z];
201  //DBfieldPoint_t *By0 = &Btable[index_x][index_y0][index_z];
202  //DBfieldPoint_t *By1 = &Btable[index_x][index_y1][index_z];
203  DBfieldPoint_t *Bz0 = &Btable[index_x][index_y][index_z0];
204  DBfieldPoint_t *Bz1 = &Btable[index_x][index_y][index_z1];
205 
206  DBfieldPoint_t *g = &Btable[index_x][index_y][index_z];
207  g->dBxdx = (Bx1->Bx - Bx0->Bx)/d_index_x;
208  //g->dBxdy = (By1->Bx - By0->Bx)/(double)(index_y1-index_y0);
209  g->dBxdz = (Bz1->Bx - Bz0->Bx)/d_index_z;
210 
211  g->dBydx = (Bx1->By - Bx0->By)/d_index_x;
212  //g->dBydy = (By1->By - By0->By)/(double)(index_y1-index_y0);
213  g->dBydz = (Bz1->By - Bz0->By)/d_index_z;
214 
215  g->dBzdx = (Bx1->Bz - Bx0->Bz)/d_index_x;
216  //g->dBzdy = (By1->Bz - By0->Bz)/(double)(index_y1-index_y0);
217  g->dBzdz = (Bz1->Bz - Bz0->Bz)/d_index_z;
218 
219  DBfieldPoint_t *B11 = &Btable[index_x1][index_y][index_z1];
220  DBfieldPoint_t *B01 = &Btable[index_x0][index_y][index_z1];
221  DBfieldPoint_t *B10 = &Btable[index_x1][index_y][index_z0];
222  DBfieldPoint_t *B00 = &Btable[index_x0][index_y][index_z0];
223 
224  g->dBxdxdz=(B11->Bx - B01->Bx - B10->Bx + B00->Bx)/d_index_x/d_index_z;
225  g->dBzdxdz=(B11->Bz - B01->Bz - B10->Bz + B00->Bz)/d_index_x/d_index_z;
226  }
227  }
228  }
229 
230  // fix this
231  vector<double>pair_spect_origin;
232  geom->Get("//posXYZ[@volume='PairSpectrometer']/@X_Y_Z",pair_spect_origin);
233  vector<double>pair_magnet_origin;
234  geom->Get("//posXYZ[@volume='PairMagnet']/@X_Y_Z",pair_magnet_origin);
235 
236  z_shift = pair_spect_origin[2]-pair_magnet_origin[2];
237 
238  //cout << " PS volume origin = *" << pair_spect_origin[0] << ", " << pair_spect_origin[1] << ", " << pair_spect_origin[2] << ")" << endl;
239  //cout << " PS magnet volume origin = *" << pair_magnet_origin[0] << ", " << pair_magnet_origin[1] << ", " << pair_magnet_origin[2] << ")" << endl;
240 
241  //jout << "Map Z shift = " << z_shift << endl;
242 
243  return Bmap.size();
244 }
245 
246 // Use bicubic interpolation to find the field at the point (x,y).
247 //See Numerical Recipes in C (2nd ed.), pp.125-127.
248 void DMagneticFieldMapPS2DMap::GetFieldBicubic(double x,double y,double z,
249  double &Bx_,double &By_,
250  double &Bz_) const{
251  // transform from hall coordinates to map coordinates
252  z -= z_shift;
253 
254  // table of weight factors for bicubic interpolation
255  static const int wt[16][16]=
256  { {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
257  {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
258  {-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1, 0, 0, 0, 0},
259  {2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0},
260  {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
261  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
262  {0, 0, 0, 0,-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1},
263  {0, 0, 0, 0, 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1},
264  {-3,3, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
265  {0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0},
266  {9,-9, 9,-9, 6, 3,-3,-6, 6,-6,-3, 3, 4, 2, 1, 2},
267  {-6,6,-6, 6,-4,-2, 2, 4,-3, 3, 3,-3,-2,-1,-1,-2},
268  {2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
269  {0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0},
270  {-6, 6,-6,6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1},
271  {4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1}};
272  double temp[16];
273  double cl[16],coeff[4][4];
274 
275  // Interpolate using coarse grid (fine grid not implemented yet)
276  // Get closest indices for this point
277  int index_x = (int)floor((x-xmin)*one_over_dx + 0.5);
278  if(index_x<0 || index_x>=Nx)return;
279  //if (index_x>=Nx) return;
280  else if (index_x<0) index_x=0;
281  int index_z = (int)floor((z-zmin)*one_over_dz + 0.5);
282  if(index_z<0 || index_z>=Nz)return;
283 
284  int index_y = 0;
285 
286  int i, j, k, m=0;
287  int index_x1 = index_x + (index_x<Nx-1 ? 1:0);
288  int index_z1 = index_z + (index_z<Nz-1 ? 1:0);
289 
290  // Pointers to magnetic field structure
291  const DBfieldPoint_t *B00 = &Btable[index_x][index_y][index_z];
292  const DBfieldPoint_t *B01 = &Btable[index_x][index_y][index_z1];
293  const DBfieldPoint_t *B11 = &Btable[index_x1][index_y][index_z1];
294  const DBfieldPoint_t *B10 = &Btable[index_x1][index_y][index_z];
295 
296  // First compute the interpolation for Bx
297  temp[0]=B00->Bx;
298  temp[1]=B01->Bx;
299  temp[2]=B11->Bx;
300  temp[3]=B10->Bx;
301 
302  temp[8]=B00->dBxdx;
303  temp[9]=B01->dBxdx;
304  temp[10]=B11->dBxdx;
305  temp[11]=B10->dBxdx;
306 
307  temp[4]=B00->dBxdz;
308  temp[5]=B01->dBxdz;
309  temp[6]=B11->dBxdz;
310  temp[7]=B10->dBxdz;
311 
312  temp[12]=B00->dBxdxdz;
313  temp[13]=B01->dBxdxdz;
314  temp[14]=B11->dBxdxdz;
315  temp[15]=B10->dBxdxdz;
316 
317  for (i=0;i<16;i++){
318  double tmp2=0.0;
319  for (k=0;k<16;k++) tmp2+=wt[i][k]*temp[k];
320  cl[i]=tmp2;
321  }
322  for (i=0;i<4;i++)
323  for (j=0;j<4;j++) coeff[i][j]=cl[m++];
324 
325  double t=(z - B00->z)*one_over_dz;
326  double u=(x - B00->x)*one_over_dx;
327  for (i=3;i>=0;i--){
328  Bx_=t*Bx_+((coeff[i][3]*u+coeff[i][2])*u+coeff[i][1])*u+coeff[i][0];
329  }
330 
331  // Next compute the interpolation for By
332  temp[0]=B00->By;
333  temp[1]=B01->By;
334  temp[2]=B11->By;
335  temp[3]=B10->By;
336 
337  temp[8]=B00->dBydx;
338  temp[9]=B01->dBydx;
339  temp[10]=B11->dBydx;
340  temp[11]=B10->dBydx;
341 
342  temp[4]=B00->dBydy;
343  temp[5]=B01->dBydy;
344  temp[6]=B11->dBydy;
345  temp[7]=B10->dBydy;
346 
347  temp[12]=B00->dBydxdy;
348  temp[13]=B01->dBydxdy;
349  temp[14]=B11->dBydxdy;
350  temp[15]=B10->dBydxdy;
351 
352  for (i=0;i<16;i++){
353  double tmp2=0.0;
354  for (k=0;k<16;k++) tmp2+=wt[i][k]*temp[k];
355  cl[i]=tmp2;
356  }
357  m=0;
358  for (i=0;i<4;i++)
359  for (j=0;j<4;j++) coeff[i][j]=cl[m++];
360 
361  By_=0.;
362  for (i=3;i>=0;i--){
363  By_=t*By_+((coeff[i][3]*u+coeff[i][2])*u+coeff[i][1])*u+coeff[i][0];
364  }
365 
366  // Next compute the interpolation for Bz
367  temp[0]=B00->Bz;
368  temp[1]=B01->Bz;
369  temp[2]=B11->Bz;
370  temp[3]=B10->Bz;
371 
372  temp[8]=B00->dBzdx;
373  temp[9]=B01->dBzdx;
374  temp[10]=B11->dBzdx;
375  temp[11]=B10->dBzdx;
376 
377  temp[4]=B00->dBzdz;
378  temp[5]=B01->dBzdz;
379  temp[6]=B11->dBzdz;
380  temp[7]=B10->dBzdz;
381 
382  temp[12]=B00->dBzdxdz;
383  temp[13]=B01->dBzdxdz;
384  temp[14]=B11->dBzdxdz;
385  temp[15]=B10->dBzdxdz;
386 
387  for (i=0;i<16;i++){
388  double tmp2=0.0;
389  for (k=0;k<16;k++) tmp2+=wt[i][k]*temp[k];
390  cl[i]=tmp2;
391  }
392  m=0;
393  for (i=0;i<4;i++)
394  for (j=0;j<4;j++) coeff[i][j]=cl[m++];
395 
396  Bz_=0.;
397  for (i=3;i>=0;i--){
398  Bz_=t*Bz_+((coeff[i][3]*u+coeff[i][2])*u+coeff[i][1])*u+coeff[i][0];
399  }
400 
401 }
402 
403 
404 // Find the field and field gradient at the point (x,y,z).
406  double &Bx_,double &By_,
407  double &Bz_,
408  double &dBxdx_,
409  double &dBxdy_,
410  double &dBxdz_,
411  double &dBydx_,
412  double &dBydy_,
413  double &dBydz_,
414  double &dBzdx_,
415  double &dBzdy_,
416  double &dBzdz_) const{
417  // transform from hall coordinates to map coordinates
418  z -= z_shift;
419 
420  // radial distance
421  double r = sqrt(x*x + y*y);
422  if (r>xmax || z>zmax || z<zmin){
423  Bx_=0.0,By_=0.0,Bz_=0.0;
424  dBxdx_=0.0,dBxdy_=0.0,dBxdz_=0.0;
425  dBydx_=0.0,dBydy_=0.0,dBydz_=0.0;
426  dBzdx_=0.0,dBzdy_=0.0,dBzdz_=0.0;
427  return;
428  }
429 
430  // Initialize z-component
431  Bz_=0.;
432 
433  // Get closest indices for this point
434  int index_x = static_cast<int>((x-xmin)*one_over_dx);
435  int index_z = static_cast<int>((z-zmin)*one_over_dz);
436  int index_y = 0;
437 
438  if(index_x>=0 && index_x<Nx && index_z>=0 && index_z<Nz){
439  const DBfieldPoint_t *B = &Btable[index_x][index_y][index_z];
440 
441  // Fractional distance between map points.
442  double ux = (x - B->x)*one_over_dx;
443  double uz = (z - B->z)*one_over_dz;
444 
445  // Use gradient to project grid point to requested position
446  Bx_ = B->Bx+B->dBxdx*ux+B->dBxdz*uz;
447  By_ = B->By;
448  Bz_ = B->Bz+B->dBzdx*ux+B->dBzdz*uz;
449  dBxdx_=B->dBxdx;
450  dBxdy_=B->dBxdy;
451  dBxdz_=B->dBxdz;
452  dBydx_=B->dBydx;
453  dBydy_=B->dBydy;
454  dBydz_=B->dBydz;
455  dBzdx_=B->dBzdx;
456  dBzdy_=B->dBzdy;
457  dBzdz_=B->dBzdz;
458 
459  dBydx_= dBxdy_;
460 
461  /*
462  printf("Grad %f %f %f %f %f %f %f %f %f\n",dBxdx_,dBxdy_,dBxdz_,
463  dBydx_,dBydy_,dBydz_,dBzdx_,dBzdy_,dBzdz_);
464  */
465  }
466 }
467 
468 
469 
470 //-------------
471 // GetFieldGradient
472 //-------------
473 void DMagneticFieldMapPS2DMap::GetFieldGradient(double x, double y, double z,
474  double &dBxdx, double &dBxdy,
475  double &dBxdz,
476  double &dBydx, double &dBydy,
477  double &dBydz,
478  double &dBzdx, double &dBzdy,
479  double &dBzdz) const{
480 
481  // transform from hall coordinates to map coordinates
482  z -= z_shift;
483 
484  // Get closest indices for this point
485  int index_x = (int)floor((x-xmin)*one_over_dx + 0.5);
486  if(index_x<0 || index_x>=Nx)return;
487  //if (index_x>=Nx) return;
488  else if (index_x<0) index_x=0;
489  int index_z = (int)floor((z-zmin)*one_over_dz + 0.5);
490  if(index_z<0 || index_z>=Nz)return;
491 
492  int index_y = 0;
493 
494  const DBfieldPoint_t *B = &Btable[index_x][index_y][index_z];
495 
496  dBxdx = B->dBxdx;
497  dBxdy = B->dBxdx;
498  dBxdz = B->dBxdz;
499  dBydx = B->dBydx;
500  dBydy = B->dBydx;
501  dBydz = B->dBydz;
502  dBzdx = B->dBzdx;
503  dBzdy = B->dBzdx;
504  dBzdz = B->dBzdz;
505 }
506 
507 
508 
509 //---------------------------------
510 // GetField
511 //---------------------------------
512 void DMagneticFieldMapPS2DMap::GetField(double x, double y, double z, double &Bx, double &By, double &Bz, int method) const
513 {
514  //cout << "in GetField()..." << endl;
515 
516  /// This calculates the magnetic field at an arbitrary point
517  /// in space using the field map read from the calibaration
518  /// database. It interpolates between grid points using the
519  /// gradient values calculated in ReadMap (called from the
520  /// constructor).
521 
522  //cout << " hall coords: x = " << x << " z = " << z << endl;
523 
524  // transform from hall coordinates to map coordinates
525  z -= z_shift;
526 
527  //cout << " local coords: x = " << x << " z = " << z << endl;
528 
529  Bx = By = Bz = 0.0;
530 
531  if(Ny>1){
532  _DBG_<<"Field map appears to be 3 dimensional. Code is currently"<<endl;
533  _DBG_<<"unable to handle this. Assuming y=0."<<endl;
534  }
535 
536  if (x<xmin || x>xmax || z>zmax || z<zmin){
537  return;
538  }
539 
540  // Interpolate using coarse grid
541  // Get closest indices for this point
542  int index_x = static_cast<int>((x-xmin)*one_over_dx);
543  if (index_x<0 || index_x>=Nx) return;
544 
545  int index_z = static_cast<int>((z-zmin)*one_over_dz);
546  if(index_z<0 || index_z>=Nz)return;
547 
548  int index_y = 0;
549 
550  //cout << " indexes: x = " << index_x << " z = " << index_z << endl;
551 
552  const DBfieldPoint_t *B = &Btable[index_x][index_y][index_z];
553 
554  //cout << " mag field: Bx = " << B->x << " By = " << B->y << " Bz = " << B->z << endl;
555 
556  // Fractional distance between map points.
557  double ux = (x - B->x)*one_over_dx;
558  double uz = (z - B->z)*one_over_dz;
559 
560  // Use gradient to project grid point to requested position
561  Bx = B->Bx+B->dBxdx*ux+B->dBxdz*uz;
562  By = B->By;
563  Bz = B->Bz+B->dBzdx*ux+B->dBzdz*uz;
564 }
565 
566 //---------------------------------
567 // GetField
568 //---------------------------------
570 {
571  /// This calculates the magnetic field at an arbitrary point
572  /// in space using the field map read from the calibaration
573  /// database. It interpolates between grid points using the
574  /// gradient values calculated in ReadMap (called from the
575  /// constructor).
576 
577 
578  double Bz=0.,Bx=0.0,By=0.0;
579 
580  if(Ny>1){
581  _DBG_<<"Field map appears to be 3 dimensional. Code is currently"<<endl;
582  _DBG_<<"unable to handle this. Assuming y=0."<<endl;
583  }
584 
585  // radial position and angles
586  double x = pos.x();
587  //double y = pos.y();
588  double z = pos.z();
589  // transform from hall coordinates to map coordinates
590  z -= z_shift;
591  if (x<xmin || x>xmax || z>zmax || z<zmin){
592  return;
593  }
594 
595  // Interpolate using coarse grid
596  // Get closest indices for this point
597  int index_x = static_cast<int>((x-xmin)*one_over_dx);
598  if (index_x<0 || index_x>=Nx) return;
599 
600  int index_z = static_cast<int>((z-zmin)*one_over_dz);
601  if(index_z<0 || index_z>=Nz)return;
602 
603  int index_y = 0;
604 
605  const DBfieldPoint_t *B = &Btable[index_x][index_y][index_z];
606 
607  // Fractional distance between map points.
608  double ux = (x - B->x)*one_over_dx;
609  double uz = (z - B->z)*one_over_dz;
610 
611  // Use gradient to project grid point to requested position
612  Bx = B->Bx+B->dBxdx*ux+B->dBxdz*uz;
613  By = B->By;
614  Bz = B->Bz+B->dBzdx*ux+B->dBzdz*uz;
615 
616  // Set result
617  Bout.SetXYZ(Bx,By,Bz);
618 }
619 
620 
621 
622 
int ReadMap(string namepath, int32_t runnumber=1, string context="")
DApplication * dapp
TVector3 DVector3
Definition: DVector3.h:14
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
#define y
DMagneticFieldMapPS2DMap(JApplication *japp, int32_t runnumber=1, string namepath="Magnets/PairSpectrometer/PS_1.8T_20150513_test")
double zmax
JApplication * japp
int Nz
Definition: bfield2root.cc:27
DGeometry * GetDGeometry(unsigned int run_number)
#define _DBG_
Definition: HDEVIO.h:12
void GetFieldGradient(double x, double y, double z, double &dBxdx, double &dBxdy, double &dBxdz, double &dBydx, double &dBydy, double &dBydz, double &dBzdx, double &dBzdy, double &dBzdz) const
double sqrt(double)
void GetFieldAndGradient(double x, double y, double z, double &Bx, double &By, double &Bz, double &dBxdx, double &dBxdy, double &dBxdz, double &dBydx, double &dBydy, double &dBydz, double &dBzdx, double &dBzdy, double &dBzdz) const
Double_t ymin
Definition: bcal_hist_eff.C:89
void GetField(const DVector3 &pos, DVector3 &Bout) const
Double_t ymax
Definition: bcal_hist_eff.C:91
void GetFieldBicubic(double x, double y, double z, double &Bx, double &By, double &Bz) const
double zmin
TH2F * temp
union @6::@8 u