11 #include <JANA/JParameterManager.h>
20 #define T2(x_2) (2*x_2-1)
21 #define T3(x_1,x_3) (4*x_3-3*x_1)
22 #define T4(x_2,x_4) (8*x_4-8*x_2+1)
23 #define T5(x_5,x_3,x_1) (16*x_5-20*x_3+5*x_1)
24 #define T6(x_6,x_4,x_2) (32*x_6-48*x_4+18*x_2-1)
25 #define T7(x_7,x_5,x_3,x_1) (64*x_7-112*x_5+56*x_3-7*x_1)
26 #define T8(x_8,x_6,x_4,x_2) (128*x_8-256*x_6+160*x_4-32*x_2+1)
27 #define T9(x_9,x_7,x_5,x_3,x_1) (256*x_9-576*x_7+432*x_5-120*x_3+9*x_1)
32 #define dT2dx(x) (2*2*pow(x,2-1))
33 #define dT3dx(x) (3*4*pow(x,3-1)-3)
34 #define dT4dx(x) (4*8*pow(x,4-1)-2*8*pow(x,2-1))
35 #define dT5dx(x) (5*16*pow(x,5-1)-3*20*pow(x,3-1)+5)
36 #define dT6dx(x) (6*32*pow(x,6-1)-4*48*pow(x,4-1)+2*18*pow(x,2-1))
37 #define dT7dx(x) (7*64*pow(x,7-1)-5*112*pow(x,5-1)+3*56*pow(x,3-1)-7)
38 #define dT8dx(x) (8*128*pow(x,8-1)-6*256*pow(x,6-1)+4*160*pow(x,4-1)-2*32*pow(x,2-1))
39 #define dT9dx(x) (9*256*pow(x,9-1)-7*576*pow(x,7-1)+5*432*pow(x,5-1)-3*120*pow(x,3-1)+9)
46 int32_t runnumber = 1;
47 jcalib = japp->GetJCalibration(runnumber);
49 JParameterManager *jparms = japp->GetJParameterManager();
50 jparms->SetDefaultParameter(
"BFIELD_MAP", namepath);
52 Init(jcalib, namepath);
60 Init(jcalib, namepath);
76 this->jcalib = jcalib;
79 vector<map<string, string> > sections_map;
80 jcalib->Get(namepath, sections_map);
83 for(
unsigned int k=0; k<sections_map.size(); k++){
84 map<string, string> &vals = sections_map[k];
88 section.
Bi = vals[
"Bi"];
89 section.
section = atoi(vals[
"sec"].c_str());
90 section.
zmin = atof(vals[
"zmin"].c_str());
91 section.
zmax = atof(vals[
"zmax"].c_str());
92 section.
zmid = atof(vals[
"zmid"].c_str());
93 section.
znorm = atof(vals[
"znorm"].c_str());
94 section.
rmid = atof(vals[
"rmid"].c_str());
95 section.
rnorm = atof(vals[
"rnorm"].c_str());
99 section.
order1 = section.
pp.size() - 1;
100 section.
order2 = section.
pp[0].size() - 1;
102 cout<<
"---- parameterized B-field section ---------------"<<endl;
103 cout<<
"namepath = "<<section.
namepath<<endl;
104 cout<<
" Bi = "<<section.
Bi<<endl;
105 cout<<
" section = "<<section.
section<<endl;
106 cout<<
" zmin = "<<section.
zmin<<endl;
107 cout<<
" zmax = "<<section.
zmax<<endl;
108 cout<<
" zmid = "<<section.
zmid<<endl;
109 cout<<
" znorm = "<<section.
znorm<<endl;
110 cout<<
" rmid = "<<section.
rmid<<endl;
111 cout<<
" rnorm = "<<section.
rnorm<<endl;
112 cout<<
" order1 = "<<section.
order1<<endl;
113 cout<<
" order2 = "<<section.
order1<<endl;
124 unsigned int &order1 = section.
order1;
125 DMatrix C1(order1+1, order1+1);
128 case 9: C1(9,1)= +9; C1(9,3)=-120; C1(9,5)=+432; C1(9,7)=-576; C1(9,9)=+256;
129 case 8: C1(8,0)= +1; C1(8,2)= -32; C1(8,4)=+160; C1(8,6)=-256; C1(8,8)=+128;
130 case 7: C1(7,1)= -7; C1(7,3)= +56; C1(7,5)=-112; C1(7,7)= +64;
131 case 6: C1(6,0)= -1; C1(6,2)= +18; C1(6,4)= -48; C1(6,6)= +32;
132 case 5: C1(5,1)= +5; C1(5,3)= -20; C1(5,5)= +16;
133 case 4: C1(4,0)= +1; C1(4,2)= -8; C1(4,4)= +8;
134 case 3: C1(3,1)= -3; C1(3,3)= +4;
135 case 2: C1(2,0)= -1; C1(2,2)= +2;
139 TMatrixD C1_t(TMatrixD::kTransposed, C1);
142 unsigned int &order2 = section.
order2;
143 DMatrix C2(order2+1, order2+1);
146 case 9: C2(9,1)= +9; C2(9,3)=-120; C2(9,5)=+432; C2(9,7)=-576; C2(9,9)=+256;
147 case 8: C2(8,0)= +1; C2(8,2)= -32; C2(8,4)=+160; C2(8,6)=-256; C2(8,8)=+128;
148 case 7: C2(7,1)= -7; C2(7,3)= +56; C2(7,5)=-112; C2(7,7)= +64;
149 case 6: C2(6,0)= -1; C2(6,2)= +18; C2(6,4)= -48; C2(6,6)= +32;
150 case 5: C2(5,1)= +5; C2(5,3)= -20; C2(5,5)= +16;
151 case 4: C2(4,0)= +1; C2(4,2)= -8; C2(4,4)= +8;
152 case 3: C2(3,1)= -3; C2(3,3)= +4;
153 case 2: C2(2,0)= -1; C2(2,2)= +2;
160 for(
unsigned int i=0; i<=order1; i++){
161 for(
unsigned int j=0; j<=order2; j++){
162 P(j, i) = section.
pp[i][j];
167 section.
Q.ResizeTo(P);
168 section.
Q = C1_t*P*C2;
171 section.
cc = section.
pp;
172 for(
unsigned int i=0; i<=order1; i++){
173 for(
unsigned int j=0; j<=order2; j++){
174 section.
cc[i][j] = section.
Q(j,i);
179 string type = vals[
"Bi"];
180 if(type==
"Bx")sections_Bx.push_back(section);
181 if(type==
"Bz")sections_Bz.push_back(section);
188 GetField(pos.x(),pos.y(),pos.z(),Bx,By,Bz);
189 Bout.SetXYZ(Bx,By,Bz);
201 double r =
sqrt(x*x + y*y);
204 for(
unsigned int i=0; i<sections_Bx.size(); i++){
205 if(sections_Bx[i].IsInRange(z)){
206 Br = sections_Bx[i].Eval(r,z);
212 for(
unsigned int i=0; i<sections_Bz.size(); i++){
213 if(sections_Bz[i].IsInRange(z)){
214 Bz = sections_Bz[i].Eval(r,z);
231 double r =
sqrt(x*x + y*y);
234 for(
unsigned int i=0; i<sections_Bz.size(); i++){
235 if(sections_Bz[i].IsInRange(z)){
236 return (sections_Bz[i].Eval(r,z));
255 double r_prime = (r-rmid)/rnorm;
256 double z_prime = (z-zmid)/znorm;
260 double &r_1 = r_prime;
261 double r_2 = r_1*r_1;
262 double r_3 = r_2*r_1;
263 double r_4 = r_3*r_1;
264 double r_5 = r_4*r_1;
265 double r_6 = r_5*r_1;
266 double r_7 = r_6*r_1;
267 double r_8 = r_7*r_1;
268 double r_9 = r_8*r_1;
271 double &z_1 = z_prime;
272 double z_2 = z_1*z_1;
273 double z_3 = z_2*z_1;
274 double z_4 = z_3*z_1;
275 double z_5 = z_4*z_1;
276 double z_6 = z_5*z_1;
277 double z_7 = z_6*z_1;
278 double z_8 = z_7*z_1;
279 double z_9 = z_8*z_1;
283 for(
unsigned int i=0; i<order1; i++){
288 const double *my_pp = &pp[i][0];
289 double p = my_pp[0]*
T0();
291 case 9: p += my_pp[9]*
T9(r_9,r_7,r_5,r_3,r_1);
292 case 8: p += my_pp[8]*
T8(r_8,r_6,r_4,r_2);
293 case 7: p += my_pp[7]*
T7(r_7,r_5,r_3,r_1);
294 case 6: p += my_pp[6]*
T6(r_6,r_4,r_2);
295 case 5: p += my_pp[5]*
T5(r_5,r_3,r_1);
296 case 4: p += my_pp[4]*
T4(r_4,r_2);
297 case 3: p += my_pp[3]*
T3(r_3,r_1);
298 case 2: p += my_pp[2]*
T2(r_2);
299 case 1: p += my_pp[1]*
T1(r_1);
304 case 0: Bi += p*
T0();
break;
305 case 1: Bi += p*
T1(z_1);
break;
306 case 2: Bi += p*
T2(z_2);
break;
307 case 3: Bi += p*
T3(z_3,z_1);
break;
308 case 4: Bi += p*
T4(z_4,z_2);
break;
309 case 5: Bi += p*
T5(z_5,z_3,z_1);
break;
310 case 6: Bi += p*
T6(z_6,z_4,z_2);
break;
311 case 7: Bi += p*
T7(z_7,z_5,z_3,z_1);
break;
312 case 8: Bi += p*
T8(z_8,z_6,z_4,z_2);
break;
313 case 9: Bi += p*
T9(z_9,z_7,z_5,z_3,z_1);
break;
319 for(
unsigned int i=1; i<=order2; i++){
320 const double *col_ptr = &cc[i][0];
321 double v = *col_ptr++;
323 for(
unsigned int j=1; j<=order1; j++, col_ptr++){
335 double rr= R(0,0) = 1;
336 for(
unsigned int i=1; i<=order2; i++)R(i,0)=(rr*=r_prime);
340 double zz= Z(0,0) = 1;
341 for(
unsigned int i=1; i<=order1; i++)Z(0,i)=(zz*=z_prime);
357 double &dBxdx,
double &dBxdy,
359 double &dBydx,
double &dBydy,
361 double &dBzdx,
double &dBzdy,
375 GetField(x,y,z,Bx, By, Bz);
396 GetField(x,y,z,Bx, By, Bz);
397 GetFieldGradient(x, y, z, dBxdx, dBxdy, dBxdz, dBydx, dBydy, dBydz, dBzdx, dBzdy, dBzdz);
virtual 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
virtual void GetFieldBicubic(double x, double y, double z, double &Bx, double &By, double &Bz) const
vector< vector< double > > pp
virtual 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 GetBz(double x, double y, double z) const
virtual ~DMagneticFieldMapParameterized()
double Eval(double &r, double &z) const
DMagneticFieldMapParameterized(jana::JApplication *japp, string namepath="Magnets/Solenoid/solenoid_1500_poisson_20090814_01_params")
vector< vector< double > > cc
void Init(jana::JCalibration *jcalib, string namepath)