[1208] | 1 | \chapter{Electromagnetic Fields} |
---|
| 2 | |
---|
| 3 | \section{Creating a New Type of Field} |
---|
| 4 | {\sc Geant4} currently handles magnetic and electric fields and, |
---|
| 5 | in future releases, will handle combined electromagnetic |
---|
| 6 | fields. Fields due to other forces, not yet included in |
---|
| 7 | {\sc Geant4}, can be provided by describing the new field and the |
---|
| 8 | force it exerts on a particle passing through it. For the |
---|
| 9 | time being, all fields must be time-independent. This |
---|
| 10 | restriction may be lifted in the future. |
---|
| 11 | |
---|
| 12 | In order to accommodate a new type of field, two classes must |
---|
| 13 | be created: a field type and a class that determines the force. |
---|
| 14 | The {\sc Geant4} system must then be informed of the new field. |
---|
| 15 | |
---|
| 16 | \paragraph{A new Field class} |
---|
| 17 | A new type of Field class may be created by inheriting from |
---|
| 18 | G4Field |
---|
| 19 | \begin{verbatim} |
---|
| 20 | class NewField : public G4Field |
---|
| 21 | { |
---|
| 22 | public: |
---|
| 23 | void GetFieldValue( const double Point[3], |
---|
| 24 | double *pField )=0; |
---|
| 25 | } |
---|
| 26 | \end{verbatim} |
---|
| 27 | and deciding how many components your field will have, and what |
---|
| 28 | each component represents. For example, three components are |
---|
| 29 | required to describe a vector field while only one component is |
---|
| 30 | required to describe a scalar field. |
---|
| 31 | |
---|
| 32 | If you want your field to be a combination of different fields, |
---|
| 33 | you must choose your convention for which field goes first, |
---|
| 34 | which second etc. For example, to define an electromagnetic field we |
---|
| 35 | follow the convention that components 0,1 and 2 refer to the magnetic |
---|
| 36 | field and components 3, 4 and 5 refer to the electric field. |
---|
| 37 | |
---|
| 38 | By leaving the GetFieldValue method pure virtual, you force |
---|
| 39 | those users who want to describe their field to create a |
---|
| 40 | class that implements it for their detector's instance of |
---|
| 41 | this field. So documenting what each component means is required, |
---|
| 42 | to give them the necessary information. |
---|
| 43 | |
---|
| 44 | For example someone can describe DetectorAbc's field by creating |
---|
| 45 | a class DetectorAbcField, that derives from your NewField |
---|
| 46 | \begin{verbatim} |
---|
| 47 | class DetectorAbcField : public NewField |
---|
| 48 | { |
---|
| 49 | public: |
---|
| 50 | void MyFieldGradient::GetFieldValue( const double Point[3], |
---|
| 51 | double *pField ); |
---|
| 52 | } |
---|
| 53 | \end{verbatim} |
---|
| 54 | |
---|
| 55 | They then implement the function GetFieldValue |
---|
| 56 | \begin{verbatim} |
---|
| 57 | void MyFieldGradient::GetFieldValue( const double Point[3], |
---|
| 58 | double *pField ) |
---|
| 59 | { |
---|
| 60 | // We expect pField to point to pField[9]; |
---|
| 61 | // This & the order of the components of pField is your own |
---|
| 62 | // convention |
---|
| 63 | |
---|
| 64 | // We calculate the value of pField at Point ... |
---|
| 65 | } |
---|
| 66 | \end{verbatim} |
---|
| 67 | |
---|
| 68 | \paragraph{A new Equation of Motion for the new Field} |
---|
| 69 | Once you have created a new type of field, you must create an |
---|
| 70 | Equation of Motion for this Field. This is required in order to |
---|
| 71 | obtain the force that a particle feels. |
---|
| 72 | |
---|
| 73 | To do this you must inherit from G4Mag\_EqRhs and create your own |
---|
| 74 | equation of motion that understands your field. In it |
---|
| 75 | you must implement the virtual function EvaluateRhsGivenB. Given |
---|
| 76 | the value of the field, this function calculates the |
---|
| 77 | value of the generalised force. This is the only function that |
---|
| 78 | a subclass must define. |
---|
| 79 | \begin{verbatim} |
---|
| 80 | virtual void EvaluateRhsGivenB( const G4double y[], |
---|
| 81 | const G4double B[3], |
---|
| 82 | G4double dydx[] ) const = 0; |
---|
| 83 | \end{verbatim} |
---|
| 84 | |
---|
| 85 | In particular, the derivative vector dydx is a vector with six |
---|
| 86 | components. The first three are the derivative of the position |
---|
| 87 | with respect to the curve length. Thus they should set equal to |
---|
| 88 | the normalised velocity, which is components 3, 4 and 5 of y. |
---|
| 89 | \begin{verbatim} |
---|
| 90 | (dydx[0], dydx[1], dydx[2]) = (y[3], y[4], y[5]) |
---|
| 91 | \end{verbatim} |
---|
| 92 | |
---|
| 93 | The next three components are the derivatives of the velocity |
---|
| 94 | vector with respect to the path length. So you should write the |
---|
| 95 | "force" components for |
---|
| 96 | \begin{verbatim} |
---|
| 97 | dydx[3], dydx[4] and dydx[5] |
---|
| 98 | \end{verbatim} |
---|
| 99 | |
---|
| 100 | for your field. |
---|
| 101 | |
---|
| 102 | \paragraph{Get a G4FieldManager to use your field} |
---|
| 103 | In order to inform the {\sc Geant4} system that you want it to use |
---|
| 104 | your field as the global field, you must do the following steps: |
---|
| 105 | |
---|
| 106 | \begin{enumerate} |
---|
| 107 | \item Create a Stepper of your choice: |
---|
| 108 | \begin{verbatim} |
---|
| 109 | |
---|
| 110 | yourStepper = new G4ClassicalRK( yourEquationOfMotion ); |
---|
| 111 | // or if your field is not smooth eg |
---|
| 112 | // new G4ImplicitEuler( yourEquationOfMotion ); |
---|
| 113 | \end{verbatim} |
---|
| 114 | \item Create a chord finder that uses your Field and Stepper. You |
---|
| 115 | must also give it a minimum step size, below which it |
---|
| 116 | does not make sense to attempt to integrate: |
---|
| 117 | \begin{verbatim} |
---|
| 118 | yourChordFinder= new G4ChordFinder( yourField, |
---|
| 119 | yourMininumStep, // say 0.01*mm |
---|
| 120 | yourStepper ); |
---|
| 121 | \end{verbatim} |
---|
| 122 | \item Next create a G4FieldManager and give it that chord finder, |
---|
| 123 | \begin{verbatim} |
---|
| 124 | yourFieldManager= new G4FieldManager(); |
---|
| 125 | yourFieldManager.SetChordFinder(yourChordFinder); |
---|
| 126 | \end{verbatim} |
---|
| 127 | \item Finally we tell the Geometry that this FieldManager is |
---|
| 128 | responsible for creating a field for the detector. |
---|
| 129 | \begin{verbatim} |
---|
| 130 | G4TransportationManager::GetTransportationManager() |
---|
| 131 | -> SetFieldManager( yourFieldManager ); |
---|
| 132 | \end{verbatim} |
---|
| 133 | \end{enumerate} |
---|
| 134 | |
---|
| 135 | \paragraph{Changes for non-electromagnetic fields} |
---|
| 136 | If the field you are interested in simulating is not electromagnetic, |
---|
| 137 | another minor modification may be required. The |
---|
| 138 | transportation currently chooses whether to propagate a particle |
---|
| 139 | in a field or rectilinearly based on whether the particle is |
---|
| 140 | charged or not. If your field affects non-charged particles, you |
---|
| 141 | must inherit from the G4Transportation and re-implement the |
---|
| 142 | part of GetAlongStepPhysicalInteractionLength that decides whether |
---|
| 143 | the particles is affected by your force. |
---|
| 144 | |
---|
| 145 | In particular the relevant section of code does the following: |
---|
| 146 | \begin{verbatim} |
---|
| 147 | // Does the particle have an (EM) field force exerting upon it? |
---|
| 148 | // |
---|
| 149 | if( (particleCharge!=0.0) ){ |
---|
| 150 | |
---|
| 151 | fieldExertsForce= this->DoesGlobalFieldExist(); |
---|
| 152 | // Future: will/can also check whether current volume's field is Zero or |
---|
| 153 | // set by the user (in the logical volume) to be zero. |
---|
| 154 | } |
---|
| 155 | \end{verbatim} |
---|
| 156 | and you want it to ask whether it feels your force. If, for the sake |
---|
| 157 | of an example, you wanted to see the effects of gravity on a |
---|
| 158 | heavy hypothetical particle, you could say |
---|
| 159 | |
---|
| 160 | \begin{verbatim} |
---|
| 161 | |
---|
| 162 | // Does the particle have my field's force exerted on it? |
---|
| 163 | // |
---|
| 164 | if (particle->GetName() == "VeryHeavyWIMP") { |
---|
| 165 | fieldExertsForce= this->DoesGlobalFieldExist(); // For gravity |
---|
| 166 | } |
---|
| 167 | \end{verbatim} |
---|
| 168 | |
---|
| 169 | After doing all these steps, you will be able to see the effects of |
---|
| 170 | your force on a particle's motion. |
---|
| 171 | |
---|
| 172 | \section{Status of this chapter} |
---|
| 173 | |
---|
| 174 | 10.06.02 partially re-written by D.H. Wright \\ |
---|
| 175 | 14.11.02 spell check by P. Arce \\ |
---|
| 176 | |
---|
| 177 | |
---|
| 178 | |
---|
| 179 | |
---|
| 180 | |
---|
| 181 | |
---|