Geant4  10.01.p01
EMField.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708
30 // The Geant4-DNA web site is available at http://geant4-dna.org
31 //
32 // If you use this example, please cite the following publication:
33 // Rad. Prot. Dos. 133 (2009) 2-11
34 
35 #include "EMField.hh"
36 #include "G4SystemOfUnits.hh"
37 
38 EMField::EMField()
39 {
40 }
41 
42 void EMField::GetFieldValue(const double point[4], double *Bfield ) const
43 {
44  // Magnetic field
45  Bfield[0] = 0;
46  Bfield[1] = 0;
47  Bfield[2] = 0;
48 
49  // Electric field
50  Bfield[3] = 0;
51  Bfield[4] = 0;
52  Bfield[5] = 0;
53 
54  G4double Bx = 0;
55  G4double By = 0;
56  G4double Bz = 0;
57 
58  G4double x = point[0];
59  G4double y = point[1];
60  G4double z = point[2];
61 
62 // ***********************
63 // AIFIRA SWITCHING MAGNET
64 // ***********************
65 
66  // MAGNETIC FIELD VALUE FOR 3 MeV ALPHAS
67  // G4double switchingField = 0.0589768635 * tesla ;
68  G4double switchingField = 0.0590201 * tesla ;
69 
70  // BEAM START
71  G4double beamStart = -10*m;
72 
73  // RADIUS
74  G4double Rp = 0.698*m;
75 
76  // ENTRANCE POSITION AFTER ANALYSIS MAGNET
77  G4double zS = 975*mm;
78 
79  // POLE GAP
80  G4double D = 31.8*mm;
81 
82  // FRINGING FIELD
83 
84  G4double fieldBoundary, wc0, wc1, wc2, wc3, limitMinEntrance, limitMaxEntrance, limitMinExit, limitMaxExit;
85 
86  limitMinEntrance = beamStart+zS-4*D;
87  limitMaxEntrance = beamStart+zS+4*D;
88  limitMinExit =Rp-4*D;
89  limitMaxExit =Rp+4*D;
90 
91  wc0 = 0.3835;
92  wc1 = 2.388;
93  wc2 = -0.8171;
94  wc3 = 0.200;
95 
96  fieldBoundary=0.62;
97 
98  G4double ws, largeS, h, dhdlargeS, dhds, dlargeSds, dsdz, dsdx, zs0, Rs0, xcenter, zcenter;
99 
100 // - ENTRANCE OF SWITCHING MAGNET
101 
102 if ( (z >= limitMinEntrance) && (z < limitMaxEntrance) )
103 {
104  zs0 = fieldBoundary*D;
105  ws = (-z+beamStart+zS-zs0)/D;
106  dsdz = -1/D;
107  dsdx = 0;
108 
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;
114 
115  By = switchingField * h ;
116  Bx = y*switchingField*dhds*dsdx;
117  Bz = y*switchingField*dhds*dsdz;
118 
119 }
120 
121 // - HEART OF SWITCHING MAGNET
122 
123  if (
124  (z >= limitMaxEntrance)
125  && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS)) < limitMinExit*limitMinExit))
126  )
127 {
128  Bx=0;
129  By = switchingField;
130  Bz=0;
131 }
132 
133 // - EXIT OF SWITCHING MAGNET
134 
135 if (
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)
139 
140  )
141 {
142 
143  xcenter = 0;
144  zcenter = beamStart+zS;
145 
146  Rs0 = Rp + D*fieldBoundary;
147  ws = (std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter)) - Rs0)/D;
148 
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));
151 
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;
157 
158  By = switchingField * h ;
159  Bx = y*switchingField*dhds*dsdx;
160  Bz = y*switchingField*dhds*dsdz;
161 
162 }
163 
164 // **************************
165 // MICROBEAM LINE QUADRUPOLES
166 // **************************
167 
168  // MICROBEAM LINE ANGLE
169  G4double lineAngle = -10*deg;
170 
171  // X POSITION OF FIRST QUADRUPOLE
172  G4double lineX = -1295.59*mm;
173 
174  // Z POSITION OF FIRST QUADRUPOLE
175  G4double lineZ = -1327*mm;
176 
177  // Adjust magnetic zone absolute position
178  lineX = lineX + 5.24*micrometer*std::cos(-lineAngle); // 5.24 = 1.3 + 3.94 micrometer (cf. DetectorConstruction)
179  lineZ = lineZ + 5.24*micrometer*std::sin(-lineAngle);
180 
181  // QUADRUPOLE HALF LENGTH
182  G4double quadHalfLength = 75*mm;
183 
184  // QUADRUPOLE SPACING
185  G4double quadSpacing = 40*mm;
186 
187  // QUADRUPOLE CENTER COORDINATES
188  G4double xoprime, zoprime;
189 
190 if (z>=-1400*mm && z <-200*mm)
191 {
192  Bx=0; By=0; Bz=0;
193 
194  // FRINGING FILED CONSTANTS
195  G4double c0[4], c1[4], c2[4], z1[4], z2[4], a0[4], gradient[4];
196 
197  // QUADRUPOLE 1
198  c0[0] = -5.;
199  c1[0] = 2.5;
200  c2[0] = -0.1;
201  z1[0] = 60*mm;
202  z2[0] = 130*mm;
203  a0[0] = 10*mm;
204  gradient[0] = 3.406526 *tesla/m;
205 
206  // QUADRUPOLE 2
207  c0[1] = -5.;
208  c1[1] = 2.5;
209  c2[1] = -0.1;
210  z1[1] = 60*mm;
211  z2[1] = 130*mm;
212  a0[1] = 10*mm;
213  gradient[1] = -8.505263 *tesla/m;
214 
215  // QUADRUPOLE 3
216  c0[2] = -5.;
217  c1[2] = 2.5;
218  c2[2] = -0.1;
219  z1[2] = 60*mm;
220  z2[2] = 130*mm;
221  a0[2] = 10*mm;
222  gradient[2] = 8.505263 *tesla/m;
223 
224  // QUADRUPOLE 4
225  c0[3] = -5.;
226  c1[3] = 2.5;
227  c2[3] = -0.1;
228  z1[3] = 60*mm;
229  z2[3] = 130*mm;
230  a0[3] = 10*mm;
231  gradient[3] = -3.406526*tesla/m;
232 
233  // FIELD CREATED BY A QUADRUPOLE IN ITS LOCAL FRAME
234  G4double Bx_local,By_local,Bz_local;
235  Bx_local = 0; By_local = 0; Bz_local = 0;
236 
237  // FIELD CREATED BY A QUADRUPOOLE IN WORLD FRAME
238  G4double Bx_quad,By_quad,Bz_quad;
239  Bx_quad = 0; By_quad=0; Bz_quad=0;
240 
241  // QUADRUPOLE FRAME
242  G4double x_local,y_local,z_local;
243  x_local= 0; y_local=0; z_local=0;
244 
245  G4double vars = 0;
246  G4double G0, G1, G2, G3;
247  G4double K1, K2, K3;
248  G4double P0, P1, P2, cte;
249 
250  K1=0;
251  K2=0;
252  K3=0;
253  P0=0;
254  P1=0;
255  P2=0;
256  G0=0;
257  G1=0;
258  G2=0;
259  G3=0;
260  cte=0;
261 
262  G4bool largeScattering=false;
263 
264  for (G4int i=0;i<4; i++)
265  {
266 
267  if (i==0)
268  { xoprime = lineX + quadHalfLength*std::sin(lineAngle);
269  zoprime = lineZ + quadHalfLength*std::cos(lineAngle);
270 
271  x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
272  y_local = y;
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;
275 
276  }
277 
278  if (i==1)
279  { xoprime = lineX + (3*quadHalfLength+quadSpacing)*std::sin(lineAngle);
280  zoprime = lineZ + (3*quadHalfLength+quadSpacing)*std::cos(lineAngle);
281 
282  x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
283  y_local = y;
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;
286  }
287 
288  if (i==2)
289  { xoprime = lineX + (5*quadHalfLength+2*quadSpacing)*std::sin(lineAngle);
290  zoprime = lineZ + (5*quadHalfLength+2*quadSpacing)*std::cos(lineAngle);
291 
292  x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
293  y_local = y;
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;
296  }
297 
298  if (i==3)
299  { xoprime = lineX + (7*quadHalfLength+3*quadSpacing)*std::sin(lineAngle);
300  zoprime = lineZ + (7*quadHalfLength+3*quadSpacing)*std::cos(lineAngle);
301 
302  x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
303  y_local = y;
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;
306  }
307 
308 
309  if ( z_local < -z2[i] )
310  {
311  G0=0;
312  G1=0;
313  G2=0;
314  G3=0;
315  }
316 
317  if ( z_local > z2[i] )
318  {
319  G0=0;
320  G1=0;
321  G2=0;
322  G3=0;
323  }
324 
325  if ( (z_local>=-z1[i]) & (z_local<=z1[i]) )
326  {
327  G0=gradient[i];
328  G1=0;
329  G2=0;
330  G3=0;
331  }
332 
333  if ( ((z_local>=-z2[i]) & (z_local<-z1[i])) || ((z_local>z1[i]) & (z_local<=z2[i])) )
334  {
335 
336  vars = ( z_local - z1[i]) / a0[i] ;
337  if (z_local<-z1[i]) vars = ( - z_local - z1[i]) / a0[i] ;
338 
339 
340  P0 = c0[i]+c1[i]*vars+c2[i]*vars*vars;
341 
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];
344 
345  P2 = 2*c2[i]/a0[i]/a0[i];
346 
347  cte = 1 + std::exp(c0[i]);
348 
349  K1 = -cte*P1*std::exp(P0)/( (1+std::exp(P0))*(1+std::exp(P0)) );
350 
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))
355  );
356 
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)
361  );
362 
363  G0 = gradient[i]*cte/(1+std::exp(P0));
364  G1 = gradient[i]*K1;
365  G2 = gradient[i]*K2;
366  G3 = gradient[i]*K3;
367 
368  }
369 
370  // PROTECTION AGAINST LARGE SCATTERING
371 
372  if ( largeScattering )
373  {
374  G0=0;
375  G1=0;
376  G2=0;
377  G3=0;
378  }
379 
380  // MAGNETIC FIELD COMPUTATION FOR EACH QUADRUPOLE
381 
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);
385 
386  Bx_quad = Bz_local*std::sin(lineAngle)+Bx_local*std::cos(lineAngle);
387  By_quad = By_local;
388  Bz_quad = Bz_local*std::cos(lineAngle)-Bx_local*std::sin(lineAngle);
389 
390  // TOTAL MAGNETIC FIELD
391 
392  Bx = Bx + Bx_quad ;
393  By = By + By_quad ;
394  Bz = Bz + Bz_quad ;
395 
396  } // LOOP ON QUADRUPOLES
397 
398 
399 } // END OF QUADRUPLET
400 
401  Bfield[0] = Bx;
402  Bfield[1] = By;
403  Bfield[2] = Bz;
404 
405 // *****************************************
406 // ELECTRIC FIELD CREATED BY SCANNING PLATES
407 // *****************************************
408 
409  Bfield[3] = 0;
410  Bfield[4] = 0;
411  Bfield[5] = 0;
412 
413  // POSITION OF EXIT OF LAST QUAD WHERE THE SCANNING PLATES START
414 
415  G4double electricPlateWidth1 = 5 * mm;
416  G4double electricPlateWidth2 = 5 * mm;
417  G4double electricPlateLength1 = 36 * mm;
418  G4double electricPlateLength2 = 34 * mm;
419  G4double electricPlateGap = 5 * mm;
420  G4double electricPlateSpacing1 = 3 * mm;
421  G4double electricPlateSpacing2 = 4 * mm;
422 
423  // APPLY VOLTAGE HERE IN VOLTS (no electrostatic deflection here)
424  G4double electricPlateVoltage1 = 0 * volt;
425  G4double electricPlateVoltage2 = 0 * volt;
426 
427  G4double electricFieldPlate1 = electricPlateVoltage1 / electricPlateSpacing1 ;
428  G4double electricFieldPlate2 = electricPlateVoltage2 / electricPlateSpacing2 ;
429 
430  G4double beginFirstZoneX = lineX + (8*quadHalfLength+3*quadSpacing)*std::sin(lineAngle);
431  G4double beginFirstZoneZ = lineZ + (8*quadHalfLength+3*quadSpacing)*std::cos(lineAngle);
432 
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);
435 
436  G4double xA, zA, xB, zB, xC, zC, xD, zD;
437  G4double slope1, cte1, slope2, cte2, slope3, cte3, slope4, cte4;
438 
439  // WARNING : lineAngle < 0
440 
441  // FIRST PLATES
442 
443  xA = beginFirstZoneX + std::cos(lineAngle)*electricPlateSpacing1/2;
444  zA = beginFirstZoneZ - std::sin(lineAngle)*electricPlateSpacing1/2;
445 
446  xB = xA + std::sin(lineAngle)*electricPlateLength1;
447  zB = zA + std::cos(lineAngle)*electricPlateLength1;
448 
449  xC = xB - std::cos(lineAngle)*electricPlateSpacing1;
450  zC = zB + std::sin(lineAngle)*electricPlateSpacing1;
451 
452  xD = xC - std::sin(lineAngle)*electricPlateLength1;
453  zD = zC - std::cos(lineAngle)*electricPlateLength1;
454 
455  slope1 = (xB-xA)/(zB-zA);
456  cte1 = xA - slope1 * zA;
457 
458  slope2 = (xC-xB)/(zC-zB);
459  cte2 = xB - slope2 * zB;
460 
461  slope3 = (xD-xC)/(zD-zC);
462  cte3 = xC - slope3 * zC;
463 
464  slope4 = (xA-xD)/(zA-zD);
465  cte4 = xD - slope4 * zD;
466 
467 
468  if
469  (
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
475  )
476 
477  {
478  Bfield[3] = electricFieldPlate1*std::cos(lineAngle);
479  Bfield[4] = 0;
480  Bfield[5] = -electricFieldPlate1*std::sin(lineAngle);
481 
482  }
483 
484  // SECOND PLATES
485 
486  xA = beginSecondZoneX + std::cos(lineAngle)*electricPlateWidth2/2;
487  zA = beginSecondZoneZ - std::sin(lineAngle)*electricPlateWidth2/2;
488 
489  xB = xA + std::sin(lineAngle)*electricPlateLength2;
490  zB = zA + std::cos(lineAngle)*electricPlateLength2;
491 
492  xC = xB - std::cos(lineAngle)*electricPlateWidth2;
493  zC = zB + std::sin(lineAngle)*electricPlateWidth2;
494 
495  xD = xC - std::sin(lineAngle)*electricPlateLength2;
496  zD = zC - std::cos(lineAngle)*electricPlateLength2;
497 
498  slope1 = (xB-xA)/(zB-zA);
499  cte1 = xA - slope1 * zA;
500 
501  slope2 = (xC-xB)/(zC-zB);
502  cte2 = xB - slope2 * zB;
503 
504  slope3 = (xD-xC)/(zD-zC);
505  cte3 = xC - slope3 * zC;
506 
507  slope4 = (xA-xD)/(zA-zD);
508  cte4 = xD - slope4 * zD;
509 
510  if
511  (
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
517  )
518 
519  {
520  Bfield[3] = 0;
521  Bfield[4] = electricFieldPlate2;
522  Bfield[5] = 0;
523  }
524 
525 //
526 
527 }
a0
const G4double a0
Definition: G4IonCoulombCrossSection.cc:68
EMField::GetFieldValue
void GetFieldValue(const double Point[4], double *Bfield) const
Definition: EMField.cc:42
c2
static c2_factory< G4double > c2
Definition: G4ScreenedNuclearRecoil.cc:128
P1
static const G4double * P1[nN]
Definition: G4ElectroNuclearCrossSection.cc:2176
z
G4double z
Definition: TRTMaterials.hh:39
volt
static const double volt
Definition: G4SIunits.hh:223
EMField::EMField
EMField()
Definition: EMField.cc:38
G4int
int G4int
Definition: G4Types.hh:78
P0
static const G4double * P0[nN]
Definition: G4ElectroNuclearCrossSection.cc:2175
deg
static const double deg
Definition: G4SIunits.hh:133
G4bool
bool G4bool
Definition: G4Types.hh:79
c1
static const G4double c1
Definition: G4BetheHeitlerModel.cc:80
micrometer
static const double micrometer
Definition: G4SIunits.hh:90
c0
static const G4double c0
Definition: G4BetheHeitlerModel.cc:80
P2
static const G4double * P2[nN]
Definition: G4ElectroNuclearCrossSection.cc:2177
m
static const double m
Definition: G4SIunits.hh:110
tesla
static const double tesla
Definition: G4SIunits.hh:247
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
mm
static const double mm
Definition: G4SIunits.hh:102
EMField.hh