source: trunk/documents/UserDoc/UsersGuides/ForApplicationDeveloper/html/Detector/electroMagneticField.html @ 1358

Last change on this file since 1358 was 1208, checked in by garnier, 15 years ago

CVS update
File size: 23.2 KB
Line 
1<!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN">
2<!--Converted with LaTeX2HTML 96.1-g (July 19, 1996) by Nikos Drakos (nikos@cbl.leeds.ac.uk), CBLU, University of Leeds -->
3<HTML>
4<HEAD>
5<TITLE>Getting started with Geant4: the Field</TITLE>
6<!-- Changed by: Katsuya Amako,    30-Nov-1998 -->
7<!-- Changed by: John APOSTOLAKIS,  7-Dec-1998 -->
8<!-- Proof read by: Joe Chuma,     29-Jun-1999 -->
9<!-- Changed by: Dennis Wright, 12-Dec-2002 -->
10<!-- Corrections for field changes: John Apostolakis, 17-Jun-2005 -->
11<!-- Corrections, changes by: John APOSTOLAKIS,   7-Jul-2005 -->
12</HEAD>
13
14<BODY>
15<P>
16
17<TABLE WIDTH="100%" >
18<TR>
19<TD>
20</A>
21<A HREF="index.html">
22<IMG SRC="../../../../resources/html/IconsGIF/Contents.gif" ALT="Contents" HEIGHT=16 WIDTH=59></A>
23<A HREF="material.html">
24<IMG SRC="../../../../resources/html/IconsGIF/Previous.gif" ALT="Previous" HEIGHT=16 WIDTH=59></A>
25<A HREF="hit.html">
26<IMG SRC="../../../../resources/html/IconsGIF/Next.gif" ALT="Next" HEIGHT=16 WIDTH=59></A>
27</TD>
28
29<TD ALIGN="Right"><FONT COLOR="#238E23"><FONT SIZE=-1>
30<B>Geant4 User's Guide</B> <BR>
31<B>For Application Developers</B> <BR>
32<B>Detector Definition and Response</B> </FONT></FONT> </TD>
33</TR>
34</TABLE>
35<P><BR>
36
37
38<CENTER><FONT COLOR="#238E23"><FONT SIZE=+3>
39<B>4.3 Electromagnetic Field</B> </FONT></FONT>
40</CENTER>
41<BR>
42<BR>
43
44<HR ALIGN="Center" SIZE="7%">
45<p>
46
47<a name="4.3.1">
48<H2>4.3.1 An Overview of Propagation in a Field</H2></a>
49
50<P>
51Geant4 is capable of describing and propagating in a variety of
52fields. Magnetic fields, electric fields and electromagnetic, uniform
53or non-uniform, can specified for a Geant4 setup. The
54propagation of tracks inside them can be performed to a
55user-defined accuracy. </P>
56
57<P>
58In order to propagate a track inside a field, the equation of motion
59of the particle in the field is integrated.  In general, this is done
60using a Runge-Kutta method for the integration of ordinary differential
61equations.  However, for specific cases where an analytical solution is
62known, it is possible to utilize this instead.  Several Runge-Kutta
63methods are available, suitable for different conditions.  In specific
64cases (such as a uniform field where the analytical solution is known)
65different solvers can also be used.  In addition, when an approximate
66analytical solution is known, it is possible to utilize it in an
67iterative manner in order to converge to the solution to the precision
68required.  This latter method is currently implemented and can be used
69particularly well for magnetic fields that are almost uniform. </P>
70
71<P>
72Once a method is chosen that calculates the track's propagation in
73a specific field, the curved path is broken up into linear chord
74segments.  These chord segments are determined so that they closely
75approximate the curved path.  The chords are then used to interrogate
76the Navigator as to whether the track has crossed a volume boundary.
77Several parameters are available to adjust the accuracy of the
78integration and the subsequent interrogation of the model geometry. </P>
79
80<P> 
81How closely the set of chords approximates a curved trajectory is governed
82by a parameter called the <I>miss distance</I> (also called the <I>chord
83distance</I> ).
84This is an upper bound for the value of the sagitta -  the distance between the
85'real' curved trajectory and the approximate linear trajectory of the chord.
86By setting this parameter, the user can
87control the precision of the volume interrogation. Every attempt has been made
88to ensure that all volume interrogations will be made to an accuracy within
89this <I>miss distance</I>. </P>
90<P>
91<img src="electroMagneticField.src/MissDistance.jpg" width="629" height="168">
92</P>
93
94<P>
95  Figure A: The curved trajectory will be approximated by chords, so that the
96  maximum estimated distance between curve and chord is less than the the
97  <i>miss distance</i>. </P>
98<P>
99  In addition to the <I>miss distance</I> there are two more parameters which
100  the user can set in order to adjust the accuracy (and performance) of
101  tracking in a field. In particular these parameters govern the accuracy of
102  the intersection with a volume boundary and the accuracy of the integration
103  of other steps. As such they play an important role for tracking. </P>
104
105<P> 
106  The <I>delta intersection</I> parameter is the accuracy to which an
107  intersection with a volume boundary is calculated. If a candidate boundary
108  intersection is estimated to have a precision better than this, it is
109  accepted. This parameter is especially important because it is used to limit
110  a bias that our algorithm (for boundary crossing in a field) exhibits. This
111  algorithm calculates the intersection with a volume boundary using a chord
112  between two points on the curved particle trajectory. As such, the
113  intersection point is always on the 'inside' of the curve. By setting a
114  value for this parameter that is much smaller than some acceptable error,
115  the user can limit the effect of this bias on, for example, the future
116  estimation of the reconstructed particle momentum. </P>
117<P>
118<img src="electroMagneticField.src/IntersectionError.jpg" width="237" height="498"></P>
119
120<P>
121 Figure B: The distance between the calculated chord intersection point C and
122 a computed curve point D is used to determine whether C is an accurate
123 representation of the intersection of the curved path ADB with a volume
124 boundary. Here CD is likely too large, and a new intersection on the chord AD
125 will be calculated.</P>
126<P>
127 The <I>delta one step</I> parameter is the accuracy for the endpoint of
128 'ordinary' integration steps, those which do not intersect a volume boundary.
129 This parameter is a limit on the estimated error of the endpoint of each
130 physics step. It can be seen as akin to a statistical uncertainty and is not
131 expected to contribute any systematic behavior to physical quantities. In
132 contrast, the bias addressed by <I>delta intersection</I> is clearly
133 correlated with potential systematic errors in the momentum of reconstructed
134 tracks. Thus very strict limits on the intersection parameter should be used
135 in tracking detectors or wherever the intersections are used to reconstruct a
136 track's momentum. </P>
137
138<P>
139<I>Delta intersection</I> and <I>delta one step</I> are parameters of
140the Field Manager; the user can set them according to the demands of
141his application.  Because it is possible to use more than one field
142manager, different values can be set for different detector regions. </P>
143
144<P>
145Note that reasonable values for the two parameters are strongly
146coupled: it does not make sense to request an accuracy of 1 nm for
147<I>delta intersection</I> and accept 100 &#956m for the
148<I>delta one step</I> error value.  Nevertheless
149<I>delta intersection</I> is the more important of the two.  It is
150recommended that these parameters should not differ significantly -
151certainly not by more than an order of magnitude. </P>
152
153<hr>
154<a name="4.3.2">
155<H2>4.3.2 Practical Aspects</H2></a>
156
157<B>Creating a Magnetic Field for a Detector</B>
158<P>
159The simplest way to define a field for a detector involves the
160following steps:
161 <P>
162 <ol>
163 <li>create a field:
164 <PRE>
165    G4UniformMagField* magField
166      = new G4UniformMagField(G4ThreeVector(0.,0.,fieldValue));
167 </PRE>
168 <li>set it as the default field:
169 <PRE>
170    G4FieldManager* fieldMgr
171      = G4TransportationManager::GetTransportationManager()
172        ->GetFieldManager();
173       fieldMgr->SetDetectorField(magField);
174 </PRE>
175 <li>create the objects which calculate the trajectory:
176 <PRE>
177    fieldMgr->CreateChordFinder(magField);
178 </PRE>
179 </ol>
180 <P>
181 To change the accuracy of volume intersection use the <tt>SetDeltaChord</tt>
182 method:
183 <PRE>
184    fieldMgr->GetChordFinder()->SetDeltaChord( G4double newValue);
185 </PRE>
186
187<P>
188<B>Creating a Field for a Part of the Volume Hierarchy</B> 
189<P> It is possible to create a field for a part of the detector. In particular
190  it can describe the field (with pointer fEmField, for example) inside a
191  logical volume and all its daughters. This can be done by simply creating a
192  <tt>G4FieldManager</tt> and attaching it to a logical volume (with pointer,
193  logicVolumeWithField, for example) or set of logical volumes.
194 <PRE>
195  G4bool allLocal = true;       
196  logicVolumeWithField->SetFieldManager(fEmField->GetLocalFieldManager(),allLocal);
197 </PRE>
198</P>
199<P>&nbsp;</P>
200<P>You can then choose whether the logical volume's daughter volumes
201  will also be given this new field.  If so, the field will be 'pushed' into
202  its daughters, unless they already have a field manager.</P>
203
204<B>Creating an Electric or Electromagnetic Field</B>
205<P> 
206 The design and implementation of the <i>Field</i> category allows and enables
207 the use of an electric or combined electromagnetic field.
208 These fields can also vary with time, as can magnetic fields.
209
210 <P>
211 Source listing 4.3.1 shows how to define a uniform electric field for the
212 whole of a detector.
213<p>
214<center>
215<table border=2 cellpadding=10>
216<tr>
217<td>
218 <PRE>
219 // in the header file (or first)
220  #include "G4EqMagElectricField.hh"
221  #include "G4UniformElectricField.hh"
222
223  ...
224    G4ElectricField*        fEMfield;
225    G4EqMagElectricField*   fEquation;
226    G4MagIntegratorStepper* fStepper;
227    G4FieldManager*         fFieldMgr;
228    G4double                fMinStep ;
229    G4ChordFinder*          fChordFinder ;
230
231 // in the source file
232
233  {
234    fEMfield = new G4UniformElectricField(
235                 G4ThreeVector(0.0,100000.0*kilovolt/cm,0.0));
236
237    // Create an equation of motion for this field
238    fEquation = new G4EqMagElectricField(fEMfield);
239
240    G4int nvar = 8;
241    fStepper = new G4ClassicalRK4( fEquation, nvar );       
242
243    // Get the global field manager
244    fFieldManager= G4TransportationManager::GetTransportationManager()->
245         GetFieldManager();
246    // Set this field to the global field manager
247    fFieldManager->SetDetectorField(fEMfield );
248
249    fMinStep     = 0.010*mm ; // minimal step of 10 microns
250
251    fIntgrDriver = new G4MagInt_Driver(fMinStep,
252                                     fStepper,
253                                     fStepper->GetNumberOfVariables() );
254
255    fChordFinder = new G4ChordFinder(fIntgrDriver);
256    fFieldManager->SetChordFinder( fChordFinder );
257
258  }
259</PRE>
260</td>
261</tr>
262<tr>
263<td align=center>
264 Source listing 4.3.1<BR>
265 How to define a uniform electric field for the whole of a detector, extracted
266 from example in examples/extended/field/field02 .
267</td>
268</tr>
269</table></center>
270
271
272<P> An example with an electric field is examples/extended/field/field02, where
273the class F02ElectricFieldSetup demonstrates how to set these and other
274parameters, and how to choose different Integration Steppers. </P>
275
276The user can also create their own type of field, inheriting from
277<tt>G4VField</tt>,
278and an associated Equation of Motion class (inheriting from <tt>G4EqRhs</tt>)
279to simulate other types of fields.
280
281<P>
282<B>Choosing a Stepper</B></P>
283<P>
284Runge-Kutta integration is used to compute the motion of a charged
285track in a general field.  There are many general steppers from which
286to choose, of low and high order, and specialized steppers for pure
287magnetic fields.  By default, Geant4 uses the classical fourth-order
288Runge-Kutta stepper, which is general purpose and robust.  If the field
289is known to have specific properties, lower or higher order steppers
290can be used to obtain the same quality results using fewer computing
291cycles. </P>
292
293<P>
294In particular, if the field is calculated from a field map, a lower
295order stepper is recommended.  The less smooth the field is, the lower
296the order of the stepper that should be used.  The choice of lower
297order steppers includes the third order stepper <tt>G4SimpleHeum</tt>, the
298second order <tt>G4ImplicitEuler</tt> and <tt>G4SimpleRunge</tt>,
299and the first order
300<tt>G4ExplicitEuler</tt>.  A first order stepper would be useful only for
301very rough fields.  For somewhat smooth fields (intermediate), the
302choice between second and third order steppers should be made by
303trial and error.  Trying a few different types of steppers for a particular
304field or application is suggested if maximum performance is a goal. </P>
305
306The choice of stepper depends on the type of field: magnetic or general.
307A general field can be an electric or electromagnetic field,
308it can be a magnetic field or a user-defined field (which requires
309a user-defined equation of motion.)
310
311For a general field several steppers are available as alternatives to the
312default (<tt>G4ClassicalRK4</tt>):
313<PRE>
314   G4int nvar = 8;  // To integrate time & energy
315                    //    in addition to position, momentum
316   G4EqMagElectricField* fEquation= new G4EqMagElectricField(fEMfield);
317
318   fStepper = new G4SimpleHeum( fEquation, nvar );         
319            //  3rd  order, a good alternative to ClassicalRK
320   fStepper = new G4SimpleRunge( fEquation, nvar );
321            //  2nd  order, for less smooth fields
322   fStepper = new G4CashKarpRKF45( fEquation );     
323            // 4/5th order for very smooth fields
324</PRE>
325
326<P>
327Specialized steppers for pure magnetic fields are also available.
328They take into account the fact that a local trajectory in a slowly varying
329field will not vary significantly from a helix.  Combining this in
330with a variation the Runge-Kutta method can provide higher accuracy at lower
331computational cost when large steps are possible. </P>
332
333<PRE>
334   G4Mag_UsualEqRhs*
335      fEquation = new G4Mag_UsualEqRhs(fMagneticField);
336   fStepper = new G4HelixImplicitEuler( fEquation );
337      // Note that for magnetic field that do not vary with time,
338      //  the default number of variables suffices.
339 
340      // or ..
341     fStepper = new G4HelixExplicitEuler( fEquation );
342     fStepper = new G4HelixSimpleRunge( fEquation );   
343</PRE>
344
345<P> 
346You can choose an alternative stepper either when the field manager
347is constructed or later. 
348At the construction of the ChordFinder it is an optional argument:
349
350<PRE>
351 G4ChordFinder( G4MagneticField* itsMagField,
352                G4double         stepMinimum = 1.0e-2 * mm,
353                G4MagIntegratorStepper* pItsStepper = 0 );
354</PRE>
355
356To change the stepper at a later time use
357<PRE>
358 pChordFinder->GetIntegrationDriver()
359             ->RenewStepperAndAdjust( newStepper );
360</PRE>
361<BR>
362
363<p>
364<B>How to Adjust the Accuracy of Propagation</B></p>
365<p>In order to obtain a particular accuracy in tracking particles through an
366 electromagnetic field, it is necessary to adjust the parameters of the field
367 propagation module. In the following section, some of these additional
368 parameters are discussed.</p>
369<p>
370 When integration is used to calculate the trajectory, it is necessary to
371 determine an acceptable level of numerical imprecision in order to get
372 performant simulation with acceptable errors. The parameters in Geant4 tell
373 the field module what level of integration inaccuracy is acceptable.</p>
374<p> In all quantities which are integrated (position, momentum, energy) there
375  will be errors. Here, however, we focus on the error in two key quantities:
376  the position and the momentum. (The error in the energy will come from the
377  momentum integration). </p>
378<p> 
379 Three parameters exist which are relevant to the integration accuracy.
380 DeltaOneStep is a distance and is roughly the position error which is
381 acceptable in an integration step. Since many integration steps may be
382 required for a single physics step, DeltaOneStep should be a fraction of the
383 average physics step size. The next two parameters impose a further limit on
384 the relative error of the position/momentum inaccuracy. EpsilonMin and
385 EpsilonMax impose a minimum and maximum on this relative error -
386 and take precedence over DeltaOneStep. (Note: if you set
387 EpsilonMin=EpsilonMax=your-value, then all steps will be made to this
388 relative precision.</p>
389
390<center>
391<table border=2 cellpadding=10>
392<tr>
393<td>
394<PRE>
395  G4FieldManager *globalFieldManager;
396
397  G4TransportationManager *transportMgr=
398       G4TransportationManager::GetTransportationManager();
399
400  globalFieldManager = transportMgr->GetFieldManager();
401                            // Relative accuracy values:
402  G4double minEps= 1.0e-5;  //   Minimum & value for smallest steps
403  G4double maxEps= 1.0e-4;  //   Maximum & value for largest steps
404
405  globalFieldManager->SetMinimumEpsilonStep( minEps );
406  globalFieldManager->SetMaximumEpsilonStep( maxEps );
407  globalFieldManager->SetDeltaOneStep( 0.5e-3 * mm );  // 0.5 micrometer
408
409  G4cout << "EpsilonStep: set min= " << minEps << " max= " << maxEps << G4endl;
410</PRE>
411</td>
412</tr>
413<tr>
414<td align=center>
415 Source listing 4.3.2<BR>
416 How to set accuracy parameters for the 'global' field of the setup.
417</td>
418</tr>
419</table></center>
420<P>
421
422<p>
423 We note that the relevant parameters above limit the inaccuracy in each step.
424 The final inaccuracy due to the full trajectory will accumulate!</p>
425<p> 
426 The exact point at which a track crosses a boundary is also calculated with
427 finite accuracy. To limit this inaccuracy, a parameter called
428 DeltaIntersection is used. This is a maximum for the inaccuracy of a single
429 boundary crossing.  Thus the accuracy of the position of the track after a
430 number of boundary crossings is directly proportional to the number of
431 boundaries.</p>
432
433<p> 
434<b>Choosing different accuracies for the same volume</b></p>
435<p>
436 It is possible to create a FieldManager which has different properties for
437 particles of different momenta (or depending on other parameters of a track).
438 This is useful, for example, in obtaining high accuracy for 'important' tracks
439 (e.g. muons) and accept less accuracy in tracking others (e.g. electrons). To
440 use this, you must create your own field manager which uses the method</p>
441<p> 
442 void ConfigureForTrack( const G4Track * );</p>
443<p>
444 to configure itself using the parameters of the current track. An example of
445 this will be available in examples/extended/field05.</p>
446<p>
447<b>Parameters that must scale with problem size</b></p>
448<p>
449 The default settings of this module are for problems with the physical size
450 of a typical high energy physics setup, that is, distances smaller than about
451 one kilometer. A few parameters are necessary to carry this information to
452 the magnetic field module, and must typically be rescaled for problems of
453 vastly different sizes in order to get reasonable performance and
454 robustness.  Two of these parameters are the maximum acceptable step and
455 the minimum step size.</p>
456
457<p>
458 The <b>maximum acceptable step</b> should be set to a distance larger than
459 the biggest reasonable step. If the apparatus in a setup has a diameter of
460 two meters, a likely maximum acceptable steplength would be 10 meters. A
461 particle could then take large spiral steps, but would not attempt to take,
462 for example, a 1000-meter-long step in the case of a very low-density 
463 material. Similarly, for problems of a planetary scale, such as the earth
464 with its radius of roughly 6400 km, a maximum acceptabe steplength of a few
465 times this value would be reasonable.</p>
466<p>
467 An upper limit for the size of a step
468 is a parameter of <tt>G4PropagatorInField</tt>, and can be
469 set by calling its <tt>SetLargestAcceptableStep</tt> method.</p>
470
471<p>
472 The <b>minimum step size</b> is used during integration to limit the
473 amount of work in difficult cases. It is possible that strong fields or
474 integration problems can force the integrator to try very small steps;
475 this parameter stops them from becoming unnecessarily small.</p>
476<p>
477 Trial steps smaller than this parameter will be treated with less accuracy,
478 and may even be ignored, depending on the situation.</p>
479<p>
480 The minimum step size is a parameter of the MagInt_Driver, but can be set in
481 the contstructor of G4ChordFinder, as in the source listing above.</p>
482
483
484<p><B>Known Issues</B> </p>
485<P>
486Currently it is computationally expensive to change the <I>miss distance</I>
487to very small values, as it causes tracks to be limited to curved sections
488whose 'bend' is smaller than this value. (The bend is the distance of the
489mid-point from the chord between endpoints.) For tracks with small curvature
490(typically low momentum particles in strong fields) this can cause a large
491number of steps
492<ul>
493<li>even in areas where there are no volumes to intersect
494   (something that is expected to be addressed in future development,
495    in which the safety will be utilized to partially alleviate this
496    limitation)
497
498<li>especially in a region near a volume boundary 
499    (in which case it is necessary in order to discover
500     whether a track might intersect a volume for only a short distance.)
501</ul>
502Requiring such precision at the intersection is clearly expensive, and new
503development would be necessary to minimize the expense.
504
505<P>
506By contrast, changing the intersection parameter is less computationally
507expensive. It causes further calculation for only a fraction of the steps,
508in particular those that intersect a volume boundary. </P>
509
510<BR>
511
512<hr>
513<a name="4.3.3">
514<H2>4.3.3 Spin Tracking</H2></a>
515
516The effects of a particle's motion on the precession of its spin angular
517momentum in slowly varying external fields are simulated. The
518relativistic equation of motion for spin is known as the BMT equation.
519The equation demonstrates a
520remarkable property; in a purely magnetic field, in vacuum, and
521neglecting small anomalous magnetic moments, the particle's spin
522precesses in such a manner that the longitudinal polarization remains a
523constant, whatever the motion of the particle. But when the particle
524interacts with electric fields of the medium and multiple scatters, the
525spin, which is related to the particle's magnetic moment, does not
526participate, and the need thus arises to propagate it independent of the
527momentum vector. In the case of a polarized muon beam, for example, it
528is important to predict the muon's spin direction at decay-time in
529order to simulate the decay electron (Michel) distribution correctly.
530<P>
531In order to track the spin of a particle in a magnetic field, you need to code
532the following:
533
534<P>
535<ol>
536<li> in your DetectorConstruction
537<PRE>
538#include "G4Mag_SpinEqRhs.hh"
539
540  G4Mag_EqRhs* fEquation = new G4Mag_SpinEqRhs(magField);
541
542  G4MagIntegratorStepper* pStepper = new G4ClassicalRK4(fEquation,12);
543                                                       <B>notice the 12 </B>
544</PRE>
545<P>
546<li> in your PrimaryGenerator
547<PRE>
548  particleGun->SetParticlePolarization(G4ThreeVector p)
549</PRE>
550
551for example:
552
553<PRE>
554  particleGun->
555  SetParticlePolarization(-(particleGun->GetParticleMomentumDirection()));
556</PRE>
557
558or
559
560<PRE>
561  particleGun->
562  SetParticlePolarization(particleGun->GetParticleMomentumDirection()
563                                         .cross(G4ThreeVector(0.,1.,0.)));
564</PRE>
565
566where you set the initial spin direction.
567
568</ol>
569</P>
570
571While the G4Mag_SpinEqRhs class constructor
572
573<PRE>
574  G4Mag_SpinEqRhs::G4Mag_SpinEqRhs( G4MagneticField* MagField )
575     : G4Mag_EqRhs( MagField )
576  {
577      anomaly = 1.165923e-3;
578  }
579</PRE>
580
581sets the muon anomaly by default, the class also comes with the public method:
582
583<PRE>
584  inline void SetAnomaly(G4double a) { anomaly = a; }
585</PRE>
586
587with which you can set the magnetic anomaly to any value you require.
588
589<P>
590For the moment, the code is written such that field tracking of the spin is
591done only for particles with non-zero charge. Please, see the Forum posting:
592
593<PRE>
594http://geant4-hn.slac.stanford.edu:5090/HyperNews/public/get/emfields/88/3/1.html
595</PRE>
596
597for modifications the user is required to make to facilitate neutron spin
598tracking.
599</P>
600
601<BR>
602
603<HR><A HREF="../../../../Authors/html/subjectsToAuthors.html">
604<I>About the authors</i></A>   <p></P>
605</BODY>
606</HTML>
607
608
Note: See TracBrowser for help on using the repository browser.