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

Last change on this file since 1282 was 1230, checked in by garnier, 14 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.