ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EMField.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 }