9 #include <evioFileChannel.hxx>
10 #include <evioUtil.hxx>
25 jcalib = japp->GetJCalibration(runnumber);
26 jresman = japp->GetJResourceManager(runnumber);
30 JParameterManager *jparms = japp->GetJParameterManager();
31 jparms->SetDefaultParameter(
"PSBFIELD_MAP", namepath);
33 int Npoints = ReadMap(namepath, runnumber);
35 _DBG_<<
"Error getting JCalibration object for magnetic field!"<<endl;
46 this->jcalib = jcalib;
49 if(ReadMap(namepath)==0){
50 _DBG_<<
"Error getting JCalibration object for magnetic field!"<<endl;
78 jerr <<
"ERROR: jcalib pointer is NULL in DMagneticFieldMapPS2DMap::ReadMap() !" << endl;
82 jerr <<
"ERROR: jresman pointer is NULL in DMagneticFieldMapPS2DMap::ReadMap() !" << endl;
86 jerr <<
"ERROR: geom pointer is NULL in DMagneticFieldMapPS2DMap::ReadMap() !" << endl;
91 jout<<
"Reading Pair Spectrometer Magnetic field map from "<<namepath<<
" ..."<<endl;
92 vector< vector<float> > Bmap;
95 jresman->Get(namepath, Bmap);
97 jout<<Bmap.size()<<
" entries found (";
106 map<float, int> xvals;
107 map<float, int> yvals;
108 map<float, int> zvals;
111 for(
unsigned int i=0; i<Bmap.size(); i++){
112 vector<float> &a = Bmap[i];
133 cout<<
" ) at 0x"<<hex<<(
unsigned long)
this<<dec<<endl;
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);
142 dx = (xmax-xmin)/(
double)(Nx-1);
143 dy = (
ymax-
ymin)/(
double)((Ny < 2)? 1 : Ny-1);
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);
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);
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);
187 _DBG_<<
"Field map appears to be 3 dimensional. Code is currently"<<endl;
188 _DBG_<<
"unable to handle this. Exiting..."<<endl;
196 double d_index_x=double(index_x1-index_x0);
197 double d_index_z=double(index_z1-index_z0);
207 g->
dBxdx = (Bx1->
Bx - Bx0->
Bx)/d_index_x;
209 g->
dBxdz = (Bz1->
Bx - Bz0->
Bx)/d_index_z;
211 g->
dBydx = (Bx1->
By - Bx0->
By)/d_index_x;
213 g->
dBydz = (Bz1->
By - Bz0->
By)/d_index_z;
215 g->
dBzdx = (Bx1->
Bz - Bx0->
Bz)/d_index_x;
217 g->
dBzdz = (Bz1->
Bz - Bz0->
Bz)/d_index_z;
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;
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);
236 z_shift = pair_spect_origin[2]-pair_magnet_origin[2];
249 double &Bx_,
double &By_,
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}};
273 double cl[16],coeff[4][4];
277 int index_x = (int)floor((x-xmin)*one_over_dx + 0.5);
278 if(index_x<0 || 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;
287 int index_x1 = index_x + (index_x<Nx-1 ? 1:0);
288 int index_z1 = index_z + (index_z<
Nz-1 ? 1:0);
312 temp[12]=B00->dBxdxdz;
319 for (k=0;k<16;k++) tmp2+=wt[i][k]*temp[k];
323 for (j=0;j<4;j++) coeff[i][j]=cl[m++];
325 double t=(z - B00->z)*one_over_dz;
326 double u=(x - B00->x)*one_over_dx;
328 Bx_=t*Bx_+((coeff[i][3]*u+coeff[i][2])*u+coeff[i][1])*u+coeff[i][0];
347 temp[12]=B00->dBydxdy;
354 for (k=0;k<16;k++) tmp2+=wt[i][k]*temp[k];
359 for (j=0;j<4;j++) coeff[i][j]=cl[m++];
363 By_=t*By_+((coeff[i][3]*u+coeff[i][2])*u+coeff[i][1])*u+coeff[i][0];
382 temp[12]=B00->dBzdxdz;
389 for (k=0;k<16;k++) tmp2+=wt[i][k]*temp[k];
394 for (j=0;j<4;j++) coeff[i][j]=cl[m++];
398 Bz_=t*Bz_+((coeff[i][3]*u+coeff[i][2])*u+coeff[i][1])*u+coeff[i][0];
406 double &Bx_,
double &By_,
416 double &dBzdz_)
const{
421 double r =
sqrt(x*x + y*y);
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;
434 int index_x =
static_cast<int>((x-xmin)*one_over_dx);
435 int index_z =
static_cast<int>((z-
zmin)*one_over_dz);
438 if(index_x>=0 && index_x<Nx && index_z>=0 && index_z<
Nz){
442 double ux = (x - B->
x)*one_over_dx;
443 double uz = (z - B->
z)*one_over_dz;
474 double &dBxdx,
double &dBxdy,
476 double &dBydx,
double &dBydy,
478 double &dBzdx,
double &dBzdy,
479 double &dBzdz)
const{
485 int index_x = (int)floor((x-xmin)*one_over_dx + 0.5);
486 if(index_x<0 || 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;
532 _DBG_<<
"Field map appears to be 3 dimensional. Code is currently"<<endl;
533 _DBG_<<
"unable to handle this. Assuming y=0."<<endl;
536 if (x<xmin || x>xmax || z>
zmax || z<
zmin){
542 int index_x =
static_cast<int>((x-xmin)*one_over_dx);
543 if (index_x<0 || index_x>=Nx)
return;
545 int index_z =
static_cast<int>((z-
zmin)*one_over_dz);
546 if(index_z<0 || index_z>=
Nz)
return;
557 double ux = (x - B->
x)*one_over_dx;
558 double uz = (z - B->
z)*one_over_dz;
578 double Bz=0.,Bx=0.0,By=0.0;
581 _DBG_<<
"Field map appears to be 3 dimensional. Code is currently"<<endl;
582 _DBG_<<
"unable to handle this. Assuming y=0."<<endl;
591 if (x<xmin || x>xmax || z>
zmax || z<
zmin){
597 int index_x =
static_cast<int>((x-xmin)*one_over_dx);
598 if (index_x<0 || index_x>=Nx)
return;
600 int index_z =
static_cast<int>((z-
zmin)*one_over_dz);
601 if(index_z<0 || index_z>=
Nz)
return;
608 double ux = (x - B->
x)*one_over_dx;
609 double uz = (z - B->
z)*one_over_dz;
617 Bout.SetXYZ(Bx,By,Bz);
int ReadMap(string namepath, int32_t runnumber=1, string context="")
DMagneticFieldMapPS2DMap(JApplication *japp, int32_t runnumber=1, string namepath="Magnets/PairSpectrometer/PS_1.8T_20150513_test")
DGeometry * GetDGeometry(unsigned int run_number)
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
virtual ~DMagneticFieldMapPS2DMap()
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
void GetField(const DVector3 &pos, DVector3 &Bout) const
void GetFieldBicubic(double x, double y, double z, double &Bx, double &By, double &Bz) const