grasshopper

This application is based on the geant4 development toolkit, which allows to program complex particle tracking (e.g. gamma, electrons, protons, etc) and particle-matter interaction MC simulations. Grasshopper is a simple geant4 application, where all the geometries and even generators parameters are defined in a gdml file, with the purpose of setting up quick and simple simulations. The goal is to allow users with no C++ and Geant4 knowledge quickly set up and run simulations.

Author: Areg Danagoulian

Creation time: 11/2015
Last update: continous

For copyright and licensing see files COPYRIGHT and LICENSE.

To install

The user is required to have the following

* xerces.  This will allow the GDML parser capability.
* For CMake Builds the User MUST Have CMake version 3.17 or higher.
* Built and installed geant4 libraries.  Also, in the cmake stage, the following flag needs to be passed:
`-DGEANT4_USE_GDML=ON`.  In some cases you also have to also add the location for Xerces with flags such as     
    `-DXERCESC_INCLUDE_DIR=/usr/local/include/ -DXERCESC_LIBRARY=/usr/local/lib/libxerces-c.so`.  See the geant4 instructions on how    to add Xerces for more detail on the paths.  Here's an example of what a Geant4 cmake command looks like:
    `cmake -DCMAKE_INSTALL_PREFIX=../geant4.10.05-install -DGEANT4_INSTALL_DATA=ON -DGEANT4_USE_GDML=ON ../geant4.10.05`
* ROOT -- optional.  Has been tested with version 6.22/07.  If you do not have ROOT the make process will recognize that and exclude        it from the build.
* When building grasshopper, the compiler might not find the GDML header files.  In that case just determine the actual file directories, and add them to the include list by appending `-I/directory_to_headers` to the `CPPFLAGS` env variable.

Version: grasshopper and been built against and tested with Geant4 10.6.

Important note: these days geant4 primarily works via the cmake framework. However grasshopper can also use the older Makefile framework. It is important that you source the appropriate shell script in geant4 directories to enable all the env. variables that are necessary for Makefile to work correctly. In my particular case I have the following line in my .bashrc file, please modify this accordingly for your build/configuration:

. /usr/local/geant4/geant4.10.05-install/share/Geant4-10.5.0/geant4make/geant4make.sh

If all the regular geant4 installations and configurations are ready, then the user can get the code by

> git clone https://github.com/ustajan/grasshopper.git

To build with GNUmake

> cd grasshopper

> make

To build with CMake

> mkdir RunGrasshopper && cd RunGrasshopper

> cmake ../../grasshopper && make -jN

To run

> grasshopper input.gdml output.root

Tutorial

The best way to learn how to use grasshopper is by using the tutorial on the project wiki, here. You can also get there by going to Wiki link above.

Input

input.gdml -- the geometry definition markup language file. This file defines
- the geometries of the objects and detectors, as well as their materials and the positioning
- the particle type, energy, and direction
- various computational optimization options
- output formats
input_spectrum.txt -- this is currently a fixed file name. The user can define the energy of the particle. If the energy is set negative, the code defaults to reading the energy spectrum from "input_spectrum.txt" and samples the energy from that spectrum.

Output

The code will generate three files:

output.root -- The code will always generate a root file, as specified in the input arguments.
output.dat -- depending on the settings in the gdml file, it can also generate a simple ascii file, output.dat, which can then be read into python/matlab/etc.
g4_00.wrl. This is the VRML visualization file. The settings for this file's content are defined in the gdml file. The wrl file can be loaded in paraview (http://www.paraview.org/download/), thus allowing for a rendering of the geometries and some particle tracks.

The root output structure (which is identical to the ascii output) can be somewhat complex. The output of the simulation can be thought of as a table -- where every line is an entry corresponding to a particular set of information within an event. Below is a set of entries in the ASCII output from a simulation of a 40meV neutron beam undergoing capture inside two 3He ionization chamber (numbered 5 and 37):

alt text

By enabling the BriefOutput flag in the gdml, it's is possible to get the shorter version of the above:

alt text

There are three type of entries:

IsSurfaceHitTrack. These entries are registered when a track enters the detector volume. To filter out these events, simply search for all entries where IsSurfaceHitTrack==1
regular tracks. At the end of an event all the tracks inside the detector can be recorded in the form of individual entries. All entries for which IsSurfaceHitTrack!=1 && IsEdepositedTotalEntry!=1 are regular track entries. These entries are useful for debugging and diagnostic purposes, as well as for general studies of the interaction physics of the detector.
IsEdepositedTotalEntry. At the end of an event all the relevant tracks' deposited energies is summed into E_deposited variable, and a dedicated entry is made in the tree with this information. By filtering IsEdepositedTotalEntry==1, all the deposition energy entries can be selected. This would be useful for studying the detector response function. These entries have the additional feature of including a concatenation of all the processes which contributed to the deposited energy. For example, for EventID=0 the energy deposition was achieved via the main track (EventGenerator), a compton scatter and a photoelectric effect (see the last line above).

The variables in the tree/table are as following:

E_beam -- the energy of the particle produced by the event generator
E_incident -- only for IsSurfaceHitTrack==1 entries this is the energy of the particle at the detector entrance. Good way to study the flux of the particles along a surface.
E_deposited -- the deposited energy in the detector
x_incident -- for IsSurfaceHitTrack==1 or SaveTrackInfo==1 the hit position on the detector, in mm
y_incident -- same
z_incident -- same
theta -- same (radians)
Time -- the (flight) time from the inception of the event
EventID -- this is the # of the event in the simulation history
TrackID -- for a particular event, many tracks may be produced. They will all have the same EventID and different TrackIDs
ParticleID -- for a particular track, this is the ParticleID, based on LLNL Particle Data Group definitions.
ParticleName -- the actual name (which is more useful than the ParticleID)
CreatorProcessName -- the name of the process which created this track
IsEdepositedTotalEntry -- the flag for total deposited energy entries
IsSurfaceHitTrack -- the flag for surface hit entries
detector# -- the number of the detector. The code assumes that all detector volumes are named det_phys#, where # is the detector number. The code trunkates the det_phys and thus extracts the detector number.

The CreatorProcessName for IsEdepositedTotalEntry==1 entries contains a concatenation of all the processes which contributed to the deposited energy. This is a good way to study the various effecst which, for example, contribute to the photopeak and the Compton continuum in a gamma detector.

Finally, it is also possible to request a brief output, in the form of EventID, Energy, ParticleName, and CreatorProcessName. To do this set the BriefOutputOn flag to 1, in the gdml input.

The output can be modified to be only limited for just one or two of these entry types. The gdml file allows to do this using the SaveSurfaceHitTrack, SaveTrackInfo, SaveEdepositedTotalEntry variables. For example, for a simulation where the flux is to be determined, only SaveSurfaceHitTrack needs to be set. For studying the energy deposition distribution for a particular type of detector the SaveEdepositedTotalEntry needs to be set. The use of these variables can significantly simplify the analysis of the MC output.

General status

In its current form, the code has been tested with <12MeV photons, neutrons, gammas, and electrons. Very general checks indicate that most processes are being simulated correctly.

A few known problems:

When running with neutrons >16MeV the code will hang.
The code is optimized to run in two configurations
- SurfaceHit analysis - providing energy/position/momentum information on the flux of particles through a "detector" surface. It will generate multiple entries for multiple tracks crossing a surface for a single MC event. However, if one of the tracks then backscatters and crosses the surface again it will be ignored.
- EnergyDeposited analysis - providing deposited energy for every event. Here the code is less sharp - if multiple tracks enter the detector, it will combine their deposited energies. As long as you specify which particle type's energy deposition is of interest (e.g. electrons for photons, or protons for neutrons) this is sufficiently accurate as most detectors will not be able to time-resolve the multiple hits from tracks that originate from a single track. This, however can fail for a) very fast detectors or b) when you have a process which produces a low energy electron and a photon. If you have to deal with complex processes like (b) then grasshopper is probably not the best thing to use (least you want to modify the code or see if various cuts on particle type allow you to isolate individual responses).
More importantly, the above circumstance results in some (mostly minor) issues where the user tries to determine the energy deposition for a whole range of energies of incident particle. If the incident particle creaters an electromagnetic or hadronic shower, then an E_incident vs. E_deposited comparison may yield strange results, as E_deposited>E_incident in some cases - this is due to the fact that E_deposited may have captured deposited energy from multiple tracks, while E_incident is the energy of an individual track.

To do

Below is a prioritized list of future tasks.

* OUTPUT
  + Implement "simple" ascii text output, along with ROOT  -- DONE
  + Make the code be able to reliably switch between ROOT (when ROOT is available) and ASCII output -- DONE
* Write the python/javascript front end to the gdml? <<===== need a UROP
* Check the physics, by comparing to data and validated geant4 simulations -- DONE
  + Did some basic comparisons of gamma transmission to exp(-mu*x)
* automatic wrl generation (even in batch mode), limited to 300 tracks -- DONE
* Implement the keeper list, which should allow to dramatically speed up the simulations -- DONE
* Clean up the code:
  + Replace the char arrays in TTree with actual strings -- DONE

Requirements for the front end

The front end will read in a simple ASCII text file with "human readable" descriptions, parse it, populate a gdml file based on some template, and run the geant4 simulation. The ASCII input file needs to have the following fields and sections

* Particle Event Generator parameters
  + Particle type definition.  Ideally, this should include ions.  A list should be provided in a separate .md document
  + Particle energy
  + Directionality.  For simplicity, limit it to:  parallel;  isotropic.
  + Beam source x,y,z
  + Beam width
* Geometries
  + Number of geometries -- N (default:1).  No daughter geometries.
    If people want to do that, they can modify the gdml directly.
  + Materials: *only* limited to NIST materials, provided by a (long) list.
    If people want custom materials, they can modify the gdml directly.
  + Geometry1:
    - Sensitive?  If yes, its energy deposition will be measured and added to the output.
      An additional branch needs to be added, which will list the detector name.
    - volume -- shape
    - logical -- add material.
    - physical placement in the world.  If this is a sensitive volume, then the name needs to
      end with _number.  This number will go into the tree branch.
  + Geometry2, etc...
* Keeper list and cuts
  + This list should contain the following:  particle type; energy threshold E_th; volume volname
  + Now, for every particle type, only tracks with E>E_th *inside* volname will be tracked
  + Everything else will be killed
  + Example of how this will look:
    - neutron:  0.1 MeV,  everywhere
    - gamma  :  0.1 MeV,  everywhere
    - e-     :  0.1 MeV,  detector_volume
    - proton :  1.0 MeV,  detector_volume
  + The code will generate a vector, and for every step will loop over the vector, checking if these
    conditions are met.  If not, it will kill the track.

ustajan / grasshopper

readme