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