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 );
73 G4cout <<
" [ Number of values x,y,z: "
74 <<
fNx <<
" " <<
fNy <<
" " << fNz <<
" ] "
82 for (ix=0; ix<
fNx; ix++) {
86 for (iy=0; iy<
fNy; iy++) {
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 ) {
117 G4cout <<
"\n ---> ... done reading " << endl;
120 G4cout <<
" ---> assumed the order: x, y, z, Bx, By, Bz "
121 <<
"\n ---> Min values x,y,z: "
123 <<
"\n ---> Max values x,y,z: "
129 G4cout <<
"\n ---> Dif values x,y,z (range): "
131 <<
"\n-----------------------------------------------------------" << endl;
135 for (ix=0; ix<
fNx; ix++)
137 for (iy=0; iy<
fNy; iy++)
139 for (iz=0; iz<
fNz; iz++)
158 double *Bfield )
const
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;}
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]
330 G4double Grad1, Grad2, Grad3, Grad4, Grad5;
348 if ( (z>=-3900*
mm && z<-3470*
mm) || (z>=-490*
mm && z<100*
mm) )
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);
static c2_factory< G4double > c2
vector< vector< vector< double > > > fXField
static const G4double * P1[nN]
void GetFieldValue(const double Point[4], double *Bfield) const
vector< vector< vector< double > > > fZField
static const G4double * P0[nN]
TabulatedField3D(G4float gr1, G4float gr2, G4float gr3, G4float gr4, G4int quadModel)
G4GLOB_DLL std::ostream G4cout
vector< vector< vector< double > > > fYField
const G4double x[NPOINTSGL]
static const G4double * P2[nN]
static const double tesla