Geant4 /

Geant4 Simulation of the Fly-Past and Alpha Correction at the ISIS MSR Facility

Introduction

Spin polarised muons are used as probes of the local fields in the solid state at the ISIS MSR facility at the Rutherford Appleton Laboratory in the UK. Muons are implanted into a sample where they thermalise and then form some well defined Quantum State. This state could involve diamagnetic muon or be a paramagnetic form of muonium, it could also be a more complex chemically correlated form, involving the intrinsic host atoms or defects an impurities in the host. The quantum state could be related to a particular lattice site, or be in association with a defect, or it could also involve some form of diffusion (classical, quantum, band-like motion). The observation of the statistics of the spatial distribution of the decay positrons allows a spectroscopy of the local fields, and this yields geometrical, chemical and dynamical information of the well defined Quantum State under consideration.

For the case of MSR in the diamond host, one particular advantage is the ability to study muonium as a hydrogenic analogue atom. In fact, the bulk of the current understanding of the simple hydrogen impurity in diamond is derived from MSR experiments.

At the ISIS EMU facility, diamond samples are small compared to the muon beam cross-section. They are therefore in a special class of samples, where a correction must be made for a systematic false asymmetry in the measurement, due to a flypast component of the muon beam which does not stop in the sample. This asymmetry correction, termed “alpha”, is usually characterised using a silver sample of identical dimension to the small diamond sample. The alpha parameter, is field dependent.

There are several effects which could lead to the field dependent variation of alpha. At higher fields (in the EMU Facility), there is a focussing of the incident muons onto the sample due to the solenoidal nature of the field, the flypast fraction therefore changes. Also, there is an effect of the field on the transport properties of the decay positrons, should they come from the sample or the flypast beam dump, altering the count rate in the forward-backward decay positron detectors.

We have developed a detailed Monte Carlo simulation of the ISIS-EMU facility which includes the details on the muon beam (including distribution of rays in the muon beam ensemble), the small sample, the flypast component of the beam, the different histories of the muons, the transport properties of the decay positrons, the accurate description of the applied double Helmholtz field and the active detectors. All other materials are also modelled. The project has been implemented on the UJ Research Cluster.

Project Description – J Hartman

The project of Jonathan Hartman has been to implement the description of the applied double Helmholtz field in the GEANT simulation. The introduction of a magnetic field in GEANT is not well described in the manuals Jonathan had to reverse engineer the technology of incorporating fields in GEANT using various codes developed from other applications, and by gradually building up the complexity of his own code, starting from toy examples, to the full Simulation. At each sage, he had to run tests on the simulation using scenarios where one could either calculate the result analytically, or one had intuition about the qualitative behaviour of the simulation. Once this code was working, he had to integrate it into the existing simulation of the ISIS-EMU system (which had been a zero-field simulation), already developed by out group. As the simulation would need to be run many times for different field values automatically in batch mode on a Research Cluster, he also had to develop the machinery to initiate each simulation automatically on several cores with different parameter values for the field, and to assemble all the results together after the individual calculations were completed. The Simulation must therefore be adapted for parallel execution. Jonathan also had to document his modifications as well as produce a didactic note on the use of Magnetic Fields in GEANT.

The final stage will be to run the simulation in a parameter space of magnetic field and sample size, determining the alpha value for the Field dependent asymmetry correction. This simulation will be benchmarked against the experimental data for the Silver calibration sample. This is hoped to lead to insight into the alpha correction, especially at higher fields. This insight will then be used to interpret some high field data obtained for the behaviour of muonium in diamond. This data seems to indicate the ionisation of muonium at higher temperatures. It is then expected that a better analysis of this experiment can be made. The experiment has impact, as it means one would be determining the level of muonium in the band-gap of diamond. This is a fascinating concept, as the possibility of doping diamond using hydrogen is current. For example, it has been used successfully in wide band gap oxide material. In diamond, the situation is expected to be different though related, as hydrogen has a slightly different behaviour (mobile T-state and deeply trapped BC state, and very high probability of complexation with impurities).

It is expected that this work will then be publishable.

Status

  1. The B-Field can now be reliably described in GEANT for the ISI-EMU system.
  2. The code is running on the UJ-Research Cluster, in parameter search mode.
  3. Tests are being done on the Cluster version of the code
  4. Once the tests are done and we are sure that the results are satisfactory, we will do production runs – essentially, large scale number crunching of the simulation, to map out the field dependent behaviour of the systematic transport issues at the ISIS-EMU system.
  5. Finally, this data will be used to provide insight into the various components in the “alpha” correction, and to assist with removing systematic dependencies that affect the interpretation of data on the ionisation of the H-species in diamond.

The GEANT model of the ISIS-EMU facility. The double Helmholtz coils are not shown. The beam enters from the left axially, the diamond sample is in the centre of the forward and backward decay positron detectors.

Report from J Hartman on Results so far for B-Fields in GEANT

Implementation of Magnetic Fields

In Geant4 propagation of tracks in a field is simulated by numerically integrating the equation of motion of a particle in that field. In order to define a field for your application, you first need to create it by declaring it somewhere in the detector-construction user class. For example by adding the following line somewhere in the DetectorConstruction.cc source file:

G4UniformMagField* fEmField = new G4UniformMagField(G4ThreeVector(B[0],B[1],B[2]));

Don't forget to include the header file for the relevant field class (in the above example, this would be G4UniformMagField.hh). Then create a global field manager, for example:

G4FieldManager* FieldManager = G4TransportationManager :: GetTransportationManager() -> GetFieldManager();

Now set the previously defined field to the newly created field manager:

FieldManager -> SetDetectorField(fEmField);

The above lines would generate a global field defined in the whole world volume although it is possible to create a field that's only defined in a given logical volume. This is useful for computational efficiency if the calculation of the trajectories is resource consuming and you only need the field effects in a certain region of the detector. To do this, first create the logical volume in which the field is to be defined and then attach the field manager to this volume with the SetFieldManager command. If the logical volume is called LogicalVolume, then the code to attach the field manager would be: ''' LogicalVolume -> SetFieldManager(FieldManager, true);'''

The second parameter in SetFieldManager determines whether or not the field is also calculated in the daughter volumes of LogicalVolume. If you don't want the field to be calculated in the sub-volumes contained in LogicalVolume then set the second parameter of SetFieldManager to false.

Pre-defined Field Classes

Geant4 only has a few pre-defined field classes that can be declared when creating a field for the application. These classes are:

  • G4UniformElectricField
    • Defines a uniform electric field. This class is derived from the class G4ElectricField.

The rest below are all derived from G4MagneticField:

  • G4DELPHIMagField
    • This is the magnetic field from the DELPHI experiment.
  • G4HarmonicPolMagField
  • G4LineCurrentMagField
    • Magnetic field generated from a line current.
  • G4QuadrupoleMagField
    • Magnetic quadrupole.
  • G4UniformMagField
    • Uniform magnetic field.

In most cases however, the field that needs to be included will not be pre-defined. In that case you have to define your own field class.

User-defined Fields

In order to determine a user defined field, your new field class must contain a GetFieldValue method and your header file for the class must declare it as a public member function with the line:

inline void GetFieldValue(Point[4], *Field);

In the source code file where the field is defined, the GetFieldValue function should be used to define a field map with respect to time as position in the detector . The meaning of the parameters passed into the function is as follows: the first three elements of the one-dimensional array Point are the position coordinates the detector setup, (x = Point[0], y = Point[1], z = Point[2]), and the fourth component is time, t = Point[3]. The parameter Field is a pointer to the first element of an array that stores the vector components of the field in rectangular coordinates which are supposed to be the output of the GetFieldValue function. For a pure magnetic field B, the output field values are Bx = Field[0], By = Field[1] and Bz = Field[2]. For a combined electromagnetic field, the output field values are Bx = Field[0], By = Field[1] and Bz = Field[2] for the magnetic field components and Ex = Field[3], Ey = Field[4] and Ez = Field[5] for the electric field components. The GetFieldValue function is called at each step of the integration of the equation of motion for each charged particle in the field. The position and time, Point, are input values for the GetFieldValue function while the field, Field, is taken as output. In a user defined field class, the GetFieldValue function must contain lines that calculate the field values (passed in as the Field array) in terms of position (the first three elements of the Point array). The field values should only be dependent on Point[3] if the field is time varying. The user defined field class must be a derivative of G4MagneticField if its a pure magnetic field, a derivative of G4ElectricField if its a pure electric field or a derivative of G4ElectroMagneticField if its a combined electromagnetic field. If the field is anything other than an electromagnetic field (by electromagnetic, I mean including pure magnetic and electric fields), then you'd have to create your own kind of field with its own associated equation of motion. The new kind of field must be in the form of a class derived from the G4Field class.

Integrity checks by Visualisation Field lines and by Calculation of trajectories

In a simulation, tracks showing curved trajectories of particles are approximated by a series of straight lines called chords. The accuracy to which the set of chords approximate the track is determined by a parameter called the miss distance. The default value for the miss distance is 0.25 mm and is the maximum distance between the chord and the real track, see figure 4.6 in the book for application developers. The miss distance can be adjusted by the application developer. In addition to the miss distance, other important parameters include the delta intersection (the accuracy to which the track is calculated at the intersection of the boundary and the delta one step parameter (the accuracy of the endpoint of a given integration step. The last to parameters mentioned are set by SetDeltaIntersection and SetDeltaOneStep respectively, both derived from the globalFieldManager class. It's also possible to set the maximum and minimum step sizes. The maximum step size is set using SetLargestAcceptableStep, derived from GetPropagatorInField in the G4TransportationManager class. The minimum step size is passed as a parameter in the G4ChordFinder function when the chord finder is declared. Using SetLargestAcceptableStep just to visualise smooth tracks is not recommended because they are already calculated as smooth Geant 4. What you need to do to visualise smooth tracks is to use the command /vis/scene/add/trajectories smooth (either type this in at the command prompt after starting your program before doing a run or put this line in the macro. Also use a filter to get rid of any unwanted tracks you don't want to display. Here is an example of a smooth trajectory:


figure 1: Wireframe drawing of the ISIS-EMU facility with the spiral trajectory of a muon in a double Helmholtz magnetic field. The beam enters from the left axially, the diamond sample is in the centre of the forward and backward decay positron detectors.

Visualisation of Field Lines

If you are using the version 9.3 and later of Geant 4, a tool called G4BlineTracer may exist for the purpose of visualising field lines. For all earlier versions though, the code for this tool may be found in the field directory in the extended examples in the folder with the name BlineTracer. So if you want to use this tool in earlier versions of Geant 4, you will have to incorporate this code into application since it has not been put into Geant 4 kernel yet. To do so, you can just follow the instructions in the readme file thats included in the BlineTracer folder. However, in order for the code to work properly, there are a few changes that you need to make to Bline Tracer code itself which are not documented. Firstly if you look in the source code file G4BlineTracerMessenger.cc, you will see that in the first if statement, in the SetNewValue function (lines 128 and 129) is written the following:

if (command == BlineCmd) theBlineTool->ComputeBlines(1).

The number 1 in ComuteBlines(1) is a problem because according to the documentation in the readme file, the number you enter after the /vis/blineTracer/computeBline command should indicate the number of field lines to be calculated. However, the messenger class is passing the integer 1 to the function, no matter what number you type in. Despite what the readme file says, only one field line will be calculated. To fix this, change theBlineTool->ComputeBlines(1);, in the source code to theBlineTool->ComputeBlines(BlineCmd->GetNewIntValue(newValues));. After this change, there will be another problem. The field lines calculated are not independent and the the tool just draws the same field line multiple times. When a field line is drawn with this tool, two halves of the field line are drawn separately. As such, you wouldn't want the two halves to follow different trajectories. The way they are made to run along the same trajectory is by setting a flag called FirstPartOfBline (line 51) which you will find in the G4BlinePrimaryGeneratorAction.cc. The trouble is, this flag is not turned off once the second part of the field line is calculated. There are two more changes to the code that need to made before the it will work properly. First, in the header file G4BlinePrimaryGeneratorAction.hh go to line 68 and press enter to skip a line. Then add inline void SetFirstPartOfBline(G4bool* flag) and press enter to go onto the next line. Now add {FirstPartOfBline=flag;} then save and exit. For the last change, open the file G4BlineTracer.cc and notice the for loop starting at line 176. This is the loop that calculates the trajectories for the field lines. Move your curser to the end of line 197 and press enter to insert a line between the second theRunManager->BeamOn(1); and the right brace } which indicates the closure of the loop (the new line should be line 198). Now type SetFirstPartOfBline(false);, then save ,exit and compile your application with the modified Bline Tracer files in the correct directories. The tool should now work the way it's supposed to. SetFirstPartOfBline was the the command we recently created, by editing the header file G4BlinePrimaryGeneratorAction.hh, in order to reset the FirstPartOfBline flag after the second half of a give field line had been drawn. It is crucial that this is done after the second run is initiated by theRunManager->BeamOn(1); but before the closure of the main for loop.

The image below shows an example of how the visualisation of a field, using the Bline Tracer tool, should look. The specific field shown was a user defined magnetic one I created and is generated from a loop of electric current located at the centre of the figure.


Figure 2: Magnetic field lines of a field generated by a current loop.

These tests have been implemented over a variety of situations, for different field types, intensities, orientations and different particle energies. Only a subset of the results are shown.

Killing of tracks

In order to kill tracks in a simulation, you need to make use of the optional user class G4UserStackingAction. If this class is used, then its ClassifyNewTrack() method is called by G4StackManager every time a G4Track object is put in a particular stack by G4EventManager. The ClassifyNewTrack() method returns an object ClassificationOfNewTrack which indicates which stack the G4Track object is placed in. ClassificationOfNewTrack has four possible values, namely fUrgent, fWaiting, fPostpone and fKill. The fUrgent stack is the stack where the G4Track object is placed when that track is the next one to be calculated. The fWaiting is the stack where the tracks are stored if they don't immediately need to be simulated until the fUrgent one's empty, at which point they will be moved to that stack. The fPostpone stack is a stack where tracks can be placed if they are only to be simulated after the current event. Finally, if a value of fKill is returned for a particular track, then it is not placed in any stack but deleted or killed so that the track does not get simulated.

With the above information in mind, I have added the classes UserStackingAction and UserStackingActionMessenger classes to the musr application in to be able to kill the secondary positron tracks depending on whether or not they originated from the diamond. Note that you are only able to use the above method to kill secondary tracks rather than primary tracks since there is no memory address available to indicate the starting position of a track that can be passed to the G4EventManager class before the primary track is simulated.

Merging the Magnetic Field into the ISIS-EMU small sample simulation

This step has been completed. The possibility to initiate each simulation with a different B-Field has also been implemented, using the mechanism of the messenger class.

Developing the Cluster version of the Single Core code

This task has also been completed, and tests of the results with respect to the single core code are being done. There are issues about the versions of GEANT and other software over the cluster, as compared to the single test machine. This is leading to work needing to be done in on the code that facilitates the reports and collates the results from different cores.

Production runs on the cluster

This component will be initiated once we are satisfied with the integrity of the GEANT model.