<!-- ******************************************************** -->
<!--                                                          -->
<!--  [History]                                               -->
<!--    Changed by: Katsuya Amako,    30-Nov-1998             -->
<!--    Converted with LaTeX2HTML 96.1-g (July 19, 1996) by   -->
<!--      Nikos Drakos (nikos@cbl.leeds.ac.uk),               -->
<!--      CBLU, University of Leeds                           -->
<!--    Changed by: John APOSTOLAKIS,  7-Dec-1998             -->
<!--    Proof read by: Joe Chuma,     29-Jun-1999             -->
<!--    Changed by: Dennis Wright, 12-Dec-2002                -->
<!--    Corrections for field changes: John Apostolakis,      -->
<!--      17-Jun-2005                                         -->
<!--    Corrections, changes by: John APOSTOLAKIS, 7-Jul-2005 -->
<!--    Converted to DocBook: Katsuya Amako, Aug-2006         -->
<!--                                                          -->
<!-- ******************************************************** -->


<!-- ******************* Section (Level#1) ****************** -->
<sect1 id="sect.EMField">
<title>
Electromagnetic Field
</title>


<!-- ******************* Section (Level#2) ****************** -->
<sect2 id="sect.EMField.Overview">
<title>
An Overview of Propagation in a Field
</title>

<para>
Geant4 is capable of describing and propagating in a variety of
fields. Magnetic fields, electric fields and electromagnetic,
uniform or non-uniform, can specified for a Geant4 setup. The
propagation of tracks inside them can be performed to a
user-defined accuracy.
</para>

<para>
In order to propagate a track inside a field, the equation of
motion of the particle in the field is integrated. In general, this
is done using a Runge-Kutta method for the integration of ordinary
differential equations. However, for specific cases where an
analytical solution is known, it is possible to utilize this
instead. Several Runge-Kutta methods are available, suitable for
different conditions. In specific cases (such as a uniform field
where the analytical solution is known) different solvers can also
be used. In addition, when an approximate analytical solution is
known, it is possible to utilize it in an iterative manner in order
to converge to the solution to the precision required. This latter
method is currently implemented and can be used particularly well
for magnetic fields that are almost uniform.
</para>

<para>
Once a method is chosen that calculates the track's propagation
in a specific field, the curved path is broken up into linear chord
segments. These chord segments are determined so that they closely
approximate the curved path. The chords are then used to
interrogate the Navigator as to whether the track has crossed a
volume boundary. Several parameters are available to adjust the
accuracy of the integration and the subsequent interrogation of the
model geometry.
</para>

<para>
How closely the set of chords approximates a curved trajectory
is governed by a parameter called the <emphasis>miss distance</emphasis> (also
called the <emphasis>chord distance</emphasis> ). This is an upper bound for the
value of the sagitta - the distance between the 'real' curved
trajectory and the approximate linear trajectory of the chord. By
setting this parameter, the user can control the precision of the
volume interrogation. Every attempt has been made to ensure that
all volume interrogations will be made to an accuracy within this
<emphasis>miss distance</emphasis>.

<figure id="fig.EMField_1">
<title>
  The curved trajectory will be approximated by chords,
  so that the maximum estimated distance between curve and chord is
  less than the the <emphasis>miss distance</emphasis>.
</title>
<mediaobject>
  <imageobject role="fo">
    <imagedata fileref="./AllResources/Detector/electroMagneticField.src/MissDistance.jpg" 
               format="JPG" contentwidth="12.0cm" align="center" />
  </imageobject>
  <imageobject role="html">
    <imagedata fileref="./AllResources/Detector/electroMagneticField.src/MissDistance.jpg"
               format="JPG" align="center" />
  </imageobject>
  <textobject>
    <phrase>Miss Distance</phrase>
  </textobject>
</mediaobject>
</figure>
</para>

<para>
In addition to the <emphasis>miss distance</emphasis> there are two more
parameters which the user can set in order to adjust the accuracy
(and performance) of tracking in a field. In particular these
parameters govern the accuracy of the intersection with a volume
boundary and the accuracy of the integration of other steps. As
such they play an important role for tracking.
</para>

<para>
The <emphasis>delta intersection</emphasis> parameter is the accuracy to which
an intersection with a volume boundary is calculated. If a
candidate boundary intersection is estimated to have a precision
better than this, it is accepted. This parameter is especially
important because it is used to limit a bias that our algorithm
(for boundary crossing in a field) exhibits. This algorithm
calculates the intersection with a volume boundary using a chord
between two points on the curved particle trajectory. As such, the
intersection point is always on the 'inside' of the curve. By
setting a value for this parameter that is much smaller than some
acceptable error, the user can limit the effect of this bias on,
for example, the future estimation of the reconstructed particle
momentum.

<figure id="fig.EMField_2">
<title>
  The distance between the calculated chord intersection
  point C and a computed curve point D is used to determine whether C
  is an accurate representation of the intersection of the curved
  path ADB with a volume boundary. Here CD is likely too large, and a
  new intersection on the chord AD will be calculated.
</title>
<mediaobject>
  <imageobject role="fo">
    <imagedata fileref="./AllResources/Detector/electroMagneticField.src/IntersectionError.jpg" 
               format="JPG" contentdepth="5.0cm" align="center" />
  </imageobject>
  <imageobject role="html">
    <imagedata fileref="./AllResources/Detector/electroMagneticField.src/IntersectionError.jpg" 
               format="JPG" align="center" />
  </imageobject>
  <textobject>
    <phrase>Figure: The distance between the calculated chord intersection
      point and a computed curve point.</phrase>
  </textobject>
</mediaobject>
</figure>
</para>

<para>
The <emphasis>delta one step</emphasis> parameter is the accuracy for the
endpoint of 'ordinary' integration steps, those which do not
intersect a volume boundary. This parameter is a limit on the
estimated error of the endpoint of each physics step. It can be
seen as akin to a statistical uncertainty and is not expected to
contribute any systematic behavior to physical quantities. In
contrast, the bias addressed by <emphasis>delta intersection</emphasis> is
clearly correlated with potential systematic errors in the momentum
of reconstructed tracks. Thus very strict limits on the
intersection parameter should be used in tracking detectors or
wherever the intersections are used to reconstruct a track's
momentum.
</para>

<para>
<emphasis>Delta intersection</emphasis> and <emphasis>delta one step</emphasis> are
parameters of the Field Manager; the user can set them according to
the demands of his application. Because it is possible to use more
than one field manager, different values can be set for different
detector regions.
</para>

<para>
Note that reasonable values for the two parameters are strongly
coupled: it does not make sense to request an accuracy of 1 nm for
<emphasis>delta intersection</emphasis> and accept 100 &amp;#956m for the
<emphasis>delta one step</emphasis> error value. Nevertheless <emphasis>delta
intersection</emphasis> is the more important of the two. It is
recommended that these parameters should not differ significantly -
certainly not by more than an order of magnitude.
</para>

</sect2>


<!-- ******************* Section (Level#2) ****************** -->
<sect2 id="sect.EMField.Pract">
<title>
Practical Aspects
</title>

<!-- ******************* Section (Level#3) ****************** -->
<sect3 id="sect.EMField.Pract.MagField">
<title>
Creating a Magnetic Field for a Detector
</title>

<para>
The simplest way to define a field for a detector involves the
following steps:

<orderedlist spacing="compact">
  <listitem><para>
    create a field:
    <informalexample>
    <programlisting>
    G4UniformMagField* magField
      = new G4UniformMagField(G4ThreeVector(0.,0.,fieldValue));
    </programlisting>
    </informalexample>
  </para></listitem>
  <listitem><para>
    set it as the default field:
    <informalexample>
    <programlisting>
    G4FieldManager* fieldMgr
      = G4TransportationManager::GetTransportationManager()
        -&gt;GetFieldManager();
       fieldMgr-&gt;SetDetectorField(magField);
    </programlisting>
    </informalexample>
  </para></listitem>
  <listitem><para>
    create the objects which calculate the trajectory:
    <informalexample>
    <programlisting>
    fieldMgr-&gt;CreateChordFinder(magField);
    </programlisting>
    </informalexample>
</para></listitem>
</orderedlist>
</para>

<para>
To change the accuracy of volume intersection use the
<literal>SetDeltaChord</literal> method:

<informalexample>
<programlisting>
    fieldMgr-&gt;GetChordFinder()-&gt;SetDeltaChord( G4double newValue);
</programlisting>
</informalexample>
</para>

</sect3>

<!-- ******************* Section (Level#3) ****************** -->
<sect3 id="sect.EMField.Pract.PartVol">
<title>
Creating a Field for a Part of the Volume Hierarchy
</title>

<para>
It is possible to create a field for a part of the detector. In
particular it can describe the field (with pointer fEmField, for
example) inside a logical volume and all its daughters. This can be
done by simply creating a <literal>G4FieldManager</literal> and attaching it
to a logical volume (with pointer, logicVolumeWithField, for
example) or set of logical volumes.

<informalexample>
<programlisting>
  G4bool allLocal = true;       
  logicVolumeWithField-&gt;SetFieldManager(localFieldManager, allLocal);
</programlisting>
</informalexample>
</para>

<para>
Using the second parameter to <literal>SetFieldManager</literal> you choose
whether  daughter volumes of this logical volume will also be given this new
field.  If it has the value <literal>true</literal>, the field will be
assigned also to its daughters, and all their sub-volumes.  
Else, if it is <literal>false</literal>,  it will be copied only to those
daughter volumes which do not have a field manager already.
</para>

</sect3>

<!-- ******************* Section (Level#3) ****************** -->
<sect3 id="sect.EMField.Pract.Electric">
<title>
Creating an Electric or Electromagnetic Field
</title>

<para>
The design and implementation of the <emphasis>Field</emphasis> category
allows and enables the use of an electric or combined
electromagnetic field. These fields can also vary with time, as can
magnetic fields.
</para>

<para>
Source listing <xref linkend="programlist_EMField_1" /> shows how to 
define a uniform electric field for the whole of a detector.

<example id="programlist_EMField_1">
<title>
How to define a uniform electric field for the whole of a detector, 
extracted from example in examples/extended/field/field02 .
</title>

<programlisting>
 // in the header file (or first)
  #include "G4EqMagElectricField.hh"
  #include "G4UniformElectricField.hh"

  ...
    G4ElectricField*        fEMfield;
    G4EqMagElectricField*   fEquation;
    G4MagIntegratorStepper* fStepper;
    G4FieldManager*         fFieldMgr;
    G4double                fMinStep ;
    G4ChordFinder*          fChordFinder ;

 // in the source file

  {
    fEMfield = new G4UniformElectricField(
                 G4ThreeVector(0.0,100000.0*kilovolt/cm,0.0));

    // Create an equation of motion for this field
    fEquation = new G4EqMagElectricField(fEMfield); 

    G4int nvar = 8;
    fStepper = new G4ClassicalRK4( fEquation, nvar );       

    // Get the global field manager 
    fFieldManager= G4TransportationManager::GetTransportationManager()-&gt;
         GetFieldManager();
    // Set this field to the global field manager 
    fFieldManager-&gt;SetDetectorField(fEMfield );

    fMinStep     = 0.010*mm ; // minimal step of 10 microns

    fIntgrDriver = new G4MagInt_Driver(fMinStep, 
                                     fStepper, 
                                     fStepper-&gt;GetNumberOfVariables() );

    fChordFinder = new G4ChordFinder(fIntgrDriver);
    fFieldManager-&gt;SetChordFinder( fChordFinder );

  }
</programlisting>
</example>
</para>

<para>
An example with an electric field is
examples/extended/field/field02, where the class
F02ElectricFieldSetup demonstrates how to set these and other
parameters, and how to choose different Integration Steppers.
</para>

<para>
The user can also create their own type of field, inheriting from
<literal>G4VField</literal>, and an associated Equation of Motion class
(inheriting from <literal>G4EqRhs</literal>) to simulate other types of
fields.
</para>

</sect3>

<!-- ******************* Section (Level#3) ****************** -->
<sect3 id="sect.EMField.Pract.Stepper">
<title>
Choosing a Stepper
</title>

<para>
Runge-Kutta integration is used to compute the motion of a
charged track in a general field. There are many general steppers
from which to choose, of low and high order, and specialized
steppers for pure magnetic fields. By default, Geant4 uses the
classical fourth-order Runge-Kutta stepper, which is general
purpose and robust. If the field is known to have specific
properties, lower or higher order steppers can be used to obtain
the same quality results using fewer computing cycles.
</para>

<para>
In particular, if the field is calculated from a field map, a
lower order stepper is recommended. The less smooth the field is,
the lower the order of the stepper that should be used. The choice
of lower order steppers includes the third order stepper
<literal>G4SimpleHeum</literal>, the second order <literal>G4ImplicitEuler</literal>
and <literal>G4SimpleRunge</literal>, and the first order
<literal>G4ExplicitEuler</literal>. A first order stepper would be useful
only for very rough fields. For somewhat smooth fields
(intermediate), the choice between second and third order steppers
should be made by trial and error. Trying a few different types of
steppers for a particular field or application is suggested if
maximum performance is a goal.
</para>

<para>
The choice of stepper depends on the type of field: magnetic or
general. A general field can be an electric or electromagnetic
field, it can be a magnetic field or a user-defined field (which
requires a user-defined equation of motion.) For a general field
several steppers are available as alternatives to the default
(<literal>G4ClassicalRK4</literal>):

<informalexample>
<programlisting>
   G4int nvar = 8;  // To integrate time &amp; energy 
                    //    in addition to position, momentum
   G4EqMagElectricField* fEquation= new G4EqMagElectricField(fEMfield); 

   fStepper = new G4SimpleHeum( fEquation, nvar );         
            //  3rd  order, a good alternative to ClassicalRK
   fStepper = new G4SimpleRunge( fEquation, nvar ); 
            //  2nd  order, for less smooth fields
   fStepper = new G4CashKarpRKF45( fEquation );     
            // 4/5th order for very smooth fields 
</programlisting>
</informalexample>
</para>

<para>
Specialized steppers for pure magnetic fields are also
available. They take into account the fact that a local trajectory
in a slowly varying field will not vary significantly from a helix.
Combining this in with a variation the Runge-Kutta method can
provide higher accuracy at lower computational cost when large
steps are possible.

<informalexample>
<programlisting>
   G4Mag_UsualEqRhs* 
      fEquation = new G4Mag_UsualEqRhs(fMagneticField); 
   fStepper = new G4HelixImplicitEuler( fEquation ); 
      // Note that for magnetic field that do not vary with time,
      //  the default number of variables suffices.
 
      // or ..
     fStepper = new G4HelixExplicitEuler( fEquation ); 
     fStepper = new G4HelixSimpleRunge( fEquation );   
</programlisting>
</informalexample>
</para>

<para>
You can choose an alternative stepper either when the field
manager is constructed or later. At the construction of the
ChordFinder it is an optional argument:

<informalexample>
<programlisting>
 G4ChordFinder( G4MagneticField* itsMagField,
                G4double         stepMinimum = 1.0e-2 * mm,
                G4MagIntegratorStepper* pItsStepper = 0 );
</programlisting>
</informalexample>
</para>

<para>
To change the stepper at a later time use

<informalexample>
<programlisting>
 pChordFinder-&gt;GetIntegrationDriver()
             -&gt;RenewStepperAndAdjust( newStepper );
</programlisting>
</informalexample>
</para>

</sect3>

<!-- ******************* Section (Level#3) ****************** -->
<sect3 id="sect.EMField.Pract.Adjust">
<title>
How to Adjust the Accuracy of Propagation
</title>

<para>
In order to obtain a particular accuracy in tracking particles
through an electromagnetic field, it is necessary to adjust the
parameters of the field propagation module. In the following
section, some of these additional parameters are discussed.
</para>

<para>
When integration is used to calculate the trajectory, it is
necessary to determine an acceptable level of numerical imprecision
in order to get performant simulation with acceptable errors. The
parameters in Geant4 tell the field module what level of
integration inaccuracy is acceptable.
</para>

<para>
In all quantities which are integrated (position, momentum,
energy) there will be errors. Here, however, we focus on the error
in two key quantities: the position and the momentum. (The error in
the energy will come from the momentum integration).
</para>

<para>
Three parameters exist which are relevant to the integration
accuracy. DeltaOneStep is a distance and is roughly the position
error which is acceptable in an integration step. Since many
integration steps may be required for a single physics step,
DeltaOneStep should be a fraction of the average physics step size.
The next two parameters impose a further limit on the relative
error of the position/momentum inaccuracy. EpsilonMin and
EpsilonMax impose a minimum and maximum on this relative error -
and take precedence over DeltaOneStep. (Note: if you set
EpsilonMin=EpsilonMax=your-value, then all steps will be made to
this relative precision.

<example id="programlist_EMField_2">
<title>
How to set accuracy parameters for the 'global' field of the setup.
</title>

<programlisting>
  G4FieldManager *globalFieldManager;

  G4TransportationManager *transportMgr= 
       G4TransportationManager::GetTransportationManager();

  globalFieldManager = transportMgr-&gt;GetFieldManager();
                            // Relative accuracy values:
  G4double minEps= 1.0e-5;  //   Minimum &amp; value for smallest steps
  G4double maxEps= 1.0e-4;  //   Maximum &amp; value for largest steps

  globalFieldManager-&gt;SetMinimumEpsilonStep( minEps );
  globalFieldManager-&gt;SetMaximumEpsilonStep( maxEps );
  globalFieldManager-&gt;SetDeltaOneStep( 0.5e-3 * mm );  // 0.5 micrometer

  G4cout &lt;&lt; "EpsilonStep: set min= " &lt;&lt; minEps &lt;&lt; " max= " &lt;&lt; maxEps &lt;&lt; G4endl;
</programlisting>
</example>
</para>

<para>
We note that the relevant parameters above limit the inaccuracy
in each step. The final inaccuracy due to the full trajectory will
accumulate!
</para>

<para>
The exact point at which a track crosses a boundary is also
calculated with finite accuracy. To limit this inaccuracy, a
parameter called DeltaIntersection is used. This is a maximum for
the inaccuracy of a single boundary crossing. Thus the accuracy of
the position of the track after a number of boundary crossings is
directly proportional to the number of boundaries.
</para>

</sect3>

<!-- ******************* Section (Level#3) ****************** -->
<sect3 id="sect.EMField.Pract.DiffAcc">
<title>
Choosing different accuracies for the same volume
</title>

<para>
It is possible to create a FieldManager which has different
properties for particles of different momenta (or depending on
other parameters of a track). This is useful, for example, in
obtaining high accuracy for 'important' tracks (e.g. muons) and
accept less accuracy in tracking others (e.g. electrons). To use
this, you must create your own field manager which uses the
method

<informalexample>
<programlisting>
  void ConfigureForTrack( const G4Track * );
</programlisting>
</informalexample>

to configure itself using the parameters of the current track.
An example of this will be available in
examples/extended/field05.
</para>

</sect3>

<!-- ******************* Section (Level#3) ****************** -->
<sect3 id="sect.EMField.Pract.ParaScal">
<title>
Parameters that must scale with problem size
</title>

<para>
The default settings of this module are for problems with the
physical size of a typical high energy physics setup, that is,
distances smaller than about one kilometer. A few parameters are
necessary to carry this information to the magnetic field module,
and must typically be rescaled for problems of vastly different
sizes in order to get reasonable performance and robustness. Two of
these parameters are the maximum acceptable step and the minimum
step size.
</para>


<para>
The <emphasis role="bold">maximum acceptable step</emphasis> should be set 
to a distance larger than the biggest reasonable step. If the apparatus 
in a setup has a diameter of two meters, a likely maximum acceptable
steplength would be 10 meters. A particle could then take large
spiral steps, but would not attempt to take, for example, a
1000-meter-long step in the case of a very low-density material.
Similarly, for problems of a planetary scale, such as the earth
with its radius of roughly 6400 km, a maximum acceptabe steplength
of a few times this value would be reasonable.
</para>

<para>
An upper limit for the size of a step is a parameter of
<literal>G4PropagatorInField</literal>, and can be set by calling its
<literal>SetLargestAcceptableStep</literal> method.
</para>

<para>
The <emphasis role="bold">minimum step size</emphasis> is used during 
integration to limit the amount of work in difficult cases. It is 
possible that strong fields or integration problems can force the 
integrator to try very small steps; this parameter stops them from 
becoming unnecessarily small.
</para>

<para>
Trial steps smaller than this parameter will be treated with
less accuracy, and may even be ignored, depending on the
situation.
</para>

<para>
The minimum step size is a parameter of the MagInt_Driver, but
can be set in the contstructor of G4ChordFinder, as in the source
listing above.
</para>

</sect3>

<!-- ******************* Section (Level#3) ****************** -->
<sect3 id="sect.EMField.Pract.Known">
<title>
Known Issues
</title>

<para>
Currently it is computationally expensive to change the <emphasis>miss
distance</emphasis> to very small values, as it causes tracks to be
limited to curved sections whose 'bend' is smaller than this value.
(The bend is the distance of the mid-point from the chord between
endpoints.) For tracks with small curvature (typically low momentum
particles in strong fields) this can cause a large number of
steps

<itemizedlist spacing="compact">
  <listitem><para>
    even in areas where there are no volumes to intersect
    (something that is expected to be addressed in future development,
    in which the safety will be utilized to partially alleviate this
    limitation)
  </para></listitem>
  <listitem><para>
    especially in a region near a volume boundary (in which case it
    is necessary in order to discover whether a track might intersect a
    volume for only a short distance.)
  </para></listitem>
</itemizedlist>
</para>

<para>
Requiring such precision at the intersection is clearly expensive,
and new development would be necessary to minimize the expense.
</para>

<para>
By contrast, changing the intersection parameter is less
computationally expensive. It causes further calculation for only a
fraction of the steps, in particular those that intersect a volume
boundary.
</para>

</sect3>
</sect2>


<!-- ******************* Section (Level#2) ****************** -->
<sect2 id="sect.EMField.Spin">
<title>
Spin Tracking
</title>

<para>
The effects of a particle's motion on the precession of its spin
angular momentum in slowly varying external fields are simulated.
The relativistic equation of motion for spin is known as the BMT
equation. The equation demonstrates a remarkable property; in a
purely magnetic field, in vacuum, and neglecting small anomalous
magnetic moments, the particle's spin precesses in such a manner
that the longitudinal polarization remains a constant, whatever the
motion of the particle. But when the particle interacts with
electric fields of the medium and multiple scatters, the spin,
which is related to the particle's magnetic moment, does not
participate, and the need thus arises to propagate it independent
of the momentum vector. In the case of a polarized muon beam, for
example, it is important to predict the muon's spin direction at
decay-time in order to simulate the decay electron (Michel)
distribution correctly.
</para>

<para>
In order to track the spin of a particle in a magnetic field,
you need to code the following:

<orderedlist spacing="compact">
  <listitem><para>
    in your DetectorConstruction

    <informalexample>
    <programlisting>
#include "G4Mag_SpinEqRhs.hh"

  G4Mag_EqRhs* fEquation = new G4Mag_SpinEqRhs(magField);

  G4MagIntegratorStepper* pStepper = new G4ClassicalRK4(fEquation,12);
                                                       <emphasis role="bold">notice the 12 </emphasis>
    </programlisting>
    </informalexample>
  </para></listitem>
  <listitem><para>
    in your PrimaryGenerator

    <informalexample>
    <programlisting>
  particleGun-&gt;SetParticlePolarization(G4ThreeVector p)
    </programlisting>
    </informalexample>

    for example:

    <informalexample>
    <programlisting>
  particleGun-&gt;
  SetParticlePolarization(-(particleGun-&gt;GetParticleMomentumDirection()));

  // or
  particleGun-&gt;
  SetParticlePolarization(particleGun-&gt;GetParticleMomentumDirection()
                                         .cross(G4ThreeVector(0.,1.,0.)));
    </programlisting>
    </informalexample>

    where you set the initial spin direction.
  </para></listitem>
</orderedlist>
</para>

<para>
While the G4Mag_SpinEqRhs class constructor

<informalexample>
<programlisting>
  G4Mag_SpinEqRhs::G4Mag_SpinEqRhs( G4MagneticField* MagField )
     : G4Mag_EqRhs( MagField ) 
  {
      anomaly = 1.165923e-3;
  }
</programlisting>
</informalexample>

sets the muon anomaly by default, the class also comes with the
public method:

<informalexample>
<programlisting>
  inline void SetAnomaly(G4double a) { anomaly = a; }
</programlisting>
</informalexample>

with which you can set the magnetic anomaly to any value you
require.
</para>


<para>
For the moment, the code is written such that field tracking of
the spin is done only for particles with non-zero charge. Please,
see the Forum posting:

<literal>
http://geant4-hn.slac.stanford.edu:5090/HyperNews/public/get/emfields/88/3/1.html
</literal>

for modifications the user is required to make to facilitate
neutron spin tracking.
</para>

</sect2>
</sect1>