Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
EMField.cc
이 파일의 문서화 페이지로 가기
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 // Based on purging magnet advanced example.
36 //
37 
38 #include "EMField.hh"
39 #include "G4Exp.hh"
40 #include "G4SystemOfUnits.hh"
41 
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43 
44 
45 EMField::EMField()
46 {}
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 
50 void EMField::GetFieldValue(const double point[4], double *Bfield ) const
51 {
52  // Magnetic field
53  Bfield[0] = 0;
54  Bfield[1] = 0;
55  Bfield[2] = 0;
56 
57  // Electric field
58  Bfield[3] = 0;
59  Bfield[4] = 0;
60  Bfield[5] = 0;
61 
62  G4double Bx = 0;
63  G4double By = 0;
64  G4double Bz = 0;
65 
66  G4double x = point[0];
67  G4double y = point[1];
68  G4double z = point[2];
69 
70 // ***********************
71 // AIFIRA SWITCHING MAGNET
72 // ***********************
73 
74  // MAGNETIC FIELD VALUE FOR 3 MeV ALPHAS
75  // G4double switchingField = 0.0589768635 * tesla ;
76  G4double switchingField = 0.0590201 * tesla ;
77 
78  // BEAM START
79  G4double beamStart = -10*m;
80 
81  // RADIUS
82  G4double Rp = 0.698*m;
83 
84  // ENTRANCE POSITION AFTER ANALYSIS MAGNET
85  G4double zS = 975*mm;
86 
87  // POLE GAP
88  G4double D = 31.8*mm;
89 
90  // FRINGING FIELD
91 
92  G4double fieldBoundary, wc0, wc1, wc2, wc3, limitMinEntrance, limitMaxEntrance, limitMinExit, limitMaxExit;
93 
94  limitMinEntrance = beamStart+zS-4*D;
95  limitMaxEntrance = beamStart+zS+4*D;
96  limitMinExit =Rp-4*D;
97  limitMaxExit =Rp+4*D;
98 
99  wc0 = 0.3835;
100  wc1 = 2.388;
101  wc2 = -0.8171;
102  wc3 = 0.200;
103 
104  fieldBoundary=0.62;
105 
106  G4double ws, largeS, h, dhdlargeS, dhds, dlargeSds, dsdz, dsdx, zs0, Rs0, xcenter, zcenter;
107 
108 // - ENTRANCE OF SWITCHING MAGNET
109 
110 if ( (z >= limitMinEntrance) && (z < limitMaxEntrance) )
111 {
112  zs0 = fieldBoundary*D;
113  ws = (-z+beamStart+zS-zs0)/D;
114  dsdz = -1/D;
115  dsdx = 0;
116 
117  largeS = wc0 + wc1*ws + wc2*ws*ws + wc3*ws*ws*ws;
118  h = 1./(1.+G4Exp(largeS));
119  dhdlargeS = -G4Exp(largeS)*h*h;
120  dlargeSds = wc1+ 2*wc2*ws + 3*wc3*ws*ws;
121  dhds = dhdlargeS * dlargeSds;
122 
123  By = switchingField * h ;
124  Bx = y*switchingField*dhds*dsdx;
125  Bz = y*switchingField*dhds*dsdz;
126 
127 }
128 
129 // - HEART OF SWITCHING MAGNET
130 
131  if (
132  (z >= limitMaxEntrance)
133  && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS)) < limitMinExit*limitMinExit))
134  )
135 {
136  Bx=0;
137  By = switchingField;
138  Bz=0;
139 }
140 
141 // - EXIT OF SWITCHING MAGNET
142 
143 if (
144  (z >= limitMaxEntrance)
145  && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS))) >= limitMinExit*limitMinExit)
146  && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS))) < limitMaxExit*limitMaxExit)
147 
148  )
149 {
150 
151  xcenter = 0;
152  zcenter = beamStart+zS;
153 
154  Rs0 = Rp + D*fieldBoundary;
155  ws = (std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter)) - Rs0)/D;
156 
157  dsdz = (1/D)*(z-zcenter)/std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter));
158  dsdx = (1/D)*(x-xcenter)/std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter));
159 
160  largeS = wc0 + wc1*ws + wc2*ws*ws + wc3*ws*ws*ws;
161  h = 1./(1.+G4Exp(largeS));
162  dhdlargeS = -G4Exp(largeS)*h*h;
163  dlargeSds = wc1+ 2*wc2*ws + 3*wc3*ws*ws;
164  dhds = dhdlargeS * dlargeSds;
165 
166  By = switchingField * h ;
167  Bx = y*switchingField*dhds*dsdx;
168  Bz = y*switchingField*dhds*dsdz;
169 
170 }
171 
172 // **************************
173 // MICROBEAM LINE QUADRUPOLES
174 // **************************
175 
176  // MICROBEAM LINE ANGLE
177  G4double lineAngle = -10*deg;
178 
179  // X POSITION OF FIRST QUADRUPOLE
180  G4double lineX = -1295.59*mm;
181 
182  // Z POSITION OF FIRST QUADRUPOLE
183  G4double lineZ = -1327*mm;
184 
185  // Adjust magnetic zone absolute position
186  lineX = lineX + 5.24*micrometer*std::cos(-lineAngle); // 5.24 = 1.3 + 3.94 micrometer (cf. DetectorConstruction)
187  lineZ = lineZ + 5.24*micrometer*std::sin(-lineAngle);
188 
189  // QUADRUPOLE HALF LENGTH
190  G4double quadHalfLength = 75*mm;
191 
192  // QUADRUPOLE SPACING
193  G4double quadSpacing = 40*mm;
194 
195  // QUADRUPOLE CENTER COORDINATES
196  G4double xoprime, zoprime;
197 
198 if (z>=-1400*mm && z <-200*mm)
199 {
200  Bx=0; By=0; Bz=0;
201 
202  // FRINGING FILED CONSTANTS
203  G4double c0[4], c1[4], c2[4], z1[4], z2[4], a0[4], gradient[4];
204 
205  // QUADRUPOLE 1
206  c0[0] = -5.;
207  c1[0] = 2.5;
208  c2[0] = -0.1;
209  z1[0] = 60*mm;
210  z2[0] = 130*mm;
211  a0[0] = 10*mm;
212  gradient[0] = 3.406526 *tesla/m;
213 
214  // QUADRUPOLE 2
215  c0[1] = -5.;
216  c1[1] = 2.5;
217  c2[1] = -0.1;
218  z1[1] = 60*mm;
219  z2[1] = 130*mm;
220  a0[1] = 10*mm;
221  gradient[1] = -8.505263 *tesla/m;
222 
223  // QUADRUPOLE 3
224  c0[2] = -5.;
225  c1[2] = 2.5;
226  c2[2] = -0.1;
227  z1[2] = 60*mm;
228  z2[2] = 130*mm;
229  a0[2] = 10*mm;
230  gradient[2] = 8.505263 *tesla/m;
231 
232  // QUADRUPOLE 4
233  c0[3] = -5.;
234  c1[3] = 2.5;
235  c2[3] = -0.1;
236  z1[3] = 60*mm;
237  z2[3] = 130*mm;
238  a0[3] = 10*mm;
239  gradient[3] = -3.406526*tesla/m;
240 
241  // FIELD CREATED BY A QUADRUPOLE IN ITS LOCAL FRAME
242  G4double Bx_local,By_local,Bz_local;
243  Bx_local = 0; By_local = 0; Bz_local = 0;
244 
245  // FIELD CREATED BY A QUADRUPOOLE IN WORLD FRAME
246  G4double Bx_quad,By_quad,Bz_quad;
247  Bx_quad = 0; By_quad=0; Bz_quad=0;
248 
249  // QUADRUPOLE FRAME
250  G4double x_local,y_local,z_local;
251  x_local= 0; y_local=0; z_local=0;
252 
253  G4double vars = 0;
254  G4double G0, G1, G2, G3;
255  G4double K1, K2, K3;
256  G4double P0, P1, P2, cte;
257 
258  K1=0;
259  K2=0;
260  K3=0;
261  P0=0;
262  P1=0;
263  P2=0;
264  G0=0;
265  G1=0;
266  G2=0;
267  G3=0;
268  cte=0;
269 
270  G4bool largeScattering=false;
271 
272  for (G4int i=0;i<4; i++)
273  {
274 
275  if (i==0)
276  { xoprime = lineX + quadHalfLength*std::sin(lineAngle);
277  zoprime = lineZ + quadHalfLength*std::cos(lineAngle);
278 
279  x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
280  y_local = y;
281  z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
282  if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
283 
284  }
285 
286  if (i==1)
287  { xoprime = lineX + (3*quadHalfLength+quadSpacing)*std::sin(lineAngle);
288  zoprime = lineZ + (3*quadHalfLength+quadSpacing)*std::cos(lineAngle);
289 
290  x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
291  y_local = y;
292  z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
293  if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
294  }
295 
296  if (i==2)
297  { xoprime = lineX + (5*quadHalfLength+2*quadSpacing)*std::sin(lineAngle);
298  zoprime = lineZ + (5*quadHalfLength+2*quadSpacing)*std::cos(lineAngle);
299 
300  x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
301  y_local = y;
302  z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
303  if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
304  }
305 
306  if (i==3)
307  { xoprime = lineX + (7*quadHalfLength+3*quadSpacing)*std::sin(lineAngle);
308  zoprime = lineZ + (7*quadHalfLength+3*quadSpacing)*std::cos(lineAngle);
309 
310  x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
311  y_local = y;
312  z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
313  if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
314  }
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 > z2[i] )
326  {
327  G0=0;
328  G1=0;
329  G2=0;
330  G3=0;
331  }
332 
333  if ( (z_local>=-z1[i]) & (z_local<=z1[i]) )
334  {
335  G0=gradient[i];
336  G1=0;
337  G2=0;
338  G3=0;
339  }
340 
341  if ( ((z_local>=-z2[i]) & (z_local<-z1[i])) || ((z_local>z1[i]) & (z_local<=z2[i])) )
342  {
343 
344  vars = ( z_local - z1[i]) / a0[i] ;
345  if (z_local<-z1[i]) vars = ( - z_local - z1[i]) / a0[i] ;
346 
347 
348  P0 = c0[i]+c1[i]*vars+c2[i]*vars*vars;
349 
350  P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a0[i]/a0[i];
351  if (z_local<-z1[i]) P1 = -c1[i]/a0[i]+2*c2[i]*(z_local+z1[i])/a0[i]/a0[i];
352 
353  P2 = 2*c2[i]/a0[i]/a0[i];
354 
355  cte = 1 + G4Exp(c0[i]);
356 
357  K1 = -cte*P1*G4Exp(P0)/( (1+G4Exp(P0))*(1+G4Exp(P0)) );
358 
359  K2 = -cte*G4Exp(P0)*(
360  P2/( (1+G4Exp(P0))*(1+G4Exp(P0)) )
361  +2*P1*K1/(1+G4Exp(P0))/cte
362  +P1*P1/(1+G4Exp(P0))/(1+G4Exp(P0))
363  );
364 
365  K3 = -cte*G4Exp(P0)*(
366  (3*P2*P1+P1*P1*P1)/(1+G4Exp(P0))/(1+G4Exp(P0))
367  +4*K1*(P1*P1+P2)/(1+G4Exp(P0))/cte
368  +2*P1*(K1*K1/cte/cte+K2/(1+G4Exp(P0))/cte)
369  );
370 
371  G0 = gradient[i]*cte/(1+G4Exp(P0));
372  G1 = gradient[i]*K1;
373  G2 = gradient[i]*K2;
374  G3 = gradient[i]*K3;
375 
376  }
377 
378  // PROTECTION AGAINST LARGE SCATTERING
379 
380  if ( largeScattering )
381  {
382  G0=0;
383  G1=0;
384  G2=0;
385  G3=0;
386  }
387 
388  // MAGNETIC FIELD COMPUTATION FOR EACH QUADRUPOLE
389 
390  Bx_local = y_local*(G0-(1./12)*(3*x_local*x_local+y_local*y_local)*G2);
391  By_local = x_local*(G0-(1./12)*(3*y_local*y_local+x_local*x_local)*G2);
392  Bz_local = x_local*y_local*(G1-(1./12)*(x_local*x_local+y_local*y_local)*G3);
393 
394  Bx_quad = Bz_local*std::sin(lineAngle)+Bx_local*std::cos(lineAngle);
395  By_quad = By_local;
396  Bz_quad = Bz_local*std::cos(lineAngle)-Bx_local*std::sin(lineAngle);
397 
398  // TOTAL MAGNETIC FIELD
399 
400  Bx = Bx + Bx_quad ;
401  By = By + By_quad ;
402  Bz = Bz + Bz_quad ;
403 
404  } // LOOP ON QUADRUPOLES
405 
406 
407 } // END OF QUADRUPLET
408 
409  Bfield[0] = Bx;
410  Bfield[1] = By;
411  Bfield[2] = Bz;
412 
413 // *****************************************
414 // ELECTRIC FIELD CREATED BY SCANNING PLATES
415 // *****************************************
416 
417  Bfield[3] = 0;
418  Bfield[4] = 0;
419  Bfield[5] = 0;
420 
421  // POSITION OF EXIT OF LAST QUAD WHERE THE SCANNING PLATES START
422 
423  G4double electricPlateWidth1 = 5 * mm;
424  G4double electricPlateWidth2 = 5 * mm;
425  G4double electricPlateLength1 = 36 * mm;
426  G4double electricPlateLength2 = 34 * mm;
427  G4double electricPlateGap = 5 * mm;
428  G4double electricPlateSpacing1 = 3 * mm;
429  G4double electricPlateSpacing2 = 4 * mm;
430 
431  // APPLY VOLTAGE HERE IN VOLTS (no electrostatic deflection here)
432  G4double electricPlateVoltage1 = 0 * volt;
433  G4double electricPlateVoltage2 = 0 * volt;
434 
435  G4double electricFieldPlate1 = electricPlateVoltage1 / electricPlateSpacing1 ;
436  G4double electricFieldPlate2 = electricPlateVoltage2 / electricPlateSpacing2 ;
437 
438  G4double beginFirstZoneX = lineX + (8*quadHalfLength+3*quadSpacing)*std::sin(lineAngle);
439  G4double beginFirstZoneZ = lineZ + (8*quadHalfLength+3*quadSpacing)*std::cos(lineAngle);
440 
441  G4double beginSecondZoneX = lineX + (8*quadHalfLength+3*quadSpacing+electricPlateLength1+electricPlateGap)*std::sin(lineAngle);
442  G4double beginSecondZoneZ = lineZ + (8*quadHalfLength+3*quadSpacing+electricPlateLength1+electricPlateGap)*std::cos(lineAngle);
443 
444  G4double xA, zA, xB, zB, xC, zC, xD, zD;
445  G4double slope1, cte1, slope2, cte2, slope3, cte3, slope4, cte4;
446 
447  // WARNING : lineAngle < 0
448 
449  // FIRST PLATES
450 
451  xA = beginFirstZoneX + std::cos(lineAngle)*electricPlateSpacing1/2;
452  zA = beginFirstZoneZ - std::sin(lineAngle)*electricPlateSpacing1/2;
453 
454  xB = xA + std::sin(lineAngle)*electricPlateLength1;
455  zB = zA + std::cos(lineAngle)*electricPlateLength1;
456 
457  xC = xB - std::cos(lineAngle)*electricPlateSpacing1;
458  zC = zB + std::sin(lineAngle)*electricPlateSpacing1;
459 
460  xD = xC - std::sin(lineAngle)*electricPlateLength1;
461  zD = zC - std::cos(lineAngle)*electricPlateLength1;
462 
463  slope1 = (xB-xA)/(zB-zA);
464  cte1 = xA - slope1 * zA;
465 
466  slope2 = (xC-xB)/(zC-zB);
467  cte2 = xB - slope2 * zB;
468 
469  slope3 = (xD-xC)/(zD-zC);
470  cte3 = xC - slope3 * zC;
471 
472  slope4 = (xA-xD)/(zA-zD);
473  cte4 = xD - slope4 * zD;
474 
475 
476  if
477  (
478  x <= slope1 * z + cte1
479  && x >= slope3 * z + cte3
480  && x <= slope4 * z + cte4
481  && x >= slope2 * z + cte2
482  && std::abs(y)<=electricPlateWidth1/2
483  )
484 
485  {
486  Bfield[3] = electricFieldPlate1*std::cos(lineAngle);
487  Bfield[4] = 0;
488  Bfield[5] = -electricFieldPlate1*std::sin(lineAngle);
489 
490  }
491 
492  // SECOND PLATES
493 
494  xA = beginSecondZoneX + std::cos(lineAngle)*electricPlateWidth2/2;
495  zA = beginSecondZoneZ - std::sin(lineAngle)*electricPlateWidth2/2;
496 
497  xB = xA + std::sin(lineAngle)*electricPlateLength2;
498  zB = zA + std::cos(lineAngle)*electricPlateLength2;
499 
500  xC = xB - std::cos(lineAngle)*electricPlateWidth2;
501  zC = zB + std::sin(lineAngle)*electricPlateWidth2;
502 
503  xD = xC - std::sin(lineAngle)*electricPlateLength2;
504  zD = zC - std::cos(lineAngle)*electricPlateLength2;
505 
506  slope1 = (xB-xA)/(zB-zA);
507  cte1 = xA - slope1 * zA;
508 
509  slope2 = (xC-xB)/(zC-zB);
510  cte2 = xB - slope2 * zB;
511 
512  slope3 = (xD-xC)/(zD-zC);
513  cte3 = xC - slope3 * zC;
514 
515  slope4 = (xA-xD)/(zA-zD);
516  cte4 = xD - slope4 * zD;
517 
518  if
519  (
520  x <= slope1 * z + cte1
521  && x >= slope3 * z + cte3
522  && x <= slope4 * z + cte4
523  && x >= slope2 * z + cte2
524  && std::abs(y)<=electricPlateSpacing2/2
525  )
526 
527  {
528  Bfield[3] = 0;
529  Bfield[4] = electricFieldPlate2;
530  Bfield[5] = 0;
531  }
532 
533 //
534 
535 }
x
Float_t x
Definition: compare.C:6
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
micrometer
static constexpr double micrometer
Definition: G4SIunits.hh:100
mm
static constexpr double mm
Definition: G4SIunits.hh:115
y
Float_t y
Definition: compare.C:6
D
double D(double temp)
Definition: G4DNAElectronHoleRecombination.cc:91
z
Double_t z
Definition: advanced/microbeam/plot.C:277
P2
static const G4double * P2[nN]
Definition: G4ElectroNuclearCrossSection.cc:2177
P0
static const G4double * P0[nN]
Definition: G4ElectroNuclearCrossSection.cc:2175
m
static constexpr double m
Definition: G4SIunits.hh:129
G4double
double G4double
Definition: G4Types.hh:76
G4bool
bool G4bool
Definition: G4Types.hh:79
c2
TCanvas * c2
Definition: plot_hist.C:75
EMField::GetFieldValue
void GetFieldValue(const double Point[4], double *Bfield) const
Definition: EMField.cc:50
deg
static constexpr double deg
Definition: G4SIunits.hh:152
EMField::EMField
EMField()
Definition: EMField.cc:45
a0
const G4double a0
Definition: G4IonCoulombCrossSection.cc:68
G4int
int G4int
Definition: G4Types.hh:78
P1
static const G4double * P1[nN]
Definition: G4ElectroNuclearCrossSection.cc:2176
volt
static constexpr double volt
Definition: G4SIunits.hh:244
tesla
static constexpr double tesla
Definition: G4SIunits.hh:268
다음에 의해 생성됨 :   doxygen 1.8.5