84 G4double fieldBoundary, wc0, wc1, wc2, wc3, limitMinEntrance, limitMaxEntrance, limitMinExit, limitMaxExit;
86 limitMinEntrance = beamStart+zS-4*
D;
87 limitMaxEntrance = beamStart+zS+4*
D;
98 G4double ws, largeS, h, dhdlargeS, dhds, dlargeSds, dsdz, dsdx, zs0, Rs0, xcenter, zcenter;
102 if ( (z >= limitMinEntrance) && (z < limitMaxEntrance) )
104 zs0 = fieldBoundary*
D;
105 ws = (-z+beamStart+zS-zs0)/D;
109 largeS = wc0 + wc1*ws + wc2*ws*ws + wc3*ws*ws*ws;
110 h = 1./(1.+std::exp(largeS));
111 dhdlargeS = -std::exp(largeS)*h*h;
112 dlargeSds = wc1+ 2*wc2*ws + 3*wc3*ws*ws;
113 dhds = dhdlargeS * dlargeSds;
115 By = switchingField * h ;
116 Bx = y*switchingField*dhds*dsdx;
117 Bz = y*switchingField*dhds*dsdz;
124 (z >= limitMaxEntrance)
125 && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS)) < limitMinExit*limitMinExit))
136 (z >= limitMaxEntrance)
137 && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS))) >= limitMinExit*limitMinExit)
138 && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS))) < limitMaxExit*limitMaxExit)
144 zcenter = beamStart+zS;
146 Rs0 = Rp + D*fieldBoundary;
147 ws = (std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter)) - Rs0)/D;
149 dsdz = (1/
D)*(z-zcenter)/std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter));
150 dsdx = (1/
D)*(x-xcenter)/std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter));
152 largeS = wc0 + wc1*ws + wc2*ws*ws + wc3*ws*ws*ws;
153 h = 1./(1.+std::exp(largeS));
154 dhdlargeS = -std::exp(largeS)*h*h;
155 dlargeSds = wc1+ 2*wc2*ws + 3*wc3*ws*ws;
156 dhds = dhdlargeS * dlargeSds;
158 By = switchingField * h ;
159 Bx = y*switchingField*dhds*dsdx;
160 Bz = y*switchingField*dhds*dsdz;
178 lineX = lineX + 5.24*
micrometer*std::cos(-lineAngle);
179 lineZ = lineZ + 5.24*
micrometer*std::sin(-lineAngle);
190 if (z>=-1400*
mm && z <-200*
mm)
204 gradient[0] = 3.406526 *
tesla/
m;
213 gradient[1] = -8.505263 *
tesla/
m;
222 gradient[2] = 8.505263 *
tesla/
m;
231 gradient[3] = -3.406526*
tesla/
m;
234 G4double Bx_local,By_local,Bz_local;
235 Bx_local = 0; By_local = 0; Bz_local = 0;
239 Bx_quad = 0; By_quad=0; Bz_quad=0;
243 x_local= 0; y_local=0; z_local=0;
262 G4bool largeScattering=
false;
264 for (
G4int i=0;i<4; i++)
268 { xoprime = lineX + quadHalfLength*std::sin(lineAngle);
269 zoprime = lineZ + quadHalfLength*std::cos(lineAngle);
271 x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
273 z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
274 if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=
true;
279 { xoprime = lineX + (3*quadHalfLength+quadSpacing)*std::sin(lineAngle);
280 zoprime = lineZ + (3*quadHalfLength+quadSpacing)*std::cos(lineAngle);
282 x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
284 z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
285 if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=
true;
289 { xoprime = lineX + (5*quadHalfLength+2*quadSpacing)*std::sin(lineAngle);
290 zoprime = lineZ + (5*quadHalfLength+2*quadSpacing)*std::cos(lineAngle);
292 x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
294 z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
295 if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=
true;
299 { xoprime = lineX + (7*quadHalfLength+3*quadSpacing)*std::sin(lineAngle);
300 zoprime = lineZ + (7*quadHalfLength+3*quadSpacing)*std::cos(lineAngle);
302 x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
304 z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
305 if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=
true;
309 if ( z_local < -z2[i] )
317 if ( z_local > z2[i] )
325 if ( (z_local>=-z1[i]) & (z_local<=z1[i]) )
333 if ( ((z_local>=-z2[i]) & (z_local<-z1[i])) || ((z_local>z1[i]) & (z_local<=z2[i])) )
336 vars = ( z_local - z1[i]) / a0[i] ;
337 if (z_local<-z1[i]) vars = ( - z_local - z1[i]) / a0[i] ;
340 P0 = c0[i]+c1[i]*vars+c2[i]*vars*vars;
342 P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a0[i]/a0[i];
343 if (z_local<-z1[i]) P1 = -c1[i]/a0[i]+2*c2[i]*(z_local+z1[i])/a0[i]/a0[i];
345 P2 = 2*c2[i]/a0[i]/a0[i];
347 cte = 1 + std::exp(c0[i]);
349 K1 = -cte*P1*std::exp(P0)/( (1+std::exp(P0))*(1+std::exp(P0)) );
351 K2 = -cte*std::exp(P0)*(
352 P2/( (1+std::exp(P0))*(1+std::exp(P0)) )
353 +2*P1*K1/(1+std::exp(P0))/cte
354 +P1*P1/(1+std::exp(P0))/(1+std::exp(P0))
357 K3 = -cte*std::exp(P0)*(
358 (3*P2*P1+P1*P1*
P1)/(1+std::exp(P0))/(1+std::exp(P0))
359 +4*K1*(P1*P1+P2)/(1+std::exp(P0))/cte
360 +2*P1*(K1*K1/cte/cte+K2/(1+std::exp(P0))/cte)
363 G0 = gradient[i]*cte/(1+std::exp(P0));
372 if ( largeScattering )
382 Bx_local = y_local*(G0-(1./12)*(3*x_local*x_local+y_local*y_local)*G2);
383 By_local = x_local*(G0-(1./12)*(3*y_local*y_local+x_local*x_local)*G2);
384 Bz_local = x_local*y_local*(G1-(1./12)*(x_local*x_local+y_local*y_local)*G3);
386 Bx_quad = Bz_local*std::sin(lineAngle)+Bx_local*std::cos(lineAngle);
388 Bz_quad = Bz_local*std::cos(lineAngle)-Bx_local*std::sin(lineAngle);
427 G4double electricFieldPlate1 = electricPlateVoltage1 / electricPlateSpacing1 ;
428 G4double electricFieldPlate2 = electricPlateVoltage2 / electricPlateSpacing2 ;
430 G4double beginFirstZoneX = lineX + (8*quadHalfLength+3*quadSpacing)*std::sin(lineAngle);
431 G4double beginFirstZoneZ = lineZ + (8*quadHalfLength+3*quadSpacing)*std::cos(lineAngle);
433 G4double beginSecondZoneX = lineX + (8*quadHalfLength+3*quadSpacing+electricPlateLength1+electricPlateGap)*std::sin(lineAngle);
434 G4double beginSecondZoneZ = lineZ + (8*quadHalfLength+3*quadSpacing+electricPlateLength1+electricPlateGap)*std::cos(lineAngle);
436 G4double xA, zA, xB, zB, xC, zC, xD, zD;
437 G4double slope1, cte1, slope2, cte2, slope3, cte3, slope4, cte4;
443 xA = beginFirstZoneX + std::cos(lineAngle)*electricPlateSpacing1/2;
444 zA = beginFirstZoneZ - std::sin(lineAngle)*electricPlateSpacing1/2;
446 xB = xA + std::sin(lineAngle)*electricPlateLength1;
447 zB = zA + std::cos(lineAngle)*electricPlateLength1;
449 xC = xB - std::cos(lineAngle)*electricPlateSpacing1;
450 zC = zB + std::sin(lineAngle)*electricPlateSpacing1;
452 xD = xC - std::sin(lineAngle)*electricPlateLength1;
453 zD = zC - std::cos(lineAngle)*electricPlateLength1;
455 slope1 = (xB-xA)/(zB-zA);
456 cte1 = xA - slope1 * zA;
458 slope2 = (xC-xB)/(zC-zB);
459 cte2 = xB - slope2 * zB;
461 slope3 = (xD-xC)/(zD-zC);
462 cte3 = xC - slope3 * zC;
464 slope4 = (xA-xD)/(zA-zD);
465 cte4 = xD - slope4 * zD;
470 x <= slope1 * z + cte1
471 && x >= slope3 * z + cte3
472 && x <= slope4 * z + cte4
473 && x >= slope2 * z + cte2
474 && std::abs(y)<=electricPlateWidth1/2
478 Bfield[3] = electricFieldPlate1*std::cos(lineAngle);
480 Bfield[5] = -electricFieldPlate1*std::sin(lineAngle);
486 xA = beginSecondZoneX + std::cos(lineAngle)*electricPlateWidth2/2;
487 zA = beginSecondZoneZ - std::sin(lineAngle)*electricPlateWidth2/2;
489 xB = xA + std::sin(lineAngle)*electricPlateLength2;
490 zB = zA + std::cos(lineAngle)*electricPlateLength2;
492 xC = xB - std::cos(lineAngle)*electricPlateWidth2;
493 zC = zB + std::sin(lineAngle)*electricPlateWidth2;
495 xD = xC - std::sin(lineAngle)*electricPlateLength2;
496 zD = zC - std::cos(lineAngle)*electricPlateLength2;
498 slope1 = (xB-xA)/(zB-zA);
499 cte1 = xA - slope1 * zA;
501 slope2 = (xC-xB)/(zC-zB);
502 cte2 = xB - slope2 * zB;
504 slope3 = (xD-xC)/(zD-zC);
505 cte3 = xC - slope3 * zC;
507 slope4 = (xA-xD)/(zA-zD);
508 cte4 = xD - slope4 * zD;
512 x <= slope1 * z + cte1
513 && x >= slope3 * z + cte3
514 && x <= slope4 * z + cte4
515 && x >= slope2 * z + cte2
516 && std::abs(y)<=electricPlateSpacing2/2
521 Bfield[4] = electricFieldPlate2;
void GetFieldValue(const double Point[4], double *Bfield) const
static c2_factory< G4double > c2
static const G4double * P1[nN]
static const G4double * P0[nN]
static const double micrometer
const G4double x[NPOINTSGL]
static const G4double * P2[nN]
static const double tesla