44 if (choiceModel==1)
G4cout<<
"-> Square quadrupole field"<<
G4endl;
45 if (choiceModel==2)
G4cout<<
"-> 3D quadrupole field"<<
G4endl;
46 if (choiceModel==3)
G4cout<<
"-> Enge quadrupole field"<<
G4endl;
60 const char * filename =
"OM50.grid";
63 G4cout <<
"\n-----------------------------------------------------------"
64 <<
"\n 3D Magnetic field from OPERA software "
65 <<
"\n-----------------------------------------------------------";
67 G4cout <<
"\n ---> " "Reading the field grid from " << filename <<
" ... " << endl;
68 ifstream
file( filename );
71 file >> fNx >> fNy >> fNz;
73 G4cout <<
" [ Number of values x,y,z: "
74 << fNx <<
" " << fNy <<
" " << fNz <<
" ] "
78 fXField.resize( fNx );
79 fYField.resize( fNx );
80 fZField.resize( fNx );
82 for (ix=0; ix<fNx; ix++) {
83 fXField[ix].resize(fNy);
84 fYField[ix].resize(fNy);
85 fZField[ix].resize(fNy);
86 for (iy=0; iy<fNy; iy++) {
87 fXField[ix][iy].resize(fNz);
88 fYField[ix][iy].resize(fNz);
89 fZField[ix][iy].resize(fNz);
94 double xval,yval,zval,bx,by,bz;
96 for (ix=0; ix<fNx; ix++) {
97 for (iy=0; iy<fNy; iy++) {
98 for (iz=0; iz<fNz; iz++) {
99 file >> xval >> yval >> zval >> bx >> by >> bz >> permeability;
100 if ( ix==0 && iy==0 && iz==0 ) {
101 fMinix = xval * lenUnit;
102 fMiniy = yval * lenUnit;
103 fMiniz = zval * lenUnit;
105 fXField[ix][iy][
iz] = bx ;
106 fYField[ix][iy][
iz] = by ;
107 fZField[ix][iy][
iz] = bz ;
113 fMaxix = xval * lenUnit;
114 fMaxiy = yval * lenUnit;
115 fMaxiz = zval * lenUnit;
117 G4cout <<
"\n ---> ... done reading " << endl;
120 G4cout <<
" ---> assumed the order: x, y, z, Bx, By, Bz "
121 <<
"\n ---> Min values x,y,z: "
122 << fMinix/
cm <<
" " << fMiniy/
cm <<
" " << fMiniz/
cm <<
" cm "
123 <<
"\n ---> Max values x,y,z: "
124 << fMaxix/
cm <<
" " << fMaxiy/
cm <<
" " << fMaxiz/
cm <<
" cm " << endl;
126 fDx = fMaxix - fMinix;
127 fDy = fMaxiy - fMiniy;
128 fDz = fMaxiz - fMiniz;
129 G4cout <<
"\n ---> Dif values x,y,z (range): "
130 << fDx/
cm <<
" " << fDy/
cm <<
" " << fDz/
cm <<
" cm in z "
131 <<
"\n-----------------------------------------------------------" << endl;
135 for (ix=0; ix<fNx; ix++)
137 for (iy=0; iy<fNy; iy++)
139 for (iz=0; iz<fNz; iz++)
142 fXField[ix][iy][
iz] = (fXField[ix][iy][
iz]/197.736);
143 fYField[ix][iy][
iz] = (fYField[ix][iy][
iz]/197.736);
144 fZField[ix][iy][
iz] = (fZField[ix][iy][
iz]/197.736);
158 double *Bfield )
const
186 gradient[0]=fGradient1*(
tesla/
m)/coef;
187 gradient[1]=fGradient2*(
tesla/
m)/coef;
188 gradient[2]=fGradient3*(
tesla/
m)/coef;
189 gradient[3]=fGradient4*(
tesla/
m)/coef;
190 gradient[4]=-fGradient3*(
tesla/
m)/coef;
192 for (quad=0; quad<=4; quad++)
194 if ((quad+1)==1) {z = point[2] + 3720 *
mm;}
195 if ((quad+1)==2) {z = point[2] + 3580 *
mm;}
196 if ((quad+1)==3) {z = point[2] + 330 *
mm;}
197 if ((quad+1)==4) {z = point[2] + 190 *
mm;}
198 if ((quad+1)==5) {z = point[2] + 50 *
mm;}
204 x>=fMinix && x<=fMaxix &&
205 y>=fMiniy && y<=fMaxiy &&
206 z>=fMiniz && z<=fMaxiz
211 double xfraction = (x - fMinix) / fDx;
212 double yfraction = (y - fMiniy) / fDy;
213 double zfraction = (z - fMiniz) / fDz;
217 double xdindex, ydindex, zdindex;
221 double xlocal = ( std::modf(xfraction*(fNx-1), &xdindex));
222 double ylocal = ( std::modf(yfraction*(fNy-1), &ydindex));
223 double zlocal = ( std::modf(zfraction*(fNz-1), &zdindex));
227 int xindex =
static_cast<int>(xdindex);
228 int yindex =
static_cast<int>(ydindex);
229 int zindex =
static_cast<int>(zdindex);
233 (fXField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
234 fXField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
235 fXField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
236 fXField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
237 fXField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
238 fXField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
239 fXField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
240 fXField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal)*gradient[quad]
244 (fYField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
245 fYField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
246 fYField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
247 fYField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
248 fYField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
249 fYField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
250 fYField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
251 fYField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal)*gradient[quad]
255 (fZField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
256 fZField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
257 fZField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
258 fZField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
259 fZField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
260 fZField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
261 fZField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
262 fZField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal)*gradient[quad]
293 if (z>=-3770*
mm && z<=-3670*
mm) G0 = (fGradient1/coef)*
tesla/
m;
294 if (z>=-3630*
mm && z<=-3530*
mm) G0 = (fGradient2/coef)*
tesla/
m;
296 if (z>=-380*
mm && z<=-280*
mm) G0 = (fGradient3/coef)*
tesla/
m;
297 if (z>=-240*
mm && z<=-140*
mm) G0 = (fGradient4/coef)*
tesla/
m;
298 if (z>=-100*
mm && z<=0*
mm) G0 = (-fGradient3/coef)*
tesla/
m;
330 G4double Grad1, Grad2, Grad3, Grad4, Grad5;
348 if ( (z>=-3900*
mm && z<-3470*
mm) || (z>=-490*
mm && z<100*
mm) )
355 G4double c0[5],
c1[5], c2[5], z1[5], z2[5], a0[5], gradient[5];
366 gradient[0] =Grad1*(
tesla/
m)/coef;
375 gradient[1] =Grad2*(
tesla/
m)/coef;
386 gradient[2] = Grad3*(
tesla/
m)/coef;
395 gradient[3] = Grad4*(
tesla/
m)/coef;
404 gradient[4] = Grad5*(
tesla/
m)/coef;
407 G4double Bx_local,By_local,Bz_local;
408 Bx_local = 0; By_local = 0; Bz_local = 0;
412 x_local= 0; y_local=0; z_local=0;
434 for (
G4int i=0;i<5; i++)
439 zoprime = lineZ + i*140*
mm;
442 z_local = (z - zoprime);
446 zoprime = lineZ + i*140*mm +(3150-40)*mm;
450 z_local = (z - zoprime);
454 if ( z_local < -z2[i] || z_local > z2[i])
462 if ( (z_local>=-z1[i]) && (z_local<=z1[i]) )
470 if ( ((z_local>=-z2[i]) && (z_local<-z1[i])) || ((z_local>z1[i]) && (z_local<=z2[i])) )
473 myVars = ( z_local - z1[i]) / a0[i];
474 if (z_local<-z1[i]) myVars = ( - z_local - z1[i]) / a0[i];
477 P0 = c0[i]+c1[i]*myVars+c2[i]*myVars*myVars;
479 P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a0[i]/a0[i];
480 if (z_local<-z1[i]) P1 = -c1[i]/a0[i]+2*c2[i]*(z_local+z1[i])/a0[i]/a0[i];
482 P2 = 2*c2[i]/a0[i]/a0[i];
484 cte = 1 + std::exp(c0[i]);
486 K1 = -cte*P1*std::exp(P0)/( (1+std::exp(P0))*(1+std::exp(P0)) );
488 K2 = -cte*std::exp(P0)*(
489 P2/( (1+std::exp(P0))*(1+std::exp(P0)) )
490 +2*P1*K1/(1+std::exp(P0))/cte
491 +P1*P1/(1+std::exp(P0))/(1+std::exp(P0))
494 K3 = -cte*std::exp(P0)*(
495 (3*P2*P1+P1*P1*P1)/(1+std::exp(P0))/(1+std::exp(P0))
496 +4*K1*(P1*P1+P2)/(1+std::exp(P0))/cte
497 +2*P1*(K1*K1/cte/cte+K2/(1+std::exp(P0))/cte)
500 G0 = gradient[i]*cte/(1+std::exp(P0));
507 Bx_local = y_local*(G0-(1./12)*(3*x_local*x_local+y_local*y_local)*G2);
508 By_local = x_local*(G0-(1./12)*(3*y_local*y_local+x_local*x_local)*G2);
509 Bz_local = x_local*y_local*(G1-(1./12)*(x_local*x_local+y_local*y_local)*G3);
void GetFieldValue(const double Point[4], double *Bfield) const
TabulatedField3D(G4float gr1, G4float gr2, G4float gr3, G4float gr4, G4int quadModel)
G4GLOB_DLL std::ostream G4cout