source: trunk/examples/advanced/microbeam/src/MicrobeamEMField.cc@ 1345

Last change on this file since 1345 was 1230, checked in by garnier, 16 years ago

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