55 if (choiceModel==1)
G4cout<<
"-> Square quadrupole field"<<
G4endl;
56 if (choiceModel==2)
G4cout<<
"-> 3D quadrupole field"<<
G4endl;
57 if (choiceModel==3)
G4cout<<
"-> Enge quadrupole field"<<
G4endl;
80 G4cout <<
"\n-----------------------------------------------------------"
81 <<
"\n 3D Magnetic field from OPERA software "
82 <<
"\n-----------------------------------------------------------";
84 G4cout <<
"\n ---> " "Reading the field grid from " << filename <<
" ... " <<
G4endl;
85 std::ifstream
file( filename );
90 G4cout <<
" [ Number of values x,y,z: "
91 <<
fNx <<
" " <<
fNy <<
" " << fNz <<
" ] "
99 for (ix=0; ix<
fNx; ix++)
104 for (iy=0; iy<
fNy; iy++)
113 double xval,yval,zval,bx,by,bz;
115 for (ix=0; ix<
fNx; ix++)
117 for (iy=0; iy<
fNy; iy++)
119 for (iz=0; iz<
fNz; iz++)
121 file >> xval >> yval >> zval >> bx >> by >> bz >> permeability;
122 if ( ix==0 && iy==0 && iz==0 )
148 G4cout <<
" ---> assumed the order: x, y, z, Bx, By, Bz "
149 <<
"\n ---> Min values x,y,z: "
151 <<
"\n ---> Max values x,y,z: "
158 G4cout <<
"\n ---> Dif values x,y,z (range): "
160 <<
"\n-----------------------------------------------------------" <<
G4endl;
165 for (ix=0; ix<
fNx; ix++)
167 for (iy=0; iy<
fNy; iy++)
169 for (iz=0; iz<
fNz; iz++)
188 double *Bfield )
const
228 for (quad=0; quad<=4; quad++)
230 if ((quad+1)==1) {z = point[2] + 3720 *
mm;}
231 if ((quad+1)==2) {z = point[2] + 3580 *
mm;}
232 if ((quad+1)==3) {z = point[2] + 330 *
mm;}
233 if ((quad+1)==4) {z = point[2] + 190 *
mm;}
234 if ((quad+1)==5) {z = point[2] + 50 *
mm;}
253 double xdindex, ydindex, zdindex;
257 double xlocal = ( std::modf(xfraction*(
fNx-1), &xdindex));
258 double ylocal = ( std::modf(yfraction*(
fNy-1), &ydindex));
259 double zlocal = ( std::modf(zfraction*(
fNz-1), &zdindex));
269 int xindex =
static_cast<int>(std::floor(xdindex));
270 int yindex =
static_cast<int>(std::floor(ydindex));
271 int zindex =
static_cast<int>(std::floor(zdindex));
275 (
fXField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
276 fXField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
277 fXField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
278 fXField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
279 fXField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
280 fXField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
281 fXField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
282 fXField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal)*gradient[quad]
286 (
fYField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
287 fYField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
288 fYField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
289 fYField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
290 fYField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
291 fYField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
292 fYField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
293 fYField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal)*gradient[quad]
297 (
fZField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
298 fZField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
299 fZField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
300 fZField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
301 fZField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
302 fZField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
303 fZField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
304 fZField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal)*gradient[quad]
372 G4double Grad1, Grad2, Grad3, Grad4, Grad5;
390 if ( (z>=-3900*
mm && z<-3470*
mm) || (z>=-490*
mm && z<100*
mm) )
408 gradient[0] =Grad1*(
tesla/
m)/coef;
417 gradient[1] =Grad2*(
tesla/
m)/coef;
428 gradient[2] = Grad3*(
tesla/
m)/coef;
437 gradient[3] = Grad4*(
tesla/
m)/coef;
446 gradient[4] = Grad5*(
tesla/
m)/coef;
449 G4double Bx_local,By_local,Bz_local;
450 Bx_local = 0; By_local = 0; Bz_local = 0;
454 x_local= 0; y_local=0; z_local=0;
476 for (
G4int i=0;i<5; i++)
481 zoprime = lineZ + i*140*
mm;
484 z_local = (z - zoprime);
488 zoprime = lineZ + i*140*mm +(3150-40)*mm;
492 z_local = (z - zoprime);
496 if ( z_local < -z2[i] || z_local > z2[i])
504 if ( (z_local>=-z1[i]) && (z_local<=z1[i]) )
512 if ( ((z_local>=-z2[i]) && (z_local<-z1[i])) || ((z_local>z1[i]) && (z_local<=z2[i])) )
515 myVars = ( z_local - z1[i]) / a0[i];
516 if (z_local<-z1[i]) myVars = ( - z_local - z1[i]) / a0[i];
519 P0 = c0[i]+c1[i]*myVars+c2[i]*myVars*myVars;
521 P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a0[i]/a0[i];
522 if (z_local<-z1[i]) P1 = -c1[i]/a0[i]+2*c2[i]*(z_local+z1[i])/a0[i]/a0[i];
524 P2 = 2*c2[i]/a0[i]/a0[i];
526 cte = 1 +
G4Exp(c0[i]);
530 K2 = -cte*
G4Exp(P0)*(
532 +2*P1*K1/(1+
G4Exp(P0))/cte
536 K3 = -cte*
G4Exp(P0)*(
538 +4*K1*(P1*P1+P2)/(1+
G4Exp(P0))/cte
539 +2*P1*(K1*K1/cte/cte+K2/(1+
G4Exp(P0))/cte)
542 G0 = gradient[i]*cte/(1+
G4Exp(P0));
549 Bx_local = y_local*(G0-(1./12)*(3*x_local*x_local+y_local*y_local)*G2);
550 By_local = x_local*(G0-(1./12)*(3*y_local*y_local+x_local*x_local)*G2);
551 Bz_local = x_local*y_local*(G1-(1./12)*(x_local*x_local+y_local*y_local)*G3);