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