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